The Mediterranean subsurface phytoplankton dynamics and their impact on Mediterranean bioregions

. Ocean bioregions are generally deﬁned using remotely-sensed sea surface chlorophyll ﬁelds, based on the assump-tion that surface chlorophyll is representative of euphotic layer phytoplankton biomass. Here we investigate the impact of subsurface phytoplankton dynamics on the characterisation of ocean bioregions. The Mediterranean Sea is known for its con-trasting bioregimes despite its limited area, and represents an appropriate case for this study. We modelled this area using a high resolution regional dynamical model, NEMO-MED12, coupled to a biogeochemical model, PISCES, and focused our 5 analysis on the bioregions derived from lower trophic levels. Validated by satellite and Biogeochemical-Argo ﬂoat observations, our model shows that chlorophyll phenology can be signiﬁcantly different when estimated from surface concentrations or integrated over the ﬁrst 300m deep layer. This was found in both low chlorophyll, oligotrophic bioregions as well as in high chlorophyll, bloom bioregions. The underlying reason for this difference is the importance of subsurface phytoplankton dynamics, in particular those associated with the Deep Chlorophyll Maximum (DCM) at the base of the upper mixed layer. 10 Subsurface phytoplankton are found to signiﬁcantly impact the bloom bioregions, while in oligotrophic regions, surface and subsurface chlorophyll are of similar importance. Consequently, our results show that surface chlorophyll is not representative of total phytoplankton biomass. Analysis of the DCM ﬁnds that its dynamics are extremely homogeneous throughout the Mediterranean Sea, and that it follows the annual cycle of solar radiation. In the most oligotrophic bioregion, the total phytoplankton biomass is almost constant along the year, implying that the summertime DCM biomass increase is not due to DCM 15 photoacclimation, nor an increase of DCM production, but instead of the "migration" − with photoacclimation − of surface phytoplankton into the DCM. Integrated Chl bioregions are primarily deﬁned by a surface to subsurface signal balance, and present a progression from the bloom bioregion to the most oligotrophic one. The results illustrate the importance of subsurface dynamics − and of the Deep Chlorophyll Maximum in particular − in the Mediterranean Sea. A bioregionalisation of the DCM revealed a very homogeneous annual cycle of its chlorophyll content, with seasonal progression that evolves with the annual solar radiation cycle. Furthermore, in oligotrophic areas, stable phytoplankton biomass implies that the summer Chl increase in the DCM results from a phytoplankton "migration" from the surface layer to the DCM. This highlights the misunderstanding that comes from surface Chl analysis. The phytoplankton biomass does not decrease drastically in summer in oligotrophic bioregions because of nutrient limitation. Instead, it mostly goes deep to the DCM, where growth conditions are adequate with higher nutrients concentration and sufﬁcient light availability.

Correspondence: Julien Palmiéri (julien.palmieri@noc.ac.uk) Abstract. Ocean bioregions are generally defined using remotely-sensed sea surface chlorophyll fields, based on the assumption that surface chlorophyll is representative of euphotic layer phytoplankton biomass. Here we investigate the impact of subsurface phytoplankton dynamics on the characterisation of ocean bioregions. The Mediterranean Sea is known for its contrasting bioregimes despite its limited area, and represents an appropriate case for this study. We modelled this area using a high resolution regional dynamical model, NEMO-MED12, coupled to a biogeochemical model, PISCES, and focused our 5 analysis on the bioregions derived from lower trophic levels. Validated by satellite and Biogeochemical-Argo float observations, our model shows that chlorophyll phenology can be significantly different when estimated from surface concentrations or integrated over the first 300m deep layer. This was found in both low chlorophyll, oligotrophic bioregions as well as in high chlorophyll, bloom bioregions. The underlying reason for this difference is the importance of subsurface phytoplankton dynamics, in particular those associated with the Deep Chlorophyll Maximum (DCM) at the base of the upper mixed layer. 10 Subsurface phytoplankton are found to significantly impact the bloom bioregions, while in oligotrophic regions, surface and subsurface chlorophyll are of similar importance. Consequently, our results show that surface chlorophyll is not representative of total phytoplankton biomass. Analysis of the DCM finds that its dynamics are extremely homogeneous throughout the Mediterranean Sea, and that it follows the annual cycle of solar radiation. In the most oligotrophic bioregion, the total phytoplankton biomass is almost constant along the year, implying that the summertime DCM biomass increase is not due to DCM 15 photoacclimation, nor an increase of DCM production, but instead of the "migration" − with photoacclimation − of surface phytoplankton into the DCM.

