MPAS‐Ocean Simulation Quality for Variable‐Resolution North American Coastal Meshes

Climate model components utilizing unstructured meshes enable variable resolution, regionally enhanced simulations within global domains. Here we investigate the relationship between mesh quality and simulation statistics using the JIGSAW unstructured meshing library and the Model for Prediction Across Scales‐Ocean (MPAS‐Ocean) with a focus on Gulf Stream dynamics. In the base configuration, the refined region employs 8 km cells that extend 400 km from the coast of North America. This coastal‐refined region is embedded within a low‐resolution global domain, with cell size varying latitudinally between 30 and 60 km. The resolution transition region between the refined region and background mesh is 600 km wide. Three sensitivity tests are conducted: (a) The quality of meshes is intentionally degraded so that horizontal cells are progressively more distorted; (b) the transition region from high to low resolution is steepened; and (c) resolution of the coastal refinement region is varied from 30 to 8 km. Overall, the ocean simulations are shown to be robust to mesh resolution and quality alterations. Meshes that are substantially degraded still produce realistic currents, with Southern Ocean transports within 0.4% and Gulf Stream transports within 12% of high‐quality mesh results. The narrowest transition case of 100 km did not produce any spurious effects. Refined regions with high‐resolution produce eddy kinetic energy and sea surface height variability that are similar to the high‐resolution reference simulation. These results provide heuristics for the design criteria of variable‐resolution climate model domains.


Introduction
Climate models based on unstructured horizontal meshes have matured in recent years. Unstructured global simulations of historical periods compare well when validated against observations and against other future climate projections (Golaz et al., 2019;Petersen et al., 2019;Scholz et al., 2019). Unstructured meshes offer great freedom in placing resolution in the areas of interest for regionally refined simulations and also suggest the possibility of improving global simulation quality with targeted areas of high resolution. However, modelers now have a dizzying array of choices to make in designing their meshes, compared to the limited

