Incorporating tides and internal gravity waves within global ocean general circulation models: A review Progress in Oceanography

Until recently, high-resolution global modeling of tides has been done separately from high-resolution global modeling of the atmospherically-forced oceanic general circulation. Here we review the emerging class of high-resolution global models that are simultaneously forced by both atmospheric fields and the astronomical tidal potential. Such models simulate barotropic (surface) tides, internal tides, near-inertial motions, the eddying general oceanic circulation, and a partially resolved internal gravity wave (IGW) continuum spectrum (Garrett- Munk spectrum) simultaneously. We review the technical aspects of such global models and their myriad ap-plications, for example, in satellite oceanography, operational oceanography, boundary forcing of regional models, tidal-cryosphere interactions, and assessment of future coastal flooding hazards in a changing climate with altered tides.


Introduction
Tides have fascinated humankind for thousands of years (Cartwright, 1999). Accurate analysis and prediction of tides in ports and harbors has been practiced since the 1800s. The advent of the TOPEX/JASON satellite altimetry series (Fu & Cazenave, 2001;Stammer & Cazenave, 2017;International Altimetry Team, 2021) revolutionized physical oceanography by allowing for highly accurate global maps of surface tidal elevations (e.g., Le Provost, 2001; and other oceanic motions. The TOPEX/JASON tidal maps build upon work demonstrating the feasibility of tidal mapping with the earlier Seasat and Geosat altimeters (e.g., Cartwright, 1983;Mazzega, 1985;Woodworth & Cartwright, 1986;Cartwright & Ray, 1990). Despite the advent of accurate tidal maps from altimetry, many questions about tides remain, including the response of tides to climate change, the details of tide-cryosphere interactions and open-ocean tidal dissipation mechanisms, the predictability of internal tides, and the impact of tides on upcoming remote sensing missions. Some of these questions are addressable with global high-resolution simulations that incorporate tidal and atmospheric forcing simultaneously, a frontier direction in ocean modeling that is the subject of this paper.

Background and definition of terms
Tides are best known for effecting a periodic rise and fall of the sea surface height (SSH). The large-scale motions that dominate SSH are variously known as surface tides or barotropic tides. The surface tides are arguably weakly resonant in the open-ocean (e.g., Wunsch, 1972;Garrett & Greenberg, 1977;Heath, 1981;Arbic et al., 2009), because the spatial structures and frequencies of oceanic normal modes (e.g., Platzman et al., 1981;Platzman, 1991;Zahel & Müller, 2005;Müller, 2007) match those of the astronomical forcing reasonably well. The weakly resonant open-ocean tides drive strong tides in coastal regions (e.g., Hudson Strait, the Bay of Fundy, the Patagonian and Northwest European Shelves, and others) that are well-shaped for resonance (e.g., Garrett, 1972;Clarke, 1991;Cummins et al., 2010).
The barotropic SSH perturbations are driven by converging and diverging barotropic (depth-independent) flows. When barotropic tidal currents flow over topographic features they initiate vertical motions, which in a stratified fluid such as the ocean implies the existence of tidal oscillations along the interfaces between layers of different densities (e. g., Bell, 1975;Wunsch, 1975;Baines, 1982;St. Laurent & Garrett, 2002;Lahaye et al., 2020;amongst many). The heaving of interfaces below the sea surface at tidal frequencies, known as the internal or baroclinic tide, is associated with temperature perturbations, as we will show later in frequency spectra of temperature, and brings about a considerable baroclinic (depth-dependent) component in tidal velocities (see, for instance, Figs. 2 and 5 of Timko et al., 2013). Internal tides can be divided into stationary (also referred to as "coherent") components, which can be defined with a predictable amplitude and phase lag, and non-stationary (also referred to as "incoherent") components, which lie within the tidal frequency band but which do not have a predictable amplitude and phase lag. The largest internal tide vertical displacement signals take place deep within the thermocline, but internal tides are also associated with an SSH perturbation signal, of smaller amplitude than the barotropic tide SSH signal, and with much smaller horizontal scales. Mitchum (1996, 1997) extracted stationary internal tide SSH signals from satellite altimetry; they separated internal tides from barotropic tides by applying a spatial high-pass filter to the total (barotropic plus baroclinic) tidal SSH signals. (In some models, internal tide SSH signals have been isolated from barotropic tide SSH signals through division of SSH into steric and non-steric components-see, for instance, Savage et al., 2017a). Ray and Mitchum's pioneering work on extraction of internal tides from altimetry has been built upon in numerous studies, including Kantha & Tierney (1997), Zhao et al. (2010), Zaron (2019), Ubelmann et al. (2022), and many others. See Carrère et al. (2021) for a recent review of empirical internal tide models based upon altimetry.
Empirical internal tide models, based upon altimetry, have demonstrated that low vertical-mode internal tides propagate over thousands of kilometers, consistent with evidence from acoustic tomography measurements (Dushaw et al., 1995) and moorings (e.g., Alford et al., 2007;Zhao et al., 2010). Low vertical mode internal tides are primarily generated by topographic features such as the Hawaiian islands, which have horizontal scales of about 100 km. High vertical mode internal tides are generated by small-scale (e.g., ~10 km) features such as abyssal hills. Vertical modes are solutions to the Sturm-Liouville problem for internal gravity waves (IGWs), wave motions in a stratified fluid with a gravitational restoring force. See, for instance, Nugroho (2017), Nugroho et al. (2018), the appendix of Arbic et al. (2018), and other sources. Higher vertical modes have more structure in the vertical direction, and shorter horizontal scales, than low vertical modes.
Spectra of kinetic energy, temperature variance, and other oceanographic properties indicate a continuum of internal wave energy extending out to high frequencies (N) and high vertical wavenumbers (Garrett & Munk, 1975;Cairns & Williams, 1976). The IGW continuum spectrum, or Garrett-Munk spectrum, is thought to emerge from nonlinear interactions between internal gravity waves, with internal tides and NIWs serving as the primary energy sources (e.g., Olbers, 1976;McComas & Bretherton, 1977;McComas & Müller, 1981a, 1981bMüller et al., 1986;Polzin, 2004). The final state of the cascades in the IGW spectrum is short vertical scale IGWs, which turn over and break, causing a substantial amount of the mixing in the interior of the open ocean-that is, away, from surface and bottom boundaries (Moum et al., 2003;Klymak et al., 2008). Indeed, observed energy dissipation rates are predicted remarkably well by theories based on observed IGW energy levels and an assumed cascade to smaller vertical scales at which breaking takes place (Gregg, 1989;Polzin et al., 1995;Whalen et al., 2015;Kunze, 2017a).
Low-vertical-mode oceanic internal gravity waves have horizontal spatial scales of order 100 km, similar to the length scales of oceanic mesoscale eddies, the oceanic dynamical equivalent of atmospheric synoptic weather systems. Mesoscale eddies contain most of the ocean's kinetic energy, and are ubiquitous throughout the ocean, but are strongest in regions of strong currents including the Antarctic Circumpolar Current and the western boundary currents such as the Gulf Stream, Kuroshio, etc. (e.g, Stammer, 1997;Wunsch, 1997;amongst many). Interactions between internal tides and mesoscale eddies will be discussed further below. The similar length scales of internal tides and mesoscale eddies means that they can be entangled in some measures, such as the SSH measurements of satellite altimetry (e.g., Ray and Byrne 2010;Shriver et al., 2012;Zaron, 2017;Nelson et al., 2019).

Motivation for simulating tides and IGWs within global ocean models
Until recently, high-resolution global modeling of tides was performed separately from high-resolution global modeling of the oceanic general circulation. Models with simultaneous tidal and atmospheric forcing can be used to examine the interactions of tides with other components of the climate system, for instance, mesoscale eddies, the oceanic general circulation, sea ice, ice streams, floating ice shelves, marine ecosystems, river outflows and estuaries, etc. As an example of Fig. 1. Schematic of open-ocean internal wave generation and dissipation processes that were considered as part of a Climate Process Team led by Jennifer MacKinnon of the University of California San Diego. Barotropic tidal flow interacts with topographic features to generate highvertical-mode internal waves (e.g., at mid-ocean ridges) and low-vertical-mode internal waves (e. g., at tall steep ridges such as the Hawaiian Ridge). Deep low-frequency flows, including mesoscale eddies and currents, over topographic features can generate lee waves (e.g., in the Southern Ocean). Storms cause near-inertial oscillations in the mixed layer, which can generate both low-and high-mode internal waves (e.g., beneath storm tracks). In the open ocean these internal waves can scatter off of topographic features and potentially interact with mesoscale fronts and eddies, until they ultimately dissipate through wave-wave interactions that generate a nonlinear cascade to small-scale turbulence. Internal waves that reach the shelf and slope can scatter, dissipate via bottom boundary layer drag, or amplify as they propagate towards shallower water. Reproduced from Fig. 1  such an interaction, in coastal areas with strong tides, tidal mixing fronts in sea surface temperature arise in the summertime (Simpson & Hunter, 1974). Tidal mixing fronts represent boundaries between stratified and tidally mixed waters, and feature enhanced biological productivity. Tides affect Arctic sea ice drift and deformation (Hibler et al., 2006), the stick-slip flow of continental ice streams (Bindschadler et al., 2003), and the horizontal flows of floating ice shelves (Doake et al, 2002). Tides in sea ice-covered oceans have been examined in several studies (Kowalik & Proshutinsky, 1994;Lyard, 1997;Albrecht et al., 2006;Holloway & Proshutinsky, 2007;Cancet et al., 2018). Seasonal changes impact openocean internal tides , barotropic tides on shelves (e. g., Pugh and Vassie, 1992;St-Laurent et al., 2008;Müller, 2012;Gräwe et al., 2014), and open-ocean barotropic tides . Tides change over centennial scales (Haigh et al., 2020), likely due to sea level rise (Müller et al., 2011;Schindelegger et al., 2018) and other climatic changes. Secular changes in tides are often comparable to changes in mean sea level (Jay, 2009) and therefore should be considered in assessments of coastal flooding hazards in a changing climate system. Ruault et al. (2020) demonstrated that inclusion of tides improves the simulation of the offshore spreading of river plumes. All of the tidal interactions with other climate system components can at least in principle be examined in high-resolution global models with simultaneous tidal and atmospheric forcing. In another application of general circulation models with embedded tides, Ray et al. (2021) constructed bounds for the small ψ 1 ocean tide, which allows for estimation of the ψ 1 earth tide, a window into geophysical processes such as free core nutation.
As was first documented in Müller et al. (2015), the internal gravity wave continuum spectrum, or Garrett-Munk spectrum (Garrett & Munk, 1975;Cairns & Williams, 1976), is partially resolved in global highresolution models with simultaneous atmospheric and tidal forcing. Such models carry the "ingredients" necessary for the development of the IGW spectrum via the classical wave-wave interaction paradigm-NIWs forced by the fast component of wind stress, internal tides forced by barotropic tidal flow over topography, and nonlinear interactions enabled by the high resolution in the models. (See Barkan et al., 2017, for an alternative view implicating the potential importance of mesoscale eddies in generating the IGW continuum spectrum). Below, we will often refer to such models as "global internal wave models" as a shorthand for "global high-resolution models with simultaneous atmospheric and tidal forcing".
Global internal wave models are important for operational and satellite oceanography. Because tides are present in many oceanic observations, including frequency spectra of SSH variance, kinetic energy, and temperature variance, it stands to reason that operational models will benefit from inclusion of tides. A major new project funded by the US Navy is using global internal wave models to examine the impact of internal waves on basin-to global-scale acoustic propagation. Another recent National Oceanographic Partnership Program (NOPP) project will focus on large arrays of in-situ instruments to test global internal wave models. Satellite altimetry is ingested into global operational models such as the US Navy system, based upon the Hybrid Coordinate Ocean Model (HYCOM; Bleck, 2002). Altimetry enables operational models to accurately place mesoscale eddies at their correct locations . The Surface Water Ocean Topography (SWOT) mission Morrow et al., 2019), scheduled for launch in late 2022, will measure SSH at finer spatial resolutions than is possible with current generation nadir altimeters, and will measure in twodimensional swaths rather than the one-dimensional tracks of nadir altimeters. In order for SWOT to fully realize its potential for the study of mesoscale and submesoscale eddies, high-frequency internal wave motions will have to be accurately removed, just as barotropic tides have been accurately removed from nadir altimetry (Egbert et al., 1994;Le Provost et al., 1994;Shum et al., 1997;Ray, 1999;Egbert & Erofeeva, 2002;Lyard et al., 2006Lyard et al., , 2021Stammer et al., 2014). As will be discussed in Section 5.4, global internal wave models have been used extensively to quantify the relative contributions of internal tides and gravity waves to the SSH wavenumber spectrum, in preparation for SWOT. Global internal wave models have also been used to quantify the kinetic energy in geostrophically balanced vs. internal wave motions, in preparation for proposed satellite missions that will measure surface ocean velocity Rodriguez et al., 2020), and the predecessor airborne Sub-Mesoscale Ocean Dynamics Experiment (S-MODE) mission (Rodriguez et al., 2020). Another application of global internal wave models, likely to grow in importance in the coming SWOT era, lies in regional modeling, which traditionally employs only barotropic tidal forcing at the lateral boundaries. Recent studies, to be discussed in Section 5.5, emphasize the importance of remotely generated internal tides and waves for the lateral boundaries of regional models.
Finally, global internal wave models are likely to be valuable in the quantification of mixing. The array of IGW production and dissipation mechanisms is complex, as illustrated in the schematic shown in Fig. 1. Quantifying the mechanisms and four-dimensional space-time geography of internal tide and IGW dissipation (e.g., Whalen et al., 2012Whalen et al., , 2015Kelly et al., 2013;Whitehouse et al., 2014;MacKinnon et al., 2017) and the impact of open-ocean tidal dissipation on the large-scale oceanic circulation (e.g., Munk and Wunsch, 1998;Wunsch & Ferrari, 2004;Simmons et al., 2004a;St. Laurent & Simmons, 2006;Ferrari & Wunsch, 2009;Kunze, 2017b, amongst many) are matters of intense current research interest in the physical oceanography and climate modeling communities. Global internal wave models cannot resolve breaking IGWs, but they can provide information about generation and propagation and can therefore contribute to the discussion of internal wave dynamics. An initial attempt to look at the mixing implied by the energetics of HYCOM internal tide simulations (Buijsman et al., 2016) will be discussed later in this paper.
Early modeling of baroclinic (internal) tides was done in limited domains. For example, internal tides in the Northwest Australian Shelf, northern British Columbia, Line Islands Ridge, East China Sea, and the mid-Atlantic Ridge were examined by Holloway (1996), Cummins & Oey (1997), Johnston et al. (2003), Niwa & Hibiya (2004), and Zilberman et al. (2009), respectively. Xing & Davies (1999 outlined several factors influencing the internal tides in a regional model of the shelf off the west coast of Scotland. Robertson performed regional internal tide studies in the Weddell Sea (Robertson, 2001a(Robertson, , 2001b, Fieberling Guyot (Robertson, 2006), and Indonesian Seas (Robertson & Ffield, 2008). The internal tides near Hawai'i have been simulated in many regional studies (e.g., Kang et al., 2000;Merrifield et al., 2001;Merrifield & Holloway, 2002;Carter et al., 2008). Carter et al. (2012) provides a review of regional internal tide studies up to that point. For an example of a more recent regional study, conducted with increased computer power, see Nugroho (2017) and Nugroho et al. (2018).
The first basin-and global-scale simulations of internal tides were performed with "tide-only" models, in which atmospheric forcing fields were absent, the stratification was fixed across the model domain, and only tidal forcing was applied. A North Pacific basin simulation of the internal tide was conducted by Niwa & Hibiya (2001), and the first global simulations followed in Arbic et al. (2004a), Simmons et al. (2004b). More recent papers based on global "tide-only" models with fixed stratification include Hibiya (2011, 2014), who examined the sensitivity of internal tide generation to model resolution, Kodaira et al. (2015), who examined the baroclinic tide contribution to Fig. 2. Global map of amplitude (cm) of M 2 surface tidal elevation in (A) TPXO7.2 (an update to the model described by Egbert et al., 1994), a barotropic tide model constrained by satellite altimetry, and (B) 32-layer HYCOM simulations forced concurrently by atmospheric fields and the tidal gravitational potential. Lines of constant phase lag plotted every 45 • are overlaid in white. Reproduced from Fig. 2 of Shriver et al. (2012), ©American Geophysical Union, Wiley Online Library, used with permission. surface tidal currents, and Timko et al. (2017), who examined the sensitivity of resolved internal tide generation to small-scale topographic features that are not well resolved in global bathymetric datasets. Coastal and regional models often include fine grid spacing (high resolution) along with combined atmospheric and tidal forcing. As an example of such work from one group, Kerry et al. (2014aKerry et al. ( , 2014bKerry et al. ( , 2016 examined the impacts of subtidal flows on internal tide generation, propagation, mixing, and non-stationarity (incoherence) in a regional model of the Philippine Sea that resolved barotropic tides, baroclinic tides, and the time-evolving atmospherically forced general circulation and associated mesoscale eddy field.
Finally, the inclusion of tidal forcing in high-resolution global ocean general circulation models is a more recent development. The earliest simulations that included both atmospheric and tidal forcing in global general circulation models (Thomas et al., 2001;Schiller and Fiedler, 2007;Müller et al., 2010) had grid spacings of order 1 • , which is not fine enough to simulate mesoscale eddies or internal tides and gravity waves. The first documented global simulations with high resolution and simultaneous tidal and atmospheric forcing were performed with HYCOM . The STORMTIDE simulations Li et al., 2015) of the Max Planck Institute Ocean Model (MPI-OM) and simulations of the National Oceanic and Atmospheric Administration (NOAA) Generalized Ocean Layered Model (GOLD; e.g., Waterhouse et al., 2014) followed shortly afterwards. Tidal elevations and currents in the HYCOM simulations are well-vetted; they have been compared to various in-situ and remotely sensed observations including satellite altimetry, moored temperature and velocity time series, tide gauges, and dissipation estimates inferred from the ARGO float array, among others. Global HYCOM tide simulations run at the US Naval Research Laboratory have been run with horizontal grid spacings of 1/12.5 • and 1/25 • , and we will refer to these simulations as "HYCOM12" and "HYCOM25" below. HYCOM employs a hybrid Arbitrary Lagrangian Eulerian vertical coordinate (Bleck, 2002;Griffies et al., 2020) which smoothly transitions from z-levels in the open-ocean mixed layer to isopycnal (density) coordinates in the open-ocean interior, where mixing is weak, to terrain-following coordinates in shallow water regions. The number of vertical levels in the HYCOM simulations reported on here varies from 32 to 41. Simulations of the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al., 1997), performed on NASA supercomputers with 90 z-levels and horizontal grid spacings of 1/12 • , 1/24 • , and 1/48 • (referred to below as MITgcm12, MITgcm24, and MITgcm48, respectively) were first reported on in Rocha et al. (2016a,b), and have subsequently been used in a number of studies, especially studies with a focus on the upcoming SWOT mission (e.g., Savage et al., 2017b;Qiu et al., 2018;Torres et al., 2018Torres et al., , 2019Wang et al., 2018Wang et al., , 2019Klein et al., 2019;Chereskin et al., 2019;Luecke et al., 2020;Yu et al., 2021). A separate global simulation of MITgcm with simultaneous tidal and atmospheric forcing has been performed in a collaboration between two Chinese institutions (Fu et al., 2021).
Other large modeling groups have followed suit in the inclusion of tidal forcing in high-resolution simulations of the basin-and global-scale general circulation. Tides have recently been inserted into global simulations of the NOAA Modular Ocean Model version 6 (MOM6) model, the Department of Energy (DOE) Model for Prediction Across Scales (MPAS), and the Nucleus for European Modelling of the Ocean (NEMO), into 1/50 • North Atlantic simulations of HYCOM, into 1/60 • North Atlantic simulations of NEMO, and into high-resolution North Pacific simulations of the Regional Modeling Ocean System (ROMS). Regional simulations, forced at their lateral boundaries by global internal wave models, are being performed by a University of Toronto group (using the MITgcm), and the University of California Los Angeles (UCLA) ROMS group. Initial results and plans for high-resolution global "wind plus tides" simulations in Australia (Callum Shakespeare, Australian National University, personal communication, 2021), the United Kingdom (Helene Hewitt, United Kingdom Meteorological Office, personal communication, 2016), Italy (Nadia Pinardi, University of Bologna, personal communication, 2019), and elsewhere are underway. It is clear that this class of models has wide appeal and will continue to grow.

Atmospheric forcing
The ocean exchanges momentum, heat, fresh water, and gases with the atmosphere (see Csanady, 2001;Josey et al., 2013, amongst many). Together with background mixing, much of which is due to IGW breaking, wind stress and buoyancy forcing set up the large-scale stratification and circulation of the ocean (e.g., Pedlosky, 1996;Vallis, 2017). Three other aspects of atmospheric forcing merit a brief mention here, due to their impacts on high-frequency oceanic motions. Nearinertial waves are generated by high-frequency variations in wind stress. Models with hourly atmospheric fields have been shown to contain NIW kinetic energy up to three times higher than models forced by 6-hourly atmospheric fields (e.g., Rimac et al., 2013;Flexas et al., 2019). High-frequency oceanic motions are also driven by atmospheric pressure loading (e.g., Ponte, 1994;Tierney et al., 2000;Stammer et al., 2000;Carrère & Lyard, 2003), an effect that is not always activated in oceanic general circulation simulations. Most of the power in frequency spectra of atmospheric pressure is broadband, but there are peaks at 12 and 24 h due to the atmospheric tides, which arise from the diurnal cycle of solar heating (e.g., Chapman & Lindzen, 1970). The small S 1 tide in the ocean has a period of 24 h and is predominantly driven by the atmospheric S 1 tide . The 12-hour S 2 tide, the second largest tidal constituent in the ocean, is predominantly driven by the solar gravitational field. However, in a global average, about 15% of the S 2 barotropic tide SSH is forced by the pressure loading of the atmospheric S 2 tide (e.g., Arbic, 2005 and references therein).
The earlier HYCOM tide simulations presented here were forced by the United States Navy Operational Global Atmospheric Prediction System (NOGAPS; Rosmond et al., 2002), and the more recent HYCOM simulations are forced by the United States Navy Global Environmental Model (NAVGEM; Hogan et al., 2014). The frequency of NAVGEM output has varied over time, such that some recent HYCOM simulations are forced by hourly outputs and some by 3-hourly outputs. The MITgcm simulations presented here were forced by six-hourly output from the 0.14 • European Centre for Medium-Range Weather Forecasts (ECMWF) operational atmospheric reanalysis, converted to surface fluxes using bulk formulae from Large and Yeager (2004). Both the HYCOM and MITgcm groups are performing coupled high-resolution ocean-atmosphere runs with frequent coupling (approximately once per hour, or more frequently). In addition, some newer HYCOM results include data assimilation acting on both oceanic and atmospheric components. However, none of the HYCOM and MITgcm results shown here are from the coupled ocean-atmosphere runs, which are relatively new and are still under intense investigation.

Bathymetry
Bathymetry is a topic of continuing interest within the tide modeling community. The accuracy of tides in both regional and global tide models depends on the accuracy of the underlying bathymetric dataset (Florent Lyard, Centre National de la Recherche Scientifique, personal communications over many years; Joannes Westerink, personal communication, University of Notre Dame, 2020). Most global bathymetries use updated versions of the Smith & Sandwell (1997) product (see, for instance, Becker et al. 2009;Tozer et al. 2019; and the SRTM products at https://topex.ucsd.edu/). These products are based upon acoustic soundings of bathymetry where available and satellite altimetry measurements in other locations. Altimetry can map bathymetry because positive and negative mass anomalies, such as undersea mountain chains and trenches, perturb sea level. Altimetry cannot resolve features smaller than about π times the ocean depth (Smith & Sandwell, 1997), which is about 10-20 km in the pelagic ocean. Altimetric sensing of bathymetry works less well in continental shelves where there is significant sediment cover. The most accurate bathymetric measurements in shelf areas are often held closely by local governments because of their military significance. Abyssal hills, common small-scale features on the ocean floor (Goff & Jordan, 1988, 1989, are not well resolved in global bathymetric databases except in locations where acoustic soundings are available (about 10% of the seafloor, and most often in coastal regions; Charette & Smith, 2010;Wessel & Chandler, 2011). As a remedy, Goff & Arbic (2010) and Goff (2010) put forth a statistical description of abyssal hills that can be used in either spectral or physical space domains to "make up" for the abyssal hills that are lacking in global bathymetry products. The Goff & Arbic (2010) and Goff (2010) products have been used in several studies. For instance, Melet et al. (2013) combined the statistical abyssal hill results with the tidal conversion scheme of Nycander (2005) to estimate the conversion resulting from the "missing" abyssal hills, and Timko et al. (2017) followed with a quantification of the extra resolved internal tide conversion in a global internal tide model. See Sandwell et al. (2022) for an updated bathymetric product that also incorporates synthetic smallscale bathymetry.

Governing equations for tides
To illustrate the mechanics of tidal insertion into general circulation models, we provide the governing equations for tides in the simplest case, namely, one-layer shallow water motion (Gill, 1982). If we include quadratic bottom drag and eddy viscosity, and assume a tensor form of the wave drag, as in Garner (2005) and Nycander (2005), the governing momentum equation is and the governing mass conservation equation is where t is time, g is gravitational acceleration, f is the Coriolis parameter, k is a unit vector in the vertical direction, H is the resting water depth, η is the perturbation tidal elevation, η EQ is the equilibrium tidal potential, η SAL is the self-attraction and loading term, u → = (u,v) is the horizontal velocity vector, with u and v representing the zonal (east-west) and meridional (north-south) components, respectively, K H is the horizontal eddy viscosity, c d is the quadratic drag coefficient, T is the topographic internal wave drag tensor, and ρ 0 ≈ 1035 kg m − 3 is the average density of seawater. The equilibrium tide, self-attraction and loading term, and drag terms are discussed in more detail below. Arbic et al. (2004a) presents one-layer governing equations that are equivalent to (1) and (2), except that they are written in a slightly different (flux-divergent) form. The governing equations for two-and multi-layer shallow water tide-only models are provided in Arbic et al. (2004a) and Simmons et al. (2004b), respectively. The "Laplace tidal equations" (Laplace 1775(Laplace , 1776) are shallow-water equations similar to (1) and (2) but without the nonlinear advective terms and the damping terms. Thus one might argue that the very first ocean model in history was a tide model. As Laplace understood, the actual tidal elevation field in the ocean (represented by η) differs from the much simpler equilibrium tide η EQ for several reasons; friction, the obstruction of continents, the Earth's rotation, the presence of the self-attraction and loading term, and the fact that the shallowwater phase speed in the open ocean-about 200 m s − 1 -is not fast enough to keep up with the Moon overhead. The shallow-water equations (1) and (2) model the dynamical response of the ocean to the equilibrium tidal forcing. The equilibrium tide, or astronomical tidal potential, arises from the difference in the gravitational tidal potential from a distant body across a body of finite size (Newton, 1687;Pugh, 1987;Cartwright, 1999;Pugh & Woodworth, 2014), and varies as the inverse third power of the distance between the Earth and the disturbing body (Moon or Sun). Thus, the tidal potential raised by the Sun on Earth is about 46% that raised by the Moon on Earth, despite the fact that the Sun is about 27,000,000 times more massive than the Moon. An illustrative diagram, and a more detailed explanation, for the equilibrium tide, can be found in Arbic et al. (2018), amongst many other sources. Here we will only discuss tides arising from the degree two tidal potential; tides arising from the degree three tidal potential are detectable but much smaller (Cartwright, 1975;Woodworth, 2019;Ray, 2020).
The equilibrium tide consists of two bulges under which the Earth rotates. Hence the period of the principal lunar semi-diurnal tide, M 2 , is one-half of a lunar day (a day measured against the Moon), or about 12.42 h. The period of the principal solar semi-diurnal tide, S 2 , is onehalf of a solar day (a day measured against the Sun), or 12 h. The unequal periods of M 2 and S 2 yield a classic "beat" pattern in the tides. Spring tides occur twice a month, during New Moon and Full Moon, when the Earth, Moon, and Sun lie in a straight line, an occurrence known as "syzygy". Spring tides feature a maximal tidal range-the difference between high and low tides. In contrast, during "neap tide", the line between the Earth and Moon is at right angles to the line between the Earth and Sun, and the difference between high and low tides is at a minimum. Tides of diurnal period arise because the canonical equilibrium tidal bulges are tilted with respect to the Earth's equator.
The solar equilibrium tidal bulge is so tilted because of Earth's obliquity-the tilt of the Earth's spin axis with respect to the plane of the ecliptic. The lunar equilibrium tidal bulge is tilted due to obliquity and also to lunar inclination-the inclination of the plane of the lunar orbit to the ecliptic by about 5.145 • . Due to the tilt of the equilibrium tidal bulges with respect to Earth's equator, the two high equilibrium tides every day are unequal, and the two low equilibrium tides are also unequal (amongst many references, see, for instance, Figure 13.11 in Arbic et al., 2018). This "diurnal inequality" gives rise to diurnal tides. The largest diurnal tide, K 1 , has a period exactly equal to a sidereal day, meaning a day measured against the distant stars (86164 s, or 3 min and 56 s less than a solar day), and includes contributions from both the lunar and solar gravitational field. The second and third largest diurnal tides, O 1 and P 1 , are forced exclusively by the lunar and solar gravitational fields, respectively. Complete treatments of the equilibrium tides (Doodson, 1921;Cartwright & Tayler, 1971;Roosbeek, 1996) are computed in terms of spherical harmonic expansions, where the semi-diurnal and diurnal terms respectively represent the "22" and "21" terms of the degree 2 spherical harmonic (see, for instance, Williams & Boggs, 2016). The spherical harmonic expansions also yield "20" terms, the long-period tides. Doodson's treatment yielded more than 400 tidal lines, and Roosbeek's yields an even greater number of tidal lines.
For a single semi-diurnal tidal line, the equilibrium tide η EQ below, is given by where t is time, t ref is a reference time, A and ω represent amplitude and frequency, φ is latitude, λ is longitude, and χ (t ref ) is the astronomical argument referenced to t ref (e.g., Schwiderski, 1980;Pugh & Woodworth, 2014). The Love numbers h 2 and k 2 account respectively for the solid earth "body-tide" deformation response to the astronomical forcing, and the perturbation to the gravitational potential resulting from the redistributed mass of the solid Earth. The nodal factors f(t ref ) and ν(t ref ) account respectively for slow modulations in the amplitude and phase lag of a major tidal line due to the action of smaller lines of nearby frequency during the 18.6 year nodal cycle (Pugh & Woodworth, 2014).
[Note that in this paper, as in many sources, the symbol f is used for both Coriolis parameter (equation (1), and elsewhere) and nodal factor for amplitude Eqs. (3, 4, 5 and 7); the reader can determine through context , and ν(t ref ) are all frequencydependent-they differ from one tidal line to another. The values of h 2 and k 2 are constant within the semi-diurnal band but vary within the diurnal band due to free-core nutation resonance within the solid Earth (Wahr and Sasao, 1981). The nodal modulations are ignored in some applications. For a single diurnal tidal line, the equilibrium tide is given by and for a single long-period tidal line, the equilibrium tide is given by The long-period constituents have relatively small amplitudes and will not be further discussed here. If one wants to include many tidal lines in a model, one simply includes equilibrium tidal terms like those in (3), (4), and (5), one term for each tidal line. Values of A, ω, period (2 π / ω), and the combination 1 + k 2 -h 2 , for the four largest semi-diurnal lines (M 2 , S 2 , N 2 , K 2 ), the four largest diurnal lines (K 1 , O 1 , P 1 , Q 1 ) and the two largest long-period lines (M m , M f ) can be found in Cartwright & Tayler (1971), Table 13.1 of Arbic et al. (2018), and many other sources. The combination 1 + k 2 -h 2, sometimes called a diminishing factor , has a value of about 0.7 for most tidal constituents, meaning that the elasticity of the solid Earth reduces the amplitude of the ocean tides by about 30% from the amplitudes that would be seen on an Earth that was completely rigid. Early HYCOM tide simulations (e.g., Shriver et al., 2012) included the four largest diurnal lines and the four largest semidiurnal lines. More recent HYCOM simulations include only the five largest tidal lines-M 2 , S 2 , N 2 , K 1 , and O 1 . The MITgcm simulations discussed here employ the full luni-solar tidal potential (Weis et al., 2008), which as discussed earlier includes hundreds of lines.
Because the solid Earth is elastic, it deforms under the load of the ocean tides as well as the direct astronomical forcing. The self-gravity of the load-deformed solid Earth, and of the ocean tides themselves, both perturb the gravitational potential. The self-attraction and loading (SAL) term (Hendershott, 1972;Ray, 1998) represents the sum of these effects, and is given by where n is a spherical harmonic index, η n is the nth spherical harmonic of η, and ρ earth ≈ 5518 kg m − 3 is average solid Earth density. The load numbers h' n and k' n (Munk and MacDonald, 1960), account respectively for solid-earth yielding and the resulting perturbation potential. Because the self-attraction and loading term involves spherical harmonics, which are computationally expensive, special procedures for handling the selfattraction and loading term have been developed over the years. In the "scalar approximation" (Accad & Pekeris, 1978), which is used in older HYCOM tidal simulations and in some of the other model results shown in this paper, η SAL is given by a simple coefficient (of order ~ 0.1) multiplied by η. The self-attraction and loading term can also be computed iteratively (e.g., Egbert et al., 2004;Arbic et al., 2004a). This iterative approach has been used in some HYCOM tide runs, but in newer HYCOM tide runs, SAL fields are read in from the TPXO (Topex/ Poseidon Cross-Over) satellite-altimeter-constrained barotropic tide model, which has continued to evolve over time (Egbert et al., 1994;Egbert & Erofeeva, 2002).
In a promising new approach for the treatment of SAL, outlined in Schindelegger et al. (2018), the spherical harmonics are computed directly inline (within the model, as it is running), taking advantage of new, computationally efficient spherical harmonic packages (see also Stepanov and Hughes, 2004;Kuhlmann et al., 2011;Vinogradova et al., 2015). The inline approach is more accurate than the scalar approximation, which is strictly incorrect (Ray, 1998) because it does not account for the horizontal scale dependence in (6). The inline approach is less tedious than the iterative approach, because the latter involves performing multiple simulations and the former involves performing only one simulation. The inline approach to SAL can be used for simulations of paleotides or tides in a future world with a different climate; in contrast, the approach of reading in SAL from established models of the present-day tides cannot be used for such applications. In addition, whereas the iterative approach is only easily applicable in the case of periodic motions such as tides, the inline approach applies the selfattraction and loading term to all motions with a mass signal, whether they are periodic or not, as Hendershott (1972) anticipated should be done.
Through analysis of tidal dissipation in the Irish Sea, Taylor (1919) determined that bottom boundary layer drag should be modeled as quadratic in the velocity, as seen in (1). Modern ocean models still use his formulation, with values of the drag coefficient c d of order 0.001-0.003. The energy dissipation equation, which is obtained by taking the dot product of (1) with u → , thus yields a bottom boundary layer dissipation term proportional to the cube of the velocity. Because barotropic tidal velocities are much larger in coastal regions (up to 1-2 m s − 1 ) than in the open ocean, where they are of order 1-2 cm s − 1 , the boundary layer dissipation term is very large in shallow seas and negligible in the open ocean (Munk, 1968(Munk, , 1997Jayne & St. Laurent, 2001). Parameterized topographic wave drag schemes have been employed in several barotropic tide models (e.g., Jayne & St. Laurent, 2001;Egbert et al., 2004;Arbic et al., 2004a;Griffiths & Peltier, 2008, 2009Green, 2010;Buijsman et al., 2015). The schemes in Egbert et al. (2004) and Arbic et al. (2004a) employ tensors, as shown schematically in (1). The tensors vary geographically according to topographic roughness and the local Brunt-Väisälä frequency N. Recent HYCOM simulations employ the simpler (non-tensor) Jayne & St. Laurent (2001) scheme.
To our knowledge, the HYCOM tidal simulations are the only published tidally forced general circulation simulations for which parameterized topographic wave drag has been inserted. Several major complications arise when inserting wave drag into circulation models. First, insertion of a wave drag may seem like "double counting" in internal wave models, which after all resolve the generation of some of the internal tides and gravity waves. In several papers (Arbic et al., 2004a, 2012, Ansong et al., 2015Buijsman et al., 2016Buijsman et al., , 2020, we have argued that parameterized wave drag is still needed in global internal wave models. Such models do not resolve the breaking of smallscale internal waves, and this energy loss must be otherwise parameterized. As will be discussed later in this paper, Ansong et al. (2015) demonstrated that the internal tides are too large in HYCOM simulations that do not include parameterized wave drag. Second, one must decide where to apply the wave drag. Thus far, we have applied it to the nearbottom flow, defined in our HYCOM tide papers as the flow averaged over the bottom 500 m. The near-bottom flow is a combination of the barotropic and low-vertical-mode baroclinic flows. Thus, application of the wave drag to the near-bottom flow implies that the drag acts on both barotropic and baroclinic flows. Third, the dissipation patterns (Egbert & Ray, 2003) and wave drag strengths for semi-diurnal and diurnal tides differ, such that it is not straightforward to optimally tune a wave drag for both semi-diurnal and diurnal tides (Skiba et al., 2013). In HYCOM we tune the wave drag for semi-diurnal tides, and therefore the wave drag is not well tuned for diurnal tides. Fourth, one must decide whether to apply the wave drag to all resolved motions in the model, including mesoscale eddy motions. Because the wave drag acting on lowfrequency and tidal motions is not equal (Bell, 1975), we have thus far decided to apply the wave drag in HYCOM only to tidal motions. Separation of tidal and non-tidal flows is achieved through filtering the near-bottom flow, as described in ; see Shriver et al. (2012) for an updated description of the wave drag procedure in HYCOM.

Data assimilation
Global operational US Navy HYCOM forecasts assimilate nadir altimetry, sea surface temperature, and available in-situ data. We are just beginning to analyze new simulations of HYCOM that simultaneously include data assimilation and tides. Preliminary analysis demonstrates that the internal tides in assimilative HYCOM runs are more accurate than in non-assimilative HYCOM runs, in part because assimilation improves the agreement of HYCOM stratification with observations (Luecke et al., 2017). Since 2016, HYCOM tidal simulations have employed an Augmented State Ensemble Kalman Filter (ASEnKF), taken from the assimilation literature, to improve the accuracy of barotropic tides. The ASEnKF applies only to the barotropic tidal motions in HYCOM. Thus, the newer HYCOM simulations that do not assimilate nadir altimeter data to forecast mesoscale eddies are still referred to here as "non-assimilative". All of the HYCOM results shown in this paper are from non-assimilative HYCOM simulations. The HYCOM results shown here are based upon both older simulations that do not employ the ASEnKF, and newer simulations that do.

Tidal analysis
Because tidal motions are periodic, tidal variables V (e.g., tidal elevations, and the components of horizontal tidal velocities) can be written in terms of an amplitude and phase lag, viz.
This decomposition is known as a tidal harmonic analysis. The amplitude and phase lag are solved for with a least-squares fitting procedure in the tidal harmonic analysis routines of Foreman (1977Foreman ( , 2004, Pawlowicz et al. (2002), and Foreman et al. (2009). An alternative approach is based on the work of Doodson (1921) which describes tidal species in terms of lunar speeds rather than solar ones. Differences between the two approaches are described in, for instance, Pugh & Woodworth (2014). The different histories of these approaches have resulted in important differences in tidal analysis software packages, and it is always good practice for consistency to use the same package for both tidal analysis and prediction (tidal prediction is the usage of the derived harmonic constants to predict tidal motions).
The squared discrepancy D 2 at a particular location, between tidal elevations η MODEL from a model and η OBS from observations, is defined as where time-averaging is denoted by < >. In terms of tidal amplitudes and phase lags, D 2 is given by where A MODEL and A OBS are respectively the model and observation amplitudes, and φ MODEL and φ OBS are respectively the model and observation phase lags. See Shriver et al. (2012) and Arbic et al. (2018) for an alternative decomposition of D 2 into an "amplitude error" due to differences in amplitudes and an "amplitude-weighted phase error" due to differences in the phase lags. An area-weighted RMS error is given by

Barotropic tide results
In this section, we describe some representative barotropic tide results from global tidally forced high-resolution general circulation models.

Barotropic tides in tidally forced general circulation models
We focus on the barotropic tidal errors in HYCOM, the most thoroughly vetted tidally forced global general circulation model. The surface tidal elevations in older HYCOM tidal simulations have been compared to tide gauges Stammer et al., 2014) and the surface elevations in both older and newer HYCOM simulations are routinely compared to TPXO. Visual comparison of M 2 surface elevation amplitudes and phase lags in the older HYCOM simulations with those in the TPXO altimeter-constrained barotropic tide model (Fig. 2) indicates that HYCOM is capturing many of the important features. However, differences between HYCOM and TPXO are also discernible, in, for instance, the amplitude patterns around the southern and southeastern coasts of Africa, and the amphidrome pattern at about 160 • W, 60 • S. The global open-ocean RMS M 2 elevation errors between early HYCOM simulations and tide gauges or altimetry, computed from (10) over waters deeper than 1000 m and latitudes equatorward of the TOPEX/ JASON turning latitude (66 • ), are at the level of about 7-8 cm, comparable to those in earlier forward (non-data-assimilative) barotropic tide models (e.g., Jayne & St. Laurent, 2001;Egbert et al., 2004;Arbic et al., 2004a). A summary comparison between hydrodynamical models (HYCOM, STORMTIDE, and several forward barotropic tide models--Stammer et al., 2014; their Table 12) and both tide gauges and satellite altimetry indicates that STORMTIDE errors are comparable to, but slightly larger than, those in HYCOM. Barotropic tidal currents in HYCOM, STORMTIDE, several forward barotropic tide models, and four altimeter-constrained barotropic tide models, were compared to barotropic tidal currents inferred from historical current meters in Section 8.1 of Stammer et al. (2014). The discrepancies between HYCOM, STORMTIDE, and the forward barotropic tide models and the current meters are comparable to each other and are larger than the discrepancies between the data-assimilative models and current meters. However, the level of the discrepancies between the purely hydrodynamical forward tide models and current meters (about 1 cm s − 1 , compared to an RMS current meter signal of 1.5 cm s − 1 ) is only slightly higher than the ~ 0.8-0.9 cm s − 1 discrepancy between the dataconstrained models and current meters. In contrast, tide gauge comparisons  demonstrate that data-constrained models have much more accurate elevations than the forward hydrodynamical tide models. Data-constrained tide models perform better with elevations, which are constrained by satellite alltimetry, than they do with tidal currents. An additional complexity in modeling tidal currents, and comparing them to observations, is the strong vertical structure of currents, due to baroclinic tidal motions. This structure is not  Fig. 6), inclusion of an iterated self-attraction and loading term, rather than the simple scalar approximation used in earlier HYCOM simulations, and a more accurate treatment of the water depths underneath floating ice shelves in Antarctica-a well-known factor in accurate modeling of global barotropic tides-yield substantially lower M 2 tidal elevation errors with respect to TPXO. Further improvements are obtained with an Augmented State Ensemble Kalman Filter (ASEnKF), in which sets of perturbations with horizontal scales comparable to those of open-ocean barotropic tides are introduced one-by-one to the tidal forcing. The ASEnKF machinery examines the spatial patterns of error introduced by the various perturbations, and generates a linear combination of perturbations that minimizes the errors. The globally averaged RMS M 2 elevation error in HYCOM with respect to TPXO, again computed from (10) in the open ocean equatorward of 66 • , is reduced from 4.4 cm in a best-case non-assimilative solution, to 2.6 cm in the ASEnKF solution. This improvement, while substantial, still does not bring the barotropic tidal errors in HYCOM down to the levels seen in data-constrained barotropic tide models.
We are seeking further improvements in HYCOM barotropic tidal accuracy, based upon the "back-effect" of coastal tides upon open-ocean tides and the known importance of bathymetry. The traditional conceptual model of coastal tides consists of a coastal region that receives information from the open-ocean, at a lateral boundary, but does not transmit information back into the open-ocean. A series of papers (Arbic et al., 2007;Arbic et al., 2009;Arbic & Garrett, 2010)  and with results from a simple coupled oscillator model in Arbic and Garrett (2010). The most important coastal region, whose removal has the largest back effect, is the Hudson Strait/Hudson Bay system in eastern Canada. Removal of this system increases the globally averaged M 2 amplitude by about 13%. Because coastal regions impact open-ocean tides, more accurate modeling of coastal tides may improve the modeling of open-ocean tides as well. Simulations with two-way nesting, and higher horizontal resolution, in regions of large coastal tides, do indeed see an improved accuracy in the modeled open-ocean tides (Jeon et al., 2019). However, the improvements come with a substantial added computational cost. In another effort, we are performing additional ASEnKF perturbations in regions of large coastal tides. The perturbations in such regions have horizontal scales similar to those of coastal tides. Finally, we are attempting to improve the underlying topography, which has been shown to improve the accuracy of global barotropic tide models (e.g., Lyard et al., 2021).
Barotropic tides embedded within a general circulation model have a measurable seasonal cycle (Fig. 3). Regions with large seasonal changes include the Labrador Sea, adjacent to the Hudson Bay/Hudson Strait system. Observations of the Hudson Bay/Hudson Strait system display a measurable seasonal cycle, which has been attributed to seasonal sea-ice cover (St-Laurent et al., 2008). Tides drag along the underside of sea ice, effectively doubling the quadratic friction in winter compared with summer. More details on seasonal tidal variations, using observations and regional models, can be found in Müller et al. (2014). Global tidally forced general circulation simulations such as the STORMTIDE results used in Müller et al. (2014) provide a global look at seasonal barotropic tide variability.
Global HYCOM reproduces known tidal mixing fronts in several shallow seas around the global ocean with a high degree of skill (Timko et al., 2019). The addition of tides into HYCOM impacts summertime sea surface temperatures in regions of strong tides by as much as 2 • C (Fig. 4).

Barotropic tides in a changing climate
Analysis of tide gauge observations indicates that tidal elevations have been changing on centennial timescales (e.g., Flick et al., 2003;Ray, 2006;Jay, 2009;Woodworth, 2010;Müller et al., 2011;Devlin et al., 2014;Feng et al., 2015;Talke et al., 2018;. See the recent review by Haigh et al. (2020) for more detail and more references on this evolving topic. Both positive and negative trends in tidal elevations are visible (Fig. 5). An early attempt (Müller et al., 2011) to model the impact of sea level rise on the centennial-scale secular changes in tides was less successful than a more recent attempt (Schindelegger et al., 2018). The latter study employed higher model resolution than the former study, and used an inline formulation of the self-attraction and loading, rather than an iterative procedure. Iterations are not only tedious, but produce "iteration jitter" from one iteration to the next, and this iteration jitter is comparable to the small secular signals. The Müller et al. (2011) andSchindelegger et al. (2018) studies employed barotropic tide models, with water column thicknesses adjusted using estimates of centennial-scale sea-level rise and glacial isostatic adjustment. The advent of tidally forced general circulation models discussed here represents a potential opportunity to study centennial-scale secular changes in tides within a more complete earth system model. Indeed, DOE and NOAA modeling groups plan to run their tidally forced general circulation simulations over centennial time scales, in order to examine secular tide changes.
The increase in open-ocean tides with the removal of shelves, predicted by the analytical models in Arbic et al. (2009) and Arbic & Garrett (2010), are consistent with the finding that the modeled deep-ocean tides of the most recent ice ages, when sea level was lower and many regions of large present-day tidal dissipation on shelves did not exist, were substantially larger (e.g., Egbert et al., 2004;Arbic et al., 2004bArbic et al., , 2008Griffiths & Peltier, 2008, 2009Green, 2010;Wilmes & Green, 2014). Tides change even more dramatically (Bills & Ray, 1999;Green & Huber, 2013;Green et al., 2017;Daher et al., 2021) over geological time scales (millions to billions of years), in part due to changes in Earth rotation rate, which is continually slowing down because of tidal friction (e.g., Williams & Boggs, 2016) and in part because of tectonically driven changes in continental configurations.

Internal tide and gravity wave results
In this section, we describe some representative results on internal tides and gravity waves from global high-resolution general circulation models that include tidal forcing.

Comparison of modeled internal tides with altimetry
An early model-altimeter comparison of stationary internal tide SSH perturbations (Fig. 6) is taken from Merrifield et al. (2001), who ran a regional model of the internal tides around Hawai'i. The regional model was forced at its lateral boundaries (Merrifield and Holloway, 2002) by the data-assimilative TPXO barotropic tide model. The TPXO boundary forcing yields accurate barotropic tides within the regional model and hence contributes to the "wiggle matching" of internal tide perturbations displayed in Fig. 6. In global internal wave models, the barotropic tides must be generated within the model, without benefit of boundary conditions from accurate altimeter-constrained models, thus making it more difficult to "wiggle-match" altimeter data with high accuracy.
The stationary internal tide SSH perturbations in HYCOM have been compared with along-track altimeter results in several papers Ansong et al., 2015;Buijsman et al., 2020). As demonstrated in Ansong et al. (2015), internal tides in HYCOM simulations with wave drag applied to the near-bottom flow (Fig. 7b,c) lie closer to altimeter results (Fig. 7a) than do internal tides in simulations with wave drag applied only to the barotropic flow (Fig. 7d) or in simulations without wave drag (Fig. 7e,f). An in-preparation paper by Joseph Ansong demonstrates that internal tides in other tidally forced generally circulation models (MITgcm, MOM6, and NEMO), none of which include wave drag Fig. 6. The M 2 baroclinic sea surface elevation amplitude around Hawai'i from a regional model (dark lines) and as estimated from TOPEX-Poseidon satellite altimeter data (Ray & Mitchum, 1996). The amplitudes have been detrended along each track to remove variations with long spatial scales that are associated with the barotropic tide. Numbers within the figure represent altimeter track numbers. Reproduced from Fig. 3 of Merrifield et al. (2001), ©American Geophysical Union, Wiley Online Library, used with permission.  Ansong et al. (2015), which can be consulted for more details of the E051-E059 simulations, ©American Geophysical Union, Wiley Online Library, used with permission.
to date, are also too large relative to observations. The importance of damping on the amplitudes of internal tides has also been demonstrated in reduced physics models (Kelly et al., 2021). Internal tides in STORMTIDE have also been compared to altimetry, in a manner similar to what is seen in Fig. 7-see Fig. 1 of Müller et al. (2012). In STORM-TIDE as in HYCOM and the models examined in the Ansong et al. inpreparation paper, internal tide generation hotspots identified by Shriver et al. (2012)-Hawai'i, east of Philippines, tropical South Pacific, tropical Southwest Pacific, and Madagascar-are all clearly visible.
Next we turn to the ability of global hydrodynamical internal wave models to match satellite altimetry observations of internal tides on a "wiggle-by-wiggle" level. Visual comparisons such as those in Fig. 7, and comparisons of amplitudes averaged over large areas, provide information about the ability of global internal wave models to match the general spatial patterns seen in altimetry. Matching the internal tides "wiggle-by-wiggle", as in the Fig. 6 regional model results but on a global scale, requires greater model skill. A first test of HYCOM's ability to wiggle match altimeter observations was performed in Carrère et al. (2021), who included non-assimilative HYCOM as an example hydrodynamical model in an inter-comparison of several empirical internal tide models, and one data-assimilative reduced-gravity model, with altimetry. The internal tide model inter-comparison in Carrère et al. (2021) is motivated by the need for internal tide correction models for existing nadir altimeter missions and for the upcoming SWOT wideswath altimeter mission. The Zaron (2019) empirical model is currently providing internal tide corrections for nadir altimeters. Carrère et al. (2021) found that non-assimilative HYCOM has some skill in reducing internal tide SSH variance from altimeter records, but not as much skill as the empirical and data-assimilative reduced-gravity models. A comparison of the internal tide amplitudes north of Hawai'i ( Fig. 8) demonstrates that HYCOM exhibits slightly larger amplitudes in this region, and some differences in its geographical patterns relative to data-constrained internal tide models, but also shares many similarities. Preliminary analysis of the internal tides in assimilative HYCOM solutions indicates that they lie closer to altimetry than the tides in nonassimilative HYCOM.
Internal tide non-stationarity is a topic of continuing interest in the community, including those using global internal tide models. Müller et al. (2012) demonstrated that seasonal changes are seen in the internal tides of the STORMTIDE model. The seasonal tidal changes to internal and barotropic tides shown in Müller et al. (2012) and Müller et al. (2014), respectively, represent a special case of tidal non-stationarity. HYCOM-based studies of non-stationary tides include Shriver et al. (2014), who examined differences between internal tide estimates computed over different time windows, Buijsman et al. (2017), who demonstrated that internal tide energy fluxes emanating from generation hotspots in the North and South Pacific pass through the equator, but become non-stationary due to seasonal changes in stratification and interactions with equatorial jets, and Duda et al. (2018), who documented strong "bending" of internal tides in the Gulf Stream. Ansong et al. (2015), who focused on the sensitivity of the stationary internal tides to parameterized topographic internal wave drag, found that the stationary internal tide amplitudes were also sensitive to the duration of the model records. We interpret this as being due to the increased ability of a tidal harmonic analysis to discriminate between a stationary internal tide and other motions (e.g., mesoscale eddies) with longer record durations. The sensitivity of stationary tide estimates to model record duration is another manifestation of internal tide non-stationarity in HYCOM.
Internal tide non-stationarity can also be estimated via the residual variance in the semi-diurnal band after the stationary component has been removed by standard tidal harmonic analysis-see, for instance, Figs. 14 and 15 in Savage et al. (2017a). Comprehensive results from this residual approach are displayed in Fig. 9, taken from Nelson et al. (2019), which displays the semi-diurnal non-stationary variance fraction (SNVF) -the ratio of residual non-stationary semi-diurnal variance to total (stationary plus non-stationary) semi-diurnal variance-in three different analyses of HYCOM and in an analysis of satellite altimetry by Zaron (2017). The uppermost subplot of Fig. 9 shows the SNVF computed from frequency spectra of five years of hourly HYCOM SSH output. The second and third subplots of Fig. 9 display the SNVF computed using wavenumber spectra computed from hourly and 9.9156-day samples of HYCOM output, respectively. Wavenumber spectra can be used to quantify non-stationarity because internal tides map to known wavenumbers which can be integrated over both before and after stationary tides have been subtracted from model time series. However, the wavenumber spectra technique requires the subtraction of a background mesoscale eddy signal, due to the fact that internal tides and mesoscale eddies have similar length scales. The wavenumber analysis is applied to 9.9156-day output as well as hourly output because 9.9156 days is the repeat period of the TOPEX/JASON nadir altimeter ground track. The 9.9156 day sampling time prevents us from using frequency spectra to compute the SNVF from altimetry. The bottom subplot of Fig. 9 displays the SNVF computed by Zaron (2017), using wavenumber spectra applied to nadir altimeter data. The third subplot of Fig. 9 thus displays HYCOM results computed using the same techniques and sampling used in the altimeter analysis. The four subplots of Fig. 9 display many similarities, including high values of SNVF (indicating near-total non-stationarity) in the equatorial Pacific, consistent with the results of Buijsman et al. (2017), and the equatorial Indian Ocean. The similarities of the HYCOM and altimeter analyses in Fig. 9 suggest that HYCOM has substantial skill in simulating the geographical variability of internal tide non-stationarity. Indeed,  used HYCOM in a study of the viability of predicting nonstationary internal tides in satellite altimeter observations.

Modeling the internal gravity wave continuum
Global internal wave models simulate a partially resolved IGW continuum spectrum, as first shown in Müller et al. (2015), who examined HYCOM results. The horizontal wavenumber-frequency (K-ω) spectrum, computed over a North Pacific region within global HYCOM, displays increasing variance in the higher-resolution HYCOM25 simulation relative to the lower-resolution HYCOM12 simulation (Fig. 10). The increased energy in HYCOM25 lies along the linear IGW dispersion curves of the first three vertical modes, defined by the equation where c e is the eigenspeed of the vertical mode in question. Eigenspeeds are computed from the IGW Sturm-Liouville problem, applied to vertical Fig. 9. Global maps of semidiurnal total SNVF (semidiurnal non-stationary variance fraction). The top three maps are computed from a 5 year, 1/12.5 • simulation of HYCOM using (top) f-space (frequency space) methodology applied to the hourly output, (second from top) kspace (wavenumber space) methodology applied to the hourly output, and (third from top) k-space methodology applied to the altimeter-sampled output. These are plotted alongside (bottom) the results from altimetry in Zaron (2017). In all subplots, results are only plotted where the total semidiurnal variance is greater than twice the error variance. Reproduced from Fig. 4 of Nelson et al. (2019), ©American Geophysical Union, Wiley Online Library, used with permission.
profiles taken from the model along the southern and northern boundaries of the North Pacific region. Eigenspeeds and f values vary from gridpoint to gridpoint. The maximum and minimum bounding curves for each vertical mode are displayed in Fig. 10. Müller et al. (2015) also performed a comparison of kinetic energy frequency spectra in models vs. mooring observations, and demonstrated that global high-resolution internal wave models capture some of the high-frequency (supertidal) kinetic energy variance seen in the observations, especially when model resolution is increased-in other words, the models simulate a partially resolved IGW continuum spectrum. The horizontal wavenumber-frequency spectra and other characteristics of the IGW continuum in global models has been examined in several papers since, including studies based on high-resolution simulations of the MITgcm (e.g., Rocha et al., 2016aRocha et al., , 2016bQiu et al., 2018;Torres et al., 2018Torres et al., , 2019Wang et al., 2018Wang et al., , 2019Klein et al., 2019;Chereskin et al., 2019), and studies that compare IGWs in HYCOM and MITgcm (Savage et al., 2017b;Luecke et al., 2020).

Comparison of modeled internal tides and gravity waves with in-situ observations
Internal tide and gravity wave fields in HYCOM (and, to a lesser extent, MITgcm) have been compared to multiple in-situ platforms. For instance, Timko et al. (2012Timko et al. ( , 2013 compared tidal currents in HYCOM at specific depth levels (meaning that both barotropic and baroclinic components were included) to currents inferred from historical current meter observations. Ansong et al. (2017) compared internal tide energy fluxes, and their variability, in HYCOM and 79 historical mooring sites. Savage et al. (2017a,b) compared the dynamic height frequency spectra in HYCOM and MITgcm to moored McLane profiler (Doherty et al., 1999) spectra. Rocha et al. (2016a) and Chereskin et al. (2019) compared wavenumber spectra in MITgcm48 and Acoustic Doppler Current Profiler (ADCP) data along ship tracks. Luecke et al. (2020) compared frequency spectra of temperature variance and kinetic energy in HYCOM12, HYCOM25, MITgcm12, MITgcm24, and MITgcm48 to spectra computed from ~ 2000 moored historical time-series records. Luecke et al. (2020) integrated over six different frequency bands ranging from low-frequency bands dominated by sub-tidal and mesoscale motions, through diurnal, semi-diurnal, and (in the case of kinetic energy only) near-inertial bands, to supertidal bands dominated by IGW continuum motions. The MITgcm48 energy levels compare most closely to observations in the supertidal band, due to the increased resolution of MITgcm48. However, HYCOM displays higher spatial correlations with observations over all frequency bands, presumably because HYCOM, as the backbone of an operational model, has undergone parameter tunings in order to place western boundary currents and other features in their correct locations. Example frequency spectra of temperature variance and kinetic energy, from the historical mooring data, the two HYCOM simulations, and the three MITgcm simulations (Fig. 11) demonstrate that MITgcm48 lies closest to the mooring observations at supertidal frequencies, and that in both models, increased resolution raises the supertidal energy levels. Tidal peaks and the low-frequency mesoscale eddy continuum, are also visible in Fig. 11, and near-inertial peaks are visible in the kinetic energy spectra displayed in that figure.
Measurements from the global drifter program offer another promising observational platform for comparison with high-frequency motions in global high-resolution internal wave models. The drifters are sampled frequently in time, such that hourly products (e.g., Elipot et al., 2016) can be produced. Drifters therefore offer a global product that can discriminate between low-and high-frequency motions, as with moored time series, and that provides global coverage, as with satellite altimetry. Model-drifter comparisons (Fig. 12) demonstrate that surface nearinertial kinetic energy in MITgcm is too weak while surface semi-diurnal kinetic energy is too strong. The former is likely due to the 6-hourly wind forcing of the MITgcm simulations, which is too infrequent to excite NIWs (Rimac et al., 2013;Flexas et al., 2019). The latter is likely primarily due to the lack of wave drag in MITgcm, which as shown in Ansong et al. (2015) yields tidal motions that are too strong. An inrevision paper by Arbic et al. indicates that HYCOM lies closer to the drifters in both the near-inertial and semi-diurnal bands, likely due to inclusion of more frequent wind forcing, and a parameterized topographic internal wave drag, respectively. The spectra E(ω) integrated over all horizontal wave numbers K, and E ITW (K), where "ITW" denotes "internal wave motions", integrated over frequencies larger than the inertial frequency at 29 • N (ω ≥ f 29N ), are shown as bottom and right-hand-side insets, respectively. In the E(ω) spectra, cyan lines indicate the inertial frequencies at the bounding latitudes 29 • N and 43 • N, while M 2 and M 4 frequencies are also indicated. Reproduced from Fig. 2 of Müller et al. (2015), ©American Geophysical Union, Wiley Online Library, used with permission.

Using global internal wave models in planning of field work and remote sensing missions
Global internal wave models have been used to plan for field work and remote sensing missions, especially SWOT. Richman et al. (2012) demonstrated, using HYCOM tide simulations, that the highwavenumber end of the SSH wavenumber spectrum, of greatest interest to SWOT, can be impacted by IGWs. IGWs flatten spectral slopes, from the − 5 or − 11/3 values expected from interior and surface quasigeostrophic turbulence, respectively (Le Traon et al., 2008;Sasaki & Klein, 2012), to slopes closer to the − 2 values expected from the Garrett-Munk IGW continuum spectrum. The slope of the SSH wavenumber spectrum in the mesoscale part of the spectrum (as defined by Xu & Fu, 2012) for the subtidal, tidal, and supertidal frequency bands, computed from HYCOM12, HYCOM25, MITgcm12, MITgcm24, and MITgcm48 over seven different regions of interest, is displayed in Savage et al. (2017b;their Fig. 10). [See Savage et al., 2017b for definitions of the subtidal, tidal, and supertidal bands and the seven regions of interest.]. Savage et al. (2017b) demonstrates that the slopes in the subtidal band generally lie close to the interior and surface quasi-geostrophic values mentioned above, while slopes in the tidal band are flatter and slopes in the supertidal band are also flatter, at least in the simulations with greater resolution. Fig. 13, also taken from Savage et al. (2017b), displays the SSH wavenumber spectra integrated over the subtidal, tidal, and supertidal frequency bands, computed from MITgcm48 over the seven regions. The tidal and supertidal bands often dominate the spectra at the high wavenumbers of greatest interest for SWOT, and the ratio of supertidal to total variance approaches unity over a wide band of high wavenumbers.
Following Richman et al. (2012), a number of papers have used global internal wave models to examine entanglement of internal waves and eddies in the high-spatial-resolution SSH fields that SWOT will measure (e.g., Rocha et al., 2016aRocha et al., , 2016bSavage et al., 2017aSavage et al., , 2017bQiu et al., 2018;Torres et al., 2018Torres et al., , 2019Wang et al., 2018Wang et al., , 2019Klein et al., 2019;Chereskin et al., 2019). The transition scale delineating the dominance of geostrophically balanced low-frequency eddying motions and unbalanced internal wave motions, computed from MITgcm48, varies from around 200-250 km near the equator to values around 50-100 km in mid-to high-latitudes (Fig. 14). Such studies further emphasize the need for accurate removal of high-frequency motions from nadir altimeter and SWOT observations before examination of lowfrequency motions takes place. Wang et al. (2018) used MITgcm48 simulations to design a ground truth "Cal/Val" (calibration and validation) field experiment for SWOT.

Boundary forcing of regional and coastal models
There is a growing realization that global internal wave models can provide important sources of information at the lateral boundaries of regional and coastal models (Reginald Beach, Office of Naval Research, Fig. 11. Sample frequency spectra of temperature variance (in the Eastern Pacific) and kinetic energy (Northeast Atlantic) for five global simulations (HYCOM12, HYCOM25, MITgcm12, MITgcm24, and MITgcm48), and moored historical observations. Instrument locations and depths are given in the subplots. The solid vertical lines show the diurnal (left) and semidiurnal (right) tidal frequencies, and the dashed vertical lines in the right-hand side subplots show the local Coriolis frequency, which is much more prominent in kinetic energy than in temperature variance. The approximate bounds of each frequency band discussed in Luecke et al. (2020) are delineated by boxes in the lower right subplot. Reproduced from Fig. 2 of Luecke et al. (2020), ©American Geophysical Union, Wiley Online Library, used with permission.
personal communication, 2020). Horizontal boundary conditions for regional models typically include only the barotropic tides (e.g., Merrifield & Holloway, 2002). However, several recent papers have shown that remotely generated internal waves can affect the flows inside regional models. For instance, Ponte & Cornuelle (2013) showed that baroclinic tidal kinetic energy in a regional model continues to increase as the domain size of the regional model increases. Kerry et al. (2013) showed that baroclinic tidal energy conversion in a regional model is affected by the presence of remotely generated internal tides. Suanda et al. (2017Suanda et al. ( , 2018 demonstrated that the addition of remote internal tides reduces the background stratification on the shelf by up to 50% and increases both the horizontal and vertical dispersion of threedimensional drifters by up to a factor of 2-3. The studies described above were able to demonstrate the importance of remotely generated internal tides and waves without using global internal wave models. Two recent studies make the same point with the use of global internal wave models. Mazloff et al. (2020) showed that the IGW continuum spectrum in a regional model run without boundary conditions that included remotely generated internal waves was deficient relative to observations. In contrast, a global internal wave model (MITgcm48) analyzed over the same region featured a more energetic IGW continuum that was closer to observations. One might describe Mazloff et al.'s result as a "null result"; a regional model run at the same resolution as a global model, and without remotely generated internal waves at the boundaries, is insufficiently energetic at high frequencies. In contrast, a regional model that is so forced at its lateral boundaries, and that also features increased vertical and horizontal resolution over what is feasible in global models, has an IGW continuum spectrum that lies closer to observations and to predictions of the Garrett-Munk spectrum . As seen in Fig. 15, the frequency spectrum of the Nelson et al. (2020) regional simulation with finer vertical grid spacing is little different from the "One-to-one" simulation, in which the vertical and horizontal grid spacings are the same as in the global MITgcm48 simulation used to force the regional model at lateral boundaries. The regional simulation with finer horizontal grid spacing has a shallower frequency spectrum, lying closer to the − 2 prediction of the Garrett-Munk spectrum, but with a flattened slope, deviating from the − 2 prediction, at the highest frequencies. The spectrum of the regional simulation having both finer vertical and finer horizonal grid spacing ("Finer both") lies close to the Fig. 12. Zonally averaged (a) low-frequency (<0.5 cpd), (b) near-inertial, (c) semi-diurnal, and (d) diurnal kinetic energy (KE) in 1 • latitude bins estimated from MITgcm48 (blue) and drogued drifter data (black). The colored shading shows the 95% confidence interval determined using a boostrapping resampling approach. Note that near-inertial KE at latitude range 10 • S to 10 • N is not included because inertial motions there are indistinguishable from low-frequency variability in the spectra. Reproduced from Fig. 3 of Yu et al. (2019), ©American Geophysical Union, Wiley Online Library, used with permission.
− 2 prediction over all resolved frequencies, without the deviation at the highest frequencies seen in the "Finer Delta x" simulation. The "Finer-Both" simulation also lies closer to the theoretically predicted "consistency relations" than does the "Finer Delta x" simulation (see Fig. 7 of Nelson et al., 2020), indicating that both vertical and horizontal grid spacing are important in the modeling of IGWs.
We expect that regional models forced at their lateral boundaries by global internal wave models will have many uses in the coming SWOT era, which will emphasize internal tides and gravity waves, their interactions with mesoscale and submesoscale eddies, and the signatures of all of these flows in SSH measured by SWOT (and in velocity measured by missions such as the airborne S-MODE mission and the proposed space-based Odysea mission). Such regional models can also be used to better quantify the mechanisms involved in the formation of the IGW continuum spectrum. Pan et al. (2020) examined the regional simulation outlined by Nelson et al. (2020) through the lens of wave turbulence theory, and determined that the "induced diffusion" mechanism (McComas & Bretherton, 1977) is primarily responsible for the IGW spectral formation in the regional simulation. The parametric subharmonic instability (PSI) and elastic scattering mechanisms are less important, at least in the particular regional simulation studied. The relative unimportance of PSI is consistent with analysis of global HYCOM simulations in Ansong et al. (2018). However, it must be kept in mind that the global (and even regional) internal wave models discussed here do not resolve all of the nonlinear wave-wave interactions present in the actual ocean. The relative importance of PSI and other mechanisms may increase as model resolutions continue to increase. Ansong et al. (2018) found a ten-fold increase in PSI energy transfer in HYCOM25 relative to HYCOM12, implying that numerical convergence has not yet been achieved in these simulations.

Summary, ongoing work, and future challenges
General circulation models forced by the astronomical tidal potential as well as atmospheric fields-often referred to as "global internal wave models" here-represent a new frontier in ocean modeling. In this paper we have focused on three early models of this class, with most emphasis on HYCOM simulations run by the US Naval Research Laboratory, secondary emphasis on MITgcm simulations run by the NASA Jet Propulsion Laboratory, and some emphasis on STORMTIDE simulations run at the Max Planck Institute. Another early simulation (Waterhouse et al., 2014) of this type was done with the NOAA GOLD model. New basinand global-scale simulations of this class are being run by the Mercator and Grenoble NEMO groups, the Florida State University HYCOM group (e.g., Chassignet & Xu, 2017-newer simulations include tides), the DOE MPAS group, the NOAA MOM6 group, the Fu et al. (2021) collaboration between Chinese institutions, the UCLA ROMS group, groups in Australia, the UK, and Italy, and other groups. This exciting class of models is here to stay.
When run with state-of-the-art fine horizontal and vertical grid spacings (e.g., 1/12th degree or finer horizontal spacing and of order 40 or more vertical levels), such models simultaneously simulate the eddying general circulation, barotropic tides, internal (baroclinic) tides, near-inertial waves (NIWs), and a partially resolved IGW continuum spectrum (Garrett-Munk spectrum). The IGW continuum spectrum continues to increase in energy as grid spacings are decreased.
Global and basin-scale internal wave models are being used to examine basin-to global-scale impacts of internal waves on ocean acoustics, and to examine internal tide non-stationarity and internal wave interactions with mesoscale and submesoscale eddies. Embedding tides into general circulation models also allows for the study of barotropic tidal interactions with other components of the climate system--for example, tidal effects in rivers and estuaries, tidal mixing fronts on continental shelves, tidal interactions with sea ice, continental ice streams, and floating ice shelves, seasonal changes in barotropic (and internal) tides, and climate-change-induced secular (centennial-scale) changes in barotropic tides. Ongoing work includes an in-preparation intercomparison of global hydrodynamical internal tide models, led by Joseph Ansong, further comparisons with drifters, which represent a promising new avenue for in-situ observational comparisons with internal wave models, and intense continuing work on planning for SWOT, S-MODE, and potential space-based missions for measurements of ocean surface velocities. The US Navy is also planning extensive field work to test global internal wave models.
Challenges for such models include the overly large internal tides in MITgcm48, as documented by Yu et al. (2019), an in-review paper by Arbic et al., and an in-preparation paper by Joseph Ansong. On a related note, the question remains whether the parameterized topographic wave drag employed in HYCOM is the best solution to the problem of damping employed in global internal wave models. The use of wave drag does yield more accurate internal tides (Ansong et al., 2015), and models without wave drag feature tidal motions that are too large (Yu et al., 2019), but internal wave breaking and turbulence at the seafloor is not the only pathway for dissipation of internal tides and waves. Global maps of dissipation inferred from finescale parameterizations applied to hydrographic casts (Kunze, 2017a) indicate substantial dissipation in the upper ocean, which the bottom-focused wave drag parameterization employed in HYCOM does not capture. Wave-wave interactions in the upper ocean represent another pathway to dissipation, one that is only beginning to be captured in global internal wave models as they progress further towards resolving the Garrett-Munk IGW continuum spectrum. Examination of IGW spectrum formation in such models is just beginning to be undertaken. The spectrum is more fully developed in MITgcm48 relative to HYCOM. Another problem in HYCOM, not touched upon earlier in this paper, is the presence of a numerical instability (Buijsman et al., 2016, likely a thermobaric instability (Sun et al., 1999;Hallberg, 2005), in the North Pacific. Thermobaric instabilities, which arise if the compressibility of seawater is not treated with sufficient accuracy, have been fixed in MOM6 (Adcroft et al., 2008), but attempts to suppress instabilities in HYCOM have not yet been successful. Finally, another ongoing challenge for all ocean models,   Nelson et al. (2020). Spectra are computed from the MMP and four regional simulations in a 6 • by 8 • patch of the North Pacific that includes the MMP (see Fig. 1 of Nelson et al. 2020). The regional simulations include a "One-to-one" simulation with the same horizontal and vertical grid spacing as the global MITgcm48 simulation used to force the regional models at their lateral boundaries, a simulation with finer grid spacing in the vertical (z) direction, a simulation with finer grid spacing in the horizontal (x) direction, and a "Finer-Both" simulation with finer grid spacing in both vertical and horizontal directions. There are marks on the x-axis for the inertial frequency f, the semidiurnal lunar tide M 2 , their sum f + M 2 , and the buoyancy frequency from the observations N (computed as the square root of the time mean of N 2 ). The dashed curve above the spectra represents the spectral slope theorized by the GM76 model (Cairns and Williams 1976). Reproduced from Fig. 5 of Nelson et al. (2020), ©American Geophysical Union, Wiley Online Library, used with permission.
including global tidally forced general circulation models, is the need for improved bathymetry, especially the finer scale bathymetric features in coastal areas and in areas of rough seafloor in the open-ocean.
The use of global internal wave models to study open-ocean mixing is just beginning. An initial result (Fig. 16) indicates similarities between inferred depth-averaged dissipation in HYCOM (upper panels) and from ARGO floats (bottom panel). See, for instance, enhanced dissipation along mid-ocean ridges in the Atlantic and Indian Ocean, and in the western Pacific (both northern and southern hemispheres). Further study will focus on the vertical distribution of the dissipation in internal wave models and in observations. The four-dimensional space-time geography of ocean mixing is one of many important problems that global internal wave models will be used for in the future.
Finally, another significant opportunity for extending global internal wave models awaits in the new class of coupled ocean-atmosphere simulations with high resolutions in both fluids (e.g., Bryan et al., 2010;McClean et al., 2011;Griffies et al., 2015;Strobach et al., 2020Strobach et al., , 2022Barton et al., 2021). Coupled models necessarily have frequently updated atmospheric fields, which as discussed above are needed for accurate simulation of near-inertial waves. Analysis of the internal wave fields in some coupled ocean/atmosphere models is just beginning.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements
I thank the reviewers, Nadia Pinardi and Phil Woodworth, for their comments, which greatly improved this manuscript; the editor, Enrique Curchitser, for handling the manuscript; Walter H.F. Smith, for comments on Section 3.2; Chris Hughes, for comments on the discussion of self-attraction and loading (SAL) in Section 3.3, and Shaun Johnston, for providing Figure 3 and pointing me towards the Merrifield & Holloway (2002) reference on the TPXO boundary forcing of that result. Jialin Lin is thanked for handling the editorship of an earlier version of this Fig. 16. Log 10 of depth-mean dissipation rates in W/kg (a) of resolved semi-diurnal low modes in global HYCOM, (b) of resolved semi-diurnal low modes plus parameterized (un-resolved) semi-diurnal high modes in HYCOM (the latter represented by a wave drag), and (c) based on Argo floats (Whalen et al. 2015). All grid cells with Argo values have incomplete coverage in the vertical. The dissipation rates shown are the layer-averaged values. They can be considered upper-bound estimates. Only grid cells that have both HYCOM and Argo data are shown. Reproduced from Fig. 8 of Buijsman et al. (2016), ©American Meteorological Society, used with permission. manuscript, originally intended to be part of an American Geophysical Union (AGU) monograph that was eventually canceled.
I express profound gratitude to the co-authors of all of the papers that we have written on surface tides, internal tides, and internal gravity waves over the years, and to the physical oceanography community for much encouragement as we have pursued this work. For this review paper I would like to single out Jorge Sarmiento for special appreciation. Jorge supported me during a critical time in my career, as this global internal tide research was in its infancy, from National Science Foundation (NSF) grants OCE-0327189 and OCE-0097316 and a Ford Motor Company Carbon Mitigation Initiative Grant.
Finally, I acknowledge with appreciation the generous institutional support that I have obtained over the years; a Jackson School of Geosciences Development Grant, during 2006-2008, faculty start-up funds at the University of Michigan, and funding from the University of Michigan Associate Professor Support Fund (Margaret and Herman Sokol Faculty Awards).

Dedication
I dedicate this paper to my beloved brother Joel Bernard Arbic, who passed away on December 2, 2017. Joel was passionate about a great many things, especially travel, adventure, family, and friends, and was tremendously gifted in endeavors mechanical, athletic, and artistic. Joel will forever be missed by his wife, two sons, parents, two brothers, inlaws, two nieces, other relatives, and numerous friends around the world.