Toward global maps of internal tide energy sinks

Internal tides power much of the observed small-scale turbulence in the ocean interior. To represent mixing induced by this turbulence in ocean climate models, the cascade of internal tide energy to dissipation scales must be understood and mapped. Here, we present a framework for estimating the geography of internal tide energy sinks. The mapping relies on the following ingredients: (i) a global observational climatology of stratification; (ii) maps of the generation of M2, S2 and K1 internal tides decomposed into vertical normal modes; (iii) simplified representations of the dissipation of low-mode internal tides due to wave-wave interactions, scattering by smallscale topography, interaction with critical slopes and shoaling; (iv) Lagrangian tracking of low-mode energy beams through observed stratification, including refraction and reflection. We thus obtain a global map of the column-integrated energy dissipation for each of the four considered dissipative processes, each of the three tidal constituents and each of the first five modes. Modes ≥6 are inferred to dissipate within the local water column at the employed half-degree horizontal resolution. Combining all processes, modes and constituents, we construct a map of the total internal tide energy dissipation, which compares well with observational inferences of internal wave energy dissipation. This result suggests that tides largely shape observed spatial contrasts of dissipation, and that the framework has potential in improving understanding and modelling of ocean mixing. However, sensitivity to poorly constrained parameters and simplifying assumptions entering the parameterized energy sinks calls for additional investigation. The attenuation of low-mode internal tides by wave-wave interactions needs particular attention.


Introduction
Ocean general circulation models (OGCMs) do not resolve the cascade of energy down to dissipation scales. They rely on parameterizations of the irreversible mixing accomplished by small-scale turbulence. This parameterized mixing is pivotal to model behaviour on climatic timescales because it is the principal source of density transformation away from boundaries (Iudicone et al., 2008). Yet ad hoc representations of irreversible mixing, such as the specification of diffusivities that are fixed in time, remain common in ocean modelling. Tuning of fixed diffusivities may allow fair reproduction of some targeted observations, but their prescription generally precludes a realistic model response to perturbations. Indeed, for an ocean model to appropriately respond to changes in external forcing, irreversible mixing must be connected to the forcing via conservative and realistic tracking of energy .
One important forcing to account for is the gravitational attraction of the Moon and Sun, which generates surface (barotropic) ocean tides. About two-thirds of the power input to surface tides is directly lost to shear-driven mixing and frictional heating at continental margins (Taylor, 1919;Jeffreys, 1920). The remaining third-about 1 TW-feeds the generation of internal (baroclinic) tides (Egbert and Ray, 2000;Nycander, 2005). Internal tides are internal waves of tidal frequency that radiate from sloping topography where barotropic tidal currents push stratified fluid up and down the slope (Munk, 1966;Bell, 1975). They are a major energy source for the internal wave field and, via instabilities of that field, for irreversible mixing in the ocean interior (Sjöberg and Stigebrandt, 1992;Munk and Wunsch, 1998;Rudnick et al., 2003).
Internal tides of small vertical scale are more prone to instability than larger-scale waves, which can propagate over longer distances and dissipate far from their generation site (Olbers, 1976 Garrett, 2002;Garrett and Kunze, 2007). To characterize the vertical scales of internal tides, it is customary to describe their vertical wavenumber spectrum by a discrete set of vertical normal modes or "equivalent" modes. The first few modes have been observed to travel up to several thousands of kilometres (Dushaw et al., 1995;Ray and Mitchum, 1996;Zhao et al., 2016). In contrast, higher modes are thought to break into small-scale turbulence within a relatively small radius of their emission. This premise, together with available maps of internal tide generation (Jayne and St Laurent, 2001;Carrère and Lyard, 2003;Nycander, 2005), has been used to map the near-field dissipation of internal tides and parameterize the associated mixing in OGCMs: it is commonly assumed that a fixed fraction of the power input to internal tides feeds local, bottom-intensified turbulence Simmons et al., 2004). This power fraction-related to the modal content of internal tide generation-is traditionally taken as onethird, based on turbulence observations from the Brazil Basin , for lack of a global map of its spatial variations (Falahat et al., 2014a).
Most models do not explicitly represent the contribution of far-field internal tide dissipation to mixing. Rather, a fixed diffusivity of about 10 −5 m 2 s −1 is commonly prescribed as representative of "background" internal wave-driven mixing (e.g., Jochum, 2009). The power effectively consumed by this mixing is then not controlled nor tied to the supply of energy by external forcings. In order to resolve this disconnect, substantial effort has been deployed to refine our understanding of the generation, propagation and dissipation of low-mode internal tides (MacKinnon et al., 2017). An estimate of the modal distribution of internal tide generation has been attained (Falahat et al., 2014b). Further estimation of the direction of emission of energy beams is ongoing (Pollmann et al., 2019). The propagation through variable stratification and geostrophic flows has been studied through a combination of in-situ measurements, satellite altimetry and high-resolution numerical experiments, and found to agree well with linear wave theory (Rainville and Pinkel, 2006;Alford and Zhao, 2007b;Arbic et al., 2010;Zhao et al., 2010;Ray and Zaron, 2016).
Wave-wave interactions have long been recognized as a conduit to wave breaking via downscale energy transfer (Olbers, 1976;McComas and Bretherton, 1977;Müller et al., 1986;Henyey et al., 1986). They are dominated by triadic wave instabilities, also known as parametric subharmonic instability (PSI), which are active equatorward of about 30°(15°) for the semidiurnal (diurnal) internal tides ). PSI appears to be a significant sink of low-mode energy but seemingly not the dominant dissipation pathway of the semidiurnal mode-1 tide (Hibiya et al., 1998;Alford et al., 2007;Hazewinkel and Winters, 2011;Xie et al., 2011;MacKinnon et al., 2013b;Sun and Pinkel, 2013). Much attention has therefore been directed to the interaction of low-mode waves with topography. Rays bouncing off smallscale seafloor roughness lose energy through scattering to higher modes (Müller and Xu, 1992;St Laurent and Garrett, 2002;Bülher and Holmes-Cerfon, 2011;Kelly et al., 2013). Beams reaching larger-scale topographic obstacles can either reflect back to deeper waters, dissipate at the slope, shoal, or undergo a combination of the latter: their fate depends in a largely predictable way on the height and shape of the obstacle (Müller and Liu, 2000;Johnston and Merrifield, 2003;Nash et al., 2004;Klymak et al., 2011Klymak et al., , 2013Mathur et al., 2014;Legg, 2014). Other dissipative processes, such as wave capture or scattering by mesoscale eddies (Bühler and McIntyre, 2005;Polzin, 2008;Dunphy and Lamb, 2014;Dunphy et al., 2017), remain less well documented or quantified.
Building on this progress, steps toward comprehensive parameterization of tidal mixing have been taken. Niwa and Hibiya (2011) constructed a horizontal map of low-mode internal tide dissipation using high-resolution OGCM experiments that include tidal forcing and a 30-day linear damping of baroclinic wave fluctuations. Such experiments allow investigation of the generation and propagation of the first few modes, but they do not resolve the full wave spectrum and the inferred dissipation is dependent on model numerics and adjustable damping terms (e.g., Buijsman et al., 2016). An alternative approach consists of calculating the horizontal propagation of depth-integrated low-mode energy through an observed or modelled stratification, using linear wave theory and parameterized energy sinks . This approach allows estimation of energy sinks due to specific processes, known to produce specific vertical structures of dissipation (e.g., Melet et al., 2016), thus paving the way for a comprehensive, three-dimensional mapping of internal tide-driven mixing. Indeed, the vertical distribution of far-field dissipation is a key ingredient of internal wave-driven mixing parameterizations (de Lavergne et al., 2016a(de Lavergne et al., , 2016bMelet et al., 2016).
Here we adapt and extend the framework developed by  to construct climatological two-dimensional maps of lowmode dissipation, broken down into four dissipative processes, for the first five vertical modes of the three main tidal constituents. We then sum the estimated low-mode dissipation with the local dissipation of higher modes to obtain a horizontal map of the total internal tide dissipation. This map shows fair agreement with full-depth dissipation estimates from observed finescale strain (Kunze, 2017a), suggesting that the framework has potential to narrow the gap toward comprehensive mapping and parameterization of tidal mixing.
The presentation proceeds as follows. The methodology, including the estimation of energy sources, energy sinks and beam trajectories, is detailed in section 2. The resultant distributions of internal tide energy loss and their sensitivity to least constrained parameters of the calculation are presented in section 3. Next, we compare our results to independent observational estimates of internal wave energy dissipation (section 4) and to previous global budgets of internal tide energy sinks (section 5). Limitations and perspectives are discussed in the concluding section. Details about the treatment of interactions with topographic slopes are provided in Appendix A. The rationale for choosing a Lagrangian-rather than Eulerian-propagation scheme is exposed in Appendix B.

Strategy
Consider internal waves of given tidal frequency and vertical mode number. Their column-integrated energy E is a function of time t, geographical position = r x y ( , ), and angle ϕ of the horizontal wavevector k h . The energy evolution is governed by ) is the horizontal energy transport by modal group velocity c g , G is the position-and angle-dependent generation rate of the considered mode, and D encapsulates energy sinks. Note that divergence of F includes both position (propagation) and angle (refraction or reflection) terms. Since c g is parallel to k h , ϕ represents the direction of energy propagation in the horizontal.
Our objective is to map the steady-state dissipation D(x, y). Because propagation and dissipation depend on the frequency and mode number of the internal tide, we separately consider each mode of the three most energetic tidal constituents-M 2 , S 2 and K 1 (Egbert and Ray, 2003). We exploit the mode-partitioned generation estimates of Falahat et al. (2014b) (section 2.2), propose simple formulations for the decay due to wave-wave interactions (section 2.3.1) and scattering off abyssal hills (section 2.3.2), and introduce a geometric treatment of interactions with large-scale topographic slopes (section 2.3.3). Propagation and refraction are modelled using linear wave theory applied to the annual mean WOCE global hydrographic climatology (Gouretski and Koltermann, 2004). We thus ignore time variations and build the steady-state dissipation field by tracking individual energy beams from sources to sinks, then summing over all beams (section 2.4).
This strategy resembles that of , who solved Eq. (1) for the first-mode M 2 internal tide using the WOCE climatology, a map of internal tide generation rates and parameterizations for energy transfers to high modes. The present approach differs from theirs in the following principal ways: (i) energy beams are tracked one by one using a Lagrangian scheme, rather than (1) being solved using Eulerian advection schemes; (ii) the first five vertical modes are treated individually; (iii) interactions with large-scale topographic slopes are explicitly modelled; (iv) parameterization choices for wave-wave interactions and topographic scattering differ. Falahat et al. (2014b) developed a semi-analytical model for the tidal generation of internal waves projected onto vertical normal modes. Using global observational climatologies of stratification (Gouretski and Koltermann, 2004), bathymetry (Smith and Sandwell, 1997) and barotropic tidal velocities (Egbert and Erofeeva, 2002), they mapped the generation rate of each of the first 10 modes of the semidiurnal M 2 tide. For the purpose of this study, we calculated analogous maps for the semidiurnal S 2 and diurnal K 1 tides. The global total power input to M 2 , S 2 and K 1 internal tides is 634, 148 and 68 GW, respectively. 1 These integral values exclude seafloor depths < 400 m, where the calculation becomes unreliable due to violation of the assumptions of small tidal excursion and subcritical topography (Falahat et al., 2014b). They are only slightly lower than equivalent values derived from the same products and the non-modal formulation of Nycander (2005): 680, 150 and 83 GW for M 2 , S 2 and K 1 , respectively.

