Assessing the performance of the Cretan Sea ecosystem model with the use of high frequency M3A buoy data set

Abstract. During the Mediterranean Forecasting System Pilot Project a buoy was deployed in the Cretan Sea and for the first time high-frequency physical and biogeochemical data were collected over an extended period, providing a unique opportunity for the evaluation of an ecosystem model. The model both tuned and validated in the Cretan Sea in the past, is explored and quantified. In addition, the optimal parameter set is determined while the effects of high-frequency forcing are explored. The model results are satisfactory, especially at the upper part of the water column, while the inability of 1-D modelling in fully exploring the hydrodynamics of the particular area is depicted and further developments are suggested. Key words. Oceanography; general (numerical modeling) – Oceanography; biological and chemical (ecosystems and ecology)


Introduction
Heavily populated coastal areas bound the Mediterranean Sea.The economies of these regions are highly dependent on fishing, transportation, recreation and other industries, which in turn depend on a healthy coastal environment.Predicting the behaviour of the marine environment is an essential part of the management of marine resources under anthropogenic stress.Therefore, It is an essential requirement for the marine science community to determine the potential time-scales of predictability of the marine ecosystem if an operational coastal ocean environmental monitoring and forecast system is to be developed.Such a system would provide estimates of the changes in both the physical and the biogeochemical marine environments and an enhanced understanding of the marine ecosystem, essential to guiding resource management.Additionally, it would form an early warning system of potentially harmful ecological events and aid the formulation of cost-effective preventive and remedial Correspondence to: G. Triantafyllou (gt@imbc.gr)measures.Towards these issues an existing ecosystem model both tuned and validated in the Cretan Sea (Triantafyllou et al., 2002), is explored and its performance is quantified.As a first step historical and recent data sets are used to exploit the model to determine the optimal parameter set.The particular ecosystem model was chosen due to its generic nature in response to the physicochemical environment within which it is placed.With all significant biological pathways in the modelled system included, the model can respond to the physical and chemical forcing in a way that is at least qualitatively correct under a wide range of conditions.
In this study, high-frequency forcing data from the POSEI-DON system (Soukissian et al., 1999) and the M3A buoy are used and compared with climatological forcing, in an attempt to investigate whether the performance of the ecosystem model is enhanced.High-frequency forcing model results are further validated with in situ data of nutrients and chlorophyll collected during maintenance trips.Finally, high-frequency data of chlorophyll, nitrate, oxygen and temperature during M3A deployment is hindcast with model results.

Study area
The area under study is the M3A mooring site located at 35 • 07 07.54 N, 24 • 59 44.66 E in the Cretan Sea, north of Heraklion, with a depth of 1100 m (Fig. 1).The open character of the area in conjunction with the significant distance from the north coast of Crete (21 miles) ensures that there is no influence from land activities in the functioning of the specific ecosystem.The simulation period is the M3A deployment period starting from 30 January 2000 until 20 April 2001.
The Cretan Sea at the eastern part of the Mediterranean is dominated by multiple scale circulation patterns.The hydrological structure is complex and is characterised by mesoscale variability.It is an area of deep-water formation with a sporadic and episodic formation of intermediate water occurring predominantly in the late winter.From spring through to late autumn/ early winter the region is thermally stratified to a depth of 50-70 m.The ecosystem of the outer Cretan Sea has two different modes of operation.During periods of stratification it exhibits a microbial food web, when small phytoplankters dominate, taking advantage of their superior surface-to-volume ratio and a considerable portion of the total organic pool is in the form of detritus.Phosphorus is thought to be the limiting nutrient for phytoplankton and bacterial growth (Becacos-Kontos, 1977;Berland et al., 1980;Krom et al., 1991Krom et al., , 1992;;Thingstad and Rassoulzadegan, 1995) with very low concentrations even below the euphotic zone at 200 m.Nutrients are recycled in the top layers with a prominent P limitation of both phytoplankton and bacteria, and high surface concentrations and vertically decreasing gradients of DOP and DOC.During this mode very little energy is passed to the higher trophic layers and to the benthos.When the thermocline erodes the system switches into the classical mode during which trapped nutrients find their way to the euphotic zone, initiating a new production bloom and biomass distribution, being linearly ordered according to trophic levels along which the organisms increase sequentially in size (Smetacek and Pollehne, 1986).
The variability observed in biogeochemical parameters is due to the complex hydrological structure characterised by eddies interconnected with jets without defined time scales (Balopoulos, 1996).These features are important because they can transport, entrap and disperse chemicals, particulate matter, small organisms, heat, etc., while they significantly modify the vertical mixing patterns in an area, creating among other phenomena zones of upwelling and downwelling that in turn can have a major effect on the biological productivity (Krom et al., 1993).Those structures, in conjunction with the presence of transient Mediterranean water (TMW), a mass rich in nutrients and poor in oxygen, can drastically alter the chemical conditions of the intermediate layers of the entire south Aegean (Theocharis et al., 1999), acting as a nutrient source to the euphotic zone or as sinks for organic material to the deeper layers.

