2.1. Individual swimmer wake model

This section introduces an analytical model to compute the flow field generated by many individual wakes in close proximity while conserving mass and momentum. This method is inspired by the approach adopted by Zong & Porté-Agel (Reference Zong and Porté-Agel2020) to superpose wind turbine wakes. Unlike previous formulations, which prescribe a drag coefficient and calculate momentum deficits, the present formulation prescribes the net force generated by the swimmers and calculates momentum excess. Importantly, this formulation did not assume a priori that the convective velocity would trend towards a plateau.

We assume that a vertical swimmer generates a downstream wake defined in the swimmer-fixed frame $u_w(x,y,z)$ to generate a net force $F_z$ that counteracts negative buoyancy and thus maintains a constant swimming speed $u_0$ through a fluid with constant density $\rho$. These assumptions allow for the simplification of the integral form of the momentum equation,

(2.1)\begin{equation} F_z = \rho \iint_{wake} {u_w(x,y,z)} (u_w(x,y,z)-u_0) \,{\rm d}\kern0.7pt x\,{\rm d}y .\end{equation}

By introducing the wake velocity surplus, $u_s = u_w-u_0$, and substituting this definition into (2.1) we obtain the following:

(2.2)\begin{equation} F_z = \rho \iint_{wake} {u_w(x,y,z)} u_s(x,y,z) \,{\rm d}\kern0.7pt x\,{\rm d}y. \end{equation}

We introduce an effective wake convection velocity, $u_c(z)$, which varies with downstream distance from the swimmer, but is constant in the spanwise directions. Consequently, the net vertical force can be rewritten as

(2.3)\begin{equation} F_z = \rho {u_c(z)} \iint_{wake} u_s(x,y,z) \,{\rm d}\kern0.7pt x\,{\rm d}y. \end{equation}

The wake convection velocity effectively represents the average speed at which the local velocity surplus is advected in the wake of the swimmer. To derive a mathematical expression for $u_c$, we substitute (2.3) into (2.2) to get

(2.4)\begin{equation} {u_c(z)} = \frac{\displaystyle\iint_{wake} {u_w(x,y,z)} u_s(x,y,z) \,{\rm d}\kern0.7pt x\,{\rm d}y}{\displaystyle\iint_{wake} u_s(x,y,z) \,{\rm d}\kern0.7pt x\,{\rm d}y}. \end{equation}

The numerical evaluation of (2.4) is described in § 3.4.1.

2.2. Wake superposition

To calculate the flow field at the aggregate scale, we define $U_{\infty }$ as the swimming speed of all organisms in the volume, or the free stream velocity in a swimmer-fixed frame. Furthermore, we introduce $U_w(x,y,z)$ as the global flow field generated by the swimmers. Lastly, $U_s$ is defined as the velocity surplus generated by the swimmers expressed as $U_s(x,y,z) = U_w(x,y,z) - U_{\infty }$. Following a procedure analogous to that in § 2.1, the effective convection velocity of the combined wakes is given by the following:

(2.5)\begin{equation} {U_c(z)} = \frac{\displaystyle\iint {U_w(x,y,z)} U_s(x,y,z) \,{\rm d}\kern0.7pt x\,{\rm d}y}{\displaystyle\iint U_s(x,y,z) \,{\rm d}\kern0.7pt x\,{\rm d}y}. \end{equation}

The force exerted by the $i$th swimmer in the streamwise direction is denoted $F_z^i$. To conserve momentum in the wake, we require

(2.6)\begin{equation} \sum_{i} F_z^i = \rho {U_c(z)} \iint U_s(x,y,z) \,{\rm d}\kern0.7pt x\,{\rm d}y. \end{equation}

Substituting the left-hand side of (2.6) with (2.3) yields the following:

(2.7)\begin{equation} \sum_{i} \rho {u_c^i(z)} \iint u_s^i(x,y,z) \,{\rm d}\kern0.7pt x\,{\rm d}y = \rho U_c(z) \iint U_s(x,y,z) \,{\rm d}\kern0.7pt x\,{\rm d}y. \end{equation}

