PhytoSFDM version 1 . 0 . 0 : Phytoplankton Size and Functional Diversity Model

Biodiversity is one of the key mechanisms that facilitate the adaptive response of planktonic communities to a fluctuating environment. How to allow for such a flexible response in marine ecosystem models is, however, not entirely clear. One particular way is to resolve the natural complexity of phytoplankton communities by explicitly incorporating a large number of species or plankton functional types. Alternatively, models of aggregate community properties focus on 5 macroecological quantities such as total biomass, mean trait, and trait variance (or functional trait diversity), thus reducing the observed natural complexity to a few mathematical expressions. We developed the modelling tool PhytoSFDM, which can resolve species discretely and can capture aggregate community properties. The tool also provides a set of methods for treating diversity under realistic oceanographic settings. This model is coded in Python and is distributed as an open-source 10 software. PhytoSFDM is implemented in a 0D physical scheme and can be applied to any location of the global ocean. We show that aggregate-community models reduce computational complexity while preserving relevant macroecological features of phytoplankton communities. Compared to species-explicit models, aggregate models are more manageable in terms of number of equations and have faster computational times. Further developments of this tool should address the caveats 15 associated with the assumptions of aggregate community models and on implementations into spatially resolved physical settings (1D and 3D). With PhytoSFDM we embrace the idea of promoting open source software and encourage scientists to build on this modelling tool to further improve our understanding of the role that biodiversity plays in shaping marine ecosystems.


Introduction
Numerical models are simplified abstractions of complex phenomena.They are engineered for the problem at hand and cannot be designed to maximize simultaneously the three key requirements of generality, precision, and realism, because one of these must be sacrificed in favour of the other two (Levins, 1966).Marine ecosystem models are no exceptions, and the scientific community has questioned the trend towards increasing model complexity in terms of large numbers of state variables and parameters (Fulton et al., 2003;Anderson, 2005;Hood et al., 2006;Anderson, 2010).Alternatives such as trait-based models have been put forward as a way to simplify overly parameterized ecosystem models (Follows and Dutkiewicz, 2011).
In the past 2 decades, trait-based models of planktonic ecosystems have become important tools for elucidating the fundamental mechanisms behind emergent patterns of community structure and diversity.Most of these models describe the phytoplankton community by a discrete representation of many species or functional groups (Baird and Suthers, 2007;Follows et al., 2007;Bruggeman and Kooijman, 2007;Barton et al., 2010;Banas, 2011;Ward et al., 2012;Smith et al., 2015).Alternatively, models have been developed that Published by Copernicus Publications on behalf of the European Geosciences Union.E. Acevedo-Trejos et al.: PhytoSFDM treat the whole phytoplankton species assemblage as a single entity (Wirtz and Eckhardt, 1996;Norberg et al., 2001;Merico et al., 2009;Bruggeman, 2009;Wirtz, 2013;Wirtz and Sommer, 2013;Terseleer et al., 2014;Acevedo-Trejos et al., 2015).These models use aggregate community properties such as total biomass, mean trait, and trait variance to describe changes in phytoplankton community composition.Hence, by approximating the full spectrum of species or functional types with just a few macroecological properties, these models present a way of reducing the complexity of natural communities (Merico et al., 2009).
The simplification of both types of trait-based models (i.e.discrete and aggregate) relies on the use of a key trait, for which relationships with other traits can be formulated.Cell size is recognized as one of the most important traits for characterizing phytoplankton communities (Litchman and Klausmeier, 2008;Litchman et al., 2010;Finkel et al., 2010;Marañón, 2015), and it has been commonly used in plankton ecosystem models (Baird and Suthers, 2007;Banas, 2011;Ward et al., 2012;Wirtz, 2013;Wirtz and Sommer, 2013;Terseleer et al., 2014;Acevedo-Trejos et al., 2015;Smith et al., 2015).This morphological trait affects trophic organization of foodwebs and the sequestration of CO 2 into the ocean interior (Chisholm, 1992).Phytoplankton size also impacts on many ecological and physiological functions and is linked to other relevant traits via trade-off relationships (see reviews by Litchman and Klausmeier, 2008;Litchman et al., 2010;Finkel et al., 2010).Therefore, studies on how cell size is associated with ecological and physiological processes and on the impact that these associations have on the structure and functioning of planktonic communities are of fundamental importance (Marañón, 2015;Andersen et al., 2015).
Here we present the new Phytoplankton Size and Functional Diversity Model (called PhytoSFDM) that allows for five different ways of describing the size composition of phytoplankton communities in the upper mixed layers of the world oceans.In the first variant, the phytoplankton community is described according to the classical approach that resolves the discrete assemblage of many different species and then we present four alternative ways of expressing aggregate community properties of phytoplankton based on four different ways of treating size diversity.We provide this model as open-source so that it can be used, modified, and redistributed freely with the aims of fostering reproducibility and encouraging investigations about the impact of environmental conditions on properties of phytoplankton community structure and diversity.

