The effect of stochastic perturbations on plankton transport by internal solitary waves

<jats:p>Abstract. Internal solitary and solitary-like waves are a commonly observed feature of density stratified natural waters, including lakes and the coastal ocean. Since such waves induce significant currents throughout the water column they can be responsible for significant transport of both passive and swimming biota. We consider simple models of moving zooplankton based on the Langevin equation. The small amplitude randomness significantly alters the nature of particle motion. In particular, passage through the wave leads to strongly non Gaussian particle distributions. When the plankton swims to return to its equilibrium photic level, a steady state that balances randomness, swimming and wave-induced motions is possible. We discuss possible implications of this steady state for organisms that feed on plankton.                     </jats:p>

section_content: Due to the changes in water density caused by the variable distribution of both temperature and salinity, the interior of the world's lakes, seas and oceans serves as a waveguide for a variety of wave motions for which gravity provides the restoring force (Gill, 1982). The focus of the present work will be on those motions which exhibit high enough frequencies so that the effects of the rotation of the Earth may be neglected. These are commonly known as internal gravity waves, and the particular internal waves we will focus on are finite amplitude internal solitary or solitary-like waves. The waves we consider propagate horizontally, largely without changing form, and their surface manifestations, alternating bands of light and dark, may often be observed from airplanes and spacecraft (Helfrich and Melville, 2006). The reason for the ease of observation is the systematic modulation Correspondence to: M. Stastna (mmstastn@uwaterloo. ca) of surface wave activity by the convergence and divergence of internal solitary wave-induced currents. Regions of high surface wave activity are visible as dark bands, while regions of low surface wave activity are visible as light bands in the visible spectrum. Internal solitary waves are generated by a variety of mechanisms. In the ocean, the most prominent is tidally induced flow over topography. As such, internal solitary waves are often visible at locations in which the ocean depth varies quickly, for example in straits (e. g. the Strait of Gibraltar) and fjords (e. g. Knight Inlet, Canada) (Helfrich and Melville, 2006; Farmer and Armi, 1999). In lakes, internal solitary-like waves can be produced by the breakdown of wind-induced, lake-scale internal seiches (de la Fuente et al., 2008). 

From a theoretical point of view, internal solitary waves (henceforth ISWs) are a translating, baroclinic vorticity distribution of permanent form (though since the streamlines are not closed, different fluid particles make up the wave at any given time) and can be described mathematically by a variety of theories. Weakly nonlinear theories explicitly balance wave steepening due to nonlinearity with dispersion in a model equation that considers horizontal and temporal variations. The model equation is derived by positing solutions that are separable in the vertical and horizontal spatial variables, then performing an asymptotic expansion in the amplitude and aspect ratio parameters. This procedure generally yields completely integrable equations in the Korteweg deVries, or KdV, family (Grimshaw, 1997). While these theories are mathematically rich, they often misrepresent the spatial structure of finite amplitude waves (Lamb, 1999). Indeed, outside of highly specialized situations this is true even for moderate amplitudes, though, to be fair, model equations are commonly used to successfully interpret various features of field data (e. g. mass transport by ISWs on the Malin shelf, Inall et al., 2001). An alternative to model equations, is to give up on the many attractive mathematical properties of the weakly nonlinear theories, such as an infinite hierarchy of conserved quantities and an inverse scattering transform, and consider fully nonlinear ISWs which are solutions of the full Euler equations that govern an inviscid fluid. In the oceanographic literature this approach is often carried out using numerical solution of the Euler equations (Vlasenko and Hutter, 2002). Semi-analytical results can be achieved by deriving a single equation for the displacement of lines of constant density (or isopycnals) and the unknown wave speed. The resulting nonlinear, elliptic eigenvalue problem is known as the Dubreil-Jacotin Long (henceforth DJL) equation. It has been shown that fully nonlinear ISWs are characterized by a variational structure in which the available potential energy is held fixed and the kinetic energy of the disturbance is minimized (Turkington et al., 1991). This facilitates the construction of a computational algorithm based on standard techniques in continuous optimization theory for the computation of ISWs. We have implemented the algorithm using standard spectral methods (Trefethen, 2000), with details of an open source, MATLAB code described in Dunphy et al. (2011). The currents induced by fully nonlinear solitary waves are often sizeable, and vary over the entire depth of the water column. They are thus important to the distribution of various types of passive and semi-passive tracers including pollutants, sediment and biota (especially plankton). Studies of deterministic internal solitary wave-induced particle transport (Lamb, 1997) have shown that horizontal, wave-induced transport is strongly dependent on initial height and wave amplitude. Indeed, Lamb showed that transport predicted by model equation-based theories (e. g. those leading to the KdV equation) can be significantly different than that due to exact solitary waves. The effect of internal waves on plankton patchiness was studied using currents predicted by linear and weakly nonlinear theory and an Eulerian formulation of plankton concentration by Lennert-Cody and Franks (1999). The authors found that wave-induced concentration changes were largest at the wave crest and concluded that measurements of patchiness should be conducted along with measurements of lines of constant density (isopycnals). These conclusions remain valid for exact ISWs, even if the quantitative predictions based on linear and weakly nonlinear theories do not. 

