Pelagic primary production in the coastal Mediterranean Sea: variability, trends and contribution to basin scale budgets

. We estimated pelagic primary production (PP) in the coastal (<200 m depth) Mediterranean Sea from satellite-borne data, its contribution to basin-scale carbon fixation, its variability and long-term trends during the period 2002-2016. Annual coastal PP was estimated at 0.041 Gt C, which approximately represents 12 % of total carbon fixation in the Mediterranean 15 Sea. About 50 % of this production occurs in the eastern basin, whereas the western and Adriatic shelves contribute with 25% each of total coastal production. Strong regional variability is revealed, from high-production areas (>300 g C m -2 ) associated with major river discharges, to less productive provinces (<50 g C m -2 ) located in the southeastern Mediterranean. PP variability in the Mediterranean Sea is dominated by interannual variations but overall trend during the study period shows notable decrease (17%) since 2012 concurring with a period of increasing sea surface temperatures in the Mediterranean Sea 20 and positive North Atlantic Oscillation and the Mediterranean Oscillation climate indices. PP declines in most coastal areas (-0.05 to -0.1 g C m -2 per decade) except in the Adriatic where PP increases at +0.1 g C m -2 per decade. Regionalization of coastal waters based on PP seasonal patterns reveals the importance of river effluents in determining PP in coastal waters that can regionally increase in up to five-fold. Our study provides insight on the contribution of coastal waters to basin scale carbon balances in the Mediterranean Sea while highlighting the importance of the different temporal and spatial scales of variability. patterns. Two main PP groups were observed: zones with strong cross-shore gradients, typically found in wider estuarine regions and homogeneous zones within narrow continental shelf areas.

