Increasing ocean wave energy observed in Earth’s seismic wavefield since the late 20^th century

Domain: earth_science
Published: 2023-11-01
Editorial outcome: minor_revision

Abstract
Ocean waves excite continuous globally observable seismic signals. We use data from 52 globally distributed seismographs to analyze the vertical component primary microseism wave field at 14–20 s period between the late 1980s and August 2022. This signal is principally composed of Rayleigh waves generated by ocean wave seafloor tractions at less than several hundred meters depth, and is thus a proxy for near-coastal swell activity. Here we show that increasing seismic amplitudes at 3σ significance occur at 41 (79%) and negative trends occur at 3σ significance at eight (15%) sites. The greatest absolute increase occurs for the Antarctic Peninsula with respective acceleration amplitude and energy trends (±3σ) of 0.037 ± 0.008 nm s−2 y−1 (0.36 ± 0.08% y−1) and 4.16 ± 1.07 nm2 s−2 y−1 (0.58 ± 0.15% y−1), where percentage trends are relative to historical medians. The inferred global mean near-coastal ocean wave energy increase rate is 0.27 ± 0.03% y−1 for all data and is 0.35 ± 0.04% y−1 since 1 January 2000. Strongly correlated seismic amplitude station histories occur to beyond 50∘ of separation and show regional-to-global associations with El Niño and La Niña events.

Introduction
Earth’s seismic wave field was revealed by the late 19th century to be incessantly excited at periods between ~8 and 30 s. It was well established by the 1960s that seafloor forces due to wind-driven ocean gravity waves are the principal source of seismic waves in this period range in the absence of earthquakes or other large transient events, and that this microseismic wave field is primarily composed of seismic surface waves1. Mid-20th century studies established phenomenological understanding of the distinct ocean wave to solid Earth source coupling processes responsible for the primary (≈14–20 s) and secondary (≈6–12 s) microseism period bands. Microseism signals occur at much longer periods than typical anthropogenic seismic noise2 and thus constitute an ocean-state proxy that is well recorded by seismographs essentially everywhere on Earth.

Global microseism signals are generated by two distinct source processes. The primary microseism between ~14 and 20 s is excited by normal and shear tractions due to dynamic pressure variations of ocean swell on the near-coastal seafloor 3–5. Ocean waves in this period range have wavelengths λ in the range of 300–600 m, and their dynamic pressure field decays exponentially with depth at a rate determined by the wavenumber-depth product. The primary microseism spectrum6 matches the swell spectrum of partially to fully developed seas excited by worldwide storm systems, and microseism amplitudes reflect near-coastal seafloor tractions influenced by swell amplitude, propagation direction, and/or wave period. The distinct and more energetic secondary microseism between ~6 and 12 s, not analyzed in this study, has been investigated more thoroughly than the primary mechanism. It arises from wave period-shifted pressure perturbations on the ocean floor generated by generalized standing (clapotic) components of the ocean wavefield resulting from swell interference. These standing wave components may be generated by wave interactions occurring within a single translating storm system, between swell arising at multiple storm centers, or between incident and reflected waves due to coasts7. Notable efforts in the theory and numerical modeling of seismic ocean source processes in recent years have sought to more thoroughly incorporate necessary elements of seafloor bathymetry, physical oceanography, and seismology7–11 to refine the understanding of the ocean swell wavefield responsible for the secondary microseism and the influences of water depth, wave propagation direction, seafloor slope, and general bathymetry on primary and secondary microseism generation12,13.

Increasing ocean basin surface wind speeds since the mid-20th century have been inferred from meteorologic, oceanographic, and satellite altimetry data. The greatest wave powers and power increases are attributed to the Southern Ocean between 40∘ and 80∘ S14,15. Relevant factors vary geographically over multi-decade timescales under the influence of troposphere and ocean warming16–21 with 14 estimating that global wave power has on average increased by 1.087 × 103 kW m−1 y−1 and correlates strongly with sea surface temperature for 1948–2008.

The microseism wave field arises from geographically distributed forces applied to the seafloor by wind-driven ocean wave activity and is thus a proxy for the ocean wave state that complements surface and remote measurements6. This indicates that secular and other trends in ocean waves state are expected to be globally reflected in seismic data. Microseism-based wave state studies spanning shorter timescales than presented here have identified spatiotemporal trends in storm intensity, duration, and tracking22–26, and have been applied to estimate high-latitude sea ice variations and associated ocean wave attenuation27,28. Wave state is modulated by inter-annual climate processes14,29–31 and extreme microseism intensity associations, particularly with El Niño (ENSO) states, have previously been noted using shorter duration data sets25. However, global evidence for widespread secular intensification of microseism amplitudes has not been previously documented.

In this work, we assess primary microseism intensity since the late 20th century using data from globally distributed research and monitoring seismographs. We demonstrate a strong prevalence of increasing seismic amplitudes and energies at rates that are consistent with ocean wave intensification estimates from other disciplines. We additionally perform a correlation and clustering analysis of seismic amplitude time series to demonstrate that the primary microseism signal is a consistent proxy for long-range ocean wave spatial and temporal variability.

Results
Calculation of microseism metrics
We calculated microseism metrics from power spectral densities (PSDs) through 1 August 2022 using the section averaging method applied for the conterminous United States in ref. 32. The PSD estimation methodology is similar to that of ref. 33 but does not implement smoothing or binning beyond the inherent frequency discretization and spectral leakage properties of the discrete Fourier transform. We calculate seismic spectra in 1-hour 50% overlapping windows from continuous long-period high-gain vertical component (channel LHZ) 1 sample per second seismic data, which is insensitive to Rayleigh wave propagation direction, recorded by primary station sensors (location code 00) retrieved from the EarthScope Data Management Center (DMC) from the IU34 and II35 networks. Additional detail on spectral estimation is included in Methods.

Global earthquakes occur with a clustered Poissonian time distribution36 and large shallow events produce strong seismic surface waves that overlap with the 14–20 s primary microseism period band37,38. We remove potential earthquake transients from global earthquakes with magnitudes ≥ 5.75 using origin times and magnitudes from the U.S. Geological Survey Comprehensive Catalog of Earthquake Events and Products (ComCat)39 to cull intervals during which the primary microseism may be obscured. ComCat event selection is based on W phase40 moment magnitude Mw or on body wave magnitude mb for smaller earthquakes, and we remove data for all stations beginning at the earthquake origin time for 3, 6, 12, 24, and 48 hours for magnitude ranges 5.75 ≤ M ≤ 7.0, 7.0 ≤ M ≤ 7.5, 7.5 ≤ M ≤ 8.0, 8.0 ≤ M ≤ 9.0, and M ≥ 9, respectively. These data removal windows are highly conservative given that the global transit time for surface waves in this period range is about 3 hours and only very large earthquakes produce significant multi-circumglobal signals37.

