# Investigating the lateral resolution of the Rayleigh wave focal spot imaging technique using two-dimensional acoustic simulations

Bruno Giammarinaro, Christina Tsarsitalidou, Gregor Hillers

► To cite this version:

Bruno Giammarinaro, Christina Tsarsitalidou, Gregor Hillers. Investigating the lateral resolution of the Rayleigh wave focal spot imaging technique using two-dimensional acoustic simulations. Comptes Rendus. Géoscience, 2024. hal-04155285v2

HAL Id: hal-04155285

https://hal.science/hal-04155285v2

Submitted on 29 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# Comptes Rendus
## Géoscience
### Sciences de la Planète

Bruno Giammarinaro, Christina Tsarsitalidou and Gregor Hillers

**Investigating the lateral resolution of the Rayleigh wave focal spot imaging technique using two-dimensional acoustic simulations**

Published online: 22 February 2024

**Part of Special Issue:** New Developments in Passive Seismic Imaging and Monitoring

**Guest editors:** Michel Campillo (University Grenoble-Alpes, Institut des Sciences de la Terre, Grenoble, France),

Andrew Curtis (School of GeoSciences, University of Edinburgh, Scotland),

Anne Obermann (Swiss Seismological Service, ETH, Zurich, Switzerland) and Nikolai

Shapiro (CNRS, Institut des Sciences de la Terre, Grenoble, France)

https://doi.org/10.5802/crgeos.254

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Research article

New Developments in Passive Seismic Imaging and Monitoring

# Investigating the lateral resolution of the Rayleigh wave focal spot imaging technique using two-dimensional acoustic simulations

Bruno Giammarinaro<sup>®,*,a</sup>, Christina Tsarsitalidou<sup>®,a</sup> and Gregor Hillers<sup>®,a</sup>

<sup>a</sup> Institute of Seismology, University of Helsinki, Helsinki, Finland

*Current address:* LabTAU, INSERM, Centre Léon Bérard, Université Lyon 1, Lyon, France (B. Giammarinaro)

*E-mails:* bruno.giammarinaro@inserm.fr (B. Giammarinaro), christina.tsarsitalidou@helsinki.fi (C. Tsarsitalidou), gregor.hillers@helsinki.fi (G. Hillers)

**Abstract.** We investigate the lateral resolution power of the seismic Rayleigh wave focal spot imaging technique. We use two-dimensional acoustics simulations in a closed cavity for the passive Green's function and focal spot reconstruction. Four different velocity distributions target different resolution aspects. The finite data range that is necessary to constrain the Bessel function model controls the lateral spreading of material contrasts, the distinction of two objects on sub-wavelength scales, and the image quality of complex random media. Good data quality from dense networks supports short range estimates and super-resolution.

**Keywords.** Numerical modeling, Wave propagation, Seismic noise, Cavity, Focal spot, Resolution, Seismic imaging.

**Funding.** Academy of Finland (decision number 322421).

Manuscript received 8 December 2022, revised 5 December 2023 and 9 February 2024, accepted 29 January 2024.

## 1. Introduction

Disadvantages of earthquake tomography associated with limited illumination can now be compensated by ambient noise tomography with its flexible virtual source and receiver configurations. Both approaches invert far-field observations of travel time differences, obtained from earthquake seismograms [Tromp et al., 2005, Liu and Gu, 2012] or from passive Green's function reconstructions [Shapiro and Campillo, 2004, Sabra et al., 2005], for a model of the velocity structure.

Modern dense seismic arrays support alternative local surface wave speed estimations from noise correlation functions [Lin et al., 2009, Lin and Ritzwoller, 2011], which includes the large scale application of the frequency domain spatial autocorrelation (SPAC) method [Aki, 1957, Ekström et al., 2009, Ekström, 2014] that is otherwise typically applied to local sparse array data [Asten, 2006]. Dense arrays can now contain on the order of 1000 sensors, which facilitates the proper sampling of the noise correlation amplitude distribution at near-field distances. At zero lag time, the time domain representation of the spatial autocorrelation field is referred to as focal spot, which contains the same information as SPAC

<footnote>* Corresponding author.</footnote>
and can be analyzed using the same mathematical tools [Cox, 1973; Yokoi and Margaryan, 2008, Tsai and Moschetti, 2010, Haney et al., 2012, Haney and Nakahara, 2014].

Focal spot analysis is used in applications that work with a high sensor density. Focal spots have first been studied in time-reversal experiments in underwater acoustics and medical imaging [Fink, 1997]. Noise correlation medical imaging transferred the noise seismology approaches to passive elastography, first using far-field surface waves along muscle fibers [Sabra et al., 2007], later using refocusing shear waves. Different properties of the shear wave focal spot have been analyzed including its cross-section width [Catheline et al., 2008, Gallot et al., 2011], two-dimensional shape [Benech et al., 2013, Brum et al., 2015], and its tip curvature [Catheline et al., 2013], which can all be reconstructed using MRI [Zorgani et al., 2015] and ultrasound [Barrere et al., 2020] speckle tracking methods. Importantly, Zemzemi et al. [2020] demonstrated that the ability to discriminate two objects is not controlled and hence limited by the shear wavelength, but instead by the ultrasonic frequency and the pixel density. In seismology, this corresponds to the array station density.

Hillers et al. [2016] first applied Rayleigh wave focal spot imaging in seismology to image lateral velocity variations in a fault zone environment. As for SPAC, the phase velocity is estimated from the focal spot shape using Bessel function models. Using numerical time-reversal experiments based on a Green's function calculator for one-dimensional layered media [Cotton and Coutant, 1997], Giammarinaro et al. [2023] demonstrated the feasibility to accurately estimate phase velocity and dispersion from noisy focal spots. Results for non-isotropic surface wave illumination showed that the wave speed bias is negligible if the sensors are isotropically distributed, compatible with the SPAC results by Nakahara [2006]. Moreover, the effect of interfering P waves can be mitigated using a data or fitting range around one wavelength. Whereas these results from Giammarinaro et al. [2023] suggest the overall robustness and utility of the focal spot method for seismic imaging applications, most notably because of the increase in depth resolution, the study could not address lateral resolution. Considering that super-resolution can be obtained with tip-curvature measurements [Zemzemi et al., 2020], it is important to assess the effect of the

data range on speed estimates from densely sampled seismic focal spots.