Background
There now exists a growing selection of unstructured mesh models that are used for various global and regionally focused forecasts and analyses. This includes MPAS , FESOM , ICON (Korn, 2017), FVCOM (Chen et al., 2003), SCHISM , SLIM (Kärnä et al., 2013), and Fluidity (Davies et al., 2011). Mesh creation tools such as Shingle 2.0 have been developed to produce high-quality reproducible meshes efficiently (Candy & Pietrzak, 2018). Models differ in the arrangement of variables on the underlying computational grid and in the numerical techniques employed, with both unstructured triangle-and polygon-based finite-volume and finite-element type discretization schemes adopted in various frameworks. As such, different approaches to the construction and optimization of the models' underlying unstructured meshes have been explored, including techniques based on Centroidal Voronoi Tessellation (CVTs) Yang et al., 2018), optimization via optimal transport (McRae et al., 2018;Weller et al., 2016), as well as triangulation-based refinement schemes (Lambrechts et al., 2008a;Remacle & Lambrechts, 2018). A mesh generation tool developed by Lambrechts et al. (2008b) refines according to bathymetry, bathymetry gradients, and distances from coasts. In the context of MPAS-Ocean, the numerical scheme requires that the mesh define a highly regular, orthogonal 10.1029/2019MS001848 tessellation, constraining grid generation choice to algorithms that can generate optimized Voronoi-type meshes .
Variable resolution is advantageous in situations where highlighting a region may help to correct a bias or resolve a dynamic condition. In many cases, the resolved region will also be the focus of the investigation, but resolution can also be placed to correct a bias that is impacting a global simulation. Variable resolution can serve as a replacement for nested grids, with the advantage that variable resolution can be applied in more complex configurations and is more integrated with the global simulation (Biastoch et al., 2018;Hagos et al., 2013). Nested grids have the advantages of more easily implemented variable time stepping and simplified grid geometry. However, nesting introduces challenges with conservation, coupling, interpolation, and noise control (Debreu & Blayo, 2008).
Meshes in which the resolution varies as a function of latitude have been used to compensate for the changing Rossby radius with latitude. This approach is used in the standard high-resolution MPAS mesh . Variable-resolution meshes are designed to improve the dynamics of a particular region or process and also to provide good global dynamics. Variable resolution meshes may refine particular regions, for example, the Arctic Ocean  or a coastal region (Androsov et al., 2019). Variable resolution has been applied in regional ocean models to capture a wide range of scales, from tens of kilometers to tens of meters. SCHISM and FVCOM have been applied in variable-resolution cases to place high resolution in estuaries and straits, where narrow channels and complex bathymetry must be properly represented, for example, the Chesapeake Bay (Ye et al., 2018) and the Canadian Archipelago . It has also been applied for a variety of coastal processes, for example, storm surge (Fernández-Montblanc et al., 2019; and nutrient distribution (Tian et al., 2014).
Resolution can also be placed based on a particular parameter. For example, FESOM uses meshes that refine to the local Rossby radius (Sein et al., 2017), a more sophisticated approach than refining based on latitude alone. FESOM also uses meshes that refine according to eddy variability. This approach is much less computationally expensive than refining based on Rossby radius but has been shown to improve deep ocean biases and Gulf Stream separation (Rackow et al., 2019). FESOM's meshes which refine based on sea surface height (SSH) variability are useful for capturing boundary currents (Biastoch et al., 2018). A configuration which used high resolution over areas of high SSH variability, areas upstream of the separation of mid-lattitude jets, and in the Nordic Seas improved Gulf Stream separation and biases in the Northwest Corner (Sein et al., 2016).
Because of the computational cost and complexity of global simulations, the majority of variable-resolution tests have been performed on idealized or simplified domains. For example, in order to eliminate the effects of continental geography, many tests have used aquaplanet configurations (Abiodun et al., 2008;Hagos et al., 2013;Lorant & Royer, 2001;Rauscher & Ringler, 2014;Rauscher et al., 2012;Zhao et al., 2016). Others have used two-dimensional domains (Düben & Korn, 2014). These simplified domains can demonstrate the effects of mesh resolution independent of other variables. Additionally, atmospheric variable-resolution simulations can inform choices in ocean domains (Abiodun et al., 2008;Düben & Korn, 2014;Park et al., 2014;Rauscher & Ringler, 2014;Zarzycki et al., 2015;Zhao et al., 2016). However, mesh resolution and design consequences on more realistic simulations are still largely unknown, even though the use of variable resolution in realistic simulations is becoming more widespread.
While mesh design is still a developing field, the literature points to several important considerations. In the past, parameter values for sub-grid scale physics were typically tuned for each resolution. Now, for variable-resolution meshes, parameterization schemes must work well across the span of grid cell sizes. Another consideration is that variable-resolution results compared against uniform high-resolution simulations may not necessarily be comparable near mesh transition regions. For example, a current flowing from a non-eddy permitting to an eddy-permitting region may not immediately develop eddies. Instead, eddies will develop downstream of the beginning of the high-resolution region once perturbations have time to evolve (Danilov & Wang, 2015). A similar result was found in atmospheric variable-resolution aquaplanet simulations, in which precipitation error was decreased in the eastern (downstream) section of the high-resolution region but not in the western (upstream) section (Hagos et al., 2013).
A high-resolution region will also have effects on the rest of the domain. Most obviously, a high-resolution region will have an effect immediately downstream, as the increased variability of the high-resolution region is carried into the low-resolution region (Danilov & Wang, 2015). Changes to dynamics within the high-resolution region can propagate to other global processes (Hagos et al., 2013;Lorant & Royer, 2001;Sakaguchi et al., 2016;Sein et al., 2017). Conversely, the impact of the global domain on the high-resolution region is also important. A high-resolution region can decrease local error but will have a limited impact on processes that are due to causes outside the high-resolution region (Zarzycki et al., 2015).

The Model for Prediction Across Scales-Ocean (MPAS-Ocean)
The Model for Prediction Across Scales (MPAS) is an open source framework that provides common functionality for climate model components on unstructured meshes. This includes a mesh specification, decomposition of variables across processors, parallel input and output specified in a run time streams file, timers, and error handling. Finite volume operators were developed for Voronoi tesselations in Ringler et al. (2010) for the shallow water equations using mimetic methods to guarantee that mass, velocity, and potential vorticity evolve in a consistent and compatible manner.
MPAS-Ocean solves prognostic equations for momentum, thickness (volume), and tracers using these operators  and can be run using both regular and unstructured meshes on Cartesian and spherical domains. The time stepping is split-explicit, where the 2D barotropic equations are sub-cycled within 3D baroclinic time steps. Both parts use a second-order predictor-corrector method based on Higdon (2005), as detailed in Appendix A5 of Ringler et al. (2013). Advection uses the flux-corrected transport scheme (Skamarock & Gassmann, 2011), which blends high-and low-order fluxes to preserve monotonicity and is second-order accurate on variable-resolution meshes . The simulations presented here use a z-star vertical coordinate, which is the standard choice for global simulations. The MPAS-Ocean vertical coordinate is designed within an Arbitrary Lagrangian-Eulerian (ALE) framework Reckinger et al., 2015). Simulations typically include 60, 80, or 100 vertical layers, which vary from 2 m thick at the surface to 150 m thick at a depth of 5,000 m.
The vertical mixing scheme is the K-Profile Parameterization (KPP) (Van Roekel et al., 2018), calculated in the CVMix library (https://github.com/CVMix/CVMix-src, https://doi.org/10.5281/zenodo.1000800) and applied implicitly. The horizontal mesoscale eddy parameterization is Gent-McWilliams thickness advection (Gent & Mcwilliams, 1990), applied to variable-resolution meshes with a coefficient of 600 m 2 s −1 at grid cells larger than 30 km and tapering linearly to zero between 30 and 20 km. Viscosity (del-2) and hyperviscosity (del-4) are applied to the momentum equation with coefficients that depend on the grid cell size as respectively, where Δx is the horizontal grid cell width. The coefficients were tuned in Petersen et al. (2019) to be as small as possible, while ensuring that dissipation is sufficient for stability and to prevent grid scale noise. The addition of hyperviscosity removes energy more strongly at the highest wave numbers and allows a smaller viscosity coefficient. No horizontal diffusion is explicitly applied to the tracers.
For this study MPAS-Ocean was run with the same choice of parameters as typical global simulations, such as those presented in Petersen et al. (2019). One exception is that this study used the stand-alone version of MPAS-Ocean, rather than the coupled E3SM code. Stand-alone mode applies idealized, constant atmospheric forcing, where wind forcing is averaged over a 65-year CORE (coordinated ocean-ice reference experiments) cycle (Griffies et al., 2009). The choice to use stand-alone mode was made for two reasons. First, it considerably simplified the required setup, streamlining the work required for a large parameter study. Achieving realistic climatological results would require a lengthy spin-up process and longer run time, neither of which were possible for this number of simulations. Second, the idealized forcing simplified the conditions of the simulations, making it easier to evaluate any numerical effects of the meshes. Because this study used both new variable-resolution meshes and a new mesh creation tool, it was important to test simplified domains before running the meshes with the complexities inherent in fully coupled E3SM simulations.
The simulation is spun up for 1 year from an initial climatology of Polar Science Center Hydrographic Climatology, version 3 (PHC3.0, Steele et al., 2001). Surface salinity and temperature restoring to yearly averaged PHC3.0 is conducted with a piston velocity of 1.37 m day −1 to represent surface fluxes. Sea-ice is not included in these simulations. Simulations with more realistic atmospheric forcing (6-hourly CORE winds and surface fluxes) and active sea ice have been run within E3SM using the coastal-refined mesh (CUSP8) are currently underway and will be presented in a future publication.

JIGSAW Mesh Generation
JIGSAW is an unstructured meshing library designed to generate high-quality grids for computational simulation, with a focus on constructing optimized Voronoi-type grids for unstructured mesh GCM's. JIGSAW is a hybrid algorithm that combines both Delaunay-refinement and Voronoi optimization type approaches to enable the rapid generation of very high quality, high-resolution Voronoi/Delaunay meshes on the sphere. The Delaunay mesh is composed of triangular cells, and the Voronoi mesh is the dual composed of polygonal cells. A key advantage of this combined strategy is efficiency and guaranteed mesh quality. Previous mesh generation methods used in MPAS  used an iterative Lloyd's method and were extremely slow.
With JIGSAW, highly optimized, large-scale variable-resolution Voronoi-type meshes can be generated in the order of minutes, allowing model users to easily create and explore a range of alternative configurations, investigate mesh quality and resolution dependence, and tailor the overall mesh and model configuration to their simulation needs. This capability was exploited in the present study to design and assess a range of coastal-enhanced MPAS-Ocean configurations and to explore various model/mesh feedbacks.
Meshes can be generated in local two-dimensional domains and over general spheroidal surfaces. Mesh resolution can be adapted to follow complex user-defined metrics, including topographic contours, solution profiles, and/or coastal features. This flexibility enables the construction of complex, variable-resolution model configurations, offering enhanced simulation fidelity in regions of interest or importance.
Given a particular geometry definition and resolution specification, JIGSAW proceeds to assemble the unstructured mesh incrementally-first creating a conforming Delaunay triangulation of the domain using a "frontal" Delaunay-refinement strategy (Engwirda & Ivers, 2016), before optimizing the resulting Voronoi/Delaunay tessellation using Optimal Delaunay Tessellation (ODT) type techniques (Chen & Holst, 2011;Engwirda, 2017). The final mesh is guaranteed to consist of high-quality triangular and polygonal cells that facilitate a locally orthogonal unstructured C-grid staggering. The final meshes are heavily optimized, typically satisfying the stringent mesh quality requirements imposed by the TRiSK (Thuburn Ringler Skamarock Klemp) discretization scheme (Ringler et al., 2010) used in MPAS-Ocean.
For TRiSK-based schemes, a complex array of geometrical and topological constraints must be satisfied (Engwirda, 2018), requiring tessellations to be orthogonal, centroidal, well-centered, and smoothly varying. These criteria require that the vertices of the triangular and polygonal grid cells lie close to the centroids of their enclosing control-volumes, that the staggered Voronoi and Delaunay edges intersect near their midpoints, that the Delaunay triangles contain their own circumcenters, and that the cell angles and edge-lengths be "nicely" distributed with respect to the desired mesh resolution constraints. Satisfying such criteria is nontrivial, and failure to do so has been shown to impact on the asymptotic accuracy and stability of the underlying numerical scheme (Peixoto, 2016) in idealized cases.
The expected accuracy of the TRiSK formulation is thus a function of both the geometry and topology of the mesh and can be quantified by considering the nature of the discrete gradient, divergence, curl, and interpolation operators used to discretize the continuous PDE's (Engwirda, 2018;Ringler et al., 2010). Based on theoretical analysis, it is expected that the accuracy of TRiSK is maximized (achieving quasi second-order scaling) only for "perfect" tessellations consisting of regular hexagons and equilateral triangles. For general unstructured meshes incorporating irregular and/or deformed polygonal and triangular cells, numerical accuracy is expected to degrade-leading to quasi first-order behavior in many practical configurations (Peixoto, 2016). The goal of mesh optimization is to construct a tessellation that serves to minimize these numerical errors, thus maximizing the quality of the resulting simulation.
A key question in the current study is to assess what impact mesh quality has on practical MPAS-Ocean simulations and to define an associated set of "best practice" guidelines for mesh generation. To this end, an "ensemble" of meshes was considered in the current work-exploring the impact of different mesh quality perturbations and variable-resolution designs on the characteristics of spun up ocean simulations.

Meshes and Simulations
All the meshes used are based on two base configurations, a global low-resolution mesh and a mesh with refinement along the coast of North America. The global low-resolution mesh, EC60to30, varies from 30 km resolution at the equator and poles to 60 km resolution at the mid-latitudes and uses 100 vertical layers.
The base EC60to30 mesh created using JIGSAW was compared against the EC60to30-E3SM-V1 mesh created using a parallel Lloyd's algorithm , which was used in previously published E3SM simulations (Golaz et al., 2019;Petersen et al., 2019). Images of the two EC60to30 meshes can be seen in the first two panels of Figures 1 and 2, which show two different metrics for measuring cell quality. Figure 2 shows the percent change between the size of neighboring cells, and Figure 1 shows close-up images of the mesh and the ratio of the smallest to largest sides of the cells. These metrics show the different strategies used by each of the mesh creation methods. In order to cover the sphere, the mesh must deviate from regular hexagons. E3SM-V1 spreads these imperfections between large numbers of cells, resulting in smooth regions of lower-quality cells. JIGSAW concentrates the imperfections into "seams" of low-quality cells separating regions of very high quality cells.
The second base mesh is the North American refined mesh, created to investigate processes affecting North American coastal regions at high resolution while avoiding the cost of running a global high-resolution model. In addition to the improvements in the dynamics of the Gulf Stream investigated in this study, using the CUSP8 mesh will allow improved simulation of a variety of coastal processes around North America. The CUSP8 mesh (Coastal United States "Plus" with 8 km coastal resolution) has high resolution along the Atlantic and Pacific coasts from Central America to the Arctic, with additional high resolution in the Caribbean and around Greenland, Hawaii, and the Bering Strait (see Figure 3). The CUSP8 mesh is built on top of a background low-resolution EC60to30 mesh. It uses 80 vertical layers.
In the CUSP8 mesh, the transition between the high-resolution region and the background mesh begins 400 km off the coast. This was chosen so that the high-resolution region encompassed the flow of the Gulf Stream along the coast and other important coastal processes. The transition region is 600 km wide and  follows the following functions: where W is the weight, D is the distance from the coast, D start is the distance from the coast where the transition region begins, and D width is the transition width. The final cell width, C, shown in Figures 3 and  4, is simply a linear combination of the coastal and background cell widths, C coast and C back .
In addition to these two base meshes, a mesh with 8 km resolution spanning the full North Atlantic basin (NA8) was created. Like the CUSP8 mesh, it was built on a background EC60to30 mesh (see Figure 3). A global high-resolution simulation was not feasible for this study, but the NA8 mesh provides high resolution within the region of interest in the North Atlantic, providing a benchmark for the performance of the CUSP8 mesh.
All resolutions are first created on a full sphere, and then continental and island land cells are culled if the cell center is within a high-resolution coastline defined by connected points(http://www.naturalearthdata. com). The bathymetry is obtained by interpolation of the ETOPO2 2-Minute Gridded Global Relief Dataset available from the National Geophysical Data Center (section 4.1; Ringler et al., 2013). Partial bottom cells are used for a better representation of the bathymetry. All domains presented in this study use this standard method of initializing coastlines and bathymetry, with the highest resolution data available. This means that regions with higher mesh resolution also have finer coastal and depth features.  (4)) for the transition width study (left) and the convergence study (right). The background resolution plotted is 60 km; however, the background resolution varies from 30 to 60 km depending on latitude. In order to ensure that all the meshes could be compared, EC60to30 simulations were run in each vertical configuration: 60, 80, and 100 layers. All three EC60to30 meshes performed similarly in terms of kinetic energy (KE), sea surface height (SSH), eddy kinetic energy (EKE), and sea surface height root mean squared (SSH RMS) (see Figure A2).
Three studies were performed to investigate mesh features and their effects on simulation quality. All the simulations were run for 10 years, with analysis performed on the last 9 years.
The first study uses the EC60to30 mesh to examine the effect of poor mesh quality on simulations. Meshes were intentionally degraded, producing poor quality cells. Variable-resolution meshes by necessity contain distorted cells within the transition regions. This is a particular concern when designing complex meshes such as the CUSP8 mesh that have large variations in resolution and relatively narrow transition regions.
Because of the difficulty of decoupling the effects of poor cell quality from the effects of a change in resolution, the effect of poor cell quality on simulations was investigated using EC60to30 meshes with cell quality degraded globally.  A mesh degradation heuristic was developed to systematically reduce the quality of meshes, perturbing the position of vertices and updating topology to effectively "de-optimize" the overall structure of a given mesh and degrade the shape of its cells. Care was taken to ensure that degraded meshes inherited the large-scale properties of their parent grids, adhering to variations in resolution and matching cell counts exactly. The kernel of the degradation operation consisted of randomly perturbing a subset of vertices toward the centroid of their largest neighboring triangle. By controlling the magnitude of the average relative vertex perturbation, the notion of a " -degraded" mesh was introduced-a 0.5-degraded mesh would re-position vertices (on average) halfway between their current position and the neighboring centroid location. Mesh topology was updated following the re-positioning of vertices to ensure that the orthogonality of the mesh was preserved. Starting from a fully optimized initial mesh, several iterations of this process were repeated to ensure that degraded grids were sufficiently randomized.
Three degraded meshes were created, EC60to30-degraded-0.25, EC60to 30-degraded-0.50, and EC60to30-degraded-0.75, with larger degradation fractions indicating a more degraded mesh. Figure 1 shows the mesh quality of the standard EC60to30 mesh and the degraded meshes.
The second study investigates the effects of the steepness of the transition function in the CUSP8 mesh (equation (3)) by varying the transition width from 100 to 900 km (Figure 4). A 10 km transition was attempted as well but failed early in the spin-up process. This study was designed to investigate how steep the transition function could be without negatively affecting the simulation quality. In addition to exploring the steepness of the transition function, this study also investigates the impact of the size of the higher-resolution region. Because the beginning of the transition region was kept fixed, the center of the transition region and the beginning of the low-resolution region were closer to the coast for steeper transitions, effectively shrinking the higher-resolution region (see Figure 5).
The third study investigates different coastal resolutions ranging from 8 km (CUSP8) to 30 km (CUSP30) in order to explore the improvements in dynamics with increased resolution. Resolutions were chosen to span the range between the highest resolution in the MPAS high-resolution model, which varies from 8 to 16 km, and low-resolution model, which varies from 30 to 60 km. The computational performance of the meshes was also examined in order to give a better sense of the trade-off between higher resolution and higher simulation costs. These meshes were compared against the EC60to30 and NA8 meshes. Ideally, the CUSP8 mesh would show dynamics comparable to the NA8 mesh within the high-resolution region with a much lower cost than a global high-resolution mesh. This study is designed to examine the degree to which we can produce the correct dynamics on variable-resolution meshes. It is likely that the CUSP8 mesh will not be able to fully recover the dynamics of a high-resolution simulation and that further modifications to the mesh will be required. Table 1 shows the parameter values used for each simulation. These values were chosen based on the highest resolution region of the simulation. The EC60to30-degraded-0.50 and EC60to30-degraded-0.75 meshes had to be run at a smaller time step than the standard EC60to30 meshes due to the smaller cell sizes introduced by the degradation process. All simulations were started from a short initialization. This initialization process maintains stability using Rayleigh damping and small time steps after the run is started from rest. It is not intended to produce an equlibrium state. All runs used a 7 day initializations except the EC60to30-E3SM-V1, EC60to30-degraded-0.50, and EC60to30-degraded-0.75 meshes. The EC60to30-E3SM-V1 mesh used a 21 day initialization. The EC60to30-degraded-0.50 and EC60to30-degraded-0.75 meshes required longer initializations and were spun up to a different point because of the smaller time step required.

Results and Discussion
The analysis focuses on the Gulf Stream because it is the most prominent feature within the high-resolution region of the CUSP simulations. The Gulf Stream also crosses out of the high-resolution region, allowing Note. See Figure 7 for plots of the data and Figure 6 for a map of the Gulf Stream transects. Observational references: Florida-Cuba (Johns et al., 2002), Florida-Bahamas (Johns et al., 2002), Cape Hatteras (Halkin & Rossby, 1985), New Jersey (Rossby et al., 2014), Drake Passage (Donohue et al., 2016), Tasmania-Ant (Ganachaud & Wunsch, 2000), and Africa-Ant (Ganachaud & Wunsch, 2000).
the effect of the transition in resolution to be investigated. The sea surface height, kinetic energy, sea surface height root mean squared, and eddy kinetic energy were analyzed for all simulations. Transport through transects along the Gulf Stream was calculated (see Figure 6 for a map of the Gulf Stream transects). Transport through Southern Ocean transects were also calculated in order to see if the high-resolution region had an impact on global dynamics (see Table 2 and Figure 7 for the transect results). SSH RMS and EKE were averaged along the Gulf Stream region (see Figure 6). These results are not expected to closely match observations, both because of the idealized forcing used and because of the differences between the sampling techniques used to calculate observational estimates and those used in our calculations. Global analysis was also run looking at global temperature, salinity, SSH, and EKE. However, because of the extremely similar results for all simulations, this paper focuses only on analysis of the areas within and around the high-resolution region. Preliminary results from simulations with realistic climatological forcing are also used to give an indication of how CUSP meshes perform in realistic climate simulations. Further results will follow in subsequent papers.
The comparison of the JIGSAW EC60to30 mesh and the EC60to30-E3SM-V1 mesh showed that they performed very similarly, confirming that the meshes created using JIGSAW produce comparable results to those used in previous MPAS studies (see Figure A1).

Study 1: Degraded Meshes
Though the degradation factor for the degraded mesh study and the transition widths for the transition width study were chosen independently, the degraded meshes were found to be a good proxy for the transition regions (see Figure 8). The cell quality in the transition region for the 100 km transition width is comparable to the cell quality in the 0.75 degraded mesh, and the cell quality in the transition region for the 900 km transition is comparable to the cell quality in the 0.25 degraded mesh. Thus, the results of the degraded mesh study should also be considered when interpreting the results within the transition regions of the CUSP meshes.
The results of the degraded mesh study are summarized in Figure 9, showing snapshots and averaged distributions of sea surface height and kinetic energy in the CUSP region. Overall, it was found that mesh degradation did not significantly effect the quality of the simulations, with the pattern and magnitude of sea surface height and kinetic energy for the set of degraded meshes and the optimized EC60to30 mesh visually near identical. The more degraded meshes were found to have slightly higher average sea surface height variability and eddy kinetic energy (see Figure 10). Transport through all transects measured showed no significant variation between the degraded meshes (see Figure 7) and the reference EC60to30 configuration.  Table 2 shows the data and Figure 6 shows the locations of the transects.
Overall, it was not found that a reduction in mesh quality had notable adverse effects on the simulations, beyond the need for the use of smaller time steps in highly degraded cases, due to the presence of smaller grid cells. While increasing computation time, the use of smaller time steps does not impact the quality of results. An EC60to30 mesh was run with a 2 min time step (the time step used for the most degraded mesh), and the results were compared against the standard EC60to30 case. There was no apparent effect of the smaller time step. While these results are encouraging-showing that the TRiSK-based numerical formulation employed by MPAS-Ocean appears to be relatively insensitive to mesh distortion-these conclusions should be tempered by the nature of simulations run. Specifically, our analysis is restricted to relatively low-resolution, eddy parameterized configurations, where it may be expected that the dissipation due to viscous mixing acts to damp down any noise and/or oscillations generated at the grid scale. It is further noted that at low resolution and with constant forcing, energy is primarily injected into the system at relatively long, well-resolved wavelengths. Future studies may expand on the results presented here, using a set of high-resolution, eddy-resolving configurations to study interactions between the discretization, mesh quality, and grid scale response in the absence of explicit viscous damping.

Study 2: Transition Width
The analysis of the transition width study can be found in Figure 11. As the transition width increases, the dynamics of the simulations improve. The simulations with wider transition regions show greater SSH RMS and EKE (see Figure 10). This is to be expected, both because the transition is less steep, leading to higher quality cells in the transition region, and because the higher-resolution area is effectively larger with a greater transition width (note the locations of the center of the transition region in Figure 11). The three widest transitions (900, 600, and 300 km) have closer average values. The 100 km transition, where the Gulf Stream is meandering into the low-resolution region, shows a more significant decline in average SSH RMS and EKE. Even within 400 km of the coast, where the resolution is 8 km in all the simulations, the dynamics were improved by a wider transition region. With a narrower transition, meanders and eddies from the Gulf Stream cross into regions of lower resolution. It appears that these features are then smoothed out and do not have time to recover even when returning to the high-resolution region. This result is consistent with that found by Danilov and Wang (2015), in which eddies did not develop at the beginning of the eddy-permitting region but instead developed only once perturbations had developed further downstream. This is clearly seen in Figure 11. It also appears possible that the transition region is affecting the path of the Gulf Stream, "trapping" it within the high-resolution region. However, there is not a wide enough spread of transition widths in this study to say anything definitive about this effect. The transport through the Gulf Stream transects increases with wider transition widths, with the exception of the Cape Hatteras transect, which shows the opposite pattern.
In addition to the results examined here, the results of the degraded mesh study should be considered as a proxy for the transition regions. Although the cell quality within the CUSP8-transition-100 transition region is comparable to that of the EC60to30-degraded-0.75 mesh, the CUSP8-transition-100 mesh did not require the smaller time steps that the EC60to30-degraded-0.50 and EC60to30-degraded-0.75 meshes did. If the reduction in time step in the case of the degraded meshes is due to the presence of smaller cells, it would make sense that it is not required here, as the degraded cells in the transition region are larger than the cells in the high-resolution region. The results of the degraded mesh study indicate that mesh quality does Figure 10. Plot of the average surface SSH RMS and EKE over the region shown in Figure 6 for years 2-10. The CUSP meshes have significantly higher average SSH RMS and EKE than the EC60to30 meshes. The variability increased as the mesh degradation increased. As the transition width was narrowed, the variability decreased, though this effect was small between CUSP8-transition-900 and CUSP8-transition-300. As the resolution decreased, the variability decreased, reaching the same values as the EC60to30 mesh for CUSP30, as would be expected.
not have a large impact on simulation results. The variation between the CUSP meshes is probably due primarily to other effects, such as the smaller region of higher resolution, rather than the cell quality within the transition region.

Study 3: Coastal Resolution
The analysis of the coastal resolution study can be found in Figure 12. The meshes with higher coastal resolution showed significantly improved dynamics, particularly in eddy kinetic energy and sea surface height variability, which were much smaller in CUSP30 (see Figure 12). The Gulf Stream within the high-resolution region in CUSP8 is similar to that of NA8. However, as noted in the transition width study, features that cross into the lower resolution transition region and then back into the high-resolution region, such as meanders and eddies, are less well resolved in the CUSP8 simulation. Figure 10 shows the effect of this as well.
The integration region is within the high-resolution region in all of the simulation. However, the average SSH RMS and EKE decreases with decreasing transition width. This indicates that the variability within the high-resolution region is affected by the adjacent low-resolution region, presumably by current passing from low to high resolution. This effect was described by Danilov and Wang (2015). Figure 13 shows the path of the Gulf Stream in the coastal resolution study. The NA8 simulation shows very little variability in the path of the Gulf Stream, while the CUSP8, CUSP12, and CUSP20 simulations show much more. The CUSP30 simulation also does not show much variability in the Gulf Stream path, but this is expected as the resolution is too low to be eddy permitting. The variation in the Gulf Stream path is also apparent in the transport through the transects along the Gulf Stream. In the southernmost transects (Florida-Cuba and Florida-Bahamas) where the flow is geographically constrained, the transport increases with increased resolution. The Cape Hatteras and New Jersey transects do not show this pattern. Figure 13 and Table 2 show that in the CUSP8, CUSP12, and CUSP20 simulations, there is significant variability in the path of the flow in the region of the New Jersey transect. Periods of very low or negative transport are probably due to North-South flow through the transect as the Gulf Stream separates from the coast. This can be seen in some of the monthly Gulf Stream paths seen in Figure 13.
This high variability is not due only to high resolution, as the NA8 simulation, which has the same coastal resolution, shows very little variability in the path of the Gulf Stream. The very low variability in the NA8 simulation is probably due largely to the idealized forcing used, as this effect did not show up in global high-resolution simulations with realistic climatological forcing. Initial results from global high-resolution simulations with climatological forcing show that the Gulf Stream had a realistic path and variability (see Figure A3). It appears that the lack of variation in the forcing or in the mesh itself prevents the NA8 mesh from developing meandering features. Variable forcing appears to resolve this problem. The variability in cell size and quality in the CUSP meshes may allow these features to develop. For all the transects, the variability increased significantly between CUSP30 and the higher-resolution meshes, as would be expected when transitioning to an eddy-permitting resolution.
The CUSP simulations are also a significant improvement on global high resolution in terms of cost (see Figure 14). CUSP8 offers an order of magnitude improvement in speed when compared to global high-resolution simulations with cell sizes ranging from 18 to 6 km. EC60to30, while an order of magnitude faster than CUSP8, lacks the improvements in coastal dynamics that motivated the creation of the CUSP8 mesh. Performance tests were run on Grizzly at Los Alamos National Laboratory. Grizzly is an Institutional Figure 12. Resolution Study: Same as Figure 11. A higher-resolution improves the dynamics of the Gulf Stream significantly, with CSUP8 approaching the dynamics of NA8 along the coast. Variability increases with increasing resolution.
Computing (IC) cluster, running on the TOSS operating system (Tri-Lab Operating System Stack) and using the Intel OmniPath interconnect. Each core is a 2.1 GHz Broadwell with a 45 MB cache, with 36 processors per node.

Conclusion
Overall, this mesh resolution case study indicates that simulations are robust to changes in the mesh.
Changes to mesh quality were found to have little impact on the simulation quality and statistics. Problems with the stability of the simulations at large time steps occurred in spinning-up the two most degraded  configurations, but with a modified time step, these simulations were found to perform similarly to the non-degraded cases. Such behavior is consistent with the expected reduction of the time step due to CFL limits associated with heavily degraded meshes that incorporate small grid cells. Despite previous theoretical analysis suggesting a strong link between mesh quality and numerical discretization error in inviscid settings (Peixoto, 2016), it was found that simulation quality was not obviously diminished with increasing levels of mesh degradation. In this sense, it appears that the TRiSK formulation used in MPAS-Ocean may outperform its theoretical bounds in many practical cases, when using the typical suite of parameters for global simulations with non-zero dissipation. Changes to the transition width were also found to have relatively little impact on the quality of the simulations.
It is likely that much of the variation in the transition width study was due to the change in the size of the higher-resolution portion of the transition region rather than the transition itself. The difference between the steady Gulf Stream path in the NA8 simulation and the variable paths in the CUSP simulations shows that the transition region has some impact on variability. It is not clear if this is due to mesh quality or to the effect of changing resolution. In this case, the CUSP meshes had more realistic variability, but it is not clear that this added variability would be desirable in a simulation with realistic forcing. This study also demonstrated that higher coastal resolution improved the dynamics of the Gulf Stream at a much lower cost than a high-resolution global model.
When designing a mesh, the effect of processes outside of the high-resolution region is essential. The transition width study showed that processes within the high-resolution region cannot be properly resolved if they interact with processes in the low-resolution region. For example, in the CUSP8-transition-100 simulation, meanders and eddies crossing into the low-resolution region had a strong impact on the dynamics present along the Gulf Stream within the high-resolution region. More broadly, it is important to evaluate the dependence of the coastal dynamics on basin scale or global dynamics. A coastal high-resolution model may be of limited use if the ultimate drivers of the coastal dynamics are not modeled accurately. For instance, flooding during a hurricane requires that offshore storm surges are modeled at appropriate resolution in order to predict accurate coastal surges.
Physical dynamics considerations appear to be much more important than mesh metrics considerations in these stand-alone ocean simulations. The cases presented here were limited to idealized surface forcing in order to conduct a large parameter study. Future studies will look in more detail at CUSP8 simulations with realistic atmospheric forcing and in coupled configurations, which may have more stringent mesh quality requirements due to cross-component feedbacks. Upcoming papers will examine the CUSP8 mesh in simulations with realistic forcing and active sea ice. Subsequent research has been done to explore the causes of the weak Gulf Stream in the CUSP meshes and has revealed that it may be linked to biases in the Labrador Sea. Work to resolve these issues in both low and variable-resolution meshes is ongoing. Based on the results of this study and further research, modifications have been made to the CUSP8 mesh. The high-resolution region has been extended to encompass the Gulf Stream extension, which will prevent the current being steered by the transition. High resolution will also extend into the Labrador Sea and the coastlines of Greenland and the Canadian Arctic, resolving other processes that are essential to the Gulf Stream, including downwelling in the Labrador sea and the Labrador current.
Similar variable-resolution meshes are in development for investigating other regions of interest, including the Arctic and Southern Oceans. Our results suggest robust capabilities inherent in the MPAS-Ocean discretization and mesh generation approaches. These provide the capability to create a diverse range of variable-resolution configurations, which will allow modelers to accurately resolve additional physical processes at lower computational costs.
Appendix a Figure A1. EC60to30-E3SM-V1 vs EC60to30: Averages are taken from years 2-10 of the simulation, snapshots from 0002-06-01.  In the North Atlantic, the high-resolution grid in (a) is similar to the NA8 grid. The Gulf Stream separation, variability, and transport are much more realistic in the CORE-forced high-resolution run (a) than in any of the climatology-forced runs ( Figure 13).