OCEAN CIRCULATION MODELING WITH K-OMEGA AND K-EPSILON PARAMETERIZATIONS OF VERTICAL TURBULENT EXCHANGE Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

Original Russian Text O The Authors, 2019 published in MORSKOY GIDROFIZICHESKIY ZHURNAL, Vol. 35, Iss. 6 (2019)

Ocean Circulation Modeling with K-Omega and K-Epsilon Parameterizations of Vertical Turbulent Exchange

V. B. Zalesny1,2,*, S. N. Moshonkin 12, V. L. Perov3, A. V. Gusev12

, 4

lMarchuk Institute of Numerical Mathematics of RAS, Moscow, Russian Federation Marine Hydrophysical Institute of RAS, Sevastopol, Russian Federation 3Hydrometeorological Research Center of Russian Federation, Moscow, Russian Federation 4Shirshov Institute of Oceanology of RAS, Moscow, Russian Federation * vzalesny@yandex.ru

Purpose. The main goal of the work is to advance the ocean general circulation model by improving description of the processes of vertical turbulent exchange of heat, salt and momentum which significantly affect quality of reproducing the ocean circulation and thermohaline structure using the models based on a system of the ocean hydrothermodynamics primitive equations. Methods and Results. The main instrument of the research is the sigma model of the oceanic and marine circulation developed at the Marchuk Institute of Numerical Mathematics, RAS. In the incompressibility, hydrostatics and Boussinesq approximations, the system of equations is supplemented with the k - ro and k - e parameterizations of the vertical turbulent exchange, the equations for which are solved by the splitting method applied to the physical processes. The equations are split into the stages describing transport-diffusion of the turbulence characteristics and their generation-dissipation. At the generation-dissipation stage, the equations for turbulent characteristics are solved analytically. At that, the stability functions resulted from application of the spectral algorithm are used. To assess quality of two parameterizations of the vertical turbulent exchange, the North Atlantic-Arctic Ocean circulation is numerically simulated and the upper ocean layer characteristics are studied.

Conclusions. It is shown that the structure of the North Atlantic-Arctic Ocean large-scale fields is sensitive to choice of the vertical turbulence models. In particular, application of the k - e parameterization is accompanied by a noticeably higher rate of involvement of the seasonal pycnocline waters in the developed turbulence zone than that resulting from application of the k - ro parameterization.

Keywords: ocean circulation, k-omega parameterization, k-epsilon parameterization, splitting method, vertical turbulent exchange, North Atlantic.

Acknowledgements: The investigation is carried out in the Marchuk Institute of Numerical Mathematics of RAS, Marine hydrophysical Institute of RAS and the Hydrometeorological Research Center of Russian Federation at support of the Russian Science Foundation (grant No. 17-77-30001) and the Russian Foundation for Basic Research (grant No. 18-05-00177).

For citation: Zalesny, V.B., Moshonkin, S.N., Perov, V.L. and Gusev, A.V., 2019. Ocean Circulation Modeling with K-Omega and K-Epsilon Parameterizations of Vertical Turbulent Exchange. Physical Oceanography, [e-journal] 26(6), pp. 455-466. doi:10.22449/1573-160X-2019-6-455-466

DOI: 10.22449/1573-160X-2019-6-455-466

© 2019, V. B. Zalesny, S. N. Moshonkin , V. L. Perov, A. V. Gusev © 2019, Physical Oceanography

Introduction

Ocean general circulation models (OGCMs) are, as a rule, based on systems of primitive equations [1-4]. It is possible to distinguish several main differences of the OGCMs from the models of classical hydrodynamics, namely: the presence of rotation (Coriolis acceleration); simplification of the third equation of motion to a hydrostatic relationship; "artificial" vertical mixing upon the achievement of stable stratification; free upper surface - the sea level and simplification of the

ISSN 1573-160X PHYSICAL OCEANOGRAPHY VOL. 26 ISS. 6 (2019) 455

The content is available under Creative Commons Attribution-NonCommercial 4.0

kinematic boundary conditions; complex coastal boundary, the presence of islands, variable bottom topography. The models of marine and oceanic dynamics are described by complex non-classical systems, including evolutionary and diagnostic equations, and require the development of special numerical algorithms [1, 2, 5-8]. The features of the marine and oceanic dynamics problems include a small horizontal scale of modeled phenomena and a small amount of observational data.

