Surface heat fluxes and ecosystem function in the Cretan Sea (eastern Mediterranean): a modelling study

Abstract. As a component of the Mediterranean Forecast System Pilot Project, a data buoy was deployed in the Cretan Sea. A 1-D ecosystem model of the site has been used to investigate the role of surface heat fluxes in determining modelled ecosystem behaviour. The method of calculation of these fluxes, the quality of the data used, and the temporal resolution of the data all had an impact upon the modelled ecosystem function. The effects of the changes in heat flux formulation were substantial, with both annually averaged properties of the system and the seasonal evolution of the biology being affected. It was also found that the ecosystem model was extremely sensitive to the accuracy of the meteorological forcing data used, with substantial changes in biology found when offsets in the forcing data were imposed. The frequency of forcing data was relatively unimportant in determining the biological function, although lower frequency forcing damped high frequency variability in the biology. During periods of mixing the biology showed an amplified response to changes in physical dynamics, but during periods of stratification the variations in the physics were found to be less important. Zooplankton showed more sensitivity to physical variability than either phytoplankton or bacteria. The consequences for ecosystem modelling are discussed. Key words. Oceanography: physical (air-sea interactions; turbulence, diffusion, and mixing processes) – Oceanography: biological and chemical (plankton)


Introduction
The effective modelling of ecosystems requires a suitable knowledge of both the governing biogeochemical equations and physical processes. The ecosystem function is influenced by physical processes, through changes in the temperature, light and mixing regimes (Huisman et al., 1999;Margalef, 1997;Sharples and Tett, 1994;Pingree et al., 1978). There-Correspondence to: J. I. Allen (jia@pml.ac.uk) fore, the physical parameterisation of the model may have an important effect upon the biological function (Chen and Annan, 2000). Additionally, the choice of forcing at the airsea interface can strongly influence modelled ecosystem behaviour (Lacroix and Nival, 1998).
This work has been undertaken as part of the Mediterranean Forecasting System Pilot Project (MFSPP), which aims to predict the marine ecosystem variability in coastal areas of the Mediterranean Sea. A forecasting system requires two parts, an observing system and a numerical modelling component. The M3A buoy, in the Cretan Sea, has been installed to supply the observational data. The forecast capability of the modelling system is dependent upon the ecosystem model's responses to variability in physical forcing, and temporal and spatial resolution of forcing functions.
The aim of this work is to investigate the role of surface heat flux, as determined by the frequency of meteorological data and the choice of heat flux formulations, in determining the biological function of a one-dimensional ecosystem model of the Cretan Sea. The primary productivity of a system is determined by both nutrient availability and the residence time of plankton in the euphotic zone, which, in turn, are both dependent upon the stability of the water column (Huisman et al., 1999). In a 1-D model this is dependent upon surface fluxes of heat and momentum.
The Cretan Sea, in the eastern basin of the Mediterranean, is a seasonally stratified oligotrophic system, and in the winter months the water column frequently overturns, mixing up bottom waters . Phosphorus is generally considered to be the limiting nutrient in this region. These overturning events are important to the ecosystem function of the Cretan Sea, since they sporadically mix up nutrients to the surface waters. Air-sea transfers of momentum and heat, which determine the timing and extent of these overturning events, are, therefore, an important factor in controlling the local biogeochemistry. The extent to which modelled ecosystem behaviour is affected by the method of modelling surface fluxes of heat and momentum is, therefore, a pertinent question. Castellari et al. (1998)  ing an OGCM of the Mediterranean Sea, that the hydrodynamics are sensitive to the heat flux formulation and meteorological data used. Lacroix and Nival (1998) have also shown that in the western Mediterranean the frequency of forcing has a strong effect upon the biological function.

