## Abstract

A key ingredient of energetically consistent ocean models is the parameterized link between small-scale turbulent mixing, an important energy source of large-scale ocean dynamics, and internal gravity wave energetics. Theory suggests that this link depends on the wave field’s spectral characteristics, but because of the paucity of suitable observations, its parameterization typically relies on a model spectrum [Garrett–Munk (GM)] with constant parameters. Building on the so-called “finestructure method,” internal gravity wave spectra are derived from vertical strain profiles obtained from Argo floats to provide a global estimate of the spatial and temporal variability of the GM model’s spectral parameters. For spectral slopes and wavenumber scales, the highest variability and the strongest deviation from the model’s canonical parameters are observed in the North Atlantic, the northwest Pacific, and the Southern Ocean. Internal wave energy levels in the upper ocean are well represented by the GM model value equatorward of approximately 50°, while they are up to two orders of magnitude lower poleward of this latitude. The use of variable spectral parameters in the energy level calculation hides the seasonal cycle in the northwest Pacific that was previously observed for constant parameters. The global estimates of how the GM model’s spectral parameters vary in space and time are hence expected to add relevant detail to various studies on oceanic internal gravity waves, deepening the understanding of their energetics and improving parameterizations of the mixing they induce.

## 1. Introduction

Internal gravity waves constitute a crucial component of the ocean’s energy cascade from large to small scales and play an important role for ocean dynamics as the mixing induced by their breaking provides some of the energy required to maintain the overturning circulation (Munk and Wunsch 1998; Wunsch and Ferrari 2004). The scale-bridging energy transfer is associated with nonlinear wave–wave interactions, which flux energy through the spectrum without changing its total amount (e.g., Müller et al. 1986). Neither this transfer nor the subsequent mixing can be resolved by state-of-the art ocean models, but because of their relevance for large-scale ocean dynamics, reliable parameterizations are indispensable. Their development requires a thorough understanding of internal gravity wave generation, interaction, and dissipation mechanisms and hence knowledge of the wave field’s spectral properties. This study presents the first global observational reference of the oceanic internal wave field’s spectral characteristics.

The dominant energy source of internal gravity waves is the interaction of barotropic currents with the rough seafloor, generating wave disturbances of tidal frequency, the so-called “internal tides” (e.g., Bell 1975; Garrett and Kunze 2007). Another process creating internal waves of specific frequency is a fluctuating wind stress, leading to resonance of the surface mixed layer at frequencies close to the local Coriolis frequency and to the emission of internal waves of near-inertial frequencies into the stratified ocean below [see, e.g., Alford et al. (2016) for a recent review]. Other sources with a somewhat lower global energy input are the generation of lee waves by geostrophic flows over rough bottom topography (Nikurashin and Ferrari 2011; Trossman et al. 2013), the resonant interaction of surface waves (Olbers and Eden 2016; Haney and Young 2017), and spontaneous generation (Vanneste 2013, and references therein). The typical imprint of a specific forcing on the internal wave field is lost after a certain distance, because processes like wave–mean flow interaction (Müller 1976), topographic scattering (Müller and Xu 1992), and wave–wave interactions (Olbers 1976) affect the local spectral energy level. The latter process is typically grouped into three classes, which dominate the energy transfer in different regions of the internal waves’ wavenumber–frequency spectrum (McComas and Bretherton 1977). At critical layers, where the waves’ frequency approaches the local inertial frequency, and through overturning induced by gravitational or shear instability (e.g., Olbers et al. 2012), internal wave energy is dissipated—sometimes hundreds of kilometers away from their generation sites. This energy is thought to be either a direct (McComas and Müller 1981a) or indirect (Kunze 2019) source of small-scale turbulence, which can create potential energy through density mixing.

Garrett and Munk (1972) combined available observations with linear wave kinematics to develop a model of the energy spectrum of a horizontally isotropic and vertically symmetric wave field. This Garrett–Munk (GM) model, including the later refinements by Garrett and Munk (1975), Cairns and Williams (1976), and Munk (1981), describes a universal spectrum, which develops after the spectral signal of a specific forcing has already been lost through the influence of the wave field’s interactions with surrounding features or other waves. In regions of strong forcing, the internal wave spectrum is hence expected to deviate from the GM model; everywhere else, however, this model is considered a good albeit somewhat smoothed representation of the ocean’s internal wave field (Briscoe 1977; Müller 1976; Olbers 1983). It is therefore widely used as a reference model, although the spectral parameters such as energy level or spectral slope are expected to depart from the model constants depending on the local conditions. The regional catalogue of Polzin and Lvov (2011), who fit the GM model to fine- and microstructure datasets from various locations, illustrates the geographic variability of these parameters. The increased awareness of the role of internal gravity waves for the ocean’s energy budget and of the need to adequately represent this in numerical models (e.g., Polzin 2009; Müller and Natarov 2003; Olbers and Eden 2013), however, also led to an increased need for a similar overview on global scales and of higher spatial and temporal resolution. By fitting the GM model to spectra obtained from vertical profiles collected by the Argo program, which since its introduction in the year 2000 has collected over two million vertical profiles in the ice-free ocean’s upper 2000 m (e.g., Riser et al. 2016), this study provides such a reference. In the following section, the theoretical background and the method are described. Global maps of the obtained results as well as a presentation of their vertical and seasonal variability are presented in section 3, while section 4 is an in-depth discussion of the uncertainties associated with this study’s approach combined with concluding remarks.

## 2. Data and methods

### a. The Garrett–Munk model and spectral relations

GM-type energy spectra in vertical wavenumber *m* and frequency *ω* space can be factorized into the local energy level *E* and shape functions *A*(*m*) and *B*(*ω*)

with

where *f* is the Coriolis frequency and *N* is the buoyancy frequency. The scaling functions *n*_{A} and *n*_{B} ensure that the functions *A* and *B*, respectively, integrate to unity and ∫*S*_{E}(*m*, *ω*) *dm dω* = *E*. Equations (1)–(3) describe the modified form of the original GM model (Garrett and Munk 1972) suggested by Cairns and Williams (1976), with no cutoff at low wavenumbers and spectral parameters slope *s*^{GM} = 2, wavenumber scale $m*GM=0.01\u2009m\u22121$, and energy level *E*^{GM} = 3 × 10^{−3} m^{2} s^{−2}. In the GM model, the energy level *E*^{GM} is furthermore scaled by *N*(*z*)/*N*_{0}, where *N*_{0} is a constant reference buoyancy frequency; for the estimation of energy levels from Argo profiles, however, no such scaling will be assumed. The vertical wavenumber bandwidth *m*_{b} is defined as

and related to the vertical wavenumber scale $m*$ by

Quantities typically observed by standard instrumentation are vertical shear **u**_{z}, where **u** is the two-dimensional horizontal velocity, and vertical strain $\xi z=(N2\u2212Nfit2)/N2\xaf$, where $Nfit2$ is a quadratic fit to the vertical profile of *N*^{2} and $N2\xaf$ the average. The polarization relations link spectra of these quantities, $Suz\u2061(m,\u2009\omega )$ and $S\xi z\u2061(m,\u2009\omega )$, to the energy spectrum *S*_{E}(*m*, *ω*):

Integration over all internal wave frequencies *f* ≤ *ω* ≤ *N* gives vertical wavenumber spectra $S\xi z\u2061(m)$ and $Suz\u2061(m)$, respectively (see also Pollmann et al. 2017):

with

### b. Determining internal wave saturation

Vertical wavenumber spectra of strain can be obtained from vertical profiles of temperature, salinity, and pressure (CTD measurements). Such observations have been collected in large numbers almost everywhere in the ice-free ocean by Argo floats, which profile the ocean’s upper 2000 m roughly every 10 days at a resolution of a few meters (Argo 2000). This implies that motions that do not necessarily represent internal gravity waves but, for example, their transition to turbulence are resolved in the Argo data. For the analysis of internal wave spectral properties it is therefore essential to define an upper wavenumber threshold, above which spectral information is discarded.

The GM model involves a cutoff at *m* = *m*_{c} = 2*π*/10 m^{−1}, up to which the vertical wavenumber shear spectrum is basically flat and contains a total shear of

