TY - JOUR AU - Gudmundsson,, Olafur AB - Summary A double-correlation method is introduced to locate tremor sources based on stacks of complex, doubly-correlated tremor records of multiple triplets of seismographs back projected to hypothetical source locations in a geographic grid. Peaks in the resulting stack of moduli are inferred source locations. The stack of the moduli is a robust measure of energy radiated from a point source or point sources even when the velocity information is imprecise. Application to real data shows how double correlation focuses the source mapping compared to the common single correlation approach. Synthetic tests demonstrate the robustness of the method and its resolution limitations which are controlled by the station geometry, the finite frequency of the signal, the quality of the used velocity information and noise level. Both random noise and signal or noise correlated at time shifts that are inconsistent with the assumed velocity structure can be effectively suppressed. Assuming a surface wave velocity, we can constrain the source location even if the surface wave component does not dominate. The method can also in principle be used with body waves in 3-D, although this requires more data and seismographs placed near the source for depth resolution. Persistence, memory, correlations, clustering, Interferometry, Volcano seismology INTRODUCTION Many Earth processes, especially in volcanic areas, produce complex seismic signals lacking a clear onset, such as emergent events or continuous tremor. We cannot pick their onsets to locate their sources. Instead, amplitude decay with distance has sometimes been used (e.g. Battaglia & Aki 2003; Di Grazia et al. 2006). This strategy requires simplistic assumptions about geometrical spreading, site magnification and attenuation and does not work well for the tremor data at Katla volcano, which we use as an example in this paper (Sgattoni 2016). Another strategy is to use array methods, such as beam forming (Rost & Thomas 2002), frequency-wavenumber analysis (Capon 1973) or semblance methods (Furumoto et al. 1992; Almendros 2003; Konstantinou & Schlindwein 2003) on station subarrays. Array methods require dense sampling of the wave field to yield reliable direction and slowness or vectorial wave number. Combining such estimates at multiple dense arrays constrains the source location. In classical beam forming (Rost & Thomas 2002), after appropriate delay that corresponds to a given trial source location, data series from all stations are summed, the result squared, and this summed over time. The result is identical to the sum of all autocorrelations and appropriately delayed cross correlations of the individual time series. The autocorrelations add noise, but contain no timing information, so excluding these and summing only the cross correlations helps suppress noise without loss of relative time information. Interstation cross correlation has become a common tool for analysis of seismic ambient noise over the past decade. Numerous authors have used cross-correlation analysis to assess the spatial distribution of ambient noise sources. Shapiro et al. (2006) located a peak near the primary microseismic peak at 26 s period in the Gulf of Guinea by back projecting noise correlations. Gudmundsson & Brandsdottir (2010), Zeng & Ni (2010), Ballmer et al. (2013) and Droznin et al. (2015) used a similar approach to locate persistent sources of geothermal noise in Iceland, microseisms around Japan and volcanic tremor in Hawaii and Kamchatka, respectively. Haney (2010) measured differential traveltimes extracted from cross correlations of recordings of very-long-period volcanic tremor at different stations and inverted those for the source location. Haney (2014) derived an analytical description of the array response of the squared envelope of stacked, time-shifted seismograms for monochromatic waves and deconvolved it from his back-projection results to resolve multiple sources and study their time evolution. All such methods require detailed information about seismic velocity in the area. Because this information is often insufficiently precise, envelopes or various integral measures of cross correlations were stacked rather than the cross correlations themselves in a number of the above studies. If we consider a tremor or noise source as a point source which is incoherent in time and a simple medium between that source and all recording seismographs, then the recorded tremor (noise) will consist of a continuous interference pattern between the same source functions at random time intervals, but at a fixed delay for one station compared to any other corresponding to the differential traveltime from the source to the two stations. Recordings of the tremor at two stations will, therefore, correlate at that time lag, cancelling the random time of each source pulse by a relative measurement. This is the principle that all of the above methods make use of. In addition, correlation is an effective tool to suppress noise, or any signal that is incoherent at the particular time lags assumed. Each differential traveltime extracted by correlation of tremor recorded at two stations does not constrain the source location to a point alone. Instead, it will constrain the location to lie anywhere on a hyperbolic surface corresponding to equal differential distance from the source when assuming a uniform velocity. For propagation in 2-D the solutions are constrained to lie on a hyperbolic curve. With information from two station pairs, we have two hyperbolic curves, which will usually intersect at one or two points. Further station pairs will provide redundant constraints on a single-point solution. Thus, we need differential traveltime measurements from at least two station pairs to locate a source in 2-D. Similarly, we need three pairs to locate a source in 3-D. When we back project correlation functions to hypothetical source locations, again using uniform velocity, the energy distribution for each of them will peak on a hyperbolic surface. Stacking two correlation functions will usually produce a peak at a point, while stacking more station pairs will enhance that peak relative to other parts of the hyperbolae. However, this suppression is rather weak and is fundamentally limited by the number of available station pairs. The method we present below is intended to enhance signal by using a modified form of correlation, which we refer to as double correlation. The double correlation calculates the cross correlation of correlations from two pairs of stations. It uses the fact that when the signals are temporarily aligned for a given source, they should correlate not only between stations, but also between cross-correlation functions. The double correlation will suppress noise correlated at time shifts that are inconsistent with the assumed velocity model. It will also suppress random noise in the signals. Stehly et al. (2008) used a similar approach to reconstruct the Green's function between two stations by correlating the coda of correlations of ambient noise. Various assumptions are implicit in all correlation approaches, including that a source is not significantly autocorrelated and that it radiates similar time series in all directions. Besides, we need to know the velocity well, which we may not. The tremor signals may contain both body and surface waves, so it becomes ambiguous which velocity to use and whether we need to solve a 2-D or a 3-D problem. The Green's functions between the source and different receivers may be complex and differ by more than a time delay from one station to the other (multipathing, phase conversion). In that case the covariograms may correlate poorly. However, variations in the tremor source activity level should still produce similar variations at each station pair in the more slowly-varying covariogram envelopes. Use of envelopes may help to stabilize the back projection at a loss of spatial resolution. METHOD Our method of location extends the previous correlation methods by applying double correlation instead of a single correlation. We make the same assumptions as are normally made in, for example, single-correlation methods, including that when multiple sources are present, they are not correlated with each other. The double-correlation location method can be divided into three aspects: (1) the computation of the double correlations of triplets of seismographs, (2) the back projection of the stacked double correlation to hypothetical sources in a geographic grid and (3) the selection of station triplets. In the following subsections we address these aspects separately. Double correlation Here we introduce the procedure to compute a double correlation of a station triplet. The correlation at a pair of stations, a and b, is \begin{eqnarray} C_{ab}(j) = \sum \limits _{i=1}^{N-j} a(i) b(i+j) , \end{eqnarray} (1) where i and j are sample indices, a(i) and b(i) are the seismograms at stations a and b, respectively, and N is the total number of sample points. Having computed a second such correlation for stations c and d \begin{eqnarray} C_{cd}(j) = \sum \limits _{i=1}^{N-j} c(i) d(i+j) , \end{eqnarray} (2) we can then define a double correlation as \begin{eqnarray} C_{abcd}(m) = \sum \limits _{j=-M}^{M} C_{ab}(j) C_{cd}(j+m) , \end{eqnarray} (3) where M = N − 1 is the maximum sample lag in the first cross correlation. The double correlation is useful for two purposes. First, the level of redundancy is increased over a single correlation, which suppresses random noise. Second, this double correlation will peak at the double differential traveltime between the two pairs. Back projecting this function to the source in 2-D will concentrate energy around a point as opposed to the hyperbola obtained by the back projection of a single correlation. However, we do not implement double correlation as described above. We note that the sum in eq. (3) will not have much of a contribution from time lags beyond a finite limit corresponding to the possible differential traveltimes of sources within the network. Therefore, we implement double correlation in the following manner. First, we convert our time series to an analytic signal \begin{eqnarray} \skew4\hat{f}(t) = f(t) + i \mathcal {H}[f(t)] , \end{eqnarray} (4) where f(t) is the original time series, |$\mathcal {H}[...]$| denotes the Hilbert transform and |$\skew4\hat{f}(t)$| is the analytic signal of f(t). Then, we divide the time span of our signals, T, into K equal sub intervals of length ΔT = T/K. Next, we select a triplet of stations, a, b and c. We choose one of those as a reference station, for example a, and correlate stations a and b on one hand and stations a and c on the other. Choosing this particular way of combining stations is not fundamentally required for our approach to work, but for our data character and geometry, this approach is effective. The issue related to the selection of station triplets will be discussed later. For the kth interval the correlations at stations a and b, and stations a and c are \begin{eqnarray} {}_k\hat{C}_{ab}(j) = \sum \limits _{i=k\Delta T+1}^{(k+1)\Delta T} \hat{a}(i) \hat{b}^{\ast }(i+j) \end{eqnarray} (5) and \begin{eqnarray} {}_k\hat{C}_{ac}(j) = \sum \limits _{i=k\Delta T+1}^{(k+1)\Delta T} \hat{a}(i) \hat{c}^{\ast }(i+j), \end{eqnarray} (6) respectively, where the time series are all complex and * indicates the complex conjugate. Finally, we correlate those two \begin{eqnarray} \hat{C}_{abc}(m) = \sum \limits _{k=0}^{K-1} {}_k\hat{C}_{ab}(i) \thinspace \thinspace {}_k\hat{C}_{ac}^{\ast }(j+m) , \end{eqnarray} (7) where i and j are the expected sample lags calculated from the theoretical differential traveltimes between stations a and b, and a and c, respectively. These theoretical differential traveltimes are obtained from, for example, ray tracing in a pre-defined velocity model. The amplitude of this complex double correlation |$\hat{C}_{abc}$| is a measure of the double correlation at a given differential time shift. It will peak at the double differential traveltime between the two station pairs. The phase of |$\hat{C}_{abc}$| measures the precise phase lag or time lag between the time series. In order to make use of the phase information in |$\hat{C}_{abc}$|⁠, a complete velocity model, which we may not know precisely, is needed to calculate the differential traveltimes between stations. We, therefore, focus our analysis on the amplitude of |$\hat{C}_{abc}$|⁠. The double correlation will suppress random noise more than single correlation and suppress non-stationary correlated noise. Stacking the back-projected double correlations for all triplets will further suppress correlated peaks that are not common for all triplets, for example, components of the Green's functions between source and stations that vary from one path to another. A further advantage of this double correlation is that each back projection provides a localized estimate of the source location (in 2-D). Therefore, when we back project to a hypothetical source and stack different such double correlations, we stack equivalent estimates of the source location. However, stacking single correlations involves stacking hyperbolic distributions that have very different localization properties. Back projection The back projection of the above double correlation is done in the following manner: For each hypothetical source location the travel time to any given station is computed using a velocity model. The time lags between stations a and b and stations a and c are thus estimated and their double difference computed. We calculate the modulus of the corresponding double correlation at that differential lag. Finally, the back projection is defined as the stack of double-correlation moduli over all station triplets. It is common that amplitude variations of seismograms are not compensated in the back projection (e.g. Shapiro et al. 2006; Haney 2014). The amplitude variation depends on a number of effects, for example, focusing effects due to velocity variations, attenuation, etc. We argue that for Katla volcano, we cannot know the amplitude variation sufficiently well since we are not able to model these effects precisely. Besides, the amplitude variation is not important, because we are exploiting coherent phase information in our method. The amplitude simply affects the relative weighting of the different covariograms or double covariograms. To optimize the use of the phase information, it is more appropriate to select such weighting based on some measure of correlation uncertainties in each case. However, such uncertainties may be difficult to estimate reliably and will vary from one data set to another. Selection of station triplets In the calculation of double correlations we choose to combine stations according to a reference-station scheme, rather than including all possible combinations of stations. We have tested this synthetically, also using the same real-data examples as presented later. We find that this reference-station scheme works better than including all possibilities, partly because by restricting the double correlations to correlations within closed triangles reduces the occurrence of nearly parallel station pairs. It, in turn, alleviates overlapping of tangential hyperbolae. As a result, hyperbolic streaking is reduced in the result. Under the reference-station scheme, a network of n stations will have |$C^n_3$| combinations of three stations and for each, three possible choices of a reference station. The order of the remaining two stations in each triplet is arbitrary and the two choices are completely redundant. There are, therefore, N independent triplets, where N is \begin{eqnarray} N = 3 C^n_3 . \end{eqnarray} (8) For a station network with n = 10 seismographs the number of triplets is N = 360. APPLICATION TO REAL DATA We apply the double-correlation method to tremor at Katla volcano, southern Iceland, during unrest in July 2011. Katla volcano is located near the central south coast of Iceland near the tip of the southward propagating Eastern Volcanic Zone (Fig. S1 in the Supporting Information). Katla is a subglacial volcano with a well-developed caldera structure (about 10 km in diameter). The geology of the area is complex, indicating that there could be strong attenuation and focusing effects. Both the ice cover and the caldera rim are outlined in Figs 1 and 2 for reference. The unrest of the volcano started in August 2010 and an increase of seismicity in the area was observed in July 2011. The tremor started at ∼19:00 GMT on 2011 July 8 and lasted for 23 hr (see Sgattoni (2016) for a detailed description of this unrest). Three main phases of tremor were observed (Sgattoni 2016). Two were related to the cauldrons (depressions on the glacier surface, Fig. 1) that collapsed during the tremor episode (Gudmundsson & Sólnes 2013). It is unclear if the cauldrons collapsed due to a hydrothermal process or a small, subglacial magmatic eruption (Sgattoni 2016). The third phase was associated with a glacial flood to the south east of the caldera. This phase only lasted for an hour and was only clearly identified at the station closest to the flood (Fig. S2 in the Supporting Information). Including the permanent stations from the Icelandic Meteorological Office and three temporary stations deployed by Uppsala University and Reykjavik University on nunataks (rock protruding the ice surface) in the ice, a total of 12 nearby stations were in operation during the tremor period. Here we use the 10 stations closest to the caldera for our analysis (Fig. S2 in the Supporting Information). Figure 1. Open in new tabDownload slide (a,b) The back projection of single correlations for a tremor at Katla volcano, Iceland in July 2011. Two strategies of amplitude normalization are used. From left to right: without any normalization and one-bit normalization of seismograms. (c,d) The same as panels (a) and (b) except showing the back projection of double correlations. The square root of back-projected double correlations is taken to allow a fair comparison to the single correlations. The colours define normalized energy, dark red for maximum energy. Red inverted triangles: seismic stations. Black solid line: glacier. Black dashed line: caldera outline. Black dots: locations of earthquakes occurred during the tremor period. White open circles: locations of cauldrons. Figure 1. Open in new tabDownload slide (a,b) The back projection of single correlations for a tremor at Katla volcano, Iceland in July 2011. Two strategies of amplitude normalization are used. From left to right: without any normalization and one-bit normalization of seismograms. (c,d) The same as panels (a) and (b) except showing the back projection of double correlations. The square root of back-projected double correlations is taken to allow a fair comparison to the single correlations. The colours define normalized energy, dark red for maximum energy. Red inverted triangles: seismic stations. Black solid line: glacier. Black dashed line: caldera outline. Black dots: locations of earthquakes occurred during the tremor period. White open circles: locations of cauldrons. First, local earthquakes are removed from the seismograms by automated clipping and manual removal (Fig. 1). Such events do not invalidate our approach, but they may not be at the same location as the tremor source. As the events are strong signals, they may thus inject strong ‘correlated noise’ to the analysis of tremor. Data are then filtered between 0.8 and 1.5 Hz using a fourth-order Butterworth filter. We choose this low-frequency band because the energy is more stable there than at higher frequencies (Sgattoni 2016). We expect that the tremor sources are located inside the station network and the amplitude of some of the closer stations may dominate. It is, therefore, sensible to normalize amplitudes in some way. We apply one-bit normalization to each seismogram before computing the covariograms and compare resulting back-projected maps with those obtained from raw covariograms (without any normalization of traces). The choices of normalization are somewhat subjective. To optimally weight the covariograms in the back-projected stack, we need to know their uncertainties, which are difficult to estimate. We assume that the tremor contains a strong surface-wave component (Sgattoni 2016). It has been shown that volcanic tremor often consists mainly of surface waves (Ferrazzini et al. 1991; Wegler & Seidl 1997). We, therefore, focus on the surface-wave component and perform the back projection in a 2-D grid with a resolution of 0.5 km in each direction. To choose velocity for the back projection, we generate back-projected maps with several feasible velocities based on the tomographic results of Jeddi (2016). We choose that velocity which gives a local minimum of the entropy (Li et al. 2016). This entropy quantifies the diffuseness of a distribution and a small entropy implies localized distribution. We find that a velocity of 1.2 km s−1 gives the smallest entropy and most focused energy distribution. Since we are working with the modulus of the complex double correlation, which is related to the surface-wave envelopes (energy), this velocity should correspond to group velocity. This velocity may appear low. However, it is comparable to the results of Haney et al. (2011) at Hekla volcano and Obermann et al. (2016) at Snæfellsjökull volcano, as well as Jeddi (2016) at Katla volcano and Benediktsdóttir et al. (private communication, 2016) at Eyjafjallajökull, all in Iceland, at similar frequencies (1 Hz). In the case when no trace normalization is applied the energy of the back-projected single and double correlations peaks along hyperbolae corresponding to constant time shifts for the three aligned stations along the northeastern caldera rim (Fig. 1). This energy is significantly better focused in the double correlation than the single correlation results. However, an elongated streak of energy persists to the south of the caldera. The background energy is also lowered in relative terms by double correlation compared to the single correlation, demonstrating the effect of noise suppression. In the case where we correlate one-bit normalized traces the hyperbolic streaking artefact is significantly suppressed. Again, the double correlation suppresses these artefacts more effectively than the single correlation. We are not able to evaluate which normalization scheme is optimal unless we know the characteristics of noise processes in the signals. In our particular examples (real and synthetic data, see also Fig. 2) using one-bit normalized traces appears to be the better choice. Correlating traces without any normalization may be suboptimal because the high-amplitude stations near the source dominate and we, therefore, do not take full advantage of all possible geometries offered by the network. Figure 2. Open in new tabDownload slide (a,b) The back projection of single correlations for synthetic examples. Two strategies of amplitude normalization are used. From left to right: without any normalization and one-bit normalization of seismograms. (c,d) The same as panels (a) and (b) except showing the back projection of double correlations, again the square root is taken for the double correlations. For better comparison the same colour scale is used as the real-data example. The symbols follow Fig. 1, except that the black cross denotes the source location. Figure 2. Open in new tabDownload slide (a,b) The back projection of single correlations for synthetic examples. Two strategies of amplitude normalization are used. From left to right: without any normalization and one-bit normalization of seismograms. (c,d) The same as panels (a) and (b) except showing the back projection of double correlations, again the square root is taken for the double correlations. For better comparison the same colour scale is used as the real-data example. The symbols follow Fig. 1, except that the black cross denotes the source location. In both cases, the maximum back-projected energy lies close to the location of two closely spaced cauldrons that collapsed in the southern part of the caldera. This lends credibility to the location method, because it is reasonable to expect tremor in association with volcanic, geothermal or glacio-hydrological processes near cauldrons. The third cauldron appears to be associated with weaker tremor (Sgattoni 2016). SYNTHETIC TESTS We have tested the double correlation approach with various synthetic tests. When using perfect data in a uniform velocity structure, the single source is well recovered with a resolution limitation controlled by the frequency content of the data. With significant velocity heterogeneity (compared to the finite frequency limitation) we lose our ability to predict travel time precisely. This results in erratic phase behaviour in the complex covariograms. However, the modulus of the back-projected double correlation remains stable, but its finite width, controlled by the frequency band-width of the data, translates into a broader location estimate. Our tests of the method include both various noise processes and propagation effects. We have added random incoherent noise to the time series, coherent noise in the form of the source impulses delayed according to a uniform body-wave velocity of 2.7 km s−1 to simulate the presence of body waves in the tremor signal, and signals corresponding to plane waves from a range of azimuths to simulate the effect of distant sources, for example, microseisms. We have also tested adding a set of random scatterers given a random strength, a random scattering width and a random orientation, and modelling the effect of heterogeneity by delaying the source impulses according to a 2-D heterogeneous velocity model for surface waves and given those random amplitude perturbations on top of simple geometrical spreading. The velocity model (Fig. S3 in the Supporting Information) is generated by adding a perturbation, which is described by a Gaussian autocorrelation function, to a homogeneous velocity model with a velocity of 1.2 km s−1 (Klimeš 2002). The tests are done with seismograms filtered in the same frequency range as the Katla tremor data (0.8–1.5 Hz). The back projection is done with the same uniform velocity as chosen for the Katla tremor data, which is also the average group velocity in the synthetic velocity structure. We have done a number of tests with different strength of the described effects added. These tests reveal that, even with a strong velocity variation of ±20 per cent (a standard deviation of 10 per cent) with a correlation length of 7 km, we still recover the source location to within a few km. Suppression of incoherent noise is very effective and the point-source location is well recovered even if this noise has the same amplitude as the signal on average (signal-to-noise ratio is 1). The body-wave component of the synthetic seismograms is coherent at time delays corresponding to distances and velocities that are very different from those for the surface waves and does not significantly affect the recovery of the point source even if given a relative amplitude as big as the surface-wave component. Despite the appearance of some hyperbolic artefacts, distant plane-wave sources have little effect on the source location estimate unless they dominate over the local source in amplitude. The most significant synthetic components are the addition of scatterers and random amplitude variations. With 70 scatterers randomly distributed within the caldera, given a random orientation and a Gaussian width of ±60° (one standard deviation) and a random scattering strength between 0 and 1, we simulate the appearance of the stacked, back-projected energy distribution obtained with the real data. By enhancing the amplitudes of synthetic recordings at stations on nunataks in the ice cap we simulate the effects of different normalization strategies on the back-projected distribution. The structure in Katla volcano is very complex so that the site magnification and attenuation effects could be strong. This would result in dramatic variations in amplitudes between different seismograms. In the testing we have also tried to simulate the features appeared in the back-projection results of the real-data (Fig. 2). The simulations include both single and double correlation, again using traces with and without one-bit normalization. In all cases we use synthetic data containing velocity heterogeneity (±20 per cent velocity variation, 7 km correlation length), added white noise (signal-to-noise ratio of 1 on average), a body-wave component equal to the surface wave component in amplitude and a plane-wave component with a relative amplitude of 0.1 on average. We have added 70 random scatterers with the same geometrical characteristics as described above. In Figs 2(a) and (c), we show the back-projected energy without any amplitude normalization. The source, which is located at coordinates (2, −4), is clearly resolved (Fig. 2). These synthetic tests show similar behaviour as the real-data examples in which the strong spurious peaks are suppressed in the double correlation. The energy is also more focused by double correlation than single correlation. Note the common structure in the peripheral distribution, which is controlled by the station geometry. One-bit normalization of traces reduces the streaking hyperbolic artefacts (Figs 2b and d). The energy is more focused around the central source region in the double- than the single correlation. DISCUSSION AND CONCLUSIONS We introduce a double-correlation method to locate tremor sources. This method back projects stacks of complex, doubly-correlated tremor records to hypothetical sources in a geographic grid for location estimates. Compared to the existing single-correlation method, the double correlation is better in suppressing both random noise and noise correlated at time shifts that are inconsistent with the assumed velocity. The back projection of each double correlation provides a localized estimate of source locations. Stacking different back-projected maps will, therefore, give equivalent estimates of the source locations, as opposed to the single correlation, where the individual maps contain strong signatures of the station geometry. Application of the method to real data from Katla volcano, Iceland reveals a stable energy source during the tremor period. The location of that source lies close to the cauldrons collapsed during that period. In the synthetic tests we are able to recover the primary source location with an accuracy of a few km despite various types of noise added to the seismograms. Although we have added significant secondary sources, such as scatterers and plane-wave sources, these are well suppressed by the double correlation. Unless they are comparable in amplitude to the primary source, they will not have a significant effect on its source-location estimate. Random noise is effectively suppressed by the double correlation. Correlated noise, such as the body-wave component, does not have a significant effect on the result because it has a time delay which is very different from the surface-wave component. Velocity perturbation has a little effect on the location estimate except broadening it to have a resolution of a few km, depending on the frequency content of the data. We are able to simulate the real-data example qualitatively with the added effects to the synthetic data. The energy distributions in the back-projected maps are similar for both data sets. In particular, hyperbolic streaking artefacts in the real data are well simulated. We, therefore, conclude that the peak in the back-projected double correlation of Katla's data (Fig. 1(d)) reveals the location of the dominating energy source during the tremor period. We need not assume that the tremor is dominantly composed of surface waves for a 2-D back-projection to work as at surface-wave velocities the body-wave components do not stack coherently. It suffices that the tremor contains a significant surface-wave component. Focusing on surface waves limits the analysis to 2-D, so source depth cannot be recovered. For a source with a finite depth, a phase shift will be introduced into the surface waves, but this will be common to all observations and, therefore, cancel in the double-differential times of the double correlation. For the higher frequencies (1.5–4 and 4–9 Hz) of the Katla data we have not been able to recover a coherent picture of the source with this method. We suspect that is caused by increasing vigour of scattering effects at increasing frequency in tune with the fact that addition of discrete scatterers for the 0.8–1.5 Hz range proved to be a significant obscuring effect in our synthetic examples. In this short article, we have not discussed all possible complications, such as that for some station geometries with multiple dominating sources some spurious correlation may be incompletely suppressed even after double correlation. Such problems may be alleviated by further processing steps, for example, with higher-order cross correlations. The specific method presented here is only one of a family of methods based on higher-level cross correlation. Slightly different approaches may be better suited depending on the characteristics of the data. However, the suggested double-correlation strategy does suppress some types of uncorrelated and correlated noise better than other alternatives and may be useful for a spectrum of array-processing problems, not only tremor analysis. Acknowledgments We would like to thank two anonymous reviewers for their constructive comments. This work was partially supported by the Standup for Energy national research program. REFERENCES Almendros J. , 2003 . Performance of the radial semblance method for the location of very long period volcanic signals , Bull. seism. Soc. Am. , 93 ( 5 ), 1890 – 1903 . Google Scholar Crossref Search ADS WorldCat Ballmer S. , Wolfe C.J., Okubo P.G., Haney M.M., Thurber C.H., 2013 . Ambient seismic noise interferometry in Hawai’i reveals long-range observability of volcanic tremor , Geophys. J. Int. , 194 ( 1 ), 512 – 523 . Google Scholar Crossref Search ADS WorldCat Battaglia J. , Aki K., 2003 . Location of seismic events and eruptive fissures on the Piton de la Fournaise volcano using seismic amplitudes , J. geophys. Res. , 108 ( B8 ), 2364 , doi:10.1029/2002JB002193 . Google Scholar Crossref Search ADS WorldCat Capon J. , 1973 . Signal processing and frequency-wavenumber spectrum analysis for a large aperture seismic array , Methods Comput. Phys. , 13 , 1 – 59 . Google Scholar Crossref Search ADS WorldCat Di Grazia G. , Falsaperla S., Langer H., 2006 . Volcanic tremor location during the 2004 Mount Etna lava effusion , Geophys. Res. Lett. , 33 ( 4 ), L04304 , doi:10.1029/2005GL025177 . Google Scholar Crossref Search ADS WorldCat Droznin D. , Shapiro N., Droznina S.Y., Senyukov S., Chebrov V., Gordeev E., 2015 . Detecting and locating volcanic tremors on the Klyuchevskoy group of volcanoes (Kamchatka) based on correlations of continuous seismic records , Geophys. J. Int. , 203 ( 2 ), 1001 – 1010 . Google Scholar Crossref Search ADS WorldCat Ferrazzini V. , Aki K., Chouet B., 1991 . Characteristics of seismic waves composing Hawaiian volcanic tremor and gas-piston events observed by a near-source array , J. geophys. Res. , 96 ( B4 ), 6199 – 6209 . Google Scholar Crossref Search ADS WorldCat Furumoto M. , Kunitomo T., Inoue H., Yamaoka K., 1992 . Seismic Image of the Volcanic Tremor Source at Izu-Oshima Volcano, Japan , in Volcanic Seismology , pp. 201 – 211 , eds Gasparini P., Scarpa R., Aki K., Springer . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Gudmundsson M. , Sólnes J., 2013 . Activity at Katla and jökulhlaup in Múlakvísl river 2011 (in Icelandic) , in Natural Hazard in Iceland: Volcanic Eruptions and Earthquakes (in Icelandic) , pp. 228 – 229 , eds Sólnes J., Sigmundsson F., Bessason B., University of Iceland Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Gudmundsson O. , Brandsdottir B., 2010 . Geothermal noise at Olkelduhals, SW Iceland , Jokull , 60 , 89 – 102 . OpenURL Placeholder Text WorldCat Haney M.M. , 2010 . Location and mechanism of very long period tremor during the 2008 eruption of Okmok Volcano from interstation arrival times , J. geophys. Res. , 115 ( B10 ), doi:10.1029/2010JB007440 . OpenURL Placeholder Text WorldCat Haney M.M. , 2014 . Backprojection of volcanic tremor , Geophys. Res. Lett. , 41 ( 6 ), 1923 – 1928 . Google Scholar Crossref Search ADS WorldCat Haney M.M. , Nies A., Masterlark T., Needy S., Pedersen R., 2011 . Interpretation of Rayleigh-wave ellipticity observed with multicomponent passive seismic interferometry at Hekla Volcano, Iceland , Leading Edge , 30 ( 5 ), 526 – 531 . Google Scholar Crossref Search ADS WorldCat Jeddi Z. , 2016 . Seismological investigation of Katla volcanic system (Iceland): 3D velocity structure and overall seismicity pattern , PhD thesis , Uppsala University . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Klimeš L. , 2002 . Correlation functions of random media , Pure appl. geophys. , 159 ( 7 ), 1811 – 1831 . OpenURL Placeholder Text WorldCat Konstantinou K.I. , Schlindwein V., 2003 . Nature, wavefield properties and source mechanism of volcanic tremor: a review , J. Volcanol. Geotherm. Res. , 119 ( 1–4 ), 161 – 187 . Google Scholar Crossref Search ADS WorldCat Li K.L. , Gudmundsson Ó., Tryggvason A., Bödvarsson R., Brandsdóttir B., 2016 . Focusing patterns of seismicity with relocation and collapsing , J. Seismol. , 20 ( 3 ), 771 – 786 . Google Scholar Crossref Search ADS WorldCat Obermann A. , Lupi M., Mordret A., Jakobsdóttir S.S., Miller S.A., 2016 . 3D-ambient noise Rayleigh wave tomography of Snæfellsjökull volcano, Iceland , J. Volcanol. Geotherm. Res. , 317 , 42 – 52 . Google Scholar Crossref Search ADS WorldCat Rost S. , Thomas C., 2002 . Array seismology: methods and applications , Rev. Geophys. , 40 ( 3 ), 2-1 – 2-27 . Google Scholar Crossref Search ADS WorldCat Sgattoni G. , 2016 . Characteristics and geological origin of earthquakes and tremor at Katla Volcano (S-Iceland) , PhD thesis , Alma Mater Studiorum - Università di Bologna . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Shapiro N.M. , Ritzwoller M.H., Bensen G.D., 2006 . Source location of the 26 sec microseism from cross-correlations of ambient seismic noise , Geophys. Res. Lett. , 33 ( 18 ), L18310 , doi:10.1029/2006GL027010 . Google Scholar Crossref Search ADS WorldCat Stehly L. , Campillo M., Froment B., Weaver R.L., 2008 . Reconstructing Green’s function by correlation of the coda of the correlation (C3) of ambient seismic noise , J. geophys. Res. , 113 ( B11 ), B11306 , doi:10.1029/2008JB005693 . Google Scholar Crossref Search ADS WorldCat Wegler U. , Seidl D., 1997 . Kinematic parameters of the tremor wave field at Mt. Etna (Sicily) , Geophys. Res. Lett. , 24 ( 7 ), 759 – 762 . Google Scholar Crossref Search ADS WorldCat Zeng X. , Ni S., 2010 . A persistent localized microseismic source near the Kyushu Island, Japan , Geophys. Res. Lett. , 37 ( 24 ), L24307 , doi:10.1029/2010GL045774 . Google Scholar Crossref Search ADS WorldCat SUPPORTING INFORMATION Supplementary data are available at GJIRAS online. Figure S1. A map of Katla volcano. Red inverted triangles show the locations of seismic stations. The white areas mark glaciers. Figure S2. Seismograms from all stations in the network. The amplitudes are normalized by a common constant for all stations. Note the differences in scale for the normalized amplitudes. The tremor was strongest between the 4th and the 10th hour. A number of short bursts occurred during this period. A distinct phase was observed at HVO around the 9th hour. This phase occurred during the glacial flood close to the station. Figure S3. The 2-D velocity model used in the simulation of Fig. 2 in the paper. The velocity perturbation is described by a Gaussian autocorrelation function with a correlation length of 7 km and a standard deviation of 10 per cent, corresponding to ±20 per cent maximum velocity variation. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Authors 2016. Published by Oxford University Press on behalf of The Royal Astronomical Society. TI - A double-correlation tremor-location method JO - Geophysical Journal International DO - 10.1093/gji/ggw453 DA - 2017-02-01 UR - https://www.deepdyve.com/lp/oxford-university-press/a-double-correlation-tremor-location-method-VH2kWfx7eY SP - 1231 VL - 208 IS - 2 DP - DeepDyve ER -