Adriatic shelf

Abstract. In this paper we show results from numerical simulations carried out with a complex biogeochemical fluxes model coupled with a one-dimensional high-resolution hydrodynamical model and implemented at three different locations of the northern Adriatic shelf. One location is directly affected by the Po River influence, one has more open-sea characteristics and one is located in the Gulf of Trieste with an intermediate behavior; emphasis is put on the comparison with observations and on the functioning of the northern Adriatic ecosystem in the three areas. The work has been performed in a climatological context and has to be considered as preliminary to the development of three-dimensional numerical simulations. Biogeochemical model parameterizations have been ameliorated with a detailed description of bacterial substrate utilization associated with the quality of the dissolved organic matter (DOM), in order to improve the models capability in capturing the observed DOM dynamics in the basin. The coupled model has been calibrated and validated at the three locations by means of climatological data sets. Results show satisfactory model behavior in simulating local seasonal dynamics in the limit of the available boundary conditions and the one-dimensional implementation. Comparisons with available measurements of primary and bacterial production and bacterial abundances have been performed in all locations. Model simulated rates and bacterial dynamics are in the same order of magnitude of observations and show a qualitatively correct time evolution. The importance of temperature as a factor controlling bacteria efficiency is investigated with sensitivity experiments on the model parameterizations. The different model behavior and pelagic ecosystem structure developed by the model at the three locations can be attributed to the local hydrodynamical features and interactions with external inputs of nutrients. The onset of the winter/spring bloom in the climatological simulations is primarily driven by local stratification conditions. During summer the major carbon-transfer pathway developed by the model is the microbial web at all the sites, indicating that a large fraction of organic matter is processed through bacteria during productive periods, as suggested by field observations. The site directly influenced by riverine inputs differs from the others, showing a more alternate shifting among trophic pathways. Applying the conceptual scheme proposed by Legendre and Rassoulzadegan (Ophelia, 41, 153-172, 1995), it can be recognized as a herbivorous spring phase tightly followed by a microbial loop development, a summer microbial web phase and a multivorous trophic web pattern during autumn with a subsequent recovery of microbial processes. Results are discussed in terms of regime shifting from transient to stable water column conditions. Key words. Oceanography: general (continental shelf processes; numerical modelling) – Oceanography: biological and chemical (biochemistry and food chains)


Introduction
This paper deals with one-dimensional simulations of the northern Adriatic Sea ecosystem dynamics.The work has been developed within the framework of the Mediterranean Sea Forecasting System Pilot Project (MFSPP), whose longterm goal is the development of reliable tools to forecast the Mediterranean marine coastal environment at the level of the primary producers.
The preliminary and necessary initial step toward this goal, is the application of state-of-the-art, one-dimensional ecosystem models at selected implementation sites, and the realization of numerical simulations under climatological forcing functions in preparation for three-dimensional numerical simulations.Such an effort allows for a proper model validation, aimed to assess the skill of the model in reproducing observed features of the marine coastal environment.
Among the coastal regions of the Mediterranean Sea, the northern Adriatic basin (Fig. 1) emerges as a very interesting area for carrying out this modeling effort.The coastal ecosystem is characterised by a wide spatial and temporal variability determined by the atmospheric forcing functions, the general circulation, the continental (riverborne) inputs of nutrients and organic matter, and the water column-sediments interactions.Trophic conditions varies from mesotrophic (with occasional dystrophic crisis) in areas directly affected by external inputs, to oligotrophic in the central part of the basin (Zavatarelli et al., 1998;Stachowisch, 1984;Giordani et al., 1992;Fonda Umani, 1996).This variability imposes a strong constraint on the ecological model to be used in this area, as the model has to be generic enough to reproduce this wide trophic range, and the associated biogeochemical interactions, dependent only on the physical and biogeochemical boundary conditions used.
The model has been implemented in three different locations of the northern Adriatic Sea that are representative of the spatial environmental variability of the basin and that are of interest to monitoring programs, so that their distinctive features can be defined and compared with the model results.Two of these locations were already the subject of a modeling exercise in the past (Vichi et al., 1998a, b).We show here new simulations for these two locations, because (as described below) the model has been updated with a more detailed description of the dynamics of the dissolved organic matter (DOM) and the associated bacterial processes, which required a re-assessment of the simulations previously carried out.Moreover, here we go further in the comparison between model results and observations, since we extended it to the model predicted biogeochemical rates (primary and bacterial production), in an effort to provide insights on the model behavior in reproducing correctly the fluxes of matter within the planktonic food web.
The analysis of the ecosystem functioning at the three implementation sites and the role of the different food chains in transferring matter from lower to higher compartments have also been the subject of investigation.Extending the work done in Vichi et al. (1998b) we have analyzed model results from the point of view of the conceptual marine ecosystem functioning scheme proposed by Legendre and Rassoulzadegan (1995) with particular emphasis on the role of DOM and microbial community processes in the dynamics of the northern Adriatic Sea.
In Sect. 2 we provide a general characterization of the implementation sites, while Sect. 3 illustrates the model used with particular emphasis on the hydrodynamical parameterizations and the description of the DOM dynamics.Results are described, discussed and validated with observations in Sect. 4. In Sect. 5 we describe some features of the pelagic food web functioning and carbon cycling at the study sites and in Sect.6 we draw conclusions.

Sites of implementation
Three different areas have been chosen (Fig. 1), each one representing some characteristics of the northern Adriatic shelf and having different interactions with the hydrodynamic features of the basin.Sites S1 and S3 are placed in the northwestern part of the Adriatic Sea, while site AA1 is located in the Gulf of Trieste in the northeastern region.The site names correspond to specific locations that have been of interest to (time-limited) multidisciplinary monitoring programs (except for AA1), but, in the context of their climatological seasonal characterization, we use them with the broader meaning of hydrographic sub-regions affected by different dynamical processes.
The hydrodynamical characteristics at S1 are directly influenced by the Po River runoff (1500 m 3 s −1 on annual average).S1 is located at about 5 km offshore from the Po delta.The bottom depth is 20 m.On the contrary, S3 reflects conditions more typical of an open-sea area, less affected by landward derived inputs.It is located at about 37 km offshore of the western Adriatic coast; the bottom depth is 30 m and is approximately in the central part of the northern Adriatic cyclonic gyre (Artegiani et al., 1997b).Both stations show a similar seasonal variability in the physical parameters, but the water column at S1 is stratified for most of the year due to the freshwater input, while the S3 site is well mixed during the winter and autumn seasons.S1 site is also characterised by the Po River input of large amounts of inorganic nutrients, organic and inorganic suspended matter (Cozzi et al., 1999;Pettine et al., 1998Pettine et al., , 1999;;Giani et al., 2000Giani et al., , 2001)), the latter

Z Z Microzooplankton
(3) Fig. 2. General overview of the biogeochemical state variables and matter fluxes implemented in the pelagic module of ERSEM.Square boxes represent functional groups defined in the model, indicated as vectors in which the subscript i = c, n, p, s represents some of the exchanged element.Continuous arrows indicate fluxes of carbon and inorganic nutrients, dashed arrows the fluxes of inorganic nutrients (N, P, Si) alone and dotted arrows the gas exchange of O 2 , N 2 and CO 2 .Small double arrows above the boxes indicate fluxes with the benthic compartment of the model system.Modified and updated from Baretta et al. (1995).mainly constituted of clay and finer silt.The hydrodynamics of site S3 instead is mainly influenced by the cyclonic general circulation system of the northern Adriatic basin, even if intermittent phytoplankton blooms can be associated with the offshore spreading of the Po River plume.
The Gulf of Trieste, where the AA1 station is located (Fig. 1), is a semi-enclosed basin with a surface area of about 600 km 2 and a maximum depth of 26 m.The whole area is strongly affected by river runoff, expecially along the shallow northwestern coast (Isonzo River), with an average freshwater input ten times lower than the Po River (150 m 3 s −1 ) (Stravisi, 1983a, b;Naudin et al., 1996).However, due to the small volume of the basin (9.5 km 3 ), the average volumespecific discharge is nearly three times higher than in the northern Adriatic as a whole (Malej et al., 1995a).The hydrological features of the basin show a very large interannual and seasonal variability due to the strong easterly winds characteristic of the northern area (Bora), particularly intense in the Gulf of Trieste.The structure of the plankton community reflects the variability of the abiotic factors: species composition and relative abundances change dramatically from year to year (Fonda Umani et al., 1995;Malej and Fonda Umani, 1995;Cataletto et al., 1995;Mozetič et al., 1998).The AA1 site is located in the central part of the basin where the bot-tom depth is 20 m, and it is an operational site of the monitoring programme of the Marine Biology Laboratory (LBM) of Trieste.

