Finite-Difference Approximation of the Potential Vorticity Equation for a Stratified Incompressible Fluid and an Example of its Application for Modeling the Black Sea Circulation Part II. Discrete Equation of Potential Vorticity in a Quasi-Static Approximation and an Example of its Application for Simulation

Full Text

See Part I: Demyshev, S.G., 2024. Finite-Difference Approximation of the Potential Vorticity Equation for a Stratified Incompressible Fluid and an Example of its Application for Modeling the Black Sea Circulation. Part I. Finite-Difference Equation of Potential Vorticity of Ideal Fluid. Physical Oceanography, 31(2), pp. 149-160.

Introduction

To study circulation in the atmosphere and ocean, analysis of potential vorticity (PV) is of fundamental importance, since PV characterizes the role of nonlinear processes in fluid dynamics. In the field of potential mass force in the absence of viscosity and density diffusion (ideal fluid), PV is an invariant [1] and therefore its structure determines the trajectory of fluid particles preserving potential vorticity. The difficulty is that in reality the abovementioned conditions are not met or only met in an approximate form, since friction, diffusion and diabaticity change the PV of seawater particles. When analyzing the equation of potential vorticity of a stratified fluid, it is possible to estimate the influence of nonlinear and diffusion factors on its evolution.The importance of Ertel’s theorem for research in the field of physical oceanography is discussed in [2]. It emphasizes that for large-scale water movements, the appropriate form of vorticity is potential vorticity, which includes such physically different elements as velocity vortex and seawater density. Ertel's theorem, which most other vorticity theorems are derived from, determines dynamic evolution of potential vorticity. In turn, the absence of potential vorticity means an inertial-gravitational mode of motion that depends on the ocean stratification. Since Ertel's theorem does not depend on the specific type of Lagrangian invariant, modified formulas for potential vorticity can be applied. As an example, the work [3] introducing “optimal” potential vorticity can be pointed out. As a result, the approach used permits to quantify the degree of disequilibrium in atmospheric climate processes.

Studies of the dynamics of currents in the atmosphere and ocean based on the potential vorticity analysis are few. Apparently, this is due to two factors. Firstly, potential vorticity according to Ertel is a kinematic quantity [4], from which it is impossible to determine intensity of the vortex structure of currents and even the rotation sign. Secondly, it is often sufficient to consider the potential vorticity in the quasi-genostrophic approximation according to Rossby [5] that is a dynamic characteristic [4] and contains the necessary information about the dynamics of currents.

For the analysis of atmospheric forecasts, the work [6] plays an important role. It analyses the possibilities of using isentropic maps of potential vorticity to represent some dynamic processes in the atmosphere. Examples from operational weather analysis and from idealized theoretical models are given to illustrate this approach and its relationship to classical synoptic concepts. The structure, reasons for the formation and stability of cyclones and blocking anticyclones, physical mechanisms of Rossby wave propagation, baroclinic and barotropic instability in space and time are discussed.

In [4], the concept of “potential vorticity” is analyzed and the basic relationships for calculations are considered, the approaches of Rossby and Ertel are studied. As an example, estimates are given from observational data of potential vorticity for the quasi-permanent anticyclonic Lofoten eddy in the Norwegian Sea. The authors showed that Ertel potential vorticity is a kinematic characteristic and Rossby potential vorticity in the quasi-geostrophic approximation is dynamic.

Analysis of the Ertel potential vorticity equation permits to estimate the contribution of nonlinear and diffusion effects to the balance of forces that determine PV evolution. This study is a continuation of the research [7]. It is aimed at obtaining a discrete equation for the potential vorticity of a stratified fluid in a quasi-static approximation as a consequence of the initial finite-difference system of equations of the Black Sea dynamics model using the example of the Black Sea circulation with realistic atmospheric conditions for 2011 and to analyze the resulting potential vorticity equation.

Equation of potential vorticity of a stratified incompressible fluid in the quasi-static approximation

In the Boussinesq approximation and quasi-statics in the Cartesian coordinate system, the fluid motion in the domain $\Omega$ with boundary $\partial \Omega$ in the Gromeka-Lamb form is described by the following system of equations:

$\frac{\partial {\overrightarrow{U}}_H}{\partial t}+\overrightarrow{\xi }x{\overrightarrow{U}}_H=-\frac{1}{{\rho }_0}\nabla {(P+E)}_{}+\overrightarrow{F},$ (1)

$\frac{\partial P}{\partial z}=g\rho ,$ (2)

$\nabla \overrightarrow{U}=0$, (3)

$\frac{\partial T}{\partial t}+\text{div}\left(T\overrightarrow{U}\right)={\left({\kappa }^T{T}_z\right)}_z-{\kappa }^H{\nabla }^{4}T$, (4)

$\frac{dS}{dt}+\text{div}\left(S\overrightarrow{U}\right)={\left({\kappa }^S{S}_z\right)}_z-{\kappa }^H{\nabla }^{4}S$, (5)

$\rho =\Phi (T,S)$. (6)

