Observed and modelled tidal bar sedimentology reveals preservation bias against mud in estuarine stratigraphy

Mud plays a pivotal role in estuarine ecology and morphology. However, field data on the lateral and vertical depositional record of mud are rare. Furthermore, numerical morphodynamic models often ignore mud due to long computational times and simplifications of mixed depositional processes. This study aims to understand the spatial distribution, formative conditions and preservation of mud deposits in the intertidal zone of bars in high‐energy sand‐dominated estuaries, and to elucidate the effects of mud on morphology, ecology and stratigraphic architecture. To meet these objectives, field data (historic bathymetry, bio‐morphological maps and sediment cores of the shoal of Walsoorden, Western Scheldt estuary, the Netherlands) were combined with complementary hydro‐morphodynamic numerical modelling (Delft3D). Based on the field observations, two types of mud deposits were distinguished: (1) mudflat deposits, which are thick (>10 cm) mud beds at the surface associated with high elevations and low accumulation rates; and (2) mud drapes, which are thin (millimetre to centimetre) buried laminae that form and preserve at a wide range of elevations and energy conditions. Model results show that deposition on mudflats occurs just after high‐tide slack water in areas shielded from high flood velocities, suggesting that mud accumulation is mostly controlled by elevation, flow velocity and flow direction. Mud accumulation increases shoal elevation, sometimes to supratidal levels. This reduces flow over the shoal, which in turn reduces chute channel formation, stabilises bar morphology and decreases local tidal prism. These effects further promote mud deposition and vegetation settling. Although observations show that mud cover at the surface is relatively high (20%–40% of the intertidal area), mud constitutes only a small percentage of the total estuary volume (ca 5%) revealing that only a small fraction is preserved in the stratigraphy. Due to this mismatch between surface and subsurface expression of mud, interpretations of estuarine stratigraphy risk underestimating the influence of mud at the surface on morphodynamics and habitats.

understand the spatial distribution, formative conditions and preservation of mud deposits in the intertidal zone of bars in high-energy sand-dominated estuaries, and to elucidate the effects of mud on morphology, ecology and stratigraphic architecture. To meet these objectives, field data (historic bathymetry, bio-morphological maps and sediment cores of the shoal of Walsoorden, Western Scheldt estuary, the Netherlands) were combined with complementary hydromorphodynamic numerical modelling (Delft3D). Based on the field observations, two types of mud deposits were distinguished: (1) mudflat deposits, which are thick (>10 cm) mud beds at the surface associated with high elevations and low accumulation rates; and (2) mud drapes, which are thin (millimetre to centimetre) buried laminae that form and preserve at a wide range of elevations and energy conditions. Model results show that deposition on mudflats occurs just after high-tide slack water in areas shielded from high flood velocities, suggesting that mud accumulation is mostly controlled by elevation, flow velocity and flow direction. Mud accumulation increases shoal elevation, sometimes to supratidal levels. This reduces flow over the shoal, which in turn reduces chute channel formation, stabilises bar morphology and decreases local tidal prism. These effects further promote mud deposition and vegetation settling. Although observations show that mud cover at the surface is relatively high (20%-40% of the intertidal area), mud constitutes only a small percentage of the total estuary volume (ca 5%) revealing that only a small fraction is preserved in the stratigraphy. Due to this mismatch between surface and subsurface expression of mud, interpretations of estuarine stratigraphy risk underestimating the influence of mud at the surface on morphodynamics and habitats.

| INTRODUCTION
Over the past few decades, estuary management has widened its focus from flood protection and navigation towards including nature preservation and ecosystem services (Barbier et al., 2011). This development led to an increased interest in mud (clay and silt <64 μm), which plays a pivotal role in ecologically valuable areas. Areas that are low in morphodynamics, with low flow velocities and muddy substrate, are favoured habitats for many estuarine species (vegetation, benthos and therefore birds) (Brückner et al., 2020;Dyer et al., 2000;Gingras et al., 1999, Singer et al., 2016. Mud affects the quality of the bed by attracting nutrients and/or pollutants. When suspended in large quantities, mud increases the turbidity of the water column , leading to reduced light penetration and primary production (Kromkamp et al., 1995). Estuaries with high mud concentrations at the turbidity maximum are prone to develop fluid mud (Baas et al., 2009;Winterwerp & van Kesteren, 2004). Large suspended concentrations also have economic effects, as mud can cause problems through siltation in harbours and navigation channels (van Kessel et al., 2011). Furthermore, mud deposits influence the large-scale morphology of estuaries (Braat et al., 2017. The latter is often neglected in morphological modelling (in this text model is used to refer to numerical simulations) of estuaries due to the need for model simplifications (Hibma et al., 2003;van der Wegen & Roelvink, 2012). Furthermore, very little field data are available on spatial and temporal mud distributions, especially on the scale of entire estuaries. Better understanding and improvement of predictive capabilities of the distribution of mud are necessary for sustainable estuary management.
Sand and mud distribution in the substrate often shows considerable horizontal and vertical variations due to differences in erosion and deposition characteristics of mud compared to sand. Mud needs low velocities to accumulate, but due to its cohesion (van Ledden et al., 2004a) a high critical shear stress for erosion reduces the erodibility of mud, that is the scour lag effect (van Straaten & Kuenen, 1957). Because of these characteristics, the import of mud in estuaries is efficient but the export of mud is relatively difficult.
Mudflats generally develop along the fringes of estuaries and on top of bars (Dalrymple & Choi, 2007;Kleinhans et al., 2021). Mud deposition increases with distance from channels (van Straaten & Kuenen, 1957), which is especially visible in tide-dominated systems with large intertidal areas. Facies descriptions of intertidal deposits from different estuaries, like the Severn (Allen, 1990), Salmon river estuary (Dalrymple et al., 1991;Dalrymple & Choi, 2007), the Bristol channel (Harris & Collins, 1988) and measurements from the Western Scheldt (McLaren, 1994;van der Wal et al., 2010) confirm this typical large-scale distribution of mud in estuaries with mud deposits predominantly at high intertidal elevations (Allen, 1990). In line with these observations, previous modelling research has shown that mud can reduce estuary surface area and width by confining the estuary over centennial to millennial timescales (Braat et al., 2017, similar to floodplains in rivers (Schuurman et al., 2016;Tal & Paola, 2007;van Dijk et al., 2013). Mudflats flank the estuary and limit lateral migration and expansion of channels. Mud deposits also increase the height of intertidal flats, which can enhance the probability of a window of opportunity for pioneer marsh establishment (Cao et al., 2017;de Haas et al., 2018). Vegetation establishment will in turn enhance estuary confinement and mud sedimentation on bars (Brückner et al., 2020;Kleinhans et al., 2018;Lokhorst et al., 2018).
The stratigraphy in estuaries contains information about past conditions of the system. Sedimentation and erosion in estuaries are both spatially and temporally highly variable. Nonetheless, the distribution of mud is predictable. Temporally, erosion and deposition of sediments is possible during one tidal cycle (de Boer et al., 1989), with spring-neap or seasonal cyclicity van der Wal et al., 2010) and during large storm events. Spatially, sedimentation and erosion rates differ largely, for example, between lateral channel migration and vertical bar accretion. This means that the stratigraphic record of estuaries only represents a small fraction of the period over which physical processes occurred and it does not completely capture all processes (Davis, 2012;Jerolmack & Paola, 2010;Paola et al., 2018). Numerical models can help verify processes deduced from the geological record and fill in the gaps. Results from numerical models and the study of currently active estuaries in combination with geological reconstructions allow us to better understand the processes behind the preservation of mud.

K E Y W O R D S
estuarine mud deposits, estuary stratigraphy, mud preservation, shoal accretion, tidal bar morphology Although it is common to include mud in hydrodynamic models, mud and its effects are often neglected in morphodynamic models despite its importance. This is because of difficult calibration, limited modern field data and long computation times, with the exception of some models that break new ground or deviate from the norm (Dam & Bliek, 2013;Le Hir et al., 2011;Sanford, 2008;van Ledden et al., 2004b;Waeles et al., 2007). Models that include mud focus mainly on lateral distribution at the surface and neglect mud present below the surface (Braat et al., 2017;Le Hir et al., 2011;van Ledden et al., 2004b;Waeles et al., 2007). While much evidence of mud beds is found in sedimentary geology (Dalrymple & Choi, 2007;Fietz et al., 2021;Ghinassi et al., 2021;La Croix & Dashtgard, 2014;Martinius & van den Berg, 2011;Shchepetkina et al., 2016;Thomas et al., 1987;van de Lageweg et al., 2018), there is a gap between highly detailed field observations at sparse locations and the relatively coarse system-wide numerical models and maps. Due to this knowledge and data gap, it is still unclear how mud beds form, preserve and affect the morphology of tidal bars. Therefore, ecologists, estuarine managers and geologists have a strong need for field data and morphodynamic models that study mud stratigraphy.
The aim of this study is to elucidate how mud deposits form and preserve in dynamic estuaries. The distribution on a tidal bar and in the subsurface is studied by combining sedimentological field data and numerical modelling. The specific objectives are to determine (1) the locations and conditions of mud deposition on a shoal in a highenergy sand-dominated estuary; (2) how, and over what time span, these mud deposits are preserved in the stratigraphy; (3) the implications for the morphodynamics of the shoal with extrapolation to estuary scale; and (4) the implications for geological interpretations and benefits for ecology. The field data provide information on spatial extent of current and recent mud deposition and mud layer thickness, while the model outcomes reveal the underlying processes that determine mud deposition and preservation. As such, this study bridges the gap between large-scale numerical model outcomes and detailed field observations.

| Site description: the Western Scheldt
The Western Scheldt was chosen as the study area ( Figure 1, the Dutch part of the Scheldt estuary), since it is one of the most well-studied and monitored estuaries in the world, and much historical data on bathymetry, sediment dynamics and hydrodynamics are available. It is the last remaining estuary on the west coast of the Netherlands that has not been (semi-) closed naturally or by the Delta Works, since it is the access point to the Port of Antwerp, the second largest harbour in Europe.
The Western Scheldt developed when the Honte tidal channel expanded landward in the Middle Ages during storm surges, eventually connecting to the Scheldt river in the 17th century Pierik et al., 2017;van der Spek, 1997). Because of this connection, the Scheldt river could drain via the Western Scheldt tidal system, instead of the Eastern Scheldt. The former therefore became deeper and wider (van der Spek, 1997). The estuary had an irregular planform, with secondary branches to which large amounts of fine sediments were imported, causing accretion of tidal flats and marshes van den Berg et al., 1996). During the last centuries, these secondary branches were embanked stepwise, making the estuary's planform smoother, causing the tidal range and average channel depth to increase (van den Berg et al., 1996;van der Spek, 1997;. This increase in tidal energy over the last centuries is clearly reflected in sedimentary deposits; secondary branches are muddy, while modern deposits are predominantly sandy (van den Berg et al., 1996).
The present-day Western Scheldt (Figure 1) is a tidedominated, semi-diurnal, well-mixed, macrotidal estuary with a tidal range at the mouth of about 4 m that increases landward to 5.5 m at Rupelmonde, 110 km from the mouth. The discharge of the Scheldt river is approximately 100 m 3 /s, which is (integrated over 12 h) less than 1% of the tidal prism (2 × 10 9 m 3 ; Wang et al., 2002). Additionally, the Scheldt estuary is alluvial (i.e. developed in loose sediment) and predominantly sandy, with a median sand grain size of about 200 μm (McLaren, 1994). The estuary has a typical exponential convergent shape (Lanzoni & D'Alpaos, 2015;Savenije, 2015) with some variations in width deviating from this trend, typical for alluvial estuaries . Although the estuary has been completely embanked since the Middle Ages and influenced by large amounts of dredging and dumping (Santermans, 2013), there are no dominant constraints on morphodynamics by subsurface geology, except for some erosion-resistant layers underlying the deepest channels (Dam, 2013).
Even though the Western Scheldt is dominantly sandy, mud is essential in this system. Before the onset of significant human interference in the form of dredging , the system imported net 0.5-1.5 million m 3 mud per year, while net 1.4-2.4 million m 3 sand was exported per year from the entire estuary (Dam & Bliek, 2017). Although mud is important to the sediment balance, the total sediment volume of cohesive sediment in the current estuary substrate is estimated to be only 5% (van de Lageweg et al., 2018). On a large-scale vertical resolution (ca 0.1 m, so very thin layers are ignored) the average thickness of cohesive layers is 1.2 m in the Western Scheldt and they are more abundant towards the flanks of the estuary and at the surface (van de Lageweg et al., 2018). Typical suspended particulate matter (SPM) concentrations are between 30 and 60 mg/L (Rijkswaterstaat, 2017; Figure S1), with median settling velocities of 0.1-0.2 mm/s up to a maximum of 0.7 mm/s (Winterwerp et al., 1993).

| Site description: shoal of Walsoorden
Within the estuary, the shoal of Walsoorden was chosen as the study site, representing an average mid-channel tidal bar in a sandy estuary with mudflats and marshes. It is located in the middle of the estuary, approximately 50 km from the coast. The shoal shows strong morphodynamics, possibly showing responses of more stable mud beds over short timescales. The shoal of Walsoorden developed about 30 years ago when two smaller shoals grew together, so the sediments studied are all recent and the morphology is self-formed by the estuarine dynamics. At Walsoorden, water level ranges from −2.5 to 3 m during spring tide and −2 to 2 m during neap tide (station Walsoorden Rijkswaterstaat, 2017). The measuring station of Hansweert, closest to Walsoorden, measures a median SPM concentration of ca 40 mg/L ( Figure S1). Wave energy is almost zero near Walsoorden because of its distance from the mouth (Chen et al., 2005). It is only influenced by locally generated wind waves, which can have a large effect on the morphology of tidal bars (Maan et al., 2018), but the fetch of the dominant wind direction is very small.
In this research, the term shoal is used in a geographic context, like the study area the shoal of Walsoorden, while the term tidal bar is used to describe the study site in a stratigraphic and geological context. The term mudflat is used to describe the smaller-scale, more surficial morphological unit consisting dominantly of mud in the intertidal zone. By this definition, mudflats can occur on top of shoals and are part of the tidal bar.

MATERIALS
This study uses field data obtained from the shoal of Walsoorden in the Western Scheldt and a complementary numerical hydro-morphodynamic model in Delft3D (Figure 2). The distribution of mud based on field data of the past 20 years is analysed and additionally used to validate the model. The model elucidates the conditions and processes of deposition and preservation, as the field data do not provide information on hydrodynamic conditions. In addition, the temporal resolution of the field data is too low to capture erosion processes on the timescales of a tidal cycle of spring-neap cycle.

| Use of existing data to map surface and subsurface mud distribution, shoal of Walsoorden
Ecotope maps were used to quantify the spatial distribution of mud over time, these are bio-morphological maps that were readily available for the years 1996, 2001, 2004, 2008, 2010, 2011, 2012 and 2016 (Table S1). These maps were constructed through visual classification of aerial photography with ground truthing and tidal zonation based on bathymetry (Paree & Burgers, 2017). The maps contain mud-rich sediment classes and vegetation classes, which were grouped in analysis as muddy sediments. The vegetation classes were included because field observation showed that mud-rich sediments occurred where vegetation grew.
Bathymetries were used of the years that correspond to the availability of the bio-morphological maps. The gridded bathymetries are a combination of echo-sounding (vaklodingen; 20 m resolution) and laser altimetry on the intertidal areas (2-5 m resolution). For 1996, the bathymetry was solely based on echo-sounding since no laser altimetry was available.
Besides the maps, three profiles of the bar above the low water line were measured regularly from 1991-2017 using Differential GPS. On these profiles, a few fixed points are located from which soil samples of the top 2 cm and top 10 cm were regularly analysed for mud content ( Figure S2). All bathymetry, profiles and bio-morphological maps were made available by Rijkswaterstaat (Dutch Water Authorities).
An interpolated stratigraphy was reconstructed from the existing field data of the last 20 years in MATLAB (similar to van de Lageweg et al., 2016). A mud and sand class were generated from the bio-morphological maps F I G U R E 2 Conceptual figure of the methodology, which combines field data and a numerical morphological model including mud and stratigraphy and coupled to the corresponding bathymetry. From old to new, these map layers with elevation, lithology and age were stored and combined into an interpolated stratigraphy (varying thickness per location). When deposition occurred, the older, underlying stratigraphy was preserved. However, when erosion occurred, the interpolated stratigraphy present above the new bed level was removed. Note that Supplemental Material Figure 2 uses all available DGPS measured transects of multiple measurements a year, while the transects in Figure 6 are based on a combination of LiDAR and DGPS elevations and are only shown for times that bio-morphological maps were available.

| Collected field data on subsurface mud distribution, shoal of Walsoorden
Although abundant coring data exist of the Western Scheldt (stored in DINO database, as shown in van de Lageweg et al., 2018), only one core, with low vertical resolution, is digitally available for the shoal of Walsoorden. Therefore, fieldwork was carried out in October 2017 to gather detailed subsurface data of the shoal of Walsoorden. These data were then combined with the interpolated stratigraphy in figures, for analysis of the mud distribution in the subsurface. During the fieldwork 36 cores (black dots in Figure 1B) up to 3.1 m depth were taken over a wide variety of environments on the shoal, including marsh, low-energy and high-energy tidal-flat environments ( Figure 3). Cores were made using a gouge auger, a Van der Staay suction corer and a larger suction corer for sampling. Outcrops on the edge of the salt marsh were cleared using a shovel and then described. The cores and transects were photographed and described lithologically in the field. The GPS coordinates (not elevation) of the sample locations were acquired by handheld GPS and phone GPS. Throughout this study, the Dutch national coordinate system (RD coordinates, Dutch: Rijksdriehoeksmetingen), expressed in kilometres, was used.
Four cores 45-70 cm in length, were collected near a large flood channel, then further analysed in the laboratory as these locations were expected to show preserved small-scale mud deposits and diatoms. These data were used to determine the conditions at the time of deposition rather than inferring the living environments (e.g. salinity) at the core location in the past. For example, a lot of marine/saline species imply that a storm or spring tide transported them far into the estuary, since they do not typically occur at this site, thereby implying decreased influence of discharge. The cores were studied in more detail by diatom analysis and sedimentological analysis from lacquer peels ( Figure 4). The lacquer peels were made using a colourless 210J Flits Coating (similar to the method used in Martinius & van den Berg, 2011)  Figure 4). The lacquer peels revealed sedimentary structures of the cores that were not visible in the photographs. For diatom analysis, five mud layer samples were obtained from two cores (white circles in Figure 4) that were hypothesised to have a different sedimentary origin based on thickness, sequence, colour and organic material. In the diatom analysis, different functional groups were identified (benthic, tychoplanktonic and planktonic) and some species and genera indicative for specific environments. The screening did not involve any counting of the diatoms. The model includes two sediment types: one sand fraction (200 μm) and two mud fractions. The mud characteristics of the fractions are the same (Table 1), but they are supplied from different boundaries (river or sea side), to present marine and riverine mud, so the origin of the mud deposits can be tracked. Sand transport is calculated using van Rijn (2007aRijn ( , 2007b: TRANSPOR2004 equation and mud erosion and deposition is calculated with the Partheniades-Krone formulations (Partheniades, 1965, Winterwerp & van Kesteren, 2004; Equations 1 and 2).

| Model description and set-up
where E is the erosion flux (kg m −2 s −1 ), D the deposition flux (kg m −2 s −1 ), M the defined erosion parameter (kg m −2 s −1 ), S(τ cw ,τ cr,e ) is the erosion step function, S(τ cw ,τ cr,d ) is the deposition step function, w s the settling velocity (m/s) and c b the average sediment concentration. Winterwerp et al. (1993) showed with flume experiments that the shear strength of a Western Scheldt sediment mixture with 70% mud is approximately 0.1 N/ m 2 just after deposition. However, the critical bed shear stress for erosion is variable in reality, most importantly due to compaction. To account for some compaction, a critical shear stress for erosion of 0.2 N/m 2 (Table 1) was chosen. A very high critical bed shear stress for mud deposition was chosen, so continuous settling occurs (Sanford & Halka, 1993) (Table 1). To keep track of sand and mud in the bed, the underlayer module was used. Although this module in Delft3D has been available for a couple of years, it has not been widely used, yet (van Kessel et al., 2012). The underlayer was used with a Lagrangian active layer with a fixed thickness and Eulerian storage layers (van Kessel et al., 2012) (Table 1). There is no interaction between sand and mud in the bed, erosion and deposition is handled separately. Bed level change is calculated from the divergence of bedload sediment fluxes and erosion-deposition differences for suspended sediment. To speed up morphodynamic calculation, the bed level changes are multiplied with a morphological acceleration factor of 20 each time step (Roelvink, 2006).
The model has a curvilinear grid that extends from Bakendorp to the Dutch-Belgian border between the dikes (approximately 25 × 15 km). It is decomposed into two domains: the outer coarser areas with a median cell size of 64 × 44 m, and the inner refined domain of Walsoorden with median cell size of 33 × 40 m ( Figure 1C). Bathymetry data from Rijkswaterstaat (vaklodingen) were used for the initial morphology.
One year of water level data from 2013 obtained from the full NeVla model (van Dijk et al., 2019) were used at the two boundaries of the model including astronomic tides, storm surges and wave set-up. The same boundary conditions were used for all model runs. The original NeVla model has already been calibrated for water levels and velocities and shows a good correlation with measurements (Schrijvershof & de Vet, 1988); the model therefore does not need any further calibration for hydrodynamics. Sediment conditions at the boundaries Morphological acceleration factor 20 -are in equilibrium for sand, meaning that the transport gradient perpendicular to the boundary is zero. Mud supply is a constant concentration of 40 mg/L (Table 1, Figure S1). The reference scenario is a long-term run of 20 morphological years, based on 1 year of hydrodynamics. The starting year, 1996, corresponds with the oldest available bio-morphological data set. The initial bed composition (fraction of sand and mud in the bed) is from 1994 (McLaren, 1994) and is vertically uniform. The mud class assigned to the initial bed is the same as that which enters at the marine boundary, defined as marine mud. Sensitivity tests were conducted to obtain the best settings for mud and morphological parameters (mainly settling velocity, critical bed shear stress for erosion, erosion parameter and mud supply concentration: Table 1) that reproduced mud patterns observed in the bio-morphological maps within the range off values suggested by literature (Table 1). About 100 tests and scenarios were conducted when preparing the reference scenario, but were not included in this paper.
Model results are visualised similar to the field data in Figure 5, to allow for straightforward comparison. Since the modelled morphology becomes increasingly less comparable to the real bathymetry over time, six short-term models were also run with different initial bathymetries of 75 hydrodynamic days (4 morphological years) to validate if the mud deposition pattern holds for different initial morphologies ( Figure S3). In addition to the longterm and short-term morphological runs, a run of one tidal cycle was done to capture potential dependencies of the mud deposition and erosion over one tidal cycle.

| Limitations
Despite both field data and modelling having their limitations, the strength of this research is their combination that allows novel insights to be generated into process and system understanding. Field interpretations are limited in terms of spatial and especially temporal resolution. Field stratigraphy provides the results of the processes of interest, while the processes themselves can only be inferred from the observed end result. In this study specifically, the mud maps lacked seasonal information because the biomorphological maps are always made at the same time of year, and the temporal resolution (1 map per ca 4 years) determined the resolution of the interpolated stratigraphy. The cores grant more subsurface detail but only for the situation at the moment of collection, as not all of the sediment may have been preserved in the cores and later erosion may change the stratigraphy of the active system in the future.
Models, on the other hand, are limited by simplifications of physical laws underlying the morphological processes. In addition, in this study, the effects of salinity variations on mud dynamics, 3D circulation, wind, waves, flocculation and soil compaction are ignored to maintain achievable computational times. It is assumed that these simplifications are reasonable because the estuary is well-mixed and tidal processes dominate sediment transport while the wind fetch is short and the study site is 50 km land inwards from the coast. Bed composition dependent roughness is also partly excluded, as a spatially variable, but time fixed roughness calibrated the NeVla model. Because physical processes influence the key attributes of mudflats much more than biological processes, this can be justified (Dyer et al., 2000). Despite the simplifications, the model results provide a complete coverage of the area, include many parameters and are of very high temporal resolution. So, when field data and morphodynamic models are combined, the model results can complement limited process understanding from field data. In turn, the field data provide much more vertical detail and corroboration of the model results.

| RESULTS
The shoal morphology is presented first, followed by the spatial mud distribution and then the bar stratigraphy. The results of the field cores are shown in the Figure S4 and added to the interpolated stratigraphy in Figure 7.

| Morphological evolution
The field data reveal morphology changes over two to three decades. The south-eastern tip of the shoal lengthened and narrowed between 1996 and 2004. After 2004, the shoal length decreased again, while narrowing continued ( Figures 5A through F and 6A,C). At the northwestern tip, a flood channel started to develop from 1990 that eroded the northern flank of the shoal but developed a partly attached parallel bar on the north side ( Figures 5A  through F and 6B). Additionally, the field data show that the southern flank of the shoal has a more irregular trend of erosion and deposition compared to the other sides of the shoal, which is shown by the crossing lines on the left side of Figure 6B The morphological evolution in the model (for 20 years, starting in 1996) shows similar initial trends as the field data, but diverges increasingly with time. The model reproduces the lengthening and narrowing of the southeastern side of the shoal ( Figure 5G through L). However, the model does not show the erosion of the south-eastern tip after 2004 (transect 3 Figure 6) nor the north-northeastern flood channel and small bar formation ( Figure 5). Despite differences in the detailed morphology, the overall morphological evolution is still realistic.
The short-term model runs (4 years; Figure S3) show a greater similarity of mud deposits with the field data compared to the reference scenario, because the underlying morphology is more similar to the field data. These shortterm runs demonstrate that the initial mud concentration of the bed is not very important for the mud layer development and the model predicts mud deposits well on different initial morphologies compared to the field map data.
Both the field data and the long-term model show continuous vertical accretion of the shoal (Figure 6). Accretion up to 0.5 m is observed for both data types in transect 4 between 1996 and 2016 ( Figure 6), which elevates certain intertidal areas up to mean high water level. Additionally, the plan-view size of the shoal fluctuates slowly within the same range for both data types (between 4.5 and 6 km 2 ; Figure 8A,B).
The differences between the model simulation and the field data can be attributed to inaccuracies in the sediment transport predictors and simplifications of mud sedimentation processes and storage in the model stratigraphy. Furthermore, human interference in the model is ignored, most importantly dredging and dumping. For example, in the large model domain the formation of a submerged bar is observed at Hansweert, but in reality this is an important dredging location (Plancke et al., 2014;van Dijk et al., 2021). The model was run without dredging, and bars formed at this location which influence the channelshoal arrangement and evolution of the river direction (also observed in van Dijk et al., 2019), including the shoal of Walsoorden. This problem at Hansweert was not dependent on initial bathymetry, as this was tested in model scenarios with other initial bathymetries. Despite these differences with the field data, the bar complex shows realistic morphology with natural slopes, channels and mud cover.

| Spatial mud distribution
Mud is mostly located at the south-eastern side of the shoal ( Figure 5A  time span. Vertically, mud deposits are only observed above mean water level (+0 m NAP; Figure 7), with most deposits occurring between neap and mean high water level. The vegetation classes in the field maps (green; Figure 5A through F) indicate that the vegetation presence on the shoal of Walsoorden increased over time. The field cores and observations showed that vegetation strongly correlates with thicker mud deposits compared to locations without vegetation (cored; Figure 7). The distribution of mud in the model simulation is comparable to the field data ( Figure 5G through L). Modelled mud deposits also occur on high intertidal elevations, dominantly at the eastern side of the shoal and are constant over time. , mud in the model also spread to the middle of the shoal before it retreats again slightly in 2012. However, in the field data, mud does not spread towards the centre, but more towards the northern side of the bar where vegetation establishes. As a result, the mud within the vegetated area at the north side is under-represented in the model ( Figure 5D,E,F compared to J,K,L), especially because the thickness of this deposit in the field data is quite significant (transect 2, Figure 7). This difference is probably attributed to vegetation roughness in the field that stimulates mud deposition (Brückner et al., 2020;Mudd et al., 2010). The difference in mud layer thickness in this particular location suggests that vegetation might have established itself before mud started to accumulate.
In general, the total mud cover in the field data is smaller than the modelled mud cover ( Figure 8C,D). This is not a real but an apparent difference, because the E model calculates mud cover as a mud fraction in the bed, whereas the mud cover by field data uses a simple present/absent rule. Small percentages over the entire shoal are therefore included in the mud cover calculations in the model, while they would not be classified as mud from aerial observations or during fieldwork. The numerical model results indicate that mud deposits originate mostly from the seaward boundary and are therefore of marine origin ( Figure 8D). In total, the model results indicate that 1.64 × 10 6 m 3 mud was deposited in the Walsoorden model domain between 1996 and 2016. Of this volume, 87% originated from the seaward boundary (or reworked initial mud) and 13% from the river boundary. The amount of mud entering the domain during flood is about 8500 m 3 compared to 2250 m 3 during ebb, so 79% of the sediment entering the domain originates from the seaward boundary. When comparing these percentages, there is a slight preferential settling of mud from the seaward boundary.
The total mud surface cover, including the vegetation class, doubled between 1996-2001 as shown by the field maps ( Figures 5A,B and 8C). After this peak, a small, slow decrease is observed. This trend is also visible in the model results ( Figures 5G,H and 8D). Since the model uses only 1 year of water level data to calculate 20 years of morphology, this trend is not a result of hydrological or ecological cycles. Therefore, this trend in mud distribution is probably coupled to the large-scale morphological change of the shoal. The increase in mud cover is attributed to the heightening of the shoal and the decrease to the narrowing of the shoal. In addition to the long-term trend, strong spring-neap cyclicity is also observed in the model ( Figure 8D), which could not be observed in the field data due to low (supra-annual) temporal resolution ( Figure 8C).

| Mud in estuarine deposits
Focussing on the vertical distribution of mud, it can be seen that most of the subsurface consists of sand (Figure 7 and Figure S4). In the morphodynamic model, mud accumulation significant enough to form layers only starts to occur in the intertidal zone and thick mud layers (>10 cm) are found only at elevations near the high water level. This mud deposition in the model is consistent with field data (Figure 7). Calculated from the model, only 5.7% of the sediment in the bed that was reworked over the 20 years at the shoal of Walsoorden is mud (specifically in Delft3D this means all layers except the base layer). Based on the cores, two types of mud deposits are distinguished: (1) Mudflat deposits, which are thicker beds (>10 cm) that occur at the surface; and (2) Mud drapes, which are thin laminae of millimetres to at most a few centimetres (Martinius & van den Berg, 2011, Figure S4). The two types of deposits are associated with different settling and preservation conditions.
The thicker mud beds, interpreted as mudflat facies, only occur at high intertidal elevations on top of the shoal, both in the field data and in the model simulation. They form by slow accumulation over time. The deposits can be as thick as 1 m ( Figure 7B,F,H). The thickest beds were found in the vegetated marsh; however, the model shows that similar thicknesses can also be reached without vegetation. No evidence was found of mudflats buried by sand.
The thickness of the mudflat deposits on top of the shoal is a good representation of the height by which the shoal has increased. When the reference model scenario is compared to a model without mud supply, a difference in shoal elevation is observed ( Figure 9B,C). Where the model with mud accreted, the model with only sand maintained its initial elevation. The increase in shoal elevation is therefore not forced by long-term changes in hydrodynamics, but by mud settling on the high intertidal areas. When these insights a transferred to the field data, which also shows an accreting shoal ( Figure 9A), it is possible to conclude that this is also caused by mud settling on the shoal. Other effects of mud on morphology, for example estuary confinement, are limited at this spatial and temporal scale. The mud drape data come exclusively from the field data, as the drapes are too thin  to be captured in the model stratigraphy. The cores and lacquer peels show fine sands, which are generally laminated with millimetre thick mud drapes. These are interpreted as tidal bar deposits in the upper part of the intertidal environment, which is confirmed by their current surface position on the shoal. Where the saltmarsh is present, thick (>10 cm) rooted mud beds occur on top of these deposits. The sand-dominated cores contain mud laminae, which are occasionally rich in diatoms. Samples 1 and 4 from cores 37 and 39 show a higher abundance of planktonic species over benthic species (Table 2). Possibly, the planktonic diatoms were transported to this location under strong marine influence. Marine inwash is common in estuarine environments (Vos et al., 1988), and mostly occurs during high-energy events. This is confirmed by the occurrence of high foraminifera abundances and coarse sand grains in these samples (and in one of the samples also broken shells). This seems to confirm that these drapes were formed during storms or exceptional high spring tide. Also diatoms preserved during more quiet conditions were found. In samples 2, 3 and 5, benthic species dominate (Table 2). Benthic diatoms are unlikely to live in high dynamic and subtidal areas, so the mud in these drapes was probably deposited under calmer hydrodynamic conditions. This is certainly true for sample 2, where root remnants were found indicating vegetation settling; sample 3 did not contain many diatoms. Sample 5 contains clay pebbles which probably represents reworked cohesive clay layers (Figure 4). This is supported by a high abundance of diatom fragments and the presence of a (sub)-aerophilous genus that is unlikely to occur in the wet tidal environment. These pebbles and its diatom assemblages suggest that these pebbles were eroded from a salt marsh environment (similar to the one now on top) and deposited in the bar. The analysis of these five mud drapes indicates that they formed during slack tide under a variety of conditions, locally sheltered estuarine and storm conditions. Burial probably occurred in all cases due to bedform migration or event-related sediment transport under dynamic conditions and greater water depths.
The formation and preservation of mud drapes and mudflat layers are related to different processes (rapid or accidental burial versus slow accumulation). Deposits with thicknesses between these two classes (e.g. a few centimetres) are uncommon and can be considered either temporary (seasonal or spring-neap related) or at an intermediate stage of becoming a thicker deposit. Herman et al. (2001) also considers the mud deposits in the Western Scheldt in the top 10 cm of the bed (dated by 7 Be radionuclide tracer) as temporary due to seasonal cyclicity. Deposits that are typically thinner than 10 cm, mostly accumulate in summer and erode in winter and are therefore seasonal.

| DISCUSSION
First, the spatial distribution of mud on the shoal of Walsoorden is discussed and compared to fluvio-tidal bars in other estuaries. Second, the preservation of mud in stratigraphy is analysed and the implications for geological interpretations are outlined. Last, the implications of these findings are discussed in terms of shoal morphology and marsh vegetation.

| Mudflat formation on tidal bars
Modelling results and field data show that mud on the shoal of Walsoorden is most successfully trapped on the eastern side of the shoal at high intertidal elevations. Figure 10 illustrates why thick deposits, that is mudflats, only form at high elevations. Mud flat deposition during low-tide slack water (below low water level) has a high chance of being re-suspended during peak flood flow. However, mud deposited at high elevation during hightide slack water falls dry during peak ebb and flood velocities and is therefore protected.
On average, the Walsoorden model domain is flood dominant (Figure 10). The short-term numerical model shows that the peak velocities occur in the flood phase, just before maximum high water. In this period, the surface area of mud decreases by ca 2% compared to the maximum surface area of mud during ebb just after high-tide slack water. The magnitude of the tidal cyclicity in erosion and deposition is exaggerated in Figure 8D by the morphological acceleration factor: the response of the morphology to the hydrodynamics is faster due to the acceleration factor. Nonetheless, these velocity fluctuations explain the cyclicity observed in the modelled mud area cover on the scale of individual tidal cycles ( Figure 8D).
The preferential settling of mud on the eastern side of the shoal cannot be explained by elevation difference alone. The velocity magnitudes on both sides of the shoal are different for equal elevations. This difference can be explained by flow direction, energy dissipation over the shoal and flow divergence at the top of the shoal. As the highest velocities occur during flood, the east side of the shoal will be largely shielded from the flood flow ( Figure 10). During the ebb phase the velocity is smaller, with its peak at lower water levels.
Because most mud settles just after high-tide slack water, when a lot of water and sediment just entered the area from the western boundary, a larger portion of the deposited mud is of marine origin. Riverine mud entering the domain from the eastern boundary only dominates at low tidal water levels when the flow is localised in the main channels.
It is expected that these findings will be comparable to other tidal bars in tide-dominated, high-energy, sanddominant estuaries under flood-dominant conditions. Small differences might occur closer to the mouth of the estuary or in systems with high-energy wind waves. With regard to the mud accumulation under the influence of waves, it is expected that mud deposition is overestimated because waves are ignored Maan et al., 2018), but mud deposition is underestimated because biota is ignored (Weerman et al., 2010).
A similar mud distribution was found on top of the shoal named Paap in the Ems estuary (van Straaten, 1960). This estuary contains more mud, is less morphodynamic than the Western Scheldt (van Maren et al., 2016) and therefore represents a good system to compare to the shoal of Walsoorden and analyse if these findings hold for muddier systems. Spatially, the middle and higher parts of the Paap shoal show higher fractions of 20%-40% mud, while the sides are sandy with only 2%-3% mud (Wiggers, 1960). In addition, sediment samples from the Rijkswaterstaat (2009, 800 m resolution of the whole Ems estuary from 1989) also indicate a higher concentration of mud on the middle and higher part of the Hond-Paap shoal. The channel on the south-western side of the shoal is filled with mud and therefore samples towards the channel will be even muddier than on the shoal. These types of filling deposits were not observed in the Western Scheldt model, but demonstrate another mechanism of the formation of thick mud deposits in estuaries, 'residual tidal channel mud fills', which does not apply to this case.
In contrast to the shoal of Walsoorden, the Hond-Paap shoal in its current state does not show a general correlation between mud and topography. It is thought that this is due to lateral sorting, which does not occur on the shoal of Walsoorden. It is hypothesised that the current shoal has grown so large after the merger of the Hond and Paap shoal (ca 4 km wide compared to ca 1.5 km for the shoal of Walsoorden and 36 km 2 compared to an estuary size of 460 km 2 ), that lateral sorting started to occur on the bar similar to floodplain sedimentation, with fining away from the channel. The Hond-Paap shoal takes up more than 50% of the estuary width (de Jonge et al., 2014). Sand settles relatively close to the shoal edge and mud deposits in the centre that now has a lower elevation than the edges of the shoal. It is hypothesised that lateral sorting might also be important for significantly larger shoals in the Western Scheldt, which would require modelling with several mud fractions.

| Preservation of mud and geological implications
Estuaries are often considered very efficient sediment traps over long timescales (10-100 kyr), as they often develop in confined and deeply incised palaeovalleys under sea-level rise (i.e. transgressive systems) (Dalrymple et al., 1992;Demarest & Kraft, 1987). However, in estuaries formed under low rates of sea-level rise (i.e. the highstand situation), such as the Western Scheldt, the preservation potential of sediments is much lower; because channel-shoal migration laterally reworks older deposits (Pierik et al., 2016), whilst relatively constant sea level limits the stacking of sediments due to limited vertical accommodation space. On top of that, embanking has increased the tidal energy in the Western Scheldt. While the preservation potential in the Western Scheldt is low compared to prograding systems like deltas or other systems that experience sea-level rise (transgressive systems), there is large spatial and temporal variability in preservation (Martinius & van den Berg, 2011). On the shoal of Walsoorden, net sedimentation on top of the shoal, on the order of centimetres per year, is mostly mud according to the model results corroborated by field data (Figure 7). These deposits are laterally eroded by the main channels migrating at several metres per year. Because the Western Scheldt is a very dynamic system, bars are constantly eroded and rebuild on timescales of decades to centuries. With only one moment of subsurface measurements or geological data, only the last channel-bar configuration can be studied. So when a system has been active over many centuries, the stratigraphic record consists mostly of sediment that was deposited during an extremely short period (years to decades) due to constant sediment reworking. Nevertheless, the stratigraphic record contains valuable information on past conditions but should be studied with care (Davies et al., 2019;Jerolmack & Paola, 2010;Paola et al., 2018), especially regarding interpretations related to mud deposits as will be explained in this section.
As a result of sand-mud segregation in tidal bars, the internal architecture is found to be predictable: dominantly built of sand with mud beds on top. Van de Lageweg et al. (2018) also found 5% of mud in the subsurface, compared to 5.7% in this study with deposits fining upwards in both studies. The model and field data presented here demonstrate that it is difficult to preserve thick mud beds (>10 cm) in deeper subtidal deposits due to resuspension after low-tide slack water ( Figure 10). Even though the mud volume in the bar and estuary is small, the surface expression is large: 20%-40% of the intertidal area is covered by mud. As such, the volume of mud might seem small when studying estuary stratigraphy, but this mud volume is not a good representation of the importance of mud in the estuary, as mud at the surface has a strong influence on morphodynamics and habitats, and is still a significant contributor to the sediment budget. Thus, a preservation bias exists. Most mud is deposited at the surface and is therefore more easily reworked than sand in the next high-energy event. In addition, on geological timescales, transgression has a strong potential to erode the top beds. Therefore, the importance of mud in the system should not be underestimated based on subsurface data.
The field data did however show very thin mud drapes (2-20 mm) in deeper stratigraphy ( Figure S4). These drapes were found at a wide range of depths and at locations with high and low energy conditions near channels. A mud drape can be deposited during every slack water phase, given that the depositional environment is water covered. However, the deposits can only be preserved if they are either very rapidly covered by sand before the drapes are eroded or if they are deposited in an area shielded from flow, for example in bedform troughs or in pools during low tide (Martinius & van den Berg, 2011). This requires a higher energy environment than on the top of the shoal, but should not be so high that net erosion occurs over a tidal cycle. These findings are also represented in the results of the diatom and lacquer peel analyses. So, counter-intuitively, the mud drapes preferentially preserve under high nonerosive energy conditions. Preservation at the shoal scale of the model simulations shows that over the 20 years 52% of the deposited sediment was reworked within the region of interest, using an observational time step of 1 day. This implies that all sediment that settled at some point during the simulation still seems to be at the same place by the end of the simulation. This percentage is an underestimation of reality, because preservation and reworking calculations are dependent on the observational time step (Figure 11). Within a timestep of 1 day, erosion and re-deposition can occur during two tidal cycles, which can be overlooked without intermediate time steps. A model simulation with shorter time intervals (0.2 min) shows that the percentage of reworked sediment was as high as 58%, which is probably close to reality. This insight is a step forward in understanding the classical problem of time and preservation in geological data (Davies et al., 2019;Paola et al., 2018;Trabucho-Alexandre, 2014). Unlike for rivers, preservation potential in estuaries cannot be calculated from sedimentation relative to the influx of sediment, because of tidal flow that reworks the sediment daily. Moreover, in the model simulation a morphological acceleration factor is used, thus creating a mismatch between the depositional volume of sediment and the actual bed level change. Vice versa, volume budgets from the geological record cannot be directly used to estimate sediment fluxes or suspended concentrations, as constant sediment recycling occurs on very short timescales.
Summarised, shoal stratigraphy of high-energy, sanddominated estuaries reflects the spatial distribution of mud deposited over the last 10-100 years. However, mud contents or occurrences do not directly relate to hydrodynamic processes. It is crucial to consider the preservation mechanisms and preservation bias related to stored mud volumes. First, the overall preservation in estuaries is limited compared to for example rivers and deltas, due to constant reworking of sediment, and the temporal resolution of data influences the interpretations regarding preservation. In addition, because mud is predominantly deposited on the top of the shoal, mud cover and its morphological and ecological importance at the surface are underestimated if based solely on cores and outcrops.

| Implications for shoal morphology and ecology
The largest influence of mud on the estuary morphology is shoal accretion (Figures 7 and 9). Mud increases the elevation of the shoal by accumulating at high intertidal elevations where sand is hardly transported due to low flow velocities and short inundation times. By the increase in elevation the flow over the shoal is further reduced and starts a positive feedback for mud deposition (Braat et al., 2017. The accretion and reduction of flow over the shoal increases resistance to the tide propagating through the estuary, which can lead to local reduction in tidal prism , further reducing velocities and increasing mud deposition. The tidal prism reduction was not observed in the Walsoorden model domain, because it was masked by deepening of the main channel which increases tidal prism. In addition, flow reduction over the shoal leads to a reduction in small tidal channels forming on top of the shoal. Due to less tidal channels and velocity reduction, the bar becomes more morphologically stable and is less probable to be wiped clean of mud, or even split by cross-cutting channels. The positive feedback on mud deposition is only limited by decreasing accommodation space and inundation times resulting in less time for deposition. Because tidal bar elevation increases with mud deposition, the low-dynamic mud-dominated surface area increases. These are very valuable ecological areas for benthic species, birds and fish (Bouma et al., 2001;Gingras et al., 1999;Meire et al., 2005). If the elevation of the shoal increases up to high water level, these low-dynamic intertidal areas might be lost to the formation of supratidal salt marshes, accommodating different species. Established vegetation species experience less inundation stress and new seeds will establish more easily due to lower flow velocities (Bouma et al., 2014;Brückner et al., 2020;Cao et al., 2017;Lokhorst et al., 2018). However, vegetation patterns are not limited to muddy substrates, possibly due to species specific preferences (Brückner et al., 2020). As discussed previously, at some locations of the shoal the mud deposits were not well-reproduced in the model, because in reality these were locations where vegetation was probably established before mud accumulation. Nonetheless, once vegetation is established, it reduces flow velocities and increases mud layer thickness (Brückner et al., 2020). Inclusion of biota in models would further improve predictions of mud distributions in the future models (Brückner et al., 2020). It is expected that mud deposition is currently underestimated by ignoring biostabilising effects and overestimated by ignoring waves, and the net effect is expected to be seasonal and of minor influence to the morphology (van der Wal et al., 2010). Understanding the spatial and vertical distribution of mud layers will improve future modelling studies of estuaries. This study gives a methodology to account for mud and a range of parameters that can be used to obtain realistic results verified with field data. Especially for the Western Scheldt and similar estuaries, it is shown that although the total volume of mud in the bed is very small (5% of the sediment), it has significant effects on shoal elevation, ecological area and the sediment balance as the surface cover is much higher than 5% (namely 20%-40% of intertidal area). The successful modelling of these mud deposits will also allow assessment of possible effects of changes in tidal range, freshwater influx from the river and sea-level rise under climate change and human interference. The detailed sedimentological data of mud could have relevance for outcrop interpretation, permeability of the substrate and ecosystem development in high-energy estuaries. A better understanding of the processes that control deposition and preservation is provided, which supports future geological and ecological studies on estuaries.

| CONCLUSIONS
The objective of this study was to better understand the locations and conditions of mud deposition and preservation on a shoal in a high-energy, sand-dominated, laterally constrained estuary.
The morphodynamic model showed that the locations of mud deposits are mainly determined by a combination of elevation, flow velocity and flow direction. High intertidal areas that are shielded from the peak velocities during high water form optimal conditions in high-energy systems for mud deposition. Most mud on the shoal of Walsoorden therefore accumulated at high intertidal elevations on the south-eastern side of the shoal and originated from the seaward boundary.
Field data showed that almost no mud is preserved in deeper stratigraphy. Two types of mud deposits were distinguished that are associated with different locations and processes: (1) mudflat deposits, which are thick (>10 cm) mud beds on top of the shoal formed by slow (multi-year) accumulation over time; and (2) mud drapes, which are thin (millimetre to centimetre) deposits that form and preserve more rapidly (presumably within one tidal cycle) at a wider range of elevations and energy conditions. The thin drapes have negligible influence on the morphology, while the thick beds at the surface increase shoal elevation and are a prerequisite to rise to supratidal levels. Deposits of intermediate thicknesses are considered temporary (seasonal or spring-neap related) or in transition to become mudflat. An increase in shoal elevation by mud increases the amount of low-dynamic high intertidal muddy areas that are valuable for ecology and stimulate salt marsh growth on the shoal raising the shoal to supratidal levels. Further, it reduces cross-cutting by small tidal channels and reduces the tidal prism locally and further upstream.
Stratigraphic mud content or occurrences can only be related to hydrodynamics, morphodynamics or surface ecology with consideration of the preservation mechanisms and biases. Because most mud is deposited at high elevations, it is more easily reworked than sand in highenergy events and rarely preserved deeper in the bed. While only a small fraction of the stratigraphy consists of mud (here ca 5%), the surface mud cover is much higher (here 20%-40% of the intertidal area). Therefore, the effects of mud on the ecology and morphological evolution of the system are probably underestimated if these are based on limited geological outcrop or subsurface data.

ACKNOWLEDGEMENT
This research was funded by the Domain of Applied and Engineering Sciences TTW (grant Vici 016.140.316/13710 to MK) of the Netherlands Organisation for Scientific Research (NWO) and is part of the PhD project of LB. HJP and MB were funded by the ERC Consolidator project 647570 awarded to MK. The authors thank Rijkswaterstaat for the sharing of data and their support during the fieldwork. In particular Marco Schrijver, Gert-Jan Liek, Robert Jentink and Edwin Paree. They also thank Rike Wagner-Cremer, Marcio Boechat Albernaz, Christian Schwarz, Anne Baar and Jasper Leuven for helpful discussions. Diatom screening of the samples was conducted by Koeman and Bijkerk. Joep Storms and an anonymous reviewer are gratefully acknowledged for helpful feedback on the manuscript. They also thank Bram van Prooijen, Ton Hoitink, Ashish Mehta and two anonymous reviewers for comments on an earlier draft.

DATA AVAILABILITY STATEMENT
The Delft3D model software is open source and the code is available from the Deltares website (https://oss.deltares.nl/web/delft3d). Input files for the long-term reference model are included in a supporting information zip file. All field data from Rijkswaterstaat is publicly available from a variety of web portals or via the service desk (https://www.rijkswaterstaat.nl/zakelijk/opendata). The raw field core descriptions from Walsoorden are also included in the supporting information as an excel file.