Development of a Mathematical Model for Simulation of Macroalgae Farming in the Coastal Areas

A mathematical model consisting of a system of three coupled partial differential equations (PDEs) was proposed to estimate the concentrations of nitrogen, phosphorous and macroalgae biomass in coastal open waters. However, some simplifying assumptions were used in the model to cope with the complexity of real conditions. For the macroalgae biomass, the system works as a batch mode, while input and output were accounted for nitrogen and phosphorous. The MATLAB pdepe feature, applying the finite element method was used in model solving and the simulation of model equations. The program was split into four functions that included the solver and post-processing of the results, a function containing the PDEs, a function setting the initial conditions, and one setting the boundary conditions. For model validation, the experimental measurement of nitrogen, phosphorous and macroalgae biomass concentrations of Bandar Abbas coastal open waters were analyzed by standard methods at three depths of 1, 5 and 10 m. The predictive values of the developed model demonstrated its applicability for the management of coastal macroalgae cultivation systems by assessing the impact of nitrogen and phosphorous strategies on the farming system.


Introduction
acroalgae growth typically occurs in natural and pond systems.The productivity of marine ecosystems is characterized by various biological, physical-chemical and meteorological variables such as nutrient availability, water temperature, salinity, grazers, phytoplankton, etc. [1,2].The sustainability of macroalgae cultivation in coastal open waters is essential from environmental and economic perspectives, and for human communities [3].Independent and integrated cultivation of macroalgae in both open waters and fish ponds have been investigated [4,5].

M
The results of these evaluations demonstrate the crucial influence of the availability of various forms of nitrogen and phosphorous, and their effects on biomass productivity in different seasons and at various water depth levels [6].For example, nitrate concentrations in the range of 0.1 -1.0 µM were evaluated in the cultures of some species and found to have a significant effect on biomass productivity [5].In another study nitrogen sources (sum of ammonium and nitrate) in the range of 5 to 30 mg l -1 were studied [7].These findings suggest the requirements for managing the complex behavior of water nutrients, which are critical in macroalgae cultivation technologies [7].
Modeling of macroalgae growth is difficult as the behavior of biotic and abiotic macroalgae cultivation systems are complex, and have defied precise description using mathematical models [8,9].Furthermore, it is currently impossible to reliably predict how various modifications of system parameters might affect cultivation performance, because many parameters are not included in the model equations.However, mathematical models can be effective tools for predicting macroalgae farming productivities, by analyzing and simulating nutrient removal and the macroalgae biomass growth rate [1,2,10,11].A MARS-Ulves model was proposed to predict nitrate threshold values to control macroalgae blooms, due to anthropogenic nitrogen loading [7].The Ulva growth and nutrient assimilation were predicted by an ecological model with 8 biotic variables and 2 abiotic variables, of temperature and lighting.However, ecological models such as MARS-Ulves mostly consider hydrodynamic features which should be derived experimentally [12,13].In addition, they are complex and require various input data [14][15][16].Other classes of models have been used in natural and artificial wastewater treatment technologies such as treatment pond systems (e.g. the Reed model) which might be used for macroalgae cultivation [17][18][19].In this regard, several removal pathways have been investigated for the reduction of nitrogen and phosphorous including nitrification-denitrification, algal uptake, and sedimentation, etc. [18].These variations complicate the development of a robust mathematical expression to account for all these processes simultaneously.Another applicable class of equations for macroalgae farming were developed for fall in modeling of photobioreactors.These models are usually focused on micro and macroalgae biomass productivity and optimum nutrient availability for maximizing productivities in artificial open and closed algae cultivation systems [20].The typical models in these category are those which postulate time variations of variables such as lighting, gas concentrations, flow characteristics, etc. [18].
All currently applied models assume a homogeneous farming system in depth layers.This approximation is less valid in depths of greater than 3 meters.They are also less applicable when light and nutrients have strong gradients with depth [11,21].For such cases, we need a model which considers both spatial and temporal changes during cultivation.In other words, the biomass productivity can be better predicted by developing a depth dependent model.This research provides a greatly simplified but effective model which takes time and depth variations during cultivation into account, linking the biomass growth model to nitrogen and phosphorous concentrations.In other words, we will analyze the biomass productivity and nitrogen and phosphorous assimilation using PDEs and then examine the graphs generated by the model.These variables will affect both farm productivity and the safety of coastal waters.These simulations include the production of macroalgae maps at 3 depths of 1, 5 and 10 m against time.