ERSEM
The European Seas Ecosystem Model (ERSEM) is a generic ecosystem model, with a proven record of use in the Mediterranean Sea (Allen et al., 2002a;Zavatarelli et al., 2000;Allen et al., 1998;Vichi et al., 1998). ERSEM is a modelling framework in which an ecosystem is represented as a network of physical, chemical and biological processes that display coherent system behaviour . A "functional group" approach is used to describe the biota. The ecosystem is subdivided into three functional types: producers (phytoplankton), decomposers (bacteria) and consumers (zooplankton), and subdivided on the basis of trophic links and/or size (Table 1 and Fig. 1). Physiological (ingestion, respiration, excretion and egestion) and population (growth, migration and mortality) processes are included in the descriptions of functional group dynamics. Physiological processes and population dynamics are described by fluxes of carbon or nutrients between functional groups. Each functional group, therefore, has a number of components, each of which is explicitly modelled. These include carbon, nitrogen, and phosphorus for all functional groups and, in the case of diatoms (P1), silicon. Detailed descriptions of ERSEM and its sub-models can be found in , Baretta-Bekker et al. (1995;1998) and Ebenhöh et al. (1997). Sensitivity analyses of parameterisations of ERSEM have been undertaken by Ebenhöh et al. (1997) and Varela et al. (1995). The model used in this study is a version of ERSEM as described above, with the main adaptation being the inclusion of dynamically varying carbon to chlorophyll ratios in the primary producers, following methods described in Geider et al. (1996) (Allen, 2002b).

POM
The ERSEM code is coupled with a 1-D version of the Princetown Ocean Model (POM) (Blumberg and Mellor, 1987). The POM code calculates Richardson number dependent eddy diffusion coefficients for momentum (K M ) and scalar variables (K H ), using the Mellor-Yamada 2.5 turbulence closure model (Mellor and Yamada, 1982), modified after Galperin et al. (1988) and using a prescribed turbulence length scale following Bakhmetev (1932). These coefficients are used to transport variables in both the physical and ecosystem sub-models. The physics is driven at the air-sea interface by surface fluxes of heat and momentum, calculated from meteorological boundary data, with salinity held at the surface to climatological values . Horizontal velocity is derived from surface and bottom stresses, and transported in the same way as other parameters. These velocity gradients are used to determine the shear production.
The vertical resolution of the model is a metre at the surface, increasing a metre at a time to five metres. The resolution stays at five metres for the rest of the 250 m deep water column.

Surface heat flux formulations
The net flux of heat (Q T ) across the air-sea interface is given by: where Q S is the solar radiation flux, Q E is the latent heat of evaporation flux, Q H is the sensible heat flux and Q B is the long-wave radiation flux. Q S is calculated at every time step Table 2. Long-wave radiation formulations; (a) Rosati and Miyakoda (1998) and (b) Budyko (1974) Formulation using astronomical calculations of solar radiation, modified by cloud cover (Dobson and Smith, 1988). The sensible (Q H ) and latent heat (Q E ) fluxes are calculated using standard formulae; the net flux of heat (Q T ) across the air-sea interface is given by:

MEDITERRANEAN SEA
where p A is the air density, q A is the specific humidity of air, q S is the saturation specific humidity of air, C p is the specific heat capacity of water, L E is the latent heat of vaporisation, T S and T A are the surface and air temperatures, respectively, and C E and C H are the exchange coefficients for latent and sensible heat, respectively. The calculations for the latent and sensible heat are dependent upon the calculations of the coefficients C E and C H . Two formulations taken from Castellari et al. (1998) are used; the "neutral" formulation (Rosati and Miyakoda, 1988), which sets both coefficients equal to 1.1 × 10 −3 , and the Kondo formulation (Kondo, 1975), where the coefficients are a function of air-sea temperature difference, wind speed and a stability criterion for the surface waters.
Two formulations for the long-wave radiation flux (Q B ) are used ( Table 2). The formulations for the exchange coefficients and longwave radiation are combined to give three heat flux models (Table 3). These were found by Castellari et al. (1998) to model the hydrodynamics of the Mediterranean most successfully, although it should be noted that the model used in the Castellari et al. (1998) work is different from that used here and, therefore, direct comparisons cannot be made.
European Centre for Medium-range Weather Forecasts (ECMWF) 6-hourly data, at a 2.5 • horizontal resolution, was extracted for a position close to the M3A buoy site to provide the meteorological forcing for the model. Air temperature at 2 m (T A ), 10 m winds, and cloud cover were available. Dew point temperature (T D ) and mean sea level pressure (P A ) data were also available and were used in the calculation of the relative humidity (RH ) based on a formulation in Wallace and Hobbs (1977), and a calculation for the saturation vapour pressure (E S ) (Tetens, 1930):  Rosati and Miyakoda (1998), (b) Kondo (1975) and (c) Temperatures are in degrees Celsius and pressures in Pascals.
The buoy was deployed to collect three-hourly data for a number of physical and biological variables, including temperature and chlorophyll, at various depths. In situ meteorological information was also collected. Further details of the instrumentation used and calibration of data sets can be found in Nittis et al. (2003). Other data sets are also available for the Cretan Sea, and some relevant data are summarised in Tables 4 and 5. Table 5. Nutrient data in the Cretan Sea, and annually averaged nutrient concentrations (depths 0-250 m) for the model simulations PA1, PA2 and PA3: (a) Tselepides et al. (2000); for stations in the Cretan Sea in depths of more than 500 m; annual averages of samples taken from depths 0-200 m (b) Gotsis-Skretas et al. (1999)