Here we study systematically the lateral focal spot resolution using numerical experiments. We perform two-dimensional acoustics simulations to reconstruct the Green's function from reverberating wave fields (Section 2). The ambient field generated in a chaotic closed cavity [Draeger and Fink, 1999] yields results that are equivalent to results from open media noise correlation [Derode et al., 2003]. The obtained Green's function is identical and is here therefore taken as a proxy for seismic vertical-vertical component Rayleigh wave correlations [Sanchez-Sesma, 2006, Haney et al., 2012]. We work with a constant number of grid points and a fixed reference frequency. We implement four test cases and vary the data range or fitting distance $r_{fit}$ to investigate the effect on the resolution of the velocity structure. These cases include a homogeneous control experiment (Section 3.1), an interface between two half-spaces (Section 3.2), circular inclusions (Section 3.3), and heterogeneous or random velocity distributions (Section 3.4). In Section 4 we discuss the different resolution aspects that are investigated with the variable configurations for a comprehensive evaluation of the seismic Rayleigh wave focal spot imaging performance.

## 2. Method

### 2.1. Synthetic experiments

This study is based on synthetic diffuse wave fields generated in a closed cavity. The resulting correlation functions are equivalent to noise correlations in open media [Derode et al., 2003]. Simulations are performed using the function *kspaceFirstOrder2D* from the MATLAB toolbox *kWave* [Treeby et al., 2018]. This function solves a system of first-order acoustic equations for the conservation of mass and momentum using a wavenumber *k*-space pseudospectral method. The two-dimensional medium is composed of 500 × 500 grid points that are spaced in the *x* and *y* direction by $dx = dy = 0.1$ km (Figure 1a). The background wave speed in the cavity is $V_0 = 2$ km/s. The closed cavity is implemented using different densities outside ($\rho_{out} = 59$ kg/m³) and inside ($\rho_{in} = 2950$ kg/m³) the cavity. Choosing $\rho_{out}$ to be 2% of $\rho_{in}$ creates strongly reflecting boundaries from
impedance contrast. It supports a homogeneous stability regime of the kWave simulations through a constant Courant-Friedrichs-Lewy number. This number depends on the wave speed, hence varying only the density mitigates potentially problematic stability conditions associated with a large change in wave speed. We cannot exclude the occurrence of weak numerical dispersion in some of the heterogeneous case studies, but the overall consistency of the synthesized Green's functions and focal spots in the different experiments suggests that this effect does not govern our results. Inside the cavity we select a square target domain consisting of 151 × 151 grid points where we record the solution. In this region we define the different velocity distributions introduced in Section 2.2. Results for each of the four cases discussed in Section 3 are obtained by averaging 11 independent wave field simulations. Each simulation starts with a point source at a different position inside the cavity (Figure 1a). The source term is defined as a time-varying mass source (Figures 1b, c) emitting a 1 s long pulse. The pressure is a differentiation of this signal with the frequency range centered at 1 Hz. The wave field is recorded for 300 s inside the square sub-domain with a 100 Hz sampling frequency.

## 2.2. The acoustic medium in the target domain

The acoustic wave speed distribution $V(\mathbf{x})$ is defined as a spatial function in the square target domain

$$V(\mathbf{x}) = V_0(1 + \xi(\mathbf{x})), \quad (1)$$

where $V_0 = 2$ km/s is the background wave speed in the cavity, $\mathbf{x}$ is the position, and $\xi$ is the relative change in wave speed, i.e., the parameter that controls the medium heterogeneity. In the following sections, the 2 km background wavelength at 1 Hz is denoted $\lambda_0$. In Section 3 we begin with a control experiment of a homogeneous medium with $\xi = 0$ to study the overall system response.

## 2.3. Lateral spreading across two welded half-spaces

We modify the homogeneous control experiment and replace the 2 km/s velocity in the left half of the target domain by an increased 2.2 km/s value. This creates two half-spaces and allows us to investigate the lateral resolution as the imaging method induces spreading or averaging across the sharp

interface. Such sharp lateral velocity contrasts can occur across bimaterial interfaces in fault zone environments [Weertman, 1980, Ben-Zion, 1989], or in the contact region of intrusions with the host rock [Chamarczuk et al., 2019].

## 2.4. Resolution of circular inclusions

Next we perform a classic resolution test and study the power of the method to separate two individual entities in an image. For this we impose three pairs of circular inclusions separated by $0.25\lambda_0$, $0.5\lambda_0$ and $1\lambda_0$. The inclusions have a diameter of $1\lambda_0$. We test two different sets, where one set of inclusions is stiffer than the background with $\xi = 25\%$, and the other set has more compliant inclusions compared to the background with $\xi = -25\%$. Such a “two-body problem” is typically studied in gel phantom experiments performed in medical ultrasound imaging or medical imaging for tumor detection [Catheline et al., 2013, Zemzemi et al., 2020]. It is a less common configuration in seismology where so-called checker-board tests are typically employed to quantify the resolution of a tomography configuration. Circular cross-section features can occur in the context of magmatic intrusions or conduits. However, as said, with this configuration we can study the lateral resolution defined as the minimum distance between objects that the method allows to discriminate in an image. This and the other cases are studied using different data ranges $r_{fit}$ discussed in Section 2.6.

## 2.5. Randomly distributed wave speeds

Last we consider random media using a functional form that is often used to parameterize the heterogeneous distributions of variables such as wave speed, stress, or frictional properties in Earth materials [Frankel and Clayton, 1986, Holliger and Levanader, 1992, Mai and Beroza, 2002, Ripperger et al., 2007, Hillers et al., 2007, Sato et al., 2012, Obermann et al., 2016]. We define $\xi(\mathbf{x})$ through a spatial 2D inverse Fourier transform as

$$\xi(\mathbf{x}) = \mathrm{FT}^{-1}\left[\sqrt{P(\mathbf{k})}e^{i\phi(\mathbf{k})}\right], \quad (2)$$

where $\mathbf{k}$ is the spatial wavenumber of the 2D distribution, $P(\mathbf{k})$ is the power spectral density, and $\phi(\mathbf{k})$ is a random distribution of the phase between 0 and $2\pi$. The random distribution is calculated using
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Configuration of the numerical experiments. (a) Representation of the closed cavity with a background velocity $V_0 = 2$ km/s. The gray area represents the outer part of the medium with the same velocity but with an impedance contrast to trap the waves in the cavity. The red line is the boundary of the closed cavity. The black square indicates the target area where the results are recorded. Inside the cavity the color corresponds to the input wave speed. Red dots indicate the source positions for the 11 realizations. (b) Times series and (c) normalized power spectrum of the emitted pulse used for each simulation. (d) Example of a full time-series recorded in the cavity.</figcaption></figure>

the Python `randn` function coupled with a seed fixed to 3. Fixing the seed allows to randomly generate the phase and to keep the same distribution to observe the effect of the control parameters. The power spectral density $P(k)$ follows the von Karman probability function [Sato et al., 2012]

