Ecoregions in the Mediterranean Sea Through the Reanalysis of Phytoplankton Functional Types and Carbon Fluxes

In this work we produced a long ‐ term reanalysis of the phytoplankton community structure in the Mediterranean Sea and used it to de ﬁ ne ecoregions. These were based on the spatial variability of the phytoplankton type fractions and their in ﬂ uence on selected carbon ﬂ uxes. A regional ocean color product of four phytoplankton functional types (PFTs; diatoms, dino ﬂ agellates, nanophytoplankton, and picophytoplankton) was assimilated into a coupled physical ‐ biogeochemical model of the Mediterranean Sea (Proudman Oceanographic Laboratory Coastal Ocean Modelling System ‐ European Regional Seas Ecosystem Model, POLCOMS – ERSEM) by using a 100 ‐ member ensemble Kalman ﬁ lter, in a reanalysis simulation for years 1998 – 2014. The reanalysis outperformed the reference simulation in representing the assimilated ocean color PFT fractions to total chlorophyll, although the skill for the ocean color PFT concentrations was not improved signi ﬁ cantly. The reanalysis did not impact noticeably the reference simulation of not assimilated in situ observations, with the exception of a slight bias reduction for the situ PFT concentrations, and a deterioration of the phosphate simulation. We found that the Mediterranean Sea can be subdivided in three PFT ‐ based ecoregions, derived from the spatial variability of the PFT fraction dominance or relevance. Picophytoplankton dominates the largest part of open ocean waters; microphytoplankton dominates in a few, highly productive coastal spots near large ‐ river mouths; nanophytoplankton is relevant in intermediate ‐ productive coastal and Atlantic ‐ in ﬂ uenced waters. The trophic and carbon sedimentation ef ﬁ ciencies are highest in the microphytoplankton ecoregion and lowest in the picophytoplankton and nanophytoplankton ecoregions. The reanalysis and regionalization offer new perspectives on the variability of the structure and functioning of the phytoplankton community and related biogeochemical ﬂ uxes, with foreseeable applications in Blue Growth of the Mediterranean Sea.