Sensitivity to variations in heat flux formulation
The model was run using the three heat flux formulations, labelled PA1, PA2 and PA3 (Table 3), forced using six-hourly ECMWF meteorological data for the year 2000. A repeatingyear forcing was applied for five years to spin the models up to a quasi-steady state for each formulation. The simulations presented here were initialised for both biological and physical variables using the results from these spin-up runs.

Temperature
A comparison of these simulations with the M3A data shows that the surface temperature was underestimated by up to 3 • by all three formulations. The PA1 and PA2 runs were most seriously affected, with PA3, although still significantly underestimating the surface temperature, showing the closest fit to the data (Fig. 3). Simulated temperatures at lower depths (not shown) pick up the correct seasonal trends, but not the high-frequency variability shown by the M3A data.
The inability of the model to reproduce the surface temperature could be related to horizontal advection of heat to the site. It is well documented that the Mediterranean has a negative heat budget (Castellari et al., 1998), compensated for by transport of heat across the Gibralter Sill, and, therefore, a 1-D model with no compensation for this would be expected to underestimate the temperature. An underestimation of the solar inputs to the system may also have an effect, although this has been discounted, since the modelled solar heat flux was substantially higher than literature estimates for the Mediterranean; Garrett et al. (1993) gives the long-term mean solar heat flux as 202 W m −2 , compared to 219 W m −2 for this model.
Other potential causes of this temperature underestimation were also investigated. The ECMWF forcing data used was checked against the meteorological data taken from the M3A buoy (Nittis et al., 2003), and the two data sets showed good agreement for most variables. However, the humidity, as calculated from the ECMWF data, showed substantial differences to the measured M3A data. The importance of the accuracy of the forcing data is investigated further in Sect. 3.3.
The three runs, as mentioned above, showed substantial differences in their calculations of the hydrodynamic properties of the M3A site. As an indicator of the physical behaviour of each model run, strength of stratification, as given by the maximum Brunt-Vaisala frequency and temperature difference across the thermocline, was calculated (Fig. 4).
The onset of stratification occurred around day 80 for all three runs, with the strongest stratification not occurring until after day 100. Differences between the runs were evident, with PA1 giving a deeper thermocline in the initial stages and stronger stratification in the latter parts of the year. PA1 also appeared to give a more diffuse thermocline. The PA3 simulation resulted in the shallowest thermocline depth, and subsequently the extent of the stratification was the greatest.