Introduction
Phenology is the study of the occurrence and characteristics of recurrent biological natural phenomena (Menzel et al., 2006).
In marine science, phenology typically refers to the study of the annual cycle in phytoplankton biomass (Sverdrup, 1953;Longhurst, 2007), i.e. phytoplankton bloom starting date, the date of its maximum, bloom amplitude, bloom length, etc. These phytoplankton dynamics or regimes are strongly constrained by − and reflect − the complex interplay between ocean dynamics 5 (currents, temporal variability of the mixed layer depth, etc.), chemical environment (nutrients availability), and light forcing (Longhurst, 2007;D'Ortenzio et al., 2014). Therefore, phytoplankton phenology is representative of its surrounding ecosystem, and is considered as a good ecological indicator: a metric that provides objective information about ecosystem health, vigour and resilience (Platt and Sathyendranath, 2008), one that is particularly useful in regional comparisons and for characterizing the temporal evolution of the whole pelagic ecosystem. 10 Phytoplankton phenology is usually defined from the temporal evolution of chlorophyll-a concentration (Chl), the latter being a good proxy of phytoplankton biomass. But because of the scarcity of in situ measurements, basin-scale analyses of phytoplankton phenology are almost exclusively carried out by using Chl estimated by satellite ocean color sensors (Chl sat ). Indeed, the time frequency and spatial resolution of Chl sat are well adapted for basin-scale analysis. In particular, Chl sat is commonly used to identify ocean regions where patterns of phytoplankton phenology are similar − defined as bioregions − 15 representing areas of the ocean with (assumed) comparable ecosystem functioning (Platt and Sathyendranath, 2008;D'Ortenzio and Ribera d'Alcalà, 2009;Racault et al., 2012;D'Ortenzio et al., 2012;Sapiano et al., 2012;Mayot et al., 2016). The concept represents useful means for marine ecosystem state assessment (Platt and Sathyendranath, 2008;Racault et al., 2014).
However, satellite-based estimates have a number of limitations, principally the remote sensing only sees the surface layer 20 of the ocean. The first optical depth (i.e. the approximate depth seen from a space sensor) depends on water components (organic and inorganic), and it is generally limited to the first few meters (Gordon and McCluney, 1975;André, 1992). While satellite sensors can provide sufficient information for investigating phytoplankton variability in non-oligotrophic waters, in oligotrophic regions (or periods of oligotrophic conditions), the vertical distribution of the phytoplankton may not be directly related to the observed surface signature. In particular, the presence of deep chlorophyll maxima (DCM; Cullen (1982); Berland 25 et al. (1987); Varela et al. (1992); Cullen (2015)), ubiquitously observed in oligotrophic regimes, have been shown to coincide in the Mediterranean Sea with maximum phytoplankton biomass that is unseen in surface Chl measurements (Mignot et al., 2014). However, some indirect methods exist that permit the vertical profile of phytoplankton biomass to be estimated (Berthon and Morel, 1992;Morel and Maritorena, 2001;Uitz et al., 2006).
Hence, phytoplankton phenology − as well as ocean bioregions − estimated from Chl sat only reflects surface phytoplankton 30 dynamics, and may not necessarily be representative of phytoplankton biomass dynamics, especially in oligotrophic conditions.
In this paper, we investigate the phytoplankton dynamics of the oligotrophic Mediterranean Sea (Figure 1), because it is considered a hot spot for climate change (Giorgi, 2006;The MerMex group: Durrieu de Madron et al., 2011;Diffenbaugh and Ribera d'Alcalà, 2009;Mayot et al., 2016) that are well-known for their DCM (Berland et al., 1987;Crombet et al., 2011;Mignot et al., 2014;Lavigne et al., 2015). Moreover, the recent release of Biogeochemical-Argo floats (BGC-Argo floats) (Figure 2-a; Schmechtig et al. (2015)) that measure biogeochemical as well as physical variables into the Mediterranean (including Chl), presents a new avenue for insights into the occurrence and evolution of DCM in this sea (Mignot et al., 2014;Lavigne et al., 2015). In particular, DCM are observed across the whole Mediterranean sea in summer (Figure 2-b), even in the north-15 western region where they disappear in Winter. In general, DCM deepen eastward, from ∼40 m up to 150 m in the Levantine sub-basin, and vary on an annual cycle with a deepening in spring-summer, followed by a shallowing in autumn that tracks the depth change of an isolume (the level where the daily integrated photon flux is constant; Letelier et al. (2004);Mignot et al. (2014); Cullen (2015); Lavigne et al. (2015)). Furthermore, variations of chlorophyll concentration within the DCM are attributed to a change of phytoplankton biomass, and not simply to a photoacclimation process that only changes chlorophyll 20 (Mignot et al., 2014). Overall, DCM in the Mediterranean Sea are characterised by patterns that are already observed at global scale . The context, and these key features − the high diversity in phytoplankton dynamics in a small area; the consistent presence of a well documented DCM − serve to make the Mediterranean Sea an interesting domain to investigate the importance of subsurface dynamics for phytoplankton phenology, and its influence in the definition of the Mediterranean ecosystems. 25 Nonetheless, in situ observations of the subsurface habitat and the processes there remain sparse in both space and time at this time, and preclude a comprehensive regional analysis. Recent ecological modelling studies in the Mediterranean Sea have begun to successfully represent DCM (Lazzari et al., 2012;Macías et al., 2014). Modelling then represents an attractive alternative to explore questions concerning DCM. 30 In this study, we use a high resolution coupled dynamical-biogeochemical model (Aumont and Bopp, 2006), in a specific high resolution configuration for the Mediterranean Sea (Palmiéri, 2014), to investigate the influence of subsurface dynamics on the phenologies and bioregions of this sea. This study contributes to identifying the different ecosystems of the Mediterranean Here, we first present an observational validation of the model's performance, including the application of the bioregionalization procedure developed by D'Ortenzio and Ribera d'Alcalà (2009) to compare it with satellite estimates. Then, we use our simulation to investigate the characteristics of the DCM, and evaluate the impact on Mediterranean phytoplankton phenologies 5 and bioregions.

The Mediterranean biogeochemical model: PISCES-MED12
The model used is the coupled dynamical biogeochemical NEMO-PISCES, in a high resolution (1/12 • ) regional configuration 10 developed for the Mediterranean Sea, MED12 (Palmiéri, 2014;Richon et al., 2018c, b). PISCES is a biogeochemical model originally developed to represent the global-scale biogeochemistry, though it has previously been adapted for use in more regional studies (e.g. the Arabian Sea; (Resplandy et al., 2009(Resplandy et al., , 2011). PISCES is an intermediate complexity biogeochemical model (Aumont and Bopp, 2006), that represents 2 size class of phytoplankton, zooplankton, and particulate matter, for a total of 24 tracers. The phytoplankton growth in PISCES is limited by the availability of nutrients (nitrate, ammonium, phosphate, 15 silicic acid, and iron). Also, uptake and release of nutrient proceeds accordingly to the Redfield ratio (O/C/N/P = 172/122/16/1; Takahashi et al. (1985)).
The NEMO-PISCES simulations here have been generated using specific Mediterranean datasets for initial and boundary conditions: i) The model has been set-up on a 1/12 • resolution grid (eddy permitting), that covers the entire Mediterranean Sea, plus a 20 buffer zone west of the Strait of Gibraltar (see Figure 1). The dynamical model NEMO-MED12 (Beuvier et al., 2012b) was driven at the surface, using 53 years of atmospheric forcing from the ARPERA model (Herrmann and Somot, 2008), based on the ERA40-ERAInterim reanalysis for the 1958-2011 period (Uppala et al., 2005). Riverine water fluxes were prescribed from the Ludwig et al. (2009) inter-annual database. Initial temperature and salinity fields for the Mediterranean Sea were derived from the MEDATLAS II climatologies (MEDAR/MEDATLAS-Group, 2002). In the buffer zone west of the Strait of Gibralter, 25 these variables, as well as sea surface height, were initialised by -and then strongly relaxed toward (with a timescale of 3 days) -the World Ocean Atlas (WOA) climatologies (Locarnini et al., 2006;Antonov et al., 2006). The resulting circulation fields have previously been evaluated using simulated CFC (Palmiéri et al., 2015) and 3 He-3 H (Ayache et al., 2015) concentrations.
ii) The biogeochemical variables have been initialised with the Mediterranean SEADATANET climatologies (Schaap and Lowry, 2010), and from observational values (e.g. vertical cruise sections or single average values) where a spatially-resolved 30 climatology was not available (i.e. dissolve inorganic carbon, DIC; alkalinity; iron). In the buffer zone, WOA (Garcia et al., 2006a, b) and GLODAP (Key et al., 2004)  same river mouths that for the riverine fresh water in the dynamical model, and are adapted from the same dataset (Ludwig et al., 2009). No atmospheric input of nitrogen or iron (as dust) was prescribed in this simulation. An 1/8 • version of this model (PISCES-MED8), also exists for climate change studies on the Mediterranean sea (Richon et al., 2018a).
The biogeochemical model was run offline following the same procedure as Palmiéri et al. (2015) and Palmiéri (2014). In those, previously simulated ocean dynamics of the NEMO-MED12 model (53 year period from 1958 to 2011 as previously

Remote sensing fields
Evaluation of modelled surface chlorophyll and calculated phenology and ecoregions makes use of satellite-based chlorophyll estimates from Bosc et al. (2004). This product is a correction of the SeaWiFS surface chlorophyll estimates specifically 15 developed to better represent the low surface concentrations in the Mediterranean Sea (Bosc et al., 2004). However, while decreasing considerable SeaWiFS biases, this product still over-estimates in situ chlorophyll measurements, in particular in the most oligotrophic parts of the Mediterranean Sea. This dataset is available for a 8 years period, from November 1997 to October 2005. 20 We used the recent Biogeochemical-Argo data (BGC-Argo; Schmechtig et al. (2015)) to help evaluate the modeled chlorophyll.

Biogeochemical-Argo floats
BGC-Argo are floats transported by the current like lagrangian particles. They sample vertical profiles every 8 days, down to 1000 m depth, of temperature, salinity, nitrate, dissolved oxygen, chlorophyll, Photosynthetic Active Radiation (PAR), Coloured Dissolved Organic Carbon (CDOM), and particulate optical backscattering (bbp). In this study, only chlorophyll is used to validate the model. The observed period with the BGC-Argo (2013-2018) being different from the one simulated 25 by the model , and also because we are not interested in single profiles − but rather to validate our model with climatology on an area basis − we preferred not to perform point to point comparison for now. Instead BGC-Argo data have been gathered within 5 regions (see fig. 2), where 1) there is enough data for a reasonable data-model comparison, and 2) the regions have been chosen to be relatively close to surface chlorophyll derived clusters from D'Ortenzio and Ribera d'Alcalà (2009). For instance: 1-Liguro-Provençal area corresponds to the Bloom-Intermittently cluster; 2-Algerian area corresponds 30 to the #3 No-Bloom cluster; 3-Ionian area corresponds to the Intermittently cluster; 4-Levantine area to the #1 No-Bloom cluster; and 5-Tyrrhenian area that is a mix of No-Bloom -Intermittently clusters. All chlorophyll data have been averaged per area into monthly climatological profiles, to be easily compared to the model outputs.

