In Figure 2, the model simulated results are shown against the observations for the three mesocosm treatments (control without any addition of dust, single addition/SA and repetitive addition/RA). These include dissolved inorganic nutrients, Chlorophyll-a (Chl-a), net primary production, biomass of picophytoplankton and nanophytoplankton that represents more than 93% of the total phytoplankton, bacterial production and biomass, as well as the biomass of heterotrophic nanoflagellates (HNAN), microzooplankton and mesozooplankton.

As expected, the dust addition in SA and RA treatments provoked an increase in both dissolved inorganic nitrogen (DIN = nitrate + ammonium) and phosphorus (phosphate, PO₄) observed concentrations (Figures 2A,B), as compared to the control experiment. This difference was larger for DIN, as the amount of nitrogen added with the dust is much higher (+100% of the initial DIN concentration) compared to phosphorus (+30% of the initial PO₄ concentration). A decreasing trend can be seen for the observed DIN concentration of the control treatment and particularly for PO₄ (all treatments), which is related to the nutrient uptake by phytoplankton and bacteria. Interestingly, the measured DIN concentration in the addition treatments shows an increasing trend, suggesting a potential nitrogen excess due to a stronger phosphorus limitation. The model captured the observed decreasing trend of PO₄ and control DIN concentration, but with slightly lower concentrations by day-10 (Figures 2A,B). The higher increase in DIN (as compared to PO₄) in the addition treatments (SA/RA) was also simulated. However, the observed increasing trend of DIN in SA/RA was not simulated, even though an increasing DIN/DIP ratio from ~13 on day-1, indicative of nitrogen and phosphorus co-limitation, to ~35 on day-10, suggesting a stronger phosphorus limitation, was simulated (not shown) in agreement with the observations. One can also notice that the model PO₄ in the addition treatments presented a slightly stronger decrease after day-5, as compared to the control, which can be explained by the larger nutrient uptake from the increased phytoplankton biomass (Figures 2E,F), resulted from the nutrient enrichment by the dust additions. This stronger simulated decrease in PO₄ led to a slightly lower mean PO₄ in SA/RA over the 10-day period (Figure 3), in contrast with the observations that show an overall PO₄ increase in the addition treatments. We should note, however, that the observed PO₄ in the control treatment appears occasionally (day-5, -9) higher than SA/RA after day-4, which is in agreement with the model results.

Following the enrichment with dissolved inorganic nutrients, the simulated net primary production (netPP) appears enhanced in the dust addition treatments (Figure 2C), in agreement with the observations, even though their simulated decline starts a bit sooner (~day-3) compared to the observed (~day-6). The overall 10-day simulated increase in netPP in SA/RA (~+25%) is slightly smaller than the observed one (Figure 3). The SA netPP presents a slightly higher increase, as compared to RA, in the model simulation. This result seems reasonable, considering that in SA, added nutrients are available to phytoplankton from the start of the experiment, giving more time to stimulate the production of biomass, as compared to RA, where the same added amount is portioned on a 3-day period. Interestingly, the observations show the opposite, with RA showing a slightly higher overall increase of primary production, as compared to SA. The observed difference between SA and RA, however, might be considered relatively small (<10%, Figure 3), considering their standard deviation (from the three replicates), making it difficult to conclude safely that RA presents a higher production.

The simulated evolution of picophytoplankton biomass, showing a peak at ~day-3 and a secondary increase after day-8 that is mostly related to the decrease of its predator, HNAN (Figure 2I), is in good agreement with the observations, except a slight overestimation in the control treatment (Figure 2E). Both observations and model results show a higher picophytoplankton biomass in SA/RA treatments. The simulated nanophytoplankton biomass appears slightly overestimated (Figure 2F), showing a similar pattern to picophytoplankton, but with a slightly steeper decline after day-3, which may be attributed to the predation pressure exerted by increasing microzooplankton (Figure 2J). The same pattern is depicted from observations in the SA experiment. The observed evolution of nanophytoplankton biomass appears less clear in the control and RA experiments, with the latter showing the lowest biomass values. The observed Chl-a, dominated by picophytoplankton (Lagaria et al., 2017), is reasonably well-reproduced by the model, except an overestimation in the control treatment and a slight time-lag of its peak (Figure 2D). This deviation might be related to the observed time variability of phytoplankton Chl-a:C ratio, showing higher values between day-4 and -8, in the control treatment (Lagaria et al., 2017), which explains the measured Chl-a increase during this period. This observed variability cannot be captured by the model, which adopts a fixed Chl-a:C ratio.

