Estuarine biofilm patterns: Modern analogues for Precambrian self‐organization

 This field and laboratory study examines whether regularly patterned biofilms on present‐day intertidal flats are equivalent to microbially induced bedforms found in geological records dating back to the onset of life on Earth. Algal mats of filamentous Vaucheria species, functionally similar to microbial biofilms, cover the topographic highs of regularly spaced ridge–runnel bedforms. As regular patterning is typically associated with self‐organization processes, indicators of self‐organization are tested and found to support this hypothesis. The measurements suggest that biofilm‐induced sediment trapping and biostabilization enhance bedform relief, strength and multi‐year persistence. This demonstrates the importance of primitive organisms for sedimentary landscape development. Algal‐covered ridges consist of wavy‐crinkly laminated sedimentary deposits that resemble the layered structure of fossil stromatolites and microbially induced sedimentary structures. In addition to layering, both the morphological pattern and the suggested formation mechanism of the recent bedforms are strikingly similar to microbialite strata found in rock records from the Precambrian onwards. This implies that self‐organization was an important morphological process in times when biofilms were the predominant sessile ecosystem. These findings furthermore emphasize that self‐organization dynamics, such as critical transitions invoking ecosystem emergence or collapse, might have been captured in fossil microbialites, influencing their laminae. This notion may be important for paleoenvironmental reconstructions based on such strata. © 2019 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd

Earth's geological record boasts ample evidence for longlived sediment biostabilization. Microbialites (Burne and Moore, 1987) represent a substantial part of the sedimentary rock record from the Archean onwards (e.g. Bosak et al., 2013;Noffke and Awramik, 2013). These deposits are formed by microbial communities, thin biofilms that can develop further into cohesive microbial mats, which interact with water flow and sediment transport by excreting glue-like extracellular polymeric substances (EPS) or forming a dense mesh of filaments (Stal and Caumette, 1994;Stal, 2010). In siliciclastic depositional environments, these biophysical interactions can result in often laminated, microbially induced sedimentary structures (e.g. Noffke, 1999;Noffke et al., 2001). In carbonate-rich environments, microbial mats may also precipitate minerals, forming microbialites such as the pronouncedly laminated stromatolites, with layers typically micrometres to millimetres thick (e.g. Riding, 2000Riding, , 2011Bosak et al., 2013;Noffke and Awramik, 2013). There is mounting evidence for the important contribution of microbial mats on sediment stability (e.g. Microcoleus chtonoplastes literally means 'soil former') and hence the formation of such structures (e.g. Noffke, 1999;Schieber, 1999;Noffke et al., 2007), although lamination, a key aspect of many microbialites, might also arise due to abiotic sedimentary processes (e.g. Grotzinger and Rothman, 1996;Grotzinger and Knoll, 1999;Brasier, 2011). Microbialites could grow for up to thousands of years, range in spatial scale from millimetres to kilometres and are often clearly spatially patterned, with pattern wavelengths ranging from the order of 1 cm to 10 m (e.g. Bosak et al., 2013;Noffke and Awramik, 2013;Mariotti et al., 2014). This patterning is often regular, e.g. for the Holocene stromatolites of Shark Bay (Australia) and the Bahamas (e.g. Logan, 1961;Andres and Reid, 2006;Bosak et al., 2013), but not always (e.g. Noffke, 1999Noffke, , 2008. Fossil microbialites, given their well-preserved layer sequences and abundance in the geological record, are valuable for paleoenvironmental reconstructions (Logan et al., 1964;Jerzykiewicz and Sweet, 1988;Harzhauser et al., 2014;Lettéron et al., 2017). Despite the valuable information about formation mechanisms and environmental conditions that can be deduced from spatial patterns (e.g. Pascual and Guichard, 2005;Rietkerk and van de Koppel, 2008), microbialite patterns and their implications for paleo-reconstructions have received relatively little scientific interest.
Present-day biofilms and their bedforms may serve as analogues to Precambrian and Phanerozoic microbialites (e.g. Noffke, 2008;. Although stromatolites declined rapidly after the Precambrian, likely due to the evolution of burrowing and grazing animals (e.g. Walter and Heys, 1985) and changes in sea water chemistry (e.g. Stal, 2010), microbially induced sedimentary structures did not experience such a drastic drop in numbers and sediment-sculpting biofilms are still ubiquitous on present-day tidal flats (e.g. Noffke, 2008). These communities of primitive unicellular organisms can potentially cover and stabilize large intertidal areas rapidly, typically within one growth season (e.g. Noffke and Krumbein, 1999;de Brouwer et al., 2000). Just like their ancestors, these modern biostabilizers commonly exhibit spatial patterning (e.g. Noffke, 1999;Blanchard et al., 2000;de Brouwer et al., 2000;Lanuru et al., 2007;Weerman et al., 2010Weerman et al., , 2012. Regular spatial patterns can be induced by scale-dependent feedbacks, whose sign and magnitude depend on distance (Rietkerk and van de Koppel, 2008). Regular patterns in biofilms (Weerman et al., 2010) and salt marsh vegetation (Temmerman et al., 2007) are thought to result from scaledependent feedbacks in which water flow diverges away from topographic elevations. These elevations are then preferentially colonized and biostabilized, creating a local positive feedback. Flow convergence in surrounding topographic lows leads to increased inundation and shear stress and hence reduced biostabilization. Thus, this feedback is 'scale-dependent', as it is locally positive but has 'reversed sign' more distantly.
Although fairly similar morphologies are found in microbial mat-covered erosional remnants and uncovered erosional pockets, such bedforms are shaped by a runaway erosional disturbance, causing progressive pocket erosion in a formerly intact mat (Noffke, 1999(Noffke, , 2010Noffke and Krumbein, 1999). Although this is a positive feedback process as well, it is markedly different from the previously explained scale-dependent feedback. This might explain why erosional remnants and pockets appear not to be clearly regularly patterned. In contrast, it is likely that regularly spaced intertidal ridge-runnel bedforms do arise from a scale-dependent feedback (although this term is not used in the ridge-runnel literature): spontaneous flow instabilities are amplified by flow convergence above and net erosion of poorly drained runnel sediments and net sedimentation related to flow divergence and higher sediment stability above well-drained ridges (e.g. Allen, 1987;Gouleau et al., 2000;O'Brien et al., 2000;Williams et al., 2008;Carling et al., 2009). Although this feedback can be purely geophysical, ridge colonization by biostabilizers might amplify the stability contrast (ridge vs runnel) and hence the feedback (e.g. Blanchard et al., 2000;Lanuru et al., 2007;Williams et al., 2008;Weerman et al., 2010;Fagherazzi et al., 2013).
The emergence of large-scale structure or patterning from initial disorder due to local (scale-dependent) interactions is known as spatial self-organization (Rietkerk and van de Koppel, 2008); a Turing instability (Turing, 1952) is one of the mechanisms behind pattern emergence. Self-organization dynamics are typically nonlinear, possibly invoking bistability (e.g. degraded vs restored ecosystem state), with drastic, critical transitions between states (e.g. May, 1977;van de Koppel et al., 2001van de Koppel et al., , 2005Rietkerk et al., 2002;Weerman et al., 2010;Scheffer et al., 2012;Liu et al., 2014b). Such critical transitions might also be more stepwise in multi-stable systems, allowing a range of different pattern wavelengths (Bastiaansen et al., 2018). Thus, it is essential, for example for paleo-reconstructions based on microbialites, to understand how their patterns are formed, as self-organized regularly patterned systems may behave very differently from disturbance-driven irregularly patterned systems (e.g. Noffke, 1999;Guichard et al., 2003;Pascual and Guichard, 2005;Weerman et al., 2012).
Present-day self-organized biofilm bedforms (Weerman et al., 2010) show strong resemblance, both in morphology and underlying processes, to regularly patterned, dome-shaped stromatolites (e.g. Logan, 1961). Yet, although self-organization in stromatolites has been suggested before (e.g. Petroff et al., 2010;Brasier, 2011;Tice et al., 2011), the possible implications thereof have not been broadly explored yet (Purkis et al., 2016). Possibly, the hypothesis of self-organization has not received much attention because an essential and defining aspect of many microbialiteslaminationis not always found or evident in modern mudflat bedforms (e.g. Weerman et al., 2010). Lamination, principally a consequence of abiotic sedimentary processes (e.g. Gouleau et al., 2000;Williams et al., 2008;Carling et al., 2009), can be enhanced by long-term biostabilization with temporally variable binding rates (Williams et al., 2008). These biogenic effects can be preserved as wavy-crinkly layering (Horodyski, 1982;Schieber, 1999;Riding, 2011;Noffke et al., 2007;Sarkar et al., 2016). Yet, instead of long-term biostabilization, many previous studies have focused on environments with strong seasonality, thus emphasizing the short-term aspects of modern biofilms Le Hir et al., 2007;Weerman et al., 2010Weerman et al., , 2011Weerman et al., , 2012. As microbial biofilms form a protective 'skin layer' on the sediment, the underlying sediment regains its higher, abiotic erodibility as soon as the biofilm is eroded by grazing, hydrodynamics or both van de Koppel et al., 2001;Le Hir et al., 2007;Weerman et al., 2010Weerman et al., , 2011Weerman et al., , 2012, such that sedimentary deposits characterized by well-defined laminae cannot develop. Sediment lamination has been reported on mudflats where ridge and runnel bedforms persist multiannually, but without consideration of the influence of biofilms or similarities with stromatolites (Gouleau et al., 2000). To answer the question of whether the formation and regular patterning of many layered microbialites could be self-organized, modern analogues are needed which are morphologically and mechanistically similar, but which also withstand biophysical seasonality enough to construct long-lasting, laminated deposits.
Mats composed of simple multicellular biostabilizers such as macroalgae might have an even stronger potential to form such layered deposits. Primitive filamentous macroalgae such as Enteromorpha sp. (e.g. Raffaelli et al., 1999;Romano et al., 2003;Le Hir et al., 2007) and Vaucheria sp. (e.g. Simons, 1975;Gallagher and Humm, 1981;Wilcox, 2012) have been reported to persist multiannually (Romano et al., 2003). Moreover, all biophysical interactions characterizing microbially induced sedimentary structures (Noffke et al., 2001) have also been described for Vaucheria mats. Vaucheria baffles, traps and binds sediment (Black, 1933;Skowroński et al., 1998) in its dense mesh of filaments that can penetrate the soil down to 4-6 cm and grow upwards when buried (Gallagher and Humm, 1981), effectively biostabilizing (Webber, 1967;Paterson, 1994) and forming thick and patchy organosedimentary mats (Simons, 1975;Gallagher and Humm, 1981;Wilcox, 2012). However, the impact of such mats on long-term morphodynamics, sediment lamination or bedform patterning has not been explored in detail. Although these mats are not biofilms in their strictest definition (i.e. they are not microbial and do not excrete large amounts of EPS; e.g. Stal, 2010), their primitive form, rapid and sheet-like growth, and distribution along the entire salinity range (Simons, 1975) makes them similar to microbial biofilms. Multicellular algae like Vaucheria were rare in the Precambrian, but not absent (Butterfield, 2004). For the purpose of this study, the definition of biofilms will therefore be slightly broadened (following Fagherazzi et al., 2013) to include functionally similar mats composed of multicellular macroalgae. Do such mats create laminated and patterned bedforms, contributing to long-term estuarine morphodynamics just like cyanobacteria and other microbialite-building organisms did from the Precambrian onwards (Noffke, 1999(Noffke, , 2010Noffke and Krumbein, 1999;Riding, 2011)? This study reports on present-day intertidal ridge-runnel bedforms that seem regularly spaced and persistent for multiple years, with long-lived Vaucheria mats covering topographic highs ( Figure 1). These algal mats will also be referred to as 'biofilms' throughout this study. These metre-scale bedforms are observed on intertidal flats in the Schelde estuary in Belgium, close to the Dutch border (see Figures SI1 and SI5). Vaucheria bedforms are sometimes more dome-shaped (Figure 2), broadly resembling the morphology of many stromatolites (Logan, 1961;Andres and Reid, 2006). These observations raise multiple questions. Are the Vaucheria bedforms indeed regularly patterned and could self-organization explain this? How long do the algal and topographic patterns persist? Are the bedforms laminated, resembling layered microbialites? To address these questions, specific indicators of self-organization due to biophysical scale-dependent feedbacks will be tested, including bedform regularity and orientation, pattern formation history and local biostabilization by algae. Moreover, pattern persistence will be quantified using remote sensing data. Finally, the bedforms will be examined for wavy-crinkly lamination, a key property of many microbialites.

Study area
Field measurements were done on the intertidal flat of Ketenisse (51.285063°N, 4.312082°E), in the far upstream mesohaline zone of the Schelde estuary, with water salinity fluctuating around 3.6 ppt in 2006 (van den Neucker et al., 2007). This coastal realignment site, consisting of formerly embanked agricultural land, was de-embanked in 2002 to restore it as a mudflat and has been accreting estuarine sediments ever since. The mudflat is currently embanked on the south and west side, with salt marsh vegetation (e.g. Bolboschoenus maritimus, Phragmites australis) in front, and is bordered by the estuarine channel along the north-to-east border ( Figure SI1). Early 2003, MLW, MHW and MHWS were at respectively À2.3, +3.0 and +3.4 m relative to mean sea level (MSL); at that time the restored mudflat was elevated from about 1.9 to 2.3 m above MSL (van den Neucker et al., 2007). From visual inspection in the field between 2015 and 2017, it seems that these ridges and runnels are elongated along the direction of downhill drainage flow. It also appears that some of the runnels are somewhat deeper and streamwise more continued (i.e. more similar to drainage creeks bordered by creek levees instead of runnels bordered by ridges). This is also observed elsewhere (Gouleau et al., 2000;O'Brien et al., 2000;Whitehouse et al., 2000). The ridges and creek levees are covered by algae, primarily Vaucheria, while runnels and creeks are barely covered. Contrary to diatom-stabilized bedforms, which tend to be spatially regular but are often seasonally ephemeral Weerman et al., 2010Weerman et al., , 2012, or to more persistent but less regular erosional remnants and pockets (Noffke and Krumbein, 1999), it appears that the Vaucheria bedforms are spatially regular ( Figure 1) as well as persistent for many consecutive years, taking aside periods of extreme cold, heat or drought ( Figure SI9).

Remote sensing techniques
In order to quantify the spatial and temporal characteristics of algal pattern and bedforms on the Ketenisse mudflat, an area of approximately 50 × 50 m (around 51.285063°N, 4.312082°E; Figure SI1) was marked with six small poles in August 2015. The poles were firmly anchored in the sediment and georeferenced with an RTK-dGPS (Leica Viva GNSS GS12 receiver with CS15 controller). This delimited area was left untouched and was analysed repeatedly over the course of time, using several remote sensing techniques.
Bedform geometry and temporal dynamics were measured with a tripod-mounted 3D terrestrial laser scanner (RIEGL VZ400). By emitting a laser beam, the rotating scanner creates a three-dimensional point cloud of the sediment surface. To georeference this point cloud, reflectors (5 cm diameter, RIEGL) were mounted on the six poles surrounding the study area. The 50 × 50 m area was scanned from at least three different vantage points to avoid shadowing effects due to surface topography. Each scan position was chosen such that at least four reflectors were captured; at least three overlapping reflectors were used to merge separate scans. Raw scanner data was processed using RiScan Pro (RIEGL) to obtain georeferenced x, y, z-coordinates (x, y horizontal; z vertical). All subsequent post-processing steps were achieved using Matlab (MathWorks, Inc.). Using linear interpolation, this x, y, z-point cloud was gridded onto a regular horizontal grid with 10 cm grid spacing to obtain a digital elevation model (DEM). To accelerate the computation time while retaining precision, interpolation was conducted with random subsamples of, on average, 50 x, y, z-points per grid cell. Bed-level elevation is expressed in metres above NAP (Amsterdam Ordnance Datum, roughly coinciding with mean sea level in Amsterdam, The Netherlands). To obtain relative bed-level elevation, the background elevation profile S N (x, y) of the form with polynomial order N and constants a was calculated by regression and subtracted from the total elevation profile. Here, N = 4 was chosen as this gave higher R 2 than N = 2 or 3. As the delimited 50 × 50 m area was not exactly rectangular, a 44 × 49.5 m rectangular area was selected within the 50 × 50 m area (44 m south-to-north, 49.5 m west-to-east). The northward direction can be considered 'downslope' within this area, as the areamean background surface gradient ∇S N was oriented 359.05°f rom the north in August 2015. For a comparison of bedform to algal pattern (discussed below), three DEMs (11 August 2015, 13 December 2016 and 19 September 2017) were used. Aerial photos were taken, to analyse algal mat coverage and correlate this with bed level. To make a year-on-year comparison, aerial photos acquired in 2015 (11 August), 2016 (13 December) and 2017 (19 September) are considered. The photos were captured with a DJI Inspire 1 drone with Zenmuse X3 gimbal and camera, with sufficient overlap to make an orthophoto. The photos were stitched together using Agisoft PhotoScan Pro to obtain an orthomosaic with a resolution of 1.75 cm (2015), 0.98 cm (2016) or 1.3 cm (2017). Four white reference plates ( Figure SI9b) were placed in the field adjacent to the corner poles and georeferenced with the RTK-dGPS. In Matlab, either these four reference plates or the outermost four reference poles were used as control points to co-register the non-referenced orthomosaic to the georeferenced DEM. Each orthophoto mosaic was geometrically transformed and aligned with the DEM (i.e. with 10 cm grid spacing) to allow for a grid cell-wise comparison. The algal pattern within the 44 × 49.5 m focus area was quantified using the visible-band difference vegetation index (VDVI) (Du and Noguchi, 2017), based on the image's green (G), blue (B) and red (R) pixel values: The relative VDVI profile was computed by regression of the background profile S 4 (x, y) and subtraction from the total VDVI profile; this was done to correct for spatially large-scale spectral imbalances in the orthomosaic, hence to retain the smallerscale algal pattern. Algal presence farther back in time than the available laser scan and drone data was quantified from aerial images of the 44 × 49.5 m area obtained with Google Earth Pro (© 2016 Google Inc.). These images were acquired for 2004 (8 June), 2009 (31 August), 2012 (4 September), 2013 (7 July), 2015 (1 October) and 2017 (24 September). The original photos had a resolution of approximately 1.8 cm and were also geometrically transformed and aligned with a rectangular 10 cmresolution grid. VDVI values were computed and a fourth-order background profile was regressed and subtracted to obtain relative VDVI values.
To test if the spatial algal pattern overlaps with the topographic pattern and whether this correlation persists in time (possibly contributing to lamination), elementwise linear correlation (Pearson's r) was performed, using Matlab, for the DEMs and drone-derived VDVI maps acquired in 2015, 2016 and 2017. Since the residuals were not homogeneously distributed, p-values could not be computed directly from the linear correlation itself. Instead, p-values were computed by randomly permuting the elements of each DEM or VDVI matrix 1000 times, recomputing the correlation coefficient at each permutation and calculating the probability of obtaining an equal or more extreme r-value than the 'true' r by chance.
Two-dimensional auto-correlation and cross-correlation matrices of the DEMs and VDVI maps of 2015, 2016 and 2017 were computed to quantify pattern regularity (van de Koppel et al., 2005;Weerman et al., 2010Weerman et al., , 2012Cornacchia et al., 2018) and persistence over time. To quantify spatial regularity (auto-correlation) and displacement (cross-correlation) in west-to-east orientation (resp. north-to-south orientation), the row vector running eastward from each matrix's centre (resp. the column vector running southward from each matrix's centre) was analysed. Wavelength is quantified as the correlation distance where spatial auto-correlation is maximal, after increasing from negative to positive values. To assess the statistical significance of the auto-and cross-correlations, random permutation (as described above) was applied; each DEM and VDVI matrix was permuted 100 times. The p-value was given by the probability of finding an equally or more extreme correlation for that specific correlation distance by random chance. Since the DEM and VDVI matrices had very different units, the correlation matrices were normalized, for plotting purposes, by division by their maximal value. The same procedure was applied to VDVI matrices derived from Google Earth images ( Supplementary Information, Figure SI7).
To estimate how many of the runnels are slightly deeper (i.e. more similar to creeks), a histogram (1000 bins) was made of relative bed level in 2015. This frequency distribution is bimodal (representing ridges and runnels), with a longer tail (representing creeks) developing towards lower elevations, distinguished from the bimodal distribution by a 'bump'. This classification allows creek/runnel/ridge distinction in the DEM. From the same histogram (and similarly for 2016 and 2017), modal (i.e. most frequent) ridge-to-runnel elevation differences (a measure for bedform topographic amplitude) were computed as the difference between most frequent ridge elevation and most frequent runnel elevation.

Sediment lamination
Four sediment cores, each approximately 30-35 cm long and 7 cm in diameter, were taken on the Ketenisse mudflat (September 2016) to test the hypothesis that biofilm-covered bedforms consist of laminated deposits, while bare tidal flat does not. Two cores (referred to as 'SE'; 51.284800°N, 4.312494°E) were taken about 11 m southeast of the southeastern corner of the 44 × 49.5 m area. The other two (referred to as 'NE'; 51.285267°N, 4.312572°E), meanwhile, were cored about 13 m east of the northeastern corner of the study area. At both locations, one core was taken on an algal-covered creek levee, the other one in the bare runnel behind (as seen from the creek) the levee. In the lab, the cores were sliced in half lengthwise and an X-radiograph was taken (HP Faxitron Cabinet X-ray model 43805 N). The cores were radiated for 2.5 min with a 50 kV source at 61 cm distance between source and lightsensitive film (Konica Minolta Regius RC-300 cassette). These X-radiographs (resolution 0.085 mm) were digitalized and converted to grey-value matrices (pixel intensity ranging from 0 to 255) and a subset (25 cm from the surface downwards, 6 cm wide) was extracted to exclude boundary effects near the sides and bottom of each core. Relative X-radiograph profiles were computed by regression and subtraction of background profiles S 4 (x, y). In this way, larger-scale light intensity changes (e.g. due to the semicircular shape of the sliced core) were filtered out.

Sediment biostabilization
To test the assumption of sediment strengthening by Vaucheria, penetration resistance, a measure of bed compaction (Williams et al., 2008) or algal mat binding, was measured for all the sediment cores taken on the Ketenisse mudflat. Nine sediment cores (PVC tube, diameter 12 cm, length 20 cm) were taken in June 2015; three cores at three sites. Per site, one core was taken on an algal-covered creek levee, one core nearby where the same creek levee was not visibly algal-covered and one core in the bare runnel behind it (i.e. away from the creek). The three sites were approximately spaced at 10 m intervals along the same creek, with the central coring site at 51.285131°N, 4.311331°E ( Figure SI1). In the lab, penetration resistance was measured with an INSTRON Testing System controlled with INSTRON Bluehill 3 software. Due to technical issues, this lab analysis was delayed for 2 weeks. In the meantime, and so as to preserve their integrity and physical properties, the cores were kept in controlled conditions ( Supplementary Information, SI6). As this analysis was done to study relative differences between the substrates, rather than absolute values, the influence of this delay is not considered important. Penetration resistance as a function of depth (upper 100 mm, increments of 0.17 mm) was measured, in the centre of the sediment cores, with a metal cone (frontal area 2 cm 2 , tip angle 60°). Median penetration resistance was computed for each profile; further calculations were done with these medians. After testing for normality (Shapiro-Wilk normality test) and homoscedasticity (Levene's test) of the nine residuals, a two-way analysis of variance (ANOVA) was performed of the fitted linear model, with both substrate (runnel, bare levee, algal-covered levee) and location (sites 1-3) as fixed factors. Analyses and statistics were performed using R (http://www.Rproject.org).

Self-organization of algal pattern and bedforms
A regularly spaced and dashed algal pattern emerged seemingly spontaneously, as can be seen on aerial photos of the 44 × 49.5 m focus area on the Ketenisse mudflat (Figure 3). Although Vaucheria had already colonized this mudflat by 2003 (van den Neucker et al., 2007) (i.e. shortly after this area was de-embanked to become intertidal), the first Google Earth image from 2004 shows no signs of colonization in the focus area yet. Aerial images do not exhibit regular patterning before 2012. Some irregular clustering in a west-to-east and northwest-to-southeast orientation is observed in 2004, however. In the next aerial image available in the time series (2009), the entire focus area is covered by algal patches, but they show no apparent regularity. By contrast, in 2012, 2013, 2015 and 2017, a regularly spaced algal pattern was present, with dominant wavelengths of 2.4, 2.6, 2.5 and 3.0 m, respectively, in an east-to-west direction, as quantified by the spatial auto-correlation (p < 0.01; Figure SI7a). No clear wavelength was detected in the north-to-south direction for 2012-2017. Background bed level S 4 (x, y) in 2015 sloped down towards the north. The algal pattern visible between 2012 and 2017 is therefore oriented streamwise. Thus, emergence of a flowaligned regular algal pattern, apparently without pre-existing topographic structuring, is consistent with the hypothesis of self-organization due to scale-dependent feedbacks, with or without active biotic involvement.

Vaucheria bedform morphology
The algal pattern is in good accordance with the underlying topographic pattern (i.e. the ridges have algal cover and the runnels are bare). Apart from this being visible in the field (Figures 1 and 2), the association follows from elementwise linear correlation of laser scan-derived relative bed level with drone-derived relative VDVI within the focus area (Pearson's r = 0.488 for 11 August 2015, 0.552 for 13 December 2016, 0.384 for 19 September 2017; p < 0.001 for all three dates).
Spatial auto-correlations (p < 0.01) of both relative bed level and relative VDVI (Figure 5a) confirm that bedform ridges are elongated along the direction of downstream drainage flow and are regularly spaced in the cross-flow (spanwise) direction. The dominant wavelength (i.e. distance of maximal correlation) in the cross-flow direction is 2.6, 2.6, 2.8 m for DEM 2015, 2016, 2017 respectively and 2.3, 2.4 m for VDVI 2015, 2016, respectively. When the algal pattern gradually starts broadening and fragmenting towards 2017 (Figures 3 and SI9a), a dominant wavelength can no longer be computed (i.e. the autocorrelation of VDVI 2017 has a maximum at 2.5 m but does not change sign). No clear wavelength (at least none smaller than 10 m) exists in the northward (downslope) direction for either the DEMs or drone-derived VDVI matrices from 2015, 2016 and 2017. Roughly every 6-8 m in the spanwise direction (in 2015), ridges are longer (more continuous) in the streamwise direction (i.e. more similar to creek levees) and the runnels they border are slightly deeper (i.e. more similar to drainage creeks). The modal ridge-runnel elevation difference is 5.3 cm in 2015, 5.2 cm in 2016 and 2.3 cm in 2017. In winter and during hot and dry summer periods, algal cover is substantially reduced, or sometimes not even visible. Between May and July 2017, it was visually observed that algal cover gradually shifted from the ridges to the (sides of the) runnels ( Figure SI9a). By September 2017, however, the algae had started to recolonize the ridges. Further testing of the algal growth potential in bare sediments (van Belzen, 2018) revealed that Vaucheria spores were omnipresent and could grow everywhere in the field, including bare areas, given the right growth conditions ( Supplementary Information, Figure SI3). Thus, the observations of regularly spaced ridge-runnel bedforms in combination with stronger Vaucheria growth potential on ridges are in agreement with the hypothesized biophysical scale-dependent feedback.

Multi-year pattern persistence
For several consecutive years, the regularly spaced algal mat pattern overlaps with the bedform pattern and remains at the same location on the Ketenisse mudflat. Between 2015 and  Figure SI9a) are still measurable in the 2017 VDVI data, such that this year correlates less strongly with 2015 and 2016. The persistence of algal pattern and bedform wavelength can also be seen from the auto-and cross-correlation profiles ( Figure 5). In the cross-flow direction, cross-correlations follow a similar profile as auto-correlations, emphasizing that between subsequent years, the regular topographic and algal patterns mostly stay in the same location. The 0.3 m shift in VDVI cross-correlation between 2016 and 2017 reflects the after-effect of the temporal algal relocation ( Figure SI9a). Between 2015 and 2017, the entire 44 × 49.5 m focus area gradually accretes sediment over the course of time and rises about 2.9 cm year À1 along its northern edge ( Figure SI2) and about 4.0 cm year À1 along its southern edge ( Figure 4). Bedforms are best preserved on the southern edge of the focus area (as also observed for the algal pattern, Figure 3) and seem to gradually flatten out on the northern edge. Long-term persistence is also evident from the autoand cross-correlations computed for the Google Earth time series (Figures 3 and SI7): the algal pattern does not change much in wavelength between years and stays mostly in the same location. Hence, overall, the regularly spaced algal pattern and underlying bedforms are persistently present for multiple consecutive years, despite significant accretion rates.

Linking sediment lamination and biostabilization
Algal-covered creek levees are horizontally wavy-crinkly layered and have more stable sediments than bare runnels, which do not exhibit such lamination. Layering is clearly visible in the X-radiograph of the algal core taken around the southeastern (SE) edge of the focus area ( Figure 6). Broad dark and light bands alternate with about 4.5 cm vertical spacing, as shown by the spatial auto-correlations ( Figure SI4). Finer millimetrescale wavy-crinkly laminae are also visible, mostly within the broad dark bands of this core ( Figure SI8). Although some of these undulations appear to coincide with burrows (a, b in Figure SI8), crinkly textures are also observed away from bioturbation tracks (c, d). The ca. 1 cm wavelengths of the crinkles are of roughly the same scale as the tufted, wavy Vaucheria mat texture (Figure 1). Assuming 4 cm year À1 accretion, the broad dark bands with finer crinkles, at between 12 and 20 cm depth, date from between 2011 and 2013, coinciding with the Google Earth images that exhibit most pronounced algal cover (Figure 3). In the northeastern (NE) algal core, broad (4.6 cm spacing, Figure SI4) and fine (millimetre-scale) lamination are vaguely visible, yet less evident than in the SE-algal core. Broad-and fine-scale lamination are even more vague in the sediment cores taken in bare runnels (both on the NE and SE side). Penetration resistance (Figure 7) of Vaucheriacovered creek levees (194 kPa, i.e. least square mean of the three profile medians) was significantly higher than that of bare creek levees (89.3 kPa) and bare runnels (13.7 kPa), as follows from pairwise comparisons of the least square means (all p < 0.0001, adjusted with the Tukey method). All substrates have significant effect (two-way ANOVA, p < 0.001), whereas sites do not. Thus, the data show concurrence of algal cover, stable sediments and wavy-crinkly laminae.

Discussion
Present-day intertidal flats are often typified by regularly patterned bedforms shaped by biofilms, highlighting that spatial self-organization is an important intertidal landscape-building process. Strikingly, similarly regularly spaced biofilm bedforms are abundantly found in the fossil record as laminated microbialites, dating back to the Precambrian. Yet, modern biofilm patterns are often ephemeral, which inhibits the build-up of sedimentary lamination. This leaves the question of to what extent self-organization has shaped shallow-marine landscapes throughout Earth's geological history. This study highlights that biofilm-forming, primitive biostabilizers can stimulate selforganization of modern mudflats into regularly patterned, multi-annually persistent and hence internally laminated ridge and runnel bedforms. Flow deceleration and trapping, binding and biostabilization of ridge sediments by mats of the filamentous macroalgae Vaucheria sp. (functionally similar to microbial biofilms) likely boosts the formation and multi-year persistence of metre-scale ridge-runnel topography and enhances layered deposition, forming typical wavy-crinkly laminae. Both the crinkled layers and regular bedform morphologies are strikingly similar to modern and ancient microbialites, ranging in morphology from planar microbially induced sedimentary structures to dome-shaped stromatolites. This study therefore suggests that spatial self-organization, due to biophysical scale-dependent feedbacks, might be an important process shaping inter-and subtidal ecosystems, not only in the present, but dating back to the Precambrian, when biofilms were the predominant benthic ecosystem on the planet.

Algal/topographic pattern regularity and indications of self-organization
Multiple measurements indicate that the observed algal bedforms are likely self-organized, as explained below. Regularly spaced algal-covered ridges and bare runnels are observed from 2012 onwards (Figures 3, 4, 5 and SI7). Yet, such structuring appears absent prior to 2012. The irregular clustering observed in 2004 (Figure 3), with different orientation than the pattern from 2012 onwards, probably reflects drainage topography, given the underlying terrain slope in 2004 (see figure 5.14 in van den Neucker et al., 2007). In 2009, algal cover (and hence also topography, given the positive correlation found between cover and elevation) was irregular ( Figure 3). Hence, spatial regularity arises from seemingly disordered initial conditions, which is consistent with self-organization principles (e.g. Rietkerk and van de Koppel, 2008). This suggests that internal feedbacks are important structuring processes. This study indicates how such a scale-dependent feedback loop might work. Algal mats are preferentially found on elevated ridges (Figure 1), despite algal spore presence across the tidal flat ( Figure SI3). Ridge sediments, in turn, are more stable (Figure 7), most likely due to biostabilization (Webber, 1967;Paterson, 1994), combined with better drainage and lower shear stress compared to runnels (e.g. Gouleau et al., 2000;O'Brien et al., 2000;Williams et al., 2008;Carling et al., 2009). This mechanism, as well as the observed flow-alignment of the Vaucheria bedforms ( Figures 5 and SI1) is consistent with the biophysical scale-dependent feedback described by Temmerman et al. (2007) and Weerman et al. (2010). The idea of a feedback loop is further supported by the fact that algal bedforms retain their wavelength and amplitude for several years, despite rapid sediment accumulation (Figures 3-5 and SI7): without site-dependent sediment erosion and stabilization, runnels (or creeks) would be expected to fill up and/or ridges (or creek levees) to flatten out (e.g. Williams et al., 2008). Indeed, the observed bedforms are significantly flattened out in autumn 2017, following a hot summer in which it was visually observed that algal cover temporarily shifted from the ridges to (the sides of) the runnels ( Figure SI9a), which is likely the result of heat or desiccation stress. It is tempting to assume a causal relation, yet a multi-annual time series (encompassing both topography and algal distribution) with higher temporal resolution is required to more firmly support the notion that biofilm cover directly affects bedform geometry. In-depth study of the seasonal changes in algal cover and bedform morphology, as well as zonation along the tidal range, could for example be quantified with the biotic modification index introduced by Noffke and Krumbein (1999). While purely geophysical mechanisms without biological involvement can also explain the formation of ridge-runnel bedforms (Whitehouse and Mitchener, 1998;Blanchard et al., 2000;O'Brien et al., 2000;Lanuru et al., 2007;Williams et al., 2008;Carling et al., 2009), evidence for the importance of stabilizing biofilms on the formation and preservation of sedimentary bedforms has accumulated in the past years (e.g. de Boer, 1981;Blanchard et al., 2000;Friend et al., 2008;Noffke, 2010;Weerman et al., 2011;Mariotti et al., 2014). Although the abiotic mechanisms for ridge-runnel formation proposed in previous studies are also essentially scale-dependent feedbacks, the possibility and potential consequences of ridge-runnel self-organization (Yang et al., 2012) have received considerably less scientific attention in the geophysical context than self-organization typically enjoys in the ecosystem sciences. This study shows that regularly patterned algal bedforms are likely formed by self-organization due to scale-dependent feedbacks. This implies that tidal flats can potentially undergo catastrophic shifts between alternative states (e.g. May, 1977;van de Koppel et al., 2001;Scheffer et al., 2012;Weerman et al., 2012;Kéfi et al., 2016;Bastiaansen et al., 2018) and emphasizes the need for the development of models (e.g. Rietkerk et al., 2002Rietkerk et al., , 2004Guichard et al., 2003; Liu et al., 2012Liu et al., , 2014aBastiaansen et al., 2018) to better interpret biogeomorphic changes observed in the field and to help predict possible tipping points.

Persistence of the Vaucheria bedforms
Despite morphological and mechanistic similarities between the observed algal bedforms and self-organized diatom patterns (Weerman et al., 2010), the observed Vaucheria patterns prove to be long-lived. Both the bedforms and the coinciding algal mats remained in place for multiple years, at least from 2012 to 2016 (Figures 3-5 and SI7), but their development might have already started between 2009 and 2012 and possibly continued until 2017. Although algal cover is strongly reduced in winter and/or during hot and dry periods (visually determined, see Figure SI9), their stabilized topography does persist year-round with little change. This observation might imply that bedform stability is maintained by purely geophysical processes when Vaucheria cover is low (e.g. Gouleau et al., 2000;O'Brien et al., 2000;Williams et al., 2008;Carling et al., 2009), or by other (less prominent) algae (Gallagher and Humm, 1981). Viable Vaucheria spores are found even in sediment without visible algal cover ( Figure SI3). In any case, the observation of flattened bedforms (topographic relief halved in September 2017, as compared to August 2015 and December 2016), seemingly following a temporal shift of Vaucheria from ridges to runnels between May and July 2017 ( Figure SI9a), seems to support the expectation that biofilms can significantly contribute to bedform formation and maintenance (e.g. de Boer, 1981;Blanchard et al., 2000;Lanuru et al., 2007;Friend et al., 2008;Weerman et al., 2011;Fagherazzi et al., 2013;Mariotti et al., 2014). This disruption of the spatial biostabilization pattern might be short-term, due to disintegration of Vaucheria filaments at high temperatures (Gallagher and Humm, 1981), or announcing a critical point in bedform self-organization which is reached when the accretionary mudflat is no longer inundated with sufficient frequency (as observed for stromatolites; Logan, 1961). Vaucheria's known potential to form dense filamentous networks which deeply penetrate the sediment makes it less susceptible to erosion than, say, microbial biofilms that form a 'skin' layer Le Hir et al., 2007;Weerman et al., 2011Weerman et al., , 2012. Such a 'memory' effect is reflected in the wavy-crinkly laminae observed deep in the sediment (Figures 6 and SI8), and likely contributes to bedform persistence. Long-term biotic effects on the development of sedimentary landscapes have previously been associated mostly with more complex organisms, such as salt marsh plants (e.g. Temmerman et al., 2007;Schwarz et al., 2014). Hence, improved understanding of the long-term biogeomorphic effects of biofilm-like colonization, which is typically much faster and spatially more extensive than plant colonization, is desirable to understand how widespread and relevant this process is in shaping sedimentary landscapes.

Biofilm influence on bedform lamination
This study suggests that ridge lamination becomes more pronounced and microbialite-like due to multi-annual ridge biostabilization by Vaucheria, which promotes persistent ridge accretion while introducing a seasonal signal to the Figure 6. X-radiographs of the sediment cores (6 × 25 cm selection) taken in September 2016 in bare runnels ('mud') and algal-covered creek levees ('algae'), both taken around the northeastern edge (NE) and the southeastern edge (SE) of the Ketenisse focus area (see Figure SI1 for core locations). Background profiles S 4 (x, y) have been subtracted to obtain these relative grey intensity profiles. Outliers were determined manually for the collective (four cores) data; all cores are displayed with the same resulting range (i.e. excluding outliers: relative grey intensity between À62 and +40). Broad scale (~4 cm wavelength) and fine scale (less than centimetre scale) horizontal lamination is visible in the algal cores, but less evident in the mud cores. Darker vertical lines are signs of bioturbation. [Colour figure can be viewed at wileyonlinelibrary.com] sedimentary deposit. X-radiographs of algal-covered creek levee deposits reveal broad light and dark layers, alternating with a spacing roughly equal to annual accretion (Figures 6  and SI4). This probably reflects seasonality in biofilm growth, environmental conditions or both (e.g. Williams et al., 2008). Fine-scale wavy-crinkly laminae (Figures 6 and SI8), indicators of biogenic activity (Horodyski, 1982;Schieber, 1999;Noffke et al., 2007;Riding, 2011;Sarkar et al., 2016), are most evident in darker layers. These crinkle-comprising layers, whose dark colour might reflect lower density associated with porous, felt-like Vaucheria mats, are most abundant at depths that likely coincide with the years 2011-2013, in which the strongest algal cover is observed (Figures 3 and  6). Layering is much less evident in bare runnel sediments, despite similar year-averaged accretion rates in runnels and on ridges (or creeks and levees), which is in agreement with earlier studies showing that ridge/levee sediments are more gradually deposited and more (a)biotically stabilized, whereas runnel/creek sediments have higher turnover rates (Blanchard et al., 2000;Gouleau et al., 2000;O'Brien et al., 2000;Lanuru et al., 2007;Williams et al., 2008;Carling et al., 2009;Fagherazzi et al., 2013). Indeed, algalcovered creek levees were found to have significantly higher penetration resistance than bare creek levees and runnels ( Figure 7). Although lamination was measured at different moments in time and different locations on the mudflat than penetration resistance, and although the stability contrast might be partly related to the generally higher elevation (and thus lower water content) of algal-covered substrate, these combined measurements strongly suggest a substantial biotic contribution to ridge lamination and stabilization. The fact that lamination is less evident under algal mats closer to the waterfront ('Algae NE' in Figure 6; Figure SI1) might relate to stronger current-and wave-driven disturbances at this lower elevation. Indeed, around the latter location, overall bedform persistence seems shorter ( Figure SI2) than around the location of the clearly laminated algal core ('Algae SE' ; Figures 4 and 6). As sedimentation patterns are inherently more repetitive in bedforms shaped by a consistent feedback than in bedforms shaped by random disturbances, scale-dependent feedbacks might hence stimulate the formation of laminated deposits.

Resemblance to microbialites and wider implications
This study, combined with existing literature, supports the idea that the observed modern Vaucheria bedforms behave in a similar mechanistic and morphologic way to ancient microbialites found in the geological record. The scale-dependent feedback proposed in the current study is comparable to the mechanism described (but not pointed out as feedback) for stromatolites (e.g. Logan, 1961). The observed metre-scale, regular alternation of bare runnels and mat-covered, wavy-crinkly laminated ridges is similar to many layered microbialites (e.g. Logan, 1961;Andres and Reid, 2006;Riding, 2011;Noffke and Awramik, 2013). Vaucheria shares many characteristics with organisms known to form microbialites. Flow deceleration, grain trapping and incorporation, and biostabilization (Noffke et al., 2001;Noffke and Awramik, 2013) are known to occur in Vaucheria mats (and other biofilms) as well (Black, 1933;Webber, 1967;Simons, 1975;Gallagher and Humm, 1981;Paterson, 1994;Skowroński et al., 1998;Wilcox, 2012). Although these organisms are mainly found in siliciclastic environments, calcification has been observed in freshwater Vaucheria mats as well (e.g. Freytet and Verrecchia, 1998;Golubić et al., 2008;Pentecost et al., 2014). Further, fossil algae seemingly morphologically and ecologically equivalent to Vaucheria have been found in rocks of Precambrian age (Butterfield, 2004). Finally, Vaucheria-covered bedforms can resemble the morphology of dome-shaped stromatolites (Figure 2; but note that this photo dates from before the start of the current study and so no detailed measurements could be done). The photo in Figure 2 suggests that biological processes have contributed strongly to shaping these bedforms: they seem to have steep slopes (suggesting biotic importance; Noffke and Krumbein, 1999) and more complex (not purely flow-shaped) morphologies. Moreover, the domed bedform surfaces suggest that ridges formed by accretion (Carling et al., 2009), different from the tabular elevations that remain when a flat, mat-covered surface is erosively disturbed to form erosional remnants and pockets (Noffke, 1999(Noffke, , 2010Noffke and Krumbein, 1999). Therefore, the proposed mechanism of biotically stimulated self-organization due to biophysical scale-dependent feedbacks might well explain regularly patterned microbialites formed by Vaucheria mats and other biofilms, from the Precambrian to present day. Despite earlier reports on microbialite self-organization (e.g. Petroff et al., 2010;Brasier, 2011;Tice et al., 2011;Bosak et al., 2013 and references cited therein), the potential implications of these findings have not been elaborated on thus far (Purkis et al., 2016). If the bedforms and lamination presented in the current study were preserved in the rock record and used for a geological reconstruction, the feedbacks and selforganization mechanism proposed here may significantly affect such reconstructions (e.g. Purkis et al., 2015). The observed ridge flattening (2017) might be preserved in the rock record as horizontally wider, but vertically thinner (or even absent), ridge layers. Layer thinning might traditionally be explained as the result of an environmental change (e.g. reduced largescale sediment availability); a missing layer might complicate layer sequence dating and ridge widening might be interpreted as a wavelength increase related to rising water levels (Franca and Lemmin, 2006;Carling et al., 2009). However, the current study suggests that such layer changes could as well be triggered by a disruption of the internal biophysical feedback, leading to bedform broadening while water depth is shoaling (Figures 3, 4, 5, SI2 and SI7), and without the need for changes in external conditions. That is, the changes in layer geometry can be controlled by autogenic as well as allogenic processes (e.g. Purkis et al., 2016). Further study is required to state with more confidence whether the broadening auto-correlation signals ( Figures 5 and SI7) truly reflect an increase in the entire bedform wavelength or, instead, reflect ridge widening (runnel narrowing) or the disappearance of single ridges while the rest of the pattern stays intact. System-internal processes such as lateral bedform migration or pattern wavelength multi-stability (Bastiaansen et al., 2018) could also be preserved as a seemingly abrupt change in rock lamination. When using laminated microbialites to interpret, date or correlate between rock records (e.g. Logan et al., 1964;Jerzykiewicz and Sweet, 1988;Ndiaye et al., 2012;Harzhauser et al., 2014) it is therefore important to combine lamination analysis with the use of other proxies (such as stable isotopes) for environmental changes in sea water chemistry, ambient temperature, and so on (Dabkowski et al., 2015;Lettéron et al., 2017). Moreover, the investigated stratigraphic sections should ideally be sampled along a transect of equidistant vertical cores to capture at least one entire microbialite wavelength, as possible changes in bedform wavelength or spanwise location provide valuable information. Finally, combining stratigraphic analyses with mathematical models for (self-organized) microbialite dynamics, such as existing abiotic (Grotzinger and Rothman, 1996) and biotic (Weerman et al., 2010) models, could help to solidify the interpretation.

Conclusions
This study reports on present-day biofilms of filamentous Vaucheria algae that stimulate the build-up and maintenance of algal-covered ridges and bare runnels on intertidal flats. These bedforms plausibly self-organize due to local interactions of preferred ridge stabilization and runnel erosion; ridge biostabilization by biofilms most likely amplifies this scaledependent feedback. Flow deceleration and sediment trapping, binding and stabilization by algal filaments promote wavy-crinkly laminated ridge deposition. The proposed formation processes and observed laminae inand regular morphologies ofthese modern biofilm bedforms may mirror the microbialitic deposits that dominated the shallow-marine landscape as early as the Precambrian. This suggestion of microbialite self-organization may also be important for geological reconstructions based on organosedimentary laminae in Earth's early rock record. Figure SI2 -Bottom elevation profiles along the northern edge of the 44 × 49.5 m focus area on Ketenisse. Elevations derived from in-situ laserscan measurements; NAP roughly is mean sea level. Different lines denote different moments in time at which the elevation was measured. The bedform pattern is less persistent than on the southern edge of the focus area ( Figure 4). Figure SI3 -Vaucheria growth potential test for Paardenschor and Ketenisse (see Figure SI1 for study sites). Each bar shows the percentage of samples (5 replicates each) tested positive for Vaucheria growth. On Paardenschor, samples were taken in visually bare sediment at three different elevation ranges. On Ketenisse, samples were taken within bare creeks and on bare mudflat. Samples for positive controls were taken in clearly algal-covered sediments. Negative controls (sterilised sediments) always showed 0% growth potential, as required. Figure SI4 -Spatial auto-correlations along the vertical axis of the relative X-radiographs (Figure 6), after regridding the data to 1 mm resolution. White markers indicate non-significant data points (p≥0.01), coloured markers are statistically significant (p<0.01) based on 100 random permutations. Each correlation profile is afterwards normalised by dividing it by its maximum value, for visibility.  White markers indicate non-significant data points (p≥0.01), coloured markers are statistically significant (p<0.01) based on 100 random permutations. Each correlation profile is afterwards normalised by dividing it by its maximum value, for visibility. Figure SI8 -Close-up image between 11.6 and 16.3 cm depth in the "Algae SE" core in Figure 6, revealing wavy-crinkly laminae (examples indicated by arrows). This image is derived from the unprocessed X-radiograph, i.e. with the original 0.085 mm resolution and without subtraction of the background profiles (which was done to obtain Figures 6 and SI4). Data S1 Supporting information