A Fluid Dynamics Approach for Assessing the Intelligent Geomorphic Design of the Japanese Pufferfish Nest

Abstract

Research into the geometric nests built by white-spotted pufferfish indicated the nest’s potential for flow control and reduction in flow velocity. However, studies to date have only focused on the construction process and behaviour of the male pufferfish. Hence, the form and functions of the unique features of the nest remain unclear. The present study aims to explore the flow features most useful in understanding the habitat conditions of the nest through a combination of photogrammetric reconstructions of the nest features and two-dimensional (2D) computational fluid dynamic simulations. The findings show the role of the nest structure in reducing the flow velocity and shear stress within the nesting site. Analysis of shear stress indicates that male pufferfish build the outer zones of the nest with coarser material that improves the overall shear strength of these areas. The study identified the function of the nest structure in the protection of the eggs through reduction in flow variations and improved aeration. The addition of shell fragments to the nest peaks by the male pufferfish contributes to the resiliency of the nest structure and ensures a stable bed surface at the central zone.

1. Introduction

Through millions of years of evolution, nature has produced ingenious and practical designs that solve many of today’s technical problems. From drag reduction inspired by sharkskin riblet microstructures [1] to intelligent structures inspired by bee honeycombs [2], researchers are increasingly reliant on bioinspired designs for innovations. Among these designs, one particular structure resembling “crop circles” on the sandy bottoms off the coast of Amami-oshima Island in Japan has baffled native divers for decades. A recent chance encounter led to the discovery of the white-spotted pufferfish (Torquigener albomaculosus) [3], whose adult males, measuring, on average, 109 mm in total length, have been observed constructing the approximately two-meter-wide circular geometric nests [3,4]. These nests have been found at depths of 10 to 27 m and are characterised by an outer zone of radially aligned ridges and valleys, bordering a central zone of shallow maze-like ditches.

The entire structure, shown in Figure 1a, which takes 7 to 9 days to build (Figure 1b-d) [4], is maintained until mating takes place (Figure 1d), and thereafter is left to erode and ultimately flattened by ocean currents (Figure 1e). Based on the observations of Kawase et al. [4], the duration of egg care is approximately 5 days, at which point the nest has undergone significant erosion (Figure 1e). Despite the sheer amount of time and energy it takes to build the structure, and unlike other species that reuse nesting sites [5,6], T. albomaculosus males abandon the nesting site in search of other areas to rebuild the nest. The reason for this behaviour is not yet fully understood, however, Kawase et al. [4] put forward a plausible explanation based on the correlation between the availability of suitable construction material; in this case, fine sand particles and the choice of nesting site.

During their on-site observations, Kawase et al. [4] noticed the deposition of fine sand particles, measuring less than 0.5 mm in grain sizes (Figure 2b), onto the valley surfaces during excavation of the outer zones as a result of the stirring of the surrounding sand (Figure 2a). Subsequently, the current flow and the males’ movement through the valleys resulted in these fine sand particles flowing into the central zone. Hydrodynamic experiments carried out using a half-scale model of the nest showed flow convergence at the central zone through the upstream valleys and divergence of flow through the valleys downstream [4]. These sand particles formed the maze-like pattern in the central zone into which the eggs were released during spawning. However, the water-saturated fine sand particles within the central zone disperse downstream, hence consuming the available particles within a given nest site (Figure 2b,c). This movement of suspended fine sand particles from the ridges and valleys towards the central zone coupled with the observed behaviour of the relocation of the shell/coral fragments by the pufferfish onto the nest peaks results in a redistribution of the particle size inside the nest. The outer zones of the nest are primarily composed of coarser sand particles ranging from 1 to 3 mm in size. Based on these observations, Kawase et al. [4], hypothesised that the nests not only played a vital role in female’s mate choice but also acted as a mechanism to collect fine sand particles at the spawning point. However, the rationale or the characteristics of the nest that ultimately determine the choice of a mate remain unclear.

