Simulation of the Arctic—North Atlantic Ocean Circulation with a Two-Equation K-Omega Turbulence Parameterization

Abstract

The results of large-scale ocean dynamics simulation taking into account the parameterization of vertical turbulent exchange are considered. Numerical experiments were carried out using $k-\omega$ turbulence model embedded to the Institute of Numerical Mathematics Ocean general circulation Model (INMOM). Both the circulation and turbulence models are solved using the splitting method with respect to physical processes. We split $k-\omega$ equations into the two stages describing transport-diffusion and generation-dissipation processes. At the generation-dissipation stage, the equation for $\omega$ does not depend on k. It allows us to solve both turbulence equations analytically that ensure high computational efficiency. The coupled model is used to simulate the hydrophysical fields of the North Atlantic and Arctic Oceans for 1948–2009. The model has a horizontal resolution of 0.25${}^{\circ }$ and 40 $\sigma$-levels along the vertical. The numerical results show the model’s satisfactory performance in simulating large-scale ocean circulation and upper layer dynamics. The sensitivity of the solution to the change in the coefficients entering into the analytical solution of the $k-\omega$ model which describe the influence of some physical factors is studied. These factors are the climatic annual mean buoyancy frequency (AMBF) and Prandtl number as a function of the Richardson number. The experiments demonstrate that taking into account the AMBF improves the reproduction of large-scale ocean characteristics. Prandtl number variations improve the upper mixed layer depth simulation.

1. Introduction

The simulation of momentum, heat and salt vertical turbulent exchange is very important for ocean general circulation models (OGCMs). In OGCMs, vertical mixing is parameterized using a second-order differential operator with variable exchange coefficients $K_U$, $K_T$ and $K_S$ [1,2,3,4,5]. The two basic approaches are used to determine these coefficients: (1) they are defined as functions of stratification and velocity shift or the local Richardson number [1,6,7,8]; (2) they can be found using additional turbulence models [3,9,10,11,12,13,14,15]. The first approach is more simple and efficient from the computational point of view. This approach is traditionally used in numerical simulation of large-scale ocean circulation [6]. The second one gives better results from the physical point of view [15], however, the computation cost increases. Therefore, when using the second approach, special attention should be paid to the construction of an effective numerical technique for solving turbulence equations.

To simulate large-scale ocean circulation we use the traditional Reynolds averaged equations with parameterized vertical turbulent fluxes $\overline{u^{\prime }{w}^{\prime }}$, $\overline{v^{\prime }{w}^{\prime }}$, $\overline{T^{\prime }{w}^{\prime }}$, $\overline{S^{\prime }{w}^{\prime }}$, where $u^{\prime }$, $v^{\prime }$, $w^{\prime }$, $T^{\prime }$ and $S^{\prime }$ are turbulent disturbances for the velocity components, temperature and salinity. With the framework of K-hypothesis, the turbulent fluxes are related to the local gradients of the corresponding mean quantities by eddy viscosities and diffusivities. With this parameterization, the turbulent coefficients depend on the both turbulence characteristics and the stability functions, which are also derived from the equations for turbulent fluxes [15,16,17,18].

The turbulence models used for parameterization of turbulent viscosity and diffusivity in OGCMs are based on the equations for the turbulent kinetic energy (TKE) and the characteristic needed for closure of the system. Such characteristics may be following ones: $ϵ$ is specific dissipation rate of TKE by the water viscosity; $\omega$ is frequency characteristic of the turbulence decay process [19] or, briefly, the turbulent dissipation frequency (TDF); and l is turbulence spatial scale. The corresponding two-equation turbulence models or closure schemes are referred to as $k-ϵ$, $k-\omega$ and $k-kl$. In respect to formulation, they are similar and can be combined as a “generalized scale method” [15,20]. The characteristics $ϵ$, $\omega$ and l are connected by algebraic relations, including the TKE k and the stability function for neutral stratification $c_S^0$: ${\displaystyle ϵ={\left(c_S^0\right)}^3\frac{k^{3/2}}{l}}$, ${\displaystyle \omega =\frac{ϵ}{{\left(c_S^0\right)}^{4}k}}$, ${\displaystyle l=\frac{k^{3/2}}{c_S^0\omega }}$ [17,19]. The $k-kl$ Mellor-Yamada level 2.5 model (MY) [13,21,22] and its modifications [16,23] are most used for parameterization of turbulent viscosity and diffusivity in OGCMs [13,15]. Usage of the $k-ϵ$ and $k-\omega$ models for the parameterization of $K_U$, $K_T$ and $K_S$ in OGCMs has some specific features [9,24].

The new closure for definition of $K_U$, $K_T$ and $K_S$ on the basis of Reynolds stress equations is proposed in [25]. The coefficients are defined for developed turbulence, low turbulence in the pycnocline, and in the bottom turbulence zone. In this parameterization [25], the two-equation turbulence model $k-ϵ$ is used to compute dissipation in the developed turbulence zone as well.

In our research, we couple an OGCM with the $k-\omega$ turbulence model [9,11]. Our choice is mainly justified by the mathematical advantage of the latter. Novelty consists in applying the splitting technique for solving the $k-\omega$ equations [5,26]. We split $k-\omega$ equations into the two main stages describing transport-diffusion and generation-dissipation processes. At the generation-dissipation stage, the equation for $\omega$ does not depend on k. It allows us to solve (at this stage) the $k-\omega$ equations analytically that ensures high computational efficiency.

The splitting algorithms with respect to physical processes and spatial coordinates are applied in the both turbulence model and INMOM. It allowed to make the approach to the mixing parameterization with using the $k-\omega$ model more convenient. The aims of the paper are as follows:

• to present a new splitting algorithm for solving $k-\omega$ turbulence equations that allows to reduce the complete system to the stages of transport-diffusion and generation-dissipation;
• to find an analytical solution of $k-\omega$ equations at the generation-dissipation splitting stage that is impossible for $k-kl$ and $k-ϵ$ closures;
• to demonstrate possibilities of controlling the obtained analytical solution of $k-\omega$ equations through its coefficients by means of different physical factors. We take into account such physical factors as climatic annual mean buoyancy frequency (AMBF) and Prandtl number as function of the Richardson number to goal this aim;
• to compare results of numerical experiments with the INMOM + $(k-\omega )$ to data on the hydrographic structure of the North Atlantic and Arctic Ocean to demonstrate physical effects of the accounting AMBF and variations of the Prandtl number.

2. Model and Methods

2.1. Equations of the Ocean General Circulation Model

The INMOM is based on a primitive equation system written in a generalized orthogonal coordinates on sphere under incompressibility, hydrostatics, and Boussinesq approximations [5,27]. The vertical coordinate is dimensionless variable $\sigma \in [0,1]$:

$$\sigma =\frac{z-\zeta }{H-\zeta },$$

(1)

where z is the usual vertical coordinate directed downward, H is the ocean depth, and $\zeta$ is the sea surface height, which is positive when under the undisturbed surface. Model equations written in symmetrized form are [5]:

$$\begin{array}{c}{\displaystyle D_{t}u-Z_{\sigma }\left(f+\eta \right)v=-\frac{Z_{\sigma }}{{\rho }_0{r}_x}\left[\frac{\partial }{\partial x}\left(p-\frac{g}{2}Z\rho \right)-\frac{g}{2}\left(\rho \frac{\partial Z}{\partial x}-Z\frac{\partial \rho }{\partial x}\right)\right]+\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }{K}_U\frac{\partial u}{\partial \sigma }+D_{u}u,}\hfill \end{array}$$

(2)

$$\begin{array}{c}{\displaystyle D_{t}v+Z_{\sigma }\left(f+\eta \right)u=-\frac{Z_{\sigma }}{{\rho }_0{r}_y}\left[\frac{\partial }{\partial y}\left(p-\frac{g}{2}Z\rho \right)-\frac{g}{2}\left(\rho \frac{\partial Z}{\partial y}-Z\frac{\partial \rho }{\partial y}\right)\right]+\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }{K}_U\frac{\partial v}{\partial \sigma }+D_{v}v,}\hfill \end{array}$$

(3)

$$\frac{\partial }{\partial \sigma }\left(p-\frac{g}{2}Z\rho \right)=\frac{g}{2}\left(\rho \frac{\partial Z}{\partial \sigma }-Z\frac{\partial \rho }{\partial \sigma }\right),$$

(4)

$$-\frac{\partial \zeta }{\partial t}+\frac{1}{r_x{r}_y}\left[\frac{\partial }{\partial x}\left(Z_{\sigma }{r}_{y}u\right)+\frac{\partial }{\partial y}\left(Z_{\sigma }{r}_{x}v\right)\right]+\frac{\partial w_{\sigma }}{\partial \sigma }=0,$$

(5)

$$D_{t}T=\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }{K}_T\frac{\partial T}{\partial \sigma }+D_{T}T-\frac{\partial R}{\partial \sigma },$$

(6)

$$D_{t}S=\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }{K}_S\frac{\partial S}{\partial \sigma }+D_{S}S,$$

(7)

$$\rho =\rho \left(T,S,\widehat{p}\right)\equiv \tilde{\rho }\left(T+\overline{T},S+\overline{S},\widehat{p}\right)-\tilde{\rho }\left(\overline{T},\overline{S},\widehat{p}\right).$$

(8)

Here, x and y are generalized orthogonal coordinates in horizontal subspace; $r_x$ and $r_y$ are metric coefficients, in spherical coordinates computes as: $r_x=R_{E}cosy$, $r_y=R_E$, $R_E$ is the Earth radius; $Z=(H-\zeta )\sigma +\zeta$, $Z_{\sigma }=H-\zeta$, ${\displaystyle Z_x\equiv \frac{\partial Z}{\partial x}}$, ${\displaystyle Z_y\equiv \frac{\partial Z}{\partial y}}$:

$$f=2\tilde{\mathsf{\Omega }}sin\phi ,\eta =\frac{1}{r_x{r}_y}\left(\frac{\partial r_y}{\partial x}v-\frac{\partial r_x}{\partial y}u\right),$$

(9)

l is the Coriolis parameter, $\phi$ is geographical latitude and $\tilde{\mathsf{\Omega }}$ is the Earth angular velocity. $D_t$ is the transport operator written in symmetrized form:

$$\begin{array}{c}{\displaystyle D_t\psi \equiv D_t\left(\overrightarrow{u}\right)\psi =\frac{1}{2}\left(Z_{\sigma }\frac{\partial \psi }{\partial t}+\frac{\partial Z_{\sigma }\psi }{\partial t}\right)+\frac{1}{2r_x{r}_y}\left[\frac{\partial }{\partial x}\left(Z_{\sigma }{r}_{y}u\psi \right)+Z_{\sigma }{r}_{y}u\frac{\partial \psi }{\partial x}\right.}\hfill \\ {\displaystyle \left.+\frac{\partial }{\partial y}\left(Z_{\sigma }{r}_{x}v\psi \right)+Z_{\sigma }{r}_{x}v\frac{\partial \psi }{\partial y}\right]+\frac{1}{2}\left[\frac{\partial }{\partial \sigma }\left(w_{\sigma }\psi \right)+w_{\sigma }\frac{\partial \psi }{\partial \sigma }\right],}\hfill \end{array}$$

(10)

where $\psi$ is transported variable, $\overrightarrow{u}=\left(u,v,w_{\sigma }\right)$ is the velocity vector in the $\sigma$-coordinate system; $w_{\sigma }$ is the vertical velocity in the $\sigma$-coordinate system defined as

$$w_{\sigma }=w-\left[\left(1-\sigma \right)\frac{\partial \zeta }{\partial t}+\frac{u}{r_x}\frac{\partial Z}{\partial x}+\frac{v}{r_y}\frac{\partial Z}{\partial y}\right],$$

(11)

T and S are deviations of potential temperature and salinity from their means $\overline{T}$, $\overline{S}$; $\widehat{p}={\rho }_{0}gZ$ is an approximate function of hydrostatic pressure; R is the penetrative solar radiation flux; $\rho$ is the density computed according to [28]; and $K_U$, $K_T$ and $K_S$ are the coefficients of vertical turbulent viscosity and diffusivity.

The operators $D_T$, $D_S$ describing lateral mixing of heat and salt along isopycnal surfaces are $D_T=D_S\equiv D_{\varphi }$:

$$\begin{array}{c}{\displaystyle D_{\psi }\psi =\frac{1}{r_x{r}_y}\frac{\partial }{\partial x}\left[\mu Z_{\sigma }\frac{r_y}{r_x}\left(\frac{\partial \psi }{\partial x}-\frac{{\rho }_x}{{\rho }_{\sigma }}\frac{\partial \psi }{\partial \sigma }\right)\right]-\frac{1}{r_x{r}_y}\frac{\partial }{\partial \sigma }\left[\mu Z_{\sigma }\frac{r_y}{r_x}\frac{{\rho }_x}{{\rho }_{\sigma }}\left(\frac{\partial \psi }{\partial x}-\frac{{\rho }_x}{{\rho }_{\sigma }}\frac{\partial \psi }{\partial \sigma }\right)\right]}\hfill \\ {\displaystyle +\frac{1}{r_x{r}_y}\frac{\partial }{\partial y}\left[\mu Z_{\sigma }\frac{r_x}{r_y}\left(\frac{\partial \psi }{\partial y}-\frac{{\rho }_y}{{\rho }_{\sigma }}\frac{\partial \psi }{\partial \sigma }\right)\right]-\frac{1}{r_x{r}_y}\frac{\partial }{\partial \sigma }\left[\mu Z_{\sigma }\frac{r_x}{r_y}\frac{{\rho }_y}{{\rho }_{\sigma }}\left(\frac{\partial \psi }{\partial y}-\frac{{\rho }_y}{{\rho }_{\sigma }}\frac{\partial \psi }{\partial \sigma }\right)\right],}\hfill \end{array}$$

(12)

where $\psi$ is diffused variable, $\mu$ is lateral diffusion coefficient, ${\displaystyle {\rho }_x\equiv \frac{\partial {\rho }_{pot}}{\partial x}}$, ${\displaystyle {\rho }_y\equiv \frac{\partial {\rho }_{pot}}{\partial y}}$, ${\displaystyle {\rho }_{\sigma }\equiv \frac{\partial {\rho }_{pot}}{\partial \sigma }}$. Turbulent viscosity operators $D_u$, $D_v$ are a combination of laplacian and biharmonic operators in plain form performing mixing along $\sigma$-surfaces [5].

