Object-Based Modeling of Marine Phytoplankton and Seaweeds

Abstract

The aim of this work is to simulate the dissolved oxygen deficiency in the coastal zone that sometimes occurs during the summer water stagnation. We consider the main components of the marine ecosystem that play a major role in such processes—concentrations of nitrogen, phosphorus and sulfur compounds in water, dissolved and particular organic matter, oxygen, biomass of phytoplankton and macroalgae. We use the object-based modeling technique to simulate the spatio-temporal variability of the ecosystem in a 2D domain. In comparison with the traditional approach, it gives several advantages: more precise parametrizations of the biological components’ functionality; higher spatial resolution; possibility to account for the individual variability of hydrobionts; easy inclusion of an arbitrary number of species in the model. Our model included three species of phytoplankton and seven macroalgae. Individual-based modules control their functionality. Species of phytoplankton and seaweeds chosen for simulations are typical for the coastal zone of Crimea. In the simulations, we study the contribution of micro- and macroalgae to the processes of self-purification of a semi-enclosed basin in case of a sharp increase in nutrient concentration in water.

1. Introduction

Coastal marine ecosystems are often exposed to oxygen deficiency, especially in summer under conditions of the stable water column stratification. The reason for this is the imbalance between the sources and sinks of dissolved oxygen (DO). The sources of DO in coastal waters are macro and microalgae production, inflow at the air–sea border, horizontal and vertical exchange with the waters rich in DO. The sinks are respiration, costs on decomposing of organic matter in water and sediments, and oxidation of reduced forms of nitrogen and sulfur. Since the 1960s, hypoxia has become an increasingly acute environmental problem. The number of regions with “dead zones” grew exponentially, doubling about once every 10 years. Now, such zones are noted in 400 coastal ecosystems, and their area exceeds 245 thousand square kilometers [1]. The worst situation is observed in continental seas with limited water exchange with the ocean. These include, for example, large sections of the Baltic Sea [2], the Kattegat Strait, the Black Sea [3,4], the Gulf of Mexico [5,6], river estuaries of the East China Sea [7] and South China Sea [8].

High nutrients loads (eutrophication) lead to large production of microalgae organic matter, and as a result, great demand of DO for its decomposition. However, not all areas receiving high amount of nutrients experience oxygen deficiency. Hydrophysical processes can also induce hypoxia in coastal waters. The formation of low DO areas in the bottom layer is a complicated process caused by a number of physical and biochemical factors. Ecological modeling is a good tool for studying the various processes leading to hypoxia separately and highlighting the main ones. A comprehensive review of modeling approaches to investigation of eutrophication in marine coastal areas has been recently published by [9]. The authors analyzed 291 articles devoted to different aspects of modeling aspects: the aims of the modeling approach (statistical as well as deterministic), target sites and periods covered by the models, modeling tools, space and time resolution, strengths and weaknesses of the models. They conclude that the future models should take into account the long-term changes in benthic nutrient stocks and the possible induced changes in biodiversity. The latter could be achieved moving from a single “many-species” variable to multiple species variables, each with specific ecophysiological parameters. The introduction of more complex and adaptive descriptions of the trophic web could improve the description of the eutrophication process, including top-down control factors. In the review of the modeling studies of eutrophication in lakes [10], the authors came to similar conclusions. One of the proposed directions of development was to merge different modeling approaches, for example, empirical models built by means of the neural network, and physical dynamical models or super individual-based model, and coupled hydrodynamic-ecological model [11,12].

Most of contemporary DO dynamic models are based on advanced hydrodynamic 3D models coupled with ecological blocks. In [6], a coupled physical–biogeochemical model was used to describe the processes at the Louisiana continental shelf. The authors applied a nitrogen cycle model and a simple model for phytoplankton growth [13]. The state variables include phyto- and zooplankton biomass, two forms of inorganic nitrogen, phosphate, two pools of detritus representing large, fast-sinking particles and small, suspended particles. The simulations proved that the key controlling processes driving oxygen dynamics at the Louisiana shelf are physical transport processes and sediment respiration. In summer, stratification isolates the autotrophic surface from the lower heterotrophic waters. This isolation results in the significant outgassing of oxygen across the air–sea interface, increasing the risk of hypoxia.

In [14], a coupled physical-hydrochemical model was developed to simulate the seasonal cycle of DO in Chesapeake Bay and investigate the processes causing summer hypoxia in the estuary. The equation for DO includes advection and diffusion terms, phytoplankton production, water column and sediment respiration. Biological and chemical terms were parameterized by empirical formulas based on observation data. The budget analysis of DO for the lower layer of the water column showed that the wind speed is a more important factor compared to the river run-off. The contribution of the water column respiration to the total oxygen consumption was almost three times greater than sediment respiration.

Wang et al. [7] used a mixing model to distinguish the effects of physical mixing vs. biological processes on the DO dynamics. The difference between the calculated concentration and the measured value was accounted as caused by biological processes. Negative and positive differences correspond to biological uptake and respiratory release of the element. The calculations made it evident that the hypoxia occurred mainly due to a surface diatom bloom and the associated downward fluxes of organic carbon. More than 70% of the organic carbon decomposed within the bottom layer was associated with this bloom. Hypoxia and acidification developed for 6 days after the bloom driven by typhoon storms.

