On the variability of sources in ambient seismic noise source inversion

SUMMARY

We present a method for ambient noise cross-correlation modelling and source inversion, which accounts for spatiospectral variability in noise source distributions. It is based on numerical wavefield simulations in 2-D acoustic media. The source power spectral density is parametrized by a sum of a small number of spatial source distributions, each with a corresponding frequency spectrum held fixed during the inversion. Algorithmically, this is an extension of our previous work which assumed spatially homogeneous source spectra. In this paper, we use it to study the impact of incorrectly estimating source spectra from observed data. This is done using synthetic tests involving sources with closely spaced frequency spectra. The tests demonstrate that when the spatial variability of sources is either partially or wholly unaccounted for, the recovery of true source locations is compromised.

inverse theory, seismic interferometry, Seismic noise, waveform inversion

1 INTRODUCTION

Ambient noise interferometry seeks to extract a deterministic signal (such as a Green function) from quasi-random recordings of seismic ambient noise. The seismological literature abounds with theoretical studies establishing the equivalence between the cross-correlation of ambient noise recorded at two receivers, and the Green function between them, provided the noise field is diffuse and equipartitioned (Lobkis & Weaver 2001; Derode et al. 2003; Snieder 2004; Wapenaar 2004; Weaver 2005). On the other hand, Tromp et al. (2010) showed that the ‘ensemble-averaged’ cross-correlation (its expectation value over many realizations of the random noise sources) between two seismic sensors can be modelled in terms of the power spectral density of arbitrarily distributed noise sources. In other words, a deterministic signal can be obtained from quasi-random noise recordings, even when the noise field is bereft of special properties like diffusivity and equipartitioning. The only requirement is that the noise sources be spatially uncorrelated over the timescale of ensemble averaging. As with any geophysical observation, if it can be modelled, it can be inverted for the governing model parameters. Therefore, the work of Tromp et al. (2010) paved the way for a new line of research in ambient noise interferometry: imaging of noise sources (Datta et al. 2019, 2023; Xu et al. 2019, 2020; Ermert et al. 2021; Igel et al. 2021; Xu & Mikesell 2022; Igel et al. 2023) as well as Earth structure (Tromp et al. 2010; Fichtner et al. 2017; Sager et al. 2018, 2020), by cross-correlation modelling and inversion, that is, interferometry without Green’s function retrieval (see Fichtner & Tsai 2019). Recently, the literature highlighting the importance of this rigorous approach to ambient noise seismology, has evolved from pointing out errors in estimated Green’s functions under realistic scenarios of ambient noise sources (Halliday & Curtis 2008; Tsai 2009; Yao & van der Hilst 2009; Froment et al. 2010; Kimman & Trampert 2010) to explicitly quantifying errors in tomographic models built using the assumption of Green’s function retrieval (Valero Cano et al. 2024).

Within the realm of ambient noise cross-correlation modelling and inversion, this paper is concerned exclusively with ‘ambient noise source inversion’. Our use of this term refers specifically to the estimation of noise source power spectral density (PSD) via tomographic inversion of cross-correlation data, under the assumption of a fixed Earth structure model. This may not be a direct priority for many ambient noise seismologists, but it is an integral part of full waveform ambient noise tomography because, as with earthquake tomography, one needs a noise source model to invert for the Earth structure (and vice versa).

In general, the source PSD is space- as well as frequency-dependent, so, sensu stricto, ambient noise source inversion must ascertain the source spectra at every point in space. To keep the inverse problem tractable however, all practitioners of the technique make simplifying assumptions. The zeroth order simplification is the so-called fixed-spectrum inversion—noise sources are assumed to be spatially homogeneous with a known frequency spectrum, and only the space-dependent ‘source strength’ distribution is inverted for (Ermert et al. 2017; Datta et al. 2019, 2023; Xu et al. 2019). There are two limitations of this approximation. First, one needs a way of estimating the source spectrum. Secondly, it is not applicable in scenarios where multiple sources, having different frequency spectra (e.g. motorway + railway line + factory, for local scale studies), may be simultaneously active. If the different sources are well separated in frequency, one can filter the data in different frequency bands and hope to invert separately for the source strength in different bands. If not, the inversion could be contaminated by frequency-smearing, although these effects have not been explicitly studied. To overcome these limitations, Ermert et al. (2020) allowed for spatially varying source PSDs, and Ermert et al. (2021) developed a technique for joint spatiospectral inversion to recover the space- as well as frequency-dependence of noise sources. To our knowledge, there are no other studies that have developed multifrequency ambient seismic source inversion schemes. In fact, the fixed-spectrum assumption has been used as recently as in Valero Cano et al. (2024).

