Comparison of Methods for Determining Erosion Threshold of Cohesive Sediments Using a Microcosm System

Erosion of cohesive sediments is a ubiquitous phenomenon in estuarine and intertidal environments. Several methods have been proposed to determine the surface erosion threshold (τc0), which are still debatable because of the numerous and uncertain definitions. Based on erosion microcosm experiments, we have compared three different methods using (1) eroded mass (EM), (2) erosion rate (ER), and (3) suspended sediment concentration (SSC), and suggested a suitable method for revealing the variation of erodibility in intertidal sediments. Erosion experiments using a microcosm system were carried out in the Muuido tidal flat, west coast of South Korea. The mean values of τc0 for three methods were: 0.20 ± 0.08 Pa (EM); 0.18 ± 0.07 Pa (ER); and (3) 0.17 ± 0.09 Pa (SSC). The SSC method yielded the lowest τc0, due to the outflow of suspended sediment from the erosion chamber of the microcosm. This was because SSC gradually decreased with time after depleting the erodible sediment at a given bed shear stress (τb). Therefore, the regression between SSC and applied τb might skew an x-intercept, resulting in the underestimation (or “not-determined”) of τc0. The EM method yielded robust and accurate (within the range of τb step at which erosion begins) results. The EM method represents how the erodible depth thickens as τb increases and therefore seems better suited than the SSC and ER methods for representing depth-limited erosion of cohesive sediments. Furthermore, this study identified the spatiotemporal variations of τc0 by EM method in an intertidal flat. The τc0 in mud flat was about two times higher than that in mixed flat. Compared to the end of tidal emersion, the sediment was 10–40% more erodible at the beginning stage.


INTRODUCTION
The sediment erodibility plays a key role in the evaluation of resistance to erosion, and is quantified as the eroded mass (EM) and erosion rate (ER). Understanding the erodibility of cohesive sediments is important for many ecological (e.g., primary production, benthicpelagic coupling, and toxicity problems) and engineering applications (e.g., siltation of harbors, dredging of navigation channels, and coastal protection). The erodibility of non-cohesive sediments is primarily determined by their physical properties (e.g., grain size and bulk density) (Roberts et al., 1998). As cohesive sediments have highly variable parameters in space and time, however, it is difficult to directly measure or predict them. Furthermore, complex interactions occur between the physical, biological, and chemical parameters that impact the erodibility of cohesive sediment (Black et al., 2002;Grabowski et al., 2011). Thus, the site-and time-specific field experiments are often required to obtain the information on the erodibility of cohesive sediments.
The surface erosion threshold (τ c0 , also known as an initial critical shear stress for erosion) is largely used as one of parameters representing for sediment erodibility. It is defined as the bed shear stress (τ b ) at which the sediment erosion begins . Currently, several methods have been proposed to estimate a specific threshold. Amos et al. (2003), for instance, proposed three practical methods to define τ c0 : (1) the surface intercept of the failure-envelope on a plot of eroded depth vs. τ b ; the extrapolation to the ambient level through the regression of (2) ER and (3) suspended sediment concentration (SSC) vs. τ b , respectively.
Despite such methodological efforts, the determination of τ c0 is still debatable because of the investigator's subjective definitions used to identify the initial sediment motion on the seabed (Sutherland et al., 1998). A portable and easy-to-use erosion microcosm system capable of deployment in the field has been developed. This type of device has been widely used in littoral environments, particularly intertidal flats, to quantify spatiotemporal variations of erodibility (Dickhudt et al., 2009;Wiberg et al., 2013;. To date, many researchers have visually defined a τ c0 based on the abrupt increase in SSC (or ER) (e.g., Maa and Kim, 2002), however, this approach is quite subjective. For example, τ c0 was defined as the τ b when an abrupt increase of SSC (>1.5 mg l −1 ) was first observed . Thus, there is a need to develop an accurate and objective method for determining τ c0 using the extrapolation method (Widdows et al., 2007;Seo et al., 2020).
A series of erosion experiments were conducted on the intertidal flat, the west coast of South Korea to meet the aforementioned technical needs in the cohesive sediment community. The collected data were used as input in the three extrapolation methods using EM (the amount of erodible sediment at the bed), ER (the mass loss per unit time and area), and SSC, respectively. The main objectives of this study were (1) to compare methods for estimating τ c0 , and (2) to suggest a method suitable for revealing the spatiotemporal variation of erodibility in intertidal environments.