M3A data set
One of the innovations of the MFSPP project and the deployment of the M3A buoy was the collection for the first time of high-frequency physical and biogeochemical data over an extended period in the outer Cretan Sea.Previous studies in this region were conducted on a seasonal time-scale, which, in the case of oligotrophic systems, may prove to be inadequate in revealing the ecosystem dynamics since the temporal scales must be considered according to the biological time scales of growth (or generation time), which are rather short (Abbott, 1993).The M3A buoy provides 3-h measurements of temperature, chlorophyll-a, nitrate and oxygen, at various depths giving a unique opportunity for hindcasting modelling of the Cretan Sea ecosystem.The locations of the CTD and CT sensors were at 40, 60, and 90 m at line 2 and at 150, 250, 350 and 500 m at line 1.The chlorophyll-a and oxygen sensors were deployed at 40, 65, 90 and 115 m and the nitrate sensor at 45 m.A detailed description of the structure and setup of the three moorings comprising the M3A array are presented in Nittis et al. (2003).Chlorophyll-a, dissolved oxygen and nitrate data were post-calibrated using in situ bottle measurements.The calibration coefficients for the M3A sensors were calculated separately for each period between maintenance cruises, using the reference data collected during each re-deployment when all sensors had been cleaned (Nittis et al., 2003).

Ecosystem model
In this study the model used is a highly portable coupled physical biogeochemical water column model (Allen et al., 1998) adapted for the oligotrophic Cretan Sea (Triantafyllou et al., 2002).The biogeochemistry is based on the European Regional Seas Ecosystem Model (ERSEM) (Baretta et al., 1995), while the physics are provided by a one-dimensional version of the Princeton Ocean Model (Blumberg and Mellor, 1978).
A turbulence closure model (Mellor and Yamada, 1982) determines the vertical temperature, turbulent kinetic energy, and diffusion coefficient profiles generated by a surface heat flux, salinity and wind stress.Heat transfer is assumed to be via vertical diffusion processes.A background viscosity parameterises other mixing processes (e.g.internal waves).The physical sub-model, in addition, to transport describes irradiance, temperature, and sea surface boundary conditions.The solar radiation is provided to the biogeochemical model as photosynthetically available radiation (PAR) penetrating the water column and becoming extinct with increas- ing depth, according to an extinction coefficient calculated at each time step from the concentrations of primary producers, detritus, bacteria and inorganic suspended matter (Zavatarelli et al., 2000).
The basis of the biogeochemical model consists of state variables with a functional group, approach where organisms are grouped according to their trophic level subdivided according to size classes or feeding method.Two sub-models describe the ecosystem, one for the pelagic and another one for the benthic ecosystem, coupled in terms of carbon and nutrients.The pelagic sub-model includes all the important components, such as nutrient pools, oxygen, phytoplankton, mesozooplankton, microzooplankton, bacteria and detritus.Although the full benthic sub-model depicts a rather diverse system, in this study, due to the depth (> 1000 m) and the oligotrophic nature of the system, benthic processes are described by the simple first order benthic returns module.According to this the detritus falling into the benthos is added in both the particulate and dissolved benthic pools from which a proportion of the total carbon (5% of the particulate and 10% of the dissolved fraction) is remineralised and returned to the water column as CO 2 .The nutrient components of these forms of detritus follow the same proportions with 10% of the remineralised nitrogen added to the nitrate and the rest 90% to the ammonium pools.The remineralised phosphorus is fluxed to the inorganic phosphate pool.
The model food web is modified from the standard ERSEM version 11 (Baretta et al., 1995), according to the characteristics of the ecosystem under study as they have been revealed from literature and analysis of in situ data (Azov, 1991;Stergiou et al., 1997;Tselepides and Polychronaki, 1996) (Fig. 2).Thus, phytoplankton is com-posed of four groups with diatoms and flagellates being grazed by microzooplankton and omnivorous mesozooplankton, nanoalgae grazed by microzooplankton, while picoalgae and bacteria are grazed by heterotrophic nanoflagellates.To take into account the importance of bacterial nutrient limitation and the role of dissolved organic material in this oligotrophic system, a modified version (Allen et al., 2002) of the standard ERSEM pelagic sub-model bacteria (Baretta-Bekker et al., 1998;Baretta-Bekker et al., 1995) was used.
The 1-D Cretan Sea ecosystem model (Triantafyllou et al., 2002) which was used has a vertical resolution of 40 layers and total depth of 1100 m, with a finer resolution at the euphotic zone (0-150 m every 10 m).It was forced with six hourly wind speed, humidity, cloud cover and air temperature data obtained from ECMWF.Incident sea surface radiation is calculated from the latitude modified by the cloud cover data using the methods of Patsch (1994), while the extinction coefficient was parameterised according to the measured depth of 1% incident solar irradiance (80-90 m) (Gotsis-Skretas et al., 1999;Ignatiades et al., 1995;Psarra et al., 2000).In this study the model was also forced with high-frequency wind speed and air temperature data for the computation of wind stresses (Hellermann and Rosenstein, 1983) and relaxed to mean monthly sea surface temperature and salinity obtained from a POSEIDON buoy (Soukissian et al., 1999) located 10 miles south of M3A during the period of the M3A deployment (January 2000-April 2001).