Practical simulations reveal the fact that the increase of spatial resolution and consideration of observational data are the key questions for improving the modeling [8-10]. An important resource for increasing the OGCM adequacy is the presence of parametrizations of subgrid processes [11-14]. Models should be constantly enriched with physically justified parameterizations, since it is hardly possible to explicitly describe ocean dynamics in the near future. This is especially true of meso- and sub-mesoscale variability, the dynamics of coastal zones and water areas at high latitudes, where the Rossby radius is about 1-5 km. The parametrization of subgrid processes remains one of the most important tasks of modeling and forecasting the hydrophysical and meteorological fields [14-24].

A feature of the ocean general circulation model of Marchuk Institute of Numerical Mathematics (INM) of the Russian Academy of Sciences is the formulation of its equations in evolutionary form and a solution algorithm based on the splitting method with respect to physical processes [1, 2]. The development of the model - its physical enrichment, the inclusion of a new factor or parameterization is also carried out within the framework of the splitting method. A new model subsystem is predetermined in evolutionary form, and its operator is represented as the sum of suboperators of a simpler structure [11, 12]. k - ra and k - e vertical turbulence models also satisfy this condition.

The parameterization of vertical turbulent exchange of momentum, heat and salt is an important factor in the physical development of the model. The spatial-temporal scales of vertical turbulent exchange are much smaller than the ones of large-scale circulation. They are approximately 1-10 km horizontally, 1-10 m vertically and from seconds to several hours in time. For primitive models, this is a subgrid process to be parameterized.

Turbulent mixing is often described in the OGCMs by a second-order operator with coefficients of turbulent exchange of momentum vu, heat vT and salt vS [11-15]. In order to determine the exchange coefficients, the models based on two equations are used [17]. The first equation is written for the turbulent kinetic energy (TKE) k, the second - either for the turbulence scale l (k - kl model) or for the TKE specific dissipation rate e (k - e model) or for the frequency ra (k - ra model) . The characteristics l, e, ra are related algebraically. These models describe developed turbulence and separate layers with developed and underdeveloped intermittent turbulence. Mixing in the underdeveloped turbulence layers is performed with implying double diffusion, destruction of internal and tidal waves and other effects [18, 19].

We develop an effective algorithm for solving k - ra and k - e turbulence equations. Like the OGCM equations, the turbulence equations are solved using the splitting method with respect to physical processes. The equations are split into the two main stages: transport - diffusion and generation - dissipation of the corresponding functions. Moreover, at the second stage of splitting the turbulence 456 PHYSICAL OCEANOGRAPHY VOL. 26 ISS. 6 (2019)

equations are solved analytically. The paper describes the results of using these parametrizations in the modeling of large-scale fields of the North Atlantic and the Arctic Ocean.

Ocean general circulation model

The INM RAS OGCM is formulated in a sigma coordinate system with a free surface o = (Z - Z)/(H - Z), where Z, Z, H are a geopotential vertical coordinate, sea surface height (SSH) and ocean depth [11, 12, 25, 26]. The model is based on primitive equations written in a bipolar orthogonal coordinate system on a sphere. The poles are located at the geographical equator at 120° W and 60° E [24]. The prognostic variables of the model are horizontal current velocities, SSH, potential temperature and salinity. Sea ice is simulated according to [27].

The OGCM equation set is split into the two main subsystems: transport -diffusion of substances and adaptation of current and density fields. The equations are approximated by finite differences on the grid "C". More detailed description of the model is given in [12, 25]. In the OGCM simulations the time step rocmis 1 hour. Vertical turbulence processes are described within the framework of the "diffusion" approach. The exchange coefficients are calculated using k - e or k - ra model.

Turbulent exchange equations

The k - s model is described by a system of two equations for the turbulent kinetic energy k and its specific dissipation rate s [17, 19]. The equations of the k - s model in o-coordinate system (a is directed downward) have the following form [17]:

dk 1 d v dk .

------f---Ak = vG2-v N2 -s,

dt H2 do o! do " p

k (1)

ds 1 d vu ds r s 2 - s "I s

-------As = lc.vG -c3v N -c2s I—,

dt H2 do os do L 1 " 3 P 2 J k

where G, N are the velocity shear and buoyancy frequencies calculated in the OGCM:

