The effect of vertical coordinates on the accuracy of a shelf sea model

The vertical coordinates (VC) are one of the most important set of configuration options of an ocean model. Optimisation is, however, a non-trivial exercise. We compare nine configurations to investigate different VC options and contrast the Vanishing Quasi-Sigma (VQS), partial step z-level, s-z hybrid and Multi-Envelope (MEs) approaches. Using NEMO model simulations, a hierarchy of experiments are conducted, including: unforced simulations, multi-year climatological simulations with comparisons against tracer profile observations, and tide-only simulations. Hydrostatic pressure gradient errors on the continental slope in the VQS coordinates are found to be consistent with reduced domain-averaged accuracy in both unforced and realistic simulations. Reduced accuracy on the continental shelf is associated with larger advective tracer transports at the shelfbreak. Accuracy is improved by using separate definitions of the computational surfaces on the shelf and slope using the MEs and s-z hybridisation approaches. MEs configurations employing VQS on the continental slope with a computational slope steepness parameter, 𝑟𝑚𝑎𝑥 , of 0.04–0.07, perform comparably with s-z hybrid configurations. Restrictions on the tilt of computational surfaces on the shelf and upper slope appear less important. In contrast, tide-only experiments without stratification show that tidal simulation quality is linked with accurately representing the shelf bathymetry, which favours terrain-following systems. The experiments support transitioning the vertical coordinates across the shelfbreak using either a MEs or hybrid s-z approach as a flexible route to improving accuracy in regional and global models.