The velocity experienced by the $i$th organism is denoted $u_0^i$ and defined as $U_w(x^i,y^i,z^i)$ based on upstream swimmers. The wake velocity induced by the $i$th organism is $u_w^i$, and the wake velocity surplus for the $i$th organism, $u_s^i$, is expressed as $u_w^i-u_0^i$. The rearrangement of these terms and the subsequent application of the analysis across the entire volume lead to the derivation of an expression for the global wake surplus,

(2.8)\begin{equation} U_s(x,y,z) = \sum_{i}\frac{u_c^i(z) }{U_c(z)} u_s^i(x,y,z). \end{equation}

In light of (2.5), which describes $U_c$ as a function of $U_s$ and (2.8), which characterizes $U_s$ as a function of $U_c$, an iterative methodology is used to solve for $U_s$ and $U_c$. The procedure begins with the assumption $U_c=U_\infty$, where $U_\infty$ denotes the velocity of the free stream. This is an underestimate, as $U_c$ will increase from the free stream velocity with added momentum provided from the swimmers. Thus, in the first iteration, (2.8) is used to evaluate $U_s$, and will result in an overestimate since its value is inversely related to that of $U_c$. As the iterative process continues, this overestimate of $U_s$ is used in (2.5) to refine the calculation of $U_c$, increasing the estimate from the initial guess. In this way $U_c$ will continue increasing from the initial guess and $U_s$ will continue decreasing until the value of $U_c$ converges, satisfying condition $|U_c-U_c^*|/U_c^* \leq \epsilon$, where $U_c$ is calculated from the preceding iteration, $U_c^*$ is calculated from the ongoing iteration and $\epsilon =0.01$. The iterative development of local and global estimated convective velocities captures inherent nonlinearity and ensures the conservation of momentum in the establishment of the final 3-D flow field.

3.1. Individual swimmer response to light stimulus

The response of brine shrimp to different light intensities was investigated in a controlled environment. A 1.2 m high tank with a cross-section of $0.3~{\rm m}\times 0.3$ m (figure 1) was filled with 35 parts per thousand of salt water using Instant Ocean Sea Salt (Spectrum Brands). To reduce the influence of swimmer wakes on one another, the tank was populated with less than 1500 swimmers, or 0.015 animals per cm$^{3}$. All experiments were carried out within 24 h of animal acquisition.

Figure 1.

Schematic of experimental protocol for characterization of phototactic response of brine shrimp. (a) Brine shrimp (i) were gathered at the bottom of a 1.2 m tall tank using a flashlight positioned at the base. (b) To induce vertical migration, the bottom flashlight was turned off, and the top flashlight was turned on. The light intensity of the top flashlight was varied using neutral density filters. Recording was manually initiated once the swimmers entered the field of view of the high-speed camera (ii). (c) Example frames from the $x$–$z$ plane captured during vertical migrations under different light intensities, adjusted with neutral density filters. All trials were conducted with an infrared tank illuminator, which was used exclusively in the 0 lux condition. Three trials were performed for each light intensity: 0, 800, 1500, 2300 and 4000 lux.

To ensure consistency between trials, the animals were gathered at the bottom of the tank using a flashlight, and a minimum settling time of 30 min was allowed between each trial. Initiating a vertical migration involved turning off the flashlight at the tank's bottom and activating a target flashlight (PeakPlus LFX1000, 1000 lumens) positioned above the tank. The light intensity of this upper flashlight was adjusted using three different neutral density filters: 1/2, 1/4 and 1/8 transmittance (Neewer 52 mm ND Filter Kit). A light intensity meter (TEKCOPLUS Lux Meter with Data Logging) was used to measure the illumination at the bottom of the tank for each filter.