Macroalgae are a very important component of coastal ecosystems due to their photosynthetic activity, and it should be taken into account when simulating the DO dynamics, especially in shallow, photic ecosystems. Notwithstanding this, simulation studies with an included macrophytes submodule constitute only a small part of the growing amount of papers devoted to eutrophication processes. For example, [15] implemented a module of macroalgae dynamics within a water-quality model developed to investigate the reasons for a change in the ecosystem of the Venice lagoon, followed by sharp collapses in population and resultant releases of hydrogen sulphide. The simulations demonstrated that Ulva, the most dominant species, succeeded in the competition with the phytoplankton provided that PAR levels at the bottom were sufficient for growth. The authors emphasize that the model can give useful qualitative estimates of possible long-term changes of the ecosystem, evolving under different scenarios of environment dynamics.

In [16], Ulva lactuca and Gracilaria tikvahiae modules were developed and included into the ecosystem model of a shallow estuary as a part. They predict the biomass dynamics in units of carbon, nitrogen, and phosphorus. Primary production, respiration, grazing, decay, and physical exchange are empirical functions estimated on the local observation data. Simulations showed a strong competition for nutrients between phytoplankton and macroalgae, greater than that for light.

In [17], seagrass and macroalgae modules were implemented in the European Regional Seas Ecosystem Model to simulate the seasonal and interannual trends in nutrient concentrations and macrophyte biomasses in a coastal lagoon. The one-dimensional equations for macroalgae biomass were based on a balance between uptake, release and structural costs of the nutrients. The sensitivity analysis of macroalgae and seagrass to PAR and nitrogen concentration demonstrated high dependence of seagrass on light intensity in comparison with greater influence of nutrients concentration on macroalgae dynamics.

Complex ecosystem models are required to simulate the dynamics of biological and hydrochemical components in coastal areas where bottom macroalgae contribute to primary production and compete with phytoplanktonic community for nutrients and light. The target of this article is to propose a new tool for implementation of macroalgae and phytoplankton components to the two-layer ecological model describing dynamics of dissolved gases, nitrogen, phosphorus, and sulfur compounds. The model can reproduce oxic and anoxic conditions in the bottom layer and transitions between them. Therefore, it can be applied to simulate hypoxia occurrence and reconstruction of the oxidizing conditions.

2. Materials and Methods

First, consider a brief introduction into the object-based modeling (OBM) technique used for the development of the proposed ecological model on the Matlab platform. OBM allows a modeler to group the specific characteristics (properties) of similar entities with the dynamic procedures describing their behavior (methods) into a single definition named class. All that is needed to calculate the future state of these entities is defined in this class. An object is an instance of a class. At the beginning of the program run, the objects are created within a modeling domain (2D or 3D) in accordance with their class definition. The properties of each object represent their state, and the changes of this state depend on the specific location of this object and corresponding environmental influences.

Classes can be organized in the hierarchical structure. Common data and functionality should be placed in a superclass (parent-class), which we then use to derive subclasses (child-class). The subclasses inherit the data and functionality of the superclass and can add some other properties or functions that are unique to their particular purposes. A subclass can have two or more parents inheriting their properties and methods. In our model, the subclass “Algae” has two parents: “Bio-object” defining the common properties and methods of all biological objects, and “Species” determining the species-specific parameters of the corresponding subclass. Physical interpretation of the object is a water column, the nearest living space of a group of identical plants or animals where all the processes of exchange with the environment occur. Therefore, properties of an object include its location, radius of the horizontal projection of a water column, its height, total mass, and biomass as well as individual characteristics of the identical entities inside it. The dynamic models for the micro and macroalgae subclasses describing their evolution in the changing environment are the methods of the appropriate subclasses. These are individual-based models (IBMs) [18,19]. In comparison with the grid representation of biological components in the spatial ecological models, this approach provides the following advantages:

• IBMs using parametrizations based on laboratory experiments on defining the physiological functions of marine hydrobionts promises more precise simulation of the population dynamics than models based on grid representation of plankton fields.
• Higher spatial resolution within the aggregation of any hydrobionts due to adjustable size of the objects.
• An opportunity to account for the individual variability of organisms and plants, which we parameterize in a set of empirical coefficients used in the IBMs.
• Natural and easy inclusion of additional biological components (in particular, aquaculture objects) into a spatial ecological model.

Figure 1 shows the block diagram of the proposed model and the general scheme of classes used for the description of the biological components (hydrobiological module HB). According to our conception, all biological components (discrete entities) are described as sets of objects distributed in the 2D spatial domain (Figure A1), whereas continuous hydrophysical and hydrochemical variables are presented as fields on a regular grid. We consider organismal and population levels of the ecosystem biological components. The first one is described in the frame of the object properties. The second is a classical representation characterized by average biomass and abundance, taxonomic structure, spatial variability, etc. Variables at this level are calculated by averaging the objects’ properties. Dynamics of the chemical compounds (nitrates, nitrites, ammonia, dissolved oxygen, dissolved organic matter (DOM) and particular organic matter (POM) in nitrogen units, hydrogen sulfide, thiosulfate, sulfur, sulfate, and phosphate) are described within the hydrochemistry module (HC). The two-layer integral hydrodynamics module (HD) of the ecological model controls current velocities and water temperature. The vertical structure of the calculated domain contains the upper mixed and bottom layers divided by a seasonal thermocline. During spring and summer warming, the upper layer is deepening and the structure becomes vertically homogeneous. In the periods of summer stagnation caused by low wind speed, a new mixed layer can appear. During low wind periods, mixing cannot reach the bottom layer and under certain conditions, oxygen deficiency can occur. The mathematical description of this process is included in the HC module. Oxygen gets inside through the exchange with the atmosphere and production of the autotrophs during the photosynthesis. The sinks of oxygen are oxidation costs, respiration of plants (in the dark), and an outflow to the atmosphere.

