Modeling of the Black Sea circulation using equations of heat and salt advection–diffusion having discrete nonlinear invariants

全文:

Введение

Для воспроизведения и исследования трехмерных гидрофизических полей морей и океанов, непрерывных по времени и пространству, преимущественно используются численные модели динамики морской среды. Достоверность результатов расчетов зависит, в том числе, и от свойств разностных схем, на основе которых строятся модели диагноза и прогноза морских течений. Одним из методов уменьшения ошибки моделирования, особенно при длительных интегрированиях, является построение и использование разностных схем, удовлетворяющих законам сохранения, что по аналогии с непрерывным случаем должно быть связано с инвариантностью системы дискретных уравнений (например, [Самарский и др., 1997; Cheviakov et al., 2020]). Дискретными инвариантами гиперболической системы уравнений динамики океана в адиабатическом приближении и при отсутствии трения, диффузии и внешних источников являются интегралы по объему от плотности, полной (кинетическая плюс потенциальная) энергии, температуры, солености в степени 2 и более. При нелинейной зависимости плотности от температуры T и солености S для сохранения дискретного интеграла от плотности целесообразно использовать такие аппроксимации нелинейных слагаемых, чтобы наряду с T и S сохранялись T^K и S^L, где K и L – целые положительные числа. Наличие дискретных инвариантов T^K и S^L при нелинейном виде уравнения состояния позволяет при указанных условиях сохранять плотность и полную энергию, что повышает устойчивость разностного решения и его близость к точному.

Убедительным практическим примером важности инвариантности дискретных уравнений для атмосферных и океанических моделей служит известная аппроксимация Аракавы–Лэмба [Arakawa et al., 1981], которую удалось обобщить на схему дискретную по пространству и по времени [Sorgentone et al., 2015]. Она обладает двумя квадратичными инвариантами, что обеспечивает сохранение среднего волнового числа, взвешенного по энергии. Этим самым предотвращается передача энергии по спектру в сторону одних масштабов, например, в малые масштабы. Поскольку в численных моделях наиболее короткие волны – двухшаговые, групповая скорость которых нулевая, то используемые схемы с двумя квадратичными законами сохранения препятствуют накоплению энергии на этих масштабах.

Продолжением этих исследований являются работы [Scott et al., 1994; Palha et al., 2017; Головизнин и др., 1998; Goloviznin et al., 2022; Демышев, 2023]. Для линейных и нелинейных одномерных волновых уравнений на пяти- и девятиточечных шаблонах в [Cheviakov et al., 2020] получены конечно-разностные уравнения, сохраняющие симметрию и законы сохранения. В случае линейного волнового уравнения выписана явная пятиточечная схема, сохраняющая основные геометрические симметрии и шесть соответствующих законов сохранения. Для нелинейных волновых уравнений построена девятиточечная неявная схема, сохраняющая четырехточечные симметрии и три локальных закона сохранения. В работе [Капцов, 2019] для одномерных уравнений мелкой воды в лагранжевых координатах разработана инвариантная конечно-разностная схема, которая обладает локальными законами сохранения энергии, массы, центра масс и импульса.

Ряд работ посвящен выбору дискретных моделей для уравнений Навье–Стокса, обладающих законами сохранения или инвариантности относительно преобразования координат. В [Scott et al., 1994] рассматриваются двумерные уравнения Навье–Стокса для несжимаемой жидкости, выводятся необходимые и достаточные условия, при которых дискретная аппроксимация удовлетворяет закону сохранения в пространстве сеточных интегрируемых функций. Представлены численные результаты, демонстрирующие способность схемы точно разрешать формирующийся пограничный слой. В работе [Palha et al., 2017] предложена специальная дискретизация двумерных уравнений Навье–Стокса для несжимаемой жидкости, которая в пределе исчезающей диссипации точно сохраняет массу, кинетическую энергию, энстрофию и полную завихренность. Приведены доказательства сохранения точных дискретных свойств и демонстрационные численные тестовые примеры на нерегулярных треугольных сетках.

Уравнения адвекции–диффузии тепла и соли являются неотъемлемой частью моделей динамики морей и океанов, что требует их аккуратной аппроксимации. В оригинальных работах [Головизнин и др., 1998; Goloviznin et al., 2022] для одномерного уравнения переноса разработана и апробирована схема КАБАРЕ, которая представляет собой трехслойную явную разностную схему второго порядка аппроксимации и является точной при числах Куранта 0.5 и 1. Представлены результаты по применению новой малодиссипативной многослойной гидростатической модели, которая описывает динамику жидкости с переменной плотностью и свободной поверхностью. Для решения системы гиперболических уравнений в каждом слое используется явная схема КАБАРЕ. Результаты численных расчетов удовлетворительно согласуются с экспериментальными данными. В работе [Демышев, 2023] для уравнения адвекции непрерывного по времени и разностного по пространственным переменным получены схемы, которые обеспечивают сохранение одновременно температуры в первой и в K-й (K > 1), солености в первой и L-й (L > 1) степени. Найденные аппроксимации температуры и солености при полиномиальной зависимости плотности от температуры и солености также приводят к дивергентному виду уравнения адвекции плотности. Такой его вид обеспечивает выполнение закона сохранения полной энергии в дискретной постановке.

Цель настоящей работы – оценить точность результатов моделирования при использовании схем, которые обеспечивают сохранение температуры и солености в первой и в третьей, пятой степенях в сравнении с традиционной аппроксимацией нелинейных слагаемых в уравнениях адвекции–диффузии тепла и соли. На основе сопоставления данных контактных и спутниковых наблюдений с результатами прогностических расчетов циркуляции Черного моря в 2016 г. рассчитаны и проанализированы средние и среднеквадратические ошибки полей температуры и солености.

Постановка задачи

