Development and Application of an Empirical Dune Growth Model for Evaluating Barrier Island Recovery from Storms

.


Introduction
Numerous approaches exist to estimate the response of dunes to extreme events, including storm-impact scales [1][2][3], relatively simple 1D models [4,5], probabilistic approaches [6][7][8][9], and complex 2D numerical models that can resolve longshore variability in dune erosion, overwash, and breaching [10].Approaches to describe the recovery and rebuilding of dunes following major disturbances have also been developed [11,12].The timescale of dune recovery, which is long compared to storm-response timescale, is particularly relevant to decisions regarding restoration and management of ecologically and economically valuable coastal systems [13].Techniques are required for predicting dune formation and growth that can be used in conjunction with storm-response models to understand the evolution of the beach-dune system on decadal timescales.
Coastal dunes build when Aeolian transport carries sediment from the subaerial beach landward or seaward until a sufficient perturbation (ripples, vegetation, sand fencing, etc.) disturbs the flow.A positive feedback develops, resulting in sediment deposition and the formation of dunes that may vary spatially and temporally in height, width, and form [14][15][16][17][18]. The morphology of dunes is a key factor in the response of barrier islands to storms and, ultimately, their long-term evolution [1,19,20].Several formulations have been established to relate sediment flux via Aeolian transport to wind speed [21,22].Intermittent transport begins when a critical threshold of wind speed is reached, then increases with wind speed until continuous transport (but variable mass flux) occurs [23].Sediment flux in beach environments rarely reaches the maximum potential transport predicted by these formulations [24][25][26].For example, flux is reduced when the fetch, i.e., the length over which the wind blows, is below a critical distance [24,25,27].In the coastal zone, fetch length is controlled by the width of the beach and the wind angle [14,28].Fetch length at a given site varies with time due to variation in beach width and shoreline orientation driven by short-and long-term processes of coastal erosion and accretion (tides, storms, sea-level rise, etc.) [29][30][31][32].Aeolian mass flux also depends on sediment characteristics, such as particle size and density, and bed characteristics such as sorting or the presence of surface crusts that inhibit suspension [25].In addition, topographic features (berms, dunes, etc.) alter the local wind field and sediment entrainment [14,25,33].Flux is also reduced below theoretical maximum levels if the sediment is moist, which increases the threshold wind speed necessary to induce transport [27,34].Many of the factors controlling Aeolian transport are highly variable in space and time, are difficult to measure in nearshore environments, particularly on regional scales (e.g., wetting and drying, sediment grain size, and wind gustiness), or both [23,25].For this reason, probabilistic approaches to estimate Aeolian sediment flux have been considered as an alternative to strictly deterministic methods [23,35].Integrating models of Aeolian sediment flux up to the formation and building of coastal dunes on the landform scale remains a challenge [36], and attempts to do so on the landscape scale have occurred only recently and with limited testing in real-world environments [37,38].
The complexity of controlling factors has led to conceptual models being developed to describe the formation and evolution of coastal dunes without explicitly resolving the underlying process of Aeolian transport.Coastal dune formation following destruction by storms has been categorized into four states: (1) rapid sediment deposition at the foreshore, (2) accretion of the back beach, (3) incipient (also known as embryonic) dune formation, and (4) dune and vegetation establishment [39,40].One of the earliest frameworks used to evaluate the terminal state of dune systems is the surf-beach-dune model, which extended the Wright and Short beach state model [41,42] to the evolution and form of coastal dunes [18,43].Dissipative beaches, characterized by wide, shallow-sloping beaches resulting from low-energy wave climates, are associated with higher and more persistent dune features.In contrast, reflective beaches are associated with higher energy environments, narrower beaches, and lower/less persistent coastal dunes.Recent testing with computational fluid dynamics (CFD) modelling has supported this conceptual model [33].However, there are many challenges to extending these conceptual models into numerical frameworks that are broadly applicable at the landscape scale [12].For example, the rebuilding of dunes following a storm event relies first on the recovery of the lower and upper beach [44], and robust and accurate quantitative models for capturing this process are limited [45].
Multiple additional factors have been identified as driving the location and development of dunes [46], which also confounds the development of predictive models.A strong link exists between dune morphodynamics and vegetation, the evolution of which is complex [14,37,47].Vegetation type and density vary considerably in space and time, driven by prevailing environmental conditions (rainfall, temperature, etc.) and feedback with dune evolution [48,49].Sediment supply in part determines dominant vegetation species and thus affects bio-physical feedbacks [50].Increased frictional roughness over vegetation results in sediment deposition and incipient dune formation, which can also occur along wrack and other perturbations on the beach [37,51].The sediment trapping characteristics of vegetation control the rate of sediment deposition, while the burial and saltwater sensitivity of individual species controls vegetation survival.As a result of these feedbacks, some plant species are associated with higher, narrow dunes, while others are associated with shorter and wider landforms [37,49,52].Because of the link between dunes and vegetation and the influence of sediment moisture on Aeolian transport, dune growth may exhibit seasonal variability [53].
As an alternative to process-based methods of predicting dune growth, Houser et al. [11] proposed an empirical approach for estimating coastal dune crest elevation increase based on a formulation developed for modeling dune building in interior (non-coastal) plains [54].In this approach, the maximum elevation of the foredune is predicted using a sigmoid growth curve following growth rate patterns observed for dune-building vegetation.Houser et al. found that parameterized sigmoid growth curves were capable of capturing dune elevation change for two barrier islands in the northern Gulf of Mexico (Santa Rosa Island, Florida and Galveston Island, Texas) [11].
An empirical dune growth model (EDGR) has been developed to predict the evolution of barrier islands foredunes based on the sigmoid growth tendency of dunes.The model was calibrated using light detection and ranging (lidar) data from Dauphin Island, Alabama, a barrier island in the northern Gulf of Mexico that has seen extensive storm damage and recovery in previous decades (Figure 1).In the previous work, EDGR was coupled with the XBeach storm impact model [10] and the Delft3D hydrodynamic and morphodynamic model [55] to hindcast [56] and forecast [57,58] the decadal scale evolution of the island.The coupled model framework was shown to have good skill in evolving characteristics of the island; however, EDGR was not run or assessed as a standalone model in this prior work.In addition, the relative influence of different processes driving dune recovery were not explored.In the current work, the EDGR model is evaluated under multiple model configurations and used to explore processes driving dune evolution at Dauphin Island.The mechanics of the model are included in the next section, along with a description of the configurations used in the current work.The results section assesses model skill under these varying configurations through comparison to observational data collected between 2005 and 2015.The discussion section illustrates how the model can be used to investigate spatial and temporal variability in dune growth, including evaluating the relative contribution of erosion and accretion processes.This section also considers the scope of model applicability and discusses its underlying assumptions.