Hydrophysical state variables of the HD module provide external inputs to the HC and the HB module shown in Figure 1. Chemical and biological components of the ecosystem interact by the exchange of matter; therefore, corresponding state variables are connected through the special procedures in the HB module. We calculated available resources for a biological object in the proportion to the intersection area of a horizontal projection of an object and a grid cell. Concentrations of the chemical compounds change in the same proportion due to the uptake and release of these compounds by the macro and microalgae objects. The objects can overlap each other and form an aggregation. Therefore, they compete for light and nutrients. Light intensity for macroalgae objects located at the bottom depends on the biomass of the planktonic objects at the same location. At each time step, we resolve the objects at a random order. If an aggregation is so tight that there is deficit of resources for growth, the production rate of some objects can be zero or negative. The objects with the biomass below the threshold are eliminated from the system. Under better trophic conditions, the objects can bisect forming new ones that are included in the system.

2.1. Microalgae Objects

Properties of the microalgae object define its position in the calculated domain, size, biomass, species, number of cells, intracellular content of nutrients, and some other variables. Methods of the object determine how it can interact with other objects and the environment. The main method describes the process of growth, which is supplemented by the procedures for estimation of available resources, consumption and release of matter by the objects due to algae physiology. These three procedures are started at each time step for each object during the model integration.

For the microalgae object, we use a model describing the phytoplankton growth as a balance between the incoming energy and structure maintenance costs [19]. The specific growth rate and nutrient uptake rates are functions of the intracellular content of nitrogen and phosphorus:

$$\frac{d{B}_p}{dt}=\left({\mu }_p\left(1-k_p\right)-m_p\right)B_p,$$

(1)

$${\mu }_p={\mu }_{\max }\min \left(f\left(Q_{\mathrm{N}}\right),f(Q_{\mathrm{P}})\right)C_w,$$

(2)

$$f\left(Q_{\mathrm{N}}\right)=\frac{Q_{\mathrm{N} }-Q_{\mathrm{N}}^{\min }}{Q_{\mathrm{N}}^{\max }-Q_{\mathrm{N}}^{\min }},$$

(3)

$$f\left(Q_{\mathrm{P}}\right)=\frac{Q_{\mathrm{P} }-Q_{\mathrm{P}}^{\min }}{Q_{\mathrm{P}}^{\max }-Q_{\mathrm{P}}^{\min }},$$

(4)

$$\frac{d{Q}_{\mathrm{N}}}{dt}=V_{{\mathrm{NO}}_3}\frac{\left[{\mathrm{NO}}_3\right]}{K_{\mathrm{NO}3}+\left[{\mathrm{NO}}_3\right]}+V_{{\mathrm{NH}}_4}\frac{\left[{\mathrm{NH}}_4\right]}{K_{\mathrm{NH}4}+\left[{\mathrm{NH}}_4\right]}-{\mu }_p{Q}_{\mathrm{N}},$$

(5)

$$\frac{d{Q}_{\mathrm{P}}}{dt}=V_{\mathrm{P}}\frac{\left[\mathrm{P}\right]}{K_{\mathrm{P}}+\left[\mathrm{P}\right]}-{\mu }_p{Q}_{\mathrm{P}},$$

(6)

$$V_{{\mathrm{NO}}_3}=V_{{\mathrm{NO}}_3}^{\max }\left(1-f\left(Q_{\mathrm{N}}\right)\right),$$

(7)

$$V_{{\mathrm{NH}}_4}=V_{{\mathrm{NH}}_4}^{\max }\left(1-f\left(Q_{\mathrm{N}}\right)\right),$$

(8)

$$V_{\mathrm{P}}=V_{\mathrm{P}}^{\max }\left(1-f\left(Q_{\mathrm{P}}\right)\right).$$

(9)

Variables of this system of equations are listed in

Table A1. Definitions of the variables in Equations (1)–(9).

Parameter	Unit	Definition
Bp	mg DW/m3	biomass of the phytoplankton object
μp	1/day	specific growth rate
Z	103 cell/m3	phytoplankton abundance
V	μm3	cell volume
QC, QN, QP	pmol/cell	carbon, nitrogen, and phosphorus intracellular content
	pmol/cell	maximum nitrogen quota
$Q_N^{min}$	pmol/cell	minimum nitrogen quota
$Q_P^{max}$	pmol/cell	maximum phosphorus quota
$Q_P^{min}$	pmol/cell	minimum phosphorus quota
$K_{N{O}_3}$	μmol/L	nitrates half-saturation constant
$V_{N{O}_3}^{max}$	pmol/(cell day)	maximum uptake rate for nitrates
$K_{N{H}_4}$	μ mol/L	ammonium half-saturation constant
$V_{N{H}_4}^{max}$	pmol/(cell day)	maximum uptake rate for ammonium
$K_P$	μ mol/L	phosphates half-saturation constant
$V_P^{max}$	pmol/(cell day)	maximum uptake rate for phosphates
[NO3], [NH4], [P]	μmol/L	concentrations of nitrates, ammonium and phosphates in water
a	mg DW/m3	fitting parameter
b	mg DW/m3	fitting parameter
mp		mortality rate
kp		exudation rate
k1	mg/L	threshold concentration of hydrogen sulfide

. The mortality rate is estimated as $m_p=e^{-{\left(\frac{Z-a}{b}\right)}^2}$ to provide a fit to actual data on phytoplankton abundance in the simulated region. In case hydrogen sulfide appears in the water, the mortality rate increases: $m_p{}^{\prime }=m_p\left(1+\frac{\left[H_{2}S\right]}{k_1}\right)$.

