WL2020R13_131_6.pdf by Vlaanderen-be

13_131_6 FHR reports

Integraal plan Boven-Zeeschelde Sub report 6 – Scaldis Mud a Mud Transport model for the Scheldt Estuary

www.flandershydraulicsresearch.be

Integraal plan Boven-Zeeschelde Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Smolders, S.; Bi, Q.; Vanlede, J.; De Maerschalck, B.; Plancke, Y.; Mostaert, F.

Cover figure © The Government of Flanders, Department of Mobility and Public Works, Flanders Hydraulics Legal notice Flanders Hydraulics Research is of the opinion that the information and positions in this report are substantiated by the available data and knowledge at the time of writing. The positions taken in this report are those of Flanders Hydraulics Research and do not reflect necessarily the opinion of the Government of Flanders or any of its institutions. Flanders Hydraulics Research nor any person or company acting on behalf of Flanders Hydraulics Research is responsible for any loss or damage arising from the use of the information in this report. Copyright and citation © The Government of Flanders, Department of Mobility and Public Works, Flanders Hydraulics Research 2020 D/2020/3241/201 This publication should be cited as follows: Smolders, S.; Bi, Q.; Vanlede, J.; De Maerschalck, B.; Plancke, Y.; Mostaert, F. (2020). Integraal plan Boven-Zeeschelde: Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary. Version 4.0. FHR Reports, 13_131_6. Flanders Hydraulics Research: Antwerp. Reproduction of and reference to this publication is authorized provided the source is acknowledged correctly. Document identification Customer: Keywords (3-5): Knowledge domains

Flanders Hydraulics Research Ref.: WL2020R13_131_6 Mud, Scheldt Estuary, Modelling Hydraulics and sediment > Sediment > Cohesive sediment > Numerical modelling

Text (p.): Confidentiality:

77 ‫ ܈‬No

Author(s):

Sven Smolders; Qilong Bi

Appendices (p.): ‫ ܈‬Available online

Control Name Yves Plancke Reviser(s):

Signature Getekend door:Yves Plancke (Signature) Getekend door:Bart De Maerschalck (Sign Getekend op:2021-02-02 09:12:47 +01:00 Getekend op:2021-01-07 13:54:36 +00:00 Reden:Ik keur dit document goed Reden:Ik keur dit document goed

Bart De Maerschalck

Project leader:

Joris Vanlede

Getekend door:Joris Vanlede (Signature) Getekend op:2021-01-22 11:17:32 +01:0 Reden:Ik keur dit document goed

Approval Head of Division:

F-WL-PP10-2 Version 7 Valid as from 3/01/2017

Frank Mostaert

Getekend door:Frank Mostaert (Signature Getekend op:2021-01-07 08:51:31 +00:0 Reden:Ik keur dit document goed

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Abstract This report is number 6 in a report series for the project integrated plan Upper Sea Scheldt. This report presents the Scaldis Mud model setup and calibration results. The mud model is a sediment transport model built in SEDI3D, the 3D sediment transport module within TELEMAC-3D. The mud model is built on top of the 3D hydrodynamic model, Scaldis 2013 (Smolders et al., 2016).

Final version

WL2020R13_131_6

III

F-WL-PP10-2 Version 7 Valid as from 3/01/2017

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Contents Abstract ............................................................................................................................................................ III Contents ............................................................................................................................................................ V List of tables.................................................................................................................................................... VIII List of figures .................................................................................................................................................... IX 1

Introduction ............................................................................................................................................... 1 Integral Plan Upper Sea Scheldt ........................................................................................................ 1 Goal of the mud model...................................................................................................................... 1 B-Alternatives .................................................................................................................................... 2 1.3.1

The VaG-alternative ................................................................................................................... 2

1.3.2

The Chafing-alternative (Schaaf) ............................................................................................... 2

1.3.3

The VaH-alternative ................................................................................................................... 2

The mud model in the modelling chain ............................................................................................. 3

1.4.1

SSC in the Ecosystem Model (UA) ............................................................................................. 3

1.4.2

SSC in the fish (Twait Shad) model (INBO; Vanoverbeke et al., 2019a) .................................... 4

1.4.3

SSC in the Bird model (INBO; Vanoverbeke et al., 2019b) ........................................................ 5

1.4.4

Using Delta’s in the Modelling Chain......................................................................................... 6

1.4.5

Conclusions on handling delta’s of SSC in the modelling chain ................................................ 7

Literature Review of sediment transport modelling ................................................................................. 8 Maintenance of Deurganckdok (2004) .............................................................................................. 8 3D mud transport model Sea Scheldt (2006) .................................................................................... 9 2.2.1

Computation of bed shear stress ............................................................................................ 10

2.2.2

Continuous sedimentation ...................................................................................................... 10

2.2.3

Flocculation ............................................................................................................................. 10

2.2.4

Consolidation ........................................................................................................................... 11

Development zone Saeftinghe (2011) ............................................................................................. 11 LTV Mud model (Delwaq) ................................................................................................................ 12 Stortstrategie Beneden Zeeschelde (2016) – mud model ............................................................... 13 3

Measurement Data on SSC ...................................................................................................................... 16 Continuous SSC data ........................................................................................................................ 16 Grain Size Distribution in the Sea Scheldt ....................................................................................... 19 Ebb tide longitudinal transects ........................................................................................................ 21

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

3.3.1

Measuring Points ..................................................................................................................... 21

3.3.2

Ebb tide longitudinal transect: campaign 2012 ....................................................................... 21

3.3.3

Ebb tide longitudinal transect: campaign 2013 ....................................................................... 23

3.3.4

Ebb tide longitudinal transect: campaign 2014 ....................................................................... 23

Longitudinal transects around high and low water slack ................................................................ 26 Longitudinal transect derived from tide independent datasets ..................................................... 26 Ebb and flood asymmetry in SSC ..................................................................................................... 28 Spring/Neap variation in SSC ........................................................................................................... 29 Dredging and Disposal ..................................................................................................................... 29 4

The Hydraulic Model: Scaldis 3D ............................................................................................................. 33 General introduction ....................................................................................................................... 33 Computational Grid and Bathymetry .............................................................................................. 33 Boundary Conditions ....................................................................................................................... 33 Salinity ............................................................................................................................................. 34

Cohesive Sediment Transport Model: Scaldis Mud ................................................................................. 35 Introduction ..................................................................................................................................... 35 Conceptual Model ........................................................................................................................... 35 5.2.1

Advection-Diffusion ................................................................................................................. 35

5.2.2

Bed Shear Stress ...................................................................................................................... 35

5.2.3

Erosion ..................................................................................................................................... 36

5.2.4

Deposition................................................................................................................................ 36

5.2.5

Bed evolution........................................................................................................................... 37

5.2.6

Bed Structure ........................................................................................................................... 37

5.2.7

Density ..................................................................................................................................... 37

Sediment properties ........................................................................................................................ 38 Boundary conditions........................................................................................................................ 39 5.4.1

Hydrodynamic boundary conditions ....................................................................................... 39

5.4.2

Sediment input from the boundaries ...................................................................................... 41

Initial conditions .............................................................................................................................. 44 Handling tidal flats ........................................................................................................................... 44 Disabling bottom update in SEDI-3D ............................................................................................... 45 Decoupling the effect of SSC on relative density ............................................................................ 46 Using a constant Manning Bottom Roughness for sediment transport.......................................... 47 Reduced settling velocity on shallow areas using a logistic function .............................................. 48 Reduced settling velocity on upstream boundaries and entrance of locks .................................... 49 Dredging and Disposal Flux ............................................................................................................. 50

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Results ..................................................................................................................................................... 51 Spin-up time .................................................................................................................................... 51 6.1.1

Set up harmonic run ................................................................................................................ 51

6.1.2

Mass balance information from the listing file ....................................................................... 52

6.1.3

Time scales to reach equilibrium: conceptual model .............................................................. 56

Calibrated model results ................................................................................................................. 57 6.2.1

Mass balance of the model ..................................................................................................... 57

6.2.2

Mass transfer maps ................................................................................................................. 59

6.2.3

Dredging mass balance ............................................................................................................ 60

6.2.4

Averaged ETM.......................................................................................................................... 61

6.2.5

Ensemble analysis SSC ............................................................................................................. 63

6.2.6

Sedimentation on intertidal areas in the Upper Sea Scheldt .................................................. 65

Simulation speed-up times .............................................................................................................. 70 7

Conclusions .............................................................................................................................................. 72

Recommendations for further improvements ........................................................................................ 73

References ............................................................................................................................................... 74

Appendix A: Software choice for mud modelling ............................................................................................ A1 Delwaq ......................................................................................................................................................... A1 Comments in the .cli file (resolved) ......................................................................................................... A1 Number of exchanges in the .hyd file (resolved) .................................................................................... A1 Aggregation for Telemac not implemented (unresolved) ....................................................................... A1 Parallelisation (unresolved) ..................................................................................................................... A2 Conclusion ............................................................................................................................................... A2 SISYPHE ........................................................................................................................................................ A2 Error advection schemes ......................................................................................................................... A2 Appendix B: treatment of tidal flats in Telemac 3D ........................................................................................ A7 Variables involved in treatment of tidal flats .............................................................................................. A7 LECDON_TELEMAC3D.f................................................................................................................................ A8 TELEMAC3D.f ............................................................................................................................................... A9 Appendix C: Formulas for dredging amounts ................................................................................................ A14 Notation ..................................................................................................................................................... A14 Volume concentration ............................................................................................................................... A14 Reduced Volume V’ ................................................................................................................................... A14 Conversion of 𝑽𝑽′𝒃𝒃 [m³] to 𝒎𝒎𝒎𝒎 [TDS]......................................................................................................... A15

Appendix D: Steering file calibrated model ................................................................................................... A16

Final version

WL2020R13_131_6

VII

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

List of tables Table 1 – Independent variables in the Fish (twait shad) model (INBO)........................................................... 5 Table 2 – Independent variables in the bird model (INBO). .............................................................................. 6 Table 3 – Cohesive sediment parametrization in IMDC (2004) ......................................................................... 8 Table 4 – Parameters LTV Mud model – changes from 2006 to 2010 ............................................................ 13 Table 5 – Stations with turbidity measurements ............................................................................................ 16 Table 6 – Total dredging mass of mud from Lower Sea Scheldt in 2007 - 2015 [ton/y]. ................................ 30 Table 7 – Dredged mud mass from different zones in the Lower Sea Scheldt [ton/y]. .................................. 31 Table 8 – Keywords TELEMAC-3D for sediment properties ............................................................................ 39 Table 9 – Sediment loads per subcatchment .................................................................................................. 41 Table 10 – Total discharge and average sediment concentration per subcatchment .................................... 42 Table 11 – Keywords for tracer handling in Telemac3D.................................................................................. 43 Table 12 – Keywords for handling tidal flats in Telemac3D ............................................................................ 45 Table 13 – Possible ways of reducing sedimentation on intertidal flats ......................................................... 68 Table 14 – list of all relevant variables for treatment of tidal flats in Telemac 3D ......................................... A7

VIII

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

List of figures Figure 1 – Modelling chain ................................................................................................................................ 3 Figure 2 – Network of the state variables in the UA water column model, involved in the key pelagic processes of primary production and nutrient cycling. .................................................................................... 4 Figure 3 – Habitat suitability for twait shad larvae in function of SSC .............................................................. 5 Figure 4 – Using Delta's in the modelling chain ................................................................................................ 7 Figure 5 – Model domain and initial mud concentration in IMDC (2004) ........................................................ 8 Figure 6 – model domain in the 2006 study ...................................................................................................... 9 Figure 7 – Observed shear stresses and settling velocities, fitted to the flocculation model of Winterwerp et al. (2006).......................................................................................................................................................... 10 Figure 8 – model domain for the study of

development zone Saeftinghe.................................................. 11

Figure 9 – Model grid of the Delft3D model (red) in the NEVLA domain (green) ........................................... 13 Figure 10 – 4 different disposal scenario’s considered in (Coen et al., 2016b)............................................... 15 Figure 11 – Map of the Scheldt estuary with indication of location names and km from estuary mouth ..... 17 Figure 12 – Ensemble Analysis at Bouy 84. Median, P25 and P75 .................................................................. 18 Figure 13 – Ensemble Analysis at Oosterweel. Median, P25 and P75 ............................................................ 18 Figure 14 – Ensemble Analysis at Driegoten. Median, P25 and P75 ............................................................... 18 Figure 15 – Grain size distributions for samples taken with Delft Bottle on Frame (DBF), Delft Bottle Suspended (DBH) and Suspension Samples (SUSP) during flood .................................................................... 19 Figure 16 – Grain size distributions for samples taken with Delft Bottle on Frame (DBF), Delft Bottle Suspended (DBH) and Suspension Samples (SUSP) during ebb ...................................................................... 20 Figure 17 – Median grain sizes LISST-100X, Delft bottle on a frame and pump samples. .............................. 21 Figure 18 – Ebb tide longitudinal transect: measured SSC (water samples), campaign of August 2012 ........ 22 Figure 19 – Ebb tide longitudinal transect: campaign of November 2012...................................................... 22 Figure 20 – Measured longitudinal profiles of SSC at 1m below the surface in 2013..................................... 23 Figure 21 – Measured longitudinal profiles of SSC at 1m above the bottom in 2014 .................................... 24 Figure 22 – Measured longitudinal profiles of SSC at 1m below the surface in 2014..................................... 24 Figure 23 – Fresh water inflow in the Upper Sea Scheldt in 2014 .................................................................. 25 Figure 24 – median surface SSC (and P25 and P75) around LWS (blue) and HWS (red) ................................. 26 Figure 25 – Median surface SSC with error bars between P25 and P75 for tidally independent dataset over period 1971-2015, clustered per 10 km .......................................................................................................... 27 Figure 26 – Median SSC over different depths in the water column with error bars between P25 and P75 from OMES dataset over period 2001-2010 ................................................................................................... 27 Figure 27 – Median Surface SSC with error bars between P25 and P75 ......................................................... 28 Figure 28 – Spring Neap variation in Surface SSC, calculated over the period 1996-2015 ............................. 29

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Figure 29 – Disposed volumes of mud converted to TDS in the Lower Sea Scheldt 2007 - 2015. .................. 30 Figure 30 – Location of most important dredging (red) and disposal (green) locations................................. 31 Figure 31 – Distribution (2007-2015) of dredged material over different dredging zones. ........................... 32 Figure 32 – Representation of the bed structure in SEDI-3D .......................................................................... 37 Figure 33 – Plots of floc diameter versus settling velocity. ............................................................................. 38 Figure 34 – The synthetic boundary for the calibrated model ........................................................................ 40 Figure 35 – Discharges QN boundary condition .............................................................................................. 40 Figure 36 – New location of Q boundary at Merelbeke .................................................................................. 41 Figure 37 – Annual mean SSC in the North sea with the location of the Scaldis model sea boundary .......... 43 Figure 38 – Manning bottom roughness coefficient of Scaldis HD 2013 along the estuary axis. ................... 47 Figure 39 – Effect of Manning coefficient on bed shear stress; varying flow velocities U and water depths. ......................................................................................................................................................................... 48 Figure 40 – sigmoid (logistic) curve alpha ....................................................................................................... 49 Figure 41 – Location of point source of sediment in the Scaldis model ......................................................... 50 Figure 42 – Model domain of Scaldis for Harmonic simulations ..................................................................... 52 Figure 43 – Mass balance plot for harmonic simulations HAR_5 and HAR_6 ................................................. 54 Figure 44 – An example of intertidal areas (red) and subtidal channels (blue) indicated by the polygons used in Vandenbruwaene (2016). ............................................................................................................................ 55 Figure 45 – Averaged bed evolution for tidal flats and inside the main channel for the Upper Sea Scheldt . 55 Figure 46 - The conceptual model showing the time scale for reaching the equilibrium between erosion and deposition. ....................................................................................................................................................... 57 Figure 47 – Mass balance from listing file for the whole model domain ........................................................ 58 Figure 48 – Mass balance of the Upper Sea Scheldt of last 20 days of simulation (calculated from masses in UA polygons) ................................................................................................................................................... 58 Figure 49 – Mud transport calculated from the model compared with calculated mud transport rate from Vandenbruwaene et al. (2017) ........................................................................................................................ 59 Figure 50 – Location of polygons (red) from where the total amount of sediment on the bed is calculated 60 Figure 51 – Estimation of sedimentation in dredging polygons for last 20 days ............................................ 61 Figure 52 – Averaged location of ETM with P25 and P75 ............................................................................... 62 Figure 53 – Location of ETM; visualization of all time steps ........................................................................... 62 Figure 54 – Effect of dumping of sediment near Oosterweel on the average SSC values along the estuary . 63 Figure 55 – Ensemble analysis Buoy 84 ........................................................................................................... 64 Figure 56 – Ensemble analysis Oosterweel ..................................................................................................... 64 Figure 57 – Ensemble analysis Driegoten ........................................................................................................ 65 Figure 58 – Bed thickness after 40 days of simulation for the region Temse and Rupelmonde as an example ......................................................................................................................................................................... 66 Figure 59 – Average bed evolution of the Upper Sea Scheldt......................................................................... 66 Figure 60 – excessive sedimentation in shallow areas in SA_033 with sedimentation rates up to 7.5 m/yr. 67

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Figure 61 – Averaged erosion/deposition flux at a transect just downstream Rupelmonde for different solutions .......................................................................................................................................................... 69 Figure 62 – Zoom on left bank of Figure 61 - Averaged erosion/deposition flux at the left bank of the transect just downstream Rupelmonde for different SA runs ........................................................................ 69 Figure 63 – Scaling performance of TELEMAC-3D with SEDI-3D on different queues at FHR ........................ 70

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

1 Introduction Integral Plan Upper Sea Scheldt The implementation of the Seine-Scheldt connection will result in increased shipping traffic between France and Flanders. Waterwegen en Zeekanaal NV (W&Z) wants to improve the navigability of the Upper Seascheldt in order to prevent a bottle neck effect on the Canal Ghent – Terneuzen and the Western Scheldt. Additionally, an improved navigability of the Upper Seascheldt would also promote the shipping between the Scheldt basin and the Albert Canal (towards the river Meuse). Within this framework, an integrated plan is being developed, in which navigability, safety and nature are the key elements. At the moment, the upstream part of the Upper Seascheldt is a Class IV fairway (ships up to 85m long and 9.5m wide) and forms a bottleneck in the European network. The questions that are to be answered within the integrated plan pertain to the measures that need to be taken to upgrade the Upper Seascheldt to a Class Va fairway suitable for ships up to 2250 tons (ships up to 110m long and 11.4m wide), taking into account the other functions (safety, nature and recreation). The outcome of a feasibility study was that with relatively small measures (see further) a balance between cost and benefit can be found, but allowing navigability up to Class Va while increasing safety for ships of class IV and lower. The integrated plan is aimed at further developing the conclusions from the feasibility study towards Class Va shipping. It is of the utmost importance that the design of this enlargement leads to a multifunctional Scheldt with assets for navigability, guarantees for protection against flooding and a sustainable natural system. The mud model presented in this report is part of a model chain to evaluate the different geometrical scenarios. Fine sediments play a key role in aquatic ecosystems like the Scheldt estuary. They determine light penetration into the water column and hence affect the primary production. They determine the layers of the bed supporting benthic life and the sediment’s organic content forms food supply to filter feeders. Therefore the behavior of these fine/cohesive sediments is important in the assessment of the impact of changes in bathymetry or management of the estuary and for this project, the Upper Sea Scheldt in particular.

Goal of the mud model The goal of the mud model within this project is to estimate the effect of different possible measures in the Upper Sea Scheldt on cohesive sediment concentrations (mud) in the estuary. Modeled effects on SSC provide input for subsequent modelling done by project partners: • • •

Ecological Modelling [UA] Fish model [INBO] Water Bird model [INBO]

This modelling chain and the underlying assumptions are described in §1.3. Besides the specific goals for this project, the mud model has some general goals of a mud model for the Scheldt estuary to fulfill in order to achieve a certain quality and trustworthiness.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

B-Alternatives The B-alternatives have been defined in a preceding study (Bevaarbaarheid van de Boven-Zeeschelde en Zuidelijk vak Ringvaart voor klasse Va-schepen). The B-alternatives consist of three different potential designs. For a more detailed description of the alternatives, the reader is referred to the scenario analysis report of the hydrodynamic model (Hassan et al., 2019). 1.3.1 The VaG-alternative In the VaG-alternative, different bends are cut off and the navigation channel is also widened at a number of locations: the reclamation at Wijmeers, GOG-GGG Bergenmeersen (Uitbergen), GOG Scheldebroek, The Castle of the Dike and the Pottelbergse Schorren (Grembergen), the Groot Schoor near Vlassenbroek, the Vlassenbroekse Schorren, Uiterdijk, the Cramp and Sint-Amands. Additional modifications to the bathymetry of the Upper Sea Scheldt and the Ringvaart result in additional expansion of the navigation channel. 1.3.2 The Chafing-alternative (Schaaf) In the Chafing-alternative the bathymetry of different curves/bends is adjusted in the Sea Scheldt. For this alternative, the land boundary remains unchanged. 1.3.3 The VaH-alternative The VaH-alternative combines the above alternatives with interventions by VaG up Wichelen and changes in the bathymetry of the chafing alternative down Wichelen. This hybrid alternative makes maximum use of the existing channel.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

The mud model in the modelling chain Within the project, modelled effects of scenario’s on SSC will be used by other project partners in a modelling chain. This chain is illustrated below and further described in IMDC (2014) Figure 1 – Modelling chain

1.4.1 SSC in the Ecosystem Model (UA) SSC concentrations in the ecosystem model are important for the estimation of the primary production. SSC determines how much the photosynthetically active radiation can reach into the water column for primary production. In the ecosystem model light dependence is included using a modified Platt model (Platt et al., 1980). The depth dependence of the light climate is assumed to follow an exponential decrease (Lambert-Beer law) with light extinction coefficient Kd. The calculation of Kd from SPM concentrations is based on data. There are different fitted relations for different estuarine zones (Upper Seascheldt, Lower Seascheldt, Western Scheldt and Rupel). Primary production is integrated over photic depth (irradiance >= 1% of surface irradiance) by means of Gauss quadrature. This is followed by integration over the surface area of the specific pelagic box divided by its volume, which gives the concentration change in state variables related to primary production (e.g. oxygen, diatom C, nitrate) .”SSC is a forcing of the ecosystem model and is not a state variable in this model. Figure 2 shows the relationships of the state variables in the water column in the ecosystem model. For more detailed information about the ecosystem model the reader is referred to van Engeland et al. (2018).

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 2 – Network of the state variables in the UA water column model, involved in the key pelagic processes of primary production and nutrient cycling.

The ecosystem model of the current state of the estuary is fed by monthly averaged SSC data (interpolated daily) (and 5-day averaged discharges). These SSC values are derived from the subset of surface concentrations from the OMES data. It is the surface SSC that is limiting the light climate and thus the primary production. Therefore, the delta’s that are calculated from the Scaldis mud model for this application are derived from the calculated surface SSC. 1.4.2 SSC in the fish (Twait Shad) model (INBO; Vanoverbeke et al., 2019a) The fish model determines the habitat suitability for the twait shad and is based on the model of Stevens et al. (2011). The model contains two components: the habitat suitability for spawning and the habitat suitability for nursing the fish in their larvae state. Full grown twait shad migrates to open sea. According to Stevens et al. (2011) turbidity, which is directly related to SSC, has only a direct effect on the larvae nursery habitat suitability. Indirectly the SSC could have an effect on primary production in the ecosystem model and so on minimal oxygen concentrations. All the independent variables in the twait shad model are listed in Table 1. The habitat suitability for the twait shad larvae in function of the SSC is given in Figure 3. This figure was made based on SSC tolerance limits found in literature and is thus independent of input data. Therefore, SSC is thus not a calibration parameter for the twait shad model. The twait shad model is not calibrated. It is only validated based on recent data from 2007-2012.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Table 1 – Independent variables in the Fish (twait shad) model (INBO). Direct dependence on SSC in dark green, indirect dependence in lighter green

Parameter Water quality Average temperature

Unit

Model

Temporal resolution

Spatial resolution

°C

100 m 100 m 100 m

Average salinity (Chloride) Minimum Oxygen

‰

spawn larvae spawn larvae

mg/l

spawn larvae

Turbidity (SSC)

mg/l

larvae

May May - June May May - June May May - June May - June

mg C/l

larvae

May – June

100 m

year round

100 m

spring

100 m

Biotic variables Zooplankton

100 m 100 m

Habitat quality water depth (maximum depth at low water) average flow velocity

m/s

spawn & larvae spawn

maximal flow velocity

m/s

larvae

spring

100 m

km/day

spawn

spring

100 m

Net water displacement

Figure 3 – Habitat suitability for twait shad larvae in function of SSC

1.4.3 SSC in the Bird model (INBO; Vanoverbeke et al., 2019b) The bird model is a Generalized Linear Mixed Model (GLMM). This means that it is assumed that there is a linear relationship between the explanatory variables in the model and the number of birds. Table 2 gives an overview of the explanatory variables in two bird models: teal (wintertaling) and shelduck (bergeend). SSC is directly linked to the number of teal birds. It has no direct effect on the shelducks. SSC contributes indirectly to both bird models by its influence on oxygen. It also indirectly affects the phytoplankton, which is an direct explanatory variable in the teal bird model.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Table 2 – Independent variables in the bird model (INBO). Direct dependence on SSC in dark green, indirect dependence in lighter green

Teal (wintertaling) Shelduck (bergeend) Habitat Surface area mud flat Spread in exposure time Interaction surface area x spread Water Quality Oxygen Phytoplankton SSC

X X X

All the explanatory variables are scaled with an average of 0 and a standard deviation of 1. This how the variables are used in the two bird models. SSC is thus rescaled to a SSCbeta : 𝑆𝑆𝑆𝑆𝑆𝑆𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 =

𝑆𝑆𝑆𝑆𝑆𝑆 − 𝑎𝑎𝑎𝑎𝑎𝑎 𝑠𝑠𝑠𝑠𝑠𝑠

with avg the average SSC value and std the standard deviation on the SSC values. 1.4.4 Using Delta’s in the Modelling Chain In this project, we use a data-assimilation approach to predict SSC under different scenario’s (reference year 2050) under three new bathymetries (see §1.3). The data-assimilation is illustrated in Figure 4. The measurements of SSC serve a double purpose. First of all they are used to calibrate the mud model. The mud model is then used to calculate the difference between the reference situation and a scenario. This difference is expressed as a ratio (termed delta, or relative difference). Delta’s or relative differences in SSC between model runs are calculated as: ∆= �〈

𝑆𝑆𝑆𝑆𝑆𝑆𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 − 𝑆𝑆𝑆𝑆𝑆𝑆𝑟𝑟𝑟𝑟𝑟𝑟 〉� 𝑆𝑆𝑆𝑆𝑆𝑆𝑟𝑟𝑟𝑟𝑟𝑟

With { } an averaging operator over space and 〈 〉 an averaging operator over time. Delta is dimensionless, but has a sign. The definition of the spatial and temporal averaging is done together with the project partners. Delta’s are calculated over the spatial boxes of the ecosystem model, and over the duration of the simulation after spinup. This is further elaborated in the scenario analysis report of the mud model (Bi et al., in prep). If further down the modelling chain, the expected SSC is needed in mg/l, they can be calculated using the delta’s and the measured value. This is illustrated in Figure 4. The underlying hypothesis is that the model is better capable of predicting the sign and amplitude of effects, rather than predicting the absolute values of SSC.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 4 – Using Delta's in the modelling chain

1.4.5 Conclusions on handling delta’s of SSC in the modelling chain The effects that need to be modelled in this project are twofold: • •

What is the effect of a change in boundary conditions from 2013 (calibration) to 2050 (reference for scenario analysis)? What is the effect of changing the geometry of the system in the Upper Sea Scheldt in the reference year 2050?

Effects are calculated from calculated effects on surface SSC (top interface in the 3D Scaldis Mud model), as all models further in the modelling chain are calibrated based on the surface SSC (schepstalen) in the OMES dataset.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

2 Literature Review of sediment transport modelling In this section a short review is given of other models handling cohesive sediments.

Maintenance of Deurganckdok (2004) This study estimated the maintenance dredging volumes in the (then still to-be built) Deurganckdok. The study was meant to optimize the dredging strategy of Deurganckdok. A 3D Hydrodynamic and sediment transport model was built in Delft3D. The model boundaries are located at Waarde and Schelle. The model counts 10 sigma layers, with a higher vertical resolution close to the bottom. The model grid was based on the 2D Scalwest model. The grid resolution in the zone of interest is ~40m. The hydrodynamic boundary conditions were derived from 2D Scalwest. Mud availability in the bottom is initialized based on the maps of Wartel, Parker et al. (1999) Figure 5 – Model domain and initial mud concentration in IMDC (2004)

Cohesive sediment is modelled using a standard Krone-Partheniades with one fraction. The following parametrization is used: Table 3 – Cohesive sediment parametrization in IMDC (2004)

Parameter ws τcrit, depo τcrit, ero M Ρdry 8

Value 2 mm/s 0.2 Pa 0.4 Pa 5*10-5 kg/m²/s 244 kg/m²

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

The model is calibrated on average ebb and flood concentrations on 4 locations: Meestoof, Liefkenshoek, Drempel van Zandvliet and Bath, based on 13h measurements performed 1990 and 1997. The spring-neap variation is checked against the analysis of the long-term measurements at Prosperpolder. Model runs typically last 1 tidal cycle (although the report isn’t particularly clear on the exact run duration and how spin-up was handled). Based on the model, the sedimentation of DGD was estimated at 1975 tons of dry matter per tide, or 1.4MT/yr.

3D mud transport model Sea Scheldt (2006) The 3D mud transport model was developed to simulate the mud transport in the Lower Sea Scheldt, with special attention paid to exchange processes between a harbor basin like Deurganckdok and the estuary. To this purpose, the functionalities of Delft3D were extended for this study. The new developments include improved formulations for bed friction, erosion, flocculation, consolidation and a low Reynolds turbulence model. In order to accurately simulate consolidation, a general sigma vertical grid system was also implemented. The details are described in WL | Delft Hydraulics (2007). The hydrodynamic model was calibrated against measurements of the exchange flow at the Kallo lock. The sediment transport model was calibrated against sediment concentration measurements carried out in the main channel of the Scheldt close to Deurganckdok during the HCBS field campaign in 2005. The model is also validated against sedimentation patterns observed in Deurganckdok in 2005-2006. The model grid is based on the based on the existing 2004 model described above and in more detail in IMDC (2003, 2004). The boundaries are kept at Waarde and Schelle. The grid resolution is refined locally using domain decomposition. The grid resolution around DGD was refined down to ~20m. The grid resolution around Kallo lock is ~10m. In the vertical, 40 sigma layers are distributed logarithmically. Boundary conditions were calculated using the NEVLA model of the entire Scheldt estuary. The HLES model (Uittenbogaard,1998) is used to compute the horizontal component of the sub-grid eddy viscosity and the sub-grid eddy diffusivity. Figure 6 – model domain in the 2006 study

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

The fluxes between the water phase and the bed are calculated with the Partheniades-Krone formulations (Partheniades, 1965). The calibration was started using the mud parametrization used in IMDC (2004). For the initial mud concentration on the bottom, the same values are used as in IMDC (2004). The technical report notes high sediment concentrations in the lower depth cells of drying/flooding cells (up to 2g/l). The calculation remains stable however, and a sensitivity analysis on the drying/flooding parametrization in Delft3D gave no substantial differences. 2.2.1 Computation of bed shear stress Previous formulations to compute the bed shear stress were based on the depth-averaged flow velocity and user-specified Chézy coefficient. This Chézy coefficient incorporates both the form drag of bed forms and the friction of individual grains, and used to calibrate the hydrodynamics. However, the form drag does not physically contribute to the erosion of mud beds. Therefore the skin friction formulations by Soulsby and Clarke (2005) are used when the thickness of the mud bed exceeds a critical value. This critical thickness is set to 0.01 m. This process is parametrized with a Nikuradse roughness of 0.001m. 2.2.2 Continuous sedimentation Continuous sedimentation is assumed, and the erosion rate is based on basic soil properties that can be measured in a laboratory. This method is described in detail in Winterwerp et al. (2006). Note that continuous deposition instead of deposition during slack tide only results in a gross tide-average sedimentation flux that is in the order of 10 times higher. To maintain a similar balance between sedimentation and erosion (i.e. to maintain a similar equilibrium concentration), also the tide-average erosion flux should become in the order of 10 times higher. The critical erosion threshold and the erosion constant can be related to physical parameters following Winterwerp and van Kesteren (2004). 2.2.3 Flocculation The flocculation model of Winterwerp et al. (2006) was calibrated against field data of the HCBS campaign. This model describes flocculation as the result of turbulence- induced aggregation and floc break-up. The latter two processes work continuously, and at equilibrium, they balance. The model has one characteristic floc size, which can be regarded as the median floc size. Figure 7 – Observed shear stresses and settling velocities, fitted to the flocculation model of Winterwerp et al. (2006)

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Flocs that have formed in the Scheldt river are transported into the low energy docks, where they rapidly settle from suspension. However, the default flocculation model would predict low sediment settling rates because the degree of turbulence is low. Therefore an alternative flocculation model setting is used in the final calibrated model were the fall velocity does not decrease as much at low turbulence levels. 2.2.4 Consolidation The original aim of the new mud model was to additionally include consolidation. The numerical difficulties associated with the consolidation module forced the project to be carried out without the consolidation module.

Development zone Saeftinghe (2011) van Maren et al. (2011) describe the calibration of the sediment transport model that was used in the design of a new tidal dock in the Scheldt (Saeftinghedock). The numerical model is based on the model that was developed earlier for the research of Deurganckdok, and that is described in the previous paragraph. Figure 8 – model domain for the study of development zone Saeftinghe

The model grid has a horizontal resolution of ~50m around Saeftinghedock and Deurganckdok. In the vertical, 20 sigma layers are distributed logarithmically. Boundary conditions are imposed at Schelle and Waarde. The current deflecting wall at Deurganckdock is implemented. The model is recalibrated for hydrodynamics (van Maren et al., 2010) and mud transport (van Maren et al., 2011)

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

LTV Mud model (Delwaq) In the LTV mud model the bed of the estuary is represented by two layers. The first layer is a thin fluffy fine sediment layer deposited during slack water. At high current velocity, most or nearly all of this layer is resuspended into the water column. The critical shear stress for resuspension τcrit1 of this layer is low and its erosion constant M1 is high. If less than a certain mass m1 per unit area of fine sediment is available in layer 1, it may well be assumed that the surface coverage of the underlying bed forms is not complete. In this case, the resuspension constant M will become dependent on the percentage of surface coverage. A transition between zeroth order and first order resuspension behavior occurs. The expression for the erosion flux from layer 1, Fero1, now reads: Fero1 = min(M0,mM1) x max(0,(τ/τcrit1 – 1)) where M0 and M1 are the zeroth and first order resuspension constants, respectively and m the available sediment mass per unit area in layer 1. By definition, M0 = m10 M1. The second layer represents the sand bed, which prevails in the Scheldt estuary. There is an erosion flux of fine sediment, present in the pores of the sand bed: D* = D50 ((s-1)g/ν2)1/3 Fero2 = pM2ρs((s-1)gD50)0.5 D*0.3 (τ/τcrit2 – 1)1.5 = pM’2(τ/τcrit2 – 1)1.5 where M’2 = 333 M2 for the LTV mud model. The value 333 follows from the standard values of ρs = 2600 kg/m³, s = ρs/ρw = 2.5, D50 = 3*10-4 m. The power 1.5 in the expression is inherited from the Van Rijn type of erosion function used. The erosion rate increases linearly with the mud fraction p. The sedimentation flux is split into two fractions and parameter α steers the sedimentation towards layers 1 and 2: Fsed1 = (1-α) ws C Fsed2 = α ws C As α << 1, the rate of sediment exchange between the water column and the first layer is much higher than the rate of exchange with the second layer. In combination with a much higher typical sediment mass per unit area in layer 2 compared with layer 1, the residence and response times of sediment in layer 2 are much longer. Whereas layer 1 responds on the tidal time scale (hours), layer 2 responds on the seasonal scale (months to years). The neap spring tidal cycle (14 days) may influence both layers. There are two mud fractions present in the model: a fluvial and a marine fraction. No sand transport is modelled. The bed level remains unchanged in this model. No morphological changes are modelled. The model does not take the following physics into account: fluid mud, flocculation, biological influence, consolidation, sand/mud behavior and sediment-water interaction. (Van Kessel et al., 2006). In 2007 the LTV mud model was further improved by prolonging the hydrodynamic simulation period to one year (previous spring neap tidal cycle). A shift in dumping locations by a few 100 meters towards the navigation channel was performed. New boundary conditions were given to the coastal boundary. The number of layers of the model was doubled from 5 to 10, but gave no significant improvements. In the 2006 version, wave effects were included by applying a constant, year-averaged wind speed of 7 m/s and assuming a constant fetch length of 25 km. Subsequently, the wave-induced bed shear stress was computed from the equilibrium wave height as a function of the local water depth throughout the model domain. In this 2007 version a variable wind forcing on waves was calculated with a Swan model and the simulated wave maps were exported to Matlab to be converted to wind induced bed shear stress maps and are as such applied in the model. The effects of biology on stabilizing or destabilizing sediment were investigated, but were not taken into account in the final version. The influence of shipping traffic was investigated by increasing the bed shear stress locally by 0.5 Pa and adding a probability of p = 1%. This had significant effects, but needs to be implemented in much bigger detail to represent real shipping traffic. (Van Kessel et al., 2007) 12

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

In 2009 the model was further improved by introducing a deposition probability mimicking the effect of hindered settling. This should increase near-bed concentrations. The hypothesis is that for the original settings for settling velocity and deposition, near-bed concentration peaks during slack do not form sufficiently. The explanation is that for a vertically uniform concentration distribution, the deposition flux towards the bed Fdep equals the sedimentation flux Fsed through the water column near the bed: Fsed = Fdep = ws C. By introducing a deposition probability p<1, the near-bed concentration peaks do form during slack, as Fdep < Fsed and the excess sedimentation flux accumulates near the bed. As a second improvement the excessive sediment deposition in the upper reaches of the tributaries was suppressed (Rupel basin) (Van Kessel and Vanlede, 2009). Table 4 shows the parameter value evolution of the LTV mud model from 2006 to 2010. Table 4 – Parameters LTV Mud model – changes from 2006 to 2010

parameter ws α

τcrit1 M1 d τcrit2 M2 Psedmin

fraction 1

1 mm/s 0.1 0.2 Pa 2.3*10-5 s-1 0.05 m 0.5 Pa 3.5*10-7 kg/m²/s -

2006

fraction 2

1 mm/s 0.05 0.2 Pa 2.3*10-5 s-1 0.05 m 0.5 Pa 3.5*10-7 kg/m²/s -

fraction 1

2007

1 mm/s 0.1 0.1 Pa 2.3*10-5 s-1 0.05 m 0.5 Pa 3.5*10-7 kg/m²/s -

fraction 2

1 mm/s 0.05 0.1 Pa 2.3*10-5 s-1 0.05 m 0.5 Pa 3.5*10-7 kg/m²/s -

fraction 1

2010

1 mm/s 0.1 0.1 Pa 1.16*10-5 s-1 0.05 m 0.5 Pa 1.75*10-7 kg/m²/s 0.1

fraction 2

1 mm/s 0.1 0.1 Pa 1.16*10-5 s-1 0.05 m 0.5 Pa 1.75*10-7 kg/m²/s 0.1

Stortstrategie Beneden Zeeschelde (2016) – mud model This study was performed by Flanders Hydraulics in collaboration with Antea, UA and INBO within the framework of the “Agenda voor de Toekomst” research agenda. Within the study possible disposal strategies in the Sea Scheldt were investigated (for both sand and mud). A 2D model was used in this study that is a cut-out and 3x3 refinement of the NEVLA grid. The detailed model is implemented in Delft3D ‘SEDonline’. The downstream boundary is located downstream at Ossenisse and upstream at Tielrode and Terhagen. Figure 9 – Model grid of the Delft3D model (red) in the NEVLA domain (green)

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Coen et al. (2016a) describe the set-up of the mud model. They chose to disable all interactions from the mud transport model back to the hydrodynamics. Mud transport is modelled in two fractions. Bed roughness is modelled with a uniform Manning value of 0.025. The model is calibrated against continuous SPM data at Boei 84, Oosterweel and Driegoten. Parameter

Mud Fraction 1

Mud fraction 2

Dry Density [kg/m³]

550

Initial mud concentration [mg/l]

150

0.2

τcr, ero [Pa]

0.4

0.8

τcr, depo [Pa]

1000

Settling velocity [mm/s]

The initial sediment layer thickness at bed is input as a spatially varying field, derived from lithological information combined with the thickness of the erodible layer

In the scenario calculations (Coen et al., 2016b), mud disposal is implemented in the model as a solid discharge of 1212 TDS per disposal over 6 minutes. Four different disposal scenario’s (reference and 3 variations) are tested in the model

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 10 – 4 different disposal scenario’s considered in (Coen et al., 2016b). Reference scenario is top left

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

3 Measurement Data on SSC Continuous SSC data There are three locations in the Sea Scheldt where turbidity was measured continuously for a long time period: Bouy84, Oosterweel and Driegoten. The turbidity measurements are transformed into SSC values using a relationship between turbidity and SSC derived from measured water samples taken at these locations. The locations are listed in Table 5 with their characteristics and are shown on a map in Figure 11. For the three location data was available for 2013 and 2014 with a 10 minute interval (data provide by HIC Vlaanderen). Table 5 – Stations with turbidity measurements

Location

Height

Buoy 84 top

3.3 m above the bottom

Buoy 84 bottom

0.8 m above the bottom

Oosterweel top

4.5 m above the bottom

Oosterweel bottom

1 m above the bottom

Driegoten

3 m below the water surface

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 11 – Map of the Scheldt estuary with indication of location names and km from estuary mouth

For these three locations an ensemble analysis is performed on the 2 years of available SSC data. The time series of SSC concentrations are split up into individual tidal cycles, and grouped in 3 tidal amplitude classes with reference to tidal amplitude at Antwerp. The highest one third of tidal cycles is labeled as spring tide. The middle third is labeled as average tide and the lowest third of amplitudes are labeled as neap tide. Note that this division does not correspond to the astronomical definition of spring and neap tide (which has 2 spring and neap tidal cycles per lunar month), but still provides a meaningful subdivision along tidal amplitude, with sufficient tidal cycles in each tidal class for statistical ensemble analysis. The results are shown below for the three stations (Figure 12 - Figure 14). In the ensemble analysis the tidal signal is divided in individual tides and every tide is subdivided in time relative to the high water level. For a specific time to high water the average parameter value (in this case SSC) is shown along with its confidence interval.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 12 – Ensemble Analysis at Bouy 84. Median, P25 and P75

Figure 13 – Ensemble Analysis at Oosterweel. Median, P25 and P75

Figure 14 – Ensemble Analysis at Driegoten. Median, P25 and P75

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Grain Size Distribution in the Sea Scheldt Plancke et al. (2014b) reports measurements of sediment transport performed using a Delft Bottle (both on frame and suspended) executed specifically for this project. Because the Delft Bottle is designed to capture only sand and allows finer particles (< 63 µm) to flow through the bottle, those measurements will be used later on in the project to calibrate the sand transport. Grain sizes were determined on the sediment trapped in the Delft Bottle, and on sediment in suspension that was pumped up. Grain size is determined at the Sediment Laboratorium of Flanders Hydraulics Research (FHR) using the Mastersizer 2000, which uses laser diffraction. Samples are sieved at 2000 µm before analysis. Organic matter was not removed. As expected, grain sizes obtained from the Delft Bottle are more course than grain sizes observed in suspension, because the Delft Bottle acts as a filter that captures the sandy fraction (Figure 16 and Figure 16). The names used in Figure 15 and Figure 16 can be traced back on the map in Figure 11. If we focus on the characteristic grain size for the sediment in suspension, we have a d50 in the range of 15-25 µm. Note that this is a measure of the size of primary particles. Figure 15 – Grain size distributions for samples taken with Delft Bottle on Frame (DBF), Delft Bottle Suspended (DBH) and Suspension Samples (SUSP) during flood (Plancke et al., 2014b)

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 16 – Grain size distributions for samples taken with Delft Bottle on Frame (DBF), Delft Bottle Suspended (DBH) and Suspension Samples (SUSP) during ebb (Plancke et al., 2014b)

Thant et al. (2016) report on a measurement campaign at Ketelplaat where grain (floc) sizes where determined both from a LISST-100X, samples taken with a Delft bottle and pump samples. The Delft bottle was mounted on a frame, at 35 cm and 47 cm above the bottom. Sampling time was 3–5 minutes, after which the sediment load was defined. Figure 17 shows the grain sizes determined with different techniques. The higher median grain sizes of the LISST-100X can originate from the initial presence of sediment flocs in the system. These flocs brake up by pumping up the water, during transport or by handling of the water samples in the lab, leading to a higher concentration of fine grained particles and lower median grain sizes. Note that the sampling technique of the Delft bottle implies a selection of particle grain size (only sand is trapped in the bottle, fines (<63μm) will pass through). Median grain sizes of pump samples at Ketelplaat are considerably lower compared to the d50 values of the LISST-100X. Lab conditioned measurements by the LISST-100X on calibrated sand (105μm) also showed an overestimation of the expected median grain size (130μm–150μm). It is therefore not clear how much of the deviation in d50 between LISST-100X and water samples is due to measurement configuration or effective change in sediment properties. Over- and underestimations of median grain sizes by the LISST are also encountered in other studies.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 17 – Median grain sizes LISST-100X, Delft bottle on a frame and pump samples. US stands for ultrasonic treatment of the sample

Ebb tide longitudinal transects 3.3.1 Measuring Points In order to capture the longitudinal profile of SSC in the upper Sea Scheldt, profiles were sailed and samples were taken every 2 km. Samples were taken near the surface (approximately 1m below surface) and near the bottom (approximately 1m above the bottom). These measurements campaigns were already performed in the 1970’s with sailed transects at slack high and low water. Since the 2000’s these measurements were restarted, and since 2009, extra campaigns were sailed at half-tide ebb. Initially these transects covered only the Lower Sea Scheldt, from 2012 they were extended towards the Upper Sea Scheldt. These measurements were performed during ebb tidal conditions, because of the relatively low variation of velocities during this tidal phase. This means that the entire transect was measured within 3 days. One day covers the transect from Melle to Dendermonde another day covers the transect from Dendermonde to Rupelmonde (these locations can be found in Figure 11), and the third day covers the zone between the Belgian-Dutch border and Rupelmonde. The results for the year 2012, 2013 and 2014 are discussed in the following sections. 3.3.2 Ebb tide longitudinal transect: campaign 2012 The 2012 campaign is discussed in detail in Plancke et al. (2012). In 2012 there was a summer and autumn measurement campaign. In the summer there was more variation between SSC values at surface and near the bottom (Figure 18). During the autumn campaign the bottom and surface SSC values lie closer together (Figure 19), indicating a better mixed water column during ebb.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 18 – Ebb tide longitudinal transect: measured SSC (water samples), campaign of August 2012

Figure 19 – Ebb tide longitudinal transect: campaign of November 2012

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

3.3.3 Ebb tide longitudinal transect: campaign 2013 The campaign of 2013 is discussed in detail in Plancke et al. (2014a). Figure 20 shows the seasonal variation of surface (1 m below the surface) SSC in the Sea Scheldt. The measurement in march is considered to be representative for winter conditions, because it was preceded with a period characterized by higher fresh water inflow in the Upper Sea Scheldt. The measurement in September is considered to be representative for a summer period because it was preceded by low fresh water inflow. A clear difference between summer and winter conditions can be seen. In winter, higher sediment concentrations are found in the Lower Sea Scheldt, with a maximum of 0.5 g/l around Antwerp (km 88). The upper Sea Scheldt, i.e. upstream of km 125, is characterized by low values of less than 0.1 g/l, with value of 0.05 g/l up-esuary of KM 140. In summer, the highest sediment concentration is found around Temse (km 113) and Driegoten (km 120). Also upstream the values are higher with concentrations around 0.15 g/l. Figure 20 – Measured longitudinal profiles of SSC at 1m below the surface in 2013

3.3.4 Ebb tide longitudinal transect: campaign 2014 The campaign in 2014 is summarized here and for more in depth information we refer to Plancke et al. (2015). Figure 21 (bottom) and Figure 22 (surface) show the seasonal variation of SSC in the Sea Scheldt. The measurements in February are considered to be representative for winter conditions, as it is preceded by a period of higher fresh water inflow (Figure 23). Since there was a peak discharge in August, the measurements in September are not considered to be representative for summer conditions. The measurement in November (autumn in Figure 13 and Figure 22) is believed to be more representative of summer conditions (concerning the fresh water inflow regime; Figure 23). There is a clear difference in the longitudinal profile of SSC between summer and winter conditions. In winter, the highest sediment concentrations are in the Lower Sea Scheldt, with a maximum of 0.4 g/l (surface) or 0.6 g/l (bottom) around Antwerp. The upper Sea Scheldt upstream of Dendermonde (km 137) is Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

characterized by values lower than 0.1 g/l (both surface and bottom). In autumn (which is considered a representative for summer because of the fresh water inflow conditions), the highest sediment concentrations (around 0.4 g/l surface, peaks to 1 g/l bottom) are found around Driegoten (km 120). In the Lower Sea Scheldt, lower values around 0.15 g/l are found. Figure 21 – Measured longitudinal profiles of SSC at 1m above the bottom in 2014

Figure 22 – Measured longitudinal profiles of SSC at 1m below the surface in 2014

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 23 – Fresh water inflow in the Upper Sea Scheldt in 2014 (red lines indicate the measurement campaigns)

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Longitudinal transects around high and low water slack Vandenbruwaene et al. (2016) analyzed measurements that were gathered around high water slack (HWS) and low water slack (LWS). The longitudinal stretch of the data goes from Buoy 79 downstream (around Bath km 63) to Rupelmonde (km 106). The data cover the period of 2001-2003, 2005, 2008 and 2010-2012, for LWS and 2002-2003 and 2008-2012 for HWS. Both measured values for high and low water slack are in the same range. Only at both the downstream end and upstream end of the measured stretch, the SSC values of low water slack are higher than those at high water slack. Figure 24 – median surface SSC (and P25 and P75) around LWS (blue) and HWS (red) (Vandenbruwaene et al., 2016)

Longitudinal transect derived from tide independent datasets Vandenbruwaene et al. (2016) analyzed the available data in the DONAR, CEME and OMES datasets (all three containing tidally independent SSC measurements). Figure 25 shows median surface SSC derived from the combination of these 3 datasets. It shows surface SSC around 0.025 g/l in the Western Scheldt, rising to 0.075 g/l around Rupelmonde (km 106). The highest median surface SSC values (0.1 g/l) were measured between Temse (km 113) and Dendermonde (km 137). It is noted that the measurements during ebb tide (section 3.3) in the Upper Sea Scheldt are substantially higher than the range presented in Figure 25. No clear direct reason can be given. The OMES dataset also contains data of SSC measured in the water column. Vandenbruwaene et al. (2016) clustered that data in four depth classes (expressed in % of the water column/ water depth) (Figure 26). In general there is an increase of SSC with depth, as expected.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 25 – Median surface SSC with error bars between P25 and P75 for tidally independent dataset over period 1971-2015, clustered per 10 km

Figure 26 – Median SSC over different depths in the water column with error bars between P25 and P75 from OMES dataset over period 2001-2010

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Ebb and flood asymmetry in SSC Vandenbruwaene et al. (2016) calculated the tidal phase of every measurement in the DONAR, CEME and OMES database. Even though the measurements were performed tidally independently, the intertidal variation could still be quantified. Figure 27 shows concentrations in the Sea Scheldt (upstream from km 60) that are slightly higher during ebb than during flood phase. However notice that the variability is clearly bigger than the average difference. The difference is smaller and less uniform in the Western Scheldt. Figure 27 – Median Surface SSC with error bars between P25 and P75

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Spring/Neap variation in SSC Vandenbruwaene et al. (2016) also determined the tidal amplitude for each measurement in the DONAR, CEME and OMES database to calculate the spring-neap variation in SSC. In general, the SSC during neap tide is 0.8 times the SSC during average tide conditions. For spring tide conditions, the ratio lies around 1.2 (Figure 28). Similar ratios were found for water column SSC in the Sea Scheldt. Figure 28 – Spring Neap variation in Surface SSC, calculated over the period 1996-2015

Dredging and Disposal Santermans (2013) gives an overview of the yearly dredged and disposed sediment volumes in the Lower Sea Scheldt from 1981 - 2011. Santermans (2013) expresses the reported volumes in reduced volumes (V’). A reduced volume means that the original wet volume of mud is converted to the volume that the mud would have if its density was 2000 kg/m³. Because we use mass concentrations in the model, the reported reduced volumes were converted to mass using the following formula: 𝑀𝑀𝑑𝑑 = 1.606 ∙ 𝑉𝑉 ′

with 𝑀𝑀𝑑𝑑 the dry matter mass [tons] of the dredged volume. The derivation of this formula and the underlying assumptions are given in appendix C. The converted volumes to mass units from Santermans (2013) are presented in Figure 29. This figure shows also the distribution of the disposed sediment over the disposal sites in the Sea Scheldt. The most important sites (in magnitude of dredged and disposed sediment) are shown in Figure 30. Until 1999 most of the mud was deposited on the Plaat van Boomke. From 2000 on Punt van Melsele and Oosterweel are the most important deposition sites. The three

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

disposal zones are located in the same vicinity. Oosterweel is located on the right bank (Figure 30) and is used for disposal during ebb. Plaat van Boomke (Figure 30) is located on the right bank and is used for disposal during flood. Punt van Melsele (Figure 30) is located on the left bank and is used for disposal during flood. Table 6 gives the dredged amounts of sediments for each year in tons/year. Table 7 gives an average dredged mud mass over the period 2007-2015 for specific locations in the lower Sea Scheldt. Figure 31 shows the distribution of the total dredged volume of sediment over the different sites for the years 2007-2015. Figure 29 – Disposed volumes of mud converted to TDS in the Lower Sea Scheldt 2007 - 2015. Presented values are derived from reduced volumes of mud (m³ V’) using a multiplication factor of 1606 kg/m³ V’. Data source: Maritime Access division.

Table 6 – Total dredging mass of mud from Lower Sea Scheldt in 2007 - 2015 [ton/y]. Presented values are derived from reduced volumes of mud (m³ V’) using a multiplication factor of 1.606 ton/m³V’. The dredging in harbors is not included in this table. Data source: Maritime Access division.

year

Total mass dredged from lower Sea Scheldt [ton/year]

2007 2008 2009 2010 2011 2012 2013 2014 2015 Average (2007-2015)

2,879,608 3,062,342 3,587,897 4,602,115 7,670,508 5,513,671 4,989,945 4,769,268 3,838,599 4,545,995

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 30 – Location of most important dredging (red) and disposal (green) locations

Table 7 – Dredged mud mass from different zones in the Lower Sea Scheldt [ton/y]. Presented values are derived from reduced volumes of mud (m³ V’) using a multiplication factor of 1.606 ton/m³V’. Data source: Maritime Access division

Dredging location

Average (2007-2015)

Deurganckdok

1,473,357

Drempel van Frederik

1,169,005

Toegang Zandvliet- en Berendrechtsluis

824,186

Kallosluis

519,208

Drempel van Lillo

192,443

Drempel van Zandvliet

155,515

Toegang Boudewijn- en Van Cauwelaertsluis

135,441

Rest

76,841

Haven Linkeroever

30,631

Noordzeeterminal

26,654

Ketelplaat

18,253

Drempel van Krankeloon

1,304

Total 4,545,995

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 31 – Distribution (2007-2015) of dredged material over different dredging zones. Data source: Maritime Access division

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

4 The Hydraulic Model: Scaldis 3D General introduction The Scaldis model was developed within the framework of the project 13_131 Integrated Plan for Upper Sea Scheldt. There was a need for a numerical model with higher mesh resolution in the upstream part of the Scheldt estuary. An unstructured mesh has the advantage to be very flexible in changing mesh resolution in specified areas of interest. Therefore the Scaldis model was developed in the TELEMAC software, which is based on a finite element method.

Computational Grid and Bathymetry The computational grid includes the Belgian coastal zone, extended to Dunkirk in France and Westenschouwen in the Netherlands. The grid includes the Eastern Scheldt as well. The mesh resolution varies between 400 and 150 m in this part of the model. In the Western Scheldt and estuary mouth area the grid resolution increases towards 120 to 80 m. The grid includes all tributaries reaching as far as the tidal influence. The grid resolution keeps increasing all the way till the upstream boundary in Merelbeke where it is 5-7 m. All flood control areas of the Sigma Plan are included in the grid. The 2D grid consists of 459.692 nodes. The 3D mesh consists of prisms eventually cut into tetrahedrons and is automatically constructed from the 2D mesh. A sigma transformation is used for the vertical location of these 5 layers. Layer 1 is the bottom layer and the following layers are always situated on 0.12, 0.30, 0.60 and 1.00 fraction of the water depth. The bathymetry is interpolated from measured data from 2013 or the closest date available. For more detailed information we refer to the calibration report of the Scaldis model from Smolders et al. (2016).

Boundary Conditions There are 9 liquid boundaries in the Scaldis model. The downstream boundary is located in the North Sea. A water level is imposed on this boundary. This boundary contains 469 nodes. On every node a water level time series with 10 minute interval is imposed. These time series are extracted from the regional ZUNO model of the southern North Sea. A correction of the harmonic components was calculated based on the comparison of the harmonic components of the ZUNO results and measurements over a period of 1 year (Maximova et al., 2015). There are 8 upstream liquid boundaries with prescribed discharges. Measured daily average discharges are available as upstream boundary conditions for Merelbeke (ADCP Melle), Dender, Zenne, Dijle, Kleine Nete, Grote Nete, and channel Ghent – Terneuzen. The values for the channel Ghent-Terneuzen were taken downstream the weir of Evergem. This location was used as a proxy. For the channel of Bath hourly discharge measurements were available. Wind influence is not included in the hydrodynamic model. It is assumed to be incorporated in the boundary conditions downstream.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Salinity Salinity is applied as an active tracer. This means that density effects are on in the model. Time series with 10 minute interval are imposed on the downstream boundary. The salinity values are generated from the CSM-ZUNO model train. Salinity boundary values in the Scaldis model were corrected based on the comparison of the simulated and measured salinity time series at Vlakte van de Raan (located in the larger mouth area of the Scheldt Estuary). No tracer (Salinity) values are set for the upstream boundary conditions. No salinity is entering the domain through these boundaries in the Scaldis HD model. The model starts from an initial salinity field: a map is made based on a combination of salinity measurements (Western Scheldt) and corrected model results from ZUNO (coastal area). This initial salinity map is read at the start of a simulation by a modified subroutine FONSTR.f. The values of the 2D map are copied to the other four layers in the model.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

5 Cohesive Sediment Transport Model: Scaldis Mud Introduction SEDI-3D is a suite of subroutines inside the TELEMAC-3D module of the TELEMAC-MASCARET system, developed for modelling suspended sediment transport processes. SEDI-3D consists of two parts, one for calculating the suspended sediment movements, and the other one for bottom evolution based on the assumptions of a bed layer model. The suspended sediment is considered an active tracer. There is the possibility to model cohesive sediment, non-cohesive or a mixture of both. For the modeling of only cohesive or only non-cohesive transport only one size class of sediments can be used at this time (V7P2). In this report the focus is on the cohesive sediment modelling. Delwaq was first preferred to do cohesive sediment modelling within an offline approach. Although offline modeling is assumed to be faster since it only involves modeling a single linear advection-diffusion equation, still simulation times were too long, merely due to the lack of the ability to parallel computing. Therefore, Delwaq was considered as not ideal within this project. An offline approach using SISYPHE was also left after problems with the boundaries. A more elaborate explanation on the problems of both Delwaq and SISYPHE is given in Appendix A. And thus SEDI-3D was chosen to model cohesive sediment transport.

Conceptual Model 5.2.1 Advection-Diffusion In nature, cohesive sediment transport occurs in fluids (water column) through the combination of advection and diffusion. In SEDI-3D, a 3D advection-diffusion equation is solved by considering the cohesive sediment particles moving at the same velocity as the fluid: 𝜕𝜕 𝜈𝜈𝑡𝑡 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 � + 𝑈𝑈𝑗𝑗 = + 𝑤𝑤𝑠𝑠 𝐶𝐶𝛿𝛿𝑖𝑖3 � 𝜕𝜕𝑥𝑥𝑗𝑗 𝜕𝜕𝑥𝑥𝑖𝑖 𝜎𝜎𝑡𝑡 𝜕𝜕𝑥𝑥𝑖𝑖 𝜕𝜕𝜕𝜕

In this equation U is the mean flow velocity [m/s], t is the time [s], xj represents the components of the coordinate vector [m], vt is the eddy viscosity [m2/s], σt is the turbulent Prandtl-Schmidt number (i.e. the ratio of vt to the eddy diffusivity of the sediment particles), C is the sediment concentration [g/L or kg/m³], ws is the representative mean settling velocity [m/s], and δij is the Kronecker delta. 5.2.2 Bed Shear Stress The bed shear stress is given by 𝜏𝜏𝑏𝑏 = 𝜌𝜌𝑤𝑤 𝑢𝑢∗ |𝑢𝑢∗ |

with 𝜌𝜌𝑤𝑤 the density of the water and 𝑢𝑢∗ the friction velocity.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

In SEDI-3D, a quadratic friction law is used with a drag coefficient CD to compute τb in a rough regime. When a Manning coefficient is used the equations look as follows: 1 𝜌𝜌 𝐶𝐶 𝑈𝑈�𝑈𝑈� 2 𝑤𝑤 𝐷𝐷 𝑔𝑔 𝐶𝐶𝐷𝐷 = 2𝑛𝑛2 1/3 ℎ

𝜏𝜏𝑏𝑏 =

Where 𝑈𝑈 is the depth-averaged velocity (which is also calculated in SEDI-3D), n is the Manning coefficient, g is gravitational constant and h is the water depth. After the calculation of this shear stress, the shear velocity is calculated and is then imposed at the bottom as a boundary condition for solving the momentum conservation equations of the flow. 5.2.3 Erosion At the interface between the water column and the bed layer, erosion processes happen due to the shear motion of the flow. The erosion flux is computed with the Partheniades formula. The erosion flux is the product of an erosion rate multiplied with a probability factor as a function of the shear stress in excess of a critical erosion shear stress: 𝜏𝜏𝑏𝑏 𝑀𝑀 � − 1� 𝑖𝑖𝑖𝑖 𝜏𝜏𝑏𝑏 > 𝜏𝜏𝑐𝑐𝑐𝑐 𝐸𝐸 = � 𝜏𝜏𝑐𝑐𝑐𝑐 0 𝑜𝑜𝑜𝑜ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒

with M the Krone-Partheniades erosion constant [kg/m²/s] and the value of M can be set by the user via the keyword EROSION COEFFICIENT, τb the bed shear stress and τce the critical bed shear stress for erosion. The latter can also be set by the user by using the keyword CRITICAL EROSION SHEAR STRESS OF THE MUD LAYERS. So erosion only occurs when the bed shear stress is higher than the critical bed shear stress for erosion set by the user. The erosion constant M determines the intensity of the erosion. A larger value will mean more erosion if erosion occurs. 5.2.4 Deposition The empirical deposition law from Krone is implemented in SEDI-3D to estimate mud deposition. Here the deposition flux is approximated by the product of local sediment concentration with the settling velocity, multiplied with a deposition probability: 𝜏𝜏𝑏𝑏 � 𝑖𝑖𝑖𝑖 𝜏𝜏𝑏𝑏 < 𝜏𝜏𝑐𝑐𝑐𝑐 𝑤𝑤𝑠𝑠 𝐶𝐶 �1 − 𝐷𝐷 = � 𝜏𝜏𝑐𝑐𝑐𝑐 0 𝑜𝑜𝑜𝑜ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒

Where 𝜏𝜏𝑐𝑐𝑐𝑐 is the critical shear stress for mud deposition (keyword CRITICAL SHEAR STRESS FOR DEPOSITION), ws is the settling velocity [m/s] (keyword CONSTANT SEDIMENT SETTLING VELOCITY), and C is the sediment concentration in suspension [g/L] or [kg/m³]. If the bottom shear stress is smaller than the critical bottom shear stress for deposition, sediment is settling.

Within this project we choose to model deposition D as a shear stress independent flux, following Sanford and Halka (1993) and Winterwerp (2007). This is also in line with recent applications in modelling cohesive sediment transport (Le Hir, 2011; Van Maren, 2015). This is done by setting 𝜏𝜏𝑐𝑐𝑐𝑐 to a very large value of 1 000 000 Pa. The formula for the deposition flux then simplifies to: 𝐷𝐷 = 𝑤𝑤𝑠𝑠 𝐶𝐶

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

5.2.5 Bed evolution Bed evolution is implemented in SEDI-3D via the Exner equation (without considering the bedload transport at bottom): (1 − 𝜆𝜆)

𝜕𝜕𝑧𝑧𝑏𝑏 + (𝐸𝐸 − 𝐷𝐷) = 0 𝜕𝜕𝜕𝜕

where λ is the bed porosity and zb is the bed level. 5.2.6 Bed Structure

SEDI-3D provides a multilayer bed model for simulating complex domains with different properties in different bed layers (Figure 32). In this way, the sediment bed is represented by a fixed number of layers. Each layer is characterized by its thickness, sediment concentration and resistance to erosion, i.e. a critical shear stress for erosion. Usually the first layer (contains fresh deposits) has the lowest mud concentration and a small critical shear stress for erosion. Figure 32 – Representation of the bed structure in SEDI-3D (after Tassi, 2017)

When the bed shear stress exceeds the critical value, the first layer starts to erode. Depending on the thickness of the first layer, only after it is completely eroded, the second layer becomes erodible. For the deposition at the bottom, the deposits always have the same properties as the material in the first (= upper) layer. If the first layer is completely gone, then it will restore the first layer when new deposits are added to the bed. If needed a consolidation model is present to refill the lower layers. 5.2.7 Density SEDI-3D is capable of including both sediment and salinity in the model, and the water density will be affected by sediment concentration, as well as salinity. This is done through a relative density ∆𝜌𝜌 defined in SEDI-3D as: 𝜌𝜌𝑠𝑠𝑠𝑠𝑙𝑙 − 𝜌𝜌0 𝜌𝜌𝑠𝑠𝑠𝑠𝑠𝑠 − 𝜌𝜌0 � ∆𝜌𝜌 = + 𝐶𝐶 � 𝜌𝜌0 𝜌𝜌𝑠𝑠𝑠𝑠𝑠𝑠 𝜌𝜌0

Where ρ is the density affected by the salinity and suspended sediment, ρ0 is the original density of water, ρsal is the density of the saline water, ρsed is the density of sediment, and C is the concentration of the sediment in suspension.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

There are several density laws provided in SEDI-3D and the following one is used when salinity is presented: 𝜌𝜌𝑠𝑠𝑠𝑠𝑠𝑠 = 𝜌𝜌𝑟𝑟𝑟𝑟𝑟𝑟 (1 + 750 ∗ 𝑆𝑆 ∗ 10−6 ) 𝜌𝜌𝑟𝑟𝑟𝑟𝑟𝑟 = 999.972 𝑘𝑘𝑘𝑘/𝑚𝑚3

with S the salinity [kg/m³] and ρref the reference density.

Then the relative density is used to calculate the buoyancy source terms arising from the density gradient in momentum conservation equations of the flow in horizontal directions: 𝐹𝐹𝑥𝑥𝑥𝑥 = −𝑔𝑔

𝑍𝑍𝑠𝑠 𝜕𝜕 ∆𝜌𝜌 �� 𝑑𝑑𝑑𝑑� 𝜕𝜕𝑥𝑥𝑖𝑖 𝑧𝑧 𝜌𝜌0

where Fxi is the buoyancy source term, z is the water elevation and Zs is the elevation of the free surface. In this sense, the salinity and sediment concentrations will have an influence on the flow field, and this will in turn affect the sediment movement again.

Sediment properties For the sediment transport, only one class of mud is considered at this moment. Based on the relevant literature (van Leussen, 1994; Winterwerp, 2002; Lee et al., 2011), the following mud properties are adopted and considered to be sufficient in terms of reproducing the general patterns of mud transport in the Scheldt estuary. The characteristic diameter of mud particles used in the model is set to 50 µm. Lee et al. (2011) summarized the settling velocity of mud particles related to its size (Figure 33). Combining the data from settling column tests (van Leussen, 1994) and the modified Stokes equation (Winterwerp, 2002), the settling velocity in this reference case is set to 0.5 mm/s. The density is set to 2650 kg/m³. Flocculation and hindered settling are not active. Figure 33 – Plots of floc diameter versus settling velocity. The diamond symbols represent the measured diameter and settling velocity of a floc in settling column tests (van Leussen, 1994) and the lines represent the simulated data with the modified Stokes equation (Winterwerp, 2002).

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Another important group of physical parameters is for the erosion/deposition process at bottom. In the current model, one bed layer is assumed at the bottom. This bed layer is initially empty. The density of mud in this layer is set to 500 kg/m³. Since these simulations are only concerning mud transport and no morphology, consolidation is not taken into account. For freshly deposited sediment at the bottom the critical erosion shear stress is 0.05 Pa (cfr. Mitchener and Torfs, 1996). The critical shear stress for deposition is set to a very large value (1 000 000 Pa) to ensure deposition of sediment is happening constantly. The erosion coefficient is set to 1.0E-4 kg/m²/s. These values are given to specific keywords in the Telemac3D steering file. The list of keywords describing these sediment and bed layer properties can be found in Table 8. The values for settling velocity, erosion coefficient and critical erosion shear stress were reached after calibration of the model. Table 8 – Keywords TELEMAC-3D for sediment properties

keyword

Value

Explanation

SEDIMENT

YES

Activate sediment transport in TELEMAC-3D

COHESIVE SEDIMENT

YES

Choose between cohesive or non-cohesive sediment

MEAN DIAMETER OF THE SEDIMENT

0.00005 m

The diameter D50 for sediments

DENSITY OF THE SEDIMENT

2650.0 kg/m³

Sediment density

CONSTANT SEDIMENT SETTLING VELOCITY

0.5E-3 m/s

The settling velocity of the sediment

NUMBER OF SEDIMENT BED LAYERS

The total number of sediment bed layers

INITIAL THICKNESS OF SEDIMENT LAYERS

Sediment layers thickness for initialization

MUD CONCENTRATIONS PER LAYER

500 kg/m3

Dry density of the mud-bed layers

EROSION COEFFICIENT

1.0E-4 kg/m²/s

The erosion coefficient used in Partheniades formula

CRITICAL EROSION SHEAR STRESS OF THE MUD LAYERS

0.05 N/m²

Critical erosion shear stress of the mud per layer

CRITICAL SHEAR STRESS FOR DEPOSITION

100,000 N/m²

Value of the critical bottom shear stress under which deposition occurs.

Boundary conditions 5.4.1 Hydrodynamic boundary conditions The downstream hydrodynamic boundary is water level driven. Water level time series are extracted from a harmonic ZUNO run representing the year 2013. For the mud model a period of 42 days is chosen. 2 days for hydrodynamic spin-up; 20 days for the spin-up of the sediment (see section 6); and a 20 day production run (Figure 34). The harmonic time series are corrected after comparison with harmonic analysis of measurements over a period of one year (Maximova et al., 2015). The 42 day period starts at 29/07/2013 and stops at 07/09/2013. This period is chosen because it corresponds to the calibration period of the hydrodynamic model.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 34 – The synthetic boundary for the calibrated model

Production run 20 days

Spin-up period 22 days

Upstream tributaries feed the model with a discharge. The discharges are representative discharges for a summer period and remain constant with an event of five days starting after 34 days. This discharge event is synthetic and is representative for an event with return period 1/6 (synthetic boundary discharge values determined by IMDC, 2015). This boundary condition is called the Normal Discharge scenario (QN). The discharges are shown in Figure 35 for all tributaries and upstream boundaries. The names and location of all these boundaries can be found in Figure 11 in section 3.1. The choice to include a discharge event was made by the steering committee of this project. Figure 35 – Discharges QN boundary condition

The location were discharge enters the model upstream near Merelbeke has changed a little compared to the calibrated hydrodynamic model described in Smolders et al. (2016). The new location correspond to the real location of the weirs at Merelbeke and ensures a better introduction of the sediment in this part of the model. It is shown in Figure 36. 40

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 36 – New location of Q boundary at Merelbeke

5.4.2 Sediment input from the boundaries In the model a sediment concentration is imposed to every liquid boundary. Every boundary has a different but constant sediment concentration imposed. Van Hoestenberghe et al. (2014) calculated the average annual total sediment load at Schelle for the period 1971 – 2009 and estimated this to be 287.45 tons according to the interpolation method they present. Table 9 is made from data coming from the interpolation method from Van Hoestenberghe et al. (2014) and shows the yearly averaged sediment load for the different subcatchments of the Scheldt. The sediment loads of the subcatchments Rupel, Nete, Lower Scheldt and Durme, which do not have a liquid boundary in the Scaldis model, was redistributed over the other subcatchments to maintain the total upstream sediment load of the estuary. The sediment load of the Nete was redistributed over the Kleine Nete and Grote Nete according to their weight in the contribution to the total load. In the same way the sediment load of the Rupel was redistributed over Kleine Nete, Grote Nete, Dijle and Zenne. The sediment load of the Lower Sea Scheldt and Durme was redistributed over the Dender and Merelbeke. The calculated yearly averaged redistributed sediment loads are given in the third column of Table 9. Table 9 – Sediment loads per subcatchment

subcatchment Dender Zenne Kleine Nete Dijle Grote Nete Merelbeke (Upper Scheldt) Rupel Nete Lower Scheldt Durme

Final version

Total (Schelle)

sediment load (tons/year) 31.95 32.80 12.29 81.61 10.94 91.20 6.11

2.72 10.76 7.36 287.45

redistributed sed. load (tons/year) 36.65 32.80 14.29 85.24 12.38 104.62 to Kleine and Grote Nete, Dijle and Zenne to Kleine and Grote Nete to Merelbeke and Dender to Merelbeke and Dender

WL2020R13_131_6

sediment load (tons/388.375 days) 38.60 34.54 15.05 89.77 13.03 110.17

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

To get the concentrations for every boundary in the Scaldis mud model, the total loads are divided by the total discharge volumes of the simulation period. The simulation period will be the same as for the hydrodynamic scenarios, meaning this will be four times a QN period followed by one QE period. Together these statistical boundary conditions represent an average year. For more specific information about these QN and QE boundary conditions we refer to Smolders et al. (2017). These statistical boundary conditions are a few days longer than a normal year, i.e. 388.375 days to be exactly. For each boundary (subcatchment) the total annual discharge volume was calculated from the synthetic boundary conditions (4*QN + QE) given in the second column of Table 10. Because this boundary condition is not exactly one year, but 388.375 days, the sediment loads were recalculated to represent the total sediment load over 388.375 days (last column of Table 9). Dividing the total loads by the total discharge volume gives the mass concentrations (Table 10) per catchment. These values are used as the sediment boundary condition for the mud model. They remain constant in time. The sediment concentrations for the other two discharge boundaries, i.e. Bath and Terneuzen, are zero in the model. Table 10 – Total discharge and average sediment concentration per subcatchment

Subcatchment

total discharge volume (m³)

mass concentration (g/L)

Dender

3.9E8

0.098

Zenne

5.6E8

0.062

Kleine Nete

3.6E8

0.041

Dijle

1.2E9

0.074

Grote Nete

2.9E8

0.045

Merelbeke

1.2E9

0.094

For the downstream boundary satellite images from Fettweis et al. (2007) were used. These satellite images were already used for a mud model for Zeebrugge (Nguyen et al., 2019). Figure 37 – Annual mean SSC in the North sea with the location of the Scaldis model sea boundary the annual mean SSC values (in mg/L) in the North Sea. The downstream boundary of the Scaldis model is drawn as a grey line. On this line the values for SSC could be extracted. These values can be imposed to the model varying in space like the water level time series vary in space along the sea boundary. But because this boundary is far from our zone of interest (i.e. the upstream part of the estuary) and simulation times will be limited (weeks instead of years), an average value is taken to represent SSC values at the downstream boundary, i.e. 0.013 g/L. This value is kept constant over space and time for the downstream boundary.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 37 – Annual mean SSC in the North sea with the location of the Scaldis model sea boundary (source: KBIN – OD Natuur)

Sediment in Telemac3D is treated like an active tracer. Because of hard coded software choices, sediment must always be the last tracer in the list. So in the steering file at the keywords for tracers a second tracer (after salinity) should be added. These values are summarized in Table 11. Table 11 – Keywords for tracer handling in Telemac3D

Keyword

Value

Explanation

NUMBER OF TRACERS

Cohesive sediment is added as the last tracer. Salinity was already present as an active tracer.

NAMES OF TRACERS

'SALINITY 'MUD

INITIAL VALUES OF TRACERS

0.0; 0.5

The second value initializes the mud concentration (kg/m³ or g/L)

PRESCRIBED TRACERS VALUES

0.1; 0.0; 0.1; 0.094; 0.1; 0.098; 0.1; 0.062; 0.1; 0.074; 0.1; 0.045; 0.1; 0.042; 0.1; 0.0; 0.1; 0.012

Values are given for the boundaries. For boundary 1 the value for tracer 1 and 2 is given, then for boundary 2 the value of tracer 1 and 2 are given, and so on.

Final version

PSU G/L

'; '

Define the names and units of the tracers

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Specifies the type of profiles of tracer concentration on the vertical. Possible choices are: 0: user defined, TRACERS VERTICAL PROFILES

1;1;1;1;1;1;1;1;1

1: constant, 2: Rouse equilibrium, constant (diluted tracer) or Rouse (sediment), 3: Rouse (normalized) and imposed concentration. 4: Rouse modified with molecular viscosity

COEFFICIENT FOR HORIZONTAL DIFFUSION OF TRACERS

1.0E-06;1.0E-06

The first value is for the salinity as set before for the Scaldis 3D calibrated model. The second value is for the sediment.

COEFFICIENT FOR VERTICAL DIFFUSION OF TRACERS

1.0E-06;1.0E-06

The first value is for the salinity as set before for the Scaldis 3D calibrated model. The second value is for the sediment.

Initial conditions Hydrodynamics are initialized with a 2 day simulation without sediment. Then a new simulation with sediment is started, starting from the last time step of this initialization run. In this new simulation sediment is initialized in the water column. Different values were tried and 0.5 g/L was chosen to initialize mud in the model. Initializing sediment in the water column or on the bed has no influence on the final equilibrium solution of the model as will be shown in section 6. (and in the assumption that the same amount of sediment is used in both types of initialization). It was chosen to start from sediment in the water column. If an unlimited supply of sediment was present on the bed, this would give better results in terms of SSC, but erosion rates would be unrealistically high. It was chosen not to use this approach and stay closer to reality. Salinity is initialized together with hydrodynamics (Smolders et al., 2016). Measured values are interpolated on the model grid. In the 2050 model later or in the scenarios, the same salinity map will be used to initialize salinity in all model runs together with hydrodynamics.

Handling tidal flats To handle tidal flats in TELEMAC there are a few options to choose from. Depending on the treatment of tidal flats several subroutines make different choices in handling tracers and water level on these tidal flats. The most important is to keep the choices consistent over all simulations. The different keywords handling these choices are listed in Table 12 and the values for the mud model are also given here. Appendix B gives a list of the most common locations in the code where these choices for tidal flats can have an influence.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Table 12 – Keywords for handling tidal flats in Telemac3D

Keyword

Value

Explanation

TIDAL FLATS

YES

Activate if tidal flats are present in the model

Algorithms of treating tidal flats: 1: equations solved everywhere with correction on tidal flats, 2: dry elements frozen

TREATMENT OF NEGATIVE DEPTHS

It is only used if the keyword OPTION FOR THE TREATMENT OF TIDAL FLATS = 1. Possible choices are: 0: no treatment, 1: smoothing, 2: flux control

MINIMAL VALUE FOR DEPTH

0.1 m

the minimum water depth, below which tidal flats are considered as dry.

THRESHOLD FOR SEDIMENT FLUX CORRECTION ON TIDAL FLATS

0.1 m

Below this limiting depth, all sediment erosion rates are set to zero.

Treatment of tidal flats at the diffusion step for tracers. 0: forced to zero, 1: value before masked

OPTION FOR THE TREATMENT OF TIDAL FLATS

TREATMENT ON TIDAL FLATS FOR VELOCITIES

TREATMENT ON TIDAL FLATS FOR TRACERS

Disabling bottom update in SEDI-3D In this model approach the bottom is not updated in SEDI-3D for the hydrodynamics. The simulations are short and morphodynamics are not the focus of this study. By preventing bottom elevation updating, there is no feedback of changing bottom elevation to the modelled hydrodynamics. The bottom is however updated for the sediment module, so it keeps track of the amount of sediment in the bottom, which has a feedback tot the erosion and sedimentation rates. In SEDI-3D (TELEMAC v7p2r1), there are several variables involved in the computation of bed evolution. The general process is as follows: 1. Initialize the bed with certain thickness; 2. Compute the erosion/deposition flux at bottom, determine the gain/loss of bed material at each point; 3. The mass changes at each location is translated in to the changes of bed thickness 4. The changes of bed thickness is added to the original bottom elevation, therefore the bottom is updated.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

The bed evolution is calculated in the subroutine fonvas.f. The bottom changes in the steps 1 to 3 are always calculated for the sediment module SEDI-3D. This kind of information can be saved in output variables like bed thickness and bed evolution. At step 4 we prevent the bed changes to be added/subtracted to/from the original bottom elevation. This is done by commenting the following line (at line 310 in V7P2r1 for mixed sediment and at line 409 for a single class sediment case) in the subroutine fonvas.f: CALL OV( 'X=Y+Z ', ZF , ZR,HDEP,C,NPOIN2)

Decoupling the effect of SSC on relative density In SEDI-3D, a density law can be used to model the influence of salinity and sediment concentration on water density. The changes in water density by salinity and SSC could induce changes in the flow field, through buoyancy terms in the momentum conservation equations. However, the sediment concentration is considered to be small and cannot alter the flow field significantly. In order to simplify the physics implemented in the model the coupling between the suspended sediment and the hydrodynamics throughout density changes will be switched off. A relative density Δρ is calculated in the subroutine drsurr.f as mentioned earlier in section 5.2.7. To decouple the effects of SSC on relative density the following lines in the subroutine drsurr.f (from line 197 to line 207 in V7P2r1) are commented. This prevents SSC effects from being added to the relative density. IF(SEDI) THEN IF(MIXTE) THEN CALL OS('X=X+CY ',X=DELTAR,Y=TA%ADR(NTRAC-1)%P, &

C=(RHOS-RHO0)/(RHO0*RHOS)) CALL OS('X=X+CY ',X=DELTAR,Y=TA%ADR(NTRAC)%P,

C=(RHOS-RHO0)/(RHO0*RHOS)) ELSE CALL OS('X=X+CY ',X=DELTAR,Y=TA%ADR(NTRAC)%P,

C=(RHOS-RHO0)/(RHO0*RHOS)) ENDIF

ENDIF

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Using a constant Manning Bottom Roughness for sediment transport The hydrodynamic model was calibrated using spatial variable Manning bottom roughness coefficients. Because these coefficients are used to calibrate and they are correcting for more than only the differences in bottom roughness, e.g. a turbulence model that might be too dissipative in more tortuous parts of the estuary, it would give undesired patterns in the results if they were used for the calculation of the bottom shear stress used to calculate erosion rates. Figure 38 shows the values of the Manning bottom roughness coefficient for the Scaldis HD 2013 model along the estuary axis. Figure 38 – Manning bottom roughness coefficient of Scaldis HD 2013 along the estuary axis.

The Fortran code of the CLSEDI.f subroutine was changed so that for sediment transport only a fixed Manning coefficient of 0.02 s/m1/3 is used for the entire model domain. DO IPOIN=1,NPOIN2 !

COMPUTES THE FLUID DENSITY DENSI(IPOIN) = (DELTAR(IPOIN)+1.D0)*RHO0

COMPUTES THE STRESS AT THE BOTTOM FRICCO=0.02D0 HC = MAX(HN%R(IPOIN),1.D-1) CFS = 2*FRICCO**2*GRAV/HC**(1.D0/3.D0) UETCARS = (U2D%R(IPOIN)**2+V2D%R(IPOIN)**2)*0.5D0*CFS TOB%R(IPOIN) = DENSI(IPOIN)*UETCARS ENDDO

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

The effect of a different Manning coefficient on the calculated bed shear stress in the model is shown in Figure 39. For the same depth and the same flow velocity, a higher Manning coefficient will yield in a higher bed shear stress in the model. Figure 39 – Effect of Manning coefficient on bed shear stress; varying flow velocities U and water depths.

Reduced settling velocity on shallow areas using a logistic function When τ critical deposition is set to a very high number (1000000 Pa), there will always be sediment deposition in the model and the equation will be rewritten to: D = ws C In this equation the settling velocity is constant for the entire model domain. A parameter alfa, α, is added to this equation. The hydrodynamic model does not capture some physical processes in shallow water (like wave action). The parameter alfa is a logistic function that reduces the settling velocity in shallow areas. The formula for deposition will then look like: D = α ws C where 𝛼𝛼 =

1 + 𝑒𝑒 −𝑘𝑘(𝑑𝑑−𝑑𝑑0 )

where d is the water depth, d0 is the water depth below which a significant reduction will take place and k determines the steepness of the slope in reducing α from 1 to 0. With k = 5, d0 = 1.5, and k = 5, d0 = 3.0, this sigmoid function is plotted in Figure 40.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 40 – sigmoid (logistic) curve alpha

In the calibrated model the threshold for depth (d0) is set to 1.5 m. This value was chosen based on the results of preliminary model runs that showed that beyond this threshold excessive sedimentation became a problem. Alfa was added to the settling velocity in the subroutine vitchu.f. This alfa is applied to every node in the model. The line of code is given here below. WCHU is the settling velocity in every point at the current time step. WCHUO is the settling velocity given by the user in the steering file. HN is the water depth at point IPOIN at the current time step. The threshold value for depth is given in red, i.e. 1.5 m. WCHU%R(IP) = WCHU0/(1+EXP(-5.0D0*(HN%R(IPOIN)-1.5D0)))

Reduced settling velocity on upstream boundaries and entrance of locks The bathymetry at the upstream boundaries is artificially deepened to create a more stable hydrodynamic inflow boundary. By making the river much deeper in these small sections the inflow velocity is low and the hydrodynamic run will be more stable. The low flow velocities are a problem for the sediment transport as these sections act as sediment traps. Therefore the settling velocity at these locations was set to 10-7 m/s. Some locations like the entrance to the old and new locks of Wintam and Duffelsluis are subjected to large sedimentation rates in reality. These locations are regularly swept, i.e. the excess of sediment is swept to a location just outside the entrance zone into the main navigation channel. Therefore the settling velocity of the mud was also set to 10-7 m/s at these locations. This was done in the subroutine vitchu.f where the settling velocity for all nodes that are located inside a polygon (polygon includes the areas mentioned above) is set to 10-7 m/s.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Dredging and Disposal Flux As a first approximation of dredging and disposal of sediment, the total disposal flux of sediment is added as a point source of sediment to the simulation. The magnitude of the sediment concentration of this point source is determined based on reported disposals in recent years (2007-2015). Figure 29 gives the average (over the period 2007-2015) total dredged sediment in the lower Sea Scheldt. On average 4.5 million tons dry solids (TDS) is deposited back in the estuary each year. Figure 31 gives insight in the origin of this material. It is important to mention that dredging through alternative dredging techniques such as the sweep beam are not included in the statistics mentioned in Figure 31. In the Scaldis model a point source is added with coordinates (RD Paris): x=83430 m and y=361424 m. The sediment is release with a discharge of 0.1 m³/s and a concentration of 1441.53 g/L at -6 m TAW. This corresponds to a release of 4.5 million tons TDS per year. The location of the source point in the model is shown in Figure 41. Real dredging and disposing subroutines are available within the TELEMAC modeling suite as the module NESTOR, but works currently only with SISYPHE. Figure 41 – Location of point source of sediment in the Scaldis model

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

6 Results Spin-up time How to initiate sediment in the model and how much time does it need to reach a quasi-equilibrium state? To answer these questions a simplified model simulation was done using a simple harmonic signal as boundary condition. 6.1.1 Set up harmonic run The model domain of the Scaldis model was changed: the downstream boundary was set at the entrance of the estuary (between Cadzand and Westkapelle) and can be seen in Figure 42. On this boundary it is easier to impose a harmonic water level signal. This signal was coded in the SL3.f subroutine, which was specially designed to be modified by end users for this purpose: PI = 4.D0*ATAN(1.D0) SL3 = 1.89D0*SIN(AT*(2.D0*PI/43200.D0)+(PI/2))+2.68D0 The harmonic signal resembles the characteristics of a signal near the estuary mouth (Vlissingen) in 2013. The amplitude is 1.89 m and the mean water level is set to 2.68 m TAW. There are two tidal cycles in 24 hours. For the upstream boundaries, each tributary has a prescribed constant discharge (in the steering file) for the fresh water inflow, the same as in the Scaldis hydrodynamic model, i.e. 23 m³/s for the channel GhentTerneuzen; 34.7 m³/s for Merelbeke; 11.1 m³/s for Dender ; 15.92 m³/s for Zenne; 34.6 m³/s for Dijle; 8.3 m³/s for Grote Nete; 10.38 m³/s for Kleine Nete and 35 m³/s for channel near Bath. These discharges remain constant during the simulation.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 42 – Model domain of Scaldis for Harmonic simulations

A time step of 4 seconds is used and the simulation period is 40 days. All the keywords in the steering file, i.e. the characteristics of the simulation, remained unchanged compared to the hydrodynamic Scaldis model. After some testing two harmonic simulations were defined: - HAR_5: zero sediment in the water column, all initialized on the bottom. A layer of 0.05 m of sediments was available on the bed. - HAR_6: zero sediment on the bed, all sediment initialized in the water column. An initial concentration of 3.38 g/L was given in the water column. With this concentration more or less the same amount of sediment is present in the model domain as when 0.05 m of sediment was put on the bottom. On the boundaries SSC concentrations are given as described for the Scaldis Mud model in Table 10 – Total discharge and average sediment concentration per subcatchment in section 5.4.2. These concentrations remain constant during the entire simulation. Both harmonic simulations have a two day hydrodynamic spin-up period before the sediment is introduced in the model domain. 6.1.2 Mass balance information from the listing file Information about the mass balance inside the model domain is provided by the TELEMAC code in the listing file. The listing file is a text file where information about the simulation can be stored. For both harmonic simulations information about the solvers and the mass balance and the boundary fluxes are recorded every ten minutes of the simulation period. With a Matlab script the information about the mass balance is extracted from this file and plotted. An example is given here below where the parameters: Mass of bed, Mass of the present time step and mass leaving the domain during this time step were extracted and plotted over time. Figure 43 shows the mass of the bed, mass in suspension and the mass leaving the

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

domain for both harmonic runs in function of time. The total mass is the sum of the mass on the bed, in suspension and the mass that left the domain. This line needs to remain constant if mass is conserved in the model. MASS BALANCE SEDIMENT BED MASS BALANCE AT TIME =

1531800

MASSE OF BED

TOTAL ERODED MASS SEDIMENT BED MASS BALANCE

10124452433.0264

4091964537.03698

(GAIN>0 LOSS<0): -5.421638488769531E-004

WATER MASS AT THE PREVIOUS TIME STEP

0.5864336E+10

MASS AT THE PRESENT TIME STEP

0.5864680E+10

MASS LEAVING THE DOMAIN DURING THIS TIME STEP :

-344306.6

ERROR ON THE MASS DURING THIS TIME STEP

0.1054432E-05

FLUX BOUNDARY

23.00000

M3/S

( >0 : ENTERING

<0 : EXITING )

FLUX BOUNDARY

34.70000

M3/S

( >0 : ENTERING

<0 : EXITING )

FLUX BOUNDARY

11.10000

M3/S

( >0 : ENTERING

<0 : EXITING )

FLUX BOUNDARY

15.92000

M3/S

( >0 : ENTERING

<0 : EXITING )

FLUX BOUNDARY

34.60000

M3/S

( >0 : ENTERING

<0 : EXITING )

FLUX BOUNDARY

8.300000

M3/S

( >0 : ENTERING

<0 : EXITING )

FLUX BOUNDARY

85914.03

M3/S

( >0 : ENTERING

<0 : EXITING )

FLUX BOUNDARY

35.00000

M3/S

( >0 : ENTERING

<0 : EXITING )

FLUX BOUNDARY

0.000000

M3/S

( >0 : ENTERING

<0 : EXITING )

TRACER

1: SALINITY

, UNIT : PSU

* M3)

ADVECTIVE FLUX THROUGH BOUNDARIES OR SOURCES

-2749265.

DIFFUSIVE FLUX THROUGH THE BOUNDARIES

0.000000

MASS AT THE PREVIOUS TIME STEP

MASS AT THE PRESENT TIME STEP

0.1124493E+12

MASS EXITING (BOUNDARIES OR SOURCE)

-0.1099706E+08

ERROR ON THE MASS DURING THIS TIME STEP

14503.11

ADVECTIVE FLUX THROUGH THE BOUNDARIES

-1034.469

DIFFUSIVE FLUX + DEPOSITION

-85.29186

MASS AT THE PREVIOUS TIME STEP

0.3985857E+09

0.1124383E+12

SEDIMENT IN SUSPENSION

MASS AT THE PRESENT TIME STEP

MASS LEAVING THE DOMAIN DURING THIS TIME STEP : ERROR ON THE MASS DURING THIS TIME STEP

Final version

0.3985903E+09 -4479.044 -69.58301

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 43 – Mass balance plot for harmonic simulations HAR_5 and HAR_6

Figure 43 shows that the model reacts quickly (2 days) on the sediment input and that both simulations (sediment on bottom and sediment in suspension) evolve towards the same equilibrium condition as proven by Schramkowski et al. (2019). This graph is made considering the entire model domain. The total mass remains constant and shows that mass is conserved in the model. From this graph a spin up time of 20 days for the mud model is derived. The largest changes happen in the first two days, but after 20 days the simulation reaches values that change only slightly in the next simulation days. This is further shown in Figure 45 where the bed evolution is plotted for both harmonic simulations and with respect to the Upper Sea Scheldt part of the model domain. Distinction is made between subtidal or channel area and intertidal area. The division into these two parts is done based on polygons determined in Vandenbruwaene et al. (2016). The name intertidal area in Figure 45 refers to both the intertidal flats and marsh area of the Upper Sea Scheldt. This area is marked red in Figure 44 (figure given as an example). The name channel refers to deep, average and shallow sub-tidal areas or the blue area in Figure 44. In Figure 45 the bed evolution from all nodes in the channel and on the intertidal areas are averaged and plotted in function of time for both harmonic simulations. Inside the channel the first reaction of the sediment in suspension is to sink to the bed from where it is eroded again. This erosion rate is very similar to the erosion rate from the simulation where the sediment was initially put on the bottom. There the bed evolution becomes negative because of the erosion. On the intertidal areas the bed evolution is also very similar between both harmonic simulation. The difference lies in the beginning of HAR_6 where the excessive amount of sediment in the water column sinks to the bottom and realizes a first initial bed layer on the intertidal areas. Also in this figure the rates of bed evolution don’t change drastically anymore after 20 days.  conclusion: The way sediment is initialized in the model has no effect on the final solution, nor on the time scale to reach a kind of equilibrium situation in the model.  conclusion: A spin up time of 20 days should suffice for the sediment in the model to reach a kind of equilibrium condition.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 44 – An example of intertidal areas (red) and subtidal channels (blue) indicated by the polygons used in Vandenbruwaene (2016).

Figure 45 – Averaged bed evolution for tidal flats and inside the main channel for the Upper Sea Scheldt

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

6.1.3 Time scales to reach equilibrium: conceptual model The evolution towards equilibrium of the model shows two time scales. There is the horizontal equilibrium which is shown in the previous section. And there is a shorter timescale to reach vertical equilibrium between erosion and deposition. This is determined by the settling velocity, erosion coefficient (M = Partheniades constant) and the critical shear stress for erosion. In order to investigate this timescale, a conceptual model is used. This is a model that describes the erosion-deposition process in one point using the values and volumes of the entire TELEMAC model. The basic code lines are given: total_water_vol = 9.94e10 m³ total_area = 5.32e09 m² mud_layer_density = 500 kg/m³ initial_conc = 1.338 g/L bed_thickness = 0.0 m total_bed = total_area*bed_thickness*mud_layer_density tidal_waves = 1.89*sin(t*(2.0*pi/43200.0)+(pi/2.0))+2.68; tau_ce = 0.05 Pa M = 1.0e-5 Ws = 0.0005 m/s tau_b = 0.1*(tidal_waves-0.78) Pa erosion = max(0.0, M*(tau_b./tau_ce-1.0)*total_area) (in [kg/s]) deposition = concentration * ws * total_area (in [kg/s]) concentration(t+1) = concentration(t) - deposition(t)/total_water_vol + erosion(t)/total_water_vol; updating bed:

bed(t+1) = max(0.0, bed(t) + deposition(t) - erosion_1(t));

The bed shear stress tau_b in this conceptual model is dependent on the water depth and this equation is constructed so that it replicates the most common observed bed shear stresses (between 0 and 0.38 Pa) from the harmonic model. Simulations are done for three values of settling velocity (0.25 mm/s, 0.5 mm/s, 0.75 mm/s) and three values of erosion parameter M (0.5E-5 kg/m²/s, 1E-5 kg/m²/s, 1.5 kg/m²/s). The results are shown in Figure 46. This figure shows the effect of both parameters on the vertical time scale of the sediment in the model. The higher the settling velocity, the faster a vertical equilibrium is reached. The erosion parameter M has no big effect on the time scale but increases the amount of sediment in suspension. The conceptual model shows that time scales to reach an equilibrium condition in the vertical for sediment take less than 2 to 4 days dependent mainly on the settling velocity.  conclusion: a sediment equilibrium condition in the vertical is reached faster than in the horizontal, and needs only 2 to 4 days dependent mainly on the settling velocity.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 46 - The conceptual model showing the time scale for reaching the equilibrium between erosion and deposition.

Calibrated model results 6.2.1 Mass balance of the model The mass balance of the entire model domain was calculated from data (written during the simulation) in the listing file with mass balance output listed every ten minutes. The results are given in Figure 47 where the mass in suspension, mass on the bed, the mass leaving the domain and the total mass are plotted over the entire simulation time. The mass that is leaving the model domain is given a negative sign. After an initial two days the sediment behavior reaches a stable pattern in function of the tide. Figure 48 shows the mass balance of the Upper Sea Scheldt only for the last 20 days of the 40 day simulation. This figure shows that part of the sediment that enters the domain through the upstream boundaries is trapped on the bottom and part is flushed out of the Upper Sea Scheldt. The mass in suspension remains stable depending on the tide and spring/neap tidal cycle. The stronger increase in total mass in the last six days is caused by a higher input of sediment caused by a peak in upstream discharge. As a constant concentration of sediment is used on the upstream boundaries the amount of discharge determines the total load of sediment that enters the model domain.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 47 – Mass balance from listing file for the whole model domain

Figure 48 – Mass balance of the Upper Sea Scheldt of last 20 days of simulation (calculated from masses in UA polygons)

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

6.2.2 Mass transfer maps The mass transfer in the Upper Sea Scheldt has been studied before by Vandenbruwaene et al. (2017). The amount of sediment mass transfer in the OMES polygons was calculated based on differences in bathymetry and lithological data. The original values denote a transport rate over 10 years (2001-2011) but in Figure 49 these values are given in a value of transport per year (so divided by 10). For the current study, the same OMES polygons have been used to calculate the similar mass transfer map (Figure 49). The values coming from the mud model however were values of transport over a spring/neap tidal cycle. These values were extrapolated to transport rates per year so they can be compared with the values calculated by Vandenbruwaene et al. (2017). The directions of sediment transport on the OMES polygon boundaries are the same and the values are in the same order of magnitude in most of the polygons (Figure 49). But due to the different methods used, the directions are more meaningful. Figure 49 – Mud transport calculated from the model compared with calculated mud transport rate from Vandenbruwaene et al. (2017)

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

6.2.3 Dredging mass balance The mass silting up in access channels and docks is calculated and compared to the fixed source term that represents the disposal flux (§3.8). Because the model does not use a dredging and disposal functionality, these amounts should correspond as closely as possible. Figure 50 shows the sedimentation in the model at the end of the calibration run, together with the polygons that are used to calculate the total siltation in the access channels and docks. Note that the model does not reproduce the siltation that is observed on the location “Drempel van Frederik”. This dredging location is located in the fairway, north of Deurganckdok. It is the second most important dredging location with a dredged amount of about 1.1 MTDS/yr. The reason why the model does not reproduce siltation in this location is not well understood and requires further research. The siltation in Deurganckdok in the model is concentrated in the first half of the dock. The end of the dock is also dredged, however. Multiple processes that can cause internal redistribution of sediment are not included in the model, such as maneuvering ships that stir up the sediments, or gravitational circulation of high concentration sediment suspensions. Figure 50 – Location of polygons (red) from where the total amount of sediment on the bed is calculated

Figure 51 shows that the model predicts a siltation in the access channels of 175 kton in 20 days, or 3.2 MT/yr. In reality, the access channels have a total average siltation of 3.0 Mton/yr. In the model however, the total reported disposal flux of 4.5 Mton/year is introduced as a point source. The difference can be attributed to the fact that the model does not reproduce the siltation on the sills in the fairway. This means that an artificial source of 1.5 Mton/year is added in the model domain. This needs to be taken into account when analyzing the results.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 51 – Estimation of sedimentation in dredging polygons for last 20 days

175,000 ton in 20 days = 3.2 Mton/year

6.2.4 Averaged ETM The average (over spring/neap tidal cycle) SSC was calculated over the cross sections every 10 km along the estuary. The results are given in Figure 52. The peak in SSC around Antwerp (km 91; for locations and km see Figure 11) is mainly caused by the disposal of sediment near Oosterweel. The model uses yearly averaged discharges upstream and these discharges don’t represent a good summer (dry) or winter (wet) season in discharge. In the current model settings the sediment net transport direction is downstream and the model does not seem to be able to represent a natural ETM, which is located around km 120 in summer conditions and around Antwerp in winter conditions (Figure 20). In Plancke et al. (2019) an effort was made to reproduce a natural ETM in the Scheldt estuary in a Mike 11 1D model of the Scheldt estuary. Only with zero upstream discharge an ETM was visible. When the upstream discharge was introduced all sediment was flushed out of the system. In this model sediment was initialized with 0.1 m on the bed and this gives initially good results, but when this layer is depleted, the estuary is flushed of sediments. In Lanckriet and Conin (2018) the Delwaq cohesive sediment model of Deltares was used in the framework of the evaluation of external effects on the siltation of Deurganckdok. Also in this model, increased SSC concentration around Antwerp are strongly influenced by the disposal of sediments near Oosterweel and this also controls the location of the ETM in the model. In Brouwer et al. (2018) an idealized width-averaged process-based model, called iFlow, was used to investigate the role of a time-varying river discharge on the trapping of fine sediment in an ETM. Two different states for sediment on the bed were distinguished: in the availability-limited state, the SSC is limited by the amount of erodible sediment at the bed. Over time, under constant forcing conditions, the estuary evolves to morphodynamic equilibrium; in the erosion-limited state, there is an abundant amount of sediment at the bed so that sediment pickup occurs at the maximum possible rate. Under availabilitylimited conditions, periods of high river discharge push estuarine turbidity maxima (ETMs) downstream, while drier periods allow ETMs to move upstream. However, when the estuary is in an erosion-limited state during low river discharge, a bottom pool of sediment is formed. When the discharge then increases, it takes time to deplete this pool, so that an ETM located over a bottom pool moves with a significant time lag relative to changes in the river discharge. Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 52 – Averaged location of ETM with P25 and P75

Figure 53 shows the movement of the location of the peak in SSC for the different time steps along the estuary. The bottom panel shows the tide and the upstream discharge at Merelbeke. When the discharge increases near the end of the simulation the peak in SSC is pushed downstream. When the peak in discharge stops the peak in SSC moves quickly back upstream to its original location. Figure 53 – Location of ETM; visualization of all time steps

Figure 54 shows the average SSC values when the sediment source at Oosterweel is not active. This figure shows that the disposal of sediments near Oosterweel is causing the local peak in SSC. The model with its

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

current settings is not able to reproduce a natural ETM. In Smolders et al. (2019) different combinations of upstream discharges, critical erosion shear stress, critical disposal shear stress, Partheniades constant, settling velocity show no formation of a natural ETM in the Scheldt estuary. When sediments, mud, is put on the bottom as an unlimited source of sediments then better results are possible, but a continued high erosion rate in the navigation channel is not realistic and also does not represent the physical behavior of the sediment. Figure 54 – Effect of dumping of sediment near Oosterweel on the average SSC values along the estuary

6.2.5 Ensemble analysis SSC Both measured and modeled SSC values are timed with respect to the high water level at Antwerp tide station. Three locations with continuous measurements are discussed here: Buoy 84, Oosterweel and Driegoten (locations see Figure 11). Figure 55 shows the ensemble analysis of Buoy 84 location. The model seems to capture the average SSC, but clearly underestimates the peak after HW and especially during spring tide. Missing this peak during ebb tide might be explained by the fact that concentrations in this part of the estuary are very location dependent. During ebb the peak concentrations caused by the downstream transport (because of the ebb tide) of the disposed sediment near Oosterweel are better represented more towards the middle of the navigation channel. Figure 56 shows the ensemble analysis results for Oosterweel. Here the model tends to overestimate the SSC. Both model and measurement are highly influenced by the disposal of dredged sediments at this location. The way the sediment is added in the model as a point source might influence the order of magnitude of the SSC at this location. The model is able to represent both the ebb and flood peak in SSC for neap, average and spring tide. Figure 57 shows the ensemble analysis at Driegoten. This is the most upstream location and closest to the area of interest in this project. However, because the model with the current settings is not able to reproduce the natural ETM concentrations it delivers poor results at Driegoten. The flat line in model results shows that a kind of constant concentration of cohesive sediments is passing this location and is transported downstream. The measured variation in the signal is most probably due to the movement of a natural ETM in this area.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 55 – Ensemble analysis Buoy 84

Figure 56 – Ensemble analysis Oosterweel

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 57 – Ensemble analysis Driegoten

6.2.6 Sedimentation on intertidal areas in the Upper Sea Scheldt The siltation on the intertidal area in the Upper Sea Scheldt is analyzed. In order to avoid the influence from the spin-up period, the results of the last 20 days are used. The information about the intertidal areas in the Upper Sea Scheldt is from Vandenbruwaene et al. (2016), in which the intertidal areas are categorized as tidal flats and tidal marshes. The amount of sediment silted on the intertidal areas was calculated as follows: 𝑛𝑛

𝑛𝑛

𝑖𝑖=1

𝑆𝑆 = �� 𝑏𝑏𝑖𝑖 ∙ 𝐴𝐴𝑖𝑖 �� 𝐴𝐴𝑖𝑖

where, S is the averaged sedimentation rate during the last 20 days(cm/20d), b is the nodal value of bed evolution (cm) after 20 days, A is the nodal area derived from the domain discretization, i = 1, 2, …, n represents all the nodes in the intertidal areas indicated in in the Upper Sea Scheldt. Later this values is converted to the yearly-sedimentation rate (cm/yr) by multiplying by 18.25. Sedimentation rate on intertidal areas is on average 4.5 cm/year. But sedimentation in shallow areas was reduced by reducing strongly the settling velocity in these areas, so this average result is misleading. It is the average of a few locations with a lot of sedimentation and a lot of areas with few sedimentation; as is shown in Figure 58. This figure shows the small colored dots representing sedimentation along the navigation channel and shows no sedimentation inside the navigation channel. This model does not well represent the vertical dynamics of the sediment and is enhanced by interfering in the settling velocity in these areas. Figure 59 shows the average bed evolution in the model for the channel and the intertidal area of the Upper Sea Scheldt separately. This figure shows only results averaged over the entire domain of the Upper Sea Scheldt and specific locations can have sedimentation rates that differ from the averaged values.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 58 – Bed thickness after 40 days of simulation for the region Temse and Rupelmonde as an example

Figure 59 – Average bed evolution of the Upper Sea Scheldt

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Excessive sedimentation rates have been found in shallow areas in the model. Three possible solutions were considered to solve the problem: 1. the critical erosion shear stress is reduced to try to get more sediment back in suspension; 2. the settling velocity is reduced and made dependent on depth; 3. the bottom roughness was increased by setting a uniform Manning roughness coefficient of 0.02 s/m1/3. For a reference simulation the bed evolution shows that shallow areas have a positive bed evolution (=sedimentation) and the channel has a negative bed evolution (= erosion) (Figure 60). Only 20 days were simulated and sedimentation rates up to 7.5 m/year were found, whereas normal sedimentation rates for intertidal flats range from 0 to a few cm per year. Figure 60 – excessive sedimentation in shallow areas in SA_033 with sedimentation rates up to 7.5 m/yr.

The different solutions and reference are presented in Table 13. This table shows all the parameters that have changed between these possible solution simulations and the reference run. All other parameters were kept the same as in the calibrated mud model described above. The most important difference is that all these solution simulations start with sediment on the bed (this problem was tackled before the choice was made to start with sediment in suspension). In Solution 1 the critical erosion shear stress was lowered by adding a second bed layer that acts as a ‘fluffy layer’. In Solution 2 the settling velocity in shallow areas is reduced by adding a sigmoid function (note that this solution is used in the final calibrated model described above). In solution 3 a uniform bed roughness was introduced where especially for the Upper Sea Scheldt the local Manning bottom roughness value was increased and thus the calculated bed shear stress increased. Results are shown by means of sedimentation and erosion fluxes on a transect, located just downstream of Rupelmonde (i.e. the black transect line in Figure 60).

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Table 13 – Possible ways of reducing sedimentation on intertidal flats

Reference

Solution 1

Solution 2

Solution 3

Settling velocity

0.5 mm/s

sigmoid (k=5; d=3.5m)

0.5 mm/s

number of bed layers

thickness of bed layer

0.6 m

top 0.005 m

0.6 m

bottom 0.6 m

mud concentration per layer

500 kg/m3

100; 800 kg/m3

500 kg/m3

erosion coefficient

0.00001

0.00005

0.00001

critical erosion shear stress

0.05 Pa

top 0.01 Pa

0.05 Pa

Manning coefficient

like HD

uniform

Method

Run

bottom 0.15 Pa like HD

0.02 s/m1/3

reference run

reducing critical erosion by adding ‘fluffy layer’

reducing settling V with sigmoidal function

increasing bottom roughness

SA33

SA36

SA47

SA49

Figure 61 shows the results of the different possible solutions for the excessive sedimentation in shallow areas on the transect just downstream Rupelmonde. Figure 62 shows the same results but then zoomed in on the left bank of the estuary. Orange lines indicate sedimentation fluxes and blue lines indicate erosion fluxes. For the reference (SA33) it is shown that the sedimentation flux exceeds the erosion flux and the sedimentation flux reaches also further towards the river banks. In these parts the erosion flux is even zero. Increasing the bottom roughness by using the uniform Manning coefficient, solution 3 (SA49) increases the erosion flux, but does not create more erosion in the shallow areas. There the erosion flux remains zero. By decreasing the critical erosion shear stress, solution 1 (SA36), both the erosion flux and the sedimentation flux increase. The erosion flux is extended to the more shallow areas but also the sedimentation flux increases here, so this is not a good solution either. Finally, by decreasing the settling velocity with a sigmoid function (SA47), as described in section 5.10, the sedimentation flux reduces in the shallow areas and goes to zero. The erosion flux remains the same as in the reference simulation. Still the erosion flux is smaller than the sedimentation flux, but this can be solved by increasing the erosion parameter M.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary Figure 61 – Averaged erosion/deposition flux at a transect just downstream Rupelmonde for different solutions: orange color indicates deposition fluxes (D) and blue color indicates erosion fluxes (E).

Figure 62 – Zoom on left bank of Figure 61 - Averaged erosion/deposition flux at the left bank of the transect just downstream Rupelmonde for different SA runs - Orange color indicates deposition fluxes (D) and blue color indicates erosion fluxes (E).

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Simulation speed-up times Speed-up was tested with a simulation that was as similar as possible to the calibrations simulations. The simulation duration was two tides ( = 22350 time steps of 4 seconds). The computation was not started from a previous computation file. Culverts were included. Wind was not included. Cohesive sediment was included (SEDI-3D). Constant discharges were used for the upstream boundaries. A tidal signal was used for the downstream boundary. The Scaldis mesh and model domain was used (Smolders et al., 2016, 2017). Salinity was included as an active tracer. There was no update of the bottom for the sediment and the density of the water was also not affected by the sediment. The sediment was initialized as a concentration in the water in the beginning of the simulation and there was no sediment on the bottom. The model speed-up is defined as Tref/Tn (after Moulinec et al., 2011) where Tref is the simulation time with the smallest number of cores and Tn is the simulation time with n cores. 89400 seconds where simulated corresponding with two tides. The reference time was taken from a simulation on one Reynolds 1 node with 16 processors. Several simulations were started increasing the number of nodes used from the Reynolds queue. Figure 63 shows the scaling performance of a TELEMAC-3D simulation with the sediment module SEDI-3D on the Reynolds queue at Flanders Hydraulics. The black dashed line indicates the optimal speedup. Except for five Reynolds nodes (which corresponds to 80 cores) the curve keeps increasing, showing that using more processors will keep increasing the speed-up of the computations. Figure 63 – Scaling performance of TELEMAC-3D with SEDI-3D on different queues at FHR

The total computational capacity of the FHR Linux cluster is subdivided into different queues. Each queue has his own queuing subsystem and consist of a number of identical machines. The Reynolds queue consists of 12 nodes with each node holding 16 processors (2 x Intel Xeon E5-2640 v3 @ 2.60GHz 8 Core), 196 in total. 1

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

For comparison the same simulation was also performed on other queues at the Flanders Hydraulics computational cluster. Each queue contains nodes with the same hardware architecture, but hardware can differ between queues. Some of the queues had nodes with 12 processors and some with 16 so to compare simulations with the same number of processors, 48 processors were chosen. This means that on queues with nodes of 12 processors, four nodes are used instead of three. For the speed-up the Reynolds queue simulation on one node (16 processors) was chosen as the reference. Figure 63 shows that the Stokes queue is the best performing queue at Flanders Hydraulics closely followed by the Reynolds queue. The Navier and Stevin queues both perform less, but both have 12 processors per node instead of 16 like Reynolds and Stokes. They might be slower because of an increased communication between those nodes. The Bernoulli queue has only one node with 64 processors, of which 48 were used for this test. It scores terrible, being just a bit faster than reality. At the moment of writing this report, it was assumed this queue encounters a hardware problem. As a test, the same simulation was done on the Stokes queue, but the sediment module SEDI-3D was turned off. The simulation speed increased from six to almost nine time faster than reality. This shows that the sediment module weighs on the entire simulation.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

7 Conclusions The Scaldis Mud model is a simple cohesive sediment transport model built on top of a calibrated 3D hydrodynamics model. The most important simplifications are: • • • • •

Only one fraction is taken into account. This is a limitation of the software platform (Sedi3D). Details on the choice of the platform is given in appendix A. Hydrodynamics and Sediment transport were decoupled. There is no bottom update, no effect of SSC on the density or turbulence characteristics The settling velocity on the intertidal flats had to be artificially reduced in order to keep the sedimentation on those areas in check. Dredging and disposal is imposed as a flux, not calculated in a dredging/disposal module. This is a limitation of the software platform (Sedi3D). Flocculation and hindered settling were not taken into account

Considering the simplicity of this model, the Scaldis Mud model produces acceptable results. From the modest goals that were set in the introduction, most were reached: • • • • • • •

The model is able to reproduce a spatial SSC distribution with a local ETM, however this ETM is caused by local disposal of sediments. The natural ETM further upstream is not reproduced with the current settings of the model. There is a clear spring-neap tidal cycle in the SSC. The global mass balance reaches an equilibrium. The model is responsive to changes in SSC values at the upstream boundaries. The average sedimentation rate on the intertidal corresponds to measured rates. Sedimentation rates in the dredging polygons correspond to average dredging amounts. The model shows an acceptable tidal variation of SSC variation at Buoy 84 and Oosterweel.

The model however does not reproduce the natural ETM dynamics in the estuary. Note that currently, no complex model is able to reproduce the measured seasonal ETM dynamics. The advective transport however of suspended sediments the Upper Sea Scheldt from the upstream boundary to downstream of cohesive sediments introduced into the Scheldt estuary at the upstream boundaries is well modelled. The net fluxes are in the same order of magnitude as fluxes calculated based on bathymetrical data. This makes this model still usable for scenario analysis. Expert judgement is needed to interpret the results taking the limitations of this model into account.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

8 Recommendations for further improvements In this report the high sedimentation rate in the shallow areas is solved by reducing the settling velocity in these areas; in the LTV mud model in Delwaq this is tackled by introducing an additional bed shear stress caused by wind waves. A better understanding is necessary of the processes in the software code. Why is this problem occurring in the model and not in reality? Secondly, the model should be able to reproduce a natural ETM. Additional research already showed that with the current model settings, no parameter combination was able to get an upstream transport of cohesive sediments. The fact that also other models and software platforms are having problems modelling the seasonal ETM dynamics that can be observed in nature, demonstrates that this is not a trivial problem to solve. Finally, when the two problems above are solved the model should be improved locally on the smaller scale. The interaction with the hydrodynamic model should be investigated. At this moment the hydrodynamic model is calibrated separately from the sediment model, but it seems that both are more entangled; shown already by the effect of the bottom friction coefficient on the bed shear stress. The model can be made more complicated by adding more physical processes like consolidation and flocculation or adding more sediment classes.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

9 References Bi, Q.; Smolders, S.; Plancke, Y.; De Maerschalck, B.; Vanlede, J. (in preparation). Integraal Plan Bovenzeeschelde: Sub report 9 – Effect of B-alternatives on Mud Transport. Version 4.0. FHR Reports, 13_131_9. Flanders Hydraulics Research: Antwerp. Brouwer, R. L., G. P. Schramkowski, Y. M. Dijkstra, and H. M. Schuttelaars, 2018: Time Evolution of Estuarine Turbidity Maxima in Well-Mixed, Tidally Dominated Estuaries: The Role of Availability- and Erosion-Limited Conditions. J. Phys. Oceanogr., 48, 1629–1650, https://doi.org/10.1175/JPO-D-17-0183.1. Coen, L.; De Maerschalck, B.; Plancke, Y.; Verwaest, T.; Mostaert, F. (2016a). Sedimentstrategie BenedenZeeschelde: Deelrapport 1 – Opzet en validatie slibmodel. Versie 4.0. WL Rapporten, 14_025. Waterbouwkundig Laboratorium: Antwerpen, België. Coen, L.; De Maerschalck, B.; Plancke, Y.; Verwaest, T.; Mostaert, F. (2016b). Sedimentstrategie BenedenZeeschelde: Deelrapport 3 – Slibscenario’s. Versie 4.0. WL Rapporten, 14_025. Waterbouwkundig Laboratorium: Antwerpen, België. Dams, J.; Vanlede, J.; Plancke, Y.; Verwaest, T.; Mostaert, F. (2016). Slibbalans Zeeschelde: Deelrapport 3 Literatuurstudie. Version 4.0. WL Rapporten, 00_029. Flanders Hydraulics Research & Vrije Universiteit Brussel: Antwerp, Belgium. Fettweis, M., Nechad, B., & Van den Eynde, D. (2007). An estimate of the suspended particulate matter (SPM) transport in the southern North Sea using SeaWiFS images, in situ measurements and numerical model results. Continental Shelf Research, 27(10-11), 1568-1583. Hassan, W.; Maximova, T.; Bi, Q.; Smolders, S.; Vanlede, J.; Mostaert, F. (2019). Integraal Plan Bovenzeeschelde: Sub report 7 – Effect of the B-alternatives on the Hydrodynamics. Version 2.0. FHR Reports, 13_131_7. Flanders Hydraulics Research: Antwerp. IMDC (2003). Optimalisatie van onderhoudsbaggerwerken Deurganckdok; Deelrapport 1. Het hydraulisch model, I/RA/11239/03/049/EST. International Marine & Dredging Consultants iov. Gemeentelijk Havenbedrijf Antwerpen. IMDC (2004). Optimalisatie van de onderhoudsbaggerwerken Deurganckdok: deelrapport 2. Sedimentologisch en morfologisch model onderzoek. International Marine & Dredging Consultants iov. Gemeentelijk Havenbedrijf Antwerpen. IMDC (2014). Modelling instruments for Integrated Plan Upper Seascheldt. I/NO/11448/14.165/DDP IMDC (2015). Bepalen van randvoorwaarden voor de referentiesituatie (2050) m.b.t. debiet en waterstand. I/NO/11448/15.128/TFR/ Lanckriet T. and Cronin K. (2018). Evaluation of external effects on the siltation in Deurganckdok (20142017) Report 1.25: Delwaq model study of transient and external effects. IMDC and Deltares, Oktober 2018. Le Hir, P.; Cayocca, F.; Waeles, B. (2011). Dynamics of sand and mud mixtures: A multiprocess-based modelling strategy. Cont. Shelf Res. ISBN 0278-4343 31(10 SUPPL.). doi:10.1016/j.csr.2010.12.009

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Lee, B.J., Toorman, E., Molz, F.J. and Wang, J., 2011. A two-class population balance equation yielding bimodal flocculation of marine or estuarine sediments. Water research, 45(5), pp.2131-2145. Maximova, T.; Vanlede, J.; Verwaest, T; Mostaert, F. (2015). Vervolgonderzoek bevaarbaarheid Bovenzeeschelde: Subreport 4 – Modellentrein CSM – ZUNO: validatie 2013. WL Rapporten, 13_131. Flanders Hydraulics Research: Antwerp, Belgium Ministerie van de Vlaamse Gemeenschap (1991). De baggerwerken in de Schelde en de kwaliteit van water en bodem: stand van zaken. Ministerie van de Vlaamse Gemeenschap. Departement Leefmilieu en Infrastructuur: Brussel. 67 + appendices pp. Mitchener, H. and Torfs, H., 1996. Erosion of mud/sand mixtures. Coastal engineering, 29(1-2), pp.1-25. Moulinec, C., Denis, C., Pham, C. T., Rougé, D., Hervouet, J. M., Razafindrakoto, E., Barber, R.W.,Emerson, D.R., Gu, X.-J. (2011). TELEMAC: An efficient hydrodynamics suite for massively parallel architectures. Computers & Fluids, 51(1), 30-34. Nguyen, D.; Vanlede, J.; Scanu, A.; De Maerschalck, B.; Mostaert, F. (2019). Sediment Transport Model for the Port of Zeebrugge: Model set-up, Calibration, Validation. Version 1.0. FHR Reports, 00_067_3. Flanders Hydraulics Research: Antwerp. Partheniades, E. (1965). Erosion and deposition of cohesive soils. J. Hydraul. Div. 91(1): 105–139 Plancke, Y.; Taverniers, E.; Vanlede, J.; Verwaest, T.; Mostaert, F. (2012). Slibbalans-Zeeschelde: Deelrapport 1 - Metingen halftij-eb Boven-Zeeschelde 2012. Versie 2.0. WL Rapporten, 00_029. Waterbouwkundig Laboratorium: Antwerpen, België. Plancke, Y.; Vereecken, H.; Vanlede, J.; Verwaest, T.; Mostaert, F. (2014a). Slibbalans-Zeeschelde: Deelrapport 5 - Metingen halftij-eb Boven-Zeeschelde 2013. Versie 4.0. WL Rapporten, 00_029. Waterbouwkundig Laboratorium: Antwerpen, België. Plancke, Y.; Vereecken, H.; Vanlede, J.; Verwaest, T.; Mostaert, F. (2014b). Integraal plan BovenZeeschelde: Deelrapport 3 – Factual data-rapport sediment transport metingen in de Zeeschelde in 2014. Versie 4.0. WL Rapporten, 13_131. Waterbouwkundig Laboratorium: Antwerpen, België. Plancke, Y.; Vereecken, H.; Vanlede, J.; Verwaest, T.; Mostaert, F. (2015). Slibbalans-Zeeschelde: Deelrapport 7 - Metingen halftij-eb Boven-Zeeschelde 2014. Versie 4.0. WL Rapporten, 00_029. Waterbouwkundig Laboratorium: Antwerpen, België. Plancke, Y.; Coen, L.; Kaptein, S.; Mostaert, F. (2019). AvdT – Sedimenttransport op verschillende tijdschalen: Deelrapport 2 – Verkenning 1D-modellering slibtransport in het Schelde-estuarium. Versie 1.0. WL Rapporten, 17_088_2. Waterbouwkundig Laboratorium: Antwerpen.(in preparation) Sanford, L.P.; Halka, J.P. (1993). Assessing the paradigm of mutually exclusive erosion and deposition of mud, with examples from upper Chesapeake Bay. Mar. Geol. 114: 37–57 Santermans, J. (IMDC). (2013). Instandhouding Vaarpassen Schelde Milieuvergunningen terugstorten baggerspecie. LTV - Veiligheid en Toegankelijkheid. Baggeren en storten Achtergrondrapport A-31 Schramkowski, G.P.; Kaptein, S.; Vanlede, J.; Mostaert, F. (2019). Integraal plan Boven-Zeeschelde: Sub report 12 – Proof of uniqueness of the sediment solution for offline Telemac simulations. Version 1.0. FHR Reports, 13_131_12. Flanders Hydraulics Research: Antwerp, Belgium

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Smolders, S.; Maximova, T.; Vanlede, J.; Plancke, Y.; Verwaest, T.; Mostaert, F. (2016). Integraal Plan Bovenzeeschelde: Subreport 1 – SCALDIS: a 3D Hydrodynamic Model for the Scheldt Estuary. Version 5.0. WL Rapporten, 13_131. Flanders Hydraulics Research: Antwerp, Belgium. Smolders, S.; Maximova, T.; Vandenbruwaene, W.; Coen, L.; Vanlede, J.; Verwaest, T.; Mostaert, F. (2017). Integraal Plan Bovenzeeschelde: Deelrapport 5 – Scaldis 2050. Version 4.0. FHR Reports, 13_131_5. Flanders Hydraulics Research: Antwerp. Smolders, S.; De Maerschalck, B.; Plancke, Y.; Mostaert, F. (2019). Slibmodel Schelde-estuarium: Deelrapport 5 – Verbeteringen Scaldis slib met betrekking tot ETM karakteristieken. Versie 0.1. WL Rapporten, 17_088_5. Waterbouwkundig Laboratorium: Antwerpen.(in preparation) Soulsby, R. L., & Clarke, S. (2005). Bed shear-stress under combined waves and currents on smooth and rough beds (TR 137). Stevens M., Van den Neucker T., Gelaude E., Baeyens R., Jacobs Y., Mouton A., Buysse D. & Coeck J. (2011). Onderzoek naar de trekvissoorten in het Schelde-estuarium. Voortplantingsen opgroeihabitat van rivierprik en fint. Rapporten van het Instituut voor Natuur- en Bosonderzoek 2011 (14). Instituut voor Natuur- en Bosonderzoek, Brussel, Belgium. Tassi, P., 2016. SISYPHE User Manual – Version 7.2 (December 2016) Thant, S.; Plancke, Y. (2016). Sediment transport measurements in the Schelde-estuary: How do acoustic backscatter, optical transmission and direct sampling compare?, in: Erpicum, S. et al. (Ed.) Proceedings of the 4th IAHR Europe Congress, Liege, Belgium, 27-29 July 2016: Sustainable Hydraulics in the Era of Global Change. pp. 193-200 Uittenbogaard, R.E., 1998. Model for eddy diffusivity and viscosity related to sub-grid velocity and bed topography. Note, WL | Delft Hydraulics. Van Engeland T., T.J.S. Cox, K. Buis, S. Van Damme, P. Meire (2018). 1D Ecosystem model of the Schelde estuary: model calibration and validation. Report 018-R217. Van Hoestenberghe, T.; Ferket, B.; De Boeck, K.; Vanlierde, E.; Vanlede, J.; Verwaest, T.; Mostaert, F. (2014). Slibbalans Zeeschelde: Deelrapport 2 – Sediment load for the river Scheldt and its main tributaries (1972 -2009). Versie 5.0. WL Rapporten, 00_029. Waterbouwkundig Laboratorium & Antea Group. Antwerpen, België. van Kessel, T.; Vanlede, J.; Bruens, A. (2006). Development of a mud transport model for the Scheldt estuary in the framework of LTV: phases 1 and 2. Versie 1.0. Delft Hydraulics/Flanders Hydraulics Research: Delft. 79 + appendices pp. van Kessel, T.; Vanlede, J.; Kuijper, K.; de Kok, J. (2007). Further development and first application of a mud transport model for the Scheldt estuary: in the framework of LTV. Phase 2. Delft Hydraulics/Waterbouwkundig Laboratorium: Delft. 139 pp. van Kessel, T., J. Vanlede (2009). Impact of harbour basins on mud dynamics Scheldt estuary in the framework of LTV. Deltares & FHR report 1200253, Delft, The Netherlands. van Kessel, T.; Vanlede, J.; de Kok, J. (2011). Development of a mud transport model for the Scheldt estuary. Cont. Shelf Res. 31(10 SUPPL.). doi:10.1016/j.csr.2010.12.006 van Leussen, W., 1994. Estuarine Macroflocs: Their Role in Finegrained Sediment Transport. Universiteit van Utrecht, the Netherlands. 76

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

van Maren, D.S.; van der Werf, J.; van Holland, G.; Sas, M. (2010) Ontwikkelingszone Saeftinghe Hydraulisch, sedimentologisch en morfologisch onderzoek. Deelrapport 5: Calibratie Hydrodynamica Saeftinghedokmodel. Gemeentelijk Havenbedrijf Antwerpen, Belgium van Maren, D.S.; van Holland, G.; Sas, M. (2011) Ontwikkelingszone Saeftinghe - Hydraulisch, sedimentologisch en morfologisch onderzoek. Deelrapport 6: Calibratie Sediment Saeftinghedokmodel. Gemeentelijk Havenbedrijf Antwerpen, Belgium van Maren, D.S.; Winterwerp, J.C.; Vroom, J. (2015). Fine sediment transport into the hyper-turbid lower Ems River: the role of channel deepening and sediment-induced drag reduction. Ocean Dyn. ISBN 1023601508 65: 589–605. doi:10.1007/s10236-015-0821-2 Vandenbruwaene, W.; Vanlede, J.; Plancke, Y.; Verwaest, T.; Mostaert, F. (2016). Slibbalans Zeeschelde: Deelrapport 4 – Historische evolutie SSC. Versie 6.0. WL Rapporten, 00_029_4. Waterbouwkundig Laboratorium & Antea: Antwerpen. Vandenbruwaene, W.; Levy, Y.; Plancke, Y.; Vanlede, J.; Verwaest, T.; Mostaert, F. (2017). Integraal plan Boven-Zeeschelde: Deelrapport 8 – Sedimentbalans Zeeschelde, Rupel en Durme. Versie 4.0. WL Rapporten, 13_131_8. Waterbouwkundig Laboratorium: Antwerpen. Vanoverbeke J., Van Reyckegem G., Van Braeckel A. & Van den Bergh E.(2019a). Modelinstrumentarium voor het voorspellen van habitatgeschiktheid van de Zeeschelde voor fint (Alosa fallax) - Deelrapport voor het Integraal plan Boven-Zeeschelde. Rapporten van het Instituut voor Natuur- en Bosonderzoek 2019 (18). Instituut voor Natuur- en Bosonderzoek, Brussel. DOI: doi.org/10.21436/inbor.14517836 Vanoverbeke J., Van Reyckegem G., Van Braeckel A. & Van den Bergh E. (2019b). Modelinstrumentarium voor het voorspellen van overwinterende aantallen wintertaling (Anas crecca) in de Boven-Zeeschelde Deelrapport voor het Integraal plan Boven-Zeeschelde. Rapporten van het Instituut voor Natuur- en Bosonderzoek 2019 (15). Instituut voor Natuuren Bosonderzoek, Brussel. DOI: doi.org/10.21436/inbor.14517871 Wartel, Parker et al. (1999), acoustic reflection techniques and sediment analysis Winterwerp, J., 2002. On the flocculation and settling velocity of estuarine mud. Continental Shelf Research 22, 1339-1360. Winterwerp, J.C., W.G.M. van Kesteren, 2004, Introduction to the physics of cohesive sediments in the marine environment, Elsevier, Developments in Sedimentology, 56. Winterwerp, J. C., A.J. Manning, C. Martens, T. de Mulder, J. Vanlede, (2006), “A heuristic formula for turbulence-induced flocculation of cohesive sediment”, Estuarine, Coastal and Shelf Science, Vol 68, pp 195-207 Winterwerp, J.C. (2007). On the sedimentation rate of cohesive sediment. Proc. Mar. Sci. ISBN 9780444521637 209–226. doi:10.1016/S1568-2692(07)80014-3 WL | Delft Hydraulics (2007) New developments in the mud transport module of Delft3D. Report II: Implementation, sensitivity analysis, calibration and validation.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Appendix A: Software choice for mud modelling The cohesive sediment transport model was implemented in Sedi3D. Two other software platforms (Delwaq and Sysiphe) were also considered, but not chosen in the end. This appendix argues why the other platforms were not selected.

Delwaq The first plan to model mud in this project was to use the existing LTV-MUD model running in Delwaq (van Kessel et al., 2011). In TELEMAC-2D and TELEMAC-3D the option exist to generate output from the hydrodynamic simulations to serve as input for Delwaq. The files created by TELEMAC could then be used directly as input for the Delwaq user interface. Several problem were encountered trying this options. The most important issues are listed here below: Comments in the .cli file (resolved) The Delwaq input files were created from a run with the calibrated Scaldis 2013 3D model. A first problem occurred by not being able to visualize the model area in the WAQ GUI. If in TELEMAC the .cli or conlim-file has comments behind the comment sign at the end of each line, the WAQ GUI can’t read it properly. An example of such a comment at the end of the line is: 4 5 5 0.000 0.000 0.000 0.000 5 0.000 0.000 0.000

87161

4544 # Terneuzen (86972 - 87495)

to indicate that this boundary is the boundary at Terneuzen in the model. This problem was fixed by the developers so that the WAQ GUI could cope with such comments. Number of exchanges in the .hyd file (resolved) A second problem was that in the .hyd file for Delwaq the number of exchanges was not given by a number. Here the problem was that the variable to be written by TELEMAC did not have sufficient significant digits. The Scaldis 3D model has a large model domain and the Fortran code was not adapted to cope with numbers larger than 6 digits. In the subroutines gredelhyd_autop.f (line 343) and wrihyd.f (line 224) I6 was changed to I7 in the Fortran code and this solved this problem. The changes were also incorporated into the next release of TELEMAC (V7P2r0) (TELEMAC forum post #24471; http://www.opentelemac.org/index.php/kunena/21-telemac-3d/10280-possible-bug-in-coupling-withdelwaq#24471). Aggregation for Telemac not implemented (unresolved) The next problem occurred when the TELEMAC grid was opened with the DIDO, where you can aggregate meshes to increase the computation speed. Although a description is given in the Delwaq documentation on how to aggregate a TELEMAC mesh, the code itself could not execute this. Apparently, the code to do so was never developed. This means that it is not possible at this moment (January 2017) to aggregate a TELEMAC mesh. The amount of detail that is present in the TELEMAC mesh (Scaldis model) is not necessary to perform the Delwaq simulations. The smallest elements determine the global time step of the simulation and since grid aggregation is not possible this meant that the time step for the Delwaq simulations had to be set to five seconds.

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Another reason to set the time step to five seconds was the fact that for several areas in the model very small water volumes were present. One of the input files for Delwaq coming from TELEMAC is a file containing the exchanged volumes of water between the different elements. For some elements these volumes are almost zero (~1.0e-20). the implicit solvers of Delwaq (15/16 and 21/22) can’t handle such small volumes and the explicit solvers can’t handle the very low residence times. An experimental test with a new explicit solver that uses local smaller time steps per segment and time step to fulfill the CFL criterion ((Q*dt)/V < 1) managed to get through a small test period, but time steps of 0.05 seconds were necessary and this makes the simulation go even slower than reality. Even a small program that changed the minimum volumes (but also maintaining the global mass balance, by adding a small volume to all segments) could also not improve the situation. Parallelisation (unresolved) The Delwaq processes can run partially in parallel, but with restrictions and only on 1 node. So parallelization could not improve computation time. Conclusion In conclusion: A Delwaq model is possible on top of a hydrodynamic TELEMAC model, but the computation time is a bottleneck. The best performance was a simulation that was as fast as reality. This was not feasible to simulate mud dynamics and the Delwaq option was left.

SISYPHE After the setback with Delwaq, the idea was to keep the offline approach of Delwaq by using SISYPHE offline on top of a hydrodynamic simulation of TELEMAC-3D. SISYPHE is the existing sediment transport module of the Telemac software suite. SISYPHE itself is 2D and will only use the 2D result file of a 3D simulation. SISYPHE had the advantage that, although a small time step would be necessary for the current mesh resolution (the same restriction like Delwaq), it can run in parallel mode. First results showed performances of a speed up of 360 (one year simulated in one day). Problems arose when looking into detail to sediment transport in the upstream river reaches of the model. Advection schemes 13 or 14 in SISYPHE, to be used in combination with the presence of tidal flats, did not work. There is no advective sediment transport in the SISYPHE model, because advection scheme 14 that is used in the calibrated HD model does not write the necessary fluxes to the output file (fluxes are not a standard output parameter) to be used in an offline SISYPHE simulation. On a TELEMAC developers meeting (March 2017, Grenoble) we further learned that currently, nobody uses SISYPHE offline as we intended to and that the vision of the developers was to use SISYPHE in the future only for bed load transport and bed updating. Suspended transport will be removed from SISYPHE and included in TELEMAC-3D in a module called SEDI-3D. For the reasons stated above, it was decided to stop development of a mud transport model in SISYPHE, and move on to the SEDI-3D module. Error advection schemes following keywords: /---------------------------------------------------------------------/

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary /

SISYPHE V7P2 STEERING FILE

TEST CASE BERGENMEERSEN 3D

/---------------------------------------------------------------------/ /

AUTHORS: SVEN SMOLDERS

DATE : 03/03/2016

/ /

/---------------------------------------------------------------------/ / COMPUTER INFORMATIONS

/---------------------------------------------------------------------/ GEOMETRY FILE

: GEO_V20_2013_T_CRIT_SETL.SLF

BOUNDARY CONDITIONS FILE

: CONLIM_V17_SIS2.CLI

RESULTS FILE

: RES_SIS2_RUN01.SLF

FORTRAN FILE

: NOER6.F

HYDRODYNAMIC FILE

: MYLINK /R2D_V19_041_0.SLF

NUMBER OF PRIVATE ARRAYS

NAMES OF PRIVATE VARIABLES : TAU_CRIT_L1;SETTLING COMPUTATION CONTINUED

: NO

/PREVIOUS SEDIMENTOLOGICAL COMPUTATION FILE = PARALLEL PROCESSORS

:16

/---------------------------------------------------------------------/ / GENERAL INFORMATIONS

/---------------------------------------------------------------------/ TITLE : 'SUSPENDED SEDIMENT TRANSPORT IN SCHELDT ESTUARY' TIME STEP : 1 NUMBER OF TIME STEPS : 5 NUMBER OF SUB-ITERATIONS : 10 GRAPHIC PRINTOUT PERIOD : 6 /288 TIMES THE GRAPHIC PRINTOUT TIMESTEP IN THE HYDRODYNAMIC FILE VARIABLES FOR GRAPHIC PRINTOUTS : 'U,V,S,H,B,E,*ES,M,TOB,R,KS,TOB,MU,QS1,CS1,QSSUSP' LISTING PRINTOUT PERIOD : 6 STEADY CASE : NO TIDE PERIOD : 45000 STARTING TIME OF THE HYDROGRAM : -46000 NUMBER OF TIDES OR FLOODS : 20 MASS-BALANCE : YES CRITICAL EVOLUTION RATIO : 15000000000000000000.0 EFFECT OF WAVES : NO / FLOW-SEDIMENT INTERACTIONS /---------------------------------------------------------------------/ MIXED SEDIMENT: NO COHESIVE SEDIMENTS: YES NUMBER OF SIZE-CLASSES OF BED MATERIAL : 1 INITIAL FRACTION FOR PARTICULAR SIZE CLASS = 1 SEDIMENT DIAMETERS = 0.00006 SEDIMENT DENSITY = 2560 WATER VISCOSITY = 1.0E-6 WATER DENSITY = 1025

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary GRAVITY ACCELERATION = 9.81 LAW OF BOTTOM FRICTION FRICTION COEFFICIENT

:4 : 0.022

BED ROUGHNESS PREDICTION

: NO

BED ROUGHNESS PREDICTOR OPTION : 1 /

BED-LOAD TRANSPORT

BED LOAD : NO SLOPE EFFECT : YES FORMULA FOR SLOPE EFFECT : 1 =KOCH AND FLOKSTRA BETA : 1.3 / ONLY USED WITH FORMULA 1 = DEFAULT SECONDARY CURRENTS : YES /DEFAULT ONLY NECCESSARY BY COUPLING WITH 2D MODEL SECONDARY CURRENTS ALPHA COEFFICIENT : 0.75 /DEFAULT OPTION FOR THE TREATMENT OF NON ERODABLE BEDS :3 TIDAL FLATS : YES MINIMAL VALUE OF THE WATER HEIGHT : 0.01 OPTION FOR THE TREATMENT OF TIDAL FLATS : 2 /SUSPENDED LOAD TRANSPORT SUSPENSION : YES SETTLING VELOCITIES : 0.001 DIFFUSION : YES OPTION FOR THE DIFFUSION OF TRACER : 1 OPTION FOR THE DISPERSION : 2 TYPE OF ADVECTION : 4 / TO KEEP THE MASS LOSS DUE TO NUMERICS TO A MINIMUM WE NEED AN ADVECTION SCHEME / THAT IS MASS CONSERVATIE. SO WE NEED TO CHOOSE 4 OR 5 OR WITH TIDAL FLATS 14 SCHEME OPTION FOR ADVECTION:1 SOLVER FOR SUSPENSION : 3 SOLVER ACCURACY FOR SUSPENSION : 1.E-10 MASS CONCENTRATION : YES REFERENCE CONCENTRATION FORMULA : 1/ZYSDERMAN AND FREDSOE, EQUILIBRIUM FORMULA CORRECTION ON CONVECTION VELOCITY : YES / DEFAULT OPTION INITIAL SUSPENSION CONCENTRATIONS : 0/MASS CONCENTRATION EQUILIBRIUM INFLOW CONCENTRATION : NO /IF YES THE KEYWORD ABOVE WILL NOT BE USED TREATMENT OF FLUXES AT THE BOUNDARIES : 1;1;1;1;1;1;1;1;1 CONCENTRATION PER CLASS AT BOUNDARIES: 0;3.546E-5;3.714E-5;2.341E-5;2.799E-5;1.693E-5;1.564E-5;0;4.87E-6 / 1=TERNEUZEN (NO SEDIMENT) 2=MERELBEKE; 3=DENDER; 4=ZENNE / 5=DIJLE; 6=GR NETE; 7=KL NETE; 8=BATH; 9=ZEERAND / EROSION, DEPOSITION AND CONSOLIDATION

FORMULATION FOR DEPOSITION AND EROSION : 1 /1:KRONE ET PARTHENIADES PARTHENIADES CONSTANT : 1.0E-03 CRITICAL SHEAR VELOCITY FOR MUD DEPOSITION : 1000 MUD CONSOLIDATION = YES CONSOLIDATION MODEL = 1 NUMBER OF LAYERS OF THE CONSOLIDATION MODEL = 2 MUD CONCENTRATION PER LAYER : 80;200

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary CRITICAL EROSION SHEAR STRESS OF THE MUD = 0.1;0.5 /;0.19862961;1.2 / NUMERICAL SCHEMES

FINITE VOLUMES : NO MASS-LUMPING : YES / MASS MATRIX IS CONDENSED TO ITS DIAGONAL TO SPEED UP AND STABILIZE THE COMPUTATION, BUT RESULTS ARE SMOOTHED TETA SUSPENSION : 0.5 MAXIMUM NUMBER OF ITERATIONS FOR SOLVER FOR SUSPENSION : 50 /DEFAULT = 50 PRECONDITIONING FOR SUSPENSION : 2 /DIAGONAL PRECONDITIONING SOLVER ACCURACY FOR SUSPENSION : 1.E-8 MATRIX STORAGE : 3 /DEFAULT =1 MATRIX-VECTOR PRODUCT : 1 /DEFAULT ZERO : 1.E-10 TETA : 0.5 / &ETA

gave the following error: … ================================================================================ ITERATION

0 TIME: 0.0000 S

INITIAL QUANTITY IN SUSPENSION FOR CLASS 1 :

0.000000

CVTRVF: 200 SUB-ITERATIONS REQUIRED FOR THE DISTRIBUTIVE SCHEME. DECREASE THE TIME-STEP PLANTE: PROGRAM STOPPED AFTER AN ERROR RETURNING EXIT CODE:

RUNCODE::MAIN:

|RUNCODE: FAIL TO RUN |/OPT/OPENMPI/INTEL/64/1.6.3/BIN/MPIEXEC

--BIND-TO-SOCKET

/PROJECTS/13_131_BEVAARBOZS/SISYPHE/RUN_002/SIS_SUS2_SCALDIS_001.CAS_2017-05-04-11H13MIN28S

-WDIR -N

/PROJECTS/13_131_BEVAARBOZS/SISYPHE/RUN_002/SIS_SUS2_SCALDIS_001.CAS_2017-05-04-11H13MIN28S/OUT_NOER6 |~~~~~~~~~~~~~~~~~~ |1

subroutine CVTRVF.f is located in the BIEF. It mentions in the following part of code: IF(NIT.GE.NITMAX.AND.OPTADV.NE.4) THEN  IF NUMBER OF ITERATIONS > MAX NUMBER OF ITERATIONS AND SCHEME OPTION FOR ADVECTION IS NOT EQUAL TO 4

IF(LNG.EQ.1) WRITE(LU,900) NIT IF(LNG.EQ.2) WRITE(LU,901) NIT 900

FORMAT(1X,'CVTRVF : ',1I6,' SOUS-ITERATIONS DEMANDEES POUR LE'

& ,/,1X, '

Final version

SCHEMA DISTRIBUTIF, REDUIRE LE PAS DE TEMPS')

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary 901

FORMAT(1X,'CVTRVF: ',1I6,' SUB-ITERATIONS REQUIRED FOR THE'

& ,/,1X, '

DISTRIBUTIVE SCHEME. DECREASE THE TIME-STEP')

CALL PLANTE(1) STOP

This means that there is a problem with the keyword and our choice for the advection scheme. It iterates too much. Message Subject : advection schemes 13 and 14 don't work Category : SISYPHE Posted by : jmhervouet URL : http://www.opentelemac.org/index.php/kunena/17-sisyphe/10738-advection-schemes-13-and-14-don-t-work/26339 Message : ----Hello Sven, Sure it is the reason. To run the NERD scheme (13-14) you need to have the fluxes between points computed by Telemac-2D or 3D. Generally speaking to have a conservative advection when solving the non-conservative form of the advection equation you need to have the discrete continuity equation available (old and new depths, fluxes between points). This is not possible in stand-alone mode, especially if you do not keep the same time-step. Furthermore, in standalone mode Sisyphe updates the velocities with somewhat wrong assumptions (e.g. considering that the product h*u will remain constant). This latter approximation may produce strange patterns of sedimentation, which are totally artificial. This is why we think that we should stop the standalone mode and do not recommend it. This does not explain why you seem to have good results with scheme 4. With this latter scheme, which does not work with tidal flats I programmed an option in the case where the velocity field does not obey the continuity equation. This saves mass conservation but spoils the monotonicity. This option is also used when there is a velocity correction in Sisyphe, to take into account the fact that the sediment "sees" the velocity near the bottom. With best regards, Jean-Michel Hervouet

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Appendix B: treatment of tidal flats in Telemac 3D Variables involved in treatment of tidal flats The treatment of tidal flats and variables that set minimum depths can have a significant influence on the results obtained in the simulations. Because Telemac is open source, a user can check inside the code where, when and how the settings for treatment of tidal flats are being used in Telemac 3D. First a list of all the relevant keywords is given, with the names of these variables of how they are used in the code (= MNEMO). Then all relevant pieces of code are listed with some explanation on what is going on. Table 14 – list of all relevant variables for treatment of tidal flats in Telemac 3D

keyword TIDAL FLATS OPTION FOR THE TREATMENT OF TIDAL FLATS

Value model YES 1

TREATMENT OF NEGATIVE DEPTHS

MINIMAL VALUE FOR DEPTH THRESHOLD FOR SEDIMENT FLUX CORRECTION ON TIDAL FLATS TREATMENT ON TIDAL FLATS FOR VELOCITIES

0.1 m

TREATMENT ON TIDAL FLATS FOR TRACERS

THRESHOLD FOR VISCOSITY CORRECTION ON TIDAL FLATS

0.2

0.1 m

Explanation Activate if tidal flats are present in the model Algorithms of treating tidal flats: 1: equations solved everywhere with correction on tidal flats, 2: dry elements frozen It is only used if the keyword OPTION FOR THE TREATMENT OF TIDAL FLATS = 1. Possible choices are: 0: no treatment, 1: smoothing, 2: flux control the minimum water depth, below which tidal flats are considered as dry. Below this limiting depth, all sediment erosion rates are set to zero.

Treatment of tidal flats at the diffusion step for tracers. 0: forced to zero, 1: value before masked Treatment of tidal flats at the diffusion step for tracers. 0: forced to zero, 1: value before masked Below the threshold, viscosity will be progressively cancelled. See subroutine VISCLIP

MNEMO

default value

BANDEC OPTBAN

YES 1

OPT_HNEG

HMIN

-1000

HSED

0.2

TRBAVI

TRBATA

HLIM

0.2

With every new simulation the keywords of the steering file are read (using lecdon_telemac3d.f) and parameters get a value. After that, the main subroutine telemac3D.f is called to start the calculations. Code out of both subroutines is discussed here:

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

LECDON_TELEMAC3D.f - line 1643: ! SPECIAL TREATMENT IF PARALLELISM IF(NCSIZE.GT.1.AND.BANDEC.AND.OPTBAN.EQ.2) THEN OPTBAN=1 IF (LNG.EQ.1) WRITE(LU,121) IF (LNG.EQ.2) WRITE(LU,122) 121 FORMAT(/,'ATTENTION: VOUS AVEZ CHOISI LE MODE PARALLELE,', & /,'========= L''OPTION DE TRAITEMENT DES BANCS', & /,' DECOUVRANTS EST MISE A 1') 122 FORMAT(/,'ATTENTION: YOU HAVE CHOSEN PARALLEL MODE,', & /,'========= THE TIDAL FLATS TREATMENT IS SET TO 1') ENDIF This piece of code sets the OPTION FOR THE TREATMENT OF TIDAL FLATS back to 1, even when chosen option 2 by the user, if used in parallel. For Scaldis this has no influence because option 1 was chosen - line 1853: IF(.NOT.BANDEC) THEN OPTBAN = 0 OPT_HNEG = 0 ENDIF IF(OPTBAN.EQ.2) THEN MSK = .TRUE. ! WITH A NEGATIVE HMIN, MASKBD WILL FAIL HMIN = MAX(HMIN,0.D0) ELSEIF(MSKUSE) THEN MSK = .TRUE. ELSE ! NOTE JMH : MASKING BY THE USER DOES NOT APPEAR TO HAVE BEEN ENVISAGED MSK = .FALSE. ENDIF A new variable is introduced: MSK. For the choices made in Scaldis MSK=.FALSE. - line 1690: ! WITH SOME TIDAL FLATS VERSIONS OF DISTRIBUTIVE ADVECTION SCHEMES ! POSITIVE DEPTHS MUST BE TREATED WITH OPTION 2 IF(BANDEC.AND.OPTBAN.EQ.1.AND.OPT_HNEG.NE.2) THEN This piece of code just ends the simulation and gives an error message if the TREATMENT OF NEGATIVE DEPTHS is not equal to 2 when tidal flats are present and OPTION FOR TREATMENT OF TIDAL FLATS = 1. The code does not change the variable itself.

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

TELEMAC3D.f To run a 3D simulation the subroutine Telemac3D.f is the main script to follow. In this section this subroutine is read and pieces of code handling tidal flats are taken out and are discussed. Calls to other subroutines will be noted and these will be handled in the next sections. - line 744: ! NOTE : HMIN = -1000.0 IN DICTIONARY BUT HMIN IS AT LEAST 0.0 IF OPTBAN=2 IF(OPTBAN.EQ.2) THEN CALL CLIP (H, HMIN, .TRUE., 1.D6, .FALSE., 0) ENDIF This part of code has no influence for Scaldis - line 777: CALL VERMOY(U2D%R,V2D%R,U%R,V%R,2,Z, & T3_01%R,T3_02%R,T3_03%R,1,NPLAN,NPOIN2,NPLAN,OPTBAN) The subroutine VERMOY is called to calculated the average of a 3D variable on the vertical. The OPTBAN variable chooses the way this averaging is done. If OPTBAN = 1 the averaging is done different compared to OPTBAN = 0 or 2. - line 858: CALL FSGRAD(GRADZS,ZFLATS,Z(NPOIN3-NPOIN2+1:NPOIN3), & ZF,IELM2H,MESH2D,MSK,MASKEL, & UNSV2D,T2_01,NPOIN2,OPTBAN,SVIDE) The subroutine FSGRAD calculates the free surface gradient taking the treatment of tidal flats into account. The difference between treatments is that this subroutine needs to know if it need to take the elements of tidal flats into account or not. - line 922: CALL TFOND(AUBORF%R,CF%R,U2D%R,V2D%R,U%R,V%R,W%R,KARMAN, & LISRUF,DNUVIV,Z,NPOIN2,KFROT,RUGOF%R,UETCAR%R, & NONHYD,OPTBAN,HN%R,GRAV,IPBOT%I,NPLAN) The subroutine TFOND calculates AUBOR, the coefficient for the log law at the bottom. If OPTBAN = 1 in this subroutine a minimum velocity (U*)² = g*H is imposed (if water depth is smaller than 1 mm) to have some friction on the tidal flats to oppose the free surface gradient.  TFOND.f ! FINAL COMPUTATION OF AUBOR = - (U*)**2 / U(BOTTOM) ! ON TIDAL FLATS A MINIMUM VELOCITY IS IMPOSED (U*)**2=GH ! TO OPPOSE SOME FRICTION TO FREE SURFACE GRADIENTS ! IF(OPTBAN.EQ.1) THEN DO N=1,NPOIN IF(HN(N).LT.1.D-3) THEN UETCAR(N)=MAX(GRAV*MAX(HN(N),1.D-7),UETCAR(N))

Final version

WL2020R13_131_6

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

ENDIF ENDDO ENDIF

- line 986 - line 986: IF(OPTBAN.EQ.1) THEN CALL VISCLIP(VISCVI,VISCTA,H,NPLAN,NPOIN3,NPOIN2,NTRAC,HLIM) ENDIF At the end of the section for calculating the viscosities the subroutine VISCLIP is called if OPTBAN = 1. This subroutine is called to limit turbulent viscosity on tidal flats.  VISCLIP.f DO I=1,NPOIN2 IF(H%R(I).LT.HLIM) THEN COR=(MAX(H%R(I),0.D0)/HLIM)**2 DO IPLAN=1,NPLAN VISCVI%ADR(1)%P%R(I+(IPLAN-1)*NPOIN2)= & VISCVI%ADR(1)%P%R(I+(IPLAN-1)*NPOIN2)*COR VISCVI%ADR(2)%P%R(I+(IPLAN-1)*NPOIN2)= & VISCVI%ADR(2)%P%R(I+(IPLAN-1)*NPOIN2)*COR VISCVI%ADR(3)%P%R(I+(IPLAN-1)*NPOIN2)= & VISCVI%ADR(3)%P%R(I+(IPLAN-1)*NPOIN2)*COR ENDDO ENDIF ENDDO If the water depth (H%R) is smaller than HLIM (THRESHOLD FOR VISCOSITY CORRECTION ON TIDAL FLATS; default = 0.2 m), a correction term COR for the turbulent viscosity will be calculated and applied. The smaller the water depth the smaller this correction term becomes. The correction term is squared and a squared number smaller than 1 becomes even smaller. So the turbulent viscosity on tidal flats becomes progressively smaller. the same is done for the viscosity for the tracers. - line 1451: & & & &

CALL VITCHU(WCHU,WCHU0,TURBWC,U,V,W,H,RUGOF,LISRUF, TURBA,TURBB,T3_01,T3_02,T3_03,SVIDE,MESH3D, IELM3,NPOIN2,NPOIN3,NPLAN,NTRAC,MSK,MASKEL, UETCAR,TA,HN,HSED,FLOC, FLOC_TYPE, HINDER,HIND_TYPE,CGEL,CINI)

When sediments are active VITCHU is called to calculate the settling velocity. HSED (THRESHOLD FOR SEDIMENT FLUX CORRECTION ON TIDAL FLATS) is used only when flocculation is used and when the flocculation type = 2. Then the subroutine SOULSBYFLOC3D is called and HSED is used. In VITCHU HSED is called HMIN, which can be confusing for the user.  VITCHU.f CALL SOULSBYFLOC3D(WCHU,TRAV1%R,MESH3D,NPOIN2, & NPOIN3,NPLAN,HN,HMIN,UETCAR%R)

A10

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

- line 1490: IF(LT.GT.1) THEN CALL PLANE_BOTTOM(IPBOT%I,ZPROP%R,NPOIN2,NPLAN,SIGMAG,OPTBAN) ENDIF Call to subroutine PLANE_BOTTOM to find variable IPBOT (= plane number of last crushed plane). For every 2D point, this subroutine finds the last plane with no normal height above, i.e. delta(z) equal to zero and if no problem arises, IPBOT = 0. If tidal flats are present, IPBOT = NPLAN-1 (=total number of planes -1). So plane IPBOT+1 always exists and has the first free point, unless there is no depth.  PLANE BOTTOM.f IF(SIGMAG.OR.OPTBAN.EQ.1) THEN ! DO IPOIN2=1,NPOIN2 IPBOT(IPOIN2)=0 DO IPLAN=1,NPLAN-1 IF(Z(IPOIN2,IPLAN+1)-Z(IPOIN2,IPLAN).LT.1.D-4) THEN IPBOT(IPOIN2)=IPLAN ENDIF ENDDO ENDDO ELSE DO IPOIN2=1,NPOIN2 IPBOT(IPOIN2)=0 ENDDO ENDIF - line 1721: CALL CLSEDI & (ATABOF%ADR(NTRAC),…,HSED) Call to CLSEDI to express the boundary conditions for sediment at the bottom and the surface. HSED (THRESHOLD FOR SEDIMENT FLUX CORRECTION ON TIDAL FLATS; default = 0.2 m) is used inside CLSEDI as HMIN. Within this subroutine another subroutine is called using HMIN: ERODC and FLUSED:  ERODC.f IF(HN(IPOIN).LT.HMIN) THEN QERODE=0.D0 In the subroutine ERODC the erosion is modelled (within a multi-layer consolidation model). The user can choose the law defining the critical erosion velocity as a function of the concentration. By default the Partheniades formulation is used. If the water depth becomes lower than HSED (HMIN letters used in code), then the erosion flux becomes zero.  FLUSED.f ! COMPUTES THE DEPOSITION PROBABILITY IF(SEDCO) THEN ! COHESIVE SEDIMENT (Here FLUDPT >0) IF(SIGMAG.OR.OPTBAN.EQ.1) THEN DO I=1,NPOIN2 IF(IPBOT%I(I).NE.NPLAN-1) THEN ! DEPOSITION ON THE FIRST FREE PLANE WITH LOCAL VELOCITY

Final version

WL2020R13_131_6

A11

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

…

I3D=I+IPBOT%I(I)*NPOIN2 FLUDPT(I) = WC(I3D)*MAX(1.D0-TOB(I)/MAX(TOCD,1.D-6),0.D0) ELSE ! TIDAL FLAT FLUDPT(I) = 0.D0 ENDIF ENDDO ELSE DO I=1,NPOIN2 FLUDPT(I) = WC(I)*MAX(1.D0-(TOB(I)/MAX(TOCD,1.D-6)),0.D0) ENDDO ENDIF ENDIF ! COMPUTATION OF THE TRACER FLUX ON THE BOTTOM IF(SETDEP.EQ.1) THEN ! (SETDEP is choice of advection diffusion scheme for settling velocity; ! default = 0) ! USING HMIN TO CLIP EROSION (DIFFERENT FROM USING IPBOT) DO I=1,NPOIN2 IF(HN(I).LE.HMIN) THEN FLUER(I) = 0.D0 ENDIF ENDDO ELSEIF(SIGMAG.OR.OPTBAN.EQ.1) THEN DO I=1,NPOIN2 ATABOF(I) = 0.D0 BTABOF(I) = 0.D0 IF(LITABF(I).EQ.KLOG) THEN ! NO EROSION AND DEPOSITION ON TIDAL FLATS IF(IPBOT%I(I).NE.NPLAN-1) THEN ATABOF(I) = -FLUDPT(I) BTABOF(I) = FLUER(I) ENDIF ENDIF ENDDO

HMIN is in the subroutine FLUSED only used in a very specific case when the keyword ADVECTIONDIFFUSION SCHEME WITH SETTLING VELOCITY (default = 0) is set to 1. In other cases with tidal flats in case OPTBAN = 1, the erosion and deposition flux on tidal flats are set to zero. Be aware that this is only when the 3D planes are crushed together (IPBOT=0) and there is no or almost no water on the tidal flats. In other situation erosion and deposition are calculated as in all other cases. - line 1852: ! CLIPS HBOR IF(OPTBAN.EQ.2) THEN CALL CLIP(HBOR,HMIN,.TRUE.,1.D6,.FALSE.,0) ENDIF HBOR is the prescribed depth on lateral boundaries. It is not used in Scaldis.

A12

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

- line 2078: call to CVDF3D.f which solves the advection-diffusion step by call the chosen scheme. If OPTBAN = 1 tidal flats are solved in a piece wise linear way. In CVDF3D.f for solving the 3D diffusion and SUPG advection steps is referred to DIFF3D.f. At the end of this subroutine the value of tidal flats that where masked or if required so, was set to the value of the previous time step: IF(MSK.AND.TRBAF.EQ.1) THEN DO IPOIN3 = 1,NPOIN3 IF(MASKPT%R(IPOIN3).LT.0.5D0) FD%R(IPOIN3) = FN%R(IPOIN3) ENDDO In this piece of code TRBAF represents TRBATA (TREATMENT ON TIDAL FLATS FOR TRACERS) - line 2236: IF(OPT_HNEG.EQ.2) THEN CALL FLUX3DLIM(FLODEL%R,FLULIM%R,NPLAN,MESH2D%NSEG,NPOIN2,1) ENDIF The subroutine FLUX3DLIM limits the 3D horizontal edge by edge fluxes on points.

Final version

WL2020R13_131_6

A13

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Appendix C: Formulas for dredging amounts This paragraph is taken from Appendix A of Dams et al. (2016). In literature different units are being used to report dredging amounts. Some texts mention formulas to convert in between units. Because of differences in notation, these formulas are sometimes difficult to compare. This appendix contributes to that discussion by deriving conversion formulas from known basic definitions, and by making all assumptions explicit.

Notation Symbol

𝒎𝒎𝒈𝒈 𝒎𝒎𝒃𝒃 𝑽𝑽𝒃𝒃 𝑽𝑽𝒈𝒈 𝝆𝝆𝒃𝒃 𝝆𝝆𝒘𝒘 𝝆𝝆𝒈𝒈 𝑪𝑪

Description Mass of dry fraction, mass of grains Mass of mixture; bulk mass Volume of Mixture water sediment (Bulk Volume) Volume of dry fraction; volume of grains Bulk density Density of water Grain density Volume concentration

Unit kg, TDS kg m³ m³ kg/m³ kg/m³ kg/m³ -

Volume concentration Volume concentration is a powerful concept for conversions. It is defined as the ratio [-] of volume of grains to the volume of the mixture.

One can link this to densities as:

𝐶𝐶 = 𝐶𝐶 =

Reduced Volume V’

𝑉𝑉𝑔𝑔 𝑉𝑉𝑏𝑏

𝜌𝜌𝑏𝑏 − 𝜌𝜌𝑤𝑤 𝜌𝜌𝑔𝑔 − 𝜌𝜌𝑤𝑤

[definition in Dutch, taken from Ministerie van de Vlaamse Gemeenschap (1991)] “Onderhoudsbaggerwerken op de drempels in de vaargeul Er werd onderscheid gemaakt tussen specie met een densiteit groter of gelijk aan 1.6 en specie met een lagere densiteit.[…] In de tabel wordt bij een densiteit groter of gelijk aan 1.6 het volume aangegeven dat rechtstreeks in de baggerschepen werd opgemeten. Bij een densiteit kleiner dan 1.6, dit is het geval van slibhoudende specie of zelfs zuiver slib, wordt het volume aangegeven dat eenzelfde hoeveelheid vaste specie zal aannemen bij een densiteit gelijk aan 2.

A14

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Dit komt er in feite op neer dat het volume slibhoudende specie herleid wordt naar het volume dat eenzelfde gewichtshoeveelheid met water verzadigd zand zou innnemen.” This is the definition of the reduced volume V’. For reasons of consistency, V’ is noted in full as 𝑉𝑉′𝑏𝑏 in the formulas below. The formula for V’ is derived based on the assumption of invariance of 𝑚𝑚𝑔𝑔 and 𝑉𝑉𝑔𝑔 .

𝑚𝑚𝑔𝑔 is invariant because of the last paragraph in the definition from Ministerie van de Vlaamse Gemeenschap (1991). Invariance of 𝑉𝑉𝑔𝑔 follows, because of the fact that 𝜌𝜌𝑔𝑔 doesn’t change. Based on the invariance of 𝑉𝑉𝑔𝑔 under the conversion from 𝑉𝑉𝑏𝑏 to 𝑉𝑉′𝑏𝑏 , one can work out the conversion relation. Normal

Reduced (for 𝝆𝝆𝒃𝒃 <1.6) 𝑉𝑉𝑔𝑔 = 𝐶𝐶′ ∙ 𝑉𝑉′𝑏𝑏 𝜌𝜌′𝑏𝑏 − 𝜌𝜌𝑤𝑤 𝐶𝐶′ = 𝜌𝜌𝑔𝑔 − 𝜌𝜌𝑤𝑤 𝜌𝜌′𝑏𝑏 = 2 kg/m³ (definition) 𝐶𝐶 ∙ 𝑉𝑉𝑏𝑏 = 𝐶𝐶′ ∙ 𝑉𝑉′𝑏𝑏 𝜌𝜌𝑏𝑏 − 𝜌𝜌𝑤𝑤 𝑉𝑉′𝑏𝑏 = 𝑉𝑉𝑏𝑏 2 − 𝜌𝜌𝑤𝑤

𝑉𝑉𝑔𝑔 = 𝐶𝐶 ∙ 𝑉𝑉𝑏𝑏 𝜌𝜌𝑏𝑏 − 𝜌𝜌𝑤𝑤 𝐶𝐶 = 𝜌𝜌𝑔𝑔 − 𝜌𝜌𝑤𝑤

In case we assume 𝜌𝜌𝑤𝑤 = 1, the conversion formula above to calculate 𝑉𝑉′𝑏𝑏 from 𝑉𝑉𝑏𝑏 simplifies to: 𝑉𝑉′𝑏𝑏 = 𝑚𝑚𝑏𝑏 − 𝑉𝑉𝑏𝑏

Both 𝑚𝑚𝑏𝑏 and 𝑉𝑉𝑏𝑏 are readily available from automatic measurements on a dredging vessel. 𝑉𝑉𝑏𝑏 is derived from the level the dredged material in the hull, combined with the geometry of the hull. 𝑚𝑚𝑏𝑏 is derived from the sinking-in (immersion) of the dredging vessel.

Conversion of 𝑽𝑽′𝒃𝒃 [m³] to 𝒎𝒎𝒈𝒈 [TDS]

Combining the basic relation 𝑚𝑚𝑔𝑔 = 𝑉𝑉𝑔𝑔 ∙ 𝜌𝜌𝑔𝑔 with 𝑉𝑉𝑔𝑔 = 𝐶𝐶′ ∙ 𝑉𝑉′𝑏𝑏 (derived above), one can easily work out that

the

relation

under

reduced

volume

𝑚𝑚𝑔𝑔 = 𝜌𝜌𝑔𝑔 ∙ 𝐶𝐶′ ∙ 𝑉𝑉′𝑏𝑏

The table below gives this conversion factor under typical values: 2.65 2.65

𝝆𝝆𝒈𝒈

1 1.025

𝝆𝝆𝒘𝒘

Conversion factor 𝝆𝝆𝒈𝒈 ∙ 𝑪𝑪′ 1.606 1.59

Note that because the assumption 𝜌𝜌𝑤𝑤 = 1 is typically used in the derivation of 𝑉𝑉′𝑏𝑏 from 𝑚𝑚𝑏𝑏 and 𝑉𝑉𝑏𝑏 (derivation above), it is recommended to use the conversion factor 1.606 to convert reduced volume to mass.

Final version

WL2020R13_131_6

A15

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary

Appendix D: Steering file calibrated model TELEMAC3D V7P2r1 steering file calibrated mud model: /---------------------------------------------------------------------/ /

INPUT-OUTPUT, FILES

/---------------------------------------------------------------------/ FORTRAN FILE ='Scaldis_SEDI3D_CAL007_V7P2r1.f' LIQUID BOUNDARIES FILE

='Up_QN_2013_sedi3d.txt'

GEOMETRY FILE

='geo_v23_2013_sedi.slf'

BOUNDARY CONDITIONS FILE

='conlim_v23_sedi.cli'

2D RESULT FILE

='r2D_CAL_007_sedi3d.slf'

3D RESULT FILE

='r3D_CAL_007_sedi3d.slf'

FORMATTED DATA FILE 2

='Seaboundary_QN_2013.txt'

CULVERTS DATA FILE

='Tubes_2013_C.txt'

NUMBER OF CULVERTS

=35

OPTION FOR CULVERTS

TYPE OF SOURCES

MAXIMUM NUMBER OF SOURCES

=150

/----------------------------------------------------------------------/ /

RESTART FILE

/----------------------------------------------------------------------/ PREVIOUS COMPUTATION FILE

='HOTSTART_CAL_007'

PREVIOUS COMPUTATION FILE FORMAT = 'SERAFIND' COMPUTATION CONTINUED

=YES

INITIAL TIME SET TO ZERO

=NO

/---------------------------------------------------------------------/ /

INPUT-OUTPUT,TIME STEP, GRAPHICS AND LISTING

/---------------------------------------------------------------------/ TIME STEP

=4.0

NUMBER OF TIME STEPS

=864000

PARALLEL PROCESSORS

=160

GRAPHIC PRINTOUT PERIOD

=900

LISTING PRINTOUT PERIOD

=150

/40 days

VARIABLES FOR 2D GRAPHIC PRINTOUTS = 'U,V,S,H,EF,DF,US,HD,TA*' VARIABLES FOR 3D GRAPHIC PRINTOUTS = 'Z,U,V,W,TA*' ORIGINAL DATE OF TIME

= 2013;07;31

ORIGINAL HOUR OF TIME

= 22;20;00

MASS-BALANCE

=YES

INFORMATION ABOUT MASS-BALANCE FOR EACH LISTING PRINTOUT = YES /---------------------------------------------------------------------/ /

FRICTION

/---------------------------------------------------------------------/ FRICTION COEFFICIENT FOR THE BOTTOM

=0.02

LAW OF BOTTOM FRICTION

A16

/ Manning

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary LAW OF FRICTION ON LATERAL BOUNDARIES

/=Nikuradse Law

FRICTION COEFFICIENT FOR LATERAL SOLID BOUNDARIES

=0.054848

/---------------------------------------------------------------------/ /

EQUATIONS, BOUNDARY CONDITIONS

/---------------------------------------------------------------------/ VELOCITY PROFILES

=1;1;1;1;1;1;1;1;1

PRESCRIBED FLOWRATES

=23;34.7;11.1;15.92;34.6;8.3;10.38;35;0

PRESCRIBED ELEVATIONS

=1;1;1;1;1;1;1;1;1

OPTION FOR LIQUID BOUNDARIES =1;1;1;1;1;1;1;1;1 TREATMENT OF FLUXES AT THE BOUNDARIES =1;1;1;1;1;1;1;1;1 /---------------------------------------------------------------------/ /

EQUATIONS, INITIAL CONDITIONS

/---------------------------------------------------------------------/ NUMBER OF HORIZONTAL LEVELS

= 5

MESH TRANSFORMATION

= 2 /sigma transformation default=1

INITIAL ELEVATION

INITIAL CONDITIONS

='CONSTANT ELEVATION'

/---------------------------------------------------------------------/ /

NUMERICAL PARAMETERS

/---------------------------------------------------------------------/ MATRIX STORAGE

TREATMENT OF NEGATIVE DEPTHS

MASS-LUMPING FOR DEPTH

MASS-LUMPING FOR VELOCITIES

MASS-LUMPING FOR DIFFUSION

/---------------------------------------------------------------------/ /

TURBULENCE

/---------------------------------------------------------------------/ VERTICAL TURBULENCE MODEL

COEFFICIENT FOR VERTICAL DIFFUSION OF VELOCITIES

=1.E-2

MIXING LENGTH MODEL

HORIZONTAL TURBULENCE MODEL

COEFFICIENT FOR HORIZONTAL DIFFUSION OF VELOCITIES

=1.E-2

PRECONDITIONING FOR DIFFUSION OF K-EPSILON

ACCURACY FOR DIFFUSION OF VELOCITIES

= 1.E-5

MAXIMUM NUMBER OF ITERATIONS FOR DIFFUSION OF VELOCITIES = 201 CORIOLIS

=YES

CORIOLIS COEFFICIENT

=1.13522E-04

/----------------------------------------------------------------------/ /

ADVECTION

/----------------------------------------------------------------------/ ADVECTION STEP

=YES

SCHEME FOR ADVECTION OF VELOCITIES

SCHEME FOR ADVECTION OF DEPTH

SCHEME FOR ADVECTION OF TRACERS

=13

NUMBER OF SUB ITERATIONS FOR NON LINEARITIES

Final version

WL2020R13_131_6

A17

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary SOLVER FOR VERTICAL VELOCITY

/----------------------------------------------------------------------/ /

DIFFUSION

/----------------------------------------------------------------------/ SCHEME FOR DIFFUSION OF VELOCITIES

/default implicit (0 value cancels diffusion)

SCHEME FOR DIFFUSION OF TRACERS

SCHEME FOR DIFFUSION OF K-EPSILON

IMPLICITATION FOR DIFFUSION

=1. /default

SOLVER FOR DIFFUSION OF VELOCITIES

OPTION OF SOLVER FOR DIFFUSION OF VELOCITIES

/---------------------------------------------------------------------/ /

PROPAGATION

/---------------------------------------------------------------------/ PROPAGATION STEP

=YES

ACCURACY FOR PROPAGATION

=1.E-5

MAXIMUM NUMBER OF ITERATIONS FOR PROPAGATION

=500

SOLVER FOR PROPAGATION

FREE SURFACE GRADIENT COMPATIBILITY

=0.9

IMPLICITATION FOR DEPTH

=0.55

IMPLICITATION FOR VELOCITIES

INITIAL GUESS FOR DEPTH

OPTION OF SOLVER FOR PROPAGATION

/----------------------------------------------------------------------/ /

TRACERS

/----------------------------------------------------------------------/ NUMBER OF TRACERS = 2 NAMES OF TRACERS = 'SALINITY

PSU

';'MUD

G/L

INITIAL VALUES OF TRACERS = 0.0;0.5 PRESCRIBED TRACERS VALUES = 0.1;0.0;0.1;0.094;0.1;0.0984; 0.1;0.062;0.1;0.0742;0.1;0.045;0.1;0.042;0.1;0.0;0.1;0.012 DENSITY LAW

COEFFICIENT FOR VERTICAL DIFFUSION OF TRACERS

=1.E-6;1.E-6

COEFFICIENT FOR HORIZONTAL DIFFUSION OF TRACERS =1.E-6;1.E-6 PRECONDITIONING FOR DIFFUSION OF TRACERS

DAMPING FUNCTION

/----------------------------------------------------------------------/ / SEDIMENT TRANSPORT

/----------------------------------------------------------------------/ SEDIMENT = YES COHESIVE SEDIMENT = YES MEAN DIAMETER OF THE SEDIMENT = 0.00005 DENSITY OF THE SEDIMENT = 2650.0 CONSTANT SEDIMENT SETTLING VELOCITY = 0.5E-3 HINDERED SETTLING = NO BED LAYERS THICKNESS = 0.0

A18

WL2020R13_131_6

Final version

Integraal plan Boven-Zeeschelde - Sub report 6 – Scaldis Mud: a Mud Transport model for the Scheldt Estuary NUMBER OF SEDIMENT BED LAYERS = 1 INITIAL THICKNESS OF SEDIMENT LAYERS = 0.0 MUD CONCENTRATIONS PER LAYER = 500 EROSION COEFFICIENT = 1.0E-4 CRITICAL EROSION SHEAR STRESS OF THE MUD LAYERS = 0.05 CRITICAL SHEAR STRESS FOR DEPOSITION = 1000000.0 THRESHOLD FOR SEDIMENT FLUX CORRECTION ON TIDAL FLATS = 0.1 ADVECTION-DIFFUSION SCHEME WITH SETTLING VELOCITY = 0 /----------------------------------------------------------------------/ / SEDIMENT POINT SOURCES

/----------------------------------------------------------------------/ ABSCISSAE OF SOURCES = 83429.961 ORDINATES OF SOURCES = 361424.594 ELEVATIONS OF SOURCES = -6.0 WATER DISCHARGE OF SOURCES = 0.1 VALUE OF THE TRACERS AT THE SOURCES = 0.0;1441.53 /----------------------------------------------------------------------/ /

TIDAL FLATS

/----------------------------------------------------------------------/ TREATMENT ON TIDAL FLATS FOR VELOCITIES

TREATMENT ON TIDAL FLATS FOR K-EPSILON

TREATMENT ON TIDAL FLATS FOR TRACERS

TIDAL FLATS

=YES

OPTION FOR THE TREATMENT OF TIDAL FLATS

Final version

WL2020R13_131_6

A19

DEPARTMENT MOBILITY & PUBLIC WORKS Flanders hydraulics Research Berchemlei 115, 2140 Antwerp T +32 (0)3 224 60 35 F +32 (0)3 224 60 36 waterbouwkundiglabo@vlaanderen.be www.flandershydraulicsresearch.be