# Ambient noise tomography with non-uniform noise sources and low aperture networks: case study of deep geothermal reservoirs in northern Alsace, France

Maximilien Lehujeur, Jérôme Vergne, Alessia Maggi, J. Schmittbuhl

## ► To cite this version:

Maximilien Lehujeur, Jérôme Vergne, Alessia Maggi, J. Schmittbuhl. Ambient noise tomography with non-uniform noise sources and low aperture networks: case study of deep geothermal reservoirs in northern Alsace, France. Geophysical Journal International, 2017, 208 (1), 10.1093/gji/ggw373. hal-01482004

HAL Id: hal-01482004

https://hal.science/hal-01482004

Submitted on 19 Oct 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ambient noise tomography with non-uniform noise sources and low aperture networks: case study of deep geothermal reservoirs in northern Alsace, France

Maximilien Lehujeur, Jérôme Vergne, Alessia Maggi and Jean Schmittbuhl

IPGS, Université de Strasbourg/EOST, CNRS, 5 rue Renéé Descartes, F-67084 Strasbourg Cedex, France. E-mail: lehujeur@unistra.fr

Accepted 2016 October 5. Received 2016 September 29; in original form 2016 April 25

## SUMMARY

We developed and applied a method for ambient noise surface wave tomography that can deal with noise cross-correlation functions governed to first order by a non-uniform distribution of the ambient seismic noise sources. The method inverts the azimuthal distribution of noise sources that are assumed to be far from the network, together with the spatial variations of the phase and group velocities on an optimized irregular grid. Direct modelling of the two-sided noise correlation functions avoids dispersion curve picking on every station pair and minimizes analyst intervention. The method involves station pairs spaced by distances down to a fraction of a wavelength, thereby bringing additional information for tomography. After validating the method on synthetic data, we applied it to a set of long-term continuous waveforms acquired around the geothermal sites at Soultz-sous-Forêts and Rittershoffen (Northern Alsace, France). For networks with limited aperture, we show that taking the azimuthal variations of the noise energy into account has significant impact on the surface wave dispersion maps. We obtained regional phase and group velocity models in the 1–7 s period range, which is sensitive to the structures encompassing the geothermal reservoirs. The ambient noise in our dataset originates from two main directions, the northern Atlantic Ocean and the Mediterranean Sea, and is dominated by the first Rayleigh wave overtone in the 2–5 s period range.

**Key words:** Surface waves and free oscillations; Seismic tomography.

## 1 INTRODUCTION

Cross-correlation of ambient seismic noise records is now commonly used to estimate the Green’s function of a medium between pairs of seismic receivers (Shapiro & Campillo 2004; Snieder 2004; Roux et al. 2005a). It has been used in various contexts and at different scales (e.g. Sabra et al. 2005; Shapiro et al. 2005; Yang et al. 2007; Lin et al. 2009). Although several studies have successfully isolated and interpreted body waves in this manner (e.g. Roux et al. 2005b; Poli et al. 2012; Lin & Tsai 2013; Nakata et al. 2015), surface waves remain the primary wave-types extracted from the cross-correlation of ambient seismic noise (Campillo et al. 2011).

Recovering the true Green’s function between a pair of receivers requires either a diffuse wave field (e.g. Lobkis & Weaver 2001; Campillo & Paul 2003) or a uniform distribution of ballistic noise sources (Snieder 2004; Roux et al. 2005a). A non-uniform distribution of noise sources in a non-diffusive medium can induce errors in the station-to-station travel-time measurements. These errors remain negligible if the interstation distance is longer than a few wavelengths and if the spatial distribution of the noise sources is smooth (Tsai 2009; Weaver et al. 2009; Yao & Van der Hilst 2009; Froment et al. 2010), however some studies have reported significant

biases on the surface wave velocity measurements caused by directive seismic noise, especially in the microseismic band (Shapiro et al. 2006; Pedersen & Krüger 2007). Several methods have been proposed to correct such biases with no prior knowledge of the spatial distribution of the noise sources. The ‘directional balancing’ technique (Curtis & Halliday 2010) constructs a virtual source with near-uniform directionality by decomposing the noise field in space at one of the receivers using array processing. The ‘multi-dimensional deconvolution’ technique (e.g. Wapenaar et al. 2008, 2011) computes a point-spread function from the correlation of the noise field over large and regular arrays. This function contains the signature of the noise source spatial distribution, which can be efficiently removed from the observed correlation functions. The double-beam-forming technique (Boué et al. 2014) uses array-processing methods to isolate the wave fields from sources located within the Fresnel zones that contribute constructively to the Green’s function.

Dispersion measurements, especially group-speeds, are assumed to be less accurate for station pairs closer than 2–3 wavelengths (Bensen et al. 2007). Most studies based on the correlation of ambient noise for tomography therefore exclude such station pairs (e.g. Yao et al. 2006; Yang et al. 2007; Lin et al. 2008, 2009;
Stehly et al. 2009; Young et al. 2011), thereby reducing the lateral resolution of tomographic images at the periods and wavelengths required to probe deep structures. Recent case-studies have suggested that the distance threshold for dispersion measurements may be reduced down to 1 wavelength (e.g. Luo et al. 2015; Zigone et al. 2015), which extends application of the correlation method to small-aperture networks. The reliability of travel time measurements obtained from the shortest station pairs would still need to be verified in each case.

We propose a way to determine the Rayleigh wave velocity model and the spatial distribution of the noise sources simultaneously, by interpreting the noise correlation function (NCF) waveforms, including pairs with interstation distances below 2 or 3 wavelengths. We developed the method in order to perform local-scale ambient noise surface wave tomography while avoiding the traveltime errors caused by an unknown non-uniform distribution of noise sources. We applied it to the deep geothermal field near Soultz-sous-forêts and Rittershoffen (Alsace, France). Ambient seismic noise is particularly suited to this context because the region lacks natural seismicity and its geothermal operators require imaging methods that are cheaper than active seismic surveys.

To image the first 5 km of the crust using surface waves, we need to focus on their dispersive behaviour in the 1–7 s period range. This requires working with ambient noise dominated by the secondary microseismic peak that originates from oceanic sources (e.g. Gutenberg 1936; Longuet-Higgins 1950; Kedar et al. 2008; Sergeant et al. 2013) located in coastal or pelagic zones (e.g. Beucler et al. 2015). Local seismic stations are installed around geothermal sites primarily for the purpose of monitoring induced seismicity. Their limited aperture is sub-optimal for imaging purposes and requires us to exploit station pairs closer than 2 wavelengths.

The paper focuses on the methodology proposed to handle the case of directive noise recorded by a limited aperture network. The specificities of this case are explained using actual data recorded near the two geothermal sites. We describe our procedure for determining group and phase dispersion maps together with the azimuthal variation of the noise energy from the two-sided NCF waveforms. We explain the forward problem, compare it to other studies, test the inversion procedure on synthetic data and apply it to the observed NCFs.

## 2 DATA AND PROCESSING

The data used in this study result from a combination of several permanent and temporary networks deployed near the geothermal sites of Rittershoffen and Soultz-sous-forêts (Fig. 1, black stars) between August 2009 and end of 2014. The two short-period permanent networks (SZ and RT, Fig. 1) were the first to be installed in the area and are managed by the Ecole et Observatoire des Sciences de la Terre (EOST) and ES-Gothérmie (ES-G). In addition, a network of 16 short-period temporary stations provided by the Karlsruhe Institute of Technology (KIT) and the GeoForshungZentrum (GFZ) was deployed within a 5 km radius around the Rittershoffen site between May 2013 and December 2014 (KIT1 in Fig. 1).

The RT, SZ and KIT1 networks were originally designed to monitor micro-seismic activity associated with the installation, stimulation and exploitation of the two geothermal sites. Their long-duration continuous recordings motivated us to perform an ambient noise tomography study of the deep geothermal reservoir, which is located at the interface between the basement and the sedimentary cover at about 2.5 to 3 km depth (Baujard 2015). However,

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Map of the area. Soultz and Ritt. indicate the two deep geothermal sites of Soultz-sous-Forêts and Rittershoffen and squares represent either permanent or temporary seismic stations from various networks (see insert and text for details). The thin white lines correspond to roads, the white dashed line is the French–German border and the black line is the Rhine, which runs along a portion of the border.</figcaption></figure>

the narrow aperture of these networks hampered attempts to image the medium at reservoir depth with surface waves extracted from NCFs. We extended this aperture in 2014 by deploying 7 temporary broad-band velocimeters on a 30 km diameter circle centred on Rittershoffen (Fig. 1, yellow squares).

All the possible pairs of sensors for the whole network yielded about 600 vertical-vertical NCFs computed over time periods extending from 1 to 5 yr. The pre-processing steps applied to the hourly sequenced noise records and intended to homogenize the seismic wavefield include instrument response correction, a 0.1 to 10 s spectral whitening and a 1-bit digitization. The resulting NCFs form the dataset of this study and are presented in Fig. 2.

## 3 NOISE DIRECTIVITY AND INFLUENCE ON THE NCFs

Despite the pre-processing steps, some characteristics of the ambient noise still show up in the NCFs. We observe asymmetry between causal and acausal amplitudes (Fig. 2a), known to be an effect of noise directivity (Stehly et al. 2006). The larger amplitudes observed on the causal NCFs indicate that most of the noise energy originates from the west, according to the convention chosen to orient the station pairs.

We also observe an azimuthal variation of the NCF phases that can be highlighted by grouping the NCFs having similar interstation distances and sorting each group as a function of station pair azimuth (Fig. 2b). This two-level sorting of the NCFs reveals shorter phase arrival times for azimuths around 0°N and 180°N. The azimuthal variation of the NCF phases in the Rittershoffen and Soultz-sous-forêts data was first described by Lehujeur et al. (2015), who ascribed it to strong directivity of the ambient noise in the period band dominated by the secondary microseismic peak.

