In this paper we report the results of tests of the ambient-noise surface wave seismic tomography method at Stillwater Canada Inc.’s Platinum Group Element-Copper (PGE-Cu) deposit at Marathon in Ontario. This survey demonstrated the potential of the technique to help determine the geological setting of the deposit and, in addition, located a previously unknown exploration target. The method uses cross-correlations of recorded ambient seismic noise to measure velocities of Rayleigh surface waves at specific frequencies that are then inverted to produce a shear wave velocity model of the subsurface.

Introduction: Ambient-noise surface wave seismic tomography

Seismic imaging has been used for many years to probe the subsurface of the earth. The standard active-source method employs a controlled source of seismic energy - explosives or a Vibroseis truck - and uses reflected seismic waves for the imaging. This method is potentially useful for exploration of new ore bodies, but is relatively expensive and commonly impractical because of cost, access issues, and its environmental impact.

In the past decade, Michel Campillo and co-workers have developed a method using the normal background vibrations (ambient seismic noise) to image the subsurface. Correlations of seismic noise recorded at two stations (sensors) are used to reconstruct the response of the Earth as if one of the receivers acted as a source and the other as a receiver (Figure 1). Using this approach, it is possible to image the subsurface without an active seismic source, in an inexpensive and non-invasive way.

Fig. 01
Figure 1. Example of correlation between two sensors A and B that recorded seismic noise (i.e. ground vibrations indicated by the arrows) during a given period of time). The correlation of the record, symmetrical in this case, is represented on the right. The left side represents the response as if B were the source and A the receiver (B → A); on the right, the roles are reversed (A → B).

Ambient seismic noise has been studied for decades and it appears that its origin is different for different period bands (Gutenberg, 1958). In general, noise with a period greater than 1 second (i.e. lower than 1 Hz) is due to meteorological conditions such as waves generated by storms on lakes or ocean swell. Noise with a period less than 1 second is mainly due to human activities such as cars and machinery. Seismic noise is dominated by surface waves (e.g. Campillo & Paul, 2003; Shapiro & Campillo, 2004) of which two types can be observed: Rayleigh waves and Love waves. Both have the same propagation direction but the ground motion is different: Rayleigh waves induce ellipsoidal motion in a vertical plane (like ocean waves, Figure 2a), while Love waves induce horizontal motion (Figure 2b). Hence a triaxial seismometer records Rayleigh waves on its vertical and radial components and Love waves on its transverse component. Love waves are shear waves and Rayleigh waves comprise both shear and compressional waves.

Fig. 02
Figure 2. Surface wave propagation and particle motion of a) Rayleigh wave (anonymous (https://commons.wikimedia.org/wiki/File:Rayleigh_wave.jpg), "Rayleigh wave“, https://commons.wikimedia.org/wiki/Template:PD-US) and b) Love wave (Nicoguaro (https://commons.wikimedia.org/wiki/File:Love_wave.svg), https://creativecommons.org/licenses/by/4.0/legalcode).

An important property of surface waves is dispersion: since the subsurface generally consists of zones of different seismic velocities, waves of different wavelengths travel through different zones and thus travel at different velocities. Since wavelength and frequency are related, dispersion is the effect by which waves of different frequencies travel at different velocities. Typically, seismic velocities increase with depth as the rock becomes less weathered and more competent, and so typically the low frequency waves (travelling through both shallow slow zones and deeper faster zones) are faster than high frequency waves which travel through the slow zones close to the surface.

Fig. 03
Figure 3. a) Example of difference between group (denoted as U) and phase (denoted as c) velocities. Group velocity is the velocity of the envelope and phase velocity is the velocity of the carrier. Time passing is denoted by t1, t2 and t3; b) relationship between phase and group velocities: the group velocity has a minimum where the phase becomes steep. Figures from Stein & Wysession (2003).

