Generation of secondary microseism Love waves: effects of bathymetry, 3-D structure and source seasonality

S U M M A R Y Secondary microseisms are ubiquitous ambient noise vibrations due to ocean activity, dominating worldwide seismographic records at seismic periods between 3 and 10 s. Their origin is a heterogeneous distribution of pressure fluctuations along the ocean surface. In spherically symmetric earth models, no Love surface waves are generated by such a distributed surface source. We present global-scale modelling of three-component secondary microseisms using a spectral-element method, which naturally accounts for a realistic distribution of surface sources, topography and bathymetry, and 3-D heterogeneity in Earth’s crust and mantle. Seismic Love waves emerge naturally once the system reaches steady state. The ergodic origin of Love waves allows us to model the horizontal components of secondary microseisms for the first time. Love waves mostly originate from the interaction of the seismic wavefield with heterogeneous Earth structure in which the mantle plays an important role despite the short periods involved. Bathymetry beneath the source region produces weak horizontal forces that are responsible for a weak and diffuse Love wavefield. The effect of bathymetric force splitting into radial and horizontal components is overall negligible when compared to the effect of 3-D heterogeneity. However, we observe small and well-focused Love-wave arrivals at seismographic stations in Europe due to force splitting at the steepest portion of the North Atlantic Ridge and the ocean–continent boundary. The location of the sources of Love waves is seasonal at periods shorter than about 7 s, while seasonality is lost at the longer periods. Sources of Rayleigh and Love waves from the same storm may be located very far away, indicating that energy equipartitioning might not hold in the secondary microseism period band.


INTRODUCTION
Wind-driven ocean waves are responsible for the generation of the background seismic wavefield-referred to as 'seismic ambient noise' (Nakata et al. 2019)-recorded continuously by seismometers worldwide at periods between about 3 and 300 s. The background seismic wavefield carries imprints of the atmosphere and the ocean, where the signal is generated, as well as of the solid Earth, where seismic waves propagate. Seismic ambient noise is a unique source of information about the energy exchange between these systems (e.g. Lay et al. 2009;Gualtieri et al. 2018), as well as about the structure of the Earth (e.g. Shapiro et al. 2005). Presentday ambient noise tomography reconstructs empirical Green's functions based on cross-correlation of noise records, assuming a diffuse and isotropic distribution of sources. However, the sources of the background seismic wavefield are not isotropically distributed (e.g. Hillers et al. 2012;Stutzmann et al. 2012), and therefore empirical Green's functions may differ from theoretical predictions (e.g. Tsai 2009;Tromp et al. 2010;Fichtner & Tsai 2019). Improving our understanding of the sources and the generation mechanisms of the background seismic wavefield is therefore crucial both for studying the energy transfer between the different earth systems and for overcoming limitations and errors in imaging Earth's interior.
The spectrum of the background seismic wavefield exhibits two main peaks at periods of about 14 and 7 s, called the primary and secondary microseisms, respectively. Microseisms are due to ocean gravity waves interacting with the seafloor (primary microseisms, e.g. Hasselmann 1963) or with each other (secondary microseisms, e.g. Longuet-Higgins 1950). At longer periods (about 50-300 s), ocean infragravity waves interacting with the seafloor are responsible for the seismic 'hum' (e.g. Rhie & Romanowicz 2006;Ardhuin et al. 2015). For a review, see Ardhuin et al. (2019).  Secondary microseisms are the most energetic background vibrations, showing the largest spectral amplitude worldwide. Observations of the link between microseisms and weather and ocean storms date back to the end of the 19th century (Bertelli 1872). The generation mechanism of secondary microseisms was first investigated by Longuet-Higgins (1950) and Hasselmann (1963) who showed that second-order ocean wave motion due to interactions between trains of ocean gravity waves of equal frequency and nearly opposite propagation directions produces a pressure at the surface of the ocean which is not attenuated with ocean depth. Ordinary ocean gravity waves (periods between 1 and 30 s) are driven by wind and have a wavelength of tens to hundreds of meters (Kinsman 1984). Therefore, their interaction occurs just below the ocean surface. This second-order pressure fluctuation generates compressional waves that propagate within the ocean and are transmitted and converted at the seafloor (Gualtieri et al. 2014).