Model description
PhytoSFDM is developed from the study of Acevedo-Trejos et al. (2015), which used a size-based model of aggregate community properties to investigate the phytoplankton size structure and size diversity in two environmentally contrast-ing regions of the Atlantic Ocean.In this model, the phytoplankton community self-assembles according to a trade-off emerging from relationships between cell size and (1) nitrogen uptake, (2) zooplankton grazing, and (3) phytoplankton sinking.In PhytoSFDM we have extended this work by providing four ways of treating size diversity using a momentbased approximation (see Smith et al., 2011;Bonachela et al., 2015, and Sect. 2.1.3 in this study).In addition, we include a discrete version of the model (hereafter referred to as the full model) to better illustrate the potential of using aggregate models as compared to the equivalent discrete version.
In the following, we present the mathematical equations, a description of the code structure, and easy-to-follow examples of how to use the model.

Mixed-layer scheme
The zero-dimensional physical set-up consists of two vertical layers, the upper mixed layer containing the pelagic ecosystem and the abiotic bottom layer with nitrogen concentration as forcing.Following Evans and Parslow (1985) and Fasham et al. (1990), we describe material exchange between the two layers (K) as a function of the mixed-layer depth (M), where κ is a constant that parameterizes diffusive mixing across the thermocline and h + (t) is a function that describes entrainment and detrainment of material.The latter is given by h + (t) = max[h(t), 0], with h(t) = dM(t)/dt.Zooplankton are considered capable of maintaining themselves within the upper mixed layer; thus, their mixing term simplifies to K Z = h(t).

