This section is devoted to show results and to highlight eventual effects of the interplay between the nonlinearity characterizing the system dynamics and the presence of noisy fluctuations for the irradiance variable.

### Analysis of experimental data

The need of taking into account noisy fluctuations of such an environmental variable is well demonstrated in Fig. 1. In the first panel (a) the experimental time behaviour of the irradiance is shown. This noisy curve is based on the experimental data (purple points) of the Boussole buoy located in the Gulf of Lion, collected over a period of nine years, precisely from 2004 to 2013. The time series of the experimental data presents quite a few gaps in time due to the malfunction of the buoy. This aspect has been remedied by merging the experimental data with those of the OASIM model validated for the Boussole site61 (yellow points). The latter is a multispectral atmospheric radiative transfer model that is in turn forced by experimental-model data based on ECMWF ERAINTERIM reanalyses which provide, for example, cloud cover data. The radiative model is partly stochastic since it considers the effects stemming from the presence of clouds, averaged along a single day (this explains why the yellow points are slightly less scattered). We see that the OASIM model accurately reproduces the profile which emerges from the experimental data. Further, we stress that the experimental data are only used in this initial analysis. In the biogeochemical simulations the irradiance signal is fully reconstructed starting from a realistic seasonal cycle combined with a range of different random fluctuations, and the information from OASIM is not used. In the second panel (b) the daily (black points) as well as the three-month (red points) running mean of the experimental series are plotted. Figure 1c shows the irradiance noisy fluctuations (INF) which have been obtained by subtracting the three-month running mean curve (3MRM, red curve in Fig. 1b) from the daily running mean one (DRM, black curve in Fig. 1b) and normalizing with respect to the mean of the 3MRM (\(\overline{3MRM}\)), namely \(INF = (DRM – 3MRM) / \overline{3MRM}\). We see that a seasonal overall trend with higher oscillations during the winter time can be seen, implying that the characteristics of the noise may change over the year. Moreover, a slight imbalance between positive and negative values of the noisy fluctuations (that is, different values of the maximum fluctuation intensity) is present. The physical reason for the occurrence of such an aspect can be ascribed to the fact that the maximum value of solar irradiance corresponds to that measured during a sunny day. Conversely, the minimum level tends to zero corresponding to a dense darkness. While the former is close to the mean value of the solar irradiance (most of all in summer), the latter is much further away and then a natural asymmetry arises in the random fluctuations. However, it should be noted that, apart from the intense spikes, the asymmetry is not so pronounced, as proved by the mean value (red line in Fig. 1c) which is practically zero, namely \(0.4\%\) of the \(\overline{3MRM}\). Therefore, basing on this last observation, to model the noise affecting the irradiance dynamics, as a first approximation we consider a symmetric Gaussian autocorrelated noise as described in the next subsection.

On the basis of such experimental results, we postulate the hypothesis that random fluctuations of light cannot be neglected, most of all in the study of ecological systems where light profoundly determines the system dynamics, governing fundamental processes at the basis of of the food web.

Figure 1

(a) Experimental data (purple points) of the stochastic solar irradiance collected by the Boussole buoy in a time-window of 9 years (2004-2013); the yellow points are the data generated by the OASIM model used to fill the gaps present in the experimental time-series due to malfunctioning of the buoy. (b) Daily (black points) and three-month (red points) running mean of the light curve in panel (a). (c) Irradiance noisy fluctuations (INF), obtained by subtracting the three-month running mean curve (3MRM) from the daily running mean one (DRM) and normalizing with respect to the mean value of 3MRM (\(\overline{3MRM}\)), namely \(INF = (DRM – 3MRM) / \overline{3MRM}\); the red line represents the mean value of such fluctuations. Data already presented and validated in61.

### Solar irradiance

The solar irradiance forcing is derived considering a deterministic seasonal oscillation combined with an Ornstein-Uhlenbeck process. The coefficient of variation (CV) of simulated light forcing, Fig. 2, \(CV=\sigma / \mu\) (\(\mu\) and \(\sigma\) being mean value and standard deviation calculated over both time and numerical realizations), is shown for 231 \(D-\tau\) pairs. D and \(\tau\) represent the intensity of a Gaussian noise source and the auto-correlation time of the fluctuations, respectively (see Eqs. (2) and (3)).

Each pixel represents the mean value on time of CV calculated with respect to 1000 different stochastic realizations.

Figure 2

Coefficient of variation (\(CV=\sigma / \mu\)) of irradiance resulting from numerical integration of model equations for 231 \(D-\tau\) different scenarios.

It is easy to see the agreement between the results obtained from the numerical integration and the theoretical ones derivable from Eq. (5) by putting \(\text {var}\{F_L(0)\}=0\) and \(t \gg 1\), getting \(\sigma ^2_L=D / 2\tau\). In Fig. 2, indeed, the maximum values of \(\sigma\) lie in the upper left part of the plot corresponding to small (high) values of \(\tau\) (D). As it is clear the values of D have been chosen in order to obtain a relative standard deviation ranging from \(5\%\mu\) to \(60\%\mu\). We underline that, in this case, it is possible to interchangeably consider \(\sigma\) and CV since the dependence of CV on D and \(\tau\) does not differ from that of \(\sigma\) (meaning that the dependence of \(\sigma\) is not altered by dividing by \(\mu\)) (results not shown).

### Effects on population dynamics

In this section the noise-induced effects on the population dynamics are examined. The nine planktonic populations present a different qualitative behaviour of the CV, compared to that of the irradiance. In this case, the CV is characterized by a strong non-monotonic dependence on the parameter \(\tau\). This aspect can be appreciated in Fig. 3 where different curves of CV versus the time correlation parameter are shown for different fixed values of D.

Figure 3

Coefficient of variation (\(CV=\sigma / \mu\)) of the nine planktonic populations resulting from numerical integration of model equations plotted versus the considered values of \(\tau\); the different curves are related to different values of the noise intensity D.

The existence of a maximum value for CV can be appreciated for each species. Although the qualitative behaviour is the same for all strains, particular attention has to be payed on diatoms and nanoflagellates. All the other species, indeed, present a percent variation of standard deviation between \(2\%\) and \(15\%\). In the case of nanoflagellates, instead, the D-dependent range is \(20-90\%\), while diatoms reach values over the \(100\%\) for the highest values of D. Therefore, these two species, in particular, and the whole system, in general, are extremely sensitive to the auto-correlation time which characterizes the noise.

We note that the different curves related to the different selected values of D approach the horizontal axis, tending asymptotically to vanish as \(\tau\) increases. Such a behaviour can be explained by the fact that high values of \(\tau\) give rise to a more correlated dynamics, so that \(\tau \rightarrow \infty\) implies fully correlated time-behaviours corresponding to the deterministic case. In this instance, then, all the different realizations give the same results, making the standard deviation vanish. The same happens, independently of the value of \(\tau\), for low values of noise intensity for which the corresponding curves approach the same almost vanishing value (see orange, gray and yellow lines). Differently from the previous case, when \(\tau \rightarrow 0\) the noise tends to a delta-correlated noise, that is a white noise; for \(\tau \ne 0\), instead, the noise spectrum is not flat, being characterized by a Cauchy-Lorentz distribution. The strong nonmonotonicity of CV with respect to \(\tau\), emerging when there are relatively high values of CV, implies a greater variability of the system biomass. Lower values of CV indicate that the system dynamics is less influenced by the presence of noise where very little or no differences with respect to the deterministic case are present. Conversely, high values of CV clearly demonstrate the remarkable signature of the presence of an impacting noise source. It is interesting to note that the noise influence on the ecosystem strongly depends on both \(\tau\) and D, that is, just an intense noise is not enough to generate a greater response of the ecosystem. In particular, experimental data are characterized by a CV approximately equal to 0.361, which corresponds to values of D and \(\tau\) lying on the diagonal strip in Fig. 2 ranging from \((\tau ,D)=(0.5,10^4)\) to \((\tau ,D)=(365,10^7)\). Finally we note the presence of a noise suppression effect. High values of D, indeed, can generate slight effects when the correlation time \(\tau\) does not take on suitable values.

The results shown here are an extension of the previous work by Benincà et al.56. There, the authors analyse a simpler, less realistic model of two interacting populations, whose dynamics is affected by a randomly fluctuating temperature. In that case, moreover, the deterministic oscillations of the temperature are suppressed, and the system exhibits intrinsic Lotka-Volterra oscillations whose frequency match with the characteristic one(s) of the noise. On the contrary, here, the observed maximum response (see Fig. 3) cannot be interpreted as a synchronization effect, since our model does not present intrinsic Lotka-Volterra-like oscillations and the periodic population variability is only due to the deterministic forcing(s).

The nonmonotonic behaviour of the CV can be then interpreted as the signature of the intimate interplay between the ecological system and the noise. This interplay, indeed, has a pivotal role in both determining the dynamics of the populations and defining the characteristics of the ecosystem.

In Fig. 3 it can be observed that the value of \(\tau\) for which CV is maximum strongly depends on the noise intensity D. In particular, it is possible to note that the peaks in Fig. 3 move towards higher values of \(\tau\) as the noise intensity increases. Thus, Fig. 3 demonstrates that the maximum-response effect to the random fluctuations is sensitive to the noise intensity D.

However, it is important to underline that the response of the system to the noisy signal does not depend on the yearly oscillations induced by the deterministic forcings. Indeed, by considering constant the deterministic part of all external forcings (temperature, irradiance, wind and salinity), the non monotonic behaviour of CV with respect to both \(\tau\) and D is still present, provided that the populations are not extinct (plot not shown). In this scenario indeed, besides dinoflagellates, diatoms and nanoflagellates are practically extinct as well, exhibiting thus a constant vanishing variance. All the other strains, instead, present qualitatively the same nonmonotonicity with only slight differences (shift of the peaks and different mean values of the CV curves), probably due to the extinction of diatoms and nanoflagellates which causes relevant differences in the system dynamics. More specifically, the system’s response seems to depend on both the noise intensity and the correlation time (see Fig. 3).

In this scenario (absence of seasonal driving) we have studied the dependence on both parameters D and \(\tau\) of the probability density functions (PDFs) of the non-vanishing populations. In Fig. 4, the PDFs of bacteria (B1), picophytoplankton (P3), microzooplankton (Z5) and etherotrophic nanoflagellates (Z6) are plotted for \(\tau =0.5\) and eight different values of the parameter D.

Figure 4

Dependence of the probability density functions of non-vanishing populations on the parameter D for \(\tau =0.5\). The curves are normalized within the interval taken into account. For this reason the relative peaks of the curves in the bottom panels have different values compared to those of the top panels. However, the figure aims at showing the existence of the value of the noise intensity for which the system is more sensitive as well as the generation of a stationary out-of-equilibrium state induced by the noise.

We see that the mean value and the variance of these populations are strongly affected by the presence of random fluctuations in the irradiance. Specifically, as the noise intensity increases the mean values of picophytoplankton and bacteria concentrations exhibit a shift. In particular, the results indicate that picophytoplankton is disavantaged by the presence of a noisy component in the irradiance, which indeed tends to inhibit its ability to absorbe the solar light, slowing down its growth. As a consequence, since phytoplankton and bacteria compete for the same resources, as the former declines the latter are favoured, with a compensation mechanism which allows their predators (zooplankton populations) to be almost not affected by the noisy behaviour of the irradiance. Further, we note that for intermediate values of the noise intensity (\(D = 10^4 – 10^5\)) a maximum of the variance occurs (the PDFs are clearly spread on a wider range of values). Such an effect indicates that the noisy behaviour of irradiance strongly influences the whole ecosystem dynamics. Moreover, the nonmonotonic behaviour of the variance (its PDFs become larger and then tighter again as the noise intensity increases) indicates that the noise pushes the ecosystem away from equilibrium, driving it towards a non-equilibrium steady state. Finally, we note that the nonmonotonic behaviour of CV as a function of the noise intensity remains also in the presence of seasonal driving.

Figure 5

Coefficient of variation (\(CV=\sigma / \mu\)) of nine planktonic populations resulting from numerical integration of model equations plotted versus the considered values of D; different curves correspond to different values of the correlation time \(\tau\).

Figure 5 shows indeed the nonmonotonic response of the ecosystem to the change of D when the deterministic seasonal cycling of the four environmental parameters (temperature, irradiance, wind and salinity) is present. It is easy to observe that also in this instance the major noise-induced effect appears in nanoflagellates and diatoms with a percent standard deviation of 50\(\%\) and 100\(\%\), respectively. The coalescence of different curves (related to different values of \(\tau\)), as D decreases, is due to the fact that for \(D \rightarrow 0\) the impact of the noise is negligible and the evolution of the system practically resembles the deterministic one. On the contrary, for higher values of D remarkable differences arise and clear peaks of CV appear in the considered range of variation.

These plots show that, for a fixed value of \(\tau\), there exists a value of the noise intensity for which the planktonic concentrations are maximally spread around their mean values (corresponding to the maximum value of CV and then of the variance). Moreover, such a nonmonotonic behaviour suggests the presence of a resonance, which can be interpreted as the effect of the interplay between the nonlinearity of the system and the environmental random fluctuations.

Also in this case, the interplay between the two parameters D and \(\tau\) in determining and characterizing the dynamics of the ecosystem transparently emerges. The value of D corresponding to the maximum value of CV, indeed, basically depends on the specific value of \(\tau\).

Finally, we point out that the different dynamic scenarios identified by the D-\(\tau\) couples can be experienced by the system during the year, since the two parameters may seasonally vary depending on the different weather conditions. In other words, a seasonally varying noise (see Fig. 1c) may cause the nine populations explore different regions of the D-\(\tau\) space during the year. Therefore, the results reported in this paper can highlight the detectable yearly variability of a marine ecosystem which does not stem from the deterministic seasonal variation of environmental parameters.

### Effects on the organic carbon

In this subsection the effects of the irradiance noise on the biogechemistry are analysed. In Fig. 6 the dependence on \(\tau\) of both the CV [panel (a)] and the mean value concentration [panel (b)] of detritus, labile dissolved organic carbon (L-DOC), semi-labile dissolved organic carbon (SL-DOC) and gross primary production (GPP) are shown. All these biogeochemical properties are correlated with carbon cycling. Gross primary production is related to the amount of carbon entering in the ecosystem, and is related to the maximum energy available in the ecosystem progressively dissipated in the trophic web. Gross primary production is directly affected by light fluctuation and its CV shape is very similar to that of the irradiance, Fig. 2. We selected also detritus and DOC because they are important indicators for the carbon cycling dynamics and are related to the cycling of chemicals like heavy metals62. The different curves, related to different values of D, approach the same (vanishing) value for large \(\tau\). As previously discussed for the CV [Fig. 6(a)] of biomass concentrations, this circumstance is due to the fact that, in this case, the system dynamics tends to the deterministic case, characterized by a unique possible realization implying a vanishing standard deviation. For high correlation times thus the system is insensitive to the noise intensity. On the contrary, for small values of \(\tau\), different values of D lead to significant differences of the variance. In particular, detritus, L-DOC and SL-DOC exhibit a clear non-monotonic behaviour whose maximum value depends on the combined values of D-\(\tau\). Only the GPP presents a decreasing monotonic behaviour.

The dependence of the mean value concentration on \(\tau\), instead, is qualitatively the same for all the four parameters. Also in this case we can note a diversification with respect to D occurring at small \(\tau\) and a (deterministic) constant value arising for low (high) values of D (\(\tau\)).

These results manifest that not only the population dynamics, but also all the biogeochemical processes are profoundly affected by the presence of stochastic environmental variables. The values and the behaviour of the examined quantities are indeed determined by the intimate interplay between the intensity and the time correlation of the noise fluctuations.

Figure 6

(a) Coefficient of variation (\(CV=\sigma / \mu\)) and (b) mean value concentration (\(\mu\)) of detritus, labile dissolved organic carbon (L-DOC), semi-labile dissolved organic carbon (SL-DOC) and gross primary production (GPP) resulting from numerical integration of model equations plotted versus the considered values of \(\tau\); the different curves are related to different values of the correlation time D.