Sponge Density and Distribution Constrained by Fluid Forcing in the Deep Sea

Understanding the distribution and density of sponge grounds in the deep sea is key to appreciating the ecological importance, vulnerability, and role in ocean biogeochemistry of these important habitats. A novel combination of clustering analysis and fluid flow finite element modelling was used to study distribution of four sponge grounds from the Labrador Sea, which were chosen for their distinct assemblages of sponges. Significant small-scale clustering patterns of sponges were found within each sponge ground, measured using the Ripley K function. A new approach using finite element modelling of fluid flow was then applied to understand the drivers of the observed sponge distribution, with detailed numerical experiments providing insights into the flow patterns that could not be obtained from field measurements. Simulations using idealised sponge shapes suggested that sponge wakes could interact and influence the mean flow conditions at the spacings observed within the sponge grounds. Simulations of flow over topographic models of the sponge grounds showed that these interacting wakes generated a boundary layer of slowed flow, which may be beneficial to sponge development. The boundary layer may be acting as a protective layer, affecting larval recruitment, increasing particle fall out and increasing the effective radius of pumping. This observation has important implications for the impact of anthropogenic damage on sponge grounds, which changes the sponge distribution and may reduce the boundary layer thickness, affecting potential recovery rates. This study illustrates the power of incorporating mathematical modelling with field observations in the deep sea to increase our understanding of anthropogenic and climate change impacts on sponge ground ecosystem structure and function.


INTRODUCTION
Sponges (Phylum Porifera) are an often overlooked as important components of benthic marine ecosystems. They are known to perform several functional roles notably as reef building organisms in both the shallow and deep-sea (Miller et al., 2012;Knudby et al., 2013;Roberts et al., 2018). Sponge grounds support areas of high biodiversity (Klitgaard and Tendal, 2004;Murillo et al., 2012;Hawkes et al., 2019), act as nurseries for many fish species (Klitgaard, 1995;Freese and Wing, 2003), and can influence primary production by controlling nutrient availability (Bell, 2008;Maldonado et al., 2012;Murillo et al., 2012;Kutti et al., 2013;Maldonado, 2016). This ecological importance is complimented by their economic potential from bioengineering (Wu et al., 2019) and novel chemicals for pharmaceutical applications (Koopmans et al., 2009;Anjum et al., 2016). However, despite their ecological and socioeconomic importance, sponges are vulnerable to anthropogenic stressors such as trawl fishing (Kenchington et al., 2014;Pusceddu et al., 2014;Victorero et al., 2018) and oil and gas activities (Vad et al., 2018). Sponge ground locations and densities need to be known, in order to quantify their distribution and function, and-ultimatelyprotect the critical habitats that they form. In the deep-sea, observational data are very limited due to the difficulties in conducting science at depth, such as the cost and limited field-of-view of remotely operated vehicles (Henry and Roberts, 2014a,b;Howell et al., 2014). Computational models can be used to interpolate and quantitatively predict sponge distribution, helping to circumvent the challenges involved in understanding the factors that control sponge distribution. The large (basin) scale distribution of sponge grounds in the North Atlantic has been studied. However, sponge distribution and density within deep-sea sponge grounds has been understudied and remains poorly constrained (Klitgaard and Tendal, 2004;Cathalot et al., 2015;Kazanidis et al., 2019).
Studies on sponge density and distribution in shallow water environments have been carried out using standard ecology survey techniques, with multiple studies across environments and latitudes (Wilkinson and Evans, 1989;Roberts and Davis, 1996;Xavier and van Soest, 2012). However, the major drivers of these distribution patterns such as wave damage, water opacity, and light levels do not operate in the deep-sea. Empirical analysis of localized distribution in deep-sea sponge grounds is very limited with most data being purely observational or on species occurrence only (Henry and Roberts, 2014a). Studies of sponge density have been carried out in the Faroe-Shetland Channel where fisheries pressure and substrate type (sand, gravel, cobble) were statistically significant factors (Kazanidis et al., 2019). However, a host of other chemical oceanographic variables (silicic acid, temperature, etc.) also exhibited strong controls over sponge density (Kazanidis et al., 2019). The influence of multiple forcing mechanisms results in significant challenges in scaling-up from a small number of observational data points to whole-basin sponge ground densities using oceanographic variables (Howell et al., 2016). Creating mathematical simulations of current conditions in the deep sea allows us to explore whether the observed distribution of sponges is related to prevailing currents. This in turn provides possible explanation of the patterns seen which could be scaled up to infer large-scale distributions. The simulations of the observed sponge grounds also allow a more detailed investigation of the small-scale flow patterns that would not be possible to measure the field, with the potential to increase our understanding of this important and understudied habitat.
Patterns of deep-sea sponge distribution and density have been attributed to several factors in large scale biogeographic studies. The primary drivers of sponge distribution are thought to be silicic acid concentration (Beazley et al., 2015;Howell et al., 2016) temperature, depth (Rice et al., 1990;Bett, 2001;Murillo et al., 2012), and particulate carbon availability (Barthel et al., 1996;Howell et al., 2016), in addition to topography, internal waves, and current speeds which control flow conditions at the seafloor (Klitgaard and Tendal, 2004;Howell et al., 2016). It would be logical to assume these factors would also play a role at the smaller scale and that sponge distribution would potentially be optimized to improve their localized conditions (Kutti et al., 2013). The arrangement of organisms can affect the flow around individuals, either spacing themselves the maximum distance from each other to receive the greatest unprocessed flow potential, or clustering to slow the flow for potentially greater particle dropout or protection from damage. These patterns in distribution are seen abundantly in nature such as in corals (Chindapol et al., 2013;Cresswell et al., 2017) and terrestrially in trees (Holtmeier and Broll, 2017).
Here, we will investigate whether the patterns of density and distribution at small spatial scales in four different sponge grounds in the Labrador Sea are linked to the prevailing flow conditions, and will explore possible explanations for the patterns observed. We discuss the implications and possible explanations of our results, as well as their potential impacts on the greater understanding of the role and sensitivity of sponge grounds in the deep-sea.