The observed bacterial production (BP) presented a significant increase in the addition treatments, showing peaks between day-2 (SA) and day-3 (RA), a decline on day-6 and an increasing trend afterward (Figure 2G). The same growing trend at the end could also be seen in the control BP that was otherwise relatively constant until day-7. A similar pattern was followed by the observed bacterial biomass, except for an initial decrease in the control, until day-7 (Figure 2H). Therefore, bacterial production and biomass appeared closely coupled to phytoplankton biomass, particularly picophytoplankton that showed a very similar evolution. This was to a point expected, as both picophytoplankton and bacteria are prey for HNAN. As the latter decreased throughout the 10-day period (Figure 2I), its predation on picophytoplankton and bacteria was relaxed, allowing for their biomass increase at the end of the experiment. This top-down control, however, cannot explain the significant stimulation of BP in SA/RA that appeared to drive also the bacterial biomass increase. Another strong coupling mechanism is the production of dissolved organic carbon by phytoplankton that bacteria rely on for their growth. As later discussed in the model sensitivity section, on the short time scale (~days) of the experiment, one may assume that bacterial production is mostly affected by the labile organic matter in situ produced/excreted by plankton, rather than the semi-labile organic carbon that is decomposed on a longer time-scale. The increase in phytoplankton due to the nutrient enrichment with dust additions may thus explain the observed increase in BP. Under the assumption of bacteria utilizing mainly in situ produced dissolved organic carbon (see later discussion in section Bacterial dynamics), the model was able to reproduce the observed BP variability and particularly its increase in SA/RA experiments (Figure 2G). An overall BP increase of about +40% was simulated in SA/RA relative to the control, which was slightly smaller, as compared to the observed increase (+65%) in BP (Figure 3). The simulated bacterial biomass has a similar pattern with the observed, but shows a much weaker variability (Figure 2H). We should note however that the observed bacterial biomass initial increase (~1 mgC/m³/day) would probably require about twice the observed BP (~0.5 mgC/m³/day), even ignoring bacteria mortality and predation losses. This suggests that the conversion factor (20 mgC/cell, Lee and Fuhrman, 1987) used to calculate carbon biomass from bacteria cell abundance might be a bit overestimated. Indeed, the bacteria variability was much better reproduced by the model (not shown), when the observed bacteria biomass (and its initial value used in the model simulations) was decreased by some factor.

The model reasonably captured the evolution of other heterotrophs (HNAN, microzooplankton, mesozooplankton). A continuous decrease in HNAN with time was simulated in agreement with the observations (Figure 2I). The model HNAN was slightly higher in SA/RA, as compared to the control, while observations suggested a stronger HNAN overall increase in SA (Figure 3). The observed increase in both microzooplankton and mesozooplankton was well-reproduced by the model that also simulated a decrease of microzooplankton mainly after the end of the experiment, related to the predation by mesozooplankton (Figure 2J). Again, as in most cases, the model simulated an overall higher biomass in SA, followed by RA, for both microzooplankton and mesozooplankton (Figure 3).

In order to investigate the overall effect of the dust addition in the two different treatments, both in terms of observations (DATA-SA & DATA-RA) and model simulations (MODEL-SA & MODEL-RA), the fractional change of the Addition with respect to the Control treatment [(Addition − Control)/Control] was computed (Figure 3). Considering the standard deviation from the three replicates, overall there are no significant differences in the observations between the single (SA) and repetitive (RA) addition treatments, both showing a positive change, with the exception of nanophytoplankton, where the repetitive addition results in a negative fractional change in contrast to the single addition, which is positive and microzooplankton that shows a negative change in both SA/RA. In contrast with the observations, the model results show a positive change in all cases. The model produces a smaller positive change for autotrophs, bacteria and grazers, compared to the observations, with the exception of HNAN in RA, showing a slightly higher positive change than the observed.

