CFD simulations using our idealized arrays were set up using a hexahedron flow domain (1.75 m × 0.71 m × 0.25 m; W x D x H) that was sufficiently large to reconstruct flow patterns without influence from boundaries. Following previous setups (e.g., Rahman, 2017), our seafloor and Ernietta individuals were assigned no-slip conditions (e.g., fixing the velocity at 0 m/s), the upper face and flow-parallel faces were assigned slip conditions, and the inlet and outlet conditions were assigned on the upstream and downstream faces, respectively. For the inlet boundary condition, a fully-developed, depth-averaged velocity was prescribed (Gibson et al., 2021), and a zero-pressure condition was assigned to the outlet. We chose U = 0.1, 0.2, 0.5 m/s following previous Ernietta CFD analyses (Gibson et al., 2019). Liquid water material properties [density ρ = 1000   k g / m 3 , dynamic viscosity μ = 0.001   k g / ( s ∗ m ) ] were assigned to the flow domain. The domain was meshed using free tetrahedral elements with two layers of additional hexahedral elements proximal to the no-slip conditions to better resolve the viscous sublayer. Because our system Reynolds number was R e = 74,826   (Supplementary 1), turbulent flow was solved using the Spalart-Allmaras one-equation eddy viscosity RANS turbulence model. This turbulence model was chosen to more accurately resolve flow separation from adverse pressure gradients (Bardina et al., 1997; Blazek, 2001), such as between the smaller and larger individuals in our populations.

Recent analyses have demonstrated that LES can produce fluid flow results that are more accurate than the widely-employed RANS approach (Gibson et al., 2021). This improvement stems from solving the Navier-Stokes equations without implementing a Reynolds stress tensor (Blazek, 2001). For this reason, we used COMSOL’s LES residual-based variational multiscale (RBVM) turbulence closure on the Hansburg population CFD simulation. For computational ease, we initiated our LES simulations with a previously computed stationary RANS solution following the protocol presented in Idealized Populations CFD Protocol (Reynolds-Averaged Navier-Stokes). We used a hexahedron flow domain (150 cm × 80 cm × 15 cm; W × D × H) that was sufficiently large to reconstruct flow patterns with minimal influence on the population from boundaries in the initial stationary simulation. Following Idealized Populations CFD Protocol (Reynolds-Averaged Navier-Stokes), we applied the same boundary conditions, fluid material properties, and performed a mesh sensitivity analysis (Supplementary 2). Following the idealized population methods, we used the simpler Spalart-Allmaras turbulence closure to solve the flow fields that were used to speed the convergence of the LES simulation. Using the mesh-independent solution (“Normal”; Supplementary Figure S2), we integrated pressure (p) on the upstream and downstream faces to calculate the pressure differential ( Δ p ) that was required to drive flow at the desired U = 0.5   m / s depth-averaged velocity. As there were no paleocurrent markers preserved within the slab, we chose the direction of flow based off of our hypothesis framework. While flow could feasibly have come from other directions, this orientation was chosen to better reconstruct the flow dynamics of smaller individuals behind larger individuals while simultaneously reconstructing flow over isolated individuals both at the leading edge of the population and downstream (see Supplementary Figure S3). Because LES is so computationally intensive, only the 0.5 m/s depth-averaged velocity was solved for in this study.

Using the same flow domain and mesh procedures described in Idealized Populations CFD Protocol (Reynolds-Averaged Navier-Stokes), we modified our design for the LES analysis by changing the flow-parallel hexahedron slip faces that are orthogonal to the no slip seafloor to be periodic boundaries with a pressure differential Δ p = 0   P a . In a physical sense, this allowed fluid that intersected these boundaries to exit the domain on a face and be reintroduced on the opposing face while preserving momentum and continuity. We then changed the inlet and outlet faces to also be periodic boundaries, but additionally prescribed a pressure differential of Δ p = 2   P a to approximate flow at the desired depth-averaged velocity 0.5 m/s. A pressure point constraint was placed on the no slip seafloor away from Ernietta individuals (see Gibson et al., 2021), which provided a known pressure point for calculating variations within the pressure field. We used the solution derived from the Spalart-Allmaras preliminary simulation as the first time step (e.g., initial conditions, t = 0 ) for the flow field, and ran the simulation for 30 s using COMSOL’s dynamic time stepping to ensure adequate time resolution was solved between time steps (but only output flow fields every 0.01 s). Additionally, due to non-linearities from the fully-coupled flow field turbulent solver, we used an Automatic Fully Non-linear Newton solver instead of the default Constant Newton Iterations solver settings. Finally, we verified our flow fields were in agreement with theoretical constraints by computing the time-averaged velocity profile throughout multiple points in the flow domain. As described in theory (Prandtl, 1905; Gibson et al., 2021), these profiles follow the logarithmic law of the wall.

