Cyclone Signatures in the South-West Indian Ocean from Two Decades of Microseismic Noise

Tropical Cyclones (TC) represent the most destructive natural disaster affecting the islands in the South-West Indian Ocean (SWIO) each year. Monitoring ocean activity is therefore of primary importance to secure lands, infrastructures and peoples, but the little number of oceanographic instruments makes it challenging, particularly in real time. Long-term seismological records provide a way to decipher and quantify the past cyclonic activity by analyzing microseisms, seismic waves generated by the ocean activity and propagating through the solid Earth. In the present study, we analyze this microseismic noise generated by cyclones that develop in the SWIO basin between 1999 and 2020, using broadband seismic stations in La Réunion. The power spectral density (PSD), together with the root mean square (RMS) analyses of continuous seismic data recorded by the permanent Geoscope RER seismic station, indicate the intensification of the microseismic noise amplitude in proportion to the cyclone intensity. Thus, we establish a relationship between the cyclone intensity and the PSD of the Secondary Microseisms (SM) in frequency band ∼0.14 to 0.25 Hz (4 to 7 s period). The Pearson coefficient between the observed and estimated TC intensity are >0.8 in the presence of a cyclone with mean wind speeds >75 km/h and with a seismic station distance-to-storm center D < 3000 km. A polarization analysis in the time and frequency domains allows the retrieval of the backazimuth of the SM sources during isolated cyclone events and well-polarized signal, i.e., CpH > 0.6. We also analyzed the RMS of the Primary Microseisms (PM frequency between ∼0.05 and 0.1 Hz, i.e., for 10 to 20 s period) for cyclones passing nearby La Réunion (D < 500 km), using the available temporary and permanent broadband seismic stations. We also found high correlation coefficients (>0.8) between the PM amplitude and the local wave height issued from the global hindcast model demonstrating that the PM amplitude can be used as a robust proxy to perform a real-time wave-height monitoring in the neighboring ocean. Transfer functions are calculated for several cyclones to infer wave height from the seismic noise amplitude recorded on land. From the analysis of two decades of data, our results suggest that it is possible to quantify the past ocean activity for as long as continuous seismic archives are available, emphasizing microseismic noise as a key observable for quantifying and understanding the climate change.