$$P(k) = \frac{4\pi\Gamma[\kappa + 1]\epsilon^2 a^2}{\Gamma[\kappa](1 + a^2\|k\|^2)^{\kappa+1}} \quad (3)$$

The correlation length of the modeled parameter is $a$, $\epsilon$ governs the contrast in the medium, $\kappa$ defines the sharpness of the spectral decay, and $\Gamma$ is the Gamma function. We use Equation (3) to generate nine different media with variable $a$ and $\kappa$ (Table 1), and constant $\epsilon = 150$ km$^{-1}$.

## 2.6. Data processing and wave speed estimation

Data processing is performed in a Python3.8 environment. For each simulation, we analyze the diffusive parts of the wave field (Figure 2a) by focusing on simulated time series data between 100 s and 300 s (Figure 1d). Considering the distances between

the sources and chaotic cavity borders, this corresponds to a range between 10 and 30 mean-free times. We follow Giammarinaro et al. [2023] and filter the traces with a Gaussian filter centered at 1 Hz and with a width of 3.2% to estimate the central frequency dependent phase velocity. In the limit, the analysis could be performed for a harmonic case, however, the narrow filter stabilizes the results. As demonstrated by Giammarinaro et al. [2023], a wider frequency-band filter leads to group velocity estimates, and it biases the focal spot reconstruction by averaging over frequencies which leads to an apparent attenuation. We compute the normalized cross-correlation between each sensor pair to extract the spatial autocorrelation amplitude fields at zero lag time. Results are stacked over the 11 realizations for each case. This yields at short distances around a reference station the large-amplitude feature referred to as focal spot (Figure 2c). The shape is linked to the wave speed through the imaginary part of the Green's function. From each focal spot we estimate the free parameter wave speed $V = \omega/k$ by fitting
<table><caption>**Table 1.** Parameters for the generation of different heterogeneous media using Equation (3)</caption><thead><tr><th>Medium</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><th><i>a</i> (km)</th><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>3</td><td>3</td><td>3</td></tr><tr><th><i>κ</i></th><td>0.1</td><td>0.3</td><td>0.6</td><td>0.1</td><td>0.3</td><td>0.6</td><td>0.1</td><td>0.3</td><td>0.6</td></tr></tbody></table>

Where *a* is the correlation length and *κ* governs the sharpness of the spectral decay. The corresponding wave speed distributions are displayed in the top row of Figure 8.

the azimuthal averaged and normalized data to the $J_0(kr)$ model using a nonlinear least squares regression algorithm [Hillers et al., 2016, Giammarinaro and Hillers, 2022, Giammarinaro et al., 2023]. Again, the 2D acoustic configuration yields results that are equivalent to the lateral propagation of Rayleigh surface waves. The $J_0(kr)$ model equally describes the vertical–vertical component of the Rayleigh wave focal spot [Haney et al., 2012, Haney and Nakahara, 2014], and *V* is thus equivalent to the Rayleigh wave phase velocity $c_R$.

This process is performed for the 22,801 focal spots, and each obtained *V* estimate is associated with the location of the reference station. This imaging concept thus compiles velocity distributions across dense arrays without solving a tomography inverse problem. Importantly, we choose two different data ranges $r_{fit}$ of 1 km and 2 km associated with $0.5\lambda_0$ and $1\lambda_0$ at 1 Hz (Figure 2d). Away from edges, this corresponds to 80 and 314 samples, respectively. We vary $r_{fit}$ because it is a critical tuning parameter, and values around one wavelength yield overall stable results [Giammarinaro et al., 2023]. Larger values can stabilize a regression for noisy signals, but the short distance focal spot imaging concept essentially invites the limitation of $r_{fit}$ for improved resolution. This refers to improved depth resolution as $c_R$ values can be estimated at wavelengths that cannot be studied with tomography [Tsarsitalidou et al., 2021, Giammarinaro et al., 2023], but also to the lateral resolution investigated here.

## 2.7. Error estimation

An important advantage of the focal spot approach is the local assessment of the measurement uncertainty *c* that is estimated following Aster et al. [2019]

$$ \epsilon_k = \sqrt{\frac{RSS}{N-2}} C_k, \quad (4) $$

where RSS is the residual sum of squares, *N* is the number of samples, and $C_k$ is the diagonal element of the parameter covariance matrix *C* associated with the wavenumber *k*. Rules of error propagation apply to yield the wave speed error $\epsilon_V$

$$ \epsilon_V = V\epsilon_k/k. \quad (5) $$

In addition to this estimate we can compare imaged with ground truth input values to discuss resolution.

# 3. Results

## 3.1. The homogeneous reference case

The first test is the homogeneous control experiment which illustrates basic features of the approach. The simulations yield diffuse wave fields that can be used for focal spot imaging. Figure 2(a) shows a snapshot of a diffuse wave field, Figure 2(b) a time-space representation of the time-reversed correlation wave field, Figure 2(c) shows a focal spot of one realization, and Figure 2(d) shows the results of the nonlinear regression following the data processing described in Section 2.6. In Figure 2(b) we can see the refocusing and diverging waves of the Green's function after cross-correlation. However, the reconstruction is not perfect which is indicated by the small-amplitude fluctuations. The time domain autocorrelation field shows the focal spot at small distances around the origin (Figure 2c). The irregularity of the white zero-crossing contours again illustrates fluctuations and imperfect reconstruction after stacking over 11 realizations. We attribute this to non-perfectly diffuse wave fields associated with the modes in the cavity. However, studies of non-diffuse wave fields [Nakahara, 2006, Giammarinaro et al., 2023] demonstrate that isotropically distributed sensors yield unbiased wave speed estimates. This condition is not met along the boundaries of the domain which yields
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Workflow stage examples for the homogeneous medium in the target area. (a) Snapshot of a diffuse wave field. (b) Time-reversed space-time wave field obtained by cross-correlation. (c) Spatial autocorrelation with the focal spot at the center obtained at 1 Hz after stacking over 11 realizations. (d) Simulated 1 Hz focal spot data (gray) and nonlinear regression results for $r_{\text{fit}} = 0.5\lambda_0$ (orange) and $r_{\text{fit}} = 1\lambda_0$ (red). In this and all subsequent figures the unit $\lambda$ is equal to the 1 Hz wavelength $\lambda_0$ for the reference $V_0 = 2$ km/s.</figcaption></figure>

larger fluctuations in the estimates. Figure 2(d) displays results from the nonlinear regression using the two data distances $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$. The spread in the data with increasing distance corresponds to the fluctuations from Figure 2(c).