Seismic data contain calibration pulses that produce values that greatly exceed naturally occurring levels, and may also contain periods of nonphysically low or high amplitude values corresponding to local physical disturbances or technical malfunctions. We remove these nonphysical outliers by comparing integrated signal power to the global New High-/Low-Noise models (NHNM/NLNM) of ref. 41, applying an acceptance corridor for the primary microseism period band between five times that is predicted by the (only rarely reached) low-noise model and 50 times that are predicted by the high-noise model.

Characterization of stationary seasonal variations
Seismic data commonly exhibit large annual variability in microseism intensity corresponding to the winter development of extratropical cyclonic storms, high-latitude formation of sea ice, and other seasonal factors24,27,28. We characterize stationary annual variations as a standard Fourier series representation H(t) in which the Fourier coefficients for the tropical year fundamental and its first three harmonic coefficients are estimated by projecting microseism amplitude A(t) or energy (calculated from the square of seismic velocity) E(t) time series onto orthogonal trigonometric basis functions at periods T0, T0/2, T0/3, and T0/4 to obtain corresponding coefficients ai and bi

H(t) = a1 sin(t/T0) + b1 cos(t/T0) + a2 sin(2t/T0) + b2 cos(2t/T0) + a3 sin(3t/T0) + b3 cos(3t/T0) + a4 sin(4t/T0) + b4 cos(4t/T0) (1)

where T0 = 365.242 days. Aperiodic sampling and data gap tolerant Lomb-Scargale PSD analysis42 of daily sampled median acceleration amplitude time series and examination of decaying coefficient (ai, bi) amplitudes confirmed that four Fourier terms were sufficient to characterize all significant spectral lines for our purposes.

Robust trend estimation
To assess the influences of removing ComCat catalog earthquake windows, noise model culling, and the presence or subtraction of station-specific stationary annual harmonic (equation (1)) functions, we applied robust (ℓ-norm minimizing; see Methods) trend estimation to progressively processed datasets for each station. The removal of short data intervals using (1) earthquake catalog and (2) noise model determined outlier procedures, as described above, produces trend estimates with lesser uncertainties but does not appreciably affect overall assessment of global seismic amplitude and energy trends nor the conclusions of this study (e.g., Supplementary Fig. 1).

Hourly time series as described above are smoothed using a two-month (61-day) moving median in daily steps to produce 1 sample per dayseries signals for ℓ-norm minimizing linear trend determination. We obtain trend estimates for signals within the primary microseism band (14–20 s) at global broadband seismic stations selected solely for operational histories exceeding 20 years and data completeness exceeding 80%. These selection criteria allow for the estimation of primary microseism trends at 52 globally distributed stations (Table 1) with earliest data ranging from 1988 to 1999 (Supplementary Fig. 2).

We fit seismic vertical acceleration amplitudes A(t) (Fig. 1) with ℓ-norm minimizing linear functions to estimate annualized rates of amplitude change RA (Fig. 2a, b; Table 2) for each station. Seismic wave energy is proportional to the square of velocity amplitude, and we estimate these time series E(t) (Supplementary Fig. 3) to obtain annualized rates of change, RE (Fig. 2c; Table 3). Corresponding proportional annual percentage change rates PA and PE are calculated by normalizing rates by corresponding station median amplitudes (Supplementary Fig. 4).

PA = RA / med(A(t)) × 100% (2)

PE = RE / med(E(t)) × 100%: (3)

For both acceleration A(t) and velocity squared E(t) trend results (seasonal harmonics signals removed), 41 stations (79%) show positive and 8 (15%) show negative slopes at 3σ significance (Figs. 2, 3; Tables 2, 3). Trend estimates obtained after subtraction of the seasonal harmonic signal (equation (1)) had lesser uncertainties reflecting the reduced signal mean average deviation from the linear function fit, but showed similar results to those obtained directly from A(t) and E(t) (Fig. 3, Supplementary Fig. 1).

Stations with significant positive amplitude and energy trends have wide geographically distribution and show regional magnitude correlations that are particularly well resolved for northeastern North America, western Europe (historically the most densely instrumented regions), and for sites in the extratropical southern hemisphere. Stations with significant negative trends are restricted to the northern and western Pacific Ocean regions and to two sites (TUC, Tucson, Arizona, and HKT, Hockley, Texas) in the southern United States. Seismic acceleration and velocity squared trends (RA, RE; Tables 2, 3) with associated 3σ uncertainties range from (0.037 ± 0.008 nm s−2 y−1; 4.157 ± 1.071 nm2 s−2 y−1) at PMSA (Palmer Station, Antarctic Peninsula) to (−0.018 ± 0.002 nm s−2 y−1; −1.271 ± 0.195 nm2 s−2 y−1) at HKT. Trends exhibit a modest (correlation coefficient of 0.213) proportionality in that stations with high historical median primary microseism amplitudes (Supplementary Fig. 4) tend to also exhibit greater amplitude and energy increases (Fig. 4), as inferred for wave energy increase in the Southern Ocean43.

The highest proportional rates PE of microseism energy increase are observed at eastern North America station SSPA (Standing Stone, Pennsylvania) and at central North America station RSSD (Black Hills, South Dakota) at 0.837 ± 0.090 and 0.810 ± 0.088% y−1, respectively (Table 3). The highest absolute rate of amplitude and energy increase is observed at PMSA (Palmer Station, Antarctica), but its corresponding proportional rate of energy increase (0.576 ± 0.148% y−1) is not as high as at other (e.g., eastern North America) stations (Table 3), reflecting very high historical median primary microseism levels and perhaps the influence of variable Antarctic sea ice28. Similarly relatively moderate proportionate rates of energy increase are found for other high-energy southern hemisphere stations (TAU; Hobart, Tasmania, and PTCN; Pitcairn Island) which are expected to have sensitivity to the state of the Southern Ocean (Fig. 4).

Methods
Spectral estimation
Hourly PSD estimates were calculated via Welch’s section averaging method using eleven 1024-s subwindows with 75% overlap. Additional stations with shorter contiguous recording intervals from the Federation of Digital Seismograph Networks (FDSN) networks IC54, G55, SR56, AS57, DW58, and CU59 were identically examined and display similar trends but are not included in this study due to shorter operational intervals. Each estimate requires at least 90% data completeness within each 1-hour window. We zero-pad any data gaps within the 1-hour segments after linearly detrending each contiguous data segment to avoid broadband discontinuity-induced spectral artifacts60. For each subwindow we apply a Hann taper and calculate their average, removing the instrument response to obtain physical units as parameterized by authoritative EarthScope Data Management Center metadata, and returning acceleration PSD estimates in dB relative to 1 (ms−2)2/Hz. We estimate microseism band square root integrated power by integrating the hour-long PSD estimates in microseism period bands using the trapezoidal rule and then obtain a time series of the band-limited root power in acceleration units (ms−2). Acceleration spectra are integrated in the frequency domain and squared to obtain seismic energy metrics. Integration to velocity to estimate the velocity squared energy proxy does not have a large effect on proportional amplitude trend statistics but makes metrics more sensitive to lower-frequency signals and produces a secondary effect on station trend ranking (Tables 2, 3).

