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

Во введении раскрывается актуальность диссертационной работы, пред- ставлены основные характеристики объекта исследования, сформулированы цель и задачи исследования, основные положения, выносимые на защиту, обосновыва- ется научная новизна, практическая значимость и достоверность полученных ре- зультатов, изложена методология исследования.
Озеро Шира представляет собой бессточное озеро без островов, в которое впадает речка Сон. В силу малости притока все влияние реки сосредоточено в при- устьевой зоне, основным внешним фактором, определяющим течение в озере, яв- ляется ветровое воздействие (рис. 1).
Рис. 1 – Сила и направление ветра на озере в 2014 г. Данные с метеостанции пос. Шира
Основные характеристики: озеро имеет овальную форму, длина – 9,4 км, ши- рина – 5 км, площадь водного зеркала – 34,7 км 2 , средняя глубина – 11,2 м, мак- симальная глубина – 25 м. Батиметрия озера представлена на рис. 2, а.
абв
Рис. 2 – Батиметрия озера Шира (а), температура – сплошная синяя линия, соленость –
штриховая красная линия (б), график плотности по формуле Буссинеска4 (в)
На рис. 2, в представлен график зависимости плотности от глубины, полученный на основании обработки локальных измерений температуры и кондуктивности, который наглядно демонстрирует сильную стратификацию воды
4 Boussinesq, J. Théorie analytique de la chaleur mise en harmonic avec la thermodynamique et avec la théorie mécanique de la lumière: Tome I — Gauthier-Villars, 1903.
8
в озере Шира в течение всего летнего периода (измерение 13.07.2015 в точке глу- биной 22 м).
Весной и осенью озеро обычно является меромиктическим, т.е. в циркуляции участвует только поверхностный слой жидкости, что приводит к специфическому составу донных отложений. При смене меромиктического режима на голомикти- ческий (хотя бы один раз за год происходит перемешивание воды до дна) может также происходить нарушение баланса биоты. Процесс перемешивания воды в озере до дна прошел в озере Шира осенью 2014 – зимой 2015 г., что было опреде- лено по косвенным признакам весной 2015 г. Пока сложно сказать почему меро- миктическое озеро изменило свой тип функционирования на другой, одно из воз- можных объяснений через анализ глубины залегания весеннего термоклина при- ведено в книге5.
В первой главе диссертации представлена математическая модель ветро- вого движения неоднородной жидкости, приведены новые аналитические решения задачи ветрового движения однородной и на их основе определено взаимодействие между моделями различного уровня сложности.
В параграфе 1.1 рассматривается математическая модель трехмерного тече- ния неоднородной жидкости в приближении Буссинеска и предположении о гидростатичности давления в замкнутом водоеме, на основе которой выполнено аналитическое и численное моделирование течений в озере Шира.
Система уравнений содержит следующие уравнения: уравнения движения –
u u u u  u  u  u  g0 (1) tuxvywzlv=zKz zxKh xyKh ygx xdz,
 0z
v v v v  v  v  v  g0 (2) tuxvywzlu=zKzzxKhxyKhygy ydz;
 0z уравнение неразрывности несжимаемой жидкости –
u  v  w = 0; (3) x y z
уравнения переноса тепла и солей –
u  u T  v T  w T =   K T T     K T T     K T T  , (4)
t x y z x h x y h y z z z 
S  u S  v S  w S =   K S S     K S S     K S S ; (5) t x y z x h x y h y z z z

уравнение состояния в общем случае –
 = (T,S). (6)
