Exploring Viscosity Space in an Eddy‐Permitting Global Ocean Model: Is Viscosity a Useful Control for Numerical Mixing?

A generic shortcoming of constant‐depth (or “z‐coordinate”) ocean models such as MOM5 and Nucleus for European Models of the Ocean (NEMO) is a tendency for the advection scheme to produce unphysical numerical diapycnal mixing, which may exceed the explicitly parameterized mixing based on observed physical processes. Megann (2018, https://doi.org/10.1016/j.ocemod.2017.11.001) estimated the effective diapycnal diffusivity in the Global Ocean Version 5.0 (GO5.0) 0.25° global implementation of the NEMO model and showed that this was up to 10 times the explicit diffusivity used in the model's mixing scheme and argued that this was at least partly caused by large transient vertical velocities on length scales comparable to the horizontal grid scale. The current UK global NEMO configuration GO6, as used in the Global Coupled Model version 3.1 (GC3.1) and UK Earth System Model (UKESM1), is integrated in forced mode at 0.25° resolution with a range of viscosity parameterizations. In the present study, the effective diffusivity is evaluated for each integration and compared with the explicit value from the model mixing scheme, as well as with that in the control (using the default viscosity). It is shown that there is a strong correspondence between lower viscosity and enhanced numerical mixing and that larger viscosities lead to a marked reduction in the unrealistic internal temperature drift seen in the control configuration, without incurring excessive damping of the large‐scale circulation, mixed layer depths, or sea ice cover. The results presented here will inform the choices made in global ocean configurations used in climate and Earth System models following the sixth Coupled Model Intercomparison Project (CMIP6).

One might naively expect that reducing the viscosity in a given ocean configuration should lead to unambiguous improvements in its realism. This may hold for non-eddy-permitting resolutions: for example, Jochum et al. (2008) investigated the sensitivity of the ocean component of the CCSM3 coupled model to the choice of parameters in the laplacian viscosity scheme, which yields a viscosity that depends on both grid size and the flow deformation. At the relatively coarse resolution of this model (1.125° east-west, and from 0.25° north-south at the Equator to around 0.7° at high latitudes), this study found that reducing viscosity generally resulted in improvements in both large-scale circulation and coastal currents. Where there is significant energy at the mesoscale, however, the situation is not so clear. G. Madec (personal communication) suggests that the default value of −1.5 × 10 11 m 4 s −1 of the biharmonic isopycnal viscosity A Lm in 1/4° NEMO simulations (e.g., GO5.0, Megann et al., 2014, andGO6, Storkey et al., 2018), which was originally selected to give a surface eddy kinetic energy (EKE) field of similar level to that estimated from satellite observations, is in fact too low and that this tends to lead to variability in the surface elevation in models that is more noiselike than eddy-like in regions where the baroclinic Rossby radius is unresolved or poorly resolved. This flies against the generally accepted practice in ocean modeling to set the viscosity to as low a value as permitted by numerical stability. Megann (2018) showed that the standard deviation of the vertical velocity in GO5.0 was 2 orders of magnitude larger than the mean, characterizing four distinct regimes for variability in the vertical velocity: tropical instability waves (TIWs), found within a few degrees of the Equator in the Atlantic and Pacific, with wavelength of around 1,100 km and period of around 30 days; near-inertial wavelike features generated close to the western boundaries, which may be associated with near-inertial gravity waves (NIGWs, certainly present in both GO5.0 and GO6 at ¼° resolution, but not resolved in the time-averaged fields available for the present analysis); resolved mesoscale eddies at midlatitudes with length scales of 50-100 km in the vicinity of the separated western boundary currents and the Antarctic Circumpolar Current (ACC); and finally poorly resolved or noise-like patterns at close to the grid scale, seen at higher latitudes and close to coastlines, particularly in the North Atlantic subpolar gyre. The last of these is certainly unphysical, and the aforementioned study ascribes much of the high numerical mixing to them, so a logical approach to reducing the spurious mixing would be to suppress the small-scale motions.
Higher-order viscosity, or "hyperviscosity," schemes such as the biharmonic viscosity used in the standard 0.25° NEMO configurations, or the Smagorinsky deformation-dependent viscosity (Smagorinsky, 1963), act preferentially to dissipate grid-scale horizontal shear, unlike the laplacian (or harmonic) viscosity, which is generally considered (e.g.,  to lead to large-scale flows in ocean models of eddy-permitting or eddy-resolving resolution that are too viscous. The choice of the momentum closure scheme in an ocean model clearly has the potential to change the pathways and the strength of the flow. For example, Chassignet and Garraffo (2001) showed that, even at the relatively high spatial resolution of 1/12°, the separation of the Gulf Stream and the strength of the subpolar gyre in an isopycnic model of the North Atlantic both showed a significant sensitivity to the viscosity parameters. Over longer time scales, changes in viscosity will therefore lead to changes in the tracer fields through advection. At the same time, however, we might expect that the effect of a higher-order viscosity in suppressing horizontal shear at the grid scale will also act as a control for the transient vertical motions at this length scale, and hence on the rate of numerical mixing. Although the viscosity does not act directly on the vertical velocity, grid-scale noise in the horizontal flow will be accompanied by comparable noise in the divergence, and hence in the vertical velocity.
Over time scales of a decade and longer, unrealistically high diapycnal mixing in the ocean interior will inevitably lead to drifts in water properties relative to the initial state (normally based on climatology). It should be noted that model biases from other sources, such as those due to the choice of surface forcing, may dominate over those from numerical mixing, while deeper waters may also be affected by slow changes in the overturning circulation. Such a bias (i.e., cold near the surface and warm at depth) was seen in the comparison study of two coupled models by Megann et al. (2010): the model with a depth-coordinate ocean (HadCM3) had an excessively diffuse thermocline, while another with an isopycnal-coordinate ocean (HYCOM) coupled to the same atmosphere had a much sharper thermocline, along with temperature drifts of the opposite sign to those in HadCM3; indeed, the surface warming and interior cooling in HYCOM suggested that this ocean model might in fact be underdiffusive. Dunne et al. (2012) carried out a similar comparison, replacing the depth-coordinate MOM4p1 with the GOLD isopycnic ocean model in a coupled model and noted a significantly higher (and more realistic) thermocline stratification in the isopycnal-coordinate ocean. Other climate models using other depth-coordinate ocean components, including HadG-EM1 (T. C. Johns et al., 2006) and HadGEM3 (Hewitt et al., 2011), have also shown a temperature bias that reverses sign with depth, similarly consistent with excessive mixing, although in the latter model this was offset by a strong warm surface bias in the Southern Ocean.
This paper addresses the question of whether changes in the viscosity parameterization can lead to an overall improvement in the performance of an eddy-permitting ocean configuration of a z-coordinate ocean model, with particular focus on the numerical mixing in the model. In Section 2, the model is briefly described, and details of the integrations are given. In Section 3, we summarize the analysis used in Megann (2018) and describe the numerical method used to implement it. In Section 4, the effects of increasing the viscosity on the grid-scale noise in the vertical velocity and on the diagnosed diffusivity are presented. We describe the associated reduction in the global deep temperature drifts, the sensitivity of the large-scale circulation indices, Equatorial Pacific metrics, mixed layer depth and sea ice cover to the viscosity in Section 5, and quantify the robustness of the correlations of model metrics with the viscosity. Finally, Section 6 is a summary of results and a discussion.

Model Description
The results presented here are based on integrations of the GO6 configuration (Storkey et al., 2018) of the NEMO ocean model, which is derived from the GO5.0 configuration described by Megann et al. (2014). GO6 is based on version 3.6 of NEMO (Madec & The NEMO Team, 2017). The sea ice is version 5.1.2 of the Los Alamos National Laboratory sea ice model CICE (Hunke & Lipscomb, 2010) in the Global Sea-Ice version 8 configuration (GSI8; Ridley et al., 2018). This configuration forms the ocean and sea ice component of the UK's contributions to the sixth climate model intercomparison project (CMIP6): namely, the Global Coupled Model version 3.1 (GC3.1, Williams et al., 2018) and the UK Earth System Model (UKESM1, Sellar et al., 2019). The GO6 model grid is an extension of the ORCA025 global 1/4° resolution grid (Barnier et al., 2006) that was used in GO5.0 and is now referred to as eORCA025. The model bathymetry is DRAKKAR v3.3, derived from the ETOPO1 data set (Amante & Eakins, 2009) with additional data in coastal regions from GEBCO (IOC et al., 2003).
GO6 uses a total variance dissipation (TVD) scheme (Zalesak, 1979), also referred to as Flux Corrected Transport (FCT), for both horizontal and vertical advection of tracers. The vertical mixing scheme for tracers and momentum is a modified version of the Turbulent Kinetic Energy (TKE) scheme (Gaspar et al., 1990;Madec & The NEMO Team, 2017), with a background vertical eddy diffusivity of 1.2 × 10 −5 m 2 s −1 , decreasing linearly from ±15° latitude to a value of 1.2 × 10 −6 m 2 s −1 at ±5° latitude (Gregg et al., 2003). The mixing effect of tides is parameterized using the Simmons et al. (2004) scheme.
Lateral diffusion of tracers is along isoneutral surfaces using laplacian mixing, with a coefficient of 150 m 2 s −1 ; the mesoscale eddy parameterization scheme of Gent and McWilliams (1990) for baroclinic instability is not used. Free-slip lateral boundary conditions are applied. A biharmonic viscosity is used: the control integration and two sensitivity experiments use a fixed viscosity which has a reference value A Lm at the Equator, and with values reduced poleward as the cube of the maximum grid cell dimension in order to avoid instability in the momentum diffusion equation (see Griffies, 2004), so at 60°N the horizontal viscosity is approximately 1/8 of its value at the Equator. We shall refer to this scheme as the "fixed" viscosity. In an additional three experiments, the biharmonic viscosity is deactivated and replaced with the deformation-dependent viscosity Smagorinsky, 1963), again biharmonic, for which the viscosity B mSmag has the form: where L is a local grid scale, given by and the parameter C Msmag is used to control the size of the viscosity. In a further experiment, the Smagorinsky scheme is applied together with a fixed background value.

Experimental Design
We shall describe seven integrations in this paper, which are summarized in Table 1. These are a control, identical to the GO6 configuration described by Storkey et al. (2018); two experiments with the biharmonic viscosity A Lm increased by a factor of 2 and 3, respectively, with respect to the reference; three with the fixed biharmonic viscosity replaced with the biharmonic Smagorinsky formulation with coefficients C Msmag = 2.0, 3.0, and 4.0; and one with both Smagorinsky and fixed background viscosity applied. In the latter experiment, mixed, A Lm was identical to that in the control but was augmented by Smagorinsky viscosity with C Msmag = 3. We note that the latter experiment perforce used a modified version of the NEMO v3.6 code: in the standard code release, selection of the Smagorinsky viscosity introduces a constant background viscosity, and this was changed so that the fixed term scales as the cube of the grid size, for consistency with the default.  Depth-and Area-Weighted Averages Over Years 1996-2005 Each of the configurations was integrated for a total of 30 years, with the initial conditions and forcing identical to those described by Storkey et al. (2018). The initial conditions for temperature and salinity were an average over years 2004-2008 of the EN3 monthly objective analysis (Ingleby & Huddleston, 2007) and the model was started from a state of rest, with monthly mean output saved for the whole integration and 5-day output for the final decade. The surface boundary conditions were derived from the CORE2 surface forcing data set (Large & Yeager, 2009). Surface salinity is restored to climatology at a rate of −33.33 mm day −1 psu −1 . The two integrations with increased fixed viscosity alone were found to be unstable if started from rest with the increased viscosity, so this was ramped up gradually from the default value of A Lm over the first decade of the respective run. The standard time step of 1,350 s for the 1/4° GO6 configuration was also found to lead to instability with the higher fixed viscosity values, so that this was reduced for the entire production integration as given in Table 1 for these two experiments (bilap_2 and bilap_3). To check that the numerical mixing is not strongly affected by the choice of time step itself, the model was run for 10 years with the standard viscosity, but with half the default time step, and it was confirmed that the temperature evolution in this case was much closer to that in the control experiment than to that in any of the other experiments. Figure 1 shows the absolute value of the biharmonic viscosity in the surface level in the seven experiments (in the cases with the Smagorinsky scheme, the viscosity was averaged over the final 10 years of integration), as well as profiles of the global mean viscosity (Figure 1a). The effect of the Smagorinsky contribution can be seen to be strongly surface enhanced, as expected from the more energetic circulation in the upper ocean; the viscosity in the Smagorinsky-only runs is reduced even more at depths below 1,000 m relative MEGANN AND STORKEY 10.1029/2020MS002263 5 of 29 to the cases with fixed viscosity. The apparent dependence on depth of the fixed viscosity is an artifact, resulting from the combination of the latitude dependence of the viscosity and the spatial variation of the water depth. In the experiments using the Smagorinsky viscosity, it is clear that the viscosity is lower in the centers of the gyres than in the experiments with the fixed viscosity, but larger along the western boundaries and in the vicinity of the separated boundary currents. The effect of increasing the Smagorinsky coefficient is obvious, with the viscosity in smag_2 being lower than that in the control almost everywhere, and in smag_4 producing values on the western boundary a factor of 8 or more larger than that in the control. The viscosity in the mixed experiment can be seen to combine enhanced values in the boundary currents with a higher fixed background viscosity. For simplicity, we shall use the global mean viscosity at selected depth levels as a metric in the following analysis; in the experiments using the Smagorinsky viscosity, this of course conceals the variations with latitude, in particular the low values in the interior of the gyres and the high values in the boundary currents and close to the Equator. Nevertheless, we note that in Smag_2 the viscosity is almost everywhere smaller than in the control; in the experiments with increased fixed viscosity, the global mean viscosity increases exactly as the viscosity in individual cells; in mixed, the viscosity is everywhere larger than in the control; and in the experiments with pure Smagorinsky viscosity, the global mean viscosity increases commensurately with the local values. The global mean is therefore representative of local changes, with a caveat needing to be applied in comparison between pure Smagorinsky cases and the rest of the experiments.

Summary of Analysis
The analysis is identical to that described in Megann (2018). For clarity, we shall define again in this section the main quantities that we shall enumerate in Section 4.
As in Lee et al. (2002), a zonally integrated water mass transformation streamfunction G(Θ, ρ) is defined as where Ψ(Θ, ρ) is the overturning streamfunction at latitude Θ and potential density ρ, and V(Θ, ρ) is the volume below the isopycnal surface ρ and south of Θ. We consider here those processes that change the density below the directly surface-forced surface layer, excluding from discussion those regions in density space with potential density lower than the maximum monthly surface density. If we assume that the density transformation is entirely due to diffusive processes, we can define an effective diffusivity κ eff as where y is the northward spatial dimension, and the zonal integration in the denominator is performed on density surfaces. In reality, as discussed in Lee et al. (2002) and Megann (2018), there are additional contributions to the transformation rate from the nonlinearity of the equation of state, namely cabbeling and thermobaricity, but these tend to be spatially quite localized (mainly in the Southern Ocean and the subpolar North Atlantic: see Klocker & McDougall, 2010). The analysis we present here neglects these contributions, which tend to create negative values of the effective diffusivity in the latter regions. A further approximation is the fact that zonal integrals are present in both the numerator and the denominator of Equation 4, and there is an assumption that the meridional volume transport and the density stratification have little variation, or both vary in a similar way, across each zonal section. This approximation is necessary for the derivation of a practical effective diffusivity, but the associated uncertainties make direct comparisons between κ eff and the explicit diffusivity κ exp less reliable. Nevertheless, the majority of the comparisons we shall make in this paper will be between κ eff in different experiments, and variations of the stratification across the ensemble are very small.
To express the potential effect of transient vertical motions on diapycnal mixing, Megann (2018) defined an effective diapycnal advective velocity g adv as where w is the vertical velocity, the s operator denotes a standard deviation, and the overbar denotes zonal and temporal means on density surfaces. This can be seen to be equal to the standard deviation of the vertical velocity, weighted by the instantaneous local stratification, and may be understood as the potential of the transient vertical motions to advect the density, and hence to perform diapycnal mixing.
In a case where temperature drifts are dominated by diffusive, downgradient vertical heat flux, we can define (similarly to that used by LeClair and Madec) an effective downward heat flux Q eff , averaged globally and over the final year of each integration, as We evaluate Q eff , here from a time series of the annually and horizontally averaged mean temperature.

Numerical Solution
A full description of the numerical solution is again in Megann (2018). In summary, the above expressions were evaluated in the global domain in latitude bands 2° wide. The density coordinate was chosen to be potential density, referenced to a pressure of 2,000 dbar (σ 2 ), and 72 density classes were used (values are listed in the above paper), spanning the range 30.0 < σ 2 < 37.2 kg m −3 : a linear density scale was used for σ 2 < 35.0 kg m −3 , and a logarithmic mapping was used at densities higher than 35.0 kg m −3 . The transformation streamfunction, and thence the effective diffusivity, were averaged from 5-day mean outputs from the last 10 years of each of the 30-year integrations: this was the shortest sampling period available. We will present additional results for the RMS vertical velocity on model levels: for GO6, in contrast to GO5.0, the mean square of the vertical velocity was calculated during integration and saved in the monthly output files and so allows the true RMS velocity to be calculated, which takes account of transients unresolved in the 5-day means discussed in Megann (2018).

Eddy Kinetic Energy and Vertical Velocity
The EKE is the field that is most directly affected by changes in viscosity, since the biharmonic viscosity in the 1/4° GO6 configuration acts preferentially on the horizontal velocity shear at small length scales. In regimes where eddies are poorly resolved (or not resolved at all), the fact that the horizontal velocity is dominated by noise at the grid scale, rather than smooth eddies, means that the horizontal velocity is likely to be significantly more divergent than in the real ocean; in this case, we might expect the eddies to be associated with variability of the vertical velocity on horizontal length scales of around the grid scale and the first baroclinic Rossby radius. Here, we define the time-averaged-specific EKE (in J kg −1 ) from the mean and mean squared velocity components (updated at each time step during runtime) as where the mean specific kinetic energy MKE (also in J kg −1 ) is defined as and where the overbars denote a time mean. In order to avoid including the annual cycle of the mean flow in the EKE, we evaluate the MKE and EKE separately for the 12 individual months of the year over the 10year period, and then average over all months. Figure 2 shows the global means of the mean (MKE) and eddy (EKE) components of the kinetic energy and of the RMS vertical velocity w RMS , as functions of depth in the upper 2,000 m of the ocean, and averaged over the final decade of the integrations, along with their respective ratios to the same quantity in the control. Both MKE and EKE (Figures 2a-2d) decrease steadily with depth, and although the profiles themselves are similar, evaluating the fractional change relative to the control reveals a clear hierarchy in both MKE and EKE: as expected, increasing the viscosity reduces both, with the effect on the EKE a little larger, particularly in the upper 500 m. For comparison, we note that Delworth et al. (2012) compared the surface EKE in three GFDL coupled model configurations using the MOM ocean model with horizontal resolution of approximately 1°, 1/4°, and 1/10°: they found that in the highest resolution case the EKE was comparable in magnitude with observed estimates from satellite data. Plotting the surface EKE in the GO6 experiments on the same logarithmic scale (not shown) reveals that the sea surface variability is intermediate between that in the 1/4° and 1/10° GFDL experiments, and hence significantly closer to observations than in the corresponding 1/4° (CM2.5) GFDL model. This is consistent with the proposition that the standard 1/4° NEMO parameter set gives a model that is underviscous, relative to comparable model configurations.
In contrast to the kinetic energy in the horizontal flow, the spatial mean of the RMS vertical velocity (Figures 2e and 2f) increases with depth to a maximum at around 3,000 m, then reducing below that, consistent with a first vertical mode structure. The ordering of the profiles of the vertical velocity is also a little different from that for the horizontal kinetic energy: as with the KE, increasing the Smagorinsky term unambiguously reduces w RMS , but while doubling the fixed viscosity from the control actually increases w RMS , tripling it does decrease it at all depths. The mixed and smag_4 experiments give consistently lower w RMS than the other experiments, as is the case with the EKE. We note that there is substantial interannual variability of w RMS , even at 2,000 m depth, which reflects the fact that the transient vertical velocity field is mainly windforced, both locally (Ekman pumping) and remotely (NIGWs).
As demonstrated by Megann (2018), the instantaneous vertical velocity field contains much more energy in the 0.25° NEMO configuration than is captured even in a 5-day mean, because of the considerable inertial-time scale variability in this field, which can have amplitudes of over 100 m day −1 , even in the absence of tidal forcing. In this section, we discuss in a qualitative fashion the sensitivity of the small-scale structure in the vertical velocity field to the viscosity in two energetic regions; as shown in Figure 1 and Table 1, smag_2 has the lowest viscosity of the ensemble over most of the domain, while smag_4 (along with bi-lap_3 and mixed) has viscosities that are typically significantly higher than in the control. Figure 3 shows the instantaneous vertical velocity at 2,000 m depth in the North Atlantic, evaluated from the final restart file of the respective integrations, while Figure 4 shows the same quantity in the south-west Atlantic and Drake Passage. In both regions, there is a clear differentiation between the resolved features equatorward of 50°N and 50°S, respectively, and the complex, marginally resolved, or unresolved features poleward of these latitudes. We note further the almost plane wavefronts of the NIGWs propagating southward south of 35°N in the subtropical North Atlantic (Figure 3), which are not seen in the southwest Atlantic ( Figure 4). The mesoscale features are more energetic in the Drake Passage than in the North Atlantic (note the different color scales). The amplitudes of the near-inertial waves reduce by 10%-20% in the experiments with Smagorinsky viscosity but are relatively insensitive to increases in the fixed viscosity. By contrast, the eddies around the Gulf Stream and the smaller features along the western boundaries respond strongly to changes in the viscosity parameters, with increases in viscosity leading to reductions in the vertical velocity at the smallest length scales, at the same time as an apparent "clumping" of the features toward larger scales. The amplitude of spatial variability of the vertical velocity can clearly be seen to be reduced as the viscosity is increased from the control, with the largest changes east of Grand Banks. Although we should not necessarily expect a direct correspondence between instantaneous fields and long-time scale mean variance, the snapshots in Figures 3 and 4 are not inconsistent with the ordering of the time-mean profiles in Figures 2e and 2f: the control and smag_2 have the noisiest and highest-amplitude vertical velocity field, while mixed and smag_4 generally give the smoothest and lowest-amplitude fields.

The Cell Reynolds Number in the Experiments
The cell Reynolds number (Re) is a commonly used parameter in numerical fluid mechanics that represents the ratio of advective to viscous tendencies and is related (see e.g., Bryan et al., 1975) to numerical stability: in regimes where values are significantly above a critical threshold, unstable numerical modes tend to manifest at the grid scale. In the case of Laplacian viscosity, Re is generally defined as VΔ/ν, where V is the speed, Δ the cell horizontal dimension, and ν the viscosity; for stable flow, Re ≤ 2.  generalize the stability criterion for biharmonic viscosity (as used in the present model) as Re = VΔ 3 /ν ≤ 16. Although a high Re is neither a sufficient nor a necessary condition for numerical mixing, Ilıcak et al. (2012) suggest that "The use of a grid Reynolds number below 10 is a sensible and universal recommendation to avoid the saturation level of spurious mixing, but this would not guarantee that the remaining spurious mixing would be insignificant." A configuration in which Re exceeds the critical value in a considerable fraction of the cells may be characterized as underviscous, and their behavior at length scales around the cell dimension will be noisy. It is informative, therefore, to enumerate Re for the experiments described in this paper. Figure 5 shows the cell Reynolds number at 300 m depth, evaluated from a 5-day mean, in all the experiments, along with a cumulative histogram of this quantity in the global domain.  alone all have more grid points with Re > 16 than the control; of these, smag_4 is close to the control in this respect but has fewer cells with much larger Re. Figure 6a shows the ratio of effective to explicit diffusivity in the global domain in the GO6 control experiment, and this is quantitatively very close to that in GO5.0 (Figure 8a of Megann, 2018), as are also the ratios in the Atlantic and Pacific domains (not shown). The black contours are of the overturning streamfunction, while the bold black dashed lines show the maximum winter surface density, effectively delimiting the surface-forced portion of density space, which is masked out. The white areas in the interior indicate negative values, which are likely to result from nonlinearities in the equation of state. If the suppression of noise in the vertical velocity at the horizontal grid scale resulting from the increased viscosity does indeed lead to reduced numerical mixing, this will be evident as a clear reduction of the effective diffusivity in those runs which have a viscosity higher than the control. Figures 6b-6g show the ratios of the global κ eff to that in the control: with the enhanced viscosities, the effective diffusivity in most of the unventilated interior water masses is reduced by between 20% and 40% relative to that in the control, while in smag_2, by contrast, κ eff is 10%-20% higher than in the control over much of the interior. We note that these reductions extend up through the base of the subtropical cells, which can be seen as the near-surface recirculations between 40°S and 40°N; these typically extend to a depth of between 300 and 500 m, indicating that increasing the viscosity reduces the numerical mixing at the depth of the seasonal thermocline.  To more quantitatively examine the effect of increasing the viscosity on the numerical mixing, we evaluate κ eff ; the ratio of κ eff to the explicit diffusivity κ exp ; and the ratio of the effective diffusivity to that in the control integration, in each density class as a horizontal average over the global ocean in each integration ( Figure 7). We have removed the traces for water masses with potential densities higher than σ 2 = 37.1, as the relatively small volumes in these density classes add excessive noise to the ratio. We note that the profiles of the diagnosed mean diffusivity (Figure 7a) all vary in a very consistent manner with density, as do the profiles of the ratio of effective to explicit diffusivities (Figure 7b), confirming the robustness of the method. Looking at the ratio of κ eff to the control (Figure 7c), it is clear that there is a tendency for increasing viscosity to correspond to a reduction in the effective diffusivity: increasing the fixed viscosity by factors of 2 and 3 (red curves) gives commensurably lower diffusivity, as does increasing the Smagorinsky coefficient (green curves). The experiment smag_2 with the lowest Smagorinsky coefficient (dashed green line) has consistently the highest κ eff , with the control (black line) below it, while the lowest diffusivities occur with bilap_3, with bilap_2, smag_4, and mixed all having diffusivity lower than in the control in every density class.

Global Drifts and Biases
We would expect numerical mixing, diagnosed above in terms of enhanced values of the diapycnal diffusivity in intermediate to deep water masses, to contribute to unrealistic drifts in temperature and salinity in these water masses. If our hypothesis that the mixing is due to transient vertical velocities at close to the grid scale, which are controlled by increasing viscosity, is valid, we expect that drifts in the model tracer fields will reduce as the viscosity increases. We will only present here results for temperature drifts, since the near-monotonic temperature reduction with depth provides a conceptually simpler picture of diffusive mixing than is the case for salinity; in the latter case, the globally averaged salinity profile is dominated by a fresh halocline between 500 and 1,500 m depth, with relatively low vertical gradients of salinity below that, so it is difficult to interpret changes in terms of diffusion. Of course, other processes besides mixing (physical and numerical) will affect the time evolution of the water masses, in particular modes tend to exhibit a slow, long-term adjustment of the overturning circulation, and biases are furthermore sensitive to errors in surface forcing; drifts are therefore the sum of several tendencies, whose relative magnitudes will vary with depth and location, so attribution to a single source is difficult. This is further compounded by the observation that eddies can transport heat upward in the ocean along isopycnal surfaces (Griffies et al., 2015), contribution, of course, is to suppress the eddies and thereby to reduce the upward eddy transport of tracers. We shall, therefore restrict this section to a presentation of the model drifts and biases, and even if we are able to associate with confidence any improvements in realism to increased viscosity, it is not straightforward to attribute this unambiguously to the reductions in mixing presented in Section 4.3.  . Global means of (a) effective diffusivity κ eff , (b) κ eff divided by the explicit diffusivity κ exp , and (c) κ eff divided by κ eff in the control, and with potential density σ 2 as the vertical axis. Lines are suppressed for densities with σ 2 > 37.1 for clarity. but small changes below 2,000 m; by contrast, all the model runs have a strong cooling between about 200 and 2,000 m, which persists throughout all of the integrations. On top of this, though, there is again a clear ordering of the changes that is strongly correlated with the viscosity. Increasing the biharmonic viscosity has the clear effect of reducing the cooling at intermediate depths but enhances the net warm bias centered at around 200 m (although we shall show below that the latter results from a reduction in the cold bias at midlatitudes being compensated by a substantial warm bias in the Southern Ocean that is almost unaffected by viscosity). Conversely, decreasing viscosity leads to enhanced cold errors in the intermediate depth range: in smag_2, there is in fact near-uniform mean cooling at all depth levels from 300 m down to 4,500 m. Consistent with the results presented in the previous section for the effective diffusivity, smag_2 shows the largest rate of drift in each basin, followed by the control, while bilap_3 has the lowest drift, followed by bilap_2 and smag_4.
In with highest viscosity (Smag_4 and bilap_3) have the lowest downward heat flux, while Smag_2 has the highest, and over this depth range the experiments with higher viscosity lie closer to the profile for the EN4 climatology.
To illustrate the regional structure of the global mean biases discussed above, Figure 10a shows the temperature bias in the global domain in years 1996-2005 of the control experiment with respect to the EN4 climatology at 300 m depth, which is close to the depth at which the model has its maximum global mean biases (Figure 9), and is in the depth range typical of the seasonal thermocline. Significant basin-scale errors can be seen: a cold bias of up to −2.5 K across the tropics and the subtropical Atlantic, as well as in the Arctic; and a warm bias of comparable magnitude in the subpolar North Atlantic, the South Pacific and in the Southern Ocean. The warm mean biases at 300 m seen in Figures 9a and 9c can be seen to be associated with the latter large-scale warm features. The other panels show the temperature differences from the control at the same depth in the other experiments. Tripling the fixed viscosity (Figure 10c biases in the regions of the separated Kuroshio becomes worse as the viscosity is increased; the cold bias in the separated Gulf Stream strengthens, as does the warm bias in the retroflected Agulhas Current. There is therefore a trade-off between a reduction in the gyre-scale biases and a poorer representation of the separated boundary currents, where increasing viscosity may have an adverse effect on the dynamics of the model. The salinity biases at 300 m depth ( Figure 11) are overall similar in distribution to those in the temperature, with fresh biases of up to 0.5 psu in the tropics, subtropics and the Arctic, and salty errors across the northern subpolar regions and around the path of the ACC. Increasing viscosity generally again leads to reductions in salinity bias, but the effect is about half that on temperature and less consistent: as stated earlier, reducing mixing is not expected to robustly lead to reductions in salinity biases.

Circulation Indices and Mixed Layer Depths
For the reductions in numerical mixing and in the interior drift in the model that arise from increasing the viscosity to be considered as improvements in the model performance, these should not be countered by adverse effects on the large-scale circulation. We present results for three key volume transports: first the transport through the Florida Strait, which represents the main northward component of the upper branch of the Atlantic meridional overturning circulation (AMOC); the ACC as represented by the flow through Drake Passage; and the AMOC itself, evaluated at 26°N. We shall also discuss the meridional ocean heat transport (MOHT) at 26°N, and the volume transport and maximum velocity in the Pacific Equatorial Undercurrent (EUC), both of the latter two indices evaluated at 155°W.  annual metrics in the last 10 years of each experiment, while Figure 12 shows time series of each of these, along with the ratios of the same quantity in each of the experiments to that in the control.
Although there is an interannual variability of 1-2 Sv in the simulated Florida Strait throughflow (Figure 12a), the means are within 1 Sv of one another, and within 2% of the transport of 35.52 Sv in the control. We do not consider, therefore, that there is a strong sensitivity of the Florida Strait transport to the viscosity. Cable measurements of the transport through Florida Strait carried out under the Rapid program (https:// www.rapid.ac.uk/rapidmoc/) estimate 31.4 ± 3.1 Sv; this is not inconsistent with the simulated transport in the mid-2000s in the GO6 integrations, although the mean in the simulations is generally 2-3 Sv higher than the observations.
In the control experiment, the volume transport in the ACC (Figure 12b

Table 3 Mean Volume transports in Sv Through the Florida Strait (FS), Drake Passage (DP), and in the Atlantic Meridional Overturning Circulation (AMOC) at 26°N in the last 10 Years of the Model Experiments, and the Ratios of these to the Values in the Control Experiment
with changes in the fixed viscosity having a stronger effect than changing the Smagorinsky coefficient. This general tendency is consistent with current understanding of Southern Ocean dynamics (e.g., Hallberg & Gnanadesikan, 2006), where the transport of momentum by mesoscale eddies is thought to reduce the ACC transport; damping the eddies will therefore tend to increase the transport. It is worth mentioning that the mechanism of Johnson and Bryden (1989), in which eddy form drag transfers momentum downward from the wind stress, would by contrast predict a reduction in the ACC transport if the eddies were suppressed by elevated viscosity; however, Figure 4 suggests that substantial eddy activity remains in the Drake Passage region even in the Smag_4 case, and this may well be enough to maintain the isopycnal deformations required to generate form drag. The sensitivity to viscosity in our ensemble is opposite in sign to that reported by Jochum et al. (2008), but we note that the model used in the latter study does not have a resolution high enough for eddies to be represented in this region, so the first-order effect of increased viscosity in damping the mean flow is likely to dominate. We also note that the effect of reducing numerical mixing on the overturning circulation is likely to be progressive and become more visible in a longer simulation.
We might expect changes in the mixing at midlatitudes to have at least a small effect on the overturning circulation: Ito and Marshall (2008) and Nikurashin and Ferrari (2013) hypothesize that the depletion of deep waters by diabatic processes exerts a strong control on the global overturning circulation. This would predict that the higher numerical mixing seen in the experiments with lower viscosity would correspond to stronger overturning. The overturning at 26°N in the control experiment, peaking at almost 26 Sv in 1995, and averaging 22.59 Sv in 1996-2005, is unrealistically strong when compared with the observations made with the Rapid array at the same latitude (Smeed et al., 2018), which gives an estimate for the overturning strength of around 17 Sv between 2004 and 2016. In the GO6 experiments, however, there is only weak overall correlation between the Atlantic overturning circulation and the viscosity: the smag_4 experiment has consistently the weakest overturning at around 3 Sv, or around 12%, weaker than in the control, but the overturning strengths in the other experiments are not significantly different from that in the control. This means that the above hypothesis is not robustly supported by our results: the numerical contribution to the diapycnal mixing, at least, is not a good predictor of the overturning strength.
We note that the contours of the Atlantic overturning streamfunction against potential density and latitude on each panel of Figure 6 confirm that the overturning circulation is only weakly sensitive to the viscosity. First of all, the peak value at around 50°N and σ 2 = 36.75 lies between 24 and 25 Sv in all the experiments, with little change in its latitude and density or in the magnitude of the upper cell: only in Smag_4, is there a suggestion that there is any reduction in overturning strength. There is a hint that the lower cell strengthens weakly with reducing viscosity: in Smag_2, there is an increase in its strength by about 2 Sv, and in bilap_3 it is noticeably weaker by about 1 Sv.
The MOHT is clearly directly related to the volume transport, but also to the temperature stratification, so the sensitivity of the MOHT to the viscosity, through its effect on the numerical mixing, is of interest. Figure 12d shows that the sensitivity of the MOHT to the viscosity is qualitatively very similar to that of the AMOC (Figure 12 We conclude that, while the model has an unrealistically strong overturning circulation, its heat transport is weaker than that observed and that increasing the viscosity worsens this bias, albeit only by a few percent. In the context of the Atlantic overturning circulation, it is of interest to discuss the sensitivity of the wintertime dense water formation in the North Atlantic to the viscosity. In Figure 13, we show the March mixed layer depth, averaged over the last 10 years of integration, in the control, as well as the difference between this quantity in selected experiments and that in the control. We may compare the MLD in the control with climatological March mixed layer depth from the data set of Montégut et al. (2004) Table 4. The standard deviations of all three indices are between 4% and 7% of the respective mean, but the standard deviation of the total mixing volume is only 2.5% of the mean, suggesting that the first-order response to increasing viscosity is for the mixing to transfer from the Labrador Sea to the Nordic Seas, without a large change to the total mixing volume. This is consistent with the low sensitivity observed above for the AMOC strength. The reason for the change in spatial pattern is not immediately obvious, although it may be related to the reduction in grid-scale noise in the Labrador Sea seen in Figure 3 with increased viscosity: if these features lead to more numerical mixing, this will increase the rate of effective convection; alternatively, increased eddy activity at the edges of Labrador Sea may transfer water mass properties toward the center of the sea where most of the mixing occurs. At the same time, we note that these changes are quite small compared with the variability of mixing from one winter to another. The deepest mixed layers in the control, located west of Drake Passage as far as 150°W, are deeper than those observed, but elsewhere the model reproduces the austral winter mixed layer quite acceptably. Increasing the viscosity generally shallows the mixed layer on the flanks of the ACC but deepens it in the middle of the ACC path, and in the aforementioned region of the south-east Pacific the effect of increasing the viscosity is to mitigate the excessively deep winter mixing. The effect of increased viscosity on the boreal summertime mixed layer is not clear from the color scale in Figure 14, but narrowing the range (not shown here) reveals that the summer mixing in the subtropical gyre is less than 20 m in the model, in contrast to between 25 and 45 m in the climatology, while in the subpolar gyre the simulation mixes to over 70 m, about 50% deeper than in the climatology. Increasing viscosity tends to reduce the latter bias by up to 10 m but does not significantly affect the too-deep summer mixing in the former region.

Equatorial Pacific Metrics
A region where lower values of viscosity might be expected to improve the realism of simulations at this resolution is the Equatorial Pacific; Tropical Instability Waves (TIWs) typically have wavelength of around 1,000 km but are associated with fronts and other features in the surface temperature field down to below 1 km scale (e.g., Marchesiello et al., 2011). The filaments generated by the TIWs stir the cool upwelled waters on the Equator with the warmer surface waters at higher latitudes and tend to reduce the cold equatorial bias seen in lower-resolution models (e.g., Jochum et al., 2008). Megann and New (2001) found that increasing the harmonic Smagorinsky viscosity coefficient in an ensemble of regional isopycnic-coordinate model simulations with resolutions between 0.3125° and 1.25° led to a noticeable sharpening of the equatorial thermocline at 155°W and to a narrower and faster core of the EUC but did not result in a significant reduction in the surface cold bias. It is therefore of interest to ask whether increasing the viscosity in the present ensemble causes any reduction of realism in this region.
As can be seen in Figure 1, the Smagorinsky scheme substantially increases the viscosity in the tropics, where there is strong shear along the energetic zonal current systems; in Smag_4, the viscosity in the eastern Pacific within a degree or so of the Equator is more than 4 times that in  the control. The TIWs themselves are reasonably well represented in the ¼° GO6 simulations, and there is little noise in the SST field at the grid scale, even in the experiment with the lowest viscosity. It can be argued, therefore, that the elevated viscosities around the Equator produced by the Smagorinsky scheme are unnecessary and may lead to reduced realism in this region. To assess the variability in the simulations without aliasing the seasonal cycle, we quantify their amplitude by evaluating the mean variance of the SST on the Equator between 130°W and 110°W in December, when the variability in this region in the simulations is strongest, over the final decade of the integrations (Table 5). There is a weak but robust dependence of the amplitude on the viscosity: the standard deviation of the ensemble is 2.2% of the mean, with the lowest variability seen in bilap_3 and the highest in smag_2. To enable comparison with the results of Megann and New (2001), we define the thermocline stratification as the mean separation at 155°W of the 15°C and 26°C isotherms on the Equator. This again is quite insensitive to the viscosity, having a standard deviation across the ensemble of 1.2% of the ensemble mean, although it is largest in bilap_3 and smallest in smag_2, confirming that higher viscosities indeed tend to sharpen the thermocline, even if the effect is marginal. Finally, we evaluate the mean SST bias with respect to the CCI climatology ( [1996][1997][1998][1999][2000][2001][2002][2003][2004][2005] as between 1.5°S and 0.5°N and between 135°W and 100°W. In the control, the cold tongue is 0.56°C too warm; the error increases quite reliably with the fixed viscosity to 0.64 in bilap_3 and also increases with the Smagorinsky coefficient, although starting from a lower level in Smag_2. We evaluate the annual mean volume transport in the Pacific EUC at 155°W (Figure 12e), and the maximum eastward velocity in the core of the EUC at the same longitude (Figure 12f), which in this model occurs at about 110 m depth. The experiments with enhanced fixed viscosity, as well as the mixed experiment, lead to small but consistent increases in both transport (1.3% of the mean) and maximum velocity (0.9% of the mean), but in the experiments with Smagorinsky viscosity alone neither metric is distinguishable from that in the control. Given that the EUC transport varies over the annual cycle by more than a factor of 2, we conclude that the EUC speed and transport are not significantly affected by changing the viscosity parameters. Storkey et al. (2018) showed that in the ¼° GO6 configuration (identical to our control experiment), the seasonal cycles of the hemispheric sea ice extent lay within the error bounds of the HadISST1.2 observational estimates (Rayner et al., 2003) with the exception that the maximum summer sea ice retreat was more extensive than those observed: for 3 months, the minimum extent in both the Arctic and the Antarctic was only about 50% of that in the observations. The Arctic sea ice volume was substantially lower than that in the PIOMAS reanalysis (Zhang & Rothrock, 2003) in all months, and in the summer months the volume was less than 25% of that in PIOMAS. In Figure 15, we show the annual cycle of Arctic and Antarctic ice extent, defined, as in Storkey et al. (2018), as the total area in each monthly mean of cells in which the ice concentration is higher than 15% estimates already mentioned (we do not include observational estimates for the Antarctic ice volume, as the data coverage there does not justify this). It can be seen that the annual cycles in the seven experiments are almost indistinguishable, with the exception that in the Arctic the excessively low ice extent and volume are mitigated, albeit by a very small extent, by increasing the bilaplacian viscosity (bilap_2 and bilap_3): in both the latter experiments, the minimum extent increases by 12%-14% but is still below the lower error bound of the HadISST data. In these two experiments, the low bias in the Arctic summer sea ice volume is reduced by about 5%, but the minimum volume is still well outside observational limits. We conclude that the sensitivity of the sea ice cover in GO6 to the viscosity is rather weak, although increasing the fixed viscosity reduces the bias.

Quantitative Correlations
To confirm the robustness of the effects of changing the viscosity, we evaluate correlations of model metrics based on the fields discussed above with the mean viscosity (see Table 6). Several definitions of the latter quantity were calculated, including the global depth-weighted mean (as listed in Table 1); the mean at 2,000 m and the surface layer mean viscosity, and the largest correlations with all metrics except for the AMOC and the Drake Passage transport were obtained with the surface viscosity, the latter two being most MEGANN AND STORKEY  strongly correlated with the viscosity at 2,000 m. We will use a Pearson critical value of 0.669 (for N − 2 = 5 degrees of freedom) to judge the significance of each correlation.
The global effective diapycnal advective velocity, the effective diapycnal diffusivity, and the Florida Strait transport are strongly correlated with the viscosity, with magnitudes of correlation higher than 0.85. The EKE at 2,000 m and the mean effective diffusivity are both strongly negatively correlated with viscosity, as is the diapycnal advective velocity g adv (whose values are negative): all of these metrics would be expected to reduce in magnitude when the grid-scale noise is suppressed. The RMS vertical velocity at 2,000 m is negatively correlated with the viscosity, but less strongly than the aforementioned metrics: although increasing the viscosity unequivocally reduces the mean EKE, its effect on the vertical velocity variance appears to be sensitive to the choice of viscosity scheme. As shown in Figure 2f, increasing the viscosity through the Smagorinsky coefficient, or by combining the two schemes in the mixed experiments, leads to unambiguous reductions in the w RMS , but increasing the fixed viscosity appears to have little coherent effect. As noted in Section 4.1, though, the considerable interannual variability in the vertical velocity is likely to contribute to the lower correlation coefficient for the RMS vertical velocity. The much more robust correlation between the viscosity and g adv implies that weighting w by the stratification is indeed rather more informative on the potential for vertical motions to lead to mixing. Of the Equatorial Pacific metrics, the SST variance and the SST bias in the cold tongue are both strongly correlated with viscosity: as viscosity is increased, the variance reduces (in other words the TIWs become weaker) and the bias increases, and the signs of both these correlations infer that the realism of the simulations reduces with larger viscosities. The stratification in the equatorial thermocline is not significantly affected by viscosity.
The correlation of the temperature change in the 500-1,500 m depth range of 0.624 is marginally significant, but if the outlier point corresponding to the mixed experiment is removed, the correlation increases to more than 0.90. Although the Florida Strait throughflow has a strong negative correlation with the mean viscosity, it has a very narrow range (±3% of the mean), confirming the remark made in Section 5.2 that it is only weakly dependent on the viscosity. The spin-down of the Drake Passage transport in the control, as discussed in Section 5.2, is significantly reduced as viscosity is increased, with a positive correlation of 0.93. There is weak sensitivity of the AMOC at 26°N to the viscosity, and removing the outlier data point (the mixed experiment) does little to change this. Similarly, there is no significant correlation with the total winter mixing volume in the North Atlantic, although the transfer of mixing from the Labrador Sea to the Nordic Seas is significant, albeit weak. The heat transport at 26°N is, by contrast, strongly negatively correlated with the viscosity, although the standard deviation of this metric over the ensemble is less than 2% of the mean transport, and so the sensitivity may be considered negligible in the context of the much stronger interannual variation seen in Figure 12d. There a weak correlation between the RMS 300 m temperature bias and the viscosity, but the RMS salinity bias at the same depth has a strong positive correlation with the viscosity; despite this, the differences in SSS bias between experiments are still only of the order of 10% of the bias in the control, so the sensitivity of the SSS to the viscosity may be judged to be weak.

Summary and Discussion
We have used an ensemble of the GO6 global NEMO model configuration at 1/4° resolution to investigate the sensitivity of the effective numerical mixing to the viscosity, as well as the effects on the large-scale behavior of the model. The ensemble includes three experiments using the default biharmonic viscosity scheme, with the coefficient (which tapers as the inverse cube of the grid size) increased by factors of 2 and 3 relative to that of the control; three experiments with the fixed viscosity replaced by the Smagorinsky deformation-dependent viscosity, again with three different values for the respective parameter; and, finally, an experiment combining the two schemes. The effect of the viscosity was evaluated as a set of metrics, evaluated over the final decade of a 30-year integration of each of the experiments: these included the horizontal EKE; the RMS vertical velocity; the effective advective diapycnal velocity (the RMS vertical advection of density, weighted by the vertical density gradient, used as a measure of diapycnal fluxes from the transient vertical flow); the cell Reynolds number Re, a measure of the noisiness of the horizontal flow at the grid scale, and an indicator of the potential for computational instabilities; the effective diapycnal diffusivity as defined by Lee et al. (2002); and temperature drifts and circulation indices.
Increasing the viscosity through either scheme reduces the eddy KE with a high level of predictability, and also the RMS vertical velocity, albeit with a lower correlation. In the control experiment, the cell Reynolds number is found to exceed the critical value in around a third of the grid cells on the surface level; at the lowest value of the Smagorinsky coefficient, this increased to 86%, while tripling the fixed viscosity reduced the fraction to 13%. There is a robust correlation between both the effective diffusivity and the effective diapycnal velocity and the mean viscosity, with higher viscosities corresponding to lower diffusivity and a reduced rate of density advection by the vertical velocity.
In contrast to climatological data which suggest a near-surface warming over the last few decades, with no significant interior cooling, all of the experiments cool markedly between 250 and 2,000 m depth, with a maximum cooling at around 500 m of between 0.15°C and 0.4°C. The magnitude of the cooling in the ensemble is negatively correlated with the mean viscosity, with increasing viscosity giving a reduced cooling, and hence a more realistic trend. While the reductions in the temperature drift below 2,000 m are consistent with a reduction in downward effective heat flux, the situation at shallower depths is more complicated: increasing viscosity tends consistently to reduce the biases at 300 m over most of the ocean, but the largescale mean effective heat flux is upward. We may, therefore, conclude that increasing viscosity tends to improve the realism of the simulations, but it is difficult to interpret these sensitivities unambiguously as being due to a reduction in mixing. In some regions, for instance the subtropical and tropical Atlantic below the seasonal thermocline, it is possible to associate reduced mixing with reductions in bias and drift, but in many other regions a combination of processes, such as eddy heat transport, biases in the surface forcing or spin-up of the overturning circulation, may compete and at least partially cancel. Griffies et al. (2015) describe the net effect of transient mesoscale eddies in transporting heat upward in the ocean, opposing the tendency for the processes of advection by the mean flow and vertical mixing to move heat downward. We expect this to be the case in the GO6 simulations described here, since at low latitudes and midlatitudes eddies are reasonably well represented. It is likely that changes to the viscosity settings will change the balance between the advection by the mean flow and the advection by eddies by damping one or other preferentially, and the results presented here do not distinguish between the two contributions to heat transport. In addition, while the results we have presented here show that changes in viscosity have a robust influence on the global diapycnal mixing and hence on the large-scale climatically relevant fields, we admit that there are questions that the application of the mixing analysis to a global model cannot address satisfactorily. An idealized model in a small domain can simulate specific flow and forcing regimes and is a valuable tool for addressing more precisely targeted questions. For example, such an approach would illuminate the effect of the grid-scale features along the western boundaries poleward of 40°N, as described in Section 4.1, in modifying the local water mass properties.
If increasing the viscosity led to useful reductions in internal temperature and salinity drifts, but at the cost of an inferior simulation of other ocean metrics such as the large-scale circulation and increased surface biases, it would be considered unacceptable. We have demonstrated that the surface temperature and salinity biases, as well as the transports though Florida Strait, the strength of the AMOC and the Atlantic Ocean meridional heat transport, are only weakly sensitive to the viscosity, while increasing the viscosity in fact reduces the rate of spin-down of the ACC, which we ascribe to a suppression of the marginally resolved eddies that tend to transfer momentum away from the ACC. In addition, increased viscosity reduces the unrealistic summer sea ice retreat in both the Arctic and the Antarctic, while having only a weak effect on ice cover over the rest of the year. The question then arises of whether there are significant trade-offs for the improvements described above.
Comparison of Figures 2c and 7c suggests that, in the experiment with the largest Smagorinsky coefficient, a typical reduction of 15% in the effective diffusivity comes with a commensurate attenuation of EKE in the 1,000-2,000 m depth range. Raising the viscosity beyond the cases described in this paper will therefore continue to damp the eddy variability toward what is seen in non-eddy-permitting configurations, and in a coupled model configuration a reduction in surface variability will certainly have an adverse effect on ocean-atmosphere coupling. Jochum et al. (2008) show that the use of no-slip boundary conditions leads to excessive values of the Smagorinsky viscosity close to western boundaries, which adversely affect the circulation in their lower-resolution simulation. The Smagorinsky formulation is this sense a blunt instrument; Figure 1d shows that in the C Msmag = 4 simulation the largest viscosity values occur mainly on the western boundaries between 40°N and 40°S and along the Equator, where flows are already quite smooth at this resolution. Given that the highest viscosities in this case can locally be more than 5 times that in the control, it may be that dynamics in these regions are unrealistically damped, without being compensated by further reduction in the cell Reynolds number at higher latitudes. Indeed, we have shown that the variance of the surface temperature associated with TIWs, as well as the realism of the cold tongue, is reduced-albeit weakly-as the viscosity increases. We therefore conclude that the viscosity provides a useful control on numerical mixing that can be used in ocean-only and coupled configurations and-provided the reduction in eddy KE and the perhaps unnecessarily large local viscosities near western boundaries and along the Equator are acceptable for a given application-the considerable advantages of increasing viscosity within the range investigated here outweigh the potential disadvantages.
The question inevitably arises as to which combination of parameters is optimal: of the experiments discussed here, those with the biharmonic viscosity A Lm tripled with respect to the default; with the Smagorinsky parameter C Msmag set to 4.0; and the combination of the default background value of A Lm = −1.5 × 10 11 m 4 s −1 and C Msmag = 3.0 give comparable improvements over the control. Increasing the fixed viscosity alone is impractical, in that it demands reductions in the time step by up to a third, and requires a carefully controlled spin-up strategy in order to be numerically stable, while we note that the mixed scheme is not available in the standard v3.6 or v4.0 releases of NEMO. Using the higher value of the Smagorinsky coefficient is therefore the most practical and efficacious way to use viscosity to control the numerical mixing. We may conclude from Figures 2 to 5 that the choice of C Msmag = 4.0 offers a good compromise between reducing the numerical mixing and excessive damping of EKE, while retaining a comparable spread of cell Reynolds number to that in the control. We can speculate as to whether further increases in viscosity (assuming numerical stability) would continue to improve the simulations: we note that an integration of the GO6 configuration on a 1° global grid with laplacian viscosity (which may be considered comparable to an integration at higher resolution but with viscosity high enough to suppress all mesoscale variability) showed a much reduced temperature drift, consistent with a lower rate of numerical mixing. In the latter case, of course, the improvement comes at the cost of much poorer representation of many of the featuresparticularly western boundary currents and strait throughflows-that are acceptably simulated at 1/4°.
The approach described here to reducing numerical mixing is based on suppressing the poorly resolved horizontal and vertical motions at around the model grid scale, which we have shown form a strong contribution to numerical mixing. We have proposed that their effect on mixing is chiefly through advection of tracers across isosurfaces, and this is supported by the strong negative correlation between the RMS vertical advection and the viscosity. The mechanism for this is via the diffusive component of the vertical advection scheme, by which transient vertical velocities mix the tracers downgradient, so if the advection scheme were replaced by a less diffusive scheme, we might expect the noise-like flow features to have a smaller effect on numerical mixing. Hofmann and Morales Maqueda (2006) and Hill et al. (2011) have confirmed that using higher-order vertical advection in a global z-coordinate ocean model can lead to significantly reduced numerical mixing. Other schemes, such as fourth-order TVD (Gottlieb & Shu, 1998) and piecewise parabolic (James, 2000), promise to improve the model performance further but have not yet been tested in a global NEMO, and their computational overheads have not yet been assessed. An alternative approach is to allow the coordinate surface themselves to be displaced vertically in response to wavelike or eddy-like vertical motions on time scales up to a few days, so reducing the numerical advection across coordinate surfaces. This is a fundamental property of density-coordinate or isopycnal models, which have very low numerical mixing by construction, and in the hybrid-coordinate HYCOM model (Bleck, 2002) such a mechanism (called "slack") is also implemented in the near-surface constant-depth layers. In z-coordinate models, this approach has only rarely been used: LeClair and Madec (2011) proposed a similar modification called  z to the vertical coordinate in NEMO and demonstrated its effectiveness in an idealized channel domain in suppressing by a factor of 5 the numerical mixing from simulated internal waves. Although this has yet to be comprehensively tested in a global domain in NEMO, it offers a promising approach to reducing numerical mixing in this model.
Although it is risky to extrapolate results from an ocean-only model to predict the behavior of a coupled ocean-atmosphere model, we expect that numerical mixing in the ocean component will have consequences for the rate of heat uptake simulated by such models under anthropogenic warming. We note that the GC3.1/UKESM1 platform (Sellar et al., 2019) submitted to the CMIP6 intercomparison project uses a GO6 ocean with a nearly identical parameter set to our control experiment (and, specifically, the same fixed viscosity) and we have shown here that the default viscosity scheme in GO6 has numerical mixing that is toward the higher end of the present ensemble. Several other CMIP6 model configurations have z-coordinate oceans with comparable horizontal resolutions to the ¼° of GO6, so it would be surprising if the rates of numerical mixing observed here were not typical of this generation of coupled models. Although the ocean configurations for the climate and Earth System models participating in CMIP6 have already been frozen, the results presented here offer a challenge to the next generation of climate simulations.

Data Availability Statement
The data used for the analyses and plots are available at https://zenodo.org/record/3933125.