Ecological Phytoplankton acclimation to changing light intensity in a turbulent mixed layer: A Lagrangian modelling study

A new individual-based plankton model is used to test the hypothesis that the timescale of photoacclimation of phytoplankton within the surface mixing layer of the ocean is slow relative to mixing, in which case the chlorophyll to carbon (Chl:C) ratio of individual cells shows little adjustment in response to changes in light environment driven by vertical displacement. Rates of photoacclimation are shown to be a strongly non-linear function of light intensity that depends on the balance of intrinsic chlorophyll synthesis at low irradiance versus increasing growth rate at high irradiance. Predicted photoacclimation was negligible for cells experiencing rates of turbulent mixing typical of the open ocean surface boundary layer (10 − 3 to 10 -1 m 2 s -1 ), in which case Chl:C is set by mean light intensity. The model was extended to incorporate a simple ecosystem of nutrient, phyto- plankton, zooplankton and detritus and, using two-layer slab physics, used to study photoacclimation in a more realistic setting, the seasonal cycle of plankton dynamics at Ocean Weather Station India in the North Atlantic (59°N, 20°W). Results were remarkably similar when compared with an equivalent ecosystem model that used an Eulerian representation of phytoplankton, reinforcing our conclusion that mixing rates within the surface mixed layer of the ocean are typically too fast to permit photoacclimation by phytoplankton to ambient light.