Здесь x, y, z – прямоугольная система координат, ось Ox направлена на восток, ось Oy – на север, ось Oz – вертикально вверх; (u,v,w) – вектор скорости течения, u=u(x,y,z,t), v=v(x,y,z,t),w=w(x,y,z,t) – восточная, северная и вертикальная
5 Рогозин, Д. Ю. Меромиктические озера Северо-Минусинской котловины: закономерности стратификации и экология фотофторных серных бактерий — Красноярск: ИФ СО РАН, 2019.
компонентысоответственно; =(x,y,t)–отклонениесвободнойповерхностиот невозмущенного положения; T – температура воды, отсчитываемая от некоторого
среднего значения T0 ; S – соленость воды; Kh , Kz – коэффициенты горизонтального и вертикального турбулентного обмена для скорости; KhT , KhS ,
K zT , K zS – коэффициенты горизонтальной и вертикальной диффузии для температуры и солености соответственно; l – параметр Кориолиса; g – ускорение свободного падения;  – плотность воды в точке с координатами (x,y,z); 0 – некоторое среднее значение плотности.
Система уравнений дополняется граничными условиями.
Рис. 3 – Границы водоема На свободной поверхности при z = (x,y,t) (рис. 3):
 K Vh   w , w    u   v  , T = Q , S = Q . 0zz ztxyz Tz S
z z
Здесь V (u,v) – вектор горизонтальной скорости; w (w,w)– напряжение
hxy
ветра на свободной поверхности; Q , Q – потоки тепла и солености через TS
свободную поверхность.
На дне водоема (рис. 3) возможно несколько постановок граничных условий
для скорости, в частности используется нелинейное или линейное условие про- скальзывания.
Для вертикальной составляющей скорости на дне всегда выполняется: w zH = u zH H x  v zH H y.Для температуры и солености задаются усло-
вия: T z  Q , S z = Q . TS
На вертикальных границах 1, 2 (рис. 3) ставятся условия прилипания u  v  0 или непротекания: равенство нулю нормальной составляющей вектора го-
ризонтальной скорости Vh. Для температуры и солености полагаем T N =Q , S N =Q , где N – нормаль к вертикальным границам.
TS
Уравнения (1) – (6) и граничные условия конкретизируются в зависимости от решаемой задачи.
В параграфе 1.2 рассматриваются аналитические решения для модели Экмана ветрового движения жидкости с условием прилипания и проскальзывания на дне. Представлены полученные автором решения в случае более часто используемого в численных расчетах условия проскальзывания на дне при постоянном и переменном коэффициентах вертикального турбулентного обмена. Соответствующие решения сравниваются с известными решениями для модели
Экмана с условиями прилипания на дне. В параграфах 1.3 и 1.4 приведены аналитические решения для модели с учетом горизонтального обмена в двумерном (течение в вертикальной плоскости) и трехмерном случаях. Решения найдены при постоянном и переменном коэффициентах вертикального турбулентного обмена. В трехмерном случае решение при переменном коэффициенте вертикального турбулентного обмена найдено для дрейфовой составляющей течения.
Вторая глава посвящена численным расчетам для озера Шира с использованием системы моделирования ROMS. В параграфе 2.1 приведено краткое описание численной модели, реализованной открытым кодом ROMS. В параграфе 2.2. представлена адаптация численной модели, реализованной откры- тым кодом ROMS, к расчету течений в озере Шира, которая включала в себя: со- здание цифровой модели рельефа дна по пространственно-разбросанным данным измерений батиметрии с использованием метода интерполяции бигармониче- скими сплайнами; задание начального профиля температуры и солености; подклю- чение метеорологических данных для учета ветрового воздействия и потоков тепла. Представлены результаты верификации численных алгоритмов ROMS на полученных аналитических решениях.
Для численного моделирования гидродинамики водоема необходимо опре- делить его рельеф дна в виде непрерывно-дифференцируемого поля. В настоящем параграфе были рассмотрены ранее известные методы построения полей на регу- лярных сетках из пространственно-разбросанных данных – в нашем случае это ба- тиметрия, измеренная эхолотом с катера, оборудованного GPS – и разработана ме- тодика построения цифровой модели рельефа дна6, которая включала в себя: пред- варительную обработку данных батиметрии для выявления и исключения проти- воречивых данных; интерполяцию данных батиметрии на равномерную сетку с ис- пользованием бигармонической сплайн-интерполяции7; бигармоническое сглажи- вание функций Грина в методе бигармонической сплайн-интерполяции. Результаты применения разработанной методики для озера Шира представлены на рис. 4.
аб
Рис. 4 – Данные измерений глубины в озере Шира (а), поле данных батиметрии, построенное
на регулярной сетке ∆ = ∆ = 10 м (б)
Изложены результаты верификации модели ROMS на полученном в Главе 1
аналитическом решении для модели с учётом горизонтального турбулентного об-
6 Shchepetkin, A.F., Volodko, O.S. A method for reconstructing bottom topography for an enclosed basin from scattered, sparse measurement data. Океанологические исследования, 46(3), 2018, P.67–84.
7 Sandwell, D.T. Biharmonic spline interpolation of GEOS‐3 and SEASAT altimeter data. Geophysical research letters, 1987, 14(2), P.139–142.
11

