← News
Pollution Wins Pollution Wins Planet

The Hidden Order of Plankton: How Noise and Mixing Shape Ocean Life

Gulf of Naples
The Hidden Order of Plankton: How Noise and Mixing Shape Ocean Life
Correlation Length Key parameter
DNA, Microscopy, Satellites Data sources used
Clustering, Diversity, Biomass Patterns explained

Across the world’s oceans, from the sunlit surface to the twilight zone, life pulses in invisible waves. A single liter of seawater can contain tens of thousands of plankton—microscopic organisms that form the foundation of marine food webs, produce half the planet’s oxygen, and absorb vast quantities of carbon dioxide. Yet despite their global importance, plankton communities are maddeningly patchy: dense blooms appear and vanish like storms, separated by stretches of relative emptiness. For decades, ecologists have struggled to explain this patchiness—why some regions teem while others starve, why species distributions follow certain statistical patterns, and how local chaos gives rise to global regularity.

Now, a new study reveals that two simple forces—random fluctuations in growth rates and the slow mixing of ocean waters—can explain a surprising range of plankton patterns, from local diversity to large-scale patchiness. The researchers show that these dynamics generate a hidden length scale, a correlation length, that governs how far plankton patches extend before dissolving into the background. This single parameter, emerging from the interplay of biology and physics, unifies seemingly unrelated ecological laws and offers a minimal framework for understanding marine biodiversity.

The Science

The team, led by Giorgio Vittorio Visco and Samir Suweis at the University of Padua, developed a mathematical model that combines stochastic population dynamics with spatial diffusion. At its core is a deceptively simple equation describing how plankton abundance changes over time and space:

Here, $n_x$ is the density of plankton at location $x$, $\bar{\mu}(n_x)$ is the deterministic growth rate, and $\sigma \xi_x$ represents random environmental fluctuations—noise in growth caused by shifting temperatures, nutrient pulses, or viral outbreaks. The term $D \partial_x^2 n_x$ captures spatial diffusion, modeling how ocean currents and mixing spread plankton across space.

The deterministic growth rate $\bar{\mu}(n_x)$ itself has three components:

  • $\beta / n_x$: immigration or rare bloom events that sustain populations at low densities,
  • $g$: intrinsic growth at intermediate densities,
  • $-n_x / c$: density-dependent regulation (e.g., competition or grazing) that limits growth at high densities.

This structure produces a sigmoidal per-capita growth curve—positive at low and medium abundances, negative at high ones—mirroring real ecological constraints.

Crucially, the model treats environmental noise as multiplicative: fluctuations scale with population size, leading to a characteristic $\Sigma(n) \sim n^2$ scaling in abundance variance. This is a signature of environmental stochasticity, distinct from demographic noise, and it predicts log-normal-like distributions over time—a pattern long observed in plankton time series.

To test their theory, the researchers analyzed three independent datasets:

  1. LTER-MC: A 20-year time series (1995–2015) of plankton abundance from the Long-Term Ecological Research station in the Gulf of Naples, tracking diatoms, dinoflagellates, and coccolithophores.
  2. GRUMP: The Global rRNA Universal Metabarcoding Plankton database, containing genetic data from over 1,000 oceanic samples, used to reconstruct species abundance distributions (SADs) for prokaryotes and eukaryotes.
  3. Tara Oceans expeditions: High-resolution spatial transects of chlorophyll concentration—a proxy for total phytoplankton biomass—across thousands of kilometers in the Pacific and other oceans.

These datasets span temporal, taxonomic, and spatial scales, allowing the team to probe whether a single mechanism could explain patterns across them all.

What They Found