Since plankton are microscopic creatures, the extent of their motion is not large in its spatial extent. Compared to fluid motions occurring on larger length scales (e. g. internal waves) this motion is erratic in time, and can be described by various complex models from the literature (Genin et al., 2005; Rubjakov, 1970). The large number of such creatures in a cubic meter of water, along with their irregular individual trajectories, suggests that their interaction with waveinduced currents can be idealized using a stochastic Lagrangian particle description. We will pursue such a description in this manuscript. An illustrative example of a point release experiment, with the initial position of the 10 000 member particle ensemble upstream of a rightward propagating wave of depression and at the center of the region of rapid density change (or pycnocline), is shown in the upper panel of Fig. 1. The wave in the figure induces horizontal currents that are oriented in (against) the direction of wave propagation above (below) the pycnocline. While the vertical currents induced by the wave are not negligible, they are typically less than a quarter of the horizontal currents. In the upper panel of Fig. 1 we show the shaded contours of density and four instances of the particle cloud (shown every 16. 4 dimensionless time units) as it passes through the wave. Length is scaled by the fluid depth H, and time by the advective time scale H /c where c is the wave propagation speed. The aspect ratio of the upper panel in Fig. 1 is somewhat misleading since the particles initially diffuse isotropically. In the lower four panels we account for the "long" wave nature of the solitary wave by showing the particle clouds with a one to one aspect ratio, following the same labeling for the particle clouds as in the upper panel. It can be clearly seen that the dominant wave-induced currents are horizontal currents, with a shear layer across the pycnocline (note the change in scale from panel 1 to panel 4). 

The remainder of the manuscript is organized as follows. We begin by discussing the simplest, most symmetric stochastically perturbed motion, that due to a stationary vortex with circular streamlines. This example allows us to introduce the theoretical and numerical methods employed. Using both analytic and geometric arguments we show that stochastically perturbed motion in the circular vortex must lead to a net (or deterministic) drift that is radially outward and that the rate of drift decays with radial distance from the center of the vortex. We proceed to discuss how stochastically perturbed motion due to ISWs breaks the symmetry exhibited by the circular vortex and show the unexpected consequences of stochasticity in this case, namely the strongly non Gaussian particle distribution during and after passage through the wave. Since purely stochastic motion does not allow for a nontrivial steady state, we introduce the effects of combined random and systematic swimming. While swimming behaviour is complex, and many varieties are discussed in the biology literature (e. g. Genin et al., 2005; Rubjakov, 1970), plankton generally seeks to maintain its position at a fixed photic (light) level. We thus consider a combination of stochastic, wave-induced and return to photic level motions. In this case a non-trivial steady state is possible. We present simulations of this steady state, and discuss the implications of the resulting plankton distribution for plankton feeders.
section_name: Introduction

