Satellite Evidence for Strengthened M2 Internal Tides in the Past 30 Years

Satellite altimetry sea surface height measurements from 1993 to 2022 are used to show the strengthened mode‐1 M2 internal tides in the past 30 years. Two mode‐1 M2 internal tide models M9509 and M1019 are constructed using the data in 1995–2009 and 2010–2019, respectively. The results show that the global mean M2 internal tides strengthened by 6% in energy. However, the internal tide strengthening is spatially inhomogeneous. Significantly strengthened internal tides are observed in a number of regions including the Aleutian Ridge and the Madagascar‐Mascarene region. Weakened internal tides are observed in the central Pacific. On global average, M1019 leads M9509 by about 10° (20 min in time), suggesting that the propagation speed of M2 internal tides increased. M9509 and M1019 are evaluated using independent altimetry data. The results show that M9509 and M1019 perform better for the data in 1993–1994 and 2020–2022, respectively.

satellite altimetry data into two data sets of 48 satellite years long each, from which two mode-1 M 2 internal tide models are constructed using a newly developed mapping procedure (Zhao, 2022a(Zhao, , 2022b)).The resulting models have errors as low as 1 mm, allowing us to examine the internal tide changes by contrasting the two models.We find that mode-1 M 2 internal tides strengthened and their propagation speeds increased, and that the internal tide strengthening is spatially inhomogeneous.Our finding is consistent with a recent study by Opel et al. (2023), who report that the M 2 barotropic tide weakened in the past 30 years, due to stronger barotropic-to-baroclinic tidal conversion.
The rest of this paper is arranged as follows.Section 2 describes the satellite altimetry data and the analysis methods.Section 3 presents the strengthened mode-1 M 2 internal tides by contrasting the two internal tide models.Section 4 contains conclusions and discussion.