Energy sources
The modal partitioning of global generation rates is similar for the two semidiurnal constituents (Table 1): modes 1, 2 and 3 receive about 35%, 25% and 15% of the total energy flux, respectively. Together, modes ≥6 represent only about 10% of the total power. However, the latter fraction is a lower bound since it does not incorporate modes ≥11, which are not resolved by the calculation. In particular, abyssal hill topography, which is absent from the 1/30°-resolution 'etopo2v2' bathymetry product (Smith and Sandwell, 1997), is responsible for significant high-mode (≳50) internal tide generation (Melet et al., 2013;Lefauve et al., 2015). Adding the M 2 internal tide generation by abyssal hills as estimated by Melet et al. (2013), the power received by modes ≥6 rises to 23% of the M 2 total (Table 1). Generation over abyssal hills has not been estimated for the other tidal constituents. Compared to M 2 and S 2 , K 1 internal tide generation projects less strongly onto modes 1-2, which receive no more than 40% of the total flux (Table 1). This different modal distribution relates to the different spatial distribution of diurnal conversion, which is dominated by hotspots in Luzon Strait and Indonesian seas (not shown).
Total (modes 1-10) generation rates integrated over the Pacific, Atlantic and Indian basins are shown as a function of seafloor depth in Fig. 1. For all three constituents and all three basins, generation rates decrease with depth. The decrease is sharper in the Pacific, which hosts 62%, 58% and 82% of the M 2 , S 2 and K 1 barotropic-to-baroclinic conversion, respectively. Profiles of conversion predicted by the non-modal formulation of Nycander (2005) are very similar (compare thin and thick curves in Fig. 1). However, spatial maps of conversion reveal important differences: as opposed to the non-modal estimate (Fig. 2c), the calculation of Falahat et al. (2014b) often predicts negative conversion rates (Fig. 2a). Negative rates could be interpreted as baroclinic-to-barotropic conversion. However, it is not known whether such conversion actually occurs, and it is unclear whether these apparent internal tide energy sinks can feasibly connect with the sources.
To circumvent the issue, we apply a correction to each modal conversion field to remove negative values while preserving key integral properties. Specifically, we set all negative values to zero, and multiply the resulting field by a depth-dependent correction factor to recover the original depth-dependent generation rate integrated over each of the Pacific, Atlantic and Indian oceans. This procedure is performed on each field at the original 1/30°-resolution using a vertical bin size of 200 m. We finally average generation rates onto the 0.5°-resolution grid of the WOCE hydrographic climatology. The resulting conversion rates of the M 2 tide are shown in Fig. 2. Summed over modes 1 to 10, the total corrected field (Fig. 2b) closely resembles its nonmodal counterpart (Fig. 2c): grid cells where the two fields agree within a factor of 2 cumulate 87% of the total power input of 634 GW. We expect the slightly smoother patterns of the corrected field to be more realistic, reflecting superior accuracy of the modal formulation demonstrated in test cases (Falahat et al., 2014b).
Corrected generation maps for each mode also display realistic patterns ( Fig. 2d-g). The western Pacific (including Indonesian seas) stands out as the main source region for the first few modes. The western Indian Ocean and Mid-Atlantic Ridge are the next two most prominent generation regions. As mode number increases, power input gradually spreads from localized hotspots at abrupt topography (such as Hawaii) to more distributed sources above flatter but corrugated topography (such as ridge flanks), in accord with theory and observations (St Laurent and Nash, 2004;Falahat et al., 2014a). The estimated power input to modes ≥6 ( Fig. 2g) is dominated by the contribution of abyssal hills, and therefore by regions centred around slow-spreading ridges (Goff, 2010). Similar trends are found for the conversion of S 2 and K 1 tides (not shown). In summary, the applied correction for negative conversion rates introduces important uncertainties but produces plausible fields that are suitable for the present mode-by-mode global mapping exercise.
The generation of internal waves by the principal eight tidal constituents has been estimated by Nycander (2005). Excluding regions shallower than 400 m, M 2 accounts for 67% of the global energy flux summed over all constituents; M 2 , S 2 and K 1 together account for 90%. Tidal constituents with proximate frequencies cause internal wave generation rates that differ in magnitude but little in distribution (Egbert and Ray, 2003). Using M 2 , S 2 and K 1 as our best proxy for N 2 , K 2 and the next three diurnal constituents, respectively, we estimate the mode-by-mode generation across all constituents as where appropriate power ratios have been introduced. Including the contribution of abyssal hills, we thus calculate that mode 1 takes up about 30% of the total internal tide generation ( Table 1). The Table 1 Global power input to internal tides per mode and per constituent. Percentage and total power values derive from the semi-analytical model of Falahat et al. (2014b); they exclude seafloor depths < 400 m and include positive and negative local rates. Values in bold incorporate the contribution of abyssal hill topography to high-mode M 2 internal tide generation, as estimated by Melet et al. (2013). 'All' incorporates the eight most energetic tidal constituents following Eq. (2). percentage decreases monotonically with mode number, reducing to 5% at mode 5. These global percentages mask pronounced spatial contrasts ( Fig. 3): internal tides generated at continental slopes are dominantly mode 1, whereas those radiated off rough and less steep topography are dominantly mode 2 and higher. Modes ≥6 dominate over the global mid-ocean ridge system. The angular distribution of energy sources has not yet been mapped. As a surrogate for such a map, we use bathymetry information to direct energy sources away from large-scale topographic slopes. We fit planes to the 'etopo2v2' bathymetry over half-degree grid squares, find the vector normal to each of these planes, obtain the angle ϕ g of that vector in the horizontal (with respect to due east), and distribute energy sources around ϕ g with a cosine weight: x y G x y x y ( , , ) 1 2 ( , ) max ( 0, cos ( ( , )) ) .
This distribution ensures that most of the energy radiates away from bathymetric slopes, with no energy directed upslope. Most internal tide sources are focused in one direction . However, propagation through the time-variable, eddying ocean tends to rapidly disperse time-average energy fluxes (Rainville and Pinkel, 2006;Alford and Zhao, 2007a;Zaron and Egbert, 2014;Vic et al., 2018), motivating the choice (3) of angularly distributed sources. We examined sensitivity to this choice by running an additional experiment, referred to as BEAM, with In the numerical code, Eqs. (3) and (4) are discretised at the angular resolution of π/30. Eq. (4) means that local energy sources are placed within a single angle bin, that which contains ϕ g (x, y).

