Московский Государственный Университет

им. М.В. Ломоносова

Географический факультет

Кафедра метеорологии и климатологии

Институт вычислительной математики РАН

Дипломная работа

Численное моделирование процессов

тепловлагопереноса в системе «водоём – грунт»

Степаненко В.М. (5 курс)

Научный руководитель:

чл.-корр. РАН, д.ф.-м.н. Лыкосов В.Н.

Рецензент:

доцент, к.г.н. Архипкин В.С.

Москва – 2004 г.

Содержание Введение ……………………………………………………………………. 3 1. Физическая постановка задачи и описание модели…………….. 5 1.1. Тепловлагоперенос в теле водоёма……………………………….. 1.1.1. «Эмпирическая» параметризация………………………………… 1.1.2. Параметризация на основе уравнения для кинетической энергии турбулентности с использованием пути смешения (« E l » параметризация)………………………………………….. 1.1.3. Параметризация на основе уравнений для кинетической энергии турбулентности и скорости её диссипации (« E » параметризация)………………………………………….. 1.2. Тепловлагоперенос в слое снега…………………………………… 1.3. Теплоперенос в слое грунта………………………………………… 1.4. Тепловой баланс на подстилающей поверхности………………… 2. Численная реализация модели

2.1. Аппроксимация уравнения теплопроводности……………………. 2.2. Аппроксимация и решение уравнения теплового баланса……….. 2.3. Аппроксимация уравнений для турбулентной энергии и скорости её диссипации…………………………………………….. 3. Верификация модели……………………………………………… 3.1. Эксперименты с данными наблюдений в Колпашево……………. 3.2. Эксперименты с данными наблюдений на озере Сырдах……….. 3.2.1. Моделирование сезонного промерзания и испарения…………… 3.2.2. Моделирование талика…………………………………………….. 3.2.3. Моделирование теплового режима водных масс………………… 3.2.4. Характеристики турбулентности………………………………….. Заключение…………………………………………………………………. Список литературы………………………………………………………… Введение В связи с прогрессом, достигнутым к настоящему времени в конструировании вычислительной техники и разработке систем параллельного программирования, современный этап развития математических моделей климатической системы характеризуется постоянным совершенствованием пространственного разрешения и отказом (пока на региональном уровне) от гидростатического приближения.

Указанные тенденции порождают новые проблемы в детализации описания физических процессов подсеточных масштабов, среди которых важное место занимает взаимодействие атмосферы с различными типами подстилающей поверхности на суше. Одним из ключевых вопросов здесь является описание процессов взаимодействия атмосферы с сетью мелких гидрологических объектов, важнейшую часть которой составляют небольшие озёра и болота.

Это особенно важно для северных территорий Евразии (Западно-Сибирская низменность, Карелия, Финляндия) и Северной Америки (большая часть территории Канады), где данная сеть представлена наиболее плотно и где, как показывают эксперименты с климатическими моделями, региональные температурные изменения в связи с глобальным потеплением проявляются наиболее ярко. Для адекватной параметризации процессов взаимодействия атмосферы и суши в этих условиях необходимо, чтобы соответствующий блок климатической модели учитывал эффекты "гидрологической неоднородности" подстилающей поверхности. Особое значение при этом приобретают исследования, направленные на сравнение различных подходов к решению этой задачи на основе использования рядов длительных метеорологических наблюдений (приземные температура, влажность, давление, скорость ветра и т.п.) в отдельных регионах в качестве входных параметров. Примером таких исследований является проект сравнения схем параметризации происходящих на поверхности суши процессов PILPS-2(e) с особым акцентом на гидрологические процессы высоких широт (Bowling et al., 2003 [13]).

Значительный интерес представляет моделирование процессов тепловлагопереноса в системе "водоём - почва", где под водоёмом понимается мелкое озеро или болото. До настоящего времени представление болот в климатических моделях осуществляется как соответствующая спецификация того или иного участка подстилающей поверхности без учета происходящих в их толще термодинамических процессов, роль которых в процессах обмена энергией и массой между атмосферой и сушей все еще недостаточно изучена. Вместе с тем к настоящему времени накоплен большой опыт по решению задач, связанных с отдельными аспектами взаимодействия атмосферы и болот (см., например, работы по расчету водного режима болот при их осушении (Балясова, Пакутин, 1988 [1]), по гидравлике болот (Кашеваров, 2002 [4]), по исследованию процессов генерации и переноса метана (Крылова, Крупчатников, 2002 [5])).