The first test was whether growth-rate fluctuations follow the predicted $n^2$ scaling. Using the LTER-MC time series, the researchers computed the variance of abundance changes over a 7-day interval, conditioned on initial population size. Across all three plankton groups—diatoms, dinoflagellates, and coccolithophores—they found a clear power-law relationship: $\Sigma(n) \propto n^2$. This held not only for individual species but also for total community abundance (

Figure 3: Population distributions. In all panels, solid curves denote the best-fit Generalized Inverse Gaussian (GIG) distribution, whereas dashed curves show the best-fit log-normal distribution. (A): Species abundance distributions (SADs) for representative eukaryotic and prokaryotic samples from the GRUMP dataset. The inset summarizes the transition between the GIG and log-normal regimes as a function of λ\lambda. Triangles indicate the median λ\lambda values, and the gray region marks the interval where the two distributions are statistically indistinguishable according to the Bayesian Information Criterion (|ΔBIC|<5|\Delta_{\mathrm{BIC}}|<5). (B): Time-aggregated SADs for diatoms, dinoflagellates, and coccolithophores in the LTER-MC dataset. For coccolithophores, we have ΔBIC=0.08\Delta_{\mathrm{BIC}}=0.08, suggesting that lognormal and GIG are equivalent. (C and D): Examples of total abundance distributions. Panel C shows the distribution of chlorophyll concentration from the Tara Microbiome/Tara Pacific dataset, whereas panel D shows the total abundance distributions (TADs) for the LTER-MC dataset.
Figure 3: Population distributions. In all panels, solid curves denote the best-fit Generalized Inverse Gaussian (GIG) distribution, whereas dashed curves show the best-fit log-normal distribution. (A): Species abundance distributions (SADs) for representative eukaryotic and prokaryotic samples from the GRUMP dataset. The inset summarizes the transition between the GIG and log-normal regimes as a function of λ\lambda. Triangles indicate the median λ\lambda values, and the gray region marks the interval where the two distributions are statistically indistinguishable according to the Bayesian Information Criterion (|ΔBIC|<5|\Delta_{\mathrm{BIC}}|<5). (B): Time-aggregated SADs for diatoms, dinoflagellates, and coccolithophores in the LTER-MC dataset. For coccolithophores, we have ΔBIC=0.08\Delta_{\mathrm{BIC}}=0.08, suggesting that lognormal and GIG are equivalent. (C and D): Examples of total abundance distributions. Panel C shows the distribution of chlorophyll concentration from the Tara Microbiome/Tara Pacific dataset, whereas panel D shows the total abundance distributions (TADs) for the LTER-MC dataset. Source: Giorgio Vittorio Visco, Kobe Simoens

D,E), confirming that environmental noise dominates over demographic stochasticity.

Next, they examined species abundance distributions (SADs). In the absence of diffusion, the model predicts a Generalized Inverse Gaussian (GIG) distribution—a flexible form that can resemble power laws or log-normals depending on parameters. With diffusion, the effective noise is rescaled, altering the shape of the distribution.

The key parameter is $\lambda = 2(1 - g / \sigma^{\star 2})$, where $\sigma^\star$ is the effective noise strength after accounting for spatial mixing. The theory predicts a transition: when $\lambda \gtrsim 1.5$, the SAD follows a power-law-like GIG; when $\lambda \lesssim 1.5$, it becomes log-normal-like.

