SPEAD 1.0 – Simulating Plankton Evolution with Adaptive Dynamics in a two-trait continuous fitness landscape applied to the Sargasso Sea

Diversity plays a key role in the adaptive capacity of marine ecosystems to environmental changes. However, modelling the adaptive dynamics of phytoplankton traits remains challenging due to the competitive exclusion of sub-optimal phenotypes and the complexity of evolutionary processes leading to optimal phenotypes. Trait diffusion (TD) is a recently developed approach to sustain diversity in plankton models by introducing mutations, therefore allowing the adaptive evolution of functional traits to occur at ecological timescales. In this study, we present a model called Simulating Plankton Evolution with Adaptive Dynamics (SPEAD) that resolves the eco-evolutionary processes of a multi-trait plankton community. The SPEAD model can be used to evaluate plankton adaptation to environmental changes at different timescales or address ecological issues affected by adaptive evolution. Phytoplankton phenotypes in SPEAD are characterized by two traits, the nitrogen half-saturation constant and optimal temperature, which can mutate at each generation using the TD mechanism. SPEAD does not resolve the different phenotypes as discrete entities, instead computing six aggregate properties: total phytoplankton biomass, the mean value of each trait, trait variances, and the inter-trait covariance of a single population in a continuous trait space. Therefore, SPEAD resolves the dynamics of the population's continuous trait distribution by solving its statistical moments, wherein the variances of trait values represent the diversity of ecotypes. The ecological model is coupled to a vertically resolved (1D) physical environment, and therefore the adaptive dynamics of the simulated phytoplankton population are driven by seasonal variations in vertical mixing, nutrient concentration, water temperature, and solar irradiance. The simulated bulk properties are validated by observations from Bermuda Atlantic Time-series Studies (BATS) in the Sargasso Sea. We find that moderate mutation rates sustain trait diversity at decadal timescales and soften the almost total inter-trait correlation induced by the environment alone, without reducing the annual primary production or promoting permanently maladapted phenotypes, as occur with high mutation rates. As a way to evaluate the performance of the continuous trait approximation, we also compare the solutions of SPEAD to the solutions of a classical discrete entities approach, with both approaches including TD as a mechanism to sustain trait variance. We only find minor discrepancies between the continuous model SPEAD and the discrete model, with the computational cost of SPEAD being lower by 2 orders of magnitude. Therefore, SPEAD should be an ideal eco-evolutionary plankton model to be coupled to a general circulation model (GCM) of the global ocean.

1 Introduction

Phytoplankton are a polyphyletic group of microscopic primary producers widespread in aquatic environments. They are mainly single-celled, although colonial or multicellular species exist in most phytoplankton phyla (Beardall et al., 2009). Despite accounting for only 1 % of the global photosynthetic biomass, phytoplankton perform more than 45 % of Earth's net primary production (Field et al., 1998; Falkowski et al., 2004). They are the basis of all oceanic food webs and play key roles in biogeochemical cycles (Falkowski, 2012). In particular, they have a large impact on global climate through the export of detritic carbon from the surface to the ocean interior, sequestrating carbon in the deep ocean for timescales from a few years to more than a millennium depending on the depth they reach (DeVries and Primeau, 2011; DeVries et al., 2012). This process, called the “biological carbon pump”, regulates the concentration of carbon dioxide in the atmosphere (Volk and Hoffert, 1985; Falkowski et al., 1998).

Phytoplankton are highly diverse and live in many different environments. They differ in their ecological interactions and the processes through which they mediate biogeochemical cycles. For instance, some species can fix atmospheric nitrogen and enrich oligotrophic regions; some produce ballast minerals (mainly silica and calcium carbonate) and sink faster to the deep ocean, and others are mixotrophic, being able to both photosynthesize and feed on organic sources (Le Quéré et al., 2005). Most species are denser than seawater and eventually sink, but some are buoyant (Lännergren, 1979; Villareal, 1988). Phytoplankton size ranges from less than 1 µm for cyanobacteria like Prochlorococcus (Chisholm et al., 1988) to more than 1 mm for the giant diatom Ethmodiscus rex (Swift, 1973; Villareal and Carpenter, 1994). Their half-saturation constants for the main limiting nutrients range over 3 orders of magnitude (Edwards et al., 2012), reflecting adaptation to different nutrient supply levels. They are also adapted to very different temperatures: some diatoms can grow within sea ice (Ackley and Sullivan, 1994), whereas some hyperthermophilic cyanobacteria can grow at up to 75 ^∘C in hot springs (Castenholz, 1969). Most oceanic species have an optimal temperature for growth between 0 and 35 ^∘C (Thomas et al., 2012). Even within the same species or genus, wide variability has been observed for key traits such as iron requirements (Strzepek and Harrison, 2004), light requirements (Biller et al., 2015), and resistance to predation (Yoshida et al., 2004). Given that any change in the abundance and composition of phytoplankton has far-reaching consequences for other organisms and for the Earth's climate, it is important to understand the factors driving the dynamics of such communities.

Numerical modelling studies can address this issue by finding the mechanistic equations and parameters that most correctly account for the observations, and thereby provide invaluable insights into the general rules controlling ecosystems. Models can also be used to make predictions of how phytoplankton will impact or be impacted by future environmental changes (Norberg et al., 2012; Irwin et al., 2015). Mathematical models of phytoplankton growth as a function of nutrient concentration, temperature, and radiation have been developed since the 1940s (Riley, 1946; Steele, 1958; Riley, 1965), leading to the now common NPZD (nutrient, phytoplankton, zooplankton, detritus) models representing one or several compartments of nutrients, phytoplankton, zooplankton, and detritus (Fasham et al., 1990). Since the early 1990s (Maier-Reimer, 1993), the increase in computational power has allowed biogeochemical models to be fully coupled with ocean circulation (Aumont et al., 2003; Follows et al., 2007). However, representing all the phytoplankton diversity in models is neither feasible nor desirable. The computational cost would be high, and even if computationally feasible, the existing observations would not suffice to constrain the many free parameters.

Instead, ecological models account for biodiversity through a few key traits representing physiological characteristics or adaptation to different environments. The most widely investigated phytoplankton traits are cell size, nutrient niche, optimal temperature, optimal irradiance, and resistance to predation. Some trait-based models divide the phytoplankton community into discrete entities or “boxes” with different traits. The boxes can be as simple as diatoms and small phytoplankton groups (Aumont et al., 2015), with diatoms having higher nutrient concentration niches, or include more complex divisions into functional groups (Baretta et al., 1995; Le Quéré et al., 2005; Follows et al., 2007). In this article, we will call these models “discrete” or “multi-phenotype”. The other approach, which further reduces the number of equations while still allowing communities to adapt to changes in their environments, is to consider one or several continuously distributed traits and to compute only the dynamics of aggregate properties, such as community biomass, mean trait values, and trait variances (Wirtz and Eckhardt, 1996; Norberg et al., 2001; Bruggeman and Kooijman, 2007; Merico et al., 2009; Acevedo-Trejos et al., 2016; Smith et al., 2016; Chen and Smith, 2018). In this method trait variance can be used as a quantitative index of biodiversity. A community with a higher trait variance is considered to be more diverse because it has a wider spread and more even distribution of trait values (Li, 1997), although it does not necessarily have a higher number of taxonomic species (richness) (Vallina et al., 2017). This second type of model will be called “aggregate” or “continuous trait”.

One weakness induced by the simplification of phytoplankton communities in both aggregate and discrete models, however, is that competitive exclusion (Hardin, 1960; Hutchinson, 1961) often leads to a collapse of the modelled diversity (Merico et al., 2009), making adaptation impossible unless trait variance is artificially imposed (Norberg et al., 2012; Wirtz, 2013) or a mechanism is explicitly added to sustain it. One way to sustain biodiversity is through immigration from a distant community (Norberg et al., 2001; Savage et al., 2007). Yet, immigration does not explain the diversity observed in closed laboratory experiments, including continuous cultures (Fussmann et al., 2007; Kinnison and Hairston, 2007; Beardmore et al., 2011). Biodiversity can also be sustained by viruses (Thingstad and Lignell, 1997) or predators (Murdoch, 1969; Kiørboe et al., 1996) if they specialize in a narrow range of prey or switch their preference to the most common phytoplankton species. This is the idea behind the “kill the winner” theory (Thingstad, 2000; Vallina et al., 2014b), whereby predation concentrating more on the most dominant species maintains diversity because then each prey species persists at the abundance at which the predation rate equals its growth rate.

An alternative approach recently introduced to sustain diversity in models is to allow the simulated phytoplankton to mutate their functional traits (Kremer and Klausmeier, 2013; Merico et al., 2014), as mutations are the ultimate cause of all diversity and adaptation. Due to their short generation times of around 1 d (Marañon et al., 2013), phytoplankton are known to evolve at the timescale of a few years (Schlüter et al., 2016). For phytoplankton, ecological timescales featuring successions of dominant species in reaction to changes in the environment and selection of the fittest overlap evolutionary timescales on which species can also evolve genetically to adapt to their new environment (Irwin et al., 2015). As far as we know, the first aggregate phytoplankton model allowing a phytoplankton trait to randomly mutate through subsequent generations, before being selected by the environment, was developed by Merico et al. (2014). They called their scheme trait diffusion (TD), where “diffusion” is a mathematical term referring to the spreading of a property, in this case the trait value, not to physical transport. Trait diffusion of a single physiological trait was recently introduced in a model coupled with oceanic circulation (Chen and Smith, 2018). Upgrading the trait diffusion framework to several traits requires more complex equations and the introduction of a new class of state variables: the covariances between traits. However, multi-trait models are more realistic, and conceptual modelling studies have shown that the dynamics of correlated traits sometimes differ from those of single-trait models (Savage et al., 2007).

There are several other types of plankton ecological (without mutations) or eco-evolutionary (with mutations) models, such as individual-based models, resident–mutant models, and models with stochastic mutations. Each represents and sustains diversity in its own unique way. All the modelling approaches mentioned in this paper and others are reviewed with their assumptions, costs, and benefits by Ward et al. (2019). For instance, the aggregate approach corresponds to the parametric trait distribution model in their Fig. 3 (eco-v), and the multi-phenotype approach encompasses the “everything is everywhere” (eco-iii), functional groups (eco-iv), and deterministic mutations (evo-iii) cases depending on the number of phenotypes and the presence or absence of mutations. What makes aggregate models particularly advantageous is their cost efficiency and their applicability to a spatial context. Adding adaptive evolution to models has been identified as a key challenge for the near future, as it will allow researchers to answer several unresolved ecological questions more fundamental than sustaining variance that cannot be answered by purely ecological models. For instance, eco-evolutionary models can be used to assess the relative roles of natural selection and neutral evolution in explaining the observed diversity patterns. They could also serve to explore hypotheses on the emergence of new species and new environments after a mass extinction. Finally, representing evolution is necessary to estimate the prevalence of extinction, adaptation, and migration in response to future environmental changes.

Here we present a new aggregate phytoplankton model called SPEAD (Simulating Plankton Evolution with Adaptive Dynamics), an eco-evolutionary model using the trait diffusion framework for two key phytoplankton traits: the nitrogen half-saturation constant and optimal temperature for growth. The SPEAD model is based on an NPZD model (Vallina et al., 2014a, 2017), wherein the phytoplankton compartment is represented by the community biomass, mean trait values, trait variances, and covariance. SPEAD is embedded in a 1D (water column) physical setting simulating the Sargasso Sea using data from the Bermuda Atlantic Time-series Studies (BATS). We chose the 1D rather than a simpler 0D setting because vertical turbulent diffusion (not to be confused with trait diffusion) is the main source of covariance by mixing communities from different depths. Since the trait diffusion equations can easily be adapted to a discrete model, we have also built a discrete version of SPEAD wherein the phytoplankton community consists of 625 different phenotypes (i.e. 25 half-saturation constants and 25 optimal temperatures), each characterized by its own fixed set of trait values. The discrete version is more intuitive and easier to programme, and it provides a useful control experiment.

In the following sections, we first describe our ecological model, the differential equations controlling the growth of phytoplankton, and the adaptive evolution of their trait distribution, as well as the physical model setting. Then, we present the model outputs. In order to validate SPEAD and to highlight its novelties, we will focus on answering the following four questions.

1. How well does SPEAD represent the bulk properties of phytoplankton communities observed in the Sargasso Sea?

2. Do the aggregate and discrete approaches agree?

3. How are phytoplankton dynamics changed by the value of the mutation rates?

4. Can the mean value and variance of each trait be represented independently by a one-trait model wherein only nitrogen half-saturation or optimal temperature varies between phenotypes?

Finally, we discuss the reach of our modelling framework, focusing on three aspects: the performance of aggregate models, the choice of phytoplankton traits, and the relationship between trait diffusion and evolution.

2 Methods

2.1 A phytoplankton community model with two traits

Our phytoplankton community model SPEAD extends an existing nitrogen-based NPZD model (Vallina et al., 2017). Nitrogen is partitioned into four pools, all expressed in millimoles of nitrogen per cubic metre (mmol N m⁻³): phytoplankton (P in the equations), zooplankton (Z), dissolved inorganic nitrogen or DIN (N), and particulate organic nitrogen or PON (D as in “detritus”). Phytoplankton increase biomass by taking up DIN (V_p). Zooplankton increase biomass by grazing phytoplankton (G_z). The non-predation mortalities of phytoplankton (M_p) and zooplankton (M_z) and the nitrogen exudation by zooplankton (E_z) are divided between DIN and PON. Given that nitrogen is the limiting nutrient for phytoplankton growth, we do not consider nitrogen exudation by phytoplankton. ω_p, ω_z, and ϵ_z are constants representing the respective proportions of M_p, M_z, and E_z going to DIN. PON is remineralized to DIN (M_d). The fluxes from one pool to another are controlled by the pool concentrations and by two environmental forcings: temperature (T; ^∘C) and photosynthetically active radiation or PAR (I; W m⁻²). The main state variables of the model and their relationships are shown in Table 1 and Fig. 1a, and the expressions of the biogeochemical fluxes between the pools (mmol N m⁻³ d⁻¹) are given by the following equations, with their dependencies.

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}=V_{\mathrm{p}}(P,N,T,I)-M_{\mathrm{p}}(P,T)-G_{\mathrm{z}}(P,Z,T)\text{(2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}=G_{\mathrm{z}}(P,Z,T)-E_{\mathrm{z}}(P,Z,T)-M_{\mathrm{z}}(Z,T){\displaystyle }& {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}}={\mathit{ϵ}}_{\mathrm{z}}{E}_{\mathrm{z}}(P,Z,T)+{\mathit{\omega }}_{\mathrm{p}}{M}_{\mathrm{p}}(P,T)+{\mathit{\omega }}_{\mathrm{z}}{M}_{\mathrm{z}}(Z,T)\\ \text{(3)}& {\displaystyle }& {\displaystyle }+M_{\mathrm{d}}(D,T)-V_{\mathrm{p}}(P,N,T,I){\displaystyle }& {\displaystyle \frac{\mathrm{d}D}{\mathrm{d}t}}=(\mathrm{1}-{\mathit{ϵ}}_{\mathrm{z}})E_{\mathrm{z}}(P,Z,T)+(\mathrm{1}-{\mathit{\omega }}_{\mathrm{p}})M_{\mathrm{p}}(P,T)\\ \text{(4)}& {\displaystyle }& {\displaystyle }+(\mathrm{1}-{\mathit{\omega }}_{\mathrm{z}})M_{\mathrm{z}}(Z,T)-M_{\mathrm{d}}(D,T)\end{array}$$

Table 1 State variables of the ecosystem model.

Figure 1 NPZD (nutrient, phytoplankton, zooplankton, detritus) model within its physical setting (a). The phytoplankton pool is represented by a discrete set of species with different traits (b) or by moments of the trait distributions, assuming a bivariate normal distribution (c). Colours in (b) and (c) represent phytoplankton concentration.

Zooplankton, DIN, and PON are generic pools characterized by a single variable: their concentration. Phytoplankton and zooplankton mortality, zooplankton exudation, grazing, and the particle remineralization rate have simple expressions as a function of the nitrogen pool concentrations and temperature.

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{z}}(P,Z,T)=g_{\mathrm{0}}{e}^{{\mathit{\alpha }}_{\mathrm{h}}(T-T_{\mathrm{0}})}{\displaystyle \frac{P^{\mathrm{2}}}{P^{\mathrm{2}}+K_{\mathrm{p}}^{\mathrm{2}}}}Z\text{(6)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{z}}(P,Z,T)=(\mathrm{1}-{\mathit{\beta }}_{\mathrm{z}})G_{\mathrm{z}}(P,Z,T)\text{(7)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{z}}(Z,T)={\mathit{\psi }}_{\mathrm{z}}{e}^{{\mathit{\alpha }}_{\mathrm{h}}(T-T_{\mathrm{0}})}{Z}^{\mathrm{2}}\text{(8)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{p}}(P,T)={\mathit{\psi }}_{\mathrm{p}}{e}^{{\mathit{\alpha }}_{\mathrm{h}}(T-T_{\mathrm{0}})}P\text{(9)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{d}}\left(T\right)={\mathit{\psi }}_{\mathrm{d}}{e}^{{\mathit{\alpha }}_{\mathrm{h}}(T-T_{\mathrm{0}})}D\end{array}$$

The constants appearing in this equation (α_h, T₀, K_p, β_z, ψ_z, ψ_p, and ψ_d) are described in Table 2. Zooplankton mortality depends on the square of zooplankton concentration in order to prevent an explosion of zooplankton concentration. This stabilizing quadratic mortality term represents consumption by animals higher on the trophic chain, which is expected to increase faster than a linear function of Z biomass. Grazing is formulated as a Holling type III function (Holling, 1959) of phytoplankton concentration, with a niche at low concentrations to prevent the whole phytoplankton community from going extinct, even when they have very low growth rates. Grazing, mortality, and remineralization are considered heterotrophic processes and as such increase exponentially with temperature. The exponential factor is α_h = 0.092 ^∘C⁻¹. This is equivalent to multiplying the speed of all these processes by 2.5 when the temperature increases by 10 ^∘C, as in a Q₁₀=2.5 formulation. This value of Q₁₀ is close to measured values for zooplankton grazing (Hansen et al., 1997) and to the theoretical predictions of the metabolic theory of ecology for respiration (Gillooly et al., 2001; Allen et al., 2005).

Table 2 Parameters of the ecosystem model.

