SANS: Publicly Available Daily Multi‐Scale Seismic Ambient Noise Source Maps

Seismic ambient noise sources have received increased attention recently, creating new possibilities to study the Earth's subsurface and the atmosphere‐ocean‐solid Earth coupling. In efforts to locate such noise sources using nonlinear finite‐frequency inversions, methodological developments such as pre‐computed wavefields and spatially variable grids were necessary. These make inversions feasible for the secondary microseismic sources in a frequency range up to 0.2 Hz on a daily basis. By obtaining a starting model for the inversion using Matched Field Processing (MFP) we are able to further steer the inversion toward acceptable global noise source models and improve the final result. Analysis of 1 year of daily inversions shows the seasonal variations of the secondary microseisms and their dependence on the atmosphere‐ocean‐solid Earth coupling due to storm‐induced ocean waves. We present a web framework, seismic ambient noise sources (SANS, sans.ethz.ch), where daily regional‐ to global‐scale SANS maps for the secondary microseisms are made available to the public. This enables the implementation of time‐variable noise source distributions into full‐waveform ambient noise tomography and time‐dependent subsurface monitoring methods. Additionally, it can act as a reference for other studies that locate secondary microseismic sources on a larger scale.

open ocean when two opposing ocean waves overlap to create a vertical pressure wave that causes tiny displacements at the ocean bottom (Longuet-Higgins, 1950;Longuet-Higgins & Ursell, 1948). In our study we focus on secondary microseismic sources.
In theory, cross-correlation functions approach the Green's functions for a homogeneously distributed, uncorrelated, random noise source distribution and an equipartitioned wavefield. Many ambient noise interferometry studies assume that the noise source distribution is sufficiently homogeneous for the cross-correlations to converge to Green's functions (e.g., Nakata et al., 2019;Snieder & Wapenaar, 2010;Wapenaar, 2004;Wapenaar & Fokkema, 2006;Weaver et al., 2009). However, several studies have shown that the omni-present ambient noise wavefield changes on a daily basis (e.g., Ardhuin et al., 2015;Bertelli, 1872;Longuet-Higgins, 1950) and the cross-correlation and Green's function diverge if more realistic constraints-such as global or local energy and directionality constraints-are implemented into the modeling (Halliday & Curtis, 2008;Tsai, 2009;Tsai & Sager, 2022). The heterogeneity of noise source distributions can have a significant effect on travel times, particularly for monitoring applications (Delaney et al., 2017;Zhan et al., 2013). Prior knowledge of the noise source locations can help distinguish if arrival time changes in the cross-correlations are due to changes in the noise source distribution or subsurface velocities.
Inspired by work in helioseismology (Woodard, 1997), recent works introduced the direct numerical modeling of noise cross-correlations for any heterogeneous noise source distribution on Earth (Datta et al., 2019;Ermert et al., 2017;Fichtner, 2014;Hanasoge, 2013b;Tromp et al., 2010). This has resulted in several studies using adjoint techniques (e.g., Fichtner et al., 2006) and sensitivity kernels (e.g., Fichtner, 2014;Tromp et al., 2010) to invert for the seismic ambient noise source (SANS) distribution for different frequency ranges on various scales Igel et al., 2021;Xu et al., 2019Xu et al., , 2020. Expanding on these adjoint and sensitivity kernel techniques, Xu et al. (2019) and Bowden et al. (2021) showed that certain beamforming algorithms are mathematically similar to the first iteration of nonlinear finite-frequency inversions.
The direct forward-modeling of ambient noise cross-correlations allows us to circumvent common assumptions in ambient noise studies-for example, wavefield equipartitioning and a quasi-random noise source distribution-that are necessary for Green's function retrieval (e.g., Sánchez-Sesma & Campillo, 2006;Shapiro & Campillo, 2004;Shapiro et al., 2005;Wapenaar & Fokkema, 2006). Additionally, full-waveform ambient noise tomography methods are capable of directly implementing information about the noise source distribution . Recent developments have made the computation of cross-correlations for ambient noise source inversions more efficient by using spatially variable grids and pre-computed Green's function databases (Ermert et al., 2020;Igel et al., 2021). This allows us to rapidly invert for the noise source distribution of the secondary microseismic sources between 0.1 and 0.2 Hz on a regional to global scale with reasonable computational cost by taking advantage of high performance computing (HPC) resources.
As with other inversion methods, a good initial model is advantageous to ensure convergence to an appropriate final model. In this paper, we introduce the combination of a beamforming method, namely Matched Field Processing (MFP) (Bowden et al., 2021;Gal et al., 2018;Xu et al., 2019), with nonlinear finite-frequency inversions to improve the final inversion result. MFP is capable of locating ambient noise sources by considering sources anywhere in a given domain and estimating the travel times for any station pair to obtain a map of the noise source "power" from the cross-correlations. By producing an initial model with MFP we introduce additional information to the inversion workflow and help to avoid local minima during the subsequent iterations. Since both methods, MFP and nonlinear finite-frequency noise source inversions, have previously been described individually in detail, we refer the interested reader to earlier publications for more in-depth derivations and explanations (Bowden et al., 2021;Ermert et al., 2020;Igel et al., 2021).
Building on these various developments, we present a web framework to make daily SANS maps available to the public (sans.ethz.ch). In-depth knowledge of the ambient noise source distribution has the potential to help improve ambient noise tomography and imaging methods; particularly to ensure that changes in the subsurface are not confused with the spatio-temporal variations of the microseismic noise source distribution. The provision of daily noise source maps is a step toward the incorporation of time-variable noise sources into full-waveform ambient noise inversions, and it adds a new component to near real-time analysis of atmosphere-ocean-solid Earth coupling.