This prediction was borne out in the data. In the GRUMP dataset, eukaryotic plankton (e.g., diatoms, dinoflagellates) had $\lambda = 1.9$ and were best fit by a GIG distribution, while prokaryotes (e.g., bacteria) had $\lambda = 1.1$ and followed a log-normal (

Figure 4: Spatial patterns. (A): Spatial Taylor’s law for prokaryotic GRUMP samples from transect KM1906 near Hawaii. Triangles represent the data, and the solid line is the theoretical prediction from Eq. 7 using the fitted correlation length and local variance. The vertical dashed line marks the crossover scale L=2​ΓL=2\Gamma, separating the correlated and de-correlated regions. The inset shows the spatial two-point correlation function, with triangular markers denoting the median correlation within distance bins and the red line the exponential fit. (B): Comparison across all transects between the observed Taylor’s law exponent α\alpha and the theoretical prediction derived from the fitted correlation length. The inset compares the correlation lengths Γ\Gamma estimated from the empirical correlation function and from Eq. 7. (C): Conditional variance V​a​r​(n|L)Var(n|L) as a function of LL for chlorophyll concentration in a representative Tara Microbiome–Tara Pacific transect. Dots represent the data, and the solid line is the fit of Eq. 8. The inset shows the corresponding spatial correlation function. (D): Comparison between the observed patchiness exponent pp and the theoretical prediction from Eq. 8 over all samples. Gray dots denote individual samples, blue squares are median values within bins, and the blue line is the best linear fit. The inset compares the correlation lengths Γ\Gamma estimated from Eq. 8 and from the empirical spatial correlation function.
Figure 4: Spatial patterns. (A): Spatial Taylor’s law for prokaryotic GRUMP samples from transect KM1906 near Hawaii. Triangles represent the data, and the solid line is the theoretical prediction from Eq. 7 using the fitted correlation length and local variance. The vertical dashed line marks the crossover scale L=2​ΓL=2\Gamma, separating the correlated and de-correlated regions. The inset shows the spatial two-point correlation function, with triangular markers denoting the median correlation within distance bins and the red line the exponential fit. (B): Comparison across all transects between the observed Taylor’s law exponent α\alpha and the theoretical prediction derived from the fitted correlation length. The inset compares the correlation lengths Γ\Gamma estimated from the empirical correlation function and from Eq. 7. (C): Conditional variance V​a​r​(n|L)Var(n|L) as a function of LL for chlorophyll concentration in a representative Tara Microbiome–Tara Pacific transect. Dots represent the data, and the solid line is the fit of Eq. 8. The inset shows the corresponding spatial correlation function. (D): Comparison between the observed patchiness exponent pp and the theoretical prediction from Eq. 8 over all samples. Gray dots denote individual samples, blue squares are median values within bins, and the blue line is the best linear fit. The inset compares the correlation lengths Γ\Gamma estimated from Eq. 8 and from the empirical spatial correlation function. Source: Giorgio Vittorio Visco, Kobe Simoens

A). In the LTER-MC dataset, all three groups had $\lambda < 1.5$ and $\Delta_{\textrm{BIC}} \leq 0$, indicating log-normal-like SADs (

Figure 4: Spatial patterns. (A): Spatial Taylor’s law for prokaryotic GRUMP samples from transect KM1906 near Hawaii. Triangles represent the data, and the solid line is the theoretical prediction from Eq. 7 using the fitted correlation length and local variance. The vertical dashed line marks the crossover scale L=2​ΓL=2\Gamma, separating the correlated and de-correlated regions. The inset shows the spatial two-point correlation function, with triangular markers denoting the median correlation within distance bins and the red line the exponential fit. (B): Comparison across all transects between the observed Taylor’s law exponent α\alpha and the theoretical prediction derived from the fitted correlation length. The inset compares the correlation lengths Γ\Gamma estimated from the empirical correlation function and from Eq. 7. (C): Conditional variance V​a​r​(n|L)Var(n|L) as a function of LL for chlorophyll concentration in a representative Tara Microbiome–Tara Pacific transect. Dots represent the data, and the solid line is the fit of Eq. 8. The inset shows the corresponding spatial correlation function. (D): Comparison between the observed patchiness exponent pp and the theoretical prediction from Eq. 8 over all samples. Gray dots denote individual samples, blue squares are median values within bins, and the blue line is the best linear fit. The inset compares the correlation lengths Γ\Gamma estimated from Eq. 8 and from the empirical spatial correlation function.
Figure 4: Spatial patterns. (A): Spatial Taylor’s law for prokaryotic GRUMP samples from transect KM1906 near Hawaii. Triangles represent the data, and the solid line is the theoretical prediction from Eq. 7 using the fitted correlation length and local variance. The vertical dashed line marks the crossover scale L=2​ΓL=2\Gamma, separating the correlated and de-correlated regions. The inset shows the spatial two-point correlation function, with triangular markers denoting the median correlation within distance bins and the red line the exponential fit. (B): Comparison across all transects between the observed Taylor’s law exponent α\alpha and the theoretical prediction derived from the fitted correlation length. The inset compares the correlation lengths Γ\Gamma estimated from the empirical correlation function and from Eq. 7. (C): Conditional variance V​a​r​(n|L)Var(n|L) as a function of LL for chlorophyll concentration in a representative Tara Microbiome–Tara Pacific transect. Dots represent the data, and the solid line is the fit of Eq. 8. The inset shows the corresponding spatial correlation function. (D): Comparison between the observed patchiness exponent pp and the theoretical prediction from Eq. 8 over all samples. Gray dots denote individual samples, blue squares are median values within bins, and the blue line is the best linear fit. The inset compares the correlation lengths Γ\Gamma estimated from Eq. 8 and from the empirical spatial correlation function. Source: Giorgio Vittorio Visco, Kobe Simoens

B).

For total abundance, chlorophyll data from the Tara expeditions showed log-normal distributions across most transects (

Figure 4: Spatial patterns. (A): Spatial Taylor’s law for prokaryotic GRUMP samples from transect KM1906 near Hawaii. Triangles represent the data, and the solid line is the theoretical prediction from Eq. 7 using the fitted correlation length and local variance. The vertical dashed line marks the crossover scale L=2​ΓL=2\Gamma, separating the correlated and de-correlated regions. The inset shows the spatial two-point correlation function, with triangular markers denoting the median correlation within distance bins and the red line the exponential fit. (B): Comparison across all transects between the observed Taylor’s law exponent α\alpha and the theoretical prediction derived from the fitted correlation length. The inset compares the correlation lengths Γ\Gamma estimated from the empirical correlation function and from Eq. 7. (C): Conditional variance V​a​r​(n|L)Var(n|L) as a function of LL for chlorophyll concentration in a representative Tara Microbiome–Tara Pacific transect. Dots represent the data, and the solid line is the fit of Eq. 8. The inset shows the corresponding spatial correlation function. (D): Comparison between the observed patchiness exponent pp and the theoretical prediction from Eq. 8 over all samples. Gray dots denote individual samples, blue squares are median values within bins, and the blue line is the best linear fit. The inset compares the correlation lengths Γ\Gamma estimated from Eq. 8 and from the empirical spatial correlation function.
Figure 4: Spatial patterns. (A): Spatial Taylor’s law for prokaryotic GRUMP samples from transect KM1906 near Hawaii. Triangles represent the data, and the solid line is the theoretical prediction from Eq. 7 using the fitted correlation length and local variance. The vertical dashed line marks the crossover scale L=2​ΓL=2\Gamma, separating the correlated and de-correlated regions. The inset shows the spatial two-point correlation function, with triangular markers denoting the median correlation within distance bins and the red line the exponential fit. (B): Comparison across all transects between the observed Taylor’s law exponent α\alpha and the theoretical prediction derived from the fitted correlation length. The inset compares the correlation lengths Γ\Gamma estimated from the empirical correlation function and from Eq. 7. (C): Conditional variance V​a​r​(n|L)Var(n|L) as a function of LL for chlorophyll concentration in a representative Tara Microbiome–Tara Pacific transect. Dots represent the data, and the solid line is the fit of Eq. 8. The inset shows the corresponding spatial correlation function. (D): Comparison between the observed patchiness exponent pp and the theoretical prediction from Eq. 8 over all samples. Gray dots denote individual samples, blue squares are median values within bins, and the blue line is the best linear fit. The inset compares the correlation lengths Γ\Gamma estimated from Eq. 8 and from the empirical spatial correlation function. Source: Giorgio Vittorio Visco, Kobe Simoens

C), consistent with spatial variability driven by fluctuating growth rates.

The most striking validation came from spatial patterns. The model predicts that the two-point correlation function—the similarity in abundance between two locations—decays exponentially with distance:

where $\Gamma$ is the correlation length. Fitting this to GRUMP transects, the researchers found excellent agreement (

Figure 5: Simulations and analytical predictions. Numerical simulations of Eq. 1 are run on top of a one-dimensional lattice with L=170L=170 sites, u=1u=1 and periodic boundary conditions. We fix β=0.3\beta=0.3, c=105c=10^{5} and D=1D=1, while the other parameters are chosen to explore different values of λ\lambda. (A): we show the two-point correlation function with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. The predicted correlation length is Γ=11.2\Gamma=11.2 . The solid line represents the theoretical curve given by Eq. 4, while dots results from the simulation. (B): comparison between the fitted correlation length and the theoretical one for 140 simulations performed for different values of gg and σ\sigma. These combinations are selected to have 20 different values of λ\lambda between 1 and 2, and 7 values of Γ\Gamma spanning the range from 3 to 14. (C and (D)): spatial Taylor’s law (Eq. 7) and empirical variance (Eq. 8) in log-log scale with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. Solid lines are the analytical prediction and dots represent the simulation outcomes. (E): stationary abundance distribution for the two population regimes. In red, we report the log-normal regime for which λ=1.15\lambda=1.15 and ΔBIC<0\Delta_{\textrm{BIC}}<0. The parameters are g=3.5×10−3g=3.5\times 10^{-3} and σ=0.43\sigma=0.43, and the correlation length is Γ=11.2\Gamma=11.2. Red dots are the simulated data, while the curve follows Eq. 58. In blue, we show the GIG regime with λ=1.84\lambda=1.84 and ΔBIC>0\Delta_{\textrm{BIC}}>0. For this case, we choose g=7.4×10−4g=7.4\times 10^{-4} and σ=0.46\sigma=0.46. The value of Γ\Gamma is identical to that used in the log-normal regime. Here as well, the continuous line denotes the theoretical prediction, whereas the dots indicate the outcomes of the numerical simulations. Panel F: ΔBIC\Delta_{\textrm{BIC}} as a function of the fitted slope λ\lambda over 140 simulations. Triangles are median values of ΔBIC\Delta_{\textrm{BIC}}, and the gray band is the neutral interval |ΔBIC|<5|\Delta_{\textrm{BIC}}|<5.
Figure 5: Simulations and analytical predictions. Numerical simulations of Eq. 1 are run on top of a one-dimensional lattice with L=170L=170 sites, u=1u=1 and periodic boundary conditions. We fix β=0.3\beta=0.3, c=105c=10^{5} and D=1D=1, while the other parameters are chosen to explore different values of λ\lambda. (A): we show the two-point correlation function with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. The predicted correlation length is Γ=11.2\Gamma=11.2 . The solid line represents the theoretical curve given by Eq. 4, while dots results from the simulation. (B): comparison between the fitted correlation length and the theoretical one for 140 simulations performed for different values of gg and σ\sigma. These combinations are selected to have 20 different values of λ\lambda between 1 and 2, and 7 values of Γ\Gamma spanning the range from 3 to 14. (C and (D)): spatial Taylor’s law (Eq. 7) and empirical variance (Eq. 8) in log-log scale with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. Solid lines are the analytical prediction and dots represent the simulation outcomes. (E): stationary abundance distribution for the two population regimes. In red, we report the log-normal regime for which λ=1.15\lambda=1.15 and ΔBIC<0\Delta_{\textrm{BIC}}<0. The parameters are g=3.5×10−3g=3.5\times 10^{-3} and σ=0.43\sigma=0.43, and the correlation length is Γ=11.2\Gamma=11.2. Red dots are the simulated data, while the curve follows Eq. 58. In blue, we show the GIG regime with λ=1.84\lambda=1.84 and ΔBIC>0\Delta_{\textrm{BIC}}>0. For this case, we choose g=7.4×10−4g=7.4\times 10^{-4} and σ=0.46\sigma=0.46. The value of Γ\Gamma is identical to that used in the log-normal regime. Here as well, the continuous line denotes the theoretical prediction, whereas the dots indicate the outcomes of the numerical simulations. Panel F: ΔBIC\Delta_{\textrm{BIC}} as a function of the fitted slope λ\lambda over 140 simulations. Triangles are median values of ΔBIC\Delta_{\textrm{BIC}}, and the gray band is the neutral interval |ΔBIC|<5|\Delta_{\textrm{BIC}}|<5. Source: Giorgio Vittorio Visco, Kobe Simoens

A inset). Moreover, $\Gamma$ could be estimated independently from spatial Taylor’s law—the scaling of local variance with mean abundance across different spatial scales.

Taylor’s law in ecology usually takes the form $\text{Var}(n) \propto \langle n \rangle^\alpha$. In spatial contexts, the exponent $\alpha$ changes with the size $L$ of the sampling region. The model predicts a crossover: for $L \ll \Gamma$, $\alpha \approx 2$ (correlated regime); for $L \gg \Gamma$, $\alpha \approx 1$ (uncorrelated regime). The transition occurs at $L = 2\Gamma$.

Testing this on prokaryotic data from the KM1906 transect near Hawaii, the researchers observed exactly this crossover (

Figure 5: Simulations and analytical predictions. Numerical simulations of Eq. 1 are run on top of a one-dimensional lattice with L=170L=170 sites, u=1u=1 and periodic boundary conditions. We fix β=0.3\beta=0.3, c=105c=10^{5} and D=1D=1, while the other parameters are chosen to explore different values of λ\lambda. (A): we show the two-point correlation function with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. The predicted correlation length is Γ=11.2\Gamma=11.2 . The solid line represents the theoretical curve given by Eq. 4, while dots results from the simulation. (B): comparison between the fitted correlation length and the theoretical one for 140 simulations performed for different values of gg and σ\sigma. These combinations are selected to have 20 different values of λ\lambda between 1 and 2, and 7 values of Γ\Gamma spanning the range from 3 to 14. (C and (D)): spatial Taylor’s law (Eq. 7) and empirical variance (Eq. 8) in log-log scale with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. Solid lines are the analytical prediction and dots represent the simulation outcomes. (E): stationary abundance distribution for the two population regimes. In red, we report the log-normal regime for which λ=1.15\lambda=1.15 and ΔBIC<0\Delta_{\textrm{BIC}}<0. The parameters are g=3.5×10−3g=3.5\times 10^{-3} and σ=0.43\sigma=0.43, and the correlation length is Γ=11.2\Gamma=11.2. Red dots are the simulated data, while the curve follows Eq. 58. In blue, we show the GIG regime with λ=1.84\lambda=1.84 and ΔBIC>0\Delta_{\textrm{BIC}}>0. For this case, we choose g=7.4×10−4g=7.4\times 10^{-4} and σ=0.46\sigma=0.46. The value of Γ\Gamma is identical to that used in the log-normal regime. Here as well, the continuous line denotes the theoretical prediction, whereas the dots indicate the outcomes of the numerical simulations. Panel F: ΔBIC\Delta_{\textrm{BIC}} as a function of the fitted slope λ\lambda over 140 simulations. Triangles are median values of ΔBIC\Delta_{\textrm{BIC}}, and the gray band is the neutral interval |ΔBIC|<5|\Delta_{\textrm{BIC}}|<5.
Figure 5: Simulations and analytical predictions. Numerical simulations of Eq. 1 are run on top of a one-dimensional lattice with L=170L=170 sites, u=1u=1 and periodic boundary conditions. We fix β=0.3\beta=0.3, c=105c=10^{5} and D=1D=1, while the other parameters are chosen to explore different values of λ\lambda. (A): we show the two-point correlation function with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. The predicted correlation length is Γ=11.2\Gamma=11.2 . The solid line represents the theoretical curve given by Eq. 4, while dots results from the simulation. (B): comparison between the fitted correlation length and the theoretical one for 140 simulations performed for different values of gg and σ\sigma. These combinations are selected to have 20 different values of λ\lambda between 1 and 2, and 7 values of Γ\Gamma spanning the range from 3 to 14. (C and (D)): spatial Taylor’s law (Eq. 7) and empirical variance (Eq. 8) in log-log scale with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. Solid lines are the analytical prediction and dots represent the simulation outcomes. (E): stationary abundance distribution for the two population regimes. In red, we report the log-normal regime for which λ=1.15\lambda=1.15 and ΔBIC<0\Delta_{\textrm{BIC}}<0. The parameters are g=3.5×10−3g=3.5\times 10^{-3} and σ=0.43\sigma=0.43, and the correlation length is Γ=11.2\Gamma=11.2. Red dots are the simulated data, while the curve follows Eq. 58. In blue, we show the GIG regime with λ=1.84\lambda=1.84 and ΔBIC>0\Delta_{\textrm{BIC}}>0. For this case, we choose g=7.4×10−4g=7.4\times 10^{-4} and σ=0.46\sigma=0.46. The value of Γ\Gamma is identical to that used in the log-normal regime. Here as well, the continuous line denotes the theoretical prediction, whereas the dots indicate the outcomes of the numerical simulations. Panel F: ΔBIC\Delta_{\textrm{BIC}} as a function of the fitted slope λ\lambda over 140 simulations. Triangles are median values of ΔBIC\Delta_{\textrm{BIC}}, and the gray band is the neutral interval |ΔBIC|<5|\Delta_{\textrm{BIC}}|<5. Source: Giorgio Vittorio Visco, Kobe Simoens

A). The theoretical prediction, using $\Gamma$ estimated from the correlation function, matched the data without free parameters. Across all transects, the observed and predicted Taylor exponents were tightly correlated (

Figure 5: Simulations and analytical predictions. Numerical simulations of Eq. 1 are run on top of a one-dimensional lattice with L=170L=170 sites, u=1u=1 and periodic boundary conditions. We fix β=0.3\beta=0.3, c=105c=10^{5} and D=1D=1, while the other parameters are chosen to explore different values of λ\lambda. (A): we show the two-point correlation function with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. The predicted correlation length is Γ=11.2\Gamma=11.2 . The solid line represents the theoretical curve given by Eq. 4, while dots results from the simulation. (B): comparison between the fitted correlation length and the theoretical one for 140 simulations performed for different values of gg and σ\sigma. These combinations are selected to have 20 different values of λ\lambda between 1 and 2, and 7 values of Γ\Gamma spanning the range from 3 to 14. (C and (D)): spatial Taylor’s law (Eq. 7) and empirical variance (Eq. 8) in log-log scale with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. Solid lines are the analytical prediction and dots represent the simulation outcomes. (E): stationary abundance distribution for the two population regimes. In red, we report the log-normal regime for which λ=1.15\lambda=1.15 and ΔBIC<0\Delta_{\textrm{BIC}}<0. The parameters are g=3.5×10−3g=3.5\times 10^{-3} and σ=0.43\sigma=0.43, and the correlation length is Γ=11.2\Gamma=11.2. Red dots are the simulated data, while the curve follows Eq. 58. In blue, we show the GIG regime with λ=1.84\lambda=1.84 and ΔBIC>0\Delta_{\textrm{BIC}}>0. For this case, we choose g=7.4×10−4g=7.4\times 10^{-4} and σ=0.46\sigma=0.46. The value of Γ\Gamma is identical to that used in the log-normal regime. Here as well, the continuous line denotes the theoretical prediction, whereas the dots indicate the outcomes of the numerical simulations. Panel F: ΔBIC\Delta_{\textrm{BIC}} as a function of the fitted slope λ\lambda over 140 simulations. Triangles are median values of ΔBIC\Delta_{\textrm{BIC}}, and the gray band is the neutral interval |ΔBIC|<5|\Delta_{\textrm{BIC}}|<5.
Figure 5: Simulations and analytical predictions. Numerical simulations of Eq. 1 are run on top of a one-dimensional lattice with L=170L=170 sites, u=1u=1 and periodic boundary conditions. We fix β=0.3\beta=0.3, c=105c=10^{5} and D=1D=1, while the other parameters are chosen to explore different values of λ\lambda. (A): we show the two-point correlation function with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. The predicted correlation length is Γ=11.2\Gamma=11.2 . The solid line represents the theoretical curve given by Eq. 4, while dots results from the simulation. (B): comparison between the fitted correlation length and the theoretical one for 140 simulations performed for different values of gg and σ\sigma. These combinations are selected to have 20 different values of λ\lambda between 1 and 2, and 7 values of Γ\Gamma spanning the range from 3 to 14. (C and (D)): spatial Taylor’s law (Eq. 7) and empirical variance (Eq. 8) in log-log scale with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. Solid lines are the analytical prediction and dots represent the simulation outcomes. (E): stationary abundance distribution for the two population regimes. In red, we report the log-normal regime for which λ=1.15\lambda=1.15 and ΔBIC<0\Delta_{\textrm{BIC}}<0. The parameters are g=3.5×10−3g=3.5\times 10^{-3} and σ=0.43\sigma=0.43, and the correlation length is Γ=11.2\Gamma=11.2. Red dots are the simulated data, while the curve follows Eq. 58. In blue, we show the GIG regime with λ=1.84\lambda=1.84 and ΔBIC>0\Delta_{\textrm{BIC}}>0. For this case, we choose g=7.4×10−4g=7.4\times 10^{-4} and σ=0.46\sigma=0.46. The value of Γ\Gamma is identical to that used in the log-normal regime. Here as well, the continuous line denotes the theoretical prediction, whereas the dots indicate the outcomes of the numerical simulations. Panel F: ΔBIC\Delta_{\textrm{BIC}} as a function of the fitted slope λ\lambda over 140 simulations. Triangles are median values of ΔBIC\Delta_{\textrm{BIC}}, and the gray band is the neutral interval |ΔBIC|<5|\Delta_{\textrm{BIC}}|<5. Source: Giorgio Vittorio Visco, Kobe Simoens

B).

Finally, they tested “patchiness”—how variance scales with spatial scale. The model predicts:

This yields $p \approx 2$ for small $L$ and $p \approx 1$ for large $L$. Empirical patchiness exponents matched this curve across all samples (

Figure 5: Simulations and analytical predictions. Numerical simulations of Eq. 1 are run on top of a one-dimensional lattice with L=170L=170 sites, u=1u=1 and periodic boundary conditions. We fix β=0.3\beta=0.3, c=105c=10^{5} and D=1D=1, while the other parameters are chosen to explore different values of λ\lambda. (A): we show the two-point correlation function with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. The predicted correlation length is Γ=11.2\Gamma=11.2 . The solid line represents the theoretical curve given by Eq. 4, while dots results from the simulation. (B): comparison between the fitted correlation length and the theoretical one for 140 simulations performed for different values of gg and σ\sigma. These combinations are selected to have 20 different values of λ\lambda between 1 and 2, and 7 values of Γ\Gamma spanning the range from 3 to 14. (C and (D)): spatial Taylor’s law (Eq. 7) and empirical variance (Eq. 8) in log-log scale with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. Solid lines are the analytical prediction and dots represent the simulation outcomes. (E): stationary abundance distribution for the two population regimes. In red, we report the log-normal regime for which λ=1.15\lambda=1.15 and ΔBIC<0\Delta_{\textrm{BIC}}<0. The parameters are g=3.5×10−3g=3.5\times 10^{-3} and σ=0.43\sigma=0.43, and the correlation length is Γ=11.2\Gamma=11.2. Red dots are the simulated data, while the curve follows Eq. 58. In blue, we show the GIG regime with λ=1.84\lambda=1.84 and ΔBIC>0\Delta_{\textrm{BIC}}>0. For this case, we choose g=7.4×10−4g=7.4\times 10^{-4} and σ=0.46\sigma=0.46. The value of Γ\Gamma is identical to that used in the log-normal regime. Here as well, the continuous line denotes the theoretical prediction, whereas the dots indicate the outcomes of the numerical simulations. Panel F: ΔBIC\Delta_{\textrm{BIC}} as a function of the fitted slope λ\lambda over 140 simulations. Triangles are median values of ΔBIC\Delta_{\textrm{BIC}}, and the gray band is the neutral interval |ΔBIC|<5|\Delta_{\textrm{BIC}}|<5.
Figure 5: Simulations and analytical predictions. Numerical simulations of Eq. 1 are run on top of a one-dimensional lattice with L=170L=170 sites, u=1u=1 and periodic boundary conditions. We fix β=0.3\beta=0.3, c=105c=10^{5} and D=1D=1, while the other parameters are chosen to explore different values of λ\lambda. (A): we show the two-point correlation function with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. The predicted correlation length is Γ=11.2\Gamma=11.2 . The solid line represents the theoretical curve given by Eq. 4, while dots results from the simulation. (B): comparison between the fitted correlation length and the theoretical one for 140 simulations performed for different values of gg and σ\sigma. These combinations are selected to have 20 different values of λ\lambda between 1 and 2, and 7 values of Γ\Gamma spanning the range from 3 to 14. (C and (D)): spatial Taylor’s law (Eq. 7) and empirical variance (Eq. 8) in log-log scale with g=2.3×10−3g=2.3\times 10^{-3} and σ=0.45\sigma=0.45. Solid lines are the analytical prediction and dots represent the simulation outcomes. (E): stationary abundance distribution for the two population regimes. In red, we report the log-normal regime for which λ=1.15\lambda=1.15 and ΔBIC<0\Delta_{\textrm{BIC}}<0. The parameters are g=3.5×10−3g=3.5\times 10^{-3} and σ=0.43\sigma=0.43, and the correlation length is Γ=11.2\Gamma=11.2. Red dots are the simulated data, while the curve follows Eq. 58. In blue, we show the GIG regime with λ=1.84\lambda=1.84 and ΔBIC>0\Delta_{\textrm{BIC}}>0. For this case, we choose g=7.4×10−4g=7.4\times 10^{-4} and σ=0.46\sigma=0.46. The value of Γ\Gamma is identical to that used in the log-normal regime. Here as well, the continuous line denotes the theoretical prediction, whereas the dots indicate the outcomes of the numerical simulations. Panel F: ΔBIC\Delta_{\textrm{BIC}} as a function of the fitted slope λ\lambda over 140 simulations. Triangles are median values of ΔBIC\Delta_{\textrm{BIC}}, and the gray band is the neutral interval |ΔBIC|<5|\Delta_{\textrm{BIC}}|<5. Source: Giorgio Vittorio Visco, Kobe Simoens

D), and the $\Gamma$ values derived from patchiness agreed with those from correlation functions.

λ Values in LTER-MC Plankton Groups

Fitted λ values for three plankton groups in the LTER-MC dataset, all below 1.5, indicating log-normal-like species abundance distributions.

λ Values in LTER-MC Plankton Groups
LabelValue
Diatoms1.3
Dinoflagellates1.2
Coccolithophores1.1

Source: LTER-MC dataset analysis (Visco et al., 2026).

Why This Changes Things

For decades, plankton ecologists have treated patterns like SADs, Taylor’s law, and patchiness as separate phenomena, each requiring its own explanation. Neutral theory explained SADs but ignored space. Turbulence models explained patchiness but not abundance distributions. This study bridges that gap.

The key insight is that a single mechanism—stochastic growth coupled with diffusion—generates a correlation length $\Gamma \approx 4D / \sigma^2$, which sets the scale of ecological coherence in the ocean. Below $\Gamma$, populations are correlated; above it, they fluctuate independently. This length scale emerges from the balance between mixing (which homogenizes) and environmental noise (which diversifies).

This has profound implications. First, it suggests that much of marine biodiversity can be understood as a “null model”—a baseline shaped by physics and chance, against which deviations (e.g., from predation, symbiosis, or adaptation) can be measured. If a region’s plankton patterns don’t fit this model, it may signal the action of additional biological forces.

Second, it offers a way to quantify ecosystem resilience. The correlation length $\Gamma$ depends on $D$ (diffusivity) and $\sigma$ (noise strength). Climate change is altering both: warming oceans may reduce vertical mixing (lowering $D$), while increased weather variability may amplify $\sigma$. According to the model, this could shrink $\Gamma$, leading to smaller, more fragmented plankton patches—a potential early warning of ecosystem destabilization.

Third, it reframes our view of marine “hotspots” of diversity. Rather than being solely the product of favorable conditions or evolutionary history, patchiness may arise simply from the physics of mixing and noise. This doesn’t diminish their importance, but it suggests that conservation efforts should consider not just where diversity is high, but how it is structured—whether as isolated oases or as part of a connected seascape.

The model also resolves a long-standing puzzle: why some plankton communities follow power laws while others are log-normal. The answer lies in the ratio of growth to noise. When environmental fluctuations are strong relative to intrinsic growth ($\lambda > 1.5$), power-law SADs emerge; when growth dominates, log-normals prevail. This explains why prokaryotes—fast-growing, generalist bacteria—tend to be log-normal, while eukaryotes—slower, more specialized—often follow power laws.

What’s Next

The model is minimal by design, and the authors acknowledge its limitations. It assumes one-dimensional dispersal, ignores advection (directed flow), and treats all species as statistically equivalent. Seasonality is not explicitly modeled, though the framework holds when data are analyzed in monthly blocks.

Future work could extend the model to two or three dimensions, incorporate seasonal forcing, or include species interactions like competition or predation. The framework could also be applied beyond plankton—to coral reefs, forests, or even microbial communities in the human gut—where stochastic growth and dispersal shape diversity.

Perhaps most importantly, this work opens the door to real-time monitoring of ocean health. If $\Gamma$ can be estimated from satellite chlorophyll data or autonomous ocean sensors, it could serve as a dynamic indicator of ecosystem state. A shrinking correlation length might signal increased fragmentation; a shift in $\lambda$ could reflect changing noise regimes.

In a world where marine ecosystems face unprecedented stress—from warming, acidification, and overfishing—having a simple, predictive theory for biodiversity is not just academically satisfying. It’s a tool for stewardship. By revealing the hidden order beneath the ocean’s apparent chaos, this study offers a compass for navigating the uncertain future of the seas.

λ Values in GRUMP Dataset by Domain

Comparison of λ values between prokaryotic and eukaryotic plankton in the GRUMP dataset, showing a clear divergence in abundance distribution shapes.

λ Values in GRUMP Dataset by Domain
LabelValue
Prokaryotes1.1
Eukaryotes1.9

Source: GRUMP dataset analysis (Visco et al., 2026).

View on map →