Dune Dynamics Drive Discontinuous Barrier Retreat

Barrier islands and spits tend to migrate landward in response to sea‐level rise through the storm‐driven process of overwash, but overwash flux depends on the height of the frontal dunes. Here, we explore this fundamental linkage between dune dynamics and barrier migration using the new model Barrier3D. Our experiments demonstrate that discontinuous barrier retreat is a prevalent behavior that can arise directly from the bistability of foredune height, occurring most likely when the storm return period and characteristic time scale of dune growth are of similar magnitudes. Under conditions of greater storm intensity, discontinuous retreat becomes the dominant behavior of barriers that were previously stable. Alternatively, higher rates of sea‐level rise decrease the overall likelihood of discontinuous retreat in favor of continuous transgression. We find that internal dune dynamics, while previously neglected in exploratory barrier modeling, are an essential component of barrier evolution on time scales relevant to coastal management.


Introduction
Coastal barrier systems (barrier islands and spits), which account for ∼10% of the world's continental coastlines (Stutz & Pilkey, 2011), are of significant social, economic, and ecologic importance, yet the future of these landforms remains uncertain in the face of projected accelerated relative sea-level rise (RSLR; Sweet et al., 2017) and changes in tropical storm activity, including increasing intensity (e.g., Knutson et al., 2020). Barriers tend to migrate upward and landward in response to RSLR and storms, thereby maintaining subaerial exposure. Barrier migration occurs primarily through the process of overwash, where sediment eroded from the front of a barrier (shoreface and beach) is transported landward of the foredune crest to the barrier interior and back-barrier bay during storms (e.g., Dolan & Godfrey, 1973;Donnelly et al., 2006). Mounting evidence suggests that discontinuous (or "punctuated") barrier retreat, whereby a barrier oscillates between periods of migration and relative immobility over time scales of decades to centuries, is a prevalent behavior (Ashton & Lorenzo-Trueba, 2018;Ciarletta et al., 2019;Lorenzo-Trueba & Ashton, 2014;Mellett & Plater, 2018). This behavior contrasts with continuous transgression, where a barrier consistently retreats over time scales greater than the typical return period of overwash events. Sedimentological analyses of barriers and relict barrier deposits indicate that many barrier systems worldwide have experienced fluctuations in their rate of migration, which have been attributed to changes in external forcing such as the rate of RSLR (Cooper et al., 2016;Mellett et al., 2012); sediment supply (Forbes et al., 1991;Rampino & Sanders, 1982); shelf morphology and back-barrier accommodation (Nordfjord et al., 2009;Storms et al., 2008); antecedent Abstract Barrier islands and spits tend to migrate landward in response to sea-level rise through the storm-driven process of overwash, but overwash flux depends on the height of the frontal dunes. Here, we explore this fundamental linkage between dune dynamics and barrier migration using the new model Barrier3D. Our experiments demonstrate that discontinuous barrier retreat is a prevalent behavior that can arise directly from the bistability of foredune height, occurring most likely when the storm return period and characteristic time scale of dune growth are of similar magnitudes. Under conditions of greater storm intensity, discontinuous retreat becomes the dominant behavior of barriers that were previously stable. Alternatively, higher rates of sea-level rise decrease the overall likelihood of discontinuous retreat in favor of continuous transgression. We find that internal dune dynamics, while previously neglected in exploratory barrier modeling, are an essential component of barrier evolution on time scales relevant to coastal management.
Plain Language Summary Barrier islands and spits tend to move landward in response to sea-level rise and storms when sediment from the beach is washed into the barrier interior and beyond, but tall dunes at the front of a barrier can prevent this process from happening for all but the largest storms. Here, we explore the interactions between dunes, storms, and barrier migration with a new model (Barrier3D). Our model shows that barriers migrate when dunes are low but are stationary when dunes are tall. Over decades to centuries, repeated cycles of dune loss and regrowth lead to sporadic (or "discontinuous") barrier retreat. Additionally, we find that sporadic behavior may become less common in the future in response to rising sea levels, while increasing storm intensity will change which environments are most likely to experience discontinuous retreat. These findings emphasize the importance of dune dynamics in controlling barrier evolution. Barrier migration over decadal timescales -timescales relevant to coastal communities and managers -is commonly modeled as a constant background process, but our results suggest that is not always the case.