Introduction
The composition of the phytoplankton community has a profound influence on the functioning of marine ecosystems, impacting ecosystem health (e.g., eutrophication and harmful algae blooms, e.g., Hallegraeff, 2003), services (e.g., aquaculture and fishery, e.g., Cury et al., 2008), and climate feedbacks through carbon sequestration (e.g., Bopp et al., 2005). Phytoplankton impact marine biogeochemical fluxes and trophic networks, in relation to their size (e.g., picophytoplankton (PIC), nanophytoplankton (NNP), and microphytoplankton (MIC) have different affinities for nutrients and only the latter are palatable to fishes) and to their functional characteristics (e.g., diatoms, typically large in size, can assimilate silicate, while nanosized coccolitophores impact the cycling of calcium; Chisholm, 1992). Phytoplankton functional types (PFTs) are defined as groups of phytoplankton that have a distinct function in the ecosystem (Baretta et al., 1995).
PFTs can be investigated by means of ocean color observations (IOCCG, 2014). A variety of algorithms have been proposed to derive PFTs biomass and production from satellite observations (Kostadinov et al., 2017). However, these data sets are limited to the first optical depth of the ocean, and the derivation of vertical structure and of fluxes other than primary production is highly uncertain (Lee et al., 2015). State-of-the-art ecosystem models can describe PFTs and the related biogeochemical and trophic fluxes, but the skill of model simulations is often poor, mainly due to the difficulty of formulating and parameterizing PFT processes (Shimoda & Arhonditsis, 2016). Merging PFT ocean color and ecosystem models through data assimilation was proposed as a novel approach to estimate better the phytoplankton community structure and related fluxes in marine ecosystems (Ciavatta et al., 2018). Better simulations of diatoms, dinoflagellates, NNP, and PIC biomass, as well as estimates of related biogeochemical fluxes, were delivered in both a reanalysis simulation (Ciavatta et al., 2018) and (pre)operational predictions (Skákala et al., 2018) in the North West European Shelf Sea ecosystem, through the assimilation of a regional ocean color PFT product .
In the Mediterranean Sea, regional ocean color algorithms and products were used to characterize the spatialtemporal variability of PFTs (Di Cicco et al., 2017;Navarro et al., 2014Navarro et al., , 2017Sammartino et al., 2015;Uitz et al., 2012;Volpe et al., 2007). Observation-based analysis showed that small-sized groups, either PIC (Magazzù & Decembrini, 1995) or NNP (Uitz et al., 2012), can dominate the primary production in the Mediterranean Sea, although MIC production can also be relevant in some spring-blooming Mediterranean areas (Uitz et al., 2012). Modeling studies have investigated the variability of primary production (e.g., Di Biagio et al., 2019;Lazzari et al., 2012Lazzari et al., , 2016, but they only focused on the phytoplankton community as a whole. Model and ocean color total chlorophyll were used to define ecoregions in the Mediterranean Sea (Ayata et al., 2018;Basterretxea et al., 2018;d'Ortenzio & Ribera d'Alcalà, 2009;Mayot et al., 2016). Ecoregions can be defined as areas characterized by distinctive and rather homogeneous ecosystem features, such as the phenology of the total phytoplankton community, as conducted in the above cited papers. The above total chlorophyll-based regionalizations, along with regionalizations based on hydrodynamics and biogeochemistry, have been recently integrated as layers of a "consensus" map, aiming to support the understanding of the ecosystem functioning as well as the implementation of the Marine Strategy Framework Directive in the Mediterranean Sea (Ayata et al., 2018). To the authors' knowledge, a regionalization based on PFTs from either ocean color or models has not been carried out yet, despite PFTs having profound effects on the ecosystem features, such as biogeochemical fluxes and trophic webs.
The objective of this work is to define PFT-based ecoregions in the Mediterranean Sea and to characterize their carbon fluxes, by merging information on PFTs from model simulations and ocean color observations. To achieve this objective, we assimilated a regional ocean color product for PFTs (Di Cicco et al., 2017) into an assimilative modeling system for the Mediterranean Sea. This system uses the coupled physical-biogeochemical model Proudman Oceanographic Laboratory Coastal Ocean Modelling System (POLCOMS)-European Regional Seas Ecosystem Model (ERSEM) of the Mediterranean Sea developed by Holt et al. (2008) and updated by Kay and Butenschön (2018), and we applied the PFT assimilation framework developed by Ciavatta et al. (2018). Previous applications of ERSEM to investigate the Mediterranean Sea biogeochemistry, in both assimilative and not assimilative simulations, include the works by Vichi et al. (1998), Allen et al. (1998Allen et al. ( , 2002Allen et al. ( , 2003, Petihakis et al. (1999Petihakis et al. ( , 2002Petihakis et al. ( , 2007, Triantafyllou et al. (2000Triantafyllou et al. ( , 2003, Zavatarelli et al. (2000), Blackford (2002), Pinardi et al. (2003), Kalaroni et al. (2016), and Tsiaras et al. (2017). A skill assessment was performed to evaluate the reanalysis adequacy in representing the average spatial variability of the phytoplankton community structure. The ecoregions have been defined by computing multiannual average of reanalyzed surface PFT fractions, assessing spatial patterns of the PFT dominance or relevance, and evaluating the spatial variability of vertically integrated carbon fluxes. Holt & James, 2001) describes the hydrodynamics. This provides the physical forcing to the pelagic biogeochemical submodel, namely, the ERSEM (Baretta et al., 1995;Butenschön et al., 2016). The third submodel is the ERSEM benthic biogeochemical model (Blackford, 1997;Butenschön et al., 2016). The two ERSEM submodels are coupled at the same temporal and spatial resolution as the physical model to capture the effects of the three-dimensional hydrodynamics on the biogeochemical cycles (Holt et al., 2004). The model grid has a horizontal resolution of 0.1°(approximately~11 km in the meridional direction and between 9.3 and 7.8 km in the zonal direction at the latitude of the study region) and 42 sigma coordinate levels in the vertical.

The Physical Submodel: POLCOMS
The physical model POLCOMS (Holt & James, 2001) is a three-dimensional, finite difference, primitive equation model formulated in spherical-polar coordinates on an Arakawa B-grid. The model includes an advection scheme with stability and conservation properties (James, 1996), a vertical turbulence model (General Ocean Turbulence Model; Burchard et al., 1999), and calculation of horizontal pressure gradients. POLCOMS was configured for the Mediterranean Sea by Holt et al. (2008).

The Pelagic Biogeochemical Submodel: ERSEM
The biogeochemical processes are described by the ERSEM (Baretta et al., 1995), version 15.06 , as in the PFT assimilative model by Ciavatta et al. (2018). ERSEM uses a functional type approach to model the dynamics of the low trophic levels of the ecosystem. Primary producers are split into four PFTs, including two size-based types (PIC and NNP), dinoflagellates, and large diatoms (which, summed up, represent the model MIC). Each of these PFTs is represented in terms of its content of chlorophyll, carbon, nitrogen, phosphate, and (for diatoms only) silicate. Three functional types of zooplankton (mesozooplankton, microzooplankton, and heterotrophic nanoflagellates) prey on the phytoplankton, bacteria, and particulate organic matter as a function of their size. A bacterial functional type drives the microbial loop, the production and recycling of dissolved organic matter in labile and recalcitrant forms, and the regeneration of inorganic nutrients in the water column (Hansell, 2013;. In the functional types, the stoichiometric ratios of nutrients-to-carbon and chlorophyll-to-carbon (in the PFTs) vary dynamically (Baretta-Bekker et al., 1997;Geider et al., 1997). The model simulates the dynamics of five inorganic dissolved nutrients (carbon, nitrate, ammonium, phosphate, and silicate) and dissolved oxygen. Nitrogen undergoes nitrification in the water column, as a function of available ammonium, temperature, oxygen, and pH (Blackford & Gilbert, 2007;Butenschön et al., 2016). Denitrification is modeled in the sediments only .
ERSEM coupled to POLCOMS was used in the Mediterranean Sea by Kay and Butenschön (2018) and , by using the ERSEM parameterization of Butenschön et al. (2016). However, here the maximal PFT carbon-to-chlorophyll ratios were set equal to the values in the comparable Mediterranean model by Lazzari et al. (2012Lazzari et al. ( , 2016 and Teruzzi et al. (2014Teruzzi et al. ( , 2018; i.e., 0.02 mg Chl/mg C). In fact, the original ERSEM values were calibrated for the global ocean and led to ocean color overestimation in the Mediterranean Sea in our preliminary sensitivity simulations (not shown).

The ERSEM Benthic Submodel
The benthic submodel is the ERSEM benthic model (Blackford, 1997), as described in Butenschön et al. (2016). In the configuration applied here, the submodel includes 35 biogeochemical variables, subdivided into five living functional types (including three zoobenthos types, and aerobic and anaerobic bacteria), along with particulate matter and dissolved organic and inorganic nutrients. The fluxes at the sediment-water interface are determined by sedimentation of organic matter and diffusion of inorganic nutrient across the seabed. 2.1.4. Atmospheric Forcing, Boundary Constrains, and Initial Conditions Atmospheric forcing was provided by the ERA-Interim data set (Dee et al., 2011), using 6-hourly values of surface winds, sea level pressure, temperature and relative humidity, and daily accumulation of precipitation and shortwave and longwave radiation fluxes. The resolution of the surface data was 0.75°. Ocean boundary conditions at the western boundary were taken from a run of the same model configured for a wider domain  using the same surface forcing: This provided daily values of temperature, salinity, currents, nitrate, phosphate, silicate, dissolved inorganic carbon (DIC), and bioalkalinity. Climatological values for river inputs of fresh water and nutrients were taken from GlobalNEWS2 (Mayorga et al., 2010). River nutrient concentrations were set constant; riverine freshwater discharges varied seasonally using historic data from Dai et al. (2009) to set the variation. Climatological inputs for the North Adriatic (NAd) rivers, including the Po River, were derived from data in Cozzi and Giani (2011). Atmospheric depositions of bioavailable nitrogen and phosphorus were represented as in Lazzari et al. (2016), that is, temporally constant and spatially different in the eastern and western basins, with values derived from Ribera d' Alcalà et al. (2003).
The initial conditions for nitrogen were derived to match the average subbasin biogeochemical-ARGO profiles in de Fommervault et al. (2015). Phosphate initial conditions were set equal to 1/25 of the nitrogen values, that is, setting the initial N:P ratio equal to an average Mediterranean value of 25 (comparable, e.g., to Pujo-Pay et al., 2011). The initial conditions for the other biogeochemical and physical variables were derived from a model simulation on a wider model domain . A spin-up of 20 years was performed before starting the reference and assimilative simulations, having verified that this was sufficient to stabilize the patterns of the nutrients profiles.

The Assimilation System
The assimilation system is described in full in Ciavatta et al. (2016Ciavatta et al. ( , 2018. It uses the ensemble Kalman filter (EnKF; Evensen, 1994). This is a sequential assimilation method, which starts by randomly sampling an ensemble of N state vectors x h i links the four state variables representing the ERSEM phytoplankton types to the respective ocean color products (Ciavatta et al., 2018). Observations and model states are log-transformed prior to the analysis, to guarantee positivity of the solutions (Janjić et al., 2014), as in the applications by Nerger and Gregg (2007), Ciavatta et al. (2016Ciavatta et al. ( , 2018, Ford and Barciela (2017), and Skákala et al. (2018).
The PFT assimilation scheme was set up as in Ciavatta et al. (2018). We used the EnKF with an ensemble size of N = 100 members. We included in the analyzed state vector 35 out of the 51 ERSEM pelagic state variables, to keep the analysis computationally affordable, in analogy with Ciavatta et al. (2018).
Coherently with the paper's objective to reanalyze the plankton community structure, we included in the state vector all the 27 variables related to the phytoplankton, zooplankton, and bacteria functional types, as well as eight variables representing DIC, labile, and refractory dissolved organic matter and small particulate matter. The remaining 16 variables were updated through the model equation during the simulation runtime ("forecast" step; e.g., Pradhan et al., 2019). The radius of the localized analysis was set to be spatially variable as a function of the bathymetry. In particular, we increased the "resolution" of the analysis from oceanic toward coastal waters, by setting a radius of 100 km for grid points where the bathymetry is deeper than 2,000 m, 50 km for bathymetry between 50 and 2,000 m, and 25 km for bathymetry shallower than 50 m (Ciavatta et al., 2016).
The model error is accounted for by random perturbations of the model forcing, namely, the surface solar irradiance, thus inducing fluctuations in the underwater light field that drives photosynthesis (Ciavatta et al., 2016). A Gaussian perturbation with standard deviation equal to 20% of the irradiance value is added during the model forecast step. Furthermore, at the first assimilation step of each year, the model error is added to all the variables undergoing the analysis, as white noise drawn from a distribution of pseudo-random fields with error equal to 10% of the value of the variables. The error is lowered to 1% for those variables that have relatively high mean values (DIC, small particulate matter), to avoid divergences in the concentrations of the largest pool in the model (Ciavatta et al., 2011).
assimilation. The perturbation applied Gaussian pseudo-random fields with error equal to 30% of the value of the variables.
The observational error was derived from values of the bias and root-mean-square deviation (RMSD) of the PFT ocean color product (see section 2.3). The variance of the unbiased PFT concentrations regridded onto the surface model domain and was computed following the procedure described in the Appendix of Ciavatta et al. (2016).The PFT assimilation approach was developed in the framework of the EU Copernicus Marine Environment Monitoring Service (CMEMS) -Service Evolution Project TOSCA, and was applied with the operational CMEMS model of the North West European Shelf-Sea by Skákala et al. (2018).
The reanalysis simulation was run on the U.K. Natural Environment Research Council (NERC) High Performance Computing facility "ARCHER," using 7200 CPUs and~10 mega allocation units and producing 10 Tb of output.

Data
The assimilated ocean color PFT chlorophyll products were developed by Di Cicco et al. (2017) using regional PFT algorithms for the Mediterranean Sea. These algorithms are empirical relationships developed following a global approach (Hirata et al., 2011) based on the relation between phytoplankton abundance and the trophic status of the environment (Brewin et al., 2010Chisholm, 1992;e.g., Devred et al., 2006;Uitz et al., 2006). These algorithms resulted from the existing covariability between in situ PFT groups (expressed as chlorophyll fraction) and the corresponding log 10 -transformed in situ total chlorophyll a concentration, taking account of the lognormal distribution of this pigment (Chisholm, 1992;Hirata et al., 2011).
Following previous works carried out for the global ocean (Barlow et al., 1993;Gieskes et al., 1988;Uitz et al., 2006), the contribution of each PFT to the Mediterranean assemblage in terms of chlorophyll concentration was computed by performing a multiple regression analysis between in situ total chlorophyll a and seven diagnostic pigments, markers for the main algal groups (Vidussi et al., 2001). The final in situ PFT estimation formulas and the regional satellite algorithms, with their references, are fully described in  2003) was used for the exact identification of the two different water types. The time series for this work includes daily values of chlorophyll a concentration quantified for four phytoplankton groups (PIC, NNP, diatoms, and dinoflagellates) and their associated errors (bias and RMSDs) at 2-km resolution. The ocean color PFT errors were estimated through model validation using an independent PFT validation in situ data set (Di Cicco et al., 2017) derived from the SeaBASS data set (Werdell & Bailey, 2005).
Using the same methods applied in Ciavatta et al. (2018), we computed unbiased concentrations and associated standard errors of the PFT concentrations. These were then projected onto the model grid. Finally, we computed 5-day composites centered on the last day of each month of years 1997-2015. These composites were assimilated on a monthly time step, and the results of the reanalysis are presented in this work for years 1998-2014.
To corroborate the simulated sea surface temperature (SST), we used data of the Group for High Resolution Sea Surface Temperature (GHRSST) global Level 4 AVHRR_OI product. This product uses optimal interpolation (OI) by interpolating and extrapolating SST observations from different sources, which include satellite (Advanced Very High Resolution Radiometer [AVHRR]) and in situ platforms, i.e., ships and buoys (Reynolds et al., 2007). The original GHRSST global Level 4 AVHRR_OI data at daily temporal resolution and 0.25°spatial resolution were averaged to monthly values and interpolated to the model grid for comparison with the model output. The in situ biogeochemical observations used to evaluate the reanalysis skill were obtained from publicly available data sets. Temperature, salinity, nitrate, phosphate, and chlorophyll data were provided through OGS/NODC infrastructure (http://nodc.ogs.trieste.it); we used data of biogeochemical-ARGO profiles for nitrate and oxygen (http://doi.org/10.17882/42182#58265) as well as for salinity, temperature, and chlorophyll (BOPAD-prof, Barbieux et al., 2017).
In situ observations of PFTs were obtained from Di Cicco et al. (2017), in turn, derived from the SeaBASS data set (Werdell & Bailey, 2005).

Skill Assessment and Similarity Evaluation
The skill of the assimilation of the ocean color PFTs (reanalysis for 1998-2014) was assessed in comparison with a reference simulation without assimilation and in comparison with assimilated ocean color products and independent in situ observations, by using quantitative statistical metrics, as in Ciavatta et al. (2018). The performance of the reference and reanalysis output (y) in matching the assimilated composites of ocean color PFT concentrations (y′ c ) was evaluated by computing spatial correlation and RMSD at each month of the simulation years and then averaging in time. The skill of the reanalysis was evaluated for the output of both the forecast and analysis steps of the assimilation algorithm (see section 2.2 for the definition of these steps). In particular, the assessment of the analysis skill used the same ocean color products that were assimilated monthly, while the assessment of the forecast used ocean color data that had not been assimilated yet. Although this is a common procedure in the evaluation of assimilative systems (see, e.g., Gregg et al., 2009), we stress that the not-assimilated ocean color data used for skill assessment cannot be considered independent from the assimilated ones, because of temporal correlations among subsequent ocean color fields retrieved from satellites. Similarly, also the nonassimilated in situ PFT observations cannot be considered fully independent because they were used to calibrate the assimilated ocean color PFT products (Di Cicco et al., 2017). The remaining in situ physical and biogeochemical observations are independent from the assimilated products.
The similarity between the reanalysis and ocean color average phytoplankton community was synthetized in a map scaled from 0 to 1. This map was derived by (1) computing and summing the absolute differences between the reanalysis and ocean color maps of the four PFT average fractions; (2) scaling the sums between 0 and 1, by dividing by the maximum sum; (3) computing the differences in (1) above minus the scaled sums.
In the map, the value 1 represents the maximum similarity, while 0, the lowest.
Quantitative metrics were used to evaluate the skill of the reanalysis in matching independent in situ observations. Daily values of the variables in the reanalysis data set (y) were matched up point to point in space and time with the data (o). In particular, we computed and then showed in robust skill diagrams Ciavatta et al., 2018) (a) the bias calculated as the median value of the reanalysisto-observation mismatch, bias* = median(y − o), normalized by the interquartile range of the data (IQR o ); (b) the unbiased median absolute error, MAE′ = median{abs[y − o − bias*]}, normalized by IQR o , and taken with the algebraic sign of the differences between the interquartile range of the output and the data, sign (IQR − IQR o ); and (c) the correlation coefficient. In these diagrams, the closer the point is to the center, the lower the error of the simulation (i.e., the median bias and the unbiased absolute error).
The coherence between in situ observations and model output was evaluated also by comparing their spatial distributions in bar plots, representing median and interquartile range. For discussion purposes, the matchups were subdivided in broad zones, encompassing coastal waters (bathymetry less than 200 m), open surface ocean and open deep ocean (with bathymetry more than 200 m and surface/deep thresholds of 10 m for chlorophyll and 200 m for nutrients). The in situ chlorophyll observations have a detection limit of 10 −3 mg/L, and the same threshold was set to the matchup reanalysis data used in the skill assessment, for the sake of consistency.

Definition of the PFT Ecoregions and Carbon Fluxes
Three PFT ecoregions were defined by assessing the spatial variability of the reanalyzed MIC, NNP, and PIC fractions computed at the surface of the model domain, averaged in the years 1998-2014.
The PFT fractions were computed as ratios of each PFT chlorophyll to total chlorophyll. Diatoms and dinoflagellates fractions were summed up to define a region of dominance of MIC as a whole. The NNP and PIC regions were defined by weighting the ratios of these two PFTs (see Table 3 for details).
In order to characterize the different regions of the Mediterranean Sea, we have selected three simulated ecosystem fluxes that can be of interest for applications such as aquaculture or fisheries: (1) the ratio of gross primary production to community respiration (GPP:CR), which defines an autotrophic system when greater than 1; (2) the trophic efficiency, here approximated by the ratio between mesozooplankton and phytoplankton biomass (mesozooplankton is the highest trophic level represented in the model), which indicates how much of the primary production is transferred to the upper trophic levels (see, e.g., Kwiatkowski et al., 2019); and (3) the sedimentation efficiency, here defined as the flux of particulate organic carbon (POC) reaching the sediments, normalized by the GPP. All the fluxes, except the sedimentation efficiency, were computed by integrating the model output from the surface to the maximum depth of 200 m, which on average incorporates the entire euphotic zone, encapsulating the vertical distribution of phytoplankton in the Mediterranean Sea (Lavigne et al., 2015). Sedimentation efficiency was computed at the bottom of the water column and provides an indication of the POC available to the benthic filter feeders, for example, commercially relevant shellfish.

Skill in Matching the Ocean Color Products
The reanalysis reproduced the distributions of the PFT chlorophyll concentrations as inferred from ocean color observations in the surface waters of the Mediterranean Sea ( Figure 3). On average, PIC and NNP reanalysis had the highest chlorophyll concentrations and displayed east-west differences, though the gradients were more evident in the satellite products. The MIC group (i.e., both diatoms and dinoflagellates) had lower average concentrations but high peak values and steep gradients in coastal areas (e.g., in the Adriatic Sea and Gulf of Gabès). The model overestimated ocean color PIC, leading to an overall overestimation of total chlorophyll. This holds in particular for the eastern basin, where NNP and diatoms were overestimated as well.
Crucially, the simulated average distribution of the phytoplankton community structure matched well the one inferred from the ocean color products, as shown in the maps of The PFT community structure has evident gradients in both the model and data fields (see PFT fractions in Figure 4 and Figure S1 in the supporting information) that reflects the gradients of the simulated surface nutrients (see section 3.2). PIC ratio is higher in the more oligotrophic eastern Mediterranean than in the western Mediterranean. Diatoms (and dinoflagellates in the reanalysis) dominate in the nutrient-rich coastal waters of the Adriatic Sea, the Gulf of Gabès, and the area surrounding the Nile estuary. NNP relevance is intermediate among the other PFTs ( Figure 6): It is lower than PIC fractions in both the west and East Mediterranean according to both the satellite and reanalysis products (see Figures 4 and S1). However, the east-to-west and coast-to-open sea gradients are more evident in the satellite product that in the model simulation of the NNP and PIC ratios. Figure 5 shows that the reanalysis and ocean color representation of the plankton community structure tends to be more alike in the open ocean than in coastal areas, where the normalized similarity index might be lower than 75%. In coastal areas, often MIC (i.e., diatoms plus dinoflagellates) dominate in both the reanalysis and ocean color product but with rather higher fraction values in the ocean color product (e.g., in the NAd, Gabès gulf and Nile delta, Figure 4). Similarity is lower than 75% also in the North West Mediterranean, North Aegean and Ionian Sea, where the reanalysis NNP fractions are noticeably lower than in the ocean color products.

Comparison With the Reference Simulation
The assimilation of ocean color PFT improved the reference model simulation of the phytoplankton community structure (i.e., PFT fractions, Figure 7a). This holds for the diatoms and dinoflagellates fractions in particular, which halved their normalized bias. The bias of the PIC also decreased slightly. NNP skill did not 10.1029/2019JC015128 Journal of Geophysical Research: Oceans improve significantly, and their reference correlation remained negative, since the west-east gradient of the ocean color fractions remained unrepresented (Figure 4).
The reanalysis did not change significantly the reference skill for the PFT concentrations, excluding dinoflagellates that were slightly deteriorated (Figure 7b). This resilience of the skill is not related to a malfunctioning of the assimilative system but to feedbacks of the simulated ecosystem to the analysis improvements. In fact, the analyzes had better skill than the 1-month forecasts, as expected from a functional assimilative system, but the forecast skill remained comparable to the reference one ( Table 1). The EnKF analysis increments shifted the model toward the ocean color gradients, by increasing strongly the biomass in the West Mediterranean (NNP) and along the coasts (dinoflagellates and diatoms) and decreasing the biomass in the east (NNP and PIC; see Figure S2 in the supporting information). However, these changes led also to proportional EnKF increments of nutrient content  Figure S2), which became a source of nutrients to the system. This source of nutrients exacerbated the model nutrient bias and sustained basin-wide growth of total phytoplankton, despite localized assimilative corrections of single PFTs (e.g., MIC increase in coastal areas). Therefore, the analysis was able to improve the reference spatial

Skill in Matching In Situ Observations
The reanalysis was relatively skilled in fitting the not-assimilated observations of physical and biogeochemical variables (Figure 8 and Table 2) collected at the sampling sites shown in Figure 9. Median bias and absolute errors were smaller than the variability of the observations for all the variables, except for oxygen, as indicated by the circles within the unit radius in Figure 8. The interquartile ranges of data and reanalysis overlapped for all the variables (Table 2) and the correlations for temperature, salinity, and chlorophyll were as high as 0.92, 0.66, and 0.44, respectively (Figure 8). Correlations for the other variables were relatively low. Nitrate reanalysis and data overlapped, though the simulations overestimated the observed distribution (Table 2). Similarly, phosphate observations were overestimated in the reanalysis: The positive bias and rightside position in Figure 8 reflects the higher median value and larger interquartile range in Table 2. The phosphate overestimation likely explains the overall model overestimation of chlorophyll concentrations (Figure 3). In fact, phosphate remains on average the limiting factor of primary production in our simulation, in agreement with previous experimental and modeling findings (e.g., Krom et al., 1991;Lazzari et al., 2016), implying overproduction where this nutrient was overestimated. However, we note that phosphate data were available for model assessment only in limited areas of the Mediterranean Sea.
The reanalysis captured the highest observed PFTs concentrations well (see the comparable values of IQ3 in Table 2, with the exception of the underestimated NNP) but the simulated distributions tended to be skewed toward low values (Table 2), leading to an overall negative median bias (Figure 8). We note that the number of in situ observations available for the assessment of PFTs was relatively low compared to the other variables.
The spatial patterns of the errors were different for the different variables, suggesting that the coupled physical-biogeochemical model did not have systematic regional failures (Figure 9). Errors of simulated temperature were higher in the area of the Gulf of Lion and Ligurian Sea (Figure 9). In the same area, also a comparison with the GHRSST SST confirms a larger model uncertainty ( Figure S3 in the supporting information). This can be related to the model underestimation of the velocities of the Ligurian and Balearic currents and excessive summer stratification of the water column . The model simulated poorly other complex hydrodynamic features, such as the gyres in the North Tyrrhenian and South Adriatic and currents around the Sicily Island: See Figure 9 and Figure S3 in the supporting information. This uncertainty is probably a consequence of the coarse resolution of the model grid (0.1°) and of the meteorological forcing (0.75°) and has likely contributed to the bias of the model simulation of nutrients and PFT concentrations.
The error patterns of salinity were less systematic, since high and low error points coexisted in the same area; see, for example, the Algerian Basin and the Tyrrhenian Sea in Figure 9. An exception is the NAd, in proximity of the Po River, where the model simulated fresher waters than observed. Nitrate, chlorophyll, and dissolved oxygen also reached the highest errors in the NAd (Figure 9); this was a direct consequence of using approximate climatological riverine discharges.

Journal of Geophysical Research: Oceans
The model was skilled in simulating phosphate observations in the NAd (Figures 9 and 10). This mitigated the impact of nitrate overestimation on primary production in the NAd, given that phosphate was the limiting factor in our simulation. The model tended to overestimate the phosphate data in other coastal areas and in the open ocean ( Figure 10), but too few data were available for a robust assessment (Figure 9).
Chlorophyll error did not have clear spatial patterns, since high and low error data points coexisted in the same area (e.g., in the Algerian Sea, the Gulf of Lion, and the NAd, Figure 9). Figure 10 shows that the reanalysis slightly overestimated the in situ observations in the open surface water (depths 0-10 m), where the model tended to overestimate also the ocean color product (Figure 3). At the same time, the reanalysis distribution in the deep open ocean (depths 10-200 m) overlapped with the lowest tail of the observed distribution, meaning that the reanalysis tended to underestimate the deep chlorophyll maximum values. In the coastal areas, the simulated and observed distributions are comparable and both reached the highest values in the North Adriatic Sea. Here the model did not produce systematic chlorophyll outliers despite the nitrogen overestimation because simulated phytoplankton production was limited by phosphate ( Figure 10). If chlorophyll data from the NAd are removed from our comparison, the reanalysis and observation distributions in coastal areas were 1 order of magnitude lower ("No NAd" in Figure 10).
The simulated PFT concentrations overlapped with the observations, though the reanalysis was skewed toward low values (Tables 2 and S4). However, the relatively low number of sampling stations affects the robustness of the assessment (Figure 9).

Comparison With the Reference Simulation
The reanalysis skill in matching the in situ observations of the biogeochemical variables was not significantly different from the skill of the reference model without assimilation (Figure 8 and Table 2). For the PFT concentrations, assimilation reduced the reference bias versus in situ observations (in particular that of PIC), but it also increased their MAE′ slightly ( Figure 8). The reanalysis skill overlapped with that of the reference for nutrients and dissolved oxygen, while phosphate was degraded (Figure 8). The limited impact of PFT assimilation on fitting the in situ data of these biogeochemical variables was expected, because none of them was included in the EnKF analysis (see section 2.2.1), but could only be updated in the forecast model integration (see to this regard the discussion in Pradhan et al., 2019). The relatively low frequency of the assimilative analyses and the relative Figure 7. Skill of the reanalysis (circles) in representing the spatial distributions of (a) the ocean color phytoplankton functional type fractions and (b) concentration. The skill of the reference simulation is shown for comparison (squares). IQR = interquartile range of the model estimates; IQRo = interquartile range of the observations; MAE′ = unbiased median absolute error; Bias = median value of the model-to-data deviations. The notation of the variables is explained in Table 2. high error of the ocean color PFT data can also help explain this outcome (see Ciavatta et al., 2018). The degradation of phosphate can be related to the analysis increments of phosphorus content of phytoplankton, which acted as a source of nutrients to the system, exacerbating further the reference bias (see section 3.1.1 and Figure S2 in the supporting information). We note that deterioration of nutrients can occur when assimilating ocean color, in particular when issues in the model lead to systematic biases in both the phytoplankton biomass and the nutrient concentrations (see discussion in Gregg et al., 2009), such as overestimation of nutrients inputs and chlorophyll concentration in our case. For example, assimilation of ocean color chlorophyll worsened estimates of nitrate in the NAd in the recent simulations of both Tsiaras et al. (2017) and Teruzzi et al. (2018), in relation to the systematic underestimation of both chlorophyll and nutrients in that region. Though we consider the model suitable to assess the PFT data assimilation strategy in the context of this work (see the relatively good match of model and ocean color PFT fractions in Figures 4-7 and with PFT in situ observations in Figure S4), we highlight that a comprehensive assessment of both nitrate and phosphate fluxes in the model needs to be performed in future studies.

PFT Ecoregions in the Mediterranean Sea
The good skill of the reanalysis in estimating the PFT fractions inferred from ocean color (Figures 4, 5, and 6) and the improvements with respect to the reference simulation (in particular for MIC ratios, Figure 7a) supported the use of the reanalysis output in Figure 4 to define PFT-based ecoregions in the Mediterranean Sea (Table 3 and Figure 11).
The region of PIC dominance (blue in Figure 11) is the largest in the Mediterranean Sea by far. It covers the offshore waters from the East Mediterranean to the Balearic Islands. It also includes the Central Aegean Sea and the central and southern offshore areas of the Adriatic Sea.
The region of NNP relevance (yellow in Figure 11) covers large bands of the Mediterranean coasts, in particular in the Western Mediterranean, Adriatic Sea, and Levantine.
The region of MIC dominance (in Figure 11) is comparatively small and associated with coastal areas. It covers mainly the NAd and the Gulf of Gabès. However, it spreads also along the Italian Adriatic coast, the Libyan coast, and the Nile Delta.
In general terms, the simulated PFT distribution is largely driven by the gradients of surface dissolved inorganic nutrients available for growth (see the phosphate and nitrate distributions in Figure 9). MIC dominates  Table 2. Reference and reanalysis skill for T and S were identical because of the one-way coupling of biogeochemistry to physics in the model applied here. 7.9 (6.5-9.2) 7.9 (6.6-9.2) 5.7 (4.7-6.9) 84,537 Nitrate (mmol/m 3 ) NO3 5.2 (3.0-7.2) 4.9 (2.8-6.8) 2.9 (0.6-6. Note. "Not" is the notation of the variables used in Figure 8; Q1 and Q3 are the first and third quartiles, respectively; N is the number of matchups.

10.1029/2019JC015128
Journal of Geophysical Research: Oceans the areas near the main river mouths, where nutrient abundance can sustain relatively large cells, which require high intracellular nutrient concentration. When nutrient are not limiting, these PFTs tend to accumulate high amount of nutrients (luxury uptake) reducing their availability for smaller phytoplankton. Large phytoplankton is also less heavily grazed with respect to small species: In our model MIC is only preyed on by meso and microzooplankton, while small phytoplankton is preyed on by all the zooplankton functional types, including heterotrophic nanoflagellates. Conversely, the fast-growing, small PIC dominate the large areas of oligotrophic open waters of the Mediterranean. This PFT has a higher surface to volume ratio (due to their smaller size), and therefore, they are able to take up nutrients at lower concentration than large cells.
In the model, this feature is described through a parameter representing the affinity for nutrients during uptake ; for example, the PIC affinity for phosphorus is double that in diatom. Finally, NNP, which has intermediate characteristics with respect to the other two groups, is relevant The regionalization is based on a reanalysis that matched well the PFT structure inferred from ocean color, but additional considerations might be needed for those areas where the consensus between reanalysis and ocean color was lower. These areas are exemplified by the 75% similarity line in Figure 11, which was reported from the results in Figure 5.
We have confidence in our definition of the MIC region, although the region is included in the coastal areas where the similarity was lower than 75% ( Figure 11). In fact, the MIC-dominated areas were overlapping substantially in the reanalysis and ocean color products. The two products differed in the magnitude of the MIC domination, which was noticeably higher from the satellite algorithm than from the model, for example, in the MIC spots in the NAd and in the Gulf of Gabès (Figure 4).
Portions of the PIC region are also bounded by the 75% similarity line, most notably the highly productive Ligurian Sea/Gulf of Lion and the North Aegean Sea (Figure 11). Here NNP was dominant in the ocean color product, while PIC fraction was more than 50% in the reanalysis product ( Figure 4). Given the difference in the two products, the consideration of information from in situ PFT observations and local-scale modeling might be needed to confirm the classification of the Ligurian Sea/Gulf of Lion and the North Aegean Sea.
The Ligurian Sea/Gulf of Lion was previously classified as a "blooming" area, based on clustering analysis of the phenology of total chlorophyll from ocean color (d 'Ortenzio & Ribera d'Alcalà, 2009;Mayot et al., 2016). The ocean color product applied in our work indicates a climatological NNP blooming or dominance in the region from November to May, and PIC dominance from June to October; on the other hand, our reanalysis simulated NNP dominance only in the month of December (climatological analysis), and PIC dominance in the rest of the year (not shown). The reanalysis and ocean color dissimilarity in representing the PFT seasonal succession and average dominant PFT in the Ligurian Sea/Gulf of Lion is possibly related to our inadequate simulation of the Rhone river inputs (see the highest dissimilarity at the river mouth in Figure 5), also to the Table 3 Definition of the phytoplankton functional type (PFT) ecoregions in Figure 11 PFT ecoregion Definition MIC Microphytoplankton is dominant: Here microphytoplankton chlorophyll (sum of diatoms and dinoflagellates) ratio to total chlorophyll is higher than the remaining PFT chlorophyll ratios. NNP Nanophytoplankton is relevant: Here the nanophytoplankton ratio is the second highest ratio (after picophytoplankton) and summed with micro ratio is higher than the dominant pico ratio. PIC Picophytoplankton is absolutely dominant: Here the pico ratio is higher than 50%.
Note. The ecoregion were computed using the surface reanalysis PFT fractions shown in Figure 4.

10.1029/2019JC015128
Journal of Geophysical Research: Oceans misrepresentation of the complex hydrodynamics of the area (see section 3.2). In this region, higher-resolution models might be needed to resolve the hydrodynamic impact on the nutrient transport and zooplankton grazing on the phytoplankton community structure (Auger et al., 2014). However, we note that the satellite model used here was parameterized using high-performance liquid chromatography (HPLC) pigment data (Di Cicco et al., 2017). In the Atlantic Ocean, differences have been seen between size-fractionated chlorophyll estimated through HPLC and through size fractionated filtration , with implications for satellite size-fractionated chlorophyll and primary production algorithms (see . For Mediterranean waters, future efforts should be placed on collecting in situ data using other techniques (e.g., flow cytometry, microscopy, and size-fractioned filtration) alongside HPLC, to aid in calibration of future ocean color phytoplankton community structure algorithms.
In the North Aegean Sea, field data suggested that the area is overall dominated by PIC (>60%, Ignatiades et al., 2002;Siokou-Frangou et al., 2002), supporting our reanalysis-based classification ( Figure 11). Also, a local-scale ecosystem model confirmed the dominance by PIC in the North Aegean Sea (44%, Tsiaras et al., 2014), while MIC dominance and NNP dominance were driven by nutrient river discharges in confined coastal spots of the North Aegean Sea, for example, the Thermaikos Gulf identified as NNP also in our regionalization ( Figure 11).
Finally, the NNP region is bounded by the 75% similarity line in the Alborán Sea, Italian Adriatic coast, and South East Levantine ( Figure 11). Here the reanalysis and the ocean color are in general coherent in pointing out NNP in relevant fractions, supporting our classification, though the absolute magnitudes of the fractions can be quite different in the two products (see Figure 4). The different magnitudes are likely related to both the coarse resolution and approximate riverine inputs in our model. Because of that, our model cannot resolve the PFTs as accurately as locally tuned models, which could capture the ecosystem effect of eddies near the Gibraltar Strait (Ramirez-Romero et al., 2014), of the Po input transport  and of the Nile plume (e.g., Petihakis et al., 2009). Some in situ observations available in the above areas seem to support their classification in the NNP region. For example, relevant autotrophic NNP fractions were observed in the surface waters of the South East Levantine (comparable to the PIC biomass, Tanaka et al., 2007), in the Alborán Sea (up to 40% of the phytoplankton community, Arin et al., 2002), and in the west coast of the Adriatic Sea (e.g., the largest abundance fractions in Zoppini et al., 1995).
It is worth noting that the MIC and NNP ecoregions in Figure 11 are larger in space than the ecoregions obtainable from the reference simulation ( Figure S5 in the supporting information). If compared to ecoregions estimated from the ocean color product, the reanalysis area of NNP relevance is less extended ( Figure S6 in the supporting information) and comparable to the area where NNP is dominant in the ocean color product ( Figure S7 in the supporting information).

Regional Plankton-Driven Biogeochemical and Trophic Flows
The magnitude and variability of selected, normalized carbon fluxes, and trophic flows in the Mediterranean PFT ecoregions are presented in Figure 12. Our regionalization approach is supported by these results,  Table 3. MIC = microphytoplankton; NNP = nanophytoplankton; PIC = picophytoplankton. The ecoregions were derived from the surface reanalysis PFT fractions shown in Figure 4. The red line represents the similarity value of 75% derived in Figure 5.

10.1029/2019JC015128
Journal of Geophysical Research: Oceans because the vertical fluxes were relatively homogeneous within the regions and distinct between regions. In fact, the spatial interquartile ranges were relatively small and not overlapping in Figure 12, with some exceptions.
MIC is the most distinct region, having the maximal and more variable values of all the normalized fluxes.
Crucially, the NNP region also diverges from the PIC region in terms of trophic and sedimentation efficiency. That means "relevant" fractions of NNP can make a difference in fluxes in areas dominated by PIC. Apparently, this does not hold for the GPP normalized by the total respiration, since this ratio is not significantly different in NNP and PIC.
The regional differences of the simulated fluxes in Figure 12 can be explained based on our model formulation of the PFT traits. The differences are also determined by "emergent properties" of the ecosystem. These are features that are not coded explicitly in the model but arise from nonlinear synergies of biogeochemical and physical states simulated by the model, such as spatial shifts of the plankton trophic web in relation to nutrient and organic matter availability.
The community production-to-respiration ratio (GPP:CR) is higher than 1 on average, characterizing the epipelagic Mediterranean as an overall autotrophic system ( Figure 12). The ratio is higher in the nutrient-rich coastal areas dominated by MIC. In the model, large phytoplankton is preyed by mesozooplankton, supporting the herbivorous food chain and leading to higher transfer efficiency toward upper trophic levels than in regions where only smaller zooplankton are active. This trophic structure implies production of relatively large POC (from fast-sinking large organisms, fecal pellets, grazing, and mortality detritus), which leads to the highest sedimentation efficiency in the MIC region.
At the other extreme, in the oligotrophic PIC region, the median primary production is low and detritus is mainly directed to the dissolved pool (or to small, slowly sinking particles), explaining the smallest sedimentation efficiency. In this region small grazers are advantaged with respect to mesozooplankton, leading to the low efficiency of transfer toward upper trophic levels. These results are coherent with the findings of Allen et al. (2002) and by Lazzari et al. (2016), who estimated that 50% of the potential production is channelled toward the dissolved pool in the Eastern Mediterranean, almost entirely covered by our PIC region.
Normalized fluxes and flows in the NNP region sit between the two other regions (e.g., trophic and sedimentation efficiency), consistent with the intermediate traits of the NNP type (e.g., intermediate parameter values for the mesozooplankton predation and sinking rates in the model).
We argue that the simulated ecoregion efficiencies are meaningful for the real Mediterranean ecosystem because of the substantial correspondence between the reanalysis and ocean color PFT fraction distributions (Figures 4-7), despite some model biases and lack of data to corroborate the carbon flux estimates. In fact, the flows and fluxes in Figure 12 are all normalized quantities; thus, they should be affected only at the second order level by the biases in the phytoplankton biomasses shown in Figure 3. In addition, although flux data were not available to corroborate our novel PFT ecoregions, there is some consistency between flux estimates from our reanalysis and previous works at the scale of the whole Mediterranean Sea. Table 4 Figure 12. Normalized carbon fluxes in the reanalysis phytoplankton functional type ecoregions in Figure 11 (average values in 1998-2014). The height of the box represents the median value of the region, while the error bar represents the spatial variability within the region (i.e., the interquartile range). Fluxes are defined in section 2.5 and were integrated from the surface to the maximum depth of 200 m. The sedimentation efficiency is an exception, because it was computed at the bottom of the water column. The distributions of the carbon fluxes used in the computation of the efficiencies are shown in Figure S8 in the supporting information. CR = community respiration; GPP = gross primary production; MIC = microphytoplankton; NNP = nanophytoplankton; PIC = picophytoplankton.
shows that our estimates of the net primary production by the total phytoplankton community matches well the one by Uitz et al. (2012), who used observations from ocean color and in situ campaigns. The confidence bound for our estimates was derived from the range of our 100-member ensemble simulation, and it overlaps with the range provided by Uitz et al. (2012). We note that these estimates are lower than the ones derived from sparse in situ observations (Sournia, 1973) and from other model simulations (Di Biagio et al., 2019).
Our MIC net production is comparable to the one estimated by Uitz et al. (2012) based on their PFT ocean color algorithm, but NNP and PIC production diverge by~100%. This difference is possibly related to our overestimation of the PIC biomass, in particular in the North West Mediterranean area ( Figure 3). However, the PIC production might be underestimated by the ocean color algorithm of Uitz et al. (2012). The Uitz et al. (2012) algorithm is parameterized using pigment data. Differences have been seen in the Atlantic when using ocean color algorithms of size-fractionated primary production tuned using size-fractionated filtration-based methods , with higher estimates of PIC production (see Table 3 of . In addition, in situ studies have reported a 65% PIC contribution to the total Mediterranean production (Magazzù & Decembrini, 1995), comparable to 60% we estimated here. Reducing uncertainty in phytoplankton community structure in the Mediterranean Sea will ultimately require more in situ campaigns collecting information on phytoplankton community structure using a range of techniques (e.g., flow cytometry, microscopy, and size-fractioned filtration, alongside HPLC).
Carbon export to the sediment represents 8% of the net primary production in our simulation (Table 4). This value needs to be corroborated in future studies, since other estimates are not available for comparison at the scale of the Mediterranean Sea, to the authors' knowledge.

Conclusions
The Mediterranean Sea can be subdivided into three PFT ecoregions, which are homogeneous in the surface phytoplankton community structure as well as in selected, vertically integrated carbon and trophic fluxes: 1. a relatively small region of microphytoplankton dominance, in the most productive coastal areas, where the flows of food and POC toward the upper trophic levels and the sediments, respectively, reach the highest efficiency; 2. a large region of picophytoplankton dominance, covering mainly oligotrophic open ocean waters, which is the less efficient in transferring food and POC; 3. a region where nanophytoplankton is relevant, in coastal and Atlantic-influenced waters, where the efficiency of the carbon fluxes is intermediate.
The regionalization was based on a skilled reanalysis of the phytoplankton community structure of the Mediterranean Sea (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014). The reference simulation of the community structure (i.e., PFT fractions) was improved by the assimilation of a regional, quality-assessed ocean color PFT product, confirming the usefulness of this approach for the reanalysis of marine ecosystems (Ciavatta et al., 2018). However, the reanalysis did not improve significantly the simulation of in situ data of not-assimilated biogeochemical variables, which were not included in the multivariate analysis, and deteriorated the simulation of phosphate, as a consequence of the analysis leading to an increase in the phytoplankton nutrient contents. The

10.1029/2019JC015128
Journal of Geophysical Research: Oceans application of mass-conserving assimilative schemes preventing nutrient divergences could be useful to mitigate the issue (e.g., Hemmings et al., 2008), but model improvements to reduce the reference nutrient biases and thus increase the assimilation efficacy are certainly a priority. Such possible improvements include the coupling of higher-resolution physical models or reanalysis (e.g., Adani et al., 2011;Teruzzi et al., 2018), more realistic representation of riverine discharges and improved parameterization of the model PFT chlorophyllto-carbon ratios.
We acknowledge the limits of our regionalization, which is based on time averaged, surface PFT chlorophyll estimates, and therefore cannot represent the seasonal, interannual, and vertical variability of the plankton community structure. These choices were functional to provide a concise regionalization of the Mediterranean, corroborated by surface satellite data, and characterized by homogenous, vertically integrated ecosystem fluxes.
Our regionalization can add a useful layer of information for basin-scale planning of ecosystem studies, management actions, and sustainable exploitation of the Mediterranean Sea, for example, in the context of the European Marine Strategy Framework Directive (Ayata et al., 2018). For example, one could speculate that the MIC region identified here is the most suitable for fisheries or aquaculture activities sustained by large and efficient particulate organic matter flows (e.g., shellfish production). Then, local ecosystem, high trophic level, and aquaculture models (e.g., Brigolin et al., 2014;Petihakis et al., 2012;Politikos et al., 2015;Tsiaras et al., 2014) could be nested off-line to our reanalysis for fine-scale analysis and aquaculture planning within the sustainable macroregions identified in this work. Such kind of application is the object of our current research in the framework of the European H2020 TAPAS Project.