G2 =| 1 — ] + \ — — | N2 — ; A is a horizontal diffusion operator; v is

^ H do J ^ H do J Hp0 do

a coefficient of vertical turbulent viscosity, vp is a coefficient of potential density

vertical turbulent exchange. In the developed turbulence layer where k > kmin = 0.03 cm2/s2, the temperature and salinity diffusion coefficients are also assumed vT = vp and v^ = vT [28]. Dimensionless turbulent Schmidt

numbers for TKE and dissipation rate are ok = 1, os = 1,3. The remaining

s s s [-0,4, N2 > 0,

parameters are equal to: ci =1,44, c2 =1,92, c =< 2 Here and

[ 1,0, N < 0.

further it is assumed that k ^ 0, s ^ 0, œ ^ 0.

The solution of turbulence equations depends on the buoyancy and velocity shear frequencies calculated by the OGCM. In its turn, the coefficients of vertical turbulent exchange used in the OGCM depend on the turbulence characteristics. The coefficients of vertical turbulent viscosity vu and diffusion vT in the OGCM are calculated using similarity relations. For k - e models we have:

v=( c0 )3 Vt =( c0 )3 cp £ • (2)

Here cs, cp = cs are dimensionless stability functions for a vector and scalar; c°° = 0,5544 is a value of stability function at a neutral stratification. In weak turbulence layers the background values vu = 1, vs =vT = 0,05 cm2/s are used. The k - ro model equations have the following form [17]:

dk 1 d V dk 2 Ar2 i 0\4 ,

-M = VuG -vpN2 -(c] ) rok,

dt H2 da aak da dro 1 d vu dro = dt H2 da a da

cVG2 -cfyN2 -cro2 (cs )4 kro

(3)

Dimensionless turbulent Schmidt numbers for TKE and dissipation rate are = 2, aa= 2. According to [17], the rest of parameters are

c? = 0,555, c? = 0,833, < =i 0,6 N2 > 0 .

^ 2 3 1 1, N2 < 0

When calculating the coefficients of vertical turbulent viscosity vu and diffusion vT , we use the following relations

cu k cp k

vu = \ ?, Vt = % ?, (4)

where c^ = 0,5562. Like in k - e model, the background values vu = 1, vT =vS = 0,05 cm2/s are used in weak turbulence layers.

Splitting of turbulent exchange equations

The systems of equations (1) and (3) are solved using the splitting method with respect to physical processes. We split the procedure for solving equations into the two stages describing the transport - diffusion and generation - dissipation processes. At each time step tJ < t < t]+1, simpler split subsystems are solved. The solution obtained at the current stage is used as an initial condition at the next stage.

The transport - diffusion and generation - dissipation processes have different characteristic times. Transport - diffusion is a slow three-dimensional evolution of k, e, ? fields similar to the one of the OGCM background fields. Generation -dissipation is a fast spatial process describing the dynamics of turbulent disturbances.

At the stage of three-dimensional transport - diffusion, writing down the equation for TKE in one line, we have:

dk dt

ds dt

dro _ dt

1 d v dk 1 d v ds

doo lb

1 d v dro H2 do o do

+ Ak,

+ As.

- Aro.

(5)

On the upper boundary of the ocean for o = 0 we set the conditions:

V 1 dk i ,\3

tHdo=-c*(u*) ,

V 1 ds os H do

=Q0

(6)

v _1 dro=

oro H do~^

where u* is a friction velocity in the water at the ocean surface; cg is a dimensionless parameter depending on wind and waves (wind generation parameter): cg»10 [21] or cg « 40 [22]; Qs0,are the surface fluxes of specific

dissipation rate and dissipation frequency. One of the methods for setting the "dissipation flux" is given at p. 22, 23 in [28]. In our experiments it was supposed that = 0 . At the ocean bottom for o = 1 the normal fluxes k, s, ro are equal to zero.

At the generation - dissipation stage, substituting (2) in (1) and (4) in (3), we have:

ak dt '

ds ~dt~

~dt~

dro "dt

(c0 )3

(c0 )3

G - N

v

Pr

k s s k

2 A

cG2 - c3s

N Pr

c