Model includes interactive module which describes sea ice dynamics and thermodynamics [29].

2.2. The Two-Equation K-Omega Turbulence Model

The OGCM was subsequently coupled with the two-equation $k-\omega$ turbulence model. This model involves the solution of evolutionary equations for the turbulent kinetic energy k (TKE) and the turbulence dissipation frequency $\omega$ (TDF) [15,19]. The coefficients of turbulent viscosity $K_U$ and turbulent diffusivities for heat $K_T$ and salt $K_S$ are given as functions of TKE and TDF:

$$K_U=\frac{C_S^U}{c_S^0}\frac{k}{\omega }.$$

(13)

Here $C_S^U$ is the dimensionless stability function, $c_S^0=0.5562$ is its value for neutral stratification [15]; $k=0.5\overline{\left(u^{\prime }{u}^{\prime }+v^{\prime }{v}^{\prime }+w^{\prime }{w}^{\prime }\right)}$ is TKE, the bar above denotes averaging over the ensemble (or time), and the values with strokes are turbulent pulsations relative to the mean flow $\overrightarrow{u}=\left(u,v,w\right)$; ${\displaystyle \omega =\frac{ϵ}{{\left(c_S^0\right)}^{4}k}}$ is TDF, $ϵ$ is rate of TKE transformation to internal thermal energy of water due to dissipation.

For the temperature turbulent diffusivity, ${\displaystyle K_T=\frac{C_S^T}{c_S^0}\frac{k}{\omega }}$. Suppose that $C_S^U=c_S^0$ and ${\displaystyle C_S^T=\frac{C_S^U}{Pr}}$ in the elementary $\sigma$-layer, where Pr is Prandtl number (see, e.g., equations (2.73)–(2.74) from [30]). As well, suppose that $K_S=K_T$ for fully developed turbulent flows. As a result, we use the following dependencies: ${\displaystyle K_S=K_T=\frac{K_U}{Pr}}$.

The standard two-equation $k-\omega$ turbulence model is [15,19]:

$$\begin{array}{c}{\displaystyle D_{t}k=\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }\left(\frac{K_U}{{\sigma }_k}\frac{\partial k}{\partial \sigma }\right)+D_{k}k+Z_{\sigma }\left(K_U{G}^2-K_{\rho }{N}^2-{\left(c_S^0\right)}^4\omega k\right),}\hfill \end{array}$$

(14)

$$\begin{array}{c}{\displaystyle D_t\omega =\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }\left(\frac{K_U}{{\sigma }_{\omega }}\frac{\partial \omega }{\partial \sigma }\right)+D_{\omega }\omega +Z_{\sigma }\frac{\omega }{k}\left(c_1^{\omega }{K}_U{G}^2-c_3^{\omega }{K}_{\rho }{N}^2-c_2^{\omega }{\left(c_S^0\right)}^{4}k\omega \right).}\hfill \end{array}$$

(15)

Here ${\displaystyle G^2={\left(\frac{1}{Z_{\sigma }}\frac{\partial u}{\partial \sigma }\right)}^2+{\left(\frac{1}{Z_{\sigma }}\frac{\partial v}{\partial \sigma }\right)}^2}$ and ${\displaystyle N^2=\frac{g}{{\rho }_0}\frac{1}{Z_{\sigma }}\frac{\partial \rho }{\partial \sigma }}$ are the velocity shift and buoyancy frequency in the second degree ($N^2>0$ corresponds to stable stratification), $D_k$ and $D_{\omega }$ are mixing operators computed by the Equation (12).

Note: relations (14)–(15) are used only in layers where the developed turbulence presents and $k>0.03$ cm${}^2/$s${}^2$ [30]. Otherwise, small background values of exchange coefficients (intermittent turbulence) $K_{Umin}=1.0$ and $K_{Tmin}=0.05$ cm${}^2$/s are used.

2.3. Boundary and Initial Conditions

Boundary conditions for (2)–(8), (14) and (15) are as follows. At the sea free surface $\sigma =0$:

$$\begin{array}{c}{\displaystyle -\frac{K_U}{Z_{\sigma }}\frac{\partial u}{\partial \sigma }=\frac{{\tau }_{ax}}{{\rho }_0},-\frac{K_U}{Z_{\sigma }}\frac{\partial v}{\partial \sigma }=\frac{{\tau }_{ay}}{{\rho }_0},w_{\sigma }=0,}\hfill \\ {\displaystyle -\frac{K_T}{Z_{\sigma }}\frac{\partial T}{\partial \sigma }+{\alpha }_T\left(T-T_{obs}\right)=q_T,}\hfill \\ {\displaystyle -\frac{K_S}{Z_{\sigma }}\frac{\partial S}{\partial \sigma }+{\alpha }_S\left(S-S_{obs}\right)=q_S,}\hfill \\ {\displaystyle -\frac{K_U}{Z_{\sigma }{\sigma }_k}\frac{\partial k}{\partial \sigma }=C_g{\left(u_{\ast }^S\right)}^3,-\frac{K_U}{Z_{\sigma }{\sigma }_{\omega }}\frac{\partial \omega }{\partial \sigma }=0,}\hfill \end{array}$$

(16)

where ${\tau }_{ax}$, ${\tau }_{ay}$ are the wind stress components, ${\alpha }_T$, ${\alpha }_S$ are the relaxation coefficients, $q_T$ and $q_S$ are the total heat and salt fluxes, $C_g$ is wind generation parameter [13,31,32], ${\displaystyle u_{\ast }^S={\left(\frac{\sqrt{{\tau }_{ax}^2+{\tau }_{ay}^2}}{{\rho }_0}\right)}^{1/2}}$ is friction velocity in the upper ocean layer.

At the sea bottom $\sigma =1$:

$$\begin{array}{c}{\displaystyle \frac{K_U}{Z_{\sigma }}\frac{\partial u}{\partial \sigma }=\frac{{\tau }_{bx}}{{\rho }_0},\frac{K_U}{Z_{\sigma }}\frac{\partial v}{\partial \sigma }=\frac{{\tau }_{by}}{{\rho }_0},w_{\sigma }=0,}\hfill \\ {\displaystyle \frac{K_T}{Z_{\sigma }}\frac{\partial T}{\partial \sigma }=0,\frac{K_S}{Z_{\sigma }}\frac{\partial S}{\partial \sigma }=0,}\hfill \\ {\displaystyle \frac{K_U}{Z_{\sigma }{\sigma }_k}\frac{\partial k}{\partial \sigma }=C_g{\left(u_{\ast }^B\right)}^3,\frac{K_U}{Z_{\sigma }{\sigma }_{\omega }}\frac{\partial \omega }{\partial \sigma }=0,}\hfill \end{array}$$

(17)