[for a wavenumber dependence of $\u2061(m+m*)\u22122$, the right-hand side amounts to 0.7*N*^{2}; see e.g., Polzin et al. 1995]. This indicates that the internal wave field is associated with Richardson numbers $Ri\u22121=\u2329\u2061(uz)2\u232a/N2~O\u2061\u2009(1)$ (Munk 1981) and that the upper wavenumber cutoff *m*_{c} can be defined as that wavenumber up to which the shear spectrum needs to be integrated for its variance to reach a critical value (e.g., Polzin et al. 1995). The corresponding threshold for the strain variance can be found by specifying the shear-to-strain ratio, which is defined as

it amounts to $R\omega GM=3$ for the GM model and is observed to be *R*_{ω} = 7 ± 3 in the global analysis of Kunze et al. (2006). In consequence, the maximum strain variance is typically set to 0.1 (Kunze et al. 2006; Waterman et al. 2014) or 0.2 (Whalen et al. 2012, 2015). It has been argued (e.g., Munk 1981; Gargett et al. 1981; Gargett 1990) that if the wavenumber cutoff is to be associated with a fixed Richardson number, the product *Em*_{c} should be constant so that the cutoff *m*_{c} shifts to lower wavenumbers *m* < 2*π*/10 m^{−1} as the energy increases above *E*^{GM}. For higher wavenumbers *m*_{c} ≤ *m* ≤ *m*_{u} ≈ 10*m*_{c}, Munk (1981) and Gargett et al. (1981) suggest a different spectral shape proportional to (*m*)^{−1} (see also, e.g., Lindborg 2006; Kunze 2019). Based on simple stability arguments applicable to both atmospheric and oceanic internal gravity waves, such a wavenumber dependence is associated with spectral saturation (Orlanski and Bryan 1969; Munk 1981; Smith et al. 1987), that is, the situation when the waves are only just stable and added energy will trigger instability and thus enhanced energy loss (Thorpe 2007). For the GM model with an average shear spectral level of 0.5*N*^{2}/*m*_{c}, the wavenumber spectrum in this wavenumber range is given by (see also Kunze et al. 2006):

Analogously, for an observed variance ⟨(**u**_{z})^{2}⟩ distributed over a wavenumber range *m*_{2} − *m*_{1}, the wavenumber spectrum follows

Extending this wavenumber dependence to lower wavenumbers, i.e., the internal wave wavenumber range, produces so-called “saturation spectra” and provides an upper bound for the fitted Argo spectra.

These concepts introduce a variety of related but different criteria to determine spectral saturation and to identify an unsaturated wavenumber range in which the spectral analysis should be performed. A suitable saturation criterion for the purpose of this study is selected based on sensitivity tests of spectral parameters to these criteria combined with different, a priori defined wavenumber ranges for a subset of Argo data. Due to its high level of technicality, this sensitivity analysis is presented in the appendix. Its main outcome can be summarized as follows: the differences between the various test cases are rarely significant within the uncertainty of the method [Pollmann et al. (2017) find the finestructure method’s uncertainty to underlying assumptions and parameter settings can reach a factor of 5]. For practical reasons, a combination of wavenumber range and saturation criterion was chosen that ensures (i) a relatively high number of individual estimates, (ii) a parameter variance not too far above or below that of previously used and evaluated settings (Whalen et al. 2012, 2015), (iii) a heuristic upper level for internal wave energy, and (iv) not too broad fit confidence intervals. The last point is the reason why the previous approach of considering a lower wavenumber bound *m*_{1} = 2*π*/100 m^{−1} and an upper wavenumber bound *m*_{2} given by either 2*π*/10 m^{−1} or the wavenumber for which a maximum strain variance $\u2329\xi z2\u232a$ of 0.1 is reached (Kunze et al. 2006; Waterman et al. 2014; Pollmann et al. 2017) was not followed. Instead, the wavenumber range is *m*_{1} = 2*π*/100 m^{−1} to *m*_{2} = 2*π*/100 m^{−1} and only spectra with at least three points within this range and a corresponding shear variance below *N*^{2} are considered [i.e., saturation is associated with an inverse Richardson number Ri^{−1} = ⟨(**u**_{z})^{2}⟩/*N*^{2} = 1].

The shear variance needed for this calculation is estimated from the internal wave energy level *E*, which is obtained from observed strain spectra using Eq. (8) and is then used to estimate shear spectra based on Eq. (9). This procedure relies on the assumptions that the local internal wave field is adequately represented by the observed strain spectrum and that its form, in particular its frequency dependence, does not depart too much from that of the GM model, that is, similar assumptions as those underlying the approach of defining an upper strain variance limit to identify the transition to saturation. Shear variance is obtained by integrating these spectra over the same wavenumber range as chosen before in the calculation of strain variance and internal wave energy from Eq. (8), that is, *m*_{1} = 2*π*/100 m^{−1} to *m*_{2} = 2*π*/10 m^{−1}.

### c. Spectral fits and TKE dissipation rates

Using data collected from January 2006 to December 2017 and following the procedure detailed in Pollmann et al. (2017), vertical wavenumber spectra of strain, $S\xi z\u2061(m)$, are calculated from half-overlapping 200-m segments of strain profiles with an average vertical resolution not coarser than 10 m (see also, e.g., Kunze et al. 2006; Whalen et al. 2012, 2015). From each individual spectrum, the spectral slope *s* and the wavenumber scale $m*$ are estimated from a least squares fit using MATLAB’s trust-region-reflective algorithm. The starting points are given by the GM values *s*^{GM} = 2 and $m*GM=0.01\u2009m\u22121$, and the upper and lower bounds for the fitted parameters are a factor of 20 around these starting points except for the lower bound of *s*, which is *s*_{min} = 1.001. The fitting iteration is continued until convergence within numerical precision is achieved. The internal wave energy level *E* is then calculated by inserting the fitted $m*$ and *s* into Eq. (8) and integrating between *m*_{1} = 2*π*/100 m^{−1} and *m*_{2} = 2*π*/10 m^{−1}. This integration produces the strain variance $\u2329\xi z2\u232a$ on the left-hand side, which is obtained by summing the squared Fourier amplitudes; the integration on the right-hand side is done numerically. Inserting these energy levels *E* into Eq. (9) and integrating over the same wavenumber range gives an estimate of the associated shear variance, which is used to identify whether the spectrum in question covers the unsaturated part of the internal wavenumber spectrum only, assuming saturation at an inverse Richardson number Ri^{−1} = ⟨(**u**_{z})^{2}⟩/*N*^{2} = 1. If the spectra pass this test, the results are kept and sorted into depth bins based on the segments’ average depth.

The strain variance $\u2329\xi z2\u232a$ is a key component of the finestructure method proposed by Gregg (1989) to estimate the dissipation rate of turbulent kinetic energy (TKE) *ε*_{TKE} from observations at scales of $O \u2061(10)$ m. The expression used here is the result of extensive evaluation and advancement (e.g., Wijesekera et al. 1993; Polzin et al. 1995, 2014) and is given by

with *ε*_{0} = 6.73 × 10^{−10} W kg^{−1} and *N*_{0} = 5.24 × 10^{−3} s^{−1} (e.g., Kunze et al. 2006). The GM model variance $\u2329\xi zGM2\u232a$ is calculated by integrating Eq. (8) over the same wavenumber range as chosen for the observed strain spectrum. The function *L*_{f}(*f*, *N*) is a latitudinal correction

where the subscript “30” denotes the value at a latitude of 30°. The function *h*(*R*_{ω})

accounts for the frequency content of the wave field and reduces to unity when the shear-to-strain ratio *R*_{ω} is set to the GM value of 3. This is typically done in cases when only strain information is available, such as for Argo CTD profiles (Whalen et al. 2012; Pollmann et al. 2017). Assuming the spectral forms and relations described in the previous section, in principle allows the calculation of the shear-to-strain ratio *R*_{ω}—since it is not based on actual shear measurements and relies on integrated GM frequency spectra, the notation $R\omega I$ is introduced. To investigate the influence of $R\omega I$, which mainly represents latitude and stratification effects, TKE dissipation rates are calculated from Eq. (16) in two ways, using either $h\u2061(R\omega GM)=1$ or $h\u2061(R\omega I)$. This is meant to be a theoretical exercise only; whether $h\u2061(R\omega I)$ can have similar effects as *h*(*R*_{ω}) requires a thorough evaluation of shear and strain measurements under various dynamical conditions and will be addressed in a different study.

## 3. Results

Following the procedure detailed above, individual spectral fits are calculated and sorted into 100-m-depth bins and 1° × 1° horizontal bins. In each bin, the results are averaged over these individual estimates. Derived quantities like the internal wave energy and the TKE dissipation rate are calculated from individual results and then averaged. Figure 1 shows two examples of observed and fitted spectra as well as the geographic distribution of the number of individual estimates obtained in this study. Of the 913 bins with more than 400 estimates, 754 are located in the Pacific and 139 in the Indian Ocean.

### a. Global variation of internal wave spectral properties

Figures 2 and 3 show the global variation of the bin-averaged fitted quantities and derived internal wave energy levels in three vertical depth ranges (300–500, 500–1000, 1000–2000 m). The variation of the spectral slopes (Figs. 2a,c,e) is rather weak, with the highest slopes in the Southern Ocean, in the northwest Pacific, and at high latitudes in the North Atlantic. Large areas of intermediate to low values are confined to low latitudes, but individual values well below the GM model value can also be found elsewhere. The geographic differences are most pronounced above 500 m and weakest in the middle depth range, where the imprint of the wave field’s main forcing mechanisms (winds, tides) is weakest. The location of slope maxima persists with depth, but these are sometimes less distinct (the northwest Pacific, for example, can only be identified as a region of higher spectral slopes in the upper and the lower depth range). The geographic variation of the wavenumber scale $m*$ (Figs. 2b,d,f) is similar, with maximum values observed in the North Atlantic, the northwest Pacific and in the Southern Ocean. The horizontal gradients are again strongest in the upper and weakest in the middle depth range. With a few exceptions in particular in the Indian Ocean and around the equator, the wavenumber scale at low latitudes is relatively small and increases with depth.

The internal wave energy levels (Figs. 3a,c,e) agree with the results published in Pollmann et al. (2017), who assumed a spectral slope of *s*^{GM} = 2 and a wavenumber scale of $m*GM=0.01\u2009m\u22121$: there is a strong equator-to-pole gradient of two orders of magnitude, the maxima coincide with the subtropical gyres of the Pacific Ocean, and the North Atlantic Current as well as parts of the Antarctic Circumpolar Current (ACC) feature elevated energy levels. The band of high energy at the equator, however, is not seen in the present analysis and particularly in the central ocean basins, the energy is somewhat reduced. This is likely due to the nonconstant spectral parameters and the choice of a different saturation criterion, which discards spectra with increased levels at high wavenumbers more often than the specification of an upper strain variance limit would (note that spectral slopes lower than *s*^{GM} = 2 indicate increased strain spectral levels at higher vertical wavenumbers, see also Fig. 1a). The bandwidth *m*_{b} (Figs. 2b,d,f) is derived from the wavenumber scale $m*$ and the spectral slopes *s* according to Eq. (5) and varies by two orders of magnitude between equator and poles in the upper ocean. This latitudinal gradient is reduced with depth and not symmetric around the equator: in the east Pacific, elevated values are found at all latitudes, while in the west Pacific, these are confined to middle to high latitudes.

Figure 4 illustrates the agreement with the GM model values. Wavenumber scales $m*$ within a factor of 2 of the corresponding GM reference are set to $m*GM=0.01\u2009m\u22121$, which affects about 40% of the shown averages. Since about 98% of the bin-averaged spectral slopes agree with *s*^{GM} = 2 within a factor of 2, only values between 1.8 and 2.2 were marked as the GM reference, which concerns 50% of the averages. Energy levels, too, agree very well with the buoyancy scaled GM reference in the upper depth range, so that in Fig. 4a, values within 30% of 3 × 10^{−3}*N*/*N*_{0} m^{2} s^{−2} are set to *E*^{GM} = 3 × 10^{−3} m^{2} s^{−2}, which affects about half of the shown averages. Considering all available estimates, 15%–20% of the wavenumber scales, 50%–60% of the energy levels, and about 25% of the spectral slopes agree with the GM model reference within a factor of 2, within 30%, or within 10%, respectively, with the range denoting ocean basin differences (for wavenumber scales, the lower bound represents the Atlantic and the upper bound the Pacific, while the opposite holds for energy levels).

The figure highlights that the buoyancy scaled GM energy level (Fig. 4a) is a good reference equatorward of approximately 50° (depending on ocean basin and hemisphere), with only a few exceptions in parts of the Gulf Stream and around the equator (too high) as well as near 20° in the Pacific (too low). Poleward of approximately 50°, the buoyancy scaled GM energy level is rarely reached. Spectral slopes (Fig. 4b) in very good agreement with the GM model value can be observed at all latitudes, although the North Atlantic, the ACC, and the northwest Pacific as well as the equator emerge as regions with higher or lower slopes, respectively. The GM wavenumber scale $m*GM$ is a good representation of the observed $m*$ at high and middle latitudes. In particular in the northwest Pacific and in the Indian Ocean, good agreement can also be found closer to the equator. There is no clear relation between how well spectral slopes and wavenumber scale on the one hand and energy levels on the other agree with the respective GM model values, underlining that the local forcing, wave processes, and the surrounding environment influence these spectral parameters differently. In the lower depth range (not shown), the agreement with the GM reference is notably reduced for energy levels, with only 3% of the results reaching within 30% of *E*^{GM} = 3 × 10^{−3}*N*/*N*_{0} m^{2} s^{−2} below 1000 m (40% reach within a factor of 2). For the spectral slopes, the agreement slightly improves with depth, while the most estimates agreeing with the GM wavenumber scale $m*GM$ within a factor of 2 are found between 500 and 1000 m. The estimate distributions (not shown for all data, but for a subset of data in Fig. 8) demonstrate that estimates not agreeing with the respective GM reference are mostly lower. The distributions’ positive skew and taking the mean of a large number of individual spectra produce the observed agreement, which is in line with the model’s aim of representing the average internal wave field under no direct influence of forcing. It also underlines that to investigate the variation of the model parameters as done here, it is necessary to average over many individual estimates, which is made possible by the wealth of profiles collected by the Argo program.

### b. Global variation of TKE dissipation rates

TKE dissipation rates are calculated from Eq. (16) based on the spectra used for the spectral fits presented in the previous section. The geographic variation (see Fig. 5a) generally agrees with that published in previous studies (Whalen et al. 2012, 2015; Pollmann et al. 2017, which are based on fixed spectral parameters from different GM model versions) with high TKE dissipation rates over rough topography (e.g., Hawaii or Indonesia), in western boundary currents (e.g., Kuroshio, Gulf Stream), in the wind-driven gyres of the Pacific Ocean, and in regions of high eddy kinetic energy (e.g., Agulhas retroflection, ACC). Local deviations are likely caused by (i) the larger database of Argo profiles used here and (ii) methodological differences (see also section 4b).

A striking, overarching difference is that with the exception of the dissipation hot spots, the magnitude of TKE dissipation is strongly reduced compared to these previous studies and that only few regions of elevated dissipation rates are present below 500 m (see Fig. 5b). The difference between relating internal wave saturation to a maximum shear variance of *N*^{2} (Fig. 5a) and to that of 0.5 *N*^{2} (Fig. 5d) is in many areas much smaller than the deviation from previously published results. This suggests that an increase of the potentially too strict upper shear variance limit would not result in a better agreement with these previous estimates (while remaining realistic). Figure 5e shows the TKE dissipation rates obtained when using the same wavenumber range and saturation criterion as in Pollmann et al. (2017), i.e., *m*_{1} = 2*π*/100 m^{−1} to *m*_{2} = 2*π*/10 m^{−1} while reducing *m*_{2} if necessary to satisfy $\u2329\xi z2\u232a\u22640.1$, contrasting the results with those obtained by Pollmann et al. (2017) for the same input data and in the same depth range in Fig. 5f. Hence, the difference between Figs. 5e and 5f lies in the constraint, that the unsaturated part of the spectrum considered for the integration and the fits is required to contain at least three wavenumber points—a criterion found necessary to produce reliable spectral fits (see section 4) but not applied in the analysis shown in Fig. 5f. Because of the smaller number of input data (to save computational time, only the years 2006–11 were considered), the spatial coverage is reduced compared to that in the other subfigures, but there are enough data points to illustrate that apart from the hot spots near Hawaii, Indonesia, and in the northwest Pacific, TKE dissipation rates are much lower in Fig. 5e. It is noteworthy that this strong sensitivity to *n*_{min} is only observed for TKE dissipation rates, not for energy levels (not shown). This could be explained by the TKE dissipation rates’ dependence on the GM model strain variance, whereas energy levels are estimated based on the variable spectral parameter fits (see section 4b and the appendix for further discussion). The additional constraint *n*_{min} = 3 does not affect the results’ sensitivity to topographic roughness: TKE dissipation rates of both Pollmann et al. (2017) and of this study are on average higher above rough^{1} topography than above smooth topography, with the largest differences observed equatorward of 30° and at 65°–70° latitude (not shown) similar to the behavior reported by Whalen et al. (2012, Fig. 3a).

Eden et al. (2019) evaluate resonant and nonresonant nonlinear interactions between internal gravity waves of GM-type form and find that their energy dissipation *ε*_{IW} depends on the spectral slope according to

which extends the parameterization of *ε*_{IW} suggested by McComas and Müller (1981b) (see also Olbers 1976; Henyey et al. 1986) by a factor (*s* − 1)^{−3}. Without this slope dependence, Eq. (19) forms the basis of the fine-structure method (Gregg 1989). Replacing the shear production term in the TKE balance by *ε*_{IW} and assuming a fixed balance of shear production, turbulent buoyancy flux, and TKE dissipation (Osborn 1980) *ε*_{TKE} can be estimated from Eq. (19), when furthermore a fixed ratio of turbulent buoyancy flux and TKE dissipation is assumed (Olbers and Eden 2013):

where Γ ≈ 0.2 is the so-called “mixing efficiency.”

The TKE dissipation rates calculated from Eqs. (19) and (20) for *μ*_{0} = 0.6 are shown in Fig. 5c. There is a noteworthy increase of TKE dissipation at middle and high latitudes and a minor decrease at the hot spots (e.g., northwest Pacific, Indonesia) compared to Fig. 5a. This leads to a better agreement with the results of Whalen et al. (2012, 2015) and Pollmann et al. (2017), although discrepancies remain with respect to the TKE dissipation rate magnitude, especially near the equator, and with respect to the spatial extent of dissipation hot spots (e.g., around Hawaii or in the South Pacific Gyre, elevated dissipation rates cover a substantially larger area in these previous studies).

### c. Seasonal variation

Whalen et al. (2012) showed that TKE dissipation rates have a distinct seasonality in the northwest Pacific. This section addresses the question whether this holds true for the fitted spectral properties. Figure 6 shows the time series of monthly mean energy levels, wavenumber scales, spectral slopes, and TKE dissipation rates in 5° latitude bands between 150°W and 150°E and between 300 and 400 m depth. To avoid biases due to different numbers of estimates in different parts of the domain, the parameters are first averaged onto a regular latitude–longitude grid before computing the average time series for the entire domain. As expected from previous work (e.g., Jing and Wu 2010; Whalen et al. 2012), the TKE dissipation rates (Fig. 6b) feature a clear seasonal cycle and a noteworthy latitudinal dependence with values increasing equatorward. A similar behavior is seen in the internal wave energy time series (Fig. 6a), with increased variability as well as higher amplitudes closer to the equator. There is, however, no clear cycle in the annual variability. This holds true for the spectral slopes (Fig. 6c), but not for the wavenumber scales (Fig. 6d), which in particular during the first part of the considered period feature clear summer minima and winter maxima south of 40°. Note that the time series of wavenumber scales and spectral slopes are offset in the figure for clarity as the amplitudes are comparable in all latitude bands.

Seasonal variability in TKE dissipation rates and diapycnal diffusivities has previously been linked to similar variations of near-inertial energy in the mixed layer (Jing and Wu 2010; Whalen et al. 2012). This implies that the seasonal variability in the wind forcing causes variations of near-inertial internal wave generation, which in turn lead to corresponding variations in wave energy dissipation and mixing. The lack of a distinct seasonal cycle in the internal wave energy levels is at odds with this picture. It is also at odds with the time series produced from the energy levels of Pollmann et al. (2017), which were derived assuming *s* = 2 and $m*=0.01\u2009m\u22121$ and feature a distinct seasonal cycle in the domain in question (not shown). Taking the variability of spectral slope *s* and wavenumber scale $m*$ into account instead of using constant GM model values hence removes the seasonality of the derived energy levels. It also explains why time series of TKE dissipation rates (which only depend on GM model spectral constants) behave differently from those of energy levels (which depend on variable spectral parameters): there is a pronounced seasonal dependence in the time series of vertical strain, which forms the basis of both TKE dissipation rate and energy level estimates. This seasonality remains discernible after scaling by the GM model variance and multiplying by $N2/N02$ to estimate the TKE dissipation rates. It is, however, almost completely erased when the vertical strain is divided by the integrated wavenumber-dependent part of Eq. (8), that is, $\u222bm1m2m2/{m*\u2061[1+\u2061(m/m*)s]}\u2009dm$, in the process of estimating energy levels (not shown). In conclusion, estimates of how the internal waves’ energy levels vary in time depend substantially on whether the wave spectrum is assumed to be fixed or allowed to vary. As illustrated in Figs. 6c and 6d, the vertical wavenumber slope *s* and the wavenumber scale $m*$ of the internal wave energy spectrum vary on seasonal and other time scales. Also in frequency space, a temporal variability of the spectral parameters seems possible, as variations in frequency spectral slope have been linked to stratification changes and the associated modifications of the coupling between near-inertial and tidal motions (van Haren et al. 2002). Such variations of the spectral shape in turn affect mechanisms like the nonlinear energy transfer through the spectrum, a prerequisite for wave-induced turbulent mixing. Its temporal variability and that of the internal wave energy levels might be understood in terms of seasonal changes in wave generation alone when only a narrow internal wave frequency band is considered [Alford and Whitmont (2007) observe a pronounced seasonality of near-inertial kinetic energy consistent with local wind forcing variability]. Owing to the lack of frequency information in the vertical Argo profiles, this study represents all internal wave frequencies present when and where a profile was collected. These include internal tides, for which a seasonal cycle is observed in some diurnal constituents but barely in the most energetic, semidiurnal M_{2} tide (Müller et al. 2012; Xu et al. 2013), and the higher-frequency continuum generated for example by nonlinear or kinematic interactions. The disappearance of the seasonal cycle in the internal wave energy levels for variable spectral parameters underlines, that the seasonality of the forcing can only partly account for the temporal variability of the waves’ energy when a broad frequency range is considered. In addition, variations in the spectral form and in wave interaction processes (with possible reciprocal effects; van Haren et al. 2002) as well as in the environmental parameters that shape them have to be recognized.

All time series in Fig. 6 display a notable reduction in variability and average amplitude in all latitude bands beginning in the year 2012. This result is robust to methodological modifications: when using only delayed-mode data, which involves a higher degree of quality control and for example flags profiles with strong salinity biases, or computing strain from potential temperature rather than buoyancy frequency, the time series barely change (not shown). Figure 6b includes a time series of eddy kinetic energy (EKE) at 30°–35°N in the same domain, which features higher intra-annual variability and an increasing trend until early 2012 and lower intra-annual variability combined with a decreasing trend afterward, analogous to the behavior of energy and TKE dissipation rate levels in this latitude band. Jing and Wu (2014) and Whalen et al. (2018) demonstrate that strong mesoscale activity increases the seasonal variability of turbulent mixing. The observed long-term trends in internal wave energy and TKE dissipation rates agree with this argument. It is likely that other processes play a role, too, seeing that the trends in EKE prior to and after 2012 are comparable, which is not the case for most time series shown in Fig. 6. As the period of available Argo observations only just reaches two decades, hypotheses about the impact of large-scale, lower-frequency phenomena like El Niño–Southern Oscillation on the local wave field and on turbulent mixing remain, however, pure speculation.

### d. Variation with depth

The vertical variation of spectral properties and TKE dissipation rates is depicted in Fig. 7, averaged in a 4° latitude band through Luzon Strait (18°–22°N) in the west Pacific. The black line in Fig. 7a denotes the minimum water depth in the averaging region calculated from Becker et al. (2009) and marks Taiwan and Luzon around 120°E, the Northern Mariana Islands near 145°E, and the individual seamounts east of the Mariana Trench. The internal wave energy Fig. 7a is unaffected by these topographic variations and smoothly decreases with depth, which implies a strong influence of the stratification. Spectral slopes *s* (Fig. 7b) and wavenumber scales $m*$ (Fig. 7c) exhibit a rather high variability around intermediate and high values near the surface, low values with small vertical and horizontal differences between 400 and 1000 m, and again higher values and gradients below 1000 m. TKE dissipation rates (Fig. 7d) generally decrease with depth, but with somewhat higher horizontal variation than seen in the internal wave energy (Fig. 7a).

A noteworthy feature is the patch of high spectral slopes, wavenumber scales, and TKE dissipation rates near 160°E around 1000–1400 m depth. While the line of minimum water depth (Fig. 7a) suggests high topographic variability not only in this area but also farther east and west, maps of topographic roughness (e.g., Decloedt and Luther 2010) show that this parameter is indeed patchy in the northwestern Pacific, implying a similarly irregular influence on internal wave generation, dissipation, and wave-induced mixing. Near 160°E, the seamount density is particularly high in the latitude band in question and, contrary to most regions of the global ocean, the majority of the Argo profiles represent the entire water column. Near the ocean bottom (in most cases not deeper than 1400 m) they show enhanced TKE dissipation rates and often also enhanced spectral parameters. This is in line with numerous previous studies reporting enhanced mixing over rough topography, indicative of internal wave generation through tide–topography interaction and turbulent mixing induced by their subsequent instability and breaking (e.g., Polzin et al. 1997; Klymak et al. 2008; Tian et al. 2009). It is the binning and averaging over four latitude bins of the individual results that produces the square shape of the anomaly and is responsible for concealing the bottom enhancement observed in the vicinity of other sea mounts in the west Pacific in bins including also a substantial fraction of smooth sea floor.

Near 120° and near 140°, intermediate TKE dissipation rates are observed deeper in the water column than in the surrounding areas, but without a corresponding signal in $m*$ and *s*. The pattern of increased values between 300 and 700 m near 120° and lower values west of it is mirrored in the internal wave energy (Fig. 7a), and energy levels higher than 10^{−3} m^{2} s^{−2} can be found deeper in the water column near 140° and 165° than in the surrounding regions. In particular the increase near 165° is, however, insignificant within the uncertainty of the method (contrary to the increase of TKE dissipation rates). This underlines, that the wave field’s energy level, its spectral form, and its energy dissipation are influenced differently by local characteristics such as seafloor roughness owing for example to scattering processes (Müller and Xu 1992) or the interaction with remotely generated internal waves.

## 4. Discussion

Identifying suitable strain spectra to estimate internal wave spectral properties as is done here is based on a variety of assumptions, mainly that 1) the observed variability is caused by internal gravity waves and that 2) the observed wave field does not depart too much from the GM-type form. Assumption 2 implies that fits to this spectral form, the application of the polarization relations, and the application of saturation criteria, each of which involves uncertainties of its own, can be justified within a reasonable uncertainty.

### a. Uncertainty caused by nonwave sources

To meet assumption 1, profiles (partially) spanning the mixed layer and mode waters are discarded by applying the variable temperature criterion (de Boyer Montégut et al. 2004) twice following Whalen et al. (2012). Moreover, profiles with too strong variations in buoyancy frequency are discarded to avoid overestimating the strain variance when a segment’s buoyancy profile is poorly represented by a quadratic function [see also Whalen et al. (2012, 2015) for a more detailed motivation]. There remains a risk that nonwave processes of similar scales as those assumed for the internal wave field (10–100 m wavelengths) add to the observed variability and incorrectly increase the spectral estimates. A prominent example is double-diffusive convection (DDC), which arises in a hydrostatically stable density profile due to the different molecular diffusivities of heat and salt and takes the form of salt fingers or diffusive layers depending on the sign of the vertical gradient of temperature and salinity (Stern 1960; Schmitt 1994). Its strength is described by the Turner angle (Tu; Ruddick 1983),

where *α* and *β* are the coefficients of thermal expansion and haline contraction, respectively, Θ is Conservative Temperature, and *S*_{A} is Absolute Salinity. DDC is characterized by Turner angles between −45° and −90° (diffusive regime) and between +45° and +90° (saltfinger regime). For further distinction, Tu = ±60° and Tu = ±75 mark the transition from weak to intermediate and from intermediate to strong DDC. For −45° < Tu < 45°, the water column is stable to this type of instability. A global atlas of the Turner angle was presented by You (2002).

The influence of DDC on the estimated internal wave spectral properties is investigated based on Turner angle classifications for a subset of data (years 2011–15). Like buoyancy frequency, the Turner angle is computed using the Gibbs SeaWater Oceanographic Toolbox based on TEOS-10 standards (McDougall and Barker 2011). Figure 8 shows the estimate distribution of spectral slopes and internal wave energy levels for this subset of data, distinguishing four different cases: 1) using all estimates independent of Tu as done in the full data analysis, 2) discarding segments with at least 50% of their Tu estimates in the DDC range, 3) considering only segments with at least 50% of their Tu estimates in the DDC stable range, 4) considering only segments with at least 50% of their Tu estimates in the DDC range. The last case illustrates the influence of DDC on the spectral estimates, the first three serve to evaluate the uncertainty introduced by not correcting for DDC in the results. Since Turner angles characterize DDC favorability, not its actual occurrence, and the threshold of 50% to categorize individual segments is somewhat arbitrary, this has to be understood as a guideline rather than an error quantification.

In all ocean basins and for both energy levels and spectral slopes, the estimate distribution is comparable in cases 1, 2, and 3, which holds true independent of depth range and for wavenumber scales and TKE dissipation rates (not shown). In case 4, the maxima of the energy estimate distributions in all ocean basins, and, to a small degree, that of spectral slopes in the Indian Ocean are shifted toward higher values compared to cases 1–3. For spectral slope estimates in the Atlantic and Pacific Oceans as well as for TKE dissipation rates and wavenumber scales, the distributions remain comparable in all four scenarios (not shown). Especially in the Atlantic Ocean, the number of estimates per bin is strongly decreased when profiles favorable to DDC are discarded; in fact, the majority of the estimates (76%) represents profiles in the DDC favorable range (but only 6% of the estimates are obtained from profiles favorable to strong DDC). This is mirrored in the strongly reduced coverage of TKE dissipation rate estimates in the associated global map (Fig. 8h) compared to the reference scenario of using all estimates (Fig. 8g). For the remaining data points, there is no substantial change in the geographic variability: the northwest Pacific remains a hot spot, the ACC and the North Atlantic still feature intermediate to high values, and the equator, at least in the Pacific and Indian Oceans, remains an area of very low TKE dissipation rates. This is supported by the relative difference shown in Fig. 8i, where the majority of the differences is well below the uncertainty of the method. Removing profiles with at least 50% of the Tu values in the DDC favorable range mostly leads to reduced TKE dissipation rate estimates, but there are also exceptions (e.g., the Indonesian Seas or individual points in the Southern Ocean). The pattern is similar for $m*$, *s*, and *E*, albeit with slightly increased differences for $m*$ and reduced differences for *E* and *s*, much so for the latter (not shown). The overall decrease as well as the shifted maxima in the estimate distributions of energy when only DDC favorable profiles are considered indicates that DDC typically increases finestructure estimates of spectral properties and derived TKE dissipation rates, in line with the notion that nonwave variability is incorrectly assigned to internal waves. There are also cases, when the opposite is true, underlining that 200-m vertical segments can easily include parts of different degrees of DDC favorability and also different stages of DDC activity, such that the overall influence on the total variability is diminished. In conclusion, the effect of DDC on finestructure estimates could hence be relevant for individual profiles, but for the abundance of Argo profiles used in this study, the effect is judged insignificant within the uncertainty of the method.