Our objective in this study is two-fold. First, we demonstrate a method for noise cross-correlation modelling and inversion, with a source PSD that incorporates general spatiospectral variability (Section 2). It is an extension of the work by Datta et al. (2023) – hereafter DSK23 – so the various features and benefits of DSK23 are retained. These include relatively high tolerance for velocity model inaccuracies, improved resolution of noise sources outside the receiver array and a fast-converging, Hessian-based optimization scheme. The chief limitation of DSK23, that it is a local-scale technique which does not account for Earth’s sphericity or distant noise sources (outside the computational domain), is also inherited. Valero Cano et al. (2024) reported that the inability to account for distant sources, which is an inherent feature of local or regional scale ambient noise studies, produces errors in structure models that are comparable to those caused by (implicitly) assuming a homogeneous noise source distribution. However, no non-global study can truly account for all possible ambient seismic sources, and local-scale analyses are pertinent when the study area is situated sufficiently far away (relative to frequencies considered) from major distant sources such as ocean microseisms. For example in cryoseismology, one may be interested in studying sources of ambient seismicity within a glacier, such as water tremor from englacial or subglacial drainage channels (e.g. Labedz et al. 2022).

Our second objective in this paper, with a view towards practical applications, is to examine the effects of inaccuracies in assumed source spectra, on ambient noise source inversion. Noise source spectra are often estimated from observed cross-correlation data (e.g. Datta et al. 2019; Xu et al. 2020). This can lead to inaccurate estimates which do not (fully) account for the spatiospectral variability of sources, because one may fail to recognize different source signatures. Section 3 presents synthetic tests that systematically consider such inaccurate source spectra.

2 MODELLING METHODOLOGY

Ambient noise source inversion, as defined in Section 1, is founded on the following model for the ensemble cross-correlation (Tromp et al. 2010; Fichtner & Tsai 2019):

$$\begin{eqnarray} C(\mathbf {x}_A, \mathbf {x}_B;\omega ) = \int G(\mathbf {x}_A, \mathbf {x}; \omega ) G^{*}(\mathbf {x}_B, \mathbf {x}; \omega ) S(\omega , \mathbf {x})\, {\rm d}^{2}\mathbf {x}, \end{eqnarray}$$

(1)

where |$\mathrm{\mathbf {x}_A}$| and |$\mathrm{\mathbf {x}_B}$| are two receiver locations, G is the medium Green’s function and |$S(\omega , \mathbf {x})$| is the source PSD. As in DSK23, we limit ourselves to the modelling of scalar (e.g. acoustic) wavefields, but the ideas presented in this paper also directly apply to the elastic scenario, where C, G and S would each be a second-order tensor. The variable-separable form of |$S(\omega , \mathbf {x})$| used in DSK23, is now generalized to:

$$\begin{eqnarray} S(\omega , \mathbf {x}) = \sum _{l=1}^{N_s} {P_l(\omega ) \sigma _l(\mathbf {x})}. \end{eqnarray}$$

(2)

This means that the source PSD is a superposition of different strength distributions |$\sigma _l(\mathbf {x})$|⁠, each weighted by its own power spectrum |$P_l(\omega )$| (see Fig. 1). The number of contributions in this sum, |$N_s$|⁠, is expected to be small, because it represents different types of sources active in different frequency bands. The different |$P_l(\omega )$|⁠, for different values of l, are independent of each other inasmuch as their shapes or functional forms can be different. This is a slight departure from the work of Ermert et al. (2021), who used a well-defined spectral basis of orthogonal sinc functions. We do not present a real data application in this paper, but envision that in actual applications, the choice of the functions |$P_l(\omega )$|⁠, will vary on a case-to-case basis.

