Numerical Simulations of Flow around Copepods: Challenges and Future Directions

Abstract

Copepods are small aquatic creatures which are abundant in oceans as a major food source for fish, thereby playing a vital role in marine ecology. Because of their role in the food chain, copepods have been subject to intense research through different perspectives from anatomy, form-function biology, to ecology. Numerical simulations can uniquely support such investigations by quantifying: (i) the force and flow generated by different parts of the body, thereby clarify the form-function relation of each part; (ii) the relation between the small-scale flow around animal and the large-scale (e.g., oceanic) flow of its surroundings; and (iii) the flow and its energetics, thereby answering ecological questions, particularly, the three major survival tasks, i.e., feeding, predator avoidance, and mate-finding. Nevertheless, such numerical simulations need to overcome challenges involving complex anatomic shape of copepods, multiple moving appendages, resolving different scales (appendage-, animal- to large-scale). The numerical methods capable of handling such problems and some recent simulations are reviewed. At the end, future developments necessary to simulate copepods from animal- to surrounding-scale are discussed.

1. Introduction

Copepods (from the Greek word for oar feet) are among the most diverse animals in the world with more than 14,000 described species. They can be found in almost any kind of aquatic environment, from subzero waters of Arctic to hot springs, from the top of Himalayas to 10,000 meters down the deep sea, in mud, subterranean groundwater, lakes, seas and oceans [1]. Some copepods are parasitic and live off other animals such as fish—e.g., salmon louse Lepeophtheirus salmonis Kroyer is a devastating pest for salmon farms [1]. They are among the most abundant animals on the planet and a major food source for fish. They play a vital role in the food chain and marine ecology which have made them the subject of intense research through different perspectives from anatomy, form-function biology, to ecology. In this manuscript, the numerical methods capable of supporting such research are reviewed and their limitations and strengths are discussed.

The anatomy of copepods is well-studied and a general body plan is consistent across myriad orders of these small crustaceans. Their size typically ranges between 0.5 to 5 mm. They have a segmented body with several appendages attached to it (see Figure 1). A copepod has a pair of first and a pair of second antennae, the first pair is long and the other is short, mandibular palps, two pairs of maxillas, and a pair of maxilliped. There are four or five pairs of swimming legs attached to the abdomen and the urosome with setae at the caudal area. They have one simple eye in the middle of the head (at least in the larval stage), which can only tell the difference between light and dark. Most of the sensory ability of copepods comes from the chemoreceptors on the mouthparts [2] and mechanoreceptors over appendages especially the first antennae [3]. Overall, simulating a copepod with an anatomically realistic shape and moving appendages poses a great challenge even to the most advanced numerical tools—see Section 2.

The behavior of copepods and the function of their appendages has been typically investigated through observations using high-speed cameras. Copepods feed on a wide variety of prey, ranging from algae of a few micrometers to metazooplankton and fish larvae [5]. They use their mouth appendages i.e., the antennae, maxillia and maxillipeds to create a feeding current toward the mandibular palps (jaws) [6,7,8]. The flow field created by the feeding copepod is laminar with a low Reynolds number based on length and swimming speed of order 1 to 10, where the viscous forces are dominant. Copepods spend most of their time swimming and feeding at this low Reynolds number environment. However, if a threat is detected they respond with rapid jumps to escape from harm’s way. Such escape maneuvers give copepods a much greater chance of survival than the zooplanktons which do not exhibit such predator-avoidance mechanism [9]. The Reynolds number during jumps may reach as high as several thousands, which may even transition to turbulence.

The appendages move differently during an escape maneuver from the normal cruising mode. When cruising, only the mouth appendages move to create a laminar feeding current while in the escape the power strokes of antennae and the legs are mainly used to create large acceleration. The power-stroke starts with the beating of the first antennae followed by multiple metachronal beating of the legs while other mouth appendages stay in the retracted position [10]. During the return stroke swimming legs minimize their surface area and move synchronously to the original position [10,11,12]. Such power-strokes can give copepods the incredible speed of 50 to 500 body-length per second during escape [12,13,14]. To understand how copepods achieve such high velocities, the hydrodynamic forces produced by the appendages movements should be determined.

