An expanded Eulerian model of phytoplankton environmental response

Ecological Modelling 118 (1999) 237–247

An expanded Eulerian model of phytoplankton environmental response Gerald S. Janowitz *, Daniel Kamykowski Department of Marine, Earth and Atmospheric Sciences, Box 8208, North Carolina State Uni6ersity, Raleigh, NC 27695 -8208, USA Accepted 25 January 1999

Abstract An Eulerian approach to modeling plankton physiological responses to environmental factors is developed wherein the time history of cell exposure to two external environmental fields over specified time intervals are utilized as independent variables along with position and time to help characterize the cell population. We seek to find the concentration of cells per unit volume as a function of depth, time, and the time histories of exposure to PAR (photosynthetically active radiation) as it influences internal cellular carbon through phototsynthesis and to nitrate as it influences internal cellular nitrogen through nutrient assimilation. The response under consideration here, vertical swimming, is taken to depend on historical exposure to the external PAR and nitrate fields. The model can be readily extended to other external fields and to more than the one historical time scale here associated with each external field. This type of model joins Lagrangian models as most beneficial when phytoplankton physiology responds to environmental factors in a nonlinear fashion, i.e. when the mean response does not depend on the mean exposure. A simple example is discussed and the impact of wind-driven mixing is explored. © 1999 Elsevier Science B.V. All rights reserved. Keywords: Eulerian; Model; PAR; Nitrate; Swimming

1. Introduction The natural way of modeling the response of a plankton cell to its environment is to track the response of the cell to a range of external fields to which the cell has been exposed over time as it moves through the environment, a Lagrangian * Corresponding author. Tel.: + 1-919-5157837; fax: +1919-5157802. E-mail address: [email protected] (G.S. Janowitz)

description. Woods and Onken (1982) popularized this approach that subsequently was applied to a variety of problems including the deep chlorophyll maximum (Wolf and Woods, 1988), phytoplankton motility in the upper mixed layer (Yamazaki and Kamykowski, 1991), diatom demography (Woods and Barkmann, 1993), and phytoplankton photoresponse in the upper mixed layer (Franks and Marra, 1994; Kamykowski et al., 1994). Alternatively, traditional Eulerian models focus on predicting cell concentrations as a

0304-3800/99/$ - see front matter © 1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 4 - 3 8 0 0 ( 9 9 ) 0 0 0 3 7 - X

238

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

function of space and time variables and physiological responses of these populations via the use of mean properties of these cells, e.g. Janowitz and Kamykowski (1991). A review of other theoretical approaches to ecosystem modeling is given in Muller (1997). Lande and Lewis (1989), McGillicuddy (1995) explored the relationship between Lagrangian and Eulerian points-of-view in biophysical models. While many Lagrangian models do exist, it is difficult to compare field results to Lagrangian predictions except with the use of cell specific techniques like epifluorescence microscopy or flow cytometry (Kamykowski, 1995). Field observations made at points in space and based on the more common bulk population techniques (Kamykowski, 1995) will contain a large number of plankton with different time histories, i.e. the observations are Eulerian. The concentration of cells in some region in space at some time can be decomposed into categories of cells with different exposure histories. A straightforward Eulerian model poses the following question: if a population is present at some location at some time, what is the average exposure history of this population? Here we pose a more complicated question: at some location and time and set of exposure histories, what portion of the population represents each possible exposure history? The answer to this question can yield information on the central tendency, the variance, and higher order statistics of the local population’s exposure history. A straightforward Eulerian model is useful if the physiological responses in which we are interested depends only on the average exposure history of all the cells present at some location. However, if the physiological response of interest depends on more than one exposure history in a nonlinear manner, then an expanded Eulerian model developed here is more useful; by expanded we mean an increase in the number of independent variables utilized. In the next section, we demonstrate the construction of an expanded Eulerian model in a simple context. In what follows we use exposure histories along with location and time as independent Eulerian variables. We shall construct a simple one dimensional (in space) time-dependent

model using two exposure histories, one to PAR modeled as a cellular carbon controller and the other to nitrate modeled as a cellular nitrogen controller along with the vertical variable z and time as four independent variables. We seek to determine the concentration of cells as a function of these four variables and predict the response of cell concentration to a swimming behavior specified as a function of the exposure variables. Turbulent mixing, parmeterized via a constant eddy diffusivity, will also be included down to a specified depth below which a stable stratification in mass density is assumed present to eliminate turbulent mixing. This simplified turbulent model is utilized to demonstrate the effects of winddriven mixing. Results of this mixed case will be contrasted to the unmixed, swimming only, case.