In contrast, the phytoplankton pool is composed of diverse organisms responding to environmental conditions in different ways. The diversity of phytoplankton is represented by variations in the values of two traits: the logarithm of the half-saturation constant for nitrogen uptake (x) and the optimal temperature for growth (y). From now on, we will refer to each set (x,y) of trait values as a “phenotype”. Nutrient uptake by phytoplankton depends on the trait distribution. The bivariate trait distribution is represented by a density p(x,y) (mmol N m^{−3  ∘}C⁻¹) so that the biomass of phytoplankton (mmol N m⁻³) with trait values between x₁ and x₂ and between y₁ and y₂ is equal to ${\int }_{x_{\mathrm{1}}}^{x_{\mathrm{2}}}{\int }_{y_{\mathrm{1}}}^{y_{\mathrm{2}}}p(x,y)\mathrm{d}x\mathrm{d}y$, and by extension the total phytoplankton biomass P is equal to the density integrated over the whole trait domain. Any phenotype has its own uptake rate u_p(x,y) (d⁻¹). The uptake rate is the product of a constant ($u_{\mathrm{p}}^{\mathrm{0}}$) and three dimensionless growth factors: a nutrient factor (γ_n(N,x)), a temperature factor (γ_T(T,y)), and an irradiance factor (γ_I(I)). Two of these factors, γ_n(N,x) and γ_I(I), represent limitations by resources. The third factor, γ_T(T,y), represents the kinetic effect of temperature on growth. In this study, we use the Monod approach (Monod, 1949) so that cells do not store nutrients and the uptake rate is equal to the reproduction or growth rate. We assume that all phytoplankton are unicellular and we do not consider changes in their cell volumes, which allows us to use the words “growth” and “reproduction” interchangeably. All phenotypes share the same rates of mortality and grazing.

The last term in the equation of trait density (Eq. 10) is trait diffusion (TD), as defined by Merico et al. (2014). Trait diffusion represents the fact that offspring can exhibit different trait values than their parents due to mutations or otherwise heritable plasticity. In our numerical model, we assume only that these mutations are heritable and random. They can represent genetic mutations and other, e.g. epigenetic, phenotypic plasticity. We assume that mutations on x and y are independent of each other. In the limit of small but frequent mutations, stochasticity can be neglected (Dieckmann and Law, 1996; Champagnat et al., 2006), and this process can be represented as a deterministic diffusion depending on diffusivity parameters ν_x and ν_y and on the second derivatives of the growth rate (u_p(x,y)) relative to each trait. Note that in TD, diffusion is a mathematical term referring to the spreading of a property, in our case trait values, not to a physical mixing process. It should therefore not be confused with vertical turbulent diffusion, which is also present in our model (see Sect. 2.3). To avoid ambiguity, from now on, we will refer to the trait diffusivity parameters as “mutation rates”. ν_x and ν_y are mutation rates per generation, not per unit of time; therefore, time does not appear in their units. They have the same units as trait variances. The derivation of the trait diffusion term is explained in Appendix A. The differential equations followed by a given phenotype (x,y) are as follows.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial p(x,y,t)}{\partial t}}=\\ {\displaystyle }& {\displaystyle }\left[u_{\mathrm{p}}(N,T,I,x,y)-{\displaystyle \frac{M_{\mathrm{p}}(P,T)}{P}}-{\displaystyle \frac{G_{\mathrm{z}}(P,Z,T)}{P}}\right]\\ \text{(10)}& {\displaystyle }& {\displaystyle }p(x,y,t)+{\mathit{\nu }}_x{\displaystyle \frac{{\partial }^{\mathrm{2}}(u_{\mathrm{p}}\cdot p)}{\partial x^{\mathrm{2}}}}+{\mathit{\nu }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}(u_{\mathrm{p}}\cdot p)}{\partial y^{\mathrm{2}}}}\text{(11)}& {\displaystyle }& {\displaystyle }{u}_{\mathrm{p}}(N,T,I,x,y)=u_{\mathrm{p}}^{\mathrm{0}}{\mathit{\gamma }}_{\mathrm{n}}(N,x){\mathit{\gamma }}_T(T,y){\mathit{\gamma }}_I\left(I\right)\end{array}$$

Like all biodiversity models, SPEAD must not allow a phenotype to outcompete all other phenotypes in all environments because any such Darwinian demon would drive all its sub-optimal competitors to extinction and trait variance would collapse to zero. In order to make competition for nutrients possible, we have defined two uptake traits so that either low or high values are advantageous in certain environments and disadvantageous in others. The shape of the two trade-offs and the three growth factors are presented in Fig. 2.

Figure 2 Phytoplankton growth factors γ_n (nutrient-dependent), γ_T (temperature-dependent), and γ_I (PAR-dependent). Panels (a) and (b) represent the growth factor as a function of nutrient concentration and temperature, respectively, for different phenotypes. Panels (c) and (d) represent the growth factor as a function of the corresponding trait for different values of the environmental parameter. The maximum of each curve corresponds to the phenotype most adapted to a given environment. Panels (e) and (f) are normalized versions of (c) and (d), respectively, so that their maximum is always 1. Panel (g) is the PAR-dependent growth factor, which is common to all phenotypes in this version of SPEAD.

The first trait allowed to mutate in SPEAD, x, is the logarithm of the half-saturation constant that controls the nutrient limitation factor γ_n(N,x). The half-saturation constant K_n(x) is the DIN concentration at which the nitrogen uptake rate is equal to one-half of the maximum uptake rate for the same temperature and solar irradiance. As the concentrations are always positive and span several orders of magnitude, we use a natural logarithmic scale and define $x=\log \left(\frac{K_{\mathrm{n}}}{K_{\mathrm{n}}^{\mathrm{0}}}\right)$ as our first trait axis, with $K_{\mathrm{n}}^{\mathrm{0}}=\mathrm{1}$ mmol N m⁻³ as a reference value for K_n. The half-saturation constant can be linked to the well-known trade-off between the affinity for a nutrient and the maximum uptake rate, also known as the gleaner–opportunist trade-off (Frederickson and Stephanopoulos, 1981). For a given phenotype, γ_n(N,x) follows a Michaelis–Menten function of nutrient concentration.

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{\mathit{\gamma }}_{\mathrm{n}}(N,x)={\mathit{\gamma }}_{\mathrm{n}}^{\mathrm{\infty }}\left(x\right){\displaystyle \frac{N}{K_{\mathrm{n}}\left(x\right)+N}}={\mathit{\gamma }}_{\mathrm{n}}^{\mathrm{\infty }}\left(x\right){\displaystyle \frac{N}{K_{\mathrm{n}}^{\mathrm{0}}{e}^x+N}}\end{array}$$

The nutrient limitation factor at a high nitrogen concentration, ${\mathit{\gamma }}_{\mathrm{n}}^{\mathrm{\infty }}\left(x\right)$, is set so that each phenotype has an ecological niche wherein it grows faster than its competitors. To obtain the expression of ${\mathit{\gamma }}_{\mathrm{n}}^{\mathrm{\infty }}\left(x\right)$, we first introduce two variables that represent the competitive ability of phytoplankton in two types of environments: the growth rate in the limit of a high nitrogen concentration ($u_{\mathrm{p}}^{\mathrm{\infty }}(T,I,x,y)$; d⁻¹) and the affinity for nitrogen ($f_{\mathrm{p}}(T,I,x,y)$; d⁻¹ mmol⁻¹ m³). These two variables fully define the Michaelis–Menten function for nutrient uptake. The affinity represents the initial slope of the curve at zero resource concentration, while the maximum growth rate represents the horizontal asymptote of the curve that is reached at an infinite resource concentration (see Fig. 2 in Vallina et al., 2019). The affinity is closely related to the growth rate at a low nitrogen concentration, which is equal to $f_{\mathrm{p}}(T,I,x,y)N$. The maximum growth rate is related to the biomass-specific handling rate of nutrient ions at high nitrogen concentrations. The values of $u_{\mathrm{p}}^{\mathrm{\infty }}(T,I,x,y)$ and $f_{\mathrm{p}}(T,I,x,y)$ are given by introducing Eq. (12) into Eq. (11) and taking the limit of high or low concentrations.

$$\begin{array}{ll}\text{(13)}& {\displaystyle }& {\displaystyle }{u}_{\mathrm{p}}^{\mathrm{\infty }}(T,I,x,y)=u_{\mathrm{p}}^{\mathrm{0}}{\mathit{\gamma }}_{\mathrm{n}}^{\mathrm{\infty }}\left(x\right){\mathit{\gamma }}_T(T,y){\mathit{\gamma }}_I\left(I\right){\displaystyle }& {\displaystyle }{f}_{\mathrm{p}}(T,I,x,y)={\displaystyle \frac{u_{\mathrm{p}}^{\mathrm{0}}{\mathit{\gamma }}_{\mathrm{n}}^{\mathrm{\infty }}\left(x\right){\mathit{\gamma }}_T(T,y){\mathit{\gamma }}_I\left(I\right)}{K_{\mathrm{n}}\left(x\right)}}\\ \text{(14)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{u_{\mathrm{p}}^{\mathrm{\infty }}(T,I,x,y)}{K_{\mathrm{n}}\left(x\right)}}\end{array}$$

In order to prevent phenotypes from dominating all nutrient niches (i.e. from being Darwinian demons), we assume that the product of maximum growth and affinity $f_{\mathrm{p}}{u}_{\mathrm{p}}^{\mathrm{\infty }}={\left(u_{\mathrm{p}}^{\mathrm{\infty }}\right)}^{\mathrm{2}}{\left(K_{\mathrm{n}}\right)}^{-\mathrm{1}}$ is constant with (independent of) x (Meyer et al., 2015). Therefore, $u_{\mathrm{p}}^{\mathrm{\infty }}$ is proportional to the square root of K_n and to a factor independent of x, yielding the following expression for ${\mathit{\gamma }}_{\mathrm{n}}^{\mathrm{\infty }}\left(x\right)$:

$$\begin{array}{}\text{(15)}& {\mathit{\gamma }}_{\mathrm{n}}^{\mathrm{\infty }}\left(x\right)=\sqrt{e^x}=e^{x/\mathrm{2}}.\end{array}$$

With this “maximum growth rate–nutrient affinity” trade-off, phenotypes that grow faster than their competitors at high nutrient concentrations (high $u_{\mathrm{p}}^{\mathrm{\infty }}$), called opportunists, are disadvantaged at low concentrations (low f_p) by the same factor, and those which grow comparatively fast at low nutrient concentrations (high f_p), called gleaners, are at a disadvantage under high concentrations (low $u_{\mathrm{p}}^{\mathrm{\infty }}$). This constraint has been observed in bacteria and phytoplankton (Healey and Hendzel, 1980; Button et al., 2004; Elbing et al., 2004; Cermeño et al., 2011; Vallina et al., 2019), has been used in models (Dutkiewicz et al., 2009; Vallina et al., 2014b; Smith et al., 2016; Vallina et al., 2017), and is the simplest way to discriminate the bloom-forming opportunist phytoplankton (such as diatoms) from the ubiquitous gleaners (such as Prochlorococcus) without explicitly representing the complex effects of cell size or the way uptake sites are packed at the cell surface. However, the thermodynamic bases of the trade-off are still unclear (Wirtz, 2002). For any nutrient concentration N, we note that the phenotype corresponding to the largest growth rates is $x=\log \left(\frac{N}{K_{\mathrm{n}}^{\mathrm{0}}}\right)$. This is why, under the assumption of the gleaner–opportunist trade-off defined above, K_n defines the optimal nutrient concentration of each x phenotype at which they are competitively superior (Vallina et al., 2017). This result, however, is dependent on the specific model assumption that $f_{\mathrm{p}}{u}_{\mathrm{p}}^{\mathrm{\infty }}$ is a constant.

The second phytoplankton trait that is allowed to mutate in SPEAD is the optimal temperature. Temperature affects microbes in two ways. One is generic and applies to the whole plankton community. An increase in temperature increases the speed of both primary production and heterotrophic processes for thermodynamic reasons. This effect is often assumed to be exponential. In our model, the exponential factor for autotrophic primary production is α_a=0.056 ^∘C⁻¹, which corresponds to a Q₁₀ of 1.75, slightly lower than the classical value of 1.88 from Eppley (1972) but higher than the values based on the metabolic theory of ecology for photosynthesis (López-Urrutia et al., 2006). The second effect of temperature is phenotype-specific. Each phenotype has an optimal temperature for growth, which is the second trait axis and is denoted by y. The effect of temperature on a given phenotype (x,y) is asymmetric: at temperatures more than 5 ^∘C above y growth ceases, but temperatures below y merely slow growth. We defined our temperature multiplicative growth factor to be as close as possible to the species-specific curves of Eppley (1972):

$$\begin{array}{}\text{(16)}& {\mathit{\gamma }}_T(T,y)=e^{\frac{(T-y)}{\mathrm{\Delta }T}}\left({\displaystyle \frac{y+\mathrm{\Delta }T-T}{\mathrm{\Delta }T}}\right)e^{{\mathit{\alpha }}_{\mathrm{a}}(y-T_{\mathrm{0}})}.\end{array}$$

The temperature tolerance ΔT is set to 5 ^∘C. T₀ is a reference temperature with no ecological meaning. For a fixed value of y, γ_T(T,y) has a maximum at T=y with a value of $e^{{\mathit{\alpha }}_{\mathrm{a}}(y-T_{\mathrm{0}})}$. At $T=y+\mathrm{\Delta }T$ and warmer, growth is impossible. For a given value of the environment temperature T, the phenotypes with the largest growth rates have an optimal temperature y around 2 ^∘C larger than T. This apparent mismatch, whereby the dominant phenotype at temperature T can grow even faster at temperatures a few degrees warmer, is both coherent with other models (Beckmann et al., 2019) and observed in nature (Thomas et al., 2012; Irwin et al., 2012).

In this study, the PAR limitation factor γ_I(I) is the same for all phenotypes. It includes an optimal PAR (I_opt) of 25 W m⁻² and photoinhibition above this level. Our value for I_opt is in the middle of the range considered by Follows et al. (2007), and our expression for γ_I(I) is equivalent to theirs.

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_I\left(I\right)={\mathit{\gamma }}_I^{\mathrm{0}}\left(\mathrm{1}-e^{-\ln (\mathrm{1}+\mathit{\chi })\frac{I}{I_{\mathrm{opt}}}}\right)e^{-\frac{\ln (\mathrm{1}+\mathit{\chi })}{\mathit{\chi }}\frac{I}{I_{\mathrm{opt}}}}\text{(18)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_I^{\mathrm{0}}={\displaystyle \frac{\mathit{\chi }+\mathrm{1}}{\mathit{\chi }}}{e}^{-\frac{\mathrm{1}}{\mathit{\chi }}\ln \left(\frac{\mathrm{1}}{\mathit{\chi }+\mathrm{1}}\right)}\end{array}$$

In the above equation, ${\mathit{\gamma }}_I^{\mathrm{0}}$ is a normalization factor (to ensure that γ_I(I) cannot exceed 1) and χ is an inhibition factor. The higher the inhibition factor, the less photoinhibition there is at irradiances larger than I_opt. In this study, we use χ=12, which is the average value in Follows et al. (2007) for large phytoplankton and corresponds well to published photoinhibition curves (Platt et al., 1980; Whitelam and Codd, 1983; Walsh et al., 2001).

For comparison with data, two additional variables can be estimated from the model: primary production and chlorophyll a concentration. Primary production (PP) is expressed in milligrams of carbon per cubic metre per day (mg C m⁻³ d⁻¹). Our model-based estimate is calculated by multiplying the phytoplankton concentration and the uptake rate, normalizing from nitrogen to carbon with the 106:16 Redfield molar ratio (Redfield, 1934) and then converting from the amount of substance to mass using the molar mass of carbon (12 g mol⁻¹). The chlorophyll a concentration (Chl; mg Chl m⁻³) is obtained by dividing the phytoplankton concentration in mass of carbon by a variable carbon to chlorophyll mass ratio (C : Chl). The C : Chl ratio is estimated as in Vallina et al. (2008) using a function of depth and time developed by Lefèvre et al. (2002), with parameter values calibrated with the observations of Goericke and Welschmeyer (1998). At the surface, C : Chl is a sinusoidal function of the day of year, varying between a maximum of 160 mg C mg Chl⁻¹ at the summer solstice and a minimum of 80 mg C mg Chl⁻¹ at the winter solstice. From the depth at which $I(z,t)=\mathrm{25}$ W m⁻² to the bottom, C : Chl decreases linearly with I(z,t) down to a value of 40 mg C mg Chl⁻¹ when light is absent.

2.2 Aggregate and multi-phenotype models

Traits x and y have an infinity of possible values. In order to solve the equations numerically, the problem needs to be simplified. Two approaches are considered here. In the multi-phenotype or discrete model approach (Fig. 1b), the trait space is discretized, and only a finite number of phenotypes with fixed trait values are simulated. Phenotypes with intermediate trait values are neglected. In the aggregate or continuous model approach (Fig. 1c), the state variables are total phytoplankton concentration, the mean trait values, the trait variances, and the inter-trait covariance. In the continuous trait model, a specific shape of the trait distribution must be assumed a priori (Wirtz and Eckhardt, 1996; Bruggeman and Kooijman, 2007). In the discrete trait model, the trait distribution is an emergent property, and thus it does not need to be assumed beforehand.

In the multi-phenotype model, only n_x values of x and n_y values of y are allowed. The phytoplankton community is divided into n_x×n_y phenotypes. The values of both traits are explicitly bounded by x_min, x_max, y_min, and y_max. Each phenotype is separated from its immediate neighbours by a trait interval $\mathrm{\Delta }x=\frac{x_{max}-x_{min}}{n_x-\mathrm{1}}$ on x or a trait interval $\mathrm{\Delta }y=\frac{y_{max}-y_{min}}{n_y-\mathrm{1}}$ on y. Mutation fluxes at the boundaries (i.e. mutations of the phenotypes with the highest or lowest trait values leading out of the domain) are set to zero. In the interior of our trait domain, the concentration of the phenotype with the jth value of x and the kth value of y, noted P_jk, is controlled by the following equation, where $a_{j,k}=u_{j,k}-g_{j,k}-m_{j,k}$, called the net growth rate, is the sum of all growth and death terms affecting phytoplankton concentration.

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{P}_{j,k}}{\mathrm{d}t}}=& {\displaystyle }{a}_{j,k}(N,T,I)P_{j,k}+{\displaystyle \frac{{\mathit{\nu }}_x}{{\left(\mathrm{\Delta }x\right)}^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }\left(P_{j-\mathrm{1},k}{u}_{j-\mathrm{1},k}+P_{j+\mathrm{1},k}{u}_{j+\mathrm{1},k}-\mathrm{2}{P}_{j,k}{u}_{j,k}\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\mathit{\nu }}_y}{{\left(\mathrm{\Delta }y\right)}^{\mathrm{2}}}}\left(P_{j,k-\mathrm{1}}{u}_{j,k-\mathrm{1}}+P_{j,k+\mathrm{1}}{u}_{j,k+\mathrm{1}}\right.\\ \text{(19)}& {\displaystyle }& {\displaystyle }\left.-\mathrm{2}{P}_{j,k}{u}_{j,k}\right)\end{array}$$