The 2D velocity distributions obtained with $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$ are displayed in Figures 3(a) and (b), respectively. The key feature in these images are the speckle patterns. They illustrate that the imperfect reconstruction affects the measurement process to produce these fluctuations around the average reference value. Using $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$ leads to $V = 2.014 \pm 0.034$ km/s and $V = 2.008 \pm 0.021$ km/s, respectively, so the $V_0 = 2$ km/s reference value is well recovered with an uncertainty level of 1.5% and 1%. Figures 3(c) and (d) display the spatial error distribution obtained from

Equation (5), where the values indicate a similar error range. Both illustrations demonstrate that an increasing data range improves accuracy and precision. Another strong feature are the boundary effects in Figures 3(a) and (b). Recall that we only use observations inside the square target area, similar to an array deployment. The estimation error increases towards the boundaries and the affected area appears to depend on the $r_{\text{fit}}$ length. However, the biasing effects are not homogeneously distributed along the boundaries but are largest along the southern edge. This is likely explained by the centered, low position of the target area in the cavity, together with the excitation of specific modes in the cavity that sustain a non-isotropic energy flux in the scattered wave field. This also explains the corresponding pattern of larger uncertainties in the lower half of Figure 3(c).
The boundary effects emerge in all our discussed cases but do not affect our conclusions. These effects can be mitigated by a different configuration or geometry, or by averaging over more sources located at more diverse locations.

## 3.2. Resolution of the interface between half-spaces

Figure 4 displays the input wave speed distribution (Figure 4a) and the focal spot-based images of two homogeneous half-space media for the data ranges $r_{\text{fit}} = 0.5\lambda_0$ (Figure 4c) and $r_{\text{fit}} = 1\lambda_0$ (Figure 4d). As for the control experiment, the distributions show small residual fluctuations around the well resolved input values (the distributions in the left half in panels (c) and (d) yield $V = 2.200 \pm 0.099$ km/s and $V = 2.196 \pm 0.095$ km/s, respectively), and again larger edge effects along the lower boundary.

More interesting are the profiles across the domain shown in Figure 4(b). These profiles are averages along the N-axis. The imaged distributions across the interface do not follow the blue input step function but are spread out. To quantify the resolution we estimate the width $\Delta L$ of the transition. This is the distance between the samples where the amplitude of the estimated profile equals the 5% and 95% values of the 0.2 km/s velocity jump, i.e., when the values are 2.01 km/s and 2.19 km/s. Figure 5(a) enlarges the area around the interface located at $3.75\lambda_0$. It shows four profiles obtained with $r_{\text{fit}}$ values ranging from $0.5\lambda_0$ to $2\lambda_0$ in $0.5\lambda_0$ intervals. It confirms the previous observation that the overall velocity in each half-space is correctly estimated for every data range. As for the results in Figures 3 and 4 we can perhaps discern a weak trend to overestimate the reference values. We attribute this to the imperfect simulation configuration since similar effects are not observed in numerical time-reversal experiments, even in the presence of noise [Giammarinaro et al., 2023]. The width $\Delta L$ is shown in Figure 5(a) as dotted lines at the bottom, with the vertical dashed lines being the 5% and 95% boundaries of the transition zone. The calculated widths are compiled in Figure 5(b). For $r_{\text{fit}} = 0.5\lambda_0$, the transition zone has a width of $\Delta L = 0.4\lambda_0$. For $r_{\text{fit}} = 1\lambda_0$ we obtain $\Delta L = 1.2\lambda_0$, and for $r_{\text{fit}} = 2\lambda_0$, it is $\Delta L = 2.2\lambda_0$. The spreading effect quantified as the transition width between the two media can hence well be approximated to scale with

the data range, $\Delta L \approx r_{\text{fit}}$. Figure 5(c) shows the average $r_{\text{fit}}$ dependent $\epsilon_V$ profiles. The lack of features at the position of the interface indicates that the 10% velocity contrast does not influence the reconstruction. The profiles indicate the tendency of an overall better reconstruction in the right slower domain. This is likely controlled by the proportionally smaller data range in the left faster domain, which is a consequence of the constant $r_{\text{fit}}$ value tied to the reference wavelength in the right domain. As expected, this trend weakens with increasing $r_{\text{fit}}$. Figures 5(a) and (c) together show that an increasing $r_{\text{fit}}$ leads to a generally increased precision as $\epsilon_V$ decreases, but to a loss in accuracy around the position of the interface.

## 3.3. Resolution and contrast of circular inclusions

We now examine the resolution of pairs of $1\lambda_0$-wide circular inclusions separated by variable distances. The velocity in the inclusions increases and decreases by 0.5 km/s with respect to the background reference $V_0 = 2$ km/s. Nonlinear regression of the focal spot data are performed using data ranges $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$. Figure 6 collects in the left two columns the input wave speed distributions, and the focal spot obtained images for $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$. As a general observation, every inclusion is visible on the images for both data ranges, and for the three different distances separating the inclusions. This highlights the effectiveness of the Rayleigh wave focal spot method for “discovery mode” applications [Tsai, 2023]. However, the contrast—the difference in “velocity amplitude”—depends on the data range. Using $r_{\text{fit}} = 1\lambda_0$ increases the area of the biased velocity estimates. The quantitative aspect of the method decreases with an increase of the data range, which appears as an averaging effect, consistent with the observations across the interface in the half-space case. Panels in the right two columns in Figure 6 show the data from Figures 6(a) to (f) normalized by the input wave speed maps (Figures 6a, b). The overall neutral or white backgrounds illustrate the well resolved reference value. The fact that circular features are visible in Figures 6(i) to (l) demonstrates that the velocity estimates around the edges are biased. The width of the halo or $\Delta L$ scales again with the data distance,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The control experiment. Focal spot obtained images for (a) $r_{fit} = 0.5\lambda_0$ and (b) $r_{fit} = 1\lambda_0$ for a homogeneous medium with wave speed $V = 2$ km/s. The reference wavelength $\lambda_0$ is indicated by the black line and the $r_{fit}$ range used for the nonlinear regression is the radius of the black circle. The indicated wave speeds $V$ are the domain average values and the standard deviation quantify the fluctuations across the domain. Panels (c) and (d) illustrate the spatial distribution of wave speed errors obtained with Equation (5) for $r_{fit} = 0.5\lambda_0$ and $r_{fit} = 1\lambda_0$, respectively.</figcaption></figure>

the reconstruction benefits from smaller $r_{fit}$ values (Figures 6i, j).

These interpretations are supported by wave speed profiles across the inclusions (Figure 7), which demonstrate again that the focal spot image quality, i.e., resolution and contrast, depends on the data range. Inclusions separated by $0.25\lambda_0$ (Figures 7a, b) are well imaged for small $r_{fit} = 0.5\lambda_0$. The wave speed value away from an interface is well approximated. This applies to the stiff and the compliant inclusions. For $r_{fit} = 1\lambda_0$ the contrast cannot be recovered, which is linked to the data range dependent transition zone width that is here similar to the inclusion diameter. The same mechanism applies to the results obtained with $r_{fit} = 1.5\lambda_0$. Importantly, two inclusions can always be discriminated when the separation distance is $0.25\lambda_0$, which is yet more obvious from Figure 6.

For a separation distance of $0.5\lambda_0$ this quality of the imaging approach increases, and inclusions are completely discriminated for a distance of $1\lambda_0$. Taken together, the data range dependent sensitivity of the resolution mostly affects the contrast estimate, but less so the ability to discriminate two objects. The power to discriminate inclusions separated by only $0.25\lambda_0$ suggests that the property of super-resolution applies, similar to focal spot medical imaging results obtained with ultrasound wavelengths [Zemzemi et al., 2020]. Focal spot imaging resolution thus benefits from high station density across short data ranges. The wave speed estimates are, however, sensitive to noisy data at short distances. In turn, longer data distances reduce fluctuations, which leads to a trade-off between accuracy or contrast and resolution or discrimination [Giammarinaro et al., 2023].
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Results for the two half-spaces. (a) The input model. (b) Reference (blue) and focal spot based velocity profiles obtained with two data ranges (red, black) that are averaged along the *N*-axis. Panels (c) and (d) show focal spot based images of the velocity. The reference wavelength $\lambda_0$ is indicated by the black line and the $r_{\text{fit}}$ range used for the nonlinear regression is the diameter of the black circle.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The dependence of the transition width $\Delta L$ on the data range $r_{\text{fit}}$. (a) Averaged input and imaged wave speed profiles using four different data ranges. Averaging is performed along the *N*-axis. The vertical dashed lines indicate the location where the amplitudes of the empirical profiles equal the 5% and 95% values of the 0.2 km/s velocity jump. The dotted lines indicate the estimates of the transition width $\Delta L$. (b) The transition width $\Delta L$ as a function of the data distance. (c) The averaged error estimate $\epsilon_V$. Colors in all panels correspond to the same $r_{\text{fit}}$ values.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Results of the circular inclusion resolution test. The left two columns panels (a–f) show absolute values, the right two columns panels (g–l) show scaled values. (a,b) Input models of the three circular inclusion pairs separated by $1\lambda_0$, $0.5\lambda_0$, and $0.25\lambda_0$. The inclusion diameter is $1\lambda_0$. Focal spot images obtained with (c,d) $r_{\text{fit}} = 0.5\lambda_0$ and (e,f) $r_{\text{fit}} = 1\lambda_0$. (g,h) Normalized input wave speed map. (i,j) Images in panels (c,d) scaled by the corresponding input distributions in panels (a,b). (k,l) Images in panels (e,f) scaled by the corresponding input distributions in panels (a,b). In the lower two rows, the reference wavelength $\lambda_0$ is indicated by the black line and the $r_{\text{fit}}$ range used for the nonlinear regression is the diameter of the black circle.</figcaption></figure>

### 3.4. Imaging random media

The last experiment explores focal spot imaging results of heterogeneous media that are characterized by a von Karman spectral density probability function (Equations (2), (3)). Figure 8 shows input wave speed distributions and the images obtained with data ranges $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$ together with the input-normalized images and the formal uncertainties obtained with Equation (5). For each case, we calculate a coefficient of correlation *R* between the input distribution and the image. As in the previous experimental configurations, the overall pattern of

the velocity variations is well retrieved. Peak *R* correlation values are obtained at zero-lag, which implies no systematic phase change and hence the effectiveness of focal spot “inference-mode” imaging [Tsai, 2023]. The smallest coefficient of correlation is *R* = 0.66 for the small scale, high contrast wave speed variation Medium 1 (*a* = $0.5\lambda_0$, *κ* = 0.1) and the longer data range $r_{\text{fit}} = 1\lambda_0$. The best estimation is obtained for Medium 9 (*a* = $2\lambda_0$, *κ* = 0.6) imaged with $r_{\text{fit}} = 0.5\lambda_0$, which corresponds to analysis using shorter ranges than the correlation length *a*. This leads to a high coefficient of correlation *R* = 0.97. Increasing $r_{\text{fit}}$ to $1\lambda_0$ yields *R* = 0.95. This again
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Cross-sections through the circular inclusions. Profiles of input (blue) and estimated wave speeds for $r_{fit} = 0.5\lambda_0$ (orange), $r_{fit} = 1\lambda_0$ (green), and $r_{fit} = 1.5\lambda_0$ (red) for each pair of the circular inclusions. Results are obtained along the E-axis passing through the center of the inclusions. (a,b) Inclusions separated by $1\lambda_0$. (c,d) Inclusions separated by $0.5\lambda_0$. (e,f) Inclusions separated by $0.25\lambda_0$.</figcaption></figure>

confirms our previous results, the images become smoother with increasing $r_{fit}$, which is synonymous with a loss in contrast and hence details. As an example, in row 6, the green 45-degree trending feature around position $x = 3\lambda$, $y = 2.5\lambda$ illustrates a narrow low-velocity zone that is overestimated by the averaging, longer $r_{fit} = 1\lambda_0$. The other way around, the smoother the input distribution for larger $a$ and $\kappa$ (Medium 3, 6, 9), the better is the estimate. The corresponding normalized results in rows 3 and 6 illustrate the imaging quality in a complementary style. Predominantly neutral distributions indicate overall good reconstructions, the red-green pattern amplitude range indicates the loss in accuracy, and the color-speckle size is related to $r_{fit}$. As an example, in row 6, the green 45-degree trending feature around position $x = 3\lambda$, $y = 2.5\lambda$ illustrates a narrow

low-velocity zone that is overestimated by the averaging, longer $r_{fit} = 1\lambda_0$. The spatial variations of the formal uncertainty estimates in row 4 indicate a medium and wave field dependence. As for the control experiment (Figure 3c) the generally larger error in the lower region is governed by the configuration (Figure 1). The resolution of the comparatively large error spot at $x \approx 6\lambda$, $y \approx 3\lambda$ in the short correlation length media 1 to 6 highlight the advantage of the local error estimation. Again as for the control experiment, the larger data range in row 7 reduces the uncertainty associated with fluctuations significantly.

We estimate the relative wave speed change $\xi$ from the wave speed images by inverting Equation (1)

