*2.1. Ambient Noise Data*

The selected ambient noise data were recorded at 33 broadband seismic stations (Figure 1, see IRIS Data Management Center for further details of the instruments) of the Transportable Array (TA) network along parts of the middle and southeastern North Atlantic coastal area and the Shenandoah Valley. The original time series of amplitude in vertical (V), north-south (N), and east-west (E) directions archived with a sampling rate of 40 samples-per-second were first parsed into 1-h segments, followed by removing the mean, linear trend, and instrument response in each segment [44]. Then each segment was processed following a 14-step procedure summarized below, to estimate the power spectral density (PSD) in vertical (V) direction (steps 1–4), and the primary vibration direction by the radial-to-transverse spectral ratio (*Ra*) method (steps 5–9) as well as the polarization analysis method based on [45] (steps 10–13) in the three DF bands (DF1, 0.1–0.2 Hz; DF2, 0.2–0.3 Hz; and DF3, 0.3–0.4 Hz). The two methods used to estimate the primary vibration directions are both based

on an assumption that DF microseisms propagate dominantly as a fundamental mode Rayleigh wave [34,46–48]. The *Ra* method searches the direction of the largest ratio of radial to transverse components on the horizontal plane which is expected to be along the propagation direction of a Rayleigh wave [29,44]. The polarization analysis method represents the particle motion within a short time range (when the propagation direction is reasonably stable) as an ellipsoid with three axes perpendicular to each other [45], from which the back azimuth of the major axis with highest probability for the whole recording period can be calculated. If the primary energy source is stable and strong enough, and significantly larger than secondary sources, the primary vibration directions obtained by the two methods would agree because the secondary sources do not alter the direction of the major axis but increase the *Ra* values in directions other than the major axis. based on an assumption that DF microseisms propagate dominantly as a fundamental mode Rayleigh wave [34,46–48]. The *Ra* method searches the direction of the largest ratio of radial to transverse components on the horizontal plane which is expected to be along the propagation direction of a Rayleigh wave [29,44]. The polarization analysis method represents the particle motion within a short time range (when the propagation direction is reasonably stable) as an ellipsoid with three axes perpendicular to each other [45], from which the back azimuth of the major axis with highest probability for the whole recording period can be calculated. If the primary energy source is stable and strong enough, and significantly larger than secondary sources, the primary vibration directions obtained by the two methods would agree because the secondary sources do not alter the direction of the major axis but increase the *Ra* values in directions other than the major axis.

*J. Mar. Sci. Eng.* **2020**, *8*, x FOR PEER REVIEW 3 of 21

**Figure 1.** Study area and locations of transportable array (TA) stations (triangles) and buoys (yellow circles) of National Oceanic and Atmospheric Administration, recordings of which are analyzed in this study. The table lists the original numbers and a simplified numbering scheme of the ocean buoys grouped according to their four distinct locations: deep ocean (DO), the continental slope-deep ocean side (SlDO), the continental slope-shelf side (SlSh), and the continental shelf (Sh). The color-encoded relief base map is from [49]. The black lines A-A' and B-B' are the transects presented in Figure 9a,b. The white box indicates the area where the underground shear velocity model (Figure 9c,d) was computed based on [15]. **Figure 1.** Study area and locations of transportable array (TA) stations (triangles) and buoys (yellow circles) of National Oceanic and Atmospheric Administration, recordings of which are analyzed in this study. The table lists the original numbers and a simplified numbering scheme of the ocean buoys grouped according to their four distinct locations: deep ocean (DO), the continental slope-deep ocean side (SlDO), the continental slope-shelf side (SlSh), and the continental shelf (Sh). The color-encoded relief base map is from [49]. The black lines A-A0 and B-B0 are the transects presented in Figure 9a,b. The white box indicates the area where the underground shear velocity model (Figure 9c,d) was computed based on [15].

The steps of the data analysis are: The steps of the data analysis are:

method with a bandwidth coefficient of 40 [51].