section_content: In the absence of perturbations (e. g. those due to particle swimming, treated in Sect. 2. 3), pathlines followed by a particle in the fluid can be found by integration of the fluid velocities (u,w). These may be specified analytically (as for the irrotational vortex case below) or numerically (as for the ISW case below). For steady flow, such as a circular vortex x ig. 1. The motion of a 10,000 member ensemble of particles initially located upstream of the ISW and at the center of the pycnocl uperimposed on a shaded plot of density contours. Subsequent times evolve from right to left and are separated by 16. 4 time units. In pper panel the entire horizontal extent of the ISW is shown, however only a portion of the computational domain is shown. Since ISWs ong waves, panels (1)-(4) show a 1:1 aspect ratio of the particle clouds with the same labels as in the upper panel. Note the increase in sc s time increases. or that induced by the ISW in a frame moving with the wave speed, the pathlines match the shape of the streamlines. We compute pathlines for a flow in which the stochastic perturbations are taken to effect the velocity of individual particles, but not the fluid velocities themselves. The governing equations in the case of coloured noise are given by

where the correlation time of the two noise terms ξ i is specified by γ i > 0 and the variance by σ i > 0. These must be specified as part of the model. Here dW represents the increment of the Wiener process, and essentially reflects the fact that stochastic differential equations require a careful mathematical interpretation due to the non-smooth nature of the noise (Gardiner, 1990). There is a considerable literature on the mathematics, approximation and numerical solution of stochastic differential equations (Kloeden and Platen, 1992; Ottinger, 1996; Gardiner, 1990). We will adopt the so-called Langevin approach. In other words, we consider individual integrations of Eqs. ( 1 ) and ( 2 ). In the following we apply the Milstein (Ottinger, 1996) method for the discretization of the stochastic differential equation and the Mersenne Twister algorithm (Matsumoto and Nishimura, 1998) to generate realizations of the stochastic perturbation. The numerical algorithm employed in the integration of the equations is of order t in time in both the deterministic and stochastic parts of the governing equation and converges strongly (Ottinger, 1996). The strong convergence implies we can make mathematically valid conclusions regarding individual paths (for example we can plot accurate approximations of individual paths). The Mersenne Twister algorithm allows for the rapid construction of large ensembles while maintaining confidence in the temporal characteristics of the stochastic perturbations. The second approach to dealing with stochastic differential equations forgoes the individual paths of the Langevin approach in favour of the Fokker-Planck equation, a partial differential equation, governing the probability density. Under certain mathematical conditions, such as the case of additive white noise (see Ottinger, 1996 or Gardiner, 1990 for details), the two approaches are equivalent. However, in many instances, the former approach is preferable from a computational viewpoint (see the discussion in Ottinger, 1996) and this is especially true for exotic types of noise in biological applications such as the run and tumble model (Bearon, 2007). While the present model allows for a variety of correlation time-variance combinations, throughout the following we will assume that both noise terms have the same correlation time and variance. The coloured noise model approaches the classical white noise limit, though the white noise limit requires an appropriate choice of stochastic calculus (see Gardiner, 1996 for example). While we have performed white noise simulations corresponding to many results reported on below, for numerical simulations we report only on the coloured noise case which we believe to be more relevant to the motions of plankton and other biota.
section_name: Results

section_content: We consider a steady flow with closed, circular streamlines. The flow is along streamlines and a particular type of vortex is distinguished by the radial structure of the azimuthal velocity field. In polar coordinates we write

where êθ is the unit vector in the azimuthal direction. We have considered several standard models of vortices (irrotational, Rankine, etc. ) but have found no qualitative differences in behaviour. We thus focus on the simple irrotational vortex for which the velocity field is given by

Corresponding Cartesian expressions for the velocity field u i = (u(x,z),w(x,z)) follow from a standard change of variables. 

The study of vortices in two dimensions is a well developed discipline with many topics covered in the classical review by (Aref, 1983). Dynamical systems theory has been applied to chaotic mixing by vortices in (Rom-Kedar et al., 1990), and more recently by (Vassilicos, 2002). Our approach is, in comparison, much simpler, serving to motivate the subsequent work on ISW-induced motion. 

In Fig. 2a the solid curve shows the effect of noise on the mean radial position of an ensemble of 10 000 particles initially found at r = 0. 2. It is clear that the mean radius gradually increases, with the rate of increase slowing as the mean radius increases. In the same panel we demonstrate that increasing the correlation time (dashed curve) while holding the intensity of the noise fixed leads to a faster rate of increase of the mean radius. In Fig. 2b we repeat the solid curve from 2a and show the effect of a one dimensionless unit increase in initial radius (dashed line). It is clear that the rate of increase of the mean radius is much lower for the ensemble of particles that begins farther from the center of the vortex. We observed similar behaviour in simulations with various vortex types (e. g. Rankine). 

In the white noise limit, if we assume the Itô stochastic calculus, we can derive an analytical estimate for the outward drift. Gardiner (Sect. 4. 4. 5, notation follows the reference) shows that the completely noise driven motion

with isotropic noise of amplitude can be written in polar coordinates (a(t),φ(t)) as

where

As Gardiner points out, ( 8 ) is an orthogonal transformation, so that dW a (t) and dW φ (t) can be taken as increments of independent Wiener processes. Thus the stochastic motion, through the Itô-calculus chain rule, is responsible for a deterministic outward drift term (i. e. the first term on the right hand side of Eq. 7). The case of any purely azimuthal motion, including the irrotational vortex, follows in the same manner (with somewhat more complicated algebra). 

In Fig. 3a we demonstrate the geometric reason for the stochastically induced outward drift. The unperturbed flow around the grey vortex is indicated by arrows. A circle, the interior of which represents all points within one standard deviation of the mean particle position is shown along with the bisector of the circle formed by the tangent line to the circle. It can be seen that more than one half of the area of the circle lies outside of the grey region. Thus the particles are more likely to move out of the grey circle, and the tangential speed of motion along the circle is irrelevant to the rate at which this occurs. As the radius of the grey circle (and hence the radius of curvature) increases, the fraction of the area of the circle corresponding to the standard deviation of the particle inside the grey circle tends closer and closer to one half, and hence the rate of outward drift decreases. The same general geometric principle holds for more complicated paths, as illustrated in Fig. 3b (for such cases an analytical argument analogous to that above is impossible to carry out). As the concavity of the path changes, so does the direction of stochastically induced drift. Moreover, without the symmetry of the circular path, the speed with which the particle moves along the curve becomes important, for example if the particle slows down in a region of small radius of curvature. Internal solitary waves, which are a baroclinic vorticity distribution of permanent form, do not have circular, or even closed, streamlines (essentially due to symmetry breaking by the stratification) and as such can be expected to exhibit more complicated behaviour than that found for the irrotational vortex.
section_name: Stochastically perturbed motion in a vortex

section_content: We performed computations of 10 000 particles initially released from the same location upstream (at the center of the pycnocline) of a large solitary wave, as shown in the upper panel of Fig. 1. The center of the pycnocline is chosen as the site of particle release since it is the region of maximum shear in the wave-induced currents. The solitary wave is a solution of the DJL equation,

an elliptic, nonlinear eigenvalue problem for the isopycnal displacement η and the wave propagation speed c that is equivalent to the full Euler equations governing an invisicid fluid. In all cases considered we make the Boussinesq approximation. The particular stratification profile, as a function of depth, is given by the functional form of N 2 (z), the so-called buoyancy frequency squared (Gill, 1982) which is defined from the background density profile ρ(z) as

where g is the acceleration due to gravity and ρ 0 is a reference density. Throughout we consider a simple, single pycnocline profile for which

and

where z 0 specifies the pycnocline center (taken to be 20 % of the total depth below the surface) and d specifies the pycnocline thickness (taken to be 5 % of the total depth). Once η and c are known the density field is given by ρ = ρ(z -η) and the wave-induced velocities in a frame moving with the wave are given by (u,w) = c(η z,-η x ), where subscripts denote partial derivatives. The vertical velocities are antisymmetric across the wave crest (x = 0), however the horizontal velocities do not have a simple line of symmetry, The axis labels, 'x' and 'z' have been removed to make larger sub-panels.
section_name: Stochastically perturbed motion through a solitary wave

section_content: Particle clouds with the mean position subtracted (white), the ellipse corresponding to the covariance matrix and 10 000 draws from a bivariate Gaussian whose covariance matrix matches that of the particle cloud (black). (a) t = 5 and t increases by 5 units for each panel. The axis labels, "x" and "z" have been removed to make larger sub-panels. 

