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

HAL
open science

# Ray-theoretical modeling of secondary microseism P waves

V. Farra, E. Stutzmann, L. Gualtieri, M. Schimmel, F. Ardhuin

► To cite this version:

V. Farra, E. Stutzmann, L. Gualtieri, M. Schimmel, F. Ardhuin. Ray-theoretical modeling of secondary microseism P waves. Geophysical Journal International, 2016, 206, pp.1730-1739. 10.1093/gji/ggw242. insu-03581297

HAL Id: insu-03581297

https://insu.hal.science/insu-03581297

Submitted on 19 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Ray-theoretical modeling of secondary microseism *P* waves

V. Farra,<fnref n="1" /> E. Stutzmann,<fnref n="1" /> L. Gualtieri,<fnref n="2" /> M. Schimmel<fnref n="3" /> and F. Ardhuin<fnref n="4" />

Accepted 2016 June 23. Received 2016 June 13; in original form 2015 November 4

## SUMMARY

Secondary microseism sources are pressure fluctuations close to the ocean surface. They generate acoustic *P* waves that propagate in water down to the ocean bottom where they are partly reflected and partly transmitted into the crust to continue their propagation through the Earth. We present the theory for computing the displacement power spectral density of secondary microseism *P* waves recorded by receivers in the far field. In the frequency domain, the *P*-wave displacement can be modeled as the product of (1) the pressure source, (2) the source site effect that accounts for the constructive interference of multiply reflected *P* waves in the ocean, (3) the propagation from the ocean bottom to the stations and (4) the receiver site effect. Secondary microseism *P* waves have weak amplitudes, but they can be investigated by beamforming analysis. We validate our approach by analysing the seismic signals generated by typhoon Ioke (2006) and recorded by the Southern California Seismic Network. Backprojecting the beam onto the ocean surface enables to follow the source motion. The observed beam centroid is in the vicinity of the pressure source derived from the ocean wave model WAVEWATCH III<sup>R</sup>. The pressure source is then used for modeling the beam and a good agreement is obtained between measured and modeled beam amplitude variation over time. This modeling approach can be used to invert *P*-wave noise data and retrieve the source intensity and lateral extent.

**Key words:** Seismic interferometry; Body waves; Seismic noise; Wave propagation.

## 1 INTRODUCTION

Seismic noise is the continuous oscillation of the ground recorded by seismometers worldwide (Gutenberg 1936; Webb 1998; Stutzmann *et al.* 2000, 2009; Berger 2004). Ocean waves generate most of the noise energy in the period band 1–300 s but depending on the period, different source mechanisms and/or ocean waves are involved (Longuet-Higgins 1950; Hasselmann 1963; Ardhuin *et al.* 2015). Seismic noise can be divided in three period bands and it is called secondary microseisms, primary microseisms and hum in the period bands 1–10, 10–20 and 20–300 s, respectively. In the period band 1–10 s, the secondary microseism sources are pressure fluctuations close to the ocean surface. These pressure sources result from the so-called secondary mechanism, that is non-linear interactions of ocean gravity waves having similar frequencies and coming from nearly opposite directions (Longuet-Higgins 1950; Hasselmann 1963; Ardhuin *et al.* 2011). At longer periods (10–20 s), primary microseism sources are pressure fluctuations located only close to the coasts. These pressure fluctuations occur at the ocean bottom and result from the so-called primary mechanism, that is, the direct coupling between ocean gravity waves and the bathymetry in shallow water (Hasselmann 1963; Ardhuin *et al.* 2015). At even longer periods (20–300 s), the hum source location and mechanism

have long been debated (e.g. Fukao *et al.* 2002; Rhie & Romanowicz 2004; Tanimoto 2007; Webb 2008; Nishida 2013; Traer & Gerstoft 2014). Recently, Ardhuin *et al.* (2015) showed that the primary mechanism applied to infragravity waves is presumably the source mechanism of Rayleigh waves in the hum period band. In that case, a pressure source at the ocean bottom is generated by the interaction of infragravity waves with the continental shelf. Long-period Love waves have also been detected (e.g. Kurrle & Widmer-Schnidrig 2008; Nishida *et al.* 2008) but they cannot be generated by pressure sources and complementary source mechanisms have to be considered. Fukao *et al.* (2010) and Nishida (2013) proposed that Love wave sources are random shear tractions at the ocean bottom due to the non-linear topographic coupling of infragravity waves.

In this paper, we focus on noise *P* waves in the period band 1–10 s, that is, on noise body waves generated by ocean gravity waves interacting through the secondary mechanism. The corresponding ocean wave phenomena can be classified into three sea-state configurations (Ardhuin *et al.* 2011). The first configuration (class I) occurs when a storm has a wide angular distribution, with ocean gravity waves coming from many different azimuths. In this case, the interacting waves are in the vicinity and/or inside the storm. For what concerns the second sea-state configuration (class II), ocean gravity waves reaching the coast are reflected back and interfere with

<footnote n="1">¹Institut de Physique du Globe de Paris, PRES Sorbonne Paris Cité, CNRS UMR 7154, 1 rue Jussieu, 75005 Paris, France. E-mail: stutz@ipgp.fr</footnote>

<footnote n="2">²Lamont-Doherty Earth Observatory, Columbia University, 61 Route 9W, Palisades, NY 10964, USA</footnote>

<footnote n="3">³Institute of Earth Sciences Jaume Almera, CSIC, Lluis Sole i Sabaris, s/n, E-08028 Barcelona, Spain</footnote>

<footnote n="4">⁴Ifremer, Univ. Brest, CNRS, IRD, Laboratoire d'Océanographie Physique et Spatiale (LOPS), IUEM, F-29200 Brest, France</footnote>
the upcoming ocean gravity waves. The corresponding interaction area is confined close to the coast. The third sea-state configuration (class III) relates to interactions of ocean gravity waves coming from different storms.

Whatever source class, ocean gravity wave interactions can be approximated by pressure fluctuations close to the ocean surface (Hasselmann 1963; Gualtieri et al. 2015). These pressure fluctuations generate acoustic P waves which propagate in the ocean. When the P waves reach the ocean bottom, they are partly reflected and partly transmitted into the crust and they continue their propagation through the Earth (Ardhuin & Herbers 2013; Gualtieri et al. 2014). Both body waves and surface waves are generated at the ocean bottom but body waves have a much smaller amplitude than surface waves and cannot be directly observed on seismograms (Stutzmann et al. 2012; Gualtieri et al. 2013).