To further highlight the effects of ambient noise directivity on the NCFs, we performed a frequency–wavenumber analysis (FK) over a time period during which all the sensors were available
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Noise correlation functions (NCFs) filtered between 1 and 5 s. Each line of the plot corresponds to an NCF from one single station pair, and the lines are sorted in two different ways. (a) NCFs are sorted by increasing interstation distance; by convention, the positive part of the NCF corresponds to noise sources located west to the station pairs. (b) NCFs are grouped by similar distance ranges (delimited by the black horizontal solid lines); each group is sorted as a function of azimuth of the station-pair (angle from 0° to 180° measured at the western station, clockwise from north).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Dominant noise propagation directions and resulting effects on the near-field NCFs. (a) Frequency–wavenumber (FK) transform obtained at 2.5 s period from the continuous vertical noise records acquired for 27 d in May 2014. During this period all the sensors of the networks were available simultaneously. Amplitudes are expressed in decibels. The white (resp. red) dashed circle indicates the locations in the ($k_x, k_y$) domain where the fundamental mode (resp. 1st overtone) of the Rayleigh waves is expected according to the local 1-D velocity model. A positive $k_x$ (resp. $k_y$) value corresponds to a plane wave propagating toward the east (resp. toward the north). The black dashed circle indicates a horizontal phase speed of 7 km s$^{-1}$. The white arrows indicate two dominant patterns in the $k_x - k_y$ domain corresponding to the noise coming from back-azimuths ~310°N and ~150°N. (b) Common Shotpoint Gather: NCFs filtered at 2.5 s of period and displayed in the space domain relative to a common virtual origin ($x = 0, y = 0$; white star) at time $t = 5.25$ s. The black and white dashed lines indicate the two dominant wavefronts that can be identified in the NCFs. The large and small black arrows indicate the back-azimuths ~310°N and ~150°N respectively.</figcaption></figure>

simultaneously. Fig. 3(a) shows the FK transform of the noise records at 2.5 s period. The different modes of the Rayleigh waves are expected to show up on circles whose radii depend on the phase-speed. The FK analysis has a marked peak for positive $k_x$ (i.e. waves propagating eastward, Fig. 3a, large white arrow), corresponding to nearly plane waves with a 310°N back-azimuth. A second, less en-

ergetic peak appears for positive $k_y$ values (i.e. waves propagating northward, Fig. 3a, small white arrow) corresponding to a 150°N back-azimuth. These two dominant peaks indicate surface waves with a phase speed of about 2.7 km s$^{-1}$, which is in good agreement with the theoretical phase speed of the first Rayleigh wave overtone (Fig. 3a, red dashed circle) computed using the programs
from Herrmann & Ammon (2002) and a local 1-D velocity model constrained by borehole measurements (Beauce et al. 1991; Cuenot et al. 2008). Application of the FK transform at different epochs and period ranges indicated that this noise directivity pattern was very stable over time between 1 and 7 s period (see supplementary materials S1).

As the NCF should tend towards the Green's function between each receiver pair, we can theoretically reconstruct the wave field emitted by each station and recorded by all others (e.g. Lin et al. 2009). In Fig. 3(b), we display a snapshot of all these spatial domain wave fields centred on a virtual common origin located at Cartesian coordinates $x = 0$ and $y = 0$ and filtered at 2.5 s. We refer to this kind of representation as a Common Shotpoint Gather (CSG). In the hypothetical case of a perfectly uniform noise distribution and laterally invariant velocity, a CSG snapshot would look like a circular wavefield centred on the virtual source. The part of Fig. 3(b) at $x > 0$ corresponds to the causal part of the NCFs (i.e. the part of the correlation functions affected by eastward propagating noise). We interpret this figure as the interference between two dominant wave fronts (Fig. 3b, black and white dashed lines) propagating away from two common sources as time increases (black arrows). Movie S2 in the supplementary materials shows the time history of the Fig. 3(b) CSG.

This observation confirms that our NCFs are controlled to first order by the main noise propagation directions. The significant phase distortions we observe are probably caused by the small extent of the network compared to the wavelength, as corroborated by Tsai (2009), Yao & Van der Hilst (2009) and Weaver et al. (2009) who find that the effect of a non-uniform noise distribution increases with decreasing interstation distance and increasing period. For small aperture networks in which most station pairs are shorter than a few wavelengths, these phase distortions become prominent and may induce significant errors in the resulting tomographic maps. The errors depend on the proportion of short and long station-pairs and the orientation of the station pairs with respect to the dominant source directions. Using the inner networks (RT, SZ and KIT1, Fig. 1), Lehujeur et al. (2015) showed that the strong phase distortions can cause errors on tomographic maps of up to 30 per cent at 2 s period. In the following we address the feasibility of correcting this first order effect to extract the correct spatial variation of Rayleigh wave velocity.

Rejecting station pairs with an interstation distance shorter than 2 or 3 wavelengths would reduce the effect of the phase distortions on the tomographic maps, but would severely decrease the amount of available data. In our case, rejecting station pairs closer than two wavelengths would remove 80 per cent of the dataset at 3 s of period (wavelength ~7.3 km) and 96 per cent at 4 s (wavelength ~10.4 km).

Instead of rejecting data, we propose to interpret all the NCFs by taking the noise directivity into account. To do so, we need to determine the spatial distribution of the noise and release the hypothesis that the NCFs are equal to the Green's function of the medium. This procedure is similar to how Yao & Van der Hilst (2009) inverted for phase speed anomalies, azimuthal noise intensity and azimuthal anisotropy. Their study was restricted to station pairs longer than the two-wavelength limit and they worked in a period band between 10 and 30 s, where the ambient noise has different characteristics than in the 1–7 s range used here. They concluded that the noise distribution had only minor effects on their phase dispersion maps. In our case, however, the azimuthal distribution of the noise energy dominates the NCFs (Figs 2 and 3) because of the small aperture of the network.

# 4 INVERSION PROCEDURE

In this study, we invert the two-sided noise correlation waveforms to determine simultaneously the spatial noise distribution and the Rayleigh wave group and phase dispersion maps at different periods. This procedure ensures safe velocity measurements by separating the effects of a heterogeneous medium from those caused by a non-uniform source distribution. It also avoids dispersion curve computation and picking on each station pair, which in turn reduces analyst involvement in the data processing and mitigates the risk of misinterpretation.

## 4.1 Forward problem

### 4.1.1 Statement of the forward problem

This section describes how we model the correlation waveforms for each station pair of a network, for a known velocity model and spatial distribution of noise sources.

To solve this problem we assume that all noise sources are located at the surface and far from the network and that the vertical component of the seismic noise is dominated by unimodal Rayleigh waves in the far-field domain. Under these assumptions - the validity of which is discussed in Section 5—a simplified expression of the Green's function can be written as:

$$G(r_A, r_S, v) = \frac{e^{-\gamma(v) SA} e^{-i(2\pi v \frac{S_A}{c(v)} + \varphi_0)}}{\sqrt{2\pi h(v) SA}} \quad (1)$$

where the subscripts S and A denote the source and receiver located at coordinates $r_S$ and $r_A$ respectively and separated by a distance $SA = \|r_S - r_A\|$, $v$ is the frequency, $c(v)$ is the apparent phase dispersion velocity along the SA path, $\gamma(v)$ is a positive real function controlling the intrinsic attenuation, $h(v)$ represents the thickness of the layer affected by the surface wave (the denominator corresponds to the geometrical attenuation of a cylindrical wave), and $\varphi_0$ is the phase at the source.

From this expression for the Green's function, we derive a simplified expression of the correlation function (see Appendix A for intermediate steps):

$$C_{AB}(v) = \frac{1}{2\pi h(v)} \int_{0}^{\pi} \mathcal{H}(\xi, v) e^{4i\pi v a \frac{\cos \xi}{c(v)}} d\xi \quad (2)$$

where $a$ is half the distance between the receivers A (eastern station) and B (western station) and $\xi$ is the hyperbolic coordinate (ranging from 0 to $\pi$ and measured counterclockwise from the A side of the station pair) controlling the curvature of a hyperbola whose focus points are the two stations A and B ($\xi = \pi/2$ is the perpendicular bisector of the AB segment, $\xi = 0$ (resp. $\xi = \pi$) corresponds to the hyperbola whose branches are merged with the path joining the two stations on the A (resp. B) side of the station pair, see Fig. A1 in Appendix A). The exponential term of eq. (2) is referred to as the interferometric phase. $\mathcal{H}(\xi, v)$ is called the hyperbolic contribution and represents the integrated amount of energy contributed to the correlation function from sources located along the hyperbola of parameter $\xi$. In the case of a uniform noise distribution, Snieder (2004) used the stationary phase integration technique to approximate eq. (2) by its behaviour near the angular domains where the interferometric phase is stationary (Fresnel zones). In this study, we need to account for the azimuthal distribution of the noise energy that is carried by the $\xi$-dependency of function $\mathcal{H}$, so we chose to solve eq. (2) by numerical integration over $\xi$ using a classical trapezium approach.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Synthetic cross-correlation functions obtained by numerical integration of the weighted interferometric phase. Left part: Imposed uniform (black) and non-uniform (red) distributions of the noise energy as a function of the hyperbolic parameter $\xi$. Central part: Interferometric phase weighted by the non-uniform noise distribution, displayed in the time domain and derived temporally (i.e. along the horizontal axis). The green bars correspond to the angular width of constructive interference area at different frequencies (Fresnel zones to within a quarter of period). Lower part: Synthetic correlation functions obtained for the two noise energy distributions by integrating the interferometric phase over the $\xi$ axis. Insert: Phase dispersion curve used for this example. The interstation distance is $2a = 15$ km and the period is bounded to the 1–5 s interval. In this specific example, the hyperbolic contribution is assumed to be independent of the frequency.</figcaption></figure>

As we assume that the noise sources are located far from the network, the spatial distribution of the noise energy can be reduced to a single azimuthal function (a similar assumption is proposed by Yao & Van der Hilst (2009) and referred to as a plane wave approximation). Under this assumption, we can write the hyperbolic contribution $H_\alpha$ for a pair oriented with an angle $\alpha$ measured at the western station clockwise from north as (see Appendix A for further details)

$$H_{\alpha}(\xi, \nu) = A(\alpha - \xi, \nu) + A(\alpha + \xi, \nu) \quad (3)$$

where function $A(\theta, \nu)$ is the azimuthal noise distribution in a north reference system and $\theta$ is the backazimuth measured clockwise from north.