Wave-wave interactions
Wave-wave interactions constitute a prominent but poorly quantified dissipation pathway of low-mode internal tides. Satellite altimetry allows tracking of low-mode energy beams that remain coherent with the astronomical tide (Ray and Mitchum, 1996;Zhao et al., 2016). However, loss of coherence means that satellite observations provide only an upper-bound on the open-ocean dissipation rate of low-mode beams (Rainville and Pinkel, 2006;Buijsman et al., 2017;Zaron, 2017). Other constraints are called for. Theory indicates that the decay rate due to wave-wave interactions increases with mode number and is elevated equatorward of the "PSI latitude" where the tidal frequency ω is twice the Coriolis frequency f (McComas and Müller, 1981;Hibiya et al., 2002;. Observations of a beam travelling northeastward from Hawaii suggest that PSI barely attenuates the firstmode semidiurnal internal tide in this region (Alford et al., 2007;Zhao et al., 2010;MacKinnon et al., 2013b). This is supported by dedicated numerical experiments (Hazewinkel and Winters, 2011), in which the energy loss of the mode-1 M 2 tide was inferred to be 15% over 760 km. Using an appropriate local group speed of 2.5 m s −1 , the latter estimate translates into an e-folding decay time of 21 days. We choose 20 days as a representative timescale over which mode-1 tides are damped by PSI-noting that, in reality, this timescale varies in space and time and decreases with the amplitude of the participating waves (MacKinnon et al., 2013b;Sun and Pinkel, 2013;Ansong et al., 2018).
Attenuation of low modes by other wave-wave interactions is thought to be an order of magnitude slower (Hibiya et al., 1998;. In-situ observations northeast of Hawaii suggest a decrease of diffusivity by a factor of ∼4 north of the M 2 PSI latitude (28.8°), occurring over four latitude degrees, with little concomitant change in stratification (MacKinnon et al., 2013a). We take the efolding decay time of mode 1 to increase fourfold over the same latitude span, plateauing at 80 days poleward of the transition (Fig. 4). For lack of independent constraints, an analogous latitudinal variation is employed for S 2 and K 1 mode-1 tides, taking into account their different PSI latitudes (29.9°and 14.5°, respectively). Note that linear free waves do not exist poleward of the turning latitude at which their frequency equals the inertial frequency: M 2 , S 2 and K 1 horizontal energy fluxes are thus confined equatorward of 74.5°, 85.8°and 30.0°, respectively.
We must also specify the dependence on mode of these attenuation rates. Evaluation by Olbers (1983) of the transfer integral within a Garrett-Munk wave spectrum suggests that decay times scale with the −2 power of mode number (cf. his Fig. 33) This approximate dependence might not be representative of all types of wave-wave interactions nor of all types of wave fields . Lacking better knowledge of the relative efficiency of wave-wave interactions as a function of mode, we set the e-folding decay time of mode number n as τ wwi (n) = τ wwi (1)/n 2 . Since the group speed is inversely proportional to mode number, this choice means that wave-wave interactions cut horizontal energy fluxes F by a factor 1/e over a length c g τ wwi proportional to 1/n 3 . This decay length is shown in Fig. 5 for the first four modes of the M 2 tide. The group speed has been computed from the annual mean stratification of the WOCE climatology as (Llewellyn Smith and Young, 2002) where N is the depth-mean buoyancy frequency and H the ocean thickness. Shortest decay lengths occur within regions of weak depthintegrated stratification, such as shelves and polar seas, where slow  Fig. 1. Conversion rate into modes 1-10 as a function of seafloor depth, integrated over the (dark blue) Pacific, (pale blue) Atlantic and (orange) Indian oceans, including both positive and negative local rates, for each of (a) M 2 , (b) S 2 and (c) K 1 tidal constituents. Thick curves are estimates of Falahat et al. (2014b); thin curves are estimates of Nycander (2005). The three ocean masks have been extended south to Antarctica as in Falahat et al. (2014b); the Arctic is included in the Atlantic. group velocities leave more time for wave-wave interactions to act (Fig. 5). Away from these regions, the decay length of the first-mode M 2 tide is typically a few 1000 km in the latitude range where PSI is active, and 10,000 km or more at higher latitudes. By contrast, the decay length of mode 2 is generally < 1000 km; that of mode 4 rarely exceeds 100 km.
The chosen latitude-and mode-dependence of the e-folding time over which wave-wave interactions erode low-mode beams is broadly compatible with observations but remains largely a place-holder awaiting refined representation of this energy cascade. More complex and more realistic parameterizations of the rate at which low-mode internal tides lose energy to the background wave field can be considered in future implementations. As a first step, in order to gauge the sensitivity of the calculation to the chosen timescales, we performed two sensitivity experiments, referred to as FAST and SLOW, in which efolding decay times are uniformly halved and doubled, respectively.

Scattering by abyssal hills
Bühler and Holmes-Cerfon (2011) studied the decay of a low-mode internal tide as it bounces over random bottom topography. They showed that the energy flux decays as exp (−λ sca r), where the decay rate per unit horizontal distance r is a function of topography statistics only: Topographic corrugations h(r) enter Eq. (6) through the probabilistic expectation of their squared height and of their squared spatial derivative. Here we apply Eq. (6) to abyssal hill statistics mapped by Goff (2010). The decay rate can be approximated in terms of the local rms height h rms and local mean horizontal wavenumber κ of abyssal hills as The corresponding e-folding decay length λ sca −1 is shown in Fig. 6.
Scattering is efficient near ridge crests, where the decay length is typically several 100 km, but is weak over old or sediment-rich seafloor. The predicted decay rates do not depend on the frequency and mode number of the low-mode internal tide, and are a factor of 3 to 10 smaller than those estimated by . By contrast, they are roughly one order of magnitude stronger than deep-ocean scattering rates of the first-mode M 2 tide calculated by Kelly et al. (2013) using the model of Nycander (2005) and bathymetric variations resolved in 'etopo2v2'. This would suggest that scattering by resolvedscale topographic roughness in the abyss is negligible. The role of topographic slopes resolved in current global bathymetric maps is considered next.

Interactions with topographic slopes
The fate of a low-mode internal tide impinging onto a large-scale topographic obstacle depends on the ratio of the topographic slope to the ray-path slope (Fig. 7). The ray-path slope or wave slope s follows from the dispersion relation as independent of mode number. Topographic slopes smaller than s, termed subcritical, allow the internal tide to bounce forward and shoal (e.g., Müller and Liu, 2000;Legg, 2014). Those approximately equal to s, termed critical, cause a transfer of energy to very high modes and boundary turbulence (e.g., Eriksen, 1982;Legg and Adcroft, 2003;Martini et al., 2013). Topographic slopes significantly steeper than the wave slope, termed supercritical, mainly result in the backward reflection of the internal tide (e.g., Klymak et al., 2013;Johnston et al., 2015). These general principles, allied with knowledge of the ocean's topography and stratification, can be used to estimate a priori the Percentage generation rate transmitted, reflected and dissipated fractions of incident energy fluxes. The rationale is illustrated in Fig. 7 in the simplified situation of a threesegment obstacle of height h immersed in constant stratification. Consider a horizontal energy beam or energy flux F b encountering shoaling topography. This beam is made of individual rays travelling up and down the water column with slope s, meeting the surface every 2H/s of travelled horizontal distance. We assume rays to be randomly distributed, a reasonable assumption given decoherence of internal tides by the fluctuating ocean environment. The fraction of F b dissipated at critical topography is first estimated as the proportion of rays encountering a near-critical slope, i.e. as the ratio of the horizontal length shown in red to the distance between consecutive bounces 2H/s (Legg, 2014). We define near-critical slopes as those ranging between 0.8 and 1.5 s, the approximate range over which critical reflection dominates (Muller and Liu, 2000;Mathur et al., 2014;Legg, 2014). Next, we calculate the reflected fraction as the portion of rays ultimately impinging on supercritical topography (slopes exceeding 1.5 s). The remaining fraction of the energy flux is transmitted. Transmitted internal tides undergo shoaling, which may cause their direct breaking via increase of their energy density and shear (Legg, 2014). We further parameterize the fractional dissipation of the transmitted flux due to shoaling as (h/H) 2 , a simplifying approximation based on the model results of Legg (2014;cf. her Figs. 6 and 14).
As detailed in Appendix A, this geometric treatment of interactions with bathymetric slopes can be generalized to depth-varying stratification and arbitrary topography. Dissipated, reflected and transmitted fractions were estimated within each grid column of the 0.5°-resolution WOCE climatology, for each of the four cardinal directions and each of the three tidal constituents, using depth-varying wave slopes and 1/30°resolution bathymetric variations. Dissipation and reflection are zero in downslope directions, except at crests, whose specific treatment is outlined in Appendix A. With non-uniform stratification, the fractional dissipation due to shoaling of the transmitted energy flux becomes the square of the WKB-scaled bathymetric rise, where z is depth, H is the local ocean thickness and h is the thickness of a shallower neighbouring grid column. This simple parameterization neglects the dependence on Froude number of shoaling-driven wave 10,000 km 1,000 km 100 km 10,000 km 1,000 km 100 km breaking (Legg, 2014), likely underestimating (overestimating) fractional dissipation of most (least) energetic transmitted waves. A low bias of shoaling-induced dissipation could be implied on average, though uncertainties make a more detailed assessment premature. Fig. 8 shows the estimated fractions for M 2 . There are four values (for four directions) at each grid point: we show the maximum of the four, usually corresponding to the direction along which topography shoals most rapidly. Energy loss due to shoaling is substantial at the shelf break and shoreward, but quasi-negligible elsewhere. Critical and supercritical topography primarily consists of continental slopes. In accord with previous studies (e.g., Kelly et al., 2013;Johnston et al., 2015;Klymak et al., 2016), supercritical reflection generally dominates, implying that much of the internal tide energy reaching continental slopes turns back to the open ocean. Nonetheless, the estimated fractional dissipation at critical slopes is significant around most of the open ocean's perimeter, and is non-negligible near ridge crests and major fractures or seamounts. The S 2 tidal constituent has a frequency close to that of M 2 , and therefore has similar ray paths and similar transmission/reflection/dissipation characteristics as M 2 . In contrast, the lower frequency of K 1 internal tides means that their ray paths are less steep, implying fewer but more often supercritical encounters with the bottom. As a result, K 1 energy reflection at continental slopes is proportionately large (not shown).

Energy propagation
To quantify and map the modelled energy sinks, we propagate energy fluxes from generation sites through climatological stratification, accounting for refraction and for reflection at supercritical slopes. We adopt a Lagrangian approach: beams are tracked one by one throughout their lifecycle. The procedure for a given beam is the following. The energy flux F b carried by the beam emerges from the centre of the generation grid cell (x, y), at an angle ϕ with respect to the east, with magnitude G(x, y, ϕ). We first determine the exact position it reaches along the grid-cell perimeter, and save the energy lost to wave-wave interactions within the generation grid cell. The beam tracing algorithm is then identical from grid cell to grid cell: 1. Given the angle of F b , we find the next grid-cell edge and position along that edge reached by F b , and determine crossed distance Δr. 2. The dissipation by wave-wave interactions and by scattering off abyssal hills within the traversed grid cell are saved. F b is reduced accordingly. 3. The angle of F b is updated due to refraction. 4. The magnitude and emission angle of reflected energy flux are saved. The dissipation due to critical slopes is saved. F b is updated. 5. The dissipation due to shoaling is saved. F b is updated.  Fig. 7. Schematic illustrating the geometric treatment of interactions with topographic slopes. To simplify the illustration, stratification (N) is taken to be constant, so that ray-path slopes (s) are also constant. Rays that encounter subcritical slopes reflect forward; they either shoal or bounce against a supercritical slope to eventually reflect backward. Rays that encounter a critical slope dissipate along the slope. Rays that directly encounter a supercritical slope reflect backward. The fraction of the incident energy flux that shoals, dissipates or reflects backward is given by the proportion of rays undergoing each scenario: it is the ratio of the yellow, red or blue horizontal length to the distance between consecutive bounces (2H/s). Ratios are evaluated within each half-degree grid square, for each tidal constituent and each of four directions, using the 1/30°-resolution bathymetry and the climatological stratification (see Appendix A for details).
Interactions with topographic slopes are performed last because they apply to bathymetric variations from the traversed grid cell to the next. The emission angle of reflected fluxes is determined by specular reflection against plane-fits to the 1/30 c -resolution bathymetry. The remainder of slope interactions do not depend on the exact incident angle but only on the appropriate cardinal direction-a practical but reasonable simplification see Appendix A). Specular reflection at the land-ocean boundary ensures no energy is lost. The latter reflected fluxes are found to be very small. The sequence above is applied to all source triplets (x, y, ϕ) characterised by G(x, y, ϕ) > 0. Next, the source term G is replaced by the saved field of reflected fluxes and the whole procedure repeated. Five such iterations suffice to reach a total dissipation within 1% of the total generation. We perform the calculation for each of M 2 , S 2 and K 1 and each of the first five modes; we find that barely any mode-6 internal tide energy escapes generation grid cells given the modelled attenuation rates (see propagation distances in Table 2). Using a π/30 angular resolution for G, the calculation for a single mode takes between a few minutes (mode 5) and a few hours (mode 1) to complete on a laptop computer. This makes the framework well-suited for sensitivity studies.
By treating modes 1 to 5 independently, we have implicitly assumed that there is no exchange of energy between these low modes. In other words, we have assumed that: (i) dissipative processes transfer energy directly to locally-dissipating modes or to small-scale turbulence, and (ii) propagation and supercritical reflection do not alter the modal content of low-mode beams. Assumption (i) is justified for the modelled interactions with topography, thought to generate small-scale motions (Bühler and Holmes-Cerfon, 2011;Nash et al., 2004;Legg, 2014). PSI also tends to generate relatively high vertical wavenumbers (Hibiya et al., 1998(Hibiya et al., , 2002. Nonetheless, wave-wave interactions-or interactions with vortical flows, not modelled here-may cascade energy from one low mode to another (McComas and Bretherton, 1977;Dunphy and Lamb, 2014). Such transfers presently lack quantification and are thus not considered here. Assumption (ii) is unlikely to hold everywhere. Open-ocean propagation can be assumed to preserve mode number (Rainville and Pinkel, 2006). However, reflection at supercritical slopes generally alters the modal content of energy fluxes (Müller and Liu, 2000;Klymak et al., 2016). To explore the impact of wavenumber modification upon reflection, we carried an additional experiment, termed EXCH, where only half the magnitude of reflected fluxes conserves the mode number of incident fluxes. The other half is equally distributed among higher propagating modes (n ≤ 5): it is saved as an additional energy source for these modes. Experiment EXCH thus demonstrates the feasibility of including a downscale energy transfer between low modes.

Energy transport and content
We begin by illustrating the energy redistribution accomplished by the first-mode M 2 internal tide. Example beam trajectories, originating from three source grid cells with three different emission angles each, are depicted in Fig. 9. Relatively slow open-ocean attenuation, particularly poleward of 30°, allows the beams to travel over basin scales. Most beams lose significant energy at continental slopes and give birth to one or several reflected beams. Refraction acts to bend trajectories (i) toward the equator and (ii) toward weaker depth-integrated stratification. Effect (ii) steers beams into shallower regions, increasing the total mode-1 dissipation due to critical reflection and shoaling by a few percent (not shown). Effect (i) diminishes the energy transport to high latitudes. For example, if refraction is unaccounted for, the predicted mode-1 dissipation south of 60°S increases by 30%. Horizontal energy fluxes accumulated across all beams, shown in Fig. 10a,b, highlight prominent regions of first-mode semidiurnal internal tide activity. Hotspots include waters around Madagascar, Indonesia, New Zealand, French Polynesia, Hawaii, east of China and Japan, and west of Morocco and Europe.
To assess the realism of our framework, we compare M 2 mode-1 energy transports of the REF and BEAM experiments (Fig. 10a,b) to those extracted from satellite observations of sea surface height ( Fig. 10e; Zhao et al., 2016). All hotspots listed above also stand out in the satellite-derived map. Conversely, only one observed hotspot is not captured here: the continental slope of northeastern Brazil. Colocated hotspots indicate that major sites of M 2 mode-1 generation are accurately placed in the corrected conversion field (Fig. 2d). However, transports detected in satellite observations are much weaker than their  Kelly et al., 2013;Zhao et al., 2016), so that substantial overestimation of energy sources is unlikely. Alternative explanations for the discrepancy include: (i) limited detection of internal tides propagating perpendicularly to satellite tracks; (ii) absence of incoherent internal tide energy fluxes in the satellite map. Because satellite tracks are generally oriented in the north-south direction, the high-pass filter applied to sea surface height data tends to suppress east-west oriented waves Zhao, 2018). This suppression is conspicuous in the aggregate angular distribution of energy transports, largely north-south in the satellite product (Fig. 10f,  red), but largely isotropic and dominantly eastbound in the REF and BEAM experiments (Fig. 10f, blue and black). To mimic the suppression of eastbound/westbound waves, we multiply local angle-resolved transports of REF and BEAM by the red curve in Fig. 10f, that is, by the aggregate angular distribution of satellite-derived transports normalized by its maximum. The resulting maps of energy transports (Fig. 10c,d) are in better agreement with the observational map (Fig. 10e), both in terms of magnitude and spatial patterns. In particular, the global energy content in REF (BEAM) decreases from 165 (172) to 64 (69) PJ, a level more comparable to the 36 PJ reported by Zhao et al. (2016). Masking out regions unmapped in Zhao et al. (2016) further reduces by 10 PJ the energy content in REF and BEAM. The remaining mismatch likely reflect loss of coherence, particularly in the equatorial band (Shriver et al., 2014;Buijsman et al., 2017;Zaron, 2017), where satellite-derived transports are much weaker than those modelled here.
Comparison with satellite-observed internal tides also illustrates biases in the angular distribution of sources in REF and BEAM. The REF energy emission is too diffuse, whereas the BEAM experiment underestimates angular spreading at and away from sources. The source  . Panel (f) shows the angle distribution of globally integrated energy transports in (black) REF, (blue) BEAM and (red) satellite observations. Each angle distribution is normalized by its maximum; the dashed circle marks a normalized value of 1. In (c) and (d), angle-resolved transports have been multiplied by the red curve in (f) prior to summation over all angles; this selects the northbound/southbound orientations that dominate in the satellite product. Energy transport is the sum over all beams of F r b , per unit area. Energy content is the same quantity divided by the group velocity. The globally integrated energy content in PetaJoules (PJ) is given for each of (a-e). orientation ϕ g deduced from topography appears to be a decent approximation at some but not all prominent generation sites. Future work should strive to constrain the angle-dependence of generation and to model the dispersion of energy beams by fluctuating currents and stratification. Here only a preliminary assessment of the sensitivity of dissipation to angular spreading will be performed, by contrasting REF and BEAM experiments. We henceforth focus on the REF experiment; sensitivity will be examined in section 3.3. Table 2 summarizes global characteristics of M 2 low-mode energy propagation. On average, energy in the first mode travels 1440 km, at a speed of 2.3 m s −1 . The corresponding mean residence time is 9 days. Because the e-folding length of attenuation by wave-wave interactions is proportional to 1/n 3 , travel distance and energy content decrease sharply with mode number. The distribution of the energy content among modes 1 to 5 resembles local estimates from mooring observations (Zhao et al., 2010;Vic et al., 2018). Cumulated over the first five modes, the global energy content of M 2 internal tides approaches 300 PJ, in the range of extrapolations from in situ observations (Alford and Zhao, 2007a;Kelly et al., 2013). Mode-5 internal tide energy travels only 68 km on average, so that its dissipation is near local.

Energy dissipation
The contrasting propagation characteristics of modes 1 to 5 translate into contrasting distributions of dissipation (Fig. 11). The first mode is the only one capable of redistributing energy on planetary scales, and is therefore the main contributor to dissipation in regions remote from important generation sites, e.g., at southern high latitudes. Despite this redistribution, the dissipation of the first-mode M 2 tide remains focused around major source regions-namely the western Indian, western Pacific and mid-Atlantic-and equatorward of 30° (  Fig. 11a). We estimate that wave-wave interactions and scattering off abyssal hills in the open ocean contribute 45% of the dissipation of the mode-1 M 2 tide, leaving substantial energy loss at continental slopes and shelves (Table 2). Dissipation becomes increasingly clustered around energy sources (Fig. 11b-f) and increasingly dominated by the contribution of wave-wave interactions (Table 2) as mode number increases. Whereas interactions with topography account for the bulk (72%) of the dissipation of the first mode, this balance is reversed at higher modes: wave-wave interactions monopolise 74%, 90%, 96% and 99% of the dissipation of mode 2, mode 3, mode 4 and mode 5, respectively. A small fraction of the energy loss by wave-wave interactions occurs on shelves (defined here by H < 400 m), so that the percentages of dissipation by open-ocean wave-wave interactions given in Table 2 are slightly lower.
Dissipation summed over modes 1 to 5 is shown in Fig. 12 for each tidal constituent. The spatial pattern of energy loss by low-mode S 2 internal tides is very similar to that of M 2 low modes. However, S 2 dissipation extends further into polar oceans and is about four times weaker. Diurnal K 1 internal tides are less energetic still, and their dissipation is focused in the tropical western Pacific (Fig. 12c). In particular, there is substantial energy loss near the two dominant K 1 generation sites, Luzon Strait and the Indonesian archipelago, but almost no activity in the Atlantic basin. Low-mode dissipation patterns also reflect the uneven distribution of decay rates: dissipation tends to be enhanced over rough or steep topography and equatorward of PSI latitudes. Even though PSI decay times are uniform, there is a slight enhancement of dissipation near the PSI latitudes (∼30°for S 2 and M 2 , 15°for K 1 ) relative to lower latitudes, due to energy travelling from higher latitudes into the region of active PSI.
The decomposition by process of the total dissipation of modes 1-5 is broadly similar across the three tidal constituents (Table 3), though K 1 has somewhat larger percentage dissipation from shoaling and wavewave interactions. Because they act everywhere and powerfully on modes ≥2, wave-wave interactions dominate, causing about 60% of the overall low-mode dissipation. Only 5% of the dissipation occurs in marginal seas shallower than 400 m, implying that most of the lowmode energy approaching continental margins is reflected or dissipated before crossing the shelf break. Shoaling and critical slopes together account for about 25% of the low-mode dissipation.
Guided by the similar dissipation distributions of M 2 and S 2 , we estimate the total dissipation across the eight most energetic tidal constituents as analogous to Eq.
(2). The global integral of D All is 1044 GW, equal to that of G All . Of this total internal tide energy loss rate (Fig. 13f), we estimate that 217 GW or 21% is the local dissipation of modes ≥6 (Fig. 13e). The remaining 827 GW transit through lower modes, but eventually power small-scale turbulence via four modelled routes of the energy cascade (Fig. 13a-d and Table 3): wave-wave interactions (521 GW or 63%), scattering by small-scale seafloor roughness (83 GW or 10%), near-critical reflection (128 GW or 15%) and shoaling (95 GW or 12%). The obtained global maps of internal tide energy sinks possess plausible, distinct patterns that integrate present knowledge of the topography, stratification, internal tide generation, beam propagation and beam degradation.

Sensitivity experiments
Multiple assumptions and approximations underlie the maps of internal tide energy sinks. The angular distribution of sources, attenuation rates by wave-wave interactions and preservation of mode number upon reflection rank among key choices we have made. To explore sensitivity to these choices, we examine four variants of the reference (REF) calculation: • BEAM: a single beam emanates from each source grid cell, in the direction normal to topographic slope.
• FAST: e-folding decay times due to wave-wave interactions are halved.
• SLOW: e-folding decay times due to wave-wave interactions are doubled.
• EXCH: mode n horizontal energy fluxes that reflect at supercritical slopes become 50% mode n and 50% equally partitioned into modes n + 1 to 5.
Changes to the partitioning by process of low-mode dissipation are summarised in Table 3 and changes to the maps of total dissipation illustrated in Fig. 14. The BEAM experiment leads to a more heterogeneous spatial structure of dissipation, but has virtually unchanged process contributions to the overall dissipation. Larger changes occur under modified decay times: wave-wave interactions in the open ocean contribute 51% and 71% of the modes 1-5 dissipation in SLOW and FAST, respectively, bracketing the REF contribution of 61%. Opposite tendencies are implied for the dissipation due to interactions with topography. In particular, the summed percentage contribution of critical slopes and shoaling climbs from 18% in FAST to 31% in SLOW, representing an increase from 150 to 253 GW. Over the factor of 4 range of decay times explored, sensitivity of the process distribution of dissipation is thus sizeable but linear. This overall sensitivity is reflected in the difference maps (Fig. 14e,f): relative to its REF and FAST counterparts, internal tide energy loss in SLOW is less focused around the main source regions and is elevated along basin perimeters.
Reflection of first-mode energy is substantial: in REF, mode-1 reflected fluxes total 88, 24 and 13 GW for M 2 , S 2 and K 1 , respectively. Transferring half of reflected fluxes to higher modes accelerates the degradation of low-mode beams, effectively augmenting damping by wave-wave interactions and reducing energy fluxes impinging on continental slopes. As a result, the EXCH experiment entails modifications analogous-albeit smaller-to those obtained in FAST (Table 3 and Fig. 14d,f). Faster beam degradation also lowers the energy content of C. de Lavergne, et al. Ocean Modelling 137 (2019) 52-75 low-mode internal tides. For example, first-mode M 2 internal tide energy totals 198 PJ in SLOW but only 141 PJ in EXCH and 127 PJ in FAST. Refined constraints on energy content, if attainable, could therefore help to reduce uncertainties in bulk attenuation rates.

Comparison to finestructure observations
Estimates of the dissipation rate of turbulent kinetic energy from measured microscale velocity fluctuations are relatively sparse (Waterhouse et al., 2014). Estimates of the dissipation rate of internal wave energy from the finescale strain observed in CTD profiles are much more abundant (Whalen et al., 2015;Kunze, 2017a). The latter estimates rely on a parameterization Gregg, 1989;Polzin et al., 1995) whose underlying assumptions and implementation choices imply larger uncertainties . For individual strain-based estimates, the uncertainty tied to parameter choices is thought to be a factor of 3 to 5 (Kunze, 2017a;Pollmann et al., 2017). Here we employ the global dataset of Kunze (2017a), who estimated dissipation from vertical strain contained in ship-based hydrographic casts. Four characteristics of this dataset make it most suitable for comparison with present maps of column-integrated internal tide energy loss: (i) wide geographical coverage, (ii) full-depth extent of most profiles, (iii) (theoretical) restriction to energy lost by the internal wave field and (iv) degree of spatio-temporal averaging implicit in the parameterization (Whalen et al., 2015).
The dataset, downloaded from ftp.nwra.com/outgoing/kunze/ iwturb, contains dissipation rates ϵ [W kg −1 ] and diffusivities K ρ [m 2 s −1 ] applying to half-overlapping 256-m segments spanning each hydrographic profile. We first calculated an average dissipation profile within each sampled half-degree grid cell of the WOCE climatology. Values at depths < 400 m are often unreliable due to presence of the surface mixed layer (Kunze, 2017a). Taking a similar approach to Kunze (2017a), we approximated dissipation at depths z < 400 m as ϵ (z) = K ρ (z = 400 m)N 2 (z)/R f , where N is the local climatological buoyancy frequency, K ρ the finestructure-inferred diffusivity, and R f = 0.2 the mixing efficiency assumed in the finescale parameterization. This conservative prolongation of dissipation profiles to the surface rests upon the scaling ϵ(z) ∝ N 2 (z), typical of internal wave C. de Lavergne, et al. Ocean Modelling 137 (2019) 52-75 breaking in the upper ocean (Gregg, 1989), but could underestimate the shallow dissipation. In particular, breaking of near-inertial internal waves generated by atmospheric storms (Pollard, 1970) is expected to be underrepresented in the profiles-an expectation reinforced by the bias of ship-based sampling toward moderate winds. Next, we dismissed profiles that cover < 90% of the total water column stratification (i.e., such that < N dz N dz 0.9 2 2 R ). This criterion excludes 6% of profiles. We finally integrate dissipation ϵ (times the seawater density ρ) over the whole water column to obtain depth-integrated internal wave energy dissipation rates.
The spatial pattern of the finestructure-inferred dissipation exhibits strong similarities to the mapped internal tide energy loss (Fig. 15): dissipation in the Indian Ocean is focused in the western half of the basin, generally decreasing away from Madagascar; dissipation in the Pacific is focused in the western tropics, and relatively weak along a crescent joining the northern, eastern and southern portions of the deep basin; dissipation in the Atlantic is relatively high above the Mid-Atlantic Ridge but tends to have weaker large-scale contrasts than in the Indo-Pacific; dissipation in the Southern Ocean is generally weak, except for hotspots near the Southwest Indian Ridge, Macquarie Ridge and Drake Passage. A more quantitative comparison is presented in Fig. 16a for the global ocean and in Fig. 17 for separate basins. The qualitative similarities outlined above translate into a global correlation of 0.41 between the finestructure and REF dissipations (when computed C. de Lavergne, et al. Ocean Modelling 137 (2019) 52-75 using the logarithm of dissipation values, the global correlation coefficient becomes 0.55). The correlation is strongest in the Pacific and Indian oceans, and weakest in the Atlantic. All correlation coefficients have high statistical significance (p-value indistinguishable from zero). Encouraging agreement in magnitudes is also observed (Figs. 15c and 16a): 65% (80%) of values agree within a factor of 3 (factor of 5). Clear discrepancies also emerge from the comparison. First, the mapped tidal dissipation displays a slight pattern amplification compared to the observational dataset: the REF experiment generally predicts somewhat larger rates in the most dissipative regions, and somewhat weaker rates in more quiescent open-ocean areas (Fig. 15c). This could indicate a lack of energy redistribution in our calculation, possibly related to overestimated transfers to dissipation. Biases in the finestructure estimates, notably in regions of strong forcing (Hibiya et al., 2012;Waterman et al., 2014;Thurnherr et al., 2015;Bouruet-Aubertot et al., 2018), could also contribute to the pattern of differences. Superimposed on this broad pattern is a second, more significant mismatch: the REF dissipation is substantially lower than that estimated by Kunze (2017a) in five important regions, namely the western boundary of the Atlantic, the Agulhas region, the region west and south of Australia, the Kuroshio region, and the southeast Pacific sector of the Southern Ocean (Fig. 15c). These discrepancies are most pronounced in the Atlantic and Southern Oceans, leading to weaker overall agreement compared to the Indo-Pacific (Fig. 17).
The large underestimates by REF of the finestructure-inferred dissipation consistently occur in the regions of atmospheric and/or oceanic storm tracks. This suggests that non-tidal forcings or non-modelled interactions dominate in these areas. Possible contributors to the observed dissipation include: (i) breaking of wind-powered near-inertial waves, perhaps catalysed by mesoscale ocean eddies (Kunze, 1985;Zhai et al., 2007;Jing and Wu, 2014); (ii) high-mode internal wave generation by eddies interacting with topography, notably along western boundaries (Zhai et al., 2010;Nikurashin and Ferrari, 2010;Clément et al., 2016;Pollmann et al., 2017); (iii) scattering of low-mode internal tides by mesoscale eddies (Dunphy and Lamb, 2014). To evaluate the ability of the energy source (i) to bridge the gap between the REF and finestructure dissipations, we employed an estimate of the atmospheric power input to near-inertial motions (Rimac et al., 2013), and assumed that a 10%-20% fraction of this power fuels near-inertial waves that break in the near-local water column (Furuichi et al., 2008;Zhai et al., 2009;Cuypers et al., 2013;Jouanno et al., 2016). This additional dissipation barely modifies the overall correlation with finestructure data, but raises the agreement within a factor of 3 (factor of 5) to 72-74% (88%) of all values (Fig. 16c,d). Order-of-magnitude underestimates only persist along the western boundary of the Atlantic (not shown), calling for additional contributions there.
Importantly, the REF mapping of internal tide energy dissipation has not been tuned to match the finestructure data. Other experiments compare similarly with the observational dataset (Table 4); agreement within a factor of 3 ranges from 63% in BEAM and FAST to 67% in SLOW. Hence, comparison lends confidence that both the present approach and the finescale parameterization have skill in mapping the energy flow from internal waves to small-scale turbulence. To put this skill in context, we mapped the dissipation implied by a standard OGCM parameterization of internal-wave-driven mixing (Fig. 14b), and compared it to the finestructure dataset (Fig. 16b). As in the NEMO model (Mignot et al., 2013), near-field mixing is parameterized using one-third of the local power input to internal tides (Nycander, 2005), and far-field mixing is represented by a fixed diffusivity, equal to 10 −5 m 2 s −1 except for a reduction around the equator and under sea ice. Applied to the WOCE climatology, this OGCM parameterization implies a total power consumption of 920 GW, and generally yields order-of-magnitude agreement with the observational dissipation rates (Fig. 16b). However, the depth-integrated dissipation implied by the OGCM parameterization varies mostly with latitude and misses the large spatial contrasts of far-field mixing evident in the observations.

