CITRATE 1 . 0 : Phytoplankton continuous trait-distribution model with one-dimensional physical transport applied to the North Pacific

Diversity plays critical roles in ecosystem functioning, but it remains challenging to model phytoplankton diversity in order to better understand those roles and reproduce consistently observed diversity patterns in the ocean. In contrast to the typical approach of resolving distinct species or functional groups, we present a ContInuous TRAiT-basEd phytoplankton model (CITRATE) that focuses on macroscopic system properties such as total biomass, mean trait values, and trait variance. This phytoplankton component is embedded within a nitrogen–phytoplankton-zooplankton– detritus–iron model that itself is coupled with a simplified one-dimensional ocean model. Size is used as the master trait for phytoplankton. CITRATE also incorporates “trait diffusion” for sustaining diversity and simple representations of physiological acclimation, i.e., flexible chlorophyllto-carbon and nitrogen-to-carbon ratios. We have implemented CITRATE at two contrasting stations in the North Pacific where several years of observational data are available. The model is driven by physical forcing including vertical eddy diffusivity imported from three-dimensional general ocean circulation models (GCMs). One common set of model parameters for the two stations is optimized using the Delayed-Rejection Adaptive Metropolis–Hasting Monte Carlo (DRAM) algorithm. The model faithfully reproduces most of the observed patterns and gives robust predictions on phytoplankton mean size and size diversity. CITRATE is suitable for applications in GCMs and constitutes a prototype upon which more sophisticated continuous trait-based models can be developed.


Introduction
Phytoplankton are a polyphyletic group of oxygenic organisms that account for nearly half of the global primary production (Field et al., 1998) and also play indispensable roles in other biogeochemical cycles in the Earth system (Falkowski, 2012).They have astonishingly high diversity, with several thousand species already documented and many remaining to be explored (Sournia et al., 1991;Moon-van der Staay et al., 2001).Their equivalent spherical diameter (ESD) can range from less than 1 µm for cyanobacteria such as Prochlorococcus (Chisholm et al., 1988) to more than 1 mm for some giant diatoms (Villareal, 1993).Furthermore, physiology differs substantially even within the same genera or species and the role of intraspecific variability in population dynamics and biogeochemical cycles remains to be investigated (Strzepek and Harrison, 2004;Johnson et al., 2006;Palenik et al., 2006;Kooistra et al., 2008;Biller et al., 2015).The roles of phytoplankton diversity in marine ecosystem functioning have not been understood as thoroughly as those of plant diversity in terrestrial ecosystems (Tilman et al., 1997(Tilman et al., , 2014)).
Although various ocean models have been developed by accounting for different functional groups or categories of phytoplankton (e.g., Le Quéré et al., 2005;Hashioka et al., 2013), the finite number of such distinct types included limits their ability to resolve the vast diversity of trait values.Some pioneering studies have considered greater numbers of species, each of which is defined by a particular set of multivariate trait axes that constitute a hyper-volume niche space (Follows et al., 2007;Barton et al., 2010;Follows and Dutkiewicz, 2011;Masuda et al., 2017).It is worth noting that these diversity models usually focus on "functional traits", which are the key to linking phytoplankton diversity, environmental conditions, and ecosystem functioning.Important phytoplankton traits include maximal growth rate, the light absorption and nutrient uptake affinities, optimal growth temperature, and edibility (i.e., susceptibility to grazing; Litchman et al., 2007;Litchman and Klausmeier, 2008;Edwards et al., 2011Edwards et al., , 2012Edwards et al., , 2015;;Merico et al., 2009;Barton et al., 2010, Thomas et al., 2012;Chen, 2015).The total species pool in these modeling studies should ideally cover the entire multi-dimensional trait space constrained by tradeoffs (Barton et al., 2010;Smith et al., 2011), although computational limits make this impossible in practice.As a compromise, only a limited set of trait combinations is sampled from the entire trait space.Although this approach has effectively generated large-scale patterns of plankton diversity, it generally underestimates local diversity for two reasons: (1) lack of appropriate mechanisms for sustaining diversity (but see Vallina et al., 2014) and ( 2) insufficient trait resolution so that fitness differences between species are too large to allow coexistence (i.e., insufficient equalizing effect; see Chesson, 2000).In any case, a substantial proportion of the idealized species so modeled cannot survive under realistic oceanic conditions, and therefore the models do not capture the functions associated with many species.
Continuous trait-based models have been developed to address the above questions (Wirtz and Eckhardt, 1996;Norberg et al., 2001;Bruggeman, 2009;Merico et al., 2009Merico et al., , 2014;;Terseleer et al., 2014;Acevedo-Trejos et al., 2015, 2016;Smith et al., 2016).Instead of modeling the dynamics of individual species, continuous trait-based models or socalled "adaptive dynamics" models focus on macroscopic or aggregate properties of a community such as total biomass, average trait, and trait variance by assuming that phytoplankton traits follow some distribution (usually Gaussian; Smith et al., 2011).These models do not have the problem of inadequate trait resolution because they have infinitesimally fine trait resolution.The trait variance, treated as a tracer in the model, serves as a measure of trait diversity.Thus, the continuous trait-based model has the advantage that the factors controlling diversity can be directly quantified and better understood because the sources (e.g., speciation or immigration) and sinks (e.g., resource competition) for diversity are specified explicitly.Although the size variance cannot be simply equated to species richness, it can be converted to other diversity metrics such as the continuous entropy (Quintana et al., 2008).Moreover, the diversity of functional traits is arguably a better diversity index than species richness relating to ecosystem functioning (Loreau et al., 2001).In addition, these models are computationally much more efficient than classic discrete species approaches.For example, assuming two independent traits for the phytoplankton community, a continuous trait-based model only requires 1 (biomass) + 2 × 2 (trait mean and variance) = 5 tracers for the phytoplankton community, while a discrete species-based model requires 2 × 10 = 20 tracers if assuming 10 discrete values in each trait dimension, which still provides only coarse trait resolution.Furthermore, this difference increases linearly with trait dimension.
Relatively few continuous trait-based models have been coupled with physics transport and calibrated against oceanic observations.Here we describe a new one-dimensional model, CITRATE 1.0, built upon the classic nitrogenphytoplankton-zooplankton-detritus (NPZD) model with a phytoplankton community represented using a continuous distribution of size taken as a master trait (Fig. 1).In this way, not only total phytoplankton biomass, but also phytoplankton mean size and size variance, are explicitly modeled.The distributions of other important functional traits are implicitly modeled via well-established scaling power laws.Although this approach might overlook some other important traits that are not related to size and thereby underestimate trait diversity to some extent, it serves as a starting point for later development of more comprehensive diversity models that can include more traits or be integrated with the discrete functional group approach.For the model to be implemented in the subarctic North Pacific, a well-known high-nitrate lowchlorophyll (HNLC) region, CITRATE also incorporates an iron limitation module.We optimized the model parameters against the extensive observational data at two contrasting stations (K2: 160 • E, 47 • N; S1: 145 • E, 30 • N) in the North Pacific (Fig. 2a).The station K2 is located within the western subarctic North Pacific gyre and is characterized by low temperature, high nitrate, and high carbon export (Matsumoto et al., 2014;Wakita et al., 2016).Iron limitation on phytoplankton growth has been suggested at this station (Fujiki et al., 2014).The station S1 is located within the western subtropical North Pacific and is characterized by high sea surface temperature, low levels of nitrate, and carbon export ef- ficiency (Matsumoto et al., 2016;Sasai et al., 2016;Wakita et al., 2016).To independently validate the model, we also use the optimized model parameters from stations K2 and S1 to run the model for station ALOHA (158 • W, 22.75 • N) and compare the model outputs with the extensive observational data collected there.
In the following sections, we first describe the details of the model structure and the parameter optimization subroutine.Then we show the results of parameter optimization and modeled patterns of nutrients, phytoplankton biomass, mean size, and size diversity.We also discuss the merits and limitations of the model and of the continuous trait-distribution approach.CITRATE is intended as a prototype for later incor-poration into three-dimensional (3-D) general ocean circulation models (GCMs) and for further development of more comprehensive trait-based models.

Model description
The aim of the present study is to design and implement a continuous trait-based model (CITRATE 1.0) at two representative stations in the North Pacific.The overall goal of this model is not only to simulate the phytoplankton size diversity but also to faithfully reproduce the seasonal and vertical dynamics of other important quantities, such as nu- trients, Chl a, and productivity, for later investigations of the roles of phytoplankton diversity in biogeochemical cycles in different oceanic regions (using 3-D regional and/or global simulations).Therefore, these two contrasting stations were used to provide a single set of parameter values by fitting the model to observations before the obtained model was validated against data from another independent station (ALOHA).Hence, CITRATE 1.0 consists of the following key features.
1.It models the mean and variance of a continuous phytoplankton size (i.e., log cell volume; µm 3 ) distribution and incorporates "trait diffusion" to sustain size diversity (Merico et al., 2014).
2. It contains an iron cycle in addition to the nitrogen cycle because in the subarctic and equatorial Pacific iron instead of nitrogen should be the main limiting nutrient for phytoplankton growth (Behrenfeld et al., 2006;Fujiki et al., 2014).
3. The phytoplankton cells have variable chlorophyll-tocarbon (θ ) and nitrogen-to-carbon (Q N ) ratios that respond to light and nutrient conditions in a realistic fashion.

4.
A single set of model parameters are optimized against field observational data at two time-series stations in the northwestern Pacific.

Description of the ecosystem model
CITRATE 1.0 contains nine tracers in total: dissolved inorganic nitrogen (DIN, abbreviated as N in all the equations; unit: µmol N L −1 ), phytoplankton biomass (P ; µmol N L −1 ), microzooplankton biomass (MIC; µmol N L −1 ), mesozooplankton biomass (MES; µmol N L −1 ), detritus in terms of nitrogen (D; µmol N L −1 ) and iron (D Fe ; nmol Fe L −1 ), dissolved iron (fer; nmol Fe L −1 ), and the products of P l and P (v+l 2 ) where l (ln µm 3 ) is the phytoplankton mean log cell volume and v ((ln µm 3 ) 2 ) is the log volume variance (Fig. 1).We assume that phytoplankton size is the master trait that determines all physiological functions (Litchman et al., 2007;Finkel et al., 2010;Edwards et al., 2011Edwards et al., , 2012Edwards et al., , 2015;;Marañón, 2015).We also assume that phytoplankton size follows a lognormal distribution, which is supported by some observational data (Finkel, 2007;Quintana et al., 2008Quintana et al., , 2016)).Since l and v are not real standing stocks that can be directly transported in hydrodynamic models but are emergent properties of phytoplankton size structure, we follow Bruggeman (2009) to use the raw moments of biomass probability (i.e., P l and P v + l 2 for mean and variance) as independent tracers involved in transport.All the assumptions made here will be discussed later in Sect. 4. Below we describe the equations for each tracer.For simplicity, phytoplankton cells are assumed not to excrete inor-ganic nitrogen or to have any natural mortality to be converted into detritus.Phytoplankton are eaten by both microzooplankton and mesozooplankton: where µ com is the phytoplankton specific growth rate (d −1 ) of the whole community (i.e., integrated over the whole size spectra).The equation of µ com , along with those of l and v, will be described later in Sect.2.2.E z is the activation energy (in electron volts [eV], 1 eV = 96.49kJ mol −1 ) for heterotrophic processes; g max,i (i = 1 for microzooplankton and 2 for mesozooplankton) is zooplankton maximal grazing rate (d −1 ).K P ,i is the grazing half-saturation constant of zooplankton.Here we have assumed that zooplankton grazing follows a Holling type III functional response.P T,i is total palatable prey concentration for zooplankton (µmol N L −1 ), the details of which will be given later in Sect.2.3.If zooplankton grazing has no size selectivity on phytoplankton, then P T = P .We assume that microzooplankton preferably feed on small phytoplankton, while mesozooplankton prefer large phytoplankton (Table 1).Mesozooplankton also feed on microzooplankton.More descriptions of zooplankton size-dependent grazing will be given later; z is the depth of the model grid (m).K v is the vertical eddy diffusivity (m 2 s −1 ).
The total amount of food ingested by zooplankton is divided among three fates: zooplankton net growth, excretion into the inorganic nitrogen pool, and defecation of unassimilated food into the detritus pool (Buitenhuis et al., 2010).Mesozooplankton mortality is set to be proportional to the squares of its biomass and is also converted into detritus pool.As such, the dynamics of microzooplankton and mesozooplankton as described by the following.Maximal microzooplankton specific ingestion rate for 1.
Detritus is converted to DIN at a rate (R dn , d −1 ) that has the same temperature sensitivity with zooplankton grazing.Detritus is also assumed to have a constant sinking rate (W d , d −1 ).
where unass i represents the fraction of unassimilated food by zooplankton.DIN is taken up by phytoplankton and is replenished by zooplankton excretion, detritus regeneration, and diffusion from the depth: The sources and sinks of fer largely follow DIN with an additional source (atmospheric deposition; Fe depo ) and sink (scavenging; fer scav ; Aumont et al., 2003;Buitenhuis et al., 2010;Nickelsen et al., 2015).
To translate between nitrogen and iron in phytoplankton and zooplankton, a constant fer : N ratio (R fer_N ) of 0.0265 is assumed.The data of monthly atmospheric deposition of total soluble iron are extracted from Scenario III in Luo et al. (2008).Following Nickelsen et al. (2015), the iron scavenging rate (fer scav ) is composed of both the linear scavenging rate (k scm ) and the particle absorption rate (k sc ): in which Fe prime is the concentration of free iron.
where k eq is the equilibrium constant between free iron and ligands and is assumed to depend only on temperature: Note that T is in absolute temperature (K) and l fe is the total iron ligand concentration that is assumed constant (0.6 nM).
The equation for D Fe is as follows.

Continuous trait-based phytoplankton model
Following the moment closure techniques in Merico et al. (2009) and the introduction of "trait diffusion" (Merico et al., 2014), the equations for µ com , l, and v can be written as follows.
Here, µ(l) is the phytoplankton growth rate (d −1 ) at mean size l and u is the trait diffusion parameter, which describes the probability of the parental size l(i) changing to adjacent size values l(i − 1) or l(i + 1) in offspring cells (Merico et al., 2014).Equation (7a)-(7c) are approximations because the higher-order moments, such as the skewness and kurtosis, have been ignored and a Gaussian distribution needs to be assumed for l; dg i (l) dl and d 2 g i (l) dl 2 are the first and second derivatives of zooplankton clearance rate (d −1 ) against phytoplankton size and will be described in detail in Sect.2.3.
The equations of P l and P (v + l 2 ) are the following.
dz Following previous studies (Flynn, 2003;Geider et al., 1997;Follows et al., 2007;Chen and Laws, 2017), phytoplankton growth rate (µ) depends on temperature (T , K), light (I , W m −2 ), DIN, and fer: in which µ m is a function of T : The trait parameters µ m , K N , K fer , and α c are all dependent on cell size l.
µ m = µ 0,m e α µ l+β µ l 2 (10a) K fer = K 0,fer e α fer l (10c) Equation (10a) follows that the maximal phytoplankton growth rate is a unimodal function of phytoplankton size (Chen andLiu, 2010, 2011;Marañón et al., 2013).It is worth noting that the light term of phytoplankton growth (the right side of Eq. 8) is usually modeled as 1 − e − αcI µm (Flynn, 2003), in which both α c and µ m are dependent on size.We use α I to represent the net effect of size on α c µ m for mathematical convenience, which leads to Eq. ( 8).
Following Flynn (2003), we have derived equations to directly calculate phytoplankton chlorophyll-to-carbon (θ, g Chl (mol C) −1 ) and nitrogen-to-carbon (Q N , mol N (mol C) −1 ) ratios from ambient light and nutrient levels: where θ min and θ max are minimal and maximal Chl : C ratios, respectively.Q min and Q max are minimal and maximal N : C ratios, respectively.The total Chl a concentrations (Chl, µg L −1 ) and net primary production (NPP, µg C L −1 d −1 ) integrated over the whole size spectra can be calculated as follows.
To calculate the fractions of Chl within a size range (i.e., < 1, 1-3, 3-10 and > 10 µm), we had to discretize the size spectra into 60 even size classes between l − 6 √ v and l + 6 √ v and calculated the µ, α c , K N , Q N , θ , and eventually the Chl of each size class following Eq.(11a)-(11c).This is because the distributions of Chl do not follow the lognormal distribution of cell volume and an analytic solution is not yet available for calculating only a fraction of Chl.Fortunately, this approach only adds a minor computational cost because we only need to calculate the size-fractionated Chl once per day when saving model outputs.