Dynamics of the full phytoplankton community
The description of the phytoplankton community is a trait-based variant of the classical nutrient-phytoplanktonzooplankton-detritus (NPZD) model (Fasham et al., 1990).We consider only one nutrient, nitrogen, which constitutes the currency of our model, one zooplankton population (composed of individuals assumed to be identical), and a single detritus pool.We define n morphologically distinct phytoplankton types (hereafter referred to as morphotypes), and we consider n equal to either 10 or 100.Each morphotype is characterized by a biomass P i and a cell size S i , in units of µm equivalent spherical diameter (ESD).The distribution of biomass along the size dimension is known to be positively skewed (i.e. an asymmetrical size distribution with a pronounced right tail compared to its left tail), due to physiological, morphological, and ecological constraints that limit phytoplankton from a minimum size of around 0.15 µm ESD to a maximum size of about 575 µm ESD (Marañón, 2015;Andersen et al., 2015).Consequently, we assume a log-normal distribution of size to represent the size of each morphotype, thus transforming the cell size S i as follows: L i = ln(S i ).The net growth rate of the whole phytoplankton community (P) is then given by dP dt where f i (L i , E) is the net growth rate of size class i, which we assume to be a proxy for fitness (Smith et al., 2011).Hence f i accounts for the gains and losses of each morphotype as a function of cell size L i and environment E. The latter includes changes in nitrogen, irradiance, temperature, and grazing.The equation describing the fitness functions of each size class i is thus given by where µ P indicates the maximum growth rate and F (T ) = e 0.063•T is Eppley's formulation for temperature-dependent growth (Eppley, 1972).The light-limiting term, H (I ), represents the total light I available in the upper mixed layer.
According to Steele's formulation (Steele, 1962), where I s is the light level at which photosynthesis saturates and I (z) is the irradiance at depth z.The exponential decay of light with depth is computed according to the Beer-Lambert law with a generic extinction coefficient k w : The current version of our model does not specify any size dependence for light absorption, although we provided suggestions on how this could be done (Sects.4 and 6).
The nutrient-limiting term U in Eq. ( 3) is determined by a Monod function with a half-saturation constant K N , which scales allometrically with phytoplankton cell size L (Litchman et al., 2007), with β U and α U , respectively, intercept and slope of the K N allometric function (i.e. the power law β U •S α U ).This empirical relationship is based on observations of different phytoplankton groups (see Fig. 3b in Litchman et al., 2007), with the regression parameters rescaled from cell volume to ESD.The loss term G(L i , P i ) in Eq. (3) represents zooplankton grazing.As mentioned above, here we consider a single zooplankton population, which is assumed to be an assemblage of identical individuals with a size-selective feeding preference given by where α G is the slope for size-dependent grazing (or the power law 1 • S α G ) and K P is the half-saturation constant.This formulation is inspired by meta-analyses of laboratory data (Hansen et al., 1994(Hansen et al., , 1997) ) and reflects a grazing preference of zooplankton for smaller phytoplankton cells.For demonstration purposes, we use here a simple formulation for zooplankton grazing; however, other functional relationships can be implemented and tested in future versions of PhytoSFDM (see also Sects. 4 and 6).
The loss term V (L i , M) in Eq. ( 3) represents the sinking of phytoplankton as a function of size and depth of the mixed layer, where the constants α V and β V are the parameters of the function relating phytoplankton cell size to sinking velocity according to Stokes' law (Kiørboe, 1993) or the power law β V •S α V .These parameters are expressed here in units of metres per day.Our model formulation does not specify an explicit size dependence for the phytoplankton maximum growth rate (µ P ).Various compilations of data from laboratory experiments reveal different size scalings for µ P , either as a power law of cell volume (Litchman et al., 2007;Edwards et al., 2012) or as a unimodal distribution in terms of cell size (Wirtz, 2011;Ward et al., 2012;Marañón et al., 2013).Therefore, we adopted an approach similar to that of Smith et al. (2015), who reproduced the unimodal distribution of realized growth rate over size using two physiological trade-offs.We specified our trade-off in terms of three allometric relationships, and this results in an indirect size dependence of phytoplankton growth rate.
The loss term m P in Eq. ( 3) accounts for all phytoplankton losses other than those from grazing and mixing.
Differential equations for the nutrient (N), zooplankton (Z), and detritus (D) complete the model system: www.geosci-model-dev.net/9/4071/2016/Geosci.Model Dev., 9, 4071-4085, 2016 where δ D is the mineralization rate and N 0 is the nitrogen concentration below the upper mixed layer.µ Z , δ Z , and m Z are, respectively, maximum growth rate, prey assimilation coefficient, and mortality rate of zooplankton.All parameter values and their units are reported in Table 1.