мена (рис. 5). Параметры данного расчета: длина и ширина прямоугольного бас- сейна 10000 м, ветер восточный  0 104 м2 с2 , условие прилипания на дне, ко-
эффициент вертикального турбулентного обмена K  0,02 м2 с, коэффициенты z
 1 м2 с , глубина H  50 м.
аб
Рис. 5 – Годографы скоростей в центральной точке бассейна: а – модель с учетом
горизонтального турбулентного обмена (ромбики), решение для модели Экмана (непрерывная линия), б – результаты расчета в модели ROMS
Масштабы скоростей совпадают, различие в поведении на дне объясняется тем, что в модели ROMS скорости выводятся в точке, находящейся на расстоянии половина шага сетки в вертикальном направлении. Верификация показала доста- точно хорошее совпадение результатов численных расчетов с аналитическим ре- шением.
В параграфе 2.3 представлены результаты исследования свойств разностных схем, применяемых для расчета в численной модели ROMS.
При численном решении системы уравнений (1) – (6) методом конечных раз- ностей основные проблемы возникают при аппроксимации адвективных членов. В модели ROMS для аппроксимации имеется возможность выбора этих слагаемых. Сравнительный анализ различных разностных аппроксимаций адвективных чле- нов, применяемых для решения уравнения переноса и диффузии тепла и солей (4)- (5), был проведен на аналитических решениях линейного уравнения конвекции – диффузии в одномерном и двумерном случаях.
В общем случае решение для уравнения конвекции-диффузии с постоянными коэффициентамиводномерномслучае C tuC xK 2C x2 можетбыть
горизонтальной вязкости K
h
exp{
где С0 (x) – начальные данные. Для четырех начальных конфигураций: непрерыв-
ных – функция Гаусса, треугольник и разрывных – «полочка» и «помост» удается взять интеграл, который во всех случаях выписывается через экспоненциальную функцию и/или функцию ошибки erf. Это дает возможность провести сравнитель- ный анализ разностных аппроксимаций на решениях различной степени гладкости.
Для аппроксимации диффузионного члена использовалась обычная трехто- чечная аппроксимация. Численные расчеты были проведены для двенадцати вари- антов явной аппроксимации адвективных членов, в том числе для используемых в
получено с использованием формулы С(x, t)  0
 2(Kh)
}d,
 C (x)
(xut)2
h
12
1/2
4Kht

