Model Sensitivity Analysis for Coastal Morphodynamics: Investigating Sediment Parameters and Bed Composition in Delft3D

Abstract

Numerical simulation of sediment transport and subsequent morphological evolution rely on accurate parameterizations of sediment characteristics. However, these data are often not available or are spatially and/or temporally limited. This study approaches the problem of limited sediment grain-size data with a series of simulations assessing model sensitivity to sediment parameters and initial bed composition configurations in Delft3D, leading to improved modeling practices. A previously validated Delft3D sediment transport and morphology model for Dauphin Island, Alabama, USA, is used as the benchmark case. A method for the generation of representative sediment grain sizes and their spatially varying distributions is presented via end-member analysis of in situ surficial sediment samples. Derived sediment classes and their spatial distributions are applied to two sensitivity case simulations with increasing bed composition complexity. First, multiple sediment classes are applied in a single fully mixed layer, regardless of sediment type. Second, multiple sediment classes are applied in a thin, fully mixed transport layer with underlayers containing only the non-cohesive sediment classes below. Simulations were carried out in a probabilistic, Delft3D MorMerge configuration to capture long-term morphology change for 10 years. We found there is sensitivity to the inclusion of additional sediment classes and sediment distribution made evident in bed level and morphology change. Inclusion of highly mobile fine sediments altered model results in each sensitivity case. The model was also found to be sensitive to initial bed composition in terms of bed level and morphology change, with notable differences between sensitivity cases on decadal timescales, indicating an armoring effect in the second sensitivity case, which used the transport and underlayer bed configuration. The results of this study offer guidance for numerical modelers concerned with sediment behavior in coastal and estuarine environments.

1. Introduction

Sediment transport processes impact key coastal management interests [1,2], including restoration decisions [3,4,5], dredging impacts [6,7], habitat suitability [8], aquaculture [9,10], and submerged aquatic vegetation [11]. Assessing how natural and human drivers may impact sediment dynamics in coastal systems requires detailed observations and complex predictive models. Process-based numerical models such as Delft3D can be used to simulate hydrodynamics, sediment transport, and resultant morphological evolution in estuarine and coastal environments. Accurately representing wave- and current-driven sediment (re)suspension, transport, and deposition is critical for simulating morphological evolution. Delft3D model errors may arise from inaccurate initial conditions (e.g., bathymetry, topography), boundary conditions [12], or input parameters such as bottom roughness [13]. Input parameters are often “tuned” or adjusted during the calibration process to improve the accuracy of model output in comparison with observations. The need for calibration may be reduced when sufficient input data are available to accurately describe the physics of the system. However, these data are not always available or are spatially and/or temporally limited. As increased sensitivity to model input leads to increased uncertainty in model outputs [14], understanding the sensitivity of the model to input parameters and initial conditions improves model assessments and design. Further, an analysis of model flow component celerity values [15] suggests that in contrast to uncertainties in flow velocities or water depths, uncertainties and errors in both the bed elevation and the bed sediment composition are retained in the model domain for the duration of the simulation. In other words, while initial error in hydrodynamics may dissipate over a sufficient simulation time, errors in bed composition and bed level will persist.

Delft3D (D3D) [16,17] is a widely used hydrodynamic, sediment transport, and morphology model applied in a range of coastal, estuarine, and fluvial environments across a range of timescales. Sediment transport models, such as Delft3D, XBeach, and the fluvial surface-based formulation model presented by Hu et al. [18], are capable of modeling sediment transport when multiple grain sizes are present. The improvement of such computationally intense models, either by efficient time-stepping schema [19] or accelerated morphological activity [20,21], has been a decades-long pursuit. Model efficacy may also be improved by deep understanding of model sensitivity. Evaluation of D3D sediment transport and morphology module sensitivity to various model inputs, parameters, and coefficients has formed a broad, yet non-exhaustive knowledge base for numerical modelers. Welsch [22] evaluated Manning’s friction coefficient and asymmetry coefficient, finding little impact on rip channels, but sensitivity of nearshore bar formation to asymmetry coefficients. Another study [23] assessed sensitivity to morphological acceleration factor, active layer thickness, and wind waves, and critically showed sensitivity (overprediction) of sediment loads for spatially uniform bed composition during model spin-up. Model parameters of bed porosity [24] and the constitutive relation for bed slope [25] have also been shown to have strong controls on sediment processes and morphological evolution. Of particular interest in barrier island settings, a sensitivity analysis of the threshold minimum water depth parameter (minimum depth for sediment computations) [26] showed an influence of this parameter on the interaction of barrier islands and adjacent inlets, with a pronounced influence on the erosion of barrier island termini. Further, it was shown that island tip erosion influenced the formation and size of secondary inlet channels [26]. This finding has potential implications for the sensitivity of overwash channel and breach formation in higher-energy regimes and highlights the potential for second-order effects of model sensitivity to inputs and parameterization. Boechat-Albernaz et al. [27] focused on the sensitivity of long-term, decadal, morphology-to-wave orbital velocity parameterization (necessary in 2D depth-averaged models) finding that a skewness–asymmetry parameterization led to more stable and realistic shoreline change, opening the possibility of improved long-term morphological modeling of complex coastlines.

Sediment grain size is a known control for increasing beach-face slope trends with increasing sediment grain size and decreasing beach exposure to wave energy [28], shoreline variability [29], and riverbed stratigraphic evolution [19]. Beach-face slope, dune heights, and other important subaerial morphological thresholds are not typically the target outputs modeled by D3D, which does not by default include swash or aeolian processes needed to move sediment onshore to model dry points, but such morphological thresholds can be handled by a model such as XBeach [30], specifically with respect to storm impacts. A combined field and modeling study [31] evaluated the need for spatial accuracy in grain-size distributions for the XBeach model, a factor that is as yet unevaluated for D3D and which this study aims in part to address. Koktas [31] showed the ability of that model to sort grain sizes given accurate initial spatial distributions. Sediment grain size has also been shown to be critical to reproducing observed overwash extent [32] in the XBeach model, which is likely due in part to grain-size control on beach-face slope. Mickey et al. [33] established a high degree of XBeach model sensitivity to antecedent beach slope and topographies (i.e., dune crest height, dune toe, dune volume).

Barrier island overwash, beach slope, and shoreline position each contribute to the overall morphology (topobathy) of a coastal system. Dune dynamics also shape an island and are a critical driver of barrier island evolution [34] at timescales relevant to coastal modelers and managers. Nearshore sediment transport processes may influence subaerial island morphology, with shoal changes influencing dune dynamics [35]. Further complexities arise with barrier island response to sea-level rise (SLR) and increased storm frequency [34,36,37], which are difficult to reproduce or forecast with any single model.

The ability to forecast long-term multi-decadal evolution of coastal systems under present-day and future forecast conditions of sea-level rise (SLR) and storm frequency/intensity more often requires a “model train” framework. Such frameworks [38,39] connect disparate models, each handling different timescales or energy regimes, and loosely couple them with the outputs from each serving as inputs of the others (e.g., fair-weather shorelines and bed-level change output of D3D as initial pre-storm elevations in XBeach, or post-storm bed levels from XBeach serving as initial elevations for D3D). The need for model train frameworks further motivates this study, as D3D can produce antecedent morphologies for subsequent physical processes such as aeolian transport, dune development, and total barrier island response to storms.

This study assessed the influence of input sediment parameterization on coastal sediment transport and resultant morphological evolution using a two-dimensional Delft3D model developed for Dauphin Island, Alabama (AL) [40,41]. This region (Figure 1) was selected due to the diversity of coastal features (semi-enclosed bay, barrier islands, estuaries) and breadth of collected field data.