A few experimental studies have been conducted to measure the hydrodynamic forces during the power and return strokes. Alcaraz and Strickler [11] examined the relation of forces and appendages movement for escape mode by measuring the spring force attached to the tethered copepod while filming the appendages movement from the side. They found the force to be in the thrust (forward) direction during the power stroke and in the drag (backward) direction during the return stroke. More recently, Lenz et al. [12] performed a similar study. They reported that the peaks on the force record corresponded to the power stroke of each leg and the stroke of forth and third pairs of legs produced the largest peak. The estimated power per muscle mass by Lenz et al. [12] was found to be particularly high relative to other organisms. In the above studies, the force record was the total force produced by all the appendages and it was not completely clear how much each appendage contributed to the total force. A major contribution of numerical simulations is that they specifically provide how much each appendage contributes to the total force and how the appendages movements affect the forces produced, thereby clarify the form-function relation of anatomical structures or appendages.

Another point of interest in copepod research is the flow field that is generated by the copepods and how copepods respond to the flow in their environment [15,16]. In fact, copepods respond to hydrodynamic signals within a few milliseconds [17]. They use their mechanoreceptors to feel the hydrodynamic disturbances in the flow [18]. It is believed that they can distinguish between prey, predator, or mate from the hydrodynamic signatures of the flow [19]. In addition, the hydrodyanmic disturbance that they create is used by fish to capture them. The flow field and the coherent vortical structures has been visualized in experiments with live copepods [14,19,20,21,22]. Such flow visualizations typically capture the flow over a 2D plane. 3D flow measurements are possible, e.g., using tomographic particle image velocimetry (tomo-PIV) [23], but they are more challenging and normally have a lower resolution than 2D PIV. Numerical simulations can complement such experiments by providing the detailed 3D flow around the copepods. Nevertheless, the numerical simulation of a copepod is a challenging undertaking due to its complex shape and the thin appendages and their rapid movements. In the next section, available methods for simulating flow around copepods with different levels of complexity are reviewed and their strengths and limitations are discussed.

2. Overview of Numerical Methods for Simulations of Copepods

The methods for simulating copepods can be classified into two main categories: Simulations that resolve the shape and motion of copepod (appendage-scale methods), and the simulations that ignore the exact shape of a copepod or the motion of its appendages and model their effect as a force field on the flow (force-field method). The former is better suited for investigating form-function relation of different anatomical features of a copepod, whereas the latter is useful for investigating the general flow features generated by copepod (copepod-scale), the interaction of copepod-scale flow and the large-scale (e.g., oceanic) flow, or their collective behavior. In what follows, these two categories are further discussed.

2.1. Force-Field Simulations

When the large- or copepod-scale flow is of interest, i.e., the length scales of surrounding flow are much larger than the appendages size, the copepod shape or its appendages are ignored but its effect on the flow is modeled as a force field—see the reviews [24,25] for more details and applications of this method. In fact, the background flow is governed by the incompressible Navier-Stokes (or Stokes depending on the Reynolds number of the flow) equations and a force is added at the locations where the copepods are present [26]. This method is similar to what has been used for simulation of small particles [27] or bacteria and other microorganisms (active fluids) in the flow [28],

The 3D incompressible, non-dimensional Navier-Stokes equations are as follows:

$$\begin{array}{cc}\hfill & \nabla .\mathbf{u}=0\hfill \\ \hfill & \frac{D\mathbf{u}}{Dt}=-\nabla p+\frac{1}{Re}{\nabla }^2\mathbf{u}+\mathbf{f}\hfill \end{array}$$

(1)

where $D/Dt=\partial /\partial t+\mathbf{u}.\nabla$, $\mathbf{u}$ is the velocity vector, p is pressure, $Re$ is the Reynolds number, and $\mathbf{f}$ is the force field (force per unit volume) exerted by the copepods onto the flow. The above equations require appropriate boundary conditions to form a well-posed problem. The boundary conditions are the no-slip condition at copepod locations in flow

$$\mathbf{u}={\mathbf{V}}_{\mathbf{swimming}}$$

(2)

where ${\mathbf{V}}_{\mathbf{swimming}}$ is the swimming velocity of the copepod.

This method has been used to simulate the flow at the copepod-scale. Jiang et al. [26] developed five models for the force $\mathbf{f}$ in Stokes flow for different types of steady swimming such as hovering, freely sinking, and upward/backward/forward swimming. The force field $\mathbf{f}$ was obtained through force balance by assuming the shape of the copepod body as a sphere and using Stokes drag formulas. More recently, Jiang and Kiørboe [29] used an impulsive stresslet to model the force field generated by a jumping copepod. In addition, Jiang and Kiørboe [30] simulated the copepod jumps by assuming an spheroidal shape for the body and applying a force field to account for the beating legs. The flow field created by self-propelled jumping copepod compared well with PIV measured flow fields and spatial decay rate [30].