Eq. (2) allows us to compute synthetic correlation waveforms (i.e. the $C_{AB}$ functions for every station pair of the network) by summing the contribution of the noise sources located on hyperbolas whose focal points are the receiver pairs. It represents the forward problem of our study and provides a non-linear relation between the model parameters (i.e. the azimuthal distribution of noise energy $A(\theta, \nu)$ and the phase dispersion curve $c(\nu)$) and the data domain (i.e. the correlation waveforms along each path). Fig. 4 shows an application of eq. (2) in the time domain for synthetic parameters. We generate a two-sided correlation function using a given dispersion curve (Fig. 4, insert) and either a uniform (left side, black curve) or a non-uniform (left side, red curve) azimuthal distribution of noise energy. This example illustrates some of the effects highlighted on the observed NCFs: asymmetry of the causal-acausal amplitudes,

distortion of the spectral content and frequency-dependent shift of both phase and envelope of the correlation waveform.

### 4.1.2 Validation of the forward problem

To test the reliability of the correlation waveforms modelled by the forward problem described above, we generated a set of synthetic correlation functions from an imposed azimuthal noise distribution (Fig. 5, blue curve). We computed a time shift value by comparing each correlation waveform to a reference correlation function obtained with the same amount of energy coming from all directions. This time shift represents the error on the Rayleigh wave phase arrival time when neglecting the effect of a non-uniform noise distribution. It is measured by computing the difference in the instantaneous phase at the group arrival time. The comparison was done at 2 s period by pre-filtering the two correlation functions using a Gaussian filter applied in the frequency domain. We then compared our measurements (Fig. 5, red and black dots) to the time-shifts predicted using the formula from Weaver et al. (2009) and Froment et al. (2010) (Fig. 5, solid red and black curves). We first measured the time shift for all possible orientations of a hypothetic station pair, keeping the interstation distance constant (2.6 times the wavelength, Fig. 5a). The observed and predicted time-shifts are in good agreement for all the orientations of the station pair and remain low (<4 per cent) since the azimuthal noise distribution imposed is smooth and the interstation distance is relatively long (Weaver et al. 2009). We also analysed the influence of the interstation spacing by
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison between noise correlation time error as predicted using the formula by Weaver *et al.* (2009) (solid curve labelled as ‘predictions’) and the time-error measured using our synthetic noise correlation functions (dots labelled as ‘measurements’) for a given synthetic azimuthal noise distribution (blue curve, expressed in arbitrary units). In this application, the period and phase velocity are 2 s and 1.94 km s⁻¹ respectively. (a) Effect of the station pair orientation for a fixed interstation spacing 2a = 10 km = 2.6λ. Red dots (resp. black dots) correspond to the time shifts measured on the acausal (resp. causal) part of the correlation function for a hypothetical station pair whose orientation varies from 0° (south-north) to 180° (north-south). The time shift (δt) is measured at the wave arrival time (t). The black dashed line indicate the pair orientation used for subplot b. (b) Effect of the interstation distance relative to the wavelength (2a/λ) for a fixed orientation of the station pair (260°N). Measurements (black dots) and predictions (black solid curve) are shown for the causal side of the correlation function. The black dashed line indicate the distance used for subplot a.</figcaption></figure>

fixing the pair orientation to 260°N and testing various interstation distances from a fraction of a wavelength to 5 wavelengths (Fig. 5b). The test confirms that the phase-error due to a non-uniform source distribution tends to zero when increasing the interstation distance. Below one wavelength, our observations are not consistent with theoretical predictions. This is due to the fact that the formula by Weaver *et al.* (2009) used to predict the time shifts is not valid for distances too short relative to the wavelength (their formula have a singularity for distances going to 0). In addition, our ‘measurements’ are difficult to interpret for short distances since they are obtained by picking the phase difference near the zero correlation lag-time. In the rest of this study, we interpret the whole correlation waveforms and not only the instantaneous phase at the wave arrival time. This procedure is more robust for short interstation distances. We always apply similar filters to synthetic and observed waveforms before computing the waveform misfit function.

### 4.1.3. Parametrization

The full correlation waveform allows us to determine both group and phase velocity of the Rayleigh waves. Because both azimuthal noise distribution and Rayleigh wave velocities are frequency dependent, we chose to solve an inverse problem for each of several periods in the range 1–7 s. Each inversion is bounded to the neighbourhood of a period 1/v₀ using a Fourier domain gaussian taper (e.g. Bensen *et al.* 2007, eq. 6 with parameter α = 15) applied to both the synthetic and the observed NCF waveforms. Our model parameters are given by a vector *m* whose components correspond to the Rayleigh wave group and phase velocity in the *N* cells of a grid and to the amplitude of the noise energy in *M* directions.

(1) The first *N* components of the model vector (noted $c_0, \dots, c_{N-1}$) contain the Rayleigh wave phase velocity at period 1/v₀ in each cell of the grid.

(2) The next *N* components of the model vector (noted $l_0, \dots, l_{N-1}$) contain to the first derivative of the dispersion curve at period 1/v₀ in each cell. In the neighbourhood of frequency v₀, we assume that the phase dispersion curve ($c(v)$) can be linearly approximated in the decimal logarithmic domain ($\log(v), \log(c)$):

