Non-linear Invariants of a Discrete System of the Sea Dynamics Equations in a Quasi-Static Approximation Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

ISSN 1573-160X, PHYSICAL OCEANOGRAPHY, Vol. 30, Iss. 5, pp. 523-548 (2023)

Original Russian Text © S. G. Demyshev, 2023, published inMorskoy Gidrofizicheskiy Zhurnal, 2023, Vol. 39, Iss. 5, pp. 557-583.

Original article

Nonlinear Invariants of a Discrete System of the Sea Dynamics Equations in a Quasi-Static Approximation

S. G. Demyshev 0

Purpose. The study is purposed at obtaining the approximations providing the presence of discrete non-linear invariants for the difference system of the sea dynamics equations in the absence of external forces, friction and diffusion, and at analyzing the features of the resulting schemes at the example of calculating the Black Sea circulation for 2011.

Methods and Results. The method of undetermined coefficients at which the new unknowns are introduced is applied, that makes it possible to satisfy the additional conditions. The schemes providing simultaneous preservation of temperature in the first and the K-th (K > 1) degrees and salinity in the first and the L-th (L > 1) degrees, were obtained. The approximations of temperature and salinity found on the box faces with a polynomial dependence of density on temperature and salinity lead to a divergent form of the density advection equation. This form provides fulfilling the law of conservation both of the total energy and the sum of kinetic and dynamic potential energy in a discrete formulation. Based on the analysis of circulation in the Black Sea in 2011, it is shown that at increase of the degree of invariants, the following effects take place: the gradients in the temperature field in the frontal zones as well as the processes of the saltier water upwelling in the sea center and the fresher water downwelling along its periphery are intensified, and the intensity of small-scale features in the vertical velocity field decreases.

Conclusions. A discrete dynamical model in a quasi-static approximation was obtained. It has a number of nonlinear invariants corresponding to the continuous problem. The results of calculating the Black Sea circulation for real conditions in 2011 showed that presence of the degree invariants exceeding two made it possible to specify the circulation features on small scales.

Keywords: numerical modeling, kinetic energy, potential energy, discrete conservation laws, Black Sea, cyclonic circulation, anticyclonic eddies, sea dynamics

Acknowledgments: The work was carried out with financial support of the Russian Science Foundation grant 23-27-00141.

For citation: Demyshev, S.G., 2023. Nonlinear Invariants of a Discrete System of the Sea Dynamics Equations in a Quasi-Static Approximation. Physical Oceanography, 30(5), pp. 523-548.

© S. G. Demyshev, 2023

© Physical Oceanography, 2023

1. Introduction

Currently, the requirements for the accuracy of solving discrete equations of the dynamics of oceans and seas, which ensures the reliability of the marine environment condition forecast, are increasing. One of the aspects of work to improve numerical models of sea dynamics is due to the need to fulfill conservation laws for a finite-difference system of equations. This approach is ideologically based on the well-known Noether's theorem [1], establishing one-to-one correspondence between conservation laws and the properties of its solution for a hyperbolic system of hydrodynamic equations. It can be assumed that for a discrete system of equations, the conservation laws must also be strictly satisfied, which will increase the finite-difference solution stability and proximity to the exact one.

Marine Hydrophysical Institute of RAS, Sevastopol, Russian Federation H demyshev@gmail.com

H

Abstract

ISSN 1573-160X PHYSICAL OCEANOGRAPHY VOL. 30 ISS. 5 (2023)

523

The content is available under Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0) License

In [2], a difference analogue of Noether's theorem was constructed for the one-dimensional Euler equation, which, due to the transition to the space of grid functions, differed significantly from the original in additional conditions. The construction of the corresponding difference schemes is closely related to the property of their invariance, which definition was given in [3]. The property of invariance under transformation groups is a necessary condition for applying Noether's theorem to obtain conservation laws. Invariant difference schemes were built for second-order ordinary differential equations [4]. In [5], difference analogs for linear and nonlinear one-dimensional wave equations with symmetry conditions and conservation laws were obtained. For one-dimensional shallow water equations in Lagrangian coordinates, a new invariant finite-difference scheme was obtained. It has local laws of energy, mass, mass center, and momentum conservation [6]. When applied to computational geophysics problems, the law of total energy conservation is used as the norm of a discrete solution, which increases its stability. In [7-9], based on energy-stable schemes, the properties of the numerical solution of the Navier-Stokes equations were studied. Test problems demonstrated the accuracy and performance of the algorithms presented. A number of works were devoted to the CABARET scheme, which was introduced as a new three-layer explicit difference scheme with spatial splitting of the time derivative [10]. It was carefully developed at a high level for the one-dimensional transport equation without dissipation [11, 12], which made accurate simulation of the solution properties to the corresponding hyperbolic equation possible.

When solving oceanological problems, the main attention of researchers is focused on the calculation and analysis of difference analogs of energy characteristics, which are not always an exact consequence of the original discrete formulation and, strictly speaking, reflect model dynamics with errors.

As an example of the energy analysis use, the following works can be figured out. In [13], estimates of the spatial distribution and transformation of eddy energy in a model ocean basin are presented. The authors show that local generation and dissipation of eddy energy is compensated by so-called non-localized flows from other parts of the basin, especially along the periphery of eddy structures. In [14], based on an analysis of the distribution of the average kinetic energy of eddies calculated from satellite altimetry data, estimates of the spatial distribution of sources and sinks of eddy energy for actually observed mesoscale structures were carried out. The influence of the horizontal resolution of models on the description of energy cascades was studied in [15]. The authors prove that when calculating kinetic (KE) and potential (PE) energy budgets, forward and reverse energy cascades are more accurately described in a model that resolves structures with dimensions smaller than the baroclinic Rossby deformation radius. Model estimates of the seasonal variability of eddy kinetic energy carried out in [16], indicate that the seasonal increase in the thermal effect of the atmosphere causes a mesoscale variability intensification in the velocity field at the ocean surface. The Lorentz energy cycle analysis for the entire World Ocean based on reanalysis data was carried out in [17]. In this paper, the mean 524 PHYSICAL OCEANOGRAPHY VOL. 30 ISS. 5 (2023)

long-term estimates of energy components were obtained and the main differences between the energy cycles of the ocean and atmosphere were presented. The methodology for estimating eddy energy and its transformation mechanisms, proposed by Lorenz, is widely used for regional circulation studies. For example, the region of the Kuroshio Current is considered in [18], eddy variability in the Red Sea is studied in [19], and a detailed analysis of the mechanisms of mesoscale eddies formation in the Sea of Okhotsk is presented in [20].

For the Black Sea, a number of studies estimated the kinetic energy of eddies and currents based on observational data [21, 22]. Study of annual and interdecadal variability of available potential energy based on observational data for 19101998, was given in [23]. The energy balances of semi-enclosed seas (including the Black Sea) were calculated in [24]. It was also shown that the work of the buoyancy force was one of the main factors in the formation of mesoscale dynamics in semi-enclosed seas.