The following designations are introduced: $\overrightarrow{U}=({\overrightarrow{U}}_H,w)=(u,v,w)$ are the components of the current velocity vector along the axes (x, y, z) directed to the east, north and vertically downwards, respectively; $\overrightarrow{F}=(F^u,F^v)$, $\overrightarrow{g}$=(0, 0, g) is free fall acceleration; (T, S, P and p) are temperature, salinity, pressure and sea water density; p₀ = 1 g/cm³ (here and henceforth it is pressure and density normalized to p₀ ); $\overrightarrow{f}$=(0, 0, $f^z$) is the Coriolis parameter, where $f^z$ = 2$\omega$ sin $\phi$; $\omega$ is angular velocity of the Earth’s rotation; $\phi$ is latitude.

In equation (1), taking into account the quasi-static approximation, the absolute velocity vortex and kinetic energy of motion are introduced:

$\overrightarrow{\xi }=({\xi }^x,{\xi }^y,{\xi }^z),\text{ где }{\xi }^x=-\frac{\partial v}{\partial z},{\xi }^y=\frac{\partial u}{\partial z},{\xi }^z=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}+f^z={\xi }^r+f^z,$

$F^u={\left({\mu }^{\text{ver}}{u}_z\right)}_z-{\mu }^{\text{hor}}{\nabla }^{4}u,F^v={\left({\mu }^{\text{ver}}{v}_z\right)}_z-{\mu }^{\text{hor}}{\nabla }^{4}v,\overrightarrow{F}=(F^u,F^v),$

$E=\rho {}_0\frac{u^2+v^2}{2},$

where ${\mu }^{\text{ver}},{\mu }^{\text{hor}}$ are coefficients of vertical and horizontal exchange of momentum.

On the surface at z = 0

${\nu }_V{u}_z=-{\tau }^x,{\nu }_V{v}_z=-{\tau }^y,$ $w=-{\varsigma }_t,$ ${\kappa }^T{T}_z=Q^T,$ ${\kappa }^S{S}_z=\frac{Ev-\mathrm{Pr}}{{\rho }_1}{S}_0;$ (7)

at the bottom at z = H(x, y)

u = v = w = 0, T_z = S_z = 0. (8)

The following designations are used: (${\tau }^x,{\tau }^y)$ is tangential wind stress; $Q^T$ is heat flow; Ev is sea water evaporation; Pr is precipitation; S₀ is model salinity on the sea surface; ${\rho }_1$ is sea water density on the sea surface.

Functions ${\mu }^{ver},{\kappa }^T,{\kappa }^S$ were calculated in accordance with the Mellor–Yamada parameterization [6].

On solid side walls for meridional boundary sections:

$u={\nabla }^{2}u=v_x={\nabla }^2{v}_x=\mathrm{0,}{T}_x={\left({\nabla }^{2}T\right)}_x=S_x={\left({\nabla }^{2}S\right)}_x=\mathrm{0,}$ (9)

for zonal boundary sections:

$v={\nabla }^{2}v=u_y={\nabla }^2{u}_y=\mathrm{0,}{T}_y={\left({\nabla }^{2}T\right)}_y=S_y={\left({\nabla }^{2}S\right)}_y=0.$ (10)

In boundary areas where water flows in, the following conditions are used:

for meridional sections

$u=u^p,{\nabla }^{2}u=v_x={\nabla }^2{v}_x=\mathrm{0,}T=T^p,S=S^p,{\left({\nabla }^{2}T\right)}_x={\left({\nabla }^{2}S\right)}_x=\mathrm{0,}$ (11)

for zonal sections

$v=v^p,{\nabla }^{2}v=u_y={\nabla }^2{u}_y=\mathrm{0,}T=T^p,S=S^p,{\left({\nabla }^{2}T\right)}_y={\left({\nabla }^{2}S\right)}_y=0.$ (12)

For the Upper Bosphorus Current and for the Kerch Strait, when the current is directed from the Black Sea to the Sea of Azov

$v=v^s,{\nabla }^{2}v=u_y={\nabla }^2{u}_y=\mathrm{0,}{T}_y=\mathrm{0,}{S}_y=\mathrm{0,}{\left({\nabla }^{2}T\right)}_y={\left({\nabla }^{2}S\right)}_y=0.$ (13)

When $t=t^0$, the following initial conditions are set:

$u=u^0(x,y,z),v=v^0(x,y,z),\varsigma ={\varsigma }^0(x,y),T=T^0(x,y,z),S=S^0(x,y,z).$ (14)

Based on the system of equations (1)–(6), the Ertel equation is derived. Applying the corresponding operation from relations (1)–(2) and taking into account the continuity equation (3), the equation for $\overrightarrow{\xi }$ is obtained

$\frac{\partial \overrightarrow{\xi }}{\partial t}+\overrightarrow{U}(\nabla \overrightarrow{\xi })-\overrightarrow{\xi }\nabla \overrightarrow{U}=\nabla \times (\overrightarrow{g}\rho +{\overrightarrow{F}}^{\xi }),$ (15)

