Океанология, 2022, T. 62, № 4, стр. 558-563

ВВЕДЕНИЕ

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

Данная проблема возникла, в частности, при численном моделировании термохалинных полей на полигоне ИО РАН в районе Геленджика на Черном море и их сопоставлении с результатами CTD-съемок в сентябре 2019 г. Изначально планировалось выполнение двух последовательных термохалинных съемок вод шельфа и верхней части континентального склона 14 и 15 сентября с пространственным покрытием приблизительно 15 × 15 км. Предполагалось, что перерыв между окончанием первой съемки и началом второй составит около 12 ч. При скорости адвективного переноса порядка 10 см/с и наличии устойчивой погоды на период работ имелась возможность постановки задачи краткосрочного прогноза параметров морской среды от первой съемки ко второй как задачи Коши.

По объективным причинам вторая съемка была выполнена примерно через 40 ч (16 сентября) после первой, что сделало практически невозможным прогноз в рамках решения задачи Коши, так как вода, находившаяся на полигоне в период первой съемки, покинула полигон.

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

Для качественного анализа прогноза будем использовать уровень моря, для которого имеется расчет уровня по данным второй съемки, и данные спутникового зондирования, приведенные в работе [3].

ПОСТАНОВКА ЗАДАЧИ И АНАЛИЗ РЕШЕНИЙ

Для решения задачи Коши используем систему осредненных по Рейнольдсу уравнений геофизической гидродинамики для ньютоновской жидкости в приближении Буссинеска:

Здесь $\frac{D}{{Dt}} \equiv \frac{\partial }{{\partial t}} + {\mathbf{U}} \cdot \nabla $, ${\mathbf{U}} = (u,{v},w)$ – вектор скорости течения, где $u$, ${v}$ – горизонтальные компоненты по осям $x$ и $y$ соответственно, а $w$ – скорость по вертикальной оси $z$; $f$ – параметр Кориолиса; $p$ – гидростатическое давление; ${{\rho }_{0}}$ – среднее (постоянное) значение плотности; $g$ – ускорение свободного падения; $T$, $S$ – температура и соленость морской воды. Операторы ${{Q}_{u}}$, ${{Q}_{{v}}}$, ${{Q}_{T}}$, ${{Q}_{S}}$ параметризуют процессы турбулентного перемешивания для компонент скорости и процессы турбулентной диффузии для тепла и соли; ${{Q}_{{mT}}}$, ${{Q}_{{mS}}}$ – операторы ассимиляции натурных глубоководных термохалинных измерений для, соответственно, тепла и соли. Параметр δ здесь принимает значение 0 при решении прямой задачи и равен 1 при решении обратной задачи по времени для нахождения начальных и граничных значений прямой задачи, как это приведено в работах [7, 8 , 10].

Система уравнений (1)–(7) дополняется следующими краевыми условиями:

– на свободной поверхности моря $z = - \zeta $: ${{D\zeta } \mathord{\left/ {\vphantom {{D\zeta } {Dt}}} \right. \kern-0em} {Dt}} = - w$, ${{\nu }_{T}}{{\partial u} \mathord{\left/ {\vphantom {{\partial u} {\partial z = {{ - {{\tau }_{x}}} \mathord{\left/ {\vphantom {{ - {{\tau }_{x}}} \rho }} \right. \kern-0em} \rho }}}} \right. \kern-0em} {\partial z = {{ - {{\tau }_{x}}} \mathord{\left/ {\vphantom {{ - {{\tau }_{x}}} \rho }} \right. \kern-0em} \rho }}}$, ${{\nu }_{T}}{{\partial {v}} \mathord{\left/ {\vphantom {{\partial {v}} {\partial z = {{ - {{\tau }_{y}}} \mathord{\left/ {\vphantom {{ - {{\tau }_{y}}} \rho }} \right. \kern-0em} \rho }}}} \right. \kern-0em} {\partial z = {{ - {{\tau }_{y}}} \mathord{\left/ {\vphantom {{ - {{\tau }_{y}}} \rho }} \right. \kern-0em} \rho }}}$, ${{\partial T} \mathord{\left/ {\vphantom {{\partial T} {\partial z}}} \right. \kern-0em} {\partial z}} = {{\partial S} \mathord{\left/ {\vphantom {{\partial S} {\partial z = }}} \right. \kern-0em} {\partial z = }}0$. Здесь ${{\nu }_{T}}$ – коэффициент турбулентной вязкости, ${{\tau }_{x}}$, ${{\tau }_{y}}$ – напряжения трения ветра (потоки импульса). Граничное условие непрерывности давления на свободной поверхности, $p = {{p}_{{atm}}}$, исключено из рассмотрения ввиду малости района расчета;

– на жидкой боковой цилиндрической поверхности горизонтальные компоненты скорости течения, температура и соленость находятся из условия излучения при решении прямой и обратной задач с учетом того, что каждая из точек втока при смене направления интегрирования по времени становится точкой вытока, см. [7, 8];

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

В начальный момент времени задавалось распределение температуры и солености, уровень моря; горизонтальные компоненты скорости течений $u$ и ${v}$ полагались равными нулю.

Прогноз выполнялся от начального момента времени, который соответствует последнему по времени решению для 1-й съемки. Интервал прогностического расчета от первой до второй съемки составлял чуть больше 40 ч. Результат прогноза сравнивается с данными расчета на основе постановки задачи (1)–(7) на период 2-й съемки.

На рис. 1а приведено начальное для прогноза поле уровня, полученное на основе обработки глубоководных термохалинных данных первой съемки с помощью технологии, приведенной в работах [7, 8]. Интервал изолиний уровня моря здесь и далее – 0.1 см. На рисунке 1б показано поле уровня в начальный момент для второй съемки, полученное на основе такой же обработки натурных данных второй съемки. Отметим, что данные первой съемки отражают только антициклонический вихрь, а второй съемки, – циклонический. Указанные решения качественно согласуются с данными, полученными на основе анализа спутниковой информации [3] (рис. 2). Результат прогноза на основе решения задачи Коши приведен на рисунке 1в. Сравнивая рисунки 1б (расчет уровня по данным 2-й съемки) и 1в (результат прогностического расчета как задачи Коши), видно, что расчет не является приемлемым.

Численный эксперимент показал, что необходим другой подход к решению задачи прогноза. Было решено воспользоваться методом, предложенным Г.П. Курбаткиным [4] и реализованным для других задач в работах [2, 9 ], где рассматриваются априорные неадиабатические источники как основной инструмент.

Расчет таких источников предполагает численный анализ предыстории процессов на основе натурных данных, сводящийся, по сути, к вычислению временны́х трендов измеренных априори параметров. В нашем случае имеются данные термохалинных измерений для первой съемки. Рассмотрим уравнения переноса импульса (1)–(2), тепла и солей (4)–(5). В правой части уравнений стоят неадибатические слагаемые ${{Q}_{u}}$, ${{Q}_{{v}}}$, ${{Q}_{T}}$, ${{Q}_{S}}$ которые параметризуют процессы турбулентного перемешивания в осредненной по Рейнольдсу модели. В численной модели используется замыкание Лаундера, Рисса и Роди [6, 12, 13], включающее уравнения для интенсивности турбулентной энергии и скорости ее диссипации. Модель хорошо известна и многократно использовалась в предыдущих работах [7, 8 , 10 ]. По положениям модели, используемые константы являются квазиуниверсальными. При этом процессы турбулентного обмена теплом и солью отражены в данных измерений. Вторые слагаемые, ${{Q}_{{mT}}}$ и ${{Q}_{{mS}}}$, описывают эффект ассимиляции данных глубоководных термохалинных измерений при нахождении краевых условий при интегрировании системы уравнений (1)–(7) назад по времени и имеют вид ньютоновских источников с нелинейными весовыми функциями.

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

Предполагаем, что реальная физика неадиабатических процессов будет отражена в слагаемых ${{Q}_{{mT}}}$ и ${{Q}_{{mS}}}$. Необходимо перейти от системы (1)– (7) к квазигиперболической системе уравнений для исключения из рассмотрения процессов, связанных с численной параметризацией процессов турбулентного перемешивания. Используя гипотезу Гульдберга-Моона [5, 11] о внутреннем трении, получим:

Здесь слагаемые $ku$, $k{v}$ – так называемое внутреннее трение, где коэффициент $k$ в выполненных расчетах полагался постоянным $k = {{{{{10}}^{{ - 2}}}} \mathord{\left/ {\vphantom {{{{{10}}^{{ - 2}}}} {\Delta t}}} \right. \kern-0em} {\Delta t}}$, т.е. данное слагаемое принято на два порядка меньшим, чем нестационарное слагаемое. В океанологии трение в такой форме использовал Стоммел в задаче о западной интенсификации течений [15], Нейман в работе [14], в метеорологии использование данного слагаемого можно найти в работе Гилла [1]. Введение внутреннего трения помогает также избежать проблем с устойчивостью численного решения системы (8)–(14) в отсутствие слагаемых, описывающих турбулентный обмен и диффузию.

Слагаемые в правых частях уравнений (11) и (12) вида $Q_{{mS}}^{*}$ и $Q_{{mT}}^{*}$ имеют смысл указанных выше априорных неадиабатических источников и содержат, в том числе, и процессы турбулентного перемешивания. Эти функции рассчитываются априори на первом этапе решения задачи прогноза путем вычисления левой части уравнений (11) и (12) на основе решения системы уравнений (8)– (14) с ассимиляцией натурных данных 1-й съемки. На этом этапе в правых частях уравнений (11) и (12) слагаемые $Q_{{mS}}^{*}$ и $Q_{{mT}}^{*}$ заменяются на операторы усвоения натурных данных $\delta {{Q}_{{mT}}}$, $\delta {{Q}_{{mS}}}$ и выполняется анализ на основе данных первой съемки, после чего и вычисляются неадиабатические источники $Q_{{mS}}^{*}$ и $Q_{{mT}}^{*}$.

Рассматриваются следующие краевые условия:

– на свободной поверхности моря $z = - \zeta $: ${{D\zeta } \mathord{\left/ {\vphantom {{D\zeta } {Dt}}} \right. \kern-0em} {Dt}} = - w$;

– на жидкой боковой цилиндрической поверхности компоненты $u$, ${v}$ и $T$, $S$ задаются на основе условий излучения при решении прямой и обратной задач, см. [7, 8]. При реализации прямой и обратной задач по времени каждая жидкая граничная точка втока становится точкой вытока и будет определяться, в том числе, на основе данных измерений, т.е. будет согласована с данными измерений;

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

В начальный момент времени полагается: $\zeta = 0$, $u = {v} = 0$, а распределения температуры $T$ и солености $S$ заданы.

Решение задачи прогноза состоит из двух этапов. На первом этапе решается задача анализа данных термохалинных измерений для системы уравнений (8)–(14), где в правых частях уравнений (11) и (12) переноса тепла и солей стоят слагаемые ${{Q}_{{mT}}}$, ${{Q}_{{mS}}}$, описывающие ассимиляцию натурных данных. Выполняется обработка данных термохалинных измерений для первой съемки 14 сентября и расчет неадиабатических источников $Q_{{mS}}^{*}$ и $Q_{{mT}}^{*}$. Для этого рассчитываются полные производные в левых частях уравнений (11) и (12) на основе полученного решения. Нестационарные слагаемые вычислялись как разность между решением в конечный и начальный моменты времени для первой съемки. Расчет адвективного слагаемого в уравнениях (11) и (12) выполнялся на основе поля течений, полученного путем его осреднения за период расчета 1-й съемки. Непосредственно для описания процесса адвективного трехмерного переноса тепла и соли использовалась численная схема Годунова, предназначенная для описания негладких решений с разрывами производных.

После вычисления неадиабатических слагаемых $Q_{{mS}}^{*}$ и $Q_{{mT}}^{*}$, последние подставлялись в систему (8)–(14) для выполнения прогноза от решения в конечный момент первой съемки 14 сентября на период до конца второй съемки 16 сентября.

Идея прогноза в данном случае базируется на предположении, что временнáя тенденция процессов, отраженная в априорных источниках $Q_{{mS}}^{*}$ и $Q_{{mT}}^{*}$, сохранит свое значение в течение периода прогноза.

На рис. 3а показано прогностическое поле уровенной поверхности, полученное в результате указанного расчета. Приведенное решение качественно согласуется с рис. 1б, полученным на основе обработки с ассимиляцией натурных данных второй съемки.

При изложенной постановке расчета неадиабатического источника на основе системы (8)– (14), последний содержит изменчивость полей за счет эффекта усвоения измерений температуры и солености. Если же добавить априори к указанному эффекту изменчивость за счет адвективного, уже адиабатического процесса, то есть сложить тенденции в правой части уравнений (11)–(12), то получим решение для уровенной поверхности, приведенное на рис. 3б. Физически данное решение лучше согласуется с устойчивой схемой переноса вдоль побережья с востока на запад, в соответствии с циклоническим характером крупномасштабной циркуляции вод в восточной части моря. Кроме того, циклонический вихрь движется в последнем случае вдоль характеристики ${f \mathord{\left/ {\vphantom {f h}} \right. \kern-0em} h}$, прижимаясь ближе к берегу. Здесь $f$ – параметр Кориолиса, $h$ – глубина моря. Данная характеристика следует из уравнения квазигеострофического потенциального вихря. Однако, если считать за истину данные обработки непосредственно натурных данных второй съемки (рис. 1б), то трудно отдать предпочтение одному из вариантов, показанных на рис. 3. Здесь оба варианта прогноза приведены единственно для того, чтобы указать на важность правильных представлений о физике процессов в исследуемом районе, что необходимо для адекватной формулировки задачи прогноза.

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

ЗАКЛЮЧЕНИЕ

Использование метода неадиабатических источников при наличии данных натурных измерений может оказаться полезным при решении ряда океанологических задач, в том числе связанных с процессами взаимодействия океана и атмосферы. Например, ранее такой подход был использован для нахождения потоков тепла на поверхности океана в районе Гольфстрима в зимний период [2, 9 ]. Задачи прогноза, особенно среднесрочного и тем более долгосрочного, также могут быть сформулированы с учетом “технологии” неадиабатических источников, т.е. с использованием априорных временны́х трендов искомых характеристик.

В работе удалось применить данный подход к задаче краткосрочного прогноза. Показана возможность использования “технологии” неадиабатических источников для решения задач прогноза в том случае, когда неадиабатические временны́е тренды существенны в формировании гидрофизических полей на рассматриваемом интервале времени. Эта возможность оказалась связанной с особенностью физики процессов в рассмотренном регионе, поскольку эти особенности отражены в ассимилируемых натурных данных.

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

Источники финансирования: Работа выполнена в рамках темы госзадания № 0128-2021-0002 и при поддержке гранта РНФ № 18-11-00163 “Разработка иерархии математических моделей нового поколения для решения задач вычислительной океанологии на основе гиперболической декомпозиции и балансно-характеристического подхода”.