Administration (NASA) archive website (http://oceancolor.gsfc.nasa.gov/), were used for model calculations. Only night-time orbits were selected to avoid problems with skin temperature during daylight. Orbits with quality flag 2 in SST were included after checking their validity and accuracy in order to have a more complete dataset. Daily (24-hour averaged) 100 Photosynthetically Active Radiation (PAR, in E m -2 ) was obtained as a Level 3 product at 9 km, the best available resolution from the NASA archive from both MODIS and Medium Resolution Imaging Spectrometer (MERIS).
All satellite-derived variables were remapped onto a regular 1-km spatial grid over the study area, by averaging all available pixels within each grid cell. For each parameter, outliers were removed whenever they exceeded about 3-times the mean ± SD of the time series. For the purpose of this study, coastal areas were defined as the waters lying between 5 and 200 105 m depth. Only values at depths exceeding 5 m depth were considered in order to avoid any chlorophyll (Chl,   Bathymetric data were obtained from ETOPO1 (Amante and Eakins, 2009). The black contour indicates the 200m isobath, the limit of coastal waters as defined in the present study.

Primary production estimates
PP was estimated from satellite-derived Chl, SST and PAR values using the time, depth, and wavelength-resolved lightphotosynthesis model of Morel (1991). This model was previously used for estimating PP in the Mediterranean sea (Antoine 115 and André, 1995) and at global scale  and performs well when compared to in situ measurements (Campbell et al., 2002;Friedrichs et al., 2009) or when compared to other similar algorithms designed for https://doi.org/10.5194/bg-2020-457 Preprint. Discussion started: 11 January 2021 c Author(s) 2021. CC BY 4.0 License. use with satellite observations (Carr et al., 2006;Saba et al., 2011). Instantaneous production at depth z (m) of the water column, time t of the day, and for absorption of irradiance at wavelength , P(, z, t), is calculated as: 120 P (, z, t) = E (, z, t) Chl (z) a * (, z)  (mol C m -3 s -1 ), where E() is the spectral scalar irradiance (mol photons m -2 s -1 ), a * () is the spectral chlorophyll-specific absorption coefficient of phytoplankton (m 2 mg Chl -1 ), and  is the quantum yield of photosynthesis for carbon fixation (mol C mol photons -1 ; its possible spectral changes are ignored). Note that neither Chl, a* and  are made variable with time. 125 The triple integration of (1) w.r.t. wavelength, depth and time gives the daily column-integrated primary production, PP: = 12 ∫ ∫ ∫ P(λ, z, t) dλ dz dt 700 400 min ( / ) 0 0 (g C m -2 ), (2) 130 where the factor 12 is the conversion from moles to grams of carbon, D is the day length (h), Zp (m) is the depth where the photosynthetically available radiation (PAR) falls to 0.1% of its value just below the sea surface (so approximately 1.5 times the euphotic depth), and Zb is the bottom depth taken from the ETOPO1 data base (Amante and Eakins, 2009). The time integration used intervals equal to 1/30 of the day length (about 20 to 30 min depending on season). The depth integration used 135 intervals equal to 1/50 of Zp and goes down to whichever is shallower between Zp and Zb. The spectral integration is performed over the visible range (400 to 700 nm) with a 5 nm resolution.
The spectral irradiance at a given depth z, E (, z, t), is calculated as (starting from just below the sea surface): 140 E(λ, z, t)=E(λ, z − dz, t)e [− (λ,z) The model was operated both for clear sky conditions and for the actual MODIS PAR values, in which case a 150 reduction of the clear-sky irradiance is uniformly applied across the entire day, as being the ratio of the satellite to clear-sky PAR daily values. Chl is assumed to be uniformly distributed with depth, and equal to the satellite-derived value. This simplification was considered more appropriate for the generally shallow and well-mixed waters of coastal areas than the use of global parameterization of the shape of the vertical profile as a function of the surface Chl value (e.g. Morel and Berthon, 1989;Uitz et al., 2006), whose validity outside of open ocean waters is not established. 155 From PP estimates, new (PPnew) and regenerated (PPreg) production were calculated using the ratio of export production to total production (i.e., ef-ratio) (Laws et al., 2000;2011). Indeed, assuming a steady state, the export production must equal the new production fuelled by new nutrients brought to the surface layers. The ef-ratio as a function of satellitederived temperature and production can be obtained from the empirical relationship obtained by Laws et al. (2011): where T is temperature in degrees Celsius (°C) and PP is the daily production (mg C m -2 ).
We report annual PP estimates (Gt C) for the entire Mediterranean coastal areas (ΣPPcoast) and separately for the Western, Eastern and Adriatic basins (ΣPPbasin). While some authors include the Adriatic in the eastern basin (e.g. Bosc et al., 170 2004), we treated it separately because its peculiarities (i.e. bathymetry, influence of rivers, eutrophic character) differentiate it from the rest of the Mediterranean Sea (Cushman-Roisin et al., 2001). Most of the Adriatic has a shallow (<200 m) bathymetry and it collects some 30% of the freshwater flowing into the Mediterranean, acting as a dilution basin for the nutrients discharged by the Po and other Adriatic rivers and becoming one of the most human-impacted regions of the Mediterranean Sea (Ludwig et al., 2009;Micheli et al., 2013;Raicich et al., 2013). 175

Coastal regionalization
We used a two-step classification procedure to define coastal regions along the Mediterranean based on their temporal PP patterns. First, 9 regions (R1 to R9) were identified using a classification technique based on an unsupervised learning neural network (Self-Organizing Maps or SOM; Kohonen, 1982;2001). Then, 18 alongshore marine ecoregions were obtained considering the most relevant cross-shore limits of the SOM-derived regions (Z1 to Z18).
SOM is an unsupervised neural network method that reduces the high dimensional feature space of the input data to a lower dimensional network of units called neurons. SOM is especially suited to extract patterns in large datasets of satellite data (Ben Mustapha et al., 2014;Charantonis et al., 2015;Farikou et al., 2015). Unlike other classification methods, like kmeans, SOM tends to preserve data topology (i.e. preserves neighbouring regions) and, therefore, it is particularly suited for pattern recognition (Liu and Weisberg, 2005). It allows adequate classification of areas with high spatial complexity and strong 185 gradients. Similar neurons are mapped close together on the network facilitating the visualization of patterns and a topological ordination of the classified areas and the relative distance among neurons is obtained as results of the analysis.
For typical satellite imagery, SOM can be applied to both space and time domains. Here, we have addressed the analysis in the time domain of the datasets, which allows regionalizing the studied area on the basis of similitudes in the time variation of PP. We chose a map size of (3 x 3), with 9 neurons (for further details, see Basterretxea et al. 2018). We used a 190 hexagonal map lattice in order to have equidistant neighbours and to avoid introducing anisotropy artefacts. For the algorithm initialization, we opted for linear mode, batch training algorithm, and 'ep' type neighbourhood function since this parameter configuration produces the lower quantitative and topological error and computational cost (Liu et al., 2006). These SOM computations were performed using the MATLAB toolbox of SOM v.2.0 (Vesanto et al., 2000a(Vesanto et al., , 2000b provided by the Helsinki University of Technology (http://www.cis.hut.fi/somtoolbox/). 195

Climate data
To identify possible drivers of long-term PP variability we searched for correlations with two climate indices, the North Atlantic Oscillation index (NAO) and the Mediterranean Oscillation Index (MOI). The corresponding data were downloaded from the Climate Research Unit at the University of East Anglia (https://crudata.uea.ac.uk/cru/data/). Climate indices are defined either as anomalies of a climate variable, using the difference between two geographical points, or as principal 200 components (Hurrell, 1995).
NAO is the central mode of climate variability of the Northern Hemisphere atmosphere. It is based on the pressure difference between the middle of the North Atlantic Ocean and Iceland, which affects winter conditions in the North Hemisphere (Hurrell and Van Loon, 1997;Marshall et al., 2001). Positive NAO results in a relatively dry winter in the Mediterranean but a warmer and wetter winter in northern Europe, and vice versa. Because of its influence on precipitation, 205 Mediterranean river inflows are generally anti-correlated with the NAO (Trigo et al., 2006). MOI is the most widely used teleconnection index for the Mediterranean basin. It reflects differences in temperature, precipitation, circulation, evaporation and other parameters between the eastern and western basin. There are different versions depending on the points of reference (Criado-Aldeanueva and Soto-Navarro, 2013). We used the version obtained as the normalized pressure difference between Gibraltar and Israel (Palutikof, 2003). Positive MOI phases are associated with 210 increased atmospheric pressure over the Mediterranean Sea that promotes a shift of the wind trajectories toward lower latitudes leading to milder winters (Criado-Aldeanueva and Soto-Navarro, 2013). Under these conditions, reduced precipitation is observed in the southeastern Mediterranean region (Törnros, 2013). With some regional differences, NAO and MOI express https://doi.org/10.5194/bg-2020-457 Preprint. Discussion started: 11 January 2021 c Author(s) 2021. CC BY 4.0 License. relatively similar climate patterns over the Mediterranean Sea. They are highly positively correlated in winter, and weakly but still significantly correlated in summer (Efthymiadis et al., 2011;Martínez-Asensio et al., 2014). 215

Statistical analyses
Linear temporal trends in the PP series were calculated using Theil-Sen slope adjustment (Sen, 1968) of the residuals of the deseasonalized series. Only pixels with a trend statistically significant at the 95% level were considered. Correlation analyses were performed using the Pearson Product Moment correlation. Differences between means were tested using the Kolmogorov-Smirnov test (Massey et al., 1951). 220

Coastal primary production
Annual primary production in coastal waters of the Mediterranean Sea (ΣPPCoast) is estimated to be 0.041±0.004 Gt C, which represents some 12% of total carbon fixation in the Mediterranean Sea (see Table 1 and 2). Approximately, 80% of this ΣPPCoast is sustained by recycling processes and, the rest, PPnew, is exported to the seafloor or to nearby areas. Although average surface 225 Chl concentration is 3-fold higher in coastal areas (0.3 mg m -3 ) than in open areas (0.11 mg m -3 ), the annual carbon fixation per surface area (PPannual) over the shelf is, on average, 26% lower than in the open ocean (100±91 and 136±40 g C m -2 respectively; see Table 1 and 2). We would have expected that Mediterranean coastal annual PP would also be higher (about 1.7-fold) than oceanic annual PP. This hypothesis would have been observed if depth integration in coastal areas would have been down to Zp. However, depth integration in coastal areas is quite often stopped at a much shallower depth, i.e. Zb, hence 230 the lower PP per unit area. The surface volumetric PP has been also estimated with a mean value of 2.93±9.60 g C m -3 (Table   1). This mean value of volumetric productivity could not have been compared with oceanic productivity in the Mediterranean Sea considering that previous works presented only integrated production (e.g. Azov, 1986;Bosc et al., 2004;Bricaud et al., 2002a;Coll et al., 2010;Krom et al., 1991;Macias et al., 2015).
Some differences in PP are observed between the more productive shelf waters in the western basin and those in the 235 eastern basin (PPannual =98±55 g C m -2 and 92±96 g C m -2 , p<0.001), the Adriatic shelf being by far the most productive (PPannual= 123±106 g C m -2 , Table 1). Annual carbon fixation is 97% higher in the eastern (ΣPPeast=0.021±0.002 Gt C) than in the western shelf (ΣPPwest=0.011±0.001 Gt C), which is due to greater eastern surface area (about twice the western surface area; Table 1). PPannual varies spatially from 90 to 250 g C m -2 in the western shelf, and from 50 to 400 g C m -2 in the eastern basin where lowest values (<75 g C m -2 ) are found mainly along the Gulf of Sirte. Contrastingly, PPannual exceeds 100 g C m -2 240 in the Adriatic basin reaching values above 400 g C m -2 in the north western coast (Fig. 2b). The most productive coastal regions (>150 g C m -2 ) are mainly located along the European coasts and seem to be related with the outflow regions of major rivers. Indeed, the highest values of the coefficient of variation of primary production (CVPP) are observed in the mouth of the Ebro, Rhone, Tiber, Po, Neretva or Nestos/Evros rivers. Along the western African coast, PPannual displays values >150 g C m -https://doi.org/10.5194/bg-2020-457 Preprint. Discussion started: 11 January 2021 c Author(s) 2021. CC BY 4.0 License. 2 ; however, since the shelf is narrow, its contribution to ΣPPCoast is marginal ( Fig. 1 and Fig. 2). The annual volumetric 245 productivity follows a similar pattern than the annual integrated production with most values varying between 1-3 g C m -3 but reaching up to > 10 g C m -3 in the most productive coastal regions of the Adriatic Sea and the Gulf of Gabes (Fig. 3a). Surface area, annual mean surface Chl, annual average PP (PPannual), annual integrated PP (ΣPP) and annual average 250 productivity per unit volume PP (PPVOLannual).  Western basin 1980-1985 120 In situ (Oxygen) (Bethoux, 1989(Bethoux, ) 1981 157,7 Satellite (CZCS) (Morel and André, 1991b) 1979-1983 157 -197 a Satellite (CZCS) André, 1995) 1996 140 -150 In situ ( 14 C data) (Conan et al., 1998(Conan et al., ) 1997(Conan et al., -1998 Satellite (SeaWiFS) (Bricaud et al., 2002a(Bricaud et al., ) 1991(Bricaud et al., -1999 In situ ( 14 C method) (Marty et al., ) 1996 175

Long-term variability and trends
As shown in Fig. 4, variability in annual PP is dominated by short-scale variations (i.e. subdecadal). ΣPPCoastal exhibits moderate interannual variability (up to 25%) whereas basin scale interannual variations range from 26% in the Adriatic basin, up to 28% 290 in the western basin and 29% in the eastern basin. When considering the whole basin, positive anomalies in coastal PP extended between 2004 and 2011 (mean 0.044±0.001 Gt C y -1 ; Fig. 4a). Conversely, year 2012 was particularly unproductive in all three basins (specific annual mean PP for 2012 were 0.037 Gt C y -1 for the whole basin, 0.010 Gt C y -1 for the western,0.019 Gt C y -1 for the eastern and 0.009 Gt C y -1 for the Adriatic basin). This negative anomaly marked the beginning of a less productive period, particularly noticeable in the eastern basin (Fig. 4c). 295 Long-term trends in PP at 95% of confidence level are significant at basin scale and also in the western and in the eastern basins (p<0.05; Fig. 4a-c) whereas the Adriatic Sea does not display significant trends. However, while a slight negative tendency in seen in the western coast (-10.73 TC per decade. Fig. 4b), a more dramatic tendency, driven by a shift in year 2012, is observed in the eastern basin (-25.39 TC per decade; Fig. 4c). As revealed by Fig. 5a, some regionally coherent patches 305 of significant trend in PP are observed along the coast. Most of these regions presented declining PP trends, particularly along the north African coast where SST temperature is increasing at a higher rate (Fig. 5b). Typical PP trend magnitudes observed along the Spanish Mediterranean and the North African coast from the Gulf of Gabes range from -0.05 to +0.05 g C m -2 per decade. Some positive PP trends, exceeding +0.1 g C m -2 per decade, can be determined in some coastal regions of the north of the Adriatic Sea. 310

Coastal regionalization
The nine characteristic temporal PP patterns, their corresponding spatial distribution obtained from the SOM analysis and the 18 zones in which the coastal region was classified are shown in Fig. 7. Generally, wider shelves present higher spatial complexity manifested as a larger number of SOM patterns. About 78% of the shelf waters include R1, R2, R3 and R4 patterns.
In particular, R1 and R2 characterize areas of low production with scarce seasonality (PPannual=44±17 and 69±22 g C m -2 , 330 respectively; Table 3) typically occurring in the southern and eastern Mediterranean (12 and 29% of the total surface area).
They are representative of the productivity patterns in vast shelf regions in the Gulf of Gabes and Sirte and in the central Aegean (Z17, Z16, Z14 and Z11). R3 and R4 correspond to higher production and a wider range of variation (PPannual=90±32 and 98±39 g C m -2 ; Table 3    As shown in Table 3, areas of low production with seasonal patterns R1 to R4 contribute to more than 62% of total pelagic carbon fixation in Mediterranean coastal areas. In contrast, systems of higher production (PPannual>280 g C m -2 d -1 ) barely contribute to 17% of total production. These regions of enhanced production are generally constrained to Regions of 355 Freshwater Influence (ROFIs; Simpson, 1997) where terrestrial nutrients fuel coastal production. Indeed, R8 is almost exclusively restricted to the river mouths and it presents elevated PP values (1.29±0.50 gC m -2 d -1 ) and a wide range of variation of 0.67 to 2.14 gC m -2 d -1 . R7 pattern is exclusively located in the shallowest inner shelf of the Gulf of Gabes and it is bounded by R5, a transition region between the inner and outer shelf. Unlike the other regions, where PP peaks in late winter-spring, maximum PP in R7 occurs in fall. Finally, R5 and R6 patterns correspond to transition regions accounting for 20.4% of the 360 total production and 14.8% of the Mediterranean coast. While R5 mainly occurs in deltas areas, R6 is characteristic of the western Mediterranean shelf, including the North African coast (0.36-0.08 gC m -2 d -1 ; Fig. 7).
The SOM-based regionalization reveals two groups of coastal waters: those with low cross-shore variability and including only one or two SOM regions (i.e. Z1, Z13, Z16 and Z18) and those with strong cross-shore gradients including several SOM-regions (i.e. Z4, Z9, Z12, Z15). The first pattern is typically observed in narrow continental shelf areas with low 365 influence of river inputs whereas the second group is found in regions with wider continental shelf such as ROFIs (the Rhone delta, the north and western coastline of the Adriatic Sea and the Nile Delta) and in the Gulf of Gabes. The western Adriatic (Z9) and the Gulf of Gabes (Z17) are the largest contributors to ΣPPCoast, contributing together to 35.9% of shelf production in the Mediterranean Sea but, in the case of Z17, it is mainly due to its large extension (Table 4). PP is also high in the northern Alboran Sea (Z1), Nile delta (Z15), the western Adriatic (Z9), and Gulf of Lions (Z4; Table 4). With the exception of Z1, 370 influenced by the entrance of waters from the Atlantic Ocean and by local coastal upwelling, these zones receive important riverine fluxes (Q).

Coastal primary production
To our knowledge, this is the first study focused on the contribution of coastal waters to the overall pelagic PP in the Mediterranean Sea. While the mean coastal values for the Mediterranean (100±91 g C m -2 ) are somewhat lower that the mean 380 values over the continental shelves of the World ocean (160±40 g C m -2 ; Smith and Hollibaugh, 1993), the impact of coastal pelagic PP to total basin production (12%) is in the high range of the estimations for other Seas (Muller-Karger et al., 2005).
This estimation is subject to the uncertainties inherent to using satellite ocean colour, which is limited to the upper ocean (down to 20 m at best in clear waters) and has poor performance in some areas (i.e. Case-2 waters). It nevertheless provides an assessment of net rates of carbon fixation in coastal areas that is consistent with global estimations of the contribution of coastal 385 areas to oceanic production (Gattuso et al., 1998;Ducklow and McCallister, 2004;Muller-Karger et al., 2005). Bias in coastal Chl estimations is mainly due to the presence of non-phytoplankton components such as coloured dissolved organic matter (CDOM) or other terrestrial substances . These compounds originate from coastal erosion, resuspension in shallow areas, river inputs or anthropogenic effluents. Likewise, they affect the propagation of photosynthetic radiation through the water column (Morel, 1991). However, the possible uncertainties and biases caused by Chl estimation through satellite data might have alter very weakly our estimation of coastal PP. Indeed, Case-1 waters are largely predominant in the coastal Mediterranean regions whereas Case-2 waters are reduced to less than 5% of the whole basin. In particular, they are confined to the north Adriatic Sea, Gulf of Gabes and around Nile delta (Antoine et al., 1995;Morel and André, 1991;Bosc et al., 2004) where our PP estimations may present larger uncertainties. However, PPannual values off the Nile river delta, >100 g C m -2 estimated here, are only slightly higher than those reported by Antoine et al. (1995) (80-100 g C m -2 ). Highest values have 395 been reported for this region (>300 g C m -2 ) but, as shown in Fig. 2, they are restricted to a narrow coastal band. In the case of the Adriatic Sea, Umani (1996) reported values of PP from 50 to 200 g C m -2 y -1 , while Zoppini et al. (1995) estimated PP rates from 210 to 260 g C m -2 y -1 in the northern coastal areas. Our estimations range between 100 and >350 (with mean values of 123±106 g C m -2 ).
Because of its extension, the eastern basin contributes more than the western basin to overall coastal production (50% 400 and 25% respectively; Table 1). Furthermore, carbon fixation from the Adriatic Sea represents 24% of total coastal production, which is significant considering the area of this Sea (19% of Mediterranean coastal waters). The relevance of the contribution of the Adriatic Sea in overall coastal PP relies in two main characteristics; (1) coastal waters (<200 m) constitute a large part of the Adriatic Sea and (2) about one third of the river discharge in the Mediterranean is concentrated in the Adriatic Sea (see Table 4). Indeed, patterns in the northern Adriatic Sea reflect a variation in the drivers of PP with respect to other regions. For 405 example, while internal processes (i.e. vertical diffusion and mixing) and, less so, atmospheric deposition, drive PP in most coastal waters, production in the north Adriatic would be mainly driven by fluvial sources of carbon and regeneration through bacterial pathways (Umani et al., 2007). Moreover, distinctive dynamics in this sea is driven by the influence of river outflows on stratification and general circulation patterns Giani et al., 2012).
Here, from ef-ratios, we estimated that on average only 22±20 % of the production in the coastal Mediterranean Sea 410 is new and the rest is sustained via regenerated sources (Fig. 8a-c). This PPnew value (Fig. 8b) is comparable to the mean organic carbon that sinks to the sea floor (28%) estimated from Muller-Karger et al. (2005) and Pace et al. (1987) but higher than PPnew estimations provided by Vidussi et al. (2001) for oceanic waters in the eastern basin (15% of total production).
Contrarily to Vidussi et al. (2001) who estimated PPnew in the eastern basin, the coastal PPnew average here includes both eastern and western basins of the Mediterranean, but also the highly productive areas in the northern Adriatic. This could explain that 415 the PPnew observed here is higher than the one observed for oceanic waters in the eastern basins. Additionally, high ef-ratios (> 0.3) are observed in our case in the areas where nutrient inputs from the Atlantic and river effluents significantly enhance PPnew (Fig. 8a). Furthermore, ef-ratios present significant seasonality, varying between 0.26±0.04 in the most productive winter-spring season and 0.15±0.02 in summer, when the water column is strongly stratified and food web shifts to a more recycling dominated system. 420 https://doi.org/10.5194/bg-2020-457 Preprint. Discussion started: 11 January 2021 c Author(s) 2021. CC BY 4.0 License.

Long-term variability and trends
Available satellite ocean colour data span about 20 years, so that temporal trends derived from their analysis are highly 425 depending on decadal variability (Henson et al., 2010). Despite these limitations, satellite observations of ocean colour over the past two decades suggest relationship between warming and reduced productivity in permanently stratified areas (Behrenfeld et al., 2006). Since clear tendencies of warming are observed in the Mediterranean Sea (Nykjaer, 2009;Pastor et al., 2018), intensification of stratification would decrease nutrient supply to phytoplankton and, thus, decrease PP (Behrenfeld et al., 2006;Stambler, 2014). 430 Barale et al. (2008) observed a general decrease in Chl biomass in the Mediterranean Sea over the period 1998-2003.
However, some coastal areas in their study displayed the opposite tendency. Macias et al. (2015) anticipated no future global changes of integrated PP in the Mediterranean Sea from modelling results. They predicted a tendency to oligotrophication in the western basin and increase in the productivity of the eastern basin. Our study reveals that ΣPPCoastal in the Mediterranean Sea varies nonlinearly and a reduction of carbon fixation rates is observed since 2012 (Fig. 4a). Overall negative trends are 435 reported here in both the Western and in the Eastern basin (-10.70 and -25.39 TC per decade; Fig. 4b-c). A spatial analysis of the long-term decadal variability reveals weak but spatially coherent and significant tendencies (p<0.05; Fig. 5 Font et al., 2007;Schroeder et al., 2008;Šolić et al., 2008;Viličić et al., 2012). 445 Long-term decadal variations in the eastern and western basins are mostly coupled, suggesting that they share the same PP drivers at this basin scale (Fig. 4b-c). A major feature in the interannual pattern is a global decrease in production in 2012 that is extended to the following years in the eastern basin. Durrieu de Madron et al. (2013) reported peculiar atmospheric conditions in the Mediterranean Sea during 2012 that triggered a massive formation of dense water on the continental shelf and in the deep basin of the Gulf of Lions. A similar anomaly was described in the Adriatic shelf where unprecedented dense 450 water generation was preconditioned by a dry and warm year resulting in a significant reduction of coastal freshwaters and basin-wide salinity increase (Mihanović et al., 2013;Raicich et al., 2013). Additionally, Pastor et al. (2018) observed an anomalously temperature increase in the Mediterranean Sea during summer 2012. From our analysis, we infer that this climaterelated event had strong influence on the global coastal PP of the Mediterranean Sea.
Several studies have reported influence of climate variations in the coast (Belgrano et al., 2008;Cloern et al., 2007;455 Tiselius et al., 2016). In agreement, we observed an influence of climate scale variability on coastal productivity as suggested by the inverse correlations between ΣPP and SST and, more loosely, with NAO and MOI (Fig. 6). While these correlations emphasize the pre-eminent role of climate variability in the regulation of interannual to decadal scale coastal productivity, the pathways through which this control of the atmosphere over coastal productivity is exerted are complex and may regionally differ (Grbec et al., 2009). Climate can influence phytoplankton growth by the direct effect of temperature on algal metabolism, 460 by changes in basin scale circulation (including exchanges with adjacent seas), by regulating nutrient supply through variations in the thermocline intensity, by changes in wind patterns affecting mixing and dust deposition pathways or through changes in precipitation that have direct influence on wet deposition and on river runoff. These effects are modulated by changes in the biota and in the interaction between organisms (e.g. Molinero et al., 2005). The relative importance of climate-driven processes relative to other productivity enhancing processes depends on regional characteristics and may be seasonally varying. For 465 example, variations in dust deposition, which may sustain up to 50% of new production in the Levantine basin ( Kress and Herut, 2001;Herut et al., 2002), are expected to be more important in the eastern and southern Mediterranean coasts because of their proximity to the Saharan dust sources. Likewise, variations in cooling and vertical mixing are expected to be more effective during late winter when PP peaks and when diatoms dominate in the Mediterranean Sea (Lacroix and Nival, 1998;Marty, 2002;Marty and Chiavérini, 2010). 470 Our results reveal that, in contrast to other regions like the North Sea (Capuzzo et al., 2018) or the Arctic Ocean (Gregg et al., 2003), the coastal Mediterranean Sea did not globally display a marked decline in PP during the last decades.
We suggest that in some coastal areas, a decrease in vertical nutrient supply though the thermocline may be compensated by other nutrient sources. Variations in atmospheric deposition, groundwater and river outflows together with the influence of human activities through changes in landscape use and nutrient management are important sources of nutrient in the ecosystem 475 and thus, act as major drivers of PP in these waters (e.g. Paerl et al., 1999). As a consequence of human activities, both terrestrial and coastal ecosystems have experienced progressive nutrient enrichment (Conley et al., 2009;Deegan et al., 2012).
However, while this effect is evidenced in shallow nearshore waters, its influence in the ocean is estimated to be minimal (Wang et al., 2018). In the Mediterranean Sea, high coast population growth rates and concomitant food demand have resulted in dramatic increase of water demand for irrigation farming and fertilizer use (Ryan, 2008). Indeed, while the freshwater 480 discharge of Mediterranean rivers has significantly reduced during recent decades (~20%), the corresponding total nitrogen inputs to coastal seas are estimated to have increased by a factor up to 5, fuelling PP in river influenced areas (Ludwig et al., 2009). While the importance of groundwater in the Mediterranean Sea could be comparable to that of rivers (Rodellas et al., 2015) and generalized nitrification of Mediterranean coastal aquifers is acknowledged (EEA, 1995;Zalidis et al., 2002), general trends in groundwater discharges remain largely unknown. 485

Coastal regionalization
Coastal regionalization reveals marked differences in coastal waters PP in the Mediterranean Sea. Annual values range from 215±124 g C m -2 in the north Alboran Sea (Z1) to 48±17 g C m -2 along the coasts of Egypt and Libya (Z16). These values are globally lower than published data; yet, literature values in coastal waters are highly variable depending on methodology, depth and/or sampling date. For example, García-Gorriz and Carr (2001) estimated annual PP of 300-900 g C m -2 for Z1 but 490 Morán and Estrada (2001) narrowed this range to mean values between 121 and 366 g C m -2 depending on distance from the coast. Pugnetti et al. (2008) reported mean values of 150 g C m -2 that are almost twice higher than our values at Z8. In the lower range, Sournia (1973) estimated 30-60 g C m -2 in Z16, which is in accordance with our values for this zone. In this sense, despite the limitations inherent to satellite data, the present work provides estimations based on a long data record (14 years) and a homogeneous methodology. 495 The nine characteristic temporal patterns obtained from the SOM analysis (Fig. 7) reveal small differences in PP among the different regions. Most variations are due to changes in the magnitude of annual carbon fixation, although seasonality little varies. Exceptions are R7, R8 and R9 which represent the dynamics in coastal regions regulated by terrestrial inputs. Likewise, interannual variations are highly coherent among regions, following the basin-scale pattern shown in Fig. 4, including the remarkable decline in productivity during 2012. Exceptions are R1, representing the dynamics in the Gulf of 500 Sirte and R7 in Gulf of Gabes where a different interannual variability suggests alternative sources of PP variability in this region. Indeed, the Gulf of Gabes is a peculiar region displaying consistently high Chl and PP in most studies (e.g. Bosc et al., 2004;Barale et al., 2008). Drira et al. (2008) reported high biomass and toxic dinoflagellate blooms in the inner shelf of the Gulf of Gabes where surface nitrate concentration often exceeded 1µM. This enrichment is associated with degradation of the water quality attributed to industrial and urban activities (Hamza-Chaffai et al., 1997;Zairi and Rouis, 1999). However, even though these waters may suffer from eutrophication, satellite-borne data overestimates Chl within these waters, as revealed by Katlane et al. (2011) who observed constant high turbidity and suspended matter of industrial origin affecting these waters but also, reflection from the bottom affecting MODIS data. This suggests that general Chl algorithms may be particularly inaccurate in this region.
The magnitude of coastal PP has been often related to both shelf width and magnitude of river discharge (Liu et al., 510 2010). Our data does not display a general relationship between shelf width (Q) and PPannual (Table 4 and Fig. 9). Indeed, wide shelves with important river discharge flux from the Po, the Rhone and the Nile rivers display high productivity (Z4, Z9 and Z15 > 170 g C m -2 ) whereas production is low in narrowest shelves like Z2, Z5 and Z16. However, PP in some regions with important river inflows, like Z10, are significantly lower (89±37 g C m -2 ). In other regions like Z1 and Z8 PP is high despite the lack of important freshwater sources. 515 The role of river discharges depends both on Q and on nutrient loads. Most Mediterranean rivers have lost their natural flows and their discharges to the sea are strongly regulated by dams and water abstractions. Consequently, their outflow to the Mediterranean Sea is highly uncoupled from weather and climate variability. For instance, some rivers flowing into the Adriatic and Ionian Seas, like the Acheloos, Nestos or Aliakmon nowadays present high to maximum discharge in July due to peak hydropower production (Skoulikidis et al., 2009). 520 Figure 9: Relation between primary production, shelf width and river discharge flow (Q). Bubbles sizes are proportional to the annual mean PP (PPannual) of each of the 18 defined zones (see Fig. 8 and Table 4).
Nutrient and organic matter loads have globally increased during the last century (Beusen et al., 2016). Concentrations 525 exported to nearby seas depend on the combined effects of lithology, urban effluents, industry and agriculture in catchment basins that are often difficult to quantify. The Land Use and Land Cover (LULC) data collection provide indices of the threat https://doi.org/10.5194/bg-2020-457 Preprint. Discussion started: 11 January 2021 c Author(s) 2021. CC BY 4.0 License. of potential development for setting land and water quality policies. Rivers like the Rhone and Po with important influence on coastal productivity flow through extensive areas and therefore accumulate the impact from anthropogenic activities.
Agricultural practices and urban effluents can strongly determine the concentration and molar ratios of the nutrients flowing 530 into coastal waters. For example, despite the flux of the Nile river has been drastically reduced after the operation of the Aswan Dam (from 47 to 17 km 3 y -1 , Ludwig et al. (2009)), the coastal region is still highly productive. A remarkable increase in the concentrations of nitrate derived from fertilizers and sewage is responsible for this sustained productivity (Turley, 1999;Nixon, 2003Nixon, , 2004. Conversely, pollution pressures in the western Balkan basins are relatively low and the Neretva (Z10), running through a karstic region in Croatia, displays low nutrient levels (Ludwig et al., 2009;Skoulikidis et al., 2009). 535 Finally, other oceanographic processes determine the productivity of coastal regions. In particular, Z1 and Z18, in the Alboran Sea are comparatively more productive than other areas. The influence of winds and circulation patterns favouring subsurface water upwelling higher productivity in the northern Alboran Sea where described by García-Gorriz and Carr (2001). Also, localized patterns of relatively high primary production were found in persistent deep water density fronts resulting from the interaction of MAW and Mediterranean water by Lohrenz et al. (1988). 540

Conclusions
In summary, pelagic PP in coastal shelves of the Mediterranean Sea during the period 2002-2016 was estimated in this study for the first time using available satellite ocean colour product. We estimated that 12 % of PP of the Mediterranean Sea is attributable to coastal pelagic production and from that, about 80% of this carbon fixation is sustained by regenerated pathways.
High PP spatial variations were observed among the different regions, as mainly driven by major river effluents. Our analysis 545 also reveals a weak negative PP trend, which cannot be qualified as climate-driven because most of the temporal variability is dominated by interannual or sub-decadal variations and the satellite record is only 14-year long. Finally, we identify 18 alongshelf zones based on their temporal PP patterns. Two main PP groups were observed: zones with strong cross-shore gradients, typically found in wider estuarine regions and homogeneous zones within narrow continental shelf areas.

Competing interests 550
The authors declare that they have no conflict of interest.

Author contribution
AR and GB designed the study in collaboration with DA. PMS and DA conducted most of the analyses. PMS wrote the manuscript, with substantial contributions from all co-authors. https://doi.org/10.5194/bg-2020-457 Preprint. Discussion started: 11 January 2021 c Author(s) 2021. CC BY 4.0 License.