The force field method can also be used to simulate the interaction of copepod with the large-scale flow, e.g., turbulence. In the pioneering work of Yamazaki et al. [31], direct numerical simulations of an isotropic turbulence were performed and copepods were simulated as passively drifting (inertialess) particles by the turbulent flow plus a random walk. Later, Squires and Yamazaki [32] simulated marine particles with inertia in a similar setup and showed the preferential concentration (rather than a homogeneous distribution) generated by such particles. Lewis and Pedley [33] used a probablity density function of predator and prey velocities, rather than direct numerical simulations, to model the contact of microorganisms in turbulent flow. Recently, Ardeshiri et al. [34] carried out direct numerical simulations of isotropic turbulence and simulated the motion of copepods as particle similar to Yamazaki et al. [31] but instead of the random walk they developed and used a Lagrangian model based on the recorded velocities of copepods in an experimental setting. In all of the above studies [31,32,34], a non-uniform body force was applied to the largest scales of motion in order to maintain a statistically stationary flow, but the effect of copepods on the turbulent flow was ignored and no force from the copepods were added to the right hand side of Equation (1). In fact, the copepods’ motion were determined by the flow but the existence of the copepod did not affect the flow, i.e., the copepods’ motion was coupled one-way to the flow. Nevertheless, it is known that the existence of particles in the turbulent flow affects the turbulent characteristics as well as particle dispersion, which requires two-way coupling [35,36]. In two-way coupling, the force exerted by the copepods on the flow needs to be added to the right hand side of Equation (1).

2.2. Appendage-Scale Simulations

For simulations of copepods with realistic body and appendages shapes, the Navier-Stokes Equation (1) are solved over a grid with appropriate no-slip conditions on the body (Equation (2)). Nevertheless, handling the motion of copepod appendages require special techniques because the boundary at which the no-slip conditions should be applied is moving. The numerical methods for moving boundary problems are classified into two main categories [37]: (a) boundary-conforming methods in which the grid moves the moving boundary and (b) non boundary-conforming methods in which the grid is fixed and the moving boundary moves over a fixed background mesh.

In boundary-conforming methods, also known as the arbitray Lagrangian-Eulerian (ALE) methods [38], the Navier-Stokes Equation (1) are modified to account for the grid motion. The ALE methods can efficiently keep a high-resolution near the moving boundaries. For large deformations, however, the grid can become highly skewed or stretched, i.e., the grid quality degrades, which may cause convergence issues for the numerical method. Remeshing may solve the grid quality, but it is computationally expensive and difficult.

Non boundary-conforming, i.e., fixed-grid, methods do not need any modification to the Navier-Stokes equations as the grid is fixed and does not move with the boundary. These methods are more suitable for large deformations because they do not create highly skewed grids [39]. However, additional complexities arise in terms of identifying the grid nodes adjacent to the moving boundaries and transferring the effects onto those nodes. There are different fixed-grid methods, e.g., the immersed boundary method [39], cut-cell methods [40], and fictitious domain method [41], among others. Out of these methods, only the immersed boundary has been used for simulations of copepods [4]. Consequently, only the immersed boundary is briefly described below.

The original immersed boundary, pioneered by Peskin [42,43,44], smeared the effect of boundaries over several grid nodes, which required additional resolution near the boundary. Since then, a number of sharp-interface methods have been developed [39]. In the method used for copepods [4], the sharp-interface is maintained by reconstructing the boundary conditions at the nodes that are exterior to, but adjacent to the immersed-boundary surface using a quadratic interpolation along the local normal to the boundary [45]. The background nodes at each time step are classified into fluid, wall, and immersed boundary nodes using an efficient ray-tracing algorithm [46]. The method uses curvilinear grids as the background grid to efficiently stretch the grid in regions of interest, i.e, the curvilinear-immersed boundary (CURVIB) method [47]. The method is fully parallelized using Message Passing Interface (MPI) and Portable, Extensible Toolkit for Scientific Computation (PETSc) [48]. It has been fully validated [46,49,50,51] and applied to a wide range of applications: vortex induced vibrations [52,53], aquatic locomotion [54,55,56,57,58,59,60,61,62,63], cardiovascular flows [64,65,66,67,68,69,70], sediment transport [71], large-eddy simulations and flow control [72,73,74,75], rheology [76,77], and copepods [4]. In the next section, some of the copepod simulations using this method are presented and discussed.