Recent work by Mizuuchi et al. [8], analysed the behaviour of the pufferfish during the early and middle stages of nest construction. The algorithm produced during this study provided evidence that simple repetitive action on the part of the fish results in the distinct radial patterns observed in the completed nests. Furthermore, they observed that the fish morphology affected the structural features of the nest and these differences in nest topology may act as extended phenotypic signals that influence mate choice.

The research to date has focused mainly on understanding the construction process of the nest and the male behaviour [4,8]. Although some research has been carried out on flow dynamics within the nest [4], no single study exists which investigates this in detail. Central to understanding the flow field is a robust description of the physical model of the nest that enables reliable hydrodynamic experiments and numerical simulations [9,10]. There are no published data on the topology of the nest apart from measurements of nest depth, diameter, size of the central zone, number of ridges and valley pairs and average current velocity near the nest [3,4,8].

Taken together, existing studies indicate that the nest serves as a flow control mechanism that gathered fine sand particles at the nesting site. However, little is known about the interrelations between the various topological features of the nest and its efficacy as a viable nest site. This presents us with some interesting questions. If we consider the upstream ridge, how does the height of the ridge affect the flow within the nest? Are there any recirculation zones within the nest, and if so, how does the ridge height relate to the size of recirculation zones? How do these ultimately affect the transport of sediment particles within and around the nest site, safeguarding the optimal conditions for spawning and hatching its eggs?

In this paper, we use photogrammetric reconstruction and 2D computational fluid dynamics to understand the flow interactions within the pufferfish nests. The geometric nature of the nest is such that the flow field in and around the nest cannot be axisymmetric and the unstable nature of the nest structure requires a complex three-dimensional model. However, the uncertainty in the dimensions of the nest and flow field warrant the use of a reduced dimensional model to capture the critical flow features of the nest that are most useful in understanding of the ecological significance of the nest. The male pufferfish are observed carrying out maintenance of the nest topology until mating takes place. This behaviour indicates the importance of the topology during the early days of egg care; hence, a flow simulation may provide insight into the significance of the structure. 2D simulations are particularly useful in studying the critical flow features of a habitat, as demonstrated by Crowder and Diplas [11,12,13]. Through 2D simulations, we show the formation of a recirculation zone at the centre of the nest that could provide protection from varying flows for the eggs as well as improved aeration rate. Using the distribution of wall shear stress across the nest and surrounding seabed, we demonstrate the significance of the shell/coral fragment placement on top of the nest peaks. We show the low shear stress within the nesting zone and its significance in providing a stable bed surface that is suitable for laying eggs.

2. Materials and Methods

2.1. Nest Reconstruction

Digital photogrammetry techniques such as "Structure from Motion", able to calculate the relative three-dimensional positions of points on an object by monitoring their motion across a number of overlapping images from a single moving camera [14], are increasingly used within academia to process survey images into three-dimensional spatial data such as textured polygonal models. Numerous studies have shown its effectiveness in the reconstruction of topology [15,16,17,18,19], archaeological sites [20], animal morphology and tracks for ecological monitoring [21]. The procedure follows three main steps: (1) matching feature points are mapped across multiple overlapping images to form tie points. During this process, if the data are not provided as input during setup, the intrinsic (such as focal length, camera type, resolution) and extrinsic (camera orientation) camera parameters for each image are estimated by the algorithm; (2) camera parameters are used to calculate the depth information for each image and combined into a single dense point cloud; (3) a 3D polygonal mesh is generated based on the dense point cloud information [22,23].

Previous studies have based their methodology on collecting images and videos from pre-planned surveys alongside direct measurements of the structure being studied. The benefit of this approach is that it provides control over the quality of available data and results in an accurate scaled model that can be verified and validated against direct measurements. The current study utilises existing video footage (see Figure 3 for a summary of workflow), thereby limiting control over the choice and quality of input data. Despite these limitations, Young et al. [24] showed that the Agisoft PhotoScan algorithm (now Agisoft Metashape) is capable of producing accurate and precise models of underwater structures without the need for camera calibration. The paper identifies a similar technique that could be used for historical footage, which provided the inspiration for the model reconstruction in this study. Three-dimensional reconstruction of the nest and surrounding terrain was completed using the standalone photogrammetry software Agisoft Metashape Professional Edition (hereafter Metashape) [25].