Methodology
In the following section, we will explain the main steps of the two methods: nonlinear finite-frequency inversions and MFP. Despite the differences in the approaches taken, Bowden et al. (2021) show that these methods are well connected. Both have their advantages and disadvantages: constant velocity MFP is an efficient, data-driven approach that works on any cross-correlation data but does not model wave propagation properly and does not allow for an iterative optimisation. An inversion is computationally more expensive but-in contrast to MFPmodels the wave propagation more accurately and allows us to account for the nonlinearity by using an iterative approach.
More importantly, an inversion allows for prior knowledge to be implemented. Hence we use the more efficient, data-driven MFP algorithm to compute a starting model for the nonlinear finite-frequency inversion, to avoid local minima and accelerate the convergence toward an acceptable model. Both methods rely on the fact that vertical-component seismic ambient noise data in the frequency range of 0.1-0.2 Hz are dominated by surface waves. As we solely use vertical-component data, we focus on the Rayleigh wave component of the secondary microseisms.

Nonlinear Finite-Frequency Inversion
The inversion method is based on a concept from helioseismology (Woodard, 1997) which enables the direct modeling of cross-correlations for any noise source power-spectral density (PSD). The work has been adapted for applications to Earth by several authors (e.g., Ermert et al., 2017;Fichtner, 2014;Hanasoge, 2013a;Tromp et al., 2010) with some additional implementations of pre-computed wavefields (Ermert et al., 2020) and spatially variable grids (Igel et al., 2021) to improve efficiency and make inversions feasible for higher frequencies up to 0.2 Hz. In the following section we will provide a short overview of the gradient-based iterative inversion method. For more details, the reader is referred to the aforementioned publications.

Cross-Correlation Modeling
The following equation allows us to forward model the cross-correlation wavefield  , for two stations at locations x 1 and x 2 , for an arbitrary noise source PSD S nm , at points ξ on the Earth's surface ∂⊕, using the Green's functions G, in the frequency domain Fichtner, 2014;Igel et al., 2021;: We imply the Einstein summation convention for repeated indices, and * indicates the complex conjugate. Each noise source is modeled as a vertical point source and the frequency spectrum of each noise source PSD S nm is represented by a Gaussian with a centre frequency of 0.15 Hz and a standard deviation of 0.05 Hz. To reduce the computational cost we pre-compute the Green's functions G using the time-domain spectral-element codes for spherically symmetric Earth models AxiSEM to model wave propagation (Nissen-Meyer et al., 2014), and Instaseis to extract seismograms (van Driel et al., 2015) with 1-D isotropic PREM (Dziewonski & Anderson, 1981) as underlying velocity, attenuation, and density structure model. The Green's function database is then re-used for subsequent iterations and inversions. Additionally, we implement spatially variable grids with regional dense areas and no grid points on land, to reduce the number of possible noise sources, and thus the modeling parameters for regional applications (Igel et al., 2021). The combination of pre-computed wavefields and spatially variable grids allows us to efficiently invert for the noise source distribution of the secondary microseisms in a frequency range of 0.1-0.2 Hz (Ermert et al., 2020;Igel et al., 2021).

Inversion
Once we have modeled the synthetic cross-correlations, we measure the difference to the observed cross-correlations using the logarithmic energy ratio (e.g., Ermert et al., 2017). This measurement quantifies the asymmetry of the cross-correlation which arises from a heterogeneous noise source distribution as illustrated in Figure 1. The logarithmic energy ratio computes the ratio of the energies E + and E − of the expected surface wave arrival window w(τ) in the causal and acausal parts of the cross-correlation C(τ): In contrast to full-waveform misfits, the logarithmic energy ratio aims to match the energy in the causal and acausal expected surface wave arrival windows. For our frequency range we use a surface wave speed of 2,900 and an adaptive window size starting at 200 s. By linearly increasing our window size with inter-station distance we take surface wave dispersion into account and ensure the main wavelets are within the measurement windows. Although this measurement contains less information, it is much more robust, and relative to other measurements mostly insensitive to unknown 3-D Earth structure . This is because we only compare energies in the expected surface wave arrival time windows instead of the entire waveform. Independent of the frequency content of the cross-correlation, a time shift of the direct wave arrival introduced by a more complex velocity model does not have a significant effect on the measurement as long as the main wavelet is still within the window. Consequently, this allows us to use a simple 1-D PREM (Dziewonski & Anderson, 1981) velocity model to compute synthetic cross-correlations.
A disadvantage of the measurement is that we solely focus on direct wave arrivals and thus ignore any spurious arrivals. Different measurements, for example, a full waveform or envelope measurement, would take spurious arrivals into account but would not allow us to use 1-D velocity models. During the inversion we aim to minimize the squared L 2 -norm, that is, the misfit χ, of the measurements A i and on the synthetic and observed cross-correlations, respectively: where N is the number of measurements.
Adjoint techniques (e.g., Fichtner et al., 2006) allow us to compute source sensitivity kernels (e.g., Fichtner, 2014;Hanasoge, 2013b;Tromp et al., 2010) which provide a spatial reference of where an increase or decrease in noise source strength should decrease the misfit. By compiling the gradient, that is, the sum of all sensitivity kernels, we can update the noise source distribution and continue with the next iteration by re-computing the cross-correlations, misfits, and sensitivity kernels. To minimize the misfit we adopt a gradient-based iterative scheme using the steepest descent method, including regularization and step-length tests. Several synthetic and real-data tests have shown that there are usually no significant improvements in the noise source distribution after roughly five iterations independent of the initial model. Hence, we run eight iterations of the inversion to ensure that we have converged to a model that explains the data based on our measurement.
In previous research (Igel et al., 2021), we used a homogeneous distribution in the ocean as the initial noise source distribution. For a gradient-based iterative inversion method like ours, a good initial model can be helpful in steering the inversion toward an acceptable global noise source model and avoid local minima. By introducing a different method to locate noise sources-namely MFP-we are able to efficiently create a more realistic initial model from the same observed cross-correlations. This is similar to full-waveform inversions, where starting models are often constructed with more efficient methods such as ray-based travel time tomography or dispersion curve analysis (e.g., Teodor et al., 2021;Virieux & Operto, 2009).

Matched Field Processing
Matched Field Processing, in this context, may be considered similar to beamforming and backprojection methods, where time-shifts are applied to the data and rays are backprojected to obtain a source location. However, whereas beamforming generally assumes plane waves arrive at an array of sensors, MFP directly considers sources anywhere within a computational domain and estimates travel times accordingly. This makes it very suitable for global ambient noise source studies, where stations from all over the globe may be used. Additionally, it is able to map noise sources on any grid, which allows us to use the same source grid for MFP and the inversion.
MFP algorithms of varying complexity have been developed, for example,: to locate hydrothermal acoustic sources (Cros et al., 2011); microseismic sources in exploration geophysics (Corciulo et al., 2012); glacial tremors (Umlauft et al., 2021); or applied to three-component seismic array data for microseisms (Gal et al., 2018). The algorithm could also be adapted to be nearly identical to full-waveform methods by including synthetic Green's functions (Bowden et al., 2021;Schippkus & Hadziioannou, 2022). Although there may be some value to more complex MFP implementations, we prefer the computationally efficient version described below, as the subsequent inversion iterations will add further complexity.

Constant Velocity MFP
Our MFP algorithm is based on the assumption that a point source at a proposed noise source location for a set surface wave group velocity will lead to signal in the cross-correlation at a certain lag. The first step is to compute all cross-correlations and create a grid of possible noise sources. Subsequently, we iterate over all possible noise sources and compute the travel times to the stations based on a constant surface wave speed v of 2,900 . This surface wave speed is roughly the average Rayleigh wave group velocity in the 0.1-0.2 Hz frequency range in PREM (Dziewonski & Anderson, 1981) and has provided good results for synthetic and real-data applications. For such a narrow frequency band we can consider the group velocity to be roughly constant. The travel time difference Δt ij between arrivals t i and t j at the two receiver locations x i and x j determines the lag at which the current noise source location x would result in a signal in the cross-correlation: Note that ‖.‖ denotes the vector norm for a 2-D example as illustrated in Figure 2. For our applications we extend the vector norm from 2-D to a sphere where it is adapted to be the great circle distance between the noise source locations and stations.
Finally, the corresponding value of the cross-correlation-or in our case the value of the square envelope of the cross-correlation as explained in Section 2.2.2-is added to the "power" of that grid point, and we repeat the process for the next possible noise source location. This is equivalent to applying phase shifts to the raw signals and then measuring coherencies, as MFP or beamforming is often described. An illustration of this algorithm can be seen in Figure 2. By using a larger array of stations, we are able to spatially restrict the locations of noise sources and obtain a map of noise source "power."

