Ocean Dynamics (2008) 58:237–246 DOI 10.1007/s10236-008-0145-6 On the mathematical stability of stratified flow models with local turbulence closure schemes Eric Deleersnijder · Emmanuel Hanert · Hans Burchard · Henk A. Dijkstra Received: 3 June 2008 / Accepted: 22 August 2008 / Published online: 19 September 2008 © Springer-Verlag 2008 Abstract Occasionally, numerical simulations using lo- cal turbulence closure schemes to estimate vertical tur- bulent fluxes exhibit small-scale oscillations in space, causing the eddy coefficients to vary over several orders of magnitude on short distances. Theoretical develop- ments suggest that these spurious oscillations are essen- tially due to the way the eddy coefficients depend on the vertical gradient of the model’s variables. An instability criterion is derived based on the assumptions that the artefacts under study are due to the development of small-amplitude, small time- and space-scale pertur- bations of a smooth solution. The relevance of this criterion is demonstrated by applying it to a series a clo- Responsible Editor: Tal Ezer E. Deleersnijder Centre for Systems Engineering and Applied Mechanics, Louvain School of Engineering, Université catholique de Louvain, Louvain-la-Neuve, Belgium e-mail: ericd@uclouvain.be E. Hanert (B ) Department of Meteorology, University of Reading, Reading, UK e-mail: emmanuel.hanert@uclouvain.be H. Burchard Leibniz Institute for Baltic Sea Research Warnemünde, Warnemünde, Germany e-mail: hans.burchard@io-warnemuende.de H. A. Dijkstra Institute for Marine and Atmospheric research Utrecht, Utrecht University, Utrecht, The Netherlands e-mail: dijkstra@phys.uu.nl sure schemes, ranging from the Pacanowski–Philander formulas to the Mellor–Yamada level 2.5 model. Keywords Marine modelling · Vertical mixing · Turbulence closure · Stability 1 Introduction A number of present-day models of geophysical and en- vironmental fluid flow have recourse to a Fourier–Fick parameterisation of the vertical fluxes of momentum, heat, and dissolved constituents: the flux of the relevant quantity is expressed as the product of its vertical gradi- ent and a suitably defined eddy coefficient. The simplest closure assumption consists in assuming that the eddy coefficients are constant, an approach which was shown to be inappropriate for studying stratified marine and oceanic flows (e.g. Ruddick et al. 1995; Goosse et al. 1999). To consider flows with vertical density contrasts, parameterisations were suggested in which the eddy coefficients are functions of the Richardson number (Munk and Anderson 1948; Pacanowski and Philander 1981). Clearly, this option is better, but is not always able to properly capture the main processes governing the evolution of turbulent motions. This is why more complex models are often preferred, such as the k, k − ε or Mellor–Yamada models (Mellor and Yamada 1982; Rodi 1987; Burchard 2002a). The original Mellor–Yamada level 2.5 closure scheme (Yamada 1977) is known to be prone to in- stability. As mentioned by Mellor and Yamada (1982), “for some model simulations, a discontinuity in velocity could develop and persist” causing large-amplitude,