Rheological analysis and rheological modelling of mud sediments: What is the best protocol for maintenance of ports and waterways?

Abstract Natural mud sediments display complex rheological behaviour like thixotropy, viscoelasticity and yield stress. These rheological characteristics can significantly vary over depth, from one mud layer to another, as each layer can have a different density and composition. Fast and reliable measurements of yield stresses of mud samples are important for maintenance operations in ports and waterways. These protocols, performed in the laboratory, should give a rheological fingerprint which is representative of the in-situ behaviour of the mud. In this article, we show that our recently developed stress ramp-up rheological protocol is a time-efficient and well-grounded protocol to determine the yield stresses of natural mud samples by comparing with other existing well-grounded protocols. In this study, we also refine the stress ramp-up protocol such as to reduce the experimental time for different mud layers based on their densities. The protocol was tested on a large number of mud samples obtained from different locations/depths of the Port of Hamburg, Germany. An empirical model is proposed to fit the two-step yielding behaviour that the mud samples exhibit. The model captures the two-step yielding phenomenon in mud samples quite well, within the density range of 1050–1200 kg. m−3. This two-step yielding is a feature of mud samples as found in various harbours and estuaries worldwide in rheometry.


Introduction
Mud sediments consist of water, clay minerals, sand, silt, and organic matter. Usually, (fluid) mud exhibits a complex rheological response including viscoelasticity, yield stress and thixotropy. It is already known that the rheological and cohesive properties of mud sediments are dependent on their density and the amount of organic matter (Malarkey et al., 2015;Parsons et al., 2016;Paterson et al., 1990;Paterson and Hagerthey, 2001;Schindler et al., 2015;Shakeel et al., 2019;Tolhurst et al., 2002;Wurpts and Torn, 2005). Knowledge of these rheological properties of mud sediments is very important to predict density currents and fluid mud flows, which in turn impact turbulence, navigation in channels, maintenance of ports (dredging activities) and general quality control of water in coastal areas. Being able to quantify the properties of fluid mud allows to set the appropriate boundary conditions for computational modelling of sediment transport, which are required to facilitate dredging operations (Whitehouse et al., 2000) and the maintenance of navigational channels (Kirichek et al., 2018;May 1973;Parker and Kirby, 1982).
The presence of fluid mud in a port is linked to the criterion for the nautical bottom, which is defined as the level beyond which the manoeuvrability and controllability of the ship becomes difficult due to the contact of mud layer with the ship's keel. For practical reasons, the nautical bottom was until recently defined as a critical fluid density (McAnally et al., 2007). However, there is another criterion (i.e., yield point) which is quite important for the ports where the organic matter content is significantly varying as a function of location within the port. Therefore, for these ports, only density may not be enough to predict the rheological characteristics of mud and the nautical bottom can be defined on the basis of yield stresses of mud (Wurpts and Torn, 2005). Thanks to the technical progress in in-situ monitoring, it is now possible to assess the in-situ yield stresses (Kirichek et al., 2020), even though the techniques still require further improvements due to the complexity of the fluid mud systems. In order to calibrate and improve the in-situ equipment, extensive laboratory tests are still needed. Claeys et al. (2015) reported a protocol to analyse the rheological properties of mud sediments in the laboratory using a vane-type rheometer. Their main objectives were to attain the equilibrium flow curves (EFC) of the samples with a good repeatability and to deduce the dynamic yield stress of the disturbed samples. Their protocol starts with a stress growth test at a very small shear rate (1 s − 1 ) to obtain the undrained shear strength of "undisturbed" samples, followed by a pre-shearing step at a very high shear rate (1000 s − 1 ), to completely disturb the sample. The stress growth test is highly dependent on the applied shear rate and, therefore, the selection of suitable shear rate for the samples of different consistencies (i.e., fluid mud, consolidated, etc.) is very critical (Rogers et al., 2010;Stokes and Telford, 2004;Yuan et al., 2017). This protocol is, therefore, not very straightforward to analyse the yield stress of "undisturbed" mud samples. Moreover, this protocol is based on nine cycles or more (depending upon the selected shear rates) of applying/removing shear rate and the total time of an experiment is about 15-20 min/sample. Therefore, this protocol is not suitable to measure yield stresses of a large number of samples. Recently, we investigated another rheological method (stress ramp-up test) to assess the rheological properties (in particular the yield stresses) of mud samples (Shakeel et al., 2020b).
Different empirical rheological models have been used in literature for the cohesive mud sediments (Coussot, 1997;Van Kessel and Blom, 1998;Bai et al., 2002;Babatope et al., 2008;Huang and Aode, 2009). Van Kessel and Blom (1998) used a thixotropic model developed by Toorman (1997), which is an extension of the Moore model (Moore, 1959) for characterizing the rheological behaviour of estuarine mud. Fonseca et al. (2019) analysed the rheological properties of mud sediments from Port of Santos, the Port of Rio Grande, the Port of Itajaí and the Amazon south navigation channel in Brazil, by fitting the data with a Bingham model. Yang et al. (2014) investigated the rheological behaviour of mud sediments from shoal of the Hangzhou Bay, Yangtze River and Yangcheng Lake, China with the help of the Herschel-Bulkley model. Huang and Aode (2009) studied the use of a Dual-Bingham model and Worrall-Tuliani model (Worrall and Tuliani, 1964) to fit the rheological data of mud samples obtained from two different locations of Hangzhou Bay, China. Similarly, Xu and Huhe (2016) used a Dual Herschel-Bulkley model to examine the rheological data of estuarine mud obtained from Lianyungang, China. All these reported empirical models are usually effective for single-step yielding in the sample. However, no empirical/theoretical model is yet available in literature to capture the behaviour of two-step yielding. This two-step yielding feature has already been reported for mud samples from other locations (Huang and Aode, 2009;Nie et al., 2020;Yang and Yu, 2018). These two yielding points are typically associated to the breakage of interconnected network of flocs (first yield point) followed by the collapse of individual flocs into smaller flocs or particles (second yield point) (Nie et al., 2020). Two-step yielding is also typical for wall slip in rheometry, specifically at low shear stress/shear rate (Barnes, 1995).
In the present study, a comparison between stress ramp-up (Shakeel et al., 2020b), Claeys protocol and conventional rheological protocols (i. e., equilibrium flow curve, shear rate ramp up and ramp down, etc.) is carried out in order to investigate the differences between protocols. Stress ramp-up protocol is further refined such as to reduce the experimental time for the different mud layers based on their densities. Moreover, an empirical model is also proposed to capture the dual yielding behaviour in mud samples.

Experimental
In this study, 1 m core sampler was used to collect the mud samples from different locations of Port of Hamburg (PoH), Germany (Fig. 1, Table S1). The obtained samples were then subsampled into different layers on the basis of their consolidation stage (which is a function of depth). The dry density of the mud samples was considered to be 2650 kg. m − 3 (Coussot, 1997), which is actually the density of solid particles/quartz. The bulk density of the mud samples was estimated by the oven drying method (Coussot, 1997). The mud samples obtained from different locations of port were analysed with the aim to derive an empirical rheological model for the observed two-step yielding behaviour. In order to compare different rheological protocols, only mud samples from location 'L6' (see Fig. 1, Table 1) were considered. Particle size distribution (PSD) for different mud layers, collected from location 'L6', was measured by using static light scattering technique (Malvern MasterSizer, 2000MU). The mud samples were diluted enough, without any pre-treatment, in order to achieve the required range of obscuration (i.e., less than 20%) in the Malvern equipment (Shakeel et al., 2019(Shakeel et al., , 2020d. The obtained results are presented in supplementary material, (Fig. S1). The properties of the mud samples from location 'L6' are presented in Table 1. All the mud samples were homogenized by mild hand stirring, prior to the rheological experiments.
HAAKE MARS I rheometer (Thermo Scientific, Germany) was used to perform rheological experiments with Couette (CC25 DIN Ti) and vane (FL 22) geometries in rheometer cup (diameter = 27 mm). Peltier controller system was used to maintain the temperature of 20 • C for each experiment. The repeatability of the experiments was checked by performing each experiment in duplicate and the repeatability error was always less than 2%. Grooved bob geometry (grooves of 0.5 mm depth, 1 mm width and spaced by 2 mm) was also used to analyse the effect of wall slip (Barnes, 1995). The similar results obtained from smooth and grooved bob geometries (data not shown) indicated the absence of wall slip. Following types of rheological tests (Sec. 2.1-2.5), already reported in literature, were performed to analyse the mud samples:

Stress ramp-up test
The controlled stress mode of the rheometer was used to perform stress ramp-up tests. A linear increasing stress is applied at a rate of 0.1-10 Pa/s, till the shear rate reaches 300 s − 1 (Fig. 2a). The resultant  rotation of geometry is recorded, which eventually used to calculate the shear rate and apparent viscosity. (Claeys et al., 2015) First, a stress growth test is performed at a shear rate of 1 s − 1 for 30 s followed by a pre-shearing step at a particular shear rate (300 s − 1 for Couette and 150 s − 1 in case of vane) for 15 s. The desired shear rate is then set for a period of 100 s to get the steady state value of stress, followed by a low shear rate step (0.001 s − 1 ) for 6 s. The same cycle is repeated for each selected shear rate from 100 to 0.01 s − 1 by performing 11 cycles (Fig. 2b).

Equilibrium flow curve (EFC) test
The increasing or decreasing equilibrium flow curves are performed by linearly increasing or decreasing the shear rates (within the range of 0.01 and 100 s − 1 ) after the time period of 100 s to get the steady state stress values, without performing any pre-shearing, stress growth and low shear rate steps ( Fig. 2c and d). The duration of 100 s was enough to attain the steady state values of stresses for the applied shear rates (see Fig. S2).

Shear rate ramp up and ramp down test (CSRT)
Controlled shear rate ramp up and ramp down (CSRT) tests consist of shear rate ramping up and down by linearly increasing the shear rate from 0.01 to 100 s − 1 and then linearly decreasing from 100 to 0.01 s − 1 . These experiments are performed without giving enough time to attain steady state values. However, at lower shear rates (from 0.01 to 5 s − 1 ) a slower increasing or decreasing rate is used as compared to the higher shear rate range (from 5 to 100 s − 1 ), as can be seen in Fig. 2e.

Pre-shear test
In this protocol, similar to Claeys protocol, a pre-shearing step is performed at a shear rate of 100 s − 1 for 100 s, followed by a linear shear rate ramp down step from 100 to 0.01 s − 1 . The change from one shear rate to another during the ramp down is too short for the system to reach steady-state stress values. The pictorial representation of the protocol is presented in Fig. 2f.

Existing rheological models
Several rheological models have been reported in the literature to fit the flow curves of the yield stress materials such as Bingham model (Bingham, 1922), Dual-Bingham plastic model (Huang and Aode, 2009), Papanastasiou model (Papanastasiou, 1987), Herschel-Bulkley model (Herschel and Bulkley, 1926), Worrall-Tuliani model (Worrall and Tuliani, 1964), Toorman model (Toorman, 1997), etc. However, Bingham, Herschel-Bulkley and Worrall-Tuliani models have been found in the literature to well describe the flow curves of different types of sediment. These models are given by: Where K is the consistency index (K = μ ∞ when n = 1), n is the flow behaviour index, μ ∞ is the viscosity at higher shear rate, τ 0 is the yield stress and τ B is the Bingham yield stress, as at high shear, Eq. (3) reduces to τ(γ →∞) = τ B + μ ∞γ .

Proposed rheological model for two-step yielding
From experimental evidence (i.e., rheometry), it appeared that the mud samples collected from the Port of Hamburg exhibit a two-step yielding behaviour (Shakeel et al., 2020a(Shakeel et al., , 2020b. This two-step yielding is associated with two characteristic shear stresses that we (Shakeel et al., 2020b) have termed static yield stress τ s and fluidic yield stress τ f . These names were chosen as τ s represents the transition between a (solid) cohesive structure and a structure consisting of mobile clusters of particles and τ f depicts the transition between this structure and a structure consisting of even smaller mobile clusters of particles. The shear stress τ(γ) is, at lower shear rate (γ), depending on τ s and can be fitted by the function τ stat (τ s ,γ) and, at higher shear rate (γ), depending on τ f and fitted by the function τ fluid (τ f ,γ 0 ,γ), where the transition between the two functions occurs at a shear rate defined as γ 0 .
We propose the following equation that enables to fit τ(γ) over the whole range of shear rates by introducing a step function α(γ 0 ): where the step function α is given by: This function gives α(γ <γ 0 ) = 1 and α(γ >γ 0 ) = 0 with a transition at γ =γ 0 , whose sharpness depends on the value of k. Different values of k (i.e., 1, 10 and 100) were tested for the experimental data fitting of mud sediments (data not shown). The best fitting was achieved by using k = 10 s, which was then used for the complete data fitting. The function τ stat is given by: This function can be seen as an adaptation of the Worrall-Tuliani model, whereby τ s /γ s can be identified as Δμ and 1/γ s as β, see Eq.
(3). The Worrall-Tuliani model includes an initial shear stress τ 0 which does not appear in Eq. (6). The reason is linked to the accuracy of the measurements: with modern devices, it is possible to estimate the shear stresses at very low shear rates and, hence, the value of τ 0 would be lost in the "noise" of the data at low shear. We found it, therefore, preferable to eliminate an unnecessary parameter by setting τ 0 = 0, which amounts to state that τ(γ = 0) = τ stat (γ = 0) = 0. From the analysis of this function, we find that: The shear rate γ s represents, therefore, the shear rate for which the stress τ s is half of its value. The curvature of the first function can be changed by varying γ s . The Worrall-Tuliani model assumes that τ(γ →∞) = τ ∞ + μ ∞γ and, hence, reduces to a Bingham model. This part has been left out of our model, since for higher shear rate, the stress does not follow a Bingham behaviour but display a second yield behaviour.
This second step in yielding behaviour is captured by the function τ fluid given by: This model is an adaptation of the Worrall-Tuliani model, whereby the coefficient d enables to tune the "sharpness" of the curvature. For mud sediments, d = 1 was used, and hence the function can be seen as a Worrall-Tuliani model with as origin γ =γ 0 , which can be written: From Eq. (3), we identified τ 0 as τ s , Δμ as τ f /(γ f −γ 0 ) and β as 1/ (γ f − γ 0 ). This function τ fluid does include a Bingham part, since where τ s + τ f can be seen as a (pseudo) Bingham yield stress for γ =γ 0 .
Furthermore, we find that: where τ s is the "true" yield stress at γ =γ 0 . We also have: where τ s + τ f 2 can be seen as a (pseudo) Bingham yield stress for γ =γ f . In short, the two-step yielding model contains six fitting parameters: γ 0 , τ s , γ s , τ f , γ f and μ ∞ (see Fig. 3). Two parameters, τ s (static yield stress) and τ f (fluidic yield stress) are quite important for practical applications.
The correspondence between these terminologies of yield stresses (static and fluidic) and the terminologies used in literature is presented in Table S2.

Comparison of different protocols for flow curve
Different protocols are available to measure the rheology (particularly the flow curve) of mud sediments, as detailed in the experimental section. In order to compare these different rheological protocols, mud sample 'L6-3' was selected. The two representations of the flow curves (shear stress vs shear rate and viscosity vs shear stress) for mud sample are shown in Fig. 4a and b, using different protocols with Couette geometry. The static and fluidic yield stress values were estimated from the viscosity declines in Fig. 4b. The mathematical method of yield stress determination is explained in more detail in Shakeel et al. (2019). Due to the curvature of the shear stress behaviour, yield stress ranges are given in Table 2, where necessary. It is clear from Fig. 4 that large differences are observed in shear stress/viscosity behaviour at lower shear rates/shear stresses. However, the differences at shear stress/shear rate values above the static yield point are also quite significant. This suggests that the mud sample 'L6-3', in the partially disturbed state, still has some structure depending on its shear history (i.e., protocol).
It is evident from Fig. 4b and Table 2 that "CSRT-ramp up", "EFCincreasing" and "stress ramp-up" displayed higher static and fluidic yield stress values as compared to the rest of the methods. All these methods started with a low shearing action and this implies a less disturbed state of the sample resulting in higher shear stressesbut also higher yield stress values. For the characterization of mud for harbour applications, it is important that the measurements are done in conditions reproducing in-situ situations. This implies that a method should be favoured in which the sample undergoes a transition from an (almost) undisturbed state to an unstructured state. All three methods cited above ("CSRT-ramp up", "EFC-increasing" and "stress ramp-up") correspond to this requirement, and as can be seen in Table 2, the values found for the fluidic yield stress are very close and in the range of 38-40 Pa. The fluidic yield stress is of importance for harbour applications, as it corresponds to the structural state most encountered in harbours. The authorities of Port of Emden for example use 100 Pa of yield point as a criterion for their nautical bottom (Wurpts and Torn, 2005), which corresponds to the definition of fluidic yield stress in our case. These three methods are, therefore, suitable to estimate the fluidic yield stress of mud for harbour applications. However, the EFC-increasing protocol is a longer test and from CSRT-ramp up the determination of the static yield stress is not very straightforward. Therefore, the stress ramp-up protocol is chosen for further analysis and empirical data fitting.
The fluidic yield stress values found for the other methods are also very similar, i.e., in the range of 26-29 Pa. These values are lower than the ones found from the other protocols, as the samples are in a less structured state. However, the static yield stress values found for the 'L6-3' mud sample are very different and lies in the range of 3.7-11 Pa. The lowest value is obtained by using Claeys protocol (Claeys et al., 2015), which is again due to the destruction of structure during measurement. This protocol is quite lengthy which implies that for liquid-like samples, settling can also occur in the samples during the time interval required for the measurement (∼ 20 min).
The comparison of the different protocols was also performed for mud sample using vane geometry, as the Claeys protocol was proposed with the vane type geometry. The results of different protocols for mud sample (L6-3), using vane geometry, are displayed in Fig. 5a. The differences between Couette and vane geometry results are discussed more into detail in Shakeel et al. (2020c). For Claeys protocol, the    conventional yield stress values are also reported in Table 3, in the form of an "undrained shear strength" obtained from the first step in the protocol (see Fig. 5b) and Bingham yield stress. The "Bingham yield stress" is determined from the extrapolation of shear stress vs shear rate curve for shear rate approaching to zero. In stress ramp-up test, an undisturbed mud sample is sheared and at the first viscosity decline the structure is partially broken or disturbed and the sample is more or less completely disturbed above second yield point. However, in Claeys protocol, the first step of applying shear rate of 1 s − 1 gives the peak stress of about 46.7 Pa (Table 3), which represents the fluidization or structural breakup for undisturbed mud sample. Therefore, this value corresponds to the fluidic yield stress (i.e., 40 Pa) in Couette, which also shows the structural breakup for undisturbed mud sample (blue dashed arrow in Fig. 6a). On the other hand, Bingham yield stress (i.e., estimated by extrapolating the flow curve to zero shear rate, as shown by the orange dashed line in Fig. 6b) for Claeys protocol represents the stress required to break the residual structure of mud sample (i.e., 26.5 Pa) or the minimum stress required to keep the material in flow condition. This value, therefore, corresponds to the difference between fluidic yield stress and the point where the static part is completely disturbed (i.e., 40-15 = 25 Pa), as shown by the orange dashed arrow in Fig. 6a. Furthermore, the undrained shear strength value is highly dependent on the applied shear rate, which is another limitation of the Claeys protocol. The differences in undrained shear strength values between Couette and vane geometry is also striking. This illustrates the importance of the measuring geometry to get information on the systems (Shakeel et al., 2020c).
In order to further verify this correlation of static and fluidic yield stresses with the conventional yield stresses (i.e., undrained shear strength and Bingham yield stress), stress ramp-up test followed by a constant high shear test (at 300 s − 1 for 500 s) and then stress ramp-down test was performed for mud samples having two different densities (1151 and 1256 kg. m − 3 ) obtained from location 'L8'. Claeys protocol was also carried out for the same samples to get conventional yield stresses. The results are presented in Fig. S3 and Fig. S4 with the corresponding yield stress values in Table S3 and Table S4. It is clearly evident that the undrained shear strength corresponds to the fluidic yield stress and the Bingham yield stress matches with the difference between the fluidic yield point and the static yield point. The fluidic yield stress obtained from stress ramp-down step also resembles with the Bingham yield stress, due to the fact that the sample was extensively sheared in the constant shear rate step performed before the stress rampdown step. However, this fact needs to be verified by systematic investigation of larger range of subsamples from PoH and mud samples from other sources.
The analysis of undisturbed consolidated mud samples (with higher densities) is not possible with Couette geometry due to the difficulty in experimental procedure (inclusion of the bob in the cup). Another problem, with consolidated samples, is the wall/sample interaction that can lead to unwanted slippage, which disqualifies the Couette method. However, Couette geometry with roughened surface or with grooves can be used to try to avoid this slippage problem. In order to estimate the fluidic yield stress of very dense mud samples, a pre-shear method with vane geometry can be used instead of Claeys protocol. The pre-shear method gives similar results as the Claeys protocol (Fig. 5a) but the experimental time is very much reduced (i.e., less than 5 min as compared to ~ 20 min).

Optimization of stress ramp-up test
In order to optimize the experimental time of stress ramp-up protocol for mud samples of different consistencies/densities, stress ramp-up tests with varying sweep rates (0.1-10 Pa/s) were performed for different mud layers (i.e., samples L6-1, L6-2 and L6-4) using Couette geometry. The results are shown in Fig. 7a-c. It is clear from Fig. 7a that for the sample with lower densities (i.e., fluid mud), slower sweep rates (0.1-0.2 Pa/s) are more appropriate to obtain both viscosity declines, since the transition for static yield stress lies in the range of very low shear stress. For the consolidated sample (with higher density), higher stress ramp-up rates (3-10 Pa/s) are possible to determine the static and fluidic yield stress values (Fig. 7c). The preferred stress ramp-up rates for different mud layers along with their approximate experimental times are presented in Table 4. Based on this analysis, it is obvious that stress ramp-up test is very fast, reliable and robust for measuring the yield stresses (both static and fluidic) of mud samples.

Empirical fitting of flow curves with existing models
The results presented so far (see for instance Fig. 4) show that the shear rate/shear stress function display a two-step yielding behaviour for mud. This feature has already been reported for mud samples from other locations (Huang and Aode, 2009;Nie et al., 2020;Yang and Yu, 2018). In order to fit the shear rate/shear stress function, these authors either used a Dual-Bingham model or a Worrall-Tuliani model. Due to its formulation, see Eq. (3), the Worrall-Tuliani model can only be used to fit the higher shear rate region quite effectively. The Dual Bingham model or Dual Herschel-Bulkley model can capture the characteristic features of two-step yielding by fitting the basic model (i.e., Bingham or Herschel-Bulkley) to two different regions of the same flow curve. However, as shown in supplementary information (Fig. S5), it requires to define a critical shear rate which delimit each of the yielding regions (See Sec. 3 of supplementary information for detailed information about empirical fitting with existing models). Therefore, for practical purpose, it is quite useful to develop a model that can capture both yielding steps Fig. 6. (a) Apparent viscosity as a function of shear stress for mud sample L6-3 using stress-ramp up test with Couette geometry; (b) shear stress as a function of shear rate for mud sample using Claeys protocol with vane geometry; solid lines are just the guide for the eye.
in a single equation.

Proposed model for two-step yielding
Though the Dual-Bingham and Dual-Herschel-Bulkley models can be adapted to fit each of the two yielding regions of flow curves, the selection of the critical shear rate which delimits these regions is quite difficult and time consuming to find, particularly in case of significant amount of experimental dataset. Therefore, we derived an empirical model (described in Sec. 3.2) which can be used to fit the data over the whole range of shear rates.
The experimental data of stress ramp-up tests, performed for the mud samples obtained from different locations/depths of port (Fig. 1), were fitted with our model, as shown in Fig. 8a for one mud sample. The model displays a good agreement with the experimental data and the values of the fitting parameters are also shown in Fig. 8a for that sample. It is already known that the mud samples obtained from different locations of port have quite different rheological fingerprint due to the variation in organic matter content, sand content, etc. (Shakeel et al., 2019(Shakeel et al., , 2020b. Therefore, after fitting all the experimental data with the model, the fitting parameters were correlated with the density and organic matter content (TOC), in order to generalize the model for the whole port. Fig. 8b shows the fluidic yield stress (τ f ) as a function of excess density for the samples from all locations. The experimental data above the excess density (ρ − ρ w ) of 200 kg .m − 3 shows deviation from the linear behaviour on a semi-log scale. This deviation from the linear trend could be related to the existence of a plateau behaviour at high density values, i.e., yield stress becomes independent of density (within the large experimental deviation) after that critical value of density. Moreover, the mud samples with such a high density are not important for defining nautical bottom in ports and waterways. Consequently, the data above this density value was not considered further for the correlation between fitting parameters and density/TOC.
The experimental data of fluidic yield stress (τ f ) as a function of excess density (ρ − ρ w ) for different locations was fitted with a power law model, given as: where 'a' and 'b' are two fitting parameters. It is evident that the parameter 'b' of different lines is quite similar while the other parameter ('a') is significantly different for different locations, as shown in Fig. 9a and b. It can also be seen that the organic matter content (TOC) shows a decreasing trend from location L1 to L10 (from the river side to the sea side). As the parameter 'b' was not varying significantly, a fixed value of 2.4 was used for parameter 'b'. The power law fitting was again performed with just one fitting parameter 'a', as shown in Fig. 9c. This fixed value of the fitting parameter 'b' is slightly lower than the one reported in literature for mud samples by Nie et al. (2020) who claimed 'b' to be   within the range of 5.2-5.8, due to the higher density of their mud samples.
A similar approach was also used for the correlation between the static yield stress (τ s ) and excess density (see Fig. S6). The correlation between fitting parameter 'a' (both static and fluidic) and TOC for different locations is shown in Fig. 9d. It is clearly evident that a strong correlation exists between the fitting parameter 'a' and the TOC for different locations. A same approach was used for the correlation between infinite viscosity (μ ∞ ) and TOC where the fixed value of parameter 'b' was 1. The corresponding values for parameter 'a' is given in Fig. 10b. The correlation of the static shear rate γ s , as defined in Sec. 3.2, with excess density for location 'L3' is shown in Fig. 10a. For γ s , it can be seen that the parameter shows a decreasing trend with the excess density, except at very small densities, where a cluster of data is observed. This cluster of data was not included in the fitting because such low density samples either behave like a Newtonian fluid (i.e., have no yield stress) or exhibit single-step yielding, and therefore, fitting with our twostep yielding model can produce inaccurate values of the parameters.
Discarding this data, a power law fitting was applied for this parameter (γ s ) as well, by fixing the value of parameter 'b' equal to − 2.4 (Fig. 10a). Fig. 10b shows a strong correlation between the fitting parameter 'a' for different model parameters and the TOC for different locations. The remaining two model parameters, γ f and γ 0 showed no correlation with the density of the samples. These parameters were observed to vary within the range of 2-22 s − 1 for γ f and 0.1-13 s − 1 in case of γ 0 for all the studied locations. However, the values of these parameters may depend on the stress sweep rate, type of TOC, etc. which needs further investigation.
The upper and lower limits of the parameter 'a' for different model parameters were also estimated for all the locations, in order to cover the scattering of the data for each location (see Table S8). These limits of parameter 'a' were estimated by fitting the higher and lower values of yield stresses with the same power law. Fig. 11 shows the power law fitting line (black line) along with the upper and lower limit lines (red lines) for location 'L6' in case of fluidic yield stress. The samples having excess densities higher than 200 kg. m − 3 are also plotted in Fig. 11, but  these samples were not included in the empirical fitting. These samples also lie within the upper and lower limits. However, some of these samples show a clear deviation from the predicted increase, as already explained. The proposed model may not be suitable for other investigated protocols, particularly, where two-step yielding is not observed. For those protocols, simple Bingham model or Worrall-Tuliani model can be used to fit the experimental data.

Conclusions
Natural mud sediments can exhibit complex rheological behaviour like thixotropy, viscoelasticity and yield stress. In this study, rheological analysis was performed on the natural mud samples obtained from different locations of Port of Hamburg, Germany. Fast and easily repeatable measurements of yield stresses are important for practical applications. Therefore, different rheological protocols were compared with our recently developed, time-efficient, stress ramp-up rheological protocol to obtain the flow curves and yield stress values.
The results showed that "CSRT-ramp up", "EFC-increasing" and "stress ramp-up" displayed higher static and fluidic yield stress values as compared to the rest of the methods. All these methods started with a low shearing action, implying starting from an undisturbed state of the sample, which resulted in higher values of the yield stresses. These methods are, therefore, suited to measure the yield stresses of mud sediments for harbour applications, where mud sediments usually exist in undisturbed state (i.e., structured state). However, the EFC-increasing protocol is still a quite longer test and from CSRT-ramp up the determination of static yield stress is not very straightforward. Therefore, the stress ramp-up protocol is the most suited and fastest method to analyse the yield stresses of mud sediments for harbour applications. Furthermore, the yield stress values from stress ramp-up test corresponds well with the conventional yield stress values obtained from Claeys protocol. The optimization of stress ramp-up test enabled to reduce the experimental time for different mud layers (∼ 10-200 s).
The important result from the empirical fitting was that the fluidic yield stress corresponds to the undrained shear strength while the difference between fluidic and static yield stresses corresponds to the Bingham yield stress. However, this comparison between different yield points needs further investigation. Our model was quite accurate in capturing the two-step yielding behaviour of mud samples from Port of Hamburg, within the density range of 1050-1200 kg. m − 3 . The applicability of stress ramp-up protocol and our empirical model need further investigation for the mud sediments from different ports.