Введение

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

В работе [1] представлены результаты исследования влияния марикультуры мидии Mytilus galloprovincialis на морскую среду, полученные при помощи математического моделирования. Авторы использовали основанную на теории динамического бюджета энергии популяционную модель, дополнив ее натурными данными по содержанию сестона и химическому составу тканей мидий. Данные были получены в 2006–2007 гг. на трех мидийных фермах, расположенных на территории Италии и Словении в Адриатическом море. Динамику потоков углерода, азота и фосфора оценивали в течение 10-месячного срока выращивания. Расчеты показали, что количество азота и фосфора, выделяемое мидиями в процессе жизнедеятельности в виде фекалий и псевдофекалий, в два и пять раз больше зафиксированного в структуре моллюсков соответственно. Полученные результаты свидетельствуют о необходимости комплексного подхода к вопросам размещения марикультуры, учитывающего влияние морского хозяйства на окружающую среду.

Исследование [2] посвящено оценке влияния стока р. Апалачикола (Флорида, США) на популяцию гигантской устрицы Magallana gigas. Авторы использовали трехмерную гидродинамическую модель циркуляции и натурные данные, касающиеся жизненного цикла устрицы, для построения стохастических моделей. Результатом работы является определение ряда факторов среды, влияющих на динамику роста устриц. Было установлено, что главным фактором является гидрологический режим р. Апалачикола. Несмотря на детальное изучение процессов влияния морской среды на популяцию устриц, в исследовании остался неизученным вопрос о влиянии самой устричной фермы на экосистему акватории. Добиться более ясного понимания процессов, происходящих между марикультурой и средой, можно при объединении гидродинамической модели с биологической.

Работа [3] посвящена исследованию влияния устричной фермы на экосистему приливной зоны зал. Маренн-Олерон, Франция. Основное внимание уделено структуре пищевой цепи залива. Для исследования влияния устричной фермы на залив используется углеродная модель пищевой цепи. Анализ результатов показал, что присутствие устриц в экосистеме изменяет связи между бентосными и пелагическими организмами, что делает бентосных животных главными потребителями органического углерода. Увеличение площади устричной фермы в два раза привело к возрастанию вторичной продукции, которая стала пищей для молоди нектона. При помощи математического моделирования авторам удалось показать, как наличие марикультуры двустворчатых моллюсков в экосистеме может повлиять на видовой состав фауны приливной зоны. В процессе роста устрицы выделяют достаточно большие объемы аммония и фосфатов, однако в указанной работе не исследовалось прямое влияние этих соединений на экосистему.

Эффективность выращивания марикультуры часто оценивают при помощи показателя экологической емкости акватории. В работах [4–6] приведен ряд моделей, используемых для определения экологической емкости акватории, где расположены объекты марикультуры. В [4] представлена двухмерная физико-биогеохимическая модель, реализованная для зал. Сунго, Китай. Результаты [4, 7] и других работ были использованы при размещении морского хозяйства, выращивающего поликультуру в зал. Сунго. Несмотря на активную эксплуатацию залива в течение 60 лет, его экосистемы по-прежнему здоровы благодаря реализации научных рекомендаций по ведению морских хозяйств [8].

Анализ литературы показал широкое применение комплексных экологических моделей для исследования экосистем, включающих в себя марикультуру. Работы [1, 3, 6, 8] свидетельствуют о возможности применения таких моделей для оценки потоков нутриентов и углерода в экосистемах. Также комплексные экологические модели позволяют определять влияние гидрохимического режима на продуктивность культивируемых видов [2, 7, 9], прогнозировать экономическую прибыль [6, 8]. Анализ рассмотренных ранее моделей показал, что все они имеют географическую привязку и результаты моделирования существенно зависят от локальных условий. Перспективным является направление, при котором в качестве ядра биологического блока модели используется индивидуум-ориентированная модель культивируемого вида.

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

Материалы и методы

Предмет и объект исследования. Европейская устрица Ostrea edulis в прошлом веке являлась промысловым видом, однако ввиду большого сокращения численности утратила свое промысловое значение . Культивирование O. edulis представляет интерес в связи с актуальностью проблемы восстановления численности дикорастущей устрицы и выращивания в качестве марикультуры [10]. Объектом имитационного моделирования был выбран зал. Донузлав, в котором находится ряд действующих морских хозяйств по выращиванию мидии и устрицы [11]. Наличие натурных данных о динамике морфометрических характеристик плоской устрицы в этом районе ³ способствовало выбору зал. Донузлав в качестве района размещения модельной плантации.

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