Comparison of biological time-series data
The physical differences between the modelled water columns have an effect upon the biological function of the model, through the timing of stratification, the water temperature and the extent of mixing. The timing and magnitudes of the spring blooms were quite different for the three simulations. The PA1 simulation bloomed the latest and with the smallest peak (87 days and 0.51 mg-Chl m −3 ). PA2 gave a rather diffuse bloom, with substantial secondary blooms at most depths; these peaks occurred earlier than for the PA1 run (at days 62 and 71), and with a similar magnitude (0.54 and 0.53 mg-Chl m −3 ). The PA3 simulation resulted in an early bloom (day 61) which was of a greater magnitude (0.68 mg-Chl m −3 ) than either the PA1 or PA2 blooms. The comparisons of simulated chlorophyll concentrations with data show that they are of the right order of magnitude (Fig. 5). However, the relatively large peaks in chlorophyll in the model runs are not consistent with the apparent lack of any significant spring bloom in the measured data. There is unfortunately little data available for the period up until day 65, and none between day 40 and 65, so it is possible that some elevated chlorophyll concentrations were present in this period.
The three runs showed remarkable similarity in the summer and autumn, and substantial differences in the winter and spring (Fig. 6). The winter and spring periods coincide with overturning events, and hence, the nutrient supply to, and residence times of, phytoplankton in the euphotic zone are controlled by the physical properties of the water column. In the summer and autumn stratification disconnects the surface waters from the nutrient rich bottom waters. Hence, the limit to primary production is the in situ recycling of nutrients in the surface layers, and the ecosystem comes under biological control. Turley et al. (2000) found a highly significant relationship between bacterial and primary production in the Cretan Sea, and this relationship is well replicated in the simulations (Fig. 7). No obvious differences were evident between model runs, with all showing relative bacterial and primary production of the same order as the Turley data.

Annually averaged biology
The comparison of the three runs with data for phytoplankton production show good agreement (Table 4 and Fig. 8), although the scatter in the data does not allow a judgement to be made on the relative merits of them. The phytoplankton biomass also agrees well with published data for the region, although possibly slightly on the high side in all cases. Both bacterial biomass and production lie within the range of published data. Similarly, modelled nutrient levels show good agreement with observations ( Table 5).
The annual average production and biomass of the three runs showed some significant differences (Fig. 8). The main differences between the three runs lie in the production data, with the zooplankton being particularly affected (the difference between the PA1 and PA2 annually averaged production is greater than 300%). The biomass varied by less than 10% for both the phytoplankton and bacteria, although the zooplankton showed a much more marked variation (a 60% difference between the PA1 and PA2 runs).

Sensitivity to the frequency of meteorological forcing
To investigate the role of surface forcing frequency, the three simulations were each rerun with the meteorological data read in at 12-hourly, daily, weekly and monthly intervals. Annually-averaged simulated zooplankton were more sensitive than either the bacteria or the phytoplankton to changes in forcing frequency (Fig. 9); the bacteria and phytoplankton biomass values were surprisingly unaffected by the use of extremely low temporal resolution surface forcing data. The PA3 run was particularly insensitive to changes in frequency of forcing, with little change in any of the annually averaged production values for the 6-, 12-or 24-hourly forced simulations, and relatively small changes when using weekly or monthly forcing. The production was more sensitive to frequency of forcing than the biomass, with the PA1 run showing by far the most sensitivity.
The chlorophyll concentrations showed little dependence upon the frequency of forcing. The differences were mainly in the fine scale detail, with the simulations that used lower frequency forcing showing damping of the short-term fluctuations in chlorophyll concentrations, but the main chlorophyll peak remaining of largely the same amplitude and timing.