$$ \xi(\mathbf{x}) = V/V_0 - 1, \qquad (6) $$
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Results from the random media case. The top row shows the heterogeneous velocity distributions that are synthesized using Equations (2) and (3). Values of the corresponding tuning parameters are given in Table 1. Rows 2–4 and rows 5–7 show focal spot based results obtained with $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$, respectively. For each case, the top rows 2 and 5 show the estimated wave speed distributions, the center rows 3 and 6 show the ratio between the estimated and the input wave speed, and the bottom rows 4 and 7 show the local wave speed error estimates. For each panel, the vertical axis and the horizontal axis correspond respectively to the North and East axis, scaled by the mean input wavelength. The *R* coefficient is the correlation between the input and the estimated map.</figcaption></figure>

with $V_0 = 2$ km/s. We compile histograms of the $\xi$ distributions to compare properties of the reference input distributions and of the obtained images. Figure 9 collects the $\xi$ histograms from the input and from the estimates for $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$. The similarity between reference and image is better for $r_{\text{fit}} = 0.5\lambda_0$ than $r_{\text{fit}} = 1\lambda_0$. The best result is obtained for the smoothest Medium 9. The green histograms obtained with $r_{\text{fit}} = 1\lambda_0$ are more narrow and have a higher peak value around small $\xi$ values compared to the orange $r_{\text{fit}} = 0.5\lambda_0$ results. This indicates again

the low-pass filter property of larger data ranges observed in the previous experiments.

## 4. Discussion

Resolution can mean different things in different imaging contexts, including the number or density of measuring points, the spectral sensitivity of an imaging device or method, the ability to detect or discriminate features and to accurately estimate their properties, contrast in brightness or color, or phase
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Histograms of the relative wave speed change $\xi$ for the nine media imaged in Figure 8. Blue data correspond to the input reference values, and orange and green data correspond to images obtained with $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$, respectively.</figcaption></figure>

fidelity [e.g., Smith, 2003]. We use numerical simulations of two-dimensional acoustic wave propagation in a cavity (Figures 1, 2) to investigate the lateral resolution of the focal spot imaging technique for a fixed acquisition system with a constant number of grid points. The increase in depth resolution for such a compact dense array configuration compared to measurements made on traveling seismic Rayleigh surface waves is established by Giammarinaro et al. [2023]. We implement four test cases that together allow us to study the lateral resolution power of Rayleigh wave focal spots reconstructed from vertical–vertical component noise correlation data. Lateral resolution is discussed as the ability to resolve a step function in the material properties, to discriminate and characterize two closely spaced objects, and to measure position, amplitude, and phase of random distributions.

In a first homogeneous control experiment (Figures 2, 3) the reference wave speed is estimated with an error below 2% on average, which includes, however, areas of larger bias associated with edge effects. This average error is larger compared to the focal spot results based on numerical time-reversal experiments, where noise-free synthetics lead to errors in the 0.01% range for vertical–vertical component records and data ranges of $r_{\text{fit}} = 0.5\lambda_0$ and $r_{\text{fit}} = 1\lambda_0$, and where anisotropic surface wave incidence results in biases in the 0.1% range [Giammarinaro et al., 2023]. The different error levels are associated with the different methods used to synthesize Green's functions and focal spots. In time-reversal experiments, the wave field and hence the ballistic wave correlations are fully controlled by the mirror properties. More mirror elements lead to better refocusing results. Controlled lab experiments

can stack over different space realizations. In seismic data applications an improved Green's function and refocusing reconstruction is achieved by time averaging to better conform with the decorrelated noise source assumption. Here, the focal spots are retrieved from cross-correlation of the reverberating cavity wave field, where the quality of the Green's function is controlled by the ability to excite and average a sufficiently large number of modes in the cavity [Draeger and Fink, 1999]. Hence we stack over different realizations using different source positions. This increases the number of independent modes to enhance the narrow-band refocusing, which is equivalent to using more time reversal mirror elements in a time-reversal experiment. Our approach converges towards the theoretical focal spot, but it remains sensitive to details of the implementation such as the source positions of the relatively few realizations. This explains the observed fluctuations in the focal spot reconstructions, which indicate the imperfect Green's function synthesis, even for the homogeneous case. This means that our imperfect cavity results are comparable to focal spots obtained from noisy data.

The second half-spaces experiment (Figures 4, 5) shows the feasibility to resolve the contrast between two media that have a 10% difference in wave speed. However, the interface is not perfectly resolved. Whereas the wave speed is correctly estimated away from the interface, the finite data or fitting range creates an averaging or low-pass filter effect that depends on $r_{\text{fit}}$. Tests with different $r_{\text{fit}}$ values indicate that the transition width $\Delta L$ scales with good approximation linearly with the data range, $\Delta L \approx r_{\text{fit}}$.

The third circular inclusion configuration (Figures 6, 7) shows the possibility to discriminate two
separate objects separated by $0.25\lambda_0$, even for data ranges of one wavelength. However, the data range dependent low-pass filter properties lead to averaged amplitude values at edges. We have not observed situations where the focal spot method yields biased phase properties, so the inclusion positions are accurate. This suggests that the super-resolution property demonstrated with passive elastography in soft tissues [Zemzemi et al., 2020] also applies to 2D Rayleigh surface wave propagation. Again this means that for good data quality and high station density the method has the potential to meet the formal criterion of super-resolution, i.e., the sensitivity at small scales is sufficient to discriminate objects or features that are separated by distances that are much shorter than the wavelength.

The fourth test case considers random media (Figures 8, 9, Table 1) which show the highest similarity to distributions of Earth material properties. The quality of the focal spot reconstruction as quantified by the correlation coefficient $R$ between reference and image depends on the roughness or smoothness of the distribution in relation to wavelength and data range. $R$ is small when the distributions are rough, have small-scale fluctuations compared to the probing wavelength, and the data range is large. The reconstruction is almost perfect when the distributions are smooth, have large-scale fluctuations, and the data range is small. These conclusions are further corroborated comparing histogram properties of the velocity variation parameter $\xi$ (Equation (2), Figure 9), which again demonstrate the low-pass filtering effect of large $r_{fit}$ values. Thus, positions are well estimated but the amplitudes diverge for small-scale heterogeneities. The best estimates are obtained for variations on scales larger than $r_{fit}$. Deblurring methods can potentially be developed to further improve the lateral resolution of the method considering that the focal spot Bessel function shape acts as a spatial convolution filter. If resources permit, an imaging campaign could be optimized by first detecting target features with low contrast before then a sensor re-configuration helps to improve the quantitative estimate by increasing the signal-to-noise ratio through network densification.

The present study employs two-dimensional acoustic simulations. These scalar simulations yield the same Green's function as for vertical-vertical component Rayleigh wave propagation. Our main

