Modelling of Hydrobiological Processes in Coastal Waters

In general, marine hydrodynamics is related to many different physical, chemical, geological and biological processes, which often constitute important environmental problems that need special care and detailed investigation. Among the most frequent phenomena of the afore‐ mentioned environmental issues, observed mainly in coastal basins, are the algal blooms, which are strongly related to optical, aesthetic and other important disturbances, such as eutrophication processes. The appearance of such phenomena is becoming more and more important, because algal blooms may contain potential toxic populations, called harmful algal blooms (HAB). As such episodes may correlated with eddies development [1], their study may obviously furnish a potentially effective prediction tool. The next important step in this rather preliminary process of tracing the potential locations of harmful algal blooms is the develop‐ ment and application of a transport model, which computes the matter transfer in space and time. A further step is the integration of a biological model into a transport model. It is expected that mathematical models that include physicochemical and biological processes will be able to describe the marine hydrodynamics and the relevant matter transfer issues, such as episodes of the spread of phytoplankton cells, contributing to the better understanding and more effective investigation and management of similar phenomena. Such models have been applied and discussed in literature [2-6].


Introduction
It is well known that the dispersion of organic and inorganic matter in a coastal basin is closely related to the circulation of the seawater masses. The more detailed, accurate and reliable the knowledge of the hydrodynamic circulation is, the more successful the initial tracing of matter distribution in a coastal basin is expected to be.
In general, marine hydrodynamics is related to many different physical, chemical, geological and biological processes, which often constitute important environmental problems that need special care and detailed investigation. Among the most frequent phenomena of the aforementioned environmental issues, observed mainly in coastal basins, are the algal blooms, which are strongly related to optical, aesthetic and other important disturbances, such as eutrophication processes. The appearance of such phenomena is becoming more and more important, because algal blooms may contain potential toxic populations, called harmful algal blooms (HAB). As such episodes may correlated with eddies development [1], their study may obviously furnish a potentially effective prediction tool. The next important step in this rather preliminary process of tracing the potential locations of harmful algal blooms is the development and application of a transport model, which computes the matter transfer in space and time. A further step is the integration of a biological model into a transport model. It is expected that mathematical models that include physicochemical and biological processes will be able to describe the marine hydrodynamics and the relevant matter transfer issues, such as episodes of the spread of phytoplankton cells, contributing to the better understanding and more effective investigation and management of similar phenomena. Such models have been applied and discussed in literature [2][3][4][5][6].
The models to be used depend on the available data, the specific scientific questions to be answered and the particular characteristics of the ecosystem. The mathematical models used for hydrobiological processes in coastal waters can be focused on diagnostic matters, like the trophic state and the water quality of a specific coastal zone, or on prognostic matters, like the dispersion of HAB in a coastal basin.
Consequently, the use of mathematical simulation is of great significance for the diagnosis as well as the prognosis and prevention of hazardous situations that may threaten the marine environment. From this point of view, the description and application of three different types of models, presented in the following sections, constitutes the objective of the present chapter.

Description of the method
The trophic index (TRIX) proposed by Volleinweider et al. [7] is estimated as a linear combination of the logarithms of four variables, referring to primary production [chlorophyll-a (Chla), dissolved inorganic nitrogen (DIN), total phosphorus (TP) and the absolute percentage deviation from oxygen saturation (aD%O)], in the form: Using data from the Northern Adriatic Sea, Volleinweider et al. [7] suggested the following values for the parameters k and m: k=1.5 (to fix the differences between the upper and lower range of the variables to 3 log units) and m=1.2 (to fix the scale range from 0 to 10). Low TRIX values (2)(3)(4) indicate poor productivity in waters, and thus imply "high quality" waters, while TRIX values from 4 to 5 indicate "good quality", 5 to 6, "mediocre quality" and 6 to 8, "bad quality" waters [8].
In this form, TRIX has been applied for the classification of coastal waters in the Mediterranean Sea, the Marmara Sea, the Northern European seas, the Black Sea, the Caspian Sea, the Baltic coastal waters, the seas of Southeast Mexico, and the Persian Gulf [9][10][11][12][13][14][15][16][17]. However, when the concentrations of the variables appear with different upper and lower ranges in respect to those of the Northern Adriatic Sea, the TRIX index needs modification by re-estimating the k and m parameters according to the new upper and lower limits [7,18].
Two classification procedures for water were proposed. The first is based on a simple comparison between box and whisker plots for UNTRIX data from: 1) the investigation site and 2) the respective data from the reference area. The second procedure is based on the TQR TRIX trophic index, defined according to equation (3) where 50 th is the median and 75 th percentile is a value indicating that 75 % of the values of the variable fall below that value.