2.1.1. Image Acquisition

The video footage used for 3D reconstruction of the nest was acquired from the 2017 documentary series "Big Pacific" by Public Broadcasting Service (PBS) [7]. The duration of the video clip pertaining to the pufferfish nest was approximately nine minutes (29.99 frames per second). Still images were extracted using VLC media player (VLC) [26] at a recording ratio of 30 to obtain the highest possible number of images. Extracted image sequences were manually sorted to closely match the capture requirements given in the user manual. Scenes that reduced the number of "blind zones" and provided good overlap between consecutive images were chosen. A sequence of 23 images was obtained as input for nest reconstruction and a sequence of 21 images was selected for the reconstruction of the surrounding terrain.

2.1.2. Nest Dimensions

Due to the lack of direct measurements of the nest, an image showing the top-down view of the nest containing a pufferfish was chosen for digital measurement using Fiji, a distribution of ImageJ2, an open source image processing package [27,28]. Taking the total length of the fish in the image as measurement of a male specimen given by Matsuura [3] (109 mm), a global scale to measure the nest dimensions can be set (

Table 1. Nest dimensions measured from the top down image.

CZ (mm)	FP (mm)	SP (mm)	D (mm)
709.798	983.674	1476.653	1739.118

). Instead of measuring lengths of individual ditches, the feret diameter of the entire zone was taken since the nest shape was distorted due to the camera angle (Figure 1). To see if this method gave measurements that are representative of the real nest, the ratio of central zone (CZ) to total nest diameter (D) was compared with direct measurements provided by Kawase et al. [4]. Direct nest measurement ratios ranged between 0.38 and 0.50 whilst the digital measurements provided a ratio of 0.40.

2.1.3. Reconstruction of Nest and Surrounding Terrain

Reconstruction was based on the workflow guidelines provided in the software user manual using the recommended parameters and through trial-and-error. The three main steps in the reconstruction are: photo alignment, building dense cloud points, and mesh generation (see Figure 4 for sample outputs from the workflow). The images were imported into Metashape without any further processing. The imported images did not contain sufficient EXIF data for estimating focal length, therefore Metashape classified the images as not calibrated (NC) and assumed the photos were taken using a 50mm lens. Prior to commencing photo alignment, image quality was estimated for each input image. This gave a quality value of over 0.7 for all of the input images, which is above the recommended minimum limit of 0.5.

Photo alignment was carried out at the highest accuracy which increases the computation time but is required for accurate estimation of camera positions. To offset the increased processing time, generic preselection was chosen, which initiates the point matching for overlapping images with lower accuracy settings first. All input images were successfully aligned, giving a sparse point cloud of 162,596 points with a reprojection error of 0.411 pixels (Figure 4a,b). Agisoft defines reprojection error as the average root mean square reprojection error for all input image tie points [22]. Once the sparse point cloud was generated, it was visually inspected to remove outliers. Following the point cloud editing, camera alignment optimisation was performed to improve the model accuracy.

Metashape supports the use of ground control points (markers) to set a coordinate system. Defining a coordinate system improves camera alignment optimisation and sets up an accurate scaled model that can be used for measurements. Pixelation on zooming and lack of unique features for confident marker placement meant that manual marker placement was omitted. The lack of scaling in the reconstruction process is circumvented by proportional scaling of the two-dimensional slices to match the digital measurements in later stages of the study.

