Spatio-temporal dynamics of a planktonic system and chlorophyll distribution in a 2D spatial domain: matching model and data

Field data on chlorophyll distribution are investigated in a two-dimensional spatial domain of the Mediterranean Sea by using for phytoplankton abundances an advection-diffusion-reaction model, which includes real values for physical and biological variables. The study exploits indeed hydrological and nutrients data acquired in situ, and includes intraspecific competition for limiting factors, i.e. light intensity and phosphate concentration. As a result, the model allows to analyze how both the velocity field of marine currents and the two components of turbulent diffusivity affect the spatial distributions of phytoplankton abundances in the Modified Atlantic Water, the upper layer of the water column of the Mediterranean Sea. Specifically, the spatio-temporal dynamics of four phytoplankton populations, responsible for about 80% of the total chlorophyll a, are reproduced. Results for phytoplankton abundances obtained by the model are converted in chlorophyll a concentrations and compared with field data collected in twelve marine sites along the Cape Passero (Sicily)- Misurata (Libya) transect. Statistical checks indicate a good agreement between theoretical and experimental distributions of chlorophyll concentration. The study can be extended to predict the spatio-temporal behaviour of the primary production, and to prevent the consequent decline of some fish species in the Mediterranean Sea.

shape and width of the DCM.
The chemical analysis performed on the bottle samples confirms that the Mediterranean Sea is a phosphorus limited basin due to the high N:P ratio measured along the water column in all sites investigated [1,3]. Specifically, the phosphate concentration takes on very low values due to the absorption of Saharan dust, while the higher concentration of nitrate, nitrite and ammonium are strictly connected to: (i) the excess of nitrogen fixation in the Mediterranean basin, (ii) the nutrient flow coming from the coasts.

