Численный расчет продольно-вертикальной термической структуры Телецкого озера в годовом цикле

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

Том 5, № 3, 2000

ЧИСЛЕННЫЙ РАСЧЕТ ПРОДОЛЬНО-ВЕРТИКАЛЬНОЙ ТЕРМИЧЕСКОЙ СТРУКТУРЫ ТЕЛЕЦКОГО ОЗЕРА В ГОДОВОМ ЦИКЛЕ *

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

Барнаул, Россия e-mail: kvon@ad-sbras.nsc.ru

The lake Teletskoye ice-thermal regime was studied on the basis of longitudinal-vertical (lateral-averaged) mathematical model. The model takes into account water compressibility, meteorological factors, wind and inflow discharges currents. The latter are determined by the river inflows, including lateral and the lake outflow. Allowing for the water compressibility makes it possible to describe more accurately the hydrothermal conditions in the deep part of the lake. The additional equation for the ice heat balance was included into the model for the description of the thermal regime in winter. The numerical simulations of Lake Teletskoye annual cycle thermal structure were performed. The comparison of the calculated vertical temperature profiles for different moments of time and for number of places and also of ice cover thickness with field measurement data is given.

Введение

Первые работы по математическому моделированию ледотермического режима Телецко-го озера опираются на нестационарные одномерные вертикальные модели. В работе [4] моделируются вертикальные термический режим и турбулентная структура озера в периоды осенне-зимнего выхолаживания и весеннего прогревания, а также динамика толщины льда. В работе [6] изучается годовой термический режим Телецкого озера на основе вертикальной одномерной модели. Для описания гидротермического режима озера в ней использовалась методика осреднения уравнений баланса тепла, импульса и параметров турбулентности по горизонтальным сечениям водоема. При моделировании процессов роста и таяния ледового покрова учитывалось теплосодержание слоя льда.

Однако следует отметить, что из-за значительной протяженности в продольном направлении и относительно малой ширины водоема существенные изменения в сезонной динамике термической структуры озера происходят не только в вертикальном направлении, но и вдоль его продольной оси. Поэтому при моделировании этого озера в работах [9, 10] используется предложенный в [20] подход с осреднением по ширине озера уравнений

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №99-01-00547.

© Д. В. Квон, В. И. Квон, А.Н. Семчуков, 2000.

импульсов, сохранения массы, переноса тепла и уравнений для энергии турбулентности и скорости ее диссипации. В них проведено исследование термического режима Телецкого озера в весенне-летний период без учета сжимаемости воды. При этом предполагается, что плотность воды в уравнении состояния зависит только от температуры. Вместе с тем недавние теоретические исследования (см., например, [17, 8, 22]) глубоких озер, в особенности озера Байкал, показали, что сжимаемость воды, точнее, зависимость плотности воды от давления в уравнении состояния, играет важную роль в поведении их глубинных вод. Это главным образом связано с тем, что температура воды максимальной плотности убывает с возрастанием глубины (с увеличением давления воды) от 4 °С на водной поверхности на 0.2 °С каждые 100 м.

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

1. Постановка задачи

Уравнения осредненной по ширине озера модели гидротермических процессов [10] с учетом сжимаемости воды имеют следующий вид:

/ г0 \

дЬи дЬпп дЬши ,д 1 Г , \ д , ди д , ди . . ,„.

—¡тт + + = -ЗЬ^- П + — \Ф\ +ттЬКХ- + — ЬКг— - 8т1и1и, (1) дъ дх дг дх у ро 7 I дх дх дг дг

дЬТ дЬиТ + ^— +

дъ

дх

дЬи дЬ-ю дх дг

дЬ-Т д 1ГЛ дТ д 1ГЛ дТ

— = т- ЬОх— + — ЬПг— + дТк дг дх дх дг дг

р = рт/(1.0 - р/кр),

2

где ^ = ^

г=1

1 +

1/2

ъ

время, ось х направлена вдоль озера по его

дъу + (дЬи

дх) \дг