(1) Apply an anti-triggering algorithm based on a prescribed range of short (1 s) to long (30 s) term average amplitude ratios (0.2 < STA/LTA < 2.5) to filter each segment for avoiding occasional energy bursts [4,50]. (2) Apply fast Fourier transform with a 10% cosine taper on the filtered segments in three (1) Apply an anti-triggering algorithm based on a prescribed range of short (1 s) to long (30 s) term average amplitude ratios (0.2 < STA/LTA < 2.5) to filter each segment for avoiding occasional energy bursts [4,50].

directions to calculate spectra ((), (), and ()) and then smooth them using Konno–Ohmachi

(3) Compute the PSDs in the vertical direction [52] in the unit of (m/s2)2/Hz dB:


$$PSD(f) = 10 \log \left[ \frac{1}{0.825} \cdot \frac{2\Delta t}{N} \cdot V(f)^2 \right] \tag{1}$$

where ∆*t* is the sample interval (0.01 s); *N* is the number of samples in each selected time-series segment; the constant 1/0.825 is a scale factor to correct for the 10% cosine taper applied [53].


$$
\begin{bmatrix} R(t, F, \varphi) \\ T(t, F, \varphi) \end{bmatrix} = \begin{bmatrix} -\cos(\varphi) & -\sin(\varphi) \\ -\sin(\varphi) & -\cos(\varphi) \end{bmatrix} \begin{bmatrix} N(t, F) \\ E(t, F) \end{bmatrix} \tag{2}
$$

in which ϕ is defined as the back-azimuth angle between the north and the radial direction from the recording station toward the source.


$$\mathbf{M}\_{ij} = \begin{bmatrix} \text{cov}(V\_{ij\prime}, V\_{ij}) & \text{cov}(V\_{ij\prime}, N\_{ij}) & \text{cov}(V\_{ij\prime}, E\_{ij}) \\ \text{cov}(\mathbf{N}\_{ij\prime}, V\_{ij}) & \text{cov}(\mathbf{N}\_{ij\prime}, N\_{ij}) & \text{cov}(\mathbf{N}\_{ij\prime}, E\_{ij}) \\ \text{cov}(\mathbf{E}\_{ij\prime}, V\_{ij}) & \text{cov}(\mathbf{E}\_{ij\prime}, N\_{ij}) & \text{cov}(\mathbf{E}\_{ij\prime}, \mathbf{E}\_{ij}) \end{bmatrix} \tag{3}$$

(11) Calculate the three eigenvalues λ*ij* and associated eigenvectors <sup>→</sup> *x ij* of the covariance matrix **M***ij* by solving:

$$\mathbf{M}\_{\rm ij} \overrightarrow{\mathbf{x}}\_{\rm ij} = \lambda\_{\rm ij} \overrightarrow{\mathbf{x}}\_{\rm ij} \tag{4}$$

and define the maximum eigenvalue <sup>λ</sup>*ij*1, and associated eigenvectors *<sup>x</sup>ij*<sup>11</sup> *<sup>x</sup>ij*<sup>12</sup> *<sup>x</sup>ij*<sup>13</sup> <sup>|</sup> . (12) Find the back azimuth angle ϕ*ij* corresponding to the major axis of the polarized ellipse:

$$\begin{cases} \boldsymbol{\varrho}\_{ij} = \arctan \frac{\boldsymbol{\chi}\_{ij\_{13}}}{\boldsymbol{\chi}\_{ij\_{12}}} & \text{if } \boldsymbol{\chi}\_{ij\_{11}} > 0 \\ \boldsymbol{\varrho}\_{ij} = \arctan \frac{\boldsymbol{\chi}\_{ij\_{13}}}{\boldsymbol{\chi}\_{ij\_{12}}} + 180 & \text{if } \boldsymbol{\chi}\_{ij\_{11}} < 0 \end{cases} \tag{5}$$

Compute the probability of the back azimuth angle within 0◦–360◦ range with a 10◦ bin width. The back azimuth of highest probability is considered as the primary vibration direction by the polarization method.