Introduction
It can be time-consuming to choose a suitable vertical coordinate system for an ocean model, yet it has been argued that the vertical coordinates are the single most important configuration option because they steer the representation of modelled processes (Griffies, 2004). A significant factor is the effect on the model's representation of processes at regions of steeply sloping bathymetry, such as at the continental slope, which can be important in terms of transports and dissipation (e.g. Zhai et al., 2010;Marsh et al., 2017;Dukhovskoy et al., 2006;Desbruyères et al., 2020). The vertical coordinates are particularly important for shelf sea models where there is a need to accurately represent shallow water dynamics as well as include the adjacent deeper off-shelf regions. In such situations the vertical coordinates play an important role in how on-and off-shelf processes are represented and therefore, in some sense, how the shelf seas are forced via cross-shelf exchange (Huthnance, 1995;Brink, 2016). The three main coordinate options are: (i) vertical coordinates defined by depth, giving horizontally uniform ''z-levels'', (ii) terrain-following sigma-coordinates following the bottom topography, (iii) coordinates that follow isopycnals. Each has their strengths and weaknesses, and it can be possible to hybridise the coordinates in density driven flow that contributes to shelf-ocean exchange (Luneva et al., 2020), although the larger scale significance of intermittent cascading is hard to quantify (Ivanov et al., 2004). The representation of bottom topography at ocean boundaries is also important for the simulation of ocean boundary dynamics, such as boundary currents (Ezer, 2016) and boundary waves (Dukhovskoy et al., 2006;Wise et al., 2020a), which can influence coastal sea level (Wise et al., 2018(Wise et al., , 2020bFukumori et al., 2015;Calafat et al., 2012) and communicate changes in ocean circulation around ocean basins (Johnson and Marshall, 2002;Roussenov et al., 2008). The artificially extreme depth gradients found in z-level models can be reduced, though not removed, by introducing partial steps in the bottom cells (Pacanowski and Gnanadesikan, 1998).
It is possible to define an arbitrary vertical coordinate = ( , , , ), where ( , , ) are the spatial coordinates and is the temporal coordinate, with the restriction that is a monotonic function of height above the sea surface, . The primitive equations can then be transformed into the ( , , , ) coordinate system ( -coordinates) such that a dependent variable can be described as a function of the new coordinates, i.e. = ( , , ( , , , ), ). A particular and common application of this transformation is to allow the spatial variation to fit with the changes in the bottom topography. This gives terrain-following coordinates that are commonly referred to as sigma ( )-coordinates, i.e. = = ∕ ( , ), where is the ocean bottom depth. Terrain-following coordinates naturally improve on the noted deficiencies in z-level models by representing topographic features more faithfully, and they can do so with a coarser overall vertical resolution since computational levels naturally compress where the bathymetry is shallower. A well studied problem that arises in numerical models with -coordinates, however, is that the hydrostatic pressure gradient is transformed to become the sum of two terms that may be large, of comparable magnitude and of opposite sign; in finite difference form this can introduce a truncation error (Haney, 1991;Beckmann and Haidvogel, 1993;Mellor et al., 1994Mellor et al., , 1998. Reducing the hydrostatic pressure gradient (HPG) error is an active area of research and model development and a range of approaches have been investigated to achieve this including: increasing eddy viscosities, preconditioning by removing the background stratification, interpolation of the density field back onto z-levels and various different HPG algorithms and higher order approximations (see Berntsen and Oey, 2010 and references therein for further discussion). This issue makes clear that in addition to the vertical coordinates affecting the representation of physical processes, the accuracy of a numerical model depends on how well the model's numerical schemes cater to the different vertical coordinates. Haney (1991) originally noted that a sufficient condition for maintaining hydrostatic consistency in a finite difference scheme is where is the bottom depth, is the horizontal change in of adjacent cells, = ∕ , is height measured from the undisturbed sea surface (negative in the direction of the ocean bottom), and is the vertical grid spacing. However, since Mellor et al. (1994) showed that Eq. (1) and the discretisation error are not strictly aligned, the concept of hydrostatic consistency has been considered less helpful. A more useful approach was adopted by Beckmann and Haidvogel (1993), who studied the local error due to the HPG error by defining a ''slope parameter'' that is the ratio of the slope of the computational surface and the mean local bottom depth,̄. Defining the slope parameter, , by gives 0 < < 1 with the error tending to zero as → 0. To limit the truncation error a restriction on the maximum allowable slope parameter, , can be imposed such that ≤ . Following Dukhovskoy et al. (2009), we refer to such coordinate systems as Vanishing Quasi-Sigma (VQS)-coordinates, since the computational levels are no longer strictly terrain-following. As described in O'Dea et al. (2012) this can be achieved by creating a smoothed version of the bathymetry with the restriction that ≤ (by deepening relative to the real bathymetry). The terrain-following coordinate system is then generated based on this new smoothed bathymetry surface (envelope) and the grid cells that are deeper than the true bathymetry are masked. For example, see the zonal cross-section of the vertical coordinates in Fig. 1a. A disadvantage of this approach is that vertical resolution can be reduced on the continental shelf and slope. In addition, a saw-tooth representation of the bottom can be created when there is a jump in the number of active levels between adjacent points, e.g. see inset panel in Fig. 1b. This 'in cropping' of sigma layers into the topography potentially provides a physical barrier to any downslope flow. Bruciaferri et al. (2018Bruciaferri et al. ( , 2020 allowed for the definition of a spatially varying by introducing multiple arbitrarily defined surfaces (envelopes), which vertically divide the ocean domain into sets of levels, each potentially possessing an individual 2D varying parameter (envelopes can also be optimised according to other criteria). These coordinate systems have been referred to as Multi-Envelope scoordinates (MEs), and an aspect that will be explored in this paper is their flexibility to allow a stricter (smaller) to be applied on the continental slope and a less strict (larger) on the shelf. In both and (like) coordinates the computational levels are in general misaligned with density surfaces (isopycnals). This misalignment introduces spurious diapycnal mixing via the numerical diffusion that accompanies advection of the tracer fields. Over steep topography this misalignment can be more extreme in terrain-following coordinates (Griffies et al., 2000). Furthermore, a spurious circulation can be generated in terrain-following coordinates as this mixing can reinforce horizontal density gradients; in z-level models the same gradients tend to be weakened. These issues are removed in isopycnic models where the vertical coordinate surfaces follow isopycnals; although this in itself creates complications in representing water masses that are well mixed (Griffies, 2004). At the time of writing isopycnic coordinates are not an option in NEMO.
For our experiments we use NEMO version 4.0.2 and the 7 km horizontal resolution, 51-level, Atlantic Margin Model (AMM7), which is a shelf model spanning the northwest European shelf. A version of this configuration is used as part of the UK Met Office's operational forecasting capabilities for the northwest European shelf (O'Dea et al., 2017) and is developed by both the National Oceanography Centre and the Met Office through the Joint Marine Modelling Programme. The Met Office's operational model uses VQS vertical coordinates with a constant parameter that, as of version CO6 (Coastal Ocean version 6), is relatively large, i.e. weak restriction on the slope of computational surfaces. We will compare the performance of the model using the current CO6 vertical coordinates against 8 other vertical coordinates described in detail in the next section, but encompassing (i) smaller (ii) variable via the MEs-coordinates (iii) different stretching functions, (iv) s-z hybridisation and (v) z-level with partial steps. Our focus is on the coordinates' effect on overall model accuracy, and this implicitly includes the overall effect of the coordinates on the model's numerical approximations in addition to the effect of the coordinates on representing physical processes. As part of the Joint Marine Modelling Programme, this case-study will help inform future development of the vertical coordinate system used in the operational 1.5 km Atlantic Margin Model (AMM15) (Graham et al., 2018;Tonani et al., 2019).
In the following section we describe the vertical coordinates in detail and the core of the model. We then discuss the experiments (unforced, multi-year realistically forced, tide-only) and results, and follow with a summary and conclusion.

Vertical coordinate configurations
Nine configurations of a 51-level vertical coordinate system are used in our experiments. A vertical-zonal plot of each is given in Figs. 1-3 (transect in Fig. 4 denotes location) showing the vertical maximum of (Eq. (2)) for each vertical column. The configurations are described below and summarised in Table 1.

Vanishing quasi-sigma (VQS)-coordinates
Configurations 1-3 (sf_r24, sf_r10, sh_r10) use VQS coordinates that are generated by creating a smoothed version of the bathymetry with a restriction on the Slope Parameter such that ≤ . The coordinates are then generated based on this new smoothed bathymetry (envelope) and then the grid cells that are deeper than the true bathymetry are masked out. This last step ensures that the bathymetry has not been smoothed.
Configuration sf_r24 (Fig. 1a) is the current AMM7 CO6 configuration and uses = 0.24. An parameter of 0.24 is relatively large and allows computational surfaces to closely follow the terrain. The vertical distribution of levels is defined by the Siddorn and Furner (2013) (sf) stretching equation. Configuration sf_r10 (Fig. 1b) is the same as configuration sf_r24 with the exception that = 0.10, hence computational slopes are generally flatter on the slope and shelf; however, there is reduced vertical resolution as well as more saw-toothing (layers that 'in crop' into the bottom topography). Configuration sh_r10 (Fig. 1c) is the same as configuration sf_r10, but with the vertical distribution of levels defined by the Song and Haidvogel (1994) (sh) stretching formulation. The Siddorn and Furner (2013) stretched configuration (sf_r10) typically gives flatter surfaces while the Song and Haidvogel (1994) stretched configuration (sh_r10) exhibits a finer resolution at the bed on the mid-slope where the levels are steepest. Note that parameters can be modified in the Siddorn and Furner (2013) formulation to increase resolution at the bed. Table 1 Description of vertical coordinates used. Stretching options used: SF (Siddorn and Furner, 2013); SH (Song and Haidvogel, 1994); MI (Madec and Imbard, 1996 Fig. 2a) the upper envelope has maximum depth of 250 m and is generated by smoothing the actual bathymetry ≤ 250 m with the Martinho and Batteen (2006) smoothing algorithm with = 0.24. The deeper envelope is a smoothed version of the actual bathymetry with minimum depth of 250 m and using = 0.07. In the case of the deeper envelope, an = 0.04 has also been applied to specific grid points where horizontal pressure gradient errors were larger than 0.1 m/s in a tuning phase. The upper part of the domain is discretised using 25 VQS levels stretched via the Siddorn and Furner (2013) formulation while the deeper part of the domain consist of 26 VQS levels distributed to ensure continuity of the Jacobian of the vertical coordinate transformation (see Bruciaferri et al., 2018 for details). This configuration has the same as configuration sf_r24 on the shelf but a lower on the slope. Configuration MEs_r10-07 ( Fig. 2b) is the same as MEs_r24-07, except that the upper envelope is constrained by an = 0.10 (the same as sf_r10 and sh_r10), which enables us to isolate the effect of on the shelf and at the shelfbreak.

S-z hybrid and z-level coordinates
Configurations 6-8 (sz_r10_s26, sz_r24_s21, sz_r10_s21) use hybrid sz coordinates and configuration 9 (zp) uses z-levels with partial steps. Configuration sz_r10_s26 (Fig. 2c) uses hybrid s-z coordinates where the bottom 25 levels are z-level with partial steps and the upper 26 levels (approximately upper 490 m) are VQS with an = 0.10. The s-z coordinates employ vertical stretching using a double tanh function based on Madec and Imbard (1996) in the z-levels, transitioning to pure sigma (all layers of equal thickness) in shallow waters. Configuration sz_r24_s21 (Fig. 3a) also uses hybrid s-z coordinates, but here the bottom 30 levels are z-level with partial steps and the upper 21 levels (approximately upper 225 m) are VQS with an = 0.24. Configuration sz_r10_s21 (Fig. 3b) is the same as sz_r24_s21, but using an = 0.10. Contrasting the performance of configurations sz_r24_s21 and sz_r10_s21 will again test the sensitivity to on the shelf and at the shelfbreak. Contrasting the performance of configurations sz_r10_s26 and sz_r10_s21 will test the sensitivity to the number of slevels just below the shelfbreak. The comparison of the MEs and s-z hybrid configurations allows us to quantify the effect of gently sloping computational surfaces on the continental slope.
Finally, configuration zp ( Fig. 3c) uses z-levels with partial steps and the same double tanh stretching used in the z-level portions of the s-z hybrid configurations, which concentrates resolution towards the surface.

Discussion and summary
With these nine configurations, summarised in Table 1, we are able to investigate a number of permutations: (i) the VQS coordinates with different allow us to examine the effect of restricting the tilt of computational surfaces at all depths; (ii) the MEs coordinates allow us to examine the effect of individually on the continental shelf and slope; (iii) the s-z coordinates allow us to examine the effect of z-levels with partial steps below the shelfbreak; (iv) the z-level with partial steps allows us to contrast sigma and z-levels on the shelf; (v) the different stretching options allow us to compare the effect of vertical resolution. Fig. 4 shows the vertical resolution of the nine configurations at points on the shelf, at the shelfbreak and in the deep ocean. The sf_r24, sf_r10 and sh_r10 configurations have notably higher resolution on the shelf than the other configurations; however, comparing them against one another makes clear the effect that different stretching formulations and have on vertical resolution. The Siddorn and Furner (2013) stretching (sf_r24, sf_r10) has in general a more uniform cell thickness throughout the water column with smaller thicknesses  Table 1) at 47 • N showing the vertical maximum of (Eq. (2)) for each vertical column. The thick black line indicates the bathymetry and the grid cells without colour shading are masked cells.  Table 1) at 47 • N showing the vertical maximum of (Eq. (2)) for each vertical column. The pink dashed line indicates the transition between envelopes. The thick black line indicates the bathymetry and the grid cells without colour shading are masked cells.

