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

Barrier islands are low-lying coastal landforms vulnerable to inundation and erosion by sea level rise. Despite their socioeconomic and ecological importance, their future morphodynamic response to sea level rise or other hazards is poorly understood. To tackle this knowledge gap, we outline and describe the BarrieR Inlet Environment (BRIE) model that can simulate long-term barrier morphodynamics. In addition to existing overwash and shoreface formulations, BRIE accounts for alongshore sediment transport, inlet dynamics, and flood–tidal delta deposition along barrier islands. Inlets within BRIE can open, close, migrate, merge with other inlets, and build flood–tidal delta deposits. Long-term simulations reveal complex emergent behavior of tidal inlets resulting from interactions with sea level rise and overwash. BRIE also includes a stratigraphic module, which demonstrates that barrier dynamics under constant sea level rise rates can result in stratigraphic profiles composed of inlet fill, flood–tidal delta, and overwash deposits. In general, the BRIE model represents a process-based exploratory view of barrier island morphodynamics that can be used to investigate long-term risks of flooding and erosion in barrier environments. For example, BRIE can simulate barrier island drowning in cases in which the imposed sea level rise rate is faster than the morphodynamic response of the barrier island.

1 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, 1969, 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, 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 parameterizations (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.

2 Background

2.1 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 and Pilkey, 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).

2.2 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 long-term 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).

2.3 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 short-lived 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 import, 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.

2.4 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 (McCall 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 and Murray, 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.

3 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).

Figure 1 Schematized model domain in (a) the plan view highlighting the three moving boundaries (the shoreface toe, shoreline, and the back-barrier (lagoon) location) and the sediment fluxes that determine their coupling. Vectors indicate the direction of potential changes, with the dot symbolizing a movement up. (b) Implementation of fine sediment dynamics into the barrier overwash model. (c) A close-up of (a), showing the littoral sediment fractionation within an inlet and its translation to barrier change. (d) A close-up of (b), showing the barrier volume deficit approach. See Table 1 for model variable names and units.

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.

Figure 2 Model structure showing the time loop in which the model updates SLR and calculates resulting barrier island change.

3.1 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).

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).

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.

3.2 Cross-shore morphodynamics

3.2.1 Cross-shore barrier model

At a minimum, barrier islands can be described in the cross-shore 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 cross-shore evolution of the barrier system can be fully determined with the rates of migration of the shoreface toe (x_t),

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{t}}}{\mathrm{\Delta }t}}=\mathrm{4}{Q}_{\mathrm{sf}}{\displaystyle \frac{H+D_{\mathrm{T}}}{D_{\mathrm{T}}\left(\mathrm{2}H+D_{\mathrm{T}}\right)}}+{\displaystyle \frac{\mathrm{2}\dot{z}}{s_{\mathrm{sf}}}},\end{array}$$

the shoreline (x_s),

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{s},\mathrm{ow},\mathrm{sf}}}{\mathrm{\Delta }t}}={\displaystyle \frac{\mathrm{2}\left(Q_{\mathrm{ow},\mathrm{b}}+Q_{\mathrm{ow},\mathrm{h}}\right)}{\left(\mathrm{2}H+D_{\mathrm{T}}\right)\left(\mathrm{1}-f\right)}}-\mathrm{4}{Q}_{\mathrm{sf}}{\displaystyle \frac{H+D_{\mathrm{T}}}{{\left(\mathrm{2}H+D_{\mathrm{T}}\right)}^{\mathrm{2}}}},\end{array}$$

the back-barrier shoreline (x_b),

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{b},\mathrm{ow}}}{\mathrm{\Delta }t}}={\displaystyle \frac{Q_{\mathrm{ow},\mathrm{b}}}{H+D_{\mathrm{lagoon}}}},\end{array}$$

and the barrier height above sea level (H),

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{\Delta }{H}_{\mathrm{ow}}}{\mathrm{\Delta }t}}={\displaystyle \frac{Q_{\mathrm{ow},\mathrm{h}}}{W_{\mathrm{b}}}}-\dot{z},\end{array}$$

where Q_sf is the sediment flux at the shoreface, $\dot{z}$ 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,

$$\begin{array}{}\text{(5)}& Q_{\mathrm{ow},\mathrm{h}}=Q_{\mathrm{ow},\max }{\displaystyle \frac{V_{\mathrm{d},\mathrm{h}}}{max(V_{\mathrm{d},\mathrm{b}}+V_{\mathrm{d},\mathrm{h}},H_{\mathrm{crit}}\cdot W_{\mathrm{b},\mathrm{crit}})}},\end{array}$$