Ultra-high-quality dense cloud was built using mild depth filtering which prevents small features such as the intricate patterns in the central zone from being sorted out as outliers. This resulted in a dense point cloud of over 6.5 million points (Figure 4c–f). Metashape calculates the number of depth maps used in creating each cloud point and assigns a point confidence value to each point (Figure 4e,f). This dense point cloud confidence was used to identify areas of low confidence that can be excluded during the final model build. Information from the point cloud was then used to build a high-face-count polygonal mesh. Interpolation was used to reduce holes in the geometry since the dense point cloud does not contain enough information to provide an accurate reconstruction of the surface at certain areas. The computed model contained over 1.2 million faces, which was reduced to 200,000 by performing mesh decimation (Figure 4g,h). The decimated model is exported in OBJ file format to Blender for post-processing.

2.2. Post-Processing and 2D Slice Generation

The polygonal model was imported into the open source 3D computer graphic software Blender [29]. Photogrammetry provided satisfactory reconstruction of three ridge-valley pairs closest to the camera. These three pairs were used in the subsequent 2D slice generation. 2D slices were cut along the ridges and valleys and the resulting side view port renders were exported to Adobe Photoshop and Adobe Illustrator. This procedure was repeated to obtain three slices of the nest from the terrain model. The slices were scaled based on the digital measurements and a reference sketch was created for SolidWorks (Figure 5). SolidWorks was used to produce splines based on the reference sketches, which were subsequently used to create the computational domains, and the resultant 3D models were exported to STAR CCM+ for mesh generation (note that STAR CCM+ requires the import of a 3D CAD model from which users can define a 2D domain). During construction of the splines, the right-hand side ridge was elevated to reduce the height difference between the left-hand side ridge. The large difference in height measured for the right-hand peaks is as a result of the limited footage and camera angle used for reconstruction. It was not possible to validate the digital height measurements of the reconstructed nests due to a lack of direct measurements in the literature.

2.3. Numerical Simulation

2D simulations are carried out over the reconstructed nest geometry, which is a geometric structure with potentially complex 3D flow. 2D simulations are often used for 3D problems due to the immense computational resources required when performing 3D simulations. 2D simulations cannot capture the complexities of the axisymmetric flow field within the nest, however, they can prove beneficial as a first approximation of the flow features that characterize the pufferfish habitat. The commercial computational fluid dynamics software Simcenter STAR-CCM+ is used in this study.

The utility of the 2D simulation is based on the following assumptions:

• A steady incompressible current flow along the maximum dimension of the nest is modelled;
• The angle between two adjacent ridges is sufficiently small so that the flow remains invariant in the radial direction;
• The nest is away from surf zones and coast for the flow to be fully developed.

2D simulation of the flow along the maximum dimension of the nest ensures that the critical flow features of the nest are captured along the direction of the nest erosion. Uncertainty in the direction of the incoming flow as a result of dominant surface wind direction, wind–wave energy transfers and reflections of the waves near the coast are reduced by using a mean velocity profile.

The present study uses V2F K-Epsilon model which is a Reynolds-Averaged Navier Stokes (RANS) turbulence model. RANS simulations are a computationally efficient method of turbulent simulations since they solve for a time-averaged flow. This is done by solving mean flow equations obtained by decomposing the solution variables in the instantaneous Navier–Stokes equations into an average and a fluctuating component. The V2F K-Epsilon model tries to solve the additional stress tensor term that appears in the momentum transport equation as a result of this decomposition by modelling turbulent eddy viscosity using turbulence kinetic energy, turbulent dissipation rate, wall–normal stress component and the elliptic relaxation parameter [30]. The current study requires accurate prediction of flow separation and skin friction, which is possible through the V2F model, since it is better suited for resolving the near-wall turbulence effects compared with other K-Epsilon models [30].

2.4. Computational Domain and Boundary Conditions

The 2D computational domain along with the dimensions and boundary conditions imposed are shown in Figure 6. The size of the domain is 5D × 30D, where D is the nest diameter, with the origin located at the centre of the nest. The inlet is located at a distance of 10D from the nest centre and the outlet is located at a distance of 20D from the nest centre. This is to ensure the boundaries has no effect on the flows upstream and downstream of the nest. Similarly, a depth of 5D is taken to ensure the boundary has no effect on the flow around the nest and seabed.