As shown in Figure 2, both the small and large Ernietta individual flow patterns show patterns consistent with CFD descriptions elsewhere (Gibson et al., 2019; 2021). In the small individual simulations, flow enters the cavity from above, is directed to the downstream interior face, and is then forced upstream toward the interior and upstream face, where it is finally directed out of the cavity by mixing with the current above the cavity. A wake of slower, more turbulent water is formed directly downstream of the individual. These patterns are stronger or weaker with faster and slower far field velocities, respectively. The larger individual developed the same patterns of recirculation within the cavity and downstream with a similar relationship between recirculation speed and far field velocity, but the velocities and mixing were greater than those experienced by the smaller individual.

Figure 3 show the results for each of the idealized population CFD studies. In all simulations, there is recirculation in both large and small individuals in the population simulation, though absolute velocities vary depending on size, location, and proximity to neighbouring individuals. Cavity velocities are always slower than the far field velocities, and there is a higher-pressure region upstream of the population where fluid is directed around and over the aggregate. A wake is also formed downstream of the population in all simulations.

Figure 4 and Supplementary Figure S4 show that flow patterns within the reconstructed Hansburg slab broadly follow those described for idealized arrays, with some key distinctions. All individuals receive recirculation within their cavities, and these cavity velocities are reduced compared to the far field velocity. While there is consistent mixing within cavities, most fluid is preferentially directed around and above individuals. Slower velocities are observed in upstream individuals, while faster velocities are observed within the aggregations and within the downstream individuals’ cavities. Smaller Ernietta within the aggregate that are located directly downstream of the largest individual consistently receive fast recirculation–faster than even larger individuals located closer to the upstream face of the largest individual (Figure 4; Supplementary Figure S4). Narrow individuals have slower cavity velocities. Eddies are shed off of all individuals, but orientation to flow, individual shape, and individual location dictate how far downstream eddies persist. Narrow individuals that are oriented broad-side into flow typically exhibit rapid mixing, including upstream recirculation directly behind the individual, while the wakes of those individuals that form part of aggregates are typically disrupted by wakes of neighbouring individuals, leading to additional turbulence. Isolated individuals typically develop small vortex streets that are disrupted by individuals further downstream.

Previous CFD analyses on Ernietta individuals and idealized populations investigated how burial depth and individual spacing may have controlled feeding efficiency (Gibson et al., 2019). This work demonstrated that turbulence increased above downstream individuals in addition to generating faster recirculation within the cavities of the same individuals (see Gibson et al., 2019 Figure 3), suggesting that the faciliatory lifestyle exhibited by Ernietta may have imparted a selective advantage for feeding and waste dispersion. While the new CFD results targeting individuals of mixed sizes presented here broadly support the previous results regarding cavity recirculation and suspension feeding, these new data provide additional insights that were not available in the previous analyses. Recirculation is faster within the cavity of larger individuals than in smaller ones, regardless of far field velocity. This indicates that a greater volume of fluid should cycle through larger individuals, potentially delivering larger quantities of food. Also, these new idealized population simulation results partially agree with previous work suggesting that mixing is increased above and within the cavities of downstream individuals. This agreement is strongest for simulations where the larger individuals are downstream of the smaller ones. When smaller individuals were upstream of larger individuals, they received slower velocities than those that were placed downstream (Figure 3). Interestingly, our results demonstrate that the locations of smaller individuals in idealized populations of Ernietta directly affect flow patterns, and thus impact our interpretation of their feeding efficiency. In idealized populations constructed with smaller individuals downstream of larger individuals (Figure 3), smaller individuals experience slower cavity recirculation than the larger individuals (contra Gibson et al., 2019). Based on these comparisons, this would signal that, while vertical mixing is enhanced above aggregations of similarly sized Ernietta, those individuals that are smaller and downstream of larger ones might be flow starved compared to the larger individuals. These idealized populations would then suggest that it might be optimal for Ernietta to form populations of similarly-sized individuals, (i.e., as single cohorts that result from reproduction with spat-fall, however, see below).