Bioregionalization
Bioregionalization follows the method of D'Ortenzio and Ribera d'Alcalà (2009), in which bioregions are defined by applying a k-means clustering algorithm to normalized Chl sat annual cycles. We performed the same procedure, with some adjustments, to both modelled (Chl surf ) and satellite estimates (Chl sat ) for the period November 1997 to October 2005: i) Satellite estimates (Chl sat ) were monthly averaged to have the same frequency as model outputs. ii) Since available Chl sat are known to have artificially high chlorophyll concentration along the coast (due to high concentrations of coloured organic matter that confound the ocean colour algorithm), coastal areas are filtered from remote sensing estimates. In order to be consistent and avoid discrepancies between model and satellite estimates, filtering used model coastal area (bathymetry shallower than 160 meters) for both Chl sat and Chl surf .
Finally, as we are primarily interested in the seasonal cycle of chlorophyll, we normalized all values towards specific local 10 annual maxima (i.e. if we consider the surface chlorophyll for example, each grid cell's surface seasonal cycle is normalized to its surface annual maxima; same procedure apply when considering the DCM and the 0-300m vertically-integrated chlorophyll).
Next, the resulting normalized patterns of phenology were classified into different clusters, using a specific K-mean function ("kmeansruns" from the "fpc" R-CRAN package) to identify different clusters in the Mediterranean Sea, each of which corre-15 sponds to a specific chlorophyll annual cycle (and is characteristic of a particular Mediterranean ecosystem). The membership of each grid cell is then mapped (regionalization) to identify regions of the Mediterranean Sea with equivalent ecosystems: the so-called bioregions.
This clusterization procedure is performed initially on surface Chl, but is additionally applied on the basis of the vertical maximum of Chl and the vertically-integrated Chl. This may potentially result in the diagnosis of different bioregions. 20 The number of significant clusters is defined by analysing the stability of each cluster, and keeping the maximum number of clusters for which each resulting cluster is stable. For a fixed number of cluster, 100 "clusterization" experiments are processed.
Each experiment is based on data randomly disturbed by maximum 5%. Three different disturbing method are used: Bootstrap, Noise, or Jitter (see Hennig (2007) for further details). The Jaccard similarity coefficient of a cluster (ratio of the number of points common to both disturbed and original cluster, compare to the total number of points that belong to, at least, one of 25 the 2 clusters) is calculated for each cluster, on each experiment. A cluster is then considered stable if its averaged Jaccard coefficient over the 100 experiments is ≥ 0.75 (see for example table 2).  Morel and Gentili, 2009). This is also the case for our analysis, where low Chl values are especially overestimated (Bosc et al., 2004). Switching to temporal performance, the seasonal cycle is generally well simulated by the model (Figure 3), which is an imperative characteristic for conducting phenology analysis. In the Western basin, the bloom occurs in February in the model, slightly earlier than satellite estimates that present a maximum in March. Both model and satellite show minimal Chl values in summertime (June to September), then Chl surface concentrations rise again from September, with a slower 15 growth rate in the model. The amplitude of the seasonal cycle is less important in the Eastern basin. The model produces a broadly realistic seasonal cycle with the minima in summer, occurring however slightly later, and for a longer period, than in satellite estimates. There is also a slower and later increase toward the winter maximum values, indicating some difficulty for the model in simulating autumn phytoplankton growth. As noted previously, in this study we will pay attention to the influence of subsurface Chl. Thus, we now evaluate our model results with in situ observations on vertical profiles. In this purpose, we , what is all − qualitatively − in good agreement with the observations. However this comparison also highlights that the model underestimates the Chl concentration, not only at the surface, but also at depth by ∼60%. It suggests that the model underestimates the chlorophyll in general (see the appendix A, for further evaluation and discussion on this point). Furthermore DCM depth is always deeper than observed by approximately ∼30 to ∼50 m. Finally, despite some shortcomings compared to observations, the model generally produces a global Chl distribution with realistic temporal and spatial structures. Notwithstanding the mismatches in absolute chlorophyll identified above, the overall performance of the model is sufficient to permit an investigation of the general patterns of chlorophyll phenology in the Mediterranean Sea (see discussion, section 4.2).