Study Area
The study area is located in the Muuido tidal flat, the west coast of South Korea (Figure 1). In Gyeonggi Bay, extensive tidal flats comprising a total area of 873 km 2 have developed around the bedrock islands (Koh and Khim, 2014). The surface sediments consist of fine sand to silt with a mean grain size of 51.5 µm . Erosion experiments have been conducted at two sites: mud flat (GB01) and mixed flat (GB02) from September 11 to 17, 2020. Sediment cores were selectively collected at the neap (N1: September 11; and N2: September 13) and spring (S1: September 16; and S2: September 17) tides to examine the spring-neap tidal variation of erodibility (Figure 2A). The tides have a semi-diurnal regime, with mean tidal range of approximately 9 m (KMA, 2017), and tidal currents are ebb-dominated, reaching up to 1.5 m s −1 (Lee et al., 2013). Winds are dominated by the regional monsoon with mild south-southeasterly winds during summer and strong northnorthwesterly winds during winter (KMA, 2017). The SSC was observed ranging from about 100-1,000 mg l −1 at nearby tidal channel (Park and Lee, 2016).

Microcosm Experiments
Sediment erodibility was measured using a dual-core Gust Erosion Microcosm System (GEMS) designed by Gust and Müller (1997). A GEMS comprises an erosion chamber, electronic control box, turbidimeter (Hach, 2100AN), water pump, erosion motor, and rotating disk (Supplementary Figure 1). A series of increasing τ b is applied to the surficial of the sediment core by controlling both the rotation rate of disk and the water pumping rate. A conceptual diagram and a photograph of the GEMS could be found in  and Seo et al. (2020). Sediment cores for erosion microcosm experiments were collected using an acrylic cylinder at the beginning and end of tidal emersion. The emersion time, which indicates how long since the site emerged from tidal flooding, was determined by the pressure sensor (RBR SoloD) at GB02. This sensor was deployed at 0.1 m above bed. After the collection, sediment cores were carefully carried to the laboratory, and then installed in the GEMS. Measurement of erodibility was started within 2 h after collecting a sediment core, to minimize the effects of dewatering and self-weight consolidation. Prior to the beginning of the erosion experiment, the sediment surface within the core was positioned at 10 cm below the bottom unit of rotating disk. Seven stepwise sequences of τ b (0.01, 0.05, 0.1, 0.2, 0.3, 0.45, and 0.6 Pa) were applied to the top surface of core for approximately 10 min at each step. During the first 10 min of the experiment, a minimum τ b (i.e., 0.01 Pa) was applied to flush out pre-existing suspended sediments within the core and inlet/outlet tubes. While the rotating disk began applying a given stress, the seawater obtained from the field site was pumped into the erosion chamber. The effluent water containing the eroded or resuspended sediments was passed through a turbidimeter, which was continuously recorded in turbidity [in nephelometric turbidity units (NTU)]. Throughout the erosion experiments, the τ b and turbidity data were logged at 1 s interval. The effluent water from each stress step was sampled and then filtered through a 0.7 µm glass fiber filter (GF/F). The filters were oven dried at 105 • C for 24 h and then weighed to estimate the ground-truth SSC (in mg l −1 ). The turbidity was converted to SSC using a linear regression between recorded NTU and the SSC derived from filtering (R 2 = 0.92, not shown). This was then used to determine the EM at each τ b . An erosion formulation developed by Sanford and Maa (2001) and Sanford (2006) was used as follows: where m is the eroded mass (kg m −2 ), t is the elapsed time (s), ER is the erosion rate (kg m −2 s −1 ), M is the erosion rate parameter (kg m −2 s −1 Pa −1 ), τ c is the depth-varying erosion threshold (Pa), and λ is the rate of sediment depletion (s −1 ).

