On the strength of the phase cross-correlation in retrieving the Green’s function information in a region affected by persistent aftershock sequences

Although research on seismic interferometry is now entering a phase of maturity, earthquakes are still the most troublesome issues that plague the process in real applications. To address the problems that arise from spatially scattered and temporally transient enormous earthquakes, preference is usually given to the use of time-dependent weights. However, small earthquakes can also have a disturbing effect on the accuracy of interpretations if they are persistently clustered right next to the perpendicular bisector of the line joining station pairs or in close proximity to one of the stations. With regard to the suppression of these cluster earthquakes, commonly used solutions for dealing with monochromatic microseismic cluster events (e.g., implementing a band-reject filter around a comparatively narrow frequency band or whitening the amplitude spectra before calculating the cross-spectrum between two signals) may not have the necessary efficiency since earthquake clusters are generally a collection of events with different magnitudes, each having its own frequency and energy contents. Therefore, the only solution left in such a situation is to use stronger non-linear time-dependent weights (e.g., square of the running average or one-bit normalization), which may cause Green’s function amplitude information to be lost. In this paper, by simulating the records of a benchmark earthquake MN 5.2 with the help of empirical Green’s functions (EGF) obtained after the Ahar-Varzeghan Earthquake Doublet (MN 6.4 and MN 6.3), it is shown that the amplitude-unbiased phase cross-correlation is a relatively efficient approach in the face of the issues concerning long-standing cluster events.

Abstract Although research on seismic interferometry is now entering a phase of maturity, earthquakes are still the most troublesome issues that plague the process in real applications. To address the problems that arise from spatially scattered and temporally transient enormous earthquakes, preference is usually given to the use of time-dependent weights. However, small earthquakes can also have a disturbing effect on the accuracy of interpretations if they are persistently clustered right next to the perpendicular bisector of the line joining station pairs or in close proximity to one of the stations. With regard to the suppression of these cluster earthquakes, commonly used solutions for dealing with monochromatic microseismic cluster events (e.g., implementing a band-reject filter around a comparatively narrow frequency band or whitening the amplitude spectra before calculating the crossspectrum between two signals) may not have the necessary efficiency since earthquake clusters are generally a collection of events with different magnitudes, each having its own frequency and energy contents. Therefore, the only solution left in such a situation is to use stronger non-linear time-dependent weights (e.g., square of the running average or one-bit normalization), which may cause Green's function amplitude information to be lost. In this paper, by simulating the records of a benchmark earthquake M N 5.2 with the help of empirical Green's functions (EGF) obtained after the Ahar-Varzeghan Earthquake Doublet (M N 6.4 and M N 6.3), it is shown that the amplitude-unbiased phase cross-correlation is a relatively efficient approach in the face of the issues concerning long-standing cluster events.
Keywords Seismic interferometry . Empirical Green's function . Spatially localized and temporally persistent events . Phase cross-correlation