Case study
The three abovementioned trophic indices (TRIX, UNTRIX and modified TRIX) were applied to quantify the water quality of the coastal waters in the Kalamitsi area, Ionian Sea, Greece ( Fig. 1), before and after the operation of the Waste Water Treatment Plant (WWTP) [11]. As mentioned, this quantification allows a classification of the water quality. Before the WWPT operation, the TRIX values showed a temporal variation in the range of 1.9-4 TRIX units, while no spatial variation (between stations) was observed. The mean TRIX value w 3.5. After the WWPT operation, the TRIX values showed temporal and spatial variations (2.2-4 TRIX units). The mean TRIX value was 3.4. As can be seen, no significant differences in TRI values were detected before and after the WWPT operation (T-test, p=0.608). These values indica "high" water quality, in contrast with the TRIX values obtained for the gulfs of the Aegean Sea, lik the Thermaikos Gulf (5.0-6.0) and the Saronikos Gulf (3.7-6.2), with a mean of 5.3, indicatin "mediocre" water quality [14].
As the concentrations of Chl-a and DIN from Kalamitsi displayed a profile much lower than those the Northern Adriatic Sea [11,7], the TRIX was re-estimated (KALTRIX) according to the Kalamit data set, following this procedure: -zero values of the variables were substituted by the detection limit value of the applied analytic method -extreme values (m ± 2.5 SD) were excluded (where m is the mean value and SD is the standa deviation of the sample) -the range between upper and lower values of each variable was set to 1 log unit Before the WWPT operation, the TRIX values showed a temporal variation in the range of 1.9-4.7 TRIX units, while no spatial variation (between stations) was observed. The mean TRIX value was 3.5. After the WWPT operation, the TRIX values showed temporal and spatial variations (2.2-4.8 TRIX units). The mean TRIX value was 3.4. As can be seen, no significant differences in TRIX values were detected before and after the WWPT operation (T-test, p=0.608). These values indicate "high" water quality, in contrast with the TRIX values obtained for the gulfs of the Aegean Sea, like the Thermaikos Gulf (5.0-6.0) and the Saronikos Gulf (3.7-6.2), with a mean of 5.3, indicating "mediocre" water quality [14].
As the concentrations of Chl-a and DIN from Kalamitsi displayed a profile much lower than those of the Northern Adriatic Sea [11,7], the TRIX was re-estimated (KALTRIX) according to the Kalamitsi data set, following this procedure: • zero values of the variables were substituted by the detection limit value of the applied analytical method • extreme values (m ± 2.5 SD) were excluded (where m is the mean value and SD is the standard deviation of the sample) • the range between upper and lower values of each variable was set to 1 log unit • test of normality was performed for log transformations.
This procedure furnished equation (4): KALTRIX values displayed high temporal variability (0.8-9 units) before the WWPT operation, whereas both temporal and spatial variability occurred after the WWPT operation. In both cases, KALTRIX values were higher than those of TRIX, although the same data set was applied. The mean KALTRIX value was 5.3 before WWPT operation and 5.0 after it, with no statistically significant differences (T-test, p=0.735), indicating minor influence on water quality of WWTP operation. However, greater divergences of KALTRIX values were observed after WWTP operation. So, it could be assumed that KALTRIX is more sensitive, and is able to discriminate minor trends not easily detected by TRIX.
The same data set from Kalamitsi and data from South Evoikos used as a reference site (which has been reported as a having high-good ecological quality status, directive 2000/60/EC) were applied to the UNTRIX index (equation 2). Using the box and whisker plots ( According to the TQR TRIX index (equation 3), water quality is classified as "moderate" both before (0.56) and after (0.59) the WWTP operation. On the other hand, according to the boxplot procedure [10], water quality is classified as "moderate" before the WWTP operation (75 th of the site = 3.142>3.097=100 th of the reference), while after the WWTP operation it is classified as "good" [75 th of the Kalamitsi =2.994>2.118=75 th of the reference (S. Evoikos) and simultaneously 100 th of the Kalamitsi = 4.063>3.097=100 th of the reference (S. Evoikos)].

Description of the method
If an investigation is focused not only on the classification of the coastal waters' quality but also on the assessment of critical thresholds, then piecewise linear regression models with breakpoints may constitute useful tools [19]. Furthermore, piecewise linear regression models with breakpoints are able to express discontinuity appearing in ecosystems, associated either with habitat fragmentations (usually independent variables) or with "decisions" taken by individuals according to their situation [20,21]. The breakpoint (bt) is a value of either an independent variable (x) or a dependent variable (y) according to the specific characteristics of the ecosystems and the interests of the investigation, where two separate linear regression equations are estimated: one with y values less or equal to bt and the other one with y values greater than bt. The general form is given by equation 5: Such models have been used, for example, to estimate the production of crops, to search for thresholds in coastal mangrove forests and to examine a limpet feeding rate and its response to temperature [22][23][24].

Case study
As chlorophyll-a and nutrients are two of the most common indicators for eutrophication assessment [25], a piecewise linear regression model was applied with data (275 observations) for chlorophyll-a and nitrates (NO 3 -N) from the Thermaikos Gulf, Greece, which were collected over a period of two years [20].
Nitrogen was chosen because it is used by cells for the biosynthesis of molecules like chlorophyll-a, DNA, RNA and proteins. Thus, it could be assumed that discontinuity might have appeared because each phytoplankton individual "decides" the way that a given amount of nitrogen will be used, depending on the stage of its cell cycle and cell demands. In the fitted model, the concentration of dissolved nitrogen (NO 3 -N, μg/l) is the independent variable, while the concentration of chlorophyll-a (mg/m 3 ) is the dependent one.
Because the above model is piecewise linear, the parameters were estimated using the quasi-Newton method, an interactive method that minimizes the least squares loss function [∑ (observed-predicted) 2 ] through iterative convergence of the predefined empirical equation.
According to the present data, the breakpoint that minimizes the loss function and simultaneously gives the highest correlation coefficient (r=0.84) is 2.5 mg/m 3 of chlorophyll-a. So the model of equation ( Plotting the predicted and observed values [20], it was found that the predicted values of Chla were overestimated when observed values were lower than 0.9 (mg/m 3 ), and underestimated when observed values were higher than 4.4 (mg/m 3 ). Thus, we simulated the model starting with a lower breakpoint value (0.5 mg/m 3 Chl-a) and terminated it at breakpoint = 4.8 mg/m 3 .
The results show that when breakpoint values are lower than 0.9 or higher than 4.4 mg/m 3 Chla, no relationship could be estimated between Chl-a and NO 3 . However, when 0.9≤break-point≤4.4 mg/m 3 Chl-a, a relationship between both parameters appeared. The negative slopes b i (-0.08411 and -0.16438 respectively) indicate a "tendency of the ecosystem to be stabilized" at the values of 2.5 or 0.9 mg/m 3 of Chl-a, when its concentration in the Thermaikos Gulf varies between 2.5 and 4.4, and 0.9 and 2.5, respectively. It can be assumed that, when Chl-a values are higher than 2.5, then competition and predation may tend to stabilize the phytoplankton community, while processes like decomposition and inflows may increase NO 3 concentration at the same time.
The estimated critical values of breakpoints (0.9, 2.5 and 4.4) are very close to the criteria for trophic assessment in marine and coastal ecosystems given by the Swedish Environmental Protection Agency (<1.5, 1.5-2.2, 2.2-3.2, 3.2-5 and >5 μg/l Chl-a) and the low risk "trigger values" (2-4 μg/l Chl-a) which were proposed by the Australian and New Zealand Environment and Conservation Council [26]. So, it could be assumed that, in the case of the Thermaikos Gulf, values of Chl-a lower than 0.9 may indicate oligotrophic waters, from 0.9 to 2.5, mesotrophic, from 2.5 to 4.4, eutrophic, and above 4.4, hypertrophic waters. Furthermore, it has been reported that Chl-a concentrations between 2 and 4.44 appeared during algal blooms [27,28], so it could be suggested that values of Chl-a close to the breakpoint may indicate the start of an algal bloom, although more field and laboratory investigations are needed to confirm this hypothesis.

Description of the method
In this part of the chapter, a hydrobiological/bio-hydrodynamic model is described. This model is constituted by two parts: a first hydromechanical part, which computes the hydrodynamics and the matter transfer in a coastal geophysical basin; and a second hydrobiological part, which computes the cells' growth rate (generation of new mass) and the cells' decay (loss of mass). When shallow waters characterize a coastal basin, a two-dimensional, depth-averaged, hydrodynamic model is considered to be quite sufficient for the simulation of the seawater circulation. In this case, the hydrobiological model is also a depth-averaged model, like the one presented below.
As mentioned above, the first part of the model refers to the hydromechanical processes. These hydromechanical processes concern (a) the hydrodynamics, e.g., the computation of the velocity field and the sea surface elevation field; and (b) the matter transport, described by the advection and dispersion processes.

Hydrodynamics.
The hydrodynamic-hydromechanical model is based on the usual equations for the conservation of mass and momentum, as described in equations (7) where h is the seawater depth, U and V are the vertically averaged horizontal velocities, ζ is the sea surface elevation, f is the Coriolis parameter, τ sx and τ sy are the wind surface shear stresses, τ bx and τ by are the bottom shear stresses, ν h is the dispersion coefficient related to the local vorticity, ρ is the density of the water and g is the gravity acceleration [29,31].
The output from the hydrodynamic simulation, e.g., the velocity field, is then used as an input in the transport simulation. The specific transport model applied here is based on the tracer method which can be described as a Lagrange-Monte Carlo simulation (random walk method) [29,3]. The simulation of the transport processes based on this method has the important advantage of avoiding numerical problems. Advection and dispersion processes are simulated by means of this method [32][33][34][35][36]. In this specific work, the two parts of the hydromechanical model, i.e., the hydrodynamic and transport processes, are fully coupled, including the biological processes at the same time. In more detail, the transport part of the model uses, as already mentioned, the velocity field, produced by the hydrodynamic model at each time step. The transport model is described analytically in the following paragraphs.
Matter transport. According to the tracer method, a large number of particles that represents a particular amount of mass are released in the marine environment through a source. The location where an algal bloom episode occurs determines the position of that source. The transport and fate of each particle is traced with time. Advection of the particulate matter is computed by the local seawater velocity while turbulent diffusion is simulated by the random Brownian motion of the particles. More specifically, the tracer method is described by the following steps [29,6]: a. the velocity field is determined by a set of values of the velocity components at specific grid points computed by the hydrodynamic part of the model b. a suitable time step is selected c. the range of the random velocity ± Ur is computed from the following equation (8): where D is the local diffusion coefficient and dt is the time step d. in the case of instantaneous discharge, a specific amount of particles is released at once from the initial source; while in the case of continuous discharge, for each time step a specific number of particles is introduced in the study domain e. integration in time is executed and the new coordinates of each particle are computed; the motion of each particle is analysed considering: i) a deterministic part which concerns the advective transport, and ii) a stochastic part which concerns the transport-spreading due to diffusion processes.
The horizontal positions of the particles are computed from the superposition of the deterministic and stochastic displacements:

x t dt and x u dt rnd
where Δx i n and Δy i n are the deterministic displacements and Δx i n ΄ and Δy i n ΄ are the stochastic displacements, with u i n (x i n , t n ) and v i n (y i n , t n ) being the deterministic velocities at time t n at the location x i n and y i n of the i particle, u i n ΄& v i n ΄are the random (stochastic) horizontal velocities at time t n at the locations x i and y i , respectively, u i n ΄= v i n ΄ = √(6D h /dt), with D h being the horizontal particle diffusion coefficient, and rnd is a random variable distributed uniformly between -1 and +1. Eventually, the spatial particle distribution allows for the computation of the particle concentrations in each grid box.

If we define
• N 0 as the starting population of a phytoplanktonic species • μ as the growth rate (divisions/day) and • TL as the total losses then these three factors regulate the density of population N, which, in a continuous model, follows an exponential growth, according to equation (11): where t is the time (in days), TL is influenced by various processes such as predation, sinking, water movement and stratification, and μ is dependent on abiotic factors such as temperature, light and nutrients, since they act directly on the individuals' biochemistry [37,38].
In the biological part of the model, the maximum growth rate of phytoplankton μ (T) is considered to follow Eppley's equation [39,40]. Then light intensity and sun shine duration are considered to limit μ (T) according to Steele's equation [41]. The sunshine duration and the photosynthetic active radiation (par) were estimated using field data [42,43]. Nitrogen and phosphorus limitations (N lim and P lim respectively) are computed according to Michaelis-Menten enzyme Kinetics. Temperature and nutrients concentrations are considered to take random values, following normal distribution which was estimated from field data [27].
The abovementioned processes, with the relevant factors as functions of different elements and the relevant constants, are given in Tables 2a and 2b. In more detail, the information for the growth rate coefficients was introduced into the model as a time series of mean daily growth rates. That is, it was computed in a discrete way, instead of using the continuous equation (11). The modelling approach, concerning the increase of the particles due to division and the decrease due to particle losses, is described by the following steps [6]: a. After a time period of one day, a new number of particles equal to (μ•N 0 ) are generated.
In that product (i.e. μ•N 0 ) μ is the growth rate coefficient and N 0 is the number of particles before the division, denoting the total number of particles of the previous day. The "new" total number of particles N, after the growth process, is afterward computed from the equation N=N 0 + μ•N 0 .
b. The new particles' position is determined from the positions of other particles randomly selected.
c. Additionally, a decay coefficient (TL) is adopted for the population decrease and particle disappearance of the water column. So, at the end of a period of a day, a number of particles equal to TL•N is drawn out, where TL is the aforementioned decay coefficient and N is the total number of particles of the previous day, plus the particles added due to the growth process of the present day, e.g., TL•N = TL•(N 0 +μ•N 0 ). These particles are no longer taken into account in the computations.
d. Finally, the total number of particles Nt, after growth and decay process, is given by eq. 13: Concerning the boundary conditions for the transport model, it should be pointed out that if the particles reach the coast, they return to their previous position, while the particles reaching the open sea boundary are trapped there and no further computation is made for them; these particles are then excluded from the computational loops. The position of the source is determined from the place where an episode of an algal bloom is observed. Particles' concentrations are then computed from the number of particles, counted at each grid box. As far as the initial conditions are concerned, the numerical simulations always start from scratch, i.e., zero current velocities at the starting point of the simulation time.