Sensitivity to variations in the forcing data
Inaccuracy in the forcing data could well be a cause of the consistent shortfall in sea surface temperatures (Fig. 3). Comparison of measured (M3A) data and the ECMWF forcing data shows that the ECMWF gives good estimates of air temperature, wind speed and mean sea level pressure, but a poor estimate of humidity. The average relative humidity for days one to 160 was 51% for the ECMWF and 65% for the M3A measured data. To investigate the significance of this difference, a simulation (spun up for a year from the PA3 standard initialisation), forced with the M3A humidity, rather than ECMWF data, was run. The accuracy of the M3A humidity sensor is ±3% (Nittis, pers. comm.). In the second half of the year, and at other times where no M3A data was available, an average value of 65% was used. This run resulted in a marked increase in modelled temperatures, and a much better fit with the chlorophyll data (Fig. 10). It also showed markedly stronger stratification than the previous run using the PA3 formulation (Fig. 4). The chlorophyll peak that was being produced at between 40 and 70 days no longer appeared and a far lower, broader peak at around day 90 was found. It is apparent from this that the ecosystem function of the model is extremely sensitive to the quality of the humidity data used. Even when using rather patchy data, with average values used in periods of no data, a far better result is obtained than if using high resolution, but inaccurate, data.
The quality of the humidity data influences the behaviour of the model, which indicates that the accuracy of other meteorological forcing data should be considered. Therefore, perturbations of the forcing variables were performed to compare the influence that these data have on the surface heat flux. These were run from the standard PA3 initialisation, and the differences in model behaviour over 160 days (the period of M3A data availability) analysed. An increase in air temperature led to increases in sea surface temperature and stronger, earlier stratification, and vice versa. Similarly, increasing the humidity also acted to raise sea surface temperature (and give more pronounced stratification). The wind also had an effect upon stratification, with increased winds giving less stratification as would be expected.
The simulations that gave earlier stratification also had earlier phytoplankton blooms. These earlier blooms were more intense but shorter lived. When the onset of stratification was delayed the bloom became delayed and far less intense. The overall net effect upon the biomass and production is a balance between the longevity and intensity of the bloom, and hence there was no strict pattern as to how a change in air temperature or wind speed affected the net production or biomass (see Table 6). However, it is apparent that even relatively small offsets in the forcing data may have major effects upon the ecosystem function of the model. A comparison between a standard run and runs using perturbed forcing data show that substantial changes in heat flux can be induced (Fig. 11). Large, but realistic, perturbations over a period of a month give changes of the order of ± 50-150 Wm −2 differences in the total heat flux (compared with a heat flux of the order of −200-200 Wm −2 ). The strength of influence of the perturbation does not appear to be dependent upon its sign, Table 6. Percent change in 160 day mean biomass ( B) and production ( P ) for depth-integrated phytoplankton, zooplankton and bacteria with changes in forcing variables. Note the increase in relative humidity is scalar -a 20% increase in a humidity value of 50% gives 70%, not 60%

Phytoplankton
Zooplankton Bacteria and is significant for all three parameters investigated.