Figure 1.

An example of a noise source PSD model with |$N_s=2$|⁠. (a), (c) the two source strength distributions, in a |$50 \times 50$| km modelling region containing 20 receivers as shown by the small triangles. The receiver numbers marked in (c) are referenced in Fig. 2. (b), (d) the corresponding frequency spectra, which are both Gaussian functions peaked at 0.2 and 0.5 Hz, respectively, with standard deviations of 0.05 and 0.08 Hz, respectively.

We note that |$N_s=1$| in eq. (2) corresponds to DSK23, whereas the case |$N_s>1$| necessitates algorithmic changes, which are described in this section. We highlight the changes relative to DSK23, and point the reader to that and other past papers for details that remain unchanged.

2.1 Forward modelling

We simulate acoustic wave propagation using a time-domain (finite difference) solver (Devito, Louboutin et al. 2019). While it is theoretically possible to obtain band-limited Green’s functions using filtered delta functions in finite difference simulations, we use a band-limited source function to ensure numerical stability and accuracy. With |$N_s=1$| in eq. (2), it is easy to evaluate eq. (1) using source–receiver reciprocity and point sources only, as done by DSK23. However for |$N_s>1$|⁠, this point-source approach would require |$N_s$| simulations per receiver, one for each source, |$\sqrt{P_l(\omega )}$|⁠. It is more efficient to employ the distributed-source algorithm of Tromp et al. (2010), which is valid for any general |$S(\omega , \mathbf {x})$|⁠. This algorithm has been discussed by other authors (e.g. Sager et al. 2018), but our adaptation involves an extra deconvolution step as we do not have direct access to Green’s functions. For the sake of completeness, we list here all the steps in the recipe for evaluating (1). It is illuminating to rewrite eq. (1) as:

$$\begin{eqnarray} C(\mathbf {x}, \mathbf {x}_B) = \int G(\mathbf {x}, \xi ) \left[G(\xi , \mathbf {x}_B)^{*}\, S(\xi )\right]\, {\rm d}^{2}\xi , \end{eqnarray}$$

(3)

where source–receiver reciprocity has been invoked at the receiver location |$\mathbf {x}_B$|⁠. The |$\omega$|-dependence of each term has been omitted for brevity. With reference to eq. (3), the recipe for evaluating |$C(\mathbf {x}_A, \mathbf {x}_B)$|⁠, is:

• Inject a dummy source |$w(\omega )$| at |$x_B$|⁠, to evaluate the wavefield |$w(\omega )G(\xi ,\mathbf {x}_B)$|

• Multiply the complex conjugate of the wavefield computed above, with |$S(\xi )$|

• Use |$w^{*}(\omega )G^{*}(\xi ,\mathbf {x}_B)S(\xi )$| as a distributed source to solve for the wavefield |$C^{\prime }(\mathbf {x}, \mathbf {x}_B)=w^{*}(\omega )C(\mathbf {x}, \mathbf {x}_B)$|

• Sample this wavefield at |$\mathbf {x}=\mathbf {x}_A$| and deconvolve |$w^{*}(\omega )$|⁠, to obtain the cross-correlation |$C(\mathbf {x}_A, \mathbf {x}_B)$|⁠.

Choosing |$w(t)$| such that its bandwidth spans that of the correlograms to be modelled, reduces the artefacts invariably associated with deconvolution. This is verified in Fig. 2, where we compare the results of the above numerical recipe, with an analytical calculation in a homogeneous wave speed (⁠|$c=2$| km s⁻¹) model, assuming the source PSD of Fig. 1. The analytical calculation invokes source–receiver reciprocity at both receivers, and uses |$G(\mathbf {x}, \mathbf {x}_i) = H_0^{(1)}(\frac{\omega }{c}(|\mathbf {x}-\mathbf {x}_i|))$|⁠, |$H_0^{(1)}$| being the Hankel function of the first kind of order 0 (Hanasoge 2013; Datta et al. 2019). The minor discrepancies in the waveforms of Fig. 2 are likely due to the water-level regularization (Menke 2012) applied to stabilize the deconvolution.