Model formulation
One graphical and one conceptual schematic representations of the system Ω (z, t) studied is shown in Figure 1, (a) and (b).Domain Ω is defined on a rectangular (x, z) plane with 0 < x < 1 and 0 < z < H.The height of the domain is H in z direction.The macroalgae biomass is represented by green shapes, and the nitrogen and phosphorous by white circles.The light depletion with water depth is shown by a color gradient from yellow to gray.Boxes indicate the components represented as state variables in the mathematical model: NO, combined nitrite and nitrate, P, phosphorous, and biomass B, as an organic particles in the form of macroalgae.NO 3inlet = 1 gl (homogeneous), and initial concentrations of N and P were assumed zero in selected control volumes [3].
The system contains a mixture of seawater and suspended macroalgae B, where the biomass is grown by uptaking N and P variable compounds from water [19].Macroalgae biomass B, nitrogen N and phosphorous P concentrations were considered as state variables.The domain is exposed to air from above and no diffusible layer from bottom.For the purpose of describing the biomass growth mathematically, we consider a 1-D water column domain 0 < z < H with the z = 0 in the bottom and the bulk water interface at z = H , and set B (t, z) to be the concentration of biomass.Height H is assumed constant.Biomass B satisfies a no-flux condition B z | z=0 = 0 at the water bottom and a Dirichlet condition B z | z=h = B 0 at the water interface.A diffusive boundary layer above the interface z = H could be included to allow for mass transfer resistance but does not qualitatively affect results.We suppose the sides of the system to be insulated.We also suppose initial conditions B (z, 0) = 0.03, 0 < z < H, i.e., application beginning at the t = 0 of macroalgae biomass of 0.03 gl -1 .

Figure 1 (a). A graphical view of the entire domain Ω (t, z). (b). A conceptual model scheme of nitrogen and phosphorous input, transformation and removal processes in domain Ω (t, z).
In the processes incorporated as PDEs, B (z, t) denotes the biomass change at each time, N (z, t) denotes the nitrogen change at each point and z is the water depth, and P (x, t) denotes the phosphorous change at each point and z is the water depth.The system consists of a column of water, with variables of biomass, B, nitrogen, N and phosphorous, P concentrations (see Figure 1).Assimilation of nitrogen and phosphorous which is used for macroalgae feeding is not uniform in different water depths, as the profile of macroalgae concentration varies with water depths.In this study, standard convection-diffusion-reaction (CDR) systems of partial differential equations were used for macroalgae biomass production and nutrients in a one-dimensional static domain [22,23].

Biomass
The aim of a biomass model equation is to describe the spatial profile of biomass B in (x, y) plane as a function of time in order to predict the effects of nitrogen and phosphorous penetration barriers.Let  ≔ ( 0 ,   ) × ( 0 ,   ) × ( 0 ,   ) ⊂  3 be the system volume domain, in which biomass is growing, decaying or both.For modeling, the fact that the biomass growth µ and decay β rate is variable at any given depth location H inside of the domain Ω, we impose that the biomass accumulate/degrade equation in 3 levels with µ 1 , µ 2 and, µ 3 , in order to take into account different biomass growth rates.That is, the biomass change in the domain is not divergent.In this case, the biomass concentration B satisfies a simplified convection-diffusion-reaction (CDR) equation of the form which is a function of three variables defined on a region in Ω, where B is a vector of biomass (gl -1 ), f is a vector showing the kinetics of biomass growth and decay, and D B stands for the diffusion coefficient of the biomass (m 2 s −1 ).We set (, ) =  −   .
(2) Here B is the biomass concentration.Note that here, the specific growth rate µ depends on various other nitrogen and phosphorous factors, but broadly it depends on several others such as light intensity, temperature, etc., and the extra complications are neglected for simplicity.
Using (2), equation ( 1) can be rewritten as where the μ and β coefficients represent biomass growth and loss, respectively.In order to uniquely determine the biomass profile from Eq. ( 3), we have to add external initial and boundary conditions.In this case, we assume that no biomass goes in or out by the bottom of Ω: B z | z=h = B.Moreover, for numerical convenience, we assume that B is periodic in the x and y-directions.On the other hand, assuming that on the top of the domain there is a constant biomass concentration, and since this is always defined up to an additive constant, without loss of generality, we will impose the following boundary condition: .