Объект ООМ определяется как совокупность идентичных организмов и их ближайшего жизненного пространства (БЖП), представленного, например, скоплением фитопланктона, небольшим участком дна, занятым макроводорослями определенного вида, коллектором мидий или садком устриц; БЖП – область, в которой протекают обменные процессы организмов и среды. Геометрически модельные объекты представляются в виде цилиндров, размеры которых соответствует реальным физическим размерам выделенного БЖП для совокупности морских организмов. Предполагается однородное распределение биомассы внутри объекта. Для цилиндра задается высота и радиус основания. Начальное положение объектов в пространстве модели можно выбирать или распределять случайным образом. При моделировании морских ферм объекты, которые имитируют выращиваемую марикультуру, располагаются в соответствии со схемой фермы. Объекты фитопланктона и фитобентоса распределяются случайным образом.

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

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

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

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

ООМ морского хозяйства. На описанных выше принципах ранее была разработана двухмерная химико-биологическая объектно-ориентированная модель морского хозяйства [12]. Для имитации функционирования устричной фермы и процессов ее взаимодействия с морской средой в описание класса «Моллюски» был добавлен новый метод, разработанный на основе одномерной модели динамического энергетического баланса (ДЭБ) европейской устрицы [13], которая математически описывает процессы фильтрации, питания, дыхания, экскреции, роста и размножения устрицы. Результаты верификации модели показали хорошее соответствие значений расчетных величин линейного и весового роста устриц измеренным в ходе натурного эксперимента в зал. Донузлав в 2001–2003 гг. Рассмотрим более подробно структуру классов ООМ (рис. 1).

Р и с.  1. Схема классов химико-биологической объектно-ориентированной модели морской экосистемы. Черные стрелки обозначают наследование, голубые – включение множества объектов одного класса в свойства другого, более старшего класса

F i g.  1. Class diagram of the chemical-biological object-oriented model of marine ecosystem. Black arrows indicate heredity, blue ones – inclusion of multiple objects of one class into the properties of another higher class

«Расчетная сетка» – класс с координатами расчетной области и регулярной сеткой, в узлах которой задаются значения гидрофизических и гидрохимических переменных. Свойства класса: глубина водоема H, шаг сетки в dx (зональном) и dy (меридиональном) направлениях, размеры области Ndx ´ Mdy, границы водоема в пространстве, шаг по времени dt.

Класс «Гидродинамика» – наследник класса «Расчетная сетка». Основным методом класса является упрощенная интегральная модель верхнего слоя Черного моря [14, 15], которая описывает пространственно-временну́ю динамику температуры (T⁰, ℃) и глубины (h, м) верхнего перемешанного слоя, а также температуры нижнего слоя (T^h, ℃). Указанные слои предполагаются однородными по плотности. Учитывая, что масштаб пространственной изменчивости полей ветра и температуры воздуха над морем (особенно в летний период) превышает размер расчетной области, используется одномерный вариант модели:

h∂T0∂t=Г0-Гh-0 ,

Г0=a(T0-Ta) ,

∂h∂t=1T0-Th(Гh-0-Гh+0) ,

G-D=0.5αghρw(Г0+Гh-0) ,                                   (1)

G=CgV*3 , D=ρV*3(fhL1+a2hL|N|) ,

fhL1=a1hL1  при h≤Cga1L1Cg при h>Cga1L1 ,

где Г⁰ – поток тепла через поверхность; Г^h-0 и Г^h+0 – потоки тепла на нижней границе перемешанного слоя и верхней границе нижнего слоя соответственно; Т_а – температура воздуха (℃); a – коэффициент термического расширения воды; r_w – плотность воды (кг/м³); G, D – интегральные генерация и диссипация энергии турбулентности; V*=τ0ρw  – динамическая скорость трения в воде; t₀ – напряжение трения ветра на поверхности; L – масштаб длины Монина – Обухова (м); L₁ – масштаб толщины пограничного слоя (м); N = L/L₁ – безразмерный параметр стратификации Казанского – Монина; C_g, а, а₁, а₂ – эмпирические коэффициенты.

Уравнение (1) определяет условие переключения между двумя режимами динамики верхнего слоя – вовлечения и антивовлечения части водной массы в верхний слой. Вовлечение происходит при G – D – 0,5r_waghГ⁰ > 0 (характеристики верхнего перемешанного слоя меняются, нижнего – нет). В обратном случае происходит антивовлечение (меняются характеристики нижнего слоя, верхний остается неизменным). Однако в результате газо- и теплообмена через поверхность температура и содержание кислорода в верхнем слое могут измениться. Таким образом, двухслойная модель с переменным во времени верхним слоем вносит дополнительное слагаемое в уравнения для гидрохимических переменных верхнего слоя [16]

∂Yi0∂t+U∂Yi0∂x+V∂Yi0∂y=Kx∂2Yi0∂x2+∂2Yi0∂y2+Gi0+Ai0+Гh-0-Гh+0(T0-Th)h(Yih-Yi0)    (2)

и нижнего слоя

