Measuring the 3D wake of swimming snakes (Natrix tessellata) using volumetric particle image velocimetry

Undulatory swimming kinematics are usually classified into four main modes involving different proportions of the body and/or caudal fin (Lindsey, 1978; Sfakiotakis et al., 1999). Among these, anguilliform swimming describes the motion of elongated animals, where the kinematics consist of an undulation that increases in amplitude along the body. Since Lighthill's analytical large-amplitude elongated body theory (Lighthill, 1971), numerical investigations have been conducted (Kern and Koumoutsakos, 2006; Borazjani and Sotiropoulos, 2009, 2010; Nangia et al., 2017; Battista, 2020a,b; Khalid et al., 2020, 2021; Ogunka et al., 2020) in order to obtain an estimation of quantities such as forces, pressure fields and swimming efficiency for an idealized swimmer. Anguilliform swimming hydrodynamics have been experimentally studied in the eels Anguilla anguilla (Muller et al., 2001) and Anguilla rostrata (Tytell, 2004; Tytell and Lauder, 2004), the lamprey Petromyzon marinus (Gemmell et al., 2016; Du Clos et al., 2019; Lehn, 2019) and the catfish Plotosus lineatus (Tack et al., 2021), using particle image velocimetry (PIV) in 2D.

The elongated and limbless morphology of snakes is associated with over 10 different gaits used in various environments (Jayne, 2020), with swimming kinematics first studied quantitatively by Taylor (1952), inspired by the work of Gray (1933). Snakes' great diversity in size and shape makes them ideal to study anguilliform locomotion, and their numerous independent radiations into aquatic environments provide insight into the adaptation of elongated limbless tetrapods to this environment. However, apart from prey capture studies (Van Wassenbergh et al., 2010; Segall et al., 2019), a hydrodynamic analysis based on live swimming snake data has yet to be attempted.

Although 2D PIV helps us to understand aquatic locomotion dynamics, 3D velocimetry is necessary to fully describe the vortex dynamics produced by a moving finite body (Tytell et al., 2008), especially in out-of-plane movements as in swimming snakes. The main methods are synthetic aperture PIV (e.g. Lehn, 2019; Mendelson and Techet, 2020), tomographic PIV (Skipper et al., 2019) and digital defocusing PTV (e.g. Flammang et al., 2011; Bartol et al., 2016). These methods instantly capture the entire 3D unsteady flow field, but flow field measurements of swimming animals are scarce. Aside from unpublished work on lampreys (Lehn, 2019), the 3D wake of an anguilliform swimmer has not been characterized experimentally.

We propose a method to examine the 3D flow induced by an anguilliform swimmer using a free-swimming dice snake (Natrix tessellata) as our model. Using DDPTV (defocusing digital particle tracking velocimetry), we measured the wake and investigated the vortical structure. We compared results with previous computational fluid dynamics and 2D experimental PIV studies.

Swimming kinematics and PIV data were obtained from three (one male and two females) adult captive dice snakes, Natrix tessellata (Laurenti 1768). Dice snakes are semi-aquatic snakes and were chosen for their relatively small size.

Although conducted separately from the PIV experiments, these trials have the potential to be performed concurrently in a separate study, with appropriate laser filtering in place. Movie 1 shows an example of a swimming trial video.

In order to quantify the flow around a swimming snake, we used a volumetric three-component DDPTV setup (V3V-9000-CS system, TSI; Fig. 1B). The experiments were conducted in a water-filled Plexiglas tank (length 210 cm, width 18 cm, height 30 cm, water height 22 cm). A narrower tank was chosen in order to increase the chances of the snake swimming through the measurement volume without over-constraining the kinematics. The water was seeded with 50 μm polyamide particles (PSP-50, Dantec Dynamics) with a concentration of around 5×10⁻² particles per pixel (ppp) (Cambonie and Aider, 2014). The tank was illuminated from above by a 14 Hz pulsed 200 mJ dual head Nd:YAG laser (Quantel Evergreen) expanded with two −25 mm cylindrical lenses. The particles were filmed by three CCD arrays (Powerview Plus 4MP-HS) with 50 mm camera lenses (Nikon AF Nikkor, aperture f/16) mounted on a 170 mm equilateral triangular frame (V3V 9000-CS). The acquisition setup was piloted by the INSIGHT V3V-4G software and calibrated as explained by Troolin and Longmire (2010). The resulting measurement volume is the intersection between the fields of view of the cameras and the laser cone. Its dimensions are approximately 18×18×12 cm and it is centred on the water tank in width (z-axis) and height (y-axis).

During each trial, the snake was placed at one side of the tank and allowed to swim. When put into the water, snakes typically swam to the other side of the tank along the x direction, passing through the measurement volume while being recorded by the cameras at a frequency of 14 Hz. Preliminary trials were conducted to estimate the best time difference (ΔT) between two consecutive PIV frames so that displacement of the seeding particles did not exceed 6 pixels. The resulting ΔT ranged between 1 and 2.5 ms. As the width of the measurement volume was smaller than the width of the water tank, only those trials where the snake swam completely through it were retained. The PIV swimming trials were performed on six different days over the course of several weeks. During the first two recording days, the females were gravid, which impacted their body mass and diameter. The measurements in Table 1 were made after the eggs had been laid. A total of 112 sequences of free-swimming N. tessellata with at least one tailbeat in the measurement volume were recorded. The number of sequences was respectively 31, 52 and 29 for Natrix M, Natrix F1 and Natrix F2.