Figure 2.

Comparison of cross-correlation waveforms computed analytically and numerically, as described in the main text, for the source PSD model of Fig. 1, and assuming a homogeneous wave speed of 2 km s⁻¹. The eight plots in this figure show waveforms for eight randomly chosen receiver pairs, as labelled in the top left corner of each plot.

2.2 Inversion

In the inverse problem of estimating |$S(\mathbf {x},\omega )$| from observed cross-correlation waveforms, we only invert for the spatial components |$\sigma _l(\mathbf {x})$|⁠. |$N_s$| and each of the |$P_l(\omega )$|⁠, are fixed a priori, so there is no spectral inversion. The spatial parametrization is the same as in DSK23, with a set of |$N_b$| 2-D Gaussian basis functions, so the total number of model parameters is limited to |$N_s \times N_b$|⁠. We have

$$\begin{eqnarray} \sigma _l(\mathbf {x}) = \sum _{j=1}^{N_b} (m^l_j)^2B_j(\mathbf {x}), \end{eqnarray}$$

(4)

where |$B_j(\mathbf {x})$| are the spatial basis functions, with coefficients |$m_j$|⁠. We note that the index l does not appear on the basis |$B_j(\mathbf {x})$|⁠, which must be the same for all frequency bands l. However, there is freedom in choosing the particular form of |$B_j(\mathbf {x})$|⁠, which can be exercised to implement different types of spatial parametrizations (e.g. a ring of sources as in Datta et al. 2019).

The discrete inverse problem of estimating the model coefficients |$m^l_j$|⁠, is solved iteratively using the Gauss–Newton optimization scheme of DSK23, but with the gradient and Hessian matrices augmented to incorporate the larger model space (for |$N_s > 1$|⁠). To see this, we go back to the source sensitivity kernels (Hanasoge 2013), which form the basis for inverse modelling.

Given any measurable quantity d, its variation is related to the variation in the (cross-correlation) waveform, |$\delta C$|⁠, via the ‘adjoint source’ f (so named due to its role in adjoint methods):

$$\begin{eqnarray} \delta d = \int f^{*} \delta C\, {\rm d}\omega . \end{eqnarray}$$

(5)

Using (1) and (2) in (5), under the assumption of a fixed structure model, we have

$$\begin{eqnarray} \delta d &= \iint f^{*}\, G^{*}(\mathbf {x}_A, \mathbf {x}; \omega )\, G(\mathbf {x}_B, \mathbf {x}; \omega )\, \delta S(\mathbf {x}, \omega )\, {\rm d}\omega \, {\rm d}^{2}\mathbf {x}\\ &= \sum _l \iint f^{*}\, G^{*}(\mathbf {x}_A, \mathbf {x}; \omega )\, G(\mathbf {x}_B, \mathbf {x}; \omega )\, P_l(\omega )\, {\rm d}\omega \, \delta \sigma _l(\mathbf {x})\, {\rm d}^{2}\mathbf {x}\\ &= \int \left[ \sum _l K^{AB}_l \delta \sigma _l \right]\, {\rm d}^{2}\mathbf {x} , \end{eqnarray}$$

where

$$\begin{eqnarray} K^{AB}_l(\mathbf {x}) = \int {f^{*}\, G^{*}(\mathbf {x}_A, \mathbf {x}; \omega )\, G(\mathbf {x}_B, \mathbf {x}; \omega ) P_l(\omega )} \, {\rm d}\omega , \end{eqnarray}$$

(6)

is the (spatial) source sensitivity kernel for the |$\mathbf {x}_A$|-|$\mathbf {x}_B$| receiver pair, with respect to the l^th source strength distribution, |$\sigma _l(\omega )$|⁠. Thus there are |$N_s$| sets of sensitivity kernels (⁠|$K^{AB}$| or |$K^i$|⁠) for each receiver pair and they combine to produce |$N_s$| misfit kernels (K) for the tomographic inverse problem:

$$\begin{eqnarray} K_l(\mathbf {x}) = \sum _i \Delta d_i K^i_l(\mathbf {x}) . \end{eqnarray}$$

(7)

Here the index i runs over all receiver pairs, and |$\Delta d$| is the data misfit. Given a misfit function |$\chi$| which quantifies the difference between observed and predicted values of d, the misfit gradient is given by |$\delta \chi =-\int \sum _l K_l(\mathbf {x}) \delta \sigma _l( \mathbf {x})\, \mathrm{d}^{2}\mathbf {x}$|⁠. The particular misfit function used in our study is

$$\begin{eqnarray} \chi =\dfrac{1}{2}\sum _{i}\left( \ln \dfrac{E_{i}^{\mathrm{obs}}}{E_{i}^{\mathrm{syn}}}\right)^{2}, \end{eqnarray}$$

(8)

where the index i runs over all receiver pairs and E is the energy of the cross-correlation waveform, contained in a time window |$w(t)$|⁠:

$$\begin{eqnarray} E_{i} = \sqrt{\int w(t) C_{i}^2(t)\, {\rm d}t}. \end{eqnarray}$$

(9)

For further details, the reader is referred to the foundation papers (DSK23; Hanasoge 2013; Datta et al. 2019).

Once the misfit kernels are obtained, the gradient of the misfit g, and the Jacobian matrix G (used to approximate the Hessian), are built by ‘stacking’ the |$N_s$| contributions row-wise (vertically) and column-wise (horizontally) respectively, that is,

$$\begin{eqnarray} g = {\begin{bmatrix}g^1 \\ g^2 \\ \vdots \end{bmatrix}} \end{eqnarray}$$

(10)

$$\begin{eqnarray} G = {\begin{bmatrix}G^1 & G^2 & \cdots \end{bmatrix}} , \end{eqnarray}$$

(11)

where

$$\begin{eqnarray} g^l_j=-\int 2 K_l(\mathbf {x}) m^l_j B_j(\mathbf {x})\, \mathrm{d}^2 \mathbf {x}, \end{eqnarray}$$

(12)

$$\begin{eqnarray} G^l_{ij} = \int 2 K^i_l(\mathbf {x}) m^l_j B_j(\mathbf {x}) \mathrm{d}^2 \mathbf {x}. \end{eqnarray}$$

(13)

3 SYNTHETIC TESTS FOR INVERSION

When working with real data, it is customary to filter the cross-correlation waveforms in a few specific frequency bands, and separately invert the filtered data in each band. Such an analysis works well for noise sources that are well-separated in frequency. However, complications may arise if sources are spectrally variable within a chosen frequency band. This can happen, for example, when multiple, localized noise source regions have closely spaced, overlapping frequency spectra. In such a scenario, it can be challenging to estimate the source spectra accurately from the observed data. We present synthetic inversion tests in this section to assess the impact of inaccuracies in this estimate.

For the sake of simplicity, we choose test models of noise source PSD parametrized with |$N_s=2$|⁠, and Gaussian functions (with variable parameters) for the two |$P_l(\omega )$|⁠. The source PSD model of Fig. 1 is split up into two test models, in order to separately consider source distributions within and outside the receiver array (Fig. 3). In each case, the spectral parametrization is with Gaussian functions |$P_1(\omega )$| and |$P_2(\omega )$| peaked at the arbitrarily chosen values of 0.3 and 0.4 Hz, respectively. Given the spatial extent of our computational domain and the dominant wavelengths modelled, we note that the maximum distance of noise sources from the centre of the receiver network is about |$5\lambda$|⁠. Finally, since our stated objective is to isolate the effects of spectrally variable sources, we use a known, homogeneous, wave speed model in all cases. The inability to precisely localize more distant sources, as well as the impact of wave speed model inaccuracies, have already been demonstrated by DSK23.

Figure 3.

