Численное исследование гидрофизических процессов в стратифицированных озёрах

Белолипецкий Павел Викторович

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

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

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

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

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

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

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

                                                                           (1)

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

                                                                           (2)

Здесь T – температура воды, K_T(z) – коэффициент вертикального турбулентного обмена, F_H –теплообмен с дном, F_n – полный тепловой поток через свободную поверхность, F_I – приходящая коротковолновая радиация, b – коэффициент поглощения излучения; c_p – удельная  теплоемкость воды, r₀ – характерное значение плотности воды, H – глубина водоема.

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

                                                                                                      (3)

Здесь S – соленость воды, K_S(z) – коэффициент вертикального турбулентного обмена для солености, F_SH – массообмен с дном, F_S – поток через свободную поверхность. Необходимо также задать начальные распределения температуры и солености: .

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

                                      (4)

здесь , – напряжение трения ветра, K_min=0.02 см²/сек – минимальное значение коэффициента вертикального турбулентного обмена, , ,  ,  f – параметр Кориолиса.

                                         (5)

В описанной методике интенсивность вертикального турбулентного обмена определяется градиентом скорости и стратификацией.

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

                                                                                             (6)

где r₀=1.0254 г/см³, e₁=0,9753, e₂= – 0.00317, e₃=0,02737, T₀=17.5 ^oC,  S_o=35 ‰.

Согласно натурным измерениям, плотность соленой воды оз. Шира отличается от рассчитанной по формуле (6). Поэтому коэффициенты e₁, e₂ , e₃  были уточнены  согласно экспериментальным данным:

e₁=0,98452, e₂= – 0.007125, e₃=0,04112.

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

                                                                                                       (7)

где  F_R – суммарный радиационный теплообмен, F_e – теплоотдача испарением, F_c – конвективный теплообмен. Суммарный радиационный теплообмен F_R:

                

где  F_sc – коротковолновая солнечная радиация при ясном небе, m – облачность в долях единицы, T_a – температура воздуха на высоте 2 м  T_w  – температура водной поверхности (в градусах Цельсия).

Теплоотдача испарением:

          

где   f(W_z) – ветровая функция, e_m – давление насыщенного пара при температуре водной поверхности, W_z – модуль скорости ветра на высоте  z метров, e_z  –  давление пара на высоте z над поверхностью (обычно z=2м). Для определения ветровой функции можно использовать формулу Браславского–Нургалиева, учитывающую не только скорость ветра, но и соотношение между температурой воздуха и водной поверхности:

          

          

Расчет теплообмена конвекцией обычно производят, исходя из соотношения Боуэна, полученного в предположении подобия процессов тепло- и массопереноса

          

Напряжение  ветра определяется по формуле Давтян Н.А. [8]:

                                                                        (8)

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

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

           Математическая модель двумерных в вертикальной плоскости ветровых течений в замкнутых стратифицированных водоемах основывается на уравнениях турбулентных течений неоднородной жидкости в приближениях Буссинеска и "твердой крышки" [5]:

                

                                (9)

                                                                                                                 (10)

          

                                                    (11)

                

Здесь – составляющие скорости течения воды в направлениях Оx и Оz, ось Оz направлена вниз, – время, - плотность воды, - характерная плотность воды, – давление, – ускорение свободного падения, – температура воды, – соленость воды, - коэффициенты турбулентного обмена.

Введем функцию :

           ,                                                                                      (12)

тогда уравнения (9) приводятся к виду

                

                         (13)

Уравнения (10), (11), (13) дополняются начальными и граничными условиями. Граничные условия на водной поверхности (z=0)

                              

                                                                             (14)

на дне и на боковой поверхности (z=H или x=x_i)

                                              

                                 

                                                                  (15)

Здесь - внешняя нормаль к поверхности, – напряжение трения ветра, – полный поток тепла через свободную поверхность, – теплообмен с ложем водоема, – поток соли через свободную поверхность, – массообмен c ложем водоема.

Для определения коэффициента вертикального турбулентного обмена применяется формула Прандтля-Обухова:

                                    (16)

где - фоновое значение коэффициента вертикального турбулентного обмена, h – глубина верхнего квазиоднородного слоя, которая определяется по первой от поверхности расчетной точке, в которой выполняется условие

                

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

Коэффициенты горизонтального турбулентного обмена считаются постоянными. где - турбулентные числа Прандтля.

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

                 .                                                                                              (17)

С учётом этого приближения гидростатики уравнения для исследования двумерных в вертикальной плоскости течений  имеют следующий вид:

                                                          (18)

                 ,                                                                                                  (19)

                 ,                                                                                                  (20)