Two types of velocities can be measured: the velocity of the envelope of the signal, called the group velocity, and the velocity of the carrier, called phase velocity (Figure 3). The comparison of a dispersive wave at different travel times shows that the envelope of the signal propagates at a different speed from the carrier. The relationship between the group velocity and the phase velocity implies that if a wave is not dispersive then the phase and group velocity are equal (Figure 3b). In this paper, group velocities are determined using the Frequency Time Analysis (FTAN) introduced by Dziewonski et al. (1969) and average phase velocity using a FK analysis.

Surface waves are recovered by comparing the seismic noise records registered at two different seismic stations. Once surface waves have been retrieved from seismic noise, it is possible to compute tomographic images of the subsurface. Two tools are used to compute shear wave velocity models from surface waves reconstructed with seismic noise. The method of Mordret et al. (2013) is used for 2-D spatial travel-time tomography and that of Sambridge (1999) and Mordret et al. (2014) is used for 1-D depth tomography.

Geological Setting

The Marathon deposit is one of several gabbroic to ultramafic hosted disseminated PGE-Cu deposits located around the outer margin of the Coldwell Complex. This Complex is located along the north eastern shore of Lake Superior, near the town of Marathon and is associated with the Midcontinent Rift (Figure 4). It was dated at 1108 ± 1 Ma by Heaman et al., (2007) and between 1107.7 to 1105.3 Ma by Good and Dunning (2018). The Complex, approximately 25 km in diameter, is the largest alkaline intrusion in North America. It is composed of three centres of predominantly alkaline magmatism that intruded the Archean greenstone terrane (Mitchell and Platt, 1977). Centre 1 is composed of Augite Syenite, Quartz Syenite and the Eastern Gabbro Suite. The Eastern Gabbro is composed of numerous gabbroic and ultramafic intrusions of the Metabasalt, Layered and Marathon Series. Geological and geophysical evidence suggests that the Marathon Series intruded along conduits within a large plumbing system beneath the thin flat-lying syenite intrusions.

Fig. 04
Figure 4. Geological map of the region (Good et al., 2015 after Miller and Nicholson 2013).

The mineralization consists of mainly disseminated PGE-bearing Cu-rich sulfides hosted in the Two Duck Lake Gabbro of the Marathon Series. The deposit has dimensions of approximately 3000 m strike length and a variable thickness up to 200 m. Like the intrusion, it dips to the west at about 35 degrees. There are three mineralized zones: the Footwall Zone at the lower contact contains blebby to semi-massive sulfides with low Cu/S ratios; the Main Zone, the thickest and most continuous, is characterized by disseminated interstitial sulfide minerals and the W-Horizon is a zone of high PGE, low S mineralization located at the south end of the deposit. Estimated total proven and probable reserves of 91 Mtonnes at 0.83grams/tonne Pd, 0.24grams/tonne Pt, 0.085grams/tonne Au and 0.247% Cu, containing 3.39 Mounces of combined Pd+Pt+Au and 497 Mpounds of Cu (Murahwi, 2010).

Fig. 05
Figure 5. Simplified geological map of the Marathon deposit showing the position of the seismic array.


a) Phase I - Noise test

A noise test using an array of 31 vertical component stations was performed to determine the Power Spectral Density and azimuthal distribution of local ambient seismic noise. Use of vertical component stations was dictated by cost and availability of recording equipment for the Phase II survey and allowed only Rayleigh wave analysis to be performed. The station distribution and beamforming results (Figure 6) show that the source of the high-frequency noise lies to the SW and is probably waves on Lake Superior. Low-frequency noise comes from the NE, probably ocean swell on the North Atlantic.

Fig. 06
Figure 6. Station distribution and beamforming analysis averaged over the 26 days of acquisition, in five different period bands. The panel on the lower right shows the stations used in the beamforming analysis.

b) Phase II – demonstration of the ambient noise surface-wave technique

The dataset consists of 26 days of ambient seismic noise recorded over 90 vertical component stations (Figure 7). The nominal distance between the sensors is 300 m.

The Power Spectral Density (PSD) and spectrograms for each individual station show for most stations relatively stable seismic noise with the distinct pattern of oceanic storms at periods above 1 s. In the period band used for the tomography (0.15 - 1.5 s), the noise is made of discrete bursts, probably related to local weather.

