## Abstract

A nonlinear model for biological and physical dynamical interactions in a laminar upwelling flow field in parts I and II of this study is extended to turbulent flow. In the previous studies, a prescription for obtaining quadrature solutions to the fundamental biodynamical equations was developed. In this study, we use a probability density function approach on these solutions to obtain statistics of the biodynamical state variables and their self-interaction for the case of turbulent advection. To illustrate the theory, a simple nutrient (*N*), phytoplankton (*P*) problem is considered, that of upwelling into a surface turbulent layer. Biological interaction is modelled as bilinear, representing the uptake of *N* by *P* in a uniform light euphotic zone. A random walk model is used to obtain the appropriate probability density function for the advective turbulent field. The mean quantities, , , as well as the biological interaction term are calculated. The term has two contributions, , and the turbulence-induced interaction term, . It is shown that the often neglected turbulence-induced coupling term is of the order and opposite in sign. This results in, over a wide range of Peclet numbers, the mean interaction term being significantly smaller than either of its constituent terms, and .

## 1. Introduction

There is great interest in understanding the role that turbulence plays in affecting the dynamics of upper ocean plankton ecosystems. Recent research has focused on the relative dominance of physical forcing and biological processes in controlling coupled biological–physical dynamics. Yamazaki *et al*. (2002) provide a general overview and cite new evidence for small-scale turbulence influencing the population dynamics of plankton. Random small-scale turbulence and associated small-scale coherent structures can affect primary production, feeding, predation, aggregation and mating. Spectra of turbulence and patchiness of plankton have been reviewed by Okubo & Mitchell (2001). The spectra of the plankton follow that of the turbulence for shorter time scales (1–10^{2} min), but not for longer scales where growth and grazing become important. Denman *et al*. (1977) derived theoretical spectra for both scale ranges. Cullen *et al*. (2002) developed the classification structure of pelagic ecosystems in terms of high/low turbulence and high/low nutrients, a concept first introduced by Margalef (1978).

Modelling coupled physical/biogeochemical-ecosystem dynamics in the sea is complex and highly nonlinear. Nonlinearities arise not only from physical advection of the biological state variables, but also from a number of interactions among the biological state variables representing nutrient uptake by phytoplankton, grazing, predation, etc. Fundamental modelling ambiguities and issues occur because: (i) there is no counterpart for biological dynamics to the Navier–Stokes equations for physical fluid dynamics, (ii) the continuum limit for concentration distributions of plankton involves finite-sized individuals, and (iii) the number of biological state variables in the real ocean is very large (e.g. numerous species with age classes, many macro- and micro-nutrients, etc.). Nonetheless, coupled dynamical modelling research is being vigorously pursued in terms of: four-dimensional Eulerian field equations; individual- (or cohort) based Lagrangian models; and hybrid models. (See reviews by Hofmann & Friedrichs 2002; Robinson & Lermusiaux 2002; Yamazaki *et al*. 2002; Runge *et al*. 2005.)

Introduction of turbulence in coupled ecosystem models further complicates an already very complex problem. Owing to the nonlinearities noted above, turbulence affects plankton biodynamics in two general ways. These are by: (i) advective stirring and mixing and (ii) inducing zero average fluctuations about the mean biological fields in the nonlinear biodynamical interactions. Products of such fluctuations generally correlate and are non-zero in the mean. The first effect is typically modelled by a type of eddy diffusivity for the diffusion of the biological state variables, such as *K* profile parameterization (Large *et al*. 1994) or Mellor/Yamada (Mellor & Yamada 1982). The second effect represents some very complex interactions induced on the biodynamical coupling of the ecosystem. To model this latter effect, an additional assumption beyond that of the eddy diffusivity assumption for the turbulent field has to be made in order to close the coupled physical–biological mean field equations. However, this is rarely done formally and explicitly.

Commonly, only mean-field/mean-field interactions are taken into account and fluctuation correlation effects ignored. Donaghay & Osborn (1997) do recognize this problem and recommend treating these difficult terms (in their case, mean nonlinear birth and mortality rates of plankters) as input information or forcings. However, they proceed only qualitatively thereafter. Okubo *et al*. (2001) adopt a deterministic viewpoint to study the nonlinear instability of a diffusive predator–prey system in terms of finite amplitude perturbations of an equilibrium state.

