High-Resolution Numerical Modeling of Barotropic Global Ocean Tides for Satellite Gravimetry

The recently upgraded barotropic tidal model TiME is employed to study the influence of fundamental tidal processes, the chosen model resolution, and the bathymetric map on the achievable model accuracy, exemplary for the M 2 tide. Additionally, the newly introduced pole-rotation scheme allows to estimate the model’s inherent precision (open ocean rms: 0.90 cm) and enables studies of the Arctic domain without numerical deviations originating from pole cap handling. We find that the smallest open ocean rms with respect to the FES14-atlas (3.39 cm) is obtained when tidal dissipation is carried out to similar parts by quadratic bottom friction, wave drag, and parametrized eddy-viscosity. This setting proves versatile to obtaining high accuracy values for a diverse ensemble of


Plain Language Summary
We introduce the upgraded computer model TiME, that simulates ocean tides originating from the gravitational attraction of the sun and moon. The model relies on the physics of relevant processes without incorporating actual observations of water level variations. Formerly unconsidered effects that strongly impact tidal dynamics are now included. We discuss the individual impact of these effects on the model accuracy, which is estimated relatively to local measurements from tide gauges. We further compare our results to external tidal models, that employ satellite observations for increased accuracy. Here we find that the upgraded model performs well in the open ocean, and has a reduced accuracy in shallow and coastal waters. The final model setting can simulate tides that recur once or twice per day at a similar level of accuracy. This weak dependence on the excitation amplitude renders TiME especially suited for studying minor tides. Due to their low amplitudes, these tides make up a smaller part of tidal dynamics and are hard to determine with satellite data, thus rendering solutions by our model being free of data constraints valuable. Comparing our solutions with routinely used, empirically motivated estimates of minor tides we show that an increased accuracy is obtained.
SULZBACH ET AL. and other altimetry missions into hydrodynamic models (Carrere et al., 2015;Egbert & Erofeeva, 2002;Taguchi et al., 2014) or using these data to construct empirical corrections to an adopted model (Cheng & Andersen, 2011;Fok, 2012;Ray, 1999;Savcenko et al., 2012), those models are extensively used for the processing of unrelated observations, as e.g., satellite gravimetry missions. Presently, all 34 tidal constituents given by the FES14 tidal atlas (Carrere et al., 2015;Lyard et al., 2006) are directly removed from Gravity Recovery And Climate Experiment (GRACE) and GRACE Follow-on (GRACE-FO) data, and more than 300 additional minor constituents inferred by admittance methods are also subtracted (Kvas et al., 2019). The existing weaknesses in present-day admittance methods, however, have been discussed extensively in the past (Ray, 2017), so that explicit tidal simulations with unconstrained numerical ocean tide models provide potentially valuable information on tidal lines less well constrained by satellite altimetry.
The sensitivity of satellite gravimetry to periodic mass re-distributions in the Earth system is expected to increase even further when the full potential of the satellite-to-satellite tracking by means of laser ranging interferometry (Ghobadi-Far et al., 2020) is also used for gravity field processing. Employing endto-end satellite simulations, Flechtner et al. (2016) found that ocean tide errors are among the top three factors that limit the accuracy of global mass distribution estimates from GRACE-FO. Various concepts of multi-satellite constellations are currently being evaluated by space agencies in Europe, the US, and China for possible implementation as a next-generation gravity mission (e.g., Hauk and Wiese [2020]). Scientific requirements and user demands for such new missions almost always request higher spatial resolution and greater accuracy (see Pail et al., 2015). Equivalently, in order to re-process the already existing data record from GRACE and GRACE-FO into more precise time-series of terrestrial water storage and ocean bottom pressure suited for climate monitoring (Tapley et al., 2019), better ocean tide models are critically important.
While data-constrained tidal models provide highly accurate estimates of tidal constituents in regions where altimetry data is dense (open ocean residuals below 1 cm), model accuracy decreases as the data quality decreases (minor tides, polar, and shelf areas). In effect, the ratio of model uncertainty to signal typically increases considerably for tidal excitations with smaller amplitudes . Even more, additional errors can be introduced, when estimating minor tidal excitations with admittance methods. These deviations might be reduced by the explicit numerical modeling of minor tides.
In this contribution, we present efforts toward extending a hydrodynamic model of ocean tidal dynamics particularly suited to study minor and compound tides. Our work is based on the Tidal Model forced by Ephemerides (TiME; Weis et al., 2008) introduced in Section 2. We describe various improvements to the numerics of the model including the rotation of the poles (Section 3), an extension of the physical model by implementing explicitly the effects of self-attraction and loading (Section 4) and the incorporation of topographic wave drag as a new dissipation mechanism (Section 5). Exemplary for the principal semidiurnal lunar tide M 2 , we will report about the accuracy of the simulated tidal heights both with respect to tide gauge data and the state-of-the-art global tide solution FES14 that is constrained by observations. Various sensitivity experiments are presented documenting the individual contributions of the various changes made to TiME in terms of achieved accuracy (Section 6). The paper is augmented with an assessment of energy dissipation patterns of the model and additional simulations of partial tides in the diurnal and semidiurnal tidal bands (Section 7). Building on the results of previous chapters, we focus in Section 8 on selected minor tides that can be simulated with higher accuracy than solutions constructed from linear admittance estimates on data-constrained models. Finally, the article is closed with a summary (Section 9).

Tidal Model Forced by Ephemerides (TiME)
The barotropic ocean model presented in this contribution is based on decade-long work to simulate global ocean tides in Hamburg, Germany. Starting from the fundamental work of Wilfried Zahel (Zahel, 1977, 1978, unconstrained hydrodynamic models were used to quantify the contributions of ocean tides to Earth rotation (Seiler, 1991), the evolution of tides since the last glacial maximum, and its consequences for oceanic torques acting on the solid Earth (Thomas & Sündermann, 1999), interactions among SULZBACH ET AL. 10.1029/2020JC017097 2 of 21 ocean tides and the general circulation (Thomas et al., 2001), and the identification of free barotropic normal modes in the world's ocean under the influence of friction and sea-bottom deformations caused by surface loading (Zahel & Müller, 2005).
We start our work with the Tidal Model forced by Ephemerides (TiME) as described by Weis et al. (2008), which simulates global barotropic tidal dynamics by solving the nonlinear shallow water equations (e.g., Pekeris, 1974) parametrizes horizontal turbulence as eddy viscosity. Here the horizontal eddy-viscosity coefficient A h and the earth-radius R earth = 6,371 km were introduced. The principal forcing term g∇ ζ eq = ∇V tid can either be an individual partial tide excitation or the full luni-solar tidal potential as quantified by the ephemeridies of sun and moon (Bartels, 1957). Since TiME considers nonlinear accelerations arising mainly from bottom-friction (−r/H|v|v), advection (−(v·∇)v) and the wave-drift-term (−∇·(ζv)), full ephemeridic forcing enables interactions between individual partial tides (compare also Einšpigel and Martinec [2017]). Time-variable surface loading from ocean tidal elevations causes elastic deformation of the sea-bottom that induces a secondary gravity potential (Hendershott, 1972). Weis et al. (2008) approximated the effect as a fraction ϵ of the local tidal elevation, setting ζ SAL = ϵζ as introduced before by Accad and Pekeris (1978). Typical values for ϵ are between 0.08 and 0.12.
When considering forcing by a single tidal constituent of frequency ω, the tidal solution converges to  describes the part of the field oscillating at frequency nω, comprised of in-phase , consequence of tidal rectification (e.g., Pérenne & Pichon, 1999), as well as the overtides with n ≥ 2 are the result of said nonlinear interaction, representing minor contributions to the tidal flow field. Equations 1 and 2 are solved on a global, regular latitude-longitude grid at a resolution of  1 12 employing a semi-implicit finite-difference algorithm as described in Backhaus (1983) and Backhaus (1985). Since the smallest zonal grid-cell size limits the time-step, the zonal resolution is reduced at two certain latitude circles toward the North Pole and finalized by a spherical cap to avoid a polar coordinate singularity. Numerical experiments were based on global Gebco (GEBCO Compilation Group, 2019) and Etopo1 (NOAA, 1988) bathymetries, where the best results presented in Weis et al. (2008) were obtained with Gebco data. Results from this model configuration will be presented in Section 6 in comparison to the latest results that incorporate various extensions to the model physics as outlined in the following.
SULZBACH ET AL.

Pole Rotation Scheme and Tide-Raising Potential
Positioning a pole of the numerical grid in the middle of the Arctic ocean requires the application of a spherical cap in order to work around the singularity. This can be avoided by shifting the pole locations to dry grid cells that are ignored by the solver. In order to adapt the numerical scheme to the rotated pole location, we transform the explicitly space-dependent parts of the shallow water Equations 1 and 2, comprising bathymetric function H(ϕ, λ), tide-raising potential V tid = gζ eq (ϕ, λ) and Coriolis vector f(ϕ, λ) into a new coordinate system which is spanned by the rotated coordinates ϕ′(ϕ, λ), λ′ (ϕ, λ). This coordinate transformation where ϕ p , λ p is the location of the new North Pole on the unrotated grid, can immediately be used to transform the Coriolis factor f(ϕ) = f sin(ϕ). Due to the different local definitions of north and east on a rotated grid, special attention has to be spent when transforming vectorial quantities (e.g., the tidal velocity field v).
In total, 4 different pole locations (see Figure 1) are implemented for experiments: The original latitude-longitude grid with the pole in the Arctic at 90° N (following denoted Arc), and two realizations with land-covered poles at Argentina and China (Chi) and Antarctica and Greenland (Gre), with dry areas around the poles having a radius of at least 4°. An additional configuration with the South Pole in Australia and the North Pole in the North Atlantic (Aus) was also implemented to highlight the potential errors that might be introduced from a poorly selected grid configuration. In the following, the consequences of the pole rotation scheme for bathymetry and tide-raising potential are discussed.
Previously, TiME was run with Gebco and Etopo1 bathymetries. More recently, the high-resolution Rtopo2 data (Schaffer et al., 2016) became also available. Based on Gebco data, the Rtopo2 grid uses additional data sources in polar latitudes with a special focus on accurately representing sub-ice-shelf cavities as the position of the ice-shelf grounding line is known to have a strong effect on tidal oscillations (Arbic, Karsten, & Garrett, 2009;Wilmes & Green, 2014). We included these areas by considering the water-column height as the difference between the ice base and bedrock depth. Since the resolution of our tidal model (5 arc minutes) is well below the resolution of the available bathymetric data (30 arc seconds), special attention has been paid to how the necessary resolution reduction is optimally performed. We employ an interpolation strategy motivated by perturbation theory, that uses conservative interpolation of inverse depth (1/H), to directly derive bathymetric maps for the respective pole-orientations (Arc/Chi/Gre/Aus) from high-resolution data (see supporting information). The minimum water depth was set to 15 m as by Weis et al. (2008), while values of 10 or even 5 m have a negligibly small influence on tidal dynamics. The Caspian Sea Level (CSL) that is subject to climatic rapid changes compared to the open ocean (Chen et al., 2017;Naderi Beni et al., 2013) has been fixed to −26.5 m.
The luni-solar tide-raising potential acting on the oceans at epoch t SULZBACH ET AL.
10.1029/2020JC017097 4 of 21 can be decomposed into orthonormal, real-valued spherical harmonics Y lm with time-dependent coefficients. These coefficients can further be decomposed into temporally harmonic constituents at discrete frequencies ω i with amplitudes , lm lm i i s c (Hartmann & Wenzel, 1995). Here lm P are orthonormal Legendre-Functions (e.g., Heiskanen & Moritz, 1967) and is a combination of Love numbers (Love, 1909) that encrypts the action of (elastic) body tides of the solid Earth on oceanic tides. The frequency dependence of  originates from a coupling of tidal motions to the near diurnal free wobble (NDFW), one of the fundamental eigenmodes of the solid Earth (Wahr, 1981;Wahr & Sasao, 1980), which has prominent effects on admittances in the K 1 tidal group (Ray, 2017).
We use the tidal catalog published by Hartmann and Wenzel (1995) considering contributions up to degree l = 3, so that minor ocean tides even of third degree origin (e.g., M 3 ) can be simulated.  ( ) b l i a up to degree-3 has been taken from Spiridonov (2018) (model 9), where diurnal Love numbers are calculated by linear interpolation from neighboring data-points (linear Admittance, see Table 1). Please note that due to the frequency-dispersion and degree dependence introduced by  ( ) b l i a , the ocean tide-raising potential differs from the initial gravitational potential created solely by sun and moon.
Since V tid is composed of few low-degree spherical harmonics, transformation to rotated-pole grids can be performed efficiently by utilizing Wigner-D-functions initially derived in the framework of quantum theory of angular momentum (e.g., Warschalowitsch et al., 1988). Dependent on ϕ p and λ p a single component transforms into a combination of spherical harmonics of the same degree. The rotation is therefore entirely described by a transformation   )( Here ( ) lm lm i i c s describes a set of 2l + 1 coefficients of the same degree. This transformation is implemented by using the rotation algorithm by Gooding and Wagner (2010) following routines published by Risbo (1996).
The maximum dry radius around the North and the South Pole of 4° defines the smallest zonal resolution on this equally distributed lat-lon grid. Since this distance restricts the minimum possible time-step, it was necessary to introduce additional zones of reduced zonal resolution to the Southern Hemisphere. In Weis et al. (2008), this was solely implemented for the Northern Hemisphere since former bathymetries provided a much bigger dry radius around the South Pole (compare Figure 5). This modification allows the enlargement of the time-step by the factor of 4 and thus provides a considerable speed up. For all experiments shown in this paper, zonal resolution is first halved at ϕ = ±60° and afterward at ϕ = ±75° ensuring that zonal resolution varies between the margins of  1 12 and     1 1 4 sin (4 ) 12 43 on the gre-grid (    1  38 on chi). The SULZBACH ET AL. 10.1029/2020JC017097

of 21
Tide

Self Attraction and Loading
The elastic deformation of the sea-bottom in response to a water load induces a secondary gravity potential that itself excites a response in ocean dynamics. The inclusion of this so-called self attraction and loading (SAL) potential V SAL = gζ SAL has been found to be important many years ago (Hendershott, 1972). Since the additional local forcing depends on the entire global distribution of tidal elevations, it can be computed by a convolution integral. Such computations are rather costly, so that local parametrization of SAL (e.g., Accad & Pekeris, 1978) were often applied. However, present-day accuracy demands on global ocean tide models make a rigorous consideration of SAL indispensable (Ray, 1998). We consider SAL in the spectral domain, where the computation translates to a set of algebraic operations on spherical harmonic coefficients ζ nm , reading Here, α l = 1 + k l − h l is a combination of degree dependent Load-Love numbers (LLNs) (Munk & MacDonald, 1960) that encode the effects of additional gravitation (1), as well as seafloor deformation (h l ) and potential shift (k l ) caused by the varying water mass. The factor α l /(2l + 1) decreases with rising degree l, within ensuring slow but steady convergence, thus avoiding Gibbs Phenomenon at coastal load discontinuities (Agnew, 2007;Hewitt & Hewitt, 1979). We employ LLNs of Wang et al. (2012) (PREM Earth model) with a correction to represent low degree LLNs in the frame of figure (Blewitt, 2003). Further we set l max , the maximum considered degree in Equation 10, to 1024, thereby guaranteeing sub-mm accuracy in SAL-representation apart from some shelf areas and estuaries as also incorporated and discussed before by Schindelegger et al. (2018). The computational burden is thereby shifted to repeated transformations between the spectral and the spatial domain, which are efficiently handled with the highly optimized SHTNS-package by (Schaeffer, 2013).

Topographic Wave Drag
A realistic representation of dissipative forces is critically important for modeling the system's resonant coupling to oceanic normal modes and has been the main issue in numerical tidal modeling for a long time.
Besides turbulent friction of the bottom boundary-layer, which is strongest in shallow water, the inversion of satellite altimetry data indicated that significant dissipation is happening in the deep ocean located at prominent topographic features (Egbert & Ray, 2000, 2001. The reason for this energy loss from barotropic tides is the excitation of baroclinic tidal motion, known as internal tides (e.g., St. Laurent et al., 2012;Wunsch, 1975). This phenomenon is hard to observe with remote sensing techniques since the induced changes in the ocean state are nearly entirely internal. However, it was possible to detect the minuscule ocean surface fingerprint of internal tides with space-geodetic techniques (e.g., Zhao, Alford, Girton, Johnston, & Carter, 2011;Zhao, Alford, Girton, Rainville, & Simmons, 2016). Further dissipation is expected from shallow water processes, horizontal turbulent friction and friction at ice-water boundaries in polar latitudes.
The explicit simulation of internal tides requires a depth-resolving, baroclinic ocean general circulation model (Arbic et al., 2012) and is hence much more complex than barotropic tidal modeling as pursued here. Nonetheless, parametrizations for baroclinic processes are available from theoretical considerations and have shown to be able to capture the dissipation patterns accurately ( Nycander, 2013). We implement a parametrization proposed by (Nycander, 2005) that builds on previous considerations of Bell (1975) and Llewellyn Smith and Young (2003). It is based on linear wave theory and produces dissipation patterns that closely match available observations. The wave drag information is described by the second rank tensor . The drag then enters the tidal partial differential equation (PDE) as an additional dampening acceleration We use depth-resolved hydrographic data from the World Ocean Database for salinity  and temperature (Locarnini et al., 2019) with the TEOS-10 equation of state (McDougall & Barker, 2011) to compute a global map of the buoyancy-frequency N(ϕ, λ, z), where z is a depth coordinate. The excitation of internal waves is a strongly frequency-dependent process that differs for diurnal (ω d = Ω), semidiurnal tidal species resulting in different wave drag tensors for each species ). As Buijsman et al. (2015), we follow the approach introduced and developed by Nikurashin and Ferrari (2011), Melet et al. (2013), and Scott et al. (2011) to reduce potentially overestimated wave dragstrength at supercritical slopes. To achieve this, the drag strength is normalized at supercritical slopes by its criticality squared to compensate for overestimated dissipation. Further, we allow for a tuning parameter κ to adjust the overall dissipation strength and introduce a cutoff depth of 150 m below which shelf oceans are assumed to be well-mixed. For more detail about calculation and properties of  consider the supporting information.
With the discussed wave drag-parametrization, TiME now includes three dominant dissipation mechanisms comprised in the dissipation operator . This further encompasses parameterized turbulent horizontal eddy-viscosity R (Zahel, 1977) and quadratic bottom friction. Alltogether, Equation 1 can be rewritten as By tuning the parameters κ, A h and r model dissipation channels can be weighted individually.

Tidal Elevations for M 2 from TiME
In order to highlight the importance of individual model changes to TiME as given in the previous sections, we now report results from a number of sensitivity experiments for the principal lunar tide M 2 as outlined in Table 2. The model performance will be benchmarked against a data set comprised of 151 tide gauge stations compiled by Ray (2013) as well as the global state-of-the-art tidal atlas FES2014 (Carrere et al., 2015;Lyard et al., 2006), that was produced by Noveltis, Legos and CLS and distributed by Aviso +, with support from CNES(https://www.aviso.altimetry.fr/). Misfits will be reported in terms of time averaged rms that can be further averaged over a certain ocean domain D o with area A o yielding the space-averaged . We calculate averages for shallow water if the ocean depth is smaller than H = 1,000 m (10.4% ocean surface) or open ocean if the depth exceeds this limit (83.0% ocean surface). Both areas are restricted to latitudes with |ϕ| < 66° as altimetry data in these region is dense and guarantees a high quality of derived tidal atlases. When mentioned in the following sections, the results for As an additional benchmark for our model we monitor planetary dissipation conducted by the M 2 tide. Being laid out by Platzman (1984), the theory of planetary dissipation was employed by Egbert and Ray (2000) and Egbert and Ray (2001) to derive estimates of M 2 -tidal dissipation utilizing altimetry data-constrained tidal models. Herein, the planetary dissipation field d was derived by evaluating the relation that uses the mean tidal energy consumption field w = ρ sw 〈a c ·v H〉 (work done by tidal forces) and the energy flux field ∇p = ρ sw g∇〈ζHv〉.  Primarily the introduction of topographic wave drag and an improved bathymetric map allowed to decrease global M 2 amplitudes and dissipation rates. Tuning experiments concentrated on finding an optimal ratio between damping by eddy viscosity and topographic wave drag, while the bottom friction was left constant at r = 0.003. Even though the original Nycander scheme does not contain a free, tunable parameter, experiments with κ = 100% showed that additional dissipation is necessary to obtain optimum results as was also found by Buijsman et al. (2015). This can either be provided by increasing A h or κ (altering r worsened the accuracy). Several tuning experiments with κ = 100…225% and   The resulting mean tidal power consumption field w as well as the energy dissipation field by wavedrag-acceleration d wd (see Figure 2) match the results derived with altimetric data (compare Egbert and Ray [2001]).
While the achieved accuracy of experiment RE was the highest in our ensemble, this model setup lacks a solid physical foundation due to the excessive dissipation mediated by parameterized eddy-viscosity. A physically more reasonable setting can be obtained when minimizing the dissipation by eddy-viscosity as pursued by most modern barotropic models (e.g., Egbert et al., 2004;Schindelegger et al., 2018). The obtained open ocean rms-values of 4.03/6.54 cm for W0 increase with respect to setting RE, while the dissipation-overshoot is reduced to 120 GW (5%). On the other hand, the shallow-water accuracy is not altered considerably to 17.86 cm. This can be seen as a trade-off between maximized model accuracy and well-founded model physics that will be beneficial with respect to sensitivity studies (e.g., paleo simulations, climatic impacts). On the other hand, this trade-off is undesirable for high-precision applications such as satellite gravimetry. As the accuracy that has to be sacrificed with setting W0 increases for minor tides (see Section 7), we decide to favor setting RE for the present study.
The amphidromic system and global rms-data for experiment RE is shown in Figure 3. In comparison to FES2014-data, M 2 oscillation systems are predicted precisely with exception of some features around Antarctica (compare Figure 5, top). It is worth noting that the reproduction of large-scale features (e.g., tidal phases defining amphidromic systems) was also possible on a similar level of detail by experiment WE. Exceptions were mainly constituted of bathymetry-induced aberrations around Antarctica (compare Section 6.3). The principal accuracy gain is attributed to a more realistic representation of dissipation, bathymetry, and SULZBACH ET AL.
10.1029/2020JC017097 9 of 21 SAL effects. The remaining critical regions are shelf and coastal areas and especially concentrated around Antarctica and in the North Atlantic shelf areas. This suggests possible origins for these discrepancies in tide-ice interaction as well as in possible bathymetric inaccuracies or insufficient representation of (nonlinear) shallow water effects. With respect to similar modern barotropic tidal models, TiME produces solutions on the same level of accuracy, while the shallow water accuracy is moderately decreased (e.g., +4.4 cm to Schindelegger et al. (2018).

Impact of Pole Location
As a first finding we note similarly accurate results when performing experiments on alternative grids with land-covered poles (compare experiment P2 in Table 2). The accuracy obtained on the Gre-grid could be increased further to 3.63/5.00 cm open ocean rms, when additional tuning was applied (experiment P2b). Since the zonal model-resolution increases toward the numerical poles, the bathymetric information contained in two differently oriented grids differs slightly. This also impacts the wave drag tensor that depends on H(ϕ, λ  . We conclude that optimized parameters vary for experiments on differently oriented grids, but similar accuracy is achievable. Overall, the rotated pole scheme did not improve the global accuracy level of present-day tides significantly, since tidal elevations in the Arctic (near to the former pole cap) are diminutive (experiment P1). Nevertheless, deviations induced by a pole cap situated in an area of presumably high tidal elevation are a source of imprecision (P3) and should be avoided. The non-optimally placed pole cap on the "Aus"-grid resulted in altered and diminished dissipation as well as in an increased rms  4 75 5 57 . / . cm.
An additional benefit of the implemented pole-rotation scheme lays in its versatility: The mitigation of the large pole-cap grid cell allows the unbiased study of historical situations in which tidal elevations in the Arctic might have been significant, as proposed by Griffiths and Peltier (2008) and more recently by Velay-Vitow et al. (2020) for the Last Glacial Maximum. Additionally, alternative grid orientations can be used to guarantee approximately equal aspect ratios for grid-cells in the Arctic (e.g., when using the Chigrid), further recommending TiME to be used for studies of Arctic tides.

Impact of Bathymetry
Additional experiments were performed with bathymetries constructed from Etopo1 (Amante & Eakins, 2009) and Gebco data (GEBCO Compilation Group, 2019) treating sub-ice-shelf cavities as dry gridcells. This configuration resembles the bathymetric maps used in Weis et al. (2008). Interpolation to the model's resolution was done using first-order conservative interpolation.
The results (experiment B1, B2) show that ignoring the effects of Antarctic sub-shelf cavities on ocean tide resonances leads to large scale deviations of the displayed amphidromic systems, especially in the southern ocean (Wilmes & Green, 2014) (compare Figure 5). The most striking deviation hereby occurs in the Weddell-Sea. As the Rtopo2-bathymetry is mainly based on Gebco-data, the model setup for experiment B1 can be seen as a blocking experiment for the Antarctic shelf regions. Blocking experiments are useful to investigate the back-action of shelf-tides on open ocean tides. Arbic, Karsten, and Garrett (2009) conducted such simulations for a number of shelf areas (e.g., Patagonian Shelf, Hudson Bay) and also considered analytical solutions. Both approaches predict that blocking a near-resonant shelf-region enhances the amplitude of the open ocean tide, as it is shown in Figures 5a and 5c for the Southern Ocean. As the shelf-blocking increases the open ocean amplitudes, tidal dissipation is also increased and contributes to the overestimated dissipation in experiment WE (B1: +150 GW, B2: +350 GW). It is possible to reduce this overestimation slightly by enhancing the dampening forces. However, this only leads to minor improvements and cannot rectify the imprecisely represented oscillation systems as depicted in Figure 5, it is therefore not further investigated.
SULZBACH ET AL.
10.1029/2020JC017097 11 of 21 In summary, the conducted experiments highlight the irreplaceability of constructing a realistic bathymetric map. Further, the results point out that poorly represented areas can have strong near-and far-field effects on tidal dynamics even if they have only a small spatial extent.

Impact of SAL
In this section, we discuss the impact of a number of SAL representations on the simulation results. When evaluating Equation 10 up to l max = 100 we find that the open ocean rms does not increase indicating sufficient handling of this effect with regard to open ocean tidal dynamics (experiment S1). The shallow water rms is also not altered considerably, which might be due to a generally less precise model performance in shallow waters. When further decreasing the maximum degree to l max = 10 the open ocean rms increases to 3.99/5.63 cm (S2) which is still a significantly more precise result than S3, that was obtained by the local SAL-parametrization by (Accad & Pekeris, 1978) (5.41/6.69 cm). However, this result can be seen as a valuable improvement when comparing it to completely neglecting the effect (experiment S4). This leads to profound misrepresentation of tidal phases and strongly increased misfit. Overall, the estimated planetary dissipation increases with a less precise SAL representation (S3: +240 GW, S4: +500 GW). Thus, the local SAL-parametrization contributes significantly to the overestimated dissipation in experiment WE.
As discussed by Müller (2007) the inclusion of SAL primarily leads to a phase-shift of oscillation systems by altering properties of the underlying normal modes. This phase shift is well approximated in first order by applying an effectively reduced gravity factor (S3). However, to represent the precise far-field effects of the SAL a treatment in terms of spherical harmonics with l max ≥ 100 is necessary. To precisely represent near field effects the maximum degree l has to be extended to higher values (Schindelegger et al., 2018). Since the SULZBACH ET AL.
10.1029/2020JC017097 12 of 21 efficient handling of SAL-transformation by the SHTNS-package (Schaeffer, 2013) does not considerably increase computation time we treat the SAL effect up to degree l max = 1,024 in our experiments.

Impact of Topographic Wave Drag
The only dissipation channel that leads to considerable deep ocean dissipation included in the model is the excitation of internal tides by topographic effects. Since deep ocean dissipation is an experimental matter of fact, it does not surprise that complete neglection of this effect has severe impacts on the achievable accuracy (experiment W3) causing an rms increase of +4.64/+4.91 cm and a surplus dissipation of an additional 260 GW. Finding an optimum interplay between the present dissipation channels on the other hand is more complicated. When abstaining from tuning wave drag-strength, e.g., setting κ = 100%, best results are achieved by allowing significant dissipation by eddy-viscosity (W1). If one instead decides to increase wave drag substantially (W2), as for example, done by Buijsman et al. (2015), the open ocean results slightly worsened in comparison to experiment RE, without improving shallow water rms. As discussed in Section 6.1 a minimization of the obtained wavedrag dissipation leads to the physically well-founded model setting W0 while the open ocean accuracy moderately deteriorates (+0.64/1.71 cm). This setting should always be favorable with respect to sensitivity experiments that benefit from realistically represented tidal physics.
Overall the tuning of the dissipation channels suggested that best results are obtained when wavedrag dissipation contributes about 900 GW to planetary dissipation, which is close to the expected value. The slight tuning of the wave drag tensor (factor 1.25 to the original Nycander-tensor) stresses that it is based on a reliable theoretical basis and can be expected to provide ad-hoc precise results. This is a valuable result when it comes to adapting the model to other tidal groups or paleo settings.

Impact of Spatial Resolution
were designed to represent a similar physical situation as chosen for experiment RE. Therefore, dissipation channels were tuned to achieve comparable dissipation ratios leading to an effective increase in A h . We emphasize, that altering the model resolution renders repeated model tuning necessary. Parameters cannot be transferred directly without altering tidal dynamics.
We observe that, while overall dissipation decreases, open ocean rms-values are considerably increased to 5.21/7.95 cm for R2/R1 (compare Egbert et al. [2004]). The cause for this might originate from the model-inherent resonant behavior of oceanic tides. With reduced spatial resolution, the geometry of the ocean basins determining oceanic normal-modes cannot be properly represented. The resulting slight shifts in resonance frequencies then strongly impact tidal dynamics, especially in shallow waters, where oscillation systems reside on smaller spatial scales. Thus, to further increase the precision and accuracy of TiME, an increased resolution beyond  1 12 should be considered.

Additional Tidal Excitations
In this section, we present simulation results for additional partial tides. We selected partial tides that differ in excitation amplitude, frequency as well as in excitation pattern from M 2 to test the sensitivity of the preferred model setting RE. The overall aim is to demonstrate model setting RE as robust for simulating partial tides of differing character on a comparable level of accuracy.
As discussed in Section 3, excitation patterns relate to the spatial dependence of the partial tide forcing that is proportional to the spherical harmonic functions Y lm . Within, l defines the degree, m the order of the spherical harmonic function, where m = 0, 1, 2, 3 further enumerates the tidal species (0: long period, 1: diurnal, 2: semidiurnal, 3: terdiurnal). Since the tidal forcing strength spans several scales for different partial tides, the level of accuracy obtained for different partial tides cannot be compared directly to each other without considering the overall signal amplitude.
SULZBACH ET AL.

10.1029/2020JC017097
13 of 21 To facilitate this comparison we introduce the admittance function that relates the tidal response, expressed by its elevation   i , to its g-normalized forcing amplitude (compare supporting information). Hereby, Z lm is only evaluated at discrete tidal frequencies for partial tides with forcing pattern Y lm . Since the tidal PDEs are only weakly non-linear and tidal frequencies within one tidal band only differ slightly, Z lm takes a related shape for each partial tide. Hence, it can be used to compare the response strength and especially the relative level of accuracy compared to the excitation strength for individual excitations by considering x (compare Equation 13). Please refer to the supporting information where the tidal potential catalog used for this study can be obtained.

Semidiurnal Tides
The K 2 -excitation is a semidiurnal partial tide of second degree origin (l = 2), thus representing another evaluation point of admittance Further, we considered the tidal response to ν 2 tide (approx. 4% M 2 -forcing) as a third evaluation point of Z 22 . Here, nonlinear effects will play an even more important role, due to the diminished excitation amplitude, while frequency difference to M 2 is reduced  The validation of ν 2 resulted in a rms of 0.19/0.18 cm, revealing a moderately enhanced level of inaccuracy compared to K 2 -results. The reason for this could be found in an imperfect representation of non-linearities in tidal dynamics. Nonetheless, the results demonstrate that TiME is able to perform simulations on a similar scale of accuracy within one tidal band (in this case Z 22 ) without the need to adapt model parameters for each simulation. Results for these partial tides are shown in Figure 6.
On the other hand, monitoring tidal dissipation reveals increasingly altered weights for individual dissipation channels. While K 2 dissipates 34.7 GW, partitioning as (D bf , D wd , D ed = 4.9/13.6/16.3 GW) the distribution for ν 2 (D bf , D wd , D ed = 0.34/1.33/2.21 GW) is even more shifted toward a dominant eddy-dissipation. Not surprisingly, the dissipation by quadratic bottom friction ∼|v| 3 is strongly reduced compared to dissipation by linear forces ∼|v| 2 . Remarkably, the overall dissipation lost by reduced bottom friction is transferred to D ed , while deep ocean dissipation by D wd amounts to a comparable fraction (33.6/38.4/34.2% for M 2 /K 2 /ν 2 ) of dissipation.
On the other hand, simulation results obtained with model setting W0 (compare Table 3) reveal decreasingly accurate results for minor tides. While for M 2 the decrease of open ocean accuracy was at a moderate level of 19% the toll for implementing setting W0 increases by 54 (84)% for the K 2 (ν 2 ) excitation. The reason for this could reside in a possible overestimation of quadratic (non-linear) shelf-dissipation mechanisms when adapting the setting W0 for the M 2 tide. The importance of a nonlinear dissipation channel (quadratic bottom friction) reduces strongly for minor tides, which could result in distorted ratios between deep ocean and shelf-dissipation. As the accentuation of an alternative linear dissipation mechanism (eddy-viscosity) in model setting RE improves the results for minor tides drastically, it could be beneficial to consider novel linear shelf-dissipation mechanisms for the precise prediction of minor tides. In spite of the poor physical justification of dominant eddy-dissipation it might thus be its linear nature that benefits the accuracy of hydrodynamical tidal simulations.
SULZBACH ET AL.

Diurnal Tides
The K 1 tide is the principal excitation in the diurnal band with a magnitude of 58% M 2 -forcing strength. It is important to note that the resulting forcing applied on ocean masses is enhanced by 6.2%, compared to an equivalent forcing at semidiurnal frequencies due to the NDFW-resonance (compare Table 1). In contrast to Z 22 , the tidal excitation pattern is proportional to Y 21 and the strongly dispersive wave drag-parametrization is further limited to low latitudes with |Φ| <30°. This causes the admittance function Z 21 to take a different shape compared to Z 22 . Tidal elevations are now concentrated in the North Pacific, Indian, and Southern Ocean. Validation performance yields an rms of 0.90/1.32 cm (compare Figure 7). A possible explanation for the overall higher accuracy might be the larger scales of diurnal oscillations systems, that are less sensitive to detailed bathymetric information. Moreover, the overall smaller admittance in the diurnal tidal band (Z 21 ) indicates less resonant tidal behavior, and thus, less sensitivity to slight changes. We further display a second diurnal oscillation system (Q 1 ) to demonstrate the achieved accuracy over multiple scales of excitation amplitude. While the forcing-strength is reduced by 86.4% compared to K 1 , validation accuracy is on a similar level, yielding 0.19/0.25 cm open ocean rms. Due to the shift in excitation frequency by  1.70 h , the admittance-function Z 21 changes notably. This also leads to an altered concentration of uncertainty in (shelf)-regions.
Concluding this chapter, we want to stress that the achieved high accuracy for 5 partial tides of diverse character proves the model setting RE as suitable and favorable over other settings for high accuracy applications.

Tidal Solutions for Satellite Gravimetry
In the previous section, it was shown that with model setting RE it is possible to simulate minor tides at a similar level of relative accuracy for a wide range of tidal frequencies, excitation amplitudes, and excitation patterns. On the other hand, the precision of satellite data-constrained partial tide solutions depends on the available data quality. This quality is not constant but depends on the respective frequency and domain of tidal observations. Typically, polar tides are less accurately known since many SULZBACH ET AL.

10.1029/2020JC017097
15 of 21  satellite orbits are limited to |ϕ| < 66°. This leads to prominent GRACE residuals (Ray, Loomis, et al., 2019;Ray, Luthcke, & Boy, 2009;Wiese et al., 2016) in polar seas. The same is true for minor tidal constituents that are routinely considered for gravity field de-aliasing (Savcenko & Bosch, 2010), but are often not provided explicitly by data-constrained tidal atlases as the data quality is poor. As this might change at sometime with continuously extended altimeter time-series (Ray, 2020), minor tides are currently routinely derived by admittance assumptions (Gérard & Luzum, 2010), that are prone to reduced accuracy especially in shallow waters and ice shelf regions (e.g., Pedley et al., 1986) that are governed by nonlinear processes. In this framework, it is natural to ask if purely hydrodynamic solutions can perform more accurate than data-constrained solutions. While this is certainly not true for major constituents (M 2 , K 1 , …) we want to take a closer look at minor tidal constituents, that are yet relevant for gravity field de-aliasing.
Most promising results can be expected by tides at the edges of tidal bands, as these contain the potentially largest errors when utilizing linear admittance theory (Ray, Luthcke, & Boy, 2009). Thus, we choose the diurnal 2Q 1 and OO 1 tides (1.0%/1.8% M 2 -forcing strength) as first test cases. The validation is performed with a set of tide gauges stations of predominantly coastal character (TICON, Piccioni et al., 2019). Additionally, we probe Q 1 -results in the Antarctic region by validation with a data set of Antarctic tide gauges stations (Howard & Padman, 2020). The respective solutions are either directly included in the FES2014-atlas (Q 1 ), or derived via linear admittance supported by Q 1 , O 1 -and K 1 tide, where we consider perturbations in the tidal potential height by the NDFW-resonance (compare Tab. 1). Hereby we assume ∂ ω Z lm = c 1 + c 2 ω (compare Equation 15), evaluate the constants c 1 , c 2 by the two closest supporting points and use the result to extrapolate the results linearly to 2Q1 (support Q 1 , O 1 ) and OO 1 (support O 1 , K 1 ).
As the distribution of the respective rms-values is considerably askew, especially for the TICON-dataset (compare Figure 8), we decided to utilize the distribution median as an effective validation metric. The median will be listed in the following for the (TiME/FES2014)-distributions. TiME performs on a similar level of accuracy as FES2014 for Q 1 in the Antarctic domain (1.04/0.83 cm). While this is already quite remarkable for an unconstrained model, it proposes that the local, nonlinear particularities of tides below ice shelves must be considered explicitly to obtain more accurate results. On the other hand, the accuracy for 2Q 1 (0.07/0.11 cm) and even more for OO 1 (0.17/0.34 cm) is clearly increased when TiME solutions are employed. From this improvement of validation accuracy with respect to linear admittance solutions, we draw the conclusion that the utilization of TiME-solutions for certain partial tides will result in a reduction SULZBACH ET AL.
10.1029/2020JC017097 16 of 21 Figure 7. Dimensionless admittance-function Z 21 for K 1 -(a) and Q 1 tide (b) and corresponding rms-values for Z 21 (c + d). Rescaling Z 21 to real ocean elevations amounts to 88.14/12.00 cm for K 1 /Q 1 . Note the overal lower response level compared to Z 22 (Figures 3 and 6).
of the aliased tidal signal in GRACE-data. This de-aliasing potential for GRACE-data exhibited by our minor tide solutions emerges from the integrated improvement of the model and its high accuracy over a wide range of partial tide.

Summary and Outlook
In this paper, we introduced several modifications to the barotropic tidal model TiME (Weis et al., 2008), which resulted in a considerable increase of the open ocean accuracy.
First, we showed that the introduction of a comprehensive pole rotation scheme allows us to get rid of numerical artifacts potentially induced by the former pole cap handling. The realization of two "pari passu" grid orientations with land-covered numerical poles further allowed to estimate resolution-connected model uncertainties.
Second, the introduction of a non-local online treatment of the effect of SAL (Ray, 1998) and the implementation of a wave drag-parametrization (Nycander, 2005) allowed for a substantial increase in the model accuracy. We further discussed the relevance of constructing optimized bathymetric maps from different available global data sets. The resulting update on TiME's bathymetry evoked another significant increase in model accuracy, especially due to the inclusion of sub-ice-shelf cavities (Wilmes & Green, 2014). The individual contributions of these updates to model accuracy are summarized in Table 4.
Due to the removal of the numerical pole cap, TiME proved to be highly versatile in simulating arbitrary oceanic regions with the same level of accuracy. An open question is the correct representation of (nonlinear) processes in shallow water, beneath ice shelves, and coastal areas, where the model accuracy considerably drops. Though dissipation by eddy-viscosity (a linear dissipative force) increased the overall model accuracy considerably, the obtained high values for A h are hard to justify. The question of shallow-water dissipation should be readdressed, also given the persistently overestimated M 2 -tidal dissipation.
Tuning experiments of the updated model resulted in a set of model parameters, that equally distributes M 2 tidal dissipation to friction by wave drag, quadratic bottom-turbulence, and parametrized eddy-viscosity. The set of model parameters proved robust toward the simulation of diurnal tides and minor tidal excitations, where results with comparable relative accuracy were obtained. On the other hand, parameters had to be adjusted slightly, when changing the model's resolution. While the discussed setting is favorable for the gravimetric applications we are envisioning, a second, physically better-founded setting was derived, that is favorable for sensitivity studies, or paleo experiments.
The achieved model performance qualifies TiME as a purely hydrodynamic tidal model for simulation of modern-day tides. While absolute model deviations from tide gauge data are considerably bigger than results obtained by data-constrained tidal models for major tides, we could show that the accuracy for minor tides can be improved. The same might be possible for polar tides of major origin if crucial polar tidal processes as sea-ice-friction are considered. This potential arises from TiME's SULZBACH ET AL.
10.1029/2020JC017097 17 of 21   independence off satellite data and allows for an almost constant relative model accuracy over multiple scales in tidal forcing strength. To fully benefit from the de-aliasing potential of the obtained solution a comprehensive study focusing on the accuracy improvement of all relevant minor tides by unconstrained simulations should be conducted and augmented with direct estimates of GRACE-gravity field solutions.

Data Availability Statement
The tidal data presented in this study can be obtained from https://doi.org/10.5880/GFZ.1.3.2021.001 in the form of Stokes coefficients including load-tide mass variations (Sulzbach et al., 2021).