Biogeochemical control of marine productivity in the Mediterranean Sea during the last 50 years.

The temporal dynamics of biogeochemical variables derived from a coupled 3-D model of the Mediterranean Sea are evaluated for the last 50 years (1960-2010) against independent data on fisheries catch per unit effort (CPUE) for the same time period. Concordant patterns are found in the time series of all of the biological variables (from the model and from fisheries statistics), with low values at the beginning of the series, a later increase, with maximum levels reached at the end of the 1990s, and a posterior stabilization. Spectral analysis of the annual biological time series reveals coincident low-frequency signals in all of them. The first, more energetic signal peaks around the year 2000, while the second, less energetic signal peaks near 1982. Almost identical low-frequency signals are found in the nutrient loads of the rivers and in the integrated nutrient levels in the surface marine ecosystem. Nitrate concentration shows a maximum level in 1998, with a later stabilization to present-day values, coincident with the first low-frequency signal found in the biological series. Phosphate shows maximum concentrations around 1982 and a posterior sharp decline, in concordance with the second low-frequency signal observed in the biological series. That result seems to indicate that the control of marine productivity (plankton to fish) in the Mediterranean is principally mediated through bottom-up processes that could be traced back to the characteristics of riverine discharges. The high sensitivity of CPUE time series to environmental conditions might be another indicator of the overexploitation of this marine ecosystem. KEY POINTS: Biogeochemical evolution of the Mediterranean over the past 50 yearsRiver nutrient loads drive primary and secondary productionsStrong link between low trophic levels and fisheries.

1. Introduction

The Mediterranean Sea is the largest semienclosed water body on Earth [Azov, 1991] and is surrounded by heavily populated areas. Anthropogenic pressure on this system is high, from fisheries to tourism, land occupation, and inputs of substances and energy into the ecosystems. An adequate understanding of the consequences of human actions is thus a prerequisite for management plans aiming at securing sustainability while maximizing ecosystem services.

Given its limited connection with the open ocean and its horizontal and vertical dimensions, the Mediterranean basin has been termed as a natural laboratory, where conducting oceanographic research is possible under quite controlled conditions. Its strategic position as the main southern European sea has facilitated a long tradition of oceanographic research in this basin. Hence, there is a considerable amount of information on its hydrodynamic and biogeochemical characteristics (for a comprehensive review, see Siokou-Frangou et al. [2010]) but with scattered spatiotemporal resolution.

The main drivers of biological productivity in the basin still need to be identified. Recent works have questioned the limiting role of different nutrients for primary production [Huertas et al., 2012] and disputed the main sources of nutrients to the basin [Siokou-Frangou et al., 2010; Ribera d'Alcalá et al., 2003; Ludwig et al., 2009] and what is driving the temporal dynamics of nutrient loads into the marine environment [Ludwig et al., 2010]. The eastern Levantine basin is mainly limited by phosphate, while nitrate is the limiting nutrient in the western basin [Krom et al., 2004]. Water interchange through the Strait of Gibraltar has been identified as a key controlling mechanism for the nutrient budget of the basin [Macías et al., 2007; Huertas et al., 2012] while main sources of allochthonous nutrients are atmospheric deposition [Ribera d'Alcalá et al., 2003] and riverine runoff [Ludwig et al., 2009].

The limited period covered by the instrumental record on the biogeochemical characteristics of the basin makes it difficult to perform a comprehensive study on the links between nutrient loads and biological productivity for an extended time period. An alternative method that complements the observational record is the use of numerical coupled models that allow for hindcast simulations of the ecosystem's characteristics, considering the main external forcings. Such models are very useful tools for studying the dynamics of complex systems, as they provide mechanistic explanations and quantitative estimations for observed patterns. These models can also be used to study the ecosystem with high spatial and temporal resolutions that are impossible to obtain by other means, and they could provide a glimpse into unobserved periods through hindcast and forecast simulations.

Model validation is a fundamental step in testing the validity of the assumptions included in the modeling system [Jorgensen and Fath, 2011]. Unfortunately, the availability of independent data on the main characteristics of an ecosystem is constrained by the instrumental record. For a basin-wide hydrodynamic-biogeochemical coupled model, the observational period with sufficient data to allow for a proper model validation is mainly limited to the satellite era, spanning from 1985 for sea surface temperature [e.g., Macías et al., 2013] and from 1996 for ocean color [e.g., Macías et al., 2014]. Field data are much scarcer and scattered and usually cannot be used for long-term, basin-wide model validation.