and an overwash flux widening the barrier,

$$\begin{array}{}\text{(6)}& Q_{\mathrm{ow},\mathrm{b}}=Q_{\mathrm{ow},\max }{\displaystyle \frac{V_{\mathrm{d},\mathrm{b}}}{max(V_{\mathrm{d},\mathrm{b}}+V_{\mathrm{d},\mathrm{h}},H_{\mathrm{crit}}\cdot W_{\mathrm{b},\mathrm{crit}})}}.\end{array}$$

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:

$$\begin{array}{}\text{(7)}& {\displaystyle }{V}_{\mathrm{d},\mathrm{b}}=max\left[\mathrm{0},\left(W_{\mathrm{b},\mathrm{crit}}-W_{\mathrm{b}}\right)\left(H+D_{\mathrm{lagoon}}\right)\right],\text{(8)}& {\displaystyle }{V}_{\mathrm{d},\mathrm{h}}=max\left[\mathrm{0},\left(H_{\mathrm{crit}}-H\right)W_{\mathrm{b}}\right].\end{array}$$

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,

$$\begin{array}{}\text{(9)}& Q_{\mathrm{sf}}=k_{\mathrm{sf}}\left(s_{\mathrm{sf},\mathrm{eq}}-s_{\mathrm{sf}}\right).\end{array}$$

3.2.2 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 $z-x_{\mathrm{b}}\cdot s_{\mathrm{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$\le D_{\mathrm{lagoon}}\le z-x_{\mathrm{b}}\cdot s_{\mathrm{background}}$. We then compute the fine sediment thickness in the back barrier (Fig. 1) as

$$\begin{array}{}\text{(10)}& D_{\mathrm{fines}}=z-x_{\mathrm{b}}{s}_{\mathrm{background}}-D_{\mathrm{lagoon}}.\end{array}$$

In turn, we can geometrically define f as follows:

$$\begin{array}{}\text{(11)}& f={\displaystyle \frac{D_{\mathrm{fines}}}{D_{\mathrm{fines}}+D_{\mathrm{lagoon}}+H}}.\end{array}$$

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., $\mathrm{\Delta }{x}_{\mathrm{s},\mathrm{ow},\mathrm{sf}}\left(\mathrm{2}H+D_{\mathrm{T}}\right))$ deposits on top and/or back of the barrier (Fig. 3).

Figure 3 Effect of lagoon depth (fine sediment fraction) on (a) the overwash flux, (b)  the barrier width, and (c) the shoreline location. More overwash flux is needed to maintain a barrier with a deep lagoon, resulting in barrier drowning for a lagoon depth of > 8 m, where the required barrier overwash flux is greater than the maximum potential overwash flux Q_ow,max.

3.2.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₅₀, 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₀ 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),

$$\begin{array}{}\text{(12)}& k_{\mathrm{sf}}=(\mathrm{3600}\cdot \mathrm{24}\cdot \mathrm{365}){\displaystyle \frac{{\int }_{z_{\mathrm{0}}}^{D_{\mathrm{T}}}\mathit{\varphi }\frac{H(z{)}^{\mathrm{5}}{g}^{\mathrm{5}/\mathrm{2}}}{\mathrm{32}{w}_{\mathrm{s}}^{\mathrm{2}}{z}^{\mathrm{5}/\mathrm{2}}}\mathrm{d}z}{D_{\mathrm{T}}-z_{\mathrm{0}}}},\end{array}$$

where $\mathit{\varphi }=\mathrm{16}{e}_{\mathrm{s}}{C}_{\mathrm{s}}\mathit{\rho }/\left(\mathrm{15}\mathit{\pi }\right({\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho }\left)g\right)$ 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, $H\left(z\right)=H_{\mathrm{s}}\sqrt{\sqrt{g}\cdot T/\mathrm{4}\mathit{\pi }/\sqrt{z}}$. We derive a simple analytical expression of the integrated shoreface response rate,

$$\begin{array}{}\text{(13)}& k_{\mathrm{sf}}=(\mathrm{3600}\cdot \mathrm{24}\cdot \mathrm{365}){\displaystyle \frac{e_{\mathrm{s}}{c}_{\mathrm{s}}{g}^{\frac{\mathrm{11}}{\mathrm{4}}}{H}_{\mathrm{s}}^{\mathrm{5}}{T}_{\mathrm{p}}^{\frac{\mathrm{5}}{\mathrm{2}}}}{\mathrm{960}R{\mathit{\pi }}^{\frac{\mathrm{7}}{\mathrm{2}}}{w}_{\mathrm{s}}^{\mathrm{2}}}}\cdot {\displaystyle \frac{\frac{\mathrm{1}}{\frac{\mathrm{11}}{\mathrm{4}}{z}_{\mathrm{0}}^{\frac{\mathrm{11}}{\mathrm{4}}}}-\frac{\mathrm{1}}{\frac{\mathrm{11}}{\mathrm{4}}{D}_{\mathrm{T}}^{\frac{\mathrm{11}}{\mathrm{4}}}}}{D_{\mathrm{T}}-z_{\mathrm{0}}}},\end{array}$$

