Formation of the Zeebrugge coastal turbidity maximum

Formation of the Zeebrugge coastal turbidity maximum

The role of uncertainty in near-bed exchange processes

van Maren, D.S.; Vroom, J.; Fettweis, M.; Vanlede, J.

DOI

10.1016/j.margeo.2020.106186

Publication date

2020

Document Version

Final published version

Published in

Marine Geology

Citation (APA)

van Maren, D. S., Vroom, J., Fettweis, M., & Vanlede, J. (2020). Formation of the Zeebrugge coastal

turbidity maximum: The role of uncertainty in near-bed exchange processes. Marine Geology, 425,

[106186]. https://doi.org/10.1016/j.margeo.2020.106186

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

Contents lists available atScienceDirect

Marine Geology

journal homepage:www.elsevier.com/locate/margo

Formation of the Zeebrugge coastal turbidity maximum: The role of

uncertainty in near-bed exchange processes

D.S. van Maren

_{, J. Vroom}

_{, M. Fettweis}

_{, J. Vanlede}

c_{Royal Belgian Institute of Natural Sciences, Operational Directorate Natural Environment, Vautierstraat 29, 1000 Brussels, Belgium} d_{Department of Mobility and Public Works, Flanders Hydraulics Research, Berchemlei 115, 2140 Antwerp, Belgium}

A R T I C L E I N F O Keywords: Uncertainty Sand-mud interaction Turbidity Maximum Fine sediments Zeebrugge A B S T R A C T

Despite availability of a large amount of observational data and modelling studies, the mechanisms maintaining the Turbidity Maximum in the Belgian-Dutch coastal zone around the port of Zeebrugge (Belgium) are in-sufficiently understood. In order to better understand the dynamics of this turbidity maximum we examine the role of baroclinic (salinity and sediment-induced) processes and local sediment sources on the formation and persistence of the turbidity maximum through two different numerical model approaches. One model approach allows erosion of the highly compacted muddy seabed, serving as a sediment source, in line with observations of bed level change over several decades. The other approach reduces the exchange between the bed and the water column, to mimic the formation of highly concentrated near-bed suspensions with concentrations of several g/l observed around the port of Zeebrugge. Both model approaches are calibrated to various sources of available data (in situ sediment concentration observations, satellite image, bed level changes, mud content and dredging data), which they reproduce comparably well. However, reducing the water-bed exchange strengthens sediment convergence in the turbidity maximum, whereas the sediment source leads to sediment export. With the available data, it is difficult to determine which of the approaches is more realistic. Apparently, the lack of knowledge on near-bed exchange processes introduces an important source of uncertainty which cannot be adequately addressed with currently available observations. This work therefore shows that more quantitative knowledge on water-bed exchange processes in turbid marine environments is needed. It is further hypothesized that the large-scale erosion of the muddy seabed following the extension the port of Zeebrugge in the early 1980's brought such a large amount of sediment in suspension (50–100 million ton) that sediment convergence was strengthened. This increasing sediment convergence introduces a positive feedback mechanism that maintains sediment in the Turbidity Maximum, or even strengthens it. The high sediment concentrations ob-served today may therefore be a long-term effect of port construction carried out decades earlier.

1. Introduction

Estuarine Turbidity Maxima (ETM's) are regions of elevated sus-pended sediment concentration within an estuary, see e.g.de Nijs and Pietrzak (2012),Ralston et al. (2012),McSweeney et al. (2016),Grasso et al. (2018),Burchard et al. (2018), andHesse et al. (2019)for recent examples and detailed references. ETM's are the result of converging sediment pathways generated by estuarine circulation and lag effects (Dyer, 1994). Estuarine circulation is the combined effect of gravita-tional circulation (Postma, 1967), internal tidal asymmetry (Jay and Musiak, 1994) due to tidal straining (Simpson et al., 1990), lateral tidal residual flows (Lerczak and Geyer, 2004) and river flow; the relative

importance of each component is highly variable. Residual transport by time lag effects, such as settling lag and scour lag, are the result of the sediment properties (settling velocity, critical shear stress for erosion) in combination with asymmetries in the hydrodynamics (Allen et al., 1980) and topographical effects, such as divergence or convergence of channel cross sections (Friedrichs et al., 1998). Most ETM's are set within their estuary, although some are pushed into the coastal zone by high river discharge, such as the Turbidity Maxima (TM's)s of the Amazon (Kineke et al., 1996) and Yangtze river (Beardsley et al., 1985). Coastal areas with elevated suspended sediment concentration such as the Wadden Sea have been explained by the occurrence of estuarine circulation (Burchard et al., 2008). However, not all coastal TM's show

https://doi.org/10.1016/j.margeo.2020.106186

Received 17 December 2019; Received in revised form 13 March 2020; Accepted 20 March 2020

⁎_{Corresponding author at: Deltares, Dept. of Coastal and Marine Systems, P.O. Box 177, 2600, MH, Delft, the Netherlands.} E-mail address:bas.vanmaren@deltares.nl(D.S. van Maren).

Available online 08 April 2020

0025-3227/ © 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/).

a clear correlation with estuarine trapping processes, one of them being the Zeebrugge turbidity maximum. The Zeebrugge TM is located on the Belgian-Dutch continental shelf, with the Scheldt estuary as the primary fresh water source. The river discharge is low (20–600 m3_{/s with an}

average of 100 m3_{/s) compared to the tidal discharge (~50.000 m}3_/s),

and therefore insufficient to create a seaward-located ETM (as in the above examples). We therefore refer to a Coastal Turbidity Maximum (TM) rather than an ETM. However, despite generally well mixed conditions, the proximity of a fresh water source may generate salinity-driven residual currents which may play a role in the formation of the Zeebrugge TM, similar as for ETM's.

First studies on the Zeebrugge TM date back to the end of the 19th century when the planning works for the port of Zeebrugge started; they indicated a trend to deposition of muddy sediment in the near-shore area south of Zeebrugge (Van Mierlo, 1899). Later,Nihoul (1975) postulated that the Zeebrugge TM was generated by a large-scale gyre trapping fine-grained sediments transported north-eastward by the re-sidual currents. Fettweis and van den Eynde (2003) argued that the northward decrease in fine sediment transport capacity was the main mechanism, resulting in sediment deposition and subsequent strength-ening of the bed through consolidation processes during neap tides. They also hypothesize that sediment not only originates from the Strait of Dover (Irion and Zollmer, 1999) but also from local marine clay deposits deposited during the Holocene transgression. The erosion of these clay deposits (hereafter referred to as consolidated Holocene mud) accelerated in recent times as a result of maritime access works of the port of Zeebrugge (Fettweis et al., 2009a), providing a source for the TM. Bathymetric changes derived from charts indeed indicate yearly erosion up to 2.4 million tonnes of the consolidated Holocene mud bed (Bastin, 1974). Quantitative clay mineral composition of po-tential source areas showed that the paleo-estuary was the main pro-venance of the clay deposits and of the fine-grained sediments in the TM, rather than the present-day Scheldt river and estuary (Adriaens et al., 2018). Additionally, high SSC gradients (> 3 g/l difference be-tween observations 0.2 and 2.2 m above the bed) observed near-bed around Zeebrugge (Fettweis et al., 2010) may lead to sediment-induced suppression of turbulence (Winterwerp, 2001) which may generate and enhance high concentration near-bed suspensions and fluid mud layers and may also contribute to siltation in the port of Zeebrugge (Winterwerp, 2006). Clear evidence for turbulence suppression of such near-bed suspensions is provided in Fettweis's data by the observations that the concentration in the upper sensor (2.2 m above the bed (mab)) is minimal during a storm, when concentrations in the lowest sensor (0.2 mab) reach peak values. Fluid mud exists in the approach channel to the port and in the port itself, but outside these deeper channels high near-bed concentrations occur irregularly and depend on hydro-dynamic conditions. The location of the high concentration near-bed suspensions therefore also varies. And finally, disposal of fine-grained sediment from maintenance dredging works in the port of Zeebrugge and the navigation channels may strengthen the Zeebrugge TM and act as a semi-permanent source (Fettweis et al., 2011, 2016).

Mechanisms that have been demonstrated to contribute to the for-mation of Estuarine Turbidity Maxima (density-driven effects of both salinity and sediments) have to our knowledge not yet been in-vestigated for the Zeebrugge TM. Important research questions are therefore to what extent the Zeebrugge TM is formed and maintained by (1) erosion of the consolidated Holocene mud bed, (2) salinity-in-duced density effects, and (3) sediment-insalinity-in-duced density effects. In order to answer these questions a process-based numerical 3D model is set up in which sediment- and salinity-induced density effects as well as a consolidated Holocene mud source are implemented to study the effect of individual processes. Since sediment transport mechanisms in this area (tide- and wave induced resuspension and salinity-driven flows in combination with highly concentrated near-bed suspensions, in com-bination with anthropogenic dredging and dumping operations) are extremely complex, reproducing the observed sediment concentration

observations in detail is beyond the present state-of-the-art of numerical modelling. Yet, our model reproduces the most essential transport mechanisms and can therefore be applied to execute numerical ex-periments to better understand the dynamics of the TM. Of particular importance in the numerical model setup are the exchange processes between the water and the bed. Near the bed, the suspension increases in concentration from a dilute suspension to a highly concentrated benthic suspension, and within the bed from a poorly consolidated mud at the surface to -over time - a well consolidated soil at greater depth. The resulting strong variability in sediment concentration and asso-ciated transport processes requires a vertical resolution that cannot yet be resolved in numerical 3D models due to computational limitation (especially when modelling dynamic morphologic equilibrium, which requires simulation periods of several years). In order to advance our knowledge on the formation mechanisms of the TM, despite the low vertical model resolution and simple parameterizations, we apply two alternative concepts of water-bed exchange with the numerical model: (1) sediment may be supplied by erosion of consolidated Holocene mud deposits, and (2) the high near-bed sediment concentration may influ-ence sediment dynamics of the TM. Despite inevitable shortcomings in the model set-up, the results provide valuable insight into the role of model assumptions on the dynamics of turbidity maxima, and on the mechanisms responsible for the formation and maintenance of the turbidity maximum near Zeebrugge in particular.

2. Study area

The port of Zeebrugge is situated on the Belgian coastline (Fig. 1), just south of the Dutch-Belgian border and the Western Scheldt estuary (the downstream part of the Scheldt estuary). The Scheldt River is a fresh water source with an average discharge around 100 m3_{/s, ranging}

from 20 to 600 m3_{/s (Fettweis et al., 1998). The salinity in the Belgian}

coastal zone are also influenced by the larger Rhine and Meuse rivers entering the North Sea some 100 km northeast of Zeebrugge (Lacroix et al., 2004). The Zeebrugge TM is located on the Belgian-Dutch con-tinental shelf in a long band of elevated sediment concentration origi-nating from the English Channel on the French coast and extending onto the Dutch continental shelf. Residual currents along the French, Belgian, and Dutch coastlines transport sediment from the Straits of Dover northward resulting in an alongshore north-east-directed sedi-ment flux of about 30 million tonnes/year (Fettweis et al., 2007). The yearly averaged sea surface sediment concentration in the TM is around 50 mg/l (Fig. 1). Near the bed, the sediment concentration becomes several 100 mg/l during peak tidal flow conditions, and may even be several g/l during storm conditions (Fettweis et al., 2010). The highest observed suspended sediment concentrations occur near Zeebrugge (regularly exceeding 3 g/l;Fettweis et al., 2010), although high near-bed sediment concentrations (up to 2 g/l) have been episodically ob-served along the Dutch coast as well (van der Hout et al., 2017) during and after storm conditions. Trapping of sediments close to the Dutch coastline under the influence of density effects and cross-shore currents has been studied in greater detail (de Boer et al., 2009;Van der Hout et al., 2015). The Dutch coastline more towards the north is strongly influenced by the fresh water discharge from the Rhine and Meuse Rivers, generating strong salinity-induced currents. However, the Zeebrugge TM occurs in a zone where the salinity is higher (28 to 34 psu) and vertical and horizontal salinity gradients are smaller. Salinity-driven currents will therefore be less important than in the coastal area along the Dutch coast where salinity-induced trapping of sediment was extensively investigated.

The seabed of the Belgian-Dutch nearshore is composed of fine sands with variable mud content. The recent sediments in the seabed around Zeebrugge predominantly consists of sand-dominated cohesive sediment (following the classification ofvan Ledden et al., 2004), with pockets of clay-dominated cohesive sediments (Fettweis et al., 2009b). The cohesive sediments occur mainly in the Zeebrugge TM area and are

characterised by consolidated to more recently deposited layers from various ages (Fettweis et al., 2009a;Adriaens et al., 2018). Note that the stiff Holocene clays underneath the recent deposits have a much higher clay content. This consolidated mud is outcropping in-between the main navigation channel and the port of Zeebrugge (Fig. 2) and covered with recently deposited mud with critical erosion strengths of 0.5 (top layer) to 2.5 Pa (Fettweis et al., 2010). These recent deposits become several dm to > 1–2 m thick in the channels and the port, re-quiring regular maintenance dredging works. Approximately 10 million tonnes of dry sediments are annually dredged from the port of Zeeb-rugge and its approach channels (De Maerschalck and Vanlede, 2013) and disposed northwest and northeast of Zeebrugge (seeFig. 2).

The location of the mud deposits is very comparable to the location of the TM, suggesting that sediments in the TM may have been eroded from the subsoil. This relation is also suggested by the clay mineralogy of the consolidated Holocene mud layers and in the suspended matter in the TM (Adriaens et al., 2018). An analysis of frequent bathymetrical soundings data from the period 1975–2011 suggests that on average 2.8 million m3_{is eroded annually outside the maintained fairways (Fig. 3a),}

in the area with outcropping consolidated Holocene mud (compare Fig. 2withFig. 3b). Note that the spatial extent of the consolidated Holocene mud is much larger than the area for which frequent bathy-metric surveys are available and therefore the total amount of eroded sediment may be larger. A major enlargement of the port took place between 1979 and 1986, with the construction of two breakwaters extending 4 km into the North Sea. Flow contraction led to erosion of on average 3–6 million m3_{/y throughout the 1980's (Fig. 3a) and}

ap-proximately 50 million m3_{over the period 1975–2011. The bulk}

den-sity of the consolidated Holocene mud is 1500–1800 kg/m3_(Fettweis

et al., 2006), equivalent to 800–1300 kg/m3_{dry weight density.}

As-suming a dry weight density of 1000 kg/m3_{implies that the eroded dry}

mass directly resulting from the port expansion is approximately 50 million tonnes over 37 years. Erosion in front of the jetties gradually declined significantly and reached an equilibrium depth by 2005.

Fig. 1. Near surface sediment concentrations in the English Channel and southern North Sea, revealing elevated SSC near Zeebrugge computed from satellite images

covering the period 2003–2011 and using the algorithms ofNechad et al. (2010).

Fig. 2. Dredging locations (red), disposal locations

(yellow), and the area with consolidated Holocene mud in the bed (blue; based onFettweis et al., 2006). The main navigation channels connecting the North Sea with the Port of Zeebrugge and the Western Scheldt estuary cross the consolidated Holocene mud deposits. (For in-terpretation of the references to colour in this figure le-gend, the reader is referred to the web version of this article.)

3. Model setup 3.1. General model layout

The model is setup in Delft3D, solving the unsteady shallow water equations in three dimensions under the hydrostatic pressure assump-tion. Sediment transport and morphological updates are computed si-multaneously with the flow (seeLesser et al. (2004)for a description and validation). The generation, transport, and dissipation of turbu-lence are resolved with a k-ε model, in which turbulent mixing is modified by sediment-induced buoyancy effects through the equation of state (seeWinterwerp and van Kessel, 2003).

The model domain covers an 80 km long stretch along the Belgian-Dutch coast (extending 30 km in the offshore direction) and a part of the Scheldt estuary, seeFig. 4. The mesh resolution is approximately 150 × 150 m near the port of Zeebrugge and increases to 1500 × 1500 m in the North Sea. The model domain consists of 19,000 horizontal cells and has 10 vertical σ-layers with a thickness logarith-mically increasing from the bed (2% of the water depth) to the surface (20% of the water depth). The bathymetry of the model domain has

been derived from bathymetric data collected by the Dutch Ministry of Public Works and completed with lower resolution EMODnet data in the North Sea. The bed roughness is defined with a Chézy formulation using a Manning roughness which value has been obtained by cali-bration of the hydrodynamics. The model has been set up for 2014 because of availability of observations. The hydrodynamic model is coupled online with a SWAN wave model (Booij et al., 1999), using the same computational grid and bathymetry. Timeseries of wave height, period and direction measured at a wave buoy close to the model boundary are prescribed as boundary conditions, whereas wind-gen-erated waves are computed with wind observed at a coastal meteor-ological station. All boundary conditions and initial conditions are summarized inTable 1, and for conciseness not repeated in the text.

A multi-fraction morphodynamic sediment transport model has been setup with 1 sand fraction and 3 mud fractions, in which the bed level is updated at each hydrodynamic time step. Multiple mud frac-tions are needed to represent the vertical variation of the sediment concentration typically observed in nature, with the faster settling particles leading to large intratidal variability near the bed, and the slower settling particles to more constant sediment concentrations

near-Fig. 3. Annual volume change dV (panel (a), positive values indicate erosion) based on bed level observations over the Belgian shoreface covering the consolidated

Holocene mud patch (Fig. 2) but excluding the fairways. The blue line provides volumetric changes in the area of the consolidated Holocene mud patch, the red line the volumetric changes in front of the port of Zeebrugge (black polygon in panel b). Harbour extension works were carried out during the period indicated with the grey box (1979–1986). The dashed lines indicate the average volumetric changes in the period from 1975 to 2011 (2.8 and 1.6 million m3_{/y for the Holocene mud}

area and the Zeebrugge area, respectively). The right panel (b) depicts the bathymetric change (in m, red is erosion) between 1976 and 2014, with names of the two main fairways (Scheur and Zeebrugge approach channel). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Bathymetry in the numerical domain, and

location of observation stations. Black triangles are water level observations (red when combined with flow velocity observations; blank when combined with wave observations). Red circles are stations with only flow velocity observations; blank circles only wave observations; orange squares only SSC observations; the orange star refers to a station where water levels, waves, and SSC is collected. The water level and flow velocity station names are used for the aggregate bias and RMSE values and are not discussed individually. Station I is MOW1, II is Blankenberge, and III is VvdR. The wave station letters correspond to panels in Fig. S1. (For inter-pretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

surface. Sand is transported following the van Rijn (2007a, 2007b) formulations, prescribed as particles with a diameter of 200 μm and an initial layer thickness of 30 cm. Sand transport is reduced with a factor (1-p2), where p2is the fractional mass of mud in the buffer layer (see

hereafter). Transport of fine sediments (mud) is modelled with the al-gorithms developed byVan Kessel et al. (2011a), which have also been applied previously to the North Sea (van Kessel et al., 2011a), the Western Scheldt estuary (Van Kessel et al., 2011b), and the Ems Estuary (van Maren et al., 2015a). The model distinguishes two bed layers: a typically thin, fluffy upper layer of mud only (S1) which rapidly

accu-mulates and erodes, and a deeper layer (S2), which typically represents

a sandy bed in which fine sediments accumulate during calm conditions and supplies sediment during high energy events. The consolidated Holocene mud deposits (when used) and recent mud deposits are also part of this buffer layer – this is explained in the following section. The computed thickness of the fluff layer depends on erosion parameter settings and local hydrodynamics (also see hereafter).

Dredging and disposal is implemented following actual dredging strategies: sediment that deposits in navigation channels and the port of Zeebrugge is dredged when the bed level exceeds a certain intervention level and disposed on the allocated disposal areas (see Fig. 2 for dredging and disposal locations). The amount of dredged sediment can therefore be used as a validation parameter.

3.2. Erosion of mud

The erosion rate E1of layer S1linearly depends on the amount of

sediment in this layer below a user-defined threshold M0/ M1and is

independent of the amount of sediment above this threshold:

= < E mM m M M 1 cr 1 1 ,1 0 1 (1a) = > E M m M M 1 cr 1 0 ,1 0 1 (1b)

Here m (kg/m2_{) is the mass of sediment per unit surface in layer S} 1,

M0is the standard zero-order erosion parameter (kg/m2/s) and M1(1/

s) is the erosion parameter for limited sediment availability. Reducing the erosion rate of fine sediment when the availability is limited, has the benefit of a smoother and more realistic model behaviour in mixed sand-mud environments (m < M0/M1). For completely muddy areas

(m > M0/M1), the buffer model switches to the standard

Krone-Partheniades formulations for erosion of bed layer S1.

The erosion E2of S2scales linearly with the excess shear stress:

=

E p M 1

cr

2 2 2

,2 (2)

Here, p2is the fraction of mud in S2as computed by the model and

M2is the erosion parameter for S2(kg/m2/s). The erosion of sand and

mud is scaled proportionally with the available mass of sand and mud in S2. E2becomes zero if the mass of mud in layer S2becomes zero (i.e.

layer S2consists of pure sand).

Physically, when the mud fraction in layer S2is larger than the sand

fraction, layer S2represents recently deposited mud with critical

ero-sion strengths of 0.5–2.5 Pa (Fettweis et al., 2010). When sand is dominant, layer S2represents a sandy layer in which mud fills the

pores. Even though the initiation of motion for sand with D50 = 200 μm is slightly below 0.2 Pa (e.g. van Rijn, 2007a), mud within a sandy matrix is only released when large amounts of sand are brought in suspension. For such conditions, earlier model work suggests a value for τcr, 2 around 1 Pa (Van Kessel et al., 2011a). This

con-veniently corresponds to the critical shear strength of the recent mud deposits (see above). Sediment which does not or only marginally consolidates (represented by the fluff layer) typically has a critical shear stress for erosion, τcr, of several 0.01 to several 0.1 Pa (e.g.Widdows Table 1

Initial and boundary conditions.

Boundaries Initial conditions

Hydrodynamics Water levels and Riemann invariants (a combination of water level and currents allowing waves to propagate out of the model domain) at sea; tidal discharge in the Scheldt estuary; derived from a numerical model (Zijl et al., 2015)

– Salinity Offshore: measurements at various observations stations operated by the Dutch Ministry of Public

Works, close to the model boundaries. Salinity in the Scheldt boundary is from an existing Scheldt model.

30 ppt with 3 year spin-up Waves Observed significant wave heights and peak periods at a station close to the offshore model boundary,

assuming a JONSWAP spectrum (Hasselman et al., 1985). The prescribed wind field is based on observations at station Vlissingen (coastal station 30 km NE of Zeebrugge)

–

Sand Local equilibrium 30 cm with 3 year spin-up

Mud 5 mg/l offshore to 30 mg/l (summer) to-60 mg/l (winter near the coast, based on satellite images (Fettweis et al., 2007). In the Scheldt estuary 45 mg/l is prescribed. All concentrations are equally distributed over the three fractions.

0 mg/l in suspension and on the bed (both the fluff and buffer layer); 3 year spin-up

Table 2

Parameter settings of the sediment transport model.

Parameter Description Fraction

1 2 3

ws,x Settling velocity sediment fraction x [mm/s] 0.5 2 4

τcr, 1 Critical bed shear stress layer S1[Pa] 0.2 0.2 0.2

τcr, 2 Critical bed shear stress layer S2[Pa] 1 1 1

M0,x Zero-order erosion rate layer S1for fraction x [kg/m2_{/s] 10}−3 ₁₀−3 ₁₀−3

M1,x First-order erosion rate layer S1for fraction x [/s] 10−4 ₁₀−4 ₁₀−4

M2,x Erosion rate layer S2for fraction x [kg/m2_{/s] 10}−3 _2*10−3 _4*10−3

M2,x,source Erosion rate layer S2for fraction x kg/m2/s] for marine sediment source 5*10−5 5*10−4 1*10−3 General

D50 Median grain size sand [] 200 μm (ws= 22 mm/s)

et al., 2007). Therefore the critical shear stress for the fluffy layer was set low (τcr, 1=0.2 Pa for all fractions – see Table 2). Except for the

model with a consolidated Holocene mud source, all erosion parameters are prescribed as uniform values. The erosion parameters M0, M1, and M2 (Table 2) are obtained through calibration.

3.3. Settling and deposition

The deposition flux D is the product of the near-bed settling velocity ws, band the near-bed sediment concentration Cb and is divided

be-tween layers S1and S2with a burial parameter α:

=

D1 (1 ) w Cs b b, (3b)

=

D2 w Cs b b, (3c)

The value for α is based on calibration, and is typically between 0.05 and 0.2 (Van Kessel and van Maren, 2013). A low value for α implies a slow exchange with buffer layer S2. In combination with settings for M2

and τcr, 2it determines the residence time of fines in the buffer layer.

The dimensionless β is the ‘reduced deposition’ factor, introduced by Van Kessel and Vanlede (2009)to approximate physical processes that occur near the bed, but are not part of the model setup or formulations. Without reduced deposition (i.e. β = 1) the deposition rate into the bed may be overestimated because of 4 mechanisms:

1. The vertical discretization with 10 logarithmic layers is insufficient to properly reproduce high near-bed SSC. This has two implications. First, turbulent mixing is modified by sediment-induced buoyancy effects through the equation of state. An underestimation of con-centration gradients (due to low vertical resolution) leads to over-estimation of mixing (thereby further reducing near-bed con-centrations). Secondly, the settling velocity scales non-linearly with the sediment concentration. Hindered settling is modelled through an adaptation of theRichardson and Zaki (1954)formula usingws, b= ws, 0(1 − c/cg)n. Herein is ws, 0the clear water settling velocity,

cg a reference concentration (for cohesive sediments equal to the

gelling concentration, in our case set to 100 g/l) and n a power depending on the particle Reynolds number (for low particle Rey-nolds numbers n = 5). Underestimation of near-bed SSC due to the vertical discretization will therefore lead to overestimated settling velocities.

2. The critical shear stress for erosion τcr is not attained

in-stantaneously, but after deposition a particle only gradually gains strength through consolidation. A large amount of the particles that deposit on the bed, are therefore immediately re-entrained (even by fairly weak currents). A model with β = 1 will therefore always overestimate deposition rates during a phase of net deposition, especially since values for τcrare often based on the erosion phase

and therefore fairly high. As far as we know, there is no field ob-servational data on deposition / erosion processes on such small time and spatial scales.

3. Flocs may be broken up in the near-bed boundary layer, resulting in lower settling velocities (e.g.Hill et al., 2001).

4. Irregularities such as bed forms on a spatial scale smaller than the grid cell size may result in areas with low deposition rates (e.g. top of ripples) further reducing deposition. On the other hand, mud may at the same time more easily deposit on the lee side of ripples – the net effect is therefore not known.

With the current state of knowledge, it is not clear how especially mechanism 2 and 3 operate under dynamic field conditions, and con-sequently which of the 4 mechanisms is more important. However, their combined effect is a reduced sediment flux into the bed, leading to higher near-bed sediment concentrations. In this study we primarily focus on the effect of the resulting high near-bed sediment concentra-tions on ETM dynamics, but do not quantify the individual contribution

of the various mechanisms to β. We therefore evaluate the effect of highly concentrated near bed layers and TM dynamics by comparing model results without reduced deposition (β = 1, default value) to si-mulations with pronounced reduced deposition (β = 0.1). A more elaborate discussion on β follows in the discussion.

The sediment settling velocity wsis computed from current velocity

and SSC observations, measured from 2005 to 2009 at stations Blankenberge and MOW1 (seeFig. 4for location) at two positions in the water column (0.2 and 2.2 m above the bed). Under the assumption that the sediment concentration profile is in equilibrium during peak flow conditions and that the current velocity profile is logarithmic, the set-tling velocity wscan be computed from observations of SSC and current

velocities using the Rouse profile: = c c a h a h z z a P (4) In which c is the sediment concentration (g/l), cais the sediment

concentration at z = a, h is the local water depth (m), z is the vertical position in the water column (M), =P w

u

T s_{is the Rouse number with σ} T

being the Prandtl-Schmidt number (equal to 0.7), κ the von Karman constant (equal to 0.41), and u∗the bed shear velocity (m/s). Under the

assumption that a = 0.2 m (and therefore cais equal to the near-bed

measurements), P can be fitted to the data by rewriting Eq.4to = c c a h a h z z P log log a (5)

Both u∗and z0are calculated by fitting a logarithmic profile through

the current velocity observations using =

u u ln( )z u ln( )z0 ₍₆₎

For every pair of SSC observations at two depths during peak flow velocity conditions, one value for wscan be computed. Note that P does

not include the effect of turbulence damping, because turbulence damping is negligible during peak flow conditions. The computed set-tling velocities are converted into a frequency distribution of setset-tling velocity (Fig. 5). The median settling velocity is slightly > 2 mm/s.

We realise that fitting a Rouse profile to two observation points introduces inaccuracies, that the profile might not be in equilibrium at peak flow velocities, and as a consequence the higher and lower range of the computed settling velocity have high uncertainties. Also the settling velocities seem rather large compared to quite a few modelling

10−1 100 101 102 0 10 20 30 40 50 60 70 80 90 100 w_s [mm/s] % lower Blankenberge MOW1

Fig. 5. Cumulative frequency distribution (over time) of the settling velocity,

computed from a Rouse profile fit during peak flow conditions to continuous SSC measurements at stations Blankenberge and MOW1.

studies carried out in the past (e.g.Van Kessel et al., 2011a, 2011b;van Maren et al., 2015a, 2015b;Grasso et al., 2018;Toublanc et al., 2016). However, we believe that such large settling velocities are realistic for several reasons. First, such large settling velocities are not unusual in areas with concentrations up to several g/l (resulting in fast-settling flocs, see e.g. observational data by Manning and Dyer, 2007) and proved necessary for modelling the ETM of the Weser (Hesse et al., 2019). Secondly, the settling velocity of particles of sand-mud mixtures (as is the Zeebrugge area) increases relative to pure muds (Manning et al., 2011). In their experiments (with a typical range of the settling velocity of the flocs of 0.5–4 mm/s), microflocs settle much faster than macroflocs whereas in pure muds the fragile macroflocs settle faster. This may be the result of increased bonding potential between clay particles when sand is added to a muddy matrix (Mitchener et al., 1996; Manning et al., 2010). Thirdly, the good agreement in the settling ve-locity distribution between the two stations suggests that methodology (despite its shortcomings) has predictive value. And fourthly, settling velocities of several mm/s were needed to reproduce the observed SSC in the Zeebrugge TM (see hereafter).

Based onFig. 5and model calibration results, the settling velocities of the three fine sediment classes implemented in the model have been chosen as ws= 0.5, 2 and 4 mm/s (Table 2), corresponding to ~15%,

45%, and 60% of exceedance, respectively. The slowest settling fraction represents poorly flocculated mud, fraction 2 moderately flocculated mud, and fraction 3 strongly flocculated mud.

The parameter values inTable 2 are based on an extensive cali-bration procedure which involved variation of the erosion parameters, bed roughness and settling velocity. Results of this calibration proce-dure are not presented here. The model results described hereafter as the standard model results are based on the best results obtained by varying the parameters inTable 2.

4. Model results 4.1. Hydrodynamics

The hydrodynamic model has been calibrated for the year 2014 by varying the Manning's roughness parameter n. Best agreement with model and observations was obtained using a Manning's value (throughout the model domain) of 0.022 s/m1/3_{(equivalent to a Chézy}

value of 67 and a Nikuradse roughness height of 2.4 cm at a typical water depth of 10 m). Such a value represents a moderately smooth bed, indicating that bed roughness is dominated by sand particles and small-scale bedforms (the majority of the model domain is not domi-nated by mud, this is only around Zeebrugge). A detailed evaluation of the hydrodynamic model is beyond the scope of the present paper, and therefore we only present a concise summary of model performance. We compute the model bias B and the root-mean-square-error (RMSE) for 16 water levels and 9 flow velocity observation stations, using the following definitions: = M D B (6) = = RMSE sign M M D D N (M D) n [( ) ( )] N n n 1 2 (7) with M and D as the time-averaged modelled and observed variable respectively, σMand σD the standard deviation of the modelled and

observed variable, and N the number of observations. The model bias B is a measure for the mean reproduction of a particular variable whereas the RMSE quantifies the difference between model and data variability. The bias B of the astronomical water levels is close to 0 whereas the B for water levels including meteorological effects is up to 11 cm (gen-erally positive, implying that the model slightly overestimates surges). The modelled bias for flow velocity is typically 1–2 cm/s, implying hydrodynamic residuals are well reproduced. The RMSE of flow velo-cities is typically 5–10 cm/s, meaning that the model overestimates

peak flow velocities.

The modelled wave height decreases in the shoreward direction, corresponding to the observational data (Fig. S1). With a correlation coefficient typically between 0.85 and 0.9 the agreement between ob-served and modelled wave height can be considered good. The salinity typically varies between 30 and 34 ppt, with a tidal variation in the same order of magnitude as low frequency salinity variations (Fig. S2). Both the low-frequency variation in salinity (in response to wind-driven flows and river discharge – see panel a) and tidal variability (panels b and c) are largely represented by the model. The offshore salinity (imposed on the model boundaries) is 33–34 ppt whereas the salinity ranges between 14 and 25 ppt at the Scheldt boundary location. The modelled tidal variability, with values well below those imposed on the offshore boundaries, therefore demonstrates the importance of the Scheldt River discharge on the salinity distribution around Zeebrugge. 4.2. Sediment transport

The sediment transport model has been setup and calibrated in two steps. The first step is the generation of the seabed composition, which is initially only composed of sand. Tidal currents, and the superimposed northeast-directed residual flow, bring fine-grained sediments into the model domain at concentrations determined by the boundary con-centrations (5 to 60 mg/l, three fine sediment fractions, seeTable 1). As a result, the sandy seabed fills in with fine sediments until the model achieves a dynamic equilibrium between hydrodynamic forcing, off-shore boundary conditions, sediment availability and parameter set-tings. Dynamic equilibrium is defined as a state where the long-term variation (timescale of years) in the amount of sediment in the bed and in suspension is small compared to short-term variations (intratidal and spring-neap variation), and sediment fluxes into the model are equal to sediment fluxes out of the model.

Initially, the sediment flux into the model (mostly the southern boundary) is much larger than the sediment flux out of the model (the northern boundary) as sediment deposits within the model domain. Additionally, the sediment influx at the southern boundary is larger during initialisation than in equilibrium, as outflowing tidal currents (which are an order of magnitude larger than residual currents) are low in mud content during outflow. After 2 years the flux into the model is equal to the flux out of the model (around 30 million ton/year) and the influx no longer decreases. These values are close to earlier estimates by Eisma and Irion (1988)andFettweis et al. (2007), between 1.4 and 1.7 times higher than values reported byVan Alphen (1990),Lafite et al. (1993)andVelegrakis et al. (1997), and about 1.5 times lower than those given byMcManus and Prandle (1997). After 2 years, the mud content in the bed is up to 30 kg/m3_{near the harbour, and decreasing in}

offshore direction. After three years, this amount has become quasi-stationary, and this equilibrium bed provides the initial bed sediment composition for the various model runs executed as part of this paper (seeFig. 10below, computed after 4 years).

The second step is the calibration against observations of the sus-pended sediment concentration and the mass of dredged sediments. Note that for each set of model parameters, step 1 (generation of the initial condition) is done anew. Available data are (1) satellite-based maps with distribution of near-surface SSC, (2) time series of in situ observed SSC, (3) dredged masses of both mud and sand; and (4) the amount of mud in the bed. The best agreement between data and model (using realistic model settings and boundary conditions) is referred to as the standard simulation (seeFig. 6andFig. 7), using the settings of Table 2. The model-data comparison reveals that although the com-puted spring-neap variation is in line with observations, the general level of SSC is too low. Traditional calibration methods for the sediment transport module are insufficient to raise concentration levels. In-creasing erosion rates through the parameter settings leads to higher sediment concentration over short timescales, but also depletes the bed of sediment, resulting in lower SSC over longer timescales. Reducing

the erosion rate simply leads to a reduction in SSC (over short but also longer timescales). The standard simulation is an optimum setting ob-tained through a calibration procedure not shown here: apparently essential processes are missing.

Important processes cited in literature (seeSection 2) that poten-tially influence the sediment dynamics around the study area are (1) sediment is eroded from the bed, and (2) high near-bed sediment concentrations influence turbulence damping and thereby residual transport and deposition. Both these aspects are therefore further evaluated with the model.

Highly concentrated near-bed suspensions (observed during field surveys after storm conditions, see e.g.Fettweis et al. (2010)and also Fig. 6 andFig. 7) profoundly influence hydrodynamics and sediment dynamics. The resulting high vertical concentration gradients reduce vertical mixing of sediment, further strengthening the near-bed sedi-ment concentration. Also a local sedisedi-ment source (not implesedi-mented in the standard model settings) leads to higher near-bed concentrations and therefore better agreement with data. High near-bed sediment concentrations are not properly represented in the reference model because of vertical resolution limitations and poor understanding of near-bed exchange processes related to flocculation and consolidation (as explained in detail inSection 3.3); both are represented by the bulk parameter β. Using β, the near-bed concentrations are higher in areas where sediment converges, which better corresponds to observations but also influence mixing properties within the water column.

To test the effect of a local source and the reduction of near-bed

exchange of sediment, two alternative model configurations were setup and compared against the standard simulation: one with a local sedi-ment source and one with reduced bed exchange. For model alternative 1, mud is prescribed as a net sediment source on the location indicated inFig. 2. This mud replaces the sand-mud mixture and is thicker (3 m) than can be eroded within the timeframe of this model study. As this mud is consolidated, the erodibility is lower compared to more recently deposited mud. The erodibility parameters M0, M1, and M2 of the mud patch were therefore varied (seeTable 2) until the net loss from the consolidated Holocene mud area in the model (2.6 million m3_{/y) was}

comparable to the observed loss of 2.79 million m3_{/year (Fig. 3). Note}

that better agreement between observed and modelled sediment con-centrations can be achieved by either increasing erosion rates from the consolidated Holocene mud area (by decreasing τcr, 2, or increasing the

erodibility parameters M0, M1 and M2) or increasing the local rough-ness. This has not been done, however, since it would result in un-realistically large erosion rates of the consolidated Holocene mud. For model alternative 2, the near-bed sediment concentration is enlarged through reduced deposition (as explained above), parameterized by the parameter β. We have not used combinations of both approaches in order to be able to differentiate between their individual impact.

Both adaptations to the standard simulation lead to higher sus-pended sediment concentrations throughout the water column. At 2.2 m above the bed (mab) the computed sediment concentration is comparable with observations, especially the spring-neap variation (see Fig. 6c). The intra-tidal SSC variability is reasonably reproduced for −2 0 2 [m ] (a) ζ Hs 0 0.4 0.8 1.2 u [m/s] (b) 0 0.3 0.6 0.9 SSC [g/l] (c) Standard Holocene source Reduced deposition Data

13Apr0 20Apr 27Apr 04May 11May 18May 25May

0.4 0.8 1.2 SSC [g/l] [2014] (d)

Fig. 6. Computed water level ζ and significant wave height Hs (a), Depth-averaged flow velocity u (b), observed (black) and computed (coloured) SSC at 2.2 mab (c)

and 0.2 mab (d) for the standard simulation (grey), model alternative 1 (consolidated Holocene mud source, in red) and model alternative 2 (reduced deposition, in blue) during summer conditions at station MOW1 in 2014. The grey shaded areas indicate two distinct periods (storm-dominated and tide-dominated) for which more detailed results are provided inFig. 7. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

both model adaptations during tide-dominated conditions (Fig. 7f): the concentrations are typically around 0.2 g/l during peak flow conditions and reduce to several 10's of mg/l during flow reversal. During storm conditions (Fig. 7e), the model with a consolidated Holocene mud source overestimates SSC levels, but the standard model and reduced deposition approach show reasonable agreement. Again, the con-centrations are typically highest during periods with high flow velocity. The data-model comparison suggests that the additional shear gener-ated by the wave model leads to an overestimation of the sediment erosion rates.

Closer to the bed the general levels of SSC are reasonably re-produced for both model adaptations during tide-dominated conditions, but the intra-tidal phasing is poorly represented. The observations suggest that the concentrations are highest in-between peak flow and slack water, whereas the modelled concentrations are highest during peak flow (standard approach and consolidated Holocene mud source approach) or during slack tide (reduced deposition). During storm conditions, the level of near bed sediment concentration is reasonably reproduced with the consolidated Holocene mud source, but the phasing is completely wrong (Fig. 7g). The data suggests that sediment concentration variations result from vertical mixing, with a vertical uniform concentration around 0.2 g/l during high shear. During slack tide conditions, the sediment settles, leading to highly concentrated suspensions of 1.5 g/l near the bed. This vertical mixing behaviour is conceptually best represented with reduced deposition. Using this ap-proach, the near-bed concentration is indeed lowest during periods with maximum hydrodynamic energy, in agreement with observations. This is the result of vertical mixing, which will be elaborated in more detail hereafter.

The vertical structure of the computed suspended sediment con-centration (Fig. 8) reveals that during periods of high flow velocities (shortly before high and low waters) the sediment concentration fol-lows a Rouse profile. The sediment concentration is higher for the

simulation with a consolidated Holocene mud source and lower for the standard simulation. However, the concentration profiles differ most during flow reversal. With reduced deposition, suspended sediment settles towards the lower computational layers, but remains suspended in the water column. This leads to strong near bed gradients in the concentration profiles, influencing the water density and hence turbu-lence damping. These high-concentration near-bed layers are regularly observed at the observation stations in the Zeebrugge TM, supporting the reduced deposition approach.

On a larger spatial scale, the computed spatial distribution of sur-face sediment concentration in the TM is compared to satellite images (Fig. 9) and computed bed sediment distribution to observations (Fig. 10). The general pattern of elevated SSC is in line with the satellite observations, although the peak in computed SSC is slightly north of Zeebrugge whereas the data suggests SSC is highest just south of Zeebrugge (Fig. 9g,h). The spatial distribution of surface sediment is better reproduced with the additional consolidated Holocene mud source (Fig. 9c, d), especially over the triangle-shape sand bar situated at the mouth of the Scheldt estuary. The fact that the consolidated Holocene mud source is located south of estuary mouth implies that sediment eroded from the consolidated Holocene mud is transported towards this sand bar. Reduced deposition only marginally influenced the magnitude of the surface sediment concentrations but interestingly leads to a southward shift of the TM (compareFig. 9a, b withFig. 9e, f) which is more in line with the observations. Apparently, reducing the exchange of sediment between the water and the bed generates a south-westward directed transport component and/or reduces the north-eastward transport compared to standard water bed exchange simula-tions.

The seasonal variability in sediment concentration is under-estimated by the model for all scenarios. Seasonal variations in SSC are typical for mid-latitude shelf seas such as the North Sea and are related to the seasonal patterns in wind forcing and wave heights (Howarth −2 0 2 [m] (a) ζ Hs (b) 0 0.4 0.8 1.2 u [m/s ] (c) (d) 0 0.3 0.6 0.9 SSC [g/l] (e) Standard (f) Holocene source Reduced deposition Data 00 03 06 09 12 15 18 21 00 0 0.4 0.8 1.2 SSC [g/l] (g) 11−May−2014 00 03 06 09 12 15 18 21 00 18−May−2014 (h)

Fig. 7. Details of water levels and wave

height (a-b), depth-averaged flow ve-locity (c-d) and sediment concentra-tions at station MOW1 2.2 m.a.b. (e-f) and 0.2 m.a.b. (g-h) for the grey shaded periods inFig. 7. The left panels re-present a period with wave-dominated conditions (high waves, neap tide), the right panels a period with tide-domi-nated conditions (low waves, spring tide).

et al., 1993;van der Hout et al., 2017), density effects caused by higher fresh water discharges in winter, and biological effects (e.g.Jago et al., 2007). The model accounted for the seasonal variability in wind forcing (and its effect on wave height) and discharge, but not for biological effects on the settling velocity or bed erosion. This introduces three potential reasons why the observed sediment concentration in summer deviates more from that in winter, compared to the model. First, par-ticles settle faster in summer due to the greater availability of organic material strengthening flocculation (Mietta et al., 2009), especially during the spring and summer phytoplankton blooms (Jago et al., 2007; Fettweis and Baeye, 2015). Secondly, in muddy areas affected by or-ganic exopolymers the erosion rates are lower in summer due to strengthening of the bed (Kornman and De Deckere, 1998;Paterson and Hagerthey, 2001), leading to lower summer sediment concentrations. And thirdly, biological cohesion reduces bedform size and therefore bed roughness (Malarkey et al., 2015; Parsons et al., 2016); a lower bed roughness leads to lower erosion rates.

The computed amount of sediment in the bed differs strongly among the different model simulations. Obviously, the simulation with pre-scribed bed sediment shows close agreement with observed bed sedi-ment. Without sediment prescribed to the bed, the highest mud content is observed in the channels. These deposits best agree with the recent, unconsolidated muds described by Fettweis et al. (2006). The model alternative with reduced deposition predicts the lowest amount of mud in the bed, which is an obvious consequence of reducing the sediment flux to the bed. Apparently, despite a fairly good agreement of surface concentrations (as inFig. 12, for all model realisations), the predicted amount of sediment in the bed may strongly vary.

Computed dredging volumes serve as a last model validation para-meter. The dredging quantities (Table 3) correspond most to simula-tions with a consolidated Holocene mud source. This model alternative overestimates the siltation in the navigation channels, but under-estimates siltation in the Port of Zeebrugge with a factor 2. This sug-gests that mechanisms transporting sediment through the channels into the port itself are underestimated. This may be the result of vertical processes (as indicated by the low computed near-bed sediment con-centrations at MOW1 in winter) but also by the underestimated SSC

west of the port of Zeebrugge (Fig. 9- bringing relatively low SSC water into the port during flood). Furthermore, the entrance of the port is represented by 5 horizontal grid cells, which is insufficient to accu-rately resolve horizontal circulation cells (eddies), while horizontal exchange accounts for 43% of net sediment transport into the port (Vanlede and Dujardin, 2014). The computed siltation rates are overall slightly lower using reduced deposition; lowest siltation rates are computed for the standard simulation.

Data-model comparison reveals that the large-scale sediment con-centration levels are reasonably eproduced by reducing the exchange between the water column and the bed (reduced deposition approach), or by prescribing a local sediment source (consolidated Holocene mud source approach). Both approaches also resemble the intratidal varia-tion of SSC observed several meters above the bed, but fail at capturing the intratidal variation near the bed. But even though the variability is wrong, the levels of suspended sediment concentration are captured by applying reduced deposition, supporting the reduced deposition ap-proach. The bathymetric changes and the dredging quantities are more in line with the modelling approach prescribing a consolidated Holocene mud source.

Apparently, sediments in the water column resemble observational data best using reduced deposition, whereas sediments in the bed cor-respond better with data using a prescribed mud source. Which of the two approaches better describes sediment transport processes in the Zeebrugge TM cannot be established with the presently available data and process understanding – a combination of both approaches is most likely. The model does therefore not conclusively reveal which of the two investigated mechanisms is responsible for the formation of the TM (and if both are – which one is more important). In order to further understand the role of sediment supply and bed exchange, we explore the sensitivity of both model approaches to salinity and sediment-in-duced density effects.

4.3. Model scenarios

Both salinity and sediments influence the density of water through the equation of state, generating density currents and damping

00:00 03:00 06:00 09:00 12:00 -2 0 2

Water le

v

e

l [m]

F

low veloci

ty

[m/s]

(a)

SSC [g/l]

De

pth [m

a

bove bed

]

(b)

SSC [g/l]

(c)

SSC [g/l]

(d)

Holocene sediment source Hindered deposition

Fig. 8. Computed suspended sediment

con-centration at station MOW1 at peak flood flow (high water, b); slack tide (halfway ebb; c) and peak flood flow (high water; d) – see top panel (a) for corresponding water level (with red dots indicating the timing of computed SSC profiles) and depth-averaged flow velocity (positive is south-westward ebb flow). (For interpretation of the refer-ences to colour in this figure legend, the reader is referred to the web version of this article.)

turbulence mixing (Winterwerp, 2001, 2006). The effect of the two model alternatives on transport mechanisms are further explored by individually switching off the effect of salinity (by assuming all water is saline) and sediment (by switching of the contribution of the sediment concentration on the density in the equation of state) on the hydro-dynamics. Both salinity and SSC influence the hydrodynamics in two ways: (1) horizontal gradients in the salinity or SSC generate a hor-izontal pressure gradient that drives a near-bed flow towards the area of lower pressure, and (2) vertical gradients in SSC or salinity dampen vertical mixing (as computed by the k-ε turbulence model).

The difference in the computed SSC for the various model alter-natives increases towards the bed (Fig. 8), and therefore the impact on SSC is evaluated here for the near-bed layer. The impact of salinity and

sediment-induced density-effects on SSC differs substantially for both model alternatives. For both model approaches, including salinity-driven currents leads to an increase of turbidity in the nearshore northeast of Zeebrugge and a decrease in turbidity in the estuary mouth: the TM is pushed towards the coast. The effects are much more pronounced and stretch further in westward direction, however, for simulations with reduced deposition. The sediment-induced density coupling results in both cases in a lower turbidity in the estuary mouth. With a reduced deposition flux, both sediment-induced suppression of turbulence and salinity effects lead to > 100 mg/l increase in SSC (right panels inFig. 11). The average near-bed sediment concentration is many 100 mg/l, and therefore the relative near-bed increase in SSC is only several 10's %. However, sediment introduces a positive feedback

Fig. 9. Modelled and observed spatial distribution of surface SSC in summer (from 1 April to 1 October; a, c, e, and g) and in winter (from 1 October to 1 April; b, d, f,

and h) computed from the standard simulation (a, b); with a consolidated Holocene sediment source (c, d); with reduced deposition (e,f); and from satellite images (g, h). The model has been run for the year 2014, the sediment concentration is based on observations covering the period 2003–2011 (using the algorithms byNechad et al., 2010).

mechanism, where an increase in the sediment concentration may re-duce vertical mixing, further strengthening vertical concentration gra-dients. The computed near surface sediment concentrations (not shown) are indeed lower, meaning the vertical sediment concentration gradients become larger.

The impact of salinity and sediments on hydrodynamics is illu-strated with yearly average flow velocity profiles in the Zeebrugge ac-cess channel and in the approach channel to the Western Scheldt (Scheur) (Fig. 12). Horizontal salinity gradients generate salinity-driven currents which are directed towards the freshwater source near the bed and in opposite direction near the surface. As sediment concentrations increase towards the bed, salinity gradients generate a transport com-ponent directed towards the freshwater source (the Scheldt River, but also the port of Zeebrugge where local drainage water is discharged (up to 10 m3_{/s during peak discharge conditions). The overall effect is that}

sediment remains closer attached to the coast, as illustrated inFig. 11c and d.

Sediment-induced density effects only have a limited impact on the residual current (Fig. 12). The impact of sediment-induced density ef-fects is therefore primarily a vertical redistribution of sediments, with more sediments close to the bed and less up in the water column. However, when combined with salinity-driven currents (generally di-rected landward closer to the bed), sediment-induced density effects

lead to enhanced landward transport of suspended sediment. As the vertical gradients in SSC are stronger when using reduced deposition, salinity-driven currents have a more pronounced impact on the sedi-ment redistribution.

When SSC is sufficiently high, self-organizing mechanisms appar-ently exist which maintain or strengthen the TM: the higher the sedi-ment concentration, the stronger the sedisedi-ment-induced stratification, leading to larger salinity-induced landward sediment transport. The impact hereof on the formation of the Zeebrugge TM will be discussed in more detail inSection 5.3.

5. Discussion

Both model approaches presented here (reduced deposition and a local sediment source) reproduce the basic sediment dynamics in the Zeebrugge TM: the spring-neap variation and spatial distribution of SSC as well as typical dredging volumes. A number of shortcomings remain, however, such as the underestimation of thick, high concentration suspensions near the bed and the north-eastward shift of the TM. Despite these shortcomings, our model underlines that (1) our knowl-edge on near-bed sediment dynamics has significant deficiencies, lim-iting predictive modelling of fine sediments in turbid environments, (2) multiple model approaches exist which can reproduce the dynamics of

Fig. 10. Amount of mud in the bed computed by the model for the reference model (a), model alternative 1 (mud source, b) and model alternative 2 (reduced

deposition, c), in kg/m2_{(sum of layer 1 and 2, and all mud fractions). Panel d provides the measured clay fraction (fromFettweis et al., 2009b) in %.}

Table 3

Observed and modelled dredged sediment mass (in million tonnes/year) and the mud fraction (by mass, in %) for the three main dredging areas in the TM: the port of Zeebrugge, the approach channel to Zeebrugge, and Scheur (seeFig. 3for locations). Observed dredging numbers are based on personal communications with the Flemish ministry of Public Works, compiled inVroom and Schrijvershof (2015). The model scenarios with the closest agreement are in bold.

Simulation Zeebrugge Approach channel Zeebrugge Scheur

Total (106_t/

year) Obs./comp.(%) Mud fraction(%) Total (10 6_t/

year) Obs./comp.(%) Mud fraction(−) Total (10 6_t/

year) Obs./comp.(%) Mud fraction(−)

Observed 5.96 / 77 1.06 / 66 0.68 / 43

Standard 1.17 19 64 0.68 64 16 1.32 194 6

Consolidated Holocene mud

source 2.99 50 98 1.13 107 91 1.45 213 85

the turbidity maximum to a similar level of accuracy, and the under-lying assumptions are important when interpreting modelled turbidity maximum dynamics; and (3) both approaches provide new insights into the formation and maintenance of the Zeebrugge Turbidity Maximum. These will be discussed in more detail below.

5.1. Near-bed sediment dynamics

The reduction in the deposition flux β is to some extent similar to the well-known and often-used Krone equation (Krone, 1962). The Krone equation only allows sedimentation below a critical user -defined parameter, and the closer the bed shear is to the critical bed shear for deposition, the lower the deposition flux (if we express reduced de-position as β = (1 − τb/τb, cr) for τb< τb, cr, we obtain the Krone Fig. 11. Increase in near-bed SSC by (a, b) sediment-induced density effects on turbulence suppression (simulation with sediment and salinity effects minus a

simulation without sediment effects but with salinity effects) and (c, d) salinity-induced density effects (simulation with sediment and salinity effects minus a simulation with sediment effects but without salinity effects). The left panels (a, c) depict simulations with the consolidated Holocene sediment source whereas the right panels (b, d) depict simulations with reduced deposition.

-0.02 0 0.02

u [m/s]

z/

h

[-]

(a)

u [m/s]

(b)

Fig. 12. Modelled flow velocity profile (averaged over one

year) in the Zeebrugge entrance channel (a) and in the approach channel to the Western Scheldt (Scheur; panel b). Positive values are directed towards the port of Zeebrugge in (a) and towards the Western Scheldt in (b). The thin solid lines are velocities computed with a barotropic model excluding salinity, whereas a baroclinic model including salinity was run with deposition efficiency DE (thick solid line) and without DE (dashed line).

equation). However, there are several arguments against the Krone equation. As demonstrated byWinterwerp (2007), there are no physical arguments (or proper data) to exclude mutual erosion and sedimenta-tion as in the Krone formulasedimenta-tion as a funcsedimenta-tion of shear. Also the as-sumption of constant settling velocities (as implicitly in the Krone for-mula) is not realistic for high sediment concentrations (Mehta et al., 2014).

Reduced deposition is hypothesized to be the combined effect of vertical discretization, strength development, floc dynamics in the bed boundary layer, and local irregularities (Section 3.3). In its present form β is independent of shear, sediment concentration and settling velocity. However, there are three arguments suggesting that reduced deposition should become more pronounced at higher sediment concentrations.

First, reduced deposition leads to a steeper sediment concentration profile (Fig. 8 andFig. 13b) which in turns leads to suppression of turbulent mixing at the lutocline (Winterwerp, 2001, 2006). With less turbulent mixing, sediment above the lutocline rapidly settles from the upper part of the water column, further steepening the concentration profile (Fig. 13c). Reduced deposition therefore favours the formation of highly concentrated near-bed layers. In the Rouse concentration profile, the sediment concentration gradients scale linearly with the sediment concentration SSC, and therefore the contribution of me-chanism (1) to β scales with SSC.

Secondly, the sediment concentration influences β through con-solidation. Consolidation times scale with the square thickness of the initial deposits (as follows from the diffusion term in theGibson et al. (1967) consolidation equation). At low sediment concentrations, the deposition flux is small, and therefore all sediment on the bed con-solidates instantaneously. However, sediment does not consolidate in-stantaneously when settling rates are large, leading to a longer period of low critical bed shear stresses. And since deposition scales linearly with the sediment concentration (as in Eq. (3)), the contribution of strength development by consolidation to β scales with SSC2_.

Thirdly, the sediment concentration influences deposition through hindered settling. It was argued inSection 3.3that as models with low vertical resolution underestimate the near bed SSC, they underestimate hindered settling and therefore overestimate deposition. However, an underprediction of SSC also reduces the deposition flux wsC, which

would imply β > 1. Which of these two is more important, depends on the near bed concentration and the gelling concentration cg(Fig. 14) via

de formula ofRichardson and Zaki (1954). Typically, cgvaries between

20 and 100 kg/m3_{. For c}

g= 20 kg/m3, the deposition flux wsC

de-creases with increasing SSC when SSC > 3 kg/m3_{, hence hindered}

settling is dominant, requiring β < 1; for cg= 50 kg/m3 hindered

settling dominates when SSC > 8 kg/m3_{. The higher the gelling}

con-centration, the higher the near bed sediment concentration needs to be

before hindered settling becomes dominant.

When the concentrations are in the hindered settling range, un-derprediction of the SSC (due to model resolution) leads to over-prediction of the deposition flux (therefore β < 1). At low sediment concentrations (SSC less than several kg/m3_{), the near-bed vertical}

discretization is more likely to underestimate the deposition flux (as in Spearman and Manning, 2008). This means that β should be smaller at higher sediment concentrations. In the concentration range near Zeebrugge (SSC several g/l near-bed) the effects of increased deposition and hindered settling are comparable, motivating a β ≈ 1.

The settling velocity wsinfluences β in several ways. At very low

settling velocity,β = 1 because the vertical gradients are absent whereas at very high settling velocity β = 1 because particles deposit instantaneously. Hence, an optimum in wsmust exist where it

maxi-mally influences β. The settling velocity is, in turn, related to the se-diment concentration through flocculation: typically, the floc size and therefore wsincreases with the sediment concentration. However, ws

also depends on shear.Dyer et al. (2002)andManning et al. (2007) observe that the settling velocity in a highly concentrated near-bed suspension in the Tamar Estuary increases, due to damping of turbulent

Fig. 13. Illustration of the effect of reduced deposition (RD) on a Rouse sediment concentration profile (a). Reduced deposition fluxes result in higher near-bed

sediment concentrations and a lutocline-shaped profile (b), which is enhanced by turbulence suppression (DT; c).

Fig. 14. Relation between deposition flux wsC as a function of concentration C,

for a gelling concentration cgof 20, 50 and 80 kg /m3(black lines) and at

in-crements of 5 kg/m3_{in-between 10 and 100 kg/m}3_{(grey). Ws is computed}