Phytoplankton properties
In this work, the contribution of each picophytoplankton group to the total amount of chlorophyll is based on the experimental estimation of cellular chlorophyll a content obtained by performing the high-performance liquid chromatography (HPLC) analysis on the seawater samples collected between Cape Passero and Malta during a previous oceanographic survey [4,5].
According to the study carried out in the Mediterranean Sea by other authors [4][5][6][7], the phytoplankton community can be divided into three main size fractions: pico-(< 3µm), nano-(3−20µm) and micro-phytoplankton (> 20µm). Specifically, in the marine ecosystem studied (Strait of Sicily), the analysis is focused on the picophytoplankton fraction, which accounts, on average, about for 80% of the total chl a and Dvchl a [4,5,[8][9][10], and consists of two main domains: picoprokaryotes and picoeukaryotes. The former is composed of two genera of cyanobacteria, i.e. Synechococcus and Prochlorococcus. The latter is dominated by Prymnesiophytes, Pelagophytes and green algae. On the other hand, the nano-and micro-phytoplankton fraction accounts for about 20% of the total chl a and Dvchl a on average, and is mainly represented by the eukaryotes domain, i.e. Prymnesiophytes, Pelagophytes, Dinophytes and diatoms. This fraction is poorly present in DCM, and is almost uniformly distributed along the water column.
In this study, we analyze the behaviour of four picophytoplankton populations, i.e. Synechococcus, Prochlorococcus (HL-ecotype), Prochlorococcus (HL-ecotype) and the whole picoeukaryotes domain, which are located at different depths along the water column. Field observations indeed indicate a prevalence of Synechococcus close to the water surface, while Prochlorococcus and picoeukaryotes dominate the intermediate deeper layers [4,5].
Synechococcus is mostly present close to the coasts, where its cell concentration reaches usually the maximum value. Moreover, the analyses carried out along the water column showed the presence of nine different phylogenetic groups (clades) of Synechococcus in the Mediterranean Sea, most of which are observed in the Strait of Sicily [9].
The Prochlorococcus concentration is characterized by a bimodal distribution in the intermediate layers of the water column [5,9], indicating the coexistence of two ecotypes of this genus: high light-adapted (HL-) ecotype and low light-adapted (LL-) ecotype. The former ecotype is localized in the upper part of the MAW between the surface water and 90 m of depth. The latter is mostly present at depths greater than 50 m [5,8,9]. Moreover, recent studies showed that the Prochlorococcus HL-ecotype prevails close to the Sicilian coast [4,5,9], while Prochlorococcus LL-ecotype dominates the marine ecosystem between the Sicilian -Maltese shelf and the Libyan coast [8,9]. However, the biological features of both Prochlorococcus ecotypes can change along the Cape Passero -Misurata transect. In our study, the prevalence of different clades, ecotypes and/or groups inside the four picophytoplankton populations, whose dynamics is modeled in this work, is taken into account for each hydrological station of the Cape Passero -Misurata transect. In particular, we use different settings for some biological parameters, such as the half-saturation constants for light intensity and nutrient concentration, in order to consider the phylogenetic modifications among the marine stations investigated.
The analysis of seawater samples shows that Synechococcus contributes on average to more than 20% of the total chlorophyll concentration in the Strait of Sicily [4,5]. In particular, the average concentration of Synechococcus is 5.8 × 10 3 cell ml −1 , while its concentration peak (1.05 × 10 4 cell ml −1 ) is localized at 50 m of depth. The chl a cellular content of Synechococcus has never been estimated in the Strait of Sicily. Therefore, assuming the oligotrophic conditions, we chose to use the content measured by Vaulet and Courties in the English Channel waters, whose value was fixed equal to 1.18 fg chl a cell −1 [11].
In the Strait of Sicily, Prochlorococcus and picoeukaryotes contribute equally to the primary production in terms of chl a and Dvchl a concentrations in intermediate and deeper layers of the MAW. However, picophytoplankton is numerically dominated by Prochlorococcus with an average concentration of 5.2 × 10 4 cell ml −1 . Specifically, this group is mainly localized in DCM, where can reach the mean value of 12.5 × 10 4 cell ml −1 . The marker of Prochlorococcus is divinil chlorophyll a, whose molecular structure is very similar to that of chlorophyll a. The Dvchl a cellular content of total Prochlorococcus ranges between 0.25 and 2.20 fg Dvchl a cell −1 along the water column, with a mean value exponentially increasing with the depth [5].
The analysis performed on the seawaters shows that the average picoeukaryotes concentration in the DCM is 0.6 × 10 3 cell ml −1 [4,5,7], while the mean value of chl a cell −1 ranges between 10 and 660 fg chl a cell −1 along the water column, with a significant exponential increase with the depth [5]. Therefore, the concentration of chl a per picoeukaryotes cell changes within the MAW, assuming significantly higher values in the DCM respect to the upper layer [12,13].
On the basis of these findings, we convert the picophytoplankton abundances into chlorophyll concentration by using the two conversion curves obtained by Brunet et al. [5] for Prochlorococcus and picoeukaryotes (see Fig. 2

Setting of parameters
In this work, we reproduce the two-dimensional distribution of the total chlorophyll concentration experimentally observed in the Strait of Sicily, by setting the values of the environmental and biological parameters so that the monostability condition in the marine ecosystem is obtained. In particular, initially the environmental parameters are estimated by the experimental data collected along the Cape Passero -Misurata transect. Afterwards, on the basis of the environmental conditions observed along the water column, the biological parameters are fixed in such a way to guarantee the coexistence of all four planktonic populations [14][15][16] in the intermediate and deeper layers. By this way, we mimic the presence of a deep chlorophyll maximum (DCM) in all hydrological stations, although the peak of abundance for each group is localized at a different depth [14][15][16][17][18].
In accordance with the field observations, the hydrodynamical variables are fixed in the model at constant values during the whole sampling period investigated (15-30 July 2008), even if they can change during the year. Since the only field data available in this study were those acquired during the MedSudMed-08 oceanographic survey (summer season 2008), we did not have any data to reproduce the effects of seasonal changes in hydrological variables during the whole period investigated by the model (approximately two years and three months). However, in previous works [16,19] it has been shown that, in Mediterranean Sea, the seasonal changes of environmental parameters do not modify significantly the steady spatial distributions of phytoplankton abundance obtained in summer season, since the hydrological variables remain almost constant between late spring and late autumn. For these reasons, the stationarity assumption can be considered valid, even though the steady solu- In this work, the vertical turbulent diffusivity is reproduced according to the method by Pacanowski and Philander [20,21], by using the experimental profiles of geostrophic velocity components and density collected along the Cape Passero -Misurata transect. In particular, to obtain the vertical turbulent diffusivity as a function of the depth, D v (z), in each sampling site we use the following expression based on the empirical studies [20,22,23]: where ν 0 = 36m 2 /h, α = 5.0, and n = 1.0 are adjustable parameters chosen for lower intensities of wind stress, ν b = 0.36m 2 /h is the background dissipation parameter, Ri(z) is the gradient Richardson number, depending on the depth, given by Here, the buoyancy frequency N (z) is calculated directly by the vertical profile of water density as follows where g is the gravity acceleration, ρ w is the density of the sea water, and ∂ρ(z) ∂z is the vertical density gradient. On the other hand, the shear of the mean currents, shear(z), depends on the vertical profiles of the horizontal velocity components according to the following expression [20,21] shear 2 (z) = ∂u(z) ∂z where u(z) and v(z) are the zonal and the meridional geostrophic current components mea-  [4,5,11,12,18,[26][27][28][29][30][31][32][33][34][35].
In particular, the maximum specific growth rates are fixed in accordance with experimental results given in Refs. [28][29][30], while the specific loss rates are estimated on the basis of experimental findings given in Refs. [28,[30][31][32]. On the other hand, the swimming velocity and nutrient recycling coefficients are chosen in agreement with the theoretical results obtained by other authors [27,31]. Specifically, the magnitudes of swimming velocities of the four planktonic populations are set to the values obtained by Raven [27], while nutrient recycling coefficients are calculated by considering the assimilation efficiencies estimated by Thingstad [31].
In previous works [14][15][16], as well as in our preliminary analysis, the half-saturation constants of the picophytoplankton groups are fixed to the average values observed experimentally (see parameters of the full model in Table I Conversely, on the basis of the field observations, the core of this study exploits halfsaturation constants set to different values in the twelve hydrological stations investigated (see parameters of the reduced model in Table I). Indeed we recall that in the Strait of Sicily the picoeukaryotes domain includes several groups [4,5], while Synechococcus, Prochlorococcus HL and Prochlorococcus LL are characterized by the presence of different clades [8,9].
For each phytoplankton population this biodiversity is a marker of the skills of adaptation at the different environmental conditions observed within the 2D domain of the marine ecosystem. As a consequence, depending on the marine site analyzed, different clades, ecotypes and/or groups prevail inside each population, and the half-saturation constants have to change accordingly [36]. Thus, the half-saturation constants, K R i and K I i , for the four populations are set so that the production layers and the peaks of phytoplankton abundance are placed at depths compatible with those observed by Brunet et al. [4,5] in the Strait of Sicily. In particular, according to Ryabov et al., for the i-th population and for each hydrological station, we consider the following expressions where R * i and I * i are the critical values, defined as the values of the resource availability at the boundaries of the production layer, r i is the maximum specific growth rate, and m i is the specific loss rate [15,37]. As a preliminary step, we calculate the values of half-saturation constants by using in Eqs. (5), (6) Table I). Specifically, the values of these parameters are estimated for Synechococcus and picoeukaryotes according to the experimental findings obtained in previous works [34,35], while no data are available for any Prochlorococcus ecotype (neither HL, nor LL). Therefore, in order to get phytoplankton abundances in accordance with the experimental results, we fix the nutrient content of the Prochlorococcus (both ecotypes) in such a way to respect the ratios of the average concentrations of the planktonic groups experimentally observed in the Strait of Sicily [4,5].
Finally, the chl a-normalized average absorption coefficients are calculated on the basis of the light absorption spectra obtained for phytoplankton cultures by Hickman et al. [16,18].
In general, the values estimated are in agreement with the absorption coefficients measured by Brunet et al. in the coastal zones of the Mediterranean Sea [12,13].

Statistical analysis
Since some parameters (seven half-saturation constants) are freely estimated in the model, it is worth to investigate the real goodness of the fit obtained by performing a comparison, based on the χ 2 test, between the reduced model and the full model [39]. Specifically, in order to establish which of them reproduces better the field data, we applied for both models the Akaike Information Criterion (AIC) [38,39] defined as follows where χ 2 df denotes a chi-squared with a number of degrees of freedom, df , equal to the number of independent constraints of the model, H. The minimum value of AIC indicates the model with the best fit.
In Table II, we show the AIC statistics calculated for each hydrological station and the whole transect. Here, the results indicate that the best AIC in the most sites (eleven hydrological stations over twelve) is obtained by the full model, while the χ 2 test show that the best fit is usually observed by using the reduced model. However, for the whole transect we observe that the best AIC statistics is obtained by using the reduced model, in accordance with the result of the χ 2 test.
Therefore, on this basis the reduced model seems to be the best tool to reproduce the overall two-dimensional distribution of chlorophyll concentration. On the other hand, the full model is able to well reproduce the experimental chlorophyll profiles in the most sites of the transect.
This conclusion can be further checked by using an other statistical tool, which defines a magnitude difference between the chi-square for the reduced model and the chi-square for the full model. In particular, we estimate this magnitude difference by using the Cohen's effect size measure [39] w = ( χ 2 df /(N · df )) 0.5 , where N is the sample size, equal for both models, χ 2 df is the difference in chi-square and df the difference in degrees of freedom between the two models.
According to the Cohen's suggested standard for a small effect (w ≤ 0.1) [39], the value obtained for the index w indicates that using the reduced model affects weakly (see Table II of Supplementary Tables) the result of the χ 2 test in the whole transect (w = 0.08).
Specifically, according to the values of w, the decrease of χ 2 is negligible in the sites from M 4 to M 12, while it is more significative, although still small, in the three hydrological stations localized close to the Sicilian coast. On the whole, the reduced model therefore allows to improve the results of the χ 2 test compared to the full model, even if this improvement is not actually very significative.