The usual procedure to validate a hindcast simulation extending several decades back in time is to assess model performance against available data during the observational period and assume that model-data relationships hold during the whole simulated period. Additionally, a secondary proxy for biological production covering the period of interest could be pursued. One option is to use fisheries data on catches and efforts because both are usually available in National Statistics for the last few decades.

However, in that case, caution is necessary because fisheries yields are usually determined by a combination of bottom-up and top-down processes, the so-called bot-top control [Ruiz et al., 2009]. Catch per unit effort (CPUE) data are often used in the assessment of fish populations as an index of relative abundance, assuming proportionality between CPUE and abundance [Myres and Worm, 2003; Watson et al., 2013]. This proportionality is particularly important in ecosystems such as the Mediterranean Sea, where CPUE data are often the only form of data available to assess temporal evolution and patterns of distribution and abundance of exploited populations.

In the present work, we used a coupled 3-D hydrodynamic-biogeochemical model of the entire Mediterranean Sea to study the time evolution of its biochemical characteristics during the last 50 years. Model results were compared with independent fisheries statistics to assess its performance, and a detailed analysis of the possible causes of the observed pattern is presented.

2. Material and Methods

A fully coupled 3-D hydrodynamic-biogeochemical model of the entire Mediterranean basin (Figure 1) was run, covering the period from 1957 to 2012. The first 3 years (1957–1959) was not used in the analysis to avoid any influence of the initial conditions (spin-off period). The coupling of hydrodynamic and biogeochemical models was performed using the Framework for Aquatic Biogeochemical Model (FABM) [Bruggeman and Bolding, 2014]. The FABM is a two-way coupled model system in which hydrodynamics modify biogeochemistry by water movement, substance transport, light availability, and temperature-dependent process rates and in which biogeochemistry influences water column properties through light attenuation modifications by phytoplankton shelf-shading [Burchard et al., 2006].

Figure 1

Model domain. Background color indicates bottom depth as included in the model configuration (isobaths of 100, 1000, and 3000 m are also indicated). Blue stars along the coast indicate the position of the 37 rivers included in the model setup. The three main rivers of the Mediterranean basin are indicated by larger red dots.

The model integration time step was 27 s, and monthly mean values of all hydrodynamics and biogeochemical variables were stored throughout the entire simulation time span. Details of the two models used are provided in sections 2.1 and 2.2.

2.1. Hydrodynamic Model

We used the 3-D General Estuarine Transport Model (GETM) 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]. Detailed descriptions of the GETM equations can be found in Stips et al. [2004] and at http://www.getm.eu. The configuration of the Mediterranean Sea is the same as described in Macías et al. [2013] and has a horizontal resolution of 5' × 5', using 25 vertical layers. ETOPO1 [Amante and Eakins, 2009] was used to build the bathymetric grid by averaging to the corresponding horizontal resolution of the model grid. The salinity and temperature climatologies required at the start of the model integration were obtained from the Mediterranean Data Archaeology and Rescue-MEDAR/MEDATLAS database (http://www.ifremer.fr/medar/). The Strait of Gibraltar is described as an open boundary, and the Dardanelles inflow is treated as a riverine inflow within the basin. The current configuration of the model includes 37 rivers discharging along the Mediterranean coast (Figure 1). The corresponding river discharges (water flow) are derived from the Global River Data Center (Germany) database. Nutrient content (nitrate and phosphate) of freshwater runoff were obtained from Ludwig et al. [2009].

Detailed information for the water fluxes and nutrient content in each river and every time step is not available for all of the areas. Climatologic values were used to fill in the existing gaps within the data series. The best available time series of nutrient concentrations corresponded to the three major rivers, namely, the Po, the Rhone, and the Ebro, which represent 76% of the total freshwater discharges into the Mediterranean Sea [see Ludwig et al., 2009, Table 3]. It is typically assumed that these rivers also drive the overall characteristics of freshwater discharges into the basin. Therefore, only the information from these three rivers was used in the time series dynamic analysis shown below.

Table 3

Fitting Parameters for the Statistical Comparisons of Signals 1 of the Different Biological Variables and the Mean Marine Nitrate Concentration Low-Frequency Signal

