Major fertilization sources and mechanisms for Mediterranean Sea coastal ecosystems

The Mediterranean Sea is considered an oligotrophic basin with a few biological production hotspots. Some of these productive regions are under the influence of large rivers’ discharges and have been described to suffer from eutrophication indications. In order to quantify the impact of the river‐borne inorganic nutrients on the pelagic production levels and bottom oxygen conditions in the northern coast of the Mediterranean Sea, two hindcast simulations using a coupled hydrodynamic‐biogeochemical model are compared. Both simulations cover a 60 yr period and are identical except for the rivers’ conditions, one is performed considering realistic inorganic nutrient levels in the rivers’ waters while the second has no inorganic nutrients coming from the freshwater discharges. Comparing these two simulations we found that inorganic nutrients in the rivers are responsible for a large fraction of the pelagic primary production and for the totality of bottom hypoxic zones in both the Adriatic and Aegean basin. On the contrary, river fertilization is of marginal importance for the coastal ecosystem of the NW Mediterranean Sea, where mesoscale processes and associated vertical transport of nutrients are the main processes responsible for the majority of the production. We, finally, discuss the consequences of these fundamental differences in fertilization for freshwater management and ecosystems status evaluation in the light of different policy options and within the context of future climate change.

The Mediterranean Sea is considered an oligotrophic basin with a few biological production hotspots. Some of these productive regions are under the influence of large rivers' discharges and have been described to suffer from eutrophication indications. In order to quantify the impact of the river-borne inorganic nutrients on the pelagic production levels and bottom oxygen conditions in the northern coast of the Mediterranean Sea, two hindcast simulations using a coupled hydrodynamic-biogeochemical model are compared. Both simulations cover a 60 yr period and are identical except for the rivers' conditions, one is performed considering realistic inorganic nutrient levels in the rivers' waters while the second has no inorganic nutrients coming from the freshwater discharges. Comparing these two simulations we found that inorganic nutrients in the rivers are responsible for a large fraction of the pelagic primary production and for the totality of bottom hypoxic zones in both the Adriatic and Aegean basin. On the contrary, river fertilization is of marginal importance for the coastal ecosystem of the NW Mediterranean Sea, where mesoscale processes and associated vertical transport of nutrients are the main processes responsible for the majority of the production. We, finally, discuss the consequences of these fundamental differences in fertilization for freshwater management and ecosystems status evaluation in the light of different policy options and within the context of future climate change.
The Mediterranean Sea is typically described as an oligotrophic basin with some isolated "hotspots" of enhanced biological productivity (Sournia 1973;Siokou-Frangou et al. 2010). One of the main originators of this idea have been the ocean color images from remote sensing (as the one shown in Fig. 1) which typically show a very blue, oligotrophic sea with "red," productive areas in some specific regions (e.g., Barale and Gade 2008;Uitz et al. 2012).
Basin-wide, around six "productive" regions could be easily isolated from the satellite information (see boxes in Fig.  1). If we move from west to east along the northern coast we find the Alboran Sea, the Gulf of Lion, the Adriatic Sea, and the Aegean Sea as mesotrophic regions. Along the African coast, the Gulf of Gabes and the Nile delta are the most productive zones (e.g., D'Ortenzio and Ribera d'Alcal a 2009).
The reasons for the production detected in each area are diverse. For example, the production in Alboran is related with the water interchange through Gibraltar (Garcia-Gorriz and Carr 1999; Ruiz et al. 2001;Sole et al. 2016), the tidal forcing (S anchez-Garrido et al. 2015) and the meandering Atlantic Jet (Sarhan et al. 2000;Macias et al. 2006). All these mesoscale processes promote nutrient upwelling into the upper water column, enhancing biological production in this sub-basin . Tidal movements, which promote vertical mixing (Lozano and Candela 1995;Sammari et al. 2006), or satellite bias related with water transparency issues (Bosc et al. 2004) have been referred as potential causes of the satellite-detected production in the Gulf of Gabes.
In the rest of the boxes ("NW Mediterranean," "Adriatic," "Aegean," and "Nile," Fig. 1), there are relatively large riverine discharges. It is, thus, assumed that in the last four regions a large fraction of the biological production happening nearby the coastal areas is triggered by allochthonous nutrients supplied by the freshwater flow (e.g., Bosc et al. 2004;Siokou-Frangou et al. 2010).
River-induced fertilization of coastal waters is a natural process, but it has accelerated worldwide in recent times due to human activities (Diaz 2001) leading to eutrophication. Eutrophication and its consequences are, thus, an increasing threat to coastal ecosystems through harmful algal blooms, alteration of plankton community structure, and hypoxia in bottom waters (Cloern 2001). The excess of nutrients delivered to the coastal areas by the riverine outflow is the main cause of eutrophication (Bouwman et al. 2013), but the imbalance between the different macro-nutrients (mainly nitrogen and phosphorus compared to silica) is also a crucial driver for the development of harmful algal blooms (Billen and Garnier 2007).
A major manifestation of eutrophication is, thus, the development of micro-or macro-algal blooms (Heisler et al. 2008) in the river plume and adjacent shelf water (Lohrenz et al. 1997;Dagg and Breed 2003). A fraction of this organic matter (OM) produced in the plume, along with terrestrial OM inputs, settle to the seafloor where it is decomposed by benthic fauna and bacteria using oxygen (Rabouille et al. 2008). When water exchange and mixing are limited, the oxidation of OM may reduce oxygen below 2 mg/L producing a state of hypoxia in the bottom layer (Rabalais et al. 2002). Under these conditions benthic fauna become stressed or die due to the lack of oxygen (Diaz and Rosenberg 1995).
Globally, the number and size of anoxic and hypoxic areas has grown dramatically in recent years (Diaz and Rosenberg 2008), including several regions in the northern coast of the Mediterranean Sea (Rabouille et al. 2008). For this basin, it has been described how the temporal evolution of the Mediterranean Sea overall productivity is linked to freshwater quality (inorganic nutrients content) during the last 50 yr (Macias et al. 2014a) and how important Po discharges are for the Adriatic Sea environment (e.g., Justić et al. 1987;Solidoro et al. 2009). However a direct assessment of the exact importance of inorganic nutrient loads from the rivers for the coastal eutrophication status is still missing.
In Europe, freshwater quality falls under the regulation of the Water Framework Directive (WFD, 2000/60/EC) but this quality has also a clear impact on the status of nearby marine ecosystems located downstream , and which are regulated by a different piece of the European Union (EU) legislation, the Marine Strategy Framework Directive (MSFD, 2008/56/EC). WFD and MSFD have different objectives and methodologies to achieve "Good Environment Status" (GES) but it is quite obvious, given the existing connections, that a strong harmonization between legislations and a holistic approach to ecosystems management are necessary (COM 2014).
Being eutrophication one of the major MSFD descriptors (D5), it concerns the EU Member States the obligation to assess the environmental status of their marine areas and to establish GES in their jurisdictional waters. Thus, quantifying the impacts of river-borne inorganic nutrients on the production and status of nearby marine ecosystems in the northern Mediterranean sea is not only a scientifically relevant issue but also a necessary exercise to understand how and to what extent management options in the rivers' catchment could influence those marine ecosystems downstream.
An integrated modeling approach covering the landocean continuum in the Mediterranean Sea is used in the present contribution to perform this assessment by comparing the results of two different simulations; one performed including nitrate and phosphate in the rivers' waters and the second executed considering no inorganic nutrients in the freshwater inputs. Results of both simulations are then compared in terms of pelagic primary production and seabottom oxygen levels in the vicinity of the estuarine areas.