Empirical Dune Growth (EDGR) Model
The EDGR dune model is briefly reviewed here with a complete description available in Table A2 of the work in [56].The foundation of EDGR is a novel mechanism for decomposing beach or barrier island elevation profiles into the sum of a set of shape functions representing landform features, referred to hereafter as the parameterized island Gaussian fit (PIGF) method.Each crossprofile (ZE(x,y,t)) is represented as the sum of Gaussian curves (ZG,n(x,y,t1)), where "n" enumerates each curve) representing the underlying island platform and, if present, one or more berm and/or dune features as (Figure 2 The fit deviation (error) is captured in ε(x,y).Each Gaussian is defined by height (Ghigh,n(x,t1)), full width at half maximum (Gwide,n(x,t1)), and location (GY,n(x,t1)) as Dune height is predicted following an equation presented in Houser et al. [11].The total maximum elevation of the foredune (Dhigh(x,t)), equivalent to the height of the foredune Gaussian curve (Ghigh(x,t)) plus the height of the underlying island platform, is calculated from the longshorevariable initial dune elevation (Dhigh,0), terminal dune elevation (Dhigh,F), growth rate (r), and time of incipient dune formation (t0) as A parameterized linear model is used to calculate the foredune Gaussian width (Gwide(x,t)) from the foredune Gaussian height (Ghigh(x,t)) using prescribed values for slope (mdWdH) and y-intercept (bdWdH) as (Figure 2) After the foredune is evolved, EDGR reconstructs the cross-shore profiles by replacing the foredune in the island profile from the prior time step with the predicted foredune curve (Figure 2).In addition to allowing foredune growth to be modeled, the PIGF method allows characteristics of individual dunes to be rapidly identified and quantified [59].EDGR includes formulations for incipient dune formation, with dunes forming either at a prescribed location (GY(x)) for use in hindcasting or at a location based on the distance (DSL,0(s)) to a post-storm shoreline (SLY(x)) to allow for barrier island retreat when forecasting: The model can be run standalone or in conjunction with an external model that predicts storm impacts or other changes to the beach and dune not associated with foredune growth.
Model parameters (Dhigh,F, GY, r, mdWdH, bdWdH) can be derived in multiple ways depending on the application of the model and available data.In cases where EDGR is used to fill gaps between elevation surveys, Gaussian basis set decomposition of the initial and final elevation data can be used.Otherwise, template (proxy) data from other time periods or sites can be used under the assumption that growth will be similar for the modelled dunes.When template data are used and the final configuration of the dunes is unknown, the longshore distribution of Dhigh,F for the template data is decomposed via fast Fourier transform (FFT) [60], the phases of the wavenumbers are randomized, and an inverse FFT is applied to reconstruct a new longshore distribution of Dhigh,F.This approach enables sensitivity testing of model predictions by allowing multiple scenarios of longshore varying dune height to be modeled, each of which retains the statistical characteristics (mean, standard deviation, and spatial scales of variability) of the template data.

Study Site and Model Configuration
The focus of the current study was on the western half of Dauphin Island, an uninhabited region that is separated from the community of Dauphin Island by Katrina Cut.This breach originally formed during Hurricanes Ivan (2004) and Katrina (2005).Although it was subsequently closed with a rock wall and has begun to infill with sediment, it provides a natural barrier between areas to the west of the Cut that evolve predominantly under natural processes and areas to the east that are heavily driven by anthropogenic influence.The island has been frequently inundated or overwashed during storm events over the last decades, resulting in multiple cycles of dune erosion and loss followed by recovery and growth.Additional information on the geophysical characteristics of Dauphin Island can be found in [61][62][63][64][65].
The island was divided into four analysis regions based on observed variability in dune evolution over the study period: R1, the wide western terminus of the island; R2, a narrower portion of the island associated with overwash fans in 2005 that approached or reached the bay shoreline; R3, a narrow portion of the island that incurred numerous small breaches during Katrina; and R4, the immediately adjacent to the region of Katrina Cut (Figure 3).Model parameters for Dauphin Island were previously derived from two lidar elevation surveys taken in 2005 and 2015 for each region as well as the complete model domain (Table 2 [56]).In the current study, the sensitivity of model predictions to these parameters is explored.The model framework was applied to predict the recovery of the dunes at Dauphin Island following Hurricane Katrina over the period of August 2005 to January 2015.Data for model assessment were taken from 12 topographic lidar surveys collected over this period (Figure 4).Additional information on the processing of these data and citations for data access may be found in [56,66].For the EDGR model application, data were interpolated to cross-shore transects with 5 m longshore resolution and 2.5 m cross-shore resolution.There were seven tropical storm events and hurricanes that passed within 200 km of the study area during this period (Figure 4; https://coast.noaa.gov/hurricanes/),resulting in dune erosion and occasional overwash and inundation of lower-lying areas of the island.Additional information on historical storms and their impacts at Dauphin Island can be found in [62,[67][68][69][70]. EDGR was applied at Dauphin Island to predict dune growth from 2005 to 2015 using six model configurations (Table 3) that evaluated the model's skill and capacity for analyzing dune evolution, including testing the sensitivity of the model to the quantity of input data and applicability in predicting dune growth where local, historical dune growth data are unavailable.In the first configuration, denoted E1, EDGR operated standalone with Dhigh,F and GY derived from 2015 data.This configuration benchmarks the mode skill in predicting dune width, cross-shore dune profile shape, and dune height between initial and terminal values.Three standalone experiments were run using FFT randomization to generate Dhigh,F and DSL,0 parameters from the 2015 lidar using different sections of the island as the terminal dune template: the entire site, denoted E2; data from region R3, denoted E3; and data from region R4, denoted E4.These configurations evaluated the skill of the model when parameterized with nonlocal data, such as in cases where historical data are unavailable for the site of interest.An additional experiment, denoted E5, was conducted using the observed dune locations and the FFT-randomized dune height values from configuration E2; this configuration allowed the error associated with a randomized terminal dune height to be isolated from errors in dune position.In the final configuration, denoted E6, EDGR was coupled to the lidar observations as a proxy erosion model to identify if elevation loss of the dune crest exceeded a defined threshold to restart dune growth.A threshold of 0.25 m was used, which was estimated to be greater than the error in the bias-corrected lidar [66].In this case, Dhigh,0 was reset as the elevation of the lidar dune crest with no other data assimilation.Observed dune position and FFT-generated dune heights from experiment E5 were used so that comparison between E5 and E6 quantifies the model improvement if dune erosion is accounted for.This configuration evaluates and illustrates the use of EDGR in analyzing the relative contribution of erosion and accretion processes to dune evolution.L2005 is a "no-change" model taken as the 2005 lidar, and the other six models are different configurations of EDGR.Dhigh,F is terminal dune height and DSL,0 is distance to shore.In all experiments, the growth rate ® is fixed at 0.53, and the slope and y-intercept of the linear model used to estimate foredune Gaussian width from height (mdWdH and bdWdH) are fixed at 9.5 and 13.9.

Experiment
ID.The 0.5 m contour in the 2005 lidar survey was taken as the shoreline location in calculating DSL,0 so that profiles with incomplete coverage of the intertidal zone could be included.For experiments using FFT randomization (E2-E6), 100 different sets of terminal dune characteristics were generated and run for each template configuration and model performance is presented as the mean and standard deviation of error statistics calculated individually for each of the 100 initializations.As a benchmark of EDGR's skill, an alternate "no-change" model (the 2005 lidar, hereafter referred to as L2005) was also evaluated.The L2005 model assessment represents the error that would result from using increasingly outdated dune topography.Dune height and position were extracted from L2005 as the point of maximum island elevation.EDGR time steps were chosen such that the model generated output at the times of the observed lidar survey (Figure 4, time step of several months to upwards of two years).The longshore and cross-shore resolution of the model matched the gridded lidar (2.5 m cross-shore resolution and 5 m longshore resolution).

Prescribed Terminal Dunes (Experiment E1)
EDGR formulations for identification of ZR(x,y,t1) and growth of the foredune as a Gaussian form were robust in application at Dauphin Island (Figure 5).Across R1-R4, the Gaussian fitting method for isolating the residual island profile in the 2005 initial DEM converged for 1912 out of 1918 of the cross-shore profiles (99.7%).A portion of these profiles were subaqueous in 2005 and/or a foredune was not identifiable in the 2015 lidar survey.The Gaussian foredune building method was successfully applied to all other profiles (100%, 94%, 62%, and 96% of profiles in R1-R4).The largest errors in the predicted DEM were seaward and shoreward of the EDGR-evolved foredune, i.e., in areas where the initial elevations were not updated by the model (Figure 5).ZE(x,y,t) was compared to the corresponding lidar (ZL(x,y,t)) over the area (Zdune(x,y,t)) that EDGR updated elevation (|ZE(x,y,t)−ZL(x,y,t0)| > 0.05 m) (Figure 6).EDGR's RMSE for Zdune increased slightly over the first 1-4 years post-storm before leveling out (0.31-0.50 m in 2015; lowest RMSE in R1).Bias was relatively consistent over time (−0.35 to 0.14 m).R2 decreased somewhat over the first 1-4 years in R2-R4 before leveling out to 0.51-0.61,while remaining consistently high for R1 (0.66-0.77).These predictions were benchmarked against L2005 over the same cross-shore area; the no-change model prediction skill initially deteriorated before leveling off in 2012-2015.3).Dhigh (left column) was taken as the maximum dune elevation for each cross-shore profile.The same statistics are shown for the 2005 lidar.Statistics for Zdune(x,y,t) (right column) are calculated at all cross-and longshore locations where EDGR updated the elevation for a given lidar survey.

Template Terminal Dunes (Experiments E2-E6)
Using R1-R4 as the EDGR template (experiment E2) improved predictions of DSL,0 over L2005 in R1-R3 while performing similarly to L2005 in R4 (Figure 7; Tables A1 and A2).EDGR had a negative bias in ΔSL,0 for the western regions (−15 ± 10 m and −8 ± 8 m for R1 and R2) and a slight positive bias for the eastern regions (3 ± 8 m and 7 ± 7 m for R3 and R4).RMSE in ΔSL,0 was similar throughout the study site, decreasing slightly from east to west (R1, 36 ± 6 m; R2, 27 ± 4 m, 27 ± 3 m, and 25 ± 3 m).R2 exhibited no spatial pattern and ranged from 0.06 to 0.18, with low values reflecting the high amount of scatter associated with randomizing terminal dune position.Spatial variation in the post-Katrina shoreline with a length scale greater than the 100 m Hanning filter is evident in the error in ΔSL,0 (Figure 8).Both EDGR and L2015 have spikes in error at locations where there are small discontinuities in the foredune, such as where overwash or breaching resulted in the dune reforming farther inland than the primary foredune at adjacent locations.3) referenced as distance to the post-Katrina shoreline (ε,ΔSL,0) for 100 independently initialized FFT randomizations (colored lines represent the mean, gray represents the standard deviation).The same statistics are shown for the 2005 lidar, which was used as the initial elevation in EDGR.
Experiment E2 had similar bias to L2005 for Dhigh in the first 1-2 years following Hurricane Katrina, after which time EDGR began to outperform L2005 (Figures 9; Tables A1 and A2).EDGR bias was positive for Dhigh in R1 (0.16 ± 0.26 m in Jan 2015) and R3 (0.44 ± 0.23 m in Jan 2015) and negative in R2 (−0.30 ± 0.24 m in Jan 2015) and R4 (−0.18 ± 0.11 m in Jan 2015).RMSE was similar for regions R1, R2, and R4 (0.69 ± 0.10 m, 0.74±0.12m, and 0.74 ± 0.10 m in Jan 2015) and lower than for L2005 after the first 2-4 years post-Katrina.RMSE in EDGR was somewhat higher for R3 (0.93 ± 0.16 m in Jan 2015) and was of a similar magnitude to RMSE for L2005 until 7-8 years post-Katrina.R2 for Dhigh in EDGR was similar across all regions (mean of 0.03-0.05for R1-R4 in 2015) and lower than for E1 and L2005.
Temporal patterns in model performance in experiment E3 (data from R3 as the dune template) and experiment E4 (data from R4 as the dune template) were similar to experiment E2 (full island used as the dune template), but use of a more spatially confined template with less longshore variability reduced the range of variability in predictions and error (Figures 7 and 9).Experiment E3 had similar error statistics to E2 in ΔSL,0 (mean RMSE varied by less than 3 m; mean bias varied by less than 3 m; and mean R2 varied by less than 0.03).A larger effect was seen in experiment E4 (mean RMSE varied by up to 5 m and mean bias became more negative by up to 10 m).For Dhigh, use of R3 as the template (experiment E3) increased RMSE in R1, R2, and R4 by up to 0.22 m; created a more negative bias; and reduced R2.As expected, RMSE and magnitude of bias in R3 were reduced slightly (by 0.07 and 0.10 m).RMSE for configuration E4 decreased in R1, R2, and R4 and increased in R3, resulting in an RMSE for this region of ~0.25 m more than the rest of the study site.Bias became more positive compared to experiment E2 (up to 0.31 m in R3), while R2 increased somewhat for each region.The most notable effect of using prescribed dune locations (experiments E5 vs. E2) on Dhigh was that bias became slightly more positive in all regions (change of 4-6 cm).The coupled configuration (E6) allowed evaluation of dune crest erosion through time.Elevation loss (ΔDhigh) throughout the study domain occurred more frequently prior to 2010 (Figure 4).The most widespread losses occurred following Hurricane Gustav in 2008, with ΔDhigh > 0.25 m in 3-13% of profiles across R1-R4.The following year, ΔDhigh > 0.25 m occurred in 22% of profiles in R3, with minimal exceedance of this threshold (<2%) in R1, R2, and R4.After 2010, ΔDhigh > 0.25 m was minimal (<2% for all regions and surveys).Comparing experiment E6 to the identical standalone run (experiment E5), RMSE decreased in all regions with the largest improvement (0.17 m) in R3 (Figure 7, Figure 9).This region also saw a slight improvement in magnitude of model bias (0.05 m) and R2 (0.09), while error statistics in the other regions negligibly improved or degraded slightly (e.g., increase in magnitude of negative bias in R2).

Discussion
EDGR relies on several assumptions that define the scope of model applicability, including that (1) the cross-shore structure and alongshore variability of a barrier island can be characterized by spatial decomposition with a basis set of Gaussian curves, and (2) this spatial decomposition supports a temporal process decomposition, separating dune growth from other coastal changes.In the following sections, we discuss the validity of these assumptions and draw conclusions on the generalizability of the modeling approach for broader scientific, engineering, and coastal management applications.

Dune Characterization
The EDGR model may be applied when a scale separation allows the forms of the barrier island platform and foredune to be isolated from each other.In its present formulation, the model uses the PIGF method and assumes Gaussian forms for the platform and any dune or berm features.Additional testing at other sites is required to test the validity of this approach at other barrier islands and to determine if alternate approaches can be used at sites where features are not well approximated by a Gaussian form.For example, a linear approximation for the residual elevation could be applied in application to mainland coasts, or a half-waveform model could be used to represent the ocean-facing half of the domain [71].The PIGF method can also be modified to use a piecewise function approach to capture different regions of the cross-shore profile, such as using a linear fit to approximate the beach [59].
A benefit of spatial decomposition is that it enables EDGR as a mechanism for characterizing dunes in addition to predicting growth.Applying a dune template and benchmarking against historical data allows evaluation of the similarity of the study site to the template and sensitivity testing to dune drivers and response.For example, the bias in ΔSL,0 became more positive moving from west to east in all experiments at Dauphin Island using FFT-randomization, reflecting formation of the dunes somewhat closer to the post-Katrina shoreline where the island narrows in width.Dhigh did not exhibit the same trend and was overpredicted in R1 and R3 when R1-R4 or R4 was used as the template, and underpredicted in the other two regions.Dhigh bias became more positive in all four regions when known terminal dune locations were prescribed in the model (experiment E5) rather than using template locations (experiment E2) despite using an identical distribution of Dhigh,F.This variation was found to have occurred because the 2005 elevation tended to be higher at the location of actual foredune formation than at the template location, resulting in larger values of Dhigh,I.However, the foredune did not typically form at the overall elevation maximum in each post-storm cross-shore profile (captured in the error in L2005 dune location).Incipient dunes are known to form at the site of flow perturbations, which include local elevation highs.The results at Dauphin Island suggest there may be a site-specific window of potential dune formation constrained by distance to shore.
Decomposition of lidar profiles into a Gaussian basis set for parameterization of mdWdH and bdWdH (Equation 4) allowed rapid and automated analysis of the relationship between dune width (Gwide) and height above the residual island (Ghigh).There was considerable profile-to-profile variation in the relationship of Gwide to Ghigh, but a temporal pattern emerged wherein low dunes in the first 1-2 years after Katrina had a relatively steep relationship between Gwide and Ghigh (Figure 10).As the dunes grew over time, they also became narrower for a given height (decrease in mdWdH).The relationship between Gwide to Ghigh followed this trend over the first 6-8 years following Katrina, after which the relationship stabilized in most of the analysis regions.This result is consistent with prior findings suggesting that the dune ramp-a relatively shallow sloping feature at the base of the dune that is incorporated in the Gaussian form-is a source of sediment for dune crest growth [72].

Dune Evolution
A requirement for EDGR application is the availability of sufficient data for model initialization, which can be a challenge for coastal models in general [73].The skill of the model using R3 and R4 as the template (experiments E3 and E4) suggests that the minimum longshore spatial scale of data needed to parameterize the model is on the order of 2-4 km.Data on terminal dune state and the time scale of post-storm dune recovery are required, while testing of the influence of the growth rate parameter suggests that the model is relatively less sensitive to the availability of interim data points during post-storm recovery.In cases where local dune recovery data are unavailable, the limited impact on prediction error for a given region of the study site when using a different region as the template indicates that dune data may be transferred effectively between similar sites using the FFT method.In the application at Dauphin Island, EDGR outperformed the no-change model (L2005) within each region in terms of bias and RMSE even when R3, the region most dissimilar to the rest of the study site, was used as the template (experiment E3).Data from a prior time period are also likely to be an appropriate template for future growth at any given site under the assumption that recovery from different storms will be broadly similar.It may be possible to generate template values from a more distant location when dune recovery data at a target site are limited or unavailable, thus increasing the general applicability of the model to include sites with insufficient historical data for model parameterization.Gutierrez et al. [74] found that Dhigh is correlated to factors such as the shoreline change rate and beach width, suggesting that such information can be used to constrain the choice of template at sites where long-term dune evolution data are unavailable.Factors known to control dune building rate and form (vegetation type, dominant wind direction and speed, sediment grain size, etc.) should also be assessed and compared between the study area and potential EDGR template sites.
Where sufficient data are available, EDGR can be used to evaluate the relative contribution of processes of erosion and accretion to dune evolution.Using R3 as the template (experiment E3) increased the negative bias in Dhigh for R2 and R4 and introduced a negative bias in R1, while excluding R3 from the template (experiment E4) improved model performance in all the other regions.These results suggest that either R3 has the lowest terminal dune heights of the study site (i.e., Dhigh,F is constrained by dune-building processes) or frequent overwash inhibits dunes from approaching their terminal state (i.e., Dhigh is constrained by erosion).The relative influence of these processes can be considered by comparing otherwise identical standalone (experiment E5) to lidarcoupled (experiment E6) model predictions.From 2005 to 2010, R3 had the highest percentage of profiles with dune elevation loss between lidar surveys exceeding 0.25 m (Figure 4).In particular, R3 exhibited widespread erosion in the 2008-2009 time period that was not observed in the other regions.Coupling to lidar reduced the RMSE and bias, and increased R2 for R3 over the period of 2005-2012 (Figure 9).From 2010 to 2015, few profiles in any of the regions exhibited loss of Dhigh greater than the 0.25 m, indicating limited overwash during this time.However, positive bias in R3 predictions began to grow in the coupled model after 2012 and exceeded the bias in the standalone run using R3 as the template (experiment E3).This result suggests that reduced capacity for dune building is also contributing to the relatively low dune heights found in R3.Saltwater inundation impacts vegetation type, resulting in a secondary effect of storms on dune evolution [75].Therefore, one possible explanation may be that the higher frequency of overwash during the 2005-2010 period altered the local vegetation.Alternately, sediment supply and/or other factors controlling Aeolian flux may differ in this region.The coupled model is also informative as to what types of events are driving significant dune erosion at Dauphin Island.Two periods of widespread elevation loss (in 2008 and 2009) corresponded to passage of tropical storms within 200 km of Dauphin Island, while elevation loss in 2007 suggests that extratropical storms and/or far-field tropical storms have also resulted in overwash over portions of the island.
The application of EDGR to investigate dune growth also highlighted potential model modifications.In the current version, incremental changes to the dune profile from an external erosion model (those which do not reset Dy or Dhigh,0) are discarded in a coupled framework.Dune elevation change in EDGR could alternately be formulated based on differential dune growth [11] as Positive elevation change (from EDGR) and negative elevation change (from an erosion model) could then be explicitly resolved in a coupled framework.Doing so would require iterative calibration of the growth parameter to describe uninhibited dune elevation gain, which cannot be derived from observational field data that includes times of overwash and inundation.Dune width change could similarly be reformulated and calibrated within EDGR as dDwide/dt, which may resolve some of the challenges found in identifying a robust relationship between Ghigh and Gwide.Reformulation in this manner would further segregate the impacts of Aeolian and hydrodynamic sediment transport on dune evolution, thus facilitating explicit examination of these processes and the factors that control them in addition to potentially improving model performance.
The EDGR model could also be modified to replace the Houser dune growth formulation with a formulation that included time-varying environmental conditions that influence dune growth [14,[21][22][23][24][25][26][27][28].For example, growth rate (r) could be formulated as a function of predicted Aeolian sediment flux, which could either be derived from an existing process-based model (see, e.g., in [37,38]) or parameterized as a function of wind speed, fetch length, and sediment bed characteristics (grain size, moisture content).Vegetation growth and sediment trapping [14,37,47] could similarly be explicitly modeled or characteristics of dune form (width, height) parameterized as a function of dominant vegetation type [53].Modification of the model formulation to include explicit linkage to environmental characteristics or processes correlated to dune growth would replace the need for historical dune growth data with a requirement for, e.g., site-specific wind and vegetation data to predict dune growth parameters.

Broader Model Application
When initial and terminal dune characteristics are known and prescribed (experiment E1), EDGR predicts the interim evolution of dune height and the cross-shore profile of the dune.In this case, the model has the highest accuracy in predicting spatially and temporally variable dune elevations.Even when the precise spatial distribution of terminal dune characteristics is unknown, EDGR captures spatially averaged characteristics in dune parameters and had lower region-wide bias than the no-change alternative model (Figure 7).Capturing island-scale dune statistics (experiments E2-E6) is of similar value to predicting the precise longshore distribution (experiment E1) in some cases.For example, dune habitat may be evaluated based on regional or island-wide acreage [65].Similarly, estimates of the probability of island-scale overwash and inundation [2] may be improved by updating estimates of dune maximum elevation.In cases where the longshore distribution of dune state is required [8,76], EDGR's computational efficiency enables extensive sensitivity testing and ensemble prediction to constrain possible future scenarios.For reference, the 1918 cross-shore profiles initialized in under 30 s and the model completed simulations for 2005 to 2015 in less than 5 min while running on a Windows desktop computer with 2.4 GHz Intel Xeon CPU.This efficiency enables a wide range of dune outcomes and associated uncertainty to be generated and evaluated.Last, EDGR may be coupled with other models that evolve the storm impacts and other changes to the coast not associated with the foredune as part of a model framework for predicting decadal-scale coastal evolution [56][57][58].Therefore, the EDGR method is readily applicable to coastal vulnerability and restoration studies and associated management applications.

Conclusions
In the current study, an empirical dune growth model (EDGR) was developed and used to predict the recovery of the primary foredune of a barrier island following erosion during storms.Cross-shore profiles are first approximated as the sum of a set of Gaussian curves representing dune and berm features atop an underlying island form using a novel parameterized island Gaussian fit (PIGF) method.The foredune Gaussian is evolved with maximum height following the Houser et al. [11] growth model and width calculated from a new linear formulation relating dune width to height.Terminal dune location and height can be prescribed if known, such as when EDGR is being used to fill temporal data gaps between post-storm and recovered lidar surveys.Alternatively, EDGR can create a longshore distribution of terminal dune parameters via fast Fourier transform randomization of available template data from a region that is assumed to be similar to the study site.Bulk characteristics (mean and standard deviation) and length-scales of spatial variability from the template are preserved, while EDGR's high computational efficiency allows multiple potential realizations to be rapidly generated and tested to constrain the range of possible outcomes.EDGR can operate as a standalone model or be coupled with an external model that predicts dune erosion and/or evolution of other characteristics of the barrier island, such as the shoreline.
EDGR was used to predict foredune evolution of the western portion of Dauphin Island over the 10-year period following widespread dune destruction by Hurricane Katrina (2005) and compared to 12 topographic lidar surveys.Model predictions were benchmarked against 2005 lidar, representing elevation values that would be used if updated lidar were unavailable.EDGR improved predictions of maximum dune elevation over the no-change model, with highest performance when known terminal dune characteristics were prescribed.This application also illustrated the potential of the model to be used in investigating dune growth processes.Use of the Gaussian basis set to extract dune width and height for model parameterization revealed that the relationship between dune height and width evolved as the dunes grew; over time, dunes became relatively narrower for a given height.In addition, comparison between a standalone model configuration and one wherein EDGR was coupled to lidar surveys as a proxy for a dune erosion model suggests that dune heights in a lower elevation portion of the study site were constrained by erosion in the first few years following Katrina.Relatively limited elevation gain compared to other portions of the study site in more recent years, when dune crest erosion was not identified in the lidar, suggests that the area has lower terminal heights, and is thus now growth-limited.This application illustrates that EDGR is a Table A1.Root mean square error (RMSE), bias, and coefficient of determination (R 2 ) of dune position referenced as distance to the 2005, post-storm shoreline (ΔSL,0) for six configurations of EDGR and an alternate, no-change model taken as the 2005 lidar (Table 3).

Figure 1 .
Figure 1.(A) Location of Dauphin Island within (B) the northern Gulf of Mexico.

Figure 2 .
Figure 2. (A) Diagram illustrating the Gaussian decomposition method and empirical dune growth model (EDGR) formulations for growth in (B) dune height and (C) width.(D) The evolving foredune is incorporated into the rest of the cross-shore profile as it evolves.Model formulations may be found in Equations 1-5 with descriptions of model parameters provided inTable 1. Figure from the work in [56].

Figure 3 .
Figure 3. Elevation of Dauphin Island from light detection and ranging (lidar) surveys taken in (A) 2005 and (B) 2015; also shown are analysis regions used in the current study (dashed black lines),

Figure 4 .
Figure 4. (A) Lidar surveys (black lines) taken over the study time period of August 2005 to January 2015, along with tropical storms or hurricanes that passed within 200 km of Dauphin Island (black dots).Survey data are not available for R1 and R2 for the March 2006 survey.(B) Percentage of crossshore profiles for each analysis region for which the loss in elevation of dune height (ΔDhigh) was great than 0.25 m compared to the previous survey.Taken from the work in [56].

Figure 5 .
Figure 5. (A) EDGR 2015 model prediction for experiment E1 (Table3) and (B) difference between the prediction and 2015 lidar.From within each region R1-R4 (demarcated by dashed black lines), a representative cross-shore transect has been selected (shown as colored lines in panels (A,B) and denoted as T1-T4).The initial (gray line) and 2015 lidar elevation profiles (black line) are shown for each of these transects (panels (C-F)), along with the EDGR prediction for 2015 (colored dashed line, where the color for each transect is the same as in panels (A,B).

Figure 6 .
Figure 6.(A,D) Root mean square error (RMSE), (B,E) bias, and (C,F) coefficient of determination (R2) for experiment E1 (Table3).Dhigh (left column) was taken as the maximum dune elevation for each cross-shore profile.The same statistics are shown for the 2005 lidar.Statistics for Zdune(x,y,t) (right column) are calculated at all cross-and longshore locations where EDGR updated the elevation for a given lidar survey.

Figure 7 .
Figure 7. (A,D) Root mean square error (RMSE), (B,E) bias, and (C,F) coefficient of determination (R2) for model predictions using six experiment configurations of EDGR (E1-E6; Table3) of Dhigh (A-C) and ΔSL,0 (D-F).Statistics are calculated over all lidar surveys from 2006 to 2015.For EDGR configurations using template dunes, the value shown is the mean over the 100 realizations, with the vertical bar indicating plus and minus one standard deviation.Error statistics for ΔSL,0 are not shown in cases (gray) where dune location was prescribed using historical data.

Figure 9 .
Figure 9. (A-E) Root mean square error (RMSE), (F-J) bias, and (K-O) coefficient of determination (R2) of Dhigh for EDGR applied to predict dune growth from 2005 to 2015.Each row represents a different configuration of EDGR (Table3), with model predictions benchmarked against a no-change model represented by 2005 lidar (L2005).

Figure 10 .
Figure 10.Relationship between the height (Ghigh) and width (Gwide) of the foredune identified via sumof-Gaussians fitting of observed cross-shore profiles (Figure 2).Shown are (A) September 2006 (~1 year post-Katrina) and (B) January 2015 (~10 years post-Katrina).Boundaries of the analysis regions (R1-R4) are found in Figure 3.The relationship between Ghigh and Gwide was found to vary in time (linear model slope, mdWdH, and intercept, bdWdH, shown in panels (C,D).This temporal variability was beyond the scope of this initial study and a single relationship (black line in panels (A,B)) was used in all EDGR applications.

Table 1 .
Figure from the work in [56].

Table 1 .
Description of EDGR model parameters.Model formulations may be found in Equations 1-5 with a visual representation of model application given in Figure2.
mdWdH Slope of the linear model used to predict Gwide(x,t) from Ghigh(x,t) bdWdH Y-intercept of the linear model used to predict Gwide(x,t) from Ghigh(x,t)

Table 2 .
EDGR model parameters calculated for each analysis region (Figure3) based on lidar surveys in 2005 and 2015.Model parameters include initial elevation, Dhigh,0; terminal elevation, Dhigh,F; growth rate, r; distance to the 2005 shoreline, ΔSL,0; and mdWdH and bdWdH, the slope and y-intercept of the linear model used to estimate foredune Gaussian width from height.

Table 3 .
Model configurations used to predict dune growth at Dauphin Island from 2005 to 2015.