# Development of a 3D Coupled Physical-Biogeochemical Model for the Marseille Coastal Area (NW Mediterranean Sea): What Complexity Is Required in the Coastal Zone?

Marion Fraysse, Christel Pinazo, Vincent Martin Faure, Rosalie Fuchs, Paolo Lazzari, Patrick Raimbault, Ivane Pairaud

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marion Fraysse, Christel Pinazo, Vincent Martin Faure, Rosalie Fuchs, Paolo Lazzari, et al., Develop- ment of a 3D Coupled Physical-Biogeochemical Model for the Marseille Coastal Area (NW Mediter- ranean Sea): What Complexity Is Required in the Coastal Zone?. PLoS ONE, 2013, 8 (12), pp.e80012. 10.1371/journal.pone.0080012. [hal-01232380](https://hal.science/hal-01232380v1)

HAL Id: hal-01232380

https://hal.science/hal-01232380v1

Submitted on 17 Sep 2018

HAL is a multi-disciplinary open access archive for the deposit and dissemination of sci- entific research documents, whether they are pub- lished or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Development of a 3D Coupled Physical-Biogeochemical Model for the Marseille Coastal Area (NW Mediterranean Sea): What Complexity Is Required in the Coastal Zone?

Marion Fraysse¹,²,³*, Christel Pinazo²,³, Vincent Martin Faure²,³, Rosalie Fuchs²,³, Paolo Lazzari⁴, Patrick Raimbault²,³, Ivane Pairaud¹

¹ IFREMER, Laboratoire Environnement Ressources Provence Azur Corse, La Seyne sur Mer, France, ² Aix Marseille Université, CNRS/INSU, IRD, Mediterranean Institute of Oceanography (MIO), UM 110, Marseille, France, ³ Université de Toulon, CNRS/INSU, IRD, Mediterranean Institute of Oceanography (MIO), UM 110, La Garde, France, ⁴ Dept. of Oceanography, Istituto Nazionale di Oceanografia e di Geofisica Sperimentale (OGS), Trieste, Italy

## Abstract

Terrestrial inputs (natural and anthropogenic) from rivers, the atmosphere and physical processes strongly impact the functioning of coastal pelagic ecosystems. The objective of this study was to develop a tool for the examination of these impacts on the Marseille coastal area, which experiences inputs from the Rhone River and high rates of atmospheric deposition. Therefore, a new 3D coupled physical/biogeochemical model was developed. Two versions of the biogeochemical model were tested, one model considering only the carbon (C) and nitrogen (N) cycles and a second model that also considers the phosphorus (P) cycle. Realistic simulations were performed for a period of 5 years (2007–2011). The model accuracy assessment showed that both versions of the model were able of capturing the seasonal changes and spatial characteristics of the ecosystem. The model also reproduced upwelling events and the intrusion of Rhone River water into the Bay of Marseille well. Those processes appeared to greatly impact this coastal oligotrophic area because they induced strong increases in chlorophyll-a concentrations in the surface layer. The model with the C, N and P cycles better reproduced the chlorophyll-a concentrations at the surface than did the model without the P cycle, especially for the Rhone River water. Nevertheless, the chlorophyll-a concentrations at depth were better represented by the model without the P cycle. Therefore, the complexity of the biogeochemical model introduced errors into the model results, but it also improved model results during specific events. Finally, this study suggested that in coastal oligotrophic areas, improvements in the description and quantification of the hydrodynamics and the terrestrial inputs should be preferred over increasing the complexity of the biogeochemical model.

**Citation:** Fraysse M, Pinazo C, Faure VM, Fuchs R, Lazzari P, et al. (2013) Development of a 3D Coupled Physical-Biogeochemical Model for the Marseille Coastal Area (NW Mediterranean Sea): What Complexity Is Required in the Coastal Zone? PLoS ONE 8(12): e80012. doi:10.1371/journal.pone.0080012

**Editor:** Moncho Gomez-Gesteira, University of Vigo, Spain

**Received** May 24, 2013; **Accepted** September 27, 2013; **Published** December 4, 2013

**Copyright:** © 2013 Fraysse et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Funding:** This work was supported by the PACA region (http://www.regionpaca.fr/), IFREMER grant (www.ifremer.fr/), GIRAC (http://www.polemerpaca.com/ Environnement-et-amenagement-du-littoral/Gestion-de-l-eau-en-zone-cotiere/GIRAC) and PNEC-EC2CO MASSILIA (http://www.insu.cnrs.fr/actions-sur-projets/ ec2co) projects, sustained by the water agency (AERMC, www.earmc.fr/) and the PERSEUS European FP7 project (http://www.perseus-fp7.eu/). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing Interests:** The authors have declared that no competing interests exist.

* E-mail: marion.fraysse@univ-amu.fr

## Introduction

Coastal regions, located at the interface between oceanic and terrestrial systems, play a crucial role in earth system functioning [1]. Greater than 60% of the world's population lives less than 60 km from the sea, increasing the human pressure on these systems [2]. Understanding the fate of anthropogenic inputs from major cities and their impacts on the adjacent marine ecosystems is essential for the protection and management of coastal waters. Although it is clear that coastal systems are locally strongly impacted by human activities, it remains difficult to distinguish between climatic and anthropogenic forcing [1-3]. Modeling approaches can assist with this difficulty and are useful tools for studying such a complex coastal environment.

Marseille is the second largest city in France, and the metropolitan area of Marseille extends beyond the city limits with a population of 1 038 940 and a density 1 718 people per km² [4]. The density of contaminant-generating industries in the city of

Marseille and the quantity of sewage are highly representative of large, modern Mediterranean cities. Hydrodynamic and biogeochemical processes affect the transport and form of chemical compounds; for example, contaminant speciation often depends on suspended matter and on organic compounds (Particulate Organic Carbon (POC) and Dissolved Organic Carbon (DOC)) [5]. Marseille was thus chosen for the development of a numerical tool (a chain of models) to assess the raw inputs (from city to sea) and exports (from mid-sea to open sea) of chemical contaminants. This tool was developed based on the coupling of a hydrodynamic model [6], a sedimentary model, a biogeochemical model and a model of chemical contamination. This paper presents the coupled hydrodynamic-biogeochemical compartment.

The Marseille coastal area is located in the eastern part of the Gulf of Lions (GoL) in the western Mediterranean Sea (Figure 1). The GoL is one of the most productive areas of the Mediterranean Sea [5], even if the Mediterranean Sea remains oligotrophic. The biogeochemical functioning of the GoL is strongly impacted by
inputs from the Rhone River as it is the most significant source of
freshwater and nutrients in the Mediterranean Sea; these inputs
have a direct influence on the primary production. On an annual
basis, approximately 50% of the primary production in the GoL
can be attributed to terrigenous inputs [7,8]. The hydrodynamics
of the GoL are complex and highly variable [9] because there is
strong temporal and spatial variability in the forcing occurring at
the eastern part of the GoL [10].

The biogeochemical functioning of the Bay of Marseille (BoM)
is complex and highly driven by hydrodynamics. The hydrody-
namics of the Marseille coastal area were studied in details using a
modeling approach [6]. The primary forcing components are the
two dominant winds (north-northwesterly winds, which favor
upwelling, and southeasterly winds, which favor downwelling) and
the oligotrophic Northern Current (NC), which flows along the
continental slope toward the west [11] and occasionally intrudes
on the shelf [6,12–14]. North-northwesterly wind gusts induce
upwelling zones off the “Cote Bleue” and off “Calanques” [6,11].
During an upwelling event, cold, rich waters are advected upwards
[15], which can induce (in favorable cases) an increase in primary
production. A case study undertaken in the BoM [16] showed that
primary production tripled at a coastal station influenced by
upwelling compared with a reference offshore station. Observation
[17] and modeling studies of the Region Of Freshwater Influence
(ROFI) [18–20] show the predominant westwards direction of the
Rhone River plume. A less common orientation of the Rhone
River plume, toward the East up to 40 km from the Rhone River
mouth and offshore in the BoM, was recently observed [21]. The
presence of water from the Rhone River in the BoM was
established in a modeling study by Pairaud et al. [6]. The effects of
the eastward intrusion events of the Rhone River plume were

observed on biological production, local phytoplankton blooming
and chromophoric dissolved organic matter (CDOM) production
was measured [22]. The BoM is also impacted by urban rivers, a
Wastewater Treatment Plant (WWTP) and atmospheric deposi-
tion, but their impacts on the biogeochemical functioning of the
Bay and nutrient limitation remain to be studied.

The knowledge available on nutrient limitation is at the scale of
the Mediterranean Sea and its Western Basin. Thus, it is essential
to understand how photosynthetic production is or is not limited
by nitrogen (N) and phosphorus (P) because the cycles of key
nutrient elements, such as N and P, have been massively altered by
anthropogenic activities [23]. In the Mediterranean Sea, N or P
concentrations are generally considered to be limiting factors for
algal production [24]. Van Wambeke et al. [25] confirm that P
limitation of bacterioplankton is a generalized phenomenon in the
Mediterranean Sea. Indeed, with a NO₃:PO₄ ratio of 65–80
[26,27], the Rhone River contributes to the relative P deficit of the
Mediterranean Sea [5]. Atmospheric inputs also provide N in
excess relative to P (mean DIN/DIP ratio of 60 in the Western
Mediterranean Basin) and contribute to the P limitation [28]. In
addition, in the northwestern Mediterranean Sea, P limitation for
both phytoplankton and bacteria was suggested by Thingstad et al.
[29]. Therefore, a modeling approach could aid in the
understanding of the limitation functioning of the Marseille
coastal area.

3D physical and biogeochemical modeling approaches were
used in the GoL to study upwelling [30], nitrate fluxes between the
margin and the open sea [31], the functioning of the planktonic
ecosystem during spring and its impact on particulate organic
carbon deposition [32] and the biogeochemical functioning of
eddies [33]. Recently, Fontana et al. [34] performed assimilation

<figure><img src="image_3.png" /><figcaption>**Figure 1. Map of the Gulf of Lion and the model domain.** The RHOMA model domain (red dotted lines), the Bay of Marseille (green rectangle), the “Cote Bleue” and “Calanques” upwelling spots (blue circles), the Rhone River plume and intrusion zone (in brown) and the Northern current (NC) are represented.</figcaption></figure>

doi:10.1371/journal.pone.0080012.g001
of chlorophyll-a remote sensing data in a 3D coupled physical-biogeochemical model for the area near the mouth of the Rhone River. This study represents the first 3D, high-resolution coupled physical-biogeochemical modeling study of the BoM. Coastal zones are highly variable both hydrodynamically and biogeochemically; therefore, the modeling approach cannot be the same as used in the open sea. Franks [35] strongly suggested that the model should be designed around the scientific question being addressed. Indeed, the choice of the complexity of the biogeochemical model usually depends on the study to be conducted [36]. In our case, a model light in computational cost has many advantages for operational applications and for coupling it with other models (sediment and contamination models). Moreover, the complexity of 3D coupled models can make them particularly difficult to evaluate. Arhonditsis and Brett [37] performed a systematic analysis of 153 biological models that incorporated plankton. Only 47% of the models assessed had any performance-related accuracy assessment, and only 30% determined some measure of goodness of fit with respect to observed values. In this paper, we attempted to develop the simplest model to address our questions and then to assess its ability to reproduce both temporal and spatial variability and the primary processes of this coastal ecosystem.

The objective of this study was to develop and to assess a physical-biogeochemical coupled model to characterize the biogeochemical functioning of the Marseille coastal zone. An initial version of the coupled model was developed, and then the biogeochemical model was improved by adding the P cycle. Then, we assessed both versions of the model and whether the addition of the P cycle yields an improvement in the model simulation. In the final section, we discuss the complexity required in coastal ecosystem modeling and the contribution of model results to the understanding of the Marseille coastal area.

# Materials and Methods

## 1. Models

**1.1. The hydrodynamic model: MARS3D RHOMA.** The hydrodynamic model used for this study was the free surface, three-dimensional MARS3D model (3D hydrodynamic Model for Applications at regional Scale, IFREMER). The high resolution MARS3D-RHOMA configuration was applied and validated to the forecast of the oceanic circulation off Marseille [6], with a horizontal resolution of 200 m and 30 sigma vertical levels. The time step was fixed at 30 s. Atmospheric forcing, hydrodynamic open boundary conditions and the numerical schemes of the model were described by Pairaud et al. [6]. The model grid resolution used in this work is a downgrade of the version described in Pairaud et al. [6], with a 400-m horizontal resolution. Similar validation was performed on the 400-m version, and the processes were well reproduced.

**1.2. Coupled model.** The coupled model domain covered an area of 100 km × 48 km between the mouth of the Rhone River and Cape Sicié. The study area was discretized horizontally using a uniform mesh of 252 × 120 cells at 400-m resolution. The vertical direction was divided into 30 sigma levels (refining the resolution close to the surface and bottom).

The hydrodynamic model and the biogeochemical model were coupled online (Equation 1). The biogeochemical model calculated all variables tendency (Source-Minus-Sinks, SMS) for each grid point every 20 minutes, whereas the physical model performed the advection-diffusion of the biogeochemical concentrations (C) with a 30-second time-step. The analytical formulations of the advection and diffusion terms are the same as those for the

temperature and salinity in the physical model.

$$ \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} + w \frac{\partial C}{\partial z} = \frac{\partial}{\partial x} \left( K_x \frac{\partial C}{\partial x} \right) + \frac{\partial}{\partial y} \left( K_y \frac{\partial C}{\partial y} \right) + \frac{\partial}{\partial z} \left( K_z \frac{\partial C}{\partial z} \right) + SMS \quad (1) $$

The coupling was introduced as the extinction light coefficient calculated by the biogeochemical model used in the function of the chlorophyll-a concentration, which induced a retro-action on the temperature calculation of the physical extinction light coefficient; thus, it was truly a coupled online model. Total radiance was read from the meteorological file (MM5) by the hydrodynamic model.

**1.3. The biogeochemical model: ECO3M MASSILIA.** The biogeochemical model was implemented using the Eco3M (Ecological Mechanistic and Modular Modeling) modeling platform [38,39]. A new biogeochemical model (ECO3M-MASSILIA) was developed for this study. The model structure used is primarily based on the pelagic plankton ecosystem model published by Faure et al. [40,41] but without the P cycle.

The biogeochemical model was split into 5 compartments (phytoplankton, heterotrophic bacteria, dissolved and particulate organic matter and dissolved inorganic matter). This model allowed for a variable intracellular content of phytoplankton and bacteria. Another particularity was that the zooplankton were not represented as a state variable of the model; rather, their physiological fluxes (e.g., grazing or excretion) were considered as explicit functions. The closure formulation of the model was based on the assumption that all of the matter grazed by the zooplankton and the higher trophic levels was returned to one of the pools of organic or inorganic matter (dissolved and particulate organic matter and dissolved inorganic matter).

In this paper, we present and discuss only the modifications applied to the initial version of the model. First, the ECO3M-MASILLIA-noP was developed with new parameterization adapted for the Mediterranean Sea and the addition of grazing limitation by temperature. Then, the ECO3M-MASILLIA-noP was modified with the addition of the P cycle (ECO3M-MASILLIA-P). Thus, the pelagic ecosystem was summarized in 12 or 17 state variables (Table 1), and the cycles of C and N were described by the ECO3M-MASSILIA-noP, with the addition of P for the ECO3M-MASSILIA-P (Figure 2). We refer the reader to Annex S1 for the detailed equations of the biogeochemical model.

The parameterization of the biogeochemical model was adapted for the Marseille coastal area (Annex S2). Previous studies of the phytoplankton community in the Marseille - Rhone area concluded that diatoms dominated the phytoplankton community [42,43]. When possible, phytoplankton parameters were chosen to represent diatoms. For the purpose of parameter refinement, we implemented zero-dimensional models, which allowed us to perform many simulations with different parameters sets.

A coastal study in the northwestern Mediterranean Sea demonstrated that the grazing rate had strong seasonality: the grazing is lower in winter and higher in summer [44]. Another study reported that the development of the zooplankton was impacted by the water temperature [45]. This seasonality was introduced into the model as a function linking the grazing rate to the temperature (T) (Equation 2). The limitation function f(T) [46] was applied to the grazing rate.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The pelagic biogeochemical model. The ECO3M-MASSILIA-noP (only white-colored variables) and ECO3M-MASSILIA-P (white- and purple-colored variables) versions of the biogeochemical model are represented.
doi:10.1371/journal.pone.0080012.g002</figcaption></figure>

$$f(T) = 0.0498 \cdot \exp(0.13 \cdot T) \quad (2)$$

The literature on the limitation by P [25,29] and the potential deficiency in P of diverse inputs into the study area [26–28] highlight the importance of considering this element. For the addition of the P cycle in ECO3M-MASSILIA-P, five supplementary states variables were added (Figure 2): dissolved inorganic phosphorus (PO₄), labile dissolved organic phosphorus (LDOP), detrital particulate organic phosphorus (DPOP), phosphorus phytoplankton biomass (Pb) and phosphorus bacterial biomass (Pba). For coherence, the representation of the P cycle mimicked that of the N cycle. The structural modifications are described hereafter for the P equations that differed from those of the N. The full model equations are detailed in Annex S1.

Concerning the phytoplankton compartment, the phytoplankton biomass in phosphorus (Pb) was introduced as a state variable. Pb depends on the PO₄ uptake and zooplankton grazing (Equation 3).

$$\frac{\partial Pb}{\partial t} = +uptake - grazing \quad (3)$$

The equations of the phytoplankton biomass in carbon and chlorophyll were modified. The carbon biomass equation was modified through the term for the limitation of the phytoplankton maximal growth rate (P<sub>m</sub><sup>C</sup>) by nutrient (Q*) (Equation 4). The co-limitation of phytoplankton by nutrients was introduced into the equation for phytoplankton carbon biomass as an independent nutrient co-limitation of type I [47]. The co-limitation between N and P was calculated using the Liebig's Law of the minimum [48,49]. The phytoplankton limitation by P was calculated in the same manner as by N: with the intracellular quota (Q) and the

<table><caption>**Table 1.** Biogeochemical model variables (ECO3M-MASILLIA-P version).</caption><thead><tr><th>Variable</th><th>Definition</th></tr></thead><tbody><tr><td>1</td><td>C<sub>B</sub> Phytoplankton Carbon</td></tr><tr><td>2</td><td>N<sub>B</sub> Phytoplankton Nitrogen</td></tr><tr><td>3</td><td>Chl.a Phytoplankton Chlorophyll-a</td></tr><tr><td>4</td><td>CBA Bacterial Carbon</td></tr><tr><td>5</td><td>NBA Bacterial Nitrogen</td></tr><tr><td>6</td><td>DPOP Detrital Particulate Organic Carbon</td></tr><tr><td>7</td><td>DPON Detrital Particulate Organic Nitrogen</td></tr><tr><td>8</td><td>LDOC Labile Dissolved Organic Carbon</td></tr><tr><td>9</td><td>LDON Labile Dissolved Organic Nitrogen</td></tr><tr><td>10</td><td>NO3 Nitrates</td></tr><tr><td>11</td><td>NH4 Ammonium</td></tr><tr><td>12</td><td>O Oxygen</td></tr><tr><td>13</td><td>PB Phytoplankton Phosphorus</td></tr><tr><td>14</td><td>PBA Bacterial Phosphorus</td></tr><tr><td>15</td><td>DPOP Detrital Particulate Organic Phosphorus</td></tr><tr><td>16</td><td>LDOP Labile Dissolved Organic Phosphorus</td></tr><tr><td>17</td><td>PO4 Phosphorus</td></tr></tbody></table>

Only variables from 1 to 12 for ECO3M-MASILLIA-noP version. Units in µmol.L⁻¹ except Chlophyll-a which is in µg.L⁻¹.
doi:10.1371/journal.pone.0080012.t001

Geider formulation [50]. Consequently, phytoplankton production is controlled by the most limiting internal ratio (Q*).

$$Q^* = \min \left( \left( \frac{Q_C^N - \min Q_C^N}{\max Q_C^N - \min Q_C^N} \right), \left( \frac{Q_C^P - \min Q_C^P}{\max Q_C^P - \min Q_C^P} \right) \right) \quad (4)$$

The phytoplankton chlorophyll-a concentration was a diagnostic variable related to the phytoplankton biomass by carbon (C<sub>B</sub>) and the variable ratio Chla:C (Q<sub>C</sub><sup>Chla</sup>) (Equation 5). Equation 6 describes the Q<sub>C</sub><sup>Chla</sup> of phytoplankton depending on the phytoplankton N-to-C ratio and the limitation by the most limiting internal ratio (Q*). As Q* decreases toward zero, Q<sub>C</sub><sup>Chla</sup> also decreases, as does the limitation on Chla. In Smith and Tett [51] and Faure et al. [52], the variable Q* only represented the internal N:C ratio. To also consider the limitation of chlorophyll-a production when phytoplankton was limited by P, we modified the equation for Q* (Equation 4).

$$[Chla] = Q_C^{Chla} \cdot C_B \quad (5)$$

$$\text{With } Q_C^{Chla} = Q_C^N \cdot (\max Q_N^{Chla} \cdot Q^* + \min Q_N^{Chla} (1 - Q^*)) \quad (6)$$

As for phytoplankton, heterotrophic bacteria were described in terms of carbon, nitrogen and phosphorus. To include the P cycle in the microbial loop, the bacterial biomass in phosphorus (Pba) was added as a new state variable. Pba varied with the uptake and grazing of organic matter (DPOP and LDOP) and inorganic matter (PO₄) (Equation 7).
$$ \frac{\partial P_{ba}}{\partial t} = + \text{uptake} (\text{DPOP}, \text{LDOP}, \text{PO}_4) - \text{grazing} \quad (7) $$

We simulated the temporal evolution of heterotrophic bacteria according to the cell quota theory. This concept [53] considers that growth rate (μ<sup>BA</sup><sub>max</sub>) is a function of limiting nutrients within the cell. A new system of co-limitation by C, N and P was needed to control the bacterial production (BP) (Equation 8). We choose to represent the co-limitation (f) of the bacterial production by the Liebig's law (Equation 9), similarly to the representation of phytoplankton production. However, even if Liebig's law was used for the co-limitation of the phytoplankton, the limitation of the maximum growth rate (μ<sup>BA</sup><sub>max</sub>) of bacteria was calculated differently. Q<sup>C</sup><sub>BA</sub>, Q<sup>N</sup><sub>BA</sub> and Q<sup>P</sup><sub>BA</sub> represents the carbon, nitrogen and phosphorus cell quota of bacteria (see Annex S1 for more details). As the ratio (min Q<sup>X</sup><sub>BA</sub> : Q<sup>X</sup><sub>BA</sub>) increased, the limitation of the bacteria by the element X also increased. The most restrictive element corresponds to the higher ratio (min Q<sup>X</sup><sub>BA</sub> : Q<sup>X</sup><sub>BA</sub>).

$$ BP = \mu_{\max}^{BA} (1 - f) \cdot NBA \quad (8) $$

$$ \text{With } f = \max \left( \frac{\min Q_{BA}^{C}, \min Q_{BA}^{N}, \min Q_{BA}^{P}}{Q_{BA}^{C}, Q_{BA}^{N}, Q_{BA}^{P}} \right) \quad (9) $$

**1.4. Modeling strategy.** To determine the spin-up period of the model, tests were conducted for the year 2007 by replacing winter initial conditions with summer initial conditions. We considered the maximum spin-up period to be completed when all state variables were equal between the two runs. The period lasted 90 days; therefore, the model results can be studied from April 2007 onwards (Annex S3).

The hydrodynamic open boundary conditions and initial conditions of the MARS3D-RHOMA configuration are described in Pairaud et al. [6]. The biogeochemical open boundary conditions and initial conditions were built from a coupled model applied to a larger area (GoL) at a horizontal resolution of 1.2 km and 30 vertical sigma levels. The hydrodynamic model (MARS 3D-GOL configuration) was based on the MARS3D-MENOR configuration validated by Nicolle et al. [54]. The MARS 3D-GOL configuration was coupled online with the biogeochemical ECO3M-MASSILIA-P [55]. The biogeochemical open boundary conditions of this larger model were provided by the OPATM-BFM pre- and operational runs 2007–2011 performed in the framework of the Mersea and MyOcean projects [56].

Atmospheric dry deposition (dry and wet deposition) was sampled weekly by the national MOOSE program (Mediterranean Oceanic Observing System on Environment) and the “Service d'Observation” of the Mediterranean Institute of Oceanography (S.O. MIO). Inputs from rain (wet deposition) were measured for each rain event. The device (collecteur MTX-Italia) was installed on an island located in the bay off Marseille. Atmospheric samples were analyzed for soluble components (nitrate (NO<sub>3</sub>), ammonium (NH<sub>4</sub>), phosphate (PO<sub>4</sub>), dissolved organic carbon (DOC), dissolved organic nitrogen (DON) and dissolved organic phosphorus (DOP)) and insoluble particulate matter (particulate organic carbon (POC), particulate organic nitrogen (PON) and particulate phosphorus (POP)). Wet atmospheric concentrations were applied to the rainfall represented by the atmospheric model MM5. Dry

deposition was applied as a mean flux between each pair of sample dates.

Between 2007 and 2011, daily averaged Rhone River discharges at the Beaucaire station were available and provided by the “Compagnie Nationale du Rhone”. The “Grand Rhone” located in the modeling domain represents only 90% of total Rhone River discharge [57,58]. Marseille has 4 main Urban Rivers (Aygalade, Belvedère-Figuière, Huveaune-Jarret and Bonneveine) and a Wastewater Treatment Plant (WWTP), which flows in at Cortiou (Annex S4). The Marseille Urban Rivers discharges were available with a time step of 6 minutes, and the daily discharge from the Marseille WWTP was available (Data provided by the DEA-MPM (Direction de l'Eau et de l'Assainissement-Marseille Provence Métropole)). The Berre Lagoon, which is a shallow semi-confined ecosystem, is connected to the Mediterranean Sea via the Caronte channel. The latter was represented in the coupled model as a river with a constant discharge of 20 m<sup>3</sup>·s<sup>-1</sup>. River and WWTP concentrations.

Biogeochemical concentrations in the Rhone River (NO<sub>3</sub>, NH<sub>4</sub>, PO<sub>4</sub>, POC, PON, POP, DOC, DON and DOP) were measured daily and provided by the national MOOSE program and the SO MIO. As our model considered only the labile fraction of organic matter, a percentage representing the labile fraction of organic matter was applied to the daily Rhone River concentrations. However, studies of the lability of Rhone River organic matter were rare. The labile fraction of particulate organic carbon (DPOC) was considered to be 18% of POC [59]. Déliat [60] estimated that approximately 20% of the Rhone River DOC was biodegradable; thus, we computed LDOC as equal to 20% of DOC. The labile fraction of DON (LDON) was also estimated to be 20% of DON [60]. Experiments in Loch Creran (Scotland) showed that bioavailable DOP (BDOP) accounted for 88±8% of DOP (average±SD) [61]. The same percentages were applied to the particulate organic matter concentrations to obtain DPON and DPOP.

The concentrations in the urban rivers were not always available from in-situ data for all of the variables and hence were derived from empirical relationships (Table 2). The Rhone River concentrations (described above) were applied to the concentrations of nutrients in other Marseille Urban Rivers. Organic matter concentrations in the Marseille Urban Rivers and from the WWTP were deduced from suspended particulate matter (SPM) concentrations [62] with constant ratios. The same hypotheses on the lability of the Rhone River organic matter were applied to the urban river concentrations. The concentrations of NH<sub>4</sub> and NO<sub>3</sub> from the Marseille WWTP were available daily (Data provided by the DEA-MPM (Direction de l'Eau et de l'Assainissement-Marseille Provence Métropole)). The PO<sub>4</sub> concentration was fixed to 13.4 µmol (Faure, comm. Pers.). Marine phytoplankton and bacteria species were considered to be absent from all of the rivers and WWTP inputs. The Caronte concentrations were fixed at constant values (Table 2) obtained by averaging the measured concentrations [63].

## 2. Observational data sets

**2.1. In-situ data.** To validate the coupled model, we used a long time series of hydro-biogeochemical data collected twice monthly at the Somlit station (43°14.30′N; 5°17.30′E) located in the BoM (Annex S4). High vertical resolution profiles of temperature, salinity and oxygen were obtained between 0 and 55 m using a conductivity temperature-depth-oxygen profiler (CTDO, Seabird 19+). Water samples were collected at 3 depths with hydrological Niskin bottles for the determination of inorganic and organic nutrients concentrations. Samples for nitrate, nitrite
<table><caption>**Table 2.** Concentrations of River inputs (µmol.L⁻¹).</caption><thead><tr><th>River</th><th>NO<sub>3</sub></th><th>NH<sub>4</sub></th><th>PO<sub>4</sub></th><th>LDOC</th><th>LDON</th><th>LDOP</th><th>LPOC</th><th>LPON</th><th>LPOP</th></tr></thead><tbody><tr><td>Rhone</td><td>data</td><td>data</td><td>data</td><td>0.2*DOC</td><td>0.2*DON</td><td>0.88*DOP</td><td>0.18*POC</td><td>0.2*POC</td><td>0.88*POP</td></tr><tr><td>Caronte</td><td>0.75</td><td>2.72</td><td>0.4</td><td>45</td><td>3.73</td><td>0.11</td><td>1.8</td><td>1.36</td><td>0.078</td></tr><tr><td>WWTP</td><td>data</td><td>data</td><td>13.4</td><td>135</td><td>17</td><td>3</td><td>8403</td><td>864</td><td>372</td></tr><tr><td>Bonneveinee</td><td>RR</td><td>RR</td><td>RR</td><td>52.4</td><td>6.62</td><td>1.17</td><td>667</td><td>68.60</td><td>29.6</td></tr><tr><td>Huveaune</td><td>RR</td><td>RR</td><td>RR</td><td>38.5</td><td>4.86</td><td>0.86</td><td>869</td><td>89.4</td><td>38.57</td></tr><tr><td>Emiss 1</td><td>RR</td><td>RR</td><td>RR</td><td>553</td><td>69.9</td><td>12.37</td><td>2430</td><td>250</td><td>107.8</td></tr><tr><td>Emiss 2</td><td>RR</td><td>RR</td><td>RR</td><td>14.4</td><td>1.82</td><td>0.32</td><td>f(RD)</td><td>f(RD)</td><td>f(RD)</td></tr><tr><td>Aygalade</td><td>RR</td><td>RR</td><td>RR</td><td>52.4</td><td>6.62</td><td>1.17</td><td>f(RD)</td><td>f(RD)</td><td>f(RD)</td></tr><tr><td>Belvedere</td><td>RR</td><td>RR</td><td>RR</td><td>52.4</td><td>6.62</td><td>1.17</td><td>f(RD)</td><td>f(RD)</td><td>f(RD)</td></tr></tbody></table>

RR: Rhone River daily measured concentrations; f(RD) : concentrations linked to urban rivers discharges.
doi:10.1371/journal.pone.0080012.t002

and phosphate were collected into 60 ml polyethylene nutrients flasks and were frozen at -20°C until analysis at laboratory according to Aminot and Kerouel [64]. Samples for acid silicic were collected in 60 ml polyethylene flasks and stored at 5°C until analysis at laboratory according to [64]. Samples for ammonium determination were collected in triplicate in 60 ml polycarbonate tubes. The reagent was immediately added to the tubes and ammonium level was determined by fluorometry according to Holmes et al. [65]. Particulate organic carbon (POC) and particulate organic nitrogen (PON) in suspended matter collected on Whatman GF/F glass micro-fibre filters pre-combusted for 4 h at 450°C, were determined by using the high combustion method (1000°C) on a CN Integra mass spectrometer [66]. Chlorophyll concentrations were estimated by fluorometry [67] on suspended matter collected on Whatman GF/F filter. The data were provided by the SOMLIT network (Service d'Observation en Milieu Littoral, http://somlit.epoc.u-bordeaux1.fr).

**2.2. Satellite remote sensing data.** We also compared the model-predicted fields of surface chlorophyll-a concentrations with the corresponding remote sensing-derived concentrations. The MODIS (Moderate Resolution Imaging Spectroradiometer) and MERIS (MEdium Resolution Imaging Spectrometer) ocean color sensors have a spatial resolution of approximately 1 km, which is more course than the spatial resolution of our model (400 m). For the comparisons, we used the remotely sensed chlorophyll-a concentration processed with the algorithm OC5 [68,69] from the IFREMER (Institut Français de Recherche pour l'Exploitation de la Mer) database. The OC5 method [68] is empirical and derived from the OC4 algorithm of NASA (or OC3M-547 for MODIS and OC4E for MERIS). This method gives results similar to OC4 in open waters but provides more realistic values over the continental shelf [70]. The quality of the remote sensing data was evaluated in Annex S5. The remote sensing data contained observational errors, which was also the case for the in-situ data. Nevertheless, the large number of existing data, their availability throughout the year and their spatial coverage make remote sensing data essential to the consideration of spatial patterns and gradients. Thus, the comparison with the model results had to consider the observational error.

### 3. Statistical comparison with observations

To assess the accuracy of the 3D coupled online model, we used statistics highlighting correspondences between the model results (M) and in-situ and remote sensing observations (O). The study of model performance often requires multiple stages of analysis, the

evaluation of the ability of the model to reproduce instantaneous station values (and trends) and the ability of the model to recreate spatial characteristics (and trends) of processes [71]. Stow et al. [72] explained that model start to have skill when the observational and predictive uncertainty halos overlap; in the ideal case, the halos overlap completely. Thus, accuracy assessment requires a set of quantitative metrics and procedures for comparing model output with observational data in an appropriate manner for the particular application. Here, we used different statistical indicators [72,73]:

- The percentage model bias (PB, model error normalized by the observations) measures whether the model underestimated or overestimated the observations [74]. The closer to zero the value is, the better the model.

- The cost function (CF) gives a non-dimensional value that is indicative of the “goodness of fit” between two sets of values [73]. Radach and Moll [75] proposed an interpretation of the values of the cost functions adapted from the OSPAR Commission [76] rating. Model results are classified as very good (cost function between 0 and 1), good (1 and 2), reasonable (2 and 3) and poor (higher than 3) [77].

- The correlation coefficient (R) is a measure of the strength and direction of the linear relationship between two sets of values [74].

- The average absolute error (AAE) quantifies the magnitude rather than the direction of each discrepancy (Equation 10). The closer to zero the value is, the better the model.

$$AAE = \frac{\sum_{i=1}^{n} |O_i - M_i|}{n} \qquad (10)$$

- The Root Mean Square Deviation (RMSD) is also a measure of the difference between values predicted by the model and the values actually observed from the environment [74]. The RMSD has the same units as the quantity being evaluated.

To compare the variables, we used signed, normalized, unbiased RMSD (snuRMSD) to quantify the magnitude direction of each discrepancy and the normalized bias to quantify the direction of each discrepancy [78]. Those indicators are used to construct target diagrams [78].
Model developments and parameterization were performed for the years 2007 and 2008, and then the observations available for the years 2009, 2010 and 2011 were used only to assess the accuracy of both versions of the model.

## Results

We assessed the ability of the coupled model to reproduce the main characteristics of the Marseille coastal area for the two versions of the biogeochemical model, with and without the P cycle.

An evaluation of a 3D ecosystem model starts with the hydrodynamics. Pairaud et al. [6] performed the hydrodynamics study for the years 2007 and 2008. They computed model/ observation statistics and showed the ability of the model to capture monthly to seasonal variability in the thermal structure and to reproduce the observed features over the shelf, such as the warming or cooling of the sea due to upwelling and downwelling events and the punctual extension of the Rhone River plume.

Realistic 3D simulations were performed using the coupled model for 2007 through 2011. First, we present the model-observation comparisons, where we examined the temporal dynamics in the surface and bottom layers at the Somlit station and computed statistics. Then, we present spatial maps to assess the representation of the spatial gradient and their variability during the year. Finally, we focused on short time scale shelf processes and their impacts on the biogeochemical functioning of the pelagic coastal ecosystem.

# 1. Observations/model comparisons at the Somlit station

**1.1. Temporal dynamics.** To estimate the reliability of the model in qualitative terms, we described the agreement between simulated and observed monthly means, in-situ measurements of temperature, chlorophyll-a, NO₃, particulate organic matter and PO₄ (Figure 3-4). At the Somlit coastal station (Annex S4), both model versions (with and without the P cycle) captured the annual cycle, with strong seasonal oscillations in temperature, NO₃ and chlorophyll-a. The error bars showed the variability at a smaller time scale than the month of the simulation. If the winter temperatures were estimated well, the maximum summer temperatures seemed to be significantly underestimated. Nitrate concentrations and seasonal evolution were well quantified by the model. Higher chlorophyll-a variability was observed over shorter time scales in spring and summer than in winter and autumn (Figure 3). Both versions of the model reproduced the winter period well when the vertical mixing induced high NO₃ concentrations and low temperature in the surface layer. The spring was characterized by high chlorophyll-a concentrations (spring bloom) associated with a decrease in NO₃. Peaks in the chlorophyll-a concentrations during this season were slightly overestimated by both model versions (Figure 3). During summer, nutrient depletion was reproduced well, while the chlorophyll-a concentration was overestimated. Some isolated values of nitrate and chlorophyll-a were captured by the model.

Both versions of the model demonstrated their abilities to capture inter-annual variability, with a higher spring peak in chlorophyll-a (Figure 3), POC and PON (Figure 4) in 2010 than in 2011. In 2009, the chlorophyll-a concentration increased in November and December, which was not the case for 2007, and the model reproduced this inter-annual variability well.

The model suggests a clear and regular seasonal variation for phosphate, similar to those for nitrate, which is not the pattern observed in the field by bi-weekly sampling. The concentrations in PO₄ remained very low (lower than 0.2 µmol.L⁻¹) throughout the

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Surface model/data time series comparisons for temperature, chlorophyll-a and NO₃. Time series of the monthly mean temperature (°C)(A), chlorophyll-a concentration (µg.L⁻¹) (B) and NO₃ concentration (µmol.L⁻¹) (C) from the model with the P cycle (grey) and the model without P (black), with one standard deviation (error bar), compared with the in-situ data (blue) at the Somlit station.
doi:10.1371/journal.pone.0080012.g003</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Surface model/data time series comparisons for POC, PON and PO₄. Time series of the simulated monthly mean POC (µmol.L⁻¹) (A), PON (µmol.L⁻¹) (B) and PO₄ (µmol.L⁻¹) (C) concentrations from the model with the P cycle (grey) and the model without P (black), with one standard deviation (error bar), compared with the in-situ data (blue) at the Somlit station.
doi:10.1371/journal.pone.0080012.g004</figcaption></figure>
year, which made it difficult to distinguish a seasonal signal in the in-situ PO₄ data. The two model versions had rather similar results over the five years. However, there were a few discrepancies between them in chlorophyll-a, POC and PON, particularly in summer and spring when the model version without the P cycle overestimated the chlorophyll-a concentration more than the model version with the P cycle (Figure 3, Figure 4). Thus, the addition of the P cycle improved the results during spring and summer.

**1.2. Statistical model/observations comparisons.** This qualitative comparison was completed through the calculation of statistical indicators between the Somlit time series and the results of the simulations (only for the years 2009 to 2011). The surface and bottom results comparisons are presented in Table 3 and Table 4, respectively. Although the Somlit station is a coastal station near Marseille, the low mean concentrations of NO₃, chlorophyll-a and PO₄ (both in model and in observations) indicated that oligotrophic conditions were predominant. The standard deviation of the in-situ data was frequently equal to or very close to the mean concentration at the surface for NO₃ (0.68±0.64 µmol.L⁻¹), PO₄ (0.03±0.03 µmol.L⁻¹) and chlorophyll-a (0.4±0.32 µg.L⁻¹), which highlights that this oligotrophic ecosystem had a strong variability, possibly due to pulsed nutrients. This variability made the model/observations comparisons difficult, so that different statistical indicators were used.

First, we computed the cost function. According to the OSPAR classification [76], our results were very good or good except for bottom chlorophyll-a, POC, PON and PO₄. Then, we computed the average absolute error (AAE, Equation 10) and the RMSD, which measure the average magnitude of the errors in a dataset, without considering their direction. AAE and RMSD were expressed in the same units as the variable. In the surface layer, the addition of the P cycle induced a decrease or stagnation in the AAE for all of the variables, except for the POC, for which the increase in the AAE was not significant. The model version with the P cycle performed better at representing the chlorophyll-a concentration, with an AAE of 0.33 µmol.L⁻¹ versus 0.38 µmol.L⁻¹ for the model version without the P cycle. In the bottom layer, the AAEs were higher than in the surface layer, except for the NH₄, salinity and O₂. The addition of the P cycle induced slight variations or stagnations in the AAE for all the variables. The model with the P cycle performed slightly worse for

the bottom chlorophyll-a concentrations, with an AAE of 0.52 versus 0.49 µmol.L⁻¹ for the model without the P cycle.

Target diagrams [78] provide summary and visual information regarding the bias and the error between a model and observations. In order to compare variables with each other normalized statistics were used as described in Jolliff et al. [78]. In the Target diagram presented (Figure 5), the Y-axis corresponds to the normalized bias (Bias*) and the X-axis corresponds to the normalized unbiased RMSD (uRMSD*). Variables located in the upper part of the diagram (Y>0) were overestimated by the model (chlorophyll-a, PON and PO₄). For the X-axis, the model standard deviation is larger than the standard deviation for the observations when variables are located in the right part of the target (X>0). The distance between any point and the origin is then the value of the total RMSD (see Jolliff et al. [78] for more details). As shown by the time series and previous comparisons, the discrepancies between the two versions of the biogeochemical model are small both at the surface (Figure 5 A–B) and at the bottom (Figure 5 C–D). The main differences were in the chlorophyll-a concentrations, for which the model version with the P cycle was more efficient in the surface layer with a lower normalized bias and a smaller uRMSD* than the model version without the P cycle. Both model versions overestimated the amplitude of the variations in the state variables except for TEMP, NO₃ and NH₄. The target diagrams also demonstrate that models were less efficient in representing the bottom layer, with an increased bias for PON and POC.

To summarize, the timing and magnitude of the surface chlorophyll-a and NO₃ concentrations were generally well matched. The shortcomings of the simulation were an overestimation of the chlorophyll-a concentrations in spring and a slight overestimation in summer. The model also had more difficulty in simulating the bottom than the surface at the Somlit station. However, using the in-situ data to evaluate the ability of the model as discussed above refers to data from a single station sampled every fortnight and with undoubted local bias and particularity. We thus used remote sensing data for the validation of the representation of spatial processes.

## 2. Evaluation of spatial processes representation

Ocean color observations derived from remote sensing data were used to assess the abilities of the model to capture spatial gradients and seasonality over the whole study area.

<table><caption>**Table 3.** Surface model and data comparison at the Somlit station from 2009 to 2011.</caption><thead><tr><th></th><th>TEMP</th><th>SAL</th><th colspan="2">CHL</th><th colspan="2">NO<sub>3</sub></th><th colspan="2">NH<sub>4</sub></th><th colspan="2">O<sub>2</sub></th><th colspan="2">POC</th><th colspan="2">PON</th><th colspan="2">PO4</th></tr><tr><th></th><th></th><th></th><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th></tr></thead><tbody><tr><td>n</td><td>78</td><td>78</td><td>78</td><td>78</td><td>73</td><td>73</td><td>77</td><td>77</td><td>74</td><td>74</td><td>73</td><td>73</td><td>73</td><td>73</td><td>64</td></tr><tr><td>Mean (in situ)</td><td>17.15</td><td>37.97</td><td>0.40</td><td>0.40</td><td>0.68</td><td>0.68</td><td>0.31</td><td>0.31</td><td>241.12</td><td>241.12</td><td>5.13</td><td>5.13</td><td>0.79</td><td>0.79</td><td>0.03</td></tr><tr><td>Mean (model)</td><td>16.59</td><td>37.99</td><td>0.52</td><td>0.58</td><td>0.58</td><td>0.56</td><td>0.13</td><td>0.13</td><td>231.55</td><td>233.46</td><td>3.59</td><td>4.29</td><td>0.88</td><td>0.93</td><td>0.04</td></tr><tr><td>Std (in situ)</td><td>3.63</td><td>0.17</td><td>0.32</td><td>0.32</td><td>0.64</td><td>0.64</td><td>0.68</td><td>0.68</td><td>14.19</td><td>14.19</td><td>2.09</td><td>2.09</td><td>0.35</td><td>0.35</td><td>0.03</td></tr><tr><td>Std (model)</td><td>2.68</td><td>0.36</td><td>0.52</td><td>0.64</td><td>0.57</td><td>0.59</td><td>0.20</td><td>0.26</td><td>18.66</td><td>20.89</td><td>2.85</td><td>3.71</td><td>0.59</td><td>0.67</td><td>0.05</td></tr><tr><td>CF</td><td>0.27</td><td>1.24</td><td>1.01</td><td>1.16</td><td>0.80</td><td>0.81</td><td>0.40</td><td>0.43</td><td>1.30</td><td>1.37</td><td>1.20</td><td>1.19</td><td>1.17</td><td>1.19</td><td>1.39</td></tr><tr><td>Bias (%)</td><td>3.21</td><td>-0.05</td><td>-28.75</td><td>-44.41</td><td>14.17</td><td>17.05</td><td>59.95</td><td>59.65</td><td>3.97</td><td>3.18</td><td>30.03</td><td>16.30</td><td>-10.96</td><td>-16.84</td><td>-15.35</td></tr><tr><td>AAE</td><td>0.99</td><td>0.21</td><td>0.33</td><td>0.38</td><td>0.51</td><td>0.51</td><td>0.27</td><td>0.29</td><td>18.39</td><td>19.47</td><td>2.49</td><td>2.48</td><td>0.41</td><td>0.41</td><td>0.04</td></tr><tr><td>RMSD</td><td>1.37</td><td>0.33</td><td>0.52</td><td>0.63</td><td>0.7</td><td>0.7</td><td>0.73</td><td>0.75</td><td>22.56</td><td>23.94</td><td>3.17</td><td>3.38</td><td>0.55</td><td>0.6</td><td>0.05</td></tr><tr><td>R</td><td>0.96</td><td>0.44</td><td>0.34</td><td>0.34</td><td>0.34</td><td>0.36</td><td>0.00</td><td>-0.01</td><td>0.24</td><td>0.20</td><td>0.40</td><td>0.47</td><td>0.41</td><td>0.46</td><td>0.10</td></tr></tbody></table>

n: the number of in-situ data available for comparison, std: the standard deviation.

doi:10.1371/journal.pone.0080012.t003
<table><caption>**Table 4.** Bottom model and data comparison at the Somlit station from 2009 to 2011.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">TEMP</th><th rowspan="2">SAL</th><th colspan="2">CHL</th><th colspan="2">NO<sub>3</sub></th><th colspan="2">NH<sub>4</sub></th><th colspan="2">O<sub>2</sub></th><th colspan="2">POC</th><th colspan="2">PON</th><th colspan="2">PO4</th></tr><tr><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th><th>No P</th><th>P</th><th>No P</th></tr></thead><tbody><tr><td>n</td><td>74</td><td>74</td><td>74</td><td>74</td><td>71</td><td>71</td><td>73</td><td>73</td><td>67</td><td>67</td><td>69</td><td>69</td><td>69</td><td>69</td><td>63</td></tr><tr><td>Mean (in situ)</td><td>14.56</td><td>38.04</td><td>0.36</td><td>0.36</td><td>1.03</td><td>1.03</td><td>0.14</td><td>0.14</td><td>243.53</td><td>243.53</td><td>4.16</td><td>4.16</td><td>0.66</td><td>0.66</td><td>0.03</td></tr><tr><td>Mean (model)</td><td>14.78</td><td>38.16</td><td>0.81</td><td>0.77</td><td>0.5</td><td>0.52</td><td>0.11</td><td>0.11</td><td>233.38</td><td>233.5</td><td>6.03</td><td>6.3</td><td>1.5</td><td>1.47</td><td>0.09</td></tr><tr><td>Std (in situ)</td><td>1.75</td><td>0.09</td><td>0.22</td><td>0.22</td><td>0.89</td><td>0.89</td><td>0.11</td><td>0.11</td><td>12.38</td><td>12.38</td><td>1.6</td><td>1.6</td><td>0.28</td><td>0.28</td><td>0.03</td></tr><tr><td>Std (model)</td><td>1.25</td><td>0.11</td><td>0.56</td><td>0.52</td><td>0.58</td><td>0.58</td><td>0.15</td><td>0.15</td><td>17.41</td><td>17.92</td><td>3.3</td><td>3.55</td><td>0.66</td><td>0.65</td><td>0.06</td></tr><tr><td>CF</td><td>0.42</td><td>1.75</td><td>2.39</td><td>2.23</td><td>0.88</td><td>0.87</td><td>1.08</td><td>1.08</td><td>1.27</td><td>1.29</td><td>1.89</td><td>2.03</td><td>3.21</td><td>3.09</td><td>2.32</td></tr><tr><td>Bias (%)</td><td>-1.55</td><td>-0.32</td><td>-122.52</td><td>-111.72</td><td>51.96</td><td>50.23</td><td>20.15</td><td>16.66</td><td>4.17</td><td>4.12</td><td>-44.98</td><td>-51.49</td><td>-128.24</td><td>-123.06</td><td>-164.91</td></tr><tr><td>AAE</td><td>0.74</td><td>0.15</td><td>0.52</td><td>0.49</td><td>0.78</td><td>0.77</td><td>0.12</td><td>0.12</td><td>15.68</td><td>15.94</td><td>3.02</td><td>3.25</td><td>0.89</td><td>0.85</td><td>0.07</td></tr><tr><td>RMSD</td><td>0.95</td><td>0.17</td><td>0.71</td><td>0.65</td><td>1.09</td><td>1.08</td><td>0.18</td><td>0.18</td><td>19.49</td><td>19.76</td><td>3.84</td><td>4.15</td><td>1.06</td><td>1.03</td><td>0.08</td></tr><tr><td>R</td><td>0.86</td><td>0.29</td><td>0.25</td><td>0.26</td><td>0.21</td><td>0.22</td><td>0.06</td><td>0.05</td><td>0.41</td><td>0.41</td><td>0.19</td><td>0.2</td><td>0.24</td><td>0.26</td><td>0.04</td></tr></tbody></table>

doi:10.1371/journal.pone.0080012.t004

**2.1. Climatology of surface chlorophyll-a patterns.** The maps of time-averaged surface chlorophyll-a concentrations permit the evaluation of the ability of the model to reproduce horizontal patterns of chlorophyll-a. The outputs from both versions of the model were compared with observed maps of surface chlorophyll-a concentrations from MODIS and MERIS during the 4 seasons of the year (Figure 6). The Rhone River plume area was characterized by high chlorophyll-a concentrations, in contrast with the eastern portion of the study area, where concentrations in chlorophyll-a remained very low. The models captured the west-east gradient in surface chlorophyll-a concentration. The BoM seems to be a transition zone between the rich water of the Rhone River plume and the oligotrophic water located eastward. The models reproduced well the spatial extension of the Rhone River intrusion plume, except very close to the mouth, where the concentrations were underestimated by the model. This underestimation was caused by the inability to distinguish terrestrial and fluvial chlorophyll-a from marine chlorophyll-a in the remote sensing data. The marine chlorophyll-a is the only one included in the model.

The chlorophyll-a remote sensing data were highly variable in both the spatial distribution and the concentrations; this variability was well captured by both versions of the model. In winter, both version of the model reproduced well the chlorophyll-a concentrations in the BoM, but the concentrations associated with the Rhone River plume were underestimated. During spring and summer, the model fields deviate from the observations in overestimating the remotely sensed chlorophyll-a concentrations in the Rhone River plume. Nevertheless, the spring remains the season during which the surface layer was the most productive and the chlorophyll-a concentration was at its maximum in the BoM. The summer period had lower concentrations than in the spring in the surface layer, particularly in the BoM. The fall was characterized by the lowest chlorophyll-a concentrations of the year in the entire study area; the signal of the Rhone River was weak in comparison with the three other seasons.

The spatial comparisons exhibited little discrepancy between the two versions of the model. During spring and summer, the chlorophyll-a signal associated with the Rhone River plume was lower and more similar to the observations in the model with the P cycle both spatially and in intensity. Therefore, the addition of the P cycle causes an improvement in the representation of the effects of the Rhone River on chlorophyll-a concentrations.

**3. Influence of short time scale shelf processes forcing on co-limitation functioning**

In the evaluation of the representation of short time scale shelf processes in the coupled models, we focused on the year 2008, for which hydrodynamic shelf processes off Marseille were analyzed in detail by Pairaud et al. [6]. The Marseille coastal area was under strong forcing influences (upwelling events, intrusion of Rhone River plume in the BoM), which induced significant daily variability in the biogeochemical variables. Figure 7 presents the observations and daily model results for 2008 used to evaluate the ability of the model to reproduce processes over short time scales. In early April, a peak in chlorophyll-a was present in the model results and in the remote sensing and in-situ data, and we noticed strong discrepancies between the in-situ data and the remote sensing data. The model values were between those of the remote sensing and in-situ data, leading to the assumption that chlorophyll-a concentrations could be underestimated from remote sensing data in spring. In summer, the BoM experienced enrichment events, when modeled chlorophyll-a concentrations increased from 0.5 to 1.5 µmol.L⁻¹. Those events matched the observations very well, except for the event of late August 2008. The chlorophyll-a concentrations remained low during fall with slight variations in intensity for both the models and observations. Weak discrepancies were observed at short time scales between the two versions of the model, except in summer, when enrichment events were clearly less intense in the model with the P cycle. To explain these discrepancies, we focused on phytoplankton limitation by P and N.

Figure 8 presents phytoplankton limitation by nutrients calculated using the cell quota formulation in the surface and bottom layers at the Somlit station. As the limitation increases in value toward 1, the nutrient limitation on phytoplankton increases. Temperature was also presented on the same graph (Figure 8) to evaluate the links between nutrient co-limitation and temperature. The results of the model showed that phytoplankton limitation by nutrients was greater in the surface layer (Figure 8A) than in the bottom layer (Figure 8B), as was the variability in limitation. As shown in Figure 3, high concentrations of nutrients were available from November to March due to winter vertical mixing, so that phytoplankton limitation by nutrients remained very low in the surface and bottom layers. An increase in limitation occurred from March to the beginning of May, corresponding to the spring bloom. The end of the spring bloom was marked by the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Target diagrams of the model/data comparison. Target diagrams of a set of variables for the model with the P cycle (A–C) and without the P cycle (B–D). The model/data comparison is made using in-situ data for the years 2009 to 2011 at the surface (A–B) and at 55 m (bottom C–D) at the Somlit station. The normalized bias (dimensionless) was computed versus the sign unbiased normalized RMSD (dimensionless).
doi:10.1371/journal.pone.0080012.g005</figcaption></figure>

exhaustion of nutrient supplies in the surface layer, the stabiliza-
tion of limitation at high levels and the stabilization of the
chlorophyll-a concentration at approximately 0.3 µg.L⁻¹ at the
surface (Figure 7). Short, strong decreases in limitation happened
during summer at the surface. Then, the fall was associated with a
decrease in limitation by nutrients. Therefore, phytoplankton
limitation by nutrients had an apparent seasonal signal, and
increases in the chlorophyll-a concentration in the surface layer
were linked with short decreases in phytoplankton nutrient
limitation. However, phytoplankton nutrient limitation was not
directly correlated with chlorophyll-a because temperature and
light limitation were also considered to affect chlorophyll-a
seasonality.

The model results also showed that phytoplankton were more
limited by P (grey points) in the surface layer and more by N (green
points) in the bottom layer. However, switches between P and N
limitation occurred several times during the year, highlighting the

strong variability of ecosystem functioning most likely due to
physical shelf processes.

**3.1. Influence of upwelling.** The Marseille coastal area is strongly impacted by upwelling events associated with strong upwards vertical velocity (maximum of 5 cm.s⁻¹ during an upwelling event in November 2008), which led to an important decrease in surface temperature [6]. Significant upwelling events occurred during the summer of 2008, and these strong temperature decreases were always associated with N limitation (Figure 8A). This result was consistent with the modeled N limitation of phytoplankton development in bottom layer (Figure 8B).

The maps in Figure 9 present a comparison of surface
temperature, chlorophyll-a concentration and phytoplankton
nutrient limitation during an upwelling event in July 2008. During
this event, the area impacted by the ascent of deep water was
characterized by low temperature (T<20°C) (Figure 9A) and
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Mean surface chlorophyll concentration for winter, spring, summer and fall. Mean surface chlorophyll concentration for winter (December–January–February), spring (March–April–May), summer (June–July–August) and fall (September–October–November) for the years 2007 to 2011 from the remote sensing data (MODIS and MERIS) and the simulation (with and without the P cycle). The model results presented were averaged across the 10 first meters of the water column.
doi:10.1371/journal.pone.0080012.g006</figcaption></figure>

limitation by N rather than P (Figure 9B). The lower temperature was localized to the two primary upwelling points (“Cote Bleue” and “Calanques”), but the majority of the Marseille coastal area was impacted by upwelling as described by Pairaud et al. [6]. The remote sensing data and modeled chlorophyll-a concentrations were in good agreement (Figure 9 C–D); the concentrations were

high in the two primary upwelling points in response to the
nutrient enrichment induced by the upwelling. A weak increase in
chlorophyll-a was observed far from the two main upwelling spots
but still located in the temperature-impacted upwelling area. This
could be explained by the rapid uptake of nutrients near the
upwelling points. The nutrients were consumed by phytoplankton

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Surface chlorophyll concentrations for the year 2008. Evolution of the surface chlorophyll concentrations from the remote sensing data of MODIS (green point) and MERIS (blue point), in-situ data (blue point) and the simulation with the P cycle (black line) and without the P cycle (blue line) at the Somlit station.
doi:10.1371/journal.pone.0080012.g007</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Evolution of temperature and phytoplankton nutrient limitation for the year 2008. Phytoplankton limitation by P (grey point) and N (blue point) and model temperature (blue line) for the year 2008 at the Somlit station at the surface (A) and at 60 m (B). doi:10.1371/journal.pone.0080012.g008</figcaption></figure>

before being transported away from the upwelling points due to
the strong nutrient limitation of phytoplankton in the summer and
in the surface layer. These results explained why the surface
temperature was strongly impacted and the surface chlorophyll-a
concentration more weakly impacted at the Somlit station during
upwelling events (Figure 7–8A).

**3.2. Influence of the intrusion of water from the Rhone River in the BoM.** The intrusion event of water from the Rhone River into the BoM that was studied by Pairaud et al. [6] occurred in the second part of June 2008 and was detected in the model results at the Somlit station as a decrease in salinity (Figure 10A). The study of the summer decreases in phytoplankton nutrient limitation revealed that they were linked to low salinity events. Rich Rhone River water induced an increase in nutrients and thus a decrease in phytoplankton nutrient limitation (Figure 10A). Both versions of the model were in agreement with the observed chlorophyll-a concentrations (Figure 10B), highlighting the strong biological response to the Rhone River water intrusion events. The intensity of the increases in chlorophyll-a were lower in the version of the model with the P cycle during most of the low salinity events. Therefore, according to the results of the model, phytoplankton was limited by P during most of the low salinity events. This finding was in agreement with spatial comparisons of chlorophyll-a concentrations in the Rhone River plume (Figure 6), which were less extended and better represented in the model version with the P cycle. Therefore, the P cycle slightly improved the results, such as during the low salinity event in mid-July where a net improvement was noticed.

Discussion

Due to a worldwide movement toward ecosystem-based
management [79], the demand for quantitative tools to support
ecosystem-based management initiatives increases [79,80]. There
is now a proliferation of end-to-end ecosystem models which
attempt to represent the entire ecological system and the
associated abiotic environment [81]. This kind of models has
been developed in various areas [82,83] and some are listed in the
MEECE European project (http://www.meece.eu). However this
kind of models has lots of parameters, needs a lot of input data,
human and material resources to be developed, which make them
difficult to transpose in numerous areas.

In this study, we preferred a different modeling approach and
we developed a model adapted to targeted questions: study of the
impact of the physical processes on low trophic levels and study of
the inputs of chemical contaminants in the trophic chain. So, we
chose to develop a physical biogeochemical coupled model with a
high spatial and temporal resolution and with an on-line coupling
technique in order to describe finely the link between physics and
biogeochemistry. This coupled model focused on the low trophic
level in order to study the key processes which permit the input of
contaminant in the trophic chain. The biogeochemical model with
a low number of parameters relatively to end-to-end models is
easier to adapt to other oligotrophic coastal area influenced by
large river inputs like the Gulf of Gabes in Tunisia and the Nile
delta in eastern Mediterranean Sea. Indeed, the biogeochemical
model ECO3M-MASSILIA-P was modified from an initial model
[40,41,84].
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Upwelling event of 15/07/2008. Simulated temperature (A), simulated phytoplankton limitation by P (red) and N (blue) (B), chlorophyll-a concentrations from MODIS and simulated chlorophyll-a concentration (average across the first 10 meters of depth) compared for 15/07/2008. doi:10.1371/journal.pone.0080012.g009</figcaption></figure>

During the biogeochemical model development, the following question of the complexity level of the biogeochemical model was asked: What was the level of detail necessary and sufficient to

<figure><img src="image_12.png" /></figure>

describe this coastal ecosystem? The Mediterranean Sea is generally known to be P limited; thus, we decided to add the P cycle to the model. In choosing the most reliable model, particular attention was paid to the accuracy of the model and the assessment of the shortcomings of the two new versions of the biogeochemical model.

**1. Model accuracy assessment**

The models simulated 5 years, from 2007 to 2011. The model development and parameterizations were performed for 2007 and 2008, whereas model accuracy assessment of both final model versions was performed for 2009 to 2011, which allowed a robust evaluation of the model results.

The results showed that both versions of the model had significant ability to simulate the biogeochemical functioning of the Marseille coastal area. Indeed, the models reproduced with relatively good accuracy the seasonal and the inter-annual dynamics of chlorophyll-a and nutrient concentrations. Nevertheless, the chlorophyll-a concentration was slightly overestimated by the model during the spring bloom and in summer. The models reproduced the spatial pattern of this area well, with a good representation of the west-east gradient in chlorophyll-a concentration associated with the Rhone River plume extensions. The general ability assessment of model results also showed that there was strong variability in biogeochemical concentrations at short time scales (from a few hours to several days), mostly in spring and summer. The model was able to reproduce different types of shelf events, such as upwelling or the intrusion of Rhone River water into the BoM, and their associated biogeochemical response. The use of target diagrams allowed a visual comparison of statistical indicators. The results appeared more accurate at the surface than at the bottom, and the discrepancies between the two versions of the biogeochemical model were rare.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Impact of the intrusion of Rhone River water on phytoplankton nutrient limitation and chlorophyll-a concentrations. Evolution from May to mid-August 2008 at the Somlit station of the simulated phytoplankton limitation by P (grey point) and N (blue point) and the simulated salinity (red line) (A), the surface chlorophyll-a concentrations from remote sensing observations of MODIS (green point) and MERIS (blue point), the in-situ data (bleu point) and the simulation with the P cycle (black line) and without the P cycle (blue line) (B). doi:10.1371/journal.pone.0080012.g010</figcaption></figure>
Nevertheless, the implementation of the P cycle in the biogeochemical model improved the agreement between the model predictions and the surface chlorophyll-a observations for a given event. In addition, the model with the P cycle performed well at reproducing the plume and processes associated with the Rhone River. However, in the bottom layer, its use increased the error in predicting the chlorophyll-a concentration due to the overestimation of the PO₄ concentration. We attributed the improved fit in the surface chlorophyll-a concentrations to the ability of the model to reproduce the P limitation in the spring and summer periods. In winter, phytoplankton production was principally limited by light because the strong vertical mixing brought nutrients in excess to the surface layer. During this period, the P cycle did not impact the surface chlorophyll-a concentration.

## 2. Model complexity

The general assessment of the results of both versions of the model showed that the accuracies of the models with and without the P cycle were similar. The addition of the P cycle to the biogeochemical model improved the description of the ecosystem functioning, but it also injected errors.

The errors associated with the P cycle came from different sources. The primary difficulty was that the measured PO₄ concentrations in the Marseille coastal area were often lower or very close to the detection limit. At the Somlit station in the surface layer, the concentration of PO₄ ranged from 0 to 0.3 µmol.L⁻¹, with a mean value of 0.03 µmol.L⁻¹ for the years 2009 to 2011, which indicated that concentrations of PO₄ were generally lower than or close to their detection limits, as the precision of the PO₄ concentrations was ±0.02–0.03 µmol.L⁻¹. Therefore, comparisons between the model and the observations were not possible. Then, as previously stated by several authors [85–87], there is a critical need for new experimental data examining N and P dynamics under different and extreme input ratios and growth rates [48]. Thus, the need for improved biogeochemical understanding of N/P co-limitation increased the difficulty of modeling the P cycle. Furthermore, Poggiale et al. [88] demonstrated that models were sensitive to uptake formulations; even if two formulas provided similar values, large numerical differences in the stability criteria may occur. Thus, errors were certainly introduced by the choice of the parameter values and the mathematical functions representing the P cycle.

Even if it appears plausible from a biological a priori standpoint that more complex models should mirror reality better, other recent studies have suggested that it is not always the case. A meta-analysis of mechanistic aquatic biogeochemical models also found that increased model complexity did not improve the fit [37]. Kriest et al. [89] suggested that increasing complexity of unturned, unoptimized models that were simulated with parameters commonly used in large-scale model studies did not necessarily improve performance. Los et al. [90] also indicated that adding more complexity did not necessarily improve the quality of the model results in terms of their ability to reproduce measurements and hence their applicability as prognostic tools. Instead, they argued that there should be a balance between ecological and physical resolution in relation to the specific question to be addressed. Crout et al. [91] evaluated models using model reduction for three models. In all cases, they identified reduced models that outperformed the more complicated ones, “suggesting some over-parameterization has occurred during model development” [91].

In the study of the Marseille coastal area, we followed an inverse approach, with the simplest model (model without P cycle) at the beginning and consecutively increasing model complexity with the

P cycle. The addition of the P cycle increased the computation time by approximately 30% and did not automatically improve the fit with the observations. Therefore, the model version (with or without the P cycle) should be chosen depending on the question to be addressed and by considering the general accuracy assessment achieved in this work.

## 3. The Marseille Coastal area: a challenging zone for a modeling approach

The region of interest is very complex (biologically and dynamically) and was thus a challenge to model. Principally, 4 difficulties were encountered: (i) there was a strong west-east gradient in oligotrophy, namely the eastern portion of the studied area was mainly oligotrophic, while the western portion was more productive; (ii) few observations were available for comparisons with the model results; (iii) the studied region was a coastal area highly influenced by large river inputs (the Rhone River and the Marseille urban rivers); and (iv) the hydrodynamic shelf processes were extremely variable at short time scales, and they strongly impacted most of the coastal domain.

The BoM was generally oligotrophic and characterized by low concentrations of chlorophyll-a and nutrients, with mean concentrations of chlorophyll-a, PO₄ and DIN of 0.4 µg.L⁻¹, 0.03 µmol.L⁻¹ and 1 µmol.L⁻¹, respectively. In comparison, in the surface layer of the Baltic Sea, the mean concentrations of chlorophyll-a, PO₄ and DIN were approximately 3.5 µg.L⁻¹, 0.33 µmol.L⁻¹ and 3.6 µmol.L⁻¹ [92]. As in the Baltic Sea, a modeling study in the Channel and Southern Bight of the North Sea showed mean concentrations of chlorophyll-a of 5.35 µg.L⁻¹, mean PO₄ of 0.82 µmol.L⁻¹ and mean DIN of 17.57 µmol.L⁻¹ [77]. Therefore, in the BoM, the concentrations were ten times less than the concentration range in the Channel and Southern Bight of the North Sea and in the Baltic Sea. The BoM contrasted with the Rhone River plume, which was characterized by high concentrations of nutrients and chlorophyll-a. Thus, the large range of trophic conditions in this small model domain put the model to the test.

The scarcity of observations to compare with the model results was a problem during the model development phase. The Somlit station was sampled only twice monthly, and as shown by the model results, many processes occurred at shorter time scales. One risk to avoid was to tune the model to obtain good results in comparison with the in-situ data at the Somlit station. We also evaluated the spatial pattern of the processes reproduced by the model and compared them with remote sensing data, which provided spatial patterns and gave daily information (during non-cloudy conditions). Nevertheless, the errors associated with the remote sensing data were high.

Coastal areas are highly influenced by cross-shore and earth/sea exchange. The Marseille coastal area has two large open boundary conditions (OBC). Thus, a good representation of the OBC was essential to accurately describe the cross-shore exchanges. Coastal modeling also requires long-term monitoring of inputs to represent the discharge and concentrations associated with each river well. Although the Rhone River was well-studied and the daily discharge and concentrations of the primary state variables were available, it was necessary to construct from the literature or from the few available observations the concentrations for the missing variables for the Marseille urban rivers and the Rhone River. Progress could be made in the description of the concentrations of the rivers by monitoring each river to better describe the organic matter and nutrient inputs from the rivers.

In addition, the low impact of the addition of the P cycle to the model suggested that physics drives the majority of the coupled
model results. Hydrodynamic shelf processes appeared extremely variable at short time scales, and they strongly impacted most of the coastal domain. In coastal areas, hydrodynamics were of primary importance for coupled modeling. An overestimation of the plume extension of the Rhone River by the hydrodynamic model induced a strong error in the chlorophyll-a concentration at the Somlit station, which demonstrated that the choice of the physical model is very important for biogeochemistry, as stated in many studies [93,94].

## 4. Contribution of the model results to the understanding of the Marseille coastal area

As previously discussed, although this coastal ecosystem was complex and thus a challenge to model, the coupled model performed well at reproducing the primary characteristics and processes. Therefore, the coupled model provided interesting information on the functioning of this coastal ecosystem.

First, the temporal evolution of surface chlorophyll-a clearly distinguished the four main seasons. Spring was characterized by high concentrations of chlorophyll-a (spring bloom) and was the most productive season. In contrast, the fall appeared to be the least productive period. This finding is not in accordance with general knowledge of the Northwestern Mediterranean Sea as a small autumnal bloom is supposed to occur during this period. Our coupled model reproduced the autumnal bloom only for 2009.

From a spatial point of view, the Marseille coastal area was characterized by a strong west-east oligotrophy gradient. Outside the ROFI, the coastal domain near Marseille generally remained oligotrophic, whereas at the mouth of the Rhone River, the phytoplankton biomass was always high [95]. Nevertheless, the Rhone River was limited by P because the NO₃·PO₄ ratio of the Rhone River inputs was approximately 65–80 [26,27]. This P limitation was also detectable in the model results because the chlorophyll-a concentrations in the Rhone River plume were higher in the model without the P cycle than in the model with the P cycle, better mimicking the observations.

The urban rivers and the Marseille WWTP delivered high concentrations of nutrients and organic matter, but the eutrophication risk remained extremely low in the BoM because the Urban River inputs and sewage from the WWTP were rapidly diluted by strong hydrodynamic events [6] due to the short residence time. Nevertheless, the Gulf of Fos was characterized by shallow depths and higher chlorophyll-a concentrations (up to 10 µmol·L⁻¹ [34]), which indicated that the eutrophication risk was higher. Previous studies [96] demonstrated that the greatest risk of algal blooms occurs during periods of calm weather because the Gulf of Fos is flushed intermittently but strongly by wind-driven lateral circulation [34].

The coupled model permitted the evaluation of the eutrophication risk in the BoM, but the model was also very useful in contextualizing the in-situ data. Indeed, the Somlit station is impacted by many shelf events (e.g., upwelling and the intrusion of Rhone River water). For example, an increase in nutrient and chlorophyll-a concentrations associated with a decrease in salinity could be caused by the intrusion of Rhone River water, heavy rainfall events or an extension of the urban river plumes.

The coupled model presented a significant advantage over satellite images and in-situ measurements by allowing the study of the processes at short time scales and in three dimensions. It was proven to be a pertinent tool in the study of the biological response associated with physical processes. The study of upwelling events in summer 2008 showed that a large part of the Marseille coastal area was impacted by upwelling events. However, the chlorophyll-

a response to the nutrient enrichment caused by upwelling was mostly noticeable at the two main upwelling points ("Cote Bleue" and "Calanques"). Another interesting event was the intrusion of Rhone River water into the BoM. The influence of the Rhone River water in the BoM has been observed hydrodynamically [6,21], but the only impacts on the ecosystem that had previously been investigated were in the CDOM concentrations [22]. The coupled model revealed that the intrusion of Rhone River water into the BoM induced an obvious increase in chlorophyll-a concentrations.

Finally, further studies are necessary to better explore and characterize the biogeochemical functioning of this coastal ecosystem. The coupled model allowed the quantification of the impact of atmospheric inputs (organic matter and nutrients) on the biogeochemical functioning of this coastal ecosystem. Preliminary results (not shown) suggested that atmospheric deposits might be negligible, except during episodes of heavy rainfall. The initial results concerning the intrusion of Rhone River water also raised new questions: What caused the eastward transport of Rhone River water? What quantities of nutrients are transported into the BoM during Rhone River water intrusion events? Thus, the coupled model could be a useful tool in quantifying the relative impact of the Rhone River and the Urban Rivers on the biogeochemical functioning of the BoM using a mass balance approach.

## Conclusion

In this paper, we presented the development and the evaluation of a 3D coupled physical-biogeochemical model for the Marseille coastal area. The high-resolution coupled model demonstrated its ability to recreate realistic situations. It underlined that the biogeochemical functioning of the BoM was very complex because the physical processes and river inputs were important drivers of the biogeochemical model results. Indeed, the ecosystem of the BoM quickly switched between oligotrophic conditions and important enrichment events due to external forcing.

The P cycle was added to the biogeochemical model because this element was well known to limit biological production in the Mediterranean Sea. The addition of P improved the description of ecosystem functioning and in some cases improved the model results, but it also introduced errors into the model. In using the model results, awareness of its shortcomings and positive features is important. Therefore, the general accuracy assessment in this paper and the questions addressed should be considered in choosing the version of the biogeochemical model (with and without the P cycle).

The coupled model appeared to be a useful tool for the analysis of processes and the estimation of budgets in a very dynamic environment where it is difficult to extrapolate from discrete measurements. Thus, despite its imperfections, large quantities of information are available in the results of the coupled model. The model data could help researchers to design field campaigns by better anticipating the biogeochemical front and gradient. In addition, the model results allow better evaluation of the impact of urban inputs on the coastal area and could thus help policy managers propose solutions.

In conclusion, this study suggested that improvements in the description of hydrodynamics and terrestrial inputs should be preferred over increasing the complexity of the biogeochemical model in this coastal oligotrophic area.
## Supporting Information

**Annex S1 Equations of the biogeochemical model of the
Marseille coastal area (ECO3M-MASSILIA-P).**
(DOCX)

**Annex S2 Parameters.**
(PDF)

**Annex S3 Spin up of the MARS3D-RHOMA/ECO3M-MASSILIA coupled model for the year 2008 at the Somlit station at the surface.** "Model 1" corresponding to the summer initial conditions and "Model 2" corresponding to the winter initial conditions.
(DOCX)

**Annex S4 Maps of the Marseille city inputs.**
(TIF)

**Annex S5 Quality of the remote sensing data.**
(DOCX)

## Acknowledgments

The authors acknowledge the staff of the Somlit national network for the littoral observations, (INSU-CNRS) for providing the data time series,

## References

1. Le Tissier MD, Buddemeier R, Parslow J, Swaney DP, Crossland CJ, et al. (2006) The role of the coastal ocean in the disturbed and undisturbed nutrient and carbon cycles A management perspective. LOICZ, Geesthacht, Germany.

2. Goberville E, Beaugrand G, Sautour B, Tréguer P, Somlit T (2010) Climate-driven changes in coastal marine systems of western Europe. Marine Ecology Progress Series 408: 129–147.

3. Behrenfeld MJ, O'Malley RT, Siegel DA, McClain CR, Sarmiento JL, et al. (2006) Climate-driven trends in contemporary ocean productivity. Nature 444: 752–755.

4. INSEE (2009) Résultats du recensement de la population. Available: http://www.recensement.insee.fr/chiffresCles.action?zoneSearchField=MARSEILLE&ecoZone=241300391-GFP&idTheme=3. Accessed 17 April 2013.

5. Durrieu de Madron X, Guieu C, Sempéré R, Conan P, Cossa D, et al. (2011) Marine ecosystems' responses to climatic and anthropogenic forcings in the Mediterranean. Progress in Oceanography 91: 97–166.

6. Pairaud IL, Gatti J, Bensoussan N, Verney R, Garreau P (2011) Hydrology and circulation in a coastal area off Marseille: Validation of a nested 3D model with observations. Journal of Marine Systems 88: 20–33.

7. Coste B (1974) Role des apports nutritifs mineraux rhodaniens sur la production organique des eaux du Golfe du Lion. Tethys 6: 727–740.

8. Morel A, Bricaud A, André JM, Pelaez-Hudlet J (1990) Spatial-temporal evolution of the Rhone plume as seen by CZCS imagery. Consequences upon the primary production in the Gulf of Lions. Water Pollution Research Reports In: Martin: 45–62.

9. Petrenko A, Dufau C, Estournel C (2008) Barotropic eastward currents in the western Gulf of Lion, north-western Mediterranean Sea, during stratified conditions. Journal of Marine Systems 74: 406–428.

10. Allou A, Forget P, Devenon J-L (2010) Submesoscale vortex structures at the entrance of the Gulf of Lions in the Northwestern Mediterranean Sea. Continental Shelf Research 30: 724–732.

11. Millot C (1990) The Gulf of Lions' hydrodynamics. Continental Shelf Research 10: 885–894.

12. Millot C, Wald L (1980) The effect of mistral wind on the ligurian current near Provenc. Oceanologica Acta 3: 399–402.

13. Petrenko A (2003) Variability of circulation features in the Gulf of Lion NW Mediterranean Sea. Importance of inertial currents. Oceanologica Acta 26: 323–338.

14. Gatti J (2008) Intrusions du courant nord méditerranéen sur la partie est du plateau continental du golfe du Lion University Aix-Marseille II.

15. El Sayed MA, Aminot A, Kerouel R (1994) Nutrients and trace metals in the northwestern Mediterranean under coastal upwelling conditions. Continental Shelf Research 14: 507–530.

16. Minas HL (1968) A propos d'une remontée d'eaux "profondes" dans les parages du golfe de Marseille (oct. 1964), conséquences biologiques. Cahiers Océanographiques 20: 647–674.

17. Broche P, Devenon J-L, Forget P, de Maistre J-C, Naudin J-J, et al. (1998) Experimental study of the Rhone plume. Part I: physics and dynamics. Oceanologica Acta 21: 725–738.

IFREMER for providing the ocean color remote sensing data computed with the OC5 algorithm and the Compagnie Nationale du Rhone for the data on Rhone River discharges and Meteo France. The Rhone concentration data were provided by the national MOOSE program (Mediterranean Oceanic Observing System on Environment) and the Service d'Observation of the Mediterranean Institute of Oceanography (MIO). The authors gratefully acknowledge N. Garcia, V. Lagadec and M. Fornier for analytical and field assistance. This study was part of the "MERMEX WP3-C3A" and international "IMBER" projects. We are grateful to DEA-MPM (Direction des Eaux et Assainissement- Marseille Provence Metropole) and SERAM (Société d'Exploitation du Réseau d'Assainissement de Marseille) for providing the Marseille urban river and WWTP discharge. Many thanks also to Benedicte Thouvenin, Pierre Garreau, Francis Gohin, Romaric Verney, Cassandre Jany and Cedric Garnier for fruitful discussion. We thank American Journal Experts for English corrections and the two anonymous reviewers for their helpful comments to improve this paper.

## Author Contributions

Analyzed the data: MF CP IP PR RF VF PL. Wrote the paper: MF CP IP RF VF PR PL. Conceived and designed modelling configurations: MF CP IP RF VF PL PR.

18. Estournel C, Broche P, Marsaleix P, Devenon JL, Auclair F, et al. (2001) The Rhone River plume in unsteady conditions: numerical and experimental results. Estuarine, Coastal and Shelf Science 53: 25–38.

19. Arnoux-Chiavassa S, Rey V, Fraunie P (2003) Modeling 3D Rhône river plume using a higher order advection scheme. Oceanologica Acta 26: 299–309.

20. Reffray G, Fraunie P, Marsaleix P (2004) Secondary flows induced by wind forcing in the Rhone region of freshwater influence. Ocean Dynamics 54: 179–196.

21. Gatti J, Petrenko A, Devenon J-L, Leredde Y, Ulses C (2006) The Rhone river dilution zone present in the northeastern shelf of the Gulf of Lion in December 2003. Continental Shelf Research 26: 1794–1805.

22. Para J, Coble PG, Charrière B, Tedetti M, Fontana C, et al. (2010) Fluorescence and absorption properties of chromophoric dissolved organic matter (CDOM) in coastal surface waters of the northwestern Mediterranean Sea, influence of the Rhône River. Biogeosciences 7: 4083–4103.

23. Elser JJ, Bracken MES, Cleland EE, Gruner DS, Harpole WS, et al. (2007) Global analysis of nitrogen and phosphorus limitation of primary producers in freshwater, marine and terrestrial ecosystems. Ecology letters 10: 1135–1142.

24. Moutin T, Raimbault P, Golterman HL, Coste B (1998) The input of nutrients by the Rhône river into the Mediterranean Sea: recent observations and comparison with earlier data. Hydrobiologia 373–374: 237–246.

25. Van Wambeke F, Christaki U, Giannakourou A, Moutin T, Souvemergoglou K (2002) Longitudinal and vertical trends of bacterial limitation by phosphorus and carbon in the Mediterranean Sea. Microbial ecology 43: 119–133.

26. Ludwig W, Dumont E, Meybeck M, Heussner S (2009) River discharges of water and nutrients to the Mediterranean and Black Sea: Major drivers for ecosystem changes during past and future decades? Progress in Oceanography 80: 199–217.

27. Ludwig W, Bouwman AF, Dumont E, Lespinas F (2010) Water and nutrient fluxes from major Mediterranean and Black Sea rivers: Past and future trends and their implications for the basin-scale budgets. Global Biogeochemical Cycles 24: GB0A13.

28. Markaki Z, Loçe-Pilot MD, Violaki K, Benyahia L, Mihalopoulos N (2010) Variability of atmospheric deposition of dissolved nitrogen and phosphorus in the Mediterranean and possible link to the anomalous seawater N/P ratio. Marine Chemistry 120: 187–194.

29. Thingstad TF, Zweifel UL, Rassoulzadegan F (1998) P limitation of heterotrophic bacteria and phytoplankton. Limnol Oceanogr 43: 88–94.

30. Pinazo C, Marsaleix P, Millet B, Estournel C, Vehil R (1996) Spatial and temporal variability of phytoplankton biomass in upwelling areas of the northwestern mediterranean: a coupled physical and biogeochemical modelling approach. Journal of Marine Systems 7: 161–191.

31. Tusseau-Vuillemin M-H, Mortier L, Herbaut C (1998) Modeling nitrate fluxes in an open coastal environment (Gulf of Lions): Transport versus biogeochemical processes. Journal of Geophysical Research 103: 7693.

32. Auger PA, Diaz F, Ulses C, Estournel C, Neveux J, et al. (2011) Functioning of the planktonic ecosystem on the Gulf of Lions shelf (NW Mediterranean) during spring and its impact on the carbon deposition: a field data and 3-D modelling combined approach. Biogeosciences 8: 3231–3261.
33. Campbell R, Diaz F, Hu Z, Doglioli A, Petrenko A, et al. (2013) Nutrients and plankton spatial distributions induced by a coastal eddy in the Gulf of Lion. Insights from a numerical model. Progress in Oceanography 109: 47–69.

34. Fontana C, Grenz C, Pinazo C, Marsaleix P, Diaz F (2009) Assimilation of SeaWiFS chlorophyll data into a 3D-coupled physical–biogeochemical model applied to a freshwater-influenced coastal zone. Continental Shelf Research 29: 1397–1409.

35. Franks PJS (2009) Planktonic ecosystem models: perplexing parameterizations and a failure to fail. Journal of Plankton Research 31: 1299–1306.

36. Ourmières Y, Brasseur P, Lévy M, Brankart J-M, Verron J (2009) On the key role of nutrient data to constrain a coupled physical–biogeochemical assimilative model of the North Atlantic Ocean. Journal of Marine Systems 75: 100–115.

37. Arhonditsis G, Brett M (2004) Evaluation of the current state of mechanistic aquatic biogeochemical modeling. Marine Ecology Progress Series 271: 13–26.

38. Baklouti M, Diaz F, Pinazo C, Faure V, Quéguiner B (2006) Investigation of mechanistic formulations depicting phytoplankton dynamics for models of marine pelagic ecosystems and description of a new model. Progress in Oceanography 71: 1–33.

39. Baklouti M, Faure V, Pawlowski L, Sciandra A (2006) Investigation and sensitivity analysis of a mechanistic phytoplankton model implemented in a new modular numerical tool (Eco3M) dedicated to biogeochemical modelling. Progress in Oceanography 71: 34–58.

40. Faure V, Pinazo C, Torreton J-P, Jacquet S (2010) Modelling the spatial and temporal variability of the SW lagoon of New Caledonia I: a new biogeochemical model based on microbial loop recycling. Marine pollution bulletin 61: 465–479.

41. Faure V, Pinazo C, Torreton J-P, Douillet P (2010) Modelling the spatial and temporal variability of the SW lagoon of New Caledonia II: realistic 3D simulations compared with in situ data. Marine pollution bulletin 61: 480–502.

42. Harmelin-Vivien M, Loizeau V, Mellon C, Beker B, Arlhac D, et al. (2008) Comparison of C and N stable isotope ratios between surface particulate organic matter and microphytoplankton in the Gulf of Lions (NW Mediterranean). Continental Shelf Research 28: 1911–1919.

43. Beker B, Romano J-C, Arlhac D (1999) Le suivi de la variabilité des peuplements phytoplanctoniques en milieu marin littoral: le golfe de Marseille en 1996-1997. Oceanis 25: 395–415.

44. Gutiérrez-Rodríguez A, Latasa M, Scharek R, Massana R, Vila G, et al. (2011) Growth and grazing rate dynamics of major phytoplankton groups in an oligotrophic coastal site. Estuarine, Coastal and Shelf Science 95: 77–87.

45. Halsband-Lenk C, Hirche H-J, Carlotti F (2002) Temperature impact on reproduction and development of congener copepod populations. Journal of Experimental Marine Biology and Ecology 271: 121–153.

46. Rose JM, Caron DA (2007) Evidence Does low temperature constrain the growth rates of heterotrophic protists? and implications for algal blooms in cold waters. Limnology and Oceanography 52: 886–895.

47. Saito MA, Goepfert TJ, Ritt JT (2008) Some thoughts on the concept of colimination: Three definitions and the importance of bioavailability. Limnology and Oceanography 53: 276–290.

48. Bougaran G, Bernard O, Sciandra A (2010) Modeling continuous cultures of microalgae colimated by nitrogen and phosphorus. Journal of theoretical biology 265: 443–454.

49. Rhee GY (1978) Effects of N:P atomic ratios and nitrate limitation on algal growth, cell compositions, and nitrate uptake. Limnology and Oceanography 23: 10–25.

50. Geider RJ, MacIntyre HL, Kana TM (1998) A Dynamic Regulatory Model of Phytoplanktonic Acclimation to Light, Nutrients, and Temperature. Limnology and Oceanography 43: pp. 679–694.

51. Smith CL, Tett P (2000) A depth-resolving numerical model of physically forced microbiology at the European shelf edge. Journal of Marine Systems 26: 1–36.

52. Faure V, Pinazo C, Torréton J-P, Douillet P (2006) Relevance of various formulations of phytoplankton chlorophyll a:carbon ratio in a 3D marine ecosystem model. Comptes rendus biologies 329: 813–822.

53. Droop MR (1968) Vitamin B12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri. Journal of Marine Biology Association 48: 689–733.

54. Nicolle A, Garreau P, Liorzou B (2009) Modelling for anchovy recruitment studies in the Gulf of Lions (Western Mediterranean Sea). Ocean Dynamics 59: 953–968.

55. Pinazo C, Fraysse M, Fuchs R, Thouvenin B, Pairaud I, et al. (2012) Modelling the Gulf of Lions coastal area with a 3D coupled physical and biogeochemical approach. 50th ECSA Conference. Venice, Italy.

56. Lazzari P, Solidoro C, Ibello V, Salon S, Teruzzi A, et al. (2012) Seasonal and inter-annual variability of plankton chlorophyll and primary production in the Mediterranean Sea: a modelling approach. Biogeosciences 9: 217–233.

57. Ibanez C, Pont D, Prat N (1997) Characterization of the Ebre and Rhone Estuaries: A Basis for Defining and Classifying Salt-Wedge Estuaries. Limnology and Oceanography 42: pp. 89–101.

58. Pont D, Simonnet J-P, Walter A V (2002) Medium-term Changes in Suspended Sediment Delivery to the Ocean: Consequences of Catchment Heterogeneity and River Management (Rhône River, France). Estuarine, Coastal and Shelf Science 54: 1–18.

59. Sempere R, Charriere B, Van Wambeke F, Cauwet G (2000) Carbon inputs of the Rhone River to the Mediterranean Sea: Biogeochemical implications. Global Biogeochemical Cycles 14: 669–681.

60. Déliat G (2001) Matière organique dissoute des zones côtières: sources, distribution et biodégradabilité Université Paris VI.

61. Lønborg C, Davidson K, Álvarez-Salgado XA, Miller AEJ (2009) Bioavailability and bacterial degradation rates of dissolved organic matter in a temperate coastal area during an annual cycle. Marine Chemistry 113: 219–226.

62. Jany C, Thouvenin B (2012) Récapitulatif des hypothèses utilisées pour les rejets de MES dans les simulations pour METROC, Rapport Ifremer, RST.ODE/LER/PAC/12-19.

63. Gouze E (2008) Bilan de matière de l'étang de Berre : influence des apports des tributaires et des processus de régénération dans le maintien de l'eutrophisation Université d'Aix marselle II.

64. Aminot A KR (2004) Hydrologie des écosystèmes marins. Paramètres et analyses. Quae.

65. Holmes MR, Aminot A, Kerouel R, Hooker BA, Peterson JB (1999) A simple and precise method for measuring ammonium in marine and freshwater ecosystems. Can J Fish Aquat Sci 56: 1801–1808.

66. Raimbault P, Garcia N, Cerutti F (2008) Distribution of inorganic and organic nutrients in the South Pacific Ocean – evidence for long-term accumulation of organic matter in nitrogen-depleted waters. Biogeosciences 5: 281–298.

67. Raimbault P, Neveux J LF (2004) Dosage rapide de la chlorophylle a et des phaeopigments a par fluorimétrie après extraction de méthanol. Comparaison avec la méthode classique d'extraction à l'acétone. Oceanis 30: 189–205.

68. Gohin F, Druon JN, Lampert L (2002) A five channel chlorophyll concentration algorithm applied to SeaWiFS data processed by SeaDAS in coastal waters. International Journal of Remote Sensing 23: 1639–1661.

69. Gohin F, Loyer S, Lunven M, Labry C, Froidefond J-M, et al. (2005) Satellite-derived parameters for biological modelling in coastal waters: Illustration over the eastern continental shelf of the Bay of Biscay. Remote Sensing of Environment 95: 29–46.

70. Gohin F (2011) Annual cycles of chlorophyll-a, non-algal suspended particulate matter, and turbidity observed from space and in-situ in coastal waters. Ocean Science 7: 705–732.

71. Shutler JD, Smyth TJ, Saux-Picart S, Wakelin SL, Hyder P, et al. (2011) Evaluating the ability of a hydrodynamic ecosystem model to capture inter- and intra-annual spatial characteristics of chlorophyll-a in the north east Atlantic. Journal of Marine Systems 88: 169–182.

72. Stow C, Jolliff J, McGillicuddy DJ, Doney SC, Allen JI, et al. (2009) Skill assessment for coupled biological/physical models of marine systems. Journal of Marine Systems 76: 4–15.

73. Allen JI, Holt JT, Blackford J, Proctor R (2007) Error quantification of a high-resolution coupled hydrodynamic-ecosystem coastal-ocean model: Part 2. Chlorophyll-a, nutrients and SPM. Journal of Marine Systems 68: 381–404.

74. Allen JI, Somerfield PJ, Gilbert FJ (2007) Quantifying uncertainty in high-resolution coupled hydrodynamic-ecosystem models. Journal of Marine Systems 64: 3–14.

75. Radach G, Moll A (2006) Review of three-dimensional ecological modelling related to the North Sea shelf system—Part 2: model validation and data needs. Oceanography and Marine Biology; An Annual Review 44: 1–60.

76. OSPAR Commission (1998) Report of the ASMO Modelling Workshop on Eutrophication Issues (5–8 November 1996, The Hague, The Netherlands).

77. Lacroix G, Ruddick K, Park Y, Gypens N, Lancelot C (2007) Validation of the 3D biogeochemical model MIRO&CO with field nutrient and phytoplankton data and MERIS-derived surface chlorophyll images. Journal of Marine Systems 64: 66–88.

78. Jolliff JK, Kindle JC, Shulman I, Penta B, Friedrichs MAM, et al. (2009) Summary diagrams for coupled hydrodynamic-ecosystem model skill assessment. Journal of Marine Systems 76: 64–82.

79. Rose Ka, Roth BM, Smith EP (2009) Skill assessment of spatial maps for oceanographic modeling. Journal of Marine Systems 76: 34–48.

80. Pikitch EK, Santora C, Babcock EA, Bakun A, Bonfil R, et al. (2004) Ecosystem-Based Fishery Management. science 305: 346–347.

81. Fulton Ea (2010) Approaches to end-to-end ecosystem models. Journal of Marine Systems 81: 171–183.

82. Travers M, Shin Y-J, Jennings S, Machu E, Huggett Ja, et al. (2009) Two-way coupling versus one-way forcing of plankton and fish models to predict ecosystem changes in the Benguela. Ecological Modelling 220: 3089–3099.

83. Libralato S, Solidoro C (2009) Bridging biogeochemical and food web models for an End-to-End representation of marine ecosystem dynamics: The Venice lagoon case study. Ecological Modelling 220: 2960–2971.

84. Fuchs R, Dupouy C, Douillet P, Caillaud M, Mangin A, et al. (2012) Modelling the impact of a La Niña event on a South West Pacific Lagoon. Marine Pollution Bulletin 64: 1596–1613.

85. Flynn KJ (2008) The importance of the form of the quota curve and control of non-limiting nutrient transport in phytoplankton models. Journal of Plankton Research 30: 423–438.

86. Klausmeier Ca, Litchman E, Levin Sa (2007) A model of flexible uptake of two essential resources. Journal of theoretical biology 246: 278–289.

87. Leonardos N, Geider RJ (2005) Elemental and Biochemical Composition of Rhinomonas Reticulata (Cryptophyta) in Relation To Light and Nitrate-To-Phosphate Supply Ratios1. Journal of Phycology 41: 567–576.

88. Poggiale J-C, Baklouti M, Queguiner B, Kooijman SALM (2010) How far details are important in ecosystem modelling: the case of multi-limiting nutrients in phytoplankton-zooplankton interactions. Philosophical transactions of the Royal Society of London Series B, Biological sciences 365: 3495–3507.
89. Kriest I, Khatiwala S, Oschlies A (2010) Towards an assessment of simple global marine biogeochemical models of different complexity. Progress in Oceanography 86: 337–360.

90. Los FJ, Villars MT, Van der Tol MWM (2008) A 3-dimensional primary production model (BLOOM/GEM) and its applications to the (southern) North Sea (coupled physical–chemical–ecological model). Journal of Marine Systems 74: 259–294.

91. Crout NMJ, Tarsitano D, Wood a T (2009) Is my model too complex? Evaluating model formulation using model reduction. Environmental Modelling & Software 24: 1–7.

92. Wan Z, She J, Maar M, Jonasson L, Baasch-Larsen J (2012) Assessment of a physical-biogeochemical coupled model system for operational service in the Baltic Sea. Ocean Science 8: 683–701.

93. Skogen MD, Moll A (2005) Importance of ocean circulation in ecological modeling: An example from the North Sea. Journal of Marine Systems 57: 289–300.

94. Edwards KP, Barciela R, Butenschön M (2012) Validation of the NEMO-ERSEM operational ecosystem model for the North West European Continental Shelf. Ocean Science 8: 983–1000.

95. Younes W, Bensoussan N, Romano J-C, Arlhac D, Lafont M-G (2003) Seasonal and interannual variations (1996–2000) of the coastal waters east of the Rhone river mouth as indicated by the SORCOM series. Oceanologica Acta 26: 311–321.

96. Tett P, Gilpin L, Svendsen H, Erlandsson CP, Larsson U, et al. (2003) Eutrophication and some European waters of restricted exchange. Continental Shelf Research 23: 1635–1671.