Nonlinear Internal Tides in a Realistically Forced Global Ocean Simulation

1. Introduction

Internal tides are formed when barotropic tides propagate along steep slopes where barotropic energy is transferred to baroclinic (internal) tides at a rate of about 1 TW globally (Egbert and Ray 2000; Egbert and Ray 2001). The baroclinic motions induced by the tides play an important role in the vertical mixing of the deep ocean, affecting the thermohaline overturning circulation (Wunsch and Ferrari 2004; Melet et al. 2013; Waterhouse et al. 2014). While globally about half of the energy converted to internal tides is dissipated locally as high vertical modes (Vic et al. 2019), the remaining low modes radiate away. In both satellite altimetry (Zhao et al. 2016) and realistically forced global ocean models (Shriver et al. 2012; Muller et al. 2012; Buijsman et al. 2020) low-mode internal tides are observed to propagate across ocean basins. The dissipation of these low-mode internal tides is attributed to several mechanisms: 1) scattering at topography (Johnston and Merrifield 2003; Lamb 2004; Mathur et al. 2014), 2) scattering by mesoscale eddies (Kelly et al. 2016; Dunphy et al. 2017; Savage et al. 2020), 3) breaking at topography and shelves (Legg 2014; Klymak et al. 2016) and 4) wave-wave (triad) interactions (Olbers 1976; Varma and Mathur 2017; Eden et al. 2019; Baker and Sutherland 2020; Sutherland and Dhaliwal 2022).

In this paper we consider the decay of low-mode internal tides, generated at tall ridges or shelves, due to wave-wave interactions that drive a forward energy cascade to higher frequencies and wavenumbers. Linear internal tides manifest in the ocean as plane propagating wave modes that form horizontal wave beams. These modes can either self interact (Baker and Sutherland 2020; Sutherland and Dhaliwal 2022) or interact with other modes (Varma and Mathur 2017) in resonant triad interactions, generating modes at higher frequencies and wavenumbers. A commonly studied triad interaction, in which an internal tide interacts with an inertial mode of half the tidal frequency, referred to as Parametric Subharmonic Instability (PSI; e.g., Mackinnon et al. 2013; Ansong et al. 2018), is not considered in this study.

When the background rotation is weak and the stratification is surface intensified, a low-mode internal tide with significant energy may undergo a “superharmonic cascade" (Sutherland and Dhaliwal 2022), in which the mode disperses into a train of nonlinear internal waves (NLIWs) (Helfrich and Grimshaw 2008; Helfrich and Grimshaw 2008), also referred to as solitons or solitary waves. These waves have higher wavenumbers and frequencies than the mother wave they evolved from. This process, in which a long sinusoidal internal tide transforms into a short wavelength solitary wave of $\mathcal{O}\left(1\right)$ km, is also referred to as “nonlinear steepening". Analytical solutions for these NLIW result from the modified Korteweg-de Vries (KdV) equations (Ostrovsky and Stepanyants 1989), in which the advective term balances nonhydrostatic and Coriolis dispersion. Due to its higher frequency, the semidiurnal internal tide is more affected by nonlinear steepening than the diurnal internal tide. Equatorward of the diurnal turning latitude, Coriolis dispersion inhibits the steepening of the diurnal internal tide because its frequency is closer to the local Coriolis frequency than the frequency of the semidiurnal internal tide (Helfrich and Grimshaw 2008; Farmer et al. 2009).

NLIW may also result from other mechanisms than the steepening of the low-mode internal tide. These mechanisms include lee waves (Maxworthy 1979), tidal beams impinging on a pycnocline (New and Pingree 1990; New and Pingree 1992; Gerkema 2001; Grisouard et al. 2011; Mercier et al. 2012) and topographic scattering (Lamb 2004; Gerkema et al. 2006). Lee waves are generated over a ridge or bank when the tidal barotropic flow arrests upstream propagating internal wave modes. As the tidal flow relaxes the lee-wave may evolve into an internal solitary wave train. Solitary wave trains are also formed when a tidal beam, generated at a ridge or shelf, impinges on a moderately stratified pycnocline. Laboratory experiments confirm the generation of higher-harmonics at multiples of the primary frequency and solitary waves at even higher frequencies by tidal beams impinging on a pycnocline (Mercier et al. 2012). In topographic scattering, the internal tide interacts with supercritical topography to generate higher-harmonics.

Low-mode nonlinear internal waves with large enough amplitudes in the interior are easily detectable with Synthetic Aperture Radar (SAR) and optical sunglint imagery from satellite (Alpers 2018; Jackson et al. 2012; Jackson et al. 2013) because the waves’ surface currents affect ocean surface roughness. Using optical sunglint imagery, Jackson (2007) compiled a global map of NLIWs that shows they are ubiquitous globally in coastal shelf areas and in the (sub)tropics in the open ocean, e.g., in the South China Sea (Simmons et al. 2011; Alford et al. 2015), near the Mascarene Ridge (Konyaev et al. 1995; da Silva et al. 2015), in the Bay of Bengal and Andaman Sea (Osborne and Burch 1980; Jensen et al. 2020; Magalhaes et al. 2020), and offshore of the Amazon Shelf (Ivanov et al. 1990; Magalhaes et al. 2016). The higher harmonic frequencies associated with these nonlinear internal waves have also been observed with moorings, e.g., in the Bay of Biscay (van Aken et al. 2007) and globally with surface drifters (Yu et al. 2019).

Regardless of the generation mechanism, the presence of NLIW represents a potential energy sink for low-mode internal tides and a pathway to dissipation. Strong vertical mixing associated with large amplitude low-mode NLIW has been observed in areas such as the Oregon shelf (Moum et al. 2003), the Monterey Bay shelf (Carter et al. 2005), the New Jersey shelf (Shroyer et al. 2010), and the South China Sea (Alford et al. 2015; Zhang et al. 2023).

