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

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

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

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

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

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

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

Предполагается, что весь азот в водохранилище можно представить как неорганический (где N₁ – переменная) и органический (где P, B, Z, N₀ – переменные). Органический азот, в свою очередь подразделяется на живой (где P, B, Z – переменные) и неживой (где N₀ – переменная).

Теперь попытаемся описать процессы, происходящие между этими динамическими переменными.

Сбросы в водоем органического и неорганического азота в сутки (равны С, где С – некоторая константа).

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

Фиксация свободного азота. Запас связанного азота в водоеме может пополниться путем усвоения молекулярного азота воздуха свободноживужими фиксаторами азота (разными видами Azotobacter, Clostridium pasteurianum) и некоторыми сине-зелеными водорослями. Поскольку последние в водоемах распределены очень широко, то фиксация азота из воздуха за счет их развития может достигать в природных водах очень большой величины.

Здесь рассматриваются процессы азотфиксации фито- и бактериопланктона. Для наших модельных условий запишем:

• – азотфиксация фитопланктоном за сутки равна kP;
• – азотфиксация бактериопланктоном за сутки равна kB.

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

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

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

Минерализация органического азота. Образовавшиеся в водоеме белковые вещества растений и животных в дальнейшем, после отмирания организмов, подвергаются минерализации с помощью бактерий. Участие бактерий в азотном цикле является необходимым условием круговорота. Азот распадающихся белковых веществ выделяется в первую очередь в виде аммиака. Таким образом, в результате этого процесса водоем пополняется минеральными соединениями и обедняется органическими. Заметим, что потребность фитопланктона в биогенных элементах обеспечивается за счет регенерации в умеренных широтах и составляет 18–25 %. При моделировании этот процесс разбивается на две стадии:

• минерализация органического азота бактериопланктоном за сутки равна kBN₀ (по закону взаимодействующих масс);
• поступление неорганических соединений азота в водоем за счет минерализации за сутки равно kB (линейная зависимость).

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

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

Этот процесс называется косвенной динитрификацией .

Потери азота возможны также в результате деаммонификации, NH₄ → N₂, а также при образовании летучих окислов азота с участием мочевины, некоторых аминокислот и разнообразных нестойких азотсодержащих соединений.

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

Разделим и этот процесс на две стадии:

• потери неорганического азота в водоеме за счет денитрификации за сутки равны kBN₁ (по закону взаимодействующих масс);
• переход азота в газообразную форму (выход из системы) в сутки равен kB (линейная зависимость).

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

Экспериментальным путем показано, что взятые для опыта животные-детритофаги (пресноводные планктонные рачки-фильтраты, дафнии, пресноводные моллюски лужанки и насекомые жучки-кожееды) обеспечивают достаточно высокий уровень поступления биогенных элементов, доступных для усвоения звеном первичных продуцентов. Так, в эвтрофных водоемах (чем является и наш объект) выделения аммиака зоопланктоном почти на 56 % обеспечивают ежегодные потребности фитопланктона в азоте. Все это свидетельствует о необходимости учета этого процесса в круговороте азота.

Таким образом, для модели выделения зооплактона за сутки равны kZ (линейная зависимость).

Смертность. При моделировании учитывается также и смертность живых компонентов (P, Z, B), т.е. переход в мертвую органику (N₀):

• отмирание фитоплактона за сутки равно kP (линейная зависимость);
• отмирание бактериопланктона за сутки равно kB (линейная зависимость);
• отмирание зоопланктона за сутки равно kZ + k₁Z (линейная зависимость), где под k₁Z подразумевается выедание зооплактона хищниками более высокого трофического уровня, другими словами, выход из системы.

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

• выедание зоопланктоном фитопланктона за сутки равно kZP;
• выедание зоопланктоном бактериопланктона за сутки равно kZB;
• потребление зоопланктоном детритного азота за сутки равно kZN₀.

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

Рассмотрим простейший случай круговорота биогенного элемента в водоеме, где в качестве фазовых переменных приняты следующие величины: P – количество фитопланктона; N₀ – количество органического азота; N₁ – количество неорганического азота в водоеме.

Блок-схема круговорота азота в Куйбышевском водохранилище (три фазовые переменные) представлена на рис. П.1.

На ее основе составлена модель, состоящая из трех обыкновенных дифференциальных уравнений:

(1)

где k₀ – коэффициент азотфиксации фитопланктона;
k₁ – коэффициент потребления фитопланктоном неорганического азота;
k₂ – смертность фитопланкона;
k₃ – коэффициент трансформации органического азота;
k₄ – выедание фитопланкона зоопланктоном;
k₅ – сток неорганического азота через плотину;
k₆ – сток органического азота через плотину; C₁ – поступление неорганического азота в водоем; C₂ – поступление органического азота в водоем.

Нахождение реальной стационарной точки модели. Объем Куйбышевского водохранилища равен 58 км³. Средняя концентрация неорганического азота (N₁) ~ 0,4 мг/л или в переводе на все водохранилище N₁ = 23200 т. Для органического азота средняя концентрация ~ 0,657 мг/л, т.е. в пересчете на весь водоем N₀ = 38106 т.