For the inlet boundary conditions in the simulations, a precursor simulation was used to produce a fully developed flow. The outlet data (Turbulent Intensity, Turbulent Viscosity Ratio and Velocity Magnitude) from the precursor simulation were then used as the input for the main simulations. The inlet boundary for the main simulations were set as a velocity inlet. A mean wave–current velocity profile is used at the inlet for the precursor. Zhang et al. [31] provides an extensive list of investigations into wave–current interactions using analytical approaches, laboratory experiments and numerical methods carried out over the years. For this study, a mean vertical velocity profile for wave–current interactions near Amami-oshima was based on an analytical model proposed by You [32], where a three-layer streamwise velocity profile is given by

$$\overline{u}=\frac{{\overline{u}}_{\ast }}{\kappa }\left(\frac{{\overline{u}}_{\ast }}{u_w^{\ast }}\right)\ln \frac{z}{z_0}{z}_0\le z\le {\delta }_1$$

(1)

$$\overline{u}=\frac{{\overline{u}}_{\ast }}{\kappa }\left(\frac{{\overline{u}}_{\ast }}{u_w^{\ast }}\right)\left\{\frac{z}{{\delta }_1}+\ln \frac{z}{z_0}-1\right\}{\delta }_1\le z\le \mathsf{\delta }$$

(2)

$$\overline{u}=\frac{{\overline{u}}_{\ast }}{\kappa }\ln \frac{z}{z_1}\mathsf{\delta }\le z\le h$$

(3)

Here, $\overline{u}$ is the current velocity, ${\overline{u}}_{\ast }$ is the friction velocity, $\kappa \approx 0.4$ is the von Karman constant, $u_w^{\ast }$ is the wave friction velocity, $z$ is the height above seabed, $z_0$ is the height above seabed where velocity is zero, and ${\delta }_1$ is the wave boundary layer thickness. The benefit of this approach is that most of the input parameters for the model are readily available (wave height $H$, wave period $T$, water depth $h$, bed roughness $K_s$ and reference velocity ${\overline{u}}_r$), the calculations do not require iterations and the authors’ comparison of the model shows sufficient accuracy and agreement with other models. However, the accuracy greatly depends on the measured reference velocity. Kawase et al. [4] noted an approximate current velocity of $0.8 \mathrm{m}{\mathrm{s}}^{-1}$, but the height and method of measurement is not available. Therefore, in this study we assume that the measurement was recorded one meter above the seabed, which is the upper limit for easy and accurate current measurements, as cited by You [32] with reference to Coffey [33]. The input parameters for calculating the wave number and wave friction velocity [34], the mean wave period and wave height around Amami-oshima, were obtained from He et al. [35] for the month of July. The bed properties were assumed to be on the higher end of fine sand, which gives a bed roughness length of 0.3 mm ($K_s=9$) [34]. Sufficient data are not available to determine an accurate distribution of particle size across the nest. However, based on the available literature and photo documents, the particle distribution appears relatively uniform except for the outer and inner regions. A uniform bed roughness that is on the higher end of fine sand is assumed to provide a sufficient representation of the particle distribution. To obtain the density and dynamic viscosity for the physics model, mean salinity at 20 m depth was obtained from Japan Oceanographic Data Center to be 34.46 PSU [36].

