Biophysical Modeling of Mangrove Seedling Establishment and Survival Across an Elevation Gradient With Forest Zones

Mangrove forest development critically depends on the establishment and survival of seedlings. Mechanistic insights into how water levels, waves and bed level dynamics influence the establishment process of individual mangrove seedlings are increasing. However, little is known about how spatial and temporal changes in water levels, waves and bed level dynamics across elevation gradients in mangrove forests facilitate/limit seedling dynamics. For this study, a new seedling establishment and growth model was integrated into a process‐based hydrodynamic and morphodynamic numerical model. This biophysical model was applied to a fringing mangrove forest located in the southern Firth of Thames, Aotearoa, New Zealand. This study quantifies the increasing establishment density and survival probability of mangrove seedlings from the lower‐elevated unvegetated intertidal flat toward the higher‐elevated mature mangrove forest. Three cross‐shore zones with distinctive seedling dynamics were identified: (a) a zone with daily tidal inundation where seedling dynamics are episodic and limited by the dispersal of individual propagules that rapidly anchor to the substrate by root growth, (b) a zone with daily to bi‐weekly tidal inundation where seedling dynamics respond to variations in spring‐neap tidal cycles and, (c) a zone with less than bi‐weekly inundation where seedling dynamics are governed by high propagule supply and seedling survival probability. The seedling establishment density and survival probability are dominated by annual extremes in tidal hydroperiod and bed shear stresses, respectively. The obtained parameterizations can be used to incorporate seedling dynamics in decadal‐timescale mangrove forest development models that are instrumental for mangrove management and restoration.


Introduction
Seedling establishment and survival are key to the development and persistence of mangrove forests (Friess et al., 2012;Krauss et al., 2008).Establishing seedlings enables mangrove forest recovery from storms and adaptation to environmental change.However, mangrove seedlings (stem height <0.5 m) are substantially more vulnerable to hydrodynamics (i.e., tidal flows and waves) and morphodynamics (i.e., sedimentation and erosion) than adult mangrove trees (Balke et al., 2013a(Balke et al., , 2013b)).Successful establishment and growth of mangrove seedlings depends on periods with environmental forcing conditions below species-specific tolerances (Balke et al., 2014).The water levels, waves and bed level dynamics that facilitate/limit seedling establishment and growth vary with cross-shore elevation.Moreover, mangrove trees influence these processes and may facilitate seedling establishment and survival (Huxham et al., 2010).While mechanistic knowledge on the establishment and survival tolerances of individual mangrove seedlings has increased over the past decade (Balke et al., 2013a(Balke et al., , 2013b;;Van Hespen et al., 2022), little is known about how seedling dynamics (i.e., the temporal variability of seedling numbers per area) vary with cross-shore gradients in elevation and vegetation in mangrove forests, and what these variations imply for mangrove forest dynamics.
Understanding seedling dynamics is of crucial importance for mangrove forest management and restoration projects.The restoration of mangrove forests that have been largely lost due to human activities (e.g., land reclamation, infrastructure, aquaculture) is rapidly gaining attention as a "Nature-based Solution" for climate change mitigation (e.g., through carbon sequestration) and adaptation (e.g., through improved coastal safety) (Leal & Spalding, 2022).However, many of the initial restoration attempts failed due to a lack of knowledge on the required environmental conditions for specific mangrove species to establish and grow (Ellison, 2000;Kodikara et al., 2017;Primavera & Esteban, 2008;Wodehouse & Rayment, 2019).The implementation and upscaling of mangrove forest restoration requires knowledge as well as unified methods for monitoring and modeling of seedling dynamics.The benefits of assessing and possibly creating the required environmental conditions for seedling establishment and growth are widely emphasized (Balke & Friess, 2016;Lewis, 2005;Van Bijsterveldt et al., 2022;Van Cuong et al., 2015;Winterwerp et al., 2020) as well as the relevance of monitoring and evaluation (Beeston et al., 2023;Lewis & Brown, 2014).However, mangrove forest development models with predictive skills on seedling establishment and survival are yet to be developed (Gijsman et al., 2021a;Van Hespen et al., 2023).
The process of mangrove seedling establishment starts with propagules falling from the parent tree (i.e., propagule abscission).Propagule abscission occurs during the fruiting season of mangrove trees.Mangrove propagules are buoyant and can be transported by tides and waves (Di Nitto et al., 2013;Tomlinson, 1986;Van der Stocken et al., 2019).After dispersal, propagules require a period free of tidal inundation to rapidly anchor their roots in the soil and establish as seedlings, a period referred to as Window of Opportunity 1 (WoO 1) (Balke et al., 2011).Once established, the seedlings require a period with limited hydrodynamic disturbance and bed elevation change to survive and grow stronger, a period referred to as Window of Opportunity 2 (WoO 2) (Balke et al., 2015).Inundation-free periods and the increase in survival tolerances of individual mangrove seedlings have been quantified in experimental studies for several mangrove species (Balke et al., 2011(Balke et al., , 2013a(Balke et al., , 2013b;;Van Hespen et al., 2022), and were used to explain the absence of mangrove seedling establishment on an unvegetated mudflat (Balke et al., 2015).Hu et al. (2015) showed that spatial variations in environmental forcing conditions also influence the seaward extent of marsh vegetation establishment.For mangroves, however, only limited attention has been given to quantifying the role of spatial variations in environmental forcing conditions for (a) the seaward extent of mangrove seedling establishment and survival and (b) mangrove seedling dynamics across an elevation gradient and with mangrove forest zonation (Gijsman et al., 2021b).
Even though the seedling establishment and survival process is mechanistically depending on environmental conditions, mangrove growth models typically include seedling establishment probabilistically (e.g., Berger & Hildenbrandt, 2000;Chen & Twilley, 1998).Mangrove growth models simulate tree growth based on interactions between individual neighboring trees (e.g., considering competition for space, nutrients and light), and seedlings are introduced individually with a pre-defined probability of occurrence.This probabilistic approximation is based on seedling dynamics being characterized by natural variability in for example, propagule abscission and propagule dispersal.The imposed probability of occurrence in mangrove growth models is often independent of the local environmental processes (e.g., Van Maanen et al., 2015) or calibrated with parameterizations based on the hydroperiod (P; the relative period of inundation per tide) (Bessely et al., 2023;Buffington et al., 2021) and/or the bed shear stress (τ; the surface shear stress induced by tidal currents and/or waves) (Xie et al., 2020(Xie et al., , 2022)).
The aim of this study is to gain new insights into the role of ecological and physical conditions that influence mangrove seedling dynamics, that is, (a) the natural variability in propagule abscission and seedling physiology, (b) the effects of variations in water levels, waves and sediment dynamics, and (c) the reduction of environmental disturbances by mangrove tree presence.We investigated seedling dynamics by performing (manipulated) simulations with a newly developed process-based biophysical numerical model that incorporates both the probabilistic characteristics of propagule abscission and seedling physiology (e.g., stem height, root length) as well as the hydrodynamic and morphodynamic processes that can impede their establishment and survival.The model was applied to the southern Firth of Thames estuary in Aotearoa New Zealand, where a mangrove forest comprising of Avicennia marina var.australasica (the only mangrove species in New Zealand) exhibits a distinct forest zonation.The model results were used to explore parameterizations between the hydroperiod, bed shear stress and seedling establishment and survival.

Study Site Characteristics and Field Monitoring Data
The Firth of Thames estuary (hereafter Firth) is located on the northeast coast of the North Island of Aotearoa New Zealand (Figure 1a).The estuary is mesotidal, wave-dominated and characterized by an extensive intertidal flat (Figure 1a; Swales et al., 2015).Similar to other estuaries in New Zealand, the Firth intertidal flat has been expanding due to increased sediment inputs from river runoff (Horstman et al., 2018).Since the early 1960s, Avicennia marina var.australasica mangroves have been colonizing the expanding mudflat.At present, the upper intertidal flat is colonized by an approximately 1 km wide fringing mangrove forest (Swales et al., 2015).The cross-shore forest characteristics are divided in distinct zones, from the most seaward bare intertidal mudflat, to an area where seedlings are establishing in front of the forest, the dynamic forest fringe with trees up to 4 m in height, the scrub forest with a denser tree cover of trees up to 2 m in height and, most landward, the relict forest fringe with sparser trees up to 5 m in height (Figure 1b; Swales et al., 2015Swales et al., , 2019;;Horstman et al., 2018).These forest zones are based on a cross-shore survey of profile elevation and vegetation characteristics from November 2016 (Horstman et al., 2018;Montgomery et al., 2018) and the mangrove forest edge was defined at the seaward front of the dynamic forest fringe (Figure 1b).
Mangrove colonization in the Firth primarily occurred during four major seedling establishment events in which the forest fringe expanded seaward by several tens of meters (Swales et al., 2015).These events occurred during El Niño years that resulted in lower water levels and a calmer wind and wave climate (Lovelock et al., 2010;Swales et al., 2007Swales et al., , 2015)).Other years were characterized by a lack of (major) forest expansion events because of wave and bed level disturbances in excess of the tolerances for seedling establishment (Balke et al., 2015;Swales et al., 2007).Since the last major seedling establishment event in the 1990s, however, the mangrove forest in the Firth has steadily expanded seaward (Figure 2).From aerial images, it was observed that the gradual expansion of the mangrove forest edge between 2007 and 2022 was about 50 m (i.e., averaging 3.3 m/year), with an increasing number of trees observed up to approximately 100 m seaward from the contiguous forest edge (Figure 2).In December 2022, field observations were performed showing that a substantial number of seedlings had established in the seedling establishment zone during the previous fruiting season (i.e., November 2021-February 2022) (Figure 3).Environmental forcing conditions in the Firth have been monitored extensively, with water levels recorded at the Tararu tidal gauge (Waikato Regional Council, site nr.1033.1)since 1992, and wind and temperature recorded at the NIWA climate station (station 38619) since 2010 (Figure 1a).This study used these field data as boundary conditions for the model scenarios.Moreover, this study used data from field measurement campaigns conducted in November 2016, from May to June 2017, from January 2020 to January 2021, and in December 2022 to setup, calibrate and validate the model (Figure 1a; Gijsman et al., 2024).

Model Description and Parameter Selection
The model developed for this study is process-based (i.e., mechanistic) and includes mangrove propagule and seedling dynamics (abscission, dispersal, establishment, growth and survival) based on interactions between hydrodynamics (i.e., tidal currents and waves), morphodynamics (e.g., sediment erosion and deposition) and mangrove tree characteristics (e.g., stem density, stem diameter and stem height).The seedling dynamics model is based on field and flume data presented in Balke et al. (2015) and the hydrodynamics and morphodynamics are resolved in Delft3D Flexible Mesh (Deltares, 2022a, 2022b).

Propagule Abscission and Dispersal and Seedling Establishment and Growth
Establishment and growth of the Avicennia marina var.australasica mangrove seedlings are simulated on an individual basis and include the processes of propagule abscission and dispersal, and seedling establishment, growth and survival.The seedling establishment process starts with the delivery mechanism of new propagules by the existing parent trees (i.e., abscission).Propagule abscission at the study site occurs during the fruiting season from November until February (Balke et al., 2015).The maximum probability of each tree releasing a propagule (P abscission ) is calibrated with calibration coefficient β.P abscission is assumed to vary throughout the fruiting season and approximated with a sinusoidal relationship: with t [s] the timestep during the fruiting season and T f [s] the duration of the fruiting season.
The released propagules have an obligate dispersal period (ODP) before the propagule can initiate root growth, and a buoyancy period (BP) during which the propagule remains buoyant ( Van der Stocken et al., 2019).Values of 7 days (ODP) and 15 days (BP) were adopted in this study based on the characteristics of Avicennia marina propagules (reviewed in Van der Stocken et al. (2019)).Propagules are dispersed by the tide through inundation during their BP.Tidal dispersal of propagules is imposed to be offshore directed because of the study focus on offshore propagule transport from the forest onto the bare intertidal flat, a process that has been observed in the Firth (Swales et al., 2007).This offshore transport of propagules corresponds to the ebb-dominant tidal flow characteristics measured in the seedling establishment zone and dynamic forest fringe (Gijsman et al., 2024;Lovett, 2017;Tablada Torres, 2020).The net distance of offshore propagule transport by each tide is equal to a calibration factor (D i ) multiplied by the maximum tidal inundation depth at the source location of each propagule.When a location does not inundate, the propagules remain at that location.Propagules stop moving when their BP ends.
Propagules start to anchor their roots in the soil after their ODP lapses and when they get stranded during or at the end of their BP.Seedling establishment (the propagule takes root and turns into a seedling) occurs when a cumulative required inundation-free period (T IF ) has occurred without the rooting propagule being disturbed by tides (i.e., WoO-1; Balke et al., 2014).The required T IF is propagule-dependent and is randomly selected from a normal (i.e., Gaussian) distribution.In this study, the normal distribution of the required inundation-free period was characterized by a mean μ T IF ) of 3 days and a standard deviation (σ T IF ) of 1 day (Balke et al., 2015).Upon initial establishment, the seedlings also obtain the following individual characteristics: (a) an initial shoot length (H i,0 ), (b) an initial root length (R i,0 ) and (c) a growth limitation/acceleration factor (M i ).In this study, these characteristics were randomly selected from normal distributions that were characterized by μ H i,0 = 2 cm and σ H i,0 = 0.5 cm for the initial shoot length, μ R i,0 = 1.4 cm and σ R i,0 = 1 cm for the initial root length, and μ M i = 1 and σ M i = 0.2 for the growth limitation/acceleration.These values are selected following previously reported Avicennia marina seedling characteristics from the Firth (Figures 3a-3c in Balke et al., 2015) and following Van Hespen et al. (2022).The initial shoot length was set equal to the propagule size.Please note that a negative initial root length (R i,0 ) parameterizes a delayed start of root growth/germination after stranding of the propagule.
After establishment, the seedling shoot length H i (Equation 2) follows a growth curve approximated by a powerlaw function.The shoot growth speed (G H ) was set equal to 3.425•10 4 m per day (0.03425 cm per day) following Balke et al. (2015).The root length also follows a power-law function (Equation 3).The root growth speed (G R ) was set equal to 3.266•10 3 m per day (0.3266 cm per day) following Balke et al. (2015).In this study, powers c 1 and c 2 were set to 0.75 to obtain realistic shoot/root lengths of 1 year old seedlings (Balke et al., 2015).A i is the age of a seedling in days in Equations 2 and 3: (2) Established seedlings are dislodged by hydrodynamic and morphodynamic disturbances when one (or more) of the following thresholds is (are) exceeded: (a) if the bed shear stress exceeds their bed shear stress tolerance (τ cr,i ), (b) if the erosion depth exceeds their critical erosion depth (E cr,i ), or (c) if the deposition exceeds their shoot height (H cr,i ).The survival tolerances of seedlings increase with seedling growth following Equations 4-6, respectively (Balke et al., 2015;Hu et al., 2015): In Equation 4, τ cr,i,0 [N/m 2 ] is the initial bed shear stress tolerance at establishment and G τ i [N/m 3 ] is a growth coefficient for the bed shear stress tolerance.τ cr,i,0 was used as a model calibration parameter and G τ i was set to 41.35 N/m 3 following Balke et al. (2015).In Equation 5, E cr,i,0 [m] is an initial erosion depth tolerance at establishment and G E i [ ] is a growth coefficient for the erosion depth tolerance.E cr,i,0 was set equal to 0.1108 m and G E i to 0.03508 following Balke et al. (2015).The temporal growth and the development of the survival tolerances of the seedlings are displayed in Figures 4a-4d.

Hydrodynamic Forcing and Morphodynamic Changes
The environmental forcing of hydrodynamics (water levels and waves) and morphodynamics (sediment erosion and deposition) are solved with Delft3D Flexible Mesh (DFM; Deltares, 2022a, 2022b).DFM solves wave dynamics by a coupling with a SWAN based wave model that computes wave properties based on the wave action balance with energy sources and sinks (Booij et al., 1999;Deltares, 2022c).
Wave energy losses due to the presence of vegetation (ϵ v ) are included through the equation of Mendez and Losada (2004), using the mangrove parameters stem density (N), stem diameter (D), stem height (H), and stem drag coefficient (C D ).Tidal flow dynamics are computed with the unsteady shallow water equations (continuity and momentum conservation; Deltares, 2022a).The loss in flow momentum is computed through vegetation drag, bed friction and turbulent fluctuations.The presence of the emergent mangrove stems is included through an additional flow resistance parameter in the momentum balance related to vegetation drag (M v ), depending on N and D, H and C D (following Baptist et al., 2007).A Manning roughness coefficient (n) was imposed as bed friction and set to 0.018 s/m 1/3 , following Gijsman et al. (2024).Momentum loss due to turbulent fluctuations is included through an artificially increased eddy viscosity (ν e ), set to 5 m 2 /s (Deltares, 2022a).
The bed shear stresses induced by flow velocity (τ c ) (Equation 7) and waves (τ w ) (Equation 8) are combined to a bed shear stress by flow and waves (τ cw ) with the equations of Fredsoe (1984) (Deltares, 2022a).
Herein, ρ [kg/m 3 ] is the density of water, g [m/s 2 ] is the gravitational acceleration, u [m/s] is the depth-averaged flow velocity, f w [ ] is the wave friction factor and u orb [m/s] is the near-bed wave orbital velocity.
The combined bed shear stress τ cw can dislodge mangrove seedlings when in excess of their tolerance to bed shear stress (τ cr,i ).In this study, seedling dislodgement was evaluated by comparing the tidal hydroperiod-averaged bed shear stress (i.e., the cumulative bed shear stress averaged over the hydroperiod during each tide; τ cw,TA ) to the tolerance of the seedlings at each preceding low tide.The hydroperiod was approximated with the maximum tidal (e) Assumed variability in probability of propagule abscission throughout the fruiting season.Model offshore (Station A, Figure 1a) boundary conditions for the reference season (Table 1) with (f) measured water level, (g) spectral wave height derived from wind velocity measurements and (h) suspended sediment concentrations derived from spectral wave height.
water depth (h M ) in the particular grid cell in relation to the tidal range, assuming a sinusoidal tidal wave.Only grid cells with a hydroperiod longer than 0.01 were considered (i.e., submergence longer than 1% of the tidal period).
In addition to the potential for dislodging mangrove seedlings, τ cw also influences the erosion rate (ER) and the deposition rate (DR) following the equations of Partheniades and Krone (Partheniades, 1965): Erosion occurs with an erosion parameter (M e , set to 0.00004 kg/m 2 s) when τ cw exceeds the critical bed shear stress for erosion (τ cr,e ).Deposition occurs at a rate of the sediment fall velocity (w s , set to 0.0004 m/s) and the depth-averaged sediment concentration (SSC).The critical bed shear stress for deposition (τ cr,d ) was set to 1000 N/m 2 .Erosion and deposition are computed when the local water depth exceeds 0.05 m.Suspended sediments are transported in the model domain following an advection-diffusion equation.Sediment is advected by the computed depth-averaged flow velocity and diffused by a sediment diffusion parameter modeled with an artificial eddy diffusivity (ε s , set to 5 m 2 /s).The abovementioned model parameters are quantified following the DFM model of the Firth of Gijsman et al. (2024).The cumulative bed level change due to erosion and deposition can erode or bury mangrove seedlings when the cumulative erosion (∑(ER-DR)) exceeds the seedling tolerance to erosion (E cr,i ) or the cumulative deposition (∑(DR-ER)) exceeds the shoot length (H i ).

Model Setup
For the DFM computations, a cross-shore oriented depth-averaged (2DH) spatial grid was setup for the Firth, spanning the subtidal and intertidal mudflats, the forest fringe and mature mangrove forest (Figure 1a).Waves were resolved on a 25m × 25m spatial grid with a cross-shore length of 7,300 m (i.e., the x-direction) from the NIWA climate station until the most seaward measurement station (station A), and an alongshore width of 600 m (i.e., the y-direction) to avoid boundary effects within the area of interest.Flow velocities and morphodynamics were resolved on a 20mx20m spatial grid centrally placed inside the wave domain, of the same length as the wave domain and with an alongshore width of 300 m.The spatial grid was refined to 10mx10m in the area of interest around the forest fringe.The wave spectra were resolved with 24 frequency bins between 0.05 and 1 Hz, while directions were resolved in bins of 2°.
Bathymetry of the model was interpolated from a cross-shore survey conducted in December 2022 with an RTK-GPS (Trimble GNSS).The vegetation cover for the model was converted from a field survey conducted in November 2016 (Montgomery et al., 2018).The tree density observations were linearly interpolated to determine the number of individual trees in each 10 m cross-shore section.Tree locations inside these cross-shore sections were randomly allocated, and tree heights were selected from a normal distribution with the mean and standard deviation of the observed tree heights (Montgomery et al., 2018).A fixed relation between tree height and tree diameter was used based on the conducted survey (H i = a D i b , with a = 14.72 and b = 0.58, p < 0.05, r 2 = 0.37) following Gijsman et al. (2024).Trees were allocated a C D value of 3 after Gijsman et al. (2024).Trees were treated individually for propagule abscission.For the interactions with hydrodynamics and morphodynamic changes, a grid-cell-dependent tree density was determined from the number of trees in each cell.Tree height and diameter and stem drag were grid-cell averaged.
The model was forced at the offshore model boundary with the water levels measured at the Tararu tidal station (Figures 1a and 4f).In addition, wave spectra and suspended sediment concentrations were imposed at the offshore model boundary (Figures 1a, 4g and 4h) based on a derived relation with the wind velocity and direction measured at the NIWA climate station for the monitoring period from May-June 2017 (Gijsman et al., 2024).The wave spectrum was a Pierson-Moskowitz spectrum with spectral wave height (H m0 ) and mean spectral wave period (T m01 ) and a limited directional spreading (power function cos 50 (θ θ p )).

10.1029/2024JF007664
Tides were solved on a Courant-number limited timestep of maximum 30 s.After every 15 min of flow and morphodynamic computations, a new wave field was generated with the updated water depths (i.e., 15 minute online flow-wave coupling).Propagule and seedling dynamics were resolved at each low water occurrence of the semi-diurnal tide at the offshore boundary (Figure 4f).For an overview of the most important DFM model parameters and additional details about the model setting, we refer to Gijsman et al. (2024).

Model Calibration
The hydrodynamic and morphodynamic processes are calibrated by Gijsman et al. (2024).The reference scenario was setup from 1st of November 2021 until 31st of October 2022 (Table 1), with a fruiting season from 1st of November 2021 until 1st of March 2022 (Figure 4e).Three parameters were calibrated: (a) the propagule abscission probability factor β, (b) the net tidal dispersal distance factor D i and (c) the initial bed shear tolerance of seedlings at establishment τ cr,i,0 .The reference scenario was repeated with increasing values of D i and τ cr,i,0 , and evaluated for the most seaward propagule and surviving-seedling location.One year old seedlings were observed up to 75 m seaward of station E in December 2022 (Figure 3) and it was assumed that propagules were dispersed further offshore.Lastly, β was increased until the maximum seedling density was approximately 4-5 stems/m 2 , as was approximated in the forest in December 2022.The model calibration yielded values of 0.0025, 4, and 1.4 N/ m 2 for β, D i, and τ cr,i,0 , respectively, resulting in a most seaward propagule location of 101 m w.r.t.station E, a most seaward seedling location of 95 m w.r.t.station E and a maximum seedling density of 4.2 stems/m 2 .

Model Scenarios and Output
The reference scenario with the calibrated model settings was repeated to explore the influences of (a) probabilistic propagule abscission and allocation of seedling physiology properties, (b) variations in environmental forcing conditions, and (c) mangrove forest presence on seedling dynamics.Firstly, the reference scenario was repeated 25 times each with a new random propagule abscission process as well as a new random allocation of seedling physiology (sim.#2-25 in Table 1) to quantify the influence of the probabilistic seedling properties and study the importance of ecological variability.With these 25 repetitions, a constant mean and standard deviation of surviving seedling density were obtained.Secondly, the reference scenario was repeated with different environmental forcing conditions (i.e., water levels, waves and suspended sediment concentrations) (sim.#26-36 in Note.The timeframe of the simulated environmental forcing conditions is presented.Simulations with variations in the probabilistic propagule/seedling input (1-25) were performed to study the influence of ecological variability.The influence of mangrove tide and wave attenuation was studied by enabling/disabling tide and wave attenuation (37).
Table 1) to study the influence of inter-annual variability in environmental forcing.For these runs, the boundary conditions from 1st of November 2021-31st of October 2022 were replaced by those from previous years, back to the year 2010-2011.Lastly, the reference scenario was repeated without the tide and wave attenuation caused by mangrove drag to study the influence of the mangrove trees present (sim.#37 in Table 1).
Model results regarding the seedling density and the physiology of the surviving and non-surviving seedlings were evaluated in 14 predefined sections, each spanning 40 m in cross-shore and 100 m in alongshore direction.As a result, the intertidal flat was represented by approximately two cross-shore sections ( 1,200 m < x < 1,120 m), the seedling establishment zone by three sections ( 1,120 m < x < 1,000 m), the dynamic forest fringe by three sections ( 1,000 m < x < 880 m) and the scrub forest by six sections ( 880 m < x < 640 m).The results were also averaged per forest zone.Additionally, seedling density was studied over time throughout and after the fruiting season.Lastly, for each computational grid cell, differences in annual-timescale forcing conditions (i.e., yearly median and 90th percentile hydroperiod and bed shear stress) were compared to the relative seedling establishment density (R E : the number of seedlings that established (S E ) relative to the cross-shore maximum (S E,max )) and seedling survival probability (P s : the fraction of established seedlings that survived (S S ) relative to S E ) of each year.Relations between R s , P s and the 90th percentile of the annual timescale hydroperiod (P) and tidal-hydroperiod averaged bed shear stress (τ cw,TA ) were explored with fit coefficients a, b, c and d following Equations 11 and 12:

R s = S E S E,max
= a • e b • P 90th percentile (11) 3. Results

Density and Physiology of Surviving and Non-Surviving Seedlings
The seedling density at the end of the reference simulation varied considerably across the elevation gradient of the transect through the Firth's mangroves (Figure 5).In the seedling establishment zone, seaward of the present forest, the simulated surviving seedling density was on average 0.025 seedlings/m 2 (i.e., 1 seedling per 40 m 2 ).Seaward of the seedling establishment zone on the intertidal flat, the average surviving seedling density was simulated to be smaller with 0.004 seedlings/m 2 (i.e., 1 seedling per 250 m 2 ).Landward of the seedling establishment zone in the dynamic forest fringe, the seedling density also decreased to 0.022 seedlings/m 2 (i.e., 1 seedling per 47 m 2 ) in the lower elevated area but increased rapidly toward the center (0.17 seedlings/m 2 or 1 seedling per 6 m 2 ) and higher elevated areas (0.87 seedlings/m 2 ).The dense scrub forest was characterized by the largest densities of surviving seedlings, exceeding 4.2 seedlings/m 2 .Further landward, the seedling density decreased toward the relict forest fringe to densities of approximately 1 seedling/m 2 (Figure 5).
The density of surviving seedlings was dependent on the probabilistic occurrence of propagule abscission and attribution of seedling physiology as well as on annual variability in environmental forcing (Figure 5).These dependencies varied across the elevation gradient.The repeated model simulations with different probabilistic propagule abscission and probabilistic allocation of seedling physiology (sim #1-25, Table 1) revealed that these probabilistic factors influenced surviving seedling density in the seedling establishment zone with a standard deviation (std.) of 0.0027 seedlings/m 2 (Figure 5b).In the dynamic forest fringe, the std.was larger with 0.0078 seedlings/m 2 , and in the scrub forest, the std.increased to 0.0288 seedlings/m 2 (Figure 5a).The model simulations with different environmental forcing (sim.#1 and 26-36, Table 1) produced std.values of 0.0039 seedlings/m 2 , 0.0580 seedlings/m 2 , and 0.0510 seedlings/m 2 in the seedling establishment zone, dynamic forest fringe and scrub forest, respectively.The relative variability of seedling survival due to probabilistic propagule input, calculated by the std.divided by the mean at each cross-shore section, was 0.11 in the seedling establishment zone, 0.022 in the dynamic fringe and 0.0082 in the scrub forest.The relative variability due to environmental forcing was 0.22, 0.17, and 0.014 for these respective forest zones.
The surviving seedling density decreased considerably when disregarding the damping of tidal flows and attenuation of wind waves by the present mangrove vegetation (Figure 5a).In the seedling establishment zone, minimal variability in surviving seedling density was simulated.In the dynamic forest fringe and the scrub forest, surviving seedling density dropped to an average of 0.31 seedlings/m 2 and 1.70 seedlings/m 2 , respectively.In the relict forest fringe, at approximately 400 m from the forest edge, the established seedling density remained similar regardless of the tide and wave attenuation by the mangrove trees.

Journal of Geophysical
The initial seedling root length was found to be a clear indicator of the survival probability of seedlings (Figure 6).The required inundation-free period, initial shoot height and growth acceleration/limitation factor appeared to have limited influence.The differences between the surviving and non-surviving seedlings varied across the elevation gradient (Figure 6).In the seedling establishment zone in particular, seedlings that survived had a mean initial root length (μ R i,0 ) of 0.50 cm that was substantially less negative than the seedlings that did not survive (μ R i,0 = 1.58 cm).Here, a negative initial root length indicated a delayed start of the root growth/sprouting processes and the root growth rate (G R ) was 0.33 cm/day.In the dynamic forest fringe, the difference in initial root length was smaller, with μ R i,0 = 0.85 cm and μ R i,0 = 1.64 cm for the surviving and non-surviving seedlings, respectively.In the lower dynamic forest fringe ( 1,000 m < x < 960 m), where the lowest surviving seedling density was simulated, seedlings also required a substantially faster growth rate with an average of μ M i = 1.08 to survive, in contrast to the remainder of the transect where μ M i of the surviving and non-surviving seedlings was almost the same.

Intra-Annual Seedling Dynamics
On an intra-annual timescale, the seedling density generally varied substantially during the fruiting season and in the month afterward (until approximately 150 days from the start of the fruiting season; Figure 7).The seedlings that survived until 2 months after the fruiting season ended (approximately 180 days from the start of the fruiting season) typically developed a large enough tolerance to bed shear stress, erosion and deposition to survive the rest of the year.The intra-annual development of seedling density was distinctly different between the seedling establishment zone, the dynamic forest fringe and the scrub forest (Figure 7).In the seedling establishment zone and in the scrub forest, the seedling density increased gradually throughout the fruiting season (Figures 7a-7c and 7g-7i).In contrast, in the dynamic forest fringe, particularly in the lower parts of this zone, the seedling density showed sharp fluctuations during the fruiting season (Figures 7d and 7e).

Variability in Seedling Survival Probability and Relation to Environmental Forcing
The variations in the surviving seedling density across the cross-shore elevation gradient and forest zones originate from cross-shore differences in the number of establishing seedlings (e.g., due to propagule supply) as well as differences in seedling survival probability (e.g., due to disturbances).Local propagule supply is imposed by the location of the present trees and tidal dispersal, and thus by elevation and hydroperiod.Seedling survival probability is influenced by the bed shear stresses occurring due to tides and waves, the local erosion and deposition, as well as the seedling physiology.With hydroperiods and bed shear stresses varying across the crossshore elevation gradients and forest zones, seedling establishment density and seedling survival probability accordingly show an inverted pattern (Figure 8).
While the interannual variability in seedling establishment density (S E ) was minimal (no visible spread in Figure 8a), the variability in seedling survival probability (P S ) was larger (Figure 8b).In the seaward intertidal flat and seedling establishment zones, the seedling survival probability varied with a standard deviation exceeding 10% between the years 2011-2022, while this was only 5% in the dynamic forest fringe (Figure 8b).The cross- shore and interannual variability in seedling establishment density can be explained by an exponential decay function in relation to extremer hydroperiods.The relation between the variation of the relative seedling density at the end of each simulated year (relative to the cross-shore maximum) and the yearly 90th percentile of the hydroperiod (i.e., P exceeded during 10% of the tides annually) has an r 2 -value of 0.82 and a root-mean-squarederror (rmse) of 0.15 (Equation 11; Figure 9a).This relation explains the variability in relative seedling density better than the relation with the yearly median hydroperiod (i.e., P exceeded during 50% of the tides; r 2 = 0.82; rmse = 0.28).The variability in seedling survival fraction can be attributed to variations in extremes of the tidalhydroperiod averaged bed shear stress.The relation between the seedling survival fraction and the 90th percentile of the yearly tidal-hydroperiod averaged bed shear stress (i.e., τ cw,TA exceeded during 10% of the tides) shows an r 2 -value of 0.90 (rmse = 0.12) (Equation 12; Figure 9b), while the annual median produced a poorer fit (r 2 = 0.82, rmse = 0.

Seedling Dynamics Across Elevation and Forest Zones in the Firth
Established seedling density and survival probabilities greatly varied between the lower-elevated intertidal flat and higher-elevated mature mangrove forest (Figure 5).In addition, the processes facilitating and/or limiting seedling dynamics varied across elevation and forest zones (Figure 7).Three characteristic cross-shore zones   The seedling survival fraction related to the 50th and 90th percentile tidal-hydroperiod averaged bed shear stress.Fits are based on alongshore-averaged results of simulations with mangroves for the period 2011-2022 (sim #1 and 26-36) and without mangroves for 2022 (sim #37; Table 1).r 2 values remain relatively high for all fits due to the larger concentration of data points with similar hydroperiod and bed shear stress.
Journal of Geophysical Research: Earth Surface 10.1029/2024JF007664 Seaward of the forest edge, on the intertidal flat and in the seedling establishment zone, both seedling establishment density (<5% of maximum) and survival probability (10%-50%) were low (Figure 8).In these zones, the seedling density increased gradually throughout the fruiting season at times when sufficiently strong propagules were dispersed from the forest, anchored their roots, and managed to survive in the occurring conditions (Figure 7) (Balke et al., 2014).Surviving seedlings typically showed a significantly shorter delay of the start of root growth (less negative initial root length) than non-surviving seedlings, showcasing the importance of rapid root anchoring and growth (Figure 6) (Balke et al., 2015;Van der Stocken et al., 2019).The survival probability increased further seaward at lower elevations, where water depths are larger during high tide conditions, the slope is milder (∼1:1000 slope) and wave breaking is less severe, all contributing to reduced bed shear stresses (Figure 8b).
In the dynamic forest fringe, the seedling establishment density increased (5%-50% of maximum) and their survival probability varied from low (<5%) to high (>80%) (Figures 5 and 8).Seedling density varied abruptly during the fruiting season in response to inundation-free periods during neap tide conditions, followed by seedling mortality due to wave disturbances during spring tide conditions (WoO-1 and WoO-2, respectively, in Balke et al., 2015) (Figure 7).The spring-neap tidal cycles caused large variations in propagule and seedling density, as observed by Swales et al. (2007) in the Firth.Reduced seedling survival probability is associated with the high bed shear stresses at the steepest part of the profile (∼1:100 slope), where a large portion of the waves break during high tide conditions (Figure 9).
The higher-elevated areas of the dynamic forest fringe and scrub forest are characterized by an increase in seedling establishment density and survival probability (Figure 5).Here, the supply of propagules is large and the area is sheltered from waves by the presence of trees in combination with a reduced hydroperiod, giving rise to a gradually increasing seedling density throughout the fruiting season (Figure 7).The attenuation of waves by mangrove trees contributes to the high seedling survival probability (Figure 8b).Further into the forest, high seedling survival probability occurred independently of the wave attenuation by mangrove trees (Figure 8b).The seedling establishment density reduced in that area due to the smaller tree density and hence reduced propagule supply (Figure 8a).
The distinct zonation in seedling dynamics in the Firth is dominated by water levels and waves and to a lesser extent by ecological variability of seedling properties.In other biophysical settings, ecological variability may increase in dominance and the transitions between forest zones may be more gradual for example, when more than one species is present, each having different traits and tolerances (Van Hespen et al., 2022;Yu et al., 2023).The presented results show that process-based biophysical models that incorporate both physical and ecological processes are essential for developing a mechanistic understanding of seedling dynamics across a wide range of biophysical settings.

Model Limitations and Application
The newly developed model represents the complexity of biophysical interactions between propagules/seedlings and hydrodynamics and morphological processes in the field, and also parameterizes seedling growth and its response to environmental conditions.Herein, a number of assumptions were made.Ecological factors that may influence propagule/seedling establishment and survival and growth such as the presence of macroalgae, availability of sunlight or nutrients, temperature, salinity, or consumption by crabs are not considered (Balke et al., 2015;Clarke & Kerrigan, 2002;Clarke & Myerscough, 1993;Krauss et al., 2008;Osunkoya & Crees, 1997), whilst these influences can also vary in time and space across elevation gradients and forest zones.
The presence of local-scale topographic features such as mud mounds and bed ridges that increase the spatial variability in soil strength and propagule/seedling dynamics are not considered.The presence of such features in the Firth could explain the large difference with reported seedling densities exceeding 36 seedlings/m 2 during the 2006 season (Swales et al., 2007).Moreover, the model is calibrated with observations of seedling densities approximately 1 year after the 2021 fruiting season.The reported seedling density in 2006 also reduced exponentially to below 1 seedling/m 2 after 180 days (Swales et al., 2007).
This study found that annual-timescale seedling survival probability is mostly hampered by bed shear stresses, being a proxy for sediment entrainment and thus local-scale morphodynamics.The derived relations for seedling establishment density and seedling survival probability can be implemented in longer-term morphodynamic mangrove development models that assume probabilistic establishment dynamics (Figure 9).In this way, longer-term models can incorporate annual seedling dynamics without the need to resolve hydrodynamics and morphological processes for the entire year.The resulting reduction in computational time can be invested in the consideration of other morphological and ecological processes relevant at these longer timescales, such as the influence of surface elevation and associated hydroperiod changes, competition for space and resources (Chen & Twilley, 1998;Clarke, 1995), or adaptive seedling growth in response to environmental forcing (Van Hespen et al., 2022).

Implications for Mangrove Forest Development, Recovery and Restoration
Variations in mangrove seedling dynamics across elevation gradients and forest zones are expected to influence mangrove forest recovery and restoration success.The identified cross-shore zones with distinctive seedling dynamics also comply with observations of forest zonation in the Firth with (a) low tree densities on the intertidal flat and in the seedling establishment zone, (b) higher tree densities with varying characteristics in the forest fringe and (c) high tree densities with rather uniform characteristics in the mature forest (Horstman et al., 2018;Lovelock et al., 2010;Montgomery et al., 2018;Swales et al., 2015Swales et al., , 2019)).This study shows that, in addition to cross-shore variations in and associated species competition (Van Loon et al., 2016;Watson, 1928), cross-shore and inter-annual variations in seedling dynamics also influence long-term forest characteristics and zonation.This becomes particularly clear in the Firth' mangroves that are exposed to waves, have relatively small cross-shore variations in hydroperiod and are composed of a single mangrove species.
While mangrove recovery from damage may be rapid in the mature forests, in the forest fringe, seedling recruitment and forest recovery may take years to decades.Rapid forest recovery in a storm-damaged fringe is dependent on the occurrence of inundation-free periods and limited disturbances (i.e., Windows of Opportunity) in subsequent spring-neap tidal cycles (Figure 7).In the Firth, past forest expansion was dominated by large establishment events initiated during calm El Niño years (Lovelock et al., 2010;Swales et al., 2015).More recently, however, the forest expansion and fringe development is more gradual and the inter-annual variability in established seedling densities is estimated at only ±44% in the seedling establishment zone and ±34% in the forest fringe.The recent gradual forest expansion rate of approximately 3.3 m/year (Figure 2) corresponds to the seaward expansion of the profile elevation (Gijsman et al., 2024), showing the positive feedback between profile elevation and forest development in the forest fringe.As a result, mangrove forest fringes may develop highly non-linearly, as has also been observed in mangroves in French Guyana (Proisy et al., 2009) and Vietnam (Nardin et al., 2016).
When mangrove restoration is attempted, the restoration design could be optimized in terms of location and timing.Restoration of dynamic intertidal flats may benefit from increasing the supply of propagules and temporarily supporting the root anchoring process, for instance by planting (Van Bijsterveldt et al., 2022).In addition, the placement of artificial permeable structures could enhance the seedling survival probability by reducing wind-wave driven bed shear stresses (Van Cuong et al., 2015;Winterwerp et al., 2020).In the Firth, the hydroperiod on the unvegetated mudflat is relatively low in comparison to other mangroves (Balke & Friess, 2016;Van Loon et al., 2016;Watson, 1928), but it is paramount that the physiological survival limits due to for example, increased periods of inundation are considered (Balke et al., 2015).The timing of restoration efforts that involve planting should consider the spring-neap tidal cycle and aim at optimizing the use of calm neap tide conditions, particularly in the forest fringe.Importantly, in the determination of the success of restoration efforts, the natural inter-annual variability in seedling survival probability due to variations in environmental forcing and ecological variability should be considered.Overall, it is paramount that (temporary) mangrove restoration measures should only be attempted in response to mangrove degradation and that they require careful consideration of potential negative side effects on the intertidal morphology and biodiversity.

Conclusions
This study identified cross-shore zones governed using biophysical processes that control mangrove seedling establishment and survival.On the intertidal flat, only seedlings that rapidly anchor were able to survive.In the mangrove forest fringe, the timing of seedling establishment with regard to the spring-neap tidal cycle was a key control.In the mature forest, the sheltered conditions provided an opportunity for high establishment and survival rates.The resulting differences in seedling establishment density and survival probability imply that the limitations to and rate of mangrove forest development and recovery vary between the daily inundated intertidal flat, daily to biweekly inundated mangrove fringe and less than biweekly inundated mature forest.The derived parameterizations of the relations between annually prevailing hydrodynamic conditions (i.e., hydroperiod and bed shear stresses), seedling establishment density and seedling survival probability provide a way to incorporate threshold-limited seedling dynamics into decadal-timescale schematized mangrove development models.The presented results and developed model can be used to inform mangrove management and restoration projects through the quantification and simulation of seedling establishment success in response to management strategies and/or environmental changes.

Figure 1 .
Figure 1.The mangrove forest in the southern Firth of Thames estuary, Aotearoa New Zealand.(a) Locations of field measurement stations where data were collected for model development and scenarios.The model domain is delineated by the black rectangle.(b) Cross-shore bathymetry of the upper intertidal flat with measurement stations, vegetation characteristics and vegetation zones (data from Horstman et al., 2018).Figure adapted from Gijsman et al. (2024).

Figure 2 .
Figure 2. The gradual expansion of the mangrove forest in the southern Firth of Thames between 2007 and 2022.The locations of the measurement stations E-H are indicated.The dashed lines indicate the borders between the distinctive forest zones (see Figure1b).The forest edge location is indicated.Aerial images obtained from Google Earth.

Figure 3 .
Figure 3. Field observation in the Firth of Thames (20 December 2022) showing seedlings established in the previous fruiting season from November 2021-February 2022.Aerial image from March 2023 (source: Google Earth).Measurement stations E-G and forest edge are indicated.Photo credit: R. Gijsman, 20 December 2022.

Figure 4 .
Figure 4. Growth curves for seedlings and their tolerances as well as model boundary conditions for the reference scenario (November 2021-October 2022).(a) Seedling shoot length growth curve.(b) Seedling root length growth curve with time.(c) Seedling tolerance to erosion growth curve with time.(d) Seedling tolerance to bed shear stress growth curve with time.Mean growth curves (solid lines) and standard deviations (dashed lines) are based on field measurements by Balke et al. (2015).(e) Assumed variability in probability of propagule abscission throughout the fruiting season.Model offshore (Station A, Figure1a) boundary conditions for the reference season (Table1) with (f) measured water level, (g) spectral wave height derived from wind velocity measurements and (h) suspended sediment concentrations derived from spectral wave height.

Figure 5 .
Figure 5. (a) Cross-shore density variations of the surviving and non-surviving seedlings with and without the sheltering provided by mangroves' wave attenuation.The variability in surviving seedling density caused by ecological variability (i.e., probabilistic propagule/seedling input, sim #1-25) and environmental forcing (sim #1 and 26-36) is included.(b) A close-up of panel (a) focusing on the seedling establishment zone.The cross-shore elevation profile, forest zones and characteristic water levels are presented in panels (a) and (b).

Figure 6 .
Figure 6.Distributions of the physiological characteristics of the surviving and non-surviving seedlings in all simulations combined, except for simulation without mangroves (sim.#37): (a) the required inundation-free period, (b) the initial shoot length, (c) the initial root length, and (d) the growth acceleration/limitation factor.
24).The calibration coefficients a, b, c and d in Equations 11 and 12 were fit with values of a = 3.395, b = 20.38,c = 1.045 and d = 0.408.

Figure 7 .
Figure 7.The temporal variability in seedling density in different cross-shore zones throughout and after the fruiting season for the simulations with mangroves for 2022 (sim.#1) and 2019 (sim.#28) and without mangroves for 2022 (sim.#37).(a-i) The seedling density in the different cross-shore zones, note the differences in scale between the horizontal axes (j) The tidal maximum water levels for 2019 and 2022.(k) An overview of the cross-shore coastal zones and characteristic water levels.

Figure 8 .
Figure 8. Alongshore-averaged seedling establishment density and seedling survival probability as well as 50th and 90th percentile of annual hydroperiods and bed shear stresses across forest zones.(a) Seedling establishment density of a simulation with mangroves for 2022 (sim.#1) and without mangroves (sim #37) as well as its variability due to environmental forcing (sim # 1, 26-36).Metrics of hydroperiods are presented on the right-hand axis.(b) Seedling survival fraction of simulations stated in (a).Metrics of tidal-hydroperiod averaged bed shear stresses are presented on the right axis.

Figure 9 .
Figure 9. Relative seedling establishment density and seedling survival probability at the end of each simulated year in relation to annual forcing conditions.(a) The relative seedling establishment density (relative to maximum seedling density along the cross-shore profile) related to the 50th and 90th percentile hydroperiod.(b)The seedling survival fraction related to the 50th and 90th percentile tidal-hydroperiod averaged bed shear stress.Fits are based on alongshore-averaged results of simulations with mangroves for the period 2011-2022 (sim #1 and 26-36) and without mangroves for 2022 (sim #37; Table1).r 2 values remain relatively high for all fits due to the larger concentration of data points with similar hydroperiod and bed shear stress.

Table 1
Model Simulations With Simulation Number and Simulation Name