Fig. 07
Figure 7. Map of the 90 1-C stations on top of the local topography. The blue lines show the surface extension of the gabbro intrusion (Good et al., 2015).

Noise cross-correlations and noise source distribution

Correlation functions are computed from daily ambient noise records. Signals are first down-sampled from a sample rate of 4 ms to 20 ms and normalized in the frequency band [0.1 - 20] Hz. Figure 8 shows correlations between a reference station (00035 and 00110 acting as a virtual source) versus the other stations in the array. Correlations are sorted as a function of distance from the reference station. Propagation of Rayleigh waves is clearly visible across the array. Most of the noise energy is coming from the southwest (acausal, or negative time, energy on station 00035 correlations and causal, or positive time, energy on 00110 station correlations) in the direction of the Great Lakes (Xu et al., 2017). Some low frequency energy is coming from the northeast. This is confirmed by the beamforming analysis performed in different frequency bands (Figure 6). In the following we will work on the symmetric components of the correlations, i.e., the average of the causal and acausal parts.

Fig. 08a Fig. 08b
Figure 8. Noise cross-correlations for stations a) 00035 and b) 00010 acting as a virtual source, versus the other stations for the frequency-band [(0.1-20) Hz].

Average dispersion curves over the entire network

We now build an average seismic section by summing acausal and causal portions of the correlations and binning the correlations in fixed distance intervals (every 50 m). Doing so, we assume that the underlying velocity model is mostly 1D. By stacking the large amount of data in each distance bin, the signal-to-noise ratio (SNR) is highly increased, which allows us to extract the fundamental mode (i.e. the first harmonic) average phase velocity dispersion curve (Figure 9).

Fig. 09
Figure 9. Stacked section of noise cross-correlations filtered between (5-10) Hz. The fundamental mode of the Rayleigh waves is clearly apparent with no apparent interference from overtones of the Rayleigh wave.

A frequency-wavenumber (FK) analysis of the stack section allows us to accurately pick the average phase velocity dispersion curve for the entire network (Figure 10a). A mean group velocity dispersion curve is then calculated as a derivative of the mean phase velocity dispersion curve according to

Equation 01

where U is the group velocity, c is the phase velocity and ω the frequency. These mean dispersion curves are presented on Figure 10b. We verified that the mean group velocity dispersion curves calculated from FK plots are coherent with the mean group velocity dispersion curves from FTAN picking described later on. The derivation of group velocity from phase velocity involves a differentiation; we smoothed the phase velocity prior to the computation to obtain a stable group velocity dispersion curve. We therefore get an average group velocity dispersion curve which exhibits the correct average velocity value but lacks some variation along the different periods (see Figure 11a).

Fig. 10
Figure 10. a): FK analysis of stack sections of noise-correlations with picked phase velocity. b): Mean dispersion curves calculated from FK plots: Phase velocity dispersion curve and group velocity dispersion curve of Rayleigh waves (fundamental mode).

Group velocity dispersion curve picking and quality control

The Rayleigh wave group velocity dispersion curves for the 4005 possible vertical-vertical (ZZ) correlations have been manually picked using the FTAN algorithm. FTAN is a set of narrow Gaussian band-pass filters applied to the correlation. The maximum of the envelopes of the band-pass filtered correlation defines the group velocity dispersion curve for this correlation. We then rejected all the dispersion curves for station pairs separated by less than 900 m to avoid larger uncertainty measurements due to near source effects. Finally, after computing the statistics of all remaining dispersion curves (Figure 11a) we rejected all the dispersion curves for which at least one point falls outside the boundaries defined in Figure 11b, i.e. the values larger (or smaller) than the most probable velocity value (called mode) ± 1 km/s. In Figure 11b, the mode value is shown by the blue curve and the 1 km/s limits by the dashed white curves. The background of the figure shows the probability density function of the Rayleigh wave group velocity dispersion curves, before 11a and after the 11b quality control.