A high-speed camera (Edgertronic SC1) was set up with a $20~{\rm cm}\times 25$ cm ($1024\times 1280$ pixel) field of view, 60 cm above the bottom of the tank. For each test, a recording was manually triggered once the first swimmer entered the camera field of view and captured for 30 s at 40 frames per second (f.p.s.). Four trials were carried out with each of the three filters (800, 1500, 2300 lumens per square metre (lux)), without a filter present (4000 lux) and without the target light (0 lux). An infrared (850 nm) light was used to illuminate the tank and collect control data for the case in which no visible illumination was present. For consistency across all tests, this infrared illumination remained on for all tests.

3.2. Individual swimmer response to background flow

To simulate the vertical flows induced during collective vertical migration of brine shrimp, water was drained from a 2.4 m tall tank with a cross-section of $0.5~{\rm m} \times 0.5$ m, producing uniform flows in the range of 0.05–0.5 cm s$^{-1}$ (figure 2). These flow speeds correspond to those observed in vertical migrations of brine shrimp, with animal number densities between 100 000 and 600 000 animals per cubic metre (Houghton & Dabiri Reference Houghton and Dabiri2019). The uniformity produced by this set-up reflects the uniform jet produced by steadily moving dilute swarms in discrete swimmer simulations (Ouillon et al. Reference Ouillon, Houghton, Dabiri and Meiburg2020).

Figure 2.

Schematic of the experimental protocol for characterizing the flow response of brine shrimp. (a) Brine shrimp (i) were initially gathered at the bottom of a 1.2 m tall tank using a flashlight positioned at the base. (b) Brine shrimp (i) were initially gathered at the bottom of a 1.2 m tall tank using a flashlight positioned at the base. (c) Once the swimmers entered the field of view of the high-speed camera (ii), a flow valve was opened, introducing bulk flow in the opposite direction of the swimmers’ motion. Recording was manually initiated at this point. The flow rate was controlled using a system of two flow valves and a flow meter arranged in series. Three trials were conducted for each target flow speed: 0, 0.07, 0.14, 0.21 and 0.3 cm s$^{-1}$.

Before testing, the tank was filled with 10 $\mathrm {\mu }$m silver-coated glass spheres (CONDUCT-O-FIL, Potters Industries, Inc.) to facilitate imaging of the flow field with a laser sheet. To confirm the quiescence of the tank, particle image velocimetry (PIV) was employed after introducing the animals with a 15 ml centrifuge tube. The tank was considered quiescent when the maximum time-averaged streamwise velocity was below 0.02 cm s$^{-1}$. Flow rate control was achieved using two series connected flow valves (1 in. NPT PVC Ball Valve), one for flow control and one for shut-off, and an inline flow meter (FLOMEC Flow meter/Totalizer 5–50 gpm).

Once the tank was confirmed to be quiescent with PIV, a migration was induced with the same procedure explained in § 3.1. A high-speed camera (Edgertronic SC1) was set up with a field of view of $21~{\rm cm} \times 26$ cm ($1024 \times 1280$ pixel), 90 cm up from the bottom of the tank. Once the first swimmer entered the camera field of view, the shut-off valve was manually opened to initiate the flow, and the camera was manually triggered to record for 30 s at 15 f.p.s. Three trials were carried out for each of the five target speeds: 0, 0.07, 0.14, 0.21 and 0.3 cm s$^{-1}$. The trials were carried out on different days using different animals, considering the limited number of trials achievable with the volume of the tank.

To identify potential sources of measurement uncertainty, two significant factors were addressed. First, before the onset of the flow, a streamwise velocity variation of 0.02 cm s$^{-1}$ was allowed. Second, the flow speed was manually set using a ball valve and changed during each test due to variations in the height of the water column during draining. To address these uncertainties, an additional camera recorded flow rates displayed on the inline flow meter during each test. Both of these sources of variability were accounted for in the error bars of all flow measurements.