Estimation of Surface Erosion Thresholds
In microcosm experiments, sediment erosion can be determined by the responses of SSC under increasing τ b , which is expressed as the EM and ER. The τ c0 , a point of initial erosion of the bed, was estimated by three methods: (1) EM method: τ c0 was determined by an x-intercept of the regression of the eroded mass vs. τ b ; (2) ER method: τ c0 was determined, at background level of ER, through the linear regression of the time-averaged erosion rate vs. τ b ; and (3) SSC method: τ c0 was determined, at background level of SSC, through the linear regression of the time-averaged SSC vs. τ b . All the three have been used to compare the τ c0 from various erosion devices (e.g., Neumeier et al., 2007;Widdows et al., 2007;Thompson et al., 2013;Seo et al., 2020) (Supplementary Table 1). In case of the ER and SSC methods, using time-averaged erosion rate and SSC from each τ b , the linear regression was computed to remove the horizontal heterogeneity in τ b and erosion parameters (Schoellhamer et al., 2017). Many researchers (e.g., Sanford, 2006;Jacobs et al., 2011) have interpreted the erosion rate for each applied τ b as the average over the time interval. Detail explanations of each method are given in Widdows et al. (1998) and Amos et al. (2003).  The type of erosion was determined from the trends in erosion rate through time (Amos et al., 1992). Depth-limited erosion (type I), which is characterized by an increase in bed resistance and a decrease in erodibility with depth, is identified by an erosion rate that peaks rapidly when τ b is increased, but then declines exponentially with time. Steady-state erosion (type II) is identified as having a near-constant erosion rate for a given τ b because the τ c does not change with depth in the sediment bed. Erosion of type I/II is a transitional form between type I and type II. The process of type I began with entrainment of the organicrich fluff layer (type Ia), which then gave way to erosion of surface bed materials (type Ib) (Amos et al., 2003). Type Ib is considered more significant threshold concerning bed erosion (Widdows et al., 2007). The erosion amount of the fluff layer was negligibly small because it was a thin layer; and its threshold for erosion was small because this material was newly deposited and did not have sufficient time to develop self-weight consolidation (Ha and Maa, 2009). In this study, therefore, τ c0 was determined by the onset of type Ib. If the linear regression between the erosion variables (EM, erosion rate, and SSC) and τ b is divided into two parts by an inflection point, thereby transitioning from type Ib to type I/II (Seo et al., 2020). This partitioning of erosion type would be better suited for threshold estimation than no partitioning of this. Figure 2B shows the time series of SSC under the applied τ b 's. At GB01, SSC decreased with time until τ b = 0.05 Pa, except the sediment core at the beginning of the emersion during S2 period. At GB02, SSC increased slightly in the transition period at which τ b increased stepwise from 0.01-0.05 Pa, and then it quickly decreased. This small amount of erosion is representative of the removal of the fluff layer. The highest SSC generated among all the sediment cores was 24.45 mg l −1 at GB01, and it was 170.18 mg l −1 at GB02, indicating the presence of more erodible sediments at GB02. In GB01 and GB02, there was initial erosion of sediment at the τ b of 0.3-0.45 Pa and 0.1-0.2 Pa, respectively. τ c0 was derived using linear regression between the erosion variables and τ b . Representative cases of the three methods for estimating τ c0 are shown in Figures 3, 4. For example, the time series of SSC at each τ b in GB02 during N2 period (see Figure 2B) was converted to the data in Figure 3F indicated by the light red squares. The regression of ER and SSC method was fit through the time-averaged ER and SSC at each τ b , respectively. At GB01, FIGURE 5 | Relationship between eroded mass and erosion threshold (τ c ) in GB01 (green circles) and GB02 (red squares). Data at other sites are from Sanford (2006); Dickhudt et al. (2009), Wiberg et al. (2013), and Xu et al. (2014Xu et al. ( , 2016. See Table 2 for site information in legend.

Comparison of Different Surface Erosion Thresholds
Frontiers in Marine Science | www.frontiersin.org the relationship between erosion variables and τ b showed a linear regression without an inflection point corresponding to type Ib (Figures 3A,C,E). At GB02, the relationship between the EM and τ b revealed two phases (except the sediment cores at the end of emersion during N1 and S2): the first was a region of steady increase in EM with increasing τ b from 0.1-0.3 Pa; the second was a distinct inflection in the regression line at the τ b between 0.3 and 0.45 Pa (Figure 3B). These phases might correspond to type Ib and type I/II, respectively. There was also an inflection in the regression line between the ER (and SSC) and τ b (Figures 3D,F). At the end of emersion during N1 and S2 periods, a linear regression between the EM and τ b reflected type Ib (i.e., a single regression line without inflection) in GB02 (Figure 4B).
A summary of τ c0 is provided in Table 1. At GB01, the mean values of τ c0 for three methods were: (1) 0.28 ± 0.03 Pa (EM); (2) 0.21 ± 0.06 Pa (ER); and (3) 0.20 ± 0.10 Pa (SSC). The τ c0 in ER and SSC methods were similar. The τ c0 by EM method was consistently higher than the results from ER and SSC methods, due to the steeper slope in the regression line (Figures 4A,C,E). At GB02, the τ c0 by ER and SSC methods partly yielded a negative x-intercept caused by the extremely low slope in regression (Figures 4D,F), resulting in the "not-determined (N.D.) τ c0 " ( Table 1). During N2 and S1 periods, the τ c0 by EM method was about 25% higher than the results from ER and SSC methods. During S2 period, it was about 10% lower than the τ c0 by ER and SSC methods. Because the EM representing the erodible depth was accumulated with increasing τ b , the EM method is suitable for reflecting the depth-limited erosion of cohesive sediments. After the erodible sediment was removed at each step of τ b , SSC and ER gradually decreased with time. Thus, the regression of the SSC and ER vs. τ b might skew an x-intercept. Especially, the SSC method yielded the lowest τ c0 due to the outflow water containing eroded sediments from the erosion chamber.
The data from the GEMS experiments could be represented as the plots of the EM vs. τ c (Figure 5). Based on the data from all sediment cores at each site, the equation was obtained using a best-fit power law regression in the form of m = a(τ c -τ c0 ) b , where m is the eroded mass (kg m −2 ) (Sanford and Maa, 2001;Dickhudt et al., 2009). This regression equation could be used in sediment transport models (Fall et al., 2014). The results of other published erodibility data in tidal flat/marsh environments FIGURE 6 | The surface erosion threshold (τ c0 ) at GB01 (upper panel) and GB02 (lower panel) sites during N1, N2, S1, and S2 periods (see Figure 2A). measured by GEMS were similar to the results presented here from the Muuido tidal flat (Figure 5). Based on the best-fit curve, the τ c0 in GB02 was determined as 0.03 Pa, which corresponds to an erosion threshold of fluff layer (i.e., type Ia) (Jago et al., 2002;Ha and Maa, 2009). However, the τ c0 in GB01 could not be determined because of a negative x-intercept of the regression curve caused by low erodibility. In this case, the previous studies (see references in Table 2) suggested either 0.01 Pa (the minimum τ b applied for GEMS) or 0 Pa for τ c0 to ensure that the fit produced reasonable data.

Spatiotemporal Variation in Erodibility
The erodibility of cohesive sediment has been quantified as the τ c0 and EM (Dickhudt et al., 2009;. Our results indicated a spatial variation of erodibility between mud flat and mixed flat. The τ c0 by EM method in GB01 (mud flat, mean = 0.28 ± 0.03 Pa) was about two times higher than that in GB02 (mixed flat, mean = 0.12 ± 0.03 Pa). The EM in GB01 (mean = 4.89 ± 3.80 g m −2 ) also was about one order lower than that in GB02 (mean = 47.81 ± 18.03 g m −2 ). As seen in the erosion curves (Figure 5), the mixed flat was more erodible than the mud flat. Previous studies suggested that the spatial variations of erodibility of intertidal flats were caused by two main factors: (1) the sediment composition; and (2) the extracellular polymeric substances (EPS) concentration (Panagiotopoulos et al., 1997;Andersen et al., 2010). Even small amounts of mud content can create the cohesive attachment when mixed in non-cohesive sediment. Proportions of clay greater than 10% is sufficient for sediment stabilization (van Ledden et al., 2004).  found that the sediment erodibility in mud flat was about seven times lower than that in mixed flat due to increasing clay content on the Yeochari tidal flat, northern Gyeonggi Bay (for location, see Figure 1A). Therefore, the field-based measurements of spatial variability in erodibility and sediment composition are needed to understand the regional-scale sediment transport and morphology evolution. The EPS, which is secreted out of microphytobenthos (MPB), plays a key role in the "biostabilization" (Andersen et al., 2010). The EPS concentration and MPB biomass, which have a semiannual to annual variability, are at their lowest levels from September to November (Park et al., 2014). In this study, it seems unlikely that the sediment erodibility was affected by EPS. Temporal variation of sediment erodibility has been found on various timescales such as intra-tidal and spring-neap tidal cycles (Grabowski et al., 2011). As shown in Figure 6, the τ c0 at the beginning of tidal emersion was about 10-40% lower than that at the end, indicating that the subaerial emersion time might be a factor to determine the intra-tidal variation of τ c0 (Fagherazzi et al., 2017). In addition, the sediment in the spring tide was about 80% less erodible than in the neap tide at GB02 (Figure 6). As the spring tides have stronger currents compared to neap tides, more sediments are resuspended into the water column, and the consolidated (i.e., low erodibility) sediment beds are more exposed.

CONCLUSION
Three methods for determining the τ c0 of cohesive sediments in an intertidal flat were compared using an erosion microcosm system. The conclusions drawn from this study could be summarized as follows: (1) The τ c0 by EM method was slightly higher than the results from the ER and SSC methods. The EM method is more robust and accurate (within the range of τ b step at which erosion begins) than the ER and SSC methods. The τ c of the consolidated sediment bed increased with sediment depth, and thus the eroded mass representing the erodible depth was accumulated with increasing τ b . This suggested that the EM method is suitable for reflecting the depth-limited erosion of cohesive sediments.
(2) The SSC method produced the lowest τ c0 , owing to the outflow of suspended sediment from the erosion chamber. After depleting the erodible sediment at each step of τ b , the SSC gradually decreased with time. Therefore, the regression between SSC and applied τ b might skew an x-intercept, leading to a lower τ c0 . (3) This study identified the spatiotemporal variations of sediment erodibility in an intertidal flat. The τ c0 in mud flat was about two times higher than that in mixed flat. Compared to the end of tidal emersion, the sediment was 10-40% more erodible at the beginning stage. The τ c0 during spring tide was about 80% higher than that during neap tide.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
HJH designed the study, performed the experiments, and wrote the manuscript. HKH identified research goals, performed the analysis, and supervised the study. Both authors contributed to the article and approved the submitted version.

FUNDING
This study is supported by the National Research Foundation of Korea (NRF-2018R1D1A1A02085804). This study is also supported by the project entitled "Marine Ecosystem-Based Analysis and Decision-Making Support System Development for Marine Spatial Planning" (20170325), and "Development of Advanced Science and Technology for Marine Environmental Impact Assessment" (20210427), funded by the Ministry of Oceans and Fisheries of Korea (MOF).