Introduction
Primary production by phytoplankton drives the ocean food web, leading to export and sequestration of carbon in deep waters (Lutz et al., 2007). Much of this production occurs in the surface mixing layer, the thickness of which can vary seasonally from only a few metres during summer stratified conditions, to greater than 500 m over winter in some parts of the ocean. However, production can only occur within the euphotic zone, which typically extends to an average of 50 m in the North Atlantic (Ducklow, 1999). Turbulence, with velocities typically 5−10 mm s −1 , gives rise to rapid (time scales of minutes to hours) vertical movement of water in the surface mixing layer, e.g. as shown by the trajectories of neutrally buoyant floats (Harcourt et al., 2002(Harcourt et al., , 2010D'Asaro et al., 2014). It also means that phytoplankton cells experience vertical displacements and thereby a fluctuating light environment given the exponential attenuation of irradiance with depth (Gardner et al., 1999), with potentially important consequences for primary production.
Phytoplankton acclimate in response to changing light intensity, leading to variation in chlorophyll to carbon (Chl:C) ratios. Chlorophyll synthesis is maximised under low irradiance in order to improve light harvesting and photosynthetic efficiency, whereas there is no such need to produce high levels of light-harvesting pigment when exposed to saturating levels of irradiance (van de Poll et al., 2006). The extent to which cells become fully acclimated depends on the relative timescales of mixing and chlorophyll synthesis and it may be the case that phytoplankton do not undergo significant photoacclimation in the surface layer of the ocean given the high mixing rates. Previous studies, both observational (Lewis et al., 1984;Anning et al., 2000) and modelling (Falkowski and Wirick, 1981;Franks and Marra, 1994;McGillicuddy, 1995), have differed in their conclusions on this matter, which has important consequences for the calculation of depth-integrated photosynthesis within the mixed layer of the ocean. Individual-based (Lagrangian) modelling approaches should be most appropriate for calculating primary production if photoacclimation is prevalent, in which the phytoplankton assemblage is represented as the sum of individual cells, each with its own unique position in the water column and associated life history, including the photophysiological response (Ross et al., 2011a). Standard Eulerian models should otherwise suffice, in which phytoplankton are treated as a single state variable with uniform Chl:C.
Here, we present a new Lagrangian model and use it to calculate primary production in the surface mixing layer of the ocean, focusing on the associated sensitivity to photoacclimation as seen in the predicted life histories of individual cells in a turbulent environment. Two model setups are examined. The first incorporates an individual-based representation of phytoplankton into an ecosystem model, embedded within a 2 vertical layer slab physics in which the depth of the upper (surface mixed) layer, as well as surface irradiance, vary seasonally. It is applied to a realistic ocean setting, namely the seasonal cycle at Ocean Weather Station India (OWSI) in the northern North Atlantic (59°N, 20°W), using turbulence forcing based on observed winds. The second model version provides a theoretical analysis of the interaction between photoacclimation and rate of mixing by representing phytoplankton as the sole biological entity, embedded within a mixing layer of fixed depth and exposed to constant irradiance. Using these two analyses, our aim is to test the following hypothesis: the timescale of photoacclimation of phytoplankton within the surface mixing layer of the ocean is slow relative to mixing such that individual cells undergo little adjustment of their Chl:C ratio in response to changes in light environment driven by vertical displacement. If true, this would mean that most, if not all, cells within the mixing layer would exhibit approximately uniform Chl:C ratios.

Model description
We first introduce the Eulerian model, detailing all of the equations and the processes they represent. The two Lagrangian model setups, i.e. the full ecosystem model, and the simplified fixed slab model, are then described in turn. Lists of model variables and parameters are provided in Tables 1 and 2, respectively.

The Eulerian model
The Eulerian model assumes that turbulence occurs within a single surface layer, and that mixing is sufficient to ensure all components are fully mixed at all times. Concentrations of nutrients, phytoplankton, zooplankton, and detritus (NPZD) are therefore necessarily homogeneous throughout this surface mixing layer. Primary production is calculated numerically, using an average value for Chl:C throughout the surface mixing layer, but allowing changing light intensity with depth (at 1 m resolution).
Phytoplankton biomass is represented by two state variables that permit flexible stoichiometry in terms of nitrogen (N) and chlorophyll (Chl) content: P (mmol N), and Chl (g Chl). Carbon biomass is calculated assuming a fixed Redfield ratio, parameter ζ, of 12.6 mmol N g C −1 (equivalent to 6.625 mol C mol N −1 ). Photosynthesis, V (g C g Chl −1 h −1 ), is described by a Smith function: where V max is the maximum chlorophyll-specific rate of photosynthesis (g C g Chl −1 h −1 ), α is the initial slope of the photosynthesis-irradiance curve (g C g Chl −1 h −1 (W m -2 ) −1 ) and I(z) is irradiance at depth z (W m -2 ). A constant irradiance, I 0 is imposed at the ocean surface which is then attenuated with depth according to Beer's law: The vertical attenuation coefficient, k PAR , depends on chlorophyll concentration in a piecewise calculation in the water column: 0-5 m, 6-23 m, and > 23 m (Anderson, 1993;Anderson et. al., 2015). Attenuation is calculated in 1-metre intervals based on the summed chlorophyll content in each interval. The specific growth rate, μ P (d −1 ), is then: where θ is chlorophyll to carbon ratio (g Chl g C −1 ; equal to ζChl/P). ρ Chl (g Chl mmol N −1 ) is the scaling factor which specifies instantaneous chlorophyll synthesis relative to N and is an inverse function of irradiance such that phytoplankton acclimate at low light intensities by producing extra chlorophyll: where θ max is the maximum achievable chlorophyll to carbon ratio. This equation is based on the dynamic photoacclimation model of Geider et al. (1997), which predicts the dependencies of Chl:C and growth rate on irradiance, temperature, and nutrient availability. The equation for the rate of change of N biomass is:  [1] Adjusted to fit the seasonal cycle of chlorophyll as OWSI (Fig. 2). References: FDM90 Fasham et al. (1990); Y13 Yool et al. (2013); A15 Anderson et al. (2015).
Growth (the first term) depends on photosynthesis and nutrients. The maximum rate of photosynthesis, V max , is temperature-dependent according to Eppley (1972): where V max 0 is V max at 0°C and k N is the half saturation constant for nutrient uptake. There are three biological loss terms in Eq. (5), namely grazing and linear and non-linear non-grazing losses (Anderson et al., 2015). The linear term represents a respiration rate of 0.02 d −1 (parameter m P2 ). The non-linear term represents density-dependent processes such as viral lysis and takes a Michaelis-Menten form governed by parameters m P3 (maximum rate, mmol N −1 d −1 ) and k P (the half saturation constant, mmol N m -3 ). This loss term is calculated as a function of the average phytoplankton concentration in the mixing layer, P (mmol N m -3 ). A Holling Type III equation is used for grazing: where I max is the maximum grazing rate (d −1 ) and k g is the half-saturation constant for grazing (mmol N m -3 ). Phytoplankton are lost from the surface layer to the layer below via a constant mixing rate, k mix (d −1 ), as well as detrainment as the mixing layer shoals (term h ' in Eq. (5)). The corresponding equation for the chlorophyll content of phytoplankton is: The equation for zooplankton is (Anderson et al., 2015): where β is absorption efficiency, k Z is net production efficiency and m Z is the mortality rate (mmol N −1 d −1 ). Zooplankton are subject to mixing, as well as detrainment and dilution due to changes in the depth of the mixing layer (term h ' ). The equation for detritus is: There is also an input from zooplankton mortality, and detritus is subject to a temperature-dependent (Eppley function: Eq. (11)) remineralisation rate, m D (d −1 ). As well as losses due to mixing and changes to the depth of the mixing layer, detritus is also subject to a sinking rate of 10 m d −1 (parameter V D ).
Finally, the equation for dissolved inorganic nitrogen is: A corresponding Lagrangian model was formulated based on the same equations, designed to investigate two scenarios: firstly, a 'fixed slab', phytoplankton-only scenario, for a detailed analysis of photoacclimation in the model in response to mixing rate and light intensity and, secondly, development of a full ecosystem model to provide a direct comparison between the Eulerian and Lagrangian formulations.

Fixed slab model
Phytoplankton are the only biological entity in this version of the model and are represented as "super-individuals" (Scheffer et al., 1995) within a vertical slab of water with a fixed depth, H (m). The slab is subject to a constant mixing rate and constant irradiance at the surface. It is computationally impossible to represent each individual phytoplankton cell in a real ocean water column, hence the use of superindividuals, each of which represents a group of cells. Each super-individual has its own unique vertical location in the water column which changes over time according to turbulent mixing and particle movement rules.
Now, I i becomes the irradiance experienced by a super-individual at a given point in time (W m −2 ). The specific growth rate of a superindividual, μ P,i (d -1 ), is then: The equations for rate of change of N and Chl biomass in superindividuals are now: Mortality is calculated as a simple quadratic function with parameter m P = 0.3 (mmol N) −1 d −1 . This simple mortality function does not explicitly represent processes such as grazing, but acts as a control on phytoplankton population as our main focus is rates of phytoplankton growth and acclimation, rather than the full ecosystem. Superindividuals are transported passively in the water column by turbulent motion. A simple random walk model is used based on the specification of turbulence via a diffusivity coefficient, K v (Falkowski and Wirick, 1981;Visser, 1997;Ross and Sharples, 2004). Turbulent mixing is calculated using the non-local K-profile Parameterization (KPP) (Large et al., 1994). The eddy diffusivity K v is calculated from: where S is the vertical shape function and the vertical velocity scale of turbulence, w (m s −1 ), is calculated from: where κ = von Kármán's constant (0.4) and u* is the wind friction velocity approximated by u* = 0.001U (Large, 1994), where U is surface wind speed (m s −1 ). Here, we use a single value for S, therefore assuming a constant turbulence throughout the mixing layer. This simplification enables a more straightforward analysis of the model results. The equation for the random walk is (Visser, 1997): Here, z n is current depth, z n+1 is depth at the next time step and R represents a random process, with mean of 0 and variance r (we use a uniform distribution for R between +1 and -1, such that r = 1/3). Note that this equation is only strictly valid if K v is assumed to remain constant throughout the displacement from z n to z n+1 . The distance moved each time step is therefore a function of surface wind speed and depth of the mixing layer. The model has reflecting boundary conditions such that super-individuals bounce off the surface and base of the water column. The constant fixed depth in the fixed slab model means that, unlike the mixing layer model with ecosystem (below), there is no detrainment or dilution of super-individuals. The model was run to steady state using a time step of 1 min. It was initialised with 20,000 super-individuals, each with P i = 0.015 mmol N, Chl i = 0.012 mg Chl and Chl:C ratio of 0.01 g Chl g C −1 . Model results were shown to be insensitive to the initial Chl:C ratio and so a midrange value was chosen.

Mixing layer model with ecosystem
This version of the model has two layer slab physics (Anderson et al., 2015) in which a seasonally-varying upper mixing layer that incorporates an ecosystem made up of nitrate-phytoplankton-zooplankton-detritus (NPZD) is positioned above a deep layer that contains only nitrate in a fixed vertical profile (Fig. 1). State variables are exchanged between the two layers via mixing across the base of the mixing layer, entrainment and detrainment as the mixed layer deepens or shallows, and by sinking of detritus. Z, N and D are defined as standard Eulerian state variables representing average concentrations in the surface layer, i.e., with a single differential equation for each. Phytoplankton, on the other hand, are represented in Lagrangian mode, with state variables P i and Chl ,i . The average phytoplankton nitrogen biomass per metre (P conc ), calculated by summing across all P i and divided by the depth of the surface layer, is then used as the input for zooplankton, detritus, and dissolved inorganic nitrogen. The equations for P i are: Super-individuals are lost from the surface layer to the layer below via a constant mixing rate, k mix (d −1 ), as well as detrainment as the mixing layer shoals. The same reflecting boundaries were used as in the fixed slab model, such that super-individuals bounce off the surface and base of the water column.
The corresponding equation for the chlorophyll content of superindividuals is: In order to compensate for mortality and maintain a representative number of super-individuals in the water column, super-individuals divide (split in two) when they reach a threshold size (P div , mmol N). This use of splitting is solely a means of controlling the number of particles within the simulation and is not related to the biological process of cell division.
The equation for zooplankton is (Anderson et al., 2015): The equation for detritus is: Finally, the equation for dissolved inorganic nitrogen is: Fig. 1. The structure of the ecosystem (adapted from Ross and Sharples, 2007) The slab physics and forcing were set up for a climatological scenario for Ocean Weather Station India (OWSI) in the northern North Atlantic (59°N, 20°W). This is an ideal site for testing the model due to the strong seasonality in the depth of the surface mixing layer and characteristic spring bloom. The annual cycle of mixing layer depth was taken from World Ocean Atlas (Antonov et al., 2010). Irradiance at the ocean surface was calculated using standard trigonometric equations (Brock, 1981;Anderson et al., 2015) in conjunction with attenuation by clouds according to Reed (1977). Photosynthetically active radiation (PAR) as a fraction of the total was 0.43 (Anderson et al., 2015) and, as for the fixed-slab model, attenuation in the water column is a piece-wise approximation (Anderson, 1993;Anderson et al., 2015). Irradiance varies throughout the day according to a sinusoidal function. Daily wind speeds were used to calculate turbulent mixing (Eq. (19)), based on a sinusoidal function fitted using non-linear least squares to an annual cycle obtained by averaging data from 1997 to 2001, taken from the ECMWF ERA-40 atlas (Uppala et al., 2005). Finally, nutrient was assumed to be present in the lower model layer as a linear gradient (Anderson et al., 2015): where N(z) is nutrient at depth z. The regression coefficients ( aN = 0.0074, b N = 10.85) were fitted for OWSI from World Ocean Atlas data (Levitus et al., 2010). Model parameters were mostly taken from previous modelling studies, including those focusing on the North Atlantic (Fasham et al., 1990;Yool et al., 2013;Anderson et al., 2015). The magnitude of acclimation to light in the model phytoplankton is controlled by the ratio θ θ / max i and thus θ max , the maximum achievable Chl:C ratio, is an important parameter although assigning a value is not straightforward. Maximum ratios of Chl:C can range over an order of magnitude, from ∼0.007 to over 0.07 g Chl (g C) −1 (Geider et al. (1997). Previous models of photoacclimation have focused only on single species. Flynn et al. (2001), for example, used data for the marine diatom Skeletonema costatum when comparing different models of photoacclimation in phytoplankton and estimated values of between 0.054 and 0.079 g Chl (g C) −1 . Oceanic phytoplankton comprise a multitude of species and we therefore tuned this parameter for OWSI, resulting in a value of 0.02 g Chl (g C) −1 . This value is consistent with satellite-derived values of Chl:C ratio in the open ocean (Behrenfeld et al., 2005) that are considerably lower than those seen in laboratory studies. The threshold for particle division, P div , was chosen so as to ensure that the number of super-individuals in each simulation was greater than 2000.

Ocean Weather Station India
The model was first run for OWSI to examine photoacclimation in a realistic setting, comparing results for Lagrangian and Eulerian representations of phytoplankton. This analysis involved tuning the value of parameter θ max , providing a value of 0.02 g Chl g C −1 that is representative of open ocean phytoplankton and which is subsequently used in the theoretical analysis of the interplay between acclimation and the timescales of mixing (Section 3.2). The model was compared with chlorophyll data which are 8-day averages for a characteristic year, 2006 (Anderson et al., 2015) and with climatological nitrate data from World Ocean Atlas (Antonov et al., 2010). The predicted seasonal cycles for these two variables, using the Lagrangian representation of phytoplankton and θ max = 0.02 g Chl g C −1 , is shown in Fig. 2 (simulation length was 3 years, with the third year shown by which time a repeating annual cycle occurs). The general agreement is good, similar to other slab modelling studies of this site (Fasham, 1995;Anderson et al., 2015). The timing of the spring chlorophyll bloom is about a month late as compared to the data which may be a consequence of the climatological mixing layer forcing, whereas the data are for a particular year. The underestimation of overwintering stock of phytoplankton at high latitudes is a common feature of slab models (Fasham, 1995;Anderson et al., 2015) which may be due to the assumption of complete mixing within what is a deep mixing layer. The predicted timing and magnitude of the seasonal nitrate drawdown are good, which is encouraging given that the nitrate data themselves represent an average climatology.
The Eulerian version of the model generates results that are nearly identical to those when phytoplankton are modelled in a Lagrangian manner (Fig. 3). Predicted Chl:C varied between about 0.012 and 0.17 g Chl (g C) −1 in both model versions, higher in the winter months because of low irradiance, but not showing any signs of a diel cycle. The fact that the Lagrangian and Eulerian versions of the model generated predictions that match closely for phytoplankton, nutrients and Chl:C ratio suggests that the timescales of mixing are too fast for photoacclimation to occur, a phenomenon that is examined in more detail in the next section.

Fixed slab model
The interaction between rates of mixing and photoacclimation was next investigated using the fixed slab model with phytoplankton as the sole state variable (no explicit ecosystem or limitation by nutrients), with parameters as for OWSI (notably, θ max = 0.02 g Chl g C −1 ). Six simulations were carried out, each with 20,000 super-individuals in a fixed slab of 100 m and with constant forcing in terms of irradiance at the ocean surface (200 W m -2 ). Wind speed was varied for the different model runs, between 0 and 100 m s −1 , generating turbulent mixing, K v , of between 0 and 0.4 m 2 s −1 . The timescales (T m ) associated with the turbulent mixing in the surface layer can be approximated using the layer depth (H) and the turbulent diffusivity (K v ) (Taylor and Ferrari, 2011), and are given in the legend for Fig. 4: The simulations were run to steady state and snapshots taken that show Chl:C versus incident irradiance for 10,000 randomly chosen super-individuals (Fig. 4). A sensitivity analysis was carried out to investigate how the predicted spread of Chl:C values was affected by changes to θ max and V max (default values ± 10 %). The results (not shown) indicate that, while changes to these parameters affect the fully acclimated Chl:C ratio for super-individuals at each irradiance, the spread of ratios shown over the whole population showed low sensitivity to changes in the parameterisation of θ max or V max . Fully acclimated Chl:C, θ A , occurs when there is zero mixing such that superindividuals have fixed depths, distributed evenly throughout the 100 m slab (Fig. 4a). The red lines in Fig. 4 show θ A : Acclimated Chl:C thus decreases with increasing irradiance, as expected. Photoacclimation was, however, incomplete when even a low rate of mixing is introduced (4 × 10 −5 m 2 s -1 ; Fig. 4b). Higher Chl:C ratios were seen in super-individuals experiencing low irradiance, but there is a large spread about the red line. This spread of points in the plots diminishes as the rate of mixing is increased although not towards the fully acclimated state but, rather, towards uniformity whereby, at the highest mixing rate of 4 × 10 -1 m 2 s -1 , all the super-individuals ended up with much the same Chl:C ratio close to 0.011 g Chl g C -1 (Fig. 4f). Super-individuals were unable to photoacclimate to incident irradiance in the fluctuating light environment at high mixing. Rather, they acclimated to the average irradiance in the euphotic zone (surface to depth of 1 % irradiance; red spots in Fig. 4).
Although irradiance determines fully acclimated Chl:C, as described above, the rate of chlorophyll synthesis (photoacclimation) depends also on growth rate, μ P , which is an increasing function of irradiance: where μ Chl sp is the specific rate of chlorophyll synthesis (d −1 ). This relationship is shown in Fig. 5, with data points from Fig. 4a superimposed. There is thus an optimum irradiance, in this case 32 W m -2 , at which the rate of chlorophyll synthesis, i.e., photoacclimation, is maximised, with μ Chl sp declining on either side of the optimum. Photoacclimation is greatest when super-individuals experience sustained exposure to irradiances that maximise μ Chl sp . The key factor is the timescales involved, i.e., to quantify what is meant by "sustained" in this regard. Working from the simulations shown in Fig. 4, we identified, for each mixing rate, the super-individual that experienced the closest average irradiance to 32 W m -2 during the last 5 days and show its trajectory in terms of Chl:C, depth and irradiance during this period (Fig. 6a). Zero mixing leads to fully acclimated Chl:C is 0.014 (Eq. 19), the super-individual occupying a depth of 6.6 m. When mixing is low (≤ 4 × 10 -4 m 2 s −1 ), the super-individuals in question maintain a depth of between 3 and 10 m throughout the 5-day period, with an average depth of 7.0 m. An increase of mixing to 4 × 10 -3 m 2 s −1 , however, means that irradiance is no longer sustained over several days and photoacclimation is minimal. Higher mixing rates lead to rapidly fluctuating vertical position and irradiance, resulting in little or no photoacclimation. Similar results are seen when a lower average irradiance, 10 W m -2 , is selected (Fig. 6b). Photoacclimation is seen for mixing rates ≤ 4 × 10 -4 m 2 s −1 , with super-individuals maintaining a depth of between 4 and 15 m and Chl:C tending towards 0.019 (the acclimated value for 10 W m -2 ). In contrast, Chl:C tended towards 0.011 at higher mixing rates, corresponding to the average irradiance in the water column, 50.9 W m -2 .

Discussion
Rates of photoacclimation are a strongly non-linear function of light intensity, peaking at an optimum irradiance that maximises the balance  between increasing Chl:C at low irradiance with increasing growth rate (and associated synthesis of chlorophyll) at high irradiance (Fig. 5). A new Lagrangian model was used to test the hypothesis that rates of photoacclimation in the surface mixed layer of the ocean, due to changing irradiance as a consequence of vertical displacement, are slow relative to mixing in which case little variation occurs in Chl:C between phytoplankton cells in the water column. The hypothesis was supported by our model analyses. The theoretical results of the fixed slab model indicate that significant photoacclimation only occurs when mixing rates are ≤ 4 × 10 −4 m 2 s -1 , which corresponds to an approximate mixing timescale of nearly 30 days (phytoplankton Chl:C ranged from 0.003 to 0.02 g Chl g C -1 ). These slow mixing rates allowed phytoplankton cells to remain at the depth of optimum irradiance for a period of several days or more. Mixing rates in the ocean surface boundary These studies imply diffusivities ranging from (10 -3 to 10 -1 m 2 s -1 ) in the surface boundary layer and that phytoplankton experience rapid vertical displacements on time scales of hours. These magnitudes and timescales of displacement are evident immediately in the trajectories of Lagrangian floats as shown in D' Asaro (2001). When mixing rates of this order were used in our model (≥ 4 × 10 -2 m 2 s -1 , corresponding to an approximate mixing timescale of 7 h), results showed that phytoplankton photoacclimate to average light intensity in the euphotic zone, with Chl:C ratio showing a relatively narrow range between 0.008 and 0.017 g Chl g C -1 . Moving beyond the theoretical analysis, we also examined the performance of the Lagrangian phytoplankton model as part of an NPZD ecosystem model embedded within 2-layer slab physics, run for Station India in the northern North Atlantic. The rate of turbulent mixing in the surface mixing layer was determined from observed wind speeds, varying between 0.01 and 0.29 m 2 s −1 . Phytoplankton growth was driven by nutrients and seasonally varying irradiance at the ocean surface. Results were compared with an equivalent Eulerian model in which phytoplankton were treated as a single homogeneous state variable, with a single Chl:C ratio calculated from depth-averaged light intensity. Predicted seasonal cycles of N, P, Z and D were remarkably similar for the two model versions, showing the characteristic bloom of phytoplankton and associated nitrate drawdown (Fig. 2). Chl:C also exhibited strong seasonality, with highest values in winter when incident irradiance is low and the mixed layer is deep. Predicted Chl:C was remarkably similar between the Lagrangian and Eulerian model versions, being lower in the former, but only marginally so (Chl:C varied 0.0103 -0.0165 and 0.0109 -0.0166 g Chl g C −1 in the   Tomkins, et al. Ecological Modelling 417 (2020) 108917 Lagrangian and Eulerian versions of the model, respectively (Fig. 3)). These marginal differences are not a result of individual photoacclimative properties, but rather result from differences in the growth rates. The mixing rates in the Lagrangian model are not always fast enough to completely overcome the growth rates in the upper part of the mixing layer, resulting in a slight gradient in the phytoplankton concentration. Enforcing the vertical phytoplankton profile predicted by the Lagrangian model onto the Eulerian model results in a perfect match between their respective predictions (not shown). However, the aim of the study was not to just test when Eulerian and Lagrangian models diverge for all scenarios, which would have required the addition of a diffusive term to the Eulerian model, but to test when they diverge in the special case of a well-mixed surface layer. In this case, we find that divergence occurs at mixing rates characteristically much lower than those of the surface mixing. The fact that the two model versions produce the same results at the high mixing rates that are characteristic of the surface mixing layer of the ocean (10 -3 to 10 −1 m 2 s −1 ) is sufficient for us to conclude that mixing in the natural environment is high enough to prevent variation in Chl:C (from photoacclimation in response to ambient light). These results re-emphasise those of the theoretical analysis, indicating that photoacclimation to ambient irradiance is of little consequence in the surface mixed layer of the ocean because rates of mixing are too fast. Furthermore, they indicate that, at least in terms of areas of the ocean where primary production occurs primarily in the surface mixed layer, there may be limited benefit in using a Lagrangian representation of phytoplankton in models, rather than a standard Eulerian approach. Of course, photoacclimation does occur on longer timescales as mean light intensity changes, e.g., seasonally or due to changing mixed layer depth. Our results are in agreement with field studies that have examined timescales of photoacclimation. Lewis et al. (1984) took samples from the surface and near bottom in the Bedford Basin (a coastal inlet of Canada's Atlantic coast) and found no evidence of photoacclimation when mixing occurred on time scales less than ∼10 h. Anning et al. (2000) subjected cultures of the marine diatom Skeletoneme costatum to irradiance shifts over 15 days. Results showed that it took 3-5 days for the diatoms to fully acclimate to changes in irradiance, with the chlorophyll content of the cells changing 2.5 fold. Our results likewise show that acclimation takes several days in order that significant changes in Chl:C occur, and that such changes are negligible on timescales of less than a day (Fig. 6). Previous modelling studies have been equivocal about the interaction between rates of mixing and photoacclimation, and the resulting consequences for primary production. Falkowski and Wirick (1981) found that mean and variance of Chl:C ratios in phytoplankton, along with photosynthetic rates, showed little sensitivity to mixing rates ranging from 0 to 0.001 m 2 s −1 when using a random walk Lagrangian model within an idealised 20 m mixed layer. On the other hand, Franks and Marra (1994), using a similar idealised model with a 15 m mixed layer, found that depth-integrated photosynthesis increased with increasing mixing. Their model included photo-inhibition which may have negatively impacted on cells near the surface, especially at low mixing rates. McGillicuddy (1995) compared the performance of two photoacclimation models, Wolf and Woods (1988) and Lande and Lewis (1989), at two mixing rates (0.01 m 2 s −1 , and 4.64 m 2 s −1 ) and found little or no difference in predicted phytoplankton growth rate for the two model implementations under weak mixing, whereas a significant reduction in mean growth rate occurred in the Wolf and Woods (1988) model under strong mixing. The advantage and novelty of our approach is two-fold. First, we used the fixed slab model to examine a wide range of mixing rates, both representing those in the ocean mixed layer and slower rates, using a contemporary model of photoacclimation (Geider et al., 1997) and a maximum Chl:C (parameter θ max ) appropriate for ocean phytoplankton (obtained by tuning). Second, we modelled a full seasonal cycle for a realistic ocean setting (Ocean Weather Station India), using a complete NPZD ecosystem model and seasonal forcing for wind and irradiance, and compared the performance of Lagrangian and Eulerian formulations for this scenario.
An interesting application of Lagrangian models has been to study whether measurements of primary production from bottles suspended in the water column, in which algae experience more or less constant irradiance to which they acclimate, are compromised given that cells in the surrounding water experience vertical mixing and thereby fluctuating irradiance. Barkmann and Woods (1996) constructed and parameterised a 1D model for this purpose and found that the primary production calculated of freely moving phytoplankton at a given depth in the water column was 40 % lower compared to cells that were confined to the same depth without being able to move. Similar results were found by Ross et al. (2011a, b) when undertaking a Lagrangian modelling experiment in which calculated growth rates were compared for phytoplankton cells subject to mixing versus cells that were "frozen" for 24 h, for two ocean settings: shallow, coastal regions (surface mixed layer =6 m) and the open ocean (surface mixed layer =60 m). Our results suggest that photoacclimation, leading to significant changes in Chl:C, takes place over several days. Incubation experiments to measure primary production are usually carried out over 24 h (Mingelbier et al., 1994), in which case errors related to photoacclimation may not be large.
Our findings do not, however, mean that there is no place for Lagrangian approaches to modelling phytoplankton and their role in marine ecosystems. Previous studies have shown that models formulating growth as a nonlinear function of an internal nutrient quota (eg. Droop kinetics) give differing predictions when represented in a Lagrangian, as opposed to Eulerian formulation (Hellweger and Kianirad, 2007;Baudry et al., 2018). Additionally, much of the primary production in the ocean does not occur in the surface mixed layer but, rather, in stratified waters such as the deep chlorophyll maximum and elevated chlorophyll content may be expected as a consequence of the low irradiances at depth (Cullen, 1982). Individual differences in functional traits may be important when it comes to understanding phenomena such as the deep chlorophyll maximum and ephemeral features such as blooms caused by eddies, in which case the spatial patterns and dynamics may not be fully captured using traditional Eulerian techniques (DeAngelis and Yurek, 2017). Moreover, Lagrangian models can be integrated with geographical information systems (GIS) data, such as satellite observations of sea surface temperature, in order to better understand the forces driving phytoplankton evolution and which determine ecological niche. The results of our study nevertheless indicate that standard Eulerian approaches should be more than sufficient for many contemporary marine ecosystem model applications.

Conclusions
The results of our individual-based (Lagrangian) modelling study indicate that photoacclimation in phytoplankton is negligible for cells experiencing rates of turbulent mixing that are typical of the open ocean surface boundary layer (10 -3 to 10 -1 m 2 s −1 ). The individualbased representation of phytoplankton was incorporated with an ecosystem model that also contains nutrient, zooplankton and detritus and used to simulate the seasonal cycle of plankton dynamics at Ocean Weather Station India in the North Atlantic (59°N, 20°W). Results were compared with a separate version of the model in which all state variables, including phytoplankton, were represented in Eulerian mode. The two simulations showed remarkable similarity, again indicating that rates of mixing in the surface ocean are too fast to permit photoacclimation of phytoplankton to ambient light. There may therefore be little or no benefit in using Lagrangian models to calculate primary production, rather than a standard Eulerian approach, at least for ocean sites where production is predominantly in the surface mixed layer.

Declaration of Competing Interest
None.