Barotropic tides in MPAS-Ocean (E3SM V2): impact of ice shelf cavities

. Oceanic tides are seldom represented in Earth sys-tem models (ESMs) owing to the need for high horizontal resolution to accurately represent the associated barotropic waves close to coasts. This paper presents results of tides implemented in the Model for Prediction Across Scales–Ocean or MPAS-Ocean, which is the ocean component within the U.S. Department of Energy developed Energy Exascale Earth System Model (E3SM). MPAS-Ocean circumvents the limitation of low resolution using unstructured global meshing. We are at this stage simulating the largest semidiur-nal (M 2 , S 2 , N 2 ) and diurnal (K 1 , O 1 ) tidal constituents in a single-layer version of MPAS-O. First, we show that the tidal constituents calculated using MPAS-Ocean closely agree with the results of the global tidal prediction model TPXO8 when suitably tuned topographic wave drag and bottom drag coefﬁcients are employed. Thereafter, we present the sensitivity of global tidal evolution due to the presence of Antarctic ice shelf cavities. The effect of ice shelves on the amplitude and phase of tidal constituents are presented. Lower values of complex errors (with respect to TPXO8 re-sults) for the M 2 tidal constituents are observed when the ice shelf is added in the simulations, with particularly strong improvement in the Southern Ocean. Our work points towards future research with varying Antarctic ice shelf geometries and sea ice coupling that might lead to better comparison and prediction of tides and thus better prediction of sea-level rise and also the future climate variability.


