1 Introduction

A comprehensive operational oceanography service today should consist of both observation data (remote sensing data, in situ measurements) and model data (forecasts and reanalysis) of both physical and biogeochemical parameters (She et al. 2016). Such a service exists both at European level (CMEMS – Le Traon et al. (2017)) and mostly also at national level with a greater focus on the immediate coastal waters (e.g. in Germany, Brüning et al. (2014); in the Netherlands, De Kleermaeker et al. (2012); in Portugal, Mateus et al. (2012)). In particular, at national level, there is often an additional downstream drift prediction component of great importance in crisis situations such as oil spills (Broström et al. 2011; Maßmann et al. 2014) or search and rescue operations (Breivik et al. 2013). Obviously the current data of the operational circulation models are the most important input for the drift prediction. Due to the lack of direct current measurements, operational current forecasts are regularly validated indirectly by comparing drift model results with observations from surface drifters (Callies et al. 2017). Unfortunately, such drifter experiments usually take place only under normal conditions, while drift predictions often have their highest relevance in extreme situations, whereas the quality of current forecasts during normal conditions cannot easily be transferred to extreme situations. Accidentally, instabilities resulting in unrealistic currents in the CMEMS-Baltic MFC NRT-forecast product (Le Traon et al. 2017) during storm events limited in time and space have been noticed, which were not detected by the very extensive, usual calibration and validation procedures of the Baltic MFC including the validation of extreme sea level events (CMEMS (n.d.); Golbeck et al. (2015)). More detailed analyses of these events showed that beneath current profiles especially diffusivity profiles were unrealistic, so that the turbulence model was considered more closely. It was found that despite explicit stability checks like a check for the CFL-criteria and checks for (linear) stability in the von Neumann sense, in certain situations, numerical instabilities occurred, which could only be avoided by additional non-linear stability and realizability checks within the turbulence closure scheme. Those checks are already described for turbulence closure schemes not using double diffusion in Umlauf and Burchardt (2005), so that these checks were extended to our turbulence model including double diffusion.

The paper is organized as follows: Section 2 describes the used model HBM and the model setups in detail. Section 3 gives an explanation of realizability in general (Section 3.1) and documents the added stability and realizability criteria (Section 3.2). In Section 4, both the physical and also the technical (ε-tests) validation of model results with and without the new criteria are shown, and finally the impact of the enhanced stability on the results of a downstream drift model is presented in Section 5. Main conclusions are provided in Section 6.

2 Model and setups

2.1 Numerical model

The physical model used in this study is HBM (HIROMB-BOOS model) described in Berg and Poulsen (2012) which is used both for operational forecasts (e.g. by the CMEMS Baltic MFC or the BSH), for reanalysis (Fu et al. 2012) and for research projects at various institutes, especially for the North Sea and Baltic Sea region (e.g. MeRamo:Neumann et al. (2018, in review), CLAIM (n.d.)).

HBM is a three-dimensional baroclinic ocean circulation model using Boussinesq approximations. The model is a further development of the operational circulation model BSHcmod (Dick et al. 2001). Like in BSHcmod, advection and diffusion are realized by a flux-corrected transport scheme, and the horizontal viscosity is parametrized by Smagorinsky (1963). In HBM the user has the choice between z-coordinates with free surface and so-called dynamical or generalized vertical coordinates (Dick et al. (2008), Kleine (2004)) and the possibility of a fully dynamical two-way nesting with any number of grids. The vertical mixing is realized by a two-equation k-ω turbulence model accounting for buoyancy-affected geophysical flows (Umlauf et al. 2003). By parametrization, the shear due to internal waves (Axell 2002), an estimate production in the surface layer from below and unresolved bottom shear due to tides (Canuto et al. 2010) are also taken into account. A detailed description of the used parameters resp. the parameter-making could be found in Berg (2012). In this study, the turbulence model is coupled to an algebraic second order closure scheme either based on Canuto et al. (2002) or based on Canuto et al. (2010). Both closure schemes consider double diffusion which is relevant in the Baltic Sea – one of the main application areas. For example, at the Baltic Sea Science Congress 2017, it was shown that double-diffusive instabilities may constitute a key mixing process in this region (Gillner et al. 2017). The latter scheme finally has been extended by additional stability and realizability checks. A detailed description of this extension can be found in Chapter 3.