ROMS аппроксимаций вида предиктор-корректор со схемами вида Leap-Frog в ка- честве предиктора и параметром  , отвечающем за метод интерполяции значений на трех слоях по времени в корректоре. В предикторе первая производная аппрок- симировалась по пространству центральной разностью 2-го и 4-го порядка и раз- ностями против потока 2-го и 3-го порядка.
В случае 2-х пространственных переменных рассматривается уравнение C tuC xvC yK 2C x2 2C y2, где K –постоянныйкоэффи-
h
циент диффузии вдоль осей Оx и Оy. Скорости переноса u   py,v  px, p  const.
h
В начальный момент времени задается распределение C0(x,y,0). В этом случае
осуществляется два процесса: процесс вращения начальной конфигурации вокруг начала системы координат и диффузия.
В двумерном случае аналитическое решение уравнения диффузии строится как произведение решений соответствующих одномерных задач. На рис. 6, а, б представлен график, когда начальные данные – пирамида высоты b, в основании которой лежит квадрат со стороной 2l, с центром в точке x0 , y0 . Расчет прово-
дился при следующих параметрах: расчетная область со стороной квадрата L  10 с шагом по обоим направлениям x  y  0,5, шаг по времени t  0,01, на всех
границах ставилось условие свободного прохода, p  2, Kh  0,05. На рис. 6 (в – е)
приведены сечения аналитического и разностного решения вертикальной плоско- стью, проходящей через биссектрису первого и третьего квадрантов. Все решения показаны на момент времени t  3,14, когда начальное распределение осуществ- ляет полный оборот вокруг начала системы координат.
аб
вгде
Рис. 6 – Начальные данные пирамида при b  1, l  1 с центром в точке x0  1, y0  1 и
аналитическое решение соответственно (а, б). Результаты расчетов для схем вида предиктор-
корректор с Leap-Frog в качестве предиктора и   1 12 в корректоре с различными
аппроксимациями первой производной по пространству: в –центральная разность 2-го порядка, г – центральная разность 4-го порядка, д – разность против потоков 2-го порядка, е – разность против потоков 3-го порядка
В результате этого исследования по оценке величины абсолютной ошибки было определено преимущество аппроксимации четвертого порядка для адвектив- ных членов в предикторе Leap-Frog, которая в дальнейшем использовалась в рас- четах.
13

В параграфе 2.4 представлены результаты численных расчетов на основе мо- дели ROMS c модельным ветром и учетом реальных метеорологических данных в летний период. Расчеты проведены с гидрофизическими характеристиками озера Шира, в том числе распределения температуры и солености, типичной для летнего периода, нулевых начальных скоростях и возвышении свободной поверхности, различных ветровых сценариях при следующих значениях параметров: шаг по го-
ризонтали x 15 м, y  20 м,
по вертикали z  1 м
раллеливания OpenMP, MPI. Начало расчета – покой, коэффициент придонного
, применялись методы распа- трения Kb  104 м с , коэффициент горизонтального турбулентного обмена
 2 м2 с , линейное уравнение состояния. Использовалась KPP-параметризация
турбулентного вертикального перемешивания.
Были рассмотрены различные сценарии, в том числе мгновенно начинаю-
щийся постоянный ветер, который мгновенно заканчивается. В этом случае после остановки восьмичасового ветра северо-восточного, юго-восточного, восточного и южного направлений со скоростью 8–10 м/с визуально наблюдаются колебания свободной поверхности с периодом 11 и 6 часов. Движение происходит против ча- совой стрелки и максимальные значения достигаются на противоположных бере- гах бассейна, сама свободная поверхность также подвергается колебаниям мень- шего масштаба (рис. 7).
Рис. 7 – Колебания свободной поверхности после прекращения действия восьмичасового северо-восточного ветра со скоростью 8 м/с
В северной и центральной географических точках озера наблюдаются незна- чительные колебания изотерм и изогалин, в южной, восточной и западной колеба- ния достигают 3-х метров по глубине (рис. 8).
аб
Рис. 8 – Колебания изотерм (а) и изогалин (б) в западной части озера после воздействия
северо-восточного ветра со скоростью 8 м/с в течение 8 часов, далее 24 часа – без ветра
В случае западного, юго-западного, северо-западного и северного направле- ний ветра возникают сейшевые колебания свободной поверхности без вращения. В изотермах, построенных по натурным наблюдениям (рис. 9), во время прохож- дения внутренней волны также наблюдаются колебания с размахом от трех до пяти метров.
K
h
14