Nitrogen and phosphorous
For the purpose of simplicity, we assume that µ depends only on the concentration of limiting substrates, e.g., N and P. In particular, we assume that µ = µ (t,U(N,P)), which is the so called substrate uptake rate function, which indicates the reaction rate of substrate usage.We used Monod kinetics for N and P uptake by macroalgae.Here we define µ by the equation where µ max is the maximum biomass growth rate.The presence of macroalgae in the domain influence the N and P by assimilation.The equation which describe the diffusion and assimilation of N is   =    2   2 + ()  +  1 () −   . ( Equation 5 has two boundary conditions: one at top water surface and the other at the bottom water.The concentration of nitrogen and phosphorous were assumed constant at the top water surface or N| z = L Z = S m (the substrate concentration is at its maximum level at the top of Ω) and the flux of nitrogen and phosphorous were assumed zero at the bottom surface or ∂N/∂t| z = 0 (no nitrogen goes in or out by the bottom of Ω).Similar equation and boundary conditions to N were also used for the phosphorous profile as follow: The top boundary condition is P| z = L Z = S m (the substrate concentration is at its maximum level at the top of Ω) and ∂P/∂t| z = 0 (no phosphorous goes in or out by the bottom of Ω).In the above equations ( 5) and ( 6), D N and D P stand for the diffusion coefficients of the nitrogen and phosphorous substrates, respectively (m 2 s −1 ).Over longer durations, (5 and 6) equilibrate to C (N) = k N zz and C (P) = k P zz .Note that U and µ do not explicitly depend on time or space in the model presented here.In low water layers, where B is small and C (B) = C 1 B, nitrogen concentration is depleted exponentially with length scale l = k/C 1 .Thus for those nutrients like nitrogen and phosphorous, a permanent reactive penetration layer results (at least, permanent on the time scales considered here).
At the top of the water column, if the concentration of biomass is saturating, then (3) simplifies to, approximately, where the constant D B is biomass diffusivity, and the constant C 0 is the saturation level of C (B), so that nitrogen and phosphorous show a decreasing profile.
Gathering the above PDE equations, we have the following model, with incorporated initial and boundary conditions.
The model has been translated to be accessible by the input solver format of the MATLAB pdepe feature, which is used in the solving and simulation of model equations.In the pdepe feature, the ordinary differential equations (ODEs) resulting from discretization in space are integrated to obtain approximate solutions at specified times to solve initial-boundary value problems for parabolic-elliptic PDEs in 1-D.The model parameters of maximum growth rate, coefficient of biomass decay, nitrogen and phosphorous generation results of biomass decay, diffusion coefficients, decomposition of biomass based nitrogen compounds d 1 , decomposition of biomass based phosphorous compounds d 2 , and yield were obtained from literature for Gracilariopsis persica in the Persian Gulf or by experimental analysis as shown in Table 1.
To use the MATLAB pdepe feature for solving the model PDEs, we split the program into four part functions that included the solver and post-processing of the results, a function containing the PDEs, a function setting the initial conditions, and one setting the boundary conditions.The equations were translated into the MATLAB format.The outputs were set to obtain profiles of state variables.

Materials and methods
The nitrogen, phosphorous and macroalga Gracilariopsis persica concentrations in Bandar Abbas coastal open waters were determined in September 2016 and January 2017 according to the standard methods for 3 depths of 1, 5 and 10 m.Experimental data used for model validation [5,27] were obtained from analysis of three water depths, i.e., n = 1, 5 and 10 m.

Results and discussion
We have performed numerical experiments with the corresponding initial values for N and, P and for different initial biomass concentrations.The biomass specific growth rate µ and decay β can affect the final biomass B during cultivation time.It also predicts higher biomass productivity in different depth layers.The greater the specific growth rate µ, the higher resulting biomass productivity in a shorter time.Some results and analysis of macroalgae growth and nutrient assimilation are presented below.

Time-and depth-dependent macroalgae growth
The solution to the model for single state variable of biomass is shown in Figure 2. The simulated results revealed a time-and depth-dependent growth pattern, decreasing from top down, corresponding to the faster growth empirically observed in cultures at depths of 1 and 5 m than at 10 m depth cultures (Figure 2).The simulated results suggest utilizing a cultivation depth of ≥ 5 meters.The model estimate should be refined for differing macroalgae species and geographical locations.

Time-and depth-dependent nitrogen and phosphorous
The solution to the model for N and P variations is shown in Figure 3.The simulated results revealed a time-and depth-dependent growth pattern with differing concentrations of nitrogen and phosphorous (Figure 3).The purpose of this research was to develop a mathematical model to assist in optimizing the cultivation performance and design of macroalgae cultivation systems in open waters.Nitrogen N and phosphorous P are important elements for the growth of photosynthetic organisms in marine environments including seagrasses, macroalgae, and microalgae.Their availability can improve macroalgae farms' productivity, but at the same time may cause eutrophic phenomena which contribute to harmful conditions such as coastal red tide blooms [11,28].Therefore, control of inorganic water nitrogen and phosphorous concentrations are important in macroalgal farming activities [10,[29][30][31].The proposed model is different from typical existing models, since it simulates macroalgae productivity as a function of spatial and temporal variations by means of spreading biomass reaction, instead of convective biomass under the influence of substrate diffusion.Indeed, the conditions describing the productivity in various water levels are a key issue in the presented model that have been extensively used in macroalgae productivity problems.The novel a c b approach to biofilm modeling followed in the paper is, as far as we know, applied for the first time to the cultivation of Gracilariopsis in the Persian Gulf area, and hence introduces local communities to a method that is likely to find application [3].Using this model, numerical simulations were performed that predict the behavior of state variables in domain in a range of different conditions.As the nutrient availability increases, there is a gradual shift towards more productivity in water depth layers.This is consistent with the fact that the greater the surface light and nutrient availability (and so the higher the growth rate is), the better suited it is for commercial macroalgae production systems [7].
The model simulation results suggest that compensation for lower than required nitrogen and phosphorous sources in coastal waters might be supplied by local wastewater, to increase macroalgal farms' productivity [32,33].However, the seasonal and depth variation of nitrogen and phosphorous sources in natural water are also dependent on other variables such as light and temperature, and should be considered for more accurate simulation of macroalgae cultures.Integration of macroalgae cultivation with fish farming is another promising option for higher nutrient sources and economic applicability [5,34].The model could consider two or more substrates interacting with each other and biomass productivity.This amounts to introducing new equations to the model, taking into account transport and interaction of the various substrates.The model could also consider systems with multiple water species.
The most relevant characteristic of the present novel and more rigorous model approach, based on existing mathematical analysis, is that it can simulate and follow the behavior of a range of previously described models that simulate practical macroalgae cultivation behavior.Detailed comparisons with experimental data including the effect of seasons and light irradiation require future research.

Conclusion
This study provides an estimate of macroalgae productivity as a function of inorganic nitrogen and phosphorous concentrations by a model of three PDE equations.The important features of the model for improving understanding of the macroalgae cultivation system are its (1) parameterization of the link between nitrogen, phosphorous and biomass in terms of N and P-to-biomass ratio, and its (2) linking macroalgae growth to nitrogen-phosphorous assimilation.Simulations were in reasonable agreement with the experimental data obtained from analysis of different depths of water, and illustrated the strong dependence of biomass productivity on N and P concentrations.The values predicted by the developed model are useful in providing insight on experimental data.The model may be used for further investigation of coastal macroalgae farming for the development of local communities.We are now working on refining the model by including fluid hydrodynamic and growth related parameters, in order to simulate for higher macroalgae productivity, as well as for sustainability of the coastal waters.

Figure 2 .
Figure 2.Estimated concentrations of macroalgae biomass influenced by water depth.

Figure 3 .
Figure 3.Estimated concentrations of (a) nitrogen and (b) phosphorous with water depth.

Figure 4
Figure 4 compares the measured values of biomass, nitrogen and phosphorous concentrations for three depths of 1, 5 and 10 m with the corresponding simulated values.Correlation coefficients (R 2 ) with 95% confidence intervals between the model predictions and the measured biomass, nitrogen and phosphorous for different depths (Figure 4) were 0.97, 0.95, and 0.96, respectively.These correlation coefficients indicate that the mathematical model successfully fitted the experimental data.

Table 1 .
Model parameters used for simulation of Gracilariopsis persica.