Известия РАН. Физика атмосферы и океана, 2022, T. 58, № 5, стр. 493-503
Вариационный метод решения задачи о квазигеострофической циркуляции в двухслойном океане
В. Б. Залесный *
Институт вычислительной математики им. Г.И. Марчука РАН
119333 Москва, ул. Губкина, 8, Россия
* E-mail: vzalesny@yandex.ru
Поступила в редакцию 09.03.2022
После доработки 29.04.2022
Принята к публикации 09.06.2022
- EDN: XHNZYU
- DOI: 10.31857/S0002351522050133
Аннотация
Рассматривается новая вариационная постановка и метод решения задачи квазигеострофической динамики в двухслойном периодическом канале. Область имитирует лежащую в Южном океане зону Антарктического кругового течения. Особенностью задачи является двусвязность области ее решения (периодичность по широте). Используя разложение в ряды Фурье, задача сводится к нелинейной системе обыкновенных дифференциальных уравнений (ОДУ) по времени. Двусвязность области приводит к тому, что вместе с решением ОДУ требуется выполнить стационарное интегральное соотношение, определяющее полный расход течения. Предлагается вариационный численный алгоритм решения задачи, близкий к технике четырехмерного усвоения данных (4DVAR). Основой функции ценности является стационарное интегральное соотношение. С помощью серии вычислительных экспериментов изучаются стационарные режимы течений, зависящие от модельных параметров. Расчеты показывают, что наличие в рельефе дна высоких гармоник может вызывать формирование противотечения в нижнем слое. Противотечение устойчиво к небольшим вариациям возмущений рельефа и коэффициента турбулентной вязкости.
ВВЕДЕНИЕ
Квазигеострофические модели циркуляции в атмосфере и океане активно используются и развиваются на протяжении многих лет [1–13]. Можно выделить основные направления их развития. Это – математические вопросы, включая существование и единственность классических решений и задач ассимиляции данных [3–7]; методы исследования и оценка аттракторов атмосферной и океанической динамики [2, 3, 7, 9, 10], [6]; разработка численных алгоритмов [2, 6, 13]; исследование физических свойств их решений [8, 11, 12]. Во всех указанных направлениях получен ряд интересных результатов.
Данная работа является продолжением исследований [11, 12], развивая постановку задачи на нестационарный случай и двухслойный океан. Акцент делается не на параметризациях, а на новом вариационном методе решения задачи. Основная цель работы связана с построением эффективного численного алгоритма решения уравнений квазигеострофической динамики в двухслойном периодическом канале, имитирующем зону Антарктического кругового течения (АКТ). Математическая особенность задачи связана с тем, что область ее решения – двусвязная. Это приводит к тому, что вместе с эволюционной системой уравнений требуется выполнить стационарное интегральное соотношение. Рассматривается новый алгоритм, основанный на вариационном подходе, связанном с техникой сопряженных уравнений [2, 14–18]. С помощью серии вычислительных экспериментов изучаются некоторые стационарные режимы квазигеострофических течений, возникающие при различных модельных параметрах.
1. ИСХОДНЫЕ УРАВНЕНИЯ
Рассмотрим двухслойную квазигеострофическую модель течений в периодическом канале, параллельном оси $Ox$ [11, 12]. Система уравнений в верхнем 1-м и 2-м нижнем слое с граничными и начальными условиями имеет вид
(1.1)
$\begin{gathered} \frac{{\partial {{q}_{1}}}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {{{u}_{1}}{{q}_{1}} - {{k}_{1}}\frac{{\partial {{q}_{1}}}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{{v}_{1}}{{q}_{1}} - {{k}_{1}}\frac{{\partial {{q}_{1}}}}{{\partial y}}} \right) - \\ - \,\,{{\mu }_{1}}\Delta \left( {\frac{{\partial {{v}_{1}}}}{{\partial x}} - \frac{{\partial {{u}_{1}}}}{{\partial y}}} \right) = \frac{1}{{{{H}_{1}}}}\left( {\frac{{\partial {{\tau }_{y}}}}{{\partial x}} - \frac{{\partial {{\tau }_{x}}}}{{\partial y}}} \right), \\ \end{gathered} $(1.2)
$\begin{gathered} \frac{{\partial {{q}_{2}}}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {{{u}_{2}}{{q}_{2}} - {{k}_{2}}\frac{{\partial {{q}_{2}}}}{{\partial x}}} \right) + \\ + \,\,\frac{\partial }{{\partial y}}\left( {{{v}_{2}}{{q}_{2}} - {{k}_{2}}\frac{{\partial {{q}_{2}}}}{{\partial y}}} \right) + \left( {r - {{\mu }_{2}}\Delta } \right)\left( {\frac{{\partial {{v}_{2}}}}{{\partial x}} - \frac{{\partial {{u}_{2}}}}{{\partial y}}} \right) = 0, \\ \end{gathered} $(1.3)
${{\left. {{{q}_{i}}} \right|}_{{x = 0}}} = {{\left. {{{q}_{i}}} \right|}_{{x = {{L}_{x}}}}},\,\,\,\,{{\left. {\frac{{\partial {{q}_{i}}}}{{\partial x}}} \right|}_{{x = 0}}} = {{\left. {\frac{{\partial {{q}_{i}}}}{{\partial x}}} \right|}_{{x = {{L}_{x}}}}},\,\,\,\,i = 1,2,$(1.4)
${{v}_{i}}{{q}_{i}} - {{k}_{i}}\frac{{\partial {{q}_{i}}}}{{\partial y}} = 0\,\,{\text{при}}\,\,y = 0,L,$2. ПОСТРОЕНИЕ МОДЕЛЬНЫХ УРАВНЕНИЙ
Пусть компоненты трения ветра на поверхности океана заданы следующим образом
(2.1)
${{\tau }_{x}} = {{\tau }_{0}}\sin \left( {\frac{{\pi y}}{L}} \right),\,\,\,\,{{\tau }_{y}} = 0.$Cледуя [8, 11, 12], будем искать частное решение задачи (1.1)–(1.5) в виде
(2.3)
${{u}_{i}} \equiv - {{({{\psi }_{i}})}_{y}} = {{U}_{i}} - \frac{\pi }{L}\cos \left( {\frac{{\pi y}}{L}} \right){{\Phi }_{i}}(x),$(2.4)
${{v}_{i}} \equiv {{({{\psi }_{i}})}_{x}} = sin\left( {\frac{{\pi y}}{L}} \right)\Phi _{{ix}}^{{}}(x),$(2.5)
${{q}_{1}} = \Delta {{\psi }_{1}} + ({{f}_{0}} + \beta y) - \frac{{f_{0}^{2}}}{{g{\kern 1pt} '{{H}_{1}}}}\left( {{{\psi }_{1}} - {{\psi }_{2}}} \right),$(2.6)
${{q}_{2}} = \Delta {{\psi }_{2}} + ({{f}_{0}} + \beta y) + \frac{{f_{0}^{2}}}{{g{\kern 1pt} '{{H}_{2}}}}\left( {{{\psi }_{1}} - {{\psi }_{2}}} \right) + \frac{{{{f}_{0}}}}{{{{H}_{2}}}}B,$(2.7)
$B = sin\left( {\pi y{\text{/}}L} \right)h(x),\,\,\,\,g{\kern 1pt} ' = g\frac{{{{\rho }_{2}} - {{\rho }_{1}}}}{{{{\rho }_{0}}}}.$Подставим (2.1)–(2.4) в (1.1)–(1.2) и проинтегрируем уравнения с учетом граничных условий (1.4) по $y$ от 0 до $L$. После преобразований получим
(2.8)
$\begin{gathered} \frac{\partial }{{\partial t}}\left[ {{{\Phi }_{{1xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{1}} - \frac{\alpha }{{{{H}_{1}}}}\left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right)} \right] + \\ + \,\,{{\left[ {{{\Phi }_{{1xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{1}} - \frac{\alpha }{{{{H}_{1}}}}\left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right)} \right]}_{x}}{{U}_{1}} + \\ + \left[ {\beta + \frac{{\alpha \left( {{{U}_{1}} - {{U}_{2}}} \right)}}{{{{H}_{1}}}}} \right]{{\Phi }_{{1x}}} - \\ - \,\,{{k}_{1}}{{\left[ {{{\Phi }_{{1xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{1}} - \frac{\alpha }{{{{H}_{1}}}}\left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right)} \right]}_{{xx}}} - \\ - \,\,{{\mu }_{1}}\Delta \left( {{{\Phi }_{{1xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{1}}} \right) = 0, \\ \end{gathered} $(2.9)
$\begin{gathered} \frac{\partial }{{\partial t}}\left[ {{{\Phi }_{{2xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{2}} + \frac{\alpha }{{{{H}_{2}}}}\left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right)} \right] + \\ + \,\,{{\left[ {{{\Phi }_{{2xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{2}} + \frac{\alpha }{{{{H}_{2}}}}({{\Phi }_{1}} - {{\Phi }_{2}}) + \frac{{{{f}_{0}}h}}{{{{H}_{2}}}}} \right]}_{x}}{{U}_{2}} + \\ + \left[ {\beta - \frac{{\alpha \left( {{{U}_{1}} - {{U}_{2}}} \right)}}{{{{H}_{2}}}}} \right]{{\Phi }_{{2x}}} - \\ - \,\,{{k}_{2}}{{\left[ {{{\Phi }_{{2xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{2}} + \frac{\alpha }{{{{H}_{2}}}}({{\Phi }_{1}} - {{\Phi }_{2}}) + \frac{{{{f}_{0}}h}}{{{{H}_{2}}}}} \right]}_{{xx}}} + \\ + \,\,\left( {r - {{\mu }_{2}}\Delta } \right)\left( {{{\Phi }_{{2xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{2}}} \right) = 0. \\ \end{gathered} $Умножим (1.1), (1.2) на $\cos \left( {\pi y{\text{/}}L} \right)$ и проинтегрируем их вначале по $y$, а затем по$x$. Получим
(2.10)
$\begin{gathered} \frac{\gamma }{2}\frac{{\partial \left( {{{U}_{1}} - {{U}_{2}}} \right)}}{{\partial t}} - 2\alpha \left( {{{\Phi }_{{1x}}},{{\Phi }_{2}}} \right) + \\ + \,\,3{{L}_{x}}{{k}_{1}}\alpha \left( {{{U}_{1}} - {{U}_{2}}} \right) = {{f}_{1}}, \\ \end{gathered} $(2.11)
$\begin{gathered} \frac{\gamma }{2}\frac{{\partial \left( {{{U}_{1}} - {{U}_{2}}} \right)}}{{\partial t}} - 2\left( {\alpha {{\Phi }_{{1x}}} + {{f}_{0}}{{h}_{x}},{{\Phi }_{2}}} \right) + \\ + \,\,3{{L}_{x}}{{k}_{2}}\alpha \left( {{{U}_{1}} - {{U}_{2}}} \right) = {{f}_{2}}, \\ \end{gathered} $(2.12)
$\begin{gathered} {{f}_{1}} = \left( {\frac{{3{{\tau }_{0}}\pi }}{4} - 3{{k}_{1}}\beta {{H}_{1}}} \right){{L}_{x}},\,\,\,\,{{f}_{2}} = 3{{L}_{x}}{{k}_{2}}\beta {{H}_{2}}, \\ \gamma = \frac{{3{{L}_{x}}{{L}^{2}}\alpha }}{{{{\pi }^{2}}}},\,\,\,\,\alpha = f_{0}^{2}{\text{/}}g{\kern 1pt} ', \\ \end{gathered} $(2.13)
$\left( {{{\varphi }_{1}},{{\varphi }_{2}}} \right) \equiv \int\limits_0^{{{L}_{x}}} {{{\varphi }_{1}}{{\varphi }_{2}}dx} .$Замечание. Использование параметризации крупномасштабного турбулентного обмена со вторыми производными требует дополнительного граничного условия для ${{q}_{i}}$ на твердых стенках при $y = 0,\,\,L$. Естественными (в вариационном смысле) условиями будет обращение в ноль потоков ${{q}_{i}}$ по нормали вида (1.4). Учитывая кинематические условия ${{v}_{1}} = 0,\,\,{{v}_{2}} = 0$, выполнение (1.4) возможно, например, если положить ${{k}_{i}} = 0$. Вопрос постановки данного краевого условия обсуждается в [12].
3. ФОРМУЛИРОВКА УРАВНЕНИЙ В ТЕРМИНАХ ${{V}_{1}} = \left( {{{U}_{1}} + {{U}_{2}}} \right){\text{/}}2$, ${{V}_{2}} = \left( {{{U}_{1}} - {{U}_{2}}} \right){\text{/}}2$
Введем новые переменные
Получим
(3.1)
$\begin{gathered} \left( {\frac{\partial }{{\partial t}} + {{V}_{1}}\frac{\partial }{{\partial x}}} \right)\left[ {{{H}_{1}}\left( {{{\Phi }_{{1xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{1}}} \right) - \alpha \left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right)} \right] + \\ + \,\,{{\left[ {{{H}_{1}}\left( {{{\Phi }_{{1xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{1}}} \right) - \alpha \left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right)} \right]}_{x}}{{V}_{2}} + \\ + \,\,\left( {\beta {{H}_{1}} + 2\alpha {{V}_{2}}} \right){{\Phi }_{{1x}}} - \\ - \,\,{{k}_{1}}{{\left[ {{{H}_{1}}\left( {{{\Phi }_{{1xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{1}}} \right) - \alpha \left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right)} \right]}_{{xx}}} + \\ + \,\,{{H}_{1}}{{{\tilde {\mu }}}_{1}}\left( {{{\Phi }_{{1xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{1}}} \right) = 0, \\ \end{gathered} $(3.2)
$\begin{gathered} \left( {\frac{\partial }{{\partial t}} + {{V}_{1}}\frac{\partial }{{\partial x}}} \right) \times \\ \times \,\,\left[ {{{H}_{2}}\left( {{{\Phi }_{{2xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{2}}} \right) + \alpha \left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right) + {{f}_{0}}h} \right] - \\ - \,\,{{\left[ {{{H}_{2}}\left( {{{\Phi }_{{2xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{2}}} \right) + \alpha \left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right) + {{f}_{0}}h} \right]}_{x}}{{V}_{2}} + \\ + \,\,\left( {\beta {{H}_{2}} - 2\alpha {{V}_{2}}} \right){{\Phi }_{{2x}}} - \\ - \,\,{{k}_{2}}{{\left[ {{{H}_{2}}\left( {{{\Phi }_{{2xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{2}}} \right) + \alpha \left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right) + {{f}_{0}}h} \right]}_{{xx}}} + \\ + \,\,{{H}_{2}}(r + {{{\tilde {\mu }}}_{2}})\left( {{{\Phi }_{{2xx}}} - \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}{{\Phi }_{2}}} \right) = 0, \\ \end{gathered} $(3.3)
$\begin{gathered} - \frac{{3{{L}^{2}}\alpha }}{{{{\pi }^{2}}}}\frac{{\partial {{V}_{2}}}}{{\partial t}} + \frac{{2\alpha }}{{{{L}_{x}}}}\left( {{{\Phi }_{{1x}}},{{\Phi }_{2}}} \right) - 6{{k}_{1}}\alpha {{V}_{2}} = \\ = - 3\left( {\frac{{{{\tau }_{0}}\pi }}{4} - {{k}_{1}}\beta {{H}_{1}}} \right), \\ \end{gathered} $(3.4)
$\begin{gathered} \frac{{{{f}_{0}}}}{{{{L}_{x}}}}\left( {{{h}_{x}},{{\Phi }_{2}}} \right) - 3({{k}_{2}} - {{k}_{1}})\alpha {{V}_{2}} = \\ = \frac{3}{2}\left[ {\frac{{\pi {{\tau }_{0}}}}{4} - \beta ({{k}_{1}}{{H}_{1}} + {{k}_{2}}{{H}_{2}})} \right], \\ \end{gathered} $Интегральное соотношение (3.4) есть разность (2.10)–(2.11).
Закон сохранения. Умножим (3.1)–(3.2), соответственно на ${{\Phi }_{1}},\,\,{{\Phi }_{2}}$, проинтегрируем их по замкнутому контуру $x$ от 0 до ${{L}_{x}}$ и добавим соотношения (3.3), (3.4) умноженные, соответственно на ${{V}_{2}}$, $({{V}_{1}} - {{V}_{2}})$. Складывая все уравнения, получим энергетическое соотношение
Отсюда при ${{\mu }_{i}} = 0,\,\,\nu = 0\,$, ${{k}_{1}} = {{k}_{2}}\, = 0$ следует закон сохранения
(3.5)
$\begin{gathered} \frac{\partial }{{\partial t}}\int\limits_0^{{{L}_{x}}} {\left[ {{{H}_{1}}\left( {\Phi _{{1x}}^{2} + \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}\Phi _{1}^{2}} \right) + {{H}_{2}}\left( {\Phi _{{2x}}^{2} + \frac{{{{\pi }^{2}}}}{{{{L}^{2}}}}\Phi _{2}^{2}} \right)} \right.} + \\ \left. { + \,\,\alpha {{{\left( {{{\Phi }_{1}} - {{\Phi }_{2}}} \right)}}^{2}} + \frac{{3\alpha {{L}^{2}}}}{{{{\pi }^{2}}}}V_{2}^{2}} \right]dx = 0. \\ \end{gathered} $4. УРАВНЕНИЯ ДЛЯ КОЭФФИЦИЕНТОВ ФУРЬЕ. ВАРИАЦИОННАЯ ФОРМУЛИРОВКА НАЧАЛЬНО-КРАЕВОЙ ЗАДАЧИ
Будем искать решение (3.1)–(3.4) методом Галеркина, разлагая все функции от $x$ в ряды Фурье по полной ортогональной системе $\left\{ {sin\left( {2n\pi x{\text{/}}{{L}_{x}}} \right),} \right.$ $\left. {{\text{cos}}\left( {2n\pi x{\text{/}}{{L}_{x}}} \right)} \right\}$:
(4.1)
$\begin{gathered} {{\Phi }_{i}} = \sum\limits_n {a_{n}^{{(i)}}\left( t \right)\cos \left( {Nx} \right)} + \\ + \,\,\sum\limits_n {b_{n}^{{(i)}}\left( t \right)sin\left( {Nx} \right)} ,\,\,\,\,i = 1,2, \\ \end{gathered} $(4.2)
$\begin{gathered} h = \sum\limits_n {{{c}_{n}}\cos \left( {Nx} \right)} + \sum\limits_n {{{d}_{n}}sin\left( {Nx} \right)} , \\ N = 2n\pi {\text{/}}{{L}_{x}},n = 1,2,3, \ldots \\ \end{gathered} $Подставим (4.1)–(4.2) в (3.1)–(3.4) и умножим уравнения скалярно на базисные функции $sin\left( {Nx} \right)$, $\cos \left( {Nx} \right)$ для разных $n$. В силу ортогональности базиса получим после преобразований систему уравнений для коэффициентов разложений $a_{n}^{{\left( i \right)}},\,\,b_{n}^{{\left( i \right)}}$, ${{c}_{n}},\,\,{{d}_{n}}$
(4.3)
$\begin{gathered} \left( {\frac{d}{{dt}} + {{k}_{1}}{{N}^{2}}} \right)\left[ {\left( {{{H}_{1}}{{s}_{0}} + \alpha } \right)a_{n}^{{\left( 1 \right)}} - \alpha a_{n}^{{\left( 2 \right)}}} \right] + \\ + \,\,\left( {{{V}_{1}} + {{V}_{2}}} \right)N\left[ {\left( {{{H}_{1}}{{s}_{0}} + \alpha } \right)b_{n}^{{\left( 1 \right)}} - \alpha b_{n}^{{\left( 2 \right)}}} \right] = \\ = \,\,\left( {\beta {{H}_{1}} + 2\alpha {{V}_{2}}} \right)Nb_{n}^{{\left( 1 \right)}} - {{\mu }_{1}}{{H}_{1}}s_{0}^{2}a_{n}^{{\left( 1 \right)}}, \\ \end{gathered} $(4.4)
$\begin{gathered} \left( {\frac{d}{{dt}} + {{k}_{1}}{{N}^{2}}} \right)\left[ {\left( {{{H}_{1}}{{s}_{0}} + \alpha } \right)b_{n}^{{\left( 1 \right)}} - \alpha b_{n}^{{\left( 2 \right)}}} \right] - \\ - \,\,\left( {{{V}_{1}} + {{V}_{2}}} \right)N\left[ {\left( {{{H}_{1}}{{s}_{0}} + \alpha } \right)a_{n}^{{\left( 1 \right)}} - \alpha a_{n}^{{\left( 2 \right)}}} \right] = \\ = - \left( {\beta {{H}_{1}} + 2\alpha {{V}_{2}}} \right)Na_{n}^{{\left( 1 \right)}} - {{\mu }_{1}}{{H}_{1}}s_{0}^{2}b_{n}^{{\left( 1 \right)}}, \\ \end{gathered} $(4.5)
$\begin{gathered} \left( {\frac{d}{{dt}} + {{k}_{2}}{{N}^{2}}} \right)\left[ {\left( {{{H}_{2}}{{s}_{0}} + \alpha } \right)a_{n}^{{\left( 2 \right)}} - \alpha a_{n}^{{\left( 1 \right)}} - {{f}_{0}}{{c}_{n}}} \right] + \\ + \,\,\left( {{{V}_{1}} - {{V}_{2}}} \right)N\left[ {\left( {{{H}_{2}}{{s}_{0}} + \alpha } \right)b_{n}^{{\left( 2 \right)}} - \alpha b_{n}^{{\left( 1 \right)}} - {{f}_{0}}{{d}_{n}}} \right] = \\ = \left( {\beta {{H}_{2}} - 2\alpha {{V}_{2}}} \right)Nb_{n}^{{\left( 2 \right)}} - {{H}_{2}}{{s}_{0}}\left( {r + {{\mu }_{2}}{{s}_{0}}} \right)a_{n}^{{\left( 2 \right)}}, \\ \end{gathered} $(4.6)
$\begin{gathered} \left( {\frac{d}{{dt}} + {{k}_{2}}{{N}^{2}}} \right)\left[ {\left( {{{H}_{2}}{{s}_{0}} + \alpha } \right)b_{n}^{{\left( 2 \right)}} - \alpha b_{n}^{{\left( 1 \right)}} - {{f}_{0}}{{d}_{n}}} \right] - \\ - \,\,\left( {{{V}_{1}} - {{V}_{2}}} \right)N\left[ {\left( {{{H}_{2}}{{s}_{0}} + \alpha } \right)a_{n}^{{\left( 2 \right)}} - \alpha a_{n}^{{\left( 1 \right)}} - {{f}_{0}}{{c}_{n}}} \right] = \\ = - \left( {\beta {{H}_{2}} - 2\alpha {{V}_{2}}} \right)Na_{n}^{{\left( 2 \right)}} - {{H}_{1}}{{s}_{0}}\left( {r + {{\mu }_{2}}{{s}_{0}}} \right)b_{n}^{{\left( 2 \right)}}. \\ \end{gathered} $(4.7)
$\begin{gathered} \frac{{3\alpha {{L}^{2}}}}{\pi }\frac{{\partial {{V}_{2}}}}{{\partial t}} + 6\alpha {{k}_{1}}{{V}_{2}} = 3\left( {\frac{{\pi {{\tau }_{0}}}}{4} - \beta {{k}_{1}}{{H}_{1}}} \right) - \\ - \,\,\alpha \sum\limits_n {N\left( {a_{n}^{{(1)}}b_{n}^{{(2)}} - b_{n}^{{(1)}}a_{n}^{{(2)}}} \right)} , \\ \end{gathered} $(4.8)
$\begin{gathered} 6\alpha ({{k}_{1}} - {{k}_{2}}){{V}_{2}} = 3\left[ {\frac{{\pi {{\tau }_{0}}}}{4} - \beta \left( {{{k}_{1}}{{H}_{1}} + {{k}_{2}}{{H}_{2}}} \right)} \right] + \\ + \,\,\sum\limits_n {N{{f}_{0}}\left( {{{c}_{n}}b_{n}^{{(2)}} - {{d}_{n}}a_{n}^{{(2)}}} \right)} , \\ \end{gathered} $Таким образом, в результате преобразований исходная задача преобразована к системе шести уравнений (4.3)–(4.8) для шести неизвестных $a_{n}^{{\left( 1 \right)}}$, $b_{n}^{{\left( 1 \right)}}$, $a_{n}^{{\left( 2 \right)}}$, $b_{n}^{{\left( 2 \right)}}$, ${{U}_{1}}$, ${{U}_{2}}$. Видно, что если ${{U}_{i}}$ заданы, то уравнения (4.3)–(4.6) распадаются на несвязанные подсистемы для отдельной гармоники $n$. При конечном наборе отличных от нуля коэффициентов ${{c}_{n}} \ne 0,\,\,{{d}_{n}} \ne 0$ ряды в (4.1)–(4.2) и система уравнений (4.3)–(4.8) конечны.
Итак, задача сведена к решению эволюционных уравнений (4.3)–(4.7) и стационарного интегрального соотношения (4.8). Решением задачи является вектор-функция Ψ(t) = $ = \left( {a_{n}^{{\left( 1 \right)}}\left( t \right),\,\,b_{n}^{{\left( 1 \right)}}\left( t \right),\,\,a_{n}^{{\left( 2 \right)}}\left( t \right),\,\,b_{n}^{{\left( 2 \right)}}\left( t \right),\,\,{{V}_{2}}\left( t \right)} \right)$, и параметр ${{V}_{1}}$. Особенностью задачи, значительно осложняющей ее численное решение, является отсутствие уравнения для ${{V}_{1}}$. ${{V}_{1}}$ является параметром системы (4.3)–(4.7) и его следует подобрать так, чтобы в каждый момент времени выполнялось соотношение (4.8). Переформулируем задачу в виде, используемом при решении задач 4-х мерной вариационной ассимиляции [2, 13–16].
Следуя [2, 15], будем искать такое решение (4.3)–(4.7) с неизвестным параметром ${{V}_{1}}$, на котором достигается минимум функционала J
(4.9)
$J = \frac{\lambda }{2}{{{\rm Z}}^{2}}\left( {{{t}_{0}}} \right) + \frac{1}{{2T}}\int\limits_{{{t}_{0}}}^{{{t}_{0}} + T} \eta (t){{{\rm Z}}^{2}}\left( t \right)dt \to \min ,$то функционал (4.9) можно выбрать в виде
(4.10)
$\begin{gathered} J = \frac{\lambda }{2}{{{\rm Z}}^{2}}\left( {{{t}_{0}}} \right) + \frac{1}{{2T}}\int\limits_{{{t}_{0}}}^{{{t}_{0}} + T} \eta (t){{{\rm Z}}^{2}}\left( t \right)dtdt + \\ + \,\,\frac{\varepsilon }{2}{{\left[ {{{V}_{1}} + \frac{{{{H}_{1}} - {{H}_{2}}}}{{{{H}_{1}} + {{H}_{2}}}}{{V}_{2}}\left( {{{t}_{0}}} \right) - U_{{brt}}^{{dat}}\left( {{{t}_{0}}} \right)} \right]}^{2}}. \\ \end{gathered} $5. СОПРЯЖЕННЫЕ УРАВНЕНИЯ ДЛЯ КОЭФФИЦИЕНТОВ ФУРЬЕ. МЕТОД РЕШЕНИЯ СИСТЕМЫ ОПТИМАЛЬНОСТИ
Используя известную технику [2, 15–18], решение вариационной задачи (4.3)–(4.7), (4.9) можно свести к решению системы прямых и сопряженных уравнений, часто называемой системой оптимальности.
Перепишем (4.3)–(4.7), (4.9) в матричном виде
Прямые уравнения
(5.1)
$\begin{gathered} \frac{{3\alpha {{L}^{2}}}}{\pi }\frac{{\partial {{V}_{2}}}}{{\partial t}} + 6\alpha {{k}_{1}}{{V}_{2}} + \\ + \,\,\alpha \sum\limits_n {N\left( {a_{n}^{{(1)}}b_{n}^{{(2)}} - b_{n}^{{(1)}}a_{n}^{{(2)}}} \right) = 3\left( {\frac{{\pi {{\tau }_{0}}}}{4} - \beta {{k}_{1}}{{H}_{1}}} \right).} \\ \end{gathered} $Сопряженные уравнения
(5.2)
$\begin{gathered} \left[ {\left( { - \frac{d}{{dt}} + {{k}_{2}}{{N}^{2}}} \right){{s}_{2}} + {{H}_{1}}{{s}_{0}}\left( {r + {{\mu }_{2}}{{s}_{0}}} \right)} \right]b{\kern 1pt} *_{n}^{{\left( 2 \right)}} - \\ - \,\,\alpha \left( {{{V}_{1}} + {{V}_{2}}} \right)Na{\kern 1pt} *_{n}^{{\left( 1 \right)}} - \,\,\alpha \left( { - \frac{d}{{dt}} + {{k}_{1}}{{N}^{2}}} \right)b{\kern 1pt} *_{n}^{{\left( 1 \right)}} + \\ + \,\,N\left[ {\left( {{{V}_{1}} - {{V}_{2}}} \right){{s}_{2}} - \left( {\beta {{H}_{2}} - 2\alpha {{V}_{2}}} \right)} \right]a{\kern 1pt} *_{n}^{{\left( 2 \right)}} + \\ + \,\,\alpha \sum\limits_n {Na_{n}^{{(1)}}V_{2}^{*} = - {{f}_{0}}\eta \left( {\sum\limits_n {N{{c}_{n}}} } \right)\bar {J},} \\ \end{gathered} $(5.3)
$\begin{gathered} \bar {J} = \frac{1}{T}\int\limits_t^{t + T} \eta \left\{ {{{f}_{0}}\sum\limits_n N \left( {a_{n}^{{\left( 2 \right)}}{{d}_{n}} - b_{n}^{{\left( 2 \right)}}{{c}_{n}}} \right)} \right. - \\ \left. { - \,\,6({{k}_{2}} - {{k}_{1}})\alpha {{V}_{2}} - 3\left[ {\frac{{\pi {{\tau }_{0}}}}{4} - \beta ({{k}_{1}}{{H}_{1}} + {{k}_{2}}{{H}_{2}})} \right]} \right\}dt. \\ \end{gathered} $Напомним, что сопряженные уравнения можно получить, умножая прямые уравнения, соответственно, на $a{\kern 1pt} *_{n}^{{\left( 1 \right)}},\,\,b{\kern 1pt} *_{n}^{{\left( 1 \right)}},\,\,\,a{\kern 1pt} *_{n}^{{\left( 2 \right)}},\,\,b{\kern 1pt} *_{n}^{{\left( 2 \right)}},\,\,V_{2}^{*}$, интегрируя их по времени на интервале усвоения от ${{t}_{0}}$ до ${{t}_{0}} + T$ и приравнивая нулю производные по $a_{n}^{{\left( 1 \right)}},\,\,b_{n}^{{\left( 1 \right)}},\,\,\,a_{n}^{{\left( 2 \right)}},\,\,b_{n}^{{\left( 2 \right)}},\,\,{{V}_{2}}$. Учет функционала (4.10) приводит к появлению слагаемых в правой части (5.2). В качестве управления выбираются начальные условия для прямой системы $a_{n}^{{\left( 1 \right)}},\,\,b_{n}^{{\left( 1 \right)}},\,\,\,a_{n}^{{\left( 2 \right)}},\,\,b_{n}^{{\left( 2 \right)}},\,\,{{V}_{2}}$ (в момент времени $t = {{t}_{0}}$) и параметр ${{V}_{1}}$.
Система оптимальности (5.1)–(5.2) решается итерационно на отрезке $({{t}_{0}},\,\,{{t}_{0}} + T)$, причем прямые уравнения (5.1) решаются вперед по времени, а сопряженные (5.2) в обратном времени. Решением задачи являются $a_{n}^{{\left( 1 \right)}}\left( t \right),\,\,b_{n}^{{\left( 1 \right)}}\left( t \right),\,\,\,a_{n}^{{\left( 2 \right)}}\left( t \right),\,\,b_{n}^{{\left( 2 \right)}}\left( t \right),\,\,{{V}_{2}}\left( t \right)$, $a{\kern 1pt} *_{n}^{{\left( 1 \right)}}\left( t \right),b{\kern 1pt} *_{n}^{{\left( 1 \right)}}\left( t \right),a{\kern 1pt} *_{n}^{{\left( 2 \right)}}\left( t \right),b{\kern 1pt} *_{n}^{{\left( 2 \right)}}\left( t \right),V_{2}^{*}\left( t \right)$, плюс вектор управления $a_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}}} \right),b_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}}} \right),a_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}}} \right),b_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}}} \right),{{V}_{2}}\left( {{{t}_{0}}} \right)$, V1. Управление находится с помощью итерационного процесса: в наших расчетах использован алгоритм M1QN3 [19]. При решении прямой системы (5.1) используются начальные условия при $\,t = {{t}_{0}}$
(5.4)
$\begin{gathered} a_{n}^{{\left( 1 \right)}}\left( t \right) = a_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}}} \right),\,\,b_{n}^{{\left( 1 \right)}}\left( t \right) = b_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}}} \right), \\ a_{n}^{{\left( 2 \right)}}\left( t \right) = a_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}}} \right),\,\,\,b_{n}^{{\left( 2 \right)}}\left( t \right) = b_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}}} \right),\,\,\,{{V}_{2}}\left( t \right) = {{V}_{2}}\left( {{{t}_{0}}} \right). \\ \end{gathered} $При решении сопряженной системы (5.2) используются нулевые условия при $\,t = {{t}_{0}} + T$
(5.5)
$\begin{gathered} a{\kern 1pt} *_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}} + T} \right) = b{\kern 1pt} *_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}} + T} \right) = a{\kern 1pt} *_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}} + T} \right) = \\ = \,\,b{\kern 1pt} *_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}} + T} \right) = V_{2}^{*}\left( {{{t}_{0}} + T} \right) = 0. \\ \end{gathered} $Правые части в начальных условиях (5.4) (компоненты вектора $\,{{\Psi }_{0}}$) зависят от компонент градиента лагранжиана, которые в свою очередь связаны с решением сопряженной системы при $\,t = {{t}_{0}}$. Нумеруя компоненты производной от лагранжиана ${{(\partial L)}_{i}}$ в соответствии с компонентами ${{\Psi }_{0}} = \left( {a_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}}} \right)\,,\,\,\,b_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}}} \right),\,\,\,a_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}}} \right)\,,\,\,\,b_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}}} \right),\,\,{{V}_{2}}\left( {{{t}_{0}}} \right),\,\,{{V}_{1}}} \right)$, получим
(5.6)
$\begin{gathered} {{(\partial L)}_{5}} = V{\kern 1pt} *_{n}^{{\left( 2 \right)}}(t) + \\ + \,\,\varepsilon \left( {{{V}_{1}} + \frac{{{{H}_{1}} - {{H}_{2}}}}{{{{H}_{1}} + {{H}_{2}}}}{{V}_{2}} - {{V}_{{dat}}}} \right)\frac{{{{H}_{1}} - {{H}_{2}}}}{{{{H}_{1}} + {{H}_{2}}}}, \\ \end{gathered} $6. ЧИСЛЕННЫЕ РАСЧЕТЫ
Расчеты сделаны для периодического зонального канала, протяженностью Lx = 4 × 106 м, шириной L = 106 м, глубиной слоев H1 = 103, H2 = $ = 4 \times {{10}^{3}}$ м (общая глубина ${{H}_{1}} + {{H}_{2}} = 5 \times {{10}^{3}}$ м). Расчетная область имитирует положение Антарктического кругового течения, которое находится в Южном полушарии. В расчетах использовались следующие (либо близкие к ним) параметры
(6.1)
$\begin{gathered} {{\tau }_{0}} = {{10}^{{ - 4}}}{{m}^{2}}{{\sec }^{{ - 1}}},\,\,\,\,{{f}_{0}} = - {{10}^{{ - 4}}}{{\sec }^{{ - 1}}}, \\ \beta = 1.4 \times {{10}^{{ - 11}}}{{m}^{{ - 1}}}{{\sec }^{{ - 1}}},\,\,\,\,{{L}_{x}} = 4 \times {{10}^{6}}m, \\ {{L}_{y}} \equiv L = {{10}^{6}}m,\,\,\,\,{{Н}_{1}} = {{10}^{3}}m,\,\,\,{{H}_{2}} = 4 \times {{10}^{3}}\,m, \\ {{\mu }_{1}} = {{\mu }_{2}} = 1{{m}^{2}}{{\sec }^{{ - 1}}},\,\,\,\,{{k}_{1}} = {{10}^{3}},\,\,\,{{k}_{2}} = {{10}^{3}}{{m}^{2}}{{\sec }^{{ - 1}}}, \\ r = {{10}^{{ - 7}}}{{\sec }^{{ - 1}}},\,\,\,\,\alpha = \frac{{f_{0}^{2}\rho }}{{g\Delta \rho }} \approx {{10}^{{ - 6}}}{{m}^{{ - 1}}}, \\ {{\rho }_{1}} = {{10}^{3}}(1. + 0.026),\,\,{{\rho }_{2}} = {{10}^{3}}(1. + 0.028)\,\,kg{{m}^{{ - 1}}}. \\ \end{gathered} $Дополнительно укажем значения некоторых величин и параметры рельефа дна:
(6.2)
$\begin{gathered} N = 2\pi n{\text{/}}{{L}_{x}},\,\,\,\,{{s}_{0}} = {{N}^{2}} + {{\pi }^{2}}{\text{/}}{{L}^{2}}, \\ {{s}_{i}} = {{H}_{i}}{{s}_{0}} + \alpha ,\,\,\,\,B = sin\left( {\pi y{\text{/}}L} \right)h(x), \\ h(x) = \sum\limits_n {\left[ {{{c}_{n}}cos(2\pi nx{\text{/}}{{L}_{x}}) + {{d}_{n}}sin(2\pi nx{\text{/}}{{L}_{x}})} \right]} . \\ \end{gathered} $Алгоритм. Процесс решения задачи состоял из двух этапов. На первом этапе решалась прямая система (5.1) с фиксированным ${{V}_{1}}$ по схеме Кранка-Николсон. Интервал расчета по времени $0 \leqslant t \leqslant {{t}_{0}} + T$ был достаточно большим. На втором этапе, на небольшом интервале по времени ${{t}_{0}} \leqslant t \leqslant {{t}_{0}} + T$, находилось решение системы оптимальности (5.1)–(5.2), (5.4)–(5.5), (4.9), включая параметр ${{V}_{1}}$. Данная двухэтапная процедура итерационно повторялась до достижения (с заданной точностью) сходимости решения полной задачи. На втором этапе для решения оптимальной задачи использовалась стандартная программа M1QN3 [19]. Полное решение задачи, таким образом, включало решение прямой системы $a_{n}^{{\left( 1 \right)}}\left( t \right),\,\,b_{n}^{{\left( 1 \right)}}\left( t \right),\,\,\,a_{n}^{{\left( 2 \right)}}\left( t \right),\,\,b_{n}^{{\left( 2 \right)}}\left( t \right),\,\,\,{{V}_{2}}\left( t \right)$ при ${{t}_{0}} < t \leqslant {{t}_{0}} + T$, решение сопряженной системы $a{\kern 1pt} *_{n}^{{\left( 1 \right)}}\left( t \right),b{\kern 1pt} *_{n}^{{\left( 1 \right)}}\left( t \right),a{\kern 1pt} *_{n}^{{\left( 2 \right)}}\left( t \right),b{\kern 1pt} *_{n}^{{\left( 2 \right)}}\left( t \right),V_{2}^{*}\left( t \right)$ при t0 ≤ t ≤ $ \leqslant {{t}_{0}} + T$ и вектор управления $a_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}}} \right),b_{n}^{{\left( 1 \right)}}\left( {{{t}_{0}}} \right),$ $a_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}}} \right),b_{n}^{{\left( 2 \right)}}\left( {{{t}_{0}}} \right),{{V}_{2}}\left( {{{t}_{0}}} \right),{{V}_{1}}$.
Во всех расчетах в (4.9) было положено $\lambda = 1$, шаг по времени составлял $\Delta t = {{10}^{6}}\,с \approx \,10$ сут., интервал 1-го этапа ${{10}^{{13}}}\,\,c \approx 2.7 \times {{10}^{6}}$ лет, интервал 2-го этапа =30 сут., число полных циклов 2-х этапной итерационной процедуры равнялось 20. В серии экспериментов использовались различные значения коэффициентов трения и диффузии ${{k}_{i}},\,\,r,\,\,{{\mu }_{i}}$, глубин слоев и т.д.
Замечание. Оценим приблизительное время выхода решения на стационарный режим, которое определяют параметры ${{k}_{i}},\,\,r,\,\,{{\mu }_{i}}$. Удержим в (5.1) слагаемые, обеспечивающие экспоненциальное затухание возмущений решения по времени. Оставляя только «диссипативную» часть четвертого уравнения (5.1), видно, что затухание происходит за счет трех диссипативных процессов, соответственно, с коэффициентами ${{k}_{2}},\,\,r,\,\,{{\mu }_{2}}$
(6.3)
$\begin{gathered} \frac{d}{{dt}}b_{n}^{{\left( 2 \right)}} + .... = - \left[ {{{k}_{2}}{{N}^{2}} + {{H}_{1}}{{s}_{0}}\left( {r + {{\mu }_{2}}{{s}_{0}}} \right){\text{/}}{{s}_{2}}} \right]b_{n}^{{\left( 2 \right)}}, \\ b_{n}^{{\left( 2 \right)}}\left( t \right) = b_{n}^{{\left( 2 \right)}}({{t}_{0}}) \times \\ \times \,\,\exp \left\{ { - \left[ {{{k}_{2}}{{N}^{2}} + {{H}_{1}}{{s}_{0}}\left( {r + {{\mu }_{2}}{{s}_{0}}} \right){\text{/}}{{s}_{2}}} \right]t} \right\}. \\ \end{gathered} $Это позволяет грубо оценить время стационирования (уменьшения в ${{e}^{\kappa }}$) $n$-ой гармоники решения для всех диссипативных процессов – диффузии с коэффициентом ${{\mu }_{2}}$. С учетом (6.1)–(6.2) имеем
(6.4)
$\begin{gathered} t \leqslant \kappa {{\left[ {{{k}_{2}}{{N}^{2}} + {{H}_{1}}{{s}_{0}}\left( {r + {{\mu }_{2}}{{s}_{0}}} \right){\text{/}}{{s}_{2}}} \right]}^{{ - 1}}} \approx \\ \approx \,\,\kappa {{\left( {2.5 \times {{{10}}^{{ - 9}}} + {{{10}}^{{ - 10}}} + {{{10}}^{{ - 15}}}{{n}^{2}}} \right)}^{{ - 1}}}{{n}^{{ - 2}}}. \\ \end{gathered} $Из (6.4) видно, что самый медленный процесс – это диффузия с коэффициентом ${{\mu }_{2}}$. Время затухания возмущений за каждого диссипативного процесса (с коэффициентами ${{k}_{2}},\,\,r,\,\,{{\mu }_{2}}$) в е-раз будет при $n = 1$, соответственно, $250$, $3000$ и $3 \times {{10}^{8}}$ лет.
Эксперимент 1, базовый. В рамках базового эксперимента было проведено 7 расчетов: 1.1–1.7, отличающихся друг от друга высотой рельефа, глубиной канала и значением ${{\mu }_{i}}$. Основная цель эксперимента состояла в подборе параметров и оценке особенностей численного решения оптимальной задачи. Во всех расчетах кроме 1.2 было положено $n = 1$, в первых шести расчетах принято ${{H}_{1}} = 1$ км, ${{H}_{2}} = 4$ км, ${{H}_{1}} + {{H}_{2}} = 5$ км, в последнем глубины уменьшены в 5 раз. Основные параметры модели базового эксперимента приведены в таблице 1. Здесь ${{V}_{{ACC}}}$ – расход в свердрупах (аналог расхода АКТ): ${{V}_{{ACC}}} = {{U}_{1}}{{H}_{1}} + {{U}_{2}}{{H}_{2}},(1\,\,sv \equiv {{10}^{6}}\,{{{\text{m}}}^{{\text{3}}}}{\text{/s)}}$; знак ∞ означает стационарность решения, ${{J}_{0}}$, ${{J}_{{opt}}}$ – умноженные на ${{10}^{{10}}}$, начальное и конечное значение минимизируемого функционала или функции ценности (4.9). Различные значения ${{J}_{0}}$ связаны с тем, что почти во всех начальных приближениях было положено ${{V}_{1}} = 0.1,\,\,{{V}_{2}} = {{10}^{{ - 4}}}$.
Таблица 1.
№ | ${{c}_{n}}$ | ${{d}_{n}}$ | ${{k}_{1}} = {{k}_{2}}$ | $r$ | ${{\mu }_{1}} = {{\mu }_{2}}$ | $n$ | ${{U}_{1}},\,\,{{U}_{2}}$${\text{cm/s}}$ | ${{V}_{{ACC}}}$, $sv$ | ${{a}_{1}}$, ${{b}_{1}}$, | ${{a}_{2}}$, ${{b}_{2}}$ | ${{J}_{0}},{{J}_{{opt}}}$$ \times {{10}^{{10}}}$ |
---|---|---|---|---|---|---|---|---|---|---|---|
1.1 | 100 | 0 | 103 | 10–7 | 1 | 1 | 23 17 |
919 ∞ |
39 114 2257 |
26 760 1631 |
1 10–7 |
1.2 | 100 | 0 | 103 | 10–7 | 1 | 2 | 12 6.5 |
383 ∞ |
19 306 1220 |
8320 815 |
$4 \times {{10}^{{ - 2}}}$ ${{10}^{{ - 9}}}$ |
1.3 | 200 | 0 | 103 | 10–7 | 1 | 1 | 10 5 |
290 ∞ |
30 975 1559 |
9814 815 |
0.3 $3 \times {{10}^{{ - 9}}}$ |
1.4 | 200 | 0 | 103 | 10–7 | 103 | 1 | 10 4 |
264 ∞ |
29 188 1602 |
8226 816 |
0.9 $5 \times {{10}^{{ - 8}}}$ |
1.5 | 100 | 0 | 103 | 10–5 | 1 | 1 | 12 0.05 |
339 ∞ |
1257 4553 |
408 1631 |
0.3 ${{10}^{{ - 9}}}$ |
1.6 | 100 | 0 | 103 | $5 \times {{10}^{{ - 7}}}$ | 1 | 1 | 13 0.06 |
383 ∞ |
18 174 3683 |
7339 1629 |
0.2 $5 \times {{10}^{{ - 7}}}$ |
1.7 | 100 ${{Н}_{i}}/5$ |
0 | 103 | 10–7 | 1 | 1 | 27 22 |
232 ∞ |
235 390 14606 |
193 813 12284 |
92 $8 \times {{10}^{{ - 4}}}$ |
Расчеты показывают, что использованные параметры модели дают разумную величину зональной скорости в верхнем слое $ \sim 10 - 20$ см/с. Первый вариант 1.1 приводит к завышенной скорости в нижнем слое, что связано, как показывают расчеты 1.2, 1.3, с номером гармоники $n = 1$ и небольшой амплитудой рельефа. Во всех случаях вариационный алгоритм показал высокую скорость сходимости и точность. Минимизируемый функционал (4.9) уменьшался в процессе счета для 6-и расчетов, примерно, на 7–8 порядков, а для 1.7 – на 5 порядков.
Оценим использованные значения коэффициента придонного трения $r$, сравнивая его с величиной ${{k}_{2}}$. Если предположить, что $r$ лежит в пределах ${{\left( {2\pi n{\text{/}}{{L}_{x}}} \right)}^{2}}{{k}_{2}} < \varepsilon < {{10}^{3}}{{\left( {2\pi n{\text{/}}{{L}_{x}}} \right)}^{2}}{{k}_{2}}$, то при ${{k}_{2}} = {{10}^{3}} - {{10}^{4}},\,\,n = 1$, можно получить, что 10–9 < $ < r < {{10}^{{ - 6}}}$. Это явилось основанием выбора его значений в вариантах 1.5, 1.6. Сравнение вариантов 1.1, 1.5, 1.6, показывает разумность величины $r \approx 5 \times {{10}^{{ - 7}}}$, при которой полный расход течений ближе к наблюдаемому расходу АКТ. Отметим, что если глубины уменьшить в три раза: ${{Н}_{i}}{\text{/}}3$, то получим ${{U}_{1}} = 34$ ${\text{см/с}}$, ${{U}_{2}} = 28$ ${\text{см/с}}$, а расход будет ${{V}_{{ACC}}} = 492$ св. Для этого случая приблизительно выполняется реальная пропорция отношения длины и глубины АКТ.
Эксперименты 2, 3. Зависимость решения от номера гармоники n. Цель второго и третьего экспериментов состояла в изучении зависимости решения от гармоники n. Их результаты приведены в таблицах 2, 3. Расчеты проводились для различных параметров придонного трения и турбулентной вязкости КГВ. Выбирая в качестве основного ${{k}_{1}} = {{k}_{2}} \equiv k$, выпишем связь и зависимость остальных от n. Имеем
Таблица 2.
№ | ${{c}_{n}}$ | ${{d}_{n}}$ | ${{k}_{i}}$ | $r$ | ${{\mu }_{i}}$ | $n$ | ${{U}_{1}},\,\,{{U}_{2}}$${\text{cm/s}}$ | ${{V}_{{ACC}}}$, $sv$ | ${{a}_{1}}$${{b}_{1}}$ | ${{a}_{2}}$ ${{b}_{2}}$ |
${{J}_{0}},{{J}_{{opt}}}$$ \times {{10}^{{10}}}$ |
---|---|---|---|---|---|---|---|---|---|---|---|
1.1 | 100 | 0 | 103 | 10–7 | $1$ | $1$ | 23 17 |
919 ∞ |
39 114 2257 |
26 760 1631 |
1 10–7 |
2.1 | 100 | 0 | 103 | 10–7 | 1 | 2 | 12 6.5 |
383 ∞ |
19 306 1220 |
8320 815 |
$4 \times {{10}^{{ - 2}}}$ 10–9 |
2.2 | 100 | 0 | 103 | 10–6 | 1 | 3 | 65 0.7 |
93 ~ |
217 10–3 |
||
2.3 | 100 | 0 | 103 | 10–6 | 1 | 4 | 6 0.2 |
69 ∞ |
4795 1363 |
–813 408 |
1022 $4 \times {{10}^{{ - 6}}}$ |
2.4 | 100 | 0 | 103 | 10–6 | 1 | 5 | 55 –0.05 |
52 ∞ |
3895 2174 |
–1046 326 |
3921 $4 \times {{10}^{{ - 8}}}$ |
2.5 | 100 | 0 | 103 | $3 \times {{10}^{{ - 6}}}$ | 1 | 6 | 6 0.3 |
73 ∞ |
2231 2129 |
–499 271 |
1280 $3 \times {{10}^{{ - 9}}}$ |
Таблица 3.
№ | ${{c}_{n}}$ | ${{d}_{n}}$ | ${{k}_{i}}$ | $r$ | ${{\mu }_{i}}$ | $n$ | ${{U}_{1}},\,\,{{U}_{2}}$${\text{cm/s}}$ | ${{V}_{{ACC}}}$, $sv$ | ${{a}_{1}}$ ${{b}_{1}}$ |
${{a}_{2}}$ ${{b}_{2}}$ |
${{J}_{0}},{{J}_{{opt}}}$$ \times {{10}^{{10}}}$ |
---|---|---|---|---|---|---|---|---|---|---|---|
3.1 | 50 | 0 | 103 | 10–7 | 1 | 5 | 63 57 |
2914 ∞ |
25 $5 \times {{10}^{{ - 6}}}$ |
||
3.2 | 50 | 0 | 103 | 10–7 | 1 | 10 | 23 17 |
887 ∞ |
–5460 487 |
–5072 326 |
126 $3 \times {{10}^{{ - 11}}}$ |
3.3 | 50 | 0 | 103 | 10–7 | 1 | 15 | 12 7 |
386 ~ |
0.04 0.01 |
||
3.4 | 75 | 0 | 103 | 10–7 | 1 | 15 | 17 11 |
628 ∞ |
4 $2 \times {{10}^{{ - 10}}}$ |
||
3.5 | 85 | 0 | 103 | 10–7 | 1 | 15 | 20 14 |
752 ∞ |
8 | ||
3.6 | 85 | 0 | 103 | 10–6 | 1 | 15 | 4.6 –0.41 |
295 ∞ |
–1001 1507 |
–1317 128 |
250 $2 \times {{10}^{{ - 5}}}$ |
Таким образом, при ${{k}_{1}} = {{k}_{2}} = {{10}^{3}}$ и предположении, что вклад остальных диссипативных процессов меньше, верхнюю грань соответствующих коэффициентов можно оценить как
Отсюда видно, что от номера волны существенно зависит только коэффициент $r$.
Отметим два результата, полученных в эксперименте 3. Во-первых, увеличение номера гармоники (см. (3.1)–(3.3)) приводит к уменьшению скоростей течений в обоих слоях и уменьшению расхода. Во-вторых, вариант 3.6 показывает, что при некоторых условиях, (например, при увеличении коэффициента придонного трения) в нижнем слое может формироваться противотечение. Для проверки чувствительности решения к изменению модельных параметров было проведено две серии дополнительных расчетов. В первой серии варьировалась амплитуда рельефа дна (табл. 4). Расчеты показывают, что решение остается стационарным и устойчивым. Противотечение, слегка изменяясь по величине, в нижнем слое присутствует, его скорость $ \approx {\kern 1pt} - {\kern 1pt} 0.4$ см/с (табл. 4).
Таблица 4.
${{c}_{{15}}}$ | 81 | 82 | 83 | 84 | 85 | 86 | 87 | 88 | 89 | 90 |
${{V}_{{ACC}}}$ | 34.16 | 32.95 | 31.75 | 30.60 | 29.5 | 28.33 | 27.25 | 26.2 | 25.14 | 24.13 |
${{U}_{1}}$ | 4.78 | 4.73 | 4.68 | 4.63 | 4.6 | 4.54 | 4.50 | 4.45 | 4.40 | 4.36 |
U2 | –0.34 | –0.36 | –0.38 | –0.39 | –0.41 | –0.43 | –0.44 | –0.46 | –0.47 | –0.49 |
Во второй серии дополнительных расчетов варьировался коэффициент турбулентной вязкости КГВ (${{k}_{1}} = {{k}_{2}} = k$). Результаты представлены в таблице 5. Видно, что противотечение в нижнем слое возникает лишь при повышении $k$: $k > 920$$\,{{{\text{м}}}^{2}}{{c}^{{ - 1}}}$. Примерно, при $k \leqslant 930$ в обоих слоях течения совпадают с направлением ветра. Увеличение $k$ от $k \approx 930$$\,{{{\text{м}}}^{2}}{{c}^{{ - 1}}}$ приводит к уменьшению расхода и скоростей в верхнем и нижнем слоях. Скорость в нижнем слое направлена против ветра (отрицательна). Уменьшение коэффициента турбулентной вязкости от $k \approx 920$ приводит к увеличению расхода и скоростей в обоих слоях. Скорость в нижнем слое положительна.
Таблица 5.
$k \times {{10}^{3}}$, ${{{\text{м}}}^{2}}{{c}^{{ - 1}}}$ | 0.99 | 0.98 | 0.96 | 0.95 | 0.93 | 0.92 | 0.91 | 0.90 | 0.89 |
${{V}_{{ACC}}}$(св) | 32.86 | 36.32 | 43.53 | 47.19 | 54.91 | 58.93 | 63.21 | 67.46 | 71.82 |
${{U}_{1}}$, ${{U}_{2}}$, ${\text{м}}{{c}^{{ - 1}}}$ | 4.71 –0.36 |
4.83 –0.30 |
5.11 –0.19 |
5.25 –0.13 |
5.55 –0.015 |
5.72 0.04 |
5.90 0.01 | 6.07 0.17 |
6.26 0.23 |
Дополнительные расчеты показывают, что решения описывают стационарные предельные точки [4] и “непрерывно” зависят от возмущений рельефа и коэффициента турбулентной вязкости. Заметим, что в этой работе изучаются только стационарные режимы, наблюдаемые не для всех модельных параметров. Дополнительные расчеты подтверждают наличие устойчивых стационарных точек, описывающих режим противотечения.
7. ЗАКЛЮЧЕНИЕ
Сформулирована новая вариационная формулировка и построен численный алгоритм решения уравнений, описывающих нестационарную квазигеострофическую динамику в двухслойном периодическом канале, имитирующем АКТ. Алгоритм решения задачи основан на технике сопряженных уравнений и состоит в итерационном решении системы оптимальности.
Рассчитан и изучен ряд стационарных режимов квазигеострофических циркуляций в океане, возникающих при различных модельных параметрах. Показано, что увеличение номера гармоники, генерируемой гармоникой рельефа дна, приводит к уменьшению скоростей течений в обоих слоях и уменьшению полного расхода АКТ. Наличие в рельефе дна высоких гармоник может приводить к формированию противотечения в нижнем слое. Противотечение устойчиво к небольшим вариациям возмущений рельефа дна и коэффициента турбулентной вязкости.
При некоторых характеристиках рельефа дна (наличие высоких гармоник) и значениях коэффициента турбулентной вязкости, в нижнем слое может формироваться противотечение. Оно стационарно, направлено против ветра и устойчиво к небольшим изменениям входных модельных параметров.
Работа выполнена при финансовой поддержке РНФ, грант № 18-11-00163.
Список литературы
McWilliams J.C. Fundamentals of Geophysical Fluid Dynamics. Cambridge University Press, 2006. 299 p.
Дымников В.П., Залесный В.Б. Основы вычислительной геофизической гидродинамики. Москва: Геос, 2019. 448 с.
Дымников В.П., Грицун А.С. Ляпуновские показатели и размерность аттрактора двуслойной бароклинной модели атмосферы // ДАН СССР. 1996. Т. 347. № 4. С. 535–538.
Дымников В.П., Филатов А.Н. Основы математической теории климата. Москва: ВИНИТИ, 1994. 252 с.
Ипатова В.М. Задача инициализации для модели общей циркуляции атмосферы. Москва: Труды МФТИ. 2012. Т. 12. № 2. С. 121–130.
Agoshkov V.I., Ipatova V.M. Convergence of solutions to the problem of data assimilation for a multilayer quasigeostrophic model of ocean dynamics // Russ. J. Numer. Anal. Math. Modelling. 2010. V. 25. № 2. P. 105–115.
Bernier Ch. Existence of attractor for the quasi-geostrophic approximation of the Navier-Stokes equations and estimate of its dimension // Adv. Math. Sci. Appl. 1994. V. 4. № 2. P. 465–489.
Charney J.G., Shukla J., and K.C. Mo. Comparison of a barotropic blocking theory with observation // J. Atmos. Sci. 1981. V. 38. P. 762–779.
Gritsun A., Branstator G., Dymnikov V. Construction of the linear response operator of a atmospheric general circulation model to small external forcing // Rus. J. Num. Anal. Math. Modelling. 2002. V. 17. № 5. P. 399–416.
Gritsun A.S. Unstable periodic trajectories of a barotropic model of the atmosphere // Russ. J. Num. Anal. Math. Modelling. 2008. V. 23. № 4. P. 345–367.
Ивченко В.О., Залесный В.Б. Диффузионно-ротационная параметризация вихревых потоков потенциального вихря: баротропное течение в зональном канале // Изв. РАН. Физика атмосферы и океана. 2019. Т. 55. № 1. С. 3–16.
Ivchenko V.O., Zalesny V.B., Sinha B. Is the coefficient of eddy potential vorticity diffusion positive? Part 1: barotropic zonal channel // J. Phys. Oceanogr. 2018. V. 48. № 6. P. 1589–1607.
Dong-wook Shina, Younghun Kangb, Eun-Jae Parkb. C0-discontinuous Galerkin methods for a wind-driven ocean circulation model: Two-grid algorithm // Comput. Methods Appl. Mech. Engrg. № 328. 2018. P. 321–339.
Марчук Г.И. Избранные научные труды. Том II. Сопряженные уравнения и анализ сложных систем. Москва: РАН, 2018. 500 с.
Шутяев В.П. Методы усвоения данных наблюдений в задачах физики атмосферы и океана // Изв. РАН. Физика атмосферы и океана. 2019. Т. 55. № 1. С. 17–34.
Data assimilation for the Earth System. Swinbank R., Shutyaev V., Lahoz W.A. (eds.). Kluwer Academic Publishers, Dordrecht/Boston/London. 2003. 377 p.
Shutyaev V. Control operators and fundamental control functions in data assimilation. In R.Swinbank , (eds.). Data assimilation for the Earth System. 2003. P. 55–64.
Zalesny V., Agoshkov V., Shutyaev V., Parmuzin E., Zakharova N. Numerical Modeling of Marine Circulation with 4D Variational Data Assimilation / Journal of Marine Science and Engineering. 2020. 8. 503. 19 p.
Gilbert J.C., Lemarechal C.L. The modules M1QN3 and N1QN3. Version 2.0c (June 1995).
Дополнительные материалы отсутствуют.
Инструменты
Известия РАН. Физика атмосферы и океана