Fig. 3.
Depth-longitude slice of vertical coordinates 7-9 sz_r24_s21, sz_r10_s21, zp (see Table 1) at 47 • N showing the vertical maximum of (Eq. (2)) for each vertical column. The thick black line indicates the bathymetry and the grid cells without colour shading are masked cells. in the upper and lower levels. A more uniform surface layer may give more consistent air-sea fluxes, however the Song and Haidvogel (1994) stretching (sh_r10) has thinner cells around 2500 m depth compared to the Siddorn and Furner (2013) formulation, which we see from Fig. 1c appears to give increased resolution just above the mid region of the continental slope. On the shelf the MEs, s-z and zp configurations have lower resolution. This is more pronounced in the lower half of the water column for the s-z and zp configurations. This is also the case at the shelf break. From Figs. 1-3 it is clear that the trade-off for increased resolution is greater potential for HPG errors, with the sf_r24 and zp configurations at opposite extremes. This is particularly notable from the mid-continental slope region upwards onto the shelfbreak. The MEs and s-z configurations can be thought of as compromises (of varying degrees) between these two extremes; indeed this is the general basis for the presented ordering of the configurations.
Bottom boundary processes, which can have vertical scales smaller than in the interior ocean have been demonstrated to be particularly sensitive to vertical resolution. In experiments investigating the effect of resolution on gravity currents down a slope, it has been shown that increased resolution in the Ekman layer improves representation of dense water cascades (Laanaia et al., 2010;Wobus et al., 2011). Laanaia et al. (2010) showed that when applying a no-slip bottom condition, realistic gravity current dynamics could be captured with 3 computational surfaces in the bottom 40 m, whereas a single surface resulted in a reduced descent rate. Finer resolution in the bottom 40 m-200 m range improved performance, but to a lesser degree. Berntsen et al. (2018) showed that the application of a quadratic drag law friction parameterisation and ''law of the wall'' scaled drag coefficient gave good representation of Ekman drainage and spiral with 4 to 5 grid cells in the bottom Ekman layer. While our focus is necessarily on overall model performance, clearly bottom boundary processes are an important consideration for processes such as Ekman transport off-shelf and the associated downwelling circulation ) that can connect shelf-ocean carbon exchange . Panel (c) in Fig. 4 shows that the sf_r24 configuration has a resolution at the bottom of less than 10 m at 300 m bathymetry with finer resolution at 200 m bathymetry. Fig. 1 panel (a) shows that increased resolution in the Ekman layer is maintained on the upper slope. From there the resolution will typically reduce in the bottom as the number of overall levels becomes a limitation. The MEs configurations also have about 3 layers in the bottom 40 m at 300 m depth (more at shallower depth). The other configurations typically either have 2 or 1 layers. While even the sf_r24 configuration has reduced bottom resolution a great depth, our ensemble does include a range of resolution characteristics including in the bottom boundary region.

Model description
The AMM7 regional configuration of NEMO encompasses the North West European shelf and at the time of writing the Met Office's operational version uses NEMO version 3.6. As part of the transition to version 4 (Madec and NEMO System Team), a review of the vertical coordinates has taken place and we therefore use NEMO version 4.0.2. Before moving on to describe the specific configurations used for the individual experiments, we will first outline the core model set-up, which is based on O' Dea et al. (2017).
The curvilinear horizontal grid extends from 20 • W, 40 • to 13 • E, 65 • and has a meridional resolution of 7.4 km and a zonal resolution  that varies from 9.4 km along the southern boundary to 5.2 km along the northern boundary with a 7.4 km mean. The bottom topography uses bathymetry derived from the North West European Shelf Operational Oceanographic System (NOOS), which is a modified GEBCO 1-arcmin dataset. The bathymetry, Fig. 4a, extends to a maximum depth of order 5000 m to the south west. The model employs a non-linear split-explicit free surface, with variable volume layers, and a 10 s barotropic and a 300 s baroclinic time step. For momentum, lateral viscosity is applied with a bilaplacian operator along computational levels with a coefficient of 10 10 m 4 /s -note that for an idealised two layer reduced gravity system, the Internal Rossby Radius on the NW European shelf would be in the region of 2-5 km (Holt and Proctor, 2008), and as such would not be resolved. The vertical diffusion is via the Generic Length Scale scheme with Neumann conditions at the boundaries and a quadratic friction parameterisation is applied with ''law of the wall" scaled drag coefficient. The Coriolis and momentum advection terms are computed using the vector invariant formulation with the enstrophy and energy conserving (EEN) scheme. Finally, the NEMO specific 'PRJ' option is used for the hydrostatic pressure gradient scheme. This uses a cubic spline to reconstruct the density, with a subsequent analytical integration to calculate the pressure and pressure Jacobian method to obtain the horizontal gradient. It should be emphasised that this is a NEMO specific option and, as such, conclusions drawn in this paper are limited to that scope since it is not clear how the specifics of the numerical scheme transfer to other models.
For tracers, lateral diffusion is applied with a geopotential laplacian operator with a coefficient of 50 m 2 /s, vertical diffusion is via the Generic Length Scale scheme, advection is applied with an FCT scheme that is 4th order in the horizontal and 2nd order in the vertical. The EOS80 equation of state is used. Supporting code and configuration files used in the experiments can be found in the software repository (Wise et al., 2021).