The captured videos of the vertical migration of brine shrimp were analysed using FIJI (Schindelin et al. Reference Schindelin2012) and the wrMTrck plugin (Husson et al. Reference Husson, Costa, Schmitt and Gottschalk2018). The resulting swimming trajectories were fitted in MATLAB with a smoothing spline algorithm, which minimizes a combination of squared residuals and curvature penalties utilizing cubic smoothing spline interpolation to fit a curve to the provided data points. A smoothing parameter of 0.95 (figure 3) was selected to prioritize the reduction of oscillations, effectively smoothing out the fitted curve while preserving the overall trend of the data.

Figure 3.

Plot of brine shrimp swimming trajectories generated with ImageJ's wrMTrck plugin over a 30 s interval on top of the final frame in the image sequence. Colour transitions from blue to yellow represent progression in time, illustrating the trajectory of each swimmer within the tank. Brine shrimp present in the final frame can be identified as the white silhouettes at the end of the trajectories. Gradations in the background shading are due to the illumination used to induce phototaxis.

3.3. Characterizing collective swimming

To reconstruct the 3-D swimming trajectories of brine shrimp during induced vertical migration, a 3-D particle tracking velocimetry (PTV) method was used (figure 4), using scanning optics and a single high-speed camera (Photron FASTCAM SA-Z). A 671 nm continuous wave laser (5 W Laserglow LRS0671 DPSS Laser System) was directed through a condenser lens (370 mm back focal length) and a sheet-forming glass rod by a mirror to ensure parallel beams. The distance between the laser plane and the high-speed camera was adjusted with a galvanometer (Thorlabs GVS211/M) controlled by a voltage signal from an arbitrary function generator (Tektronix AFG3011C). This experiment was validated on a smaller scanning volume in the same facility with the same equipment by Fu et al. (Reference Fu, Houghton and Dabiri2021).

Figure 4.

Schematic of the scanning system, modified from Fu et al. (Reference Fu, Houghton and Dabiri2021). The laser beam was directed at a mirror whose angle was controlled by a galvanometer. Following reflection by the mirror, the beam was passed through a condenser lens and glass rod in order to generate a laser sheet parallel to the camera field of view.

Each laser sheet sweep covered 6.6 cm of tank depth and took 0.1 s to complete. The high-speed camera captured 300 two-dimensional (2-D) $22~{\rm cm} \times 22$ cm (1024 x 1024 pixel) slices during this period. The scanned volume was centred on the tank cross-section and positioned 0.9 m from the tank floor. The swimmers move at approximately 1 cm s$^{-1}$; therefore, during the length of a scan (0.1 s), the swimmers will have moved approximately 0.1 cm. Given their body size of 1 cm, we effectively treated each scan as a still frame for the purposes of the current analysis. In addition, this combination of scanning rate and camera frames per second results in approximately 48 image sheets per centimetre in the scanning direction ($y$) and 46 pixels per centimetre in the camera plane ($x$ and $z$). This level of resolution for detecting swimmers 1 cm long results in well-formed and easily identifiable swimmers, as shown in figure 5(a). A movie of 3-D reconstructed swimmers for 6 s is in the supplementary materials available at https://doi.org/10.1017/jfm.2024.1102.

Figure 5.

Volume reconstruction and animal tracking during an induced vertical migration in the positive $z$ direction. (a) A 3-D scanned reconstruction of 100 animal bodies. A close-up of an individual animal shown for reference. (b) Six-second 3-D swimming trajectories, with colour transition from blue to yellow represent progression in time, illustrating the trajectory of each swimmer within the tank.

The experiments were carried out in the tank described in § 3.1. The brine shrimp were added to the tank in densely packed 0.25 teaspoon increments (approximately 125 swimmers) and three vertical migrations were induced as explained in § 3.1 for each increment of swimmers added. Throughout vertical migration, the centre of the tank was scanned for 6 s every 50 s, totalling 7.5 min, to assess configuration changes over time.