Square Envelope Measurement
The simplest MFP method uses the value of the cross-correlation waveform to obtain the "power" for each noise source location. This is equivalent to applying a phase shift and measuring a zero-lag correlation coefficient. However, since this often results in strong fluctuations of the noise source power due to the oscillatory nature of the waveforms and struggles with low signal-to-noise ratios we instead take the value of the square envelope S(C(τ)) of the cross-correlation C(τ).
where  indicates the Hilbert transform. Additionally, we compute the standard deviation σ(S(C(τ))) of the square envelope and set all values below twice the standard deviation to 0, that is, we do not add any "power" for those noise source locations. We iterate over all possible noise sources, calculate the travel time difference Δt and finally obtain the value from the cross-correlation, for example, red dot for the actual waveform value or blue dot for the value of the square envelope. Additionally, we calculate the standard deviation of the envelope and set everything below twice the standard deviation to 0. This increases the signal-to-noise ratio and improves the final MFP power map. This process is repeated for every possible noise source location and every station pair.
Besides increasing the signal-to-noise ratio when a signal is present, this also smooths the resulting noise source and avoids the fluctuations of noise source power. This does mean that for cross-correlations with no clear signal, that is, where the square envelope is nearly constant, nothing is removed and the signal-to-noise ratio can not be increased. This has little effect on the final MFP power distribution as it would add a near constant value.
where P(C(τ i )) is the MFP power for the time lag τ at source location i and cross-correlation C(τ i ). Synthetic and real data tests show that using the square envelope with a cut-off threshold greatly increases the contrast of the final MFP maps and ensures that we mainly use the signal from the cross-correlations.
To account for geometric spreading we multiply each value of the square envelope of the cross-correlations with an amplitude decay factor D i as introduced by previous studies (e.g., Bowden et al., 2021;Corciulo et al., 2012), which depends on the surface wave group velocity v, the average frequency of our bandpass filter f, and the average distance of the station pair to the proposed source location r i : Hence, we ensure that the contribution of any potential source will be scaled according to the average distance to the station pair. This process is repeated for all possible noise source locations and cross-correlations, and the values of the square envelope of the cross-correlations are added up as illustrated in Figure 2. Of course, more sophisticated methods to model either the travel times or amplitude decays and attenuation can be implemented in MFP (Bowden et al., 2021;Schippkus & Hadziioannou, 2022). Such modeling is precisely the point of subsequent iterations of the full-waveform approach, whereas the MFP is only intended to give a computationally efficient initial model.
Whilst other array-based beamforming methods require a more local array to justify a plane wave assumption, MFP is advantageous for spatially large arrays that surround the dominant noise source location. We illustrate this in Figure 3 by running a synthetic example, where we forward model cross-correlations using the pre-computed Green's function database and cross-correlation model code from the inversion with a dominant noise source blob within the domain, and a frequency content of 0.1-0.2 Hz. We apply the MFP algorithm to two sets of stations: 35 stations in closer proximity and 35 stations spread out in the whole domain. If the dominant noise source is outside the array we see strong smearing and MFP is only able to give us a direction of the dominant noise source.
On the other hand, if the dominant noise source is surrounded by stations, MFP is able to constrain the spatial extent of the dominant noise sources more accurately.

MFP Starting Model
Thanks to MFP being most capable when the dominant noise sources are within the array, it is a useful method to locate noise sources on a regional to global scale. MFP and the finite-frequency inversion use slightly different information from the cross-correlation to obtain a noise source distribution. The logarithmic energy ratio is largely insensitive to unknown Earth structure, as it only takes the energy in a given window but ignores the actual waveform. On the other hand, MFP with the square envelope measurement uses more information from the waveform itself but does not properly account for wave propagation. Additionally, the resulting MFP maps are harder to interpret in terms of physical units as they are not an actual model of a physical quantity but rather an image of the noise source distribution.
To combine the two methods we normalize a smoothed MFP noise source map and set it as the initial PSD model for the finite-frequency inversion. In contrast to the previously used homogeneous starting model, this greatly reduces the presence of inversion artifacts. Synthetic and real-data tests have shown that this improves the final noise source maps without significantly increasing the computational cost. Figure 4 shows a regional synthetic comparison of two inversions with a homogeneous and an MFP starting model. The synthetic cross-correlations were modeled using the noise source distribution on the left with added random noise to make them more realistic. The random noise is introduced by normalizing a random time series, multiplying it by the maximum amplitude of the cross-correlation and a scaling factor of 1.5, applying a bandpass filter between 0.1 and 0.2 Hz, and finally adding it to the cross-correlation. Comparisons show that this resembles our real ambient noise cross-correlations more closely with a similar signal-to-noise ratio. Note that the initial coherent signal is still the strongest arrival after adding the noise as it is superimposed on the initial signal.
The inversion with the homogeneous starting model does contain the most dominant noise sources but shows a strong tendency to move noise sources closer to the coast, especially for the large dominant noise source area off the European coast. In contrast, the inversion with an MFP starting model does not lead to strong coastal sources, and better represents the spatial distribution of the dominant noise sources in the actual model. This is particularly useful for global inversions where MFP can help to avoid inversion artifacts due to lack of data by increasing the probability of noise sources in certain areas before the first iteration.
MFP introduces new information to the inversion, as it actually uses the cross-correlation waveforms, as opposed to the finite-frequency inversion where we measure the energy in the expected surface wave arrival time windows. Hence, we expect this to help avoid local minima by steering the nonlinear inversion in the right direction and thus produce a more accurate noise source map. Despite the clear differences in the resulting inversion models and a lower initial misfit for the MFP starting model, the misfits of the final iterations shown in Figure 4 are very similar. However, it is clear that the final model of the inversion with the MFP starting model is visually more similar to the target model than the inversion with a homogeneous starting model. This indicates that including the additional waveform information via the MFP starting model does help to steer the inversion in a direction that is more closely aligned with the actual noise source distribution.

