The sediment samples used for the erosion experiments were collected from two sites of the Elbe in the area of the Port of Hamburg. Upstream of the port area, the Elbe divides into two parts, the Norderelbe (NE) and the Suederelbe (SE), which subsequently reunite in the center of the port and surround the island of Wilhelmsburg (see Figure 1 ). The sediment samples were collected during two measurement campaigns, one each on the NE and SE, at known sedimentation hotspots of cohesive sediments. The SE-campaign was conducted in June 2023 and located in a ship turning point. As the ship turning point provides a widening of the cross-section, sediments accumulate due to relatively low flow velocities and evolving flow shadows. The sediment samples at this site were collected during flood slack water. The NE-campaign was carried out in November 2023 in the harbor entrance to Baakenhafen. Reduced flow velocities lead to fine-sediment accumulation in this area as well, whereas at this location the extension of the accumulation tends to reach further into the main flow area. The NE-samples used for the experiments in this paper were collected during ebb slack water.

The sampling device used was a sediment corer developed at the Institute of River and Coastal Engineering (Patzke et al., 2019, 2022). The corer penetrates the upper sediment layers by its own weight and extracts cores with a diameter of 20 cm and a maximum height of 1.20 m. The penetration depth of the corer varies depending on the attached extra weights and the present bed densities. In both campaigns, the corer penetrated about one meter into the bed. In Witt et al. (2023) it was shown that the device is capable of collecting naturally stratified bed samples, thus the properties (grain size distribution, loss on ignition) of the sediment used in this study can be seen as an average of the top one-meter layer of the bed.

The sediment samples from the two sampling areas exhibit a comparable grain size distribution (GSD) (see Table 1 ). The GSD was derived by a combined sieve and hydrometer analysis according to DIN EN ISO 17982-4. The clay content (d < 2 μm) of both samples with around 24% is way above the threshold of 5-15% found in literature for dominant cohesive behavior of the sediment (e.g. van Rijn (1993), clay defined as d < 4 μm). In terms of silt (2 – 63 μm) and sand (> 63 μm) content the samples vary slightly with a three percent difference in each case, with the sample from SE showing a higher silt and correspondingly a lower sand content. For both samples, the sand fraction is in the range of fine sand (< 200 μm) only. Consequently, the values for the median diameter D₅₀ and D_{50, mud}, which is defined as the median diameter of the mud-fraction (clay and silt), are comparable. The loss on ignition is higher for sampling site SE (11.5%) than for site NE (8.4%) indicating a higher organic content.

For all conducted experiments the SSC evolution in the measuring chamber was captured and the corresponding erosion rates were derived (see sections 2.3.2 and 2.3.3). In Figure 3 the results for both the SSC and the erosion rates for three different exemplary bulk densities (1100, 1150 and 1200 kg/m³) for the two sampling sites SE and NE are shown. The chosen densities represent the range of tested densities for which data from both sampling locations is available and illustrate the observed trends and relationships. The figure gives an impression of the generally high repeatability of the SSC measurements and erosion rate calculations, especially when considering the rapid change of erodibility with increasing bulk densities of the sediment mixture.

Focusing on the SSC plots against experimental time (panel (A)-(C) for SE and (D)-(F) for NE) the strong relationship between ρ_b of the samples and the observed SSC evolution becomes evident. For sediment from SE at ρ_b =1100 kg/m³ the first increase in SSC on the applied linear scale becomes visible at around minute 50 and the curve rises quickly over the following load steps to reach values of 20 g/l after 120 minutes duration. The SSC evolutions for ρ_b =1150 kg/m³ and ρ_b =1200 kg/m³ exhibit a consistent shift to the right, meaning a longer experimental duration and thus higher applied bed shear stresses are needed to cause a similar SSC-increase. For ρ_b =1200 kg/m³ the SSC reaches values of < 1 g/l after the maximum duration of 3.25 h. The decreasing erodibility with rising bulk densities of the slurry is caused by an increased amount of cohesive particles per volume element and therefore more and stronger bonds between the single particles. Since the samples tested are homogenized and have uniform erosion characteristics over the depth of the sample, the gradient of the SSC curves is theoretically constant for each load step and increases with the applied bed shear stress. This relationship is generally well represented in the measured data, especially for higher applied bed shear stresses.