Noise body wave characteristics can be extracted by stacking signals recorded by an array of stations. Beamforming techniques have been extensively used for determining noise wave types and source location. The first noise P-wave source detections were reported more than 50 yr ago (Haubrich et al. 1963; Lacoss et al. 1969). Vinnik (1973) used secondary microseism P waves recorded in central Russia to estimate the force source associated with medium power cyclones. More recently, several studies demonstrated that a significant amount of P-wave microseism energy is generated far from the coast in deep oceans (e.g. Gerstoft et al. 2008; Koper et al. 2009, 2010; Landès et al. 2010; Obrebski et al. 2013). Sources of body waves have also been associated with specific storms (e.g. Schulte-Pelkum et al. 2004; Gerstoft et al. 2006; Koper & de Foy 2008; Zhang et al. 2010a,b; Davy et al. 2014). Gal et al. (2015) further showed that P waves recorded in Australia are mostly generated close to the coast at periods shorter than 2 s and in deep water at longer periods.

In this paper, we focus on the amplitude of P waves in the secondary microseism period band. In Section 2, we present the theory for computing the power spectral density (PSD) of the P-wave displacement recorded by an array of seismic stations and due to secondary microseism sources. In order to validate our computation, we analyse in Section 3 the well-documented typhoon Ioke which occurred in 2006 (Zhang et al. 2010a,b; Gualtieri et al. 2014) and generated seismic signals recorded by the California network (CI) over several days. We compute the beam per day and backproject the beam centroid onto the ocean surface in order to make a comparison with the pressure PSD model derived from the ocean wave model WAVEWATCH III<sup>R</sup> (Tolman 2009; Ardhuin et al. 2011; Rascle & Ardhuin 2013). We follow the P-wave source over several days and investigate the related intensity variations. We compare observed and modeled beam PSD amplitudes and lateral extent. The good agreement between them validates our modeling approach.

## 2 MODELING NOISE P-WAVE AMPLITUDE

In this section, we present the modeling of P waves recorded in the far field and generated by the secondary mechanism, that is by ocean wave interactions close to the ocean surface (Longuet-Higgins 1950; Hasselmann 1963). The P-wave amplitude recorded at receivers in the far field depends on: (1) the source magnitude, (2) the source site effect that takes into account the propagation within the water layer, (3) the propagation from the ocean bottom to the receiver and (4) the receiver site effect.

In this paper, we use the Fourier transform sign convention used by engineers for signal processing. Therefore, the Fourier transform of the function p(t) is given by:

$$P(f) = \int_{-\infty}^{\infty} p(t)e^{-2i\pi ft} dt.$$

### 2.1 Secondary microseism sources

Secondary microseism sources in the period band 1–10 s are due to non-linear interactions between ocean gravity waves and can be represented as a random pressure field, P(x, f), acting at the ocean surface (Hasselmann 1963). The pressure field varies with the seismic frequency f, which is twice the frequency of the interacting ocean waves. The PSD of this random pressure field, $F_p(x, f)$ (Pa$^2$ m$^2$ s), as a function of the 2-D spatial coordinate vector x and frequency f, is given by (Hasselmann 1963; Ardhuin et al. 2011):

$$F_p(x, f) = [2\pi]^2 [\rho_w g]^2 f E^2(f/2) \times \int_0^{\pi} M(f/2, \phi) M(f/2, \phi + \pi) d\phi, \quad (1)$$

where $\rho_w$ is the water density, g is the gravity acceleration, f/2 is the ocean wave frequency, E(f/2) is the sea surface elevation PSD (in m$^2$ s) and M(f/2, φ) is the non-dimensional ocean gravity-wave energy distribution as a function of ocean wave frequency f/2 and azimuth φ.

Eq. (1) is derived from eq. (2.15) in Hasselmann (1963), which gives the pressure PSD $F_p(K, \omega)$ as a function of angular frequency $\omega = 2\pi f$ and horizontal wavenumber K. K is the sum of the wavenumbers of two ocean gravity waves with frequency f/2 and traveling in opposite directions. Using statistical hypotheses about the ocean wave spectrum, Hasselmann (1963) showed that the pressure PSD $F_p(K, \omega)$ is independent of the horizontal wavenumber K; therefore, the random pressure field P(x, f) is white in the space domain, which means that pressure sources at two different locations are independent and

$$\langle \bar{P}(\mathbf{x}, f) P(\mathbf{x}', f) \rangle = F_p(\mathbf{x}, f) \delta(\mathbf{x}' - \mathbf{x}), \quad (2)$$

where cornered brackets denote ensemble mean, x and x' are 2-D position vectors and $\bar{P}(\mathbf{x}, f)$ denotes the complex conjugate of P(x, f).

### 2.2 Source site effect

Pressure sources at the ocean surface generate only P waves that propagate down to the ocean bottom. In the ocean layer, upgoing P waves are generated by reflection at the seafloor and downgoing P waves by reflection at the free surface. Gualtieri et al. (2014) computed the source site effect as the constructive interference of multiply reflected plane P waves in the ocean that are converted to either P or SV plane waves at the ocean–crust interface. They showed that the ocean site effect on the transmitted plane P-wave potential can be written as:

$$C_P(i_w, h, \omega) = \frac{T_P(i_w)}{1 + R_P(i_w) e^{-i\phi_w(h, \omega, i_w)}}, \quad (3)$$

where $i_w$ is the plane P-wave take-off angle in the ocean, h is the water depth at the source location and $\omega$ is the seismic angular frequency. $T_P(i_w)$ and $R_P(i_w)$ are the P-wave potential transmission and reflection coefficients at the seafloor, respectively. The phase
<figure><img src="image_3.png" /><figcaption>Figure 1. Map of the source site effect modulus $|C_p|$ at a period of 4 s (left), 5 s (middle) and 6 s (right). The black line corresponds to the typhoon Ioke track estimated as the maximum significant wave height. Black dots indicate time every 12 hr from 2006 September 2 00:00 to 6 00:00. The white area on the lower left corner corresponds to epicentral distances larger than 90° with respect to the CI seismic network.</figcaption></figure>

shift $\phi_w(h, \omega, i_w)$ is due to propagation within the water layer and is defined as:

$$ \phi_w(h, \omega, i_w) = 2\omega \frac{\cos(i_w)}{\alpha_w} h, \quad (4) $$