In this submodule, we use ecological archive presented in [20] to evaluate their values for three phytoplanktonic species typical and predominant in the Sevastopol bay: Sceletonema, costatum (Bacillariophyta), Gymnodinium catenatum (Dinoflagellata) и Emiliania huxleyi (Coccolithophyceae). These parameters are presented in

Table A2. Species specific parameters for the phytoplankton growth model.

Parameter	Sceletonema costatum	Gimnodinium catetum	Emiliania huxleyi
μp	0.7	0.36	0.76
V	113	74.000	47
QC	1.50	7388	1.5
$K_{N{O}_3}$	2.12	7.6	0.2
$V_{N{O}_3}^{max}$	0.29	156	0.25
$K_{N{H}_4}$	0.20	33.6	0.2
$V_{N{H}_4}^{max}$	0.70	80.9	0.724
$K_P$	0.58	3.4	0.44
$V_P^{max}$	0.48	34.1	0.552
$Q_N^{max}$	0.56	150	1.312
$Q_N^{min}$	0.05	31	0.04
$Q_P^{max}$	0.035	9	0.082
$Q_P^{min}$	0.0028	1.83	0.0026
a	500	500	500
b	300	300	300
kp	0.15	0.15	0.15
k1	0.005	0.005	0.005

.

Light and dissolved inorganic compounds of carbon, nitrogen and phosphorus are resources for growth. We do not use any parameterization of photoinhibition bearing in mind that microalgae cells can use all the space within the water column and focus at the level where light intensity is comfortable for growth. Self-shading was parameterized by the coefficient in Equation (2):

$$C_w={\left(\frac{D_{ph}}{D_{th}}\right)}^{-0.75},$$

(10)

where D_ph is the average phytoplankton biomass in the nearest surroundings of the object (on the population level); D_th is the threshold value of the biomass indicated the beginning of the self-shading process.

2.2. Macroalgae Objects

Bottom seaweeds community is described as a set of modeling objects at the lower boundary of the calculated domain. Each object has a cylinder shape with a height of up to 2 m and simulates the nearest living space of a group of identical plants. Objects of this class perform the array “Phytobenthos” (Figure 1). The main method of this class is the model describing physiology of the macrophytes [18]. The dynamical equations are similar to Equations (1)–(9) with the only difference in the description of the specific growth rate. We can simulate the photosynthesis of the macroalgae in more detail due to better scientific knowledge, understanding, and parameterization of these components of the ecosystem. Therefore, we use the following formulas for modeling the photosynthesis rate P_m, and the specific growth rate of the macroalgae µ_m:

$$\frac{d{B}_m}{dt}=\left({\mu }_m\left(1-k_m\right)-m_m\right)B_m,$$

(11)

$$P_m=\left(P_{max}\tanh \left(\alpha \frac{I}{P_{max}}\right)-R_d\right)f\left(T\right);$$

(12)

$$I=I_0\exp \left(-\gamma z\right); \gamma =0.0104D_{chl}+0.0015; D_{chl}=0.012D_{ph}$$

(13)

$$f\left(T\right)=\sin \left(\frac{2\pi T}{c_1}+c_2\right);$$

(14)

$${\mu }_m=P_m{c}_3\min \left(f\left(Q_{\mathrm{N}}\right),f(Q_{\mathrm{P}})\right)C_w,$$

(15)

Variables of these equations are listed in

Table A3. Definitions of the variables in Equations (11)–(14).

Parameter	Unit	Definition
Bm	g DW/m3	biomass of an object
$P_{max}$	mg O2/(g DW h)	maximum photosynthesis rate
Pm	mg O2/(g DW h)	photosynthesis rate
µm	1/h	specific growth rate
$R_D$	mg O2/(g DW h)	dark respiration rate
α (ΔP/ΔI)	–	slope of the PI-curve at small values of radiation flux
Dph	mg DW/m3	average phytoplankton biomass above the bottom macroalgae object
Dchl	mg/m3	chlorophyll a concentration at the location
T	°C	water temperature
I	μmol photons/(m2 s)	photosynthetically active radiation
c1, c2		empirical coefficients
c3, c4		conversion factors, depending on molar ratio of C and O2 and QC
DW	g	dry weight of the algae tissue
km		exudation rate
mm		mortality rate
z	m	depth of the top of a macroalgae object

. We use an average empirical ratio between chlorophyll concentration and dry weight of phytoplankton biomass following the study [21]. The light attenuation coefficient describing the influence of phytoplankton growth on the light intensity at the bottom layer is given as a function of chlorophyll concentration according to [22]. Self-shading (coefficient C_w in Equation (14)) is parameterized by analogy with Equation (10).

To simulate the production processes in the bottom phytocenosis in the coastal zone of the Crimea, we chose seven macroalgae species typical for the rocky substrate: Cystoseira barbata, Ulva lactuca, Ceramium tenuicorne, Cladophora glomerata, Polysiphonia nigrescens, Phyllophora truncata, Enteromorpha prolifera. The model for this macroalgae community was studied and verified on annual data in [23] for a one-dimensional case.

Table A4. Species specific parameters for the macroalgae growth model.