3. Results from Previous Simulations of Copepods

The earliest numerical simulations of copepods considered copepods as particles drifting by a turbulent flow [31,32]. In their pioneering work Yamazaki et al. [31] and Squire and Yamazaki [32] tracked the trajectories of 512 and 165,888 copepods (particles), respectively, drifting in an isotropic turbulent flow. The results showed that the contact was increased due to turbulence [31] and clustering can occur with densities of 10 to 60 times of the mean average population density [32]. Recently, Ardeshiri et al. [34] simulated 256,000 copepods with a Lagrangian model of jumping behavior in an isotropic turbulent flow and found clustering of jumping copepods. Furthermore, they used their simulations to estimate the contact rate of jumping copepods and found that it can be increased by a factor up to ${10}^2$ compared to that experienced by passively transported fluid tracers of the same size [78]. These simulations investigated the role of large-scale flow on the copepod motion, but ignored the copepod-scale flow generated by the copepods.

The first attempt to simulate copepod swimming numerically was reported in a series of papers by Jiang et al [26,79,80,81]. In their theoretical work, the force field approach (Section 2.1) was used by approximating the shape of the copepod as a sphere and modeling the force for different swimming conditions [26,29]. In their numerical simulations [79,81], a realistic body shape was used but the multiple moving appendages were neglected and their effect was accounted for via a distribution of body forces acting on the water. The first simulation of an anatomically realistic copepod with all of its swimming appendages was reported by Gilmanov and Sotiropoulos [45] who employed a sharp-interface Cartesian method (Section 2.2) to a model the swimming of a tethered copepod. Later, Borazjani et al. [4] built on that work by carrying out high-resolution simulations over a range of governing parameters with more biologically accurate kinematics prescribed from high-resolution experiments.

Borazjani et al. [4] used 3D numerical simulations of a tethered copepod to investigate the role of the first antennae in production of hydrodynamic force during hopping (Figure 1). The details of body shape and the kinematics of the appendages were prescribed from experimental observations using a dual digital holography setup with 200 frame per second high-speed camera to record the motion of hopping copepods [21]. By analyzing the high-resolution experimental holographic movies frame-by-frame, it was observed that the copepod antennae deform distinctively during the return stroke to the fully open position [4]. Consequently, two sets of simulations were conducted to examine the importance of the realistic motion of the antennae on the hydrodynamics of during hopping as observed from the top view in Figure 2 [4]: one treating each antenna as a rigid, oar-like structure moving symmetrically during power and return stroke; and the other using an asymmetrical motion with a deformed shape during the return stroke which was prescribed from experiments.

The force produced by each individual appendage was computed directly from the numerical simulations and its relative significance on the total hydrodynamic force was quantified [4]. The computed results: (1) show that for both cases the antennae are major contributors to the net thrust during hopping; and (2) clearly demonstrate the significant hydrodynamic benefit in terms of thrust enhancement and drag reduction due to biologically realistic, asymmetric antenna motion. In addition, the antennae were found to produce the largest drag- and thrust-type forces among all appendages regardless of whether the antennae were treated as rigid or deformable. This conclusion, however, could be partly affected by the fact that the leg kinematics in the present model does not exhibit the required degree of biological realism. Nevertheless, the large drag of the antennae provides a plausible explanation why copepods have been observed to use multiple leg strokes before re-deploying their antennae. Finally, the calculated net force time history over a swimming cycle explains why copepods swimming in the ocean move in bursts (jumps) during escape. The high thrust (produced by all appendages) during the power stroke, which rapidly accelerates the copepod, is followed by a large drag force (produced by the body and the appendages) during the return stroke, which decelerates the copepod, yielding an intermittent, burst-like behavior during the swimming cycle.

