The 100^{th}, 75^{th} and 50^{th} box and whisker plots of UNTRIX values from S. Evoikos (reference site) and Kalamitsi before and after the WWTP operation

## 1. 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-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.

## 2. Methods and case studies

### 2.1. Trix — Estimating the trophic status of a coastal ecosystem

#### 2.1.1. 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 (Chl-a), 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-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-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].

An unscaled trophic index, UNTRIX, was introduced by Pettine et al. [10] following the European Water Framework Directive (WFD, 2000/60/EC):

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.

#### 2.1.2. 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.

Data were collected seasonally, from three fixed stations, from March 2001 until January 2003 (before WWTP operation) and from August 2004 until July 2006 (after WWTP operation).

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 (Table 1), we defined the 100^{th}, 75^{th} and 50^{th} plots of the Southern Evoikos and Kalamitsi areas before and after the WWPT operation.

S. Evoikos(reference site) | Kalamitsi(before WWTP operation) | Kalamitsi(after WWTP operation) | |

100^{th} | 3.097 | 4.091 | 4.063 |

75^{th} | 2.118 | 3.142 | 2.994 |

50^{th} | 1.768 | 2.673 | 2.452 |

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 box-plot 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)].

### 2.2. Piecewise linear regression models with breakpoints

#### 2.2.1. 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-24].

#### 2.2.2. 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 (6) assumes the form:

For all the estimated parameters, the t-test was implemented with:

In all cases, H_{0} is rejected as p<0.05. The fitted model explains about 71 % of the observed variability (r^{2}=0.706).

Plotting the predicted and observed values [20], it was found that the predicted values of Chl-a 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} Chl-a, no relationship could be estimated between Chl-a and NO_{3}. However, when 0.9≤breakpoint≤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.

### 2.3. A numerical dynamic hydrobiological model

#### 2.3.1. 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) [29, 30, 6]:

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-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]:

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

a suitable time step is selected

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

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

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:

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.

**The second part of the model refers to the hydrobiological processes**

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.

Function | Definition | Formula | References |

μ | Phytoplankton growth rate (day^{-1}) | μ_{(T)} . L_{lim} . N_{lim} . P_{lim} | [44, 37] |

μ_{(T)} | Temperature-dependent growth rate (°C) | 0.81 e^{0.0631 T} | [39, 40] |

L_{lim} | Light limitation | [41] | |

N_{lim} | Nitrogen limitation | (Nt: μmole N) | Michaelis–Menten enzyme kinetics |

P_{lim} | Phosphorus limitation | Michaelis–Menten enzyme kinetics | |

(a) | |||

Constants | Definition | In model | In literature |

TL | Total losses (grazing + sinking + natural death) | ~ 0.33 (day-1) estimated from numerical tests and bibliographical references | 0.07-0.67 (season dependent) [45-48] |

I_{opt} | Optimum light irradiance | 110 W. m^{-2} | 100-170 W. m^{-2}[37, 49] |

K_{N} | Half-saturation for nitrogen uptake | 3 μΜ | 0.4-10 μΜ [44, 49-51] |

K_{P} | Half-saturation for phosphorus uptake | 0.2 μΜ | 0.09-3.4 μΜ [52-54] |

(b) |

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]:

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}.The new particles’ position is determined from the positions of other particles randomly selected.

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.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.

**2.3.2. Case studies**

**The hydrobiological model**

The application of the mathematical hydrobiological model described in the previous section is presented here.

**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.

Fig. 3a depicts the study area of the Thermaikos basin, with A being the initial location of the appearance of an HAB episode, B the control position-location with reported measured values of HAB concentrations, T the city of Thessaloniki, P the port of Thessaloniki, E the area called Megalo Emvolo, and M1, M2 and M3 the mussel culture areas lying at the west coasts of the gulf. Fig. 3b depicts the grid described above, as the platform for the numerical solution of the equations of the model.

The present application to the real coastal basin of Thermaikos Gulf tried to simulate the fate of harmful algal cells after a sudden appearance of HAB, using *Dinophysis* abundances 10^{4}-10^{5} cells/lit recorded in north-eastern Thermaikos on January 10^{th}_{,} 2000 (position A in Fig. 3a). For the simulation, a number of 1000 particles was initially used at the location (A), and 619,000 particles were finally counted 21 days later (this latter number included particles that may have dispersed in the coastal basin of the gulf or escaped out of the gulf). The model runs led ultimately to the distribution of *Dinophysis* populations in the coastal basin of the gulf. Two 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), *Dinophysis* 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 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].

**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:

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.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.