Parameter	Cystoseira barbata	Ulva lactuca	Ceramium tenuicorne	Cladophora glomerata	Polysiphonia nigrescens	Phyllophora truncata	Enteromorpha prolifera
$P_{max}$	3.58	28.00	12.31	27.50	11.12	2.45	41.00
$R_D$	0.80	2.97	2.88	2.00	2.59	0.74	1.80
α (ΔP/ΔI)	0.12	0.45	0.13	0.13	0.11	0.04	0.15
$K_{N{H}_4}$	16.09	20.00	16.93	32.68	23.90	7.93	13.40
$V_{{\mathrm{NH}}_4}^{\max }$	37.15	211.00	143.43	327.83	138.18	9.71	188.00
$K_{N{O}_3}$	9.08	34.00	3.91	5.27	14.01	9.21	13.30
$V_{{\mathrm{NO}}_3}^{\max }$	8.65	116.00	18.67	115.72	49.89	1.69	169.00
$K_P$	0.29	4.41	1.00	0.33	4.47	0.37	1.51
$V_{\mathrm{P}}^{\max }$	0.06	9.31	0.78	3.06	7.39	0.12	4.59

shows the photosynthetic and kinetic parameters used in our study. The extremum nutrients contents in tissue were set equal for all macroalgae: $Q_{\mathrm{P}}^{\min }$ = 28 μmol P/g DW; $Q_{\mathrm{N}}^{\min }$= 850 μmol N/g DW; $Q_{\mathrm{P}}^{\max }$ = 80 μmol P/g DW; $Q_{\mathrm{N}}^{\max }$ = 2570 μmol N/g DW.

Macroalgae objects interact with the environment in the same manner as the phytoplankton objects. Concentrations of nitrogen and phosphorus compounds (both inorganic and organic) and oxygen changed under the influence of uptake and release of matter by the macrophytes. Objects containing algae of different species can exist in the same place or intersect each other simulating coexistence of the epiphytes and host plants and their competition for nutrients.

2.3. Hydrochemistry Module

Dynamics of water chemistry were described by the model for nitrogen and sulfur cycles in the Black Sea [24]. Despite the fact that this model was initially developed for the redox zone of the Black Sea, it can be applied to the simulation of anoxic conditions that are possible in the eutrophic coastal waters during the summer stagnation. Even in a non-eutrophied system, restricted vertical mixing can lead to anoxia. High concentrations of nutrients can dramatically worsen the situation. The model is able to mathematically describe the transition of normal oxic conditions to anoxic ones. It parameterizes the cycles of nitrogen and sulfur including the processes of ammonification, nitrification, nitrate reduction, denitrification, three stages of hydrogen sulfide oxidation, two stages of sulphate reduction, and thiodenitrification (Figure 2). The concentration of phosphorus was added to the model using the Redfield ratio to simulate macro- and micro algae physiological processes more accurately (uptake, internal stocks, possible limitation of growth by nitrogen and phosphorus).

Hydrochemical fields were represented as arrays of discrete values in the cells of a regular grid. Coupling with the hydrodynamic module was performed using a special procedure. The concentration of the chemical element in each grid cell was considered as a sum of changes caused by dynamical factors, chemical processes, and biologically induced changes:

$$\frac{dC}{dt}=\Delta C_{dyn}+\Delta C_{chem}+\Delta C_{bio}.$$

(16)

The vertical distribution of water temperature is determined by heat flux and wind stress at the surface in accordance with the simple integral model [25]. Two layers with a homogeneous vertical distribution of temperature and chemical compounds can exchange water volumes under certain influences at the air–sea boundary. General equations for the chemical components of the model include advective and diffusive terms, exchange between layers, chemically and biologically induced terms:

$$\frac{\partial X_i^u}{\partial t}+U^u\frac{\partial X_i^u}{\partial x}+V^u\frac{\partial X_i^u}{\partial y}=K_x\left(\frac{{\partial }^2{X}_i^u}{\partial x^2}+\frac{{\partial }^2{X}_i^u}{\partial y^2}\right)+G_i^u+A_i^u+Q_h\left(X_i^b-X_i^u\right)+Q_0,$$

(17)

$$\frac{\partial X_i^b}{\partial t}+U^b\frac{\partial X_i^b}{\partial x}+V^b\frac{\partial X_i^b}{\partial y}=K_x\left(\frac{{\partial }^2{X}_i^b}{\partial x^2}+\frac{{\partial }^2{X}_i^b}{\partial y^2}\right)+G_i^b+A_i^b+Q_h\left(X_i^u-X_i^b\right),$$

(18)

where U and V are the averaged current velocities in the corresponding layers; indices u and b indicate the upper and bottom layers; X_i are chemical components of the model; G_i are the terms induced by chemical processes; and A_i are the inputs and uptake by biological components. Q₀ and Q_h are fluxes at the surface and intermediate boundaries.

A detailed description of the chemical model and numeric values of the model parameters can be found in [25]. Here, we show the equations for dissolved oxygen in the upper and bottom layers as an example (Equations (18) and (19)).

$$\begin{array}{c}\frac{dD{O}^u}{dt}=\Delta D{O}_{hor}^u+\beta \left(D{O}_{sat}-D{O}^u\right)f\left(T,W^2\right)+Q_h\left(D{O}^b-D{O}^u\right)+P_p{B}_p-\\ p_{1}Am-p_{2}N{f}_1-p_{3}N{f}_2 -p_{4}H{s}_1-p_{5}H{s}_2-p_{6}H{s}_3\end{array}$$

(19)