The outlet boundary was set as a pressure outlet, the top boundary was set as symmetry plane, and finally, a wall boundary condition was used on the bottom. Wall type essentially applies no slip conditions ($u=v=w=0$) at the boundary. The turbulent boundary layer is divided into three sub-layers defined by the non-dimensional wall distance $y^{+}$, requiring different modelling approaches to resolve the flow. The V2F model allows the use of low-$y^{+}$ and all-$y^{+}$ wall treatment [37]. The low-$y^{+}$ wall treatment resolves the transport equations up to the wall cell, hence a sufficiently fine mesh with cell centroid at $y^{+}~1$ is required [38]. This increases the computational resources, especially for a high Reynolds number flow requiring significantly finer mesh near the wall. On the other hand, all-$y^{+}$ wall treatment tries to cater for both high-$y^{+}$ and low-$y^{+}$ wall treatment subject to the mesh resolution [38]. Computational cost for this method is lower than that of low-y⁺ wall treatment. Ultimately, the decision on wall treatment depends on the available computational resources, the accuracy required in resolving the boundary layer, the nature of the flow problems. The current study uses the low-$y^{+}$ ($y^{+}<1$) to resolve the flow separation that may occur as a result of adverse pressure gradients within the nest.

A polyhedral mesh with prism layers near the bottom wall is used to ensure the required $y^{+}$ value is achieved. Mesh refinement is performed in a 2D radius to achieve a fine mesh resolution around the nest area. To ensure domain independence, the upstream vertical velocity profiles were plotted for a ridge and valley section to check if the flow is fully developed (see Figure 7 for details of the locations monitored and results obtained for one of the ridge sections). The criterion for assessing fully developed flow was based on the idea that the vertical velocity profiles of a fully developed flow will have a constant velocity distribution in the flow direction. For both sections, velocity fluctuations near the seabed were observed. As the flow approached the nest, these fluctuations were reduced and were fully developed at a distance of 1 m from the nest centre.

Furthermore, the outlet was extended to 30D and maximum velocity magnitudes of 84 points were compared with the corresponding points on the 20D domain to evaluate the change. The average change across the 84 points between two coarse mesh was 0.00251 (0.251%), whilst the difference between medium and fine meshes was 0.00243 (0.243%). The difference between the meshes is low, hence the flow is considered independent of the domain.

Grid independence was tested by taking one ridge and valley section. The variations in height and overall shape between the two ridge sections and the two valley sections is small, hence it was assumed that testing one from each section would be sufficient to establish the required mesh resolution. The mean velocity magnitude was compared across 84 points for three different meshes. Due to available computational resources, the limit for mesh base size and cell count was determined through trial-and-error to be 0.1 m with a cell count not exceeding 600,000. Having determined the maximum refinement, the base value was increased by a factor of 4 and halved for each successive mesh. This improved the computational time as the solution from the coarser grid is used as initial solutions for the finer meshes which reduced the convergence time. A minimum convergence criterion for residuals were set at 1.0 × 10⁻⁶.

For the ridge section, the average change between coarse and medium mesh was 0.0003697 (0.037%), whilst the difference between medium and fine mesh was $2.413\times {10}^{-6}$ ($2.413\times {10}^{-4}$%). For the valley section, the average change between coarse and medium meshes was 0.000179 (0.017%), whilst the difference between the medium and fine meshes was $5.967\times {10}^{-5}$ (0.006%). Hence, the fine meshes for both ridge and valley sections were considered to give sufficient mesh resolution since the difference between medium and fine mesh is small.

3. Results

3.1. Recirculation Zones

Plots of Line Integral Convolution and the velocity scalar field were used to visualise the streamwise velocity distribution in and around the nest and identify the recirculation zone on the lee of the upstream ridge. Figure 8a shows the mean streamwise velocity field in the region. Acceleration of the flow is observed as it goes over the first peak of the upstream ridge, forming a dome over the top of the nest with a decrease in the mean flow velocity inside the nest central zone. A point probe placed at the centre of the nest ($x=0$) at a height of 0.03D above the nest measured a maximum mean velocity magnitude of $0.016 \mathrm{m}{\mathrm{s}}^{-1}$s, whilst a probe placed 0.86D upstream of the nest centre measured a value of $0.464 \mathrm{m}{\mathrm{s}}^{-1}$.

