## 1 Introduction

It has been shown that correlating long time series of seismic ambient noise recorded at two stations makes it possible to retrieve the Earth response (i.e. the Green's function) between these two stations (e.g., Campillo, 2006; Gouédard et al., 2008; Larose et al., 2006; Sabra et al., 2005a; Shapiro and Campillo, 2004). This technique has provided useful results for studies of tomographic imaging (e.g., Sabra et al., 2005b; Shapiro et al., 2005; Yang et al., 2007; Yao et al., 2006, 2008) and monitoring (Brenguier et al., 2008a, 2008b; Sens-Schönfelder and Wegler, 2006; Wegler and Sens-Schönfelder, 2007).

In practice, since the ambient seismic noise is dominated by the surface-wave contributions, only the direct arrivals of surface waves are clearly visible in the noise-correlation signals. However, theory says that in the case of a perfectly isotropic noise field, the noise correlation $C\left({x}_{A},{x}_{B},\tau \right)$ provides the complete Green's function $G\left({x}_{A},{x}_{B},\tau \right)$, including all types of waves that propagate between the two stations in x_{A} and x_{B} (e.g., Gouédard et al., 2008; Lobkis and Weaver, 2001; Roux et al., 2005; Sánchez-Sesma et al., 2006; Snieder, 2007; Wapenaar, 2004). If we let u(x_{A},t) and u(x_{B},t) denote time-dependent random wavefields recorded at two sensors, A and B, we can write:

$$C\left({x}_{A},{x}_{B},\tau \right)={\u2329u\left({x}_{A},t\right),u\left({x}_{B},t\right)\u232a}_{T}$$ | (1) |

_{∞}as the correlation when T→ ∞. A theoretical guess for the reconstruction of the Earth response is:

$$\frac{\partial}{\partial \tau}{C}_{\infty}\left({x}_{A},{x}_{B},\tau \right)\propto {G}^{+}\left({x}_{A},{x}_{B},\tau \right)-{G}^{-}\left({x}_{A},{x}_{B},-\tau \right)$$ | (2) |

^{+}and G

^{−}represent the causal and anticausal Green's function, respectively. Note that in practical applications, most studies choose to deal with a normalized correlation function:

$$C\left({x}_{A},{x}_{B},\tau \right)={\u2329u\left({x}_{A},t\right),u\left({x}_{B},t\right)\u232a}_{T}^{norm}=\frac{{\u2329u\left({x}_{A},t\right),u\left({x}_{B},t\right)\u232a}_{T}}{\sqrt{C\left({x}_{A},{x}_{A},0\right)C\left({x}_{B},{x}_{B},0\right)}}$$ | (3) |

_{A}and x

_{B}. In this case, the maximum of the normalized correlation function corresponds to the coherence level between the two stations.

In the ideal case, the later part of the noise correlation should contain the coda part of the Green's function, and it might be possible to correlate these waves as we correlate standard earthquake codas, to reconstruct the Green's function (Campillo and Paul, 2003). This is the concept of the so-called C^{3} method, i.e. re-correlating the coda of the noise-correlation functions to reconstruct the Green's function (Stehly et al., 2008).

For microseism records dominated by surface noise, the required spatial isotropy of the recorded wavefield is mainly produced by the distribution of sources and by the scattering of seismic waves in the medium. However, usually the distribution of the noise sources does not provide a perfect random noise field. In this case, the theoretical requirements are not completely fulfilled, and the noise correlation does not perfectly reconstruct the surface-wave Green's function, as some fluctuations remain in the correlation results. This, in practice, prevents the identification of other amplitude arrivals that are smaller than the predominant fundamental mode of the surface waves. This is the main limitation in the implementation of the C^{3} method, i.e. does the later part of the noise correlation contain the coda part of the surface-wave Green's function? Note that this is also a central issue for noise-based temporal monitoring of seismic velocity, for which measurements are performed on noise correlation codas.

In their recent study, Stehly et al. (2008) demonstrated successful results for the C^{3} function, clearly showing the direct surface-wave part of the Green's function. The success of this technique is an indication that the later part of the ambient noise correlations is meaningful and contains multiply scattered waves with properties of regular earthquake codas. Stehly et al. (2008) showed that the C^{3} method helps to improve the time symmetry of correlation functions, thereby increasing the quality of the data used in noise-based tomographic studies. Indeed, since the noise correlation functions should converge to the causal and anticausal Green's function (Eq. (2)), time symmetry in noise correlations may constitute a way of controlling the quality of the Green's function reconstruction from the correlation process. In particular, travel-time estimations must be the same in both causal and anticausal parts. That is why some tomographic studies use only time-symmetric noise correlations with good signal-to-noise ratios. Since Stehly et al. (2008), some studies have proposed a formalization of the C^{3} method (de Ridder et al., 2009; Garnier and Papanicolaou, 2009). In particular, Garnier and Papanicolaou (2009) proved the validity of this method based upon the stationary-phase analysis of the C^{3} function leading terms, and they confirmed that this method can enhance the quality of travel-time estimates in the case of anisotropic source distributions.

As indicated above, the noise-correlation technique has provided useful information for seismological applications, and this in spite of the imperfect reconstruction of the Green's function in correlation functions. In particular, an anisotropic source distribution leads to biased travel-times extracted from correlation functions (Froment et al., 2010; Weaver et al., 2009). One of the main issues in the improvement of noise-based measurements is thus to understand and eliminate these source effects in correlation results. In this study, we show the improvement of noise-based results in seismological applications through the use of the C^{3} method.

In Section 2, we present the data and the processing used in this study to compute the C^{3} function. In Section 3, we examine the time symmetry of the C^{3} function and its relationship with the station distribution in the network. Then, we discuss the interest of combining the information from the different correlation functions. Finally, we show the correlation function that is obtained by iterating the correlation process on the C^{3} coda. In particular, we investigate the evolution of the reconstruction of the Green's function throughout the iteration process.

## 2 Data and C^{3} computation

Throughout this study, we have used pairs of stations that are located in the Alps as part of a regional network of 150 broad-band stations. We consider the vertical motions of seismic ambient noise recorded continuously over one year. The noise processing is the same as in Stehly et al. (2008).

For each station pair, we correlate the noise records, which are pre-filtered in the 0.1–0.2-Hz frequency band, and the noise correlation (C^{1}) is obtained by stacking all of the day-by-day correlations. We then apply the C^{3} method between two stations A and B, according to the following successive steps.

1. We take a third station, S, in the network located at x_{S}, and we consider the result of the A-S and B-S noise correlations (C^{1}) as the signal recorded at A and B, respectively, if a source was present at x_{S}. The receiver in S thus has the role of a virtual source.

2. We then select a 1200-s-long time window (T_{coda}) in the later part of the noise correlation (C^{1}), beginning at twice the Rayleigh wave travel-time (i.e., in the same way as Stehly et al. (2008)). Note that the coda in both the positive and negative parts of the noise correlation (C^{1+} and C^{1−}) are selected.

3. We cross-correlate the codas in the A-S and B-S noise correlations. Note that we whiten the codas between 0.1 and 0.2 Hz before correlation; the motivation for this processing is discussed later. We compute the correlation functions between the positive-time and negative-time codas to finally obtain two correlations that we denote as C^{3++}_{S} and C^{3−−}_{S}, respectively.

4. We average over the N stations of the seismic network to obtain C^{3++} and C^{3−−}:

$${C}^{3++}\left({x}_{A},{x}_{B},{\tau}^{\prime}\right)=\frac{1}{N}\sum _{j=1}^{N}\left[{\u2329{C}^{1+}\left({x}_{A},{x}_{Sj},\tau \right),{C}^{1+}\left({x}_{B},{x}_{Sj},\tau \right)\u232a}_{{T}_{coda}}^{norm}\right]$$ | (4) |

$${C}^{3--}\left({x}_{A},{x}_{B},{\tau}^{\prime}\right)=\frac{1}{N}\sum _{j=1}^{N}\left[{\u2329{C}^{1-}\left({x}_{A},{x}_{Sj},\tau \right),{C}^{1-}\left({x}_{B},{x}_{Sj},\tau \right)\u232a}_{{T}_{coda}}^{norm}\right]$$ | (5) |

5. We stack C^{3++} and C^{3−−} to obtain the actual C^{3} correlation function between A and B:

$${C}^{3}\left({\tau}^{\prime}\right)=\frac{1}{2}\left({C}^{3++}\left({\tau}^{\prime}\right)+{C}^{3--}\left({\tau}^{\prime}\right)\right)$$ | (6) |

The steps described above highlight the difference in the correlation technique between noise correlations (C^{1}) and the C^{3} method, that points out advantages of this iterative method. For instance, since the C^{3} method correlates signal that has been extracted from the C^{1} computation, it is expected to be more coherent than seismic ambient noise. Furthermore, we can distinguish coda and ballistic contributions in C^{1} functions (which is not possible in seismic noise); this provides scattering averaging in the correlation process when only coda waves are selected. Finally, we point out that, by construction, the virtual sources S are uncorrelated in the C^{3} method, which may not be the case for noise sources in C^{1} computation.

Fig. 1 illustrates the C^{1} and C^{3} correlation functions for one station pair. The good agreement between the pulse shapes in the C^{3} and C^{1} correlation functions shows that the C^{3} correlation process has extracted from the coda of the C^{1} correlation the necessary phase information to reconstruct the direct surface wave part of the Green's function.

## 3 Network stations as virtual sources

### 3.1 Source signature in correlation functions

In the case of noise sources that are homogeneously distributed, the positive and negative parts of the noise-correlation functions are expected to be symmetric, as they represent the causal and anticausal parts of the Green's function, respectively. In practice, these two parts may not have the same amplitude, which reflects a difference in the energy flux that propagates between the two stations in both directions (Paul et al., 2005; van Tiggelen, 2003). Thus, if the source density (or virtual source power) is larger for one direction than for the other, a time asymmetry will be seen for the correlation function (Fig. 2). This means that the level of symmetry in the correlations, up to first order, is an indication of the source distribution.

What makes C^{3} appealing is that, in theory, the “noise” sources that contribute to C^{1} and C^{3} are different. Considering the two stations A and B, the noise correlation (C^{1}) is directly affected by the noise-source distribution. In contrast, the C^{3} method is based on the use of the network stations as virtual sources (Stehly et al., 2008). In this case, the symmetry in the C^{3} function should no longer reflect the noise-source distribution, but the distribution of the network stations around the path A-B.

### 3.2 Time symmetry of C^{3} and C^{1}

To verify these theoretical expectations, we compare the time symmetry of both correlation functions (C^{1} and C^{3}) for the particular case of the station pair HAU-BOURR (Fig. 3). In the C^{3} correlation process, we chose to select only the 55 network stations south of HAU-BOURR (see map in Fig. 3a). This means that the theoretical C^{3} virtual sources are localized south of the station pair of interest and opposite to the incident-noise direction, as the main source of noise is located in the North Atlantic ocean or on the northern European coast (e.g., Friedrich et al., 1998; Kedar et al., 2008; Landès et al., 2010; Stehly et al., 2006). Fig. 3b shows the results of both the C^{1} and C^{3} correlation functions for this geometry. The time symmetry of both signals shows two opposite peaks of maximum amplitude, which reflect a main energy flux propagating in opposite directions: from north to south for C^{1}, and from south to north for C^{3}. Note that the same convention is used for all of the correlation functions, i.e. the signal in negative (resp. positive) correlation times corresponds to waves propagating from north to south (resp. south to north). This observation indicates that the time symmetry of the two correlation functions is consistent with the location of their expected contributing source distributions.

However, the C^{3} function has a better time symmetry than C^{1}. This indicates that the coda waves used to compute the C^{3} function constitute a more isotropic field than the ambient seismic noise. This can be explained by the nature of the coda waves, which correspond to waves scattered from the heterogeneities in the Earth subsurface (e.g., Aki and Chouet, 1975). This scattering process is expected to make the coda wavefield more isotropic, depending on the distribution of the scatterers in the medium and the lapse-time considered in the coda (Paul et al., 2005). To confirm the influence of the scattering in our observations, we compare the C^{3} results for the same station pair (HAU-BOURR) computed for two different time windows T, in the C^{1} correlation (Fig. 4). In the first case, we consider a 1200-s time window in C^{1}, including the direct surface waves (Fig. 4a). This time window begins 15 s before the Rayleigh wave travel-time (i.e., twice the dominant period of the signals). We denote it as T_{all}, as it takes into account both the ballistic and the scattered waves. Due to the strong amplitude of the direct waves in the time window, we expect that their contribution will dominate in the C^{3} process. The resulting C^{3}_{all} function is shown in Fig. 4b. In the second case, we use C^{3}_{coda} for the C^{3} function calculated (as in Section 2) from a 1200-s time window T_{coda} of the correlation located entirely after the direct arrivals (Fig. 4c).

In view of the nature of the dominant waves involved in both of these computations, we can infer the origin of the energy flux observed in C^{3}_{all} and C^{3}_{coda}. The energy flux is expected to be totally controlled by the distribution of the virtual sources in the first case, whereas it must also be influenced by scattering (and the distribution of the scatterers) in the second case. The station distribution considered in this case helps in the separation of the different contributions to the time symmetry in the correlation process.

On the one hand, Fig. 4b confirms the role of network stations as virtual sources in the C^{3} method, as C^{3}_{all} shows complete time asymmetry, which reflects a directive energy flux in the opposite direction to the seismic noise direction that is consistent with the distribution of the network stations. On the other hand, the time symmetry is preserved in C^{3}_{coda} (Fig. 4c), which indicates the link between the scattering events and the wavefield isotropy in the correlation process. Finally, we notice the lower coherence observed in C^{3}_{coda} compared to C^{3}_{all}. This likely reflects the higher level of fluctuations in the coda of C^{1} than in its ballistic part. Note that the coherence in C^{3}_{all} is about twice as high as that in C^{1} in this geometry particularly favorable to C^{3}_{all} (i.e., a dense distribution of virtual sources in the alignment of receivers).

This analysis shows that the time symmetry of the C^{3} function does not depend any more on the distribution of the noise sources, but is instead essentially controlled by the distribution of the stations in the network when direct waves are selected, and is also influenced by the contribution of the scatterers when the C^{3} process is performed from the coda. This means that the sources in the C^{3} method are controllable to a certain extent, which is a noticeable advantage compared to the C^{1} process, which usually suffers from uneven noise-source distribution.

## 4 Combining information from the different correlation functions

The C^{3} method is an alternative method for the reconstruction from seismic noise of the Green's function between two points, and we have shown that it can improve noise-based measurements, especially by making correlation functions independent of the noise source. However, this iterative method can also be viewed as a way to obtain several estimates of the Green's function between two points that can be used to complement each other for tomography purposes. Indeed, noise-based tomography requires unbiased travel-time measurements between station pairs. In the case of a directive noise, it has been shown that the noise-correlation functions carry the footprint of the noise-source direction, which can lead to biased travel-time estimations. In general, when both the causal and anticausal parts of the noise-based correlation functions are reconstructed, the travel-time estimation is said to be unbiased. Further bias might come from a poor signal-to-noise ratio in the correlation functions, which also prevents satisfactory measurements of the travel-time.

We show in Fig. 4 for a given station pair, that the combination of the C^{1} and C^{3} results may alleviate these bias. In this case, we have seen that both the C^{3}_{coda} and C^{3}_{all} functions provide information with respect to C^{1} since they exhibit also the time-reversed counterpart of the ballistic signal reconstructed in C^{1}. Note that in the general case, this information is provided by C^{3}_{coda}, since the time symmetry of the C^{3}_{all} function is strongly controlled by the station distribution. There are therefore configurations for which the information in C^{1} and C^{3}_{all} is redundant. On the other hand, in terms of coherence, C^{1} and C^{3}_{all} can provide an advantage with respect to C^{3}_{coda} for better-quality measurements to be obtained. Besides, Fig. 4 shows a higher level of coherence in C^{3}_{all} than in C^{1}.

Several aspects have roles in the improvement of the correlation results. In the present study, we have mentioned in particular the noise-source effects, time symmetry and coherence. The C^{3} functions complement C^{1} by contributing to one or several of these aspects. The different correlation functions can thus be viewed as complementary information that can be considered together to improve noise-based tomography analysis.

## 5 Iterating the correlation process

The results presented in the previous sections have shown that it is possible to reconstruct (at least partially) the Green's function by iterating the correlation process from C^{1} to C^{3}. These results emphasize that the late arrivals in the C^{1} function still carry information about the propagation medium since the correlation of the C^{1} coda averaged over the station networks allows the reconstruction of the direct arrivals.

We can then ask about the next iterative step, i.e. is it possible to use the coda of the C^{3} function to reconstruct part of the Green's function?

For Fig. 5, we computed the C^{5} function^{1} for a particular station pair from the coda of the C^{3} functions averaged over the network station pairs. Note that in Fig. 5 we compare correlation functions instead of Green's function. The relation between the Green's function and the correlation function (Eq. (2)) refers implicitly to a time derivative that is written as $\mathrm{Im}\left(G\right)\propto i\omega \times C$ in the frequency domain. When dealing with correlation of correlations, an extra ω^{2} − term is introduced into the counterpart of Eq. (2), which becomes $\mathrm{Im}\left(G\right)\propto i{\omega}^{3}\times {C}^{3}$. This implies a difference in the spectral content of the Green's function and the correlation function that becomes more pronounced during the iteration process. However, the spectral whitening used in our processing maintains the spectrum as almost constant in a narrow band during the iteration. This allows us to consider that the spectrum variations remain negligible throughout the iteration process.

The good agreement between the direct arrivals of the C^{5} and the C^{1} functions shows that some coherent signal is still present in the coda of the C^{3} function. This result supports the idea that the C^{5} function could also be used in noise-based tomography or for monitoring techniques. However, we observe in Fig. 5 that the signal-to-noise ratio in the C^{5} function is smaller than that in the C^{3} function, which was also smaller than that in the C^{1} function (note the decrease in the coherence level in Fig. 5b, c and d). This might limit the practical interest of this iterative correlation process in geophysical inversion techniques.

Sabra et al. (2005c) proposed a theoretical expression for the level of fluctuations in the C^{1} function that helps in the understanding of the Green's function reconstruction during the iterative correlation process. They considered a white-noise model and a finite bandwidth B. In the case of normalized correlations (Eq. (3)), the expected level of fluctuations, n, is given by:

$$n=\frac{1}{\sqrt{2BT}}$$ | (7) |

The stack over N virtual sources, as performed in the C^{3} method, leads to an expected fluctuation level in the C^{3} function of (e.g., Larose et al., 2008; Stehly et al., 2008):

$$n=\frac{1}{\sqrt{2NBT}}$$ | (8) |

Measurements of the level of fluctuations carried out for the iterative correlation functions have shown good agreement with Sabra et al. (2005c). This means that the fluctuation levels observed (Fig. 5) are indeed controlled by the evolution of B, T and N during the iterative correlation process. In the following, we examine the evolution of these different parameters in this iteration process:

1. The frequency bandwidth B: in the frequency domain, the correlation consists of multiplying the spectra. If the frequency content is peaked, the correlation process heightens the contrast between the frequencies, which, in practice, amounts to having a narrower spectral content (B is decreasing). This corresponds to a loss of information, which in Eq. (7) leads to higher fluctuation levels in the correlations. In our case, the noise spectrum between 0.1 Hz and 0.2 Hz is peaked near 0.14 Hz, and this peak will be accentuated throughout the successive iterative steps. To prevent this loss of spectral information, we whiten the selected coda between 0.1 Hz and 0.2 Hz, so as to correlate a flat spectrum signal in the same frequency band at each iteration step.

2. The duration T: it is well established that correlation results are highly dependent on the duration of the ambient noise time series. In the present study, C^{1} has been calculated from one year of seismic ambient noise, and C^{3} results from the average of about 200 correlations (C^{3++} and C^{3−−} averaged over a hundred stations) performed on a time window of T = 1200 s extracted from the C^{1} coda. Considering Eq. (8), this results in larger fluctuations in the C^{3} function than in C^{1}. As the fluctuation level is higher in C^{3} than in C^{1}, the C^{5} function was computed from the selection of a shorter time window in the C^{3} coda, which results in a stronger fluctuation level in the C^{5} function. This explains why the level of fluctuations will always increase through the correlation iterations after the C^{3} processing.

3. The number of virtual sources N: Eq. (8) gives a level of fluctuations that is decreasing as $1/\sqrt{N}$. This means that the stacking over many virtual sources improves the C^{3} processing. However, N remains the same after the C^{3} function, and therefore this parameter does not influence the fluctuation level in the higher-order correlation functions.

This analysis shows that the level of fluctuations can be improved in the C^{3} correlation (with respect to the C^{1} correlation) by considering large networks (and thus more virtual sources), although it is limited by the decrease in the time interval T after the C^{3} iteration.

## 6 Conclusions

The present study shows that the later part of the noise correlation contains coherent signals that can allow the reconstruction of the Green's function by the so-called C^{3} method. Indeed, the C^{3} correlation function clearly shows the direct surface-wave part of the Green's function, and the success of the next iterative step, which is C^{5}, shows that coherent signals are also present in the coda of the C^{3} function.

This method can be useful to improve noise-based measurements, especially by suppressing source effects that are caused by non-isotropic source distributions. Two points should be emphasized. First, the C^{3} method is computed from scattered waves that have already been extracted in the C^{1} function. Secondly, we have shown that the C^{3} correlation function loses its dependency on the noise-source spatial distribution, and instead highly depends on the network stations that have the role of virtual sources. Thus, both scattering and an even spatial distribution of stations help in converging toward a two-sided Green's function with this technique. The C^{3} method can therefore help in the resolution of problems that arise from uneven source distributions in the primary noise correlations (C^{1}).

However, fluctuations in C^{3}_{coda} are higher than in C^{1}, and this might mitigate the practical interest in this higher-order correlation function for tomography purposes. We can adopt a somewhat different approach by combining information from the different correlation functions. Indeed, the C^{1} and C^{3} correlation functions provide several travel-time estimations that can be viewed as complementary information for noise-based tomography inversions.

## Acknowledgements

All the seismic data used in this study have been obtained from the IRIS DMC (http://www.iris.edu/), the ORFEUS database (http://www.orfeus-eu.org/), the ETH Zürich, the CEA (Commissariat à l’Énergie Atomique, France), and the LGIT (Laboratoire de Géophysique Interne et Tectonophysique de Grenoble, France). We thank N. Shapiro, U. Wegler and D. Draganov for their helpful comments and suggestions. One of the authors (B.F.) acknowledges the support of Shell Research. This research was supported by European Research Council (Advanced grant Whisper).

^{1} As for C^{3}, the exponent indicates the number of C in the correlation function description. The C^{5} function corresponds to the Correlation of Coda of the C^{3} function.