Numerical simulations of the effect of turbulence in biodynamical coupling problems is typically approached by using some type of random walk model for the turbulent fluid element displacement (Visser 1997; Yamazaki *et al*. 2002). Care must be taken when using a spatially varying diffusivity for the random walk (Hunter *et al*. 1993; Ross & Sharples 2004) to insure that a spurious build-up of particles does not occur.

Alternative analytical approaches to solving the general advective–diffusive problem have recently been explored by Pope (1994) and Wilson & Sawford (1996). These approaches involve obtaining the probability density function, which is then applied to the field variables for their statistics. However, development of analytical solutions for advective–diffusive–reactive problems such as those modelling the effect of turbulence on primary production depends upon the ability to be able to formulate this problem both in Lagrangian as well as Eulerian coordinates. The diffusive aspect of the problem is essentially Lagrangian, while the measurement aspect and much of the theoretical underpinning are essentially Eulerian.

In this study, we develop an analytical approach to obtaining the statistics of the biodynamical variables and their coupling by extending a formalism developed in parts I and II of this series, Robinson (1997, 1999); hereafter known as R97, R99. Robinson uses the method of characteristics to obtain a set of symbolic solutions to the advective–reactive equations associated with a general reaction of *n*-tuple biodynamical interactions in the presence of general nonlinear advective flow. Then, he details a number of particular solutions relevant to the *NP* and *NPZ* (nutrient–phytoplankton–zooplankton) problems for a variety of flow field cases including coastal and open-ocean upwelling into an upper ocean euphotic zone. Here, we use the formalism developed in R97, R99 for the case of purely laminar-advective upwelling as a starting point and extend the theory to include the superposition of turbulence on the advecting flow field.

A key factor in allowing the extension of the Robinson theory to include turbulence arises from the fact that the fundamental set of equations imposes minimal restrictions on the form of the advection velocity. A constraint on the form of the advection velocity arises when the solutions for the biological state variables need to be expressed in terms of Eulerian coordinates. This restriction is handled for the turbulence case by using a probability density function approach to obtain ensemble averages of the biodynamical variables and their coupled interaction.

Section 2 of this study develops the theory. Section 3 applies that theory to a simple example problem, *NP* upwelling into a uniform light turbulent euphotic zone. Section 4 discusses the results obtained for the ensemble averaged mean quantities, , , and the biodynamical coupling with detailed emphasis on the significance of the turbulence-induced coupling term . A summary and conclusions are presented in §5.

## 2. Fundamental theory

### (a) Robinson model

A nonlinear *NPZ* model for laminar upwelling into a euphotic zone has been developed in the series of studies by R97 and R99. The fundamental equations used in the model are(2.1)where repeated indices are summed. The continuity equation for the flow field has been invoked. The variables *ϕ*_{i} are the biological state variables and *F*_{i} represents the interactions between these biological state variables.

Note thatcan be expressed as a function of both ** x** and

*t*in Eulerian coordinates and as a function of in Lagrangian coordinates. Here,(2.2)describes the Lagrangian trajectory of a fluid element at position

**at time**

*x**t*, which was at position

**at time**

*X**t*=

*t*

_{0}. Robinson employs the method of characteristics for a variety of laminar advective flow fields to obtain solutions. The subscript L on a spatial variable indicates a Lagrangian functional dependence. An upper case letter is used to indicate the initial position of a fluid element trajectory, i.e. . For the case of no interaction,

*F*

_{i}=0, the

*ϕ*

_{i}are not coupled and (2.1) represents the advection of a passive scalar.

It is straightforward to see that the formal solution to (2.1) can be written in terms of Lagrangian coordinates, whence(2.3)where are the values of the biological state variables at *t*=*t*_{0} at position ** X**. The Eulerian solution to (2.1)(2.4)is obtained by using the Lagrangian trajectory (2.2) in (2.3).

### (b) Formalism for adding turbulence