Finally, the combination of the following cases was considered:

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,

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,

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.

**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.

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.

### 2.4. A stochastic markovian model

#### 2.4.1. 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-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:

To states | To states at time (t+1) | ||||

From states | time (t) | 1 | 2 | . . | k |

1 | n_{1}(t) | n_{11}(t+1) | n_{12}(t+1) | . . | n_{1k}(t+1) |

2 | n_{2}(t) | n_{21}(t+1) | n_{22}(t+1) | . . | n_{2k}(t+1) |

. . . | . . . | . . . | . . . | . . . . . . | . . . |

k | n_{k}(t) | n_{k1}(t+1) | n_{k2}(t+1) | . . . | 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}, _{0}), spring to spring (P_{1}), summer to summer (P_{2}) and autumn to autumn (P_{3}) are given by:

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 (^{..} =lim(autumn to winter), and give the probabilities of the system in “statistical equilibrium” for each season [65-67].

**2.4.2. 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-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.

From winter to spring | From spring to summer | ||||||||||||||

A_{1} | A_{2} | A_{3} | A_{4} | A_{1} | A_{2} | A_{3} | A_{4} | ||||||||

A_{1} | 0.5 | 0.5 | 0 | 0 | 0.584 | 0.25 | 0.083 | 0.083 | |||||||

A_{2} | 0.714 | 0.143 | 0.143 | 0 | 0.112 | 0.444 | 0.333 | 0.111 | |||||||

A_{3} | 0.461 | 0.461 | 0.078 | 0 | 0 | 0.333 | 0.667 | 0 | |||||||

A_{4} | 0 | 0.5 | 0.25 | 0.25 | 0 | 0 | 0.5 | 0.5 | |||||||

From summer to autumn | From autumn to winter | ||||||||||||||

A_{1} | A_{2} | A_{3} | A_{4} | A_{1} | A_{2} | A_{3} | A_{4} | ||||||||

A_{1} | 0.5 | 0.5 | 0 | 0 | 0.584 | 0.25 | 0.083 | 0.083 | |||||||

A_{2} | 0.714 | 0.143 | 0.143 | 0 | 0.112 | 0.444 | 0.333 | 0.111 | |||||||

A_{3} | 0.461 | 0.461 | 0.078 | 0 | 0 | 0.333 | 0.667 | 0 | |||||||

A_{4} | 0 | 0.5 | 0.25 | 0.25 | 0 | 0 | 0.5 | 0.5 | |||||||

(a) | |||||||||||||||

From winter to spring | From spring to summer | ||||||||||||||

B_{1} | B_{2} | B_{3} | B_{4} | B_{1} | B_{2} | B_{3} | B_{4} | ||||||||

B_{1} | 0.833 | 0.167 | 0 | 0 | 0.786 | 0.071 | 0.143 | 0 | |||||||

B_{2} | 0.712 | 0.224 | 0.064 | 0 | 0.4 | 0.2 | 0.2 | 0.2 | |||||||

B_{3} | 0.222 | 0.222 | 0.334 | 0.222 | 0.2 | 0.4 | 0.2 | 0.2 | |||||||

B_{4} | 0.571 | 0.143 | 0.286 | 0 | 0 | 0 | 0.5 | 0.5 | |||||||

From summer to autumn | From autumn to winter | ||||||||||||||

B_{1} | B_{2} | B_{3} | B_{4} | B_{1} | B_{2} | B_{3} | B_{4} | ||||||||

B_{1} | 0.571 | 0.143 | 0.214 | 0.072 | 0.333 | 0.167 | 0.25 | 0.25 | |||||||

B_{2} | 0.4 | 0.2 | 0.2 | 0.2 | 0.25 | 0.25 | 0.25 | 0.25 | |||||||

B_{3} | 0.25 | 0.25 | 0.25 | 0.25 | 0 | 0.4 | 0.4 | 0.2 | |||||||

B_{4} | 0 | 0 | 0.333 | 0.667 | 0 | 0 | 0.6 | 0.4 | |||||||

(b) |

The matrices P_{i}, i={0, 1, 2, 3} and their limits (_{i} and B_{i}) in statistical equilibrium could be estimated as follows:

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) 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.

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].

## 3. 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 phytoplanktonic 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.

## Acknowledgments

We owe many thanks to the Forestry Research Institute for providing us with wind data for the period of January 2000. We also wish to thank Dr O. Gotsis-Skretas, Dr Ch. Kontoyiannis and Dr A Pavlidou from the Hellenic Center for Marine Research for providing the data from S. Evoikos Gulf.