conclusions thus hold for seismic Rayleigh wave imaging, including the advantageous estimation of the local error that can help improve the uncertainty management of shear wave inversions. However, this set-up does not allow to study the biasing effect of interfering body wave energy. Giammarinaro et al. [2023] shows that the error on the vertical-vertical component increases in the presence of P waves, but this can be compensated for by increasing the data range. Together with our observations here this implies that an increase in $r_{fit}$ improves the results if refocusing P wave energy distorts the surface wave focal spot, but that this remedy negatively affects lateral resolution power. This trade-off situation would benefit from efficient focal spot filters. Alternatively, Rayleigh wave phase speed estimates obtained from radial-vertical component data are much less sensitive to P waves [Giammarinaro et al., 2023], which offers independent constraints for the improvement of vertical-vertical results. The application of radial-radial and transversal-transversal component focal spots requires first an efficient separation of the combined Rayleigh and Love wave energy [Haney and Nakahara, 2014]. Spatial autocorrelations of surface wave strain and rotational data can further provide additional constraints on the local velocity structure [Nakahara et al., 2021, Nakahara and Haney, 2022].

## 5. Conclusion

We investigate the lateral resolution power of the Rayleigh surface wave focal spot imaging method using two-dimensional numerical experiments of reverberating wave fields in a cavity. Most importantly the resolution depends on the data range. This means that focal spot imaging exhibits super-resolution properties provided the data quality supports sub-wavelength data ranges. Longer data ranges still allow imaging of small-scale features at super-resolution albeit with a loss in contrast. Seismic Rayleigh wave focal spot imaging shows convincing resolution properties that make it suitable for a wide range of imaging applications ranging from feature detection to accurate wave speed estimates. There are hence no fundamental disadvantages compared to established passive surface wave tomography methods. Here as there, the station configuration can be tuned to support image quality and properties
for different goals, and in both cases the data quality or signal-to-noise ratio ultimately has the largest impact on the resolution.

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## Open research/data availability

No observational data were used in this study.

## Dedication

The manuscript was written with contributions from all authors. All authors have given approval to the final version of the manuscript.

## Funding

This work was sponsored by a research grant from the Academy of Finland (decision number 322421).

## Acknowledgments

We thank Julien de Rosny, Stefan Catheline, Michel Campillo, Léonard Seydoux, Laurent Stehly, Alexander Meaney, and Markus Juvonen for helpful discussions. We thank the editor N. Shapiro and two anonymous reviewers for comments that helped to improve the manuscript.

## References

* Aki, K. (1957). Space and time spectra of stationary stochastic waves, with special reference to microtremors. *Bull. Earthq. Res. Inst.*, 35, 415-456.

* Asten, M. W. (2006). On bias and noise in passive seismic data from finite circular array data processed using SPAC methods. *Geophysics*, 71(6), V153-V162.

* Aster, R. C., Borchers, B., and Thurber, C. H. (2019). *Parameter Estimation and Inverse Problems*. Elsevier, Amsterdam, Oxford, Cambridge, 3rd edition.

* Barrere, V., Melodelima, D., Catheline, S., and Giammarinaro, B. (2020). Imaging of thermal effects during high-intensity ultrasound treatment in liver by passive elastography: a preliminary feasibility in vitro study. *Ultrasound Med. Biol.*, 46, 1968-1977.

* Ben-Zion, Y. (1989). The response of two joined quarter spaces to SH line sources located at the material discontinuity interface. *Geophys. J. Int.*, 98(2), 213-222.

* Benech, N., Brum, J., Catheline, S., Gallot, T., and Negreira, C. (2013). Near-field effects in Green's function retrieval from cross-correlation of elastic fields: experimental study with application to elastography. *J. Acoust. Soc. Am.*, 133, 2755-2766.

* Brum, J., Catheline, S., Benech, N., and Negreira, C. (2015). Quantitative shear elasticity imaging from a complex elastic wavefield in soft solids with application to passive elastography. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control*, 62(4), 673-685.

* Catheline, S., Benech, N., Brum, J., and Negreira, C. (2008). Time reversal of elastic waves in soft solids. *Phys. Rev. Lett.*, 100(6), article no. 064301.

* Catheline, S., Souchon, R., Rupin, M., Brum, J., Dinh, A. H., and Chapelon, J.-Y. (2013). Tomography from diffuse waves: Passive shear wave imaging using low frame rate scanners. *Appl. Phys. Lett.*, 103(1), article no. 014101.

* Chamarczuk, M., Malinowski, M., Nishitsuji, Y., Thorbecke, J., Koivisto, E., Heinonen, S., Juurela, S., Mężyk, M., and Draganov, D. (2019). Automatic 3D illumination-diagnosis method for large-N arrays: Robust data scanner and machine-learning feature provider. *Geophysics*, 84(3), Q13-Q25.

* Cotton, F. and Coutant, O. (1997). Dynamic stress variations due to shear faults in a plane-layered medium. *Geophys. J. Int.*, 128, 676-688.

* Cox, H. (1973). Spatial correlation in arbitrary noise fields with application to ambient sea noise. *J. Acoust. Soc. Am.*, 54, 1289-1301.

* Derode, A., Larose, E., Tanter, M., de Rosny, J., Tourin, A., Campillo, M., and Fink, M. (2003). Recovering the Green's function from field-field correlations in an open scattering medium (L). *J. Acoust. Soc. Am.*, 113(6), 2973-2976.

* Draeger, C. and Fink, M. (1999). One-channel time-reversal in chaotic cavities: Theoretical limits. *J. Acoust. Soc. Am.*, 105(2), 611-617.

* Ekström, G. (2014). Love and Rayleigh phase-velocity maps, 5-40 s, of the western and central USA from
USArray data. *Earth Planet. Sci. Lett.*, 402(C), 42–49.

Ekström, G., Abers, G. A., and Webb, S. C. (2009). Determination of surface-wave phase velocities across USArray from noise and Aki's spectral formu-lation. *Geophys. Res. Lett.*, 36(18), 5–9.

Fink, M. (1997). Time reversed acoustics. *Phys. Today*, 50(3), 34–40.

Frankel, A. and Clayton, R. W. (1986). Finite differ-ence simulations of seismic scattering: Implica-tions for the propagation of short-period seismic waves in the crust and models of crustal hetero-geneity. *J. Geophy. Res.*, 91(B6), 6465–6489.

Gallot, T., Catheline, S., Roux, P., Brum, J., Benech, N., and Negreira, C. (2011). Passive elastography: shear-wave tomography from physiological-noise correlation in soft tissues. *IEEE Trans. Ultrason. Ferroelectr. Freq. Control*, 58(6), 1122–1126.