Sponge Ground Imaging
Sponge ground images were recorded on the DY081 RRS Discovery cruise summer 2017 at three locations: Orphan Knoll, Nuuk, and Southwest Greenland (SWGL) (Figure 1) each chosen due to their differing water mass interactions. High-definition video was recorded via Remotely Operated Vehicle (UK National Marine Facilities ROV Isis) using a fixed view Scorpio camera. Images were recorded from a typical height of 1 m above the sea floor, with the field of view varying in width between 1 and 3 m. Lasers were attached to the Scorpio camera at 10 cm spacing so scale could be accurately resolved from images. The study areas were randomly selected within each location. A single study area was chosen for Nuuk (Figures 1, 2A) and SWGL Localities (Figures 1, 2B), which were uniform in depth and have no dominant species. Two study areas were assigned at Orphan Knoll, due to the presence of two distinct sponge grounds at separate depths characterized by different dominant groups. These are (i) a Geodia dominated sponge ground (GOK; Figures 1, 2C) and (ii) a hexactinellid dominated sponge ground (HOK; Figures 1, 2D) at greater depth. Still images were the taken from the video at 2 s intervals which was sufficient to ensure that typically 30 individual, overlapping images were available to create a digital surface model of between 3 and 5 m in length.

Digital Surface Model
Surface mesh models were created from sets of images for each study area using photogrammetry. Agisoft PhotoScan  software was used to align the images and create a 3D surface mesh using structure-form-motion algorithms (Westoby et al., 2012); the SGM algorithm (Hirschmüller, 2005) and scale invariant feature transform (SIFT) operator (Lowe, 2004).
A scaled 3D surface mesh was exported as a.stl file and a 2D point cloud map of sponge locations exported as a point pattern dataset. The 3D surface mesh was converted using AutoCAD Software to AutoCAD Drawing format to FIGURE 3 | Plots of the Ripley K function clustering analysis (K) for the four study areas (A) HOK study area, (B) Nuuk study area, (C) GOK study area, and (D) SWGL study area. Observed value for K [K obs (r)] theoretical random distribution value for K [K theo (r)] < 95% confidence limit of random distribution between upper confidence bound [K Hi (r)] and lower confidence bound [K lo (r)].