∂Yih∂t+U∂Yih∂x+V∂Yih∂y=Kx∂2Yih∂x2+∂2Yih∂y2+Gih+Aih+Гh-0-Гh+0(T0-Th)(H-h)(Yih-Yi0) , (3)

где Y_i⁰ и Y_i^h – переменные блока «Гидрохимия» в верхнем и нижнем слоях (мкМ); U, V – средние в слое скорости течения (м/с); G_i, A_i – функции источников (стоков), определяемые химическими и биологическими процессами соответственно; верхние индексы «0» и «h» определяют принадлежность члена верхнему или нижнему слою; H – глубина бассейна (м); К_х – коэффициент горизонтальной диффузии.

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

– одномерная химико-биологическая модель, совместно описывающая циклы азота и серы в аэробно-анаэробных условиях и разработанная для моделирования процессов в редокс-зоне Черного моря (обеспечивает задание функций G_i в уравнениях (2) и (3)). Уравнения химических реакций и значения параметров, использованных при расчетах, представлены в работе [17];

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

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

Класс «Фитопланктон» является наследником классов «Расчетная сетка» и «Вид водорослей». В настоящем варианте модели было использовано представление фитопланктона в виде пассивной примеси, что упрощает модель и экономит время счета при допустимой ошибке [18]. Основной метод этого класса – модель эволюции поля фитопланктона под воздействием внешних условий и взаимодействия с объектами марикультуры моллюсков. Она является двухмерным вариантом ранее разработанной и апробированной одномерной модели [19]:

∂B∂t+U∂B∂x+V∂B∂y=Kx∂2B∂x2+∂2B∂y2+Pr(1-kr)-ke-mB-Ai,       (4)

       ∂QN∂t+U∂QN∂x+V∂QN∂y=Kx∂2QN∂x2+∂2QN∂y2+VNO3max1-fQN[NO3]KNO3+[NO3]+       + VNH4max1-fQN[NH4]KNH4+[NH4]- PrQN ,                                                        (5)

∂QP∂t+U∂QP∂x+V∂QP∂y=Kx∂2QP∂x2+∂2QP∂y2+VPmax1-fQP[P]KP+[P]- PrQP ,   (6)

Pr=μmin⁡(fQN,fQP) ,