Now, we show the prescription for formally extending the Robinson model to include a turbulent field superposed on a laminar advective field. Equations (2.3) and (2.4) remain solutions to (2.1) for both the turbulent and laminar cases. For the turbulence case, the problem with using the fluid trajectory given by (2.2) to obtain the Eulerian solution (2.4), is twofold. First, equation (2.2) cannot, in general, be expressed in closed form. Second, for the turbulence case, the trajectory is a random variable.

Our approach is to obtain ensemble averages of the solution of the biodynamical variables (equation (2.4)) using Eulerian statistics. In the Eulerian frame, defined at a fixed point ** x** and time

*t*, we need to obtain the statistics of the initial position

**in terms of the final fixed point position,**

*X***, time,**

*x**t*, and initial time,

*t*

_{0}. That is, we need to obtain the probability density function (PDF) for(2.5)which is just the inverse of (2.2). We introduce the subscript E to indicate an Eulerian or, fixed point, coordinate. Thus, for the turbulent case, we will use the general solution (2.4) along with the conditional PDF(2.6)to obtain ensemble averaged statistics of

*ϕ*

_{i}=

*ϕ*

_{i}(

**,**

*x**t*). The variables to the right of the slash in (2.6) are the conditioned variables.

The mean of the state variables are given by(2.7)and the cross variance by(2.8)where is the fluctuating component. The second term on the r.h.s. of (2.8) is given by(2.9)We will use both the overbar as well as bracket notation to indicate an ensemble average.

To obtain *p*_{E}, some assumption needs to be made about the statistics of the turbulent field itself. This is accomplished here by assuming that at each point a fluid element undergoes a random increment or ‘walk’. This results in determining the Lagrangian PDF(2.10)Bayes' theorem (Papoulis 1965) is used to relate *p*_{E} to *p*_{L}, namely(2.11)Equation (2.11) is also seen as the ratio of the fractional number of ensembles of paths that arrive at ** x** at time

*t*from the initial position

**divided by the total number of such paths.**

*X*Finally, it is interesting to contrast this approach of explicitly using the PDF (2.11) on the solution (2.4) with the standard approach of taking the ensemble average of (2.1), which results in(2.12)The eddy diffusivity approximation (Okubo & Mitchell 2001) is typically invoked, whence(2.13)where *κ*^{*} is the diffusivity, with the asterisk indicating that it is in dimensional form. If , then the usage of (2.13) in (2.12) results in the advective–diffusion equation for .

However, in order to proceed in the case of , some additional assumption has to be made. Often this consists of assuming that the interaction of the total state variable can be approximated by the interaction of only its ensemble mean (Donaghay & Osborn 1997), i.e.(2.14)Usage of (2.13) and (2.14) in (2.12) results in the advective–diffusive reaction equation for . However (2.14) is, in general, not valid and involves neglecting terms of the form of (2.9), i.e. setting(2.15)

In §3, we illustrate the application of the PDF approach to a very simple example problem considered by Robinson, namely that of a two-component *ϕ*_{1}=*P*, *ϕ*_{2}=*N* nutrient, phytoplankton (*NP*) bilinear interaction in a linear strain upwelling field. We will explicitly evaluate the mean interaction term given by (2.14), which for this example becomes(2.16)We will show that, not only is but that, in fact, it is comparable in magnitude over a wide range of turbulence values to . Thus, the turbulence interaction term plays a significant role in the biodynamics of *P*,*N* interaction and cannot, in general, be neglected.

## 3. Example problem: bilinear *NP* interaction in linear strain upwelling

In order to focus on turbulence-induced biological fluctuation interactions, we will illustrate and apply the theoretical approach developed in §2 to a very simple *NP* problem considered in R99. Figure 1 provides an illustration of the problem. Seed phytoplankton and nutrients contained in an infinitely deep and wide reservoir are upwelled into a euphotic zone of depth *Z*_{e}. There, the phytoplankton and nutrients undergo a simple bilinear interaction. No zooplankton grazing occurs and there is no mortality. The light field in the euphotic zone is taken to be uniform in the vertical, which in R99 is shown to be qualitatively similar to attenuated light.