Giammarinaro, B. and Hillers, G. (2022). Resources for “Seismic surface wave focal spot imaging: numerical resolution experiments. University of Helsinki.

Giammarinaro, B., Tsarsitalidou, C., Hillers, G., de Rosny, J., Seydoux, L., Catheline, S., Campillo, M., and Roux, P. (2023). Seismic surface wave focal spot imaging: numerical resolution experiments. *Geophys. J. Int.*, 232(1), 201–222.

Haney, M. M., Mikesell, T. D., van Wijk, K., and Nakahara, H. (2012). Extension of the spatial autocorrelation (SPAC) method to mixed-component correlations of surface waves. *Geophys. J. Int.*, 191(1), 189–206.

Haney, M. M. and Nakahara, H. (2014). Surface-wave Green's tensors in the near field. *Bull. Seism. Soc. Am.*, 104(3), 1578–1586.

Hillers, G., Mai, P. M., Ben-Zion, Y., and Ampuero, J.-P. (2007). Statistical properties of seismicity of fault zones at different evolutionary stages. *Geophys. J. Int.*, 169, 515–533.

Hillers, G., Roux, P., Campillo, M., and Ben-Zion, Y. (2016). Focal spot imaging based on zero lag cross-correlation amplitude fields: Application to dense array data at the San Jacinto fault zone. *J. Geophy. Res. Solid Earth*, 121(11), 8048–8067.

Holliger, K. and Levander, A. R. (1992). A stochastic view of lower crustal fabric based on evidence from the Ivrea Zone. *Geophys. Res. Lett.*, 19(11), 1153–1156.

Lin, F.-C. and Ritzwoller, M. H. (2011). Helmholtz

surface wave tomography for isotropic and az-imuthally anisotropic structure. *Geophys. J. Int.*, 186(3), 1104–1120.

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

Liu, Q. and Gu, Y. (2012). Seismic imaging: From classical to adjoint tomography. *Tectonophysics*, 566–567, 31–66.

Mai, P. M. and Beroza, G. C. (2002). A spatial random field model to characterize complexity in earth-quake slip. *J. Geophy. Res. Solid Earth*, 107(B11), ESE 10–1–ESE 10–21.

Nakahara, H. (2006). A systematic study of theoretical relations between spatial correlation and Green's function in one-, two- and three-dimensional ran-dom scalar wavefields. *Geophys. J. Int.*, 167, 1097–1105.

Nakahara, H., Emoto, K., and Nishimura, T. (2021). Extending the formulation of the spatial autocorrelation (SPAC) method to strain, rotation and tilt. *Geophys. J. Int.*, 227(1), 287–302.

Nakahara, H. and Haney, M. M. (2022). Connection between the cross correlation and the Green's func-tion: strain and rotation of surface waves. *Geo-phys. J. Int.*, 230(2), 1166–1180.

Obermann, A., Planès, T., Hadziioannou, C., and Campillo, M. (2016). Lapse-time-dependent coda-wave depth sensitivity to local velocity perturba-tions in 3-D heterogeneous elastic media. *Geo-phys. J. Int.*, 207(1), 59–66.

Ripperger, J., Ampuero, J.-P., Mai, P. M., and Giardini, D. (2007). Earthquake source characteristics from dynamic rupture with constrained stochastic fault stress. *J. Geophys. Res.*, 112(B4), article no. B04311.

Sabra, K. G., Conti, S., Roux, P., and Kuperman, W. A. (2007). Passive in vivo elastography from skele-tal muscle noise. *Appl. Phys. Lett.*, 90(19), article no. 194101.

Sabra, K. G., Gerstoft, P., Roux, P., Kuperman, W. A., and Fehler, M. C. (2005). Surface wave tomography from microseisms in Southern California. *Geophys. Res. Lett.*, 32(14), 1–4.

Sanchez-Sesma, F. J. (2006). Retrieval of the Green's function from cross correlation: the canonical elas-tic problem. *Bull. Seism. Soc. Am.*, 96(3), 1182–1191.

Sato, H., Fehler, M. C., and Maeda, T. (2012). *Seis-
*mic Wave Propagation and Scattering in the Heterogeneous Earth: Second Edition.* SpringerLink.
Springer, Berlin, Heidelberg, 2nd edition.

Shapiro, N. M. and Campillo, M. (2004). Emergence of broadband Rayleigh waves from correlations of the seismic ambient noise. *Geophys. Res. Lett.*, 31, article no. L07614.

Smith, S. W. (2003). *Digital Signal Processing: A Practical Guide for Engineers and Scientists*. Demystifying Technology Series. Newnes, Amsterdam, Boston.

Treeby, B. E., Budisky, J., Wise, E. S., Jaros, J., and Cox, B. T. (2018). Rapid calculation of acoustic fields from arbitrary continuous-wave sources. *J. Acoust. Soc. Am.*, 143(1), 529–537.

Tromp, J., Tape, C., and Liu, Q. (2005). Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels. *Geophys. J. Int.*, 160(1), 195–216.

Tsai, V. C. (2023). The future of Earth imaging. *Seismol. Res. Lett.*, 94(5), 2119–2128.

Tsai, V. C. and Moschetti, M. P. (2010). An explicit relationship between time-domain noise correlation and spatial autocorrelation (SPAC) results. *Geo-phys. J. Int.*, 182(1), 454–460.

Tsarsitalidou, C., Boué, P., Hillers, G., Giammarinaro, B., Campillo, M., Seydoux, L., and Stehly, L. (2021). Seismic imaging with focusing surface waves obtained from USArray noise correlation functions. In *EGU General Assembly Conference Abstracts*, pages EGU21–9045. EGU.

Weertman, J. (1980). Unstable slippage across a fault that separates elastic media of different elastic constants. *J. Geophy. Res. Solid Earth*, 85(B3), 1455–1461.

Yokoi, T. and Margaryan, S. (2008). Consistency of the spatial autocorrelation method with seismic interferometry and its consequence. *Geophys. Prospect.*, 56, 435–451.

Zemzemi, C., Zorgani, A., Daunizeau, L., Belabhar, S., Souchon, R., and Catheline, S. (2020). Super-resolution limit of shear-wave elastography. *Euro-phys. Lett.*, 129(3), article no. 34002.

Zorgani, A., Souchon, R., Dinh, A.-H., Chapelon, J.-Y., Ménager, J.-M., Lounis, S., Rouvière, O., and Catheline, S. (2015). Brain palpation from physiological vibrations using MRI. *Proc. Natl. Acad. Sci. USA*, 112, 12917–12921.