The present paper reviews invariants known for the hyperbolic system of differential equations of hydrodynamics of an ideal fluid (absence of friction, diffusion, external forces, and under the assumption of adiabaticity), but not obtained in discrete form as an exact consequence of the finite-difference formulation. In [25], difference equations for the rate of change of kinetic and potential energy corresponding to a numerical model of dynamics on the C-grid were written and analyzed. In a discrete formulation, to comply with the law of total energy conservation, it is necessary to adequately describe the work of the buoyancy force in the kinetic and potential energy budget equations. Its approximation depends on the difference type of density, which is generally determined by temperature, salinity, and pressure. This paper considers the case of a polynomial dependence on temperature and salinity. In [26], corresponding approximations of density were obtained. They ensure compliance with the total energy law in the case of representation of density in the form of a differentiable functional for a number of parameters. The present paper, as a continuation of these studies, considers discrete invariants obtained for a numerical model of sea dynamics under the condition of adiabaticity, absence of friction, diffusion, and external forces.

2. Conservation laws in differential formulation

Let us write down the equations of the ideal fluid model in the Boussinesq approximation and the incompressibility of sea water:

u, -(£ + f)v + wuz = -gqx-~(P'+ E)x, (1)

p0

v, + + f)u + wv2 =-gqy - — (P'+ E)y, (2)

^ P0 ^

u +v +w = 0, (3)

x y z ' v ;

z

P = gPo Ç + S \ pdv = gPo Ç + P' (4)

o

T + (uT ) x + (vT ) y + ( wT ) z = 0, (5)

St+(uS ) x + (vS )y + ( wS ) z = 0, (6)

P = G (T S). (7)

When z = 0, w = -çt ;

at the bottom at z = H(x, y) u = v = w = 0, Tz = Sz = 0. (8)

On the side walls: for meridional -

u = vx = 0, Tx = Sx = 0, (9a)

for zonal boundary sections -

v = uy = 0, Ty = Sy = 0, (9b)

Initial conditions: at t = t0

(T, S) = (T0, S0), u = u0, v = v0, Z = Z0 (10)

Here: Ç = vx -uv, E =

1 u2 + v2

In system (1)-(10), u, v, w are the velocity vector components directed along axes x, y, z, respectively; T, S, P are the sea water temperature, salinity, and pressure; f is the Coriolis parameter; Z is the reduced sea level; g is the free fall acceleration.

In (7), G(T, S) is a polynomial of arbitrary degree depending on temperature and salinity. In general, the equation of state includes pressure, but many models use relation (7), a version of which is recommended by the Intergovernmental Oceanographic Commission 1. In accordance with it, conservation schemes were built taking into account the polynomial dependence of density on temperature and salinity. For the further simplification of the recording of calculations, it is assumed that p0 = 1 g/cm3.

1 IOC, SCOR and IAPSO, 2010. The International Thermodynamic Equation of Seawater -

2010: Calculation and use of thermodynamic properties. Intergovernmental Oceanographic

Commission Manuals and Guides No. 56. UNESCO, 196 p.

2. 1. Equations of temperature, salinity, and density advection

The following notation is introduced:

h 1

< ^ >H = J^dz, < ^ >Qz = —JJfydxdy,

0 Qz Qz

1 H H

< ^ >V = — J JJ fydxdydz, V = J JJ dxdydz,

V 0 Qz 0 Qz

where Q z is the horizontal surface at z-level.

Due to the fulfillment of condition (3), equations (5) and (6) can be written in non-divergent form:

T + uTx + vTy + wTz = 0, (11)

St + uSx + vSy + wS2 = 0. (12)

From (3), (11), (12), and boundary conditions (8)-(9), it follows that for any K > 1, L > 1:

TtK + uTK + vTK + wTK = 0, (13)

t x y z 5 v '

SL + uSL + vSL + wSL = 0. (14)

Returning to the divergent form of writing relations (13), (14) and integrating over the domain, the first two invariants (K > 1, L > 1) are obtained:

< TtK >V +< c,tT0K >n° = 0, < SL >v +< qtSL >°° = 0. (15)

It should be noted that the equation (15) corresponds to:

<TK >V +<cT0K >Q0 = V(TK)dxdydz = 0,

V - Z Qz 1 Pi h

< SL >v +< cSL >Q° =--i ii (SL )dxdydz = 0.

V dt-z h\

Let G be a polynomial in T and S, then from expressions (11) and (12) the density advection equation is obtained:

p, + (u p )x + (vp ) ^ + ( w p )z = 0. (16)

The divergent form of expression (16) ensures the preservation of the density

1 d H

integral, similar to equation (15)---J JJ pdxdydz = 0.

V - z Qz

2. 2. Equations for the rate of kinetic and potential energy change

From equations (1)-(4), the equation for the kinetic energy budget is obtained: Et + [u (gq + P'+ E )]x + [v (gq + P'+ E)] , + [w( gq + P'+ E )]z = gp w. (17)

From equation (17), when integrating over the basin, it follows:

< Et > + < qtE0 + g(¿X > =< gwp > . (18) From the form of equations (17), (18), it can be seen that:

1 d h

< Et >V + < qtE0 >°0 =--1UEdxdydz = 0.

V - z nz

We assume that the potential energy has the form n = - g z p . Then, the equation for the rate of potential energy change is written as follows:

TT, + (un) x + (vn) y + (wn) z = - gpw. (19)

From equation (19), invariant III is obtained:

< Et + n, >V + < q E + g(¿X >°0 = ^^ 11|(E + n)dxdydz = 0. (20)

2 V dt--z h\

2. 3. Equations for the rate of dynamic potential energy change

An important physical characteristic of the sea and atmosphere dynamics is available potential energy, which was introduced by E. Lorenz in [27]. The integral energy, just like the total energy, under the condition that the processes are adiabatic, is conserved, or it is an invariant. Due to the state equation nonlinearity, ambiguity in determining the available potential energy takes place (in [27], a linear dependence of density on temperature was assumed). Therefore, in recent years, a number of works used the potential energy anomaly - dynamic potential energy - to analyze the results of numerical experiments [28].

In the continuous case, we assume that:

p( x, y, z, t ) = p* ( x, y, z, t ) + ps ( z ),

1 ( 1 ^

p* ( z) = -1 —1| p( x, y, z, t )dxdy

where r

T-

— ^ 4z —

dt, (21)

t is the integration time. In contrast to the definition given in [29], equation (21) assumes that p* is independent of time.

Let us introduce the notations: Dpe = - g z p .

Taking into account the continuity equation, substituting expansion (21) into expression (16), we obtain:

* * * * * /-v

p* + u px + vp + w pz + w p* = (22)

Transforming expression (22), the dynamic potential energy equation is obtained:

Dpe +(uDpe ^ + (vDpe) +(wDpe ^ -(gzwpS ) + gp*w = 0. (23)

Due to the absence of flows through the basin sides, taking into account expression (23), invariant IV is obtained:

<Et + <;tEo + Dpe >V + < g£■), >°° = -- f f[ (E + Dpe)dxdydz = 0. (24)

2 V dtJz

3. Conservation laws in a discrete model

We approximate a basin with an uneven bottom by boxes with centers corresponding to integer index values i, j, k (i = ii,..., iN, j = j!,..., jM, k = 1,..., Kj), edges - i + 1/2, j + 1/2, k + 1/2. The horizontal dimensions of boxes (hx,hy) are

constant; vertically, uneven approximation

( K = zk+1/2 - zk-1/2 , hzk+1/2 = zk+1 - zk ) is used.

Difference operators have the form (similarly for j, k )

£ = <U,jk + Wjk s ^ k = <U, jk - km,jk V2 ^ = g^ + g2^

T), j, k 2 xn j, k h x, y i, J, k xTi, j, k yTj, j, k>

x

(25)

(*>* =T~"ZtjAhy, (klF = ^Z Zki,jkfthxhy, V = Z Zh,.

i2k i, j V i, j k=1 i, j k=1

At horizons zk, temperature, salinity, and horizontal velocity components are

calculated, at horizons zk+1/2 - vertical velocity [25]; Qk is surface area at horizon

zk. The distribution of variables is shown in Fig. 1.

Let us consider equations that are continuous in time. A finite-difference system of model equations with second-order accuracy in spatial variables (up to a non-uniform step) was written out repeatedly [25].

3. 1. Invariants of differential-difference equations of temperature, salinity, and density advection

Let us write the state equation at point (i, j, k) as a polynomial:

N M

Pi,j,k = G(Tuk, suk ) = Guk =ZZa«,Tj,kS:j,k, (26)

n=0 m=0

where n > 0, m > 0, anm are the constants. Note that N, M do not have to be integers.

F i g. 1. Distribution of variables in box (i, j, k)

To preserve the spatial integral of the density under the condition of adiabaticity and the absence of external sources, it is necessary to obtain an approximation of the nonlinear terms in such a way that, along with T and S, TK and SL are preserved, respectively:

L JL

V dt

H i is H

11| (TK ) dxdydz = 0, — — 11| (SL L dxdydz = 0

-z

(27)

-Z

As a result, the volume integral of each member of series (26) is limited due to the positive temperature and salinity.

The continuity equation in discrete form at point (i, j, k), taking into account notation (25), has the form:

5 xui, j ,k + 5 yvi, j ,k + 5 zw,, j ,k =

(28)

To derive the corresponding approximations that ensure the preservation of two moments, the first and the highest ones, the following method is used.

Let us introduce new unknowns TMI2jk ,SMI2jk ,p+1I2jk, at points (i+1/2, j, k), at points (i,j + m, k) - Tj+1|U ,Sij+1|2)k ,Pj+1|2)k, at Points (h j k + 1/2) -

Ti, j,k+1I2' Si, j,k+1I2' pi, j,k+1I2 ■

Let us write down the equations of temperature and salinity advection at point (i, j, k):

^+) + 8, ) + 5, (ju.k) = 0, (29)

dS, n

-¡tT + 8, (Ui, jkSuk ) + 8, (Vi, j,kSi,j,k ) + 8, (Wi, j,kSi,j,k ) = 0. (30)

Further, all reasoning and calculations carried out for temperature will be

identical for salinity. Note that Ti+j,2, j,k , T .+H2,k , T, j,k+1I2

and Si+1,2, j,k, Si , j+1I2,k , Si, j,k+1I2 haVe not been determined yet-

Taking into account discrete continuity equation (28), equation (29) is rewritten in the following form:

Ui+1I2, j ,k (T+1I2J, k — Ti ,j ,k ) — Ui-1I2, j ,k (Ti-1I2, j ,k - Ti, j ,k )] h- +

i, j ,k

dt

-1

Vi, j+1I2,k (T j+1I2,k Ti, j,k ) Vi, j-1I2,k (T j-1I2,k Ti, j,k )]

(Ti, j ,k+1I2 - T,j ,k )- Wi, j,k-1I2 (Ti, j,k-1I2 - T,j ,k )]( hk y = 0. (31)

w,j,k+1I2 \Ti,j,k+1I2

In a similar way, we proceed with equation (30). Let functional Q, for example, on temperature and its derivative be written as follows:

Qi,j,k=Qi,j,k (x,j,k x dQj=Qi,j¿¡¡if-, where Qi,j,k=d^.

We multiply equation (31) by QQ.ijk and transform the result taking into account continuity equation (28). We get the equation for Q. k :

ddQj + {Ui+1I2, j,k [Q, jk + (T'+H2, j,k - Ti,j,k )] - Ui-1I2, j,k [Qij - Uj,k (Ti,j,k - T'-H2, j,k )]} +

+ {Vi, j +1I2,k [Qi,j,k + Qi,j,k (Ti, j+1I2,k - Ti,j,k )] - Vi, j-1I2,k [Qi, j,k - Qi,j,k (Ti, j,k - Ti, j-1I2,k )]} hy +

+ W j,k+1I2 [Qi,j,k + Q j,k (T,j,k+1I2 - T,j,k )] - W j.k- 1I2 [Qi,j,k - Q j,k (T;j,k - T,j,k- 1I2 )]}(K ) 1 = 0

(32.1)

Assume that Qi j k satisfies the advection equation, which is written in the form:

dQt

i, j ,k

dt

[(Qi

i+1/2,j ,k

- Qi, j,k )u

i, j,k ' i+1/2, j,k

- (Qi

i -1/2, j ,k

- Qi, j,k )u

i, j,k ' i-1/2, j,k J x

h"1

+[(Qi+ , j+1/2,k i, j,k / i, j+1/2,k - (Qi ,j-1/2,k i^i,j,k/ i,j-1/2,k h + +[(Qi, j,k+1/2 - Qi,j,k M, j,k+1/2 - (Qi,j,k-1/2 - Q:,jk K, j,k-1/2](hz ^ +

+Qi,j ,k (gxUi, j,k + gyVi,j,k + g ) = 0.

Then from equations (32.1) and (32.2), an expression, e.g., for T

(32.2)

i+1/2, j ,k

(similarly for Ti

Ti.i.k+1/2 ) is obtained:

i, j+1/2,k' * i, j,k+1/2

T

(a

T

i+1,j,k i+1,j,k

Qi +1, j,k ) (Qi,j,kTi, j,k Qi, j,k )

i+1/2,j ,k

Qi+1, j ,k Qi, j ,k

Assume that Qi j k = TKjk. Then temperature (salinity in the same way) is approximated on the sides of box (i, j, k) as follows:

T

K

¿+1/2, j,k

K-1

K

K

K

K

T+1,j,k ~Ti,j,k fTrK-1 fTrK-1

V i+1,j,k i,j,k

TK

Ti, j ,k+1/2

K -1

K

T, j+1/2,k K

K-1

K

K

K

K

T,j+1,k - t/,j,k rpK-1 rpK-1

V i,j+1,k i,j,k J

K

T& - tK \

Ti, j,k+1 _ Ti, j,k TK-1 - tK-1 V i,j,k+1 i, j,k J

in case Q. ,, = SL,, SL,,„ ,

j ,k i, j ,k /+1/2, j ,k

L -1

SL

j+1/2,k

L -1

-S

L \

i+1,j ,k i',j ,k sL-1 - SL-1 V^+U,k Sij,k J

^ SL

i, j+1,k

-S

l ^

i,j,k

SL-1 - sL-1

Vs/,j+1,k ^ij j

SL

5 °i, j,k+1/2

L -1

^ SL

L

i, j ,k+1

-S

i, j ,k

sL-1 - sL-1

V Si, j,k+1 Si, j,k J

(33)

In case of hydrostatic instability of the fluid column, convective mixing procedure is used. It leads to temperature and salinity equalization at adjacent horizons. The relationships on the vertical sides of the box, which follow from equation (33), can be used:

T

K-1

i, j ,k+1/2

where

K

Tjk +

T

K-1 i, j,k+1

®(T/,k) J

S

' ,k+1/2

K-2

L-1

Sj +

sL-1 j,k+1

j,k+1, Si, j,k )

■j,k/ /

(34)

L-2

® (Ti, j ,k+1, Ti, j ,k ) = i Ti J ,k+lTi "j ,k , ^ ( Si, j, k+1, Si, j ,k ) = i Sf,j !k-1 Si

i, j ,k'

n=0

m=0

Thus, from the approximations of expressions (33) and (34), it follows that when integrating over the entire domain, T, S and TK, SL (K > 2, L > 2) are preserved, that is, the following conservation laws are satisfied:

f dTKk 1 \f < j^K i dSi,j, k 1 F.f

{~dTJ "{~drT,K1/2J "0, {^TJ Sj,1/2

These expressions correspond to integrals (27) (invariants I and II). To ensure the total energy conservation law, it is necessary to describe the buoyancy force work adequately in the equations for the rate of kinetic and potential energy change. For this purpose, it is necessary to obtain an equation for density from heat and salt advection equations (29), (30), taking into account the fulfillment of relations (33) and (34), which should have a divergent form.

Let the density advection equation be written in the form in which

P+1/2, j,k ,Pi, j+1/2,k ,Pi, jM1/2 are the unknowns:

dPi, j,k dt

A , jk P,, jk ) + ^ (v,, jk P i, jk ) + 5 , (w,, jk p i, jk ) = (35)

Density in point (i, j, k) has the form pi,jkk = Gt,j,k (T j,k, S;, j,k). We believe the following to be fulfilled:

lpr=( Gi,M )T dJf-+(Gi,j,k l^f-,where (Gi,M )T=, Km )S ==j

We multiply equations (29) by (G^t and (30) by (Gj)S. Carrying out the appropriate transformations, taking into account continuity equation (28) and the requirements of the divergent form of equation (35), we obtain the relation, for example for pi+1/2,. ,k (for p

i, j+1/2,k ,Pi, j k+1/2 - similarly)

P,

i+1/2, j k

T

i+1/2, j ,k

(Gi+1,j k ))t + (Gi, Jk ))t

A

(Gi+1, j,k )' T+1, j k + (Gi, jk )' Ti, jk

s..

¿+1/2, jk

(G+ jk ))s + (Gi, jk) )s

(G+1, j,k )ss S+1, j,k + (G, j,k )sss, j,k

GM,j,k + G,, jk

(36)

The divergent form of the density advection equation, as a consequence of equations (29), (30), is provided by approximating the density on the box sides in the form of formula (36).

In genera1, Pi+1/2, j,k , j+1/2,k ,Pi,j,k+1/2 should have the form:

k -1 t k1,, - tk

i+1, j,k i, j,k a nTn-1 Sm-1

1 Z_lZ_l n ,m-1'" i+1/2, j ,k°i+

i+1/2, j ¿-if rpK-1 rpK-1 /1 / / n ,m-1"^ i+1/2, j ,k i+1/2, j ,k K Ti+1, j,k - Ti, j,k n=1 m=1

L -1 sL+j - sf NM

___ijL ii a mTn-1 sm-1

T SL-1 - sL-1 -1,m i+1/2, j,k i+1/2,j,k

L Si+1,j,k Si,j,k n=1 m=1

ii (an,m-1nTi+1/2,j,kSi+l/2,j,k n=1 m=1

+an-1,mmTi +-/2, j ,kSi+1/2, j ,k ) + iian,mTi ?1/2, j ,kSi+1/2, j ,k > (37a)

TK - TK N M _y

D = K 1.+1, J ,k ., J ,k iia nTn-1 Sm-1 y

Ki, j+1/2,k „ rpK-1 y,K-1 Z_lZ.J n,m-1'" i, j+1/2,k i, j+1/2,k

K Ti+1, j,k - Ti, j,k n=1 m=1

L -1 Sf+1^ - St

N M

UM_iJL ii a mTn-1 Sm-1

T SL-1 - SL-1 n-1,m"lli,j+1/2,ki,j+1/2,k

L Si,j+1,k Si,j,k n=1 m=1

y

-iiKm-1nTt:+1/2,kS™-1

j+1/2,k i, j+1/2,k

N M

+an-1,mmTi ,j+1/2,kSi,j+1/2,k Ti nj+1/2,kSi, j+1/2,k .

(37b)

n=1 m=1

K - 1 TK - TK N M

., j ,k+1 ., j n"T n-1 S™-1

1 Z_lZ_l n,m-1'" i, j,k+1/2°i,

i, j,k+1/2 „ rpK-1 rpK-1 ' ' / ' n,m-1"-1 /, j,k+1/2i, j,k+1/2

K Ti, j,k+1 - Ti, j,k n=1 m=1

L -1 Sf ,^+1 - Sf NM

_JL iia mTn-1 Sm-1

1 AuAu Un-1,m'nli, j,k+1/2°i,

T <?L-1 ^L-1 ' ' n-1,m i, j ,k+1/2 i, j ,k+1/2

L Si, j,k+1 - Si, j,k n=1 m=1

N M

ii (an,m-1nTi,j,k+1/2Si,j,1+1/2 n=1 m=1

+an-1,mmTinjjк+1/2 SiJ ,k+1/2 ) + iian,mTi "j ,k+1/2 Si'"j ,k+1/2 • (37c)

n=1 m=1

From expressions (33), (34), it follows that when integrated over the entire region, Tt j k and TtKjk, S.k and SL k remain unchanged. Note that relations

(37) do not depend on the conservation scheme for temperature and salinity in a certain sense. Namely: formulas (33) and (37a)-(37c) can contain different values of N, M, K, L . For example, the approximation of the temperature advection

equation ensures conservation of TL (expression (33)), expressions (37) lead to a divergent form of the density advection equation, and the equation of state has the highest degree of T equal to N.

D

Thus, approximations (37a)-(37c) provide a divergent form of density advection equation (35) and satisfy

I F r , > U

dp■ k[ i dc. . ] 0 1 d H CC

1,', 1,' p1 j 1/2 I = 0, corresponding to--J JJ (pp)dxdydz = 0

J',j,k I + )"^,j

dt | [ dt r1,j,1/2J ' ' ~ F dt _z U

(invariant III).

Note that, in contrast to the CABARET scheme [11, 12], approximations of temperature, salinity, and density ensure the fulfillment of the corresponding conservation laws and have a "nonlinear" form of (30), (37a)-(37c). On the edges of the box, relations (33) and (37) involve the T and S values at only two points.

Analysis of the derivation of relations (36) shows that the obtained result is easily generalized to the general case of dependence of G on r (r = 1,..., R, where R is an integer) of Fr functions:

Gk=G (Fjk,Fjk). (38)

Let G be a differentiable functional for each Fr

dG

G' 1 jk

dF.r k 1, j, k

and corresponding derivatives are limited. It is assumed that each Fr at point (1, j, k) satisfies the advection equation:

- + 5 u .tFr.t + 5 v .tFr.t + 5 w. .tFr.t = 0, (39)

dtt x 1,j,k 1,j,k y i,j,k i,j,k z 1,j,k 1,j,k ' v '

where r = 1, ..., R. For each equation (39) (r = 1, ..., R), carrying out the corresponding transformations, the following is obtained:

dF

r

_LLL + iF -Fr )u h -(Fr - Fr )u h~l +

d^ ^ 1+1/2, j,k i, j,k '1+1/2, j,k x V^ 1-1/2, j,k 1 1, j,k'1* 1-1/2, j,k"x

+(F,j+1/2,k - Fj )V1,j+1/2,kh'h - (Fij-1/2,k - Fj,k )V1,j-1/2,kK + (40)

+^+1/2 -Fjk )W, )-1 - (F,r;,k-1/2 -Fjk )W, ,k_1/2(hh )-1 + +Fjk (5 u + 5 v + 5zWj) = 0.

dG,. k _ , dF,1 r k The following relation holds: —— = > (G ,) r—'—. Then, after

* dt r 1,j,k'F dt

transforming equation (40), taking into account definition (38), multiplying it by G, r and summing the attained result, we have:

dG

—^ + y (G k) r[(Fi+1/2 . k-Fr k)u.+1/2 .h - (F' . k-F.'k)u.1/2 . kh +

dh ' 1 v i, J,k' Ar LV- i+1/2,j,k i,j,ks i+1/2,j,k x V i-1/2,j,k i,j,ks i-1/2,j,k x

+(F'k -fr.k)V, +1/2 h - (F r 1/2 k -f r k)V,. 1/2 ,hl +

V i,j+1/2,k i,j,k' i,j+1/2,k y V i,j-1/2,k i,j,k' i,j-1/2,k y

-KkW ,k+1/2(hk)-1 - (Fjk-1/2 -Fjk)wt j,k-1/2(hZ)-1 + +Fj (5xU,. k + SyVi . k + 5. k)] - 0. (41)

Assume that the advection equation holds for functional G. Taking into account expression (41), a recurrence relation for G on sides (i + 1/2, j, k) of box (i, j, k) (for Gij+ink and Gijk+1/2- similarly) is obtained:

Gi+1/2 ,j,k y \_Fi+1/2,j,k

rG.„ ., +G.., ^ (G.+1 .kYF (G k)' rFr

' V i +1,j,k' Fr i+1,j,k V i,j,k'Fr i,i

i+1,j,k i,j,k

Jf'

G... .k + G .k

+ i+i,j,k i,j,k (42)

2 '

This relation can be interpreted as a finite-difference analogue for functional G on the sides of box (i, j, k). Approximation (42) leads to a divergent form of the advection equation of functional G:

dG, . k

——— + 5 u. ,G ., + 5 v. ,G ., + 5 w. ,G ., = 0

dt x i,j,k i,j,k y i,j,k i,j,k z i,j,k i,j,k

and, thereby, ensures that the volume integral of G is equal to zero. The resulting property exactly corresponds to the continuous formulation. Note that difference relation (42) does not depend on specific form F++ y2JJc, Fi"j+1/2,k, F*JMm.. As an example, the equation of state from work 1 is given in the form:

p. .. =y y a Tn.kSm.k. (43)

"i,j k / ' / ' nm i,j,k i,j,k v '

n=0 m=0

In accordance with equation (43), it is required that in the adiabatic approximation and the absence of external forces, the volume integrals of T 5 k and

Sfjk do not change with time. Then expression (37a), taking into account expression (43), is written in the form (for j, k similarly):

4 T4 + T3 T + T2 T2 + T T3 + T4 5 f -x -x

411+1, j,k+1i+1, j k1!, j,k+1i+1, jMi, j ,k+1i+1, jMi, j,k+1i, ja nT x + a nT S x) +

H+1/2,j,k 5 T3 + T2 T + T T2 1 T3 n,0 i+1/2,j,k ~'~an,1nl i+1/2,j.k^ij,k J

5 +1, j ,k + Ti+1, j,kTi,j,k + +1, j ,kT i, j,k + Ti, j ,k n=1 V J

S;,. : U + S; : , ^ ~ x

+ ' r. y an-11Ti+1/2, j,k y I an,0nTi+1/2, j,k + an-1,1Ti+1/2, j,kSi+1/2, j,k +an,1Ti+1/2, j,kSi+1/2, j,k

2 n=1 n=1 V

3. 2. Discrete total energy conservation equation

In the absence of friction and external forces, the discrete kinetic energy equation has the form [25]:

dE,

1, j,k

dt

+ 5.

1, j ,k 1, j ,k

+ gC,, + (K ,k )1

+ 5,

1,j ,k 1 ,j ,k

+ gc, ■■ + (K,;-,k)

+5,

Wj ,k |g c,j + (p!/,k)

= g

W, j,k P, j,khk hk

1, j,k>

+ 5.

w.

1, j ,k

^EU )1,j,k + WUky (EV )1,j,k

(44)

Equation (44) uses the following notations:

1, j ,k+1/2

U- ,,u

( e ) = 1,j ,k+1 1,j ,k ( F )

V^u/1, j,k+1/2 _ „ ' V v /

1, j,k+1 1, j,k j JJ, ) = 1, j,k+1 1, j,k

v /1, j,k+1/2 _ _

E1, j ,k =

(U1, j,k )2 + (V1, j,k )2 2

The equation for density in the absence of diffusion and in the adiabatic

approximation has form (35X where p. j,k , p1+1/2, j,k , p1, j+1/2,k > p1, j,k+1/2 satis^ to relations (37a)-(37c).

Then the equation for the rate of potential energy change n 1 jk = -gz k p1 j k is written in the following way [25]:

d n

1, j ,k

dt

5x (u1, j,kn 1, j,k ) + 5y (V, j ,kn 1, j,k ) + 5z (W1, j,kn 1, j,k ) = -g

w. .p ,h

. (45)

In equation (40), n 1+1/2, jk =-gzk P1+1/2, j ,k ,n 1, j+1/2,k =-gzk P1, j+1/2,k ,

n,

*-1, j,k+1/2 = -gz k+1/2 p, j,k+1/2.

Note that the difference analogue of the buoyancy force work in equation (45) is identical to the term in equation (44) for Et j k in the case of a divergent form of

equation for density (35).

Let us integrate equations (44) and (45) in a difference sense over space, then the discrete analogue of the total energy satisfies the equation (invariant IV):

dEj dP1,j,k . , j -gz,,-— I

dt

dt

-y

-y

dCu (u1,j,1) (v1,j,1) +gc dCu

dt

dt

31, j

dt

= 0. (46)

The total energy conservation law (46), being an analogue of continuous invariant (20), is satisfied provided that the approximation of the buoyancy force

work (the right side of equations (44) and (45)) is exactly consistent in both equations.

3. 3. Discrete equation for the rate of dynamic potential energy change

Assume that (^k+1/2 is the surface square at horizon z k+1/2):

1 ^2 1 j ^ Pk+1/2 = —— y^- ypi,j,k+1/2hxhy h ' where h is the time step.

(t2 t1) k k+1/2

V ' ,J

Then suppose:

-z -z -z

P, j ,k = p,j ,k + Pk , p+1/2, j ,k = p+1/2, j ,k + Pk , j+1/2, k = j+1/2,k + Pk Pi, j,k+1/2 = p, j,k+1/2 + P/fc+1/2 • (47)

Taking into account continuity equation (28), we obtain:

d Pj

dt

" 5* (u-jk Pj ) + 5, ( v,j,k Pj ) + Sz ( wt,j,k ftj ) +

+PkZ (5 Am + 5 Am + 5 zW,M ) + Z 5 zPk =

(48)

Let DipeJk = -gzk p* jk . Then the consequence of notations (47) and equation

(48) is the equation for dynamic potential energy (Dpej k) in the adiabatic approximation, in the absence of diffusion and external sources:

dDpj

i, j ,k

dt

+

5 ( ,k )+5 ( v^z ,k )+5z ( w,,k<; ,k )] + ( ^^ ^^ S zft- ) =

=- g

i k j hzW:,j,k Pj

Equation (49) uses the following notations:

pj _ * pj _ * pj _ *

+1/2, j,£ = — gZ k pi+1/2, j,£, ai, j+1/2,/t = — gZ k pi, j+1/2,i, ai, j,£+1/2 = — gZ k+1/2 Pi, j,£+1/2. Then the following property holds:

(49)

dEi,j,k - z dP*,j,k y dt k dt

("i,;i) , d?i,; (Vn)

dt

. -.,j v i,j,i/ =i,j

+ —- —— +gÇi j

dt

*,j dt

= °. (50)

Approximations (47) ensure the fulfillment of equation (49) and, therefore, (50) (invariant V, which corresponds to its continuous counterpart (24)). Note that conservation law (50) can also be obtained in a more general case - the dependence

of the average density over time pS = pS ( z, t ).

4. Numerical experiments

Let us consider examples of the Black Sea circulation calculation for 2011 [29] using some of the approximations obtained.

In the numerical experiments, a uniform step along the horizontal coordinates was set to 1.6 km; 27 horizons were used vertically with condensation in the upper layer of the sea. The runoff of the Black Sea rivers was taken into account in accordance with [30] and amounted to about 340 km, of which the rivers of the northwestern part of the sea (the Danube, the Dniester, the Dnieper, and the Southern Bug) account for approximately 78%, the rivers of the Caucasus (the Rioni and smaller rivers) - 13%, and the rivers of Turkey (the Yeshil-Irmak, the Kyzyl-Irmak, and the Sakarya) - 5%. The rivers of the Caucasus, due to their large number, were represented along the Caucasian coast in the form of three sources. The salinity at river mouths was assumed to be zero. The water temperature at the mouths of rivers, except for the rivers of Turkey, was set from [31]. The temperature of the Turkish rivers was set equal to the temperature of the coastal sea waters.

It was assumed that in the Upper Bosphorus current the temperature and salinity are the same as in the sea. In the Lower Bosphorus stream, the temperature was taken to be 16°C and salinity - 35%o.

To set the atmospheric forcing, SKIRON data for 2011 were used; vertical mixing was described based on the Mellor-Yamada theory [31]. The initial conditions in this calculation corresponded to January 1, 2011.

Three numerical calculations were carried out. They varied in difference schemes for approximating the temperature and salinity advection equations. In the first (I) calculation, a "traditional" scheme was used, which ensured the preservation of T jk, T2j k, St jk, Sf j k, in the second (II) one -

T .., T5.,, S , S3.,, and in the third (III) one - T .., Tw,, S , Sj .

1, j, k' 1, j, k' 1, j, k' 1, j, k' v ' 1, j, k' 1, j, k' 1, j, k' 1, j, k

The choice is based on the following considerations. The second experiment corresponds to the maximum degree of temperature and salinity in the modern equation of state, the third one demonstrates the effect of the approximations used.

The initial conditions for the second and third calculations correspond to February 1, 2011 according to the results of the first experiment.

The integration time for the second and third calculations was 20 days. All figures are as of February 21, 2011.

Fig. 2 demonstrates that the basin was covered by an extensive cyclonic gyre with two centers in its western and eastern parts. Mesoscale eddies ranging in size from several kilometers to tens appeared and evolved between the shore and the Black Sea Rim Current. There is a qualitative agreement in the level fields (Fig. 2), and quantitative differences in the level structure among three calculations are insignificant and observed in areas of intense temporal variability of the Black Sea Rim Current. The qualitative agreement of the results indicates the correctness of the approximations used.

b

F i g. 2. Reduced sea level in experiments I (a), II (b), and III (c)

During this period, on the sea surface (Fig. 3), a zone of cold water is observed on the northwestern shelf, and that of warm water - in the southwestern anticyclone area. Upwelling is observed near the Kerch Strait and a thin alongshore strip in the southeastern corner of the basin. These features were more clearly obtained in experiments II and III (Fig. 3, a, b). In the southeastern anticyclone area, narrow frontal zones, which are more clearly observed in calculations II and III, are distinguished.

t, °c

28 °E 30°E 32 °E 34°E 36°E 33°E 40 °E b

t °c

M°E 30°E 32°i 34°E 3ffi 38°E 40°E

c

F i g. 3. Temperature at the 3-m horizon in experiments I (a), II (b), and III (c)

The greatest difference in temperature fields is observed in the upper layer of the sea in coastal areas (Fig. 4). The areas identified are the northwestern shelf, the upwelling zone near the coast of Crimea, and the alongshore strip near

the Anatolian and Caucasian coasts. The maximum difference between the temperature values in experiments I and II, I and III (Fig. 4) in absolute value is 2 °C for both options, and the average is 0.14 °C.

46°N-

44°N-

42 °N-

AT," C

1---7---1-1---1---1-----T"

28°E 30°E 32°E 34°E 36°E 38°E 40°E

.

a

46°N-

44°N-

42 N

AT", °C

28°E SCE 32°E 34°E SS^E 38°E 40°E

b

F i g. 4. Difference between the temperature values at the 3-m horizon in calculations I and II (a) and I and III (b)

The salinity field has a structure that corresponds to an extensive cyclonic circulation covering the deep-water part during this period, and contains areas of less saline water in the coastal zone. A clearly expressed tendency is observed: when using schemes with invariants of a higher degree (calculations II and III - see Fig. 5, b, c), the effect of deep waters upwelling in the central part and downwelling of less saline waters along the cyclonic circulation periphery is enhanced. This is evidenced by the structure of isoline 18.35%o, which serves as a conditional marker. Its relatively smooth structure demonstrates that increasing the order of the invariant does not increase the computational noise in the model, manifested in the form of a small-scale component.

S, %0

23 E 30°E 32 E 34 E 36°E 38°E 40"E

16.4 18 3 18 2 18 1 18 17. S 17. s 17.7 17.6 17 5 17 4 17 3 17.2 17.1 17

44°N-

¡"E 30°E 32°E 34° E 36°E 33°E 40°E

1 a.A

183 18.2 18.1 1 8 17.9 1 7.S 1 7.7 1 7.6 1 7. 5 1 7. 4 1 7.3 17.2 17,1 17

a

b

c

F i g. 5. Salinity at the 30-m horizon in experiments I (a), II (b), and III (c). Isoline 18.35 %o is highlighted in red

In the salinity field, the greatest difference is observed in the river inflow area and in the Bosphorus region. In the central part of the sea, its negative values are

PHYSICAL OCEANOGRAPHY VOL. 30 ISS. 5 (2023) 543

observed. It corresponds to the excess of salinity in experiments II and III compared to calculation I, and positive values at the continental slope of the northwestern shelf and along the coastal area of the northern and western parts (Fig. 6). The extreme values are approximately 1 and 2%o for both calculation options.

a

b

F i g. 6. Difference between the salinity values at the 30-m depth in calculations I and II (a) and I and III (b)

When implementing schemes with nonlinear invariants, comparative analysis reveals one peculiarity, which is the sharpening of transverse gradients in the frontal zones. Therefore, it can be assumed that in calculations II and III, the vertical velocity is reproduced more adequately. And first of all, it is necessary to consider the presence of a small-scale component, which can contain a large error.

From a comparison of Fig. 7, a and b, it follows that when using patterns with invariants of a degree greater than 2, the density of two to four step perturbations decreases. For any numerical model, the variability of fields on these scales is distorted due to the incorrect value of the group velocity.

45°N-

44°N-

43°N-

42°N-

41°N-

28 °E 29 °E 30 °E 31°E a

F i g. 7. Vertical velocity at the 400-m depth in calculations I (a) and III (b) (western part of the sea)

A similar picture is observed in the eastern part of the sea (Fig. 8). Large-scale features in the vertical velocity structure in all three calculations are approximately the same, but the intensity of small-scale variability is lower in calculations II and III.

45° N

43°N"

42° N-

41°N

I

H

£ CJ

£ CJ

_o

-0.005

36 E 37 E 38 E 39 E 40 E 41"E

F i g. 8. Vertical velocity at the 400-m depth in calculations I (a) and III (b) (eastern part of the sea)

5. Conclusion

The results obtained demonstrate the possibility of constructing discrete analogues of a continuous problem that have a set of nonlinear invariants. The method of introducing an overdetermined grid used in the present paper allows obtaining necessary approximations to ensure a number of conservation laws. In the absence of friction and external forces in the adiabatic approximation, the expressions obtained ensure the conservation of

TK (K > 1), SL (L > 1), p, (E + n), (E + Dpe).

It is easy to see that expressions (31.1), (31.2), and (32) can be generalized to the case of an arbitrary form of differentiable functional Qt k depending on Tt jk, and provided that the continuity equation is satisfied. When using patterns that

b

b

a

have TiKj k, Sfj k invariants, where K and L can be large, at least two questions arise.

First, it is necessary to estimate the calculation error, since the number of operations grows correspondingly as the polynomial degree increases, which can lead to an increase in the order of the approximation error. Secondly, it is not clear how well these approximations apply when describing the inflow of rivers, exchange through straits, and action of external forces (taking into account the atmospheric effect on the sea surface).

When deriving the difference equation for density advection from functional Gi jk, it should be differentiable from Tt j k and Sf Jk. It is easy to generalize

the results obtained to the case when Gi j k depends on R functions (formula (42)). Assuming that each variable of the functional satisfies an advection equation of type (35) in a velocity field satisfying expression (28), the equation for Gt k is

obtained in divergent form.

An attractive feature of the obtained approximations is the ability to use independently patterns that ensure conservation of TiKjk, Sfj k at K > 2, L > 2, and

the total energy in the equation of state, which has the form of a polynomial of temperature and salinity of degree N > 1, M > 1.

The analysis of the results obtained indicates that with an increase in K and L, three effects take place: gradients in the frontal zones in the temperature field become more acute, processes of water upwelling in the center of the sea and descent along the periphery intensify, and the density of small-scale features in the vertical velocity field decreases. How connected they are and how they affect the accuracy of calculations are the next questions that require a separate study.

REFERENCES

1. Noether, E., 1918. Invariante Variationsprobleme. Nachrichten von der Gesellschaft der Wissenschaften zu Gottingen, 2, pp. 235-257 (in German).

2. Dorodnitsyn, V.A., 1993. A Finite-Difference Analogue of Noether's Theorem. Doklady Mathematics, 38(2), pp. 66-68.

3. Samarskij, A.A., Mazhukin, V.I. and Matus, P.P., 1997. Invariant Difference Schemes for Differential Equations with the Transformation of Independent Variables. Doklady Mathematics, 55(1), pp. 140-143.

4. Dorodnitsyn, V.A. and Kaptsov, E.I., 2014. Invariant Difference Schemes for Second Order Ordinary Differential Equations Possessing Symmetries. Keldysh Institute Preprints, (16), 42 p. (in Russian).

5. Cheviakov, A.F., Dorodnitsyn, V.A. and Kaptsov, E.I., 2020. Invariant Conservation Law-Preserving Discretization of Linear and Nonlinear Wave Equations. Journal of Mathematical Physics, 61(8), 081504. doi:10.1063/5.0004372

6. Kaptsov, E.I., 2019. Numerical Implementation of an Invariant Scheme for One-Dimensional Shallow Water Equations in Lagrangian Coordinates. Keldysh Institute Preprints, (108), 28 p. doi: 10.20948/prepr-2019-108 (in Russian).

7. Lin, L., Ni, N., Yang, Zh. and Dong, S., 2020. An Energy-Stable Scheme for Incompressible Navier-Stokes Equations with Periodically Updated Coefficient Matrix. Journal of Computational Physics, 418, 109624. doi:10.1016/j.jcp.2020.109624

8. Chen, H., Sun, S. and Zhang, T., 2018. Energy Stability Analysis of Some Fully Discrete Numerical Schemes for Incompressible Navier-Stokes Equations on Staggered Grids. Journal of Scientific Computing, 75(1), pp. 427-456. doi:10.1007/s10915-017-0543-3

9. Lin, L., Yang, Z. and Dong, S., 2019. Numerical Approximation of Incompressible Navier-Stokes Equations Based on an Auxiliary Energy Variable. Journal of Computational Physics, 388, pp. 1-22. doi:10.1016/j.jcp.2019.03.012

10. Goloviznin, V.M. and Samarsky, A.A., 1998. Finite Difference Approximation of Convective Transport Equation with Space Splitting Time Derivative. Matematicheskoe Modelirovanie, 10(1), pp. 86-100 (in Russian).

11. Afanasiev, N.A., Goloviznin, V.M. and Solovjev, A.V., 2021. CABARET Scheme with Improved Dispersion Properties for Systems of Linear Hyperbolic-Type Differential Equations. Numerical Methods and Programming, 22(1), pp. 67-76. doi: 10.26089/NumMet.v22r105

12. Goloviznin, V.M., Maiorov, P.A. and Solovjev, A.V., 2022. Validation of the Low Dissipation Computational Algorithm CABARET-MFSH for Multilayer Hydrostatic Flows with a Free Surface on the Lock-Release Experiments. Journal of Computational Physics, 463, 111239. doi: 10.1016/j.jcp.2022.111239

13. Grooms, I., Nadeau, L.-P. and Smith, K.S., 2013. Mesoscale Eddy Energy Locality in an Idealized Ocean Model. Journal of Physical Oceanography, 43(9), pp. 1911-1923. doi:10.1175/JP0-D-13-036.1

14. Shang, X., Xu, C., Chen, G. and Lian, S., 2013. Review on Mechanical Energy of Ocean Mesoscale Eddies and Associated Energy Sources and Sinks. Journal of Tropical Oceanography, 32(2), pp. 24-36.

15. Kjellsson, J. and Zanna, L., 2017. The Impact of Horizontal Resolution on Energy Transfers in Global Ocean Models. Fluids, 2(3), 45. doi: 10.3390/fluids2030045

16. Rieck, J.K., Boning, C.W., Greatbatch, R.J. and Scheinert, M., 2015. Seasonal Variability of Eddy Kinetic Energy in a Global High-Resolution Ocean Model. Geophysical Research Letters, 42(21), pp. 9379-9386. doi:10.1002/2015GL066152

17. Von Storch, J.-S., Eden, C., Fast, I., Haak, H., Hernandez-Deckers, D., Maier-Reimer, E., Marotzke, J. and Stammer, D., 2012. An Estimate of the Lorenz Energy Cycle for the World Ocean Based on the 1/10° STORM/NCEP Simulation. Journal of Physical Oceanography, 42(12), pp. 2185-2205. doi:10.1175/JPO-D-12-079.1

18. Yang, Y. and San Liang, X., 2018. On the Seasonal Eddy Variability in the Kuroshio Extension. Journal of Physical Oceanography, 48(8), pp. 1675-1689. doi:10.1175/JPO-D-18-0058.1

19. Zhan, P., Subramanian, A.C., Yao, F., Kartadikaria, A.R., Guo, D. and Hoteit, I., 2016. The Eddy Kinetic Energy Budget in the Red Sea. Journal of Geophysical Research: Oceans, 121(7), pp. 4732-4747. doi:10.1002/2015JC011589

20. Stepanov, D.V., 2018. Mesoscale Eddies and Baroclinic Instability over the Eastern Sakhalin Shelf of the Sea of Okhotsk: A Model-Based Analysis. Ocean Dynamics, 68(10), pp. 13531370. doi:10.1007/s 10236-018-1192-2

21. Kubryakov, A.A. and Stanichny, S.V., 2015. Mesoscale Eddies in the Black Sea from Satellite Altimetry Data. Oceanology, 55(1), pp. 56-67. doi:10.1134/S0001437015010105

22. Menna, M. and Poulain, P.-M., 2014. Geostrophic Currents and Kinetic Energies in the Black Sea Estimated from Merged Drifter and Satellite Altimetry Data. Ocean Science, 10(2), pp. 155-165. doi:10.5194/os-10-155-2014

23. Suvorov, A.M. and Shokurova, I.G., 2004. Annual and Interdecadal Variability of the Available Potential Energy in the Black Sea. Physical Oceanography, 14(2), pp. 84-95. doi:10.1023/BTOCE.0000037872.25674.ac

24. Cessi, P., Pinardi, N. and Lyubartsev, V., 2014. Energetics of Semienclosed Basins with Two-Layer Flows at the Strait. Journal of Physical Oceanography, 44(3), pp. 967-979. doi:10.1175/JPO-D-13-0129.1

25. Demyshev, S.G., 2004. Energy of the Black Sea Climatic Circulation. Part I: Discrete Equations of the Rate of Change of Kinetic and Potential Energy. Metjorologiya i Gidrologiya, (9), pp. 65-80 (in Russian).

26. Demyshev, S.G., 1998. Buoyant-Force Approximation in a Numerical Model of Baroclinic Ocean Currents. Izvjstiya, Atmosphjric and Ocjanic Physics, 34(3), pp. 362-369.

27. Lorenz, E.N., 1955. Available Potential Energy and the Maintenance of the General Circulation. Tjllus, 7(2), pp. 157-167. doi:10.3402/tellusa.v7i2.8796

28. Roquet, F., 2013. Dynamical Potential Energy: A New Approach to Ocean Energetics. Journal of Physical Ocjanography, 43(2), pp. 457-476. doi:10.1175/JPO-D-12-098.1

29. Demyshev, S.G. and Dymova, O.A., 2018. Numerical Analysis of the Black Sea Currents and Mesoscale Eddies in 2006 and 2011. Ocjan Dynamics, 68(10), pp. 1335-1352. doi:10.1007/s10236-018-1200-6

30. Altman, E.N. and Simonov, A.I., eds, 1991. [Hydromjtjorology and Hydrochjmistry of Sjas in thj USSR. Vol. IV. Black Sja. Issuj 1. Hydromjtjorological Conditions]. Saint-Petersburg: Gidrometeoizdat, 428 p. (in Russian).

31. Mellor, G.L. and Yamada, T., 1982. Development of a Turbulence Closure Model for Geophysical Fluid Problems. RjVijws of Gjophysics, 20(4), pp. 851-875. doi:10.1029/RG020i004p00851

About thj author:

Sergey G. Demyshev, Chief Researcher, Head of the Wave Theory Department, Marine

Hydrophysical Institute of RAS (2 Kapitanskaya Str., Sevastopol, 299011, Russian Federation), Ph.D.

(Math-Phys ), Scopus Author ID: 6603919865, SPIN-Kog: 1848-2350, IstinaResearcherID

(IRID): 17369115, ResearcherID: C-1729-2016, ORCID ID: 0000-0002-5405-2282,

demy shev@gmail. com

Thj author has rjad and approVjd thj final manuscript.

Thj author djclarjs that hj has no conflict of intjrjst.