Additionally, the setting of Dauphin Island and the adjacent coastal region provide areas known to be highly morphologically active with well-documented historical landform changes, allowing a broad evaluation of sensitivity of the morphological response of the seafloor and subaerial landforms (e.g., sand spits and shoals). The seafloor surrounding Dauphin Island is spatially dynamic [42], with tidally dominant tidal inlets (Main Pass and Petit Bois Pass, Figure 1) [43] reworking the seafloor and a morphologically active Ebb-Tidal Shoal (Figure 1) worked by the joint action of waves and currents [44,45]. Historically, subaerial shoals have formed on the Ebb-Tidal Shoal [43] and migrated in the direction of net sediment transport [46,47] (i.e., north-north-west, and westward [48,49,50], reforming or joining Pelican Island [43] and ultimately supplying sediment to Dauphin Island. The shoreface of Dauphin Island is also morphologically active, being subject to alongshore littoral sediment transport, wave action, overwash, and breaching by storms [50,51,52]. We present a novel methodology that uses collected field data and grain-size analysis to establish a series of appropriate sediment parameters for the model. The model was run in probabilistic MorMerge [20] mode to assess the influence of the input sediment parameterization and bed configurations on decadal-scale cohesive and non-cohesive sediment transport and resulting morphological change. The results of this sensitivity analysis can be used to improve the understanding of sediment parameterization and bed-level dynamics and to recommend optimal strategies for developing representative sediment inputs. Further, the use of the P-J model for Dauphin Island as a benchmark combined with in situ sediment data-derived model input allows for a novel holistic analysis for both a narrow quantitative and broad qualitative assessment of system benefits and effects produced by common modeling practices when ample sediment data are available.

2. Materials and Methods

2.1. Model Formulation

Delft3D (D3D) (developed by Deltares; see [16]) is an integrated modeling system capable of simulating hydrodynamics, waves, sediment transport, and morphological evolution using a suite of modules that can be coupled in a variety of configurations [16]. D3D is commonly applied in coastal, estuarine, and riverine environments.

Sediment classes are the ranges of sediment types and grain sizes supplied by the user to the transport module. Each sediment class is parameterized by sediment type (cohesive “mud” or non-cohesive “sand”), median sediment diameter (d50), specific density, and dry bed density. Additional parameters are required in the case of cohesive sediments, including saline and fresh water settling velocities. The initial thickness of each sediment class at the bed must also be specified and may be uniform or spatially varying. The model may accept any number of cohesive and non-cohesive sediment classes and solves the mass balance equation for each. Increasing the number of sediment classes leads to longer wall-clock times for model simulations if computational resources are unchanged.

Delft3D solves unsteady shallow-water equations and shallow-water assumptions for vertical momentum for two-dimensional (2D) depth-averaged simulations [16]. The transport portion, as described by [16], applies the advection–diffusion equation in the finite-volume approximation form (Equation (1)).

$${\displaystyle \frac{\partial \left[hc\right]}{\partial t}}+{\displaystyle \frac{\partial \left[hUc\right]}{\partial x}}+{\displaystyle \frac{\partial \left[hVc\right]}{\partial y}}+{\displaystyle \frac{\partial \left[\omega c\right]}{\partial \sigma }}=h \left[{\displaystyle \frac{\partial }{\partial x}}\left(D_H{\displaystyle \frac{\partial c}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial x}}\left(D_H{\displaystyle \frac{\partial c}{\partial x}}\right)\right]+{\displaystyle \frac{1}{h}}{\displaystyle \frac{\partial }{\partial \sigma }}\left[D_V{\displaystyle \frac{\partial c}{\partial \sigma }}\right]+hS$$

(1)

For further definition of the standard form above, see [16]. Of particular interest to this analysis is the final term ($hS$), in which $S$ captures all source and sink terms per unit area. Source and sink terms depend on settling velocity, vertical diffusivity ($D_V$), bed shear stress, bottom orbital velocity, and sediment concentration, and vary with sediment type (i.e., cohesive or non-cohesive), sediment grain size (i.e., sediment diameter, d50), and sediment availability (i.e., sediment distribution and relative abundance).

For cohesive sediments (mud), the source/sink term for a given sediment class, $S^{\left(l\right)}$, is described as per Parthenides–Krone [53] as the balance of deposition and erosion (Equation (2)).

$$S^{\left(l\right)}=E-D$$

(2)

Erosion ($E$), given by Equation (3), is a function of bottom shear stress due to currents and waves (${ \tau }_{cw}$), critical shear stress for a given sediment class (${\tau }_{cr\left(l\right)}$), a user-tunable erosion parameter ($M$), and tunable erosion power (n) (default = 1).

$$E=M max(0,{\displaystyle \frac{{ \tau }_{cw}}{{\tau }_{cr\left(l\right)}}}-1{)}^n$$

(3)

Deposition ($D$), given by Equation (4), is the product of settling velocity (${\omega }_s$), sediment concentration at the bed ($c_b$), and a dimensionless reduction factor (Γ). Γ is itself a function of modeled bed shear stress and critical shear stress for a given sediment diameter.

$$D={\omega }_s{c}_b\Gamma$$

(4)

For non-cohesive sediments (sand), transport is calculated in accordance with the formulations of Van Rijn [54]. Suspended- and bed-load transports are distinguished by a reference height a, with suspended load being above and bed load being below that height, respectively. Source and sink terms for suspended non-cohesive sediments (sand) are given by Equations (5) and (6):

$$\mathrm{S}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}=c_a{\displaystyle \frac{D_v}{∆z}}$$

(5)

$$\mathrm{Sink}=c_{kmx}\left({\displaystyle \frac{D_v}{∆z}}+{\omega }_s\right)$$

(6)

where $D_V$ is the vertical eddy diffusivity [16], $∆z$ is the distance (vertically) from the Van Rijn reference level to the center of reference cell, and ${ c}_{kmx}$ is the mass concentration of the sediment fraction in question in the reference cell. The calculation of $c_a$ is given by Equation (7).

$$C_a=f_{sus}\eta 0.015{\rho }_s{\displaystyle \frac{d50}{a}}{\displaystyle \frac{T_a}{D_{\ast }}}$$

(7)

The terms $T_a$ and $D_{\ast }$ are dimensionless shear and sediment diameter, respectively, ${\rho }_s$ is the sediment density, $f_{sus}$ is a user-tunable calibration term (default = 1), and $\eta$ is the relative availability of the sediment fraction at the bed.

Finally, the sediment source/sink term for bed-load transport, $S_b$ [54], is given by Equation (8).

$$\left|S_b\right|=\eta 0.006{\rho }_s{\omega }_s\left(d50\right)M^{0.5}{M}_e^{0.7}$$

(8)

The trailing terms $M$ and $M_e$ are the mobility due to waves and currents and excess mobility, respectively. Each is a function of sediment diameter, sediment density, water density, depth-averaged velocity, and orbital velocity. In addition to the aforementioned factors, there is also a factor of critical depth-averaged velocity for initiation of motion. Equation (8) gives the rate of bed-load transport, while the direction vector is derived from the part of transport due to waves and the part due to currents.

2.2. Benchmark Model

In the present study, we utilized and built upon a previously published and evaluated Delft3D FLOW/WAVE model developed for the barrier island platform surrounding and including Dauphin Island, AL [40,41], hereafter referred to as the P-J model. The spatial extent covers the area from Fort Morgan Peninsula, AL to the east, Petit Bois Island, Mississippi (MS) to the west, Mobile Bay to the north, and the 30 m depth contour to the south (Figure 1). As the open coast setting is considered well mixed, the P-J model was developed in a two-dimensional (2D), depth-averaged configuration to reduce computational cost. The P-J model was developed as a framework to evaluate sediment transport and morphological changes under both deterministic and probabilistic oceanographic conditions. Simulating littoral sediment transport on monthly timescales and decadal morphological change under fair-weather, quiescent conditions was the purpose. The latter decadal-scale simulation was the focus of this study, and thus P-J model output for a ten-year MorMerge simulation, modeled year-by-year, served as the benchmark simulation for this analysis: model output from this simulation is hereafter referred to as the Benchmark case. The model bed level was initialized using a 2015 digital elevation model (DEM) [40,41]. Morphology was allowed to evolve continuously, though sediment availability (bed composition) was reset after each year of the simulation period to ensure numerical stability, as no fluvial or other sediment source was applied as a boundary condition. A morphological acceleration factor of 365 was used, with each year of morphology change being simulated by a 30 h period (6 h of spin-up without bed updating plus 24 h of simulation time with bed updating). A time step of 0.375 min was used.

The P-J model used a single sediment class with a single representative d50 of 200 microns (μm), with associated reference density, dry bulk density, and settling velocity. This single sediment class was applied in one fully mixed layer. An initial bed-layer thickness of 10 m was chosen to provide ample sediment supply on annual timescales. This single non-cohesive sediment class in a fully mixed bed layer is among the simplest parameterizations of sediment in D3D and serves as the benchmark against which sensitivity to sediment parameterization and bed configuration is measured. A detailed description of various initial and boundary conditions, as well as details on grid resolution, input elevations, and coupling configuration implemented in the P-J model are available in [41], and morphological and sediment initial conditions are summarized in [42]. Additional details of model setup and forcing are provided in Appendix A. Validation and model skill information for the P-J model are reiterated in Appendix B.

In addition to the ten-year MorMerge model output, which served as the benchmark in this study, deterministic P-J model output was also available from a one-month simulation (simulation period from 19 June 2015 through 19 July 2015), which included morphological evolution modeled in real time.