### b. Uncertainty caused by underlying assumptions and methodological choices

The validity of assumption 2 can to some extent be assessed by the confidence intervals of the fit—the larger these are (see Figs. 9c,d), the harder it is to force an observed spectrum into the shape given by Eqs. (1) and (2). Narrower confidence intervals can be achieved by requiring that fitted spectra consist of a minimum number of points *n*_{min} (Figs. 10c,f), but in combination with an upper variance limit this constraint introduces an artificial bias toward spectra in which the variance is spread over a larger number of wavenumbers, i.e., the less energetic spectra. The strong decrease in TKE dissipation rates with increasing *n*_{min} in the r1a scenario (the upper wavenumber bound *m*_{2} is given by the largest wavenumber satisfying $\u2329\xi z2\u232a\u22640.1$ but at most 2*π*/10 m^{−1}, Fig. 10e) reveals that the maximum strain variance of 0.1 is typically reached after integrating between very few wavenumbers only. The difference of the TKE dissipation rates’ magnitude to those published in Pollmann et al. (2017), which do not rely on a minimum spectral length in the integration range, supports this conclusion (see Figs. 5e,f). That energy levels still agree well with previous results despite the specification of *n*_{min} = 3 suggests that the scaling with the GM model reference in the calculation of *ε*_{TKE} reinforces the bias introduced by this additional constraint (for *n*_{min} = 10, the spectral slopes depart only little from *s*^{GM} = 2 and the energy levels, too, decrease notably, see Fig. 10).

