- Email: [email protected]
Minimum length scale for growth-limited oceanic plankton distributions
Ecological Modelling 158 (2002) 111 /120 www.elsevier.com/locate/ecolmodel
Minimum length scale for growth-limited oceanic plankton distributions P. McLeod *, A.P. Martin, K.J. Richards School of Ocean and Earth Science, Southampton Oceanography Centre, European Way, Southampton SO14 3ZH, UK Received 12 October 2001; received in revised form 17 May 2002; accepted 28 May 2002
Abstract A patch of phytoplankton in a shear environment will be drawn out into a long ribbon-like structure or filament. Phytoplankton are typically growth-limited due to nutrient limitation or grazing pressure. For such a reactive tracer it is demonstrated that there is a minimum filament width representing a balance between a diffusion and growth fuelled outward-propagating wave and the opposing shear flow. The minimum width is given by sffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k m max 1; 0:4 W 8:6 1:6 1 l l where k is the effective diffusivity for scales small compared to the filament width, l the straining rate of the flow, and m the net tracer growth rate. Two regimes are apparent with behaviour governed by the ratio of biological growth rate to the strain rate. Counter-intuitively, for m /l /2.5, the minimum width exceeds that for an unlimited exponentiallygrowing tracer by a factor of 0.4â(m /l )/0.36. For m /l B/2.5 the behaviour mimics that of an exponentially-growing or inert tracer with filament width independent of growth rate. The results demonstrate that predictions of the finest scale of structure seen in patchy phytoplankton distribution are very sensitive to the biological processes active and their representation. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Plankton; Filaments; Logistic growth; Minimum length scale
1. Introduction Patchiness or heterogeneity is an undisputed, ubiquitous feature of oceanic plankton distributions (e.g. Bainbridge, 1957; Fasham, 1978; Mackas et al., 1985). Scales of observed structure
* Corresponding author. Tel.: /44-23-8059-6631 E-mail address: [email protected] (P. McLeod).
range from centimetres to hundreds of kilometres (Haury et al., 1978). The characteristics of a patch, such as spatial extent, patch intensity and biodiversity, can be a critical factor in determining a population’s dynamics, crucial to the creation and persistence of ecological stability, or the ability to absorb perturbations (Simberloff and Wilson, 1969; Steele, 1974; Smith et al., 1996). A plankton patch under the influence of a straining flow will be drawn out into a long ribbon-like structure or
0304-3800/02/$ - see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 4 - 3 8 0 0 ( 0 2 ) 0 0 1 7 6 - X
112
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
filament. Because of their transient and convoluted nature there have been few studies of the role of planktonic filaments in larger scale ecosystem dynamics. However the expected behavioural regimes of an inert tracer in a mesoscale shear environment have been studied theoretically and applied to in situ observations. Here we show that a biological tracer can exhibit a dramatically different structure to that of an inert tracer. The evolution of an inert tracer field in a straining flow goes through three distinct regimes (Garrett, 1983). For small length scales (L /â(k / l )), where k is the effective diffusivity for small scales and l the straining rate, a patch falls below the influence of the shearing field and disperses through diffusive processes alone. Once the patch is sufficiently large it is influenced by the mesoscale field, which tends to strain the patch into a filament. The filament is expected to be of variance s2 /k /l . For larger scales streaks merge and again diffusive processes become primarily responsible for the growth of the patch. The predicted time taken for merging or homogenisation of the tracer field depends upon both the assumed large scale and small-scale diffusive parameters and the local straining rate. Inert tracer release experiments (Ledwell et al., 1998; Sundermeyer and Price, 1998) used the expected variance of filaments, s2 /k /l, and observed straining rate to provide an estimate of the effective diffusion coefficient. The first study to address the competition between growth and diffusion in a pure straining flow was Martin (2000). An exponentially-growing tracer in a shear environment converges to a Gaussian of the same variance as an inert tracer, s2 /k /l , (Martin, 2000). This prediction of a minimum width for phytoplankton filaments has been verified by recent iron fertilisation experiments (Boyd et al., 2000) which produced a plankton patch that was drawn into a filament, visible in satellite imagery of the area (Abraham et al., 2000). Assuming the phytoplankton behaved as an exponentially-growing tracer the same approach as Ledwell et al. (1998) was used to calculate a value for effective diffusivity for the region (Abraham et al., 2000). Caution must be urged with this approach as the derived quantities
depend greatly upon the assumed growth process and the consequent final structure. Unlimited growth results in an infinitely increasing phytoplankton concentration whereas in reality processes such as nutrient limitation and predator grazing will limit the magnitude of any plankton community. We investigate the influence of this limitation. As we will show a minor change to the equation describing growth of the reactive constituent can have a dramatic impact upon the observed tracer structure and its time-varying behaviour.
2. Methods Consider the evolution of the distribution of a reactive tracer, such as phytoplankton, in a fluid environment. We describe this by: DP P k92 PmP 1 (1) DT P0 where P is the tracer concentration, D/DT is the total derivative, k is an effective diffusivity which parameterises physical processes at scales smaller than the model resolution, and m is the maximum net growth rate of the tracer. The final term in Eq. (1) is a greatly simplified parameterisation of biological processes. Otherwise known as logistic growth it has a threshold limit, P0, to which the population will tend and above which the population is unable to sustain itself. This is the simplest parameterisation of growth limiting factors such as grazing and nutrient limitation. The growth rate, m, is considered homogeneous, i.e. any variation, either spatially or temporally, is assumed to be at scales larger than the structure produced by the flow. It should be noted that we only take ambient waters to have zero concentration of phytoplankton for convenience. Simply replacing the growth term P mP 1 P0 with
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
mPb
P 1 Pb
1
P P0
allows any background concentration, Pb, without changing any of the dynamics or following results. We will assume the flow is two-dimensional and divergence free, a reasonable assumption away from strong fronts and at scales greater than /1 km. Any such flow can be resolved locally into a pure strain and a rotation. Straining processes are primarily responsible for the formation of filaments (Ottino, 1989). Hence we restrict attention to a pure straining flow, u /(/lx , ly ), where l is the straining rate. In such a flow field the tracer will be drawn out in the y direction into a ribbonlike structure or filament. Hence a one-dimensional approximation can be justified, as the across-filament gradients (in the x direction) will far exceed those along the filament. To further aid analysis Eq. (1) can be non-dimensionalised firstly by normalising the concentration with respect to the population limit, C P=P0 ; and secondly by pffiffiffiffiffiffiffiffiffiffiffi substituting t(1=l)t; and x (k=l)f:/ Thus @C @C @ 2 C f bC(1C): @t @f @f2
(2)
The remaining parameter, b/m/l is an inverse Damko¨hler number (Damko¨hler and Heumann, 1982) being a measure of the ratio of reaction to flow time scales. Typical net growth and strain rates in the ocean both vary in the range 0.01 /1 day1 (Sundermeyer and Price, 1998; Abraham et al., 2000). We therefore choose to examine b in the range 0.01 /100. From now on an asterisk will be used to denote dimensional parameters. Our model is a combination of two well-studied systems. As such one may expect the solution of Eq. (2) to tend towards the solution of each component in certain limits. Without advection Eq. (2) reduces to the wellknown Fisher’s equation (Murray, 1993). Fisher’s equation, which combines diffusion with a logistic growth term, possesses travelling wave solutions, which can take the form of propagating fronts (Fisher, 1937). For an initially localised ‘seed’ population these fronts have a minimum velocity,
113
Umin /2âb (U min /2â(km ), in dimensional form) (see e.g. Murray, 1993). The velocity of these fronts increases with decreasing steepness of the front. We conjecture that with the inclusion of an opposing advective velocity the fronts will stop where the opposing flow is of similar magnitude. Considering the flow field u /f(u* /lx ) and frontal velocity UF /âb , and assuming the fronts propagate symmetrically from the origin this results in a final expected non-dimensional width of the filament: W
pffiffiffiffi b
(3)
pffiffiffiffiffiffi km W l The functional form of this width differs significantly from that for an exponentially-growing (or inert) tracer, â(k /l ). The final width here is dependent on the growth rate of the tracer and is more sensitive to straining rate. For small C Eq. (2) tends to the Martin (2000) model. There the use of a growth term without logistic correction results in a Gaussian of width â(k /l) (/1 in model dimensions), the magnitude of which increases or decreases exponentially with time depending of the relative magnitude of the growth and shear rates (Martin, 2000). For small b stretching should dominate, preventing the concentration reaching its limiting value of 1, and the solution is expected to tend towards that for an inert or exponentially-growing tracer i.e. a Gaussian of width â(k /l). For large b biological effects should dominate and the solution is expected to tend to that for Fishers equation. To find the limits on the behaviour of the system we solve Eq. (2) numerically. The chosen advection scheme is flux-corrected upstream differencing (Smolarkiewicz, 1984). To avoid stability problems and the question of at which point(s) in a time step growth occurs the physical (advection and diffusion) and reactive (growth) components are treated separately (Rood, 1987). The reactive component is solved by a 4th order Runge /Kutta scheme.
114
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
3. Results The inclusion of the factor (1/C ) in the growth term radically changes the behaviour of the system compared to one with unlimited exponential growth. As expected the behaviour of the system is very dependent on the value of b. Fig. 1 shows the tracer distribution for four different values of b which characterise the different regimes of the system. Henceforth we define the width, W , as the distance between the ‘noses’ of the fronts. We define the ‘nose’ as 0.01% of the maximum tracer concentration as it gives the most consistent results. In Fig. 2(b) we show the variation of final width with parameter value.
For values of b 0/1 (Fig. 1(a)) the distribution decays to zero as a Gaussian of variance 1. Consequently the width of the distributions is constant for b B/1 (Fig. 2(b)) i.e. the width is independent of growth rate, reproducing the result of Martin (2000) for an exponentiallygrowing tracer in a convergent flow field. In this regime the suppression of biological effects by the flow is sufficient that the tracer concentration never reaches the threshold where fronts begin to form. For b E/10 (Fig. 1(c) and (d)) a propagating front develops. The nose of the front stops where the advective current balances the propagation velocity (Fig. 2(a)). As expected the asymptotic
Fig. 1. Evolution of a reactive tracer, in a convergent flow field, from an initial ‘tanh’ function. Results for (a) b /0.01, (b) b /2, (c) b /10, (d) b /50. The arrows show direction of movement with time, t . The solution at large time is shown in bold.
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
115
Fig. 2. Asymptotic variation with growth rate of (a) Velocity of fronts, /4.3*max (1, 0.4âb/0.36), (b) Width, /8.6*max (1, 0.4âb/0.36), (c) Gradient of fronts.
filament width varies as the square root of the parameter b (Fig. 2(b)). However there is an additional constant correction to Eq. (3) suggesting that the frontal behaviour is modified by the
advective field beyond simple additive effects. The asymptotic frontal gradient is seen to increase with b, Fig. 2(c). We suggest this is because for a higher growth rate the population is better able to sustain
116
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
itself closer to the stopping point. For a lower growth rate the current has a greater effect on the population ‘pushing’ it away from the equilibrium point resulting in a gentler front. There is a transition region where the asymptotic solution is neither exactly Gaussian nor a propagating front. We choose to differentiate between the Gaussian and propagating front regimes by using Kurtosis, K m4 =s4 ; a measure of ‘peakiness’ of a distribution, where m4 is the fourth moment about the mean and s2 is the distribution variance. The Kurtosis of a Gaussian is 3. In its extreme a frontal solution may be thought of as a box or double step function, which has K /1.79. Here a frontal solution is defined by its Kurtosis being within 10% of that of an extreme frontal solution, i.e. K B/1.97. The Kurtosis of the steady state solution is shown in Fig. 3(b) as a function of b. The three zones may also be distinguished through their maximum asymptotic height, Fig. 3(a). A frontal solution has a limiting concentration $/1 whilst Gaussian solutions tend to zero. Referring to Fig. 3(a) and (b) the transition regime therefore spans 10/b 0/10 and 1.8 B/Kurtosis B/2.7, where the maximum concentration of the distribution neither decays to zero nor reaches its limiting value. In its steady state the concentration reaches a limit where, at the centre of the filament, net growth is matched by losses due to diffusion. For large values of b we find a fourth regime. The steady state solution is the same as for lower values of b (/10). However, there exists a bc such that for b /bc instead of stopping directly at the steady state point fronts ‘overshoot’ before subsequently converging. The precise value of bc is initial condition dependent although it increases as the difference between the initial condition and the steady state solution decreases. The magnitude of overshoot increases generally with b for a range of given initial condition. Results (not shown) suggest that for a sufficiently high growth rate gentle fronts propagate outwards faster than the advective field is able to adjust the frontal gradient and hence propagation velocity. The fronts eventually become steeper, reducing the frontal velocity with the consequence that the filamental width de-
creases until the expected width, at the velocity balance, is reached. The velocity and width distributions, Fig. 2(a) and (b), have a two-part linear regression: pffiffiffiffi (4) Vf 4:3 max(1; 0:4( b 1:6)1) pffiffiffiffi W 8:6 max(1; 0:4( b 1:6)1) (5) sffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k m W 8:6 max 1; 0:4 1:6 1 l l A very good fit to data is observed. The nonunity coefficients are dependent on the definition of width. Our formula for the width changes from constant to increasing as âb at the point where b/1/o , and o :/1.5. In the transition regime the asymptotic width is under predicted, but the maximum error is less than 10% of the predicted width. The time taken to reach steady state is also of interest. The behaviour of the concentration at the centre of the filament can be modelled by: t C(t)C0 C?exp (6) h where C0 is the final concentration and h the decay timescale is given to good approximation by h 1=(jb1j) (Fig. 3(c)). The apparent singularity at b/1 shows the transition between an exponentially decaying solution and a non-zero steady state being achieved. The time for the filament width to reach its equilibrium value is more difficult to quantify. For b 0/1 the time scale is 0.5, or a decay rate of 2. This agrees once again with that found by Martin (2000) for an exponentially-growing tracer. For b /1 the decay timescale is initial condition dependent. The time taken for the distribution to achieve its steady state solution decreases with increased proximity of the initial condition to the asymptotic solution. We are interested in the longest expected convergence time for a plankton distribution. We therefore choose to describe convergence times for the smallest practicable initial condition (width two boxes, Cmax /0.01). For b B/5 (approximately Gaussian) steady state is achieved within /0.28 model time units. Assuming l /1 /105 s 1 this corresponds to /0.5
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
117
Fig. 3. (a) Variation of maximum asymptotic tracer concentration, Cmax(T ), with parameter value, (b) Kurtosis, Gaussian/3, propagating front/1.8, (c) Variation of (concentrationwise) decay timescale, h , with parameter value. Solid line indicates 1 : Note log scale. approximation h jb1j
118
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
day. For the propagating front solutions the mean convergence time is 0.22 model time units. Assuming l /1/10 5 s 1 again this corresponds to / 0.25 day. However l/105 s 1 is the predicted maximum straining rate. For the smallest expected strain rate, l /10 7 s 1, this would be 25 days, but a mean value is 3 days corresponding to l/ 106 s 1. Both mentioned time scales are dependent on the assumed shearing rate. As the Lagrangian integral timescale for a mesoscale ocean flow varies from 2 to 10 days (Griffa et al., 1995), this suggests that not all filaments will achieve the predicted steady state width during the life span of the shearing process.
4. Discussion The predicted width of a reactive tracer differs from that predicted for an inert or exponentiallygrowing tracer by LR/LI /max(1, 0.4âb/0.36) / G(b ), where LR is the width of a reactive tracer filament and LI the width of an inert filament. For a given l using the expected width of an inert filament to calculate k for a growth-limited reactive tracer leads to an overestimate of k by a factor of G2 for b /2.5. Choosing typical values of k /1 m2 s 1 and l /106 s 1 results in a predicted inert filament width of 8.6 km. The predicted width of a growthlimited reactive filament, in the same environment, is the greater of 8.6 km or 3.44â(m /l )/3.1 km. In practice filament width is limited above by oceanic integral length scales which are of the order 37 km zonally and 28 km meridionally (Schafer and Krauss, 1995). However for the values above to increase the reactive filament to the order of the integral length scale of the mesoscale field requires b /62. The largest discrepancy between exponentially-growing and growth-limited phytoplankton that could be expected will occur in bloom conditions in low energy regions. Here we may expect b /20 giving a factor of two difference. Although b /20 may be realistic for short periods, such values are unlikely to be sustained for periods larger than a few weeks as net growth rates rapidly decrease with the onset of decreased nutrient availability and increased predation.
It should be remembered that effective diffusivity varies with scale. If we use Okubo’s (1971) empirical formula k /0.0103L1.15 cm2 s 1, with growth and shearing rates of similar order such that b /3, and m and l are in the range 10 5 / 107 s 1 then our formula predicts filament widths in the range 0.8 /80 km accordingly. It should be noted that all estimates of filament width using Okubo’s formula have an associated error of up to an order of magnitude due to the sensitivity of the parameterisation of diffusivity to sea and weather conditions (Okubo, 1971). To date the only application of the expected width of a tracer filament to observational data has been to derive a value for the local effective diffusivity. The expected width of an exponentially-growing tracer was applied to a plankton filament formed during the recent SOIREE experiment (Boyd et al., 2000; Abraham et al., 2000). The values described by Abraham et al. suggest b/1.59/1. The formulae for growth-limited and unlimited exponentially-growing tracers do not differ for this value. However caution must be urged as for higher values of b the predicted solutions differ dramatically. Comparison of results through derived values of diffusivity is not recommended as there is no independent means of assessing the correct value. A better means of comparison of model results is through the structure or Kurtosis of the filament cross-section. However this requires high-resolution information on tracer concentration, from low-level remote sensing or fine scale in situ surveys. Although such data is currently unavailable, it should be available in the near future, due to developments in instrumentation and the increasing frequency of multi-disciplinary cruises. Subsequent to completion of this work and submission of this manuscript we were made aware of a manuscript (Neufeld et al., submitted for publication) investigating a similar problem from a more mathematical viewpoint. Neufeld et al. (submitted for publication) also predict the existence of a regime where tracer reaction (there auto-catalysis) can affect filamental width and derive an estimate for the width from a dimensional argument that is consistent with our Eq. (7).
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
5. Conclusions The structure observed in a reactive tracer field, such as plankton, is often markedly different from that of an inert tracer, (Denman and Platt, 1976; Seuront et al., 1999). We have shown that one reason may be a different functional form for the expected asymptotic width of any filaments formed. sffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k m W 8:6 max 1; 0:4 1:6 1 (7) l l For relatively large growth rate or small strain rate, such that m/l /16, the predicted structure differs significantly from that predicted for an inert or exponentially-growing tracer. Observed mesoscale filament widths are expected to be less than 1 km up to the integral length scale of the straining field. These filaments are expected to have a mean convergence time of /3 days. Numerical results suggest larger filaments may be formed but there are factors likely to disrupt filaments on this scale: grazing modifying the distribution; variation of the current field and growth rate on temporal and spatial scales less than those required for the asymptotic solution to be achieved. It should be noted that using a more sophisticated growth model, specifically explicitly modelling nitrate instead of using a logistic growth term, gives near identical results (McLeod, pers. comm.). An interesting feature of our results is that the biological processes assumed can have a marked effect on the filamental structure. This is in contrast to the case of an exponentially-growing tracer where filament width is independent of growth rate. Perhaps counter-intuitively imposing a limit to the population such that fronts of plankton propagate out from the centre of a convergent flow, resulting in a wider final distribution than if growth was unrestricted. The imminent availability of sufficiently high-resolution remote sensed and in situ data for the ecological structure of planktonic filaments will allow this hypothesis to be tested. Meanwhile explicit simulation of such tracers in two-dimensional turbulent fields and the explicit inclusion of
119
processes such as nutrient limitation and grazing provide ample scope for investigators to improve our scanty understanding of the ubiquitous plankton filament.
Acknowledgements This work was funded by a NERC studentship (GT04/99/MS/262).
References Abraham, E.R., Law, C.S., Boyd, P.W., Lavender, S.J., Bowie, M.T., 2000. Importance of stirring in the development of an iron-fertilized phytoplankton bloom. Nature 407, 727 /733. Bainbridge, R., 1957. The size shape and density of marine phytoplankton concentrations. Biological Review 32, 91 / 115. Boyd, P.W., Watson, A.J., Law, C.S., Abraham, E.R., Trull, T., Murdoch, R., et al., 2000. A mesoscale phytoplankton bloom in the polar southern ocean stimulated by iron fertilization. Nature 407, 695 /702. Damko¨hler, R., Heumann, T., 1982. Experimental studies to clarify the diffusion behavior of copper /nickel /alloys with a high copper content. Physica Status Solidi A */Applied Research 73 (1), 117 /127. Denman, K.L., Platt, T., 1976. The variance spectrum of phytoplankton in a turbulent ocean. Journal of Marine Research 34 (4), 593 /601. Fasham, M.J.R., 1978. The application of some stochastic processes to the study of plankton patchiness. In: Steele, J.H. (Ed.), Spatial Patterns in Plankton Communities. New York, Plenum Press, pp. 131 /156. Fisher, R.A., 1937. The wave advance of advantageous genes. Annals of Eugenics 7, 355 /369. Garrett, C., 1983. On the initial streakiness of a dispersing tracer in two- and three-dimensional turbulence. Dynamics of Atmospheres and Oceans 7, 265 /277. Griffa, A., Owens, K., Piterbarg, L., Rozovskii, B., 1995. Estimates of turbulence parameters from Lagrangian data using stochastic particle model. Journal of Marine Research 53, 371 /401. Haury, L.R., McGowan, J.A., Weibe, P.H., 1978. Patterns and processes in the time-space scales of plankton. In: Steele, J.H. (Ed.), Spatial Patterns in Plankton Communities. Plenum Press, New York, pp. 277 /327. Ledwell, J.R., Watson, A.J., Law, C.S., 1998. Mixing of a tracer in the pycnocline. Journal of Geophysical Research 103 (C10), 21499 /21529. Mackas, D.L., Denman, K.L., Abbott, M.R., 1985. Plankton patchiness: biology in the physical vernacular. Bulletin of Marine Science 37 (2), 652 /674.
120
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
Martin, A.P., 2000. On filament widths in oceanic plankton distributions. Journal of Plankton Research 22 (3), 597 / 602. Murray, J.D., 1993. Biological waves: single species models. In: Mathematical Biology. New York, Springer, pp. 274 /310. Neufeld, Z., Haynes, P.H., Te´l, T., Chaotic mixing induced transitions in reaction-diffusion systems, Chaos, Submitted for publication. Okubo, A., 1971. Oceanic diffusion diagrams. Deep-Sea Research 18, 789 /802. Ottino, J.M., 1989. The Kinematics of Mixing: Stretching, Chaos and Transport. Cambridge university press. Rood, R.B., 1987. Numerical advection algorithms and their role in atmospheric transport and chemistry models. Reviews of Geophysics 25 (1), 71 /100. Schafer, H., Krauss, W., 1995. Eddy statistics in the South Atlantic as derived from drifters drouged at 100 m. Journal of Marine Research 53, 403 /431. Seuront, L., Schmitt, F., Lagadeuc, Y., Schertzer, D., Lovejoy, S., 1999. Universal multifractal analysis as a tool to
characterize multiscale intermittent patterns: examples of phytoplankton distribution in turbulent coastal waters. Journal of Plankton Research 21 (5), 877 /922. Simberloff, D.S., Wilson, E.O., 1969. Experimental zoogeography of islands: the colonization of empty islands. Ecology 50, 278 /296. Smith, C.L., Richards, K.J., Fasham, M.J.R., 1996. The impact of mesoscale eddies on plankton dynamics in the upper ocean. Deep-Sea Research I 43 (11-12), 1807 /1832. Smolarkiewicz, P.K., 1984. A fully multidimensional positive definite advection transport algorithm with small implicit diffusion. Journal of Computational Physics 54, 325 /362. Steele, J.H., 1974. Spatial heterogeneity and population stability. Nature 248, 83. Sundermeyer, M.A., Price, J.F., 1998. Lateral mixing and the North Atlantic tracer release experiment: observations and numerical simulations of Lagrangian particles and a passive tracer. Journal of Geophysical Research 103 (C10), 21481 / 21497.
Minimum length scale for growth-limited oceanic plankton distributions P. McLeod *, A.P. Martin, K.J. Richards School of Ocean and Earth Science, Southampton Oceanography Centre, European Way, Southampton SO14 3ZH, UK Received 12 October 2001; received in revised form 17 May 2002; accepted 28 May 2002
Abstract A patch of phytoplankton in a shear environment will be drawn out into a long ribbon-like structure or filament. Phytoplankton are typically growth-limited due to nutrient limitation or grazing pressure. For such a reactive tracer it is demonstrated that there is a minimum filament width representing a balance between a diffusion and growth fuelled outward-propagating wave and the opposing shear flow. The minimum width is given by sffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k m max 1; 0:4 W 8:6 1:6 1 l l where k is the effective diffusivity for scales small compared to the filament width, l the straining rate of the flow, and m the net tracer growth rate. Two regimes are apparent with behaviour governed by the ratio of biological growth rate to the strain rate. Counter-intuitively, for m /l /2.5, the minimum width exceeds that for an unlimited exponentiallygrowing tracer by a factor of 0.4â(m /l )/0.36. For m /l B/2.5 the behaviour mimics that of an exponentially-growing or inert tracer with filament width independent of growth rate. The results demonstrate that predictions of the finest scale of structure seen in patchy phytoplankton distribution are very sensitive to the biological processes active and their representation. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Plankton; Filaments; Logistic growth; Minimum length scale
1. Introduction Patchiness or heterogeneity is an undisputed, ubiquitous feature of oceanic plankton distributions (e.g. Bainbridge, 1957; Fasham, 1978; Mackas et al., 1985). Scales of observed structure
* Corresponding author. Tel.: /44-23-8059-6631 E-mail address: [email protected] (P. McLeod).
range from centimetres to hundreds of kilometres (Haury et al., 1978). The characteristics of a patch, such as spatial extent, patch intensity and biodiversity, can be a critical factor in determining a population’s dynamics, crucial to the creation and persistence of ecological stability, or the ability to absorb perturbations (Simberloff and Wilson, 1969; Steele, 1974; Smith et al., 1996). A plankton patch under the influence of a straining flow will be drawn out into a long ribbon-like structure or
0304-3800/02/$ - see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 4 - 3 8 0 0 ( 0 2 ) 0 0 1 7 6 - X
112
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
filament. Because of their transient and convoluted nature there have been few studies of the role of planktonic filaments in larger scale ecosystem dynamics. However the expected behavioural regimes of an inert tracer in a mesoscale shear environment have been studied theoretically and applied to in situ observations. Here we show that a biological tracer can exhibit a dramatically different structure to that of an inert tracer. The evolution of an inert tracer field in a straining flow goes through three distinct regimes (Garrett, 1983). For small length scales (L /â(k / l )), where k is the effective diffusivity for small scales and l the straining rate, a patch falls below the influence of the shearing field and disperses through diffusive processes alone. Once the patch is sufficiently large it is influenced by the mesoscale field, which tends to strain the patch into a filament. The filament is expected to be of variance s2 /k /l . For larger scales streaks merge and again diffusive processes become primarily responsible for the growth of the patch. The predicted time taken for merging or homogenisation of the tracer field depends upon both the assumed large scale and small-scale diffusive parameters and the local straining rate. Inert tracer release experiments (Ledwell et al., 1998; Sundermeyer and Price, 1998) used the expected variance of filaments, s2 /k /l, and observed straining rate to provide an estimate of the effective diffusion coefficient. The first study to address the competition between growth and diffusion in a pure straining flow was Martin (2000). An exponentially-growing tracer in a shear environment converges to a Gaussian of the same variance as an inert tracer, s2 /k /l , (Martin, 2000). This prediction of a minimum width for phytoplankton filaments has been verified by recent iron fertilisation experiments (Boyd et al., 2000) which produced a plankton patch that was drawn into a filament, visible in satellite imagery of the area (Abraham et al., 2000). Assuming the phytoplankton behaved as an exponentially-growing tracer the same approach as Ledwell et al. (1998) was used to calculate a value for effective diffusivity for the region (Abraham et al., 2000). Caution must be urged with this approach as the derived quantities
depend greatly upon the assumed growth process and the consequent final structure. Unlimited growth results in an infinitely increasing phytoplankton concentration whereas in reality processes such as nutrient limitation and predator grazing will limit the magnitude of any plankton community. We investigate the influence of this limitation. As we will show a minor change to the equation describing growth of the reactive constituent can have a dramatic impact upon the observed tracer structure and its time-varying behaviour.
2. Methods Consider the evolution of the distribution of a reactive tracer, such as phytoplankton, in a fluid environment. We describe this by: DP P k92 PmP 1 (1) DT P0 where P is the tracer concentration, D/DT is the total derivative, k is an effective diffusivity which parameterises physical processes at scales smaller than the model resolution, and m is the maximum net growth rate of the tracer. The final term in Eq. (1) is a greatly simplified parameterisation of biological processes. Otherwise known as logistic growth it has a threshold limit, P0, to which the population will tend and above which the population is unable to sustain itself. This is the simplest parameterisation of growth limiting factors such as grazing and nutrient limitation. The growth rate, m, is considered homogeneous, i.e. any variation, either spatially or temporally, is assumed to be at scales larger than the structure produced by the flow. It should be noted that we only take ambient waters to have zero concentration of phytoplankton for convenience. Simply replacing the growth term P mP 1 P0 with
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
mPb
P 1 Pb
1
P P0
allows any background concentration, Pb, without changing any of the dynamics or following results. We will assume the flow is two-dimensional and divergence free, a reasonable assumption away from strong fronts and at scales greater than /1 km. Any such flow can be resolved locally into a pure strain and a rotation. Straining processes are primarily responsible for the formation of filaments (Ottino, 1989). Hence we restrict attention to a pure straining flow, u /(/lx , ly ), where l is the straining rate. In such a flow field the tracer will be drawn out in the y direction into a ribbonlike structure or filament. Hence a one-dimensional approximation can be justified, as the across-filament gradients (in the x direction) will far exceed those along the filament. To further aid analysis Eq. (1) can be non-dimensionalised firstly by normalising the concentration with respect to the population limit, C P=P0 ; and secondly by pffiffiffiffiffiffiffiffiffiffiffi substituting t(1=l)t; and x (k=l)f:/ Thus @C @C @ 2 C f bC(1C): @t @f @f2
(2)
The remaining parameter, b/m/l is an inverse Damko¨hler number (Damko¨hler and Heumann, 1982) being a measure of the ratio of reaction to flow time scales. Typical net growth and strain rates in the ocean both vary in the range 0.01 /1 day1 (Sundermeyer and Price, 1998; Abraham et al., 2000). We therefore choose to examine b in the range 0.01 /100. From now on an asterisk will be used to denote dimensional parameters. Our model is a combination of two well-studied systems. As such one may expect the solution of Eq. (2) to tend towards the solution of each component in certain limits. Without advection Eq. (2) reduces to the wellknown Fisher’s equation (Murray, 1993). Fisher’s equation, which combines diffusion with a logistic growth term, possesses travelling wave solutions, which can take the form of propagating fronts (Fisher, 1937). For an initially localised ‘seed’ population these fronts have a minimum velocity,
113
Umin /2âb (U min /2â(km ), in dimensional form) (see e.g. Murray, 1993). The velocity of these fronts increases with decreasing steepness of the front. We conjecture that with the inclusion of an opposing advective velocity the fronts will stop where the opposing flow is of similar magnitude. Considering the flow field u /f(u* /lx ) and frontal velocity UF /âb , and assuming the fronts propagate symmetrically from the origin this results in a final expected non-dimensional width of the filament: W
pffiffiffiffi b
(3)
pffiffiffiffiffiffi km W l The functional form of this width differs significantly from that for an exponentially-growing (or inert) tracer, â(k /l ). The final width here is dependent on the growth rate of the tracer and is more sensitive to straining rate. For small C Eq. (2) tends to the Martin (2000) model. There the use of a growth term without logistic correction results in a Gaussian of width â(k /l) (/1 in model dimensions), the magnitude of which increases or decreases exponentially with time depending of the relative magnitude of the growth and shear rates (Martin, 2000). For small b stretching should dominate, preventing the concentration reaching its limiting value of 1, and the solution is expected to tend towards that for an inert or exponentially-growing tracer i.e. a Gaussian of width â(k /l). For large b biological effects should dominate and the solution is expected to tend to that for Fishers equation. To find the limits on the behaviour of the system we solve Eq. (2) numerically. The chosen advection scheme is flux-corrected upstream differencing (Smolarkiewicz, 1984). To avoid stability problems and the question of at which point(s) in a time step growth occurs the physical (advection and diffusion) and reactive (growth) components are treated separately (Rood, 1987). The reactive component is solved by a 4th order Runge /Kutta scheme.
114
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
3. Results The inclusion of the factor (1/C ) in the growth term radically changes the behaviour of the system compared to one with unlimited exponential growth. As expected the behaviour of the system is very dependent on the value of b. Fig. 1 shows the tracer distribution for four different values of b which characterise the different regimes of the system. Henceforth we define the width, W , as the distance between the ‘noses’ of the fronts. We define the ‘nose’ as 0.01% of the maximum tracer concentration as it gives the most consistent results. In Fig. 2(b) we show the variation of final width with parameter value.
For values of b 0/1 (Fig. 1(a)) the distribution decays to zero as a Gaussian of variance 1. Consequently the width of the distributions is constant for b B/1 (Fig. 2(b)) i.e. the width is independent of growth rate, reproducing the result of Martin (2000) for an exponentiallygrowing tracer in a convergent flow field. In this regime the suppression of biological effects by the flow is sufficient that the tracer concentration never reaches the threshold where fronts begin to form. For b E/10 (Fig. 1(c) and (d)) a propagating front develops. The nose of the front stops where the advective current balances the propagation velocity (Fig. 2(a)). As expected the asymptotic
Fig. 1. Evolution of a reactive tracer, in a convergent flow field, from an initial ‘tanh’ function. Results for (a) b /0.01, (b) b /2, (c) b /10, (d) b /50. The arrows show direction of movement with time, t . The solution at large time is shown in bold.
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
115
Fig. 2. Asymptotic variation with growth rate of (a) Velocity of fronts, /4.3*max (1, 0.4âb/0.36), (b) Width, /8.6*max (1, 0.4âb/0.36), (c) Gradient of fronts.
filament width varies as the square root of the parameter b (Fig. 2(b)). However there is an additional constant correction to Eq. (3) suggesting that the frontal behaviour is modified by the
advective field beyond simple additive effects. The asymptotic frontal gradient is seen to increase with b, Fig. 2(c). We suggest this is because for a higher growth rate the population is better able to sustain
116
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
itself closer to the stopping point. For a lower growth rate the current has a greater effect on the population ‘pushing’ it away from the equilibrium point resulting in a gentler front. There is a transition region where the asymptotic solution is neither exactly Gaussian nor a propagating front. We choose to differentiate between the Gaussian and propagating front regimes by using Kurtosis, K m4 =s4 ; a measure of ‘peakiness’ of a distribution, where m4 is the fourth moment about the mean and s2 is the distribution variance. The Kurtosis of a Gaussian is 3. In its extreme a frontal solution may be thought of as a box or double step function, which has K /1.79. Here a frontal solution is defined by its Kurtosis being within 10% of that of an extreme frontal solution, i.e. K B/1.97. The Kurtosis of the steady state solution is shown in Fig. 3(b) as a function of b. The three zones may also be distinguished through their maximum asymptotic height, Fig. 3(a). A frontal solution has a limiting concentration $/1 whilst Gaussian solutions tend to zero. Referring to Fig. 3(a) and (b) the transition regime therefore spans 10/b 0/10 and 1.8 B/Kurtosis B/2.7, where the maximum concentration of the distribution neither decays to zero nor reaches its limiting value. In its steady state the concentration reaches a limit where, at the centre of the filament, net growth is matched by losses due to diffusion. For large values of b we find a fourth regime. The steady state solution is the same as for lower values of b (/10). However, there exists a bc such that for b /bc instead of stopping directly at the steady state point fronts ‘overshoot’ before subsequently converging. The precise value of bc is initial condition dependent although it increases as the difference between the initial condition and the steady state solution decreases. The magnitude of overshoot increases generally with b for a range of given initial condition. Results (not shown) suggest that for a sufficiently high growth rate gentle fronts propagate outwards faster than the advective field is able to adjust the frontal gradient and hence propagation velocity. The fronts eventually become steeper, reducing the frontal velocity with the consequence that the filamental width de-
creases until the expected width, at the velocity balance, is reached. The velocity and width distributions, Fig. 2(a) and (b), have a two-part linear regression: pffiffiffiffi (4) Vf 4:3 max(1; 0:4( b 1:6)1) pffiffiffiffi W 8:6 max(1; 0:4( b 1:6)1) (5) sffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k m W 8:6 max 1; 0:4 1:6 1 l l A very good fit to data is observed. The nonunity coefficients are dependent on the definition of width. Our formula for the width changes from constant to increasing as âb at the point where b/1/o , and o :/1.5. In the transition regime the asymptotic width is under predicted, but the maximum error is less than 10% of the predicted width. The time taken to reach steady state is also of interest. The behaviour of the concentration at the centre of the filament can be modelled by: t C(t)C0 C?exp (6) h where C0 is the final concentration and h the decay timescale is given to good approximation by h 1=(jb1j) (Fig. 3(c)). The apparent singularity at b/1 shows the transition between an exponentially decaying solution and a non-zero steady state being achieved. The time for the filament width to reach its equilibrium value is more difficult to quantify. For b 0/1 the time scale is 0.5, or a decay rate of 2. This agrees once again with that found by Martin (2000) for an exponentially-growing tracer. For b /1 the decay timescale is initial condition dependent. The time taken for the distribution to achieve its steady state solution decreases with increased proximity of the initial condition to the asymptotic solution. We are interested in the longest expected convergence time for a plankton distribution. We therefore choose to describe convergence times for the smallest practicable initial condition (width two boxes, Cmax /0.01). For b B/5 (approximately Gaussian) steady state is achieved within /0.28 model time units. Assuming l /1 /105 s 1 this corresponds to /0.5
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
117
Fig. 3. (a) Variation of maximum asymptotic tracer concentration, Cmax(T ), with parameter value, (b) Kurtosis, Gaussian/3, propagating front/1.8, (c) Variation of (concentrationwise) decay timescale, h , with parameter value. Solid line indicates 1 : Note log scale. approximation h jb1j
118
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
day. For the propagating front solutions the mean convergence time is 0.22 model time units. Assuming l /1/10 5 s 1 again this corresponds to / 0.25 day. However l/105 s 1 is the predicted maximum straining rate. For the smallest expected strain rate, l /10 7 s 1, this would be 25 days, but a mean value is 3 days corresponding to l/ 106 s 1. Both mentioned time scales are dependent on the assumed shearing rate. As the Lagrangian integral timescale for a mesoscale ocean flow varies from 2 to 10 days (Griffa et al., 1995), this suggests that not all filaments will achieve the predicted steady state width during the life span of the shearing process.
4. Discussion The predicted width of a reactive tracer differs from that predicted for an inert or exponentiallygrowing tracer by LR/LI /max(1, 0.4âb/0.36) / G(b ), where LR is the width of a reactive tracer filament and LI the width of an inert filament. For a given l using the expected width of an inert filament to calculate k for a growth-limited reactive tracer leads to an overestimate of k by a factor of G2 for b /2.5. Choosing typical values of k /1 m2 s 1 and l /106 s 1 results in a predicted inert filament width of 8.6 km. The predicted width of a growthlimited reactive filament, in the same environment, is the greater of 8.6 km or 3.44â(m /l )/3.1 km. In practice filament width is limited above by oceanic integral length scales which are of the order 37 km zonally and 28 km meridionally (Schafer and Krauss, 1995). However for the values above to increase the reactive filament to the order of the integral length scale of the mesoscale field requires b /62. The largest discrepancy between exponentially-growing and growth-limited phytoplankton that could be expected will occur in bloom conditions in low energy regions. Here we may expect b /20 giving a factor of two difference. Although b /20 may be realistic for short periods, such values are unlikely to be sustained for periods larger than a few weeks as net growth rates rapidly decrease with the onset of decreased nutrient availability and increased predation.
It should be remembered that effective diffusivity varies with scale. If we use Okubo’s (1971) empirical formula k /0.0103L1.15 cm2 s 1, with growth and shearing rates of similar order such that b /3, and m and l are in the range 10 5 / 107 s 1 then our formula predicts filament widths in the range 0.8 /80 km accordingly. It should be noted that all estimates of filament width using Okubo’s formula have an associated error of up to an order of magnitude due to the sensitivity of the parameterisation of diffusivity to sea and weather conditions (Okubo, 1971). To date the only application of the expected width of a tracer filament to observational data has been to derive a value for the local effective diffusivity. The expected width of an exponentially-growing tracer was applied to a plankton filament formed during the recent SOIREE experiment (Boyd et al., 2000; Abraham et al., 2000). The values described by Abraham et al. suggest b/1.59/1. The formulae for growth-limited and unlimited exponentially-growing tracers do not differ for this value. However caution must be urged as for higher values of b the predicted solutions differ dramatically. Comparison of results through derived values of diffusivity is not recommended as there is no independent means of assessing the correct value. A better means of comparison of model results is through the structure or Kurtosis of the filament cross-section. However this requires high-resolution information on tracer concentration, from low-level remote sensing or fine scale in situ surveys. Although such data is currently unavailable, it should be available in the near future, due to developments in instrumentation and the increasing frequency of multi-disciplinary cruises. Subsequent to completion of this work and submission of this manuscript we were made aware of a manuscript (Neufeld et al., submitted for publication) investigating a similar problem from a more mathematical viewpoint. Neufeld et al. (submitted for publication) also predict the existence of a regime where tracer reaction (there auto-catalysis) can affect filamental width and derive an estimate for the width from a dimensional argument that is consistent with our Eq. (7).
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
5. Conclusions The structure observed in a reactive tracer field, such as plankton, is often markedly different from that of an inert tracer, (Denman and Platt, 1976; Seuront et al., 1999). We have shown that one reason may be a different functional form for the expected asymptotic width of any filaments formed. sffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k m W 8:6 max 1; 0:4 1:6 1 (7) l l For relatively large growth rate or small strain rate, such that m/l /16, the predicted structure differs significantly from that predicted for an inert or exponentially-growing tracer. Observed mesoscale filament widths are expected to be less than 1 km up to the integral length scale of the straining field. These filaments are expected to have a mean convergence time of /3 days. Numerical results suggest larger filaments may be formed but there are factors likely to disrupt filaments on this scale: grazing modifying the distribution; variation of the current field and growth rate on temporal and spatial scales less than those required for the asymptotic solution to be achieved. It should be noted that using a more sophisticated growth model, specifically explicitly modelling nitrate instead of using a logistic growth term, gives near identical results (McLeod, pers. comm.). An interesting feature of our results is that the biological processes assumed can have a marked effect on the filamental structure. This is in contrast to the case of an exponentially-growing tracer where filament width is independent of growth rate. Perhaps counter-intuitively imposing a limit to the population such that fronts of plankton propagate out from the centre of a convergent flow, resulting in a wider final distribution than if growth was unrestricted. The imminent availability of sufficiently high-resolution remote sensed and in situ data for the ecological structure of planktonic filaments will allow this hypothesis to be tested. Meanwhile explicit simulation of such tracers in two-dimensional turbulent fields and the explicit inclusion of
119
processes such as nutrient limitation and grazing provide ample scope for investigators to improve our scanty understanding of the ubiquitous plankton filament.
Acknowledgements This work was funded by a NERC studentship (GT04/99/MS/262).
References Abraham, E.R., Law, C.S., Boyd, P.W., Lavender, S.J., Bowie, M.T., 2000. Importance of stirring in the development of an iron-fertilized phytoplankton bloom. Nature 407, 727 /733. Bainbridge, R., 1957. The size shape and density of marine phytoplankton concentrations. Biological Review 32, 91 / 115. Boyd, P.W., Watson, A.J., Law, C.S., Abraham, E.R., Trull, T., Murdoch, R., et al., 2000. A mesoscale phytoplankton bloom in the polar southern ocean stimulated by iron fertilization. Nature 407, 695 /702. Damko¨hler, R., Heumann, T., 1982. Experimental studies to clarify the diffusion behavior of copper /nickel /alloys with a high copper content. Physica Status Solidi A */Applied Research 73 (1), 117 /127. Denman, K.L., Platt, T., 1976. The variance spectrum of phytoplankton in a turbulent ocean. Journal of Marine Research 34 (4), 593 /601. Fasham, M.J.R., 1978. The application of some stochastic processes to the study of plankton patchiness. In: Steele, J.H. (Ed.), Spatial Patterns in Plankton Communities. New York, Plenum Press, pp. 131 /156. Fisher, R.A., 1937. The wave advance of advantageous genes. Annals of Eugenics 7, 355 /369. Garrett, C., 1983. On the initial streakiness of a dispersing tracer in two- and three-dimensional turbulence. Dynamics of Atmospheres and Oceans 7, 265 /277. Griffa, A., Owens, K., Piterbarg, L., Rozovskii, B., 1995. Estimates of turbulence parameters from Lagrangian data using stochastic particle model. Journal of Marine Research 53, 371 /401. Haury, L.R., McGowan, J.A., Weibe, P.H., 1978. Patterns and processes in the time-space scales of plankton. In: Steele, J.H. (Ed.), Spatial Patterns in Plankton Communities. Plenum Press, New York, pp. 277 /327. Ledwell, J.R., Watson, A.J., Law, C.S., 1998. Mixing of a tracer in the pycnocline. Journal of Geophysical Research 103 (C10), 21499 /21529. Mackas, D.L., Denman, K.L., Abbott, M.R., 1985. Plankton patchiness: biology in the physical vernacular. Bulletin of Marine Science 37 (2), 652 /674.
120
P. McLeod et al. / Ecological Modelling 158 (2002) 111 /120
Martin, A.P., 2000. On filament widths in oceanic plankton distributions. Journal of Plankton Research 22 (3), 597 / 602. Murray, J.D., 1993. Biological waves: single species models. In: Mathematical Biology. New York, Springer, pp. 274 /310. Neufeld, Z., Haynes, P.H., Te´l, T., Chaotic mixing induced transitions in reaction-diffusion systems, Chaos, Submitted for publication. Okubo, A., 1971. Oceanic diffusion diagrams. Deep-Sea Research 18, 789 /802. Ottino, J.M., 1989. The Kinematics of Mixing: Stretching, Chaos and Transport. Cambridge university press. Rood, R.B., 1987. Numerical advection algorithms and their role in atmospheric transport and chemistry models. Reviews of Geophysics 25 (1), 71 /100. Schafer, H., Krauss, W., 1995. Eddy statistics in the South Atlantic as derived from drifters drouged at 100 m. Journal of Marine Research 53, 403 /431. Seuront, L., Schmitt, F., Lagadeuc, Y., Schertzer, D., Lovejoy, S., 1999. Universal multifractal analysis as a tool to
characterize multiscale intermittent patterns: examples of phytoplankton distribution in turbulent coastal waters. Journal of Plankton Research 21 (5), 877 /922. Simberloff, D.S., Wilson, E.O., 1969. Experimental zoogeography of islands: the colonization of empty islands. Ecology 50, 278 /296. Smith, C.L., Richards, K.J., Fasham, M.J.R., 1996. The impact of mesoscale eddies on plankton dynamics in the upper ocean. Deep-Sea Research I 43 (11-12), 1807 /1832. Smolarkiewicz, P.K., 1984. A fully multidimensional positive definite advection transport algorithm with small implicit diffusion. Journal of Computational Physics 54, 325 /362. Steele, J.H., 1974. Spatial heterogeneity and population stability. Nature 248, 83. Sundermeyer, M.A., Price, J.F., 1998. Lateral mixing and the North Atlantic tracer release experiment: observations and numerical simulations of Lagrangian particles and a passive tracer. Journal of Geophysical Research 103 (C10), 21481 / 21497.