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 the Black Sea circulation in 2011

Full Text

Часть II. Дискретное уравнение потенциальной завихренности в квазистатическом приближении и пример его использования для расчета циркуляции Черного моря 2011 года¹

Введение

Для исследования циркуляции в атмосфере и океане анализ потенциальной завихренности (англ. potential vorticity, сокр. PV) имеет принципиальное значение, так как PV характеризует роль нелинейных процессов в динамике жидкости. В поле потенциальной массовой силы при отсутствии вязкости и диффузии плотности (идеальная жидкость) PV является инвариантом [1], и поэтому его структура определяет траекторию движения частиц жидкости, сохраняющих потенциальную завихренность. Трудность заключается в том, что в реальности вышеперечисленные условия не выполняются или выполняются в приближенном виде, поскольку трение, диффузия и диабатичность меняют PV частиц морской воды. При анализе уравнения потенциальной завихренности стратифицированной жидкости можно оценить влияние нелинейных и диффузионных факторов на ее эволюцию.

Важности теоремы Эртеля для исследований в области физической океанографии посвящена работа [2]. В ней подчеркивается, что для крупномасштабных движений воды подходящей формой завихренности является именно потенциальная завихренность, которая включает в себя такие физически разные элементы, как вихрь скорости и плотность морской воды. Теорема Эртеля, из которой следует большинство других теорем о завихренности, определяет динамическую эволюцию потенциальной завихренности. В свою очередь, отсутствие потенциальной завихренности означает инерционно-гравитационный режим движения, который зависит от стратификации океана. Поскольку теорема Эртеля не зависит от конкретного вида лагранжева инварианта, то могут использоваться модифицированные формулы потенциальной завихренности. В качестве примера можно указать на работу [3], в которой вводится «оптимальная» потенциальная завихренность. В результате используемый подход позволяет количественно оценить степень неравновесности атмосферных климатических процессов.

Исследования динамики течений в атмосфере и океане на основе анализа потенциальной завихренности немногочисленны. Видимо, это связано с двумя факторами. Во-первых, потенциальная завихренность по Эртелю – это кинематическая величина [4], по которой нельзя определить интенсивность вихревой структуры течений и даже знак вращения. Во-вторых, часто достаточно рассмотреть потенциальную завихренность по Россби [5] в квазигенострофическом приближении, которая является динамической характеристикой [4] и содержит необходимую информацию о динамике течений.

Для анализа атмосферных прогнозов большую роль играет работа [6], в которой анализируются возможности использования изоэнтропических карт потенциальной завихренности для представления некоторых динамических процессов в атмосфере. Приведены примеры из оперативного анализа погоды и из идеализированных теоретических моделей для иллюстрации этого подхода и его связи с классическими синоптическими концепциями. Обсуждаются структура, причины формирования и устойчивости циклонов и блокирующих антициклонов, физических механизмов распространения волн Россби, бароклинной и баротропной неустойчивости в пространстве и во времени.

В работе [4] анализируется понятие «потенциальная завихренность» и рассматриваются основные соотношения, по которым проводятся ее расчеты, исследуются подходы Россби и Эртеля. В качестве примера приводятся оценки по данным наблюдений потенциальной завихренности для квазипостоянного антициклонического Лофотенского вихря Норвежского моря. Авторами показано, что потенциальная завихренность по Эртелю является кинематической характеристикой, а потенциальная завихренность по Россби в квазигеострофическом приближении – динамической.

Анализ уравнения потенциальной завихренности по Эртелю позволяет оценить вклад нелинейных и диффузионных эффектов в баланс сил, определяющих эволюцию PV. Данная работа является продолжением исследований [7]. Цель ее – получить дискретное уравнение потенциальной завихренности стратифицированной жидкости в квазистатическом приближении как следствие исходной конечно-разностной системы уравнений модели динамики Черного моря и на примере черноморской циркуляции с реалистическими атмосферными условиями для 2011 г. провести анализ полученного уравнения потенциальной завихренности.

Уравнение потенциальной завихренности стратифицированной несжимаемой жидкости в квазистатическом приближении

В приближении Буссинеска и квазистатики в декартовой системе координат движение жидкости в области Ω с границей ∂Ω в форме Громеки – Лэмба описывается следующей системой уравнений:

$\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}(T\overrightarrow{U})={({\kappa }^T{T}_z)}_z-{\kappa }^H{\nabla }^{4}T$, (4)

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

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

Введены обозначения: $\overrightarrow{U}=({\overrightarrow{U}}_H,w)=(u,v,w)$ – компоненты вектора скорости течения по осям (x, y, z), направленным на восток, север и вертикально вниз соответственно; $\overrightarrow{F}=(F^u,F^v)$, $\overrightarrow{g}$= (0, 0, g) – ускорение свободного падения; (T, S, P, ρ) – температура, соленость, давление и плотность морской воды; ρ₀ = 1 г/см³ (здесь и в дальнейшем полагаем давление и плотность нормированными на ρ₀ ); $\overrightarrow{f}$= (0, 0, $f^z$) – параметр Кориолиса, где $f^z$ = 2ω sin ; ω  – угловая скорость вращения Земли; $\phi$ – широта.

В уравнении (1) с учетом квазистатического приближения введены абсолютный вихрь скорости и кинетическая энергия движения:

$\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={({\mu }^{\text{ver}}{u}_z)}_z-{\mu }^{\text{hor}}{\nabla }^{4}u,F^v={({\mu }^{\text{ver}}{v}_z)}_z-{\mu }^{\text{hor}}{\nabla }^{4}v,\overrightarrow{F}=(F^u,F^v),$

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

где ${\mu }^{\text{ver}},{\mu }^{\text{hor}}$ – коэффициенты вертикального и горизонтального обмена количеством движения.

На поверхности при 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)

на дне при z = H(x, y)

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

Использованы следующие обозначения: (${\tau }^x,{\tau }^y)$ – касательное напряжение трения ветра; $Q^T$ – поток тепла; Ev – испарение морской воды; Pr – осадки; S₀ – модельная соленость на поверхности моря; ${\rho }_1$ – плотность морской воды на поверхности моря.

Функции ${\mu }^{ver},{\kappa }^T,{\kappa }^S$ рассчитывались в соответствии с параметризацией Меллора – Ямады [6].

На твердых боковых стенках

для меридиональных участков границы:

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

для зональных участков границы:

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

На участках границы, где вода втекает, используются следующие условия:

для меридиональных участков

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

для зональных участков

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

Для верхнебосфорского течения и для Керченского пролива, когда течение направлено из Черного моря в Азовское,

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

При $t=t^0$ задаются следующие начальные условия:

$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)

Исходя из системы уравнений (1)–(6), выведем уравнение Эртеля. Применяя соответствующую операцию из соотношений (1)–(2), с учетом уравнения неразрывности (3) получаем уравнение для $\overrightarrow{\xi }$

$\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)

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

Следствием уравнений (4), (5) является уравнение для плотности

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

в котором введены обозначения

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

Потенциальная завихренность несжимаемой жидкости в квазистатическом приближении имеет вид

$\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}+[(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y})+f^z]\frac{\partial \rho }{\partial z}.$ (18)

Тогда из уравнений (15), (16) с учетом выражений (3), (17) и (18) следует уравнение Эртеля в квазистатическом приближении для вязкой жидкости

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

где

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

Дискретное уравнение потенциальной завихренности в квазистатическом приближении

В соответствии с введенными в работе [7] разностными операторами выпишем дифференциально-разностные уравнения модели (1)–(6) (дифференциальные по времени):

$\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(u_{i,j,k}^{}{T}_{i,j,k}^{})+{\delta }_y(v_{i,j,k}^{}{T}_{i,j,k}^{})+{\delta }_z(w_{i,j,k}^{}{T}_{i,j,k}^{})=$

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

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

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

${\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}.$

В квазистатическом приближении компоненты вихря скорости (рис. 1) имеют вид

${\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)

Из представления (28) следует, что в вершинах бокса ($i,j,k$) (рис. 1) выполняется важное соотношение

${\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.$

Полагаем, что слагаемые в квадратных скобках в левой части уравнений (21)–(22) записываются в виде, обеспечивающем сохранение энстрофии и энергии [9] для приближения мелкой воды, и соответствуют формулам (32) работы [10].

Уравнения для компонентов абсолютного вихря скорости (аналог уравнения (15)) – для ${\xi }^x$ в точке $i,j+1/\mathrm{2,}k+1/2$, для ${\xi }^y$ в точке $i+1/\mathrm{2,}j, k+1/2$ и для ${\xi }^z$ в точке $i+1/\mathrm{2,}j+1/\mathrm{2,}k$ – с учетом вязкости следуют из соотношений (21)–(24) и имеют вид 

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

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

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

Рис. 1. Распределение переменных в боксе (i, j, k). В вершинах бокса, обозначенных звездочкой, определена PV (ω), на его ребрах – компоненты абсолютного вихря скорости ${\xi }_{}^x,{\xi }_{}^y,{\xi }_{}^z$

Fig. 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

Введем обозначения:

${\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.$

Уравнение плотности в точке (i, j, k) является следствием соотношений (25)–(27) и записывается следующим образом:

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

где

$\begin{array}{l}{D}_V^{\rho }={\alpha }_1^T{\delta }_z[{\kappa }^V({\delta }_{z}T)]+{\alpha }_1^S{\delta }_z[{\kappa }^V({\delta }_{z}S)]+2{\alpha }_2^T[{\delta }_z[{\kappa }^V{\overline{T^{}}}^z({\delta }_{z}T)]+\\ +{\alpha }^{TS}\left[T{\delta }_z[{\kappa }^V({\delta }_{z}S)]+S{\delta }_z[{\kappa }^V({\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].$

Проводя аналогичные [5] преобразования с уравнениями (29), (30), получаем в точке $i+1/\mathrm{2,}j+1/\mathrm{2,}k+1/2$ уравнение потенциальной завихренности в квазистатическом приближении

$\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({\overline{\rho }}^{yz})+$$+{\overline{{\Lambda }^y}}^y{\delta }_y({\overline{\rho }}^{xz})+{\overline{{\Lambda }^z}}^z{\delta }_z({\overline{\rho }}^{xy})+{\overline{{\xi }^x}}^x{\delta }_x({\overline{D^{\rho }}}^{yz})+{\overline{{\xi }^y}}^y{\delta }_y({\overline{D^{\rho }}}^{xz})+{\overline{{\xi }^z}}^z{\delta }_z({\overline{D^{\rho }}}^{xy}),$ (31)

где правая часть уравнения (31) является разностным аналогом выражения (20). С учетом квазистатического приближения и уравнений (29) приняты следующие обозначения:

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

${\Upsilon }_{i+1/\mathrm{2,}j,k+1/2}^y={\delta }_z{([w,{\xi }_{}^y])}_{i+1/\mathrm{2,}j,k+1/2}-{\delta }_z{([v,{\xi }_{}^z])}_{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(\right[w,{\xi }_{}^z\left]\right)}_{i+1/\mathrm{2,}j+1/2,k}+{\delta }_y{\left(\right[v,{\xi }_{}^z\left]\right)}_{i+1/\mathrm{2,}j+1/2,k}-\\ -{\delta }_x{\left(\right[w,{\xi }_{}^x\left]\right)}_{i+1/\mathrm{2,}j+1/2,k}-{\delta }_y{\left(\right[w,{\xi }_{}^y\left]\right)}_{i+1/\mathrm{2,}j+1/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}$

Разностный аналог потенциальной завихренности Эртеля имеет вид

${\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)

где обозначения ${\varpi }_{}^x,{\varpi }_{}^y,{\varpi }_{}^z$ очевидны.

Отличие уравнения (31) от уравнения (45) в [7] состоит не только в том, что учтена вязкость и диффузия, но и в том, что компоненты абсолютного вихря скорости имеют вид (28).

Численный анализ компонентов уравнения потенциальной завихренности по результатам расчета циркуляции с атмосферными условиями для 2011 г.

В численных прогностических экспериментах задавались следующие параметры. Расчеты проводились с равномерным шагом по горизонтальным координатам 1,6 км, по вертикали использовались 27 горизонтов со сгущением в верхнем слое моря. Учет стока причерноморских рек и расходы через Босфорский и Керченский проливы (условия (9)–(13)) соответствовали данным работы [11]. Температура воды в устьях рек (условия (11)–(12)), кроме рек Турции, задавалась из [11]. Предполагалось, что температура рек Турции равна температуре прибрежных вод моря. В верхнебосфорском течении температура и соленость принимались теми же, что и в море, в соответствии с условиями (11). В нижнебосфорском потоке соленость принималась равной 35 ‰, а температура 16 °С.

Для задания атмосферного воздействия в уравнениях (7), (8) использовались данные SKIRON за 2011 г. [12], вертикальное перемешивание описано на основе теории Меллора – Ямады [8]. Начальные условия (14) в этом расчете соответствовали 1 января 2011 г. Расчет был проведен на год модельного времени, его параметры и результаты подробно описаны в [13].

В качестве примера рассмотрим PV (формула (32)) для двух моментов времени – зимнего (рис. 2, а) и летнего (рис. 2, b) периода, в которые структура циркуляции заметно различается.

Рис. 2. Приведенный уровень моря на 1 февраля (а) и 1 августа (b) 2011 г.

Fig. 2. Spesified/Stated sea level on February 1 (a) and August 1 (b), 2011

На 1 января 2011 г. приведенный уровень моря представлял собой обширный циклонический круговорот с двумя синоптическими вихрями – Севастопольским и юго-западным антициклонами. В отличие от зимней циркуляции, в летний сезон (рис. 2, b) циклонический круговорот распадается на два. В западной части бассейна вдоль свала глубин Основное Черноморское течение распространяется в виде узкого струйного течения (рис. 2, а). Зимой вдоль северной периферии северо-западного шельфа образуются мощные меандры, а в юго-восточном углу моря формируется интенсивный антициклон.

Зимней циркуляции соответствует потенциальная завихренность, представленная на рис. 3, а – d. Наблюдаются качественные различия в ее структуре по глубине. В верхнем 30-метровом слое большие значения PV концентрируются в двух областях (рис. 3, а, b). Первая – северо-западная часть моря, ограниченная приблизительно координатами 44°–46° с. ш., 29°–31° в. д., в которой динамика вод во многом определяется стоком рек, в первую очередь Дуная. Поэтому имеют место большие пространственные градиенты в поле плотности, что и обусловливает значения  в этой зоне моря. Вторая область – свал глубин, где выделяется Анатолийское побережье и южная периферия северо-западного шельфа (рис. 3, b).

Рис. 3. PV на горизонтах 3,75 (а), 35 (b), 106,25 (с) и 350 м (d) на 1 февраля 2011 г.

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

В области Севастопольского и юго-западного вихря экстремальных значений $\varpi$ не наблюдается. Объясняется это тем, что $\varpi$ является скалярной величиной, равной произведению градиента плотности и абсолютного вихря скорости, и поэтому большая его величина не означает обязательного увеличения завихренности и, более того, его знак не определяет знак вращения вихря [2]. Ниже верхнего 50-метрового слоя наибольшие значения  наблюдаются в прибрежной полосе (рис. 3, с, d). Локальные максимумы сосредоточены в относительно малых зонах (~ 10 км), которые отчетливо проявляются на глубине 100 м (рис. 3, с). В центральной части моря структура PV довольно однородна по пространству.

Показательной иллюстрацией для анализа PV является его расчет для летнего периода, когда циркуляция носит менее регулярный характер и ее вихревая структура выражена более ярко (см. рис. 2, b). Две особенности проявляются в структуре потенциальной завихренности в августе 2011 г. Полоса малых значений  в верхнем слое (рис. 4, а, b) соответствует области наибольшего перемешивания по вертикали поля плотности, то есть значение $\delta {}_z\overline{\rho }{}_{}^{xy}$ очень мало – на два – три порядка меньше, чем вертикальные градиенты плотности окружающей воды. Вторая особенность заключается в однородной структуре и малых значениях PV в областях, примерно соответствующих ядрам юго-восточного антициклона, юго-западного и восточного круговоротов. Эти особенности также определяются структурой $\delta {}_z\overline{\rho }{}_{}^{xy}$, которая незначительно меняется по пространству. Такой вид потенциальной завихренности в центральных частях круговоротов согласуется с выводами работы [4], в которой PV, восстановленный по данным наблюдений, в области Лофотенского круговорота имеет аналогичную структуру. На нижних горизонтах (рис. 4, с, d) вдоль границы области вследствие перепадов в рельефе дна имеет место узкая полоса неоднородных значений $\delta {}_z\overline{\rho }{}_{}^{xy}$, в центральной части наблюдается малая пространственная изменчивость этой величины.

Рис. 4. PV на горизонтах 3,75 (а), 35 (b), 106,25 (c) и 350 м (d) на 1 августа 2011 г.

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

Основной вклад в структуру вихря, как правило, вносит составляющая ${\varpi }_{}^z$ [4], величина которой определяется квазигеострофическим характером движения и вертикальной стратификацией морской воды. В областях, где имеет место вток пресных (устья рек) или соленых (проливы) вод моря, преобладающее значение могут иметь горизонтальные составляющие потенциальной завихренности. В качестве примера на рис. 5 приведены значения ${\varpi }_{}^x,{\varpi }_{}^y,{\varpi }_{}^z$ на горизонтах 3,75 и 106,25 м.

Сопоставление рис. 5, а, c, e и 3, а показывает, что горизонтальная составляющая ${\varpi }^x$ (рис. 5, а) вносит основной вклад в зоне втока пресных вод Дуная в северо-западной части моря, в остальной области – структуру вихря определяет ${\varpi }_{}^z$ (рис. 5, e).

Рис. 5. Составляющие PV: ${\mathit{\varpi }}_{}^{\mathbf{x}}$ на горизонтах 3,75 (а), 106,25 м (b), ${\mathit{\varpi }}_{}^{\mathbf{y}}$ на горизонтах 3,75 (с), 106,25 м (d), ${\mathit{\varpi }}_{}^{\mathbf{z}}$ на горизонтах 3 (e), 106,25 м (f) на 1 февраля 2011 г.

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

На горизонте 106,25 м (рис. 5, b, d, f) вертикальная составляющая на порядок больше ${\varpi }^x$ (рис. 5, b), ${\varpi }^y$ (рис. 5, d), поэтому она (рис. 5, f) довольно точно определяет вид потенциальной завихренности на горизонте 106,25 м (см. рис. 3, c).

Прямыми вычислениями устанавливается, что в верхнем слое моря вид составляющей ${\varpi }_{}^z$ качественно соответствует ${\delta }_z{\overline{\rho }}^{xy}$, но она по абсолютной величине на несколько порядков меньше, чем ${\overline{{\xi }^z}}^z$. В свою очередь структура абсолютного вихря довольно однородна и положительна, поэтому при умножении ${\delta }_z{\overline{\rho }}^{xy}$ на ${\overline{{\xi }^z}}^z$ структуру потенциальной завихренности характеризует ${\delta }_z{\overline{\rho }}^{xy}$, а ее количественное значение зависит от ${\overline{{\xi }^z}}^z$. Величина ${\xi }^z$ определяется двумя слагаемыми – относительной и планетарной завихренностью. Если оценивать вклад $f^z$ в абсолютный вихрь, то $f^z$ сравнимо по величине с относительным вихрем и увеличивает значения ${\xi }^r$. По величине ${\xi }^z$ в среднем больше на два порядка, чем ${\delta }_z{\overline{\rho }}^{xy}$. В этот период года в результате зимней конвекции сомножитель, обусловленный градиентом плотности по вертикали, в верхнем слое моря мал, за исключением области стока рек, где его величина может быть значима.

На горизонте 106,25 м (рис. 5, f) оба сомножителя положительные ${\overline{{\xi }^z}}^z$ и в среднем меньше на несколько порядков. Изменчивость PV наблюдается во вдольбереговой полосе, в центральной части моря PV однородна. Отметим, что, во-первых, на нижних горизонтах (приблизительно ниже глубины 50 м) относительный вихрь по абсолютной величине на порядок меньше $f^z$. Во-вторых, так как интеграл по горизонтальной поверхности от ${\xi }^r$, отличие от нуля которого определяется стоком рек и обменом воды через проливы, мал, то в структуре ${\xi }^r$ имеют место зоны циклонического и антициклонического вращения вод. В то же время планетарная завихренность положительная и больше ${\xi }^r$, и поэтому именно она определяет количественные значения PV с поправками, вносимыми относительным вихрем, на качественную структуру потенциальной завихренности на глубинных горизонтах.

Рассмотрим вклад нелинейных сил в эволюцию $\varpi$. Введем обозначения

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

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

Основной вклад во временную эволюцию PV нелинейные силы в верхнем слое вносят в прибрежной области моря (рис. 6, а, b). Их вклад для различных областей неодинаков: больше в северо-западной части (рис. 6, а) и вдоль Анатолийского побережья (рис. 6, b). Оценки показывают, что их количественные различия по абсолютной величине между центральной частью моря и его периферией составляют несколько порядков. На нижних горизонтах (рис. 6, с, d) наибольшие значения нелинейных слагаемых в уравнении Эртеля концентрируются в виде узкой вдольбереговой полосы с более ярко выраженным характером у южного берега моря.

Рассмотрим вклад слагаемых отдельных слагаемых $C^x,C^y,C^z$ в $C^S$.

Рис. 6. ${\mathit{C}}^{\mathbf{S}}$ на 1 февраля 2011 г. на горизонтах 3 (а), 106.25 м (b) и на 1 августа 2011 г. 3,75 (c), 106.25 м (d)

Fig. 6. ${\mathit{C}}^{\mathbf{S}}$ at horizons 3 (a) and 106.25 m (b) on February 1, 2011, and at horizons 3.75 (c) and 106.25 m (d) on August 1, 2011

В верхнем слое (рис. 7, а, c, e) зоны больших по абсолютной величине значений $C^x,C^y,C^z$ имеют сходную структуру. В юго-восточном углу бассейна и в северо-восточной части моря, ограниченной координатами 42°– 44° с. ш., 37°–39° в. д., наблюдаются области значений $C^x,C^y,C^z$, близких к нулю. Вычисленные средние и максимальные значения $C^x,C^y,C^z$ (рис. 7, а, c, e) в сопоставлении с  $C^S$ (рис. 6, а) свидетельствуют о том, что экстремальные значения различаются на порядок, средние – на два порядка. Это означает, что $C^x,C^y,C^z$ взаимно компенсируются и в результате получается структура, представленная на рис. 6, а. Прямыми вычислениями устанавливается, что основной вклад в нелинейные слагаемые $C^x,C^y,C^z$ в верхнем слое дают $C_1^x,C_1^y,C_1^z$ соответственно, то есть ${\delta }_x({\Upsilon }^x{\overline{\rho }}^{yz}),{\delta }_y({\Upsilon }^y{\overline{\rho }}^{xz}),{\delta }_z({\Upsilon }^z{\overline{\rho }}^{xy}).$ Оценка порядка величин показывает, что $C^S$ по величине в среднем на два порядка меньше каждого из слагаемых $C^x,C^y,C^z$.

Рис. 7. Составляющие ${\mathit{C}}^{\mathbf{S}}$: ${\mathit{С}}_{}^{\mathbf{x}}$ на горизонтах 3,75 м (а), 106.25 м (b), ${\mathit{C}}^{\mathbf{y}}$ на горизонтах 3.75 м (c), 106.25 м (d), ${\mathit{С}}_{}^{\mathbf{z}}$ на горизонтах 3,75 м (e), 106,25 м (f) на 1 февраля 2011 г.

Fig. 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 February 1, 2011

Аналогичная ситуация имеет место и для рассчитанных полей на 1 августа 2011 г.

Заключение

Для системы дискретных уравнений динамики моря в приближении гидростатики и с учетом вязкости, диффузии, втока рек, водообмена через проливы и атмосферного воздействия как ее следствие получено конечно-разностное уравнение потенциальной завихренности стратифицированной несжимаемой жидкости. Так же, как и в более общем случае, оно имеет дивергентный вид и отличается от своего дифференциального аналога. Поскольку используется нелинейное состояние для расчета плотности, полученное дискретное уравнение для PV не является точным следствием конечно-разностных уравнений модели. Для оценки влияния нелинейного характера уравнения состояния на полученные результаты необходимы дополнительные исследования.

Анализ величины самой потенциальной завихренности подтвердил полученные ранее результаты, что ее вертикальная компонента является основной. Горизонтальные составляющие вносят заметный вклад в областях стока рек, водообмена через проливы и в зонах резких градиентов поля плотности. Качественный вид PV в верхнем слое моря имеет аналогичные особенности, зафиксированные по данным наблюдений. Однородная структура для центральной части вихревых образований и интенсивный характер в области больших поперечных градиентов в поле плотности во фронтальных зонах определяют структуру потенциальной завихренности. В глубинных слоях моря ее наибольшие значения концентрируются в виде узкой вдольбереговой полосы, в остальной части моря значения PV малы.

Расчет составляющих ${\varpi }^z$ в зимний период на верхних горизонтах показал, что в верхних слоях моря, за исключением прибрежных зон стока рек, определяется ${\xi }^z$, который является суммой относительного вихря и $f^z$ (величины примерно одного порядка). В нижних слоях моря количественные значения PV в большей степени определяет планетарная завихренность, а ее качественные особенности обусловлены структурой относительного вихря.

Из анализа нелинейных слагаемых в уравнении PV следует, что в верхнем слое моря основной вклад в адвекцию потенциальной завихренности $C^S$ вносится в северо-западной части и вдоль Анатолийского побережья. На нижних горизонтах наибольшие значения $C^S$ наблюдаются вдоль береговой полосы с более ярко выраженным характером у южного берега моря, что соответствует структуре PV.

Расчет слагаемых $C^x,C^y,C^z$ для зимнего и летнего периодов позволил установить два факта. Во-первых, величина каждого $C^x,C^y,C^z$ определяется ${\delta }_x({\Upsilon }^x{\overline{\rho }}^{yz}),{\delta }_y({\Upsilon }^y{\overline{\rho }}^{xz}),{\delta }_z({\Upsilon }^z{\overline{\rho }}^{xy}),$ то есть дивергенцией от произведения нелинейных слагаемых в уравнениях движения и плотности. Во-вторых, $C^S$ на полтора – два порядка по величине меньше, чем каждая из составляющих $C^x,C^y,C^z$, то есть локально полная сумма вертикальной и горизонтальной адвекции потенциальной завихренности на два порядка меньше, чем каждая по отдельности. Возможное объяснение такого результата заключается в следующем. Предположим, что конечно-разностные аналоги нелинейных слагаемых в уравнении PV близки к дифференциальному виду $\text{div}(\overrightarrow{U}\omega ).$ Тогда, представляя $\omega ={\omega }^S+{\varpi }^{*},$ где ${\omega }^S$ – величина, осредненная по пространству, в каждой точке области имеем $\text{div}(\overrightarrow{U}{\omega }^S)={\omega }^S\text{div}(\overrightarrow{U})=0$ или близко к нулю. Поскольку изменчивость во времени потенциальной завихренности в преобладающей степени зависит от $\text{div}(\overrightarrow{U}{\omega }^{*})$ и происходит взаимная компенсация нелинейных составляющих по x, y, z при расчете адвекции PV.

Насколько полученные результаты имеют общий характер – вопрос дальнейших исследований.

¹ См. Часть I: Демышев С. Г. Конечно-разностная аппроксимация уравнения потенциальной завихренности для стратифицированной несжимаемой жидкости и пример его использования при расчете циркуляции Черного моря. Часть I. Дифференциально-разностное уравнение потенциальной завихренности идеальной жидкости // Морской гидрофизический журнал. 2024. Т. 40, № 2. С. 165–179. EDN BCKKBN.

About the authors

S. G. Demyshev

Author for correspondence.
Email: demyshev@gmail.com
ORCID iD: 0000-0002-5405-2282
SPIN-code: 1848-2350
Scopus Author ID: 6603919865
ResearcherId: C-1729-2016

зав. отделом теории волн, главный научный сотрудник, доктор физико-математических наук

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 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 (1MB)

Indexing metadata

3. Fig. 2. Spesified/Stated sea level on February 1 (a) and August 1 (b), 2011

Download (349KB)

Indexing metadata

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

Download (699KB)

Indexing metadata

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

Download (638KB)

Indexing metadata

6. Fig. 5. PV components: at horizons 3.75 (a) and 106.25 m (b); at horizons 3.75 (c) and 106.25 m (d); at horizons 3 (e) and 106.25 m (f) on February 1, 2011

Download (1MB)

Indexing metadata

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

Download (1MB)

Indexing metadata

8. Fig. 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 February 1, 2011

Download (5MB)

Indexing metadata