SANS: Daily Seismic Ambient Noise Sources
In light of the recent developments that have significantly decreased the computational cost of ambient noise source inversions for the secondary microseisms on a global scale (Ermert et al., 2020;Igel et al., 2021), we introduce a new web framework, SANS, where daily seismic ambient noise sources are made available to the public (sans.ethz.ch). Currently, we run two inversions every day: one for a global station list and one more regional with stations surrounding the North Atlantic. A regional inversion allows for a higher spatial resolution of the noise source distribution in that area.
Users can obtain the inversion results by either directly looking at a plot of the noise source distribution maps online or downloading the inversion output and analyzing it themselves, for example, for implementation in other studies.

Data Selection and Processing
We download and process the seismic ambient noise data automatically every morning at 4 a.m. (CET) using ObsPy . All stations within a chosen station list are checked for available data. The station lists are based on globally available broadband sensors but limit the minimum distance between stations to 1°(≈111 km) to avoid smaller arrays. Dense station arrays would lead to high local singularities of sensitivity that would distort the final noise source distribution. This results in 414 stations for the global and 153 stations for the North Atlantic station list. The global distribution of stations is illustrated in the resolution analysis in Figure 6 and both station lists can be downloaded from the website. The data availability changes on a daily basis, leading to roughly 70% of these stations having data available on average.
After downloading all available data, we remove the instrument response, downsample to 1 Hz, segment the data into 2 hr windows, and remove any windows containing earthquakes that are in the GCMT catalogue (Ekström et al., 2012) with a minimum magnitude of 5.6. Occasionally this can lead to all windows being removed if there was one strong or several smaller earthquakes in a day. Subsequently, we compute the daily cross-correlations of the windowed seismic ambient noise data by stacking the individual cross-correlations of the 12 windows. This helps to increase the signal-to-noise ratio of the final daily ambient noise cross-correlations.
Similar to Igel et al. (2021), we ignore cross-correlations with a signal-to-noise ratio below 3.5 as this is just above the threshold of the average signal-to-noise ratio of a random time series. The signal-to-noise ratio is determined by dividing the maximum amplitude within the expected surface wave arrival window by the standard deviation of the whole time series. Hence, we define a signal as a clear surface wave arrival within the expected window. This is usually the case if the dominant noise source is in-line with the station pair. Due to our chosen measurement being the logarithmic energy ratio, cross-correlations with little signal-that is, no clear surface wave arrival-would not contribute much to the final gradient and update of the noise source distribution. Besides improving the final result of the inversion, ignoring cross-correlations below a signal-to-noise threshold also decreases the computational cost, as fewer cross-correlations have to be modelled during the inversion. During the inversion we apply a band-pass filter between 0.1 and 0.2 Hz as we focus on the secondary microseismic noise sources.