Robust trend estimation
The robust parameter estimates in this study solve for intercept and trend (slope) parameters for the linear function that minimizes the outlier-resistant ℓ-norm metric

min ||d - Gm||1 = abs(di - (Gm)i)

where m is the 2-parameter model vector specifying the linear trend intercept and slope, G is an n by 2 system matrix consisting of a first n-length column where each element Gi,1 = 1 and a second n-length column Gi,2 = ti where ti is the time of the ith data point, and d is an n-element acceleration amplitude or seismic energy time series vector. We do not interpolate across data gaps. The slope and y-intercept parameters and their covariance matrix C are estimated using iteratively re-weighted least squares (IRLS)61 implemented via the robustfit function of MATLAB®. C is estimated from the convergent least-squares solution of the weighted IRLS equations.

Dendrogram calculation
Dendrogram associations between the microseism acceleration and climate index time series (Fig. 5) were calculated with MATLAB® from the correlation coefficients between the 54 detrended and demeaned time series shown in Fig. 6. The associated hierarchical cluster tree was generated using single-linkage agglomerative clustering implementing the Ward objective function62 in which each iterative linking step minimally increases the total within-cluster variance of the time series. The associated separation metric used in the MATLAB® linkage function is Dij = 1 - Cij where Cij is the correlation between time series i and j.