Fig. 11
Figure 11. a) Probability density distribution of all Rayleigh wave group velocity dispersion curves. The red, black and blue curves show the mean, median and most probable dispersion curve, respectively. The green curve is the average group velocity dispersion curve derived from the FK analysis. b) Probability density distribution of the Rayleigh wave dispersion curves left after the rejection procedure. The red, black and blue curves show the mean, median and most probable dispersion curve, respectively. The dashed blue curves are the limits used to reject outlier dispersion curves.

The number of dispersion curves left at each period and used for the tomography is shown in Figure 12. It goes up to 1380 dispersion curves for the period band [0.15 - 1.5] seconds. The number of pairs is limited due to the large (900 m) minimum inter-station distance rejection criterion.

Fig. 12
Figure 12. Number of measurements kept at each period after the quality control.

The azimuthal distribution of the Rayleigh wave paths at each inverted period is shown in Figure 13. We can see that the azimuths in the southwest/northeast direction are dominating because the dispersion curves kept during the QC are mostly oriented in this direction, the direction of the main noise sources.

Fig. 13
Figure 13. Azimuthal distribution of the Rayleigh wave paths used for the tomography at each frequency. The outer circle defining the number of Rayleigh wave paths is given by the number on the right. The inner circle value is half this number.

Tomography at different frequencies

The dispersion curves are inverted in the period band [0.15-1.5] seconds (with a regular step of 0.1 s) and regionalized into a regular grid of 100 m by 100 m rectangular cells (in longitude and latitude directions) using the Mordret et al. (2013) approach. Figure 14 shows the results for the Rayleigh wave group velocity maps.

Fig. 14
Figure 14. Rayleigh wave group velocity maps inverted at different periods. The blue dots show the station locations, the black contours show the surface extent of the gabbro intrusion.

Figure 15 shows the density of rays in each cell of the model at each period for Rayleigh waves. A density of more than 20 rays per cell permits us to achieve a nominal lateral resolution of about 250 m for most of the area, with a minimum resolution of about 500 m at the edges of the area. The density is higher (up to 120 rays per cell) in the middle of the array.

Rayleigh wave group velocity maps show evolving structures with increasing periods. The map at 0.25 s correlates particularly well with the surface extent of the gabbro intrusion.

Fig. 15
Figure 15. Ray density maps at each period for the Rayleigh wave group velocity. Depth inversion and estimation of the model robustness.

The full set of group velocity maps describes local dispersion curves at each single cell of the map. These local dispersion curves are then inverted at depth to obtain local 1D shear wave velocity (Vs) profiles. The ensemble of every best local 1D model constitutes the final 3D Vs model. We used a Monte-Carlo approach to invert the dispersion curves at depth because of the strong non-linearity of the problem and the absence of an accurate a priori starting velocity model. This Monte-Carlo approach is based on a Neighborhood Algorithm (Sambridge, 1999) and is described in detail by Mordret et al. (2014). This algorithm has the advantage of exploring entirely the model space using a Voronoï cell discretization while increasing the focus on its most promising areas as the number of iterations increases. The ensemble of best models can be used to estimate some uncertainties on the model parameters. We jointly inverted both Rayleigh group velocity dispersion curves and the mean phase velocity dispersion curve to stabilize and complete the inversion. In total, we inverted 2 dispersion curves for each single cell. These curves are defined in Table 1.

Wave type Velocity Mode Period (s) Measurements
Table 1: Dispersion curves used in depth inversion for each point of the grid.
Rayleigh Phase Fundamental 0.1 - 1.2 FK
Rayleigh Group Fundamental 0.15 - 1.5 FTAN

For each dispersion curve presented in Table 1, we define error bars. The error bars of the mean dispersion curves calculated with FK analysis are twice as high as the error bars calculated with FTAN analysis (FK: +/- 200 m/s and FTAN: +/- 100 m/s). This equalizes the contribution of the mean dispersion curves in the inversion. The depth inversion searches for the optimum parameters that minimize the difference between the synthetic dispersion curves and the observed ones. The 1D model is discretized with 3 homogeneous layers with various thicknesses. During the inversion, the P-wave velocity is scaled to Vs using a Vp/Vs ratio of 1.73 and the density is scaled to Vp using the empirical relationship of Brocher (2005). In total, we invert for 5 parameters: velocities in 3 layers and the thicknesses of the 1st and the 2nd layer (Table 2).

Parameters Units Min Max
Table 2: Parameterization of the local 1D Vs models.
Vs1 m/s 2500 3500
Vs2 m/s 2800 4000
Vs3 m/s 2800 3800
Depth1 m 5 700
Depth2 m Depth1+100 Depth1+1000

We sampled a total of 21,000 models. The final model is the average of the 1000 best models (with the lowest misfits). This Monte-Carlo technique also allows us to estimate uncertainties for our model as the standard deviation of the distribution of the 1000 best velocity models at each depth. We interpolate the 3-layer 1D models every 1 meter at depth to obtain the final velocity cube.

Fig. 16
Figure 16. Local depth inversion. The data are in red, the sampled models and associated synthetic dispersion curves are in grey (the shade depends on the misfit value) and the best model is in green. a) The point corresponding to the green star. b) The point corresponding to the red star. The fit to the average phase velocity curve is poor because in this area the 1D assumption we took to estimate the average dispersion curve does not hold. c) Map of the average normalized misfit for the 1000 best models at each geographical point.

Figure 16 shows two examples of local depth inversion along with the final misfit map. We see that there is a band of higher misfit in the middle of the model (Figure 16c). In this area, the high misfit comes from the poor fit of the average phase velocity that we measured assuming a 1D average velocity model. The dipping slab in this area violates this assumption: the data are sensitive to the thickness of the gabbro intrusion slab and to the material below it, introducing a lower velocity layer which is poorly constrained in terms of velocity. However, the local group velocity is sensitive to the slab geometry. We chose to keep narrow ranges for the inversion parameters to reduce the non-uniqueness of the solutions and avoid over-fitting our results. Although the misfit is larger in these areas, the inverted model parameters are still robust. Given the periods that we used for the inversion, the model is resolved between ~100 m and 1500 m depth.

Final 3D Vs model

In the following figures we present horizontal (Figure 17) and vertical slices (Figure 18) into the final 3D cube. We observe that the gabbro intrusion slab is well resolved as a higher velocity layer, dipping toward the west. The contact between the upper syenite and lower gabbro intrusion is well resolved but the bottom of the gabbro intrusion is less clear and seems to be picked up by the data only in the middle of the array. However, in this area, the misfit is higher so the interpretation of the slab bottom should be made cautiously.

Fig. 17
Figure 17. Slice at 250 m depth in the velocity model. The locations of the vertical slices presentedin Figure 18 are shown by the black lines. The red contours show the surface trace of the gabbrointrusion.
Fig. 18
Figure 18. Vertical slices in the Vs model. The locations of the slices are presented in Figure 17.

A high velocity anomaly seen on Figure 18 on A and D vertical slices at ~700 m depth is associated with larger uncertainties (Figures 19 and 20) and its true nature is unknown. In order to have a better view of the slab upper surface, we plot the iso-velocity surface at 3250 m/s in Figure 21.

Fig. 19
Figure 19. Slice at 250 m depth in the uncertainty model. The locations of the vertical slices presented in Figure 20 are shown by the black lines.
Fig. 20
Figure 20. Vertical slices in the Vs uncertainty model. The locations of the slices are presented in Figure 19.
Fig. 21
Figure 21. Iso-velocity surface at 3250 m/s. The colour indicates the depth of the surface. We clearly image the dipping of the slab toward the west.


We used 26 days of passive seismic records to retrieve the fundamental mode of the Rayleigh waves propagating across the Marathon area. Most of the noise is coming from Lake Superior.

We jointly inverted local group velocity dispersion curves and the average phase velocity dispersion curve to obtain a robust 3D shear-wave velocity model of the studied area resolved down to 1.5 km. Preliminary interpretation, to be developed in a subsequent paper, is that we were able to image 1) the contact of the gabbroic intrusion which dips to the west, 2) a near-surface layer that could in part be the syenite in the centre of the complex; and 3) a high-velocity zone at ca. 700 m depth that could represent mineralized, olivine-rich gabbro. This interpretation generally agrees with results from the previous core hole drilling program in the survey area.