The numerical simulations also provide the 3D flow around the copepod. The antennae create complex vortical structures that can be viewed from the top (Figure 2). Figure 2A shows the antennae at the end of the its power stroke, which creates a strong shear layer at the tips of each antenna. Figure 2B,C shows the wake at the middle and end of the return stroke, respectively. It can be observed that the downstream wake structures have been destroyed by the rigid antennae while still moving downstream in the deformable case i.e., more thrust-type wake for the deformable antennae. The numerical simulations also elucidate for the first time the origin of the distinct toroidal vortices observed in flow visualization experiments in the wake of escaping copepods [19]. Visualizing the coherent structures by the iso-surfaces of q-criterion (defined as the difference between the norm of vorticity and strain-rate [82]; consequently, q > 0 are regions where the rotation rate dominates the strain rate, i.e., a vortex) in the wake of the copepod showed the three-dimensional structure of the toroidal vortices (Figure 3). It was found that the toroidal vortices are formed mainly from the leg strokes by forming two columnar vortices, which are attached to each other at the end forming an arch-like vortex loop [4].

4. Discussions and Future Directions

Numerical simulations have come a long way from simplified geometries to realistic anatomic ones. Numerical simulations can provide the 3D flow around copepods with realistic anatomies and appendages motion, which is challenging to measure experimentally, and provide hydrodynamic forces generated by individual appendages or body. In addition, one of the main benefits of numerical simulations is to test hypothetical cases which might not be possible in real world. Such numerical simulations can play an important role in investigating the form-function relations of biological features. In fact, Borazjani et al. [4] used numerical simulations to investigate the form and function of copepod antennae during hopping. This was possible because the force generated by each appendage can be computed separately from others in numerical simulations which is quite challenging, if not impossible, to do in experiments. In addition, this was also possible through a comparative analysis between a realistic antenna motion and a hypothetical symmetric motion (not observed in nature).

Appendage- and copepod-scale flow simulations need experiments to provide the accurate anatomical geometries and motions for their results to be biologically relevant. The results of such simulations should be validated against available measurements before applying them to hypothetical cases for comparative analysis. Such complementary simulations and experiments are powerful tools to advance knowledge and test hypotheses on form and function of different biological features or behaviors. There are many possibilities in this line of research as many features of copepods has yet to be simulated. The current simulations were mostly tethered [4] and self-propelled simulations, which have been performed for other aquatic swimmers [55,59,61,62,83], has yet to be carried out for appendage-scale flow of copepods. Such simulations can investigate the role of multiple metachronal beat of legs during multiple hops for feeding or jumps to avoid predators. In addition, the role of hairy appendages versus non-porous ones is an interesting feature to investigate. Nevertheless, simulations with hairy or flexible appendages need: (a) realistic material properties measured experimentally for flexible appendages; and (b) a finite element solver coupled with a fluid solver [49] to solve for the deformations due to hydrodynamic forces.

Apart from investigating the form and function of different biological features or appendages, investigating how the copepods interact with the large-scale surrounding flow is an exciting new avenue for research. It is interesting to see how different flow regimes, e.g., turbulence [15] or vortices [16,84], affect the motion (either cruising or hopping) of copepods. In addition, the effect of vortices or flow gradients on the deformation of hairy and deformable appendages may shed light on the mechanosensory mechaisms of copepods. Appendage- and copepod-scale simulations can help in such studies.

Appendage- and copepod-scale simulations, however, cannot investigate the effect of presence of many copepods on the large-scale surrounding flow. It is well-known that the presence of small microorganisms [28] and even passive particles [77] can affect the large-scale flow and its rheology. Previous simulations of copepods motion in large-scale flow [31,32,34] did not consider the effect of copepods on the large-scale flow, i.e., one way coupled (Section 2.1). To investigate the effect of copepods on the large-scale surrounding flow a multi-scale, two-way coupled simulation is required. One possible way for such multi-scale simulations is to couple appendage-scale (Section 2.2) to large-scale (Section 2.1) simulations. Each appendage-scale simulation will provide the force at the copepod location ($\mathbf{f}$ in the right-hand-side of Equation (1)) for the large-scale simulation while the large-scale simulation will provide the flow conditions (incoming velocity magnitude and direction) at copepod locations for appendage-scale simulations, i.e., two-way coupled. Note such simulations will be computationally expensive and using parallel computing is a must. Finally, such multi-scale simulations do not need modeling (random walks [31] or Lagrangian modeling of jumps [34]) to calculate the trajectories of copepods’ motion in a turbulent flow, thereby providing a more accurate contact rate and clustering of copepods. Contact rates and clustering of copepods are directly related to the probability of encountering food, mates, or predators. Better approximations of contact rate through two-way, multi-scale simulations enable answering important ecological question involving the three major survival tasks, i.e., feeding, predator avoidance, and mate-finding.