Surface chlorophyll phenologies and bioregions 5
Applying the clusterization procedure to both Chl sat and model Chl surf results in different total cluster numbers (see Table   2). Model Chl results in 4 clusters, while satellite Chl results in 5. For parsimony, we then use the common 4 clusters in the rest of the study for both model Chl surf and Chl sat in order to enable further comparisons.
The 4 clusters identified from the remote sensing estimates (based on a monthly dataset; see When applied to the model outputs, the clusterization procedure produces analogous results to those of the satellite estimates 30 (see Figure 5-b). It generates the main structures − one Bloom-Intermittently bioregion (red), and 3 No-Bloom − with similar annual cycles that have realistic winter-spring maxima and summer minima in chlorophyll. However, despite these robust similarities, some differences obviously exit. For instance, the amplitudes of the seasonal cycles tend to be more elevated than those deduced from Chl sat , due to overly oligotrophic conditions in summer in the model (Palmiéri, 2014 Chl surf increases later, due to a nutricline in the model that is too weak, especially in the western basin (see the appendix A).
Model blooms also occur slightly earlier (February-March) than those observed (April).
The clusters' spatial structures also reveal some differences. The modelled Bloom-Intermittently bioregion is satisfyingly located in the Gulf of Lions, but it extends to the east, taking in the Northern Ionian, the Southern Adriatic, and the Levantine sub-basins, while it is excluded from the Tyrrhenian sub-basin. The No-Bloom bioregion observed off the Maghrebin coasts in In summary, despite some unavoidable differences, when applied to the model surface chlorophyll, the clusterization pro-10 cedure generates results that are consistent with those from satellite estimates. The model produces coherent bioregions, with realistic chlorophyll annual cycles/phenology from surface chlorophyll fields, even if their respective geographical location may differ from those of satellite estimates. Bearing in mind the limitations associated with these model-observation discrepancies, using the model to study the phenology and bioregions of the Mediterranean Sea remains plausible.

15
Surface chlorophyll is commonly used as a proxy of phytoplankton biomass, a key facet of interest for biogeochemists. However, to our knowledge, very few studies link the surface annual Chl signal to the integrated Chl for the whole epipelagic layer − and these are still early stage studies (Uitz et al., 2006;Alvain et al., 2008;Uitz et al., 2012). In part, as implied above, this stems from the limited availability of appropriate observational data. Here, we take advantage of complete model fields to address the issue of the relationship between surface and deep chlorophyll. Applying the clusterization approach to chlorophyll 20 integrated across the 0 − 300 m deep layer (Chl tot ; a layer wide enough to include the euphotic layer). produces bioregions and corresponding annual cycles ( Figure 6) that are very different from those obtained using Chl surf . The procedure generates the same number of stable clusters (Table 3), but there are a number of significant differences. In particular, the temporal structure of all clusters is significantly modified, with annual cycle amplitudes largely reduced, and the appearance of 2 local maxima.
The first one occurs in February-March, as with Chl surf , but an additional maxima now manifests in summer. Such large mod- 25 ification can only be attributed to subsurface chlorophyll. Furthermore, unlike surface chlorophyll, for which annual minima happen in summer, total chlorophyll minima now occur in November-December across all bioregions. The annual cycles of all Chl tot regimes present a progression from bioregions where the seasonal amplitude is relatively high, and where the late winter-early spring signal is dominant (in red), through to bioregions where the amplitude is very low, and where the summer signal is slightly dominant (in blue). 30 Furthermore, the locations of the Chl tot bioregions differ from those of Chl surf . These locations also present the progression highlighted in Chl tot phenology, from the red to the blue cluster. Interestingly, the locations of these only 2 extreme bioregions are very similar to the 2 Chl surf extremes: the red Bloom-Intermittently and the yellow No-Bloom (Figure 5-b). But the red Chl tot bioregion's area is less expended than from the surface, and the blue Chl tot does not include the Aegean sub-basin.  (see table 4). One (in 20 green) covers the western basin plus the areas of the eastern basin where Chl surf has its red Bloom-Intermittently bioregion (see the red cluster in Figure 5-b); and the second one (yellow) extends across most parts of the eastern basin. The annual cycles differ completely from those of the corresponding surface bioregions. The dynamics of Chl max follow the seasonal solar radiation cycle, with a minimum in winter (January or February) and a maximum in summer (June or July), suggesting that phytoplankton growth in DCM is mainly controlled by solar energy rather than nutrient availability. Differences between 25 the "Western" and "Eastern" clusters reside mainly in the amplitude of the Chl max annual cycle, which is less important in the "Eastern" cluster (yellow), due to its more oligotrophic conditions. The amplitude of the western cluster is more important because this region benefits from higher nutrient supply, due to mixing with rich intermediate and deep water, and to the inflow of nutrient-rich Atlantic Water at the surface. These higher nutrient concentrations enable a shallower DCM and a higher receipt of solar energy by the phytoplankton, both of which result in higher chlorophyll concentrations.

Bioregions and Mixed Layer Depth
The Bloom-Intermittently cluster made from Chl Surf (Figure 5) is the most recognizable due to its specific phenology and locations. Modelled phenology is correct for this cluster, but its locations differ from that derived from remotely sensed Chl Surf .
Spring-blooms are induced by the high amount of nutrients introduced to the euphotic layer by mixing, from deep and interme-5 diate waters. Furthermore, the interplay between the Mixed Layer Depth (MLD) and euphotic depth, can induce a temporal lag between the establishment of conditions favouring blooms, and bloom occurrence (Sverdrup, 1953;Dutkiewicz et al., 2001;Lavigne et al., 2013;D'Ortenzio et al., 2014). As MLD is a key factor for spring bloom occurrence, the modelled MLD may explain the differences between the satellite and the model Bloom-Intermittently cluster. Applying a regionalization procedure (with 4 clusters like with chlorophyll) to both a Mediterranean climatology (Houpert et al., 2015) and NEMO-MED12 MLD (see Figure 9), results in a cluster (in red) located exactly in the same regions that the 15 Bloom-Intermittently chlorophyll cluster. This red MLD cluster is associated with the Mediterranean blooms in both model and satellite estimates, and presents exactly the same differences between the model and the observations. It confirms that the differences seen between the model and the satellite Bloom-Intermittently bioregions are caused by differences in the circulation model dynamics.
Differences within the other Chl surf clusters between satellite and model results ( Figure 5) are more difficult to investigate. 20 At first glance, the physical dynamics do not appear to be the main determinant factor. MLD clusters in the model and in the climatology (except the red one discussed immediately above) are very similar, and hence cannot explain the differences in the No-Bloom clusters. Thus, these differences may be the effect of other factors highlighted in PISCES-MED12 results like a too weak nutricline in western basin, or too low phosphate concentrations over the first 100 m depth of the Mediterranean sea (see Appendix A and Palmiéri (2014)). 25

Modelled chlorophyll maximum and total chlorophyll phenologies
In order to evaluate in more detail the simulated Chl phenology, we compare it to the novel information gained from the recently released BGC-Argo floats (Schmechtig et al., 2015)). We averaged BGC-Argo Chl profiles to get a monthly averaged annual cycle in the 5 areas defined in figure 2, and normalized them by their respective annual maximum. In each area, we selected column - Figure 11). First, across all selected areas it produces a surface chlorophyll signal (in green) that is comparable with that estimated from satellite ( Figure 5). It gives confidence to our procedure for treating the information from the Argo floats.
The difference with simulated Chl surf confirms the model difficulties to reproduce the autumnal surface chlorophyll increase (Figure 3), which appears too late, and too weak, because of a generally low PO 4 surface concentration (see the appendix A).
The vertical Chl maxima obtained from the Argo float observations (in red) is globally similar to that produced by the model.

5
Its phenology presents more variability, but the general Chl max dynamics confirm a maximum in summer and minimum in winter, with the exception of the Algerian sb, where Chl max is quite flat all over the year at the exception of a low value in December (but this might be due to the small number of observations in this sub-basin, or reflects sub-basin peculiarities such as its intense meso-scale eddy activity). The second exception is in the Gulf of Lions, where the spring bloom significantly dominates the Chl max phenology with an amplitude of 0.8 in the Liguro-Provençal sub-basin (what is not seen in the model, 10 as − even if the bloom is modeled, and the bloom cluster appears − its amplitude is largely underestimated), that compares to 0.4−0.5 in other areas. These apparent Chl max discrepancies could possibly be attributed to the non-homogeneity in BGC-Argo data, itself caused by the relatively short time records to date of the Lagrangian floats, and/or by the non homogeneity of the lagrangian floats sampling method itself. However, the elevated Chl max during the bloom period in the Liguro-Provençal sub-basin is most likely a real signal.

15
Despite the BGC-Argo float data being at an early stage, it still provides some interesting information on phenology that is reproduced by our simulation. For instance, Chl tot phenology is different everywhere from that of Chl surf . Also, the sub-basins have a range of different seasonal regimes, with a progression from a surface-dominated regime with high phenology amplitude (0.7) in the Liguro-Provençal sub-basin through slightly surface-dominated situations in the Algerian, Tyrrhenian and Ionian sub-basins (with interesting differences concerning winter to summer local maxima), to a more equilibrated surface/subsurface 20 regime with low phenology amplitude (0.2) in the Levantine sub-basin (i.e: in the most oligotrophic area). As already noted, all of these characteristics derived from the BGC-Argo profiles are in agreement with our modelling results. Only differences compare to model Chl tot are likely to be related to the underestimated surface chlorophyll in the western basin. The shape is good in the Liguro-provençal sub-basin, but the dominance of the bloom is underestimated. As well, the maximum Chl tot is in summer in the Algerian and Tyrrhenian sub-basins, with a too flat Chl tot phenology (amplitude of ∼0.2 compare to ∼0.4), 25 when surface signal should be the slightly dominating one.

Underestimated chlorophyll and DCM depth.
This study has shown that the model underestimates the chlorophyll at both the surface, and the DCM, and that it simulates a DCM that is too deep (see discussion of the later issue in the appendix A). How can this affect the Chl tot phenologies and resulting bioregions? 30 First, we use normalized chlorophyll, as done by D'Ortenzio and Ribera d'Alcalà (2009) and Mayot et al. (2016), in part because we are interested in the chlorophyll variation along the year, but also to avoid chlorophyll bias to impact the results.
Regarding Chl tot , it won't have any consequences, as far as the bias is the same in the whole water column, that Chl surf and Chl max are given the same weight in their contribution to Chl tot , than observed. That is the case in the Eastern basin: There, Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-423 Manuscript under review for journal Biogeosciences Discussion started: 8 November 2018 c Author(s) 2018. CC BY 4.0 License.
Chl surf /Chl max is 0.3 on average in BGC-Argo data, and 0.28 in the model (table 1). In the Western basin, we know it is not the case. Chl surf contribution is underestimated, with an observed Chl surf /Chl max of 0.45, and of 0.28 in the model (table 1).
It explains why the modelled Chl tot 's surface signal is not dominant in the Algerian and Tyrrhenian sub-basin, and not enough dominant in the Liguro-Provençal. That may have an impact on the bioregion definition of the Western basin.
The too deep DCM is less problematic for Chl tot . As we are interested of the monthly variation of the vertical chlorophyll 5 sum, as far as the Chl surf / Chl max balance is respected (as discussed above), the depth of Chl max will not have any impact on Chl tot phenology, or the induced bioregions. The definition of Chl tot is independent of the chlorophyll depth. It is still important to notice and investigate as it highlights a default in the simulation (here the nutricline is too smooth, see the appendix A), but a better CHL max depth would not improve Chl tot .
As the number of BGC-Argo floats is still regularly increasing in the Mediterranean Sea, the spatial and temporal coverage 10 of the float database will − probably within a few years − be large enough to derive an appropriate 3D climatology of the depth distribution of Mediterranean chlorophyll. This potentially represents an important alternative for investigating ecosystems of the Mediterranean Sea compared to both in situ observations and satellite-estimated surface chlorophyll.

Phytoplankton dynamics in the oligotrophic bioregion.
The preceding sections highlighted that the bioregions and phenology defined from Chl surf and Chl tot differ, and that these 15 differences are strongly impacted by the behaviour of the DCM. The DCM results from the interplay between light and nutrient availability. It is situated at the top of the nutricline, where the balance of downward photon flux and upward nutrients flux is optimal (Mignot et al., 2014;Cullen, 2015;Lavigne et al., 2015). The DCM experiences a seasonal vertical displacement, following the isolumes (Letelier et al., 2004;Mignot et al., 2014;Cullen, 2015;Lavigne et al., 2015) that deepen in summer with solar radiation annual cycle, and also with the summertime disappearance of the surface chlorophyll-induced shadow 20 effect. Mignot et al. (2014) suggest that the DCM is a result of photoacclimation processes (an increase of the phytoplankton Chl/C ratio when available light decreases) that result in a higher Chl/C ratio in the DCM than in surface waters (this has been confirmed by a modelling approach Ayata et al. (2013)). However, they also revealed that, inside the DCM, the seasonal Chl max concentration variations do not result from a further photoacclimation process (as DCM is supposed to follow an isolume) but rather result from changes in the phytoplankton biomass inside the DCM. 25 This photoacclimation process occurs in our simulation since the biogeochemical model PISCES includes the photoadaptative model of Geider et al. (1997). It is illustrated when analyzing the phytoplankton biomass bioregionalisation ( Figure   10), which results in 4 clusters ( Table 3) that are extremely similar to those from Chl tot (Figure 6). The only difference is due to phytoplankton photoacclimation. In particular, in the most oligotrophic region (blue cluster), the phytoplankton biomass is almost constant from February to August, while Chl tot reveals a Chl increase on the same period. Hence, this difference does 30 not result from any global increase of the phytoplankton biomass in this bioregion, but rather from a vertical "migration" − not in sense of individual movement, but rather understood as the "movement" of favourable growth conditions − and gathering of the phytoplankton − with its related photoacclimation − into the DCM.

Surface versus total chlorophyll bioregionalization
Unsurprisingly, surface and total chlorophyll phenologies and bioregions are different. In particular, phenologies highlight the fact that both surface − maximum peak in late winter to early spring − and subsurface chlorophyll dynamics − maximum peak in summer − are discernible and have a significant impact on Chl tot phenology across all bioregions. Surface Chl is characterized by a minimum in summertime and a maximum in late winter-early spring that are driven by physical dynamics 5 and the related winter vertical mixing that brings nutrients to the surface. In contrast, the subsurface signal is mainly controlled by the annual cycle of solar radiation that generates maximum Chl concentrations in the summer season. The integrated signal (Chl tot ) reveals characteristics that results mainly in a combination of these two modes of variability, with characteristics of both. This result was expected for the most oligotrophic part of the Mediterranean Sea, but the model shows that Chl tot phenology is also different from the Chl surf one in the Bloom-Intermittently bioregion. 10 Most of the ocean biogeographical studies are based on Chl surf , with the hypothesis that the surface chlorophyll is a good proxy of the total phytoplankton biomass (D'Ortenzio and Ribera d' Alcalà, 2009;Racault et al., 2012;D'Ortenzio et al., 2012;Sapiano et al., 2012) the latter being of greatest interest for biogeochemistry. To our knowledge, only one study (Platt et al., 2009(Platt et al., , 2010 compared surface chlorophyll phenology to phytoplankton biomass, with both being estimated by remote sensing (the latter calculated using an empirical model estimating phytoplankton biomass from Chl sat ). It highlighted that new features, 15 unseen by the Chl sat appeared on the phytoplankton phenology, in particular in low latitudes (Platt et al., 2009). However, the comparison has not to date been extended to biogeography comparisons.
When using Chl tot , the bioregions definition changes from bioregions defined on surface characteristics, to bioregions defined on a surface-to-subsurface Chl balance. The main reason for these differences being that subsurface productivity is shown to be important, principally in oligotrophic regions (as expected) but also in the bloom areas of the Mediterranean sea. Biore-20 gions based on Chl tot generate a more continuous bioregionalization, with an intermittently bioregion (yellow) now marking the transition between the bloom bioregion (red; strong surface signal) and the No-Bloom bioregions (green and blue, the most oligotrophic; equivalent amplitude in both surface and subsurface Chl signal). Our results show that in the Mediterranean Sea − a domain with a wide range of phytoplankton dynamics − the total chlorophyll is uniformly much closer to phytoplankton biomass ( Figure 10) than it is to surface chlorophyll, with significant consequences for the results and induced interpretations. 25 Interestingly, the 2 extreme bioregions (Bloom-Intermittently and the most oligotrophic) are the only ones that present reasonably close clusters location, within both Chl surf and Chl tot . This could have different reasons. 1-The surface signal is sufficiently dominant to keep being the main signal on Chl tot . This might be plausible regarding the Bloom-Intermittently bioregion, but not for the most oligotrophic one. Or 2-in these bioregions, surface and subsurface evolve together, and are more strongly correlated, getting a coherent and geographically identifiable cluster from either surface or vertically integrated 30 chlorophyll signal.
Thus, our results show that the surface chlorophyll alone is not a good proxy of total phytoplankton biomass. This does not, of course, mean that the surface chlorophyll is incorrect, but merely that it is incomplete: subsurface features that can completely change the interpretation of phytoplankton dynamics, and hence of the epipelagic ecosystem, are not taken into account (we 14 Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-423 Manuscript under review for journal Biogeosciences Discussion started: 8 November 2018 c Author(s) 2018. CC BY 4.0 License. know since Reygondeau et al. (2017) that epipelagic, mesopelagic and bathypelagic bioregions, and hence ecosystems, are different). This is obvious and already known, but marine ecosystem studies using remote sensing surface chlorophyll estimates already interpret their results as if they were representative of the whole epipelagic ecosystem. Phenologies and bioregions defined from surface chlorophyll are characteristic of the surface phytoplankton dynamics only, and cannot be said to represent the whole epipelagic ecosystem.

5
Nonetheless, recent remote sensing products may still help with marine ecosystem analysis, for instance the community structure products of Alvain et al. (2006) and Uitz et al. (2006). Another solution could be using remote sensing derived primary production products (see for example for the Mediterranean sea Bricaud et al. (2002); Bosc et al. (2004); Uitz et al. (2012). But one must be very prudent using these primary production products, as it has been shown that depending on the algorithm, the resulting primary production can have up to a factor 2 of difference (Carr et al., 2006;Saba et al., 2010). In particular, the Bloom-Intermittently cluster overlaps that derived from satellite estimates, but has differences in its spatial Finally, we show that phenology and bioregions based on Chl surf are very different from those based on integrated phytoplankton biomass. It appears that surface chlorophyll alone is insufficient to describe epipelagic ecosystem functioning, and hence cannot be considered as a good proxy of the whole phytoplankton biomass (at least in the Mediterranean). So phenology and bioregional studies that are based on remotely sensed surface chlorophyll estimates are only representative of surface plankton dynamics and not the whole epipelagic ecosystem.         contribution to Chl tot was too unbalanced compare to observations, to perform correct Chl tot analysis. This simulation shows us anyway that the improved intermediate water ventilation helps improving the nutricline, which in turn improved the DCM concentration and depth. It might also helped increasing the surface Chl minimum in summer, but that run also includes atmosphere deposition what is more likely to be the determinant factor that boosted the summer surface primary production (as shown in Richon et al. (2018c)). Both simulations (this study and Richon et al. (2018c)) present interesting differences, that 5 would benefit from more detailed comparison.