High-frequency nonlinear (solitary) internal waves have generally been studied in high-resolution process and regional nonhydrostatic model simulations with horizontal resolutions of $\mathcal{O}\left(100\right)$ m (e.g., Simmons et al. 2011; da Silva et al. 2015; Magalhaes et al. 2020; Buijsman et al. 2010). It is not yet clear how well they are simulated in realistically forced hydrostatic global ocean model simulations with horizontal resolutions of $\mathcal{O}\left(1\right)$ km. The appearance of solitary-like waves in hydrostatic models (Jensen et al. 2020) has been attributed to numerical dispersion that mimics the nonhydrostatic dispersion (Vitousek and Fringer 2011). Global ocean models do simulate internal waves at supertidal frequencies (e.g., Müller et al. 2015; Savage et al. 2017; Luecke et al. 2020; Arbic et al. 2022). However, these latter analyses have not specifically focused on tidally generated low-mode nonlinear internal waves, which are the target of this paper. While the depth-integrated energy budget of the semidiurnal internal tide has been well investigated in realistically forced global ocean simulations (Buijsman et al. 2020; Müller 2013; Buijsman et al. 2016; Li and von Storch 2020), such energy budget is lacking for low-mode nonlinear internal waves that have evolved from the primary frequency internal tides.

The overarching goal of this paper is to understand the decay of low-mode internal tides due to wave-wave interactions that drive a forward energy cascade to supertidal frequencies in a global Hybrid Coordinate Ocean Model (HYCOM) simulation with a 4-km horizontal resolution (Raja et al. 2022). Specific goals of this study are to 1) quantify the kinetic energy in the supertidal band, 2) quantify the energy transfer rate from the tidal to supertidal frequencies, 3) identify the mechanisms responsible for these energy transfers, and 4) compare simulated hotspots of supertidal internal wave energy with observations of NLIW.

The paper layout is as follows: the ocean model and methods used for this study are briefly described in section 2. Global internal tide energetics are presented in section 3, and discussed in more detail for the main internal tide beam radiating from the Amazon Shelf. In section 4 we correlate the occurrence of simulated supertidal energy flux with satellite imagery of solitary internal waves globally. We then estimate the nonlinear energy transfer rate from primary to higher-harmonic frequencies based on the coarse-grained kinetic energy transfer, and attribute the spatial patterns of the energy transfers to constructive interference between modes. Finally, we summarize our findings and make suggestions for future work in section 5.

2. Model and Methods

2.1. Model

HYCOM is a primitive equation ocean general circulation model (Bleck 2002) and the operational global ocean forecast model of the United States Navy (Metzger et al. 2014). The model uses a hybrid vertical coordinate that is isopycnal in the open ocean and transitions to terrain-following in shallow water, with z-coordinates used to resolve the surface mixed layer. The simulation presented in this paper is run in a forward (non-assimilative) mode on a tripolar grid at 1/25^∘ (4 km) nominal horizontal resolution with 41 nonuniform layers in the vertical. In the mixed layer, vertical grid points are located at fixed depths most of which are 8 m apart, while in the deep ocean they move vertically with their assigned isopycnals. This grid resolution allows HYCOM to resolve at least four semidiurnal internal waves modes in deep water (Buijsman et al. 2020).

The model is forced with realistic atmospheric fields from the Navy Global Environmental Model (Hogan et al. 2014) and tidal forcing for the M${}_2$, S${}_2$, K${}_1$, O${}_1$, and N${}_2$ tidal constituents. To account for numerical errors in the tidal solution due to imperfect topography and damping terms, an Augmented State Ensemble Kalman Filter is applied to optimize the spatially varying Self Attraction and Loading term used in the simulations (Ngodock et al. 2016). The sea surface height root mean squared error with TPXO8-atlas, an altimetry constrained model, is 2.6 cm. The surface tidal error is minimized by tuning the strength of the linear wave drag (Jayne and St. Laurent 2001), which is a parameterization of the surface-tide energy conversion to the unresolved baroclinic wave modes.

Model output from global HYCOM simulations has been validated extensively against moorings (Luecke et al. 2020), altimetry (Shriver et al. 2012; Buijsman et al. 2020; Buijsman et al. 2016), drifters (Arbic et al. 2022) and other models (Savage et al. 2017; Arbic et al. 2022). Subtidal kinetic energy is lower in HYCOM compared to surface drifters, especially near the equator. Near inertial kinetic energy shows good agreement with drifters, where the high (3-hr) wind forcing frequency plays an important role. Compared to altimetry, HYCOM is more energetic in the semidiurnal band, about 23% in hot-spot regions and 9% in global-mean variance (Buijsman et al. 2020). Diurnal and semidiurnal frequency band integrated steric sea surface height (SSSH) variances in HYCOM match reasonably well with tide gauges, although the model is clearly deficient at supertidal frequencies (Savage et al. 2017). Comparison with higher horizontal grid resolution models suggests that the variance in the supertidal band in HYCOM has not yet converged, and is significantly lower than observations.

In this paper we demonstrate that our HYCOM simulation can resolve the initial stages of the superharmonic cascade associated with the low-mode internal tide. We demonstrate in Appendix A) that the simulation does not resolve NLIW modes with wavelengths $<8$ km that can result from this cascade. In Appendix A), we compare the generation of NLIW in the hydrostatic HYCOM simulation with a high-resolution nonhydrostatic MITgcm simulation at the Mascarene Ridge (da Silva et al. 2015). We find that HYCOM underestimates the supertidal energy by about 50% at this location.

2.2. Filtering