where ${\overrightarrow{F}}^{\xi }=(F^u,F^v\mathrm{,0}).$

A consequence of equations (4), (5) is the equation for density

$\frac{\partial \rho }{\partial t}+\text{div}\left(\overrightarrow{U}\rho \right)=\Phi (D^T,D^S),$ (16)

where the following designations are introduced

$D^T={\left({\kappa }^T{T}_z\right)}_z-{\kappa }^H{\nabla }^{4}T,D^S={\left({\kappa }^S{S}_z\right)}_z-{\kappa }^H{\nabla }^{4}S.$ (17)

The potential vorticity of an incompressible fluid in the quasi-static approximation has the following form:

$\omega =\overrightarrow{\xi }\nabla \rho =-\frac{\partial v}{\partial z}\frac{\partial \rho }{\partial x}+\frac{\partial u}{\partial z}\frac{\partial \rho }{\partial y}+\left[\right(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y})+f^z]\frac{\partial \rho }{\partial z}.$ (18)

Then, from equations (15), (16), taking into account expressions (3), (17) and (18), the Ertel equation follows in the quasi-static approximation for a viscous fluid

$\frac{\partial \omega }{\partial t}+\text{\hspace{0.33em}div}\left(\overrightarrow{U}\omega \right)={\text{\hspace{0.33em}Ф}}^{\omega },$ (19)

where

${\text{Ф}}^{\omega }=(\overrightarrow{\xi }\nabla )\Phi +\left[\right(\nabla \times \overrightarrow{F^{\xi }})\nabla ]\rho .$ (20)

Discrete potential vorticity equation in a quasi-static approximation

In accordance with the difference operators introduced in [7], the following differential-difference equations of model (1)–(6) (differential in time) are written out:

$\frac{d{u}_{i+1/\mathrm{2,}j,k}}{dt}-{\left[v,{\xi }^z\right]}_{i+1/\mathrm{2,}j,k}+{\left[w,{\xi }^y\right]}_{i+1/\mathrm{2,}j,k}=$$-{\delta }_x(E_{i+1/\mathrm{2,}j/k}+P_{i+1/\mathrm{2,}j,k})+F_{i+1/\mathrm{2,}j,k}^u,$ (21)

$\frac{d{v}_{i,j+1/\mathrm{2,}k}}{dt}+{\left[u,{\xi }^z\right]}_{i,j+1/\mathrm{2,}k}-{\left[w,{\xi }^x\right]}_{i,j+1/\mathrm{2,}k}=$$-{\delta }_y(E_{i,j+1/\mathrm{2,}k}+P_{i,j+1/\mathrm{2,}k})+F_{i,j+1/\mathrm{2,}k}^v,$ (22)

${\delta }_z{P}_{i,j,k+1/2}=g{\rho }_{i,j,k+1/2},$ (23)

${\delta }_x{u}_{i,j,k}^{}+{\delta }_y{v}_{i,j,k}^{}+{\delta }_z{w}_{i,j,k}^{}=\mathrm{0,}$ (24)

$\frac{d{T}_{i,j,k}}{dt}+{\delta }_x\left(u_{i,j,k}^{}{T}_{i,j,k}^{}\right)+{\delta }_y\left(v_{i,j,k}^{}{T}_{i,j,k}^{}\right)+{\delta }_z\left(w_{i,j,k}^{}{T}_{i,j,k}^{}\right)=$

$={\delta }_z\left({\kappa }_{i,j,k}^T{\delta }_z{T}_{i,j,k}\right)-{\kappa }^H{\nabla }_{xy}^2\left({\nabla }_{xy}^2{T}_{i,j,k}\right),$ (25)

$\frac{d{S}_{i,j,k}}{dt}+{\delta }_x\left(u_{i,j,k}^{}{S}_{i,j,k}^{}\right)+{\delta }_y\left(v_{i,j,k}^{}{S}_{i,j,k}^{}\right)+{\delta }_z\left(w_{i,j,k}^{}{S}_{i,j,k}^{}\right)=$

$={\delta }_z\left({\kappa }_{i,j,k}^S{\delta }_z{S}_{i,j,k}\right)-{\kappa }^H{\nabla }_{xy}^2\left({\nabla }_{xy}^2{S}_{i,j,k}\right),$ (26)

${\rho }_{i,j,k}^{}={\rho }_0+{\alpha }_1^T{T}_{i,j,k}+{\alpha }_1^S{S}_{i,j,k}+{\alpha }_2^T{T}_{i,j,k}^2+{\alpha }_{}^{ST}{S}_{i,j,k}{T}_{i,j,k}.$ (27)

In the quasi-static approximation, the components of the velocity vortex (Fig. 1) have the following form:

${\xi }_{i,j+1/\mathrm{2,}k+1/2}^x=-{\delta }_z{v}_{i,j+1/\mathrm{2,}k+1/2},$${\xi }_{i+1/\mathrm{2,}j,k+1/2}^y={\delta }_z{u}_{i+1/\mathrm{2,}j,k+1/2},$