Previous studies have noted that individual Ernietta fossils can vary in shape, with authors attributing these shape variations to a combination of individual growth response to the local environment and/or taphonomy (Elliott et al., 2016; Ivantsov et al., 2016; Smith et al., 2017). In our idealized population models, individuals were identical in shape but were binarily scaled. While cavity velocity magnitudes differed depending on location, the overall recirculation patterns were still comparable between simulations. In contrast, our Hansburg individuals were asymmetrically scaled to better approximate the spatial reconstruction of the Farm Hansburg in-situ population (Supplementary Figure S1). Those results demonstrate cavity and downstream flow patterns are strongly affected by overall individual shape (e.g., slender vs round). Slender individuals receive slower and more asymmetric recirculation compared to rotund individuals. These patterns were present within aggregated and isolated individuals, indicating that overall shape exerts a strong control over internal cavity flow. There may, therefore, have been strong selective pressure for individuals to grow with a more uniform cavity diameter, as evident in many Ernietta deposits (see figures and descriptions from Jenkins et al., 1981; Ivantsov et al., 2016; Elliott et al., 2016; Smith et al., 2017), while cavity filling by deposited sediment would serve to retain the overall oval to circular shape of the cavity.

Our results also show that an individual’s location (or proximity) to neighbors is an equally important control on the generation of flow patterns thought to be important for feeding. Isolated individuals within our Hansburg population experience slower cavity recirculation even when having a more uniform cavity diameter (Figure 4; Supplementary Figure S4), while aggregated individuals experience comparatively faster cavity recirculation. Flow patterns around our individual Ernietta and idealized arrays (Individuals and Idealized Populations vs Previous Work) suggest that smaller individuals experience slower cavity velocities when downstream of larger individuals. However, our Hansburg population results demonstrate that even when smaller individuals were placed immediately downstream of larger individuals, the smaller individuals were constantly exposed to faster velocity flow when compared to isolated forms. At first impression this appears contradictory, but it may highlight the importance of modeling natural systems as opposed to overly simplified (i.e., idealized) ones. Furthermore, our reconstructed population simulations agree with previous results that demonstrated increased mixing above and within downstream cavities (Gibson et al., 2019). This suggests that smaller individuals within aggregates are not hampered when behind a larger individual, but in fact are likely to receive more and better-mixed fluid flow. Thus, while it is less advantageous for an individual to be near the edge of the aggregation, smaller individuals would likely gain access to additional resources that they might not otherwise receive should they be on the edge of the population, or entirely isolated. While little is known about the ontogenetic history of Ernietta, it is reasonable to assume that younger individuals would be smaller; our results suggest that there would thus have been an advantage to settling, or living near, larger individuals within the greater population. In this fashion, aggregations of adult Ernietta would likely have performed a “nursery” function, creating localized conditions favoring the settlement and growth of younger individuals.

This complex interplay of controls on individual feeding success has some parallels in extant benthic marine ecosystems. Studies focused on intraspecific competition arising from spatial density within populations of sessile suspension feeders demonstrate that, in multi-sized aggregations, intraspecific competition can affect both individual size and spatial distribution (Werner and Gilliam, 1984). When different sizes preferentially exploit disparate resources, the competition between adults and juveniles is reduced (Werner and Gilliam, 1984). However, in such cases where exploited resources overlap as is plausible with Ernietta–and size does determine the outcome of competition, younger individuals are often less successful due to metabolic constraints (Peters, 1983). Importantly, the magnitude of this general relationship is greatly influenced by resource availability. While a paucity of resources does typically lead to increased competition, when resources are abundant this can decrease levels of exploitative competition (Hazlett et al., 1975; Werner and Gilliam, 1984). This distinction is important when applying these models to Ediacaran-aged communities in Namibia; recent geochemical work on the Nama Group suggests large scale unstable redox conditions accompanied by decreasing nutrient delivery to inner ramp habitats where Ernietta are typically found (Bowyer et al., 2017; 2020). In combination with the evidence presented here for relatively low-levels of intraspecific competition, this would indicate that populations of Ernietta on Farm Hansburg either: 1) locally received short-lived, adequate nutrient levels to offset multi-sized intraspecific competition; or, 2) Ernietta was an ecological generalist not endemic to oxygenated, nutrient-saturated environments (Darroch et al., 2018b).