2. Model description The most important external fields which determine physiological responses include PAR, nutrients, temperature and salinity. Any response of a cell to these fields should depend not only on the present state of these fields but also on the past history of exposure of a cell to these fields. If Xi (t) represents any one of these external fields at time t, we define the exposure of the cell to this field over a time scale T, EXi, by EXi (t,T) =(1/T)0/t e( − (t − t%)/T)Xi (t%) dt% +EXi (0,T) e − t/T,

(1)

and hence, d EXi /dt =(Xi −EXi )/T.

(2)

The subscript i ranges from 1 to the number of relevant external fields.We note that this definition is equivalent to a first order kinetic response (Cullen and Lewis, 1988). We also note that using several different time scales for the same external field yields different exposure histories and that Xi (t) may be either a single external field or a function of one or more external fields. It is clear from Eq. (2) that if EXi (t,T) is known, we can determine the Xi (t). It is convenient, though not essential, to scale the exposure to a value between

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

0 and 1. If we know the minimum (EXimin) and maximum (EXimax) possible exposures, we can define the scaled exposure, EX*i = (EXi −EXimin)/(EXimax − EXimin).

(3)

Henceforth, we assume all exposures have been so scaled and drop the asterisk. Here we shall utilize only two external fields, the PAR and nitrate fields. We consider one spatial variable, the depth. The variable z is zero at the surface and decreases downwards (z B 0), with the vertical domain decomposed into M layers each of thickness dz; here M=50 and dz = 1 m. The integer k refers to the layer number, with fluid in layer k ranging from a depth of (k −1)*dz to k*dz below the surface (15k 5M). The external nitrate and PAR fields are specified as follows: N(z) =NIT/NITmax =0.5*(1.− tan h((z + ZN)/2),

(4a)

P(z,t)=PAR/PARmax = −e(z/ZP)sin (2*y*t/86 420 s),

(4b)

All negative values of P(z,t) are set to 0. The fields, N(z) and P(z,.75 days), are given in Fig. 1 with ZN= ZP=10 m. This simple concentration independent form of the PAR distribution is utilized to most simply study the effects of swimming behavior on the population distribu-

239

tion. We note that the nitrate field will be held fixed with time in our experiments. Our first experiment will have no wind mixing so this is reasonable, but our second will have mixing down to a specified depth. Although the time variation of the nitrate field in this mixing case could be readily calculated independently of the remainder of the calculation, we fix this field in time to evaluate the effects of mixing on the cell density distributions subject to the same nitrate and PAR fields. The nitrate field increases monotonically with depth. The PAR field is 0 for the first 12 h and positive for the next 12 and has a period of 1 day. The exposures to NIT and PAR, now denoted as EN and EP, can be calculated from Eq. (1). The minimum values of these variables are both zero and the maximum values are NITmax and PARmax/SC, where SC can be calculated by inserting equation Eq. (4b) at z=0 into Eq. (1), requiring that the value at t =1 day equals the value at t=0, and noting the maximum value of EP. This numerical calculation shows that SC ranges from 1 for T =0 to y for very large values of T and equals 2.453 for T=1 day. The equation governing the dimensionless EN and EP are thus: d EN/dt/(N(z)−EN)/TN,

(5)

and d EP/dt =(SC*P(z,t)−EP)/TP,

Fig. 1. Initial scaled PAR and nitrate–nitogen profiles. Initial values of nn and np are 10× value.

(5a)

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

240

Both EP and EN vary continuously from 0 to 1. For mathematical convenience we decompose the EN field into ten equal ranges specified by the integer nn which specifies a value of EN between (nn −1)/10. and nn/10. The EP field is broken up into ten parts specified by the integer np, though not of equal lengths. The value np = 1 specifies values of EP B0.01; np= 2 refers to 0.01 B EP B 0.05; np =3 refers to 0.05B EP B 0.125; and np 4 – 10 refers to 0 .125*(np − 3)B EP B 0.125*(np − 2). The length of each interval is specified by LN(nn)=0.1 for all nn; LP(1)= 0.01,LP(2)= 0.04,LP(3) =0.0725, and for np]4,LP(np)= 0.125. We utilize the nonuniform spacing for EP since most of the population will have values of EP B0.4 and this spacing gives better resolution for the weakly exposed population. Here we shall use TN=TP=1 day. The primary dependent variable in our model is the number of plankton cells per unit volume at a given depth and some time with exposures EN and EP, i.e. C(z, EN, EP, t), or in our discretized notation C(k, nn, np, t). The total concentration of cells in each layer, Ctotal(k, t), is the sum of C(k, nn, np, t) over np and nn from 1 to 10. We wish to compute how C and Ctotal change with time from some initial distribution, due to the effects of time and depth varying EP and EN, a vertical swimming speed W(nn, np) to be specified, and turbulent mixing. All plankton cells with the same values of k, nn, and np are said to be in the same ‘box’, see Fig. 2, and are assumed to be uniformly distributed in the box. In a time interval dt each cell undergoes a ‘displacement relative to the width of the box’ Dn, Dp, Dk, in the directions of increasing nn, np, k of: Dn = (dEN/dt)*dt/LN = (N(z*) −EN(nn*))*dt/LN(nn),