Introduction
Historically speaking, seismologists have long dreamed of turning ambient seismic noise into useful signals. https://doi.org/10.1007/s10950-021-10008-1 Article highlights • Spatially localized and temporally persistent small events have a destructive impact on correlation results in the seismic interferometry.
• Conventional pre-processing methods cannot completely eliminate these adverse impacts.
• Phase cross-correlation partially addresses concerns regarding these persistent cluster events. Since the beginning of the current century, a new field of research called seismic interferometry has provided an overarching theoretical framework to translate this visionary ambition into a concrete reality. In brief, theoretical studies have shown that the ensemble-averaging of cross-correlations of random seismic noises records over available time windows is a straightforward way to retrieve the Empirical Green's functions (EGFs) between a pair of stations (for a comprehensive review, see Wapenaar et al. 2010a, b;Campillo et al. 2014;Nakata et al. 2019 and references therein), provided that certain conditions are fulfilled. The first precondition for the successful implementation of this method is that noise sources are uncorrelated such that the crosscorrelations of signals coming from simultaneously acting sources located in different places (the "crossterms") nullify each other reciprocally (e.g., Snieder 2004;Weemstra et al. 2014;Medeiros et al. 2015). Another prerequisite is the availability of a diffusive seismic noise. In strictly theoretical terms, a perfect diffuse system is defined under the theorem of the modal energy equipartitioning (e.g., Hodgson 1996;Hennino et al. 2001;Shapiro et al. 2000;Weaver and Lobkis 2006;Sánchez-Sesma et al. 2008). Fundamentally, a wavefield can be assumed to be equipartition when, in a given frequency band, all the modes of the system are statistically equiprobable, and the share of energy between these modes is balanced in such a way that the net flux becomes zero (e.g., Margerin et al. 2009;Margerin 2017). In general, the equipartitioning required for performing the seismic interferometry is deemed to be fulfilled in a finite medium with an irregular bounding surface or at the stage of the isotropically developed multiple scattering processes in a heterogeneous medium (e.g., Margerin et al. 2009;Pilz and Parolai 2014).
In the absence of a fully scattered wavefield, it is imperative that studies be carried out at a long lapse of time to make sure station pairs are sufficiently surrounded by a spatially homogeneous distribution of random noise sources (shown in black circles in Fig. 1). The rationale behind this idea is that if the azimuthal homogeneity of the incoming noise is met, the coherent energy sent from the sources located in the Fresnel zone, i.e., the endfire lobes along with the line joining a pair station, contributes constructively in the crosscorrelation at a given lag-time. On the other hand, in such conditions, most of those energy produced by the sources outside the Fresnel zone cancel one another out by destructive interferences (see, e.g., Snieder and Hagerty 2004;Roux et al. 2005;Lin et al. 2008;Harmon et al. 2010).
The most important challenge that researchers face when attempting to do a study over a longer period of time is that it makes the correlation results more vulnerable to the disturbing influence of spatially scattered and temporally transient (SS-TT) high-energy impulsive sources (e.g., earthquakes, explosive) (shown in blue circles in Fig. 1). These events by imposing too much energy along certain azimuths prevent the energy equipartitioning condition from being met. Matters become even more complicated when clusters of continually repeating events are grouped in clusters along the isophase hyperbolas with foci at two stations in the off-Fresnel zone (shown in green circles in Fig. 1). In such cases, spatially localized and temporally persistent (SL-TP) small events may create a series of unwelcome bumps, sometimes called precursory noise (see Yanovskaya and Koroleva 2011;Ritzwoller and Feng 2018), in the earlier time-lags of cross-correlation results if they occur right next to the perpendicular bisector of the line joining two stations (shown in perpendicular orange dash line in Fig. 1). As the cluster events shift from the center toward one of the stations, these cluster events manifest their problematic effect in the form of a frequency-dependent asymmetry in the correlation results (Frank et al. 2009, see Fig. 8;Zheng et al. 2011). Contrary to the common view that taking the average of the causal and anticausal parts of correlation functions is a way to homogenizing incoming noise azimuths (e.g., Lin et al. 2008), Zheng et al. (2011, Fig. 5) showed that this averaging not only does not help to mitigate such frequencydependent asymmetry but also may cause a significant bias in the group and phase velocity measurements.
In this paper, we simulated a benchmark earthquake M N 5.2 using the empirical Green's functions (EGFs) obtained from the amplitude-unbiased phase crosscorrelation (PCC) analysis (Schimmel 1999;Schimmel et al. 2011b;D'Hour et al. 2016;Schimmel et al. 2018;Ventosa et al. 2019) in a three-month seismic active timeframe after the Ahar-Varzeghan Earthquake Doublet (M N 6.4 and M N 6.3) to examine the operational effectiveness and efficiency of this method in regards to mitigating the precursory noise created by cluster events when extracting coherent signals from the uncorrelated ambient noise data. A comparison of the consistency of the results provides evidence to support this claim that the PCC appears to be an optimal option in the face of the issues concerning the retrieval of Green's function information, even in the presence of long-standing cluster events.

Correlation techniques
The seismic interferometry offers an innovative basis to obtain an approximate assessment of the empirical Green's function (i.e., EGF) between station pairs, simply by the ensemble average (i.e., stack) of cross-correlations of received uncorrelated random seismic noise. These correlations can be obtained by calculating the inverse Fourier transform of the cross-spectrum of two signals, Þand u k m X B ; ω ð Þare respectively fixed-length signals (e.g., hourly data) recorded by the component of l at the station X A and the component of m at the station X B in the angular frequency domain, ω, also k = 1, …, K denotes the index of the hourly time windows and * stands for the complex conjugate. When it comes to the actual implementation, earthquakes are major barriers on the way to achieving the satisfactory results in this method. So far, various time-domain pre-processing approaches have been proposed to lessen the disturbing influence of SL-TP events (Ritzwoller and Feng 2018). In the case of monochromatic microseismic SL-TP events, it is also a common practice to use a band-reject filter around a comparatively narrow frequency bandwidth (Schimmel et al. 2011b, see Subsection 3.1.2) or to whiten amplitude spectra, before calculating the cross-spectrum between two signals (Shapiro et al. 2006, see Discussion Section;Ritzwoller and Feng 2018). Advantage of spectral whitening is that it broadens the bandwidth of the estimated empirical Green's functions (Seats et al. 2012). It should be borne in mind, however, that the aftershock clusters cannot be viewed as monochromatic signals, since they are generally composed of a collection of events with different magnitudes, each having its own frequency and energy contents. Therefore, the conventional pre-processing techniques may not be of much help as regards to aftershock clusters that may affect a broad frequency range of the correlation results.
To make the frequency domain whitening as adaptable as possible, attention may be drawn to use frequencydependent whitening approaches as e Φ k lm X A ð ; X B ; ωÞ≡e u k l X A ; ω ð Þe u k m X B ; ω ð Þ * , where e u ω ð Þ ¼ u ω ð Þ=N ω ð Þ is the spectrum of ground motion after normalizing with a frequency-dependent real-valued function, N(ω). The particular case that N k l X A ; ω ð Þ¼N k m X B ; ω ð Þ¼ u k l ω ð Þ È É is called the deconvolution method (e.g., Prieto et al. 2011;Viens et al. 2014Viens et al. , 2015Viens et al. , 2016Viens et al. , 2017 Fig. 1 A schematic drawing illustrating the spatially scattered and temporally transient enormous events (the blue solid circles), the spatially localized and temporally persistent small events (a cluster of red solid circles), a line joining two stations (the horizontal black dash line), the azimuthally homogeneous distribution of uncorrelated noise sources (the black solid circles), the perpendicular bisector of the line joining two stations (the perpendicular orange dash line), and endfire lobes. The series of temporally persistent events that are spatially localized along the isophase hyperbolas with foci at two stations in the off-Fresnel zone are shown by the green solid circles where the operator || indicates the absolute value and the operator {} represents the smoothing of the spectra using a moving average over predefined points, δ DEC is the regularization parameter that has to be added to prevent numerical instability. This regularization partly relaxes the dependence of the correlation analysis to the features of excitations (Prieto and Beroza 2008;Prieto et al. 2011).
, we end up in the crosscoherency method as (e.g., Prieto and Beroza 2008;Prieto et al. 2011 where δ CCOH is the regularization parameter for the cross-coherency method. By introducing the amplitude-unbiased PCC method, Schimmel (1999) contributed substantially to the new ways of thinking about calculating the coherence of two signals in terms of their instantaneous phase (Schimmel et al. 2011b). The PCC is obtained according to the following equation: where θ l (X A ) and θ m (X B ) represent the phases belonging to the stations A and B that their elements are the instantaneous phase of three components at the instantaneous time t. It is important to note that in this equation, the results have been normalized to θ k lm ≤ 1. Therefore, θ k lm ¼ 1 indicates that there is a perfect correlation between two signals and θ k lm ¼ −1 implies a complete anticorrelated feature between two signals.

Region of study and data
In an attempt to evaluate the robustness of the correlation methods in addressing the issues concerning the retrieval of the true Green's function information in the presence of long-standing cluster events, we limit our focus to the northwestern Iran where the region was in continuous vibration by sustained aftershock sequences, since 11 August 2012 when a runaway release of a stored elastic strain of the Ahar-Varzeghan Earthquake Doublet (M N 6.4 and 6.3, M N represents the Nuttli's magnitude scale, Nuttli 1973) triggered a succession of localized aftershock sequences. In total, 869 aftershocks with magnitudes greater than M N 2.5 have jolted this region at a period from the first of August 2012 to the end of August 2013. The yearly pattern of events occurring in a rectangular box bounded by latitude (min): 36.50°, latitude (max): 39°longitude (min): 45.50°and longitude (max): 49°and their magnitude distribution are shown in Fig. 2.
The data used in this study were recorded by threecomponent broad-band and the short-period sensors in different-length Nanometrics Y-File and the GCF fragments. Stations used in this study belong to Iranian Seismological Centre, Institute of Geophysics, University of Tehran (IRSC-IGUT, http://irsc.ut.ac.ir), and International Institute of Earthquake Engineering and Seismology (IIEES, http://www.iiees.ac.ir/en/). We applied a baseline correction by the removal of mean and trends from all fragments, before using them. Thereafter, following applying a cosine taper, fragments were corrected to the ground velocity (for a polarity analysis, see Section 4.1) and ground displacement (for an interferometry analysis, see Section 4.2) by removing instrumental responses. The different-length fragments were merged to obtain a 1-h-long times-series for all three-component records. All windows were visually inspected to determine whether there has been any indication of failure to record or poor quality data. If so, these windows were excluded from the processing cycle. Finally, data are bandpass-filtered between 0.01 and 2 Hz and then decimated up to a Nyquist frequency of 2 Hz.
4 Method: comparison of EGFs obtained in two seismic active and quiescent timeframes

Background noise analysis
We selected the ZNJK station in Northwestern Iran as our virtual source, and tried to make the pairwise comparison of the EGFs of the hourly three-component displacement noise data of this station and the other six stations. To investigate the effect of aftershock clusters on cross-correlation results, we performed this comparison for two different timeframes: a threemonth-length seismic active timeframe and a 3-monthlength quiescent timeframe. The seismic active timeframe was selected from 2 days after the occurrence of the doublet earthquake to 13 November 2012 (Fig.  3a). During this interval, stations have recorded 623 earthquakes with a magnitude larger than M N 2.5 of which 37 earthquakes had a magnitude larger than or equal to M N 4. Another 3-month-length timeframe was selected from a seismic quiescent time period: 01 May to 01 August 2013 in which only 26 events with a magnitude of M N 2.5-3.5 occurred (Fig. 3b).
In this paper, the timescales of post-seismic healing after the occurrence of each event were considered very short. Additionally, it is assumed that the accumulation of stress before triggering an event does not affect the nature and content of the noise. Therefore, with such assumptions, the background noise received by the stations in the seismic active timeframe is expected to have the same pattern as the seismic quiescent timeframe, provided that the seasonal differences are negligible. We evaluated the validity of this assumption by the information obtained from instantaneous back azimuths (BAZs) of the elliptical motions as a function of time and frequency using frequency-dependent polarization analysis (FDPA) in the region (Schimmel andGallart 2003, 2004;Schimmel et al. 2011a;Davy et al. 2015;Berbellini et al. 2019;Carvalho et al. 2019). In line with this strategy, an adaptive-length Gaussian-shaped sliding window moves along the three-component data vector. The length of the Gaussian window is inversely proportional to the frequency and used to construct a frequency-dependent cross-spectral matrix. Then, this matrix is decomposed into its eigenvalues and corresponding eigenvectors, helping us determine semimajor and semi-minor axes of an ellipse that best fits the ground motion. Using these results, in addition to the instantaneous BAZs, the variation of the instantaneous planarity vectors of polarized signals with respect to the mean direction of their vectors can be measured by a parameter called the degree-of-polarization (DOP). The planarity vector is the perpendicular vector to the instantaneous plane of the ground motion ellipse. It should be noted that the DOP changes in a range between 0 and 1. In the cases in which the semi-major axis direction is strongly varying, it is expected to have a smaller value for the DOP, while for a perfectly polarized motion the DOP will be 1.
We performed polarity analysis on both the seismic active and quiescent timeframes by setting the DOP in the range (0.7-1). This constraint allowed us to downweight less polarized signals and to increase the intensity of weakly coherent Rayleigh waves. We obtained the results by accumulating the outputs of polarity analysis for 180-s time segments in which no earthquake occurred. Bootstrap resample analysis was used before performing analysis on the time segments to reduce the effect of very small cluster earthquakes in space that were not visually recognizable, especially for seismic active intervals. Figure 4 shows the variations of the BAZ, as a function of frequency, at each seismic station. The number of signals has been normalized to 1 per day, and the color scale has been saturated to 0.5.
Considering the results (Fig. 4), it is particularly evident that coherent noise in the lower-frequency band (0.1-0.5 Hz) entered the station from all directions; nevertheless, a degree of directivity is also evident in the patterns of their azimuth-frequency distribution. At most stations, this directivity is observed along the northeast and southwest, which may be attributed to secondary microseismic waves. However, the SRB station shows a predominant direction of elliptically polarized seismic waves in the southeast direction, which may be considered the effect of regional structures, especially the directivity due to the large-scale topography response (Burjánek et al. 2014;Massa et al. 2014). Despite the complexities noted above, what matters is that the patterns are almost unchanged in two seismic active and quiescent timeframes. By this token, it can be assumed that any difference in the results of correlation methods is most likely due to the presence of cluster earthquakes.
4.2 Efficiency of applying a time-domain binary filter to eliminate the precursory noise When performing the correlation, we first selected hourly time intervals randomly from both the seismic active and quiescent timeframes so that the gradual decline of the rate of occurrence of aftershocks following the trigger of the mainshock does not introduce a systematic bias on the results. After removing the effect of the enormous spikes, glitches, and SS-TT events with a time-domain binary filter on the data (i.e., setting 0 for time periods having peaks larger than a certain threshold value and 1 for other time periods) (e.g., Ritzwoller and Feng 2018), the adverse effect on the correlation results is only due to the cluster of small earthquakes (Pedersen and Krüger 2007). To make it easier to extract the EGFs of Rayleigh and Love waves, we rotated the coordinate system from northeast-down coordinates to radial transversedown coordinates before calculating the correlation. To accelerate the convergence of Green's function, we overlapped the time windows by 30 min (Seats et al. 2012).
In Fig. 5, the EGFs obtained for the vertical-tovertical cross-correlations are shown in the microseismic frequency range (0.01-0.3 Hz). The causal window is an interval between times corresponding to group speeds of 2 and 5 km/s, and its counterpart is shown by an acausal window. A comparison of the SNR (i.e., the ratio of the maximum amplitude in the causal window to the RMS obtained for precursory noise) at the seismic quiescent (5a) and active timeframes (5b) for the ZNJN-MRD station pair decreased from 2.09 to 0.27. As such, the time-domain binary filter does not seem to be doing very much in mitigating the harm posed by the aftershock clusters, particularly because the amplitude of the signals generated by some of the earthquakes was much smaller than the selected threshold, so they were able to pass through the time-domain binary. In other words, a time-domain binary filter assigns 0 for time periods having peaks larger than a certain threshold value and 1 for other time periods. If the amplitudes of the aftershocks are smaller than this threshold, they are safe from this filter and become further magnified in the ensemble average (stacking) operation.
Since the events are mostly clustered close to the perpendicular bisector of the line joining the ZNJN-MRD station pair, the precursory noise is more pronounced in the initial lags of their crosscorrection function. As shown in Fig. 3, however, the location of the cluster events is inclined towards one of the pairs of stations in the cases of the ZNJN-TBZ and ZNJN-AZR station pairs. Because of that, the precursory noise in these two cases is seen in the later lags, and the SNR obtained for these station pairs shows a decrease in the seismic active timeframe compared to the seismic quiescent timeframe (with a drop from 2.17 to 0.35 for the ZNJN-TBZ pair and from 2.15 to 0.66 for the ZNJN-AZR station pair). In the cases of ZNJN-HRS, ZNJN-BST, and ZNJN-SRB, the clusters are mostly located outside the space between the station pairs. Therefore, the effect of cluster events shows themselves as a clear asymmetry in the correlation results rather than precursory noise. Examples of such asymmetries have already been reported by studies conducted over microseismic cluster events at the Gulf of Guinea (Shapiro et al. 2006), near the Kyushu Island (Zeng and Ni 2010), and in shallow waters of the Juan de Fuca plate Ritzwoller 2015, 2017). As such, aftershock clusters are expected to bear a close resemblance to the microseismic clusters, because these events continuously shake a limited area of space, months, or years after the main shock. 4.3 Efficiency of applying a frequency-dependent whitening approach to eliminate the precursory noise Now, we try to study the efficiency of the frequencydependent whitening approaches in mitigating the earthquake clusters. To this end, we choose the regularization parameter δ DEC to be 1% of the average spectral power of u l (X A , ω) when applying the deconvolution method (Nakata et al. 2011). In the cross-coherency method, when the spectral amplitude is small, the numerator and denominator are both small. Therefore, the regularization parameter δ CCOH can be selected much lower than δ DEC . According to Fig. 5c and Fig. 5d, the deconvolution and cross-coherency methods have both a positive effect on increasing the SNR of the correlation results for ZNJN-HRS, ZNJN-BST, and ZNJN-SRB station pairs. However, the use of these methods was Fig. 4 We performed the polarity analysis on both the seismic quiescent (first column) and active (second column) timeframes by setting the DOP in the range (0.7-1) for different stations in the frequency range of 0.1 to 0.5 Hz. The number of signals has been normalized to 1 per day, and the color scale has been saturated to 0.5. As can be seen, the coherent incoming waves enter the station from all directions; nevertheless, a degree of directivity is also evident in the patterns of their azimuth-frequency distribution.
Additionally, in most stations, this directivity is observed along the northeast and southwest, which may be attributed to the microseismic waves. However, the SRB station shows the predominant direction of elliptically polarized seismic waves in the southeast direction. Despite the complexities noted above, what matters is that the patterns are almost unchanged in two seismic active and quiescent timeframes  Fig. 5 The vertical-to-vertical correlation results between ZNJN and six stations located at northwest of Iran obtained by the crosscorrelation analysis at seismic (a) quiescent and (b) active timeframes. (c) Results obtained by the deconvolution method at the seismic active timeframe, (d) results obtained by the crosscoherency method at the seismic active timeframe, and (e) results obtained by the PCC method at the seismic active timeframe. The signal-to-noise ratio (SNR) is inserted at the top left of each plot. In calculating the SNRs, we have considered the signal window to be a causal window and the noise window to be the precursory window (the gray-colored box). A comparison of the results shows that the PCC method is better than the other methods to eliminate the impact of precursory noises at active time interval (for details see the main text) not able to eliminate the impact of the precursory bumps and the asymmetry of the ZNJN-TBZ, ZNJN-TBZ, and ZNJN-AZR results.

Evaluating the robustness of the PCC method in retrieving the true Green's function by simulating a benchmark earthquake
Looking at the results in Fig. 5e, it seems that the PCC method can better overcome the impact of cluster events than other correlation methods because the vertical-tovertical PCC results are less affected by precursory noise, and the asymmetry between the causal and acausal parts of the correlation results is insignificant. To more accurately assess the ability of this method, we try to simulate the records of a benchmark earthquake with the help of EGFs obtained from the PCC noise analysis. To this end, we first compute the pairwise PCC between the three components of the station pairs in random hourly time intervals in the active timeframe. If simulated signals are fitted with an acceptable degree of accuracy to the original records of the earthquake, it can be concluded that the PCC can retrieve a reliable estimate of Green's function, even in the presence of cluster aftershocks.

Setting a benchmark
To put the idea described above into practice, it is necessary to establish an appropriate benchmark earthquake. Generally speaking, a benchmark event is referred to as a candidate earthquake that its epicenter locates very close to the location of the virtual source. As a rule of thumb, a benchmark earthquake should not be an event smaller than M W 4 since smaller events usually do not have an adequate signal-to-noise ratio at longer periods. Furthermore, it should not be an event stronger than M W 5, since the source time function of strong events is so complex (i.e., a M w 4-5 earthquake is that the earthquake can be regarded as a point source at long periods > 3s). In line with these requirements, we selected a M N 5.2 event occurred on 2008 May 27 at 06:18:08 as our benchmark earthquake (shown with a red star in Fig. 3). Table 1 shows information about the epicenter of this event reported by IRSC. From the inversion of the first pulse of the phase onset together with the sign information of the amplitudes recorded by the Northwest Iranian stations, an oblique-reverse mechanism was obtained for this earthquake. The seismic moment tensor analysis package, Hybrid (Kwiatek et al. 2016), was used to obtain this mechanism. Table 2 provides a concise description of the moment tensor solutions that we obtained for these events.

Converting EGFs obtained from noise to earthquake Earth Green's response
Before proceeding further with the analysis, a set of corrections should be made upon the noise EGF tensors to make them quite similar to a real earthquake Green's response function. To this end, we first assume that there is no energy leakage to off-diagonal entries of the EGF tensors due to the Love and Rayleigh waves coupling. Based on this assumption, we calibrate the amplitude of the non-zero entries of the EGF tensors to that of the earthquake using a set of calibration factors by the method described in Denolle et al.'s study (2013).
Comparing the epicenter of this event with the location of the ZNJK station (36.67 N, 48.68 E) indicates that these two are only 4.93 km away from each other which compared to the average interstation distance (in the order of 254 km) is somewhat negligible. Due to the proximity of the epicenter of the benchmark earthquake and the location of the virtual source, we suppose that the attenuation and geometrical spreading effects for seismic waves propagating from benchmark-to-stations and virtual source-to-stations are almost identical. Furthermore, the arrival time differences between the simulated signals and the waves received from the benchmark earthquake were ignored. However, the difference between the excitation depth of the near-surface noise sources and the centroid depth of the benchmark earthquake remains a matter of serious concern, as spectrums of seismic waves strongly depend on their excitation depth (i.e., high-frequency surface waves are more likely to be excited for shallow events than deep events). To address this difference, we pursued the course of action outlined in Denolle et al.'s work (2013).
The first thing needed for the depth correction is to determine the radial displacement eigenfunctions associated with Rayleigh waves, r 1 (ω), the vertical displacement eigenfunctions associated with Rayleigh waves, r 2 (ω), and the transverse displacement eigenfunctions associated with Love waves, l 1 (ω), for the region surrounding ZNJN station both at the centroid depth of the benchmark earthquake, h, and at the surface of Earth (see Denolle et al. 2013Denolle et al. , 2014. Figure 6 shows the ratios of these displacement eigenfunctions that were obtained by the Chebyshev Spectral Collocation algorithm (Denolle et al. 2012, VEA Package). The velocity model and density profiles in calculating these parameters were adapted from the model obtained by Donner et al. (2013) (Table 3). Although this model was achieved by analyzing the group-velocity dispersion curves of surface waves obtained at the Alborz mountains region, however, Donner et al. (2014Donner et al. ( , 2015 demonstrated that it is also applicable to a broad regional scale including the NW of Iran. In the next step, the effect of the radiation pattern of the benchmark earthquake was taken into account by inserting the moment tensor solutions in the surface-wave displacement equations (see Donner et al. 2013, Eqs. 13 to 15). This effect on the results can be seen by comparing the absolute values of these converting terms in the polar coordinate system at given periods 5, 7, and 10 s (Fig. 7).
An earthquake source is generally considered a spontaneous rupture propagation along the fault. Therefore, convolving the modified noise displacement fields with the displacement source time function of the benchmark earthquake makes it more compatible with the benchmark earthquake records (Denolle et al. 2013, Eq. 17). Frequently, the source time function of a moderate event is considered a pulse with the width of T. This pulse width is inversely proportional to the corner frequency f c as T = 1/2f c (Udias 1999). This corner frequency is a threshold frequency beyond which the falloff of the source spectrum is proportional to f 2 and can be obtained by f c = 0.491 β(Δσ/M 0 ) where Δσ is called the stress drop (Hanks and Thatcher 1972). This stress drop is typically considered Δσ = 3 MPa and shear velocity β is supposed to be β =3 km/s. We obtained f c = 0.36 Hz for the corner frequency of the benchmark event. Finally, we applied a conservative free phase shift up to 2 s upon the correlation results on account of the possible uncertainties regarding errors in the determination of the location of benchmark events, their excitation depths, their fault plane solutions, and so on.

4.4.3
The degree of consistency between the waveforms and amplitudes of simulated and observed records Figure 8 shows waveforms obtained for the records of simulated signals and the records received from the benchmark earthquake. A visual comparison of these results indicates that there is a certain degree of consistency between the amplitude and frequency content of these two records, but any claim in this regard must be made on a mathematical basis. As a metric to evaluate the accuracy of the PCC method in eliminating the disturbing impact of the persistent sequence of the cluster aftershocks, therefore, we used the following parameters to compare the degree of consistency between the waveform and amplitude information of simulated and observed records: (1) A normalized correlation coefficient function (hereinafter called CCF) that forms a standard to evaluate the consistency between the travel time information of these two sets of data Viens et al. 2014): CCF ≈ ∑ N 92:5 i¼N 2:5 x i y i ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ∑ N 92:5 where x is the Green's function and y is the earthquake record in the time domain. In this equation, N 92.5 and N 2.5 stand for the 2.5 and 92.5% of the cumulative energy of both signals, respectively.
(2) Peak ground displacement (hereinafter called PGD) differences (in percentage) that quantitatively assesses the accuracy of the amplitude retrieved as: where PGD PCC and PGD eq are respectively PGDs of the simulated and earthquake wavefield, and PGD max = max(PGD noise , PGD eq ). The CCFs for the vertical component of the simulated results in the ZNJK-AZR and ZNJK-TBZ pairs exhibit a rather low correlation (≤ 0.6) in comparison with other pairs. In our view, the low CCF values obtained for the PCC results of the ZNJK-AZR and ZNJK-TBZ pairs may be relevant to the complexities created by the low-velocity large-scale heterogeneities that have been reported by previous studies at the eastern lower crust margins of the Sahand inactive volcanic region (see Rezaeifar et al. 2016;Bavali et al. 2016). Because of these large-scale anomalies, it seems that factors affecting the wave propagation of seismic noise on the surface of the Earth, which are slightly different from factors affecting the wave propagation of the benchmark earthquake and the complexities of the propagation path, were not compensated merely by applying the simple corrections provided above. However, what is clear from the comparisons is that the PCC method has been able to simulate amplitudes and waveforms of the benchmark earthquake at the SRB, BST, and MRD stations in a relatively acceptable way. By comparing the results obtained for Fig. 6 Ratio of the radial (Rayleigh) r 1 (h)/r 1 (0), transverse (Love) l 1 (h)/l 1 (0), and vertical (Rayleigh) r 2 (h)/r 2 (0) displacement eigenfunctions at the source depth and the surface for the region surrounding ZNJN station. These ratios were obtained based on the velocity and density profiles under seismic ZNJN station (Table 2) PGD max values, it can be said that although the matching accuracy of the amplitudes is less than the phase information, the use of the PCC method has resulted in better retrieval of amplitude information in some components of the station pairs.

Discussion and conclusion
When retrieving EGFs from noise data in the microseismic frequency range, the main focus is often on eliminating the damaging effects of very large earthquakes, while less attention is paid to small earthquakes. However, Jousset and Douglas (2007) showed that small earthquakes may also make a major contribution to the amplitude spectrum at longer periods. Moreover, we observed that the disturbing effect of temporarily persistent small earthquakes may be much greater, especially in two specific cases: (1) If these earthquakes are clustered around the perpendicular bisector of the line joining station pairs, they cause precursory noise to be observed in the earlier time-lags of cross-correlation results; (2) If these earthquakes are clustered in the close proximity to one of the stations, they may lead to the frequency-dependent asymmetry in the correlation results. Experiences to date have demonstrated that stronger non-linear operations such as the square of the running average (Yanovskaya and Koroleva 2011) or one-bit normalization (Ritzwoller and Feng 2018) may be useful in mitigating the adverse effects of clusters of small earthquakes. However, a number of previous scholarly papers have provided evidence demonstrating that the efficiency of these methods is mainly limited to cases where we wanted to extract information on seismic velocities (e.g., Sabra et al. 2005;Cupillard and Capdeville 2010;Cupillard et al. 2011), while some uncertainty, even controversy, surrounds the application of these methods to address the problem of properly extracting Green's function amplitudes (Weaver 2011;Prieto et al. 2011;Lawrence et al. 2013;Weemstra et al. 2014). In general, the retrieval of Green's function amplitude information from noise correlation analysis is extremely important in applications such as ground motion simulation, attenuation estimation, evaluation of site amplification effects, and anisotropy analysis (e.g., Prieto et al. 2011). Therefore, there is widespread consensus in the literature that the use of pre-processing schemes should be limited to keep Green's function amplitude information as intact as possible.
Another important point to note is that all timedomain and frequency-domain normalization methods destroy the uniqueness of waveforms, which has two major drawbacks. First, they increase the time required to extract the signal in the stacking process. Second, they may worsen the situation regarding the destructive effect of the signal associated with spatially localized and temporally persistent small events due to the increased ambiguity between signals and noise (D'hour et al. 2016;Schimmel et al. 2018).
According to the explanation provided in this paper, it seems that the PCC can be a way to address the issues we discussed above, perhaps because in this method, waveform similarity is measured through the number of Fig. 7 The absolute values of converting terms (see Section 6.2 in the main text) are shown in polar plots at given periods T= 5 (red dashed line), 7 (orange dashed line), and 10 s (black dashed line). We obtained the factors at those periods by a flat response assumption for the moment-rate function. The moment tensors used to obtain these values are listed in Table 1. The maximum amplitude at the azimuth 15°also shown at these plots phase-coherent samples rather than the sum of amplitude products (Schimmel et al. 2011a, b). As the epicenter of the cluster earthquakes inclined towards the perpendicular bisector of the line joining two stations, the waves generated by them travel in different paths before reaching the stations. Therefore, the recorded Fig. 8 Prediction of the 2008 MN 5.2 benchmark earthquake by the modified and calibrated noise impulse responses. Waveforms were filtered between 0.01 and 0.3 Hz. The radial, transverse and vertical components of the simulated records are shown by red blue color, while records of the benchmark earthquake are shown by blue. The normalized correlation coefficient function, CCF, and the percentage of retrieved amplitude errors, PGD diff , are shown at the upper right and left part of each subplot, respectively signals from these station pairs have different coherence under the influence of different anomalies and heterogeneities along their path. This can prevent the formation of precursory noise and frequency-dependent asymmetry in the correlation results. Additionally, the amplitude unbiasedness of the PCC allows it to be used under difficult circumstances with a high-amplitude variability of signals and noise, without the need for implementing pre-processing methods in mitigating the influence of energetic features such as strong earthquakes, sensor failures, and local episodic noise (Schimmel et al. 2018). Therefore, the PCC may partially address concerns regarding the retrieval of Green's function amplitude information. Of course, just doing a study on a small number of stations cannot guarantee the PCC's efficiency in this regard. Therefore, this conclusion must be examined by other studies. In addition, the fact that the PCC method does not require pre-processing techniques can also be an advantage in the face of accelerating signal extraction. The relatively satisfactory results obtained for EGFs with the help of 3-month data in this study can be evidence of this claim.
It should be noted that the performance of the PCC may also be affected by the temporal variability in the seismic anisotropy of background seismic noise due to the alignment of the fractures and cracks that accompanies major earthquakes (e.g., Saade et al. 2017). This may reduce the efficiency of the PCC method, since in such cases, signals may be less phase-coherent than in the surrounding noise (Schimmel et al. 2018). This may have been the cause of the slight mismatch observed at some stations.