2.3. Suite of Sensitivity Tests

Multiple sediment classes (sediment type and diameter), the distribution of each sediment class, and the configuration of the bed layer may be specified in D3D. The default configuration, as applied in the Benchmark simulation, is that of a fully mixed bed in a single layer, with a single sediment class. Two sensitivity cases were developed for this study, which were designed to progressively increase model complexity with each test, first by adding multiple sediment classes and second by increasing the complexity of the bed configuration.

In the case of multiple sediment classes in a single fully mixed layer, the total bed thickness was taken as the sum of the initial thickness of each sediment class and all sediments were available for mobilization and transport at the start of the simulation in proportion to the initial ratio in the bed. No bookkeeping of sediments was performed by the model in this case. In a second bed configuration, there was an initial overlying transport layer, with bed layers underlying. The concentrations of individual sediment classes in the transport layer and underlayers were tracked throughout the simulation in this configuration.

Figure 2 shows the arrangement of these layers and the bookkeeping that took place in the second bed configuration. Sediment was only initially available for transport from the transport layer. As sediment is eroded, sediment from the underlayers replenishes the transport layer. As sediment is deposited, additional sediment is sequestered in the underlayers. In addition to the transport layer and underlayers, this configuration contains a base layer into which sediments are sequestered if a maximum number and thickness of underlayers is reached.

Figure 2 also presents the various model configurations applied in this analysis. Each case was run in the probabilistic D3D MorMerge mode for a period of ten years. The P-J model Benchmark case (denoted A in Figure 2) uses a uniformly well mixed bed with a single sediment class and a spatially uniform initial thickness of 10 m. The sediment type was set to “sand” and a d50 of 200 μm was used.

The first sensitivity case, “Single Layer, Four Class” (B in Figure 2, referred to as “SL4C”), maintains the single fully mixed layer of the Benchmark case and the overall initial thickness of 10 m in a single bed-layer. Four sediment classes were applied (development, characteristics, and spatial distribution of sediment classes are described in Section 2.4 and Section 2.5), including one 6 μm mud class. Where observationally derived sediment classes were not present in the domain, the 200 μm sand class made up the entire initial thickness. This test case was designed to evaluate overall morphological sensitivity to the inclusion of multiple spatially varying sediment classes.

The second sensitivity case, “Underlayer, Four Class” (C in Figure 2) (referred to as “UL4C”), maintained the four sediment classes and overall initial thickness of SL4C. However, the configuration of the bed was altered to include a transport layer with underlayers. A 0.5 m transport layer was applied, and the remaining thickness was arbitrarily initially divided equally into two underlayers. The 6 μm mud class was only present in the transport layer. Only sand classes were initially in the underlayers. Each class contributed to the total thickness of each layer, as previously described. Where derived sediment classes were not present in the domain, the 200 μm sand class made up the entire initial thickness. The UL4C case was designed for evaluation of model sensitivity to bed configuration (i.e., limitation of fine sediments to the transport layer) and to highlight possible armoring effects of fine sediments by coarser sediments.

2.4. Development of Sediment Classes

When multiple sediment classes of type and grain size are introduced, computational resources are taxed, as the model solves Equation (1) for each sediment class individually. An effort was made to generate sediment grain sizes, distributions, and abundance that were informed by the study area, while also limiting computational complexity. Grain-size data for over 300 surficial sediment samples (at 240 sample locations) were used to develop a model input layer for seabed characterization. The geographic location, type of sampling device, and physical setting at time of collection for the 240 surficial samples has been previously reported in [55]. The east–west extent of samples ranged from Main Pass and the Ebb-Tidal Shoal to the eastern half of Petit Bois Pass. North to south, the samples were collected within approximately 6 km (km) from either the fore- or backshore at depths between 0 m and 10 m. Each sample reflected sediment generally collected within the uppermost two centimeters of the seabed. Samples were collected with approximately 8% spatial replication.

Detailed sediment texture was measured on each sample. A summary of sediment texture is available in Ellis et al. [55], along with complete analytical procedures. Grain or particle size distribution was measured at least six times (N ≥ 6) on each sample using a laser diffraction system (Coulter LS13-320; Beckman Coulter, Inc., Brea, CA, USA). Particle size distributions were placed into 92 size-class bins (M), with each bin representing the relative volume of that size class to a unit volume of sample. The particle size distributions for each sample were imported into Matlab (version 2018a) (Mathworks, Inc., Natick, MA, USA) as an M × N matrix (92 × 6). Means and standard deviations were computed bin-wise to produce average particle size distribution (92 × 1) (n = 6). The resulting bin means were summed and normalized back to unit volume to be consistent with original analytical results. Implementing 92 size classes in the P-J model would be computationally expensive, so size-class data were reduced using an end-member analysis (EMA) method implemented with the AnalySize package designed for Matlab [56]. We utilized a parametric end-member formulation using general Weibull functions to generate and then optimize a constrained number of end-member (EM) curves (2 to 7 end-members). We evaluated the appropriate number of end-members (EMs) using goodness-of-fit metrics (R2 and angular deviation; see [56]). We determined qualitatively that four EMs provided statistically robust estimates of all sediment samples while greatly limiting the computation complexity that would be passed to the P-J model (see Supplementary Materials for further details).

2.5. Application of Sediment Classes to Delft3D

The median diameter (d50, in μm) for each of the four sediment classes derived by EMA was used to develop input size-class layers for the model. The first sediment class (hereafter referred to as EM1) was a mud class (d50 = 6 μm), while the second, third, and fourth (EM2, EM3, and EM4 respectively) were distinct sand classes (d50 = 200, 350, and 570 μm, respectively). For each sediment sample, output from the EMA provided each sediment class as a fraction of unity. However, prior to generating geospatial size-class input layers for the model, EM data were adjusted to clearly distinguish cohesive and non-cohesive samples. Dalyander et al. [57] and Soulsby and Whitehouse [58] suggest that cohesive behavior is applicable when clay content exceeds 7%. Using EM1 as a proxy for clay content, EM1 values were adjusted to 0% for all samples with initial EM1 abundance less than 7%. All samples were subsequently renormalized to ensure the sum of all EMs was 100%. The updated EM percentages were combined with geospatial data (WGS84 latitude and longitude from [55]) to produce four separate georeferenced vector datasets. Each sediment class dataset was loaded into ArcGIS Pro v2.5 (ESRI, Inc.) and interpolated using the inverse distance weighted (IDW) method (Spatial Analyst Toolbox) to produce a raster surface of the relative abundance of each sediment class (Figure 3). The IDW interpolation method was applied for the simplicity of the interpolation scheme and ability to account for alongshore variability and gradients in sediment availability along cross-shore transects.

A map of EM1 (Figure 3A) shows the bulk of the cohesive sediments are focused behind the island. A small amount of fine-grained sediment is found in other protected regions (e.g., east of Pelican Island) or offshore, likely at depths greater than what fair-weather waves influence. In contrast, sand end-members are prevalent throughout the entire sampled region. The fine-grained sand, EM2 (Figure 3B), is most prevalent in the shallow regions fronting Dauphin Island, while the medium sand, EM3 (Figure 3C), is concentrated around the flanks of modern ebb-tidal deltas (Petit Bois Pass and Main Pass). The coarse sand, EM4 (Figure 3D), is located along the shoreline of Little Dauphin Island, adjacent to flood channels that connect the Gulf of Mexico with both Mississippi Sound and Mobile Bay.

The interpolated spatial relevance percentages of the end-members derived by EMA were used to introduce spatial variability of sediment distribution in the D3D model sensitivity cases. For a given thickness of a bed layer, regardless of whether that layer were a single fully mixed, transport layer, or underlayer, the contributing thickness of each class was normalized by that end-member’s percentage relevance. For example, if for a point in a 10 m-thick layer, EM1 showed 20% relevance, that translated to a 2 m thickness of the 6 μm mud class at that point, and so on with the remaining EMs, adding up to a total thickness of 10 m. For regions outside the area where surficial samples were available for which it was not possible to derive end-member prevalence, the entire thickness was set uniformly to 200 μm, keeping the sensitivity cases analogous to the Benchmark case.

Parameters for each sediment class were prescribed and are presented in

Table 1. Sediment parameters for each sediment class. Sediment type, D50 (sediment diameter, in μm), specific density, (${\rho }_s$ in kg/m³), dry bed density (${\rho }_{dbd}$ in kg/m³), and saline settling velocity (${\omega }_s$ in m/s) were prescribed in model input.