We compute the energy terms for the filtered 3D fields from our global HYCOM simulation. We use hourly output from global HYCOM for a 30-day period, starting May 20 and ending June 20, 2019. We use a 4^th order Butterworth band-pass filter with 2-30 hr (green), 9-15 hr (blue) and 2-9 hr (red) cutoff periods to evaluate total internal tide energy (D1+D2+HH), semidiurnal (D2) and higher-harmonic (HH) frequency bands. We remove the first and last days to avoid ringing effects such that fields are time-averaged over a 28-day period. Figure 1 shows the horizontal energy frequency spectra of baroclinic surface velocities offshore of the Amazon Shelf, where tidal peaks (D1, D2) and higher-harmonic (HH) peaks can be observed. The HH (supertidal) peaks occur at multiples of the primary tidal frequencies (e.g., at 3, 4, 5, 6, etc., cycles per day (cpd)).

Because internal tide energy is found at narrow-band peaks, energy estimates are not very sensitive to the choice of the cutoff frequencies, although some small errors may be observed, for example around 3 cpd (Figure 1). The 30-day simulation period with hourly output frequency restricts our analysis to motions with a minimum fundamental frequency of 1/30 cpd and a maximum Nyquist frequency of 12 cpd, which captures our lowest frequency of interest (D1) and at least the first few supertidal peaks (D3-D6). The simulation period also allows for two spring-neap cycles, which is important for the generation cycle of nonlinear internal tides.

2.3. Internal Wave Energetics

The depth-integrated internal tide energy balance provides insight into energy sources and sinks in the ocean, and can be written as (Buijsman et al. 2016; Kang and Fringer 2012)

$$\frac{\partial \overline{E}}{\partial t}+{\nabla }_H·\overline{\mathbf{F}}=\overline{C}-\overline{R},$$

(1)

where the overbars indicate depth-integration, $\frac{\partial \overline{E}}{\partial t}$ is the tendency term, ${\nabla }_H·\overline{\mathbf{F}}$ is the horizontal divergence (${\nabla }_H=\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y}\right)$) of the depth integrated baroclinic energy flux $\overline{\mathbf{F}}$, $\overline{C}$ is the barotropic to baroclinic energy conversion, and $\overline{R}$ is the residual. In the absence of wave-mean flow and wave-wave energy exchange energy terms (Kelly 2016; Pan et al. 2021; Barkan et al. 2021), the residual term is interpreted as baroclinic energy dissipation.

The total baroclinic wave energy ($\overline{E}$) in equation (1) is the sum of the horizontal wave kinetic energy ($\overline{KE}$) and available potential energy ($\overline{APE}$). For a water column height h the horizontal wave kinetic energy is computed as

$$\overline{KE}=\frac{1}{2}{\rho }_0{\int }_{-h}^0{\left|{\mathbf{u}}_H^{\prime }\right|}^{2}dz,$$

(2)

where ${\rho }_0$ is the reference density and ${\mathbf{u}}_H^{\prime }=\left(u^{\prime },v^{\prime }\right)$ the baroclinic horizontal velocity vector along the x and y coordinates. The linear form of available potential energy is computed as (Kang and Fringer 2010)

$$\overline{APE}=\frac{1}{2}{\rho }_0{\int }_{-h}^{0}N{\left(z\right)}^2\xi {\left(z\right)}^{2}dz,$$

(3)

where $\xi$ is the vertical isopycnal displacement and $N\left(z\right)$ is the buoyancy frequency. Although maps of time-mean $\overline{APE}$ feature similar patterns as time-mean $\overline{KE}$, we only present maps of time-mean $\overline{KE}$ because maps of time-mean $\overline{APE}$ are noisy.

The depth-integrated hydrostatic baroclinic energy flux $\overline{F}$ is a measure of the horizontal internal wave energy transport. It can be written as

$$\overline{\mathbf{F}}=\overline{{\mathbf{u}}_H^{\prime }E}+\overline{{\mathbf{u}}_H^{\prime }{p}^{\prime }},$$

(4)

where ${\mathbf{F}}_{\mathbf{a}}=\overline{{\mathbf{u}}_H^{\prime }E}$ represents nonlinear self-advection, ${\mathbf{F}}_{\mathbf{p}}=\overline{{\mathbf{u}}_H^{\prime }{p}^{\prime }}$ represents hydrostatic pressure work, and $p^{\prime }$ the perturbation pressure. The perturbation velocities and pressure have depth-integrals that are equal to zero, and are computed as described in Nash et al. (2005). Although advective energy fluxes can be as large as pressure fluxes in strongly nonlinear waves in high resolution Large Eddy Simulations (Venayagamoorthy and Fringer 2005), we will show that the contribution from energy advection is small in our HYCOM simulation.

The barotropic to baroclinic conversion is a source term in equation (1), and in the absence of nonhydrostatic effects can be computed as

$$\overline{C}=-p\left(z=-h\right)\overline{{\mathbf{u}}_H}·{\nabla }_{H}h,$$

(5)

where ${\nabla }_{H}h$ is the horizontal gradient of the ocean seafloor depth.

Finally, after time averaging eq. (1), the internal wave energy budget for both the tidal and supertidal bands is

$$\langle \overline{C}\rangle ={\nabla }_H·\langle \overline{\mathbf{F}}\rangle +\langle \overline{R}\rangle .$$

(6)

3. Results

In this section we discuss the geographical distribution of internal tide energy in the abyssal ocean, neglecting areas shallower than 250 m, and then explore energy sources and sinks. We limit our analysis to low-latitudes, between ${25}^{\circ }$S and ${25}^{\circ }$N, to provide more detail where higher-harmonic energy is important and to avoid aliasing of high frequency motions due to mesoscale eddies and near-inertial waves.

3.1. Global internal tide energetics

