Permanent loss of barrier island resilience due to a critical transition in dune ecosystems

The study investigates what controls the transition of barrier islands from a vegetated, dune-backed 'high' state to a dune-less 'barren' state, and the implications for coastal protection, ecosystem integrity, and island evolution under sea level rise (SLR). Barrier islands buffer inland coasts from storms, with dunes playing a central role by elevating island topography and reducing overwash frequency. Without dunes, barriers are prone to frequent flooding, accelerated transgression, potential exposure of carbon-rich backbarrier deposits, and reduced protection for infrastructure and ecosystems. Barrier elevation thus serves as a proxy for barrier state and resilience, with dune formation and recovery after storms being key indicators. Prior modeling work ranges from simple mass-conservation approaches that omit dune dynamics to process-based models that capture storm-scale impacts but often simplify long-term dune behavior. Because overwash occurrence is primarily set by barrier elevation, a consistent description of barrier migration and resilience requires resolving dune dynamics. A recent stochastic point model captured the competition between wind-driven dune growth and stochastic water-driven overtopping/erosion, defining barrier elevation states via a probability density function (PDF) parameterized by overtopping event frequency and size. Observations from the Virginia Barrier Islands (VBI) show three qualitative states—'barren,' 'mixed,' and 'high'—raising the central research question: what controls transitions among these states, particularly from 'high' to 'barren'?

Previous approaches to barrier island dynamics include: (1) Simple, planform-averaged models using sediment mass conservation that do not resolve dune dynamics; (2) Process-based models (e.g., XBeach, CShore) that simulate storm-scale impacts and overwash but are computationally intensive and not typically used to predict long-term stochastic elevation states; (3) Large-scale morphodynamic and probabilistic models focused on long-term response to sediment supply, storms, and SLR that often simplify dune dynamics. A recent stochastic point model (Vinent et al., 2021) provides a physics-based, analytically tractable framework for dune growth versus stochastic overtopping, yielding elevation PDFs governed by overtopping frequency and size. Empirical and numerical studies show dunes grow toward a maximum height H controlled by vegetation and shoreline-related stresses, often following an exponential approach with characteristic formation time Ta. Prior work has also characterized extreme total water level statistics as marked Poisson processes with exponential marks, supporting the stochastic representation of high-water events (HWEs). Together, these studies motivate extending the point model alongshore to compare with spatially extensive elevation datasets and to identify control parameters governing transitions between barrier states.

- Study area and data: Virginia Barrier Islands (VBI). USGS LiDAR-derived DEMs for 2014, 2016, 2017 were referenced to a representative beach elevation (1.5 m NAVD88 shift based on HWE/runup analysis). Homogeneous island segments were selected; Metompkin was split into North and South due to differing dune presence. Cross-shore transects spaced at 10 m alongshore were used to extract barrier elevation h(y,t) as the maximum elevation per transect, capturing berms, overwash fans, foredunes, and secondary dunes. - Stochastic point model: Barrier elevation h is bounded between base elevation h0 (post-overwash elevation, i.e., beach berm/washover surface) and maximum dune height H (set by vegetation and shoreline stress). Dune growth follows an exponential saturation to H with characteristic formation time Ta and maximum growth rate Ga; water-driven erosion during HWEs instantaneously lowers dunes to h0 when overtopped. The stochastic elevation change is dh = (H − h) dt − Δh(h,t), with Δh > 0 during overtopping events. HWEs are modeled as a marked Poisson process with exponentially distributed sizes; frequency of events overtopping h0 is λ0 = λ1 exp(−h0/S), where λ1 and S are the frequency and average size of HWEs overtopping the reference beach elevation. Control parameters at a point are A = λ0 Ta and S̃ = S/(H − h0), with ξ = (h − h0)/(H − h0) the normalized elevation. The steady-state point PDF fξ(ξ; A,S̃) is known analytically; it is bimodal when 0 < ξmin < 1. - Alongshore extension: Assume alongshore randomness in H(y) and h0(y) is approximately normal (Gaussian) and quasi-stationary over timescales ~Ta. The alongshore elevation PDF fh(h) is the integral of the point PDF conditioned on local H and h0 over their alongshore distributions NH and N0. For tractability, global control parameters A and S̃ are evaluated at mean H and h0, reducing dependency to NH and N0 and two control parameters A, S̃. - Parameter estimation: From elevation change rates G(y) = [h(y,2017) − h(y,2014)]/3 yr versus initial elevation h(y,2014), identify equilibria where G ≈ 0 and G′ < 0. Estimate: (i) Base elevation distribution N0 by fitting low-change elevations with a Gaussian (mean h0, stdev σh0). (ii) Maximum dune height distribution NH using the high-elevation mode when present (H, σH); if no clear dune mode, approximate from Hog Island or equilibrium from G(h). (iii) Maximum dune growth rate Ga as the average of the top ten (excluding the maximum) growth rates for initial elevations > ~0.4 m. (iv) HWE parameters for the region: average size S ≈ 0.3 m and frequency λ1 ≈ 18 events/yr overtopping the reference beach elevation. Compute Ta = (H − h0)/Ga, A = λ0 Ta with λ0 = λ1 exp(−h0/S), and S̃ = S/(H − h0). - Model evaluation: Numerically integrate the alongshore PDF formulation to predict fh(h) and compare to empirical PDFs for 2014, 2016, 2017 across islands (North Metompkin, Smith, Cedar, Parramore, South Metompkin, Hog). Also compute the rescaled mean dune recovery time Tr (mean excursion time below ξmin divided by Ta) as a diagnostic of barrier state and resilience. - SLR scenario and critical transition: Define a phase curve by varying mean h0(t) under relative SLR (RRSLR = 10 mm/yr, 2020–2050) while keeping other parameters as in Hog Island to simulate approach to a tipping point. Track fh(h) and Tr along the curve to identify critical slowing down and loss of the dune-mode basin of attraction.