Comparison with previous budgets
Whether continental slopes and shelves host an important proportion of low-mode internal tide energy loss remains a matter of debate. Kelly et al. (2013) argued that the first-mode semidiurnal tide loses the bulk of its energy at large-scale topographic obstacles. In contrast,  estimated that wave-wave interactions and scattering off abyssal hills in the open ocean contribute 90% of the dissipation of the first-mode M 2 tide. Our estimate lies in between: interactions with large-scale topographic obstacles account for about half of the dissipation of the mode-1 M 2 tide ( Table 2). The lesser role attributed to basin margins by  stems from different parameterization choices: they estimated the open-ocean attenuation timescale to be O (1 day), compared to O (10 days) here. The latter timescale seems to better accord with the observed long-range propagation of the first-mode M 2 internal tide (Zhao and Alford, 2009;Dushaw et al., 2011;Waterhouse et al., 2018).
However, we emphasize that dissipation characteristics of the first vertical mode are not representative of the overall non-local internal tide dissipation. In agreement with mooring observations (Zhao et al., 2010;Vic et al., 2018), we find that modes 2-5 also accomplish important horizontal energy redistribution. Yet because of stronger damping by wave-wave interactions, this redistribution occurs over distances that are one to two orders of magnitude smaller than that associated with mode 1. As a result, dissipation in the basins' interior, away from continental slopes and shelves, accounts for 65-80% of the total modes 1-5 dissipation of 827 GW and is dominated by wave-wave interactions ( Table 3). The presence of substantial low-mode dissipation in the open ocean concurs with inferences from high-resolution OGCM simulations (Ansong et al., 2015;Buijsman et al., 2016).
According to present calculations, the proportion of low-mode dissipation occurring in marginal seas shallower than 400 m is only 4-6% (Table 3). A comparison between available microstructure measurements of water-column dissipation rates and estimated internal wave generation rates suggested a higher proportion of ∼30% (Waterhouse et al., 2014). Nevertheless, it was recognized that this inferred residual may also incorporate under-sampled continental slopes. The higher proportion could also reflect in part a bias of microstructure sampling Table 3 Low-mode dissipation per process, across tidal constituents and sensitivity experiments. The percentage contribution of each process to the dissipation of low-mode (modes 1-5) internal tides is given for each constituent in the REF experiment, and for all constituents across sensitivity experiments. Shelves, defined as areas where the bottom depth is < 400 m, have been separated out. 'All' incorporates the eight most energetic tidal constituents following Eq. (10) C. de Lavergne, et al. Ocean Modelling 137 (2019) 52-75 toward regions of strong forcing, where the ratio of water-column dissipation to generation tends to be low (Kunze, 2017a). Fig. 18 shows the distribution of the ratio of internal tide dissipation to generation, D/ G, according to the REF calculation. Because sources are more localised than sinks, the ratio exceeds 1 over the vast majority of the ocean area.
The bulk (64%) of the dissipation also occurs at ratios above 1. In contrast, 72% of the total internal tide generation occurs in locations where generation exceeds dissipation. These distributions of D/G reflect the substantial energy redistribution achieved by low-mode beams, and illustrate the difficulty of predicting dissipation rates from generation rates alone.