2.2 Setups

During this study, two setups both running in operational mode and both covering the entire North- and Baltic Sea were used. Both setups were forced by atmospheric data from the operational atmospheric model of the German Weather Service (DWD) and run-off data from the operational run-off model E-hype operated at the Swedish Meteorological and Hydrological Institute. At the open boundary in the northern North Sea and in the English Channel, the water level has been set to the sum of surge data generated by BSH’s operational North Atlantic model and tides based on 19 partial constituents. Temperature and salinity at the open boundary were taken from the Janssen et al. (1999) climatology.

2.2.1 CMEMS setup

The CMEMS setup is using z-coordinates with free surface and consists of four nested grids:

  • “North Sea” with a horizontal resolution of 3 nautical miles and up to 50 vertical layers

  • “Wadden Sea” with a horizontal resolution of 1 nautical mile and up to 24 vertical layers

  • “Inner Danish Waters” with a horizontal resolution of 0.5 nautical miles and up to 77 vertical layers

  • “Baltic Sea” with a horizontal resolution of 1 nautical mile and up to 122 vertical layers

2.2.2 BSH setup

The BSH setup is using dynamical/generalized vertical coordinates and consists of two nested grids:

  • “North Sea/Baltic Sea” with a horizontal resolution of 3 nautical miles and up to 36 vertical layers

  • “German Coastal Waters” with a horizontal resolution of 0.5 nautical miles and up to 25 vertical layers

3 Realizability and stability

3.1 General explanation

Of course, models are always a simplification of reality. In particular, numerical models actually solve only mathematical equations, so unfortunately such models are not realizable in general, meaning that the quantities, which are positive by definition, may become negative due to the model implementation. Explicitly implemented realizability criteria in the form of limiters are intended to ensure that the physically relevant model state resp. the physically relevant parameter range is not left, so that the model does not predict unphysical values, such as, for example, negative variances, diffusion coefficients, length scales, tracer concentrations, masses and volumes.

In practice stability and realizability are often difficult to distinguish, both are local characteristics, which might be triggered (in different directions) by small numeric discrepancies (e.g. due to more or less aggressive compiler tuning flags). They can spread, disappear or even lead to model crashes. Therefore, the goal is to define explicit criteria that guarantee robust numeric.

3.2 Additional stability/realizability criteria

As described in Berg (2012), the three vertical diffusivities, namely, Km for momentum, Kh for heat and Ks for salinity can be expressed as

$$ {K}_i=2\frac{k^2}{\epsilon }{S}_i\ \mathrm{with}\ i=m,h,s $$
(1)

where k is the turbulent kinetic energy and ϵ is the dissipation rate of turbulent kinetic energy. The dimensionless functions Si are the structure or stability functions, which depend on three variables and can be written as

$$ {S}_i={S}_i\left({a}_m,{a}_h,{a}_s\right)\ \mathrm{with}\kern0.5em i=m,h,s\ \mathrm{and} $$
(2)
$$ {a}_m={\left(\tau \varSigma \right)}^2,{a}_h={\tau}^2{R}_h,{a}_s={\tau}^2{R}_s $$
(3)

Equation (3) describes the dimensionless shear number, heat number and salinity number with the dynamic dissipation time scale \( \tau =2\frac{k}{\epsilon } \), the mean shear Σ, \( {R}_h={\alpha}_T\frac{\partial T}{\partial z} \)and \( {R}_s={\alpha}_S\frac{\partial S}{\partial z} \). In this case, z is the water depth, T the water temperature, S the salinity and \( {\alpha}_{T,S}=-\frac{1}{\rho}\frac{\partial \rho }{\partial \left(T,S\right)} \) denotes the thermal and the haline concentration coefficient, respectively, where ρ is the water density. It is important to mention that, according to Canuto et al. (2002) and Canuto et al. (2010), the buoyancy production N2 can be expresses as follows:

$$ {N}^2={R}_h-{R}_s $$
(4)

In contrast to the structure functions (2), which are used in this study and which include double diffusion, structure functions without double diffusion are of the form

$$ {S}_i={S}_i\left({a}_m,{a}_n\right)\ \mathrm{with}\ i=m,n $$
(5)

where n stands for buoyancy and an = (τN)2 for the dimensionless buoyancy number. For the functions (5), the following stability and realisability criteria are explicitly stated in Umlauf and Burchardt (2005):

At first all vertical diffusivities and therefore all structure functions must be greater or equal a background value, which must be greater or equal zero, i.e.

$$ {K}_i\ge {c}_i\ge 0\ \mathrm{with}\ i=m,n\ \mathrm{and}\ \mathrm{constants}\ {c}_i. $$
(6)

The condition guaranteeing increasing effective vertical shear anisotropy with increasing dimensionless vertical shear number can be formulated as

$$ \frac{1}{2}{\partial}_{a_m}\left({K}_m\left({a}_m,{a}_n\right){a_m}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)\ge 0 $$
(7)

The value of an must be greater or equal the value \( {a}_n={a}_n^{\ast } \) which describes the value in shear-free convective conditions for the turbulence equilibrium, in which the buoyancy production G equals the dissipation rate ϵ, i.e.

$$ {a}_n\ge {a}_n^{\ast } $$
(8)

To prevent oscillation between two mathematically possible values of an, monotonicity of \( \raisebox{1ex}{${K}_n\Big({a}_{m,}{a}_{n\Big)}$}\!\left/ \!\raisebox{-1ex}{${a}_n$}\right. \) with respect to an must be insured for negative an:

$$ -{\partial}_{a_n}\left(\raisebox{1ex}{${K}_n\Big({a}_{m,}{a}_{n\Big)}$}\!\left/ \!\raisebox{-1ex}{${a}_n$}\right.\right)>0 $$
(9)

In Umlauf and Burchardt (2005), it is shown that (8) always implies (9), so that (9) is not an additional condition from a numerical point of view. Finally the velocity variances must be positive, with the only critical condition being

$$ <{w}^{\prime },{w}^{\prime }>\le 2k $$
(10)

where <w′, w′> is the vertical velocity variance.

In the original implementation in HBM, condition (6) is always fulfilled by definition. Conditions (7), (8) and (10), however, are not explicitly queried and—as it turned out—not always fulfilled. It is therefore necessary to extend these conditions to the used structure functions including double diffusion. While (10) is very easy to implement due to the numerical descriptions of Canuto et al. (2010) (<w′, w′> is listed explicitly in the equations), (7) is also relatively easy to expand to

$$ \frac{1}{2}{\partial}_{a_m}\left({K}_m\left({a}_m,{a}_h,{a}_s\right){a_m}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)\ge 0 $$
(11)

To generalize (8), we consider the shear-free convection case once for the case without vertical temperature differences (heat number = 0) and once for the case without vertical salinity differences (salinity number = 0). From this consideration, minimum values \( {a}_h^{\ast } \) and \( {a}_s^{\ast } \) result for the heat number or the salinity number. Due to reasons of symmetry, \( {a}_h^{\ast }={a}_s^{\ast }={a}_{h,s}^{\ast } \), so that (9) becomes

$$ {a}_h\ge {a}_{h,s}^{\ast },{a}_s\ge {a}_{h,s}^{\ast } $$
(12)

in the case with double diffusion.

4 Model results