REEVES ET AL.
geology (Raff et al., 2018); and grain size (Forbes et al., 1995). Recently, experiments conducted using the exploratory model of Lorenzo-Trueba and Ashton (2014) suggest that decadal-to centurial-scale punctuated retreat can occur as the result of internal system dynamics, related specifically to lags in shoreface response to overwash and the assumption that overwash flux is related to a critical barrier width and height.
Here, we suggest that discontinuous retreat can also arise from dune dynamics, a suite of processes internal to barrier systems. Barrier dune systems tend to exist in one of two stable states as a result of feedbacks between dune recovery and storm erosion processes: a high-elevation, overwash-resistant state; or a low-elevation, overwash-vulnerable state (Durán Vinent & Moore, 2015;Durán Vinent et al., 2021;Goldstein & Moore, 2016). When biophysical dune recovery processes dominate, dunes tend to recover quickly if knocked down, remaining at or near their maximum height (Durán & Moore, 2013) and withstanding all but the largest of storms. When storm processes dominate, dunes cannot recover before the next storm and become trapped in a perpetual state of low elevation. However, when the time scales of dune recovery and storm occurrence are similar, foredune height becomes bistable (Durán Vinent & Moore, 2015;Durán Vinent et al., 2021;Goldstein & Moore, 2016). In this bistable regime, foredune systems respond nonlinearly to environmental drivers and are prone to abrupt shifts from one state to the other as the result of only small changes in environmental forcing. Under bistable conditions, the tendency for foredunes to exist in a state of either high or low elevation, and to alternate between the two states over time, has the potential to interact with barrier migration processes through alterations to overwash flux that occur with changes in dune height (Houser, Barrineau, et al., 2018;Sallenger, 2000). This fundamental linkage between dune dynamics and barrier evolution has yet to be explored using a long-term (decadal to millennial) modeling framework.
To investigate the role of dune dynamics in barrier migration, we introduce the new exploratory model Barrier3D, which builds upon earlier formulations for components of barrier evolution and incorporates previously neglected interactions between dunes and overwash. Previously, exploratory barrier models have coarsely resolved only large-scale, long-term changes in barrier position, height, and width and thus average over or neglect discrete events and associated temporal and alongshore variations in dune state (e.g., Lorenzo-Trueba & Ashton, 2014;Moore et al., 2010) in an effort to discover the processes most essential to barrier evolution. In contrast, Barrier3D tackles the scale separation between event-based and long-term models by explicitly, yet efficiently, simulating dune evolution, storm overwash, and a dynamically evolving shoreface in response to individual storm events and RSLR. Using this model, which we parameterize with data from Hog Island, VA, USA (a well-studied, undeveloped barrier island; Figure 1a), we explore barrier migration behaviors for a range of dune growth rates and storms frequencies. We then vary RSLR rates and storm intensities to explore how dune-storm interactions influence barrier migration under a range of anticipated future conditions.