For this model, using (2.1) with *ϕ*_{1}=*P* and *ϕ*_{2}=*N*, it follows that(3.1)and(3.2)where *z* is the vertical coordinate, taken as positive indicating increasing depth; *Z*_{e} is the depth of the euphotic zone. The quantity *a* is the phytoplankton growth rate coefficient, which has units of volume/(time)(mass). Adding (3.1) and (3.2) yields the conservation of biomass equation(3.3)where the total biomass density is given by *M*=*P*+*N*. Note that (3.3) has the simple solution(3.4)where the subscript ‘0’ indicates its initial value at time *t*=*t*_{0}. Following R99, we normalize *P* and *N* by *M*_{0}, length scales by *Z*_{e}, time scales by *M*_{0}/*a*, the nutrient uptake rate. Solutions for equations (3.1) and (3.2) with this normalization are readily seen as(3.5)and(3.6)where *t*_{e} is the time that a fluid element is located at *z*=1, the base of the euphotic zone.

### (a) Pure advection

Following R99, to obtain the spatial dependence of *P* and *N*, we need to assume some type of flow field. Consider the purely advective linear strain case of R99(3.7)and(3.8)where *α*=*α*^{*}*M*_{0}/*a*, with *α*^{*} as the dimensional strain rate. Equations (3.7) and (3.8) have the Lagrangian trajectories(3.9)and(3.10)where *X* and *Z* are the initial positions. We will use the subscript ‘adv’ to indicate a variable associated with the laminar advective flow. Substitution of (3.10) into (3.5) and (3.6) results in(3.11)and(3.12)Solution (3.11) with *z*_{adv}=*z* is equation (4.3*b*) of R99. Note that with and , equations (3.11) and (3.12) have the Lagrangian functional dependence *P*=*P*(*Z*, *t*, *t*_{0}) and *N*=*N*(*Z*, *t*, *t*_{0}).

### (b) Combining advection and turbulence

Now, we will superpose within the euphotic zone a turbulent field on the advective velocity field (figure 1). Below the euphotic zone the flow is laminar and is given by (3.7)–(3.10). Equations (3.11) and (3.12) remain as solutions for *P* and *N* but, located at the Lagrangian vertical position,(3.13)whence, in (3.11) and (3.12) . Figure 2 shows for the turbulence case two possible Lagrangian trajectories, A_{1}B and A_{2}B. The dotted lines are the mean Lagrangian paths *z*_{adv}.

Following the theory shown in §2, we need to obtain the Eulerian PDF, *p*_{E}. Then, we will apply this to solutions (3.11) and (3.12) to calculate ensemble averaged statistics of *P* and *N*. As discussed in §2, to obtain *p*_{E} we will first obtain the Lagrangian PDF, *p*_{L}, and then use Bayes' theorem, equation (2.1), to relate *p*_{E} to *p*_{L}.

Owing to the horizontal symmetry of our problem, solutions (3.11) and (3.12) have spatial dependence only in the vertical. The vertical velocity can be decomposed into a mean and turbulent fluctuation component, namely(3.14)where(3.15)

Differentiating (3.13) with respect to time *t* and using (3.10) yields(3.16)However, using (3.14) and (3.15) gives(3.17)Comparing (3.17) and (3.16) results in(3.18)which has the solution(3.19)At this point, no assumption on the nature of the turbulent field has been made. To model the vertical turbulent velocity fluctuation *w*′, consider a one-dimensional random walk model(3.20)where *κ* is the normalized diffusivity, taken to be constant. Using the random walk model (3.20) for a vertically unbounded geometry, the solution for the perturbed random displacement *z*′ from (3.19) is seen as a sum of uncorrelated weighted random increments . Thus, by the central limit theorem (Papoulis 1965) has the Gaussian distribution(3.21)where the mean displacement is given by (3.10) and the variance is given by(3.22a)and(3.22b)Subsequent in the study, we take *t*_{0}=0. The Peclet number *Pe*=*α*/*κ*=*α*^{*}/*k*^{*} in (3.22*a*) and (3.22*b*) represents the ratio of the turbulent time scale to the advective time scale.

The difference in the form of the variances shown in (3.22*a*) and (3.22*b*) arise from the difference in the amount of time in which a fluid element has spent in the turbulent euphotic zone. In figure 2, this is illustrated by showing two paths that arrive at the same point but originate either below or above the euphotic zone, respectively, and will be each associated with variances (3.22*a*) and (3.22*b*), respectively.

For our case we have a turbulent layer of finite domain between *z*=0 and 1, which requires replacing (3.21) with the form appropriate for a bounded domain by applying the no turbulent flux boundary conditions (Papoulis 1965)(3.23)This boundary condition insures that, for each ensemble, the turbulent displacement perfectly reflects from both the top and bottom boundaries, beyond which there is no turbulent fluid.

Using the method of images on (3.21) to satisfy boundary conditions (3.23) results in(3.24)Note that, as required from the definition of a PDF, boundary conditions (3.23) result in the constraint being satisfied. It should also be noted that both (3.21) and (3.24) satisfy the advective–diffusion equation. The Eulerian PDF is obtained by the prescription given by (2.11), which for our one-dimensional case can be written as(3.25)

## 4. Results

The Eulerian PDF *p*_{E} (3.25), together with solutions (3.11) and (3.12), will now be used to obtain the vertical profiles of and , and the turbulence-induced interaction term , whence(4.1a)(4.1b)and(4.1c)withandNote that the mean of the biodynamical total interaction term *F* of (3.1) and (3.2) is given by(4.2)

In the series of figures below, we will show profiles of the biomass quantities given by (4.1*a*)–(4.1*c*) and (4.2). Owing to their dependence on turbulence, all of these quantities are functions of the Peclet number. In addition, they are also functions of the same input parameters that would arise in the no turbulence pure advective upwelling case considered in R99. These input parameters are *P*_{0}, *N*_{0}, the normalized seed nutrient and phytoplankton density, taken to be constant in the reservoir where *Z*≥1, and *α* is the normalized strain rate.

For the figures given below, we will choose one set of the input parameters used extensively in examples presented in R99, namely *P*_{0}=0.01 and *N*_{0}=0.99, which represent an initial seed density of phytoplankton of 1% of the total nutrient and phytoplankton biomass density. We will also use *α*=0.1, which, for a nutrient uptake time of 12 hours, corresponds to a linear strain rate of *α*^{*}=2×10^{−6} s^{−1}, a reasonable value of open-ocean upwelling (R99).

To determine values of the Peclet numbers to use, we note that the typical turbulent vertical eddy diffusivity *κ*^{*} in an ocean surface mixed layer ranges from *κ*^{*}=10^{−4} to 10^{−2} m^{2} s^{−1} (Large *et al*. 1994). Using a mixed-layer depth of 30 m and an upwelling strain rate of *α*^{*}=2×10^{−6} s^{−1}, we obtain expected values of *Pe*≲10. Thus, we will use the three values *Pe*=0, 1 and 10, as well as a fourth value, *Pe*=∞; the latter corresponds to the no turbulence pure advection case considered in R99. Note that the *Pe*=0 case corresponds to infinitely large turbulence relative to the advection.

### (a) Time-dependent solutions

In figure 3*a*,*b* we show the solutions of and for *Pe*=10, a moderate turbulence case, for different times *t*, taking the initial time *t*_{0}=0. Shown in these figures as a solid line is the steady-state limit *t*→∞. At the bottom boundary, and initially increase with time *t*. Note the discontinuity across *z*=1. This results from our model of upwelling into a turbulent mixing region confined to a finite depth range of *z*=0 to 1. The no turbulent flux boundary conditions (3.23) insures that both and have zero vertical gradients at *z*=0 and 1, with the result that no mean turbulent flux of biomass material occurs across these boundaries.

From the form of the solutions (3.11) and (3.12) using (3.10), it is known that, for the laminar case, *P* and *N* profiles asymptote to a steady state over a time period , where, for our example problem input parameters, . For the profile in figure 3*b* this appears to be the case, but it is clearly not the case for the profile of figure 3*a*, which asymptotes to a steady state at the much later time, *t*≈50, not shown in figure 3*a*. To understand this, note that after a time *t*≫*τ*, the nutrient composition within a fluid element becomes depleted and the fluid element contains only phytoplankton. In a turbulent field, the fluid elements are rearranged to depths different from that which would occur under pure advection. In the steady state, contributions to the nutrient profile principally come from fluid elements that have resided in the mixed layer for the residence time . The phytoplankton density has continued to increase beyond *t*=*τ* due to the continual input of the new fluid elements carrying nutrients into the turbulent layer. These nutrients are taken up by phytoplankton as turbulence mixes these new fluid elements throughout the region *z*≤1, with the degree of mixing dependent upon the Peclet number.

To illustrate further the nature of this process, consider the intense turbulence limiting case *Pe*→0, which, with the change of variable in (3.25), leads to(4.3)Using (4.3) in (4.1*a*)–(4.1*c*) with solutions (3.11) and (3.12) results in(4.4a)(4.4b)and(4.4c)where, for this limiting case of infinitely intense turbulence, the mean values shown in (4.4*a*)–(4.4*c*) are just the vertically averaged values of the purely advective case. This result, as expected, implies that the intense turbulence case of *Pe*=0 perfectly mixes *P*, *N* and the interaction term (*PN*).

In figure 4 for *Pe*=0, we show such plots of , , and versus time, whereand and asymptote to a steady state at time of order , a similar time when asymptotes to a steady state for the moderate turbulence case shown in figure 3*b*. From figure 4, we see that the quantities and asymptote to a steady state about an order of magnitude later than and , again a similar result to the moderate turbulence (*Pe*=10) case of figure 3*a*. Note the interesting results of figure 4 that for all time and for time , , while for , . This latter result follows from the fact that total biomass is conserved and that as *t*→∞ throughout the turbulent layer.

### (b) Steady-state solutions

We show in figure 5 the steady-state (*t*→∞) profiles of mean phytoplankton density as a function of the Peclet number *Pe*. Figure 5 exhibits the expected behaviour: the more intense turbulence (smaller *Pe*) mixes the phytoplankton more effectively than the less intense (higher *Pe*) turbulence does. This results in smaller gradients for the profiles with decreasing *Pe*. Figure 6*a*,*b* shows the mean interaction term and the two terms of which it is composed, and (equation (4.2)). Figure 6*a* shows the expected behaviour of as a function of *Pe*. Decreasing *Pe* decreases the gradient of . Note that, in figure 5, the depth of convergence of the solutions of is about *z*≈0.6. This depth is close to the depth where the mean interaction term, , shown in figure 6*a*, is a maximum.

In figure 6*b*, we show contributions of and to . For all of the values of *Pe*, except the non-turbulent value of *Pe*=∞, the magnitude of is of the order , and opposite in sign,(4.5)This result clearly demonstrates that neglecting the turbulent interaction term would result in an order of magnitude overestimate of .

To understand the overall effect of turbulence on net phytoplankton growth, we plot, in figure 7 as a function of *Pe*, the vertically integrated net phytoplankton density defined by(4.6)along with the vertically integrated local growth rate term , where(4.7)with *F* given by (4.1*a*)–(4.1*c*) and (4.2). In the extreme limits of *Pe*→0 and *Pe*→∞, as indicated in figure 7, there is a balance between these two terms(4.8)A maximum in is seen to occur at *Pe*≈5. To understand this behaviour of , note that the case of *Pe*=0 is characterized by uniform mixing and zero vertical gradients for and , with the result that(4.9)This can also be seen in figure 5, where the vertical integral of the light line, the *Pe*=∞ case, equals that of the vertical integral of the dark line, the *Pe*=0 case. For the intermediate values of *Pe*, vertical gradients of and are maintained and turbulent transport makes locally available more nutrients for uptake than would occur in either the non-turbulent case or the extreme turbulent case.

## 5. Summary and conclusions

A major research area of contemporary interdisciplinary ocean science research is understanding and quantifying biological–physical interactions. Such interactions occur on multiple space and time scales and are typically characterized by nonlinear feedback. Fundamental processes that arise from the nonlinear biodynamical interactions of biological state variables in the presence of turbulent flow, for the most part, have been ignored in mean field dynamical studies. To include such effects by closure in terms of mean fields and gradients presents a formidable problem. However, if the statistics of the turbulence are known or parameterized, turbulence-induced fluctuation correlations of biological variables can, in principle, be obtained by using requisite probability density functions. This approach can also deal directly with the more familiar problem of turbulent diffusion of biological state variables. In the case that the turbulence arises from physical dynamics without biological feedbacks, considerable knowledge exists for its statistical parametrization. In the case that biodynamical interactions cause, or significantly influence, turbulent or chaotic flow, much more research is required.

In this study, we have shown that the analytical formalism developed by R97, R99 for modelling the biodynamical interaction of seed plankton and nutrients upwelling into a euphotic zone can be extended to include a turbulent layer. This results from the fact that in R97, R99, the biodynamical interaction terms on the r.h.s. of (2.1) are isolated from the fluid advection terms, which appear on the l.h.s. of (2.1). The separation of these terms allows quadrature solutions for the biological state variables to be obtained in terms of Lagrangian coordinates. Using the method of characteristics, R97, R99 has shown that these solutions can also be expressed in terms of Eulerian coordinates. This latter step is critical since observations and theories on which they are based are typically made in the Eulerian reference frame.

We extend the Robinson approach by considering the total advection flow field to consist of the superposition of a laminar and turbulent component. As in R97 and R99 for the advective case, we obtain the Eulerian solution (2.4) to the fundamental equation (2.1). To obtain the solution for the mean biological state variables and their coupled interaction, we use the probability density approach given in equations (2.10) and (2.11). Note that the alternative approach of using the equations for the mean biological state variables (2.12) and then modelling the turbulent flux terms by an eddy diffusivity (2.13) does not correctly or fully model the effect of the turbulent field and requires some additional assumption similar to that of (2.14). In other words, because the mean biodynamical interaction term is, in general, nonlinear and does not necessarily have some analytical prescription that allows it to be expressed in terms of the mean biological state variables, i.e. .

For the *NP* example in this study, we use a random walk model given in (3.20) and Bayes' theorem (3.25) to obtain the Eulerian probability density function. Then, we can use the solution to the fundamental equations to explicitly obtain the mean quantities and and the turbulence-induced interaction term . These quantities are functions of the Peclet number, which represents the ratio of the turbulent time scale to the advective time scale. Small values of *Pe* occur when turbulence dominates, large values occur when advection dominates.

Figures 3 and 5 show the expected effect of turbulent mixing in smoothing property profiles. Recall that the problem has been formulated to insure horizontal independence of the vertical profiles. Turbulence plays a critical role in making more nutrients available for uptake throughout the water column than would occur for the non-turbulent case. The continual advection of new nutrients into the euphotic zone, coupled with the effect of turbulence on nutrient uptake by the phytoplankton, results in the phytoplankton profile reaching a steady state at a later time than that of the nutrient profile. This effect can be seen in figures 3 and 4.

The critical result of this study has been to show that the turbulence-induced interaction term that, for the case of our simple model, is cannot be neglected. In fact, for the example discussed here, as shown in figure 6, it is comparable and opposite in sign to the mean field interaction term . In formulating solvable advection–diffusion reaction equations from (2.1), it is the latter type of term that is often retained and the turbulence-induced interaction terms neglected (Donaghay & Osborn 1997).

To elucidate some of the basic principles of nonlinear biodynamical interactions and the importance of fluctuation correlations, we have considered here a very simple *NP* interaction model. This model allows a clear elucidation of the two roles played by turbulence, that of random property flux or diffusion and that of perturbing the growth through nonlinear interaction. This study serves as a simple first step in examining in more detail the turbulent advective effects on biophysical dynamics in the sea. Extensions of this study might include: (i) using more realistic attenuation of light, (ii) improved representation of uptake and competitive nonlinear interactions of biological state variables such as the Michaelis–Menton kinetics of nutrient uptake (∼*PN*(*K*+*N*)^{−1}), zooplankton grazing (*PZ*) and quadratic mortality (*Z*^{2}), (iii) more complicated mean advection, and (iv) more realistic turbulent statistics representing larger scale anisotropic turbulent eddies.

## Acknowledgments

It is a pleasure to thank Mr Wayne G. Leslie for his substantial contribution to the production of this manuscript. We thank the Office of Naval Research for partial support of this research under grants N00014-02-1-0989 to Harvard University and N00014-04-1-0254 to the University of Massachusetts Dartmouth.

## Footnotes

- Received May 25, 2007.
- Accepted November 21, 2007.

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

- Copyright © 2007 The Royal Society