where we estimate z₀ 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),

$$\begin{array}{}\text{(14)}& D_{\mathrm{T}}=\mathrm{0.018}{H}_{\mathrm{s}}{T}_{\mathrm{p}}\sqrt{{\displaystyle \frac{g}{{\mathrm{RD}}_{\mathrm{50}}}}}.\end{array}$$

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

$$\begin{array}{}\text{(15)}& s_{\mathrm{sf},\mathrm{eq}}={\displaystyle \frac{\mathrm{3}{w}_{\mathrm{s}}}{\mathrm{4}\sqrt{g{D}_T}}}\left(\mathrm{5}+{\displaystyle \frac{\mathrm{3}{T}_{\mathrm{p}}^{\mathrm{2}}g}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}{D}_{\mathrm{T}}}}\right),\end{array}$$

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

$$\begin{array}{}\text{(16)}& w_{\mathrm{s}}={\displaystyle \frac{Rg{D}_{\mathrm{50}}^{\mathrm{2}}}{\mathrm{18}\cdot {\mathrm{10}}^{-\mathrm{6}}+\sqrt{\frac{\mathrm{3}}{\mathrm{4}}Rg{D}_{\mathrm{50}}^{\mathrm{3}}}}}.\end{array}$$

3.3 Inlet model

Inlets can form along barrier island chains if there is sufficient potential for tidal flow between the lagoon and the open ocean (Escoffier, 1940). In turn, the potential for tidal flow is determined by factors influencing the potential tidal prism (e.g., the proximity of other tidal inlets nearby, the width and depth of the basin and the barrier, the marsh cover) and factors reducing tidal flow (e.g., tidal inlet friction, wave-driven transport into tidal inlets). Once inlets exist, they alter barrier morphodynamics by distributing sediments and enhancing storm overwash potential.

3.3.1 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. We do not open a new inlet if the flow velocity through a new inlet is insufficient (see below). If a new inlet is opened, we place the barrier volume in the flood–tidal delta by increasing the back-barrier location,

$$\begin{array}{}\text{(17)}& \mathrm{\Delta }{x}_{\mathrm{b},\mathrm{breach}}=W_{\mathrm{b}}{\displaystyle \frac{H+D_{\mathrm{inlet}}}{D_{\mathrm{lagoon}}}},\end{array}$$

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.

Figure 4 Analytical solutions to the equations that govern tidal inlet size (Eq. 21) as a function of tidal basin area and offshore tidal amplitude, illustrating that (a) inlets close if a velocity amplitude of 1 m s⁻¹ cannot be maintained, (b) inlet cross-sectional area is dependent on tidal amplitude and tidal basin area, and (c) lagoon and inlet friction reduce the volume of water transported through the inlet, which is less than the potential tidal prism (offshore tidal range multiplied by intertidal area).

3.3.2 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

$$\begin{array}{}\text{(18)}& {\stackrel{\mathrm{̃}}{A}}_{\mathrm{inlet}}={\displaystyle \frac{\mathrm{2}{\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}}{\mathrm{3}}}+{\displaystyle \frac{{\mathrm{2}}^{\frac{\mathrm{2}}{\mathrm{3}}}{F}_{\mathrm{0}}}{\mathrm{6}}}+{\displaystyle \frac{{\mathrm{2}}^{\frac{\mathrm{1}}{\mathrm{3}}}(\mathrm{3}+{\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}^{\mathrm{2}})}{\mathrm{3}{\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}^{\mathrm{2}}{F}_{\mathrm{0}}}},\end{array}$$

where F₀ is