Conclusions
Mapping the breaking of internal tides is a long-standing challenge of oceanography (Munk, 1966). Here we proposed a framework for estimating the geography of internal tide energy sinks and applied it to the contemporary ocean. The framework builds upon the pioneering work of , recent calculations of internal tide generation projected onto vertical normal modes (Falahat et al., 2014b), and the growing understanding of the downscale cascade of internal tide energy (MacKinnon et al., 2017). It consists of tracking every internal tide energy beam from generation to dissipation, through observed stratification, accounting for refraction and for reflection against supercritical topography. Using the WOCE global hydrographic climatology, generation maps for the three most energetic tidal constituents (M 2 , S 2 and K 1 ), and simplified representations of dissipation via interactions with topography and with the ambient wave field, we constructed a global horizontal map of internal tide energy loss (Fig. 13f) decomposed into contributing dissipative processes (Fig. 13a-e).
Different breaking pathways produce distinct vertical structures of dissipation. By coupling each map of column-integrated energy dissipation with an appropriate vertical structure, three-dimensional maps of dissipation and diffusivity may be built. Such maps could then serve to refine estimates of water mass transformation in the ocean interior (e.g., de Lavergne et al., 2017;Kunze, 2017b), to constrain inverse estimates of ocean circulation (e.g., Ganachaud and Wunsch, 2000;Groeskamp et al., 2017), or to parameterize internal tide-driven mixing in OGCMs (e.g., Melet et al., 2016). These applications may in turn illuminate the impact of each dissipative process on ocean circulation and the sensitivity of circulation to the assumptions and choices made in the mapping procedure.
Two principal limitations of the framework must be underlined. First, the choice to track low-mode energy beams one by one-and then sum over all beams to recover the total dissipation-disallows a straightforward treatment of nonlinear interactions between beams. For example, if the efficiency of a dissipative process at some location In (b) we illustrate the dissipation implicit in a standard OGCM parameterization of internal wave-driven mixing. This is one-third the total power input to internal tides (Nycander, 2005) plus the column-integral of 6ρK b N 2 , where K b is a fixed background diffusivity equal to 10 −5 m 2 s −1 except for a reduction around the equator and under sea ice. Dissipation is in units of 10 −3 W m −2 .
depends on the total low-mode energy at this location, such dependence cannot be easily factored in. Only the energy carried by the beam being tracked has thus been used to diagnose its contribution to dissipation. In return, the adopted Lagrangian approach avoids the numerical diffusion/dispersion of energy fluxes produced by Eulerian advection schemes (see Appendix B). This is a key advantage if the trajectories and energy levels of observed beams are to be reproduced, and if the eddy-modulated dispersion of energy fluxes is to be explicitly represented.
Second, implementation of the Lagrangian calculation within a running OGCM is fraught with challenges, related to the number of beams and to the complexity and stratification-dependence of interactions with topographic slopes. Presently, the framework is designed to produce static dissipation maps appropriate to a given mean ocean state. Several elements suggest that such static maps can still serve OGCM applications: (i) they can replace maps of internal tide generation as the power available to the model's tidal mixing parameterization; (ii) such a parameterization could obviate the need for a fixed background diffusivity, ensuring mixing is energy constrained, and obviate the computational cost of an interactive representation of low-mode energy propagation; (iii) we expect the sensitivity of the global horizontal distribution of dissipation to changes in stratification to be modest, as is the case of internal tide generation (Egbert et al., 2004); (iv) the vertical structure of dissipation, most crucial for ocean Fig. 15. Column-integrated internal wave energy dissipation predicted (a) by the present REF calculation and (b) by a finescale parameterization applied to shipbased CTD data (Kunze, 2017a). In (c) we show the log of the ratio of (a) to (b).
ventilation (de Lavergne et al., 2016a), can be coupled to the evolving model stratification; (v) the maps could be updated at distant times of a simulation if large changes in the simulated ocean state warrant doing so.
Limitations of the present application to the modern ocean are numerous. Primary among them is the ad hoc representation of low-mode energy loss through wave-wave interactions: we chose decay times that depend only-and crudely-on latitude and mode number. Since this a nonlinear process, the decay time should in reality also depend on the amplitude of the internal wave field. Further work is required to assess our choices and incorporate more realistic variations of decay rates, so as to reduce the uncertainty attached to this dominant sink of internal tide energy. The modelled interactions with topographic slopes also suffer from important approximations. In particular, we did not account for the dependence on Froude number of the near-critical slope range and of the fractional energy loss due to shoaling (Legg, 2014). Additional uncertainty originates from the estimation of energy sources, which demanded a correction for the negative values present in modal generation fields (Falahat et al., 2014b) and an arbitrary choice for the angular distribution of energy emission. Substantial uncertainty in the dissipation maps also stems from processes that have not been explicitly modelled. Not least are interactions between low-mode beams and geostrophic flows, whose role in the spreading and downscale cascade of energy (Rainville and Pinkel, 2006;Dunphy and Lamb, 2014;Buijsman et al., 2017) deserves more attention.
In spite of these uncertainties, the estimated maps of internal tide energy sinks appear to be plausible and broadly compatible with available observations. Most of the energy generated at Hawaii radiates  Kunze (2017a); the power input to near-inertial motions is from Rimac et al. (2013). away and feeds dissipation across the Pacific (Klymak et al., 2006;Alford et al., 2007). Most of the energy travelling from Macquarie Ridge toward Tasmania reflects or dissipates at the Tasman continental slope (Johnston et al., 2015;Waterhouse et al., 2018). Most of the dissipation in the eastern Brazil Basin originates from locally generated high modes (Polzin, 2004). Indonesian seas host 10% of both the total source and total loss of internal tide energy (Koch-Larrouy et al., 2007). Patterns and magnitudes of the mapped internal tide dissipation compare well with observational internal wave dissipation estimates along hydrographic sections of the World Ocean (Kunze, 2017a). Regional discrepancies with the latter dataset are consistent with significant local dissipation of wind-generated near-inertial waves (Jing and Wu, 2014) and are suggestive of mesoscale eddy energy loss to internal waves near western boundaries (Clément et al., 2016). On the whole, these results bolster our understanding of tidal dissipation.