General features
The biogeochemical flux model used in this work derives from the European Regional Seas Ecosystem Model (ERSEM, Baretta et al., 1995).ERSEM is a biomass-based biogeochemical flux model which was originally constructed to simulate the dynamical cycling of carbon, oxygen and the macronutrients N, P and Si over the seasonal cycle in temperate marine systems; the model consists of an interlinked set of differential equations, describing the biological and chemical processes both in the water column and in the benthic system, as forced by external environmental conditions (light, temperature, water hydrodynamics and allochthonous nutrient sources).One of the main interesting features that made it suitable for the MFSPP framework is that, having all significant biochemical pathways implemented in the system, the model can respond to the physical and chemical forcing in a way that is at least qualitatively correct under Background turbulent diffusivity for T and tracers (χ 1 , m 2 s −1 ) 2.5 10 −5 1.3 10 −5 1.0 10 −5 Background turbulent diffusivity for S (χ 2 , m 2 s −1 ) 1.3 10 −6 9.0 10 −6 1.3 10 −6 Apportioning coefficient for infrared light (ε i , %) 90 50 80 Attenuation coefficient for visible light (λ v , m −1 ) 0.17 0.17 0.17 the wide range of trophic conditions observed in the Mediterranean Sea.
The trophic web structure implementation follows the functional group approach and the concept of standard organism introduced in Baretta and Ruardij (1988).The pelagic functional groups defined in the present model formulation are shown in Fig. 2, which also gives a schematic view of the functional interactions among groups and other state variables.Since the benthic system dynamics are not the primary subject of this paper -while being a fundamental part of the model system and used in the simulationswe refer to Vichi (2002) and to the original publications of Ebenhöh et al. (1995), Blackford (1997) and Ruardij and van Raaphorst (1995) for a complete description of benthic processes and dynamical variables.Each functional group exchanges matter in units of carbon and macronutrients with the other system components; therefore, all the state variables in the model are mathematically treated as vectors (see Sect. 3.4).A functional group constitutes an implicit sizeclass, i.e. there is no internal size structure in the biological constituents.All the organisms considered to be in a particular group (i.e.diatoms P (1) and picophytoplankton P (3) in phytoplankton) share the same functional properties in the ecosystem and have the same trophic level, but the two groups can have different physiological parameter values and relationships with the other components, such as, for instance, specific basal respiration rates or predator's preference factors.
The biogeochemical flux model is coupled with the one-dimensional version of the Princeton Ocean Model (POM, Blumberg and Mellor, 1987), which is essentially a vertically-resolved boundary layer model based on the second moment turbulence closure scheme derived by Mellor and Yamada (1982) and Galperin et al. (1988).The model determines the dynamical vertical structure and the actual turbulent diffusive vertical transport in the water column, as driven by boundary condition exchange fluxes at the air-sea interface.The suitability of the coupling between ERSEM and POM in the Adriatic Sea has been already tested in other one-dimensional applications (Allen et al., 1998;Vichi et al., 1998a, b) and in a three-dimensional idealized simulation experiment in Zavatarelli et al. (2000).

Vertical discretization
The water column at the implementation sites has been discretized in the vertical by applying a logarithmic distribution of layer depths at the surface and at the bottom, and a uniform depth distribution in the interior.The vertical grid at S1 and AA1 is composed of 30 levels with 6 logarithmic levels both at the surface and at the bottom (the depth of the first layer is 0.1 m, the central layers are 1 m deep).Site S3, which is deeper, has been divided in 40 vertical levels, in order to have the same resolution of 1 m in the central layers (surface layer is 0.05 m deep).