The time mean of the depth-integrated baroclinic kinetic energy (eq. 2) is shown in Figure 2 for the (a) diurnal, (b) semidiurnal and (c) higher-harmonic (supertidal) frequency bands. To highlight regions where supertidal waves are important, the ratio of higher-harmonic to total tidal energy (D1+D2+HH) is plotted in Figure 2d. The kinetic energy in the diurnal band contains energy from both diurnal internal tides and equatorward propagating near-inertial waves generated at the critical latitudes ($\sim \pm {30}^{\circ }$). Diurnal internal tides are most prominent at the Luzon Strait and Indonesian Archipelago (rectangles 2 and 3 in Figure 2d), while semidiurnal energy is more ubiquitous in the global ocean. In the HYCOM simulation analyzed here, supertidal kinetic energy is 1-2 orders of magnitude lower than tidal energy almost everywhere, except along strong internal tide beams at low-latitudes, most notably at the Bay of Bengal and Amazon Shelf (rectangles 1 and 6 in Figure 2d), where supertidal energy makes up 30-50% of the total energy.

The spatial distribution of internal tide energy yields some insights into their generation. Kinetic energy at diurnal and semidiurnal frequencies is highest near shallow and steep ridges, where energy is transferred from barotropic to baroclinic modes, and then decays rapidly with distance from the source. Conversely, supertidal kinetic energy can peak away from the generation sites of the internal tides, suggesting that the source of supertidal kinetic energy is the propagating internal tide. For example, equatorward of Hawaii (rectangle 4 in Figure 2d) and equatorward of the French Polynesian Islands (FPI; rectangle 5 in Figure 2d), supertidal kinetic energy emerges several hundred kilometers away from the generation sites along the path of the primary internal tide. The supertidal $KE$ fraction is much smaller for beams radiation poleward from Hawaii and the FPI (Figure 2d). The asymmetry between equatorward and poleward radiating supertidal beams from Hawaii and FPI suggests that the Coriolis and buoyancy frequencies play an important role in the generation of superharmonics (Sutherland and Dhaliwal 2022).

The magnitudes of the depth-integrated and time-averaged horizontal tidal and HH pressure and advective fluxes are presented in Figure 3. The energy flux is dominated by pressure work and resembles the distribution of $KE$ in both frequency bands (Figure 3a,b). Higher-harmonic fluxes in Figure 3b feature beams of $\mathcal{O}\left(1\right)$ kWm${}^{-1}$ appearing some distance away from the generation sites and remaining steady for several hundred kilometers. The tidal band advective flux in Figure 3c is $\mathcal{O}\left(10\right)$ Wm${}^{-1}$, which is negligible compared to the tidal pressure flux. In contrast, the advective flux in the supertidal band in Figure 3d is larger and correlates well with the supertidal pressure flux and kinetic energy beams in the equatorial region. In these regions, internal tides have transferred a significant portion of their energy to waves at higher wavenumbers and frequencies, which contribute more to the nonlinear advective fluxes.

To better understand the latitudinal dependence, we zonally integrate various energy terms for $1^{\circ }$ bins from ${25}^{\circ }$S to ${25}^{\circ }$N (Figure 4). On average, supertidal $KE$ accounts for 5% of the total internal tide energy at low-latitudes. It gradually increases from the mid to low latitudes reaching the largest value of about 10% at $4^{\circ }$N (Figure 4a), likely due to contributions from the Bay of Bengal and the Amazon shelf (Figure 2d). The peak in diurnal $KE$ near ${20}^{\circ }$N is attributed to the internal tide generation in Luzon Strait (rectangle 2 in Figure 2d). Diurnal $KE$ increases poleward towards the diurnal turning latitudes due to adulteration by near-inertial waves. Globally integrating between these latitudes, we estimate that $K{E}_{tot}=113.23$ PJ, where $K{E}_{D1}=31.67$ PJ, $K{E}_{D2}=76.14$ PJ and $K{E}_{HH}=5.42$ PJ. Figure 4b shows that most of the energy converted from the surface to the baroclinic tide is input at semidiurnal frequencies (440 GW) as opposed to diurnal frequencies (60 GW), with a negligible amount converted to supertidal frequencies. The residuals of the tidal and supertidal energy budgets (eq. 6) in Figure 4c are positive for the tidal band, mostly due to dissipation, and negative for the supertidal band, suggesting a source term is missing. To better understand how internal tide energy at the primary frequencies is transferred to supertidal frequencies, we analyze the energetics along the main internal tide beam generated at the Amazon Shelf in the next section.

3.2. Nonlinear internal tides at the Amazon shelf

The Amazon Shelf region (rectangle 6 in Figure 2d) features an energetic internal tide beam, in which up to half of the internal tide energy is in the supertidal band. In this beam, nonlinear internal solitary waves appear about 500 km from the coastline, coinciding with the rapid decay of the primary tide energy fluxes (Magalhaes et al. 2016). The direction of the main wave beam is modulated by the boundary currents along the shelf, namely the North Brazil Current (NBC) and North Equatorial Counter Current (NECC) (Magalhaes et al. 2016; Tchilibou et al. 2022). In this section we analyze the energy budget terms along the main internal tide beam generated at the Amazon Shelf to indicate potential sources of supertidal energy.

The main internal tide beam is generated in at least 3 hot spots along the steep Amazon Shelf (Figure 5a), in agreement with previous studies (Magalhaes et al. 2016; Tchilibou et al. 2022). Energy flux at tidal frequencies is high close to the generation sites and along the propagation path of the internal tide and reaches values of over 14 kWm${}^{-1}$. The supertidal energy flux in Figure 5b is much lower at the generation site but increases significantly between 200 and 500 km from the shelf break and then remains relatively steady at ∼2 kWm${}^{-1}$ until the Mid-Atlantic Ridge (top right contours in Figure 5b).

To understand the spatial distribution of energy sources and sinks, we analyze the baroclinic energy budget (eq. 6) averaged over rectangular bins of 0.5^o by 1.5^o (55.5 km by 166.5 km) along the propagation direction of the main internal tide. The bins are shown in Figure 6e superimposed over supertidal flux divergence, where sources (sinks) are represented by the red (blue) colors. The bin-averaged energy flux in Figure 6a is mainly in the tidal frequency band. It reaches a maximum value of ∼9 kWm${}^{-1}$ at the shelf slope, which can be attributed to the barotropic to baroclinic energy conversion in the tidal band in Figure 6b. On the contrary, the supertidal energy flux is initially low (0 Wm${}^{-1}$) and increases offshore until it peaks at ∼2 kWm${}^{-1}$ about 500 km from the source. This increase does not coincide with the supertidal conversion, which is negligible in Figure 6b.

The decay of the tidal energy flux and the growth of supertidal flux is more easily distinguished with the bin-averaged horizontal flux divergence in Figure 6c. Along the propagation path, the energy flux divergence at the primary and higher-harmonic frequencies appears to be anti-correlated, suggesting that energy is transferred from tidal to supertidal frequencies (and vice versa) at intervals spaced roughly 150 km apart, e.g., near 220 and 370 km (near 130 and 300 km).

The residual $\langle \overline{R}\rangle$ (eq. 1) of the energy balance for the tidal band is generally positive in Figure 6d, in agreement with the overall decay of the flux in Figure 6a. The tidal residual is negative near 300 km, suggesting that terms may be missing, e.g., due to wave-mean flow advection (Kelly et al. 2016) or energy exchanges between different frequency bands (Kelly et al. 2016; Barkan et al. 2021). Although wave-mean flow interactions are not the focus of this study, they may be a source or sink of energy for internal tides in the presence of mesoscale currents, such as the NECC and NBC. As in Figure 4, the supertidal residual is mostly negative, implying that an important source term is not accounted for. The mostly negative supertidal residual term and the spatial correlation between the tidal and supertidal flux divergence and residual terms suggest a net energy flux from the tidal to the supertidal band. This energy exchange will be investigated in more detail in the discussion section.

3.3. Solitary-like waves in HYCOM

In HYCOM, the internal tide at the Amazon Shelf steepens into solitary-like waves near 400 km (Figure 7), coinciding with the second patch of positive supertidal flux divergence (Figure 6e). At this location, the steric sea surface height (SSSH) in Figure 7a and the thermocline displacement in Figure 7c feature the largest amplitudes and steepest slopes. Simultaneously, the vertical velocities at depth in Figure 7b and the supertidal kinetic energy at the surface in Figure 7b are enhanced. The separation between the peaks of SSSH is about 110 km, which is the mode 1 wavelength at this location, implying that it is the mode 1 wave that steepens.

To validate the accuracy of the simulated solitary-like waves near the Amazon Shelf, we compare SSSH fields from HYCOM with crests extracted from SAR imagery (magenta lines) in Figure 8 during the same time-period (May-June). The internal wave crests in HYCOM steepen in regions where the solitary NLIW are observed. The satellite derived solitary waves emerge around $4^{\circ }$N and $44.5^{\circ }$W, roughly 400 km from the shelf break (see also Figure 7). These figures show that in HYCOM, large amplitude internal tides steepen into solitary-like waves as they propagate offshore in regions where supertidal energy is relatively large, in good agreement with satellite observations.

4. Discussion

In this section we further discuss the correlation between supertidal energy in the global HYCOM simulation and the observations of solitary waves in the global ocean, we estimate the nonlinear energy transfers to supertidal frequencies using a coarse-graining approach, and we attribute the spatio-temporal modulation of the nonlinear energy transfers to constructive interference between modes.

4.1. Solitary NLIW in global HYCOM

For the Amazon region we have shown that a strong supertidal energy flux of $\mathcal{O}\left(1\right)$ kWm${}^{-1}$ coincides with NLIW that are both simulated in HYCOM and observed in SAR images. To confirm this correlation globally, we compare the supertidal energy flux with the locations of NLIW observed by the 250 m Moderate-Resolution Imaging Spectrometer (MODIS) satellite over a four-year period (Jackson 2007) in Figure 9. Most regions in the abyssal ocean with supertidal energy fluxes $>1$ kWm^-1 (purple shading in Figure 9), such as the Bay of Bengal, Andaman Sea, South China Sea, Indonesian Archipelago, Amazon Shelf, and Mascarene Ridge, agree well with the locations of the observed solitary waves. However, solitary waves are also detected in regions with relatively low supertidal energy flux, such as over large areas in the eastern Atlantic and Pacific Oceans (rectangles 8 and 9 in Figure 9b) and on the continental shelves. A recent study (da Silva et al., in review) demonstrates that the solitary waves in the eastern Pacific are caused by non-tidal sources, possibly explaining the discrepancies in the eastern Pacific and Atlantic Oceans. The discrepancy between HYCOM and MODIS in rectangles 8 and 9 in Figure 9b may also be explained by weaker surface tides and/or unresolved bathymetric features in the HYCOM simulation, which may cause weaker primary and higher-harmonic internal wave flux beams. NLIW on shelves generally feature smaller wavelengths and energy fluxes ($<1$ kWm^-1) than in the open ocean. Hence, these waves are not well resolved in our HYCOM simulations, as opposed to the satellite images, which horizontal resolution is about eight times finer.

4.2. Nonlinear energy cross-scale transfer

The negative residuals of supertidal energy at the Amazon Shelf in Figure 4c and Figure 6d-e, point to a missing source term in the higher-harmonic energy budget. Therefore, we diagnose the nonlinear energy wave transfers as a function of time and space in HYCOM with a coarse-graining approach (Eyink and Aluie 2009). The coarse-graining method has been used in physical oceanography to map the energy cascade in the North Atlantic (Aluie et al. 2018) and to quantify the energy transfer between mesoscale eddies and internal waves (Barkan et al. 2021). To our knowledge, the coarse-graining method has not yet been used to study the energy pathways of internal tides in the global ocean. Our use of the coarse-graining method here is therefore novel.

Coarse-graining can be used to quantify energy transfers across an arbitrary length or time scale $\tau$. The coarse-grained kinetic energy transfer (Wm^-3) can be written in index notation as:

$${\Pi }_{\tau }(x,t)=-{\rho }_0(\widetilde{u_i{u}_j}-\widetilde{u_i}\widetilde{u_j})\frac{\partial \widetilde{u_i}}{\partial x_j}$$

(7)

where $i=1,2,3$ and $j=1,2$ are indices indicating directions $(x,y,z)$ and $\widetilde{(·)}$ denotes the low-pass filter of quantity $(·)$ with cutoff period $\tau$. We use $\tau =9$ hr (2.67 cpd) as our cutoff period (frequency), such that ${\Pi }_{\tau }$ represents the energy transfer rate between (sub)tidal and supertidal frequencies. However, we conclude from our results that the majority of this energy transfer likely comes from the semidiurnal internal tide.

In the absence of topographic conversion we now can approximate the HH energy balance as

$$\langle \overline{{\Pi }_{\tau }}\rangle ={\nabla }_h·\langle \overline{{\mathbf{F}}_{HH}}\rangle +\langle \overline{R}\rangle$$

(8)

where the residual $\langle \overline{R}\rangle$ includes terms such as the dissipation and wave-mean flow interaction terms.

We compare the nonlinear cross-scale energy transfer to the divergence of the supertidal flux in three areas with strong HH energy fluxes in Figure 10a-f: the Bay of Bengal, the Amazon Shelf and Mascarene Ridge. To reduce noise due to aliasing of high frequency motions, we average the supertidal flux divergence over a nine grid point stencil in Figure 10d, which yields a smoother appearance than the flux divergence in Figure 6e. The positive regions of cross-scale energy transfer indicate a forward cascade from tidal to supertidal frequencies, and they coincide with estimates of supertidal energy flux divergence. The energy transfer at all 3 regions in Figure 10a-c is mostly forward and form a banding pattern along the propagation path of the internal tides at regularly spaced intervals separated by a distance longer than the mode 1 internal wavelength. Note that while in some areas internal tides transfer significant energy to higher-harmonics close to the generation sites (e.g., on the east side of the Mascarene Ridge), in other areas, such as the Amazon Shelf, most of their energy is transferred away from topography.

One goal of this paper is to estimate the energy scale transfer rate from tidal to supertidal frequencies. The time-mean and area integrated coarse-grained $KE$ between $\pm {25}^{\circ }$ is 44.5 GW forward and 7.5 GW inverse transfer from (sub)tidal to higher-harmonic frequencies. This implies that about 7% of the energy converted from barotropic to baroclinic tides ($\sim 500$ GW) between $\pm {25}^{\circ }$ cascades to higher-harmonics in the deep ocean, assuming that most of these energy transfers come from the primary tides.

4.3. Spatial modulation of nonlinear energy transfer

In this section, we investigate the mechanisms that cause the regularly spaced banding patterns in the supertidal flux divergence and the nonlinear cross-scale energy transfers. A similar banding pattern is also visible in tidally band-passed surface kinetic energy (Figure 10g-i), while it is not present in the depth-integrated tidal $KE$ (Figure 10j-l). Although we compute $KE$ for the tidal (D1+D2) band in Figure 10g-l, the banding patterns of time-mean $KE$ are mostly due to the semidiurnal tides (not shown), which are the most dominant at these internal tide hotspots (Figure 2a,b). In this banding pattern, the patches of large tidal surface $KE$ also coincide with large nonlinear energy transfers. This is analogous to internal-wave dissipation scaling with energy density (e.g., Pollmann et al. 2017) because shear-driven turbulent internal wave dissipation is the end result of a forward energy cascade due to wave-wave interactions. The spatial patchiness of $KE(x,y,z)$ is attributed to the horizontal and vertical shear in the tidal baroclinic velocities, which drive the nonlinear energy transfers in eq. (7). The difference in banding patterns between surface $KE$ and depth-integrated $KE$ suggests that the banding patterns are caused by constructive and destructive interference between semidiurnal mode 1 and mode 2 waves, which are generally the most energetic waves generated at tall ridges and shelves in 4-km HYCOM simulations see Figure 8f-g of (Buijsman et al. 2020).

We decompose the baroclinic fields into vertical modes to demonstrate that the banding pattern in Figure 10 can be attributed to the interference patterns between semidiurnal mode 1 and mode 2 waves. We adopt the modal amplitude time series for baroclinic velocity ${\widehat{\mathbf{u}}}_n(x,y,t)$ and pressure ${\widehat{p}}_n(x,y,t)$, the normalized horizontal velocity eigenfunctions ${\mathcal{U}}_n(x,y,z)$, and eigenspeeds $c_n$ as computed by Raja et al. (2022) for the same HYCOM simulation used in this study.

The kinetic energy density for any frequency band at a given depth and time for the sum of the mode 1 and 2 velocities reads (omitting z and t on the left-hand side and x and y everywhere)

$$\begin{array}{cc}\hfill KE(z,t)& =\frac{1}{2}{\rho }_0\left[{({\widehat{u}}_1{\mathcal{U}}_1+{\widehat{u}}_2{\mathcal{U}}_2)}^2+{({\widehat{v}}_1{\mathcal{U}}_1+{\widehat{v}}_2{\mathcal{U}}_2)}^2\right]\hfill \\ \hfill & =\frac{1}{2}{\rho }_0\left[{\widehat{u}}_1^2{\mathcal{U}}_1^2+{\widehat{u}}_2^2{\mathcal{U}}_2^2+{\widehat{v}}_1^2{\mathcal{U}}_1^2+{\widehat{v}}_2^2{\mathcal{U}}_2^2+2{\widehat{u}}_1{\widehat{u}}_2{\mathcal{U}}_1{\mathcal{U}}_2+2{\widehat{v}}_1{\widehat{v}}_2{\mathcal{U}}_1{\mathcal{U}}_2\right],\hfill \end{array}$$

(9)