The hydrobiological model
The application of the mathematical hydrobiological model described in the previous section is presented here.
a. Application to a real episode in the Thermaikos Gulf. The model runs were based on meteorological data for a particular period in which a real episode of algal bloom took place in the Thermaikos Gulf, starting from January 10 th , 2000. The Thermaikos Gulf is a semi-enclosed coastal area, located in the north-western corner of the Aegean Sea, in the east of the Mediterranean Sea, as depicted in Figs. 2 and 3. The Axios, Aliakmon and Loudias are the three main rivers that discharge into the western and north-western coasts of the gulf, while the River Pinios outflows much further south, in the south-western coast of the outer gulf. During recent decades, there has been a major reduction in the total freshwater input to the gulf, mainly as the result of the extraction of river water for irrigation, and the construction of a series of hydroelectric power dams on the Aliakmon [55,56].
The hydrodynamic model was fully coupled with the transport and biological models described in the previous sections. The equations of the hydrodynamic model were numerically solved by the finite difference method on an orthogonal, staggered grid (Arakawa C grid), where U and V components refer to the nodes' sides and ζ refers to the interior of each mesh. The field of the gulf was discretized using 32×36 grid cells. The spatial-discretization step was dx = 1000 m. control points were used corresponding to the locations A and B, as depicted in Fig. 2. The hydrodynamic circulation of the gulf was activated by variable winds prevailing over this period (between January 10 th and 31 st , 2000, with a predominance of northern winds that correspond to the most frequent conditions in the Thermaikos Gulf). Concerning the effect of tide on the general pattern of dispersion, it is considered negligible, since the mean tidal signal of 25 cm is quite small and the residual currents due to the tide are negligible. Moreover, riverine water inputs and seawater density differences in general were not taken into account, because the total freshwater input of the rivers has been dramatically reduced in the last decades. Finally, the present numerical simulation was performed for wind-generated circulation, under the winds blowing over the area during the period between January 10 th and 31 st , 2000. The time series of wind data was recorded from a meteorological station of the Institute of Forestry Research at the location of Sani, Chalkidiki. As mentioned, the transport model was applied considering the position A as the local source of particulate matter and instantaneous appearance of algal bloom, where large concentrations of Dinophysis (HAB) were recorded in the field.
The daily variable estimated growth rate used in the biological model has a range between 0.176 and 0.489 divisions per day [6]. These values are comparable with those estimated using the long-term (one year) deterministic model for dinoflagellates in the same area [37,38]. In the literature, the values of dinoflagellates' growth rates vary within a range of 0.04-0.7 divisions per day, depending on time and local conditions [1,48,57].
The model results concerning the distribution of the harmful phytoplankton cells for time periods of five days, 10 days and 21 days after the most intense bloom of January 10 th are given in Figs. 4a, 4b and 4c [6]. According to the abovementioned patterns of the Dinophysis dispersion, it may be observed that five days after the bloom at location A (Fig. 4a), Dinophy- The present application to the real coastal basin of Thermaikos Gulf tried to simulate the fate of harmful alg sudden appearance of HAB, using Dinophysis abundances 10 4 -10 5 cells/lit recorded in north-eastern T January 10 th , 2000 (position A in Fig. 3a). For the simulation, a number of 1000 particles was initially used (A), and 619,000 particles were finally counted 21 days later (this latter number included particles t dispersed in the coastal basin of the gulf or escaped out of the gulf). The model runs led ultimately to the Dinophysis populations in the coastal basin of the gulf. Two control points were used corresponding to t and B, as depicted in Fig. 2. The hydrodynamic circulation of the gulf was activated by variable winds p this period (between January 10 th and 31 st , 2000, with a predominance of northern winds that correspon frequent conditions in the Thermaikos Gulf). Concerning the effect of tide on the general pattern of d considered negligible, since the mean tidal signal of 25 cm is quite small and the residual currents due negligible. Moreover, riverine water inputs and seawater density differences in general were not taken because the total freshwater input of the rivers has been dramatically reduced in the last decades. Final numerical simulation was performed for wind-generated circulation, under the winds blowing over the a period between January 10 th and 31 st , 2000. The time series of wind data was recorded from a meteorolog the Institute of Forestry Research at the location of Sani, Chalkidiki. As mentioned, the transport mode considering the position A as the local source of particulate matter and instantaneous appearance of algal large concentrations of Dinophysis (HAB) were recorded in the field.
The daily variable estimated growth rate used in the biological model has a range between 0.176 and 0.489 day [6]. These values are comparable with those estimated using the long-term (one year) determini dinoflagellates in the same area [37,38]. In the literature, the values of dinoflagellates' growth rates vary of 0.04-0.7 divisions per day, depending on time and local conditions [1,48,57].
The model results concerning the distribution of the harmful phytoplankton cells for time periods of five and 21 days after the most intense bloom of January 10 th are given in Figs. 4a, 4b and 4c [6]. Acc abovementioned patterns of the Dinophysis dispersion, it may be observed that five days after the bloom (Fig. 4a), Dinophysis masses were transported and concentrated in the area of Megalo Emvolo (position E of position A; 10 days after the initial bloom (Fig. 4b), higher concentrations of Dinophysis reached some c the west and north-west Thermaikos Gulf, where the largest mussel farms of Greece lie; finally, 21 days a (Fig. 4c), the present results show that concentrations of Dinophysis have been dispersed to a relatively la outer Thermaikos Gulf, and mainly along the west coast of the gulf [6]. Of course, because a large part of well as the eastern coastline consists of bathing coasts, e.g., tourist sea beaches, situations of large algal sis masses were transported and concentrated in the area of Megalo Emvolo (position E in Fig.  2) west of position A; 10 days after the initial bloom (Fig. 4b), higher concentrations of Dinophysis reached some coastal areas of the west and north-west Thermaikos Gulf, where the largest mussel farms of Greece lie; finally, 21 days after the bloom (Fig. 4c), the present results show that concentrations of Dinophysis have been dispersed to a relatively large area of the outer Thermaikos Gulf, and mainly along the west coast of the gulf [6]. Of course, because a large part of the western as well as the eastern coastline consists of bathing coasts, e.g., tourist sea beaches, situations of large algal concentrations are not only aesthetically undesirable but also hazardous on some occasions. In more detail, the simulation shows that 21 days after the bloom at the north-eastern side of the inner The with an initial concentration of 57,000 cells/lit at position A, and under the influence of prevailing northe eastern winds, the main mass was spread to the central and eastern coasts of the inner Thermaikos, as western coasts of the gulf. Concentrations higher than 200 cells/lt were computed for the region nor Aliacmon river-mouth and south of the Axios river-mouth (position Β in Fig. 2). These concentrations ar those reported by [27]. Concerning the north-western coasts of the gulf, e.g., the region of mussel culture the Axios delta  fig. 2). Along the coast of Thessaloniki City (position T, Fig. 2), significant Din were not found [6]. In more detail, the simulation shows that 21 days after the bloom at the north-eastern side of the inner Thermaikos Gulf, with an initial concentration of 57,000 cells/lit at position A, and under the influence of prevailing northern and north-eastern winds, the main mass was spread to the central and eastern coasts of the inner Thermaikos, as well as to the western coasts of the gulf. Concentrations higher than 200 cells/lt were computed for the region north-east of the Aliacmon river-mouth and south of the Axios river-mouth (position Β in Fig. 2). These concentrations are very close to those reported by [27]. Concerning the north-western coasts of the gulf, e.g., the region of mussel cultures north-east of the Axios delta (north of position B, area M1), the model resulted to values of up to 400 cells/lt. Furthermore, concerning the west and south-western coasts of the gulf, hosting the mussel cultures south of the Loudias river-mouth (mussel culture area M2, Fig. 2) and west and south of the Aliakmon river-mouth (mussel culture area M3, fig. 2), the of the model resulted to values of up to 900 cells/lt. The Dinophysis concentrations reached similar values at the northern area of Megalo Emvolo (position E, fig. 2). Along the coast of Thessaloniki City (position T, Fig. 2), significant Dinophysis masses were not found [6].