Теперь переведем биомассу фитопланктона P в азотные единицы. Средняя концентрация фитопланктона ~ 9.334 мг/л. В переводе на сухой вес эта цифра принимает значение 2,3335 мг/л. Общий азот составляет в среднем 9 % от сухого веса, т.е. 0,21 мг/л (азота), или для всего водохранилища P = 12181 т.

Таким образом, получаем стационарную точку модели (1) (в тоннах азота):

(2)

Поиск коэффициентов. За год в Куйбышевское водохранилище поступает 426900 т общего азота, т.е. C₁ = 1185,8 – C₂ (т/сут). (3)

За год сток общего азота – 385900 т, т.е. т или 1071,944 т/сут, отсюда

k₅ = 0,046 – (сут^–1).(4)

В Куйбышевском водохранилище за год аккумулируется 41000 т азота, т.е. 113 т/сут. Следовательно, отсюда

k₁ = 2,856 · 10^–8 (С·т·сут)^–1. (5)

Коэффициент прижизненных выделений и смертность фитопланктона составляет

k₂ = 0,35/360 = 0,972 · 10^–3 сут^–1.. (6)

Коэффициент азотфиксации фитопланктона k₀ = 0,9458 год^–1 или 2,63 · 10^–3 сут^–1.

Оставшиеся коэффициенты C₂, k₆, k₄, k₃ верифицируем по условию стационарности (2) модели (1), используя соотношениями (3)–(6):

k₄ = 0,0109 сут^–1, (7)

k₃ = (8)

Относительно двух коэффициентов k₆ и C₂ сделаем два предположения: k₆ = k₄ – k₀ = 8,304 · 10^–3 сут^–1, тогда из (4) получим k₅ = 0,0324 сут^–1; C₁ = C₂, тогда C₂ = 592,9 т/сут, а по (8) получим k₃ = 5,404 · 10^–4 (С·т·сут)^–1.

Из третьего уравнения системы (1) (стационарный вариант), уточняя коэффициент C₁ и используя (5)–(8), получаем C₁ = 576,4 т/сут.

Список коэффициентов:

Исследование на устойчивость. Исследуя на устойчивость стационарную точку модели (1) (приравнивая к нулю ее правые части), получаем следующие необходимые условия устойчивости:

(9)

Условия (9) для численных значений коэффициентов выполняются.

Система (1), а также все последующие модели были исследованы на устойчивость численным образом с помощью определителей Гурвица.

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

На рис.П.2 представлена блок-схема процесса круговорота азота в Куйбышевском водохранилище (четыре фазовые переменные).

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

(10)

где: k₇ – смертность зоопланктона; k₉ – коэффициент потребления зоопланктоном органического азота;
k₁₀ – коэффициент выедания зоопланктона хищниками более высокого трофического уровня; – измененный коэффициент выедания фитопланктона зоопланктоном.

Нахождение реальной стационарной точки модели. Стационарное значение для зоопланктона обозначим через Z⁰. Средняя биомасса зоопланктона составляет 448 мг/м³. Предположим, что процентное соотношение по азоту в зоопланктоне то же, что и в бактериях, тогда значение средней концентрации зооплактона в азотных единицах – 22,4 мг/м³. Для всего водохранилища – 1299,2 т.

Таким образом, получаем стационарную точку модели (10) (в тоннах азота).:

(11)

Поиск коэффициентов. Приравнивая к нулю правые части системы (9), из первого уравнения, используя (11), получаем: 6,012 10^–7 (С·т·сут)^–1. В эвтрофном водоеме пища зоопланктона содержит 5–15 % фитопланктона и 10–20 % детрита, причем Из третьего уравнения системы получим k₇ = 0,32 сут^–1. Из второго уравнения – k₁₀ = 0,743 сут^–1.

Список коэффициентов:

На рис. П.3 представлена блок-схема процесса круговорота азота в Куйбышевском водохранилище при наличии пяти фазовых переменных.

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

(12)

где k₈ – коэффициент выделения бактериопланктоном минерального азота; k₁₁ – коэффициент потребления зоопланктоном бактериопланктона; k₁₃ – коэффициент азотфиксации бактериопланктона; k₁₂ – коэффициент выделения бактериопланктоном молекулярного азота; k₁₅ – смертность бактериопланктона; k₁₆ – коэффициент денитрификации; – измененный коэффициент трансформации органического азота.

Нахождение реальной стационарной точки модели. Найдем стационарное значение для бактериопланктона, средняя биомасса которого составляет 1,9 мг/л. 10% от сухого веса составляет азот, а 75% – вода, т.е. получим 0,471 мг/л. Для всего водохранилища B = 2731,8 т азота.

Таким образом, получаем (в тоннах азота) стационарную точку модели (12):