длине от р. Бия в сторону р. Чулышман, ось г — вертикально вверх; и и ю — компоненты скорости по х и г соответственно; Т — температура воды; д — распределение удельного расхода боковых притоков; Тп — температура воды боковых притоков; Ь(х, г) — ширина озера; Ь = Ь2 + Ь1; — Ь1 и Ь2 — ординаты боковых поверхностей озера; д — ускорение силы тяжести; Кх и Кг (Ох и Ог) — коэффициенты турбулентной вязкости (температуропроводности) соответственно в горизонтальном (х) и вертикальном (г) направлениях; П — отклонение уровня воды от его невозмущенного значения го; т — коэффициент сопротивления трения боковой поверхности озера; р — плотность воды и ро — характерное значение этой величины; рт (Т,Б) — плотность воды при стандартном атмосферном давлении, Б — минерализация воды, принимаемая здесь постоянной и равной 75 мг/л [16]; кр(Т, Б,р) — объемный модуль упругости воды; р — давление воды, определяемое по формуле гидростатики. Для определения функций рт(Т, Б) и кр(Т, Б,р) используются формулы, рекомендуемые Объединенной комиссией ЮНЕСКО по океанографическим таблицам и стандартам (см., например, [2]).

К системе уравнений (1)-(4) необходимо присоединить начальные и граничные условия. В начальный момент времени задаются распределения скорости, температуры и уровня водной поверхности. На водной поверхности при г = п задаются кинематическое условие, напряжение ветра и поток тепла:

дП дп „ ди Т™ п дТ

аь + идх = т Кдг = 7' ^дг =^

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

На дне водоема при г = г0(х) задаются условие непроницаемости, квадратичный закон сопротивления трения и равный нулю поток тепла:

^ ди , . . дТ т = 0, — = кьЫМ' -г— = 0.

дг дг

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

Напряжение ветра и поток тепла через водную поверхность в (5) вычислялись с использованием потоков солнечной радиации и среднесуточных метеоданных (скорость ветра, температура и влажность воздуха, атмосферное давление и облачность) по следующим формулам:

Т™ = Ра С™ Ж2,

где Си = 1.1  10-3 при Ж < 6 м/с и с™ = 1.0  10-3(0.72 + 0.063Ж) при W > 6 м/с, ра = 1.3 кг/м3 — плотность воздуха, Ж — скорость ветра (см., например, [1]);

Яп = Язг + Яаг - ЯЬг - Яе - Я'

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

Для вычисления этих потоков (ккал/(ч-м2)) используются формулы [19, 23]:

Язг = Язе (1 - «и) (1 - 0.65 С2)'

Яаг = 4.46 х 10-13 (Та + 273.15)6 (1 + 0.17 С2), Яы = 4.74 х 10-8 (Тз + 273.15)4'

Яе = /

Я = 0.459 /

где Язе — поток солнечной радиации (ккал/(ч-м2)), падающий при ясном небе на единицу площади горизонтальной поверхности за единицу времени; С — облачность в долях единицы; аи — альбедо воды; Та — температура (°С) воздуха на высоте 2 м; Тз — температура (°С) водной поверхности; Ж2 — скорость ветра (м/с) на высоте 2 м; е2 — парциальное

давление водяного пара (мм рт. ст.) на высоте 2 м. Давление насыщенного водяного пара (мм рт. ст.) при температуре водной поверхности е3 и ветровая функция f (Ш2) вычисляются по следующим формулам:

( 5278

в, = 25.4 ехр 17.62 -

Ts + 273.15

f

4.3 W2, если A0v < 0.0148 W3, 3.54 W2 + 3.09 (A0v)1/3, если A0v > 0.0148 W23,

A0v Tsv Tav,

Tsv = (Ts + 273.15)/(1 - 0.378es/pa),

Tav = (Ta + 273.15)/(1 - 0.378 ea/pa),

где Tsv — виртуальная температура (°C) тонкого слоя пара в контакте с водной поверхностью, Tav — виртуальная температура (°C) воздуха, pa — атмосферное давление (мм рт. ст.).

Распределение удельного расхода q боковых притоков по длине озера в уравнении неразрывности (2) задавалось в виде дельта-функции в точках втекания притоков. Данные измерений температуры воды боковых притоков использовались при задании Tin в уравнении переноса тепла (3).

Коэффициенты вертикального турбулентного обмена определялись с использованием уравнения для энергии турбулентности e и скорости ее диссипации е [15]:

dbue dbwe д 7 ^ де д 7 ^ de 1 . . .

"ТТ" + IT- = 7ГbKxe^- + JT-bKze- + b

дх Oz дх дх Oz dz

dbe дbue дbwe д де д де

-тгг + -z--+ —z— = т- bKX£— + — bKzsT-

Ot дх дz дх дх дz дz

-,-bKx^ + 7-bKZ£— + ci^b

Здесь Кхе (Кхе) и Кге (Кге) — коэффициенты диффузии энергии турбулентности (скорости ее диссипации) соответственно в горизонтальном (х) и вертикальном (г) направлениях; о1е = 1.44, о2е = 2.0[1.0 - 0.3ехр(-Ие^)], Ие^ = в2/(ие), с3е = 0.8, V — кинематическая вязкость воды.

Для уравнений модели турбулентности (9) и (10) здесь приняты граничные условия, подробно описанные в [9, 10].

Члены, описывающие генерацию турбулентности в этих уравнениях, при учете сжимаемости воды и поперечного ветрового напряжения имеют вид:

P=K

ди OZ

+

Ov Oz

OT

G = -датKfi [ — + Ya

где коэффициент термического расширения определяется по формуле

в = -

1 др POT

2

2

2

p

7а — адиабатический градиент температуры, равный 10 °С-м (значение для пресной озерной воды, по данным D. Farmer; см., например, [3]). Для параметризации влияния течения в поперечном направлении озера на генерацию турбулентности P, следуя идее Архипова и Солбакова [1], в выражение для генерации введен поперечный сдвиг скорости dv/dz, вычисляемый с использованием математической модели ветрового течения [7].

Коэффициенты горизонтального турбулентного обмена определяются с использованием формулы Ричардсона [14]. Конкретные их значения приведены в [10].

Модель формирования и развития ледового покрова в работе представлена в виде уравнения для толщины льда 0(t), выражающего баланс тепла в нем [21]:

- fи=- qw'

где pi — плотность льда, Lf — скрытая теплота плавления льда, qni — поток тепла, поступающий в ледовый покров через его верхнюю поверхность, qw — поток тепла, поступающий из ледового покрова в воду через нижнюю поверхность льда (через границу раздела лед — вода).

Поток тепла изо льда в воду определяется по формуле

qw = -kw

где kw — интегральный коэффициент теплопроводности воды, Tw — температура воды у водной поверхности, Tf — температура замерзания воды.

Поток тепла через верхнюю поверхность льда qni вычисляется по формуле (6) для qn, в которой температура водной поверхности заменена на температуру верхней поверхности льда, а альбедо воды — на альбедо льда, вычисленное по формуле [21]:

0.45' если Ta < 0 °С,

ai = <

0.25 + 0.16 exp-a07Ta, если Ta > 0 °С.

Кроме того, в случае наличия ледового покрова в формулах (7) и (8) необходимо ввести поправочные коэффициенты, а именно:

qe = kef (W2)(es - б2),

qc = ke • 0.459 • f (W2)(Ti - Ta),

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

Для нахождения температуры Ti используется уравнение баланса потоков тепла на верхней поверхности льда:

где qi вычисляется по формуле

qni - qi = 0'

Tf - T-

qi = -ki^^,

в которой ki

— коэффициент теплопроводности льда.

Алгебраическое уравнение (13) решается относительно Т итерационным методом Ньютона. Если полученное значение Т больше температуры замерзания Tf, что соответствует таянию льда и образованию воды на его поверхности, то полагаем Т = Tf. В данной постановке задачи, как и в [21], пренебрегаем толщиной воды на поверхности льда.

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

При наличии ледового покрова на водной поверхности (на границе раздела лед — вода) в условии (5) вместо потока тепла из воздуха в воду Яп задается поток тепла изо льда в воду Я™, а именно:

_ дТ

СррДг — = Яи

Кроме того, предполагается, что на водной поверхности отсутствует ветровое воздействие (ти = 0 в условии (5)). (В проведенных ниже расчетах на основе вычислительного эксперимента было принято, что указанные условия задаются при толщине льда > 2 мм). Заметим, что при наличии ледового покрова исчезает источник генерации турбулентности, обусловленный сдвиговым течением, порождаемым ветром. Из результатов работы [18] следует, что важным источником турбулентности при наличии ледового покрова являются сейшевые колебания в озере. В данной работе порождаемая такими колебаниями турбулентность параметризуется через коэффициент турбулентной вязкости, равный 10-4 м2/с [18].

2. Результаты расчетов

Модель реализована с использованием полунеявной конечно-разностной схемы [10]. При численном решении уравнений движения используется алгоритм, основанный на схеме расщепления по физическим процессам. На первом дробном шаге осуществляется перенос импульса адвекцией и диффузией, а на втором дробном шаге моделируется адаптация гидродинамических полей [11]. Уравнение для толщины льда (11) решается методом Рун-ге — Кутта.

Напомним, что Телецкое озеро — глубокий проточный водоем вытянутой формы (рис. 1). Длина озера 77.8 км, максимальные ширина 5.2 км, глубина 325 м. Основная часть притока в озеро поступает через р. Чулышман в южном его конце, и отток происходит в противоположном, северо-западном, его конце через р. Бию [16].

Расчеты термической структуры Телецкого озера и анализ их результатов проведены для 1968-1969 гидрологического года, данные измерений по которому наиболее полно представлены в литературе и в материалах наблюдений [12, 13, 16]. Измерения показывают, что к 20 мая 1968 г. в озере устанавливается состояние гомотермии [16]. Эта дата была принята за начало обсуждаемых здесь расчетов. В качестве начальных данных приняты состояние покоя и равномерное распределение начальной температуры, равное в соответствии с измерениями 2.3 градуса.

Для анализа влияния сжимаемости воды на гидротермический режим Телецкого озера в его сезонной динамике были также проведены расчеты без учета сжимаемости (в уравнении состояния воды р = 0, что соответствует атмосферному давлению) при идентичных входных данных. Результаты расчетов показывают (рис. 2, б), что пренебрежение сжимаемостью усиливает интенсивность перемешивания воды в глубинных зонах озера. В этих зонах, согласно численному решению, под изотермой 4 °С имеет место однородное

- Рейдовые вертикали у \

19 \

Чулышман (

Рис. 1. Карта-схема Телецкого озера (В. В. Селегей, Т. С. Селегей [16]).

распределение температуры со значением около 4 °С. Это связано с потерей устойчивости и возникновением проникающей турбулентной конвекции, механизм которой хорошо известен (см., например, [5]). Суть его заключается в появлении турбулентности при гидростатической неустойчивости. В модели этот механизм реализован через генерацию турбулентности за счет силы плавучести.

При учете сжимаемости воды (рис. 2, а) глубинные области озера меньше подвергаются конвективному перемешиванию, связанному с потерей гидростатической устойчивости. Хотя различие по температуре небольшое, процессы турбулентного обмена отличаются

р. Бия р. Чулышман

О 10000 20000 30000 40000 50000 60000 70000

Длина, м

р. Бия р. Чулышман

О 10000 20000 30000 40000 50000 60000 70000

Длина, м

Рис. 2. Изотермы, вычисленные на 20 августа 1968 г. с учетом (а) и без учета (б) сжимаемости воды (продольно-вертикальное сечение вдоль оси озера).

качественно, что может привести к изменениям в процессах теплообмена и газообмена (в том числе переноса кислорода) в глубинных зонах озера.

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

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

Результаты расчетов обнаруживают перемещение фронтов весенне-летнего термобара к центральной части озера с последующим их слиянием (рис. 3 и 4). По данным наблюдений, к 15 июля фронты смыкаются и озеро становится однородным в продольном направлении [16]. По результатам численных расчетов, смыкание фронтов термобара и наступление однородности озера по длине (по крайней мере в значительной части озера) имеет место примерно с 13 - 15 июля.

На рис. 5 представлены продольно-вертикальные распределения температуры (изотермы), характеризующие поведение термического бара в период осеннего охлаждения озера. Фронты осеннего термобара сначала возникают также в северо-западной и южной частях озера и затем перемещаются к его центру.

На рис. 6 приводится сопоставление вычисленных и измеренных профилей температуры на рейдовой вертикали в центральной части озера (вертикаль №26) для 20 июня, 20 июля, 20 августа, 20 сентября, 20 октября и 20 ноября. Согласно результатам расчетов, на 20 августа значения температуры воды в придонной зоне озера составляют 3.6 — 3.8 °С, а на 20 ноября устанавливается равномерный вертикальный профиль температуры со значением около 4.0 °С, что соответствует данным измерений [16]. На рис. 7 сопоставлены вычисленные и измеренные профили температуры на 20 августа для разных рейдовых вертикалей, расположенных вдоль озера по его длине (см. рис. 1).

На рис. 8 дано сопоставление данных измерений толщины льда в пункте измерения у пос. Яйлю [13] и результатов расчета. И наконец, на рис. 9 приведено распределение толщины льда вдоль озера в различные моменты времени.

Метеоданные были измерены на суше, на береговой станции Яйлю. Так как скорость ветра при переходе с суши на озеро обычно возрастает, в проведенных расчетах были приняты следующие значения отношения скоростей ветра над озером и над сушей: 1.18 в июне и июле, 1.36 в августе и 1.5 в сентябре [16]. Для остального периода года значение отношения скоростей ветра было принято равным 1.5. Это связано с тем, что по имеющимся среднемноголетним данным указанное отношение в октябре — декабре равно примерно 1.5 [16]. Кроме того, учитывался тот факт, что из-за наличия ледового покрова воздействие ветра на термический режим озера в зимних условиях играет малую роль.

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

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

Предварительные результаты работы были доложены на XVI Международной школе-семинаре по численным методам механики вязкой жидкости (Новосибирск, 13-18 сентября 1998 г.). Авторы выражают благодарность О.Ф. Васильеву за полезные обсуждения, а также В. В. Селегею, В. В. Кириллову и С. В. Вахненко за предоставление дополнительных исходных данных по Телецкому озеру.

Рис. 3. Динамика продвижения фронтов весенне-летнего термобара: распределение температуры воды в Телецком озере на 5 (а) и 8 (б) июля 1968 г. (продольно-вертикальное сечение вдоль оси озера).

р. Бия р. Чулышман

О 10000 20000 30000 40000 50000 60000 70000

Длина, м

р. Бия р. Чулышман

О 10000 20000 30000 40000 50000 60000 70000

Длина, м

Рис. 4. Динамика продвижения фронтов весенне-летнего термобара: распределение температуры воды в Телецком озере на 11 (а) и 15 (б) июля 1968 г. (продольно-вертикальное сечение вдоль оси озера).

р. Бия р. Чулышман

О 10000 20000 30000 40000 50000 60000 70000

Длина, м

р. Бия р. Чулышман

О 10000 20000 30000 40000 50000 60000 70000

Длина, м

Рис. 5. Динамика продвижения фронтов осеннего термобара: распределение температуры воды в Телецком озере на 20 (а) и 25 (б) ноября 1968 г.

Рис. 6. Сопоставление вычисленных (сплошные линии) и измеренных (квадратики) вертикальных профилей температуры в Телецком озере на 20.06.68, 20.07.68, 20.08.68, 20.09.68, 20.10.68 и

20.11.68 (а - е) для рейдовой вертикали №26.

Рис. 7. Сопоставление вычисленных и измеренных вертикальных профилей температуры в Телецком озере на 20.08.68 для рейдовых вертикалей №22, 24, 6, 26, 13 и 27 (а - е).

Рис. 8. Изменение толщины льда во времени в Телецком озере вблизи пос. Яйлю (вертикаль №6). Сплошная линия — вычисленные значения толщины льда, точки — данные наблюдений.

Рис. 9. Распределения вычисленных значений толщины льда по длине Телецкого озера (с севера на юг) на 20-е числа ноября (XI), декабря (XII), января (I), февраля (II), марта (III), апреля (IV), мая (V). Сплошные линии относятся к четным месяцам, пунктирные — к нечетным.

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

[1] Архипов Б. В., СолБАков В. В. Расчет термодинамического режима водоема по двумерной модели. Изв. АН. Физика атмосферы и океана, 30, №5, 1994, 671-685.

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

[3] Голубев В. А. Геотермия Байкала. Наука, Новосибирск. Сиб. отд-ние, 1982.

[4] Дмитриев Н. В. Моделирование термодинамики приповерхностного слоя водоема при образовании и таянии льда. Метеорология и гидрология, №7, 1995, 68-76.

[5] ЗилитинкЕвич С. С. Проникающая турбулентная конвекция. Валгус, Таллинн, 1989.

[6] Зиновьев А. Т. Моделирование ледотермического режима Телецкого озера. Материалы международного симпозиума "Гидрологические и экологические процессы в водоемах и их водосборных бассейнах": Тез. докл. Новосибирск, 199б.

[Т] Игнатова Г. Ш., Квон В. И. Одномерная модель сезонного термоклина в озерах. Водные ресурсы, №б, 19Т9, 118-12б.

[8] Квон В. И., Квон Д. В. Численный анализ механизма глубокого проникновения поверхностных вод в прибрежной зоне озера в период весенне-летнего термобара. Вычисл. технологии, 2, №5, 199Т, 4б-5б.

[9] Квон Д. В. Численное моделирование сезонных изменений температуры воды в Те-лецком озере. Там же, 1, №1, 199б, 48-5б.

[10] Квон Д. В., Квон В. И. Численный расчет стоковых и термогравитационных течений в Телецком озере. Метеорология и гидрология, №б, 1998, б8-Тб.

[11] Марчук Г. И. Численные алгоритмы решения уравнений прогноза погоды. Докл. АН СССР, 156, №2, 19б4.

[12] Материалы наблюдений на озерах и водохранилищах. Дополнение к Гидрологическому ежегоднику, б, вып. 0. Главное управление гидрометеорологической службы при СМ СССР, Западно-Сибирское управление гидрометеорологической службы, озерная станция Яйлю, 19б8.

[13] Материалы наблюдений на озерах и водохранилищах. Дополнение к Гидрологическому Ежегоднику, б, вып. 0. Главное управление гидрометеорологической службы при СМ СССР, Западно-Сибирское управление гидрометеорологической службы, озерная станция Яйлю, 19б9.

[14] Озмидов Р. В. Горизонтальная турбулентность и горизонтальный обмен в океане. Наука, М., 19б8.

[15] Роди В. Модели турбулентности окружающей среды. В "Методы расчета турбулентных течений". Мир, М., 1984, 227-322.

[16] Селегей В. В., Селегей Т. С. Телецкое озеро. Гидрометеоиздат, Л., 1978.

[17] Цветова е. А. Численная модель термобара в озере Байкал. Метеорология и гидрология, №9, 1997, 58-б8.

[18] Lars Bengtsson. Mixing in ice covered lakes. Abstracts of the first international lake Ladoga simposium Î993. Joensuu, 1995.

[19] Ryan P. J., Harleman D. R. F, Stolzenbach K. D. Surface heat loss from cooling ponds. Water Resources Research, 10, No. 5, 1974.

[20] Vasiliev O.F., KvON V. I., Chernyshova R. T. Mathematical modelling of the thermal pollution of a water body. Froc. of the XV JAHR Congress, Istanbul, 2, 1973, 129-137.

[21] WAKE A., RuMER R. R. Modeling ice regime of lake Erie. J. Hyd. Div. ASCE, NHY7, No. 105, 1979, 827-844.

[22] Walker S. J., Watts R.G. A three-dimensional numerical model of deep ventilation in temperate lakes. J. Geophysical Research, 100, No. C11, 1995, 22,711-22,731.

[23] Wunderlich W. O. Heat and mass transfer between a water surface and the atmosphere. Laboratory Report TVA Engineering laboratory No. 14. Norris, 1972.

Поступила в редакцию 5 мая 1999 г., в переработанном виде 21 июня 1999 г.