since wave-induced flow above (below) the deformed pycnocline is in (against) the direction of wave propagation. As the lower panels 1-4 of Fig. 1 show the primary effect of the wave-induced currents is to greatly enhance the horizontal shear. As the particle cloud passes through the wave, it is initially stretched in both the horizontal and vertical directions. However, as the particle cloud passes beyond the wave crest (x = 0) it begins to be compressed in the vertical direction. Thus particles which find themselves below the pycnocline may begin their upward motion while particles which find themselves above the pycnocline are still moving downward. This leads to the tilde-shaped cloud labeled as 3 in Fig. 1. 

Figure 4 shows further details of this process, by identifying the outliers of the particle cloud after passage through the wave and tracing them backward in time. Outliers found farthest upstream are marked white, while those farthest downstream are marked black. It can be seen that the departure from the unperturbed location is quite significant (leftmost cloud of particles and solid circle). Particles that are farthest upstream (downstream) are found above (below) the center of the pycnocline. For early times (first cloud from the right) it can be seen that the particles that end up farthest upstream are those that have been perturbed to lie above the center of the pycnocline before they have moved into the wave. When they do move into the wave, both clouds of particles are advected downward to roughly the same degree. However the particles lying above (below) the center of the pycnocline are also advected in (against) the direction of wave propagation. This effects the time spent in the wave. Indeed, for the third cloud from the right, the outliers found below the pycnocline are moving up the rear face of the wave, while the outliers found above the pycnocline are at approximately the same height, two units of length upstream. 

The above described results are also reflected in the statistical properties of the particle distribution during, and after, passage through the wave. The standard deviation of the vertical position is non-monotonic (not shown) in time, with a local minimum occurring roughly at a time corresponding to cloud (3) in Fig. 1. Note, however that the cross-correlation between vertical and horizontal positions increases rapidly during passage through the wave indicating that the probability distribution is no longer aligned with the x and z axes. In Fig. 5 we show the particle clouds with the mean position subtracted off during passage through the wave (white), along with the quadratic form (an ellipse) corresponding to the covariance matrix of each particle cloud. Panels in the same row have the same axes. It is particularly interesting how the ellipses' axes are not aligned with the x and z axes after the particles pass completely through the wave (panels g-i). Were the particle distribution bivariate Gaussian, all the statistical information to describe it would be contained in the means and the covariance matrix. The black clouds of particles in each panel correspond to 10 000 draws from a bivariate Gaussian distribution matching the covariance matrix of the particle cloud at each particular time. It can be seen that from panel (d) onward the particle cloud is clearly non Gaussian. Indeed standard statistical tests (Mardia, 1970, Trujillo-Ortiz and Hernandez-Walls, 2003, the latter modified by D. Graham, personal communication, 2010) indicate that from panel (c) onward, the particle distribution is non Gaussian using either a skewness (third moment) or kurtosis (fourth moment) statistic, to a significance level smaller than 10 -4. Figures 1, 4 and 5 were all produced with the same ISW, as well as standard deviation and correlation time of red noise. In order to ensure that the presented results are broadly representative of parameter space we have varied the standard deviation of the noise, the noise memory and ISW amplitude. Figure 6 shows particle clouds as they pass the wave crest (ac) and as they leave the wave (d-f). The first column (a and d) corresponds to Fig. 1, the second to a 15 % reduction in wave amplitude (33 % reduction in wave available potential energy), and the third to a 55 % reduction in wave amplitude (83 % reduction in available potential energy). It can be seen that all the cases yield qualitatively unchanged behaviour, especially surprising in the case of the smallest wave. In particular, the particle distribution after passage through the wave has a sharp cutoff on the downstream side. Since exact ISWs are computed by a priori specifying the available potential energy, and not the wave amplitude, a very broad range of wave available potential energy values will yield similar particle behaviour. 

Figure 7 repeats Fig. 6 for changes to the stochastic parameters for the particle. The first column (a and e) corresponds to Fig. 1, the second (b and f) to noise with a tenfold increase in correlation time accompanied by a tenfold decrease of standard deviation, the third to a tenfold increase in correlation time with standard deviation (one fifth of the value in the first column) chosen by trial and error so that the particle cloud matches the first column, and the fourth to a one hundredfold increase in correlation time and a tenfold decrease in standard deviation. As was the case with varying wave amplitude, the figure shows that, by and large, the qualitative nature of the results is consistent over a very broad region of parameter space. Indeed, from the third column we note that even very different noise properties can yield essentially identical particle distributions. There are some differences evident in panels (b) and (f), and these indicate that for small standard deviations even strongly correlated noise does not yield a great deal of spreading of the particle cloud. Were we to allow the particle cloud to continue evolving, its spread in the vertical would continue unchecked, eventually tending to an unrealistic uniform distribution in the vertical. Systematic swimming by the plankton can modify this situation, and we discuss this case next.
section_name: Fig. 5.