Optimal parameter set
The optimal parameter set for the autotrophic and heterotrophic groups in the model are given in Tables 1 and 2, respectively.The maximum specific uptake rate in all phytoplankton groups was lowered to avoid excess production of dissolved organic carbon (DOC) and subsequent growth of bacteria, which will compete for nutrients with phytoplankton.Faster rates of small cells are, however, retained.Oligotrophic phytoplankton was shown to grow at rates of 0.5-2 d −1 , significantly faster than bacterioplankton (Laws et al., 1987), as expected if one considers that in such systems the heterotrophic biomass exceeds that of autotrophs.In oligotrophic systems phytoplankters exhibit a dormant behaviour where photosynthesis rates are minimum, relying on a chance encounter with nutrient patches to revive their photosynthetic mechanism (Smetacek and Pollehne, 1986).In addition the maximum N:C and P:C ratios for phytoplankton have been altered to restrict the luxury uptake of nutrients.
In the microzooplankton module some food preference factors have been changed to account for the fact that dinoflagellates, unlike standard ERSEM V.11 where they were considered inedible, are parameterised as a food source for omnivorous microzooplankton.The bacterial specific uptake rate was significantly lowered, in order to reduce the bacterial growth rate in the range of 0.006-0.086(d −1 ), as suggested by Turley et al. (2000) and Pedros-Alio et al. (1999).Finally the assimilation efficiency was also lowered, simulating an average bacterial growth efficiency of 5-20% (Carlson and , 1996;Kirchman et al., 1991;Pedros-Alio et al., 1999;Turley et al., 2000).(Nittis et al., 2003).

Results and discussion
At 65 m the simulated temperature follows the general trend in the data, but fails to reproduce the observed variability.This variability is less pronounced than at 40 m and can be ascribed to similar processes.At both 90 m and 150 m there is no significant seasonality in the data and the model reproduces the observations.

In situ validation
The in situ data collected during the maintenance trips was compared with model results for nutrients and chlorophyll.
Figure 4   Field data (phosphate, nitrate, ammonia and silicate) indicate the presence of a weak nutracline at 80 m, which is successfully reproduced, although more pronounced in the model results.In the top 100 m the model phosphate values are close to the minimum measured concentrations, with the data fit improving below 100m.It should be noted that compared to the CINCS data set for the same region (Tselepides et al., 2000) (April 1994-July 1996), phosphate values during the M3A deployment period are significantly higher.Nitrate exhibits a similar picture to phosphate with model concentrations close to the lower end of the range of in situ data.Once again, M3A nitrate concentrations are higher than the CINCS data and exhibit significantly higher variability at all depths.Simulated silicate and ammonia concentrations lie within the variability of measured concentrations.The maximum measured chl-a values at 80m are indicative of the deep chlorophyll maximum (DCM), a characteristic of the area (Tselepides et al., 2000), which is reproduced by the model at a slightly shallower depth.This is due to the fact that in the model the DCM coexists with the maximum biomass, while the analysis of the in situ data shows that the DCM is deeper compared to the biomass or cell count maximum (Psarra et al., 2000), indicating an elevated chlorophyll per cell content at deeper layers.The model fails to reproduce chlorophyll in the deeper parts of the water column below 10 m, where small chlorophyll levels are detected even down to 200-400 m due to wind mixing, as has been suggested by Tselepides and Polychronaki (1996).