where ${\tau }_{bx}=-{\rho }_0{C}_D\sqrt{u^2+v^2+e_b^2}u$ and ${\tau }_{by}=-{\rho }_0{C}_D\sqrt{u^2+v^2+e_b^2}v$ are bottom stress components, $C_D=2.5\times {10}^{-3}$, $e_b=5$ cm/s are empirical constants, ${\displaystyle u_{\ast }^B={\left(\frac{\sqrt{{\tau }_{bx}^2+{\tau }_{by}^2}}{{\rho }_0}\right)}^{1/2}}$ is friction velocity in the bottom ocean layer.

Thus, at the top and bottom surfaces we apply the corresponding boundary conditions for momentum, heat and salt fluxes, turbulent kinetic energy, and dissipation frequency [7,9]. At the lateral boundary, normal velocity, normal derivative of tangent velocity, heat, salt, turbulent kinetic energy and dissipation frequency fluxes are assumed to be equal to zero. The corresponding initial conditions for u, v, $\zeta$, T, S, k and $\omega$ are specified at $t=0$.

2.4. Numerical Algorithm

The numerical solution of the OGCM + $k-\omega$ model is based on the splitting method applying for evolutionary equations [33]. We formulate the system of governing equations in evolutionary form and represent the operator of the differential problem as a sum of more simple nonnegative suboperators. Just for this purpose it is convenient to use $\sigma$-coordinate system [5,27]. It allows us (1) to transform 3D computational domain to 3D volume with unity depth; (2) to exclude time derivative from boundary condition at the sea surface (${\displaystyle w=\frac{d\zeta }{dt}}$ for z-coordinates); (3) to represent transport operator in convenient nonnegative form acting on $D(x,y)\times [0,1]$; (4) rewrite, finally, both OGCM and turbulence model operators in evolutionary form. The key points of the approach proposed are as follows. When using the splitting method the form of a differential problem is of great importance. The most convenient form of equations is their symmetrized form. By the symmetrized form we mean the form of equations, which satisfies the conditions:

• The symmetrized form gives the form of the adjoint operator, which is close to the original one.
• This form leads to the finite difference approximation retaining the main properties typical of original differential operators (symmetry, skew-symmetry, nonnegativeness).
• From the form naturally follows the splitting of the problem operator into the sum of simple nonnegative operators.

Formulation Equations (2)–(8), (14) and (15) in symmetrized form makes it possible to use the splitting algorithm with respect physical processes and spatial coordinates x, y and $\sigma$ [27,33]. The equations for u, v, T, S, $\zeta$ at each time interval $t^j<t<t^{j+1}$ are split with respect to physical processes into two macro-stages: transport-diffusion of u, v, T, S:

$$D_{t}u-Z_{\sigma }\eta v=\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }{K}_U\frac{\partial u}{\partial \sigma }+D_{u}u,$$

(18)

$$D_{t}v+Z_{\sigma }\eta u=\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }{K}_U\frac{\partial v}{\partial \sigma }+D_{v}v,$$

(19)

$$D_{t}T=\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }{K}_T\frac{\partial T}{\partial \sigma }+D_{T}T-\frac{\partial R}{\partial \sigma },$$

(20)

$$D_{t}S=\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }{K}_S\frac{\partial S}{\partial \sigma }+D_{S}S,$$

(21)

and the adaptation of velocity and density fields:

$$\begin{array}{c}{\displaystyle \frac{\partial u}{\partial t}-fv=-\frac{1}{{\rho }_0{r}_x}\left[\frac{\partial }{\partial x}\left(p-\frac{g}{2}Z\rho \right)-\frac{g}{2}\left(\rho \frac{\partial Z}{\partial x}-Z\frac{\partial \rho }{\partial x}\right)\right],}\hfill \end{array}$$

(22)

$$\begin{array}{c}{\displaystyle \frac{\partial v}{\partial t}+fu=-\frac{1}{{\rho }_0{r}_y}\left[\frac{\partial }{\partial y}\left(p-\frac{g}{2}Z\rho \right)-\frac{g}{2}\left(\rho \frac{\partial Z}{\partial y}-Z\frac{\partial \rho }{\partial y}\right)\right],}\hfill \end{array}$$

(23)

$$\frac{\partial }{\partial \sigma }\left(p-\frac{g}{2}Z\rho \right)=\frac{g}{2}\left(\rho \frac{\partial Z}{\partial \sigma }-Z\frac{\partial \rho }{\partial \sigma }\right),$$

(24)

$$-\frac{\partial \zeta }{\partial t}+\frac{1}{r_x{r}_y}\left[\frac{\partial }{\partial x}\left(Z_{\sigma }{r}_{y}u\right)+\frac{\partial }{\partial y}\left(Z_{\sigma }{r}_{x}v\right)\right]+\frac{\partial w_{\sigma }}{\partial \sigma }=0,$$

(25)

If needed, within the transport-diffusion macro-stage the equations can be secondary split with respect to separate coordinates x, y and $\sigma$. At the adaptation macro-stage, a representation

$$u=\overline{u}+u^{\prime },v=\overline{v}+v^{\prime },\overline{u}=\underset{0}{\overset{1}{\int }}ud\sigma ,\overline{v}=\underset{0}{\overset{1}{\int }}vd\sigma$$

(26)

is used and implicit time scheme for the treatment of the depth averaged velocities and sea surface height $\zeta$ (linear shallow water equations) is applied. With respect to spatial coordinates, staggered grid [34], known in meteorology as “C”-grid [35], is used. Note, that the numerical algorithm for the solution of (2)–(8) is described in more detail in [5,27].

To solve (14) and (15) on the time interval $t_j<t<t_{j+1}$ we assume that three velocity components in the left-hand side of (14) and (15) are known. Splitting of turbulence Equations (14) and (15) is performed in terms of physical processes into two stages. At the first splitting stage we solve three-dimensional transport-diffusion equations:

$$D_{t}k=\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }\left(\frac{K_U}{{\sigma }_k}\frac{\partial k}{\partial \sigma }\right)+D_{k}k,$$

(27)

$$D_t\omega =\frac{1}{Z_{\sigma }}\frac{\partial }{\partial \sigma }\left(\frac{K_U}{{\sigma }_{\omega }}\frac{\partial \omega }{\partial \sigma }\right)+D_{\omega }\omega ,$$

(28)

with initial conditions $k=k^j$, $\omega ={\omega }^j$ at $t=t^j$. Note that within the transport-diffusion stage the equations can be secondary split with respect to separate coordinates x, y and $\sigma$ [5,26,33]. As a result we have $k^{j+1}\equiv k\left(t^{j+1}\right)$, ${\omega }^{j+1}\equiv \omega \left(t^{j+1}\right)$.

At the second splitting stage, on the same time interval $t_j<t<t_{j+1}$, we solve equations describing the generation-dissipation process [9,10,11]:

$$\begin{array}{ccc}\hfill \frac{\partial \omega }{\partial t}& =& B-C{\omega }^2,\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill \frac{\partial k}{\partial t}& =& \left(\frac{A}{\omega }-D\omega \right)k.\hfill \end{array}$$

(30)

with initial conditions equal to the solution of the first splitting stage: $k=k^{j+1}$, $\omega ={\omega }^{j+1}$ at $t=t^j$.

Here we can write separately the coefficients A, B, C and D for (29) and (30), using both the stability function for diffusivity $C_S^{\rho }=C_S^T$, and the Prandtl number Pr:

$$\begin{array}{ccc}\hfill A& =& \frac{1}{c_S^0}\left(C_S^U{G}^2-C_S^{\rho }{N}^2\right)=\frac{C_S^U}{c_S^0}\left(G^2-\frac{1}{Pr}{N}^2\right),\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill B& =& \frac{1}{c_S^0}\left(c_1^{\omega }{C}_S^U{G}^2-c_3^{\omega }{C}_S^{\rho }{N}^2\right)=\frac{C_S^U}{c_S^0}\left(c_1^{\omega }{G}^2-\frac{c_3^{\omega }}{Pr}{N}^2\right),\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill C& =& c_2^{\omega }{\left(c_S^0\right)}^4,D={\left(c_S^0\right)}^4.\hfill \end{array}$$

(33)

Note that $B>0$, $C>0$ and $D>0$ while A can change the sign. On the interval $t^j<t<t^{j+1}$ during the calculating k and $\omega$, variables ${\displaystyle N^2}$ and ${\displaystyle G^2}$ do not vary with respect to time.

As an example, we choose the coefficients A and B as functions of the Prandtl number. The latter is a function of $N^2$ and $G^2$ (see further) rather than a function of time at the interval $t_j<t<t_{j+1}$. Stability function for momentum is chosen as a constant $C_S^U=c_S^0$ (see, e.g., Equation (2.73) from [30]). In this case, the system (29) and (30) has a clear analytical solution at the time step of the circulation model. The analytical solution is available for $k-\omega$ model because at the generation-dissipation stage equation for $\omega$ does not depend upon k. The analytical solution of (29) and (30) has the form:

$$\begin{array}{c}{\displaystyle \omega =-r_d(r_m+r_p{a}_r)/(r_m-r_p{a}_r),}\hfill \\ {\displaystyle k=k^0{\left(\frac{E_1}{E_2}\right)}^{\frac{A}{2B}}{\left(\frac{E_3}{E_4}\right)}^{\frac{D}{2C}},}\hfill \end{array}$$

(34)

where

$$\begin{array}{c}{\displaystyle r_d=\sqrt{B/C},r_m={\omega }^0-r_d,}\hfill \\ {\displaystyle r_p={\omega }^0+r_d,a_r=exp\left(2\sqrt{BC}·t\right),}\hfill \\ {\displaystyle E_1={(r_m+r_p{a}_r)}^2,E_2=4{\left({\omega }^0\right)}^2{a}_r,}\hfill \\ {\displaystyle E_3=4r_d^2{a}_r,E_4={(r_m-r_p{a}_r)}^2.}\hfill \end{array}$$

We refer to the algorithm using the analytic solution (34) as the “splitting turbulence algorithm” (STA).

It follows from (34) that $\omega \to \sqrt{B/C}$ for $t\to \infty$ [9]. The expression $\omega =\sqrt{B/C}$ is the asymptotic of the solution (34). Using it, we can also find the asymptotic behavior of the TKE. Then the asymptotic solution of (29) and (30) has the form:

$$\omega =\sqrt{B/C},k=k^{0}exp\left(\gamma t\right),$$

(35)

where ${\displaystyle \gamma =\frac{A}{\omega }-D\omega =\frac{AC-BD}{\sqrt{BC}}}$.

3. Scenarios of Numerical Experiments

A set of experiments EX1–EX4 is performed as summarized in

Table 1. Model parameters corresponding to the numerical experiments EX1–EX4.