Clustering Analysis
A point pattern map for each locality was input into RStudio, such that each point represents the position of an individual sponge. Multi-distance spatial cluster analysis was carried out using Ripleys K functions (Ripley, 1976(Ripley, , 1977(Ripley, , 1981. The K-function is defined by: Where d ij is the Euclidean distance between ith and jth points in a dataset of n points, t is the search radius, λ is the average density of points (generally estimated as n/A, where I is the area which contains all the points), and πt 2 is the indicator function (1 if its operand is true, otherwise 0). K d approximates πt 2 if the points are homogenous with the region. This radial spatial analysis can become unrepresentative when applied close to the edges of the domain, because of the truncation of information at the edge. We used periodic boundary conditions (such that an edge of the domain is joined to its opposite edge and analysis treats this as a continuous surface) to eliminate the edge effect, which may have been significant given the narrow transects.
Upper and lower bounds of confidence of statistical significance were produced using Monte Carlo simulations. Monte Carlo simulations randomly generate a distribution of points equal to the number of input points. Multiple simulations were run and used to ensure clustering patterns were not due to random variation (at the 95% confidence level).

Fluid Flow Finite Element Analysis
We modeled the flow around idealized sponge shapes and around the 3D surface models of the sponge ground study areas using the finite-element package Comsol Multiphysics 5.1. The dimensions of the idealized sponge shapes were based on average sponge sizes of geodia and Euplectella species collected on the DY081 Cruise. The flow field was described using the Reynoldsaveraged Navier-Stokes (RANS) equations coupled with the K-epsilon (K-ε) Turbulence model (Hanjaliae and Launder, 1972). The flow was computed using the 3D surface meshes generated from photogrammetry, with a no-slip basal boundary condition, open boundaries (zero viscous stress condition) for those perpendicular to flow, downstream and upper boundaries, and a uniform inflow velocity at the upstream end. Initial flow conditions within the model area for the simulation calculations were set as a steady flow equal in velocity to that of the inflow. Surface models were orientated so the simulated flow matched the observed current direction over the sea floor recorded at the time of the ROV survey. A mesh was created using free tetrahedral elements for the geometries using the finite element mesher by Comsol. A coarse mesh was used to enable an efficient exploration of model parameter space, with some simulations being run at finer mesh resolution to ensure grid scale independence for the simulation results. The inbuilt material properties for water were used with fixed solid geometry for the idealized simulations. In all cases, three-dimensional, incompressible (constant density) flow of water was simulated with idealized structures and sea floor reconstructions were held stationary.
Inlet and initial flow speed were varied between 0.1 and 1 m/s in 0.1 m/s intervals (Reynolds numbers for all simulation fall between approximately 1.53 × 10 5 and 7.5 × 10 6 ). The volumetric flow rate used was matched to the mean bottom flow rates observed using the sonde on the DY081 RRS Discovery cruise (Hendry, 2017) and the range based on yearround observations from Orphan Knoll (Greenan et al., 2010). A stationary solver was used to compute the steady-state flow patterns and a Turbulent flow k-ε model was used to solve the RANS equations and conservation of mass.
An initial simulation of three-dimensional flow around a sphere was used to benchmark the main Comsol simulations.
The benchmarking method laid out in Nita and Allaire (2017) was used with drag coefficients for the flow around the sphere being in the range of 0.46-0.5 to verify the robustness of the numerical method used. Results were visualized as twodimensional cross sections of flow velocity magnitude and wake length, and diameters were plotted against height and diameter for the idealized simulations.

Clustering Analysis
All study areas exhibited non-random (p < 0.05) sponge distribution patterns, with variation in the presentation of distribution between sites (Figure 3). The SWGL ( Figure 3D) and GOK (Figure 3C) sites exhibited similar patterns of distribution. Individual sponges cluster together at distances of 0.1-1.2 m and the distribution becomes more uniform at distances greater than 1.4 m indicating that they form dense, sparse clusters. The HOK (Figure 3A) site shows a very strong clustering pattern within the study area, although there is a random distribution at distances less than 0.25 m. The Nuuk ( Figure 3B) study site in contrast exhibits a consistent nonrandom but weak dispersion patterning.

Models of Flow Around Idealized Sponge Morphologies
We initially used two simplified shapes to represent sponge morphology: a sphere and an ellipsoid. For both the spherical structure (Figure 4) and the ellipsoidal structure (Figure 5), the same patterns of wake (downstream regions of reduced flow velocity) are produced. There is a direct correlation between structure height and length of the initial wake, as well as the diameter of the wake and diameter of the structure (Figure 6). Wake length and diameter are mostly unaffected by flow speed with the changes of flow speeds within the wake being directly proportional to the changes to the ambient flow speed (Figure 6). The length of the effective slowing due to the wake is approximately 10 times the height of the structure, which is similar to wake lengths of submerged hemispherical objects (Shamloo et al., 2001). Beyond this wake region, the near-wall flow velocities are similar to those upstream. A second feature of both models is the presence of an area of lower velocity upstream and around the sides of the sponge structure caused by flow separation, which increases in size at high velocities (Figures 4, 5). These initial tests indicate that, for the sponge spacings (up to 1.4/m 2 ) and heights (up to 0.5 m) in the four study areas, we expect the wakes of the sponges to be interacting and influencing the flow conditions.

Models of Flow at the Study Area
For each of the study areas, we ran simulations with sponges in situ and removed. Individual sponges created wakes that were significantly deeper than the boundary layer resulting from flow over the sea floor topography. The structures in the Nuuk study area (Figure 7) produced a boundary layer of similar height to that of the seafloor topography boundary layer but with a significantly slower velocity, which appears reduced by approximately 50% from the same point within the boundary layer. The area of wake influence from each sponge structure extends horizontally in a narrow cone pattern but the distance between each structure is such that a coherent and continuous boundary layer is produced, with properties set by the sponge height and spacing. Both SWGL (Figure 8) and HOK (Figure 9) sites show similar flow patterns as the Nuuk site but with slightly thinner boundary layers. Both SWGL and HOK have clusters of sponge structures sitting within small depressions causing the sheltering effect of the topography to make a greater contribution to the overall flow pattern. The GOK site (Figure 10) differs due to the prominence of a single large sponge > 0.5 m in diameter and height. The flow patterns around the large sponge result in a significantly sheltered area where smaller sponges are clustering.
There are other small clusters of medium-sized sponges causing

Do Sponges Cluster?
Patterns of distributions within sponge grounds have been noted before, but there has not been an analysis to date of the clustering patterns themselves (Kutti et al., 2013;McIntyre et al., 2016). Clustering analysis gives a more robust approach to studying these patterns, by removing observational bias and allowing robust statistical tests to be performed to validate observations and hypotheses (Baddeley et al., 2014). As the sponge grounds chosen as study sites are defined by distinct assemblages of sponges, it is interesting but not surprising that different functional groups follow distinct clustering patterns. The ability to scale up the distributions described here are limited due to the small study area. Their size is limited by the diameter of the areas surveyed by the ROV (Hendry, 2017;Hendry et al., 2019). However, our flow-modeling method could be extended to a broader geographical area, and combined with localized environmental data and predictions, to make inferences about environmental driving factors for sponge distribution in sponge grounds on wider spatial and temporal scales. Uniquely, a modeling approach also allows us to simulate a great range of environmental conditions and therefore can be used as a predictive tool for both potential hypotheses and expected optimal conditions. Why Are Sponges Clustering?
The variability and the driving factors behind these distributions could be wide-ranging, from dispersion of gametes (Abelson and Denny, 1997) to food availability (Robertson et al., 2017). One factor that can be studied within limited geographical areas is localized flow around the sponges. Is this driving distribution patterns in deep-sea sponge grounds? This new application of modeling techniques to localized flow around sponge allows us to address such questions, previously not possible from observational approaches, especially in deepsea environments.
We need to quantify the pattern of flow around a single structure to understand how fluid flow could be affecting the distribution of individuals within a sponge ground. The creation of both idealized single organism models and study area models give us a better understanding of this flow. The idealized models give us the predicted size and length of the wake given average heights of the sponges. In our study areas, the length of the wake is greater than the spacing between individuals in the sponge grounds, so the wake effect of other sponges is predicted to influence organism structures. The wake length-scales found in this study suggest that sponges that appear isolated may be influenced by the wake of other individuals to a greater extent than previously expected. When the flow around the sponge structures within each study area is modeled and compared with the same seafloor surface without these structures, a significantly thickened boundary layer is observed as predicted by the idealized models (Figures 4, 5). This slowing effect is due to the overlapping of the wakes of each individual sponge to create a relatively uniform area of slowed flow. Directly linking the slowed boundary layer as the driver behind patterns of distribution is difficult and would require many more and varied study sites. However, the observation of this slowed boundary layer and its potential effect as a factor in sponge distribution is important in its own right. Are the sponges distributing themselves to benefit from the slowed flow in this layer and, if so, what advantage does it provide?

Advantages of Manipulating the Boundary Layer
The effects of flow on an individual sponge organism broadly fall into three categories: gamete dispersion (Abelson and Denny, 1997), nutrient supply (Kahn et al., 2018), and mechanical forcing (Palumbi, 1986). The presence of a slowed boundary layer whose thickness is governed by the velocity of the flow and roughness of the surface has the effect of reducing any mechanical forcing. This lowering in the forces applied to a sponge will reduce the likelihood of damage or attachment failure due to high flow. Average flow velocities in deep sea sponge grounds are low in our study sites (∼0.1 m/s) (Hendry, 2017;Hendry et al., 2019) so this protective effect is likely negligible. However, records of much higher flow velocities during storm events > 1 m/s are recorded the depths of these study areas (Greenan et al., 2010). The high mechanical forcing at these velocities would mean the slowed boundary layer would offer significant protection for both sponges and other organisms associated with sponge grounds (Klitgaard, 1995;Freese and Wing, 2003;Miller et al., 2012).
The most obvious impact of a reduced flow around the sponge on nutrient supply appears to be deleterious: reduced flow would result in a reduced supply of particulate and dissolved organic and inorganic (e.g., dissolved silicon) nutrients to the organism. However, flow speed and nutrient supply may be a more nuanced association due to particulate fallout rate and the active pumping effect of the sponges themselves. The slowing of the boundary layer may result in a greater proportion of the particulates carried by the flow above to be deposited around the sponges. The flow above the boundary layer is not slowed so the vertical supply of particulates to the sponges may increase even if the flow acting on the sponges themselves is reduced. Imposing a sediment flux into the model would be interesting to investigate as a next step in the development of this research.
The reduction in flow speed at the sea floor due to the boundary layer could potentially lead to an increase in larvae deposition and recruitment (Abelson and Denny, 1997). The deposition and retention of larvae in the sponge ground should increase with lower flow speeds, and larvae attachment rates should also increase due to reduced flow stress (Mariani et al., 2006; Guillas et al., 2019). An increase in larvae deposition and retention impacts the recruitment of new individuals within the sponge grounds and should increase its ability to grow and recover from damage.
The importance of ambient flow rate compared to pumping rates for sponges is challenging to determine (Leys et al., 2011;Ludeman et al., 2017). Sponges appear to stop pumping at higher flows but whether this is to protect the canal system from being clogged and damaged (Ludeman et al., 2017;Strehlow et al., 2017), because it is unnecessary for flow though the sponge (Leys et al., 2011), or due to another underlying mechanism, is currently unclear. Regardless, the presence of a slowed boundary layer is influential for each of these mechanisms. A slowed boundary layer may be advantageous as the sponge may be able to pump for longer before it must shut down for protection, or disadvantageous due to higher energy expenditure.
A noticeable feature of the idealized model is the boundary effect of the structures themselves (Figures 4, 5). A future avenue of exploration would be to investigate how this impacts the potential radius of intake around the sponge due to pumping. At what point does the pumping rate become negligible compared to the dividing flow and increased turbulence around the structure, effectively isolating it from the surrounding flow? Does this result in a slower rate of flow becoming advantageous, as the pumping rate with lower turbulence would be more effective and could draw in particulates that would otherwise be lost to the current? Modeling would be particularly useful in investigating the behavior of particles and turbulent currents around sponges, both of which are challenging to measure or observe in situ because of interference from any recording device on the flow.
Regardless of the mechanism by which a slowed boundary layer is created, its presence has additional and important consequences for other organisms associated with sponge grounds (Freese and Wing, 2003;Miller et al., 2012). It also has potential corollaries for conserving sponge grounds under threat of damage by anthropogenic activities. Removal of the sponge structures (e.g., by trawling) would result in the boundary layer velocities not being reduced as much, which-if this feature is beneficial-will have implications for any future recruitment of larvae and recovery of sponges, and so the sponge ground habitat itself.

CONCLUSION
This study shows that there are significant small-scale patterns in distribution of sponges within sponge grounds in the Labrador Sea. Further, we have developed a new approach to studying the drivers of these observed patterns using fluid flow modeling. The models produced allowed detailed quantitative experiments to be conducted, producing data that would be impossible to obtain from field observations alone. By applying these models to sponge structures taken from four study areas in the North Atlantic, we identified a slowed boundary layer caused by sponge aggregations that may account for some of the distribution patterns. This slowed boundary layer may be acting as a protective layer, increasing larval recruitment or promoting increased particulate organic matter uptake due to increases in both particulate fallout from flow and effective radius of pumping. These observations increase the importance of protecting sponge grounds from anthropogenic damage as the removal of this slowed layer coupled with the relatively slow growth rates of sponges in the deep-sea and decreased larval recruitment could result in the inability of sponge ground to recover once damaged. This study highlights the importance of coupling mathematical modeling with field observations in the deep sea to provide a much greater understanding of anthropogenic and climate change impacts on sponge ground ecosystem structure and function, and biogeochemical cycling by sponges.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
TC, JP, and KH conceived the study. TC carried out the study, wrote the manuscript, and prepared the figures. All authors analyzed the data, revised drafts of the manuscript, and gave final approval for the publication.