Chlorophyll
Figure 5 shows the simulated chl-a against the in situ fluorometer derived concentrations.Model results are very close to the measured concentrations at 40 m, with the exception of the period from 22 August 2000 until 29 October 2000, which can be attributed to biofouling of the sensor and the vertical displacement of the mooring line (Nittis et al., 2003).During this period the comparison of model and data is meaningless.This argument is strengthened by the fact that after 31 October 2000, the measured chl-a concentrations are once more very close to the model simulation.
At 65 m the model overestimated chlorophyll during March and April 2000, producing a peak at the end of March due to intense mixing.For the rest of the deployment period the simulation is very good, closely following the measured At 90 m field data was unavailable for the first seven months due to sensor problems.For the rest of the deployment period the model, although following the measured trend, systematically produced lower concentrations.Analogous model behaviour is observed at 115 m, where there is an underestimation of Chl-a, although differences appear significant only during the first seven months.This poor behaviour of the model may partly be due to the use of a constant C:Chl ratio at all depths.It is well documented that there is a common trend of an increase in chlorophyll-a and other light harvesting pigments as growth irradiance decreases (Dubinsky, 1992;Estrada, 1985;Estrada et al., 1993;Gasol et al., 1997;Venrick, 1990).However, there is a rather wide range of C:Chl values suggested by researchers for different environ- ments, which confuses the parameterisation issue.A variable C/Chl ratio of decreasing values with increasing depth would improve the simulations at the deeper layers.An additional further source of discrepancy may arise from the calibration of the data, since there is a wide variation even within three days in the in situ measurements.

Oxygen
The dissolved oxygen sensors only gave reliable data during the first six months of the buoy deployment and then failed at the beginning of July 2000 (Fig. 6).At 40 m model results generally underestimate the oxygen concentration, except at the end of March 2000.Since the region is one of low productivity, this discrepancy can probably be ascribed to the fact that the model overestimates temperature at that depth.A  better agreement between simulation and field data is found at 65 m and 115 m, where the temperature simulations have a better fit with data.At 65 m the maximum oxygen concentrations are found during April-May, which may be attributed to increased primary production, as evidenced by the Chl-a concentrations.

Nitrate
Although model results (Fig. 7) are within the right range of measured values, they do not closely follow the in situ variability.The extreme high values measured at the beginning of the data set are due to an unsuccessful deployment during which the analyser was almost at the surface.The marked decrease in nitrate during mid-May, suggesting an increased activity of primary producers, cannot be explained by the chl-  area.

Sensitivity analysis, the response of the model to changes in forcing
The effects of changes in the atmospheric forcing upon the model behaviour were investigated by making two separate runs.In the first run, climatological wind forcing was used and the model was relaxed to climatological sea surface temperature, and salinity obtained from the Mediterranean Ocean Data Base database (Brasseur et al., 1996).For the second run, high-frequency wind forcing (half hour), for the period from 1 January 2000 until 3 February 2001, for the specific area, was obtained from the POSEIDON system, while the model was relaxed to mean monthly in situ sea surface temperature and salinity.The simulations of chloro-phyll at the depths of 10, 30, 50 and 70 m from the two runs are shown in Fig. 8.In the deeper layers, at 50 and 70 m, there are no significant differences between the two runs, and the model produces a spring and autumn bloom.At 50 m the spring bloom is observed during March, while at 70 m the spring bloom is observed during the end of April/beginning of May.In both layers the autumn bloom is shown during September-October.
However, in the upper layers (10 and 30 m), there are significant differences between the two runs, with more pronounced minimum and maximum chlorophyll concentrations in the run with the high-frequency data forcing, due to the diffusion of nutrients caused by the vertical mixing processes, which are much better described.
Thus, the model results indicate that the 1-D physical model only responds to sea surface forcing in the near surface layers and that the effect of high-frequency forcing data in the simulation of the ecosystem is only significant in the mixed layer.As mentioned above, the appearance of the Transition Mediterranean Water (TMW) at approximately 300-500 m, characterized by high nutrient and low oxygen concentrations in conjunction with the gyral features populating the area, significantly affects the intermediate layers of the Cretan Sea.Thus, in the deeper layers of the euphotic zone, the ecosystem functioning is significantly affected by the advective processes, which cannot be described by the present model requiring a 3-D model.