Barrier3D
Barrier3D is an exploratory model (Murray, 2003(Murray, , 2013 that resolves cross-shore and alongshore topographic variations to simulate the morphological evolution of a barrier segment over time scales of years to centuries. The model operates with a 1-year time step over a 10-by-10 m grid. The barrier segment consists of one or more alongshore rows of dune cells at the front (ocean side) of the grid, backed by an interior domain with a predetermined constant alongshore length and dynamically changing cross-shore width ( Figure 1b). As explained below, dune and interior cells follow different sets of rules. Due to the complexities in modeling beach morphodynamics, and to focus on dune-storm interactions, Barrier3D does not simulate beach evolution and instead assumes invariant beach width and slope. Explanations of previously developed formulations and additional details not critical for understanding the fundamental workings of the model or its results are provided in the supporting information (Text S1) for completeness.

Storms and RSLR
RSLR is treated in the model using a Lagrangian frame of reference by reducing all elevations at the beginning of each modeled year relative to a sea level fixed at zero. Within each year, a probabilistically determined number of storms occur by randomly sampling from a normal distribution. We vary the mean of the distribution for our experiments and derive the standard deviation from a 35-year empirical storm record for Hog Island developed using tide gauge and wave hindcast data. Each individual storm in Barrier3D is described by three statistics: (a) duration; (b) R high , the highest elevation of the landward margin of runup; and (c) R low , the lowest runup elevation. The statistics for each storm are determined by randomly selecting from a list of 10,000 synthetic storms ( Figure 1c). We generate the synthetic storms using the multivariate sea-storm model from Wahl et al. (2016), which identifies interdependencies among the most relevant seastorm variables using a vine copula model with the 35-year empirical storm record as input (Text S1).

Dunes
The elevation of a dune cell (Z d ) is taken as the sum of the dune toe elevation (Z t ), which remains fixed over time relative to sea level, and the height of the dune crest (H d ) above its toe: d At the beginning of each modeled year, dune cells grow in a logistic manner following Houser et al. (2015) up to a maximum dune height (H d_max ; Durán & Moore, 2013): where r is the growth rate controlling the shape of the logistic curve. Within the model, the dune growth rate (r) varies randomly for each dune cell alongshore between a user-defined r min and r max , following observations from Houser et al. (2015). The characteristic dune growth rate (r ) is calculated as the mean of r min and r max . When varying r min and r max , we hold the range between the two values constant at 0.5 (Houser et al., 2015). Additionally, the static width of the dune domain (W d ) determines the number of rows (crossshore) that make up the active dune field, with d H  decaying exponentially for each row in the landward direction.
Where water elevations exceed dune elevations during storms, dune heights are reduced using a predictor of dune erosion as a function of total water level developed by Goldstein and Moore (2016). This predictor was created using data from Long et al. (2014) passed through a machine learning technique to yield a smooth nonlinear function relating the change in dune elevation ( d Z  ) to the elevation of water above the prestorm dune elevation (i.e., R high ), both normalized by the pre-storm dune elevation: REEVES ET AL.
Barrier3D does not simulate dune recovery during interstorm periods that fall within one model year, assuming that most dune growth likely occurs outside the peak storm season when the density of dune vegetation is higher and disturbance is relatively rare (e.g., Delgado-Fernandez & Davidson-Arnott, 2011;Montreuil et al., 2013;Ruz & Meur-Ferec, 2004;van Puijenbroek et al., 2017).
When the ocean shoreline erodes one cell length, the front row of the dune domain is removed and the first row of the barrier interior functionally becomes the back row of the active dune field (thereby maintaining W d ). When the ocean shoreline progrades one cell width, the most landward dune row functionally becomes part of the barrier interior (while retaining its elevation) and a new row is added to the front of the dune domain, with dune cell elevation set to Z t .

Overwash
To simulate overwash processes landward of the dune crest, Barrier3D uses a version of the cellular flow routing scheme developed by Paola (1994, 1997), in which water introduced at overtopped dune cells is carried landward row-by-row over the barrier interior, transporting sediment with it (Text S1). When all the water has been routed over the barrier interior, the model updates the barrier elevation according to the difference between the volume of sediment lost and gained from each cell. This process occurs iteratively for each hour of storm duration.
Sediment transported across the barrier interior to the back-barrier is distributed across the bay following exponential decay and contributes to landward extension of the back-barrier shoreline. The bottom of the back-barrier bay is flat (except where coarse, immobile overwash sediment has been deposited) and set to a constant depth (here 3 m), following the assumption that depositional and erosional processes for fine back-barrier sediment are able to maintain an equilibrium depth over the course of the simulation (Marani et al., 2007).

Shoreline Change
Shoreline change, applied uniformly across the barrier segment, follows the equations of Lorenzo-Trueba and Ashton (2014; Text S1). Shoreline erosion or accretion results from a combination of: (a) RSLR; (b) the cumulative volume of sediment removed from the upper shoreface by overwash and dune growth; and (c) net sediment exchange between the upper and lower shoreface. Shoreface sediment flux depends on the shoreface slope, which tends to dynamically adjust towards an equilibrium in response to perturbations (i.e., overwash and dune growth).

Experiments and Results
By examining barrier behavior across wide ranges of several key input parameter values, our simulations investigate barrier dynamics relevant to most barrier systems. To provide a common starting point for all simulations and to ground the model in reality, we parameterize the initial morphology and storm climatology for Hog Island within the Virginia Coast Reserve (VCR), VA, USA (Figure 1a). The initial morphology for the interior domain in our experiments consists of a digital elevation model (DEM) for a 500-m-long segment of Hog Island discretized into 10-by-10 m cells. H d_max and Z t are based on values extracted from a Hog Island DEM, and the synthetic storms are derived from wave hindcast data offshore Hog Island (USACE Wave Information Studies) and the nearest tide gauge in Wachapreague, VA. For a list of all input variables and the values used in these experiments, see Table S1.
Each simulation runs for 1,000 years or until the barrier drowns. Punctuated retreat is determined to occur if the simulation includes two or more alternating periods of both transgression and relative immobility; REEVES ET AL.

10.1029/2021GL092958
4 of 11 otherwise, the behavior is considered continuous. Periods of immobility (transgression) are defined as 30 years or more of shoreline change rates (SCR) under (over) 0.5 m/yr (see Text S1 for additional details).
Varying the duration and SCR thresholds used to define punctuated retreat produces predictable differences in our results but does not change the fundamental conclusions we draw from them ( Figures S1 and S2).

Dune-Driven Discontinuous Retreat
To examine how dune-storm interactions influence barrier migration, we run simulations across a range of storm frequencies, given by the average number of storms per year ( s N ), and characteristic dune growth rates (r ) (Figure 2). RSLR is held at 4 mm/yr. To account for storm stochasticity within each unique simulation, we run each of these parameter combinations 100 times and report the mean of these 100 simulations for our results.
We find that punctuated retreat is most probable diagonally across the center of the parameter space, where the storm return period and characteristic time scale of dune growth are of similar magnitudes (Figure 2a). Here, dunes are more likely to exhibit bistable behavior (Figure 2d), transitioning between stable states of tall, overwash-resistant dunes and low, overwash-vulnerable dunes. This dune behavior is consistent with previous models of dune-storm interactions (Durán Vinent & Moore, 2015;Durán Vinent et al., 2021;Goldstein & Moore, 2016). The bistability of dunes directly impacts overwash flux and, consequently, barrier migration, evidenced by instances of taller dunes occurring predominantly within periods of immobility and shorter dunes within periods of transgression (Figure 2d). Additionally, the average durations of immobile and transgressive periods are of similar magnitude where dunes are bistable (Figures 2b and 2c), suggesting that the barriers tend to spend relatively equal amounts of time stationary as actively retreating. During periods of immobility, when overwash fluxes into the barrier interior are limited, the barrier tends to lose elevation relative to sea level and become narrower from passive inundation of the back-barrier shoreline; the barrier tends to aggrade and widen during periods of transgression ( Figure S3).

REEVES ET AL.
10.1029/2021GL092958 5 of 11 When the storm return period and characteristic dune-growth time scale are of dissimilar magnitudes, punctuated retreat is less likely. The combination of fast dune growth and low storm frequency tends to result in the formation of tall dunes that limit overwash flux and reform quickly if knocked down (Figure 2e). This typically precludes periods of barrier transgression during the simulations in favor of continuous immobility and, ultimately, passive drowning via RSLR. As such, transgressive periods tend to be considerably shorter-lived than immobile periods (Figures 2b and 2c). Conversely, the combination of slow dune growth and high storm frequency tends to prevent dunes from ever recovering (Figure 2f), resulting in continuous barrier transgression. Therefore, the average duration of transgressive periods tends to be much longer than the average duration of immobile periods (Figures 2b and 2c). Punctuated retreat can occur within these corners of the parameter space (Figure 2a), sometimes as the result of anomalous periods of high or low storm activity or potentially time lags in the shoreface response to overwash ( Figure S4).

Sea-Level Rise
Next, we explore the impacts of RSLR on discontinuous retreat by running the same 25 combinations of s N and r presented in Section 3.1 (RSLR = 4 mm/yr) at two additional RSLR rates: 2 and 8 mm/yr ( Figure 3). Overall, the likelihood of discontinuous barrier migration decreases with higher rates of RSLR (Figures 3a-3c). Additionally, periods of transgression (Figures 3g-3i) tend to be much longer than periods of immobility (Figures 3d-3f) at higher rates of RSLR.
These results stem in part from the model implementation of the Bruun (1962) rule. Greater background shoreline retreat rates resulting from higher rates of RSLR make it less likely for periods of immobility to occur even in the absence of significant overwash. However, increased background shoreline retreat rates alone do not fully explain the trends we observe ( Figure S5). Higher rates of RSLR also typically lead to shorter dunes ( Figure S6), the result of increased landward migration of the dune crest in response to faster shoreline change rates (Keijsers et al., 2016;Rodriguez et al., 2018), and reduced interior elevations, which increases accommodation space for overwash. Both of these effects allow for greater overwash fluxes, and therefore greater barrier migration rates, tending to reduce the likelihood of immobile periods and therefore punctuated retreat in favor of continuous barrier transgression and, if sea level rises faster than overwash fluxes can maintain the barrier elevation, barrier drowning (though to maintain our focus on migration behaviors influenced by dune dynamics, we do not explore high enough RSLR rates in our simulations to produce this particular drowning behavior).

Storm Intensity
To explore the impacts of changing storm intensity, we run the same 25 combinations of s N and r presented in Section 3.1 with both higher and lower storm intensities (Figure 4), with storm intensity represented by the TWL. Instead of drawing randomly from the 10,000 simulated storms as before, we bin them by R high , fit a beta distribution to the histogram, and draw from the distribution shifted +0.15 (−0.15) m such that bins containing storms with higher (lower) R high values are preferentially chosen; storms are then selected randomly from the specified bin. We refer to the three levels of storm intensity as "low," "observed" (i.e., the same simulations presented in Figure 2), and "high." Increasing storm intensity shifts the maximum potential for discontinuous migration towards barriers under conditions of lower storm frequency (Figures 4a-4c). Dune bistability (and therefore discontinuous retreat) becomes more likely at lower storm frequencies because a greater proportion of storms overtop the dunes when storms are more intense. Barriers that typically experience continuous immobility as a result of tall, resilient dunes at low storm intensity transition to discontinuous retreat at high storm intensity as dunes become bistable; barriers that typically exhibit punctuated retreat at low storm intensity transition to continuous transgression at high storm intensity as low-elevation dunes become the single attracting state of the dune system. The same shifts occur for the average durations of immobile and transgressive periods (Figures 4d-4i), as more barriers experience longer transgressive periods and shorter immobile periods with higher storm intensities.
REEVES ET AL.

Discussion and Conclusions
The goal of this work is to explore and explain large-scale barrier behavior and the key feedbacks that give rise to it, rather than to numerically predict the evolution of any particular location or environment. The standard exploratory modeling approach we adopt here -reducing the number of processes involved and treating them in considerably simplified ways to enhance the clarity of potential insights -supports this goal, but bears several limitations worth considering. For instance, we do not model a beach domain, and instead assume that beach characteristics remain invariant over time and space. In reality, beach characteristics (e.g., slope and width) can vary greatly over hourly to decadal time scales (e.g., O'Dea et al., 2019) and 10 2 -10 3 m length scales (Vos et al., 2020) in ways that can affect dune growth (Durán & Moore, 2013), dune erosion (Beuzen et al., 2019), and overwash (Donnelly et al., 2006). Additionally, the model does not include tidal inlet processes, which can account for a significant proportion of overall transgressive flux in a barrier system (e.g., Leatherman, 1979;Nienhuis & Lorenzo-Trueba, 2019). The current model formulation also lacks some of the sediment pathways that can contribute to deflation of the barrier (i.e., narrowing and flattening; Passeri et al., 2020), such as gradients in alongshore transport, breaching, or storm-driven seaward transport, which would all tend to increase the vulnerability of the barrier to drowning in our model. Nevertheless, although simplifications may limit the quantitative precision of our results, our model experiments provide valuable qualitative insights into barrier behavior and how it may change in the future. We REEVES ET AL.
10.1029/2021GL092958 7 of 11 find that punctuated retreat is a common behavioral regime, consistent with recent modeling (Ashton & Lorenzo-Trueba, 2018;Ciarletta et al., 2019) and observational (e.g., Mellett & Plater, 2018;Raff et al., 2018) studies. In contrast to previous work, however, the discontinuous retreat in our experiments is driven principally by dune dynamics, though this does not invalidate other explanations for discontinuous behavior or suggest that multiple mechanisms cannot operate in conjunction. In fact, intrinsic couplings between dune dynamics and other proposed forcing mechanisms can alter the potential for dune state change. For example, RSLR can increase the vulnerability of dunes to erosion and therefore overwash (a dynamic we capture in the model; see also Rodriguez et al., 2018), and changes in sediment supplied to beaches (e.g., from the shoreface) can limit or bolster dune recovery (Houser, Barrineau, et al., 2018). In the discontinuous retreat of our model, dune dynamics lead to interior elevation loss (relative to MSL) during periods of immobility. This contrasts with the shoreface-driven discontinuous retreat from the model of Lorenzo-Trueba and Ashton (2014), in which periods of immobility are characterized by barrier aggradation (i.e., the barrier maintains elevation relative to sea level; see Text S2 for a comprehensive comparison of discontinuous retreat in the two models).
In response to projected accelerated RSLR, our experiments suggest that dune-driven discontinuous retreat will become less common overall in favor of continuous transgression as immobile periods become shorter and more infrequent. Meanwhile, increasing storm intensity will shift the range of storm frequencies under which punctuated retreat is most likely to occur, becoming the dominant behavior of barriers that were REEVES ET AL.
10.1029/2021GL092958 8 of 11 previously stable. Barrier migration behaviors could vary latitudinally, as a result of latitudinal gradients in storm characteristics and, conceivably, dune growth rates. Therefore, projected climate-driven shifts in maximum storm intensity (e.g., Knutson et al., 2020) and the distribution of dune grass species (e.g., Goldstein et al., 2018) may result in latitudinal shifts in barrier behavior; these are potentially fruitful topics for future research.
Dune-driven punctuated retreat is supported by empirical evidence from Durán Vinent and Moore, (2015) relating lower dune heights to faster shoreline change rates from the VCR. Houser, Wernette, et al. (2018) observe a similar association between dune height and storm-scale (<10 years) shoreline change from barrier islands in Texas and Florida, but raise important questions about length scales given that a single barrier can exhibit km-scale reaches of both high and low dunes (Zinnert et al., 2016). On barriers exhibiting kmscale alongshore variability in dune state, alongshore transport would tend to redistribute sediment from sections with tall dunes towards overwashing areas without dunes (Houser, Wernette, et al., 2018), tempering local effects of dune-driven punctuated retreat. The extent of punctuated behavior thus depends on the km-scale alongshore continuity of dune state, as well as the strength of alongshore diffusion of shoreline change rates.
These results demonstrate the importance of internal controls and nonlinear state changes in barrier evolution. Barrier rollover is commonly perceived (and modeled) as a constant background process, but over time scales of even a few decades to several centuries, our results suggest that is not always the case. For modeling barriers on time scales relevant to coastal communities and managers, we find that explicitly including dune dynamics -erosion, recovery, and interaction with storm overwash processes -is essential in determining both the rate and style of barrier retreat. Further, dune-driven nonlinearity in barrier retreat raises questions about the reliability of linear regression for calculating long-term historical shoreline change, which might mischaracterize the historical behavior or perceived stability of a barrier. Without contextualization for dune-storm dynamics, such historical observations could fail to provide insights necessary for accurately forecasting coastal change. Future work is needed to address the potential for cascading state changes across the entire barrier-marsh-bay system triggered by dune-driven discontinuous barrier retreat.

Data Availability Statement
Barrier3D is available for download from the online Community Surface Dynamics Modeling System model repository at https://csdms.colorado.edu/wiki/Model:Barrier3D. Data for model experiments are in repository at https://doi.org/10.6073/pasta/d2b133e539f87f38f241b97e64ddaa06.