The initial two phases were carried out by a partnership between Stillwater, UGA and Sisprobe and were supported by a Booster Grant from the EIT Raw Materials. The final phase III was one of the main activities of PACIFIC (grant agreement No 776622), a project in the European Union’s Horizon 2020 program. We also acknowledge support from the European Union’s Horizon 2020 research and innovation program for this project.

About the Author(s)

Dan Hollis is a seismologist and CEO of Sisprobe SAS, a company that specializes in ambient seismic noise imaging. Dan received a Bachelor of Science in Earth Science from California State University at Long Beach. He is a member of SEG, EAGE, AGU, AAPG, and is a registered geophysicist.


Brocher, T. M., 2005. Empirical relations between elastic wavespeeds and density in the Earth’s crust, Bulletin of the seismological Society of America, 95(6), 2081–2092.

Campillo, M. & A. Paul, 2003. Long-range correlations in the diffuse seismic coda. Science, 299, 547-549.

Dziewonski, A., Bloch, S. Landisman, N. 1969. A technique for the analysis of transient seismic signals. Bull. Seism. Soc. Am. 59, 427-444.

Ekström, G., Abers, G.A. Webb, S.C. 2009 Determination of surface‐wave phase velocities across USArray from noise and Aki's spectral formulation. Geophysical Research Letters, 36, doi.org/10.1029/2009GL039131.

Good, D.J., Epstein, R., McLean, K., Linnen, R., & Samson, I.M., 2015. Evolution of the Main Zone at the Marathon Cu-PGE sulfide deposit, Midcontinent Rift, Canada: spatial relationships in a magma conduit setting, Economic Geology, 110(4), 983–1008.

Good, D.J. and Dunning, G.R., 2018. Geochronology of mafic magmas of the Coldwell Complex, Midcontinent Rift, and age constraints on magma evolution in relation to Cu-PGE mineralization. Abstract in RFG 2018, Vancouver.

Gutenberg, B., 1958, Two types of microseisms: Journal of Geophysical Research, 63, 595–597.

Heaman, L.M., Easton, M., Hart, T.R., Hollings, P., Macdonald, C.A. and Smyk, M. 2007. Further refinement to the timing of Mesoproterozoic magmatism, Lake Nipigon region, Ontario, Canadian Journal of Earth Sciences 44: 1055-1086.

Mitchell, R.H., and Platt, R.G., 1977, Field guide to the aspects of geology of the Coldwell alkaline complex: Institute on Lake Superior Geology, Technical Report.

Mordret, A., Landès, M., Shapiro, N., Singh, S., Roux, P., & Barkved, O., 2013. Near-surface study at the Valhall oil field from ambient noise surface wave tomography, Geophysical Journal International, 193(3), 1627–1643.

Mordret, A., Landès, M., Shapiro, N., Singh, S., & Roux, P., 2014. Ambient noise surface wave tomography to determine the shallow shear velocity structure at Valhall: depth inversion with a Neighbourhood Algorithm, Geophysical Journal International, 198(3), 1514–1525.

Murahwi, C. Shoemaker, S. Gowans, R., Lemieux, J., and Jacobs, C., 2010, Technical report on the updated feasibility study for the Marathon PGM-Cu Project. Micon International Ltd. NI43-101 report.

Sambridge, M., 1999. Geophysical inversion with a neighbourhood algorithm—I. searching a parameter space, Geophysical Journal International, 138(2), 479–494.

Shapiro, N. M. and M. Campillo, 2004, Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, Geophysical Research Letters, VOL. 31, L07614, doi:10.1029/2004GL019491.

Shapiro, N. M., M. Campillo, L. Stehly, and M. H. Ritzwoller, 2005, High resolution surface wave tomography from ambient seismic noise: Science, 307, 1615–1618.


Join the Conversation

Interested in starting, or contributing to a conversation about an article or issue of the RECORDER? Join our CSEG LinkedIn Group.

Share This Article