${\xi }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k}^z={\delta }_x{v}_{i+1/\mathrm{2,}j+1/\mathrm{2,}k}-{\delta }_y{u}_{i+1/\mathrm{2,}j+1/2}+f_{j+1/2}^z.$ (28)

F i g. 1. Distribution of variables in box (i, j, k). At the box vertices indicated by asterisks, PV (ω) is determined and on its edges – the components of absolute vortex velocity

Representation (28) follows that at the vertices of the box () (Fig. 1) an important relation is satisfied

${\delta }_x{\xi }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^x+{\delta }_y{\xi }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^y+{\delta }_z{\xi }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^z=0.$

We believe that the terms in square brackets on the left side of equations (21)–(22) are written in a form that ensures conservation of enstrophy and energy [9] for the shallow water approximation and correspond to formulas (32) in [10].

Equations for the components of the absolute velocity vortex (analogous to the equation (15)) – for ${\xi }^x$ at the point $i,j+1/\mathrm{2,}k+1/2$, for ${\xi }^y$ at $i+1/\mathrm{2,}\text{\hspace{0.33em}j, }k+1/2$ and for ${\xi }^z$ at $i+1/\mathrm{2,}j+1/\mathrm{2,}k$ – taking into account viscosity, follow from relations (21)–(24) and have the following form:

$\frac{d{\xi }_{}^x}{dt}+{\delta }_z\left(\right[w,{\xi }_{}^x\left]\right)-{\delta }_z\left(\right[u,{\xi }_{}^z\left]\right)=g{\delta }_y{\overline{\rho }}^z-{\delta }_z{F}^v,$

$\frac{d{\xi }_{}^y}{dt}+{\delta }_z\left(\right[w,{\xi }_{}^y\left]\right)-{\delta }_z\left(\right[v,{\xi }_{}^z\left]\right)=-g{\delta }_x{\overline{\rho }}^z+{\delta }_z{F}^u,$ (29)

$\frac{d{\xi }_{}^z}{dt}+{\delta }_x\left(\right[u,{\xi }_{}^z\left]\right)+{\delta }_y\left(\right[v,{\xi }_{}^z\left]\right)-{\delta }_x\left(\right[w,{\xi }_{}^x\left]\right)-{\delta }_y\left(\right[w,{\xi }_{}^y\left]\right)={\delta }_x{F}^v-{\delta }_y{F}^u.$

The following designations are introduced:

${\Lambda }_{i,j+1/\mathrm{2,}k+1/2}^x=-{\delta }_z{F}_{i,j+1/\mathrm{2,}k+1/2}^v,$

${\Lambda }_{i+1/\mathrm{2,}j,k+1/2}^y={\delta }_z{F}_{i+1/\mathrm{2,}j,k+1/2}^u,$

${\Lambda }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k}^z={\delta }_x{F}_{i,j+1/\mathrm{2,}k+1/2}^v-{\delta }_y{F}_{i+1/\mathrm{2,}j+1/\mathrm{2,}k}^u.$

The density equation at the point (i, j, k) is a consequence of relations (25)–(27) and is written as follows:

$\frac{d\rho }{dt}+{\delta }_x\left(u\rho \right)+{\delta }_y\left(v\rho \right)+{\delta }_z\left(w\rho \right)=D_V^{\rho }+{\kappa }^H{D}_H^{\rho }=D^{\rho },$ (30)

where

$\begin{array}{l}{D}_V^{\rho }={\alpha }_1^T{\delta }_z\left[{\kappa }^V\right({\delta }_{z}T\left)\right]+{\alpha }_1^S{\delta }_z\left[{\kappa }^V\right({\delta }_{z}S\left)\right]+2{\alpha }_2^T\left[{\delta }_z\right[{\kappa }^V{\overline{T^{}}}^z\left({\delta }_{z}T\right)]+\\ +{\alpha }^{TS}\left[T{\delta }_z\left[{\kappa }^V\right({\delta }_{z}S\left)\right]+S{\delta }_z\left[{\kappa }^V\right({\delta }_{z}T)\right.]\left.\right],\end{array}$

$D_H^{\rho }={\alpha }_1^T{\nabla }_{xy}^{2}T+{\alpha }_1^S{\nabla }_{xy}^{2}S+2{\alpha }_1^{T}T{\nabla }_{xy}^{2}T+{\alpha }_{}^{TS}[T{\nabla }_{xy}^{2}S+S{\nabla }_{xy}^{2}T].$

Carrying out similar [5] transformations with equations (29), (30), the equation of potential vorticity is obtained at the point $i+1/\mathrm{2,} j+1/2,k+1/2$ in the quasi-static approximation

$\frac{d\omega }{dt}+{\delta }_x({\Upsilon }^x{\overline{\rho }}^{yz}+{\xi }^x{R}^x)+{\delta }_y({\Upsilon }^y{\overline{\rho }}^{xz}+{\xi }^y{R}^y)+{\delta }_z({\Upsilon }^z{\overline{\rho }}^{xy}+{\xi }^z{R}^z)={\overline{{\Lambda }^x}}^x{\delta }_x\left({\overline{\rho }}^{yz}\right)+$

$+{\overline{{\Lambda }^y}}^y{\delta }_y\left({\overline{\rho }}^{xz}\right)+{\overline{{\Lambda }^z}}^z{\delta }_z\left({\overline{\rho }}^{xy}\right)+{\overline{{\xi }^x}}^x{\delta }_x\left({\overline{D^{\rho }}}^{yz}\right)+{\overline{{\xi }^y}}^y{\delta }_y\left({\overline{D^{\rho }}}^{xz}\right)+{\overline{{\xi }^z}}^z{\delta }_z\left({\overline{D^{\rho }}}^{xy}\right),$ (31)

where the right side of equation (31) is a difference analog of expression (20). Taking into account the quasi-static approximation and equations (29), the following notations are adopted:

${\Upsilon }_{i,j+1/\mathrm{2,}k+1/2}^x={\delta }_z{\left(\right[w,{\xi }_{}^x\left]\right)}_{i,j+1/\mathrm{2,}k}-{\delta }_z{\left(\right[u,{\xi }_{}^z\left]\right)}_{i,j+1/\mathrm{2,}k},$

${\Upsilon }_{i+1/\mathrm{2,}j,k+1/2}^y={\delta }_z{\left(\right[w,{\xi }_{}^y\left]\right)}_{i+1/\mathrm{2,}j,k+1/2}-{\delta }_z{\left(\right[v,{\xi }_{}^z\left]\right)}_{i+1/\mathrm{2,}j,k+1/2},$

$\begin{array}{l}{\Upsilon }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k}^z={\delta }_x\left(\left[u,{\xi }^z\right]\right)+{\delta }_y{\left(\left[v,{\xi }^z\right]\right)}_{i+1/\mathrm{2,}j+1/\mathrm{2,}k}-\\ -{\delta }_x{\left(\left[w,{\xi }^x\right]\right)}_{i+1/\mathrm{2,}j+1/\mathrm{2,}k}-{\delta }_y{\left(\left[w,{\xi }^y\right]\right)}_{i+1/\mathrm{2,}j+1/\mathrm{2,}k},\end{array}$

${\overline{{\overline{\rho }}_{i,j+1/\mathrm{2,}k+1/2}^{yz}}}^x={\overline{{\overline{\rho }}_{i+1/\mathrm{2,}j,k+1/2}^{xz}}}^y={\overline{{\overline{\rho }}_{i+1/\mathrm{2,}j+1/\mathrm{2,}k}^{xy}}}^z={\overline{{\rho }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}}}^{xyz},$

$\begin{array}{l}{\overline{R_{i,j+1/\mathrm{2,}k+1/2}^x}}^x={\overline{R_{i+1/\mathrm{2,}j,k+1/2}^y}}^y={\overline{R_{i+1/\mathrm{2,}j+1/\mathrm{2,}k}^z}}^z={\overline{{\delta }_x\left(u_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^{}{\rho }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^{}\right)+}}^{xyz}\\ +{\overline{{\delta }_y\left(v_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^{}{\rho }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^{}\right)+{\delta }_z\left(w_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^{}{\rho }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^{}\right)}}^{xyz}.\end{array}$

The difference analogue of the Ertel potential vorticity has the following form:

${\varpi }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}={\overline{{\xi }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^x}}^x\delta {}_x\overline{\rho }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^{yz}+{\overline{{\xi }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^y}}^y\delta {}_y\overline{\rho }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^{xz}+$

$+{\overline{{\xi }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^z}}^z\delta {}_z\overline{\rho }_{i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2}^{xy}={\varpi }_{}^x+{\varpi }_{}^y+{\varpi }_{}^z,$ (32)

where the designations ${\varpi }_{}^x,{\varpi }_{}^y,{\varpi }_{}^z$ are obvious.

The difference between equation (31) and equation (45) in [7] is not only that viscosity and diffusion are taken into account, but also that the components of the absolute velocity vortex have the form (28).

Numerical analysis of components of the potential vorticity equation based on calculation results of circulation with atmospheric conditions for 2011.

In numerical forecasting experiments, the following parameters were set. Calculations were carried out with a uniform step along horizontal coordinates of 1.6 km; 27 horizons were used vertically with condensation in the upper sea layer. Accounting for the runoff of the Black Sea rivers and flows through the Bosporus and Kerch Straits (conditions (9)–(13)) corresponded to the data in [11]. The water temperature at the mouths of rivers (conditions (11)–(12)), except for the rivers of Turkey, was set from [11]. It was assumed that the temperature of the Turkish rivers is equal to the temperature of the coastal sea waters. In the Upper Bosphorus Current, the temperature and salinity were assumed to be the same as in the sea, in accordance with conditions (11). In the Lower Bosporus stream, the salinity was taken to be 35‰ with the temperature of 16 °С.

SKIRON data for 2011 [12] were used to set the atmospheric effect in equations (7), (8); based on the Mellor–Yamada theory [8], vertical mixing was described. The initial conditions (14) in this calculation corresponded to 1 January 2011. The calculation was carried out for a year of model time, its parameters and results are described in detail in [13].

As an example, the PV is considered (formula (32)) for two points in time – winter (Fig. 2, a) and summer (Fig. 2, b) periods, when the circulation structure is noticeably different.

F i g. 2. Specified sea level on 1 February (a) and 1 August (b) 2011

As of 1 January 2011, the specified sea level was an extensive cyclonic gyre with two synoptic eddies – the Sevastopol and the southwestern anticyclones. In contrast to the winter circulation, the cyclonic circulation splits into two in the summer season (Fig. 2, b). In the western part of the basin along the deep slope, the Black Sea Rim Current spreads in the form of a narrow jet current (Fig. 2, a). In winter, powerful meanders form along the northern periphery of the northwestern shelf and an intense anticyclone forms in the southeastern corner of the sea.

Winter circulation corresponds to the potential vorticity shown in Fig. 3, a – d. There are qualitative differences in its structure in depth. In the upper 30-meter layer, large PV values are concentrated in two areas (Fig. 3, a, b). The first is the northwestern part of the sea, limited approximately by coordinates 44°–46°N, 29°–31°E, where water dynamics is largely determined by the river runoff, primarily the Danube. Therefore, there are large spatial gradients in the density field, which determinevalues in this sea area. The second area is a deep slope, where the Anatolian coast and the southern periphery of the northwestern shelf stand out (Fig. 3, b).

F i g. 3. PV at horizons 3.75 (a), 35 (b), 106.25 (c) and 350 m (d) on 1 February 2011

In the area of the Sevastopol and the southwestern vortex, no extreme values of are observed. The reason is that it is a scalar quantity equal to the product of density gradient and absolute velocity vortex and, therefore, its large value does not necessarily mean an increase in vorticity and, moreover, its sign does not determine the vortex rotation sign [2]. Below the upper 50-meter layer, the highest values of in the coastal strip are observed (Fig. 3, c, d). Local maxima are concentrated in relatively small zones (~ 10 km), which are clearly visible at a depth of 100 m (Fig. 3, c). In the central part of the sea, the PV structure is quite homogeneous in space.

An illustration for the PV analysis is its calculation for summer period, when the circulation is less regular and its vortex structure is more pronounced (see Fig. 2, b). Two features appear in the structure of potential vorticity in August 2011. The band of small values $\varpi$ in the upper layer (Fig. 4, a, b) corresponds to the region of greatest mixing along the vertical density field, that is, the value $\delta {}_z\overline{\rho }_{}^{xy}$ is very small – two to three orders of magnitude less than vertical density gradients of the surrounding water. The second feature is a homogeneous structure and low PV values in areas approximately corresponding to the cores of the southeastern anticyclone, southwestern and eastern gyres. These features are also determined by the structure $\delta {}_z\overline{\rho }_{}^{xy}$, which varies slightly over space. This type of potential vorticity in central parts of gyres is consistent with the conclusions of the work [4], where the PV reconstructed from observational data in the Lofoten gyre area has a similar structure. On the lower horizons (Fig. 4, c, d) along the area boundary, there is a narrow band of heterogeneous values due to differences in the bottom topography; a small spatial variability of this value is observed in the central part.

РF i g. 4. PV at horizons 3.75 (a), 35 (b), 106.25 (c) and 350 m (d) on 1 August 2011

The main contribution to the vortex structure, as a rule, is made by the component ${\varpi }_{}^z$ [4]. Its magnitude is determined by the quasi-geostrophic nature of the movement and the vertical stratification of sea water. In areas where there is an inflow of fresh (river mouths) or salt (straits) waters, horizontal components of potential vorticity may be of predominant importance. As an example, Fig. 5 shows the values at horizons of 3.75 and 106.25 m.

Comparison of Fig. 5, a, c, e and 3, a shows that horizontal component ${\varpi }_{}^z$ (Fig. 5, a) makes the main contribution in the zone of the Danube fresh water inflow in the northwestern part of the sea, in the rest of the area ${\varpi }_{}^z$ determines vortex structure (Fig. 5, e).

F i g. 5. PV components: ${\varpi }_{}^x$ at horizons 3.75 (a) and 106.25 m (b); ${\varpi }_{}^y$ at horizons 3.75 (c) and 106.25 m (d); ${\varpi }_{}^z$ at horizons 3 (e) and 106.25 m (f) on 1 February 2011

At the horizon 106.25 m (Fig. 5, b, d, f), the vertical component is larger than ${\varpi }^x$ (Fig. 5, b) and ${\varpi }^y$ (Fig. 5, d), so it (Fig. 5, f) quite accurately determines potential vorticity type at the horizon 106.25 m (see Fig. 3, c).