Озёра также существенно воздействуют на структуру приземного слоя атмосферы и тем самым на потоки тепла, водяного пара и импульса. В большинстве моделей прогноза погоды и в климатических моделях эффекты, связанные со сравнительно небольшими и мелкими озёрами, либо совсем не учитываются, либо параметризуются очень грубо, например, в предположении, что водоём хорошо перемешан по глубине. Это, фактически, означает, что озеро рассматривается как элемент подстилающей поверхности. В реальности же озера в умеренных и высоких широтах большую часть года вертикально стратифицированы по плотности (Mironov et al., 1991 [17]). Вместе с тем, описание эффектов вертикальной стратификации на основе современных теорий турбулентного переноса (см., например, (Tsuang et al., 2001 [22])) всё ещё требует значительных вычислительных затрат. Это в особенности существенное ограничение, если принять также во внимание необходимость рассмотрения на больших временных масштабах процессов тепловлагопереноса в расположенном под водоёмом слое грунта.

В данной работе представлен компромиссный подход к разработке параметризации эффектов болотно-озёрных объектов, сочетающий достаточную полноту физического описания процессов тепловлагопереноса в системе "водоём-грунт" и вычислительную эффективность соответствующих алгоритмов реализации. Работа построена следующим образом. Вначале (раздел 1) формулируется физическая постановка задачи и приводится описание модели. В разделе 2 излагается численная реализация модели.

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

1. Физическая постановка задачи и описание модели Существенной особенностью мелкого водоёма (см., например, экспериментальные работы (Malm et al., 1997 [16], Pavlov, 1995 [18])) является сравнительно небольшая горизонтальная неоднородность в распределении его термодинамических параметров. Это позволяет в первом приближении рассматривать только процессы вертикального переноса тепла в теле водоёма. Очевидно, что на рассматриваемых временных масштабах одномерное приближение применимо и к расположенному под водоёмом грунту, а также к формирующимся в зимний период снежному и ледяному покрову.

термодинамического режима водоёма положено одномерное уравнение теплопроводности. Предполагается, что начало направленной вниз вертикальной координаты z совмещено с уровнем свободной поверхности водоёма, так что область, в которой ищется решение, представляет собой отрезок [0, h], где h=h(t) – толщина водоёма, а t - время. Удобно перейти от исходной вертикальной координаты z к новой независимой переменной =. В переменных (, t ) уравнение теплопроводности принимает вид:

(1) Здесь с – теплоемкость воды, – её плотность, – коэффициент теплопроводности, Т – температура, - скорость приращения слоя воды на верхней границе, с которой связано начало системы координат, S – поток солнечной радиации в толще водоёма, M – слагаемое, учитывающее влияние вертикального плотностного перемешивания на температуру. Для расчета потока солнечной радиации в толще водоёма применяется широко используемая в различных исследованиях экспоненциальная зависимость где е – коэффициент экстинкции.

вертикального обмена тепла и массой в водоёме. Его параметризация в предлагаемой модели осуществлена в трех вариантах.

1.1.1. «Эмпирическая» параметризация. Входящие в неё параметры определялись в ходе численных экспериментов и подобраны так, чтобы было наилучшее согласие модельных и экспериментальных данных. При этом параметризуется процесс вертикального перемешивания в случае неустойчивой плотностной стратификации (слагаемое М), а при вычислении коэффициента теплопроводности учитывается ветровое воздействие на турбулентный режим водоёма.

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

причём температура T выражена в градусах Цельсия.

Учет ветрового воздействия при вычислении выражен следующей зависимостью:

где s - коэффициент теплопроводности на поверхности водоема, b m, m - коэффициент молекулярной некоторый коэффициент, теплопроводности, V – скорость ветра, V0 – скорость ветра, при которой коэффициент теплопроводности достигает максимального значения max.

Очевидно, что описанная параметризация не лишена недостатков:

отсутствует, в частности, теоретическое обоснование задания временного интервала мгновенного перемешивания, а также выбора параметров формулы для расчёта. Другими словами, эти параметры могут сильно различаться в разных географических регионах, а также между отдельными водоемами. Главное же достоинство рассмотренной параметризации – простота её реализации. В литературе (см., например, (Simoes 1998[21])) турбулентной теплопроводности, однако и они обладают перечисленными выше недостатками, и потому в данной работе не используются.

1.1.2. Параметризация на основе уравнения для кинетической энергии турбулентности с использованием пути смешения («Е - l»

параметризация). В рамках данного подхода используется формула:

где l - путь смешения (масштаб турбулентности), E = (u ')2 + (v ')2 + ( w ') 2 кинетическая энергия турбулентности (черта сверху – знак осреднения, величины со штрихами – отклонения от среднего значения). Кинетическая энергия турбулентности рассчитывается с помощью следующего уравнения (Монин, Яглом, 1965 [10]):

где k = l E - коэффициент турбулентности, E - безразмерная константа, турбулентности за счет сдвига скорости и за счет эффекта плотностной стратификации:

а - скорость диссипации кинетической энергии турбулентности.

Для вычисления зональной u и меридиональной v компонент скорости течения воды в общем случае предполагается использовать уравнения Экмана, в которых учтен турбулентный перенос импульса и отклоняющая сила вращения Земли:

где f - параметр Кориолиса. Однако в текущей версии модели сила Кориолиса не учитывается, что может быть обосновано сравнительно небольшими размерами рассматриваемых озёр. Скорость диссипации энергии турбулентности рассчитывается по формуле Колмогорова:

где Ce = 0.07.

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

1.1.3. Параметризация на основе уравнений для кинетической параметризация). В рамках данной получившей широкое признание и распространение (Shnaydman, 2002 [20]) параметризации коэффициент турбулентной теплопроводности определяется следующим образом:

где k - коэффициент турбулентности, который вычисляется по следующей формуле:

Для расчёта турбулентной энергии E и скорости ее диссипации используются следующие уравнения (Лыкосов, 1992 [7]):

где E и - безразмерные константы, C1 - функция числа Рейнольдса Re :

- молекулярная вязкость воды, а константа C0 принята равной 1.9.

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

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

Граничные условия для энергии турбулентности и скорости её диссипации записываются в следующем виде. На свободной поверхности = 0 (граница вода-атмосфера) предполагается, что где C - безразмерная константа (в текущей версии модели её значение принято равным 1), а a - модуль напряжения трения ветра. На твёрдых границах (вода-лед, вода-грунт) аналогичным образом принимается, что Граничные значения для определяются из граничных значений E по формуле (3), где путь смешения l принят равным 0.4h.

На границе раздела «водоём-атмосфера» задается условие теплового баланса (см. ниже), а на границе раздела «водоём-грунт» используются условия непрерывности температуры и потока тепла. Для определения h записывается уравнение водного баланса водоёма:

где r – интенсивность осадков, Es – скорость испарения, Rs – поверхностный сток, Rb – водообмен тела водоёма с нижележащим грунтом. Горизонтальный грунтовый сток в модели не рассматривается.

поверхности водоёма опускается ниже 0 °С, то образуется первичный слой льда толщиной 1 см. В этом слое также решается уравнение (1), но при этом полагается, что суммарная солнечная радиация (за вычетом отраженной) полностью поглощается на поверхности льда и вниз не проникает. Это означает, что в уравнении (1) не учитывается четвертый («радиационный») член в правой части, а также не используется «конвективное» слагаемое M.

На границе раздела лёд-вода в качестве граничного условия задаётся температура фазового перехода 0 °С. Таяние льда на границе с атмосферой происходит при температуре фазового перехода и его скорость определяется тепловым балансом на этой границе. При этом считается, что стаявшая вода мгновенно добавляется к водному слою.

1.2. Тепловлагоперенос в слое снега. В зимнее время на замерзший водоём могут выпадать твердые осадки, образуя снежный покров, который в модели характеризуется распределением двух основных параметров – температуры и удельного содержания жидкой влаги. Их эволюция рассматривается в координатах (z, t) и описывается следующей системой уравнений (Володина и др., 2000 [3]):