Discussion and conclusion
The behaviour of the biological model was shown to be sensitive to variations in the heat flux formulations used. It is difficult to judge the relative merits of these formulations, since they all produced broadly acceptable ecosystem responses. However, the combination of the Kondo and May formulations (simulation PA3) gave the best temperature validation, which is in agreement with Castellari et al. (1998). Annual average biomass and production estimates from the three simulations showed distinct differences, although all showed agreement with literature data. Changes in production were far greater than the changes in biomass, and zooplankton seemed to be particularly sensitive. The timing and amplitude of the spring bloom showed quite substantial differences for the different simulations. All three simulations had peaks of chlorophyll in the period between 60 and 90 days, yet stratification did not occur in any of the simulations until about day 80, with significant levels of stratification not occurring before day 100. In addition, the stratification occurred at similar times in the three runs, although with some difference in intensity, and yet the phytoplankton blooms occurred at quite different times. This does not fit with the classic theory of phytoplankton bloom development, which requires stratification to take place before blooming occurs (e.g. Mann and Lazier, 1996). Huisman et al. (1999) developed a model where blooming may occur without stratification if the turbulence is below a critical level and the light penetration is sufficient for growth to occur. They suggest that in clear waters the critical turbulence theory would become important. The results of these simulations support this; the clear waters of the oligotrophic Cretan Sea allow the turbulence to dictate the onset of phytoplankton blooming, and hence, the timing of the spring bloom is sensitive to changes in turbulence induced by small changes in surface heat flux. This is illustrated in Fig. 12, which shows the temporal evolution of scalar eddy diffusion coefficients, primary productivity and phosphate concentration for the first 90 days of the simulation. In addition, this shows that in the PA2 and PA3 simulations strong mixing events (K H > 0.3 m 2 s −1 ) coincide with net respiration. The reduced exposure to light and poor adaptation to ambient light conditions due to the mixing of phytoplankton out of the euphotic zone is likely to be the cause. The same is not true for the PA1 simulation, where the mixing is less intense. Similarly, there appears to be a threshold of phosphate concentrations below which there is little primary production, as can be seen by the late winter lack of primary productivity in the PA1 simulation, in contrast to the PA2 and PA3 simulations. The PA2 and PA3 simulations have higher average vertical mixing constants than PA1, resulting in enhanced mixing of phosphate into the euphotic zone. A simple rule of thumb for this system seems to be that threshold values of phosphate of approximately 0.004 mmol m −3 must be exceeded and eddy diffusion values of approximately 0.3 m 2 s −1 must not be exceeded for primary productivity greater than 100 mg-C m −2 d −1 to occur. The surface heat flux, which, to a large extent, determines the turbulence of the surface waters, can, therefore, be seen to heavily influence biological behaviour through, first, determining the transport of phosphate up from deeper waters and second in determining the residence time of phytoplankton in the surface, euphotic waters.
The influence of the forcing frequency upon the biology was investigated. The main differences with changes in frequency forcing were found in the high-frequency phenomenon, which were reduced when using lower frequency forcing data. The annual mean biological properties showed relatively little dependence upon the frequency of forcing data used, although the level of dependence varied for the different flux formulations. As with other simulations the major changes in the annual mean values were generally found in the production values, and the zooplankton data were most affected. This is consistent with the observation of Colebrook (1985) that the interannual variability in zooplankton in the North Sea is dependent upon phytoplankton abundance early in the year. It is felt that the model behaves similarly; Figure 8 shows the zooplankton biomass and production to be substantially higher in the simulations where the bloom occurs the earliest (PA3) or persists (PA2). Even though changing the heat flux formulation had significant impacts upon the model function, it still did not significantly alter the fact that all three formulations significantly underestimated the surface heating. The model sim-ulates the M3A site in the Cretan Sea, which is close to two gyres, to the east cyclonic and to the west anticyclonic (Georgopoulos et al., 2000). One-dimensional water column models are incapable of simulating the horizontal processes associated with these gyres. However, the sporadic advection to the gyres is thought unlikely to be responsible for the consistent underestimation of temperature. It is more likely that it is due to the well documented continuous input of heat to Mediterranean waters across the Gibralter Sill.
Other possible sources of error in the model were investigated, and the hydrodynamic properties of the simulations were found to be greatly improved when using humidity forcing based upon M3A data, even though the data set was incomplete and averages had to be used throughout much of the year. Improving the physical performance of the model was matched with a better simulation of the chlorophyll; the chlorophyll maximum was delayed significantly, and reduced in amplitude, to give a dramatically improved validation. Accurate meteorological forcing data, therefore, seems to be extremely important in ecosystem modelling.
In summary, the accuracy of meteorological data is of paramount importance in determining ecosystem behaviour in the 1-D ERSEM/POM model of a Cretan Sea site. The frequency of the forcing data is of only relatively minor importance, and has very little effect on any of the seasonal properties of the system. The heat flux formulation has a significant effect upon the biological function of the system. In the winter and early spring, when the system is overturning and hence, mesotrophic, the biomass is sensitive to changes in heat flux. The system is said to be under physical control, with nutrient availability being primarily dependent upon the strength and depth of mixing; small variations in the physical regime are amplified by the modelled biological system. Later in the year, when the water column stratifies and the system becomes oligotrophic, the variations in heat flux have little effect on the biomass estimates, and the system is under biological control; old production dominates, and the system becomes dependent upon bacterially mediated cycling of nutrients.
Therefore, in systems where this physical control is likely, great care must be taken to effectively model the physics if there is to be any chance of effectively predicting biological behaviour. During times of biological control, the physics becomes less important, and emphasis must be placed on paramaterizing the biogeochemical model. This has implications for data assimilation systems; ideally, both biological and physical parameters would be assimilated, but knowledge of the properties of the system would allow for judgement to be made on the relative importance of biological or physical parameters to the ecosystem function of the model.
Acknowledgements. This work was partly funded by the Mediterranean Forecasting System Pilot Project (EU-MAST Project -MA53-CT98-0171) and partly funded by the UK Natural Environment Research Council through the Plymouth Marine Laboratory core strategic research programme. Thanks are extended to the people involved with the M3A buoy, without which this work would not have been possible. Thanks also go to George Petihakis for the  Fig. 12. Comparison of (a) vertical diffusion coefficient, depth averaged over whole water column (250 m), (b) phosphate concentration, depth averaged over the euphotic zone (surface 120 m) and (c) the depth integrated primary production over the whole water column, for the PA1 (large dashes), PA2 (small dashes) and PA3 (solid line) runs. maps of the Mediterranean and Cretan Seas.
Topical Editor N. Pinardi thanks K. Nittis and another referee for their help in evaluating this paper.