For load-steps of lower bed shear stresses most of the SSC curves approach a plateau asymptotically. This observation is confirmed by the corresponding plots of ε, which are shown under the SSC-plots (panel (G)-(I) for SE and (J)-(L) for NE). For low τ_b the erosion rates tend to peak at the beginning of each load-step and decrease over the remaining load step duration. This trend in ε-evolution is commonly observed in experiments with density-stratified samples and increasing erosion resistance over the depth of the sample (type I erosion, see section 1). The reason why the described trend is also apparent in the presented experiments, despite working with homogeneous samples, is seen in the fact that during the first load steps only the topmost part of the sample (< 1 mm depth) is eroded. In this top-layer particles are incorporated into the grain structure to varying degrees and therefore exhibit a differing erosion resistance, meaning that the erosion resistance (as also the applied load) follows a probability distribution (Winterwerp et al., 2012). This leads to decreasing erosion rates over the load-step duration, since the number of particles that might get eroded declines. With increasing load and erosion depth, the erosion rates become more and more constant over each load-step and reach values of up to 5 x 10^-3 kg/(m²s).

The setup 2 of the utilized C-GEMS generally records an earlier increase in SSC and ε compared to setup 1, as expected due to the slightly higher τ_b applied. Comparing the results of the two sampling sites SE and NE, the evolution of SSC- and ε indicates a higher erodibility for material from site NE. For example, for ρ_b =1200 kg/m³ the measured final SSC for site SE is < 1 g/l, whereas values of ~4-5 g/l are measured for site NE. For the other densities, this relationship is reflected by faster increasing SSC and erosion rates as well.

From the SSC and ε-evolutions, the mean values of the parameters per load-step have been calculated (see section 2.3.4) and are shown in Figure 4 (SE) and Figure 5 (NE) against the applied bed shear stress. Focusing on the exemplary bulk densities for sampling site SE ( Figure 4 ), the diagrams illustrate that for both methods to derive τ_c , the erosion threshold increases together with the bulk density. Additionally, the gradient of the regression lines consistently decreases with increasing bulk density. These fundamental relationships are in line with the observations described in section 3.2 and theory and indicate a lower erodibility for higher bulk densities. For the illustrated densities the erosion rate method leads to values for τ_c of 0.053 – 0.430 N/m². While for bulk densities of 1100 kg/m³ and 1150 kg/m³, this method yields a small semi-range (SR) for calculated τ_c in repetition tests of < 0.005 N/m², the spread of the derived values rises drastically for ρ_b =1200 kg/m³ (SR=0.086 N/m²). Also, the calculated coefficients of determination (R²; calculated based on log-transformed data) indicate a strong linear relation between the applied bed shear stress and ε for bulk densities 1100 kg/m³ and 1150 kg/m³ (R² ≥ 0.94), but drops for ρ_b =1200 kg/m³. The SSC-method yields values for τ_c in the range of 0.060 – 0.159 N/m² and results in small semi-ranges of τ_c in repetition tests (SR < 0.008 N/m²) and high R² values (R² ≥ 0.91) for all three shown bulk densities. Compared to the ε-method the SSC-method leads to slightly higher τ_c values for ρ_b =1100 kg/m³ and ρ_b =1150 kg/m³, but significantly lower values for ρ_b =1200 kg/m³. The reason for the inconsistent results of the ε-method for ρ_b =1200 kg/m³ is seen in the fact that this bulk density is close to the application limit of the utilized C-GEMS devices for sediment from this specific site. Only a small amount of sediment is eroded over the experimental duration, leading to low erosion rates and large scatter in the erosion rate data (see Figures 3 , 4 ). The large scatter in erosion rate data in turn yields low R² values and a large spread in derived erosion thresholds. Since, in contrast to ε, the SSC is accumulating over time, the SSC data exhibits much less scatter and allows a more reliable determination of τ_c even at the application boundary. In the upper part of Table 3 , the results for all experiments conducted with material from sampling site SE are summarized. The data derived in the additional experiments underlines that both applied approaches yield comparable values for τ_c , SR and R², except at the upper application limit where the SSC-method is superior to the ε-method.

The described trends and relationships are also valid for the experiments conducted with sediment from site NE. The results for the three exemplary bulk densities are illustrated in Figure 5 and the data for all experiments is summarized in the lower part of Table 3 . As expected from the data shown in Figure 3 , the erosion thresholds derived for site NE are generally lower than for site SE. The higher erodibility of the NE-sediment also leads to a shift of the application limit of the GEMS device to higher bulk densities. While for SE sediment for ρ_b =1200 kg/m³, only the SSC-method yielded reasonable results, for sediment from site NE the ε-method was still applicable for this bulk density as well. The application limit was reached at ρ_b =1250 kg/m³ for site NE, where the ε-method again led to unreasonably high values for τ_c , a high spread of the derived values and low R².

For NE-sediment with ρ_b =1100 kg/m³, ρ_b =1150 kg/m³ and ρ_b =1250 kg/m³ additional experiments have been conducted to check and quantify the repeatability of the experimental setup and approaches to derive τ_c . The additional repetitions were conducted around one month after the initial experiments, which were carried out in the first two weeks after the measurement campaign. In the additional experiments, slightly lower values for τ_c were derived. While the SSC-method for ρ_b =1150 kg/m³ in the initial experiments yielded a mean τ_c of 0.081 N/m² (SR=0.008 N/m², n=2), for the additional repetitions a mean of 0.075 N/m² (SD=0.005 N/m², n=5) was derived, leading to an overall value of 0.077 N/m² (SD=0.007 N/m², n=7) (see Table 3 ). The measured change might be caused by a change in sediment properties evoked by biological activity. However, due to the relatively small database and small measured changes in τ_c , this is only an assumption yet and needs further investigation.

The two methods used to derive τ_c delivered comparable results for most of the bulk densities tested, with a small spread in derived values in repetition tests and strong linear correlation (in log-space) between applied bed shear stress and measured SSC/ε (slightly higher for SSC-method). Especially because of the better performance of the SSC-method at the upper limits of application, this method should be preferred though and the data derived by this method is used to fit models for τ_c in section 3.3. To investigate the influence of the applied values of the thresholds for SSC (SSC_T) and ε (ε_T ) on the derived values for τ_c , a sensitivity analysis was conducted using the example of site SE ( Figures 6A, B ). The chosen thresholds were varied by ±50%. As expected from the decreasing slope of the regression lines with increasing bulk density ( Figure 4 ), the total variation of the derived τ_c also increases with ρ_b . Especially for the lower and middle values of the tested density range, the variation is relatively small, e.g. for SSC-method and ρ_b =1150 kg/m³ the variation of SSC_T of ±50% leads to a change in the derived τ_c of +10.6% and -15,7%, respectively. Figure 6C additionally shows an analysis of the influence of different time periods (t_av ) over which the average value of ε is calculated for each load step when using the ε-method. A value of t_av =5 min means only the first five minutes of each load step were used for averaging. This analysis can be interpreted as the expected sensitivity of the calculated τ_c with respect to a corresponding change in load step duration. Since the magnitude of ε decreases significantly over the applied load step duration for load steps of lower τ_b (see Figure 3 ), a shortening of t_av leads to higher calculated averages of ε and consequently to a decrease of the derived values of τ_c . For the SSC-method, a similar analysis is not feasible based on the dataset, since this parameter accumulates over the duration of the experiment und thus the actually applied load step duration directly influences the SSC in subsequent load steps.

Summaries of τ_c derived from field and laboratory experiments in other studies can be found in Debnath and Chaudhuri (2010) and Houwing (1999). As stated in section 1, the derived values are strongly influenced by the erosion device used, the specific sediment properties and the applied definition of the erosion threshold. There is only a limited amount of data available for bulk densities in the range that was tested in this study. However, the studies from Righetti and Lucarelli (2007), who tested benthic lake sediments with bulk densities from 1025 – 1175 kg/m³ in a Sedflume and derived values for τ_c in the range of 0.01 – 0.23 N/m², and Amos et al. (1997), who worked with mud flat sediments with ρ_b from 787 – 1273 kg/m³ in a Sea Carousel and derived values for τ_c of 0.15 – 0.5 N/m², show that the findings in the present study are well covered by the observations in previous studies.

The data on τ_c derived by SSC-method in section 3.2 is used to fit two mathematical models. The first model is a power-law relation between ρ_db and τ_c (Equation 1) and the second model is a more complex and physically based model shown in Equation 3, proposed in Chen et al. (2021). The application of both models requires the calculation of ρ_db as shown in Equation 2. From the assumed density of the particles of the mud-fraction ρ_pm =2575 kg/m² (cf. Malcherek (2010)) and sand particles (ρ_pps =2650 kg/m³) the weighted particle density for the sediment mixture is calculated as ρ_s = ρ_ps + ρ_pmP_m . The derived values for ρ_db are shown in Table 3 . In Figure 7 the fitted models are illustrated. Both models have been fitted to the data obtained from the two sampling sites separately and additionally to the whole dataset using least-squares method. For site NE only the derived data from the initial two experiments (one in case of ρ_b =1250 kg/m³) was used, even if more repetitions have been carried out, due to the unclear effect of the time span between the experiments (see section 3.2). The obtained fit parameters and the resulting root mean squared errors (RMSE) are summarized in Table 4 .

The power law-fit for site SE suggests an almost linear increase of τ_c over the considered range of bulk densities (exponent n=1.15) while for site NE the power-law relation is more distinct. However, the derived model for NE underpredicts the erosion thresholds in the lower range of ρ_b =1100 – 1150 kg/m². Fitting the model on the whole dataset (SE+NE) yields a reasonable model for the whole considered range of densities and both sampling sites. The derived fit-parameters for this model are m=1.83 x 10^-6 and n=1.96. Earlier investigations, which were carried out by Owen (1975) and Thorn and Parsons (1980) (see see Zhu et al. (2008)) with remolded sediment samples from four different estuaries, led to values for m of 6.85 x 10^-6 resp. 5.42 x 10^-6 and values for n of 2.44 resp. 2.28. An extensive overview of derived fit-parameters from other studies carried out with different muds and erosions devices can also be found in Schweim (2005), with m ranging between 2.5 x 10^-8 – 1.5 x 10^-2 and n between 0.73 – 2.44. Although the fit-parameters derived in this study are comparable to earlier studies, their values vary significantly with the underlying dataset both for this study and earlier works.

The model proposed in Chen et al. (2021) includes only a single fit-parameter, denoted as A (units J/m²). Additional sediment properties used to fit the model are shown in Table 1 . The mud-fraction is calculated as the sum of clay and silt fractions and for parameters ρ_ps and ρ_pm the same assumptions are made as described for the power-law model. The derived values for the fit parameter A are 8.56 x 10^-6 J/m² based on SE-data and 8.38 x 10^-6 J/m² based on NE-data. The fit on the whole dataset yielded A=8.43 x 10^-6 J/m². It stands out that the fit-parameters show only a marginal spread. In Chen et al. (2021) the authors applied the model to multiple natural muds from different locations and concluded that the value of A is generally in the order of 10^-6 or 10^-5 J/m² and specifically for natural mud and sand-mud mixtures in the range of 2.86 x 10^-6 – 1.04 x 10^-5 J/m² (if no field data is available the authors propose a general value of A=3.97 x 10^-6 J/m²). The values for A derived in this study are covered well by the findings of the authors, lying close to the mean of A-values determined by them. The small spread in A suggests that the consideration of the physical parameters describing the mud-fraction of the sediment leads to a more profound model than the power-law model also applied, relating the evolution of τ_c to the actual sediment properties and leaving a lower degree of freedom for the empirical fit. Figure 7B additionally shows that the slightly lower erosion thresholds for NE-sediment might be at least partly caused by the (also slightly) different grain size distribution, since the same fit-parameter (solid lines) leads to lower τ_c values. The lower spread of A is also reflected in the calculated RMSE. While for the power-law fit the RMSE rise considerably for both individual sites when the model is fitted on the whole dataset and not on the site-associated data only, for this model, the RMSE is almost independent of the chosen database for the model fit (see Table 4 ). A comparison of the RMSE of the two different models fitted on the whole dataset additionally reveals a slightly better adaption of the Chen-model to the measured erosion threshold data overall. A sensitivity analysis conducted for the value of the fit-parameter A with respect to the applied values of SSC_T and ε_T , using the example of site SE (fitting lines shown in Figures 6A, B ), shows that a variation of the values of ±50% leads to a range of A of 7.30 x 10^-6 – 1.0 x 10^-5 J/m² for ε-method and 6.54 x 10^-6 – 1.01 x 10^-5 J/m² for SSC-method. The analysis of different time periods t_av used for averaging ε for each load step when using the ε-method ( Figure 6C ) shows that the applied reduction of t_av from 15 to 5 min leads to a range of A of 6.55 x 10^-6 – 8.91 x 10^-6 J/m². Despite the applied variation of SSC_T, ε_T and t_av the derived values for A are well covered by the expected range of A, accordingly.

Some additional considerations regarding the results obtained and the experimental procedure are discussed in this section. Since the C-GEMS is a closed system and the eroded sediment accumulates throughout the duration of the experiment, the question arises, as with all closed erosion devices, whether deposition occurs simultaneously with erosion. In general, two concepts to describe the onset of deposition of cohesive sediment exist: i) the existence of a critical shear stress for deposition τ_d , which is lower than τ_c , meaning deposition and erosion are mutually exclusive (Krone, 1962; see Whitehouse et al. (2000)) and ii) the assumption that deposition occurs continuously also at higher shear stresses (see e.g. Winterwerp et al. (2021)). While modeling of in-situ data may yield better results under the assumption of continuous deposition, many laboratory studies have demonstrated the existence of a critical shear stress, below which deposition does not occur (see Sanford and Halka (1993)). This discrepancy may be due to differences in scale and turbulence conditions and the absence of still water/flocculation periods in laboratory experiments. In the C-GEMS, the turbulent mixing in the erosion chamber is supported by the constant suction and feeding of the suspension and flocculation of eroded particles is hardly possible because the particles are exposed to high shear rates in the different pumping cycles, resulting in very small floc sizes and (theoretical) settling rates. Both of these factors would significantly harm deposition. Even if these flocs/particles briefly were to come into brief contact with the sediment surface it is not unlikely that they would be able to withstand the applied bed shear stress, so they would be immediately eroded again. Nevertheless, we cannot completely rule out the presence of deposition during the experiments, so the derived erosion rates can be interpreted as net erosion flux.

In section 3.3 it was shown that the Chen-model can explain the observed higher erodibility of sediments from site NE at least partly with the sediment composition (cf. Figure 7B ), slightly lower predicted τ_c for NE-sediment for the same value of A). To make a statement whether this is the only reason for the noticed difference in erodibility or if other reasons e.g., a difference in the kind of organic matter, influence the erosion behavior, further investigation is needed. Focusing on the fit-parameter A, in Chen et al. (2021) it is stated that the value of A is supposed to be site-specific because it describes the bed surface roughness and the cohesion strength, which is, among others, influenced by the mineral composition of the sediment and the pore water environment. Since the sampling sites investigated in this study are only a few kilometers apart and comparable in characteristics, it seems plausible to assume a broad agreement in characteristics. This assumption in turn leads to the expectation, that the derived parameters for A are in the same range for both sites, as demonstrated in this study.

Another noteworthy point is that the lowest bulk densities tested for site SE (ρ_b =1050 kg/m³ and 1075 kg/m³) are below the structural density (as stated in section 2.3.2). Despite the adjusted, shortened time between sample preparation and the beginning of the erosion experiment for these densities, the onset of the settling process presumably led to higher densities below the lutocline than initially set. The increased densities might have led to an overestimation of the erosion threshold from the conducted measurements and explain the underprediction of τ_c for these bulk densities by the models.

Looking at the applied techniques to derive τ_c , erosion rate and SSC method, the chosen values for the thresholds of ε and SSC are obviously debatable. However, no other method to derive τ_c is available that does not include either defined threshold values, the differentiation of erosion-modi or the definition of erosion/no-erosion states. Indeed, the implementation of a kind of threshold (that is always arguable) is imperative, since with highly sensitive sensors even at very low applied bed shear stresses measurable erosion rates can be detected (as shown in this study), which rather result from particle/background erosion than from surface erosion. This observation rises the question if a single value for τ_c is a reasonable concept to describe the initiation of motion of cohesive sediment. However, since most the erosion models in use depend on this type of data, the practical benefits are beyond doubt. Some models (e.g. Parchure and Mehta (1985)) also incorporate the unsharp initiation of motion. The benefit of the tested methods in this work in deriving these data is, that only one single SSC-/ε-threshold is used and that the thresholds are transferable to other devices straightforward, which is hoped to contribute to better comparability of derived erosion thresholds of future investigations.

In this work, an improved version of a closed microcosm system (C-GEMS) proposed in Patzke et al. (2022) is presented. The device is used to perform erosion experiments with freshly deposited cohesive sediments collected from two sites on the River Elbe in the area of the Port of Hamburg. For the erosion experiments, the sediments were homogenized and diluted to bulk densities from 1050 – 1250 kg/m³, which is the observed range of bulk densities in earlier investigations over the top layer of the sediment bed (multiple ten-centimeters to about one meter). During the experiments a maximum of 13 load steps were applied with increasing bed shear stresses from 0.028 to 0.33 N/m² (setup 1) resp. 0.034 – 0.40 N/m² (setup 2). The evolution of SSC in the measuring chamber of the C-GEMS was recorded using a Hach^® Solitax sc probe and the corresponding erosion rates were calculated from the change in SSC over time. The results of the SSC- and ε-evolutions indicate a distinct relationship between the bulk density of the cohesive sediment-mixture and the erodibility and exhibit good repeatability. The further investigations focused on the derivation and quantification of the surface erosion thresholds. Two methods were applied for this purpose: i) interpolation of the log-log-regression of mean erosion rates to the threshold value of ε=1 x 10^-5 kg/(m²s) and ii) interpolation of the log-log-regression of mean SSC to the threshold value of 40 mg/l. The derived values for τ_c for both methods showed a good agreement for the considered range of bulk densities except at the upper application limits of the utilized C-GEMS device, where only the SSC-method yielded reasonable results due to significantly lower scatter in the SSC-data compared to ε-data. The reason for the lower scatter of the SSC-data is seen in the accumulating nature of this parameter, as the C-GEMS is a closed system. Comparable studies with other erosion devices support the conclusion that deriving erosion thresholds based on an accumulating parameter leads to more robust results (Amos et al., 2003; Ha and Ha, 2021). Because of the better performance at the application limits it’s suggested that the SSC-method should be preferred over the ε-method when working at the application limits of the device. The derived mean τ_c for location SE (ρ_b =1050 – 1200 kg/m²) with SSC-method are in the range of 0.037 – 0.151 N/m² and for location NE (ρ_b =1100 – 1250 kg/m²) in the range of 0.063 – 0.305 N/m². For NE-site with ρ_b =1150 kg/m³, a total of seven repetitions has been carried out yielding a mean τ_c of 0.077 N/m² with a standard deviation of 0.007 N/m², underlining the reproducibility of the utilized experimental procedure. In a last step the derived τ_c -data was used to fit two mathematical models, a power-law relationship between dry bulk density and τ_c and a model proposed in Chen et al. (2021). The fitting of the power-law model showed a strong dependency of the derived fit-parameters on the respective data used for the fit (SE, NE or SE+NE), while the fitting of the Chen-model yielded practically the same value for the fit-parameter A independent of the chosen database. Using the whole dataset to fit the models, the power-law fit yielded values of the prefactor m=1.83 x 10^-6 and the exponent n=1.96. The fitting of the Chen-model on the whole database led to A=8.43 x 10^-6. All derived fit-parameter are well-covered by the range of values reported in previous studies. Since the Chen-model relates the change in τ_c to the physical properties of the sediment mixture, shows a better adaption to the measured data and the value of the fit-parameter is stable independently of the underlying data subset, it’s concluded that the model is more profound and should be preferred over the highly empirical power-law model.

In further investigations, different models for the erosion rate will be applied to the collected data. Taking the derived relationship between ρ_b and ε for the examined samples into account, the often used linear model based on the work of Partheniades (1965) (see McAnally and Mehta (2000)) reading ε = m(τ_b -τ_c ) for τ_b > τ_c , with m as erosion constant, presumably allows only an unsatisfactory adaption to the collected data. In combination with in-situ measured density profiles the derived relations of τ_c and ε enable the generation of depth profiles of the regarding values for the upper layers of the sediment bed. Additionally, it is planned to utilize the presented C-GEMS device for erosion experiments with (almost) undisturbed naturally stratified sediment samples as well. These experiments will be conducted i) with sediment cores extracted from the navigation channel of the Elbe right after withdrawal directly on the vessel and ii) by applying the C-GEMS in-situ in the tidal flats of the Elbe by mounting the system on a piercing cylinder.