Rapid finite-frequency microseismic noise source inversion at regional to global scales


 Ambient noise cross-correlations can be used as self-consistent observables, opening novel possibilities for investigating ambient noise sources. To optimize the forward-modelling of global ambient noise cross-correlations for any given distribution of noise sources in the microseismic frequency range up to 0.2 Hz, we implement (i) pre-computed wavefields and (ii) spatially variable grids. This enables rapid inversions for microseismic noise sources based on finite-frequency source sensitivity kernels. We use this advancement to perform regional and global gradient-based iterative inversions of the logarithmic energy ratio in the causal and acausal branches of microseismic noise cross-correlations. Synthetic inversions show promising results, with good recovery of the main dominant noise sources of the target model. Data inversions for several consecutive days at the beginning of October 2019 demonstrate the capability of inverting for the spatio-temporal variations of the sources of secondary microseisms in the ocean. This paves the way for daily ambient noise source inversions which could help improve full-waveform ambient noise tomography and subsurface monitoring methods.

or matched-field processing (e.g. Gal et al. 2018) which require different assumptions, such as plane wave propagation, or a sufficiently dense array. To properly represent ambient wave propagation through the 3-D Earth, several authors adopted the direct numerical modelling of noise correlations for heterogeneous noise sources (e.g. Nishida & Fukao 2007;Tromp et al. 2010;Hanasoge 2013a;Ermert et al. 2017;Sager et al. 2018a;Datta et al. 2019;Xu et al. 2019), based on concepts originally developed in helioseismology (Woodard 1997;Gizon & Birch 2002). This approach naturally yields noise source sensitivity kernels that may be used to infer the spatial distribution of sources, while honouring the physics of wave propagation through a 3-D heterogeneous medium (Ermert et al. 2017(Ermert et al. , 2020. Another approach that has been shown to be effective to model ambient noise correlations is normal mode summation (e.g. Nishida & Fukao 2007;Gualtieri et al. 2013).
To reduce computational cost, especially when higher frequencies are involved, Ermert et al. (2017) proposed an implementation based on pre-computed wavefields from numerical wavefield solvers. This is particularly effective in an inversion framework, because the wavefield only has to be computed once. Due to the chosen measurement, namely the logarithmic energy ratio which is largely insensitive to unknown 3-D Earth structure (Sager et al. 2018a), we can justify using simpler 1-D models like PREM (Dziewonski & Anderson 1981) as underlying velocity model.
We aim to efficiently forward model noise cross-correlations on a global scale for any noise source distribution for frequencies up to 0.2 Hz. To achieve this, we implement a spatially variable grid alongside pre-computed wavefields to reduce the number of model parameters and thus the computational cost. We assess the accuracy of the reduced parameterisation by carrying out benchmark tests with a globally homogeneous parameterisation and a different wave propagation solver. Furthermore, the aforementioned sensitivity kernels allow us to easily implement a gradient-based iterative method to invert for the noise source distribution. Following a brief review of the underlying equations, we perform synthetic inversions to test the method under idealistic circumstances. Finally, we present real-data inversions for several consecutive days.

FORWARD MODELLING
The concept of modelling cross-correlations for arbitrary noise source distributions originated in helioseismology (Woodard 1997) and has been modified for the use on Earth by several authors (Tromp et al. 2010;Hanasoge 2013b;Fichtner 2014). To provide the necessary context, we give a short derivation of the cross-correlation wavefield equations, similar to Ermert et al. (2017Ermert et al. ( , 2020.
Subsequently, we describe our approach to make the computation feasible for global, high-frequency problems using pre-computed wavefields and spatially variable grids.
Rapid finite-frequency microseismic noise source inversion at regional to global scales 5

Continuous and discretised modelling equations
The ambient noise field is a superposition of elastic waves excited by a distributed source N n (ξ, ω) with n−component where ξ denotes the location on the Earth's surface ∂⊕. We henceforth omit the frequency dependency ω in the interest of condensed notation. The i-component of ground motion, u i (x), at receiver position x can therefore be expressed as a convolution of the Green's function G in (x, ξ) for some suitable Earth model and the noise sources N n (ξ): (1) The Einstein summation convention for repeated indices is implied. The ensemble-averaged frequency-domain correlation of two wavefield recordings at positions x 1 and x 2 is given by In Eq. (3), · denotes the averaging over several time intervals and * the complex conjugate. Thus, to model a cross-correlation wavefield, we need to model the Green's functions G in (x, ξ) and parameterise the noise sources N n (ξ). The temporal averaging has to be sufficiently long to ensure that purely stochastic fluctuations (e.g., instrumental noise) become irrelevant compared to the ambient noise signal that we wish to exploit. For modern broadband seismometers, as used in the application part of this work, instrumental noise levels are typically insignificant compared to ambient noise levels, especially within the microseismic frequency band. It is therefore safe to assume that instrumental noise is practically irrelevant when averaging time scales are an order of magnitude longer than the longest period considered in the correlation analysis. In practice, the averaging length is more controlled by the need to obtain correlation time series that are sufficiently simple, e.g., close to a Green's function, in order to make intuitively interpretable measurements that allow us to formulate a well-behaved inverse problem.
The correlation of the two noise sources can be described by their location-dependent power-spectral density (PSD) S nm (ξ), where we made the common assumption that neighbouring noise sources at ξ 1 and ξ 2 are temporally 6 Jonas K.H. Igel, Laura A. Ermert, Andreas Fichtner uncorrelated (e.g. Tromp et al. 2010;Hanasoge 2013b;Fichtner 2014). Upon inserting (4) into (3) the ensemble correlation condenses to Eq. (5) allows us to evaluate a cross-correlation wavefield C ij at position x 1 with a reference station at position x 2 for any noise source distribution S nm (ξ).

Pre-computed Wavefields
Since we are solely interested in the noise source distribution S nm (ξ) and not in changes of Earth structure, we can take advantage of pre-computed Green's function databases to reduce the computational cost of the correlation wavefield modelling. This approach can be further justified by choosing a suitable correlation waveform measurement that is largely insensitive to unknown 3-D Earth structure. This topic will be discussed in more detail in section 3.
Ambient noise propagation is an inherently global phenomenon, even when stations used for data analysis cluster at regional scale (e.g. Hillers et al. 2012;Retailleau et al. 2017;Sager et al. 2018a).
Several wave propagation solvers have been developed for global wave simulations -for instance, the finite-difference implementation of Igel & Weber (1996), SpecFEM (Komatitsch & Tromp 2002 To model ambient noise we compute a Green's function database with the help of reciprocity where each receiver is taken as a point source and the corresponding seismograms are recorded at all grid points. We then approximate the frequency-domain correlation wavefield by a sum over Green's functions multiplied with their corresponding noise source weight S(ξ k ), where ∆ξ k is an approximation of the finite surface element, and n g is the number of grid points (discrete noise sources). The computational cost of (6) scales with the number of grid points, Rapid finite-frequency microseismic noise source inversion at regional to global scales 7 suggesting that the noise source parameterisation should be optimised for efficient noise correlation modelling. A possible optimisation is the topic of the following section.

Spatially Variable Grid
We reduce the number of grid points by implementing a spatially variable grid. This is possible thanks to the shape of the noise source sensitivity kernels, which become broader linearly with increasing distance from the station pair (Hanasoge 2013b), thereby requiring less spatial resolution for the parameterisation. Furthermore, attenuation causes distant noise sources to be of lesser importance, additionally avoiding the refocusing of sensitivity at station antipoles. To obtain a parameterisation that is both flexible and simple, we implement spatially variable Gaussian grids, some examples of which are visualised in Fig. 1. The radial distance dφ between concentric circles of grid points on a sphere is described by where the minimum and maximum incremental grid point distance is represented by dφ min and dφ max , the radius of the area of high grid point density by σ, and the steepness of the increase in grid point distance by β. The iteratively increasing variable i, representing the index of each circle, is used to create the distribution up to the point where the sum of all dφ is 180 • , that is, the whole globe is accounted for by the concentric radial circles. We set the distance between grid points on each circle to roughly the radial distance, to ensure approximate homogeneity of the grid. Additionally, the geographical centre of the grid can be shifted to any location. We implement further options to increase adaptivity to heterogeneous noise sources. To account for multiple strong noise sources, we allow multiple areas with higher spatial resolution. This is achieved with a coarse background grid from which grid points are removed, and smaller-scale denser grids added. If prior information about the noise source distribution is available (for instance from ocean surface pressure maps that are related to secondary microseimic sources, see Ardhuin et al. 2011) the spatial resolution can be increased in areas of higher PSD. We choose the threshold value of the PSD above which additional sub-grids with higher spatial resolution are added in a given radius around the sources. This ensures that the structure of the strong noise sources is considered, whereas areas of weak noise sources have a sparser grid to save computational cost. An example of such an adapted grid is shown in Fig. 2. Here, the input model used for correlation computations is an ocean surface pressure map by Ardhuin et al. (2011).
To obtain suitable surface elements ∆ξ k for spatially variable grids, we compute the Voronoi cell surface area for each grid point using the Python package SciPy (Virtanen et al. 2020). This ensures that an accumulation of grid points does not lead to stronger noise sources in that particular area.
Combining pre-computed wavefields and spatially variable grids enables us to efficiently model global cross-correlations for any given noise distribution.

Benchmarks
To quantify the potential computational savings of the adaptive gridding, we compare correlations computed with spatially variable and globally homogeneous grids. Considering only oceanic sources, Rapid finite-frequency microseismic noise source inversion at regional to global scales 9 the total number of points ranges from 1,000 to 80,000. Since the 24 stations selected for the benchmarks are located around the North Atlantic, we increase resolution of the spatially variable grid in that region. Two example grids for roughly 20,000 grid points and the station locations can be seen in Fig The underlying Earth model is PREM (Dziewonski & Anderson 1981).
As a measure of correlation waveform fit we employ, in the interest of simplicity, the L 2 distance between the correlations computed using the test grids C test i and the reference grid C ref i : where N denotes the number of station pairs. Fig. 3 shows how the L 2 distance evolves with an increasing number of grid points. The spatially variable grid initially has a lower misfit and also converges faster than the homogeneous grid. For instance, for the homogeneous grid we require around 12,000 grid points to reach a misfit of 0.1, whereas the spatially variable grid only requires 4,000 grid points. Hence, in this case, we may reduce the number of model parameters by a factor of grid points for a homogeneous grid but only 4,000 grid points with a spatially variable grid, meaning that we may reduce the number of grid points by a factor of 3.
3 and still expect results to be nearly identical. Depending on the chosen threshold, this parameter reduction varies between 2.5 and 4.5.
While this example illustrates the potential of the variable grid approach, its efficiency is clearly application-specific and needs to be assessed carefully on a case-by-case basis prior to a large-scale inversion that requires numerous forward simulations. Particularly global inversions poorly benefit from such a spatially variable grid as most areas are high-sensitivity regions. However, the grid could still be adapted if strong noise source areas are known beforehand.
An alternative to the direct discretisation of Eq. (6) combined with a Green's function database is on-the-fly computation of the correlation wavefield source G jm (x 2 , ξ)S nm (ξ) on the grid of a numerical wave propagation solver, and its use in a subsequent forward simulation (Tromp et al. 2010;Fichtner 2014;Sager et al. 2020). This approach trades the larger storage requirements of the Green's function database for larger computational requirements in an iterative inversion, where correlations must be computed for a sequence of noise source updates. Comparing these independent approaches provides another useful plausibility check on the quality of the correlation functions. This is described in more detail in supplementary section S1.
To constrain the distribution of noise sources, we require information on how changes in the noise source PSD, S nm (ξ), affect the fit between observed and synthetic correlations. Applying adjoint techniques allows us to reduce the computational cost to obtain the source sensitivity kernels by using the pre-computed Green's function databases. More details and in-depth derivations for finite-frequency sensitivity kernels can be found in Tromp et al. (2010), Fichtner (2014) for the components i and j of a station pair with source PSD components n, m is where f ij (x 1 , x 2 ) is the measurement-dependent adjoint source which excites the adjoint wavefield at the receiver position x 1 . Since we have already pre-computed the Green's function database, the kernel computation does not require any additional simulations.
While there are numerous possible measurements that may be used to quantify differences between observed and synthetic correlation waveforms, we limit ourselves to the causal-acausal logarithmic energy ratio, defined as with the time window w(τ ) at lag τ centred around the expected surface-wave arrival time of the correlation C(τ ). The observed energy ratio A obs is defined accordingly. The outstanding advantage of the asymmetry measure A is its relative insensitivity to unknown 3-D Earth structure (Sager et al. 2018a), which adds another justification to the pre-computed wavefield approach where we keep the Green's functions constant during the inversion process. Discrepancies in the surface wave arrival times due to unknown 3-D Earth structure can be counteracted by increasing the expected surface wave arrival time window in the measurement. This ensures that the energy calculation still includes the main wavelets even when the given surface wave velocity is not entirely accurate. Due to this characteristic of the measurement we currently compute the Green's functions using the 1-D PREM (Dziewonski & Anderson 1981) as underlying velocity model.
Based on Eq. (10), we define the L 2 misfit as a sum over asymmetry differences measured at all station pairs, Finally, we derive the equation for the measurement-dependent adjoint source from Eq. 10 and Eq.
11 to obtain A more detailed derivation of the adjoint source for this specific choice of measurement and misfit can be found in Ermert et al. (2016).
An illustration of the signal energy measurement and the corresponding sensitivity kernel can be seen in Fig. 4. To obtain the example sensitivity kernel we computed an artificial observed correlation for one station pair with a dominant noise source south of Iceland. The synthetic cross-correlation is based on a homogeneous noise distribution in the oceans. As expected, the sensitivity kernel is largely negative in the North Atlantic, suggesting that an increase in noise source PSD in that region will decrease the misfit.

MICROSEISMIC NOISE SOURCE INVERSION
We perform several synthetic and real-data ambient noise source inversions using vertical-vertical-component correlations and a gradient-based iterative scheme to minimise the misfit in Eq. (11), similar to the method applied to the Earth's hum by Ermert et al. (2017). Since vertical-component data mostly constrain vertical-component sources, we restrict all inversions to S zz , setting the remaining components of the PSD tensor to 0. Prior to showing results we describe the overall setup, which is the same for both synthetic and real-data applications.
The underlying Green's function database is the result of AxiSEM (Nissen-Meyer et al. 2014) simulations with frequencies up to 0.2 Hz. This allows us to invert for secondary microseismic sources between 0.1 and 0.2 Hz. Velocity, attenuation, and density structure are based on PREM (Dziewonski & Anderson 1981), and the noise source model is parameterised using a spatially variable grid with denser grid point spacing in the North Atlantic. Due to the shape of source sensitivity kernels, it is beneficial to use stations around the area of interest to obtain good coverage.
Thus, we choose stations surrounding the North Atlantic in North America and Europe.
An illustration of ray coverage and station sensitivity for a set of station locations is presented in Fig. Rapid finite-frequency microseismic noise source inversion at regional to global scales 13 To further decrease the computational cost, we choose the observed cross-correlations for the inversion based on a modified signal-to-noise ratio. This is computed by dividing the maximum amplitude in the expected arrival time windows by the standard deviation of the whole time series (see Fig. S2 for visualisation). Since the computational cost directly scales with the number of station pairs, we set a minimum signal-to-noise threshold of 3 to filter out cross-correlations with little information in the surface wave arrival time windows. This is usually the case if the dominant noise sources are not in-line with the station pair.
Since the secondary microseisms are solely generated by ocean waves, we take a spatially homogeneous noise distribution in the oceans as the initial model. We parameterise the noise sources using a spatially variable grid with roughly 14,000 grid points, removing all grid points on land.
Since our stations surround the North Atlantic we increase the spatial resolution there. Furthermore, we parameterise the frequency spectrum of each noise source by the peak amplitude of a Gaussian with centre frequency of 0.15 Hz and standard deviation of 0.05 Hz. Observed and synthetic correlations are filtered accordingly. surrounding the Northern Atlantic. Due to the shape of the sensitivity kernels, the rays are not between but behind the station pair (top right inset). The station sensitivity is computed by taking the sum of the absolute values of all sensitivity kernels. Regions where station sensitivity is below 1% of the maximum are masked in order to indicate that information to constrain noise sources is likely to be insufficient. Since we take attenuation into account, station sensitivity does not refocus at the station antipoles, in contrast to the simplified rays shown in the top panel.
Rapid finite-frequency microseismic noise source inversion at regional to global scales 15 To minimise the asymmetry misfit, we employ an iterative steepest-descent method, where correlations, misfits, and kernels are re-computed for every model update. Prior to any update, we apply widely used preconditioning, including 95th percentile clipping and Gaussian smoothing of the kernels. The standard deviation of the Gaussian smoothing filter decreases during the inversion from 4 • to 1.5 • in order to initially avoid local minima and later permit more spatial detail. These values are largely the result of trial and error based on a series of synthetic inversions. Since PSD cannot be negative, we introduce a non-negativity constraint by setting all sources with negative PSD to 0 after each model update.
In the following paragraphs we present synthetic inversions intended to test the inversion scheme for different noise distributions and station locations under idealised conditions. Furthermore, we perform real-data inversions for several consecutive days.

Synthetic inversions
As a first synthetic application of the inversion method we use a significant wave height map from the WaveWatch III model (Tolman & Chalikov 1996) as the input noise source PSD distribution for the computation of artificial observed correlations. We choose 217 stations surrounding the North Atlantic to perform the inversion, resulting in roughly 12,500 measurements excluding all correlations with a signal-to-noise ratio below 3. Executing the approach outlined above, the misfit drops by ≈ 75% after the first iteration, and by 92% after ten iterations. The observation that the misfit plateaus quickly is typical also for other inversions presented later, and it leads us to generally stop the misfit minimisation after 10 iterations.
A summary of this first inversion is presented in Fig. 6, where we do not plot areas where the station sensitivity (seen in Fig. 5) is below 1%. Comparing the artificial distribution with the final inversion model, we see that the dominant noise sources off the coasts of Ireland, France and Spain can be recovered. Additionally, the spatially smaller sources on the southern coast of Greenland and between Iceland and Norway are included, even though the spatial resolution is evidently limited by the station geometry.
Similar to the first inversion, we use data from ocean surface pressure maps by Ardhuin et al. (2011) as input noise distribution for a second synthetic test. As can be seen in Fig. 7, these maps have smaller-scale noise sources compared to the significant wave height maps. Though the inversion is, as expected, unable to resolve such small details, it provides a qualitative assessment of the smearing length scales in our region of interest. After 10 iterations and a misfit reduction of 84%, the noise sources east of Iceland and north of Scandinavia can be distinguished, whereas other features are

Real-data inversions
Finally, we apply our method to real-data correlations of ambient noise in order to infer the spatio-temporal behaviour of secondary microseism sources. To download and cross-correlate the noise data we use the mass downloader from the Python-based toolbox ObsPy .
Motivated by the presence of Hurricane Lorenzo in the North Atlantic, we choose the first six days of October 2019 as time frame. Furthermore, there were no notable earthquakes during that time period, thereby allowing us to apply only basic linear processing that avoids potentially nonphysical effects . First, we segment the data into 1-hour windows and remove the instrument response to obtain true displacement data. Subsequently, we downsample the data to 1 Hz and compute daily correlations by stacking individual correlations of the 24 time windows. Once the observed correlations are computed, we perform the inversions for six consecutive days.
Since we only use one day of ambient noise data, the signal-to-noise ratio in the correlations is small, particularly compared to the synthetic correlations. To ensure that we are actually inverting for a signal and not just a random time series, we ignore correlations with a signal-to-noise ratio below 3.
Additionally, we add random noise to the synthetic correlations during the inversion to introduce a water level for the measurement which ensures we do not have 0 energy in a window and make the synthetics resemble observed data more closely. The noise level added is based on the maximum amplitude of each individual correlation and is kept constant throughout the inversion. Synthetic tests (see Fig. S6 and Fig. S7) show that this approach considerably improves the inversion model if the signal-to-noise ratio in the observed data is low. More analysis of different noise levels for the synthetic correlations and its effect on the inversion framework can be found in the supplementary material.
The final inversion models after 10 iterations are summarised in Fig. 8, where the models are normalised with the maximum PSD of the 6 day period. The models show the spatio-temporal variations of the source distribution of secondary microseisms in the wider North Atlantic region.
The average misfit reduction is approximately 28%.
Theoretically, the ocean surface pressure maps by Ardhuin et al. (2015) should correspond to the sources of the secondary microseisms when the ocean site effect is taken into account (Longuet-Higgins 1950;Gualtieri et al. 2013Gualtieri et al. , 2014Ardhuin et al. 2019). Fig. 9 shows these maps for the six consecutive days in the same frequency range of 0.1 -0.2 Hz, with the same smoothing filter applied as for the noise correlation inversion. In these maps, the path of Hurricane Lorenzo is visible as it moves north towards Greenland.
Comparing the final inversion models with the ocean surface pressure maps, we observe that the areas of dominant noise sources are similar throughout the 6 day period. The source strengths in the North  Rapid finite-frequency microseismic noise source inversion at regional to global scales 19 Atlantic gradually increase until the 4th October 2019. This corresponds to the day that Hurricane Lorenzo dissipates, after which the source strengths start to decrease again. We also observe some dominant noise sources in marginal seas like the Caspian Sea. However, we hypothesise that these sources have no real physical meaning as they originate from the parameterisation of the inverse problem. Correlations with a low asymmetry ratio correspond to a homogeneous noise distribution, i.e. the same source strength on both sides. Due to the lack of grid points on land this results in the sources in the marginal seas to overcompensate with very high source strengths.
A representative collection of correlations is shown in Fig. 10  attempts to match causal and acausal energy, which contains less information but is more robust.
Correlations with this lack of asymmetry automatically have less weight as the measurement is small compared to correlations with a clear asymmetry. By adding noise to the synthetics we avoid big discrepancies between the observed and synthetic measurements which would cause unwanted effects due to the parameterisation of the noise sources.

DISCUSSION
In the previous sections, we presented a global-scale finite-frequency ambient noise source inversion in the secondary microseismic frequency range. This is made possible by optimising the parameterisation and reusing pre-computed Green's function to greatly reduce the computational cost. Compared to other studies, we are able to circumvent the most common assumptions like the reduction to simple rays, the quasi-random nature of noise sources, and the plane wave assumption commonly used in beamforming, by using finite-frequency sensitivity kernels. By doing so, we properly account for global wave propagation with visco-elastic attenuation. Additionally, we only apply linear processing steps to the data reducing the risk of numerical artifacts in the cross-correlations. Thus, this method can be applied to any existing, linearly processed ambient noise cross-correlation data set. In the following, we provide a more detailed discussion on several key issues, such as measurements, 3-D structure effects, and the effective nature of the results. the causal-to-acausal energy ratio is matched (bottom). This behaviour is expected as we do not attempt to fit waveforms but only the less informative but more robust energy ratio.

Measurements and time windows
Our measurement, logarithmic energy ratio, relies on the existence of surface wave arrivals in the expected time window. However, when the dominant noise sources are not in-line with the receiver pair, spurious (but still physically meaningful) arrivals may appear outside the expected time window.
These spurious arrivals are not taken into account so far. A more flexible measurement window, e.g., where we correlate for the time shift as well, could further improve the method and contribute more measurements to the inversion. Additionally, we only consider the vertical component of ground motion. Taking measurements on the horizontal components has been shown to improve results (Xu Rapid finite-frequency microseismic noise source inversion at regional to global scales 21 et al. 2019) and should be implemented in the future. Improving the measurement and using more information from the correlations is part of current research.

3-D Earth structure, oceans and effective noise sources
While the proposed method is in principle independent of the Earth model used to compute Green's functions, we work with radially symmetric Earth models in order to reduce computational cost to a feasible level. Even though it has been shown that the measurement of logarithmic energy ratios is robust against unknown or diregarded 3-D Earth structure (Sager et al. 2018a), this simplification may still affect the inversion results. In the same context, it is important to keep in mind that all inversion results -not only for the specific inverse problem considered here -a relative to generally unavoidable simplifications of forward modelling physics. In this sense, the source distributions that we obtain are an effective physical quantity. A well-known and related example are earthquake moment tensors, which are effective point-localised simplifications of a more complicated finite-rupture process, and which are mostly estimated on the basis of 1-D Earth models without fluid ocean layers. Despite being effective, the resulting sources are useful, e.g., for tectonic studies.

Temporal averaging and temporal resolution
Since the estimation of the source distribution is based on temporal averaging, we have to ensure that the length of the finite-time cross-correlations is sufficient to neglect stochastic fluctuations such as instrumental noise. Though we did not investigate the lower limit of temporal resolution, the low levels of instrumental noise in modern broadband seismometers, ensure that averaging over time scales that exceed the longest period by an order of magnitude, eliminates instrumental noise for all practical purposes.
The actual limitation of temporal resolution comes from the measurement process and our intuition for the inverse problem. These rest on the assumption that approximations of surface waves appear within the measurement windows. Hence, the averaging should be sufficiently long to ensure that surface waves become visible. Though this is technically not required, it still allows us to interpret the measurements and to set up the inverse problem accordingly. In our specific application, averaging over one day typically produces surface wave approximations that can be used in the logarithmic energy measurements.

Coverage and station locations
The station coverage is of great importance for ambient noise studies as this determines the area of sensitivity. If we include continental grid points (see Fig. S3) we see inversion artifacts on the continents due to the lack of station coverage (and thus sensitivity) where we do not expect any sources. Choosing the station locations based on a more elaborate scheme using optimal design could further improve the resolution of the method whilst possibly decreasing the computational cost.

Computational requirements
For a large number of stations and global inversions, access to high performance computing facilities is useful as it significantly decreases the compute time due to the embarrassingly parallel nature of the code. A wavefield conversion with an AxiSEM wavefield including all frequencies up to 0.2 Hz for 150 stations and about 13,000 grid points takes between 1 and 4 hours per station depending on the sampling rate and correlation length. If we keep the station locations and the grid the same, this computation only has to be done once. A consequent inversion with 10 iterations using a subset of 3,100 cross-correlations from the 150 stations takes about 50 minutes on 600 cores. Smaller inversions on a local scale can easily be run on small clusters using publicly available pre-computed AxiSEM wavefields from, e.g., the Syngine repository (IRIS 2015; Krischer et al. 2017).

CONCLUSIONS AND OUTLOOK
We efficiently invert for the seismic noise distribution of the secondary microseisms by computing finite-frequency sensitivity kernels using a logarithmic energy ratio measurement of noise cross-correlations. The computational cost is reduced by implementing a spatially variable grid to optimise the parameterisation, and by pre-computed wavefields from which Green's functions for forward modelling and kernel computations can easily be extracted. Benchmark tests show that these simplifications have no effect on the measurement we use. Several advantages arise from our approach: (i) we properly account for global wave propagation with visco-elastic attenuation; (ii) no assumptions on the wavefield (e.g. equipartioning or plane waves) have to be made; (iii) our Rapid finite-frequency microseismic noise source inversion at regional to global scales 23 measurement is robust to unknown 3D Earth structure; (iv) the method can be applied to any noise cross-correlation data set; and is ( Synthetic inversions show promising results and we are able to resolve the dominant noise sources of the target model. The spatial resolution is limited to a few hundred kilometres due to the spreading of the sensitivity kernel as we move away from the station pairs. Data inversions for several consecutive days demonstrate that we are able to see the spatio-temporal variations in the sources of secondary microseisms with our method. Further analysis should be done on more elaborate window picking and choosing the station locations based on an optimal design scheme. A synthetic inversion with globally distributed station locations (see Fig. S4) additionally demonstrates the potential for efficient noise source inversions with sensitivity all around the globe.
Due to the efficiency of the method, we pave the way for global daily ambient noise source inversions. Publicly available daily noise source maps could help to improve methods in full waveform ambient noise tomography and near real-time subsurface monitoring.

ACKNOWLEDGMENTS
The authors would like to thank Fabrice Ardhuin for providing the ocean surface pressure maps and

Data and Code Availability
The code and data are available from the authors upon request.