In all our discrete simulations, we impose $x_{min}=-\mathrm{2.5}$ (K_n=0.082 mmol N m⁻³), $x_{max}=+\mathrm{1.5}$ (K_n=4.48 mmol N m⁻³), $y_{min}=\mathrm{18}$ ^∘C, and $y_{max}=\mathrm{30}$ ^∘C. All model values of temperature and DIN concentrations are within these boundaries. We set n_x=25 and n_y=25 in order to ensure that in most cases Δx and Δy are less than the standard deviations of x and y, respectively. Thus, the total number of discrete phenotypes (x,y) is $\mathrm{25}\times \mathrm{25}=\mathrm{625}$.

In the aggregate model, the trait distribution is assumed to be continuous. In this case, and contrary to the multi-phenotype case, the trait axes are formally unbounded, although phenotypes with extreme trait values always have low net growth rates, making them extremely rare. The prognostic variables are six statistical moments of the trait distribution: the total phytoplankton concentration P(t), the mean trait values $\stackrel{\mathrm{‾}}{x}\left(t\right)$ and $\stackrel{\mathrm{‾}}{y}\left(t\right)$, the trait variances V_x(t) and V_y(t), and the inter-trait covariance C_xy(t). They are defined as follows.

$$\begin{array}{ll}\text{(20)}& {\displaystyle }& {\displaystyle }P\left(t\right)=\int \int p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\text{(21)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{x}\left(t\right)={\displaystyle \frac{\mathrm{1}}{P\left(t\right)}}\int \int x\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\text{(22)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{y}\left(t\right)={\displaystyle \frac{\mathrm{1}}{P\left(t\right)}}\int \int y\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\text{(23)}& {\displaystyle }& {\displaystyle }{V}_x\left(t\right)={\displaystyle \frac{\mathrm{1}}{P\left(t\right)}}\int \int {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{\mathrm{2}}p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\text{(24)}& {\displaystyle }& {\displaystyle }{V}_y\left(t\right)={\displaystyle \frac{\mathrm{1}}{P\left(t\right)}}\int \int {\left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)}^{\mathrm{2}}p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y{\displaystyle }& {\displaystyle }{C}_{xy}\left(t\right)={\displaystyle \frac{\mathrm{1}}{P\left(t\right)}}\int \int \left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)\left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)p(x,y,t)\\ \text{(25)}& {\displaystyle }& {\displaystyle }\cdot \mathrm{d}x\mathrm{d}y\end{array}$$

The second-order moments (V_x, V_y and C_xy) are difficult to interpret directly due to their dimensions. In the analyses, we thus transform variances into standard deviations (σ_x and σ_y) and covariance into correlation (R_xy) as follows:

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_x\left(t\right)=\sqrt{V_x\left(t\right)},\text{(27)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_y\left(t\right)=\sqrt{V_y\left(t\right)},\text{(28)}& {\displaystyle }& {\displaystyle }{R}_{xy}\left(t\right)={\displaystyle \frac{C_{xy}\left(t\right)}{{\mathit{\sigma }}_x\left(t\right){\mathit{\sigma }}_y\left(t\right)}}.\end{array}$$

These three diagnostic variables, along with P, $\stackrel{\mathrm{‾}}{x}$, and $\stackrel{\mathrm{‾}}{y}$, are also computed for the multi-phenotype model for comparison. The standard deviations have the same dimensions as the mean traits and can thus be compared to them. Ecologically, they represent trait diversity. Inter-trait correlation is a dimensionless number between −1 and +1, which is easier to interpret than the covariance. A correlation of −1 means above-average values of x always coincide with below-average values of y and vice versa. A correlation of +1 means above-average values of x always coincide with above-average values of y. A correlation of 0 means all combinations are equally possible (i.e. the two traits are independent).

We follow the method developed by Norberg et al. (2001), based on Taylor expansions of the uptake and net growth rates to derive the differential equations for the moments of the trait distribution. We assume a bivariate normal distribution of traits, which is a generalization of the 1D Gaussian function. Normal distributions are observed in nature for the logarithm of size (Cermeño and Figueiras, 2008; Quintana et al., 2008; Downing et al., 2014) and are convenient assumptions for models because they produce the simplest forms for the equations (Wirtz and Eckhardt, 1996; Merico et al., 2009). The derivation is explained in detail in Appendix B. In the absence of trait diffusion, our equations are a particular case of the general equations derived by Bruggeman (2009) for multivariate normal trait distributions. In the single-trait case, they are simpler than the original equations of Merico et al. (2014) and identical to the more recent formulation of Coutinho et al. (2016). In the absence of physical transport, the differential equations followed by the prognostic variables are as follows.

$$\begin{array}{ll}\text{(29)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}=P\left(a(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},t)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_x{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}+C_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}\right)\text{(30)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{x}}{\mathrm{d}t}}=V_x{\displaystyle \frac{\partial a}{\partial x}}+C_{xy}{\displaystyle \frac{\partial a}{\partial y}}\text{(31)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{y}}{\mathrm{d}t}}=V_y{\displaystyle \frac{\partial a}{\partial y}}+C_{xy}{\displaystyle \frac{\partial a}{\partial x}}{\displaystyle }& {\displaystyle \frac{\mathrm{d}{V}_x}{\mathrm{d}t}}=V_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}+\mathrm{2}{V}_x{C}_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}+C_{xy}^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}+\mathrm{2}{\mathit{\nu }}_x\\ \text{(32)}& {\displaystyle }& {\displaystyle }\left(u(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},t)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_x{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial y^{\mathrm{2}}}}+C_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x\partial y}}\right){\displaystyle }& {\displaystyle \frac{\mathrm{d}{V}_y}{\mathrm{d}t}}=V_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}+\mathrm{2}{V}_y{C}_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}+C_{xy}^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}+\mathrm{2}{\mathit{\nu }}_y\\ \text{(33)}& {\displaystyle }& {\displaystyle }\left(u(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},t)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_x{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial y^{\mathrm{2}}}}+C_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x\partial y}}\right){\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{xy}}{\mathrm{d}t}}=V_x{C}_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}+(V_x{V}_y+C_{xy}^{\mathrm{2}}){\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}\\ \text{(34)}& {\displaystyle }& {\displaystyle }+V_y{C}_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}\end{array}$$

The net growth rate $a=u-g-m$ and its derivatives with respect to traits are in all cases taken near the mean trait values ($\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$) and for the values of N, T, and I at time t. The growth rate of the whole phytoplankton community depends first on the net growth rate of the most abundant phenotype (the “winner” of the competition), $a(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},t)$, with correction terms $\left(\frac{\mathrm{1}}{\mathrm{2}}{V}_x\frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}+\frac{\mathrm{1}}{\mathrm{2}}{V}_y\frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}+C_{xy}\frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}\right)$ for the less abundant and generally less fit phenotypes (“losers”). Mean traits increase when larger trait values are associated with larger net growth rates $\left(\frac{\partial a}{\partial x}>\mathrm{0}\mathrm{or}\frac{\partial a}{\partial y}>\mathrm{0}\right)$ and decrease in the opposite case $\left(\frac{\partial a}{\partial x}<\mathrm{0}\mathrm{or}\frac{\partial a}{\partial y}<\mathrm{0}\right)$. The change is faster when trait variances (V_x and V_y) are high. As a consequence, the overall effect of trait diversity on primary production depends on the environmental conditions. In a stable environment, high trait variances diminish the primary production because phenotypes with low growth rates are present. Under frequent disturbances, however, high trait variances increase the short-term adaptive capacity, allowing the community to maintain mean traits close to the optimum and thereby increasing primary production (Smith et al., 2016). We note that when traits covary (C_xy≠0), the change in each trait depends on both nutrient concentration N $\left(V_x\frac{\partial a}{\partial x}\mathrm{for}x\mathrm{and}{C}_{xy}\frac{\partial a}{\partial x}\mathrm{for}y\right)$ and water temperature T $\left(V_y\frac{\partial a}{\partial y}\mathrm{for}\mathrm{y}\mathrm{and}{C}_{xy}\frac{\partial a}{\partial y}\mathrm{for}\mathrm{x}\right)$. Variances decrease due to competition when mean trait values are close to the values that maximize the net growth rate $\left(\frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}<\mathrm{0}\mathrm{and}\frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}<\mathrm{0}\right)$. This is most often the case, since phenotypes that are not optimal tend to be outcompeted. Trait diversity must therefore be maintained by some other process: this is the role of trait diffusion. In these equations, trait diffusion is a (positive) source of variance equal to $\mathrm{2}{\mathit{\nu }}_{x/y}\left(u(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},t)+\frac{\mathrm{1}}{\mathrm{2}}{V}_x\frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}+\frac{\mathrm{1}}{\mathrm{2}}{V}_y\frac{{\partial }^{\mathrm{2}}u}{\partial y^{\mathrm{2}}}+C_{xy}\frac{{\partial }^{\mathrm{2}}u}{\partial x\partial y}\right)$ but does not affect the equations for phytoplankton concentration, mean traits, or covariance. This is coherent with the fact that mutations are symmetrical (no effect on $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$) and neither create nor remove biomass (no effect on P). Trait diffusion does not affect covariance because mutations of the two traits are independent of each other. There is no mechanistic relationship between optimal temperature and half-saturation. Mutations can create all combinations: cold-water gleaners, warm-water gleaners, cold-water opportunists, and warm-water opportunists. However, by increasing variances, trait diffusion decreases the absolute value of correlation. Only the environment can correlate the traits by favouring some combinations over others. Although correlation is defined as a local quantity, for a given depth and time, it is expected to be influenced by the spatio-temporal patterns of environment variations because local communities always contain remnants of past communities and migrants from other locations.

2.3 Physical setting

SPEAD 1.0 has one spatial dimension: the vertical. A depth-resolved simulation is the minimum physical setting in the ocean to resolve the different temperature and nutrient niches and the decisive effect of the vertical mixing on the variances and the covariance. The model is divided into 20 vertical levels from the surface to 200 m deep, with a uniform vertical step of 10 m.

Two processes can transport matter from one vertical level to another and thus need to be added to the differential equations of the 0D model presented in Sect. 2.1 and 2.2. First, PON sinks at a speed of w=1.2 m d⁻¹. Sinking is represented by an extra term in the PON equation (Eq. 4) equal to $-w\frac{\partial D}{\partial z}$, where z is depth and is always positive. The vertical derivative $\frac{\partial D}{\partial z}$ is solved by a first-order upwind scheme. Second, tracers are vertically mixed by turbulent diffusion. Vertical turbulent diffusion (called vertical diffusion from now on, which is unrelated to trait diffusion) tends to homogenize the spatial distribution of each tracer. It is controlled by the vertical diffusivity parameter κ_z, expressed in square metres per day (m² d⁻¹). The vertical diffusion of a conserved tracer A is $\frac{\partial }{\partial z}\left({\mathit{\kappa }}_{\mathrm{z}}(z,t)\frac{\partial A(z,t)}{\partial z}\right)$. This expression applies to all concentrations of the discrete model and to concentrations N, P, Z, and D of the continuous trait model. As a consequence, their full local (Eulerian) time derivatives at a given depth z are as follows.

$$\begin{array}{}\text{(35)}& {\displaystyle }& {\displaystyle \frac{\partial Z}{\partial t}}={\displaystyle \frac{\mathrm{d}Z}{\mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{z}}{\displaystyle \frac{\partial Z}{\partial z}}\right)\text{(36)}& {\displaystyle }& {\displaystyle \frac{\partial N}{\partial t}}={\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{z}}{\displaystyle \frac{\partial N}{\partial z}}\right)\text{(37)}& {\displaystyle }& {\displaystyle \frac{\partial D}{\partial t}}={\displaystyle \frac{\mathrm{d}D}{\mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{z}}{\displaystyle \frac{\partial D}{\partial z}}\right)-w{\displaystyle \frac{\partial D}{\partial z}}\text{(38)}& {\displaystyle }& {\displaystyle \frac{\partial P}{\partial t}}={\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{z}}{\displaystyle \frac{\partial P}{\partial z}}\right)\end{array}$$

In contrast, the vertical turbulent diffusions of the statistical moments $\stackrel{\mathrm{‾}}{x}$, $\stackrel{\mathrm{‾}}{y}$, V_x, V_y, and C_xy do not follow a simple analytical expression because they are not material quantities and thus are not conserved during mixing. For instance, the mixing of two phytoplankton communities with different mean traits creates additional trait variance. Fortunately, vertical mixing conserves the sums of phytoplankton trait values $\left(\mathrm{equal}\mathrm{to}P\stackrel{\mathrm{‾}}{x}\mathrm{and}P\stackrel{\mathrm{‾}}{y}\right)$ and the sum of squared trait values $\left(P(V_x+{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}),P(V_y+{\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}),\mathrm{and}P(C_{xy}+\stackrel{\mathrm{‾}}{x}\stackrel{\mathrm{‾}}{y})\right)$. Following the method of Bruggeman (2009), we therefore use these conserved moments as tracers involved in physical transport. Their full Eulerian derivatives are computed as follows using the product rule.

$$\begin{array}{ll}\text{(39)}& {\displaystyle }& {\displaystyle \frac{\partial \left(P\stackrel{\mathrm{‾}}{x}\right)}{\partial t}}=P{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{x}}{\mathrm{d}t}}+\stackrel{\mathrm{‾}}{x}{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{z}}{\displaystyle \frac{\partial \left(P\stackrel{\mathrm{‾}}{x}\right)}{\partial z}}\right)\text{(40)}& {\displaystyle }& {\displaystyle \frac{\partial \left(P\stackrel{\mathrm{‾}}{y}\right)}{\partial t}}=P{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{y}}{\mathrm{d}t}}+\stackrel{\mathrm{‾}}{y}{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{z}}{\displaystyle \frac{\partial \left(P\stackrel{\mathrm{‾}}{y}\right)}{\partial z}}\right){\displaystyle }& {\displaystyle \frac{\partial \left(P\left(V_x+{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}\right)\right)}{\partial t}}=P{\displaystyle \frac{\mathrm{d}{V}_x}{\mathrm{d}t}}+\mathrm{2}P\stackrel{\mathrm{‾}}{x}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{x}}{\mathrm{d}t}}+\left(V_x+{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}\right)\\ \text{(41)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{z}}{\displaystyle \frac{\partial \left(P\left(V_x+{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}\right)\right)}{\partial z}}\right){\displaystyle }& {\displaystyle \frac{\partial \left(P\left(V_y+{\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}\right)\right)}{\partial t}}=P{\displaystyle \frac{\mathrm{d}{V}_y}{\mathrm{d}t}}+\mathrm{2}P\stackrel{\mathrm{‾}}{y}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{y}}{\mathrm{d}t}}+\left(V_y+{\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}\right)\\ \text{(42)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{z}}{\displaystyle \frac{\partial \left(P\left(V_y+{\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}\right)\right)}{\partial z}}\right){\displaystyle }& {\displaystyle \frac{\partial \left(P\left(C_{xy}+\stackrel{\mathrm{‾}}{x}\stackrel{\mathrm{‾}}{y}\right)\right)}{\partial t}}=P{\displaystyle \frac{\mathrm{d}{C}_{xy}}{\mathrm{d}t}}+P\stackrel{\mathrm{‾}}{y}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{x}}{\mathrm{d}t}}+P\stackrel{\mathrm{‾}}{x}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{y}}{\mathrm{d}t}}\\ \text{(43)}& {\displaystyle }& {\displaystyle }+\left(C_{xy}+\stackrel{\mathrm{‾}}{x}\stackrel{\mathrm{‾}}{y}\right){\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}+{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{z}}{\displaystyle \frac{\partial \left(P\left(C_{xy}+\stackrel{\mathrm{‾}}{x}\stackrel{\mathrm{‾}}{y}\right)\right)}{\partial z}}\right)\end{array}$$

The above Eulerian derivatives are the derivatives used in SPEAD to proceed from one time step to the next. Then, at each time step, we recover the necessary statistical moments $\left(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y},V_x,V_y,\mathrm{and}{C}_{xy}\right)$ by back-computing them.

$$\begin{array}{}\text{(44)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{x}={\displaystyle \frac{\left(P\stackrel{\mathrm{‾}}{x}\right)}{P}}\text{(45)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{y}={\displaystyle \frac{\left(P\stackrel{\mathrm{‾}}{y}\right)}{P}}\text{(46)}& {\displaystyle }& {\displaystyle }{V}_x={\displaystyle \frac{\left(P\left(V_x+{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}\right)\right)}{P}}-{\displaystyle \frac{{\left(P\stackrel{\mathrm{‾}}{x}\right)}^{\mathrm{2}}}{P^{\mathrm{2}}}}\text{(47)}& {\displaystyle }& {\displaystyle }{V}_y={\displaystyle \frac{\left(P\left(V_y+{\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}\right)\right)}{P}}-{\displaystyle \frac{{\left(P\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}{P^{\mathrm{2}}}}\text{(48)}& {\displaystyle }& {\displaystyle }{C}_{xy}={\displaystyle \frac{\left(P\left(C_{xy}+\stackrel{\mathrm{‾}}{x}\stackrel{\mathrm{‾}}{y}\right)\right)}{P}}-{\displaystyle \frac{\left(P\stackrel{\mathrm{‾}}{x}\right)}{P}}{\displaystyle \frac{\left(P\stackrel{\mathrm{‾}}{y}\right)}{P}}\end{array}$$

The derivatives with respect to z used in the vertical diffusion terms are estimated by an implicit scheme in order to avoid numerical instability. The depth-resolved model is solved in time with a fourth-order Runge–Kutta numerical scheme.

Three environmental forcings are necessary to run the model: temperature, PAR, and vertical diffusivity. All three depend on depth and time and have been set to values from the Sargasso Sea. The forcings are seasonal. Inter-annual variations and the day–night cycle are neglected. Temperature and surface PAR (I₀(t)) directly affect the rates of plankton growth and death. They are set for each day using observations collected during the Bermuda Atlantic Time-series Study (BATS) (Steinberg et al., 2001). PAR availability is assumed to decrease exponentially with depth $\left(I(z,t)=I_{\mathrm{0}}\left(t\right)e^{-k_{w}z}\right)$, with a PAR vertical attenuation coefficient (k_w) of 0.04 m⁻¹. Self-shading by phytoplankton is neglected. The vertical diffusivity κ_z is the third forcing. Contrary to temperature and PAR, it has not been directly observed. Therefore, the turbulent diffusivity comes from the physical model GOTM for the Sargasso Sea (Bruggeman and Bolding, 2014). All three forcings as functions of time and depth are shown in Fig. 3.

Figure 3 Distribution in depth and time of three environmental variables: (a) temperature, (b, d) photosynthetically active radiation (PAR), and (c) vertical turbulent diffusivity. The black curve in (d) represents the lower limit of the mixed layer.

2.4 List of simulations

The simulation of the aggregate two-trait model with mutation rates ν_x=0.001 and ν_y=0.01 ^∘C² is our standard simulation for this study. ν_x is expressed without a unit because the trait axis x is in logarithmic scale, but like ν_y it is a variance increase per generation. Most of the results presented in Figs. 4, 5, 6, 7, and 9 come from this standard simulation. The bulk properties of SPEAD 1.0 (total primary production, total phytoplankton biomass, nutrient and detritus concentrations) are validated using observations from the BATS station in the Sargasso Sea (Steinberg et al., 2001; Vallina et al., 2008; Vallina, 2008). The multi-phenotype discrete version of SPEAD is used to validate (i) the assumption made in the aggregate continuous model that traits are normally distributed and (ii) the simulated values and tendencies of the moments of the continuous trait distribution. In order to better understand the behaviour of the model, the standard simulation is also compared to simulations with different mutation rates and to simulations with adaptive dynamics for only one trait, keeping the other trait unable to mutate but at its optimal value.

Trait diffusion is a relatively recent concept, and the values of the mutation rates are not yet well calibrated by observations. To obtain a qualitative idea of the ecosystem model behaviour, we tried a wide range of values for ν_x (from 0.00001 to 0.1). The largest value was chosen for its similarity to the trait diffusivity parameter used by Merico et al. (2014) and Chen and Smith (2018) to account for the observed trait variance. However, ν_x=0.1 allows the phytoplankton to reach a variance of V_x=4 in only 20 generations, since 2ν_x is added to the phytoplankton population variance at each generation. This variance is the maximum allowed in the discrete model and corresponds to having half the community at each extreme of the trait axis ($x=-\mathrm{2.5}$ and $x=+\mathrm{1.5}$). However, laboratory experiments based on single clones show significant evolution only on timescales of hundreds to thousands of generations (Schlüter et al., 2016). For this reason we also conducted simulations with mutation rates as low as 0.00001 and a control simulation without trait diffusion at all (ν_x=0). As the mutation rate has the dimension of trait squared and as the range of temperature is around 3 times larger than the range of nutrient concentration logarithms, we fixed the same ratio of mutation rates, $\frac{{\mathit{\nu }}_y}{{\mathit{\nu }}_x}=\mathrm{10}$ ^∘C², for all simulations. We checked that departing from this ratio did not qualitatively affect our results. In total, we conducted simulations for 10 sets of mutation rates, including the control case.

A two-trait model is not simply the superposition of two one-trait models for at least two reasons. First, when two environmental factors limit biomass growth but only one is included in the model, the simulation is likely to overestimate the phytoplankton growth rate. Second, when there is a strong inter-trait correlation, each environmental factor impacts both traits. For instance, if the ambient DIN concentration (N) is below the (geometric) mean half-saturation constant ($e^{\stackrel{\mathrm{‾}}{x}}$), the competition for nutrients will select for phenotypes with a lower half-saturation constant. If at the same time the half-saturation is negatively correlated with optimal temperature (i.e. if phenotypes with low half-saturation constants also tend to have high optimal temperatures), the competition for nutrients will also increase the concentrations of phenotypes with high optimal temperatures, in addition to the effect of environment temperature. In the conceptual model of Savage et al. (2007), the inter-trait correlation in a two-trait model led to higher variances and to a considerable improvement in the ability of the mean phytoplankton traits to track optimal values controlled by environmental conditions compared with one-trait models. In order to know whether these results also apply to our model, we compare the dynamics of traits x and y in SPEAD to the dynamics of simplified one-trait models wherein either x or y varies between phenotypes and the other trait is optimized instantaneously (i.e. set to the optimal value at each location and time given the environmental conditions).

The time step for our simulations is 6 h. At the first time step and at all vertical levels, the DIN concentration is initialized to 1.8 mmol N m⁻³, phytoplankton and zooplankton concentrations to 0.1 mmol N m⁻³, and the PON concentration to 0.0 mmol N m⁻³. The total amount of nitrogen in the water column is conserved, and every loss below 200 m due to PON sinking is compensated for by an equivalent gain of DIN, also at 200 m. Mean logarithm of half-saturation and mean optimal temperature are initialized at −0.5 (corresponding to K_n=0.61 mmol N m⁻³) and 24 ^∘C, respectively, with initial standard deviations of 0.1 and 0.3 ^∘C. Each simulation is run for at least 3 years and until convergence is reached. Our convergence criterion is that, for every day of year and every depth level, the difference between the last 2 years should be less than 0.1 % for P, V_x, and V_y, less than 0.1 % of the modelled range for $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$, and less than 0.001 for R_xy. In other words, convergence is achieved when the seasonal cycle of the model state variables is repeated from year to year. The results shown are in all cases from the last simulated year. We checked that total nitrogen was the only feature in the initial conditions that affected the results.

3 Results

3.1 Bulk modelled properties and comparison with observations

The first step to validate the SPEAD model is to compare some bulk properties with observations from the Sargasso Sea in order to assess the realism of the trait-independent biogeochemical parameters (Table 2). In Fig. 4, the primary production, chlorophyll, and DIN and PON concentrations of the aggregate model are compared to a 10-year climatology of monthly observations of primary production, chlorophyll, nitrate, and PON concentrations (Vallina et al., 2008). The distribution pattern of bulk properties in SPEAD does not depend much on the mutation rates and is very similar to those of 0-trait simulations without trait diversity (see the Supplement). Differences between these simulations are small compared with the differences between the model and observations. With carefully chosen values for the parameters of the ecosystem model, the model and observations agree well, although some minor discrepancies exist due to the simplicity of our parameterizations.

Figure 4 Distribution in depth and time of (a) model primary production, (c) chlorophyll concentration, (e) dissolved inorganic nitrogen concentration, and (g) particulate organic nitrogen concentration for ν_x=0.001 and ν_y=0.01 ^∘C². Each variable is compared with equivalent observations in the Sargasso Sea (b, d, f, h).

Primary production is the state variable best reproduced by the model, with a maximum around 10 mg C m⁻³ d⁻¹ in the first 50 m in February and March in both the model and the observations. From May to September, primary production spreads slightly more in depth but is overall around half its maximum value. Primary production is negligible deeper than 80 m. The chlorophyll concentration is reproduced with the right order of magnitude, an absolute maximum correctly located in February between 50 and 80 m deep at about 0.25 mg Chl m⁻³, and a deep chlorophyll maximum around 100 m in summer. However, the model chlorophyll concentration is lower than the observations in spring and higher in late autumn. The chlorophyll concentration begins to increase in December in the model but only in February in the observations. Both the chlorophyll concentration and primary production are proportional to the phytoplankton concentration. The reason for the temporal mismatch between SPEAD and the observations for the chlorophyll concentration, but not in primary production, must then be related to the temporal variability of the other factors affecting these two quantities: the nitrogen uptake rate and the C : Chl ratio. The relatively high primary production and very low chlorophyll concentration observed in December might be better accounted for if the uptake rate were faster in December than in February, despite the lower availability of nutrients and light, or if turbulence were included in the estimation of the C : Chl ratio so that it reaches its lowest value in February, when the waters are best mixed and phytoplankton cannot stay close to the surface (Taylor et al., 1997), rather than in December, when the surface light intensity is minimal (Lefèvre et al., 2002; Jakobsen and Markager, 2016).

Dissolved inorganic nitrogen is compared to the observed nitrate concentration, knowing that this form of DIN dominates at depth but co-occurs with nitrite and ammonium, which are also components of DIN. The modelled DIN and the observed nitrate concentrations share the same range, with a maximum of 2.8 mmol m⁻³ in the observations and 3.3 mmol m⁻³ in the model. Another common point between the model and the observations is that both reach a maximum at the bottom of our setting at 200 m, with concentrations between 2.1 and 2.8 mmol m⁻³ during most of the year. Because of the strong vertical mixing, from February to April, the concentrations are lower, between 1.5 and 2.0 mmol m⁻³, but still maximal. However, from June to January, the modelled DIN concentration exhibits a second maximum that is absent from the observations. This second maximum is located just below the euphotic layer between 100 m and 150 m deep, with values as high as 3.3 mmol m⁻³ in late November. Both modelled and observed concentrations are minimal at the surface due to the nitrogen uptake by phytoplankton. However, their values diverge by more than 1 order of magnitude. Modelled DIN concentrations at the surface vary between 0.18 and 1.34 mmol m⁻³, with a mean of 0.67 mmol m⁻³, whereas observed nitrate concentrations vary between 0 and 0.11 mmol m⁻³, with a mean of 0.03 mmol m⁻³. We assume that these discrepancies are due to the contribution of ammonium and possibly nitrite, since the few studies reporting measured concentrations of ammonium in the Sargasso Sea (Menzel and Spaeth, 1962; Brzezinski, 1988) have shown that ammonium was more homogeneously distributed in the upper 200 m than nitrate and was the dominant form of dissolved nitrogen from the surface to 100 m deep.

Particulate organic nitrogen distributions from the model and observations are relatively similar, with a maximum around 0.5 mmol m⁻³ (0.45 mmol m⁻³ in the observations, 0.51 mmol m⁻³ in the model) in April and May at depths between 30 m and 80 m. However, the seasonality and vertical gradients are much larger in the model, wherein particles are very rare in autumn and nearly absent at depths greater than 100 m, whereas observed PON concentrations are never below 0.08 mmol m⁻³. The observations might be better explained if a minority of particle production went to a slowly remineralizing refractory pool, enabling them to stay during the whole year and to reach greater depths (Aumont et al., 2017), but we did not increase the complexity of the particle parameterization because this is not the focus of our study.

3.2 Trait distribution of SPEAD and comparison with a multi-phenotype model

The second step to validate the aggregate SPEAD model and the only validation of its bivariate trait distribution is comparing it to the multi-phenotype model. Although both are models and thus simplify reality in similar ways, the multi-phenotype model is used as a reference for two reasons. First, it is more intuitive than the aggregate model, with birth and death processes and mutations to the nearest neighbours as the only terms in the equations. Therefore, the moments of the trait distribution in the multi-phenotype model can be used as a control to confirm that the equations of the aggregate model are correct. Second, the multi-phenotype model does not assume any particular trait distribution shape and can be used to validate the a priori assumption of the aggregate model that the trait distribution is a bivariate normal distribution.

The spatial and temporal patterns of the phytoplankton concentration, mean traits, trait standard deviations, and inter-trait correlation for the standard simulation with mutations rates of ν_x=0.001 and ν_y=0.01 ^∘C² are shown in Fig. 5. The value of $\stackrel{\mathrm{‾}}{x}$ varies between −0.83 (K_n=0.44 mmol N m⁻³) and +0.6 (K_n=1.82 mmol N m⁻³), with standard deviations between 0.31 and 0.77. The value of $\stackrel{\mathrm{‾}}{y}$ varies between 22.0 and 26.1 ^∘C, with standard deviations between 0.81 and 1.92 ^∘C. By comparison, the modelled DIN concentration varies between 0.18 and 3.31 mmol N m⁻³, and the water temperature varies between 18.5 and 27.8 ^∘C. As expected, the mean trait values remain consistently within the range of the environmental drivers to which they adapt. Because the best competitor at a given time and depth needs tens of generations to become dominant after having been a rare phenotype, the mean traits react with a delay of 1 to 2 months and with a lower amplitude than their drivers. Cold-water opportunists (high x trait, low y trait) dominate in winter and spring throughout the water column. In summer, they are slowly replaced by warm-water gleaners (low x trait, high y trait) in the upper 70 m but retain dominance at greater depths at which their half-saturation constants continue increasing and their optimal temperatures continue decreasing. The coexistence of two very distinct communities in summer and early autumn is made possible by the intense stratification, which creates a physical barrier between the different depth levels. In late autumn, the two communities are rapidly mixed by vertical turbulent diffusion, producing a peak in the standard deviation of each trait: in other words, a peak in the local (alpha) diversity. Then, as the water column becomes more homogeneous, competition selects for a single dominant phenotype, reducing the trait diversity until the next autumn. Inter-trait correlation is negative at all times and depths due to the negative correlation of the environmental drivers. High DIN concentrations generally coincide with low temperatures, favouring cold-water opportunists. This happens during winter because turbulent mixing brings nutrient-rich cold waters from the deep layers up to the surface. During summer, the consumption of nutrients by primary producers leads to a coincidence of warm temperatures with low DIN concentrations at the surface. In late autumn, the negative correlation reaches its maximum absolute value when the two main communities are suddenly mixed. During the rest of the year, trait diffusion progressively reduces the inter-trait correlation.

Figure 5 Distribution in depth and time of trait distribution moments for ν_x=0.001 and ν_y=0.01 ^∘C²: (a) phytoplankton concentration, (b) (geometric) mean half-saturation, (c) optimal temperature, (d) inter-trait correlation, (e) half-saturation logarithm standard deviation, and (f) optimal temperature standard deviation. For readability, the mean value of trait x is transformed into a nitrogen concentration in (b). To speak properly and in contrast to other means present in this study, the mean half-saturation is a geometric mean, not an arithmetic mean.

In Fig. 6, the state variables of the aggregate model are compared to the trait distribution moments of the discrete model. The discrete model is considered as a “truth” and the difference between the two models as an “error”, which is positive if the value is higher in the aggregate model. The aggregate model reproduces P and $\stackrel{\mathrm{‾}}{x}$ very precisely, with linear determination coefficients (R²) of 0.998 and 0.988, respectively. The biases (mean errors) are very low: +0.0005 mmol N m⁻³ on P and −0.04 on $\stackrel{\mathrm{‾}}{x}$. The bias of $\stackrel{\mathrm{‾}}{y}$ is larger at −0.50 ^∘C, but the coefficient of determination is still very high at R²=0.862. The error is largest in the deep community in summer and early autumn, reaching a maximum of −1.51 ^∘C in early September around 100 m. The most likely reason why the aggregate model underestimates $\stackrel{\mathrm{‾}}{y}$, but not $\stackrel{\mathrm{‾}}{x}$, is that the response of phytoplankton to temperature is asymmetrical. Increases in the environment temperature put more selective pressure on the phytoplankton community than decreases. This feature is poorly taken into account by the aggregate SPEAD model because of its assumption that traits are normally distributed. There is a mismatch between the asymmetrical shape of temperature niches and the imposed symmetrical shape of the distribution of optimal temperatures (y trait) in the aggregate model. This mismatch does not happen for the x trait. As is typically the case in aggregate models, there are more errors in the higher-order moments, in our case the standard deviations and the inter-trait correlation. The coefficients of determination for σ_x, σ_y, and R_xy are R²=0.813, R²=0.462, and R²=0.896, respectively. Their biases are −0.04, +0.09 ^∘C, and −0.001, respectively, which is around 10 % of the mean value for σ_y and σ_x and negligible for R_xy. All three variables decrease much faster in early winter in the aggregate model than in the multi-phenotype model. Additionally, in summer, there is a strong discrepancy for σ_y. In the discrete model, σ_y can reach as low as 0.57 ^∘C, but in the aggregate model it is never less than 0.81 ^∘C and very rarely less than 1.0 ^∘C.

Figure 6 Comparison at all depths and times of the aggregate and multi-phenotype model state variables for ν_x=0.001 and ν_y=0.01 ^∘C²: (a) total phytoplankton concentration, (b) (geometric) mean half-saturation, (c) mean optimal temperature, (d) inter-trait correlation, (e) half-saturation logarithm standard deviation, and (f) optimal temperature standard deviation. Blue is winter, green is spring, red is summer, and black is autumn.

The main errors on σ_x, σ_y, and R_xy are caused by the aggregate model's assumption of a multivariate normal distribution, which is not strictly correct based on the results of the discretely resolved model. In Fig. 7, we show in a 2D colour plot how the two traits are distributed in the discrete model at three different depths (surface, 50 m, and 100 m) at the end of each season (March, June, September, December). This distribution is compared with the bivariate normal distribution of the aggregate model, represented by ellipses. In March, when the waters are well mixed, and in June, when the stratification has just begun, the traits are normally distributed and the two models agree. There is only a small error on the distribution of optimal temperature. In March at all depths and in June at 100 m, the normal distribution of the aggregate model contains more phenotypes with low optimal temperatures than the distribution of the discrete model. In summer, the traits are also normally distributed near the surface. However, the distribution of optimal temperature is markedly right-skewed deeper in the water column. Optimal temperatures below that of the most common phenotypes are extremely rare, whereas those larger than this level are more common. The aggregate model has its $\stackrel{\mathrm{‾}}{y}$ below the peak of the multi-phenotype model and has a much larger σ_y. This is the largest error for both a mean trait value and a trait standard deviation in this study. Since nothing similar occurs with half-saturation, this error must be linked to the right skew in the temperature-dependent growth factor when expressed as a function of optimal temperature (Fig. 2d). In stable environments with little change in temperature with time and little vertical mixing, the distribution of optimal temperature tends to become naturally right-skewed. However, our results show that re-mixing (in late autumn), fast environmental change (near the surface), and trait diffusion can reduce or eliminate this skew so that the trait distribution is often close to normality. In December during the re-mixing phase, the trait distribution completely deviates from normality and becomes bimodal, with a community of warm-water gleaners and a community of cold-water opportunists co-occurring throughout the water column. At this time the standard deviations and the inter-trait correlation are at their annual maxima. The moments of the trait distribution at that time are very well captured by the aggregate model. However, assuming a normal trait distribution is not only wrong in terms of ecological description but also leads to incorrect dynamics during winter. In winter, the ecological selection in the now mixed waters reduces trait diversity, and trait diffusion reduces the inter-trait correlation. These processes occur faster in the aggregate model than in the multi-phenotype model (see Fig. 6) because the selective pressure is larger for a normal distribution than for a bimodal one. From a mathematical point of view, this can be shown in a simplified one-trait model. Selection through competition reduces the trait variance at a speed equal to $\frac{\mathrm{1}}{\mathrm{2}}(M_{\mathrm{4}}-V_x^{\mathrm{2}})\frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}$ (this is a one-trait unskewed version of Eq. B7), where M₄ is the fourth-order moment or “kurtosis”. In a Gaussian distribution, $M_{\mathrm{4}}=\mathrm{3}{V}_x^{\mathrm{2}}$. Bimodal distributions have a lower kurtosis; therefore, they are affected more slowly by ecological selection. By design, the aggregate model cannot account for this effect because it assumes a unimodal Gaussian distribution. From a more ecological point of view, it can be noted that in order to replace a bimodal distribution by a unimodal one with a smaller variance, a previously rare intermediate phenotype must rise to prominence and previously dominant phenotypes must become rare, which is a dramatic change. By comparison, in an already unimodal Gaussian distribution, reducing the variance only means making rare and extreme phenotypes even rarer.

Figure 7 Concentrations of each phenotype of the multi-phenotype model (colour) compared with lines of equal density of the aggregate model (ν_x=0.001 and ν_y=0.01 ^∘C²). Subplots correspond to days 71, 161, 251, and 341 and depths of 0, 50, and 100 m. Dashed lines indicate the optimal competitor.

3.3 Trait dynamics with different mutation rates

In this section, we compare the results of simulations conducted with nine different sets of mutation rates, from ν_x=0.00001 and ν_y=0.0001 ^∘C² to ν_x=0.1 and ν_y=1.0 ^∘C². A control simulation with no trait diffusion is also included. That amounts to a total of 10 simulations. This comparison highlights the unique role played by trait diffusion in SPEAD, even at very low mutation rates. The ratio of mutation rates, $\frac{{\mathit{\nu }}_y}{{\mathit{\nu }}_x}=$ 10 ^∘C², is the same in all simulations presented in this study. We do not assume that this ratio always has this value in nature. Simulations with other values were also performed, and the ratio was found to have little effect on the results. The mean value and variance of each trait are very similar in all simulations in which they have the same mutation rate, whereas primary production and correlation are impacted in similar proportions by both mutation rates. Therefore, the variability in $\frac{{\mathit{\nu }}_y}{{\mathit{\nu }}_x}$ was overlooked for simplicity.

For each simulation, Fig. 8 shows the values of depth-integrated primary production per year and the yearly averaged values and ranges of $\stackrel{\mathrm{‾}}{x}$, $\stackrel{\mathrm{‾}}{y}$, σ_x, σ_y, and R_xy. Additional diagnostics are presented in Table 3. The number of years to converge to a steady state, beyond being just a numerical issue, can also serve as an ecological indicator of the time needed to damp a perturbation or to adapt to a new physical setting, although, by design, our convergence times cannot be less than 3 years. For either trait, the maximum value of $\frac{V}{\mathrm{2}\mathit{\nu }}$ is the number of generations required to reach the highest trait variances of the simulation in the absence of ecological selection and with trait diffusion as the only source of variance. Although highly idealized, this number is a proxy for the timescale of evolutionary processes. Table 3 also assesses whether bimodality is present in the discrete model and whether the mean traits come within 1 standard deviation of the discrete model boundaries at some point in the year. Bimodality was assessed visually based on the trait distributions of the discrete model during the December mixing event.

Figure 8 Primary production and trait distribution moments for different mutation rates. The ratio $\frac{{\mathit{\nu }}_y}{{\mathit{\nu }}_x}$ is kept constant and equal to 10 ^∘C². Moment ranges are represented by error bars and their mean by dots. The mean traits are compared to the extreme values of their environmental drivers (dissolved inorganic nitrogen concentration and temperature). The trait standard deviations are compared to their values in a uniform distribution within the boundaries of the discrete model.

Table 3 Convergence time and various properties of SPEAD 1.0 simulations with different mutation rates.

Primary production lies between 146.8 and 147.0 g C m⁻² yr⁻¹ for all mutation rates between ν_x=0 (control) and ν_x=0.001, then decreases at higher trait diffusivities to finally reach 130.8 g C m⁻² yr⁻¹ for ν_x=0.1. This result agrees with the model of Chen et al. (2019) as applied to the North Pacific. Their phytoplankton community was characterized by one trait, cell size, which is somehow related to our half-saturation trait x. They found that primary production was diminished when ν_x increased, but they only considered relatively high mutation rates between 0.01 and 0.1, as well as a control simulation. Under relatively stable conditions, we find that fast mutation rates (ν_x>0.01) are a drawback for primary production because they promote large trait variances, allowing non-competitive phenotypes (i.e. under-performers) to proliferate. However, phytoplankton mutating very fast could be invaded by phenotypes mutating more slowly. Therefore, we do not expect them to be common in nature.

In the simulations with ν_x=0.01, ν_x=0.03, and ν_x=0.1, the mean trait values remain close to their environmental drivers, and their range over the year is as wide as that of the DIN concentration and temperature, respectively. On average, the community adapts nearly instantaneously to its environment. However, the cost for this apparent success in fast-tracking the environmental conditions is that the standard deviations of the trait distribution are very high and close to that of a uniform distribution between the trait boundaries of the discrete model. Given that the trait domain of the discrete model is already wider than the ranges of DIN concentration and in situ temperature, this result suggests that either (1) phenotypes that are maladapted at all depths and throughout the year are common (explaining the low primary production) or (2) that we have reached the limit of validity of the aggregate approach. In these simulations, skewed or bimodal distributions are extremely rare and smoothed out in a few days, and correlations are negligible, even in December when the stratification is broken, because trait diffusion is a symmetrical process that constantly replenishes all rare phenotypes, including warm-water opportunists and cold-water gleaners that are maladapted at all depths and during the entire year. The simulations with large mutation rates converge in 3 or 4 years and can sustain their variances in fewer than 60 generations, which is in less than a year. They use mutations to follow the seasonal cycle of their environment faster than the usual timescales of evolution, even for phytoplankton (Schlüter et al., 2016).

When the mutation rates are lower, the mean traits still vary during the year but not as much or as fast as the physical environment, and no phenotypes are found outside the trait domain of the discrete model. With ν_x at 0.001 or lower, several years are required to sustain the variance and to converge to a seasonally stable state. In this case, the mutations create variance over the long term, facilitating the ecological successions of phenotypes seasonally and adaptive evolution inter-annually. However, low mutation rates do not allow the community to evolve seasonally. Bimodality is present in the discrete version, at least during the late autumn mixing, and lasts longer as the mutation rates decrease. The variances increase when the trait diffusivity parameters increase, which is what trait diffusion was designed for. We note that, contrary to chemostat models (Merico et al., 2009), SPEAD 1.0 does not require trait diffusion to sustain a positive trait diversity: the trait variances do not collapse to zero even in the absence of trait diffusion. The late autumn mixing is a source of variance in its own right, avoiding the collapse of trait diversity even in the control simulation. However, the trait standard deviations in the control case are very low, between 0.13 and 0.43 for x and between 0.49 and 1.64 ^∘C for y. Trait diversity appears even lower when accounting for the fact that correlation between x and y is blocked at −1. The x and the y traits totally determine each other, as if there were only one trait and no extra degree of freedom. The only active phenotypes are located on a straight line. Trait diffusion is not necessary to sustain variance, but it is necessary to allow the model to explore the entire trait space and to adapt to entirely new sets of environmental conditions.

Increasing trait diffusivity to ν_x=0.00001, ν_x=0.00003, or ν_x=0.0001 does not produce a large increase in trait variance, which keeps being extremely low and mainly controlled by the December mixing, but it reduces the inter-trait correlation to moderate values, mainly between −0.75 and −0.5. Simulations with the above mutation rates also share high primary production and timescales of a few years (approximately hundreds to thousands of generations) to adapt to their environments, which are features coherent with our expectations of the effect of inter-generational mutations.

3.4 Trait dynamics compared with one-trait models

In Figs. 9 and 10, the trait distribution moments of SPEAD at the surface are compared with the environmental drivers (DIN concentration and water temperature) and with the outputs of two single-trait aggregate models, wherein only the half-saturation constant or only optimal temperature is allowed to vary between phenotypes, subject to trait diffusion. Figure 9 shows the comparison for the standard values of the mutation rates: ν_x=0.001 and ν_y=0.01 ^∘C². However, the differences between two-trait and one-trait models are likely to be larger when correlations between x and y are large. Therefore, in Fig. 10, we compare the two-trait and one-trait model distributions obtained with the lowest non-zero mutation rates of ν_x=0.00001 and ν_y=0.0001 ^∘C², which lead to inter-trait correlations between −0.8 and −1 in the two-trait model.

Figure 9 SPEAD state variables at the surface for ν_x = 0.001 and ν_y=0.01 ^∘C² (standard) compared with the state variables of one-trait models: (a) total phytoplankton concentration, (b) (geometric) mean half-saturation, (c) mean optimal temperature, (d) inter-trait correlation, (e) half-saturation logarithm standard deviation, and (f) optimal temperature standard deviation. The dashed lines represent the environmental drivers.

Figure 10 SPEAD state variables at the surface for ν_x=0.00001 and ν_y=0.0001 ^∘C² (low mutation rates) compared with the state variables of one-trait models: (a) total phytoplankton concentration, (b) (geometric) mean half-saturation, (c) mean optimal temperature, (d) inter-trait correlation, (e) half-saturation logarithm standard deviation, and (f) optimal temperature standard deviation. The dashed lines represent the environmental drivers.

With standard mutation rates, the trait dynamics are very similar in all three models. The two-trait model has slightly lower standard deviations than the one-trait models during some parts of the year, but the difference is always within 10 %. The seasonal patterns are very similar in both timing and amplitude. The greatest differences are found in summer from mid-June to mid-August, when the one-trait model with adaptive dynamics for half-saturation has a greater phytoplankton concentration by as much as 29 % and a lower nutrient concentration by as much as 24 % compared to the other two models, with concentrations that are very similar to each other. This result means that at the onset of summer, the most important factor decreasing the ability of phytoplankton to grow is not the lack of nutrients, but temperature itself. In other words, the phenotypes that dominated in spring decline, not because they are not adapted to oligotrophic conditions, but because they are not adapted to the high temperatures of summer and the growth rate of a phenotype declines sharply when temperature exceeds its optimal value. This effect is negligible at the highest mutation rates because in this case the community is able to evolve and adapt very quickly to the summer warming, but ti becomes more important as the mutation rates and hence the optimal temperature variances decrease.

The differences between models are larger at low mutation rates. With the smallest non-zero mutation rates, the summer difference in phytoplankton biomass increases. The one-trait half-saturation model has a phytoplankton biomass as much as 57 % greater and a DIN concentration as much as 30 % smaller than in the other models. Trait variances are again lower in the two-trait model during most of the year but sometimes exceed the one-trait variances during the autumn mixing. However, the most notable change is that the seasonal amplitude of mean half-saturation ($\stackrel{\mathrm{‾}}{x}$) is 56 % higher in the two-trait model than in the one-trait half-saturation model. Having a second trait allows the ecosystem to adapt faster to environmental changes. This effect is even more notable when considering that both the half-saturation variance and the nutrient-mediated selective pressure are lower in the two-trait model. This effect does not extend to the other trait, although the mean optimal temperature of the two-trait model and that of the one-trait optimal temperature model sometimes show a slight departure from each other in a seasonally dependent way.

The effects described above are related to inter-trait correlation, which is driven by correlated environmental conditions and becomes very large in the case of low mutation rates. Equations (29) to (34) can help us understand the effect of trait correlation on the seasonality of mean traits and trait variances. In the mean trait equations (30 and 31), correlation implies that both temperature and the DIN concentration drive changes in both mean trait values. The covariance term can either accelerate or slow down the response of each mean trait, but generally the sign of covariance is such that the change is accelerated. This is what occurs from December to March, when the environment selects for higher half-saturation constants and lower optimal temperatures, and the environmentally induced negative correlation between traits further accelerates this adaptation. From June to October, the same effect occurs but is significant only for half-saturation. During these months, the two-trait model actually experiences a slower increase in optimal temperature than the one-trait model because it has a smaller variance and because the selective pressure of high temperatures is much sharper than the nutrient-mediated pressure conveyed by the correlative term. In November, correlation has the opposite effect: as the temperature decreases while the DIN concentration remains low, the environment at that time selects for both low optimal temperature and low half-saturation, and the negative correlation prevents the optimal temperature from decreasing.

The effect of correlation on variance is even more convoluted. In Eqs. (32) and (33), inter-trait correlation adds a second variance-reducing competition term ($C_{xy}^{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}$ and $C_{xy}^{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}$, respectively, both very likely to be negative). This is why variances are smaller in the two-trait model during most of the year. However, one source of variance is not accounted for in these equations: vertical mixing. Trait variance is not a conservative tracer. Indeed, mixing two communities with different mean trait values “creates” additional variance. As phytoplankton adapt better to their environment in the two-trait model than in the one-trait models, the difference between surface and subsurface communities when the water column is stratified is larger in the two-trait model, and therefore the late autumn mixing event adds more trait variance in the two-trait model. This is why the variances are higher in the two-trait model than in the one-trait models in December.

4 Discussion

4.1 Strengths and weaknesses of aggregate models

SPEAD is an aggregate phytoplankton model. Aggregate models, used as far as we know since Wirtz and Eckhardt (1996), do not compute the abundance of each phytoplankton species as discrete entities, but represent the phytoplankton community by its total biomass together with the mean values, variances, and covariances of a few key traits controlling its competitive ability along different environmental gradients. Aggregate models are known to reduce the computational cost of ecosystem models by at least 1 order of magnitude. In a 0D physical setting (i.e. not spatially resolved), the single-trait aggregate model of Acevedo-Trejos et al. (2016) was found to be 18 times faster to run than its alternative discrete model with as few as 10 phenotypes. SPEAD 1.0 is 70 times faster in its aggregate version than in its alternative discrete version with $\mathrm{25}\times \mathrm{25}=\mathrm{625}$ phenotypes, despite the two model versions having to compute the same number of non-phytoplankton variables. The exact factor of cost reduction depends on the number of phenotypes used in the discrete model, and it is likely to be even larger if the number of traits increases. The computational cost of a multi-phenotype model with phenotypes covering the entire trait domain is an exponential function of the number of traits. For instance, allowing 25 values of optimal irradiance would multiply the cost of our discrete model by 25. By contrast, the cost of an aggregate model is a quadratic function of the number of traits. In our case, adding a third trait would simply increase the number of phytoplankton state variables from 6 to 10 (total phytoplankton concentration, three mean traits, three trait variances, three inter-trait covariances) and add a few new terms in each equation, which is far less computationally demanding. This makes aggregate models promising tools to explore high-dimensional trait spaces of the ecology and evolution of microbial ecotypes (Vallina et al., 2019). We note that some discrete models use an approach different from the one described in our study: they run several simulations with a relatively low number of randomly sampled phenotypes and then make an ensemble mean of all their simulations. For instance, in Follows et al. (2007), each member of the ensemble contains 78 random phenotypes, which is not a particularly large number given that they have four functional types and that their traits include half-saturation constants for several nutrients, optimal temperature, and optimal irradiance. Still, running 10 such simulations to compute an ensemble mean is computationally expensive. Aggregate models do not need ensemble means for sampling the trait space since they are continuous by design: their phytoplankton communities fill the trait space completely without the need for any arbitrary sampling.

Aggregate models are very efficient because their state variables are the quantities that make the most ecological sense. In thermodynamics, computing the trajectory of each atom or molecule is not only infeasible, but also of little use. The collection of trajectories does not provide more information on the macroscopic behaviour of a thermodynamic system than aggregate properties such as temperature, pressure, and density. Equally, modelling the dynamics of thousands of species would be incredibly costly, and sufficient observational data would not be available to validate the models. Furthermore, the results would also be extremely difficult to interpret (Levins, 1966). Given that the most important quantities for understanding a community of species with similar niches are biomass, mean trait values, and trait diversity, the aggregate model focuses computational power where it is most needed without much loss of information. Aggregate models also explicitly quantify the factors controlling biodiversity, such as the second derivatives of the net growth rate and the trait diffusivity parameters (Chen et al., 2019).

However, the aggregate approach has one major weakness: a specific shape for the trait distribution must be assumed a priori, with only as many degrees of freedom as there are free parameters (Wirtz and Eckhardt, 1996; Bruggeman and Kooijman, 2007). There is no universal distribution shape for phytoplankton traits, which is why the equations describing their dynamics are not as precise as the equations of thermodynamics. In this study, we assumed that the optimal temperature was normally distributed and that the half-saturation constant was lognormally distributed. This Gaussian closure was chosen because of its simple moment equations and low number of free or arbitrary parameters. Normally distributed temperature and irradiance niches have been observed by Irwin et al. (2015). The half-saturation constant is strongly correlated with cell size (Litchman et al., 2007; Edwards et al., 2012), and lognormal distributions have been observed in nature for size (Cermeño and Figueiras, 2008; Quintana et al., 2008; Schartau et al., 2010; Downing et al., 2014; Marañon, 2015), although not in all cases. Another very common distribution is the power (or log-linear) law (Rodríguez, 1994; Cermeño and Figueiras, 2008; Huete-Ortega et al., 2012), although the power law must be truncated on at least one side. We note that a power-law distribution with a cutoff on the left might be better able to represent the right-skewed size distribution (even in logarithmic scale) of oligotrophic environments wherein Prochlorococcus, the smallest known phytoplankton, dominates and coexists only with larger species (Marañon, 2015). What neither a unimodal nor a power-law distribution can capture is bimodality, which is known to occur at least in lakes due to common herbivores, in particular daphnids, feeding optimally on prey of intermediate size (Gaedke and Klauschies, 2017).

Normal distributions are symmetrical, unimodal, and unbounded. If the real trait distribution deviates from these three properties, errors will arise in aggregate models based on a normal distribution. Not only is information lacking by not including higher-order moments such as skewness and kurtosis, but the dynamics of mean traits and trait variances could be significantly altered. If the trait distribution is skewed, the community will respond faster to a certain type of perturbation than to the opposite perturbation. For instance, if optimal temperature is right-skewed, the phytoplankton community will adapt faster to warming than to cooling environmental conditions. Phytoplankton with larger optimal temperature will need less time to become dominant in the case of warming than cold-water phenotypes in the case of cooling because they will start from a higher concentration. To express the above in terms of moment equations, the variance of a right-skewed distribution increases when the environment favours larger trait values, thus facilitating the adaptation, but it decreases when the environment favours smaller trait values (see Appendix B and the neglected term M₃₀ in Eq. B7). If the trait distribution is multimodal, the reduction in trait diversity induced by competitive exclusion (Hardin, 1960) will be slower. This is because replacing all pre-existing communities by intermediate and previously rare phenotypes takes more time than making the most abundant phenotype even more abundant and the rarest even rarer, as in a unimodal distribution (see Appendix B and the neglected kurtosis term M₄₀ in Eq. B7, knowing that multimodal distributions have low kurtosis).

Normal distributions are unbounded, with the assumption that extreme values are rare and ecologically meaningless. The consequence of this apparently reasonable assumption is that model phytoplankton can adapt to any environmental change if they are given enough time, irrespective of the intensity of that change. This contrasts with expectations of the behaviour of real phytoplankton communities, as explained in the following example. If a closed (i.e. without immigration) phytoplankton community experiences temperatures between 15 and 25 ^∘C, the local phenotypes should be adapted to temperatures between 15 and 25 ^∘C, and no single individual should be optimized for temperatures out of the boundaries, since it would be outcompeted at all places and times. If the environment suddenly warms, the phenotypes with an optimal temperature closer to 25 ^∘C should come to dominate. However, if the temperature reaches 30 ^∘C, no phenotype with an optimal temperature of 30 ^∘C can rapidly come to dominate, since no such phenotype pre-exists in the system. Adaptation to temperatures higher than 25 ^∘C can only occur through mutations or immigration. In an aggregate model, however, none of these processes are required. Phytoplankton will be able to adapt to any warming because an extremely small but non-zero biomass of phenotypes adapted to very high temperatures is always present by model design and can become dominant if the environment selects them.

In the present study, the aggregate (continuous trait) model agrees very well with a multi-phenotype (discrete trait) model, wherein no distribution shape is imposed but the trait distribution is spontaneously close to normality during most of the year. Skewness and kurtosis occur during some times of the year, only to be removed later, and do not strongly impact our estimates of the lower-order moments. The assumed normal trait distribution is symmetrical and unimodal; therefore, some errors occur when the trait distribution is skewed or bimodal. The mean optimal temperature is slightly underestimated and its variance is overestimated because SPEAD does not account for the slight right skew of optimal temperatures distributions. The other main error is that variance tends to decrease too fast in winter, after the re-mixing of the previously stratified water column, because SPEAD cannot account for bimodality. The seasonal cycle and the orders of magnitude, however, are accurate. Our results are similar to that of Acevedo-Trejos et al. (2016). However, other studies show much larger errors (Coutinho et al., 2016; Klauschies et al., 2018) and consider Gaussian-based aggregate models to be inaccurate.

Whether the trait distribution of a model ecosystem is normal or not depends on the ecological processes included by the modeller. At least two factors in SPEAD play in favour of a normal distribution. The first factor is trait diffusion. In a fluid, the diffusive movement of a tracer follows a Gaussian law provided that the diffusivity coefficient is constant (Einstein, 1905). Trait diffusion plays the same role here for traits and tends to erase skewness and bimodality. The second factor is the simplicity of our ecological model. All our phytoplankton phenotypes compete for the same resource. In a given environment defined by nutrient concentration and temperature, there is a single most competitive phenotype, and the phytoplankton net growth rate decreases continuously when moving away from this optimum. Grazing and mortality rates do not act against this trend because we impose them to be identical for all phenotypes. Conversely, two factors in SPEAD play against normality, but their reach is relatively minor. The asymmetry of the temperature response curve promotes right-skewed distributions of optimal temperature. However, and despite the naive expectation, this right skew is generally not large in the discrete model because the standard deviation of optimal temperature is always smaller than the temperature tolerance (ΔT=5 ^∘C). The second effect is the alternation of stratification and mixing during the year. In summer, stratification leads to the formation of two distinct communities in the surface and subsurface. When the vertical mixing strengthens again in late autumn, the two communities mix into a temporary bimodal distribution.

Other ecological settings yield more widespread multimodality. Multimodality can be induced by immigration (Norberg et al., 2001), resting stages (Beckmann et al., 2019), fast environmental oscillations (Beckmann et al., 2019), spatially heterogeneous environments (Wickman et al., 2019), “convex” trade-offs favouring extreme phenotypes (Coutinho et al., 2016), and zooplankton prey selectivity (Wirtz, 2013; Klauschies et al., 2018). In particular, evolutionary branching of both phytoplankton and zooplankton into tens of size clusters can occur when each zooplankton grazes only on a small size range (Sauterey et al., 2017). However, it is important to point out that only fixed zooplankton preferences cause disruptive selection. By contrast, active switching by grazers (“kill the winner”) is insufficient to promote evolutionary branching as by design it promotes uniform distributions, flattening the peaks and filling the gaps in trait distributions. Using a Gaussian-based aggregate model in a setting promoting branching would of course be inappropriate: SPEAD would be unable to simulate evolutionary branching of phenotypes in these kinds of ecological scenarios

Alternatives to Gaussian closures have been proposed since early in the development of aggregate models. Norberg et al. (2001) and Norberg (2004) used more complex closures to estimate skewness and kurtosis. However, these closures had free parameters varying from ecosystem to ecosystem, and a discrete model was required to compute them, cancelling the advantage of aggregate models in terms of computational cost. Klauschies et al. (2018) replaced the normal distribution by a beta distribution. The beta distribution is bounded, allows bimodality, and was proven to increase the realism of trait-based aggregate models in a bounded trait scenario. However, applying this method requires defining fixed boundaries for phytoplankton traits. Phytoplankton have a minimum cell size (and half-saturation) at 0.5 µm, which is the size of Prochlorococcus, but they do not have a well-defined maximum cell size. Also, optimal temperature at local scales does not have clear boundaries. Therefore, any set of boundaries would be arbitrary and might prevent further adaptation to changing environments beyond those limits.

A more practical approach to account for non-Gaussian distributions would be to divide the community into several functional groups, each one having a normal trait distribution of its own (Terseleer et al., 2014; Chen and Laws, 2017). The sum of these communities can have a skewed or multimodal trait distribution. Trait variances (in particular size variance) could be high within phytoplankton as a whole, without bolstering the adaptive capacity of each functional group (see Sect. 4.3). All phenotypes within a given functional group must feed on the same nutrients, be subject to the same trade-offs, and should not be subject to processes promoting evolutionary branching. Ideally, and in order to prevent the convergence of all functional groups on the same trait values, each group should have distinct qualitative properties or trade-offs. Functional groups could include diatoms, mixotrophs, diazotrophs, or Prochlorococcus, among others. Two communities defined by the same parameters could coexist and avoid merging if an intermediate trait range is permanently disadvantageous due, for instance, to a convex trade-off or to a size-specific grazer. The multi-Gaussian approach would combine the moderate computational cost of Gaussian aggregate models with the more thorough description of planktonic ecosystems allowed by discrete models.

4.2 Traits in phytoplankton community models

In nature, many different traits define phytoplankton niches: nitrogen, phosphorus or iron uptake abilities, requirements for other nutrients (for instance, silica or calcite), stoichiometry, optimal temperature, optimal irradiance, mixotrophy, diazotrophy, motility, buoyancy, resistance to predation, toxicity, and many others. In many trait-based models, this complexity is reduced to one trait. The most common trait is cell size (Terseleer et al., 2014; Acevedo-Trejos et al., 2016; Smith et al., 2016; Chen and Smith, 2018). Cell size is used as a master trait because it is the most observable trait and correlates strongly with many other phytoplankton traits, such as light requirements (Taguchi, 1976; Edwards et al., 2015; Álvarez et al., 2017) and resistance to predation (Kiørboe, 1993; Thingstad et al., 2005). Some other commonly modelled traits include resistance to predation (Norberg et al., 2001; Merico et al., 2009), optimal temperature (Norberg et al., 2012; Beckmann et al., 2019), and optimal irradiance (Follows and Dutkiewicz, 2011).

The first trait included in SPEAD, half-saturation for nutrients, is known to be strongly correlated with cell size (Litchman et al., 2007; Edwards et al., 2012). Small species, such as cyanobacteria, have low half-saturation constants and can thrive in oligotrophic waters. They correspond to the gleaners of our model. Large species are more likely to be involved in blooms but require more nutrients. They correspond to the opportunists of our model. In this study, we chose not to use size but to impose only a simple gleaner–opportunist trade-off. This choice was made for two reasons: to maintain compatibility with the Darwin model (Follows et al., 2007; Vallina et al., 2014a) and to facilitate the analysis of the outputs of our otherwise complex model, as in our setting the best competitors in a given environment have a half-saturation equal to the dissolved inorganic nitrogen concentration. Our modelling framework can, however, be easily adapted to use size instead of half-saturation as the first trait (Smith et al., 2016).

The trait dynamics of models with two traits differ from those of simpler and less realistic single-trait models. Savage et al. (2007) obtain larger trait variances and much larger adaptive capacities when two traits are modelled together rather than in separate models. In our study, we also find a larger adaptive capacity, although it is conveyed by inter-trait correlation only. We actually find decreased variances caused by stronger competition during most of the year, except during the late autumn re-mixing of the water column following the summer stratification. This discrepancy might have been caused by the presence of an immigration term to sustain variance in Savage et al. (2007) but not in SPEAD, since our re-mixing of the water column is most analogous to a dispersal or migration process and we simulate it explicitly. SPEAD simulations coupled with a realistic 3D circulation model wherein phytoplankton are explicitly allowed to migrate in all directions (ideally in a patchy environment) would finally tell us if variance is increased or decreased by including more traits.

The low computational cost of aggregate models allows increasing the number of modelled traits, provided that sufficient observational data are available to constrain the corresponding trade-offs. Since the environmental drivers, such as nutrient concentrations, temperature, and light, are correlated with each other, the traits are likely to be correlated unless some processes erasing the correlations are introduced. Regardless of the effect of interactions between traits on variance, multi-trait models will be able to adapt to their environments faster without the need for large and unrealistic mutation rates or other terms sustaining large variances, such as immigration or “kill the winner” grazing.

4.3 Trait diffusion, variance, and evolution

Trait diffusion is a key process in SPEAD. Indeed, SPEAD is the first model to include diffusion of multiple traits, providing insights into how both mutations and selection can impact phytoplankton communities. For each modelled trait, a diffusivity parameter, or mutation rate, has to be set. The chosen values of these parameters decisively affect trait dynamics. However, the mutation rates remain poorly constrained. The most appropriate rate depends on what the modeller intends to represent by trait diffusion. To understand why, we will need to discuss the notions of ecological and evolutionary timescales, as well as the notions of adaptation and species.

A first interpretation of trait diffusion is that it is the most conservative way to add variance when the exact processes sustaining trait variance are unknown or too complex to be implemented in models. Indeed, trait diffusion simply adds new variance (aggregate approach) or disperses phytoplankton in a trait space (discrete approach), leaving little room to arbitrary parameters. By contrast, immigration (Norberg et al., 2001) requires assumptions about both immigration rates and the composition of immigrant communities. Likewise, “kill the winner” (Vallina et al., 2014b) requires assumption about whether it is mediated by viruses or zooplankton, whether consumers are specialized or can actively switch their preference between prey, and whether or not there is a time lag in their response.

This interpretation of trait diffusion as a generic source of variance is implicitly followed when the diffusivity parameter is set by an optimization algorithm in order to account for the observed trait variance and no other mechanism sustaining variance is included. This way, Chen and Smith (2018) found a diffusivity for the logarithm of size of 0.1, equal to our largest diffusivity for the logarithm of half-saturation. Even in the homogeneous environments of mesocosms, Wirtz (2013) reports a logarithmic size variance of 0.2 to 0.5, which corresponds to standard deviations between 0.45 and 0.7. In SPEAD, reaching these high values of trait variance is only possible with diffusivities higher than or equal to 0.003, despite the fact that our physical setting creates its own variance by mixing phenotypes adapted to the environmental conditions of different depths. This interpretation of trait diffusion is also coherent with the use of trait diffusion as a variance treatment (Chen et al., 2019) to study the effect of diversity on primary production and with the original goal of trait diffusion, which was to sustain trait variance in 0D settings (Merico et al., 2014; Acevedo-Trejos et al., 2016). In these models, trait diffusion is never run combined with other variance-sustaining mechanisms. In real ecosystems, however, mutations are expected to occur at the same time as migration, “kill the winner” grazing, and many other mechanisms that may promote diversity, including multiple convex trade-offs (Beardmore et al., 2011) and mixotrophy (Ward and Follows, 2016). The way trait diffusion is derived opens a second interpretation of what it represents. Trait diffusion is symmetrical: mutations occur at the same rate towards higher and lower trait values. Trait diffusion is also heritable: mutants transfer their mutations to their offspring. These properties correspond to the evolutionary process of random mutations and selection of the fittest by the environment. They do not correspond to environmentally induced non-heritable variations such as phenotypic plasticity (Ghalambor et al., 2007) or to any selective ecological process driven by the environment, even if some promote variance.

According to Fussmann et al. (2007), “Evolution is the change of genotype frequencies within populations or species, whereas community dynamics represent the change of abundances of different species”. This definition depends on the notion of species. Like many phytoplankton models, SPEAD lacks the notion of species and does not distinguish between intraspecific and interspecific trait variance. Our trait space is continuous by design. By mutating, phytoplankton can cross the boundaries between phenotypes as if they were all of the same species. The only distinction in SPEAD is between mutations and selection. Mutations are a strictly evolutionary process represented in SPEAD by trait diffusion through subsequent generations. Selection encompasses both intraspecific selection, which is a second evolutionary process, and interspecific selection, which is an ecological process. In SPEAD, selection has two effects: it drives the mean traits towards their optimum values (see Eqs. 30–31) and decreases trait variances by eliminating the rarest and least fit phenotypes (see Eqs. 32–34). Some cases of adaptive evolution to environmental changes in only a few generations have been reported (Fussmann et al., 2007; Kinnison and Hairston, 2007) but are not necessarily caused by mutations occurring at these timescales. They can also be driven by ecological selection from previously existing intraspecific diversity. In order to choose correct mutation rates these two evolutionary processes must be distinguished.

An alternative interpretation is that ecological processes are particular cases of eco-evolutionary processes whereby the phenotypes of the offspring are identical to those of their parents (Doebeli et al., 2017). This definition emphasizes the lack of a fundamental difference between the two types of processes from a mechanistic point of view, since both are rooted in the same birth–death dynamics. Under this alternative interpretation, ecological timescales are simply the timescales at which the effect of mutations is small, and evolutionary timescales are those at which the effect of mutations is large. Therefore, the ecological selection timescales overlap the adaptive evolution timescales, and the difference between the two is diffuse at intermediate timescales. Schlüter et al. (2016) showed that an algal culture starting with a single clone of the abundant coccolithophore Emiliania huxleyi could evolve new traits in response to ocean acidification in a time measurable in the laboratory. Their experiment lasted for 2100 generations, and the changes after only 100 generations were small. These numbers agree with seminal studies on Escherichia coli wherein bacteria were shown to adapt to temperature increases or to changes in nutrient availability in 100 to a few thousand generations (Bennett et al., 1990; Lenski et al., 1991; Travisano et al., 1995). These are the timescales a mutation rate should reflect. In the case of phytoplankton this means a few years, which falls under the category of “contemporary evolution” but does not allow each species to adapt easily to a seasonal cycle or to faster perturbations. The corresponding trait diffusivity parameters in our study are in the middle of our range, between $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ and 10⁻³ for half-saturation and between $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ and 10^−2 ∘C² for optimal temperature (see Table 3).

The trait distributions in SPEAD provide additional insights. In the absence of trait diffusion, the two traits become almost totally correlated: one cannot vary without the other. This is a soft version of the diversity collapse observed in 0D models of a one-trait fitness landscape. In our 1D model of a two-trait fitness landscape, trait variances do not collapse to zero, but a bi-dimensional trait space becomes unidimensional and phytoplankton lose their ability to adapt in other directions. Nutrient and temperature niches are known to be correlated in nature, but their correlation is never perfect (Irwin et al., 2012). A small trait diffusivity is sufficient to avoid the collapse of the dimensional trait space into a unidimensional one and likely limits such trait correlations in natural ecosystems, given that mutations affecting half-saturation and optimal temperature are likely independent and hence able to freely fill the full trait space. With very fast trait diffusion, the mean phenotypes adapt instantaneously to their environment but at the cost of keeping a large pool of maladapted phenotypes, which regularly represent more than 15 % of the community even if they have very low fitness because mutations are continuously creating them. These maladapted phenotypes explain the decrease by up to 10 % of the modelled annual primary production. In nature, the optimal niches of species do not cover all the variability of nutrient concentration and water temperature (Irwin et al., 2012). Furthermore, most real mutations having an effect on the phenotype are deleterious (Timofeeff-Ressovsky, 1940). A mutation rate able to permanently sustain maladapted phenotypes despite strong selection against them would imply a large number of deleterious mutations not represented in our model and hence massive mortality. Moderate mutation rates are therefore more likely.

Modelling several communities with their own trait distributions and their own mutations might relax the contradiction between the use of trait diffusion to explain trait variance and the use of trait diffusion to represent evolutionary processes. The variance within a species or a group is lower than the total community variance and can be sustained with lower mutation rates. This approach can also be used to separate the adaptive evolution of each species from the ecological successions (i.e. inter-group competition) in response to environmental change (Norberg et al., 2012), finally disentangling all components of eco-evolutionary processes.

4.4 Future directions

SPEAD 1.0 is the first step of the SPEAD project, whose aim is to simulate plankton evolution with adaptive dynamics in the ocean. In this first version of the model, we kept the complexity manageable, with only one spatial dimension (the vertical) and two functional traits, in order to facilitate the validation of our aggregate approach and to diagnose the effects of trait diffusion. The equations of SPEAD 1.0 for the mean trait, trait variance, and covariance resolving two functional traits can be used as a starting point to build more comprehensive trait-based models in multi-dimensional continuous trait spaces, with or without mutations. Three axes of potential future improvement have already been identified: (1) coupling SPEAD with a general circulation model, (2) increasing the number of traits, and (3) dividing the community into several functional groups, which implies combining the continuous trait distribution approach with the discrete ecotype approach.

More concretely, our goal for the near future is to include optimal solar irradiance as a third functional trait and implement the aggregate approach with trait diffusion in a 3D trait space into the Darwin model (Follows et al., 2007; Dutkiewicz et al., 2009; Barton et al., 2010; Follows and Dutkiewicz, 2011; Ward et al., 2012; Dutkiewicz et al., 2013). Darwin is a versatile model that allows many discrete ecotypes to be resolved along several environmental gradients and can be coupled with the MIT general circulation model (Marshall et al., 1997). Optimal irradiance has been present as a trait in the Darwin model since its origin. Indeed, light can be a limiting resource for phytoplankton in the ocean in mixed water columns, where plankton cannot stay close to the surface, and near the deep chlorophyll maximum of stratified water columns. Optimal irradiances are known to cover more than 2 orders of magnitude (Edwards et al., 2015) and to determine the niches of many ecotypes, even within the same species (Biller et al., 2015). The inclusion of optimal solar irradiance as a mutating trait in SPEAD will be a key step to better capture ecological successions. One challenge will be to allow phenotypes to adapt to their environment while at the same time accounting for the mechanistic correlation between nutrient and irradiance niches, since both are related to cell size. The phytoplankton in Darwin were originally divided into four functional groups: Prochlorococcus analogues, other small phytoplankton, diatoms, and other large phytoplankton. Representing each functional group by its own normal distribution is a possible starting point to develop the multi-Gaussian approach.

Improved versions of SPEAD should be able to address various ecological issues related to community assembly and responses to climate change that current models cannot address. The response of the phytoplankton community to environmental changes in a simple NPZD model can only be an increase or a decrease in primary production. Estimates of changes in primary production from NPZD models are likely to be inaccurate because the parameters of NPZD models are validated against observations of present environments. Under environmental changes, the composition of the plankton community is likely to change and the model parameters have no reason to remain valid. Models based on plankton functional types are more accurate, as they can account for ecological selection: a plankton type (e.g. diatoms) can be replaced by another type (e.g. cyanobacteria) in response to a perturbation (e.g. increased stratification). However, they do not represent adaptive evolution through which groups or species can change their traits and maintain local dominance. Therefore, they might overestimate the extinction rate and the shift in community composition. Eco-evolutionary models like SPEAD do not explain why groups that diverged millions of years ago exist, but a multi-Gaussian version of SPEAD, wherein each group follows its own adaptive dynamics, could account for their contemporary evolution in new environments. By including both trait diffusion and ecological selection of the fittest phenotypes competing in a given environment, SPEAD can potentially be used to disentangle the role of ecological and evolutionary processes in shaping diversity patterns in phytoplankton. In particular, it can be used to determine the conditions under which species or functional groups may survive climate change by evolving new traits or may be replaced by other species or functional groups from other regions. Predicting changes in phytoplankton composition is particularly important as species perform different functions and have different impacts on their environment. For instance, they contribute differently to the carbon pump (e.g. by sinking more or less fast) and to the nitrogen cycle (e.g. by fixing atmospheric nitrogen or not). Therefore, changes in community composition might dramatically impact global climate and should be included in climate prediction models. The effect of environmental changes, such as warming and increased stratification, on plankton size structure and the effect of biodiversity – controlled by trait diffusion among other processes – on primary production and ecosystem functioning are other examples of contemporary ecological issues that SPEAD might contribute to addressing.

5 Conclusions

In this article, we present a new aggregate model of a phytoplankton community called SPEAD (Simulating Plankton Evolution with Adaptive Dynamics), wherein different phenotypes competing for dissolved inorganic nitrogen are characterized by two traits: their half-saturation constants for nitrogen uptake (in logarithmic scale) and their optimal temperature for growth. The phytoplankton community is represented by the six lowest-order moments of its trait distribution: total concentration, the mean value of each trait, the variance of each trait, and the inter-trait covariance. The dynamics of these state variables are driven by three environmental factors: nutrient concentration, temperature, and solar irradiance. The physical setting represents a water column down to 200 m. The seasonal alternation of stratification and vertical mixing also has a strong effect on the trait distribution. Trait diffusion through subsequent generations is included to represent heritable mutations and hence sustain trait diversity. To our knowledge, SPEAD is the first aggregate model to simultaneously include two traits (with a proper representation of inter-trait correlation) and trait diffusion.

The ecological parameters of SPEAD were set to reproduce the observed primary production as well as the chlorophyll, nitrate, and particulate organic nitrogen concentrations observed by BATS in the Sargasso Sea. Despite its strong assumption that traits are normally distributed, SPEAD was shown to agree precisely with a discrete model explicitly representing all phenotypes, with only minor deviations at depth in summer when optimal temperature is underestimated and in early winter when trait variances decrease too fast. This good agreement is made possible by trait diffusion and by the simplicity of our ecological setting, and it might not be extendable to all ecosystem models. The trait dynamics depend strongly on the imposed trait diffusivity parameters. With very high diffusivities, primary production is low, variances are high, and the two traits are independent, filling the entire trait space. With very low diffusivities, variances are low (albeit non-zero) and the two traits are very strictly correlated: only warm-water gleaners and cold-water opportunists can survive. We think that intermediate values of mutation rates are more realistic, but the precise value depends on whether trait diffusion is meant to sustain the trait diversity of a whole community or to represent the mutations occurring within a given species.

SPEAD has a computational cost 2 orders of magnitudes lower than a full discrete model, and its variables are readily interpretable in ecological terms. This effectiveness makes it possible to increase the number of traits. As optimal irradiance is key to explaining phytoplankton distribution in the water column and is already present in the Darwin model, the next step of the SPEAD project will be to include it as a third dynamic trait. In agreement with Savage et al. (2007), we showed that adding traits accelerated the response of pre-existing traits to environmental changes. Other avenues of future improvement include representing various functional groups, each with their own distinct normal distributions, and coupling SPEAD with a general circulation model. Future versions of our multi-trait framework may address ecological questions related to the impact of selection, mutations, and biodiversity on community dynamics and to the response of phytoplankton to climate change.

Appendix A: Why mutations can be represented as trait diffusion

In this study, we represented phytoplankton mutations as a trait diffusion term, following the work of Merico et al. (2014). In this Appendix we show how the expression for trait diffusion is derived and discuss its conditions of validity.

Let us consider the dynamics of an isolated phytoplankton community wherein each individual is characterized by the values of two traits: x (in trait unit x or tux) and y (in tuy). Phytoplankton are distributed in the trait space with a density $p(x,y,t)$ (mmol N m⁻³ tux⁻¹ tuy⁻¹). The mass concentration of phytoplankton cells with values of the first trait between x and x+dx and values of the second trait between y and y+dy is $p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y$ if dx and dy are small.

If traits are strictly inherited, the equation governing $p(x,y,t)$ for a given phenotype (x,y) depends on the reproduction ($u(x,y,t)$; d⁻¹) and death ($d(x,y,t)$; d⁻¹) rates.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial p}{\partial t}}(x,y,t)=\left(u\right(x,y,t)-d(x,y,t\left)\right)p(x,y,t)\\ {\displaystyle }& {\displaystyle \frac{\partial p}{\partial t}}(x,y,t)=a(x,y,t)p(x,y,t)\end{array}$$

In the above equation, $a(x,y,t)=u(x,y,t)-d(x,y,t)$ is the net growth rate. In our model, the reproduction rate is identical to the nitrogen uptake rate because nitrogen, the limiting nutrient, is not exuded, cells do not modify their nitrogen content, and all nitrogen taken up is used for reproduction. However, genetic mutations or phenotypic plasticity can produce offspring with trait values different from that of their parents. For simplicity, we will consider mutations increasing or decreasing the traits to be equally probable. We assume that the offspring of a parent with trait value x will have trait value x−δ_x with a probability $a_x={\mathit{\nu }}_x{\left(\mathit{\delta }x\right)}^{-\mathrm{2}}$ and trait value x+δ_x with the same probability, where ν_x is a diffusivity parameter expressed in tux² that is considered independent of trait value, mutation step (δ_x), and time. Mutations also occur on trait y, with a mutation step δy and a y-diffusivity parameter ν_y. Mutations on both traits are assumed independent. The probability of having a mutation on both x and y is just the product of the probabilities of each mutation. Hence, the time derivative of $p(x,y,t)$ with mutations is as follows.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial p}{\partial t}}(x,y,t)=\left[(\mathrm{1}-\mathrm{2}{a}_x)(\mathrm{1}-\mathrm{2}{a}_y)u(x,y,t)-\mathrm{d}(x,y,t)\right]\\ {\displaystyle }& {\displaystyle }p(x,y,t)\\ {\displaystyle }& {\displaystyle }+(\mathrm{1}-\mathrm{2}{a}_y)a_x\left[u(x-\mathit{\delta }x,y,t)p(x-\mathit{\delta }x,y,t)\right.\\ {\displaystyle }& {\displaystyle }\left.+u(x+\mathit{\delta }x,y,t)p(x+\mathit{\delta }x,y,t)\right]\\ {\displaystyle }& {\displaystyle }+(\mathrm{1}-\mathrm{2}{a}_x)a_y\left[u(x,y-\mathit{\delta }y,t)p(x,y-\mathit{\delta }y,t)\right.\\ {\displaystyle }& {\displaystyle }\left.+u(x,y+\mathit{\delta }y,t)p(x,y+\mathit{\delta }y,t)\right]\\ {\displaystyle }& {\displaystyle }+a_x{a}_y\left[u\right(x-\mathit{\delta }x,y-\mathit{\delta }y,t\left)p\right(x-\mathit{\delta }x,y-\mathit{\delta }y,t)\\ {\displaystyle }& {\displaystyle }+u(x+\mathit{\delta }x,y-\mathit{\delta }y,t)p(x+\mathit{\delta }x,y-\mathit{\delta }y,t)\\ {\displaystyle }& {\displaystyle }+u(x-\mathit{\delta }x,y+\mathit{\delta }y,t)p(x-\mathit{\delta }x,y+\mathit{\delta }y,t)\\ {\displaystyle }& {\displaystyle }+u(x+\mathit{\delta }x,y+\mathit{\delta }y,t)p(x+\mathit{\delta }x,y+\mathit{\delta }y,t)]\end{array}$$

In the limit of small but frequent mutations, this equation can be simplified by making a second-order approximation of u⋅p.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial p}{\partial t}}(x,y,t)=\left[(\mathrm{1}-\mathrm{2}{a}_x)(\mathrm{1}-\mathrm{2}{a}_y)u(x,y,t)-d(x,y,t)\right]\\ {\displaystyle }& {\displaystyle }p(x,y,t)\\ {\displaystyle }& {\displaystyle }+(\mathrm{1}-\mathrm{2}{a}_y)a_x\left[\left(u\cdot p-{\mathit{\delta }}_x{\displaystyle \frac{\partial (u\cdot p)}{\partial x}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+\left(u\cdot p+{\mathit{\delta }}_x{\displaystyle \frac{\partial (u\cdot p)}{\partial x}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}\right)\right]\\ {\displaystyle }& {\displaystyle }+(\mathrm{1}-\mathrm{2}{a}_x)a_y\left[\left(u\cdot p-{\mathit{\delta }}_y{\displaystyle \frac{\partial (u\cdot p)}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+\left(u\cdot p+{\mathit{\delta }}_y{\displaystyle \frac{\partial (u\cdot p)}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}\right)\right]\\ {\displaystyle }& {\displaystyle }+a_x{a}_y\left[u\cdot p-{\mathit{\delta }}_x{\displaystyle \frac{\partial (u\cdot p)}{\partial x}}-{\mathit{\delta }}_y{\displaystyle \frac{\partial (u\cdot p)}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}\right.\\ {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}+{\mathit{\delta }}_x{\mathit{\delta }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x\partial y}}\right]\\ {\displaystyle }& {\displaystyle }+a_x{a}_y\left[u\cdot p+{\mathit{\delta }}_x{\displaystyle \frac{\partial (u\cdot p)}{\partial x}}-{\mathit{\delta }}_y{\displaystyle \frac{\partial (u\cdot p)}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}\right.\\ {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}-{\mathit{\delta }}_x{\mathit{\delta }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x\partial y}}\right]\\ {\displaystyle }& {\displaystyle }+a_x{a}_y\left[u\cdot p-{\mathit{\delta }}_x{\displaystyle \frac{\partial (u\cdot p)}{\partial x}}+{\mathit{\delta }}_y{\displaystyle \frac{\partial (u\cdot p)}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}\right.\\ {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}-{\mathit{\delta }}_x{\mathit{\delta }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x\partial y}}\right]\\ {\displaystyle }& {\displaystyle }+a_x{a}_y\left[u\cdot p+{\mathit{\delta }}_x{\displaystyle \frac{\partial (u\cdot p)}{\partial x}}+{\mathit{\delta }}_y{\displaystyle \frac{\partial (u\cdot p)}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}\right.\\ {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}+{\mathit{\delta }}_x{\mathit{\delta }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x\partial y}}\right]\end{array}$$

The sum of terms in $\frac{\partial (u\cdot p)}{\partial x}$, $\frac{\partial (u\cdot p)}{\partial y}$, and $\frac{\partial (u\cdot p)}{\partial x\partial y}$ is zero, and the sum of factors before u⋅p is 1. Hence, the above equation simplifies to the following.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial p}{\partial t}}(x,y,t)=\left[u(x,y,t)-d(x,y,t)\right]p(x,y,t)\\ {\displaystyle }& {\displaystyle }+a_x{\mathit{\delta }}_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}+a_y{\mathit{\delta }}_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}\\ {\displaystyle }=& {\displaystyle }a(x,y,t)p(x,y,t)+{\mathit{\nu }}_x{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}+{\mathit{\nu }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}\end{array}$$

This second-order approximation is valid in the limit as mutations become small (δ_x and δ_y tend to zero) and frequent ($a_x={\mathit{\nu }}_x{\left(\mathit{\delta }x\right)}^{-\mathrm{2}}$ and $a_y={\mathit{\nu }}_y{\left(\mathit{\delta }y\right)}^{-\mathrm{2}}$ tend to infinity) so that higher-order terms can be neglected. The mathematical expression of the mutation term $\left({\mathit{\nu }}_x\frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}+{\mathit{\nu }}_y\frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}\right)$ is somewhat analogous to that representing the diffusion of a tracer in physical space. In this analogy, the trait space replaces the physical space, ν_x⋅u and ν_y⋅u replace the diffusivity, and the density $p(x,y,t)$ is the diffused tracer. This analogy is the origin of the phrase “trait diffusion”. We note that in a multi-trait space, whatever the number of traits, the diffusion of each trait has the same expression as in the one-trait space.

Appendix B: Differential equations of a multi-trait aggregate model with trait diffusion

Phytoplankton community models can be discrete or aggregate. In a discrete model, the phytoplankton community is divided into a finite number of phenotypes, each characterized by a different set of trait values. Mutations are discrete with steps equal to the difference between a phenotype and its nearest neighbours. The differential equation for a discrete phenotype is intuitive and depends only on its net growth rate and a trait diffusion term.

The variables of aggregate models and the differential equations they follow are less intuitive. In an aggregate model, a general shape must be assumed for the trait distribution with some degrees of freedom, and the prognostic variables are the moments of the trait distribution that are free to vary. In a single-trait model, the most commonly assumed shape is the normal (or Gaussian) distribution (Wirtz and Eckhardt, 1996; Bruggeman and Kooijman, 2007), wherein the phytoplankton density p(x,t) (mmol N m⁻³ tux⁻¹, with tux or the trait unit x) is equal to

$$\begin{array}{}\text{(B1)}& p(x,t)={\displaystyle \frac{P\left(t\right)}{\sqrt{\mathrm{2}\mathit{\pi }{V}_x\left(t\right)}}}{e}^{-\frac{{\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{\mathrm{2}}}{\mathrm{2}{V}_x\left(t\right)}}.\end{array}$$

In this distribution, the three free parameters are the phytoplankton concentration P(t), the mean trait value $\stackrel{\mathrm{‾}}{x}\left(t\right)$, and the trait variance V_x(t). They are respectively equal to $\int p(x,t)\cdot \mathrm{d}x$, $\int x\cdot p(x,t)\cdot \mathrm{d}x$, and $\int {\left(x-\stackrel{\mathrm{‾}}{x}\right)}^{\mathrm{2}}p(x,t)\cdot \mathrm{d}x$. In a multi-trait space, we will assume that the traits follow a multivariate normal distribution, which is a generalization of a normal distribution. If N is the number of dimensions, then

$$\begin{array}{}\text{(B2)}& p(x,t)={\displaystyle \frac{P\left(t\right)}{{\left(\mathrm{2}\mathit{\pi }\right)}^{\frac{N}{\mathrm{2}}}\det (\mathrm{\Sigma }{)}^{\frac{\mathrm{1}}{\mathrm{2}}}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}}\left(t\right)\right)}^T\mathbf{\Sigma }(t{)}^{-\mathrm{1}}\left(\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}}\left(t\right)\right)}.\end{array}$$

In this case, x is the vector containing all traits, $\stackrel{\mathrm{‾}}{\mathit{x}}\left(t\right)$ is the vector of mean trait values, and Σ(t) is the matrix of variances and covariances. There are $\frac{(N+\mathrm{1})(N+\mathrm{2})}{\mathrm{2}}$ free parameters in total: one phytoplankton concentration, N mean trait values, N trait variances, and $\frac{N(N-\mathrm{1})}{\mathrm{2}}$ covariances. In the following, we will show how to derive the differential equations for each type of variable. To this end, we use the method developed by Norberg et al. (2001) based on Taylor expansions of the rates of reproduction and net growth. The assumption of a normal trait distribution is only required to compute the time derivatives of variance and covariance, but not for the equations of total biomass and mean trait values.

The trait space is considered to be unbounded, with the implicit assumption that extreme values are extremely rare and ecologically meaningless. This is expressed in the fact that p and all products including p or any of its derivatives tend to 0 when a trait tends to (plus or minus) infinity. For simplicity, in the following part of this section, we will not show the dependencies on the environmental factors and will limit ourselves to two traits, but our method can be extended to derive the equations for any given number of traits. The reproduction and net growth rates of the phenotype defined by trait values (x,y) are denoted u(x,y) and a(x,y), respectively. Integrals are over the whole bi-dimensional domain. We will first derive the equations in the absence of trait diffusion and then discuss what terms are added by the trait diffusion scheme. The net growth of a given phenotype is

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t}}=a(x,y,t)p(x,y,t).\end{array}$$

As a consequence, the equation controlling P(t) is

$${\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}=\int \int a(x,y,t)p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y.$$

We will use the notation a_jk (d⁻¹ tux^−j tuy^−k) for $a(x,y,t)$ derived j times with respect to x and k times with respect to y, normalized by the factorials of j and k, and the notation M_jk (tux^j tuy^k) for the central moment of order j with respect to x and k with respect to y.

$$a_{jk}\left(t\right)={\displaystyle \frac{\mathrm{1}}{j\mathrm{!}k\mathrm{!}}}{\displaystyle \frac{{\partial }^{j+k}a}{\partial x^j\partial y^k}}\left(\stackrel{\mathrm{‾}}{x}\right(t),\stackrel{\mathrm{‾}}{y}(t),t)$$

$$M_{jk}\left(t\right)={\displaystyle \frac{\mathrm{1}}{P\left(t\right)}}\int \int (x-\stackrel{\mathrm{‾}}{x}(t){)}^j\cdot (y-\stackrel{\mathrm{‾}}{y}\left(t\right){)}^k\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y$$

We note that M₀₀=1 (by definition of the total concentration), M₁₀=0, M₀₁=0 (by definition of the mean traits), M₂₀(t)=V_x(t), M₀₂(t)=V_y(t) (by definition of the variance), and M₁₁(t)=C_xy(t).

A Taylor expansion of the net growth rate on x and y around $\left(\stackrel{\mathrm{‾}}{x}\left(t\right),\stackrel{\mathrm{‾}}{y}\left(t\right)\right)$ yields

$$\begin{array}{ll}{\displaystyle }a(x,y,t)& {\displaystyle }=\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{\displaystyle \frac{\mathrm{1}}{j\mathrm{!}k\mathrm{!}}}{\displaystyle \frac{{\partial }^{j+k}a}{\partial x^j\partial y^k}}\left(\stackrel{\mathrm{‾}}{x}\left(t\right),\stackrel{\mathrm{‾}}{y}\left(t\right),t\right)\\ {\displaystyle }& {\displaystyle }\cdot {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^j\cdot {\left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)}^k\\ {\displaystyle }& {\displaystyle }=\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)\cdot {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^j\cdot {\left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)}^k.\end{array}$$

The time derivative of P(t) depends on the time derivative of p(x,t). Unless explicitly indicated otherwise, all derivatives with respect to traits are taken at the current mean trait values.

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}& {\displaystyle }=\int \int {\displaystyle \frac{\partial p}{\partial t}}(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\int \int a(x,y,t)p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\int \int \sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)\cdot {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^j\cdot {\left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)}^k\\ {\displaystyle }& {\displaystyle }\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)\int \int {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^j\cdot {\left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)}^k\\ {\displaystyle }& {\displaystyle }\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=P\left(t\right)\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)M_{jk}\left(t\right)\end{array}$$