Introduction
Ocean tides are generated due to the gravitational force from celestial bodies -the Sun, Moon and the Earth. Such gravitational dynamics causes the sea surface height of the global ocean to oscillate at well-defined amplitude and frequency (Mellor, 1996). These oscillations are dominated by diurnal (daily) and semidiurnal (twice daily) constituents or frequencies. The tide plays an important role in a broad range of global oceanic processes including oceanic mixing and heat, salt and biogeochemical (BGC) fluxes (Vic et al., 2019;Ledwell et al., 2000;Munk, 1997), sedimentation (Kidder, ranging from space geodetic observations to the measurement of ocean currents. In this paper, we focus on the mechanical interactions of the ocean tides with the Antarctic ice shelf cavities. Needless to say, accurately representing such tidal interactions with ice shelves is indispensable in the current scenario of a warming climate, melting ice shelves and constantly evolving ice shelf geometries (De Kleermaeker et al., 2017;Blakely et al., 2022;Wang et al., 2021;Richter et al., 2022;Verfaillie et al., 2022). Ice shelves are permanent floating sheets of ice that are connected or "grounded" to a landmass (see Fig. 1). Most of the world's ice shelves hug the coast of Antarctica. However, ice shelves can also form wherever ice flows from land into cold ocean waters, including some glaciers in the Northern Hemisphere. The point at which the ice shelf loses contact with the bed is called the "grounding line". This grounding line can migrate with the tide, due to the elastic interactions between the tides and ice or the tidal flexure on the surface of an ice shelf (Brunt et al., 2010;Padman et al., 2018). The region over which the ice shelf loses contact with the bedrock and just floats in open ocean is called the "grounding zone", and in this region the ice shelves directly interact with oceans and thus ocean tides. There are a number of detailed works on ocean-ice-shelf interactions (Jacobs et al., 1992;Straneo et al., 2013;Warburton et al., 2020;Turner et al., 2017). In this work, we present the impacts of the Antarctic ice shelf cavities on the dynamics of global tides, simulated using the Los Alamos National Laboratory (LANL)led global ocean model known as the Model for Prediction Across Scales-Ocean or MPAS-Ocean (Ringler et al., 2013;Petersen et al., 2015). We have modified MPAS-Ocean to work in a barotropic or "single-layer" framework to increase computational efficiency.
The barotropic, or single-layer, global tide model within MPAS-Ocean is supplemented with specific schemes to dissipate tidal energy. The total energy of the tidal motion in the oceans across the world is about 3.5 TW (Egbert and Ray, 2003), which was originally believed to be dissipated only through bottom friction in shallow coastal shelf seas (Taylor, 1920). Thus, in the early barotropic tide models the tidal dissipation term was either linear or quadratic bottom drag in the coastal oceans, with the drag coefficients appropriately tuned to match observations (Schwiderski, 1980;Le Provost et al., 1994). However, works of Egbert and Ray (2000) showed that about a third of the barotropic tidal dissipation occurs in the deep ocean, thus debunking previous assumptions that all tidal energy is dissipated in shallow water. Works of Munk (1966) and Munk and Wunsch (1998) established that baroclinic tides contribute strongly to tidal energy dissipation. Tidal currents flowing over topography in a stratified ocean can give rise to tidal-period oscillations, known as internal or baroclinic tides. In fact, energy is transferred from the barotropic to the baroclinic tide in rough topography (Munk and Wunsch, 1998). The energy conversion from barotropic to the baroclinic tide in rough topography is rep-resented through a linear wave drag term (Stigebrandt, 1999;Jayne and St. Laurent, 2001). Such wave drags are commonly known as "internal wave drag" in the literature (Egbert et al., 2004;Arbic et al., 2004) and have been shown to improve tidal elevations by reducing root-mean-square elevation error with observations (Egbert et al., 2004;Arbic et al., 2004;Green and Nycander, 2013;Lyard et al., 2006). The internal wave drag scheme can be implemented with or without a tunable parameter. In our work, which is an early effort at introducing a linear wave drag scheme in a tide model implemented in MPAS-Ocean, we use a scalar wave drag formulation with a tunable parameter. Our scheme, originally implemented by Jayne and St. Laurent (2001), is based on a linear scaling relationship and strongly depends on a tunable parameter. As we show later, our tuning parameter essentially sets the length scale associated with dissipation.
Some drag schemes (Egbert et al., 2004) are derived using linear transformations of the bottom topography (Bell, 1975) and do not have a tunable parameter. Such schemes have been shown to perform better when applied to barotropic models, with global dissipation rates close to TPXO7.2 (Green and Nycander, 2013). However, the drawback of schemes without a tunable parameter is that the dissipation rates and root-mean-square errors strongly depend on the global bathymetry and stratification databases and thus do not guarantee an optimal prediction of elevation and dissipation rates. For example, bathymetric databases might not resolve all abyssal features. The analytically computed Nycander scheme has been shown to increase regional and global dissipation rates due to the inclusion of abyssal hill roughness on ocean spreading ridges (Melet et al., 2013). On the other hand, high-resolution bathymetric databases increase the linear wave drag strength indicating a need for more tuning (Nycander, 2005;Zilberman et al., 2009). The reason for such an increase in wave drag has been attributed to a supercritical topography (Nikurashin and Ferrari, 2011;Scott et al., 2011), but a tunable parameter is still needed to correct the effect.
The particular scalar drag scheme that we utilize in this study is the Jayne and St. Laurent (2001) scheme. This drag acts on a barotropic global ocean model forced with five major tidal constituents: M 2 , S 2 , N 2 , K 1 and O 1 . We show that our model performs very well with respect to rootmean-square tidal elevation comparisons against TPXO8 data (https://doi.org/10.5281/zenodo.7084790; Pal, 2023). We present a sensitivity study using meshes of varying resolution and show that the root-mean-square error values improve when we use high-resolution meshes. In the last section of the paper, we focus on the mechanical interactions of the ocean tides with the Antarctic ice shelf cavities and show that inclusion of ice shelf cavities improves the tidal elevation comparison results against TPXO8 data.
One other important factor that must be accounted for in global tide models is the self-attraction and loading (SAL). SAL accounts for a combination of effects: the deforma-tion of the Earth's crust due to mass loading and the selfgravitation of the load-deformed Earth as well as of the ocean tide itself (Hendershott, 1972). Self-attraction and loading can change tidal amplitudes to first order (up to 20 % in some regions) and also significantly impact tidal phases and amphidromic points (Gordeev et al., 1977). A full calculation of SAL calls for the convolution of tidal elevation with a proper Green's function or a multiplication with load Love numbers in the spectral -i.e., spherical harmonic -domain (Ray, 1998). However, due to computational efficiency, we have used a scalar SAL term in the major part of this study (Accad and Pekeris, 1978). A full calculation of SAL has been recently investigated in detail in Barton et al. (2023) and Brus et al. (2023).
The paper is organized as follows. In Sect. 2, we present a brief overview of our model components, with a focus on new modifications and developments within the MPAS-Ocean framework. In Sect. 3, we show some verification and validation results for our barotropic single-layer MPAS-Ocean ocean framework against analytical solutions from a planar two-dimensional ocean test case. We show that our barotropic ocean model reproduces known theoretical results. In Sect. 4, we present some background on TPXO8 data and the MPAS-Ocean mesh used for the simulations. In Sect. 5, we present the MPAS-Ocean model's sea surface elevation root-mean-square error results against observational TPXO8 data, the sensitivity of the results to different meshes, the impact of adding Antarctic ice shelf cavities, and the effect of adding a full (inline) self-attraction and loading term in our simulations. We finish our paper with discussions, conclusions and future directions in Sect. 6.

Model formulation and equations
This work is a systematic analysis of tides implemented in MPAS-Ocean, including the addition of ice shelf cavities. MPAS-Ocean is the ocean component of the DOE's Energy Exascale Earth System Model (E3SM) (Golaz et al., 2019;Petersen et al., 2019) (we use E3SM version V2 here). MPAS-Ocean is solved using finite-volume methods (FVMs) on an unstructured grid, and this is where our work has an extra edge. The MPAS-Ocean code has the ability to run on unstructured meshes, which helps resolve sharp bathymetry critical for tides. This multi-resolution approach is built upon two key components: a variable-resolution mesh with exceptional mesh-quality characteristics and a finite-volume method that maintains all of its conservation properties even when implemented on a highly nonuniform grid. The unstructured variable-resolution mesh is required to resolve coastal processes in greater detail, while maintaining a relatively low resolution for the open-ocean processes which are relatively simpler to parameterize. Such variable resolution is achieved using the spherical centroidal Voronoi tessellations or SCVTs (Ringler et al., 2008). Numerical algorithms specifically designed for these grids guarantee that volume and energy are conserved (Ringler et al., 2010). The timestepping scheme is the fourth-order Runge-Kutta method.
We present results of global tidal evolution over a number of horizontal mesh configurations, ranging from uniform 64 to 8 km configurations as well as variable-resolution designs. Such single-layer framework is ideal for a seamless implementation and testing of new features (e.g., the topographic wave drag and bottom friction schemes that we use). The governing momentum equations for MPAS-O that we use in our study are the shallow-water equations, written in so-called vector-invariant form (Ringler et al., 2010): p(x, y, z) = p s (x, y) + ρ 0 gh, where u represents the depth-averaged horizontal velocity, t is the time coordinate, f is the Coriolis parameter, k is the vertical unit vector, K = |u| 2 /2 is the kinetic energy, g is the gravitational acceleration constant, η is the sea surface height relative to the unperturbed ocean surface, η EQ is the equilibrium tide, η SAL is the perturbation of tidal elevations due to SAL, χ is a tunable scalar dimensionless wave drag coefficient, C h is a topographic wave drag timescale, H is the resting depth of the ocean, C D is the bottom friction coefficient, and h is the total ocean thickness such that H + η = h. The full form of the drag terms in Eq. (1) would use the total thickness h, but our implementation uses the linearized version with the resting depth H . For the linear drag scheme, we use the formulation due to Jayne and St. Laurent (2001), with the inverse timescale H C computed according to Buijsman et al. (2016). Further details are provided in Sect. 2.3. Here p s represents the surface pressure on the ocean due to the floating ice shelf cavities (see Sect. 2.5), and p refers to the total pressure on the vertical column of water.
We have the following important modification to MPAS-Ocean in our model: (a) we introduce an external forcing due to astronomical tides (Barton et al., 2023;Brus et al., 2023); (b) we introduce a new drag scheme due to friction with the seabed representing the log law of wall turbulence and a new drag scheme called topographic wave drag due to scattering with rough topography; (c) dynamics due to Antarctic ice shelf cavities, together with a full inline SAL.

Tidal forcing
The tidal forcing is obtained from an astronomical potential among celestial bodies (see Barton et al., 2023, for details).
Here we just state the governing equations for the sea surface height perturbation due to tides. The perturbation to the sea surface height η EQ is a combination of principal semidiurnal and diurnal tides: where these terms are valid for semidiurnal (sd) and diurnal (d) tidal constituents (Arbic et al., 2018). The total forcing comes from summing over each of the constituents, c. Here A and ω are the forcing amplitude and frequency, respectively, dependent on the tidal constituent; t ref is a specified reference time; t is time; φ is latitude; λ is longitude; is an astronomical argument dependent on tidal constituent; and f (t ref ) and ν(t ref ) are amplitude and phase nodal factors accounting for small known astronomical modulations in the tidal forcing. L = 1 + k 2 − h 2 represents the body tide Love numbers accounting for changes in the gravitational potential (k 2 ) due to deformation of the Earth's crust and mantle from tidal forcing (h 2 ).

Drag based on ocean bottom roughness
An ocean bottom-depth-based friction coefficient is important when we want to use a smaller value in deep water (for non-tidal energetic tuning, say) and a high value for tides in shallow regions. We implement a bottom friction coefficient based on the ocean bottom roughness, in the MPAS-Ocean model code. The friction coefficient C D in Eq. (1) is evaluated according to the following formula (Oey, 2006): where κ = 0.4 is the von Kármán constant (Von Kármán, 1931), z 0 is the roughness parameter and z b is the ocean layer thickness (i.e., the sea surface height + the ocean bottom depth). C D,min is the minimum quadratic bottom friction coefficient. Typically, in previous studies, z 0 = 10 mm and C D,min = 0.0025. The derivation of C D is given in Oey (2006) and is obtained by matching the near-bottom modeled velocity to the law of the wall (Schlichting and Gersten, 2016).

Topographic wave drag
As we discuss in the Introduction, we utilize the Jayne and St. Laurent (2001) scheme to calculate the topographic wave drag in this study. This is a scalar scheme, in which the wave drag strength is independent of the flow direction. Works by Egbert et al. (2004) show that the scalar scheme produces reasonable agreement with observations after being appropriately tuned. This scalar scheme is easier to implement than the more sophisticated tensor-based scheme (Nycander, 2005). In the tensor-based scheme, the tensor components are functions of the directionality of topographic roughness, and hence, the wave drag strength depends on the flow direction relative to the rough topography.
The term χ Cu H represents the topographic wave drag on the velocity u. The factor C H is an inverse timescale, and χ is a tunable drag coefficient. As we discussed in the Introduction, this is a scalar linear wave drag scheme, or the Jayne and St. Laurent formulation for wave drag (Jayne and St. Laurent, 2001;Buijsman et al., 2015): HereĤ is the bottom roughness, N b is the buoyancy frequency at the bottom, and L is the wavelength of the bottom topography typically set as 10 km (Jayne and St. Laurent, 2001;Buijsman et al., 2015). We implement the wave drag term as χCu representing a linear drag on the velocity u. The factor C (=Ĥ 2 N b ) is a data field we obtain following the procedure described in Buijsman et al. (2015), and χ is a tunable drag coefficient. It is interesting to note here that after tuning, the optimum value of χ settles at a value between 1.08 and 1.44 for the high-resolution meshes, which is in the same range as the pre-factor π L in Eq. (7).

Self-attraction and loading (SAL)
In the first part of the work, SAL is implemented in MPAS-Ocean via the scalar approximation (Accad and Pekeris, 1978;Ray, 1998): where η is the sea surface height prior to alterations and β = 0.09 is a scalar parameter used to approximate the influence of SAL. This approximation is a computationally inexpensive method that is sufficiently accurate for many cases. However, this does not capture the spatial dependence and large-scale smoothing of the full calculation.
To capture the full spatial dependence, in the latter part of the paper, SAL is implemented as additional body force via the sea surface height gradient term in Eq. (1). As detailed in Barton et al. (2023), we express the inline SAL for tides in terms of the spherical harmonic decomposition of the sea surface height (Hendershott, 1972): where each spherical harmonic sea surface height term η n is multiplied by a scalar coefficient (Wang et al., 2012).
Here ρ 0 is the average density of seawater, ρ earth is the average density of the solid Earth, and the multiplicative term (1+k n −h n ) represents load Love numbers corresponding to physical effects of SAL. The "1", k n and h n terms account for gravitational self-attraction of the ocean, gravitational selfattraction of the deformed solid Earth, and deformation due to loading of the solid Earth, respectively. However, the usage of sea surface height for calculating SAL is only appropriate for tides and wind-driven barotropic motions. For other motions one must use bottom pressure anomalies.

Ice shelf cavities and their connection to pressure
Within MPAS-Ocean the ocean domain is extended to include cavities under Antarctic ice shelves (Fig. 1, Jeong et al., 2020;Comeau et al., 2022). In our work, we do not study the dynamic ice-sheet-ocean coupling, so the ice-shelf-ocean interface is essentially static in time, adjusting only to the relatively small changes in dynamic pressure from the ocean. We use a constant density, salinity and temperature in our study and thus do not evolve these tracers separately. Since we have a single vertical ocean layer, we use a hydrostatic formulation for pressure with a linear dependency on the height of the vertical column (see Eq. 3). The implementation of ice shelf cavities require introducing an extra pressure term to Eq. (1) to account for the depression of sea surface height below the Antarctic ice sheet and above the continental shelf (see Fig. 1). This extra pressure p s represents the weight of floating ice shelves. The vertical coordinate below ice shelves in multi-layer MPAS-O is described in Comeau et al. (2022). For the single-layer version, the pressure is simply computed as p s = ρgD, where D is the ice shelf draft below the open-ocean sea surface datum and ρ is the density of the displaced water (Fig. 1b).

Idealized experiments
In this section we show that Eqs. (1)-(3), with the tidal forcing set to zero, reproduce known theoretical results for the ideal test cases: (a) the aquaplanet gravity wave test case and (b) Stommel's solution on a planar ocean basin. These two test cases are used to verify to operation of the single-layer shallow-water equations on the sphere, including the spatial discretization of the advection, Coriolis and pressure gradient terms and the time-stepping scheme.
An aquaplanet is a planet without any continents. If initialized with a Gaussian hump sea surface height test case, the Gaussian hump spreads at the speed of gravity waves in the ocean. The speed of gravity waves would be c = √ gH , where g is the acceleration due to gravity = 9.8 m s −2 and H is the bottom depth of the ocean (here 1000 m). We have verified that our MPAS-Ocean simulations do indeed produce a gravity wave speed of c = √ gH (details in Sect. A in the Appendix).
Next we check the Stommel solutions for a planar, flatbottomed ocean basin. Stommel assumed a rectangular ocean with the origin of a Cartesian coordinate system at the south-west corner. The y axis points northward; the x axis is eastward. The shores of the ocean are at x = 0, λ and y = 0, b. The values of these variables are listed in Table B1 in the Appendix. The ocean is regarded as a homogeneous layer of constant depth D when at rest. We observe that MPAS-Ocean simulation results agree very well with the analytical functions derived originally in Stommel (1948).
The system is initialized with a sinusoidal wind forcing as shown in the Appendix, Fig. B1b, and we simulate Eq. (1) for different values of Coriolis parameter. Below we show simulation results and also the corresponding analytical solutions, for a typical f plane test case, in which the Coriolis parameter is a constant (f = 2.5×10 −4 ). The comparison of the MPAS-Ocean simulation results with the analytical results is shown in Fig. 2. A table with the system parameters is provided in Table B1 in the Appendix.
4 Global tidal dynamics: simulation details

Icosahedral mesh design
The horizontal meshes used throughout this study consist of quasi-uniform spherical centroidal Voronoi tessellations (SCVTs) based on a semi-structured icosahedral decomposition of the sphere. The coarsest grid is composed of 12 pentagonal Voronoi cells (primal grid) and an underlying dual triangular grid (dual grid). Refinements of this mesh are obtained by incremental bisection of the spherical triangle edges and application of an SCVT optimization procedure, which iteratively re-positions triangle vertices to lie at the center of mass of their associated Voronoi cells (Ju et al., 2011). The computational grids used in this study were generated using the JIGSAW unstructured meshing library (Engwirda, 2017). To facilitate mesh convergence analysis, we employ a set of four quasi-uniform, icosahedraltype meshes in this study, with between 7 and 10 levels of structured refinement. The icosahedron 7 configuration results in approximately 60 km global mesh spacing, while the highest-resolution icosahedron 10 configuration leads to approximately 8 km global resolution.

Variable-resolution mesh design
In addition to quasi-uniform icosahedral configurations, we also explore the use of unstructured, variable-resolution meshes to enhance the representation of coastal features and regions of high tidal dissipation through selective mesh refinement. Following previous work of Pringle et al. (2021) and consistent with Barton et al. (2023) and Brus et al. (2023), we employ a 45-to-5 km variable-resolution configuration with mesh spacing l(x) tailored to resolve barotropic wavelengths and sharp bathymetric gradients: l * (x) = max min l wav (x), l slp (x), l max , l min , Here, l wav (x) and l slp (x) are barotropic tidal length-scale heuristics, with l wav increasing mesh resolution in shallow regions to resolve the wavelength of shallow-water dynamics and l slp increasing mesh resolution in areas of sharp bathymetry. β wav and β slp are tunable "mesh-spacing" parameters, set to β wav = 1 80 and β slp = 1 4 in this study. To control grid-scale noise,H and∇ H are Gaussian-filtered (σ = 1 2 ) depths and gradients obtained from the GEBCO2021 bathymetry (GEBCO Compilation Group, 2020). l * (x) is a combined estimate of mesh spacing throughout the domain, clipping l wav , l slp to l min = 5 km and l max = 45 km globally. To control the gradation of the mesh overall, l * (x) is limited to enforce a sufficiently slow growth in length scale. We use γ = 1 8 in this study.

TPXO8 and tidal evaluation
To assess the accuracy of the implemented of tidal forcing, we compare the amplitude and phase of the tidal constituents obtained from the MPAS-Ocean simulations against those of data-assimilated TPXO8. We use a number of different meshes to assess their impact on the accuracy of our results. The meshes used are as follows: (1) a globally uniform icosahedron 7 mesh with a low resolution of 64 km; (2) a globally uniform icosahedron 8 mesh with a medium resolution 30 km; (3) a globally uniform icosahedron 9 mesh with a high resolution of 15 km; (4) a globally uniform icosahedron 10 mesh with an ultra-high resolution of 8 km; (5) a variable-resolution mesh ranging from 45 to 5 km. A singlelayer global ocean configuration was initialized for each of the above meshes, and without any external wind forcing. An initial ramp-up time of 15 d was allowed for the tidal potential forcing, after which we observe a steady oscillation in the global mean kinetic energy. We compare the five major tidal constituents, of which three are semidiurnal, i.e., M 2 , S 2 , N 2 , and two are diurnal, i.e., K 1 , O 1 as obtained from our MPAS-Ocean simulations against observational TPXO8 data. The mean harmonic decomposition of each tidal constituent over a period of 90 d was calculated from the MPAS-Ocean simulations. As a comparison metric, we choose the global complex root mean square error (RMSE) defined as follows: where A is the tidal amplitude, θ is the phase lag, and the subscripts "o" and "m" refer to the observed and modeled values, respectively. The amplitude RMSE is defined as the limit of RMSE when θ o = θ m , i.e., The quantity RMSE A is weighted by the area dA of each cell. In our study, we perform the tidal analysis calculations for the M 2 constituent of the tide only.

Results
This section is organized in the following way: in Sect. 5.1 we show the tuning of the linear wave drag coefficient, χ, for different mesh resolutions. Thereafter, in the same section we show the impact of adding ice shelf cavities to the global tidal analysis results (in particular the amplitude and phase of the M 2 tidal constituent) on the different meshes. In Sect. 5.2 we show the impact of adding inline SAL (see Barton et al., 2023, for details) to the MPAS-Ocean simulations on a variable-resolution mesh. In Sect. 5.3 we show comparisons of MPAS-Ocean simulation results against data obtained from tide gauge measurements. In Sect. 5.4 we compare the results from the different MPAS-Ocean meshes in bar chart plots. Finally, we wrap up our paper with some discussions and conclusions in Sect. 6. Figure 3. Tuning of the topographic wave drag coefficient χ for runs Icos8 n (solid red with crosses) and Icos8 is (dashed red with crosses), Icos9 n (solid green with circles) and Icos9 is (dashed green with circles), and Icos10 n (solid blue with triangles) and Icos10 is (dashed blue with triangles). See Table 1 for reference.

TPXO8 comparisons and wave drag tuning
Due to the presence of a tunable parameter, χ, in the internal wave drag formulation, we carry out a series of simulations on each mesh configuration (icosahedron 8, 9 and 10) to find the optimum χ which gives the best agreement against TPXO8 data, i.e., the lowest value of RMSE A (see Sect. 4.3 for details). These simulations were carried out with and without adding ice shelf cavities. From these runs, we conclude that the optimum χ lies between 0.5 and 2.5 (see Fig. 3). The icosahedron N mesh is henceforth referred to as IcosN is and IcosN n for cases with and without ice shelf, respectively (here N = 7, 8, 9, 10). Note that the optimum value of χ for Icos9 is (χ = 1.08) and Icos9 n (χ = 1.44) are different. Similarly for the Icos10 is and Icos10 n meshes. In particular, the optimum value of χ is lower when there are ice shelf cavities in the simulation. This fact suggests that ice shelf cavities are responsible for the extra dissipation arising from a lower wave drag in the ocean. The lowest error obtained from these icosahedron mesh simulations is 0.088 m (Icos10 is ). It must be noted that we have used a scalar SAL term in these simulations. The effects of tuning χ on the spatially varying RMSE of the M 2 constituent, for Icos8 n,is , Icos9 n,is and Icos10 n,is , are shown in Fig. 4b-g. Evidently, RMSE A improves when ice shelf cavities are added in the simulation. In particular, the global RMSE A decreases from 8.77 cm (Fig. 4b) to 8.3 cm (Fig. 4c) for the Icos8 n,is meshes. The global RMSE A decreases from 8.0 cm (Fig. 4d) to 7.6 cm (Fig. 4e) for the Icos9 n,is meshes and similarly from 7.9 cm (Fig. 4f) to 7.4 cm (Fig. 4g) for the Icos10 n,is meshes.
As Fig. 4 shows, we perform two simulations for each mesh resolution: one with no Antarctic ice shelf cavities and one with Antarctic ice shelf cavities. As is apparent from a visual comparison of Fig. 4f and g, the Antarctic ice shelf cavi-    ties have the impact of reducing the errors along the Antarctic coastline. For a quantitative estimate of the reduction in errors, in Fig. 5a we show the difference in complex RMSE as obtained from Icos10 n and Icos10 is . Figure 5b shows a South Polar projection of the difference in RMSE between Icos10 n and Icos10 is .

Effect of ice shelf cavities and SAL
In this subsection we show that adding inline SAL (described in Sect. 2.4) improves the quality of MPAS-Ocean simula-tions substantially. In addition, the RMSE values show further improvement when we add Antarctic ice shelf cavities to our simulation. We use a variable-resolution horizontal mesh for this part of our work. We show the M 2 RMSE plots for the MPAS-Ocean simulations without and with ice (VR n and VR is ) in Fig. 6. The global deep water RMSE is 0.053 m for VR45-5 n , which improves to 0.044 m (VR45-5 is ) when ice shelf cavities are added. The compounded effect of the inline SAL, ice shelf cavity forcing and the variable-resolution mesh results in the improvement in RMSE values for these two simulations. This is the lowest error achieved using the current stage of development of tides within MPAS-Ocean. As discussed later in this paper, state-of-the-art simulations by Schindelegger et al. (2018) and Shihora et al. (2022) have errors of 0.044 and 0.034 m, which are very close to the RMSE values in this study. We have shown that carefully tuned internal tide parameter, inline-SAL, ice shelves and a locally high-resolution mesh are all necessary to obtain accurate, high-quality tides simulations.

Tide gauge comparison
To further investigate the impact of ice shelf cavities on tidal dynamics, we compare the results of MPAS-Ocean to tide gauge datasets including the "ground truth" stations from the University of Hawaii Sea Level Center (UHSLC). The locations of these tide gauges are shown in Fig. 7e. Tidal harmonic data at these stations were consolidated by Pringle (2019). These data include directly provided tidal harmonics or those derived from using the Matlab Function for Unified Tidal Analysis and Prediction (UTIDE) (Codiga, 2011) on time level histories. Figure 7a-d show the model vs. tide gauge amplitudes and phases for runs VR n and VR is . These plots include gauges at latitudes south of 60 • S. Figure 7a and b show scatterplots of the mean amplitude and phase of the tidal M 2 constituent for the run VR n (i.e., no ice shelf cavity case). Similarly, Fig. 7c and d show scatterplots of the mean amplitude and phase of the tidal M 2 constituent for the run VR is (i.e., with ice shelf cavity case). For the phase data, we shifted the values so that the phase differences were all within 180 • . The RMSE when comparing against the 21 UH-SLC stations is 0.07 m for the variable-resolution mesh and 0.085 m for the 8 km quasi-uniform mesh (Table 1), which is consistent with the results seen from the TPXO8 comparison of 0.06 and 0.088 m, respectively, for the Southern Ocean.
We show values of correlation coefficient or the coefficient of determination R 2 defined as follows: Here x denotes M 2 amplitude (in Fig. 7a and c) and M 2 phase (in Fig. 7b and d) from MPAS-Ocean simulation data. y denotes M 2 amplitude (in Fig. 7a and c) and M 2 phase (in Fig. 7b and d) from tide gauge measurements. When there is no ice shelf in the simulation, the R 2 values of the amplitude and phases are 0.8 and 0.83, respectively. Including ice shelf cavities increases the R 2 values to 0.955 and 0.96, respectively. The ice shelf cavities, although static, provide an accurate boundary condition along the Antarctic coastline, and thus the tidal errors reduce appreciably in the Southern Ocean. The locations of the stations from which we record the tide gauge data are shown in Fig. 7e. Figure 7e shows a South Polar projection of the Southern Ocean, in particular, ocean between latitudes −90 • and −60 • .

Summary of all runs
Finally, Fig. 8 shows three bar chart plots exhibiting how ice shelf cavity forcing improves the RMSE in the different regions of the global ocean. In particular, the bar charts are in three categories: (a) the global ocean; (b) the Southern Ocean only (latitudes south of 60 • S); (c) ocean for latitudes north of 60 • S. Although the RMSE shows improvement with ice shelf cavities for all three categories, the strongest improvement is observed for the Southern Ocean (i.e., latitudes south of 60 • S. This result further confirms that Antarctic ice shelf cavities are indeed required to accurately capture tides in an Earth system model. The lowest error achieved is on the 45 to 5 km variableresolutions mesh, with inline SAL and ice shelf cavities (VR is with a deep M 2 RMSE of 0.044 m). As a point of comparison, Schindelegger et al. (2018) and Shihora et al. (2022) both included full inline SAL calculations in a barotropic tide model and found deep-ocean M 2 RMSEs of 0.044 and 0.034 m, respectively. M 2 RMSEs with the advanced circulation model ADCIRC were found to be 0.0287 m by Pringle et al. (2021), further lowered to 0.0194 m by Blakely et al. (2022). All four of the previous studies used a more optimized wave drag, and the last study implemented SAL by reading in SAL from a data-assimilated model. Stammer et al. (2014) includes a comparison of errors for various purely hydrodynamical, non-data assimilative models ranging from 0.0525-0.0776 m.

Conclusion
In this paper we have shown the effect of adding tides in a global ocean model and have explored the effects of drag and ice shelf cavities on global tides. We have computed the tidal analysis on meshes of different resolutions, in particular, the icosahedron meshes and also a variable-resolution mesh. In particular, we use a linear wave drag scheme, and show that a systematic tuning of the wave drag coefficient, χ, results in very accurate tides. We also found that the optimum value of χ is lower (see Fig. 3) when there are ice shelf cavities in the simulation, indicating the ice shelf cavities provide the extra dissipation. The detailed analysis of simulations with and without ice shelf cavities indicates that introducing static ice shelf cavities reduces the RMSEs appreciably along the Antarctic coastline (see Figs. 4 and 5). Finally, we have shown that using an inline self-attraction and loading term, along with ice shelf cavities, results in very accurate tides, as shown by comparison against observational TPXO8 data (see Fig. 6). Our results are further validated by comparison against tide gauge datasets (Fig. 7).
The computation of the full SAL term has been explored in detail in Barton et al. (2023) and Brus et al. (2023). In those works, the full SAL has been shown to predict better tides compared to a constant SAL coefficient. Although a detailed discussion of SAL is beyond the scope of this paper, we have shown (using one reference simulation) that inline SAL, along with ice shelf cavity forcing, reduces tidal errors appreciably.
It is clear from this work that ice shelf cavities improve tidal predictions. However, the ice shelf cavities have a fixed geometry in this work. This paper paves the way for future investigations of impacts of varying ice shelf cavity geometries on global tides. Since tides also affect the melting and freezing of ice shelf cavities as shown in Rosier et al. (2014), a two-way interaction between tides and ice shelf cavities will potentially improve tidal predictions. In this work we have primarily focused on tidal errors in the open ocean near the ice shelf cavities. However, a more detailed analysis of tides and tidal pressures near the grounding line is possible as shown by Begeman et al. (2020). Also, investigations of sea ice dynamics and floating ice shelves are expected to improve tidal errors (Sun et al., 2022;Lei et al., 2017). We would like to conduct more detailed investigations of these factors in a future effort. The gravity wave speed on an aquaplanet domain, with no continents, is used to validate the barotropic model. If initialized with a Gaussian hump sea surface height test case (see Fig. A1a), the Gaussian hump spreads at the speed of gravity waves in the ocean (Fig. A1b-d represent the Gaussian profiles at times t = 3.7, 6.1, 10 s, respectively). We calculate the velocity of propagation by tracking the position of a point at the edge of the circular patch in Fig. A1a-d. In particular, we track the maximum longitudinal coordinate where the sea surface height has a value 0.2 m. The speed of gravity waves would be c = √ gH , where g is the acceleration due to gravity = 9.8 m s −2 and H is the bottom depth of the ocean (here 1000 m). Thus the gravity wave speed is ∼ 99.8 m s −1 . The formula for gravity waves is the solution to the linear wave equation. Since we are running the MPAS-Ocean model with the nonlinear term, the speed will be slightly lower as shown in Fig. B1. The simulation parameters are shown in Table A1.

Appendix B: Stommel basin test case
Next we check the Stommel solutions for a planar, flatbottomed ocean basin. Stommel assumed a rectangular ocean with the origin of a Cartesian coordinate system at the southwest corner. The y axis points northward; the x axis is eastward. The shores of the ocean are at x = 0, λ and y = 0, b. The ocean is regarded as a homogeneous layer of constant depth D when at rest. Below we show the equations of motion according to Stommel (1948).
Here f is the Coriolis parameter, D is the bottom depth (200 m in our test case), h is the height above the sea surface, u is the x component of velocity, v is the y component of velocity and R is the Rayleigh friction coefficient. Equations (B1)-(B2) are the steady-state Navier-Stokes equations. In Stommel (1948), Stommel derived analytical expressions for the sea surface height and the horizontal and vertical velocities for a planar domain and for different values of the Coriolis parameter f . The analytical solutions can be exactly verified against the MPAS-Ocean model solutions in the steady state, provided the domain configurations are the same as in Stommel (1948). Thus in MPAS, we choose a rectangular domain with the horizontal edges at x min , x max so that λ = x max − x min and the vertical edges at y min and y max so that b = y max − y min .
Also, the coefficients in Eqs. The external force is provided by wind stress, which, for the sake of simplicity we consider to be a regular periodic function (Stommel, 1948). The wind force is given by F = F 0 cos πy b . This simple functional form of the external force allows evaluation of analytical solution of Eqs. (B1)-(B2). The wind forcing amplitude is F 0 = 1 N m 2 in MPAS, which when divided by the density ρ gives the amplitude force per unit mass as F 0 ρ = 1 × 10 −3 m 2 s −2 . In Table B1 we give the different parameters, which are the same for the simulations in this section. We show a visualization of the wind forcing in the rectangular domain in Fig. B1.
There is a wind forcing at t = 0, with the initial condition shown in Fig. B1. We simulate Eq. (1) for different values of Coriolis parameter and list the cases below. We also provide the corresponding analytical solutions and compare the simulation results with the analytical results. If the Coriolis parameter is 0, the zonal and meridional velocities and the sea surface heights have the following analytical solutions in the steady state: h(x, y) = −(F /gD)(pe Ax /A + qe Bx /B) − (b/π ) 2 (F /gD)(pAe Ax + qBe Bx ) p = e −πλ/b , q = 1.
The corresponding comparisons are shown in Fig. 2.
Code availability. Code for the Energy Exascale Earth System Model (E3SM) used in this study is publicly available: https://github.com/E3SM-Project/E3SM/ (last access: 10 February 2023). A listing of the code and model configurations used in this work can be found in the following archive: https://doi.org/10.5281/zenodo.7084857 (Pal et al., 2022).
Author contributions. NP, KNB, MRP, DE and SRB prepared the paper, designed and implemented the coding upgrades into MPAS-Ocean, designed and performed the experiments, and conducted the tidal analysis. NP led the study, conducted the majority of simulations and analysis, and wrote the first draft of this paper. NP and MRP contributed to including the ice shelf cavities, topographic wave drag and bottom drag subroutines within MPAS-Ocean. DE contributed to the coding and development of the JIGSAW software integral to the horizontal mesh creation in this study. KNB and SRB and MRP contributed to the inline SAL calculation software. MRP and DE provided project oversight and code review. All authors participated in discussions to formulate the plan and provide feedback and reviewed and revised the text.