the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Validating the spatial variability in the semidiurnal internal tide in a realistic global ocean simulation with Argo and mooring data
Gaspard Geoffroy
Jonas Nycander
Maarten C. Buijsman
Jay F. Shriver
Brian K. Arbic
The autocovariance of the semidiurnal internal tide (IT) is examined in a 32 d segment of a global run of the HYbrid Coordinate Ocean Model (HYCOM). This numerical simulation, with 41 vertical layers and $\mathrm{1}/\mathrm{25}$^{∘} horizontal resolution, includes tidal and atmospheric forcing, allowing for the generation and propagation of ITs to take place within a realistic eddying general circulation. The HYCOM data are in turn compared with global observations of the IT around 1000 dbar, from Argo float parkphase data and mooring records. HYCOM is found to be globally biased low in terms of the IT variance and decay of the IT autocovariance over timescales shorter than 32 d. Except in the Southern Ocean, where limitations in the model cause the discrepancy with in situ measurements to grow poleward, the spatial correlation between the Argo and HYCOM tidal variance suggests that the generation of lowmode semidiurnal ITs is globally well captured by the model.
 Article
(10559 KB)  Fulltext XML
 BibTeX
 EndNote
Internal tides (ITs) are internal waves generated by the interaction of tidal currents with rough bathymetry. The radiated wave beams can travel thousands of kilometers (e.g., Zhao et al., 2016; Buijsman et al., 2020), over which they undergo dissipative processes as they interact with the eddying ocean and bottom topography to eventually break (MacKinnon et al., 2017). The dissipation of ITs represents a major source of vertical mixing in the ocean interior (de Lavergne et al., 2020). As such, ITs have a key influence on the ocean state (Melet et al., 2016) and therefore on the global climate system (see Melet et al., 2022, for a review on the subject).
At any given position in a stationary medium, tidally forced waves would have a constant phase difference to the astronomical forcing. However, since they propagate within the timevarying ocean circulation, ITs are subject to a variety of mechanisms that cause their phase difference to the tidal forcing at their generation site to shift with time (Rainville and Pinkel, 2006; Shriver et al., 2014; Zaron and Egbert, 2014; Buijsman et al., 2017). In other words, ITs lose coherence by interacting with the eddying ocean. They decorrelate: the autocovariance of a time series representing the internal tide variability at a fixed position (away from the source) inevitably decays with time lag (CasparCohen et al., 2022; Geoffroy and Nycander, 2022).
Only the coherent fraction of the IT energy decays with time lag, that is the energy carried by waves with a constant phase difference to the astronomical forcing. Conversely, the incoherent fraction of the energy grows so that the total IT energy, or variance, at a given location (the sum of the coherent and incoherent components) is unaffected by the decorrelation. For very long time lags, the coherent and the complementary incoherent asymptotic limits are often called stationary and nonstationary variance, respectively. While the global field of the stationary (lowmode) IT is widely considered to be well constrained by multiyear satellite altimetry observations (Dushaw et al., 2011; Ray and Zaron, 2016; Zaron, 2019; Zhao, 2019), the nonstationary component is not.
Some current highresolution global ocean circulation models enable the estimation of the barotropic and internal tides concurrently with the ocean circulation (Arbic et al., 2010; Shriver et al., 2012; Buijsman et al., 2020). Such fully nonlinear numerical simulations also incorporate the interactions of the generated ITs with the eddying ocean and its boundaries. The model utilized in this study is the HYbrid Coordinate Ocean Model (HYCOM; Chassignet et al., 2006), with 41 vertical layers and $\mathrm{1}/\mathrm{25}$^{∘} horizontal resolution. The literature on the validation of HYCOM with observations is already vast (see Buijsman et al., 2020, and Arbic, 2022, for recent accounts). In the latest developments, surface drifters have been used to validate the geographical variability in the kinetic energy in various frequency bands at the global scale (Arbic et al., 2022). While the global coverage of drifter data is comparable to that of satellite altimetry, the contribution from the baroclinic tides to the kinetic energy observed by drifters has not been determined yet. Hence, until recently, only altimetry could unveil the geographical variability in the ITs at the global scale.
The empirical mapping of ITs from altimetry remains challenging for various reasons (Egbert and Ray, 2017). Most notably, the long sampling intervals of altimeters and the low signaltonoise ratio preclude any direct estimation in the time domain of the total IT variance at a single location. Notwithstanding, Zaron (2017) analyzed alongtrack wavenumber spectra of the sea surface height to map the total and nonstationary semidiurnal IT variance (for the baroclinic mode 1 only). The author found a globalmean ratio of nonstationary to total semidiurnal IT variance of 44 %. He also outlined the spatial inhomogeneity of the tidal variability, with this ratio being larger than 50 % in much of the equatorial Pacific and Indian oceans (Zaron, 2017). In a comparison with data from HYCOM, Nelson et al. (2019) showed the “kspace” methodology of Zaron (2017) to miss a large fraction of both the nonstationary and total variance. Nevertheless, the spatial correlation of the nonstationary fraction (i.e., the ratio of the nonstationary to total variance) between the model and altimetry suggests that the model at least qualitatively captures the generation of ITs and their interactions with the background circulation (Nelson et al., 2019).
Recently, Geoffroy and Nycander (2022) used observations from Argo parkphase data (Argo, 2000) to empirically map the variance of the semidiurnal IT at 1000 dbar. The high sampling rate of the floats captures the total variance of the IT, i.e., the autocovariance at short time lags, before any significant loss of coherence occurs and after most of the noise has decorrelated. On the other hand, the Lagrangian sampling of the drifting floats results in decorrelating effects on top of the decorrelation of the IT itself (Zaron and Elipot, 2021; CasparCohen et al., 2022; Geoffroy and Nycander, 2022). Following CasparCohen et al. (2022), we call the effects of the Lagrangian sampling on the autocovariance apparent decorrelation. This apparent decorrelation cannot be disentangled from the decorrelation of the IT (i.e., the decorrelation due to interactions with the background circulation) using Argo data only. Thus, to gain insights into the IT decorrelation around 1000 dbar, one can instead apply the methodology of Geoffroy and Nycander (2022) to Eulerian observations from moorings. In the present work we compare observations of the (total) variance and decorrelation of the semidiurnal IT around 1000 dbar from Argo floats and moorings, respectively, to data from a global HYCOM run. Contrarily to other recent validations of HYCOM with mooring data (e.g., Ansong et al., 2017; Luecke et al., 2020), the Eulerian component of our analysis is not meant as a standalone pointtopoint comparison. Rather, it is designed to bolster and extend the main analysis of the Lagrangian data.
The remainder of this work is organized as follows: in Sect. 2, the in situ datasets and the numerical simulation are presented. In Sect. 3, we review the methodology of Geoffroy and Nycander (2022) in an example location from the Lagrangian and Eulerian perspectives. In Sect. 4, the main results are presented. We compare the geographical variability in the variance and decorrelation of the semidiurnal IT obtained from the in situ data and numerical simulation. Then, the model data are used to quantify the decorrelation induced by the Lagrangian sampling. Finally, we outline the potential biases affecting the datasets. We conclude in Sect. 5 by discussing the results and giving a summary.
The different comparisons in this work are all done in terms of vertical displacement of the isotherms. The measurable variables needed to compute vertical isotherm displacement at a fixed depth are the temperature anomaly and vertical temperature gradient at that depth. In this section we briefly describe the temperature time series from each dataset.
2.1 Argo floats with Iridium communication
We use data from a global collection of Argo Iridium floats deployed by the University of Washington as part of the National Ocean Partnership Program during the period 2004–2022. Between the descending and ascending profiling phases, these floats also record temperature and pressure with an hourly resolution while adrift at 1000 dbar. This socalled park phase typically lasts 10 d. As in Hennon et al. (2014) and Geoffroy and Nycander (2022), we use the pressure records to correct for the small departures of a float from its drifting control level during a park phase.
Stitching together data from successive cycles, by filling the time between park phases (typically 6 h) with NaN (not a number) values, one can construct longer time series. In this study, we use segments of 32 d of data (i.e., the duration of the segment of numerical simulation we are using). The sampling period of the park phase can occasionally vary by more than a few seconds. To ensure evenly spaced time series, we linearly interpolate each concatenated record of 32 d onto a time axis with a constant 1 h step. Any interpolated value lying between two original records that are more than 1.5 h apart is replaced by NaN.
The position of the floats can only be determined when they reach the surface. We assume straight trajectories between two successive surfacings (typically 10 d apart). This assumption has been shown to be reasonable by Geoffroy and Nycander (2022), especially since we discard segments of data for which the mean speed of the float is larger than 0.1 m s^{−1}. This criterion is primarily intended to avoid contamination by lee waves: a flow with speed U≈0.1 m s^{−1} passing over a bathymetric feature with horizontal length scale λ≈10 km will generate lee waves with angular frequency $\mathit{\omega}=(\mathrm{2}\mathit{\pi}/\mathit{\lambda})U$ rad s^{−1} of the order of the diurnal frequency.
The Argo dataset used in this work is an updated version of the one used by Geoffroy and Nycander (2022). After the base processing and quality controls, we are left with 22 414 valid 32 d segments from 891 individual floats. A more detailed description of the Argo data processing and quality controls we employed is available from Hennon et al. (2014) and Geoffroy and Nycander (2022).
2.2 Global MultiArchive Current Meter Database (GMACMD)
The Global MultiArchive Current Meter Database (GMACMD; Scott et al., 2011) compiles tens of thousands of oceanographic time series from moorings. Previous model–data validation efforts involving both HYCOM and mooring data notably include Luecke et al. (2020), where the authors compared the temperature variance and kinetic energy over various frequency bands in realistic global ocean simulations to more than 3800 instrumental records.
Here, we extracted 331 temperature time series spanning 1972–2010 and meeting the following criteria:

The mooring lies in water deeper than 2000 m.

The record is longer than 64 d, with a sampling interval shorter than 3 h, for adequate resolution of the semidiurnal tidal signals.