A 3-D volume was constructed from 295 2-D slices for each laser sweep (figure 5a). A custom MATLAB script was used to segment the 3-D volume and identify centroids. The volumetric data was median and Gaussian filtered to reduce noise and enhance object visibility. The filtered data was then binarized and morphological operations were applied, including opening with a spherical structuring element and hole filling, to refine the binary mask. This preprocessing ensured a clean and noise-reduced dataset for object detection. Through object detection, properties such as centroid, principal axis and volume were identified for each connected component within the volume. Subsequently, a forward and backward nearest neighbour search in time was applied to centroid locations to identify and label swimming trajectories (figure 5b).

3.4. Modelling assimilation

3.4.1. Wake profile models

Two models for the individual wake structure were studied to explore the impact of local flow geometry on the aggregation-scale flow. The local flow was defined as a function of the radial distance, $r = \sqrt {x^2 + y^2}$, and the characteristic width of the wake, $\sigma (z)$, at each value of $z$. First, a Gaussian model was implemented, consistent with the wake models previously used for wind turbine modelling,

(3.1)\begin{equation} \xi_{gaussian}\left(\frac{r}{\sigma(z)}\right)= {\rm e}^{-{r^2}/{2\sigma(z)^2}}. \end{equation}

Second, a Ricker wavelet model was used to represent a local flow both in the direction of swimming and in the opposite direction of swimming,

(3.2)\begin{equation} \xi_{wavelet}\left(\frac{r}{\sigma(z)}\right)={-}\left(1 - \frac{r^2}{2\sigma(z)^2}\right) {\rm e}^{-{r^2}/{3\sigma(z)^2}}. \end{equation}

These models will be referred to as the Gaussian and wavelet models, respectively. The Gaussian model, commonly employed in wind turbine modelling (Zong & Porté-Agel Reference Zong and Porté-Agel2020), represents the flow behind the swimmer as a single downward jet. As evidenced by the qualitative match between figure 6(a) and the schematic in figure 6(c), the Gaussian wake function effectively captures the flow characteristics behind a swimmer utilizing a single main propulsor, generating a distinct single-lobed jet. The wavelet model, derived from the modified Ricker wavelet, is based on the second derivative of a Gaussian function, with an adjusted prefactor in the exponent denominator in order to create a function with a non-zero integral. As illustrated by the quantitative similarity between figure 6(b) and the schematic in figure 6(d), the wavelet wake function captures the flow distribution generated by a swimmer with two sets of propulsive appendages, producing a double-lobed jet and a region of backflow immediately behind the swimmer due to the drag created by the body.

Figure 6.

Comparative analysis of flow fields generated by different swimmer types and corresponding wake models. (a) The PIV results showing the flow field generated by a single free-swimming Pacific krill, adapted from Catton et al. (Reference Catton, Webster, Kawaguchi and Yen2011) with the permission of the Journal of Experimental Biology. The colour map represents the velocity magnitude, with arrows indicating flow direction. (b) The PIV results for a brine shrimp, highlighting the flow characteristics generated by its swimming motion, adapted from Wilhelmus & Dabiri (Reference Wilhelmus and Dabiri2014), with the permission of AIP Publishing. (c) A schematic representation of how a Gaussian distribution qualitatively captures the flow behind Pacific krill, which utilizes a single primary propulsor. (d) A schematic representation of a wavelet model that qualitatively captures the flow behind a brine shrimp, characterized by two sets of propulsors and a drag region immediately behind the body that induces backflow. The panels illustrate the differences in flow structures arising from distinct swimming mechanisms and motivate the comparison between a Gaussian and wavelet wake superposition model.

We modelled the spatial evolution of the propulsive jet as a self-similar, axisymmetric jet,

(3.3)\begin{equation} u_s = u_0 c(z)\xi\left(\frac{r}{\sigma(z)}\right), \end{equation}

where $c(z)$ is a scaling factor and the shape of the wake is determined by $\xi (r/\sigma (z))$. An expression for $c(z)$ that conserves momentum by definition was derived by prescribing $u_0$, $\xi (r/\sigma (z))$, and $F_z$ (detailed steps provided in the Appendix (A)),

