Журнал вычислительной математики и математической физики, 2020, T. 60, № 11, стр. 1915-1932

ВВЕДЕНИЕ

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

Существуют различные приближения, которые используются для задания граничных условий на внешних жидких границах. Одним из способов учета граничных условий можно считать проведение расчета по всей акватории Мирового океана – на грубой сетке вне рассматриваемой акватории, на более мелкой – внутри: при этом существенно возрастают вычислительные затраты, однако отсутствуют внешние жидкие границы. Иногда результаты расчетов по всей акватории Мирового океана на грубой сетке используются для задания граничных условий на внешней жидкой границе. Развитием данных идей является метод вложенных сеток – с обратной связью, с неполной обратной связью и без обратной связи [1], [2]. В некоторых региональных моделях (в частности, в модели ROMS) используются адаптивные граничные условия [3]. Еще одним распространенным приемом является использование осредненных данных о потоках через внешнюю жидкую границу [4] или задание потоков с помощью данных наблюдений [5]. Также перспективным методом для решения проблемы постановки граничных условий на внешней жидкой границе можно считать ассимиляцию данных. Ассимиляция данных может использоваться для коррекции результатов расчетов по данным измерений. Обзор методов ассимиляции данных в региональных моделях приведен в [6].

В задачах математического моделирования водных областей часто возникает необходимость более детального рассмотрения некоторой подобласти. Разделяя исходную область путем введения внутренней жидкой границы, можно также добиться некоторого упрощения процедур численного решения задачи в подобластях. Применение методов разделения области позволяет сводить процесс решения задачи в исходной области к поочередному (зависимому только в местах пересечения) решению задач в подобластях, возможно, имеющих более простую форму, а также применять сетки разного масштаба (т.е. характерный масштаб пространственных сеток в подобластях может отличаться). К основным направлениям развития методов разделения области можно отнести построение новых алгоритмов для различных классов задач и эффективную их реализацию на многопроцессорных вычислительных системах. В последнее время методы разделения области стали активно применяться при математическом моделировании процессов в океане. Так, для моделей, основанных на примитивных уравнениях океана, использовалась разновидность метода Шварца [7], для моделирования прибрежных потоков [8] были предложены методы соединения различных моделей, основанные на методах декомпозиции области. Некоторые специальные задачи совмещения моделей решаются также с использованием разновидностей методов разделения области в [9], [10]. Более того, методы разделения области получили распространение в задачах вариационной ассимиляции данных наблюдений [11], [12]. Хотя большинство работ, касающихся совместного применения методов разделения области и вариационной ассимиляции данных, посвящено созданию параллельных алгоритмов [11], [13], [14], методы разделения области могут также быть применимы и для случаев, в которых необходимо выделить одну из подобластей с данными наблюдений. Для таких задач процедура вариационной ассимиляции данных будет применяться только на подобласти с данными наблюдений. Приложение методов разделения области к задачам вариационной ассимиляции данных – относительно новое направление в математическом моделировании океанов и морей.

Рассматривается упрощенная модель гидродинамики открытой акватории, основанная на линеаризованной системе уравнений мелкой воды. Данная система является подзадачей метода расщепления по физическим процессам в нелинейной модели гидротермодинамики открытой акватории [15]. Модель основана на примитивных уравнениях, записанных в $\sigma $-координатах со свободной поверхностью в приближении гидростатики и Буссинеска [16]. Для аппроксимации модели по времени используется метод расщепления [17]. Для этого весь временной интервал делится на подынтервалы; на каждом подынтервале решаются следующие подзадачи, названные шагами метода расщепления.

Шаг 1. Задача о распространении тепла.

Шаг 2. Задача конвекции-диффузии для солености.

Шаг 3. Задача адаптации гидрологических полей, которая решается в три этапа: а) расчет плотности и нахождение поправок к скорости; б) решение задачи бароклинной адаптации; в) решение задачи баротропной адаптации (линеаризованная система уравнений мелкой воды).

Шаги метода расщепления сформулированы в [15]. Эта работа продолжает работы авторов по вариационной ассимиляции данных [15], [18]–[20], разделению области [21], [22] и их совместному применению [12], [23].

В данной работе рассматривается шаг 3в). Предполагается, что имеются данные о баротропной составляющей скорости на внешней жидкой границе, формулируется задача вариационной ассимиляции данных. Существует несколько способов получить эти данные. Первый (комбинированный) подход основан на использовании данных наблюдений со спутников и дрейфующих буев [24], [25]. По данным с буев можно приближенно оценить осредненное по времени трехмерное поле плотности [24], а по данным спутниковой альтиметрии вычисляется баротропная поправка к скорости. Путем интегрирования по глубине формул для расчета скорости геострофического течения можно оценить баротропные скорости [5]. Такой подход в применении к проблеме жидких границ имеет недостатки: оценки для скоростей получаются грубыми, и их можно получить только там, где проводятся частые измерения с буев и где проходят треки спутников, оснащенных альтиметрическим оборудованием. В случае, если имеются только данные спутниковой альтиметрии, для оценки относительных изменений баротропной скорости можно использовать баротропную поправку [5]. Если трек спутника проходит непосредственно по жидкой границе, можно оценить баротропную поправку к составляющей скорости, перпендикулярной этому треку; как будет показано в статье, этих данных достаточно для формулировки алгоритма ассимиляции. Второй подход основан на использовании результатов расчетов по модели большей акватории. В данном случае алгоритм ассимиляции может давать лучшие результаты по сравнению с непосредственной подстановкой данных в граничные условия. Это связано с тем, что региональные модели гидротермодинамики с более мелкой сеткой могут воспроизводить те физические явления, которые не воспроизводятся в глобальных моделях [1], в этом случае жесткое задание граничных условий по результатам расчетов глобальной модели может оказаться неэффективным. Необходимо дальнейшее исследование эффективности алгоритма ассимиляции данных на жидкой границе в рамках данного подхода, а также сравнение с известными алгоритмами вложенных сеток. В настоящей работе к задаче вариационной ассимиляции данных [26] применяется метод разделения области, сформулированный в [21]. Для этого исходная область разделяется на 2 подобласти путем введения внутренней жидкой границы. Cтавится обратная задача о восстановлении граничных функций на внешних и внутренних жидких границах, приводится доказательство однозначной и плотной разрешимости задачи, а также формулируется итерационный алгоритм ее решения. Представлены результаты численных экспериментов по использованию алгоритма в модельной задаче. Данный алгоритм может быть применен в модели гидротермодинамики морской акватории, основанной на методе расщепления. В таком случае алгоритм следует использовать на шаге метода расщепления, соответствующем задаче баротропной адаптации, на каждом временном слое модели. Общий подход, который предложен в работе, может быть расширен, однако для каждой конкретной задачи граничные условия на внешних и внутренних жидких границах следует формулировать в зависимости от выбранных упрощений, проводить соответствующее теоретическое и численное исследование, а также подбирать параметры алгоритма, влияющие на скорость сходимости и точность, по результатам исследований. Отметим также, что в реальных задачах необходимо детально рассматривать вопрос о доступности данных наблюдений на жидкой границе.

1. ПОСТАНОВКА ЗАДАЧИ

1.1. Основные обозначения и приближения, в которых записываются системы уравнений гидротермодинамики морской акватории, приведены в [15], [26]. В данной работе будем придерживаться тех же обозначений.

Рассмотрим географическую (– геодезическую) систему координат ($\lambda ,\theta ,r)$, где $\lambda \in [0,2\pi ]$ – географическая долгота, увеличивающаяся с запада на восток, $\theta \in [ - \pi {\text{/}}2,\pi {\text{/}}2]$ – географическая широта, растущая с юга на север, $r$ – расстояние точки от центра Земли, средний радиус которой принимается равным ${{R}_{E}}$. Часто вместо $r$ удобно вводить координату $z = {{R}_{E}} - r$ оси $Oz$, направленной по нормали от поверхности сферы ${{S}_{R}}$ радиуса ${{R}_{E}}$ к еe центру, т.е. по направлению силы тяжести. Единичные векторы в $\lambda $-, $\theta $- и $z$-направлениях обозначим соответственно через ${{e}_{\lambda }},\;{{e}_{\theta }},\;{{e}_{z}}$. Тогда вектор скорости в океане записывается в виде

Обозначим через $\Omega $ многосвязное многообразие на сфере ${{S}_{R}}$, которую будем называть также поверхностью отсчета [15]. Далее рассматривается случай, когда $\Omega $ не включает точки полюсов. Поверхность моря будем задавать уравнением $z = \xi (\lambda ,\theta ,t)$, где ($\lambda $,$\theta $,${{R}_{E}}) \in \Omega $, а $t$ – временная переменная, $t \in $[0, $T$], $T$<$\infty $. Функцию рельефа дна определим как $z = H(\lambda ,\theta )$ при ($\lambda $,$\theta $,${{R}_{E}}) \in \Omega $, где $H(\lambda ,\theta )$ > 0. В этих обозначениях полная глубина моря выражается следующим образом: ${{H}_{{\text{F}}}} = H - \xi $. Предположим дополнительно, что существует $\varepsilon > 0$ такое, что $H \geqslant \varepsilon $.

Введeм следующие дифференциальные операции градиента, дивергенции и полной производной в сферической системе координат (оставив за этими операциями известные обозначения из векторного анализа):

где $m = 1{\text{/}}({{R}_{E}}cos\theta )$, $n = 1{\text{/}}{{R}_{E}}$. Введем следующие обозначения для осредненных по глубине скоростей:

Обозначим через ${\mathbf{n}}$ единичный вектор внешней нормали к $\partial \Omega $. Пусть $l$ – параметр Кориолиса, $g$ – ускорение свободного падения. Рассмотрим линеаризованную задачу мелкой воды:

где ${\mathbf{f}} = {{({{f}_{1}},{{f}_{2}})}^{{\text{T}}}}$ – заданная функция, ${{R}_{f}}$ – коэффициент трения. Решая задачи на $({{t}_{{j - 1}}},{{t}_{j}})$, полагаем функции ${{{\mathbf{U}}}_{{j - 1}}}$, ${{\xi }_{{j - 1}}}$ известными. Подробное описание обозначений и приближений, используемых в работе, можно найти,например, в [15]. В настоящей работе рассматривается граничное условие на внешней жидкой границе вида [27]: где ${{m}_{{{\text{op}}}}}$ – характеристическая функция внешней жидкой (открытой) границы ${{\Gamma }_{{{\text{op}}}}}$, т.е. ${{m}_{{{\text{op}}}}} = 1$ если $(\lambda ,\theta ,{{R}_{E}}) \in {{\Gamma }_{{{\text{op}}}}}$, в противном случае ${{m}_{{{\text{op}}}}} = 0$. Заметим, что на твердой части границы (${{m}_{{{\text{op}}}}} = 0$) граничное условие принимает вид условия непротекания. Граничная функция ${{d}_{s}}$ будет рассматриваться ниже как дополнительная неизвестная. Пусть для определенности ${{d}_{s}} = 0$ на $\partial \Omega {\backslash }{{\Gamma }_{{{\text{op}}}}}$.

Весь интервал по времени $[0,T]$ разбивается на подынтервалы $({{t}_{{j - 1}}},{{t}_{j}})$, $\Delta t = {{t}_{j}} - {{t}_{{j - 1}}}$, ${{t}_{0}} = 0$, ${{t}_{J}} = T$, $j = 1, \ldots ,J$. Запишем неявную схему аппроксимации задачи (1), (2) на интервале $({{t}_{{j - 1}}},{{t}_{j}})$, по времени:

В дальнейшем будем рассматривать задачу (3) на интервале времени $({{t}_{{j - 1}}},{{t}_{j}})$ в полудискретном виде, поэтому для удобства будем опускать индексы $j$, т.е.: ${\mathbf{U}} \equiv {{{\mathbf{U}}}^{j}}$, $\xi \equiv {{\xi }^{j}}$, … Также введем вектор-функцию Тогда задача (3) запишется в виде

1.2. Пусть граница $\partial \Omega $ области $\Omega $ – кусочно-гладкая класса ${{\mathbb{C}}^{2}}$ и удовлетворяет условию Липшица. Будем считать, что кривая (простая дуга) ${{\Gamma }_{{{\text{in}}}}} \subset \Omega $ класса ${{\mathbb{C}}^{2}}$ делит область $\Omega $ на две подобласти ${{\Omega }_{1}}$ и ${{\Omega }_{2}}$ (см. фиг. 1), $\Omega \equiv {{\Gamma }_{{{\text{in}}}}} \cup {{\Omega }_{1}} \cup {{\Omega }_{2}}$, ${{\overline \Gamma }_{{{\text{in}}}}} \equiv \overline {(\partial {{\Omega }_{1}} \cap \partial {{\Omega }_{2}})} $, ${{\Gamma }_{{{\text{in}}}}} \subset \partial {{\Omega }_{1}}$ – открытое множество в $\partial {{\Omega }_{1}}$, ${{\Gamma }_{{{\text{in}}}}} \subset \partial {{\Omega }_{2}}$. Пусть границы $\partial {{\Omega }_{1}}$ и $\partial {{\Omega }_{2}}$ – кусочно-гладкие класса ${{\mathbb{C}}^{2}}$, липшицевы. Через $\partial {{\Omega }_{c}}$ обозначим твердую границу, т.е. $\partial {{\Omega }_{c}} \equiv \partial \Omega {\backslash }{{\Gamma }_{{{\text{op}}}}}$. Будем считать, что ${{\Gamma }_{{{\text{op}}}}}$ (внешняя жидкая граница) не пересекается с ${{\Gamma }_{{{\text{in}}}}}$ (внутренней жидкой границей).

Будем обозначать индексом $i = 1,2$ решение в каждой из подобластей ($i = 1,2$). Задача (3) может быть записана для функций ${{{\mathbf{U}}}^{{(1)}}},\;{{\xi }^{{(1)}}}$ в подобласти ${{\Omega }_{1}}$, ${{{\mathbf{U}}}^{{(2)}}}$, ${{\xi }^{{(2)}}}$ в подобласти ${{\Omega }_{2}}$. Пусть ${{{\mathbf{n}}}_{i}}$ – внешняя нормаль к границе $\partial {{\Omega }_{i}}$, $i = 1,2$. Потребуем выполнения условий сшивки на внутренней жидкой границе [21]:

Введем дополнительную неизвестную функцию $\nu $ так, что выполнено условие сшивки (5):

Мы получаем следующий вид граничных условий на границе области ${{\Omega }_{1}}$: а на границе области ${{\Omega }_{2}}$:

Пусть имеются данные наблюдений ${{I}_{{{\text{obs}}}}} = (H{\mathbf{U}} \cdot {\mathbf{n}})$ вдоль ${{\Gamma }_{{{\text{op}}}}}$, прошедшие предварительную обработку. Такие данные могут быть получены из наблюдений [5], а могут быть оценены по результатам расчетов по глобальной модели с более грубой сеткой. И в том, и в другом случае эти данные содержат ошибки, поэтому их непосредственное использование в качестве граничных условий не является желательным. Введем дополнительное условие (условие замыкания):

Из (4) в каждой подобласти ($i = 1,2$) мы получаем

Введем пространство $\mathbb{H}_{0}^{{(i)}}$ вектор-функций ${{{\mathbf{U}}}^{{(i)}}} \in {{({{\mathbb{L}}_{2}}({{\Omega }_{i}}))}^{2}}$ со скалярным произведением:

Введем пространство ${{\mathbb{X}}^{{(i)}}} = \{ {\mathbf{U}} \in {{({{\mathbb{L}}_{2}}({{\Omega }_{i}}))}^{2}},\;{{\left\| {div(H{\mathbf{U}})} \right\|}_{{{{\mathbb{L}}_{2}}({{\Omega }_{i}})}}} < \infty ,\;{\mathbf{U}} \cdot {\mathbf{n}} = 0\;{\text{на}}\;(\partial {{\Omega }_{c}} \cap \partial {{\Omega }_{i}})\} $. Обозначим через ${{\mathbb{W}}^{{(i)}}}$ пространство вектор-функций ${{\Phi }^{{(i)}}} \in {{\mathbb{X}}^{{(i)}}} \times {{\mathbb{L}}_{2}}({{\Omega }_{i}})$. Пусть ${{I}_{{{\text{obs}}}}} \in \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{op}}}}})$, ${\mathbf{\tilde {f}}} \in {{({{\mathbb{L}}_{2}}(\Omega ))}^{2}}$, ${{\tilde {f}}_{3}} \in {{\mathbb{L}}_{2}}(\Omega )$, $H \in {{\mathbb{C}}^{1}}(\bar {\Omega })$. Обозначим через ${{\mathbb{V}}_{1}}(\Gamma )$, ${{\mathbb{V}}_{2}}(\Gamma )$ ($\Gamma \in \partial \Omega $) пространства функций из ${{\mathbb{L}}_{2}}(\Gamma )$ со “взвешенными” скалярными произведениями: Введем гильбертово пространство ${{\mathbb{H}}_{c}}$ вектор-функций $\omega = {{({{d}_{s}},\nu )}^{{\text{T}}}}$, ${{d}_{s}} \in {{\mathbb{V}}_{1}}({{\Gamma }_{{{\text{op}}}}})$, $\nu \in {{\mathbb{V}}_{1}}({{\Gamma }_{{{\text{in}}}}})$, с нормой: где $\gamma = {\text{const}} > 0$. Также введем пространство ${{\mathbb{H}}_{{{\text{ob}}}}}$ вектор-функций $\psi = {{({{\psi }_{1}},{{\psi }_{2}})}^{{\text{T}}}}$, ${{\psi }_{1}} \in {{\mathbb{V}}_{2}}({{\Gamma }_{{{\text{op}}}}})$, ${{\psi }_{2}} \in {{\mathbb{V}}_{2}}({{\Gamma }_{{{\text{in}}}}})$, с нормой Коэффициент $\gamma $ может быть произвольным положительным числом. Как будет показано в разд. 6 настоящей работы, введение $\gamma $ позволяет увеличить скорость сходимости алгоритма решения обратной задачи.

Умножим первое уравнение в (8) на некоторую вектор-функцию ${{{\mathbf{\hat {U}}}}^{{(i)}}} = {{({{\hat {u}}^{{(i)}}},{{\hat {v}}^{{(i)}}})}^{{\text{T}}}} \in {{\mathbb{X}}^{{(i)}}}$ скалярно в $\mathbb{H}_{0}^{{(i)}}$ и проинтегрируем полученное выражение по частям с учетом граничных условий. Получим

где Сделаем замечание относительно введенных билинейных форм ${{b}_{{{\text{op}}}}}({{d}_{s}},{{\hat {\Phi }}^{{(1)}}})$, ${{b}_{i}}(\nu ,{{\hat {\Phi }}^{{(i)}}})$. Как известно [28], след $H{{{\mathbf{\hat {U}}}}^{{(i)}}} \cdot {{{\mathbf{n}}}_{i}}$ вектор-функции ${{{\mathbf{\hat {U}}}}^{{(i)}}} \in {{\mathbb{X}}^{{(i)}}}$ на $\partial {{\Omega }_{i}}$ принадлежит пространству $\mathbb{W}_{2}^{{ - 1/2}}(\partial {{\Omega }_{i}})$, следовательно, необходимо потребовать ${{d}_{s}} \in \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{op}}}}})$, $\nu \in \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{in}}}}})$, чтобы интегралы в определении данных билинейных форм имели смысл.

Для того, чтобы сформулировать обратную задачу о восстановлении граничных функций на внешних и внутренних жидких границах, провести теоретическое исследование и записать алгоритм ее решения, исключим ${{\xi }^{{(i)}}}$ из системы (8). Из второго уравнения в (8) получим

Подставив это представление в (9), (10), получим где

Сформулируем обратную задачу о восстановлении граничных функций на внешних и внутренних жидких границах: найти неизвестные граничные функции ${{d}_{s}} \in \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{op}}}}})$, $\nu \in \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{in}}}}})$ и вектор-функции решений в подобластях ${{\Phi }^{{(i)}}} = {{({{{\text{U}}}^{{(i)}}},{{\xi }^{{(i)}}})}^{{\text{T}}}} \in {{\mathbb{W}}^{{(i)}}}$ ($i = 1,2$), которые удовлетворяют уравнениям (12), (13) и дополнительным условиям (11) ($i = 1,2$), (6), (7).

Билинейные формы ${{\tilde {a}}_{i}}({{{\mathbf{U}}}^{{(i)}}},{{{\mathbf{\hat {U}}}}^{{(i)}}})$, $i = 1,2$, ограничены и положительно определены для функций ${{{\mathbf{U}}}^{{(i)}}},\;{{{\mathbf{\hat {U}}}}^{{(i)}}} \in {{\mathbb{X}}^{{(i)}}}$. Равенства (12), (13) могут быть сформулированы в следующей форме [29]:

где операторы ${{L}_{i}}:{{\mathbb{X}}^{{(i)}}} \to ({{\mathbb{X}}^{{(i)}}}){\kern 1pt} '$, ${{B}_{i}}:{{\mathbb{V}}_{1}}({{\Gamma }_{{{\text{in}}}}}) \to ({{\mathbb{X}}^{{(i)}}}){\kern 1pt} '$, $D({{B}_{i}}) \equiv \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{in}}}}})$, ${{B}_{{{\text{op}}}}}{\kern 1pt} :\;{{\mathbb{V}}_{1}}({{\Gamma }_{{{\text{op}}}}}) \to ({{\mathbb{X}}^{{(i)}}}){\kern 1pt} '$, $D({{B}_{{{\text{op}}}}}) \equiv \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{op}}}}})$ (пространство $({{\mathbb{X}}^{{(i)}}}){\kern 1pt} '$ является сопряженным к ${{\mathbb{X}}^{{(i)}}}$) вводятся с использованием билинейных форм ${{\tilde {a}}_{i}}({{{\mathbf{U}}}^{{(1)}}},{{{\mathbf{\hat {U}}}}^{{(1)}}})$, ${{\tilde {b}}_{i}}(\nu ,{{{\mathbf{\hat {U}}}}^{{(i)}}})$, ${{\tilde {b}}_{{{\text{op}}}}}({{d}_{s}},{{{\mathbf{\hat {U}}}}^{{(i)}}})$, $i = 1,2$, соответственно [30], ${{\hat {f}}_{i}} \in ({{\mathbb{X}}^{{(i)}}}){\kern 1pt} '$. Известно [30], что операторы ${{L}_{i}}$ являются изоморфизмами из ${{\mathbb{X}}^{{(i)}}}$ на $({{\mathbb{X}}^{{(i)}}}){\kern 1pt} '$, поэтому существуют ограниченные обратные операторы $L_{i}^{{ - 1}}$. Можно ввести сопряженные операторы $L_{i}^{*}$ так, что будут выполняться следующие соотношения ($i = 1,2$): Заметим, что $({{L}_{i}}{{{\mathbf{U}}}^{{(i)}}}) \in ({{\mathbb{X}}^{{(i)}}}){\kern 1pt} '$, ${{{\mathbf{U}}}^{{(i)}}} \in {{\mathbb{X}}^{{(i)}}}$ и $( \cdot , \cdot )$ обозначает их скалярное произведение, т.е. действие функционала $({{L}_{i}}{{{\mathbf{U}}}^{{(i)}}})$ на элементе ${{{\mathbf{\hat {U}}}}^{{(i)}}}$. Запишем соотношения, определяющие операторы $B_{{{\text{op}}}}^{*}$, $B_{i}^{*}$:

Введем операторы ${{C}_{{{\text{obs}}}}}:{{\mathbb{X}}^{{(1)}}} \to {{\mathbb{V}}_{2}}({{\Gamma }_{{{\text{op}}}}})$, ${{C}_{i}}:{{\mathbb{X}}^{{(i)}}} \to {{\mathbb{V}}_{2}}({{\Gamma }_{{{\text{in}}}}})$, $i = 1,2$:

Заметим, что для гладких вектор-функций ${{{\mathbf{U}}}^{{(i)}}}$ действие операторов ${{C}_{{{\text{obs}}}}}$, ${{C}_{i}}$ может быть выражено следующим образом: ${{C}_{{{\text{obs}}}}}{{{\mathbf{U}}}^{{(1)}}} = H{{{\mathbf{U}}}^{{(1)}}} \cdot {{\left. {{{{\mathbf{n}}}_{1}}} \right|}_{{{{\Gamma }_{{{\text{op}}}}}}}}$, ${{C}_{i}}{{{\mathbf{U}}}^{{(i)}}} = H{{{\mathbf{U}}}^{{(i)}}} \cdot {{\left. {{{{\mathbf{n}}}_{i}}} \right|}_{{{{\Gamma }_{{{\text{in}}}}}}}}$. Нетрудно показать, что области определения $D(C_{{{\text{obs}}}}^{*})$, $D(C_{i}^{*})$ сопряженных операторов включают в себя пространства $\mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{op}}}}})$, $\mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{in}}}}})$ соответственно.

Условия замыкания (6), (7) можно записать в виде

Сформулируем операторную постановку задачи: найти ${{{\mathbf{U}}}^{{(i)}}} \in {{X}^{{(i)}}}$, ${{d}_{s}} \in \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{op}}}}})$, $\nu \in \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{in}}}}})$, удовлетворяющие (16), (17) и (20), (21).

Если исключить ${{{\mathbf{U}}}^{{(i)}}}$, то система уравнений (16), (17), (20), (21) может быть записана в виде одного операторного уравнения:

где Таким образом, исходная задача была сведена к единому операторному уравнению (22) с оператором $A:{{\mathbb{H}}_{c}} \to {{\mathbb{H}}_{{{\text{ob}}}}}$, $D(A) \equiv \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{op}}}}}) \times \mathbb{W}_{2}^{{1/2}}({{\Gamma }_{{{\text{in}}}}})$, для которого применимы многие результаты общей теории операторных уравнений [31]. При дальнейшем рассмотрении задачи используются методы оптимального управления и сопряженных уравнений, рассмотренные в [29], а также методы решения обратных и некорректно поставленных задач [32], [33].

1.3. Включим задачу (22) в семейство задач оптимального управления [29]: найти $\omega \in D(A)$, доставляющую минимум функционалу ${{J}_{\alpha }}(\omega )$:

В терминах исходных функций условие (23) имеет вид Заметим, что функционал ${{J}_{\alpha }}(\omega )$ имеет вид функционала А.Н. Тихонова [32]–[34], где $\alpha $ – параметр регуляризации, поэтому будем далее называть $\alpha $ параметром регуляризации.

2. ЕДИНСТВЕННОСТЬ РЕШЕНИЯ ЗАДАЧИ

Покажем, что для сформулированной задачи выполняется условие единственности решений.

Теорема 1. Уравнение (22) однозначно разрешимо.

Доказательство. Покажем, что ядро оператора $A$ состоит только из нулевого элемента (${\text{ker}}(A) = \{ 0\} $). Пусть существует функция $\omega \ne 0$, удовлетворяющая однородному уравнению $A\omega = 0$. Это означает, что существует ненулевой вектор $\omega = {{({{d}_{s}},\nu )}^{{\text{T}}}}$ такой, что:

где $i = 1,2$. Заметим, что функции $\xi = \{ {{\xi }^{{(i)}}}\;{\text{в}}\;{{\Omega }_{i}},\;i = 1,2\} $, ${\mathbf{U}} = \{ {{{\mathbf{U}}}^{{(i)}}}\;{\text{в}}\;{{\Omega }_{i}},\;i = 1,2\} $ будут в данном случае являться решениями однородной задачи в области $\Omega $ [21]. Отсюда получаем, что ${{\xi }^{{(i)}}} = 0$, ${{{\mathbf{U}}}^{{(i)}}} = 0$ и ${{d}_{s}} = 0$, $v = 0$. Значит, $u = 0$ и ${\text{ker}}(A) = \{ 0\} $, что и означает однозначную разрешимость задачи (22). Таким образом, для рассматриваемой задачи выполняется условие единственности решений [29].

3. УСЛОВИЕ ОПТИМАЛЬНОСТИ

Необходимое условие (23) минимума функционала (условие оптимальности) ${{J}_{\alpha }}$ можно записать в виде

где $A{\kern 1pt} *$ – сопряженный к $A$ оператор, который имеет вид: Запишем условие оптимальности в терминах исходных операторов. Для этого введем вспомогательную (“сопряженную”) задачу Тогда условие оптимальности запишется в виде системы операторных уравнений: Для того, чтобы выражения в (27) имели смысл, потребуем ${\mathbf{U}}{\kern 1pt} {{*}^{{,(i)}}} \in D(B{\kern 1pt} *)$. Таким образом, если функция $\omega = {{({{d}_{s}},\nu )}^{{\text{T}}}}$ доставляет минимум функционалу ${{J}_{\alpha }}$, $\alpha \geqslant 0$, и $(A\omega - \phi ) \in D(A{\kern 1pt} *)$, то она необходимо удовлетворяет условиям оптимальности (27), где ${\mathbf{U}}{\kern 1pt} {{*}^{{,(i)}}}$, $i = 1,2$ – решения систем (25), (26), в которых ${{{\text{U}}}^{{(i)}}}$ являются решениями систем (16), (17).

Домножим (27) скалярно в ${{\mathbb{H}}_{c}}$ на вектор-функцию $\hat {\omega } = {{({{\hat {d}}_{s}},\hat {\nu })}^{{\text{T}}}} \in {{H}_{c}}$. Пользуясь выражениями (14), (15) и представлениями билинейных форм (18), (19), получаем следующий интегральный аналог условий оптимальности (27):

Аналогичным образом можно получить интегральные соотношения и для “сопряженных” задач (25), (26).

Сопряженные задачи (25), (26) в дифференциальной форме могут быть записаны в виде

Условия оптимальности запишутся в виде Таким образом, система вариационных уравнений состоит из задач (29), (30), (28),(8), $i = 1,2$.

4. ПЛОТНАЯ РАЗРЕШИМОСТЬ

Рассмотрим вопрос о плотной разрешимости задачи.

Теорема 2. Задача (22) плотно разрешима.

Доказательство. Для того, чтобы задача (22) была плотно разрешима, достаточно показать, что ядро оператора $A{\kern 1pt} *$ состоит только из нулевого элемента (см. [31]). Предположим, что это не так, и существуют функция $\omega {\kern 1pt} * = {{(w,h)}^{{\text{T}}}} \in D(A{\kern 1pt} *)$, $w \ne 0$, $h \ne 0$, такая, что

Условие (31) можно записать в виде системы операторных уравнений: Задачу (31) можно сформулировать следующим образом: найти ${\text{U}}{\kern 1pt} {{*}^{{,(1)}}} \in {{\mathbb{X}}^{{(1)}}}$, ${\text{U}}{\kern 1pt} {{*}^{{,(2)}}} \in {{\mathbb{X}}^{{(2)}}}$, $w$, $h$, такие, что Данная задача (31) является обобщенной постановкой следующей задачи, записанной в дифференциальной форме:

Заметим, что функции $\xi {\kern 1pt} * = \{ \xi {\kern 1pt} {{*}^{{,(i)}}}\;{\text{в}}\;{{\Omega }_{i}},\;i = 1,2\} $, ${\mathbf{U}}{\kern 1pt} * = \{ {\mathbf{U}}{\kern 1pt} {{*}^{{,(i)}}}\;{\text{в}}\;{{\Omega }_{i}},\;i = 1,2\} $ являются решениями однородной сопряженной задачи в области $\Omega $ [21] и ${\text{U}}{\kern 1pt} * = 0$ в $\Omega $. Следовательно, ${\mathbf{U}}{\kern 1pt} {{*}^{{,(i)}}} = 0\;{\text{в}}\;{{\Omega }_{i}},\;i = 1,2$, и $w = 0$, $h = 0$. Таким образом, ${\text{ker}}(A{\kern 1pt} *) = \{ 0\} $, что и означает плотную разрешимость задачи (22).

5. ИТЕРАЦИОННЫЙ АЛГОРИТМ

Приближения к функциям ${{{\mathbf{U}}}^{{(i)}}}(\alpha )$, $\nu (\alpha )$, ${{d}_{s}}(\alpha )$, доставляющим минимум функционалу ${{J}_{\alpha }}$ при некотором малом $\alpha > 0$, будем искать с помощью следующего итерационного алгоритма.

Aлгоритм

Шаг 1. Пусть найдено ${{\omega }_{k}} = \mathop {({{{({{d}_{s}})}}_{k}},{{\nu }_{k}})}\nolimits^{\text{T}} $, где $k$ – номер итерации. Решаем задачи в каждой подобласти ${{\Omega }_{i}}$, $i = 1,2$:

Шаг 2. Вычисляем функции ${\mathbf{U}}{\kern 1pt} {{*}^{{,(1)}}}$, ${\mathbf{U}}{\kern 1pt} {{*}^{{,(2)}}}$ путем решения сопряженных задач в ${{\Omega }_{i}}$, $i = 1,2$:

Шаг 3. Ищем новое приближение ${{\omega }_{{k + 1}}} = \mathop {({{{({{d}_{s}})}}_{{k + 1}}},{{\nu }_{{k + 1}}})}\nolimits^{\text{T}} $:

Поскольку задача (22) плотно разрешима, итерационный параметр ${{\tau }_{k}}$ можно выбрать по следующей формуле [29]: где $J_{\alpha }^{'}(\omega )$ – производная Гато функционала ${{J}_{\alpha }}(\omega )$. При выполнении подходящего критерия окончания итерационного процесса элементы ${{({{d}_{s}})}_{k}}$, ${{\nu }_{k}}$, $\xi _{k}^{{(i)}}$, ${\text{U}}_{k}^{{(i)}}$ принимаются в качестве приближенных решений рассматриваемой задачи в ${{\Omega }_{i}}$, $i = 1,2$. Подходящий критерий остановки следует выбирать в зависимости от параметров итерационного алгоритма, размера моделируемой области, численного метода и заданной точности. Заметим, что решение задач в подобластях может осуществляться одновременно. В следующем разделе приведено численное исследование сходимости алгоритма (32)–(35).

6. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Для численного исследования предложенного итерационного алгоритма был проведен ряд численных экспериментов. Область была выбрана квадратной $\left( {0,L} \right) \times \left( {0,L} \right)$. Она была разделена на две подобласти ${{\Omega }_{1}}$ и ${{\Omega }_{2}}$ без пересечений. Пусть ${{\Gamma }_{{{\text{op}}}}} = 0 \times \left( {0,L} \right)$ – внешняя жидкая граница, $L = 100$ м. Значения параметров $g = 9.81$ м/с², $l = 0$, ${{R}_{f}} = 0$. Глубина выбиралась $H = 1$ м. Правые части ${\mathbf{\tilde {f}}}$ и ${{\tilde {f}}_{3}}$ равны $0$.

Как было отмечено ранее в разд. 1, для аппроксимации задачи (1) по времени была использована неявная схема. Шаг по времени был выбран 0.5 с. По пространству была выбрана равномерная сетка с шагами 1 м. Алгоритм совместного применения методов разделения области и вариационной ассимиляции (32)–(35) использовался на каждом шаге по времени.

Для задания начальных условий был проведен предварительный расчет (эксперимент). Шаги по времени и по пространству были выбраны, как описано ранее, однако область была представлена прямоугольником ${{\Omega }_{p}} = \left( { - L,L} \right) \times \left( {0,L} \right)$. Для предварительного эксперимента (здесь и далее для предварительного эксперимента будем использовать индекс $p$) были выбраны следующие начальные условия:

Граничные условия для предварительного эксперимента на $\partial {{\Omega }_{p}}$ были выбраны $H({\mathbf{U}} \cdot {\mathbf{n}}) = 0$. Правые части ${\mathbf{\tilde {f}}}$ и ${{\tilde {f}}_{3}}$ равны $0$. Начальные условия для решения задачи были выбраны как результат расчета в предварительном эксперименте после 10 с. Данные наблюдений ${{I}_{{{\text{obs}}}}}$ на внешней жидкой границе ${{\Gamma }_{{{\text{op}}}}}$ были получены из предварительного эксперимента (${{I}_{{{\text{obs}}}}} = H{{{\mathbf{U}}}_{p}} \cdot {{{\mathbf{n}}}_{1}}$), который был продолжен еще на 5 с.

Для сравнения результатов расчетов с применением алгоритма ассимиляции на внешней жидкой границе и без его использования с целью исследования влияния алгоритма был проведен эксперимент, в котором область также выбиралась квадратной ${{\Omega }_{s}} = (0,L) \times (0,L)$, но на границе $x = 0$ использовалось граничное условие непротекания $H({\mathbf{U}} \cdot {\mathbf{n}}) = 0$ (здесь и далее для этого эксперимента будем использовать индекс $s$).

При реализации алгоритма совместного применения методов разделения области и вариационной ассимиляции (32)–(35) критерий остановки был выбран следующим образом:

Таким образом, итерационный алгоритм заканчивается в случае, если часть, соответствующая регуляризованным слагаемым в функционале (24), начинает преобладать. Параметр регуляризации был выбран равным $\alpha = {{10}^{{ - 5}}}$. Зависимость величины относительно величины ${\text{Res}}$ на нулевой итерации от числа итераций при различных значениях параметра $\gamma $ представлена на фиг. 2. Как видно, при выборе $\gamma = 100$ значение ${\text{Res}}$ дает лучший результат (меньшее число итераций для достижения критерия выхода из итерационного процесса), причем за 27 итераций удается подавить значение ${\text{Res}}$ почти на 6 порядков, поэтому для последующих экспериментов будет использоваться $\gamma = 100$. Также из фиг. 2 видно, что сходимость алгоритма при выбранном параметре ${{\tau }_{k}}$ немонотонна, что согласуется с результатами численных исследований схожих алгоритмов решения некорректно поставленных задач [18], [19] и связано с ошибками решения прямых и сопряженных задач на сетках.

Результаты моделирования для функции уровня ${{\xi }^{{(i)}}}$, где $i = 1,2$, а также ${{\xi }_{p}}$ и ${{\xi }_{s}}$ представлены на фиг. 3. Расположение внутренней границы ${{\Gamma }_{{{\text{in}}}}}$ обозначено белой линией. На фиг. 3а показано распределение уровня, полученное по результатам расчета с применением методов ассимиляции и разделения области в случае, когда ${{\Gamma }_{{{\text{in}}}}}$ расположена на $x = 50$; на фиг. 3б – распределение уровня, полученное по результатам расчета с применением методов ассимиляции и разделения области в случае, когда ${{\Gamma }_{{{\text{in}}}}}$ расположена на $x = 5$. Функции уровня ${{\xi }^{{(i)}}}$, где $i = 1,2$, представленные на фиг. 3а, б, можно сравнить с ${{\xi }_{s}}$ и ${{\xi }_{p}}$, изображенными на фиг. 3в и г соответственно. Отметим, что качественно фиг. 3а, б и г схожи, а для 3в наблюдается некоторое отличие, расположенное вблизи внешней жидкой границы, связанное с тем, что для данного эксперимента в качестве граничного условия на внешней жидкой границе использовалось условие непротекания.

Для более детального изучения влияния предложенного алгоритма на результаты моделирования рассмотрим также модуль разности результатов для значений функций уровня (фиг. 4) и модуль разности для величины горизонтальной компоненты вектора скорости $u$ (фиг. 5). На фиг. 4а и б (фиг. 5а и б) представлены модули разности между решением с применением методов разделения области и вариационной ассимиляции данных и предварительным расчетом, где ${{\Gamma }_{{{\text{in}}}}}$ расположена на $x = 50$ и $x = 5$. Модули разности между решением без применения методов разделения области и вариационной ассимиляции данных и предварительным расчетом представлены на фиг. 4в и фиг. 5в (заметим, что для данных изображений использована другая шкала). Использование ассимиляции данных на внешней жидкой границе значительно уменьшает различия между результатами моделирования с использованием методов ассимиляции данных и разделения области и результатами предварительного эксперимента. Вблизи внутренней жидкой границы наблюдается ошибка, которая вносится при применении метода разделения области.

Для следующих численных экспериментов расположение ${{\Gamma }_{{{\text{in}}}}}$ было выбрано как $x = 50$ м, $0 < y < 100$ м. На фиг. 6 для сравнения представлены значения $H{\mathbf{U}} \cdot {{{\mathbf{n}}}_{1}}$ по результатам моделирования и данные наблюдений ${{I}_{{{\text{obs}}}}}$, которые задавались по результатам предварительного эксперимента. Использование алгоритма ассимиляции эффективно сближает результаты моделирования и данные наблюдений.

ЗАКЛЮЧЕНИЕ

В работе был рассмотрен алгоритм восстановления граничных функций на внешних и внутренних границах для линеаризованной системы уравнений мелкой воды. Была сформулирована обратная задача, проведено теоретическое исследование (доказана однозначная и плотная разрешимость задачи) и предложен итерационный алгоритм решения задачи. Для иллюстрации результатов теоретических исследований был рассмотрен ряд численных экспериментов для области простой формы. По проведенным численным экспериментам можно сделать следующие выводы: итерационный алгоритм совместного применения методов разделения области и вариационной ассимиляции данных на жидких границах сходится при предложенном параметре ${{\tau }_{k}} > 0$, ассимиляция данных эффективно сближает данные наблюдений и результаты моделирования. Кроме того, учет граничных условий на жидких границах значительно уменьшает расхождение результатов моделирования с референтным расчетом по сравнению с расхождением результатов моделирования без учета открытых границ. Отметим, что параметры $\alpha $, $\gamma $ необходимо выбирать в зависимости от рассматриваемой задачи. К практическим недостаткам алгоритма относится то, что подбор параметров осуществляется вручную. Эти параметры влияют на точность, сходимость и скорость сходимости итерационного алгоритма. В работе был рассмотрен численный эксперимент, в котором была продемонстрирована сходимость при разных параметрах $\gamma $. Результаты экспериментов показали, что алгоритм вносит некоторую ошибку в решение вблизи внутренней границы. Исследование применимости алгоритма в практических приложениях необходимо провести в дальнейшем, а также оценить вклад ошибки в решение.

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

Авторы выражают благодарность В.П. Шутяеву за обсуждение работы и полезные замечания.