Dynamics of the aggregate phytoplankton community
The phytoplankton community comprising many distinct morphotypes (Eqs. 2 to 8) can be approximated with the socalled moment-based approach (Wirtz and Eckhardt, 1996;Norberg et al., 2001;Merico et al., 2009;Terseleer et al., 2014;Acevedo-Trejos et al., 2015).Wirtz and Eckhardt (1996), Norberg et al. (2001), andMerico et al. (2009) used a Taylor expansion together with a moment closure technique to approximate the whole community with three macroscopic properties, which correspond to the first three order moments of the approximated biomass distribution.These properties are total biomass, mean trait, and trait variance.These works (Wirtz and Eckhardt, 1996;Norberg et al., 2001;Merico et al., 2009) were inspired by earlier applications in quantitative genetics (Abrams et al., 1993) and have been reviewed by Smith et al. (2011) and more recently by Bonachela et al. (2015).
Here the whole phytoplankton community is characterized by the morphological trait cell size and by a trade-off that emerges from three allometric relationships described by Eqs. ( 6)-( 8).The equations of the respective macroscopic properties are where f is the net growth rate (or the fitness function; see Eq. 3) and f (n) is the nth derivative of the net growth with respect to the trait.Due to competitive exclusion, however, the phytoplankton community loses functional diversity over time, i.e. the variance declines to zero with time, in both full and aggregate model formulations (Merico et al., 2014).We name this standard formulation "unsustained variance".Alternatively, one can use the approximated model to focus only on changes in the mean trait, thus ignoring changes in the variance by fixing it to an arbitrary constant value: dV dt = 0. (18) While using these two formulations (i.e.unsustained and fixed variance) can be acceptable in some special cases (e.g. in experiments that lead to competitive exclusion or where diversity is being manipulated), it is clear that they fail to account for changes in the adaptive capacity of the community, which requires allowing the size variance, and thereby functional diversity, to vary over time (Merico et al., 2014).
Within our modelling tool we also provide two alternative ways of treating the size variance: immigration (following Norberg et al., 2001) and trait diffusion (following Merico et al., 2014).The treatment with immigration considers the introduction of biomass and new trait values from hypothetical adjacent communities into the resident community.The addition of incoming amount of biomass per day is named immigration I , where L I and V I are, respectively, the mean size and the size variance of the immigrating community.As implemented by Acevedo-Trejos et al. (2015), we treat I as a densitydependent process (i.e.I = δ I • P), and set L I equal to the mean size of the resident community (i.e.L I = L).Thus, we assume that phytoplankton immigrating from adjacent areas are characterized by sizes similar to the simulated community, implying that the immigrating community has been exposed to the same selection pressures as the simulated community (Terseleer et al., 2014).We also assume that the rate of immigration increases proportionally to the concentration of phytoplankton, consistent with observations of diversity patterns along the Atlantic Ocean (Chust et al., 2013).
The treatment of the size variance based on trait diffusion (Merico et al., 2014) gives where ν is the trait diffusivity parameter, r is the reproduction rate (or gross growth), and r n is the nth derivative of gross growth with respect to the trait.Note that the process of trait diffusion (last term in Eq. 24) depends on the gross growth r, via the trait diffusivity constant ν; thus, an increase in phy-toplankton gross growth causes an increase in trait variance (Merico et al., 2014).
The system of differential equations for all variance treatments is completed by equations describing gains and losses in nitrogen (N), zooplankton (Z), and detritus (D): The first term in Eq. ( 25) represents a reduction of the nitrogen pool due to phytoplankton growth, which is a function of temperature, light, nitrogen, and mean size (see the description of Eq. 3 in the previous section).The last two terms in Eq. ( 25) represent sources of nitrogen due to remineralization and mixing.The first term in Eq. ( 26) describes sizedependent grazing, while the last two terms describe losses of zooplankton due to mortality and mixing.The first term in Eq. ( 27) represents a fraction of phytoplankton biomass that is not assimilated by zooplankton and the following two terms represent the mortality of phytoplankton and zooplankton, respectively.The detritus pool is reduced by remineralization and mixing.Parameter values and their units are reported in Table 1.

Environmental forcing
We compiled monthly climatological forcing data for mixedlayer depth (MLD), photosynthetic active radiation (PAR), sea surface temperature (SST), and concentration of nitrogen immediately below the upper mixed layer (N 0 ).The MLD data were obtained from Monterey and Levitus (1997) using the variable density criterion and are openly accessible from https://www.nodc.noaa.gov/OC5/WOA94/mix.html.The PAR data were obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS), for the time period 2002-2011.This dataset is managed and distributed by the NASA's Ocean Biology Processing Group (http:// oceancolor.gsfc.nasa.gov/cms/).SST and N 0 were obtained from the World Ocean Atlas 2009 (WOA09), which is maintained and distributed by NOAA (https://www.nodc.noaa.gov/OC5/WOA09/pr_woa09.html).For consistency and efficiency, all data were transformed from their original formats (e.g.TXT and HDF) to NetCDF.All monthly forcings were spatially averaged over the selected location (square boxes in Fig. 1) and then interpolated to obtain daily values (Fig. 2).

Test-case simulation
A test-case model configuration is provided for a location of the North Atlantic Ocean at 47.5 • N 15.5 • W (Fig. 1), a region where seasonal changes in mean size and size diversity are well known (Acevedo-Trejos et al., 2015).This region presents the typical oceanographic conditions of a temperate environment (Fig. 2).The environmental conditions produce a pronounced phytoplankton bloom in spring, which stimulates secondary production and almost the full depletion of nitrogen (Fig. 3).Overturning of the water column in autumn restocks the pool of nitrogen and light limitation together with lower temperatures halt primary production (Figs. 2 and  3).