Sediment Class Name	Type	D50 (μm)	${\mathit{\rho }}_{\mathit{s}}$ (kg/m3)	${\mathit{\rho }}_{\mathit{d}\mathit{b}\mathit{d}}$ (kg/m3)	${\mathit{\omega }}_{\mathit{s}}$ (m/s) 1
EM1	Cohesive	6	2500	700	4.35 × 10−5
EM2	Non-cohesive	200	2500	1500	-
EM3	Non-cohesive	350	2500	1500	-
EM4	Non-cohesive	570	2500	1500	-

¹ Settling velocity is only mandatory to prescribe for cohesive (mud) sediment classes.

. A specific density of 2500 kg per meter cubed (kg/m³) was used across all classes. A dry bed density of 1500 kg/m³ was used for all sand classes, and a value of 700 kg/m³ was used for the mud class. It was also necessary to derive a saline settling velocity for the mud class, which was achieved the methods of [59,60]. The resultant settling velocity applied was 4.3493426 × 10⁻⁵ m/s. A reference density of 1600 kg/m³ was uniformly applied. The transport formula used for each sand class was as described by van Rijn [54] (Equations (5)–(8)), while the mud sediment class was treated with the Partheniades–Krone formulations (Equations (2)–(4)) [53].

2.6. Data Analysis Methodologies

To assess morphological sensitivity, several methods are applied. First, a simple difference between model end-state bed-levels (${DPS}_{final}$) and initial bed-level (${DPS}_{initial}$) is calculated according to Equation (9).

$$\Delta DPS={DPS}_{final}-{DPS}_{initial}.$$

(9)

Maps of $\Delta DPS$ show deposition (positive values) and erosion (negative values) for a given sensitivity case. When considering the relative difference in end-state bed levels between any two cases, this is calculated as:

$${∆Z}_f={DPS}_{final}^k-{DPS}_{final}^j$$

(10)

Maps of ${∆Z}_f$ show where the end state of one sensitivity case (k) is shallower (positive values) or deeper (negative values) than the end state of another case (j).

A series of 14 nearshore regions (following the methods of [40]) following the lines of the D3D grid were selected and further subdivided into blocks falling over the near foreshore and back-barrier regions of Petit Bois Island, Dauphin Island, and the western portion of Fort Morgan Peninsula. Little Dauphin Island and portions of the back-barrier of Dauphin Island were excluded. Also captured were the inlets between each of the land masses, where tidally driven flow velocities are expected to influence the seafloor. Figure 4 shows these regions outlined with heavy black lines, with the subdividing blocks outlined with thin gray lines.

The areas of each block were derived from model grid spacing. $\Delta Volume$ in cubic meters was then calculated for each block, as given by Equation (11):

$$\Delta Volume=area\ast \sum \Delta DPS$$

(11)

where $\Delta DPS$ (calculated as in Equation (9)) is limited to a given block. As with $\Delta DPS$, mapped $\Delta Volume$ values within blocks show sediment volume gain (deposition—positive values) and sediment volume loss (erosion—negative values) for a given case and additionally serve to highlight these differences, especially in regions where $\Delta DPS$ may be apparently minimal.

A percentage difference in volume change (${Volume}_{diff}$) was also calculated for all blocks to quantify the sensitivity of $\Delta Volume$ at the bed. This was taken as the normalized relative difference between two end volume states, as given by Equation (12).

$${Volume}_{diff}={\displaystyle \frac{{\Delta Volume}_k-{\Delta Volume}_j}{\left|{\Delta Volume}_j\right|}}$$

(12)

Maps displaying ${Volume}_{diff}$ serve to contextualize the differences between $\Delta Volume$ produced by each case and to overcome the difficulty in conceptualizing large numbers (in the order of 10⁶) by reducing the differences to a simple percentage.

Of key interest to coastal modelers is the evolution of the subaerial island shorelines and the rate at which that change occurs. This was explored using two methods. First, a quantitative analysis was undertaken of rates of shoreline change on decadal timescales. The zero-meter contour is used as a proxy for the shoreline position, following the methods of [41]. Rates of shoreline change, or shoreline linear regression rates ($LRR$) were calculated for each case by tracking the position of the foreshore (open ocean) zero-crossing on DSAS transects [61], following the methods of [41]. Second a qualitative overview was undertaken of landform changes to the subaerial island and subaqueous island platform, as the end states vary between sensitivity cases.

3. Results

3.1. Bed Level and Volume Change

Figure 5A provides a color map of $\Delta DPS$ for the Benchmark case following a 10-year MorMerge simulation: blue indicates areas of deposition, and orange indicates areas of erosion.

$\Delta DPS$ for the benchmark case showed a general range of ±2 m, with the areas of greatest morphological activity confined mainly to the Ebb-Tidal Shoal outside the mouth of Mobile Bay and in Petit Bois Pass, attributable to tidal velocities [43] through these channels. There was also significant morphological activity on the eastern flank of the deep channel through Main Pass and around the western terminus of Fort Morgan peninsula, where (ebb) tidal velocities are also significant [43] ($\Delta DPS=\pm 4 \mathrm{m}$; Figure 5A). Further, there are patterns of deposition and erosion in the nearshore surf zone of both Petit Bois Island and Fort Morgan Peninsula with $\Delta DPS$ range of ±1 m. In the nearshore region of Dauphin Island, there is also a general alongshore alternating pattern of deposition and erosion ($\Delta DPS=\pm 1\mathrm{m}$) attributable to wave action and alongshore transport [50,52]. Deposition is greatest on the westward side of Pelican Island, the Dauphin Island shoreline immediately adjoining Pelican Island, on the foreshore of Katrina Cut, and the foreshore of Dauphin Island’s western terminus (though the very tip of the western terminus is eroded).

Differences in final bed levels (${∆Z}_f$) between the Benchmark and the Single-Layer, Four-Class case (SL4C) (Figure 5B) show that the region offshore of Dauphin Island from approximately the shoreline to the 10 m depth contour experiences greater erosion when multiple sediment classes are included (${∆Z}_f$ [−2 to −8 m]). This region of intensified erosion coincides with foreshore regions where EM1 has prevalence greater than 0% (Figure 3A). Figure 6B,C also indicate that the Ebb-Tidal Shoal and the eastern flank of the deep channel through Main Pass remain morphologically active in the sensitivity cases (${∆Z}_f=0 \mathrm{m}$ generally, at depths less than 5 m). There are notable exceptions on the Ebb-Tidal Shoal, however. Southeast of Pelican Island, within the 5 m depth contour, SL4C shows greater deposition (${∆Z}_f$ up to 1 m greater) than the same region in the Underlayer, Four-Class case (UL4C) and Benchmark. In Figure 5B, ${∆Z}_f$ further shows areas of greater deposition in Petit Bois Pass relative to the Benchmark case (${∆Z}_f$ [1 to 3 m]). Taken together, these areas of increased deposition in SL4C suggest that additional EMs, particularly mobile sediments such as EM1, serve as an additional sediment source for deposition in locations adjacent to their initial placement.

The ${∆Z}_f$ between UL4C and the Benchmark (Figure 5C) shows the influence of an underlayer bed configuration in addition to multiple sediment classes on modeled morphology. A broadly erosive region is also present in the UL4C case (Figure 5C) (${∆Z}_f$ [−1 to −4 m]), again falling between the shoreline and 10 m and coinciding with greater than 0% prevalence of the EM1 sediment class. Further, in the very nearshore of Dauphin Island’s west end, very little apparent erosion takes place in the UL4C case (${∆Z}_f$ of approximately 1 m). This is significantly reduced relative to SL4C. As was true of the SL4C case, ${∆Z}_f$ between the Benchmark and UL4C reveals regions of greater deposition in Petit Bois Pass and on the Ebb-Tidal Shoal (Figure 5C; ${∆Z}_f$ of 1 m relative to the Benchmark). ${∆Z}_f$ here is reduced in the UL4C case relative to the SL4C case: 1 m shoaled compared to 2 m shoaled. Overall, ${∆Z}_f$ shows relatively less erosion taking place in the UL4C compared to SL4C, though with areas of significant erosion coinciding with greater than 0% prevalence of the highly mobile EM1 sediment class (Figure 3A).

Figure 6A–C show the calculated $\Delta Volume$ of all blocks for the Benchmark, SL4C, and UL4C cases, respectively, for all volume change analysis blocks and regions established in Figure 4. Blue indicates blocks where the volume increased overall (deposition) while orange indicates a decrease in overall volume (erosion), following the 10-year MorMerge simulation.