On the one hand, it is questionable whether spectra consisting of one or two data points only can adequately represent the local internal wave field and whether it can be determined with confidence that such spectra agree with the GM model shape. On the other hand, TKE dissipation rates estimated without a minimum spectral length (Pollmann et al. 2017) better agree with those of Whalen et al. (2015), which were shown to reproduce microstructure estimates within a factor of 2–3 in different parts of the ocean, than the TKE dissipation rate estimates of the present study. Whalen et al. (2012, 2015) indirectly set a minimum spectral length by requiring that the vertical resolution be no coarser than 15 m and by integrating the strain power spectral density between wavelengths of 100 and 40 m. If the integral (i.e., the strain variance) in that wavenumber range exceeds a threshold of 0.2, it is set to that value; if it does not, the integration is continued down to wavelengths of 10 m or until the strain variance threshold is reached (C. Whalen 2020, personal communication). The first part of this approach essentially limits the strain variance in a fixed integration range instead of, as done in Pollmann et al. (2017) and the r1 scenarios of this study (see appendix for details), restricting the integration range to a fixed strain variance. For spectral saturation at wavenumbers higher than 2*π*/40 m^{−1}, the two approaches conform, but they differ in their treatment of spectra with high strain variance at low wavenumbers: In the first case (Whalen et al. 2012, 2015), a maximum variance reached anywhere in the range between 2*π*/100 m^{−1} and 2*π*/40 m^{−1} would always be scaled by a GM model variance calculated between 2*π*/100 m^{−1} and 2*π*/40 m^{−1}. In the second case (Pollmann et al. 2017, and scenarios r1 of this study), the integration would stop at the particular wavenumber, for which the maximum strain variance is reached (which would be lower than 2*π*/40 m^{−1}), and the result would consequently be scaled by a GM reference lower than $\u222b2\pi /(100m)2\pi /(40m)S\xi zGM\u2061(m)\u2009dm$. This explains the strong sensitivity to the minimum spectral length observed only in this study and provides another reason for the differences between maps of *ε*_{TKE} presented here and in Pollmann et al. (2017) and those of Whalen et al. (2012, 2015). As illustrated in Fig. 11, TKE dissipation rates estimated by performing the spectral integration to at least 2*π*/40 m^{−1} while limiting the result to 0.2 are typically higher than those estimated by determining the upper wavenumber threshold by the strain variance only as done in this study—for Argo input data from 2006 to 2011, 15% of the estimates at 300–500 m depth differ by more than a factor of 3, mainly in the North Atlantic, the south Indian Ocean, the Kuroshio region, the Southern Ocean, and at the equator. While these deviations might be less pronounced for alternative implementations of the finestructure method (e.g., with fixed spectral parameters and with *n*_{min} = 1), they underline the relevance of apparently small methodological details and emphasize the importance of developing common standards and best practice guidelines for the community.

The largest sensitivity of *ε*_{TKE} to the constraint *n*_{min} = 3 is observed in areas that are associated with high mesoscale activity, such as Drake Passage or the Agulhas region. Here, strong interactions between mesoscale eddies and internal gravity waves are possible, which likely cause a deviation from the GM model that assumes a sufficient distance from the wave sources for a universal form to develop. Observations of much higher resolution than Argo profiles, combined with long-term temporal monitoring, would be necessary to better characterize the internal gravity wave spectra under such conditions, to explore how their agreement with the GM model is affected by the energy transfer from mesoscale eddies, and to assess the finestructure method with different parameter settings in locations where strong interactions with other dynamical phenomena are expected. Such an assessment has so far only been carried out at select locations and points in time (e.g., Sheen et al. 2013; Waterman et al. 2014; Takahashi and Hibiya 2019), owing to the scarcity of microstructure observations in the global ocean and especially in the Southern Ocean. Takahashi and Hibiya (2019, see, e.g., their Fig. 12) found that different formulations of the finestructure method all overestimated the results obtained from microstructure measurements at a site associated with strong mesoscale activity and linked this to high spectral energy levels at low wavenumbers. This agrees well with the above conclusions about the role of *n*_{min} and suggests that this additional constraint can under certain conditions lead to improved finestructure estimates of *ε*_{TKE}.

A comparison to high-resolution observations is also required to evaluate the slope correction term in the parameterization of internal wave energy dissipation suggested by Eden et al. (2019) [see Eq. (19)]. Applying this correction, the TKE dissipation rates calculated from spectra of minimum length *n*_{min} = 3 are increased particularly in the Southern Ocean and in the North Atlantic, while their maxima are decreased, leading in many regions to a better agreement with previous estimates. Only a comparison with observations of *ε*_{TKE} obtained from different instrumentation will show under what conditions a slope correction needs to be taken into consideration and when the finestructure method’s formulation in Eq. (16) is the more adequate approach. This formulation involves correction functions derived from theoretical and observational analysis, which extend the parameterization of *ε*_{IW} in Eq. (19) and might constitute an empirical form of slope correction.

The link between dissipation rates of TKE and internal wave energy through the Osborn relation [Eq. (20)] introduces further uncertainty by assuming a constant value for the mixing efficiency Γ, which has been shown to vary in numerical and observational studies (e.g., Gargett and Moum 1995; Peltier and Caulfield 2003; Mashayek et al. 2017). Moreover, replacing the shear production term in the Osborn relation by *ε*_{IW} assumes a direct connection between dissipated internal wave energy and dissipated TKE, which was recently questioned by Kunze (2019), who underlined the role of ageostrophic turbulence as an intermediate mechanism in the energy’s pathway to mixing. Further work will be needed to adequately incorporate these points in consistent parameterizations of wave-induced mixing in global numerical simulations.

Furthermore, the present analysis only covers part of the internal wave spectral characteristics, that is, its vertical wavenumber dependence. Following the approach of Hennon et al. (2014), it is possible to derive frequency spectra from Argo floats’ park phase information, if their temporal resolution is high enough. This is the case for floats equipped with Iridium Communications, which are capable of storing hourly temperature and salinity data. As the number of such floats increases, it will be possible to complement the present analysis by a global reference of the ocean’s internal wave field’s frequency spectrum. This could also shed light on how internal waves of different frequencies contribute to the observed variations in energy levels, spectral parameters, and wave-induced mixing, which here have to be considered representative of all frequencies between *f* and *N* owing to the lack of more detailed information.

Despite the plethora of associated uncertainties, the finestructure method is an invaluable tool in internal wave research. Seeing the scarcity of high-resolution observations on the one hand and the need for global-scale characterizations of the ocean’s internal wave field for the development of realistic parameterizations of wave-induced mixing in general circulation models on the other hand, the application of the finestructure method to the comparatively large database of lower-resolution observations can provide a solution to the dilemma—if the higher uncertainty of the derived results is recognized appropriately. Here the concept has been extended to estimate internal wave spectral properties and to give a first global-scale overview of their spatial and temporal variability based on the observations from the Argo network (Argo 2000). This provides a reference for many different situations, where so far the constant parameters of the GM model, representative of wintertime conditions at site *D* (Polzin and Lvov 2011), had to be used. This applies, for example, to the calculations of the internal waves’ group velocity in the internal wave model IDEMIX (Internal Wave Dissipation, Energy and Mixing; Olbers and Eden 2013), of the redistribution of internal wave energy through scattering at topography (Müller and Xu 1992), or of internal tide energy decay rates (Onuki and Hibiya 2018). Another prominent example is the finestructure method itself, which estimates TKE dissipation rates based on a reference value *ε*_{0} and a reference strain variance derived from a GM model spectrum with constant parameters. Particularly the variability of the internal wave energy level, which covers over two orders of magnitude in the global ocean and hence much more than the uncertainty of the finestructure method, could, when taken into consideration, lead to important modifications of the previously named calculations.

## Acknowledgments

This paper is a contribution to the Collaborative Research Centre TRR 181 ‘Energy Transfers in Atmosphere and Ocean’ funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Projektnummer 274762653. The data used in this study were collected and made freely available by the International Argo Program and the national programs that contribute to it (http://www.argo.ucsd.edu, http://argo.jcommops.org). The Argo Program is part of the Global Ocean Observing System. The EKE time series in Fig. 6 was calculated using Copernicus Climate Change Service Information. Discussions with C. Eden and D. Olbers, clarifications by C. Whalen of the details of her implementation, and comments by two anonymous reviewers are gratefully acknowledged.

### APPENDIX

#### Sensitivity to Saturation Criteria

This section is dedicated to a sensitivity test of fitted spectral parameters, internal wave energy, and TKE dissipation rates for a subset of Argo data (Atlantic Ocean, year 2012, 200–500 m) to different saturation criteria combined with different a priori defined wavenumber ranges in order to identify the most suitable one for the purpose of this study, that is, the estimation of internal wave spectral characteristics from Argo profiles. These wavenumber ranges are

r1a: 2

*π*/100 m^{−1}≤*m*≤ 2*π*/10 m^{−1}, while $\u2329\xi z2\u232a\u22640.1$;r1b: 2

*π*/100 m^{−1}≤*m*≤ 2*π*/10 m^{−1}, while $\u2329\xi z2\u232a\u22640.2$;r2: 2

*π*/100 m^{−1}≤*m*≤ 2*π*/10 m^{−1};r3:

*m*> 0 m^{−1}.

In each of these four cases, five different selection criteria are applied:

c1: a heuristic energy level criterion demanding

*E*< 0.01 m^{2}s^{−2};c2: a GM shear variance criterion demanding ⟨(

**u**_{z})^{2}⟩ ≤ 0.5*N*^{2}based on Eq. (12);c3: a shear variance criterion associating

*m*_{c}with Ri^{−1}= 1 demanding ⟨(**u**_{z})^{2}⟩ ≤*N*^{2};c4: a GM spectral saturation criterion demanding $Suz\u2061(m)<Suz,satGM\u2061(m)$ for all

*m*based on Eq. (14);c5: a spectral saturation criterion demanding $Suz\u2061(m)<Suz,sat\u2061(m)$ for all

*m*based on Eq. (15).

The comparison of the various scenarios and their effects on average spectral and integrated quantities is presented in Fig. 9. The effects are similar for the four quantities, albeit somewhat higher for spectral parameters and TKE dissipation rates (Figs. 9b,c,d) than for energy (Fig. 9a). This could be related to the fact that energy is an integrated quantity and that TKE dissipation rates are estimated by scaling the observed variance with that of the GM model, which accounts for departures from the canonical model parameters (relatively high in the scenarios in question, particularly for $m*$) differently than when no such scaling is included. The effects are also similar, when TKE dissipation rates are estimated with the correction term $h\u2061(R\omega I)$ instead of *h*(*R*_{ω}) = 1; the increase over the uncorrected average depends on the scenario, but is in any case much smaller than the method’s uncertainty (in fact, the differences are so small, that in the maps shown in Fig. 5 with their logarithmic color scale they cannot be distinguished by eye; not shown).

To identify a suitable integration range and saturation criterion, first the number of options is narrowed down by elimination. Focusing on the internal wave energy alone (Fig. 9a), only a few scenarios can be recognized as unsuited because of very high averages and standard deviations (r1b, c2; r1b, c3; r1b, c5; r3, c1). For practical reasons, the scenario r3 is added to that list owing to the low number of estimates in the cases c2–c5. That the criterion c1 does not belong to this group indicates that using an upper energy limit without additionally limiting the wavenumber range (r2) or defining a maximum shear or strain variance (r1a, r1b, r2c2, r2c3) is too weak a constraint to discard profiles in the saturated wavenumber range. Setting the maximum strain variance to $\u2329\xi z2\u232a=0.2$ (r1b) leads to unrealistically high energy levels when none of the criteria c1–c5 are applied. This suggests that specifying an upper strain variance limit of 0.2 cannot identify reliable data, which also holds true in combination with the criteria c2, c3, and c5, when too low spectral slopes are not filtered out and the bandwidth *m*_{b} becomes very large (Fig. 9d). The saturation criterion based on the observed spectral parameters, c5, is deemed unsuited for all wavenumber ranges tested here because of unrealistically high bandwidth values, which implies that spectra with too low slopes are not discarded.

This process of elimination leaves the scenarios r1ac1–c5, r1bc4, and r2c2–c4. For the former, the differences between the cases c1–c5 are very small or nonexistent, and when keeping all spectra the exact same results as for criteria c2–c4 obtained. This indicates that setting an upper strain variance limit of $\u2329\xi z2\u232a=0.1$ corresponds to applying any of the five criteria c1–c5 or is even the stricter constraint. When wavenumbers between 2*π*/100 m^{−1} and 2*π*/10 m^{−1} (r2) are considered and shear-based criteria are applied after the fitting and integration to identify and discard saturated spectra, the results and variability are lower (c2, c4), comparable to (c3), or higher (c1, c5) than in the r1a scenario, with some differences depending on the quantity in question. The fit confidence intervals are, however, typically narrower than when a maximum strain variance limit is specified (r1a, r1b), in particular for the spectral slopes.

It should be noted that the differences between the various scenarios are rarely significant when an uncertainty of more than a factor of 2 is assumed: comparing the reference scenario r2, c3 to scenarios r2, c2 and r2, c5 for Atlantic Ocean data from the years 2006–17 shows that internal wave energy levels and spectral slopes are changed on average by 10%–20% only when applying a different saturation criterion, while wavenumber scales and TKE dissipation rates are changed on average by a factor 2–3 below 300 m (a factor of 4–6 for TKE dissipation rates between 200 and 500 m). Pollmann et al. (2017) find a total uncertainty, including methodological and statistical variations, of up to a factor of 5 for finestructure estimates of TKE dissipation rates and energy levels obtained from following a comparable procedure to the one described here. The selection of a wavenumber range and saturation criterion is therefore mainly guided by pragmatic reasons: 1) criteria identified to be too weak, e.g., by allowing too large energy levels, in any of the scenarios are discarded; 2) the number of remaining estimates should not be too low; 3) since TKE dissipation rates estimated from the r1a or r1b scenario have been evaluated against microstructure observations, the mean and variance in any other scenario should roughly agree with them except for those discarded by either 1 or 2. This reasoning leaves the r1ac1–c4 (note that c2–c4 give the exact same results as if applying no criterion) and the r2c3 scenarios.