where $\alpha_w$ is the water P-wave velocity. Let us note that, because we do not use the same Fourier transform sign convention as in Gualtieri et al. (2014), we have changed the sign in front of the phase shift in eq. (3).

Gualtieri et al. (2014) showed that the source site effect strongly depends on frequency and ocean depth. The same result was derived with a local mode formalism by Ardhuin & Herbers (2013). The ocean site effect on P waves also depends on the take-off angle, that is, on the distance between the source and the receiver. Fig. 1 shows the modulus of the source site effect $|C_p|$ for sources in the Western Pacific Ocean and receivers in California. The source area is in the vicinity of the typhoon Ioke track (2006 September 2–5) studied in Section 3. The typhoon track (black line) is defined as the maxima of the significant wave height over time. In the vicinity of the typhoon track, the maximum source site effect is at a period of 5 s, at which it varies between 0.9 and 3.4. At shorter (4 s) and longer (6 s) periods, the modulus of the source site effect is mostly lower than 1.

## 2.3 P-wave vertical displacement in the far field

The pressure sources at seismic frequency $f$ due to non-linear interactions among ocean gravity waves at frequency $f_g = f/2$ are located within a layer at the top of the ocean. The thickness of this layer is comparable to the wavelength $\lambda_g$ of the ocean gravity waves, which satisfies the deep water dispersion relation $(2\pi f_g)^2 = 2\pi g/\lambda_g$. At the periods we are interested in, the layer thickness is much smaller than the seismic wavelength. For example, at the seismic period of 5 s, the layer thickness is around 160 m to be compared to the seismic wavelength $\lambda_s$ in water which is 7.5 km.

By using the elastodynamic representation theorem (Aki & Richards (2002), eq. 2.43), one can write the displacement in the far field as the sum of two terms: one is a volume integral over pressure sources inside the ocean and the other one is a surface integral over the vertical tractions induced by pressure fluctuations on the ocean surface. Using the Green's function satisfying the free-surface condition on the ocean surface in the elastodynamic representation theorem, one can verify that the volume integral term is of the order of $(\lambda_g/\lambda_s)^2$ smaller than the surface integral term (see also Hasselmann (1963), section 2). As only the surface integral term

contributes to the expression of the displacement given by the representation theorem, we compute the displacement in the far field only due to vertical tractions acting on the ocean surface.

Let us consider a single pressure source $P(\mathbf{x}_s, f)$ acting on a small surface $S$ centred at $\mathbf{x}_s$ on the ocean surface. By using the ray theory approach (see Appendix), one can show that the P-wave vertical displacement at a receiver $\mathbf{x}_r$ located on the free surface in the far field is given by:

$$ u_Z(\mathbf{x}_r, f) = [2 \cos(i_r)] \left[ \frac{\cos(i_c)}{4\pi \rho_c \alpha_c^2} \frac{e^{-\pi f i_r^*}}{\sqrt{J_P^E}} \hat{T}_P^E e^{-2i\pi f t_P(\mathbf{x}_s, \mathbf{x}_r)} \right] \times \left[ 2 C_P \frac{\rho_c}{\rho_w} \right] [P(\mathbf{x}_s, f) S], \quad (5) $$

where $\alpha_c$ is the crust P-wave velocity, $\rho_c$ and $\rho_w$ are the crust and water densities.

Let us describe eq. (5). For P-wave modeling in the far field, multiplying the pressure source with the source site effect, $[2 C_P(i_w, h, \omega) \frac{\rho_c}{\rho_w}] [P(\mathbf{x}_s, f) S]$, is equivalent to consider a vertical force source at the top of the crust layer and propagation of the generated P waves through an Earth model without ocean. The factor 2 is introduced to take into account that the pressure source is at the ocean surface.

Considering this vertical force source at the top of the crust layer, seismic waves propagate through the solid Earth to the receiver. In a homogeneous medium, for example a homogeneous crust, the amplitude of the P-wave displacement due to a unit vertical force at a distance $r$ is given by the expression (Aki & Richards 2002):

$$ \frac{\cos(i_c)}{4\pi \rho_c \alpha_c^2 r}, $$

where $i_c$ is the take-off angle of the ray in the crust. Considering the propagation in a spherically symmetric Earth model, the geometrical term $\frac{1}{r}$ is replaced in eq. (5) by $\frac{\exp(-\pi f t_P^*)}{\sqrt{J_P^E}} \hat{T}_P^E$, where $J_P^E$ is the geometrical spreading given by eq. (A4), and $t_P^*$ is the seismic attenuation term along the P-wave ray path. The coefficient $\hat{T}_P^E$ is the product of the normalized P-wave transmission coefficients corresponding to all the discontinuities encountered by the P-wave ray along its trajectory inside the Earth. In the phase term $e^{-2\pi i f t_P(\mathbf{x}_s, \mathbf{x}_r)}$ of eq. (5), $t_P(\mathbf{x}_s, \mathbf{x}_r)$ is the P-wave traveltime between the source at $\mathbf{x}_s$ and the receiver at $\mathbf{x}_r$.

The factor $[2\cos(i_r)]$ in eq. (5) is the receiver site effect needed to obtain the vertical displacement recorded by a receiver located
on the free surface, where $i_r$ is the ray incident angle at the receiver. The angles $i_w$, $i_c$ and $i_r$ (take-off angles of the ray in the ocean, in the crust under the seafloor, and at the receiver, respectively) are related by the Snell's law:

$$ \frac{\sin(i_w)}{\alpha_w} = \frac{\sin(i_c)}{\alpha_c} = \frac{\sin(i_r)}{\alpha_r}. \qquad (6) $$

If pressure sources are distributed along the ocean surface, eq. (5) can be generalized as:

$$ u_Z(\mathbf{x}_s, f) = \iint A(\mathbf{x}_s, \mathbf{x}_r, f) e^{-2i\pi f t_p(\mathbf{x}_s, \mathbf{x}_r)} P(\mathbf{x}_s, f) dS(\mathbf{x}_s), \quad (7) $$

where the amplitude term is given by:

$$ A(\mathbf{x}_s, \mathbf{x}_r, f) = [2 \cos(i_r)] \left[ \frac{\cos(i_c)}{4\pi \rho_c \alpha_c^2} \frac{e^{-\pi f t_p^*}}{\sqrt{J_p^E}} \hat{T}_p^E \right] [2C_p \frac{\rho_c}{\rho_w}], \quad (8) $$

and the double integral is computed over the ocean surface.