The first and largest term of this sum is the net growth rate at the mean trait values, denoted $a\left(\stackrel{\mathrm{‾}}{x}\right(t),\stackrel{\mathrm{‾}}{y}(t),t)$. This is the expected growth rate of a community without trait variance. Since M₁₀=0 and M₀₁=0, there is no term depending on $\frac{\partial a}{\partial x}$ or $\frac{\partial a}{\partial y}$. As we want to estimate the effect of trait diversity on the community, we consider the second-order terms proportional to V_x(t), V_y(t), or C_xy(t). Higher-order terms, which vanish when variance is small, are neglected so that

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}=P\left(t\right)\left(a\left(\stackrel{\mathrm{‾}}{x}\right(t),\stackrel{\mathrm{‾}}{y}(t),t)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_x\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_y\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}+C_{xy}\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}\right).\end{array}\end{array}$$

Second-order derivatives are expected to be negative if $\left(\stackrel{\mathrm{‾}}{x}\right(t),\stackrel{\mathrm{‾}}{y}(t\left)\right)$ is in the neighbourhood of the optimal trait value, and V_x(t) and V_y(t) are always positive. Therefore, the second-order terms are generally negative. This means that having large trait variances, or in other words having a large proportion of cells with non-optimal trait values, has a negative effect on the net community growth.

In the equation for the mean trait, however, having a large variance is an advantage. Let us define an intermediate variable S_x(t) (mmol N m⁻³ tux) as $P\left(t\right)\stackrel{\mathrm{‾}}{x}\left(t\right)$ or, equivalently, as $\int \int x\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y$. The time derivative of S_x(t) is as follows.

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{S}_x}{\mathrm{d}t}}& {\displaystyle }=\int \int x{\displaystyle \frac{\partial p}{\partial t}}(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\int \int x\cdot p(x,y,t)a(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\int \int \left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)\cdot p(x,y,t)a(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }+\stackrel{\mathrm{‾}}{x}\left(t\right)\int \int p(x,y,t)a(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\int \int \sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)\cdot {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{j+\mathrm{1}}\cdot {\left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)}^k\\ {\displaystyle }& {\displaystyle }\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y+{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}\stackrel{\mathrm{‾}}{x}\left(t\right)\\ {\displaystyle }& {\displaystyle }=\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)\int \int {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{j+\mathrm{1}}\cdot {\left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)}^k\\ {\displaystyle }& {\displaystyle }\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y+{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}\stackrel{\mathrm{‾}}{x}\left(t\right)\\ {\displaystyle }& {\displaystyle }=P\left(t\right)\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)M_{j+\mathrm{1},k}\left(t\right)+{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}\stackrel{\mathrm{‾}}{x}\left(t\right)\end{array}$$

As S_x(t) is a product, its derivative can also be written as

$${\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{x}}{\mathrm{d}t}}P\left(t\right)+{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}\stackrel{\mathrm{‾}}{x}\left(t\right).$$

By equating the two previous expressions, we get

$${\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{x}}{\mathrm{d}t}}=\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)M_{j+\mathrm{1},k}\left(t\right).$$

In this equation, we only consider the highest-order terms:

$$\begin{array}{}\text{(B5)}& {\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{x}}{\mathrm{d}t}}=V_x\left(t\right){\displaystyle \frac{\partial a}{\partial x}}+C_{xy}\left(t\right){\displaystyle \frac{\partial a}{\partial y}}.\end{array}$$

This equation represents the adaptation of the community to its environment. If the mean trait is not optimal, it will increase or decrease in order to maximize the net specific growth rate. The speed of this selection process is proportional to variance: biodiversity is required to track the environmental conditions. The covariance term means that if traits are correlated, the optimization of trait y will also affect trait x.

The mean value of trait y follows a similar equation:

$$\begin{array}{}\text{(B6)}& {\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{y}}{\mathrm{d}t}}=V_y\left(t\right){\displaystyle \frac{\partial a}{\partial y}}+C_{xy}\left(t\right){\displaystyle \frac{\partial a}{\partial x}}.\end{array}$$

The equations describing time changes in variances and covariance require more assumptions. As previously, we define an intermediate variable Z_x(t) (mmol N m⁻³ tux²) as P(t)V_x(t) or, equivalently, as $\int \int {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{\mathrm{2}}\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y$. In this integral, two terms depend on time: $\stackrel{\mathrm{‾}}{x}\left(t\right)$ and $p(x,y,t)$. Hence, the time derivative of Z_x(t) is as follows.

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{Z}_x}{\mathrm{d}t}}& {\displaystyle }=\int \int (x-\stackrel{\mathrm{‾}}{x}(t){)}^{\mathrm{2}}\cdot {\displaystyle \frac{\partial p}{\partial t}}(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }+\int \int {\displaystyle \frac{\partial (x-\stackrel{\mathrm{‾}}{x}(t){)}^{\mathrm{2}}}{\partial t}}p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\int \int {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{\mathrm{2}}\cdot p(x,y,t)a(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }-\mathrm{2}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{x}}{\mathrm{d}t}}\int (x-\stackrel{\mathrm{‾}}{x})\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\end{array}$$

By definition of $\stackrel{\mathrm{‾}}{x}\left(t\right)$, we have $\int \int \left(x-\stackrel{\mathrm{‾}}{x}\right)\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y=\mathrm{0}$. Thus,

$${\displaystyle \frac{\mathrm{d}{Z}_x}{\mathrm{d}t}}=P\left(t\right)\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)M_{j+\mathrm{2},k}\left(t\right).$$

As Z_x(t) is a product, its derivative can also be written as

$${\displaystyle \frac{\mathrm{d}{Z}_x}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{V}_x}{\mathrm{d}t}}P\left(t\right)+{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}{V}_x\left(t\right).$$

By equating the two previous expressions, we get

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{V}_x}{\mathrm{d}t}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{P\left(t\right)}}\left({\displaystyle \frac{\mathrm{d}{Z}_x}{\mathrm{d}t}}-{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}{V}_x\left(t\right)\right)\\ {\displaystyle }& {\displaystyle }=\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)M_{j+\mathrm{2},k}\left(t\right)-V_x\left(t\right)\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)M_{jk}\left(t\right).\end{array}$$

In this equation, the two terms proportional to $a\left(\stackrel{\mathrm{‾}}{x}\right(t),\stackrel{\mathrm{‾}}{y}(t),t)$ compensate, and it is no longer possible to neglect the third and fourth orders of the trait distribution. With only the two lowest-order non-zero terms retained, the time derivative of variance is

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{V}_x}{\mathrm{d}t}}=M_{\mathrm{30}}\left(t\right){\displaystyle \frac{\partial a}{\partial x}}+M_{\mathrm{21}}\left(t\right){\displaystyle \frac{\partial a}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(M_{\mathrm{40}}\right(t)-V_x^{\mathrm{2}}(t\left)\right)\\ & {\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(M_{\mathrm{22}}-V_x{V}_y){\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}+(M_{\mathrm{31}}-V_x{C}_{xy}){\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}.\end{array}\end{array}$$

The moments M₃₀(t) and M₄₀(t) are called the skewness and kurtosis of x, respectively. They represent the shape of the trait distribution. These moments could be described by their own equations, but their time derivatives depend on moments of even higher orders, and so on. In order to limit mathematical complexity and computation time, we do not explicitly compute moments of higher order than variance. Instead, we close our system by giving these moments the same value as in a bivariate normal distribution. In a bivariate normal distribution, odd-order moments (where j+k is odd) are zero (for reasons of symmetry) and even-order moments can be expressed as a function of variances (Isserlis, 1916).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{M}_{\mathrm{40}}\left(t\right)=\mathrm{3}{V}_x^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }{M}_{\mathrm{31}}\left(t\right)=\mathrm{3}{V}_x{C}_{xy}\\ {\displaystyle }& {\displaystyle }{M}_{\mathrm{22}}\left(t\right)=V_x{V}_y+\mathrm{2}{C}_{xy}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }{M}_{\mathrm{13}}\left(t\right)=\mathrm{3}{V}_y{C}_{xy}\\ {\displaystyle }& {\displaystyle }{M}_{\mathrm{04}}\left(t\right)=\mathrm{3}{V}_y^{\mathrm{2}}\end{array}$$

The equation for V_x(t) then simplifies to

$$\begin{array}{}\text{(B8)}& {\displaystyle \frac{\mathrm{d}{V}_x}{\mathrm{d}t}}=V_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}+\mathrm{2}{V}_x{C}_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}+C_{xy}^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}.\end{array}$$

The equation for V_y(t) is obtained by swapping the x and y indices:

$$\begin{array}{}\text{(B9)}& {\displaystyle \frac{\mathrm{d}{V}_y}{\mathrm{d}t}}=V_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}+\mathrm{2}{V}_y{C}_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}+C_{xy}^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}.\end{array}$$

These expressions represent the effect of competition on trait variance. If the mean trait is near an optimum, then $\frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}$ and $\frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}$ are negative and trait diversity tends to collapse due to competitive exclusion.

The equation for covariance is derived in a similar way. We define Z_xy(t) as P(t)C_xy(t). The time derivative of Z_xy(t) is as follows.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}{Z}_{xy}}{\mathrm{d}t}}=\int \int \left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)\cdot \left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)\cdot {\displaystyle \frac{\partial p}{\partial t}}(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{x}}{\mathrm{d}t}}\int \int \left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{y}}{\mathrm{d}t}}\int \int \left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)\cdot p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\end{array}$$

The last two terms are zero since M₁₀=0 and M₀₁=0. Thus,

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{Z}_{xy}}{\mathrm{d}t}}& {\displaystyle }=\int \int \left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)\cdot \left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)\cdot p(x,y,t)a(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=P\left(t\right)\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)M_{j+\mathrm{1},k+\mathrm{1}}\left(t\right).\end{array}$$

Since Z_xy is a product, its derivative can also be written as

$${\displaystyle \frac{\mathrm{d}{Z}_{xy}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{C}_{xy}}{\mathrm{d}t}}P\left(t\right)+{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}{C}_{xy}\left(t\right).$$

By equating the two previous expressions, we get

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{C}_{xy}}{\mathrm{d}t}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{P\left(t\right)}}\left({\displaystyle \frac{\mathrm{d}{Z}_{xy}}{\mathrm{d}t}}-{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}{C}_{xy}\left(t\right)\right)\\ {\displaystyle }& {\displaystyle }=\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)M_{j+\mathrm{1},k+\mathrm{1}}\left(t\right)\\ {\displaystyle }& {\displaystyle }-C_{xy}\left(t\right)\sum _{j=\mathrm{0}}^{\mathrm{\infty }}\sum _{k=\mathrm{0}}^{\mathrm{\infty }}{a}_{jk}\left(t\right)M_{jk}\left(t\right).\end{array}$$

In the case of a bivariate normal distribution, the only terms remaining produce the following equation.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{xy}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[M_{\mathrm{31}}\left(t\right)-C_{xy}\left(t\right)V_x\left(t\right)\right]{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[M_{\mathrm{13}}\left(t\right)-C_{xy}\left(t\right)V_y\left(t\right)\right]{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }+\left[M_{\mathrm{22}}\left(t\right)-C_{xy}^{\mathrm{2}}\left(t\right)\right]{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}\end{array}$$

Replacing M₃₁(t), M₂₂(t), and M₁₃(t) by their expressions as a function of lower-order moments yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}{C}_{xy}}{\mathrm{d}t}}=C_{xy}\left(t\right)V_x\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}+C_{xy}\left(t\right)V_y\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}\\ \text{(B10)}& {\displaystyle }& {\displaystyle }+\left[V_x\left(t\right)V_y\left(t\right)+C_{xy}^{\mathrm{2}}\left(t\right)\right]{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}.\end{array}$$

Competition tends to reduce variance, which must therefore be sustained by another process. In our model, this is the role of trait diffusion, as described in Appendix A and originally derived by Merico et al. (2014). Trait diffusion represents the effect of heritable mutations on the trait distribution and adds two new terms to the net growth in Eq. (B3) that are mathematically very similar to tracer diffusion.

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t}}=a(x,y,t)p(x,y,t)+{\mathit{\nu }}_x{\displaystyle \frac{{\partial }^{\mathrm{2}}\left(u\cdot p\right)}{\partial x^{\mathrm{2}}}}(x,y,t)\\ & +{\mathit{\nu }}_y{\displaystyle \frac{{\partial }^{\mathrm{2}}\left(u\cdot p\right)}{\partial y^{\mathrm{2}}}}(x,y,t)\end{array}\end{array}$$

The additional term for the total biomass time derivative due to the diffusion of x is

$$\int \int {\mathit{\nu }}_x{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}\cdot \mathrm{d}x\mathrm{d}y={\mathit{\nu }}_x\int {\left[{\displaystyle \frac{\partial (u\cdot p)}{\partial x}}\right]}_{x=-\mathrm{\infty }}^{x=+\mathrm{\infty }}\mathrm{d}y=\mathrm{0}.$$

This result comes from the fact that u⋅p and its derivatives tend towards zero when any trait tends towards infinity. The effect of the diffusion of y would equally be zero. This means trait diffusion is not a source of biomass. The additional term for $\frac{\mathrm{d}{S}_x}{\mathrm{d}t}$ due to the diffusion of x is as follows.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\int \int {\mathit{\nu }}_x\cdot x{\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}\cdot \mathrm{d}x\mathrm{d}y={\mathit{\nu }}_x\int {\left[x\cdot {\displaystyle \frac{\partial (u\cdot p)}{\partial x}}\right]}_{x=-\mathrm{\infty }}^{x=+\mathrm{\infty }}\mathrm{d}y\\ {\displaystyle }& {\displaystyle }-{\mathit{\nu }}_x\int \int {\displaystyle \frac{\partial (u\cdot p)}{\partial x}}\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }={\mathit{\nu }}_x\int {\left[x\cdot {\displaystyle \frac{\partial (u\cdot p)}{\partial x}}\right]}_{x=-\mathrm{\infty }}^{x=+\mathrm{\infty }}\mathrm{d}y-{\mathit{\nu }}_x\int {\left[u\cdot p\right]}_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\mathrm{0}\end{array}$$

The effects of y diffusion on $\stackrel{\mathrm{‾}}{x}$ and of both x and y diffusion on $\stackrel{\mathrm{‾}}{y}$ are equally zero. Trait diffusion has no effect on the mean trait equations. It does not favour any direction of change. However, trait diffusion is a source of variance. Indeed, a similar integration by parts adds a new non-zero term to $\frac{\mathrm{d}{Z}_x}{\mathrm{d}t}$.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\int \int {\mathit{\nu }}_x\cdot {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{\mathrm{2}}\cdot {\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\int {\mathit{\nu }}_x{\left[{\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{\mathrm{2}}\cdot {\displaystyle \frac{\partial (u\cdot p)}{\partial x}}\right]}_{x=-\mathrm{\infty }}^{x=+\mathrm{\infty }}\mathrm{d}y\\ {\displaystyle }& {\displaystyle }-\mathrm{2}{\mathit{\nu }}_x\int \int \left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)\cdot {\displaystyle \frac{\partial (u\cdot p)}{\partial x}}(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\mathrm{0}-\mathrm{2}{\mathit{\nu }}_x\int {\left[\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)\cdot u\cdot p\right]}_{x=-\mathrm{\infty }}^{x=+\mathrm{\infty }}\mathrm{d}y\\ {\displaystyle }& {\displaystyle }+\mathrm{2}{\mathit{\nu }}_x\int \int u(x,y,t)p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\mathrm{2}{\mathit{\nu }}_x\int \int u(x,y,t)p(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\end{array}$$

This integral is similar to $\frac{\mathrm{d}P}{\mathrm{d}t}$ (replacing a by u) and can be Taylor-expanded in the same way. The contribution of y diffusion on V_x is a different case. The extra term on $\frac{\mathrm{d}{Z}_x}{\mathrm{d}t}$ is zero.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\int \int {\mathit{\nu }}_y\cdot {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{\mathrm{2}}\cdot {\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }={\mathit{\nu }}_y\int {\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)}^{\mathrm{2}}\left[\int {\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial y^{\mathrm{2}}}}(x,y,t)\cdot \mathrm{d}y\right]\mathrm{d}x\\ {\displaystyle }& {\displaystyle }=\mathrm{0}\end{array}$$

Thus, the new equations of V_x(t) and V_y(t) accounting for trait diffusion are as follows.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}{V}_x}{\mathrm{d}t}}=V_x^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}+\mathrm{2}{V}_x{C}_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}+C_{xy}^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }+\mathrm{2}{\mathit{\nu }}_x\left[u\left(\stackrel{\mathrm{‾}}{x}\left(t\right),\stackrel{\mathrm{‾}}{y}\left(t\right),t\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_x\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}\right.\\ \text{(B12)}& {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_y\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial y^{\mathrm{2}}}}+C_{xy}\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x\partial y}}\right]{\displaystyle }& {\displaystyle \frac{\mathrm{d}{V}_y}{\mathrm{d}t}}=V_y^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial y^{\mathrm{2}}}}+\mathrm{2}{V}_y{C}_{xy}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x\partial y}}+C_{xy}^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial x^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }+\mathrm{2}{\mathit{\nu }}_y\left[u\left(\stackrel{\mathrm{‾}}{x}\left(t\right),\stackrel{\mathrm{‾}}{y}\left(t\right),t\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_x\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}\right.\\ \text{(B13)}& {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_y\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial y^{\mathrm{2}}}}+C_{xy}\left(t\right){\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x\partial y}}\right]\end{array}$$

Because $u(x,y,t)$ is positive for every value of x and y, the trait diffusion term is positive. Trait diffusion is able to counter the effect of competitive exclusion and sustain variance.

The effect on Z_xy(t) is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\int \int {\mathit{\nu }}_x\left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)\cdot \left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)\cdot {\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}(x,y,t)\cdot \mathrm{d}x\mathrm{d}y\\ {\displaystyle }& {\displaystyle }={\mathit{\nu }}_x\int \left(y-\stackrel{\mathrm{‾}}{y}\left(t\right)\right)\left[\int \left(x-\stackrel{\mathrm{‾}}{x}\left(t\right)\right)\cdot {\displaystyle \frac{{\partial }^{\mathrm{2}}(u\cdot p)}{\partial x^{\mathrm{2}}}}(x,y,t)\cdot \mathrm{d}x\right]\mathrm{d}y\\ {\displaystyle }& {\displaystyle }=\mathrm{0}.\end{array}$$

Thus, trait diffusion does not add nor remove covariance to the phytoplankton community. However, by increasing the variances, trait diffusion decreases the correlation. In other words, trait diffusion decorrelates the traits by making rare trait combinations more likely.

We note that, in the absence of trait diffusion, our equations are a particular case of the general equations derived by Bruggeman (2009) for multivariate normal trait distributions. In a single-trait case, they are simpler than the original equations (Merico et al., 2014) but identical to the more recent formulation of Coutinho et al. (2016) and derived using the same method.