Web Framework
After the data have been downloaded, processed, and correlated we run eight iterations of the inversion on Piz Daint, a supercomputer at the Swiss National Supercomputing Centre (CSCS). The computational cost of the inversions varies with the number of available cross-correlations for each day. However, we greatly reduce the computational cost since we have already pre-computed the wavefield and extracted the Green's function database which is re-used every day. We run both inversions, one global and one regional surrounding the North Atlantic, on 600 cores, with the usual computational times being 60 min (50 node hr) and 30 min (25 node hr), respectively. We use two different spatially variable grids, with a more homogeneous distribution of about 29,000 grid points for the global inversion and a locally dense grid in the North Atlantic with roughly 21,000 grid points for the regional inversion. Once the inversions are done, we plot the output and copy all relevant files to the ETH web server where the website is hosted. These files are then made available to the public on sans.ethz.ch.
The web framework allows users to look through the iterations of all available inversion results and compare them to significant wave height maps (Tolman & Chalikov, 1996;WAVEWATCH III, 2005) of that day. Note that the generation mechanism of the secondary microseisms requires ocean waves traveling in opposite directions to overlap (Ardhuin, Hanafin, et al., 2011;Hasselmann, 1963;Longuet-Higgins, 1950;Nakata et al., 2019); therefore the wave height maps are merely a reference as to where the areas of dominant noise sources may be and should not be directly compared. Users can download the full inversion output folder including the parameter file, station list, source grid, further plots such as the gradients, misfit reduction and other relevant files. We provide code that helps a user to plot and analyze these results themselves. In theory, the inversions are reproducible as the inversion code is made available on GitHub (https://github.com/jigel/noisi_inv). However, this does require the additional computation or download of an AxiSEM wavefield and access to HPC facilities. The global inversion requires roughly 50 node hr which includes the data download, processing, and eight iterations of the inversion but excludes the extraction of a Green's function database.
On the website, we present four different plots directly to a user: (a) iterations of the SANS inversion, (b) significant wave height map, (c) station sensitivity of the inversion, and (d) misfit evolution plot as seen in Figure 5 for the 21 November 2022. The station sensitivity is computed by taking the sum of the absolute values of all sensitivity kernels before they are weighted by the logarithmic energy ratio measurement. The resulting map is smoothed with the same smoothing parameters as the final iteration of the inversion. All values below 1% of the  (Tolman & Chalikov, 1996), (c) station sensitivity of the inversion, and (d) misfit evolution plot.
normalized sensitivity are masked out to indicate that the station sensitivity is likely to be insufficient to infer information about the noise source distribution. The empirical threshold of 1% is based on the analysis of numerous synthetic inversions. Generally, the station sensitivity can be interpreted as a proxy to see in which regions we should be able to constrain the noise sources more accurately. Due to the varying availability of data, the station sensitivity changes on a daily basis. Generally, the sensitivity is much lower in the Southern Hemisphere due to the lack of stations. Since the sensitivity kernels attenuate as we move further away from a station pair, the highest sensitivity is closest to the stations along the coastlines.
We choose this measure of the station sensitivity since we solve a nonlinear inverse problem for which a proper quantification of resolution with reasonable computational cost is not feasible. Hence, we use the above mentioned station sensitivity to provide daily approximations and perform a more detailed resolution analysis for different station distributions in the following section. It should also be noted that this is the sensitivity for the inversion only. Due to the different measurement in the MFP, where we look at the square envelope of the entire waveform, we have a bigger area of sensitivity that is not just in-line with the station pair orientation due to the expected surface wave arrival window. Hence, if the noise source distribution from the MFP contains sources in areas with low station sensitivity during the inversion, the probability of those sources being changed during the subsequent iterations is low. That does not necessarily mean that those noise sources are incorrect; it mainly indicates that the inversion is not able to provide any noise source constraints in those areas so the main method of retrieval is MFP.

Resolution Analysis
Recent efforts have estimated the resolution and covariance of noise source full-waveform inversions by treating it as a linear problem and using singular value decomposition (Xu & Mikesell, 2022). However, due to our nonlinear measurement of the logarithmic energy ratio this is not applicable to our inversion method. To show the effect of the changing data and station availability on the resolution of the inversions, we forward model cross-correlations with added random noise to more closely resemble the signal-to-noise ratio of observed cross-correlations for 414 stations around the globe and perform inversions with different station lists in Figure 6. The noise source distribution that we use to forward-model the data is an adapted significant wave height map from the WaveWatch III model (Tolman & Chalikov, 1996;WAVEWATCH III, 2005). The inversions are run with the same parameters as the daily SANS inversions. We choose different station lists from daily inversions to give a realistic station distribution that would be used for real-data inversions.
The inversion with 414 stations shows the model that we are able to recover using all potential stations. Due to the much higher station density in the Northern Hemisphere the resolution is higher and we are able to recover the dominant noise sources more accurately than in the Southern Hemisphere. As we decrease the number of stations we can see how the recovered model changes, especially when the stations are predominantly in Europe (160 stations) or North America (143 stations). However, even the inversions with fewer stations still include the most dominant noise sources from the target model. In that sense, the daily global inversions should not necessarily be seen as the global noise source distribution for each day, but rather the noise source distribution that the given station list is able to observe. Generally, the resolution in the Southern Hemisphere is lower due to the lack of station coverage and the North Atlantic usually has the highest resolution since it is surrounded by stations in Europe and North America. Additionally, the spatial resolution of our maps is limited by the wavelength of the seismic waves we are studying and the inherent shape of the sensitivity kernels and corresponding Fresnel zones that increase in size as we move away from the stations. Thus, we expect smaller, more detailed noise source areas to be smoothed as visible in the inversions in Figure 6.

Example Applications
In the following section we present two example applications of the daily SANS maps. First, we take the average of the daily inversions for Northern Hemisphere summer and winter to study the seasonal variations of the secondary microseisms. Second, we model cross-correlations for different noise source distribution models to illustrate the effect of a changing noise source distribution on the cross-correlation waveforms. Additionally, the SANS maps could serve as a reference for other studies that analyze the secondary microseismic noise source distribution.

Seasonal Analysis
Secondary microseismic sources are generated when two oceans traveling in opposite directions overlap, which in turn creates a vertical pressure wave. This induces seismic waves at the ocean bottom. The strength of these sources is directly related to the wave height of the overlapping waves (Ardhuin, Hanafin, et al., 2011;Hasselmann, 1963;Longuet-Higgins, 1950;Nakata et al., 2019).
Due to the seasonal variations in significant wave height we would expect a similar pattern for the noise source distribution of the secondary microseisms, which has already been observed over a century ago (Banerji, 1925;Burbank, 1912;Klotz, 1910). This relationship has recently been studied for various different frequency bands of microseisms ranging from the hum to secondary microseismic sources Gualtieri et al., 2021;Landés et al., 2010;Nishida & Fukao, 2007;Rhie & Romanowicz, 2006;Stutzmann et al., 2012). Thus, we would expect similar patterns to emerge if we average the daily inversions generated by the SANS workflow.
We include 335 daily inversions from the 17 May 2021 to the 2 May 2022 in the analysis and choose to define Northern Hemisphere summer (21 April to 21 October) and winter (21 October to 21 April) based on the Icelandic first day of summer in 2022. The final iterations of all inversions within those two time ranges are averaged, resulting in 164 inversions for the summer and 171 for the winter months. Before averaging, we smooth the noise source model with a 4° Gaussian smoothing filter to avoid any artifacts from small changes in the inversion parameters during that time period.
Similarly, we average the significant wave height maps from the WaveWatch III model (Tolman & Chalikov, 1996) as a comparison. This should not be taken as a direct comparison but more as a reference of where the probability of more dominant noise sources is higher. The actual mechanism of generation of secondary microseismic sources is more complicated and requires more complex modeling (Ardhuin & Herbers, 2013;Ardhuin, Stutzmann, et al., 2011;Nakata et al., 2019). Figure 7 shows the comparison of the normalized average significant wave height with the normalized PSD of the average SANS inversions for Northern Hemisphere summer and winter.
The average inversions show clear seasonal variations that are in-line with our expectations. During the Northern Hemisphere summer the more dominant noise sources are in the Southern Hemisphere, specifically the South Pacific. As shown by previous studies on seasonal noise source variations (Gualtieri et al., 2021;Landés et al., 2010;Stutzmann et al., 2012), the Northern Hemisphere winter is dominated by noise sources in the North Atlantic. This supports the result of previous studies and shows that the SANS inversions are able to observe the spatio-temporal variations of the secondary microseismic sources on various timescales. Figure 7. Comparison of the normalized average significant wave height and normalized PSD of the SANS daily inversion results for Northern Hemisphere summer (21st April to 21st Oct) and winter (21st Oct to 21st April) using a global station distribution ( ). Northern Hemisphere summer is dominated by sources in the Southern Hemisphere and Northern Hemisphere winter is dominated by sources in the North Atlantic. Note that the significant wave height maps should only be seen as a rough reference where more dominant noise sources could be possible as the actual generation mechanism of secondary microseisms is more complex.