To finally choose between the remaining options, the sensitivity to the fitting parameter settings is explored (Fig. 10), focusing on the r1a scenario without any additional criterion and on the r2c3 scenario. In both cases, energy levels react similarly to changes in fitting parameters with the noteworthy exception of doubling or reducing the starting point $x0=\u2061[sGMm*GM]$, when much higher standard deviations are observed in the r1a scenario. Apart from the fivefold increase in the starting values, TKE dissipation rates are in both cases not affected by variations of *x*_{0} or upper and lower fit bounds. Spectral slopes, however, are sensitive to such variations, in particular to those of upper and lower fit bounds. When increasing the minimum number of points *n*_{min} per spectrum, the variability, number of estimates, and for high *n*_{min} also the averages, decrease, while the fit confidence bounds become narrower. The sensitivity to this parameter is higher in the r1a scenario (lower row) than in the r3, c2 scenario (upper row), in particular for spectral slopes and TKE dissipation rates. The analysis confirms that for most settings (except very high *n*_{min}), the spectral slope confidence intervals are substantially wider in the r1a scenario than in the r3c2 scenario, reaching an upper bound of 40 in the former but only 6 in the latter case. This is why the r2c3 scenario is judged to be better suited for the present purpose of calculating spectral fits.

Since shear is not actually measured but estimated from strain observations, the r2c3 scenario constitutes a reliable alternative to the r1a scenario only if the wave field’s frequency dependence can be described by that of the GM model and the polarization relations apply. The specification of a fixed strain variance maximum to identify saturation (r1a, r1b) is based on saturation concepts applied to the GM model as well as GM model shear–strain relations (e.g., Gargett et al. 1981; Munk 1981; Polzin et al. 1995; Kunze et al. 2006) and thus subject to similar constraints. Integrating over a wavenumber range up to *m*_{2} = 2*π*/10 m^{−1} (r2) rather than adjusting *m*_{2} such that the total strain variance is limited (r1a, r1b), however, could mean that the fits and integrals are carried into the saturated part of the wavenumber spectrum such that the variance is underestimated and variance criteria to discard saturated spectra do not take effect. This seems not to be the case in the r2, c3 scenario, which gives averages and standard deviations of energy, TKE dissipation rates, and spectral properties similar to the ones obtained for the more conservative strain variance criterion r1a with $\u2329\xi z2\u232a=0.1$ and is therefore considered the best option for the aim of this study.

This conclusion remains the same when considering only data from 300 to 500 m, so that the additional depth bin between 200 and 300 m is included in the sensitivity analysis of a single year’s data to increase the number of estimates, but excluded in the full analysis to reduce the uncertainty associated with near-surface variability from nonwave sources. Also for profiles with average depth levels between 500 and 1000 m (not shown), the conclusion holds true, although in addition, the scenarios r3c1 and r1bc1 emerge as suitable alternatives. When removing profiles with at least 50% of the Turner angles in the range representative of double-diffusive convection (45° < |Tu| < 90°) or considering only profiles with at least 50% of the Turner angles in the stable range (−45° < Tu < 45°) in the 200–500 m depth range (not shown), the fit confidence intervals in the r3 scenarios are notably reduced and in particular the r3, c3 scenario gives comparable results to the chosen r2, c3 scenario. This underlines that the choice of wavenumber range and saturation criterion is not biased by depth range or double-diffusive convection, but stresses that other combinations would have been equally valid.

## REFERENCES

*J. Phys. Oceanogr.*

*Annu. Rev. Mar. Sci.*

*Mar. Geod.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*ε*due to breaking of oceanic internal waves?

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Geophys. Fluid Dyn.*

*J. Geophys. Res.*

*Annu. Rev. Fluid Mech.*

*J. Geophys. Res.*

*J. Fluid Mech.*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*Geophys. Res. Lett.*

*Sci. Rep.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Fluid Mech.*

*Geophys. Res. Lett.*

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Fluid Mech.*

*J. Phys. Oceanogr.*

*Near-Boundary Processes and Their Parameterization: Proc. ‘Aha Huliko’a Winter Workshop*, Honolulu, HI, University of Hawai‘i at Mānoa, 95–105

*Rev. Geophys.*

*Geophys. Res. Lett.*

*Evolution of Physical Oceanography*, B. Warren and C. Wunsch, Eds., MIT Press, 264–291

*Deep-Sea Res. I*

*Geophys. Res. Lett.*

*J. Fluid Mech.*

*Rev. Geophys.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Ocean Dynamics*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*Annu. Rev. Fluid Mech.*

*J. Phys. Oceanogr.*

*Ocean Modell.*

*Rev. Geophys.*

*J. Phys. Oceanogr.*

*Science*

*J. Geophys. Res. Oceans*

*Nat. Climate Change*

*Deep-Sea Res.*

*Annu. Rev. Fluid Mech.*

*J. Geophys. Res. Oceans*

*J. Atmos. Sci.*

*Tellus*

*J. Geophys. Res. Oceans*

*An Introduction to Ocean Turbulence*

*J. Phys. Oceanogr.*

*Ocean Modell.*

*Geophys. Res. Lett.*

*Annu. Rev. Fluid Mech.*

*J. Phys. Oceanogr.*

*Geophys. Res. Lett.*

*J. Phys. Oceanogr.*

*Nat. Geosci.*

*J. Phys. Oceanogr.*

*Annu. Rev. Fluid Mech.*

*J. Geophys. Res. Oceans*

*Deep-Sea Res. I*

## Footnotes

^{1}

The topographic variance is estimated in each grid box—1.5° × 1.5° and 1° × 1°, respectively—using the SRTM_30PLUS bathymetry information at a resolution of 30 arc s (Becker et al. 2009) and values larger or smaller than the global mean are considered as rough or smooth, respectively, following Whalen et al. (2012).