SQU Journal for Science, 2018, 23(1), 32-42 DOI: http://dx.doi.org/10.24200/squjs.vol23iss1pp32-42 Sultan Qaboos University 32 Development of a Mathematical Model for Simulation of Macroalgae Farming in the Coastal Areas Abdolmajid Lababpour Department of Mechanical Engineering, Faculty of Engineering, Shohadaye Hoveizeh University of Technology, Dasht-e Azadeghan 64418-78986, Iran. Email: lababpour@shhut.ac.ir ABSTRACT: 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. Keywords: Modeling algal growth; Macroalgae cultivation; Algae simulation; Coastal open water; Numerical solution of PDEs. الساحلیة المناطق في زراعة الطحالب الكبیرة لمحاكاة رياضي نموذج تنمیه عبدالمجید لباب پور الحيوية والكتلة والفوسفور النيتروجين تركيزات لتقدير ثالثة (PDEs) جزئية تفاضلية معادالت يتألف من جملة رياضي نموذج حتم اقترا :صالملخ الكتلة بالنسبة إلى .حقيقيةال الشروط تعقيدات مواجهةفي النموذج تبسيطل االفتراضات بعض استخدام تمولكن . الساحلية المفتوحة المياه في الكبيرة طحالبلل depeباستخدام مميزات الماتالب ) .فوسفورالو نيتروجينال ومخرجات مدخالت تاحتسب بينما كاملة، دفعة شكل على النظام يعمل الكبيرة، طحالبلل الحيوية MATLAB دوال تتألف من دالة الحل أربع إلى البرنامج تقسيم وتم. معادالتة النمذج ومحاكاةفي نموذج الحل المنتهية العناصر طريقة تطبيق( تم تم اختبار النموذج، صحة من لتحقق. لالحدية لشروطوأخرى ل األولية، لشروطل دالةو الجزئية، التفاضلية المعادالت على تحتويودالة النتائج، تجهيزو معيارية طرقها بتحليلفي بندر عباس وتم المفتوحة الساحليةبالتركيز على المياه الكبيرة طحالبلل الحيوية والكتلة والفوسفور النيتروجين ات لكل منقياسال خالل من الساحلية الكبيرة الطحالب زراعة نظم لتطوير تطبيقه ت إمكانيةأظهر المقترح موذجللن التنبؤية القيمإن .رامتأ 11و 5 و 1 من ثالثة قاعمأ في .رعاالمز ظامن على والفوسفور النيتروجين استراتيجيات أثر تقييم .الجزئية التفاضلية لمعادالتل العددي الحل ،ةالساحلي المفتوحة المياه طحالب،لل المحاكاة الكبيرة، الطحالب زراعة الطحالب، نمو نمذجة :مفتاحیةالكلمات ال 1. 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 DEVELOPMENT OF A MATHEMATICAL MODEL 33 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–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–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. 2. 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. NO3inlet = 1 gl -1 , NO3 outlet = 0.4 gl -1 , PO4 inlet = 0.2 gl -1 , and PO4 outlet = 0.03 gl -1 . The initial macroalgae concentration in domain Ω = 0.17 gl -1 (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 up- taking 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 Bz| z=0 = 0 at the water bottom and a Dirichlet condition Bz| 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 . A. LABABPOUR 34 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]. 2.1 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 𝜕𝐵 𝜕𝑡 = 𝐷𝐵 𝜕2𝐵 𝜕𝑧2 + 𝑓(𝑡, 𝐵) 𝑖𝑛 𝛺 (1) 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 𝜕𝐵 𝜕𝑡 = 𝐷𝐵 𝜕2𝐵 𝜕𝑧2 + 𝜇𝐵 − 𝛽𝐵, (3) NO3 PO4 Macroalgae Assimilation Assimilation F, C NO3 Harvest Decay In-flow F, C PO4 In-flow Out-flow NO3 PO4 y z x (a) (b) DEVELOPMENT OF A MATHEMATICAL MODEL 35 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 Ω: Bz| 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: B = 0 on z = H. The solution satisfies assumed initial conditions of B (t 0 ) = B 0 = 0.03 gl -1 . 2.2 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 𝜇 = 𝜇𝑚𝑎𝑥 𝑁 𝑘𝑁 + 𝑁 × 𝑃 𝑘𝑃 + 𝑃 (4) 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(𝑁) − 𝜇𝑌𝐵 𝑖𝑛 𝛺. (5) 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 = Sm (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: 𝜕𝑃 𝜕𝑡 = 𝐷𝑃 𝜕2𝑃 𝜕𝑧2 + (𝑃)𝑖𝑛 + 𝑑2(𝑃) − 𝜇𝑌𝐵 𝑖𝑛 𝛺, (6) The top boundary condition is P|z = L Z = Sm (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 Nzz and C (P) = k Pzz. 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, 𝜕𝐵 𝜕𝑡 = 𝐷𝐵 𝜕2𝐵 𝜕𝑧2 + 𝐶0, (7) where the constant DB is biomass diffusivity, and the constant C0 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. 𝜕𝐵 𝜕𝑡 = D𝐵 𝜕2𝐵 𝜕𝑧2 + 𝜇𝐵 − 𝛽𝐵, 𝜇 = 𝜇𝑚𝑎𝑥 𝑁 𝑘𝑁 + 𝑁 × 𝑃 𝑘𝑃 + 𝑃 0.0 ≤ 𝜇 ≤ 1.0, 𝐵(𝑧, 𝑡0) = 𝐵(𝑧, 0) = 0.17, 𝐵(𝐿𝑧0 , 𝑡) = 𝐵(0, 𝑡) = 0.17, 𝐵(𝐿𝑧, 𝑡) = 𝐵(𝐻, 𝑡) = 0.17, 𝒜1: 𝐵(0, 𝑧) = 𝐵(𝐿𝑥 , 𝑧) ∀𝑧 ∈ (0, 𝐿𝑧 ) 𝒜2: 𝐵(0, 𝑧) = 𝐵(𝐿𝑦 , 𝑧) ∀𝑧 ∈ (0, 𝐿𝑧 ) 𝜕𝑁 𝜕𝑡 = 𝐷𝑁 𝜕2𝑁 𝜕𝑧2 + (𝑁)𝑖𝑛 + 𝑑1(𝑁) − 𝜇𝑌𝐵, (8) A. LABABPOUR 36 𝑁(𝑧, 𝑡0) = 𝑁(𝑧, 0) = 0.7 𝑁(𝐿𝑧0 , 𝑡) = 𝑁(0, 𝑡) = 0 𝑁(𝐿𝑧 , 𝑡) = 𝑁(𝐻, 𝑡) = 1 𝒜3: 𝑁(0, 𝑧) = 𝑁(𝐿𝑥 , 𝑧) ∀𝑧 ∈ (0, 𝐿𝑧 ) 𝒜4: 𝑁(0, 𝑧) = 𝑁(𝐿𝑦 , 𝑧) ∀𝑧 ∈ (0, 𝐿𝑧 ) 𝜕𝑃 𝜕𝑡 = 𝐷𝑃 𝜕2𝑃 𝜕𝑧2 + (𝑃)𝑖𝑛 + 𝑑2(𝑃) − 𝜇𝑌𝐵, 𝑃(𝑧, 𝑡0) = 𝑃(𝑧, 0) = 0.12, 𝑃(𝐿𝑧0 , 𝑡) = 𝑃(0, 𝑡) = 0, 𝑃(𝐿𝑧 , 𝑡) = 𝑃(𝐻, 𝑡) = 0.7, 𝒜5: 𝑃(0, 𝑧) = 𝑃(𝐿𝑥 , 𝑧) ∀𝑧 ∈ (0, 𝐿𝑧 ) 𝒜6: 𝑃(0, 𝑧) = 𝑃(𝐿𝑦 , 𝑧) ∀𝑧 ∈ (0, 𝐿𝑧 ) These model equations, namely equations (8), are nonlinear and have no simple exact solution. We note that the unknowns of the model are: B, N and P. The equations for B, N and P are the biomass balance, together with their initial and boundary conditions, and N and P substrates balance together with their initial and boundary conditions, respectively. 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 d1, decomposition of biomass based phosphorous compounds d2, 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. Table 1. Model parameters used for simulation of Gracilariopsis persica. Parameter Symbol Unit Value Ref. Maximum growth rate µmax d -1 0.2 [24] Death rate β d -1 0.04 [25] Biomass based nitrogen compounds d1 d -1 0.023 [25] Biomass based phosphorous compounds d2 d -1 0.01 [25] Yield Y - 0.32 [25] Michaelis constant for nitrogen k N µM 61.5 [26] Michaelis constant for phosphorous k P µM 37 [26] Nitrogen (initial) CN mg l -1 0.7 [25] Phosphorous (initial) CP mg l -1 0.1 [25] 2.3 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. 3. 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. DEVELOPMENT OF A MATHEMATICAL MODEL 37 3.1 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. Figure 2. Estimated concentrations of macroalgae biomass influenced by water depth. 3.2 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). Figure 3. Estimated concentrations of (a) nitrogen and (b) phosphorous with water depth. 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. a b A. LABABPOUR 38 Figure 4. Measured and model simulated (a) biomass, (b) nitrogen and (b) phosphorous in three water depths. Experimental (■) and fitted (-). 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–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 DEVELOPMENT OF A MATHEMATICAL MODEL 39 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. 4. 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. Acknowledgement The author acknowledges the support of the Shohadaye Hoveizeh University of Technology, Iran. I also sincerely acknowledge the anonymous and respected reviewers for their constructive recommendations. Appendix Model equations 𝜕𝑢1 𝜕𝑡 = 𝐷𝑏 𝜕2𝑢1 𝜕𝑥 2 − 𝜇𝑢1 + 𝛽𝑢1, IC : 𝑢1(𝑥, 0) = 0.1, Left – BC : 𝑢1(0, 𝑡) = 0.1, Right – BC : 𝑢1(𝐿, 𝑡) = 0.1. 𝜕𝑢2 𝜕𝑡 = 𝐷𝑁 𝜕2𝑢2 𝜕𝑥 2 − 𝑌1 ∗ 𝜇𝑢2, IC : 𝑢2(𝑥, 0) = 0.1, Left – BC : 𝑢2(0, 𝑡) = 0.1, Right – BC : 𝑢2(𝐿, 𝑡) = 0.1. 𝜕𝑢3 𝜕𝑡 = 𝐷𝑃 𝜕2𝑢3 𝜕𝑥 2 − 𝑌2 ∗ 𝜇𝑢3, IC : 𝑢3(𝑥, 0) = 0.1, Left – BC : 𝑢3(0, 𝑡) = 0.1, Right – BC : 𝑢3(𝐿, 𝑡) = 0.1. A. LABABPOUR 40 Program code and outputs global mu beta DB DN DP y1 y2; mu = 0.4; beta = 0.1; DB = 0; DN = 1e-5; DP = 1e-6; y1 = 0.3; %nitrogen yield y2 = 0.014; %phosphorous yield L = 10; %Length of domain maxt = 10; %Maximum simulation time m = 0; t = linspace(0,maxt,10); %tspan x = linspace(0,L,5); %xmesh sol = pdepe(m,@uu,@uuic,@uubc,x,t); u1 = sol(:,:,1); u2 = sol(:,:,2); u3 = sol(:,:,3); figure (1) surf(x,t,u1) mesh (x,t,u1) title('Biomass') xlabel('Depth, m','fontsize',12,'fontweight','b','fontname','arial') ylabel('Time, d','fontsize',12,'fontweight','b','fontname','arial') zlabel('Biomass, mg/ l','fontsize',12,'fontweight','b','fontname','arial') figure (2) surf(x,t,u2) mesh (x,t,u2) title('Nitrogen') xlabel('Depth, m','fontsize',12,'fontweight','b','fontname','arial') ylabel('Time, d','fontsize',12,'fontweight','b','fontname','arial') zlabel('Nitrogen, mg/ l','fontsize',12,'fontweight','b','fontname','arial') figure (3) surf(x,t,u3) meshc (x,t,u3) title('Phosphorous') xlabel('Depth, m','fontsize',12,'fontweight','b','fontname','arial') ylabel('Time, d','fontsize',12,'fontweight','b','fontname','arial') zlabel('Phosphorous, mg/ l','fontsize',12,'fontweight','b','fontname','arial') function [c,f,s] = uu(x,t,u,dudx) global mu beta DB DN DP y1 y2; u = zeros (3,1); c = [1;1;1]; f = [DB;DN;DP].*dudx; s = [-mu*u(1) + beta*u(1); -mu*y1*u(2); -mu*y2*u(3)]; end function u0 = uuic(x) u0 = [0.05; 0.7; 0.13]; end DEVELOPMENT OF A MATHEMATICAL MODEL 41 function [pl,ql,pr,qr] = uubc(xl,ul,xr,ur,t) pl = [ul(1); ul(2); ul(3)]; ql = [1.1; 1.1; 1.0]; pr = [ur(1); ur(2); ur(3)]; qr = [1;1;1]; end References 1. Negreanu-pîrjol, B., Negreanu-pîrjol, T., Paraschiv, G., Bratu, M., Sîrbu, R., Roncea, F. and Meghea, A. Physical- chemical characterization of some green and red macrophyte algae from the Romanian Black Sea littoral. Scientific Study and Research Chemistry and Chemical Engineerin, Biotechnology, Food Industry, 2011, 12 ,173– 184. 2. Israel, A., Martinez-Goss, M. and Friedlander, M. Effect of salinity and pH on growth and agar yield of Gracilaria tenuistipitata var. liui in laboratory and outdoor cultivation. Journal of Applied Phycology, 1999,11,543–549. 3. Abreu, MH., Pereira R., Yarish, C., Buschmann, A.H. and Sousa-Pinto, I. IMTA with Gracilaria vermiculophylla: Productivity and nutrient removal performance of the seaweed in a land-based pilot scale system. Aquaculture, 2011, 312,77–87. 4. Wegeberg, S. and Felby, C. Algae Biomass for Bioenergy in Denmark Biological/Technical Challenges and Opportunities. University of Copenhagen, 2010. 5. Handå, A., Forbord, S., Wang, X., Broch, OJ., Dahle, SW., Størseth, TR., Reitan, KI., Olsen, Y. and Skjermo, J. Seasonal- and depth-dependent growth of cultivated kelp (Saccharina latissima) in close proximity to salmon (Salmo salar) aquaculture in Norway. Aquaculture, 2013,414–415,191–201. 6. Bajpai, R., Prokop, A. and Zappi, M. Algal Biorefineries. Dordrecht: Springer Netherlands; 2014. 7. Perrot, T., Rossi, N., Ménesguen, A. and Dumas, F. Modelling green macroalgal blooms on the coasts of Brittany, France to enhance water quality management. Journal of Marine Systems, 2014,132,38–53. 8. Fennel, W. and Neumann, T. Coupling biology and oceanography in models. AMBIO A Journal of Human Environment, 2001,30,232–6. 9. Van Dam, AA. Modelling Studies of Fish Production in Integrated Agriculture-Aquaculture Systems. Wageningen Agricultural University, 1995. 10. Lorenzen, K., Struve, J. and Cowan, VJ. Impact of farming intensity and water management on nitrogen dynamics in intensive pond culture: a mathematical model applied to Thai commercial shrimp farms. Aquaculture Research, 1997, 28, 493–507. 11. Coppens, J., Hejzlar, J., Šorf, M., Jeppesen, E., Erdoğan, Ş., Scharfenberger, U., Mahdy. A, Nõges, P., Tuvikene A., Baho, DL., Trigal, C., Papastergiadou, E., Stefanidis, K., Olsen, S. and Beklioğlu, M. The influence of nutrient loading, climate and water depth on nitrogen and phosphorus loss in shallow lakes: a pan-European mesocosm experiment. Hydrobiologia, 2016,778,13–32. 12. Trancoso, ARR. Modelling Macroalgae in Estuaries. University of Lisbon, 2002. 13. Mocenni, C. and Vicino, A. Modelling Ecological Competition Between Seaweed and Seagrass: a Case Study 2006, 732-737. 14. Ménesguen, A., Cugier, P. and Leblond, I. A new numerical technique for tracking chemical species in a multisource, coastal ecosystem applied to nitrogen causing Ulva blooms in the bay of Brest. Limnol. Oceanogr, 2006,51,591–601. 15. Lazure, P. and Dumas, F. An external–internal mode coupling for a 3D hydrodynamical model for applications at regional scale (MARS). Advances in Water Resources, 2008, 31,233–50. 16. Vanhoutte-Brunier, A., Fernand, L., Ménesguen, A., Lyons, S., Gohin, F. and Cugier, P. Modelling the Karenia mikimotoi bloom that occurred in the western English Channel during summer 2003. Ecological Modelling, 2008, 210,351–376. 17. Klapper, I. and Dockery, J. Mathematical description of microbial biofilms. Society for Industrial Applied Mathematics Review, 2010,52,221–265. 18. Kim, Y., Park, J., Giokas, DL. and Albanis, TA. Performance evaluation and mathematical modeling of nitrogen reduction in waste stabilization ponds in conjunction with other treatment systems. Journal of Environmental Science Health Part A, 2004,39,741–758. 19. Chapelle, A., Ménesguen, A., Deslous-Paoli, J-M., Souchu, P., Mazouni, N., Vaquer, A. and Millet, B. Modelling nitrogen, primary production and oxygen in a Mediterranean lagoon. Impact of oysters farming and inputs from the watershed. Ecological Modelling, 2000,127,161–181. 20. Fernández, I., Acién, FG., Fernández, JM., Guzmán, JL., Magán, JJ. and Berenguel, M. Dynamic model of microalgal production in tubular photobioreactors. Bioresour Technology, 2012,126,172–181. 21. Geider, RJ., Maclntyre, HL. and Kana, TM. A dynamic regulatory model of phytoplanktonic acclimation to light, A. LABABPOUR 42 nutrients, and temperature. Limnology Oceanography, 1998,43,679–694. 22. Stewart, PS., Hamilton, MA., Goldstein, BR. and Schneider, T.B. Modeling biocide action against biofilms. Biotechnology Bioengineering, 1996,49,4445–4455. 23. Gjorgjieva, J. Turing Pattern Dynamics for Spatiotemporal Models with Growth and Curvature. Harvey Mudd, 2006. 24. Smit, AJ. Nitrogen Uptake by Gracilaria gracilis (Rhodophyta): Adaptations to a Temporally Variable Nitrogen Environment. Botanica Marina, 2002,45,196-209. 25. Fatemeh, E. Agar Production by Macroalga Gracilariaopsis persica in the Coastal Waters of the Persian Gulf. Journal of Applied Environmental and Biology Sciences, 2014,4,45–52. 26. Wang, C., Lei, A., Zhou, K., Hu, Z., Hao, W. and Yang, J. Growth and Nitrogen Uptake Characteristics Reveal Outbreak Mechanism of the Opportunistic Macroalga Gracilaria tenuistipitata. PLoS One, 2014,9,e108980. 27. Buschmann, AH., Varela, DA., Hernández-González, MC. and Huovinen, P. Opportunities and challenges for the development of an integrated seaweed-based aquaculture activity in Chile: determining the physiological capabilities of Macrocystis and Gracilaria as biofilters. Journal of Applied Phycology, 2008,20,571–577. 28. Thomann, RV. and Fitzpatrick, JJ. Calibration and verification of a mathematical model of the eutrophication of the Potomac Estuary. Government of the District of Columbia, Washington, D.C.: 1982. 29. Mitchell, D a., von Meien, OF., Krieger, N. and Dalsenter, FDH. A review of recent developments in modeling of microbial growth kinetics and intraparticle phenomena in solid-state fermentation. Biochemical Engineering Journal, 2004,17,15–26. 30. Alpkvist, E., Picioreanu, C., van, Loosdrecht, MCM. and Heyden, A. Three-dimensional biofilm model with individual cells and continuum EPS matrix. Biotechnology and Bioengineering, 2006,94,961–979. 31. Lababpour, A. Open-water cultivation of seaweed genus gracilaria in the coastal waters of Qeshm island for agar production. Acta Horticulturae, 2014,1054,325–332. 32. Mwegoha, WJS., Kaseva, ME., Sabai, SMM. Mathematical modeling of dissolved oxygen in fish ponds. African Journal of Environmental Science and Technology, 2010,4,625–638. 33. James, SC. and Boriah, V. Modeling Algae Growth in an Open-Channel Raceway. Journal of Computational Biology, 2010,17,895–906. 34. Dampin, N., Tarnchalanukit, W., Chunkao, K. and Maleewong, M. Fish Growth Model for Nile Tilapia (Oreochromis niloticus) in Wastewater Oxidation Pond, Thailand. Procedia Environmental Sciences, 2012,13,513–524. Received 6 June 2017 Accepted 23 August 2017