где B – биомасса фитопланктона (мг/м³); t – время (сут); U, V – средние по вертикали скорости адвекции (м/с); m – максимальная удельная скорость валовой продукции фитопланктона (1/сут); P_r – реальная удельная скорость валовой продукции (1/сут); k_r – коэффициент, определяющий затраты энергии на дыхание; k_e – коэффициент выделения органического вещества (экссудации); m – скорость отмирания клеток фитопланктона (1/сут); Q_N и Q_P – внутреннее содержание азота и фосфора соответственно в клетках фитопланктона (пмоль/кл.); Q_N^max, Q_N^min, Q_P^max, Q_P^min – максимальное и минимальное содержание азота и фосфора в клетках фитопланктона (пмоль/кл.); [NO₃], [NH₄], [P] – концентрация нитратов, аммония и фосфора в воде (мкМ); VNO3 , VNH4 , VP  – скорости изъятия нитратов, аммония и фосфора из воды (пмоль/(кл. сут)); VNO3max , VNH4max , VPmax  – максимальные скорости изъятия нутриентов (пмоль/(кл. сут); KP , KNO3 , KNH4  – константы полунасыщения (мкМ); Ai – сумма воздействий на поле фитопланктона со стороны объектов устричной фермы. Число клеток фитопланктона в единице объема легко вычислить, зная содержание углерода в клетке Q_C, поскольку эта величина, в отличие от Q_N и Q_P, полагается постоянной. Для диатомовых Q_C = 12 пмоль/кл., B мг С/м³ = 0,083 ммоль С/м³, тогда число клеток K=0,08312BQC109=6,9∙106B(кл/м3) .

Использованные в работе числовые параметры фитопланктона вида Thalassiosira pseudo представлены в табл. 1.

Т а б л и ц а  1

T a b l e  1

Значения параметров фитопланктона

Numerical values of phytoplankton parameters

Свойствами класса «Биологический объект» являются координаты объекта в пределах области моделирования, вертикальный и горизонтальный размер, масса объекта, интервал возможных вариаций плотности и размеров объекта.

Свойства класса «Вид моллюсков» включают все необходимые параметры для описания энергетического баланса различных видов моллюсков. Класс «Моллюски» является дочерним по отношению к классам «Вид моллюсков» и «Биологический объект» и, таким образом, содержит полное описание носителя (коллектор или плот) с группой культивируемых моллюсков, включая набор параметров модели ДЭБ, которая описывает их рост на носителе и взаимодействие с окружающей средой. Основным методом класса «Моллюски» является модель ДЭБ соответствующего вида моллюсков [13]:

KddWddt=Psom+Pgen ,

P=Psom+Pgen+Psh=AeFKcC-R-Ex ,

где K_d – калорийность тканей моллюска (кал/(г сух. вес)); W_d – сухой вес мягких тканей (г); P_som – энергетические затраты на рост соматических тканей (кал/ч); P_gen – энергетические затраты на рост генеративных тканей (кал/ч); P_sh – энергетические затраты на рост раковины (кал/ч); R – энергия, затрачиваемая на дыхание моллюска (кал/ч); E_x – экскретируемая энергия (кал/ч); F – скорость фильтрации (л/ч); K_c – калорийность кормовой взвеси (кал/мг); С – концентрация кормовой взвеси (мг/л); A_e – коэффициент ассимиляции пищи.

В настоящий вариант экологической модели включены два варианта ДЭБ: для черноморской мидии Mytilus galloprovincialis и плоской устрицы Ostrea edulis. Таким образом, модель позволяет имитировать функционирование нескольких плантаций разных моллюсков в одном районе. Подробные модели ДЭБ, использованные для описания роста устриц и мидий, а также аппроксимации процессов жизнедеятельности, приведены в [13, 20].

Множество объектов класса «Моллюски» входит в состав класса «Марикультура». Этому классу присущи общие и интегральные характеристики морской плантации: ее конфигурация в пространстве, длина, объем, расположение коллекторов относительно поверхности, общая масса выращиваемых моллюсков, их продукционные характеристики, количество потребленных и выделенных веществ.

Аналогично объект класса «Моллюски» может описывать группу идентичных особей двустворчатых моллюсков, которая занимает фиксированную площадь дна с известными размерами и координатами. В этом случае множество таких объектов будет составлять массив «Зообентос».

Класс «Планктонные животные» является наследником классов «Биологический объект» и «Вид зоопланктона». Множество таких объектов («Зоопланктон» на рис. 1) вместе с массивами «Фитобентос» и «Зообентос» включается в класс «Биота», наследник класса «Расчетная сетка». В свойствах этого класса фигурируют поля всех биологических компонентов экосистемы на двухмерной сетке, а также средние или интегральные характеристики, необходимые для анализа поведения системы. Методы этого класса позволяют рассчитать все перечисленные характеристики и преобразовать представления в виде объектов в привычные двухмерные поля, удобные для анализа и визуализации.

Класс «Макрофиты» является наследником классов «Вид водорослей» и «Биологический объект» и, соответственно, в свойствах объектов этого класса содержится вся необходимая информация для расчета продукционных и метаболических процессов макроводорослей определенного вида. Основным методом класса «Макрофиты» является модель ДЭБ, разработанная в [21, 22]. Для макрофитов и микроводорослей используется практически одна и та же модель с некоторыми отличиями, обусловленными особенностями имеющихся натурных данных, на основе которых построены следующие уравнения:

dBmdt=(Pr- ke-m)Bm,

dQNdt=VNO3[NO3]KNO3+[NO3]+VNH4[NH4]KNH4+[NH4]- PrQN,

dQPdt=VPPKP+P-PrQP,

PE=PmaxtanhαEPmax-RdfT,

E=E0exp-βz,

Pr=P(E)kDWminfQN,fQPCw,

fQN=QN-QNminQNmax-QNmin,

fQP=QP-QPminQPmax-QPmin,

VNO3=VNO3max1-fQN,

VNH4=VNH4max1-fQN,

VP=VPmax(1-fQP) ,

где E₀ (мкмоль фотон/м²/с) – фотосинтетически активная радиация (ФАР) на поверхности воды, E – на глубине z (м); B_m – биомасса объекта макрофитов (г сух. вес/м³); P_r – удельная скорость роста (1/сут); P(E) – скорость фотосинтеза (мг O₂/г сух. вес/ч); P_max – максимальная скорость фотосинтеза; α – наклон PЕ-кривой при малых значениях светового потока; R_d (мг О₂/г сух. вес/ч) – темновое дыхание (характеризует количество кислорода, затраченного на дыхание в отсутствие потока ФАР); Q_N и Q_P (мкмоль/г сух. вес) – концентрации фосфора и азота в тканях водоросли; Q_N^max, Q_N^min, Q_P^max, Q_P^min (мкмоль/г сух. вес) – минимальные и максимальные концентрации фосфора и азота в тканях водоросли; VNO3 , VNH4 , VP  – скорости изъятия нитратов, аммония и фосфора из воды (мкмоль/г сух. вес/ч); VNO3max , VNH4max , VPmax  – максимальные скорости изъятия нутриентов; KP , KNO3 , KNH4  – константы полунасыщения (мкМ); m – коэффициент отмирания; k_e – скорость экссудации; k_DW – коэффициент пересчета мг О₂ в мг С.