While our simulations were designed as an extension of recent work (Gibson et al., 2019, 2021) with the goal of moving towards reconstructing fluid flow over spatially-accurate 3D populations of Ediacaran taxa, they do allow us to evaluate other hypotheses regarding the paleobiology and life habit of Ernietta. For example, Bouougri et al. (2011) suggested that living Ernietta may have been covered in a slime or mucus layer, based on the orientation of sediment detrital mica flakes associated with fossil molds. While our results cannot support or refute this suggestion, our simulations do at least show that recirculation in the wakes of individual Ernietta would provide one possible mechanism for fine-grained sediment (presumably together with suspended food particles) delivery to downstream surfaces. Given that many modern suspension feeding organisms make use of mucus nets and other structures as a means for capturing food particles (see discussion of aerosol capture mechanisms in Rubenstein and Koehl, 1977; and LaBarbera, 1984), we suggest that an external slime- or mucus-covered surface could have been used to capture suspended particulate material, and which could then have been sorted and/or transported to other feeding organs (see, for example, LaBarbera, 1978). However, without additional evidence for mucus membranes, or indeed any other evidence revealing how food particles are ingested and processed, this model remains highly speculative.

In a similar vein, Ivantsov et al. (2016) suggested that the “apical fan” of Ernietta (i.e., the portion of the organism raised into the water column) may have been used for respiration. While no previous CFD studies have interpreted this anatomy in the context of respiration, many of the flow patterns modelled around Ernietta could have been useful for bringing oxygen-rich water closer to respiratory organs. The increased turbulence modelled in these areas could plausibly have served to replenish dissolved gases, as well as diffuse any exhalents. Although we cannot directly address hypotheses surrounding respiration without substantially adjusting the methods employed here, future coupling of fluid dynamics simulations with oxygen diffusion/uptake models (see, for example, Ghisalberti et al., 2014) may offer an elegant means to test this idea. Future work using this combination of methods could plausibly be used to test the hypothesis that Ernietta was opportunistic, and able to colonize shallow marine environments that were only transiently oxygenated (see Hall et al., 2013; Bowyer et al., 2020).

Previous paleoecological work focused on the latest Ediacaran Nama interval has suggested that communities of soft-bodied Ediacara biota were relatively depauperate, possessing lower species richness than older Ediacaran communities, as well as having relatively “simple” structures consistent with either competition for a limited number of ecological resources and/or ecological stress (Darroch et al., 2015; 2018a,b). This study adds nuance to this picture, illustrating that, while species poor and relatively simple, the structure of Nama-aged erniettomorph communities was likely influenced by a wide-range of intraspecific ecological factors, including size, location, and (potentially) developmental stage. In this study, we show that individual size and location within a population can have large effects on flow within and around the body cavities, affecting their access to nutrients and well-mixed water. Given the lack of any evidence for asexual reproductive modes such as budding in Ernietta, it is reasonable to hypothesize that these factors might have been crucial during the recruitment and settlement of new individuals. In addition, while not directly assessed here, our flow patterns showing increased fluid delivery downstream of larger individuals also suggests a potential ‘nursery’ effect, with competition for space adjacent to larger Ernietta. Many of these density-dependent factors continue to influence the structure of communities in today’s oceans (Caley et al., 1996; Gilmour, 1999; Harrington et al., 2004), and so reinforce the point that, while we still have little idea where taxa like Ernietta fit (if at all) in the metazoan tree of life, in many respects they experienced similar ecological pressures–and may have utilized similar ecological strategies–as modern benthic marine organisms.