A mathematical model for simulating the biogeochemical processes in a freshwater lake.pdf В последние десятилетия под влиянием деятельности человека происходит заметное ухудшение экологического состояния крупных пресноводных озер планеты. Большой вклад в загрязнение и эвтрофирование водоемов вносят промышленные центры, сосредоточенные на их побережье и/или территории водосборного бассейна. В частности, в литоральной зоне озера Байкал в последние годы наблюдается активный рост водорослей рода спирогира, который связан с городскими бытовыми и промышленными стоками с высоким содержанием азота и фосфора. Важнейшими химическими элементами в озерной экосистеме выступают фосфор и азот, поскольку один из них является лимитирующим фактором развития биомассы водоема. Уровень концентрации фосфора в пресном озере имеет критически важное значение в предсказании вспышки цветения фитопланктона. В связи с этим проблема создания интегративных моделей, направленных на численное воспроизведение биохимических процессов фосфорного и азотного лимитирования, является актуальной для современного этапа развития математических методов в задачах лимнологии. Несмотря на длительную историю изучения первичной продукции в озере Байкал, вопрос о лимитирующих факторах развития его биомассы остается открытым. К примеру, ряд ученых [1, 2] считает, что развитие фитопланктона в Байкале лимитируется азотом. Но в то же время Голдман и др. [3] отмечают, что в летние периоды более отчетливым становится значение как фосфорного, так и азотного лимитирования. Экспериментальные исследования в Баргузинском заливе в конце августа 2002 г. [4] указывают на роль фосфора как первостепенного лимитирующего биогенного элемента, однако при поступлении фосфора в поверхностные воды лимитирующим становится азот. Очевидно, что при моделировании биохимических процессов в пресноводных экосистемах важно рассматривать не только наличие азота в воде, но и учитывать содержание фосфора, поскольку фосфатные группы являются основными струк- 1 Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 19-3160003. Б. О. Цыденов 54 турными элементами нуклеиновых кислот в растительных клетках и участвуют в энергетическом обмене [5]. Целью данной работы является разработка и описание математической модели экосистемы пресноводного озера с учетом фосфорного и азотного лимитирования для воспроизведения в нем физико-биологических процессов, а также анализ и обсуждение результатов расчетов, проведенных для условий озера Байкал в летний период. Математическая модель Математическая модель для воспроизведения физико-биологических процессов в пресноводном озере состоит из термогидродинамического [6, 7] и экосистемного блоков. Экосистемный блок основан на моделях Фашама и др. [8], Фен-нель и др. [9], Хофмана и др. [10] и включает в себя 7 прогностических переменных: нитрат (1), аммоний (2), хлорофилл а (4), фитопланктон (5), зоопланктон (6), мелкий (7) и крупный (8) нитратный детрит. Для описания фосфорного лимитирования в пресноводной среде дополнительно введены в модель 3 переменные согласно постановке Гана и др. [11]: фосфат (3), мелкий (9) и крупный (10) фосфатный детрит. Схематическая диаграмма, демонстрирующая связь между компонентами экосистемы, представлена на рис. 1. Интенсивность роста фитопланктона в модели определяется через модули расчета фосфорной и азотной составляющих согласно их роли в лимитировании биомассы водоема [4]. В качестве индикатора лимитирования часто используют молярное отношение содержания минерального азота (N) к содержанию минерального фосфора (P). Согласно Редфилду [12], при соотношении N:P < 16 первичная продукция фитопланктона лимитируется азотом, а при N:P > 16 - фосфором. Рис. 1. Концептуальное представление экосистемного блока математической модели. Стрелками показаны биологические процессы между компонентами модели [11] Fig. 1. Conceptual representation of the ecosystem module of the mathematical model. Arrows indicate the biological processes between components of the model [11] Математическая модель для воспроизведения биогеохимических процессов 55 Биохимические процессы в пресноводном озере описываются с помощью реакционно-конвективно-диффузионных уравнений д[N] + du[N] + dw[N] _ D д[N] A + d f D d[N] dt dx dz dx I x dx J dz f z dz -^maxf (I){РЬУЩ CT [ N ] 1 - + CTr [ N ] [ P] NkN + [N] 1 + [A]/kA P [N] + [A] kp + [P] d[ A] + n[ A]; (1) d[A] + du[A] + dw[A] =d_fD d[A]A + d dt dx dz dx I x dx J dz v z dz -hmaxf (I)[phyto] \\ CT [ A] [ A] [ P] JNkA + [A] +CTp [N] + [A] kp + [P] -n[A] + lBM [ Zoo] + lE-[Phyt°]2 P[Zoo] + rSDN [SDN ] + [ LDN ]; (2) kPhyto + [ Phyo d[P] + du\\P] + dw[P] _d_fD dt dx dz dx v x -hnaxf(I)[Phyto]\\ CT [ N ] 1 d[ P] V-ii 2dz *p] dx J dz' v dz [A] \\ + -+-- Г rPN +CTi [ P] PkP + [ P] P +lBMrPN [Zoo] + lErPN [ Phyto]2 kPhyto + [ Phyto] - P[Zoo] + rSDP [SDP] + rLDP [ LDP]; (3) d[Chl] du[Chl] dw[Chl] d f d[Chl] A d f d[Chl] A r^, „ --- + -^-- + -^--_-\\ Dx--- 1+-\\ D.--- \\ + pChl p]Chl ] dt dx dz dx v dx J dz v dz [Zoo] [Chl] ; (4) [Phyto] -mPhyto [Chl] - t([SDN] + [Phyto]) [Chl] - g {Phyto h kPhyto + Wo] d[Phyto] + du[Phyto] + dw[Phyto] df d d[Phyto] \\ + d fd d[Phyto]' dt dx dz dx \\ x dx J dz v z dz +|4 Phyto] - mPhy^ [ Phyto] - SDN ] + [ Phyto])) Phyto] - gn [ Phyto]2 ' kPhyto + [Phyto] -[Zoo]; (5) d[Zoo] + du[Zoo] + dw[Zoo] df d[Zoo] A+ d f d d[Zoo]' + dt dx dz dxf x dx J dz v z dz +gm [ Phyto]2 'kPhyto + [ Phyto] 7e[Zoo] - Ibm [Zoo] - Ie [ Phyto]2 ' kPhyto + [Phyto] rP[ Zoo] - mZoo [ Zoo]2; (6) d[ SDN ] + du[SDN ] + dw[ SDN ] dt dx dz _±(d^dfD]\\+±(D dfD] \\ + [ Phyto]2 -(1 -P)[ Zoo] + Sr v x Sr J dz v z dz J ~max kPhyto + [Phyto]2 +mZo0 [Zoo]2 + mPhyto [Phyto] - t(SDN] + [Phyto]) [SDN] - rSDN [SDN]; (7) d[LDP] + du[LDP] + dw[LDP] _dfD d[LDP]^ + d fD d[LDP] dt dx dz dx v x dx J dzf z dz + 56 Б. О. Цыденов _d_fD d[LDN] dx v x dx d[LDN] du[LDN] dw[LDN] _ dt dx dz + d-z f Dz ^J + t(SDN] + [Phyto] )2 - ^ldn [LDN]; (8) dt ■М 'D d[SDP] ]+-( d d dx 1 v dx J dz 1 +mZoorPN [Zoo]2 + mPhytorPN dz d[ SDP] + du[SDP] + dw[ SDP] dx + -+-, то oN = 0 и о+ = 1 (фосфорное kN + [N] 1 + [A]/kA kA + [A] k+ + [P] лимитирование), в противном случае - oN = 1 и o+ = 0 (азотное лимитирование). Функция f(I) представляет связь фотосинтеза и света: f (I) аі I 2 , 2 т-2 Mmax +а 1 (13) L I _ Is • PAR • exp -d\\ k + kchl I Chl(z)dz d (14) где IS - приходящий на поверхность воды свет; PAR - доля света для фотосинтеза = 0.43 [9]. Параметру IS в модели соответствует коротковолновая радиация [9], которая вычисляется по формуле (8) из [14]. Математическая модель для воспроизведения биогеохимических процессов 57 Участвующая в синтезе хлорофилла доля роста фитопланктона определяется как (15) р Ѳтах V-jPhyto] PchI aI[Chl ] Скорость нитрификации вычисляется [15]: n = n тах 1 - max 0, I -10 kI +1 -10 (16) где nmax - максимальная скорость нитрификации. Параметры, входящие в расчетные формулы экосистемного блока модели, приведены в таблице. Значения биогеохимических параметров модели Обозна чение Наименование Значение Цп Скорость роста фитопланктона при 0 °C 0.59 сут-1 kwater Коэффициент ослабления света для воды 0.04 м-1 kChI Коэффициент ослабления света для хлорофилла 0.025 (мг Chl)-1 • м-2 а Начальная крутизна кривой P-I 0.025 мг C (мг Chl х х Вт-м-2-сут)-1 ѳ ^max Максимальное соотношение хлорофилла к биомассе планктона 0.054 мг Chl (мг C)-1 rPN Клеточное соотношение P:N 0.0625 kN Константа полунасыщения фитопланктона по NO3 0.8 ммоль • м-3 kA Константа полунасыщения фитопланктона по NH4 0.8 ммоль • м-3 kp Константа полунасыщения фитопланктона по PO4 0.05 ммоль • м-3 mPhvto Смертность фитопланктона 0.15 сут-1 £max Максимальная скорость питания зоопланктона 0.6 сут-1 в Эффективность усвоения зоопланктона по азоту 0.75 kPhvto Константа полунасыщения питания фитопланктона 1 ммоль N • м-3 lBM Скорость экскреции зоопланктона за счет основного обмена 0.1 сут-1 h Скорость экскреции зоопланктона за счет видовой особенности 0.1 сут-1 mZoo Смертность зоопланктона 0.025 сут-1 х х (ммоль N • м-3)-1 rSDN Скорость реминерализации мелкого детрита для азота 0.03 сут-1 rLDN Скорость реминерализации крупного детрита для азота 0.01 сут-1 rSDP Скорость реминерализации мелкого детрита для фосфора 0.075 сут-1 rLDP Скорость реминерализации крупного детрита для фосфора 0.025 сут-1 T Скорость коагуляции 0.05 сут-1 nmax Максимальная скорость нитрификации 0.05 сут-1 І0 Порог ФАР для ингибирования процесса нитрификации 0.0095 Вт • м-2 kI Полунасыщение ФАР для ингибирования процесса нитрификации 0.1 Вт • м-2 Для биологических переменных модели заданы следующие начальные условия: - значения концентрации нитратов, аммония, фосфатов, зоопланктона и хлорофилла а приняты равными 5.0, 4.0, 0.4, 0.3 и 0.3 ммоль N м-3 соответственно; Б. О. Цыденов 58 - концентрация для всех детритных компонентов в начальный момент времени составляет 0.1 ммоль N • м-3; - значение биомассы фитопланктона задано исходя из соотношения 1.59 для [Chl]/[Phyto\\. Граничные условия для уравнений модели имеют вид: а) на поверхности водоема d[N] = 0 d[A] = 0; d[P] = 0; д\\СЩ = 0; д[Phyto] 0; dz ’ dz ’ dz ’ dz ’ dz ’ d[Zoo] 0, d[SDN] = 0; d[LDN] = ^ d[SDP] = ^ d[LDP] = ^ dz ’ dz ’ dz ’ dz ’ dz ’ б) на твёрдых границах ад=0; dLH=0; d!Pl=0; dCH=0; адм=0; dn dn dn dn dn d[Zoo] 0, d[SDN] = 0. d[LDN] = ^ d[SDP] = ^ d[LDP] = 0 dn dn dn dn dn где n - направление внешней нормали к области; в) на границе впадения реки в озеро [N] = 5.0 ммоль N • м-3; [A] = 4.0 ммоль N • м-3; [P] = 0.4 ммоль N • м-3; [Zoo] = = 0.3 ммоль N • м-3; [Chl] = 0.3 ммоль N • м~3; [SDN] = 0.1 ммоль N • м-3; [LDN] = 0.1 ммоль N • м-3; [SDP] = 0.1 ммоль N • м-3; [LDP] = 0.1 ммоль N • м~3; г) на открытой (правой) границе задаются условия радиационного типа: А 0 dt ф dx (ф = [N],[A], [P], [Chl], [Phyto], [Zoo], [SDN], [LDN], [SDP], [LDP]), фазовая скорость cф рассчитывается из пространственных и временных тенденций внутри области [16]. Условную запись разностной схемы для конечного объёма P для уравнений (1) - (10) можно представить в обобщенном виде: ФP -Фр At 2[-(С*)P -(C*)P +(D*)p + (D*)p +(XФ)P] + +1 [-(* )p-(C* )p+(DX Pp +(D* )p +№>V ] ]) где Ф = [N], [A], [P], [Chl], [Phyto], [Zoo], [SDN], [LDN], [SDP], [LDP]; (C*) , (C*) и (D*) , (D*) - аппроксимация конвективных и диффузионных членов соответственно; (Srt,)p - источниковый член; ФР соответствует значению Ф в центре конечного объёма на временном слое t, а ФР - на слое t+At. Видно, что нестационарный член в (17) аппроксимируется схемой Кранка - Николсон. Для аппроксимации диффузионных членов используется центрально-разностная схема, а конвективных членов - противопотоковая схема QUICK второго порядка точности [17]. Математическая модель для воспроизведения биогеохимических процессов 59 Конечная система алгебраических уравнений, соответствующая дифференциальной задаче имеет вид ap фр = Яр фe + aw Фw + aN ФN + aS ФS + b , (18) где индекс P обозначает центр текущего конечного объёма (i, j), а индексы E, W, N и S (согласно «географической» системе обозначений) - центры соседних объёмов (i+1,j), (i-1,j), (ij+1) и (i,j-1) соответственно; aE > 0, aW > 0, aN > 0, aS > 0 о ( 0 AxAzЛ и aP = aE + aW + aN + aS + aP > 0 I aP = -^-I. Решение системы (18) осуществляется методом неполной факторизации, устойчивость которого обеспечивается диагональным преобладанием [18]. Важно заметить, что линеаризация источникового члена (S0 )P осуществляется следующим образом [19]: (So )P ~ SC + SP®P, где SC - постоянная составляющая (Sф)р , а SP - коэффициент (SP не есть значение (Sф)р в точке P). Если коэффициент SP всегда отрицателен или меняет знак, то для усиления диагонального преобладания матрицы коэффициентов системы (18) переносится в левую часть уравнения (18). Детальное описание численного метода решения уравнений вида (18) приведено в работах [6, 20]. Коэффициенты вертикальной диффузии определяются как D =^т_ + ѵ , z PrT Pr где vT - турбулентная кинематическая вязкость (рассчитывается по Л-ю-модели Уилкокса [22]); v - молекулярная кинематическая вязкость воды; Pr и PrT - молекулярное и турбулентное числа Прандтля. Коэффициенты горизонтальной диффузии задаются постоянными [22]: Dx = 2.5 м2/с. Расчеты проводились на серверной платформе Supermicro с графическим ускорителем Tesla C2050. Область исследования Областью исследования является вертикальное сечение Баргузинского залива озера Байкал, начало системы координат совпадает с устьем р. Баргузин 53°25'30" с.ш., 108°59'50" в. д. [23] (рис. 2, а). Данные по рельефу дна взяты из работы [24] (рис. 2, b). Вычислительная область имеет протяженность 20 км и глубину 100 м (рис. 2, b). Глубина открытого участка речного стока (на левой границе) составляет 7.5 м. Расчётная область покрыта равномерной ортогональной сеткой с шагами hx = 25 м и hz = 2.5 м. Шаг по времени - 30 с. Исследование сходимости численного решения уравнений вида (18) по шагам сетки и времени проводилось в работах [20, 25, 26]. В качестве критерия сходимости итерационного процесса использовалось условие малости невязки (0.06 ммоль м-3) на глубине 10-30 м, причем с увеличением времени моделирования её изменение в вертикальном направлении усиливается (рис. 3, а2). Эти результаты качественно согласуются с распределением хлорофилла а, полученным методом спектрофотометрии в центральной части Баргузинского залива 24-31 августа 2002 г.: пиковое значение хлорофилла детектировано на глубине 10 м [4]. -80 -100- 0 Щ*щш\\ ^ ) і ' ' " аі ■ ■ bl : : b2 0 -20 -40 -60 -80 -100 0 -20 -40 -60 -80 -100 0 -20 -40 -60 -80 -100 Рис. 3. Распределение хлорофилла а [мг • м-3] (а), зоопланктона [ммоль N • м-3] (b), нитрата [ммоль • м-3] (c), фосфата [ммоль • м-3] (d) и аммония [ммоль • м_3] (e) на 24-е (і) и 28-е (2) сутки моделирования Fig. 3. Distribution of chlorophyll a [mg • m 3] (а), zooplankton [mmol N • m 3] (b), nitrate [mmol • m-3] (c), phosphate [mmol • m-3] (d), ammonium [mmol • m-3] (e) on the 24 (1) and 28 (2) days Зоопланктон сосредоточен в верхней 25-метровой толще водоема (рис. 3, b), его наибольшая концентрация отмечается в приповерхностном слое. В открытой части водоема (по сравнению с приустьевым пространством) рост биомассы зоопланктона происходит интенсивнее (рис. 3, b). Такая картина развития зоопланктонных популяций в конце лета имеет качественное и количественное различие от Б. О. Цыденов 62 весенней динамики, связанной с существованием термобара [31, 32]. На 28-й день концентрация зоопланктона возрастает до 0.52 ммоль N • м-3 (рис. 3, Ъ2). Содержание биогенных элементов в воде имеет обратную корреляцию с распределением биомассы планктона: на участках водоема с высокой концентрацией фито- и зоопланктона наблюдается снижение NO3, PO4 и NH4 (рис. 3, c, d, e). Результаты моделирования показывают, что уменьшение количества NO3 в присклоновой области водоема происходит быстрее, чем в пелагической зоне (рис. 3, c), при этом развитие гидробиологических процессов ведет к усилению данной тенденции (рис. 3, c2). Повышенное содержание PO4 (>0.5 ммоль • м-3) на 24-й день заметно в 30-50-метровом слое (рис. 3, d1). На 28 день PO4 (рис. 3, d2) в водах открытого озера принимает горизонтально однородное распределение, количество PO4 возрастает с глубиной. Наибольшая концентрация NH4 отмечается в 20-30-метровой толще (рис. 3, e), выше и ниже указанного слоя происходит уменьшение количества NH4. На 24-й день наблюдается образование «языка» с высоким содержанием NH4 у нижней границы рассматриваемой области на расстоянии 4.5-8.5 км от устья речного притока (рис. 3, el). Мониторинговые исследования химического состава воды в оз. Байкал в летний период также фиксируют: 1) различие в содержании биогенных элементов в мелководных и пелагических участках озера (2002-2007 гг.) [33]; 2) увеличение концентрации нитратов и фосфатов с ростом глубины в пела-гиали (1993-2011 гг.) [34]. Заключение Построенная математическая модель позволяет воспроизводить биогеохимические процессы в пресноводном озере. Проведенные на примере Баргузинского залива оз. Байкал расчеты показали, что в последней декаде августа 1) повышенное содержание хлорофилла а в пелагиали водоема локализовано на глубине 10-30 м; 2) зоопланктон сосредоточен в верхнем 25-метровом слое, причем в открытой части озера увеличение его биомассы происходит быстрее; 3) на участках водоема с высокой концентрацией фито- и зоопланктона наблюдается снижение количества биогенных элементов. Полученные результаты качественно согласуются с данными мониторинговых исследований химического состава воды в озере Байкал в летний период.

Скачать электронную версию публикации