$$c(v) \approx 10^{l(v_0)(\log(v)-\log(v_0)) + \log(c(v_0))} \quad \text{where } l(v_0) = \frac{d \log c}{d \log v} \Big|_{v_0} \\ = v_0 \frac{c'(v_0)}{c(v_0)} \qquad (4)$$

The group velocity at period 1/v₀ can be directly inferred from the pair of parameters ($c(v_0), l(v_0)$) using:

$$u(v_0) = 2\pi \frac{dv}{dk} \Big|_{v_0} = \frac{c(v_0)}{1 - v_0 \frac{c'(v_0)}{c(v_0)}} = \frac{c(v_0)}{1 - l(v_0)} \qquad (5)$$

where $k(v) = 2\pi v/c(v)$ is the wave number. The group velocity in the *j*th cell of the grid is thus $u_j = c_j/(1-l_j)$.

(3) The last *M* components of the model vector (noted $A_0, \dots, A_{M-1}$) contain a discretized version of function $A(\theta, v_0)$ for theta ranging from 0 to 2π and represent the amount of energy received by the network from each of the *M* directions.

We end up with model space of dimension 2N + *M*. We assume that all the parameters are independent and therefore that the covariance matrix on the model parameters is diagonal. The effective velocity parameters encountered along the path joining the *i*th
<figure><img src="image_7.png" /><figcaption>Figure 6. Parametrization of the dispersion maps. (a) Spatial coverage of the network obtained from the cumulated length of interstation rays crossing each pixel, in km⁻¹. (b) Azimuthal coverage of the network defined as the complementary angle to the largest azimuthal gap found in each pixel, in degrees. (c) Sum of both the normalized spatial and the azimuthal coverage maps as an indicator of the ‘network quality’. This indicator is smoothed laterally and used to place the Voronoï seeds (white dots) that control the mesh (solid lines).</figcaption></figure>

station pair of the network are obtained by summing the contributions of the structures along the interstation ray-path

$$ egin{cases} c_i(v_0) = \frac{\sum_j d_{ij}}{\sum_j d_{ij} / c_j(v_0)} \\ l_i(v_0) = \frac{\sum_j d_{ij} l_j(v_0) / c_j(v_0)}{\sum_j d_{ij} c_j(v_0)} \end{cases} \quad (6) $$

where $d_{ij}$ is the fraction of the $i$th interstation ray that crosses the $j$th cell, expressed in distance units.

To minimize the number of cells of the grid, and thus the number of parameters in the model, we use an irregular Voronoï mesh constructed from the spatial and azimuthal coverage of the network. We define spatial coverage as the cumulated length of interstation rays that cross a given pixel of the map per surface unit (expressed in km⁻¹, Fig. 6a). We define azimuthal coverage as the complementary angle of the largest azimuthal gap found in each pixel (ranging from 0 to 180 degrees, Fig. 6b). These two functions are normalized to their maxima and averaged to give us a ‘network quality map’ (Fig. 6c, colour scale) that is smoothed and used to place the Voronoï seeds (Fig. 6c, white dots) that control the size and density of Voronoï grid cells (Fig. 6c, solid lines). The outer cells are bounded by an 18 km radius circle centred on the Rittershoffen site. The resulting grid has 81 cells ($N = 81$) with surfaces ranging from 0.7 to 84.4 km². Most of the small cells are clustered in the central part of the map, where the network quality indicator is high (>0.25).

## 4.2 Inverse problem, application to synthetic data

### 4.2.1 Inverse problem

The inverse problem corresponding to the forward problem given by 2 is nonlinear. We define the a posteriori probability density function (PDF) as

$$ \sigma_M(m) = \exp \left( -\frac{1}{2} \left( \sum_{i=0}^{I-1} \left( \frac{g(m)_i - d_i^{\text{obs}}}{\sigma_i} \right)^2 + \sum_{j=0}^{N-1} \frac{S_j}{S} \left( \frac{c_j - c_j^{\text{prior}}}{\sigma_c} \right)^2 + \sum_{j=0}^{N-1} \frac{S_j}{S} \left( \frac{l_j - l_j^{\text{prior}}}{\sigma_l} \right)^2 + \sum_{k=0}^{M-1} \left( \frac{A_k - A_k^{\text{prior}}}{\sigma_A} \right)^2 \right) \right) \quad (7) $$

With $d^{\text{obs}}$ a length I data vector made of all the samples of the NCFs filtered around frequency $v_0$. $\sigma_l$ is the uncertainty associated

to the $i$th component of $d^{\text{obs}}$ taken as the standard deviation of the NCF in the coda part, $m$ is the model vector described in the previous section, $g(m)$ is a vector of synthetic data obtained from eq. (2) and filtered in the same way as the real data, $c_j^{\text{prior}}$, $l_j^{\text{prior}}$ and $A_k^{\text{prior}}$ are the components of the a priori model and $\sigma_c$, $\sigma_l$ and $\sigma_A$ are the corresponding uncertainties (their adjustment allows us to control the regularization of the inverse problem), $S_j$ denotes the area of the $j$th cell and $S$ is the total area of the grid. The quadratic misfits on the $c_j$ and $l_j$ parameters are weighted by the $S_j/S$ coefficients because the parameter space is not a normed space due to the irregular mesh (e.g. Spakman & Bijwaard 2001). We search for the maximum of $\sigma_M$ using an iterative gradient method: at each iteration we compute the gradient of $\sigma_M$ and we perform a grid search along the steepest ascent direction. This method is fast and can handle high-dimensional problems, however, the solution might depend on the model chosen to initiate the iteration (referred to as initial model in the following).

As our NCFs are mostly controlled by the spatial distribution of the noise sources (Figs 2 and 3), we decompose the inversion into two steps: in step 1, we reduce the grid to a single cell (i.e. an homogeneous model; $N = 1$), which means that the velocity model has only two parameters corresponding to the average phase and group speed, while the noise model has 36 parameters that are free to vary ($M = 36$). The solution of this step provides a preliminary noise model as well as the average phase and group velocity of the region. These parameters are used as initial model and prior information for step 2 of the inversion, in which the grid is refined into the $N = 81$ Voronoï grid cells (Fig. 6c). During step 2, both the noise parameters and the velocity parameters are adjusted to maximize the a posteriori PDF. This second inversion step estimates the spatial variations of the Rayleigh wave velocity around the average homogeneous model and refines the preliminary noise model.

### 4.2.2 Synthetic test

We tested the two-step inversion procedure at period 4.5 s on a set of synthetic NCFs generated from a synthetic noise distribution model (Fig. 7, blue curve) and a random group and phase velocity model (Figs 8a and b). These NCFs were obtained from the forward problem presented in Section 4.1 using the same network configuration as for real data (Fig. 1). A white Gaussian noise with a standard
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Synthetic noise model and results of the inversion tests A and B at 4.5 s of period. The corresponding dispersion maps are presented in Fig. 8. The corresponding data are presented in Fig. 9.</figcaption></figure>

deviation of 1 per cent of the maximum amplitude was added to them to simulate measurement inaccuracies (Fig. 9a).

The objectives of the synthetic test were (1) to estimate the error made on the velocity maps if the noise model is constrained to remain uniform, (2) to evaluate if our network can distinguish the effects caused by the non-uniform noise distribution from those controlled by the perturbations of the velocity model and (3) to test the influence of the initial model on the stability of the solution.

We performed two inversion tests (referred to as ‘A’ and ‘B’) on the same synthetic data set. In inversion ‘A’, the noise model was constrained to remain uniform in azimuth by using a single noise

parameter (i.e. *M* = 1). This procedure forced the inversion to converge toward a solution in which, the same amount of energy was coming from all directions, therefore simulating the assumption of an azimuthally uniform noise distribution. In inversion ‘B’, the noise model was subdivided into 36 parameters (*M* = 36) to account for azimuthal variations of the noise energy. This inversion test corresponded to the method developed in this study. The parametrization of the velocity model was identical for both synthetic tests A and B (*N* = 81). The noise and velocity models obtained after inversion tests A and B are presented in Figs 7–9 for comparison with the synthetic models used to generate the synthetic data set. For both inversions A and B, we tested several initial models with no significant impact on the solution. Inversion test B explained the expected noise model (Fig. 7, green histogram). In inversion test A, the amplitude of the noise converged to the average amount of energy of the expected model (Fig. 7, red histogram) since the noise model was restricted to a single-parameter.

The inner part of the group and phase dispersion maps obtained for both inversion tests A and B are presented in Fig. 8. The velocity values obtained in each cell were attributed to the centre of the Voronoï cells and smoothed laterally using a linear interpolation for graphical display. It appears that constraining the noise model to remain uniform in azimuth (inversion A, Figs 8c and d) induced severe errors in the velocity model (velocity anomalies with opposite signs or over-estimated amplitudes). Allowing azimuthal variations of the noise energy allowed us to recover the overall characteristics of the expected velocity model (inversion B, Figs 8e and f). The solution for phase speed dispersion was more accurate (Fig. 8e) than for group speed dispersion (Fig. 8f). The good quality dispersion maps from inversion B confirm that the correlation waveforms remained sensitive to the spatial variations of the Rayleigh wave

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Inner part of the synthetic dispersion maps and solutions of inversion tests A and B at 4.5 s of period. (a, b) Synthetic noise phase and group speed dispersion maps. The solid lines correspond to the Voronoï grid. (c, d) Results of inversion A in which the noise model is kept azimuthally uniform. The maps are smoothed by linear interpolation of the velocity values in each cell of the grid. The contour lines reproduce the expected velocity anomalies from maps (a) and (b) for comparison. (e, f) Results of inversion B. The red dots correspond to the stations.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Synthetic and inverted noise correlation functions at 4.5 s of period. The NCFs are sorted as a function of the interstation distance and azimuth. (a) Synthetic dataset obtained from the noise model of Fig. 7 (blue curve) and the velocity model of Figs 8(a) and (b). A white Gaussian noise is added to the waveforms to mimic measurement inaccuracies. (b, c) Data sets obtained from the solutions of inversion tests A and B respectively.</figcaption></figure>

velocity despite the highly non-uniform noise model imposed and the small scale of the network.

If we forward model NCFs from the solutions of inversion A and B (Figs 9b and c), we notice that inversion B explains the synthetic data (Fig. 9a) much more accurately. In particular, the causal-acausal asymmetry of amplitude, the phase distortions and the azimuthal variations of amplitude are well explained by inversion B, while these details cannot be explained by inversion A because it assumed a azimuthally uniform source distribution.

The synthetic tests above show that joint inversion of the noise and velocity parameters using our two-step procedure can safely estimate the group and phase dispersion maps in non-uniform noise conditions, despite the reduced aperture of the network compared to the Rayleigh wave wavelength. They also show that we need to take the non-uniformity of the noise model into account to estimate the dispersion maps correctly given our network configuration. This is particularly true for the longest periods (above 3 s). At shorter periods (synthetic tests at 2 s period are presented in the supplementary materials), the error on the dispersion maps diminished, however inverting both the speed anomalies and the azimuthal distribution of the noise energy improved the accuracy of the dispersion maps over inverting speed anomalies alone.

The synthetic test showed that there may be a trade-off between the noise model and the velocity model, meaning that an incorrect noise model may induce incorrect velocity estimations. This effect can be more or less pronounced depending on the network configuration and the targeted period. For the configuration of this study, it is strong enough to justify the need for joint inversion of the noise and velocity parameters. The solutions of the synthetic tests were only weakly dependent on the initial model, which means that the azimuthal and spatial coverage of the network are sufficient to permit unambiguous determination of phase and group speed models and azimuthal distribution of noise sources. The outer regions of the map are less well resolved because fewer station pairs were available and the cells are larger.

## 4.3 Application to observed data

We applied our procedure to the observed NCFs at 9 periods chosen between 1 and 7 s. In the first step of the inversion we estimated the average phase and group velocity as well as a preliminary noise model at each period (Fig. 10).

The noise models obtained for all the periods indicate a dominant direction of about 310°N, i.e., the Northern Atlantic Ocean (Fig. 10, polar histograms). This dominant direction varies slightly with period from about 315°N at 1 s to 280°N at 3.7 s. It is consistent with the results of the FK transform performed on the noise records (Fig. 3a and Supporting Information). Near 1 s period (Fig. 10, case i), the relative amount of energy received by the network from the northwest diminishes and the azimuthal distribution flattens out. This might be caused by the emergence of other sources such as the long period part of the local to regional anthropogenic noise (e.g. McNamara & Buland 2004; Groos & Ritter 2009) or seismic sources activated by wind (e.g. Withers et al. 1996; Bonnefoy-Claudet et al. 2006).

For periods below 2 s and above 5 s, the dispersion parameters converged toward the fundamental Rayleigh wave mode predicted by the local 1-D velocity model (see letters a, b, g, h, i and curve SZ0 in Fig. 10, left side). On the contrary, dispersion parameters obtained between 2 and 4 s were closer to the first overtone (letters c, d, e and f and SZ1 curve). A sharp velocity offset was observed near 4 s period (i.e. between b and c, Fig. 10). These results suggest that the first overtone dominated the NCFs between 2 and 4 s, as confirmed by performing a slant stack on the time domain NCFs (Fig. 10, left side, colour scale). It is likely that the first overtone dominated the NCFs in this period range because it also dominated the ambient seismic noise itself as revealed by the FK transform (see Fig. 3a and supplementary materials). We interpreted the velocity offset near 4 s period as a change in the relative weight of the two modes occurring near the cut-off period of the first overtone predicted by the 1-D local velocity model (upper bound of the SZ1
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Results of the first inversion step applied to the real NCFs at nine periods from 1 s to 7 s (denoted by letters a to i). During this step, the velocity model is assumed to be laterally homogeneous and the mesh has only one cell. Left side: Average phase speed (blue dots) and first derivative of the dispersion curve (blue lines) obtained for each period. SZ0, SZ1 and SZ2 denote the theoretical dispersion curves of the fundamental and the two first overtones of the Rayleigh waves computed from the local 1-D velocity model. The background colour scale represents the average phase speed measured via slant-stack on the NCFs at different periods. A value of 1 means good phase coherence of the NCFs. Right side: Azimuthal distribution of the noise energy obtained at each period (T) in arbitrary units, the same scale is used for all the polar histograms.</figcaption></figure>

curve). Based on these observations, we interpreted the velocity maps obtained between 2 and 4 s period (letters c, d, e and f) as corresponding to the first overtone of the Rayleigh wave, and those at other periods (letters a, b, g, h and i) as corresponding to its fundamental mode.

The synthetic wavefield predicted from this first step agrees with the observed NCF phases and relative amplitudes. The comparison between the observed and synthetic wavefields at 2.5 s period is presented in Fig. 11 using CSGs (common Shotpoint gathers) at different times. The synthetic correlation functions fit the data down to a fraction of a wavelength (see black dotted circles spaced by one wavelength). This confirms that the forward-problem theory we used in this study to model the NCFs remains valid for short paths. Additionally, we show that the overall characteristics of the NCFs can be explained during the first step of the inversion despite assuming a laterally homogeneous velocity model, which confirms that our NCFs were indeed controlled to first order by the noise directivity. Based on the results of the synthetic tests, we infer that the information concerning the spatial variations of the Rayleigh wave velocity (which will be extracted in the second step of the inversion) is contained in the second order discrepancies between the observed and the modelled datasets.

During the second step of the inversion, we released the assumption of a laterally homogeneous velocity and split the grid into 81 Voronoï cells. We set the *a priori* phase and group-velocity models to those obtained from the first step of the inversion. We regularized the inverse problem by adjusting $\sigma_c$ and $\sigma_l$ in eq. 7 (the uncertainties on the prior dispersion parameters). Optimal values of these parameters were determined using a grid-search approach (see Appendix B).

The phase and group dispersion maps obtained after the two inversion steps for periods 6.7, 2.5 and 1.1 s (letters a, e and i on Fig. 10) are presented in Fig. 12. The maps were smoothed laterally using the same parameters as for the synthetic inversion test (Fig. 8). The relative variations of the phase and group velocities are similar. A velocity gradient from the northwest to the southeast indicates decreasing velocities towards the central part of the Rhine Graben. This trend is particularly pronounced for dispersion maps dominated by the first overtone (e.g. maps e in Fig. 12). For the shortest periods dominated by the fundamental mode (e.g. maps i in Fig. 12), we observe short scale variations of the velocity, including a +5 per cent speed anomaly located under Soultz-sous-forêts and a -10 per cent located west of the Rittershoffen fault. Although no constraint is applied to the noise model parameters during the second inversion step (i.e. $\sigma_A$ is infinity in eq. 7), the final noise models (not shown) are similar to those from the first inversion step to within a constant factor.

## 5 DISCUSSION

We have developed a method that allows us to perform noise-correlation-based imaging in situations in which there is a strong sensitivity to the azimuthal distribution of noise sources, and applied it successfully to a small-aperture network. The proposed approach was motivated by the need to exploit NCFs from station pairs closer than two wavelengths. These data, although significantly influenced by the noise directivity, contained useful information for the construction of the dispersion maps. Including short station pairs in the data processing increases the lateral resolution. Furthermore, the methodology presented allows extending the upper limit of the observable period range and thus the sensitivity to the deep
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison between observed NCFs at 2.5 s of period (left columns) and inverted NCFs after the first step of the inversion (right columns). The NCFs are presented as Common Shotpoint Gather at different times. Thin dotted circles indicate the distance with respect to the source in wavelength units. The thick black dashed circle corresponds to a distance of 3 times the wavelength of the Rayleigh wave at 2.5 s.</figcaption></figure>

structures. In our case, we were able to interpret NCFs at periods up to 7 s, while a two-wavelength distance cut-off would have restricted that limit to about 3 s. We discuss in the following the specific assumptions made.

We assumed that the noise sources were located along a horizontal plane and sufficiently far from the network so that their spatial distribution could be reduced to a 1-D azimuthal function. This assumption arose because the secondary microseismic peak dominates the ambient seismic noise at periods between 1 and 10 s (e.g. Stutzmann *et al.* 2000). In northern Alsace, the main direction of the noise (~310°N) observed in both the noise records and the NCFs corresponds to the northern Atlantic, which has been identified as the origin of most of the secondary microseismic noise recorded

in Europe and north America (Gutenberg 1936; Kedar *et al.* 2008; Sergeant *et al.* 2013). The secondary incoming direction of the noise (~150°N) might be related to sources in the Mediterranean Sea or in the Indian Ocean. At periods below 1 s, the assumption of a azimuthal distribution of far-field noise sources would no longer have been valid since the noise at these periods is dominated by more local sources such as anthropogenic activity, especially in such densely populated zone (Lehujeur *et al.* 2015).

We also assumed that the Green's function was unimodal, that is, dominated either by the fundamental mode or by a single overtone depending on the frequency. This assumption seems reasonable because the ambient seismic noise and the correlation functions were clearly dominated by a single mode for all the 9 periods studied as
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Phase and group dispersion maps (upper and lower line respectively) obtained after inversion step 2 for periods 6.7 s (column labelled ‘a’ as in Fig. 10), 2.5 s (column ‘e’) and 1.1 s (column ‘i’). Maps a and i are attributed to the fundamental mode of the Rayleigh waves; maps e are attributed to the first overtone. Amplitudes are expressed in per cent relative to a reference speed ($c_{ref}$, and $u_{ref}$). The black dots correspond to the stations, the white stars marked S and R indicate the site of Soultz-sous-Forêts and Rittershoffen respectively. The horizontal white bars indicate the wavelength of the Rayleigh waves at each period. Boundaries of the Voronoï cells are shown as dotted lines. The solid black lines indicate the faults at top basement (from Baillieux *et al.* 2011). The purple line delineates the area (non-shaded), where the network quality exceeds 0.25 (Fig. 6c).</figcaption></figure>

proven by the FK and slant-stack analyses (Fig. 3a, supplementary materials and Fig. 10 left). The average phase and group velocities observed were consistent with theoretical dispersion curves computed from the local 1-D velocity model. More specifically, the velocity values were closer to the fundamental mode below 2 s and above 5 s and to the first overtone between 2 and 4 s. The sharp offset in the phase-dispersion diagram (Fig. 10, left) occurred near the expected cut-off period of the first overtone according to the local 1-D velocity model, further corroborating our modal interpretation.

We used the ray approximation to locate the velocity anomalies, that is, we attributed the observed velocity variations to the structures encountered along the interstation paths (eq. 6). This assumption is theoretically valid when the scale of the heterogeneities is large relative to the wavelength studied and applies ideally to infinite frequency signals. Here, we worked at finite frequency but the scale of the velocity anomalies observed remains large relative to the wavelength (Fig. 12). The velocity anomalies were spread over broader areas than the ray-paths because of the size of the grid-cells and the regularization of the inverse problem. These factors smoothed the lateral sensitivity, thus approximating the use of a more complex finite frequency sensitivity kernel projected onto the discretized grid.

The influence of the noise distribution decreases with increasing spacing as originally demonstrated by Weaver *et al.* (2009). As the effects of noise directivity are strongest for low interstation distances, the short aperture of our network compared to the wave-

length explains why the azimuthal biases dominate the NCFs. It would seem that the shortest station-pairs of the network mainly constrain the noise model, and that the long station-pairs are more sensitive to the velocity model. The case of station pairs with intermediate distances (between 1 and 3 wavelengths) seems to be more complex since the NCFs are sensitive to both the noise and the velocity models. For that pairs, the relative sensitivity to the noise/velocity parameters also depends on the complexity of the noise distribution and the orientation of the station pairs. For instance, in the extreme case of a plane wave emitted by a single source at infinite distance, the station pairs that are orthogonal to the propagation direction are weakly sensitive to the velocity parameters since the wavefronts reach the two stations at the same time whatever the velocity between them. The pairs that are aligned with the noise source (i.e. having the noise source within their Fresnel constructive zone) are mostly sensitive to the velocity parameters and only weakly sensitive to the location of the source. For the proportion of long station pairs to be sufficient to constrain the velocity model, the network aperture should not be too small. The stability and accuracy of the solution seem to rely strongly on the spatial and azimuthal coverage of the network, which can be examined using a synthetic approach. We have shown that our network was indeed able to determine both the noise and the velocity models unambiguously in the targeted period range. The outer regions of the network were covered by large Voronoï cells because of low spatial and azimuthal coverage, and consequently returned poorer resolution and lower relative velocity anomalies.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparison between the group dispersion maps obtained at period 3 s with the two-step inversion procedure that takes the noise distribution into account (a) and with a more standard approach (b). Group-velocity variations are given in per cent relative to a reference velocity of 1.81 km s⁻¹. The black ellipse indicate a positive velocity anomaly interpreted as a probable azimuthal effect caused by a local dominance of station pairs orthogonal to the main noise direction (~310°N).</figcaption></figure>

To further show the usefulness of taking the noise directivity into account, we compared the obtained group dispersion map to a more standard approach at period 3 s (after rejecting the pairs shorter than 2 wavelengths, Fig. 13, see supplementary materials for details). The differences between the two approaches are significant. The lower amplitudes observed in the ‘standard’ case (Fig. 13b) are interpreted as a consequence of the uniform-noise assumption in which, most of the travel time measurements cannot be explained by the inversion because of the azimuthal biases. Such biases forced us to use relatively high regularization constraints. The SW-NE positive anomaly (Fig. 13b, black ellipse) is understood as the result of a local dominance of station pairs oriented in that direction (i.e. perpendicular to the Atlantic noise). Except from that pattern, the sign of the velocity anomalies observed with both approaches is relatively consistent.

At periods below 2 s (Fig. 12i), the short scale velocity anomalies we recovered probably resulted from heterogeneities inside the sedimentary cover (e.g. local horsts and grabens) since the fundamental mode Rayleigh wave at these periods is mostly sensitive to structures within the first kilometre from the surface. The +5 per cent velocity anomaly observed under the region of Soultz-sous-Forêts might have been caused by a local Horst structure (Renard & Courrioux 1994; Genter et al. 2010; Baillieux et al. 2011; Dezayes et al. 2011). The sharp velocity contrast observed at Rittershoffen (Fig. 12i, label R) coincides with a westward dipping normal fault (Baillieux 2012; Baujard 2015) with lower velocities to the west.

The significant velocity gradient observed at periods above 2 s (Figs 12a and e; Fig. 14a) seems to be correlated to first order with the thickness of the sedimentary cover that increases from the Vosges massif to the northwest, towards the central axis of the Upper Rhine Graben to the east (Fig. 14b; Baillieux 2012). At these periods Rayleigh waves are sensitive to a thick layer extending down to the granitic basement. At Soultz-sous-Forêts (Fig. 14b, label S), the granitic basement depth is about 1.5 km (e.g. Cautru 1988; Genter & Traineau 1996), while it is about 2.5 km at Rittershoffen (Fig. 14b, label R; Baujard 2015). The two sites are located on horst structures bounded to the west by westward dipping faults

that probably contribute to positive thermal anomalies at the top of the basement (Baillieux et al. 2014, Fig. 14c). The change in the sign of the speed anomaly occurs near Rittershoffen (Fig. 14a), and coincides with an increase of the thickness of the low density Cenozoic sediments as shown from gravity survey (Rotstein et al. 2006; Fig. 14c, mark 3, late Oligocene layer).

# 6 CONCLUSIONS

We analysed the ability to perform ambient noise surface wave tomography at a local scale, in the case of directive seismic noise combined with a limited network aperture. Using data acquired near the geothermal sites of Soultz-sous-Forêts and Rittershoffen, France, we confirmed that such conditions induce significant distortion of the phase of the NCFs and prevent operators from applying the usual procedures for ambient noise surface wave tomography.

We overcame this issue by developing a two-step inversion of the NCF waveforms that could mitigate the effects of the non-uniform source distribution and produce phase and group dispersion maps of the area as well as a model of the azimuthal distribution of the noise energy. We confirmed the importance of taking the spatial distribution of the noise sources into account for ambient noise tomography on local networks using synthetic tests.

The ambient seismic noise recorded in our target region in the period range 1 s–7 s originated from two dominant directions: the Northern Atlantic (310°N) and the Mediterranean Sea (150°N). We observed that the ambient seismic noise and the resulting NCF were dominated either by the fundamental mode or by the first overtone of the Rayleigh wave, depending on the period.

We obtained a set of group and phase dispersion maps at periods between 1 and 7 s. The primary feature that emerged from these maps at periods above 2 s was a marked velocity decrease from the edge of the Upper Rhine Graben toward its central axis. We interpreted this gradient as the signature of the dipping crystalline basement. We interpreted shorter scale velocity anomalies observed
<figure><img src="image_15.png" /><figcaption>**Figure 14.** (a) Phase dispersion map obtained in this work at period 2.5 s. The purple line delineates the area, where the network quality exceeds 0.25 (Fig. 6c). (b) Top basement elevation of the area modified after Baillieux (2012). The black line and numbers from 1 to 4 correspond roughly to the section presented on sketch (c). (c) 2-D density model of the area obtained by inverting gravity measurements (modified after Rotstein *et al.* 2006) densities are expressed in g cm⁻³. (d) Temperature anomalies in °C at top basement modified after Baillieux *et al.* (2014). White dots correspond to the seismological stations used in this study. The black lines indicate the local faults at the top of the basement (Baillieux *et al.* 2011; Dezayes *et al.* 2011).</figcaption></figure>

at periods below 2 s as the signatures of heterogeneities inside the sedimentary cover such as local horsts and grabens.

We showed that our approach can be used to include NCFs from station pairs shorter than the empirical limit of 2 or 3 wavelengths in imaging problems without inducing errors in the recovered velocity models. The short station pairs were more sensitive to the azimuthal distribution of the noise intensity than to the velocity anomalies themselves. In order to ensure the usefulness as well as the reliability of tomographic images, we suggest limiting the use of our method to networks that have a significant proportion of station-pairs longer than one wavelength.

The method we propose could be applied to the high density short-period networks that are currently being deployed for passive imaging at local scales (e.g. Lin & Tsai 2013; Hand 2014; Zigone *et al.* 2015; Vergne *et al.* 2015) and should be able to lead to high resolution 3-D models, free from the bias caused by non uniformly distributed seismic noise.

## ACKNOWLEDGEMENTS

This work is supported by the LABEX ANR-11-LABX-0050.G-EAU-THERMIE-PROFONDE and benefits from a funding from the state managed by the French National Research Agency as part of the Investments for the future program. The authors thank GEIE-EMC, ECOGI and EOST for providing the data from the permanent stations. We thank the Geophysical Instrument Pool Potsdam (GFZ) for providing temporary stations of network ‘KIT1’, as well as E. Gaucher (KIT), V. Maurer (ES-Géothermie), H. Wodling, H. Jund and M. Grünberg (EOST) who deployed the stations and collected the data. This work was mainly conducted in the framework of the PhD of M. Lehujeur funded by Groupe Electricité de Strasbourg who also funded the temporary Broad-band stations. The authors thank D. Zigone as well as the two anonymous reviewers, who accepted to review the manuscript and helped to improve it.
REFERENCES

Baillieux, P., 2012. Multidisciplinary approach to understand the localization of geothermal anomalies in the Upper Rhine Graben from regional to local scale, *PhD thesis*, University of Neuchâtel, Neuchâtel.

Baillieux, P., Schill, E. & Dezayes, C., 2011. 3-D structural regional model of the EGS Soultz site (northern Upper Rhine Graben, France): insights and perspectives, in *36th Workshop on Geothermal Reservoir Engineering*, Stanford, CA, p. SGP-TR-191.

Baillieux, P., Schill, E., Abdelfettah, Y. & Dezayes, C., 2014. Possible natural fluid pathways from gravity pseudo-tomography in the geothermal fields of Northern Alsace (Upper Rhine Graben), *Geotherm. Energy*, **2**, doi:10.1186/s40517-014-0016-y.

Baujard, C., 2015. ECOGI, a new deep EGS project in Alsace, Rhine Graben, France, *Presented at the Proceedings World Geothermal Congress*, Melbourne, Australia.

Beauce, A., Fabriol, H., Le Masne, D., Cavoit, C., Mechler, P. & Chen, X.K., 1991. Seismic studies on the HDR site of Soultz-sous-Forêts (Alsace, France), *Geotherm. Sci. Technol.*, **3**, 239–266.

Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M. & Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, *Geophys. J. Int.*, **169**, 1239–1260.

Beucler, É., Mocquet, A., Schimmel, M., Chevrot, S., Quillard, O., Vergne, J. & Sylvander, M., 2015. Observation of deep water microseisms in the North Atlantic Ocean using tide modulations, *Geophys. Res. Lett.*, **42**, 316–322.

Bonnefoy-Claudet, S., Cotton, F. & Bard, P.-Y., 2006. The nature of noise wavefield and its applications for site effects studies: a literature review, *Earth-Sci. Rev.*, **79**, 205–227.

Boué, P., Roux, P., Campillo, M. & Briand, X., 2014. Phase velocity tomography of surface waves using ambient noise cross correlation and array processing, *J. geophys. Res.*, **119**, 519–529.

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

Campillo, M., Sato, H., Shapiro, N.M. & van der Hilst, R.D., 2011. Nouveaux développements de l'imagerie et du suivi temporel à partir du bruit sismique, *C. R. Geosci.*, **343**, 487–495.

Cautru, J.P., 1988. *Coupe géologique passant par le forage GPK-1 calée sur la sismique réflexion*, BRGM.

Cuenot, N., Dorbath, C. & Dorbath, L., 2008. Analysis of the microseismicity induced by fluid injections at the EGS site of Soultz-sous-Forêts (Alsace, France): implications for the characterization of the geothermal reservoir properties, *Pure appl. Geophys.*, **165**, 797–828.

Curtis, A. & Halliday, D., 2010. Directional balancing for seismic and general wavefield interferometry, *Geophysics*, **75**, SA1–SA14.

Dezayes, C., Beccaletto, L., Oliviero, G., Baillieux, P., Capar, L. & Schill, E., 2011. 3-D visualization of a fractured geothermal field; the example of the EGS Soultz site (northern Upper Rhine Graben, France), in *36th Workshop on Geothermal Reservoir Engineering*, Stanford, CA, Vol. 191.

Froment, B., Campillo, M., Roux, P., Gouédard, P., Verdel, A. & Weaver, R.L., 2010. Estimation of the effect of nonisotropically distributed energy on the apparent arrival time in correlations, *Geophysics*, **75**, SA85–SA93.

Genter, A. & Traineau, H., 1996. Analysis of macroscopic fractures in granite in the HDR geothermal well EPS-1, Soultz-sous-Forêts, France, *J. Volcanol. Geotherm. Res.*, **72**, 121–141.

Genter, A., Evans, K., Cuenot, N., Fritsh, D. & Sanjuan, B., 2010. Contribution of the exploration of deep crystalline fractured reservoir of Soultz to the knowledge of enhanced geothermal systems (EGS), *C. R. Geosci.*, **342**, 502–516.

Groos, J.C. & Ritter, J.R.R., 2009. Time domain classification and quantification of seismic noise in an urban environment, *Geophys. J. Int.*, **179**, 1213–1231.

Gutenberg, B., 1936. On microseisms, *Bull. seism. Soc. Am.*, **26**, 111–117.

Hand, E., 2014. A boom in boomless seismology, *Science*, **345**, 720–721.

Herrmann, R.B. & Ammon, C.J., 2002. *Computer Programs in Seismology: Surface Waves, Receiver Functions and Crustal Structure*, St. Louis University.

Kedar, S., Longuet-Higgins, M., Webb, F., Graham, N., Clayton, R. & Jones, C., 2008. The origin of deep ocean microseisms in the North Atlantic Ocean, *Proc. R. Soc. A*, **464**, 777–793.

Lehujeur, M., Vergne, J., Schmittbuhl, J. & Maggi, A., 2015. Characterization of ambient seismic noise near a deep geothermal reservoir and implications for interferometric methods: a case study in northern Alsace, France, *Geotherm. Energy*, **3**(3), doi:10.1186/s40517-014-0020-2.

Lin, F.C., Moschetti, M.P. & Ritzwoller, M.H., 2008. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, *Geophys. J. Int.*, **173**, 281–298.

Lin, F.-C., Ritzwoller, M.H. & Snieder, R., 2009. Eikonal tomography: surface wave tomography by phase front tracking across a regional broadband seismic array, *Geophys. J. Int.*, **177**, 1091–1110.

Lin, F.-C. & Tsai, V.C., 2013. Seismic interferometry with antipodal station pairs, *Geophys. Res. Lett.*, **40**, 4609–4613.

Lobkis, O.I. & Weaver, R.L., 2001. On the emergence of the Green's function in the correlations of a diffuse field, *J. acoust. Soc. Am.*, **110**, 3011–3017.

Longuet-Higgins, M.S., 1950. A theory of the origin of microseisms, *Phil. Trans. R. Soc. Lond.*, **A**, **243**, 1–35.

Luo, Y., Yang, Y., Xu, Y., Xu, H., Zhao, K. & Wang, K., 2015. On the limitations of interstation distances in ambient noise tomography, *Geophys. J. Int.*, **201**, 652–661.

McNamara, D.E. & Buland, R.P., 2004. Ambient noise levels in the continental United States, *Bull. seism. Soc. Am.*, **94**, 1517–1527.

Nakata, N., Chang, J.P., Lawrence, J.F. & Boué, P., 2015. Body wave extraction and tomography at Long Beach, California, with ambient-noise interferometry, *J. geophys. Res. Solid Earth*, **120**, 1159–1173.

Pedersen, H.A. & Krüger, F., 2007. Influence of the seismic noise characteristics on noise correlations in the Baltic shield, *Geophys. J. Int.*, **168**, 197–210.

Poli, P., Pedersen, H.A. & Campillo, M., 2012. Emergence of body waves from cross-correlation of short period seismic noise, *Geophys. J. Int.*, **188**, 549–558.

Renard, P. & Courrioux, G., 1994. Three-dimensional geometric modeling of a faulted domain: the Soultz Horst example (Alsace, France), *Comput. Geosci.*, **20**, 1379–1390.

Rotstein, Y., Edel, J.-B., Gabriel, G., Boulanger, D., Schaming, M. & Munschy, M., 2006. Insight into the structure of the Upper Rhine Graben and its basement from a new compilation of Bouguer Gravity, *Tectonophysics*, **425**, 55–70.

Roux, P., Sabra, K.G., Gerstoft, P., Kuperman, W.A. & Fehler, M.C., 2005a. *P*-waves from cross-correlation of seismic noise, *Geophys. Res. Lett.*, **32**, doi:10.1029/2005GL023803.

Roux, P., Sabra, K.G., Kuperman, W.A. & Roux, A., 2005b. Ambient noise cross correlation in free space: theoretical approach, *J. acoust. Soc. Am.*, **117**, 79–84.

Sabra, K.G., Gerstoft, P., Roux, P., Kuperman, W.A. & Fehler, M.C., 2005. Surface wave tomography from microseisms in Southern California, *Geophys. Res. Lett.*, **32**, doi:10.1029/2005GL023155.

Sergeant, A., Stutzmann, E., Maggi, A., Schimmel, M., Arduin, F. & Obreb-ski, M., 2013. Frequency-dependent noise sources in the North Atlantic Ocean, *Geochem. Geophys. Geosyst.*, **14**, 5341–5353.

Shapiro, N.M. & Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, *Geophys. Res. Lett.*, **31**, doi:10.1029/2004GL019491.

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

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**, L18310, doi:10.1029/2006GL027010.

Snieder, R., 2004. Extracting the Green's function from the correlation of coda waves: a derivation based on stationary phase, *Phys. Rev. E*, **69**, 046610, doi:10.1103/PhysRevE.69.046610.

Spakman, W. & Bijwaard, H., 2001. Optimization of cell parameterizations for tomographic inverse problems, *Pure appl. Geophys.*, **158**, 1401–1423.
Stehly, L., Campillo, M. & Shapiro, N.M., 2006. A study of the seismic noise from its long-range correlation properties, *J. geophys. Res.*, **111**, doi:10.1029/2005JB004237.

Stehly, L., Fry, B., Campillo, M., Shapiro, N.M., Guilbert, J., Boschi, L. & Giardini, D., 2009. Tomography of the Alpine region from observations of seismic ambient noise, *Geophys. J. Int.*, **178**, 338–350.

Stutzmann, E., Roult, G. & Astiz, L., 2000. Geoscope station noise levels, *Bull. seism. Soc. Am.*, **90**(3), 690–701.

Tsai, V.C., 2009. On establishing the accuracy of noise tomography travel-time measurements in a realistic medium, *Geophys. J. Int.*, **178**, 1555–1564.

Vergne, J., Blachet, A. & Lehujeur, M., 2015. Sources of high frequency seismic noise: insights from a dense network of ~250 stations in northern Alsace (France), *Presented at the EGU General Assembly Conference Abstracts*, Vienna, Austria, 9164 pp.

Wapenaar, K., van der Neut, J. & Ruigrok, E., 2008. Passive seismic interferometry by multidimensional deconvolution, *Geophysics*, **73**, A51–A56.

Wapenaar, K., Ruigrok, E., van der Neut, J. & Draganov, D., 2011. Improved surface-wave retrieval from ambient seismic noise by multi-dimensional deconvolution, *Geophys. Res. Lett.*, **38**, L01313, doi:10.1029/2010GL045523.

Weaver, R., Froment, B. & Campillo, M., 2009. On the correlation of non-isotropically distributed ballistic scalar diffuse waves, *J. acoust. Soc. Am.*, **126**, 1817–1826.

Withers, M.M., Aster, R.C., Young, C.J. & Chael, E.P., 1996. High-frequency analysis of seismic background noise as a function of wind speed and shallow depth, *Bull. seism. Soc. Am.*, **86**, 1507–1515.

Yang, Y., Ritzwoller, M.H., Levshin, A.L. & Shapiro, N.M., 2007. Ambient noise Rayleigh wave tomography across Europe, *Geophys. J. Int.*, **168**, 259–274.

Yao, H. & van der Hilst, R.D., 2009. Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to SE Tibet, *Geophys. J. Int.*, **179**, 1113–1132.

Yao, H., Van Der Hilst, R.D. & De Hoop, M.V., 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis – I. Phase velocity maps, *Geophys. J. Int.*, **166**, 732–744.

Young, M.K., Rawlinson, N., Arroucau, P., Reading, A.M. & Tkalčić, H., 2011. High-frequency ambient noise tomography of southeast Australia: new constraints on Tasmania’s tectonic past, *Geophys. Res. Lett.*, **38**, doi:10.1029/2011GL047971.

Zigone, D., Ben-Zion, Y., Campillo, M. & Roux, P., 2015. Seismic tomography of the Southern California plate boundary region from noise-based Rayleigh and Love waves, *Pure appl. Geophys.*, **172**, 1007–1032.

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

**Figure S1.1.** FK analysis of the ambient seismic noise recorded in May 2014 by the available network. Each line corresponds to a given period (T = 10, 5.2 and 3.4 s). The right column is the FK transform of the noise records. The dashed circles indicate various phase speed values: the black dashed circle corresponds to a phase speed of 7 km s⁻¹; the white, red and yellow dashed circles correspond to the fundamental mode and the two first overtones of the Rayleigh waves respectively. Their radius varies with the period and is computed using the local 1-D velocity model. The left column corresponds to the array transfer function displayed at the same zoom levels as the FK transforms.

**Figure S1.2.** Same figure as Fig. S1.1 at periods 2 s (upper line) and 1.5 s (lower line).

**Figure S3.1.** Synthetic and inverted noise correlation functions at period 4.5 s. The NCFs are sorted as a function of the interstation distance and azimuth. (a) Synthetic dataset obtained from the noise

model of Fig. S3.2 (blue curve) and the velocity model of Figs S3.4(a) and (b). A white Gaussian noise is added to the waveforms to mimic measurement inaccuracies. (b, c) Inverted data sets obtained from inversion tests A and B, respectively.

**Figure S3.2.** Synthetic noise model and results of the inversion tests A and B at 2 s of period.

**Figure S3.3.** Inner part of the synthetic dispersion maps and solutions of inversion tests A and B at 2 s of period. (a, b) Synthetic noise phase and group speed dispersion maps. The solid lines correspond to the Voronoï grid. (c, d) Results of inversion A in which the noise model is kept spatially uniform. The maps are smoothed by linear interpolation of the velocity values in each cell of the grid. The contour lines reproduce the expected velocity anomalies from maps (a) and (b) for comparison. (e, f) Results of inversion B. The red dots correspond to the stations. Full-scale maps are presented in Fig. S3.4.

**Figure S3.4.** Full-scale synthetic and inverted dispersion maps at period 2 s.

**Figure S4.1.** Group-velocity measurements used in the ‘standard’ approach at period 3 s. (a) Noise correlation functions between pairs longer than two wavelengths, band-pass filtered at 3 s and sorted by distance. The vertical grey bars indicate the time interval (between 0.5 and 5 km s⁻¹) in which the maximum of the envelope is searched (black dots). The red bars indicate the measurement uncertainties taken as half the difference in time between the causal and acausal travel times. (b) Group-velocity values for each pair obtained using the mean of the causal and acausal arrival times. Red bars indicate the group-velocity uncertainties computed from the time uncertainties of subplot (a). (c) Same data as panel (b) displayed as a function of the pair orientation (measured at the eastern station, clockwise from north).

**Figure S4.2.** Comparison between the group dispersion maps obtained at period 3 s with the two-step inversion procedure that takes the noise distribution into account (a) and with a more standard approach (b). Group-velocity variations are given in per cent relative to a reference velocity of 1.81 km s⁻¹. Grey paths indicate the station pairs used in the standard approach. The black ellipse indicate a positive velocity anomaly interpreted as a probable azimuthal effect caused by a local dominance of station pairs orthogonal to the main noise direction (~310°N).

(http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggw373/-/DC1).

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.

## APPENDIX A: MODELLING THE NOISE CORRELATION FUNCTIONS

We consider that the ambient seismic noise can be modelled as a superposition of noise sources distributed randomly in space and time. This conceptualization of the noise has been used to establish the link between the cross-correlation function and the Green’s function if the spatial distribution of the noise sources is uniform (Snieder 2004; Roux et al. 2005). Here we assume that the noise sources are located in the horizontal plane, so that the vertical displacement can be written in the frequency domain as

$$U(\mathbf{r}, \nu) = \int_{\mathbb{R}^2} S(\mathbf{r}_s, \nu)G(\mathbf{r}, \mathbf{r}_s, \nu)d\mathbf{r}_s \quad (A1)$$
With $v$ the frequency; $r_S$ the location of the source in the 2-D plane; $S(r_S, v)$ the source term and $G(r, r_S, v)$ the Green's function between locations $r$ and $r_S$. Following Roux et al. (2005) we call A and B two stations located on the surface at Cartesian coordinates $r_A = (+a, 0)$ and $r_B = (-a, 0)$. The cross-correlation function between stations A and B in the frequency domain is

$$C_{AB}(v) = U(r_A, v) U^*(r_B, v) \qquad (A2)$$

where * denotes the complex conjugation. We apply it to the signal recorded by station B (i.e. the western station) for consistency with the orientation convention used by Stehly et al. (2006), which implies that a wave that propagates toward the east emerges on the causal part of the cross-correlation function. Under the assumption of temporally uncorrelated noise sources (Snieder 2004; Roux et al. 2005), the cross-correlation of distinct source functions over infinite times is zero, while their autocorrelation gives the power spectral density of the source term at a given location, which can be considered as the seismic energy delivered at that point:

$$ \begin{aligned} S(r_S, v) S(r'_S, v)^* &= |S(r_S, v)|^2 \delta(r_S - r'_S) \\ &= E(r_S, v) \delta(r_S - r'_S) \end{aligned} \qquad (A3) $$

where $\delta$ is the Dirac function and $E(r_S, v)$ is the spatial distribution of the noise energy. Combining eqs (A1) and (A3) we obtain:

$$C_{AB}(v) = \int_{\mathbb{R}^2} E(r_S, v) G(r_A, r_S, v) G^*(r_B, r_S, v) dr_S \qquad (A4)$$

The elliptical coordinate system is especially suited to that problem because the sources that affect a given lag-time of the correlation functions are located on a hyperbola whose focal points are the two stations (Roux et al. 2005). In the 2-D case, we use $x = a \cos \xi \cosh \mu$, $y = a \sin \xi \sinh \mu$ with $\xi$ ranging by convention from 0 to $\pi$ and $\mu$ ranging from $-\infty$ to $+\infty$ (see Fig. A1). In that system, A4 becomes

$$ C_{AB}(v) = \int_{-\infty}^{+\infty} \int_{0}^{\pi} E(\xi, \mu, v) G(r_A, r_S, v) G^*(r_B, r_S, v) \times J(\xi, \mu) d\xi d\mu \qquad (A4') $$

where $J(\xi, \mu) = a^2 (\cosh^2\mu - \cos^2\xi)$ is the Jacobian of the coordinate system transformation.

One can show the following relations:

$$ \begin{cases} SA - SB = -2a \cos \xi \\ SA + SB = 2a \cosh \mu \\ SA \times SB = a^2 (\cosh^2 \mu - \cos^2 \xi) = J(\xi, \mu) \end{cases} \qquad (A5) $$

with $SA$ the distance between source $S$ and receiver $A$. Using eq. (1), we obtain

$$ \begin{aligned} G(r_A, r_S, v) G^*(r_B, r_S, v) &= \frac{e^{-\gamma(v)(SA+SB)}}{2\pi h(v) \sqrt{SA \times SB}} e^{-2i\pi v \frac{(SA-SB)}{c(v)}} \\ &= \frac{e^{-2a\gamma(v)\cosh\mu} e^{4i\pi v a \frac{\cos\xi}{c(v)}}}{2\pi h(v) \sqrt{J(\xi, \mu)}} \qquad (A6) \end{aligned} $$

After substituting this result into eq. (A4'), we obtain the cross correlation function between receivers A and B in the frequency domain for any source distribution:

$$ C_{AB}(v) = \frac{1}{2\pi h(v)} \int_{0}^{\pi} \mathcal{H}(\xi, v) e^{4i\pi v a \frac{\cos\xi}{c(v)}} d\xi \qquad (A7) $$

<figure><img src="image_16.png" /><figcaption>Figure A1. Coordinate systems used to locate the noise sources. Sources (red and green stars) are disposed on an ellipse (thick black line) of parameter $\mu_0 \gg 0$ that can be approached by a circle of radius $r_0$ (red dashed circle). The station pair (AB) is oriented with an angle $\alpha$ with respect to the north. The source location is expressed relative to the station pair using the elliptical coordinate system ($\xi, \mu_0$) or relative to the north using the polar coordinate system ($r_0, \theta$). Blue lines indicate constant values of $\mu$ (ellipses) or $\xi$ (hyperbolas).</figcaption></figure>

where $\mathcal{H}$ corresponds to the integrated amount of energy contributed by the sources located on a hyperbola of parameter $\xi$:

$$ \mathcal{H}(\xi, v) = \int_{\mathbb{R}} E(\xi, \mu, v) e^{-2a\gamma(v)\cosh\mu} \sqrt{J(\xi, \mu)} d\mu. \qquad (A8) $$

The case of a uniform noise distribution can be modelled by taking $E(\xi, \mu, v)$ constant relative to the spatial coordinates $\xi, \mu$. In such case, function $\mathcal{H}(\xi, v)$ is smooth over $\xi$, which has motivated the use of the stationary phase integration technique to establish the well known link between the derivative of the correlation function and the interstation Green's function (Snieder 2004). This technique is not used in the present paper and eq. (A7) is integrated numerically to account for the $\xi$-dependency of function $\mathcal{H}$.

In this work, we assume that the noise sources are located far from the network so that the noise distribution can be reduced to a 1-D azimuthal function. Such assumption implies that the distribution of the noise energy has the following form

$$ E(\xi, \mu, v) = \tilde{E}(\xi, v) \delta(\mu - \mu_0) + \tilde{E}(-\xi, v) \delta(\mu + \mu_0) \quad \mu_0 \gg 0 \qquad (A9) $$

Where $\tilde{E}$ represents the azimuthal variation of the noise energy relative to the AB station pair [similar to the function referred as B by Weaver et al. (2009) and Froment et al. (2010)]. Eq. (A9) assumes that the noise distribution can be modelled as a set of effective sources disposed along an ellipse with constant parameter $\mu_0$. Because $\mu_0$ is taken very large, this ellipse can be approached by a circle (see Fig. A1) with radius $r_0 = a \cosh \mu_0 = a \sinh \mu_0$ (because $\lim_{\mu \to \infty} \frac{\cosh \mu}{\sinh \mu} = 1$) and $r_0 \gg a$.
<figure><img src="image_17.png" /><figcaption>**Figure B1.** Grid search on the regularization parameters for inversion step 2 at 2.5 s of period. $\sigma_c$ and $\sigma_l$ denote the uncertainties on the prior dispersion parameters $c^{\text{prior}}$ and $l^{\text{prior}}$ imposed according to the results of inversion step 1. $\chi_{\text{red}}^2$ represents the data misfit, $||c - c^{\text{prior}}||^2$ and $||l - l^{\text{prior}}||^2$ represent the misfit between the inverted and prior dispersion maps. (a) $\sigma_l$ is fixed to 0 and several values of $\sigma_c$ are tested from 0 (infinite constraint on the a priori phase dispersion map) to infinity (no constraint on the phase dispersion map). (b) Once the best value of $\sigma_c$ is determined (0.03), $\sigma_c$ is fixed and $\sigma_l$ is perturbed from 0 to infinity.</figcaption></figure>

Inserting (A9) into (A8), we get

$$ \mathcal{H}(\xi, \nu) = (\tilde{E}(\xi, \nu) + \tilde{E}(-\xi, \nu)) e^{-2a\gamma(\nu)\cosh \mu_0} \sqrt{J(\xi, \mu_0)} \\ = (\tilde{E}(\xi, \nu) + \tilde{E}(-\xi, \nu)) e^{-2r_0\gamma(\nu)} \sqrt{r_0^2 - a^2\cos^2\xi}. \qquad (A10) $$

Noting that $\sqrt{r_0^2 - a^2\cos^2\xi} = r_0 - \frac{a^2\cos^2\xi}{2r_0} + o(\frac{1}{r_0^2}) \approx r_0$, we get

$$ \mathcal{H}(\xi, \nu) \approx (\tilde{E}(\xi, \nu) + \tilde{E}(-\xi, \nu)) r_0 e^{-2r_0\gamma(\nu)} \quad \text{if } r_0 \gg a. \qquad (A11) $$

Eq. (A11) shows that in the specific case where the same amount of energy reaches the network from all directions (i.e. $\tilde{E}(\xi, \nu) = \tilde{E}(-\xi, \nu) = E_0(\nu)$), then $\mathcal{H}(\xi, \nu) \approx 2E_0(\nu)r_0e^{-2r_0\gamma(\nu)}$ suggesting that the azimuthal variation of $\mathcal{H}$ can be neglected at first order when $\tilde{E}$ is constant and when the sources are very far from the network. The assumption of an azimuthally uniform noise distribution was therefore modeled using $\mathcal{H}$ constant over $\xi$ (Fig. 4, black curve).

A discretized version of function $\tilde{E}$ might be used to parametrize the azimuthal distribution of the noise energy relative to the station pair. Instead, we use a north reference system and introduce function $A(\theta, \nu)$ that corresponds to the amount of energy received at frequency $\nu$ from back-azimuth $\theta$ counted clockwise from north. The link between functions $A$ and $\tilde{E}$ for a station pair that is oriented with an angle $\alpha$ measured at the western station clockwise from north is given by

$$ \begin{cases} A(\alpha - \xi, \nu) = \tilde{E}(+\xi, \nu)r_0e^{-2r_0\gamma(\nu)} \\ A(\alpha + \xi, \nu) = \tilde{E}(-\xi, \nu)r_0e^{-2r_0\gamma(\nu)} \end{cases} \qquad (A12) $$

Combining (A11) and (A12), we get eq. (3) that provides us the link between the north referenced distribution of the noise energy (A) and the hyperbolic contribution $\mathcal{H}_\alpha$ for each station pair.

## APPENDIX B: REGULARIZATION OF INVERSION (STEP 2) OF OBSERVED DATA

Inversion step 2 is regularized by adjusting the uncertainty on the phase speed parameters ($\sigma_c$) and on the first derivative of the group dispersion curve ($\sigma_l$). We perform a grid search on these two parameters to minimize both the data misfit ($\chi_{\text{red}}^2$) and the misfits between the inverted and the prior dispersion maps ($||c - c^{\text{prior}}||^2$ and $||l - l^{\text{prior}}||^2$) defined as:

$$ \left\{ \begin{aligned} \chi_{\text{red}}^2 &= \frac{1}{l} \sum_{i=0}^{l-1} \left( \frac{g(m_i) - d_i^{\text{obs}}}{\sigma_l} \right)^2 \\ ||c - c^{\text{prior}}||^2 &= \sum_{j=0}^{N-1} \frac{S_j}{S} (c_j - c_j^{\text{prior}})^2 \\ ||l - l^{\text{prior}}||^2 &= \sum_{j=0}^{N-1} \frac{S_j}{S} (l_j - l_j^{\text{prior}})^2 \end{aligned} \right. \qquad (B1) $$

The grid search is done in two stages: first, the $\sigma_l$ value is fixed to 0 (i.e. the first derivative of the dispersion curve in each cell $j$ is constrained to the prior value $l_j^{\text{prior}}$) and we look for the best $\sigma_c$ value. Second, $\sigma_c$ is fixed and we look for the best value of $\sigma_l$. At 2.5 s of period, the regularization coefficients retained are $\sigma_c = 0.03$ and $\sigma_l = 0.09$ (Fig. B1).