The model skill in reproducing the observed variability (Figure 2) for different variables may be graphically summarized in a quantitative way using taylor (Taylor, 2001) and target (Jolliff et al., 2009) diagrams (Figure 4). In the taylor diagram, one can see the model skill for different variables regarding correlation, standard deviation (as an index of variability) and RMS error against the observations, while the target diagram also indicates the model bias and unbiased RMS error. From the two diagrams, one can identify some model variables (mesozooplankton, PO₄, Chl-a) presenting relatively good scores in all skill indexes (correlation r > 0.6, STD ~1, BIAS < 0.5). Nanophytoplankton presents a good correlation (r ~0.7) and the correct variability (STD ~0.8), but is slightly biased, as indicated in the target diagram. Picophytoplankton, bacterial production and biomass show good correlation and bias scores, but have relatively small variability, particularly the bacterial biomass. Net primary production is a bit out of phase with the observed one (Figure 2C), which results in its poor correlation and RMS error. DIN appears with a slightly negative correlation, as the model fails to capture the observed increasing trend at the end of the experiments (Figure 2B).

Models are excellent tools in exploring processes which are difficult or impossible to monitor and measure in the field. To investigate the impact of nutrient additions in terms of carbon flows within the food web and the possible differences between the two addition treatments, the corresponding carbon fluxes and their fractional change [(Addition − Control)/Control] were computed (Figure 5). In most cases, there was an initial increase in the carbon fluxes, followed by a gentler decrease. Slightly different from the dominant evolution, is the continuously decreasing flux from bacteria to heterotrophic nanoflagellates (B1Z6), related to the decreasing biomass of both Z6 and bacteria (Figure 2), while the flux from microzooplankton to the higher predator of the food web, mesozooplankton (Z5Z4) followed an increasing trend. Overall, in both treatments, there was an increase in the carbon flux from prey to predator but depending on the time scale of the process, this increase settled back to normal (in this case the control run) toward the end of the period. Smaller heterotrophs reacted faster to the dust addition, taking advantage of their relatively higher growth rates and the increasing abundance of their prey (P2Z6, P3Z6), followed by bigger animals such as microzooplankton (P2Z5, P3Z5, Z6Z5), while the 10 days experimental period did not seem to be enough for the biggest ones (mesozooplankton) to converge to the control. In all cases, the single addition treatment resulted in a small but visible initially higher carbon flux, converging after day 4 or 5 with the repeated addition treatment, which seemed to have a couple of days delay. The same dominant trend (initial increase, more gradual decrease) was visible in the fractional change of fluxes from phytoplankton to heterotrophs (P2Z5, P2Z6, P3Z5, P3Z6). This was maximized before day-5 at around +50%, with a small time-lag between SA and RA, before these started settling back to the control. This was not the case for bacteria (B1Z6) and mesozooplankton (Z5Z4), where SA/RA did not appear to converge back to the control.

From the mean food web fluxes (Figure 6), model results showed that most of the carbon flowed from small phytoplankton and in particular from picophytoplankton (P3) to the smaller heterotrophs (Z6 & Z5). Although mean fluxes in the addition treatments were higher compared to the control run, there was no change in the pattern. This supports the argument that the dust addition has not modified the food web structure i.e., causing a shift from a microbial food web to a more classical food web, dominated by large phytoplankton-zooplankton, but just increased the trophic status (less oligotrophic), with more carbon circulating in the entire food web. Thus, the increase of different fluxes (SA/RA-Control, Figure 6) appears proportional to each flux magnitude. The two major pathways of carbon, further stimulated by the dust additions, were those from the photosynthetic P3 and bacteria to the smallest heterotrophic nanoflagellates (P3Z6, BZ6), and the second and greater channel from the small phytoplankton to microzooplankton (P2Z5 & P3Z5). Besides these two main channels, smaller pathways were the predation processes within the zooplankton groups (Z6Z5, Z6Z4, Z5Z4).

To investigate the fate of the added nitrogen and phosphorus, the evolution of the (Addition − Control) difference in different nitrogen and phosphorus pools (dissolved inorganic, particulate and dissolved organic, phytoplankton, bacteria, and zooplankton) was computed (Figure 7). This nicely illustrates the transfer of added nitrogen/phosphorus as this passed from one pool to the other. Initially the difference was zero in all pools except for the dissolved inorganic (PO₄, DIN). These were quickly reduced, taken up by phytoplankton, which was the first to increase, showing a peak on day-3 and decreasing afterward due to nutrient limitation and predation by zooplankton. Nitrogen and phosphorus were then directed to zooplankton that followed phytoplankton with a time-lag (~4 days). Bacteria, stimulated by the organic matter produced by phytoplankton/zooplankton, also took a major part of nitrogen and phosphorus pools, consuming dissolved organic nitrogen and phosphorus that appear to decrease. Particulate organic pools were initially reduced, being smaller in the addition treatments, as phytoplankton mortality is a function of nutrient limitation. After day-3 they gradually built up from mortality losses of all plankton groups. At the end of the 20-day period, some of the phosphorus pools (PO₄, phytoplankton-P DOP, bacteria) appeared to be settling back to the control (Addition − Control = 0), with the added phosphorus remaining mostly in the form of POP and zooplankton-P. For nitrogen, the largest pool remained in the form of DIN, given its excess over DIP (DIN/DIP ~35), with the increase of phytoplankton-N pool being related to luxury uptake. The other pools (bacteria, PON, zooplankton-N) followed a similar pattern as in phosphorus.