b. Application of different scenarios to the extended area of the Thermaikos basin
The application of the aforementioned hydrobiological model, based on different wind conditions prevailing over the extended area of the Thermaikos Gulf, is presented in this section. To this end, the outer part of the external Thermaikos Gulf was also included in the present simulations (Fig. 6a). Furthermore, this application was focused on the investigation of the critical conditions for the appearance of an algal bloom episode and the dispersion of its population in space and time. Additionally, different positions in the Thermaikos Gulf, as sources of the initial population, and different starting population densities were examined.
The equations of the hydrodynamic model, in both aforementioned cases, were numerically solved by the finite difference method, on an orthogonal staggered grid (Arakawa C grid) similar to the previous bathymetry (without the greater area of the outer Thermaikos Gulf). The study area was now discretized with a grid of 41×42 cells, with spatial step dx = 1852 m (Fig. 5b)  In this simulation, the following two assumptions drawn from the literature were taken into account: a. Although the cell concentration (or abundance) in the peak of the bloom exceeds the number of 10 4 ce minimum of 500 to 1200 cells/lt is a threshold for restrictions in fisheries [44]. So, in this study, 2000/lt Dinophysis spp. were considered to be the minimum population for starting a bloom episode. b. In almost all the studies, before a bloom episode, Dinophysis cell abundance is usually less than 200 ce

Hydrodynamics -Concepts and Experiments
In this simulation, the following two assumptions drawn from the literature were taken into account: a. Although the cell concentration (or abundance) in the peak of the bloom exceeds the number of 10 4 cells/lt, a minimum of 500 to 1200 cells/lt is a threshold for restrictions in fisheries [44]. So, in this study, 2000/lt cells of Dinophysis spp. were considered to be the minimum population for starting a bloom episode.