section_content: Zooplankton exhibit a wide variety of swimming behaviour, going well beyond the purely random motion discussed above. Here we discuss the coupling of wave-induced currents, random motion and the simplest behaviour from the literature: return to a constant light (or photic) level (depth), modeled by modifying the vertical motion Eq. ( 2 ):

where w s represents the plankton's desire to return to its initial light level. For phytoplankton this swimming behaviour is an expression of its desire to remain in the region suitable for photosynthesis, while for zooplankton it is an expression of the organism's desire to remain at a depth optimal for zooplankton grazing and predator avoidance. 

If we consider purely vertical motion with no wave induced motion and take the white noise limit we can write down the Fokker-Planck equation (Gardiner, 1990) for the problem (subscripts denote partial derivatives)

or in steady state, 0 = -w s P z + σ 2 P zz /2. If we consider the simplest case of piecewise constant swimming behaviour, w s = w 0 -2H (zz l )w 0 where w 0 > 0 and z l is the vertical position corresponding to the desired light level we can find the steady state probability distribution,

where ζ = zz l and the decay length scale is given by ζ 0 = w 0 σ 2 /2. P 0 is a normalization constant. Clearly an increased swimming speed leads to a more localized steady state distribution, while an increase in the standard deviation of the random motion leads to a broader distribution. 

We extended the point release experiments from the previous subsection by allowing the ensemble to settle into a vertical equilibrium before being advected into the wave. We found that, provided the deterministic component of swimming (i. e. to the original light level) was below five times of the standard deviation of the random component of swimming, the results largely matched those of point release experiments, with an exact match as the speed of the deterministic component of the swimming tended to zero. More importantly, the return to light level swimming allows for the simulation of a more realistic situation in which the entire region upstream of the wave consists of a distribution of particles that are uniformly distributed in the horizontal while having a steady state distribution in the vertical that is a balance between random and systematic swimming. Such behaviour provides further symmetry breaking in the problem. Consider a particle initially at the mid-point of the pycnocline. As it begins to descend down the wave front it descends below its desired light level and the swimming behaviour causes it to swim upward. Once it passes the wave crest the wave-induced currents begin to force it upward to an equilibrium position that is above the mid-point of the pycnocline. Thus the particle overshoots its equilibrium level at the rear of the wave and only subsequently slowly returns to its equilibrium level. 

In Fig. 8 we present the results at t = 50, of a simulation with 50 000 particles initially in a vertical steady state, and uniformly distributed in the horizontal upstream of the wave (i. e. not the point release of Fig. 1 ). Panel a shows two particle clouds, the white has an equilibrium level at z l = -0. 25, or 5 % of the water column below the pycnocline center, while the black has an equilibrium light z l = -0. 15 or 5 % of the water column above the pycnocline center. It can be seen that the white cloud has moved farther through the wave. This is due to the fact that the black cloud lies completely in the region of wave-induced horizontal currents directed in the direction of wave propagation. The white cloud, on the other hand, begins in the region of wave-induced horizontal currents directed against the direction of wave propagation. However, as the white cloud is advected down the upstream face of the wave, the particles begin to swim upward toward their equilibrium level. They cross the pycnocline between x = -1 and x = 1, are partially advected in the direction of wave propagation, exit the wave above their equilibrium height, and subsequently slowly drift downward, reaching equilibrium some distance behind the wave. The black cloud experiences a similar, though smaller overshoot. Panel b shows the particle distribution as a function of x. It can be seen that the white cloud (dashed line) has an increased particle number density in the region -1. 25 < x < -0. 75, or at the rear of the wave, while the black cloud has an even larger increase in particle number density in the region -0. 1 < x < 0. 2. Both of these regions are well behind the maximum near surface convergence zone (near x = 2) where completely passive Lagrangian particle would be expected to have their maximum number density. 