Material and methods
Here we use an integrated modeling framework (MF) developed at the Joint Research Centre -Ispra of the EU Commission to specifically assess the consequences of different scenarios on the ecosystem status of regional European Seas . This MF includes the main elements of a Regional Earth System Model, i.e., the atmosphere, the hydrological basin and the oceans (Fig. 2). By considering potential scenarios, the MF allows to evaluate how political decisions and program of measures (for example regarding freshwater management) could affect the Good Environmental Status (GES) of the marine ecosystems downstream (e.g., Garcia-Gorriz et al. 2016).
At the core of the MF there is an ocean model coupled to an ecosystem model ( Fig. 2) that has already been shown to provide a reasonable representation of past and present hydrodynamic and biogeochemical conditions of the Mediterranean basin (Macias et al. 2013(Macias et al. , 2014a. In the present work, two long-term model runs covering the period 1959-2013 are performed; the first (named RR) uses the same realistic external forcing as the one presented in Macias et al. (2014a). This includes reanalysis data for the atmospheric conditions and realistic rivers' conditions in terms of water flow and nutrients' loads (from now on, when the term "nutrient" is used in the text we refer to "inorganic nutrients," specifically to nitrate and phosphate). The second run (termed RN) is identical to the first except that no nutrients are imposed in the rivers' runoff. Henceforth, differences between both runs should be attributed only to the effect of river-borne nutrients.
This two model-runs comparison will allow, thus, to evaluate the relative importance of riverine discharges for marine productivity and to examine the impact of freshwater management (in terms of pollution control) on the GES of coastal marine ecosystems in different regions of the

Macias et al. Fertilization of Mediterranean coasts
Mediterranean basin. A detailed description of the used models and of the simulation set-ups follow in the next "Ocean model" and "Run setups" section.

Ocean model
The 3-D General Estuarine Transport Model (GETM) was used to simulate the hydrodynamics in the Mediterranean Sea. GETM solves the three-dimensional hydrostatic equations of motion applying the Boussinesq approximation and the eddy viscosity assumption (Burchard and Bolding 2002). A detailed description of the GETM equations could be found in Stips et al. (2004) and at http://www.getm.eu.
The configuration of the Mediterranean Sea has a horizontal resolution of 5 0 3 5 0 and includes 25 vertical sigmalayers. A third-order Total Variation Diminishing (TVD) numerical scheme is used as recommended by Burchard et al. (2006). ETOPO1 (http://www.ngdc.noaa.gov/mgg/ global/) was used to build the bathymetric grid by averaging depth levels to the corresponding horizontal resolution of the model grid. The salinity and temperature climatologies used as initial conditions at the start of the model integration were obtained from the Mediterranean Data Archeology and Rescue-MEDAR/MEDATLAS database (http://www. ifremer.fr/medar/) while biogeochemical initial and boundary conditions were computed from the World Ocean Atlas database (www.nodc.noaa.gov/OC5/indprod.html). All model runs described below started with exactly the same initial conditions. Boundary conditions at the western entrance of the Strait of Gibraltar were also computed from the same MEDAR/ MEDATLAS dataset imposing monthly climatological vertically explicit values of salinity. No horizontal currents were imposed at the open boundary. With this boundary configuration the circulation through the Strait is established by the internally adjusted baroclinic balance provoked mainly by evaporation and deep-water formation within the basin (Macias et al. 2016b).
The GETM configuration for the Mediterranean Sea is forced at the surface every 6 h by the following atmospheric variables: wind velocity at 10 m, air temperature at 2 m, dewpoint temperature at 2 m, cloud cover and sea level pressure provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). ERA-40 dataset (at 125 km resolution) is used from 1959 to 1978 and ERAin dataset (at 80 km resolution) from 1979 onwards, in both cases atmospheric inputs are interpolated on the Mediterranean Sea GETM grid. Bulk formulae are used to calculate the corresponding relevant heat, mass, and momentum fluxes between atmosphere and ocean (Macias et al. 2013).
GETM is coupled online to the MedERGOM biogeochemical model (Macias et al. 2014a,b) by using the Framework for Aquatic Biogeochemical Models (FABM, https://sourceforge. net/projects/fabm/, Bruggeman and Bolding 2014). MedER-GOM is a modified version of the ERGOM model (Neumann 2000) specifically adapted to represent the conditions of the pelagic ecosystem of the Mediterranean Sea. It has proven useful to describe present (Mac ıas et al. 2014b), past (Macias et al. 2014a), and future  biogeochemical conditions in this semi-enclosed basin. MedERGOM includes three phytoplankton types, three macronutrients (nitrate, ammonia, and phosphate), detritus and dissolved oxygen as main state variables. The distinction between the three phytoplankton types included in the model is based on their role within the planktonic community. Hence, macrophytoplankton can grow very efficiently under nutrient-rich, lightplentiful conditions (but are not dependent on silicate), micro-, nanophytoplankton able to grow in nutrient-limited conditions and with lower light levels, and N-fixing organisms (e.g., cyanobacteria) able to fix molecular nitrogen and which are regulated by salinity (as they are only allow to grow in low-salinity environments). For a more detailed description of the model structure, equations and parameters' values the reader is referred to Mac ıas et al. (2014b). MedERGOM also has a parameterized description of the sediment layer although no diagenetic processes have been included (e.g., Neumann 2000).

Run setups
The two runs are performed for the period 1959-2013, which is the time-span covered by the ECMWF reanalysis data. The unique difference between the reference (or realistic) run (RR) and the "no nutrients" run (RN) is on the nutrients (nitrate and phosphate) concentrations in riverine waters (see below). In both model runs, a continuous small atmospheric input of nitrate, phosphate, and ammonium (equivalent to their climatological mean) is imposed in the entire model domain: nitrate 8.0e 22 mmol/m 2 d and ammonium 4.0e 22 mmol/m 2 d from EMEP (2015); phosphate 1.2e 23 mmol/m 2 d assuming a N : P in the atmospheric deposition 100 (Markaki et al. 2003).
The configuration of the ocean model used for the RR includes 53 rivers discharging along the Mediterranean coast (red stars in Fig. 3 represent the main Mediterranean rivers). Values for river discharges were derived from the Global River Data Center (GRDC, Germany) database. Reducing the numbers of rivers to just 53 underestimates the total freshwater flow into the Mediterranean. These rivers provide a mean annual freshwater flow of 309 km 3 /yr which is lower than the total flow (including all rivers and creeks) computed by a state-of-the-art hydrological model (LisFLOOD Burek et al. 2013) with a mean value of 400 km 3 /yr. However, we need to work with the reduced number of rivers in order to make the modeling system less costly to run in terms of computation time.
On the other hand, inorganic nutrient loads (nitrate and phosphate) of freshwater runoff were obtained from Ludwig et al. (2009)  to get the most comprehensive, monthly varying dataset on freshwater quality throughout Europe. However, this dataset is far from being complete, so we filled the existing gaps in the series by imposing the climatological concentration for the specific river and the concrete month computed from all available data.
No other inputs of freshwater (and nutrients) such as point sources from large coastal cities are included in our model due to the lack of proper long-term datasets for such inputs. Ignoring these point sources in our simulations could make us underestimate the fluxes provided from RR and RN runs described below.

Results
As a first step and in order to validate the capabilities of the MF to represent current-day surface chlorophyll distribution, a comparison between simulation and satellite estimates for the period 1998-2010 is made in the lower panel of Fig. 1. As already described in Mac ıas et al. (2014b) the simulated surface chlorophyll is quite similar to the satellite estimates, with a correlation coefficient of 0.82 and similar standard deviation, although slightly larger in model than in satellite data (Fig. 1).
Also the climatologic surface chlorophyll from RR for the entire hindcast period   (Fig. 3A) compares quite well with the general pattern described for the satellite estimates above (Fig. 1). Enhanced phytoplankton concentrations are found within the Alboran Sea, at the NW Mediterranean region, in the northern part of the Adriatic Sea and within the Aegean.
The model lacks phytoplankton in the Gulf of Gabes region and in the vicinity of the Nile (Fig. 3A) when compared with satellite estimates (Fig. 1). It is, however, unclear if the high surface phytoplankton concentration estimated by the satellite images at the Gulf of Gabes corresponds to real oceanic primary production or if it could be due to biases because of the very shallow depth of this region that "confuses" the algorithm used for chlorophyll computation from the water leaving radiance (Bosc et al. 2004). Also, anthropogenic pollution derived from industrial activities in the area have been proposed as probable sources for the biological production in the area (Ayadi et al. 2015). Actually, the only way in which we can achieve similar levels of biological production with our model for this region (experiment not shown) is by including pollution from point sources at the coast. On the other hand, the water flow of the Nile and its quality are highly uncertain (e.g., Gleick 1991;Conway and Hulme 1993) and still object of debate. In this case, the imposed flow for the Nile in the MF is at the lower end of recent estimates (NBI 2016), so that might be the cause of some of the missing production in this area.
In any case, our main interest lays on European rivers because these catchment areas are the ones being regulated by EU policies and potentially affected by member states' programs of measures. Henceforth, we have decided to focus our comparison on the three regions along the northern Mediterranean coast that share the presence of productivity hotspots and the discharge of main EU rivers; i.e., the NW Med, Adriatic and Aegean "boxes" in Figs. 1, 3.
When the model is run without nutrients in the rivers' waters, there is still relatively large phytoplankton biomass in the NW Mediterranean region while almost all chlorophyll disappears from the Adriatic and Aegean areas (Fig.  3B). These differences are clearly seen in Fig. 3C where up to 95% of the phytoplankton chlorophyll disappears in the northern Adriatic and in the NW Aegean (Thessaloniki Gulf) while the NW Mediterranean region shows a much moderate decrease. In order to observe the differences in greater detail, the simulations RR and RN are compared for each individual region in Figs. 4-6.
The left column of Fig. 4 shows the mean surface chlorophyll concentration (background color), the mean horizontal currents (black arrows) at 10 meters and the bathymetric contours (gray lines) for the NW Mediterranean region in the RR (A), RN (B) and the anomalies RN -RR in % (C). It could be seen that in both runs, high surface chlorophyll values are simulated nearby the two major rivers (magenta arrows) of the region, the Ebro in the south and the Rhone in the north. The maximum differences between RR and RN are roughly 15%, being simulated in the immediate vicinity of both rivers' mouths (Fig. 4C). The same pattern could be described for the surface (10 m depth) integrated primary production rate (right column in Fig. 4). In both model runs there are comparable levels of biological production with maximum differences only reaching 7% nearby the estuarine areas (Fig. 4F). It is also interesting to note that PPR is not distributed as chlorophyll in the innermost part of the Gulf of Lion but it seems to follow the main horizontal current (black arrows) along the continental shelf in an elongated path extending southeastwards.
The same representation as above is made in Fig. 5 for the Adriatic Sea. In this case, however, and as indicated above, when nutrients from the rivers are eliminated (i.e., the RN) there is a very large (almost total) reduction of both chlorophyll (Fig. 5C) and PPR (Fig. 5F) with differences reaching 95% in the northern part of this region. The currents distribution here (black arrows) shows very weak horizontal water movement in the vicinity of the major rivers (magenta arrows), indicating a relatively large residence time of the water masses in those regions.
The case for the Aegean Sea (Fig. 6) is quite similar to the Adriatic one, a very large chlorophyll (Fig. 6C) and PPR (Fig.  6F) reduction (up to 96%) is observed in the vicinity of major rivers (magenta arrows) when the nutrients are eliminated from the rivers' water. The circulation pattern in here does not show a clear picture due the presence of many islands and the irregular topography of the coast (black arrows in Fig. 6).
To obtain a more comprehensive view of the hydrodynamic processes in each region, Fig. 7 shows along with the horizontal currents (green arrows) the maximum absolute vertical velocity (background color) in the upper 100 m of the water column. It is quite clear that in the NW Med region there is a band with large vertical velocities along the main horizontal current path, where maximum PPR values are found (see Fig. 4D,E). Neither in the Adriatic nor in the Aegean Sea large vertical velocities are found in the vicinity of the major rivers or where high biological production is typically simulated in the RR simulation.
However, only looking at pelagic properties in each region will only provide a partial view of the potential eutrophication process linked to the freshwater inputs as   explained in the introduction. The sinking and oxidation of OM created in the river plume by in situ growth or supplied by the flow (i.e., allochthonous OM) could potentially create hypoxic conditions on the lower layers of the water column (Rabouille et al. 2008). Indeed, the three regions of interest in this work have been described as suffering some sort of oxygen depletion problems (e.g., Diaz and Rosenberg 2008). Although the model is not intended to thoroughly represent the oxygen dynamics, as MedERGOM includes dissolved oxygen (DO) among its state variables (Mac ıas et al. 2014b) and the oxidation and remineralization processes of OM are parameterized in the model equations, we can use the model simulations to explore the behavior of DO at the bottom of the water column in the different scenarios.
The most usual form of hypoxia (and other related eutrophication indexes) in coastal benthic systems happens at seasonal scale (Diaz and Rosenberg 2008;Romero et al. 2013), so we have computed the climatological seasonal values (i.e., for each month) of the percentage of area in shallow waters (< 40 m) presenting DO levels below 2 mg O 2 /L (Fig. 8A) and the total absolute area (in km 2 ) occupied by those hypoxic conditions (Fig. 8B). It could be seen that for RR (solid lines, Fig. 8) hypoxic conditions become more frequent at the end of the summer/beginning of fall coinciding with the maximum in stratification conditions in all studied sites. It is also clear that only in the NW Med (dotted black line, Fig.  8) hypoxic conditions are simulated for the RN scenario. For the other two regions no hypoxic areas are detected when nutrients loads are eliminated from the rivers (dotted blue and red lines in Fig. 8). In relative terms, the hypoxic areas occupy similar percentage of shallow (< 40 m) waters in both the Adriatic and Aegean in the RR simulation ( 8-9% at the peak) (Fig. 8A) although the absolute area is much larger in the Adriatic ( 4500 km 2 ) than in the Aegean ( 2400 km 2 ) (Fig. 8B). For the NW Med, even if the relative area occupied by hypoxic bottom waters is relatively large ( 2-4%) (Fig. 8A), in absolute numbers this area is much smaller (maximum 725 km 2 ) (Fig. 8B) than in the other two regions.
We can, finally, have a look at the areas where hypoxic conditions appear during the seasonal maxima identified in Fig. 8A. For the NW Med (Fig. 9A,B) hypoxic areas appear only at the exact mouth of the Rhone (in the RR simulation) and on a coastal fringe north of the Ebro (both in RR and RN). For the Adriatic (Fig. 9C,D) it is clear that the Po discharge region is where hypoxic conditions are simulated for RR while no hypoxic areas are simulated for RN. Similarly, for the Aegean (Fig. 9E,F) only in RR hypoxic regions in the Thessaloniki Gulf are simulated.

Discussion
The total annual freshwater flow into the Mediterranean Sea is not very large ( 309 km 3 /yr, Macias et al. 2016a)  compared with the basin volume, but it has already been described that the interannual variability of primary production in this Sea is related with changes in the quality of riverine waters (Macias et al. 2014a). This land-sea connection is particularly important for some coastal regions in the vicinity of the major (largest) rivers such as the Rhone and Po (Siokou-Frangou et al. 2010) where eutrophication problems could hamper achieving GES by 2020 (Cardoso et al. 2010) and, hence, the provision of services by such ecosystems (Liquete et al. 2016).
From the bioregionalization made by D'Ortenzio and Ribera d'Alcal a (2009) based on surface chlorophyll images, our three investigated boxes (i.e., NW Med, Adriatic and Aegean) where classified as "coastal" provinces driven by river runoff and continental shelf dynamics. In fact, from the pure observation of climatological chlorophyll data (either from satellite or from model simulations) it seems that all of them share the same main characteristics, with the presence of important rivers and relatively elevated biological production (Figs. 1, 3). However, a definitive quantification of the relative importance of rivers discharges on the local production in each region has not been made up to now.
Our numerical simulations allow "experiments" not feasible in the real world such as the no-nutrients scenario imposed in the RN. By performing such simulations we could quantify the relative importance of external (allochthonous) nutrient inputs vs. the local (autochthonous) nutrient supply on the marine productivity. Eliminating the nutrients from the rivers would not greatly impact the overall nutrient budget of the Mediterranean Sea, as this is mainly controlled by the water interchange and mixing  processes that happen at the strait of Gibraltar (e.g., Mac ıas et al. 2007;Huertas et al. 2012). Moreover, the inclusion of atmospheric deposition of inorganic nutrients in our modeling set-up helps to maintain somehow the overall nutrient budget of the basin.
Adriatic and Aegean share a common pattern, as in such regions nutrients borne by the rivers' waters are responsible for a very large fraction of the planktonic primary production (mean 25.13% for the whole Adriatic and 13.6% for the entire Aegean, Table 1) and of the phytoplankton biomass (mean 26.7% for Adriatic and 12.7% for Aegean, Table 1). Although, it is clear that less primary production translates into less available food for higher trophic levels, total removal of nutrients from the rivers (which is something actually not feasible in the real world, and no desirable) might not necessarily lead to negative impacts in terms of total biological production and, hence, in fisheries that are quite important in those regions (e.g., Coll et al. 2012;Piroddi et al. 2015). Indeed high non-consumable and/or toxic biomass, often at the origin of hypoxia, are a greater threat to the ecosystem, with adverse effects for the fishing and tourism industries (e.g., Lancelot et al. 2011).
Indeed, at a more local scale, excessive nutrients inputs from freshwaters sources have been described to provoke eutrophication problems (such as hypoxic areas generation) in the northern Adriatic (e.g., Justic' et al. 1987;Barmawidjaja et al. 1995) and in some coastal gulfs of the Aegean (e.g., Friligos and Zenetos 1988;Theodorou 1996). Our model simulations seem to confirm the crucial role riverine nutrients loads play in the development of coastal bottom hypoxic regions in both the northern part of the Adriatic (Figs. 8, 9, maximum hypoxic coverage 13% of shallow zones) and in the Gulf of Thessaloniki (Figs. 8, 9, maximum hypoxic coverage 11% of shallow zones). The excessive nutrient loads from the rivers and the limited water movements in the discharge areas (green arrows in Fig. 7) could together explain the development of large areas with hypoxic conditions in these two regions of the northern Mediterranean coast.
Indeed, the occurrence of hypoxic conditions in coastal regions connected to large rivers is typically due to the interplay between water column production supplying OM to bottom waters (where it is oxidized), and physical forcings (stratification, the circulation pattern in the region and residence time of the water in the mixing zone) which limit oxygenation of such bottom waters (Rabouille et al. 2008). The presence of such seasonal hypoxia in coastal systems associated with anthropogenic-produced nutrient pollution is one of the major deleterious anthropogenic influences on estuarine and marine environments (Diaz and Rosenberg 2008) and it is clearly an issue in both the northern Adriatic and the coastal Aegean as derived from our simulations.
The same, thought, might not be necessarily true for the NW Med. The effect of river fertilization is there marginally important, both in terms of planktonic production and biomass (mean difference in PPR 0.7% and in chlorophyll a (Chl a) 2.05%, Table 1) and slightly larger for the extension of benthic hypoxic zones in the coastal fringe (mean difference 23.4%, Table 1). Also, the total extension of hypoxic areas is much smaller than in the other two coastal regions (Fig. 8B) which seems to confirm previous works in this particular area (Rabouille et al. 2008). It has been proposed that oceanographic features in the coastal fringe of this region do not favor the development of bottom hypoxia, such as relatively large horizontal currents, which reduces the retention time of bottom waters and favors ventilation, or the presence of a relative weak stratification (see review in Rabouille et al. 2008) that increases the mixing efficiency of vertical water movements (e.g., Howarth et al. 2011). The currents simulated by our modeling system (both horizontal and vertical, Fig. 7A) seem to confirm that the coastal zone of this region is much more dynamic than in the other two investigated basins and, hence, that water residence time in this area should be much shorter, preventing hypoxic conditions development.
Indeed, this particular entire region of the Mediterranean basin has been described as being under the main influence of mesoscale processes associated to the intense currents and fronts typical of the area (e.g., Castellon et al. 1990;Salat 1995;Pinot et al. 2002). Here only a small percentage of nearshore nutrients are related with riverine outflows (Salat et al. 2002), with mesoscale activity exerting a quite large influence on production patterns in the region (e.g., Estrada et al. 1999;Oguz et al. 2015). The same mechanism seems to Table 1. Difference (%) between RR and RN for each of the studied regions (area shown in the different maps) regarding planktonic primary production (PPR), phytoplankton biomass (Chl a), and benthic hypoxic areas.  Fig. 7) in concordance with measurements (e.g., Pinot et al. 1995). They are, however, sensibly lower than the ones simulated in other modeling studies (e.g., Oguz et al. 2015) in which a much higher resolution model set-up was used. This difference in horizontal resolution and our use of monthly averages instead of instantaneous model outputs should contribute to the sensible dissimilarity in maximum vertical speeds described in the two works.
Independently of the possible sub-estimation of their maximum value, the vertical velocities simulated by our model provoke mixing and fertilization of the upper water column, which seems to be the main mechanism maintaining the majority of the planktonic productivity here (in accordance with Oguz et al. 2015). Such vertical movements are typically associated to geostrophic and ageostrophic cross-frontal circulation linked to the meandering of strong horizontal currents in the presence of a density front (e.g., Mahadevan and Tandon 2006) as it also happens in the Alboran Sea (e.g., Oguz et al. 2014), one of the other production 'hotspots' identified in Fig. 1.
However, we must also consider that this region is the only one where hypoxic conditions are simulated when the inorganic nutrients are removed from the rivers (Fig. 8). It is noteworthy that in RN hypoxic areas are located right at the main rivers' mouths ( Fig. 9B) where the freshwater plumes are located. The presence of these plumes increases the stratification of the water column making vertical mixing and, hence, oxygenation of the bottom layers more difficult. This enhanced stratification together with the relative large primary production of the system in this RN simulation (e.g., Fig. 4D) that supplies enough organic matter to be oxidized, creates the conditions for hypoxic areas to develop in the close vicinity of the estuaries.
Henceforth, even without the fertilization sources from the main rivers, in this NW Med region there is abundant primary production and some seasonal benthic hypoxia do occur (Fig. 4, right column and Fig. 8). It is expected that the consequences of periodic hypoxic conditions on the ecosystem are here less severe than in other regions as, this being a natural characteristic of the ecosystem, it should have evolved to be resilient to these conditions (Kemp et al. 2005).
It must be stressed that, in any case, the extension and exact location of hypoxic zones provided by our model simulations could only be considered as approximates given the following limitations of our modeling approach for the very coastal zone. First of all, the horizontal resolution of the model grid does not allow to fully resolve the submesoscale dynamics, which is of crucial importance in coastal regions (e.g., Oguz et al. 2014) and which could greatly change the general regional behavior in terms of eutrophication indexes distribution (e.g., Romero et al. 2013). We are also missing potentially important pollution sources in the form of big coastal cities delivering waste and rain water directly to the sea. Such point sources of pollution have been shown to greatly modify production seasonality and nutrients ratios (e.g., Romero et al. 2014) and to contribute to bottom hypoxic conditions nearby the sewage areas (e.g., Lajaunie-Salla et al. 2017). We also lack the direct input of OM by the river waters, which, in reaching the bottom, could contribute to oxygen depletion. Finally, in our model set-up the sedimentary layer and diagenetic processes are simply parameterized, when it has been recurrently shown that both diagenetic and biologically mediated process in the sediments could greatly affect the bottom oxygen evolution and, even, the nutrients speciation (e.g., Seitzinger 1998;Howarth et al. 2011). To fully resolve these issues and to make a better, more exact description of the bottom oxygen levels in these regions, a dedicated ad-hoc modeling exercise should be made in the future.
Another factor not included in our modeling system is the potential effect that the suspended sediments (SS) load that rivers carry could have on the limitation of primary production in the nearby marine ecosystems. It has been shown that large rivers carrying big loads of suspended solids can create quite large turbid plumes at their mouths (e.g., Caballero et al. 2014) that effectively limit primary production and nutrients' assimilation in the area (e.g., Ruiz et al. 2013). The same is true for the rivers considered in the present work as they have been described to carry larger amounts of suspended solids than equivalent (in terms of water flow) northern European rivers (e.g., Romero et al. 2013), including two major contributors such as the Rhone (SS 16 mg/L; Higueras et al. 2014) and the Po (SS 38 mg/L; Camusso et al. 1996). The presence of these high SS plumes will like make that the nutrients supplied by the riverine waters could not be used in the vicinity of the estuaries, but only after the sediments have started to sink to the bottom (e.g., Prieto et al. 2009). Not including this factor in our model could, hence, contribute to some of the mismatches between satellite-observed and model-simulated production patterns nearby the rivers' mouths.
In any case, what our modeling exercise has clearly shown is that nutrient reduction in European rivers seems to be an effective measure to control excessive production, hypoxic conditions generation and other eutrophication-related issues in both the Aegean and Adriatic basins. These nutrient-control measures fall under the WFD which impose to the EU member states the obligation to reach the GES (Ferreira et al. 2011) and the "Good Chemical Status" of their freshwater bodies. However the measures to be taken for each different catchment area need to be carefully evaluated, ideally following a cost-benefit analysis (e.g., Lancelot et al. 2011) in order to enhance the potential benefits. Furthermore, caution should be applied so that the nutrient balance (nitrate and phosphate to silicate) is not much changed. If so happens, the potential for community shifts in the estuary and coastal areas is large (e.g., Passy et al. 2013) with a further proliferation of toxic plankton blooms being more probable (Howarth et al. 2011). This has happened in other estuarine areas within Europe (e.g., Romero et al. 2013;Passy et al. 2016) with consequences for the regions under the regulation of the MSFD, which extended to all European seas the need to achieve or maintain GES in the marine environment by the year 2020.
It is evident that such freshwater pollution control measures will have a much smaller impact on the NW Mediterranean marine ecosystem than in the other two studied areas, given the relatively smaller importance of river-borne nutrients for this particular system. However, potential changes in the currents configuration and on the vertical water stratification caused by future climate (e.g., Macias et al. 2015) will have, potentially, a relatively larger impact on the production levels and ecosystems status of this NW region of the Mediterranean Sea.

Conclusion
We have, thus, evaluated within the limits imposed by the modeling system (horizontal resolution, lack of point sources, neglecting the organic matter contribution, etc.), the relative importance of river fertilization for three important Mediterranean coastal ecosystems in terms of pelagic production levels and on hypoxic conditions in the bottom water layer. It is quite clear that river-borne nutrients are much more relevant in the semi-enclosed Adriatic and Aegean Seas than in the more open NW Mediterranean region. A total removal of nutrients from the rivers' waters in the two former areas (not fully achievable in reality) will almost completely eliminate pelagic primary production in those regions, which could largely affect also secondary productivity (Nixon and Buckley 2002). On the other hand, controlling the nutrients carried by freshwater flow in these regions could be an effective way of dealing with the presence of seasonal hypoxic conditions in shallow coastal settings. Managing freshwater quality for the main rivers discharging in the NW Mediterranean region will, however, have less important consequences for their marine ecosystems as most of the primary production here is driven by internal, regional oceanographic processes. We need to be careful, however, when numerically evaluating the exact extension and location of the areas affected by eutrophication as simulated by the model, given that its horizontal resolution does not allows for a complete resolution of very coastal processes. However, the quantification and attribution exercise here presented could be a useful tool to better direct management options for EU rivers and is a clear example on how the MF could be used to evaluate different policy options.