Поиск коэффициентов. Приравняем к нулю правые части системы (13). Из второго уравнения получим k₁₂ = 0,115 сут^–1. По результатам исследований k₁₃ = 0,01 · k₀ ≈ 3,63 · 10^–4 сут^–1. Предположим, k₁₅ = k₇ = 0,3 сут^–1, где k₁₅ – смертность бактериопланктона, а k₇ – смертность зоопланктона, т.е. k₁₅ = 0,3 сут^–1. Из четвертого уравнения получим k₁₃ = 7,595 · 10^–7 (С·т·сут)^–1.

Для процесса денитрификации примем k₁₆ = k₁·10^–2, отсюда k₁₆ = 2,856 10^–10 (С·т·сут)^–1. Из пятого уравнения системы (12) получим k₈ = 0,0509 сут^–1. Из третьего уравнения следует: k₁₄ = 6,132 · 10^–5 сут^–1.

Список коэффициентов:

Блок-схема круговорота азота в Куйбышевском водохранилище уже из расчета шести фазовых переменных изображена на рис. П.4.

На основе блок-схемы составлена модель, состоящая из шести обыкновенных дифференциальных уравнений:

(13)

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

– видоизмененный коэффициент выделения бактериопланктоном аммонийного азота; – видоизмененный коэффициент денитрификации; поступления в водоем аммонийного азота; поступления в водоем нитратного азота.

Нахождение реальной стационарной точки модели. Найдем стационарное значение для аммонийного и нитратного азота. Средняя концентрация аммонийного азота ~ 0,226 мг/л, а в пересчете на все водохранилище Тогда Таким образом, получаем (в тоннах азота) стационарную точку модели (13): P⁰ = = 12181; Z⁰ = 1299,2; B⁰ = 2731,8;

Поиск коэффициентов. Приравняем к нулю правую часть системы (13).

По условию

Отсюда где – общий минеральный азот.

Пусть тогда (С·т·сут)^–1.

Должны выполняться следующие условия:

1. откуда (С·т·сут)^–1.
2. .
3. Пусть
   Тогда сут^–1 и сут^–1.
4. или
   Предположим: k₁₇ = (С·т·сут)^–1; сут^–1; т/сут (по условию 2).

Из шестого уравнения системы (13) находим k₁₆ = 0,0138 сут^–1. Из третьего уравнения получим сут^–1. В пятом уравнении уточним сут^–1.

Список коэффициентов:

По своей структуре сезонные модели не отличаются от моделей с постоянными коэффициентами: они описывают те же процессы и используют те же численные значения коэффициентов. Существенным отличием является только введение в ряд параметров (в явном виде) времени как аргумента.

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

Перепишем модель (1) уже при условии переменных параметров:

(14)

где С₁(t) и С₂(t) – сбросы в водоем минерального и органического азота (рис. П.5), определены согласно поступлениям азота в Куйбышевское водохранилище, T(t) – температура воды и I(t) – освещенность воды Куйбышевского водохранилища определены соответственно по рис. П.6.

Во всех имитационных расчетах временной интервал брался равным 360 дней, а временной шаг – 5 дней.

Аналогичным образом перепишем остальные уравнения (10), (11) и (12):

(15)

(16)

Модель из шести уравнений

(17)

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

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

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

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

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

Если водная экосистема описывается уравнениями состояния

(19)

и по компонентам векторов начальных состояний:

(20)

Начальные условия для уравнений (19) и (20) имеют вид (16) и (17) соответственно:

(21)

(22)

Уравнения (19) и (20) решаются совместно с уравнениями состояния (13), а полученные функции чувствительности применяются для изучения влияния антропогенных факторов.

Дифференцирование по параметрам и начальным условиям уравнений (14) приводит к дифференциальным уравнениям чувствительности 2-го порядка:

(23)

При начальном условии (24)

Функции чувствительности 2-го порядка (23), (24) используются для изучения синергизма антропогенных воздействий.

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

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

На рис. П.8 изображены функции чувствительности 1-го порядка по коэффициенту потребления минерального азота фитопланктоном для системы (16). Следует отметить наибольшую реакцию системы на изменение параметра потребления фитопланктоном нитратного азота.

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

Проведем теперь качественное исследование модели (17).

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

Таблица 1

Интерпретируя факторы, представленные в таблице П.2, заметим, что фактор 1 связан в основном с минеральными формами азота, а фактор 2 – с органическими формами азота и продуцентами. График проекций для 1 года (360 суток, шаг 5 суток) по первым двум главным факторам представлен на рис. П.9.

Годовой фазовый портрет модели (17) отражает почти всю информацию, которую несут в себе сложные графики на рис. П.7 – П.8. Так, с марта по июнь возрастает роль фактора 2, четко наблюдается два пика – весенний (больший) и осенний (меньший). "Прогонка" модели (17) на 5 лет и последующее применение вышеизложенного многомерного структурного подхода к качественному исследованию результатов моделирования дает возможность установить устойчивый годовой цикл колебаний обобщенных фазовых координат.