Comparison of full and aggregate models
Within PhytoSFDM, we provide a practical example of how to implement and compare phytoplankton community models that aim to describe (a) a full assemblage of species or morphotypes (see Sect. 2.1.2),and (b) an aggregate community (see Sect. 2.1.3).The aggregate community model is an approximation of the full assemblage of species or morphotypes (Wirtz and Eckhardt, 1996;Norberg et al., 2001;Merico et al., 2009).Figures 3 and 4 show the results of, respectively, the full model and the aggregate model for the unsustained variance case.N, P, Z, and D are unaffected by the type of model considered.As expected, the dynamics of P, L, and V produced by the aggregate model are good approximations of those produced by the full model.Both models exhibit competitive exclusion, as indicated by the reduction in the number of morphotypes and consequently in the loss of size variance over time (Fig. 4).The phytoplankton community evolves towards the optimal trait value, which is expressed by the fittest few morphotypes for the chosen parameterization and the prevailing environmental conditions.Although competitive exclusion is well established theoretically (Hardin, 1960), natural communities of phytoplankton are typically very diverse; hence, we will explore in the following the effects of different ways of sustaining the variance.

Comparison of variance treatments
The key aspect of trait-based models is their ability to describe the phytoplankton community in terms of mean trait and trait variance.Figures 5 and 6 show the results of 1-year simulation after an initial spin-up phase of 4 years.While the four treatments produce very similar, if not identical, dynamics for N, P, Z, and D (Fig. 5), the results for the mean size and the size variance differ considerably among treatments (Fig. 6).
As already discussed, the system loses diversity over time when variance is unsustained.The loss of diversity reduces the capacity of the community to adapt to changing environmental conditions via shifts in species composition, as a flat year-round mean trait shows (Fig. 6, grey lines).Under fixed variance, size diversity is locked at an arbitrary value.If this value is high enough, the mean size can adapt in response to changes in nutrient availability and grazing regimes (Fig. 6).This treatment can be useful for studies focusing only on the size structure of the community, but it is otherwise based on an arbitrarily fixed level of diversity and cannot offer meaningful insights, for example about biodiversity and ecosystem functioning relationships.
Trait diffusion and immigration show similar results for the mean size but not for the size variance (Fig. 6).Since the mechanism of trait diffusion depends on reproduction, i.e. gross growth (see Eq. 24), the highest diversity of the community is reached in spring under high growth rates and declines when moving towards winter.Size diversity also peaks in spring for the case of immigration because this mechanism is density-dependent (see Eq. 21), but the variances predicted in autumn and winter are, respectively, lower and higher than those obtained with trait diffusion (Fig. 6).As mentioned above, this originates from the different assumptions underlying the trait diffusion and immigration treatments, which consider, respectively, an internal or external source of phytoplankton biomass, mean trait, and trait variance.In the case of trait diffusion, such an internal source is gross growth because the size variance of the phytoplankton community is proportional to it via the diffusivity constant ν (last terms in Eqs. 22,23,and 24).In contrast, immigration represents a source of biomass (I ) and size variance (I /P[V I − V ]) external to the phytoplankton community being simulated (e.g. from an adjacent patch).Hence, during the autumn-winter transition, the size variance tends to decline in the trait dif-  fusion case as phytoplankton gross growth is reduced by growth-limiting processes.Instead, the trait variance keeps building up to values similar to the variance of the immigrating community in the case of immigration.

Sensitivity to changes in parameter values
We tested the sensitivity of the annual mean in P, L, and V to variations of ±25 % in parameter values.To quantify this sensitivity, we formulated an index S that accounts for relative changes in model results: where X(p) is the result of the state variable X obtained with the standard parameter p and X(p ) is the result of the state variable X obtained with the modified parameter p = p ± 25 %.
The four treatments of size variance respond similarly to changes in parameter values (Fig. 7).The annual means of all three state variables (P, L, and V ) are sensitive to changes in the parameters controlling zooplankton grazing (i.e.µ Z , m Z , K P , δ Z ).However, P also shows a sensitive response to parameters affecting phytoplankton gross growth, such as k w , I s , µ P , and m P .Mean size is the most robust variable, with less than 10 % relative change compared to the standard run.The size variance treatments for immigration and trait diffusion are affected by the parameters controlling the input of exogenous (i.e.δ I for immigration) or endogenous variance (i.e.ν for trait diffusion).The results of the unsustained variance model are very sensitive to changes in µ Z , and the case of fixed variance shows a sensitivity that is similar to the other cases, except for the variance itself.

