Численное моделирование явления термобара в озере Байкал Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 1, № 3, 1996

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЯВЛЕНИЯ ТЕРМОБАРА В ОЗЕРЕ БАЙКАЛ*"

О. Б. БОЧАРОВ, Т. Э. ОВЧИННИКОВА Институт водных и экологических проблем СО РАН Новосибирск, Россия

Численными методами моделируются возникновение и развитие конвективных течений под влиянием весеннего термобара в озере Байкал. Используется негидростатическая модель, в которой плотность зависит только от температуры. Анализируются различные краевые условия на открытой границе области. Показано, что возникающие глобальные циркуляции могут объяснить проникновение поверхностных вод в придонные области озера Байкал.

Термический бар был открыт Ф. Форелем более 100 лет назад [1] на Женевском озере. В 50-60-х годах нашего столетия изучение этого явления было продолжено А. И. Тихомировым на Ладожском озере и рядом американских лимнологов на Великих озерах. В основном это были натурные наблюдения и эксперименты в лотках. Теоретических исследований было крайне мало [2]. Интерес к проблеме возродили исследования последних лет на озере Байкал. С помощью современной аппаратуры было установлено, что придонные воды Байкала "моложе"вод основного глубинного ядра [3, 4]. Это означает, что должен существовать механизм глубокого перемешивания, переносящий поверхностные воды в придонные области. Одна из гипотез [4] состоит в том, что в периоды весеннего нагревания и осеннего охлаждения воды возникают циркуляционные потоки в районе термобара, образующегося в прибрежных участках озера, подверженных влиянию вод речных притоков. При этом особое внимание обращается на то, что на больших глубинах температура воды максимальной плотности падает с глубиной. Это дало толчок целому ряду исследований по проблеме термобара [2, 4-6].

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

*© О.Б.Бочаров, Т.Э.Овчинникова, 1996.

^ Работа выполнена при финансовой поддержке Международного научного фонда и Правительства России, совместный грант №КЫ1300.

1. Математическая модель и численный алгоритм

В основу исследования положена двумерная вертикальная негидростатическая модель в приближении Буссинеска, система уравнений которой имеет вид

Здесь и, т — составляющие скорости по осям х и г (направленным перпендикулярно берегу и по вертикали вверх соответственно), Т — температура, р0 — характерная плотность воды, коэффициенты х, К и Ох, характеризуют интенсивность переноса импульса и тепла в соответствующем направлении = 0.7Кд, д = х,г). Отметим, что уравнение состояния Кнудсена [7], использующееся здесь, не учитывает зависимости плотности от давления.

Коэффициент вертикальной турбулентной вязкости рассчитывался по формуле

На свободной поверхности г = Н ставится кинематическое граничное условие типа абсолютно гладкой неподвижной "твердой крышки", на дне задаются касательные напряжения, пропорциональные квадрату скорости. Поток тепла через свободную поверхность рассчитывается в зависимости от метеоусловий (потока солнечной радиации, температуры воздуха, облачности, влажности, атмосферного давления) и температуры воды. На дне ставится условие отсутствия теплообмена. Что касается граничных условий при х = 0 (у берега) и при х = Ь (открытая граница), то здесь мы исследовали три различных варианта с целью выявления их влияния на образование глобальных циркуляционных потоков.

Вариант 1. Твердая стенка при х = 0 и условие симметрии при х = Ь —

х = Ь : и = 0, тх = 0, Тх = 0. Вариант 2. Твердая стенка при х = 0 и условие свободного выхода струй при х = Ь —

х = 0 : и = 0, т = 0, Тх = 0;

х = Ь : их = 0, т = 0, Тх = 0.

Вариант 3. Теплый приток при х = 0 и условие свободного выхода струй при х = Ь (параметры притока по порядкам соответствуют р. Селенге) —

Р = Р(Т)-

(1)

х = 0 : и = 0, т = 0, Тх = 0;

