- Email: [email protected]
The Influence of Density Stratification on Particle Settling, Dispersion and Population Growth
J. theor. Biol. (1997) 187, 65–75
The Influence of Density Stratification on Particle Settling, Dispersion and Population Growth S A. C* M B CSIRO Land and Water, GPO Box 821, Canberra, ACT 2601, Australia (Received on on 24 September 1996, Accepted in revised form on 29 January 1997)
Density stratification is common throughout much of our atmosphere, oceans and inland waters. It can have a strong influence on the sinking and flotation of a wide range of particulate matter, including organisms whose life cycle incorporates a water or airborne phase. Ecologically significant examples range from the dispersion of pollen, larvae and plankton, to the dispersion and deposition of sediment bound organic solids and heavy metals. While the biological effects of settling have been studied in turbulent environments, the influence of stratification has not previously been quantified. Here we develop a theory to predict the time dependent particle concentrations throughout a fluid characterised by non-uniform stratification. Factors such as dispersion, deposition and growth of particle populations have been incorporated into the model. It is shown that stratification can significantly enhance settling rates and thereby influence the survival and dispersion of colonising populations. The theory also provides simple physical explanations for commonly observed features of particle distributions, such as subsurface maxima and the onset and decay of algal blooms. 7 1997 Academic Press Limited
Introduction The movement of small particles and organisms is often determined by the motions of the fluid medium in which they are suspended. Dispersion and survival then depend on the amount of time that the particles remain suspended, which in turn is a function of the prevailing stratification. An obvious example is colonisation of habitable regions separated by hostile environments. This includes transfer of seeds and pollen by winds and transport of marine larvae across open ocean. Similar considerations apply to the transport of particle reactive chemicals. For example, the spread of organic solids or heavy metals from industrial sources, and the resulting concentrations in aquatic sediments, are dependent on both the suspension time and the biochemical reaction time-scales. The ecological significance of these and
* Present address: CSIRO Division of Marine Research, GPO Box 1538, Hobart, Tasmania 7001, Australia. 0022–5193/97/130065 + 11 $25.00/0/jt970417
many other examples, indicates a strong need to understand settling in stratified environments. The rate at which sediments or organisms sink or float through motionless air or water depends on their size, shape and density, as well as the density and viscosity of the air or water (Okubo, 1980; Hinds, 1982). However, within a vigorously turbulent flow mixing ensures that the distribution remains homogeneous. Settling behaviour then only becomes significant near horizontal boundaries where vertical fluid velocities diminish (Martin & Nokes, 1988). It has recently been demonstrated that a similar transition occurs at a diffusive interface separating two turbulent zones (Kerr & Kuiper, 1997). Here it is proposed that a similar transition also occurs between a turbulent region and a stratified region where turbulence is suppressed by buoyancy forces. While isolated patches of turbulence are often generated by velocity shear in stratified regions (Thorpe, 1987), these entrain particles from both above and below, then collapse without making any net contribution to the vertical particle transport. It therefore seems 7 1997 Academic Press Limited
. . .
reasonable to assume that transport throughout the stratified region is controlled by the settling velocity. The differing settling behaviours in turbulent and stratified regions of flow will be demonstrated to have importance implications for the dispersion, deposition and growth of aquatic and airborne species. This is achieved by isolating the dominant processes and developing simple theoretical descriptions, which are later validated with a numerical solution of the full flow equations. We begin by considering the settling of inert particles in a system with turbulent layers above and below a stratified region. The implications for horizontal dispersion and deposition are then investigated. Finally, the theory is extended to include exponential growth of particle populations and the limitations imposed by settling and competition with other suspended organisms. Assumptions and Boundary Conditions The theory makes a number of simplifying assumptions. The first is that the density difference between the particles and the fluid is much greater than differences within the fluid. This implies that the settling velocity V is independent of depth. The second is that motions in the turbulent regions are sufficiently energetic to keep the particle distribution uniform, and conversely that turbulent energy in stratified layers is sufficiently weak that settling motions dominate the particle behaviour. This can be expressed more precisely in terms of the Peclet number, defined by Pe = Vh/k, where h is the layer depth and k is the eddy diffusivity of the particles. The turbulent regions are then assumed to be characterised by small Peclet numbers and stratified regions by large Peclet numbers. The applicability of this assumption is tested in Appendix A by comparing the theoretical results with those from a numerical model, which makes no assumptions about the flow’s Peclet number. The third assumption is that the time-scale over which the stratification changes is long compared with the settling time H/V, where H is the total depth of the system. This also requires that the particles do not contribute significantly to the fluid density (Kerr, 1991). However, this assumption can be relaxed by applying the theory repeatedly over short intervals in which the change in stratification is not significant. This was the basis for the iterative approach adopted in the surface mixed layer model of Sherman & Webster (1994). The three assumptions allow simple analytical solutions to be found which can be applied
0 hS
z/H
66
hT
hB 1 Density F. 1. Schematic representation of the three layer density structure considered in the theory.
to quite arbitrary stratification. While not valid in all situations, they should find wide applicability in natural flows. In order to incorporate a wide range of situations and conditions, we consider a three layer system consisting of a well mixed upper layer, a central stratified zone or thermocline, and a well mixed bottom layer (Fig. 1). Such conditions are very common in water bodies where the surface layer is mixed by wind or convection, the thermocline is generated and maintained by solar heating or horizontal advection, and the turbulent bottom layer is maintained by bottom stress opposing the mean flow. Analogous conditions are also found in the atmospheric boundary layer where the stratified region takes the form of an inversion layer. The analysis for the three layer system can also be easily extended to an arbitrary number of stratified layers separated by turbulent regions. The depths of the surface layer, thermocline and bottom layer will be denoted by hS , hT and hB , respectively, while their corresponding particle concentrations are CS (t), CT (z, t) and CB (t). Note that the suppression of turbulence in the thermocline allows vertical (z) gradients in both density and particle concentration to be maintained. Settling will usually be directed towards the bottom and solutions will be obtained with and without deposition. For inert particles, the one-dimensional system is vertically symmetric and solutions are the same for sinking and floating particles. However, if particle population growth is concentrated in one layer, then symmetry is lost and sinking and floating populations need to be considered separately. The initial conditions are specified in terms of a surface layer concentration C0 and a uniform concentration C* in the thermocline and bottom layer:
CS (0) = C0 , CT (z, 0) = CB (0) = C* .
(1)
The examples will use fixed layer depths of hS = 0.40H, hT = 0.35H and hB = 0.25H, and will focus on two initial conditions of particular interest. The first is an initially uniform concentration throughout the flow (C* = C0 ), encountered when a mixing event is rapidly followed by restratification. The second is an initial population concentrated entirely in the first layer (C* = 0), encountered after a sudden influx or population growth within that layer. Settling Behaviour of Inert Particles Since particles escape from the first layer at their settling speed V, the concentration within the surface layer will decay according to dCS V = − CS , dt hS
(2)
The concentration within the bottom layer CB is determined by the flux from above and settling losses below: dCB FT − VCB , = dt hB
(3)
0 1 tQ
C B = C*
=C0
6
(4)
Turbulence is much weaker within the thermocline and particle motions are dominated by settling behaviour. The particle concentration CT is therefore described by the advection equation 1CT 1C = −V T , 1t 1z
(5)
which when combined with flux continuity at z = hS has the solution CT = C0 e−(Vt − z + hS )/hS =C* ,
(6)
where z is positive downwards. The corresponding particle flux into the bottom layer is FT = VC*
0 1 h tQ T V
=VC0 e−(Vt − hT )/hS
70 1
h* e−(Vt − hT )/hB−hS e−(Vt − hT )/hS , hB−hS
DC = C0
tq
hT V
(9)
6
7
h* e−(Vt − hT )/hB − hB e−(Vt − hT )/hS . hB − hS
(10)
The discontinuity reaches its maximum value at time tD =
6
0 17
hB2 1 hS hB hT − ln V (hS − hB ) h* hS
,
(11)
before decaying towards the limit DC : CB . For deep bottom layers (hB hS ), h* = hS + hB C* /C0 and tD =
6
0 17
hB2 1 hT + hS ln h* hS V
,
(12)
which continues to increase with hB . The mean concentration over the entire water column can also be determined from eqns (3), (6) and (9),
(hS Q z Q hS + Vt)
(hS + Vt Q z Q hS + hT )
hT V
where h* = hS + (hB − hS )C* /C0 . Because incoming particles are immediately mixed through the bottom layer, a discontinuity in concentration can develop at the base of the thermocline even when the initial distribution is continuous. This is first evident at t = hT /V, from which its magnitude develops as
while the corresponding flux of particles into the thermocline is FS = VC0 e−Vt/hS .
(8)
which combined with (7) has the solution
provided the Peclet number is small (Smith, 1982; Martin & Nokes, 1988). This can be integrated to give CS = C0 e−Vt/hS ,
67
C =
C0 h* + C* (hS + hT − Vt) H
=C0
6
0 1 tQ
hT V
70 1
h* hB e−(Vt − hT )/hB − hS2 e−(Vt − hT )/hS H(hB − hS )
tq
hT . V (13)
0 1
h tq T . V
(7)
The mean concentration therefore falls linearly for a period determined by the thermocline thickness, followed by a transition to an exponential decay.
. . .
68 0
(a) t=1.1H/V t=0.5H/V
z/H
0.2
t=0.1H/V
0.4 0.6 0.8 1.0 0
0.2
0.4
0.6
0.8
1.0
C/C o 0
(b) t=1.1H/V t=0.5H/V
z/H
0.2
t=0.1H/V
0.4 0.6 0.8
begins to decay as particles first enter the bottom mixed layer (t q hT /V). It finally disappears at time tD , after which the behaviour becomes qualitatively similar to that for uniform initial concentration. However, when the bottom layer is deep (hB hS ), eqn (12) indicates that the concentration maxima can persist for long periods and therefore should be a relatively common feature in natural systems which satisfy this criteria. The size of the discontinuity calculated from (10) is shown by the dashed lines in Figs 3(a) and 3(b). In both cases, the discontinuity first appears at t = hT /V when particles from the surface layer first reach the bottom mixed layer. For an initially uniform distribution, concentrations in the thermocline fall more rapidly than those in the bottom layer, which results in an increase in DC. However, as particles continue to be lost from the bottom layer, the discontinuity begins to decay towards the limit DC : CB . The size of the DC peak increases with the
1.0 0
0.2
0.4
0.6
0.8
1.0
C/C o F. 2. Vertical profiles of particle concentration at non-dimensional times tV/H = 0.1, 0.5 and 1.1, for stratification defined by hS /H = 0.40, hT /H = 0.35 and hB H = 0.25. (a) Initially uniform particle concentration (C* = C0 ). (b) All particles initially concentrated in the surface layer (C* = 0).
(a) C T/C o
0.6
C B/C o
C/C o
C S /C o 0.2
∆C/C o
–0.2 –0.6 –1.0 0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
tV/H 1.0
(b) C T/C o
0.6
C/C o
Examples of vertical profiles calculated for the two initial conditions are shown in Fig. 2, with layer depths specified in the caption. Corresponding time series for concentrations within the two mixed layers and mid-depth in the thermocline are shown in Fig. 3. Concentrations are uniform in the well-mixed surface layer and then increase exponentially with depth through the thermocline. The distributions are identical for the two initial conditions to a depth of hS + Vt, as required by (6). For uniform initial concentration [Figs 2(a) and 3(a)], fluxes through the base of the thermocline initially balance losses to the bottom boundary layer, so that the bottom layer concentration is maintained at C0 . However, for t q hT /V a discontinuity develops at the base of the thermocline where particles are suddenly mixed through the bottom layer. Very similar results are obtained from numerical solutions of the primitive flow equations incorporating turbulent transport of settling particles (Appendix A). The development is quite different when particles are initially concentrated in the surface layer [Figs 2(b) and 3(b)]. Of particular significance is a sharp concentration maxima in the thermocline peaking at z = hS + Vt for t Q hT /V. The feature
1.0
C B/C o
0.2
C S /C o
∆C/C o
–0.2 –0.6
–1.0 0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
tV/H F. 3. Time series of particle concentrations in the surface mixed layer, bottom mixed layer, and midpoint in the thermocline. The size of the concentration discontinuity at the base of the thermocline is shown by the dashed curve. The stratification is the same as that in Fig. 2. (a) Initially uniform particle concentration (C* = C0 ). (b) All particles initially concentrated in the surface layer (C* = 0).
69
1.0
Population Dispersion: Horizontal Advection and Downstream Deposition
0.8
If populations are generated or suspended within a localised region, the horizontal flow will tend to deposit particles with some characteristic downstream distribution. Consider for example, a horizontal flow of uniform velocity u in the x direction (neglecting boundary layer shear), with populations defined by (1) at a position x = 0. The deposition flux is then given simply by FB = VCB , which after a transformation to spatial coordinates (t : x/u) in eqn (9) gives
C/C o
0.6 Fully mixed limit
0.4
0.2
0
C/C o
Stratified limit 0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
0
FB = VC*
tv/H F. 4. Time series of depth averaged particle concentration for the stratification in Fig. 2, along with the fully turbulent (hS = H) and fully stratified (hT = H) limits.
=VC0
xQ
uhT V
1
6
70
h* e−(Vx/u − hT )/hB − hS e−(Vx/u − hT )/hS hB − hS
xq
1
uhT . V (14)
For a particle cloud of finite width Dx, the final distribution of deposited particles is FB Dx/u. If the non-dimensional deposition flux, FB /VC0 , is plotted against the non-dimensional downstream distance, xV/uH, the curves are identical to those for CB /C0 in Fig. 3. For a vertically uniform upstream distribution, Fig. 3(a) suggests uniform deposition for a distance uhT /V, followed by an exponential decay. However, if particles are initially deposited within the upper layer [Fig. 3(b)], then there is no deposition over the
1.0
0.8 1.0
1.0
0.6 h T/H
ratio hB /hS , with the limits DC : 0 as hB /hS : 0 and DC : C0 as hB /hS : a. Stratification not only affects the concentration distribution, but also the depth averaged concentration. This quantity has been calculated from (13) and plotted as a function of time in Fig. 4. The corresponding trends for fully mixed and fully stratified conditions are also shown for comparison. The presence of stratification clearly produces a more rapid decline in concentration levels, which strongly impacts long-term population levels (t q H/V). In this example, there is almost an order of magnitude difference in the populations of the partially stratified and fully turbulent flows at t = 2H/V. The depth averaged concentrations in Fig. 4 can be generalised to other stratifications by using (13) to calculate the time, t90 , required for 90% of the initial concentration to be deposited. Vt90 /H is contoured in Fig. 5 as a function of the surface layer and thermocline depths. Particle suspension times decrease as the thickness of the thermocline expands, with Vt90 /H falling by a factor of 2.5 from fully mixed to fully stratified conditions. This factor grows with deposited fraction, so that the time required to deposit 99% of the initial concentration in a fully turbulent environment is 4.6 times longer than in a fully stratified environment. Figure 5 further indicates that the presence of even a thin thermocline results in a large reduction in the suspension time. This is because the e-folding time for particle decay in each turbulent mixed layer is proportional to the depth of that layer.
0.4 1.4 1.3 1.2
0.2
1.7 1.6 1.8 2.0 1.9
1.2 1.3 1.4
1.4 1.5
1.6 1.7 1.9 1.8 2.0
0 0
0.2
0.4
0.6
0.8
1.0
h S/H F. 5. Contour plot of non-dimensional time (tV/H) required for 90% of the particles initially present to settle. This is shown as a function of surface mixed layer depth and thermocline depth.
. . .
70
distance uhT /V, following by a relatively rapid increase to a peak deposition at the location utD (C* = 0).
(a) 3 2.0
1.5
2
The theory can now be extended to organisms characterised by a growth rate m, which encompasses all physiological and environmental factors except the tendency to settle through their fluid environment. Only species in which settling motions dominate over swimming or buoyancy regulation are considered.
C/C o
Population Growth
1.1 1 0.9 G=0.0 0.1
0.5
0 0
0.4
0.8
1.2
0
1
1.6
2.0
2.4
2.0
2.4
2.0 1.5 1.1 0.9
2
(15)
G=0.0
0 0
0.4
Many plankton species have a tendency to float towards higher light levels and in some cases settle along the surface as a scum. The theory can be applied directly to this situation by inverting the model so that the layer denoted by subscript B becomes the surface layer with growth. This approach neglects growth within the thermocline, so that the populations below the surface layer are identical to
0.8
1.2 tV/H
(c) 3
2.0 1.5
(16)
where G = m¯ hS /V is a measure of the relative importance of growth and settling. Clearly, populations grow for G q 1 and decay for G Q 1, as shown in Fig. 6(a). If growth below the surface layer is negligible due to light limitation, then the development of the population at these levels is equivalent to that for inert particles with a reduced settling speed (1 − G)V. These results can be applied to time dependent systems by calculating solutions incrementally while continuously updating m¯ , hS and C0 (Sherman & Webster, 1994).
0.5
0.1
1
2
1.1
1
0.9
C/C o
C S = C0 e
2.4
(b) 3
hS where m¯ = h−1 S f0 m dz is the mean growth rate for the layer (m¯ is calculated in Appendix B for the case of light limited growth in aquatic environments). Over time-scales for which m¯ and hS are approximately constant, this has the solution −(1 − G)Vt/hS
2.0
C/C o
Consider the case where surface layer populations grow exponentially and are only limited by their tendency to sink. This balance often controls phytoplankton populations in species with negative buoyancy. The surface layer population is then described by V dCS C , = m¯ − hS S dt
1.6
tV/H
G=0.0
0.1
0.5
0 0
0.4
0.8
1.2
1.6
tV/H F. 6. Time series of particle populations with growth defined by the parameter G = m¯ hS /V and stratification as in Fig 2. (a) Particles sinking out of the surface layer. (b) Particles initially uniformly distributed (C* = C0 ), floating through the surface layer and forming a scum. (c) Particles initially concentrated within the bottom layer, floating into the surface layer and forming a scum.
those of inert particles. The surface layer can be described by dCB FT − (V − m¯ hB )CB = , dt hB
(17)
71
which combined with (6) has the solution
CB = C*
6
70 1
1 − G e−(1 − G)Vt/hB 1−G
tQ
hT V
[C* hB − (C* − C0 )(1 − G)hS ] e−((1 − G)Vt − hT (1 − G))/hB /(1 − G) − C0 hS e−(Vt − hT )/hS CG − * e−(1 − G)Vt/hB hB − hS (1 − G) 1−G
=
0 1 tq
hT V
(18) where G = m¯ hB /V. In the limit of zero growth (G = 0), eqn (18) reduces to (9). Figure 6(b) shows the growth curves for the case C* = C0 and a range of G values. In all cases there is an initial increase in population at the natural growth rate m¯ . This is because losses to the surface are initially exactly compensated by an influx from the thermocline. However, as this flux diminishes in accordance with (7), growth is only maintained for G q 1, with decaying populations for G Q 1. In some circumstances, particles reaching the surface are stirred back into the layer rather than forming a scum. This usually applies to phytoplankton under relatively strong wind conditions (Hutchinson & Webster, 1994). The settling term in (17) can then be omitted and the solution takes the form CB = C* 4(1 + G−1 ) e GVt/hB − G−15 −1
0 1 tQ
hT V
−1
=(C* (1 + G ) − [C* G − C0 (G + hB /hS )−1 ] e−GhT /hB ) e GVt/hB − C0 (G + hB /hS )−1 e−(Vt − hT )/hS
0 1 tq
hT . (19) V
Since there are no loss terms, this solution describes monotonic exponential growth for all values of G. The growth rates are not shown, but are very similar to those for the scum forming particles characterised by a higher growth factor G + 1 in Fig. 6(b). Seeding populations can also rise from depth, after being resuspended from the bottom. For example, some cyanobacteria produce resting spores (akinetes) or have overwintering benthic phases. The inverted model using (18) or (19) with C* = 0 can be easily applied to this situation, with C0 representing the seed population. Figure 6(c) shows the growth curves for a range of G values. Particles first reach the surface layer at t = hT /V and then begin to multiply. However, with settling into a surface scum, populations are again only sustained for G q 1. Without scum formation (not shown), populations grow exponentially for all positive values of G.
/ The manner in which two plankton populations compete for available chemical resources (nutrients and inorganic carbon) is highly species dependent. However, a qualitative understanding can be gained from a simple Lotka–Volterra model. Once any fluxes from the thermocline have declined, the population of each species can be described by (15), with each m¯ a function of its own concentration (intraspecific competition) and the concentration of other species (interspecific competition). The simplest case is a linear dependence of growth rates on the population of two species (denoted by subscripts 1 and 2). The stability of such a system is well known (MaynardSmith, 1974) and can be expressed in terms of the parameter R=
m¯ 1 hS − V1 V1 (G1 − 1) = , m¯ 2 hS − V2 V2 (G2 − 1)
(20)
where G1 = m¯ 1 hS /V1 and G2 = m¯ 2 hS /V2 . Only species 2 will persist in the surface layer if G2 q 1 and R Q (1m¯ 1 /1C2 )/(1m¯ 2 /1C2 ); the two species will coexist if G1 q 1, G2 q 1 and (1m¯ 1 /1C2 )/(1m¯ 2 /1C2 ) Q R Q (1m¯ 1 / 1C1 )/(1m¯ 2 /1C1 ); while only species 1 will persist if G1 q 1 and R q (1m¯ 1 /1C1 )/(1m¯ 2 /1C1 ). While the growth derivatives can only be determined empirically, the analysis at least demonstrates that R is the critical parameter controlling competition in this type of system.
Implications Large peaks in concentration of suspended material are often observed to coincide with zones of strong density gradient. This is sometimes associated with favourable light or nutrient levels, such as in the case of phytoplankton maxima near the base of the ocean thermoclines (Sharples & Tett, 1994; Denman &
72
. . .
Gargett, 1995). However, similar observations have also been made of non-reproductive particles, such as larvae and marine snow (aggregates of inorganic particles, microorganisms and detritus), which are passively advected by the flow (Alldredge & Cohen, 1987; MacIntyre et al., 1995; Tremblay et al., 1994). A number of explanations have been suggested for particle accumulation in stratified regions. These include horizontal intrusions of particle laden water, reduction of settling speeds as higher densities are encountered, and increased particle aggregation in turbulent zones. However, the density conditions and particle characteristics are often not conducive to such mechanisms. Furthermore, MacIntyre et al. (1995) have used a numerical random walk model, as described by Patterson (1991), to demonstrate that accumulations occur even in the absence of these processes. The theory presented here reveals the underlying mechanism responsible for the behaviour of the random walk model. Specifically, the concentration within the overlying turbulent layer decays uniformly, while the peak in the pycnocline moves downward at the particle settling speed. There is no decrease in the magnitude of the peak until it reaches the turbulent layer below. Furthermore, if this turbulent layer is much deeper than the surface layer, then the pycnocline peak will be a very persistent feature of the concentration profile [eqn (12)]. This may explain why such features are so frequently observed in oceanic systems. The theoretical results have particularly significant implications for the survival and dispersion of species which rely on suspended seeding populations to migrate between isolated patches of habitat or recolonise an environment each season. This is due to the strong dependence of population decay rates on stratification. If the water or air column is turbulent, then small numbers of particles can remain suspended for extended periods. This may allow micro-organisms to overwinter in suspension, or assist seeds or larvae in reaching more distance habitats. However, stratification tends to constrain this mechanism and may be responsible for anomalously poor seasons or reduced colonisation ranges. Figure 6 demonstrates the strong dependence of growing populations on the parameter G for quite a diverse range of circumstances. A particular example to which this may be relevant is the relative growth rates of diatoms and cyanobacteria (blue-green algae)
within Australia’s Murray–Darling river system. The numerical model described in Appendix A has been successfully used to make detailed predictions of these populations, including the formation of blooms (Bormans & Webster, 1996). However, the discussion here will focus on testing the underlying conditions required by the theory for one of the plankton species to proliferate. In Fig. 7, a four month time series of diatom populations measured in the lower Murrumbidgee River is compared with the associated time series of the growth factor G. Initially, both G and the diatom population are relatively high. However, a period of low flow extending into December allows the development of stratification, with a corresponding sharp reduction in G. Settling then dominates over growth and the population falls. When G again increases with the flow rate, the population recovers rapidly. While generally remaining high until the next sustained reduction in G around mid-February, there are small population drops associated with shorter term (H/V 1 5 days) reductions in G. Daily oscillations in G (associated with cloudy days) do not have a significant effect, because plankton settling through the thermocline are re-entrained into the deepening mixed layer before they can reach the bottom. The strong correlation between population changes and G (over time-scales greater than H/V) appears to support the theoretical framework and suggests that G might be a useful indicator of bloom formation in other populations limited by their settling behaviour.
Conclusion A simple theory has been developed to describe particle settling in stratified systems. It is based on the premise that settling in turbulent regions obeys a well known exponential decay law based on the turbulent layer depth, while settling in stratified regions is at the particle settling speed. This can lead to a range of interesting phenomena, such as discontinuities and sharp peaks in concentration profiles. The inclusion of particle population growth in the turbulent surface layer reveals that G = m¯ hS /V is the critical parameter controlling growth of populations limited by settling. There remains a large range of related phenomena to be investigated. These include other initial conditions (selected for specific applications), time dependent stratification, z dependent settling speeds (including inhibition by large particle concentrations), and growth within stratified regions. It is hoped that this paper will provide a theoretical basis for studying
73 10 000
1.6
1.4
G
1.0 100
0.8 0.6
Diatoms (cells ml –1)
1000
1.2
10
0.4 0.2
1 06/03/95
26/02/95
18/02/95
10/02/95
02/02/95
25/01/95
17/01/95
09/01/95
01/01/95
24/12/94
16/12/94
08/12/94
30/11/94
22/11/94
14/11/94
0
F. 7. Depth averaged concentrations of the diatom Aulacoseira sp. determined from cell counts collected every three days in the lower Murrumbidgee River in southeastern Australia between 14 November 1994 and 8 March 1995 (w). This is plotted alongside the growth factor G (W), which was calculated at two hourly intervals from (B.3) and then diurnally averaged. The calculation used mixed layer depths based on two hourly averaged temperature profiles from thermistor chain measurements, along with locally measured I0 (t) and h. Details of the field measurements are given by Webster et al. (1996). From laboratory and field data, the diatom characteristics were estimated to be V = 1.0 m day−1, mmax = 0.7 day−1 and Ik = 10 mEin m−2 s−1 (Reynolds, 1984). While there are significant uncertainties in these estimates, the trends revealed by the plot are not sensitive to their specific values.
these issues within the context of more specific applications. We thank Dr Gary Jones, CSIRO Land and Water, for the diatom count data from an NRMS funded project examining cyanobacteria blooms in the lower Murrumbidgee River. REFERENCES A, A. L. & C, Y. (1987). Can microscale patches persist in the sea? Microelectrode study of marine snow, fecal pellets. Science, 235, 689–691. B, A. F. & M, G. L. (1983). Diagnostic and prognostic numerical circulation studies of the South Atlantic Bight. J. Geophys. Res. 88, 4579–4592. B, M. & W, I. T. (1996). Modelling the dynamics of the diatom Aulacoseira sp. and the blue-green Anabaena in relation to river flow. CSIRO Centre for Environmental Mechanics Technical Report No. 127.. D, K. L. & G, A. E. (1995). Biological-physical interactions in the upper ocean: The role of vertical and small scale transport processes. Ann. Rev. Fluid Mech. 27, 225–255. H, W. C. (1982). Aerosol Technology: Properties, Behaviour, and Measurement of Airborne Particles. New York: Wiley. H, P. A. & W, I. T. (1994). On the distribution of
blue-green algae in lakes: Wind-tunnel tank experiments. Limnol. Oceanogr. 39, 374–382. K, R. C. (1991) Erosion of a stable density gradient by sedimentation-driven convection. Nature 353, 423–425. K, R. C. & K, G. S. (1997). Particle settling through a diffusive-type thermocline staircase in the ocean. Deep Sea Res. 44, 399–412. MI, S., A, A. L. & G, C. C. (1995). Accumulation of marine snow at density discontinuities in the water column. Limnol. Oceanogr. 40, 449–468. M, D. & N, R. (1988). Crystal settling in a vigorously convecting magma chamber. Nature 332, 534–536. M-S, J. (1974). Models in Ecology. Cambridge: Cambridge University Press. M, G. L. & Y, T. (1974). A hierarchy of turbulence closure models for planetary boundary layers. J. Atmos. Sci. 31, 1791–1806. M, G. L. & Y, T. (1982). Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. Space Phys. 20, 851–879. O, A. (1980). Diffusion and Ecological Problems: Mathematical Models. Berlin: Springer-Verlag. P, J. C. (1991). Modelling the effects of motion on primary production in the mixed layer of lakes. Aquat. Sci. 53, 218–238. R, C. S. (1984). The Ecology of Freshwater Phytoplankton. Cambridge: Cambridge University Press. S, J. & T, P. (1994). Modelling the effect of physical variability on the midwater chlorophyll maximum. J. Mar. Res. 52, 219–238.
. . .
74
S, V. S. & W, I. T. (1994). A model for the light-limited growth of buoyant phytoplankton in a shallow, turbid waterbody. Aust. J. Mar. Freshwater Res. 45, 847–862. S, I. R. (1982). A simple theory of algal deposition. Freshwater Biol. 12, 445–449. T, S. A. (1987). Diapyanal mixing in the thermocline: A review. J. Geophys. Res. 92, 5231–5248. T, M. J., L, J. W., W, F. E., N, C. E., P, F. H. & S M. M. (1994). Drift of sea scallop larvae Placopecten magellanicus on Georges Bank: a model study of the roles of mean advection, larval behavior and larval origin. Deep Sea Res. 41, 7–49. W, I. T., J, G. J., O, R. L., B, M. & S, B. S. (1996). Control strategies for cyanobacterial blooms in weir pools. CSIRO Centre for Environmental Mechanics Technical Report No. 119.
Figure A1(b) shows the intensity of the turbulence throughout the water column in terms of the Peclet number, Pe = VH/k, based on the total water depth. This will help interpret comparisons of concentration profiles from the theoretical and numerical solutions [Fig. A1(c)]. In the surface mixed layer Pe1, so that concentrations are always homogeneous with good agreement between theoretical and numerical esti
0
(a) 1.20 0.60 0.24
0.2
0.12
An important assumption of the theory is that the turbulent and stratified regions are respectively characterised by small and large Peclet numbers. To investigate this issue within the context of more realistic flows, the theoretical results have been compared with those from a numerical model. This has been achieved by incorporating settling particles into the primitive equation circulation model of Blumberg & Mellor (1983). This solves the nonlinear hydrostatic form of the flow equations using standard finite difference methods. Characteristics of the subgrid scale turbulence were determined using the turbulence closure scheme developed by Mellor & Yamada (1974, 1982). Flow fields and particle distributions have been calculated from this model for a small idealised lake, forced by wind stress and heat flux through the free surface. The numerical runs were conducted in a basin of depth H = 3.0 m and fetch 100 m, while neglecting gradients in the cross-wind direction. The primitive equations were solved on a uniform grid with 0.1 m vertical resolution and 0.5 m horizontal resolution. The initial stratification consisted of an upper mixed layer of depth hS = 1.14 m and temperature 30°C, a thermocline of depth hT = 1.11 m and linear temperature gradient, and a bottom mixed layer of depth hB = 0.75 m and temperature 20°C [Fig. A1(a)]. There was initially a uniform distribution of particles, with uniform settling velocity V = 1 × 10−5 ms−1. The flow was forced by a constant windstress and surface heat flux. These were tuned so that the depths of the three layers remained approximately constant over the duration of the run [Fig. A1(a)], while also maintaining turbulent flow within the two mixed layers. This system could be compared directly with theoretical results for inert particles.
z/H 0.6 0.8 1.0 15
20
25
30
35
Temperature (°C) 0
(b)
0.2 0.12
0.4
0.24
z/H
Comparison With a Numerical Solution
0.4
0.60 1.20
0.6 0.8 1.0 1.0e–02
1.0e+00
1.0e+02
1.0e+04
Pe 0 0.2
(c) 1.20
0.60
0
0.2
0.24
0.12
0.4 z/H
APPENDIX A
0.6 0.8 1.0 0.4
0.6
0.8
1.0
C/C o F. A1. (a) Development of the stratification in the numerical model with non-dimensional time, tV/H. (b) Development of the rms turbulent velocity profile with non-dimensional time. (c) Comparison of the numerical and theoretical concentration profiles at four non-dimensional times.
APPENDIX B In many instances, population growth in aquatic environments is limited by the available light. This can be described in terms of an instantaneous growth rate mmax I , I + Ik
0.8 100 0.6 10 0.4 1.0 0.2 0 0 0
2
(B.1)
where mmax is the maximum growth rate, Ik is the half-saturation constant for light limited growth, and I is the photosynthetically active radiation (PAR) at depth z. Provided the particles do not contribute significantly to light attenuation, the PAR can usually be approximated by Beer’s Law, I = I0 e−hz, where I0
4
6
8
10
ηh S F. B1. Average surface layer growth rates as a function of layer depth (non-dimensionalised by the light extinction coefficient) for a range of I0 /Ik values.
is the value at z = 0, and h is the attenuation coefficient. The mean growth rate is then m¯ =
Light Limited Growth
m=
1.0
G/G max
mates. Within the thermocline Pe1, and similar concentration gradients are predicted by the two models. There is some tendency for particles to settle more rapidly in the numerical case as the high Peclet number region penetrates into the bottom mixed layer. In any case, the range of bottom layer Peclet numbers (10−2 Q Pe Q 1) encompasses settling behaviour which is significantly more rapid than the small Pe limit. This leads to quite significant discrepancies between theoretical and numerical results within this layer. In most respects, the theory provides a very good description of the settling behaviour in the numerical model. However, it does provide a cautionary reminder that a well-mixed density field does not always imply strong turbulence. Since it is usually not practical to directly measure turbulent velocities, small Peclet numbers should only be assumed where there is an obvious source of turbulent kinetic energy, such as a strongly sheared boundary layer. While this was not present in the model lake, other systems such as rivers, tidal flows and atmospheric boundary layers are much more likely to satisfy this requirement.
75
1 hS
g
hS
m dz =
0
0
1
mmax I0 + Ik ln . hhS I0 e−hhS + Ik
(B.2)
The critical parameter G = m¯ hS /V determining the population behaviour can then be expressed in the form m I0 + Ik G max ln . (B.3) hV I0 e−hhS + Ik
0
1
Figure B1 is a graphic representation of (B.3) indicating that growth is favoured by strong incident radiation (large I0 /Ik ) with low attenuation through the surface layer (small hhS ). In the weak radiation limit (I0 /Ik : 0) G : 0, while in the strong radiation limit (I0 /Ik : a) G : Gmax where Gmax = mmax hS /V is independent of I0 .
The Influence of Density Stratification on Particle Settling, Dispersion and Population Growth S A. C* M B CSIRO Land and Water, GPO Box 821, Canberra, ACT 2601, Australia (Received on on 24 September 1996, Accepted in revised form on 29 January 1997)
Density stratification is common throughout much of our atmosphere, oceans and inland waters. It can have a strong influence on the sinking and flotation of a wide range of particulate matter, including organisms whose life cycle incorporates a water or airborne phase. Ecologically significant examples range from the dispersion of pollen, larvae and plankton, to the dispersion and deposition of sediment bound organic solids and heavy metals. While the biological effects of settling have been studied in turbulent environments, the influence of stratification has not previously been quantified. Here we develop a theory to predict the time dependent particle concentrations throughout a fluid characterised by non-uniform stratification. Factors such as dispersion, deposition and growth of particle populations have been incorporated into the model. It is shown that stratification can significantly enhance settling rates and thereby influence the survival and dispersion of colonising populations. The theory also provides simple physical explanations for commonly observed features of particle distributions, such as subsurface maxima and the onset and decay of algal blooms. 7 1997 Academic Press Limited
Introduction The movement of small particles and organisms is often determined by the motions of the fluid medium in which they are suspended. Dispersion and survival then depend on the amount of time that the particles remain suspended, which in turn is a function of the prevailing stratification. An obvious example is colonisation of habitable regions separated by hostile environments. This includes transfer of seeds and pollen by winds and transport of marine larvae across open ocean. Similar considerations apply to the transport of particle reactive chemicals. For example, the spread of organic solids or heavy metals from industrial sources, and the resulting concentrations in aquatic sediments, are dependent on both the suspension time and the biochemical reaction time-scales. The ecological significance of these and
* Present address: CSIRO Division of Marine Research, GPO Box 1538, Hobart, Tasmania 7001, Australia. 0022–5193/97/130065 + 11 $25.00/0/jt970417
many other examples, indicates a strong need to understand settling in stratified environments. The rate at which sediments or organisms sink or float through motionless air or water depends on their size, shape and density, as well as the density and viscosity of the air or water (Okubo, 1980; Hinds, 1982). However, within a vigorously turbulent flow mixing ensures that the distribution remains homogeneous. Settling behaviour then only becomes significant near horizontal boundaries where vertical fluid velocities diminish (Martin & Nokes, 1988). It has recently been demonstrated that a similar transition occurs at a diffusive interface separating two turbulent zones (Kerr & Kuiper, 1997). Here it is proposed that a similar transition also occurs between a turbulent region and a stratified region where turbulence is suppressed by buoyancy forces. While isolated patches of turbulence are often generated by velocity shear in stratified regions (Thorpe, 1987), these entrain particles from both above and below, then collapse without making any net contribution to the vertical particle transport. It therefore seems 7 1997 Academic Press Limited
. . .
reasonable to assume that transport throughout the stratified region is controlled by the settling velocity. The differing settling behaviours in turbulent and stratified regions of flow will be demonstrated to have importance implications for the dispersion, deposition and growth of aquatic and airborne species. This is achieved by isolating the dominant processes and developing simple theoretical descriptions, which are later validated with a numerical solution of the full flow equations. We begin by considering the settling of inert particles in a system with turbulent layers above and below a stratified region. The implications for horizontal dispersion and deposition are then investigated. Finally, the theory is extended to include exponential growth of particle populations and the limitations imposed by settling and competition with other suspended organisms. Assumptions and Boundary Conditions The theory makes a number of simplifying assumptions. The first is that the density difference between the particles and the fluid is much greater than differences within the fluid. This implies that the settling velocity V is independent of depth. The second is that motions in the turbulent regions are sufficiently energetic to keep the particle distribution uniform, and conversely that turbulent energy in stratified layers is sufficiently weak that settling motions dominate the particle behaviour. This can be expressed more precisely in terms of the Peclet number, defined by Pe = Vh/k, where h is the layer depth and k is the eddy diffusivity of the particles. The turbulent regions are then assumed to be characterised by small Peclet numbers and stratified regions by large Peclet numbers. The applicability of this assumption is tested in Appendix A by comparing the theoretical results with those from a numerical model, which makes no assumptions about the flow’s Peclet number. The third assumption is that the time-scale over which the stratification changes is long compared with the settling time H/V, where H is the total depth of the system. This also requires that the particles do not contribute significantly to the fluid density (Kerr, 1991). However, this assumption can be relaxed by applying the theory repeatedly over short intervals in which the change in stratification is not significant. This was the basis for the iterative approach adopted in the surface mixed layer model of Sherman & Webster (1994). The three assumptions allow simple analytical solutions to be found which can be applied
0 hS
z/H
66
hT
hB 1 Density F. 1. Schematic representation of the three layer density structure considered in the theory.
to quite arbitrary stratification. While not valid in all situations, they should find wide applicability in natural flows. In order to incorporate a wide range of situations and conditions, we consider a three layer system consisting of a well mixed upper layer, a central stratified zone or thermocline, and a well mixed bottom layer (Fig. 1). Such conditions are very common in water bodies where the surface layer is mixed by wind or convection, the thermocline is generated and maintained by solar heating or horizontal advection, and the turbulent bottom layer is maintained by bottom stress opposing the mean flow. Analogous conditions are also found in the atmospheric boundary layer where the stratified region takes the form of an inversion layer. The analysis for the three layer system can also be easily extended to an arbitrary number of stratified layers separated by turbulent regions. The depths of the surface layer, thermocline and bottom layer will be denoted by hS , hT and hB , respectively, while their corresponding particle concentrations are CS (t), CT (z, t) and CB (t). Note that the suppression of turbulence in the thermocline allows vertical (z) gradients in both density and particle concentration to be maintained. Settling will usually be directed towards the bottom and solutions will be obtained with and without deposition. For inert particles, the one-dimensional system is vertically symmetric and solutions are the same for sinking and floating particles. However, if particle population growth is concentrated in one layer, then symmetry is lost and sinking and floating populations need to be considered separately. The initial conditions are specified in terms of a surface layer concentration C0 and a uniform concentration C* in the thermocline and bottom layer:
CS (0) = C0 , CT (z, 0) = CB (0) = C* .
(1)
The examples will use fixed layer depths of hS = 0.40H, hT = 0.35H and hB = 0.25H, and will focus on two initial conditions of particular interest. The first is an initially uniform concentration throughout the flow (C* = C0 ), encountered when a mixing event is rapidly followed by restratification. The second is an initial population concentrated entirely in the first layer (C* = 0), encountered after a sudden influx or population growth within that layer. Settling Behaviour of Inert Particles Since particles escape from the first layer at their settling speed V, the concentration within the surface layer will decay according to dCS V = − CS , dt hS
(2)
The concentration within the bottom layer CB is determined by the flux from above and settling losses below: dCB FT − VCB , = dt hB
(3)
0 1 tQ
C B = C*
=C0
6
(4)
Turbulence is much weaker within the thermocline and particle motions are dominated by settling behaviour. The particle concentration CT is therefore described by the advection equation 1CT 1C = −V T , 1t 1z
(5)
which when combined with flux continuity at z = hS has the solution CT = C0 e−(Vt − z + hS )/hS =C* ,
(6)
where z is positive downwards. The corresponding particle flux into the bottom layer is FT = VC*
0 1 h tQ T V
=VC0 e−(Vt − hT )/hS
70 1
h* e−(Vt − hT )/hB−hS e−(Vt − hT )/hS , hB−hS
DC = C0
tq
hT V
(9)
6
7
h* e−(Vt − hT )/hB − hB e−(Vt − hT )/hS . hB − hS
(10)
The discontinuity reaches its maximum value at time tD =
6
0 17
hB2 1 hS hB hT − ln V (hS − hB ) h* hS
,
(11)
before decaying towards the limit DC : CB . For deep bottom layers (hB hS ), h* = hS + hB C* /C0 and tD =
6
0 17
hB2 1 hT + hS ln h* hS V
,
(12)
which continues to increase with hB . The mean concentration over the entire water column can also be determined from eqns (3), (6) and (9),
(hS Q z Q hS + Vt)
(hS + Vt Q z Q hS + hT )
hT V
where h* = hS + (hB − hS )C* /C0 . Because incoming particles are immediately mixed through the bottom layer, a discontinuity in concentration can develop at the base of the thermocline even when the initial distribution is continuous. This is first evident at t = hT /V, from which its magnitude develops as
while the corresponding flux of particles into the thermocline is FS = VC0 e−Vt/hS .
(8)
which combined with (7) has the solution
provided the Peclet number is small (Smith, 1982; Martin & Nokes, 1988). This can be integrated to give CS = C0 e−Vt/hS ,
67
C =
C0 h* + C* (hS + hT − Vt) H
=C0
6
0 1 tQ
hT V
70 1
h* hB e−(Vt − hT )/hB − hS2 e−(Vt − hT )/hS H(hB − hS )
tq
hT . V (13)
0 1
h tq T . V
(7)
The mean concentration therefore falls linearly for a period determined by the thermocline thickness, followed by a transition to an exponential decay.
. . .
68 0
(a) t=1.1H/V t=0.5H/V
z/H
0.2
t=0.1H/V
0.4 0.6 0.8 1.0 0
0.2
0.4
0.6
0.8
1.0
C/C o 0
(b) t=1.1H/V t=0.5H/V
z/H
0.2
t=0.1H/V
0.4 0.6 0.8
begins to decay as particles first enter the bottom mixed layer (t q hT /V). It finally disappears at time tD , after which the behaviour becomes qualitatively similar to that for uniform initial concentration. However, when the bottom layer is deep (hB hS ), eqn (12) indicates that the concentration maxima can persist for long periods and therefore should be a relatively common feature in natural systems which satisfy this criteria. The size of the discontinuity calculated from (10) is shown by the dashed lines in Figs 3(a) and 3(b). In both cases, the discontinuity first appears at t = hT /V when particles from the surface layer first reach the bottom mixed layer. For an initially uniform distribution, concentrations in the thermocline fall more rapidly than those in the bottom layer, which results in an increase in DC. However, as particles continue to be lost from the bottom layer, the discontinuity begins to decay towards the limit DC : CB . The size of the DC peak increases with the
1.0 0
0.2
0.4
0.6
0.8
1.0
C/C o F. 2. Vertical profiles of particle concentration at non-dimensional times tV/H = 0.1, 0.5 and 1.1, for stratification defined by hS /H = 0.40, hT /H = 0.35 and hB H = 0.25. (a) Initially uniform particle concentration (C* = C0 ). (b) All particles initially concentrated in the surface layer (C* = 0).
(a) C T/C o
0.6
C B/C o
C/C o
C S /C o 0.2
∆C/C o
–0.2 –0.6 –1.0 0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
tV/H 1.0
(b) C T/C o
0.6
C/C o
Examples of vertical profiles calculated for the two initial conditions are shown in Fig. 2, with layer depths specified in the caption. Corresponding time series for concentrations within the two mixed layers and mid-depth in the thermocline are shown in Fig. 3. Concentrations are uniform in the well-mixed surface layer and then increase exponentially with depth through the thermocline. The distributions are identical for the two initial conditions to a depth of hS + Vt, as required by (6). For uniform initial concentration [Figs 2(a) and 3(a)], fluxes through the base of the thermocline initially balance losses to the bottom boundary layer, so that the bottom layer concentration is maintained at C0 . However, for t q hT /V a discontinuity develops at the base of the thermocline where particles are suddenly mixed through the bottom layer. Very similar results are obtained from numerical solutions of the primitive flow equations incorporating turbulent transport of settling particles (Appendix A). The development is quite different when particles are initially concentrated in the surface layer [Figs 2(b) and 3(b)]. Of particular significance is a sharp concentration maxima in the thermocline peaking at z = hS + Vt for t Q hT /V. The feature
1.0
C B/C o
0.2
C S /C o
∆C/C o
–0.2 –0.6
–1.0 0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
tV/H F. 3. Time series of particle concentrations in the surface mixed layer, bottom mixed layer, and midpoint in the thermocline. The size of the concentration discontinuity at the base of the thermocline is shown by the dashed curve. The stratification is the same as that in Fig. 2. (a) Initially uniform particle concentration (C* = C0 ). (b) All particles initially concentrated in the surface layer (C* = 0).
69
1.0
Population Dispersion: Horizontal Advection and Downstream Deposition
0.8
If populations are generated or suspended within a localised region, the horizontal flow will tend to deposit particles with some characteristic downstream distribution. Consider for example, a horizontal flow of uniform velocity u in the x direction (neglecting boundary layer shear), with populations defined by (1) at a position x = 0. The deposition flux is then given simply by FB = VCB , which after a transformation to spatial coordinates (t : x/u) in eqn (9) gives
C/C o
0.6 Fully mixed limit
0.4
0.2
0
C/C o
Stratified limit 0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
0
FB = VC*
tv/H F. 4. Time series of depth averaged particle concentration for the stratification in Fig. 2, along with the fully turbulent (hS = H) and fully stratified (hT = H) limits.
=VC0
xQ
uhT V
1
6
70
h* e−(Vx/u − hT )/hB − hS e−(Vx/u − hT )/hS hB − hS
xq
1
uhT . V (14)
For a particle cloud of finite width Dx, the final distribution of deposited particles is FB Dx/u. If the non-dimensional deposition flux, FB /VC0 , is plotted against the non-dimensional downstream distance, xV/uH, the curves are identical to those for CB /C0 in Fig. 3. For a vertically uniform upstream distribution, Fig. 3(a) suggests uniform deposition for a distance uhT /V, followed by an exponential decay. However, if particles are initially deposited within the upper layer [Fig. 3(b)], then there is no deposition over the
1.0
0.8 1.0
1.0
0.6 h T/H
ratio hB /hS , with the limits DC : 0 as hB /hS : 0 and DC : C0 as hB /hS : a. Stratification not only affects the concentration distribution, but also the depth averaged concentration. This quantity has been calculated from (13) and plotted as a function of time in Fig. 4. The corresponding trends for fully mixed and fully stratified conditions are also shown for comparison. The presence of stratification clearly produces a more rapid decline in concentration levels, which strongly impacts long-term population levels (t q H/V). In this example, there is almost an order of magnitude difference in the populations of the partially stratified and fully turbulent flows at t = 2H/V. The depth averaged concentrations in Fig. 4 can be generalised to other stratifications by using (13) to calculate the time, t90 , required for 90% of the initial concentration to be deposited. Vt90 /H is contoured in Fig. 5 as a function of the surface layer and thermocline depths. Particle suspension times decrease as the thickness of the thermocline expands, with Vt90 /H falling by a factor of 2.5 from fully mixed to fully stratified conditions. This factor grows with deposited fraction, so that the time required to deposit 99% of the initial concentration in a fully turbulent environment is 4.6 times longer than in a fully stratified environment. Figure 5 further indicates that the presence of even a thin thermocline results in a large reduction in the suspension time. This is because the e-folding time for particle decay in each turbulent mixed layer is proportional to the depth of that layer.
0.4 1.4 1.3 1.2
0.2
1.7 1.6 1.8 2.0 1.9
1.2 1.3 1.4
1.4 1.5
1.6 1.7 1.9 1.8 2.0
0 0
0.2
0.4
0.6
0.8
1.0
h S/H F. 5. Contour plot of non-dimensional time (tV/H) required for 90% of the particles initially present to settle. This is shown as a function of surface mixed layer depth and thermocline depth.
. . .
70
distance uhT /V, following by a relatively rapid increase to a peak deposition at the location utD (C* = 0).
(a) 3 2.0
1.5
2
The theory can now be extended to organisms characterised by a growth rate m, which encompasses all physiological and environmental factors except the tendency to settle through their fluid environment. Only species in which settling motions dominate over swimming or buoyancy regulation are considered.
C/C o
Population Growth
1.1 1 0.9 G=0.0 0.1
0.5
0 0
0.4
0.8
1.2
0
1
1.6
2.0
2.4
2.0
2.4
2.0 1.5 1.1 0.9
2
(15)
G=0.0
0 0
0.4
Many plankton species have a tendency to float towards higher light levels and in some cases settle along the surface as a scum. The theory can be applied directly to this situation by inverting the model so that the layer denoted by subscript B becomes the surface layer with growth. This approach neglects growth within the thermocline, so that the populations below the surface layer are identical to
0.8
1.2 tV/H
(c) 3
2.0 1.5
(16)
where G = m¯ hS /V is a measure of the relative importance of growth and settling. Clearly, populations grow for G q 1 and decay for G Q 1, as shown in Fig. 6(a). If growth below the surface layer is negligible due to light limitation, then the development of the population at these levels is equivalent to that for inert particles with a reduced settling speed (1 − G)V. These results can be applied to time dependent systems by calculating solutions incrementally while continuously updating m¯ , hS and C0 (Sherman & Webster, 1994).
0.5
0.1
1
2
1.1
1
0.9
C/C o
C S = C0 e
2.4
(b) 3
hS where m¯ = h−1 S f0 m dz is the mean growth rate for the layer (m¯ is calculated in Appendix B for the case of light limited growth in aquatic environments). Over time-scales for which m¯ and hS are approximately constant, this has the solution −(1 − G)Vt/hS
2.0
C/C o
Consider the case where surface layer populations grow exponentially and are only limited by their tendency to sink. This balance often controls phytoplankton populations in species with negative buoyancy. The surface layer population is then described by V dCS C , = m¯ − hS S dt
1.6
tV/H
G=0.0
0.1
0.5
0 0
0.4
0.8
1.2
1.6
tV/H F. 6. Time series of particle populations with growth defined by the parameter G = m¯ hS /V and stratification as in Fig 2. (a) Particles sinking out of the surface layer. (b) Particles initially uniformly distributed (C* = C0 ), floating through the surface layer and forming a scum. (c) Particles initially concentrated within the bottom layer, floating into the surface layer and forming a scum.
those of inert particles. The surface layer can be described by dCB FT − (V − m¯ hB )CB = , dt hB
(17)
71
which combined with (6) has the solution
CB = C*
6
70 1
1 − G e−(1 − G)Vt/hB 1−G
tQ
hT V
[C* hB − (C* − C0 )(1 − G)hS ] e−((1 − G)Vt − hT (1 − G))/hB /(1 − G) − C0 hS e−(Vt − hT )/hS CG − * e−(1 − G)Vt/hB hB − hS (1 − G) 1−G
=
0 1 tq
hT V
(18) where G = m¯ hB /V. In the limit of zero growth (G = 0), eqn (18) reduces to (9). Figure 6(b) shows the growth curves for the case C* = C0 and a range of G values. In all cases there is an initial increase in population at the natural growth rate m¯ . This is because losses to the surface are initially exactly compensated by an influx from the thermocline. However, as this flux diminishes in accordance with (7), growth is only maintained for G q 1, with decaying populations for G Q 1. In some circumstances, particles reaching the surface are stirred back into the layer rather than forming a scum. This usually applies to phytoplankton under relatively strong wind conditions (Hutchinson & Webster, 1994). The settling term in (17) can then be omitted and the solution takes the form CB = C* 4(1 + G−1 ) e GVt/hB − G−15 −1
0 1 tQ
hT V
−1
=(C* (1 + G ) − [C* G − C0 (G + hB /hS )−1 ] e−GhT /hB ) e GVt/hB − C0 (G + hB /hS )−1 e−(Vt − hT )/hS
0 1 tq
hT . (19) V
Since there are no loss terms, this solution describes monotonic exponential growth for all values of G. The growth rates are not shown, but are very similar to those for the scum forming particles characterised by a higher growth factor G + 1 in Fig. 6(b). Seeding populations can also rise from depth, after being resuspended from the bottom. For example, some cyanobacteria produce resting spores (akinetes) or have overwintering benthic phases. The inverted model using (18) or (19) with C* = 0 can be easily applied to this situation, with C0 representing the seed population. Figure 6(c) shows the growth curves for a range of G values. Particles first reach the surface layer at t = hT /V and then begin to multiply. However, with settling into a surface scum, populations are again only sustained for G q 1. Without scum formation (not shown), populations grow exponentially for all positive values of G.
/ The manner in which two plankton populations compete for available chemical resources (nutrients and inorganic carbon) is highly species dependent. However, a qualitative understanding can be gained from a simple Lotka–Volterra model. Once any fluxes from the thermocline have declined, the population of each species can be described by (15), with each m¯ a function of its own concentration (intraspecific competition) and the concentration of other species (interspecific competition). The simplest case is a linear dependence of growth rates on the population of two species (denoted by subscripts 1 and 2). The stability of such a system is well known (MaynardSmith, 1974) and can be expressed in terms of the parameter R=
m¯ 1 hS − V1 V1 (G1 − 1) = , m¯ 2 hS − V2 V2 (G2 − 1)
(20)
where G1 = m¯ 1 hS /V1 and G2 = m¯ 2 hS /V2 . Only species 2 will persist in the surface layer if G2 q 1 and R Q (1m¯ 1 /1C2 )/(1m¯ 2 /1C2 ); the two species will coexist if G1 q 1, G2 q 1 and (1m¯ 1 /1C2 )/(1m¯ 2 /1C2 ) Q R Q (1m¯ 1 / 1C1 )/(1m¯ 2 /1C1 ); while only species 1 will persist if G1 q 1 and R q (1m¯ 1 /1C1 )/(1m¯ 2 /1C1 ). While the growth derivatives can only be determined empirically, the analysis at least demonstrates that R is the critical parameter controlling competition in this type of system.
Implications Large peaks in concentration of suspended material are often observed to coincide with zones of strong density gradient. This is sometimes associated with favourable light or nutrient levels, such as in the case of phytoplankton maxima near the base of the ocean thermoclines (Sharples & Tett, 1994; Denman &
72
. . .
Gargett, 1995). However, similar observations have also been made of non-reproductive particles, such as larvae and marine snow (aggregates of inorganic particles, microorganisms and detritus), which are passively advected by the flow (Alldredge & Cohen, 1987; MacIntyre et al., 1995; Tremblay et al., 1994). A number of explanations have been suggested for particle accumulation in stratified regions. These include horizontal intrusions of particle laden water, reduction of settling speeds as higher densities are encountered, and increased particle aggregation in turbulent zones. However, the density conditions and particle characteristics are often not conducive to such mechanisms. Furthermore, MacIntyre et al. (1995) have used a numerical random walk model, as described by Patterson (1991), to demonstrate that accumulations occur even in the absence of these processes. The theory presented here reveals the underlying mechanism responsible for the behaviour of the random walk model. Specifically, the concentration within the overlying turbulent layer decays uniformly, while the peak in the pycnocline moves downward at the particle settling speed. There is no decrease in the magnitude of the peak until it reaches the turbulent layer below. Furthermore, if this turbulent layer is much deeper than the surface layer, then the pycnocline peak will be a very persistent feature of the concentration profile [eqn (12)]. This may explain why such features are so frequently observed in oceanic systems. The theoretical results have particularly significant implications for the survival and dispersion of species which rely on suspended seeding populations to migrate between isolated patches of habitat or recolonise an environment each season. This is due to the strong dependence of population decay rates on stratification. If the water or air column is turbulent, then small numbers of particles can remain suspended for extended periods. This may allow micro-organisms to overwinter in suspension, or assist seeds or larvae in reaching more distance habitats. However, stratification tends to constrain this mechanism and may be responsible for anomalously poor seasons or reduced colonisation ranges. Figure 6 demonstrates the strong dependence of growing populations on the parameter G for quite a diverse range of circumstances. A particular example to which this may be relevant is the relative growth rates of diatoms and cyanobacteria (blue-green algae)
within Australia’s Murray–Darling river system. The numerical model described in Appendix A has been successfully used to make detailed predictions of these populations, including the formation of blooms (Bormans & Webster, 1996). However, the discussion here will focus on testing the underlying conditions required by the theory for one of the plankton species to proliferate. In Fig. 7, a four month time series of diatom populations measured in the lower Murrumbidgee River is compared with the associated time series of the growth factor G. Initially, both G and the diatom population are relatively high. However, a period of low flow extending into December allows the development of stratification, with a corresponding sharp reduction in G. Settling then dominates over growth and the population falls. When G again increases with the flow rate, the population recovers rapidly. While generally remaining high until the next sustained reduction in G around mid-February, there are small population drops associated with shorter term (H/V 1 5 days) reductions in G. Daily oscillations in G (associated with cloudy days) do not have a significant effect, because plankton settling through the thermocline are re-entrained into the deepening mixed layer before they can reach the bottom. The strong correlation between population changes and G (over time-scales greater than H/V) appears to support the theoretical framework and suggests that G might be a useful indicator of bloom formation in other populations limited by their settling behaviour.
Conclusion A simple theory has been developed to describe particle settling in stratified systems. It is based on the premise that settling in turbulent regions obeys a well known exponential decay law based on the turbulent layer depth, while settling in stratified regions is at the particle settling speed. This can lead to a range of interesting phenomena, such as discontinuities and sharp peaks in concentration profiles. The inclusion of particle population growth in the turbulent surface layer reveals that G = m¯ hS /V is the critical parameter controlling growth of populations limited by settling. There remains a large range of related phenomena to be investigated. These include other initial conditions (selected for specific applications), time dependent stratification, z dependent settling speeds (including inhibition by large particle concentrations), and growth within stratified regions. It is hoped that this paper will provide a theoretical basis for studying
73 10 000
1.6
1.4
G
1.0 100
0.8 0.6
Diatoms (cells ml –1)
1000
1.2
10
0.4 0.2
1 06/03/95
26/02/95
18/02/95
10/02/95
02/02/95
25/01/95
17/01/95
09/01/95
01/01/95
24/12/94
16/12/94
08/12/94
30/11/94
22/11/94
14/11/94
0
F. 7. Depth averaged concentrations of the diatom Aulacoseira sp. determined from cell counts collected every three days in the lower Murrumbidgee River in southeastern Australia between 14 November 1994 and 8 March 1995 (w). This is plotted alongside the growth factor G (W), which was calculated at two hourly intervals from (B.3) and then diurnally averaged. The calculation used mixed layer depths based on two hourly averaged temperature profiles from thermistor chain measurements, along with locally measured I0 (t) and h. Details of the field measurements are given by Webster et al. (1996). From laboratory and field data, the diatom characteristics were estimated to be V = 1.0 m day−1, mmax = 0.7 day−1 and Ik = 10 mEin m−2 s−1 (Reynolds, 1984). While there are significant uncertainties in these estimates, the trends revealed by the plot are not sensitive to their specific values.
these issues within the context of more specific applications. We thank Dr Gary Jones, CSIRO Land and Water, for the diatom count data from an NRMS funded project examining cyanobacteria blooms in the lower Murrumbidgee River. REFERENCES A, A. L. & C, Y. (1987). Can microscale patches persist in the sea? Microelectrode study of marine snow, fecal pellets. Science, 235, 689–691. B, A. F. & M, G. L. (1983). Diagnostic and prognostic numerical circulation studies of the South Atlantic Bight. J. Geophys. Res. 88, 4579–4592. B, M. & W, I. T. (1996). Modelling the dynamics of the diatom Aulacoseira sp. and the blue-green Anabaena in relation to river flow. CSIRO Centre for Environmental Mechanics Technical Report No. 127.. D, K. L. & G, A. E. (1995). Biological-physical interactions in the upper ocean: The role of vertical and small scale transport processes. Ann. Rev. Fluid Mech. 27, 225–255. H, W. C. (1982). Aerosol Technology: Properties, Behaviour, and Measurement of Airborne Particles. New York: Wiley. H, P. A. & W, I. T. (1994). On the distribution of
blue-green algae in lakes: Wind-tunnel tank experiments. Limnol. Oceanogr. 39, 374–382. K, R. C. (1991) Erosion of a stable density gradient by sedimentation-driven convection. Nature 353, 423–425. K, R. C. & K, G. S. (1997). Particle settling through a diffusive-type thermocline staircase in the ocean. Deep Sea Res. 44, 399–412. MI, S., A, A. L. & G, C. C. (1995). Accumulation of marine snow at density discontinuities in the water column. Limnol. Oceanogr. 40, 449–468. M, D. & N, R. (1988). Crystal settling in a vigorously convecting magma chamber. Nature 332, 534–536. M-S, J. (1974). Models in Ecology. Cambridge: Cambridge University Press. M, G. L. & Y, T. (1974). A hierarchy of turbulence closure models for planetary boundary layers. J. Atmos. Sci. 31, 1791–1806. M, G. L. & Y, T. (1982). Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. Space Phys. 20, 851–879. O, A. (1980). Diffusion and Ecological Problems: Mathematical Models. Berlin: Springer-Verlag. P, J. C. (1991). Modelling the effects of motion on primary production in the mixed layer of lakes. Aquat. Sci. 53, 218–238. R, C. S. (1984). The Ecology of Freshwater Phytoplankton. Cambridge: Cambridge University Press. S, J. & T, P. (1994). Modelling the effect of physical variability on the midwater chlorophyll maximum. J. Mar. Res. 52, 219–238.
. . .
74
S, V. S. & W, I. T. (1994). A model for the light-limited growth of buoyant phytoplankton in a shallow, turbid waterbody. Aust. J. Mar. Freshwater Res. 45, 847–862. S, I. R. (1982). A simple theory of algal deposition. Freshwater Biol. 12, 445–449. T, S. A. (1987). Diapyanal mixing in the thermocline: A review. J. Geophys. Res. 92, 5231–5248. T, M. J., L, J. W., W, F. E., N, C. E., P, F. H. & S M. M. (1994). Drift of sea scallop larvae Placopecten magellanicus on Georges Bank: a model study of the roles of mean advection, larval behavior and larval origin. Deep Sea Res. 41, 7–49. W, I. T., J, G. J., O, R. L., B, M. & S, B. S. (1996). Control strategies for cyanobacterial blooms in weir pools. CSIRO Centre for Environmental Mechanics Technical Report No. 119.
Figure A1(b) shows the intensity of the turbulence throughout the water column in terms of the Peclet number, Pe = VH/k, based on the total water depth. This will help interpret comparisons of concentration profiles from the theoretical and numerical solutions [Fig. A1(c)]. In the surface mixed layer Pe1, so that concentrations are always homogeneous with good agreement between theoretical and numerical esti
0
(a) 1.20 0.60 0.24
0.2
0.12
An important assumption of the theory is that the turbulent and stratified regions are respectively characterised by small and large Peclet numbers. To investigate this issue within the context of more realistic flows, the theoretical results have been compared with those from a numerical model. This has been achieved by incorporating settling particles into the primitive equation circulation model of Blumberg & Mellor (1983). This solves the nonlinear hydrostatic form of the flow equations using standard finite difference methods. Characteristics of the subgrid scale turbulence were determined using the turbulence closure scheme developed by Mellor & Yamada (1974, 1982). Flow fields and particle distributions have been calculated from this model for a small idealised lake, forced by wind stress and heat flux through the free surface. The numerical runs were conducted in a basin of depth H = 3.0 m and fetch 100 m, while neglecting gradients in the cross-wind direction. The primitive equations were solved on a uniform grid with 0.1 m vertical resolution and 0.5 m horizontal resolution. The initial stratification consisted of an upper mixed layer of depth hS = 1.14 m and temperature 30°C, a thermocline of depth hT = 1.11 m and linear temperature gradient, and a bottom mixed layer of depth hB = 0.75 m and temperature 20°C [Fig. A1(a)]. There was initially a uniform distribution of particles, with uniform settling velocity V = 1 × 10−5 ms−1. The flow was forced by a constant windstress and surface heat flux. These were tuned so that the depths of the three layers remained approximately constant over the duration of the run [Fig. A1(a)], while also maintaining turbulent flow within the two mixed layers. This system could be compared directly with theoretical results for inert particles.
z/H 0.6 0.8 1.0 15
20
25
30
35
Temperature (°C) 0
(b)
0.2 0.12
0.4
0.24
z/H
Comparison With a Numerical Solution
0.4
0.60 1.20
0.6 0.8 1.0 1.0e–02
1.0e+00
1.0e+02
1.0e+04
Pe 0 0.2
(c) 1.20
0.60
0
0.2
0.24
0.12
0.4 z/H
APPENDIX A
0.6 0.8 1.0 0.4
0.6
0.8
1.0
C/C o F. A1. (a) Development of the stratification in the numerical model with non-dimensional time, tV/H. (b) Development of the rms turbulent velocity profile with non-dimensional time. (c) Comparison of the numerical and theoretical concentration profiles at four non-dimensional times.
APPENDIX B In many instances, population growth in aquatic environments is limited by the available light. This can be described in terms of an instantaneous growth rate mmax I , I + Ik
0.8 100 0.6 10 0.4 1.0 0.2 0 0 0
2
(B.1)
where mmax is the maximum growth rate, Ik is the half-saturation constant for light limited growth, and I is the photosynthetically active radiation (PAR) at depth z. Provided the particles do not contribute significantly to light attenuation, the PAR can usually be approximated by Beer’s Law, I = I0 e−hz, where I0
4
6
8
10
ηh S F. B1. Average surface layer growth rates as a function of layer depth (non-dimensionalised by the light extinction coefficient) for a range of I0 /Ik values.
is the value at z = 0, and h is the attenuation coefficient. The mean growth rate is then m¯ =
Light Limited Growth
m=
1.0
G/G max
mates. Within the thermocline Pe1, and similar concentration gradients are predicted by the two models. There is some tendency for particles to settle more rapidly in the numerical case as the high Peclet number region penetrates into the bottom mixed layer. In any case, the range of bottom layer Peclet numbers (10−2 Q Pe Q 1) encompasses settling behaviour which is significantly more rapid than the small Pe limit. This leads to quite significant discrepancies between theoretical and numerical results within this layer. In most respects, the theory provides a very good description of the settling behaviour in the numerical model. However, it does provide a cautionary reminder that a well-mixed density field does not always imply strong turbulence. Since it is usually not practical to directly measure turbulent velocities, small Peclet numbers should only be assumed where there is an obvious source of turbulent kinetic energy, such as a strongly sheared boundary layer. While this was not present in the model lake, other systems such as rivers, tidal flows and atmospheric boundary layers are much more likely to satisfy this requirement.
75
1 hS
g
hS
m dz =
0
0
1
mmax I0 + Ik ln . hhS I0 e−hhS + Ik
(B.2)
The critical parameter G = m¯ hS /V determining the population behaviour can then be expressed in the form m I0 + Ik G max ln . (B.3) hV I0 e−hhS + Ik
0
1
Figure B1 is a graphic representation of (B.3) indicating that growth is favoured by strong incident radiation (large I0 /Ik ) with low attenuation through the surface layer (small hhS ). In the weak radiation limit (I0 /Ik : 0) G : 0, while in the strong radiation limit (I0 /Ik : a) G : Gmax where Gmax = mmax hS /V is independent of I0 .