b.
In almost all the studies, before a bloom episode, Dinophysis cell abundance is usually less than 200 cells/lt. For this reason, we consider here N o = 150 cells/lt.
The initial model runs' simulation started with a low phytoplankton net growth rate of μ = 0.3 divisions per day. Finally, the simulation tests showed that HAB episodes appeared when the net growth rate takes values close to 1.0 divisions per day, while the losses were embodied in the net growth rate.
A number of 1000 particles was initially used at the location A, and 128,000 particles were finally counted at the end of the simulation. The total simulation time was taken as equal to seven days (time period of one week) which was sufficient for the flow to reach steady state conditions (either for N or for S winds, with speeds of 2 and 10 m/s). The patterns of seawater circulation, corresponding to steady state flow, under the influence of north and south winds of 10 m/s, are given indicatively in Fig. 6.  In this simulation, the following two assumptions drawn from the literature were taken into account: a. Although the cell concentration (or abundance) in the peak of the bloom exceeds the number of 10 4 cell minimum of 500 to 1200 cells/lt is a threshold for restrictions in fisheries [44]. So, in this study, 2000/lt c Dinophysis spp. were considered to be the minimum population for starting a bloom episode. b. In almost all the studies, before a bloom episode, Dinophysis cell abundance is usually less than 200 cell reason, we consider here No = 150 cells/lt. The initial model runs' simulation started with a low phytoplankton net growth rate of μ = 0.3 divisions pe the simulation tests showed that HAB episodes appeared when the net growth rate takes values close to 1.0 day, while the losses were embodied in the net growth rate.
A number of 1000 particles was initially used at the location A, and 128,000 particles were finally counted the simulation. The total simulation time was taken as equal to seven days (time period of one wee sufficient for the flow to reach steady state conditions (either for N or for S winds, with speeds of 2 and patterns of seawater circulation, corresponding to steady state flow, under the influence of north and sout m/s, are given indicatively in Fig. 6. Finally, the combination of the following cases was considered: Finally, the combination of the following cases was considered: a. two different positions for the starting point, the first one, position A, in the inner part of the gulf, and the second one, position B, in the outer gulf, as depicted in Fig. 5a, b. two different wind directions, the north (N) and the south (S) winds, which are generally the most frequent in the study area of the gulf, c. two different wind speeds, the first one corresponding to low winds of 2 m/s and the second one corresponding to relatively strong winds of 10 m/s. These scenarios led to eight different patterns of algal dispersion, as presented in the following sequence: four patterns corresponding to the case of initial algal bloom at location A, and four patterns corresponding to the case of initial algal bloom at location B (with A and B depicted in Fig. 5a).

Algal bloom at location A
More specifically, in the case of starting an algal bloom at location A in the inner part of the gulf, the patterns of the dispersion of the cells under the influence of north winds with speeds of 2 and 10 m/s, respectively, are given in Fig. 7; the patterns of the dispersion of the cells under the influence of south winds with speeds of 2 and 10 m/s, respectively, are given in Fig. 8.
The following characteristics [5] can be described: • If the source is located in the Thessaloniki Gulf (inner part of the Thermaikos Gulf), the dispersion of the bloom in the outer part of the Thermaikos Gulf is rather towards the rear • A bloom can be started with a low population (150 cells/lt) and a net growth rate close to 1.0 divisions per day, in a period of one week • South winds with a low speed of 2 m/s can cause higher cell concentrations than north winds with the same speed.
a. two different positions for the starting point, the first one, position A, in the inner part of the gulf, and one, position B, in the outer gulf, as depicted in Fig. 5a, b. two different wind directions, the north (N) and the south (S) winds, which are generally the most freq study area of the gulf, c. two different wind speeds, the first one corresponding to low winds of 2 m/s and the second one corre relatively strong winds of 10 m/s. These scenarios led to eight different patterns of algal dispersion, as presented in the following sequence corresponding to the case of initial algal bloom at location A, and four patterns corresponding to the case bloom at location B (with A and B depicted in Fig. 5a).

Dummy Text Algal bloom at location A
More specifically, in the case of starting an algal bloom at location A in the inner part of the gulf, the dispersion of the cells under the influence of north winds with speeds of 2 and 10 m/s, respectively, are g the patterns of the dispersion of the cells under the influence of south winds with speeds of 2 and 10 m/ are given in Fig. 8.
The following characteristics [5] can be described:


If the source is located in the Thessaloniki Gulf (inner part of the Thermaikos Gulf), the dispersion o the outer part of the Thermaikos Gulf is rather towards the rear

Algal bloom at location B
In the case of starting an algal bloom at location B in the outer part of the gulf, the patterns of the dispersion of the cells under the influence of north winds with wind speeds of 2 and 10 m/ s, respectively, are given in Fig. 9; the patterns of the dispersion of the cells under the influence of south winds with wind speeds of 2 and 10 m/s, respectively, are given in Fig. 10.
The following characteristics [5] can be described: • An algal bloom may appear only when wind speed is low (~2 m/s) • South winds may cause dispersion of phytoplankton cells in the inner gulf.

Dummy Text Algal bloom at location B
In the case of starting an algal bloom at location B in the outer part of the gulf, the patterns of the dispersion under the influence of north winds with wind speeds of 2 and 10 m/s, respectively, are given in Fig. 9; the patt dispersion of the cells under the influence of south winds with wind speeds of 2 and 10 m/s, respectively, ar Fig. 10.
The following characteristics [5] can be described:  An algal bloom may appear only when wind speed is low (~2 m/s)  South winds may cause dispersion of phytoplankton cells in the inner gulf.
More analytically, when the source is located at position B in the outer gulf, north winds do not seem to cause and dispersion of phytoplankton cells at the inner gulf during the period of one week. However, cells can be fo the eastern coasts of the outer gulf. In this case (source at position B), HABs were not observed in th calculations, except under the influence of low winds (~2 m/s). South winds may transfer and disperse phyt cells in the inner gulf and close to the areas of mussel culture after a period of one week, but phyt concentrations in all cases are shown to be lower than 1000 cells/lt.

Dummy Text Algal bloom at location B
In the case of starting an algal bloom at location B in the outer part of the gulf, the patterns of the dispersio under the influence of north winds with wind speeds of 2 and 10 m/s, respectively, are given in Fig. 9; the p dispersion of the cells under the influence of south winds with wind speeds of 2 and 10 m/s, respectively Fig. 10.
The following characteristics [5] can be described:  An algal bloom may appear only when wind speed is low (~2 m/s)  South winds may cause dispersion of phytoplankton cells in the inner gulf.
More analytically, when the source is located at position B in the outer gulf, north winds do not seem to ca and dispersion of phytoplankton cells at the inner gulf during the period of one week. However, cells can be the eastern coasts of the outer gulf. In this case (source at position B), HABs were not observed in calculations, except under the influence of low winds (~2 m/s). South winds may transfer and disperse p cells in the inner gulf and close to the areas of mussel culture after a period of one week, but p concentrations in all cases are shown to be lower than 1000 cells/lt. More analytically, when the source is located at position B in the outer gulf, north winds do not seem to cause transport and dispersion of phytoplankton cells at the inner gulf during the period of one week. However, cells can be found along the eastern coasts of the outer gulf. In this case (source at position B), HABs were not observed in the present calculations, except under the influence of low winds (~2 m/s). South winds may transfer and disperse phytoplankton cells in the inner gulf and close to the areas of mussel culture after a period of one week, but phytoplankton concentrations in all cases are shown to be lower than 1000 cells/lt.

Model description
As randomness and uncertainty appear in natural ecosystems, Markovian models, as a part of stochastic proc been applied in many ecological processes, including succession, population dynamics, energy flow and div [58][59][60][61]. In Markovian models, the definition of "states" is a critical point and depends on the specific characteri system to be modelled. A state consists of individuals, variables or number of species, while all the possible s system correspond to the "state space" C (C={1,2,3...,k}). The vector n(t)=[n1(t), n2(t), ..., nk(t)] gives the expect of the members on each one of the k "states", for each t (

Model description
As randomness and uncertainty appear in natural ecosystems, Markovian models, as a part of stochastic processes, have been applied in many ecological processes, including succession, population dynamics, energy flow and diversity [37,[58][59][60][61]. In Markovian models, the definition of "states" is a critical point and depends on the specific characteristics of the system to be modelled. A state consists of individuals, variables or number of species, while all the possible states of the system correspond to the "state space" C (C={1,2,3...,k}). The vector n(t)=[n 1 (t), n 2 (t),..., n k (t)] gives the expected number of the members on each one of the k "states", for each t ( Table 1). The so-called transition probability (p ij (t)) is the probability of a member moving from state i, at time t, to state j at time t+1. Transition probabilities can be estimated as: . . . n kk (t+1) The "state space" and time can be either continuous or discrete [62,63]. A Markov process is called stationary or homogeneous if P(t)=p ij , ∀ t ∈ ℕ, otherwise it is called non-homogeneous. A non-homogeneous Markov process undergoes cyclic behaviour if there exists d ∈ ℕ such that P(ad+s)=P(s), ∀ t ∈ ℕ and s={1,2,...,d-1} [64]. For ecosystems with seasonal cycles d=4, let P(0), P(1), P(2) and P(3) be the transition matrices from winter to spring, spring to summer, summer to autumn and autumn to winter, respectively. Then, the transition matrices from winter to winter (P 0 ), spring to spring (P 1 ), summer to summer (P 2 ) and autumn to autumn (P 3 ) are given by: P = P 0 P 1 P 2 P 3 P = P 1 P 2 P 3 P 0 P = P 2 P 3 P 0 P 1 P = P 3 P 0 P 1 P 2 The matrices P i , i={0,1,2,3} are stochastic and converge if they are regular, and the convergence is geometrically fast. Their limits (P i t ) are independent of the starting season, meaning that lim(winter to winter)=lim(spring to winter)=. .. =lim(autumn to winter), and give the probabilities of the system in "statistical equilibrium" for each season [65][66][67].

Case study
A study of zoobenthos dynamics from the Thermaikos Gulf, Greece, is presented in this section [37,67]. In this study, the following aspects are examined: the seasonal variation of zoobenthos richness (number of species/0.2 m 2 ), abundance (individuals/0.2 m 2 ) and diversity. Concerning zoobenthos richness, four states were chosen: {A 1 (up to 10 species), A 2 (from 11 to 15 species), A 3 (from 16 to 20 species) and A 4 (more than 20 species). The definition of the "state space" was based on the literature. For example, fewer than 10 species have been reported in pollutant waters with oxygen concentrations less than 1.4 ml/l [68,69], while in gulfs similar to the Thermaikos Gulf, without pollution effects, more than 20 species have been reported [70]. The definition of A 3 and A 4 was based on data from areas under eutrophication processes, where the number of species varied seasonally between 10 and 20 [70][71][72]. On abundance, four states were defined {(B 1 (less than 90 individuals), B 2 (from 91-140), B 3 (from 141 to 200) and B 4 (more than 200)}, following the same assumptions as before when considering diversity. The estimated matrices with transition probabilities from each season to the next, starting from winter to spring, are presented in Table 4.

Hydrodynamics -Concepts and Experiments 42
The matrices P i , i={0, 1, 2, 3} and their limits (P i t ) were estimated (Fig. 11). Additionally, in order to test if the number of individuals (abundance) and the number of species (richness) are independent, correlation coefficients were estimated and, furthermore, simple regression models were applied to each set of data separately. In all cases, the statistical tests confirm that the two variables are independent (p<0.05). So, the joint probability P(A i and B i ) in statistical equilibrium could be estimated as follows: P A and B = P A * P B i= 1,2,3,4 . ® 18 correlation coefficients were estimated and, furthermore, simple regression models were applied to each set of data separately. In all cases, the statistical tests confirm that the two variables are independent (p<0.05). So, the joint probability P(A i and B i ) in statistical equilibrium could be estimated as follows: P(A i and B i ) = P(A i ) * P(B i ) i={1,2,3,4}.
a) b) Figure 11. The limits of transition probabilities' matrices for each season: a) for species richness and b) for abundance (population density). The values of the probabilities are represented in the y axes.
According to Fig. 11, it could be said that winter is the best season, because the probabilities of the system being in the more "rich states" (A 3 , A 4 and B 3 and B 4 ) are 0.6916 and 0.6221, respectively. At the same time, the probability of the system to have more than 16 species and more than 141 individuals {(A 3 and B 3 ) or (A 3 and B 4 ) or (A 4 and B 3 ) or (A 4 and B 4 )}, in steady state, is higher than 0.42 (Fig. 12). On the other hand, during spring, the probabilities of A 1 (=0.46) and B 1 (=0.5221) Figure 11. The limits of transition probabilities' matrices for each season: a) for species richness and b) for abundance (population density). The values of the probabilities are represented in the y axes.
According to Fig. 11, it could be said that winter is the best season, because the probabilities of the system being in the more "rich states" (A 3 , A 4 and B 3 and B 4 ) are 0.6916 and 0.6221, respectively. At the same time, the probability of the system to have more than 16 species and  Fig. 12). On the other hand, during spring, the probabilities of A 1 (=0.46) and B 1 (=0.5221) taking their highest values among all the seasons (Fig. 11) and, furthermore, the probability of the system having fewer than 15 species and simultaneously less than 90 individuals, is higher than 0.6 ( Fig. 12). In relation to these values, it could be suggested that in spring, strong environmental influences tend to reduce both the species richness and the populations density. During summer, it can be observed ( Table 3) that transition probabilities to A 1 and A 2 still have high values, but by contrast transition probabilities to A 3 and A 4 are higher than the respective ones in spring. Concerning abundance, it could be observed (Fig. 11) that B 1 is dominant, with the highest probability value (=0.5285). The probability of the system having fewer than 15 species and fewer than 90 individuals is 0.44. So, during summer, population density seems to stabilize in low densities, similar to the behaviour observed in the spring. During autumn, the "richest state", A 4 , takes its highest probability value (=0.1823) and simultaneously A 3 and A 4 are higher than those of spring and summer. In addition, B 1 probability is lower than that observed for summer, although it continues to dominate (Fig. 11). The joint probability of the system having 11-15 species and up to 140 individuals is estimated at 0.353 (Fig. 12). So, it could be considered that a tendency towards recovery starts from autumn and continues to winter. According to Fig. 11, it could be said that winter is the best season, because the probabilities of the system more "rich states" (A3, A4 and B3 and B4) are 0.6916 and 0.6221, respectively. At the same time, the pro system to have more than 16 species and more than 141 individuals {(A3 and B3) or (A3 and B4) or (A4 and B4)}, in steady state, is higher than 0.42 (Fig. 12). On the other hand, during spring, the probabilities of A1 (=0.5221) taking their highest values among all the seasons (Fig. 11) and, furthermore, the probability having fewer than 15 species and simultaneously less than 90 individuals, is higher than 0.6 ( Fig. 12). In re values, it could be suggested that in spring, strong environmental influences tend to reduce both the specie the populations density. During summer, it can be observed ( Table 3) that transition probabilities to A1 an high values, but by contrast transition probabilities to A3 and A4 are higher than the respective on Concerning abundance, it could be observed (Fig. 11) that B1 is dominant, with the highest probability v The probability of the system having fewer than 15 species and fewer than 90 individuals is 0.44. So, du population density seems to stabilize in low densities, similar to the behaviour observed in the spring. Du the "richest state", A4, takes its highest probability value (=0.1823) and simultaneously A3 and A4 are highe spring and summer. In addition, B1 probability is lower than that observed for summer, although i dominate (Fig. 11). The joint probability of the system having 11-15 species and up to 140 individuals i 0.353 (Fig. 12). So, it could be considered that a tendency towards recovery starts from autumn and continu Winter Spring The above discussion points out that the composition of the benthic community of zoobenthos from the Thermaikos Gulf may be regulated by seasonal sequences of extinction and recolonization. The two "forcing variables" in the Thermaikos Gulf are oxygen concentration and temperature. So, species with tolerance to high temperature, low oxygen and high pollutant concentrations may regulate their births from spring to summer, while species with no tolerance to low oxygen or/and high temperatures may regulate their births from autumn to winter, when the waters are saturated with oxygen. The low values of abundance (population density) were observed during spring to autumn, and may be caused by mass mortalities due to oxygen depletion in the bottom waters [67].

Conclusions
The fact that the TRIX index has been applied in many areas suggests its use as a tool not only for the classification of marine waters but also for the comparison between data from different regions. The problem of its calibration arises when marine ecological processes of the study region differ from those of the Northern Adriatic Sea. In such cases, a re-scaled index like the KALTRIX seems to be more sensitive than TRIX to minor changes in the parameters, but as the index was re-scaled its results are not directly comparable with those of other regions. UNTRIX is also a sensitive index, but the problem of the reference site and its characteristics is evident. So, it could be said that if the interest of a manager is to compare a marine ecosystem before and after a strong human influence (like WWTP), then more sensitive indices like KALTRIX or UNTRIX should be used; otherwise, the use of the TRIX index gives the advantage of enabling direct comparisons between different marine ecosystems.
An alternative approach for marine water classification and management is the use of piecewise linear regression models with breakpoints. Such models are able to express discontinuities in ecosystems and to detect critical thresholds in coastal ecosystems.
Concerning the case of the numerical dynamic hydrobiological model, the dispersion of harmful cells after a HAB episode in a coastal basin was investigated. As expected, the wind speed and direction significantly affect the dispersion of an algal bloom. According to the results of the present work, the location of an initial outbreak in combination with the topography of a coastal basin play an important role in the appearance and the dispersion of an algal bloom, because winds with the same speed induce two different patterns of dispersion and different phytoplankton abundances (or concentrations) in space and time. The application of the dynamic hydrobiological model simulating a real episode showed that after a phytoplankton bloom at the north-east coasts of the Thermaikos Gulf (with variable wind forcing over the coastal basin and prevailing northern winds), concentrations of harmful algal cells reached (a) the north-west and west coasts of the gulf where the largest part of the mussel culture of Greece lies, and (b) the west coasts of the gulf with important bathing beaches, within a period of 10 days. For 20 days after the bloom, the mass of Dinophysis (HAB) dispersed in quite large areas of the outer Thermaikos Gulf. Concerning the application of the above model to the larger area of the extended Thermaikos Gulf, the model runs showed that, after an algal bloom, the stronger the winds are, the lower the concentration of the dispersed phytoplank-tonic cells is expected to be. This means that stronger winds more effectively disperse the masses, "diluting" the large hazardous populations. Moreover, when the source is located in the inner part of the gulf, the dispersion of the bloom does not seem to reach to the outer parts of the basin under north or south winds. From the applications of the dynamic hydrobiological models described here, it can be said that the spreading of phytoplankton cells in the large area of a gulf or a coastal basin can be investigated analytically, so that the essential measures of prevention can be applied for each different case.
Finally, if time-dependent uncertainty and randomness affect the processes under investigation, and seasonal cycles are displayed, then non-homogeneous Markovian models with cyclic behaviour may be useful tools for their description. Furthermore, these type of models may be used as an alternative approach to the study of diversity and stability of such marine ecosystems.