Two test models (a1-a4, model A, and b1-b4, model B) of noise source PSD used to perform the synthetic inversion tests with narrow-band sources. Note that the spectral parametrization is the same in both cases (a2,a4 and b2,b4), consisting of two Gaussian functions centred at 0.3 and 0.4 Hz, with a standard deviation of 0.08 Hz. The two models differ in their source strength distributions (a1,a3 and b1,b3), which are derived from the model of Fig. 1.

For each test model, test data are generated by forward modelling and adding Gaussian noise to the resulting cross-correlation waveforms. The added noise has standard deviation equal to 20 per cent of the average RMS value of all the computed cross-correlations. The test data are then subjected to two separate inversions, with |$N_s=2$| and |$N_s=1$|⁠.

In the first type of inversion (⁠|$N_s=2$|⁠), the spectral parametrization of the source PSD model is accurate (same as what is used to generate the test data), and the inversion results are satisfactory (Appendix  A, Figs A1, A2). However, this accurate parametrization is a luxury afforded by our synthetic tests. In practice, one would tend to invert such a narrow-band data set using |$N_s=1$|⁠, therefore in this paper we are more interested in the second type of inversion. In this case, the spectral parametrization is through a single Gaussian function |$P(\omega )$| which must be estimated from the cross-correlation data. The spectrum of the data itself cannot be used as |$P(\omega )$| because the cross-correlation spectra are related to the source spectra via a spatial integral containing the medium Green’s function (eq. 1). We adopt the following procedure to estimate the source spectra:

• take a few different trial |$P(\omega )$| functions guided by the mean spectrum of the data—the data spectrum provides an estimate of the peak frequency (if any) and bandwidth of the sources.

• compute synthetic cross-correlation waveforms in the inversion starting model, for each of these trial functions, by forward modelling eq. (1).

• compare the mean spectra of the synthetics to that of the observed data, and choose the trial |$P(\omega )$| function which yields the best match between computed and observed cross-correlation spectra.

For our synthetic tests, we present results with four different trial functions |$P(\omega )$|⁠, each bearing a relation to the test model parametrization. Figs 4 and 5 show the inversion results for test models A and B, respectively. In both these figures, we see that the different source regions of the test model are not faithfully recovered for any choice of |$P(\omega )$|⁠. In Fig. 4, arguably the ‘best’ inversion (b4) corresponds to the choice of |$P(\omega )$| for which the initial synthetics best match the observations (a4). We note that this |$P(\omega )$|⁠, plot (c4), is a superposition of the two spectra used to parametrize the test model. In Fig. 5, it is harder to pick the best inversion result, and correspondingly, the winning choice of |$P(\omega )$|⁠. We also note from Fig. 5, plots (a1–a4), that the observed spectra do not show a bimodal character (a consequence of the source strength distribution in test model B), thereby strongly justifying the choice |$N_s=1$|⁠.

Figure 4.

Inversion of data from test model A, performed using |$N_s=1$|⁠. The four rows numbered 1–4 correspond to four separate inversions, performed using a Gaussian |$P(\omega )$| peaked at 0.3 Hz (c1), 0.35 Hz (c2), 0.4 Hz (c3) and 0.35 Hz (c4). c2 and c4 differ in the standard deviation of the Gaussian function: in (c1)–(c3), the standard deviation is 0.08 Hz, same as the test models, but (c4) is obtained as a sum of (c1) and (c3). Each (i^th) row shows the mean spectra of the observed and synthetic cross-correlations (ai) and the source strength distribution obtained by inversion (bi).

Figure 5.

Similar to Fig. 4, but for test model B.

Figs 4 and 5 demonstrate the impact of an inherent deficiency in the model parametrization (wrong value of |$N_s$|⁠), on the quality of source inversion. The specific tests considered here demonstrate the limitations of a fixed-spectrum inversion. However the takeaways are more general because similar tests can be conceived (although we have not shown them), where the synthetic data are modelled with |$N_s \ge 3$| but the inversion is performed with |$N_s=2$|⁠, for example. Yet more severe problems with the spectral parametrization (e.g. unmodelled variability in spectral shapes for different source types l), are expected to further diminish inversion quality. The limited tests presented in this section establish the baseline for inaccuracies to be expected in source inversion due to inaccuracies in the spectral parametrization. To this end, we find that the ‘best’ estimates of the source spectra (those which lead to the best inversions) are indeed obtained by matching, in an average sense, the spectrum of the data. The procedure outlined above may be adapted for selecting |$P_l(\omega )$| from a series of trial functions. An alternative is to parametrize the source spectra using a sufficient number of basis functions, as done in Ermert et al. (2021).