In ERSEM, the nutrient uptake is constrained by a maximum uptake rate (V₀) that depends on the dissolved inorganic nutrients external concentration and phytoplankton affinity. In the case of phosphorus uptake (PO₄) this is:

with qurP being the specific affinity for phosphorus uptake (Table A1 in Supplementary Material), while this is regulated by the amount of nutrient that is necessary for the cell to address its needs for growth and intracellular storage:

where runP is the phytoplankton net production (photosynthesis-excretion-respiration), qpP_max is the maximum phosphorus internal quota (qpP_max = 2 × Redfield ratio, Table A1 in Supplementary Material) and PhytoC, PhytoP are the carbon and phosphorus phytoplankton pools respectively. The actual nutrient rate is then taken as:

With this formulation and the adopted specific affinity parameter values [qurP = 0.0025 (mgC/m³)⁻¹ day⁻¹], the simulated phytoplankton nutrient uptake was quite low, as V₀ (Equation 1) was constrained by the very low initial nutrient concentrations of the mesocosms (~0.005 mmol/m³ PO₄, ~0.1 mmol/m³ DIN). This reduced nutrient uptake resulted in a close to zero primary production (Figure 8, Run2), as it quickly leads to a sub-optimal phytoplankton internal stoichiometry that triggered an increased carbon excretion and cell lysis, as adopted in ERSEM formulation (Baretta-Bekker et al., 1997):

where sumP is the carbon uptake, seoP_max is the maximum fraction excreted, sdoP is the lysis rate and Qlim is the Droop nutrient limitation function of internal cell quotas (qpP, qnP), which in the case of phosphorus is:

with qpP_max/qpP_min being the maximum/minimum phosphorus internal quotas.

In order for simulated net primary production and phytoplankton biomass to approach the observed values in the mesocosm experiment, much higher affinity parameters [qurP × 20 ~0.05 (mgC/m³)⁻¹day⁻¹] were necessary (Figure 8, Run3). An even better agreement of simulated results, particularly in terms of dissolved inorganic nutrients that were slightly reduced toward observations, was found when nutrient uptake rates were assumed to follow the classic Michaelis-Menten kinetics as in Geider et al. (1998):

where sumP is the carbon specific maximum growth rate, qpP_red is the Redfield phosphorus quota, K_p is the half-saturation (fitted parameter) for phosphorus uptake and Qmax = (qpP_max − qpP)/(qpP_max − qpP_min), with qpP_max/qpP_min being the maximum/minimum phosphorus internal quotas (see Table A1 in Supplementary Material). Qmax approaches zero when the phytoplankton internal phosphorus quota Q is maximum (= qpP_max), resulting in the decline of nutrient uptake. When Equation (7) was used for nutrient uptake instead of Equation (3), the simulated net primary production better reproduced the observed evolution, showing a slightly more extended peak (Figure 8, Run1), while simulated inorganic nutrients got closer to the slightly lower observed values, particularly for phosphate. This is related to the use of the Michaelis-Menten sigmoid function in Equation (7) that results in a more gradual decrease of nutrient uptake and the steepest decrease of nutrients, as compared to those simulated when the linear function Equation (1) is used.

The phytoplankton affinity, defined as the initial slope of nutrient uptake at very low nutrient concentrations (PO₄ + K_p~K_p) can be calculated from Equation (7) as qurP~V_max/K_p, which gives 0.11 and 0.13 (mgC/m³)⁻¹day⁻¹ for nanophytoplankton and picophytoplankton, respectively. These are slightly higher than the fitted affinity parameters [0.05–0.0625 (mgC/m³)⁻¹day⁻¹] using the original ERSEM formulation and much higher than the previously adopted affinity parameters [0.0025 (mgC/m³)⁻¹day⁻¹].

The prey uptake by Zi zooplankton heterotrophic groups (Z4 = mesozooplankton, Z5 = microzooplankton, Z6 = heterotrophic nanoflagellates, see Figure 1) in ERSEM is described by a Holling-II type function:

where K_Zi is a half-saturation constant (where the uptake rate is half its maximum value) and FtotZi is the total available amount of food from different sources:

with su_ij being the preference of Zi group on different preys F_j and min food_Zi another half-saturation constant used to prevent from exhausting a prey food source when this is scarce, as compared to some minimum food requirement (see Tables A2, A3). The total amount of food is thus calculated from Equation (8), based on the preference and the relative availability of different preys. In order to correctly simulate the observed biomass of zooplankton groups, the half-saturation constant (K_Zi) parameters were decreased from their initial values (K_Z4 = 14, K_Z5 = 24, K_Z6 = 49 mgC/m³, see Table 1). As shown in Figure 9U, mesozooplankton was initially significantly underestimated (Run4, see Table 1). Adopting a lower half-saturation (K_Z4 = 3 mgC/m³), the simulated mesozooplankton biomass in Run5 approached the observations (Figure 9V) but the biomass of its prey, microzooplankton, was decreased (Figure 9R). Decreasing also the microzooplankton half-saturation constant (K_Z5 = 10 mgC/m³) in Run6, the microzooplankton underestimation was removed (Figure 9S), mesozooplankton was further increased (Figure 9W), but an underestimation was now found for heterotrophic nanoflagellates (Figure 9O). This was finally removed (Figure 9P) in Run1 when the respective half-saturation constant was also decreased (K_Z6 = 16 mgC/m³), showing also a slightly better agreement for microzooplankton and mesozooplankton (Figures 9T,X). This series of sensitivity experiments illustrates the prey-predator trophic relations. One may notice, for example, the increase in bacteria in Run6 (Figures 9I–L), following the decrease of their predator, HNAN (Figures 9M–P), or the increase in nanophytoplankton (Figures 9A–D) and picophytoplankton (Figures 9E–H) in Run4/Run5 due the relatively low microzooplankton biomass (Figures 9Q–T).

The observed variability in the mesocosm experiment was characterized by a significant enhancement of bacterial production and biomass in the two addition treatments (SA, RA), as compared to the control mesocosms. This increased bacterial productivity, triggered by dust additions, appeared to be closely coupled to the phytoplankton productivity increase. The ERSEM bacteria sub-model in Petihakis et al. (2002) has been revised by Petihakis et al. (2015), following Anderson and Williams (1998) and Petihakis et al. (2009), allowing for a better representation of the DOM pool, which was particularly important in the presence of significant lateral inputs (e.g., rivers, Black Sea Water) of DOM. The uptake of DOM by bacteria is described by:

where sumB is the maximum bacterial growth rate, f(T) is the growth temperature dependence, f(O2) the oxygen limitation, (N_Lim, P_Lim) is the nutrient limitation on nitrogen/phosphorus and C_Lim the limitation on available DOC. The bacteria nutrient limitation (N_Lim, P_Lim) is assumed to depend on the intracellular nitrogen (N/C) and phosphorus (P/C) bacteria quotas (qnB, qpB), compared to the maximum internal quotas (qnB_max, qpB_max):

Alternatively, the bacteria nutrient limitation may be assumed to be a function of external nutrient concentrations, as adopted by Blackford et al. (2004):

Where DON/DOP is dissolved organic phosphorus/nitrogen and K_N/PBac is a half-saturation constant.

The carbon limitation (Clim) is described by:

where DOC represents labile and semi-labile DOC and K_DOC is a half-saturation constant (Table A2).