where the last 2 terms are the cross terms, which disappear when depth integrating because of the orthogonality condition. If the mode 1 and 2 velocities are in phase at the surface the cross-terms are positive, enhancing $KE$. At the same time, the mode 1 and 2 bottom velocities are out of phase near the bottom causing negative cross terms at depth, decreasing $KE$. In the resulting horizontal $KE$ interference pattern, the amplified kinetic energy will always be largest near the surface because the ${\mathcal{U}}_n$ eigenfunctions have the largest values near the surface. In the following we will show that the phase relation between the semidiurnal mode 1 and 2 waves (and thus the sign of the cross terms) cause the banding patterns in surface energy and energy transfers.

The time-mean semidiurnal surface $KE$ computed for the mode 1 velocities, the mode 2 velocities, and the sum of the mode 1 and mode 2 velocities (eq. 9) at the Bay of Bengal, Amazon Shelf, and Mascarene Ridge is shown in Figure 11. Both mode 1 and mode 2 surface $KE$ is high along the path of the internal tides in Figure 11a-f, but the beams lack the banding patters. However, the constructive interference between mode 1 and 2 results in a banding pattern, which patches occur at the same locations as the undecomposed surface tidal $KE$ in Figure 10g-i (and the the nonlinear cross-scale energy transfer and supertidal flux divergence in Figure 10a-f).

As the semidiurnal mode 1 and 2 waves propagate away from their generation sites, the faster mode 1 waves overtake the slower mode 2 waves. This causes alternating patterns of enhanced and reduced surface $KE$ due to positive and negative cross-term values in eq. (9). If we assume that the semidiurnal mode 1 and mode 2 waves are generated at the same time every tidal cycle with period T, the first patch occurs at the generation site and the second patch at an separation distance $\Delta S=\frac{c_1{c}_{2}T}{c_1-c_2}$, where $c_1$ and $c_2$ are the M${}_2$ mode 1 and 2 phase speeds. It follows from Table 1 that the predicted separation distance agrees well with the “observed" separation distance $\Delta S_o$ for the Bay of Bengal, the Amazon Shelf and the Mascarene Ridge.

To confirm this interference pattern is surface intensified, we present the time-mean semidiurnal $KE(x,z)$ for the modes and the time-mean nonlinear energy transfer along the main internal tide beam at the Amazon Shelf in Figure 12. The patches of surface intensified nonlinear energy transfer near 50, 200, and 380 km in Figure 12d occur above the thermocline and coincide with enhanced surface tidal kinetic energy, computed from the superposed mode 1 and 2 velocities in Figure 12c. In contrast, the mode 1 and mode 2 $KE$ in Figure 12a and Figure 12b do not show these banding patterns.

4.4. Temporal modulation of nonlinear energy transfer

Finally, we correlate the spatial-temporal variability of the 30-hr low-pass filtered nonlinear energy transfers to the 30-hr low-pass filtered semidiurnal modal surface $KE$, surface density, and supertidal steric sea surface height variance in Figure 13 along the main internal tide beam at the Amazon Shelf. The regions of enhanced nonlinear forward energy transfer in Figure 13a near 50, 200, and 380 km coincide with enhanced surface $KE$ for modes 1+2 in Figure 13b with a clear modulation by the spring-neap cycle.

Moreover, the regions of nonlinear energy transfer are also modulated by the slowly varying background flow and stratification. For example the offshore migration of the density front associated with the Equatorial Currents at 300-350 km in Figure 13c coincides with an increase of the separation distance between the patches after June 10^th. At the same time, the nonlinear energy transfers and the generation of nonlinear internal waves in Figure 13d are reduced. Thus, the time varying background flow impacts the generation and propagation of tidal modes, which in turn affect the nonlinear energy transfers. The regular spaced banding patterns in SSSH are attributed to constructive supertidal mode 1 and 2 interference, which is dominated by the M${}_4$ frequency (Figure 1).

4.5. Synthesis

The regularly spaced banding patterns, e.g., in semidiurnal surface $KE$, are not unique to our simulation. Similar patterns in semidiurnal SSSH variance at the same locations and with the same separation distances have been detected in satellite altimetry and regional model simulations of the Amazon Shelf region as shown in Figure 2c,d of Tchilibou et al. (2022). As in our simulations, the patches shown in Tchilibou et al. (2022) occur near $2.5^{\circ }$N and $4^{\circ }$N along the main beam. We note that the banding patterns in semidiurnal SSSH variance in our HYCOM simulation occur at the same locations as the patterns in $KE$ (result not shown).

These regularly spaced patches in surface intensified semidiurnal energy and SSSH variance are the locations of nonlinear energy transfers to higher harmonics. In our simulations these patches are due to constructive interference between semidiurnal mode 1 and 2 waves, which are generally the dominant modes radiating from tall ridges (Buijsman et al. 2020). Hence, maps of mode 1+2 interference could be used as a proxy to identify locations of such nonlinear energy transfers, which are cumbersome to compute. Similarly, maps of constructive interference between semidiurnal and diurnal modes and/or higher modes may also provide insight in the locations of these energy transfers.

These energy transfers enhance the supertidal energy density in the equatorial region (Figure 2d), where $f/N$ is small. These findings are in agreement with the theories of Baker and Sutherland (2020) and Sutherland and Dhaliwal (2022) that predict the evolution of supertidal wave modes from mode-1 internal tides as they propagate equatorward where f is small and N is large. According to Sutherland and Dhaliwal (2022), the first step in this superharmonic energy cascade is a triad-interaction of the primary frequency internal tide with itself, generating a mode-1 wave with a frequency that is moderately different from twice the internal tide frequency. The coincidence between the energy transfers and the constructive interference between semidiurnal mode 1 and 2 internal tides at the strongest generation sites (Figure 10 and Figure 11) suggest that a triad wave interaction between semidiurnal mode 1 and 2 may be important. This interaction generates a mode 1 at nearly twice the primary frequency (Varma and Mathur 2017). It is shown in (), that of the supertidal modes, mode 1 and mode 2 are only resolved due to the horizontal grid resolution in the 4-km HYCOM simulations and that generally the quarterdiurnal (e.g., M${}_4$) higher harmonic frequency is most dominant. However, it is beyond the scope of this paper to determine what kind of superharmonic triad interaction is important here. We plan to address that in a future study.