Computational efficiency
Trait-based models that aim at resolving the complexity of natural communities by incorporating many different species or functional types can be expensive in terms of computational time (Baird and Suthers, 2007;Follows et al., 2007;Bruggeman and Kooijman, 2007;Banas, 2011;Ward et al., 2012).Alternatively, trait-based models that focus on aggregate community properties such as total biomass, mean trait, and trait variance can be more computationally efficient.In Table 2 we report a quantification of the computation time required for calculating the full and aggregate models presented here.We obtained a more than 10-fold longer computation time for the full model than for the aggregate model.In addition, when we increase the resolution of the full model from 10 to 100 morphotypes, the difference in computation time increases by more than 20-fold.Thus, increasing the realism in terms of number of species or morphotypes comes at a significant computational cost.

Strength and weakness of moment-based approximations
Models are simplifications of reality and, as such, are based on assumptions.For example, the simple exponential growth model is based on a number of assumptions that do not hold in all circumstances (many factors affect the intrinsic growth rate, which is often not time-invariant, not all individuals within a population are identical, nothing can grow indefinitely, etc.).However, this model is widely used within its range of validity.Likewise, the approximation of full models with moment-based approaches requires an assumption about the shape of the phytoplankton trait distribution (Wirtz     2) is the second derivative of the fitness function with respect to the trait; f (2) U (L, N), f (2) G(L, P), and f (2) V (L, M) are, respectively, nitrogen uptake, zooplankton grazing, and phytoplankton sinking components of f (2) .
and Eckhardt, 1996;Norberg et al., 2001;Merico et al., 2009).Typically, unimodal distributions, e.g.normal or lognormal, are assumed.However, depending on how the fitness function (i.e. the net growth rate of the phytoplankton community) is constructed and parameterized, the value of f (2) , that is the rate of change of the variance (Eqs.15, 18, 21, and 24), can be positive, implying a disruptive selection or branching.This represents an indication that the unimodality assumption does not hold (Bonachela et al., 2015).Alternatively, f (2) can remain negative over time, implying that the community continually loses variance, thus constituting a strong indication against the occurrence of disruptive selection.Therefore, models based on moment approximations require careful checks about the validity of the unimodality assumption throughout the time of the simulations.Figure 8 shows, for our test case, the predicted variance V , f (2) , and its components for the four variance treatments.In our test case, f ( 2) is negative for all treatments and its changes are jointly driven by bottom-up, f (2) U (L, N), and top-down processes, f (2) G(L, P), i.e. the second derivatives with respect to the trait for nitrogen uptake (Eq.6) and grazing (Eq.7) terms.Sinking plays a role mainly during spring, but its influence is minor compared to the effects of nitrogen uptake and grazing.
It is unclear whether unimodality in size distributions is a robust feature in the oceans.Observational evidence from recent work (Downing et al., 2014) suggests that at large temporal scales, from 5 to 20 years, unimodality of size distributions is a consistent feature of phytoplankton communities of the North Sea.By contrast, multimodality is typically observed on temporal scales of less than 1 year (Downing et al., 2014).We consider that the observational evidence available remains insufficient to identify general patterns.However, the ocean is a highly variable environment and we considered it more likely that multimodality, for example because of size-selective grazing events, is a short-term, transient situation rather than the norm, because mixing would continuously reshuffle plankton assemblages and restore homogeneous conditions.
An aspect that our model does not include in its current version is the dependency of light acquisition on phytoplankton cell size.Given that the effect of cell size on light harvesting is well understood (Augusti, 1991;Finkel and Irwin, 2000;Finkel, 2001), it could be implemented in the model.Future versions of PhytoSFDM could address this gap by considering the vertical attenuation of light as a function of both phytoplankton biomass and cell size, following the approach proposed by Baird and Suthers (2007).
Uncertainty remains about how to describe the zooplankton population, which we simplified as an assemblage of identical individuals.This has been the standard approach in plankton ecosystem modelling for decades and we based the first version of PhytoSFDM on this simple and classical formulation.In recent years, however, significant efforts have been made to increase the level of detail of the zooplankton component in ecosystem models.Approaches are numerous and include the consideration of different zooplankton functional types, different size classes, and different feeding prefwww.geosci-model-dev.net/9/4071/2016/Geosci.Model Dev., 9, 4071-4085, 2016 erences and strategies (Banas, 2011;Ward et al., 2012;Prowe et al., 2012;Wirtz, 2012;Mariani et al., 2013;Vallina et al., 2014;Ryabov et al., 2015).A trait-based description of zooplankton can help in reducing model complexity while maintaining an adequate representation of diversity.The selection of traits to consider for ecosystem models will depend on the questions under scrutiny.For example, traits that could characterize zooplankton-related processes in ecosystem models that focus on nutrient cycling are maximum growth rates, stoichiometric requirements, and the size distribution of food particles (Litchman et al., 2013).Since many zooplankton traits scale allometrically with body size, scaling laws should be considered because they are effective ways to generalize the relationships among different traits and thus to reduce model complexity (i.e.add size-related functionality without the need for discretely parameterized zooplankton classes).
Implementing such a diversity of grazing mechanisms and processes is a natural step forward in the development of ecosystem models.However, a consistent representation of different grazing strategies remains an aspect under development (Litchman et al., 2013;Smith et al., 2014).PhytoSFDM constitutes a starting model platform for gradually building model complexity at different trophic levels.

