Simulating barrier island response to sea level rise with the barrier island and inlet environment (BRIE) model v1.0

Introduction

Barrier islands are long, narrow, sandy stretches of land that occupy a significant fraction of modern coastlines around the world. Barriers are often densely populated, support diverse ecological communities, and protect bays and wetlands that provide a range of ecosystem services (Barbier et al., 2011;McLachlan, 1983). Despite their importance, there is a critical gap in our ability to predict how barriers will respond to coastal change generally and sea level rise (SLR) specifically. A necessary condition for barrier islands to migrate landwards and keep up with SLR is sufficient sediment transport from the barrier front to the top and back via overwash fan deposition and flood-tidal delta formation (Armon and McCann, 1979;Inman and Dolan, 1989;Kraft, 1971;Lorenzo-Trueba and Ashton, 2014;Mallinson et al., 2010;Moore et al., 2010). There are few constraints, however, on the potential magnitudes of these landward sediment fluxes and how these fluxes vary as a function of the coastal setting, wave climate, or SLR. Recent models (e.g., Lorenzo-Trueba and Ashton, 2014) have suggested formulations for overwash fluxes, but the potential role of tidal fluxes, their feedbacks with overwash deposition, and the resulting ability of barriers to keep pace with SLR remain unclear.

Here we present the BarrieR Inlet Environment (BRIE) model to address this fundamental knowledge gap. Transgression in the model is driven by two main processes: overwash sedimentation and flood-tidal delta deposition (Leatherman, 1979;Pierce, 1969Pierce, , 1970)). To date, models that have aimed to assess barrier island change over geological timescales typically account for only storm overwash, which is more suitable for a cross-sectional framework. Tidal inlets, however, have been suggested to contribute a large fraction of the transgressive sediment movement in a number of field studies (Pierce, 1969(Pierce, , 1970)). The BRIE model extends the formulations of Lorenzo-Trueba and Ashton (2014) (LTA14) in the alongshore direction and incorporates tidal inlet morphodynamics through Delft3D-derived parameter-Published by Copernicus Publications on behalf of the European Geosciences Union.

izations (Nienhuis and Ashton, 2016, NA16). The purpose of the model is twofold: (i) to better understand long-term barrier island morphodynamics, including the effects of, for example, sea level rise, human development (jetties, beach nourishment), or storm pattern changes, and (ii) to improve paleoenvironment reconstructions.

Section 2 of this paper provides a background on barrier island environments and recent model developments. In Sect. 3, we discuss model formulations, including overwash fluxes, alongshore sediment transport, and tidal inlet morphodynamics. Section 4 includes a model run that demonstrates the capabilities of the BRIE framework, including inlet dynamics alongshore and the generation of alongshore stratigraphic profiles. Section 5 explores model sensitivity to grid and time resolution, as well as a comparison to other barrier island models. We conclude with a few exploratory results and a discussion of potential model applications.

Barrier islands and SLR

Barrier islands are narrow strips of land formed by waves through a variety of (hypothesized) mechanisms (e.g., Gilbert, 1885;McGee, 1890;Penland et al., 1985) associated with relatively slow SLR rates and primarily passive margins (FitzGerald et al., 2008;Stutz andPilkey, 2011, McBride et al., 2013). The emergence of many barrier islands can be traced back to about 6000 years before present, when Holocene SLR slowed down (McBride et al., 2013).

However, the relationship between barrier islands and SLR is complex. Under no SLR, barrier islands are generally not observed as their associated back-barrier environments would fill completely (e.g., Beets and van der Spek, 2000). In contrast, under moderate SLR rates marshes and tidal flats generally occupy back-barrier environments. In this case, to maintain their elevation with respect to sea level, barriers migrate towards land as storm overwash and flood-tidal flows deposit sediment. Under higher SLR rates, however, it is more difficult for barriers to maintain their subaerial portion above sea level. Consequently, when onshore-directed sediment fluxes are insufficient, barrier islands drown in place and are left offshore (Rodriguez et al., 2001;Mellet and Plater, 2018). Additionally, when onshore-directed sediment flux events are very intense and frequent, barrier islands are unable to maintain their geometry as they rapidly migrate towards land, which also results in drowning (Lorenzo-Trueba and Ashton, 2014). This potentially delicate balance between SLR and barrier response, together with the current projections of future acceleration in SLR, highlight the need to better constrain onshore-directed sediment fluxes in different barriers island systems (Carruthers et al., 2013;Lazarus, 2016;e.g., McCall et al., 2010;Rogers et al., 2015).

Barrier overwash

One way for sediment to be transported across the barrier is through storm overwash. Differences in water-level setup between the ocean and the lagoon during a storm can force the flow of water and sediment through and above the subaerial portion of the barrier. Most frequently this flow is directed landward, resulting in the transport of sediment from the ocean to the bayside where it deposits as the flow spreads laterally into the lagoon (Carruthers et al., 2013;Donnelly et al., 2006). Although this process is complex and highly intermittent, individual storm events integrated over time result in a net landward sediment flux, which allows barriers to keep pace with SLR over geological timescales (Leatherman, 1983). Despite its importance in terms of future barrier island morphodynamic response and vulnerability to flooding (Miselis and Lorenzo-Trueba, 2017), this long-term landward sediment flux is generally poorly constrained, and its relationship with modern overwash fluxes is not straightforward (Carruthers et al., 2013;Donnelly et al., 2006;Lazarus, 2016;Rogers et al., 2015). This lack of constraints on longterm overwash fluxes has resulted in a suite of barrier island models that do not compute overwash processes as a function of single storm events. Instead, such models parameterize overwash volume fluxes as a function of barrier geometry and observations of barrier island migration. For example, Leatherman (1979) observed that narrow barrier islands tend to be more susceptible to overwash events than wide barrier islands. They defined a "critical barrier width" below which overwash is frequent and the barrier migrates rapidly and above which overwash and barrier migration tend to be slow. Based on these findings, overwash is often parameterized by assuming the volume is inversely proportional to island width (e.g., Jiménez and Sánchez-Arcilla, 2004;Lorenzo-Trueba and Ashton, 2014) and additionally adjusted based on local factors such as land use (e.g., Rogers et al., 2015).

Tidal inlets

Aside from storm overwash, tidal inlets have also been found to be a major contributor to barrier transgression (Inman and Dolan, 1989;Moslow and Heron, 1978;Pierce, 1969). Tidal inlets derive their transgressive potential through the deposition of flood-tidal deltas. The volume of flood-tidal delta deposits correlates with the size of the associated inlet (Powell et al., 2006). Simple equilibrium models (e.g., Stive et al., 1998) suggest that initially flood-tidal deltas grow fast but that their growth slows down as they approach an equilibrium volume and the bay fills up near the inlet. Inlet migration can therefore add to transgressive transport by exposing new bay to flood-tidal delta deposition (Nienhuis and Ashton, 2016). For these two reasons, it has been hypothesized that shortlived and rapidly migrating inlets are most efficient for barrier transgression (Pierce, 1970). However, if the migration rate and life span of tidal inlets correlate with sediment im-port, their potential for transgression should then depend on factors such as basin size, ocean waves, and tidal conditions.

Along extended barrier coastlines, barrier morphodynamics are complicated by the existence of multiple tidal inlets. Tidal inlets interact through their control on water surface elevation in the tidal basin. This interaction can cause inlets to close or change size (van de Kreeke et al., 2008;Roos et al., 2013). Observations of tidal inlet spacing (Davis and Hayes, 1984), corroborated by a recent modeling study (Roos et al., 2013), found that increasing tidal range and basin size can allow inlets to exist closer together. For a barrier coast, a greater number of inlets likely enhances their contribution to barrier transgression.

Previous numerical modeling efforts

The joint long-term effect of storm overwash and tidal inlets on barrier island evolution remains difficult to quantify. On the one hand, engineering models typically assess barrier island changes over annual to decadal timescales, which includes overwash fluxes and tidal inlet formation during storm events. For example, models such as XBeach (Mc-Call et al., 2010;Roelvink et al., 2009) resolve wave dynamics coupled with sediment transport during storm events and are able to capture barrier morphological changes, including breaching. On decadal timescales, models like Delft3D (Deltares, 2014) have been applied to study inlets, but these typically do not include the effect of storms or SLR (e.g., Tung et al., 2009;NA16). On longer timescales, models no longer use laboratory-validated sediment transport relationships but rather use various degrees of conceptual relationships between barrier geometry and barrier island movement (Cowell et al. 1995;Storms et al., 2002;Stolper et al., 2005;Masetti et al., 2008;Wolinsky and Murray, 2009;LTA14). Some of these models are morphokinematic, i.e., based upon the conservation of mass and the maintenance of barrier geometry (Cowell et al., 1995, Wolinsky andMurray, 2009;Stolper, 2005). The models developed by Storms et al. (2002), Masetti et al. (2008), and LTA14 are morphodynamic as they account for sediment fluxes along the shoreface and across the barrier island. LTA14 represents a significant simplification compared to other morphodynamic models, making it suitable for model extensions and model coupling (such as the one presented here).

Coming from a different angle, the ASMITA (Aggregated Scale Morphological Interaction between Tidal inlets and the Adjacent coast) model couples coasts to their back-barrier environment via sediment exchanges determined by the deviation of a morphological element (ebb delta, tidal flat, etc.) from an assumed equilibrium volume. ASMITA has been developed in part to understand the effects of SLR on inlets and their back-barrier environments (van Goor et al., 2003;Stive et al., 1998;Townend et al., 2016). Inlets cannot close or migrate, and the model does not account for overwash processes (Stive et al., 1998). In ASMITA, as well as other back-barrier models (van Maanen et al., 2013;Mariotti and Canestrelli, 2017), the maximum potential sediment import through tidal inlets exerts a first-order control on the ability of back-barrier environments to sustain themselves during SLR.

Here, we describe a new model (BRIE) that accounts for tidal inlet dynamics, including opening, closing, and lateral migration, combined with barrier overwash processes as described by LTA14. Within the realm of coastal geomorphological models, BRIE can be considered a large-scale coastal behavioral (LSCB) model (de Vriend et al., 1993). It seeks to represent only the main governing mechanisms of the coast at appropriate timescales, without fully resolving the mechanics of fluid and sediment transport. There is a rich body of literature concerning LSCB models, ranging from rocky coasts (Walkden and Hall, 2011), barrier islands (Stolper et al., 2005;LTA14), tidal basins (Townend et al., 2016), tidal inlets (Kraus, 2000), and sandy coastlines (Ashton et al., 2001), to aggregates of LSCB models that couple these elements (Ashton et al., 2013;Payo et al., 2017). To our knowledge, BRIE would be the first to explicitly couple barrier islands and tidal inlet morphodynamics. Despite its simplicity, BRIE provides a novel approach to study the evolution of barrier islands under decadal to millennial timescales. Moreover, it allows us to explore complex barrier dynamics across a wide range of parameter values.

Model

We developed the BRIE modeling framework to study barrier island response to SLR. The model incorporates longshore interactions by linking the cross-shore barrier island model presented by LTA14 in a series of dynamic cross-shore profiles. We apply the storm overwash and shoreface response functions independently in each cell (Fig. 1). Feedbacks between overwash dynamics alongshore arise through the coupling with alongshore sediment transport, which can adjust the shoreline location and influence the shoreface slope and barrier overwash. The model also accounts for the formation, closing, and migration of tidal inlets following the parameterizations from NA16 (Fig. 1).

To our knowledge, this is the first morphodynamic model for long-term (decadal to millennial timescales) barrier island evolution that accounts for both tidal and overwash sediment fluxes. The model is written in MATLAB, and a typical runtime for a 100 km long barrier island stretch over a 10 000-year simulation is ∼ 1 min. Its simplicity and computational speed enable us to explore model behavior for a wide range of parameter values.

General description and model setup

After initializing the environment (typically ∼ 100 km long barrier island with periodic boundaries) and determining wave climate and shoreface parameters (Table 1), we run each time step in a for loop. In each iteration, we first raise sea level (Fig. 2). SLR affects subaerial barrier volume and shoreface slope, which in turn drives overwash and shoreface fluxes (Sect. 3.2).

Next, we determine if new inlets should be formed (Sect. 3.3.1), in which case we analyze their hydrodynamics and calculate their equilibrium dimensions. For each inlet, we distribute sediments into the flood-tidal delta, the barrier island, and the shoreface (Sect. 3.3.3). Flood-tidal delta deposition changes the back-barrier location (Fig. 1). After each time step, we add the different sources and sinks to the coastal zone, including diffusive wave-driven alongshore sediment transport, and implicitly determine a new shoreline, back-barrier, and shoreface toe position. See Table 1 for an overview of all model parameters and units.

Cross-shore barrier model

At a minimum, barrier islands can be described in the crossshore dimension as composites of three regions: the active shoreface on the ocean side, the subaerial portion of the barrier island, and the back-barrier lagoon on the terrestrial side, where infrequent overwash processes determine the volume of onshore-directed sediment fluxes (LTA14). Note that the overwash model is applied independently for every alongshore cell j , from 1 to n y (Fig. 1), but we leave out these indices for clarity. Assuming an idealized geometry, the crossshore evolution of the barrier system can be fully determined with the rates of migration of the shoreface toe (x t ),

Table 1. Model variables and their dimensions. Shortened references are as follows: LTA14 (Lorenzo-Trueba and Ashton, 2014), B80 (Bowen, 1980), N15 (Nienhuis et al., 2015), R13 (Roos et al., 2013), SH85 (Sallenger and Holman, 1985), SZ09 (de Swart and Zimmerman, 2009), AM06 (Ashton and Murray, 2006), NA16 (Nienhuis and Ashton, 2016). Eq. ( 6) m 3 m -1 yr -1 overwash flux deposited in the back barrier (LTA14) Q ow,h

Eq. ( 5) m 3 m -1 yr -1 overwash flux deposited on top of the existing barrier (LTA14)

Eq. ( 8)

Eq. ( 7) m 3 m -1 barrier width deficit (LTA14) z 0 H s /γ m minimum integration depth for shoreface flux k sf Eq. (12 m 3 m -1 yr -1 shoreface response rate (LTA14) w s Eq. ( 16) m s -1 settling velocity Q sf Eq. ( 9)

Eq. ( 24) -offshore wave direction θ arctan(dx s /dy) -shoreline angle 23)

cross-shore width of the lagoon L lagoon m alongshore length of the lagoon draining to a particular tidal inlet A inlet Eq. ( 21) m 2 inlet cross-sectional area u

Eq. ( 22) m s -1 inlet flow velocity

barrier cross-sectional area updrift of an inlet α Eq. ( 26) -fraction of Q s,in deposited as new barrier (NA16) β

Eq. ( 28) -fraction of Q s,in that bypasses the inlet (NA16) δ

Eq. ( 30) -fraction of Q s,in that deposits as flood-tidal delta (NA16) α r

Eq. ( 27) -eroded barrier deposited as new barrier, fraction of Q s,in (NA16) β r

Eq. ( 29) -eroded barrier deposited in the littoral zone, fraction of Q s,in (NA16) δ r

Eq. ( 31) -eroded barrier deposited in flood-tidal delta, fraction of Q s,in (NA16)

Model output

Eq. ( 47) m 3 m -1 yr -1 transgressive flux due to tidal inlets F Eq. ( 46) -fraction of the transgressive flux from inlets the shoreline (x s ),

the back-barrier shoreline (x b ),

and the barrier height above sea level (H ),

where Q sf is the sediment flux at the shoreface, ż is the SLR rate, Q ow,h is the top-barrier overwash component, and Q ow,b is the back-barrier overwash component. Other variables and parameters are defined in Table 1. Note that Eqs.

(1)-( 4) follow the barrier island model of LTA14 except for the (1 -f ) factor in Eq. ( 2) that accounts for fine-grained sediment in the back barrier. We discuss this modification in Sect. 3.2.2. We compute overwash flux using a simple formulation that assumes the existence of a critical barrier width (W b,crit ) and a critical barrier height (H crit ) beyond which there is no overwash to the back and the top of the barrier, respectively. When the barrier width (W b ) and height (H ) are below their critical values, the overwash rates Q ow,h and Q ow,b scale with their associated deficit volumes, V d,h and V d,b , resulting in an overwash flux heightening the barrier,

and an overwash flux widening the barrier,

We define the volume deficits with respect to an equilibrium defined by the critical barrier width and height (LTA14). In this way, we can compute V d,b and V b,h as follows:

The shoreface flux (Q sf ) is controlled by the shoreface response rate (k sf ) and the deviation of the shoreface slope from its equilibrium slope,

Modifications to the LTA14 barrier model

All the above formulations are identical to LTA14 except for Eq. ( 2), which we adjust to account for fine sediments in the back barrier. LTA14 assumes a back-barrier depth geometrically determined as zx b • s background (Fig. 1), where s background is the basement slope. This depth assumes the absence of back-barrier sediment deposition (i.e., f = 0, where f is the fine sediment fraction) and therefore represents the upper-limit depth. The BRIE model accounts for fine sediment deposition by selecting a back-barrier depth D lagoon (see Eq. 3) that is within the range 0≤ D lagoon ≤ zx b • s background . We then compute the fine sediment thickness in the back barrier (Fig. 1) as

In turn, we can geometrically define f as follows:

As barriers migrate towards land, fine sediments are absorbed in the bayside and exported at the shoreface on the ocean side; this is a dynamic that can play a significant role in the total barrier sediment volume changes (Brenner et al., 2015). BRIE accounts for fine sediment export at the shoreface by assuming that the fine sediment fraction f given by Eq. ( 11) is representative of the entire cross section of the barrier and that the sediment exchange between the upper and lower shoreface Q sf is not affected by the presence of fine sediments. In this way, Eq. ( 2) accounts for the fact that the fine sediment fraction f of the overwash sediment volume extracted from the shoreface does not contribute to the total volume of the barrier. In other words, only a fraction (1 -f ) of the shoreface volume eroded (i.e., x s,ow,sf (2H + D T )) deposits on top and/or back of the barrier (Fig. 3).

Parameter estimation for the LTA14 model

The cross-shore barrier model is a function of several parameters, including the shoreface depth D T , the equilibrium shoreface slope s sf,eq , and the shoreface response rate k sf . These three parameters, although generally poorly constrained, can be estimated as a function of wave and sediment characteristics (e.g., characteristic sediment grain size D 50 , significant wave height H s ). This allows us to investigate how storm overwash, alongshore transport, and inlet dynamics covary for a particular environment.

The shoreface response rate can be viewed as the integrated cross-shore sediment transport flux between a depth z 0 below wave breaking and the shoreface depth D T (LTA14; Ortiz and Ashton, 2016). Here we integrate the shoreface flux k sf (converted from meters per second into units of meters per year),

where φ = 16e s C s ρ/(15π(ρ s -ρ)g) and H (z) is the local wave height at depth z. We solve this integral assuming H (z) is a shallow-water wave that can be estimated by the offshore wave climate and a shoaling coefficient,

We derive a simple analytical expression of the integrated shoreface response rate,

where we estimate z 0 as the breaking wave depth H s /γ , and γ is 0.4 (Sallenger and Holman, 1985). We determine the shoreface depth D T (m) using an empirical relationship based on the wave characteristics (Hallermeier, 1981),

We estimate the shoreface equilibrium slope s sf,eq as the slope at the depth of closure (Lorenzo-Trueba and Ashton, 2014),

where the settling velocity is calculated based on the empirical formulation developed by Ferguson and Church (2004),

Inlet formation

We allow the model to form new tidal inlets every T storm years at the location of minimum barrier volume A barrier , where T storm can be considered a storm return time. An inlet can only form at a distance of at least L min away from current inlets, where L min is a minimum inlet spacing (Roos et al., 2013). Although L min is likely dependent on a wide range of factors, we are not aware of field constraints on its value and therefore choose a constant L min . the flood-tidal delta by increasing the back-barrier location,

with the implicit assumption that the flood-tidal delta top is approximately at sea level. Although inlets cannot open closer than L min away from existing inlets, differences in inlet migration rates can cause inlets to exist closer to each other (and merge). Additionally, inlets can also form when a section of the barrier drowns (negative barrier cross-sectional volume, i.e., A barrier < 0), regardless of the distance to other existing inlets.

Inlet hydrodynamics

At every time step, we compute the distance among all inlets. Assuming the lagoon water drains to the nearest inlets, we determine the lagoon area per tidal inlet (the potential for tidal prism) by multiplying the water surface area (i.e., W lagoon • L lagoon ) with a predefined fraction occupied by marshes, f marsh . We compute inlet characteristics such as cross-sectional area and flow velocities based on de Swart and Zimmerman ( 2009), who in turn followed ideas established by Escoffier (1940) (Fig. 4). We solve the inlet area-velocity relationship (Escoffier, 1940;de Swart and Zimmerman, 2009) analytically for u = u e , meaning that inlets adjust to maintain an equilibrium tidal velocity amplitude whereby sediments will be neither deposited nor eroded. In this situation, the nondimensionalized equilibrium inlet cross-sectional area is given by

where F 0 is

In this formulation, A H = ωW b / √ ga 0 is a resonance nondimensional cross-sectional area, u e = u e / √ ga 0 is a nondimensional equilibrium velocity, and γ is the ratio of the potential tidal prism and the inlet friction,

where the drag coefficient

lagoon (de Swart and Zimmerman, 2009).

Based on A inlet we determine the dimensional inlet crosssectional area (m 2 ),

and the corresponding water velocity through the inlet,

The inlet area function (Eq. 21) evaluates to the largest cross-sectional area for which u = u e (1 m s -1 in all simulations) if an equilibrium inlet area exists (Fig. 4). The inlet velocity evaluates to u < u e if an equilibrium does not exist (the friction through the inlet exceeds the potential tidal prism), at which point the inlet will close. Inlets adjust instantaneously to changes. Waves do not influence the size of the inlet, but alongshore sediment transport is assumed to be present to maintain an inlet to its equilibrium size.

Alongshore sediment transport into inlets

We calculate alongshore sediment transport into inlets (Q s,in , converted from cubic meters per second to cubic meters per year) based on the CERC formula recast into deepwater wave properties (AM06),

where k is a constant that is ∼ 0.06 m 3/5 s -6/5 (Nienhuis et al., 2015), and φ 0 is the wave direction. Shoreline orientation θ is defined by x s / y. We determine the wave direction at every time step from a cumulative distribution function defined by the wave asymmetry a and wave highness h (Ashton and Murray, 2006),

where x is uniformly distributed between 0 and 1.

Note that although we estimate sediment transport into inlets based on this method, we do not calculate shoreline change based on this particular wave angle at every time step. Instead, for model stability and efficiency, we calculate shoreline change using an implicit time step nonlinear diffusion equation, with inlets, storm overwash, and cross-shore shoreface transport acting as sediment sources or sinks (see Sect. 3.5.2).

Inlet morphodynamics

After we have determined inlet cross-sectional area and wave-driven transport into the inlet, we distribute sediments between the updrift and downdrift portions of the inlet and the flood-tidal delta, following the parameterizations of NA16 (Fig. 1). Inlets can migrate and erode into a barrier and also deposit a barrier. Inlets also form flood-tidal deltas. Ebb-tidal deltas are absent from this formulation because they do not present a sink from the littoral zone. Ebb-tidal deltas, however, implicitly determine the rate of inlet migration and the size of the flood-tidal delta through their effect on waves and currents (NA16).

Inlet migration and flood-tidal delta deposition rates are dependent on the alongshore sediment transport into the inlet Q s,in and the sediment distribution fractions α, β, δ, β r , α r , and δ r (Fig. 1; NA16). These fractions are determined by Delft3D model experiments and parameterized as

where A b,updrift and A b,downdrift are W b • (D inlet + H ), the barrier cross-sectional area updrift and downdrift of the inlet, respectively. We estimate the sediment distribution fractions based on the inlet momentum balance I , which is the ratio of the tidal

For model stability, we depart from the original formulation of NA16 on two occasions.

i. Equation ( 31) is a departure of the original formulation (NA16). The new function forces both inlet flanks to migrate at the same rate, making inlet width purely a function of inlet hydrodynamics.

ii. We impose a maximum flood-tidal delta volume following Powell et al. (2006) such that V fld,max = 1 × 10 4 u e A eq /2ω 0 0.37 (m 3 ). If this maximum is reached, we limit I to 0.1 to ensure more efficient bypassing and inlet migration. In the original parameterization of NA16, flood-tidal delta deposition (δ) is not a function of flood-tidal delta size. In BRIE this would create unrealistically large flood-tidal deltas.

Based on the sediment distribution, the inlet can deposit sediment into the flood-tidal delta. Assuming that the floodtidal delta is at sea level, we can describe its rate of growth (Fig. 1c) as follows:

change the sediment budget in the littoral zone,

and migrate alongshore in the direction of the littoral drift,

Changes to the back-barrier and shoreline locations are estimated at every time step. Inlet migration, however, per time step t (∼ 0.05 years) is typically much less than the alongshore discretization y (∼ 100 m). We therefore track inlet migration by assigning a "fraction migrated" to one grid cell in each inlet. The inlet moves along the barrier if that fraction exceeds one or drops below zero. New barrier island is constructed at sea level, H inlet = -H . A second complication is that inlets are also typically (but not necessarily) wider than the alongshore discretization y. Inlets are therefore allowed to exist on multiple alongshore cells j , dependent on the inlet width, n inlet,i = W inlet / y, where n inlet is the number of alongshore cells taken up by inlet i (Fig. 1).

Shoreline change

After we have determined the various sources and sinks of sediment to the nearshore environment, we distribute sediment alongshore between the different cells based on alongshore sediment transport. We use an implicit Crank-Nicolson scheme (Crank and Nicolson, 1947) to solve for shoreline change, governed by the following nonlinear diffusion equation:

which includes the effect of wave refraction and shoaling and is therefore suitable to apply based on offshore (deepwater) wave conditions (AM06). We have added a source-sink term x s,j,ext (m) to account for cross-shore sediment movement. D j is a nonlinear term and accounts for the fact that diffusivity depends on the wave approach angle (AM06),

where k is ∼ 0.06 m 3/5 s -6/5 (Nienhuis et al., 2015) and θ (j ) = (x s,j +1 -x s,j )/ y. is the angle dependence of the diffusivity (AM06) that we compute as follows:

which we convolve with the normalized angular distribution of wave energy E(φ),

to generate a long-term, wave-climate-averaged shoreline diffusivity for every alongshore location j . We rewrite the shoreline diffusion Eq. (28) into

where n and j denote the specific time and space locations. We solve this equation by inverting this nearly tri-diagonal matrix:

where β j = D n j t/2/ y 2 . Because we use D at n instead of n + 1 this is simply a linear diffusion equation. Indices in the lower right and upper left corner indicate periodic boundary conditions. The source x s,j,ext (m) can be described by

representing offshore and onshore sediment fluxes that can erode and accrete the shoreline and flood-tidal delta deposition that acts as a littoral sink. The shoreline model is unconditionally stable and secondorder accurate in space and time. We discretize the coastline into cells with width y (typically 100 m). We use a time step t (typically 0.05 years) to ensure smooth inlet migration and reasonably accurate shoreline change. However, we note that the overwash and inlet elements of this model are not solved by Eq. ( 41) and are therefore not necessarily second-order accurate nor unconditionally stable. Section 5 presents the grid and time resolution tests.

Stratigraphy module

Aside from the usual output, such as transgressive fluxes, inlet morphodynamics, and barrier island change, the model can also compute the synthetic stratigraphy of a barrier at a certain location x strat for all grid cells j (Fig. 6). When x b exceeds x strat , the model saves the location j , lagoon depth D lagoon , the sediment deposit thickness (i.e., D lagoon -z), and the responsible process, either flood-tidal delta deposition or storm overwash. While x s < x strat < x b , we record the height of the barrier H as dune construction or erosion bounded vertically by z and H . If an inlet is present, it erodes the deposit up to a depth d inlet . Inlet migration forms sedimentary facies between d inlet and z. If an inlet is closed, it forms inlet fill facies. These barrier island facies allow us to compare model output to geological reconstructions of barrier islands (e.g., Mallinson et al., 2010). 4 Example model runs

Model without inlets

We first investigated a simulation without inlets, focusing on the effect of alongshore transport gradients and barrier overwash on barrier evolution. As we might expect, in the case of no inlets and uniform initial conditions alongshore, the barrier retreats uniformly and alongshore sediment fluxes do not affect barrier response. We also performed a model experiment with an initially variable barrier width driven by spatial changes in the bay shoreline location (Fig. 5). In this scenario, the initially narrower barrier stretches overwash more than the wider stretches and therefore transgress faster. As shoreline curvatures increases, the magnitude of the alongshore sediment fluxes directed to the narrow stretches also increases, which reduces the width of the initially wider stretches. Interestingly, we find that time lags in shoreline interconnectivity can cause the initially rapidly transgressing stretch to stay in place and eventually become landward of other portions of the coast, a phenomenon also reported by Ashton and Lorenzo-Trueba (2018). Eventually, after a few oscillations that can last for hundreds of years, the barrier approaches a spatially uniform migration rate (Fig. 5).

Model with inlets

Including tidal inlets, we see a richer set of model dynamics.

We investigated barrier change, including inlets, for an SLR rate of 2 mm yr -1 , a wave height of 1 m, and a tidal range of 1 m. After an initial spin-up phase associated with large overwash fluxes, barrier island response stays highly dynamic and does not converge to an equilibrium response, despite the imposition of constant boundary conditions (Fig. 6). Inlets open, close, interact, and migrate preferentially with the direction of the littoral drift. Inlet migration rates vary gradually, and inlet sediment distribution is initially dominated by alongshore sediment bypassing and gradually becomes more flood-tidal delta dominated (Fig. 6e). The inlet transgressive sediment flux is highest when the flood-tidal delta deposition and alongshore sediment bypassing are roughly equal (Fig. 6f). Barrier stratigraphy at that time shows that inlet migration facies make up most of the barrier, even though not all of the transgression is due to the inlet (Fig. 6b). 5 Model tests

Conservation of mass

To investigate model mass conservation, we summed the volume of the barrier and offshore deposits (Fig. 7). Comparison to an identical model without inlets shows that slight losses and gains can be attributed to inlet morphodynamics, likely inlet migration and closure (Fig. 7b). For example, we do not track the sediment lost or gained as inlets change their crosssectional area from an initial breach width. We also assume that increases in the back-barrier location can be considered small enough so that there is one depth D lagoon , whereas in reality these deposits exist on a surface with slope s background .

Regardless of these assumptions, model volume (offshore deposits and the barrier island itself) does not depend on the time step t and the grid length y; these values only deviate a few percent around their mean, with no obvious trend in time (Fig. 7b).

Comparison to the 1-D model

For model verification, we compared model results to the original cross-shore model of barrier change that only includes overwash (LTA14) (Fig. 8). Our model without inlets produces the same dynamics as the original crossshore model (LTA14), resulting in the same overwash flux (Fig. 8a). Comparing the cross-shore model to the BRIE model with inlets (forced nonuniformity) we see some clear differences. Even though the average shoreline location along the 100 km barrier follows roughly the same trajectory (Fig. 8c) and therefore has a similar transgression (erosion) rate (Fig. 8b), the individual locations vary significantly. The straight barrier reproduced by the BRIE model without inlets is now variable alongshore. Transgression rates vary from -2 m yr -1 (progradation) to +10 m yr -1 (erosion). Even though the overall trajectory is a result of the sea level history and the passive inundation of the main (nonbarrier) coast (Wolinsky and Murray, 2009) (Fig. 8c), significant deviations from this trend appear and are reflected in the overwash rates (Fig. 8a). In particular, the inlet transgressive sediment flux rates are variable.

Model evaluation

It is challenging to evaluate long-term barrier island models against natural examples. Given their erosional nature, longterm records and barrier dynamics are scarce (Mellett and Plater, 2018). Thus, instead of a direct comparison to natural examples, we evaluate our model by exploring the sensitivity of the model output to a variety of boundary conditions (Fig. 11). We find that, even though individual simulations show great variability over time (Fig. 6), longer-timescale dynamics of barrier islands present physically meaningful relationships with model boundary conditions. For example, wave height tends to increase the effect that inlets have on barrier transgression, likely by making inlets more wave dominated and by increasing their migration rates. Inlets are The greater effect of inlets for an intermediate depth could be because flood-tidal delta growth is enhanced, thereby restricting tidal flow and forcing the opening of inlets elsewhere (Fig. 11).

Discussion and conclusion

We have built a 2-D barrier island model (i.e., the BRIE model) to simulate barrier island response to SLR that couples alongshore sediment transport processes, storm overwash, and tidal inlet morphodynamics. The mathematics of the approach are verified by comparing model predictions without inlets against the LTA14 cross-shore model. We also show that sediment volume is conserved with sufficient accuracy under a wide range of scenarios. Model results demonstrate that feedbacks between shoreface dynamics, barrier overwash, and alongshore transport processes can result in a complex history of interconnected behavior between the shoreline and barrier location. Moreover, we find that the relative importance of tidal inlets and storm overwash in transporting sediments onshore during barrier landward migration can significantly vary as a function of a wide range of factors, including sea level rise rate, wave climate, barrier and inlet geometries, and antecedent topography. Overall, model results highlight the importance of the interplay between crossshore and alongshore processes, particularly tidal processes, in understanding future and past barrier response to sea level rise.

The BRIE modeling framework does not aim to reproduce the evolution of any particular field location. Instead, we focus on exploring the relative role of tidal and overwash fluxes in the response of barriers to SLR, which requires omitting processes that could also play a significant role. For instance, the BRIE model does not account for human activities and coastal protection strategies along the coast (e.g., sea walls, groins, beach nourishment), which are known to affect coastal response at different spatial and temporal timescales (Jin et al., 2013;Murray et al., 2013). Rather than accounting for marsh-lagoon dynamics in the back-barrier environment, which can potentially influence the rate of barrier landward migration under sea level rise (FitzGerald et al., 2008;Lorenzo-Trueba and Mariotti, 2017), we define a fine sediment thickness based on the lagoon depth and the basement slope. We also ignore the stochastic nature of storms, as well as the potential dynamic influence of shoreface lithology. Given its simplicity, however, the BRIE modeling framework can be extended to account for additional processes that might affect barrier evolution, including the ones mentioned above.