Secondary microseism
As expected from sources located at Earth's surface, secondary microseisms are dominated by surface waves (e.g. Haubrich & McCamy 1969;Lacoss et al. 1969). Fundamental-mode Rayleigh waves dominate the vertical components of ambient noise records (e.g. Tanimoto et al. 2006), and their spectral amplitude alone can explain the noise spectrum of the vertical-component seismic record (e.g. Gualtieri et al. 2013). Rayleigh-wave overtones have been observed in a few cases (e.g. Kimman et al. 2012).
The generation mechanism of secondary microseisms as currently accepted (e.g. Miche 1944;Longuet-Higgins 1950;Hasselmann 1963) explains only vertical-component ambient-noise data (e.g. Gualtieri et al. 2013). The Rayleigh-wave content in vertical-component records is largely understood, but Rayleigh waves alone cannot explain the ambient-noise content on the horizontal components. The generation mechanism of secondary microseisms does not justify the presence of Love waves in the seismic noise records. Yet Love waves have been widely observed in the background seismic wavefield (e.g. Darbyshire 1954;Darbyshire& Iyer 1958;Toksöz & Lacoss 1968;R i n d&D o n n1979; Friedrich et al. 1998;Nishida et al. 2008;Tanimoto et al. 2015). Rind & Donn (1979) estimated the proportion of Rayleigh-to-Love-wave energy recorded at Palisades, NY, and found that this ratio was significantly different for varying directions of arrival. They hypothesized that a pressure source acting on an inclined surface-thus imparting a force with both vertical and horizontal components-could be the main mechanism for explaining the origin of Love waves. Each pressure source at the ocean surface generates a force perpendicular to the local seafloor. This force at the ocean bottom splits into forces that are vertical and horizontal. The horizontal forces can be responsible for the generation of Love waves. Rind & Donn (1979) also pointed out that 3-D heterogeneities along the source-receiver path could contribute significantly by modifying the amplitude of Love waves differently than the amplitude of Rayleigh waves. They also reported a summary of pioneering observations of Rayleigh-to-Love-wave energy ratios around the world. These early observations were carried out by modelling microseisms as harmonic motions and statistically estimating the correlation and coherence of the phase differences between the three components of motion. More Rayleigh than Love-wave energy was documented in regions such as Japan (e.g. Ikegami 1962) and California (e.g. Byerly & Wilson 1938;Gutenberg 1958;Haubrich & Iyer 1962), and more Love than Rayleigh-wave energy in Great Britain (e.g. Darbyshire 1954)and Montana (e.g. Haubrich & McCamy 1969).
Recent observations quantified the Rayleigh-to-Love energy ratio in the secondary microseism period band. Friedrich et al. (1998) used beam-forming analysis to locate microseism sources and assess the Love-to-Rayleigh energy ratio recorded at the Gräfenberg array (Germany). They analysed seismic waves in the period band of T = 6-11 s and found a Love-to-Rayleigh energy ratio (L/R)of 0.25. Nishida et al. (2008) applied a similar technique to Hi-Net data (Japan). They analysed data between T = 5-10 s and found L/R = 0.5-0.7. Tanimoto et al. (2015) used the rotational motion recorded by a ring laser and the displacement motion recorded by a collocated seismometer in Wettzell (Germany) to assess the amount of Love-and Rayleigh-wave energy, respectively. They found a Love-to-Rayleigh ratio of L/R = 1.2 at about T = 4.5 s. Tanimoto et al. (2016a) analysed the same data to look for seasonality patterns in the generation of Love waves and found a predominance of Love waves in winter (L/R = 1.2) and a predominance of Rayleigh waves in summer (L/R = 0.8). Using the same method applied to data at Piñon Flat (California), Tanimoto et al. (2016b) found L/R = 0.5 all year round. Juretzek & Hadziioannou (2016) confirmed the seasonal behaviour of the Love-to-Rayleigh ratio oscillating between L/R = 0.6 and L/R = 1.2. They also observed a frequencydependent variation of the Love-to-Rayleigh ratio. All these observations suggest a predominance of Rayleigh waves at stations worldwide (Friedrich et al. 1998;Nishida et al. 2008;Tanimoto et al. 2016b), with a few exceptions (Tanimoto et al. 2015), a clear seasonality in the amount of generated Love waves (Tanimoto et al. 2016a), and a frequency dependence (Juretzek & Hadziioannou 2016).
Numerical simulations successfully allowed for modelling the power spectral density (PSD) of secondary microseisms recorded on the vertical component. Kedar et al. (2008) presented the first quantitative modelling of secondary microseisms using ocean wave model hindcasts in the North Atlantic Ocean. Sources were modelled employing Longuet-Higgins' generation theory in deep water. Ardhuin et al. (2011) provided the first modelling that accounts for sources due to coastal reflections and presented the first global-scale maps of secondary microseism sources. They showed that the PSD of secondary microseisms can be modelled with great accuracy assuming great-circle approximation in a half-space earth model. Using this method, Stutzmann et al. (2012) further modelled the spectrum of secondary microseisms in a variety of environments and showed global-scale maps of sources during different seasons. Gualtieri et al. (2013) showed that the pressure field over extended areas due to storms is equivalent to a distribution of vertical point forces over a regular grid at the ocean surface and modelled the vertical component of the PSD using normal-mode summation. Longuet-Higgins' theory was extended to body waves (P and SV waves) by Gualtieri et al. (2014) and Ardhuin et al. (2013). Employing the generation theory developed by Gualtieri et al. (2014), Farra et al. (2016) modelled the amplitude of P waves resulting from beam-forming analysis. Gualtieri et al. (2013) also attempted to model the PSD recorded on the horizontal components by normal-mode summation, employing a 1-D earth structure. They found a gap between data and synthetics. Vertical forces in the absence of bathymetry and 3-D heterogeneities only generate Rayleigh waves and cannot generate Love waves. At the spectral peak T = 7s,L/R = 0.65 would have been sufficient to fill the gap between data and synthetics, in line with previous observations. Ziane & Hadziioannou (2019) investigated multiple scattering of surface waves as a possible mechanism for Love-wave generation. They modelled the wavefield generated by a single vertical point source in a highly heterogeneous 3-D half-space, altering model fluctuation strength, correlation length, and layering structure. They concluded that scattering effects alone are insufficient to obtain the observed L/R ratio from one single source.
Recently, Gualtieri et al. (2020) showed that, accounting for realistic bathymetry and 3-D structure as well as a realistic distribution of sources in the ocean, it is possible to discriminate between the roles played by bathymetry and structure in the generation of secondary microseism Love waves. Through scattering and focusing/defocusing effects generated by a distribution of sources, 3-D heterogeneities are able to generate an amount of Love waves in line with observations, while bathymetry alone yields a negligible amount of Love waves. Interference of the wavefields generated by the distribution of sources is crucial to explain the observed Love-to-Rayleigh ratio.
In this paper, we employ the spectral-element method (SEM, e.g. Komatitsch & Vilotte 1998), notably the SPECFEM3D GLOBE package (e.g. Komatitsch & Tromp 2002a,b), to extend the findings of Gualtieri et al. (2020) regarding the origin of secondary microseism Love waves, and to model three-component secondary microseism PSD at 199 stations worldwide. We simulate PSDs of secondary microseisms on the three components of seismographic stations worldwide, taking advantage of a recent technique developed in rotational seismology (e.g. Igel et al. 2005Igel et al. , 2007Hadziioannou et al. 2012) to estimate the amount of Love waves at each station. We compare results obtained employing two realistic configurations of secondary microseism sources, and two 3-D earth models. In Section 2, we discuss the source configuration, and in Section 3 the numerical implementation and simulation set-up. In Section 4, we compare our synthetic seismograms and PSDs with data. In Section 5, we discuss the generation of Love waves due to bathymetry, 3-D heterogeneities in isotropic and anisotropic earth models, and the seasonality in the generation of Love waves. We discuss and summarize our findings in Section 6.