In order to avoid artefacts induced by the presence of the snake's body on the raw images during the PIV analysis, the snake was removed from the images with the help of a custom-written masking routine in MATLAB using the Image Processing Toolbox.

The images were then processed using INSIGHT V3V-4G software in order to obtain the vector fields. The main steps of the data processing are as follows. (1) 2D particle identification in the images of each camera using a Gaussian fitting (Levenberg–Marquardt algorithm). (2) 3D particle identification, i.e. the correspondence of a particle in one image to the same particle in the other two images is determined using the volume calibration. The fine search tolerance is set to 0.5 pixels and the coarse search tolerance to 1 pixel. (3) Velocity processing using the 3D displacement between two consecutive frames with a relaxation algorithm (Pereira et al., 2006). (4) Neighbour tracking reconstruction using the trajectories from the neighbouring particles to find the probable location of a missing triplet at time t+ΔT. (5) Interpolation of the randomly spaced vectors using a Gaussian-weighted algorithm (Pereira and Gharib, 2002), resulting in a 74×74×43 evenly spaced vector grid with a 2.5 mm resolution.

Vector fields were processed using code written in MATLAB. As the swimming trials are in free-swimming conditions, the movements of the snake in the measurement volume can be quite different from one sequence to another, which makes the comparison challenging. In order to tackle this problem, we changed the regular Cartesian coordinate system (x,y,z) to a new rotated coordinate system (x_rot,y,z_rot) for each swimming trial. For each recorded sequence, the water jet created by the last tailbeat was tracked with a velocity thresholding method. The trajectory of the jet was estimated by a linear regression of the superimposed position of the jet centre. A new base was then created by a rotation around the y-axis so that the new x_rot axis was co-linear to the jet's trajectory and the z_rot axis was perpendicular to it.

Although the snakes swam with the characteristic lateral undulation of an anguilliform swimmer in the frontal plane (x–z), there were also noticeable, albeit smaller, oscillations in the sagittal plane (x–y). The tail was often positioned below the body and the tailbeats were downward facing. The mean swimming speed in body lengths per second across the swimming trials was U_M=0.51 L s⁻¹, U_F1=0.37 L s⁻¹ and U_F2=0.34 L s⁻¹ (Fig. 1C).

During a typical swimming sequence, vortices were shed after each change of direction of the swimmer's body during oscillatory movements. The displacement of the body created two vortex tubes of opposite rotation (Fig. 2A,B) above and below it, which in some cases became connected. The shed vortices were tilted at an angle of approximately 30–40 deg to the swimming direction (Fig. 2C). Their trajectory was mostly linear with a tendency to curl. Body vortices with similar topologies and dynamics were also observed.

The circulation of the vortex cores rapidly increased in magnitude just after the tail beat and reached a maximum before decreasing (Fig. 2D). The positive and the negative vortex cores displayed symmetrical variations during the observed time period. During its lifetime, the vortex size gradually increased (Fig. 2E) and changed from a flattened into a more circular shape. The kinetic energy of the wake increased starting from the moment when the swimmer entered the measurement volume. After the final tailbeat of the swimmer, the perturbations dissipated and the kinetic energy decreased (Fig. 2F).

Across the swimming sequences, the dimensional maximal circulation ranged from around 18 to 51, 25 to 70 and 22 to 88 cm² s⁻¹ for Natrix M, Natrix F1 and Natrix F2, respectively. The non-dimensional circulation Γ* ranged between approximately 0.5 and 1.5. The maximum vortex size was between approximately 2.5 and 5.3 cm, 2.7 and 5.2 cm, and 2.6 and 5.4 cm for Natrix M, Natrix F1 and Natrix F2, respectively. A Spearman correlation (R=0.8, P<0.05) showed that there was an increasing monotonic trend between the maximal vortex size and the maximal circulation. As the maximal circulation is usually measured near the beginning of the vortex lifetime, this means that a high initial circulation is likely to yield a bigger vortex.

A Spearman correlation (R=0.7, P<0.05) between the maximal kinetic energy in the whole measurement volume during the entire swimming sequence and the maximal circulation of the tailbeat vortex showed an increasing monotonic trend. This means that there is a tight scaling between the local hydrodynamic parameters of the vortices (i.e. the circulation) and the global parameters (i.e. the kinetic energy of the fluid in the whole measurement field).

As the wake is the footprint of a swimmer, it is crucial to fully capture it in order to maximize the understanding of the fluid–structure interactions between the body and the water. The 2D PIV experiments on eels by Muller et al. (2001) in the frontal plane have led to the prediction that the 3D vortical structure of an anguilliform swimmer should resemble a vortex ring. Tytell and Lauder (2004) were sceptical about the true roundness of the vortex and the computational work of Borazjani and Sotiropoulos (2009) suggested that the wakes were in fact vortex loops stretched in the stream-wise direction or ‘hairpin-like’ vortices. The traditional plane to measure the swimmer's wake in the existing literature is the frontal plane (x–z) without looking at the other ones (Fig. 3A–C). Our experimental results in the frontal plane are in accordance with the 2D vortex behaviour in the existing anguilliform swimming studies (Fig. 3D) with the creation of two counter-rotating vortex cores during a tailbeat cycle. The 3D velocimetry method presented here reveals that these vortex cores observable in the frontal plane are the links between the vortex tubes stretched in the swimming direction (Fig. 3E). These links are created after each change of direction of the tail.

The observed 3D vortex topologies were hairpin like, tubular or turbulent without any clear cohesion. We observed hairpin-like and tubular structures in most of the analysed swimming sequences. A lack of cohesion was observed essentially only when the snake kinematics involved strongly out-of-plane motions, changes of direction, or acceleration/deceleration sequences. These different vortex topologies have been predicted by computational fluid dynamics (CFD) studies (Kern and Koumoutsakos, 2006; Borazjani and Sotiropoulos, 2009; Khalid et al., 2021), and associated with different kinematic parameters such as the tailbeat amplitude, frequency and body wavelength.

As DDPTV enables the observation of the full wake structure, our measurements revealed the different stages of the formation of the hairpin vortices. The shedding of two vortex tubes of opposite rotation immediately after a tailbeat (Fig. 2G) creates a bridge where the vortex tubes interact (Fig. 2H), leading to their connection and detachment before the bridge (Fig. 2I). These steps are similar to events in coherent packets of hairpin vortices in channel flow (Zhou et al., 1999). The rounding of the vortex structure after the bridging step is also reminiscent of the Ω-shaped vortex created by self-induction of the hairpin vortex. Despite the relatively slow frame rate, our 3D measurement allowed us to observe the bridging step.

The shedding of vortices was also witnessed more anteriorly on the body. The production of body vortices has been described by Gemmell et al. (2016) for lampreys, where particular attention was given to the bending of the body. In addition to transverse motion, the body has an angular velocity in the transverse plane (y–z). The body surface rotation helps to rotate the fluid and eases the shedding of adjacent vortices. The production of counter-rotating vortex tubes along the body has been suggested to generate additional thrust (Fu and Liu, 2015) during swimming. In the case of gliding snakes, a CFD study showed that the body-induced vortices produced lift along the body (Gong et al., 2022). The body vortices may therefore have a double function of propelling and stabilizing the swimmer.

The maximal circulation measured for the three swimmers was similar to those measured in smaller anguilliform swimmers: ∼20 cm eels (Tytell and Lauder, 2004) and ∼10 cm lampreys (Gemmell et al., 2016). The circulation of small hairpin vortices decreases exponentially after attaining a maximum because of the viscous interactions of the counter-rotating vortex tubes (Sabatino and Maharjan, 2015). In contrast, the circulation of laminar vortex rings slowly decreases after the maximum (Rosenfeld et al., 1998). In our results, depending on the swimming trial, the circulation either decreased or stayed relatively constant. The decreasing vortex circulation is probably due to the viscous interactions with either the counter-rotating vortex tubes or the surrounding vortices created by the other tail beats. The relatively constant or slowly decreasing circulation, in contrast, behaves like a laminar vortex ring.

Studies on vortex rings show that the diameter of the vortex rapidly increases, attains a maximum, and then asymptotically decreases to a constant value (Arakeri et al., 2004). The distance between the vortex cores observed here follows the same pattern as they either increase or stay relatively constant. This means that the vortices with an increasing distance are in the first part of their lifetime and the ones with a relatively constant distance are likely to be near the end of their lifetime. Our data further suggest a relationship between the circulation and vortex size.

As the experiments were done in a still-water tank, the displacement of the water was essentially caused by the swimming snakes. The kinetic energy of the water measured is therefore the energy that was transferred by the snake minus the dissipation. Computational studies found faster swimming kinematics result in stronger vortices (Kern and Koumoutsakos, 2006; Bhalla et al., 2013). Gravid snakes had larger body diameter and higher mass but proportionally less locomotor muscle due to eggs, which may impact locomotor performance (Seigel et al., 1987). Gravid snakes shed vortices with high circulation and size, producing the highest overall kinetic energy, suggesting they adjust swimming effort at a higher cost. Complementary results from other species may provide an explanation for the observed phenomenon of gravid snakes exhibiting less swimming activity compared with non-gravid individuals (Aubret et al., 2005).

Volumetric measurements of the wake can be used as tools for non-invasive estimations of swimming performance (Li and Mendelson, 2023). Coupled with 3D kinematic data, the results obtained using this method could be used for interspecific comparisons of the swimming performance of snakes. This would lead to a better understanding of the differences of swimming efficiency in snakes with different lifestyles, shapes and ecologies.

We thank Jake Socha and an anonymous reviewer for helpful and constructive comments on an earlier version of the manuscript. We also thank the MECADEV laboratory's Plateau d'analyse du mouvement for experimental setup assistance.