х = 0: и = и0 = 2 ■ 10-2 м/с, т = 0, Т = Т0 = 10 °С;

х = Ь : их = 0, т = 0, Тх = 0.

х

Численный алгоритм нахождения поля скорости основан на так называемом Ф-проек-ционном методе [8], который в общих чертах выглядит следующим образом. На первом шаге рассчитывается вспомогательное поле скоростей V* из уравнений с исключенным градиентом давления:

дv* д д * д д д д р

— + — uv + — wv = — Кх— V + — Кг—V + — g. дъ дх дг дх дх дг дг р0

По этому полю определяется завихренность

rotv * = rotv = и = и* — w*х■

Функция тока вычисляется из уравнения

ДФ = и,

а окончательное соленоидальное поле скоростей — из уравнений связи по функции тока:

и=

д Ф

дг '

w = —

дФ

дх

Поля v * и Т определяются по неявным противопотоковым разностным схемам 1-го порядка точности. Алгебраические системы решаются итерационным методом верхней релаксации. Используется неравномерная ортогональная сетка с измельчением шагов у свободной поверхности (г = Н) и берега (х = 0). В результате Нх меняется от 1 до 33 м, а кх — от 33 до 1350 м.

2. Численные эксперименты

В численных расчетах мы попытались воспроизвести реальные условия озера Байкал. Прибрежный профиль озера взят из работы [9], а глубина Н = 900 м примерно соответствует средним глубинам южного бассейна Байкала. Протяженность расчетной области Ь =10 км. Начальное вертикальное распределение температуры (рис. 1) взято на основе данных натурных наблюдений в окрестности мыса Лиственичного в мае 1970 г. [ 10, с. 71]. Метеоусловия, соответствующие маю—июню, взяты из работы [11]. При расчете турбулентной вязкости коэффициент А и масштаб турбулентности Ь в формуле (1) использовались для достижения достаточно близкого соответствия вычисленных значений Кх данным натурных измерений из [10]. Следует отметить, что в разных численных экспериментах вычисленные значения коэффициента турбулентной вязкости отличались весьма мало. На рис. 2 приведены графики характерного вертикального распределения значений этого коэффициента на трех различных вертикалях.

Численные эксперименты проводились для трех описанных выше вариантов граничных условий. Начальное условие для всех вариантов расчета — покой.

На рис. 3 показана картина линий тока, сформировавшаяся через 30 суток после начального момента при первом варианте граничных условий (симметрия относительно открытой границы). Эта картина качественно стабильна в том смысле, что довольно быстро образуются четыре крупные конвективные ячейки, размеры которых изменяются очень медленно. Это позволяет поверхностным водам проникать на большие глубины. Течения здесь достигают самого дна.

Рис. 1. Начальное вертикальное распределение температуры.

2 00-

□ 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0

X, КМ

Рис. 4. Линии тока через 5 суток (вариант 2). Ф = -0.02 (1), -0.1 (2), -0.2 (3), 0 (4), 0.1 (5), 0.5 (6), 1.0 (7).

□ 1.0 2 . О 3 . О 1 . О 5. О 6. О 7 . О 8 . О 9. О

X, км

Рис. 5. Линии тока через 20 суток (вариант 2). Ф = -0.07 (1), 0.07 (2).

—.—|—.—|—.—|—.—|—.—|—.—|—.—|—.—|—I—|—.—

О 1.0 З.О З.О О О 6. О . О В. О 9. О

X, км

На рис. 4-6 показана динамика развития конвективных потоков для второго варианта граничных условий (через 5, 20 и 25 суток от начального момента). Здесь также образуется несколько крупных вихрей, которые медленно продвигаются к открытой границе, а вблизи берега образуются новые. Конвективные течения достигают достаточно больших глубин (порядка 700 м).

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

—I—|—I—|—.—|—.—|—I—|—I—|—I—|—I—|-

О 1.0 2.0 3.0 4.0 5.0 6.0 7.0 3.0

X, КМ

Рис. 7. Линии тока через 7 суток (вариант 3).