Здесь L – удельная теплота плавления, Ffr – скорость замерзания, W – удельное содержание жидкой влаги, – инфильтрационный поток жидкой влаги в снежном покрове. Поток в разностном выражении выглядит так:

где h – гидравлическая проводимость, м/с (в модели h=0.01 м/с), П – пористость снега, Whc – константа, характеризующая водоудерживающую способность снега (в модели Whc=0.04), z – шаг сетки модели по вертикальной координате. Кроме процессов, описываемых системой (7), в модели также используется параметризация процесса гравитационного оседания (уплотнения) снежного покрова во времени. На нижней границе снежного покрова (поверхность раздела «снег-лёд») температура и поток тепла предполагаются непрерывными. Для определения температуры поверхности снежного покрова также используется уравнение теплового баланса.

1.3. Тепловлагоперенос в слое грунта. В основу описания процессов представленная в работе Володина и Лыкосова (1998)[2], в которой состояние почвы характеризуется температурой и содержанием жидкой, твердой и газообразной влаги. Поскольку под телом водоёма грунт должен быть насыщенным жидкой (при промерзании - твердой) влагой, то содержанием водяного пара можно пренебречь. Если, кроме того, пренебречь соответствующая система уравнений принимает вид:

Это часто принимаемое допущение, основанное на малости коэффициента термовлагопроводности.

Здесь W – коэффициент влагопроводности, а I - удельное содержание льда.

Как видно из этой системы, в почве рассматриваются процессы диффузии тепла и влаги, инфильтрация жидкой влаги, а также промерзание/таяние воды. Коэффициенты, определяющие интенсивность этих процессов, зависят от переменных состояния грунта – T, W, I.

Гравитационный поток влаги в почве определяется по формуле, приведённой в работе (Clapp and Hornberger, 1978 [14]):

где soil,max - максимальное значение гравитационного потока2, Wmax максимально возможное содержание жидкой влаги (оно достигается при заполнении ею всего объема пор, свободного ото льда), bsoil – безразмерный показатель. Суммарная теплоёмкость почвы складывается из теплоёмкостей сухой почвы, воды и льда:

а коэффициенты диффузии влаги и тепла вычисляются в соответствии с формулами:

м):

Нижним индексом max мы будем обозначать максимальное значение соответствующей величины, а индексом soil – величину, зависящее от типа почвы.

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

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

где S – суммарный поток солнечной радиации, Fa – поток уходящего длинноволнового излучения атмосферы, Fs – поток тепла, связанный с собственным излучением поверхности, Hs и LEs – потоки явного и скрытого тепла, соответственно, – альбедо поверхности. В переходные сезоны, когда осадки в твердом виде могут выпадать на незамерзший водоем, или случается дождь в период ледостава, в уравнение (4) добавляются слагаемые, учитывающие энергетический вклад этих процессов. Все участвующие в (4) метеорологических станциях.

Для расчёта потока суммарной солнечной радиации используется формула Кондратьева (Матвеев, 2000 [9]):

где 0 – эмпирическая функция высоты солнца, – оптическая толщина атмосферы для интегрального потока, принимаемая равной 0.105, n – балл облачности в долях единицы, csh = 0.5607 – эмпирический коэффициент, h высота солнца. Приход солнечной радиации на горизонтальную площадку на верхней границе атмосферы определяется следующей формулой:

где S 0 - солнечная постоянная, – широта, – склонение солнца, – часовой угол.

Уходящее длинноволновое излучение атмосферы рассчитывается как функция температуры и влажности на высоте 2 м и балла облачности (Матвеев, 2000 [9]):

В этих выражениях T2 - температура воздуха на уровне 2 м, e2 - парциальное давление водяного пара на уровне 2 м, - постоянная Стефана - Больцмана.

Собственное излучение поверхности задается известной формулой Стефана – Больцмана.

аэродинамический метод с коэффициентами, рассчитываемыми согласно Соответствующие формулы имеют вид:

где cp – теплоёмкость воздуха при постоянном давлении, а – плотность воздуха, CH, CE – коэффициенты обмена соответственно для температуры и влажности воздуха, 2 и q2 – потенциальная температура и удельная влажность на уровне 2 м, s и qs – те же величины на поверхности земли, V – скорость ветра на уровне 2 м.

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

разностный аналог уравнения (1) при условии отсутствия ледостава и в теплопроводности имеет следующий вид:

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

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

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

2.2. Аппроксимация и решение уравнения теплового баланса. Для вычисления температуры на границе с атмосферой используется уравнение теплового баланса (9), которое в разностном виде имеет вид:

величину на поверхности, ( s + 1) - величину, относящуюся к первому узлу сетки под поверхностью, а – на уровне 2 м в атмосфере, звездочками * указаны величины, которые берутся из данных натурных наблюдений, т.е.

известны. Решение уравнения (13) Ts j +1 ищется приближенно с помощью итерационной процедуры. Этот итерационный процесс содержит три уровня итераций.

Начальным приближением является Bsj ; при заданной правой части (13) из левой части находится Ts j +1, затем из (11) прогонкой рассчитывается новый профиль температуры, и по нему новое приближение Bsj +1, далее процесс повторяется снова, и т.д. Назовем этот процесс внешним. Следующий, заданном Bsj +1. При этом для гарантии сходимости процесса итераций в (13) Наконец, процесс внутренних итераций используется для определения Hs и LEs в (13) при заданном Ts j +1 по формулам аэродинамического метода с коэффициентами по теории подобия Монина – Обухова (10). Для получения нового приближения по внешнему циклу итераций необходимо осуществить промежуточного и внутреннего циклов.

Промежуточные и внутренние итерации осуществляются согласно методу последовательных приближений (Лебедев, 2000 [6]). Внешние итерации производятся чебышевским одношаговым циклическим методом;

соответствующие формулы имеют вид:

Здесь k - коэффициенты метода, k – номер итерации, N – период баланса:

а Ts ++ находится последующей прогонкой из (11).

2.3. Аппроксимация уравнений для турбулентной энергии и скорости её диссипации. В модели используются следующие разностные аналоги исходных уравнений (6):

Как видно, данная схема является полностью неявной относительно искомых величин E и. Это обеспечивает ее устойчивость. В то же время, в силу её нелинейности, решение приходится отыскивать с помощью итерационного процесса. На каждом его шаге решение первых двух разностных уравнений определяется методом прогонки. Скорости uij +1 и vij +1 к моменту решения уравнений (15) уже определены из следующих уравнений:

После определения Ei j +1 и i j +1 по формуле (5) и (4) находится коэффициент турбулентной теплопроводности ij +1.

С целью верификации модели проведены численные эксперименты с использованием рядов стандартных метеорологических измерений. Данные наблюдений на 225 метеорологических станциях бывшего СССР. В настоящей работе приведены некоторые результаты экспериментов для Колпашево (Томская область, правобережье р. Обь) и Якутска. В последнем случае результаты моделирования сравнивались с данными натурных наблюдений на оз. Сырдах (20 км к северу от Якутска, левобережье р. Лена).

Выбор оз. Сырдах для численных экспериментов объясняется тем, что на этом водоёме проводились регулярные наблюдения за водно-тепловым режимом, включая теплобалансовые измерения.

3.1. Эксперименты с данными наблюдений в Колпашево. В экспериментах с данными по станции г. Колпашево за период 1936 – 1984 гг.

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

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

Точность измерения температуры поверхности снега на ст. Колпашево составляет 0.5 °С. Как видно, рассчитанные и наблюденные кривые очень хорошо коррелируют. В то же время имеет место систематическое занижение моделью среднемесячного значения рассматриваемого параметра (на ~2 °С).

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

Температура, С Температура, С Рис. 1. Температура поверхности снега в Колпашево по данным моделирования и измерений 3.2. Эксперименты с данными наблюдений на озере Сырдах. Озеро Сырдах – сравнительно большое, находящееся в цепи таких же, по размерам и глубине, озёр, соединяющихся только в периоды редких для этих территорий многоводий ручьями типа «травяных речек». Оно вытянуто в направлении с СЗ на ЮВ и занимает большую часть площади котловины.

Размеры озера составляют: в длину - 2 км, в ширину до 1 км, площадь зеркала воды около 2 км2. Средняя глубина воды в озере составляет величину 4.5 м, максимальная достигает 12 м. Ледостав начинается в первой половине октября, сход ледяного покрова приходится на конец мая. Озеро расположено над мерзлыми породами, мощность которых составляет 280 – 320 м, и под ним существует сквозной талик (Павлов, Тишин, 1981 [11]).

нижележащих грунтов производилось за период 1965 – 1984 гг., для которого были доступны данные метеорологических наблюдений в г. Якутск и данные натурных измерений (Павлов, Тишин, 1981 [11]) на самом озере (1976 – гг.).

Эволюция водного, ледяного и снежного покровов в данном озере по данным моделирования представлена на рис. 2. По результатам качественного сравнения модельных и натурных данных (Павлов, Тишин, 1981 [11], Pavlov, 1995 [18]) можно сделать следующие выводы:

многоснежные сезоны, что соответствует реально наблюдаемым ситуациям;

- характерные глубины промерзания (толщина слоя льда) в модельном водоеме находятся в пределах 0.7 – 1.5 м, что совпадает с реальными значениями для озёр Центральной Якутии;

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

составляет в среднем 400 мм, что близко к наблюденному значению 450 мм.

3.2.2. Моделирование талика. На рис. 3 в виде распределения термоизоплет представлены результаты моделирования теплового режима талика и мерзлого грунта под озером (нулевая термоизоплета нанесена жирным пунктиром). Как видно из рисунка, талик устойчиво существует под озером в течение 20 лет (1965 – 1984 гг.). Его нижняя граница колеблется от 1.2 до 2.0 м под дном озера. В то же время, согласно данным измерениям (Тишин, 1979 [12]) глубина талика под оз. Сырдах достигает десятков метров, а в центральной части озера, по всей видимости, превышает 100 м. Такое несоответствие модельных и эмпирических данных объясняется тем, что образование глубокого талика - процесс, существенно превышающий по длительности период интегрирования модели (20 лет), так что нижняя граница модельного талика просто «не успевает» достигнуть глубины наблюдаемой границы.

Рис. 2. Промерзание/таяние и снегонакопление в оз. Сырдах по данным Глубина, моделирования 3.2.3. Моделирование теплового режима водных масс. По данным моделирования были получены среднемесячные вертикальные профили температуры в оз. Сырдах за период сентябрь 1976 – август 1977 гг. и сравнены с данными измерений (Павлов, Тишин, 1981 [11]). При этом в E параметризации вертикального турбулентного обмена. Эти результаты представлены на рис. 4 – 15.

Как видно, наблюдается хорошее согласование обоих модельных профилей с измеренным профилем в период ледостава (декабрь-апрель: рис.

7 – 11). Это естественно, поскольку в этот период физика теплопереноса в теле водоёма предельно проста: турбулентность находится в режиме коэффициент теплопроводности приближается по своему значению к молекулярному. Такая ситуация успешно моделируется, как показывают Примечательно, что вблизи границ с грунтом и слоем льда, кривая свидетельствует о повышенной интенсивности перемешивания в этих частях профиля (рис. 6 – 8). Это вызвано генерацией энергии турбулентности за счёт сдвига скорости вблизи границ. По мере полного затухания течений к концу зимы этот эффект также исчезает (рис. 9 – 11).

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

-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8. Глубина, м Рис. 5. Вертикальный профиль температуры в оз. Сырдах, Глубина, м Рис. 6.Вертикальный профиль температуры в оз. Сырдах, Глубина, м Рис 7 Вертикальный профиль температуры в оз Сырдах Глубина, м Рис. 8.Вертикальный профиль температуры в оз. Сырдах, Глубина, м Рис. 9.Вертикальное распределение температуры в оз. Сырда Глубина, м Рис 10 Вертикальный профиль температуры в оз Сырдах воспроизводится моделью. Водная масса оказывается более прогретой, чем это предсказывается моделью. Дело в том, что в период схода ледостава (май для оз. Сырдах) солнечная радиация начинает проникать под лёд и прогревать верхние слои воды, вызывая конвекцию. Именно этим объясняется хорошая перемешанность водной массы, которая видна на рис.

12. В модели проникновение радиации под лед не описано, поэтому указанные качества профиля не воспроизводятся.