where $\Delta D{O}_{hor}^u$ is the horizontal advection and diffusion input; DO_sat is the saturating oxygen concentration; β is the coefficient of air–sea gas exchange; T is the temperature; W is the wind speed above the sea; Q_h is the flow between the two layers (when the upper layer is deepening Q_h = $\frac{dh}{h dt }$, otherwise Q_h = 0); h is the upper layer depth; P_pB_p is the oxygen production during photosynthesis of phytoplankton (P_p = c₄µ_p); p₁Am is the oxygen expenses for oxidation of organic matter; p_iNf_i are oxygen expenses for nitrification; p_iHs_i are oxygen expenses for hydrogen sulfide oxidation; p_i are stoichiometric coefficients. For DO in the bottom layer, the equation is:

$$\begin{array}{c}\frac{dD{O}^b}{dt}=\Delta D{O}_{hor}^b+Q_h\left(D{O}^b-D{O}^u\right)+P_p{B}_p+P_m{B}_m- p_{1}Am-p_{2}N{f}_1-\\ p_{3}N{f}_2 -p_{4}H{s}_1-p_{5}H{s}_2-p_{6}H{s}_3\end{array}$$

(20)

where Q_h = $\frac{dh}{dt\left(H-h\right)}$ if the upper boundary of the bottom layer rises (Δh < 0), and Q_h = 0 if not; P_mB_m is the oxygen production during the macroalgae photosynthesis.

Particular organic matter sinks at a set fixed rate. Therefore, some part of it passes from the upper layer to the lower one, thereby increasing the likelihood of anoxia.

2.4. Exchange Procedure

Algae objects can interact with the environment provided by the chemistry module in two ways. They can estimate the average concentration of nutrients within the occupied space and change their uptake rate in accordance with this value and their own internal reserves (Equations (5)–(8)). Concentrations of inorganic nutrients change due to the uptake by autotrophs. The second way is the release of matter by the algae objects and corresponding alteration of dissolved and particulate organic matter and oxygen concentrations in grid cells.

Available inorganic resources $R_i^{\mathrm{N}}$ and $R_i^{\mathrm{P}}$ are calculated at each time step as the average concentration of nitrogen and phosphorus at the location of the object:

$$R_i^{\mathrm{N}}={\displaystyle \sum }_{j=1}^M{s}_{ij}\left(D_j^{\mathrm{N}}-D_0^{\mathrm{N}}\right); R_i^{\mathrm{P}}={\displaystyle \sum }_{j=1}^M{s}_{ij}\left(D_j^{\mathrm{P}}-D_0^{\mathrm{P}}\right); s_{ij}=S_{ij}/S_i$$

(21)

where D_j^N and D_j^P are the concentrations of nitrates, ammonium and phosphates at a j grid cell, which is partially occupied by the i algae object; M is the number of grid cells partially occupied by an object; $D_0^{\mathrm{N}}$ and $D_0^{\mathrm{P}}$ are the threshold concentrations; S_ij is the area of a grid cell overlap by an object in the horizontal projection; S_i is the basal area of an object. $R_i^{\mathrm{N}}$ and $R_i^{\mathrm{P}}$ were used as inputs to the growth procedure of an object, and, at the output, according to Equations (1)–(9), we obtain the values of the actual uptake of nitrates and ammonium, released oxygen, dissolved and particular organic matter.

Output values should be transferred into the corresponding fields, and the exchange procedure does it in the following manner. The actual uptake of nonorganic nitrogen and phosphorus (µmol/m³) is calculated as:

$$N_{up}=\left(B^{t+1}{Q}_{\mathrm{N}}^{t+1}-B^t{Q}_{\mathrm{N}}^t\right)+N_{ex}+N_{mort},$$

(22)

$$P_{up}=\left(B^{t+1}{Q}_{\mathrm{P}}^{t+1}-B^t{Q}_{\mathrm{P}}^t\right)+\frac{N_{ex}+N_{mort}}{16},$$

(23)

where indices t + 1 and t indicate discrete time steps; N_ex = dt·µkB is the exudation of dissolved organic matter during a time step, and N_mort = dt·mB is the particular organic matter released due to decay. N_ex enters the pool of DOM, and N_mort—the POM pool. The rate of exudation k_p in Equation (1) and k_m in Equation (11) are defined by 15% for microalgae and 25% for macrophytes in the light and half this value in the dark [26,27]. The rate of the oxygen release due to photosynthesis is defined as P_m (Equation (12)) for macroalgae and in proportion to the specific growth rate for microalgae (P_p). These additional negative or positive inputs are distributed between the grid cells occupied by an object in the same proportion as during the estimation of available resources in Equation (21).

2.5. Sensitivity Analysis

The sensitivity analysis of the submodules simulating phytoplankton and macroalgae dynamics was considered in the previously published papers [19,23]. In case of a sharp drop in the nutrients supply, the macroalgae community adapted to new conditions in approximately 20 days. Only two slowly growing species needed more time for adaptation. The new stable state was characterized by lower growth rates and higher uptake rates. The internal content of the nutrients also decreased, but relative change of this parameter varied from 2 to 4 when the external concentration dropped 10 fold.

It was found that the reaction of the system to a decrease in concentrations depends on the molar N/P ratio in water and their initial concentration [23]. With an initially high concentration of nutrients in water, the relative decrease in the specific growth rate was significantly less than with low initial concentrations. With a lower N/P ratio, the decrease in specific growth rates of most macroalgae was less, and the relative increase in the uptake rates was greater. We explained this phenomenon by the fact that phosphorus acted as a limiting element for algae growth with an increase in the N/P ratio.