The main motivation for the investigation of this problem was the occurrence of unrealistic surface currents during a storm event in the CMEMS-Baltic MFC NRT-forecast product version 1 (Fig. 1, right side), which did not occur anymore after including stability and realizability criteria described in Chapter 3 and the following change to product version 3 (Fig. 1, left side). Of course, these criteria were not the only difference between version 1 and version 3, so that the comparison of these versions is not suitable for demonstrating their influence. Therefore, this is demonstrated in the following by using the operational BSH setup during a storm event in December 2013, whereby the model runs only differ in their different turbulence closure schemes.

Fig. 1
figure 1

(taken from CMEMS (n.d.)) Surface currents for V3 and V1 products: direction are shown in arrows; magnitude is coloured; snapshot during strong wind event (13 December 2014 4 UTC)

Full size image

4.1 (Physical) model validation

By chance, temporary unrealistic surface currents during a storm event in December 2013 were also discovered in the operational setup of the BSH (Fig. 2, left column). The fact that the water levels at the tide gauges were simulated correctly led to the conclusion that the water transports as a whole were correct, but the turbulence scheme in this situation was locally not correctly solved or unstable. The assumption that this was a local instability was also confirmed by the ε-test (see Chapter 4.2).

Fig. 2
figure 2

Eastward component of surface currents before, during and after the peak of storm event on 5 December 2013. The black dot at the top left describes position 6.55° E/55.58° N

Full size image

Moreover, the assumption was confirmed by looking at the current and eddy diffusivity profiles at 6.55° E/55.58° N (Fig. 3, left column). Tests with Canuto et al. (2010)-based closure schemes (Figs. 2 and 3, middle and right column) show that all used closure schemes provide comparable current patterns before the peak of the storm and therefore in a more or less common situation (Figs. 2 and 3, first row). The diffusivity profiles of the closure schemes without explicit realizability and stability checks, on the other hand, already look partly nonphysical (Fig. 3, first row, left and middle column). In the described storm case, the diffusivity profiles then appear to become completely unstable, so that a clear stratification of the currents occurs at a depth of approximately 15–20 m. While enormous surface currents are directed towards the east, the water flows into the opposite direction at depth, so that the entire water transport fits again. Realistic profiles and current patterns in all situations are only achieved with a closure scheme that has been extended by stability and realizability criteria (10), (11) and (12) described in Chapter 3.2 (Fig. 2 and 3, right column).

Fig. 3
figure 3

Profiles of both eastward current component and eddy diffusivity at position 6.55° E/55.58° N before and during the peak of storm event on 5 December 2013

Full size image

4.2 ε-Tests/technical model validation

The ε-test is an originally technical test in which the results of short model runs (e.g. 24-h simulation) are compared, the only difference between these runs being the compiling of the source code. It is done with different compilers and/or with a different set of compiler flags. The results of these runs are compared point by point, and the maximum differences are the ε’s, where small ε’s indicate both technically and physically stable code that provides reliable and therefore portable and reproducible results. Large ε’s, on the other hand, indicate a stability problem in the code, which can have both technical and physical causes.

Within this study, the ε’s of a 12-h simulation between a model run compiled in optimization level O2 and a run compiled in optimization level O3 were analysed. This analysis showed without explicit stability and realizability checks local ε’s of the eastward current component in the order of 0.7 m/s (Fig. 4, left and middle), i.e. in the order of magnitude of the eastward current itself. This is a clear indication of (local) instability. It was also found that instability no longer occurs with the implementation of the additional stability and realizability checks (Fig. 4, right). In this case, the ε’s in the entire area were in the range 10−7.

Fig. 4
figure 4

ε’s of the eastward current component at the surface between two model runs, one with O2- and one with O3-optimized code, at the peak of the storm on 5 December 2013 at 12 UTC after 12-h simulation

Full size image

5 Results of downstream drift model