Рис. 9 – Географическая точка T1, изотермы за период 08.07.2015–11.07.2015, T=10–18C Представлены результаты проведенных впервые сезонных расчетов для озера Шира с учетом реальных метеорологических данных (температура и влаж- ность воздуха, осадки, данные о ветре получены с метеостанции пос. Шира). Рас- четы были проведены для периода 01.06.2018–31.08.2018. Картина течения на
13.07.2018 г. в 2 часа ночи, в 6 и в 8 часов утра представлена на рис. 10.
Рис. 10 – Возвышение свободной поверхности
Период внутренней волны, визуально определяемый по изменениям возвы- шения свободной поверхности, составляет 7 часов (рис. 10). В период прохождения внутреней волны наблюдалось уменьшение температуры в восточной мелководной зоне водоема, после прохождения температура увеличивается.
В параграфе 2.5 излагаются результаты спектрального анализа данных измерений скорости, соответствующие летнему периоду 2014–2015 гг. и скорости, полученной при численных расчетах с использованием модели ROMS с учетом метеорологических данных соответствующих летнему периоду 2018 г. Для нахож- дения спектральных оценок были выбраны ряды восточной и северной составляю- щих трехмерной скорости течения на глубине 7,6 м, усредненные по времени из- мерения раз в полчаса. При построении графиков спектральной плотности мощно- сти использовался метод модифицированных периодограмм Уэлча с весовой функцией Блэкмена-Харриса (рис. 11).
абв
Рис. 11 – Спектральные характеристики: а, б – по измерениям ADCP 600 кГц, в – при анализе
численных расчетов, выполненных с использованием модели ROMS для летнего периода с 13.07.2018 – 31.07.2018 г. Для восточной компоненты скорости – черная штриховая линия, для северной компоненты скорости – красная сплошная линия, глубина 7,6 м
15