Model equations and boundary conditions
The numerical model solves with a finite difference scheme the one-dimensional differential equations for (U, V , W ) momentum components, pressure P , temperature T , salinity S, the generic pelagic biogeochemical variable C and the generic benthic variable F : The variables K M and K H represent the turbulent diffusivity for momentum and tracers, respectively, calculated by the turbulence closure model.The parameters χ 1 and χ 2 (Table 1) are the background turbulent diffusivities used to parameterize unresolved processes.The last term in (5) represents the divergence of incoming solar radiation, which has been found to be an important term in the description of surface boundary layer dynamics; term ω S in Eq. ( 6) indicates the parameterized contribution of lateral advection sources, which is a necessary requirement for the maintenance of perpetual climatological vertical structure of the water column.
The most appropriate forms of the parameterization terms are discussed in Sect.3.6.Density is calculated by an adaptation of the UNESCO equation of state proposed by Mellor (1991).Finally, the last term in Eq. ( 7) accounts for the local biogeochemical source and sink terms described in detail in Sect.3.4.Benthic variables in Eq. ( 8) have no transport terms and thus their dynamics are determined by an ordinary differential equation with local reaction terms only, and boundary conditions (fluxes) from the pelagic system are directly included as source terms.Surface and bottom boundary conditions for the momentum Eqs.(1) and (2) are: where U H ≡ (U, V ), τ A is the wind stress, τ b is the bottom stress and ρ W is a reference sea water density (1024 Kg m −3 ).The bottom stress is computed using a logarithmic drag law coefficient and a bottom roughness length of 0.01 m.Boundary conditions for temperature Eq. (5), are where c p is the specific heat of water at constant pressure, Q h is the sensible heat flux, Q b the long-wave heat flux, Q e the heat lost due to evaporation (latent heat flux) and Q corr is a parameterized correction value for balancing the annual heat flux budget already tested in Vichi et al. (1998a), and further illustrated in Sect.3.6.The site-specific computation of the surface heat fluxes and the other forcing function terms introduced above are described in Sect.3.5.During the calibration phase, the seasonal vertical thermal structure of these shallow sites has been found to be highly sensitive to the parameterization of radiative heat penetration, the last term in Eq. ( 5).Downward irradiance is usually parameterized with a depth-dependent exponential attenuation, considering that the infrared portion of the light energy spectrum is completely absorbed at the surface.In shallow coastal areas it is necessary to have a finer resolution of the heating of upper layers, and thus the approximation of imposing the whole infrared radiation at the surface is insufficient.Following Paulson and Simpson (1977), we have applied a partitioning of the incoming short-wave irradiance in a portion extinguished in the first upper layers (simulating the infrared components of the spectrum) and a portion that can penetrate the water column (visible light) as follows: where Q S is the incoming solar radiation flux calculated with the Reed (1977) formula.Therefore, there are two different extinction coefficients, one for the "near-infrared" part (λ i ) and one for the visible (λ v ), with the definition of an apportioning non-dimensional coefficient ε i between the two.The parameters λ i , λ v and the different values of ε i at the three implementation sites (Table 1) have been calibrated through trial and error sensitivity experiments comparing with the observed seasonal vertical profiles of temperature at the three sites (see Table 1 for the values used in the experiments).The adopted values account for the sitespecific environmental characteristics that are not considered in the physical parameterizations.In fact, the value of ε i increases from pelagic (S3) to coastal sites (AA1, S1) because the available transmittometer measurements (F.Frascari and M. Giani, personal communication) indicate a corresponding decrease in the light transmittance related mainly to the influence of coastal waters.The need for a rather high value at S1 is justified by the larger presence of attenuating organic and inorganic suspended matter originating from the Po River (Giani et al., 2001).
The assumption of constant appropriate values for each site is sufficient to provide a satisfactory description of the radiative processes along the water column, as shown in Sect.4.2.This is not completely true for the light used in the biological model, which needs to take into account the extinction due to suspended living particles, because the selfshading effect is an important feature of the phytoplankton blooms dynamics.Thus, the irradiance used as forcing functions for the calculation of production rates is written as:  , d −1 ) 1.0 1.0 1.0 Michaelis-Menten constant for phosphorus uptake by bacteria (h p , mmol P m −3 ) 0.5 0.5 0.5 Optimal N:C ratio in bacterioplankton (n opt , mmol N (mg C) −1 0.017 0.017 0.017 Optimal P:C ratio in bacterioplankton (p opt , mmol P (mg C) −1 0.0019 0.0019 0.0019 C:Chl-a ratio in diatoms P (1)  25 25 25 C:Chl-a ratio in flagellates P (2)  50 50 50 C:Chl-a ratio in picophytoplankton P (3)  50 50 50 where ε P AR (Table 1) is the coefficient determining the portion of Photosynthetically Available Radiation and where P i is the carbon content of phytoplankton groups in the model (P c ), D the carbon content of particulate detritus (R (6) c ) and M the suspended inorganic sediments (ISM); the c constants (see Table 1) are the specific contributions to the total extinction coefficient of each suspended substance.It is important to state that the value of Q S in Eq. ( 14) is evaluated as in Eq. ( 13) as plane irradiance, while for the photosynthetic processes the use of the scalar irradiance should be preferable (Kirk, 1983).We assume that the effect of such simplification is negligible for climatological studies, although it could be of importance when resolving the photosynthetic process with higher temporal resolution.
The surface boundary condition for salinity Eq. ( 6) is a nudging/relaxation scheme, which is written in terms of salinity data at surface S * (t) and relaxation velocity α S : The chosen values of α S at the different implementation sites are given in Table 1.At the bottom, the salinity boundary condition is a simple no-flux condition: Boundary conditions for Eq. ( 7) differ according to the pelagic variable.For functional groups, such as phytoplankton or zooplankton, both bottom and surface boundary conditions are no-flux specifications.For inorganic nutrients (phosphate, nitrate, ammonium and silicate, here indicated as a generic nutrient N ) we use a climatological time series of surface concentrations N * in order to account for inputs of allochthonous nutrients at the surface.The boundary condition is written as for salinity, where α nut is the relaxation velocity given in Table 1.As nutrient bottom boundary conditions we impose the molecular fluxes of the nutrient benthic remineralization processes (Ruardij and van Raaphorst, 1995;also Sect. 3.4) as where ω remin is the dynamically-calculated sediment-water exchange flux at the bottom interface (in mmol of nutrient m −2 d −1 ) and z bot is the depth of the last layer of the vertical grid.Surface boundary conditions for dissolved oxygen (O 2 ) are derived from the air-sea gas exchange velocity (α ox ) parameterized according to Liss and Merlivat (1986) as a function of the wind velocity W A and the dissolved oxygen concentration at saturation.( O sat , Weiss (1970)) in the following way: At the bottom, the boundary condition is equivalent to Eq. ( 19) where the water-sediment flux is a function of the benthic oxygen demand (Ruardij and van Raaphorst, 1995).

The biogeochemical flux model
The ERSEM biological state variables (functional groups) such as phytoplankton, bacteria or zooplankton, are defined as vectors, where each element represents the content of C, N, P and Si in the functional group (Blackford and Radford, 1995).Thus, the dynamics of any pelagic functional group C in (7) or benthic functional group F in Eq. ( 8) is actually written, following the mathematical notation introduced by Vichi (2002), with a vector notation as C i , i = c, n, p, s, where each vector component corresponds to the intra-cellular amount of C, N, P, Si respectively.Particularly, the Si component only applies to the phytoplankton functional group of diatoms.The biogeochemical terms in Eqs. ( 7) and ( 8) are thus expressed as an algebraic sum of source and sink rates: where V indicates the other involved model variables and pr the process affecting the evolution in time of the biogeochemical variable C i , such as uptake, respiration, excretion/exudation, predation etc.Each one of the pelagic state variables listed in Fig. 2 has a dynamical equation involving the compositional elements indicated in the subscript.For instance, there is a differential equation involving the carbon component in phytoplankton and one for each of the four nutrients (nitrate, ammonium, orthophosphate and silicate);, therefore, the internal nutrients-to-carbon ratios in the functional groups have a dynamical variability depending on the external and intracellular conditions.The same is valid for all the other system components, including the benthic variables, which are not shown in Fig. 2, but their contribution to the pelagic dynamics is marked in the scheme with a double arrow on the top of the state variable boxes indicating the benthic-pelagic exchange.

The benthic model
The inclusion of benthic system dynamics is mandatory in shallow shelf areas like the northern Adriatic.The specification of each one of the biogeochemical cycles of macronutrients in the ERSEM modeling system allows one to model the deposition of organic matter with different nutrient contents and its subsequent benthic remineralization controlled by the environmental conditions.The remineralization of organic matter at the sediment-water interface has in fact a long time scale (compared with the water column remineralization time scale) and the process is mediated by the discontinuous availability of settling substrate, the redox conditions and the hydrodynamics of the water column.Therefore it is important to provide the correct buffering capacity of the sedimentary compartment, especially in late summer/autumn, in order to simulate the proper benthic nutrient fluxes that have been found to account for a large portion of nutrient availability in the basin (Giordani et al., 1992).The ERSEM benthic biogeochemical model essentially solves the early-diagenetic equations proposed by Berner (1980), using zero-order remineralization terms derived from the activity of benthic organisms (macro, meiobenthos, filter feeders, aerobic and anaerobic bacteria).The complete description of the benthic biogeochemical nutrient equations used in this work can be found in Ruardij and van Raaphorst (1995), while the description of the processes involving benthic secondary producers and decomposers is given in Ebenhöh et al. (1995) and Blackford (1997).The potential reactivity of sediments is a function of biological activity (bioturbation, bioirrigation) and also of morphological properties, such as the grain size.In early works (Vichi et al., 1998a, b), we have calibrated some of the sediment parameters at the S1 and S3 sites, focusing on the porosity φ and the non-dimensional absorption coefficient for orthophosphate p ads (Ruardij and van Raaphorst, 1995), by means of measurements extracted from the PRISMA-I data set (F. Frascari, personal communication).The same approach has been used for the AA1 site, where the sediment is composed primarily of clay and silt, and the measured porosity has been taken from Cermelj et al. (1997).All the mentioned parameter values are given in Table 2.

Forcing functions
The wind stress forcing function τ A used to evaluate boundary conditions (Eq.10) and the heat flux terms Q b , Q h , Q e for Eq. ( 12), plus short-wave incoming radiation flux Q s in Eq. ( 13), were all calculated from the 6-hours surface reanalyses of meteorological parameters from the European Center for Medium-range Weather Forecast (ECMWF) for the period 1982-93 (Gibson, 1997).The surface fluxes are computed following the procedures described in Maggiore et al. (1998) and Zavatarelli et al. (2002).Surface salinity data S * (t) in Eq. ( 16) for the S1 and S3 sites are taken from the Adriatic BiogeoChemical Data set (ABCD, Zavatarelli et al., 1998) as monthly mean climatological time series considering an adequate area centered around the implementation sites (Fig. 1).Climatological values of salinity at the AA1 site were calculated from a time series collected between the years 1986-1997 (see also Sect.4.2).Another important external forcing function that has been shown to largely affect the simulated biological properties is the Inorganic Suspended Matter (ISM) in the water column (Vichi et al., 1998a).ISM is defined as the suspended sediment fraction, obtained by subtracting the organic fraction from the total particulate collected in sediment traps (F.Frascari, personal communication; Giani et al., 2001).The concentration of suspended sediments modifies the ambient light through extinction processes and consequently decreases the productivity of phytoplankton.The ISM profiles are applied in the model as external forcing functions affecting the light extinction coefficient for biology as in Eq. ( 15); the ISM concentrations at the three sites are computed as linear interpolation of the seasonal profiles shown in Fig. 3. Unfortunately, there is scarce climatological information concerning the seasonal mean concentrations of ISM in the northern Adriatic.At the AA1 site we used observations collected monthly over the period 1997-2000, from which we have calculated the seasonal mean concentration profiles.For S1 and S3, we used the PRISMA-I data collected during only four seasonal surveys.Therefore, those data are not climatologically representative.This additional degree of uncertainty has been taken into account, especially for the S1 site in the Po prodelta area.
In the validation phase (Sect.4.2) we have also performed additional sensitivity experiments, to assess the importance of a good representation of light processes in coastal biogeochemical deterministic models.
Perpetual time series of nutrients at surface (N * in Eq. ( 18)) are climatological mean seasonal values extracted from the ABCD data set in the case of S1 and S3, and seasonal mean values from the LBM data set at the AA1 site.The seasonal frequency was only possible due to the systematic lack of observations in some months.The main differences in the three time series (not shown) are that both S1 and AA1 show high nitrate and phosphate surface concentrations, typical of river-affected coastal stations, while S3 presents lower surface values, indicating oligotrophic characteristics that are proper of open-sea areas.Besides, at S1, all nutrients show a distinct, strong peak in autumn (concentrations are about 10 times more than at the other sites), while S3 and AA1 do not show such a significant autumn increase.

Physical parameterizations
The application of 1-D vertical models for the representation of hydrodynamical processes usually requires the use of additional parameterizations to account for the absence of horizontal dynamics.In the specific case of the northern Adriatic Sea, vertical processes alone -which are determined by the specification of surface heat, momentum and water fluxesare not completely sufficient for an appropriate description of the seasonal evolution of the water column structure.Therefore, it is necessary to include a closure of the annual heat and water budgets, in order to let the model reproduce a perpetual climatogical annual cycle in the hydrodynamics.The annual heat flux budget is negative in the northern basin (Supić and Orlić, 1999), and is likely to be compensated (on a climatological time scale) by the advection of warm waters from the south (Artegiani et al., 1997b).This has been specified in the model, as proposed in Vichi et al. (1998a), by calculating the annual surface heat loss in the forcing functions and distributing this bulk value along the year in the form of an empirical heat correction function added as a surface boundary condition (Q corr in Eq. 12).The most suitable shape of this empirical function at the three sites was established by means of trial and error methods analyzing the simulated seasonal profiles of T with respect to the observed means and their range of variability.
Concerning the closure of the water flux, a crucial problem when dealing with 1-D models is that the local buoyancy losses are not compensated by long-term, basin-wide lateral advection of buoyancy, in order to maintain a perpetual dynamical equilibrium in the water column.The imposition of local net positive or negative heat and water fluxes at the surface produces a model drift; in reality, the water column budget is closed by the horizontal advection processes.In a purely one-dimensional model the horizontal advection terms are neglected, and so it is necessary to parameterize the lateral advective adjustment process.In early works we have applied a closure of the water cycle based on the imposition of a surface salt flux correction, as we did for the heat flux in Eq. ( 12).In this work we use a parameterization of local lateral advection representing the basin-scale contribution to the vertical water column stability.The method consists of the introduction of a climatological time-varying vertical profile of salinity, to which the dynamics of the water column has to adjust within a given time scale.The source term added in Eq. ( 6) is parameterized as where α adv (z) is the depth-varying relaxation frequency (in d −1 ), defined to be 0 from surface down to a given depth δ adv .Values of α adv and δ adv used at the different implementation sites are given in Table 1.The time series of climato-logical profiles S clim has a monthly frequency consistent with surface salinity forcing data S * introduced in Eq. ( 16).

Biogeochemical parameterizations
The main parameterization improvement with respect to the previous modeling work (Vichi et al. 1998a, b), in which the standard ERSEM II implementation was used (Baretta-Bekker et al., 1995, 1997), is the inclusion of dissolved detritus as a dynamical model state variable.In early implementations it was assumed that the Dissolved Organic Matter (DOM) produced at each numerical time step was directly utilized by bacteria, independently from the nutrient-content of the substrate.The implications of such a parameterization were that it was not possible to simulate any accumulation of DOC as observed in the summer period in the northern Adriatic Sea (Pettine et al., 1999(Pettine et al., , 2001) ) and pointed out in Vichi et al. (1998b).Since the modeled dynamics of carbon in phytoplankton is decoupled from the nutrient uptake processes (Baretta-Bekker et al., 1997), the DOM excreted by primary producers can present different degrees of nutrient content in the model.Generally, the utilization of DOM by marine bacterioplankton is dependent mainly on the nutrient availability in the substrate itself (Puddu et al., 2000).More recent investigations have further highlighted the importance of the qualitative character of DOM rather than the quantitative concentrations in the dynamics of this component in the northern Adriatic Sea (Pettine et al., 2001).Therefore, in this work we have implemented a detailed representation of the bacterioplankton dynamics related to DOM utilization: bacterioplankton functional processes now include the concept of refractory organic matter.The degree of refractoriness is determined by the carbon/nutrient ratios of DOM and regulate the bacterial uptake of organic substrate.Optimal uptake is achieved when the C:N:P ratios in DOM and in bacterioplankton correspond to the optimal intracellular bacterial ratio of 45:9:1 (in atoms, Goldman et al., 1987).Bacteria generally require more phosphorus for a proper assimilation into cells with respect to phytoplankton, where the optimal reference ratio is 106:16:1 (Redfield et al., 1963).These values are used in the ERSEM model as reference ratios for the functional dynamics of each biogeochemical component.This assumption does not mean that the chosen Redfield and Goldman stoichiometric models define the regeneration of nutrients in the entire model domain.The overall stoichiometric regeneration, such as the one proposed by Degobbis (1990) for the northern Adriatic, is a final result to which the interplay between the various functional dynamics should tend, and it is not embodied in the model parameterization.Thus, each functional group participates in the actual elemental ratios in the water, according to its internal physiological requirements, and, due to the high degree of connectivity between the defined groups, slight changes in the choice of these ratios have minor effects on the overall model behavior.
The formulation of carbon uptake (BCD, bacterial carbon demand) in the bacteria functional group (B) is written as (Vichi, 2002): where is the BCD regulated by the bacterioplankton physiological state and is the BCD dependent on the dissolved (R (1) c ) and particulate (R (6) c ) detritus "quality".In Eq. ( 24) f T is the temperature uptake regulating factor (expressed with a Q10 exponential formulation), f n,p is the factor determining the health state of bacterioplankton as a function of intracellular ratios f n,p = min q n , q p (26) In Eq. ( 27) p opt and n opt are the Goldman et al. (1987) P:C and N:C intracellular ratios expressed in model units (mmol nutrient/mg C, see Table 2).In Eq. ( 25) the bacterioplankton uptake on dissolved (R (1) ) and particulate (R (6) ) detritus is defined on the basis of the C:N and C:P ratios in DOM and POM, and on "characteristic" time scales for the uptake process (ν R (1) , ν R ( 6) , values are given in Table 2) .The regulating factor f n,p R (j ) defines the trophic quality of R (1) c and R (6) c following a Liebig formulation: This factor tends to 0 when the available substrate is nutrientdepleted, and to 1 in the case of optimal proportions of N and P content with respect to the reference intracellular rates p opt and n opt .It is also important to remember that bacteria are allowed to take-up inorganic nutrients directly from the waters, according to their internal needs (Zweifel, 1993; Baretta-Bekker et al., 1997;Vichi, 2002), as shown below in the case of dissolved phosphorus (N (1) ) uptake: with The utilized parameter values are given in Table 2.This flux will be further used in Sect. 5 for analyzing the microbial system functioning under different physical regimes.The models at the three implementation sites have been initialized with climatological winter profiles of temperature and salinity extracted from the data sets described in Sect.3.5.Climatological initial condition definitions for biogeochemical pelagic components are more difficult to obtain for all the model state variables and it becomes almost impracticable in the case of benthic variables.The latter are generally considered to act as reservoirs for the shallow coastal systems, especially concerning nutrient recycling and availability, and thus their appropriate initialization could be crucial for long-term climatological simulations.The usual methodology applied for the initialization of the benthic system is to start from "educate" guesses based on sparse literature values (e.g.Giordani et al., 1992) and perform a long-term simulation until an equilibrium with the pelagic dynamics is reached.Winter profiles of pelagic biogeochemical constituents in the northern Adriatic basin are usually homogenous due to the mixed conditions of the waters (Zavatarelli et al., 1998).We have performed some sensitivity tests to the initial conditions, initializing the models with different homogenous profiles for the major pelagic state variables, and we found that the pelagic and benthic systems converge to almost the same dynamical perpetual year cycle in 4-5 years of simulation in the presence of the external nutrient input forcing functions described in Sect.3.5.This result may not be a priori generalized to all marine ecosystem models in all areas, but is valid in the case of the ERSEM model in these particular implementations.Therefore, we have initialized the three 1-D models with the same vertically-homogeneous values of biogeochemical variables, and we have analyzed the results of the 6th year of simulation, considering the first 5 years as the specific spin-up time of the benthic system for reaching the equilibrium with the forcing functions.

Validation of seasonal model results
The choice was made to validate model results by computing the seasonal mean vertical profiles of selected variables and comparing them with the corresponding climatological data, as already proposed in Vichi et al. (1998a, b).This choice was motivated by the large temporal inhomogeneity in historical data distribution, which presents monthly gaps and thus, it does not allow the computation of monthly climatologies.Historical observations for S1 and S3 have been extracted from the ABCD data set (Zavatarelli et al., 1998), covering the period 1970-1990 and considering areas of about 0.4×0.4degrees roughly centered around the model locations (Fig. 1).For the AA1 station, we have used the LBM historical data sets over the periods 1986-1990 and 1995-1997, considering casts from the whole Gulf of Trieste.A depth criterium has also been applied to all data extractions, and casts with a bottom depth off the range H ± 5 m have been rejected.The resulting data distributions among seasons, stations and available variables are given in Table 3.
The seasonal means of model state variables have been plotted against the means and ranges of variability (interval between maximum and minimum values found in the time series) of all the available observations in the data sets for the three sites, and we present here only a selection of model results, focusing on temperature (T ), salinity (S), phosphate (N (1) ) and chlorophyll-a (Chla) (Figs. 4,5,6 and 7,respectively).Each figure is divided in rows corresponding to the different sites (S1, S3 and AA1, respectively) and columns that represent the season (calendar seasons).

Temperature and salinity
Figure 4 shows the seasonal profiles of temperature at the three implementation sites.The seasonal cycle is similar at all sites, with well mixed conditions in winter, a thermal stratification starting in spring, and autumn conditions with the characteristic heat storage in the bottom layers.The seasonal observed variability is captured satisfactorily by the simulated profiles.The deeper S3 site shows the best agreement with the climatological data, with the exception of 1-2 degrees discrepancies in winter at the bottom and in autumn at the surface.At S1, the model underestimates the temperature in winter and especially in spring, where the simulated profile is off the range of variability.AA1 has an intermediate behavior that is always within the range of variability, although the modelled profiles are slightly misplaced with respect to the seasonal observed means.
The seasonal profiles of salinity (Fig. 5) show the main differences that characterise each site.S1 has the largest variability in the surface salinity, due to the influence of the Po River runoff.This signal, in spite of surface boundary conditions imposed with a fast relaxation constant (see Eq. 16 and Table 1), is only partially captured by the model; this is particularly evident in autumn, where the surface value is not matched, and the surface freshwater signal is mixed down to a greater depth (with respect to observations) by momentum fluxes.It is important to note that the discrepancies mentioned above in the description of temperature profiles (Fig. 4) occur in the seasons when the salinity mean profiles are in worse agreement with the observations, indicating that the introduced parameterizations of local advection of salt (Sect.3.6) are not completely sufficient to describe the seasonal system dynamics.This is especially true in winter and autumn at S1 and AA1.In the Gulf of Trieste station, in fact, the model is not completely capable of capturing the vertical location of the halocline and the upper layers salinity structure in general.In winter, the halocline is placed at a depth of 15 m, while the data show its location at about 5 m and a similar misfit is visible in autumn.The 1-D approach and the imposition of surface salinity data appear not to be sufficient to completely contrast the wind mixing and to maintain the strong stratification observed at the AA1 site, although they do provide a reasonable representation of the seasonal variability and the coastal features of the area.Besides, the model-data discrepancies at S1 and AA1 can also be attributed to the fact that observations at these sites are less representative of a climatological situation, due to the high dynamical variability of such coastal areas.

Biogeochemical pelagic variables
The simulated seasonal profiles of phosphate in the more pelagic S3 area (Fig. 6) are within the range of variability in all seasons, except at the bottom in winter and spring.This discrepancy with the observed concentrations could be due to the inclusion of observations taken at sites with a depth greater than 30 m, where the concentration increase at the bottom is located deeper in the water column.A similar behavior is also found in the coastal station S1.The winter phosphate concentration predicted by the model at S1 shows a clear increase with depth (as found in the observations), but the concentration of orthophosphate is overestimated.The excessive bottom nutrient concentrations during winter in the model results are also observed in the other nutrient profiles (not shown).This is not likely to be caused by an overestimation of autumn remineralization processes, because at that time the overall concentrations at the bottom are comparable with observations, but rather by the lack of nutrient removal due to advective processes that are particularly intense in the Po prodelta area.
Another cause can be the underestimation of biological uptake due to improper inorganic suspended matter boundary conditions determining an unfavorable light environment.Since the contribution of advective processes cannot be taken into account in a purely 1-D model, we have tested this second hypothesis by performing a sensitivity analysis to the ISM forcing function.The ISM winter concentration at S1 (Fig. 3) is the highest with respect to all the other seasons and sites.As noted above in Sect.3.5, this profile is not derived from climatological observations, but has been collected during a single winter campaign.Therefore, we have applied the ISM seasonal profiles used as forcing functions at AA1 (Fig. 3), which are still proper of a coastal area, but are derived from a three-year monthly time series.The results of the seasonal phosphate profiles obtained with the sensitivity run at S1 have been superimposed to the standard run in Fig. 6 using a dashed line.It is clear that changing the ISM concentration affects the phosphate vertical distribution, and this is due mainly to an increase in the winter phytoplankton standing stock that utilizes more nutrients (see Fig. 7 and paragraph below).The phosphate profile falls within the observed winter range of variability at the surface and becomes much closer to the observed mean profile than in the standard run.The same behavior is observed in spring, while summer and autumn do not show any significant change.The sensitivity of the model to an appropriate representation of the vertical light distribution shows, on the one hand, the importance of having data concerning the ISM distribution and, on the other hand, the necessity to include more sophisticated parameterization of light processes in deterministic models.
Spring and summer profiles are in satisfactory agreement at surface at all sites except AA1, where the surface value is underestimated, due to the strong uptake by phytoplankton during wintertime.The winter AA1 phytoplankton biomass (Fig. 7) is in fact the largest one with respect to the other site, and this can be explained by the favorable earlier onset of the stratification occurring at AA1, as discussed in the next section.
The comparison of simulated phytoplankton biomass is shown in Fig. 7 as seasonal profiles of chlorophyll concentrations.The model state variable representing Chlorophyll-a (Chla) is a diagnostic variable in the model, which is derived from phytoplankton carbon content using the fixed ratios given in Table 2.The simulated profiles fall within the observed ranges in almost all seasons, especially in AA1, although there is a systematic underestimation in S3, particularly during spring.Winter profiles can be conveniently compared with sufficient climatological data only at AA1, where the climatological model predicts an overestimation with respect to the observed mean profile; concentrations, however, still lie within the observed ranges.
On the contrary, at S1, the limited amount of observations indicate a possible underestimation of the biomass stock during winter, which can partially justify the high nutrient concentrations given by the model.Given the model equations for phytoplankton (Baretta-Bekker et al., 1997;Vichi, 2002), two reasons can explain the low biomass in the model at S1, in spite of the high nutrient availability: (a) light limitation due to high ISM concentration; (b) unfavorable mixing conditions which move the phytoplankton in and out from the euphotic zone, limiting the adaptation to light.Actually, these two factors are tightly coupled to each other, and the occurrence of the spring bloom is strictly determined by the interaction of the vertical light distribution and the hydrodynamics of the water column, as will be discussed in detail in the next section.The sensitivity experiment on the ISM forcing function described above has confirmed that the winter phytoplankton mean concentration is strongly affected by light limitation.The chlorophyll-a profile, computed from the sensitivity run, forced with the ISM climatological profiles used at AA1 (dashed line in Fig. 7) is more in agreement with the observations with respect to the standard run (continuous line).Moreover, the subsurface maximum developed by the model in the standard run during spring is completely dependent on the light conditions; this feature is slightly visible in S3 as well and absent in AA1, while none of the observed mean spring profiles show any clear subsurface maximum.In spring, the S1 site is already stratified (see Sect. 4.3 and Fig. 9a) and the occurrence of subsurface maxima is directly correlated to low concentrations of ISM (Fig. 3).As soon as surface ISM values are increased (sensitivity test not shown), we observe the shallowing of the maximum and the outcropping to the surface, while, with a decrease in the ISM concentration, as obtained by applying the AA1 ISM spring profiles (Fig. 3, corresponding to the dashed line profile in Fig. 7), the maximum is moved deeper in the water column.Therefore, these experiments with an improved biogeochemical model confirm the importance of the knowledge of proper ISM profiles for a good representation of primary production processes, as first suggested in a previous work (Vichi et al., 1998b).

Benthic remineralization
The use of a comprehensive coupled ecosystem model allows us to investigate the role of certain components in the determination of the overall seasonal characteristics of the northern Adriatic system.In particular, we decided to focus on the benthic-pelagic remineralization fluxes, which are considered a substantial component of the nutrient cycle in this shallow shelf sea (Giordani et al., 1992;Spagnoli and Bergamini, 1997).The sensitivity of the northern Adriatic Sea nutrient dynamics to the presence of an active benthic system has been tested with the model at the AA1 site, where we have the highest availability of nutrient observations (Table 3).We have performed a simulation without the benthic model, simply prescribing a removal of sinking particles from the model domain when they are deposited at the bottom.Our intention was to show mainly the importance of having a benthic model that dynamically responds to the inputs from the pelagic system by (re)supplying nutrients through remineralization of organic matter in the sediments.We are particularly interested in phosphate dynamics, because this component is generally referred to as a limitation for phytoplankton production in the Adriatic Sea.In previous sections we have shown that simulated phosphate profiles at AA1 are generally within the observed range of variability (Fig. 6), and we want to assess how much of the dissolved concentration comes from benthic regeneration fluxes.
Figure 8 shows the results of the sensitivity experiment without benthic remineralization for state variable phosphate (N (1) ) at the bottom, compared with the standard model time evolution and the observed seasonal means from the LBM data set.The effect of benthic fluxes are visible mainly during late summer (August-October), in agreement with benthic chamber measurements (Giordani et al., 1992;Spagnoli and Bergamini, 1997) that register the maximum fluxes during this period.This gives indications that the dynamical model is capable of developing an adequate time lag between the spring deposition of organic particles and the release of inorganic components.Benthic fluxes appear to contribute to about half of the PO 4 concentration in the bottom layers during summertime, and this is also valid for all the other modeled inorganic nutrients (not shown).Nevertheless, there is a clear discrepancy in winter and autumn concentrations, where we have, respectively, an over-and underestimation of the observed mean values.During these periods, nutrient dynamics are tightly connected to the mixing state of the water column, and the comparative analysis of results hints at a possible reduced mixing during winter, leading to the segregation of nutrients in the deeper layers.However, winter overestimation could also be due to the lack of horizontal advection processes, as previously invoked in the case of the other coastal site S1.The autumn underestimation could instead be due to low buoyancy fluxes that enhance the turbulent diffusion of nutrients from the bottom layers.This is particularly evident if we look at the modeled autumn salinity profile at AA1 (Fig. 5), which shows a greatly reduced vertical gradient with respect to the observations.

Physical-biological interactions
In this section model results are comparatively analyzed in order to explain some of the different features observed at the three sites and highlighted in the previous section.Particularly, we focus on the seasonal cycle of primary producers as affected by the seasonal vertical dynamics of the water column.Figure 9 shows the depth-averaged seasonal cycle of state variable Chla at S1, S3 and AA1 sites, plotted against the corresponding cycles of the mixed layer depth and of the critical compensation depth (D cr ).The Chla means are computed over the depth-range comprised between the surface and D cr , which is defined (following Sverdrup, 1953) as the level where the vertically-integrated primary production rates equal the integrated autotrophic respiration rates.The dynamical time evolution of D cr at the 3 sites shows the largest variability at S1 (Fig. 9a) progressing from nearsurface (winter) to the bottom (summer), while has a more gentle deepening at S3 (Fig. 9b) but with a sharp shallowing in October.Site AA1 (Fig. 9c) has an intermediate behavior, with a longer winter shallow value, a steep deepening in March/April and a recovery similar to S1 in late summer/autumn.One of the main factors driving the vertical location of the D cr is the underwater light climate, which, as shown in the previous section, is strongly affected by the presence of light-attenuating suspended particles.The shallower location of the D cr in winter at S1 can thus be explained by the high concentration of ISM with respect to the other sites, which limits the penetration of light in the deeper layers.The seasonal cycle of Chla follows the typical behavior of mid-latitude temperate coastal areas, with a distinct spring bloom, a summer decay and an autumn bloom, particularly intense at S1 (Fig. 9a).The onset of the "spring" bloom is slightly shifted at the three locations: it occurs in late winter at S3 and AA1, while we find it in early spring at S1.These results are in partial agreement with the 30 year long-term analysis by Degobbis et al. (2000) of surface Chla data collected in a transect from the Po Delta to the Croatian coast.However, their open-sea site does not show any clear evidence of an earlier spring bloom with respect to coastal sites, as seen in S3 (Fig. 9b).This early winter bloom is, as mentioned above, more pronounced in AA1 (Fig. 9c), which is partially in contrast with climatological observations, although, due to the high interannual variability observed in the Gulf of Trieste, there have been observations of phytoplankton spring development from late February to June.
The seasonal primary producer cycle simulated by the model shows distinct differences in the 3 case studies, and they can be interpreted as consequences of the coupling between the environmental factors -represented mainly by the light availability -and the hydrodynamical conditions that develop at S1, S3 and AA1.If it can be demonstrated that the timing of the bloom in the basin is dependent mainly on the shallowing of the mixed layer and the deepening of the euphotic zone, according to a Sverdrup (1953) mechanism, then it is important for a good predictability skill to have appropriate predictions of the surface mixing processes driven by momentum and heat exchanges.At the same time, due to the intrinsic coupled nature of this process, it is also required to provide a good estimate of the underwater light climate as determined by the amount of suspended matter in the water.
To verify if the bloom timing is strictly coupled with the hydrodynamical conditions, we have computed the square of the Brunt-Väisälä frequency (N 2 ) from model results and outlined in Fig. 9 the time series of the depths where N 2 is maximum and also greater than 0.5 10 −4 rad 2 s −2 (below which these shallow water columns can be considered neutral).The obtained depth is assumed to be indicative of the extension of the mixed layer (the bottom of the mixed layer should always be shallower of this depth), and has been superimposed to the D cr time evolution in Fig. 9.We observe that the model is capable of reproducing the main dynamical characteristics of the three sites introduced in Sect.2, simulating a continuously stratified water column at S1 (due to the influence of the Po River freshwater input), a winter mixed layer reaching the bottom depth at S3 and an intermediate situation at the other coastal site AA1 (determined mainly by the interaction between the lower river inputs and the stronger winds).
At all three model implementation sites it can be noted that the chlorophyll increase in winter/spring is matched by the reduction of the mixed layer thickness and the deepening of D cr .The biomass peaks are achieved when the lower limit of the mixed layer is either shallower than D cr (S1 and S3) or has a comparable depth (AA1).This leads to the conclusion that at the three sites the winter/spring bloom onset is indeed controlled by the relative location of D cr (driven by light penetration processes) with respect to the lower limit of the mixed layer (driven by turbulent dynamics), as originally observed by Gran (1931) and quantitatively described by Sverdrup (1953).Given the model results at the three locations shown in Fig. 9, we can speculate that in shallow areas, it is apparently not necessary that the depth of the mixed layer be much less than the critical depth in order to stimulate the major growth in phytoplankton biomass, as suggested by Sverdrup (1953).It appears sufficient to have a small increase in the stability of the water column to enhance the phytoplankton production, as shown by the time-matching of the gradient change both in the biomass curve and in the mixed layer depth of Fig. 9.
During summer, D cr is always deeper than the lower limit of the mixed layer, but the water column is nutrient depleted (see also Sect. 5 and Fig. 12b) and phytoplankton biomass is at a minimum.In autumn, a recovery of the phytoplankton biomass occurs.The biomass increase is particularly evident at S1, where the autumn bloom has a magnitude comparable with the winter/spring one.The autumn breaking of the stratification, marked by the deepening of the mixed layer, brings bottom benthic remineralized nutrients back in the euphotic layers of the water column, enhancing the primary production.This increase of nutrient availability is amplified at S1 by the input of nutrients at the surface prescribed in the boundary conditions, giving rise to the simulated high biomass peak.
In order to further confirm the strict dependence of the spring bloom mostly on just the mechanical mixing, we have applied the nutrient surface boundary conditions used at AA1 as forcing functions of the model at S1.The two sites have the same depth and equivalent distance from the coast, but nutrient concentrations observed at AA1 are generally much lower than S1, except during the spring period.Therefore, by changing the surface boundary conditions in this way, we prescribe a higher nutrient availability during the bloom period that should lead to more favorable growth conditions.However, the timing of the winter/spring bloom (not shown) is completely unaffected by the changes in the nutrient import in the system.In contrast, the bloom duration and further the summer/early autumn developments show a different behavior with respect to the standard run.This indicates that the system has a delayed response to these conditions which shows up during the more stratified period, when the related benthic remineralization and biological interactions act as positive feedbacks that enhance the primary production.
Obviously, the results described above have general validity only in a climatological context and in the absence of biological control on the physical dynamics, as in our case where the radiative heat penetration is not affected by suspended particles.The feedback that can occur between the biological events and the physics of the water column have recently re-opened the discussion on the Sverdrup mechanism (Townsend, 1992;Stramska and Dickey, 1993) and it is important in the future to explore the effects of this cou-Table 4. Comparison of model simulated rates of primary production (PP), bacterial carbon production (BCP) and bacterial abundances (BA) with observations.Field observations are three-day averaged values integrated along the water column from Table 4 in Puddu et al. (1998).Model data are monthly averages integrated on the water column at each model implementation site (S1 and S3).Averages are given with the coefficient of variation (cv%) pling and the response of phytoplankton productivity to highfrequency changes in the mixed layer dynamics (Vichi, 2000(Vichi, , 2002)).

Comparison with available in situ rates
In recent years, the availability of in situ measured biological rates has increased, and it is essential that a biogeochemical flux model is validated not only with observed concentrations, but also in terms of fluxes of matter along the trophic web.At S1 and S3, we have used the primary production, bacterial production and bacterial abundances collected during the PRISMA-I cruises (Puddu et al., 1998) for a comparison/validation with corresponding model simulated rates and concentrations.It is important to state that the comparison can only be done as an order of magnitude, because observed rates were estimated as two-day averages of seasonal-frequency measurements from a specific year, while model rates are given as monthly climatological averages of the corresponding month of sampling.Table 4, modified from Puddu et al. (1998), shows the vertically integrated daily rates of primary production (PP, L. Alberighi, CNR-Venice, personal communication), bacterial carbon production (BCP) with the associated coefficients of variation, and mean water column bacteria abundances (BA) at S1 and S3 sites during the PRISMA-I campaign of April, July, October 1995 and January 1996.The observed PP has been compared with the mean monthly gross PP calculated by the model and integrated along the water column.BCP in the model has been computed as net production by subtracting the respiration term from the bacterial carbon demand in Eq. ( 23).This is assumed to be consistent with the rates computed by Puddu et al. (1998) in the experimental observations.Bacterial abundances have been derived by vertically-averaging the bacterial carbon concentration and applying a conversion factor of 20 fg C cell −1 (Lee and Fuhrman, 1987).
Considering the coefficients of variation given in Table 4, the model computed PP is generally of the order of magnitude of the observations, except in January at S1, where the value is underestimated by one order of magnitude.In addition, S1, besides the maximum production is found in July with a diminution in October, while the observations report an increase from July to October.There is, however, a good agreement in the relative behavior at the two sites: S1 has rates 2-3 times higher than S3, as found in the observations (except in winter).
BCP rates at S1 match the observations in April, July and October, partially indicating the validity of the new parameterizations described in Sect.3.7.The January mean in the model is 4 times lower, but with a large coefficient of variation with respect to the mean.The low winter bacterial activity in the model is clearly related to the low PP rates discussed above.In fact, the winter BCP/PP rate simulated by the model (∼ 50%) is lower than the observed value of 70%.This discrepancy might be explained with the absence of direct river-borne DOM input in the model, which should enhance the bacterial carbon demand, as it is thought to happen under the influence of the Po River (Puddu et al., 1998;Pettine et al., 1999).
Concerning the bacterial production at S3, we observe an overestimation with respect to the measured values, although the order of magnitude is comparable.The higher bacteria activity in S3 is also reflected in the values of bacterial abundances, which are about twice the observed values in winter.The vertical distributions of bacteria (not shown) are in good agreement with the observations of Puddu et al. (1998) and, as a general consideration, the model seems to partially confirm the hypothesis that the variability of bacteria is more pronounced in their activity rates than in their biomass (La Ferla et al., 1998;Puddu et al., 1998;Zaccone et al., 1999).In fact, for an approximate doubling in the model BCP from April to October in S1 (Table 4), the simulated abundances increase only by a small percentage (from 0.99 to 1.02 10 −9 cell l −1 ).
Biological rates at AA1 sites have been compared with surface PP data and vertically-integrated BCP data for the period 1998-2000.Figure 10 outlines the measured monthly rates for three distinct years plotted against the monthly averaged surface gross PP computed by the climatological model.Since the model PP is given as a daily mean, using the parameterization described in Ebenhöh et al. (1997, see also Vichi, 2002), it has been converted in the units of observations, taking into account the day-length according to the Dobson and Smith (1988) formulation.The time evolution of the simulated PP follows the behavior depicted in the observational time series, although it is evident that the high interannual variability of plankton activity in the data of the Gulf of Trieste (Fonda Umani et al. 1995;Mozetič et al., 1998) cannot be reproduced with a climatological simulation.Longer time series (at least 5 years) of observations are needed for a better estimation of the climatological variability of primary production rates and for a better comparison and validation with climatological model results.However, the model captures the mean magnitude of the spring bloom production peak and also the lower amplitude of the autumn recovery of phytoplankton production, indicating a satisfactory functional behavior of the model at AA1 as well.
The seasonally averaged BCP (over the period 1998-2000) measured at AA1 is compared with the model simulation in Table 5.Here, the simulated daily BCP has been converted to observation units dividing by the number of hours per day, because bacteria activity can be assumed to be not directly dependent on the light cycle.The annual evolution of BCP is reproduced correctly by the model with a progressive increase along the year and a stabilization in summer and autumn.However, the predicted values are twice the observations (although the coefficient of variation in the measurements is rather high), and the spring decrease observed in the data is not matched by the model.We have performed some sensitivity experiments in order to explain model discrepancies with respect to observations.The dynamical value of BCP in the model is dependent on both the bacterial carbon demand in Eq. ( 23) and the bacterial growth efficiency (BGE) parameter as detailed in Baretta-Bekker et al. (1995).The value of bacterial efficiency considered appropriate for the Gulf of Trieste is 0.3 (S.Fonda Umani and F. Azam, personal communication), while in the model we use 0.4, which is the average of the values used by Baretta-Bekker et al. (1995, 1998).A sensitivity experiment on this parameter has shown that simulated BCP can be reduced by about 17% in winter and 23% in the other seasons, with a reduction of the efficiency of 25%, as shown in the last column of Table 5.However, the general behavior of BCP is not modified, and the model is still not able to explain the decrease in BCP observed in the spring average value.Thus, we have performed an additional sensitivity experiment by applying a recently proposed empirical parameterization of the BGE dependency on the environmental temperature by Rivkin and Legendre (2001).By analyzing several observations in the world ocean, they have found a significant inverse linear relationship between BGE and temperature, and the constant BGE has been modified as Bacterial dynamics in the model already have a temperaturedependent respiration, but this is considered to be only a function of bacterial biomass (rest metabolic respiration).With Eq. ( 31), we also introduce an inverse temperaturedependence on the activity respiration rates, which are computed according to the organic substrate uptake in Eq. ( 23) (Baretta-Bekker et al., 1997;Vichi, 2002).
The simulated BCP time-evolution obtained with the improved parameterization is much more satisfactory and is shown in Fig. 11, compared with the model result obtained with the constant BGE of 0.3 and the observations from Table 5.There is a slight decrease in the winter bacterial activity, and the spring average value is still not within the observed range of variability, but the overall behavior of the modeled BCP is more in agreement with the observations compared to the model results with a constant BGE.Thus, temperature seems to be a determining factor in the formulation of bacterial efficiency during high DOC availability, and it is necessary to perform more experiments, in order to assess the consequences of such an effect on the other system components.Supplemental investigations and longer time series of BCP measurements are needed to refine the knowledge of long-term variability of bacterial activity and the related parameterization issues.Indeed, there is a clear sensitivity of the model to the parameterization of microbial dynamics, and these processes are key factors that control the evolution of DOM in the northern Adriatic, as will be further discussed in the next section.

Trophic interactions and DOM dynamics
In an early work (Vichi et al., 1998b) we analyzed the carbon pathways in an aggregated structure of the pelagic food web of the standard ERSEM II model with climatological simulations at the same S1 and S3 sites described in this research.
In that analysis we found that the role of bacteria in channeling the organic carbon was comparable to the classical herbivorous flux from phytoplankton to zooplankton.This was found to be valid not only at S3, which, having more opensea features, is likely to be driven by internal remineralization processes, but also at the S1 site that, on the contrary, should be dominated by new production and particulate matter formation.Therefore, we concluded that the multivorous food web concept proposed by Legendre and Rassoulzadegan (1995) is also a possible food web pattern in the coastal areas of the northern Adriatic with substantial allochthonous input of inorganic nutrients.
In this research we have extended the previous work in order to obtain further insights on the ecosystem functioning, focusing on the central role of DOM and the bacteriaphytoplankton coupling that has been improved by introducing the detailed dynamics of bacterial DOM utilization given in Sect.3.7.The assessment of the model capabilities in the description of microbial processes is an important step towards the development of predictive ecosystem models, especially in the light of recent findings in the microbial ecology of the Mediterranean Sea (Azam et al., 2000).The DOM dynamics in the northern Adriatic is of great concern due to the mucilage formation events that occurred in recent years (Azam et al., 1999;Degobbis et al., 1995Degobbis et al., , 1999;;Puddu et al., 2000 and references therein).Recent data on DOC distribution in the northern Adriatic Sea (Pettine et al., 1999(Pettine et al., , 2001) ) show a tendency to accumulate during the summer period, and a high concentration of DOC is expected to be a pre-conditioning of the mucilage formation.
Figure 12a shows the climatological time evolution of the modeled DOC concentration averaged within the critical compensation depth D cr (see Sect. 4.3) at the three implementation sites.The model predicts substantial DOC accumulation during the spring/summer periods at all locations.In the S1 area, the observed summer values in the period 1996-1997 ranged from 792 to 3372 mg C m −3 (Pettine et al., 2001), which are in good agreement with the simulated values.However, observations also indicate the presence of Table 5.Comparison of model simulated rates of bacterial carbon production (BCP) with observations at AA1 site (Gulf of Trieste).Model results and data are computed as seasonal averages integrated along the water column and given as a mean with the coefficient of variation (cv%).Results from simulations with a bacterial efficiency of 0.4 (standard value in the experiments) and 0.3 are shown

Season
Observed BCP Modeled BCP eff.= 0.4 Modeled BCP eff.= 0.3 a background level of DOC (about 800 mg C m −3 ) that persists throughout the winter periods.This feature is not captured by the model, which shows a complete consumption of the summer stock during the autumn period.This hints at a limited effect of the parameterization described in Sect.3.7, which indeed is capable of providing the observed summer build-up, but does not introduce an adequate aging mechanism within the DOM pool that renders a portion of it more refractory to bacterial attack.In fact, when the autumn conditions are established, new labile DOM of optimal quality is added to the pool that modifies the nutrient-content ratios, making it a good substrate for bacterial growth.Probably, there is also a contribution of refractory DOC from the Po River during this period (Pettine et al., 1998) that maintains the DOC quality at an intermediate level.This process needs to be further investigated in future works and here we focus mainly on the description of the summer accumulation mechanism.
Accumulation starts as soon as the stratification is established (see Sect. 4.3) and has a clear and different behavior at the three sites.The accumulation mechanism developed by the model is a consequence of the coupling between the parameterized biological functional processes and hydrodynamics: the onset of the stratification enhances the phytoplankton production and the depletion of dissolved inorganic phosphorous (DIP, Fig. 12b).As a consequence of the decoupling between phosphorus and carbon dynamics in phytoplankton (Baretta-Bekker et al., 1997;Vichi, 2002), the unbalanced nutrient conditions inhibit carbon assimilation into biomass, and the photosynthesized carbon is exududed by autotrophs in the form of nutrient-impoverished organic matter.A large portion of this standing DOC is probably composed of transparent exopolymers (TEP, see for example, Mari et al., 2001), although the model still does not account for this differentiation.This exudation flux leads to dissimilar DOC accumulation rates at the three sites (Fig. 12a): AA1 and S3 show an early build-up due to the earlier shallowing of the mixed layer discussed in the previous section, but AA1 reaches higher values due to the greater production rates typical of coastal zones.The maximum value at the other coastal site (S1) is comparable to AA1, but shows a more dynamical evolution towards a later summer maximum, indicating the formation of more complex biogeochemical pathways of the organic matter flow in the trophic web.
In order to thoroughly understand the biogeochemical dynamics developed by the model and their interactions with the hydrodynamical regimes, we have analyzed the model results in terms of fluxes among system components by computing a set of ratios reflecting the status of the planktonic food web, as proposed by Legendre and Rassoulzadegan (1995).We have modified the original indices slightly, in order to adjust to a phosphorus-limited system instead of a nitrogen-limited one.Our principal aim here is to investigate the model behavior to see if the deterministic nonlinear coupling of simple empirical biogeochemical relationships with the hydrodynamical transport allows the emergence of features that are indeed observed in marine ecosystems.The first considered measure (Fig. 13a) allows us to distinguish among the different trophic web patterns developed by the model.It is computed as the ratio between the transfer flow of carbon from autotrophs to zooplankters (sum of carbon grazing rates of all the micro-and mesozooplankton groups upon phytoplankton state variables) and the equivalent from bacteria to zooplankters (both computed as an average within the D cr ).The value of this ratio indicates if carbon is transferred to higher trophic levels through the classical herbivorous chain (greater than 1), the microbial loop (lower than 1), or, in the case of values close to 1, through a multivorous food web.An additional indicator is given by the "quality" of DOM as a proper substrate for bacterial growth, which we define as the DOC:DOP ratio (R (1) c /R (1) p , as described in Sect.3.7).The time series of this ratio, averaged within the D cr , is shown in Fig. 13b.The third index illustrates the role of bacteria with respect to the uptake/release of phosphorus.This is the phosphate flux between bacteria and the DIP pool given in (29), defined positive when bacteria release DIP and negative when bacteria compete with phytoplankton for dissolved inorganic phosphate.The integrated value of this flux within the D cr is shown in Fig. 13c.
The type of grazing (Fig. 13a) is a function of the available resources for zooplankton.In winter/spring, when the phytoplankton spring bloom develops at all sites, the transfer of carbon is mainly through the herbivorous food chain, as indicated by ratios higher than unity.In the decaying portion of the bloom, which is sharper at S1 and with a more gentle slope at S3 and AA1, the autotrophic phase is followed by the development of a microbial web, as also observed in many realistic meso-and eutrophic ecosystems (Legendre and Rassoulzadegan, 1995).The quality of the substrate (Fig. 13b) deteriorates after the bloom, and this determines the shifting of the microbial community from nutrient remineralizers to competitors for inorganic resources, as depicted in Fig. 13c.The microbial food web pattern (indicated by values lower than 1 in Fig. 13a) persists for the whole summer period, with a more steady behavior in S3 and AA1 than at S1. Particularly, S1 shows a series of oscillating phases which are also reflected in the behavior of the DOM quality (Fig. 13b), indicating that the input of nutrients from the river and/or more active remineralization processes give rise to small recoveries of the multivorous transfer pathway.
In autumn, at the breaking of the stratification, the three sites again show a different behavior.S1 develops a clear multivorous food web with the microbial "grazing" comparable to the herbivorous grazing, while S3 and AA1 remain in the status of a microbial through-flow of carbon.Thus, the quality of DOM released in the environment improves more rapidly at S1 than at S3 and AA1 (Fig. 13b), changing the activity of bacteria to phosphorus remineralizers (Fig. 13c).Examples of marine systems with these alternating trophic pathways can be found in Legendre and Rassoulzadegan (1995) both for coastal and open-sea areas.Model results suggest that the concurrent existence of different mass-transfer patterns takes place in the shallow northern Adriatic Sea as well, and the shifting from one to the other is essentially modulated by the physical conditions of the water column, the underwater light climate and by the external nutrient inputs.

Summary and conclusions
In this paper we have shown results from numerical simulations carried out with a 1-D complex biogeochemical fluxes model implemented at three different locations of the northern Adriatic shelf.Emphasis is put on a comparison with in situ measured rates, and on the functioning of the northern Adriatic ecosystem.The work has to be considered as preliminary to the development of three-dimensional numerical simulations of the ecosystem behavior.
The comparison with the seasonal observed means described in Sect.4.2 confirms the model capability of capturing the major seasonal local biogeochemical and physical dynamics at the three implementation sites.It is widely accepted in oceanography that much of the variance of the marine biogeochemical component behavior is directly determined by the abiotic dynamics.Indeed, we have found a more satisfactory representation of observed climatological means of biogeochemical state variables at the S3 site, where the simulated physical variables have shown the best (a) ratio between the carbon flow due to herbivorous grazing (from autotrophs to zooplankters) and the one due to microbial grazing (from bacterioplankton to zooplankters; in semi-logarithmic scale).(b) ratio between the C-component and P-component in dissolved organic matter (DOC:DOP), indicating the "quality" of substrate for bacteria growth.The optimal ratios for phytoplankton (106:1; Redfield et al., 1963) and bacteria (45:1; Goldman et al., 1987) are marked in the plot.(c) phosphorus flux between bacteria and dissolved inorganic phosphorus (DIP).agreement with observations.The climatological model behaves better in the S3 sub-region because the short-scale variability induced by coastal processes and riverine input is reduced.Discrepancies at the more coastal sites (S1 and AA3) can be attributed to the lack of this variability in the climatological boundary conditions, and to the absence of horizontal advective processes that cannot be taken into account with a 1-D model.Sensitivity analyses performed on the inorganic suspended matter forcing functions have shown how model results are substantially affected by such external boundary conditions, and that the observed seasonal signal cannot completely be reproduced by just applying climatological means of this dynamical component.This strengthens the need for coherent higher frequency data sets to improve the predictability skills of deterministic biogeochemical flux models.
The onset of the winter spring bloom at the three implementation sites appears to be driven by the local evolution of the water column stratification conditions in relation to the underwater light climate (according to a Sverdrup-like mechanism); the different timing of the bloom development also seems to condition the subsequent evolution of the microbial web dynamics.Results give indications that this coupled mechanism is the principal driver of phytoplankton blooms when using climatological forcing functions, and we suggest that the high variability observed in these coastal areas is composed of this background process modulated by the shorter scale signals of other external boundary conditions.The input of nutrients did not seem to affect the timing of the bloom, but we expect that the imposition of realistic higherfrequency inputs of nutrients from the river could have a substantial effect on this process.
Comparisons with observed biogeochemical rates, such as primary and bacterial production, are indicating -within the limit of the climatological implementation -the general qualitative agreement of model results with observations.Simulations also confirm the model skill to adjust to different trophic conditions, being capable of accommodating high productivity rates typical of areas under direct influence of river input (S1 and AA1) and more oligotrophic offshore regions (represented here by the S3 site).The partial agreement of the simulated bacterial carbon production with observations encourages the inclusion of more sophisticated deterministic parameterizations of microbial web processes (such as the differential DOM utilization proposed in this work), which are expected to have a large impact on the simulation of organic matter transfers in the ecosystem.In particular, the temperature dependence on the bacterial growth efficiency has been found to be an important factor for improving model skills in reproducing the observed microbial production rates.
Analyses of model results carried out in the framework of the conceptual scheme of ecosystem functioning proposed by Legendre and Rassoulzadegan (1995) seem to confirm that also for a coastal shelf basin such as the northern Adriatic Sea, the major carbon transfer pathway can shift from the classical herbivorous food web to the microbial one, and the development of a multivorous food web is likely, in particular at locations such as S1, to be characterised by strong nutrient external inputs.However, our simulations also indicate that after the winter spring bloom, during the summer season, the main carbon pathway is generally the microbial one.The important role of the microbial community revealed by the model is in agreement with the conclusion of Puddu et al. (1998) and Pettine et al. (2001) based on direct observations, which suggest that a large fraction of primary production is processed within the microbial loop.Therefore, the dynamics of DOM should be regarded as a central issue for a complete understanding of the northern Adriatic ecosystem functioning.

Fig. 1 .
Fig. 1.Map and bathymetry of the northern Adriatic shelf with location of the model implementation sites (depths in meters).Dashed lines delimit the maximum extension of the data extraction areas from the historical data sets.

Fig. 3 .
Fig. 3. Seasonal vertical profiles of Inorganic Suspended Matter (ISM) concentration used as external forcing functions at the three implementation sites.

Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 4. Model and data comparison for temperature as climatological seasonal profiles at S1, S3 and AA1.The continuous line is the simulated mean seasonal profile for each site, and season and climatological observations are plotted as seasonal means within the range of variability.

Fig. 7 .
Fig. 7.Model and data comparison for biogeochemical model state variable chlorophyll-a (Chla) as climatological seasonal profiles at S1, S3 and AA1.The continuous line is the simulated mean seasonal profile for each site, and season and climatological observations are plotted as seasonal means within the range of variability.The dashed line profiles at S1 (first row) derive from a sensitivity experiment on the concentration of Inorganic Suspended Matter applied in the model as external forcing functions.The continuous lines are the model results with the standard ISM profiles for the S1 site described in Sect.3.5, while the dashed lines are the seasonal means computed from a run in which the AA1 ISM climatological profiles (Fig.3) have been imposed.

Fig. 8 .
Fig. 8. Bottom concentration of phosphate at AA1 in the standard model and in the simulation without benthic remineralization compared with the LBM data.Data are shown as seasonal means with ranges of variability.

Fig. 9 .
Fig. 9. Biological-physical interactions at the three implementation sites: (a) S1, (b) S3 and (c) AA1.Thick continuous line is the model simulated chlorophyll concentration (in mg m −3 ) averaged within the critical compensation depth D cr (defined as the depth where vertically-integrated primary production rates equal autotrophic respiration rates).The dynamical evolution of D cr (dashed thin line) and of the mixed layer depth (in m, continuous thin line) simulated by the model are superimposed to the chlorophyll plot.

Fig. 10 .
Fig. 10.Comparison between time series of model simulated primary production rates and observations at AA1 site in the Gulf of Trieste.Model data (continuous thick line) are computed as a monthly climatological average of surface primary production rates.Observations are single field measurements collected once per month in the period 1998-2000.

Fig. 11 .
Fig. 11.Results of the sensitivity experiment on the formulation of bacterial growth efficiency (BGE).Comparison between time series of model simulated bacterial carbon production (BCP) rates obtained with the BGE parameterization proposed by Rivkin and Legendre (2001) and observations at AA1 site in the Gulf of Trieste.The circles (•) are the seasonally averaged observations in the period 1998-2000 plotted with the standard deviations.The dashed thick line (---) is the modeled depth-integrated BCP with constant BGE of 0.3, the continuous line (-) is the BCP with the BGE temperature-dependence and the squares ( ) are the relative seasonal means with the standard deviations.

Fig. 12 .
Fig. 12. Model simulated concentration of (a) dissolved organic carbon (DOC) and (b) dissolved inorganic phosphorus (DIP) as average within the critical compensation depth at the three implementation sites.

Fig. 13 .
Fig.13.Indices of ecosystem functioning and matter-transfer pathways.(a) ratio between the carbon flow due to herbivorous grazing (from autotrophs to zooplankters) and the one due to microbial grazing (from bacterioplankton to zooplankters; in semi-logarithmic scale).(b) ratio between the C-component and P-component in dissolved organic matter (DOC:DOP), indicating the "quality" of substrate for bacteria growth.The optimal ratios for phytoplankton (106:1;Redfield et al., 1963) and bacteria (45:1;Goldman et al., 1987) are marked in the plot.(c) phosphorus flux between bacteria and dissolved inorganic phosphorus (DIP).

Table 1 .
Values of selected physical parameters used in the numerical simulations

Table 2 .
Values of selected biological parameters used in the numerical simulations

Table 3 .
Number and distribution among variables of the available casts extracted from the historical data sets at the three selected sites