Table 4
Comparison to finestructure observations across sensitivity experiments. Coefficients of horizontal correlation between finestructure-based dissipation rates and those predicted by each experiment, computed using raw dissipation values (first row) or the logarithm of dissipation values (second row). Percentage agreement within a factor of 3 (factor of 5) is given in the third row (fourth row). See section 3.3 of the text for a definition of sensitivity experiments. The geometric treatment of interactions with topographic slopes resolved in the 'etopo2v2' bathymetry product follows the rationale of Fig. 7, accounting for topographic complexities and depth variations of the ray-path slope. Here we detail the calculation of relevant fractions considering the example of a grid column (I, J) of the WOCE climatology whose eastern neighbour (I + 1, J) is shallower, as sketched in Fig. A1. The thickness H of the column is the average of the 15×15 'etopo2v2' bathymetric values contained in the half-degree cell (I, J). The buoyancy frequency N is known at each level of the WOCE climatology and extrapolated to the bottom where necessary. The energy flux F b of a beam exiting the grid square (I, J) through its boundary with (I + 1, J) interacts with longitudinal topographic segments connecting H(I, J) to H(I + 1, J). To determine the fraction of F b encountering critical or supercritical topography, we begin by calculating the wave slope s at each level using Eq. (8). Next, we apply the same