Concluding remarks
Biological communities are complex adaptive systems (Levin, 1998) characterized by many components and interconnections that lead to emergent properties and non-linear responses.Models help us to formalize and simplify the complexity we observe in nature.This simplification allows us to render natural phenomena treatable and testable (Levins, 1966;Anderson, 2005Anderson, , 2010)).Over time, however, phytoplankton models have grown more complex, computationally more complicated, and often inaccessible to the wider scientific community, aspects that can all hamper advancements in the field.To help reverse this trend we developed PhytoSFDM as a tool to promote the use of trait-based models (whether species-explicit or aggregate models) of marine ecosystems.
A key decision in modelling is choosing an appropriate level of detail for the problem at hand.For example, a species-explicit model offers obvious advantages, which aggregate models cannot offer, when the interest lies in understanding the relative importance of particular species in providing certain ecological services or in quantifying the effect of disruptive selection.Aggregate models, instead, can be more useful at a higher level of abstraction, when the interest lies in macroecological properties.In addition, as we have shown, aggregate models present an advantage with respect to computation time when compared to full models.The advantages in terms of reducing complexity and computation time remain unproven in spatially explicit settings (e.g. in 1-D and 3-D), although preliminary applications have shown promising results (Bruggeman, 2009).
PhytoSFDM provides a set of methods, under the opensource concept, to quantify macroecological properties of phytoplankton communities, as an alternative to the traditional discrete, species-explicit approach.This effort, we hope, will foster our understanding about the role that biodiversity plays in shaping marine ecosystems.