Ф = -0.2 (1), 0.2 (2).

Во всех расчетных вариантах скорости вдоль склона имели порядок 10-2 см/с, что соответствует наблюдениям В. И. Верболова [12] и его оценкам темпов обновления глубинных вод в озере.

Приведенные результаты были получены в условиях равенства коэффициентов вертикального и горизонтального турбулентного обмена. Так как масштабы процессов по вертикали и горизонтали существенно отличаются, то фактически реализуется как бы приближение пограничного слоя, когда горизонтальная турбулентная вязкость не работает. С этим, по нашему мнению, связано количество крупных вихрей в рассматриваемой области. Данных натурных наблюдений для оценки их количества пока не хватает, а теоретический анализ еще предстоит сделать. Это же касается и анализа граничных условий на открытой границе: хотя из приведенных рисунков можно сделать вывод, что лучше всего в данной задаче работают условия симметрии (вариант 1), тем не менее окончательный ответ требует дополнительной работы с натурными данными.

Одновременно мы анализировали возможность применения в этой задаче гидростатической модели. График относительного отклонения от гидростатического приближения f =| 1 + 1/рд • дР / дх | для первого варианта граничных условий, приведенный на рис. 8, показывает, что в большей части области течения оно очень мало и только в узкой зоне нисходящего потока достигает значений порядка 10-4, причем для остальных вариантов условий это отклонение было еще меньше. Поэтому применение модели в гидростатическом приближении здесь вполне допустимо.

На основании проведенных исследований можно сделать следующие выводы. В период весеннего прогрева в результате возникновения термобара формируются крупномасштаб-

Рис. 8. График функции f =

1 +

1 dP рд dz

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

Список литературы

[1] ФОРЕЛЬ Ф. А. Руководство по озероведению (общая лимнология). Типография В. Ф. Киршбаума, СПб., 1912.

[2] ZiLiTiNKEViCH S.S., Kreiman K.D., TERZEVIC A. Yu. The Thermal Bar. J. Fluid Mech, 27, 1992, 27-42.

[3] IMBODEN D. M. Deep Water Formation: the Physical Mystery of Lake Baikal. В "Байкал — природная лаборатория для исследования изменения изменений окружающей среды и климата: Тез. докл.", ЛИСНА, Иркутск, 1994.

[4] Гранин Н. Г., Шимараев М. Н. К вопросу о стратификации и механизме конвекции в Байкале. Докл. АН СССР, 321, №2, 1991, 381-385.

[5] ЗИЛИТИНКЕВИЧ С. С., КРЕЙМАН К. Д. Теоретическое и лабораторное исследование термического бара. Океанология, 30, вып. 5, 1990, 750-755.

[6] ЦВЕТОВА е. А. Математическое моделирование Байкальского термобара. В "Мат. проблемы экологии: Тез. докл. Всерос. конф.", Ин-т матем. СО РАН, Новосибирск, 1994, 53-57.

[7] ЛАКОМБ А. Физическая океанография. Мир, М., 1974.

[8] Тимухин Г .И. О построении некоторых схем приближенного решения уравнений динамики вязкой несжимаемой жидкости. Числ. методы механ. сплошной среды, 3, №2, 1972, 77-95.

[9] Shimaraev M. N.,Verbolov V. I., Granin N.G., Sherstyankin P.P. Physical Limnology of Lake Baikal: a Review. Irkutsk—Okayama, 1994.

[10] ШиМАРАЕВ М. Н. Элементы теплового режима озера Байкал. Наука, Новосибирск, 1977.

[11] Россолимо Л. Л. Температурный режим озера Байкал. Тр. Байкальской лимнологической станции, 16, Изд-во АН СССР, М., 1957.

[12] ВЕРБОЛОВ В. И. Вертикальные присклоновые циркуляции и вертикальный водообмен. В "Течения в Байкале", Наука, Новосибирск, 1977, 122-128.

Поступила в редакцию 15 сентября 1995 г.