Ocean area
Internal tide generation Internal tide dissipation Fig. 18. Distribution of the ratio of internal tide dissipation to generation. Histograms display the proportion of (grey) ocean area, (blue) internal tide generation G All and (red) internal tide dissipation D All characterised by various ranges of the logarithm of the ratio D All /G All , according to the REF calculation. Because of energy redistribution by low-mode beams, generation occurs mostly at ratios < 1, whereas dissipation occurs predominantly at ratios in excess of 1.

Fig. A1.
Example row of topography segments spanning two grid columns of the WOCE climatology. Topography at 0.5°resolution is represented by the grey shading. Topography at 1/30°resolution is represented by the broken blackline. The latter has 29 segments, each characterised by a horizontal length Δx and a positive or negative height increment Δh. Light horizontal lines schematize vertical levels of the WOCE climatology. Points A and B are as defined in the text. • From west to east, find the first bathymetric value less than H(I, J) and the first value less than H(I + 1, J), or points A and B, respectively ( Fig. A1).
• Ignore segments west of A and east of B. This ensures F b does not interact twice with the same topography.
• Calculate the effective horizontal projection, or projected length L x = Δx + Δh/s, of each topographic segment. Here Δx and Δh are the horizontal length and the (positive or negative) eastward height increment of a single segment, respectively.
• From west to east, find the first segment that has a negative projected length (point R). If such a segment exists, sum L x cumulatively eastward until it becomes positive (point R ′ ). Set all L x values to zero between R and R ′ . Repeat eastward from R ′ . This procedure accounts for shadowing, i.e. for the presence of downsloping segments invisible to rays.
• Find topographic segments whose slope Δh/Δx lies between 0.8s and 1.5s, and sum their projected lengths L x to obtain the critical projected length of the whole row.
• From east to west, find the first segment whose slope exceeds 1.5s. Follow a ray path backward from the top of this segment (point S) until it encounters the seafloor (point S ′ ). Sum all non-critical projected lengths L x between S and S ′ . Repeat westward from S ′ . Cumulate the summed projected lengths to obtain the supercritical projected length of the whole row.
The critical (supercritical) projected length is then averaged over the 15 rows and divided by the distance between two consecutive bounces, 2∫ 0 H(I,J) s −1 dz, to yield the fraction of F b dissipated (reflected) at the slope. The remaining, transmitted flux loses another fraction of its magnitude through shoaling. The latter fraction is estimated as (∫ H(I+1,J) H(I,J) Ndz/∫ 0

H(I,J)
Ndz) 2 , a crude approximation based on the model results of Legg (2014).
Fractions are calculated for each tidal constituent, each half-degree grid cell and each of the four cardinal directions. They are set to zero in the directions along which H increases, except for crest grid cells where H increases in opposite directions. Indeed, we find that grid columns straddling the crests of mid-ocean ridges host significant critical dissipation: their relatively high bottom stratification implies relatively small wave slopes that often match those of the topography. To obtain critical and supercritical projected lengths at a crest cell, the procedure detailed above is applied to the 15 rows of 14 topography segments contained in the cell, defining point B as the shallowest topographic point in each row. The distance between consecutive bounces is calculated as previously. The fraction dissipated through shoaling becomes (∫ H B (I,J) H(I,J) Ndz/∫ 0 H(I,J) Ndz) 2 , with H B the bathymetric depth at point B averaged over the 15 rows. Premises implicit in this local geometric analysis include: (i) 0.5°is an appropriate horizontal resolution to estimate transmission/reflection/ dissipation fractions; (ii) internal tides behave as individual rays of infinitesimal width that are randomly distributed in the horizontal; (iii) effects of non-normal incidence can be neglected. Premises (i) and (ii) relate to the horizontal scale most relevant to interactions between a low-mode beam and a topographic obstacle. This scale is thought to be O (10-100 km), comparable to the wavelength of the incident low-mode internal tide Kunze and Llewellyn Smith, 2004). This suggests that 0.5°is an appropriate resolution for the analysis, but that accounting for topographic variations down to a horizontal resolution of ∼ 3 km may be superfluous. The calculated fractions should thus be interpreted as best estimates of the likelihood that transmission/reflection/dissipation dominates in any half-degree grid cell. We also note that the effective resolution of the 'etopo2v2' bathymetry map is closer to 10 km over most of the ocean (Falahat et al., 2014b). Premise (ii) also entails that we ignore complexities related to the interplay between rays and modes, which can modulate interactions with topography (Klymak et al., 2011). Premise (iii) is a practical choice dictated by the longitude-latitude grid of bathymetry and stratification climatologies. Fractions in non-cardinal directions could be estimated using bathymetric fields interpolated onto rotated grids. However, such added complexity was deemed unwarranted here given other uncertainties in the global mapping, and given evidence that oblique incidence has only a modest overall impact .