Figure 6D–F show ${Volume}_{diff}$ calculated from $\Delta Volume$ between each case. Figure 6D,E show ${Volume}_{diff}$ for SL4C and UL4C relative to the benchmark, respectively. Figure 6F gives ${Volume}_{diff}$ of UL4C relative to SL4C. Blue shades indicate positive percentage difference, while orange shades are associated with negative percent difference.

As $\Delta Volume$ is derived from $\Delta DPS$, previously described areas of morphological activity can be identified in Figure 6A, with the largest magnitudes of $\Delta Volume$ within Petit Bois Pass (regions 3 and 4) ($\Delta Volume$ = +0.97 million m³ (Mm³); $\Delta Volume$ = −0.77 Mm³, erosive) to the west, Main Pass to the east (region 11) ($\Delta Volume$ = +0.67 Mm³; $\Delta Volume$ = −0.47 Mm³), and on the outskirts of the Ebb-Tidal Shoal (region 13) ($\Delta Volume=$ +0.62 Mm³; $\Delta Volume$ = −0.42 Mm³). For areas alongshore (regions 6 and 7, Figure 4) where $\Delta DPS$ showed alternating patterns of deposition and erosion (Figure 5A,B), positive $\Delta Volume$ indicates more deposition than erosion on balance, allowing for overall shoreline accretion (Figure 6A; regions 6 and 7). The farther-offshore blocks of regions 8, 9, and 10 (adjacent to the 5 m contour and less than 10 m depth) show negative to negligible $\Delta Volume$ in the Benchmark case. Within regions 9 and 10, there is deposition of roughly equal magnitude to the north of the eroded volume blocks (±0.19 Mm³).

Within region 5, there is one volume block of note between sensitivity cases, which covers a portion of Dauphin Island’s western end and shows a small magnitude of positive $\Delta Volume$ in the Benchmark and UL4C cases (approximately +0.05 Mm³) (Figure 6A,C), while showing a significant magnitude of erosion ($\Delta Volume$ = −0.78 Mm³) in the SL4C case (Figure 6B). This coincides spatially with nearshore availability of EM1 and significant $\Delta DPS$ also in the SL4C case. That this block shows positive $\Delta Volume$ in UL4C and Benchmark indicates a sensitivity to the bed configuration when multiple sediment classes are present and again suggests an armoring of highly mobile EM1 in the UL4C configuration.

For both SL4C and UL4C, there is significant negative (erosional) $\Delta Volume$ throughout region 6, with a greater magnitude of erosion occurring in the single-layer case. The adjoining blocks of region 6 and region 7 (immediate area of Katrina Cut; Figure 6A) show positive $\Delta Volume$ for the Benchmark case, while the same blocks show positive and negative $\Delta Volume$ in Figure 6B,C, respectively. The magnitude of deposition in the westernmost blocks of region 7 is also increased in SL4C (6B) and UL4C (6C) relative to the Benchmark (6A), with the greatest positive $\Delta Volume$ occurring in SL4C in region 7 (${Volume}_{diff}=+90\%$ relative to Benchmark, ${Volume}_{diff}=-50\%$ UL4C relative to SL4C; Figure 6D and Figure 6F, respectively).

The farther-offshore blocks of regions 8, 9, 10, and 13 (adjacent to the 5 m contour and less than 10 m depth) falling on the edge of the Ebb-Tidal Shoal present distinct and different behavior in SL4C and UL4C. From SL4C (Figure 6B), where the Benchmark showed negligible $\Delta Volume$, SL4C case shows a large magnitude of volume loss (${Volume}_{diff}$ ranging −90% to −100%; Figure 6D). Further, there is deposition in several blocks, which in the Benchmark shows erosion or negligible $\Delta Volume$, particularly in blocks just offshore of Pelican Island (${Volume}_{diff}$ = ranging +90 to +100%; Figure 6D). Within the same offshore Ebb-Tidal Shoal blocks in UL4C, Figure 6C shows a lower magnitude of $\Delta Volume$ compared to SL4C, though still volume loss (${Volume}_{diff}$ ranging +40 to +60%; Figure 6F), while adjacent blocks, which had shown a gain of volume in SL4C, experience volume loss in the UL4C case (${Volume}_{diff}$ = −80%; Figure 6F).

Several blocks within region 10 and 11 within Main Pass that experienced erosion under the Benchmark case experience depositional $\Delta Volume$ under both sensitivity cases. Blocks in regions 10 and 11 covering the Ebb-Tidal Shoal that experienced relatively minimal $\Delta Volume$ under the benchmark case experience significant positive $\Delta Volume$ in SL4C (${Volume}_{diff}$ ≥ 100%; Figure 6D), though with a mixture of positive and negative $\Delta Volume$ in UL4C (${Volume}_{diff}$ +100%, and −40% to −80%; Figure 6F).

Across each of the three cases, $\Delta Volume$ in regions 1, 2, 3, and 4 (Petit Bois Island and Petit Bois Pass) and in region 12 (east of Main Pass) all show identical patterns of deposition and erosion. The same blocks show deposition, erosion, or zero $\Delta Volume$ across all cases, though with varying magnitude, the largest magnitudes of $\Delta Volume$ occurring for SL4C (${Volume}_{diff}$ of −20 to −80% in region 1, −40% to −60% in region 2, ±40% in region 3, ±40% and ±80% in region 4, relative to Benchmark; Figure 6D). ${Volume}_{diff}$ of UL4C relative to Benchmark is identical for regions 1–3, but has greater volume loss in three blocks of region 4 with ${Volume}_{diff}$ of −10% to −20%, as shown in Figure 6F. Morphological change is relatively less at Petit Bois Island (regions 1 and 2; Figure 4) for the Benchmark case, but not insignificant ($\Delta DPS$ in the range of ±1 m; Figure 5A). Similarly, region 14 (over Fort Morgan Peninsula) also shows identical, apparently negligible, $\Delta Volume$ across the three cases. Each of these regions (1, 2, 3, 4, 12, and 14) fall over areas of the model domain where only EM2 is available at the start of the run, indicating that a large prevalence of EM may be a controlling factor on overall bed morphology change. Increases in deposition for SL4C and UL4C suggest downstream availability of the other additional sediments, reworked by tidal velocities in a way analogous to the Benchmark. Reduced deposition for UL4C in region 4 again points to a limiting effect brought on by the underlayer bed configuration. We note that in region 14, as shown in Figure 6, panels D and E, is significant, though little to no morphological activity was shown by either the $\Delta DPS$ or $\Delta Volume$, and so large magnitude in region 14 are the result of dividing by small numbers (see Equation (12)).

Figure 6D shows that ${Volume}_{diff}$ between the SL4C and the Benchmark is generally significant, reaching or exceeding differences of ±100% along the foreshore of Dauphin Island and within the majority of Main Pass. The same is generally true of ${Volume}_{diff}$ between the UL4C and the Benchmark (Figure 6E). The demonstrated sensitivity of volume change to both multiple sediment classes and to initial bed configuration in the nearshore and shoreface has definite implications regarding the controls these two parameterizations may have on open-ocean shoreline evolution.

3.2. Modeled and Historical Trends in Shoreline Change Rates

Smith et al. [51] calculated rates of shoreline change as weighted linear regression rates (LRRs) on well-established DSAS shoreline analysis transects [62] across Dauphin Island. Shoreline positions used to calculate historical rates of shoreline change were extracted from observational shoreline images (aerial imagery and lidar) at Dauphin Island for the 75-year period of 1940–2015. Smith et al. [51] additionally provided shoreline LRRs for subsets of the 75-year period. For the purposes of this study, two subset historical periods are treated as proxies for general trends in recent rates of shoreline change at Dauphin Island, AL, and serve as a point of comparison to modeled shoreline change rates and trends. A comparison was made between observational LRRs for two overlapping periods (1997–2015; 1985–2006) and modeled shoreline LRRs calculated following 10-year MorMerge simulations of the Benchmark and sensitivity cases. The periods of 1997–2015 and 1985–2006 were selected for two reasons. First, though the historical periods do not match the simulation period, the wave climatology applied to the simulations was developed from a record that spanned these historical periods, meaning that both the historical periods and modeled periods experienced similar wave forcing (with the exception of storms). Second, the model was initialized with a DEM dated 2015 and simulated 10 years, as such robust historical LRRs that matched the simulation period were not available at the time of publication, meaning that the periods selected are the nearest “modern” shoreline change rates available. All LRRs from each modeled case are provided, even at points alongshore where the historical LRRs were not available. Extraction of shoreline positions calculation of LRRs from model outputs are described in Section 2.6 of this paper and [41].

Figure 7A,B, show results of this analysis for the periods of 1997–2015 and 1985–2006, respectively, in which historical shoreline LRRs and confidence intervals are reproduced and plotted alongside modeled LRRs.

The historical LRR (dashed black line) [51] and confidence interval plotted as gray shaded areas establish for both periods a northward migration of the shoreline on the western half of Dauphin Island (5 to 12 km alongshore) at a rate of approximately 5 m per year (m/yr) to greater than 10 m/yr. Northward migration rates of the open-ocean shoreline are similar farther alongshore (14 to 20 km) decreasing in magnitude as eastward distance alongshore increases. In both historical periods, southward migration of the shoreline (negative LRR) is only present at points toward the eastern end of the island between 20 km and 25 km alongshore. Both periods also display a southward migration of the shoreline on transects 4 km to 5 km alongshore.

Initial comparison of the Benchmark case (red dotted line) to the observed indicates a general agreement in the alongshore trends of shoreline LRRs, though the Benchmark broadly underpredicts northward migration present in the historical LRRs (1 m/yr in the Benchmark where the historical periods show up to 5 m/yr; 5 to 10 m/yr in the benchmark where the historical periods show 10 to 12 m/yr) and often show southward migration (−1 to −2 m/yr) where the observed shows northward migration. While perfect agreement between the modeled and observed cases is not expected due to discrepancies in time period and storm forcing, the Benchmark case falls largely outside of the confidence interval for the historical LRRs. The Single-Layer, Four-Class case (dashed blue line) shows greater agreement with both periods of historical LRRs: magnitudes of northward migration begin to approach those of the observed at multiple points alongshore and fall within the gray-shaded confidence interval more often than the Benchmark case. At the same time, the Single-Layer, Four-Class case between 6 km and 8 km alongshore shows greater southward migration rates (2–3 m/yr) than the Benchmark (both predict the incorrect direction of shoreline migration in this region). The Underlayer, Four-Class case (solid black line) falls within the one σ confidence interval to a greater extent than either the Single-Layer, Four-Class or the Benchmark case, while still reproducing southward migration where the observed describes northward migration, as with the prior two cases. The Underlayer, Four-Class case produces LRRs in near agreement with the observed LRRs between 16 and 18 km alongshore.

Root mean square error (RMSE) (m/yr) was calculated between the observational historical LRRs and model-derived LRRs of the Benchmark and sensitivity cases (

Table 2. Values of root mean square error (RMSE) between observationally derived and model-derived rates of shoreline change for two overlapping historical periods (1997–2015 and 1985–2006) and three model simulations. Units of RMSE are meters per year (m/yr).

Model Case	RMSE (1997–2015) (m/yr)	RMSE (1985–2006) (m/yr)
Benchmark	3.31	3.34
SL4C	3.29	3.21
UL4C	3.21	3.09

). As model complexity was increased, RMSE decreased, regardless of the historical period examined. The largest value of RMSE was always found between the Benchmark case LRR and the observation LRR, while the smallest value of RMSE was always found for the Underlayer, Four-Class case. Values of RMSE calculated for the 1997–2015 period were reduced by 1% and 4% relative to Benchmark RMSE for SL4C and UL4C, respectively. Values of RMSE calculated for the 1985–2006 period were reduced by 4% and 8% relative to Benchmark RMSE for SL4C and UL4C, respectively. Sensitivity in LRRs and sensitivity of bed level and volume change (nearshore and at depth) in combination indicates broad sensitivity of overall morphological evolution to sediment parameterization and bed configuration. This overall sensitivity may lead to significant changes in model output morphology and behavior of the modeled region.

3.3. Landform Change

Figure 8 shows the initial digital elevation model (DEM) for Dauphin Island, used to initialize each sensitivity case. The initial model shoreline is plotted (dashed pink line) along with the final shoreline contours produced by each case after the 10-year MorMerge simulations (Benchmark: red dotted line; SL4C: blue dashed line; and UL4C: solid black line). For simplicity, the zero-meter depth contour is used as a proxy for shoreline position in all cases.

Pelican Island and the nearby Ebb-Tidal Shoal have recorded patterns of island change over long periods. Pelican Island has been observed to widen and migrate northward to merge with Dauphin Island [43] Sandy shoals and marginal bars are known to form on the Ebb-Tidal Shoal, which then become subaerial and migrate to merge with Pelican Island [43]. Of great interest is the formation of these new small subaerial islands and shoals, which are described in the historical recordkeeping and in observation and theory on ebb shoals and function as an important factor in the sediment balance in the nearshore region [63]. The Benchmark case, with only a single sediment class described, is unable to produce these ephemeral features known to occur in this barrier island system that are critical for accurate assessments of island resilience in the long term.

In 2005, Dauphin Island was breached by Hurricane Katrina [64] completely severing the eastern end of the island from the undeveloped western end. The breach, referred to as “Katrina Cut” (Figure 1) remained open until it was closed in 2011 by a constructed rock wall. The morphological activity around Katrina Cut has remained of interest to modelers and coastal managers as a matter of island resilience and the long-term feasibility of such engineered features and solutions. Figure 8B shows the area around Katrina Cut, and a cross section of the island is shown in Figure 8E. In each of the three sensitivity cases (Benchmark, SL4C, and UL4C), the shoreline migrates in the offshore direction, as sediment accumulates in the foreshore of the rock wall structure. This offshore migration is greater in both sensitivity cases than in the benchmark. While Figure 7 shows the northward migration of the shoreline at the regions immediately east and west of Katrina Cut, note that the area of the cut is excluded from the historical analysis, and that both historical periods (1997–2015 and 1985–2006) include shoreline change due to Hurricane Katrina. This northward migration is reproduced in the Benchmark, SL4C, and UL4C, though at reduced rates, which may reflect real-world accumulation of sand in the foreshore of the rock wall since the time of construction [65], where subaerial sandy beach has reformed with a southward migrating shoreline, excluded from the observational shoreline LRR. There is alongshore variability in the final position in both sensitivity cases, though both agree regarding the direction of shoreline migration. The cross section given in Figure 8E shows that there is variability in the slope of the nearshore, with UL4C presenting the steepest slope compared to either the Benchmark or SL4C.

The left column of Figure 8 (panels A, B, and C) shows areas of morphological interest already highlighted by results presented in Section 3.2 and Section 3.3. Beginning with what does not vary between the sensitivity cases, there is no evolution in any case along the back-barrier shoreline (coincident with $\Delta DPS=0 \mathrm{m}$ and ${∆Z}_f=0 \mathrm{m};$ Figure 5A–C). This lack of change on the back barrier shoreline is to be expected, as the probabilistic runs contain no storms or wind forcing needed to induce overtopping or to intensify wave activity in Mississippi Sound and redistribute back-barrier sediments [38,41].

The benchmark case (red dotted line, 8A) shows little change in shoreline position over the 10-year MorMerge simulation, with some erosion at the southern tip, and deposition on the subaqueous platform (Figure 8D).

Figure 8C shows a transect of the west end sand spit of Dauphin Island. There is somewhat more westward progradation of the shoreline in SL4C (dash-dotted blue line, Figure 8C,F), which is more pronounced as we follow the shoreline from west to east, compared to the Benchmark and UL4C (dotted red and solid black line, respectively, in Figure 8C,F). UL4C displays slightly increased westward movement limited to between −88.33 and −88.325 degrees west. More apparent in UL4C is a seaward migration of the shoreline and expansion of the subaerial island width, forming a sandy lobe. This reduced westward movement and greater shoreline accretion again point to an “armoring” effect in the UL4C case, reducing resuspension and therefore transport, leading to sediments to be deposited here.

4. Discussion

4.1. Decadal Morphology, Shoreline Position, and Landform Change

It is the case that in terms of morphology change, the SL4C and UL4C cases show that the model is sensitive to the inclusion of multiple sediment classes, spatial availability of sediment classes, and to initial bed configuration. Those differences have been described in terms of bed-level change, volume change, and trends in rates of shoreline change, all of which together influence the overall morphological evolution of the barrier island system. These sensitivities produce interesting changes to the resultant morphology and landforms, both subaerial and subaqueous, of great interest to coastal modelers. More investigation is required to determine whether SL4C or UL4C is closer to the true morphological evolution at Katrina Cut and the implications of this analysis for the design of feasibility and longevity investigations for similar engineered barrier island features.

The model is relatively more sensitive to the inclusion of multiple sediment classes in a single fully mixed layer than to the same classes applied in the underlayer bed configuration. When included in a single fully mixed layer, all sediments are available for transport as soon as sediment transport and morphology are computed. A 10 m-thick sediment layer composed significantly of highly mobile cohesive sediment means that much of that sediment moves rapidly, leading to significant erosion of the bed where sediment in that sediment class is present. The strong patterns of erosion, particularly in the “single-layer” case, again show the sensitivity of the spatial availability of sediments, particularly highly mobile fine sediments (EM1), even in limited quantities (Figure 3A). At the same time, the reduced erosion in those same regions shown by the “underlayer” tests in both $\Delta DPS$ and $\Delta Volume$ indicate that the inclusion of sand-composed underlayers and limiting cohesive sediments to the transport layer does lead to some “armoring” behavior as the mud class fluxes out of the transport layer and the transport layer is replenished by the underlayers with sand classes. It may be possible that this effect could be increased by reducing the thickness of the initial transport layer.

In addition to bed-change sensitivity, Figure 5B,C, indicate that the shoreline position over time will be sensitive to fine sediments placed close to shore as those fine sediments rapidly mobilize and erode in the surf zone where wave-driven shear stresses are greatest. Small-magnitude settling velocities lead to longer time in suspension and thus a greater likelihood that fine sediments will be transported away from the shoreline before redeposition or be subject to resuspension at a greater rate than larger sediments. This shoreline sensitivity to nearby fine sediments is further illustrated by Figure 7, where rates or shoreline are increased in SL4C when compared to UL4C and the Benchmark. Figure 5C highlights that UL4C is less erosive overall than SL4C; however, Figure 5B also shows that SL4C end state is shallower than UL4C end state to the south of Pelican Island and nearshore to Dauphin. Therefore, while 6 μm sediments were very quickly eroded, that suspended supply meant more potentially available to settle out nearshore, allowing for the shoal formation and shoreline progradation. This effect is dampened by the presence of underlayers containing only sand sediments in UL4C, which resupply the transport layer and reduce available EM1 over time.

Analysis of volume change within blocks informs the assessment of model sensitivity, reasserting and confirming what was shown by bed-change analysis. Block regions with shared patterns of deposition and erosion (shown by Figure 6A–C) suggest that morphology change at these locations (regions 1,2, 3, 4, and 12) is controlled primarily by forces local to those regions (tidal velocities and waves), as each case has the same initial sediment available in the same prevalence (EM2). However, the differing magnitudes between cases point to a secondary control on morphology change by readily available sediments. The differing magnitudes of deposition and erosion (SL4C greatest, UL4C lesser, Benchmark least) suggest a transport of the other available sediments into these active regions, leading to differences even outside of areas where EMs are initially available, i.e., a non-local, “downstream” sensitivity to multiple sediment classes. This may be expected in similar study regions with significant alongshore sediment transport to the westward transport at Dauphin Island [48]. ${Volume}_{diff}$ in blocks 1, 2, and 3 (Figure 6F) shows no difference between SL4C and UL4C, with region 4 showing ${Volume}_{diff}$ of −20% to +40%. This suggests that the influence of the “armoring” of UL4C (and sensitivity to the initial bed configuration) is more spatially limited than the influence of multiple sediment classes, only showing an impact in regions adjacent to those where multiple classes are initially available. Overall, in both SL4C and UL4C, the largest ${Volume}_{diff}$ values are associated with areas where EM1 is available at model initialization (negative ${Volume}_{diff}$, regions 5, 6, 7, 8, 9 and 10; Figure 6D) or where mobile sediments are deposited (Main Pass shoaling) (positive ${Volume}_{diff}$, Figure 6D).

In examining the final island footprint of the 10-year probabilistic sensitivity runs (Figure 8), both cases (SL4C, and UL4C) capture island evolution behavior that is realistic compared to known historical behavior of this island system and regions like it, which were not captured by the probabilistic P-J case. Also indicated by the shoreline analysis (LRRs and final shoreline contours (Figure 8)) is that some morphological behaviors typical of certain barrier island systems and to Dauphin Island (i.e., the formation of sand shoals and the seaward formation of a lobe on the West End Spit) are influenced by the presence of additional sediment classes, and do not occur when only a single sand class is used (Figure 8). Further, the shoreline linear regression rates of change indicate that the shoreline is overall much more mobile than the Benchmark in the sensitivity cases. While model-derived Benchmark LRRs generally fall outside the one σ uncertainty band for the observed LRRs, this can partially be accounted for by the fact that the P-J model was designed to model long-term quiescent morphology change, without the pulses of island change caused by storms. Even so, the increasing subsequent complexity of SL4C to UL4C did move trend lines further towards and within the uncertainty band.

The results of this work are transferable to other similar coastal environments. It is important to note that this study was not meant to serve as a calibration to improve model accuracy compared to observations, though the increased complexity of the subsequent sensitivity cases did reduce RMSE against historical shoreline LRR trends. Rather, the aim was to improve the understanding of how applying different sediment classes (types and grain sizes) and initial bed-layer compositions affects the modeled sediment transport patterns and resulting bed levels. We did not consider how the sediment parameterizations may influence sediment transport under conditions of wind-driven currents or storm events, which may show additional sensitivities, particularly along back-barrier shorelines.

Though quiescent, fair-weather conditions may reflect the long-term average conditions of most barrier island systems, they do not reflect the full gamut of anthropogenic, atmospheric, and hydrodynamic forcing. Pulses to these systems come in the form of storms, which have historically reworked shorelines and elevations through storm waves, overwash [34], and surge, and in the form of anthropogenic engineered restoration in response to storm impacts. Changes to storm surge hydrodynamics caused by sea-level rise [36] further complicate barrier island forecasting, as do anthropogenic alterations [52] to barrier islands, which can increase vulnerability and island drowning. Further, continuous or discontinuous retreat of barriers is modulated by sea-level rise rates and storm intensity [34], with storm pulses to these long-standing barrier systems expected to increase in frequency over time. As recovery of these systems after storms and management actions take place on the timescale of decades [53], numerical modeling of barrier islands needs a considered approach to capture both the quiescent aspect and storm pulses. Multi-model model train frameworks [38,39], which aim to efficiently forecast fully forced barrier island evolution over multi-decadal timescales, can benefit from improvements to quiescent period morphology results such as those shown by this study, gained by including multiple sediment classes arranged in an armoring underlayer bed configuration.

The P-J Benchmark model itself was purpose-built to forecast and hindcast quiescent period decadal shoreline and seafloor change for feasibility studies of engineered restoration features in a model train. The morphological models developed as sensitivity tests for this study have wide-ranging applicability in large-scale coastal and estuarine modeling. The combination of climatological wave forcing, morphological tide forcing, probabilistic MorMerge morphological evolution, and the modeling/tracking of multiple sediment classes allows robust modeling of decadal-scale morphology changes in subaerial features, subaqueous features, and barrier island shorelines. Models such as the SL4C and UL4C setups can be further applied to habitat/environmental impact studies. For example, one such model was applied to a marsh–estuary system to study the response to restoration features [5] by simultaneously modeling morphology impacts, feature longevity, and impacts to suspended sediment concentrations. The models applied in this study also have potential applicability as companion models for sediment budget research in open-coast, estuarine, or tidal inlet settings by approaching sediment budgets via bed volume change and sediment class fraction at the bed [66,67].

4.2. Controlling Processes: Short-Term, Deterministic P-J Model

An additional limited analysis was carried out to investigate the sediment transport processes controlling the sensitivity seen in decadal morphology change of the seabed, shoreline, and island platform. The P-J model Benchmark, SL4C, and UL4C were each run in a deterministic mode from 19 June 2015 to 19 July 2015, recreating a portion of a 6-month validation period previously published [41]. Morphology and sediment settings were set in the deterministic mode, as they were set in the probabilistic MorMerge simulations of the Benchmark, SL4C, and UL4C cases, except for morphological acceleration and MorMerge settings, which were turned off. For full details of deterministic simulation model settings, see Jenkins et al. [41]. Sediment transport values were output at a set of six observational cross sections in the model domain during the one-month simulations. The locations of the model observation output transects are provided in Figure 9 in relation to the island shoreline and to the prevalence of sediment classes.

Final cumulative sediment fluxes (${CSF}_{final}$) in meters cubed (${\mathrm{m}}^3$) were output on six observational cross sections within the model domain for each sediment class individually, as well as for the combined bed-load (${CSF}_{final}^{bed}$) (excluding the 6 µm sediment class), combined suspended load (${CSF}_{final}^{susp.}$), the combined suspended load of sand classes (${CSF}_{final}^{sand}$), and finally the sum of ${CSF}_{final}^{bed}$ and ${CSF}_{final}^{susp.}$, referred to as ${CSF}_{final}^{total}$, which is the total sediment flux irrespective of type through a given transect over the course of a simulation. For the purposes of this study, we considered two cross sections placed over two regions of significant morphological activity and relatively greater sediment transport: Main Pass (MP), and Pelican Offshore (PO). For full details of output, ${CSF}_{final}$, see the associated Supplementary Materials (Table S1).

Looking first at the Benchmark deterministic P-J model case, a baseline of sediment fluxes is established. Through MP, a ${CSF}_{final}^{total}$ of 23,626 m³ of sand was transported southward through the inlet at a ratio of 89% suspended-load to 11% bed-load flux. These fluxes are in keeping with the known characterization of this inlet as tidally dominant [43], with transport suggesting dominance of ebb tidal velocities. PO is the second-most active, with a ${CSF}_{final}^{total}$ of 1498 m³ of sand transported northward at a ratio of 95% to 5% suspended to bed-load transport, The northward flux direction is in keeping with a well-established historical pattern of Pelican Island elongating southward over time as northward transporting sediments accumulate, while simultaneously translating northward, often interfacing with the shoreline of Dauphin Island’s east end and supplying sediment to Dauphin Island. Additionally, these fluxes likely contribute to the formation of subaqueous shoals and marginal bars, which historically develop into disconnected subaerial sand spits offshore and south of Pelican Island [43].

Where EM1 is present, even in relatively small percentages (20%, immediate foreshore of Dauphin Island; Figure 3: 0% initially at PO; Figure 9A), in the SL4C and UL4C cases, the suspended-load flux signal is dominated by the suspension and transport of fine-silt sediments. Through the PO cross section, the suspended-load transport of EM1 makes up greater than 70% of the total transport in both SL4C and UL4C. This is an expected outcome when including highly mobile sediment classes with relatively small settling velocities. Large negative values of ${∆Z}_f$, $\Delta Volume$, and ${Volume}_{diff}$, which are noted seaward from Dauphin Island and the Pelican Island shoreline, can be accounted for by suspension of EM1 and subsequent transport away before it can be reworked into the nearshore bed or sequestered into underlayers. Negative values of ${CSF}_{final}$ (−2464 m³) for EM1 through the PO cross section confirm this offshore flux of fine sediments. The offshore transport of EM1 is in keeping with transport patterns for fine sediments established by Hummel [43], who described clay and silt sediments as being transported seaward and deposited mainly on the shelf. The seaward flux of EM1 in SL4C and UL4C is supported by the results shown in Figure 5B,C, with erosion in regions where EM1 was initially present, but without equal measures of local or adjacent deposition. The direction of flux for EM1 compared to the sand classes also highlights that inclusion of various sediment classes and types influences the behavior of the system. In the 1-month deterministic case through this transect, suspended fine sediment moves offshore while suspended sand moves onshore, which may be significant for a given system and may lead to non-local, “downstream” impacts on sediment transport and morphology change. When a probabilistic (MorMerge) approach to morphological evolution is required of a modeling application, fine cohesive sediments such as EM1 may not be a good fit for a morphologically active surf-zone application. The rapid erosion and movement offshore conflicts with the periodicity of forcing and transport processes that Delft3D MorMerge is best suited to. In cases where such fine sediments must be included to capture transport or total suspended sediment concentrations while morphology change is also of interest, then the results of this study point to the use of an armoring underlayer bed configuration with mud classes limited to a thin, superficial transport layer.

The 200 µm non-cohesive sediment class (EM2) is present in each of the three sensitivity cases (Benchmark, SL4C, and UL4C); however, differences in the fluxes of EM2 between cases SL4C and UL4C indicate sensitivity to bed configuration. This again points to an armoring effect in the UL4C case. Through the MP cross section, ${CSF}_{final}^{total}$ varies between sensitivity cases by ±1%. EM3 and EM4 show reductions in both suspended-load and bed-load fluxes in the range of 7% to 10%. This variability in fluxes through MP highlights the sensitivity of individual sediment classes to bed configuration. Relative abundance of non-cohesive EMs (EM2, EM3, and EM4) is significant at MP and relatively greater than EM1. Differences in fluxes of non-cohesive sediments between SL4C and UL4C help illustrate the mechanisms at work between sensitivity cases. At each model time step, the mass-balance and transport equations are solved for each EM individually. Suspension and transport of each non-cohesive EM (EM2, EM3, and EM4) occur in proportion to the proportion of each EM at a given point. SL4C has a single uniformly well-mixed bed layer, and so the proportion of fluxed EMs is dictated by the initial relative abundance of EMs and the updated relative abundance as sediments are deposited. However, when a transport layer and bookkeeping underlayers are introduced, there is an additional pathway by which the relative makeup of the bed can be altered. If, after a given time step, there is a net outward flux of the transport layer (thus reducing transport layer thickness), the transport layer is replenished by the underlayers to maintain a constant thickness, and that resupply to the transport layer is proportionate to the makeup of the underlayers. For instance, EM2 is the most abundant of the EMs present at MP, and resupply will be made up of a greater fraction of EM2. These differences in fluxes scaled to the 10-year simulation and multiplied by a morphological acceleration factor help to explain the damping effect that UL4C has on deposition and erosion relative to SL4C.

5. Conclusions

This study assessed the sensitivity of sediment transport and resulting morphological change to input sediment parameters of sediment classes and initial bed-layer configuration. A previously developed and validated D3D model was used to simulate long-term (decadal) tides, waves and sediment transport, and morphology in the Dauphin Island, AL, USA, region. Using collected field data and an end-member analysis for grain sizes, a method for generating model input sediment grain sizes and their spatial distributions was presented. Two test cases of progressive complexity were developed to assess model sensitivity to sediment inputs: multiple spatially varying sediment classes and initial bed configuration. The first case used multiple sediment classes (cohesive and non-cohesive of varying size) in a single fully mixed bed layer. The second case used the same sediment classes with a thin transport layer and bookkeeping underlayers, which initially excluded cohesive sediments. A previously developed and validated model was used as the benchmark against which sensitivity was evaluated. The use of the P-J model for Dauphin Island, combined with in situ sediment data-derived model input, allowed for a novel holistic analysis that provides narrow, quantitative sensitivity and broad, qualitative assessment of system benefits and effects produced by common modeling practices when ample sediment data are available. In other words, this study provides insights into what is gained by the inclusion of more realistic sediment input data and model shortcomings when not using multiple sediment fractions or when opting for a simpler bed configuration.

The results indicated that the model was sensitive to the sediment parameterizations, particularly when simulating erosion/deposition volumes and rates of shoreline change. Shoreline analysis (LRRs and final shoreline contours) showed that morphological end states and phenomena typical of certain barrier island systems (e.g., the formation/migration of sand shoals on ebb shoals and the seaward formation of a lobe on the West End Spit) are influenced by the presence of additional sediment classes. The model is incapable of producing these morphological features in the same 10-year simulation period when only a single sand class is used. Although storm pulses to the barrier system needed to drive island rollover and retreat are absent from the P-J model (shoreline change due to storms is included in the historical record), inclusion of multiple sediment classes leads to improvements in shoreline LRRs toward historical trends.

There are sensitivities to the inclusion of multiple sediment classes, both local and “downstream,” especially highly mobile fine sediments, even with limited spatial availability and limited relative abundance. Sensitivity to bed-layer configuration is also significant, illustrated by analysis of bed-level change (${∆Z}_f$ and $\Delta Volume$) with reduced erosion in active regions and reduced deposition in areas downstream of the direction of net sediment transport. The “armoring” effect of the underlayer bed configuration is locally limited to areas where multiple EMs are initially present and immediately adjacent to those regions.

Based on the findings of this study, the following could be beneficial for future modeling studies interested in simulating morphology with cohesive and non-cohesive sediment transport in coastal barrier island environments.

• Applying a single sediment class may be permissible for studies interested in describing general or net sediment transport patterns and resultant morphology thereof.
• Including multiple sediment classes may be useful for resolving long-term shoreline change, particularly along features like spits and shoals, and to examine their influence on suspended transport of cohesive and non-cohesive sediments, especially within inlets and on ebb tidal deltas.
• Applying an underlayer configuration versus a single fully mixed layer can provide an “armoring” affect over long simulation times, localized to regions where multiple sediment classes are present initially.
• Only if fine, cohesive sediments are included in a superficial, thin transport layer in an underlayer bed configuration that limits the activity of non-periodic sediment transport processes through armoring will the results be accurate and consistent.