Unforced -hydrostatic pressure gradient experiment
As a first step we compare the different vertical coordinate configurations in an unforced experiment where the model is run from a state of no motion, without boundary forcing (closed lateral boundaries), and with a horizontally uniform density field. With no diffusion on the tracer fields there should be no evolution (the Smagorinsky scheme is used for lateral viscosity). This classic experiment has been used to investigate spurious mixing and currents produced as a result of the truncation error in the hydrostatic pressure gradient calculation discussed in the introduction (e.g. Beckmann and Haidvogel (1993), Mellor et al. (1998), Siddorn and Furner (2013)).
For the initial temperature and salinity we first tried four different climatologies that were representative of the on-shelf and off-shelf stratification for both summer and winter -see Fig. 5. Using each of these four different scenarios in an unforced model run using the sf_r24 configuration, we found that the summertime off-shelf temperature and salinity profiles produced the largest spurious velocities. In order to capture the largest errors, we have used these summertime off-shelf temperature and salinity profiles in the experiments in this section. Fig. 6 shows the domain maximum velocity generated with the different vertical coordinates. It shows significant spurious maximum velocities for the VQS coordinates that use a single value of the parameter (configs. sf_r24, sf_r10, sh_r10). The maximum absolute error is close to 1 m/s for the sf_r24 configuration, which constrains the slope of the computational levels the least. The sf_r10 configuration shows that tightening this restriction reduces the maximum absolute error to ∼ 40 cm/s. Comparing the sf_r10 and sh_r10 configurations shows that the error with the Siddorn and Furner (2013) stretching is approximately 15% larger (5 cm/s) than the error with the Song and Haidvogel (1994) stretching. This suggests that the increased vertical resolution over the mid-region of the continental slope produced by the Song and Haidvogel (1994) stretching might help limit the error. The MEs, hybrid s-z and z-level with partial stepping configurations perform significantly better. The MEs configurations both have maximum errors at the end of the experiment of ∼ 0.07 cm/s, with the change in in the upper 250 m making little difference. The differences between the VQS (sf_r24, sf_r10, sh_r10) and Multi-Envelope (MEs_r24-07, MEs_r10-07) configurations suggests that constraining the tilt of the computational levels on the continental slope improves performance, at least in the idealised case. This is further supported by the hybrid s-z configurations (sz_r10_s26, sz_r24_s21, sz_r10_s21). The spurious velocities for the configurations employing partial step z-levels are a further order of magnitude smaller, except for the sz_r24_s21 configuration -the hybrid case where the terrain following upper 21 computational levels are allowed to slope to a greater extent -where the maximum velocities are comparable to the MEs configurations. The fact that the sz_r24_s21 configuration produces errors an order of magnitude larger than the sz_r10_s21 configuration suggests that at least some of the errors are produced on-shelf or at the shelf break; although clearly it is the sloping of the levels on the continental slope that are having the largest impact. Vertical resolution on the shelf does not appear to affect the error to leading order.
This can be seen more clearly in Figs. 7 and 8 which show maps of the depth averaged horizontal kinetic energy at days 10 and 30 respectively (the white line denotes the 200 m isobath). Firstly, both figures show that there is relatively large spurious kinetic energy in the VQS configurations. While the errors are present both on-and off-shelf, it is evident that the erroneous energy is largest on the continental slope. Secondly, the figures show that the kinetic energy spreads spatially in time. Hence, while the growth in the maximum velocities stops or slows, the velocity errors also disperse, such that the domain wide kinetic energy continues to increase, as shown in Fig. 9a.
The kinetic energy can apparently disperse and grow because the errors in the hydrostatic pressure gradient calculation perturb the tracer fields which then drive velocities that do not decay, at least on the time scale considered here. We note, for example, that running the sf_r24 configuration without tracer advection produces a smaller maximum velocity error that is less noisy. Fig. 9b shows the magnitude of the temperature trend due to advection, volume averaged over the domain. It shows that the temperature field is evolving. Fig. 9c shows that there is an associated geostrophic flow that is larger in the configurations that produce larger temperature transport. Typically this process has been labelled as a Sigma Error of the Second Kind (SESK) defined by a vorticity error (Mellor et al., 1998). Berntsen (2002) concluded that large viscosities, unavailable to realistic ocean models, are required to prevent these errors from growing. They describe a situation where the error in vorticity produces continuous vertical exchange of water masses (light water downward) producing large real pressure gradient errors that may be more important than the initial erroneous pressure gradient error in driving a geostrophic flow. Berntsen (2002) also notes, however, that for any given model setup it is difficult to rule out other errors or limitations present in the model, for example the advection or pressure gradient schemes used, that may have a contributing effect. In any case, the relationship between the velocity and hydrostatic pressure gradient can be considered via the depth integrated momentum equations, where is the depth integrated velocity. For hydrostatic pressure = ∫ 0 ′ , where is density and the gravity constant, the relationship is where is the potential energy, the bottom hydrostatic pressure ( = − ) and is the vertical unit vector, ∇ the horizontal derivative operator, the bottom depth, the negative depth, and 0 the reference density. Hence, the spurious potential energy can be transformed into kinetic energy and dispersed throughout the domain.
In summary, the magnitude of the error increases as computational surfaces are allowed to slope relative to the geopotential surfaces. These errors can be large, particularly on the continental slope, and are associated with large advective transports perturbing the density field. As a consequence, the configurations using z-levels on the continental slope (sz_r10_s26, sz_r24_s21, sz_r10_s21, zp) perform the best, with the MEs configurations (MEs_r24-07, MEs_r10-07) performing better than the VQS configurations (sf_r24, sf_r10, sh_r10). Constraining the tilting of the surfaces on the shelf can also reduce the errors but to a lesser degree, since the bottom topography on the shelf is usually less steep. The sensitivity of the simulations to the tilt of the computational surfaces suggests that there is an opportunity for the implementation of different HPG schemes to reduce this sensitivity (e.g. Shchepetkin and McWilliams, 2003;Berntsen and Oey, 2010;Berntsen, 2011;Berntsen et al., 2015). Sensitivity testing to the configuration of momentum diffusion used did show, however, that the relative performance of the experiments with the different vertical coordinates was robust to different mixing setups.
While the unforced experiments in this section highlight an important issue with the sloping of computational levels, and suggest the mechanism by which it manifests, the results should be taken alongside more realistically forced experiments. One limitation of the unforced experiments is that they assume that density contours are horizontal, which in reality is not the case. This should be advantageous to the z-level coordinates. A second weakness of the unforced experiments is that due to the nonlinear evolution of the tracer fields, which clearly have an important effect on the spurious velocities via the nonlinear equation of state, it is difficult to assess what changes realistic forcing of tracers and momentum will have. It is therefore problematic to assume that the spurious velocities can simply be superimposed over a real climatology of the velocity field. In Section 3.2 we use a realistically forced model to diagnose differences between the vertical coordinates via differences in the modelled temperature and salinity fields to observations, and contextualise with the results from these unforced experiments.