The phytoplankton adapted to the step-wise changes of the nutrients supply much faster. A sharp increase in the nutrient uptake rates and flexible stoichiometry allow some taxonomic groups to keep almost the same growth rate during a 10-fold drop of the external concentrations. According to our estimates [19], insolation has a major effect on the phytoplankton dynamics. The specific growth rate can slightly increase and then decrease if insolation continues growing. The combined effect of irradiance and temperature rise leads to extremely fast growth of microalgae abundance and biomass. Due to the rapid adaptation of algae cells to changes in the nutrients supply, seasonal variability of phytoplankton abundance is mainly defined by seasonal changes of insolation and water temperature.

3. Results and Discussion

We have applied the model to study the case when the discharge of contaminated water sharply increases the nutrient concentration in the area with difficult water exchange. Usual nitrate concentration in the Sevastopol Bay varies from 0 to 4 μmol/L, and maximum values are observed in the surface layer of the South Bay (Figure 3). Extremum nitrate concentrations are caused by discharges of polluted water and can achieve 40–50 μmol/L at the upper layer of the South Bay. We have studied the process of self-purification of this basin by carrying out the simulations with the model including various modules. The initial conditions and the control variables were the same for each experiment. For the first simulation, we use only the main part of the ecological model—the hydrochemical module. In this experiment, the decomposition of organics and the oxidation of reduced forms could only occur due to the DO available in the system. In the second experiment, we added a set of phytoplankton objects controlled by the module of the phytoplankton growth. Finally, the third simulation was fulfilled with both phytoplankton and macroalgae sets of objects that provided additional inputs of DO. Then, we comprised and analyzed the obtained results to study the contribution of various components of the coastal ecosystem to the process of self-purification. To consider the case of limited water exchange in the top part of the South Bay, we excluded the advection terms from the left sides of Equations (16) and (17).

The hydrological conditions simulated summer stagnation with low wind and quasi-constant water temperature for one month. The time step was 0.125 day (3 h), and the space step was 50 m. The initial concentrations of the hydrochemical components are listed in

Table A5. Initial values of the chemical and biological components of the model.

Component	Initial Value
	Upper Layer	Bottom Layer
[O2]	7 mL/L	6 mL/L
[DOM]	15 μmol N/L	20 μmol N/L
[POM]	15 μmol N/L	20 μmol N/L
[NO3]	30 μmol N/L	5 μmol N/L
[NH4]	4 μmol N/L	1 μmol N/L
[NO2]	0.5 μmol N/L	0.2 μmol N/L
[PO4]	0.05 μmol P/L	0.01 μmol P/L
T	23 °C	15 °C
h	5 m
Bp	135 mg/m3
Bm	160 g/m2

.

3.1. Dissolved Oxygen

We estimated the ecological state of the basin by DO concentration in water, the rate of water purification from excessive concentrations of nutrients, and the likelihood of anoxia. Figure 4 shows the DO concentrations together with the hydrological conditions in the upper and bottom layers in all three simulations. Such conditions, with low wind speed above a semi-enclosed bay, are typical for summer.

In the first simulation without phytoplankton and macroalgae, DO in the lower layer was spent rather quickly and its content was reduced so much that traces of hydrogen sulfide appeared after 13 model days and its concentration in the lower layer increased over time (Figure A2). In the second simulation with phytoplankton objects in the ecosystem, the trigger for hypoxia occurred almost at the same time. Photosynthesis of phytoplankton was accompanied not only by the oxygen production but also by DOM and POM release. Therefore, phytoplankton growth was unable to prevent contamination of the bottom layer with hydrogen sulfide.

Inclusion of the macroalgae objects fundamentally changed the situation in the bottom layer (Figure 4b). DO content was stable during the simulation, showing only slight diurnal fluctuations, and did not drop below 4 mL/L. The stable states of the phytoplankton and macroalgae biomass under given conditions were 220 mg DW/m³ and 440 g DW/m², the initial biomass was set as shown in Table A5. We chose low initial values for simulations to show that even with such biomass of bottom macroalgae, the case of sharp deterioration of the water quality could be recovered by self-purification. At higher seaweeds biomass, the extraction of nitrogen and the saturation of the bottom layer with oxygen occurred more intensely. For example, at a macroalgae biomass of 500 g DW/m² the nitrate concentration in the lower layer (under the same initial conditions) decreased to 1 µmol during 5 model days, and DO concentration in the bottom layer stayed near the saturation level.

3.2. Nitrates

Consider nitrogen dynamics and compare graphs of nitrate concentration in these three simulations (Figure 5). In the first experiment, nitrate concentration grew due to oxidation of ammonia and organic matter whose contents in water were rather high (Figure 6 and Figure 7, Figure A3). Phytoplankton cells used nitrates and ammonia for growth; therefore, in the second experiment, we expected that the inorganic nitrogen would decrease. However, nitrate concentration in the lower layer was growing in the second experiment. This can be explained by input from the organic matter compound due to the ammonification process. This is confirmed by comparing the graphs in Figure 5, Figure 6 and Figure 7 and Figure A2 and Figure A3. DOM in the bottom layer was growing in this simulation mainly through the phytoplankton exudation (Figure 6). Additional input of the organic matter was coming from the sinking POM and its partial decomposition and dissolution (Figure 7). In the third experiment, inorganic nitrogen in the bottom layer slightly decreased in comparison with the initial state. In the upper layer, average nitrate concentration decreased due to phytoplankton uptake and water exchange between layers.

3.3. Dissolved Organic Matter