The instrument depth is within ±200 from 1000 m.
One particular mooring is used as an example in Sect. 3. It recorded during 366 d in the years 1982–1983 at the position 28.00^{∘} N, 151.95^{∘} W (north of the Hawaiian Ridge). As previously documented by Alford and Zhao (2007), we refer to the instrument at 1119 m depth as mooring no. 2. This instrument sampled temperature with a 0.25 h resolution.
2.3 HYbrid Coordinate Ocean Model (HYCOM)
This study uses 32 d, from 20 May to 20 June 2019, of hourly output at 1000 m depth from a global run of HYCOM, with 41 vertical layers and $\mathrm{1}/\mathrm{25}$^{∘} horizontal resolution. The nondataassimilative simulation, designated “GLBy190.04”, includes realistic tidal and atmospheric forcing enabling the generation and propagation of ITs within the eddying general circulation. The highresolution 2D fields are complemented by monthly mean 3D fields of temperature and salinity subsampled to 1^{∘}.
A Lagrangian analysis of the simulation is used for a direct model–data validation with Argo floats. The Argo quasiLagrangian sampling is mimicked by releasing 41 644 particles randomly across the world oceans. We let the particles be advected by the 2D velocity field at 1000 m for 32 d while sampling temperature with an hourly resolution. This Lagrangian sampling of HYCOM is achieved using the software Parcels (Van Sebille et al., 2021). The Lagrangian simulation uses a classic Runge–Kutta method for computing the advection of the particles (with 5 min integration time step). As for the Argo data, we discard particles with a mean speed larger than 0.1 m s^{−1}. We also discard any particle crossing the 1000 m isobath during the simulation.
In this section, we present a local example to introduce the methods that are used in Sect. 4 for the global comparison. We start by comparing the Argo observations to the Lagrangian sampling of the HYCOM data. Next, we investigate the effects of the drift by comparing Lagrangian and Eulerian samplings of the numerical simulation. We end the section by introducing the comparison between the Eulerian HYCOM time series and mooring observations.
3.1 Lagrangian sampling of the isopycnal displacement at 1000 dbar
As in Hennon et al. (2014) and Geoffroy and Nycander (2022), we define the vertical isotherm displacement observed by a Lagrangian particle ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{L}}$ as
where $\stackrel{\mathrm{\u203e}}{T}$ is the time average of the temperature measurements T(t) over a particle trajectory, and $(\mathrm{d}T/\mathrm{d}P{)}_{\mathrm{1000}}(t)$ is the temperature gradient at 1000 dbar. Hence, the vertical isotherm displacement is simply computed as the temperature anomaly recorded by a drifting particle at 1000 dbar divided by the vertical temperature gradient at that depth. For HYCOM particles, we compute $(\mathrm{d}T/\mathrm{d}P{)}_{\mathrm{1000}}$ as the mean, between 900 and 1100 m, of the vertical gradient calculated from the modeled monthly mean 3D temperature field. We then linearly interpolate the temperature gradient in the horizontal to the successive particle positions (hence the time dependence). To avoid any spurious displacements, we set ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{L}}$ to NaN whenever the magnitude of the temperature gradient is smaller than $\mathrm{3}\times {\mathrm{10}}^{\mathrm{5}}$ ^{∘}C dbar^{−1}.
In the case of Argo floats, Eq. (1) is evaluated over each park phase. $\stackrel{\mathrm{\u203e}}{T}$ then represents the time average of the temperature measurements T(t) over a park phase, and the vertical temperature gradient is estimated from the average of the two neighboring ascending profiles. Specifically, $(\mathrm{d}T/\mathrm{d}P{)}_{\mathrm{1000}}(t)$ is computed as the mean, within 100 dbar of the parking pressure, of the vertical gradient calculated from the average temperature profile. As for the HYCOM particles, we discarded the data from the whole park phase whenever the magnitude of the temperature gradient is smaller than $\mathrm{3}\times {\mathrm{10}}^{\mathrm{5}}$ ^{∘}C dbar^{−1}. As explained in Sect. 2.1, ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{L}}$ time series from successive park phases are stitched together to constitute 32 d time series.
For both Argo floats and HYCOM Lagrangian particles, the lowfrequency background activity is filtered out from ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{L}}$ using a fourthorder Butterworth filter with a cutoff frequency of 0.3 cpd (cycles per day). The inertial peak is not removed by this filter poleward of about ±10^{∘} latitude.
3.2 Averaging sample autocovariance series
The HYCOMderived Lagrangian time series ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{L}}$ can be analyzed in the same way as in Geoffroy and Nycander (2022). We illustrate the mapping process below. The geographical patch presented here is also described in Geoffroy and Nycander (2022). For this example area, Fig. 1 shows eight 32 d segments of Argo data and 13 HYCOM particles selected using their median position. The circular patch of 200 km radius containing the latter median positions is centered on mooring no. 2.
From a finite time series η(t) one can calculate the sample autocovariance $\widehat{R}\left(\mathit{\tau}\right)$:
where N is the total number of observations, and $\widehat{\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}$ is the sample mean of the series. Note that $\widehat{R}\left(\mathrm{0}\right)$ is the sample variance of the series. Here, N is taken as the number of nonmissing observations, accounting for gaps in the Argo time series (during the descent, main profiling, and surface transmission phases of the float cycle or because of failed quality controls). The sample autocovariance as defined in Eq. (2) is an unbiased estimator of the “true” autocovariance (i.e., $\widehat{R}\left(\mathit{\tau}\right)$ converges to the true value R(τ) for infinitely large N). In the remainder of this article we only compute the variance and autocovariance of vertical isotherm displacement time series. For the sake of brevity, we simply refer to them as the variance and autocovariance.
The sample autocovariances for all HYCOM particles within the circular patch shown in Fig. 1 are averaged to obtain a localmean autocovariance ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$. The localmean autocovariance from the Argo data, ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{Argo}}$, is computed in the same way. The standard error of the localmean autocovariance is computed as
where STD is the standard deviation of the sample autocovariances over the subset, and N_{p} is the number of particles in the subset. Figure 2a demonstrates the good agreement of the two datasets until τ≈200 h. Past this time lag limit, ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{Argo}}$ falls under its 95 % confidence interval (95 % C.I. = ±2 SEM) and thus cannot be considered significantly different from zero.
A handy tool for monitoring the evolution of the autocovariance of the semidiurnal IT is complex demodulation. Here, it consists of the least squares fitting of Acos (ω_{SD}τ)+Bsin (ω_{SD}τ), where ${\mathit{\omega}}_{\mathrm{SD}}=({\mathit{\omega}}_{{\mathrm{M}}_{\mathrm{2}}}+{\mathit{\omega}}_{{\mathrm{S}}_{\mathrm{2}}})/\mathrm{2}$ is the semidiurnal frequency, to the sample autocovariance in 48 h windows with 75 % overlap. This is equivalent to fitting Ccos (ω_{SD}τ+Φ), where C and Φ are the amplitude and phase, respectively. We then define the complex demodulate as $C=\sqrt{{A}^{\mathrm{2}}+{B}^{\mathrm{2}}}$. This positive definite amplitude follows the envelope of the modulated sinusoidal with frequency ω_{SD}. Note that this definition differs from the complex demodulation used in Geoffroy and Nycander (2022), where the authors fitted Ccos (ω_{SD}τ), i.e., with Φ=0. In practice, we compute the complex demodulate of the autocovariance following the harmonic analysis method of Cherniawsky et al. (2001). The demodulation over a given 48 h window is performed in two iterations. We first fit Acos (ω_{SD}τ)+Bsin (ω_{SD}τ) to the signal and compute the root mean square error (RMSE). We then repeat the fitting with a trimmed signal that excludes values outside of ±2 RMSE from the previously fitted curve.
Complex demodulation is just a convenient way of finding the envelope of the sample estimate of an underlying true oscillating function at a given frequency. However, as an estimate of the envelope of the true oscillating function, the complex demodulate can be shown to be biased high (see Appendix A). This bias is greatly mitigated (i) at short time lags and (ii) when the sample size is large (e.g., when demodulating regional or globalmean autocovariances). Throughout this paper we only rely on complex demodulation in one case or the other. A conservative estimate of the uncertainty in the envelope of the sample function (i.e., the uncertainty in the demodulate) is the uncertainty in the sample function itself. If this uncertainty range is larger than the envelope of the sample function, the conclusion is not that the envelope of the true function can be negative but that it is not significantly different from zero. Here, we evaluate the uncertainty affecting the complex demodulate of a mean autocovariance series over a 48 h window by computing the median of the standard error over that window, hereinafter denoted $\stackrel{\mathrm{\u0303}}{\mathrm{SEM}}$. The corresponding 95 % confidence intervals are taken as $\pm \mathrm{2}\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u0303}}{\mathrm{SEM}}$.
Following Geoffroy and Nycander (2022), the semidiurnal IT variance can be estimated from the first 48 h demodulate at the semidiurnal frequency. In Appendix B we investigate the potential contamination of the first demodulate by the nontidal variability (noise). We found that a background noise can contribute either positively or negatively (a consequence of filtering the time series) to the amplitude of the first demodulate, depending on its characteristic timescale. Notwithstanding, for a typical signaltonoise ratio and timescale of the background noise observed by Argo floats, the first demodulate was seen to remain a conservative estimate of the IT variance. The effects of the nontidal variability on the IT variance estimate as computed in Geoffroy and Nycander (2022) thus do not put their results into question.
Since HYCOM does not resolve the full spectrum of the oceanic variability, in particular the variability associated with short timescales, the corresponding stochastic noise is expected to differ from the one captured by the in situ data (in terms of both variance and characteristic timescale). The contamination of the first demodulate by this noise would therefore account for a systematic bias when comparing the simulated data with observations. To limit this, we chose to consistently subtract an estimate of the autocovariance of the nontidal variability from the sample mean autocovariance before computing the complex demodulates. The way we obtain such an estimate is described in Appendix C. Namely, we fit a model of the autocovariance of a tidal variability on top of a background noise to the sample autocovariance by nonlinear least squares. Besides the noise parameters, this also provides us with an estimate of the tidal variance. Still, in this work, we chose to perform the comparisons in terms of the first demodulate for it is a conservative and more robust estimate of the IT variance.
The result of the complex demodulation at the semidiurnal frequency applied to our two (noisecorrected) mean autocovariance series is presented in Fig. 2b. The first 48 h complex demodulate is taken as a direct estimate of the semidiurnal IT variance. From the Argo data plotted in Fig. 2b we obtain 24.9 m^{2}, with a 95 % C.I. of [16.7, 33.2] m^{2}. For this local example, the first demodulate of ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$ is almost identical to the first demodulate of ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{Argo}}$. Apart from the first couple of demodulates, the HYCOM demodulate series appears consistently smaller than the Argo one.
3.3 Eulerian perspective
The decay with time lag of the demodulates represented as red crosses in Fig. 2b mirrors the decorrelation captured by the Lagrangian sampling of the Argo floats. The motion of the instruments results in decorrelating effects that cannot be isolated from the loss of coherence of the IT. However, by analyzing the HYCOM data within a Eulerian framework, one can directly monitor the decorrelation of the IT. The Eulerian HYCOM time series can then be compared with observations from moorings.
In the Eulerian framework, our methods remain practically unchanged. We now define the vertical isotherm displacement at a given location as
where $\stackrel{\mathrm{\u203e}}{T}(x,y)$ is the time average of the temperature field $T(x,y,t)$, and $(\mathrm{d}T/\mathrm{d}P{)}_{\mathrm{1000}}(x,y)$ is the temperature gradient at 1000 m, calculated from the modeled monthly mean 3D temperature field. As for the Lagrangian time series, we apply a fourthorder Butterworth highpass filter with a cutoff frequency of 0.3 cpd on ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{E}}$.
For each HYCOM particle, we compute 32 d long time series of ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{E}}$ at the successive positions of the particle subsampled at a 12 h rate. We then calculate the sample autocovariance from each of these time series and average the results over the particle's trajectory. As in the Lagrangian framework, the resulting averages can be averaged over different particles to compute local and globalmean autocovariance series. By estimating the Eulerian autocovariance along the Lagrangian trajectories we account for the geographic variability in the IT. Also, our Eulerian sample autocovariance estimates contain many more degrees of freedom than their Lagrangian equivalents. As a result, the uncertainty affecting the Eulerian autocovariance estimates is much smaller.
Equation (4) is also used to compute vertical displacement time series from the mooring temperature records. Here, the vertical temperature gradient is computed from the annualmean climatology (WOA18; Boyer et al., 2018) as an average of the temperature gradient within 100 m of the instruments' depth. As for the HYCOM and Argo data, we discarded instruments for which the magnitude of the temperature gradient is smaller than $\mathrm{3}\times {\mathrm{10}}^{\mathrm{5}}$ ^{∘}C dbar^{−1}. Each mooring time series is split into successive 32 d segments. We then compute the sample autocovariance for each highpassfiltered segment and average them to obtain a mean autocovariance.
The Eulerian sampling serves two main purposes: (i) validating the variance of the IT measured by the Lagrangian particles and (ii) comparing the decorrelation of the IT in the HYCOM data with mooring observations. We illustrate these two aspects for the local example introduced in Sect. 3.2 in Figs. 3 and 4, respectively.
Figure 3 shows the localmean autocovariance series at 1000 dbar computed from both the Lagrangian ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{L}}$ (solid black curve) and Eulerian ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{E}}$ (solid red curve). As expected, the two autocovariance series demonstrate a close agreement at short time lags, before the motion of the particles causes the Lagrangian ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$ to decay faster: the first demodulate of the Eulerian ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$ (red crosses in Fig. 3b), reading 28.4 m^{2}, with a 95 % C.I. of [27.6, 29.1] m^{2}, is similar to the first demodulate of the Lagrangian ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$ (black crosses in Figs. 3b and in Fig. 2b). In contrast with ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$, ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$ is not affected by the apparent decorrelation, and it remains close to the mean autocovariance computed from the (Eulerian) mooring no. 2 time series ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{Moor}}$ at all time lags (see Fig. 4). In conclusion, for this local example, the model agrees very well with the in situ observations, in terms of both variance and decorrelation of the semidiurnal IT at 1000 dbar.
4.1 The semidiurnal IT variance at 1000 dbar
We bin the global collection of HYCOM particles based on their median position using circular geographical patches of 200 km radius centered on a regular 2.5^{∘} × 2.5^{∘} grid. Here, we use only particles for which (i) the mean speed is lower than 0.1 m s^{−1} (to avoid contamination by lee waves) and (ii) the variance of ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{L}}$ is lower than 1×10^{4} m^{2}. The second criterion accounts for 41 discarded particles (out of 41 644) that we consider to be outliers, most of which are located in the Labrador Basin or in shallow water close to Antarctica. Including it instead would only marginally affect our results. We did not further investigate the cause of these extreme values.
As in Sect. 3.2 and 3.3, we compute an autocovariance series for each HYCOM particle, in the Lagrangian and Eulerian framework, separately. Then, we average these autocovariances over the corresponding patches to obtain localmean autocovariance series ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$ and ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$. For each geographical patch, we get local estimates of the semidiurnal IT variance from the first 48 h complex demodulate of ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$ and ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$.
We start by checking how the Lagrangian sampling affects the semidiurnal IT variance. Figure 5 shows the first 48 h complex demodulate of ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$ plotted as a function of the first demodulate of ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$ for our collection of geographical patches. The agreement is close to perfect (with a Pearson r squared or coefficient of determination r^{2}≈0.99 and 0.76 in the log–log and linear domain, respectively). Thus, the motion of the Lagrangian particles, and therefore of the Argo floats, has no significant impact on the measured variance of the IT.
We can then map the semidiurnal IT variance (here taken as the first 48 h complex demodulate of ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$) and the associated $\stackrel{\mathrm{\u0303}}{\mathrm{SEM}}$ on a 2.5^{∘} × 2.5^{∘} grid (see Fig. 6a and b, respectively). In each figure we show only the bins which yield an IT variance larger than one $\stackrel{\mathrm{\u0303}}{\mathrm{SEM}}$. Note that the latter criterion is less binding than the one used for the global maps in Geoffroy and Nycander (2022) since the uncertainties are smaller here, by definition. Figure 6a can be directly compared with the global map of the semidiurnal IT variance computed from Argo data (see Fig. 6c and d, updated from Geoffroy and Nycander, 2022). As documented in Buijsman et al. (2020), HYCOM is known to be subject to a thermobaric instability (TBI) in the North Pacific (dashed black rectangle in Fig. 6a). Since only a few patches of Argo data are available at the edge of this TBI area, we do not exclude it from the subsequent analysis.
The main patterns visible in Fig. 6a and c broadly agree, with energy radiating away from lowmode IT generation hotspots, namely near Madagascar, Hawaii, French Polynesia, and the western South Pacific. In Fig. 7a we show the Argoderived semidiurnal IT variance plotted as a function of the corresponding one from HYCOM. The two datasets are well correlated (r^{2}=0.52 and 0.38 in the log–log and linear domain, respectively), but the semidiurnal IT variance is systematically smaller in HYCOM than in the Argo data.
We investigate this bias by looking at the geographical distribution of the HYCOMtoArgo semidiurnal IT variance ratio. For presentation purposes, instead of the latter ratio we plot the proxy ${\mathit{\sigma}}_{\mathrm{HYCOM}}^{\mathrm{2}}/({\mathit{\sigma}}_{\mathrm{HYCOM}}^{\mathrm{2}}+{\mathit{\sigma}}_{\mathrm{Argo}}^{\mathrm{2}})$ in Fig. 7b. This statistic is robust to outliers and maps a pair of variances to the closed range [0, 1]. The ratio is relatively homogeneous globally (in the range [0.25, 1.5], corresponding to [0.2, 0.6] for our proxy), except over the Southern Ocean, where the Argoinferred IT variance is significantly larger. Furthermore, we show the zonal mean of the Argo and HYCOMderived semidiurnal IT variance as a function of latitude in Fig. 7c. The zonalmean variance from Argo is generally larger, except around 40^{∘} N, where the HYCOMinferred zonalmean variance peaks at 1.6 times the Argo one. Discarding the data in the TBI area has virtually no effect on this peak (not shown). More noteworthy than this localized feature, the ratio of the HYCOM to Argo zonalmean variances decreases approaching the poles (poleward of ±50^{∘} latitude, not shown). In contrast, equatorward of ±50^{∘}, this ratio remains fairly constant, oscillating around a mean value of 0.73 with a standard deviation of 0.26. The global mean and standard deviation of the ratio are 0.61 and 0.33, respectively.
The representativity of the zonalmean variances is smaller north of 40^{∘} N because of the scarcity of available data (solid red curve in Fig. 7c). We therefore focus on the pronounced discrepancy affecting the Southern Ocean. We gather the data available over the unmasked bins in Fig. 7b into two groups, north and south of 50^{∘} S. For each group and the global collection of bins, we compute a mean autocovariance by averaging the corresponding Argo and HYCOM localmean autocovariance series (see Fig. 8). In Table 1, we summarize the semidiurnal IT variance values computed as the first 48 h demodulate of the mean autocovariance for each group. A significant fraction of the divergence visible in the globalmean autocovariance at short time lags (see Fig. 8a and b) can be explained by the larger discrepancy south of 50^{∘} S (see Fig. 8e and f). In this region, the Argoderived semidiurnal IT variance is close to 3.5 times larger than the simulated one. The agreement is better north of 50^{∘} S, where the corresponding factor is 1.5 (see Fig. 8c and d). In a separate section we discuss possible explanations for both the lower semidiurnal IT variance in HYCOM globally and the even lower simulated IT variance in the Southern Ocean.
4.2 The decorrelation of the semidiurnal IT at 1000 dbar
In contrast to Argo floats, the Eulerian sampling of moorings allows us to directly monitor the decorrelation of the IT. In a procedure similar to that in Sect. 4.1, we bin the global collection of HYCOM particles based on their median position, but this time using geographical patches (of 200 km radius) centered on the available moorings. We compute a sample autocovariance in the Eulerian framework for each particle and average them within the corresponding patches to obtain localmean autocovariance series ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$. The mooring time series in a given patch is split into successive 32 d segments. We then compute a sample autocovariance for each segment and average them to obtain a mean autocovariance ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{Moor}}$. Again, we use only particles for which (i) the mean speed is lower than 0.1 m s^{−1} (to avoid contamination by lee waves) and (ii) the variance of ${\mathit{\eta}}_{\mathrm{1000}}^{\mathrm{L}}$ is lower than 1×10^{4} m^{2} (outliers). After some additional quality controls, mainly discarding bins where either the mooring or HYCOM complex demodulates fall under one $\stackrel{\mathrm{\u0303}}{\mathrm{SEM}}$ within 15 d of time lag, we are left with 172 moored instruments and the corresponding patches of simulated data.
For all the datasets used in this study, we found the probability distribution of the global collection of localmean autocovariance at any given time lag to be skewed (not shown). This is not an issue when computing average statistics from the Argo or the corresponding simulated data, since the number of samples (i.e., the number of geographical bins) is very large, and the sample mean is therefore expected to be normally distributed, by virtue of the central limit theorem. In contrast, geographical bins where mooring data are available are fewer. Thus, the influence of the tail of the distribution on the sample mean is larger when analyzing the relatively small collection of bins where mooring data are available than when considering the global collections of Argo and HYCOM data. This precludes the use of statistics that assume a normal distribution when describing regional or even global averages of the autocovariances computed from moorings. To limit the effects of skewness, we discard bins for which either the mooring or the simulated first 48 h demodulate is above the 95th percentile of its observed distribution (here ${P}_{\mathrm{Moor}}^{\mathrm{95}}\approx \mathrm{575}$ m^{2} and ${P}_{\mathrm{HYCOM}}^{\mathrm{95}}\approx \mathrm{330}$ m^{2}, respectively). The latter criterion accounts for 12 additional discarded bins, all located in moderate to strongmesoscaleactivity areas. Even after discarding these extreme samples, the data remain highly skewed.
The moorings may not offer as much spatial coverage as the Argo floats do, but they still provide an opportunity to validate the geographical variability in the semidiurnal IT variance in HYCOM. As in Sect. 4.1, we can get local estimates of the semidiurnal IT variance from ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$ and ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{Moor}}$ in each geographical patch centered on a mooring. Figure 9a shows the scatterplot of the first 48 h complex demodulate of ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$ as a function of the first demodulate of ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{Moor}}$ for our collection of geographical bins. As with the Lagrangian data (see Fig. 7a), the correlation is good (r^{2}=0.51 and 0.40 in the log–log and linear domain, respectively), but the semidiurnal IT variance is systematically smaller in HYCOM than in the mooring data. In Fig. 9b and c we show the geographical location of the patches along with a proxy for the HYCOMtomooring semidiurnal IT variance ratio and the histogram of this proxy, respectively. The latter histogram shows that the distribution of the proxy is centered around 0.37, corresponding to a ratio of 0.6, approximately.
To measure the strength of the decorrelation affecting the Eulerian mean autocovariances, we define the semidiurnal coherent variance fraction (SCVF_{15}) as the ratio between the 48 h demodulate at τ≈15 d (i.e., one spring–neap period) and the first demodulate. As it is calculated from the demodulate at relatively long time lags, the SCVF_{15} is meaningful only when estimated from a large enough sample size (i.e., with minimal uncertainty). Therefore, we start by considering all the data available over our global collection of geographical bins.
Figure 10a displays the globalmean autocovariances calculated from the mooring (${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{Moor}}$, solid red curve) and HYCOM data (${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}$, solid black curve), as the average of the localmean autocovariance series over the bins presented in Fig. 9b. In Fig. 10b we show different statistics of the observed distribution of the demodulates of the localmean autocovariances in the form of boxplots as a function of time lag. Figure 10b suggests that the mooring records exhibit both a larger semidiurnal IT variance and a stronger decorrelation on average globally: the mean of the first demodulates and the mean SCVF_{15} are 76 m^{2} and 0.47 as well as and 41 m^{2} and 0.57 for the moorings and HYCOM data, respectively. Using median values instead leads to the same conclusions (not shown).
We investigate a potential latitudinal dependence by plotting the mean autocovariance series computed separately from the geographical bins lying north (149 instruments) and south (11 instruments) of 50^{∘} S (see Figs. 10c and d, and 10e and f, respectively). The mean semidiurnal IT variance and the mean SCVF_{15} for each group are shown in Table 2. Again, the divergence between the two datasets appears enhanced south of 50^{∘} S. Moreover, the SCVF_{15} is larger in the Southern Ocean than elsewhere for HYCOM but smaller for the moorings. This indicates a weaker, or slower, decorrelation of the IT in HYCOM in the Southern Ocean than in the global ocean.
The slower decorrelation of the IT in the simulation can be explained by some decorrelating processes, such as eddies or submesoscale variability, being weaker in HYCOM than in the real ocean. It could also be explained by the time variability in certain decorrelating processes. Our numerical simulation only spans 20 May to 20 June 2019. Therefore, it is potentially missing processes that would specifically occur or intensify at another time of the year (or in a different year). On the other hand, the mooring data span several decades. Thus, our single month of data from HYCOM may not be representative of the broader temporal sampling of the mooring data.
The spatial distribution of the moorings is sparse and tends to be denser in particular areas (e.g., the Gulf Stream region). This is all the more true south of 50^{∘} S, where much fewer moorings are available than in the northern region and where they are mostly located in the Drake Passage. For both these reasons, the mean autocovariances computed from moorings north and especially south of 50^{∘} S cannot be considered truly representative of the IT in these vast regions. Nonetheless, they remain representative of the collection of geographical bins used to construct them.
4.3 Apparent decorrelation
In Sect. 4.1 we demonstrated that the motion of the floats has no significant impact on the measured total variance of the IT. However, the Doppler effect and spatial decorrelation induced by the Lagrangian sampling of the floats both act as decorrelating processes causing the autocovariance to decay with increasing time lag (Geoffroy and Nycander, 2022). Following CasparCohen et al. (2022), we call this mechanism apparent decorrelation, as it is unrelated to the decorrelation of the IT. Geoffroy and Nycander (2022) estimated the apparent decorrelation timescale to be longer than that of the autocovariance function observed by Argo floats, on average globally. Thus, they concluded that the decay of the globalmean autocovariance observed by Argo floats is primarily due to the decorrelation of the IT.
The HYCOM data allow this to be studied by directly comparing the globalmean autocovariance computed in the Lagrangian and Eulerian frameworks (see Fig. 11a and b). As expected, the Lagrangian (${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$, solid black curve in Fig. 11a) and Eulerian (${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$, solid red curve) globalmean autocovariances are virtually identical until τ≈50 h. Past this limit, the more rapid decay of ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$ can only be caused by the apparent decorrelation, while ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$ continues to solely reflect the decorrelation of the IT. Our aim in the present section is to find estimates for the characteristic timescales of the apparent decorrelation and decorrelation of the IT in HYCOM and compare them with observations.
Taking inspiration from Geoffroy and Nycander (2022) and CasparCohen et al. (2022), we try a simple model for the complex demodulate at the semidiurnal frequency of the Eulerian autocovariance computed over a time lag window centered on τ:
Here, ${\mathit{\sigma}}_{\mathrm{SD}}^{\mathrm{2}}$ is the total semidiurnal internal tide variance; α is the stationary fraction; T is the IT decorrelation timescale; and ${\mathit{\sigma}}_{\mathrm{AM}}^{\mathrm{2}}$ and ω_{AM} are the variance and angular frequency of an amplitudemodulating sinusoid, respectively. A heuristic model for the demodulate of the corresponding Lagrangian autocovariance is then
where T_{app} is the apparent decorrelation timescale.
A constrained least squares fit of the model Eq. (5) to the complex demodulate series computed from ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$ (red crosses in Fig. 11b) yields T=94 h. The fitted model and the exponential decay due to the IT decorrelation (i.e., only the term proportional to ${\mathit{\sigma}}_{\mathrm{SD}}^{\mathrm{2}}$) are represented in Fig. 11b as the solid and dashed red curves, respectively. On average globally, 95 % of the decorrelation of the nonstationary IT in HYCOM is therefore achieved within 3T≈300 h (for $\mathrm{exp}(\mathrm{3})\approx \mathrm{0.05}$), i.e., well under the 32 d of data.
By dividing ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$ by ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{E}}$, we can isolate the effects of the apparent decorrelation on ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{HYCOM}}^{\mathrm{L}}$. In Fig. 11c we plot this ratio after applying a median filter with a window of 18 h (solid black curve). A least squares fit of a simple decaying exponential to the obtained curve yields an apparent decorrelation timescale T_{app}=209 h (dashed black curve). For verification, we compute the Lagrangian ${\mathcal{C}}_{\mathrm{SD}}^{\mathrm{L}}$ from Eq. (6) by multiplying the fitted model of ${\mathcal{C}}_{\mathrm{SD}}^{\mathrm{E}}$ (solid red curve in Fig. 11b) by the fitted exponential decay due to the apparent decorrelation (dashed black curve in Fig. 11c). The result (solid black curve in Fig. 11b) closely follows the demodulated Lagrangian autocovariance (black crosses in Fig. 11b). The different parameters obtained from the Eulerian and Lagrangian simulated data are summarized in Table 3. Note that the fitted amplitudemodulating sinusoid is close to the spring–neap cycle, here ${\mathit{\omega}}_{\mathrm{AM}}={\mathit{\omega}}_{{\mathrm{M}}_{\mathrm{2}}}{\mathit{\omega}}_{{\mathrm{S}}_{\mathrm{2}}}\approx \mathrm{0.0177}$ h^{−1}.
We can compare the globalmean values of T and T_{app} obtained from the simulation with the Argo observations. The geographical coverage of the Argo data is less than that of HYCOM (roughly 55 % in area). Hence, instead of the global collections of floats and Lagrangian particles, we now consider the intersection of the two datasets taken as the unmasked bins shown in Fig. 7b and the corresponding localmean autocovariances. These localmean autocovariances are averaged further to obtain a globalmean autocovariance for each dataset that is representative of a common area. We then fit the model Eqs. (6) and (5) to the demodulates of the Argo globalmean autocovariance (see Fig. 12, updated from Geoffroy and Nycander, 2022). The parameters fitted to the Argo observations, as well as the parameters obtained by repeating the above procedure for the corresponding simulated data, are gathered in Table 4. The values of T_{app} and T computed here, from Argo data only, are similar to those reported by Geoffroy and Nycander (2022), where the authors estimated them from a comparable collection of Argo floats and HRET (in particular, the stationary limit was determined from HRET instead of by fitting). The fitted stationary fraction α, however, is about 3 times larger than previously reported. This could be explained by a lowbiased stationary variance (obtained by projecting HRET at 1000 dbar) in Geoffroy and Nycander (2022). As from the simulated data, fitting our model to the Argo globalmean demodulates yield a T shorter than T_{app}. While the Argo T_{app} is almost identical to the HYCOM value, T is about twice as small.
Similarly to Geoffroy and Nycander (2022), we conclude that the decorrelation of the IT is more rapid than the apparent decorrelation on average globally. Yet, according to the simulated data, it is only more rapid by a factor of 2, while this factor is 4 according to the Argo data. Hence, some decorrelating processes appear to be weaker in the simulation than in the real ocean. As explained in Sect. 4.2, the slower decorrelation of the IT in HYCOM could also be explained by certain decorrelating processes being weaker than usual in the 20 May–20 June 2019 period of outputs used here. Nevertheless, the decorrelation of the IT typically is at least as important as the apparent decorrelation over the first few days of time lag. This result might not hold true everywhere, since the geographical variability in T and α is expected to be large.
Lastly, Zaron (2022) provides a valuable comparison point for the results above. Along the same lines as in this section, he uses the spatially averaged frequency spectrum of sea level observations from satellite altimetry to measure properties of both the baroclinic tidal peak and continuum (i.e., the coherent and incoherent IT variability, respectively). For the semidiurnal band, in the latitude range from 30^{∘} S to 30^{∘} N, the author found that about 62 % of the IT variance eventually becomes decorrelated within 34 d. Assuming an exponential decay as in the model Eq. (5), this would represent a stationary fraction α≈0.38 and a decorrelation timescale T≈163 h (for 99 % of the nonstationary IT to decorrelate in 5T). Repeating the fitting procedure above for data limited to the ±30^{∘} latitude range, we recover values in reasonable agreement from the Argo observations: α≈0.29 and T≈144 h (see Table 4). The HYCOMfitted α≈0.64 and T≈105 h agree less well with the results of Zaron (2022) but still are in the same ballpark. In this latitude range, the bins used for the fitting represent roughly 45 % of the oceanic area.
4.4 Potential sources of bias
Why is the IT variance lower and IT decorrelation weaker in HYCOM than in the observations, particularly in the Southern Ocean? At the time of writing we cannot think of a particular reason for either the Argo or the mooringderived IT variance to be biased high globally. In particular, the correction accounting for the nontidal variability we systematically subtract from the sample autocovariance precludes any contamination of the first demodulate by the background noise (see Appendix C).
We also investigated whether the bias in the Southern Ocean could be related to the contamination of the first 48 h demodulate at ω_{SD} by nearinertial waves as we approach the M_{2} critical latitude (where $f={\mathit{\omega}}_{{\mathrm{M}}_{\mathrm{2}}}$, at about 74^{∘} S). To check this, we can map the semidiurnal IT variance from the Argo data anew (as shown in Fig. 6c), this time fitting an additional cos (fτ), where f is the local Coriolis frequency, to the localmean autocovariances. The result of the fit becomes unstable at 74^{∘} S, but the zonal mean of the demodulates at f does reach a maximum around 60^{∘} S, while the zonal mean of the demodulates at ω_{SD} remains unaffected (not shown). We conclude that the first 48 h demodulate at the semidiurnal frequency is not affected by nearinertial waves at any latitude.
As for the model, the horizontal grid spacing limits the number of vertical modes correctly resolved to the first five modes, equatorward of ±50^{∘} latitude and for seafloor depths deeper than 1250 m (Buijsman et al., 2020). Approaching the poles, the number of model layers below the mixed layer, and hence the vertical resolution, decreases. This further limits the number of resolved modes in HYCOM south of 50^{∘} S: roughly, modes 3 and 2 are only partially resolved poleward of 60 and 65^{∘} S, respectively. Although the bulk of the IT variance at 1000 dbar is likely related to mode1 waves on average globally (Geoffroy and Nycander, 2022), the contribution from higher modes can become significant locally. In principle, Argo floats, moorings, and HYCOM data include the effect of higher modes, but these are less well resolved in HYCOM, particularly in the Southern Ocean. Together with the assumption that highermode ITs are less coherent (Egbert and Ray, 2017), this may explain the lower IT variance and weaker IT decorrelation in HYCOM than in the observations, both on average globally and more specifically in the Southern Ocean. It is also in line with the larger mean SCVF_{15} computed from HYCOM data in the Southern Ocean compared with the rest of the globe (indicating a weaker decorrelation there), whereas for moorings the SCVF_{15} is smaller in this region (see Table 2).
The modem vertical structure of the isopycnal displacement Φ_{m}(z) is obtained by solving the Sturm–Liouville problem
with the boundary conditions ${\mathrm{\Phi}}_{m}\left(\mathrm{0}\right)={\mathrm{\Phi}}_{m}(H)=\mathrm{0}$, where H is the ocean depth, N(z) is the buoyancy frequency profile, and $\mathrm{1}/{c}_{m}^{\mathrm{2}}$ is the eigenvalue corresponding to the eigenfunction Φ_{m}(z) for modem (Gill, 1982). The modal partitioning of the IT energy at a given location is mainly determined by the conversion rate (both local and remote) and lifetime of each mode (de Lavergne et al., 2019). Additionally, depending on the local stratification and ocean depth, the parking depth at 1000 dbar can be more or less close to the antinode (point of maximal displacement) and node (point of no displacement) of the different vertical modes. Therefore, the normalized contribution of modem relative to mode1 waves to the variance recorded at this depth is weighted by a coefficient γ_{m1}. In Appendix B of Geoffroy and Nycander (2022), the authors derived an expression for this coefficient:
In Fig. 13a and b we plot the global maps of γ_{21} and γ_{31} at 1000 dbar, respectively, computed from Eq. (8) after solving the eigenvalue problem Eq. (7) for the $\mathrm{1}/\mathrm{4}$^{∘} WOA18 summer climatology. Here, we extrapolated the climatological data down to the reference bathymetry from the General Bathymetric Chart of the Oceans 2019 (GEBCO, 2019) wherever the bathymetry is deeper than the deepest valid climatological record. A visual comparison with Fig. 7b in the Southern Ocean suggests that low HYCOMtoArgo IT variance ratios (darkerblue pixels in Fig. 7b) spatially coincide with large γ_{31} mostly or, to some extent, γ_{21} (especially in the Weddell Sea). South of 60^{∘} S, either γ_{31} or γ_{21} is larger than unity; hence the contribution from mode3 or mode2 waves to the IT variance recorded at 1000 dbar is magnified compared to the contribution from mode1 waves. However, modes 2 and higher are not well resolved in HYCOM at these latitudes.
The magnification of the contributions from modes 2 and 3 to the variance at 1000 dbar in the Southern Ocean only affects how propagating ITs are perceived at the Argo parking depth. This has no connection with the generation and dissipation processes that set the underlying modal partitioning of the IT energy. Additional explanations might be found by examining whether the main parameters affecting the generation of ITs (namely the bottom topography, barotropic tidal forcing, and stratification) are less accurate in this region than in the rest of the globe.
The pattern of the enhanced discrepancies between Argo and HYCOM in the Southern Ocean (darker blue in Fig. 7) might be correlated with the distribution of bathymetric features. For instance, the large values around 45^{∘} S, 105^{∘} W, in the South Pacific are centered on the Chile Rise, a known IT generation area. To date, only about 19 % of the global ocean seafloor has been mapped using shipborne techniques. Therefore, global bathymetric products widely rely on depth prediction from satellite gravimetry (Smith and Sandwell, 1994). Smallscale features such as abyssal hills are not resolved by this technique.
In the recent SRTM15+ bathymetry (Tozer et al., 2019), gravitypredicted depths were estimated to have root mean square (RMS) uncertainties and maximum error of the order of ±150 and 1800 m, respectively (based on 50 cruises distributed globally). With about 22 % seafloor coverage by direct highquality measurements, the International Bathymetric Chart of the Southern Ocean v2 (IBCSO; Dorschel et al., 2022) is the stateoftheart bathymetric chart south of 50^{∘} S. A cellbycell comparison between IBCSO v2 and SRTM15+ v2.2 showed marked disparities, in particular for water depths between −4000 and −1000 m with differences reaching up to 1700 m (Dorschel et al., 2022). Most of the bathymetric errors only affect the generation of highmode ITs, which are bound to dissipate locally. Note that HYCOM may not resolve these highmode ITs anyway, as discussed in the present section. Still, the less frequent errors of the order of thousands of meters are likely to affect the generation of lowmode ITs in the model.
Efforts have been made to improve the accuracy of the M_{2} barotropic tides embedded in HYCOM. The simulation used in this study incorporates the framework presented in Ngodock et al. (2016) aiming at minimizing the tidal elevation RMS errors with respect to TPXO, a state of the art dataassimilative tide model (Egbert and Erofeeva, 2002). However, the tidal elevations in TPXO itself also have errors. Both Stammer et al. (2014) and Zaron and Elipot (2021) point at the imperfection of the modeled tidal elevations near Antarctica (using data from GRACE and CryoSat2, respectively). On the other hand, it is not so much the sea surface height that matters here, but the tidal currents. Zaron and Elipot (2021) used surface drifter observations for evaluating TPXOpredicted tidal currents throughout much of the global oceans. Unfortunately, the spatial density of observations is too poor to evaluate the model performance around much of Antarctica.
Lastly, to assess the stratification in HYCOM, we compare the phase speed of a mode1 gravity wave in the model with the phase speed determined from climatology. The phase speed of a modem gravity wave c_{m} is directly related to the eigenvalue $\mathrm{1}/{c}_{m}^{\mathrm{2}}$ obtained by solving the Sturm–Liouville problem (Eq. 7). Having already solved Eq. (7) for the climatology, we now solve it for our simulated monthly mean 3D fields of temperature and salinity (subsampled to 1^{∘}). As for the climatology, the HYCOM data were extended down to the reference bathymetry from GEBCO 2019, wherever deeper, by appending the stored bottom value. We then linearly interpolated the climatological phase speed of a mode1 gravity wave ${c}_{\mathrm{1}}^{\mathrm{WOA}}$ at $\mathrm{1}/\mathrm{4}$^{∘} resolution onto the coarser 1^{∘} HYCOM grid.
In Fig. 14a we plot the climatologytoHYCOM phase speed ratio. Most of the visible differences fall in the range of expected interannual variability of less than 10 % (Chelton et al., 1998). However, in a few areas around Antarctica there are larger departures of the ratio from unity. Figure 14b shows the zonal mean of the mode1 phase speed from Argo and HYCOM as a function of latitude. Equatorward of ±60^{∘}, the two curves agree almost perfectly, and they also visually agree with the results of Chelton et al. (1998) (not shown). Poleward of 60^{∘} S and 60^{∘} N, however, differences steadily grow. In the Southern Ocean, the zonalmeanclimatologytoHYCOM phase speed ratio peaks at 1.7 (see Fig. 14c). Typically, weaker stratification (smaller phase speed) results in smaller energy conversion.
In this work we compare a 32 d segment of a global run of the HYCOM model, including realistic tidal and atmospheric forcing, with in situ observations of the semidiurnal IT around 1000 dbar. First, a Lagrangian sampling of the simulation was compared to parkphase data from Argo floats to validate the geographical variability in the semidiurnal IT variance in HYCOM (see Sect. 4.1). Then, the Eulerian simulation outputs were directly compared to geographically sparser mooring records, in terms of variance and decorrelation of the IT (see Sect. 4.2).
The main spatial patterns of the simulated IT variance at 1000 dbar broadly agree with Argo observations, with energy radiating away from lowmode IT generation hotspots (see Fig. 6). Nonetheless, on average globally, the HYCOM data exhibit a smaller semidiurnal IT variance than observed by Argo floats, by a factor 0.58 (see Table 1). This is in line with the results of Ansong et al. (2017) and Luecke et al. (2020), who found HYCOM values to be biased low by a similar factor when comparing the simulated M_{2} IT energy flux and temperature variance to historical mooring observations.
While the difference between the model and Argo data appears reasonably homogeneous across most of the world ocean, it steadily increases towards the poles (see Fig. 7). Because of the scarcity of Argo floats available in the northernmost region, we focused on the Southern Ocean. On average, south of 50^{∘} S, we found that the simulated semidiurnal IT variance is smaller than the variance observed by Argo by a factor of 0.29. North of 50^{∘} S, this factor is 0.69 (see Table 1).
The mooring data support the above results for the semidiurnal IT variance. Additionally, we found that the decorrelation affecting the semidiurnal IT in HYCOM over a 32 d window is weaker than observed in the mooring records, on average (see Fig. 10 and Table 2). In other words, over timescales shorter than 32 d, the semidiurnal IT in HYCOM is more coherent than in observations. This weaker decorrelation of the IT in HYCOM can be explained by some decorrelating processes being weaker in the model than in the real ocean. It could also be explained by certain decorrelating processes being unusually weak in the 20 May–20 June 2019 period of outputs used in this work. Depending on the location, the complete decorrelation of the nonstationary IT is not systematically observable in a 32 d duration. Longer time series are thus needed to accurately describe the decorrelation of the IT. Nonetheless, we found that the semidiurnal IT autocovariance in HYCOM actually reaches its stationary limit within approximately 300 h on average globally, i.e., well under our 32 d of simulated data (see Fig. 11). Put together, these results support the conclusions of Buijsman et al. (2020), who found that the stationary (i.e., the longterm coherent) M_{2} IT solution from HYCOM is too energetic compared with altimetry. Note that their comparison was based on 1year simulated time series corrected for the duration difference with 17yearlong altimetry records.
We also investigated the effects of the Lagrangian sampling inherent to the Argo floats. When comparing autocovariances computed from the HYCOM data sampled in the Lagrangian and Eulerian frameworks, respectively, we found the IT variance to be unaffected in the mean (see Fig. 5). Moreover, the simulated apparent decorrelation (the decorrelation due to the motion of a Lagrangian particle) agrees very well with the apparent decorrelation experienced by Argo floats, on average globally (see Sect. 4.3). The (Eulerian) decorrelation of the IT in the simulated data, on the other hand, typically is half as rapid as the one inferred from Argo observations. This would make the IT decorrelation at least as important as the apparent one over the first few days of time lag. However, the geographical variability in the duration and strength of the IT decorrelation is expected to be large. Limiting the comparison to the latitude range from 30^{∘} S to 30^{∘} N, the Argo and HYCOM data lead to an IT decorrelation timescale of about 4.5 and 6 d, respectively (see Table 4), in reasonable agreement with the results of Zaron (2022).
Finally, we discuss the potential sources of bias. We could not think of a particular reason for the IT variance obtained from either the Argo or the mooring data to be biased high, particularly in the Southern Ocean. However, HYCOM is subject to various limitations. First and foremost, the model can only correctly resolve vertical modes up to 5 in most of the global oceans. Approaching the poles, the reduced number of layers further limits the number of resolved modes. While mode1 ITs supposedly account for most of the tidal variability at 1000 dbar on average globally (Geoffroy and Nycander, 2022), in situ instruments also record the contribution from higher modes that can become significant locally. This may explain why the IT variance is larger in the in situ data than in HYCOM, particularly in the Southern Ocean. Insufficient model stratification also seems to be a specific problem in that very region. We have not been able to quantitatively explain the overall smaller IT variance in HYCOM than in the in situ data over the global ocean. In principle, it could be due to limitations in the bathymetry, barotropic tidal forcing, or stratification.
We consider the function R(τ) on a 48 h interval. This is fitted to the function Acos (ωτ)+Bsin (ωτ), and we then define
This procedure is a mapping from the space of functions R to the positive number C, i.e., a functional. We can write this symbolically as
where Φ is the functional.
Φ obviously has the property
where a is a constant. Φ also has the property (proof hereafter)
with equality only if R_{1}(τ) and R_{2}(τ) are proportional, i.e., perfectly correlated with positive correlation.
Let [R_{i}(τ)] be an infinite series of functions drawn from a stochastic distribution. We have access to only N of these (typically around 10). Define the sample average
${\stackrel{\mathrm{\u203e}}{R}}_{N}\left(\mathit{\tau}\right)$ is an unbiased estimate of the true average $\stackrel{\mathrm{\u203e}}{R}\left(\mathit{\tau}\right)$:
We also want an estimate of $\mathrm{\Phi}\left[\stackrel{\mathrm{\u203e}}{R}\right]$. Denote the true value by C_{∞}:
A sample estimate might be
However, Eq. (A4) implies that this is biased. For example, define
with equality only if ${\stackrel{\mathrm{\u203e}}{R}}_{N}$ and ${\stackrel{\mathrm{\u203e}}{R}}_{M}$ are proportional. But ${\stackrel{\mathrm{\u203e}}{R}}_{N}$ and ${\stackrel{\mathrm{\u203e}}{R}}_{M}$ are stochastic functions drawn from the same distribution and therefore are almost never proportional. On average we also have $\mathrm{\Phi}\left[{\stackrel{\mathrm{\u203e}}{R}}_{N}\right]=\mathrm{\Phi}\left[{\stackrel{\mathrm{\u203e}}{R}}_{M}\right]$; hence Eq. (A10) implies
so that C_{N} is a decreasing function of N. Thus, we expect C_{N} to be larger than C_{∞}. In other words, the complex demodulate of the sample function ${\stackrel{\mathrm{\u203e}}{R}}_{N}\left(\mathit{\tau}\right)$, i.e., the fitting defined in Eq. (A2), is a highbiased estimate of the envelope of the true function $\stackrel{\mathrm{\u203e}}{R}\left(\mathit{\tau}\right)$.
We now show the inequality Eq. (A4). It resembles the triangle inequality for two vectors:
For 2D vectors $\mathit{a}=({a}_{\mathrm{1}},{a}_{\mathrm{2}})$, $\mathit{b}=({b}_{\mathrm{1}},{b}_{\mathrm{2}})$, this gives
which is valid for arbitrary numbers a_{1}, a_{2}, b_{1}, and b_{2}.
Define a scalar product 〈f,g〉 between the functions f(τ) and g(τ) as the integral of fg over the 48 h interval:
Suppose that the period of sin (ωτ) exactly fits this window (note that this is the case in the present work, where $\mathrm{4}\times \mathrm{2}\mathit{\pi}/{\mathit{\omega}}_{\mathrm{SD}}\approx \mathrm{48.8}$ h). Then sin (ωτ) and cos (ωτ) are orthogonal:
Any function f(τ) can be expanded in a Fourier series on the interval. A least squares fit of f(τ) to Acos (ωτ)+Bsin (ωτ) is the same thing as projecting f(τ) on the basis functions cos (ωτ) and sin (ωτ). Thus, the best fit is given by
where $E=\langle \mathrm{cos}\left(\mathit{\omega}\mathit{\tau}\right),\mathrm{cos}\left(\mathit{\omega}\mathit{\tau}\right)\rangle =\langle \mathrm{sin}\left(\mathit{\omega}\mathit{\tau}\right),\mathrm{sin}\left(\mathit{\omega}\mathit{\tau}\right)\rangle $. Consider two functions f(τ) and g(τ), and denote
We then have
Two aspects of the processing originally used in Geoffroy and Nycander (2022) are important in mitigating a potential contamination of the first demodulate by the nontidal variability (noise). (i) The demodulation over a given 48 h window is performed in two iterations. Acos (ω_{SD}τ)+Bsin (ω_{SD}τ) is first fitted to the signal, and the root mean square error (RMSE) is computed. The fitting is then repeated with a trimmed signal that excludes values outside of ±2 RMSE from the previously fitted curve. (ii) The time series are highpassfiltered prior to computing their sample autocovariance. Suppose the nontidal variability ρ(t) is a firstorder autoregressive process (AR1) characterized by zero mean, the timescale τ_{ρ}, and the white Gaussian noise ε_{ρ}(t) with zero mean and variance ${\mathit{\sigma}}_{\mathit{\epsilon},\phantom{\rule{0.125em}{0ex}}\mathit{\rho}}^{\mathrm{2}}$. The variance and autocovariance of an AR1 process are given by
respectively. For very short τ_{ρ}, ρ(t) is close to a white noise, i.e., showing a flat power spectrum. As τ_{ρ} increases, ρ(t) resembles a red noise with increasingly steep exponential decay in spectral space. Therefore, as τ_{ρ} increases, an increasing fraction of ${\mathit{\sigma}}_{\mathit{\rho}}^{\mathrm{2}}$ is filtered out by the highpass filter.
We can investigate this further using a simple model, adapted from Geoffroy and Nycander (2022), of a tidal variability on top of a background noise:
Here, ω_{i} and A_{i} are the angular frequency and the amplitude of the tidal constituent i, respectively, where $i\in [{M}_{\mathrm{2}},\phantom{\rule{0.125em}{0ex}}{S}_{\mathrm{2}}]$. ϕ(t) and ρ(t) are AR1 processes representing random phase modulations and a nontidal variability, respectively. These AR1 processes can be defined by their variance and characteristic time. We vary ${\mathit{\sigma}}_{\mathit{\rho}}^{\mathrm{2}}$ and τ_{ρ} (the variance and characteristic time of ρ(t), respectively) with the rest of the variables set to realistic values: the tidal variance ${\mathit{\sigma}}_{\mathrm{SD}}^{\mathrm{2}}=({A}_{{M}_{\mathrm{2}}}^{\mathrm{2}}+{A}_{{S}_{\mathrm{2}}}^{\mathrm{2}})/\mathrm{2}=\mathrm{40}$ m^{2}, and ${\mathit{\sigma}}_{\mathit{\varphi}}^{\mathrm{2}}=\mathrm{1.19}$ rad^{2} and τ_{ϕ}=150 h (the variance and characteristic time of ϕ(t), respectively). Since the tidal variance is fixed, varying ${\mathit{\sigma}}_{\mathit{\rho}}^{\mathrm{2}}$ is the same as varying the signaltonoise ratio $S/N={\mathit{\sigma}}_{\mathrm{SD}}^{\mathrm{2}}/{\mathit{\sigma}}_{\mathit{\rho}}^{\mathrm{2}}$. For reference, in Appendix C we estimate the global distribution of the signaltonoise ratio observed by Argo floats (not shown). This distribution is skewed, with a median of 0.8 and a 5th percentile of 0.2, approximately. We also estimate the median τ_{ρ}≈1 h from the global collection of floats.
We start by varying τ_{ρ} for a fixed variance value of ${\mathit{\sigma}}_{\mathit{\rho}}^{\mathrm{2}}=\mathrm{200}$ m^{2}, corresponding to the extreme $S/N=\mathrm{0.2}$. For a given value of τ_{ρ}, we generate one thousand 32 d long synthetic time series h(t) following the model Eq. (B2). The time series are highpassfiltered using a fourthorder Butterworth filter with a cutoff frequency of 0.3 cpd. We then compute a sample autocovariance from each filtered time series and average these to obtain a mean autocovariance. Finally we compute the 48 h demodulates of the mean autocovariance following the complex demodulation method described in Sect. 3.2.
In Fig. B1, we show the mean autocovariance and the corresponding 48 h demodulates for τ_{ρ}=0.1, 1, and 12 h. The value of the first demodulate for τ_{ρ}=0.1 is essentially the same as if there were no ρ(t) at all, namely $C(\mathit{\tau}=[\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{48}\left]\right)\approx \mathrm{34}$ m^{2}. For longer τ_{ρ}, the main difference with the previous curve can be easily visualized from the first few demodulates. In particular, the first demodulates $C(\mathit{\tau}=[\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{48}\left]\right)\approx \mathrm{38}$ and 31 m^{2}, corresponding to τ_{ρ}=1 and 12 h, are systematically biased high and low, respectively.
We repeat the same experiment for τ_{ρ} ranging from 0 to 48 h and focusing on the effects of the noise on the first demodulate. In Fig. B2, we plot the first demodulate as a function of τ_{ρ} (solid curve). For any τ_{ρ}, the first demodulate is seen to remain lower than the true ${\mathit{\sigma}}_{\mathrm{SD}}^{\mathrm{2}}=\mathrm{40}$ m^{2}. Moreover, two regimes can be identified: (i) for τ_{ρ} shorter than 5 h, roughly, the contribution of the nontidal variability to the first demodulate is positive and peaks for τ_{ρ}=1 h, where $C(\mathit{\tau}=[\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{48}\left]\right)\approx \mathrm{38}$ m^{2}, and (ii) for longer τ_{ρ}, this contribution is negative, with a minimum $C(\mathit{\tau}=[\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{48}\left]\right)\approx \mathrm{30}$ m^{2} at τ_{ρ}=11 h. For a more typical $S/N=\mathrm{0.8}$, i.e., setting ${\mathit{\sigma}}_{\mathit{\rho}}^{\mathrm{2}}=\mathrm{50}$ m^{2} in our synthetic time series, we found a similar dependence of the first demodulate on τ_{ρ} as previously but with a much smaller span in amplitude (dashdotted curve in Fig. B2). The contribution of the nontidal variability to the first demodulate again peaks for τ_{ρ}=1 h, with $C(\mathit{\tau}=[\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{48}\left]\right)\approx \mathrm{35}$ m^{2}, and the minimum occurs around τ_{ρ}=8 h, where $C(\mathit{\tau}=[\mathrm{0},\phantom{\rule{0.33em}{0ex}}\mathrm{48}\left]\right)\approx \mathrm{33}$ m^{2}.
Hence, for the typical $S/N=\mathrm{0.8}$ and τ_{ρ}=1 observed by Argo floats, we estimate that the nontidal variability can lead to a first demodulate being biased high by roughly 3 %. For the same characteristic timescale but $S/N=\mathrm{0.2}$ (the 5th percentile of the global $S/N$ distribution sampled by Argo floats), this bias can reach 10 %. Notwithstanding, even for such an extreme $S/N$ value, the first demodulate is seen to remain a conservative estimate of the IT variance. The effects of the nontidal variability on the IT variance estimate as computed in Geoffroy and Nycander (2022) thus do not put their results into question.
It is not straightforward to understand why the contribution of the nontidal variability to the value of the first demodulate can be negative. Because of the way we defined ρ(t), i.e., an AR1 process with a positive parameter, its autocovariance should remain positive. Oscillations of R_{ρ}(τ) are actually a consequence of the highpass filter applied to the time series. We show this in Fig. B3, where we plot the theoretical autocovariance of ρ(t) with τ_{ρ}=12 h and ${\mathit{\sigma}}_{\mathit{\rho}}^{\mathrm{2}}=\mathrm{200}$ m^{2} (solid black curve) and the sample mean autocovariance from 1000 nonfiltered and 1000 highpassfiltered synthetic time series (solid and dashdotted red curve, respectively) constructed using the same parameters.
As mentioned in Sect. 3.2, we expect the stochastic background noise affecting the simulated data to be different from the one recorded by the in situ instruments. Therefore, to prevent any systematic bias in the comparisons presented in this work, we consistently correct the sample autocovariance by subtracting an estimate of the autocovariance of the nontidal variability before performing the complex demodulation (see Appendix C). We illustrate the effects of this correction in Fig. B4. Here, we proceed in the same way as for Fig. B2, but this time varying $S/N$ (from $\mathrm{4}\times {\mathrm{10}}^{\mathrm{2}}$ to 4) for fixed τ_{ρ}=1 and 12 h. While the first demodulates of the noncorrected mean autocovariances diverge for small $S/N$ values (solid curves), the first demodulates of the noisecorrected autocovariances are virtually constant (dashdotted curves).
In the present appendix, we derive a model for the autocovariance of a tidal variability on top of a highpassfiltered stochastic noise. In Appendix B we show that the autocovariance of the nontidal variability is affected by the highpass filter we apply on the original time series (see Fig. B3). This consequence of filtering the time series was not taken into account by Geoffroy and Nycander (2022) in the model of the background noise they used.
Assume that the nontidal variability ρ(t) is a firstorder autoregressive process (AR1) characterized by zero mean, the timescale τ_{ρ}, and the white Gaussian noise ε_{ρ}(t) with zero mean and variance ${\mathit{\sigma}}_{\mathit{\epsilon},\phantom{\rule{0.125em}{0ex}}\mathit{\rho}}^{\mathrm{2}}$. In the time domain, the filtered ρ^{′}(t) can be modeled as the convolution of the unfiltered ρ(t) with the impulse response of the filter h(t). We are, however, interested in the autocovariance of ρ^{′}(t), ${R}_{\mathit{\rho}}^{\prime}\left(\mathit{\tau}\right)$, which is closely related to the power spectrum (one is the Fourier transform of the other). Luckily, it is much simpler to perform this operation in the frequency domain, where the power spectrum at the output of a linear filter, ${\mathcal{P}}_{\mathit{\rho}}^{\prime}\left(z\right)$, is related to the power spectrum of the input stochastic process, 𝒫_{ρ}(z), by
Here, H(z) is the system or transfer function of the filter, and the z domain is limited to the unit circle z=exp (jωΔt); i.e., the z transform is equivalent to the discrete Fourier transform, where ω is the angular frequency, and Δt is the constant sampling interval (Kay and Marple, 1981). Moreover, the power spectrum of an AR1 process is a standard result:
The Butterworth filter is linear in amplitude but not in phase. As a workaround, we applied a secondorder filter twice, once forward and once backward, to recover a fourthorder filter with zero phase. The transfer function of a secondorder highpass Butterworth filter in the z domain, H_{B2}(z), can be obtained from its Laplace transform:
where ω_{c} is the cutoff angular frequency (Schlichthärle, 2000). We use the bilinear transform, namely setting
in Eq. (C3), to recover the expression
Hence, we can evaluate Eq. (C1) as
Since ρ^{′} is realvalued, the autocovariance ${R}_{\mathit{\rho}}^{\prime}\left(\mathit{\tau}\right)$ is given by the inverse discrete Fourier cosine transform of its power spectrum ${\mathcal{P}}_{\mathit{\rho}}^{\prime}\left(\mathit{\omega}\right)$ (computed numerically).
As mentioned in Sect. 3.2, to limit any systematic bias we consistently estimate and subtract the autocovariance of the nontidal variability from the sample autocovariance before performing the complex demodulation. We obtain such an estimate from the least squares fit of the model adapted from Geoffroy and Nycander (2022),
corresponding to the autocovariance of a tidal variability affected by a background noise:
Here, ω_{i} and A_{i} are the angular frequency and the amplitude of the tidal constituent i, respectively, where $i\in [{M}_{\mathrm{2}},{S}_{\mathrm{2}},{K}_{\mathrm{1}},{O}_{\mathrm{1}}]$. ϕ_{i}(t) is an AR1 process representing the random phase modulations affecting the constituent i. For the sake of simplicity, we assume that the decorrelating processes act in the same manner on neighboring tidal frequencies, that is, ${\mathit{\varphi}}_{{M}_{\mathrm{2}}}\left(t\right)={\mathit{\varphi}}_{{S}_{\mathrm{2}}}\left(t\right)$ and ${\mathit{\varphi}}_{{K}_{\mathrm{1}}}\left(t\right)={\mathit{\varphi}}_{{O}_{\mathrm{1}}}\left(t\right)$. Note that the highpass filter used in this work was designed so as not to affect the tidal signal. ${\mathit{\rho}}_{\mathrm{f}}^{\prime}\left(t\right)$ and ${\mathit{\rho}}_{\mathrm{s}}^{\prime}\left(t\right)$ are highpassfiltered AR1 processes accounting for a fast and a slow nontidal variability, respectively. All the AR1 processes are characterized by zero mean, the timescale τ_{X}, and the white Gaussian noise ε_{X}(t) with zero mean and variance ${\mathit{\sigma}}_{\mathit{\epsilon},\phantom{\rule{0.125em}{0ex}}X}^{\mathrm{2}}$, where $X\in [\mathit{\varphi},{\mathit{\rho}}_{\mathrm{f}},{\mathit{\rho}}_{\mathrm{s}}]$. The variance and autocovariance of an AR1 process are
respectively. The model Eq. (C7) does not take into account the effects of the drift of the floats or Lagrangian particles.
We fit the full model Eq. (C7), containing 12 parameters, to the sample autocovariance by nonlinear least squares, with weights inversely proportional to the corresponding SEM, and imposed bounds on the parameters (since we only fit positive timescales and variances). In particular, the timescales of the fast and slow nontidal variability are restricted to the [0, 5] and [5, 48] h range, respectively. The autocovariance of the nontidal variability is then computed as
using the fitted parameters. Lastly, ${R}_{\mathit{\rho}}^{\prime}\left(\mathit{\tau}\right)$ is subtracted from the sample autocovariance to correct for the effects of the filtered background noise.
The inclusion of a tidal variability in the model Eq. (C7) primarily aims at limiting the projection of the observed tidal variability onto our model for the background noise when fitting. Nonetheless, the fitted parameters can be used to compute an estimate of the semidiurnal IT variance ${\mathit{\sigma}}_{\mathrm{SD},\phantom{\rule{0.125em}{0ex}}\mathrm{LS}}^{\mathrm{2}}=({A}_{{M}_{\mathrm{2}}}^{\mathrm{2}}+{A}_{{S}_{\mathrm{2}}}^{\mathrm{2}})/\mathrm{2}$ and further to estimate the signaltonoise ratio characterizing the filtered Argo records as $S/N={\mathit{\sigma}}_{\mathrm{SD},\phantom{\rule{0.125em}{0ex}}\mathrm{LS}}^{\mathrm{2}}/{R}_{\mathit{\rho}}^{\prime}\left(\mathrm{0}\right)$. In Fig. C1 we show maps of ${\mathit{\sigma}}_{\mathrm{SD},\phantom{\rule{0.125em}{0ex}}\mathrm{LS}}^{\mathrm{2}}$, ${R}_{\mathit{\rho}}^{\prime}\left(\mathrm{0}\right)$, and $S/N$ computed from the same collection of localmean autocovariances as used in Fig. 6c. Note that ${\mathit{\sigma}}_{\mathrm{SD},\phantom{\rule{0.125em}{0ex}}\mathrm{LS}}^{\mathrm{2}}$ and ${R}_{\mathit{\rho}}^{\prime}\left(\mathrm{0}\right)$ are not correlated (r^{2}=0.13 and 0.009 in the log–log and linear domain, respectively). The leastsquaresfitted ${\mathit{\sigma}}_{\mathrm{SD},\phantom{\rule{0.125em}{0ex}}\mathrm{LS}}^{\mathrm{2}}$ appears reasonable when compared with the first demodulates presented in Fig. 6c. The coefficient of determination r^{2}=0.53 and 0.29 in the log–log and linear domain, respectively, indicates a good agreement between the two estimates. Moreover, the global mean and median of their ratio are 0.8 and 0.7, respectively, accounting for the decorrelation of the IT taking place in the first 48 h. For such localmean autocovariances, however, the uncertainties are simply to high to reliably estimate the parameters of the decorrelating process ϕ_{i}(t). Finally, we chose to perform the comparisons presented in this work in terms of the first demodulate, for it is a conservative and more robust estimate of the IT variance.
Argo data were obtained from US GDAC (ftp://usgodae.org/pub/outgoing/argo, last access: 11 January 2023, Argo, 2000). These data were collected and made freely available by the international Argo program and the national programs that contribute to it (https://argo.ucsd.edu, https://www.oceanops.org, last access: 26 May 2023). The Argo program is part of the Global Ocean Observing System. An Argo Iridium float list is maintained by Stephen C. Riser (http://runt.ocean.washington.edu/argo/heterographs/rollcall.html, Riser and Swift, 2023, last access: 11 January 2023). The code used to download and process the Argo data is available at https://doi.org/10.57669/geoffroy2022argoit1.0.0 (Geoffroy, 2022a). A global map of the total semidiurnal internal tide variance at 1000 dbar produced using the latter code is available at https://doi.org/10.17043/geoffroy2022argoit1 (Geoffroy, 2022b). There is no longterm availability plan for the HYCOM data used in this work. Climatological data were obtained from the World Ocean Atlas 2018 (https://accession.nodc.noaa.gov/NCEIWOA18). Bathymetric data were obtained from the GEBCO Compilation Group (https://www.gebco.net/data_and_products/gridded_bathymetry_data/). The Global MultiArchive Current Meter Database is not publicly available but can be obtained through request (http://stockage.univbrest.fr/~scott/GMACMD/gmacmd.html). NetCDF versions of the baroclinic tidal harmonic constants from the High Resolution Empirical Tide model are made available by Edward D. Zaron (https://ingria.ceoas.oregonstate.edu/~zarone/downloads.html).
GG: conceptualization, formal analysis, project administration, software, visualization, writing – original draft, writing – review and editing. JN: conceptualization, funding acquisition, writing – original draft, writing – review and editing. MCB: project administration, writing – review and editing. JFS: formal analysis, resources, software. BKA: writing – review and editing.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to thank Shane Elipot and one additional anonymous reviewer who greatly helped improve this paper.
Gaspard Geoffroy and Jonas Nycander are supported by grant number 201704623 from the Swedish Research Council. The computations and data handling were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre (NSC), partially funded by the Swedish Research Council through grant agreement no. 201805973. Maarten C. Buijsman and Brian K. Arbic are supported by the Office of Naval Research (ONR) grant numbers N000141912704 and N000141912712, respectively, which both fall under the project “Modeling, characterizing, and predicting effects of internal gravity waves on acoustic propagation on basin to global scales”. Jay F. Shriver was supported by ONR grant number N0001422WX01919, which falls under the project “Diagnosis and validation of the time and spatial variability of remotely generated internal waves in global ocean simulations”.
The article processing charges for this openaccess publication were covered by Stockholm University.
This paper was edited by Ilker Fer and reviewed by Shane Elipot and one anonymous referee.
Alford, M. H. and Zhao, Z.: Global patterns of lowmode internalwave propagation. Part I: Energy and energy flux, J. Phys. Oceanogr., 37, 1829–1848, https://doi.org/10.1175/JPO3085.1, 2007. a
Ansong, J. K., Arbic, B. K., Alford, M. H., Buijsman, M. C., Shriver, J. F., Zhao, Z., Richman, J. G., Simmons, H. L., Timko, P. G., Wallcraft, A. J., and Zamudio, L.: Semidiurnal internal tide energy fluxes and their variability in a global ocean model and moored observations, J. Geophys. Res.Oceans, 122, 1882–1900, https://doi.org/10.1002/2016JC012184, 2017. a, b
Arbic, B. K.: Incorporating tides and internal gravity waves within global ocean general circulation models: A review, Prog. Oceanogr., 206, 102824, https://doi.org/10.1016/j.pocean.2022.102824, 2022. a
Arbic, B. K., Wallcraft, A. J., and Metzger, E. J.: Concurrent simulation of the eddying general circulation and tides in a global ocean model, Ocean Model., 32, 175–187, https://doi.org/10.1016/j.ocemod.2010.01.007, 2010. a
Arbic, B. K., Elipot, S., Brasch, J. M., Menemenlis, D., Ponte, A. L., Shriver, J. F., Yu, X., Zaron, E. D., Alford, M. H., Buijsman, M. C., Abernathey, R., Garcia, D., Guan, L., Martin, P. E., and Nelson, A. D.: Frequency dependence of nearsurface oceanic kinetic energy from drifter observations and global highresolution models, arXiv [preprint], https://doi.org/10.48550/ARXIV.2202.08877, 22 July 2022. a
Argo: Argo float data and metadata from Global Data Assembly Centre (Argo GDAC), SEANOE [data set], https://doi.org/10.17882/42182, 2000. a
Boyer, T. P., Garcia, H. E., Locarnini, Ricardo A., Zweng, M. M., Mishonov, A. V., Reagan, J. R., Weathers, K. A., Baranova, O. K., Seidov, D., and Smolyar, I. V.: World Ocean Atlas 2018. Temperature, and salinity, NOAA National Centers for Environmental Information, [data set], https://accession.nodc.noaa.gov/NCEIWOA18 (last access: 11 October 2021), 2018. a
Buijsman, M. C., Arbic, B. K., Richman, J. G., Shriver, J. F., Wallcraft, A. J., and Zamudio, L.: Semidiurnal internal tide incoherence in the equatorial Pacific, J. Geophys. Res.Oceans, 122, 5286–5305, https://doi.org/10.1002/2016JC012590, 2017. a
Buijsman, M. C., Stephenson, G. R., Ansong, J. K., Arbic, B. K., Green, J. M., Richman, J. G., Shriver, J. F., Vic, C., Wallcraft, A. J., and Zhao, Z.: On the interplay between horizontal resolution and wave drag and their effect on tidal baroclinic mode waves in realistic global ocean simulations, Ocean Model., 152, 101656, https://doi.org/10.1016/j.ocemod.2020.101656, 2020. a, b, c, d, e, f
CasparCohen, Z., Ponte, A., Lahaye, N., Carton, X., Yu, X., and Gentil, S. L.: Characterization of internal tide incoherence: Eulerian versus Lagrangian perspectives, J. Phys. Oceanogr., 52, 1245–1259, https://doi.org/10.1175/JPOD210088.1, 2022. a, b, c, d, e
Chassignet, E. P., Hurlburt, H. E., Smedstad, O. M., Halliwell, G. R., Wallcraft, A. J., Metzger, E. J., Blanton, B. O., Lozano, C., Rao, D. B., and Hogan, P. J.: Generalized vertical coordinates for eddyresolving global and coastal ocean forecasts, Oceanography, 19, 118–129, https://doi.org/10.5670/oceanog.2006.95, 2006. a
Chelton, D. B., deSzoeke, R. A., Schlax, M. G., El Naggar, K., and Siwertz, N.: Geographical variability of the first baroclinic Rossby radius of deformation, J. Phys. Oceanogr., 28, 433–460, https://doi.org/10.1175/15200485(1998)028<0433:GVOTFB>2.0.CO;2, 1998. a, b
Cherniawsky, J. Y., Foreman, M. G. G., Crawford, W. R., and Henry, R. F.: Ocean tides from TOPEX/Poseidon sea level data, J. Atmos. Ocean. Tech., 18, 649–664, https://doi.org/10.1175/15200426(2001)018<0649:OTFTPS>2.0.CO;2, 2001. a
de Lavergne, C., Falahat, S., Madec, G., Roquet, F., Nycander, J., and Vic, C.: Toward global maps of internal tide energy sinks, Ocean Model., 137, 52–75, https://doi.org/10.1016/j.ocemod.2019.03.010, 2019. a
de Lavergne, C., Vic, C., Madec, G., Roquet, F., Waterhouse, A. F., Whalen, C. B., Cuypers, Y., BouruetAubertot, P., Ferron, B., and Hibiya, T.: A parameterization of local and remote tidal mixing, J. Adv. Model. Earth Sy., 12, e2020MS002065, https://doi.org/10.1029/2020MS002065, 2020. a
Dorschel, B., Hehemann, L., Viquerat, S., Warnke, F., Dreutter, S., Tenberge, Y. S., Accettella, D., An, L., Barrios, F., Bazhenova, E., et al.: The International Bathymetric Chart of the Southern Ocean version 2, Scientific Data, 9, 275, https://doi.org/10.1038/s41597022013667, 2022. a, b
Dushaw, B. D., Worcester, P. F., and Dzieciuch, M. A.: On the predictability of mode1 internal tides, DeepSea Res. Pt. I, 58, 677–698, https://doi.org/10.1016/j.dsr.2011.04.002, 2011. a
Egbert, G. D. and Erofeeva, S. Y.: Efficient inverse modeling of barotropic ocean tides, J. Atmos. Ocean. Tech., 19, 183–204, https://doi.org/10.1175/15200426(2002)019<0183:EIMOBO>2.0.CO;2, 2002. a
Egbert, G. D. and Ray, R. D.: Tidal prediction, J. Mar. Res., 75, 189–237, https://doi.org/10.1357/002224017821836761, 2017. a, b
GEBCO: GEBCO 2019 Grid, BODC [data set], https://doi.org/10.5285/836f016a33be6ddce0536c86abc0788e, 2019. a
Geoffroy, G.: Scripts for mapping the total semidiurnal internal tide variance at 1,000 dbar using Argo data, GitLab [code], https://doi.org/10.57669/geoffroy2022argoit1.0.0, 2022a. a
Geoffroy, G.: Global map of the total semidiurnal internal tide variance at 1,000 dbar from Argo data, Bolin Centre Database [data set], https://doi.org/10.17043/geoffroy2022argoit1, 2022b. a
Geoffroy, G. and Nycander, J.: Global mapping of the nonstationary semidiurnal internal tide using Argo data, J. Geophys. Res.Oceans, 127, e2021JC018283, https://doi.org/10.1029/2021JC018283, 2022. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag
Gill, A. E.: Atmosphereocean dynamics, Elsevier, Academic Press, ISBN 0122835204, 1982. a
Hennon, T. D., Riser, S. C., and Alford, M. H.: Observations of internal gravity waves by Argo floats, J. Phys. Oceanogr., 44, 2370–2386, https://doi.org/10.1175/JPOD130222.1, 2014. a, b, c
Kay, S. and Marple, S.: Spectrum analysis—A modern perspective, P. IEEE, 69, 1380–1419, https://doi.org/10.1109/PROC.1981.12184, 1981. a
Luecke, C. A., Arbic, B. K., Richman, J. G., Shriver, J. F., Alford, M. H., Ansong, J. K., Bassette, S. L., Buijsman, M. C., Menemenlis, D., Scott, R. B., Timko, P. G., Voet, G., Wallcraft, A. J., and Zamudio, L.: Statistical comparisons of temperature variance and kinetic energy in global ocean models and observations: Results from mesoscale to internal wave frequencies, J. Geophys. Res.Oceans, 125, e2019JC015306, https://doi.org/10.1029/2019JC015306, 2020. a, b, c
MacKinnon, J. A., Zhao, Z., Whalen, C. B., Waterhouse, A. F., Trossman, D. S., Sun, O. M., Laurent, L. C. S., Simmons, H. L., Polzin, K., Pinkel, R., Pickering, A., Norton, N. J., Nash, J. D., Musgrave, R., Merchant, L. M., Melet, A. V., Mater, B., Legg, S., Large, W. G., Kunze, E., Klymak, J. M., Jochum, M., Jayne, S. R., Hallberg, R. W., Griffies, S. M., Diggs, S., Danabasoglu, G., Chassignet, E. P., Buijsman, M. C., Bryan, F. O., Briegleb, B. P., Barna, A., Arbic, B. K., Ansong, J. K., and Alford, M. H.: Climate Process Team on internal wave–driven ocean mixing, B. Am. Meteorol. Soc., 98, 2429–2454, https://doi.org/10.1175/BAMSD160030.1, 2017. a
Melet, A., Legg, S., and Hallberg, R.: Climatic impacts of parameterized local and remote tidal mixing, J. Climate, 29, 3473–3500, https://doi.org/10.1175/JCLID150153.1, 2016. a
Melet, A., Hallberg, R., and Marshall, D. P.: Chapter 2 – The role of ocean mixing in the climate system, in: Ocean Mixing, edited by: Meredith, M. and Naveira Garabato, A., 5–34, Elsevier, https://doi.org/10.1016/B9780128215128.000098, 2022. a
Nelson, A. D., Arbic, B. K., Zaron, E. D., Savage, A. C., Richman, J. G., Buijsman, M. C., and Shriver, J. F.: Toward realistic nonstationarity of semidiurnal baroclinic tides in a hydrodynamic model, J. Geophys. Res.Oceans, 124, 6632–6642, https://doi.org/10.1029/2018JC014737, 2019. a, b
Ngodock, H. E., Souopgui, I., Wallcraft, A. J., Richman, J. G., Shriver, J. F., and Arbic, B. K.: On improving the accuracy of the M_{2} barotropic tides embedded in a highresolution global ocean circulation model, Ocean Model., 97, 16–26, https://doi.org/10.1016/j.ocemod.2015.10.011, 2016. a
Rainville, L. and Pinkel, R.: Propagation of lowmode internal waves through the ocean, J. Phys. Oceanogr., 36, 1220–1236, https://doi.org/10.1175/JPO2889.1, 2006. a
Ray, R. D. and Zaron, E. D.: M_{2} internal tides and their observed wavenumber spectra from satellite altimetry, J. Phys. Oceanogr., 46, 3–22, https://doi.org/10.1175/JPOD150065.1, 2016. a
Riser, S. C. and Swift, D.: University of Washington Argo profiling drifters, http://runt.ocean.washington.edu/argo/heterographs/rollcall.html, last access: 11 January 2023. a
Schlichthärle, D.: Digital filters: basics and design, Springer Berlin, Heidelberg, https://doi.org/10.1007/9783662041703, 2000. a
Scott, R. B., Goff, J. A., Naveira Garabato, A. C., and Nurser, A. J. G.: Global rate and spectral characteristics of internal gravity wave generation by geostrophic flow over topography, J. Geophys. Res.Oceans, 116, C09029, https://doi.org/10.1029/2011JC007005, 2011. a
Shriver, J. F., Arbic, B. K., Richman, J. G., Ray, R. D., Metzger, E. J., Wallcraft, A. J., and Timko, P. G.: An evaluation of the barotropic and internal tides in a highresolution global ocean circulation model, J. Geophys. Res.Oceans, 117, C10024, https://doi.org/10.1029/2012JC008170, 2012. a
Shriver, J. F., Richman, J. G., and Arbic, B. K.: How stationary are the internal tides in a highresolution global ocean circulation model?, J. Geophys. Res.Oceans, 119, 2769–2787, https://doi.org/10.1002/2013JC009423, 2014. a
Smith, W. H. F. and Sandwell, D. T.: Bathymetric prediction from dense satellite altimetry and sparse shipboard bathymetry, J. Geophys. Res.Sol. Ea., 99, 21803–21824, https://doi.org/10.1029/94JB00988, 1994. a
Stammer, D., Ray, R. D., Andersen, O. B., Arbic, B. K., Bosch, W., Carrère, L., Cheng, Y., Chinn, D. S., Dushaw, B. D., Egbert, G. D., Erofeeva, S. Y., Fok, H. S., Green, J. A. M., Griffiths, S., King, M. A., Lapin, V., Lemoine, F. G., Luthcke, S. B., Lyard, F., Morison, J., Müller, M., Padman, L., Richman, J. G., Shriver, J. F., Shum, C. K., Taguchi, E., and Yi, Y.: Accuracy assessment of global barotropic ocean tide models, Rev. Geophys., 52, 243–282, https://doi.org/10.1002/2014RG000450, 2014. a
Tozer, B., Sandwell, D. T., Smith, W. H. F., Olson, C., Beale, J. R., and Wessel, P.: Global bathymetry and topography at 15 arc sec: SRTM15+, Earth Space Sci., 6, 1847–1864, https://doi.org/10.1029/2019EA000658, 2019. a
Van Sebille, E., Kehl, C., Lange, M., Delandmeter, P., and contributors: Parcels, Zenodo, https://doi.org/10.5281/zenodo.7035503, 2021. a
Zaron, E. D.: Mapping the nonstationary internal tide with satellite altimetry, J. Geophys. Res.Oceans, 122, 539–554, https://doi.org/10.1002/2016JC012487, 2017. a, b, c
Zaron, E. D.: Baroclinic tidal sea level from exactrepeat mission altimetry, J. Phys. Oceanogr., 49, 193–210, https://doi.org/10.1175/JPOD180127.1, 2019. a, b
Zaron, E. D.: Baroclinic tidal cusps from satellite altimetry, J. Phys. Oceanogr., 52, 3123–3137, https://doi.org/10.1175/JPOD210155.1, 2022. a, b, c
Zaron, E. D. and Egbert, G. D.: Timevariable refraction of the internal tide at the Hawaiian Ridge, J. Phys. Oceanogr., 44, 538–557, https://doi.org/10.1175/JPOD120238.1, 2014. a
Zaron, E. D. and Elipot, S.: An assessment of global ocean barotropic tide models using geodetic mission altimetry and surface drifters, J. Phys. Oceanogr., 51, 63–82, https://doi.org/10.1175/JPOD200089.1, 2021. a, b, c
Zhao, Z.: Mapping internal tides from satellite altimetry without blind directions, J. Geophys. Res.Oceans, 124, 8605–8625, https://doi.org/10.1029/2019JC015507, 2019. a
Zhao, Z., Alford, M. H., Girton, J. B., Rainville, L., and Simmons, H. L.: Global observations of openocean mode1 M2 internal tides, J. Phys. Oceanogr., 46, 1657–1684, https://doi.org/10.1175/JPOD150105.1, 2016. a
 Abstract
 Introduction
 Data
 Methods
 Results
 Conclusions
 Appendix A: Bias of the complex demodulate of a sample mean function
 Appendix B: Effects of the nontidal variability on the complex demodulates
 Appendix C: Estimating the autocovariance of the nontidal variability
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Data
 Methods
 Results
 Conclusions
 Appendix A: Bias of the complex demodulate of a sample mean function
 Appendix B: Effects of the nontidal variability on the complex demodulates
 Appendix C: Estimating the autocovariance of the nontidal variability
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References