In Fig. 9 we compare the overshoot region for the two cases. Panel a indicates that the mean overshoot is about 0. 035 for the particle cloud with z l = -0. 15, and 0. 055 for the particle cloud with z l = -0. 25. 

In nature, ISWs are often observed in rank-ordered groups, or packets (Helfrich and Melville, 2006). Provided that the inter-wave spacing allows the particles to return to a statistical steady state, near their equilibrium level z l, the above results should be observed for each wave in the packet, with the only difference being the individual wave amplitude.
section_name: Effects of systematic swimming behaviour

section_content: We have considered the motion of stochastically perturbed, Lagrangian particles through a passing finite amplitude, exact internal solitary wave. For deterministic motion, the primary effect of the wave is to enhance horizontal spreading of an ensemble of particles via a height dependent, waveinduced transport, and to induce increased particle number densities in the wave-induced convergence regions, e. g. near the surface at the front of a wave of depression, (Lamb, 1997). This is due to the fact that internal solitary waves are long waves, with horizontal wave-induced currents that are considerably stronger than vertical wave-induced currents. Moreover, horizontal wave-induced currents are directed in opposite directions above and below the main pycnocline, while vertical wave-induced currents are antisymmetric across the vertical line passing through the wave crest. The anti-symmetry of the vertical currents means that for particles without any random motion, a particle must return to its initial height after passing through the wave. 

The combination of stochastic and wave-induced motion leads to an asymmetric spreading of the particle ensemble. This asymmetry leads to the observation that the variance of the vertical particle position is a non-monotonic function of time, which is in turn an expression of the cross-correlation between vertical and horizontal particle position during, and after, passage through the wave. Indeed the particle distribution is non Gaussian during passage through the wave and remains so after passage through the wave. In particular, a marked skewness in the horizontal particle position is evident. For an ensemble of particles released upstream of the wave and at the pycnocline mid-depth, the statistical behaviour of the particles may be understood as being due to differential transport of particles that move above and below the pycnocline. Particle variance increases faster than that of purely stochastic particles (i. e. particles with only random motion) as the particles enter the wave, and decreases as the mean particle position reaches the wave crest. Interestingly, the convergence of the particle ensemble as it passes over the downstream half of the wave is sufficient to reduce the variance of the particle cloud after all particles have passed through the wave to values well below those of purely stochastic particles. These results were robust over a broad range of wave amplitudes, noise correlation times and standard deviations. 

We subsequently extended these results to particles that had an organized swimming component to their motion. While many types of swimming behaviour are possible, we chose to focus on return to initial light level behaviour since it is the simplest situation that allows for a steady state particle distribution in the vertical. For a uniform horizontal distribution upstream of the wave that was in steady state in the vertical we found that the wave led to an overshoot of the equilibrium particle height, and a consistent increase in the particle number density over the rear half of the wave. 

As a whole, the results suggest that irregular motion during passage through a wave does not yield any clear biological advantage to the small particles (zooplankton) that are generally found in, or above the pycnocline in the nearsurface region. However, for larger animals (e. g. plankton eating fish) the particle distribution shown in Figs. 8 and 9 could be exploited. The predator would sense the downward currents on the upstream side of the wave and swim so as to pass through the center of the cloud of plankton as it overshoots its equilibrium light level over the downstream face of the wave. Given that internal waves, often in rank-ordered groups, are systematically forced by tidal flow over the continental shelf (Helfrich and Melville, 2006) as well as other topographic features, predators would have an opportunity to exploit internal wave-induced currents to increase the efficiency of their feeding. For semi-diurnal tides, for example, this opportunity would occur at least twice a day. 

Future work should consider a wider variety of swimming behaviour, as well as more biologically relevant stochastic behaviour. Preliminary experiments with the "run and tumble" model (Bearon, 2007) suggest that many of the conclusions found for red noise stochastic motion are quite robust. The swimming of fish through internal waves could also be pursued using more realistic models, for example immersed boundary methods, though this would likely require high resolution, three dimensional simulations.
section_name: Discussion

section_content: Acknowledgements. This work was supported by the Natural Sciences and Engineering Research Council of Canada. The constructive criticism of two anonymous referees greatly improved the presentation. 

Edited by: A. M. Mancho Reviewed by: two anonymous referees