In the previous chapters, it has been shown that missing realizability criteria can lead to at least temporally and spatially limited instabilities and thus partly unphysical currents in certain situations. In order to demonstrate the impact of the obviously incorrect currents on the whole application range, we have carried out two drift calculation comparisons based on fictitious cases using the BSH’s operational downstream drift model SeaTrackWeb (Maßmann et al. 2014), whereby, on the one hand, currents generated without additional realizability criteria in turbulence closure based on Canuto et al. (2002) and, on the other hand, currents generated with additional realizability criteria in turbulence closure based on Canuto et al. (2010) were used as the forcing for the drift calculations.

The first scenario describes the drift of an object (e.g. a container, a buoy, a floating boat or ship or even a human body, which has gone overboard), while in the second case, the drifting of oil after an oil spill has been simulated, where 15,000 t of oil have spilled. Both fictitious cases “occurred” at the position 6.55° E/55.58° N on 5 December 2013 at 10:00 UTC shortly before the storm peak (the position is shown in Fig. 2). The drift was calculated 72 h in each case, i.e. until 8 December 2013, 10:00 UTC.

5.1 Drift of an object

The calculated drift paths (Fig. 5) show significant deviations especially at the beginning of the calculation. Thus, the distance of the calculated positions after only a few hours is approximately 10 km (Fig. 6). Since the instabilities in the underlying currents, which were calculated without an explicit realizability check, and therefore the greatest difference between the two current data sets, occurred at the beginning of the drift calculation, this is certainly an expected result. In the further course, the distance of the calculated positions decreases again, whereby in this calculation, it was always between 5 and 8 km (Fig. 6). If one also takes into account, the uncertainty that a drift calculation entails anyway, this means that the differences (the error resp.) caused by the missing realizability checks, could mean a question of life and death in a search and rescue operation.

Fig. 5
figure 5

Drift of an object over 72 h from 5 December 2013 10 UTC to 8 December 2013 10 UTC (the starting point is marked by a blue star, the end points are marked by blue circles) forced by the two different current data sets

Full size image
Fig. 6
figure 6

Distance of the calculated object positions of the two drift calculations

Full size image

5.2 Oil drift

Also in the case of the artificial oil spill, the results of the drift simulation, which was forced by the currents generated without explicit realizability and stability checks, differ significantly from those of the drift simulation, which was forced by currents generated with explicit realizability checks. In particular, after a 72-h simulation, the area contaminated by the oil is significantly larger in the case without realizability and stability checks (Fig. 7). In an emergency, the lack of accuracy of the drift forecast caused by the missing realizability checks could lead to very high additional costs, in any case, it means a poorer basis for resource planning.

Fig. 7
figure 7

Oil drift simulation results of the artificial oil spill after 72 h on 8 December 2013 at 10 UTC (the position of the oil outlet is marked by a blue star, the black and red crosses represent oil positions) forced by the two different current data sets

Full size image

6 Conclusion

In this study, an extension of realizability and stability criteria to turbulence closure schemes with double diffusion was presented, which were previously only known or at least documented for those without double diffusion. The lack of explicit realizability and stability checks led to spatially and temporally limited instabilities with nonphysical model (especially current) forecasts during storm or strong wind events, which were even found in operational products such as the CMEMS NRT product of the Baltic Sea.

Unfortunately, this problem is not limited to a short-term obviously wrong current product. Although the instabilities found in the circulation model were only of short duration and spatially strongly limited, it could be shown that the instabilities can be of considerable relevance for a downstream drift product, since they can significantly reduce the quality of this service. For operational services in particular, this is therefore a serious issue to be taken into account in the usually very extensive validation procedures. In order to guarantee the prediction quality at a constant good level, the validation of special individual events should always be performed in addition to the statistical validation. Special events are not only storms, which are accompanied by extreme water levels and strong currents. For example, sudden cold or heat bursts have the potential to cause short-term instabilities in the models too. As shown in this study, for the validation of special events, which are only of short duration, technical validation such as the ε-tests is suitable in addition to the physical validation. In any case, it should by no means be neglected in operational use but should rather increase in scope in the course of automation.