When the DOC pool was initialized from the mesocosm measurements (~400 mgC/m³ Figure 10E, Run7), the simulated bacterial production/biomass overall mean was consistent with the measured values, but the model was not able to reproduce the observed variability and particularly the strong differentiation of the addition treatments in relation to the control (Figures 10A,I, Run7). DOC is produced by mortality, excretion and lysis of primary and secondary producers and is mostly labile that is readily available for bacteria. Semi-labile DOC takes longer to decompose by bacteria and can accumulate with time. One can assume that the initial measured DOC pool was mostly semi-labile, while bacterial production on the short (~days) time scale of the experiment was primarily affected by labile organic matter, as the one produced/excreted by plankton biomass. To test this hypothesis, the initial DOC pool was significantly decreased (~3 mgC/m³ Figures 10F–H), considering only its labile component, as also the half-saturation for DOC uptake (K_DOC = 4 mgC/m³). Indeed, in this case, bacterial production presented an increase in the addition treatments (Figure 10D, Run9), being closely coupled to the enhanced plankton biomass that released additional DOC. The simulated bacterial production presented a slightly weaker variability, as compared to the observed that showed a stronger and slightly earlier peak in the addition treatments, as well as a higher increasing trend at the end of the 10-day period. Moreover, the model failed to reproduce the observed early peak of bacterial biomass in the addition treatments (Figure 10L, Run9). The simulated bacterial production and biomass were further improved (Figures 10C,K, Run8) when bacterial nutrient limitation was assumed to be a function of external nutrient concentrations (Equation 12), rather than the bacteria internal quotas (Equation 11). In this case, bacterial production appeared to be directly stimulated by the addition of inorganic nutrients, showing an earlier peak on day-3, in agreement with observations, as well as a stronger differentiation between the addition and control treatments. The simulated bacterial biomass, despite the weaker variability, now showed a similar pattern with observations, initially increasing in SA/RA and decreasing in the control (Figure 10K). The simulated bacterial production was still slightly underestimated in the SA/RA treatments. Given its strong sensitivity to the bacterial assimilation efficiency parameter (see next section Model parameter sensitivity), a final model improvement was explored, adopting a variable bacterial assimilation efficiency, depending on nutrient limitation. Specifically, a function of nutrient limitation (NPlim, Equation 12) was introduced in the computation of bacterial respiration as:

where U_B is the bacterial growth rate (Equation 10), puB/puBo is the bacteria assimilation efficiency parameter at sufficient/low oxygen (Table A2), eO2 is the relative oxygen saturation, R_Bb is the temperature dependent basal respiration (R_Bb = srsB^*f(T), Table A2) and

The introduced function of nutrient limitation varied between 0.91 and 1.14 in the control and between 0.95 and 1.24 in SA/RA treatments and effectively resulted in a slightly stronger variability of Bacterial Growth Efficiency (BGE~(U_B − R_B)/U_B), depending on nutrient limitation. This model modification lead to a higher bacterial production and biomass in the addition treatments (Figures 10B,J) that were characterized by a relaxed nutrient limitation, improving the model fit with observations.

A minimum set of parameters/formulations of the existing model were carefully revised, as discussed in the previous section in order to obtain an optimum fit with the observations. The rest of model parameters and formulations were kept fixed to their old values, as their modification did not show any significant improvement. To test the model sensitivity to different parameters, a series of sensitivity experiments were performed, adopting a 10% increase in the value of chosen parameters (see Tables A1, A2 in Supplementary Material for their definition and reference value). The output from these sensitivity simulations with modified parameter values was then compared to the reference simulation, computing the fractional change [(sensitivity − reference)/reference] for different variables. As shown in Figure 11, dissolved inorganic nutrients were mostly affected by phytoplankton maximum growth rate (sumP) and half-saturation coefficients for nutrient uptake (K_P, K_N). An increase in sumP results in the increase of phytoplankton growth rate and thus, to the decrease of nutrients, while an increase in K_P/K_N results in a reduction of nutrient uptake rate (see Equation 7). Phytoplankton biomass (p2c, p3c) is mostly sensitive to sumP, as well as to the lysis rate (sdoP) and particularly to the grazing half-saturation constants (K_Z5, K_Z6) that affect their predation by zooplankton (z5c, z6c, see Equation 8). An increase in K_Z5 results in the increase in nanophytoplankton (p2c) and heterotrophic nanoflagellates (z6c), the main preys of microzooplankton (z5c), while increasing K_Z6 mainly affects picophytoplankton (p3c). Finally, increasing K_Z4 mainly affects mesozooplankton and to a lesser degree its preys z5c and z6c. Bacterial production is among the most sensitive model variables, particularly affected by bacterial assimilation efficiency (puB), growth rate (sumB) and basal respiration (srsB), as well as to the half-saturation for DOC uptake (K_DOC). It is also sensitive to the first order breakdown rate from POC to DOC (fr6r1, Petihakis et al., 2015) and the maximum DOC phytoplankton excretion rate under nutrient limitation (seoP_max), both related to the supply of DOC that is necessary for bacteria to grow. Bacterial biomass depends on the same parameters but appears much less sensitive.