Satellite Altimetry Data
The satellite altimetry SSH measurements used in this study are 30 years long from 1993 to 2022 (Figure 1).The measurements are made by 15 altimetry missions.The data used in this study are the along-track SSH measurements provided by the Copernicus Marine Service (https://doi.org/10.48670/moi-00146)(Taburet et al., 2019).The SSH measurements have been processed by standard corrections for atmospheric effects, surface wave biases, geophysical effects, ocean barotropic tide, and solid Earth tide (Pujol et al., 2016).In this study, we further make mesoscale correction to the along-track data using the gridded SSH product (https://doi.org/10.48670/moi-00148)following the method employed in previous studies (Ray & Byrne, 2010;Zaron & Ray, 2018;Zhao, 2022a).Two data sets in 1995-2009 and 2010-2019 are used to construct two models, which are labeled M9509 and M1019, respectively (Figure 1).The two data sets each are about 48 satellite years long, which is required to make low-noise internal tide models.Their time lengths (15 vs. 10 years) are different, because the number of concurrent satellite altimetry missions generally increases with time (Figure 1).The central times of the two data  1993-1994 and 2020-2022 are used for model evaluation.It is worth pointing out that the two data sets must be of the same size, such that the resulting models have the same noise level.Otherwise, the two models would have different noise levels, making it hard to observe the decadal changes of internal tides.

Lunar Nodal Cycle
In the 18.6-year lunar nodal cycle, the M 2 constituent is modulated by ±3.7% in amplitude and ±2.1° in phase (Haigh et al., 2011;Pugh & Woodworth, 2014;Ray, 2007).Figure 1c shows the lunar nodal modulations on amplitude and phase (out of phase by 90°).As revealed in this study, the global mean decadal change in internal tide energy is about 6%, comparable to the lunar nodal modulation; therefore, it is necessary to correct the lunar nodal modulation.We calculate their mean values in the two time windows, weighting by the number of altimetry missions.Their mean amplitudes are 0.993 and 1.021 in 1995-2009 and 2010-2019, respectively.Their mean phases are 0.5° and −0.7°, respectively.These values are used to adjust the satellite derived internal tide models M9509 and M1019 to get rid of the impact of the lunar nodal cycle.The correction can also be made by adjusting the input data with the lunar nodal factors.The two methods both work well and yield the same results.

Internal Tide Parameters
The wavelength of mode-1 M 2 internal tides is an input parameter in our mapping methods.We calculate the wavelength using hydrography in the World Ocean Atlas 2018 (WOA18; https://www.nodc.noaa.gov/OC5/woa18/) (Boyer et al., 2018).Ocean depth is based on the 1-arc-minute topography database constructed using in situ and satellite measurements (Smith & Sandwell, 1997).For a given stratification profile, the eigenvalue speeds and vertical structures are obtained by solving the Sturm-Liouville equation following (Kelly, 2016) subject to free-surface and rigid-bottom boundary conditions, where N(z) is buoyancy frequency profile, c n is the eigenvalue speed, and the subscript n is modal number.The wavelength λ n can be calculated from the eigenvalue speed c n following   = ∕ √  2 −  2 , where ω and f are the tidal and inertial frequencies, respectively.
To account for the decadal changes in ocean stratification, we compute the wavenumber using two sets of decadal hydrographic profiles in the WOA18.They are for 1995-2004 (95A4) and 2004-2017 (A5B7), respectively (Boyer et al., 2018).The wavelengths from 95A4 and A5B7 are used for mapping M9509 and M1019, respectively.Figure S1 in Supporting Information S1 shows the global maps of the wavelength of mode-1 M 2 internal tides.The change rate is calculated following (A5B7-95A4)/95A4 × 100%.One can see that the wavelengths from 95A4 and A5B7 have similar global patterns.The change rate has small-scale spatial variation, which is most likely noise due to the sparse and inhomogeneous measurements.The wavelength (phase speed) systematically increased in the Atlantic and Indian oceans, while the Pacific Ocean has a complex spatial pattern, due to its complex decadal dynamic processes.The output phase of internal tides is independent of the input wavelength.In this study, we have mapped M 2 internal tides using a set of different wavelengths from the WOA18 climatologies.
The phases of the resulting mode-1 M 2 internal tides are not affected by these slightly different wavelengths.
We calculate the depth-integrated energy E of an M 2 internal tidal wave from its SSH amplitude A following , where E n is the transfer function for mode-1 M 2 internal tides.E n is a function of ocean depth, stratification profile, latitude, tidal frequency, and baroclinic modal number.The transfer function is calculated using the WOA18 95A4 and A5B7 hydrography following Zhao et al. (2016) (Appendix A therein).Figure S2 in Supporting Information S1 shows that their global maps are very similar.The change rate of the transfer function only slightly changed over the two decades.A comparison shows that the change rates of wavelength and transfer function have similar patterns but with opposite signs (Figures S1 and S2 in Supporting Information S1).The transfer functions from 95A4 and A5B7 are used for M9509 and M1019, respectively.

Mapping Methods
Our core technique for mapping internal tides is plane wave analysis (Zhao & Alford, 2009;Zhao et al., 2016).At each grid point, internal tidal waves are determined using along-track data falling within the fitting window centered at the grid point.We fit five mode-1 M 2 internal tidal waves in the fitting window following where x and y are the east and north Cartesian coordinates, t is time, ω and k are the frequency and horizontal wavenumber of the target internal tides, respectively.To determine one wave, the amplitude and phase of a single plane wave are determined by the least-squares fit in each compass direction (with 1° increment).When the resulting amplitudes are plotted as a function of direction in polar coordinates, an internal tidal wave appears to be a lobe.The direction of the wave is thus determined from the biggest lobe.After one internal tidal wave is determined, its signal is predicted and subtracted from the original data.This procedure can be repeated to extract five internal tidal waves in arbitrary horizontal directions.
Our 3-step mapping procedure consists of two rounds of plane wave analysis and a two-dimensional bandpass filter in between (Zhao, 2022a(Zhao, , 2022b)).In the first step, mode-1 M 2 internal tides are mapped by plane wave analysis (Zhao & Alford, 2009;Zhao et al., 2016).The fitting window is 160 by 160 km, consistent with the wavelengths of mode-1 M 2 internal tides.The resulting M 2 internal tides are at a regular grid of 0.2° longitude by 0.2° latitude.In the second step, the spatially regular M 2 internal tide field is cleaned by a two-dimensional bandpass filter.The M 2 internal tide field is first converted to the wavenumber spectrum by Fourier transform in overlapping 850 by 850 km windows.The spectrum is truncated to [0.8, 1.25] times the local wavenumber.The truncated spectrum is converted back to the internal tide field by inverse Fourier transform.However, the 2D filtered internal tide field is a multiwave summed interference field.We thus call plane wave analysis the second time to decompose the cleaned internal tide field.In the third step, plane wave analysis is called again to decompose the filtered internal tide field into five internal waves at each grid point.The second-round plane wave analysis is the same as the first round, except that the input is the filtered internal tide field in the second step.Intermediate results obtained in the three steps are shown in Figures S3 and S4 in Supporting Information S1.The resulting internal tide models M9509 and M1019 are shown in Figures 2a and 2b, respectively (Zhao, 2023).

Model Errors
Model errors are estimated using background internal tides following Zhao and Qiu (2023).Background internal tides are mapped using the same altimetry data sets and following the same mapping procedure, but the tidal periods are slightly different.In this study, two periods are chosen to be 12.3706 hr (M 2 minus 3 min) and 12.4706 hr (M 2 plus 3 min), respectively.Because there are no tidal constituents at these two periods, the extracted M 2 internal tides are mainly errors.The two periods are on two sides of M 2 to avoid possible systematic bias in the internal tide field.Figures 2c and 2d show the resulting model errors.We find that large model errors are mainly associated with strong mesoscale eddies, because the models errors are mainly leaked mesoscale signals (Zhao & Qiu, 2023).Figure 2e shows that the two models both have errors as low as 1 mm, much lower than M 2 internal tides.The two models have comparable errors, because both of them are constructed using 48 satellite years of data.Therefore, the changes in M 2 internal tides can be observed from the differences between M1019 and M9509.We compare the model variance and the error variance, both of which are averaged from M9509 and M1019.We calculate the difference between the model variance  (  2

𝑀𝑀
) and twice the error variance  2f shows that the model variance is dominantly greater than twice the error variance (Δ > 0).The internal tide changes in regions of weak internal tides (Δ < 0) are masked out in this paper.

Amplitude and Energy
We first study the changes in the amplitude and energy of mode-1 M 2 internal tides.The global area-weighted mean SSH variances are calculated following Zhao (2022a).They are 9.1 and 9.8 mm 2 for M9509 and M1019, respectively.Because ocean stratification also changed during this period, we should study the internal tide changes in terms of energy, which has taken stratification into account via the transfer function.Following Zhao (2022a), the internal tide energy for M9509 and M1019 is calculated and shown in Figure S5 in Supporting Information S1.The globally integrated energies are 28.8 and 30.5 PJ, respectively, with a strengthening rate of about 6%.We thus find that the global mode-1 M 2 internal tides strengthened in terms of both SSH variance and energy over the two decades, consistent with global ocean warming.We want to stress that the strengthening rates are over 12.5 years from March 2003 to August 2015, which are the central times of the decadal windows.The real changes from 1993 to 2022 would be 2.4 times large, assuming constant change rates over the 30 years.
We next examine the spatial feature of the strengthened M 2 internal tides.Figures 3a and 3b show the changes in SSH and energy, respectively.For each case, we calculate the point-wise differences of M1019 minus M9509 and smooth the differences by 1° × 1° running mean.The results show that the changes are highly inhomogeneous, with positive changes dominating the global ocean.A number of regions with significantly strengthened M 2 internal tides are identified: (1) the Madagascar-Mascarene region, (2) the Conrad Rise and Kerguelen Plateau, (3) the western North Pacific, (4) the Tasman Sea, (5) the Kuril Strait, (6) the western South Pacific, (7) the central and eastern North Pacific, (7a) the Aleutian Ridge, (8) the Amazon continental region, (9) the Northeast Atlantic Ocean, and (10) the Walvis Ridge in the southern Atlantic Ocean.All these regions are associated with strong M 2 internal tides.The highest strengthening rate is about 30% to the south of the Aleutian Ridge.The second highest rate is about 25% in the Madagascar-Mascarene region.The third highest rate is about 20% in the western North Pacific.For both the Tasman Sea and the central and eastern North Pacific, the strengthening rate is about 15%.
There are some scattered small regions where mode-1 M 2 internal tides weakened during this period.Most of these regions are associated with weak internal tides.Among them, one remarkable region is the central Pacific Ocean (region 11), where mode-1 M 2 internal tides weakened by −10% in terms of energy.It is worth pointing out that the M 2 internal tides in the tropical Pacific Ocean are subject to significant inter-annual variation, as evidenced by the different M 2 internal tide fields in 2017-2020 (Zhao, 2022a).In particular, the M 2 internal tides in this region are affected by El Niño-Southern Oscillation (Zhao, 2016).

Phase and Speed
The phase changes from M9509 to M1019 are calculated point by point and smoothed in 1° × 1° windows.Figure 3c shows the resulting phase changes.Positive phase differences mean that M1019 leads M9509 in time, and vice versa.In other words, positive values mean that the speeds of M 2 internal tides increased over the two periods.Figure 3c shows that the phase differences are dominantly positive.The global mean phase difference is about 10° (20 min in time).We find that positive phase differences are throughout the Atlantic and Indian Oceans, suggesting that the propagation speeds of M 2 internal tides became faster.It is consistent with (a) the faster speeds in 2005-2017 observed using the WOA18 decadal hydrography, and (b) the strengthened M 2 internal tides in the Atlantic and Indian Oceans (Figures 3a and 3b).The Pacific Ocean, however, has a complex spatial pattern.The largest phase difference occurred in the southern Pacific Ocean.In region B (50°-30°S, 170°-140°W), the mean phase difference is about 40° (80 min in time).Figure 3c also shows that negative phase differences are observed in the tropical Pacific Ocean.In region A (the western North Pacific), the mean phase difference is −5° (−10 min).Note that, in the Pacific Ocean, the spatial feature in energy (Figure 3b) and the spatial feature in phase (Figure 3c) are different.It suggests that the internal tides in the tropical Pacific Ocean are affected by a complex ocean dynamic environment.

Model Evaluation
M9509 and M1019 are evaluated using independent altimetry data in 1993-1994 and 2020-2022 (Figure 1).For each SSH measurement of known time and location, the internal tide signal is predicted using the model under evaluation, and subtracted from the raw satellite data.The variance reduction is the difference before and after the internal tide correction.In the end, the variance reductions from all SSH measurements are binned into half-overlapping 1° by 1° boxes.The results show that both models cause significant variance reductions in the global ocean (Figures S6-S11 in Supporting Information S1).We further assess their performance by comparing the variation reductions they explained.In Figure 4, red regions indicate that M1019 reduces more variance, and blue regions indicate that M9509 reduces more variance.We find that M9509 overall performs better in 1993 and 1994 (Figures 4a and 4b), while M1019 performs better in 2020, 2021, and 2022 (Figures 4c-4e) (Figures 4c-4e).Note that the complex spatial feature in Figures 4a and 4b are consistent with the WOA18 decadal changes in speed.This contrast implies that M 2 internal tides systematically changed over the past 30 years.Thus, the performance of internal tide models is affected by the long-term changes of internal tides.We find that M1019 performs better for satellite altimetry data in recent years, because M1019 is constructed using recent satellite altimetry data in 2010-2019.We suggest that M1019 is a better model for the Surface Water

Figure 1 .
Figure 1.Satellite altimetry data.(a) Time coverage showing three decades of sea surface height (SSH) measurements from 1993 to 2022 made by 15 altimetry missions.(b) The yearly number of concurrent altimetry missions.Mode-1 M 2 internal tide models M9509 and M1019 are constructed using the data in 1995-2009 and 2010-2019, respectively.The central times of the two data sets are March 2003 and August 2015, respectively (triangles).The data in 1993-1994 and 2020-2022 are reserved for model evaluation.(c) The 18.6-year lunar nodal cycle of the M 2 tidal constituent.

Figure 2 .
Figure 2. Internal tide models and errors.(a) M9509.(b) M1019.(c) M9509 error.(d) M1019 error.(e) Zonal mean sea surface height (SSH) amplitudes showing that model errors are about 1 mm, much lower than M 2 internal tides.(f) Model variance minus twice error variance.Blue regions indicate weak internal tides where the internal tide changes are not studied.

Figure 3 .
Figure 3. Strengthened mode-1 M 2 internal tides.Shown are the differences of M1019 minus M9509.(a) Sea surface height amplitude.(b) Depth-integrated energy (see Figure S5 in Supporting Information S1).(c) Phase.Significantly strengthened internal tides are observed in the following regions: (1) the Madagascar-Mascarene region, (2) the Conrad Rise and Kerguelen Plateau, (3) the western North Pacific, (4) the Tasman Sea, (5) the Kuril Strait, (6) the western South Pacific, (7) the central and eastern North Pacific, (7a) the Aleutian Ridge, (8) the Amazon continental region, (9) the Northeast Atlantic Ocean, and (10) the Walvis Ridge.However, weakened internal tides are observed in central Pacific Ocean (11).In (c), positive phase differences are observed throughout the Atlantic and Indian Oceans; however, negative phase differences are observed in regions such as the tropical North Pacific (A).The largest phase differences occur in the southern Pacific Ocean (B).

Figure 4 .
Figure 4. Model performance.M9509 and M1019 are employed to make internal tide correction to six data sets from different years.Shown are the variance reduction differences explained by M1019 and M9509.M9509 performs better in 1993-1994, while M1019 performs better in 2020-2022 and January-June 2023.See Figures S6-S11 in Supporting Information S1 for variance reductions for the six data sets.