Internal tide dissipation due to wave-wave interactions is generally attributed to PSI (e.g., Pollmann et al. 2017; Eden and Olbers 2014; de Lavergne et al. 2019). Ansong et al. (2018) diagnosed PSI in 4-km and 8-km HYCOM simulations forced with either semidiurnal or diurnal tides. For the 4-km HYCOM simulation with semidiurnal tides, Ansong et al. (2018) estimated an energy transfer of 40 GW due to PSI out of the semidiurnal internal tide. The energy transfer due to PSI in an 8-km simulation with diurnal tides was only 10 GW. The 40 GW is comparable to the 45 GW of energy transfer from tidal to supertidal frequencies that we find in our study equatorward of $\pm {25}^{\circ }$. This suggests that energy transfers out of low-mode internal tides due to nonlinear steepening cannot be ignored. We note that both subharmonic and superharmonic energy transfers increase significantly with finer horizontal grid resolutions (Egbert and Ray 2000; Ansong et al. 2018; Nelson et al. 2020), pointing to the importance of small scale internal tides in resolving wave-wave interactions.

Finally, we note that the decay of the semidiurnal phase-locked (also referred to as stationary or coherent) internal tide in the equatorial Pacific as observed in satellite altimetry (Zhao et al. 2016; Zaron 2017; Buijsman et al. 2017) may not be only attributed to the time variability of the equatorial currents and tropical instability waves, which scatters energy to the non-phase locked internal tide (Buijsman et al. 2017), but also to the energy transfer to higher harmonics (rectangle 4 and 5 in Figure 2d. For example, for the equatorward beam from the French Polynesian Islands (rectangle 5 in Figure 2d, the fraction of energy scattered to the non-phase locked tide Figure 2b of (Buijsman et al. 2017) is comparable to the energy scattered to the supertidal band ($20-40\%$).

5. Summary and Conclusions

In this study we have for the first time diagnosed the depth-integrated and time-averaged supertidal energy budget in a realistically forced global HYCOM simulation. We have band-passed hourly 3D velocity and density fields over a 30-day period into diurnal (D1), semidiurnal (D2), and higher-harmonic (HH; supertidal) frequency bands and computed various energy terms. Depth-integrated kinetic energy and fluxes reveal horizontal internal tide beams at both tidal and supertidal frequencies. Tidal energy and fluxes are highest at internal tide generation sites and decay with distance from the source. However, supertidal energy and fluxes can peak a few hundred kilometers from the generation sites. The fluxes are dominated by the pressure flux, with the flux due to self advection being several orders of magnitude smaller for both the tidal and supertidal bands. At the strongest generation sites of super tidal energy, such as the Bay of Bengal, Amazon Shelf, and Mascarene Ridge, supertidal $KE$ constitutes about 25-50% of the $KE$ summed over the tidal and supertidal bands. The zonally integrated supertidal kinetic energy is highest near the equator, where it accounts for about 10% of the total internal tide energy, and decreases towards higher latitudes. We estimate that on average, in the deep ocean ($>250$ m) and equatorward of $\pm {25}^{\circ }$, supertidal $KE$ accounts for about 5.4 PJ ($\sim 5\%$) of the total internal tide energy.

As opposed to the tidal flux divergence, the supertidal flux divergence does not correlate with the barotropic to baroclinic energy conversion, which is elevated at topography in the tidal band and near zero in the supertidal band. Instead, the time-mean and depth-integrated supertidal flux divergence correlates well with the nonlinear kinetic energy transfers from (sub)tidal to supertidal frequency bands as estimated with a novel coarse-graining approach. At the largest supertidal energy generation sites, both supertidal flux divergence and nonlinear energy transfers are organized in regularly spaced banding patterns with separation distances larger than a semidiurnal mode-1 wavelength. These banding patterns are due to the constructive interference between semidiurnal mode 1 and 2 internal waves. The faster mode 1 wave overtakes the slower mode 2 wave, causing alternating patterns of reduced and elevated surface kinetic energy. Similar banding patterns have been observed in model simulations and satellite altimetry of the Amazon Shelf region by Tchilibou et al. (2022). These banding patterns are not present in $KE$ after depth-integration. We estimate that equatorward of $\pm {25}^{\circ }$ about 45 GW, or about 7% of the tidal barotropic to baroclinic energy conversion, is transferred to higher-harmonics. This number is comparable to the net energy transfer out of the semidiurnal internal tide due to PSI (Ansong et al. 2018).

Areas in the open ocean with simulated supertidal energy fluxes >1 kWm^-1 correlate well with observations of nonlinear solitary waves extracted from sunglint images (Jackson 2007). We demonstrate for the Amazon shelf, that as the tidal energy flux decays, the supertidal energy flux increases offshore. At the same time, the internal tide steepens into solitary waves, which surface expression agrees with SAR observations (Magalhaes et al. 2016). These simulated and observed solitary waves begin to appear near the second band of elevated nonlinear energy transfers at about 400 km from the shelf.

We acknowledge that our 4-km hydrostatic HYCOM simulation is only able to show the onset of the superharmonic energy cascade that affects the low-mode internal tide. Higher resolution simulations are needed to properly simulate the spectral slopes at higher frequency and wavenumbers (Nelson et al. 2020). Furthermore, these simulations need to be nonhydrostatic to accurately simulate the nonhydrostatic solitary waves (Vitousek and Fringer 2011), which is currently not yet feasible at global scales.