где - функция, описывающая положение водной поверхности.

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

           Интересным объектом для математического моделирования является озеро Шира. Озеро Шира относится к меромиктическим водоемам (в которых циркуляции охватывают только верхний слой) с высоким содержанием сульфата. Круглогодично в его воды поступает органическое вещество аллохтонного и автохтонного происхождения. Органика интенсифицирует процессы окисления и восстановления серы. Сложившийся комплекс физико-химических условий функционирования озера Шира не имеет аналогов. Знакомство с имеющейся литературой по озеру показывает, что внимание уделялось в основном гидрогеологии, меньше гидрохимии и совсем без внимания осталась гидробиология и экология озера; тем более отсутствуют знания о функциональном устройстве экосистемы (потоки вещества и энергии, сопряжения круговоротов, межпопуляционные взаимодействия, процессы на границах фаз и пр.).

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

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

Рис. 1. Характерные распределения по глубине температуры и солёности в различные периоды года.

Как видно, распределения существенно различные. Летом формируется устойчивая стратификация с двумя различными слоями верхним и нижним. Цель численных экспериментов заключалась в моделировании перехода от весеннего состояния к летнему. Как двумерная, так и одномерная модель показали хорошие результаты (рис. 2, 3). На основании чего был сделан вывод, что модели позволяют определять качественную вертикальную структуру, положение переходной зоны для температуры и солёности, значения величин в нижнем и верхнем слоях. На рис. 2, 3  приведены графики результатов расчётов одного из вариантов по одномерной и двумерной моделям и  натурные данные.

 

Рис. 2. Рассчитанные по одномерной модели и измеренные распределения температуры и солёности.

Рис. 3. Рассчитанные по двумерной модели и измеренные распределения температуры и солёности.

           В одномерной модели перемешивание происходит только под влиянием турбулентности, а в двумерной модели на перемешивание оказывают значительное влияние ещё и течения. Картина ветровых течений зависит от силы ветра, стратификации и геометрии водоема. Выполнена серия  расчетов, в результате были получены данные о характере течений, которые можно проиллюстрировать на некоторых характерных примерах. Для однородной жидкости в водоеме формируется одна циркуляционная зона, причем с увеличением длины водоема возрастает максимальное значение скорости течения воды. Влияние размеров водоема, силы ветра и стратификации было рассмотрено на следующих примерах. Пресноводное озеро длиной 500 м и максимальной глубиной 20 м, с линейным перепадом температуры от 7 до 4 градусов Цельсия в верхнем восьмиметровом слое и далее на глубине постоянном. В данном случае при скорости ветра более 9 м/с формируется одна циркуляционная зона с небольшими вихрями у дна (рис.6, 7). При скорости ветра 5 м/c в водоеме образуются две циркуляционные зоны (рис. 8, 9). Для соленой воды при том же перепаде температуры и постоянной солености и для скорости ветра 20 м/с формируются две циркуляционные зоны. Для пресноводного водоема длинной 5000 м, максимальной глубиной 20 м и тем же профилем температуры формируется одна циркуляционная зона с небольшими по размеру вихрями у дна. Для соленой воды такого же водоема формируются две циркуляционные зоны.

 

 

 

Выводы

           Предложен численный алгоритм  (на основе двумерной в вертикальной плоскости модели) для исследования ветровых течений в стратифицированных водоемах. Показано, что картина течения зависит от стратификации и силы ветра. Компьютерная модель позволяет определить к какому типу водоемов относится конкретное озеро: к димиктическому – с двумя периодами весной и осенью циркуляциями, охватывающими всю толщу воды от поверхности до дна, или к меромиктическому, в котором циркуляции охватывают поверхностный слой воды меньшей плотности. Озеро Шира является меромиктическим водоемом. Разработанная компьютерная модель позволяет прогнозировать динамику вертикальных распределений температуры, солености, плотности воды и при необходимости других компонентов водной экосистемы в зависимости от метеоданных с использованием грубых в горизонтальном направлении сеток («камерной» в горизонтальном направлении модели). 

 

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

 

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

1. Andrei G. Degermendzhy, Victor M. Belolipetsky, Tatiana A. Zotina, Ramesh D. Gulati. 2002. Formation of the vertical heterogeneity in the Shira Lake ecosystem: the biological mechanisms and mathematical model. “Aquatic Ecology”. 36:271-297.

2. Белолипецкий В.М., Белолипецкий П.В. Численное моделирование ветровых течений в стратифицированных водоемах методом расщепления. Гидростатическое приближение // Вычислительные технологии. 2006. Т.11, №5. С. 21-31.

3. Белолипецкий В.М., Генова С.Н. Одномерная модель вертикальной структуры озера. Температурный и солевой режимы озера // Труды Международной конференции «Вычислительные и информационные технологии в науке, технике и образовании». I том. – Павлодар. ТОО НПФ «ЭКО», 2006. С. 253-261.

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

5. Белолипецкий П.В. Численное моделирование двумерных в вертикальной плоскости ветровых течений в стратифицированных водоемах методом расщепления // Вычислительные технологии. 2005. Т. 10, № 5. С. 19-28.

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

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

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