Zooplankton size-dependent grazing
Following Smith et al. (2016), the ingestion rate of zooplankton on size class l can be formulated as = g max Zρ(l)P (l) where G(l) is the zooplankton ingestion rate (µmol N L −1 d −1 ) on the size class l and ρ(l) is the relative grazing preference on size class l.Z is the biomass of either microzooplankton or mesozooplankton; ε is the food other than phytoplankton (ε = 0 for microzooplankton and MIC for mesozooplankton).P T (total palatable phytoplankton food) is formulated as with P (l) as the phytoplankton concentration at size l.
Zooplankton clearance rate (g, d −1 ) on size class l can be formulated as For mathematic convenience, we parameterize ρ(l) = e bl+c , where b and c are constants.P T can be approximated as Thus, the first derivative of g(l) can then be derived as follows.

One-dimensional (1-D) model
The 1-D model focuses on the upper 150 m of the ocean.The vertical grid, a total of 30 layers, follows a stretched vertical coordinate with increasing resolution towards the sea surface (surface stretching parameter = 2.0), similar to that used in the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams, 2005).For computational efficiency, instead of explicitly solving the complete moment, temperature, and salinity equations, we imported the physics variables that are directly relevant to the ecological processes from external data products.Four types of external physics forcing data were imported into the 1-D model: vertical eddy diffusivity (K v ), surface photosynthetic available radiation (PAR 0 ), atmospheric dust deposition, and vertical temperature profiles.Vertical advection of water was neglected, which has been shown as relatively unimportant (Fernández-Castro et al., 2016).The most important physics forcing data, K v , determined the upward nutrient flux to the upper euphotic zone and were im-ported from the output of a three-dimensional (3-D) eddypermitting model targeted for the North Pacific (Hashioka et al., 2009).This 3-D model was able to faithfully simulate the Kuroshio Current and the spatial distributions of the Chl a fields.The extracted vertical profiles of K v were also consistent with the in situ estimated mixed layer depths (MLD) at the three stations (Fig. 2).PAR 0 values were imported from SeaWIFS satellite monthly climatology products.Seasonal temperature vertical profiles were imported from WOA2013 monthly climatology.
Light levels (I z ) at depth z were calculated based on PAR 0 and Chl a concentrations following the Beer-Lambert law: in which K w and K chl are the attenuation coefficients for seawater and Chl a, respectively.To realistically estimate the average light field that a phytoplankton cell should experience in a mixing water column (Franks, 2015), the ambient light level for phytoplankton within the surface mixed layer is calculated as the average light throughout the surface mixed layer, which is defined as the deepest depth with K v > 10 −3 m 2 s −1 .This calculation is based on Eq. ( 1) in Franks (2015), which gives the average time for a phytoplankton cell to move 100 m (an approximate estimate of MLD) at the local diffusivity of 10 −3 m 2 s −1 as roughly half a day.However, to compare with in situ NPP estimates that were calculated from incubation bottles without continuous mixing, phytoplankton µ, θ , and Q N are recalculated from I z based on the Beer-Lambert law (Eq.13).
The initial condition of inorganic nitrogen was set to the vertical profile of nitrate in January of the World Ocean Atlas (WOA) 2013 monthly climatology.Initial phytoplankton, microzooplankton, and detritus biomass were all set to 0.1 µmol N L −1 in each grid.Mesozooplankton biomass was initialized as 0.05 µmol N L −1 .Initial D Fe concentra-tions were set as detritus times R fer_N .Initial phytoplankton mean log size (l) and log size variance (v) were set as −2.2 log µm 3 and 0.09 (log µm 3 ) 2 , respectively.Initial dissolved iron concentration was set to the vertical profile of iron in January from a 3-D global biogeochemical model output (Aumont et al., 2003).The time step of the model was 5 min.All the fixed model parameters are shown in Table 1 and the model parameters that are optimized to match observational data are shown in Table 2.
We employed a Dirichlet boundary condition at the bottom for DIN and fer with the values predefined by the WOA2013 climatology and the model output from Aumont et al. (2003), respectively.For other tracers, we assumed no diffusive flux at the bottom.Detritus was allowed to sink out of the system with the loss of nitrogen and iron replenished by diffusion.

Delayed-Rejection Adaptive Metropolis-Hasting Monte Carlo (DRAM) algorithm
The Metropolis-Hasting Monte Carlo (MHMC) algorithm aims to find the posterior distribution (including mean and covariance matrix) of the parameter vectors given the data provided.The key here is to develop an appropriate proposal covariance matrix (P cvm ), which determines the magnitude and direction of the proposed perturbations to the parameter values, as the algorithm explores the parameter space.At each iteration of the algorithm, the newly proposed parameter set is either accepted or rejected based on the model-data mismatch, as explained below.In the classical random walk MHMC algorithm, the P cvm must be specified by the user to achieve sufficient acceptance rates for the proposed parameters, which typically requires a great deal of effort and many trials.
The adaptive MHMC (Haario et al., 2001) uses the already accepted parameters to approximate P cvm , which is periodically updated as more simulations are conducted.Specifically, the P cvm is tuned based on the covariance matrix (C vm ) of the already accepted parameter sets after a fixed number of iterations following Gelman et al. (2014; i.e., P cvm = C vm × 2.4 2 /d, where d is the length of the target parameter vector).Thus, the algorithm alters the magnitude and direction of proposed "jumps" in order to efficiently explore the parameter space.
With the delayed-rejection MCMC (Mira, 2001), when a newly proposed set of parameters is rejected, P cvm is temporarily downscaled (to 1 % of the original P cvm in our case) and a second set of parameters is proposed based on the rejected parameters and the downscaled P cvm .This approach is particularly efficient because low acceptance rates typically result when the P cvm is too large (the parameter jumps are too wide) to find the target distribution of the parameters.Temporarily reducing P cvm can substantially increase the acceptance rate.By using multiple stages of P cvm , the algorithm can also effectively deal with the problem of non-Gaussian posteriors, which can reduce the efficiency of the adaptive MHMC (Haario et al., 2006).
The DRAM algorithm, built upon the classic Metropolis-Hasting Monte Carlo (MHMC) algorithm, incorporates the merits of both the adaptive and delayed-rejection MHMC algorithm to increase the acceptance rate and thus more efficiently find the target distribution of parameter values (Haario et al., 2006;Laine, 2008).It has been shown to better explore the parameter space compared to other algorithms such as the families of simulated annealing, possibly because of its two-stage proposal covariance matrices (Villagran et al., 2008).Compared with the widely used ensemble Kalman filter, DRAM is more suitable for the nonlinear response typical of ecosystems (Annan and Hargreaves, 2007).
Here we briefly outline the DRAM algorithm.For further details and proofs, see Haario et al. (2006) and Laine (2008).
1. Initialize the parameter values and P cvm , assuming no correlation among parameters, and a standard deviation equaling one-sixth of the difference between the maximal and minimal value for each parameter, respectively (Table 1).
2. Run the model with the current parameter values (θ curr ) and calculate the likelihood (L).Note that all the parameter values must be within the boundaries shown in Table 2.
3. Propose a new set of parameters (θ pro ) based on θ curr and P cvm , rerun the model, and obtain a new likelihood (L 1 ).
4. If the ratio of L 1 /L is larger than a random number between 0 and 1, then accept θ pro (θ curr = θ pro ) and return to step 2.
5. Otherwise, propose a second set of parameters (θ pro2 ) based on θ pro and P cvm2 (= 0.01 P cvm ), rerun the model, and obtain the second likelihood (L 2 ).