If independent pressure sources at the ocean surface are considered, using eqs (2) and (7), one can write the PSD of the P-wave vertical displacement (m² s) at a given seismic frequency f, as:

$$ \langle |u_Z(\mathbf{x}_r, f)|^2 \rangle = \iint |A(\mathbf{x}_s, \mathbf{x}_r, f)|^2 F_p(\mathbf{x}_s, f) dS(\mathbf{x}_s). \quad (9) $$

## 2.4 Beamforming expression of the P-wave vertical displacement

In the following, we use the beamforming approach to detect P waves and measure their amplitudes recorded by the CI network. We define the vertical displacement beam PSD (m² s) at a given angular frequency ω and horizontal slowness vector s, as follows:

$$ B_Z(\omega, \mathbf{s}) = \frac{1}{N_s^2} \left\langle \left| \sum_{j=1}^{N_s} S_j(\omega) e^{-i\omega \mathbf{s} \cdot (\mathbf{x}_j - \mathbf{x}_c)} \right| \right\rangle, \quad (10) $$

where $N_s$ is the station number, $S_j(\omega)$ represents the vertical displacement spectrum $u_Z(\mathbf{x}_j, \omega)$ recorded at the jth station, s is the horizontal slowness vector toward the source, $\mathbf{x}_j$ is the position vector of the jth station and $\mathbf{x}_c$ is the position vector of the network centre.

Assuming that the stations forming the array are not too far away from each other, one can use the following approximations for the amplitude and traveltime of P waves generated by a point source at $\mathbf{x}_s$:

$$ \begin{gathered} A(\mathbf{x}_s, \mathbf{x}_j, f) = A(\mathbf{x}_s, \mathbf{x}_c, f), \\ t_P(\mathbf{x}_s, \mathbf{x}_j) = t_P(\mathbf{x}_s, \mathbf{x}_c) - \mathbf{s}_c \cdot (\mathbf{x}_j - \mathbf{x}_c), \end{gathered} \quad (11) $$

where $\mathbf{s}_c(\mathbf{x}_s, \mathbf{x}_c)$ is the P-wave horizontal slowness vector (pointing toward the source) at the network centre.

The beam PSD (10) of the vertical displacement due to a noise source extended within a small surface S centred at $\mathbf{x}_s$, can be written by using eqs (5), (8) and (11):

$$ B_Z(\omega, \mathbf{s}) = |A(\mathbf{x}_s, \mathbf{x}_c, f)|^2 R(\omega, \mathbf{s} - \mathbf{s}_c) F_p(\mathbf{x}_s, f) S, \quad (12) $$

where $R(\omega, \mathbf{s})$ is the so-called array response,

$$ R(\omega, \mathbf{s}) = \frac{1}{N_s^2} \left| \sum_{j=1}^{N_s} e^{-i\omega\mathbf{s} \cdot (\mathbf{x}_j - \mathbf{x}_c)} \right|^2. \quad (13) $$

Considering independent pressure sources at the ocean surface, and using eqs (7), (8) and (11), we can write the beam PSD of the

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

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Top: geographical map of the CI seismic network in California. Triangles correspond to station locations. Bottom: array response at a period of 5 s.</figcaption></figure>

vertical displacement as:

$$ B_Z(\omega, \mathbf{s}) = \iint |A(\mathbf{x}_s, \mathbf{x}_c, f)|^2 R(\omega, \mathbf{s} - \mathbf{s}_c) F_p(\mathbf{x}_s, f) dS(\mathbf{x}_s). \quad (14) $$

## 3 MODELING OF TYPHOON IOKE P-WAVES

In order to validate our theoretical computation, we analyse P-waves generated by typhoon Ioke (September 2006) in the Western Pacific and recorded by the CI network. This network consists of 48 stations and it has a lateral extent of about 600 km (Fig. 2a). Fig. 2(b) shows the array response at a period of 5 s. If we define the beam slowness resolution as the spike width at half of the beam maximum, we obtain 0.009 s km⁻¹. It corresponds to a spatial resolution of about 1390 km at a distance of 75°. If we consider the beam width at 80 per cent of the maximum, we obtain 0.005 s km⁻¹, which corresponds to a spatial resolution of 830 km.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Power spectral density obtained by beamforming analysis, in m² s, corresponding to 2006 September 4, 12:00–21:00. White circles are plotted along slownesses 0.08 and 0.04 s km⁻¹, which correspond to source-network distances 30° and 90°.</figcaption></figure>

Data processing is as follows. The vertical component seismic data with 1 Hz sampling rate are selected. They are demeaned and deconvolved with the instrument response and converted to displacement. Time-series corresponding to earthquakes of magnitude greater than 5.5 are removed from the data set and visual inspection is performed in order to remove earthquakes of smaller magnitude (if present). For each 512 s time window, the average spectrum amplitude is computed over all stations and data are discarded when their spectrum amplitude is larger than the mean plus twice the standard deviation. This enables to remove glitches and other anomalous signals. Beam PSDs are then computed using eq. (10) with non-overlapping time windows of 512 s and averaging over 9 hr.

Typhoon Ioke was first analysed by Zhang *et al.* (2010a,b) who used a beamforming approach on the CI network data to follow the *P*-wave source generated by the typhoon from 2006 September 1 to 7. Gualtieri *et al.* (2014) further showed that theory predicts that both *P* and *S* waves are generated in the vicinity of the typhoon and at the same location, but only *P* waves could be detected because the *S*-waves amplitude was below the beam noise level. Hereafter, we quantitatively compare observed and modeled *P* waves generated by typhoon Ioke.

Fig. 3 shows the 5 s period beam PSD computed on 2006 September 3 between 12:00 and 21:00. The beam PSD shows a clear single maximum having amplitude 1.928 × 10⁻¹⁵ m² s at horizontal slowness ($s_x, s_y$) = (−0.05, 0.016) s km⁻¹. We observe that the beam maximum area (yellow) is larger than that of the array response, suggesting an extended source area. Fig. 4 (left column) shows the beam PSD variation from 2006 September 2 to 5. The beam maximum amplitude moves over time, it is largest (2.36 × 10⁻¹⁵ m² s) on September 2 and then decreases down to 1.02 × 10⁻¹⁵ m² s on September 5.

Backprojecting slownesses associated with the beam maxima by ray tracing in the IASP91 model (Kennett & Engdahl 1991), we get the source location. Source maps from 2006 September 2 to 5 are plotted in Fig. 5. Red crosses correspond to the beam centroids and red contours correspond to beam amplitudes larger than 0.8 times the corresponding maximum. We observe that the