Cross-Correlation Modeling
A changing noise source distribution has a significant effect on the cross-correlations, particularly on a global scale. A common assumption is a homogeneous noise source distribution which, in theory, results in a symmetric cross-correlation. However, the noise source distribution is often strongly heterogeneous and changes constantly. We forward model cross-correlations using our modeling code from the inversion for three different noise source distributions to illustrate the changes: (a) homogeneous distribution everywhere, (b) homogeneous distribution in the ocean, and (c) the final SANS inversion model for 9 March 2022.
In Figure 8 we plot the cross-correlations for 6 station pairs and the three different models. The cross-correlations are filtered between 0.1 and 0.2 Hz. As the noise source distribution becomes more realistic, the changes to the cross-correlations become more and more significant. Especially for full-waveform ambient noise studies, the influence of a changing noise source distribution should not be ignored. In many cases the main arrival also shifts significantly which makes travel time picking more difficult. We encourage future ambient noise studies to consider including information about the noise source distribution.
The effect of a changing noise source distribution is also visible for higher periods as we see in Figure 9. Here, we compare synthetic cross-correlations modeled with a homogeneous distribution everywhere and the final SANS inversion for 9 March 2022 for different bandpass filters with a constant lower frequency of 0.017 Hz (60 s) and upper frequency ranging from 0.025 Hz (40 s) to 0.2 Hz (5 s). For lower frequencies the difference is less pronounced but still clear as we show in the zoomed in surface wave window panel for the second bandpass filter. As the noise source distribution becomes more realistic, the changes in the cross-correlations become more and more significant. This effect should not be ignored in ambient noise studies. Note that the noise source distribution for the homogeneous forward models are set to 1 everywhere, not 0 as indicated by the colourbar.
As we include higher frequencies, these changes become stronger and picking a potential travel time becomes difficult. This could bias travel time and other ambient noise tomography measurements and potentially leading to misinterpretations in terms of (temporal) changes in subsurface velocities. Further research has to be done to quantify the effect that different noise source distributions and velocity models have on cross-correlation waveforms.