Experiments	Prandtl Number	Computation of k and $\mathit{\omega }$	Accounting AMBF	Viscosity and Diffusivity
EX1	$\begin{array}{c}Pr=\left\{\begin{array}{c}1,Ri\le 0.2\hfill \\ 5Ri,0.2<Ri<2\hfill \\ 10,Ri\ge 2\hfill \end{array}\right.,\\ {\displaystyle Ri=\frac{N^2}{G^2}}\end{array}$	$\begin{array}{c}A=G^2-N^2,\hfill \\ B=c_1^{\omega }{G}^2-c_3^{\omega }{N}^2\hfill \end{array}$	$\begin{array}{c}{\alpha }_M=1,\hfill \\ {\alpha }_{Cl}=0\hfill \end{array}$	$\begin{array}{c}{\displaystyle K_U=\frac{k}{\omega },}\hfill \\ {\displaystyle K_T=K_S=\frac{K_U}{Pr}}\hfill \end{array}$
EX2	$Pr=\left\{\begin{array}{c}1,Ri\le 0.2\hfill \\ 5Ri,0.2<Ri<2\hfill \\ 10,Ri\ge 2\hfill \end{array}\right.$	$\begin{array}{c}A=G^2-N_s^2,\hfill \\ B=c_1^{\omega }{G}^2-c_3^{\omega }{N}_s^2,\hfill \\ N_s^2={\alpha }_M{N}_M^2+{\alpha }_{Cl}{N}_{Cl}^2\hfill \end{array}$	$\begin{array}{c}{\alpha }_M=0.1,\hfill \\ {\alpha }_{Cl}=0.9\hfill \end{array}$	$\begin{array}{c}{\displaystyle K_U=\frac{k}{\omega },}\hfill \\ {\displaystyle K_T=K_S=\frac{K_U}{Pr}}\hfill \end{array}$
EX3	$\begin{array}{c}Pr=\left\{\begin{array}{c}{Pr}_0,Ri\le 0.2\hfill \\ a_{P}R{i}^2+b_{P}Ri+\hfill \\ +c_P,0.2<Ri<2\hfill \\ 10,Ri\ge 2\hfill \end{array}\right.\\ \begin{array}{c}{a}_P=1.3227,b_P=2.2487,\\ c_P=0.2116,{Pr}_0=0.7143\end{array}\end{array}$	$\begin{array}{c}{\displaystyle A=G^2-\frac{1}{Pr}{N}^2,}\hfill \\ {\displaystyle B=c_1^{\omega }{G}^2-\frac{c_3^{\omega }}{Pr}{N}^2}\hfill \end{array}$	$\begin{array}{c}{\alpha }_M=1,\hfill \\ {\alpha }_{Cl}=0\hfill \end{array}$	$\begin{array}{c}{\displaystyle K_U=\frac{k}{\omega },}\hfill \\ {\displaystyle K_T=K_S=\frac{K_U}{Pr}}\hfill \end{array}$
EX4	No	No	No	$\begin{array}{c}{\displaystyle K_U=\frac{K_0}{{(1+5Ri)}^2}+K_b,}\hfill \\ {\displaystyle K_T=K_S=}\hfill \\ {\displaystyle =\frac{K_U}{1+5Ri}+K_{Tb},}\hfill \\ {\displaystyle K_0=100c{m}^2/s,}\hfill \\ {\displaystyle K_b=1c{m}^2/s,}\hfill \\ {\displaystyle K_{Tb}=0.05c{m}^2/s}\hfill \end{array}$

. The model equations are written in a spherical coordinate system. The poles are located on the geographical equator at 120${}^{\circ }$ W and 60${}^{\circ }$ E. The simulation domain includes the Atlantic Ocean north of 30${}^{\circ }$ S, Arctic Ocean and Bering Sea. Open boundaries of the area pass through 30${}^{\circ }$ S and the straits of the Aleutian Islands. As well, the domain includes the Mediterranean, Black and Baltic Seas. The grid spacing in latitude and longitude is 0.25${}^{\circ }$. 40 $\sigma$-levels are set along the vertical with refinement near the ocean surface. The bottom topography is taken from [36]. The data were smoothed in accordance with the horizontal resolution of the model so that there are no strong bottom gradients. The model is forced at the sea surface with the CORE historical atmospheric datasets for the period 1948–2009 [37].

The initial conditions are climatic average January temperature and salinity fields taken from [38], no motion, and no sea ice. All the experiments were performed for the period 1948–2009. The OGCM equations are solved with a time step ${\tau }_{ocm}=1$ h. At the transport-diffusion splitting stage the $k-\omega$ equations are solved with the same time step ${\tau }_{ocm}=1$ h. At the generation-dissipation stage, in experiments EX1–EX3, an analytical solution is used with time step ${\tau }_T=5$ min for secondary splitting with using fractional step method. Note that coefficients A and B in the analytical solution depend on Prandtl number and AMBF.

Basic experiment EX1. In the EX1, in coefficients A and B we put $Pr=1$, and for turbulent exchange coefficients, Pr is calculated according to [39] in the form:

$$Pr=\left\{\begin{array}{c}1,Ri\le 0.2\hfill \\ 5Ri,0.2<Ri<2\hfill \\ 10,Ri\ge 2\hfill \end{array}\right..$$

(36)

Here, the functions $N^2$, $G^2$, ${\displaystyle Ri=\frac{N^2}{G^2}}$, Pr are constant at the interval $t_j<t<t_{j+1}$. The parameterization (36) proposed in [39] is currently used to define turbulent diffusivity in the OGCM (see, e.g., [23]). As well, the articles appeared where Pr is computed as a function of buoyancy frequency and turbulent characteristics [12,13,22].

Experiments EX2. The climatic annual mean buoyancy frequency $N_{Cl}$ (AMBF) is used in coefficients A and B for the EX2. AMBF is reliably defined by the World Ocean Atlas data using the climatic annual mean potential density for the whole period of observations: ${\rho }_{Cl}=\rho (T_{Cl},S_{Cl})$ [40], where $T_{Cl}$ is the potential temperature [41] and $S_{Cl}$ is the salinity [42]. The simple argument is used in the favour of taking AMBF into account. Annual mean vertical density gradients are lower than instantaneous ones in the warm season and are higher in the cold season. This favors to the mixing balance in the ocean within the annual cycle. In coefficients A and B, instead of $N^2$, we used:

$$N_S^2={\alpha }_M\frac{1}{Z_{\sigma }}\frac{g}{{\rho }_0}\frac{\partial {\rho }_M}{\partial \sigma }+{\alpha }_{Cl}{N}_{Cl}^2,$$

(37)

where ${\displaystyle N_{Cl}^2=\frac{g}{{\rho }_0}\frac{1}{Z_{\sigma }}\frac{\partial {\rho }_{Cl}}{\partial \sigma }}$, ${\rho }_M$ is the potential density calculated in the OGCM at the time $t_{j+1}$, and ${\alpha }_M$ and ${\alpha }_{Cl}$ are the weights of the contributions of the OGCM buoyancy frequency in the second degree and the AMBF in the second degree to the total $N_S^2$, ${\alpha }_M+{\alpha }_{Cl}=1$. In the first experiment with ${\alpha }_{Cl}=0.1$, the stability of the model solution to the accounting of the AMBF was verified. Further, EX2 results are demonstrated and analyzed for ${\alpha }_{Cl}=0.9$, which allowed us to obtain the most adequate results.

Experiment EX3. The aim of the EX3 is to study the sensitivity of the model solution to the Prandtl number variation. For neutral or unstable density stratification, we use either ${Pr}_0=0.7143$ [43] or ${Pr}_0=1.0$ [39]. For stable stratification, the Prandtl number is usually given as a function of the Richardson number $Ri$ [39,43]. It is assumed that the Prandtl number depends on the Richardson number:

$$Pr=\left\{\begin{array}{c}{\mathrm{Pr}}_0,Ri\le 0.2\hfill \\ a_{P}R{i}^2+b_{P}Ri+c_P,0.2<Ri<2\hfill \\ 10,Ri\ge 2\hfill \end{array}\right.,$$

(38)

where $a_P$, $b_P$ and $c_P$ are constant. The two variants were considered:

• $a_P=0.6790$, $b_P=3.5062$, $c_P=0.2716$ and ${Pr}_0=1$;
• $a_P=1.3227$, $b_P=2.2487$, $c_P=0.2116$ and ${Pr}_0=0.7143$.

Both variants are chosen to reduce Pr relative to (36). We do it to increase mixing in the tropical regions with respect to the EX1 results. The experiments showed that both variants yield close results. Further, we show the results of the second variant only.

Experiment EX4. In the EX4, the vertical coefficients of turbulent viscosity and diffusivity were calculated as functions of Richardson number according to [8] (see the Table 1).

4. Discussion

The experiments showed that computational time of the INMOM with the STA did not noticeably differ from the one in the EX4. It demonstrates the computational efficiency of the proposed STA with analytical solution. At the same time, the quality of reproducing ocean characteristics with using the STA is noticeably higher than in the EX4.

Based on the results of the numerical simulations, consider the two main questions concerned with (a) model adequacy with using the STA and (b) sensitivity of simulated characteristics to the change of turbulent mixing parameterization EX1–EX4. Analysis of temperature $T\left(t\right)$ and salinity $S\left(t\right)$ time series for EX1–EX4 from Equator to Arctic does not show any significant trends after the period of 15 years, when the model solution becomes quasi-stationary in the layer 0–500 m. Therefore, we consider just the period after stabilization of solutions when comparing the experimental results with the observations.

4.1. Comparison to the Generalized Observational Data

We use the data of the World Ocean Atlas for temperature $T_{Cl}$ [41] and salinity $S_{Cl}$ [42]. Consider the spatial distribution of the model T and S deviation from climatic data $T_{Cl}$ and $S_{Cl}$ and quality of reproducing vertical structure of the ocean upper layer.

Deviation of simulated temperature and salinity from the climatic data in the ocean upper layer. Denote the differences “model minus climate” for T and S as $dT$ and $dS$ respectively for the period 1963–2009. $dT$ are in the range $\pm {(0.2-0.7) }^{\circ }$C in the predominant part of the domain for the EX1, EX2 and EX3 in 0–30 m layer (Figure 1). Temperature in the EX1 is more than a one degree lower than the climatic data in the Equatorial currents, Sargasso and Norwegian seas. Using the AMBF in the EX2 significantly reduces $dT$ in comparison to the EX1 (see Figure 1a,b). Negative $dT$ of more than 2 ${}^{\circ }$C disappear in the EX2 on the Equator and in the Tropics. Square of negative deviations more than 1 ${}^{\circ }$C are halved. As well, negative $dT$ decrease significantly in the Caribbean sea and African upwelling area. There is a significant decrease in positive $dT$ in the Gulf Stream, the North Atlantic and Labrador currents. Negative $dT$ in the EX3 at the Equator, in the Caribbean and Norwegian Seas, in the African upwelling region decrease even more (see Figure 1c). However, the positive $dT$ in the EX1 and EX2 change to the negative ones in the EX3 for the North Atlantic Current. It is due to the additional entrainment of cold water from the seasonal thermocline by decreasing the Prandtl number relative to the EX1 and EX2.

$dS$ are in the range of ±(0.02–0.1) PSU in the greater part of the domain for all the experiments in the 0–30 m layer (Figure 2). Simulated $S>S_{Cl}$ at the Equator, in the Guiana Current and the Gulf Stream, in the Canadian sector of the Arctic. $dS$ reaches 0.5 PSU in the EX3 in the region of the Subpolar Gyre of the North Atlantic (Figure 2c). The same disadvantages are in the EX4. These disadvantages of the EX4 and EX3 are substantially reduced in the EX1 and EX2. Module of $dS$ reduced for EX2 on the Equator, Tropics and Subtropics, North Atlantic, Nordic Seas (Figure 2). $T-S$ diagrams are reproduced most realistically for the EX1 and EX2. Deviations $dT$ and $dS$ are compensated in the density in all the experiments in different degree: for $dT>0$, $dS>0$ and for $dT<0$, $dS<0$. The best compensation is observed in the EX2.

The greatest positive $dS$ in 0.5–2.0 PSU and $dT$ up to 5–6 ${}^{\circ }$C occur in the frontal zone of the Gulf Stream for all experiments (Figure 1 and Figure 2). These high $dT$ and $dS$ can be related both to the relatively coarse spatial resolution of the model and to the complex dynamics in the zone of high T and S gradients. Data on the climatic variables $T_{Cl}$ and $S_{Cl}$ in this region may experience appreciable variations depending on the change in the amount of observations. The compensation of $dT$ and $dS$ in the density allows to reproduce realistic currents, front dynamics, meanders and eddies.

Vertical structure of the upper ocean layer. We choose the “climatic” period of 30 years 1980–2009 to compare simulated T and S profiles to the ones from the Ocean Atlas. The most observations were performed for this period. We compare the profiles of simulated T and S mean for thirty February and August for 1980–2009 with the climatic monthly mean $T_{Cl}$ and $S_{C1}$ for February and August [41,42]. February is period of the greatest loss of buoyancy by the ocean. August is month of the greatest inflow of buoyancy from the atmosphere to the ocean. We choose the regions with the most different mixing conditions.

Deviation of $T\left(z\right)$ for the EX2 from the climatic one decreases by 0.7–0.9 ${}^{\circ }$C in the 0–300 m layer relative to the other experiments in the Sargasso Sea in February (Figure 3a). Zone of developed turbulence increases for EX2 in August relative to the EX1, and seasonal thermocline is closer to the observations (Figure 3b). Negative 3 ${}^{\circ }$C deviation of $T\left(z\right)$ from $T_{Cl}\left(z\right)$ exists in the 0–50 m layer for EX4 in August. There is no such deviation for the EX2. Salinity deviations from climate in the EX4 and EX3 are noticeable, especially in the February (Figure 3c). The S profiles are the best reproduced in the EX1 and, especially, in EX2 (Figure 3c,d). EX2 better simulates T and S profiles relative to the other experiments both during August and February in the Sargasso Sea (Figure 3).

Significant deviations of profiles T and S from climatic values for the EX1 are eliminated for the EX2 in the central part of the North Atlantic Current (Figure 4). These deviations were not for the first 20 years of the EX1. Thermocline is very well reproduced in the August for EX2 (Figure 4b). Temperature profiles are better simulated in the EX2 relative to the other experiments in February and August (Figure 4a,b). The main drawback of the parameterization EX4 is high deviation of salinity profiles from the climatic data (Figure 4c,d). The salinity profiles are the best reproduced in the EX2 and EX3 (Figure 4c,d).

Thermal inversion is maintained by the stratification of salinity in the recirculation zone of the West Spitsbergen Current all year round (Figure 5). Vertical structure of T and S for the EX2 is in the maximum agreement with the climatic data relative to the other experiments over here. This is particularly true for the reproduction of a typical cold water interlayer in the 50–130 m layer (see Figure 5b). As well, the EX3 yields a good approximation to the climatic data. EX4 overestimates the salinity gradient in the winter halocline in the layer 50–100 m (Figure 5c) that leads to unreal great temperature inversion (Figure 5a). This drawback overcomes in the experiments using the splitting turbulence algorithm (STA).

In general, the use of STA in mixing parameterization for the OGCM significantly reduces the salinity deviation from climatic data compared to the EX4 (Figure 3c,d, Figure 4c,d and Figure 5c,d). The vertical distributions of T and S in the EX2 (using AMBF) are in the best agreement with climatic data.

4.2. Sensitivity of Ocean Model Characteristics to the Changes in Mixing Parametrization

The analysis of model sensitivity to the changes in mixing parametrization is accompanied by qualitative comparison of the upper mixed layer (UML) depth, currents and sea surface height (SSH) to the observations.

Thickness of the ocean upper mixed layer. Distribution of the maximum UML depth in the North Atlantic is shown in the Figure 6. This distribution is typical for all winters of the whole period 1948–2009. Maximum UML depth is defined as the level where water density differs from the sea surface density less than 0.15 kg/m${}^3$. UML depth is highly sensitive to the used parameterizations (Figure 6). UML depth reaches 2–3 km in the Labrador, Norwegian and Greenland seas for the EX1. Area of these regions are sharply reduced for the EX2. UML depth is nowhere greater than 3 km for the EX3. Distribution and values of the UML depth are most similar to their estimates from observations [44] for the EX2 and EX3. EX3 reproduces a deeper UML than the EX1 and EX2 in the frontal zone of the Gulf Stream (Figure 6). The reason for this is in the greater diffusivity of heat and salt due to the smaller Prandtl numbers in the zone of decaying turbulence at relatively small UML depths.

Ocean currents. Changes in the fields of temperature, salinity and density associated with the change in the mixing parameterization cover almost the entire upper half of the baroclinic layer (Figure 3, Figure 4 and Figure 5). This causes a high sensitivity of the circulation and SSH to the change of the mixing parameterization. The greatest changes in ocean currents happen by the transition from parameterization EX4 using a simple dependence of the turbulent coefficients on the Richardson number to EX1–EX3 where STA is used.

Vectors of the velocity difference between EX1–EX3 using STA and EX4 are collinear with the mean circulation vectors in the Gulf Stream and Subpolar Gyre of the North Atlantic for the 0–50 m layer during 1980–2009 (Figure 7). There is an increase in velocities for EX1-EX3 by 20 cm/s and more relative to EX4. Maximum velocity increase occurs in the Gulf Stream separation from the coast (Figure 7a). This brought the current velocities closer to real values. Differences in the velocity between the EX1, EX2 and EX3 can reach 5–10 cm/s in the 0–50 m layer in separate regions. The results show that using AMBF in the EX2 allows to reproduce thermohaline characteristics better, and it should be recommended to use AMBF for computation of circulation.

Sea surface height. SSH is also sensitive to the change of parameterizations. SSH is associated with the annual cycle of sea ice extent. Figure 8 shows the SSH simulated in the EX4 (a), EX1 (b), EX2 (c) and EX3 (d), for the 62nd year of integration in April 2009. Sea ice extent reached its maximum in this month for the Arctic Ocean and Nordic Seas (NSIDC data). Beaufort Gyre is here the most characteristic feature of the circulation and SSH field. SSH distribution is more similar to the reconstruction according to the climatic data [38] for the EX2 relative to the other experiments. Beaufort Gyre is better expressed in the SSH field for EX2. Unrealistic overestimation of SSH was obtained for the EX4 in the region of the Siberian shelf. The greatest difference in the maximum SSH between the EX2 and EX1 is 15 cm in the Beaufort Gyre. EX3 reproduces high SSH values on the periphery of the Beaufort Gyre due to the steric effect as the result of the greater entrainment of dense water from the pycnocline into the upper ocean layer (Figure 8d). So, SSH is reproduced most realistically in the EX2 with taking into account the AMBF in the equations of turbulence (Figure 8c).

Heat flux at the ocean surface. Heat fluxes at the upper boundary of the ocean are calculated using the model SST. Model SST depends on turbulent mixing in its upper layer. Thus, the forcing for the OGCM also depends on the used mixing parameterization. Typical example for the sum of the latent $LE$ and sensible $Q_T$ heat fluxes shows that the space distributions of the heat fluxes are close to each other in the EX1, EX2 and EX3 for the period of maximum heat loss by the ocean in February (Figure 9). However, quantitative differences of fluxes can reach 50–200 W/m${}^2$ in some regions. For example, sensitivity of the heat fluxes to the change of the mixing parameterization revealed in the subtropics. Low heat losses $LE+Q_T<50$ W/m${}^2$ for EX1 changes to greater ones for the EX2 (Figure 9a,b). This is due to the following process. Turbulence zone penetrates to deeper layers in the preceding warm period for the EX2. Greater accumulation of heat in the upper ocean layer occurs. More significant water-air temperature differences for the EX2 relative to the EX1 arises in the winter.

4.3. Turbulent Energy and Omega

The spatial distribution and temporal evolution of TKE and TDF, differing from one experiment to another, are generally similar in EX1, EX2 and EX3. In all the experiments, the period of turbulence temporal variability ($1/\omega$) varies from seconds in layers of high density gradients, to several minutes in well-mixed layers. Such quantities are close to the observations [45]. In the evolution of TKE and TDF, the annual cycle is well pronounced. Against this background, there is a very high synoptic variability of these characteristics over time. Figure 10 shows the TKE field for the 62nd year of integration (2009) in the 0–30 m layer (for August). Zones of the Northern Trade winds and Westerlies in the Northern Hemisphere are marked with high TKE values ∼50–100 cm${}^2$/s${}^2$ associated with high wind speeds. Low TKE values of 1–3 cm${}^2$/s${}^2$ are observed in Horse Latitudes and the Intertropical Convergence zone. South of the Equator, the TKE values of 100–150 cm${}^2$/s${}^2$ are due to the loss of heat by the ocean in the winter of the Southern Hemisphere.

4.4. Numerical Aspects, Data Assimilation

One of the rapidly developing areas of numerical modelling of the oceanic and marine circulation is the assimilation of observational data [46,47,48,49,50,51]. In the nearest future, we plan to develop a high-resolution model of the Arctic Ocean and include in it 4D VAR data assimilation procedure. A successful solution of this problem is connected with the choice of an economical numerical algorithm for solving the optimality system—a coupled system of forward and adjoint ocean dynamics equations [50]. Here we can use the developed algorithm to solve the $k-\omega$ turbulence equations. The algorithm is based on an effective splitting technique with respect to physical processes. The formulation of an individual splitting stage describing the generation-dissipation process, allows one to solve them analytically which determines the computational efficiency. Our experiments showed the computational efficiency of the algorithm for forward simulation of Arctic—North Atlantic general circulation. This algorithm, in our opinion, will also be effective in solving 4D VAR data assimilation problems.

5. Summary

New results in the problem of vertical mixing parameterization for the ocean general circulation model (OGCM) are shown at the example of the simulation of climatic characteristics for North Atlantic and Arctic Ocean for 1948–2009. Like OGCM equations, turbulence equations are solved using multicomponent splitting methods. This allowed us to rationally improve the mixing parametrization using the $k-\omega$ model. A new splitting algorithm is proposed for solving turbulence equations that allows to reduce the complete system to the stages of transport-diffusion and generation-dissipation. An analytical solution was found for the $k-\omega$ equations at the generation-dissipation splitting stage that is impossible for $k-kl$ and $k-ϵ$ closures. This opens new possibilities to control ocean characteristics of the OGCM’s numerical experiments by means of some physical factors through coefficients of the proposed analytical $k-\omega$ model solution. We demonstrate examples of possibilities for controlling the OGCM solution by such physical factors as climatic annual mean buoyancy frequency (AMBF) and Prandtl number as function of the Richardson number.

For this purpose, the OGCM INMOM coupled with a two equation $k-\omega$ turbulence model is used for simulation of hydrophysical fields in the North Atlantic and Arctic Ocean. The model has horizontal resolution 0.25${}^{\circ }$ with 40 $\sigma$-levels along the depth.

The numerical results show the model’s satisfactory performance in simulating large-scale ocean circulation and upper layer dynamics. The deviations between simulated and observed temperature and salinity in the upper mixed layer are about ±(0.2–0.7) ${}^{\circ }$C and ±(0.02–0.10) PSU respectively. The vertical structure of temperature and salinity and their gradients in the pycnocline are acceptably reproduced. The simulations show that taking into account of the AMBF data improves the reproduction of temperature and salinity in the upper ocean layer. During the warm-up period in the tropics, accounting AMBF data increases zones of developed turbulence and makes the seasonal thermocline more realistic. In the Subarctic, Greenland and Norwegian seas, the temperature inversion and cold intermediate layer are more realistically reproduced. The experiments demonstrate that taking into account AMBF improves reproduction of the large-scale thermohaline and dynamical characteristics. Prandtl number variations improve the upper mixed layer depth simulation especially for high latitudes.

Replacement of the simple relation between turbulent viscosity and diffusion and the Richardson number [8] by $k-\omega$ parameterization increases velocity of currents in the upper 50 m layer by 10–20 cm/s or more. This is noticeable for the Gulf Stream, Greenland, Labrador and North Atlantic currents. The Subpolar Gyre of the North Atlantic becomes more intense and realistic. The SSH is also sensitive to the choice of turbulent parameterization. The largest SSH difference in the center of the Beaufort Gyre between different runs is about 15 cm. In experiments with the analytical solution of the $k-\omega$ equations, SSH in the Beaufort Gyre is well reproduced, especially with taking into account the AMBF in the equations of turbulence. Experiments with the Prandtl number tuning show that the increase of pycnocline waters entrainment into the mixed layer due to the steric effect has a significant influence on the Arctic Ocean circulation.

The experiments also demonstrate the computational effectiveness of the proposed approach based on the splitting method with respect to physical processes. We split $k-\omega$ equations into the two main stages describing transport-diffusion and generation-dissipation processes. At the generation-dissipation stage, the equation for $\omega$ does not depend on k. It allows us to solve both turbulence equations analytically that ensures high computational efficiency.