$$\begin{array}{}\text{(19)}& F_{\mathrm{0}}={\left(\begin{array}{l}\mathrm{6}\sqrt{\mathrm{3}}\sqrt{\begin{array}{c}{\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}^{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{-\mathrm{2}}{\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}^{\mathrm{2}}-{\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}^{\mathrm{4}}{\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}^{-\mathrm{2}}+\mathrm{2}{\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}^{-\mathrm{4}}\\ -\mathrm{9}{\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{-\mathrm{2}}-{\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}^{-\mathrm{6}}+\frac{\mathrm{27}}{\mathrm{4}}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{-\mathrm{4}}{\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}^{\mathrm{4}}\end{array}}\\ -\mathrm{27}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}^{-\mathrm{2}}{\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}^{\mathrm{2}}-\mathrm{2}{\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}^{\mathrm{3}}+\mathrm{18}{\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}{\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}^{-\mathrm{2}}\end{array}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}.\end{array}$$

In this formulation, ${\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}=\mathit{\omega }{W}_{\mathrm{b}}/\sqrt{g{a}_{\mathrm{0}}}$ is a resonance nondimensional cross-sectional area, ${\stackrel{\mathrm{̃}}{u}}_{\mathrm{e}}=u_{\mathrm{e}}/\sqrt{g{a}_{\mathrm{0}}}$ is a nondimensional equilibrium velocity, and $\stackrel{\mathrm{̃}}{\mathit{\gamma }}$ is the ratio of the potential tidal prism and the inlet friction,

$$\begin{array}{}\text{(20)}& \stackrel{\mathrm{̃}}{\mathit{\gamma }}={\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{aspect}}\sqrt{\mathit{\omega }{L}_{\mathrm{lagoon}}{W}_{\mathrm{lagoon}}(\mathrm{1}-f_{\mathrm{marsh}})\sqrt{a_{\mathrm{0}}/g}}}{\frac{\mathrm{8}}{\mathrm{3}\mathit{\pi }}{c}_{\mathrm{d}}{W}_{\mathrm{b}}}},\end{array}$$

where the drag coefficient $c_{\mathrm{d}}=g\cdot n^{\mathrm{2}}/D_{\mathrm{lagoon}}^{\frac{\mathrm{1}}{\mathrm{3}}}$ (de Swart and Zimmerman, 2009).

Based on ${\stackrel{\mathrm{̃}}{A}}_{\mathrm{inlet}}$ we determine the dimensional inlet cross-sectional area (m²),

$$\begin{array}{}\text{(21)}& A_{\mathrm{inlet}}=\mathit{\omega }\cdot (\mathrm{1}-f_{\mathrm{marsh}})\cdot L_{\mathrm{lagoon}}\cdot W_{\mathrm{lagoon}}\cdot \sqrt{g{a}_{\mathrm{0}}}{\stackrel{\mathrm{̃}}{A}}_{\mathrm{inlet}},\end{array}$$

and the corresponding water velocity through the inlet,

$$\begin{array}{}\text{(22)}& u={\left({\displaystyle \frac{a_{\mathrm{0}}}{\mathrm{2}g}}\stackrel{\mathrm{̃}}{\mathit{\gamma }}{\stackrel{\mathrm{̃}}{A}}_{\mathrm{inlet}}^{\frac{\mathrm{1}}{\mathrm{2}}}\left[-\mathit{\xi }+{\left({\mathit{\xi }}^{\mathrm{2}}+\mathrm{4}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right]\right)}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

where $\mathit{\xi }=\stackrel{\mathrm{̃}}{\mathit{\gamma }}{\stackrel{\mathrm{̃}}{A}}_{\mathrm{inlet}}^{\frac{\mathrm{1}}{\mathrm{2}}}{\left({\stackrel{\mathrm{̃}}{A}}_{\mathrm{inlet}}-{\stackrel{\mathrm{̃}}{A}}_{\mathrm{H}}\right)}^{\mathrm{2}}$. The inlet area function (Eq. 21) evaluates to the largest cross-sectional area for which u=u_e (1 m s⁻¹ 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.

Figure 5 BRIE without inlets showing barrier response to variations in the initial back-barrier position x_b (see also Ashton and Lorenzo-Trueba, 2018), (a) x_s, and x_b at five instances; (b) barrier width as a function of time. Minima in barrier width drive faster transgression, which in turn results in wider barriers through the accumulation of alongshore sediment.

3.3.3 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),

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{Q}_{\mathrm{s},\mathrm{in}}& =(\mathrm{3600}\cdot \mathrm{24}\cdot \mathrm{365})\cdot k\cdot H_{\mathrm{s}}^{\frac{\mathrm{12}}{\mathrm{5}}}\\ & \cdot T_{\mathrm{p}}^{\frac{\mathrm{1}}{\mathrm{5}}}\cos {\left({\mathit{\varphi }}_{\mathrm{0}}-\mathit{\theta }\right)}^{\frac{\mathrm{6}}{\mathrm{5}}}\sin \left({\mathit{\varphi }}_{\mathrm{0}}-\mathit{\theta }\right),\end{array}\end{array}$$

where k is a constant that is ∼0.06 m^3∕5 s${}^{-\mathrm{6}/\mathrm{5}}$ (Nienhuis et al., 2015), and ϕ₀ 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),

$$\begin{array}{}\text{(24)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{0}}=f^{-\mathrm{1}}\left(X\right),\text{(25)}& {\displaystyle }f\left({\mathit{\varphi }}_{\mathrm{0}}\right)=\left\{\begin{array}{lll}a\cdot h\cdot \left(\frac{\mathrm{4}}{\mathit{\pi }}x+\mathrm{2}\right)& \mathrm{if}& -\frac{\mathrm{1}}{\mathrm{2}}\mathit{\pi }<{\mathit{\varphi }}_{\mathrm{0}}<-\frac{\mathrm{1}}{\mathrm{4}}\mathit{\pi }\\ a+\frac{\mathrm{4}}{\mathit{\pi }}x(\mathrm{1}-h)a& \mathrm{if}& -\frac{\mathrm{1}}{\mathrm{4}}\mathit{\pi }<{\mathit{\varphi }}_{\mathrm{0}}<\mathrm{0}\\ a+\frac{\mathrm{4}}{\mathit{\pi }}x(\mathrm{1}-h)(\mathrm{1}-a)& \mathrm{if}& \mathrm{0}<{\mathit{\varphi }}_{\mathrm{0}}<\frac{\mathrm{1}}{\mathrm{4}}\mathit{\pi }\\ (\mathrm{1}-a)h(\frac{\mathrm{4}}{\mathit{\pi }}x-\mathrm{2})+\mathrm{1}& \mathrm{if}& \frac{\mathrm{1}}{\mathrm{4}}\mathit{\pi }<{\mathit{\varphi }}_{\mathrm{0}}<\frac{\mathrm{1}}{\mathrm{2}}\mathit{\pi }\end{array}\right.,\end{array}$$

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).

3.3.4 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

$$\begin{array}{}\text{(26)}& {\displaystyle }\mathit{\alpha }\left(I\right)=\mathrm{1}-\mathit{\beta }\left(I\right)-\mathit{\delta }\left(I\right),\text{(27)}& {\displaystyle }{a}_{\mathrm{r}}\left(I\right)=\mathrm{0.6}\cdot \mathit{\alpha }\left(I\right),\text{(28)}& {\displaystyle }\mathit{\beta }\left(I\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{10}\cdot I^{\mathrm{3}}}},\text{(29)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{r}}\left(I\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{0.9}\cdot I^{-\mathrm{3}}}},\text{(30)}& {\displaystyle }\mathit{\delta }\left(I\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{3}\cdot I^{-\mathrm{3}}}},\text{(31)}& {\displaystyle }{\mathit{\delta }}_{\mathrm{r}}={\displaystyle \frac{A_{\mathrm{b},\mathrm{downdrift}}\cdot \mathit{\alpha }-A_{\mathrm{b},\mathrm{updrift}}\cdot {\mathit{\beta }}_{\mathrm{r}}}{A_{\mathrm{b},\mathrm{updrift}}}},\end{array}$$

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 and wave momentum flux M_t and M_w,

$$\begin{array}{}\text{(32)}& I={\displaystyle \frac{M_{\mathrm{t}}}{M_{\mathrm{w}}}}\cdot {\displaystyle \frac{W_{\mathrm{inlet}}}{W_{\mathrm{b}}}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}{u}_{\mathrm{eq}}^{\mathrm{2}}{A}_{\mathrm{inlet}}}{\frac{\mathrm{1}}{\mathrm{16}}{\mathit{\rho }}_{\mathrm{w}}g{H}_{\mathrm{s}}^{\mathrm{2}}{W}_{\mathrm{inlet}}}}\cdot {\displaystyle \frac{W_{\mathrm{inlet}}}{W_{\mathrm{b}}}}.\end{array}$$

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_{\mathrm{fld},max}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}{\left(u_{\mathrm{e}}{A}_{\mathrm{eq}}/\mathrm{2}{\mathit{\omega }}_{\mathrm{0}}\right)}^{\mathrm{0.37}}$ (m³). 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 flood–tidal delta is at sea level, we can describe its rate of growth (Fig. 1c) as follows:

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{b},\mathrm{inlet}}}{\mathrm{\Delta }t}}={\displaystyle \frac{Q_{\mathrm{s},\mathrm{in}}\left(\mathit{\delta }+{\mathit{\delta }}_{\mathrm{r}}\right)}{W_{\mathrm{inlet}}{D}_{\mathrm{lagoon}}}},\end{array}$$

change the sediment budget in the littoral zone,

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{s},\mathrm{inlet}}}{\mathrm{\Delta }t}}={\displaystyle \frac{\left(\mathit{\beta }+{\mathit{\beta }}_{\mathrm{r}}\right)Q_{\mathrm{s},\mathrm{in}}}{\left(H+D_{\mathrm{T}}\right)W_{\mathrm{inlet}}}},\end{array}$$

and migrate alongshore in the direction of the littoral drift,

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{\mathrm{\Delta }{y}_{\mathrm{inlet}}}{\mathrm{\Delta }t}}={\displaystyle \frac{Q_{\mathrm{s},\mathrm{in}}\left(\mathit{\alpha }+{\mathit{\alpha }}_{\mathrm{r}}\right)}{A_{\mathrm{b},\mathrm{updrift}}}}.\end{array}$$

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, $\mathrm{\Delta }{H}_{\mathrm{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_{\mathrm{inlet},i}=W_{\mathrm{inlet}}/\mathrm{\Delta }y$, where n_inlet is the number of alongshore cells taken up by inlet i (Fig. 1).

3.4 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:

$$\begin{array}{}\text{(36)}& {\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{s},j}}{\mathrm{\Delta }t}}=D_j{\displaystyle \frac{{\mathrm{\Delta }}^{\mathrm{2}}x}{\mathrm{\Delta }{y}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{s},j,\mathrm{ext}}}{\mathrm{\Delta }t}},\end{array}$$

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_{\mathrm{s},j,\mathrm{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),

$$\begin{array}{}\text{(37)}& D\left(\mathit{\theta }\right)={\displaystyle \frac{k}{(H_{\mathrm{b}}+D_{\mathrm{T}})}}{H}_{\mathrm{0}}^{\mathrm{12}/\mathrm{5}}{T}^{\mathrm{1}/\mathrm{5}}\left[E\left({\mathit{\varphi }}_{\mathrm{0}}\right)\cdot \mathrm{\Psi }\left({\mathit{\varphi }}_{\mathrm{0}}-\mathit{\theta }\right)\right],\end{array}$$

where k is ∼0.06 m^3∕5 s${}^{-\mathrm{6}/\mathrm{5}}$ (Nienhuis et al., 2015) and $\mathit{\theta }\left(j\right)=(x_{\mathrm{s},j+\mathrm{1}}-x_{\mathrm{s},j})/\mathrm{\Delta }y$. Ψ is the angle dependence of the diffusivity (AM06) that we compute as follows:

$$\begin{array}{}\text{(38)}& \mathrm{\Psi }\left({\mathit{\varphi }}_{\mathrm{0}}-\mathit{\theta }\right)={\cos }^{\mathrm{1}/\mathrm{5}}\left({\mathit{\varphi }}_{\mathrm{0}}-\mathit{\theta }\right)\left[\cos {\left({\mathit{\varphi }}_{\mathrm{0}}-\mathit{\theta }\right)}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{6}}{\mathrm{5}}}{\sin }^{\mathrm{2}}\left({\mathit{\varphi }}_{\mathrm{0}}-\mathit{\theta }\right)\right],\end{array}$$

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

$$\begin{array}{}\text{(39)}& E\left({\mathit{\varphi }}_{\mathrm{0}}\right)=\left\{\begin{array}{lll}a\cdot h& \mathrm{if}& -\frac{\mathrm{1}}{\mathrm{2}}\mathit{\pi }<{\mathit{\varphi }}_{\mathrm{0}}<-\frac{\mathrm{1}}{\mathrm{4}}\mathit{\pi }\\ a+(\mathrm{1}-h)& \mathrm{if}& -\frac{\mathrm{1}}{\mathrm{4}}\mathit{\pi }<{\mathit{\varphi }}_{\mathrm{0}}<\mathrm{0}\\ (\mathrm{1}-h)(\mathrm{1}-a)& \mathrm{if}& \mathrm{0}<{\mathit{\varphi }}_{\mathrm{0}}<\frac{\mathrm{1}}{\mathrm{4}}\mathit{\pi }\\ (\mathrm{1}-a)h& \mathrm{if}& \frac{\mathrm{1}}{\mathrm{4}}\mathit{\pi }<{\mathit{\varphi }}_{\mathrm{0}}<\frac{\mathrm{1}}{\mathrm{2}}\mathit{\pi }\end{array}\right.,\end{array}$$

to generate a long-term, wave-climate-averaged shoreline diffusivity for every alongshore location j.

We rewrite the shoreline diffusion Eq. (28) into

$$\begin{array}{}\text{(40)}& \begin{array}{rl}& {\displaystyle \frac{{x_{\mathrm{s}}}_j^{n+\mathrm{1}}-{x_{\mathrm{s}}}_j^n}{\mathrm{\Delta }t}}\\ & ={\displaystyle \frac{D_j}{\mathrm{2}}}{\displaystyle \frac{\left({x_{\mathrm{s}}}_{j+\mathrm{1}}^{n+\mathrm{1}}-\mathrm{2}{x_{\mathrm{s}}}_j^{n+\mathrm{1}}+{x_{\mathrm{s}}}_{j-\mathrm{1}}^{n+\mathrm{1}}\right)+\left({x_{\mathrm{s}}}_{j+\mathrm{1}}^n-\mathrm{2}{x_{\mathrm{s}}}_j^n+{x_{\mathrm{s}}}_{j-\mathrm{1}}^n\right)}{\mathrm{\Delta }{y}^{\mathrm{2}}}}\\ & +{\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{s},j,\mathrm{ext}}}{\mathrm{\Delta }t}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(41)}& \begin{array}{rl}& \left[\begin{array}{ccccc}\mathrm{\ddots }& \mathrm{\ddots }& \mathrm{0}& \mathrm{0}& -{\mathit{\beta }}_j\\ \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -{\mathit{\beta }}_j& \mathrm{1}+\mathrm{2}{\mathit{\beta }}_j& -{\mathit{\beta }}_j& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{\ddots }\\ -{\mathit{\beta }}_j& \mathrm{0}& \mathrm{0}& \mathrm{\ddots }& \mathrm{\ddots }\end{array}\right]\left[\begin{array}{c}\mathrm{⋮}\\ {x_{\mathrm{s}}}_{j-\mathrm{1}}^{n+\mathrm{1}}\\ {x_{\mathrm{s}}}_j^{n+\mathrm{1}}\\ {x_{\mathrm{s}}}_{j+\mathrm{1}}^{n+\mathrm{1}}\\ \mathrm{⋮}\end{array}\right]\\ & =\left[\begin{array}{c}\mathrm{⋮}\\ \mathrm{⋮}\\ {x_{\mathrm{s}}}_j^n+{\mathit{\beta }}_j\left({x_{\mathrm{s}}}_{j+\mathrm{1}}^n-\mathrm{2}{x_{\mathrm{s}}}_j^n+{x_{\mathrm{s}}}_{j-\mathrm{1}}^n\right)+\mathrm{\Delta }{x}_{\mathrm{s},j,\mathrm{ext}}\\ \mathrm{⋮}\\ \mathrm{⋮}\end{array}\right],\end{array}\end{array}$$

where ${\mathit{\beta }}_j=D_j^n\mathrm{\Delta }t/\mathrm{2}/\mathrm{\Delta }{y}^{\mathrm{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_{\mathrm{s},j,\mathrm{ext}}$ (m) can be described by

$$\begin{array}{}\text{(42)}& \mathrm{\Delta }{x}_{\mathrm{s},j,\mathrm{ext}}=\mathrm{\Delta }{x}_{\mathrm{s},j,\mathrm{ow},\mathrm{sf}}+\mathrm{\Delta }{x}_{\mathrm{s},j,\mathrm{inlet}},\end{array}$$

representing offshore and onshore sediment fluxes that can erode and accrete the shoreline and flood–tidal delta deposition that acts as a littoral sink.

Figure 6 Example model run showing (a) a barrier island including tidal inlet through time, (b) inlet facies after approximately 1300 years, (c) inlet location through time (such that the slope of the line represents the migration rate), (d) average inlet migration rate, (e) the average fraction of alongshore sediment brought into the inlet transported to the downdrift coast (β), the flood–tidal delta (δ), and the barrier itself (α), and (f) the transgressive sediment flux due to overwash and due to flood–tidal delta deposition. Model code and animated model output can be found in Nienhuis and Lorenzo-Trueba (2019) and Nienhuis (2019).

The shoreline model is unconditionally stable and second-order 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.

3.5 Other moving boundaries

At the end of each time step, we update the shoreface toe position x_t,

$$\begin{array}{}\text{(43)}& x_{\mathrm{t}}^{n+\mathrm{1}}=x_{\mathrm{t}}^n+\mathrm{\Delta }{x}_{\mathrm{t},\mathrm{ow}},\end{array}$$

back-barrier location (x_b),

$$\begin{array}{}\text{(44)}& x_{\mathrm{b}}^{n+\mathrm{1}}=x_{\mathrm{b}}^n+\mathrm{\Delta }{x}_{\mathrm{b},\mathrm{ow}}+\mathrm{\Delta }{x}_{\mathrm{b},\mathrm{inlet}}+\mathrm{\Delta }{x}_{\mathrm{b},\mathrm{breach}},\end{array}$$

and barrier height (H),

$$\begin{array}{}\text{(45)}& H^{n+\mathrm{1}}=H^n+\mathrm{\Delta }{H}_{\mathrm{ow}}+\mathrm{\Delta }{H}_{\mathrm{inlet}},\end{array}$$

independently for all alongshore locations j and run another time step.

3.6 Model output

After a model simulation (typically 10 kyr) we obtain shoreline, back-barrier, and shoreface morphodynamics for different scenarios given by, for example, SLR rates, wave climates, and tidal conditions. One aspect of particular interest, and the primary motivation for this model, is the transgressive flux due to inlet activity. We define a ratio F,

$$\begin{array}{}\text{(46)}& F={\displaystyle \frac{Q_{\mathrm{inlet}}}{Q_{\mathrm{overwash}}+Q_{\mathrm{inlet}}}},\end{array}$$

where $Q_{\mathrm{overwash}}=\mathrm{\Delta }y\cdot \sum _y{Q}_{\mathrm{ow},\mathrm{b}}$ (m³ yr⁻¹), and Q_inlet is the along-coast average transgressive sediment flux by inlet formation and flood–tidal delta deposition for all inlets,

$$\begin{array}{}\text{(47)}& \begin{array}{rl}{Q}_{\mathrm{inlet}}& ={\displaystyle \frac{\mathrm{\Delta }{x}_{\mathrm{b},\mathrm{breach}}}{\mathrm{\Delta }t}}{\displaystyle \frac{\left(H+D_{\mathrm{inlet}}\right)}{W_{\mathrm{b}}}}\\ & +\sum _{i=\mathrm{1}..n_{\mathrm{inlet}}}\left(\mathrm{1}-\mathit{\beta }-{\mathit{\beta }}_{\mathrm{r}}\right)Q_{\mathrm{s},\mathrm{in},i}.\end{array}\end{array}$$

F quantifies the fraction of the total transgressive flux due to inlets and can range from 0 to 1.

3.7 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

4.1 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).

Figure 7 Model volume (barrier volume and offshore deposits) relative to the analytically determined model volume for (a) averaged for different time steps and alongshore grid lengths, as well as (b) through time, showing mass conservation.

4.2 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⁻¹, 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

5.1 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 cross-sectional 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).

5.2 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 cross-shore 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⁻¹ (progradation) to +10 m yr⁻¹ (erosion). Even though the overall trajectory is a result of the sea level history and the passive inundation of the main (non-barrier) 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.

Figure 8 Comparison of (a) transgressive fluxes, (b) shoreline erosion rates, and (c) shoreline locations for the BRIE model without inlets to the model with inlets. Inlets induce stochasticity resulting from inlet opening, growing, migrating, and closing.

5.3 Sensitivity to grid resolution and time step

We investigated the sensitivity of the model output (Q_overwash, Q_inlet, and F) by varying the grid resolution and time resolution and holding all other parameters constant. In general, we find that these fluxes vary approximately ∼20 % between different settings (Fig. 9). These deviations appear only in simulations that include inlets and are likely caused by a sensitivity to small perturbations such as random wave angles. For example, comparing multiple simulations with equal settings, including grid and time, we obtain a variability in F (Fig. 10), with a standard deviation of 0.025. Sensitivity to grid spacing and time steps can also be caused by the discretization of inlet migration rates and distances (Eq. 35).

Figure 9 Tidal-induced transgressive sediment flux, storm-overwash-induced transgressive sediment flux, and the ratio F as a function of grid length and time step normalized by their values at a time step of 0.05 years and 100 m alongshore grid discretization.

Figure 10 Variability in the inlet fraction of the transgressive sediment flux F for varying grid lengths Δy and time steps Δt and for constant Δy and Δt, highlighting the fraction of the variability in F due to inherent model variability.

6 Model evaluation

It is challenging to evaluate long-term barrier island models against natural examples. Given their erosional nature, long-term 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 most effective for intermediate back-barrier depths, whereas overwash volumes are highest for deeper back-barrier depths. 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).

Figure 11 Fraction F for varying model parameters showing the mean and the range obtained for three different experiments. In red is the value of the specific parameter that is held constant across the other simulations.

7 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 cross-shore 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.