Discussion
With the daily computation of SANS maps we are able to study the interaction between the atmosphere, the ocean, and the solid Earth in near real-time. The daily maps show the clear heterogeneous nature of secondary microseismic noise sources and their strong spatio-temporal variations. Since the generation mechanism of secondary microseismic sources is quite well understood (Ardhuin et al., 2015;Ardhuin & Herbers, 2013;Ardhuin, Stutzmann, et al., 2011;Longuet-Higgins, 1950), these variations can also be studied by computing ocean surface pressure maps (Ardhuin, Stutzmann, et al., 2011) using significant wave height and bathymetry Figure 9. Comparison of differently filtered modeled cross-correlations for one station pair (AV.IVE-US.EGMT) and two different noise source distributions: homogeneous everywhere and final SANS inversion model for 9 March 2022. The difference in waveforms is less pronounced but still clear for lower frequencies but gets significantly stronger when higher frequencies are included. Especially for higher frequencies, picking a potential travel time because very difficult. Thus, different noise source distributions could significantly influence travel time and other measurements used for ambient noise tomography. data. A comparison in our previous research (Igel et al., 2021) shows that these ocean surface pressure maps and finite-frequency inversions coincide quite well with similar areas of dominant noise sources present in both. Similarly, we observe correlation between the location of dominant noise sources in our daily SANS maps with respect to the higher amplitude significant wave heights. Due to the generation mechanism requiring two overlapping waves traveling in opposite directions, this comparison should only be considered as a rough reference of where there is a higher probability of dominant noise sources. In addition to secondary microseismic sources, our framework would also observe marine seismic events that produce a prominent signal within our frequency range, such as submarine landslides or smaller earthquakes that are not in the catalog. Currently, we would not be able to distinguish between such events and microseismic sources.
Particularly for ambient noise tomography and monitoring, knowledge of the noise source distribution is advantageous to circumvent common assumptions like the quasi-randomness of the noise field and equipartitioning of the wavefield. For these methods, daily maps could help reduce the misinterpretation of noise distribution changes as subsurface velocity changes . Inverted for both the noise source distribution and subsurface structure at the same time. However, this comes at an increased computational cost. By already having knowledge of the noise source distribution beforehand, we could reduce the complexity and computational cost of such full-waveform ambient noise tomography methods. As we illustrate in Figures 8 and 9, a heterogeneous noise source distribution can have a significant effect on the cross-correlations which should not be neglected, especially when higher frequencies are included. This makes potential travel time measurements much more difficult and could strongly influence near real-time SANS studies.
To make our inversion process as efficient as possible we use a simple 1-D PREM Earth model to simulate the Green's functions and cross-correlations. Despite  showing that our measurement of the logarithmic energy ratio is largely insensitivity to unknown 3-D Earth structure, this simplification could have an effect on the inversion. However, seismic studies within our frequency range of 0.1-0.2 Hz are generally considered less sensitive to small heterogeneities in the crust. Future studies might incorporate more complex Earth models (e.g., Fichtner et al., 2018) by pre-computing the Green's function database using a wavefield solver like Salvus (Afanasiev et al., 2019). This would also allow the implementation of a fluid ocean layer and 3-D structure, albeit at the cost of increased computation time.
Furthermore, since the availability of ambient noise data changes daily, the number of stations included in the daily inversions can fluctuate greatly. This has an effect on the spatial sensitivity of the inversion, as dominant noise sources cannot be resolved without data from surrounding stations. In combination with the lack of grid points on land due to our parameterization, this can lead to inversion artifacts in areas where we would not necessarily expect dominant noise sources; for example, in marginal seas like Hudson Bay or the Mediterranean Sea. This also happens when there is a lack of coherent signals in the cross-correlations which are then ignored due to our data selection based on the signal-to-noise ratio. To analyze the effect of the changing data and station availability on the resolution, we perform a synthetic test where we forward model cross-correlations with added noise and invert for them using different station lists. Generally, the inversions with fewer stations still include the most dominant noise sources, albeit with a lower spatial resolution.

Conclusions and Outlook
We present a new web framework SANS (sans.ethz.ch) where daily SANS maps for the secondary microseisms on a regional to global scale are made available to the public. Two methods are combined to improve the final noise source distribution: MFP and a gradient-based iterative finite-frequency inversion. The efficient data-driven MFP approach provides a starting model to steer the inversion in the right direction. Pre-computed wavefields and spatially variable grids have decreased the computational cost of the inversions, allowing us to run the inversions every morning for the previous days' data and presenting the results shortly after. Users are able to download the inversion results and we provide code to ease the implementation of the noise source distribution maps into other workflows. Comparisons to significant wave height maps do show that areas with high waves and strong dominant noise sources often coincide, which is in-line with the generation mechanism of secondary microseisms. Furthermore, we compute the averages of the noise sources maps for Northern Hemisphere summer and winter and compare them to the averages of the significant wave height maps. These show very similar areas of stronger activity which are in-line with other studies: Northern Hemisphere summer has more dominant sources in the Southern Hemisphere and Northern Hemisphere winter is dominated by noise sources in the North Atlantic.
We hope that making the noise source distribution data readily available to the public encourages new tomographic studies and methods exploiting seismic ambient noise vibrations. The accuracy of tomographic models could be improved by implementing knowledge of the noise sources. Specifically studies that make assumptions about a homogeneous or quasi-random noise source distribution would benefit and this may lead to more accurate velocity models. Studies that focus on time-dependent velocity changes in the subsurface often try to observe changes on the order of 1% or less (e.g., Delaney et al., 2017;Zhan et al., 2013). Particularly for such monitoring purposes, it is important to verify that these changes are not a result of a changing noise source distribution. The near real-time SANS maps we present here are a first step to providing near real-time information about the noise source distribution of secondary microseisms that can be easily implemented in other studies. Future applications could also make this approach feasible for more local studies like the near real-time monitoring of avalanches and landslides. More analysis needs to be done to properly quantify the effect a changing noise source distribution has on ambient noise tomography and monitoring measurements. We hope this work encourages others to consider the effect an ever-changing seismic ambient noise field has on their research.