4 DISCUSSION

In this paper we have presented a method for ambient noise cross-correlation modelling and inversion at local scales, which accounts for spatiospectral variability in the noise sources. The forward modelling algorithm is similar to what has been used in previous studies, but is currently limited to the acoustic regime. On the other hand, from the perspective of ambient noise source inversion, the chosen parametrization of the source PSD allowed us to focus on a specific issue—that of inaccuracies in source spectra estimated from observed data. We have undertaken synthetic tests to investigate the severity of the issue. The results of these tests can be abstracted beyond the two specific (and arbitrarily chosen) test models used in this study, to draw the general conclusion that accurate spectral parametrization of the sources is important.

It may be encouraging to note that even with inaccurate spectral parametrization, one can recover true source distributions at least partially. However, it should be borne in mind that the inversion inaccuracies illustrated in Section 3 are the bare minimum of what should be expected in practice. Working with real data, we expect further complications due to various factors such as more elaborate source spectra (modelled with |$N_s > 2$| ), velocity model inaccuracies, unmodelled distant sources and also pre-processing of the raw data in order to produce the cross-correlation waveforms (which constitute the input to our source inversion algorithm).

Nonlinear processing schemes such as one-bit normalization or spectral whitening can produce artefacts in the noise correlations (Fichtner et al. 2020), so they must be used with caution when combined with inversion schemes that are based on modelling the correlations. Especially for an inversion technique like ours that makes use of cross-correlation amplitudes, one needs to be mindful of the impact of pre-processing, on noise correlation amplitudes. Fortunately, the technique relies on relative rather than absolute amplitudes, therefore the exact nature of applied processing does not matter as long as it is consistent across the receiver array (unless there is strong spatial variability in attenuation characteristics, which is unlikely at local scales). However, we concede the possibility of sources from different frequency bands being obscured by processing steps like spectral whitening. The consequences of using heavily processed correlations to perform source inversions by our method can only be fully understood through further analysis, which is beyond the scope of this study.

ACKNOWLEDGMENTS

Numerical simulations for forward modelling were performed using the open-source finite difference wave-equation solver, Devito (Louboutin et al. 2019). We thank Piero Poli, Laura Ermert and an anonymous reviewer for their insightful comments which helped to improve the manuscript and make its message clearer. This research was supported by the Department of Science and Technology, Government of India, via an INSPIRE Faculty grant to AD (IFA-19 EAS-78). PPB would like to acknowledge support from MoES grant number MoES/INDO-SWISS/4/2022-PC-I. BS acknowledges funding from the Institute of Eminence (IoE) grant at IIT Bombay.

DATA AVAILABILITY

No new data were generated or analysed in support of this research. The source inversion code used in this work is available on github: https://github.com/arjundatta23/cc_kern_inv.

APPENDIX A: ADDITIONAL INVERSION RESULTS

Figure A1.

Inversion of data from the test model of Fig. 3(a), performed using |$N_s=2$|⁠. (a) Comparison between the power spectrum of the data being inverted, |$|C^{obs}(\omega )|^2$| and the spectral parametrization chosen for inversion. The magnitude of the |$P_i(\omega )$| functions is set before inversion begins, by fitting the average energy of the data. (b) the result of inversion-inverted source strength distributions (b1 and b3) corresponding to the chosen source spectra (b2 and b4, respectively). (c) power spectrum of the data, compared with that of the initial and final synthetics from inversion.

Figure A2.

Inversion of data from the test model of Fig. 3(b), performed using |$N_s=2$|⁠. Subplots (a)–(c) are similar to those of Fig. A1.