Численные эксперименты проведены с использованием нелинейной вихреразрешающей модели Морского гидрофизического института [Демышев и др., 1992; Demyshev et al., 2022]. Модель построена на основе полной системы уравнений термогидродинамики океана в приближении Буссинеска и гидростатики в декартовой системе координат. Уравнение состояния представлено в виде полиномиальной зависимости плотности от температуры и солености во второй степени. Изменением давления мы пренебрегаем, что допустимо в Черном море вследствие малости его глубины по сравнению с океаном. Вертикальное турбулентное перемешивание описано моделью замыкания Меллора–Ямады 2.5 [Mellor et al., 1982]. В качестве граничных условий на свободной поверхности задаются напряжение ветра, потоки тепла, осадки и испарение. На твердых участках границы задано отсутствие нормальных потоков импульса, тепла и соли, на дне – условие прилипания. В модели учитываются климатический сток рек и водообмен через проливы [Гидрометеорология, 1991]. В начальный момент времени задаются высота уровня моря, температура, соленость и горизонтальная скорость. Один раз в сутки в модели производится усвоение спутниковой температуры поверхности моря [https://data.marine.copernicus.eu/product/SST_BS_SST_L3S_NRT_OBSERVATIONS_010_013]. Батиметрия бассейна построена на основе цифрового массива глубин EMODnet [https://emodnet.ec.europa.eu/geonetwork/srv/eng/catalog.search#/metadata/19f800a9-f0fd-4055-b4cd-90ed156dc7fc]. Конечно-разностная аппроксимация уравнений модели, граничных и начальных условий по пространственным переменным реализована на сетке C и подробно описана в работе [Демышев, 2012]. Для аппроксимации производных по времени используется схема “чехарда“ с периодическим подключением схемы Матсуно.

В работе проведено три численных эксперимента при одинаковых начальных и граничных условиях. Для задания потоков импульса, тепла и влаги на поверхности моря использован атмосферный реанализ ERA5 [https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5]. Начальные поля взяты из физического реанализа Черного моря CMEMS [https://data.marine.copernicus.eu/product/BLKSEA_MULTIYEAR_PHY_007_004/description]. Расчеты выполнены с горизонтальным разрешением 1.64 км, что соответствует примерно (1/48)˚ долготы и (1/66)˚ широты, левый нижний узел соответствует точке с координатами 27.34˚ в.д., 40.81˚ с.ш. По вертикали задано 27 горизонтов с глубинами 2.5; 5; 10; 15; 20; 25; 30; 40; 50; 62.5; 75; 87.5; 100; 112.5; 125; 150; 200; 300; 400; 500; 700; 900; 1100; 1300; 1500; 1700; 2100 м. Шаг по времени составляет 96 с. Коэффициенты при бигармонических операторах Лапласа, описывающих горизонтальную турбулентную вязкость и диффузию, равны 10¹⁶ см⁴/с. Коэффициенты вертикальной вязкости и диффузии, изменяющиеся по времени и по глубине, рассчитываются на каждом шаге по времени по [Mellor et al., 1982]. Эмпирические константы в параметризациях вертикального обмена, коротковолновой солнечной радиации на поверхности моря и другие постоянные коэффициенты приведены в Appendix работы [Demyshev et al., 2022]. Численные эксперименты различались аппроксимациями адвективных слагаемых в уравнениях переноса тепла и соли в соответствии с результатами работы [Демышев, 2023] (см. Приложение). В первом расчете T2S2 использовалась традиционная аппроксимация нелинейных слагаемых (П.1), при которой сохранялись T, S и T², S². Во втором эксперименте T5S3 использовалась схема, сохраняющая T, S и T⁵ , S³ (П.2), что соответствует наибольшим степеням температуры и солености в рекомендованном ЮНЕСКО уравнении состояния [IOC, 2010]. В отличие от океана, где плотность в основном определяется температурой и уравнение состояния [IOC, 2010] представлено в виде полинома с различными степенями по температуре, в Черном море плотность в преобладающей степени зависит от солености [Булгаков и др., 1984]. Поэтому третий эксперимент T3S5 проводился на основе аппроксимации адвективных слагаемых по формуле (П.3), которая обеспечивает сохранение T, S и . Все расчеты выполнены на один год с 1 января по 31 декабря 2016 г. Выходными данными модели являются среднесуточные поля уровня моря, температуры, солености и скорости течений.

Сопоставление с данными наблюдений

Для валидации результатов численных экспериментов использованы данные контактных измерений температуры и солености, полученные 14-ю буями-профилемерами ARGO [https://www.coriolis.eu.org/Data-Products/Data-selection] и в 87, 89, 91-м рейсах НИС “Профессор Водяницкий“ [Артамонов и др., 2018] в 2016 г. Натурные данные имеются за весь 2016 г. и охватывают практически всю глубоководную зону и часть зоны шельфа западнее Крыма. Данные были разделены по вертикальным слоям, районам и сезонам. На рис. 1 представлены карты, показывающие расположение станций измерений для теплого (июнь–август) и холодного (октябрь–декабрь) сезонов.

Рис. 1. Карты расположения станций измерений температуры и солености зимой (а) и летом (б) 2016 г. по данным НИС “Профессор Водяницкий” и буев-профилемеров ARGO.

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

$\begin{array}{c}mean=\frac{1}{N}{\sum }_{n=1}^N\left(q_n^{insitu}-q_n^{\mathrm{mod}el}\right),\\ RMSE=\sqrt{\frac{{\displaystyle {\sum }_{n=1}^N{(q_n^{insitu}-q_n^{\mathrm{mod}el})}^2}}{N},}\end{array}$

где q – исследуемый параметр по натурным (in situ) или модельным (model) данным, N – количество измерений. Полученные оценки показывают, что ниже верхнего 300-метрового слоя при изменении степени K и L в уравнениях (П.1)–(П.3) отклонения модельных полей температуры и солености от измеренных меняются незначительно. Поэтому ниже мы приводим результаты валидации до указанного горизонта. Если в выбранный сезон в определенном районе моря (западная, центральная, восточная части и т.п.) действовало несколько буев, то в табл. 1–4 представлены средние RMSE по буям, указанным в заголовках столбцов.

Таблица 1. RMSE температуры и профили температуры в ноябре и декабре 2016 г.

Глубина, м	R/V 91 (16 ноя–3 дек)			A 3901852, 3901854, 3901855 (22 окт–29 дек)			А 6900807, 6901831 (1 окт–30 дек)			A 6901832, 6901833 (1 окт–29 дек)
	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3
0–5	0.767	0.871	0.718	0.546	0.469	0.351	0.695	0.714	0.631
5–30	0.792	0.911	0.737	0.715	0.578	0.456	1.022	0.965	1.027	1.053	1.302	1.252
30–100	1.601	0.966	1.018	1.384	0.425	0.603	0.936	0.422	0.581	1.209	0.736	0.880
100–300	0.305	0.342	0.349	0.263	0.277	0.295	0.125	0.259	0.261	0.341	0.338	0.349

Из анализа RMSE и профилей температуры зимой (табл. 1) следует, что в глубоководной части моря (рис. 1а, буи ARGO 3901852, 3901854, 3901855, 6900807, 6901831) модельная толщина верхнего перемешанного слоя больше, чем реальная. В районе свала глубин в зоне действия Основного черноморского течения ОЧТ (рис. 1а, рейс 91, буи 6901832, 6901833) структура модельных профилей качественно похожа на измеренные, однако для всех расчетов слой скачка температуры залегает в среднем на 20 м ниже, чем по данным наблюдений. При этом разница между результатами трех расчетов наименьшая. В экспериментах T3S5 и T5S3 для всех районов глубина залегания слоя скачка воспроизводится точнее по сравнению с данными T2S2, что демонстрирует величина RMSE в слое 30–100 м. Для всех экспериментов в окрестностях слоя 100–150 м обнаруживается ядро вод пониженной температуры (менее 8.5˚С), которое в натурных профилях практически не проявляется.

В поле солености зимой (табл. 2) наибольшие RMSE для всех экспериментов наблюдается в слое 30–100 м в районе свала глубин (рис. 1а, рейсы 91, буи 6901832, 6901833) и слое 100–300 м в глубоководной части моря (рис. 1а, буи ARGO 3901852, 3901854, 3901855, 6900807, 6901831), что связано с отличиями в глубине залегания верхней границы основного галоклина: как видно из поведения профилей, приведенных в табл. 2, для всех станций измерений характерна меньшая толщина верхнего перемешанного слоя, чем по данным моделирования. При этом результаты экспериментов T3S5 и T5S3 показывают более близкие к реальности величины солености, что отражается в величине RMSE солености для верхнего 100-метрового слоя. Наименьшие отклонения модельной солености для всех экспериментов наблюдаются для юго-западной части моря (рис. 1а, буи ARGO 3901852, 3901854, 3901855) – района с интенсивной струйной динамикой и слабой вихревой изменчивостью в зимний период. Для остальных районов в слое 30–100 м RMSE солености, полученной в экспериментах T3S5 и T5S3, меньше примерно на 20%, чем по данным T2S2, а в слое 100–300 м – наоборот, на 15% больше (табл. 2).

Таблица 2. RMSE солености и профили солености в ноябре и декабре 2016 г.

Глубина, м	R/V 91 (16 ноя–3 дек)			A 3901852, 3901854, 3901855 (22 окт–29 дек)			А 6900807, 6901831 (1 окт–30 дек)			A 6901832, 6901833 (1 окт–29 дек)
	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3
0–5	0.179	0.142	0.137	0.332	0.181	0.174	0.286	0.165	0.155
5–30	0.176	0.144	0.138	0.298	0.153	0.183	0.349	0.235	0.203	0.180	0.112	0.130
30–100	0.590	0.531	0.545	0.748	0.481	0.541	0.961	0.911	0.915	0.678	0.583	0.620
100–300	0.653	0.746	0.777	0.558	0.471	0.503	0.403	0.531	0.528	0.619	0.692	0.748

В летний период в Черном море такие факторы, как прогрев вод моря, ослабление ОЧТ и интенсификация вихревой изменчивости, приводят к значительным вертикальным градиентам температуры и поэтому RMSE температуры увеличивается (табл. 3). Наибольшие RMSE до 3˚С в слое 5–30 м наблюдаются в центральной части моря (рис. 1б, рейс 87, буи № 6901900, 7900594) и связаны с неточностью воспроизведения глубины залегания верхней границы слоя скачка температуры. Как видно из структуры профилей температуры летом, толщина верхнего квазиоднородного слоя и верхняя граница слоя скачка температуры по данным измерений расположены глубже примерно на 5 м, чем по данным численных экспериментов. При этом величины RMSE в слое 5–30 м мало отличаются между собой. Также относительно большие значения RMSE температуры для всех экспериментов получены в слое 30–100 м, где залегают воды холодного промежуточного слоя (ХПС). Для всех модельных профилей в слое 50–100 м наблюдается ядро ХПС, где значения температуры менее 8.5˚С, и которое практически отсутствует в данных наблюдений. Анализ RMSE и модельных профилей температуры показывает, что в летний сезон изменение схемы аппроксимации уравнений адвекции тепла и соли мало влияет на точность воспроизведения полей температуры в верхнем 100-метровом слое.

Таблица 3. RMSE температуры и профили температуры в июне–августе 2016 г.

Глубина, м	R/V 87 (30 июн–18 июл)			A 6901900, 7900594 (1 июн–25 авг)			А 6900805, 6901866 (2 июн–25 июл)			A 6901831, 6901832, 6901833 (2 июн–29 авг)
	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3
0–5	0.662	1.158	1.379	1.023	1.187	1.419	0.708	0.967	0.912
5–30	3.026	3.005	3.089	2.970	2.082	2.605	1.143	1.005	1.304	2.862	2.117	2.485
30–100	1.112	1.105	1.079	0.902	0.972	1.028	1.323	1.206	1.177	0.999	0.824	0.883
100–300	0.253	0.254	0.255	0.227	0.249	0.253	0.269	0.295	0.292	0.280	0.279	0.283

На всех профилях солености летом (табл. 4) модельные значения ниже наблюдаемых, т.к. модельная глубина залегания верхней границы основного галоклина ниже, чем по данным измерений, в среднем на 50 м, что указывает на возможную систематическую ошибку в данных наблюдений. Наибольшие значения RMSE выявлены в слое 30–100 м, где в вертикальной структуре поля солености присутствуют наибольшие градиенты. При этом в прибрежной и центральной частях моря для экспериментов T3S5 и T5S3 RMSE солености в верхнем 100-метровом слое ниже по сравнению с данными T2S2. Наименьшая разница между натурной и модельной толщиной основного галоклина обнаружена в центральной глубоководной части моря (рис. 1б, буи 6901900, 790059, 46900805, 6901866). Для зоны свала глубин южнее Крыма и в районе побережья Северного Кавказа наблюдаются наибольшие отклонения между модельными и натурными данными. Для верхнего квазиоднородного слоя (0–30 м) результаты экспериментов T3S5 и T5S3 демонстрируют лучшее качественное и количественное соответствие с данными измерений: RMSE солености ниже для всех районов, а для станций по данным рейса 87 и буев 6901900, 7900594 модельное изменение солености с глубиной качественно соответствует вертикальной структуре реальных профилей.

Таблица 4. RMSE солености и профили солености в июне–августе 2016 г.

Глубина, м	R/V 87 (30 июн–18 июл)			A 6901900, 7900594 (1 июн–25 авг)			А 6900805, 6901866 (2 июн–25 июл)			A 6901831, 6901832, 6901833 (2 июн–29 авг)
	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3	T2S2	T3S5	T5S3
0–5	0.232	0.188	0.217	0.369	0.247	0.294	0.137	0.146	0.146
5–30	0.184	0.129	0.143	0.216	0.126	0.151	0.110	0.087	0.102	0.231	0.216	0.212
30–100	0.556	0.419	0.422	0.859	0.712	0.777	0.376	0.307	0.302	0.791	0.832	0.830
100–300	0.483	0.557	0.574	0.366	0.403	0.486	0.470	0.462	0.452	0.453	0.716	0.684

Таким образом, валидация результатов численных экспериментов показала, что для всех расчетов наблюдается качественное соответствие температуры и солености данным наблюдений. Количественная разница связана с неточностью воспроизведения толщины верхнего перемешанного слоя зимой и верхнего квазиоднородного слоя летом. Наибольшие значения отклонений от данных наблюдений обнаружены в слоях 5–30 м для температуры и 30–100 м для солености. Изменение схемы аппроксимации уравнений адвекции-диффузии тепла и соли улучшает точность воспроизведения толщины верхнего перемешанного слоя в зимний период и глубины залегания верхней границы слоя скачка температуры летом в центральной части моря.

Анализ результатов моделирования

Рассмотрим различия между результатами трех экспериментов. Приведенный уровень моря является важной интегральной характеристикой для анализа эволюции поля течений в верхнем слое, спутниковые данные о котором позволяют оценить соответствие модельных и реальных пространственных особенностей поля течений. Было проведено последовательное по времени сопоставление модельных полей уровня моря и карт альтиметрического уровня, построенных на основе алгоритма [Kubryakov et al., 2012], которые получены в отделе дистанционных методов исследований МГИ РАН. Сравнение показало качественное соответствие между данными альтиметрии и основными динамическими особенностями Черного моря (такими как ОЧТ, Севастопольский и Батумский антициклоны) в центральной и западной глубоководных частях моря. Однако в прибрежной зоне, вследствие недостаточного пространственного разрешения данных альтиметрии (1/8˚), затруднительно оценить мезомасштабную динамику вод. Поэтому для анализа также привлекались спутниковые снимки температуры поверхности моря и концентрации хлорофилла “а” [http://dvs.net.ru/mp/data/main_ru.shtml].

Анализ карт уровня показал, что отличия в результатах проявляются с конца апреля, после прекращения действия зимних ветров, и связаны они с динамикой мезомасштабных антициклонов. Например, на снимке концентрации хлорофилла за 28.04.2016 (рис. 2) наблюдается антициклонический круговорот в окрестностях 31˚ в.д., который вовлекает шельфовые воды с высокой концентрацией хлорофилла и переносит их в глубоководную зону. Судя по контрастности и форме “языка” вод, орбитальная скорость вихря достаточно высокая. Как видно из рис. 2, в экспериментах T3S5 и T5S3 расположенный в этом районе Севастопольский антициклон более интенсивен: площадь вихря и высота уровня больше по сравнению с экспериментом T2S2.

Рис. 2. Снимок концентрации хлорофилла и модельные поля уровня моря 28.04.2016.

Характерные для летнего сезона ослабление ОЧТ и усиление мезомасштабной изменчивости проявляются во всех экспериментах. Однако, в эксперименте T2S2 ширина струи больше и циклонический круговорот охватывает меньшую площадь, чем в экспериментах T3S5 и T5S3. Из рис. 3 видно более плотное расположение изолиний уровня над свалом глубин по результатам T3S5 и T5S3, что свидетельствует о более интенсивном течении. Качественное сопоставление карт уровня со спутниковыми данными показывает, что форма и расположение восточного циклонического круговорота в северо-восточной части моря по данным T3S5 и T5S3 соответствует зоне прогретых поверхностных вод между 43˚ и 45˚ с.ш. (красная область). Также важным отличием между расчетами в летний сезон является динамика вод около Анатолийского побережья. Типичным примером цепочки прибрежных мезомасштабных вихрей в юго-восточной части моря, существование которых подтверждается многочисленными наблюдениями [Akpinar et al., 2022], являются вихревые структуры на снимке температуры поверхности на рис. 3 между 34˚ и 39˚ в.д. Сравнение пространственного положения изолиний уровня и снимка на рис. 3 подтверждает, что эксперимент T5S3 дает результаты наиболее близкие к реальной картине. Также на рис. 3 в поле температуры видны два прибрежных вихря около Новороссийска, черты которых можно обнаружить в экспериментах T5S3 и T3S5. В окрестностях 40˚ в.д. в поле температуры прослеживается Батумский антициклон, воспроизводимый в результатах численных экспериментов, по периферии которого наблюдаются циклонические вихри.

Рис. 3. Снимок температуры поверхности моря и модельные поля уровня моря 02.08.2016.

Далее со временем наибольшие отличия в поле уровня проявляются в структуре ОЧТ и динамике Севастопольского и Батумского антициклонов. В октябре–ноябре ширина ОЧТ в западной и южной частях моря в эксперименте T2S2 больше, следовательно, здесь круговорот менее интенсивен. Севастопольский антициклон стационируется на кромке Северо-западного шельфа, однако, если в эксперименте T2S2 он ослабевает, то по данным T3S5 и T5S3 наоборот интенсифицируется. Во всех экспериментах Батумский антициклон сопровождается циклоническим круговоротом, который наиболее интенсивен в T2S2. При этом на спутниковых снимках интенсивный циклон в этом районе не наблюдался.

Помимо карт поля уровня об усилении интенсивности циркуляции в экспериментах T3S5 и T5S3 также свидетельствует и изменение кинетической энергии течений. На рис. 4 представлен временной ход средней в верхнем модельном слое плотности кинетической энергии. Как видно из рис. 4, начиная с мая, кинетическая энергия, а следовательно, и скорость течений в верхнем слое для эксперимента T2S2 меньше, чем в двух других расчетах. Так как постановка численных экспериментов и граничные условия для всех расчетов одинаковые, то, по-видимому, использование аппроксимаций, обеспечивающих сохранение температуры и солености в степенях больше двух, приводит к более интенсивной циркуляции, что соответствует реальным измерениям скорости течений [Артамонов и др., 2018].

Рис. 4. Изменение со временем средней плотности кинетической энергии в верхнем слое.

Как было показано выше, значимые отличия в результатах моделирования температуры и солености наблюдаются в верхнем 300-метровом слое. Анализ пространственных распределений поля солености на различных горизонтах показал, что разница в интенсивности циркуляции проявляется начиная с глубин примерно 75 м с конца весны. В экспериментах T3S5 и T5S3 в центральной части моря наблюдаются области повышенной солености по сравнению с T2S2 (рис. 5а), что свидетельствует об усилении подъема вод. Такая картина прослеживается до конца осени. Из рис. 1б видно, что летом в глубоководной части моря находились буи ARGO № 6901831 и 6901900. Для восточной части моря сопоставление с данными буя 6901831, двигавшегося на северо-восток, демонстрирует уменьшение ошибки в результатах экспериментов T3S5 и T5S3 (рис. 5б), когда буй переместился севернее 43˚ с.ш. в область с наибольшей разницей между экспериментами (рис. 5а). Так, анализ разницы между модельными и натурными данными для лета 2016 г. показал, что величины mean S составили 1.36, 1.22, 1.25‰ для экспериментов T2S2, T3S5, T5S3 соответственно. Отклонение солености вдоль траектории буя 6901900 свидетельствует, что в западной части моря (рис. 1а) эксперимент T3S5 дает наименьшую ошибку (рис. 5в): величины mean S составили 1.62, 1.28, 1.57‰ для экспериментов T2S2, T3S5, T5S3 соответственно. Профили солености вдоль траекторий буев демонстрируют, что уменьшение ошибки моделирования солености в центральной части связано с более точным воспроизведением верхней границы основного галоклина в экспериментах T3S5 и T5S3 (рис. 5г, д).

Рис. 5. Пространственное распределение поля солености 10.07.2016 на горизонте 100 м (а). Отклонение модельной солености от измеренной на горизонте 100 м летом 2016 г. для буя 6901831 (б) и буя 6901900 (в). Профили солености с указанием даты и координаты по данным моделирования и буев 6901831 (г) и 6901900 (д).

Анализ полей температуры на различных горизонтах показал, что наиболее значимые отличия наблюдаются в летний сезон в верхнем 30-метровом слое. На рис. 6 представлены пространственные распределения температуры по результатам трех экспериментов летом 2016 г. и отклонения от данных измерений вдоль траектории буев ARGO на горизонте 15 м. Как видно из рис. 6а, в западной глубоководной части моря температура в эксперименте T2S2 на 1–2˚С выше, чем в T3S5 и T5S3. Далее, со временем, западная часть прогревается, и во всех экспериментах получены качественно и количественно близкие результаты. По данным буя ARGO 6901866, действовавшего летом 2016 г. к северо-востоку от Анатолийского побережья (рис. 1б), отклонения температуры от данных измерений ниже по результатам T3S5 и T5S3 (рис. 6б): величины mean Т равны 1.36, 0.96, 0.85˚С для экспериментов T2S2, T3S5, T5S3 соответственно. Для северной части моря (рейс № 87 НИС “Профессор Водяницкий”, рис. 1б) наблюдаются наибольшие отклонения модельной температуры от измеренной и здесь существенной разницы между результатами экспериментов не выявлено (рис. 6в). Большие отклонения значений температуры связаны с неточностью восстановления глубины верхнего перемешанного слоя (табл. 3), что обусловлено высокочастотной изменчивостью физических процессов, не воспроизводимой моделью МГИ, но зафиксированной за счет высокой частоты и плотности выполнения измерений. Из рис. 6г видно, что в натурных данных имеется несколько вихревых образований с размером около 10 км (в правой и левой частях полигона), тогда как в модельных полях воспроизводиться только вихрь в центре полигона к северу от Крыма.

Рис. 6. Пространственное распределение полей температуры на горизонте 15 м на 04.07.2016 (а). Отклонение модельной температуры от измеренной на горизонте 15 м летом 2016 г. для буя 6901866 (б) и рейса № 87 (в). Поле температуры на горизонте 2 м по данным рейса № 87 и численных экспериментов на 04.07.2016 (г).

Наиболее важные на наш взгляд различия, обнаруженные в результатах моделирования поля температуры, – это толщина сезонного термоклина и толщина ХПС. В центральной части моря в экспериментах T3S5 и T5S3 с середины июня по конец сентября слой скачка температуры и слой ХПС тоньше, чем в T2S2. В качестве примера на рис. 7а показаны вертикальные разрезы температуры вдоль 43˚ с.ш. в июле. Как видно в окрестностях горизонта 25 м наблюдается сгущение изотерм в расчетах T3S5 и T5S3, а толщина ХПС в среднем менее 50 м. Анализ натурных данных за 2016–2019 гг. [Морозов и др., 2020] показывает, что в Черном море наблюдается поднятие нижней границы ХПС. По сведениям, представленным в монографии [Иванов и др., 2011], максимальная величина среднего для летнего сезона вертикального градиента температуры, характеризующего слой скачка, равна –1.4˚С/м. Средние за июль оценки величины вертикального градиента модельной температуры в слое 20–25 м для глубоководной части моря (мористее изобаты 1500 м) составляют –1.03, –1.23, –1.34˚С/м для первого, второго и третьего экспериментов соответственно. Анализ профилей температуры по данным ARGO (рис. 7б, в) и глубоководных станций рейса № 87 НИС “Профессор Водяницкий” (рис. 7г) вблизи 43˚ с.ш. указывает на то, что эксперименты T3S5 и T5S3 качественно и количественно точнее воспроизводят форму и толщину ХПС. Как видно из рис. 7б–г результаты расчетов T3S5 и T5S3 находятся в лучшем соответствии с наблюдениями в слое 25–50 м.

Рис. 7. Вертикальный разрез поля температуры вдоль 43˚ с.ш. 30.07.2016 (а). Профили температуры по данным моделирования, буев Арго № 6901831 (б), № 6901832 (в) и рейса № 87 (г).

Сравнение результатов моделирования полей температуры и солености указывает на то, что при использовании уравнений (П.2) и (П.3) наблюдается интенсификация подъема вод в центральной части моря. Анализ пространственно-временной изменчивости поля вертикальной скорости показал, что в летний сезон при близких по величине значениях скорости увеличивается площадь зоны отрицательных значений, соответствующих в модели МГИ (ось z направлена вниз) подъему вод. На рис. 8а представлены средние за лето поля вертикальной скорости на горизонте 300 м. Видно, что в центральной и восточной частях моря области, обозначенные голубым цветом (w < 0), занимают большую площадь в экспериментах T3S5 и T5S3. Также обращает внимание более регулярный характер зон чередования положительных и отрицательных значений скорости. Оценки среднего за весь период интегрирования профиля вертикальной скорости для глубоководной части моря (мористее изобаты 1000 м) демонстрируют увеличение скорости подъема вод в верхнем 600-метровом слое для экспериментов T3S5 и T5S3 (рис. 8б). Максимум вертикальной скорости в экспериментах T3S5 и T5S3 расположен в слое 27.5 –35 м, тогда как в эксперименте T2S2 в слое 35–45 м. Результаты, полученные с применением новых схем (П.2) и (П.3), согласуются с оценками глубины максимума скорости подъема вод по данным реанализа Черного моря [Дорофеев и др., 2016].

Рис. 8. Среднее за лето поле вертикальной скорости на горизонте 300 м (а). Средний для глубоководной части моря профиль вертикальной скорости по результатам трех экспериментов (б).

Заключение

В работе представлены результаты сравнительного анализа гидрофизических полей Черного моря, реконструированных при использовании в численной модели циркуляции МГИ дискретных уравнений адвекции-диффузии тепла и соли, обладающих свойством сохранения температуры и солености в различных степенях. Рассмотрены три схемы аппроксимации адвективных членов, обеспечивающие сохранение температуры и солености в степенях от двух до пяти.

Результаты расчетов демонстрируют, что в верхнем 100-метровом слое величина отклонения модельных значений солености от измеренных меньше в экспериментах Т5S3 и Т3S5 по сравнению с экспериментом Т2S2, в котором используется традиционная аппроксимация нелинейных слагаемых. Наименьшие отклонения модельной солености от измеренной наблюдаются для юго-западной части моря – зоны с интенсивной струйной динамикой и слабой вихревой изменчивостью в зимний период. Для остальных районов в слое 30–100 м RMSE солености для экспериментов T3S5 и T5S3 меньше примерно на 20%, чем по данным T2S2.

Для верхнего слоя 0–30 м использование аппроксимаций, обеспечивающих наличие инвариантов степени больше двух, приводит к лучшим качественным и количественным результатам по температуре. Оценки показали, что RMSE температуры в экспериментах T3S5 и T5S3 уменьшается приблизительно на 10%. В то же время в летний сезон изменение схемы аппроксимации мало влияет на точность воспроизведения полей температуры в слое 30–100 м.

Анализ профилей солености вдоль траекторий буев демонстрирует, что уменьшение ошибки моделирования солености обусловлено более точным воспроизведением верхней границы основного галоклина. Сравнение профилей модельной температуры с данными ARGO и глубоководных станций НИС “Профессор Водяницкий“ показывает, что эксперименты T3S5 и T5S3 по сравнению с расчетом T2S2 качественно и количественно точнее воспроизводят форму и толщину ХПС. По результатам трех экспериментов наименьшие отклонения от данных наблюдений получены при использовании аппроксимаций, обеспечивающих сохранение температуры в третьей степени и солености в пятой степени.

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

Благодарности

Авторы выражают благодарность сотрудникам ФГБУН ФИЦ МГИ С.В. Станичному и А.А. Кубрякову за предоставленные данных альтиметрии, Е.А. Скрипалевой за предоставление данных съемок рейсов НИС “Профессор Водяницкий”. Работа выполнена при финансовой поддержке гранта РНФ № 23-27-00141.

Приложение

Введем разностную сетку для бассейна с неровным дном, которая описывается боксами с целочисленными значениями в его центрах i, j, k (i = i₁, …, i_N, j = j₁, …, j_M, k = 1, …, K_i,j) и на гранях – i + 1/2, j + 1/2, k + 1/2. Горизонтальные размеры боксов (h_x, h_y) постоянные, по вертикали используется неравномерная аппроксимация

$\left(h_z^k=z_{k+1/2}-z_{k-1/2},h_z^{k+1/2}=z_{k+1}-z_k\right).$

Разностные операторы имеют вид (для j, k – аналогично):

$\begin{array}{c}{\phi }_{i,j,k}^{-x}=\frac{{\phi }_{i+1/\mathrm{2,}j,k}+{\phi }_{i-1/\mathrm{2,}j,k}}{2},\\ {\delta }_x{\phi }_{i,j,k}=\frac{{\phi }_{i+1/\mathrm{2,}j,k}-{\phi }_{i-1/\mathrm{2,}j,k}}{h_x},\\ {\nabla }_{x,y}^2{\phi }_{i,j,k}={\delta }_x^2{\phi }_{i,j,k}+{\delta }_y^2{\phi }_{i,j,k}.\end{array}$

Тогда уравнения адвекции температуры и солености в точке (i, j, k) записываются в виде

$\begin{array}{c}\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)=\mathrm{0,}\end{array}$

$\begin{array}{c}\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)=\mathrm{0,}\end{array}$

Для температуры и солености в работе [Демышев, 2023] при условии сохранения T, S и $T^K,S^L$ на гранях бокса (полуцелые значения индексов) получены следующие соотношения по оси х (для j, k – аналогично):

$\begin{array}{c}{T}_{i+1/\mathrm{2,}j,k}=\frac{K-1}{K}\left(T_{i,j,k}+\frac{T_{i+\mathrm{1,}j,k}^{K-1}}{\Phi (T_{i+\mathrm{1,}j,k},T_{i,j,k})}\right),\\ S_{i+1/\mathrm{2,}j,k}=\frac{L-1}{L}\left(S_{i,j,k}+\frac{S_{i+\mathrm{1,}j,k}^{L-1}}{\Psi (S_{i+\mathrm{1,}j,k},S_{i,j,k})}\right),\end{array}$

где

$\begin{array}{l}\Phi \left(T_{i+\mathrm{1,}j,k},T_{i,j,k}\right)=\sum _{n=0}^{K-2}{T}_{i+\mathrm{1,}j,k}^{K-n-2}{T}_{i,j,k}^n,\\ \Psi \left(S_{i+\mathrm{1,}j,k},S_{i,j,k}\right)=\sum _{m=0}^{L-2}{S}_{i+\mathrm{1,}j,k}^{L-m-2}{S}_{i,j,k}^m.\end{array}$

При K = L = 2 следует известная аппроксимация (для j, k – аналогично):

$\begin{array}{l}{T}_{i+1/\mathrm{2,}j,k}={\overline{T_{i+1/\mathrm{2,}j,k}}}^x,\\ \text{}{S}_{i+1/\mathrm{2,}j,k}={\overline{S_{i+1/\mathrm{2,}j,k}}}^x.\end{array}$ (П.1)

При K = 5 и L = 3 имеем (для j, k – аналогично)

$T_{i+1/\mathrm{2,}j,k}=\frac{4}{5}\frac{\begin{array}{l}{T}_{i+\mathrm{1,}j,k}^4+T_{i+\mathrm{1,}j,k}^3{T}_{i,j,k}+\\ +T_{i+\mathrm{1,}j,k}^2{T}_{i,j,k}^2+T_{i+\mathrm{1,}j,k}^{}{T}_{i,j,k}^3+T_{i,j,k}^4\end{array}}{\begin{array}{l}{T}_{i+\mathrm{1,}j,k}^3+T_{i+\mathrm{1,}j,k}^2{T}_{i,j,k}+\\ +T_{i+\mathrm{1,}j,k}^{}{T}_{i,j,k}^2+T_{i,j,k}^3\end{array}},$ (П.2)

$S_{i+1/\mathrm{2,}j,k}=\frac{2}{3}\frac{S_{i+\mathrm{1,}j,k}^2+S_{i+\mathrm{1,}j,k}^{}{S}_{i,j,k}^{}+S_{i,j,k}^2}{S_{i+\mathrm{1,}j,k}^{}+S_{i,j,k}^{}},$

и, соответственно, при K = 3 и L = 5:

$S_{i+1/\mathrm{2,}j,k}=\frac{4}{5}\frac{\begin{array}{l}{S}_{i+\mathrm{1,}j,k}^4+S_{i+\mathrm{1,}j,k}^3{S}_{i,j,k}+\\ +S_{i+\mathrm{1,}j,k}^2{S}_{i,j,k}^2+S_{i+\mathrm{1,}j,k}^{}{S}_{i,j,k}^3+S_{i,j,k}^4\end{array}}{\begin{array}{l}{S}_{i+\mathrm{1,}j,k}^3+S_{i+\mathrm{1,}j,k}^2{S}_{i,j,k}+\\ +S_{i+\mathrm{1,}j,k}^{}{S}_{i,j,k}^2+S_{i,j,k}^3\end{array}},$ (П.3)

$T_{i+1/\mathrm{2,}j,k}=\frac{2}{3}\frac{T_{i+\mathrm{1,}j,k}^2+T_{i+\mathrm{1,}j,k}^{}{T}_{i,j,k}^{}+T_{i,j,k}^2}{T_{i+\mathrm{1,}j,k}^{}+T_{i,j,k}^{}}.$

При интегрировании по всей области аппроксимации (П.1), (П.2), (П.3) обеспечивают сохранение T, S и $T^K,S^L$ $\left(K\ge \mathrm{2,}L\ge 2\right)$, то есть выполняются соотношения

$\begin{array}{c}\sum _{i,j,k}\frac{d{T}_{i,j,k}^K}{dt}{h}_x{h}_y{h}_z+\sum _{i,j}\frac{d{\varsigma }_{i,j}}{dt}{T}_{i,j\mathrm{,1}/2}^K{h}_x{h}_y=\mathrm{0,}\\ \sum _{i,j,k}\frac{d{S}_{i,j,k}^L}{dt}{h}_x{h}_y{h}_z+\sum _{i,j}\frac{d{\varsigma }_{i,j}}{dt}{S}_{i,j\mathrm{,1}/2}^L{h}_x{h}_y=\mathrm{0,}\end{array}$

где ζ − приведенный уровень моря и эти выражения при отсутствии диффузии и внешних сил соответствуют следующим интегралам:

$\begin{array}{c}\frac{1}{V}\frac{\partial }{\partial t}\underset{-\zeta }{\overset{H}{\int }}\underset{{\Omega }_z}{\iint }\left(T^K\right)dxdydz=\mathrm{0,}\\ \frac{1}{V}\frac{\partial }{\partial t}\underset{-\zeta }{\overset{H}{\int }}\underset{{\Omega }_z}{\iint }\left(S^L\right)dxdydz=0.\end{array}$

Полученные аппроксимации требуют дополнительного анализа при |

作者简介

S. Demyshev

编辑信件的主要联系方式.
Email: demyshev@gmail.com
俄罗斯联邦, 299011, Sevastopol, st. Kapitanskaya, 2

O. Dymova

Email: olgdymova@mhi-ras.ru
俄罗斯联邦, 299011, Sevastopol, st. Kapitanskaya, 2

参考

补充文件

附件文件

动作

1. JATS XML

下载

2. Supplementary

下载 (120KB)

索引源数据

3. Fig. 1. Maps of the location of temperature and salinity measurement stations in winter (a) and summer (b) 2016 based on data from the R/V Professor Vodyanitsky and ARGO profiling buoys.

下载 (202KB)

索引源数据

4. Fig. 2. Snapshot of chlorophyll concentration and model sea level fields on 28.04.2016.

下载 (506KB)

索引源数据

5. Fig. 3. Sea surface temperature image and model sea level fields on 02.08.2016.

下载 (556KB)

索引源数据

6. Fig. 4. Changes over time in the average density of kinetic energy in the upper layer.

下载 (217KB)

索引源数据

7. Fig. 5. Spatial distribution of the salinity field on 10.07.2016 at the 100 m horizon (a). Deviation of the model salinity from that measured at the 100 m horizon in the summer of 2016 for buoy 6901831 (b) and buoy 6901900 (c). Salinity profiles with date and coordinates according to modeling data and buoys 6901831 (d) and 6901900 (d).

下载 (463KB)

索引源数据

8. Fig. 6. Spatial distribution of temperature fields at the 15 m horizon on 04.07.2016 (a). Deviation of the model temperature from that measured at the 15 m horizon in summer 2016 for buoy 6901866 (b) and cruise #87 (c). Temperature field at the 2 m horizon according to cruise #87 data and numerical experiments on 04.07.2016 (d).

下载 (616KB)

索引源数据

9. Fig. 7. Vertical section of the temperature field along 43˚ N on July 30, 2016 (a). Temperature profiles based on modeling data, Argo buoys No. 6901831 (b), No. 6901832 (c) and cruise No. 87 (d).

下载 (402KB)

索引源数据

10. Fig. 8. Average summer vertical velocity field at the 300 m horizon (a). Average vertical velocity profile for the deep-water part of the sea based on the results of three experiments (b).

下载 (260KB)

索引源数据

11. Table 1

下载 (42KB)

索引源数据

12. Table 1

下载 (48KB)

索引源数据

13. Table 1

下载 (43KB)

索引源数据

14. Table 1

下载 (42KB)

索引源数据

15. Table 2

下载 (46KB)

索引源数据

16. Table 2

下载 (41KB)

索引源数据

17. Table 2

下载 (40KB)

索引源数据

18. Table 2

下载 (38KB)

索引源数据

19. Table 3

下载 (40KB)

索引源数据

20. Table 3

下载 (39KB)

索引源数据

21. Table 3

下载 (39KB)

索引源数据

22. Table 3

下载 (36KB)

索引源数据

23. Table 4

下载 (47KB)

索引源数据

24. Table 4

下载 (43KB)

索引源数据

25. Table 4

下载 (45KB)

索引源数据

26. Table 4

下载 (47KB)

索引源数据