SECONDARY MICROSEISM SOURCES
Analogously to earthquake seismology, we seek source time functions of secondary microseism sources to drive synthetic simulations. Sources of secondary microseisms are extended pressure sources at the ocean surface that can be discretized on a grid of equivalent point sources (Gualtieri et al. 2013).
Following Hasselmann (1963) and Ardhuin et al. (2011), the PSD of the pressure field (in units of Pa 2 m 2 s) at the surface of the ocean due to ocean wave-wave interaction can be written as where ρ w is the density of the ocean (assumed constant), g is the gravitational acceleration (assumed constant), f w is the frequency of ocean waves, f is the frequency of seismic waves (so that f = 2f w ), and k is the sum of the wavenumbers of the two ocean waves traveling in nearly opposite directions. For ocean waves moving in nearly opposite directions k ≃ 0,andF p (k, f ) only weakly depends on k.A l s o ,E(f w ) represents the PSD of the sea surface elevation (in units of m 2 Hz −1 ), which can be measured by buoys at the ocean surface, and I(f w ) is the non-dimensional ocean gravity wave energy distribution as a function of frequency, integrated over the ocean wave azimuth.
The quantity E 2 (f w )I(f w )i ne q .
( 1) can be estimated using numerical ocean wave model WAVEWATCH III (Tolman 2009), as improved by Ardhuin et al. (2010Ardhuin et al. ( , 2011 for coastal reflections and for wind-wave generation and dissipation. The model is specified globally with a spatial resolution of 0.5 • in latitude and longitude. At each grid point, the ocean state is described by 24 azimuths and 22 ocean-wave frequencies spaced between f w = 0.04 Hz (periods T w = 24.4 s) and f w = 0.30 Hz (T w = 3.29 s). The corresponding seismic periods range from T = 1.64 s to T = 12.2 s. WAVEWATCH III takes into account ocean-wave coastal reflection, empirically adjusting the amount of energy reflected from the coast for bathymetry and coastal shape (Ardhuin & Roland, 2012). Stutzmann et al. (2012) inverted for coastal reflection coefficients from observations by modelling the vertical component of secondary microseisms at stations in various environments. They found values between 4 and 10 per cent. This is in agreement with the 5 per cent global-scale average quoted by Longuet-Higgins (1950). In this work, we fix the coastal reflection to 5 per cent globally.
Following Hasselmann (1963), our eq. (1) can be written in the spatial (colatitude θ, longitude φ) and frequency domains as where dS = R 2 sin θ dθ dφ is the elementary surface, and R the radius of the Earth. The factor (2π ) 2 is needed to convert the pressure PSD from the wavenumber domain k ={k x , k y } to the spatial domain (Gualtieri et al. 2013).
In SPECFEM3D GLOBE, the ocean is incompressible and the entire ocean moves as the result of the normal displacement of the seafloor. The actual propagation of seismic waves in the ocean is neglected and the effect of the ocean on the seismic wavefield traveling in the solid Earth is taken into account by a load at the ocean floor. When the thickness of the ocean is small compared to the wavelength of the seismic waves (i.e. at periods T ≥ 20 s), this is a good approximation (Komatitsch & Tromp 2002b;Zhou et al. 2016). At shorter periods, the propagation of compressional (P) waves within the water layer has a non-negligible effect on the seismic wavefield and needs to be taken into account. In the case of secondary microseisms, sources are at the ocean surface and P waves generated at the sources are multiply reflected between the ocean surface and the seafloor. Instead of accounting for the ocean as a load at the seafloor, in our simulations the pressure sources of secondary microseisms are projected onto the seafloor by accounting for the reverberation of P waves in the ocean through an analytical correction. This 'source site effect' was first computed analytically by Longuet-Higgins (1950), and it varies with frequency, ocean depth, and seismic phase, for example, Rayleigh waves (Longuet-Higgins 1950;Gualtieri et al. 2013)orP and SV waves (Gualtieri et al. 2014). Love waves are not affected by the presence of the ocean (Komatitsch & Tromp 2002b). The source site effect facilitates moving the source from the top to the bottom of the ocean, to obtain an equivalent pressure source at the top of the crust (Gualtieri et al. 2013(Gualtieri et al. , 2014F arraet al. 2016).
Since the PSDs of secondary microseisms are dominated by surface waves, we multiply the time-varying pressure PSD in eq. (2)by frequency-and space-varying Rayleigh-wave source site effects, with c(f, θ, φ) as computed by Longuet-Higgins (1950)o v e r the fundamental mode and the first three overtones of Rayleigh waves. Fig. 1 shows the effect of the propagation of compressional waves in the ocean on the pressure PSD of secondary microseisms. Fig.  1(a) shows the pressure PSD at every source at the ocean surface given by the WAVEWATCH III model computed using eq. (2). In Fig. 1(b), we move the sources to the ocean-crust interface,  accounting for the propagation of compressional waves traveling in the ocean through the source site effect per eq. (3). Accounting for the compressibility of the ocean (Fig. 1b) dampens the sources, attenuating the overall amplitude of the pressure PSDs, with the strongest effect at periods shorter than about 6 s. At short seismic periods (T ≤ 6 s) compressibility needs to be taken into account for most ocean depths, while at long periods (T ≥ 6 s) the ocean can be considered incompressible except for deep basins. We note that the pressure PSD distribution (Fig. 1a) presents a splitting of the secondary microseism peak. As also observed by, for example, Figure 3. Pressure PSDs as a function of seismic wave period for (a) weak-amplitude and (b) strong-amplitude sources. The solid red circles denote the PSD at the frequencies specified by the ocean wave model. Amplitude and frequency content of these sources are very different. Bottom panels show the corresponding source time functions, as obtained from eq. (4). Meschede et al. (2017), the compressibility of the ocean enhances this feature (Fig. 1b), but the splitting is already present in the ocean wave-wave interaction sources (Fig. 1a). Fig. 2 shows the two actual source configurations employed in this work: pressure PSDs as computed by eq. (3), at three seismic wave periods: T = 5si nF i g s2(a) and (b), T = 7si nF i g s2(c) and (d) and T = 9si nF i g s2(e) and (f). The panels in the lefthand column show a typical source configuration during winter in the Northern Hemisphere, and the right-hand panels show a typical source configuration during winter in the Southern Hemisphere. We observe seasonality in the source configurations, with stronger amplitudes during the local winter: in the Northern Hemisphere in the left-hand panels and in the Southern Hemisphere in the right-hand panels. Seasonality is stronger at longer periods (T = 7sa n dT = 9 s). Each source is strongly frequency dependent.
We show the pressure PSD as a function of period for a weak source in Fig. 3(a), and for a strong source in Fig. 3(b). The PSD in Fig. 3(a) shows two peaks at about 5.6 and 7 s, while the PSD in Fig. 3(b) shows a nearly Gaussian shape, with a single peak around 7 s. Although the PSD of strong sources can be as much as two orders of magnitude above that of the weak ones, the number of weak sources is very high (as can be seen in Fig. 2), and therefore their contribution to the background noise level at any given seismographic station needs to be taken into account.
Ocean waves can be described as a sum of harmonics with energy at many frequencies. Since the phases ϕ i , i = 1, ..., N,o f these harmonic components generally are not correlated, we model them as independent random variables, uniformly distributed between 0 and 2π . The time-series of pressure (in Pa) associated with ocean wave-wave interaction at colatitude θ and longitude φ can be expressed in the form of a Fourier series, where N is the number of harmonics at seismic frequencies f i , f i is the frequency discretization of the ocean wave model, and t is time. The spectral amplitude is the square root of the PSD in eq. (3). The factor √ 2 ensures that the power of the pressure P(t, θ, φ) equals the integral of the PSD of the pressure F ′ p ( f,θ,φ)o v e ra l l frequencies.
Examples of source time functions computed using eq. (4)a r e shown in Figs 3(c) and (d), accompanying their power spectral densities shown in Figs 3(a) and (b). We applied a symmetric cosine taper to each end of the source time function, such that the signal is zero at the first and last data points and increases (decreases) smoothly at the beginning (end).

Theoretical framework
The spectral-element method (SEM) is based on the weak formulation of the equation of motion (e.g. Komatitsch & Vilotte 1998; Downloaded from https://academic.oup.com/gji/article/226/1/192/6166786 by guest on 20 April 2021 Igel 2017). In the absence of external forces and omitting complications due to rotation, attenuation and self-gravitation, which are incorporated but not relevant for the discussion here, instead of solving the momentum equation, where ρ is the spatial distribution of mass density, s denotes displacement, and T the stress tensor, with boundary conditions written in differential form, the SEM solves an integral of the momentum equation over Earth's volume dotted by a 'test' vector w: Stress and displacement are linked through Hooke's law T = c : ∇s, with c the fourth-order elastic tensor. Integrating by parts and using Gauss' theorem, the right-hand side of eq. (6) becomes where the first integral on the right-hand side of eq. (7)iso verthe surface of the Earth ∂ . The advantage of using the SEM for ambient noise forward modelling is that the source distribution can be taken into account naturally by modifying the surface boundary condition. Since the sources of ambient noise are pressure fluctuations, solving for the wavefield generated by those sources corresponds to modifying the customary stress-free (T ·n = 0) boundary condition at Earth's surface ∂ into where P is the time-varying pressure distribution, that is, the source time function in eq. (4). The unit vectorn accounts for the orientation of the forcing due to the bathymetry. In the SEM, the earth's volume is discretized with nonoverlapping elements (Komatitsch & Tromp 1999). The geometrical mapping in each element is based on Lagrange polynomials of degree n l defined at n l +1 Gauss-Lobatto-Legendre (GLL) points. We interpolate the pressure distribution in eq. (4) across all GLL points lying on Earth's surface bilinearly. The components of the source at each surface GLL point are determined by wheren i denotes a component of the unit outward normal vector in Cartesian coordinates, ω α , α = 0, 1, ...n l , are weights associated with the GLL quadrature and J (x(ξ α ,η β )) denotes the 2-D Jacobian. The total force applied on Earth's surface must be conserved while performing the interpolation over the GLL points: where n l is the number of GLL points on Earth's surface, x(ξ α ,η β ) represents the coordinates of the GLL points, P(t, θ , φ) is the original pressure PSD as given by WAVEWATCH III on each 0.5 • × 0.5 • patch (as in eq. 4)a n dP(x(ξ α ,η β ), t) is the interpolated pressure PSD distribution. The quantity ω α ω β J (x(ξ α ,η β )) represents the area of each surface element, and thus n l α,β=0 ω α ω β J (x(ξ α ,η β )) is equal to Earth's surface area. Fig. 4 shows that conservation of energy holds when interpolating the sources over all GLL points. The cumulative source time functions due to all sources before (blue dashed line, left-hand side of eq. 10) and after bilinear interpolation (red line, right-hand side of eq. 10)a r en e a r l y identical.

Simulation set-up
We u s e SPECFEM3D GLOBE (e.g. Komatitsch & Tromp 2002a,b) to perform three-component numerical simulations based on the theoretical framework described in the previous section. The earth model is elliptical in shape according to Clairaut's equation (Dahlen & Tromp 1998), self-gravitation is incorporated under the Cowling approximation (Dahlen & Tromp 1998) and the Coriolis effect is taken into account (Komatitsch & Tromp 2002b).
We perform global-scale simulations accurate down to T = 4s, with 960 spectral elements along each side of a chunk in the cubed sphere. For each simulation, the minimum resolution length is always smaller than the shortest wavelength. The total number of GLL points on Earth's surface is 230 400. Topography and 3-D heterogeneities are switched on and off to test their effects on the three components of the secondary microseism spectral amplitude.
In the absence of bathymetry, each pressure source is equivalent to a radial point force,n =r in eq. (9), applied to a spherical surface (Gualtieri et al. 2013). In a 1-D earth model (e.g. Dziewoński & Anderson 1981), a radial force generates P, SV and Rayleigh waves, but it does not generate any SH or Love waves (e.g. Dahlen & Tromp 1998). In the presence of bathymetry, each pressure source is decomposed into vertical and horizontal forces. The horizontal forces can generate SH and Love waves. The presence of 3-D heterogeneities can generate and enhance shear motion through focusing and defocusing of seismic waves, conversion of seismic waves, and scattering (e.g. Komatitsch et al. 2002).
We use a smoothed version of the ETOPO2 bathymetry and topography model (National Geophysical Data Center 2006), with a resolution of 4 min × 4 min (about 7.4 km × 7.4 km), which we call ETOPO4. This resolution captures the main bathymetric features, such as the continental slope and the structure of the ocean basins, without slowing down the computation or distorting the mesh excessively. The resolution of the topography and bathymetry is smaller than the minimum wavelength considered. Fig. 5 shows the unit vectors in the presence of bathymetry towards the east-west (Fig. 5a), north-south (Fig. 5b), and radial or downgoing (Fig. 5c). In the absence of bathymetry, the pressure sources are oriented vertically: the horizontal unit vectors are zero, and the radial unit vector is −1 everywhere over the ocean basins. The presence of bathymetry, and, notably, of the continental slope, generates non-zero horizontal forces (Figs 5a and b). As expected, north-south aligned bathymetric features, such as the North Atlantic ridge, mostly generate east-west oriented forces (Fig. 5a). The continental slopes around continents are responsible for the strongest horizontal forces, but many other non-zero horizontal forces are generated by abyssal topographic features. However, the tangential unit vectors are about two orders of magnitudes smaller than the radial ones. As observed by Rind & Donn (1979), ocean bottom slopes are typically gentle, and therefore the horizontal forces generated by force splitting at these slopes are small compared to the radially oriented force.
As we include continental topography in our simulations, in addition to bathymetry, a small amount of Rayleigh-to-Love converted energy may be due to continental features beneath each receiver. Cumulative source time function computed using the source discretization in WAVEWATCH III (dashed blue, left-hand side of eq. 10)a n d interpolated at all GLL points in SPECFEM3D GLOBE (red line, right-hand side of eq. 10). The inset at the bottom right shows a zoom between 1000 and 1100 s. The tangential vectors due to the topography (Supporting Information Fig. S1) have overall the same order of magnitude as the ones due to the bathymetry. However, since secondary microseisms are generated by a distribution of sources that act at the same time, the force splitting due to bathymetry contributes from all sources over the ocean at the same time, while a conversion in continental regions due to topography occurs only at the location of each receiver. Therefore, we expect bathymetry to play a more significant role than continental topography in the generation of Love waves.
We perform simulations in three different earth models, with and without bathymetry. As a 1-D earth model, we use the isotropic version of the PREM, the spherically symmetric Preliminary Reference Earth Model (Dziewoński & Anderson 1981). We consider two 3-D earth models: the isotropic model S40RTS (Ritsema et al. 2011)and the transversely isotropic model S362ANI (Kustowski et al. 2008). When a 3-D mantle model is chosen, the simulations are performed by incorporating 3-D crustal model Crust 2.0 (Bassin et al. 2000). This crustal model includes, in addition to three (upper, middle and lower) crustal layers, two layers of sediments, named 'soft' and 'hard' sediments. A summary of the simulations performed in this study is shown in Table 1.
In the secondary microseism period band, global-scale attenuation models are not easily constrained. Attenuation models at high frequencies have been developed only at the local scale (e.g. McNamara 2000). In the context of modelling secondary microseisms, the quality factor Q is usually treated as an unknown, and tuned in order to get the best fit between data and synthetics (Stutzmann et al. 2012;Gualtieri et al. 2013). We incorporate attenuation in our simulations using the PREM values. Attenuation is incorporated in SPECFEM3D GLOBE based on a superposition of standard linear solids (Komatitsch & Tromp 1999, 2002a. To reduce the computation time and to scale the problem towards high frequencies, we ran SPECFEM3D GLOBE on Graphics Processing Units (GPUs) available on the supercomputer Summit at Oak Ridge National Laboratory. The new portion of the code that allows for incorporating the distribution of pressure sources at the surface, and the corrections for bathymetry (see Section 3.1), was written in CUDA. Each forward simulation took about 8 hr on 600 GPUs.

Estimating the transverse component of microseisms
To quantify the amount of Love waves in the secondary microseism wavefield, it is necessary to rotate the east and north components of the seismogram into a transverse component, which is defined with respect to the azimuth of the seismic arrival-when it is known. Unlike typical earthquakes, sources of ambient noise are extended over the ocean floor, and the source azimuth dominating at a particular station needs to be derived from the data before the notion of a 'transverse' component applies.
The azimuth of the main arrival at a seismographic station can be estimated by measuring (or computing) the rotational motion, and notably the rotation rate around the vertical axis (Pancha et al. 2000;Igel et al. 2005Igel et al. , 2007. Only recently have rotational sensors enabled   (Dziewoński & Anderson 1981), S40RTS to the isotropic mantle model of Ritsema et al. (2011), S362ANI to the anisotropic mantle model of Kustowski et al. (2008). Crust 2.0 is the crustal model of Bassin et al. (2000) and ETOPO4 is a smoothed version of topographic model ETOPO2 (National Geophysical Data Center 2006 estimates of secondary microseism Love waves (e.g. Hadziioannou et al. 2012;Tanimoto et al. 2015Tanimoto et al. , 2016b. In particular, Hadziioannou et al. (2012) analysed rotational motion recorded with a ring laser at the Wettzell Geodetic Observatory in Germany. They found a consistently high correlation between rotation rate and transverse acceleration signals in the secondary microseism frequency band, and estimated the main direction of arrival as pointing toward a stormintheBayofBiscay .
Assuming plane wave propagation, the transverse acceleration a T is proportional toω z , the rotation rate around the vertical axis: where c is the phase velocity and ω z is the rotation around the vertical axis, defined as with s y and s x Cartesian components of the displacement vector s.
Since the transverse acceleration and the rotation rate around the vertical axis are in phase, during a seismic event their zero-lag cross-correlation coefficient approaches unity. Hence it can be used to estimate the transverse acceleration from the rotation rate.
To obtain an estimate of the transverse acceleration, we use SPECFEM3D GLOBE as a ring laser (Cochard et al. 2006) to compute the rotation rate. The estimated azimuth is very close to the theoretical azimuth computed from epicentral and station coordinates. The estimated transverse acceleration shows excellent agreement with the theoretical transverse acceleration. In the Appendix, we show that computing the rotation rate around the vertical axisω z with SPECFEM3D GLOBE, we can get a reliable estimate of the transverse acceleration a T owing to an earthquake, without knowing the event location. In Fig. A1, we observe that the cross-correlation coefficient betweenω z and a T computed over sliding time windows approaches unity during the event, most notably for Love waves. Following eq. (11), in Fig. A2, we also show that the transverse component of an earthquake can be estimated from the east and north components of the acceleration (a E , a N ) by seeking the azimuth φ = φ * that maximizes the zero-lag cross-correlation coefficient along the entire seismogram.

MODELLING SECONDARY MICROSEISMS
Following the framework described in the previous sections, we model three-component secondary microseisms at 199 stations worldwide, employing various earth models and source distributions (see Table 1).
In Fig. 6, we show seismograms at station HRV (Adam Dziewoński Observatory, Oak Ridge, Massachusetts, USA) of the Global Seismographic Network IU (IRIS/USGS). In black, we plot synthetics computed with SPECFEM3D GLOBE in 3-D earth model S40RTS with topography and bathymetry, for the Northern Hemisphere winter source configuration in Figs 2(a), (c) and (e) (simulation number 4 in Table 1). Figs 6(a)-(c) represent the displacement as computed on the east, north, and vertical components. Fig. 6(d) shows the rotation around the vertical axis, which is needed (after time differentiation, see eq. 11) to evaluate the azimuth of the main arrival and the transverse component, as detailed in Section 3.3. The grey seismograms in the background of the first three panels are the data during the same 3 hr. We observe that the amplitudes of the observed and synthetic seismograms are comparable only after about 2000 s (marked in red), the time needed for the system to reach steady state, as remarked by Gualtieri et al. (2020). This first portion of the synthetic seismograms is therefore removed before performing any analysis of the synthetic data.
In Figs 7 and 8, we show a comparison between observed (blue) and synthetic (magenta) PSDs at eight seismographic stations scattered around the North Atlantic Ocean, computed using mantle model S40RTS, topography and bathymetry. The shadow around each PSD represents the error associated with computing it, quoted as three standard deviations. To model horizontal components a realistic 3-D structure is needed: simplified models assuming a half-space (Stutzmann et al. 2012) or a 1-D earth model (Gualtieri et al. 2013) can only be used to model the vertical components of the data. The synthetic PSDs are in good agreement with data. Small deviations, of a few dB, between data and synthetic seismograms can derive from local variations in attenuation. It is worth recalling that the same radial attenuation model-PREM-is used to compute synthetic seismograms at all stations (see Section 3.2). Seismic waves in the secondary microseism period band are mostly affected by attenuation in the shallowest portion of the Earth. The shear quality factor at the surface of PREM, Q = 600, may be too high for some stations (Stutzmann et al. 2012) which could explain the small mismatch between observed and synthetic PSDs.

FACTORS AFFECTING THE E M E R G E N C EO FL O V EW A V E S
Without bathymetry and topography, pressure sources (Section 2) applied at the surface of a 1-D earth model, such as PREM, do not generate Love waves. However, Rayleigh waves are expected to have a small, but non-negligible, horizontal component (Dahlen & Tromp 1998) that affects the transverse component of seismic records. In the following, we apply the method described in Section 3.3 and the Appendix to estimate the azimuth of the main arrival, and the amount of Love waves for various earth structures and source configurations.

Effect of bathymetry
In the absence of bathymetry and 3-D heterogeneity, pressure sources do not generate SH (body) or Love (surface) waves. In Fig. 9, we show the cross-correlation between the rotation rate around the vertical axis and the transverse acceleration as a function of backazimuth and time, as computed in 1-D earth model PREM in the absence of bathymetry (simulation number 1 in Table 1). Each panel refers to a seismic station around the North Atlantic Ocean, where the maximum source power is located for this source configuration. Cross-correlation values (colour bar) oscillate between about ±0.2, with very few values exceeding 0.75 (black stars). The absence of large cross-correlation values and the absence of a clear maximum over time indicate that no Love waves are generated. Supporting Information Fig. S2 shows that the same conclusion can be drawn for stations located in the Southern Hemisphere.
The slope of the bathymetry beneath each source of secondary microseisms may be responsible for the generation of Love waves. Fig. 10 shows the cross-correlation between the rotation rate around the vertical axis and the transverse acceleration computed in 1-D earth model PREM in the presence of bathymetry (simulation number 2 in Table 1) for the same stations as in Fig. 9.A ss h o w n by Gualtieri et al. (2020), with this configuration it is possible to assess the effect of bathymetry on the generation of Love waves. Overall, the number of cross-correlation values exceeding a threshold of 0.75 is non-zero but very small. There is no clear indication of a preferred main arrival for stations located in the USA (top four panels in Fig. 10), while a more distinct arrival is identified at stations in Europe (bottom four panels in Fig. 10 Hemisphere do not show any clear indication of Love-wave arrivals (Supporting Information Fig. S3).

Effect of 3-D earth structure
In Fig. 11, the cross-correlation between the rotation rate around the vertical axis and the transverse acceleration is shown for the same stations as in Figs 9 and 10, but now considering 3-D earth model S40RTS in the absence of bathymetry (simulation number 3 in Table 1). Heterogeneous earth structure alone is able to generate Love waves in this set-up. All stations clearly show a distinct main arrival pointing toward the North Atlantic Ocean. Cross-correlation values exceeding 0.75 and approaching 1 are mostly aligned in narrow backazimuths. Compared to the case shown in Fig. 10,the increase of the number of cross-correlation values larger than 0.75 (black stars) is evident, a clear indication that Love waves are being generated. The same observation can be made at stations located in the Southern Hemisphere (Supporting Information Fig. S4). Fig. 12 shows the final scenario in which 3-D earth model S40RTS is used, while also taking into account bathymetry (simulation number 4 in Table 1). We observe small differences with respect to the previous scenario (which contained 3-D Earth structure but no bathymetry) in Fig. 11, and a similar main arrival at all stations. This confirms that bathymetry beneath the sources generates a negligible amount of Love waves relative to 3-D heterogeneity. Even at stations in Europe, where bathymetry alone was able to generate some Love waves, we observe a negligible difference (four bottom panels in Figs 11 and 12). At stations in the Southern Hemisphere (Supporting Information Fig. S5), cross-correlation values larger than 0.75 (black stars) occur at a wider range of backazimuths with respect to stations in the Northern Hemisphere (Fig. 12), indicating sensitivity to sources located at a wider range of backazimuths, and possibly a wider range of distances. We also observe small differences with respect to the previous scenario (Supporting Information Fig. S4), confirming that the amount of Love waves generated by bathymetry is negligible compared to the amount generated by the Earth's 3-D structure.

Bathymetry compared to 3-D structure
In Fig. 13, we show histograms of cross-correlations between the vertical rotation rate and the transverse acceleration exceeding the threshold of 0.75 as a function of azimuth for the four scenarios described above (black stars in Figs 9-12). The histograms reveal that bathymetry yields a small set of horizontals forces at the seafloor, which are responsible for a weak and diffuse Love wavefield (light blue histograms in Fig. 13). On the other hand, 3-D structure (black and red histograms in Fig. 13) produces a prominent main arrival at all stations. The maximum percentage of cross-correlation values exceeding the 0.75 threshold varies from station to station and does not depend on the distance from the coast (compare station WVT and station DWPF or HKT in Fig. 13). Stations in Europe show a small focused arrival also in the pres- histograms) and in Fig. 13], which may be due to locally prominent bathymetric features. In order to explore where this signal is coming from, we analyse the main directions of the Love-wave arrival at 10 seismic stations in Europe that show this small but clear arrival.
In Fig. 14 (e.g. Kolínskỳ et al. 2019), previous studies (e.g. Ardhuin et al. 2011;Stutzmann et al. 2012;Gimbert & Tsai 2015) showed that it is reliable to assume that most of the secondary microseismic energy is confined along the great-circle path. We observe several line crossings, mostly related to the North Atlantic Ridge or bathymetric features related to the ocean-continent boundary (e.g. see paths from stations ECH, BFO and GRFO crossing the Southern coast of Ireland). These features are likely responsible for the minor proportion of Love waves observed in Figs 10 and 13 (light blue histograms). With respect to the null case, where Love waves are not generated (blue histograms in Fig. 13), bathymetry is the only additional feature in this simulation and thus the only potential cause of Love-wave generation.
In Fig. 15, we show great-circle paths along Love-wave directions of arrival at 52 stations around the North Atlantic Ocean in the 3-D earth model S40RTS with bathymetry (simulation number 4 in Table 1) compared to the median pressure PSD (Fig. 15a) and to bathymetry (Fig. 15b). We observe that the region where the majority of great-circle paths cross (latitudes between about 45 • N and 60 • N, and longitudes between about −30 • Wand−50 • W) does not correspond to the maximum amplitude of the median sources (Fig. 15a). Moreover, it is located to the west of the steepest portion of the North Atlantic Ridge (Fig. 15b), as opposed to the case for which only bathymetry generated Love waves (Fig. 14).
Similar conclusions can be drawn for stations in the Southern Hemisphere (Supporting Information Fig. S6). Most of the stations around the equator and in the Southern tropical region (down to about 35 • S) are still sensitive to sources located in the Northern Hemisphere. Great-circle paths at other stations-most of which are located in the Southern extra-tropical regions-point towards the Southern Hemisphere. The few regions where great-circle paths cross in the Southern Hemisphere do not correspond to any maxima of the median source distribution (Supporting Information Fig. S6a) or maxima of the bathymetry (Supporting Information Fig. S6b).

Effect of mantle heterogeneity: does it matter?
Given the seismic periods involved, secondary microseisms are mostly sensitive to shallow structures, such as sedimentary and crustal layers. However, deep structure can have important effects in reshaping the entire wavefield, even at short periods. In order to investigate which portion of Earth's structure can be responsible for the generation of secondary microseism Love waves, and if deep structures play any role, we perform a simulation with the same source configuration used in the previous cases (Figs 2a, c and e) but with a different 3-D earth model-transversely isotropic earth model S362ANI (Kustowski et al. 2008)-in the presence of bathymetry. On top of mantle model S362ANI we superimpose crustal model Crust 2.0 (Bassin et al. 2000). This new simulation (number 5 in Table 1) differs from the previous one (simulation number 4 in Table 1) only in terms of mantle structure. In Fig. 16(a), we compare the cumulative per cent increase of cross-correlation values exceeding a threshold of 0.75 considering S362ANI with respect to S40RTS at 199 stations worldwide. The percentage increase is mostly positive in the Northern Hemisphere, where the sources have a larger amplitude (Figs 2a, c and e), and negative or negligible in the Southern Hemisphere. The transversely isotropic nature of the mantle enhances the generation of Love waves.
In addition, in the two simulations, Love waves originate from a different portion of the Earth. In Figs 16(b) and (c), we show the main direction of Love-wave arrivals at stations in the Northern Hemisphere, centred around the North Atlantic Ocean, and in the Southern Hemisphere, respectively. In both hemispheres, the main direction of arrival is different when considering a different mantle model (black and magenta arrows). For example, in Fig. 16(b), arrows point toward the North Atlantic Ocean in both cases, where the strongest source is located, but at most of the stations the arrows are oriented differently. This is an indication that, even if we expect crust and sediments to be the main source of secondary microseism Love waves, different heterogeneities within Earth's mantle also play a role in their origin.

Effect of source distribution
The generation of Love waves is mostly due to 3-D structure, with a significant contribution from the Earth's mantle in addition to crust and sediments, while bathymetry does not play a significant role. The source region of Love waves does not correspond to the location of the maximum source PSD, but it is mostly located in the Northern Hemisphere during its winter. In order to investigate seasonality in the generation of Love waves, as observed by Tanimoto et al. (2016a), we perform a simulation with a source configuration (simulation number 6 in Table 1) typical of winter in the Southern Hemisphere. Examples at three frequencies are shown in Figs 2(b), (d) and (f).
In Fig. 17, we show the cumulative per cent increase in crosscorrelation between vertical rotation rate and transverse acceleration exceeding 0.75 in the Southern Hemisphere winter scenario (simulation number 6 in Table 1) with respect to the same simulation using the source configuration typical of winter in the Northern Hemisphere (simulation number 4 in Table 1). The two simulations share the same mantle and crustal structures, the same bathymetry, and only differ in the source configuration. The cumulative percentage of cross-correlation values larger than 0.75 shows a north-south dichotomy, with a larger number of occurrences in the Northern Hemisphere during local winter (magenta in Fig. 17) and vice versa. As previously observed with displacement and rotation data recorded respectively with a seismometer and a ring laser (Tanimoto et al. 2016a), the generation of Love waves follows the seasonal pattern of the pressure sources. Gualtieri et al. (2020) showed that the Love-to-Rayleigh ratiodefined as the spectral ratio between the transverse and the vertical components-is realistic and in line with data and previous observations. In their fig. 4, they showed the Love-to-Rayleigh ratio at three seismic periods for 199 stations worldwide, using a Northern Hemisphere winter source configuration. In Fig. 18, we show the Loveto-Rayleigh ratio at the same 199 stations as Gualtieri et al. (2020) obtained using the Southern Hemisphere winter source distribution. The Love-to-Rayleigh ratio in Fig. 18 is computed at periods (a) T = 5s ,( b )T = 7sa n d( c )T = 9 s. As observed for the Northern Downloaded from https://academic.oup.com/gji/article/226/1/192/6166786 by guest on 20 April 2021 Figure 16. (a) Cumulative percentage increase of cross-correlation between vertical rotation rate and acceleration larger than 0.75 in mantle model S362ANI relative to S40RTS. (b) Main direction of Love-wave arrivals as computed in S40RTS (black) and S362ANI (magenta). In both simulations Crust 2.0 is taken into account; thus, the angle differences are due to differences in mantle structure. (c) Same as (b) but for the Southern Hemisphere.
Hemisphere winter-time configuration, in a Southern Hemisphere winter source configuration Rayleigh waves dominate at the majority of stations. The largest proportion of Love waves is recorded in the Southern Hemisphere (yellow-red colours in Fig. 18). However, we do not observe a clear north-south dichotomy. Fig. 19 shows the increase (red) and decrease (blue) of the Loveto-Rayleigh ratio considering a Northern Hemisphere winter-time source distribution (simulation number 4 in Table 1) with respect to a Southern Hemisphere winter source distribution (simulation number 6 in Table 1). Negative values represent an increase of the Love-to-Rayleigh ratio considering typical sources during winter in the Southern Hemisphere. At short periods (T = 5 s, Fig. 19a) the north-south dichotomy is very pronounced, with positive values in the Northern Hemisphere and negative values in the Southern Hemisphere. At long periods (T = 7-9 s, Figs 19b and c), the north-south dichotomy is less pronounced, although on average the Love-wave generation follows the seasonality of the sources, with a positive average value in the Northern Hemisphere (increase of Love waves with a Northern Hemisphere winter source distribution) and a negative average value in the Southern Hemisphere (increase of Love waves with a Southern Hemisphere winter source distribution). As observed for the much longer-period seismic 'hum' of the Earth (Tanimoto & Um 1999;Ekström 2001), seasonality persists in the average even for long-period secondary microseismic Love waves (T = 7-9 s) although it is less pronounced. At some stations, we observe a larger proportion of Love waves in the Southern Hemisphere during the local summer (red stations in the Southern Hemisphere in Figs 19b and c) and a larger proportion of Love waves at several stations in the Northern Hemisphere during the local summer (blue stations in the Northern Hemisphere in Figs 19ba n d19c). This behaviour is an indication that local structural effects can be strong at long periods in the secondary microseism period band. At these stations (red stations in the Southern Hemisphere in Figs 19bandc and blue stations in the Northern Hemisphere in Figs 19bandc),we observe that the role of the structure in the generation of Love waves is so strong that it overcomes source seasonality, generating Love waves very far away from the ultimate driving forces in the ocean. Such is not the case at short periods, where source distribution and seasonality are a direct proxy for Love-wave generation. Shortperiod seismic waves attenuate rapidly with distance and, therefore, Rayleigh waves can interact with heterogeneities only in the vicinity of their sources. On the other hand, long-period Rayleigh waves can travel long distances before being converted to Love waves.
Finally, an interesting clear north-south dichotomy is instead observed while looking at the direction of the main arrival of Love waves along the great-circle path. The main directions of arrival are quite stable through time in the Northern Hemisphere (Supporting Information Fig. S7a), while they show a pronounced seasonality in the Southern Hemisphere (Supporting Information Fig. S7b).

DISCUSSION AND CONCLUSION
In this paper, we show that employing a realistic 3-D earth modelincluding topography and bathymetry, ellipticity, gravity, rotation and attenuation-we are able to model all three components of secondary microseisms. Previous global-scale models (e.g. Ardhuin et al. 2011;Stutzmann et al. 2012;Gualtieri et al. 2013) only explained the vertical component of the data. To the best of our knowledge, this study presents the first modelling of threecomponent PSDs of secondary microseisms. This was possible through understanding and simulating the generation of Love waves (Gualtieri et al. 2020), which emerge ergodically due to the interaction of the seismic wavefield with 3-D wave speed heterogeneity.
As observed by Gualtieri et al. (2020), bathymetry beneath the sources has a negligible effect on the generation of Love waves. We observe a few exceptions for stations in Europe, where a significant proportion of Love waves is generated by splitting of the pressure force at the seafloor. Through triangulation of great-circle paths, we observe a potential source region corresponding to a steep portion of the North Atlantic Ridge. Other potential source regions, where great-circle paths cross each other, correspond to the ocean-continent boundary offshore the United Kingdom and Ireland. We did not take into account roughnesses of the seafloor smaller than about 7.4 km, but given that the seismic wavelengths involved here (about 12-30 km for waves traveling at 3 km s −1 at periods 4-10 s) are larger than the resolution of the bathymetry, we do not expect further effects due to higher resolution bathymetric features.
Love waves, and thus horizontal components, are strongly affected by 3-D heterogeneity. At the short periods analysed here, sediments and crustal layers play a major role. However, by simulating the generation of Love waves in different Earth structures and computing the horizontal components of secondary microseisms, we observe that the mantle affects the region where Love waves originate. The main source direction and the Love-to-Rayleigh ratio change according to the mantle model employed. This is evidence that deep structures reshape the entire wavefield even at short periods, and that the generation of secondary microseism Love waves cannot be understood by considering only the shallowest portions of the Earth.
The proportion of Love waves changes seasonally, as observed by Tanimoto et al. (2016a), and with frequency, as observed by Juretzek & Hadziioannou (2016). We also observe that changes in time and frequency are linked: at short periods (T 7 s) the location where Love waves originate follows source seasonality, while Downloaded from https://academic.oup.com/gji/article/226/1/192/6166786 by guest on 20 April 2021 Figure 18. Love-to-Rayleigh spectral energy ratio using a Southern Hemisphere winter source configuration at periods of (a) 5 s, (b) 7 s and (c) 9 s. The background colour represents the pressure PSD at the corresponding periods.
at long periods (T 7 s) the origin of Love waves is less tied to source seasonality. Love waves are generated in the vicinity of the dominant source region at short periods, while at long periods, sources can be as far away as the antipode. These observations have important consequences for the equipartitioning of seismic energy, which may still hold at short periods (T 7 s), as sources of Rayleigh and Love waves are relatively close in space, although at different depths. At long periods (T 7 s), we expect that energy equipartitioning no longer holds, given that secondary microseism Rayleigh and Love waves can be generated by sources thousands of kilometres away and located at different depths. This will affect tomographic modelling that uses three-component noise records as sources.
The framework described in this paper could be employed in future studies to better understand the generation of ambient noise in other period bands, such as primary microseisms or the seismic 'hum', or due to other disturbances, such as atmospheric pressure sources on continents. It is also suitable for being coupled Figure 19. Increase (red) and decrease (blue) of the Love-to-Rayleigh ratio considering a typical Northern Hemisphere winter source distribution with respect to a Southern Hemisphere winter source distribution at periods of (a) 5 s, (b) 7 s and (c) 9 s. with atmospheric or oceanic models to perform idealized studies on the coupling between the solid Earth and the other systems. In addition to being a first step towards ambient-noise tomography that accounts for a realistic, physics-based distribution of sources, findings in this paper are also relevant for the earthquake community, as realistic synthetic seismograms and spectra of ambient noise could be employed to denoise data and improve earthquake detection. obtained from the Data Management Center (DMC) of IRIS (Incorporated Research Institutions for Seismology), using the seismological community scientific library Obspy (e.g. Beyreuther et al. 2010;Krischeret al. 2015). SPECFEM3D GLOBE is freely available through the Computational Infrastructure for Geodynamics (CIG, https://geodynamics. org/cig/software/specf em3d globe/). This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC05-00OR22725. Additional computational resources were provided by the Princeton Institute for Computational Science & Engineering (PICSciE). We thank Tariq Alkhalifah and acknowledge funding from King Abdullah University of Science and Technology (KAUST) under grant number OSR-2016-CRG5-2970-01.

SUPPORTING INFORMATION
Supplementary data are available at GJ I online.    Fig. S2. Figure S4. Same as Fig. 11 in the main text, but for stations in the Southern Hemisphere. Clear Love-wave generation in a 3-D earth model without bathymetry. Cross-correlation between the rotation rate around the vertical axis and the transverse component as a function of backazimuth and time. Simulations are performed in 3-D earth model S40RTS without bathymetry. Layout as in Fig. S2. Figure S5. Same as Fig. 12 in the main manuscript, but for stations in the Southern Hemisphere. Clear Love-wave generation in a 3-D earth model with bathymetry. Cross-correlation between the rotation rate around the vertical axis and the transverse component as a function of backazimuth and time. Simulations are performed in 3-D earth model S40RTS with bathymetry. Layout as in Fig. S2. Figure S6. Segments of great-circle paths along the main direction of arrival for stations in the Southern Hemisphere centred around the South Atlantic Ocean, as computed in 3-D earth model S40RTS in the presence of bathymetry with a typical source distribution Downloaded from https://academic.oup.com/gji/article/226/1/192/6166786 by guest on 20 April 2021 during the Southern Hemisphere summer. The background colour represents (a) the median pressure PSD of the sources and (b) the slope of the bathymetry. Figure S7. Segments of great-circle paths along the main direction of arrival of Love waves at stations in the (a) Northern Atlantic Ocean and (b) Southern Atlantic Ocean, as computed in 3-D earth model S40RTS in the presence of bathymetry with a typical source distribution during the Northern Hemisphere winter (black) and Southern Hemisphere winter (magenta). Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

APPENDIX: ESTIMATING THE TRANSVERSE FROM ROTATION
We test the ability to estimate the proportion of Love waves in seismic records based on rotation using the 2004 February 7, M 7.3, Irian Jaya (Indonesia) earthquake as a case study. We perform simulations with SPECFEM3D GLOBE and we compute the three components of the displacement (s E , s N , s Z ) and the rotation rate around the vertical axisω z . We also compute the transverse acceleration based on the location of the epicentre and the station coordinates to make comparisons with the estimated transverse acceleration.
In Fig. A1(a), we show the transverse acceleration (blue) and the rotation rate (red) at station CAN (Canberra, Australia) of the Geoscope network. The two signals are in phase during the event and in particular for Love waves. The difference in amplitude is about four orders of magnitude, compatible with the order of magnitude of the phase speed of Love waves (eq. 11). In Fig. A1(b), we show the zero-lag cross-correlation coefficient between the two signals computed over 20 s time windows with 50 per cent overlap. To avoid maximum correlation for numerically zero values before the earthquake, Gaussian white noise has been added to the transverse acceleration. The maximum amplitude of the Gaussian random noise is 1 per cent of the maximum amplitude of the transverse acceleration. This step ensures zero correlation before the earthquake happens. The cross-correlation coefficient approaches 1 during the event and in particular for direct Love waves, while the coda is less well correlated, indicating the presence of both Rayleigh and Love waves. This method is therefore suitable for isolating the portion of the seismogram where Love waves dominate.
We can estimate the transverse acceleration by seeking the azimuth that maximizes the zero-lag cross-correlation coefficient over 20 s sliding windows. In Fig. A2(a), we show a comparison between the theoretical transverse acceleration (black) and the seismograms estimated by varying the azimuth from North in the range 0 • -360 • (colours). The cross-correlation coefficient between the rotation rate and these seismograms is shown in Fig. A2(b) (colour) as a function of azimuth and time. We observe that the cross-correlation coefficient is close to zero before the event, and, as observed in Fig. A1, the maximum of the correlation is found for Love waves. Taking the maximum value of the distribution of all azimuths for which the cross-correlation coefficient exceeds 0.9 (black crosses in Fig. A2b), we find that the estimated azimuth (grey dashed line) is very close to the theoretical azimuth (black dashed line). The difference is about 5 • , which is within the backazimuth ranges determined by Igel et al. (2007). The transverse acceleration computed with the estimated azimuth is nearly identical to the theoretical transverse acceleration (Fig. A2c). Figure A1. (a) Comparison between transverse acceleration (blue) and rotation rate around the vertical axis (red) for the M7.3 Irian Jaya earthquake in Indonesia (2004 February 7) at a station in Canberra, Australia. Seismograms were low-passed with a corner frequency of 7 s. The scaling factor betweenthese two quantities is approximately 3000 m s −1 , which is about the phase speed of fundamental-mode Love waves. (b) Zero-lag correlation coefficient between transverse acceleration and rotation rate compute over a 20 s sliding time window with 50 per cent overlap. One per cent Gaussian white noise has been added to the seismograms before computing the cross-correlation to avoid correlations equalling one for numerically zero signals. The correlation coefficient is stable around unity during the event.
Downloaded from https://academic.oup.com/gji/article/226/1/192/6166786 by guest on 20 April 2021 Figure A2. (a) Comparison between transverse acceleration as computed with SPECFEM3D GLOBE (black) and all possible transverse components computed by varying the azimuth from 0 • to 360 • . Seismograms were low-passed with a corner period of 7 s. (b) Zero-lag correlation coefficient (colour) between rotation rate and all possible transverse components as a function of azimuth (y-axis) and time (x-axis). One per cent Gaussian random noise has been added to the seismograms before computing the cross-correlation to avoid unit correlations for numerically zero signals. The black dashed line indicates the theoretical azimuth (159 • ), derived from epicentral and station coordinates. The grey dashed line denotes the estimated azimuth (154 • ), computed by finding the value that maximizes the zero-lag cross-correlation coefficient along the entire seismogram considering only values exceeding 0.9 (black crosses). The difference between the theoretical and estimated azimuths is relatively small and lies within the range determined by Igel et al. (2007). (c) Comparison between transverse acceleration as computed with SPECFEM3D GLOBE (red) and transverse acceleration computed with the estimated azimuth (blue   Northern-hemisphere winter source configuration Southern-hemisphere winter source configuration Northern-hemisphere winter source configuration Southern-hemisphere winter source configuration a) b) Figure S7. Segments of great circle paths along the main direction of arrival of Love waves at stations in the (a) Northern Atlantic Ocean and (b) Southern Atlantic Ocean, as computed in 3-D earth model S40RTS in the presence of bathymetry with a typical source distribution during the Northern Hemisphere winter (black) and Southern Hemisphere winter (magenta).