(3.4)\begin{gather} u_{s,gaussian}(x,y,z)= u_0 \left(1 - \sqrt{1 + \frac{F_{z}}{{\rm \pi} \rho \sigma(z)^2 u_0^2}} \right) {\rm e}^{-{r^2}/{2\sigma(z)^2}}, \end{gather}

(3.5)\begin{gather}u_{s,wavelet}(x,y,z)= u_0 \left( \frac{4}{5}-\frac{4}{5}\sqrt{1 + \frac{5 F_{z}}{12 {\rm \pi}\rho \sigma(z)^2 u_0^2}}\right) \left({-}1 \left(- \frac{r^2}{2\sigma(z)^2}\right) {\rm e}^{-{r^2}/{3\sigma(z)^2}}\right). \end{gather}

Substituting (3.5) and (3.4) into (2.4) and integrating in the streamwise direction over a circular cross-section with infinite radius, the expressions for $u_c$ is obtained

(3.6)\begin{gather} u_{c, gaussian}(z) = \frac{u_0}{2} \left(1+ \sqrt{1+ \frac{F_z}{{\rm \pi} \rho \sigma(z)^2 u_0^2}}\right), \end{gather}

(3.7)\begin{gather}u_{c, wavelet}(z) = \frac{u_0}{2} \left(1+ \sqrt{1+ \frac{5 F_z}{12 {\rm \pi}\rho \sigma(z)^2 u_0^2}}\right). \end{gather}

A swimmer moving vertically at constant velocity must overcome the negative buoyancy that arises from the swimmer having a greater density ($\rho_s$) than sea water. The balance of force on the swimmer is expressed as follows:

(3.8)\begin{equation} F_{thrust}= g V_s (\rho_s-\rho). \end{equation}

The thrust is introduced over some distance and not at an exact point in the flow. Therefore, we amended this expression to

(3.9)\begin{equation} F_{thrust}(z)= g V_s (\rho_s-\rho) \frac{1 + \operatorname{{\rm erf}} \left(\dfrac{z}{L_s}\right)}{2}, \end{equation}

for a more gradual development of the wake, where $L_s$ is the length of the swimmer's body length (BL). Similarly, an empirical model for the effective diameter of the wake as a function of the streamwise distance from the swimmer was used to capture the wake expansion,

(3.10)\begin{equation} \sigma (z)= 0.25 + 0.25 \log(1+{\rm e}^{(z-1.5)/L_s}). \end{equation}

The variables explicitly defined for the numerical implementation of this wake model are listed in table 1. These values were approximated to be of the order of observations made during laboratory experiments using brine shrimp. We normalize all length measurements by the BL of a swimmer, $L_s$.

Table 1.

Variables used in the wake superposition model for induced flow in vertical migration.

3.4.2. Collective flow field calculations

Model flow fields for various swimmer configurations were calculated to identify the impact of aggregate characteristics on induced flow. First, changes due to group length were examined. The length of the group was increased in each test, while animal number density and width remained constant (table 2a). To maintain constant animal number density and width of the group, the number of swimmers increased linearly with increasing length of the group. Second, to test the impact of animal number density on the resulting flow, the number of swimmers and the length of the group were kept constant while increasing the cross-sectional area in the spanwise dimensions (table 2b). This resulted in an animal number density that decreased with increasing width as $1/W^2$, where $W$ is the width of the group. For each calculation with a selected set of parameters, three iterations of swimmers were placed randomly with these specifications while maintaining a minimum nearest neighbour distance of one BL.

Table 2.

Parameters to be examined are the number of swimmers in the group, $N$, the length of the group, $L$, and the width of the group, $W$. Together, these three parameters result in a group metric that we refer to as the animal number density, measured in animals per BL$^3$ and calculated as $N/(W^2L)$. These parameters are used to examine the impact of changes in (a) group length, and (b) animal number density.