Morphotype differentiation in the Great Barrier Reef Halimeda bioherm carbonate factory: Internal architecture and surface geomorphometrics

The calcareous Halimeda bioherms of the northern Great Barrier Reef, Australia are the largest actively accumulating Halimeda deposits worldwide. They contribute a substantial component of the Great Barrier Reef neritic carbonate factory, as well as the geomorphological development of Australia's northeast continental shelf. Halimeda bioherm geomorphology is complex, expressing three distinct variations in morphotype patterns: annulate, reticulate and undulate. Similar regular and irregular geomorphological patterning often results from scale‐dependent biophysical feedback mechanisms. Therefore, a better understanding of morphotype differentiation can inform the biotic and abiotic drivers of spatial heterogeneity in the bioherm ecosystem. Here, 3D LiDAR bathymetry is integrated with 2D sub‐bottom profile datasets to investigate surface topography and internal sedimentary architecture of Halimeda bioherms through space and time. Using the ESRI ArcGIS 3D Analyst and Benthic Terrain Modeller extensions, the bioherm surface and subsurface geomorphometric characteristics were quantified for the annulate, reticulate and undulate morphotypes. Significant variation was found between the three bioherm morphotypes in their surface topography, internal structure, volume, slope gradients and terrain complexity. Therefore, their geomorphology is probably influenced by differing processes and biophysical feedback mechanisms. The complex surface topography does not appear to be inherited from the antecedent substrate, and preferred aspect orientations resulting from hydrodynamic forcing appear to be limited. It is suggested here that autogenic dynamics or biotic self‐organization similar to patterns and processes in other marine organo‐sedimentary systems modulates Halimeda bioherm geomorphology, and some hypotheses are offered towards future studies. Morphotype differentiation has implications for the development of the Halimeda bioherm carbonate factory, rates of sediment aggradation and progradation, and variable capacity to fill accommodation space. Self‐organization dynamics and morphology differentiation in Modern bioherm systems could potentially inform palaeo‐environmental interpretations of fossil bioherms and phylloid algal mounds on geological timescales.

and sedimentological surveys in the GBR, Indonesia and the Caribbean first revealed their surface topography and internal architecture. Descriptions of Holocene Halimeda bioherms come from the GBR (Davies and Marshall, 1985;Orme, 1985;Drew and Abel, 1988;Marshall and Davies, 1988;Orme and Salama, 1988;Searle and Flood, 1988); Kalukalukuang Bank (K-Bank), Indonesia (Roberts et al., 1987;Phipps and Roberts, 1988); Miskito Channel in the Nicaraguan Rise, south-west Caribbean Sea (Hine et al., 1988); and reported but not described in detail at Big Bank Shoals on the Sahul Shelf, Timor Sea (Heyward et al., 1997) (Figure 2; Table 1). On the bioherm tops, the seabed is consistently around 20-40 m depth, although the Miskito Channel bioherms are deeper at 40-50 m (Hine et al., 1988). The build-up of Halimeda sediment can be up to 50 m thick at K-Bank , but 10-30 m thickness is more common at other locations. These bioherm sites are all located on shelf-edge or platform-top positions, with a steep near-vertical slope dropping to an oceanic trough or basin around 1,000 m deep (Table 1). A common interpretation of formation involves the upwelling of cool, nutrient-rich water from below the oceanic thermocline delivering the volume of nutrients required to fuel Halimeda growth, not otherwise available from the oligotrophic tropical waters (Wolanski et al., 1988). These consistent characteristics make modern Halimeda bioherms an important analogue for the depositional environment, facies interpretation and processes of formation of lithified fossil Halimeda deposits (discussed in the following section), and potentially phylloid algal mounds, from the geological past.

| Global fossil Halimeda outcrops
Globally, several fossil Halimeda bioherm outcrops have been identified and described ( Figure 2; Table 2). The oldest reported example dates back to the Palaeocene and may be the earliest expression of a metre-scale Halimeda limestone deposit since this genus first appeared in the Cretaceous. Dragastan and Herbig (2007) describe Palaeocene and Eocene carbonate ramp deposits from the Atlas Mountains (Morocco), comprised almost exclusively of well-preserved Halimeda segments. Largely a fossil taxonomic study, the authors do not describe any mound-shaped geomorphology, and so stop short of naming these outcrops as Halimeda | 179 MCNEIL Et aL.
bioherms. However, they do make an analogous comparison to the modern GBR bioherms.
Perhaps the most spectacular and recognizable lithified fossil Halimeda bioherm deposits are those from the upper Miocene in the Sorbas Basin (Spain) (Braga et al., 1996;Martin et al., 1997), and second from Salento Peninsula (Southern Italy) (Bosellini et al., 2002). The lithology and geometry at these two sites are both described as in-situ midslope Halimeda packstones to floatstones/rudstones, comprising 80%-90% Halimeda segments in a carbonate mud matrix, forming discrete lenticular mounds and associated with underlying or adjacent Porites coral reefs. In the Sorbas Basin, Braga et al. (1996) describe mounds up to 40 m thick and 400 m long, suggesting that these deposits were at least as large, if not larger, than the modern examples of individual mound structures from the GBR and Indonesia. The coral reef facies association is remarkably similar to modern settings. The erosional discontinuity between the Porites reef and Halimeda bioherms described from Sorbas Basin is directly analogous to the GBR example, where Holocene Halimeda bioherms overly a Late Pleistocene sea-level lowstand erosional unconformity. The GBR underlying Pleistocene surface (antecedent substrate from herein) is comprised of a leached coralgal skeletal packstone/grainstone (Davies and Marshall, 1985).
In the tropical Indo-Pacific, early geologists described Halimeda-rich Pleistocene 'limestone reef' outcrops from the Solomon Islands (Guppy, 1887), Funafuti (Tuvalu) (Hinde, 1904) and New Hebrides (Vanuatu) (Chapman and Mawson, 1906), comprised almost entirely of well-preserved Halimeda  Orme et al. (1978); Davies and Marshall (1985); Drew and Abel (1985); Orme, (1985); Drew and Abel, (1988); Marshall and Davies (1988); Orme and Salama (1988) segments ( Table 2). The Funafuti Atoll coral boring geological expedition was embarked upon to test Darwin's theory of coral atoll formation (Bonney, 1898). Subsequently, Hinde (1904) described a 24 m thick Halimeda packstone unit from Funafuti. Chapman and Mawson (1906) estimated the oldest of the New Hebrides Halimeda limestones as possible Pliocene, but details of relative dating methods are not clear from the available records. Late Quaternary low mounds (2-4 m thick) are also described from Fifty Fathom Flat off the western coast of India (Rao et al., 1994; on what is now a drowned carbonate platform.

bioherms of the GBR
The GBR bioherms are the most extensive, actively accumulating Halimeda deposits in the world (Whiteway et al., 2013). They form part of the GBR World Heritage Area, and their size and extent contribute to the GBR World Heritage Outstanding Universal Values from a geological and geomorphological perspective. The bioherm complexes extend over six degrees of latitude from 10°S to 16°S and cover an area greater than 6,000 km 2 (McNeil et al., 2016) (Figure 3). Their calcium carbonate volume is estimated to be up to four times greater than that of the adjacent Holocene coral reefs (Rees et al., 2007). In addition to the expansive northern GBR bioherms, a much smaller area with the same topographic expression has been identified in the southern GBR Swain Reefs region (Searle and Flood, 1988;McNeil et al., 2016) ( Figure 3B). Additionally, seismic profiles (Hinestrosa et al., 2014) and dredged samples of Halimeda floatstone (Abbey et al., 2011) from submerged (up to 130 m depth) shelf-edge reef terraces suggest evidence of Halimeda-rich mound development during Late Pleistocene sea-level oscillations (Abbey et al., 2011;Hinestrosa et al., 2014). In the GBR, the previously known extent and geomorphological descriptions were based on sediment grabs, vibracoring, and widely spaced 2D seismic profiling and singlebeam echosounder profiles (Searle and Flood, 1988;Hopley et al., 2007;Rees et al., 2007;Davies, 2011). The first model of Halimeda bioherm morphology was constructed at the time, by a reasonable interpolation between 2D parallel seismic and singlebeam echosounder profiles, resulting in the elongate parallel ridges and troughs, and lenticular mound-like morphology descriptions that permeated the literature (Orme and Salama, 1988;Mathews et al., 2007). Previously proposed mechanisms to explain the GBR bioherm morphology include the merging of smaller individual lens-shaped mounds (Marshall and Davies, 1988;Searle and Flood, 1988); stabilization of in-situ accumulations by trapping and baffling by surface flora (Marshall and Davies, 1988); and hydrodynamic processes involving tidal jets and back-reef eddies (Marshall and Davies, 1988;Wolanski et al., 1988). Paradoxically, sediment textural descriptions from core tops suggest little to no winnowing or sorting action by currents (Orme, 1985). Searle and Flood (1988) proposed that bioherm morphology is not inherited from the antecedent substrate based on the planar nature of the underlying surface in the southern GBR. More recently, descriptions of the spatial distribution and morphology of the bioherms in the northern GBR were re-assessed (McNeil et al., 2016) with the benefit of modern 3D survey technologies such as Laser Airborne Depth Sounding (LADS) bathymetric LiDAR, and multibeam echosounders (Beaman, 2010). The resulting comprehensive new bathymetry dataset revealed that bioherm shapes and patterns were not consistent with the previous parallel ridges and troughs representation. McNeil et al. (2016) identified and described three morphological subtypes, termed reticulate, annulate and undulate (Table 3). The reticulate morphotype has a complex honeycomb-like topography reminiscent of coral patch reef morphology, sometimes described as 'cellular' reefs (Purkis et al., 2010;Blanchon, 2011;Blakeway and Hamblin, 2015;Schlager and Purkis, 2015). Perhaps the most striking recent discovery is the annulate morphotype. These 'mounds' appear lenticular or lens-like in 2D profile, but in fact form hollow centred circular ring shapes with elevated crests when viewed in 3D (Video S1). Each annulate ring is 250-300 m in diameter from crest to crest. Adjacent rings sometimes coalesce together sharing side walls (McNeil et al., 2016). The reticulate morphotype tends to be distributed towards the eastern (seaward) bioherm boundary, proximal to the adjacent coral barrier reefs and continental shelf edge that lies 3-4 km to the east. The annulate morphotype is situated more distal to the shelf edge, and annulate rings generally increase in size to westward (away from the shelf break). The undulate morphotype has lower relief, is smoother and more sinuous, and grades into and around the reticulate and annulate types (McNeil et al., 2016). The true nature of bioherm topographical and morphological expression is only revealed when viewed in 3D, and is consistent throughout the entire >6,000 km 2 distribution.
The recent discovery that the bioherms are much more extensive and morphologically complex than previously thought requires a re-evaluation of existing models describing their origin, growth and development (McNeil et al., 2016). In particular, the processes (e.g. depositional, erosional, hydrodynamic, biological) driving the development of this unusual geomorphology are not understood. It is not clear how the three different morphotypes are developed, particularly the circular annulate rings, although McNeil et al. (2016) postulated that biotic self-organization may play a role. Geomorphological differences between the three types, their relationships spatially, and with the underlying antecedent topography have not been assessed. Any future F I G U R E 3 Regional map of study area. (A) Northern Great Barrier Reef Halimeda bioherm distribution and locations of figures referred to in this study; (B) southern Great Barrier Reef Halimeda bioherm distribution in the Swain Reefs. Halimeda bioherm distribution from McNeil et al. (2016) analysis of morphological control mechanisms first requires a comprehensive understanding of the physical characteristics of each morphotype.
This study presents a quantitative geomorphological analysis of three different Halimeda bioherm morphotypes by integrating the most recently available geophysical datasets to address three main objectives: (a) Combine high-resolution 3D surface bathymetry with co-located acoustic sub-bottom profiles and GIS spatial and morphometric terrain analyses; (b) Quantify the variation in physical characteristics such as bioherm area, thickness, volume and terrain patterns for the annulate, reticulate and undulate morphotypes; (c) Better constrain the size of the Halimeda bioherm carbonate factory by quantifying the mass and volume of accumulated CaCO 3 by each morphotype; and (d) Use these analyses to determine the differences and similarities between the annulate, reticulate and undulate morphotypes, and resolve their sedimentological relationships. These analyses will provide crucial new information towards understanding the mechanisms controlling bioherm morphology and potential for biological self-organization.

| Surface geomorphic analysis
A GIS geomorphometric analysis was undertaken on the most recently available bathymetric digital elevation model (DEM) of the northern GBR, that was compiled as part of Project 3D-GBR (Beaman, 2010;2017) and available for download from the AusSeabed Marine Data Discovery portal (http://marine.ga.gov.au/#/). The input datasets used in this study are described in Table S1, with links to open access data repositories where available. The bathymetric DEM is a compiled dataset from all available multibeam echosounder data within the GBR from numerous scientific research voyages, and hydrographic surveys including airborne bathymetric LiDAR acquired by the Australian Hydrographic Office (AHO). The source bathymetry data were collected using a WGS84 horizontal datum and lowest astronomical tide (LAT) vertical datum, with point spacing varying between 6 and 30 m depending on the survey method. A vertical adjustment between LAT and MSL (mean sea level) was made to give an approximate MSL vertical datum across the dataset. The compiled DEM was generated with a 0.0003 arc degree (30 m) pixel size, and then converted to an ESRI raster grid (Beaman, 2017) for spatial analysis in ArcGIS 10.5 software.
This study concentrates on the Halimeda bioherm complexes in the northern GBR ( Figure 3A). A small area of bioherms is present in the Swain Reefs complex in the southern GBR ( Figure 3B). However, mapping in this region is incomplete and the much larger northern section provides a statistically more robust dataset.

| Morphometric analysis
From the >6,000 km 2 bioherm distribution, two target areas of interest ( Figure 3) were selected for detailed analysis that met the following criteria: the area is mapped by a single survey type (the LADS LiDAR bathymetry) for consistency of resolution from the original point data imported to the compiled DEM; covering the largest available area of all three morphotypes; and crossed by Topas sub-bottom profiles running normal (E-W) and parallel (N-S) to the shelf edge (see Section 2.3). Using the 3D Analyst toolbox and Benthic Terrain Modeler (BTM) extensions in ESRI ArcGIS 10.5 with the gridded bathymetry DEM as input data, a set of morphometric derivative outputs were generated to visualize, describe, quantify and compare the surface characteristics across a bioherm complex. These types of terrain attributes can link to physical and biological explanatory processes (Lecours et al., 2016;Walbridge et al., 2018). The geomorphometry derivatives presented in this study for target area 1 (e.g. slope, aspect, curvature, terrain ruggedness.) are comprehensively described by Walbridge et al. (2018), with mathematical algorithms and extended references, including their workflow through the ArcGIS BTM software. These methods, and the context to which they have been applied, are briefly described herein. Surface gradients (slope, aspect, curvature) calculated from the pixel-pixel relationships in bathymetric data form the foundation of the morphometric analyses (Lecours et al., 2016). Two slope derivatives are presented, one on the planar cartesian surface and an arc slope calculated on the geodesic ellipsoid. The arc slope method accounts for the angle between the input surface and the geodetic datum, applies a least squares fitting, and can produce more accurate results when applied to high resolution bathymetry with low positional uncertainty (Passalacqua et al., 2015). The two slope outputs are displayed with different symbology to highlight features of interest. To compare the annulate morphotype inner slope gradient with the outer slope gradient, a subset of annular rings were manually selected within target area 1 and used to calculate the average slope. The largest six rings were selected because larger rings are represented by a greater number of pixels (n = 206 pixels). Rings were cross-checked against the aspect azimuth to confirm inward-facing from outward-facing pixels, then the slope gradient for each inward and outward-facing pixel was recorded. Profile curvature calculates the change in slope across bioherm features and the adjacent non-bioherm seafloor. Application of the aspect computation generates a statistical visualization of surface direction, or 'easterness' and 'northerness' (Walbridge et al., 2018) in azimuth degrees for each cell.
Bathymetric position index (BPI) is a measure of relative position (elevation) between any cell and the overall seascape, by computing the difference between any cell and the mean elevation of all cells within a surrounding annulus-shaped neighbourhood (Lundblad et al., 2006). Usually computed over both broad and fine scales, only the broad-scale BPI is presented as this provided the most appropriate scale from the input data resolution and the size of the target area of interest. The default neighbourhood radii values of one and three pixels were used for the inner and outer annulus circles respectively. Resulting values are positive near crests and ridges, and negative near depressions and valley bottoms.
The ArcGIS BTM provides four methods across three tools to quantify terrain complexity or surface heterogeneity: Surface Area to Planar Area (SAPA) and slope corrected SAPA, Vector Ruggedness Measure (VRM) (previously termed Rugosity in BTM), and Arc-chord ratio. The VRM described by Sappington, Longshore and Thompson (2007) was applied in this analysis. The VRM calculates a dimensionless metric of terrain complexity using a moving 3 × 3 cell window, where values range from 0 (no terrain variation) to 1 (complete variation) (Walbridge et al., 2018). All of these geomorphometric outputs and numeric derivatives are scale and resolution dependent, and the process of creating the input DEM surface (30 m gridded resolution) will affect the accuracy and error margin of the outputs. For example, a single annulate ring of 300 m diameter would be covered by a grid of 10 × 10 pixels. A 3 × 3 cell window in the BTM morphometry computations was found to be the most appropriate scale for the input data.

| Subsurface geomorphic analysis
Subsurface acoustic data were collected using a Kongsberg TOPAS PS-18 sub-bottom profiler (frequency 1,000 to 6,000 Hz, 1,500 Hz central frequency) on the RV Southern Surveyor voyage SS09/2008 (Tilbrook and Matear, 2008). Over 210 km of seismic data in seg format were viewed in IHS Markit Kingdom seismic interpretation software (2017 version). To calculate the Halimeda bioherm thickness, the seafloor horizon (bioherm top surface) and the first prominent seismic reflector (bioherm bottom surface) were delineated. This first reflector was termed Reflector A following the notation of Orme et al. (1985) who interpreted this reflector as the Pleistocene erosional unconformity that marks the antecedent substrate of the Holocene Halimeda bioherms. Reflector A was delineated using a combination of automated and manual horizon picking to avoid interpreting multiples (seismic artefacts). The isochron thickness between these two horizons was then computed in two-way travel time (TWT; ms) with these xyz data exported as a text file. The isochron xyz data were imported into ArcGIS and clipped using the Halimeda bioherm shapefile extent (McNeil et al., 2016), to ensure that only areas confidently identified as Halimeda bioherms were included in subsequent calculations. To reduce the margin of error of the thickness calculation, any discontinuous or uncertain sections of Reflector A were excluded from the dataset, such that the isochron thickness, and subsequently the bioherm volume, was only calculated on sections where Reflector A could be confidently delineated.
Statistical calculations (ANOVA) of the differences in thickness between the reticulate, annulate and undulate morphotypes were calculated on TWT. To present a meaningful description of thickness and for use in bioherm volume calculations, TWT was converted to an estimate of thickness in metres by applying the best estimate of Halimeda sediment acoustic velocity from the literature. Skjold (1988) undertook a continuous seismic reflection (boomer, dominant frequency 3-4 kHz) and vibracoring survey close to the target area of interest 1. The seismic acoustic parameters and sediment physical properties from that survey are compiled in Table 4. Average acoustic velocity in the Halimeda lithofacies is 1,532 ± 19 m/s (±SE) calculated from 24 measurements across five Halimeda vibracores. These short cores (max 3.8 m) only capture the top section of the full bioherm sequence. Therefore, acoustic velocity may vary downcore with increased compaction or changes in sediment texture. However, these are the only seismic velocity data currently available for any GBR survey that specifically targets Halimeda bioherms, and are consistent with the summary of carbonate Quaternary P-velocities summarized by Hinestrosa et al. (2014). The bioherm volume (km 3 ) was calculated for each morphotype by multiplying average thickness of each morphotype (m) by planar area (m 2 ) and converting to km 3 . Furthermore, the grain density, porosity and CaCO 3 % measured by Skjold (1988) (Table 4) were used to calculate the volume and mass of CaCO 3 by each morphotype by adapting the equations described in Rees et al. (2007).
Statistical analyses to compare differences in measured morphometric parameters between the annulate, reticulate and undulate morphotypes were undertaken by analysis of variance (Kruskal-Wallis one-way ANOVA). The H-statistic tested the hypothesis of no statistically significant difference in median values between the three morphotypes.

| Surface geomorphometry and quantitative morphometrics
Results presented here are from two representative target areas of Halimeda bioherm morphology that express the topography and spatial organization of the three morphotypes, annulate, reticulate and undulate, from over 210 km of sub-bottom profiles combined with high-resolution LiDAR bathymetry.
The annulate morphology is characterized by sedimentary packages with high positive relief that are relatively uniform in thickness except where punctuated by annular rings. The central hollow of the rings plunge to the deepest section, at times almost reaching the antecedent substrate ( Figure 4A). The annulate surface is sinuous and relatively smooth between individual rings, compared with the reticulate surface ( Figure 4C). The crests of each ring are gently convexly curved, and the inward-facing slope is consistently steeper than the outer slope ( Figure 4E). Calculations based on the ArcGIS BTM outputs indicate that the inner slope gradient is approximately three times steeper than the outer slope (Table S2 and Figure 5C). The antecedent substrate (Reflector A) is prominent and some internal bedding is evident in the bioherm. The topmost internal bedding reflectors are generally parallel to sub-parallel to the seafloor, following the upwardly convex curvature of the mounded surfaces, and generally downlapping on underlying surfaces. These internal reflectors are indicating metre-scale features and it is not entirely certain what this bedding represents.
The reticulate morphology ( Figure 4C) is considerably more complex, with steeper topographic relief and sharper crests than the annulate rings. There is more variability in bioherm thickness, but the valleys and depressions seldom reach the antecedent substrate. The internal reflectors are also more complicated, at times showing the convex-up bedding similar to the annulate morphotype, but also showing evidence of abrupt truncations (e.g., Figure 4D), and multiple generations of growth or deposition.
The undulate morphotype ( Figure 4B) typically shows very little topographic relief. The surface is smooth to wavey, and internal bedding is parallel to both the seabed and antecedent substrate. There are two prominent 'basement' reflectors in the seismic section presented in Figure 4B. It is not certain whether the first prominent reflector is an internal bedding surface and the second reflector is Reflector A (Holocene/Pleistocene boundary), or whether the first reflector is Reflector A and the second is a pre-Holocene surface. The paucity of available core material limits cross-referencing of seismic stratigraphic features against lithostratigraphy in most areas.

MCNEIL Et aL.
The results from the ArcGIS BTM geomorphometry analysis ( Figure 5) reinforce the heterogeneity observed from 2D seismic profiles, and add further insights to the characteristics of the annulate and reticulate morphotypes. The reticulate morphotype is most common on the eastern (seaward) side of the bioherm complex, proximal to the outer barrier coral ribbon reefs and continental shelf edge and forms the eastern bioherm boundary ( Figure 5A). The annulate morphotype lies westward, further from the shelf edge and grading out to the western boundary. Coverage of the undulate morphotype increases to the west, away from the shelf edge. It is organized between and around the other two types, grading out to the surrounding inter-reef sediment.
The bioherm boundary is marked by a clear difference in topographic relief between the bioherm tops and the adjacent outer-shelf seafloor. There is an approximately 30 m difference in bathymetry between the outer-shelf seafloor and the western bioherm boundary ( Figure 5B). The tops of both the annulate and reticulate morphotypes are relatively consistent with each other, with the exception of the annulate ring hollows that plunge to the depth of the antecedent substrate. Table 5 shows the combined average bathymetry from target area 1 and target area 2, for each morphotype. The slope output highlights striking differences between the three morphotypes ( Figure 5C). First, the annulate rings clearly show a bi-modal slope distribution, expressed as a double 'halo' effect. The inner, brighter circle shows the steeper slope of the  Figure 3 for location. Inset (D) internal reflectors indicate some truncated erosional surfaces in the reticulate morphotype; (E) annulate ring crests and hollow showing proximity to the antecedent substrate and differentiation between inward-facing and outward-facing slope gradients. Ar = annulate ring. Reflector A is interpreted as the Holocene/Pleistocene erosional unconformity (Davies and Marshall, 1985;Orme, 1985;Orme and Salama, 1988). An acoustic velocity of 1,532 m/s (Table 4) was used to convert from time (TWT ms) to depth domain (m below sea level [m.b.s.l.]) on sub-bottom profiles inward-facing walls and the outer circle shows a more moderate slope on the outward-facing flanks of each annular ring. The total average annulate slope is 1.88° overall (Table 5) but the individual ring outer slope is 2.87° vs the inner slope of 8.92°; the inner slopes are approximately three times steeper than the outer slopes (Table S2). The low gradient of the annulate ring crests (black circle between the bright inner and outer circles, Figure 5C) illustrates the smooth curvature of the annulate ring crests compared to the much sharper nature of the reticulate crests. The reticulate slope inset (Figure F I G U R E 5 Morphotype and morphometric characteristics of area of interest Target area 1 (see Figure 3). Red insets highlight the annulate and reticulate morphotype surface characteristics in a 2.5 km x 2.5 km clip of each type 5C,D) highlights the very different surface topography of this morphotype compared to the annulate. Average slope is 2.71° overall (Table 5), steeper than the average annulate gradient, but lacking the double slope expression of the annulate rings. The difference in spatial organization and surface topography between these two morphotypes is remarkable. The coloured symbology of the arc slope morphometry ( Figure 5D) displays the same data but highlights: (a) the almost flat gradient of the adjacent shelf-edge seafloor and the undulate morphotype (which has an average slope of 1.38°; Table 5); and (b) the abrupt change in bathymetry is expressed by the steep slope at the bioherm eastern boundary ( Figure 5B). This is clearly illustrated in the curvature output, showing the change in slope ( Figure 5E). Here, the change in slope (curvature) at the bioherm boundary is very high (coloured red), showing that the topography steepens sharply from the seafloor to the bioherm surface, rather than a gradual gradient.
The aspect morphometry ( Figure 5F) illustrates that both the annulate and reticulate azimuth orientations are similarly random in all directions and do not appear to express any preferred orientations in 3D space.
The preceding surface gradient morphometrics combine to form the basis of the overall terrain complexity, illustrated by the VRM and the BPI. The BPI output ( Figure  5G) shows the relative differences in elevation between any cell and the surrounding 3 × 3 neighbourhood, such that the purple areas show neutral BPI and coloured pixels show relative positive and negative bathymetric positions. This result shows that the shelf-edge seafloor and the undulate morphotype share similarly flat terrain relative to the surrounding seascape, even though they are at different bathymetric depths. In contrast, the reticulate morphotype shows irregular and variable BPI compared to the regularly spaced patterning in the annulate morphotype. Similarly, the VRM terrain ruggedness ( Figure 5H) further illustrates the topographic complexity of the reticulate morphotype in particular, compared with the annulate and undulate morphotypes. The annulate complexity or ruggedness is pronounced at each individual ring but otherwise is relatively topographically smooth between the rings.

| Subsurface geomorphometry and quantitative morphometrics
Analysis of the seismic isochron thickness shows significant differences in the bioherm thickness between the three morphotypes (Kruskal-Wallis ANOVA, p < .001, Table S6). In summary, the average annulate thickness is 12.69 ± 0.02 m, reticulate is 9.80 ± 0.02 m, and undulate is 6.99 ± 0.03 m (Table 5). These thicknesses were calculated from a timeto-depth conversion based on the average seismic velocity presented in Table 4. The isochron thickness was only calculated where the seismic Reflector A could be confidently delineated, therefore the thickest sediment packages that were too seismically opaque to confidently pick Reflector A are not included in these calculations; these should be considered as minimum thicknesses. By calculating the sum of the bathymetric water depth and bioherm thicknesses the depth to the antecedent substrate was estimated (Table 5). Annulate and reticulate antecedent substrate depths are similar (37.21 ± 0.02 and 38.44 ± 0.02 m.b.s.l. respectively), but undulate depth is ca 6 m less (32.66 ± 0.04 m.b.s.l.) as a result of the shallower average bathymetry and less sediment thickness.
Two sub-bottom profiles from target area 1 (A-A′ westeast and B-B″ north-south) illustrate the typical bioherm topographic relationships with the underlying antecedent substrate ( Figure 6). The undulate morphotype ( Figure 6A) shows parallel to sub-parallel internal bedding overlying a planar substrate reflector. The eastern edge of Inset a grades into the reticulate morphotype ( Figure 6B) showing no apparent change in the planar antecedent substrate. In 2D profile, the reticulate crests and valleys ( Figure 6B) do not correspond with any antecedent features. The intersection of profiles A′ and B cut across the eastern edge of the Halimeda bioherm at a point where the ship track changes direction to north-south. This intersection clearly shows the expression of bioherm surface topography and the prominence of the seafloor and Reflector A at the eastern boundary edge and outer shelf. Reflector A remains planar and the shelf seafloor is covered by a very thin veneer of Holocene sediment. The reticulate bioherm morphotype again commences at the bioherm edge (B-B′, Figure 6C), rising abruptly from 5 to 15 m thick. Figure 6D highlights an annular ring, and the transition from reticulate to annulate also does not correspond with any change in the expression of the antecedent substrate, nor do the crests and hollows of the ring. Profile B-B′ and B′-B″ shows a gradual shallowing upwards of Reflector A from 40 to 30 m.b.s.l. running north-south. However, the corresponding bathymetry of the bioherm tops remains relatively uniform at ca 15 m.b.s.l., due to a thinning sediment package over Reflector A, rather than a 10 m shallowing in bathymetry. The gradual shallowing upwards of Reflector A is a result of the slope geometry of the continental shelf, whereby the top of the profile track at B is proximal to the shelf edge (and hence deeper Reflector A) and B' is more distal from the shelf edge because of the wider shelf here. The antecedent substrate deviates from planar where the outer shelf is incised by palaeochannels and inter-reef passages. In some sections, submerged coral reefs lie immediately proximal to Halimeda bioherms, while in others reef pinnacles protrude through Holocene sediments (Figure 7). Profile A-A′ crosses a broad channel and it is assumed that F I G U R E 6 Examples of surface profile relationships with antecedent basement topography. Reflector A is interpreted as the Holocene/ Pleistocene erosional unconformity (Davies and Marshall, 1985;Orme, 1985;Orme and Salama, 1988). Total sub-bottom profile horizontal distance is approximately 25 km. An acoustic velocity of 1,532 m/s (Table 4) was used to convert from time (TWT ms) to depth domain (m below sea level [m.b.s.l.]) on sub-bottom profiles the overlying sediment package is Halimeda, then continues into a bioherm with undulate to annulate morphology. The channel base highlights the prominent wavy erosional expression of Reflector A, and the overlying sediment is smoothly sinuous with an upwardly convex bedding that grades from parallel with the antecedent substrate, to parallel with the seafloor. The rest of profile A-A′ shows the characteristic annulate topography seen in Figure 4A, but some interesting internal features stand out. Figure 7B and C clearly show evidence of incipient mounds accumulating on Reflector A that are subsequently overlain by younger layers that drape across the low mounds and become parallel to the seafloor. Profile B-B′ crosses a broad channel intersected by a submerged coral reef. The channel has a more complicated antecedent substrate than that in Figure 7A, showing two erosional surfaces. There is evidence of overbank deposits either side of the channel, with laterally continuous younger Holocene sediment deposited on top that is truncated by the emergence of a submerged reef pinnacle ( Figure 7E). The submerged coral reef shows typical steep slopes and 'spikey' topographic F I G U R E 7 Bioherm relationships with adjoining submerged reefs and incised channels showing internal sedimentary architecture and antecedent topography. An acoustic velocity of 1,532 m/s (Table 4) was used to convert from time (TWT ms) to depth domain (m below sea level [m.b.s.l.]) on sub-bottom profiles relief in profile ( Figure 7F). Although somewhat mounded ( Figure 7F), it is unlikely that the overlying thin Holocene veneer of sediment is in-situ bioherm-forming Halimeda, because they favour unconsolidated sediment.

| The Halimeda bioherm carbonate factory
An estimate of the size and volume of the bioherm carbonate factory is obtained from the subsurface calculation of bioherm thickness (Table 6). In total, the bioherms cover an area of 6,111 km 2 with a total volume of 51.5 km 3 . Taking into account average porosity, grain density and % CaCO 3 , the total contribution of CaCO 3 to the system is 19.9 km 3 by volume, and 55.12 Gt by mass (ca 9 M tonnes CaCO 3 km −2 ). The three morphologies differ in their percentage contributions to area and volume respectively ( Table 6). The annulate morphotype contributes 15% of the planar area but 23% of the CaCO 3 volume, due to the greater sediment thickness ( Table 5). The reticulate morphotype contributes 20% of the area and 23% of the CaCO 3 volume, and the undulate morphotype contributes 65% of the area but 54% of the CaCO 3 volume because the sediment thickness is considerably less than that of the other two morphotypes. The values presented in Table 6 are estimates based on the accuracy and potential sources of uncertainty of the input parameters. The thickness (T) is dependent on the accuracy of the seismic velocity time-to-depth conversion, and the grain density, porosity and CaCO 3 % parameters are average values. The bioherm area is dependent on the accuracy of the original delineation of their spatial distribution at the 30 m pixel resolution. Nevertheless, these results represent the best available estimate of the size of the Halimeda bioherm carbonate factory. There are clear differences in area and volume between the three morphotypes which may be reflecting varying aggradation and progradation dynamics.

| DISCUSSION
Prior to 2016, the three different Halimeda bioherm morphotypes were unknown, and morphological descriptions were based on singlebeam echosounder and widely spaced 2D seismic profiles. The dataset is still limited to single lines of 2D sub-bottom profiles, albeit over substantial distances. Crucially, it is now possible to correlate sub-bottom profiles with high-resolution 3D LiDAR bathymetry, allowing the substrate beneath the annulate, reticulate and undulate morphotypes to be explored for the first time. This is a significant advance over previously available data, revealing new insights on the linkages between surface and subsurface geomorphology and similarities and differences in the development of the three morphotypes.
Based on the data available at the time, earlier interpretations of the bioherms described them as a single morphological entity. Some authors recognized a degree of variation, describing that the topographic relief is better developed and more pronounced eastward, with smaller, less complex mounds to westward (Marshall and Davies, 1988;Orme and Salama, 1988). The results presented here provide multiple lines of evidence that the annulate, reticulate and undulate morphotypes are each distinctly different, expressed by their differing geomorphometric characteristics and internal thickness and bedding. These geomorphological differences suggest some degree of heterogeneity between morphotypes with respect to the biotic and/or abiotic mechanisms that modulate their growth and development. First, the differences in physical characteristics between morphotypes based on the BTM morphometric results, and implications for the growth potential of the bioherm carbonate factory are discussed. Second, the subsurface internal architecture and relationships with the pre-Holocene antecedent topography are presented. Third, some biological observations linking living Halimeda meadows to sedimentological accretion are introduced. And fourth, the potential for biotic self-organization in mediating Halimeda bioherm geomorphology are discussed together with some suggested preliminary hypotheses and methodologies for future studies.

| Surface morphometry and quantitative metrics
The three morphotypes show significant differences in sediment accumulation (thickness), and in the depth to the T A B L E 6 Bioherm area and volume, % contribution, and CaCO 3 volume and mass calculations. CaCO 3 volume and mass calculated from measured grain density, porosity and CaCO 3 % average values from Skjold (1988; (Table 5). These metrics equate to vertical accommodation space and growth potential from the pre-Holocene flooding surface. The average depth to the antecedent substrate for annulate and reticulate morphotypes is similar (37.21 ± 0.02 and 38.44 ± 0.02 m.b.s.l. respectively), however, they differ in their average sediment thickness by almost 3 m. Additionally, the undulate morphotype is on average approximately half the thickness of the annulate, but over a much larger area (65% of the distribution). To explain these differences, the three morphotypes must have different growth rates and ability to produce living biomass and therefore sediment accumulation. Or alternatively, if their capacity for growth is the same, then they must be affected by different regimes of sediment redistribution or erosion (deposition). Previous work has established that the three morphotypes have a generalized cross-shelf distribution, with the reticulate morphology generally found closer to the shelf edge than the annulate, and the undulate grading between them but more common to westward (away from the shelf break) (McNeil et al., 2016). This implies that they are probably exposed to differing environmental conditions that control or modulate their growth and morphology. For example, is the undulate morphotype thinner because it grows further away from the shelf edge and possible resource input from shelf-edge upwelling?
The bioherms appear to be constrained to a vertical limit whereby the water depth to bioherm tops is fairly uniform (20-25 m.b.s.l.) for kilometre-scale distances across bioherms. Figure 6D illustrates an example where the antecedent substrate shallows upwards by approximately 10 m due to the gentle gradient of the underlying pre-Holocene continental shelf. But the depth to bioherm tops remains relatively consistent, resulting in a thinning sediment package. On average, the antecedent substrate beneath the undulate morphotype is 4.5 m shallower than the annulate, and the sediment package is on average 6 m thinner. But there is only a 1 m difference in the average bathymetry depth (m.b.s.l.) between undulate and annulate (Table 5). This raises the question of whether bioherm growth is constrained to within ca 20-25 m.b.s.l. Skjold (1988) proposed that the upward limit of bioherm growth in Lloyd Bay (GBR, Figure 3) is controlled by wave-base. The bathymetry data would appear to support this, but it is not possible to determine whether this is a true vertical limit, or whether bioherms still have the capacity to fill vertical accommodation space given more time at present sea level.
Published Halimeda vertical accretion rates from sediment cores in the GBR range from 0.75 to 2.33 m/kyr (Marshall and Davies, 1988;Davies, 2011). These rates imply a potential Holocene thickness of 7.5-23.3 m assuming 10 kyr since initiation (Orme, 1985;Orme and Salama, 1988). Potential vertical accumulation of 23 m is within the range of the present bioherm isochron thickness calculations (Table S5), but also brings the bioherm tops to within the depth of storm wave energy, thus confounding any conclusive statements about vertical growth limits. Laterally, the eastern extent of the bioherm boundary appears constrained by current flow from inter-reef passages at the continental shelf edge ( Figure  5B,F). The western extent may be limited by the availability of nutrients from upwelling (Marshall and Davies, 1988;Wolanski et al., 1988;Drew, 2001;Davies, 2011).
The bioherm thickness (Table 5), area and volume (Table 6) calculations have implications for the size of the bioherm carbonate factory. The size of the carbonate factory and capacity to fill accommodation space varies considerably between the three morphotypes (Table 6). For each morphotype, the relationship between area and volume are not uniform. The annulate morphotype contributes 15% of the planar surface area but 23% of the volume and mass of CaCO 3 ( Table 6). The reticulate morphotype fills the same relative volume (23%) but spread over a larger area (20%). The undulate morphotype covers 65% of the bioherm area, but 54% of the CaCO 3 volume. Although the total contribution from undulate is greatest, the comparative contribution as a proportion of area is less because of the significantly lower average thickness. The annulate morphotype has the greatest capacity to fill vertical accommodation space and highest relative carbonate production (13.58 Mt CaCO 3 km −2 ), while the undulate morphotype has the least (7.48 Mt CaCO 3 km −2 ). From this, it appears that morphotype is related to growth rate, which in turn may be influenced by cross-shelf gradients in resource availability or limitation. Unfortunately, the paucity of Halimeda bioherm sediment cores and published age data prevent further examination of vertical accumulation rate differentiation by morphotype.
A particularly striking difference between morphotypes is the measurably different slopes and terrain complexity (VRM) between annulate and reticulate ( Figure 5C,D,H). The reticulate surface topography appears more complex and less structured compared to the regular spatial organization of the annulate morphotype. An important new finding is the consistent difference between the inward and outward-facing slope gradient (8.9° vs 2.9° respectively, Table S2) of the annulate rings, as this may provide new insights into the processes forming these shapes and patterns. For example, the development of stoss and lee slopes in ripple and dune formation provides evidence of fluid flow direction and velocity (Baas, 1978;Allen, 1984), and steep slopes have been related to biological mechanisms in influencing bedform patterns (Noffke, 1999;Van de Vijsel et al., 2020). As the bioherm annulate rings are circular to sub-circular, it is unclear how these shapes and patterns develop. It is necessary to consider whether hydrodynamic flow and relationships to the antecedent substrate influence the formation and spatial patterns of these structures. The BTM Aspect output illustrates that inter-reef passages and channels influence the overall shape of the bioherm boundaries ( Figure 5F), presumably because flow velocity has a scouring effect that inhibits Halimeda growth in the channels Marshall and Davies, 1988;Drew, 2001). However, there is no evidence of any preferred aspect orientation in the shapes and patterns of annulate rings and reticulate mounds ( Figure 5F). In marine ecosystems, regular bedform patterns mediated by tidal flow can tend towards elongate parallel to flow direction (e.g., estuarine ridges and runnels Williams et al., 2008;Van de Vijsel et al., 2020), or perpendicular to flow direction (e.g., intertidal seagrass; Van der Heide et al., 2010). In the present study, the bioherm tops are fully subtidal, at ca 25 m water depth. In the GBR, flow dynamics at the individual bioherm scale have not been explored, therefore it is not possible to exclude a hydrodynamic (e.g., bottom currents, surface eddies, wave oscillations, etc.) influence on bioherm circular morphology. If hydrodynamics do play a role in shaping the annulate rings, the mechanism(s) would need to be consistent throughout >6,000 km 2 of bioherm distribution across six degrees of latitude.

| Internal architecture and antecedent topography
From the 210 km of sub-bottom profiles available to this study, no evidence was found that the antecedent substrate controls the annulate and reticulate shapes and patterns. Hopley, Smithers and Parnell (2007) describe the Pleistocene surface beneath the bioherms as 'everywhere planar or gently sloping'. Of course, there are some examples where the planar antecedent substrate is incised by palaeochannels or punctuated by submerged reefs (e.g. Figure 7) but bioherms tend not to form on these surfaces, and the observations reported here support this overall planar description beneath the bioherms. The annulate and reticulate crests, and corresponding hollows and valleys that produce up to 25 m of positive relief, do not appear to be inherited from any obvious variation in the antecedent substrate. This is consistent with the interpretation from the much smaller bioherm region in the southern GBR (Searle and Flood, 1988). The submerged coral reefs and pinnacles observed in Figure 7 do appear to initiate off antecedent highs, in contrast with the bioherms, although geomorphological independence from the antecedent substrate has also been demonstrated in some Holocene coral patch reefs (Collins et al., 1996;Blakeway and Hamblin, 2015). The parallel to sub-parallel lamellar bedding and flat surface topography observed in the undulate morphotype does closely follow the planar antecedent topography ( Figure 4B). Some authors (Orme et al., 1978;Davies and Marshall, 1985;Drew and Abel, 1985;Marshall and Davies, 1988) have previously termed these lamellar features as biostromes, sensu Cumings (1932). Davies (2011) discusses the association of Halimeda meadows, biostromes and bioherms and whether these deposits are representing a chronological succession. This discussion is useful as it links the disciplines of geology and biology, and introduces the subject of Halimeda meadows and their potential biological influence on bioherm development. Hillis (1988) describes the development of lagoonal Halimeda meadows at Enewetak Atoll that grow across unconsolidated sediment and produce a vertically 'raised terrace' above neighbouring barren sands. The author proposes a model of Halimeda meadow growth as an important transition towards Halimeda bank (bioherm) development (Hillis, 1988). The growth mechanism at Enewetak is attributed to the extensive Halimeda holdfast network and vertical Unfortunately there was no opportunity to extract samples build-up of shed segments that are subsequently stabilized by cyanobacterial algal mats. New in-situ Halimeda growth accumulates segments which are subsequently stabilized by new algal mats in a repeating cycle (Hillis, 1988). Anecdotal evidence from marine geologists who have scuba-dived on the GBR bioherms in the 1980s do suggest a 'boom and bust' growth cycle. Descriptions range from lush and 'rainforestlike' (Marshall and Davies, 1988) to a 'white desert'. Recent photographs from near Lizard Island show abundant lush green growth ( Figure 8A), however, there is evidence of a recent localized die-off and phase shift to algal mats at this location ( Figure 8C). The photographs of (cyanobacterial?) algal mats ( Figure 8C) were captured opportunistically by a remotely operated vehicle while re-visiting a previous dive site. The specific cause of Halimeda die-off is unknown, but follows 2 consecutive years of catastrophic marine heatwaves in the GBR (2016 and 2017). The apparent phase shift from Halimeda meadow to algal mats may provide some initial evidence of a similar process of cyclical stabilization and aggradation to that described from Enewetak. Recurring phases of cyclical growth/senescence may explain the bedded internal architecture observed in the GBR bioherms, but these are metre-scale thickness seismic features that would be representing centennial to millennial timescales in stratigraphy. A clear mechanism(s) to explain these bedding features remains unresolved, but may involve different biotic and/or abiotic cyclic processes interacting on different timescales (e.g., Williams et al., 2008;Van de Vijsel et al., 2020). A finer-scale seasonal growth/senescence cycle may be overprinted by environmental changes on longer timescales such as climate, cyclones and/or ocean circulation. For example, signals of 'super-cyclones' recurring in 100-300 year intervals since at least the mid-Holocene have been recorded in beach ridge sediments along the GBR, (Nott and Hayne, 2001). To resolve the potential mechanisms and timescales requires higher-resolution stratigraphic analysis combining 3D subsurface geophysical analyses with densely spaced, well-correlated sediment cores and geochemical proxies.

| Potential for biotic self-organization controlling bioherm geomorphology
An emerging area of research in marine carbonate sedimentary systems is autogenic dynamics, or spatial selforganization. Autogenic (self-generating) dynamics and self-organization are increasingly viewed as significant and important sedimentary processes, modulating the formation of biologically constructed features (Budd et al., 2016;Purkis et al., 2016). At the landscape-scale, local interactions between biological and sedimentary processes can induce geomorphological patterns (Weerman et al., 2010), ranging from regular to irregular depending on the feedback mechanism (Pascual and Guichard, 2005;Rietkerk and Van de Koppel, 2008;Weerman et al., 2012) In Modern reef-related systems, metre to kilometre-scale reticulate patterns have been shown to be biologically mediated, rather than substrate controlled (Blakeway and Hamblin, 2015;Schlager and Purkis, 2015;Purkis et al., 2016). Additionally, Purkis et al. (2015) have demonstrated similar morphometric patterns generated by the fossil alga Palaeoaplysina and fossil bryozoan Tubiphytes in the Barents Sea, suggesting that biological self-organization in sedimentary systems is not unique to Modern reef systems, and may indeed be recorded as far back as Precambrian microbialites ( Van de Vijsel et al., 2020).
The spatial patterning in this morphometric and subsurface data, and the potential boom and bust Halimeda/algal mat cycling suggests that the GBR bioherm formation and geomorphology may be influenced by biologically modulated self-organization similar to other marine organo-sedimentary systems. The emergence of large-scale structure or patterning can develop from initially homogeneous conditions Van de Koppel et al., 2008). This could explain the development of complex bioherm topography initiating from a relatively planar antecedent substrate. However, the significant differences between the geomorphology and growth dynamics of the three morphotypes suggests that the underlying mechanisms and feedbacks are probably non-uniform across the bioherm distribution in space and time. For example, the central dieback seen in other clonal plant systems (Cartenì et al., 2012;Borum et al., 2014;Ruiz-Reynés et al., 2017;Vilas et al., 2017a), may initiate ring formation in the annulate morphotype. Then ring accretion results from sediment accumulation, enhanced by Halimeda trapping, baffling and stabilization. Topographic elevation on the ring crests subsequently enhances access to nutrients (a positive feedback), while nutrient limitation in the hollows (negative feedback) would result in trough deepening. Nutrient limitation may explain the thinner, flatter, undulate morphology. Studies of biofilm-sediment systems (Weerman et al., 2010;Van de Vijsel et al., 2020) have demonstrated that topographic flattening can result when the scale-dependent feedback (in this case proposed nutrient limitation) is lost or lacking.
Different feedback mechanisms may generate more irregular spatial patterns (Pascual and Guichard, 2005;Rietkerk and Van de Koppel, 2008;Weerman et al., 2012) that are not necessarily self-organized. Irregular patterns (e.g. the bioherm reticulate morphotype) can be triggered by external disruption of the biophysical feedback Van de Vijsel et al., 2020), signalling critical regime shifts in response to environmental change (Scheffer et al., 2001;Pascual and Guichard, 2005;Weerman et al., 2012). The more irregular topography and abrupt truncations seen in the reticulate bedding could be an indication that such critical shifts have occurred. Given their proximity to the shelf edge and association with upwelling via tidal jets (Wolanski et al., 1988;Drew, 2001), one hypothesis is that the reticulate morphotype may be influenced by higher wave/current exposure, as seen in other spatially organized systems (Noffke, 1999;Guichard et al., 2003;Weerman et al., 2012). Additionally, super-cyclones recurring on centennial timescales (Nott and Hayne, 2001) may potentially trigger cyclical regime shifts in the Halimeda ecosystem.
These preliminary hypotheses are speculative from the presently available data, and require further research and testing. Future studies should aim to test the Halimeda bioherm ecosystem against specific indicators of self-organization which include oscillating consumer-resource interactions, local disturbance-recovery processes , spatial auto-correlation, and pattern persistence over time ( Van de Vijsel et al., 2020). Previous studies on biotic self-organization in marine ecosystems have utilized field-based observations and remote sensing (Temmerman et al., 2007;Weerman et al., 2012;Blakeway and Hamblin, 2015;Purkis et al., 2015;Van de Vijsel et al., 2020), laboratory experiments ( Van de Koppel et al., 2008;Borum et al., 2014), and numerical modelling to test mechanisms against field data (Klausmeier, 1999;Van der Heide et al., 2010;Weerman et al., 2010;Cartenì et al., 2012;Liu et al., 2012;Ruiz-Reynés et al., 2017;Vilas et al., 2017a). To inform such models and identify probable feedback mechanisms and processes in the Halimeda ecosystem requires a better understanding of the local oceanographic regimes at the bioherm scale and of biogeochemical processes such as nutrient cycling and productivity. Field measurements of Halimeda meadow patch size-frequency distributions may inform predictions of external stress or resource limitation (Van der Heide et al., 2010;Weerman et al., 2012), and pattern scaling relationships such as regular, random, clustered, or fractal . Additionally, cross-bioherm investigations of species-specific life history traits (Schwarz et al., 2018) would identify whether different Halimeda species variation drives variable growth rates and therefore influences geomorphology differentiation between morphotypes.
Nevertheless, whatever mechanisms are responsible for the formation of reticulate patterns and annulate rings, in the GBR these processes have commenced very soon after the marine transgression flooded the continental shelf at the end of the last ice age when Holocene Halimeda growth initiated (Orme and Salama, 1988). The proximity of the annulate ring hollows to the antecedent substrate (within 3 m) indicates that vertical aggradation had commenced within centuries (based on estimated vertical accretions rates (Marshall and Davies, 1988) of the marine transgression. Further constraining vertical accretion rates by individual morphotype certainly requires more sedimentological evidence from densely spaced cores and high-resolution radiocarbon dating not presently available.
It is unknown whether the three-dimensional annulate rings are unique to the GBR bioherms, or whether this morphology is common to other Modern Halimeda bioherm locations. Continuous reflection seismic profiles from K-Bank  and Miskito Channel (Hine et al., 1988) show almost identical (to the GBR) topographical relief and internal bedding in 2D, so it is reasonable to assume that they would also show the same 3D expression if the data were available. To recognize the build-up of bioherms requires high-resolution geophysical multibeam, LiDAR and seismic surveys to identify the characteristic surface and subsurface expression of these features. As global seafloor mapping efforts increase, such as through the Seabed 2030 Project (https:// seabe d2030.gebco.net/) it is probable that more Modern and Late Quaternary Halimeda bioherm sites will be discovered. It has been demonstrated that morphometric patterns in Modern carbonate platforms can be applied to the ancient rock record . Therefore, self-organization dynamics in Modern bioherm systems could potentially inform palaeo-environmental interpretations of fossil bioherms and phylloid algal mounds on geological timescales.

| CONCLUSIONS
Halimeda bioherm geomorphology is much more complex and diverse than previously understood, and clearly demonstrates morphotype differentiation. By integrating high-resolution geophysical datasets, combining 3D LiDAR bathymetry and quantitative morphometrics with 2D acoustic sub-bottom profiles, the surface geomorphology and internal sedimentary architecture of three Halimeda bioherm morphotypes have been characterized through space and time. This is a significant advance on previously available data. The main findings are summarized thus: 1. The annulate, reticulate and undulate morphotypes are distinctly different based on significant heterogeneity in their topography, thickness, internal structure, slope gradients and terrain complexity. By extension, their morphology is probably influenced by differing processes of development and biophysical feedback mechanisms. 2. Morphotype differentiation has implications for the nonuniform development of the Halimeda bioherm carbonate factory, rates of sediment aggradation and progradation, and capacity to fill accommodation space. 3. The lateral boundaries between the three morphotypes merge into one another with no apparent geological influence or change in the antecedent substrate, which is typically planar or gently sloping throughout. This is in contrast with adjacent submerged coral reefs that appear to initiate on topographic highs. Morphotype transitions may be reflecting cross-shelf gradients in biophysical growth limiting processes (e.g. access to nutrients and oxygenation), rather than geological mechanisms. 4. The apparent internal bedding seen in sub-bottom profile in all three morphotypes may result from a combination of biotic or abiotic factors operating at different timescales (e.g., cyclical Halimeda growth and death, stabilized by phases of algal mat development, or response to physical disturbance such as cyclones.
To resolve the precise mechanisms and timescales involved requires further analysis, in particular stratigraphic correlations with densely spaced sediment cores. 5. One hypothesis is that biological spatial self-organization similar to reticulate reef patterns and other marine organosedimentary systems modulates, or at least contributes to, Modern Halimeda bioherm geomorphology, in particular the annulate morphotype.
These findings resolve some knowledge gaps, and raise new questions for future work, particularly in developing new testable hypotheses around the mechanisms of potential spatial self-organization. This work contributes a new morphological and morphometric framework to develop field, laboratory and numerical modelling experiments to address these knowledge gaps. Additionally, these results provide geomorphological context to the role of Halimeda bioherms in the provision of complex benthic habitat terrain in the inter-reef seascape of the GBR. These results reinforce the importance of the Halimeda bioherms to the GBR's World Heritage Area status and Outstanding Universal Value from a geological and geomorphological perspective.

ACKNOWLEDGMENTS
Our thanks to Elias Samankassou and Dierk Hebbeln for convening the session titled 'Carbonate mounds through time and space' at the 20th International Sedimentological Congress in Quebec City 2018, out of which this publication arose. We acknowledge and thank Geoscience Australia and the Australian Hydrographic Office for the source bathymetry data used in this study, in particular access to the LADS LiDAR data. This research was supported by a grant of sea time on RV Southern Surveyor from the CSIRO Marine National Facility. Thanks to Emma Kennedy for permission to use photographs in Figures 1  and 8. The GIS analysis was supported by a science grant from The Ian Potter Foundation awarded to LN and MM. Fieldwork in the Great Barrier Reef Marine Park was undertaken under research permit G17-39618.1 and funded by support from the National Geographic Society (EC-190R-18), The National Parks and Wildlife Foundation, and the Great Barrier Reef Marine Park Authority 'Reef Guardians' science grants awarded to MM. The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT
The external datasets that support the findings of this study are detailed in Table S1. The bathymetry DEM and metadata may be accessed from http://pid.geosc ience.gov.au/datas et/ ga/115066. The Topas sub-bottom profile metadata may be accessed from https://www.cmar.csiro.au/data/trawl er/sur-vey_detai ls.cfm?surve y=S00809.