Forced
In this section we investigate whether the anomalies produced in the idealised setups are evident in a model with more realistic tracer fields and boundary forcing. In the following experiments the model, for each vertical coordinate setup, is initialised and forced at the boundaries with reanalysis data and run from 2005 up to and including 2015.

Forcing and initialisation
The model is initialised at the beginning of 2005 and is run with forcing applied at the lateral boundaries, the atmospheric interface and with river run-off along the land interface. Tidal forcing is also applied.
The ERA5 dataset (Hersbach et al., 2018) is used to provide hourly atmospheric surface forcing fluxes for the humidity, precipitation, long and short wave radiation, snowfall, air temperature, wind velocities, and pressure, using the CORE 3.5 Bulk formulation. Monthly diffuse attenuation coefficient data for downwelling irradiance at 490 nm are supplied with constant chlorophyll concentration. The daily river discharge timeseries are produced from an updated version of the river dataset used in Lenhart et al. (2010) combined with climatology of daily discharge data from the Global River Discharge Data Base (Vörösmarty et al., 2000) and from data prepared by the Centre for Ecology and Hydrology as used by Young and Holt (2007). The river runoff forcing frequency is daily and it is distributed throughout the water column. For the tides, 15 constituents are calculated from a tidal model of the northeast Atlantic (Flather, 1981) and applied at open lateral boundaries as the depth mean velocity and sea surface height. For the lateral boundaries the flather radiation scheme is applied to the barotropic velocities and a zero gradient condition is applied to the normal baroclinic velocities. The Global Seasonal Forecast System (GLOSEA) version 5 dataset (MacLachlan et al., 2015) is used for the ocean lateral boundary conditions with the flow relaxation scheme applied to the tracers.

Analysis approach
To assess the performance of the different vertical coordinates, we compare the monthly model output against the EN4 (EN 4.2.1 -D/L:19 Jan 21) observational dataset (Good et al., 2013) of quality controlled subsurface ocean profiles for temperature and salinity over the 10 year period 2006-2015. Note, our aim here is to compare the different set-ups rather than fine tune the model.
The tracer fields are a good comparison benchmark because the EN4 dataset provides a large observational dataset, with a spread both in depth and horizontally across the entire domain, see Fig. 10 for profile locations. The EN4 dataset includes depth profiles of salinity and temperature observations. Each profile is located at a specific longitude and latitude and has observations at different depths such that = 1, 2, … , denotes an individual observation, where subscript = 1, 2, … , denotes the profile, of which there are . The EN4 dataset groups profiles by their year and month. To produce the modelobservation difference anomalies, we first perform a nearest neighbour interpolation of the relevant monthly averaged model fields to each observation profile. We then interpolate the model vertically to the observation depths, and finally take the difference. The set of salinity and temperature difference anomalies for all the profiles are therefore defined by In order to highlight possible depth sampling biases we also define sets of depth averaged difference anomalies It should be acknowledged that there will be spatial and temporal uncertainties in any such analysis; however, by comparing the model tracer fields against over ten million observations, covering a range of depths around the domain, we obtain difference anomaly distributions that present a broad picture of how the vertical coordinates affect long term performance. Given the temporal and spatial non-uniformity of instantaneous observations, this is a pragmatic choice for assessing long term relative performance and a similar approach has been adopted in Graham et al. (2018) Fig. 11 shows the distributions of the difference anomalies for the entire domain over the entire 10 year period as boxplots. The red line denotes the median, the central box extends from the 25th to 75th percentiles and the lines extend from the 5th to the 95th percentiles. There is an overall fresh and cool bias that is consistent across the experiments with all vertical coordinates that can typically be attributed to the forcing datasets used as well as physics options chosen (for example see O'Dea et al., 2017, for discussion of AMM7 biases). From a comparative point of view, however, there is a clear  divergence between the configurations into broadly two categories. The VQS configurations with a single parameter value have larger biases than the others and exhibit greater variability about the median. This is an initial indication that the errors resulting from the hydrostatic pressure gradient errors seen in the unforced experiments (Section 3.1) might be important in the realistic models as well. Note also that the distributions have a skew and long tails, which is a possible indication of the uncertainties in the individual difference anomalies. Since large errors can undesirably influence the mean and standard deviation of the sample, we will also use the median in our diagnostics. Fig. 12 shows the salinity and temperature difference anomalies aggregated over the entire domain for the entire 10 years, giving an overall measure of relative performance. Panels 12(a) and 12(b) plot the mean absolute error and median absolute error against the standard deviation for the salinity defined by Eq. (9). Both metrics show that sf_r24 produces the largest difference and that sf_r10 and sh_r10 both give improvements. The Song and Haidvogel (1994) stretching also appears to perform favourably against the Siddorn and Furner (2013) stretching. These results are consistent with the unforced experiment results. The zp, s-z and MEs configurations all result in a further reduction in errors, with the MEs_r10-07 configuration also reducing variability about the mean error. The similarity between the MEs_r10-07 and MEs_r24-07 configuration and between the sz_r10-21 and sz_r24-21 configurations suggests that it is primarily on the continental slope where the sloping of computational levels becomes an issue. Panels 12c and 12d show the same metrics for the case where the difference anomalies have first been averaged for each profile, as defined in Eqs. (11) and (12) respectively. Here we see a similar result in terms of absolute errors, but there is a notable increase in variability in the configurations using z-levels. The difference between the MEs and s-z/zp configurations in the mean and median metrics shows that the s-z/zp models have very slightly smaller typical absolute errors but larger variability such that they perform slightly worse as measured by the mean absolute error.

Results and discussion
Panels 12(e)-12(h) show the same results for temperature. The picture presented is very similar to that of salinity, with the least constrained sf_r24 configuration showing the most variability and largest errors and the s-z and MEs showing the best performance -with the sf_r10 and sh_r10 somewhere in between. The same metrics applied to the profile averages (panels 12g-12h) give a very consistent picture.
The domain aggregated results from the reanalysis forced experiments qualitatively follow those of the unforced experiments. This suggests that the hydrostatic pressure gradient error at steeply sloping bottom topography (primarily off-shelf on the continental slope), when considered over the entire domain, might dominate the potential benefits of the improved representation of near bottom processes, such as dense water overflows. This point in particular is highlighted by the strong performance of the zp configuration compared to the sf_r24, sf_r10 and sh_r10 configurations. Indeed the difference anomalies for the zp configuration are similar to the more sophisticated s-z and MEs configurations that vary the constraint on computational slope to get the best of both z-like and terrain-following coordinates.
The greatest departure from the unforced experiment results is the performance of the MEs configurations relative to the zp and s-z configurations. In the unforced experiments the spurious kinetic energy for the MEs configuration is orders of magnitude larger than for the zp and s-z configurations, yet in the forced experiments the configurations have comparable absolute errors. There are at least two interpretations of this difference. Firstly, that the MEs configuration is more closely aligned with realistic tracer contours, resulting in smaller errors, and secondly that the hydrostatic pressure gradient errors are small enough, compared with the evolution driven by the boundary forcing, that they become less important relative to the benefits of terrain following coordinates in terms of capturing near bottom processes. In either case, it demonstrates that tuning of terrain-following coordinates using multiple envelopes can produce vastly improved results.
Despite the marginal effect of constraining the computational slopes on the shelf, i.e. above the continental slope, it can be seen that the errors originating on the continental slope coincide with errors on the shelf. This can be shown by looking at the difference anomalies on the outer shelf (region denoted OS in Fig. 13), as a function of depth. Fig. 14 shows the set of difference anomalies defined by (9) and (10), horizontally averaged, for the seasonal climatologies December-February (DJF) and June-August (JJA) for the outer shelf region denoted OS in Fig. 13. For both salinity and temperature, and both seasons, there is a clear clustering of the bias anomalies for the MEs, s-z and zp configurations, which are also relatively consistent with depth. The sf_r24 anomaly, by contrast, has a large fresh and cool bias at the bed, which decreases towards the surface (tending towards the main cluster). This is a strong indication that the errors seen in the unforced experiments remain important not only on the continental slope, but also at the shelf break, influencing processes in the lower half of the water column on the outershelf. The same can be seen to a lesser extent in the biases of both the sf_r10 and sh_r10 configurations. A similar picture emerges in the North Sea and Norwegian trench regions (not shown) (regions NS and NT in Fig. 13). In the English Channel/southern North Sea region denoted EC in Fig. 13, where there is no direct connection to the shelfbreak, the horizontally averaged difference anomalies for the seasonal climatologies (DJF) and (JJA)    Hughes et al. (2006), which are based on long-term observations of water masses in the channel described by Hansen and Østerhus (2000). In (a) the colourbar denotes the bathymetry, the white contour the 200 m isobath, and the coloured arrows represent the circulation pathways for the different water masses, which are labelled in (b). Panel (b) is the vertical slice between the Shetland and Faroe Islands denoted by the black line in (a) between the Shetland and Faroe latitude and longitude coordinates respectively (60.6668, −1.6668), (61.9334, −6.2224). present a different picture - Fig. 15. In this shallow region the VQS coordinate systems perform best, notably for salinity. The main difference may be due to the higher vertical resolution of the sf_r24, sf_r10 and sh_r10 configurations on the shelf. This might be an important detail for an improved representation of the seasonal pycnocline in shallow regions, for example.
The above aggregated results can be put into context by looking at how the different configurations represent water mass properties across the shelf, slope and deep ocean. The Shetland-Faroe channel is a challenging test case because there are a number of water masses present that originate at a range of depths. Fig. 16 (based on Figs. 1.3 and 1.4 in Hughes et al., 2006) is a schematic representation of the water mass circulation and vertical structure in the Shetland-Faroe channel based on long-term observations. The temperature and salinity ranges defining the water masses are given in Table 2 and represent typical properties of the different water masses (Hansen and Østerhus, 2000;Hughes et al., 2006). There are five water masses defined, starting from the deepest: the Norwegian Sea Deep Water (NSDW) is the coolest; the Norwegian Sea Arctic Intermediate Water (NSAIW) is slightly warmer and slightly fresher; the Modified East Icelandic Water (MEIW) is the freshest; the Modified North Atlantic Water (MNAW) is cooler and fresher than the North Atlantic Water (NAW), which is the warmest and saltiest and is part of the upper slope current in the channel. Using the average of the monthly mean temperature and salinity model data over the 10 years (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), Fig. 17 shows the five water masses across the transect shown in Fig. 16a for the 9 model configurations. The temperature and salinity ranges are defined in Table 2. To represent the water masses in the Table 2 Definition of the water masses. Obs denotes that the temperature and salinity ranges are derived from observations and are taken from Hansen and Østerhus (2000). Mod denotes the temperature and salinity ranges applied to model data.  . 17. Water masses for the 9 configurations for the transect shown in Fig. 16a extending across the Faroe-Shetland Channel. The water masses are defined in Table 2. The monthly average temperature and salinity has been averaged over the 10 year period 2006-2015. model, the ranges are in some cases slightly extended, compared to those used in the schematic, to be more continuous. In particular the temperature range to define the MNAW must extend to be significantly cooler to include the waters just above the MEIW. While this is partly due to model bias (generally too fresh and cool), this is also because there will be a degree of mixing between the water masses, which can have very different properties, that may not be reflected in a schematic representation. The black 'unassigned' areas denote mixed regions that are outside the defined ranges of the five water masses. Overall the zp, sz and MEs models do a good job at capturing the structure of the different water masses. A comparison with the VQS coordinates, which struggle to capture the opposite extremes of the NAW and NSDW, shows that the vertical coordinates are making a significant difference. This is consistent with the more extreme fresh and cool biases seen with the VQS configurations in the lower half of the water column in the outer shelf region in Fig. 14.
To demonstrate the relationship between the on-shelf tracer errors and the pressure gradient errors due to the vertical coordinates on the continental slope more explicitly, we consider the rate of change of the salinity and temperature fields due to advection at the shelfbreak. This is important because in Section 3.1 it was shown that the error induced by the steeply sloping vertical coordinates manifests via the nonlinear coupling between the hydrostatic pressure gradient and the transfer of heat and salt, characterised by increased advection of the tracers, i.e. Fig. 9. We consider the advection of the tracers at the shelfbreak, which is defined as the 200 m contour between the points (48N,7.5 W) and (61N,0.5E). This contour is adjacent to the Outer Shelf region, as denoted in Fig. 13. Advection at the shelfbreak is quantified in three steps: (i) take the magnitude of tracer advection at the shelfbreak, (ii) depth average this, (iii) take the median. In formal terms, the salinity and temperature advection in the ( , , ) directions are ∇ ⋅ ( ) and ∇ ⋅ ( ), where ∇ is the derivative operator in the ( , , ) directions and is the velocity vector ( , , ). We use the 10 yr time averages (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) of the tracer advection, denoted by ∇ ⋅ ( ) and ∇ ⋅ ( ). Taking the magnitude and depth averaging gives where is the bottom depth along the contour. To represent the typical values of and we take the median, denoted by ⟨⋅⟩. Fig. 18 shows this typical advective transport of the tracers at the shelf break (⟨ ⟩ and ⟨ ⟩) plotted against the median absolute error in the tracer fields in the Outer Shelf region (region denoted OS in Fig. 13). The black dashed line is the linear best fit. The four panels show an approximately linear relationship between the median absolute error and the magnitude of advective transport at the shelfbreak. It was shown in the unforced experiments, in Section 3.1, that a larger slope parameter, , on the continental slope, was associated with a larger magnitude domain averaged temperature advection. Fig. 18 shows that this is also partially the case in the realistic runs. The magnitude of the  advective transport at the shelfbreak for the s-z and zp configurations is very similar to the MEs configurations, which have an of either 0.04 or 0.07 on the continental slope. This suggests that the HPG error in the MEs configurations is not a leading order effect in the realistically forced experiment. However, we see that an on the continental slope that is larger than 0.07 is associated with an increase in advective transport in the realistically forced experiments. Hence for 0.07 ≤ ≤ 0.24 on the continental slope, the results suggest that the errors in the tracer fields on the outer shelf are connected to larger advective transports at the shelfbreak, which in the unforced experiments were associated with the HPG error due to the larger on the continental slope. This relationship breaks down for < 0.07 on the continental slope, suggesting that this is close to a critical value where the HPG error becomes small enough so as not to be of leading order in a realistically forced model. We note that in preliminary testing we used untuned MEs configurations with = 0.07 everywhere in the lower envelope and that the error in the tracer fields was larger than when using the tuned MEs configurations (with = 0.07 or = 0.04 in the lower envelope). The errors were, however, still smaller than those produced by the = 0.10 configurations, which further supports the interpretation above.
To summarise, the realistically forced runs broadly support the theme of the results of the unforced hydrostatic pressure gradient experiments. The erroneous velocity dominates the benefits of the terrain following coordinates, but importantly only in cases where the slope of the computational levels are allowed to deviate from the horizontal too much over the continental slope, i.e. where the slope parameter is larger than 0.07 on the continental slope. In fact in some locations it should be restricted even further i.e. fine tuning the most sensitive grids points with a smaller , as was done with the MEs configurations. These errors seem to penetrate across the shelfbreak onto the shelf, but the effect of restricting the slope parameter on the shelf itself appears far smaller. The fact that the MEs and s-z have comparatively similar results (while the sf_r10 and sh_r10 performed worse) suggests that = 0.07 on the continental slope is in the proximity of a general critical value in terms of getting the best out of the VQS coordinates for our NEMO shelf-model configuration.

Tide-only
Higher frequency dynamics are also important for shelf sea model applications, such as storm surge and tidal modelling. Calculation of the surge residuals requires tidal simulation and the quality of the tidal simulation can be a good indication of model performance at these frequencies (Graham et al., 2018). We therefore carry out a tide-only experiment to investigate the effect of the vertical coordinates on the tides. For this experiment the model is initialised from rest using a constant temperature and salinity of 10 • C and 35 PSU respectively. The tidal forcing (15 constituents) is the same as in the forced runs in Section 3.2. The simulations are run for two years and a comparison is performed between the harmonic analyses of the model sea surface height (SSH) output and tide gauges. There are 480 tide gauges and they are distributed as in Fig. 19. Here we are interested in the relative performance of the different tidal simulations. Fig. 20 shows the square root of the mean of the squared error (RMSE) against the mean error (BIAS) for the M2 and S2 amplitudes and phases, where the error is defined as the model minus the tide gauge. The results show that performance of the configurations can be separated into four groups that differ by their configurations' setups on the shelf (and around the shelfbreak). The smallest RMSE is achieved by those coordinates which allow greater tilting of the computational surfaces on the shelf ( = 0.24 for approximately upper 250 m -sf_r24, MEs_r24-07, sz_r24_s21). The next smallest RMSE is achieved by those configurations using the hybrid s-z or MEs coordinates with = 0.10 on the shelf (MEs_r10-07, sz_r10_s26, sz_r10_s21). The VQS configurations with = 0.10 (sf_r10, sh_r10) have a similar RMSE for amplitude but larger RMSE for M2 phase. Finally, the zp configuration has the poorest performance in all respects.
In these experiments, the spurious velocities investigated in the previous section do not factor, because there is no stratification. Here, the configurations with the most faithful representation of the real bathymetry on the shelf ( = 0.24) perform the best. The three configurations with = 0.24 on the shelf have different vertical resolutions on the shelf and different representations of the bottom topography on the continental slope, yet they all have a similar RMSE and bias. In shallow regions the bottom topography is an important factor affecting wave propagation and this appears to be the most important determinant of tidal simulation accuracy tested by our experiments.
The difference between the s-z and MEs group using = 0.10 on the shelf, and the sf_r10, sh_r10 group (all these configurations have = 0.10 on the shelf) is that the latter group of configurations have a notable saw-tooth representation of the bottom topography on the shelf, which appears to degrade tidal simulation. Note that the s-z and MEs configurations avoid the saw toothing because more of their computational surfaces vanish beneath the bottom topography on the continental slope where the = 0.07, in this sense the setup below the shelfbreak remains important. It is consistent with the analysis above that the zp configuration should perform worst, because it has the least accurate representation of the bottom topography on the shelf as well as a step like representation.
The results here are important for tidal simulations using a homogeneous ocean, however they diverge from the results of the unforced and realistically forced climatological experiments because there can be no hydrostatic pressure gradient error.

Summary and conclusions
The vertical coordinates are one of the most important configuration options in an ocean model, but they can require a significant investment of resources to optimise. Here we have compared 9 vertical coordinate configurations by carrying out a hierarchy of experiments using a shelfseas model with realistic bathymetry. To create the 9 different vertical coordinate systems we have used 4 different approaches, namely: Vanishing Quasi-Sigma (VQS), Multi-Envelope s-coordinate (MEs), z-level with partial steps and hybrid s-z-level. The overarching aim has been to investigate their effect on model accuracy and in the process identify key configuration options. In particular, we have investigated the importance of allowing the computational surfaces to follow the terrain on the continental slope and shelf, and we have looked at the effect of step-like and saw-tooth representations of the bottom topography. We have used the NEMO ocean model and the Atlantic Margin Model configuration at 7 km horizontal resolution (AMM7) for our shelf sea model and while our results are specific to this configuration, it is reasonable to expect that they are informative more generally.
We found that the scale of the error in the classic unforced hydrostatic pressure gradient experiments (e.g. Haney (1991), Beckmann and Haidvogel (1993), Mellor et al. (1994)) with realistic bathymetry gave a reasonable (up to a point) qualitative indication of relative performance between the nine vertical coordinates when used in a realistically forced climatological simulation and compared against observations of temperature and salinity. This is important because, as noted originally by Haney (1991), while the hydrostatic pressure gradient truncation error in s-coordinates can be investigated relatively cheaply, it is more difficult to assess the effect of the error on a realistic model. This is further complicated when trying to compare against coordinate systems that have different deficiencies, e.g. the step-like representation of bathymetry in z-level coordinates.
The tilt of the computational surfaces relative to depth, as measured by the slope parameter (Eq. (2)), on the continental slope was found to be the dominant determinant of accuracy down to a critical value of ≈ 0.07, where is the maximum allowed value of before computational surfaces vanish beneath the bathymetry. Allowing the surfaces to closely follow the terrain everywhere, as in configuration sf_r24 ( = 0.24), resulted in the generation of spurious velocities of the order of 10s of cm/s in the unforced runs and large transports of temperature and salinity. In the realistic simulations this coincided with the poorest domain aggregated performance, relative to the other configurations. Increasing the constraint on the tilt of computational surfaces on the continental slope improved the simulation of the domain-wide tracer fields. An of 0.10 improved results and an of 0.07 improved results further. By then restricting specific sensitive grid points on the slope further to = 0.04 (identified as grid points where large velocities would otherwise form in the unforced experiments in Section 3.1), we found very little difference between having partial step z-levels on the slope and VQS coordinates via the (MEs)-coordinates. Indeed the MEs-coordinates produced less error variability and a slightly smaller mean absolute error (although the median absolute error was slightly smaller for the s-z hybrid coordinates). We therefore conclude that an of approximately 0.04-0.07 on the continental slope is optimal for VQS coordinates.
In the unforced experiments we found that the magnitude of the spurious flow due to the hydrostatic pressure gradient error was associated with the magnitude of tracer advection. In the realistically forced experiments we found that, across all nine configurations, there is an approximately linear relationship between error of the tracers in the outer shelf region and the typical magnitude of tracer advection at the shelfbreak. In the unforced experiments, increased tracer advection was associated with increased on the continental slope. This appears to be true in the realistically forced experiments for >= 0.07. For < 0.07 it appears that the HPG error becomes small enough so as not to be a leading order effect.
In the unforced and realistically forced experiments we found a decreased on the shelf, from 0.24 to 0.10, made little difference -with the = 0.10 giving perhaps slightly better performance overall. We also found that in the case of the s-z hybrid coordinates (s above z), having the upper s-levels in the top 490 m (26-levels) or top 225 m (21 levels) made only a small difference, with the latter slightly preferable overall. This suggests that the larger errors originate below 490 m, although this may be domain specific.
The tide-only experiments primarily highlighted that in homogeneous models of the tides, accuracy is improved by vertical coordinates that most closely represent the shelf bathymetry. For example the sf_r24, MEs_r24-07 and sz_r24_s10 configurations, which all allow the coordinates to closely follow the terrain on the shelf, performed the best. The VQS coordinates in NEMO can lead to a saw-toothed (incropping of layers) representation of bottom topography on the shelf, which resulted in a poorer tidal simulation. Similarly, the step like representation of the shelf bathymetry in the zp configuration negatively impacted the solution. While the coordinate system used off-shelf made little direct difference to the tidal simulation, it should be noted that in the VQS system, layers that vanish on the continental slope may then reemerge on the shelf resulting in saw-toothing. These unwanted features can be minimised by either relaxing the criteria and/or defining Multiple Envelopes or using hybrid s-z coordinates on the shelf and slope.
Our experiments show that performance gains can be made by defining different VQS setups on the continental shelf and slope (transitioning at approximately 250 m depth with an < 0.07 applied on the continental slope) via the Multi-Envelope approach or by hybridising VQS coordinates above partial step z-levels. The Multi-Envelope approach offers flexibility for greater tuning to a specific domain while the hybrid s-z is perhaps simpler. There appears to be insensitivity to the precise depth location of the transition between coordinates providing it is not deeper than the typical upper slope (> 400 m), which means that this approach is well suited to many different types of model domain e.g. global and regional.
This case study has involved the use of a regional model at a fixed 7 km horizontal resolution and we expect that optimal parameter values will also depend on the model resolution and the geographic region being modelled. Such dependencies are possible avenues of future research. Furthermore, there is value in more specific process oriented studies that investigate the effect of the vertical coordinates on cascading and overflows with a realistic shelf sea model (e.g. dependence on near bed resolution), as well as exchange at the shelfbreak. Such investigations will help to inform the refinement of the vertical coordinates in more targeted studies. Finally, we note that HPG algorithms are an active area of research and model development. The implementation of different 4th order algorithms, for example the 4th order (Shchepetkin and McWilliams, 2003) method, holds the potential to enable a relaxation of the criteria, with possible benefits for the simulation of bottom boundary layers over steep topography.

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