Dk = −W(nn, np) dt/dz.

Thus, each box will, in time interval dt, lose a concentration C*Dn,p,k in the direction of increasing (decreasing), nn, np, k if those displacements are positive (negative). Each box will also gain concentration from adjacent boxes. The net change for each ‘advective’ process, for each box, DC(I), where I stands for nn, or np or k, with the other two integer variables fixed and D for Dn, Dp, or Dk in Eqs. (7a), (7b) and (7c) below, is as follows: DC(I)= − S(I)*D(I)*C(I) +0.5*(S(I−1)+1)*(D(I−1)*C(I−1)) +0.5*(S(I+1) −1) *D(I+1)*C(I+1) 1BIBN,

(6a)

(7a)

DC(1)= −0.5*(S(1)+1)*D(1)*C(1) +0.5*(S(2)−1)*D(2)*C(2),

Dp = (dEP/dt)*dt/LP =(SC*P(z*,t) − EP(np*))*dt/LP(np),

Fig. 2. Geometry of a box. Specific example for k=3, np =3, nn = 9. The arrows indicate the direction of migration of cells into adjacent boxes (not shown).

(7b)

DC(N)= − 0.5*(S(N) −1)*D(N)*C(N) +0.5*(S(N −1) +1)*D(N − 1)

(6b) (6c)

We note that the asterisk variables indicate quantities evaluated at the center of each box and that the negative sign in Eq. (6c) follows since W is positive upwards and k increases downwards.

*C(N−1)

(7c)

In Eqs. (7a), (7b) and (7c), S is the sign of the displacement D, either + 1 or −1 and N is 10 for nn and np, and 50 for k. In our algorithm the displacements, Eqs. (6a), (6b) and (6c) are first calculated for all nn, np, k. Eqs. (7a), (7b) and (7c)

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

are then utilized three times. For each value of nn, np, these equations are used for k = 1 through 50. Then for each value of np, k, we use Eqs. (7a), (7b) and (7c) for nn=1 – 10. Finally for each value of nn, k, we use Eqs. (7a), (7b) and (7c) for np = 1–10. We note this procedure is equivalent to ‘upwind’ differencing which is numerically stable. The total advective rate of change for each box is then computed by summing the three advective changes to compute a net advective change rate, Cadv(k, nn, np) = (DC(k) + DC(nn)+ DC(np))/dt, (8) Finally at each time step, each subgroup (nn, np) is turbulently mixed in the vertical to a prescribed depth, Dmax with a constant eddy diffusivity, i.e. we solve dC(z, nn, np)/dt =Av d2 C(z, nn, np)/dz 2 + Cadv(z, nn, np),

(9)

with zero flux at z= 0 and z = − Dmax, via the Crank–Nicholson scheme. The updated, in time, value for the concentration for each k, nn and np box in the mixed region is obtained from Eq. (9). Below this mixed region, if a mixed layer is taken to exist, the updated concentration is given by

241

to PAR is taken to occur even in nighttime as a proxy for the decrease in stored carbon due to respiration and biosynthesis as determined by the temporal decay since the last exposure to PAR. We use a continuous function of nn and np which reflects this behavior of the form: W(nn, np) =2.82 ×10 − 4m s − 1*( −exp(−(nn/5)2) +(1.0+tan h(0.7(nn −8))) *(tan h(0.7(8−np)))

(11)

The values of W/2.82 ×10 − 4 m s − 1 are given in Fig. 3. The contours were drawn with Golden Software’s Surfer using triangulation with linear interpolation. The maximum swimming speed is − 2.71×10 − 4 m s − 1 or −0.97 m h − 1 for nn = 1. We see here the necessity of the expanded Eulerian model for this nonlinear swimming behavior. Consider a case when a layer contains only equal cell densities of two different classes, say, one class with np=4 and nn=6 and a second class with np=8 and nn=8. Both of these classes would swim downwards according to Eq. (11). However the mean values np and nn for the total population are np=6 and nn= 7 with a

C(k, nn, np, t+ dt) = C(k, nn, np, t) + dt*Cadv(k, nn, np).

(10)

To complete the model specification, the vertical swimming speed must be determined as a function of the exposures. We specify the following hypothetical swimming behavior. A cell will first seek an adequate level of nitrate exposure to fill its nitrogen requirement and, when this level is achieved, swim to facilitate beneficial PAR exposure to fill its carbon requirement. For levels of EN B 0.6, the cell will swim in the direction of increasing N (downwards). When this level of nitrogen exposure is achieved then the cell responds to EP. For EP ] 0.75, the cell will swim downwards to decrease its PAR exposure because of phototinhibition and for EN 5 0.625 will swim upwards to increase PAR exposure; this response

Fig. 3. Swimming velocity/1.02 m h − 1 as a function of nn and np.

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

242

positive value of W. Thus, using mean properties would incorrectly specify swimming behavior. This example may appear extreme, but it does occur in our calculation. We next turn to two sample calculations, one without turbulent mixing and a second, for contrast, with mixing.

3. Results and discussion In this section we discuss the results of two 1-day experiments performed with this model to demonstrate the nature of the PAR and nitrate fields and the types of outputs that can be readily obtained. In the first experiment mixing was ‘turned off’, Av = 0, and in the second a reasonable value of Av was specified. This allows us to examine some impacts of mixing. The parameters utilized in our numerical experiments are as follows: dt =200 s, dz = 1 m, Av/0, or Av = .001m2 s − 1, for the mixed case, ZN=ZP =Dmax = 10 m, and TP = TN =1 day = 86 420 s. The chosen finite value of Av would correspond to a wind speed of perhaps 5 m s − 1 and TP and TN are chosen so that exposures can change significantly in our 1-day runs. Our initial distribution of cells was obtained as follows. First, we distribute 104 cells m − 2 uniformly over the 50 m depth of our domain yielding an initial concentration of 200 cells m − 3 in each layer. The initial distribution is then held in place without swimming or mixing subject to the nitrate and PAR fields specified in equations (3) and (4) until they reach equilibrium with these fields, i.e. EN(z)= N(z) and EP(z)= exp(z/ZP); these values can be obtained from Fig. 1 with z evaluated at the mid-point of each layer. The calculation then starts and runs for 1 day, first through 12 h of darkness and then through 12 h of light with distributions saved every 2 h. We first discuss the total concentration of cells in each layer (the sum over nn and np from 1 to 10, with k fixed) as a function of depth and time

Fig. 4. Total cell concentration (no. m − 3) for Av = 0, the unmixed case, as a function of depth (m) and time (h).

for both the non-mixed (Av =0) and turbulently mixed cases are shown in Figs. 4 and 5, respectively. The contours were drawn as in Fig. 3. In the unmixed case, Fig. 4, the cell concentration reaches a maximum in layer 11 of 817 cells m − 3 at hour 8 with strong gradients above and below this layer. This maximum gradually dissipates through the remainder of the day. This behavior is relatively straightforward. Cells initially in the upper 10 m have low N exposures (leading to low stored nitrogen reserves) and swim downwards to increase this exposure. Cells below layer 11 have high N exposures but low PAR exposures (leading to low stored carbon reserves) and swim upwards to anticipated PAR exposure resulting in an accumulation in the nutricline region. After several more hours the previously low N population which has swum downwards towards layer 11 achieves sufficient N exposure and starts an upward migration to anticipated PAR exposure to facilitate carbon fixation joining those cells already swimming upwards. After 12 h most of the

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

population has high N exposures and relatively low stored carbon reserves due to decay of PAR exposure during 12 h of darkness and thus largely respond to PAR-based carbon needs and swim upwards resulting in an accumulation in the top layer. We now discuss the mixed case, Fig. 5. Mixing is confined to the top 10 layers and transport from layer 10 to 11 is strictly due to swimming behavior. A maximum concentration in layer 11 still occurs although it is weaker with a maximum concentration of 574 cells m − 3 at hour 4. Mixing decreases the overall downwards flux of cells from layer 10 to 11 by eroding the aggregation of cells,that would otherwise have occurred in layers 7 – 10 due to downwards swimming, back up into the upper part of the mixed layer. With fewer cells present in layer 10 in the mixed case, the swimming flux to layer 11 is diminished. As the differences between the two cases are greatest during the first 12 h, we will examine the results of hour 6 in more detail. In Figs. 6 and 7, we present the