The most striking result is the two vortices that were observed within the nest. The smaller one is located between the two peaks of the upstream ridge (see Figure 8a)). The larger one is located just off and above the leeward face of the upstream ridges’ second peak and extends across the central zone, as shown on Figure 8b. Both vortices rotate in a clockwise direction (Figure 8c). Further analysis with a 50% reduction in the inlet velocity showed that the location of the recirculation zones did not change (Figure 8d).

The mean streamwise velocity distribution for the valley sections shows the presence of two small recirculation zone at each of the valleys (see Figure 9). The locations of the vortices inside the valleys correspond to a depression at the valley centres bordering the central area between peaks on adjacent ridges (see Figure 10). A significant decrease in flow velocity was not observed for the valley sections. Furthermore, the vortex structure downstream of the second valley section is observed to be significantly smaller than its counterpart on the first valley section. This difference may be accounted for by the difference in the topology of the depression at the centre of the valleys. The downstream valley depression for the first valley section is larger and steeper compared with the smaller shallow depression on the second valley section.

3.2. Velocity and Pressure Distribution

Vertical profiles of mean streamwise velocity at: x = −0.6D (upstream outside the nest), x = −0.3D (middle of upstream ridge), x = −0.1D (inside the central zone), x = 0.1D (inside central zone), x = 0.3D (middle of downstream ridge), and x = 0.6D (downstream outside the nest) for ridge and valley sections were compared to investigate changes in flow features (see Figure 11).

Figure 12a shows a fully developed flow undisturbed by the nest geometry. Flow reversals in the lower regions corresponding to the recirculation zone can be observed in Figure 12b–d. The flow is reattached by the time it reaches the downstream ridge. A similar behaviour is observed in the valley sections (see Figure 13). Slight variations in velocity profile are observed for the ridge sections. No significant difference is observed for valley section results, apart from the absence of the weaker recirculation zone in the downstream valley depression.

Pressure coefficient distribution for ridge section one across the nest surface is given Figure 14. A sharp drop in pressure is observed at the upstream ridge before gradual increase downstream. This decrease in pressure accelerates the flow over the second peak. Negative pressure is observed throughout the recirculation zones. In the centre of the nest, where the eggs are laid, no sudden change in pressure is observed.

3.3. Shear Stress Distribution

Wall shear stress distribution for ridge section one across the nest and comparison of the nest and surrounding seabed shear stress is given in Figure 15a,b, respectively. Mean seabed shear stress is approximately 0.598 Pa. The highest shear stresses are measured across the outer peaks at the upstream (3.599 Pa) and downstream ridges (2.982 Pa), with the lowest shear stress (0.002 Pa) observed downstream beyond the edge of the recirculation zone (0.2D). Mean shear stress within the central zone of the nest is approximately 0.291 Pa. This is a remarkable observation as the regions of higher shear stress at the outer peaks is where the pufferfish ensures coarser material is present during construction of the nest (see Figure 16), and the central zone is composed of fine sand particles.

The presence of coarser material within a sand mixture has been shown to improve shear strength [31,39]. Furthermore, Yagiz [39] showed that the shape of the gravel within sand mixtures increases the friction and prevents sliding of grain particles. Furthermore, both studies showed that the percentage of coarser material also improves the shear strength of the mixture. Improved shear strength ensures that the structure is protected and slows down degradation. This design is akin to the rubble-mound breakwaters with large stones on the outer regions that act as an armour to protect the inner structures from waves.

4. Discussion

Two recirculation zones were observed in the central zone of the nest as a result of the flow separation near the leeward face of the outer peak on the upstream ridge. The velocity within the central zone is significantly reduced. Two recirculation zones were observed within the central depressions of upstream and downstream valleys. The recirculation zone at the edge of the central zone and velocity reduction may prove beneficial for the eggs at the central zone, which are protected from external adverse velocity fluctuations.

Based on the presence of a recirculation zone, one can surmise that rigorous mixing of incoming flow occurs as a result of the recirculation cell momentum which results in local aeration and prevents the formation of a “dead zone” with reduced mixing and lowered oxygen. Aeration is important in maintaining a sufficient level of oxygen within the zone and prevent abrupt fluctuations in water temperature. The importance of current velocity, aeration, temperature and salinity on the development and survival of embryos has been investigated in various studies [40,41].