u (

c

G2 -

N

2A

s V u (

Pr

k s

I

ro

c

c

s V

2 N2 ^ c1G2 - c3ro — 1 3 Pr

1 -(c0 )4

ro

- c? (c0 )4

in v /

ro

ro.

(7)

(8)

(9)

(10)

Here the relations V = v = V = V / Pr, Pr = cu / cf are used, Pr is

the Prandtl number. It can be seen that unlike k - e , the system of k - ra equations (6) has a simpler form - the equation for ra does not depend on k. Since in the splitting scheme all the coefficients (7) - (10) are taken from the previous step, it is easy to write out the analytical solution (10).

An analytical solution for the k - ra model at the generation - dissipation stage is given in [15, 23 and 24]. In these works, it was noted that the selection of parameters and the use of observational data in the simulations are important factors in increasing the adequacy of modeling the large-scale ocean circulation. The possibility of increasing the OGCM's adequacy by taking into account the observational data on the climatic annual mean buoyancy frequency in A and B is also shown. From this point of view, the algorithm is quite flexible. Firstly, it allows us to vary A and B by changing the stability functions or the Prandtl number, while preserving the simplicity and stability of the computations. Secondly, it makes it possible to solve the turbulence equations with a large time step equal to the one for the OGCM. So, the time step rocm is 1 hour in numerical experiments for the OGCM. At the transport - diffusion stage the equations of the k - ra model are also solved with the time step 1 hour. At the generation -dissipation stage the k - ra equations are solved with the step tt varied within the range 1min <tt < 1h .

Analytic solution for the k - e model at generation - dissipation stage

We transform the k - e model equations at the generation - dissipation stage to a form similar to the one for the k - ra model. To do this, we introduce an auxiliary function

ra = e/k/(c0)3/< . (11)

The equations (7), (8) with regard to assumption that k ^ 0, e ^ 0 will be written down as follows:

d ln k

dt -G-N2/Pr)ra-(cs)3cura, (12)

^ = (ceG2 - c^N2 / Pr)1 - c2 (c0 )3 cU ra . (13)

dt v ' ra v '

Subtracting (12) from (13), we obtain the equations for k and ra, moreover, the equation for ra does not depend on k:

d ln k = A / ra - D®, (14)

= B / ra - C ra, (15)

dt

A = G1 - N2/ Pr, D = (c0 )3 cus ,

B = (ce -1) G2 - (c3 -1)N2 / Pr, C = (c0 )4 (c2 -1). (16)

460 PHYSICAL OCEANOGRAPHY VOL. 26 ISS. 6 (2019)

dt d ln râ

The analytical solution of (14), (15) has the form:

b ( cc0 /4bTc )+th (4BCt )

c (œ0 /4bTC) th (4BCt)+i, 1 + (4bTC/ cc0) th(VBCt)

co =

k = k0

1+( cc0/v b / c ) th(4bct )

(18)

CD0 = 8° / k° / D,

where C° and k° are the values at the initial time moment at the generation -dissipation stage. In terms of k, 8 we obtain:

k = k0

1+th(4BC t ) ( Wb / c k0 / s0 )

A/B

s =

1 + th(^/BCt)/(DyfBÏCk0 /s0)

( s0 / k0 ) + D>/BTC th (4ËCt )

c (s0/k°)th(4bcî)+d4bTC

D/C '

Dk.

Stability functions based on the spectral algorithm

Usually the form of the stability functions (when averaging the equations of geophysical hydrodynamics according to Reynolds) does not take into account scale differences and mix them [17, 28, 29]. Generally to say, the turbulence models used in the OGCMs should take into account the impact of multiscale processes. Their consideration can be introduced using the spectral algorithm [16, 3°]. In the spectral algorithm, nonlinear equations of ocean circulation are replaced by linear spectral equations with random external forcing known as the Langevin equations. After a series of transformations, the latter are reduced to the two differential equations determining new stability functions dependent only on the Richardson number Ri = N2 / G2 [16, 3°]:

1 + 8Ri2

c0 1 + 2,3Ri + 35Ri2

1,4 - 0,01Ri + 1,29Ri2

1 + 2,44Ri + 19Ri2

(19)

It should be noted that the spectral algorithm does not support the widespread notion of complete attenuation of turbulence by highly stable stratification. In this case, the first of the stability functions (19) gives the value 0.22 while the second one tends to zero asymptotically for large Ri.

c

c

c

Numerical experiments and simulation results

The main purpose of numerical experiments is to assess the influence of the two turbulence parametrizations k-e and k-ra on the structure of modeled hydrophysical fields. The simulations were carried out using the INM RAS OGCM with built-in k - e and k - ra subsystems and with new stability functions (19).

Statement of the numerical experiments. The simulations were carried out for the water area including the Atlantic Ocean northward of 30° S, the Arctic Ocean and the Bering Sea. The simulations were carried out from January 1, 1976 to December 31, 1977. The computational area has open boundaries located at 30°S and in the straits of the Aleutian Islands. The area includes the Mediterranean, Black and Baltic Seas. The grid step in latitude and longitude is 0.25°. 40 sigma levels are set in vertical with the refinement near the ocean surface. The ocean bottom topography was smoothed in accordance with the horizontal resolution of the model so that there are no strong bottom gradients. Model depth is limited to a minimum of 10 m.

Boundary conditions at the ocean surface are calculated using the atmospheric characteristics according to the CORE-II (Datasets for Common Ocean-ice Reference Experiments - Phase II) data [31]. The fluxes of sensible and latent heat, fresh water and wind stress are calculated with 1-hour discreteness according to CORE-II data on air temperature, humidity, wind speed components and sea level pressure using the model water temperature. The fluxes of long-wave and shortwave radiation (taking into account its penetrating ability) are set with 1 day discreteness [31]. The precipitation rate and river runoff are set with 1 month discreteness.

On solid coastal boundaries the conditions of no normal flow and no heat and salt flux are set. At the liquid boundaries from the surface to the bottom the climatic monthly mean values of temperature and salinity are set. The runoff of the main rivers is indirectly taken into account in the boundary condition for salinity.

January climatic fields of the ocean temperature and salinity, the absence of movement and sea ice [24] are taken as initial conditions. In the boundary condition (6) = 0.

The turbulence equations at the generation - dissipation stage are solved analytically with the time step tt = 1 h equal to the OGCM step toct .

The simulation results. We compare the simulation results with the observational data for "C" ocean weather station (OWS) with 52.75° N, 35.5° W coordinates, where the upper layer probing was carried out 8 times a day 1.

Fig. 1 presents the temperature profiles at the initial stage of free convection development in September - October 1977. Both the observational data at the OWS "C" and the average over the model cell (covering the OWS) OGCM calculations with k - e and k - ra mixing parametrizations are given. All profiles are constructed by discrete data using cubic spline interpolation [32].

1 A collection of climatological and statistical data at the ocean station "C" (52° 45' N, 35° 30' W) for 1976-1980 period. Section 1. Oceanographic and hydrochemical observations. Obninsk, 1984. 338 p.

(a) (b)

F i g . 1. Temperature profiles averaged for September, 1-10 (a) and October, 1-10, (b) 1977. Green line marks observational data from the oceanic weather station "C"1; black line marks the k - o> model results; red line marks the k- s model results

For the month under consideration the observed temperature of the upper quasihomogeneous layer (UQL) decreased by 1.75°C. The simulations indicate that the UQL deepening due to using the k - s model is noticeably larger in comparison with both the observational data and the k - < model. The UQL temperature is reproduced better when using the k - < parameterization.

Using the maximum curvature of the daily mean vertical temperature profiles we determine the position of the ocean UQL lower boundary h and relate it to we -the rate of thermocline water involvement into the UQL: we = dh / dt. Using this ratio we can assess the rate of thermocline water involvement into the UQL. The simulations show that during October 1977 the value is we ~ 0,5 m/day for the k - < model and observational data, while it is we ~ 0,8 m/day for the k - s model.

F i g. 2. UQL thickness in the model coordinates on October 10, 1977: a represents the k - o> model results; b represents the k- s model resutls. The water potential density in the UQL differs from the one at the ocean surface by less than 0.15 kg/m3

Fig. 2 presents the UQL thickness for October 10, 1977 in the North Atlantic, Norwegian and Greenland. All the qualitative features of the UQL thickness distribution coincide in both cases. However, we note that in most of the water area k - s the parametrization (compared with k - <) leads to stronger mixing. That is, in the initial period of free convection the UQL thickness, which is an important characteristic of vertical mixing, is sensitive to the choice of parameterization.

Conclusions

1. The application of the k - s and k - © parametrizations of turbulent mixing in the INM RAS ocean general circulation model is described. When calculating the parameters of turbulence models, the stability functions are used based on the spectral algorithm.

2. The method of splitting with respect to physical processes is used for the numerical solution of turbulence equations. A feature of the method is the exact solution of the split equations at the generation - dissipation stage. Previously this technique was used for k - © model [15, 23, 24]. In this case, as well, an analytical solution is written down for the k - s model and used in the simulations.

3. It is shown that the structure of large-scale fields of the North Atlantic -the Arctic Ocean is sensitive to the choice of vertical turbulence models. So, for example, in the k - s model the rate of seasonal pycnocline water involvement into the zone of developed turbulence is noticeably higher than when using the k - © model.

REFERENCES

1. Marchuk, G.I., Dymnikov, V.P. and Zalesny, V.B., 1987. [Mathematical Models in Geophysical Hydrodynamics and Methods of their Numerical Implementation]. Leningrad: Gidrometeoizdat, 295 p. (in Russian).

2. Sarkisyan, A.S., 2015. A View of Progress in the Numerical Modeling of Physical Characteristics of the World Ocean in the Light of Sixty-Year Experience. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 51(3), pp. 320-336. https://doi.org/10.1134/S000143381503010X

3. Demyshev, S.G., Korotaev, G.K., Knysh, V.V., 2004. Modeling the Seasonal Variability of the Temperature Regime of the Black Sea Active Layer. Izvestiya. Atmospheric and Oceanic Physics, 40(2), pp. 227-237.

4. Ivanov, V.A., Kubryakov, A.I., Mikhailova, E.N. and Shapiro, N.B., 1997. Modelling of Circulation in the Gulf of Izmir. Physical Oceanography, [e-journal] 8(1), pp. 47-55. https://doi.org/10.1007/BF02522565

5. Marchuk, G.I. and Paton, B.E., 2012. The Black Sea as a Simulation Ocean Model. Russian Journal of Numerical Analysis and Mathematical Modelling, 27(1), pp. 1-4. doi: 10.1515/rnam-2012-0001

6. Agoshkov, V.I., Parmuzin, E.I. and Shutyaev, V.P., 2013. Observational Data Assimilation in the Problem of Black Sea Circulation and Sensitivity Analysis of its Solution. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 49(6), pp. 592-602. https://doi.org/10.1134/S0001433813060029

7. Zalesny, V.B., Gusev, A.V. and Fomin, V.V., 2016. Numerical Model of Nonhydrostatic Ocean Dynamics Based on Methods of Artificial Compressibility and Multicomponent Splitting. Oceanology, [e-journal] 56(6), pp. 876-887. doi:10.1134/S0001437016050167

8. Zalesnyi, V.B., Gusev, A.V. and Agoshkov, V.I., 2016. Modeling Black Sea Circulation with High Resolution in the Coastal Zone. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 52(3), pp. 277-293. doi:10.1134/S0001433816030142

9. Ibrayev, R.A., Khabeev, R.N. and Ushakov, K.V., 2012. Eddy-Resolving 1/10° Model of the World Ocean. Izvestiya. Atmospheric and Oceanic Physics, 48(1), pp. 37-46. doi:10.1134/S0001433812010045

10. Mizyuk, A.I. and Korotaev, G.K., 2020. Black Sea Intrapycnocline Waves Resulted from Numerical Simulation of Basin Circulation. Izvestiya. Atmospheric and Oceanic Physics (in press).

11. Zalesny, V.B. and Gusev, A.V., 2009. Mathematical Model of the World Ocean Dynamics with Algorithms of Variational Assimilation of Temperature and Salinity Fields. Russian Journal of Numerical Analysis and Mathematical Modelling, 24(2), pp. 171-191. doi:10.1515/rjnamm.2009.012

12

13

14

15

16

17

18.

19

20

21

22

23

24

25

26

27

Zalesny, V.B., Marchuk, G.I., Agoshkov, V.I., Bagno, A.V., Gusev, A.V., Diansky, N.A., Moshonkin, S.N., Tamsalu, R.E. and Volodin, E.M., 2010. Numerical Simulation of Large-Scale Ocean Circulation Based on the Multicomponent Splitting Method. Russian Journal of Numerical Analysis and Mathematical Modelling, 25(6), pp. 581-609. doi:10.1515/rjnamm.2010.036

Volodin, E.M., Dianskii, N.A. and Gusev, A.V., 2010. Simulating Present-Day Climate with the INMCM4.0 Coupled Model of the Atmospheric and Oceanic General Circulations. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 46(4), pp. 414-431. doi:10.1134/S000143381004002X

Zalesny, V.B., Diansky, N.A., Fomin, V.V., Moshonkin, S.N. and Demyshev, S.G., 2012. Numerical Model of the Circulation of the Black Sea and the Sea of Azov. Russian Journal of Numerical Analysis and Mathematical Modelling, 27(1), pp. 95-112. doi:10.1515/rnam-2012-0006

Zalesny, V.B. and Moshonkin, S.N., 2019. Sensitivity of the Ocean Circulation Model to the k-ro Vertical Turbulence Parametrization. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 55(5), pp. 470-479. doi:10.1134/S0001433819050141 Perov, V.L., 2012. Calculation of a Coefficients of Turbulent Mixing on the Basis of the Spectral Algorithm and Its Implementation in the Model COSMO-RU. In: E. S. Nesterov, 2012. Proceedings of Hydrometcentre of Russia. Moscow: Hydrometcentre of Russia. Issue 347, pp. 74-85 (in Russian).

Warner, J.C., Sherwood, C.R., Arango, H.G. and Signell, R.P., 2005. Performance of Four Turbulence Closure Models Implemented Using a Generic Length Scale Method. Ocean Modelling, [e-journal 8(1-2), pp. 81-113. doi:10.1016/j.ocemod.2003.12.003 Large, W.G., McWilliams, J.C. and Doney, S.C., 1994. Oceanic Vertical Mixing: a Review and a Model with a Nonlocal Boundary Layer Parameterization. Reviews of Geophysics, [e-journal] 32(4), pp. 363-403. doi:10.1029/94RG01872

Canuto, V.M., Howard, A.M., Cheng, Y., Muller, C.J., Leboissetier, A., and Jayne, S.R., 2010. Ocean Turbulence, III: New GISS Vertical Mixing Scheme. Ocean Modelling, [e-journal] 34(3-4), pp. 70-91. doi:10.1016/j.ocemod.2010.04.006

Moshonkin, S.N., Zalesny, V.B. and Gusev, A.V., 2019. A Splitting Turbulence Algorithm for Mixing Parameterization in the Ocean Circulation Model. IOP Conference Series: Earth and Environmental Science, 231(1), 012038. doi:10.1088/1755-1315/231/1/012038

Zaslavskii, M.M., Zalesny, V.B., Kabatchenko, I.M. and Tamsalu, R., 2006. On the Self-Adjusted Description of the Atmospheric Boundary Layer, Wind Waves, and Sea Currents. Oceanology, 46(2), pp. 159-169. doi:10.1134/S0001437006020020

Noh, Y., Ok, H., Lee, E., Toyoda, T. and Hirose, N., 2016. Parameterization of Langmuir Circulation in the Ocean Mixed Layer Model Using LES and Its Application to the OGCM. Journal of Physical Oceanography, 46(1), pp. 57-78. doi:10.1175/JPO-D-14-0137.1

Moshonkin, S., Zalesny, V. and Gusev, A., 2018. Simulation of the Arctic-North Atlantic Ocean Circulation with a Two-Equation K-Omega Turbulence Parameterization. Journal of Marine Science and Engineering, [e-journal] 6(3), 95. doi:10.3390/jmse6030095

Moshonkin, S.N., Zalesny, V.B. and Gusev, A.V., 2018. Algorithm of the k-ro Turbulence Equations Solution for the Ocean General Circulation Model. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 54(5), pp. 495-506. doi:10.1134/S0001433818050079 Marchuk, G.I. and Zalesny, V.B., 2012. Modeling of the World Ocean Circulation with the Four-Dimensional Assimilation of Temperature and Salinity Fields. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 48(1), pp. 15-29. doi:10.1134/S0001433812010070 Zalesny, V.B., Agoshkov, V.I., Shutyaev, V.P., Le Dimet, F. and Ivchenko, B.O., 2016. Numerical Modeling of Ocean Hydrodynamics with Variational Assimilation of Observational Data. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 52(4), pp. 431-442. doi:10.1134/S0001433816040137

Yakovlev, N.G., 2009. Reproduction of the Large-Scale State of Water and Sea Ice in the Arctic Ocean in 1948-2002: Part I. Numerical Model. Izvestiya, Atmospheric and Oceanic Physics, [e-journal] 45(3), pp. 357-371. doi:10.1134/S0001433809030098 CAL OCEANOGRAPHY VOL. 26 ISS. 6 (2019) 465

28. Burchard, H., Bolding, K. and Villarreal, M.R, 1999. GOTM, a General Ocean Turbulence Model: Theory, Implementation and Test Cases. Space Applications Institute, 103 p.

29. Mellor, G.L. and Yamada, T., 1974. A Hierarchy of Turbulence Closure Models for Planetary Boundary Layers. Journal of the Atmospheric Sciences, [e-journal] 31(7), pp. 1791-1806. doi:10.1175/1520-0469(1974)031<1791:AHOTCM>2.0.CO;2

30. Sukoriansky, S., Galperin, B. and Perov, V., 2005. Application of a New Spectral Theory of Stably Stratified Turbulence to the Atmospheric Boundary Layer over Sea Ice. Boundary-Layer Meteorology, [e-journal] 117(2), pp. 231-257. doi:10.1007/s10546-004-6848-4

31. Large, W.G. and Yeager, S.G., 2009. The Global Climatology of an Interannually Varying Air-Sea Flux Data Set. Climate Dynamics, [e-journal] 33(2-3), pp. 341-364. doi:10.1007/s00382-008-0441 -3

32. Hyman, J.M., 1983. Accurate Monotonicity Preserving Cubic Interpolation. SIAM Journal on Scientific and Statistical Computing, [e-journal] 4(4), pp. 645-654. doi:10.1137/0904045

About the authors:

Vladimir B. Zalesny - Leading Research Associate, Marchuk Institute of Numerical Mathematics of RAS (8 Gubkina str., Moscow, 119333, Russian Federation), Dr. Sci. (Phys.-Math.), professor, ORCID ID: 0000-0003-3829-3374, ResearcherlD: E-4628-2014, SPIN-code: 5546-5227,^ Scopus Author ID: 6603178923, vzalesny@yandex.ru

| Sergey N. Moshonkin | - Leading Research Associate, Marchuk Institute of Numerical Mathematics of RAS (8 Gubkina str., Moscow, 119333, Russian Federation), Dr. Sci. (Phys.-Math.), ORCID ID: 0000-0001-7701-4521, Scopus Author ID: 6603090797, atarexm@himki.net

Veniamin L. Perov - Leading Research Associate, Hydrometeorological Research Center of the Russian Federation (11-13 Bolshoy Predtechensky lane, Moscow, 123242, Russian Federation), Ph. D. (Phys.-Math.), veniamin.perov@mail.ru

Anatoly V. Gusev - Senior Research Associate, Marchuk Institute of Numerical Mathematics of RAS (8, Gubkina str., Moscow, 119333, Russian Federation), Ph. D. (Phys.-Math.), ORCID ID: 0000-0002-6463-3179, ResearcherID: A-8528-2014, SPIN-code: 1387-3344, Scopus Author ID: 56828803500, anatoly.v.gusev@gmail.com

Contribution of the co-authors:

Vladimir B. Zalesny - general scientific supervision of the research, initiation of research, literature review on the research problem, formulation of the research goals and objectives, statement and formulation of the problem, development of a mathematical model, discussion of the work

result^_

I Sergey N. Moshonkin| - participation in the discussion of the paper materials, statement of the research problem, development of a methodology, computer implementation of algorithms, qualitative and quantitative analysis of the results, processing and description of the research results, preparation of the initial version of the paper text.

Veniamin L. Perov - study of the concept, advisory assistance, development of methodology, analysis of materials on the research theme, analysis of the applied research methodology, analysis and generalization of the research results.

Anatoly V. Gusev - selection and justification of numerical methods for solving equations, mathematical model correction and calculations, data collection and systematization, computer implementation of algorithms, development and debugging of a computer program for solving a problem, revision of the text, visualization/presentation of data in the text.

All the authors have read and approved the final manuscript.

The authors declare that they have no conflict of interest.