Обращает на себя внимание тот факт, что в августе (рис. 15), сентябре и октябре (рис. 4 – 5) профиль температуры, рассчитанный с использованием E параметризации, существенно отклоняется в сторону увеличения от измеренного по всей глубине озера. Это связано, по всей видимости, с тем, что E - параметризация занижает интенсивность сезонного (а, возможно, выхолаживанию нижних слоев. При интенсивной конвекции в водоёме Глубина, м Рис. 12.Вертикальный профиль температуры в оз. Сырдах Глубина, м Рис. 13.Вертикальный профиль температуры в оз. Сырдах, Глубина, м Рис. 14. Вертикальный профиль температуры в оз. Сырдах Глубина, м Рис 15 Вертикальный профиль температуры в оз Сырдах развиваются конвективные ячейки, формирующие «крупномасштабную»

турбулентность и связанный с ней противоградиентный перенос (Lykossov, 2001 [15]). Именно вклад этих вихрей не учитывается в данной параметризации, поскольку моделируются эффекты лишь мелкомасштабной турбулентности.

Характеристики турбулентности. В ходе интегрирования модели в фиксированный срок (22 ч.) в разные дни июля 1977 г.

записывались вертикальные распределения турбулентных характеристик E,, k. Пример такого распределения показан в табл. 1. В ней также приведены типичные значения этих характеристик для водоемов, согласно работе (Wuest, Lorke, 2003 [23]).

Таблица 1. Распределение характеристик турбулентности по глубине в оз. Сырдах в 22:00 31 июля 1977 г. по данным моделирования.

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

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

Анализ результатов численных экспериментов с натурными данными для оз. Сырдах в Центральной Якутии показал, что модель адекватно воспроизводит следующие основные параметры климатического режима озера: среднюю глубину зимнего промерзания, время начала и окончания ледостава, среднее испарение в теплый период. Кроме того, в численном эксперименте воспроизведен талик под данным озером, существование демонстрирует хорошее согласие с экспериментальными данными при воспроизведении теплового режима водных масс данного озера. Сравнение рассчитанных и натурных рядов температуры поверхности снежного покрова для г. Колпашево (Томская область) показывает их хорошее соответствие.

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

руководителю Лыкосову В. Н. Работа с ним – непростая, но очень содержательная научная школа. Главное, что получает автор в общении с Василием Николаевичем – это постоянный импульс к движению вперед.

Кроме того, автор искренне признателен Мачульской Е. Е. за полезные консультации в области моделирования и программирования, которыми она помогала ему в процессе создания данной работы, Н.Ф. Вельтищеву, А.В.

Павлову и М.А. Петросянцу - за ценные замечания и предложения, Н.Г.

Яковлеву - за полезные консультации при разработке модели. Работа поддержана Российским фондом фундаментальных исследований (гранты 01и Международной ассоциацией содействия сотрудничеству с учеными из стран СНГ (гранты INTAS-00-189 и INTAS-01Список литературы 1. Балясова Е.Л., Пакутин А.В. Расчет изменений максимального весеннего стока с болотных массивов под влиянием осушительных мелиораций. - Тр. ГГИ, Вып. 333. Л.: Гидрометеоиздат, 1988, 152 с.

2. Володин Е.М., Лыкосов В.Н. Параметризация процессов тепло- и влагообмена в системе растительность-почва для моделирования общей циркуляции атмосферы. 1. Описание и расчеты с использованием локальных данных наблюдений. - Известия РАН. Физика атмосферы и океана, 1998, том 34, с. 453-465.

3. Володина Е.Е., Бенгтссон Л., Лыкосов В.Н. Параметризация процессов тепловлагопереноса в снежном покрове для моделирования сезонных вариаций гидрологического цикла суши. - Метеорология и гидрология, 2000, № 5, с. 5-14.

4. Кашеваров А.А. Математическое моделирование взаимодействующих течений подземных и поверхностных вод на заболоченных территориях.

Сб.: Большое Васюганское болото. Современное состояние и процессы развития. – Томск: Изд-во Института оптики атмосферы СО РАН, 2002, с.

83 – 87.

5. Крылова А.И., Крупчатников В.Н. Глобальное моделирование потоков метана от болотных экосистем. Сб.: Большое Васюганское болото.

Современное состояние и процессы развития. – Томск: Изд-во Института оптики атмосферы СО РАН, 2002, с. 98 – 103.

6. Лебедев В.И. Функциональный анализ и вычислительная математика.

М.: ФИЗМАТЛИТ, 2000.

7. Лыкосов В.Н. О проблеме замыкания моделей турбулентного пограничного слоя с помощью уравнений для кинетической энергии турбулентности и скорости её диссипации. - Изв. АН СССР. Физ. атм. и океана, 1992, т. 28, 696-704.

8. Монин А.С., Обухов А.М. Основные закономерности турбулентного перемешивания в приземном слое атмосферы// Тр. Геофиз. ин-та АН СССР. – 1954. - №24(151). – С. 163-187.

9. Матвеев Л. Т. Физика атмосферы. С.-П.: Гидрометеоиздат, 2000, 778 с.

10. Монин А.С., Яглом А.М. Статистическая гидромеханика (механика турбулентности). Ч. 1. М.: «Наука», 1965, 639 с.

11. Павлов А.В., Тишин М.И. Тепловой баланс крупного озера и прилегающей территории в Центральной Якутии. - В кн.: Строение и тепловой режим мерзлых пород; Новосибирск, «Наука», 1981.

12. Тишин М.И. Температурный режим горных пород под крупным термокарстовым озером в Центральной Якутии. - В кн.: Региональные и криолитологические исследования в Сибири. Якутск, 1979.

13. Bowling L.C., Lettenmaier D.P., Nijssen B., Graham L.P., Clark D.B., Maayar M.E., Essery R., Goers S., Gusev Y.M., Habets F., van den Hurk B., Jin J., Kahan D., Lohmann D., Ma X., Mahanama S., Mocko D., Nasonova O., Niu G.-Y., Samuelsson P., Shmakin A.B., Takata K., Verseghy Y., Viterbo P., Xia Y., Xue Y., Yang Z.-L. Simulation of high-latitude hydrological processes in the Torne-Kalix basin: PILPS Phase 2(e). 1: Experiment description and summary intercomparisons. - Global and Planetary Change, 2003, vol. 38, pp.

1-30.

14. Clapp R. B., Hornberger M. G. Empirical equations for some soil hydraulic properties. -Water Resources Research. 1978. V. 14, N 4. P. 601 – 604.

15. Lykossov V. N. Atmospheric and oceanic boundary layer physics. - Wind stress over the Ocean (Eds. Ian S.F. Jones and Yoshiaki Toba), Cambridge University Press, 2001, pp. 54 – 81.

16. Malm J., Terzhevik A., Bengtsson L., Boyarinov P., Glinsky A., Palshin N., Petrov M. Temperature and salt content regimes in three shallow ice-covered lakes. - Nordic Hydrology, 28, 1997, 99-128.

17. Mironov D.V., Golosov S.D., Zilitinkevich S.S., Kreiman K.D., Terzhevik A.Yu. Seasonal changes of temperature and mixing conditions in a lake. - In:

"Modelling air-lake interaction. Physical Background" (ed. S.S. Zilitinkevich), Springer-Verlag, Berlin, 1991, pp. 74-90.

18. Pavlov A.V. Regularities in Thermal Regime of Lakes in Permafrost Areas.

- Russian Geocryological Research, V. 1, 1995.

19. Satyanarayana A. N. V., Lykossov V.N. and Mohanty U.C. A study on atmospheric boundary layer characterisitics at Anand, India using LSP experimental data sets. – Boundary-Layer Meteorol., 2000, v. 96, 393-419.

20. Shnaydman V.A. Two-equation turbulence closure for quantative description of boundary layers.//Research activities in atmospheric and oceanic modeling. Edited by H. Ritchie. Report No. 32. April, 2002.

21. Simoes, F. 1998. "An Eddy Viscosity Model for Shallow-Water Flows," Water Resources Engineering 98, ASCE, NY, 1858-1863.

22. Tsuang B.-J., Tu C.-J., Arpe K. Lake parameterization for climate models. Report No. 316, Max Planck Institute for Meteorology, Hamburg, 2001, 72 pp.

23. Wuest A. and Lorke A. Small – Scale Hydrodinamics in Lakes// Annu. Rev.

Fluid Mech. 2003. 35:373-412.