Direct calculations establish that in the upper layer of the sea the form of the component ${\varpi }_{}^z$ corresponds qualitatively to ${\delta }_z{\overline{\rho }}^{xy}$, but its absolute value is several orders of magnitude smaller than ${\overline{{\xi }^z}}^z$. In turn, the absolute vortex structure is quite homogeneous and positive, therefore, when ${\delta }_z{\overline{\rho }}^{xy}$ is multiplied by the potential vorticity structure ${\overline{{\xi }^z}}^z$, it characterizes ${\delta }_z{\overline{\rho }}^{xy}$ and its quantitative value depends on ${\overline{{\xi }^z}}^z$. The value ${\xi }^z$ is determined by two terms – relative and planetary vorticity. If we evaluate $f^z$ contribution to the absolute vortex, $f^z$ is comparable in magnitude to the relative vortex and increases the values of ${\xi }^r$. On average, the ${\xi }^z$ magnitude is two orders of magnitude larger than ${\delta }_z{\overline{\rho }}^{xy}$. During this season, as a result of winter convection, the factor due to the vertical density gradient in the upper layer of the sea is small, with the exception of the river runoff area, where its value can be significant.

At the horizon of 106.25 m (Fig. 5, f), both factors are positive and ${\overline{{\xi }^z}}^z$, on average, is smaller by several orders of magnitude. PV variability is observed in the coastal area; PV is homogeneous in the central part of the sea. It should be noted that, firstly, at the lower horizons (approximately below 50 m depth) the relative vortex is in absolute value smaller than $f^z$. Secondly, since the integral over the horizontal surface of ${\xi }^r$, whose difference from zero is determined by the river runoff and the water exchange through the straits, is small, then the ${\xi }^r$ structure contains zones of cyclonic and anticyclonic rotation of waters. At the same time, planetary vorticity is positive and greater than ${\xi }^r$ and, therefore, it determines quantitative values of PV with corrections introduced by the relative vortex to the qualitative structure of potential vorticity at deep horizons.

Let us consider the contribution of nonlinear forces to evolution of $\varpi$. The following designations are introduced:

$C^x={\delta }_x\left({\Upsilon }^x{\overline{\rho }}^{yz}\right)+{\delta }_x\left({\xi }^x{R}^x\right)=C_1^x+C_2^x,C^y={\delta }_y\left({\Upsilon }^y{\overline{\rho }}^{xz}\right)+{\delta }_y\left({\xi }^y{R}^y\right)=C_1^y+C_2^y,$

$C^z={\delta }_z\left({\Upsilon }^z{\overline{\rho }}^{xy}\right)+{\delta }_z\left({\xi }^z{R}^z\right)=C_1^z+C{\text{}}_2^z,C^S=C^x+C^y+C^z.$

The main contribution to the temporal evolution of PV is made by nonlinear forces in the upper layer in the coastal sea area (Fig. 6, a, b). Their contribution is different for various areas: it is greater in the northwestern part (Fig. 6, a) and along the Anatolian coast (Fig. 6, b). The estimates show that their quantitative differences in absolute value between the central part of the sea and its periphery are several orders of magnitude. At the lower horizons (Fig. 6, c, d), the largest values of the nonlinear terms in Ertel’s equation are concentrated in the form of a narrow alongshore strip with a more pronounced character near the southern coast.

Let us consider the contribution of individual terms $C^x,C^y,C^z$ into $C^S$.

F i g. 6. $C^S$ at horizons 3 (a) and 106.25 m (b) on 1 February 2011, and at horizons 3.75 (c) and 106.25 m (d) on 1 August 2011

In the upper layer (Fig. 7, a, c, e), areas of large absolute values $C^x,C^y,C^z$ have a similar structure. In the southeastern corner of the basin and in the northeastern part of the sea, limited by coordinates 42°– 44°N, 37°–39°E etc., areas of $C^x,C^y,C^z$ values close to zero are observed. The calculated average and maximum $C^x,C^y,C^z$ values (Fig. 7, a, c, e) in comparison with $C^S$ (Fig. 6, a) indicate that the extreme values differ by an order of magnitude, the average values – by two orders of magnitude. This means that $C^x,C^y,C^z$ are mutually compensated and the result is the structure shown in Fig. 6, a. The direct calculations establish that the main contribution to the nonlinear terms $C^x,C^y,C^z$ in the upper layer is made by $C_1^x,C_1^y,C_1^z$ accordingly, that is ${\delta }_x\left({\Upsilon }^x{\overline{\rho }}^{yz}\right),{\delta }_y\left({\Upsilon }^y{\overline{\rho }}^{xz}\right),{\delta }_z\left({\Upsilon }^z{\overline{\rho }}^{xy}\right).$ The estimate of the order of magnitude shows that $C^S$ is on average two orders of magnitude less than each of the terms $C^x,C^y,C^z$.

F i g. 7. Components: at horizons 3.75 (a) and 106.25 m (b), at horizons 3.75 (c) and 106.25 m (d) and at horizons 3.75 (e) and 106.25 m (f) on 1 February 2011

A similar situation occurs for the calculated fields as of 1 August 2011.

Conclusion

For a system of discrete equations of sea dynamics in the hydrostatic approximation and taking into account viscosity, diffusion, a finite-difference equation for the potential vorticity of a stratified incompressible fluid is obtained. Just as in the more general case, it has a divergent appearance and differs from its differential counterpart. Because a nonlinear state is used to calculate the density, the resulting discrete equation for PV is not an exact consequence of the finite-difference equations of the model. Additional research to estimate the influence of the nonlinear nature of the equation of state on the results obtained is needed.

Analysis of the magnitude of potential vorticity itself confirmed previously obtained results that its vertical component is the main one. Horizontal components make a noticeable contribution in river runoff areas, water exchange through straits and in zones of sharp density field gradients. The qualitative appearance of PV in the upper layer of the sea has similar features recorded from observational data. The homogeneous structure for the central part of the vortex formations and the intense nature in the area of large transverse gradients in the density field in the frontal zones determine potential vorticity structure. In the deep layers of the sea, its highest values are concentrated in the form of a narrow coastal strip; in the rest of the sea, PV values are small.

Calculation of ${\varpi }^z$ components in winter on the upper horizons showed that in the upper layers of the sea, except coastal zones of river runoff, ${\xi }^z$ is determined, which is the sum of the relative vortex and $f^z$ (value of approximately the same order). In the lower layers of the sea, the quantitative values of PV are determined to a greater extent by planetary vorticity and its qualitative features – by the relative vortex structure.

From the analysis of nonlinear terms in the PV equation, it follows that in the upper layer of the sea the main contribution to the advection of potential vorticity $C^S$ is made in the northwestern part and along the Anatolian coast. At lower horizons, the highest $C^S$ values are observed along the coastal strip with a more pronounced character near the southern coast of the sea, which corresponds to the PV structure.

The calculation of $C^x,C^y,C^z$ terms for winter and summer periods enabled to establish two facts. Firstly, the magnitude of each $C^x,C^y,C^z$ is determined by ${\delta }_x\left({\Upsilon }^x{\overline{\rho }}^{yz}\right),{\delta }_y\left({\Upsilon }^y{\overline{\rho }}^{xz}\right),{\delta }_z\left({\Upsilon }^z{\overline{\rho }}^{xy}\right),$ that is, by the divergence from the product of nonlinear terms in the equations of motion and density. Secondly, $C^S$ is one and a half to two orders of magnitude less than each of $C^x,C^y,C^z$ components, that is, locally, the total sum of the vertical and horizontal advection of potential vorticity is two orders of magnitude less than each separately. A possible explanation for this result is as follows. Assume that the finite-difference analogues of the nonlinear terms in the PV equation are close to the differential form $\text{div}\left(\overrightarrow{U}\omega \right).$ Then, representing $\omega ={\omega }^S+{\varpi }^{*},$ where ${\omega }^S$ is a quantity averaged over space, at each point of the domain we have $\text{div}\left(\overrightarrow{U}{\omega }^S\right)={\omega }^S\text{div}\left(\overrightarrow{U}\right)=0$ or close to zero. Since the time variability of potential vorticity depends predominantly on $\text{div}\left(\overrightarrow{U}{\omega }^{*}\right)$, mutual compensation of nonlinear components along x, y, z occurs when calculating PV advection.

The general nature of the results obtained is a matter for further research. 

About the authors

Sergey G. Demyshev

Author for correspondence.
Email: demyshev@gmail.com
ORCID iD: 0000-0002-5405-2282

Head of Wave Theory Department, Chief Research Associate, DSc (Phys.-Math.)

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. F i g. 1. Distribution of variables in box (i, j, k). At the box vertices indicated by asterisks, PV (ω) is determined and on its edges – the components of absolute vortex velocity

Download (390KB)

Indexing metadata

3. F i g. 2. Specified sea level on 1 February (a) and 1 August (b) 2011

Download (254KB)

Indexing metadata

4. F i g. 3. PV at horizons 3.75 (a), 35 (b), 106.25 (c) and 350 m (d) on 1 February 2011

Download (633KB)

Indexing metadata

5. F i g. 4. PV at horizons 3.75 (a), 35 (b), 106.25 (c) and 350 m (d) on 1 August 2011

Download (604KB)

Indexing metadata

6. F i g. 5. PV components: xϖ at horizons 3.75 (a) and 106.25 m (b); yϖ at horizons 3.75 (c) and 106.25 m (d); zϖ at horizons 3 (e) and 106.25 m (f) on 1 February 2011

Download (1MB)

Indexing metadata

7. F i g. 6. SC at horizons 3 (a) and 106.25 m (b) on 1 February 2011, and at horizons 3.75 (c) and 106.25 m (d) on 1 August 2011

Download (1MB)

Indexing metadata

8. F i g. 7. Components: at horizons 3.75 (a) and 106.25 m (b), at horizons 3.75 (c) and 106.25 m (d) and at horizons 3.75 (e) and 106.25 m (f) on 1 February 2011

Download (1MB)

Indexing metadata