CC BY
7
УДК 517.9
DOI: 10.14529/ mmp230202
ИЗУЧЕНИЕ КАЧЕСТВЕННЫХ ЗАКОНОМЕРНОСТЕЙ ПРОЦЕССА ЭВТРОФИРОВАНИЯ МЕЛКОВОДНОГО ВОДОЕМА НА ОСНОВЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ БИОЛОГИЧЕСКОЙ КИНЕТИКИ
Ю.В. Белова1, Е.О. Рахимбаева1, В.Н. Литвинов1, А.Е. Чистяков1, А.В. Никитина1, А.М. Атаян1
1 Донской государственный технический университет, г. Ростов-на-Дону, Российская Федерация
Статья посвящена моделированию процессов эвтрофирования мелководного водоема с использованием математической модели биологической кинетики, которая базируется на системе нестационарных уравнений конвекции-диффузии-реакции с нелинейными членами. Данная система учитывает такие составляющие, как разложение детрита, гравитационное оседание примесей, микротурбулентную диффузию, движение водного потока. Изучение продукционно-деструкционных явлений и процессов в мелководном водоеме осуществляется путем введения в математическую модель скорости роста фитопланктона и бактерий. Также есть возможность контроля динамики процессов биологической кинетики при достаточно малом притоке соединений серы и биогенных веществ. Проведена линеаризация непрерывной задачи и построен дискретный аналог на базе расщепления исходной трехмерной задачи на двумерную и одномерную. На увеличение точности моделирования исследуемых явлений и процессов повлияло применение линейной комбинации схем «кабаре» и «крест» для построения дискретной двумерной модели. В статье приведены результаты диагностического моделирования процессов возникновения сероводородного заражения, изучаются процессы самоочищения мелководного водоема. Разработанная математическая модель и ее численная реализация соответствуют современным представлениям о функционировании гидробиоценоза мелководного водоема.
Ключевые слова: Азовское море; эвтрофирование; математическое моделирование; явно-неявная 'разностная схема; погрешность аппроксимации; декомпозиция 'расчетной области.
Введение
Большое количество гидрофизических исследований водных экосистем и трофических цепей направлено на анализ трансформации вещества и энергии в них. Благодаря этим исследованиям в гидробиологии достигнуты большие успехи в понимании механизмов функционирования водных экосистем. В работах ученых Секерци Ю. и Петровского С. рассматривается процесс планктонно-кислородной динамики, учитывающий изменение физических характеристик океана, что, в свою очередь влияет на смертность не только планктона, но и других популяций [1]. Изучением воздействия человека на водную среду занимались Мюррей Дж.В. и Тайлефер М. В их работах исследуется регулирование, а также превышение стока рек и осадков над испарением; вынос загрязняющих веществ в водные объекты вследствие судоходства и стока речных вод; вертикальные плотностные градиенты и ослабление процесса перемешивания вод. Так, в работе [2] исследованы процессы межвидового взаимодействия, описаны потоки веществ, поступающие от одних элементов экосистемы к другим. Вопросами формирования структур в аэробных, субкислородных и анаэробных условиях занимался Еремеев А.П. Исследованием пространственного распределения профилей кислорода и сероводорода в поровых водах, а также расчетом потоков этих компонентов занимались Орехова Н.А. и Коновалов С.К. Изучением процесса самоочищения загрязненного участка реки занимались Давыдов А.С. и Михайлов М.Д., Виноградов М.Е. Они разработали математические модели, учитывающие скорость
поступления кислорода, потоков органических веществ, дефицит и скорость убыли кислорода за счет разложения придонных отложений и дыхания растений [3].
Под эвтрофированием зачастую понимается весь комплекс негативных явлений, приводящий к изменению трофического статуса водоема. При эвтрофировании обычно происходит накопление вещества, приводящее к нарушению функционирования лимнистических экосистем. В результате в водоеме возникает дефицит кислорода, резкие колебания биомассы водорослей, изменение таксономического состава сообществ гидробионтов.
Работа посвящена изучению процесса эвтрофирования мелководного водоема с помощью методов и средств математического моделирования. Разработанная авторским коллективом пространственно-трехмерная математическая модель биологической кинетики численно реализована в виде программно-исследовательского комплекса.
Применение равномерной прямоугольной сетки, которая является широко применяемой и многофункциональной, дает возможность дискретизировать поставленную гидробиологическую задачу [4]. Однако есть вероятность, что данный аспект приведет к образованию погрешности решения задачи на границе расчетной области. Такой недочет возможно значительно уменьшить с помощью применения метода частичной заполненности ячеек, который дает возможность увеличить точность построения задач биологической кинетики береговых систем [5]. Применение схем <кабаре> и <крест>, на основе линейной комбинации [6] дает возможность построить дискретную модель. Такой подход позволяет увеличить точность моделирования, если значения сеточного числа Пекле достаточно большие. На базе математического моделирования существует возможность формирования сценарного подхода для анализа критических явлений эвтрофикации и аноксии береговых систем мелководных водоемов [7].
1. Постановка задачи
Модель процесса эвтрофирования Азовского моря базируется на работах Ма-тишова Г.Г., Тютюнова Ю.В. [8], Ильичева В.Г., Сухинова А.И. [9], Четверушки-на Б.Н. [10], Якушева Е.В. [11], Вайнер Э.Р. [12], посвященных моделированию гидрохимических процессов. При разработке математической модели был сделан упор на то, что при разных сценариях ветрового волнения в прибрежных системах Юга России в придонных слоях увеличивается вероятность образования анаэробных условий. Исследования Трегер П., Ковет Г. показали, что восстановление поверхностного во-донасыщенного ила сопровождается высвобождением ряда химических соединений, а именно: сульфатов, двухвалентного марганца и железа, органических соединений, аммония, силикатов и фосфатов [13, 14]. Модель представлена системой уравнений для каждого с - значения концентрации ¿-й субстанции следующего вида:
-1 = #1(с1,Сз,С8, Сд) - Л1С1, = #2(¿2, С4, Сю) - А2С2,
-3 = 73«2(С4)С2С3 - #з(с1, С3) - А3С3 + В(сз - С3) + /,
04 = А1С1 + Л2С2 + Л5С5 - #4(¿2, С4, С5),
(1)
-5 = #5(¿4, С5, С8) - Л5С5, -6 = (#б(С4, С8) - Лб)Сб, -7 = #г(Сб) - Л7С7,
-8 = (#8(¿2, С7) - Л8)С8 - 78«1(С8)С1,-д = #д(¿1) - ЛдСд,
-10 = #10(¿1, ¿10) + 72«2(^0)С2 - (^1 ¿1 + 77С7)С10,
где и = (и,и,т) - вектор скорости водного потока; /ш91 - скорость гравитационного осаждения г-й компоненты; Д - двумерный оператор Лапласа; - горизонтальный и вертикальный коэффициенты турбулентного обмена соответственно; фг - химико-биологический источник (сток) или член, описывающий агрегирование (слипание-разлипание), если соответствующая компонента является взвесью, индекс г указывает на вид субстанции, г = 1,10. При построении модели параметризованы процессы биогеохимических циклов химических элементов, ответственных за трансформацию аэробных условий в анаэробные.
Функции зависимости скорости роста фитопланктона и бактерий от температуры (Т), солености (Б) (по Леману) и освещенности (I) (по Стилу) имеют вид:
Фт(Т, S, I)
I
I
-ex'p
1 - аг.
opt
(T - Topt)
T
opt
- e
(S — Sopt)
S
opt
I
I
, m
(1, 2, 5}
opt
где I(к) = 10ехр(-вк); к - глубина водоема; 10 - суммарная солнечная радиация на поверхности водоема; в - коэффициент затухания, пропорциональный концентрации взвешенных в воде веществ (фитопланктона и детрита); Т^, Борг - температура и соленость, оптимальные для данного вида гидробионтов; ат, вт - коэффициенты ширины интервала толерантности фитопланктона и бактерий к температуре и солености соответственно.
Расчетная область О представляет собой замкнутый бассейн, ограниченный цилиндрической боковой поверхностью а, невозмущенной поверхностью водоема Е0, дном Ея = Ея(х,у). Е - кусочно-гладкая граница области О, Е = Е0 и Ея и а, 0 < Ь < То. Пусть ип - нормальная составляющая поля скорости водного потока, и -вектор внешней нормали к Е.
Добавим граничные условия:
дСг
Сг = 0 на <т, если 17п < 0; -7— = 0 на <т, если 17п > 0;
дп
до' до'
-7JJ = на Ея, = f(ci) на =
■угог, г Е (3, 4, 6, 7} ; (Ki-a), г = 10;
«о Vi
(2)
.0, г Е (1, 2, 5, 8, 9}
где а0 - коэффициент аэрации; К10 - концентрация растворенного в воде кислорода при насыщении; - коэффициент поглощения г-й субстанции донными отложениями, г = 1,10. К модели добавляются начальные условия вида:
Ci\t=o = Cio(x,y,z), г = 1,10.
(3)
Для определения условий существования и единственности решения задачи (1)-(3) непрерывная модель линеаризована и аналитически исследована методом построения аналога функционала энергии. Результаты сформулированы в виде теорем.
Теорема 1. Пусть Ci(x,y, z,t), фг £ С^(Цг)Г\С(Цг), где Цг = Gx(t0,T0); fa = const > 0; U = (u,v,w — wgi),Ui(z) G C1(G); ci0 G C(G), i = 1,10. Тогда при выполнении неравенств: max {fa, щ} — ^max{|pj|} > 0 для всех i = 1,10, где фг = Pi(cj)ci +
фг, г ф j; фг = Td - Dd ~PiCi , Td = ^ + div(ciUi), Dd = faAd + Jj(^f^), A0 = n2(1/l2x + 1/l"2 + 1/l2z); lx,ly,lz - пространственные максимальные размеры расчетной области; задача вида (1) - (3) имеет решение.
2
2
Теорема 2. Пусть Ci(x,y,z,t), ^ G С2(Щ) Г\С(Щ), = const > 0; U, u^z) G Cl(G), i = 1,10. Тогда при выполнении неравенств: 2^(1/72 + 1 /72) + 2щ/1^ > "фг для всех г = 1,10 (функции -фг определяются источниками загрязняющего вещества (ЗВ)) задача биологической кинетики вида (1) - (3) имеет единственное решение.
2. Метод решения задачи
В работе [15] было проведено исследование алгоритма для решения многомерных задач диффузии-конвекции на основе схемы расщепления на явно-неявную задачи. Также была описана примененимость данного алгоритма для решения трехмерной задачи транспорта взвесей в прибрежных системах. Схема расщепления по пространственным координатам позволяет сократить обмены информацией, которые реализуются между соседними потоками при переходе с одного временного слоя на другой параллельным способом в приграничных узлах подобластей при геометрической декомпозиции трехмерной сеточной области.
Рассмотрим аппроксимацию трехмерной задачи биологической кинетики на примере уравнения диффузии-конвекции-реакции:
c't + (мс)^ + МУ + (wc)Z = (^сХ)Х + с'у )'у + (^cZ )Z + / (4)
с граничными условиями:
c'n (x,y,z,t) = a„c + вп,
где - компоненты вектора скорости водной среды, v - коэффициенты тур-
булентного обмена в горизонтальных и вертикальном направлениях, / - нелинейная функция правой части.
Введем равномерную расчетную сетку по времени:
шТ = {tn = пт, n = 0,Nt, Ntr = T}.
Проведем дискретизацию трехмерной задачи вида (4) на основе схем расщепления по физическим процессам (перенос веществ по каждому из пространственных направлений) :
„га+1/3 _ „га
-+«(о: = (мох (5)
т
cn+2/3 _ cn+1/3
+ v(cn+1/3)y = (M(cn+1/3)y)', (6)
т 'y V 4 ' у/ у
cn+1 _ cn+2/3
т
+ wcZ = (v4 )Z + f (7)
где сп - концентрация химико-биологической субстанции на текущем п-м временном слое Ьп , сп+2/3 - на промежуточном (п + 2/3) временном слое, сп+1 - на последующем (п + 1)-м временном слое, с = (сп+1 _ сп+2/3) /2.
2.1. Линейная комбинация разностных схем <кабаре> и <крест>
Рассмотрим задачу биологической кинетики на примере задачи транспорта веществ, выраженную нестационарным уравнением конвекции:
4 + исХ = 0, (8)
где t е [0, T], x е [0, L]; c (x, 0) = c0 (x); c (0,t) = 0, если u ^ 0 и c (L,t) = 0, если u < 0, u = const.
Введем равномерную по времени и пространству сетку ш = ujh х шт, где ujh = Xi = ih, i = 0, Жц, Nxh = L} , Жц - количество шагов по пространству, L -максимальный размер области вычислений; ojt = |ira| tn = пт, п = 0, Nt, Ntr = Т}, Nt - количество временных слоев, T - верхняя граница по времени.
Запишем классические схемы, используемые для численного решения задачи (8):
- схема «кабаре> для случая u ^ 0 :
cn+1 _ cn cn _ cn—1 cn _ cn
_1 i _ZtiL i vli_(a)
2r 2r + /г '
- схема <крест>:
ra+1 _ cn-i cn _ cn
В результате разложения членов разностных уравнений (9), (10) в ряд Тейлора [16] получим соответственно:
Н--^--b и-;- = — + и-- I
2т 2т h Vdt dx/ г—1/2 (11)
12 / i—1/2 cn+i_ cn—1 cn+1 _ cn—1 _ / 5c 5c \n 1_ r2 2( д3
2T + i + + U?) + (12)
c
ra
г
где г = иг/Л - число Куранта.
Линейная комбинация разностных схем «кабаре> и «крест>, взятых с весами 2/3 и 1/3 соответственно, с учетом (11), (12) запишется в виде:
„и+1 га „га „га-1 и „и-1 га I л „и с „га
„г — „г , оС»-1 - I „» — „» I „»+1 + » - П ^ П
--Ь 2----1-----Ь и—-—-= 0, если и ^ 0;
т 3т 3т 3Л (13)
„и+1 _ и „и — „и-1 и _ „и-1 и _ л и _ г и
„г „г , 0„»+1 „»+1 , г „г , „г+1 4„г 5„г-1 п , п
--Ь 2-1--1---Ь и—1---= 0, если и < 0.
т 3т 3т 3Л
Данная разностная схема имеет локальную погрешность аппроксимации
г (1 — г) иЛ2 (53с/5ж3)и-1/3/6 + О (т2 + Л3) относительно фиктивного узла (г — 1/3, п).
Получили, что порядок погрешности аппроксимации линейной комбинации разностных схем <кабаре> и <крест> с весами 2/3 и 1/3 соответственно выше, чем у исходных схем.
Найдем аналитическое решение уравнения (8), для чего представим его в виде конечной суммы ряда Фурье в комплексной форме:
N
c (x,t)= ^^ Cm (t) exp (jwmx), (14)
m=—N
L
где a; = 7x/L,m - номер моды, Cm (t) = / с (ж, i) exp (jburrix) dx, Cm - комплексная
0
амплитуда m-й моды, j2 = _1.
Подставляя (14) в (8), получим:
N \ ' / N
Ст (Ь) ехр (]штх)\ + и I ^^ Ст (Ь) ехр (/штх) I =0.
ут ( ь ) ех р ( /ш 11 (х / | | и |
\m=-N / + \m=-N
Изменяя порядок операций дифференцирования и суммирования ряда, полагая функцию ехр /штх) линейно независимой, получим:
(Ст (*))/ = - /ишт) Ст (Ь) . (15)
Из (15) получим выражение для Ст (Ь) и подставим его в (9). Таким образом, аналитическое решение уравнения (8) примет вид:
N
с (х,Ь) = ^^ Ст (0) ехр (— /иштЬ) ехр (/штх). (16)
m=-N
Исследуем точность разностной схемы вида (13). Для этого подставим выражение (13) для (хг,Ь) в узле Xi = гН в (14) (рассмотрим случай и ^ 0). Проведем несколько элементарных алгебраических преобразований, учтем, что функции ехр (/штх) линейно независимы, получим:
(17)
Сп+1_Сп Сп _Сп-1 Сп _Сп-1
_-т + -ш_ ехр + -+
т 3т 3т
ехр (/штН) + 4 _ 5 ехр (_/штН)
---з^-^^-- = 0-
Выражение (13) при т ^ 0 примет вид:
, ( ехр (¿штЬ) + 4 - 5 ехр {-¿штЬ) \ п {Ст -2к (2 + ехр {—]штК))-)
С учетом (15) решение, полученное с использованием схемы (13), соответствует решению уравнения с[ = _и*с'х, где и* = и (1 _ а1), а1 = 1 _ ехр (/штН) + 4 _ 5 ехр (_/штН)
-—---—-. Произведя несложные алгебраические преобразова-
2/штН (2 + ехр (_/штН))
ния, получим а\ = ] + О ((штК)А). Из полученного выражения видно, что схема
(5) аппроксимирует конвективный член с третьим порядком точности по пространству.
Легко показать, что при Н ^ 0 из выражения (17) следует, что разностная схема (13) имеет погрешность аппроксимации по времени О (т2). Таким образом, разностная схема (13) для задачи (8) имеет погрешность аппроксимации, равную О (т2 + Н3).
2.2. Разностные схемы для решения уравнений конвекции-диффузии-реакции
Расчетную область впишем в прямоугольник. Построим равномерную по пространству сетку:
шн = {Хг = 0гх, Уу = ]}гУ) гк = г = 0, ] = 0, к = 0, Л^;
Мх Нх 1х , Му Ну 1у , N2 Нх 1Х } ,
х
где Нх, Ну, Н2 - шаги в направлениях пространственных координат, 1х, 1у, - максимальные размеры расчетной области.
Для аппроксимации двумерного уравнения (5) в направлении горизонтальной плоскости хОу на основе линейной комбинации разностных схем <кабаре> и <крест> с учетом заполненности расчетных ячеек будем использовать схему расщепления по пространственным направлениям. Разностная схема для уравнения (5), описывающего перенос в направлении Ох, записывается в виде:
20_2 + 0_0 сп+,к _ Сп,],к , _ сСп,3,к _ сСп-1,],к .
+ 5^-1/2,^92 о 7-~ +
3 т 4-1 ' 3hx
4-t, , rnin (п n \ Ci+i,j,k-ci,j,k , 2A*ci-l,j,k<l?+Axci,j,fcgo _
+ui+i/2,j,k mill [qi,q2)-^--1--3--
y-t^ _ y-t^ y-t^ _ y-t^
_ о,, „ 4+l,j,fc i,j,k 0 i,j,k i— l,j,k
^fa+l/2,j,kQl-Г5--¿H>i-l/2,j,kQ2
(19)
h2 2 h2
I Ь rf I V x
ax ci, j, k + fix hx
x i, j,k x
- \qi ~ q2\fa,j,k-Г-,ui,j,k > 0
n+1/3
2<?i + <?o cij,fc cij,fc p. ci+i,j,fc ci,j,fc . / ч ci,j,fc ci-i,j,fc ,
+ oUi+l/2,j,kQl-777--г ui-l/2,j,k mm Щ1Л2) -7TT--r
3 т ^ 3hx ' 'J' x ' 3hx
OA na n 4- A na n na _ na na _ na
2^x ci+1,j,k q1 + ^x ci,j,kq0_ ci+1,j,k ci,j,k n ci,j,k ci-1,j,k
2^+l/2j,fc9l-^-— - 2fJLi
h 2x
a cn 1 о cn-2/3 _ cn-1
1 1 ax ci,j,k + в x \ n i,j,k i,j,k
~ \Qi ~ q2\Vi,j,k-^-,Ui,j,k < 0,где Axci>j>k =---.
Разностная схема для уравнения (5), описывающего перенос в направлении Оу, записывается аналогично. Здесь и далее коэффициент заполненности расчетной ячейки зависит от номера ячейки, то есть = 9о,г,^,к, Я1 = Ч1,г,з,к, q6 = 4&,г,з,к. Для аппроксимации уравнения (6) в вертикальном направлении используется явная схема:
сп+1 _ сп+2/3 сп+5/6 _ сп+5/6 сп+5/6 _ га+5/6 Чо---ТуГ--г дбЩ,^к-1/2-
2hz 2h
n+5/6 n+5/6 n+5/6 n+5/6
z (20)
_ „ ,, Сг,3,к+1 С^,к , _ Сг,.7> С^,к-1 , — Я5^,к+1/2-^--Н <?6-^--Н
Уравнение (20) решается методом прогонки, что значительно сокращает время расчетов.
3
2.3. Тестирование разработанной разностной схемы
Рассмотрим начально-краевую задачу для нестационарного уравнения теплопроводности:
ct = VcZx + fi, (21)
где t е [0,T] , x е [0,L]; c (0,t) = 0, c (L,t) = 0, c0 (x) = 0(20 _ x) _ 0(10 _ x) , где 9(x) - функция Хэвисайда; ^ = 1 м2/с.
Для численного решения задачи (21) будем использовать схему с весами:
cn+1 _ cn cn+a _ 2cn+CT I cn+a
s—^ = Zc:2 +c+fl, (22)
т h2
где cn+CT = acn+1 + (1 _ a) cn, a е [0,1] - вес схемы.
Необходимое условие устойчивости для явной схемы, полученное на основе метода гармоник, имеет вид неравенства [17]: 7 = тд/Н2 > 0, тогда 478т2 (На/2) ^ 47 ^ 2 ^ 7 ^ 1/2. Таким образом, получили т ^ ттах = Н2/(2д). Несмотря на то, что данная оценка является жестким ограничением для явных разностных схем, на практике шаг по времени необходимо брать еще меньше.
Параметры расчетной сетки задавались следующим образом: шаг по времени находится в диапазоне от 0,001 до 10 с, шаг по пространству Н =1 м, длина интервала по времени Т составила 60 с. На рис. 1 представлен график погрешности решения модельной задачи (21) на основе схемы (22): кривая 1 - явной схемы, 2 - явной схемы с весами (а = 0, 5).
8x1О"4
6х1(Г4
4x1О"4
0
С
Рис. 1. График зависимости погрешности аппроксимации от шага по времени т0: 1 - для явной схемы, 2 - для схемы с весами
Запишем выражение для расчета погрешности вычислений в норме пространства
L2 : Ф = (q — ci)2/ c2 , где ci - точное значение решения задачи (21) в узле i, q -V i i
численное решение, зависящее от величины шага по времени. По горизонтальной оси отложена величина шага по времени т0 , отнесенного к величине ттах (т0 = т/ттах ).
Для того чтобы относительная погрешность явной схемы была равна 0, 01%, необходимо величину т0 брать равной 0,0717, в случае использования предложенной схемы с весами параметр т0 равен 5,1858.
Рассмотрим решение начально-краевой задачи (8) при u = 0, 5 м/си начальном значении c0 (x) = 9 (70 — x) — 9 (60 — x) c0 (x) = 9(20 — x) — 9(10 — x), где 9(x) - функция Хэвисайда. Введем выражение для погрешности решения в норме пространства L : фп = ^ <0пh, = |cn — c (xi, tn)| , где c (xi, tn) - точное решение задачи (8) в узле i
i, cn - численное решение. На рис. 2 представлены графики значений погрешности аппроксимации численного решения задачи (8) с использованием схем «крест», «кабаре», а также их линейной комбинации при различных значениях чисел Куранта (r = |u| т/h, 0, 01 ^ r ^ 1). Шаг по временной сетке т менялся от 0,02 с до 2 с при величине временного интервала 100 с.
Из рис. 2 видно, что предложенная разностная схема (20) точнее описывает решение нестационарного уравнения конвекции при малых значениях чисел Куранта r ^ 0,1. Из расчетов для уравнения теплопроводности известно, что т ^ А • ттах, ттах = h2/2д, А = 0, 01 - отношение шага, при котором достигается необходимая точность расчетов и т остается в приемлемом диапазоне, к шагу, полученному из ограничения на устойчивость явной схемы ттах. Тогда А |u| т/h ^ 0,1, где сеточное число Пекле Pe = |u| h/д ^ 0, 2/А = 20. Отсюда получаем диапазон сеточных чисел Пекле, при которых предложенная разностная схема (20) будет аппроксимировать конвективный член с лучшей точностью, чем классические схемы: 2 < Pe ^ 20, в работе [16] показано, что при Pe ^ 2 более высокую точность дает центрально-разностная схема.
6
4
2
О
0.01 0.05 0.1 0.5 с 1
Рис. 2. График зависимости погрешности численного решения задачи (8) от значений чисел Куранта с использованием разностных схем: 1 - (20); 2 - «кабаре»; 3 - «крест»
3. Результаты численных экспериментов
На базе разработанного программного комплекса, ориентированного на вычислительную систему с распределенной памятью, был проведен численный эксперимент для задачи биологической кинетики (1) - (3). Также был изучен процесс эвтрофика-ции мелководного водоема в летний период, механизм самоочищения мелководного водоема с учетом влияния жизнедеятельности аэробных и анаэробных бактерий. Проведен анализ влияния поступающих в водоем биогенных веществ, пространственного распределения солености, температуры и освещенности на рост и смертность клеток фитопланктона в реальной области сложной формы (прибрежная система - Азовское море). При моделировании учитывалось влияние как абиотических, так и биотических факторов на развитие основных гидробионтов. При создании массивов входных данных использовались экспедиционные данные, электронный атлас Азовского моря [18].
Выполнен вычислительный эксперимент для моделирования превращений соединений серы (сероводорода Н2£, элементарной серы £, тиосульфатов £203, сульфатов £04) в присутствии растворенного кислорода и при его недостатке в фиксированной точке придонной области центрально-восточной части Азовского моря. На рис. 3 представлены различные сценарии превращений веществ (растворенного кислорода 02 и сероводорода Н2£), временной период моделирования - 23 дня, шаг по времени - 100 с. Численный эксперимент показывает, что при концентрации растворенного кислорода ниже 1,1 мг/л он полностью расходуется на окислительные процессы, при этом растет концентрация сероводорода. При увеличении количества гидратирован-ных молекул кислорода в воде концентрация ядовитого сероводорода снижается и постепенно достигает нулевого значения. Чем выше в воде концентрация растворенного кислорода, например, вследствие активного перемешивания вод из-за ветрового воздействия, тем быстрее запускается процесс самоочищения водоема. Результаты эксперимента согласуются с данными, полученными в ходе экспедиций, и литературными источниками.
Рис. 3. Изменение концентраций: 1 - растворенного кислорода сю, 2 - сероводорода с6 при а) недостатке кислорода (начальное значение сю = 0, 3 мг/л); б) достаточном количестве кислорода (начальное значение сю = 2 мг/л)
0.005
0
0.01
0.02 мг/л
0.015
Рис. 4. Результаты численного эксперимента: концентрации а) фитопланктона ; б) детрита с4; в) биогенного вещества с3 (нитратов)
На рис. 4а представлено трехмерное распределение концентрации синезеленых водорослей, преобладающих в Таганрогском заливе, и диатомовых водорослей, распространенных в основной акватории, которые вносят основной вклад в образование детрита (рис. 4б) вследствие экскреции и отмирания, а также минерального питания (нитратов, рис. 4в). Временной интервал моделирования - 50 дней, что соответствует временному периоду от начала вегетации до бурного цветения водорослей в летнее время. Численный эксперимент показал, что наибольшая концентрация детрита наблюдается в Таганрогском заливе, так как в нем продуцируется основная часть биомассы фитопланктона. Взвешенные частицы детрита выносятся течением в основную часть моря, где, оседая на дно, могут вызывать явление аноксии при определенных ветровых и температурных режимах. Считается, что фотосинтез пропорционален скорости роста первичных продуцентов с постоянным для каждой группы коэффициентом ассимиляции. Процесс реаэрации дает возможность учитывать концентрацию насыщения растворенным кислородом, а коэффициент реаэрации есть функция скорости ветра. Основные процессы, расходующие растворенный в воде кислород, - это окисление детрита в толще воды и придонном слое и дыхание живых существ, в том числе фитопланктонных популяций. При моделировании трансформации кислорода учитывалась его изменчивость в соответствии со стехиометрически-ми соотношениями. В описании химико-биологического источника кислорода учтены другие процессы, приводящие к потере кислорода при окислении биогенных веществ. В балансовом соотношении для кислорода учитывался его поток через границу с атмосферой. В свою очередь, этот поток сильно зависит от температуры воды. Эксперименты показали, что рассматриваемый метод является достоверным для классического формирования фитопланктонных популяций в теплый период времени в Азовском море. Также существует возможность исследовать процесс эвтрификации и метод самоочищения вследствие развития бактерий, которые задействованы в процессе разложения детрита и возобновлении сульфата до сероводорода при анаэробном дыхании микроорганизмов, которые содержат соленость, температуру и освещенность, на продукционно-деструкционные процессы и явления фитопланктона.
Заключение
В результате проведенных исследований разработана модификация пространственно-временной гидрофизической модели динамики планктона, которая содержит в себе изменение концентрации фитопланктона, аэробных и анаэробных бактерий, биогенных веществ, с целью изучения гидробиологических процессов мелководного водоема. Применение сценарного подхода в разработанном программном комплексе позволило изучить влияние отсутствия аэрации водной среды и постоянного притока кислорода на процессы эвтрофикации водоема. Представлены и изучены условия, при которых не учитывается устойчивость круговорота веществ в экосистеме, изучен сценарий структуры водоема при неблагоприятном действии природных и техногенных факторов на жизненный цикл аэробных и анаэробных бактерий, принимающих участие в процессе аммонификации детрита.
Схема расщепления исходной трехмерной задачи вида (1) - (3) по геометрическим направлениям на двумерную и одномерную по горизонтальному и вертикальному направлениям соответственно позволяет ученьшить количество обменов пакетами данных между соседними потоками при переходе от одного временного слоя к другому параллельным способом в приграничных узлах подобластей на основе геометрической декомпозиции трехмерной сеточной области по вертикальным плоскостям. Дискретизация разработанной модельной задачи водной экологии была произведена на базе линейной комбинации схем <кабаре> и <крест> с учетом частичной заполненности расчетных ячеек.
Численная реализация поставленной задачи биологической кинетики дает возможность изучить как внутри-, так и межвидовые химические связи между планктонными популяциями прибрежной системы - Азовское море в режиме ограниченного времени, что является актуальным при возникновении экологических ситуаций катастрофического характера, к которым можно отнести эвтрофикацию и процессы экзогенной гипоксии. Установлено, что существенная неоднородность структуры детрита, вызванная различными фракциями природных органических веществ, играет важную роль в регулировании глобальных процессов водной экологии.
Исследование выполнено за счет гранта Российского научного фонда № 22-71-10102, https://rscf.ru/project/22-71-10102/.
Литература
1. Petrovskii, S. Regime Shifts and Ecological Catastrophes in a Model of Plankton-Oxygen Dynamics under the Climate Change / S. Petrovskii, Y. Sekerci, E. Venturino // Journal of Theoretical Biology. - 2017. - V. 424, № 13. - P. 91-109.
2. Комилов, Ф.С. Учет гидро-климатических и физико-химических характеристик экосистемы рыбоводного пруда при ее компьютерном моделировании / Ф.С. Комилов, С.Х. Мирзоев, Ф. Акобирзода // Вестник Таджикского национального университета. Серия: Естественные науки. - 2015. - Т. 156, № 1-1. - С. 19-27.
3. Виноградов, М.Е. Влияние изменений плотности воды на распределение физических, химических и биологических характеристик экосистемы пелагиали Черного моря / М.Е. Виноградов, Ю.Р. Налбандов // Океанология. - 1990. - Т. 30, № 5. - С. 769-777.
4. Самарский, А.А. Численные методы решения задач конвекции-диффузии / А.А. Самарский, П.Н. Вабищевич. - М.: URSS, 2009.
5. Sukhinov, A.I. Numerical Realization of Three-Dimensional Model of Hydrodynamics for Shallow Water Basins on High-Performance System / A.I. Sukhinov, A.E. Chistyakov, E.V. Alekseenko // Mathematical Models and Computer Simulations. - 2011. - V. 23, № 3. -P. 3-21.
6. Головизнин, В.М. Некоторые свойства разностной схемы <кабаре>/ В.М. Головизнин, А.А. Самарский // Математическое моделирование. - 1998. - Т. 1, № 1. - С. 101-116.
7. Fennel, K. The Generation of Phytoplankton Patchiness by Mesoscale Current Patterns / K. Fennel // Ocean Dynamics. - 2001. - V. 52, № 2. - P. 58-70.
8. Tyutyunov, Yu.V. Spatiotemporal Pattern Formation in a Prey-Predator System: the case Study of Short-Term Interactions between Diatom Microalgae and Microcrustaceans / Yu.V. Tyutyunov, A.D. Zagrebneva, A.I. Azovsky // Mathematics. - 2020. - V. 8, № 7. -P. 1065-1080.
9. Sukhinov, A.I. Game-Theoretic Regulations for Control Mechanisms of Sustainable Development for Shallow Water Ecosystems / A.I. Sukhinov, A.E. Chistyakov, G.A. Ugol'nitskii, A.B. Usov, A.V. Nikitina, M.V. Puchkin, I.S. Semenov // Automation and Remote Control. - 2017. - V. 78, № 6. - P. 1059-1071.
10. Четверушкин, Б.Н. Пределы детализации и формулировка моделей уравнения сплошных сред / Б.Н. Четверушкин // Математическое моделирование. - 2012. - Т. 24, № 1. -С. 33-52.
11. Yakushev, E.V. Seasonal Changes in the Hydrochemical Structure of the Black Sea Redox Zone / E.V. Yakushev, O.I. Podymov, V.K. Chasovnikov // Oceanography. - 2005. - V. 18, № 2. - P. 48-55.
12. Weiner, E.R. Application of Environmental Chemistry: a Practical Guide for Environmental Professionals / E.R. Weiner. - Boca Raton; London; New York; Washington: CRC Press, 2000.
13. Treguer, P. Water Column Biogeochemistry Below the Euphotic Zone / P. Treguer, L. Legendre, R. Rivkin, O. Ragueneau, N. Dittert // Ocean Biogeochemistry. -Berlin: Springer-Verlag, 2003. - P. 145-156.
14. Cauwet, G. Determination of Dissolved Organic Carbon and Nitrogen by High Temperature Combustion / G. Cauwet // Methods of seawater analysis. - Weinheim: Willey-VCH Verlag, 1999. - P. 407-117.
15. Сухинов, А.И. Экономичные явно-неявные схемы решения многомерных задач диффузии-конвекции / А.И. Сухинов, А.Е. Чистяков, В.В. Сидорякина, С.В. Проценко // Вычислительная механика сплошных сред. - 2019. - Т. 12, № 4. - С. 435-445.
16. Sukhinov, A. Development and Research of a Modified Upwind Leapfrog Scheme for Solving Transport Problems / A. Sukhinov, A. Chistyakov, I. Kuznetsova, Y. Belova, E. Rahimbaeva // Mathematics. - 2022. - V. 10, № 19. - Article ID: 3564. - 11 p.
17. Samarskii, A.A. Numerical Methods for Solving Convection-Diffusion Problems / A.A. Samarskii, P.N. Vabischevich. - Education, 1998.
18. Экологический атлас Азовского моря. - URL: http://atlas.iaz.ssc-ras.ru/sitemap-ecoatlas.html (дата обращения 20.02.2023)
Юлия Валериевна Белова, кандидат физико-математических наук, кафедра «Математика и информатика:», Донской государственный технический университет (г. Ростов-на-Дону, Российская Федерация), [email protected].
Елена Олеговна Рахимбаева, кафедра «Программное обеспечение вычислительной техники и автоматизированных систем», Донской государственный технический университет (г. Ростов-на-Дону, Российская Федерация), [email protected].
Владимир Николаевич Литвинов, кандидат технических наук, кафедра «Математика и информатика», Донской государственный технический университет (г. Ростов-на-Дону, Российская Федерация), [email protected].
Александр Евгеньевич Чистяков, доктор физико-математических наук, кафедра «Программное обеспечение вычислительной техники и автоматизированных систем», Донской государственный технический университет (г. Ростов-на-Дону, Российская Федерация), [email protected].
Алла Валерьевна Никитина, доктор технических наук, кафедра «Программное обеспечение вычислительной техники и автоматизированных систем», Донской государственный технический университет (г. Ростов-на-Дону, Российская Федерация), nikitina. vm@ gmail. com.
Ася Михайловна Атаян, кафедра «Программное обеспечение вычислительной техники и автоматизированных систем», Донской государственный технический университет (г. Ростов-на-Дону, Российская Федерация), [email protected].
Поступила в редакцию 23 февраля 2023 г.
MSC 92B05 DOI: 10.14529/mmp230202
THE QUALITATIVE REGULARITIES OF THE EUTROPHICATION PROCESS OF A SHALLOW WATER RESEARCH BASED ON A BIOLOGICAL KINETICS MATHEMATICAL MODEL
Yu.V. Belova1, E.O. Rahimbaeva1, V.N. Litvinov1, A.E. Chistyakov1,
A.V. Nikitina1, A.M. Atayan1
1Don State Technical University, Rostov-on-Don, Russian Federation
E-mails: [email protected], [email protected], [email protected],
[email protected], [email protected], [email protected]
The article is devoted to modeling the processes of eutrophication of a shallow water body on a computer system with distributed memory. The proposed mathematical model of biological kinetics is based on a system of non-stationary convection-diffusion-reaction equations with non-linear terms, takes into account the movement of water flow, gravitational settling of impurities, microturbulent diffusion, decomposition of detritus as a result of the activity of aerobic and anaerobic bacteria. The introduction of a nonlinear dependence of the growth rate of phytoplankton and bacteria allows to describe the production-destruction processes in a reservoir, to control their dynamics under conditions of excessive intake of biogenic substances (nitrogen, phosphorus and silicon compounds), sulfur compounds, including hydrogen sulfide and sulfates, under various oxygen distribution modes, detritus, spatial and temporal variability of illumination, salinity and temperature, which corresponds to modern ideas about the functioning of the hydrobiocenosis of a shallow water body. The linearization of the continuous problem is carried out, its discrete analogue is constructed from the linearized model based on the splitting of the original three-dimensional problem into two-dimensional and one-dimensional. To build a discrete two-dimensional model, a linear combination of Upwind and Standard Leapfrog difference schemes was used, considering the partial filling of the calculation cells, which allowed to increase the accuracy of modeling the studied processes and phenomena. The results of diagnostic modeling of the processes of hydrogen sulfide contamination and self-purification of a shallow reservoir are presented.
Keywords: Azov Sea; eutrophication; mathematical modelling; explicit-implicit difference scheme; approximation error; computational domain decomposition; distributed memory computing system.
References
1. Petrovskii S., Sekerci Y., Venturino E. Regime Shifts and Ecological Catastrophes in a Model of Plankton-Oxygen Dynamics under the Climate Change. Journal of Theoretical Biology, 2017, vol. 424, no. 13, pp. 91-109. DOI: 10.1016/j.jtbi.2017.04.018
2. Komilov F.S., Mirzoev S.H., Akobirzoda F. [Accounting for the Hydro-Climatic and Physico-Chemical Characteristics of the Fish Pond Ecosystem in Its Computer Modeling]. Bulletin of the Tajik National University. Series: Natural Sciences, 2015, vol. 156, no. 1-1, pp. 19-27. (in Russian)
3. Vinogradov M.E., Nalbandov Yu.R. Effect of Water Density Changes on the Distribution Of Physical, Chemical and Biological Characteristics in the Pelagic Ecosystem of the Black Sea. Oceanology, 1990, vol. 30, no. 5, pp. 567-573.
4. Samarskij A.A., Vabishhevich P.N Chislennye metody reshenija zadach konvekcii-diffuzii [Numerical Methods for Solving Convection-Diffusion Problems]. Moscow, URSS, 2009. (in Russian)
5. Sukhinov A.I., Chistyakov A.E., Alekseenko E.V. Numerical Realization of Three-Dimensional Model of Hydrodynamics for Shallow Water Basins on High-Performance System. Mathematical Models and Computer Simulations, 2011, vol. 3, no. 5, pp. 562-574. DOI: 10.1134/S2070048211050115
6. Goloviznin V.M., Samarskij A.A. [Some Properties of the Difference Cabaret Scheme]. Matematicheskoe modelirovanie [Mathematical Modelling], 1998, vol. 1, no. 1, pp. 101-116. (in Russian)
7. Fennel K. The Generation of Phytoplankton Patchiness by Mesoscale Current Patterns. Ocean Dynamics, 2001, vol. 52, no. 2, pp. 58-70. DOI: 10.1007/s10236-001-0007-y
8. Tyutyunov Yu.V., Zagrebneva A.D., Azovsky A.I. Spatiotemporal Pattern Formation in a Prey-Predator System: the Case Study of Short-Term Interactions between Diatom Microalgae and Microcrustaceans. Mathematics, 2020, vol. 8, no. 7, pp. 1065-1080. DOI: 10.3390/math8071065
9. Sukhinov A.I., Chistyakov A.E., Ugol'nitskii G.A., Usov A.B., Nikitina A.V., Puchkin M.V., Semenov I.S. Game-Theoretic Regulations for Control Mechanisms of Sustainable Development for Shallow Water Ecosystems. Automation and Remote Control, 2017, vol. 78, no. 6, pp. 1059-1071. DOI: 10.1134/S0005117917060078
10. Chetverushkin B.N. Limits of Detail and Formulation of Continuum Equation Models. Mathematical Models and Computer Simulations, 2013, vol. 5, no. 3, pp. 266-279. DOI: 10.1134/S2070048213030034
11. Yakushev E.V., Podymov O.I., Chasovnikov V.K. Seasonal Changes in the Hydrochemical Structure of the Black Sea Redox Zone. Oceanography, 2005, vol. 18, no. 2, pp. 48-55. DOI: 10.5670/oceanog.2005.41
12. Weiner E.R. Application of Environmental Chemistry: a Practical Guide for Environmental Professionals. Boca Raton, London, New York, Washington, CRC Press, 2000.
13. Treguer P., Legendre L., Rivkin R.T., Ragueneau O., Dittert N. Water Column Biogeochemistry below the Euphotic Zone. Ocean Biogeochemistry. Berlin, Springer-Verlag, 2003, pp. 145-156.
14. Cauwet G. Determination of dissolved Organic Carbon and Nitrogen by High Temperature Combustion. Methods of Seawater Analysis. Weinheim, Willey-VCH Verlag, 1999, pp. 407-117.
15. Sukhinov A.I., Chistyakov A.E., Sidoryakina V.V., Protsenko S.V. Economical Explicit-Implicit Schemes for Solving Multidimensional Diffusion-Convection Problems. Journal of Applied Mechanics and Technical Physics, 2020, vol. 61, no. 7, pp. 1257-1267. DOI: 10.1134/S0021894420070159
16. Sukhinov A., Chistyakov A., Kuznetsova I., Belova Y., Rahimbaeva E. Development and Research of a Modified Upwind Leapfrog Scheme for Solving Transport Problems. Mathematics, 2022, vol. 10, no. 19, article ID: 3564. DOI: 10.3390/math10193564
17. Samarskii A.A., Vabischevich P.N. Numerical Methods for Solving Convection-Diffusion Problems. Education, 1998.
18. Ecological Atlas of the Azov Sea. Available at: http://atlas.iaz.ssc-ras.ru/sitemap-ecoatlas.html (accessed on 20.02.2023) (in Russian)
Received February 23, 2023