During their study on quantifying the spatial flows at various habitat scales, Crowder and Diplas [13] highlighted the use of vorticity and circulation metrics to understand the flow complexities of a habitat and how they can be used in habitat studies. Aquatic organisms tend to use wakes created by topology features such as submerged rocks, boulders, and corals as shelters to minimise energy expenditure, showing the importance of such features for aquatic organisms to flourish. Further on-site studies are required to determine the exact flow features that make the nest a suitable habitat for the spawning of pufferfish.

High wall shear stresses occur at the outermost peaks of the ridges. The pufferfish designs and builds the nest to ensure that the regions of high shear stress are adequately protected through an increase in shear strength. This increase in shear strength is achieved by ensuring that coarser materials are situated on the outer ridges of the nest. It may be the case, therefore, that the pufferfish’s unusual behaviour of decorating the peaks with shell fragments serves a greater purpose beyond that of aesthetics, and in fact ensures that the nest structure is retained for the required duration of spawning and hatching of the eggs. This, in turn, raises the question of whether the size of the outer ridges is relevant in the pufferfish’s quest for constructing a resilient structure through increased duration of erosion for a given flow field. It is possible that the outer ridges may also act as a barrier blocking the transport of foreign particles into the nest. This would prevent suffocation of the eggs due to deposition of fine particles on top of the eggs and damage to eggs due to larger particles. Furthermore, the average shear stress within the central zone was observed to be lower than the surrounding seabed. It is possible, therefore, that the lower shear stress ensures a stable bed surface for laying the eggs.

The findings of this research may well have implications for engineering design and provide a bridge between various fields and engineering such as ecology. Degradation of coral reefs and seagrass beds across the globe has resulted in restoration efforts through measures such as colony transplantations. Environmental factors such as local hydrodynamics play a significant role in the success of these efforts [42,43]. Similar geometric structures to the pufferfish nests could be employed in transplantation sites that are hydrodynamically exposed, since such a structure provides means of flow control, bed stability and protection from external forces. The recirculation zone presents the added advantages of improved aeration and nutrient circulation for the growth of polyps and seagrass.

5. Conclusions

This study set out to explore the critical flow features within the geometric nests built by the pufferfish, required for understanding the nest’s habitat conditions and the significance of the topological features to engineering. This study has shown how digital photogrammetry can be used alongside historical video footage to reconstruct the nest and surrounding terrain features. CFD simulations were conducted on the 2D cross-section of the nest along the flow direction with the maximum dimension, to assess the effect of the nest’s topographic features (such as valleys and ridges), on the flow field. The investigation of pufferfish nest has shown that:

• The topology of the nest contributes to the reduction of flow velocity and formation of a recirculation zone within the nest site that ensures eggs are protected from adverse velocity fluctuations, improved aeration and additionally act as a barrier blocking transport of foreign particles into the nest;
• The previously observed behavior of male pufferfish decorating peaks with shell fragments ensures a resilient structure through improved shear strength in regions of high shear stress;
• Lower shear stress is seen within the central zone of the nest, ensuring a stable bed surface for spawning of eggs.

The scope of the study was limited in terms of available computation resources. Hence, simulation of the nest geometry was limited to 2D simulations. However, the 2D simulations could prove as a good first approximation of the flow field within the nest for studying the habitat conditions of the nests. A natural progression of the current work, therefore, would be to perform a 3D simulation to validate the results obtained from the current 2D simulations and model the complex 3D flow of the nest. Hydrodynamic experiments would prove beneficial in determining the validity of the numerical simulations. However, the numerical approach and assumptions carried out in this study were aligned with the existing best practices in the literature. Coupled with wave makers, experiments could provide insight into how wave–current interactions affect the nest, and how the nest form influences its function.