If the ratio of
is larger than a random number between 0 and 1, then accept θ pro2 (θ curr = θ pro2 ) and return to step 2. Otherwise retain the current position, θ curr .Here q 1 (y, x) is the probability of proposing y given x, and q 2 (z, y, x) is the probability of proposing z given x and y.
7. After a certain interval, update P cvm based on C vm calculated from the accepted θ .
To increase the computational efficiency and avoid being trapped in local minima due to insufficient chain length, we modified the DRAM algorithm for parallel computing (Calderhead, 2014).That is, we initialize θ and P cvm simultaneously for n processes.Each process runs the above procedure from (1) to ( 7) except that at ( 7 are consolidated to update the global estimate of P cvm , which is then distributed to all subprocesses to propose new θ .Preliminary model runs suggested that from the third year, the model reached a quasi-steady state, exhibiting regular seasonal cycles under the climatological forcing (Fig. 3).As such, we ran the model for 4 years and the output of the final year was used for validation against observational data.The model outputs were linearly interpolated to the observational depths and time.To allow fair comparisons among different data types and downplay the effects of extreme values, both the model outputs and observational data were transformed to their 1/4 power and normalized between 0 and 1 to achieve a quasi-normal distribution before calculating the sum of squared errors (SSqEs).
SSqEs k,i is the sum of squared errors of data type i at station k, n k,i is the number of observations for data type i at station k, o k,i,j is the observed j th value for data type i at station k, and o k,i,min and o k,i,max are minimal and maximal observed values for data type i at station k, respectively (note that for all size fractions of Chl a, we intentionally set o k,i,min = 0 and o k,i,max = 1 to minimize the effects of the large measurement variability); m k,i,j is the value linearly interpolated from model outputs to the same depth and date of o k,i,j .
Following Laine (2008), the likelihood function was calculated as the product of the exponential of the sum of squared errors scaled by a measure of the model-data error for each data type, respectively: in which σ k,i is the standard deviation of the Gaussian errors of data type i at station k.
Following Laine (2008), we applied Gibbs sampling, which estimates the distribution of each σ k to match the ensemble distribution of model output to that of the data.This entails assuming that the prior of 1/σ k,i follows a gamma distribution, with the prior mean as S 2 0 and prior accuracy as n 0 .At each step the value of 1/σ k,i was sampled from a conditional gamma distribution . The model parameters were assumed to follow multivariate normal distributions.The likelihood function contributed by the priors of the parameters was in which n p is the number of parameters to be estimated, and γ i and η i are the prior estimates of the ith parameter and its standard deviation, respectively (Table 2).Values of η i were calculated as one-sixth of the difference between the preset maximal and minimal parameter boundaries; θ i is the current parameter value.The MCMC chain was run for an ensemble of 10 000 simulations with five processes running in parallel (i.e., a total of 50 000 parameter sets were obtained).Although the model contains more than 20 parameters, we only selected 9 parameters for optimization to minimize the possibility of parameter unidentifiability and avoid optimizing highly correlated parameters such as g max and K P simultaneously (Table 2).

Observational data
For stations K2 and S1, the observations including MLD and nine types of data (DIN, Chl, NPP, PON, fer, and four size fractions of Chl) were obtained from the K2S1 project (https: //ebcrpa.jamstec.go.jp/k2s1/en/index.html;Honda, 2016; Table 3).The observations spanned from 2010 to 2013 at seasonal sampling frequencies.Some of the data have been published in Wakita et al. (2016), Fujiki et al. (2016), Matsumoto et al. (2016), and Sasai et al. (2016).DIN was calculated as the sum of nitrate, nitrite, and ammonia, which were measured with a continuous-flow analyzer (QuAAtro 2-HR system; BL-Tech).Chl was measured using the nonacidification method following Welschmeyer (1994).NPP was measured with the technique of NaH 13 CO 3 uptake (Matsumoto et al., 2016).PON was measured by using an elemental analyzer (Wakita et al., 2016).Size fractions of Chl were measured by filtering seawater sequentially through 10, 3, and 1 µm polycarbonate membrane filters and finally a GF/F glass-fiber filter.The filters were soaked in N,N-dimethylformamide (DMF) and chlorophyll concentrations retained on the filters were measured with the same protocol as total Chl (Fujiki et al., 2016).
For station ALOHA, the data were downloaded from http: //hahana.soest.hawaii.edu/hot/.All the data were pooled to-gether to generate a quasi-climatological seasonal pattern and inter-annual variations were treated as random noise.To improve data coverage, we also included the nitrate data from the World Ocean Atlas (WOA) 2013 for observed DIN.Due to the lack of in situ observational data, the data of fer were obtained from a global biogeochemistry model (Aumont et al., 2003).To calculate MLD from depth profiles of temperature and salinity, MLD was defined as the first depth at which the seawater density exceeds surface density by 0.125 kg m −3 (Shigemitsu et al., 2012).

External physics forcing
The validity of external physics forcing data, particularly the vertical mixing that determines upward nutrient diffusive supply to the surface mixed layer, is essential for correct results and parameter optimization with the ecosystem model.Here we show in Fig. 2 a representative year of seasonal variations in K v , temperature, surface PAR, and atmospheric iron deposition.Vigorous winter mixing precedes summer water column stratification at K2 and S1, while the seasonal variations in mixing are less pronounced at ALOHA.At all three stations, the model estimates of mixed layer depths are consistent with those measured from in situ temperature and salinity profiles (Fig. 2b, f, j).Water temperatures and surface PAR values at the subarctic station K2 are significantly lower than at the subtropical stations S1 and ALOHA.The station K2 is also characterized by a pronounced spring peak in atmospheric dust deposition.

Parameter optimization and sensitivity analysis
For all five parallel subprocesses, the log likelihood continued to increase with the number of model runs and reached a plateau after 1000 iterations (Fig. 4).For most (but not all) types of data, model-data mismatches (SSqEs) consistently decreased.Comparing the two stations, the model fits to the Chl and NPP were better at station K2 than S1.The model fits to the size fractions of 1-3 µm were better at S1 than K2.
Most values of the optimized parameters fell into reasonable ranges (Table 2; Fig. 5).For example, the estimated K 0,N is close to the value (0.2 µM) given in Ward et al. (2012).For some of the parameters, such as W d and u, the final optimized value differed substantially from initial estimates, which is an expected outcome of the algorithm striving to match the nine different types of observations at both stations with contrasting environments.Below we show some preliminary results of a sensitivity analysis, particularly those differing with a priori estimates (Table 4).
The mean µ 0,m estimated from laboratory phytoplankton data is around 0.4 d −1 , half of the optimized value (Chen and Laws, 2017).Reducing µ 0,m to 0.4 d −1 mainly generated worse fits to the size fractions of < 1 µm of Chl at both sta-  tions.This is because the lower phytoplankton growth led to higher nutrient concentrations and lower estimates of < 1 µm fractions.
The estimate of W d of 20 m d −1 is a relatively high sinking speed.Reducing W d to 10 m d −1 only led to slightly worse fits to DIN data at station S1 (but better fits to DIN at K2) and overall did not deteriorate the results substantially.
The estimate of m z (0.2 (µM N) −1 d −1 ) is also at the high end of those used in the literature.We found that the model results were quite sensitive to the value of the closure term m z .Reducing m z to 0.1 (µM N) −1 d −1 led to higher meso-plankton biomass and generated much worse fits, particularly for DIN at K2.
We also tested whether we could assume that the light component of phytoplankton growth is size independent (i.e., α I = 0).The results suggested that with α I = 0, the model predicted much worse fits to the data.An optimized value of −0.26 for α I is also consistent with the size scaling relationship of light-dependent growth in Finkel (2001) and Edwards et al. (2015), suggesting that light limitation could drive phytoplankton to be small.
The optimized trait diffusion coefficient (u) was much higher than in Acevedo-Trejos et al. (2016).Reducing u to    0.05 led to worse fits to the size-fractionated chlorophyll since lower size variance failed to capture the observed size scatter.It also relates to the limitation of the model that has to assume a lognormal distribution of size (see Sect. 4.2.1).However, an abnormally high u could drive the model to unstable conditions in which the size variance kept increasing.

Comparison between best model outputs and observation
The best model outputs in terms of the highest likelihood could capture most of the observational patterns quantitatively (Figs.6-9).At both stations, the model could reproduce the vertically increasing trend of DIN with depth and the higher surface concentrations of DIN during winter than summer and autumn.It is noteworthy that the model could also successfully reproduce the relatively abundant summer DIN concentrations at the surface at station K2 due to the incorporation of iron and light limitation.The vertical and seasonal patterns of Chl a and NPP could also be well reproduced at station K2.The only problem is that at station S1 the high NPP at the surface could not be well reproduced (Fig. 7).
Validation against observed phytoplankton size data is critical for testing CITRATE 1.0 in which phytoplankton size structure is the core component.The model could reproduce most patterns of the proportions of size-fractionated Chl at both stations (Figs. 8, 9).For example, the model correctly reproduced the relative dominance of picophytoplankton (< 3 µm) at both stations, although nitrate concentration was high at station K2.The seasonal and vertical fractions of 3-10 µm were generally well simulated at both stations, except for an artificial surface peak at K2.The model could also simulate the relatively larger sizes at K2 than at S1.
We also note some deficiencies of the model.At both stations, the fractions of > 10 µm Chl were close to zero at both stations in the model in contrast to the substantial fractions of > 10 µm during summer at K2 and in the winter at S1.The model also tended to overestimate the 1-3 µm fractions at both stations and underestimate the < 1 µm fractions occasionally.All these problems relate to the assumption of a fixed trait distribution as discussed later.

Modeled seasonal patterns of nutrients, phytoplankton biomass, mean size, and size diversity
At both stations, DIN concentrations were higher during winter in the surface mixed layer due to more vigorous mixing (Figs. 10,11).Significant drawdown of DIN occurred in surface water following water column stratification, which occurred earlier at S1 than K2.At station K2, after an increase during June and July due to the peak in atmospheric deposition, dissolved iron concentration also decreased in the fall due to phytoplankton uptake.By contrast, surface iron concentrations accumulated from late summer to fall due to nitrogen limitation at station S1.
In accordance with the DIN patterns, higher concentrations of Chl a were found during winter at station S1, which resulted from both increased phytoplankton biomass and chlorophyll-to-carbon ratios (Fig. 11).Starting from spring to fall, subsurface maximal layers of Chl a formed and progressively deepened with time.By contrast, at station K2, Chl a concentrations peaked in summer and subsurface chlorophyll maximum layers were not evident (Fig. 10), suggesting light limitation played a stronger role in limiting phytoplankton growth at K2 than S1.
At both stations, in spite of the nutrient increases in winter, phytoplankton mean size peaked in spring or summer.This is mostly contributed by the light limitation on large cells, which can be reflected by the negative value of α I (Table 2).At both stations, the main periods of size increases were in spring when the light level increased and there were still nutrients left from winter mixing.The increases in light were contributed by both increases in surface PAR and shallower mixing.Nutrient (dissolved iron in the case of K2) depletion together with light decreases led to negative values of dµ(l)   dl in late spring or summer at both stations, resulting in subsequent decreases in mean size.In general, the modeled mean  sizes were significantly larger at station K2 than S1, mainly due to less severe nutrient limitation.
The modeled patterns of size variances (i.e., size diversity) are the focus of CITRATE.Within the surface mixed layer, modeled phytoplankton size diversity showed an opposite pattern with mean size, with peaks in fall at S1 and in winter at K2 (Figs. 10,11).At first glance, we also seemed to find a negative correlation between the growth rate µ com and size diversity at both stations (Fig. 12a).When growth rates were high, size variances were low, and vice versa.The paired scatterplots between µ com and size variances in surface waters suggested that these two quantities were not linearly correlated, particularly at S1. Instead, their relationships depended on the timing of the season.At station S1, during the transition from the end of winter to early spring, the phytoplankton cells experience a rapid increase in growth rate without much change in size diversity.During the rest of spring, the phytoplankton growth rate decreased from the maximum to nearly the minimum, while size diversity first underwent a phase of moderate decrease and then recovered.From the beginning of summer to mid-fall, there were no big changes in growth rate, but size diversity increased dramatically.From mid-fall to the beginning of the winter, the phytoplankton growth rate increased, but size diversity decreased to winter values.At station K2, the variability in size diversity was smaller, with high growth rates and low size diversity in summer and the opposite patterns in winter.
We decomposed the different factors affecting the dynamics of size diversity in surface waters at both stations (Eq.7c, e; Fig. 12b, c).Three points need to be mentioned.First, the calculated net combined effects, including the second derivatives of growth and grazing d 2 µ(l) dl 2 and d 2 g(l) dl 2 , trait diffusion d 4 µ(l) dl 4 and µ(l) , and vertical mixing (i.e., diffusion), were consistent with the net changes in size variances (some minor differences were because we saved the above quantities at a daily interval that could not account for the changes within 1 day), validating our computation.Sec- Figure 9.The same as Fig. 8, but for station S1.
ond, the contributions from the second derivatives of growth and trait diffusion (dominated by 2 uµ(l) with the contributions from d 4 µ(l) dl 4 being minor; Eq. 7c) were the two largest terms, which usually offset each other.It is the margin of these two terms plus vertical mixing that drove the changes in size variance.The values of d 2 µ(l) dl 2 were always negative at all times at both stations, suggesting that without "trait diffusion", size variance would decrease toward zero (Eq.7c).This highlights the importance of trait diffusion (which can be interpreted as genetic mutation or transgenerational phenotypic plasticity) to sustain diversity.The values of d 2 µ(l) dl 2 were more negative when growth rates were higher.For example, in early April at S1, the decrease in size variance was induced by a more negative d 2 µ(l) dl 2 (see also Fig. 11h).Similar situations also occurred at the end of December.Third, water column mixing also played a significant role in affecting size diversity, which was the main factor leading to the peak in size diversity in fall in surface waters at S1.The effect of mixing became important because at this time, a sub-surface maximum of phytoplankton biomass still existed below the surface mixed layer.With the deepening of the surface mixed layer, a substantial biomass of phytoplankton was entrained into surface waters and these phytoplankton communities had different trait properties than the surface ones, thereby enhancing diversity (see Sect. 4.1.1 for discussion).
The model also generated reasonable patterns of Chl : C and N : C ratios, which were largely determined by light and nutrient concentrations (Figs. 10i,j;11i,j).Both Chl : C and N : C ratios were high in winter when nutrient concentrations were high and light levels were low due to strong mixing.Both ratios were low in surface stratified waters where nutrient supply from below became diminished due to strong stratification and also light levels became strong due to both increased surface PAR and shallow mixing layers.

Validations of the model at station ALOHA
We used the optimal parameter sets obtained at stations S1 and K2 to run the model at station ALOHA.As there were www.geosci-model-dev.net/11/467/2018/Geosci.Model Dev., 11, 467-495, 2018 In models resolving a number of discrete species, the typical index for the intensity of resource competition under steady state is R * , the lowest nutrient concentration allowing positive net growth (Tilman, 1982;Litchman et al., 2007;Barton et al., 2010).Under nonequilibrium conditions, it is the maximal growth rate instead of R * that determines the outcome of competition (Huston, 1979;Barton et al., 2010).In any case, it is the realized growth rate that determines the outcome of competition.Compared to R * , the second derivative d 2 µ(l) dl 2 has two advantages as a proxy for quantifying the intensity of competition: (1) it applies under both equilibrium and nonequilibrium conditions and (2) it circumvents the problem of tracking many species.Using this approach, it is straightforward to test some ecological theories such as Huston's "general hypothesis of species diversity" (Huston, 1979).For example, the absolute magnitude of d 2 µ(l) dl 2 correlates positively with µ (Fig. 13), indicating that resource competition is more intense when growth rates are high.This is a mathematical manifestation of the verbal argument of the "dynamic equilibrium theory" proposed in Huston (1979), who emphasized that in natural environments where equilibrium is rarely achieved, fast-growing species tend to outcompete slow-growing species (see also Barton et al., 2010), and hence growth rates play a greater role in determining diversity than R * values.
Similarly, Eq. ( 7b) concisely specifies the factors affecting mean phytoplankton size.In fact, Eq. ( 7a)-(7c) can be understood as derived from a Taylor expansion representing an infinite number of discrete trait classes (Merico et al., 2009).Hence, even if a discrete version of a diversity model is used, it may also be helpful to calculate the terms in Eq. ( 7a)-(7c) in order to understand the factors affecting species diversity, biomass, and productivity.
The set of equations in Eq. ( 7) also provides an excellent platform to investigate the underlying mechanisms for the relationship between biodiversity and ecosystem functioning dl 2 .MIC and MES grazing equates to −v 2 d 2 g i dl 2 , "d4µ / dL4" equates to v 2 u d 4 µ dl 4 , and "trait diffusion" equates to 2 uµ.All the derivatives are evaluated at the mean size."Diffusion" means the contribution to the changes in size variance induced by diffusion with the underlying grid."Net effect" means the sum of the above terms."Net changes" mean the difference in size variance between adjacent days.(c) The same as (b), but at station K2.
(productivity, in this case), which have been extensively studied (Loreau et al., 2001;Tilman et al., 2014).While the negative relationship between productivity (µ com ) and diversity suggests that enhanced productivity can induce greater competition and reduce diversity (Huston, 1979), diversity can also be affected by other factors besides competition.
The incorporation of trait diffusion originally developed for continuous trait-based models (Merico et al., 2014) provides a means of representing mutation and other processes that sustain diversity, thus linking ecological and evolutionary processes (Rosenzweig, 1995).This allows for control of the level of diversity in simulation experiments such as those conducted herein to investigate diversity-productivity rela-tionships.The increasing effect of trait diffusion with growth rate is consistent with the metabolic theory of ecology in that metabolic rates, which are closely coupled with growth rates and generation time, are expected to correlate with mutation rates.Therefore, growth rates are expected to affect speciation and potentially contribute to the latitudinal diversity gradient (Rohde, 1992;Allen et al., 2006;Dowle et al., 2013).Our results have shown that trait diffusion can be the largest term counterbalancing competitive exclusion (Fig. 13).Without considering this mechanism, diversity could be underestimated in productive waters due to strong competition.
The approach of transporting trait moments across spatial grids, originally developed by Bruggeman (2009), also allows water mixing to affect diversity patterns.Although this approach is not perfect (see Sect. 4.2.2 and Fig. 14), it does allow the mixing of two communities with different mean traits to generate trait variance greater than the weighed mean variance of the two original communities.The larger difference in the mean traits, the greater the increase in trait variance upon mixing.Consider the case of mixing two communities with biomass P 1 and P 2 , mean size l 1 and l 2 , and size variance v 1 and v 2 .The biomass and mean size of the mixed community are P 1 + P 2 and P 1 l 1 +P 2 l 2 P 1 +P 2 , respectively.After some algebraic manipulation, we can derive the size variance (v ) after mixing.
Thus, it is clear from Eq. ( 15) that the difference between v and the biomass weighed mean variance P 1 v 1 +P 2 v 2 P 1 +P 2 depends on the difference in mean traits.Hence, mixing can enhance diversity to the extent that the traits of the original communities differ.Barton et al. (2010) have shown that the "hotspots" of high phytoplankton diversity are usually located along areas where mixing is strong enough to allow the coexistence of multiple populations with different traits.Our simulations are consistent with that view and show that vertical mixing can significantly enhance diversity, particularly during ocean mixed layer entrainment.

Flexible stoichiometry
We also consider realistic phytoplankton physiology and optimized model parameters guided by real data.For example, our model has incorporated some features of phytoplankton plasticity (acclimation) such as variable Chl : C and N : C ratios.Although, for the sake of simplicity, these variable ratios do not directly influence the phytoplankton specific growth rate as in Geider et al. (1997), they are able to reproduce the high Chl : C ratios in the DCM layer, thus providing a more realistic mechanism for the formation of the DCM layer than with models that assume fixed ratios (Fennel and Boss,Figure 13.The same as Fig. 6, but for station ALOHA. 2003).Similarly, the variable N : C ratio also allows phytoplankton cells to achieve higher carbon-based NPP in surface waters compared to models with fixed N : C ratios (Christian, 2005).Although cellular chlorophyll and nitrogen quotas are not calculated as independent tracers, model comparisons suggest that more complex models do not always yield better fits to the data (Flynn, 2003).

Realistic mechanisms for controlling phytoplankton size structure
In CITRATE 1.0 we have provided both bottom-up and topdown mechanisms to affect the size structure of phytoplankton.First, we employ an observation-based unimodal relationship between maximal growth rate and size to give the nanophytoplankton the advantage under nutrient-replete conditions (Chen andLiu, 2010, 2011;Marañón et al., 2013), thus allowing a trade-off between nutrient affinity and maximal growth rate within the pico-and nano-size range.Thus, bottom-up factors alone are sufficient to reproduce the observed decrease in the fraction of small phytoplankton with nutrient enrichment (Marañón et al., 2012).We also impose a size-dependent feeding preference of zooplankton based on the general understanding that smaller microzooplankton tend to prefer smaller phytoplankton, whereas larger mesozooplankton tend to prefer larger phytoplankton (Frost, 1972;Hansen et al., 1994;Liu et al., 2005;Ward et al., 2012).These top-down factors have additional effects on phytoplankton size structure.Our assumption about the preference of microzooplankton for small phytoplankton is similar to Terseleer et al. ( 2014) and Acevedo-Trejos et al. (2015), who assumed a combination of decreasing maximal phytoplankton growth rate with increasing size and a grazing preference for small phytoplankton in order to offset the growth advantage of small phytoplankton in eutrophic waters.In our case, small phytoplankton lose the advantage in eutrophic waters, where larger phytoplankton grow faster because of the imposed unimodal relationship between maximal growth rate and size.Meanwhile, in eutrophic waters, mesozooplankton dominate and preferentially feed on larger phytoplankton to balance the growth advantages of larger cells.Interestingly, counter to our intuition, field incubation experiments have often found that microzooplankton feed on diatoms faster than on picophytoplankton and that diatoms grow faster than picophytoplankton even in oligotrophic waters (Latasa et al., 1997;Zhou et al., 2015).These results raise a paradoxical question: how can diatoms grow so fast with negligible nutrients in oligotrophic waters, but without accumulating high biomass?Whether this is because of experimental bias is an open question.The feeding preference of mesozooplankton on large prey seems less disputable (Frost, 1972;Liu et al., 2005), but see Terseleer et al. (2014) for an assumption of decreasing feeding preference of copepods for large diatoms.This implies strong top-down control of large phytoplankton in eutrophic waters where mesozooplankton dominate, limiting the biomass of large phytoplankton.However, this implication is at odds with the common observation that large phytoplankton dominate total biomass in eutrophic waters (Marañón et al., 2012).Future refinements might include a unimodal feeding preference similar to the grazing kernel proposed earlier (Hansen et al., 1994;Poulin and Franks, 2010).In any case, for model calibration and validation, more and better data are needed concerning the size scaling of both phytoplankton traits and zooplankton grazing preference.

Assumption of trait distribution
To facilitate the calculation of trait moments, a certain distribution has to be assumed for the trait (Merico et al., 2009(Merico et al., , 2014)).The lognormal distribution can be fitted well to empirical data (Quintana et al., 2008(Quintana et al., , 2016)), and because of its mathematical convenience it has been widely used in continuous size distribution models (Terseleer et al., 2014;Acevedo-Trejos et al., 2015, 2016;Smith et al., 2016).For these reasons we have assumed a lognormal distribution in the present study.
However, other probability distributions can also describe phytoplankton size.In the literature, phytoplankton abundance (N , cells L −1 ) within the size interval from V to V + dV is more often modeled as a power-law function of cell volume V (unit: µm 3 ; Gin et al., 1999;Cavender-Bares et al., 2001;Cermeño et al., 2006): where N 0 represents the abundance of phytoplankton having cell volume 1 µm 3 , and α is the exponent of the power law.
Because models typically represent phytoplankton biomass instead of abundance, we can convert the abundance to biomass (B, µm 3 L −1 ): Although the power law of Eq. (16b) may seem to be a suitable alternative distribution for continuous size-based models, empirical data suggest that α tends to vary between −0.7 and −1 (Cermeño et al., 2006), which means that the exponent (α + 1) of the power law relating B and V should be between 0 and 0.3.In this case both the mean and variance of the power-law distribution as shown in Eq. ( 16b) are infinite (Newman, 2005).This problem can be solved by adding an upper cutoff via an exponential truncation (Clauset et al., 2009): where λ is a positive constant.
Whether the power law or the lognormal distribution fits better to empirical data has been widely debated in the literature, and many results show that both can fit the data equally well (Allen et al., 2001;Mitzenmacher, 2004;Clauset et al., 2009).This is not surprising given that the two distributions are intrinsically connected (Mitzenmacher, 2004;Newman, 2005).We suspect that the power law with an upper cutoff may be able to better capture the right skewness of phytoplankton size distributions, as is common in oligotrophic waters where large diatoms coexist with the dominant cyanobacteria such as Prochlorococcus (Campbell et al., 1994;Liu et al., 1997;Villareal et al., 1999).It remains to be investigated whether changing the distribution to the truncated power law can help solve the problem of underestimating the fraction of > 10 µm size in the current study.
Neither the lognormal nor the power law with an upper cutoff can capture multimodal size distributions, as exemplified in Fig. 1b of Marañón (2015) and reported by other studies (Banas, 2011;Bonachela et al., 2016;Coutinho et al., 2016).This is an inevitable consequence of aggregating the description of the entire community into only the three descriptors (i.e., total biomass, mean, and variance), which reduces the degrees of freedom, thus sacrificing detailed accuracy for generality and perspective.
One remedy for this problem might be to assign more functional groups in phytoplankton and assume a probability distribution for each group (Terseleer et al., 2014).Having a number of functional groups also circumvents the problem of size-independent functional differences among phytoplankton, such as the different maximal growth rates of diatoms and dinoflagellates despite their similar sizes (Chen and Laws, 2017).We expect that in the near future such a combination of continuous trait distributions and functional groups will likely provide more realistic representations of marine phytoplankton diversity.

Transport of moments
Another potential problem is the transport of trait moments in ocean circulation models.Unlike nutrients or plankton biomass, trait moments are not real "concentrations" that can be directly involved in advection and diffusion.In general two Gaussian curves differing in area (i.e., total biomass), mean, and variance do not sum to a perfect Gaussian curve (Fig. 14a).Bruggeman (2009) has derived that, following the assumption of normal distribution of traits, the raw moments of the biomass distribution can behave as normal tracers in GCMs.We have shown a few examples of the mixing of communities of different biomass, mean size, and size variance in Fig. 14.These examples demonstrate that when the mean sizes and size variances differ greatly and biomasses are similar, the mixed community may deviate from the assumed normal distribution, making this a poor approximation.For now, we assume that across adjacent grids, phytoplankton communities should in most cases be similar enough for this approximation to work reasonably well.

Lack of multiple traits
As a first step, we incorporated only size as the master trait that affects all physiological functions of phytoplankton.In reality, many phytoplankton functional traits, such as optimal temperature, diazotrophy, and mixotrophy, are independent of size.For example, the optimal growth temperature of phytoplankton is closely related to environmental temperature, but only weakly relates to size (Thomas et al., 2012;Chen, 2015).The optimal growth temperature and irradiances are certainly function traits that deserve to be incorporated into trait-based models (Follows et al., 2007;Norberg, 2004;Edwards et al., 2015) and are expected to strongly affect phytoplankton functional diversity on large scales.

Difficulty in modeling surface peaks in NPP at oligotrophic stations
The near-surface peak in NPP at the oligotrophic stations S1 and ALOHA during summer is not expected if we assume that the source of nutrients comes from below the euphotic zone.Even if variable N : C ratios are used in the model to allow more carbon to be fixed given the same amount of nitrogen near surface waters, surface NPP is still likely to be underestimated even with the presence of N 2 fixation because of phosphorus limitation (Christian, 2005).It is possible that other mechanisms, such as the vertical migration of phytoplankton, need to be taken into account (Villareal et al., 1999;Chavez et al., 2011).Therefore, this problem is not only restricted to CITRATE 1.0.

Optimized parameters for 3-D GCM
One purpose of optimizing a common parameter set for two stations with contrasting environmental conditions is to use this parameter set for 3-D GCMs.This is based on the expectation that a parameter set that can work for the two stations should work for other locations as well.However, our validation exercise at station ALOHA reveals that the parameter set optimized for stations K2 and S1 only succeeds in matching the DIN data well, but underestimates Chl, NPP, and PON at station ALOHA.This suggests that we might be overlooking some unique but important processes at ALOHA.Alternatively, it is also possible that the uneven sampling at K2 and S1 might bias the parameter optimization to some extent.Similar difficulties in parameter optimization have been shown previously (Ward et al., 2010).For optimizing parameters for 3-D GCMs, a better approach might be to use the "transport matrix" technique that has been successfully implemented for some biogeochemistry models (Khatiwala, 2007;Kriest et al., 2017).Nonetheless, our optimized parameters can provide a useful initial estimate for modeling other stations and for use in 3-D GCMs.

Future directions
Considering the above limitations, one future direction is to increase the number of traits in the model to generate more realistic phytoplankton diversity patterns, which requires both an "envelope" function relating the maximal growth rate with the optimal trait value and a relationship between growth rate and trait value for each species (Norberg, 2004).Another refinement as noted above is to model a continuous trait distribution for each functional group, thus combining the continuous trait-distribution and functional group approaches to better capture deviations in overall trait distri- It is relatively easy to couple the one-dimensional CIT-RATE model with 3-D global or regional ocean models in order to model the large-scale patterns of phytoplankton size and size diversity.Furthermore, it should be possible in the near future to optimize parameters for such a 3-D model using the transport matrix technique.In particular, by including both trait diffusion and competitive exclusion it may be possible to begin to untangle the relative roles of ecological versus evolutionary processes in shaping global phytoplankton diversity patterns.

Conclusions
We present a 1-D model with continuous size distribution for phytoplankton (CITRATE).The dynamics of phytoplankton mean size and size variance are directly linked to environmental factors and moments of the size distribution (Eq.7), facilitating an understanding of the underlying mechanisms controlling phytoplankton size and diversity.CITRATE 1.0 also incorporates "trait diffusion" as an eco-evolutionary process to sustain phytoplankton diversity.
We optimized the parameters of CITRATE using the DRAM algorithm, which revealed that the model can faithfully reproduce observed seasonal patterns of inorganic nitrogen, Chl a, and phytoplankton size structure at two contrasting time-series stations.The model structure and associated parameters obtained herein can be useful for 3-D regional and global ocean modeling.
The limitations of CITRATE include its assumption of a lognormal distribution for phytoplankton size as the sole master trait, which to some extent limits the precision with which it can reproduce large size classes of phytoplankton.These limitations and others may be overcome in future studies by building on CITRATE 1.0 to construct more elaborate continuous trait-distribution models capable of reproducing more realistic patterns of phytoplankton diversity.

Code and data availability
The code and data of CITRATE 1.0 are freely available at https://github.com/BingzhangChen/citrateunder the MIT license.

Tutorial.
The code for CITRATE 1.0 (https://doi.org/10.5281/zenodo.1034805) is written in Fortran 90 with the Intel Fortran compiler used.We have tested the codes on macOS Sierra 10.12.5 (i386 processor) and also a GNU/Linux cluster with x86-64 architecture.The user should be familiar with the Fortran language and have some basic knowledge of BASH.Some post-processing scripts are also written in the free software R (version 3.3.2).Before compiling the codes and running the model, the user needs to install the mpi (e.g., openmpi) library for parallel computation.Below we give some instructions and explanations of the codes and how to run the model.
1. Go to the directory in which you want to run the model (we assume that the root directory is under home directory: ∼/).
4. Type "vi run" to change the setting for the model run.
-Test = 0 means a fast run, usually for a formal model run for a large number of iterations.
-Test = 1 means running a model for debugging, which is much slower than the fast run.
The user can also modify the compiler flags depending on the purpose in the script.The user needs to specify the directory where the library of mpifort exists.
5. Type "./run" and the model will compile and an executable (CITRATE) will be generated.
6. Type "vi Model.nml", which contains two namelists.The namelist &Model contains the options for station names, the type of ecological model, the type of nutrient uptake function (one only for CITRATE), and the type for grazing function (four different grazing functions including the three Holling type functions and the Ivlev function).The station name determines the right physics files to be read and the filenames for model output.For now we only allow three possible stations: S1, K2, and HOT.Other station names will generate an error.If the user wants to add more station names, the subroutine Setup_OBSdata within MOD_1-D.f90 is the place to be modified.A number of ecological models besides CITRATE have been developed.It is beyond the scope of the present study to describe all of them in detail.Just note that the model lists are in the Fortran file bio_MOD.f90and some other details are in choose_model.f90and MOD_1-D.f90.
The namelist &MCMCrun contains the options for defining the total length of the MCMC chain, which is at least two, the number of the ensemble runs, the number of days for each model run, whether the model should start from previous runs (Readfile = 1) or start a new run (Readfile = 0), and the number of runs in the historical files (enssig and enspar).
7. After defining all the model settings, type "mpirun -np 5 citrate" and then the model will run with five parallel processes and some outputs will be shown on the screen.Type "mpirun -np 5 citrate > out" to make the model outputs stored in the "out" file.For each model run, the model saves the current parameters into the "enspar" file and the current values of σ and SSqEs into the "enssig" file.In this way, even if the model crashes, the user can pick up the current parameter position and the updated parameter covariance matrix.The model also generates the files of best parameters, best σ and SSqEs files, best model output files that correspond to observational data, and model output files at daily resolution at each grid after an ensemble run.
For each station, there are four different types of physics forcing data, including vertical profiles of eddy diffusive coefficients and temperatures, surface PAR, and atmospheric dust deposition.We already provided the relevant data for stations S1 and K2.The temporal resolution is 1 day for the vertical eddy diffusivity and 1 month for the three other types of data.

Code structure
All the source files including the makefile are stored in the src folder.Here we briefly describe the functions of the most important source files.
-Main.f90 is the main program for DRAM that calls each subroutine in serial.
-MOD_1-D.f90 is the major module that sets up and runs the 1-D model.The module also generates model output that matches the observational data.
-Interface_MOD.f90 is the module that initializes the absolute and normalized parameter vectors, the covariance matrix of the parameters, the prior parameter values, and the upper and lower parameter boundaries.
-SUB_MOD.f90 is the module that calculates the sum of squared errors (SSqEs) between model outputs and observational data.This module also contains the I/O subroutines that save the parameters, σ , and SSqEs for each iteration.It also contains the major subroutine MCMC_adapt that determines whether to accept new parameters, updates covariance matrix, proposes new parameter vectors, and calls the subroutine that runs the 1-D model with the newly proposed parameters.
-choose_model.f90 is the subroutine that defines the number and indices of tracers and the model outputs that need to be written into the output file.
-NPZD_cont.f90 is the major biological subroutine for the CITRATE model.
-bio_MOD.f90 is the module that declares most of the model names, indices for model input and output variables, and parameters.

Figure 1 .
Figure 1.Schematic description of the CITRATE model.Thick arrows indicate nitrogen flows and dashed lines indicate the simplified iron cycle.The inset denotes an example of a phytoplankton community with a lognormal distribution for cell volume.

Figure 2 .
Figure 2. (a) Locations of the three stations K2, S1, and ALOHA overlaid on annual Chl a climatology of the North Pacific.(b)-(e) Seasonal forcing of vertical eddy diffusivity (K v ), temperature, surface PAR, and atmospheric dust deposition, respectively, at station S1.The white squares are measured mixed layer depths from in situ temperature and salinity profiles.The thick tan line represents mixed layer depths calculated from a threshold of 10 −4 m 2 s −1 .(f)-(i) The same as (b)-(e), but for station K2. (j)-(m) The same as (b)-(e), but for station ALOHA.

Figure 3 .
Figure 3.An example of modeled patterns of total inorganic nitrogen (DIN), Chl a (Chl), mean size, and ln size variance for 4 years at stations K2 (a-d) and S1 (e-h).

Figure 5 .
Figure 5.Time evolution of fitted model parameters.Note that K P is the microzooplankton grazing half-saturation constant.

Figure 6 .
Figure 6.Model fittings to vertical profiles of (a)-(d) DIN, (e)-(h) Chl, (i)-(l) NPP, and (m)-(p) PON for four seasons at station K2.Black dots represent observational data and thick red solid lines represent the averaged seasonal values predicted by the model.Thin dashed lines represent the 95th percentiles of the seasonal data.

Figure 12 .
Figure 12.(a) Scatterplots of size variance versus phytoplankton community growth rate (µ com ).(b) Contributions of various factors to the dynamics of size variance in surface waters at S1.The term "competition" equates to v 2 d 2 µdl 2 .MIC and MES grazing equates to −v 2 d 2 g i dl 2 , "d4µ / dL4" equates to v 2 u d 4 µ dl 4 , and "trait diffusion" equates to 2 uµ.All the derivatives are evaluated at the mean size."Diffusion" means the contribution to the changes in size variance induced by diffusion with the underlying grid."Net effect" means the sum of the above terms."Net changes" mean the difference in size variance between adjacent days.(c) The same as (b), but at station K2.

Figure 14 .
Figure 14.Schematic diagrams for the mixing of two phytoplankton communities with different biomass, mean size, and size variance, each following a lognormal size distribution.

Table 1 .
Fixed parameters of the CITRATE 1.0 model.

Table 2 .
Parameters optimized by the DRAM algorithm.The values inside the parentheses of the initial values indicate the "hard" boundaries for the parameters.The numbers inside the parentheses of the optimized values indicate the standard deviation after the first 10 000 iterations have been removed.

Table 3 .
Observational data at stations S1 and K2.N : number of observations.Min and Max are minimal and maximal values used in data normalization (see Sect. 2.4 for details).DIN: dissolved inorganic nitrogen (µmol L −1 ).Chl a: total chlorophyll a concentration (µg L −1 ).NPP: net primary production measured by13C uptake (µg C L −1 d −1 ).PON: particulate organic nitrogen (µmol L −1 ).Fer: dissolved iron concentration (nmol L −1 ).SF Chl: percentages of four size fractions of Chl a.Note that the data of fer were from model outputs ofAumont et al. (2003)instead of real observations.

Table 4 .
Sum of squared errors between model outputs and observational data for sensitivity analysis.The standard run used the optimized parameter values in Table2.In other runs, only the value of the parameter shown was changed, while others were kept constant.