Appendix B. Lagrangian versus Eulerian propagation
Consider a beam of initial magnitude F b and orientation ϕ. In the Lagrangian framework, the beam propagates from one grid cell to the next, changing direction and magnitude as it interacts with the environment, but retaining a single orientation and infinitesimal width at any given time.
In an Eulerian framework, energy transports are not attached to beams but instead to a discrete three-dimensional (longitude, latitude, angle) grid. A numerical scheme propagates F b from the source cell to neighbouring cells, depending on the beam angle ϕ. Unless ϕ matches exactly one of the grid's underlying directions, the beam will necessarily be divided between two grid cells, each receiving a fraction of F b . Further splitting will occur at the next iteration, resulting in numerical spreading of energy transports. Refraction also causes numerical diffusion: unless the shift in angle over one time step is an exact multiple of the angular resolution, energy transports split into adjoining discrete angles. Numerical choices (advection schemes, resolution) can alleviate but not eliminate diffusion/dispersion inherent to discretisation in the Eulerian framework.
To illustrate these effects, we calculated M 2 mode-1 energy transports in the BEAM experiment using an Eulerian (Fig. B1b) instead of the Lagrangian (Fig. B1a) propagation scheme. In the Eulerian framework, all beams propagate simultaneously, feeding a three-dimensional field of energy fluxes F x y ( , , ). F is iteratively updated until steady-state is achieved, i.e. until the global dissipation rate matches the global generation rate. To mitigate numerical diffusion, the propagation scheme allows direct communication between cardinal as well as diagonal neighbours: over one iteration, a given beam feeds two among its eight contiguous grid cells, depending on its orientation ϕ. Substantial numerical spreading is nonetheless apparent in the comparatively diffuse distribution of energy transports produced by the Eulerian scheme (Fig. B1). Furthermore, spreading is not uniform over all orientations: transports that retain a beam-like structure are those oriented in cardinal or diagonal directions, along which propagation is less dispersive. This artificial segregation of beams depending on their orientation is amplified at higher angular resolution, offsetting potential gains of resolution increases.