detected sources are located in the vicinity of the typhoon track (black line), which was computed as the maximum significant wave height location over time. For a given time, the source detected by beamforming analysis is located behind the significant wave height maximum (black dots) and slightly to the left of the typhoon track. Indeed, the typhoon generates ocean waves at a given location at time $t_1$, which propagate and can meet with ocean waves generated by the typhoon at a different location at time $t_2$. The microseism source is located where these waves meet, that is behind the typhoon and in the vicinity of its track (Longuet-Higgins 1953). It is also confirmed by the source model derived from the ocean wave interactions.

Fig. 5 also shows the modeled body wave sources (Pa²m²s). They are computed as the product of the pressure PSD (eq. 1) and the squared modulus of the *P*-wave source site effect (eq. 3). The pressure PSD is obtained from the global scale ocean wave model WAVEWATCH III<sup>R</sup> (Tolman 2009; Ardhuin *et al.* 2011; Rascle & Ardhuin 2013), which has a constant resolution of 0.5° both in latitude and longitude. The ocean wave model enables to take into account all possible ocean wave interactions: within a storm, between distant storms and those generated at the coast between incident and reflected ocean waves. We checked that the pressure source in the area of interest is created only by interactions of ocean waves associated with the typhoon and that sources generated by coastal reflections are negligible. The interactions of ocean waves generated by the typhoon create a pressure source with an elongated shape located not exactly along the typhoon track but further south on September 2 to 4 and further west on September 5. The source elongated shape is related to the curvature of the typhoon track and its speed. Indeed, it is a well-known feature of tropical storms that the ocean waves are very different on the different sides of the storm, as a function of the storm displacement speed. In particular, the ‘left-rear quadrant’ (when looking in the storm displacement direction) usually has crossing seas with ocean waves coming from opposite directions (Wright *et al.* 2001; Holthuijsen *et al.* 2012). This is because the storm usually travels faster than the ocean wave group velocity and thus overtakes the waves previously generated that radiate as swell. In the left-rear quadrant, the ocean waves generated by the winds propagate in the opposite direction to swells coming from a previous position of the storm. The model predicts a maximum pressure source in the region with waves coming from opposite directions, and that maximum is thus shifted with respect to the typhoon position at the same time. This shift is consistent with the beam maximum location observed on Fig. 5.

The modeled pressure PSD is smooth. Multiplying the pressure PSD by the squared modulus of the source site effect, $|C_p|^2$, modifies the source amplitude by a factor 0.77–11.5 and introduces roughness in the source model related to the local bathymetry below the pressure source (Fig. 1). The bathymetry is taken from ETOPO1 model (Amante & Eakins 2009) smoothed on a grid of 0.5° both in latitude and longitude.

Red contours in Fig. 5 correspond to the backprojection of the beam amplitudes larger than 0.8 times the corresponding maximum. Contours are plotted only for comparison with the modeled sources (colour map in Fig. 5) and illustrate that the modeled source has a smaller lateral extent but is consistent with the observed source. The large size of the red contour, about 1000 km wide, is the result of the convolution between the array response and the extended source. The source lateral extent is about 800 km long and the beam array response resolution is about 1000 km.

We computed synthetic beams associated with the modeled source using eq. (14). The amplitude term in eq. (14) is calculated
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Observed (left column) and modeled (right column) power spectral density obtained by beamforming analysis, in m² s, from 2006 September 2 to 5. White arc circles are plotted along slownesses 0.08 and 0.04 s km⁻¹, which correspond to teleseismic epicentral distances 30° and 90°.</figcaption></figure>

by ray tracing in the IASP91 model (Kennett & Engdahl 1991). The attenuation term $t_p^*$ is computed for the attenuation model of Montagner & Kennett (1996) using the software developed by Herrmann (2013). Synthetic beams are plotted in Fig. 4 (right column) for days

2006 September 2 to 5. The observed and modeled beam maximum shows consistent location and lateral extent over time. Therefore, we compare their maximum amplitude variations over time. Error bars are computed as ±2 times the standard deviation of the beam PSD,
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Maps of *P*-wave sources, $F_p|C_p|^2$ in Pa$^2$·m$^2$·s, generated by typhoon Ioke at 5 s period. Red contours are the backprojection of beam amplitudes larger than 0.8 times the corresponding maximum. Red crosses are the beam maximum backprojection. The typhoon track is represented by black line and the typhoon position at 12:00 and 21:00 each day is shown with black dots. In colour, we show the *P*-wave source amplitude computed as the product of the ocean wave model pressure PSD and the squared modulus of the *P*-wave source site effect $|C_p|^2$.</figcaption></figure>

computed in the slowness range ±0.08 s km⁻¹ (e.g. Fig. 3). Fig. 6 shows that the trend is similar between real and synthetic beams and that the order of magnitude is the same. The amplitude of the modeled beam is underestimated by a factor that ranges from 1.2 to 1.9. The discrepancy between the real and modeled beams seems to be correlated with the noise level in Fig. 4. It is also within the uncertainty of the ocean wave model for large storms, and within the uncertainty of the seismic modeling, in which accurate crustal models under the source and under the seismic array have not been taken into account.

We also computed the beam for periods from 3 to 10 s (with a step of 1 s) and we observed a beam maximum with a high signal-to-noise ratio that could be associated with Ioke only at 5 s period. We checked that both the modeled pressure PSD and the source site effect are the largest at period of 5 s with respect to other periods. For example, on September 3, the ratio of modeled pressure PSD at 5 and 4 s (6 s, respectively) is equal to 2 (8, respectively). The squared ratio of the source site effect modulus is about 36 between periods of 5 and 4 s or 6 s. The beam maximum amplitude is proportional to both effects, which explains why the beam maximum is mostly visible with a high signal-to-noise ratio at period of 5 s.

## 4 CONCLUSIONS