Code availability
PhytoSFDM is written in Python (version 2.7.x) as a lightweight and user-friendly package to facilitate use and redistribution.We provide PhytoSFDM as free software under GNU General Public License version 2. The python package is hosted in (a) GitHUB (https://github.com/SEGGroup/PhytoSFDM), a software repository that allows for version control, (b) Zenodo (https://zenodo.org/record/49849),an open scientific repository, and (c) PyPI (https://pypi.python.org/pypi/PhytoSFDM), one of the most popular Python package repositories.To be able to install and operate the package, the user should be familiar with the Python language and should have a running Python distribution (preferably version 2.7.x) that includes the latest versions of the pip and setuptools libraries.Additional required libraries are matplotlib, numpy, scipy, and sympy.PhytoSFDM can then be conveniently installed by typing the following command from a terminal window: $ pip install PhytoSFDM or by downloading the tarball from the GitHub repository.This is installed using the source file setup.pycontained in the PhytoSFDM folder by typing

$ python setup.py install
The package consists of three main modules: Example, SizeModels, and EnvForcing.Example is the entry point: it computes and compares full and aggregate models with the four treatments of variance (unsustained, fixed, trait diffusion, and immigration) at the testing location in the North Atlantic Ocean (centred at 47.5 • N and 15.5 • W).The example is run from a terminal by typing $ PhytoSFDM_example or from an interactive python shell by typing >>> import phytosfdm.Example.exampleas exmp >>> exmp.main() The module SizeModels contains the model variants.Here the user can (a) modify the default parameters, (b) symbolically solve the derivatives with respect to the trait, and (c) log-transform mean trait and trait variance.To run the model at a specific location in an interactive Python shell, one should type >>> from phytosfdm.SizeModels.sizemodelsimport SM >>> Lat=47.5 >>> Lon=344.5 >>> RBB=2.5 >>> SM1= SM(Lat,Lon,RBB,"Imm") In the above example, the model is executed at a location in the North Atlantic Ocean centred at 47.5 • N and 15.5 • W (here transformed to a scale of 0 to 360 • ).RBB (range of the bounding box) specify the range of the bounding box (in degrees) for averaging the environmental forcing variables.The fourth argument SM1 is an object that contains the call of the function SM, which runs the size model at the specified location and with the desired treatment for the size variance, in this case immigration.
The last module, EnvForcing, consists of a class containing spatially averaged forcing data.The climatological data have monthly resolution, but we include a method to interpolate the data to a daily time step.Spatially averaged and temporally interpolated forcing at a specific location can be extracted by typing >>> MLD=ExtractEnvFor(Lat,Lon,RBB,'mld') Additional information on the usage of the package is contained in the Readme file and in the repository webpage in GitHUB.The source code of our model is fully and freely accessible.Users can modify or add new model variants.This can be done by manipulating the SizeModels module, which contains model variants as separated methods within the class SM.By using a version control system such as GitHUB, users can fork our repository, i.e. create a copy, which allows one to freely change and experiment without affecting the original code.Users can also modify the original code and submit a new version by pulling a request.More details can be found in our GitHUB repository (https://github.com/SEGGroup/PhytoSFDM;Acevedo-Trejos et al., 2016).

Figure 1 .
Figure 1.Environmental forcing variables considered in Phy-toSFDM.The data shown are the annual average of mixed-layer depth (MLD), photosynthetic active radiation (PAR), sea surface temperature (SST), and nitrogen concentration below the mixed layer (N 0 ).The square boxes mark the location of the test-case simulation.

Figure 2 .
Figure 2. Temporal variation of the environmental variables.The monthly climatology data (red dots) are spatially averaged over the test location (square boxes in Fig. 1).The interpolation (continuous line) is obtained with a third-(MLD and PAR) and a fifth-order (SST and N 0 ) polynomial.

Figure 3 .
Figure 3. NPZD dynamics of the full model (Sect.2.1.2) and of its equivalent aggregate model (Unsustained variance, Sect.2.1.3)for the last year of the simulations.The total phytoplankton in the full model corresponds to the sum of all P n .The red dots are observations of nitrogen concentrations (monthly data obtained from the World Ocean Atlas) and the green dots are remotely sensed Chl-a data (8-day composite obtained from MODIS).

Figure 4 .
Figure 4. Number of morphotypes and size variance over the first year of the simulation.Here we included the morphotypes with a biomass greater than 0.01 mmol N m −3 .Models that do not consider a mechanism to sustain variance exhibit competitive exclusion of morphotypes and a rapid decline of size diversity.

Figure 5 .
Figure 5. Nutrient, phytoplankton, zooplankton, and detritus dynamics over a seasonal cycle for the four variants of the aggregate model (see Sect. 2.1.3),named unsustained and fixed variance, trait diffusion, and immigration.

Figure 6 .
Figure 6.Dynamics of the size-structured phytoplankton community and its functional size diversity for the four variance treatments (see Sect. 2.1.3),named unsustained and fixed variance, trait diffusion, and immigration.

Figure 7 .Figure 8 .
Figure 7. Sensitivity of four variance treatments to an increase and a decrease by 25 % in the default parameter values.The values and definitions of all parameters are given inTable 1.

Table 1 .
Parameter definitions, their units, and their default values as provided in PhytoSFDM.

Table 2 .
Computation time in seconds for the full model with 10 and 100 morphotypes and the four variants of the aggregate model.