Model performance
The performance of the model is evaluated using the misfit error between modelled temperature (Fig. 9), Chl-a concentrations (Fig. 10), and in situ data.The model underestimates the temperature at 40 m by around 10-15% for most of the year, except in the autumn, when it is overestimated.At 65 m a similar trend is observed but the errors are smaller (5-10%).The model is consistently colder than the observations, which implies that the buoyancy is weaker than reality and that vertical mixing may be overestimated.The implementation of a surface heat flux formulation in the model may help to overcome these problems, as it would place the ecosystem model in a fully 3-D hydrodynamic environment.The data mismatch of chlorophyll is shown in Fig. 10.At 40 m the model consistently underestimates the data by 50-100%, except during July, while at 65 m Chl-a is consistently overestimated by 50-100%, except during the bloom period, when the error is larger.There are a number of possible reasons for this discrepancy.The attenuation of light down the water column may be incorrect, and as previously, noted the vertical mixing may be wrong.Alternatively, the use of a fixed carbon to chlorophyll ratio in the model may cause this difference, with the ratio being underestimated in the near surface water and overestimated towards the bottom of the euphotic zone.The absence of phytoplankton biomass and primary production measurements precludes us from making further conclusions.However, it seems apparent that further work is required to improve the response of the phytoplankton model to low light environments.One possible solution is to implement a variable carbon-chlorophyll model, such as that of Geider et al. (1997).
The fit with nitrate at 45 m is shown in Fig. 11.In situ data and the model show a good fit in March and June/July, but overestimates the concentrations by 200-300% in between.The excess of nitrate in the spring may arise from the lack of horizontal advective processes in the model or through problems with model phytoplankton uptake of nutrients and nutrient recycling mechanisms.Oxygen data (not shown) underestimates the data by up to 5% at 40 m, 65 m and 115 m.This reinforces the previous conclusion that in a low productivity environment, the oxygen concentration is primarily temperature dependent.
Alternatively, data fits can be substantially improved by the implementation of a data assimilation scheme, such as the Ensemble Kalman Filter.This has been demonstrated for the M3A data set by Allen et al. (2003).

Conclusions
For the first time the oligotrophic ecosystem of the Cretan Sea has been monitored over an extended period of time, providing valuable information and a unique opportunity for ecosystem model application.The Cretan Sea ecosystem model developed during the MATER project when forced with real atmospheric data simulates very satisfactorily the chl-a time series in the upper part of the water column as recorded by the fluorometers.Deeper down, close to the DCM the underestimation of the model concentrations can be attributed to the elevated cell chlorophyll content, a physiological adaptation to lower light conditions.A future modification of the ecosystem model to account for this phenomenon, as well as the implementation of heat fluxes computed for the particular area, are considered necessary.
Although the in situ nitrate data are rather incomplete as a result of instrument malfunction such data, are, however, very important in these studies.Considering the phosphate limitation of the area high-frequency phosphate data from an autonomous instrument would provide significant information towards the improvement of the model.
Laboratory calibration of the fluorometric sensors is considered crucial and necessary, in addition to the calibration with in situ data due to the significant variability of the chlorophyll concentrations in the samples in short-time periods (2-3 days).
Due to its 1-D nature, the model does not include important processes of the particular area, such as horizontal advection and vertical mixing arising from circulation patterns.Nevertheless, it provides a very successful numerical base for further development towards a 3-D forecasting ecosystem model and the implementation of assimilation techniques.
against in situ measurements for the 40, 60, 90 and 150 m depths are shown in Fig.3.Stratification started approximately at the middle of April, reaching the highest subsurface values towards the end of October, and the location of the thermocline lay between 40 and 65 m.From November onwards a gradual mixing of the thermocline into deeper layers occurred until it disappeared in early January(Nittis et al., 2003).The bottom of the seasonal thermocline was at 50-70 m and, thus, the signal of the annual cycle of heating and cooling is hardly detected below that depth.At 40 m the model results follow very closely the field values, with minimum temperatures in March and Maximum in October.The large variability in the in situ values of almost 4 • C within a few days after June 2000 may be associated with changes in the eddy structure of the region.While the variation during October 2000 and December 2000 to January 2001 was coincident, with significant movements of the mooring line shows simulated annual mean values plotted against depth, along with mean in situ concentrations for the period from 30 January 2000 until 19 April 2001, while range bars indicate minimum and maximum measured values.Model results are generally in good agreement with the observed in situ data in the whole water column, although the field data exhibits a substantial variability in almost all parameters.

Fig. 4 .
Fig. 4. Validation of model Phosphate, Nitrate, Ammonia, Silicate and Chl-a with in situ data obtained during maintenance trips.All sampling measurements at a particular depth are treated as replicates and averaged.Range bars are used to indicate variability (minimum and maximum values).

Fig. 8 .
Fig. 8. Model chlorophyll simulations at various depths with climatological and real-time forcing.

Table 2 .
Blackford and Radford (1995)on functional groups and bacteria.Z5 = microzooplankton, Z6 = heterotrophic flagellates and B1 = bacteria.The parameters that have been changed from the standard ERSEM version 11 are indicated in bold.Parameter names follow the nomenclature described byBlackford and Radford (1995)