Graphs of the DOM dynamics differ from the nitrate ones. In the first simulation, organic matter was oxidizing without any additional sources. Therefore, its content in both layers decreased. In the second and third simulations, phytoplankton and seaweeds released DOM by the exudation process. The additional input of partially dissolving POM increased DOM concentration. This led to the additional oxygen demand. Maximum DOM concentration was observed in the third experiment at the bottom layer where the conditions were insufficient for its fast oxidation.

3.4. Sources and Sinks of DO

It is of interest to analyze the integral sources and sinks of DO in these three experiments (Figure 8). Averaging was performed over the entire period of the experiment (1 month). In the first simulation, sources of DO are air–sea flow and internal DO stocks (Figure 8a). The sinks of DO are costs on the organic matter oxidation and nitrification. DO expenses on hydrogen sulfide oxidation are very small, so we do not show them on the graphs. In the second experiment, an additional source of DO is provided by phytoplankton growth (Figure 8b), and in the third experiment by phytoplankton and macroalgae growth (Figure 8c).

In the third experiment, the system state remains oxidizing and no hydrogen sulfide appears in the lower layer of the model basin (Figure A2). However, oxygen demand for organic matter oxidation increases by 43% from 276 to 395 mL/(m³ day) in comparison with the second simulation (Figure 8c). Additional organic matter is a product of seaweeds functioning (DOM release in the form of exudates and POM input due to decay).

Consider the rate of oxygen release by various types of autotrophs in a model ecosystem. Figure 9 shows the specific and total rates of DO release per day. According to the diagram, the specific oxygen production of microalgae is several times higher than that of macroalgae, which is connected with a higher rate of microalgae photosynthesis and, accordingly, its growth. In this experiment, the maximum specific growth rate of phytoplankton was set as 0.76 day⁻¹. The maximum rate of macroalgae did not exceed 0.11 day⁻¹.

Although phytoplankton has higher rates of DO production, it cannot be considered as the factor of self-purification of the ecosystem. Only if the primary production of phytoplankton is needed for the following trophic levels, the balance between the production and demand of DO can be positive. During microalgae blooms, DO demand for the oxidation of a huge amount of produced organic matter is bigger than DO production, especially since a big part of it flows to the atmosphere through the air–sea border.

3.5. Particular Organic Matter

We can draw the analogous diagrams (as in Figure 9a,b) for nitrogen and phosphorus uptake, release of DOM and POM by phytoplankton and macroalgae species. These low-level characteristics can help us to deepen our understanding of the processes on the population level. Consider the ratio of the specific rates of dissolved and particular organic matter release (Figure 10). In our model, POM is considered as organic matter consisting of two parts: dead cells or tissues, and those consumed by the animals at the next trophic level. Mortality rates are the most uncertain parameters in the model. We estimated these coefficients based on the requirements of biomass stability and correspondence to the values observed in the region [23]. At the same time, DOM release rates are proportional to the growth rates and can be roughly estimated in the laboratory experiments.

According to the diagram in Figure 10, POM release by phytoplankton is much higher in comparison with the macrophytes. These estimates may be overstated because this variant of the model does not include the next trophic levels—marine animals. Consumption of the part of the primary production prolongs the time of the organic matter decomposition, but, in a semi-enclosed basin, it cannot sufficiently influence the result. Therefore, primary production of phytoplankton is accompanied by much higher oxygen demand for the decomposition of organic matter in comparison with the macrophytes.

Compare now the diagrams in Figure 9b and Figure 10. Green macrophytes (Ulva, Cladophora and Enteromorpha) make an increased contribution to DO supply. However, their POM/DOM ratios are also higher than for other macroalgae species. Cystoseira, the basis of the bottom macrophytes community, shows the best relationship between DO production and DO demand on the decomposition of produced organic matter. This is why, by our estimates, this seaweed plays the major role in self-purification of the basin in the situation of excessive nutrient input in the coastal waters.

4. Conclusions

The experiments have shown that the proposed model is an effective tool for studying and parameterizing the influence of phytoplankton, marine seaweeds, and (in future) aquaculture facilities on the concentration of organic matter and oxygen. An object-based approach allows one to analyze biochemical mechanisms of the autotrophs’ physiology (nutrient uptake, growth, oxygen production, exudation, etc.) on both individual and population levels. Incorporation of the objects containing model plants with laboratory-estimated properties into a 3D ecological model simulating the marine environment is of great interest. Competition for nutrients and light between species and self-shading modify the rates of the main processes and can lead to unpredictable changes of the state variables of the ecosystem. The proposed modeling technique promotes a better understanding of internal ecosystem dynamics by identifying the causes of possible changes. Individual-based bottom-up modeling gives the possibility of tracking the variability of the internal state variables, such as the nutrients stored in phytoplankton cells and macroalgae tissues and the internal stoichiometric ratio of main elements.

In the test simulations discussed above, we study the reaction of different plant components of the ecosystem on sharp increase of inorganic nitrogen in water. We compare micro and macroalgae species with very different rates of the main physiological processes. Analyzing the results of experiments, we have shown why the phytoplankton community with much higher rates of oxygen production could not prevent hydrogen sulfide pollution in the conditions of low mixing and high nutrient concentration in water.

The proposed approach could explain and demonstrate the mechanism of the high internal plasticity of marine autotrophs, allowing them to continue growing even in severe trophic conditions. This approach is very convenient for expanding the model structure. Using its principles and rules, we can easily add new species and components to the model. This is also true for aquaculture facilities. Further development of the ecosystem model and embedding the next trophic level components can provide a new tool for studying the relationships between natural and artificial components of the marine coastal ecosystems including aquaculture.