Спектральный анализ показал, что внутренняя волн с периодом 7 часов наблюдается как в натурных измерениях скорости, так и в численных расчетах.
Энергонесущие максимумы восточной и северной составляющих течения в 2015 г. существенно не изменились по сравнению с 2014 г. и составляют примерно – 21, 11, 7 часов. Этот результат можно сопоставить с отмеченным в наблюдениях зимой 2014–2015 гг. явлением полного перемешивания воды, что было зафиксиро- вано измерениями химического и биологического состава воды. Таким образом показано, что спектры течения являются устойчивой характеристикой озера.
Третья глава посвящена обработке и анализу данных долговременных из- мерений (в течении 1,5 месяцев) температуры и скорости в озере Шира в летний период 2013-2018 гг. с использованием методов статистической обработки данных, в том числе метода главных компонент (эмпирических ортогональных функций (ЭОФ)).
Измерения температуры и скорости проводились в летний период 2009 – 2018 гг. с использованием термоплетей и акустических доплеровских про- филографов скорости (ADCP 600 кГц и ADCP 1200 кГц ), позволяющих проводить непрерывный мониторинг гидрофизических параметров по всей глубине.
Параграф 3.1 посвящён обработке и анализу данных измерений температуры в озере в зависимости от глубины и поверхностной температуры с помощью
вещественного метода ЭОФ. Для анализа данные измерений температуры T zi ,tk , где zi(i1,…,N) и tk k 1,…,K – переменные по пространству и времени
соответственно, представляются в виде суммы произведений двух сомножителей, один из которых зависит от времени, а второй – от пространства. Стандартным образом определяется матрица коэффициентов корреляции
Rzi,zj 1 KK Tzi,tk Tzj,tk  и находятся ее собственные значения n и k1
собственные векторы n (zi ) , которые называются ЭОФ или модами. Поле температур T представляется в виде
Tzi,tk N En tk n zi , En tk N n zi Tzi,tk , (7) n 1 i 1
где En tk  – модальные коэффициенты или главные компоненты. Слагаемые
ранжированы в соответствии с убыванием положительных собственных значений n.
С использованием разложения (7) получены следующие результаты:
1. Первые моды, содержащие около 86% общей дисперсии, для семи различных точек измерения положительны и совпадают ниже глубины 8 м. Это свидетельствует о том, что вдоль всей акватории озера, процессы, связанные с распределением температуры по глубине, имеют однотипный характер.
2. Характер изменения модальных коэффициентов En tk  позволяет определить
время возникновения внутренних волн.
3.Для всех долговременных измерений температуры (2013, 2015, 2018 гг.) выявлена сильная линейная связь между первым модальным коэффициентом
E1 tk  разложения поверхностной температуры в ряд по ЭОФ и температурой
воздуха.
В параграфе 3.2 представлен анализ долговременных измерений скорости с
помощью комплексного метода эмпирических ортогональных функций, который позволил определить структуру основного течения в озере: для всех долговремен- ных измерений скорости получены однотипные картины поведения первой моды (рис. 13) аналогичные спиралям в случае однородной жидкости, показанным на рис. 5.
Рис. 13 – Годографы первой моды для четырех длительных измерений
Трактовка первой моды как стационарного течения однородной жидкости позволяет оценить коэффициент вертикального турбулентного обмена. В резуль- тате анализа измерения скорости 01.07.2015 – 07.08.2015 получено
Kz 10,6104 м2 с.
В Заключении сформулированы основные результаты настоящей работы. 1. Получены аналитические решения для упрощенных вариантов модели трехмер- ного течения неоднородной жидкости в приближении Буссинеска и предположе- нии о гидростатичности давления для стационарного медленного течения одно- родной жидкости как для модели Экмана, так и для модели с учетом горизонталь- ного турбулентного обмена в случае двумерного течения (в вертикальной плоско- сти) и трехмерного течения с условием проскальзывания на дне при постоянном и переменном коэффициентах вертикального турбулентного обмена. Эти решения обобщают ранее известные результаты для условий прилипания на дне. На основе исследования полученных аналитических решений сделан вывод о применимости более простой модели Экмана на удалении от берегов и определен тип граничных условий на дне (условие проскальзывания) для проведения численных расчетов в озере Шира.
2. Выполнено численное моделирование гидрофизических процессов, происходя- щих в озере Шира, с использованием численной модели, реализованной открытым кодом ROMS и разработанной методике построения цифровой модели рельефа дна с использованием пространственно-разбросанных данных измерений батиметрии, позволяющей получить непрерывно-дифференцируемое поле батиметрии озера. Верификация численной модели ROMS осуществлена на полученных частных ана- литических решениях. Валидация математической модели выполнена путем каче- ственного и количественного сравнения полученных результатов численных рас- четов с данными натурных измерений в летний период, в частности, по характеру изменения температуры и оценкам спектральных характеристик скоростей тече- ния.
17

3. Выявлена зависимость характера возникающих внутренних волн от силы и направления ветра на основании проведенных модельных расчетов с учетом мо- дельного ветра и расчетов с учетом реальных метеорологических данных за весь летний период.
4. На основе численного моделирования и спектрального анализа внутренних волн в озере Шира сделан вывод об устойчивости гидрофизических характеристик в пе- риод потери озером свойства меромиктичности.
5. С использованием метода главных компонент (эмпирических ортогональных функций) для данных долговременных измерений температуры найден простой способ определения моментов возникновения внутренних волн и получена коли- чественная оценка влияния температуры воздуха на поверхностную температуру. Применение этого метода к данным долговременных измерений скоростей течения в озере позволило определить структуру течения в эпилимнионе (верхнем прогре- том слое) и получить оценку величины коэффициента вертикального турбулент- ного обмена.