Secondary microseism body waves are weak signals that cannot be directly observed on seismograms but can be investigated by beamforming analysis. In this paper, we presented the theory for modeling the secondary microseism *P*-wave amplitude recorded by an array of stations. We showed that the *P*-wave displacement in the far field can be computed in the frequency domain as the product of (1) the pressure source, (2) the source site effect that accounts for the constructive interference of multiply reflected *P* waves in the ocean, (3) the propagation term from the ocean bottom to the seismic stations, (4) the receiver site effect. We applied this approach to model the beam temporal variation corresponding to signals recorded by the CI network in California and generated by the typhoon Ioke in 2006 September. We showed that the beam centroid is located in the vicinity of the pressure source derived from the ocean wave model. Furthermore, we observed that the beam uncertainty area is about 1000 km that is wider than the ocean wave model source area. In the slowness domain, there is a good agreement between the real and modeled beam extent. Although the beam amplitude is slightly underestimated with respect to the one obtained by beamforming analysis, the amplitude difference is
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Observed (red) and modeled (black) beam PSD maximum amplitude generated by typhoon Ioke from 2006 September 2 to 5. Error bars correspond to ±2 standard deviations.</figcaption></figure>

within the uncertainty of the ocean wave model and errors due to the 1-D Earth model. This forward modeling will be used in future studies to invert P-wave noise data and retrieve source intensity and lateral extent. It can also be easily adapted to model secondary microseism S waves or other body wave phases.

## ACKNOWLEDGEMENTS

This work was supported by Agence Nationale de la Recherche grant ANR-14-CE01-0012 MIMOSA and grant ANR-10-LABX-19-01 'LabexMer'. LG acknowledges support from a Lamont-Doherty Earth Observatory Postdoctoral Fellowship and the Brinson Foundation. MS acknowledges MISTERIOS (CGL2013-48601-C2-1R). The authors used ObsPy (Beyreuther et al. 2010) for data preprocessing. Data were downloaded from IRIS-DMC. This is IPGP contribution number 3768 and Lamont-Doherty contribution number 8036. The authors thank the two reviewers for their constructive remarks.

## REFERENCES

* Aki, K. & Richards, P.G., 2002. *Quantitative Seismology*, 2nd edn, University Science Books.

* Amante, C. & Eakins, B.W., 2009. *ETOPO1 1 Arc-minute Global Relief Model: Procedures, Data Sources and Analysis*, US Department of Commerce, National Oceanic and Atmospheric Administration, National Environmental Satellite, Data, and Information Service, National Geophysical Data Center, Marine Geology and Geophysics Division.

* Arduin, F. & Herbers, T., 2013. Noise generation in the solid earth, oceans and atmosphere, from nonlinear interacting surface gravity waves in finite depth, *J. Fluid Mech.*, **716**, 316–348.

* Arduin, F., Stutzmann, E., Schimmel, M. & Mangeney, A., 2011. Ocean wave sources of seismic noise, *J. geophys. Res.*, **116**(C9), C09004, doi:10.1029/2011JC006952.

* Arduin, F., Gualtieri, L. & Stutzmann, E., 2015. How ocean waves rock the earth: two mechanisms explain microseisms with periods 3 to 300 s, *Geophys. Res. Lett.*, **42**, 765–772.

* Berger, J., 2004. Ambient Earth noise: a survey of the Global Seismographic Network, *J. geophys. Res.*, **109**(B11), B11307, doi:10.1029/2004JB003408.

* Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y. & Wassermann, J., 2010. Obspy: a python toolbox for seismology, *Seismol. Res. Lett.*, **81**(3), 530–533.

* Cerveny, V., 2001. *Seismic Ray Theory*, Cambridge Univ. Press.

* Davy, C., Barruol, G., Fontaine, F.R., Sigloch, K. & Stutzmann, E., 2014. Tracking major storms from microseismic and hydroacoustic observations on the seafloor, *Geophys. Res. Lett.*, **41**(24), 8825–8831.

* Fukao, Y., Nishida, K., Suda, N., Nawa, K. & Kobayashi, N., 2002. A theory of the Earth's background free oscillations, *J. geophys. Res.*, **107**(B9), 2206, doi:10.1029/2001JB000153.

* Fukao, Y., Nishida, K. & Kobayashi, N., 2010. Seafloor topography, ocean infragravity waves, and background Love and Rayleigh waves, *J. geophys. Res.*, **115**(B4), B04302, doi:10.1029/2009JB006678.

* Gal, M., Reading, A., Ellingsen, S., Gualtieri, L., Koper, K., Burlacu, R., Tkalčić, H. & Hemer, M., 2015. The frequency dependence and locations of short-period microseisms generated in the Southern Ocean and West Pacific, *J. geophys. Res.*, **120**(8), 5764–5781.

* Gerstoft, P., Fehler, M.C. & Sabra, K.G., 2006. When Katrina hit California, *Geophys. Res. Lett.*, **33**(17), L17308, doi:10.1029/2006GL027270.

* Gerstoft, P., Shearer, P.M., Harmon, N. & Zhang, J., 2008. Global P, PP, and PKP wave microseisms observed from distant storms, *Geophys. Res. Lett.*, **35**(23), L23306, doi:10.1029/2008GL036111.

* Gualtieri, L., Stutzmann, E., Capdeville, Y., Arduin, F., Schimmel, M., Mangeney, A. & Morelli, A., 2013. Modelling secondary microseismic noise by normal mode summation, *Geophys. J. Int.*, **193**(3), 1732–1745.

* Gualtieri, L., Stutzmann, E., Farra, V., Capdeville, Y., Schimmel, M., Arduin, F. & Morelli, A., 2014. Modelling the ocean site effect on seismic noise body waves, *Geophys. J. Int.*, **197**(2), 1096–1106.

* Gualtieri, L., Stutzmann, E., Capdeville, Y., Farra, V., Mangeney, A. & Morelli, A., 2015. On the shaping factors of the secondary microseismic wavefield, *J. geophys. Res.*, **120**(9), 6241–6262.

* Gutenberg, B., 1936. On microseisms, *Bull. seism. Soc. Am.*, **26**(2), 111–117.

* Hasselmann, K., 1963. A statistical analysis of the generation of microseisms, *Rev. Geophys.*, **1**(2), 177–210.

* Haubrich, R.A., Munk, W.H. & Snodgrass, F.E., 1963. Comparative spectra of microseisms and swell, *Bull. seism. Soc. Am.*, **53**(1), 27–37.

* Herrmann, R.B., 2013. Computer programs in seismology: an evolving tool for instruction and research, *Seismol. Res. Lett.*, **84**(6), 1081–1088.

* Holthuijsen, L.H., Powell, M.D. & Pietrzak, J.D., 2012. Wind and waves in extreme hurricanes, *J. geophys. Res.*, **117**, C09003, doi:10.1029/2012JC007983.

* Kennett, B. & Engdahl, B.N., 1991. Traveltimes for global earthquake location and phase identification, *Geophys. J. Int.*, **105**(2), 429–465.

* Koper, K.D. & de Foy, B., 2008. Seasonal anisotropy in short-period seismic noise recorded in South Asia, *Bull. seism. Soc. Am.*, **98**(6), 3033–3045.

* Koper, K.D., de Foy, B. & Benz, H., 2009. Composition and variation of noise recorded at the Yellowknife Seismic Array, 1991–2007, *J. geophys. Res.*, **114**(B10), B10310, doi:10.1029/2009JB006307.

* Koper, K.D., Seats, K. & Benz, H., 2010. On the Composition of Earth's short-period seismic noise field, *Bull. seism. Soc. Am.*, **100**(2), 606–617.

* Kurrle, D. & Widmer-Schnidrig, R., 2008. The horizontal hum of the Earth: a global background of spheroidal and toroidal modes, *Geophys. Res. Lett.*, **35**(6), L06304, doi:10.1029/2007GL033125.

* Lacoss, R., Kelly, E. & Toksöz, M., 1969. Estimation of seismic noise structure using arrays, *Geophysics*, **34**(1), 21–38.

* Landès, M., Hubans, F., Shapiro, N.M., Paul, A. & Campillo, M., 2010. Origin of deep ocean microseisms by using teleseismic body waves, *J. geophys. Res.*, **115**(B5), B05302, doi:10.1029/2009JB006918.

* Longuet-Higgins, M.S., 1950. A theory of the origin of microseisms, *Phil. Trans. R. Soc. Lond., A. Math. Phys. Sci.*, **243**(857), 1–35.

* Longuet-Higgins, M.S., 1953. Can see waves cause microseisms, *Nat. Acad. Sci., Washington Publ.*, (306), 74–86.

* Montagner, J.-P. & Kennett, B., 1996. How to reconcile body-wave and normal-mode reference earth models, *Geophys. J. Int.*, **125**(1), 229–248.

* Nishida, K., 2013. Earth's background free oscillations, *Annu. Rev. Earth Planet. Sci.*, **41**, 719–740.

* Nishida, K., Kawakatsu, H., Fukao, Y. & Obara, K., 2008. Background Love and Rayleigh waves simultaneously generated at the Pacific Ocean floors, *Geophys. Res. Lett.*, **35**, L16307, doi:10.1029/2008GL034753.

* Obrebski, M., Arduin, F., Stutzmann, E. & Schimmel, M., 2013. Detection of microseismic compressional (P-wave) body waves aided by numerical modeling of oceanic noise sources, *J. geophys. Res.*, **118**(8), 4312–4324.
Rascle, N. & Ardhuin, F., 2013. A global wave parameter database for geophysical applications. Part 2: Model validation with improved source term parameterization, *Oceanogr. Meteor.*, **70**, 174–188.

Rkie, J. & Romanowicz, B., 2004. Excitation of earth's continuous free oscillations by atmosphere-ocean-seafloor coupling, *Nature*, **431**(7008), 552–556.

Schulte-Pelkum, V., Earle, P.S. & Vernon, F.L., 2004. Strong directivity of ocean-generated seismic noise, *Geochem. Geophys. Geosyst.*, **5**(3), doi:10.1029/2003GC000520.

Stutzmann, E., Roult, G. & Astiz, L., 2000. GEOSCOPE station noise levels, *Bull. seism. Soc. Am.*, **90**(3), 690–701.

Stutzmann, E., Schimmel, M., Patau, G. & Maggi, A., 2009. Global climate imprint on seismic noise, *Geochem. Geophys. Geosyst.*, **10**(11), doi:10.1029/2009GC002619.

Stutzmann, E., Ardhuin, F., Schimmel, M., Mangeney, A. & Patau, G., 2012. Modelling long-term seismic noise in various environments, *Geophys. J. Int.*, **191**(2), 707–722.

Tanimoto, T., 2007. Excitation of normal modes by non-linear interaction of ocean waves, *Geophys. J. Int.*, **168**(2), 571–582.

Tolman, H.L., 2009. User manual and system documentation of wavewatch III™ version 3.14, Technical Note, MMAB Contribution 276.

Traer, J. & Gerstoft, P., 2014. A unified theory of microseisms and hum, *J. geophys. Res.*, **119**, 3317–3339.

Vinnik, L., 1973. Sources of microseismic p waves, *Pure appl. Geophys.*, **103**(1), 282–289.

Webb, S.C., 1998. Broadband seismology and noise under the ocean, *Rev. Geophys.*, **36**(1), 105–142.

Webb, S.C., 2008. The Earth's hum: the excitation of Earth normal modes by ocean waves, *Geophys. J. Int.*, **174**(2), 542–566.

Wright, C.W. et al., 2001. Hurricane directional wave spectrum spatial variation in the open ocean, *J. Phys. Oceanogr.*, **31**, 2472–2488.

Zhang, J., Gerstoft, P. & Bromirski, P.D., 2010a. Pelagic and coastal sources of P-wave microseisms: generation under tropical cyclones, *Geophys. Res. Lett.*, **37**(15), L15301, doi:10.1029/2010GL044288.

Zhang, J., Gerstoft, P. & Shearer, P.M., 2010b. Resolving P-wave traveltime anomalies using seismic array observations of oceanic storms, *Earth planet. Sci. Lett.*, **292**(3–4), 419–427.

# APPENDIX: RAY-THEORETICAL P-WAVE DISPLACEMENT FOR A PRESSURE SOURCE ACTING ON THE SEA SURFACE

Let us consider a pressure source $P(\mathbf{x}_s, f)$ acting on a small surface $S$ centred at $\mathbf{x}_s$ on the ocean surface. The medium below the source is a horizontal water layer over a smoothly inhomogeneous medium with crust seismic properties under the seafloor. We assume that a receiver is located in the far field. Using the ray theory approach (Cerveny 2001), we can write the vertical component of the displacement due to P waves generated by the vertical force acting on the surface $S$ and transmitted to the receiver as:

$$u_Z(\mathbf{x}_r, f) = 2[\cos(i_r) \cos(i_w)] \left[ \frac{1}{4\pi \rho_w \alpha_w^2} \sqrt{\frac{\rho_w \alpha_w}{\rho_r \alpha_r J_P}} \hat{T}_P^O \right] \times e^{-2i\pi f t_P(\mathbf{x}_s, \mathbf{x}_r)} [P(\mathbf{x}_s, f)S] \quad (A1)$$

where $\rho_w$, $\rho_c$, $\rho_c$, $\alpha_c$, $\rho_r$, $\alpha_r$ are the density and P-wave velocity in the water, in the crust under the seafloor and at the receiver, respectively. $J_P$ is the P-wave geometrical spreading. The angles $i_w$, $i_c$ and $i_r$ are the take-off angles of the ray in the water, in the medium under the seafloor and at the receiver, respectively. $t_P(\mathbf{x}_s, \mathbf{x}_r)$ is the P-wave traveltime between the source at $\mathbf{x}_s$ and the receiver at $\mathbf{x}_r$.

The term $[\cos(i_r) \cos(i_w)]$ is the product of the vertical components of the P-wave polarization vectors at the receiver and at the source, respectively. The factor 2 ahead is due to the fact that the source is on the sea surface.

The second term is the product of the source excitation term, the geometrical spreading term and the coefficient $\hat{T}_P^O$, which is the energy normalized P-wave displacement transmission coefficient at the seafloor. This coefficient is given by:

$$\hat{T}_P^O = \mathcal{T}_P \sqrt{\frac{\rho_c \alpha_c \cos(i_c)}{\rho_w \alpha_w \cos(i_w)}}$$

where $\mathcal{T}_P$ is the P-wave displacement transmission coefficient related to the P-wave potential transmission coefficient $\mathcal{T}_P$, introduced in eq. (3), through the relation $\mathcal{T}_P = \mathcal{T}_P^O \frac{\alpha_w}{\alpha_r}$.

In a spherically layered earth model, the geometrical spreading $J_P$ is given by:

$$\begin{aligned} J_P &= r_r^2 r_w \frac{\cos(i_r) \cos(i_w)}{\alpha_w \sin(i_w)} \sin(\Delta) \frac{d\Delta}{dp} \\ &= r_r^2 r_w^2 \frac{\cos(i_r) \cos(i_w)}{\alpha_w^2} \frac{\sin(\Delta)}{p} \frac{d\Delta}{dp}, \end{aligned}$$

where $r_w$ and $r_r$ are the radius at the source and at the receiver, respectively. $\Delta$ is the great circle distance between the source and the receiver, measured in radians. The ray parameter $p = r \sin(i)/\alpha$ is constant for P-wave rays in spherically layered models.

We can write the vertical displacement (A1) as:

$$u_Z(\mathbf{x}_r, f) = [\cos(i_r) \cos(i_c)] \left[ \frac{1}{4\pi \rho_c \alpha_c^2} \sqrt{\frac{\rho_c \alpha_c}{\rho_r \alpha_r J_P^E}} \right] e^{-2i\pi f t_P(\mathbf{x}_s, \mathbf{x}_r)} \times \left[ 2T_P \frac{\rho_c}{\rho_w} [P(\mathbf{x}_s, f)S] \right] \quad (A2)$$

where

$$J_P^E = r_r^2 r_E^2 \frac{\cos(i_r) \cos(i_c) \sin(\Delta)}{\alpha_c^2 p} \frac{d\Delta}{dp}$$

is the geometrical spreading corresponding to P-wave propagation in a spherically layered earth model without water layer at the top of the model (i.e. the source is located at the top of the crust layer). To obtain eq. (A2), we set $r_w = r_E$, where $r_E$ is the Earth radius, and we assume that the ray parameter does not change by considering the source at the top of the crust layer.

In order to take into account the ocean site effect, the transmission coefficient $T_P$ has to be replaced by $C_P$, given by eq. (3), in expression (A2).

If the P-wave ray path between the source and the receiver goes through discontinuities inside the earth, eq. (A2) has to be multiplied by the product of the normalized P-wave transmission coefficients corresponding to all these discontinuities. This product is denoted by $\hat{T}_P^E$. For the epicentral distance of 73°, $\hat{T}_P^E$ is around 0.96.

Considering the receiver located on the crust free surface, the first term $\cos(i_r)$ in eq. (A2) has to be replaced by $c_r(i_r)$, where $c_r(i_r)$ is the free-surface site effect at the receiver for the vertical displacement due to an incident P wave:

$$c_r(i_r) = \frac{2 \cos(i_r) \cos(2j_r)}{D},$$

with

$$D = \cos^2(2j_r) + \left(\frac{\beta_r}{\alpha_r}\right)^2 \sin(2i_r) \sin(2j_r).$$

The symbol $\beta_r$ denotes the S-wave velocity at the receiver. The angle $j_r$ is obtained from the Snell's law at the free surface, $\frac{\sin(i_r)}{\alpha_r} = \frac{\sin(j_r)}{\beta_r}$.
For small values of the angle $i_r$, $c_r(i_r)$ can be approximated by $2 \cos(i_r)$. Considering the P-wave incident angle $i_r = 18^\circ$ at the receiver, which corresponds to the epicentral distance of 73°, the approximation of $c_r(i_r)$ by $2 \cos(i_r)$ leads to the error of 1 per cent.

Finally, one can write the vertical component of the P-wave displacement at a receiver on the free surface in the far field as:

$$u_z(\mathbf{x}_r, f) = [2 \cos(i_r) \cos(i_c)] \left[ \frac{1}{4\pi \rho_c \alpha_c^2} \frac{e^{-\pi f t_p^*}}{\sqrt{J_p^E}} \hat{T}_p^E \right] e^{-2i\pi f t_p(\mathbf{x}_s, \mathbf{x}_r)} \\ \times \left[ 2C_p \frac{\rho_c}{\rho_w} \right] [P(\mathbf{x}_s, f)S], \quad (A3)$$

where we use the fact that $\rho_r = \rho_c$ and $\alpha_r = \alpha_c$ and introduce the seismic attenuation term $t_p^*$ along the ray path.

In eq. (A3), the geometrical spreading $J_p^E$ is given by:

$$J_p^E = r_E^4 \frac{\cos(i_r) \cos(i_c) \sin(\Delta)}{\alpha_c^2} \frac{d\Delta}{dp} \quad (A4)$$

Let us remark that, when modeling P wave in the far field, multiplying the pressure source with the source site effect, $[2 C_p(i_w, h, \omega) \frac{\rho_c}{\rho_w}] [P(\mathbf{x}_s, f)S]$, is equivalent to consider a vertical force source at the top of the crust layer and propagation of the generated P waves through an Earth model without ocean.