On the theory of advective effects on biological dynamics in the sea. III. The role of turbulence in biological–physical interactions

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 .


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 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 (1997Robinson ( , 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,N ,P, and the biodynamical coupling hNPi with detailed emphasis on the significance of the turbulence-induced coupling term hN 0 P 0 i. A summary and conclusions are presented in §5.

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 where repeated indices are summed. The continuity equation for the flow field has been invoked. The variables f i are the biological state variables and F i represents the interactions between these biological state variables.
Note that can be expressed as a function of both x and t in Eulerian coordinates and as a function of X; t; t 0 in Lagrangian coordinates. Here, x Z x L ðX; t; t 0 Þ ð 2:2Þ describes the Lagrangian trajectory of a fluid element at position x at time t, which was at position X at time tZt 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. X i Z ðX; Y ; ZÞ. For the case of no interaction, F i Z0, the f 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 f 0i Z f 0i ðX; tÞj XZ0;tZt 0 are the values of the biological state variables at tZt 0 at position X. The Eulerian solution to (2.1) 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), f i Z f i ðx; t; f 0i Þ is twofold. First, equation (2.2) cannot, in general, be expressed in closed form. Second, for the turbulence case, the trajectory x Z x L ðX; t; t 0 Þ is a random variable.
Our approach is to obtain ensemble averages of the solution of the biodynamical variables f i Z f i ðx; tÞ (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 X in terms of the final fixed point position, x, time, t, and initial time, t 0 . That is, we need to obtain the probability density function (PDF) for 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 p E h p E ðX; t 0 =x; tÞ ð 2:6Þ to obtain ensemble averaged statistics of f i Zf 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 and the cross variance by where f 0 Z f i Kf i is the fluctuating component. The second term on the r.h.s. of (2.8) is given by 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 p L Z p L ðx; t=X; t 0 Þ: ð2:10Þ Bayes' theorem (Papoulis 1965) is used to relate p E to p L , namely p E ðX; t 0 =x; tÞ Z p L ðx; t=X; t 0 Þ Ð dX p L ðx; t=X; t 0 Þ : ð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 X divided by the total number of such paths. 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 The eddy diffusivity approximation (Okubo & Mitchell 2001 where k Ã is the diffusivity, with the asterisk indicating that it is in dimensional form. If F i ½f j ; x; tZ 0, then the usage of (2.13) in (2.12) results in the advectivediffusion equation for f i . However, in order to proceed in the case of hF i ½f j ; x; ti s0, 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) ð2:14Þ Usage of (2.13) and (2.14) in (2.12) results in the advective-diffusive reaction equation for f i . However (2.14) is, in general, not valid and involves neglecting terms of the form of (2.9), i.e. setting hf 0 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 f 1 ZP, f 2 ZN nutrient, phytoplankton (NP ) bilinear interaction in a linear strain upwelling field. We will explicitly evaluate the mean interaction term F i given by (2.14), which for this example becomes We will show that, not only is hP 0 N 0 i s0 but that, in fact, it is comparable in magnitude over a wide range of turbulence values to PN . Thus, the turbulence interaction term hP 0 N 0 i plays a significant role in the biodynamics of P,N interaction and cannot, in general, be neglected.

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 f 1 ZP and f 2 ZN, it follows that ð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 where the total biomass density is given by MZPCN. Note that (3.3) has the simple solution where the subscript '0' indicates its initial value at time tZt 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 where t e is the time that a fluid element is located at zZ1, 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 u Z ax ð3:7Þ and w ZKaz; ð3:8Þ where aZa Ã M 0 /a, with a Ã as the dimensional strain rate. Equations (3.7) and (3.8) have the Lagrangian trajectories x Z x adv Z X exp½aðtK t 0 Þ ð3:9Þ and 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 Solution (3.11) with z adv Zz is equation (4.3b) of R99. Note that with P 0 Z P 0 ðZ; t 0 Þ and N 0 Z N 0 ðZ; t 0 Þ, equations (3.11) and (3.12) have the Lagrangian functional dependence PZP(Z, t, t 0 ) and NZN(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, z Z z adv C z 0 ; ð3:13Þ whence, in (3.11) and (3.12) z adv Z z Kz 0 . 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  At this point, no assumption on the nature of the turbulent field has been made.
To model the vertical turbulent velocity fluctuation w 0 , consider a onedimensional random walk model where k 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 0 from (3.19) is seen as a sum of uncorrelated weighted random increments ðdt 0 Þw 0 . Thus, by the central limit theorem (Papoulis 1965) z Z z C z 0 has the Gaussian distribution ð3:21Þ where the mean displacement z Z z adv is given by (3.10) and the variance b 2 Z hðz 0 Þ 2 i is given by ð3:22aÞ and b 2 Z hðz 0 Þ 2 i Z ðPeÞ K1 ½1KexpfK2atg; Z % 1: ð3:22bÞ Subsequent in the study, we take t 0 Z0. The Peclet number PeZa/kZa Ã Z 2 e /k Ã in (3.22a) and (3.22b) represents the ratio of the turbulent time scale to the advective time scale.
The difference in the form of the variances shown in (3.22a) and (3.22b) 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.22a) and (3.22b), respectively.
For our case we have a turbulent layer of finite domain between zZ0 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) vp L vz Z 0 at z Z 0; z Z 1: ð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

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 P and N , and the turbulence-induced interaction term hP 0 N 0 i, whence with P 0 Z P KP and N 0 Z N KN : Note that the mean of the biodynamical total interaction term F of (3.1) and (3.2) is given by F Z hPN i ZPN C hP 0 N 0 i: ð4:2Þ In the series of figures below, we will show profiles of the biomass quantities given by (4.1a)-(4.1c) 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 ZR1, and a 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 Z0.01 and N 0 Z0.99, which represent an initial seed density of phytoplankton of 1% of the total nutrient and phytoplankton biomass density. We will also use aZ0.1, which, for a nutrient uptake time of 12 hours, corresponds to a linear strain rate of a Ã Z2! 10 K6 s K1 , 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 k Ã in an ocean surface mixed layer ranges from k Ã Z10 K4 to 10 K2 m 2 s K1 (Large et al. 1994). Using a mixed-layer depth of 30 m and an upwelling strain rate of a Ã Z2!10 K6 s K1 , we obtain expected values of Pe(10. Thus, we will use the three values PeZ0, 1 and 10, as well as a fourth value, PeZN; the latter corresponds to the no turbulence pure advection case considered in R99. Note that the PeZ0 case corresponds to infinitely large turbulence relative to the advection.

(a ) Time-dependent solutions
In figure 3a,b we show the solutions of P and N for PeZ10, a moderate turbulence case, for different times t, taking the initial time t 0 Z0. Shown in these figures as a solid line is the steady-state limit t/N. At the bottom boundary, PðzÞj zZ1 and N ðzÞj zZ1 initially increase with time t. Note the discontinuity across zZ1. This results from our model of upwelling into a turbulent mixing region confined to a finite depth range of zZ0 to 1. The no turbulent flux boundary conditions (3.23) insures that both P and N have zero vertical gradients at zZ0 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 tZ lnðN 0 =P 0 Þ, where, for our example problem input parameters, tZ lnðN 0 =P 0 ÞZ lnð99ÞZ 4:6. For theN profile in figure 3b this appears to be the case, but it is clearly not the case for theP profile of figure 3a, which asymptotes to a steady state at the much later time, tz50, not shown in figure 3a. To understand this, note that after a time t[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 timetZ t K t e ( t. The phytoplankton density P has continued to increase beyond tZt 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 z adv Z Z expðKatÞ in (3.25), leads to p E ðZ=z; tÞ/ p E ðz adv =z; tÞ/ 1: ð4:3Þ Using ( where, for this limiting case of infinitely intense turbulence, the mean values shown in (4.4a)-(4.4c) are just the vertically averaged values of the purely advective case. This result, as expected, implies that the intense turbulence case of PeZ0 perfectly mixes P, N and the interaction term (PN ). In figure 4 for PeZ0, we show such plots of P, hðP 0 Þ 2 i 1=2 , N and hðN 0 Þ 2 i 1=2 versus time, where hðP 0 Þ 2 i Z hP 2 iKhPi 2 and hðN 0 Þ 2 i Z hN 2 iKhN i 2 : N and hðN 0 Þ 2 i 1=2 asymptote to a steady state at time of order tZ lnðN 0 =P 0 ÞZ 4:6, a similar time whenN asymptotes to a steady state for the moderate turbulence case shown in figure 3b. From figure 4, we see that the quantitiesP and hðP 0 Þ 2 i 1=2 asymptote to a steady state about an order of magnitude later than N and hðN 0 Þ 2 i 1=2 , again a similar result to the moderate turbulence (PeZ10) case of figure 3a. Note the interesting results of figure 4 that for all time hðN 0 Þ 2 i 1=2 ON and for time t( a K1 ðZ10Þ, hðP 0 Þ 2 i 1=2 OP, while for t[ a K1 ðZ10Þ, hðP 0 Þ 2 i 1=2 z hðN 0 Þ 2 i 1=2 . This latter result follows from the fact that total biomass is conserved and that as t/N P 0 zKN 0 throughout the turbulent layer.

(b ) Steady-state solutions
We show in figure 5 the steady-state (t/N) profiles of mean phytoplankton density P 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 6a,b shows the mean interaction term F and the two terms of which it is composed, PN and hP 0 N 0 i (equation (4.2)). Figure 6a shows the expected behaviour of F as a function of Pe. Decreasing Pe decreases the gradient of F. Note that, in figure 5, the depth of convergence of the solutions of P is about zz0.6. This depth is close to the depth where the mean interaction term, F, shown in figure 6a, is a maximum.
In figure 6b, we show contributions of PN and hP 0 N 0 i to F. For all of the values of Pe, except the non-turbulent value of PeZN, the magnitude of hP 0 N 0 i is of the order PN , and opposite in sign, This result clearly demonstrates that neglecting the turbulent interaction term hP 0 N 0 i would result in an order of magnitude overestimate of F. 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 with F given by (4.1a)-(4.1c) and (4.2). In the extreme limits of Pe/0 and Pe/N, as indicated in figure 7, there is a balance between these two terms This can also be seen in figure 5, where the vertical integral of the light line, the PeZN case, equals that of the vertical integral of the dark line, the PeZ0 case. For the intermediate values of Pe, vertical gradients of P and N 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.

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 F i ½f j ; x; t 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 F i 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. F i sF i ðf j Þ.
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 P and N and the turbulence-induced interaction term hP 0 N 0 i. 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 hP 0 N 0 i 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 PN . 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