Introduction
Diapycnal mixing is a natural process whereby water parcels of different temperature and/or salinity get mixed. Â This process, as opposed to advection, is non-adiabatic (non-reversible) and is governed by diffusive mechanisms. The rate of diapycnal mixing is elevated in coastal waters relative to the open ocean and is responsible for the introduction of nutrients to plankton-containing waters at upwelling zones and within the mixed layer (nutricline/thermocline co-location).
Resolution and management of diapycnal mixing was a major issue in all of the early global ocean models. Early work by Bryan (1987) as well as more recent research has found that global circulation models are quite sensitive to the diapycnal mixing rate in modeling realistic, regional ocean structure and composition.
Models using the density coordinate system ($\rho$) have been used extensively since their creation to study the importance of diapycnal mixing within global ocean models. Unlike $\sigma$ or $z$ coordinate models, $\rho$ models can leverage adiabatic transport equations so that the only diapycnal mixing within the system is that which is explicitly formulated. The ability to alter, tune and compare the impact of a range of diapycnal mixing levels reaffirmed the importance of accurately describing and controlling the level of diapycnal mixing in global ocean models.
While $\sigma$ and $z$ coordinate models can have truly adiabatic advection schemes, there has been much progress in how to limit the diapycnal flux associated with their advection schemes. Since the schemes used in these coordinate models can only be adiabatic as a limit of convergence, much effort has gone into determining how closely schemes need to be to accurately describe diapycnal mixing. These so-called “quasi-adiabatic” improvements aim to limit the spurious diapycnal mixing to a level consistent with the level measured in the ocean (by Ledwell, for example).
Methods
To study diapycnal mixing, the simplest system seemed like a logical starting point. The MITgcm is a $z$-coordinate model used extensively in both coastal reasons and global circulation studies. The domain used is a 320 x 1 x 60 ($x, y, z$) grid with $\Delta x = \Delta y = \Delta z = 1m$ for the sake of simplicity. Both sides and the bottom are closed boundaries with a no-slip condition, while the surface is an implicit free surface.
For physical parameters, there was no wind stress and $\beta$ and $f$ were both set to 0. The bottom drag was linear ($d = 1.0 \times 10^{-4}$) and the viscosity was set to 0.1 so that all the schemes could be run under the same conditions. The time step was 1 s and time stepping was staggered. The initial water structure was S = 0 everywhere and $T=0$ for the left half of the domain while $T = 2$ on the right side (lock-exchange, see Figure 1). At $t=0$ the two water masses were released.
Results
While the velocity fields in each of the models were nearly identical, the treatment of the tracer varied widely. As the colder, denser water fell beneath the warmer water, regions of Kelvin-Helmholtz type instabilities formed which lead to advective stirring of the two water masses. For all of the model runs the explicit diffusion was set to zero; but as we know of $z$-coordinate models, the advection scheme determined the rate of implicit mixing.
As we can see in Table 1, the order accuracy of the schemes vary from 1st to 4th order with stencils of either 3 or 5 points. While computationally more expensive, it is clear that flux limited schemes were superior to DST schemes which were, in turn, superior to centered schemes in minimizing the spurious diapycnal mixing in this model. The tracer used below is temperature, but analogous results should be found for any passive tracer (e.g. salinity).
The centered schemes suffered strongly due to the overshoot and undershoot caused by the rippling effect of this class of scheme. Of these two schemes (Scheme 2 and 4), the lower order scheme lead to significantly more mixing yet did keep the under/overshoots more local to the fronts (Figure 1). The fourth order centered scheme caused large regions above and below the fronts to become striped with ripples (Panel b).
The centered schemes were also the least stable schemes numerically. Schemes 2 and 4 both required many rounds of tweaking of the viscosity before the model would run successfully. This is a direct result of their their sensitivity to obligation to under/overshoot the temperature.
Moving on to the non-flux limited DST schemes (Schemes 20 and 30), the level of under/overshoot artefacts is largely corrected with small ripples still present in the advected, frontal zones at the top and bottom (Figure 3). The increased stencil of Scheme 30 leads to a further reduction in ripples at the expensive of both CPU cycles and spurious diapycnal mixing. A large stencil is inherently more diffusive and this is clearly present in the center of the mixing zone where multiple different water parcels are present.
Through the introduction of a flux limiter, some of the spurious mixing of the larger stencils can be corrected for. Both Scheme 33 and 77 show major improvements in the under/overshoot associated with the centered schemes (Figure 4). In addition, the high rate of diffusion present in Scheme 30 has been limited. The most dramatic example of this improvement can be seen in the resolution of the turbulent eddies formed along the front. Both flux limited schemes show high levels of detail within the microstructure of the eddies while all previous schemes lead to unresolved mounds due to their diffusivities.
Finally, the comparison of Scheme 2 and Scheme 33 at a later time (timestep = 500) is enlightening (figure 5). The level of noise and scheme-induced error is quite significant for Scheme 2 (panel a) while Scheme 33 (panel b) shows none of the same features. The eddies are well resolved (Panel b) and the homogenous interiors of each of the two water parcels are just that: homogenized.
