Численное моделирование годовой динамики вертикальной структуры соленого озера Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

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

Том 13, № 4, 2008

Численное моделирование годовой динамики вертикальной

>к

структуры соленого озера

В. М. Белолипецкий, С.Н. Генова Институт вычислительного моделирования СО РАН,

Сибирский федеральный университет, Красноярск, Россия e-mail: [email protected]

A 1D (in vertical direction) mathematical model of the temperature and salinity regimes of a lake is discussed. The model accounts for the ice formation process. A simplified model for estimation of the dynamics of the ice cover thickness is offered which accounts for a convective mixing layer. Process of ice thawing on the both lower and upper surfaces of an ice cover during springtime is modeled. The calculated isolines of salinity and temperature accords well with the observations for all seasons.

Введение

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

1. Одномерная модель температурного и солевого режимов водоема в период отсутствия ледяного покрова

Формирование температурного режима в непроточном стратифицированном водоеме определяется ветровыми течениями и теплообменом с атмосферой. Задача для температуры формулируется в виде [1]:

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, (грант № 07-01-00153_а), РФФИ-НВО (проект № 05-05-8902 НВО), междисциплинарного интеграционного проекта СО РАН № 24.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

дТ _ _д (К —

Ж _ дг \ т Ж"

с граничными условиями

К дТ

Кт^~

дг

2=0

1 п

Срр0

+ ав

К дТ

Кт^~

дг

Р/ е

—в*

срр0

*=Я

я

сРр0

(1)

(2)

где Т — температура воды; —т (г) — коэффициент вертикального турбулентного обмена для температуры; Ря — теплообмен с дном; Рп — полный тепловой поток через свободную поверхность; Р/ — приходящая коротковолновая радиация; в — коэффициент поглощения излучения; а — параметр, определяющий часть коротковолновой радиации, проникающей на глубину (0 < а < 1); ср — удельная теплоемкость воды; р0 — характерное значение плотности воды; Н — глубина водоема.

Аналогично ставится задача для определения вертикального распределения солености:

35 _ д ді дг К

дг

- ^ I —я £

дг = —Ря

2 = 0

К

дг

(3)

_ — ^кя

*=Я

Здесь 5 — соленость воды, К я (г) — коэффициент вертикального турбулентного обмена для солености, Ряя — массообмен с дном, Ря — поток через свободную поверхность. Необходимо также задать начальные распределения температуры и солености:

Т(0,г)_ Т0(г), 5(0,г)_ 50(г).

Для пресной воды плотность зависит только от температуры, уравнение состояния соленой воды принимается в приближении Буссинеска [2]:

р _ Р0

Т 5

Є1 + Є2^Т + Є37Г

Т0 50

(4)

где Р0 _ 1.0254 г/см3, Єї _ 0, 9753, Є2 _ —0.00317, Єз _ 0, 02737, Т0 _ 17.5 °С, 50 _ 35%о.

Существенное влияние на тепломассоперенос оказывает турбулентность. Для параметризации вертикального турбулентного обмена применяется формула, полученная на основе формулы Прандтля—Обухова и приближенного решения Экмана для ветровых течений [3]:

(0.05^1)2В + Ктіп при В > 0,

Ктіп при В < 0.

К*

(5)

Здесь В

| _____ І е—2а*

\Р0—0 /

р0

др

дг

л/т2’ +Т

напряжение трения ветра,

Ктіп _ 0.02 см2/с — минимальное значение коэффициента вертикального турбулентно-

(°.°5п)2т го обмена, К0 _ — ---------—, а

/

^ї _ —0, / — параметр Кориолиса.

2р0/ ’ У2—о V 2/

В работе [4] отмечается, что коэффициенты диффузии тепла и соли в воде различаются приблизительно в 100 раз, а коэффициент диффузии тепла меньше коэффициента вязкости в семь раз. В данной работе полагаем, что

—т _ —я _ —* при 0 < г < Л4,

Кт _ 0.1—*, —я _ 0.1—т при Л4 < г < Н.

Модуль напряжения трения ветра на водной поверхности рассчитывается по формуле [5]:

т = р„(0.69 + 0.107 ■ ИЪ) ■ 10-3И22,

где ра — плотность воздуха, И2 — скорость ветра на высоте 2 м (в м/с).

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

= (1 — а)Г1 — (Гэф + ^ (6)

где РЭф — эффективное длинноволновое излучение, — теплоотдача испарением, — конвективный теплообмен, — приходящая коротковолновая радиация. Составляющие теплового обмена через поверхность водоема вычисляются по известным соотношениям [6].

Численный алгоритм решения задач (1)-(3) основан на неявной схеме и методе прогонки.

2. Параметризация вертикальной структуры водоема в зимний период

2.1. Модель изменения толщины льда

Для определения динамики толщины ледяного покрова применяется упрощенная модель, основанная на квазистационарном температурном режиме в затвердевшей области

[6]. В этом случае из баланса тепловых потоков на границе раздела вода—лед с учетом скрытой теплоты фазового перехода получается уравнение

Ьрл^ = АлГф-^ - Я— - (О, (7)

2 [1 - (1 + К^)] еХР(-К^С) ,

/ (£) = —----------------------------------------------^-ехР(-к\/£);

к £

тепловой поток Гв—л на границе раздела вода—лед определяется по соотношению

= А Г(й) - Гф

Г в—1 - АВ *

Здесь и на рис. 1 £ — время, г — вертикальная координата (ось Ог направлена вниз), £ — толщина льда, й — глубина распространения конвекции, Н — глубина водоема, Ав — коэффициент теплопроводности воды, Ал — коэффициент теплопроводности льда, рл — плотность льда, Ь — удельная теплота плавления льда, Г — температура воды, Тл — температура поверхности льда, Гф — температура кристаллизации воды, зависящая от солености воды Б, ехр(—0*4£с), £с — толщина слоя снега, рл = ГС ехр(-кг),

к — коэффициент ослабления солнечной радиации в ледяном покрове (к = 0*17). Для определения ослабления солнечной радиации в слое льда использовалась эмпирическая зависимость [7].

Температура кристаллизации соленой воды определяется по следующей формуле [8]:

Тф = -0.0575Б + 1*710523 ■ 10-3Б3/2 - 2*154996 ■ 10-4Б2*

Рис. 1. Схема вертикальной структуры озера для периода нарастания льда

Численное решение уравнения (7) находится из соотношения

с+1 = л/(е)2 + 2(ап - Ъп)АЬ - 2^с£п/(е)Л£/(£рл),

А (т т ) рп рп

где ап = ------л), Ъп = , £п = Ш, ¿п+1 = Ьп + А1

^рл Ьрл

Если на поверхности льда имеется снег, то температура поверхности льда Тл выражается через температуру снега Тс [9]:

Тл = Тс

1 +

£сАл

£Ас

-1

Здесь £с — толщина слоя снега, Ас — коэффициент теплопроводности снега.

2.2. Оценка толщины слоя конвективного перемешивания

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

Для построения упрощенной математической модели выделим четыре горизонтальных слоя жидкости под поверхностью ледяного покрова (рис. 1):

£ < г < ¿1 — слой конвективного перемешивания;

¿1 < г < г2 — слой, в котором происходит конвективное перемешивание;

г2 < г < гз — внутренний слой;

дТ гдБ п

г3 < г < Н — придонный слой однородной жидкости, в котором —— = 0, —— = 0.

дг дг

Считаем известным распределение температуры и солености на п-м временном шаге:

г в?1 + 7П1 (г- - г1) при *1 < г < 32,

Тп(г) = в?2 + 7П2 (г- - £2) при ^2 < г < гз,

1 в?з при *з < г < Н;

г + 7П1 (*- - г1) при ^1 < г < г2,

Бп(г) = + 7П2 (г- - ¿2) при ^2 < г < *з,

1 в5з при *з < г < Н,

где = Лп — глубина распространения конвекции, Н - глубина водоема. В начальный

момент (начало образования ледяного покрова) г1 = 0, £ = 0, втх = Тф, = Б — соле-

ность на водной поверхности. Далее полагаем, что температура и соленость изменяются только в слое конвективного перемешивания.

Определим толщину слоя конвективного перемешивания (Л — £). Рассмотрим первый шаг по времени ¿1 = Д£. За время Д£ образовался ледяной покров толщиной £1. При замерзании воды высвобождается соль (лед практически “пресный”). В слое воды толщиной £1 солей содержалось:

51

J (Б + 7^г) ^ = Б£1 + 0.57^(С1)2.

о

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

ДБ1 = ая (БС1 + 0.575 (£1)2). (8)

Глубина конвекции на первом шаге по времени определяется из условия

Р1^1<г<Ь1 = Р^Н1,

где р|^1<*<^1 = р^р — плотность перемешанного слоя: р^ = р(Тс1р, Б^,).

Средние значения температуры и солености (Тр, Б^) в перемешанном слое находятся из уравнений

н1

Я 1 = 1 Г ад + ДБ1 = о/Л1 + £^ + ДБ1

Бср Л1 — £4 яаг + Л1 — С1 2 /+ Л1 — С1,

^ (9) т1 = ЛГ—Т / Т<к =

Учитывая (8) и (9), уравнение состояния (4) и предположение, что средняя плотность перемешанного слоя совпадает с плотностью подстилающего слоя воды, получаем формулу для определения Л1:

Л

1 = £ 1 + ,2взДБ1

а1Б0

^27т . ^з75 где Й1 = —----------+

То Бо

Аналогично определяется глубина распространения конвекции на (п+1)-м времен-

ном шаге:

Здесь

Й2

Лп+1 = £п+1 + ^/(£ п+1)2 — 2а2/аь

езДБ п+1 БО

— (Лп — £п+1)

Т2 (Тр — Тф) + Б^ (Бср — Б)

То Бо

+ а1 (Лп)2

ДБп+1 = а <?Бспр(£п+1 — £п)

тп+1 = ТСр (Лп — £п+1) + Тф (Лп+1 — Лп) + 0.57т ((Лп+1)2 — (Лп)2)

ср Лп+1 £п+1 :

Бп+1

ДБп+1 + Бср (Лп — £п+1) + Б (Лп+1 — Лп) + 0.57<? ((Лп+1)2 — (Лп)2)

ср Лп+1 £п+1

Таким образом, пока будет нарастать ледяной покров, будет увеличиваться и толщина слоя конвективного перемешивания.

2.3. Упрощенная модель таяния ледяного покрова весной

В весенний период температура льда повышается до температуры фазового перехода и тепловой поток из поверхностного слоя воды приводит к таянию льда снизу, которое рассчитывается по уравнению (7). Одновременно происходит таяние льда сверху. Для оценки этого процесса применим следующий подход (схема вертикальной структуры озера приведена на рис. 2).

Процесс таяния льда сверху описывается уравнением

рлД

-к

1Т„;

т

(10)

Здесь Тад (г) — температура воды, образовавшейся на поверхности льда, £ад — координата верхней поверхности льда.

Для оценки профиля температуры воды применим стационарное приближение:

к

12Т„:

а 1

т

1г2 срр0

—в(*—)

граничные условия:

К 1Т

кт -гаг

р

п

срр0

Тш |

Тф.

(11)

(12)

Рис. 2. Схема вертикальной структуры озера для периода таяния льда

Здесь Р/ — радиационный теплообмен, Рп — полный поток тепла через водную поверхность, определяемый по формуле (6), Кт=сопв1.

Решение задачи (11) и (12) имеет следующий вид:

Тад = Тф +

Срр0К

т

(Р — С)(£ад — г) + евНо (е-в5то — е-в")

(13)

где С = Рэф + Ре + Рс.

Для определения положения верхней границы ледяного покрова из (10), (13) получаем уравнение

рл^л% = —— + — (1 — е-в(5--н-)) . (14)

Срр0 Срро

Толщина образовавшегося на льду слоя воды вычисляется по формуле

п рл

Л п+1 = 7п + рл (£п+1 £п+1)

р у

рв

Для уравнения (14) задаются начальные условия:

£ I = £0 Л I = Л0

1^=0 1*=0 .

3. Примеры расчетов

С использованием подробных метеоданных были проведены расчеты термоклина в оз. Шира по одномерной модели. Константы, входящие в эмпирические соотношения, необходимо верифицировать для конкретного водоема. Для оз. Шира а = 0.6,

в = 0.004, = 0.9. Сравнение результатов расчетов с натурными данными показа-

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

По описанному вычислительному алгоритму была рассчитана динамика толщины льда пресного (оз. Иткуль) и соленых (озера Шира и Шунет) водоемов для конкретных метеоданных. По натурным данным толщина льда в глубоководной области оз. Шира 24.02.03 была равна 1 м, на оз. Шунет 28.02.03 — 0.73 м. Рассчитанная толщина льда на оз. Шира — 1.02 м, на оз. Шунет — 0.735 м. На рис. 3 приведено сравнение результатов расчетов динамики толщины льда для пресного и соленого озер. В соленом водоеме толщина льда меньше из-за наличия слоя конвективного перемешивания и вследствие более низкой температуры фазового перехода.

Ноябрь Декабрь Январь Февраль Март Апрель Май і

Iі / .У / 2

х ——

?, м

Рис. 3. Изменение толщины льда в соленом (1) и пресном (2) водоемах

Расчетная толщина льда зимой 2002 —2003 гг.

14 ноя 29 ноя 14 дек 29 дек 13 янв 28 янв 12 фев 27 фев 14 мар 29 мар 13 апр 28 апр

Дата

Рис. 4. Толщина льда и глубина конвективного слоя в глубоководной части оз. Шира, рассчитанные по метеоданным 2002-2003 гг. Точкой помечена измеренная толщина льда

Рис. 5. Профили температуры, солености и плотности для глубоководной части оз. Шира: ооо — натурные данные,----------------------------рассчитанные значения

Рис. 6. Профили температуры, солености и плотности для глубоководной части оз. Шира: ооо — натурные данные,----------------------------рассчитанные значения

Были рассчитаны толщина льда и глубина конвективного слоя (рис. 4), а также вертикальные распределения солености и температуры в оз. Шира для метеоданных 2002-2004 гг., взятых с сервера “Погода России” — архив погоды — Шира (http://meteo.infospace.ru/). На рис. 5 и 6 показаны рассчитанные и полученные в результате натурных измерений распределения по глубине температуры и солености в глубоководной части оз. Шира.

Авторы выражают благодарность Д.Ю. Рогозину за предоставленные данные натурных измерений.

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

[1] Белолипецкий П.В., Генова С.Н., Грицко В.В. Компьютерная модель вертикальной структуры водоема // Тр. Междунар. конф. “Вычислительные и информационные технологии в науке, технике и образовании”. Вестн. КазНУ им. аль-Фараби. Сер. “Математика, механика, информатика”. 2004. № 3(42). Ч. 1. С. 289-294.

[2] Кочергин В.П. Теория и методы расчета океанических течений. М.: Наука, 1978. 128 с.

[3] Belolipetskii V.M., Genova S.N. Investigation of hydrothermal and ice regimes in hydropower station bays // IJCFD. 1998. Vol. 10. P. 151-158.

[4] Математическое моделирование конвективного тепломассообмена на основе уравнений Навье—Стокса / В.И. Полежаев, А.В. Буне, Н.А. Верозуб и др. М.: Наука, 1987. 271 с.

[5] Судольский А.С. Динамические явления в водоемах. Л.: Гидрометеоиздат, 1991. 262 с.

[6] Белолипецкий В.М., Генова С.Н., Туговиков В.Б., Шокин Ю.И. Численное моделирование задач задач гидроледотермики водотоков. Новосибирск: Сиб. отд-ние РАН, ИВТ, ВЦК, 1994. 138 с.

[7] Красс М.С., Мерзликин В.Г. Радиационная теплофизика снега и льда. Л.: Гидрометеоиздат, 1990. 261 с.

[8] Гилл А. Динамика атмосферы и океана. М.: Мир, 1986. Т. 2. 415 с.

[9] Доронин Ю.П. Взаимодействие атмосферы и океана. Л.: Гидрометеоиздат, 1981. 288 с.

Поступила в редакцию 28 января 2008 г., в переработанном виде — 25 марта 2008 г.