Fitted Variable With Marine Nitrate	Equation	r2
Integrated PPR (signal 1)	PPR = 3.9 × 10−2 + 2.64 × 100 × Mar_nit	0.99 (p < 0.001)
Integrated zoo biomass (signal 1)	Zoo = −2.82 × 10−3 + 7.4 × 10−1 × Mar_nit	0.99 (p < 0.001)
CPUE (signal 1)	CPUE = -2.3 × 10−5 + 3.39 × 10−3 × Mar_nit	0.99 (p < 0.001)

The GETM run for the Mediterranean Sea was forced at the surface every 6 h with European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis products. Specifically, we used the ECMWF ERA40 reanalysis products from 1957 to 1978 and the ERA-Interim products from 1979 to 2012. The consistency of these data sets was checked to avoid spurious results. The ECMWF ERA reanalysis data are available from the ECMWF data server (http://www.ecmwf.int).

2.2. Biogeochemical Model

A modified version of the Ecological Regional Ocean Model (ERGOM) [Neumann, 2000] was used to simulate the low trophic levels of the pelagic ecosystem of the entire Mediterranean basin. ERGOM is a nitrogen-based model used to describe the transfer of nitrate between the different trophic compartments. Three phytoplankton types (diatoms, flagellates, and cyanobacteria) incorporate dissolved nitrogen (as nitrate and ammonium) by photosynthesis, while a single zooplankton state variable predates on the different phytoplankton groups. Dead particles and unassimilated grazing go into a detritus pool that is remineralized at a temperature-dependent rate. Both large phytoplankton (diatoms) and detritus, which escape water column remineralization and reach the sea floor, are subject to benthic processes. Those processes include sedimentary burial, which effectively represents nutrient removal from the model, remineralization within the sediment, and resuspension. Phosphate dynamics are linked to nitrate through the constant Redfield ratio [Redfield, 1934]. Nitrate, ammonium, and phosphate all limit primary production, depending on the relative affinity of each phytoplankton group to each particular nutrient.

Modifications to the original ERGOM code were made to better represent the characteristics of the Mediterranean ecosystem (MedERGOM) and are extensively described in Macías et al. [2014]; thus, only a brief description is included herein. Zooplankton mortality losses were modified by including two linear expressions (excretion and mortality) and a quadratic term (predation) following the recommendations for the Mediterranean Sea provided by Oguz et al. [2013]. The predation term is treated as the closure of the model and, thus, is lost from the system.

The light limitation of primary production was also changed in this implementation of MedERGOM. Instead of using a single set of light limitation values for all three phytoplankton types, specific values of the production-irradiance curves were adopted for each functional type [Follows et al., 2007]. The functional response of growth rate to light levels was also changed following Wan and Bi [2013].

2.3. CPUE

The time series of CPUE was used in the analysis to compare temporal evolutions of Mediterranean exploited populations with selected biogeochemical variables from the model simulation (e.g., primary production rate and zooplankton biomass). In particular, total annual fish catches in the Mediterranean basin from 1960 to 2012 were obtained from the Food and Agriculture Organization (FAO) database, accessible at http://www.fao.org/fishery/statistics/global-capture-production/en (data were downloaded in January 2014). Effort was estimated by taking into account the number of boats, length, and gross tonnage information obtained from the General Fisheries Commission for the Mediterranean (GFCM) database on the Authorized Vessel List and from the FAO Fishing Fleet online database [Food and Agriculture Organization (FAO), 2010a].

2.4. Singular Spectrum Analysis

Singular spectrum analysis (SSA) is designed to extract information from short and noisy time series and provides insight into the unknown, or only partially known, dynamics of the underlying system that generated the series [Ghil et al., 2002]. This methodology is analogous to applying an extended empirical orthogonal function analysis to successive lags of a univariate time series and is equivalent to representing the behavior of the system by a succession of overlapping “views” of the series through a sliding n point window [Vautard et al., 1992]. The SSA puts the decomposition of the time series into a sequence of elementary patterns of behavior that are classified as either trends or oscillatory patterns. From this decomposition into eigenvalues, it is possible to reconstruct each of the individual signals by adding the corresponding eigenvectors to the sample mean [Vautard et al., 1992].

SSA has been applied to biogeochemical variables time series from the coupled model system. Integrated values for the top 120 m of the water column were computed for nitrate, phosphate, primary production rate (PPR), and zooplankton biomass. The same analysis was performed on the fisheries data set (total fish catches and CPUE) and on the rivers' nutrient loads with the identical number of extracted eigenvectors (25) and using the same methodology [Broomhead and King, 1986] to calculate the covariance matrices (i.e., assuring a homogeneous analysis of all time series).

3. Results

Annual mean values of PPR and zooplankton biomass integrated in the upper 120 m of the water column were computed for the entire 1960–2010 time period (blue and orange lines in Figure 2a). The fisheries CPUE data series covering the same time period is plotted against time in Figure 2a. The temporal evolutions of all variables follow a similar pattern, with low values during the first 15 years (approximately 1960–1975), a substantial increase until approximately 1985 and a later stabilization until 2010. There are statistically positive and significant relationships between biogeochemical variables and the CPUE data (Figure 2b). Correlation coefficients are higher with zooplankton biomass than with PPR (Figure 2b), although in both cases R is above 0.66. For all of the correlation analyses, the probability of significance (p) was adjusted to correct for temporal autocorrelation following the recommendations by Pyper and Peterman [1998]. In all cases, correlations are significant to a 95% confidence level.

Figure 2

(a) Annual time series of biological variables; integrated primary production rate (PPR), blue line (mmol N m−2 d−1); integrated zooplankton biomass, orange line (mmol N m−2 d−1) and catches per unit effort (CPUE), black line (tons GT−1) during the analyzed period. (b) Scatterplot of PPR versus CPUE (blue dots) and zooplankton biomass versus CPUE (orange dots). Statistics of the linear fittings (blue and orange lines) are inserted.

In spite of the statistical significance of the relationships shown in Figure 2b, the contours of the temporal time series of the different variables (Figure 2a) do not have a distinct shape. These contours mainly follow the transition from a lower biological productivity at the beginning of the analyzed period to a state of larger production levels at the end. Therefore, considering that correlation does not entail causation [Sies, 1988], SSA was applied to all four time series for insight into the underlying processes shaping their temporal dynamics [Ghil et al., 2002].

Two main signals were identified as the main contributors to the total variability of the time series and explained between 69% and 90% of their total variability (Table 1). The identified signals had specific and concordant shapes (Figures 3a and 3b). Signal 1, which explains approximately 50% of the total variability of the individual series (Table 1), has a sigmoidal shape (Figure 3a) with minimum values near 1965 and maximum values near 1998. Signal 2 accounts for 25–30% of the total variability (Table 1) and has a distinctive peak during the second half of the 1980s (Figure 3b), with lower values before and after. The scatterplots of the signals in the biogeochemical variables time series against the signals in the CPUE data (Figures 3c and 3d) show very high (r2 > 0.95), significant (p < 0.001) relationships, indicating that the isolated signals could be representing common processes driving the main dynamics of the different time series.

Table 1

Energy of the Annual Series Associated to Each of the Identified Signals With the SSA

Variable	Signal 1	Signal 2	Σ Energy
Integrated PPR	43.6%	23.21%	66.8%
Integrated zoo biomass	52.27%	26.77%	79.05%
CPUE	51.6%	27.43%	83.1%
Marine nitrate	59.6%	--	59.6%
Marine phosphate	69.7%	--	69.7%
Riverine nitrate	59.13%	--	59.13%
Riverine phosphate	41.8%	--	41.8%

Figure 3

(a) Main signals found in the biological time series using SSA. (b) Secondary signals found in the biological time series using SSA. (c) Scatterplot of the main signals found in PPR versus the main signals found in the CPUE time series. Statistics of the linear fit is inserted. (d) Scatterplot of the main signals found in zooplankton biomass versus the main signals found in the CPUE time series. Statistics of the linear fit is inserted.

The nutrient loads in the three main rivers of the basin (Ebro, Rhone, and Po) [Ludwig et al., 2009] were analyzed with the same SSA methodology to test if bottom-up processes could be held responsible for the observed signals in the biological time series. In this case, total nitrate and phosphate loads show a unique signal as the main contributor to the total energy (Table 1), which accounts for 60–70% of the energy of the nutrient series. The nitrate signal peaks around the year 2000, while phosphate peaks near 1985 (Figure 4a).

Figure 4

(a) Main signals in the nutrients time series from the Ebro, Rhone, and Po rivers obtained with the SSA. (b) Main signals in the integrated (0–100 m) nutrients time series from the entire Mediterranean basin obtained with the SSA. (c) Scatterplot of the nitrate signal in the entire basin versus the first signal obtained from the different biological time series. (d) Scatterplot of the phosphate signal in the entire basin versus the second signal obtained from the different biological time series.

Mean integrated (0–120 m) nutrient concentrations in the whole Mediterranean basin computed from the model simulation showed different main signals for nitrate and phosphate (Figure 4b), which are strongly correlated with the signals found in the rivers' nutrient loads (Table 2). Again, the nitrate concentration peaks around the year 2000, while phosphate peaks during early the 1980s.

Table 2

Statistical Fit Between the Low-Frequency Signals in the Riverine Nutrient Loads and the Low-Frequency Signals in the Mean Marine Nutrient Concentrations

Equation	r2
Mean_marine_nit = 4.38 × 10−2 + 5.18 × 10−3 × Riv_nit	r2 = 0.93 (p < 0.01)
Mean_marine_pho = 3.7 × 10−7 + 3.24 × 10−5 × Riv_pho	r2 = 0.71 (p < 0.05)

Figure 4c shows that the first low-frequency signals found in the different biological variables (model and fisheries data) are strongly and significantly correlated (Table 3) with the signal found in the mean nitrate concentration of the surface layer in the whole Mediterranean basin. The second signal found in the biological records is more correlated (Table 4) to the low-frequency signal present in the phosphate concentration of the basin (Figure 4d). In all cases (Tables 3 and 4), correlation coefficients between nutrients and biological variables are above 0.94 with significance levels < 0.001.

Table 4

Fitting Parameters for the Statistical Comparisons of Signals 2 of the Different Biological Variables and the Mean Marine Phosphate

Fitted Variable With Marine Phosphate	Equation	r2
Integrated PPR (signal 2)	PPR = −1.66 × 10−1 + 5.4 × 104 × Mar_pho	0.98 (p < 0.001)
Integrated zoo biomass (signal 2)	Zoo = −3.87 × 10−2 + 1.36 × 104 × Mar_pho	0.99 (p < 0.001)
CPUE (signal 2)	CPUE = −1.24 × 10−4 + 5.37 × 101 × Mar_pho	0.98 (p < 0.001)

4. Discussion and Conclusions

The same model configuration employed here was compared and validated with satellite data on SST [Macías et al., 2013] and surface chlorophyll [Macías et al., 2014] during the observational period. In both cases, model-data agreement is acceptable in the majority of the Mediterranean basin, indicating that the model could be regarded as a reasonable descriptor of the present state of the system. From the results presented here, there are also obvious similarities between the annual time series of modeled biogeochemical variables and the fisheries records during the last 50 years in the Mediterranean basin (Figure 2a). The scatterplot of biogeochemical variables derived from the model run versus the independent fisheries statistics (Figure 2b) indicates that there is a quite strong and significant correlation between the variables, even if the strong temporal autocorrelation of the time series is taken into account [Pyper and Peterman, 1998].

The first conclusion derived from this comparison is that the proposed model setup captures the main forcings of marine productivity of the basin during the second half of the past century. The main sources of nutrients in the modeled domain are the riverine inputs, so it is reasonable to assume that river dynamics were responsible for a significant fraction of the temporal variability in marine productivity as previously proposed [Ludwig et al., 2010]. The surface distribution of plankton biomass simulated by the model [see Macías et al., 2014] also supports this hypothesis, as maximum values were usually simulated in the vicinity of the rivers' mouths, as also observed in the satellite-derived chlorophyll maps [e.g., D'Ortenzio and Ribera D'Alcala, 2009].

The analysis of the nutrient loads for the three main Mediterranean rivers for which the best information is available [Ludwig et al., 2009] reveals an already described pattern (Figure 4a). The nitrate load increased continuously during the last decades and reached maximum values at the beginning of the 21st century [Ludwig et al., 2009], with a slight, insignificant decline afterward. Phosphate, in contrast, reached maximum values during the first half of the 1980s, with a sharp decline later [Ludwig et al., 2010]. This pattern of phosphate levels is linked to the introduction of a nonphosphate regulation in detergents in the 1980s [Ludwig et al., 2009].

Curiously, it seems that the main driver of the biological production time evolution in the entire basin was the increase of nitrate in the rivers and, consequently, in the marine ecosystem (Figure 4b) up to the year 2000. The nitrate concentration is strong and significantly correlated with the first low-frequency signal of all of the biological variables (representing over 50% of total variability), with correlation coefficients over 0.97 in all cases (Table 3), prompting questions about the relative limiting role of the different macronutrients in this basin.

Usually, the observed N:P ratio in the Mediterranean Sea is approximately 28:1 in deep waters [Krom et al., 2004] and even higher (up to > 100:1) in surface waters [Huertas et al., 2012; Christodoulaki et al., 2013]. It is assumed, thus, that phosphate is the limiting nutrient, especially in the eastern basin [Siokou-Frangou et al., 2010]. Mean climatologic values reported in the World Ocean Atlas were NO3 = 0.04 µM and PO4 = 0.0011 µM for the first 125 m of the water column, yielding a value of N:P ∼36:1. Mean surface nutrient values computed from the model during the simulated period were NO3 = 0.05 µM and PO4 = 0.0009 µM. These values are reasonably similar to the observed values and provide a mean simulated N:P ∼55:1.

In spite of the observed and simulated phosphate limitation, simulated phosphate dynamics in the basin (also clearly linked to riverine loads (Table 2 and Figures 4a and 4b)) seem to be mainly driving the second (and not the first) low-frequency signal found in the biological variables time series, which are responsible for just 25% of the total variability of the series (Table 1). This discrepancy does not seem to be a problem with the biogeochemical model used, as the same was valid for the CPUE data; thus, it must be an intrinsic property of the studied system. A similar result was recently shown for the eastern Mediterranean by Christodoulaki et al. [2013] on a 1-D modeling study. In that case atmospheric nutrient inputs were considered in the model; such inputs not only further augmented P limitation in the surface water layer but they also increased the primary productivity of the region. These findings lead to new questions on the relative limitation role of the different nutrients in the Mediterranean basin [e.g., Huertas et al., 2012], but this is an issue beyond the scope of the present contribution.

It was recently suggested [Huertas et al., 2012] that the nutrient ratios and their temporal dynamics in the Mediterranean Sea are determined mainly by riverine and atmospheric inputs. Even if the antiestuarine circulation in the Strait of Gibraltar export high-nutrient waters and import impoverished surface waters [Macías et al., 2007], transformation and cycling within the basin seem to be the main driving process of nutrient dynamics in the Mediterranean [Huertas et al., 2012]. The strong links between river loads and marine nutrient concentrations found in our results (Table 2) seem to further support this hypothesis.

The second conclusion derived from our results is that the total biological productivity of the Mediterranean basin seems to be driven mainly by bottom-up processes. Nutrient availability and relative proportion mediated through riverine inputs seem to control the biological characteristics of the region. The spectral analysis of the different annual time series clearly revealed that the main low-frequency signals within each series are concordant, showing very similar time evolutions (Figures 3a and 3b). This result could be considered as an indication that the driving processes shaping the temporal dynamics of biogeochemical variables and CPUE are the same.

The slopes of the correlations between mean integrated nutrients and PPR are larger than the slopes of the correlations with the integrated zooplankton biomass (Figures 4c and 4d and Tables 3 and 4). This result indicates a stronger dependence of primary productivity on nutrient availability and enhances our confidence in the likelihood of the simulated results.

Previous reports based on field data analysis [Ludwig et al., 2010] and on modeling simulations [Lazzari et al., 2012] have indicated that river nutrient loads could be responsible for a significant fraction of biological production in the Mediterranean basin. However, this is the first time in which a strong and significant correlation holding during a substantial time period was described. Additionally, this work provides the first quantitative evidence of the close relationship between the chemical conditions, the productivity at the lower trophic levels, and the fish abundance in the basin. Our simulations of marine primary productivity evolution in the Mediterranean fit, also, quite well with a recently reported long-term (∼150 years) trend of surface Chla in the region [Boyce et al., 2014].

The Mediterranean Sea is generally considered an oligotrophic basin [Antoine et al., 1995] with some local regions of enhanced productivity [Siokou-Frangou et al., 2010]. The general oligotrophy of the basin is interrupted by isolated primary production hotspots in some coastal areas in the Alboran Sea [García-Gorriz and Carr, 2001], in the vicinity of river mouths [Meybeck et al., 2007], and in highly dynamic regions with strong frontal activity, such as the Balearic Sea [La Violette, 1990]. In spite of the general low primary productivity of the Mediterranean basin, fisheries yields are not as low as would be expected from an oligotrophic system. The Mediterranean basin represents only approximately 0.8% of the total marine area of the world [Azov, 1991], whereas its fisheries account for nearly 2.2% of the total marine fish catches during the last 50 years (FAO statistics).

In this basin, there is relatively abundant fish production in spite of low phytoplankton biomass, which has been named as the “Mediterranean paradox” [Sournia, 1973] given the assumption that the oligotrophic character of the region should normally constrain the feeding activity and production of zooplankton [Saiz et al., 2007]. The existence of a strong and very efficient link between low and high trophic levels was previously proposed to explain this paradox [Estrada, 1996], and our results seem to support this hypothesis, with fisheries data closely following the biogeochemical patterns derived from the model simulations (Figures 2 and 3).

Moreover, the high level of exploitation of fish stocks in the Mediterranean [Papaconstantinou and Farrugio, 2000] could make populations less resilient and more sensitive to environmental factors [Hsieh et al., 2008], such as food availability. This larger sensitivity to environmental conditions could be partially responsible for the extremely good correlations found between the biogeochemical characteristics of the basin and fish abundance (Figures 3c and 3d).

Catches of different fish groups as small pelagic and demersal species show low-frequency signals almost identical to the ones found in CPUE (not shown). This result is consistent with previous reports [e.g., Fiorentini et al., 1997] and indicates that a high proportion of species and groups of species exhibit increased landings over the last 45 years. Such an increase was observed both in the western and eastern basins and was previously attributed to potential changes in environmental conditions [Fiorentini et al., 1997]. Our model-data comparison also points to environmental factors (i.e., nutrient concentration) as the most probable cause of fisheries yield increases in the Mediterranean.

Still, there is clear evidence that in the Mediterranean Sea, as in other parts of the world, several stocks have been declining as a consequence of intensive fishing pressure [Lleonart and Maynou, 2003; Tudela, 2004]. The strong increase in total CPUE observed in our data (Figure 2a) could be partially attributed to a number of different factors. For example, the improvement in the fishing efficiency of vessels (e.g., fishing gears and onboard electronics) has largely increased the accessibility to fishing resources. Additionally, the geographic expansion of fishing fleets to deeper areas and to Northern African countries, where fishing resources are still abundant, is driven by declines of inshore catches and aided by improved technology and supported by subsidies [Sumaila et al., 2010; Swartz et al., 2010; Watson et al., 2011]. The use of fisheries statistics (both for catch and effort) supplied by individual countries to FAO and GFCM, which were not spatially explicit, could be another limiting factor. Several studies have confirmed that most of these statistics largely underestimate their likely true catch [Zeller and Pauly, 2007; Pauly et al., 2014], particularly in the Southern Mediterranean where mechanisms to collect fisheries data are less available [FAO, 2010b]. Finally, the high species diversity in catch composition (more than 150 important commercial species) of Mediterranean fisheries makes it difficult to detect a decline, changing patterns and trends of individual stocks.

Even taking into account all previous considerations, freshwater management reveals to be a key element controlling anthropogenic effects on this heavily populated and exploited marine basin. This is especially important for the future evolution of this region [Ludwig et al., 2010; Lazzari et al., 2012] given the relevance that river nutrient loads appear to have in determining marine primary and secondary productions. Our results clearly indicate that the final impacts of riverine discharges need to be considered and integrated into the management plans of continental waters, as any measurement taken on waste water management, agricultural plans, and navigation routes will have an impact on the total marine productivity of the Mediterranean Sea.

To perform any forecasting simulation on expected conditions in the basin, it is crucial to have well-defined and sound river management scenarios. Temporal and spatial fluctuations of productivity patterns (including potential changes of fisheries ground) are very likely to be influenced by future discharges of rivers. Henceforth, forecasting exercises should contemplate potential changes in river characteristics together with expected alterations in atmospheric and/or boundary conditions. Doing so entails using the same story lines for atmospheric and riverine conditions in order to obtain consistent scenarios of future conditions influencing the Mediterranean Sea characteristics.