The effect of various vertical discretization schemes and horizontal diffusion parameterization on the performance of a 3-D ocean model : the Black Sea case study

Results of a sensitivity study are presented from various configurations of the NEMO ocean model in the Black Sea. The standard choices of vertical discretization, viz. z levels,s coordinates and enveloped s coordinates, all show their limitations in the areas of complex topography. Two new hybrid vertical coordinate schemes are presented: the “s-on-top-of-z” and its enveloped version. The hybrid grids uses coordinates or enveloped s coordinates in the upper layer, from the sea surface to the depth of the shelf break, and z-coordinates are set below this level. The study is carried out for a number of idealised and real world settings. The hybrid schemes help reduce errors generated by the standard schemes in the areas of steep topography. Results of sensitivity tests with various horizontal diffusion formulations are used to identify the optimum value of Smagorinsky diffusivity coefficient to best represent the mesoscale activity.


Introduction
The early ocean circulation models (Sarkisyan, 1962;Bryan, 1963) all used z coordinate vertical grids and a very coarse representation of the bottom topography.However, meteorologists who started using numerical models a decade earlier (Charney et al., 1950) have found that using a z coordinate system has certain computational disadvantages in the areas of varying topography, in particular in the vicinity of mountains.To resolve this problem a terrain-following sigmacoordinate system was introduced (Philips, 1957).Similar constraints on the use of z coordinate have been found in ocean modelling (see e.g.Blumberg and Mellor, 1983) in particular where both the shelf seas and the deep ocean are included in the model domain.The z level systems have a lower number of layers in the shallow than in the deep regions, which makes it difficult to resolve the water column properties equally well.The effects of vertical discretization on the bottom boundary layer dynamics were studied by Ezer and Mellor (2004) for the idealised case of plain bottom slope.In this paper we extend such analysis for real world bathymetry which combines gentle shelves and steep continental slopes.
An adequate resolution of the near-bottom dynamics is important for the study of dense water overflows through the straits and sills (e.g.Legg et al., 2008;Berntsen et al., 2010) and exchanges between the shelf and the deep sea (e.g.Huthnance, 1995;Grégoire and Beckers, 2004;Huthnance et al., 2009).An accurate modelling of shelf-sea exchanges is particularly important for the Black Sea, which is the world's largest landlocked basin with an extensive shelf occupying about 20 % of the whole sea area (Shapiro, 2008).As a result of its separation from the world oceans by long and narrow straits, and enhanced eutrophication due to river discharges, the Black Sea has seen dramatic degradation of its ecosystem from [1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990].(Besiktepe et al., 1997).The flow of cold near-bottom waters from the shelf into the interior is a significant contributor (Stanev et al., 1998;Demyshev et al., 2002) to the maintenance of the cold intermediate layer (CIL), which separates oxygen-rich surface waters from the anoxic layers below.The shelf-sea exchanges both in the surface (Ginsburg et al., 2002;Zatsepin et al., 2003;Shapiro et al., 2010 ) and in the bottom (Shapiro et al., 2011) layers have significant impact on the state of the shelf ecosystem and support the self-cleaning of shelf waters.For these reasons, it is desirable to reproduce reasonably well the dense water cascades from the shelf to the open sea.An example of a dense water cascade from the shelf observed in May 2004 is shown in Fig. 1.
In order to resolve processes both in the surface and bottom boundary layers, the terrain-following vertical grids are widely used in modern numerical models, such as POM (Princeton Ocean Model) or POLCOMS (Proudman Oceanographic Laboratory community model system).However the terrain-following grids have their own disadvantages due to the errors in calculation of the pressure gradient force (see e.g.Rousseau and Pham, 1971;Mellor et al., 1994;Ezer et al., 2002).The error is partly caused by violating the condition for "hydrostatic consistency" (Rousseau and Pham, 1971), which is often written in the form as follows: where σ is the sigma coordinate of a numeric cell, H its depth below sea surface, δ x H the change in depth of horizontally adjacent grid cells, and δσ the vertical grid size in sigma coordinates.On a sloping topography this condition is severely restrictive.Mellor et al. (1994) have shown that, when using fine vertical resolution, a numerical scheme can be hydrostatically consistent even when R > 1 as the truncation error contributing to the "inconsistency error" is proportional to the following: and can decrease with increasing vertical resolution.There were multiple attempts to introduce alternative or generalised vertical grid systems that would increase the accuracy of numerical calculations; some of them use vertical grids which are not geometrically fixed but evolve in time.Examples include adaptive vertical grids (see Hofmeister et al., 2010, and references therein), isopycnal models such as MICOM (see Bleck et al., 1992), and hybrid models that combine isopycnal, sigma and z coordinate systems such as HYCOM (see Bleck, 2002).Each of these model types has its own advantages and disadvantages (see e.g.Hofmeister et al., 2010).However the most popular grids used in ocean modelling are based on traditional sigma and z level schemes or their combination and variants due to their computational efficiency.Ezer and Mellor (2004) compared sigma and z level grid simulations in a highly idealised dense water flow down a sloping bottom.The z level grid performed relatively poorly due to the step-like representation of topography.However increasing both the horizontal and vertical resolution in the z  (Shapiro, 2008).level model converged results toward lower-resolution sigma grid (Legg et al., 2008).Mellor et al. (1994) demonstrated, in an idealised setting, that the model is stable and the spurious currents decrease with time even when R = 13 if the vertical resolution is sufficiently fine.The modelling of a near-bottom density current over a very steep (39 >) slope using the s coordiante grid (a derivative of sigma) was stable and in a good agreement with laboratory experiments at R = 101 subject to having 4-10 s levels within the thin bottom plume (Wobus et al., 2011).The "safe" value of R also depends on the actual topography and the numerical scheme (Shchepetkin and McWilliams, 2009), so sensitivity studies are required to test vertical discretization for specific model configurations.
In this paper we present results of sensitivity tests for NEMO numerical model (Madec, 2008;NEMO, 2010) with a number of alternative vertical grids.For this study we use the latest version of NEMO-SHELF developed at the UK Met Office (O'Dea et al., 2012).Alongside the three standard grids incorporated in NEMO-SHELF, we developed two nonstandard grids which combine the features of both the z level and terrain-following coordinates.In contrast to earlier studies concerned with idealised topography, we carry out our testing for the complex topography of the Black Sea, where an accurate simulation of the dynamical processes near the shelf break is highly desirable.We also evaluate efficiency and identify optimised parameters for the Smagorinsky algorithm for horizontal diffusion, which has become the preferred method in modern grid-point models (Becker and Burkhard, 2007).

The grids
The NEMO-SHELF numerical model is a 3-D hydrostatic, baroclinic primitive equation model laid out horizontally on the Arakawa C-grid (Madec et al., 1998;NEMO, 2010;O'Dea et al., 2012).Five configurations of NEMO were generated for this study by utilising a number of vertical grids.These grids include 3 standard grids used in NEMO-SHELF: z level (ZCO), s coordiante (SCO) (see Song and Haidvogel, 1994), s coordiante with enveloping topography (Z*SIG, using designation from O'Dea et al., 2012) and two grids specifically developed for this study (Shapiro et al., 2012): the hybrid "s-on-top-of-z" (SZH) and the hybrid "senveloped-on-top-of-z" (SZHENV).Bottom topography was derived from ETOPO2 database (ETOPO2v2, 2006).In order to suppress unnecessary curvature of terrain-following coordinates in the upper layer due to variations of bottom topography below 1550 m, the deepest numerical level (w level) was forced to be horizontal at this depth level to represent an "artificial bottom".It is a common view that the currents below 1000-1200 m are negligibly small (Demyshev et al., 2002;Enriquez et al., 2005); so the waters between 1550 m and the seabed were excluded from the sim-ulation.The transition between the artificial and real bottom was slightly smoothed with a 9-point, 2-D filter.Another advantage of such configuration is that, by reducing the depth of the sea, we can increase the barotropic time step in the time splitting scheme.
The parameters of all NEMO configurations (except the grids) were kept as identical as possible.They all have 33 vertical levels, 1/12 deg zonal and 1/16 deg meridional resolution ( about 6.5 km).The model uses the total variation diminishing (TVD) advection scheme in the horizontal, the piecewise parabolic method (PPM) in the vertical (Liu and Holt, 2010;Colella and Woodward, 1984), the non-linear variable volume (VVL) scheme for the free surface, and the GLS turbulent closure scheme (NEMO, 2010).In the horizontal, the model uses the standard Jacobean formulation for the pressure gradient, 3-D Smagorinsky Laplacian diffusivity for scalar variables (Luneva and Holt, 2010), and a bi-harmonic viscosity formulation with a constant coefficient for momentum.The horizontal viscosity and diffusivity operators are rotated to be aligned with the horizontal levels in order to eliminate contamination of vertical diffusion/viscosity by large values of their horizontal counterparts.The standard grids are shown in Fig. 2. The ZCO grid generates serrated edges near the shelf break, exactly where the cold shelf waters sink to a greater depth in the spring season.The SCO shown in Fig. 2b produces a smooth terrainfollowing grid.However, it violates the "hydrostatic consistency" condition in many places around the shelf break.While such a grid will not necessarily produce wrong results (see Mellor et al., 1994;Wobus et al., 2011), the results should be carefully analysed.
The improvement on the standard SCO grid is the s-grid with the enveloping topography (Z*SIG) (see NEMO, 2010).The generation of the enveloping topography starts from the deepest part of the sea and is calculated using the following criteria: where h env ij , h env i+1, j , and h env i, j +1 are the depths of adjacent grid points in the enveloping topography.Based on common practice the threshold is set to rn cr = 0.1.
The grid cells which are below the actual bottom are masked out in Fig. 2. The enveloping algorithm improves hydrostatic consistency.However, too many numerical layers go below the bottom too early, hence reducing vertical resolution on the shelf and shelf break (see Fig. 2c-d).Such loss of vertical layers is caused by steep lower continental slopes and upper continental rise, where the currents are very small.Due to very weak stratification at these depths, the errors caused by not satisfying condition (2) are unlikely to be large, so that the loss of resolution in the Z*SIG grid is unnecessary.
In order to eliminate or reduce the loss of vertical resolution common for the Z*SIG grids, combine the advantages of both ZCO and SCO grids, and minimise their disadvantages, a new hybrid s-z coordinate grid (SZH) has been introduced.This grid consists of a standard SCO in the upper part of the sea (in this case above 100 m) and ZCO grid below 100 m.This grid is shown in Fig. 3a.The vertical levels are horizontal over the majority of the basin and become curved to follow topography only where necessary.Most places where the depth does not exceed 100 m, the shelf is sloping gently and the "hydrostatic consistency" condition (2) is satisfied even at relatively modest vertical resolution.However there are few places in the Black Sea where shelf is steep even if the depth is less than 100 m.To cater for such cases the SZHENV grid has been developed, which has the enveloping Z*SIG in the upper part of the sea and ZCO in the lower part (see Fig. 3b).There is a minimal loss of numerical levels (3 out of 18) in the shallowest grid cell in the western part of the transect.
It order to compare performance of various grids, we carried out a number of idealised simulations, for which the solution is known, some idealised cases where solution is known only approximately and real world simulations where the results are compared with satellite-derived sea surface temperature.

Idealised simulations
The first set of idealised simulations used the initial conditions produced by spreading horizontally a real temperature/salinity profile, obtained from a climatic atlas of the Black Sea for the month of January (Suvorov et al., 2004).There was no forcing either via meteorological fields, rivers or exchanges via Bosporus.In such conditions there should be neither currents nor a change in temperature and salinity, except a very small, horizontally homogeneous, vertical diffusion due to background diffusivity.However, due to computational errors, the spurious currents develop over time.Temperature stratification also becomes horizontally non-homogeneous.Results of simulation after 3.5 months are shown in Figs. 4 and 5.
Figure 4 shows that in the ZCO configuration the spurious currents reach the value of the order of 1 cm s −1 .The maximum values of spurious currents seen in the cross-section in Fig. 4 are as follows (in cm s −1 ): ZCO = 1.1, SCO = 18, Z*SIG = 6, SZH = 1.5, and SZHENV = 1.2.The spurious currents are within a reasonable range for ZCO, SZH and SZHENV; they are unacceptably high for SCO and Z*SIG.The currents result from spurious density gradients generated by uneven numerical diffusion of temperature and salinity (see Fig. 5).The spurious gradients in ZCO-coordinates seem to come from the bottom boundary condition in the GLS turbulent closure scheme, which generates increased vertical diffusion near bottom due to insufficient vertical resolution.
The second set of idealised simulations was carried out for the initial conditions consisting of a conical interface with a horizontal brim (see Fig. 6f).The temperature above and below the frontal interface was set at 15 • C and 10 • C, respectively; a constant salinity of 18 psu was applied throughout the basin.This setting mimics the domelike structure of the pycnocline in the Black Sea.The frontal interface has a constant slope, and, according to the "thermal wind" equation, the difference between horizontal velocities above and below the interface must not change horizontally.In this setting the z coordinate (ZCO) grid loses its a priori advantage of the numerical levels being parallel to the isopycnals as it was in the first set of simulations.
Figure 6 shows that the best uniformity of the velocity field and minimal spurious currents are generated by the ZCO grid.The Z*SIG and SZHENV produce similar results with very small spurious velocities.The standard s coordiante is the worst; the SZH grid also produces, however to a lesser extent, spurious currents near the steep continental slope.

Black Sea circulation
The Black Sea dynamics for the year 2007 was simulated using the same grids that were used for the idealised runs.The model was initialised as follows.Monthly climatic data for temperature and salinity from the Black Sea atlas (Suvorov, 2004) were interpolated onto the model grids.Then short convective adjustment runs were invoked for a few time steps with very high vertical diffusivity (100 m 2 s −1 ) to eliminate small density inversions due to interpolation.Then the model was spun-up for about 7 months with no meteorological forcing and river discharges following the procedure detailed in Enriquez et al. (2005), i.e. by freezing initial temperature and salinity distributions and allowing the currents to achieve an equilibrium with the density field.
The main simulations were driven by meteorological forcing from the SKIRON forecasting system (http://forecast.uoa.gr/index.php)with 10 × 10 km resolution using the CORE formulation (see Large and Yeager, 2004).The simu-lations included 8 rivers: the Danube, Dnieper, Don, Kizilirmak, Dniester, Bug, Sakarya, and Rioni.Exchange via Bosporus was described as a negative river discharge representing a balance between the upper (outward) and lower (inward) flows.Daily river discharges were interpolated from climatic data sets.River outflows were spread over a few wet grid cells adjacent to river mouths to eliminate the risk of running into negative salinity at high river discharges.All simulations started from 00:00 LT on 1 January 2007.Horizontal diffusion was calculated using a 3-D Smagorinsky formulation with a coefficient C = 0.8 (see NEMO, 2010) for the definition of this coefficient.In order to test the native performance of the model with various grids, no data assimilation or relaxation to climatology was invoked.For the analysis we selected the transect along 31 • E. Whilst we do not have observational data (the "truth") on this transect due to lack of observation, the model skill can be assessed qualitatively by comparing model results with typical partial transects and individual T /S profiles obtainable from the Black Sea databases (Suvorov, 2004).Such an approach allows us to identify significant deviation of the model from typical observations.Temperature cross-sections for 16 April 2007 are shown in Fig. 7.As expected, the z-grid has coarser vertical resolution at the location of the cold water plume penetrating from the shelf into the cold intermediate layer as compared to the terrain-following grids.This results in a "broken" plume and disappearance of bottom cold water on the shelf.Insufficient resolution could also lead to errors in estimating the downslope velocities of the near-bottom plume (Legg et al., 2008;Wobus et al., 2011) and hence the efficiency of the CIL replenishment.However the z-grid conserves well the CIL in the deep part of the sea.The sigma-grids, both SCO and Z*SIG, generate a reasonable representation of the cold water plume.However, it is "broken" on Z*SIG grid, probably due to high and patchy spurious velocities (see Fig. 7).The SCO grid does not conserve well the CIL in the open sea.The hybrid grids, SZH and SZHENV, preserve well both the cold plume at the shelf break and the CIL in the open sea, with SZH doing its job slightly better.

G. Shapiro et al.: The effect of various vertical discretization schemes/horizontal diffusion parameterization
The analysis of zonal velocities for the same transect (not shown) reveals that the standard sigma-grids, SCO and Z*SIG, locate the eastward flowing branch of the Rim Current unrealistically far away from the Turkish coast.
Generation of unrealistic currents by SCO and Z*SIG is consistent with generation of large spurious currents by these grids in the first idealised test.The ZCO, SZH and SZHENV grids produce more realistic results which are generally consistent with those reported in the literature (Stanev et al., 1998;Demyshev et al., 2002;Zatsepin et al., 2003;Shapiro, 2008).
The ability of the model with different grids to reproduce realistic sea surface temperature (SST) is examined in Fig. 8 where model results are compared with SST charts obtained from the Remote Sensing Department of the Marine Hydrophysical Institute (dvs.net.ua).As in the previous test, the model was initialised with climatic initial conditions and ran from 1 January to 7 July 2007 without data assimilation.The main feature seen in the satellite image is a cold water filament spreading from the tip of the Crimean peninsula to the SW corner of the sea and an intensive anticyclonic eddy containing cooler waters of about 23 • C (see Fig. 8f).These features are best represented by ZCO grid, and reasonably well by SZHENV.The standard sigma grids, SCO and Z*SIG, misplace the eddy and produce significant warmer areas along the western coast which are not seen on the satellite SST image.

Effect of horizontal diffusion
The values and algorithms of horizontal diffusion have significant influence on the accuracy of simulation (e.g.Wright and Stoker, 1992).The effect of variations in horizontal diffusion is evaluated here using the hybrid SZHENV vertical grid and the same forcing and initialisation as in the previous section.It order to isolate effects of horizontal diffusion, all other parameters were kept identical.For these tests we used the 3-D Smagorinsky algorithm (Luneva and Holt , 2010) as in the previous section.Both horizontal viscosity and diffusivity are coded in NEMO according to the Smagorinsky formulation, which is discussed in detail in Griffies and Hallberg (2000) for the case of horizontal viscosity.The Smagorinsky algorithm is based on the physical assumption that horizontal mixing, induced by eddies, should be larger near sharp fronts and in coarser-resolution grids, where more subgrid-scale eddies are unresolved.In this algorithm, horizontal viscosity A M and diffusivity A H = A M Pr, where Pr is the horizontal Prandtl number, depend on grid size and velocity gradients.Table 1 shows the value of Smagorinsky coefficient C as defined in NEMO (2010) and its equivalent horcon/Pr = (C/π ) 2 which used in a popular Princeton Ocean Model.All tests utilised a constant bi-harmonic horizontal viscosity with a constant coefficient of (−1 × 10 9 ) m 2 s −1 .For the inter-comparisons of various NEMO configurations with observations, we use the GODAE class 1 metrics (GO-DAE, 2008), i.e. daily averaged charts of sea surface temperature.The lowest value of the Smagorinsky coefficient (C = 0.6)(see Fig. 9a) gives the best resemblance of the re- motely sensed SST (Fig. 8f).The accuracy of model prediction for this configuration is assessed by calculating differences between 8-day SST averages for the period 1-8 July 2007 produced by the model and MODIS-AQUA gridded remotely sensed data (JPL, 2012).The metrics for the model skill are as follows: bias = −0.45• C (model is cooler), RMS = 0.46 • C, and the correlation coefficient R = 0.52.This compares well with the uncertainty in SST obtained from various satellites.For example, the differences between SST from MODIS-AQUA and AVHRR Pathfinder (JPL, 2012) for the same period are similar: bias = −0.43• C (Pathfinder is warmer), RMS = 0.25 • C, and the correlation coefficient R = 0.85.A slightly lower correlation with model data is caused by the displacement of the cold filament compared to observation, which is probably due to the fact that model was run from climate for 6 months without data assimilation.
From Fig. 9b-d, we see that the larger values of the Smagorinsky coefficient (C = 1.4,1.7, 2.2) destroy the strong anticyclonic eddy in the northwestern area of the sea.This eddy is a well-known recurrent feature commonly called the "Sevastopol eddy" (Zatsepin et al., 2003) so that the loss of this feature is not acceptable for an eddy resolving model, and hence lower values for the coefficient C are needed.It should be noted that the SST charts obtained with C = 0.8 (Fig. 8e) and C = 0.6 (Fig. 9a) are nearly identical, so that further reduction below C = 0.8 is inefficient.At even lower values of C, the numerical diffusion inherent to the TVD scheme exceeds the physical diffusion and masks out any further reduction in C. Future implementation of less diffusive horizontal advection schemes in NEMO (for example, the PPM) may help resolve this issue.
The effect of various values of C can be seen from the following example: the difference between simulations with C = 0.8 and 1.7 gives a bias of −0.04 • C, correlation R = 0.96, and an RMS = 0.19 • C. While the impact on the RMS value seems to be relatively small, larger values of C result in destroying the Sevastopol eddy in the NW Black Sea.

Discussion
Currently there are a large number of vertical grids which can be used in ocean circulation models (Ezer et al., 2002;Blumberg and Mellor, 1983;Grégoire and Beckers, 2004).However selection of an optimal grid for a specific application is not straightforward and needs specific experiments and sensitivity tests (Ezer and Mellor, 2004).Here we compare results of both idealised and real world simulations in the Black Sea using different vertical grids, with particular emphasis on the accurate reproduction of processes at the sea surface and at the shelf break.None of the grids was the winner in all tests.The spatial patterns of ocean parameters obtained with different grids differ quite significantly, so they are analysed by considering consistency of charts and crosssections alongside the standard statistical metrics, such as biases and RMS errors.
In the idealised simulations with horizontally homogenous temperature and salinity (but real topography), the main focus was on checking the spurious currents, which in all cases were generated near the shelf break (see Fig. 4).In this test the winners were the z level grid (ZCO) and the combined s-z grid either with (SZHENV) or without (SZH) enveloping topography, all of which generated currents not exceeding 1-1.5 cm s −1 .This inaccuracy seems acceptable in the Black Sea where the Rim Current flowing offshore of the shelf edge has velocities of 20-40 cm s −1 (Zatsepin et al., 2003).The s coordiante (SCO) and s coordiante with enveloping (Z*SIG) failed this test producing spurious currents as high as 6-20 cm s −1 .The SZHENV and SZH were the best in preserving the original stratification, while ZCO, the standard SCO and Z*SIG generated horizontal inhomogeneities within the cold intermediate layer.
In the second idealised test, a conical temperature front centred in the deep sea was used as the initial condition.As the horizontal pressure gradient near the front was constant everywhere, the difference between the current velocity below and above the front should be constant.Near the continental slope the frontal interface became horizontal to exclude the generation of complex currents due to the JEBAR effect.In this test the best performance was achieved by the ZCO discretization, with SZHENV being a close second (see Fig. 6).
The real world simulation was carried out for the Black Sea using high-resolution meteorological forcing, but no data assimilation was performed to separate the natural skill of the model from the quality of assimilation.In terms of reproduction of the SST, the best results were obtained using the ZCO grid, with the SZHENV being a very close second (see Fig. 8).However the cold water plume, which helps replenish the cold intermediate layer in the Black Sea, was better represented by SZH and SZHENV, while the ZCO generated an intermittent, "broken" plume highly affected by step-like structure of bottom relief.

G. Shapiro et al.: The effect of various vertical discretization schemes/horizontal diffusion parameterization
Summarizing the results of grids evaluation, we can see that the z level grid performs well in simulating surface phenomena in the Black Sea, while having difficulties in the near bottom layer on the extensive north-west shelf.This conclusion is consistent with results of the DOME project (Ezer and Mellor, 2004;Legg et al., 2008) where different model grids were tested for highly idealised dense water overflows.By modelling overflows on a plane slope, Ezer and Mellor (2004) concluded that the most apparent examples of cases where model results are affected by the model grid are simulations when the major current is a density-driven, down-slope flow, since the grid determines how topography is represented in the model.Our results confirm that a similar effect takes place in the case of real topography even when the down-slope currents are slow as in dense water cascades from the north-west shelf, or do not exist at all, as in the case of idealised conical interface.Ezer et al. (2002), suggested adding BBL schemes to z level models to allow sufficient penetration of cold shelf water into the deep ocean, a suggestion which led to the development of "z-on-top-of-sigma" numerical grids (see e.g.NEMO, 2010).In this paper an opposite construction is suggested -a hybrid SZH or SZHENV grid, which has a "sigma-on-top-of-z" design.This grid allows reducing the pressure gradient error, and it performed nearly identical to the ZCO in the surface layers in our experiments.It also provided a higher resolution on the shelf and in the near bottom layer, which is required for adequate estimates of the fluxes between the shelf and deep sea waters.The use of the SZHENV grid is costefficient as it uses the same number of grid nodes and does not require additional memory or computing power.
Numerical experiments with varying diffusivity have shown than the best match with the observed sea surface temperature is achieved by using Smagorinsky coefficient C = 0.6-0.8.Further reduction of this coefficient does not result in any visible changes due to the effect of numerical diffusivity inherent to the TVD advection scheme.

Conclusions
This sensitivity study is focused on assessment and optimisation of the performance of NEMO numerical model in the Black Sea using alternative vertical discretizations and horizontal diffusion parameterizations.Analysis of both idealised and real-world simulations shows the following: -The z level grid provides good quality results in the upper layer.As expected, it has difficulty in representation of near-bottom dynamics, which is important for quantifying shelf-deep sea exchanges and their effect on the marine ecosystem in the Black Sea.
-The s coordiante grid and s coordiante with enveloped topography grids better simulate the temperature struc-ture in the bottom layers but produce unacceptable spurious currents at the shelf break.
-The "s-on-top-of-z" hybrid scheme with enveloping topography (SZHENV) proposed in this paper combines the advantages of both z and s level schemes and minimises their drawbacks.
-In order to realistically represent mesoscale eddies, an optimised value of the Smagorinsky coefficient for horizontal diffusion should be applied.
Future developments in numerical algorithms can resolve some of the identified issues.However, implementation of a new code requires significant effort in testing, calibration and validation.This paper shows how improvements can be achieved using existing code by optimising parameters available for the tuning by a user with only modest modification of the "peripheral" model code and without interfering with the main solver.

Fig. 2 .
Fig. 2. Standard NEMO vertical grids shown at a Black Sea transect along 44 • N. Green areas represent wet points included in computations; white areas represent masked-out dry points.Dots show location of T-points within the cells.There is an extensive shelf on the western side of the transect and a very narrow shelf on the eastern part, which is not resolved by the model.(a) z coordinate grid (ZCO); (b) s coordiante grid (SCO).The undulation of the numerical levels is due to the underwater ridge which is shallower than 1550 m; (c) s coordiante grid with enveloping topography (Z*SIG).For clarity, only the upper parts of the grids are shown in subfigures (a)-(c).The Z*SIG grid for the full depth of sea is shown in subfigure (d).