Calibration
The Global Seismographic Network (GSN) is jointly operated and maintained by the EarthScope Consortium under funding from the U.S. National Science Foundation and the U.S. Geological Survey. This study utilizes authoritative metadata retrieved for GSN stations from the EarthScope data management system to convert time series counts to physical units. The network incorporates instrumentation configurations that have been increasingly and asynchronously standardized over time. Seismic sensors and data loggers in the network are subject to acceptance testing63–65 for adherence to manufacturer’s specifications (typically ±1% deviation gain and ±5° in phase relative to nominal response). After installation, responses are subject to calibration and consistency tests that have included Earth tides66, observations of normal modes such as 0S0 from great earthquakes67, and continuous use and consistency checking of the network for the location and quantification of earthquake sources, seismic tomography, and other applications. The EarthScope DMC also continuously calculates and monitors frequency-dependent quality assurance metrics for GSN and other stations in its ongoing operations65. Time-dependent gain reductions68 in the early GSN noted in a systematic study of Mw > 6.5 earthquakes were traced69 to humidity-related corrosion in feedback electronics in a few percent of Streckheisen STS-1 very broadband sensors. These issues were documented and corrected in the field, and effects on station responses occurred at longer periods than those analyzed here. Anomalous limited (small percentage of total time series) contiguous time periods associated with apparent system malfunctions or incorrect metadata were observed at ten stations after inspection of overlapping hourly PSD integral time series and were excised (Supplementary Fig. 2; Supplementary Table 1.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript presents an analysis on primary microseisms amplitudes\nfor wave periods in the range 14 to 20 s, complemenfing previous work based on double-frequency\nmicroseisms. The parficular interest for primary microseism is the expected linear and stable (but\nfrequency dependent) relafionship between ocean wave spectra in shallow to intermediate water\ndepths and the primary microseism spectra. As a result, microseism data may be used to infer ocean\nwave properfies in coastal areas from the fime scale of a single storm to interannual varaibility and\nsecular changes.\nHere the variability (fig 1) is mostly analyzed in terms trends over a fime frame that spans at least 22\nyears.\nI have found the manuscript well wriften and very interesfing, and raising quite a few quesfions. I thus\nbelieve that this paper can be accpeted after a minor revision that should help address the following\nissues.\n1) It is not clear that the trends presented here are consistent with published records. One difficulty is\nthe quesfion of calibrafion of instruments: in situ buoys are infamous for containing spurious trends (e.g.\nCollins and Jensen 2022), while satellite alfimeter require careful calibrafion and do not discriminate\nbetween different periods (e.g. Dodet et al. 2022). Some of the most careful trend analysis was\nperformed by Timmermans (2020) et al. using satellite data, and I believe they are generally consistent\nwith the present results. However, nothing is said about the calibrafion of seismometers. What is the\nabsolute calibrafion uncertainty for, say STS-1 or STS-2 instruments?\n2) Trends computed over 20-40 years are most often dominated by interannual variability, as\ndemonstrated by Hochet et al. (2023) for the North Atlanfic. Probably this issue should be menfionned.\nIn fact, rather than the trends, the interannual pafterns (North Atlanfic Oscillafions, ENSO ...) are\nprobably much more interesfing and meaningful: the comparison to Sasaki was parficularly interesfing\nwhen discussing the North Pacific.\nMinor issues:\nFigure 1: fime ficks are oddly placed: 1987, 1993, 1998, 2004, 2009, 2015 ... this is not a regular fime\ninterval and it makes detailed comparison to other datasets very challenging (beyond looking at the\nfigure other users may want to actually download the fime series). The authors should fix this and clarify\nif the verfical lines are on January 1st of the year stated below.\nthere are \"???\" in references 19 and 57", "reviewer_2": "Review comment on Aster et al., submifted to Nature Communicafions\nThis study presents that the primary microseism amplitudes have increased within 20 years at global\nscale, using spectral amplitudes of Rayleigh waves at a period of 14–20 s observed in the verfical\ncomponent of the broadband seismometers. The stafions where increases of the primary microseism\namplitudes are observed are spafially correlated. The results are very interesfing. However, I do not\nunderstand how the authors defined quanfifies used in the text, and also what kind of mechanisms\nmainly affected the increases of the primary microseism amplitude. I hope that the following comments\nare helpful when the authors improve the manuscript.\nMajor comments\n1) Definifion of quanfifies\nI do not understand how the authors defined quanfifies used in the text and figures. For example, it is\nnecessary to define the following quanfifies using equafions, such as L1-norm minimizing linear trend\nand its covariance matrix, proporfional rates of microseism amplitude and energy, and absolute rate of\namplitude and energy. How did the authors use median values to esfimate proporfional rate of\nmicroseism energy (or amplitude)? Due to the lack of this informafion, I do not understand the meanings\nof the quanfifies. It is also nice if the authors explain with equafions how L1-norm is calculated in the\nrobusffit funcfion of MATLAB.\n2) Mechanisms for the increases of primary microseism amplitudes\nAlthough the quanfifies used in this study should be defined, it seems that the increases of the primary\nmicroseism amplitudes are interesfing. However, there are less descripfions on the plausible mechanism\nthat increased the primary microseism amplitude. I would like to read descripfions on the mechanisms\nwith some quanfitafive esfimafions. Without such descripfions, it seems that the manuscript is weak for\npublicafion in Nature Communicafions.\nThe boftom paragraph in Page 6: This descripfion in very interesfing. However, when I read this\nmanuscript, I did not understand that the mechanisms stated here really change the primary microseism\namplitudes and explain the observed amplitude difference over 20 years. Therefore, I strongly\nrecommend the authors to perform some quanfitafive esfimafions. It is a good way if the authors\nconduct a three-dimensional numerical simulafion, maybe at a local scale, in which the Rayleigh waves\nare generated by the forcing due to ocean swells near coastal areas. If the sea level change or other\nfactors change the generated Rayleigh waves and they are matched with the observed quanfifies, it is\nvery interesfing. Other quanfitafive esfimafions with some theorefical approaches are also nice.\n3) Sentence structure\nIt is very hard to read the manuscript, and the sentence structure should be improved thoroughly. For\nexample, the descripfion in Secfion 1 (1.1 and 1.2) is related to method. The descripfion in Secfion 2\nshould be in the result secfion. The secfion 3 is too long for conclusion.\nOther comments.\nThe 10th lines from the boftom in Page 2\near-coastal --> near coastal\nThe 8th line from the boftom in Page 2\n“A second and more energefic”\nWhat is “second”?\nThe 9th line from the top in Page 3\n“Meteorological and climate reanalyses and modeling.”\nPlease rewrite this sentence.\nThe 14-11th lines from the boftom in Page 4\n“We calculate the least-squares...alter general conclusions (Fig. 3)”\nThis descripfion should be explained using equafions.\nAlso, Fig. 3 is cited before Figs. 1 and 2. Please rearrange the numbering the figures.\nFig 1 capfion\napplying thee-year --> applying three-year\nFigs. A5–A8\nWhy are the velocity amplitudes needed? In the spectra, velocity and accelerafion are linked with i-\nomega, so the tendencies of the amplitude variafions between the velocity and accelerafion should be\nthe same. It seems that A5–A8 are not needed.", "reviewer_3": "This manuscript examined long-term temporal variafions in microseism energy as an index of near-\ncoastal ocean wave energy and showed a mulfidecadal increase in microseism and ocean wave energy.\nThe authors used globally distributed seismic stafions with long operafional periods and showed clear\nsecular increases of microseism energy at most stafions and their regionally coherent paftern of the\ntrends. The data are rich, the analysis is generally rigorous, and the results contain interesfing findings.\nAlthough I do not hesitate to recommend the publicafion of this work, I feel the presentafion can be\nimproved, and some addifional figures and analyses can enhance the importance of this work and\nsupport the arguments in this manuscript.\nMajor comments:\n1. The differences between the authors' previous work [33] and this study are not very clearly described.\nThis study used the median values as the index of microseisms, which is different from the index\n(number of extreme microseism events) used by the previous work [33]. I think it would be befter to\ndescribe how the microseism index in each study reflects the state of ocean waves.\n2. The authors argued that “Trends exhibit a high degree of proporfionality in that stafions with higher\nmedian primary microseism amplitude strongly tend to exhibit greater amplitude and energy increases\n(Fig. 4)” in Paragraph 2 in Secfion 2. Although this proporfionality is one of the interesfing points\npresented by this work, I do not think that Fig. 4 provides direct evidence of the proporfionality between\nmedian amplitude and amplitude increasing rafio. This is because a nearly constant median amplitude\ncan produce the trend (proporfionality between percent change vs absolute change) shown in Fig. 4. To\nshow the proporfionality directly, the authors can show a plot taking the median amplitude in the\nhorizontal axis and the amplitude increasing rafio in the verfical axis. A map showing the median\namplitude is also informafive.\n3. Quanfitafive comparison with meteorological and oceanographic studies is limited. It would be useful\nto over-plot fime series of long-term wind speed data [20] and (regionally averaged) seismic data of this\nstudy and examine correlafion coefficients for the fime periods over which both data are available. Their\nregional differences also would be interesfing. Such quanfitafive comparison may enhance the\nusefulness and importance of microseism-based ocean wave state study including this paper.\nMinor comments:\n4. Abstract: I do not completely agree with “The primary microseism signal between 20–14 s period is\nprincipally composed of seismic Rayleigh waves”. The primary microseisms consist of not only Rayleigh\nwaves but also Love waves, and Love wave energy is as large as Rayleigh wave energy (Nishida et al.\n2008, doi:10.1029/2008GL034753).\n5. Abstract: It is hard to understand the difference between “the most rapidly proporfionately\nintensifying energy” and “the most rapidly intensifying primary microseism levels on an absolute basis”.\nThis may be due to the lack of the definifion of “the proporfionately intensifying energy” in the Abstract.\nI suggest that the authors add a descripfion of it in the Abstract.\n6. Paragraph 1, Introducfion: “ear-coastal sea floor” may be a typo.\n7. Paragraph 2, Introducfion: “Meteorological and climate reanalyses and modeling.” does not form a\nsentence.\n8. Paragraph 1, Secfion 1.2: The authors esfimated and removed annual periodic trends by fifting a\nperiodic funcfion to the original fime series with the secular trends and then esfimated the secular trend\nfrom the residual fime series. Strictly speaking, however, it is more accurate to fit annual and secular\ntrends simultaneously because there is no modeling error.\n9. Fig. 1: The spacing of the horizontal fick labels is not constant."}, "author_response": "Reviewer #1 (Remarks to the Author):\n\nWe appreciate the comments regarding buoy data and the referrals to additional\nreferences, which we have added in context to the Discussion. We agree\nregarding the uniqueness of period-dependent information in these data, and this\nis indeed a promising target for further study.\n\nWe added a GSN calibration subsection describing tolerances, calibration\npractices, and relevant references to the Methods section. The data and\nmetadata receive continuous attention from the seismological community for\nearthquake and other studies (because these are trend measurements, and much\nof the analysis uses historical median normalized data the absolute calibrations\nare not as important as the stability of station responses here).\n\nWe also plotted the Reguero et al [20] global trend in our (new) Fig. 7 to provide\ngraphical comparison to our energy trend estimates (please see points in (2)\nbelow). We note consistent El Nino/La Nina trends (Fig. 7a,b) with those noted by\nReguero, showing that those (mid-ocean) trends area also clearly reflected in the\nnear-coastal energy environment for these later data as indicated by primary\nmicroseism records from the western and southeastern Pacific regions,\nrespectively.\n\n2) Trends computed over 20-40 years are most often dominated by interannual\nvariability, as demonstrated by Hochet et al. (2023) for the North Atlantic. Probably this\nissue should be mentioned. In fact, rather than the trends, the interannual patterns\n(North Atlantic Oscillations, ENSO ...) are probably much more interesting and\nmeaningful: the comparison to Sasaki was particularly interesting when discussing the\nNorth Pacific.\n\nThis set of points (and the above reference) are now amplified in the Discussion\nand revised figures (see new Figs. 6 and 7, below). The revised manuscript also\nin this regard incorporates a correlation-based station history clustering (a,\nbelow) in response to this comment and to comments from other reviewers.\nAmong other points this clarifies ENSO and SOI index correlations with\nSouthwestern Pacific (cluster SWP) and Southeastern Pacific/Southwestern\nAtlantic (cluster SEPSWA) regional microseism intensities, respectively. These\nnew figures also demonstrate the imprint of ENSO/SOI excursions in Earth’s\nglobal median microseism wavefield (El Nino and La Nina events A – E noted in\nFig. 7b). The SWP and SEPSWA associations corroborate with the pre-2009\nanalysis of Reguero et al. (2019; their Fig. 7b) who noted correlations of global\nwave power with the Niño3 standardized index in which increased ocean wave\nenergy occurs in the western equatorial Pacific and southeastern Pacific during\npositive and negative excursions, respectively.\n\nMinor issues:\n\nThis (date roundoff problem in the plotting program) has been corrected, a\ncorrected figure was uploaded, and the tick mark timing was clarified in the\ncaption as suggested above.\n\nthere are \"???” in references 19 and 57\n\nThis has been corrected (the missing publisher addresses were added in the\nbibliography file).\n\n__________________________________________________________________\nReviewer #2 (Remarks to the Author):\n\nThank you for these comments. The revised manuscript clarifies our\nmethodology and the quantities calculated in the paper, as well as including\namplified discussion regarding the details of the primary microseism source\nmechanism.\n\nMajor comments\n\n1) Definition of quantities\n\nWe added (Methods) information on the robust trend fitting theory and\nmethodology employed here, including a reference in which the iterative re-\nweighted least squares (IRLS) procedure is described. We have also noted how\nthe covariance matrix is estimated from the reweighted system of least-squares\nequations at the convergence of this iterative algorithm. Proportional trends\nwere calculated by normalizing the trends by the historical median values\nobserved at each station (equations 2 and 3), and we incorporating a naming\nconvention for the various time series referred to in the manuscript to make these\nquantities clearer for readers.\n\n2) Mechanisms for the increases of primary microseism amplitudes\n\nThank you for inviting further elaboration on this central topic. Our principal aim\nin this (shorter format) paper is to report and quantitatively document the\nincreasing energy in Earth’s background seismic and causative ocean wavefield\nin a temporally and geographically extensive observational study using robust\nestimation. As we describe in the expanded Discussion, the proposed\nmechanism for widespread and systematic increase in primary microseism\namplitudes (and the only one that we ascribe to be credible) is a systematic\nincrease in the amplitude of the causative tractions on the seafloor at periods\nbetween about 20 and 14 s in at depths of less than several hundred meters (we\nadded the relevant equations 4 and 5 in this). The primary microseism\nmechanism is sufficiently well understood to contend that these time series are\nrecords of spatiotemporally integrated wave-induced ocean bottom forces across\nbroad regions of the Earth (e.g., Fig. 7c). In quantitatively associating seismic\nenergy increase with near-coastal wave energy increase, we appeal to linearity as\nelaborated in the Results section. We have additionally augmented the relevant\nreferences and elaborated more clearly outstanding challenges associated with\nfull modeling of primary microseism.\n\nWe appreciate this comment and respectfully assert that the 3-D modeling\nrecommendation in its fullest form extends beyond the scope of this paper, which\nprovides solid motivation for further advancement in this as a well warranted and\naspirational goal. At present this is a challenging undertaking even in lower-\ndimensional simulations. For example, Gualtieri et al. (2019), using the\ntheoretical approach of Hasselmann (1963) and Ardhuin et al. (2015, 2018) were\nable to demonstrate general agreement between modeling and observations for\nthe primary mechanism using a (linearized) 1-D Fourier method implementing\nwave and bathymetric spectra in combination with WAVEWATCH III estimates,\nbut noted further complexities that need focused attention including the\nnonlinear influences of steeper slopes, partitioning of energy into multimodal\nsurface and love waves, and 3-D (including wave directionality and small-scale\nbathymetry) coupling effects that necessitated the introduction of an empirical\nfitting parameter in their simulations. Empirically, this clarifies the broad source\nregions at play and highlights the extensive spatial integration involved\n(consistent with that noted by Gualtieri, 2019) in its correlations between El Nino\nand SW Pacific intensity, and between La Nina/SOI and SE Pacific/SW Atlantic\nintensity (Fig. 7a), as well as in the extensive spatial correlation lengths between\nstation time series depicted in Fig. 7c (which confirms the findings of Gualtieri\n(2019) the primary microseism signal at a particular station can be sensitive to\nsections of coast at ranges of up to 1000’s of km).\n\nMore broadly, our principal conclusion that the causative effect of the observed\ntrends is an increase in these seafloor forces over the past 3--4 decades due to\nincreasing wave energy performing elastic work on the seafloor is, we contend,\nvery well supported in this paper and in associated references. Because the\nchanges are (proportionately) small relative to the median or average excitation\nlevels, we invoke the reasonable assumption of quasi-linearity in our analysis to\nassociate seismic and wave energy rates of increase (as described in the revised\nmanuscript).\n\nEquation 5 in the revised manuscript notes the dynamic pressure period/depth\nrelationship for linearized ocean swell to illuminate the essential relationship\nbetween wave height and the (ocean floor) seismic source for the primary\nmicroseism. The reviewer’s perhaps broader point regarding sea level changes (if\nwe are interpreting this correctly) possibly affecting this signal is an intriguing\none that we previously investigated by probing our time series for evidence of\ntidal period modulations in primary microseism amplitudes using Lomb-Scargale\n(aperiodic sampling) spectral estimation. However, this only resolved the annual\nand annual harmonic periodicities in microseism amplitude and intensity\nassociated with the seasonal cycle (the large H(t) annual signal components\ndescribed in the manuscript). Since the process is excited by waves at ocean\nswell periods (with wavelengths of ~300 – 600 m) and the sensitivity footprint is\nwidespread (Fig. 7c) tidal-scale sea level changes do not appear so far to have a\ndiscernable signature.\n\n3) Sentence structure\n\nThe manuscript has been extensively rewritten reformatted in section\norganization, scientific units, and otherwise for Nature Communications (per\nhttps://www.nature.com/ncomms/submit/article).\n\nWe have also edited the manuscript throughout and received confidential internal\nreviews from colleagues within our institutions to assure general textural and\nmathematical readability, linkages with figures, figure ordering, and otherwise\naddress the reviewer’s general readability comment in mind. We explicitly\nspecified and cross-referenced relevant equations and implemented a time series\nnaming convention throughout the manuscript to further improve its accessibility\nfor a wide audience.\n\nOther comments.\n\nThe 10th lines from the bottom in Page 2\near-coastal --> near coastal\nCorrected\n\nThe 8th line from the bottom in Page 2\n“A second and more energetic”\nWhat is “second”?\nClarified -- this was in reference to the second (secondary) microseism source\nmechanism.\n\nThe 9th line from the top in Page 3\n“Meteorological and climate reanalyses and modeling.”\nPlease rewrite this sentence.\nThis sentence fragment has been removed.\n\nThe 14-11th lines from the bottom in Page 4\n“We calculate the least-squares...alter general conclusions (Fig. 3)”\nThis description should be explained using equations.\nDone. We also introduced a naming convention, with the functions A(t), E(t)\ncorresponding to the key time series, and the stationary seasonal harmonic\nprojections H(t). The relevant equation for H(t) is noted (equation 1).\n\nAlso, Fig. 3 is cited before Figs. 1 and 2. Please rearrange the numbering the figures.\nWe re-checked and corrected the figure order in the revised manuscript\n\nFig 1 caption\napplying thee-year --> applying three-year\nCorrected.\n\nFigs. A5–A8\nWhy are the velocity amplitudes needed? In the spectra, velocity and acceleration are\nlinked with i-omega, so the tendencies of the amplitude variations between the velocity\nand acceleration should be the same. It seems that A5–A8 are not needed.\nThese signals have a modest (~6 s) bandwidth. Thus, conversion from\nacceleration to velocity has some effect on corresponding metrics derived from\nthe integrated PSDs (e.g., weighting longer periods as 1/omega in the spectra).\nThe spectra are fairly narrow-band (6-s) and smooth (varying by a few dB) and the\nreviewer is correct in suspecting that relevant differences between these\nintegrals are modest in the context of this study (e.g., compare Figs. 1 and Supp.\nFig. 10 in the revised manuscript, the global acceleration history and global\nenergy history in Fig. 7a,b, or the spectra for GEOSCOPE network stations in\nGualtieri (2019). We agree that the seismic velocity amplitude plots and table in\nthe original manuscript were superfluous and they have been removed along with\nassociated references to velocity series throughout the paper.\n\n__________________________________________________________________\nReviewer #3 (Remarks to the Author):\n\nMajor comments:\n\n1. The differences between the authors' previous work [33] and this study are not very\nclearly described. This study used the median values as the index of microseisms,\nwhich is different from the index (number of extreme microseism events) used by the\nprevious work [33]. I think it would be better to describe how the microseism index in\neach study reflects the state of ocean waves.\n\nWe clarified the relationship between the previous work and these studies within\nthe Introduction section of the revised manuscript. The earlier (Aster et al.) index\nthat the reviewer refers to examined extreme outlier statistics of seismic data\n(e.g., extreme wave events in the record observed over short [a few hours] time\nintervals) while this study presents observations of characteristic (median) trends\nin amplitude and energy over a more recent and well-calibrated most recent 3-4\ndecades of the global seismic record, smoothed over 2-month and 3-year time\nscales. IN addition, the data examined here from the last 15 years or so were not\navailable at the time of publication for the previous study and show much more\nrobust evidence of global secular change than the shorter time series. We agree\nwith the reviewer that the data and methodologies in this paper indeed reveal a\nrich trove of multi-scale spatiotemporal statistical information relevant to\nevolving ocean wave state and hope that the study will inspire future work.\n\n2. The authors argued that “Trends exhibit a high degree of proportionality in that\nstations with higher median primary microseism amplitude strongly tend to exhibit\ngreater amplitude and energy increases (Fig. 4)” in Paragraph 2 in Section 2. Although\nthis proportionality is one of the interesting points presented by this work, I do not think\nthat Fig. 4 provides direct evidence of the proportionality between median amplitude\nand amplitude increasing ratio. This is because a nearly constant median amplitude can\nproduce the trend (proportionality between percent change vs absolute change) shown in Fig. 4. To\nshow the proportionality directly, the authors can show a plot taking the median amplitude in the\nhorizontal axis and the amplitude increasing ratio in the vertical axis.\n\nThank you for this comment -- we agree that the inherent correlation in plotting\ntwo dependent variables against each other was confusing in this figure. As\nsuggested by the reviewer, we replaced the original figure with a straightforward\nand we suggest more informative plot of absolute trend versus median\nproportionality as suggested (new Fig. 4; below). This plot reinforces the\ndominance of positive trends, shows where various stations lie (with the various\ndotted-line linear trends) and further displays trend variations with general\nregion. We have also modified the text accordingly in the discussion of this\nproportionality effect, which is evident but modest for some of the regional arms\nas indicated by color (correlation of 0.213) with proportionality being most\nevident for the North American and European Atlantic regions (red, blue,\nrespectively ) and South Pacific (green and teal) regions/stations.\n\nA map showing the median amplitude is also informative.\nWe added a map figure of median amplitudes to the manuscript (Supp. Fig. 11;\nbelow).\n\n3. Quantitative comparison with meteorological and oceanographic studies is limited. It\nwould be useful to over-plot time series of long-term wind speed data [20] and\n(regionally averaged) seismic data of this study and examine correlation coefficients for\nthe time periods over which both data are available. Their regional differences also\nwould be interesting. Such quantitative comparison may enhance the usefulness and\nimportance of microseism-based ocean wave state study including this paper.\n\nWe significantly fleshed out these critical points of connection in the revised\nmanuscript and graphically display these associations in time series format in the\nnew Figs. 5 and 7 (below). Specifically, the revised manuscript incorporates a\ndendrogram analysis of station histories that notes the revealed geographically\nassociations (Fig 5b). We associated time series correlations with the Southern\nOscillation and Bivariate ENSO indices, which are nearly anticorrelated (\nhttps://psl.noaa.gov/data/climateindices/list) and note ENSO influence in the\nregional and global data.\n\nThe Reguero et al. (initial manuscript reference [20]) analysis is a solid and well-\nknown touchstone for independent remote sensing estimates of global wave\nenergy but ends in 2008 (e.g., their Figs. 2 and 3) and thus does not overlap with\nnotable regional and global features shown in the (new) Fig. 7. We have included\nadditional discussion regarding regional and global sensitivity in the seismic data\nto ENSO excursions.\n\nFig. 7 in the revised manuscript now compares both our linear trends from this\nstudy with the overall (1948-2008) trend from Reguero et al. [20]. Figs. 5 and 7\nindicate that southwestern Pacific stations (WRAB, TAU, PMG, CTAO, GUMO, and\nCHTO; cluster SWP) cluster with ENSO (~ -SOI), and Eastern Pacific/South\nAtlantic/Antarctic stations (EFI, LCO, NNA, PAYG, PMSA; cluster SEPSWA)\ncorrelate and cluster with SOI (~ the additive inverse of ENSO) across the multi-\ndecadal observational period. This result is consistent as noted in the new\nmanuscript with the conclusions displayed in and drawn from Fig. 7 by the\nReguero et al. for 1948—2010. We tested for correlations with the other indices\nretrieved from the above NOAA site and did not identify significant associations\nat this time scale and temporal smoothing window.\n\nCorrelation analysis included in this revised manuscript further clarified the long-\nrange extent of primary microseism geographic averaging, with interstation\ncorrelations between time series showing significance to station separations\nbeyond 50o (Fig. 7 c) thus further demonstrating the extensive geographic\naveraging within the primary microseism signal. Importantly, this shows that the\nprimary microseism signal as processed here (Figs. 6, 7a) via the use of 2-month\ntemporal smoothing and normalization by local historical median levels provides\na consistent proxy for long-range ocean wave intensity across the substantial\nsectors of the Earth identified in our clustering algorithm and as evidenced by\ngeneral time series consistency from station to station across the clustered\nregions identified in this study and supported by the correlative associations with\nthe ENSO and SOI indices.\n\nMinor comments:\n\n4. Abstract: I do not completely agree with “The primary microseism signal between 20–\n14 s period is principally composed of seismic Rayleigh waves”. The primary\nmicroseisms consist of not only Rayleigh waves but also Love waves, and Love wave\nenergy is as large as Rayleigh wave energy (Nishida et al. 2008,\ndoi:10.1029/2008GL034753).\n\nThe reviewer’s point is correct -- we clarified in the abstract and text that this\nstudy uses vertical-component data that are insensitive to Love waves (and to\nRayleigh wave propagation direction as well) and are dominated by fundamental\nmode Rayleigh waves (e.g., Schimmel et al., 2011; DOI: 10.1029/2011GC003661).\nWe also added an introductory sentence to better recognize the foundational\nNishida et al. and other papers that address the multiplicity wave types within the\nmicroseism wavefield.\n\n5. Abstract: It is hard to understand the difference between “the most rapidly\nproportionately intensifying energy” and “the most rapidly intensifying primary\nmicroseism levels on an absolute basis”. This may be due to the lack of the definition of\n“the proportionately intensifying energy” in the Abstract. I suggest that the authors add a\ndescription of it in the Abstract.\n\nWe have clarified in the abstract and elsewhere that reported proportionate\ntrends are relative to long-term station median values and added relevant\nequations and symbolic definitions (equations 2 and 3) within the text, including\nnoting explicitly in the abstract that reported percentage trends at each station\nare relative to long-term station median values.\n\n6. Paragraph 1, Introduction: “ear-coastal sea floor” may be a typo.\nCorrected.\n\n7. Paragraph 2, Introduction: “Meteorological and climate reanalyses and modeling.”\ndoes not form a sentence.\nThis sentence fragment has been deleted.\n\n8. Paragraph 1, Section 1.2: The authors estimated and removed annual periodic trends\nby fitting a periodic function to the original time series with the secular trends and then\nestimated the secular trend from the residual time series. Strictly speaking, however, it\nis more accurate to fit annual and secular trends simultaneously because there is no\nmodeling error.\n\nBecause the trend and y-intercept and L-1 norm fitting functions are orthogonal\nto the trigonometric basis functions that characterize stationary annual and\nharmonic trends the two approaches yield effectively identical results. We also\nverified that four Fourier terms were sufficient to remove stationary seasonality\nby examining the complex coefficient amplitudes for each stations, The complex\nharmonic coefficient amplitudes for acceleration amplitude data sorted by\nlatitude are shown in the (not-included in manuscript) figure below, which also\nshows the N-S hemispherical phasing and other interesting details of this signal.\nWe also examined spectral lines within the Lomb-Scargale spectrum of the daily-\nsampled multidecadal time series to search for any additional periodicities (e.g.,\ntidal) but these were not apparent.\n\n9. Fig. 1: The spacing of the horizontal tick labels is not constant.\nThis (labeling numerical roundoff) issue has been corrected in the new Figs. 1\nand Supp. 10."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have made a thorough effort to consider all previous comments, with some interesfing new\nfigures and updated text, including now a menfion of interannual variability in the abstract. Also it is very\nimportant that defined in detail the quanfifies that are used.\nI believe the paper can be accepted as is.\nObviously the new figures lead to a few more quesfions (like: why is the KIV curve behaving the way it is\n? Could there be an issue with instrument response there? ... I've delayed my review because I wanted to\ncheck on this, using microseism spectra computed by IPGP - these are available as part of a European\nSpace Agency project on climate fime series - but could not find the fime)\nThank you again for an in depth revision.", "reviewer_2": "Comments to Aster et al., submifted to Nature Communicafions\nI reviewed the previous version of the manuscript, and I am almost safisfied with the current version of\nthe manuscript. The authors did a great job. Thank you very much. It is nice that the authors address the\nfollowing points before publicafion of this study.\nMain points\n7th-8th lines in the second paragraph in Page 4\nthat predicted --> that are predicted\nPage 7\nThe authors should add more explanafions for making the dendrogram in Fig. 5 using the correlafion\ncoefficient in Fig. 14. It is necessary to explain how the clustering is performed with cifing some\nreferences.\n7th line from the boftom in the second paragraph in Page 7\nsubsequent excursions, perhaps due to lower data\nsubsequent excursions, is perhaps due to lower data\nPlease check this part again.\n3rd line in the third paragraph in Page 7\na slope change near 50º\na slope change near a distance of 50º\nCorrected.\n5th line from the top in Page 8\nPlease finish the final part of “stafion within the“.\n4th line in the second paragraph in Page 9\nthis signals reflects\nthese signals reflect\nPlease check this part again.\nEquafion (6)\nPlease define the column and row sizes of matrices of m and G, and their elements precisely. What are\n“ones” in “a column of ones”?\nEquafion (4) and the capfion of Fig. 1\nBoth contents include the character of lambda (wavelength and longitude), which confused readers. It is\nbefter to change either.\nCapfion of Fig. 5\nThe two sentences after “Regions of ocean bathymetry with depths” should be moved for the\nexplanafions for panel (a).", "reviewer_3": "The authors have adequately addressed all comments on the previous manuscript, and the manuscript\nhas been greatly improved. The correlafion between primary microseisms and the ENSO indices, added\nin the revision, is fascinafing. I believe this research will have a major impact on a wide range of research\nfields. I would like to add just minor comments, which can be addressed by the authors before the\npublicafion of the manuscript.\nMinor comments:\n1. Fig. 7b indicates that both posifive and negafive ENSO phases tend to increase the global average\nenergy of primary microseisms (may be due to teleconnecfion). However, the SWP and SEPWA regions\nshow correlafions more directly, such as microseism energy increases in posifive ENSO phases and\ndecreases in negafive ENSO phases. The difference may be confusing. It would be befter to clarify the\ndifference in the near and distant influences of ENSO. The authors can display the average microseism\nfime series for the area excluding SWP and SEPWA in addifion to the global average and discuss the near\nand distance ENSO influences.\n2. The 23th line from the top in Page 3: … has not been been …\n3. The 4th line from the top in Page: The reference for the ENSO and SOI indices is necessary here.\n4. The 4th line from the top in Page8: stafion within the ?\n5. Figure 4: The unit of the horizontal axis should be nm/s^2.\n6. Figure 6: It is very hard to see the correlafion between microseisms fime series and ENSO/SOI indices.\nImprovements are required, such as stretching each panel horizontally."}, "author_response": "Reply to REVIEWERS' COMMENTS (rev 2)\nIncreasing near-coastal ocean wave energy observed in Earth's seismic\nwavefield since the late 20th century\nR. Aster, R. Ringler, R. Anthony\nPlease note that we have interspersed the authors’ responses in bold\nhere.\n\nReviewer #1 (Remarks to the Author):\n\nThe authors have made a thorough effort to consider all previous comments,\nwith some interesting new figures and updated text, including now a mention\nof interannual variability in the abstract. Also it is very important that defined in\ndetail the quantities that are used.\nI believe the paper can be accepted as is.\nObviously the new figures lead to a few more questions (like: why is the KIV\ncurve behaving the way it is? Could there be an issue with instrument\nresponse there? ... I've delayed my review because I wanted to check on this,\nusing microseism spectra computed by IPGP - these are available as part of a\nEuropean Space Agency project on climate time series - but could not find the\ntime)\nThank you for your thorough re-review of the paper. There are indeed\nmany other features of these types of data that warrant additional study,\nand it’s good that the ESA groups are analyzing microseism spectra in\ntheir efforts as well.\n\nReviewer #2 (Remarks to the Author):\n\nComments to Aster et al., submitted to Nature Communications\nI reviewed the previous version of the manuscript, and I am almost satisfied\nwith the current version of the manuscript. The authors did a great job. Thank\ny\nThank you very much. It is nice that the authors address the following points before\npublication of this study.\n\nMain points\n\n7th-8th lines in the second paragraph in Page 4\nthat predicted --> that are predicted\nCorrected.\n\nPage 7\nThe authors should add more explanations for making the dendrogram in Fig.\n5 using the correlation coefficient in Fig. 14. It is necessary to explain how the\nclustering is performed with citing some references.\nWe added a Methods subsection describing how the dendrogram was\ngenerated that includes appropriate seminal reference (Ward, 1963) and\nspecific mention of the MATLAB linkage function and dissimilarity\nmeasure used.\n\n7th line from the bottom in the second paragraph in Page 7\nsubsequent excursions, perhaps due to lower data\nsubsequent excursions, is perhaps due to lower data\nPlease check this part again.\nThis point has been split into two sentences to improve readability.\n\n3rd line in the third paragraph in Page 7\na slope change near 50º\na slope change near a distance of 50º\nCorrected.\n\n5th line from the top in Page 8\nPlease finish the final part of “station within the“.\nThis sentence fragment has been corrected. The associated sentence\nnow reads:\nThis correlation length scale indicates that primary microseism\namplitudes and energies provide consistent proxies for large spatial and\ntemporal scale ocean wave variability.\n\n4th line in the second paragraph in Page 9\nthis signals reflects\nthese signals reflect\nPlease check this part again.\nCorrected to: … indicate that these signals reflect…\n\nEquation (6)\nPlease define the column and row sizes of matrices of m and G, and their\nelements precisely. What are “ones” in “a column of ones”?\nThese points are now clarified in the associated Methods subsection\n(“ones” means all of the elements in the column are equal to 1).\n\nEquation (4) and the caption of Fig. 1\nBoth contents include the character of lambda (wavelength and longitude),\nwhich confused readers. It is better to change either.\nWe revised the caption to use the script lower case “l” symbol to\nindicate longitude.\n\nCaption of Fig. 5\nThe two sentences after “Regions of ocean bathymetry with depths” should be\nmoved for the explanations for panel (a).\nDone.\n\nReviewer #3 (Remarks to the Author):\n\nThe authors have adequately addressed all comments on the previous\nmanuscript, and the manuscript has been greatly improved. The correlation\nbetween primary microseisms and the ENSO indices, added in the revision, is\nfascinating. I believe this research will have a major impact on a wide range of\nresearch fields. I would like to add just minor comments, which can be\naddressed by the authors before the publication of the manuscript.\n\nMinor comments:\n\n1. Fig. 7b indicates that both positive and negative ENSO phases tend to\nincrease the global average energy of primary microseisms (may be due to\nteleconnection). However, the SWP and SEPWA regions show correlations\nmore directly, such as microseism energy increases in positive ENSO phases\nand decreases in negative ENSO phases. The difference may be confusing. It\nwould be better to clarify the difference in the near and distant influences of\nENSO. The authors can display the average microseism time series for the\narea excluding SWP and SEPWA in addition to the global average and\ndiscuss the near and distance ENSO influences.\nThis is a very interesting point that we agree is worthy of further\nelaboration. We added a supplemental figure (Supp. Figure 8) that is\nconstructed identically to Figure 7b but excludes the two station\nclusters (SWP and SEPWA) and associated 11 stations that have the\nstrongest ENSO and SOI correlations. This figure demonstrates that the\nENSO-correlated global microseism level peaks remain visible even\noutside of the most strongly correlated regions (consistent with Figure\n7c). The extensive geographic signature of ENSO cycles and other time\nseries features reflects some combination of Rayleigh wave\npropagation, long-distance swell teleconnection, and meteorological\ncorrelation in trans-basin storminess.\n\n2. The 23th line from the top in Page 3: … has not been been …\nCorrected.\n\n3. The 4th line from the top in Page: The reference for the ENSO and SOI\nindices is necessary here.\nWe have now cited the seminal Smith et al. (2000) and Ropelewski and\nHalpert (1987) references with respect to these indices and in\nassociation with the associated ENSO indices supplementary figure.\n\n4. The 4th line from the top in Page8: station within the ?\nThis sentence fragment has been corrected. The associated sentence\nnow reads:\nThis correlation length scale indicates that primary microseism\namplitudes and energies provide consistent proxies for long-range\nocean wave variability.\n\n5. Figure 4: The unit of the horizontal axis should be nm/s^2.\nCorrected.\n\n6. Figure 6: It is very hard to see the correlation between microseisms time\nseries and ENSO/SOI indices. Improvements are required, such as stretching\neach panel horizontally.\nThe widths of the referenced plots in Fig. 6 have been doubled to\naddress this issue."}]

Decision letter