- Two stable elevation modes and states: The stochastic model predicts and observations confirm two primary stable modes at the mean base elevation h0 (barren, dune-less) and at the mean maximum dune height H (high, mature dunes). A 'mixed' regime occurs between them. The alongshore PDFs for VBI islands exhibit these regimes and are reproduced by the model. - Tipping point near beach berm elevations: A critical transition occurs when barrier elevations are near the beach berm/base elevation (~0.4–0.5 m). Around this threshold, erosion dominates dune formation, leading to a stable barren state. VBI data suggest the transition to a barren state occurs at mean h0 ≈ 0.4 m, with the beach berm a strong attractor in growth-rate analyses G(h). - Control by rescaled HWE frequency A and sensitivity to h0: The barrier state is governed by control parameters A = λ0 Ta and S̃ = S/(H − h0). Transitions to barren are driven primarily by increases in A, particularly sensitive to mean h0 because λ0 ∝ exp(−h0/S). A roughly 30-fold increase in A across VBI islands maps their progression from high to barren states. - Recovery time as resilience metric and critical slowing: The rescaled dune recovery time Tr quantifies resilience. For Cedar Island, model Tr ≈ 30; with Ta ≈ 6 years, this implies ~180 years for dune recovery, effectively rendering the barrier dune-less at steady state. As h0 declines, Tr increases sharply (critical slowing) approaching the tipping point; below the threshold, the dune-mode basin disappears and the system relaxes to the barren equilibrium. - Sea-level rise effect and timescales: Because h0 is referenced to water level, SLR lowers h0, increasing A, decreasing resilience, and driving an effectively irreversible transition to the barren state. Under RRSLR = 10 mm/yr, a 0.25 m decrease in h0 (comparable to the difference between South Metompkin and Parramore) would occur in ~25 years, potentially shifting barrier state within decades (likely an underestimate given expected increases in wave-driven flooding). - Robustness: The transition depends weakly on H near the mixed-to-barren boundary (product A·S̃ is H-independent to first order). Predictions are relatively insensitive to details of extreme event probability and dune erosion representation, as they focus on post-storm recovery of proto-dunes rather than mature dune erosion. - Model–data agreement: Predicted steady-state PDFs reproduce the modes and main shapes of empirical alongshore elevation distributions across three years for multiple islands, supporting the model’s relevance despite simplifications.

The research clarifies the physical controls on barrier island state transitions by linking dune-building dynamics and stochastic water-driven erosion through two measurable control parameters. It shows that mean base elevation h0, via its exponential control on overtopping frequency, dominantly determines whether dunes can persist or whether barriers become dune-less. The model explains observed elevation PDFs and distinct island states in VBI, and provides a resilience metric (Tr) that captures critical slowing as systems approach a tipping point. This directly addresses the hypothesis that barrier elevation states emerge from non-linear stochastic dynamics with multiple equilibria. The results indicate that SLR, by lowering h0, can rapidly and irreversibly shift barriers to barren states, with cascading implications: greater flooding exposure, accelerated island transgression, potential changes in coastal carbon budgets, and reduced ecosystem and infrastructure protection. The framework offers a pathway to incorporate dune recovery dynamics and overwash-driven sediment fluxes into larger-scale coastal models that traditionally omit dune processes, improving predictions of barrier migration and hazard exposure.

This study extends a stochastic dune–overwash point model alongshore and tests it against LiDAR-derived elevation data from the Virginia Barrier Islands to reveal a critical transition from dune-backed ('high') to dune-less ('barren') states. The model identifies two stable elevation modes (near h0 and H), quantifies a tipping point near beach berm elevations (~0.4–0.5 m), and introduces a recovery-time-based resilience metric that exhibits critical slowing near the transition. Results indicate that the transition is primarily governed by the rescaled frequency of HWEs, exponentially sensitive to h0, and that SLR can drive a rapid, effectively irreversible shift to the barren state within decades. These insights support more physically grounded parameterizations of dune dynamics and overwash in large-scale coastal models and inform coastal management by identifying early-warning indicators of resilience loss. Future work should: (i) relax simplifying assumptions (e.g., all-or-nothing erosion, spatial constancy of Ga), (ii) couple shoreline change and vegetation dynamics explicitly, (iii) test the framework on wider or more complex multi-dune systems and armored coasts, and (iv) integrate human-coastal feedbacks and sediment management interventions into the stochastic resilience framework.

- Model assumptions: Dune erosion is simplified as complete lowering to h0 upon overtopping; real systems may experience partial erosion. The model assumes rapid vegetation colonization relative to dune formation. Ga is treated as temporally constant over ~annual scales. - Applicability: Not suitable where dunes are not eroded to h0 during HWEs, where dune growth rates vary strongly alongshore, for wide multi-ridge dune systems, or where complex beach armoring alters overwash/erosion processes. - Simplifications: Shoreline change and explicit vegetation dynamics are neglected in the core predictions, though they influence H and dune morphology. The probabilistic description of extremes is simplified. - Data constraints: Elevation change rates are derived from snapshots three years apart, averaging over several realizations of the stochastic process; some parameters (H, σH) for islands lacking clear dunes are approximated from representative sites (e.g., Hog), introducing uncertainty.