243

average values and standard deviations for nn and np for each layer. The values for these parameters in the mixed case are, predictably, more uniform than in the unmixed case in the upper 10 m where mixing occurs. Low EN and high EP cells initially present near the surface are mixed downwards and high EN low EP cells are moved upwards to achieve this uniformity. The effects of mixing are felt in the nn –np distribution of cells in layer 11 at hour 6 which is given in Table 1. Although there are 50% more cells present in this layer in the unmixed case (739) than in the mixed case (505), there are less than half as many low N cells (nn=1) which originated close to the surface. Mixing has brought these cells downward to layer 10 from whence they have swum into layer 11. These results can also demonstrate the failings in using mean values of the exposures versus our more exact calculation for our prescribed swimming behavior. For the unmixed case, the average values of np and nn in layer 11 at hour 6 are 4.88 and 6.69, respectively. Using these values in Eq. (11), we find the swimming speed is − 0.18 m h − 1. However the weighted swimming speed using the distribution given in Table 1, i.e. the sum over the variables nn and np from 1 to 10 of C(11, nn.np, 6 h))×W(nn,np)/Ctotal(11, 6 h) is+0.11 m h − 1. Thus, both the magnitude and sense of the swimming behavior is mis-specified using the mean state. Again, this is a consequence of using a non-linear physiological response. When the physiology is non-linearly forced by the environment, the expanded Eulerian approach should prove useful.

4. Conclusions

Fig. 5. Total cell concentration (no. m − 3) for Av =.001 m2 s − 1, the mixed case, as a function of depth (m) and time (h).

We have described the specification of an expanded Eulerian model of plankton physiological response (swimming behavior) to environmental exposures. The exposures to nitrate and PAR, each with one time scale, were utilized as independent variables, along with depth and time to support calculation of cell concentration. The model can be readily extended to include more environmental fields and/or time scales. This type

244

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

Fig. 6. The distribution of exposures with depth at hour 6 for Av = 0, the unmixed case, means and standard deviations.

of model, or a more traditional Lagrangian model which tracks large numbers of individual cells, is necessary when the physiological response under

consideration for a population does not respond linearly to population averages of the exposure variables.

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

The expanded Eulerian model presented here falls between the endpoints in the dichotomy provided by the strict Lagrangian and Eulerian approaches as presented in Lande and Lewis (1989).

245

The ease of the Eulerian calculation is maintained while information on the variety of cells contributing to the population occurring at a given depth at a given time is provided.

Fig. 7. The distribution of exposures with depth at hour 6 for Av =0.001 m2 s − 1, the mixed case, means and standard deviations.

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

246

Table 1 Cell concentrations (no. m−3) in layer 11 at hour 6 nn = 2

nn= 3

nn =4

nn = 5

nn =6

nn =7

nn = 8

nn =9

0 4 7 5 2 0 0 0 0 0

0 3 11 14 8 2 0 0 0 0

0 0 6 14 13 5 1 0 0 0

0 0 1 13 16 7 1 0 0 0

0 0 0 4 18 15 5 0 0 0

0 0 0 0 43 43 15 2 0 0

0 0 0 0 7 91 62 12 1 0

0 0 0 0 0 50 43 10 1 0

0 0 0 0 0 14 52 19 2 0

0 0 0 0 0 0 65 28 4 0

m2 s−1 1 8 18 24 21 10 2 0 0 0

0 3 6 9 12 8 2 0 0 0

0 0 1 2 5 4 1 0 0 0

0 0 0 0 5 5 2 0 0 0

0 0 0 0 2 2 0 0 0 0

0 0 0 0 0 4 3 0 0 0

0 0 0 0 0 74 52 10 1 0

0 0 0 0 0 4 5 1 0 0

0 0 0 0 0 9 24 8 1 0

0 0 0 0 0 13 100 41 6 0

nn= 1 Av =0 np= 10 np= 9 np= 8 np= 7 np= 6 np=5 np=4 np=3 np=2 np=1 Av =10−3 np =10 np =9 np =8 np= 7 np= 6 np= 5 np= 4 np= 3 np= 2 np= 1

Acknowledgements

Ex *i

This work was supported by NASA NAGW 3575-2 (D.K.and G.S.J.), NSF OCE 95-03253 (D. K.) and ECOHAB:Florida Program subcontract from Mote Marine Laboratory MML 106685 (D.K.).

EXimin, EXimax M

dz Appendix A. Nomenclature k Xi

EXi (t,T)

T

An external field of physiological importance. The index i ranges one to the final number of external fields used, here 2. The exposure of a cell to field Xi at time t over the time scale T. Defined in Eq. (1). An exposure time scale, which may vary from field to field or may have multiple values for the same field defining different EXi.

N(z) NIT, NITmax ZN P(z,t) PAR, PARmax

nn = 10

The exposure variable scaled to between 0 and 1. The asterisk is dropped after Eq. (3). The minimum and maximum values for unscaled EXi. The number of contiguous horizontal layers into which the water is subdivided. The thickness of each flouid layer, here 1 m. An index for the depth indicating fluid in the layer (k−1) dz to kdz below the surface. The scaled nitrogen field here taken to be independent of time. The unscaled nitrate field and its maximum value. A length scale used in defining N(z). Here taken as 10 m. The scaled PAR field The unscaled PAR field and its maximum value.

G.S. Janowitz, D. Kamykowski / Ecological Modelling 118 (1999) 237–247

ZP

The e-folding distance for PAR. Here taken as 10 m. EP, EN The scaled exposures to PAR and nitrate, replaces EX *1 and EX *. 2 SC The scaling for EP and equals 1/EPmax. nn An integer, from 1 to 10, denoting the range of EN from (nn−1)/10 to nn/10. np An integer, from 1 to 10, denoting a range of EP. LN(nn) The length of the nnth interval. Here taken as .1. LP(np) The length of the npth interval. Here taken as non-uniform and defined in the text. C (k, nn, np, t) The number of cells per m3 in layer k at time t with exposure ranges defined by nn and np. Dmax The maximum depth of turbulent mixing. Here taken as 10 m. Dn, p, k Displacement of cells in a box in the nn, np, or k directions in time dt relative to the width of the box in that direction. DC(I) The net change in concentration in a box due to displacements in the I direction in time dt, where I is nn, or np, or k. Cadv The net rate of change in a box due to displacements in all three ‘directions’.

W (nn, np)

247

The swimming speed (positive upwards) for a cell with exposures nn and np.

References Cullen, J.J., Lewis, M.R., 1988. The kinetics of algal photoadaptation in the context of vertical mixing. J. Plankton Res. 10, 1039–1063. Franks, P.J.S., Marra, J., 1994. A simpled new formulation for phytoplankton photoresponse and an application in a wind-driven mixed-layer model. Mar. Ecol. Prog. Ser. 111, 143–153. Janowitz, G.S., Kamykowski, D., 1991. An Eulerian model of phytoplankton photosynthetic response in the upper mixed layer. J. Plankton Res. 13, 983–1002. Kamykowski, D., 1995. Trajectories of autotrophic marine dinoflagellates. J. Phycol. 31, 200–208. Kamykowski, D., Yamazaki, H., Janowitz, G.S., 1994. A Lagrangian model of phytoplankton photosynthetic response in the upper mixed layer. J. Plankton Res. 16, 1059–1069. Lande, R., Lewis, M.R., 1989. Models of photoadaptation and photosynthesis by algal cells in a turbulent mixed layer. Deep-Sea Res. 36, 1161–1175. McGillicuddy, D.J. Jr, 1995. One-dimensional numerical simulation of primary production: Lagrangian and Eulerian formulations. J. Plankton Res. 17, 405–412. Muller, F., 1997. State-of-the-art in ecosystem theory. Ecol. Model. 100, 135–161. Wolf, K.U., Woods, J.D. 1988. Lagrangian simulation of primary production in the physical environment-the deep chlorophyll maximum and the nutricline. In: Rothschild, B.J. (Ed.) Toward a Theory on Biological–Physical Interaction in the World Ocean. Kluwer, Dordrecht, pp. 51–70. Woods, J.D., Onken, R., 1982. Diurnal variation and primary production in the ocean-preliminary results of a lagrangian ensemble model. J. Plankton Res. 4, 735–756. Woods, J., Barkmann, W., 1993. Diatom demography in winter-simulated by the Lagrangian ensemble method. Fish. Oceanogr. 2, 202–222. Yamazaki, H., Kamykowski, D., 1991. The vertical trajectories of motile phytoplankton in a wind- mixed water column. Deep-Sea Res. 38, 219–241.

.