Introduction
The South-West Indian Ocean (SWIO) is an active ocean basin in terms of Tropical Cyclones (TCs). Each year, during the cyclonic season (November-April), TCs cause extensive damages due to the effects of rain, wind and waves, not only to the islands located in the SWIO such as Madagascar, La Réunion and Mauritius, but also to the continental countries along the east coast of Africa such as Mozambique (e.g., TC Idai in 2019). Coastal areas of those countries are particularly vulnerable to the waves generated by the TCs or storms, motivating the approach of better forecasting and monitoring the ocean wave activity. The little number of direct wave observations available in the SWIO is, however, a strongly limiting factor in the present or past TCs analyses. In this general frame, our objective in this paper is to explore seismic data to extract proxies of the ocean activity, providing new observables associated with TC activities in the SWIO and hence a way to quantify present or past TC through the analysis of real-time or archived time series. In our approach, this is done by analyzing the microseismic noise recorded by seismic stations on La Réunion for TCs that occurred in the past 20 years Figure 1. Such longterm microseismic noise characterization is crucial for learning more details about the ocean wave activity in the past TCs, especially in places with limited oceanographic direct observations. Recovering past information related to TCs are also of major interest to assess long-term climate changes e.g., [1].
The ocean-induced microseismic noise also called microseisms reveals the continuous interactions between the atmosphere, the ocean and the solid Earth. Microseisms correspond to continuous vibrations propagating as elastic waves in the solid Earth, induced by the oceanic gravity waves e.g., [2]. They are recorded worldwide on land and on the ocean bottom by the broadband seismic stations in the frequency range ∼0.05 to 1 Hz. There are two types of microseisms, which differ from their frequencies and origins: the Primary Microseisms (hereafter called PM, with frequency ∼0.05-0.1 Hz) and Secondary Microseisms (SM, with frequency range ∼0.1 to 1 Hz). PM have the same period as the ocean wave e.g., [3] and have been used as a wave proxy e.g., [4,5]. They are accepted to be generated in the shallow waters through the interaction between the sloping seafloor pressurization by ocean waves e.g., [6,7]. On the other hand, SM have a frequency range twice those of the ocean waves and are generated both in the shallow and deep waters. SM are accepted to be induced by non-linear interactions between ocean waves of similar periods that propagate in opposite directions generating local standing waves e.g., [2]. SM sources are primarily located in the open oceans where two distant swells may interact e.g., [8][9][10][11] and may generate standing waves allowing pressure waves to be transmitted to the ocean bottom. In addition, SM are also created in the vicinity of storm centers due to the interactions of two waves from opposite directions, which are generated by the same storm at different times e.g., [8,9,12]. SM sources have also been detected in coastal regions where the sea waves reflected at the coastal slope may interact with the incoming swell e.g., [5,8,13,14]. This last origin of SM is also known as long period secondary microseisms (LPSM).
Several works used the microseismic noise to analyze historical TCs e.g., [4,9,13,[15][16][17][18]. Gualtieri et al. [17] analyzed 22 years of seismic data to estimate the relationship between the TC in the Pacific Ocean and the SM recorded by seismic stations located in the northwest of Pacific. Several studies involving the use of the microseismic data have shown the feasibility of tracking the TC using the microseismic noise e.g., [9,16,[18][19][20][21]. Apart from the works of Davy et al. [4,9] and Barruol et al. [13], none of these studies focused on the cyclone activity in the SWIO.
In the current work, we make use of the microseismic noise recorded by terrestrial seismic stations on La Réunion Figure 1a to study the cyclonic activity in the SWIO since 1999. We first use the seismic station RER (network, G, Figure 1a) that operated since 1986. We investigate the continuous microseisms data recorded by the RER station and correlate them with cyclone data (position, strength, wind). To further analyze the cyclonic events that passed at distance <500 km from the island, we also investigate the relatively recent (operated since ∼2010) seismic data from permanent (network code PF, Figure 1a) and temporary networks (2015-2021, network code ZF, Figure 1a) on La Réunion. In order to retrieve wave height from seismic data, and therefore to use a seismic station as a wave gauge, we establish relationships (or transfer functions) between the PM amplitude and significant wave height from the global wave model GOW2 [22]. , from the temporary experiment (blue triangles) and from the Geoscope RER station (magenta triangle). Colored stars show the offshore location of the nodes at which the significant wave heights (H S ) were extracted from GOW2 wave model Perez et al. [22]. (b) Number of tropical cyclones/storms, per season, between 1999 and 2020. (c) TC tracks (continuous grey lines) in the South-West Indian Ocean during the cyclonic period between 1999 and 2020. The color of each dot indicate the cyclone intensity on a 6-hours basis. The magenta triangle marks the RER seismic station.

Tropical Cyclones and Microseismic Noise Data
Cyclones develop in the SWIO each year between October and April. In this work, the word cyclone is used in a broad sense to refer to a tropical cyclone or to a tropical storm. We analyzed seismic signatures of cyclones that occurred in the SWIO between 1999 and 2020 (Figure 1b,c). The number of cyclones for each season is plotted in Figure 1a. The cyclonic seasons 2005-2006, 2010-2011, 2016-2017 and 2007-2008, 2018-2019 have the minimum (6) and maximum (16) number of cyclones, respectively. The information about the cyclone in the SWIO is freely available from the Météo France (MF) website (http:// www.meteo.fr/temps/domtom/La_Reunion/webcmrs9.0/anglais/index.html, accessed on 12 April 2021). Each cyclone is characterized by the intensity and pressure at the storm center, along with its geographical position. This information is made available on a 6 h basis during the cyclone lifespan. The intensity of the cyclone in this work is defined by the average of the maximum wind (km/h, see Table S1). More details about cyclone names and classifications can be found in Leroux et al. [23] and in the MF database. Figure 1c indicates that typical cyclones start at latitude ∼10 • S, move SW, reinforce down to latitude of 15 to 20 • S and then turn toward the SE to finally decrease in intensity. Strong cyclones form usually in the north, north-east and east of La Réunion. Individual information (i.e., track and intensity) for all the cyclones occurring between 1999 and 2020 are plotted in Figure 1c.
To assess the long-term cyclone climatology, we used seismic station RER to estimate the relationship between the cyclone activity and the microseismic noise. This station is in the south-eastern part of La Réunion and has been operating since 1986 (Figure 1a). RER has been instrumented with a broadband seismic sensor, with a sampling rate of 20 Hz since 1990. However, only since July 1999, the East (BHE) North (BHN) and Vertical (BHZ) components were continuously archived. Thus, to have a uniform data set, we focus our analysis on the data between 1999 and 2020.
To focus on some particular cyclones passing near La Réunion (Section 3.2), we also used data recorded by the seismic stations from the permanent network managed by the Observatoire Volcanologique du Piton de la Fournaise (code PF, OVPF/IPGP, 16 broadband seismometers, shown in red triangles in Figure 1a) and from the temporary network of the "Rivière des Pluies" experiment (code ZF, [24], 10 broadband seismometers plotted in blue triangles in Figure 1a). The corresponding seismic data are available at the French RESIF data portal center (http://seismology.resif.fr, accessed on 12 April 2021.) under their respective FDSN (Federation of Digital Sesimograph Networks) network codes. These stations are equipped with broadband seismometers with sampling rate of 100 Hz, with three components HHE (East), HHN (North) and HHZ (vertical).

Seismic Data Analyses
To estimate the relationship between the cyclone activity and the microseismic noise, we computed the continuous Power Spectral Density (PSD) of the vertical component (BHZ) of the seismic station RER. We first selected 1-h long seismic data with 50% overlap. Then each 1-h time series was divided into 13 segments with 75% overlap with neighboring segments. After that, each segment was transformed into the time-frequency domain using the method of McNamara and Buland [25]. The obtained PSD was converted into decibel (dB) with respect to acceleration (i.e., with respect to m 2 /s 2 /Hz). Finally, we averaged data every 6-h to have a similar time step as the cyclone intensity data, i.e., the average of the maximum wind. During our computation, we used SM with frequency ranges 0.14 to 0.25 Hz (period 4-7 s, short period SM), as it has proven to correlate well with a cyclone e.g., [16][17][18].
To evaluate the backazimuth (BAZ) of the incoming noise, we performed a polarization analysis on successive one hour-long seismic time series. The three components of the seismograms (BHE/HHE, BHN/HHN, BHZ/HHZ) were detrended and decimated to 1 Hz sampling rate before converting them into ground velocity by removing the station response. We obtained different parameters such as CpH and CpZ which correspond to the degree of polarization in the horizontal and vertical plane, respectively, and the azimuth of the ellipse long axis (between 0 and 360 • ). These different parameters allow us to characterize the full 3-D ground motion. A CpH of 1 corresponds to a perfectly polarized signal in the corresponding plane (i.e., a linear ground motion), while a value of 0 characterizes a random (circular) ground motion. More details about these parameters are presented in Section S1 and can be found in previous works [5,7,26,27]. In this work, we focused our attention to the azimuth of significantly polarized signal corresponding to CpH > 0. 6. To quantify the relationship between the PM and the significant wave height H S through a transfer function, we calculate the hourly Root Mean Square (RMS) amplitude of the microseismic noise on the vertical component (HHZ) of the seismograph. The vertical component was used as it is generally less affected by artefacts and has relatively larger amplitudes than the horizontal components in this frequency band e.g., [4,7]. We first converted the time series into displacement (units of µm) by removing the instrumental response and divided the data into 1-hour segments. We then filtered the data using a 2nd order Butterworth bandpass filter, with corner frequencies at 0.05 and 0.1 Hz (periods ∼10 and 20 s) and computed the RMS. After that, we calculated the Pearson correlation coefficient (P coe f ) between the PM amplitude and the significant wave height H S , extracted from the global wave hindcast model GOW2 developed by Perez et al. [22], at 17 selected nodes locations around the island Figure 1a. Finally, we computed transfer functions between PM (observed seismic) and H S (modeled ocean waves) amplitudes at 3 nodes located close to the coastal area (55. 25 Figure 2b, we plot the PSD at station RER (expressed in dB with respect to acceleration) calculated in the presence (continuous colored lines) and absence (bold black dashed-dotted line) of cyclone activity, during the cyclonic season 2018-2019. We computed the daily average PSD, using the sacpsd command [28], during the cyclone lifespan, for the days that have a maximum wind intensity (V TC ) ≥ 90 km/h (i.e., when it reaches the severe tropical storm stage, see Table S1). Then, we calculated the median of the obtained PSD for each cyclone and presented them as continuous colored lines in Figure 2b. Each spectrum represents a cyclone, the track of which is plotted on the map Figure 2a using the same color code. The black dashed-dotted line indicates the median of daily average during the quiet period in 2019, i.e., in the absence of cyclone or austral swell activities. Austral swells are indeed generated by strong winds and long-distance fetches, in general by depressions moving around Antarctica, which affected La Réunion mainly in winter season, between April and October e.g., [5,29]. Some of them occurring in September and October 2018 are indicated by the pink boxes in Figure 2c. During cyclones, earthquake signatures were removed from the data to avoid the alteration of the PSD. However, in absence of cyclone, the largeamplitude signal generated by an earthquake lasts generally less than a minute and has little to none effect on the computed PSD. Figure 2b shows increase of the PSD amplitude in the presence of a cyclone, in the range ∼5 to 35 dB, depending on the intensity and the location of the storm center. Despite the absence of direct impact of any cyclone on La Réunion during the selected cyclonic season 2018-2019, the PM peaks are also well expressed in the various spectra, indicating the effect of distant swell generated within the cyclone and impacting coastal areas of the island. The amplitude increase of the PM can reach 15 dB in case of Alcide (shown in blue color).
To evaluate the temporal evolution of the cyclone signature, we computed the continuous PSD of the microseisms. As an example, the 2018-2019 cyclonic season shown in Figure 2c,d indicates that cyclones developing at distance as far as ∼3000 km induce clear noise intensification in the SM frequency band. Occasional transient noises observed on the spectra result from either earthquakes or glitches in the data. The presence of high level of energy in the SM frequency band, in absence of a cyclone, is generally related to the continuous storm activity around Antarctica that are well recorded by the seismic station on La Réunion [29]. Such episodes occur mostly during austral winters (April to October). To avoid overcrowded of the Figure 2c, we only marked the austral swell in October, which are within the magenta boxes ( Figure 2c).   16), which makes them difficult to isolate but from March 10 to 13 (TC stage), Figure 2d shows a clear weak signal despite Idai was at its maximum power. Interestingly, TC Desmond and Eketsang, both also developing in the Mozambique Channel, are characterized by similarly weak signals at RER during their lifespans. A first hypothesis is that the weak SM amplitude observed in the Mozambique Channel may be related to a fetch of limited extent, which (and wind speed) controlled the size of significant wave height. The shape, size and location of the Mozambique Channel between the two continental masses of Africa and Madagascar could indeed represent limiting factors for long fetches development. A second, non-exclusive, factor for explaining the absence of SM recorded at RER during Idai is the absence of SM sources in this TC and all other TC travelling in the Mozambique Chan-nel. The complex geometry of the channel could explain the difficulty in generating waves propagating in the opposite directions with the same period that could create standing waves and therefore SM sources. A third hypothesis could be related to a bathymetry effect and is also found from SM source modeling in this area e.g., [8] and polarization analysis e.g., [29], according that the SM amplitude is proportional to the square of the ocean wave height e.g., [2] and induced by pressure fluctuation e.g., [8]. A model of the power spectral density of the equivalent surface pressure ( f p2s) in the Mozambique Channel could have lent evidence of the lack of SM source there. However, the f p2s parameter is not available on the GOW2 wave model that we used in the current study.

Correlation between Secondary Microseisms and Cyclone Distance
Visual comparison between Figure 3a (Cyclone intensity) and Figure 3b (SM PSD amplitude) suggests that the strength of the microseismic noise depends at first on the distance between the storm center and the seismic station. In the area surrounding Réunion Island, even the relatively weak cyclones (V TC < 100 km/h), sign with PSD higher than −115.0 dB. In contrast, between the longitude 70 • E and 90 • E, many powerful cyclones (V TC > 100 km/h) occurred, but are associated with a low level of the PSD at RER (∼−118.0 dB). Figure 3 indicates that a seismic station starts to record a weak but clear SM from a strong TC at a distance ∼4000 km, and that the SM intensified inversely proportional to distance. For each cyclone occurring from 1999 to 2020, we estimated the relation between cyclone distance and the intensification of the SM on the PSD. The obtained equations are plotted in grey lines in Figure 4a. These plots indicate logarithmic SM decays with respect to distance, such that PSD SM = A × log(D) + B, in which PSD SM is the strength of the SM, D is the distance between the storm center and the seismic station, and A and B two constants depending on the cyclone intensity and how fast it evolves from tropical disturbance to tropical cyclone/storm (see Table S1, for cyclone classification). Average values of A and B were computed from the different cyclones, with their respective standard deviations. We first determined the standard deviation from all data and removed the outlier cyclones (representing 6% of our data), which had a variance of more than three times the standard deviation. After removing the outliers, we recomputed the averages and obtained values equal to −3.9 ± 1.5 and −89.5 ± 9 (red curve in Figure 4a) for the parameters A and B, respectively.  We also estimated the relationship between the distance and the strength of the PSD in the SM frequency band, using all cyclones occurring between 1999 and 2020 and with tracks geolocalized within latitudes 5 • S-40 • S and longitudes 50 • E-90 • E Figure 4a. We excluded cyclones developing in the Mozambique Channel as the recorded SM appear to be always weaker for a given distance (and for different cyclone strengths, see Section 3.1.1). We derived the following relation: PSD SM = −4.41 × log(D) − 85.68 (green curve in Figure 4a). This equation confirms the direct influence of the distance of the storm center on the recorded SM as already proposed by Davy et al. [9].  Intensity   Figures 3 and 4 show that the intensification of the SM is not only depending on the distance from the storm center, but also on the cyclone intensity. Despite a widespread in the measurements, Figure 4a shows a general decrease of the recorded energy at the station RER as the cyclone is at larger distance (as discussed in Section 3.1.2). Alternatively, Figure 4b displays a general increase of the seismic noise amplitude as the cyclone is stronger. Distant cyclones (red colors Figure 4b) display generally low PSD whereas cyclone developing at short distance (in blue) are generally characterized by higher energies. Exceptions are observed while the cyclone has a V TC < 50 km/h. In such a case, despite the small distance station-to-storm center <500 km, the recorded SM are weak (Figure 4b). Figure 5 presents the relationship between the cyclone intensity and the SM PSD at different geographical locations. We divided the cyclone data into three groups according to their distance to RER and their locations in the SWIO basin. The latitude and longitude boundaries for each group are indicated on the corresponding plot in Figure 5 (left panels), presented on the maps (right panels) and listed in Table 1. To determine the relationship between the cyclone intensity and the microseismic noise, we correlate the cyclone intensity through the proxy of wind velocity, V TC , with the SM PSD (PSD SM ). Statistically, the distribution of the V TC can be approximated with a gamma distribution with relatively high standard deviations ( Figure 5, top right panel). As suggested by Agresti [31], it is more accurate to apply the generalized linear model (GLM) technique to determine the relationship between the two parameters. We therefore applied such GLM technique to estimate the transfer function between the V TC and PSD SM , for the cyclone intensity and microseismic noise data between 1999 and 2017. Similar technique was successfully used in Gualtieri et al. [17] to assess the relationship between the cyclone intensity and the spectral characteristic of the microseismic ambient noise in the northwest of the Pacific Ocean. Please note that if two cyclones overlap or if a cyclone occurs simultaneously with an austral swell, the computed PSD SM represents their cumulative effects and may bias the transfer function. The Pearson correlation coefficients and the transfer functions between the observed V TC and the PSD SM , for group 1,2, 3 and 4 are presented in Table 1. To validate these relations, we estimated the cyclone intensity for the cyclonic seasons 2017-2018 and 2018-2019 and compared them with the observed intensity ( Figure 6), using the corresponding transfer function for each group. The fit between the observed and predicted cyclones are generally good, especially when the cyclone is isolated in time and space from other events (e.g., Dumazile in March 2018 and Fakir in April 2018, Figure 6c). However, the TC intensity is generally over-estimated when two cyclones occur simultaneously in the studied area (e.g., Funani-Gelena in February 2019, Figure 6d). For the cyclones at a distance ≥3000 km, the intensities of the cyclones appear to be underestimated (e.g., Bouchra in November 2018, Figure 6d). Therefore, these equations are only valid for the cyclone with storm center at a distance <3000 km from the seismic station and more accurate in the occurrence of isolated cyclone and Pearson coefficient >0.3. Table 1. Correlation between the SM PSD (PSD SM ) and the cyclone intensity (V TC ).

Group
Longitude  The estimated TC intensity shown in blue, green and cyan dots were obtained from the transfer functions for a geographical location defined by group 2, group 3 and group 4, respectively (as defined in Figure 5 and Table 1). Dashed colored lines indicate the distance between the storm center and the seismic station RER (right axis). The Yellow dots (in Eketsang) color were estimated using the transfer function using all data (i.e., in Figure 5a). (d) Same as c but for the cyclonic season 2018-2019.
Delays (a few hours to a few days) between the observed and estimated V TC are seen for some cyclones (e.g., Berguitta in January 2018, Figure 6c). Gualtieri et al. [17] made similar observation and attributed it to the non-linear coupling between atmosphere and ocean, and a potentially slow wind-wave growth which may take from a few hours to a few days to develop as suggested by Hasselmann et al. [32]. The Pearson correlation coefficient between the observed and predicted cyclone intensity is higher than 0.8, when the observed V TC ≥ 75 km/h, i.e., in the presence of tropical disturbance or higher categories, apart from Cebile (January 2018, Figure 6c), while the estimated V TC is very weak. This could be associated with the presence of slow wind-wave growth. We can see that at the end of Cebile, the estimated V TC begins to increase.
As already mentioned in Section 3.1.1, for the cyclones located in the Mozambique Channel, despite their powers and their distances less than 2000 km from Réunion Island, the SM recorded by RER seismic station remains surprisingly weak (Figure 5d). Very low correlation between the V TC and the PSD SM is observed (P coe f = 0.05). As proposed in Section 3.1.1, the weakness or even the absence of SM sources associated with TCs in the Mozambique Channel may be related to its small size impeding long fetches and to its complex structure limiting the possibility of developing standing waves.

Secondary Microseisms Source Locations
In this section, we focus on the results corresponding to well-polarized signals (CpH ≥ 0.6, see Section 2.2). We computed the polarization analysis at RER seismic station, for the cyclonic seasons 2017-2018 and 2018-2019. The obtained CpH, CpZ and BAZ are plotted in Figure S5 and indicate that in the absence of cyclonic activity in the Indian ocean, the SM recorded at seismic station RER display BAZ ranging between ∼130 • and 180 • , centered at ∼150 ± 10 • . This overall polarization can be associated with the dominant long-distance storm activity, occurring in particular in the southernmost part of the Indian Ocean, well recorded by the seismic station in La Réunion all over the year, but more frequently between April and October e.g., [29]. The low CpH (i.e., poorly polarized signal) observed in the presence of some cyclones is most likely due to multiple simultaneous sources of SM in the basin, such as the simultaneous occurrence of several storms and/or to simultaneous cyclone and austral swells. It is also very likely that other SM sources from other ocean basins may be recorded at RER, arriving with different azimuths and reducing the general strength of the ground motion polarization in this frequency range. During low CpH periods, the obtained BAZ are likely to not represent the BAZ of a single source but an average of simultaneous sources, hence without any physical signification. Figure 7 shows, however, that in the presence of well-polarized seismic signals occurring during a cyclone passage, the polarization (and hence the SM sources) points toward the storm center (e.g., Dumazile, Fakir and Cilida), consistent with the MF data. This agreement occurs dominantly when the cyclone is at the south or south-east of the island. This may indicate that the cyclonic swell direction interacts with long-distance swells at ∼150 • , but further investigations and modeling would be necessary to better understand these characteristics.
Among the 22 cyclones that occurred between 2017 and 2019, five cyclones (Berguitta, Dumazile, Eliakim, Fakir and Cilida) passed close to the island. However, only Dumazile, Eliakim, Fakir and Cilida have a CpH > 0.6 during part of their lifespan ( Figure S5). For these cyclones, we computed the polarization at all available seismic stations in La Réunion, from the permanents (G and PF) and the temporary (ZF) seismic networks. The obtained BAZ, of Dumazile, Fakir and Cilida are presented in Figure 7 and point consistently towards the SSE or the SE, i.e., towards the storm center. The results for Eliakim are not presented here as they are already presented and discussed in detail in another work [Bousquet et al, in this special issue] but points toward the SW, also in good agreement with the storm center. In Figure 7 we compare the average theoretical BAZ issued from the storm center locations provided by MF (pink arrows) with the average seismically derived BAZ for Dumazile (lime arrow, Figure 7b), Fakir (Blue arrows, Figure 7c) and Cilida (green arrows, Figure 7d), obtained from the polarization analysis, for the various seismic stations in La Réunion. Figure 7 indicates indeed that in the presence of isolated cyclone and with CpH > 0.6, the SM BAZ point consistently toward the storm center.
Combining the azimuth issued from polarization analysis with the delays between the arrival time of PM and SM, it should be feasible to get an approximate location of the SM source, using a single seismic station e.g., [5]. However, polarization must be used with caution since, as explained above, the results may be biased by simultaneous SM sources from other storms worldwide, and in the case of La Réunion, by the strong depressions moving around Antarctica and generating austral swells [29]. Figure 7 indicates that a cyclone far from the seismic station (distance > 2000 km) has little to no influence on the computed BAZ (e.g., Flamboyan May 2018). Our polarization analyses show that a cyclone located in the Mozambique Channel (e.g., Idai, March 2019, Figure S5c) does not influence the recorded SM in La Réunion, which is fully consistent with the PSD results in Section 3.1.3 that showed weak amplitudes. Tracking the cyclone trajectory is beyond the scope of the present work, but should be achieved using different methods e.g., [9,[18][19][20][21].

Secondary Microseisms and Cyclone Intensity
In this Section, we analyze and discuss the seismic signatures of cyclones that passed close to Réunion Island, at a distance ≤500 km and recorded by the permanent (PF) and temporary (ZF) seismic networks on the Island. As discussed in Section 3.1, the intensification of the SM amplitude depends on the storm's intensity and the distance from the station to the storm center. Our results suggest that the seismic station RER starts to record SM visible on the PSD for a storm with V TC > 60 km/h, i.e., above a depression stage (Table S1) and at a distance as far as ∼4000 km (see Figures 1-3). Figure 8a maps cyclone tracks that passed in the proximity of La Réunion (distance ≤ 500 km) during the period 2011-2019 and Figure 8b-q present the variations of the amplitude of the SM during each cyclone, at all seismic stations running on the island at that time (cf map in the inset), together with the cyclone distance and intensity. Figure 8b-q show that during the life span of a given storm, all operating seismic stations display very homogeneous records. For a given cyclone, the temporal evolution of the SM amplitude (same as for PM amplitude, see Section 3.2.3) has a very similar pattern at all stations despite some differences in amplitudes, depending on the location of the station on the island. This observation suggests that the amplitude of the microseismic noise largely depends on the storm dynamics and is modulated by local sites effects. Figure 8 indicates that each cyclone has its own seismic signature and confirms that the maximum SM amplitude depends on the distance to the storm center (e.g., Cherono, Figure 8b) and/or on its strength (e.g., Dumile and Fakir in Figure 8d,o). In this case, the SM is likely to be generated at the storm center in the presence of broad ocean wave directional spectra, as modeled by Ardhuin et al. [8]. However, as seen in Section 3.1, a long-distance austral swell, if present, may influence the observed SM. A delay (a few hours to a few days) between the maximum cyclone intensity (and after the passage of the cyclone close to the island) and the generation of the SM is visible for some storms (e.g., Edilson, Berguitta, Figure 8h,l). We associate these SM to the possible interaction between the cyclonic sea-states and an austral swell or to the time needed for a cyclone to generate standing waves resulting from the storm itself and its displacement. Alternatively, this delay may be related to the slow wind-wave growth, which may take from a few hours to a few days [32], as mentioned previously.

Primary Microseisms and Significant Wave Height
After analyzing the microseismic noise in the SM band (i.e., in the 3-10 s period), we focus here on the PM band (i.e., in the 10-20 s period) which is accepted to be related to the interaction of incoming waves with the local coastal bathymetry. In this section, we therefore compare the seismic signal to the significant wave height (H S ) extracted from the global GOW2 hindcast wave model from Perez et al. [22] implemented with the WAVEWATCH III wave model [33] with a reduced 0.25 • grid space for the coastal locations and 0.5 • elsewhere. Their model parametrizations used a high resolution of the input forcings and of the output wave model allowing the cyclones to be well captured, hence motivating our choice for the GOW2 model.
The tracks of the cyclones used to investigate the PM are plotted in Figure 9a and the recording stations on La Réunion are shown on the inset map. We used the same cyclones and seismic stations as for the SM analysis in Section 3.2.2 ( Figure 8a). For each cyclone, the PM RMS amplitude measured at the various seismic stations are shown in Figure 9b-q (using the same color code as the cyclone tracks in Figure 9a). The corresponding H S during the life span of each cyclone, at different nodes are also presented as colored continuous lines. These figures suggest that the characteristic of the H S at each node may slightly differ in amplitude from each other and depends on the cyclone trajectories relative to the island. In the presence of a cyclone, the maximum H S is located at the nodes closest to the cyclone center. It is clear that cyclones tracked east of the island have a maximum H S computed at the eastern nodes of the wave model (e.g., Cherono, Dumile, Funani & Gelena, Figure 9b,d,q). Alternatively, when cyclones passed west from the island, the maximum H S are clearly recorded at the western nodes (e.g., Felleng, Bejisa, Dumazile, Figure 9e,f,m). In the presence of a cyclone tracked on both sides of the island (i.e., formed at the western side and that passed close to the island at the eastern side or viceversa), the nodes displaying maximum H S values change accordingly. For example, Fakir (Figure 9o) began on 20 April 2018 on the north-west of the island, inducing a maximum H S at the northern and western nodes. However, on April 24, Fakir was tracked at the eastern side of the island inducing maximums H S values at the eastern nodes. These observations confirm that the GOW2 model captured well the cyclone activity and its local signatures around Réunion Island, confirming its ability to study the cyclone activity around La Réunion and in particular the relations between the PM amplitude and the wave height H S .
For a given cyclone, the variations of the PM amplitudes recorded by the different seismic stations on the island display very similar patterns (Figure 9b-q). However, each station displays slightly different PM amplitudes, likely depending on local site effects and on the location of the station relative to the coast (similar conclusion as for the SM). Our results show that the strength of the PM amplitude also varies with the distance to the storm center and intensity. The cyclone trajectory also influences the PM amplitudes. Cyclones that passed west of the island (e.g., Dumile, Figure 9d) generated much stronger PM than those tracked east of the island (e.g., Fakir, Figure 9o). Dumile and Fakir have a comparable wind intensity (V TC = ∼125 km/h Figure 8d,o), but the PM amplitude generated by Dumile is twice stronger than Fakir (Figure 9d,o). The observed weak PM signal from Fakir is most likely due to the fact that the waves generated by the cyclone at the north-east of the island were only partly interacting with the local bathymetry and/or was rapidly attenuated before interacting with the local bathymetry.   At first sight, the PM RMS amplitude correlates well with the H S amplitude extracted from the GOW2 model (Figure 9b-q), in good agreement with the process of PM generated by oceanic waves interaction with the local bathymetry, creating elastic waves with the same frequency as the oceanic waves. Some exceptions are observed where the strongest PM is recorded a few hours (to a few days) after the strongest H S . Those PM are likely originating from long-distance swells. As an example, for Edilson (Figure 9h), the maximum PM amplitude (∼3.5 µm) is recorded one day and 15 h (from 7 to 9 February 2014) after it was tracked at the nearest of the island and the highest H S at different nodes were computed. However, a small PM peak (with a maximum value of ∼1.5 µm) is also visible in the measurements from 5 to 7 February and likely corresponds to the maximums H S (i.e., between 5 and 7 February 2014) generated in coastal areas. The recorded PM between 7 and 9 February was likely originating from long distant swells. The presence of the PM generated from simultaneous local waves and distant swells are seen for Giovanna, Felleng, Bejisa, Fernando (Figure 9c,e,f,k). Rindraharisaona et al. [5] have also associated the observed PM between 13 and 15 March 2017 for Fernando to a distant source.

Quantifying Relationship between Primary Microseisms and Significant Wave Heights
To determine the relationship between the PM amplitude and the significant H S , we present in Figure 10 two cyclones as examples: Dumile and Fakir. Figure 10a shows the tracks of cyclones in the neighborhood of La Réunion. Dumile formed on 27 December 2013, as a tropical disturbance (TD) at the north-east of the island and became a TC on 3 January 2014, during its passage west of the island and continued its journey southward. In contrast, Fakir started as disturbed weather on 20 April 2018, north-west from La Réunion and evolved rapidly and approached the island at its northern and eastern sides on April 24 as a TC. The following zoomed maps of La Réunion Figure 10b,c plot the Pearson correlation coefficients (P coe f ) between the PM amplitude observed at the seismic station MAT (located on the map) and the modeled H S at each surrounding node. We selected the station MAT because this permanent station operated since January 2011 and is characterized by a maximum PM amplitude (in most cases), suggesting a good sensitivity to microseisms. In general, for each cyclone, at least one node has a P coe f > 0.80. The P coe f for the different cyclones were computed using the data between the period within the vertical black dashed lines on each following subplot (Figure 10d,e). During these selected periods, the storm center was either at the north or alongside the island. This indicates that a simple linear regression is appropriate to estimate the relation between the near-coastal PM and the modeled H S . As mentioned previously, H S around the island strongly depends on the cyclone trajectory. The PM amplitude was simultaneously strengthened with the H S for Dumile and Fakir, as the wind speed increased to their maximum. Figure 10b shows that in the case of Dumile, the western nodes have higher P coe f than the northern or eastern nodes. The highest correlation (P coe f = 0.93) corresponds to the node W1 Table 2, i.e., at the node west of the island, close to the shore and close to the seismic station MAT, indicating that the oceanic waves interacted efficiently with the local bathymetry around this point (north-western side). In contrast, for Fakir (Figure 10c), H S modeled at the northern nodes and the node at 55.75 • E 21.00 • S correlate well with the PM amplitude, with P coe f > 0.93, which is consistent with the track of Fakir that passed north and east of the island. The significant wave heights modeled at the northern nodes continued travelling southward and likely interacted with the local bathymetry on the northern shores (hence the highest P coe f = 0.97 at the node N), which, in turn, generated PM that was recorded by all seismic stations on the island (Figure 9o). As Fakir continued its journey toward the south, H S in the northern nodes decreased rapidly whereas nodes on the eastern side displayed an increasing and a second H S peak on April 24 (for the nodes located in the eastern part of the island). We used these five nodes to compute the transfer functions presented in Figure 10f,g and Table 2.   To validate the transfer functions, we estimated H S for several cyclones that passed near La Réunion, and recorded by seismic station MAT (Figures S6-S8). In these figures are shown first (in a) the tracks of the cyclones considered and then the H S -PM Pearson correlation coefficients (P coe f ) at the wave model nodes located around the island for (b, c and d). The three bottom subplots (e, f and g) present the observed PM at station MAT, together with the estimated H S issued from the transfer function and the modeled GOW2. The vertical black dashed lines limit the period used to compute the P coe f . In general, the difference between the modeled and computed H S is less than 1 m, for any nodes that had P corr ≥ 0.8, whereas larger disagreement is seen between the two parameters with P coe f < 0.8. More details about the transfer function validation are given in Section S2.

Conclusions
We investigated microseisms generated by tropical cyclones/storms between 1999 and 2020 in the SWIO, using seismic stations in Réunion Island. The daily average of the PSD indicates the intensification of the microseisms in the presence of cyclone. Our analysis of the temporal evolution of the PSD of the SM, at periods ∼4 to 7 s, indicates a close relationship between the SM amplitude and the cyclone intensity. SM polarization analyses show that one can retrieve the BAZ of the storm center in the presence of well-polarized signals, i.e., CpH > 0.6 and isolated cyclone. The near-coastal PM amplitude (in the 10 to 20 s period) correlates well with the significant wave height from global hindcast model, with P coe f > 0.8. We have shown that PM represent a good proxy of the oceanic wave heights under cyclonic conditions. We thus computed transfer functions between the two parameters, allowing the derivation of modeled cyclonic significant wave heights from the seismic noise amplitude recorded on land, and to provide the new opportunity to use coastal and island seismic stations as terrestrial wave gauges. Continuous monitoring of microseisms provides therefore additional information on the cyclone and ocean activity, complementing satellite and other ground or marine observations.
Our results increase the number of the observational parameters related to past cyclones in the SWIO, which is crucial for better understanding the effects and signatures of cyclones (here the ocean activity) and for better forecasting the future cyclone activity and the related hazards. The present work demonstrates the feasibility of using microseismic noise in monitoring past oceanic activity. However, the non-linear coupling between atmosphere and ocean (hence the generation of microseismic noise) may limit the monitoring of the cyclone activity in real time, using microseismic data alone. Moreover, the simultaneous occurrence of the austral swell and/or two cyclones also complexifies the data analysis, in particular for determining its location. This issue could be solved by performing an f − k (frequency-wave number) analysis and/or other methods that were referenced in Section 3.2.1. Despite these limitations, the present work suggests that the use of the seismic data to quantify the sea state may bring a clear added value and even if there are still challenges and improvements that needs to be overcome beforehand, real-time monitoring of sea state from seismic data is undoubtedly a realistic opportunity.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/atmos12040488/s1. the PF seismic stations and for the easy access and the quality of their data. We appreciate the help of E. Delcher during the Riviere des Pluies fieldworks. We thank Laura Ermert and an anonymous reviewer for providing constructive comments and suggestions in their reviews and the Editor for additional careful reading of the manuscript. This is IPGP contribution number 4188.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: