As one of the most productive ecosystems on earth, mangroves provide a variety of ecosystem services with environmental and economic benefits (Nellemann et al., 2009; Barbier et al., 2011; de Groot et al., 2012). Mangroves can protect coastal communities from storm surge events by dissipating energy from waves and currents with their above-ground biomass (e.g., branches, leaves, and aerial roots) (Mazda et al., 1997; Mazda et al., 2006; Vo-Luong and Massel, 2008; Horstman et al., 2014). Energy dissipation by mangrove forests also creates shelter for many aquatic species and supports fisheries, which provide jobs and food for millions of people (Barbier et al., 2011; Hutchison et al., 2014).

On a global scale, mangroves provide an important mechanism for mitigating climate change by trapping and sequestering carbon-rich sediment in their soils (Mcleod et al., 2011; Twilley et al., 2017; Kauffman et al., 2020) at a rate of 200 g C m^-2 year^-1 (Temmink et al., 2022). Despite occupying only 0.5% of the global coastal area, mangroves store 10-15% of total coastal carbon (Alongi, 2014) and have a carbon density of 900 Mg C ha^-1 (Temmink et al., 2022). The mangrove carbon budget is comprised of carbon-rich sediment from autochthonous (i.e., produced in situ by the mangrove) and allochthonous (i.e., produced outside the forest) sources (Woodroffe et al., 2016). Allochthonous sediment enters the forest via tidal inundation or storm surge events, and it is the ability of mangroves to trap this sediment that makes them such a significant carbon sink (Jennerjahn and Ittekkot, 2002; Adame and Lovelock, 2011; Woodroffe et al., 2016).

Above-ground biomass, including mangrove pneumatophores (i.e., vertical aerial root structures, Figure 1 ), creates conditions that facilitate deposition by enhancing drag and slowing currents near the bed (Furukawa and Wolanski, 1996; Horstman et al., 2017; Mullarney et al., 2017b). However, mangrove pneumatophores also generate root-scale turbulence that enhances turbulent kinetic energy (TKE), which can promote sediment resuspension and lead to erosion (Mullarney et al., 2017a; Norris et al., 2017; Norris et al., 2019; Norris et al., 2021). Because of the competing effects of velocity reduction and turbulence enhancement, the relationship between vegetation density and sediment stability is not straightforward (Fagherazzi et al., 2017; Mullarney et al., 2017a; Xu et al., 2022a). Understanding how pneumatophore roots impact the balance of the competing processes of deposition and erosion is critical for improving the assessment of sediment retention and carbon storage in mangrove forests.

The rate of net deposition ( d m d t ) can be described in terms of a deposition probability,

in which m is the net mass deposited per bed area over time t , p is the probability that particles reaching the bed will remain deposited, w s is the settling velocity, and C is the near-bed suspended sediment concentration in a vertically mixed system. In Equation 1, the probability p captures the influence of resuspension on mass accumulation. In the absence of resuspension, there is pure deposition ( p = 1 ). When resuspension is present, p < 1 . Engelund and Fredsøe (1976) developed a model to predict the probability ( p ' ) that a particle on the bed is put in motion by the bed shear stress, and Zong and Nepf (2010) used this to describe the probability that a particle remains at the bed, p = 1 − p ' .

β is a bed friction coefficient, which we set to β = 1 . θ is the dimensionless shear stress (Shields parameter),

in which τ b is the bed shear stress, ρ s is the sediment density, ρ is the water density, g is the gravitational acceleration, and d s is the particle diameter. The Shields parameter is a ratio of destabilizing (time-mean stress) and stabilizing (grain weight) forces acting on a single grain. The critical Shields parameter ( θ c ) is defined by the critical bed shear stress ( τ b , c ) needed to initiate sediment motion. When θ < θ c , p = 1 , indicating pure deposition. When, θ ≤ θ_c   p < 1 , indicating the presence of resuspension.

In vegetated systems, vegetation-generated turbulence enhances resuspension by two means: (1) mixing momentum toward the bed, which enhances τ b (Liu et al., 2008; Conde-Frias et al., 2023) and (2) directly interacting with the bed and mobilizing sediment with enhanced instantaneous shear and normal stress (e.g., Xu et al., 2022b). Many previous studies have described the importance of instantaneous forces (both lift and stress) associated with turbulence in mobilizing sediment grains (e.g., Bagnold, 1941; Nino and Garcia, 1996; Zanke, 2003; Smart and Habersack, 2007; Diplas et al., 2008). Sediment transport models written in terms of bed shear stress (Equation 2) have yielded inaccurate predictions for vegetated systems, because they do not account for vegetation-generated turbulence (Yang et al., 2016; Tinoco and Coco, 2018; Yang and Nepf, 2018; Liu et al., 2022). Yang et al. (2016) and Tinoco and Coco (2018) found that near-bed TKE ( k t ( nb ) ) is a better predictor of sediment motion in vegetated systems. Therefore, we hypothesized that Equation 2 might better predict deposition in vegetated systems if it were recast in terms of ( k t ( nb ) ) ,

In bare channels, near-bed TKE is generated by the bed shear, such that bed shear stress and TKE are linearly related ( τ b = ρ ω k t ( nb ) ,   with ω = 0.2 , Soulsby, 1981). This relation suggests a method for redefining the critical Shields parameter (Equation 3) in terms of TKE,

with θ k t , c defined by the critical near-bed TKE ( k t , c ) needed to initiate sediment resuspension (Zhao and Nepf, 2021; Liu et al., 2022). Rewriting the Shields parameter in terms of TKE respects the original physical meaning, but expands the understanding of the destabilizing forces to include the effects of turbulence (e.g., Tinoco and Coco, 2018). Because the critical level of turbulence is the same in vegetated and bare channels (Yang et al., 2016), the critical turbulence level can be inferred from bare bed conditions. Specifically, k t , c = τ b , c / ρ ω , such that θ k t , c = θ c .

Liu et al. (2022) used Equation 4 to predict deposition within a model canopy of Phragmites australis, which has a morphology consisting of a central stem surrounded by multiple leaves. Good agreement was achieved in Case 4 (Figure 8 in Liu et al., 2022), but agreement was not as good for other cases (Figure 5 in Liu et al., 2022). The robustness of the deposition model was not discussed in a systematic way across flow conditions. Further, no bare bed conditions were examined. In contrast, the present study systematically considered paired vegetated and bare bed conditions across the same range of velocity, and also extended to higher values of solid volume fraction. This facilitated a more detailed description of the parameter range over which the model may be successfully applied. Advancing existing deposition models will help to improve the modeling of sediment transport in mangrove systems and facilitate the assessment of sediment and carbon retention.

Experiments were conducted in a recirculating Plexiglas flume with a 283 cm x 20 cm x 39 cm working section (dashed black outline in Figure 2 ). Plexiglas inserts (gray in Figure 2 ) were used to constrict the test section width to 20 cm. The water depth measured at the downstream end of the test section was H = 10   cm , which is within depth ranges typically observed in mangrove forests (Furukawa and Wolanski, 1996; Norris et al., 2021). A sharp-crested weir ( h = 5   cm ) located at the downstream end of the flume was used to fix the water depth.

Rigid vegetation, like pneumatophores, has been modeled using cylindrical dowels in several laboratory studies (Zong and Nepf, 2010; Yang et al., 2016; Tinoco and Coco, 2018; Liu et al., 2021). In this study, 0.8-cm diameter ( d ) PVC dowels were used to represent pneumatophores. The dowels were screwed into a PVC board that was inverted and inserted downward into the channel until the dowels just touched the bed. This allowed the dowels to be easily removed at the end of the experiment without disturbing the deposited sediment ( Figure 2 ). In the field, pneumatophore canopies are spatially heterogeneous, and the pneumatophores vary in diameter (0.5 to 2 cm), height (1 to 30 cm), and solid volume fraction ( ϕ = π 4 n d 2 = 0.005 to 0.04, in which n is roots per bed area) (Tomlinson, 2016; Yando et al., 2016; Norris et al., 2017; Norris et al., 2021). Three solid volume fractions were considered in this study: ϕ = 0.01, 0.02, and 0.04. The positions of the dowels within the dowel array boards were determined using a random array generator code (MATLAB).

To characterize the flow field, a Nortek Vectrino recorded instantaneous velocity components in the streamwise ( u ( t ) ), lateral ( v ( t ) ), and vertical ( w ( t ) ) directions in both the bare and vegetated test sections. Four channel-averaged velocities ( U o ) were considered ( Table 1 ). These velocities spanned a typical range of flow conditions observed within mangrove forests (Furukawa and Wolanski, 1996; Norris et al., 2021). Velocity was measured at multiple locations across the flume at z = 0.5 cm (near-bed) elevation for each channel-averaged velocity. Due to the short length of the flume, the flow was not fully developed in the bare test section, so that both wall- and bed-boundary layers were small compared to the flume width and depth. However, due to the channel constriction ( Figure 2 ), the velocity changed 10 to 20% across the bare test-section width. This variation was captured by a five-point lateral profile in the bare section. Within the dowel array, the velocity varied at the scale of the dowel, but the laterally-averaged conditions were fully developed after just a few cylinder rows. To capture the spatial heterogeneity in the dowel array, the lateral profile included 15 positions. At each position, the velocity was measured at 200 Hz for 60 s. Tests with longer records confirmed that 60 s was sufficient to capture the mean and turbulent velocity statistics.

To estimate depth-averaged velocity, a profile was constructed from measurements at 0.5-cm increments from z = 0.5 to 4.5 cm at a lateral position that was closest to the laterally averaged near-bed velocity. The vertical profile in the bare test section was used to calculate U o . Each velocity record was decomposed into time-averaged ( u ¯ , v ¯ , w ¯ ) and fluctuating ( u ′ ( t ) , v ′ ( t ) ,   w ′ ( t ) ) components and processed using the Goring and Nikora (2002) method to remove spikes. The acceleration and velocity variance threshold parameters for this method were set to λ = 1 and k = 1.5 , respectively. Turbulent kinetic energy per fluid mass is k t = 1 2 ( u ' 2 ¯ +   v ' 2 ¯ +   w ' 2 ¯ ) . Near-bed laterally averaged TKE was calculated for the bare ( k t , b ( nb ) ) and vegetated ( k t , v ( nb ) ) test sections using measurements at z = 0.5 cm ( Table 1 ). Velocity measurements were made separately from the deposition experiments to avoid disturbances to the water column that could impact deposition.

The experiments used solid glass spheres with diameter d s = 11   μm and density ρ s = 2500 kg/m³, which were selected based on grain size measured at field sites in a black mangrove forest ( O ( 10 μm ) ¹ ). Deposition experiments began by weighing 16.3 g of glass spheres, adding them to a 1 L container with water and surfactant (Windex^® Original Glass Cleaner was added to help the sediment slurry mix with the water), and shaking the container vigorously. This mass of glass spheres was chosen to achieve an initial concentration of C ≈ 20 mg/L throughout the flume, which is within the range of suspended sediment concentrations observed in mangrove forests (Furukawa et al., 1997; Horstman et al., 2017). The non-cohesive sediment mixture was poured across the width of the tail tank, and the recirculating pump mixed the sediment and water into a uniform concentration. Each deposition experiment ran for 4 hrs. The methodology for these deposition experiments was adapted from Zong and Nepf (2010) and Liu et al. (2022).

An optical backscatter sensor (OBS, Seapoint Sensors, Inc.) was used to measure the evolution of C over the duration of the experiment ( Supplementary Section 1 ). The OBS (20 Hz sampling rate) was located at the upstream end of the first bare baseboard and positioned at mid-depth. Preliminary studies confirmed that throughout an experiment C was the same at the upstream and downstream end of the flume, so only one OBS was needed to measure C . This reflected the fact that the time-scale over which deposition occurred (hours) was much longer than the time-scale of mixing (minutes), with complete mixing occurring each time water passed through the pumps. Therefore, the concentration remained uniform in the test section, even as it declined due to deposition. The OBS output voltage was calibrated using prepared concentrations ranging from 0 to 44 mg/L ( Supplementary Section 2 ).

After 4 hrs, the flume was left to slowly drain for 1 hr. Using the methodology discussed in Zhang et al. (2020), we found that the flume draining period had a negligible impact on the deposition pattern, and negligible additional deposition occurred during this time, consistent with the low C at the end of the experiment. The baseboards were left to dry in the flume for 1 day. Once the baseboards were dry, the dowels were carefully lifted off the baseboards, and the bare and vegetated test section baseboards (red outline in Figure 2 ) were carefully removed from the flume with gloves. An acetate template was placed over each board and secured with clips. This template divided the board into three 15 cm x 15 cm windows.

Three glass fiber filters (0.7-μm pore size, 47-mm diameter) were weighed in advance, lightly wet with water, and then used to wipe the sediment off the baseboard within each window (9 filters total). Tests with additional filters indicated that using three filters was sufficient. Adding a fourth filter increased the mass by only 6%. The filters were dried in a 60 o C oven for 4 hrs, which was sufficient for the filters to reach a constant weight. After drying, the filters were reweighed, and the average net deposition per bed area for the bare ( m bare ) and vegetated ( m veg ) test sections was calculated. The uncertainty in mass deposition predominantly came from the variation among the three windows of each test section.

The deposition probability was estimated by rearranging and integrating Equation 1 over the experiment duration, T ,

Using Equation 6, the deposition probability in the bare ( p bare ) and vegetated ( p veg ) test sections were estimated using the mass deposited in each section, m bare and m veg , respectively. Based on Equation 1, C ( t ) should follow an exponential decay with rate constant − p w s . So, it was reasonable to smooth the concentration record by fitting the form C = C o e − b t with initial concentration C o and constant b . Because the flume was a closed system, the temporal change in C was due only to deposition. The fitted concentration record was used in Equation 6 ( Supplementary Section 3 ).

The bed shear stress in the vegetated region ( τ b , v ) was estimated using Equation 7, developed by Conde-Frias et al. (2023), that describes the enhancement of bed shear stress by turbulence generated from rigid, emergent vegetation,

in which K = 9.5 ± 0.4   is a scale constant. R e d = U o d ν ( 1 − ϕ ) is the stem Reynolds number, in which ν is the kinematic viscosity of water. C f = 0.002 is the bed friction coefficient ( Supplementary Section 4 ). In the bare test section, τ b = ρ C f U o 2   was used to estimate bed shear stress. Conde-Frias et al. (2023) validated Equation 7 against data and simulations with ϕ = 0.016   to  0.25   and R e p up to 1300 (Table 2 in Conde-Frias et al. (2023)), which spans similar conditions examined in this study. Uncertainty was propagated for all calculations using the constant odds combination method described in Kline and McClintock (1953).

The critical Shields parameter was used to estimate the critical bed shear stress and critical turbulence threshold for resuspension (Zhao and Nepf, 2021). From Julien (2010),

in which ϕ R is the angle of repose. For the d s used in this study, ϕ R is 30 o . d * is the dimensionless particle diameter,

Equation 8 applies for d * = 0.3   to  1.9 . Using the parameters from this study, the critical bed shear stress and critical near-bed TKE were

In the bare test section, k t , b ( nb ) < k t , c (and τ b <   τ b , c ) for the lower velocity Cases 1 and 2 ( Table 1 ). Therefore, the bare test section was considered to be purely depositional ( p = 1 ) for these cases, from which Equation 6 can be used to estimate the settling velocity, w s ¯ =   0.0038 ± 0.0003 cm/s. This was consistent with Zong and Nepf (2010), who found w s = 0.004 ± 0.002 cm/s for solid glass spheres from the same manufacturer and of the same d and ρ s as used in this study. The estimated settling velocity was in reasonable agreement with Stokes’ Law (Stokes, 1851), w s = g ( ρ s ρ − 1 ) d s 2 18 ν = 0.010 cm/s, given that the manufacturer’s specifications included a range of diameters (10% finer: d s = 3 μm; 90% finer: 15 μm). Using Stokes’ Law, w s ¯ =   0.0038 cm/s suggested a mean diameter of d s = 7 μm. The value of w s ¯ estimated from Case 1 and 2 bare test section data was subsequently used in Equation 6 to solve for the deposition probability p in all other cases.

For the same channel-averaged velocity, the vegetated test section had elevated TKE compared to the bare test section ( Figure 3A ). As TKE increased, resuspension increased, which was reflected in lower net deposition in the vegetated test section relative to the bare test section at the same channel-averaged velocity ( Figure 3B ). At the lowest velocity, all three vegetation densities produced m veg / m bare   = 1 ( Figure 3B ). This was also observed for ϕ = 0.01   at the second velocity setting. For each of these cases, the near-bed turbulence was, within uncertainty, less than or equal to k t , c = 2.6 cm 2 / s 2 (marked with asterisks in Table 1 ), confirming the predicted value of k t , c (Equation 11).

For both vegetated and bare test conditions, when k t ( nb ) ≤   k t , c = 2.6 cm²/s², the deposition probability ( p ) (Equation 6) was 1 within uncertainty, and when k t ( nb ) >   k t , c = 2.6 cm²/s², p   decreased with increasing k t ( nb ) ( Figure 4A ). Furthermore, test cases with bare beds (open circles) and arrays of different solid volume fractions (triangles) collapsed to the same trend when plotted versus k t ( nb ) , consistent with the TKE adaptation of the Engelund and Fredsøe (1976) model (Equation 4, dashed black line in Figure 4A ). In contrast, although p generally decreased with increasing τ b , the trends were different between bare bed and vegetated conditions ( Figure 4B ). In the bare test section, the τ b model (Equation 2, solid black line in Figure 4B ) predicted a p value consistent with measurements (open circles). However, in the vegetated test section, the τ b model (Equation 2 with τ b , v predicted by Equation 7) overpredicted p by as much as 6-fold, because it failed to fully account for the impact of the vegetation-generated turbulence. The TKE model did best for low k t ( nb ) but underpredicted p for high k t ( nb ) ( Figure 4A ).

The original contributions presented in the study are included in the article/ Supplementary Material. Further inquiries can be directed to the corresponding author.

AD: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. EH: Formal Analysis, Investigation, Writing – review & editing. DR: Conceptualization, Formal Analysis, Supervision, Writing – original draft, Writing – review & editing. HN: Conceptualization, Formal Analysis, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing.