Результаты вычислительных экспериментов при математическом моделировании водных экосистем

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

Экологическое состояние водных объектов зависит от большого числа разнообразных факторов и процессов:

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

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

Известен ряд общих требований, предъявляемых к каждой математической модели:

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

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

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

Описание математической модели. В качестве динамических переменных модели взяты концентрации зеленых водорослей ( CA 0), сине-зеленых водорослей ( CA 1), диатомовых водорослей ( CA 2), зоопланктона ( CZ ) , бактерий ( CB ), детрита ( CD ), растворенных в воде неорганического фосфора ( PS ), неорганического азота ( ^S ), органики ( POB ) и кислорода ( О 2).

В модели учитываются следующие процессы:

• -    рост микроорганизмов;

• -    выделение продуктов метаболизма;

• -    смертность микроорганизмов;

• -    процессы оседания;

• -    переходы по трофической цепи;

• -    процессы разложения;

• -    атмосферная реаэрация (выделение кислорода из воды);

• -    денитрификация (процесс восстановления нитратов до молекулярного азота, вызываемый бактериями);

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

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

Структура модели, описывающая функционирование экосистемы, приведена ниже (рис. 1).

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

dCA⁰ = ( mA 0 - RA 0 - MA 0) ■ CA 0 + a₀ ■ CA 1 ■ CA 0, dt                                        ⁰

dCA = (mA1 - RA1 - MA1) ■ CA1 - a, ■ CA1 ■ CA0, dt dCA2 = (mA2 - RA2 - SA2 - MA2) ■ CA2 - mZ' CZ , dt                                        Y1

dCZ dCZ- = (mZ - RZ - MZ) CZ, dt dCB = (mB - RB - MB) CB dt mZ ■ CZ Y 2

dCD = MA 0 ■ CA 0 + dt

+ MA 1 ■ CA 1 + MA 2 ■ CA 2 + MZ ■ CZ + MB ■ CB -

(стрелками обозначены направления потоков вещества между компонентами экосистемы)

mB ■ CB mZ ■ CZ

Y 4

dPS dt

- ( mA 0 - RA 0) ■ PP 0 ■ CA 0 - ( mA1 - RA1) ■ PP 1 ■ CA 1 -

- ( mA 2 - RA 2) ■ PP 2 ■ CA 2 + RZ ■ CZ ■ PP 3 + RB ■ CB ■ PP 4,

— = RA 0 ■ PN 0 ■ CA 0 + dt

+ RA 1 ■ PN 1 ■ CA 1 - ( mA 2 - RA 2) ■ PN 2 ■ CA 2 +

+RZ ■ CZ ■ PN3 + RB ■ CB ■ PN4, dPOB dt

mB CB

—m —— + h 0 ■ RA 0 ■ CA 0 + h b RA1 - CA1 +

Y 5

+ h 2 ■ RA 2 ■ CA 2 + h 3 ■ RZ ■ CZ + h 4 ■ RB ■ CB , dO² = K 1( O 2 O - O 2) +

+ K acc ( mA 0 ■ CA 0 + mA1 ■ CA 1 + mA 2 ■ CA 2) -

• - alf ( RA1 - CA1 + RA 2 ■ CA 2 + RZ ■ CZ + RB ■ CB ) - B 1 - mZ ■ CZ , где mAi – функции, описывающие рост; коэффициенты: RAi – дыхания, MAi – смертности, SAi – оседания, Yi – пропорциональности; Т - температура, ° C; t - время.

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

В микробиологических системах скорость роста, как правило, лимитируется концентрацией субстратов. Для описания процесса лимитирования применялась гиперболическая зависимость, предложенная Ж. Моно.

Предполагается, что рост зеленых и сине-зеленых водорослей лимитируется фосфором, а диатомовых – азотом и фосфором. Функции роста, смертности, освещенности и температурной зависимости, а также все входные данные взяты из работ [1; 2; 3].

Построенная математическая модель представляет собой задачу Коши для системы из десяти обыкновенных дифференциальных уравнений. Для численного решения задачи Коши применяется метод Рунге–Кутта четвертого порядка аппроксимации:

rrrr y ^r n+ 1 =y ^r n + ¹ τ( K 1 + 2 K 2 + 2 K 3 +K 4 ), 6

rr

K 1 = F ( t n , У n ), rr      ^r

K 2 = F ( t n + 2, y n + т K 2¹), rr     ^r

K з = F ( t n + 2, У n + т ^), rrr K 4 = F ( t n + T, y n + т K 3 ), ⁿ = 0,1,^, где y ^r – вектор-функция неизвестных; F – правая часть системы (1); т - шаг по времени; y ₀ задано.

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

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

Рис. 2. Коэффициенты модели и начальные значения системы

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

Программный комплекс реализован в среде Visual C++ 6.0 с использованием MFC (Microsoft Foundation Classes) – одного из наиболее удобных и мощных инструментов создания Windows-приложений. Комплекс обладает дружественным интерфейсом, удобен в работе и предназначен не только для математиков, но и для исследователей, не являющихся специалистами в программировании.

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

Дальнейшие расчеты, проведенные с использованием усовершенствованной математической модели, были посвящены исследованию общих тенденций сезонного изменения динамических переменных модели на основе экспедиционных данных по Новоселовскому плесу Красноярского водохранилища за 1998–2000 гг. Отметим, что для сравнения с результатами расчетов ниже будут представлены только те натурные данные, которые точно известны в конкретные моменты времени. За остальной временной период, в силу неполноты имеющейся информации, проведен сравнительный анализ по усредненным данным, также показавший качественное совпадение расчетных и экспериментальных данных.

В частности, сезонная динамика диатомовых водорослей демонстрирует качественную картину двух пиков цветения: летний с максимальной биомассой 5,9 мг/л и осенний с максимальной биомассой 2,27 мг/л, что хорошо согласуется с данными наблюдений [5]. Суммарная биомасса диатомовых водорослей по наблюдениям в июле и августе снизилась до 1,2...2,0 мг/л. Модельный расчет также выявил падение значений биомассы в указанный период (рис. 3).

Сопоставление модельного представления сезонного хода зеленых и сине-зеленых водорослей с экспериментальными данными [6] показало, что при общей зани-женности модельных концентраций (примерно в 1,5 раза) относительный временной ход был реконструирован верно (рис. 4). Годовой ход биомассы зоопланктона имеет единственный пик и соответствует максимальным значениям биомасс зеленых и сине-зеленых водорослей, что также соответствует теоретическим представлениям.

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

Годовой ход концентрации органического вещества имеет два выраженных максимума со значениями 0,29 мг/л в начале лета и 0,54 мг/л осенью. Эти периоды соответствуют максимальным значениям биомасс фито-и зоопланктона. Для сезонной динамики бактерий и детрита характерны два пика развития: первый приходится на середину лета, второй – на начало осени. Таким образом, детрит оказывает влияние на рост бактерий и стимулирует их продуктивность, что также соотносится с натурными данными.

Рис. 3. Сравнение экспериментальных данных с численными результатами для диатомовых водорослей (треугольниками обозначены натурные данные)

Полученные модельные расчеты дают неплохое со- впадение с экспериментальными данными, что свидетельствует об адекватности рассмотренной модели.

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

r

r

а и   а и r_z m

+ V = F ( t , x , U ) d t      d x

с соответствующими начальными и краевыми условиями:

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

Расчет проводился для времени от t = 0 до t = 365, т. е. в течение года, для Новоселовского плеса Красноярского водохранилища, а также для участков реки Енисея от Красноярской ГЭС вниз по течению со средней скоростью течения 1,2 км/ч на расстояниях от 10 до 100 км. В частности, результаты расчетов для бактерий показали, что при указанных условиях концентрация бактерий в выбранной части водоема изменяется значительно (рис. 5) (по мнению авторов, большой перепад значений во многом обусловлен течением).

vr

U ( t ,0) = U ₀( t ),

vr

U (0, x) = Ui( x), где V – скорость течения водоема; x – пространственная переменная, соответствующая длине водоема; t – время. Правая часть уравнений (2) соответствует правой части системы (1). При таком моделировании предполагается, что вещество равномерно распределено по ширине потока и движется со средней скоростью потока, т. е. рассматриваются осредненные по глубине и ширине водоема характеристики. Отметим, что данную модель целесообразно использовать в том случае, если длина водоема много больше его ширины.

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

У n + 1 s n             s n + 1 S n + 1

• Уг^2г + Vy - ^y j-i = F n .

τ h

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

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

Рис. 4. Сравнение экспериментальных данных с численными результатами для сине-зеленых водорослей (треугольниками обозначены натурные данные)

Отметим, что экосистема реки Енисея на участке, примыкающем к Красноярской ГЭС, сильно обеднена из-за разрушающего действия ее турбин и низкой температуры воды. Процесс самоочищения воды на этом наиболее загрязненном участке реки ослаблен. Зоопланктон в процессе деструкции вещества играет ничтожную роль. Бактериальная масса, которая здесь накапливается, подвергается интенсивной деструкции только ниже впадения реки Ангары.

На основании проведенных расчетов можно сделать вывод, что чем больше расстояние от Красноярской ГЭС, тем больше изменяется концентрация искомых компонентов, в частности происходит сдвиг пиков для диатомовых водорослей, бактерий и детрита. При расстоянии от ГЭС свыше 100 км модель описывает динамику развития водной экосистемы в пространстве менее адекватно. Очевидно, это обусловлено тем, что при изучении характеристик качества воды следует учитывать более сложные и полные процессы: гидродинамические (перемешивание, давление, ветровые и глубинные течения), теплопе-ренос и освещенность. Поэтому необходимо использование моделей, основанных на различных вариантах уравнений механики жидкости и теплопе-реноса и соответствующих граничных условиях [7].

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

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

Авторы благодарят профессора З. Г. Гольд за полезные советы при работе над статьей.