Exploring Marine and Aeolian Controls on Coastal Foredune Growth Using a Coupled Numerical Model

Abstract

Coastal landscape change represents aggregated sediment transport gradients from spatially and temporally variable marine and aeolian forces. Numerous tools exist that independently simulate subaqueous and subaerial coastal profile change in response to these physical forces on a range of time scales. In this capacity, coastal foredunes have been treated primarily as wind-driven features. However, there are several marine controls on coastal foredune growth, such as sediment supply and moisture effects on aeolian processes. To improve understanding of interactions across the land-sea interface, here the development of the new Windsurf-coupled numerical modeling framework is presented. Windsurf couples standalone subaqueous and subaerial coastal change models to simulate the co-evolution of the coastal zone in response to both marine and aeolian processes. Windsurf is applied to a progradational, dissipative coastal system in Washington, USA, demonstrating the ability of the model framework to simulate sediment exchanges between the nearshore, beach, and dune for a one-year period. Windsurf simulations generally reproduce observed cycles of seasonal beach progradation and retreat, as well as dune growth, with reasonable skill. Exploratory model simulations are used to further explore the implications of environmental forcing variability on annual-scale coastal profile evolution. The findings of this work support the hypothesis that there are both direct and indirect oceanographic and meteorological controls on coastal foredune progradation, with this new modeling tool providing a new means of exploring complex morphodynamic feedback mechanisms.

1. Introduction

Coastal landscape evolution reflects the aggregation of the combined effects of oceanographic, meteorological, geological, ecological, and anthropogenic influences [1]. Because of the broad range of physical processes driving temporal variability of morphology in the coastal zone e.g., [2,3,4,5], the coastal profile has often been compartmentalized into discrete morphologic units (e.g., nearshore, beach, dune) [6,7] based on dominant transport processes (e.g., waves, currents, winds). Although dunes backing sandy beaches have largely been treated as wind-controlled features, marine-driven processes play an active role in the accretional and erosional development of these landforms e.g., [8,9]. In part reflecting these complex interactions, quantitative predictions of coastal dune evolution at scales beyond individual storm events are lacking e.g., [10,11] despite the important suite of ecosystem services that coastal dunes provide (e.g., coastal protection, habitat, recreation; [12,13,14,15]).

Much of the current understanding of coastal dune evolution remains largely conceptual e.g., [2,8,11,16,17,18]. Attempts that quantify dune growth and recovery from storms come primarily from scaling up small-scale wind-driven sediment dynamics e.g., [19,20] and explorations of ecomorphodynamics e.g., [21,22,23]. Although aeolian transport is the primary driver of coastal foredune growth, in many real-world systems wind climatology is not well correlated with dune growth rates [24,25,26,27], making empirical prediction of dune evolution difficult.

The marine environment poses numerous constraints on instantaneous aeolian transport rates in sandy coastal systems, contributing to the nonlinear relationships between dune growth and climatic variables. As examples, swash-induced moisture [28], salt cementation [29], and tidally modulated groundwater tables [30] each have been found to reduce or prevent local aeolian transport. These supply limiters contribute to the widely observed fetch effect [31,32,33] whereby there is a minimum length-scale required for aeolian sand transport on beaches to reach saturation. Obliquely oriented wind results in larger effective fetch distances between the water line and the dune toe [31,34]. Oblique winds therefore will transport more sediment over the beach as compared to cross-shore oriented winds of equal magnitude in fetch limited circumstances [35]. Conversely, the cross-shore component of aeolian transport is highest under shore-normal winds and reduces with increasing obliquity, leading to a potential decrease in sediment delivery to the dune under oblique wind conditions due to the cosine effect [31]. In part due to a lack of fetch limitations [36], larger beach widths associated with low-gradient, dissipative beaches promote higher rates of aeolian transport, leading to larger dune growth relative to narrower, steeper beaches [2,8,22]. The relationship between grain size and beach slope also generally encourages higher transport rates on low sloping beaches [37,38].

Since fine grained sediment is preferentially transported by wind over coarse-grained sediment, sediment sorting and armoring may subsequently occur on mixed-grain beaches e.g., [33,39,40]—leading to complex spatial patterns in instantaneous aeolian transport. In part due to armoring effects on the upper beach, it has been found that the intertidal zone is a primary sediment source area for backshore aeolian transport [41,42,43,44]. The onshore propagation of sandbars into the intertidal zone has been recognized as a particularly important source of sediment for dune growth e.g., [26,45]. In fact, it has been hypothesized that dune growth is limited by the temporal co-occurrence of intertidal sandbar welding events with the capacity to mobilize sediment by wind [26]. However, field data from a dissipative beach have demonstrated that dune growth can occur regardless of in-phase synchronization of maximum beach sediment supply and maximum dune growth [46].

Energetic wave conditions and high water levels in the collision regime of the Sallenger [9] Storm Impact Scaling Model have often been related to the erosion of coastal foredunes e.g., [47,48]. The exposure of coastal dunes to erosion is related to a variety of morphologic and environmental controls e.g., [49,50], although low sloping, wide beaches are generally less vulnerable to dune erosion than steep, narrow beaches. It has also been recently demonstrated, using both field data [46] and numerical models [51], that infragravity swash processes can contribute directly to dune growth on dissipative beaches. These results give further evidence of the direct and indirect roles of the marine system in governing both accretive and erosive dune processes.

Despite a qualitative understanding of the dominant marine and aeolian processes controlling coastal foredune evolution, quantitative models able to reliably simulate dune growth in real-world coastal systems, including the recovery of coastal dunes following storm-induced erosion, are lacking [10,11]. Recognizing that major gaps in knowledge of sediment transport and morphological change across the land-sea interface still exist, here we introduce the coupled coastal profile modeling framework Windsurf. Windsurf couples three state-of-the-art open-source models that independently account for subaqueous (XBeach; [52]) and subaerial (Aeolis; [33]; Coastal Dune Model (CDM); [22]) processes to explore the evolution of the coastal profile in response to both marine and aeolian forcing (Figure 1). The numerical models within Windsurf are capable of simulating complex coastal behaviors such as beach growth [53], swash dynamics [54], wave-driven dune growth [51], storm-induced dune erosion [55], overwashing [56], aeolian sediment sorting [33], and ecological controls on dune growth [22]. In this paper, we aim to test this new Windsurf model coupling framework. Windsurf is used to hindcast prograding and dissipative conditions for a one-year period at a data rich field site in the U.S. Pacific Northwest (PNW). This work serves as a proof of concept of Windsurf, with the model outputs additionally used to gain new insights into the relative controls of marine and aeolian processes on coastal foredune growth.

The structure of this paper is as follows. Details of the new Windsurf modeling framework and its sub-models are described in Section 2. Morphologic and environmental field datasets used to validate Windsurf are presented in Section 3. Windsurf simulations completed for the PNW field site are described in Section 4. Discussion of the results and conclusions are given in Section 5 and Section 6, respectively.

2. Windsurf Model Framework

The following section describes the details of the Windsurf model coupler, followed by details of the process capabilities of the individual numerical models within Windsurf.

2.1. Model Coupler

Windsurf couples three separate open-source numerical models [henceforth referred to as model cores] which simulate subtidal morphodynamics related to waves and currents (XBeach), subaerial morphodynamics related to wind shear and vegetation (CDM), and multi-fraction aeolian sediment transport that includes the effects of supply limiters (Aeolis) (Figure 1 and Figure 2). Together, these model cores simulate the evolution of the full coastal profile in response to both marine and aeolian forcings.

Windsurf serves as the back-end coupler to these three model cores, with functionality to generate input files, execute model simulations, exchange information between model cores, and save model output. The Windsurf framework runs the individual model cores in series, with morphological, environmental, and ecological information exchanged between the models at user-defined time steps (Figure 2). By using a loosely coupled approach in which the model cores are separated from the coupler, the individual models can continue to be developed independently and new model components can easily be incorporated into Windsurf. New capabilities can be added if the additional model core can run on or be interpolated to a central model grid. Currently, Windsurf model grids are limited to 1D, cross-shore only applications; however, extension to 2D is planned for future work.

The work presented in this manuscript uses a version of the Windsurf framework coded in MATLAB (v2017a or newer) which runs the process-based model cores in an offline mode. That is, the model cores are initialized and fully executed before moving on to the next model core. As schematized in Figure 2, Windsurf first runs XBeach to simulate marine-related processes and resulting morphology change. XBeach runs are followed by CDM which simulates spatio-temporal vegetation dynamics and the cross-shore wind field, including the effects of dune grass vegetation on bed shear stress reduction. Finally, Aeolis simulates aeolian sediment transport and aeolian-driven morphology change. A link to the Windsurf model coupling code can be found in the Supplementary Materials.

2.2. XBeach

XBeach is an open-source process-based numerical model which simulates nearshore wave transformation, subaqueous sediment transport, and morphology change [52]. XBeach has been used in numerous studies investigating wave transformation [57,58], swash processes [46,54,59], and morphology change [55,56,60,61,62]. Originally developed for assessing impacts from discrete storm events, the model has more recently been applied to simulate marine-driven morphologic changes on time scales of weeks to years [53,63,64]. An overview of relevant processes and model parameters for XBeach are described below, while a full description of the model and its formulations can be found in Roelvink et al. [52] and Roelvink et al. [65].

When using XBeach in surfbeat mode, the model does not resolve individual short waves but does resolve the wave group envelope using the time dependent, short-wave action balance equation [52,66]. The short-wave action balance equation is also coupled with the nonlinear shallow water equations to account for infragravity waves and mean currents. XBeach, which is a depth-averaged model, solves these spatially varying hydrodynamic processes in either 1DV (cross-shore profile model) or 2DH (area model).

Sediment transport is solved by an advection-diffusion equation [67]. In surfbeat mode, as XBeach is presently applied in the Windsurf framework, wave groups are resolved but individual incident (short) waves are not modeled. As the wave shape is not directly simulated, the effect of wave nonlinearity on sediment transport is parameterized via the short-wave-related velocity asymmetry (facAs) and skewness (facSk) coefficients. Many studies provide guidance on the effects of facAs and facSk on simulating erosional and accretional processes with XBeach e.g., [68,69,70]. Mean currents and long waves also contribute to sediment transport. Spatial gradients in the short-wave and current related sediment transport rates simulated by XBeach are used to calculate the resulting bed level changes using a form of the Exner equation. McCall et al. [56] found that morphology change associated with overwash and other high flow currents are often overpredicted. To address this issue, the dilatancy option accounts for the influence of pore water on bed stability, effectively increasing the local critical Shields parameter and thereby reducing transport rates during high flow conditions following the methods of van Rhee [71]. A morphologic acceleration factor (morfac) can also be implemented in XBeach according to the methods of Reniers [72]. Morfac allows morphology change to be scaled, enabling shorter duration numerical simulations without significant loss in process resolution e.g., [68,73].

2.3. Coastal Dune Model

The CDM can simulate the evolution of vegetated coastal foredunes by solving a suite of differential equations describing ecomorphodynamic processes [22]. The model originated as a 1D dry-environment saltation model [74] to simulate the formation and migration of desert dunes. Later work extended the model to a 2D area model and added the role of vegetation in stabilizing the dune form [75]. The model has been used to explore the formation mechanisms and time scales of linear [76], parabolic [75], and barchan [77] dunes. The latest iterations of CDM have been applied to simulate controls on coastal foredune ridge development [78], post-storm dune recovery on barrier islands [79], and controls on foredune hummockiness [80]. A key driver of these different morphological forms is resolving the feedback of the topographic field on the local wind field. The computationally efficient wind model of Weng et al. [81] that is implemented in CDM builds upon Hunt et al. [82], an analytical solution for turbulent shear flow over a low sloping hump, to solve for the spatially variable bed shear perturbation field ($\delta \tau$) on the beach and dune.

In addition to resolving wind and aeolian transport dynamics, CDM tracks the spatial and temporal growth of dune vegetation by assuming a linear growth rate of the vegetation cover fraction (${\rho }_{veg}$; [22]). The effective bed shear stress (${\tau }_s$) in the presence of vegetation, which is calculated by CDM, becomes:

$${\tau }_s=\frac{{\tau }_o+{\tau }_o\delta \tau }{1+\frac{m\beta }{\sigma }{\rho }_{veg}}$$

(1)

where ${\tau }_o$ is the bed shear stress in the absence of variable topography or vegetation, m is the vegetation friction coefficient, $\sigma$ is the ratio of the plant basal to frontal area, and $\beta$ is the ratio of the plant drag coefficient to bare sand [83].

CDM has the capability to simulate aeolian sediment transport using the derived bed shear stress field. The model assumes a single (spatially constant) grain size and has shown success in simulating coastal dune dynamics in a variety of morphodynamic settings e.g., [22,78,80]. When forced with cross-shore directed winds, CDM model results have suggested that LVeg, defined as the distance between the shoreline and the perennial vegetation line, can exert a primary control on the maximum height which dunes can grow [22]. Since the shoreline position is typically non-stationary in real-world settings, in Windsurf LVeg is modified to allow for vegetation growth that follows a fixed vertical contour position rather than a fixed beach width. That is, as beaches and dunes erode or prograde, the seaward-most location of vegetation in the model is assumed to follow the seaward-most location of a fixed vertical contour rather than remain at a fixed horizontal position. Furthermore, to allow for obliquely oriented winds in Windsurf, the cross-shore oriented CDM grid is rotated at each coupling time step based on the wind direction to properly resolve the $\delta \tau$ field over the foredune. Within Windsurf, CDM is used to compute the spatially varying topographically and ecologically influenced bed shear stress field. CDM is not used directly for aeolian transport calculations within Windsurf in part because only a single grain size can be implemented in the model. Instead, using the CDM calculated spatially variable bed shear stress field, multi-fraction aeolian sediment transport is simulated using Aeolis within the Windsurf framework. For additional information on CDM and its formulations see Durán and Moore [22], Moore et al. [78], Durán and Moore [79] and Goldstein et al. [80].

2.4. Aeolis

Aeolis is an aeolian sediment transport model specifically designed to simulate supply limiting processes [33,84]. As wind processes entrain and transport sediment via saltation, the model accounts for the temporal and spatial evolution of the grain size distribution at the bed surface, and with depth, on either a 1D (cross-shore) or 2D domain. Within each model grid cell, local wind properties (bed shear stress and wind velocity) are used to calculate sediment transport rates for each grain size bin ($d_n$). The saturated sediment mass-flux transport rate ($q_{sat}$) is calculated for each $d_n$ based on the formulation of Bagnold [38]:

$$q_{sat}=C_b\frac{{\rho }_a}{g}\sqrt{\frac{d_n}{D_n}}{\left(u_{*}-u_{*th}\right)}^3$$

(2)

where ${\rho }_a$ is the air density, $D_n$ is a reference grain size of 0.25 mm, $u_{*}$ is the shear velocity, $u_{*th}$ is the threshold shear velocity below which transport does not occur, and $C_b$ is an empirical coefficient. The shear velocity threshold for transport is calculated via Bagnold [38] and modified incorporating the effect of moisture following Belly [85] ($f_m$):

$$u_{*th}=f_{m}A\sqrt{\frac{{\rho }_s-{\rho }_a}{{\rho }_a}g{d}_n}$$

(3)

where ${\rho }_s$ is the sediment density and A is an empirical coefficient. For the present application moisture-related effects ($f_m$) are only used in the swash zone to limit aeolian transport in the presence of high water levels; precipitation and groundwater–related moisture are not yet implemented in Windsurf but will be accounted for in future versions. A weighted saturated sediment concentration is calculated using the local bed grain size distribution. Sediment transport gradients, resulting from local sediment characteristics, supply limiting effects (e.g., moisture, salt content), and/or the spatially variable wind field, produce local bed elevation changes.

The Aeolis model has shown skill at simulating complex sediment sorting and beach armoring processes [33]. Because of the added capability of simulating multi-fraction aeolian sediment transport and the inclusion of other supply limiting effects, Aeolis is used as the subaerial transport model within Windsurf. Windsurf provides Aeolis with spatially varying wind and bed shear stress fields from CDM and the total water level from XBeach at each time step. Aeolis is run using the full bed shear stress derived from CDM to calculate aeolian sediment concentrations which are advected using only the cross-shore component of wind field for the 1D Windsurf simulations described in this manuscript. For additional information on the details of Aeolis, see Hoonhout and de Vries [33].

3. Field Observations at a Prograding Coastal System

To test the ability of the new modeling framework to simulate realistic behavior, Windsurf is applied to a prograding coastal field site in the PNW. The following sections briefly describe the field site and the relevant datasets for model implementation.

3.1. Field Setting

The Long Beach Peninsula, WA (LBP; Figure 3) is a modally dissipative, mesotidal (2–3 m tidal range) coastal system in the PNW characterized by low-gradient (backshore slopes of ∼0.025), sandy beaches, and linear foredune ridges densely colonized with Ammophila breviligulata (American beach grass; [21,86]; Figure 3b,c). In the Oysterville section of LBP, the beach is rapidly prograding (>4 m/year) in response to gradients in longshore sediment transport and cross-shore sediment feeding from the shoreface [7]. The dune complex is also prograding, with approximately 10–15 m${}^3$/m/year of dune growth [46] which contributed to the generation of a new foredune (the seaward-most dune) since the late 1990s [78]. The dune crest elevation at Oysterville is about ∼9 m NAVD88, corresponding to about 7 m above mean high water (MHW, 2.1 m NAVD88).

Long-term dune growth at this site is driven primarily by wind-driven processes [46]. The outer Washington coast is characterized by a seasonal wind climate, with the largest average monthly wind speeds in January and December. The highest wind speed events typically co-occur with large wave events. While high total water levels (TWLs), driven in part by large wave heights, regularly impact the dune toe at Oysterville, WA, previous observational [46] and modeling [51] studies have indicated that these TWLs may also directly contribute to growth of the lower portion of the dune at Oysterville.

High TWLs at Oysterville are primarily a result of large wave runup associated with an energetic wave climate in the PNW. The average deep-water winter significant wave height (H${}_s$) offshore of the PNW is about 2.3 m, with numerous storms each year exceeding 8 m [86]. These high energy conditions seasonally erode the foreshore, leading to oscillations of the MHW shoreline of over 30 m despite a net annual progradation of the system [86,87]. Beach recovery occurs predominantly in the low energy summer period (H${}_s$∼1.5 m on average) in response to the onshore propagation and welding of intertidal sandbars [46,88].

The surf zone at the study site is typically characterized by numerous (2–4) subtidal sandbars which vary significantly on interannual time scales [89] as captured by a coastal monitoring program that has been ongoing since the late 1990s [86]. Recent work has focused on characterizing sub-annual coastal evolution at Oysterville [46,88]. The existence of morphology data at a wide range of time scales at a site with large morphology change signals makes Oysterville an ideal location to test the new Windsurf modeling framework.

3.2. Morphology Data

3.2.1. Annual Scale

Topographic and bathymetric measurements have been made at Oysterville using real time kinematic GPS surveying techniques [86]. A single cross-shore transect that extends from 9 m water depth to landward of the foredune crest is used here. The measured cross-shore transect data are used for input to the numerical model and for comparison to model output. From these data, volumetric changes were calculated between morphological units, where the nearshore is defined here as the region from −9 to 1 m NAVD88 [all vertical references henceforth are relative to NAVD88]. The upper limit of the nearshore zone corresponds approximately to local mean sea level at the field site. The beach is defined here to be between the 1 m and 4 m contours and the dune region includes all areas above the 4 m contour, where 4 m is the approximate dune toe position. Measurements between 4 August 2016 and 9 August 2017 show nearshore ($\Delta V_{nearshore}$), beach ($\Delta V_{beach}$), and dune ($\Delta V_{dune}$) volume changes of −166 m${}^3$/m (∼±300 m${}^3$/m), −5.4 m${}^3$/m (∼±13 m${}^3$/m), and 16.7 m${}^3$/m (∼±5 m${}^3$/m), respectively (Figure 4 and Figure 5). Please note that the listed volumetric change uncertainties are calculated by combing the random vertical uncertainties associated with topographic (∼8 cm; e.g., [86]) and bathymetric (∼13 cm; e.g., [90]) field measurements in quadrature.

The data show that on annual scale there is variability in the subtidal sandbar configuration (Figure 4). The net sediment loss from the nearshore that was not gained by the beach reflects either longshore gradients in transport or large errors associated with volumetric estimation from bathymetric data. Between the two profile measurements, there was also a net 10 m retreat of the MHW shoreline ($\Delta X_{MHW}$), with this apparent shoreline change occurring in part because of configurations of intertidal sandbar troughs at the time of the surveys (Figure 4). This indicates that despite measured dune growth at the site, there was effectively no net beach change over this particular year (within the error of the instrumentation).

3.2.2. Sub-Annual Scale

An additional 9 topographic surveys were completed between the August 2016 and August 2017 survey dates which are also used to inform the timing and magnitude of beach and dune morphology changes at sub-annual scale. Details of the morphology change on these time scales, using data from the same field site, are described in Cohn et al. [46]. Because MHW shoreline change is often used as a metric for coastal change e.g., [91], temporal changes in the MHW shoreline are also extracted from the data (Figure 5a). These data show that between the August 2016 and August 2017 profile there is seasonality in the beach behavior. The beach was widest in summer (July–September) and fall (October–December) and narrowest in winter (January–March) and spring (April–May), with about 40 m variability in $\Delta X_{MHW}$ over the full year. Volume changes in the beach and dune are also calculated at this monthly time scale (Figure 5b,c). $\Delta V_{beach}$ has the same seasonal behavior as $\Delta X_{MHW}$. Conversely, the dune shows gradual volumetric growth throughout the year. The dune volume increases or remains constant between each survey (within the measurement uncertainty), with the lowest rates of dune growth occurring in summer.

3.3. Environmental Data

3.3.1. Water Levels

Offshore waves and still water levels (SWL) are the primary drivers of marine-driven morphology change and therefore serve as important model boundary conditions. A Nortek Acoustic Wave and Current Profiler (AWAC) was deployed on a bottom mounted mooring at −9 m providing wave and SWL measurements for a 42-day period during summer 2016 as part of the Sandbar Aeolian Dune Exchange Experiment (SEDEX${}^2$) [88]. To temporally extend the record of oceanographic forcings at the study site beyond the SEDEX${}^2$ period, additional nearby measurements from the U.S. National Oceanic and Atmospheric Administration (NOAA) Tides and Currents database, the U.S. National Data Buoy Center (NDBC), and the Coastal Data Information Program (CDIP) are used.

The AWAC-derived SWLs show similar tidal amplitudes (Figure 6a) as compared to measurements from the NOAA Toke Point, WA tide gauge located 20 km to the northeast of Oysterville in Willapa Bay (Figure 3a). Slight (≤1 h) tidal phase offsets between the sites contribute to a root mean square error (RMSE) of 0.29 m between the AWAC and Toke Point SWLs for the SEDEX${}^2$ period. For modeling purposes, hourly tidal measurements from Toke Point, WA are used directly.

3.3.2. Waves

Offshore wave information, including significant wave height (H${}_s$), peak wave period (T${}_p$), and wave direction (D${}_{wave}$), were acquired from CDIP buoy 036, located 35 km to the northwest of the study site in 40 m water depth (Figure 3a). The SWAN model [92] was used to transform waves from the CDIP buoy to the AWAC location. A 2D SWAN grid (dx = 100 m, dy = 100 m) that encompasses the northern half of LBP (Figure 3a) was set up using a regional, high resolution NOAA [93] bathymetric digital elevation model to define the offshore water depths. Standard PNW model configurations e.g., [94] are used for this application. Local wind wave generation was neglected.

SWAN was used to transform waves to the AWAC location for each hour between August 2016 and August 2017 using the CDIP 036 measured H${}_s$, T${}_p$, D${}_{wave}$ and the tide gauge measured SWL (Figure 6b–d,h–j). SWAN performs well at modeling wave transformation to Oysterville, with minimal bias ($\Delta \mu$) and RMSEs of 0.12 m, 1.7 s, and 8.0${}^{\circ }$ for H${}_s$, T${}_p$, and D${}_{wave}$, respectively, relative to the AWAC data (

Table 1. Root mean squared error (RMSE) and bias ($\Delta \mu$) of transformed oceanographic and meteorological measurements, used to drive Windsurf boundary conditions, compared to local measurements. The local marine measurements used for comparison are 42 days and are 291 days for the wind measurements.

Environmental Parameter	Local Measurement Transformed Dataset	Locally	RMSE	$\mathit{\Delta }\mathit{\mu }$
SWL	AWAC	Toke Point	0.29 m	0.04 m
H${}_{\mathit{s}}$		CDIP 46211	0.12 m	<0.01 m
T${}_{\mathit{p}}$			1.7 s	0.03 s
D${}_{\mathit{w}\mathit{a}\mathit{v}\mathit{e}}$			8 $$	0.3 $$
u	Dyacon	Toke Point	1.13 m/s	−0.37 m/s
D${}_{\mathit{w}\mathit{i}\mathit{n}\mathit{d}}$			25.5 ${}^{\circ }$	0.2 ${}^{\circ }$

; Figure 6b–d).

3.3.3. Wind

For the period between 3 August 2016 and 21 May 2017 a Dyacon MS-140 weather station was deployed on the foredune crest at the Oysterville site and provided bulk wind speed (u) and direction (D${}_{wind}$) measurements. To extend the meteorological measurements beyond this local measurement period, available wind data from the NOAA Toke Point wind gauge were used. First, both datasets were transformed to their 10 m equivalent wind speeds. Then a linear transformation function was developed to relate u and D${}_{wind}$ from Toke Point, WA to Oysterville, WA using the overlapping time series. The transformation generally agrees well with local wind speed data (Figure 6e), with a RMSE for the locally predicted wind speed of 1.13 m/s with a small $\Delta \mu$ (−0.37 m/s). However, there are periods when wind direction is poorly resolved, likely because of the local topography within Willapa Bay. Therefore, while 77% of the locally transformed D${}_{wind}$ are within ± 5${}^{\circ }$ of the measurements, there is a fairly high overall RMSE value of 25.5${}^{\circ }$ for the wind direction transformation (Table 1). Because local wind measurements were made for 78% of the full year simulation period, the transformed data was used only when local measurements were not available.

4. Hindcast Simulations

Windsurf was applied to the Oysterville, WA field site. The following sections describe the procedure for model parameter calibration for a one-year simulation period and present the results of the calibrated case.

4.1. Windsurf Calibration

4.1.1. Model Setup

To calibrate Windsurf to Oysterville, a set of model simulations were initialized with the measured topographic and bathymetric profile from 4 August 2016 and run for a one-year period (until 9 August 2017) using the locally transformed SWLs, waves, and winds described above. A one-hour model coupling time step between the model cores was used, commensurate with the availability of oceanographic and meteorological forcing data and to adequately resolve processes on the intra-tide time scale. Within CDM the lower vertical limit of perennial vegetation was set to be 5.5 m based on local field observations (Figure 3b,c). CDM was initialized with spatially uniform vegetation (${\rho }_{veg}$ = 1) above this elevation.

A locally measured grab sample was used to inform the grain size (D50 = 0.25 mm, D90 = 0.335 mm) input to Windsurf. In the absence of detailed spatio-temporal grain size data at the field site, the grain size distribution is assumed to be spatially constant in both XBeach and Aeolis.

The 1D cross-shore model grid has variable grid spacing between 20 m (offshore) and 2.5 m and encompasses the region from –9 m to landward of the foredune crest. It is important to note that during energetic wave conditions wave breaking occurs deeper than the –9 m depth, inducing sediment transport near the offshore boundary and imposing implications on the incoming infragravity wave field. It was not practical to extend the model grid to significantly deeper water depths for this application. Therefore, consistent with the bathymetry observations (Figure 4) it was assumed that at the cells near the offshore model boundary (−9 m) there is no net bed level change on the time scales of interest. Any XBeach simulated sediment losses in this region are assumed to be replenished by either cross-shore feeding from the shoreface and/or longshore transport gradients. No modifications are made to modify the incoming bound long wave field, but it is assumed that infragravity generation within the inner surf zone from active wave breaking is of first order importance over far-field infragravity generation mechanisms for this model application.

4.1.2. Calibration Procedure

Process-based numerical models, including each Windsurf model core, typically have several tunable, site-specific parameters representing unknown model physics or site-specific transport properties e.g., [68,95]. Following preliminary model core sensitivity testing, four specific model parameters were identified for further exploration: the short-wave-related asymmetry (facAs, XBeach) and skewness (facSk, XBeach), the aeolian sediment transport coefficient (Cb, Aeolis), and the vegetation friction coefficient ((m, CDM). A choice of one year for model calibration was specifically chosen to ensure that the model is not biased towards simulating only accretional or erosional processes, as Oysterville has significant seasonal morphological variability (Figure 5).

As has been previously noted, XBeach is particularly sensitive to choices of facAs and facSk e.g., [68,69,70]—often requiring separate parameter combinations for accretional and erosional conditions e.g., [53]. Preliminary model core testing (not shown) did not yield realistic behavior when using constant facAs and facSk for the full year period. Therefore, following the general methods of Pender and Karunarathna [53], different model parameter combinations for these two variables were calibrated for low energy (assumed to be when H${}_s$≤ 2 m for our purposes) and high energy (H${}_s$ > 2 m) conditions. When H${}_s$ is below this 2 m threshold, field observations show that intertidal sandbars are typically observed to move onshore at the field site [88]. Above this wave height threshold, intertidal sandbars are generally eroded at Oysterville [46].

To calibrate Windsurf to the field site, a set of one hundred year-long Windsurf simulations were completed using random combinations of m, Cb, facAs (low energy), facAs (high energy), facSk (low energy), and facSk (high energy) (see

Table 2. Listing of range of values tested and optimum parameter settings based on the calibration routine.

Model Parameter	Mode Core	Min Value Tested	Max Value Test	Calibrated Value
m	CDM	0.005	0.20	0.104
Cb	Aeolis	0.20	1.50	0.370
facSk(low energy)	XBeach	0.10	0.50	0.123
facSk(high energy)		0.00	0.15	0.057
facAs(low energy)		0.10	0.50	0.273
facAs(high energy)		0.00	0.15	0.044

for range of values). Each of the 100 simulations were identical other than these parameter settings. Beyond this, typical model parameter settings are used for all three model cores. Please note that for XBeach a morfac of 2 was used and the dilatancy option was turned on.

Differences between the measured and modeled $\Delta V_{subaerial}$ ($\Delta V_{beach}$ + $\Delta V_{dune}$) were calculated at 10 time intervals (9 monthly surveys and the final August 2017 survey) from all 100 Windsurf simulations. The locally calibrated (best-fit) simulation was defined as the model simulation with the lowest RMSE $\Delta V_{subaerial}$ as calculated from all 10 $\Delta V_{subaerial}$ error estimates. This calibration approach was designed with the overall aim of discerning the relative controls of marine and aeolian processes on coastal foredune evolution. While calibrating the model to the monthly volumetric change data ensures that the model simulates realistic coastal profile changes at Oysterville, this exercise does not imply that Windsurf is validated to explore all coastal behaviors or for other sites e.g., [96,97].

4.2. Calibrated Model Results

The results from the best-fit calibrated simulation for Oysterville, WA between August 2016 and August 2017 are shown in Figure 4. Qualitatively, the simulated profile matches the field measurements, with growth of the dune and modest net changes to the beach modeled after one year. However, the model does not accurately reproduce the changes to the subtidal sandbars (Figure 4). Both subtidal and intertidal sandbars are smoothed during the simulations by XBeach (and potentially the other model cores).

In response to low wave energy conditions in summer, Windsurf simulates beach progradation by up to 44 m at the MHW shoreline (Figure 5). The widest beach state is modeled to occur in early November. Thereafter, the simulation shows gradual retreat of the MHW shoreline until May 2017. The shoreline remains relatively stable (±5 m) from May until the end of the simulation in August 2017. The beach volume changes track the same general temporal behavior as $\Delta X_{MHW}$. A total of 17.7 m${}^3$/m/year of dune growth is simulated compared to the field measured $\Delta V_{dune}$ of 16.7 ± 5 m${}^3$/m over the year.

Outputs from each model time step are used to investigate the relative contributions of marine and aeolian processes to coastal profile change and to explore the timing of depositional and erosional processes. The outputs are binned into 0.5 m vertical increments as shown in Figure 7, with these bins corresponding to the vertical elevations of the initial (4 August 2016) profile. In vertical cross-section, the numerical model results show reasonable agreement with the measured data over the annual scale (Figure 7a). Volumetric changes in the inner nearshore, between the 0 m and 1 m contours, are erosional in both the model and the observations. The model simulates more erosion around the 3 m contour and more deposition around the 4 m contour than revealed by the observations, although the model successfully reproduces the deposition of sediment between the 5.5 m and 7 m contours (Figure 4 and Figure 7a).

Windsurf simulates that growth of the lower portion of the beach (∼1–2 m elevations) is largest in summer, in agreement with the observations (Figure 7b). The largest upper beach growth (∼2–4 m elevations) instead occurs in fall, while the largest beach losses occur in winter. Although there is relatively little ($\Delta V_{beach}$ <30 m${}^3$/m) net change to the beach compartment over the year, the seasonal fluctuations in beach volume growth are significantly larger than volumetric variability within the dune region. The largest volume gains to the dune occur in the fall (10.3 m${}^3$/m), followed by winter (4.4 m${}^3$/m), spring (2.4 m${}^3$/m), and summer (0.6 m${}^3$/m).

As Windsurf computes dz related to marine transport by XBeach and aeolian transport by Aeolis at each coupling time step, it can be determined which volumetric changes across the profile can be attributed to marine versus aeolian processes, as shown in Figure 7c. The Windsurf simulations suggest that there are both marine (10.1 m${}^3$/m/year) and aeolian (7.6 m${}^3$/m/year) contributions to dune growth above the dune toe. Most of the marine-driven deposition to the dune is simulated near the 4 m contour—the approximate dune toe for the region. Above the 6 m contour, where swash did not frequently reach, all the simulated deposition (3 m${}^3$/m/year) is from aeolian processes.

These results suggest that marine-driven dune accumulation at 4 m or above occurs primarily during high SWL events and/or during periods of large wave energy (elevated H${}_s$ and/or T${}_p$; Figure 8a–c,f–h). Although most environmental conditions show positive $\Delta V_{dune}$, small (<0.1 m${}^3$/m) negative $\Delta V_{dune}$ are occasionally simulated (e.g., T${}_p$ of 19 s; Figure 8c) indicating that marine processes also drive dune erosion at this site under some forcing conditions. There are fewer environmental limiters on aeolian-driven dune growth (Figure 8), although the largest contribution to dune growth occurs from oblique or alongshore winds from the south (Figure 8e).

5. Discussion

5.1. Insights into Physical Processes Controlling Dune Evolution

This manuscript details the first application of the Windsurf modeling framework, demonstrating the ability for a coupled process-based numerical model to simulate exchanges between the nearshore, beach, and dune portions of the coastal profile. Windsurf predicts realistic morphologic evolution across the land-water divide with relatively limited model tuning. Despite the model framework not resolving the exact details of morphological evolution at the field site (e.g., sandbar dynamics), the ability to simulate beach accretion and erosion, wave-driven dune accretion, and wind-driven dune growth within a single model framework indicates that the model cores are resolving many of the dominant transport processes relevant for Oysterville. Therefore, model results are further explored to gain additional insights into the mechanisms contributing to beach and dune growth on sub-annual scale.

5.1.1. Aeolian Controls on Dune Growth

Although high wind velocities have the highest aeolian transport potential e.g., [38], many of the highest wind events co-occur with energetic waves and high SWLs (Figure 6). The combination of high H${}_s$ and high SWL results in TWL in the collision regime numerous times per year at Oysterville. During periods when the TWL is near the dune toe, there is effectively no source area for aeolian transport – limiting aeolian contributions to dune growth under these forcing conditions. Dune growth instead occurs predominantly from aggregated transport under moderate wind conditions (6–12 m/s) (Figure 8d,i). Intermittent dune growth occurs throughout the year under these frequent moderate wind conditions (Figure 5c) which are above the threshold velocity for saltation (3–4 m/s based on the local grain size). The largest total aeolian-driven dune growth is simulated in the winter (Figure 5c and Figure 7b). This is consistent with observations of the timing of maximum dune growth at Oysterville, WA as noted by Cohn et al. [46]. The model results also agree with field observations that documented limited dune growth occurring during the summer period (Figure 5c and Figure 7b) when the beach is wide (high sediment supply) (Figure 5a,b) but wind velocities are typically low (Figure 6k).

The largest modeled source area for aeolian sediment transport is between the 3 and 4 m contours, although some sediment is sourced from lower elevations during lower tidal stages (Figure 7c). While this primary source region is located on the upper portion of the beach compartment (where the beach is defined between 1 m to 4 m), it coincides with the approximate mean model-predicted TWL elevation during the fall and winter seasons (not shown). Thus, this modeled region of active aeolian transport generally supports previous findings that the upper intertidal zone is an important source of sediment for dune growth e.g., [41].

Wind direction plays a role in controlling the apparent fetch length, which is thought to be an important factor for aeolian transport to dunes e.g., [35]. Windsurf simulations suggest that, in Oysterville, winds from the south and south-west provide the largest contribution to total dune growth (Figure 8e), in part because these are among the most common wind directions at the field site (Figure 8j). In Figure 9a, it is shown that the largest hourly dune growth rates in the simulation occur when the wind speeds are highest, mostly independent of the wind direction. That is, similar $\Delta V_{dune}$ are simulated in the case of shore-normal or oblique winds for a given wind speed. There is some scatter in the hourly $\Delta V_{dune}$ predictions (Figure 9a), reflecting in part the cosine effect e.g., [31] and/or unsaturated transport resulting from sediment supply limitations. However, when considering higher portions of the dune, shore-normal winds provide a larger contribution to dune growth for a given wind speed (Figure 9b). This is because, in the case of shore-normal winds, saltation can extend farther horizontally past the vegetation line because of a larger cross-shore component of the wind velocity and increased wind velocities on the dune face associated with flow acceleration over the steeper apparent dune slope.

5.1.2. Marine Controls on Dune Growth

While the paucity of observed dune erosion on some dissipative beaches has been hypothesized to be from synchronization of nearshore-beach-dune exchanges [26], which were thought to mask signs of erosion, modeled dune erosion at Oysterville, WA is infrequent because TWLs in the collision regime are not always erosional [51]. Windsurf predicts that numerous wave-driven events contribute to dune accretion at Oysterville during the 2016–2017 time frame (Figure 5, Figure 7, and Figure 8), with these events typically occurring during high SWL cases and/or during energetic wave periods (Figure 8). The simulations show that there were 3 h, all in October 2016, which resulted in at least 0.5 m${}^3$/m of marine-driven dune growth. Of note, these marine-contributions to dune growth are about an order of magnitude larger than the highest modeled hourly dune growth rates from aeolian sediment transport (Figure 9). Seven additional hours were simulated throughout fall and winter when marine-derived dune growth of at least 0.1 m${}^3$/m occurred. In Cohn et al. [46], it was estimated that up to ∼ 7 m${}^3$/m/year of dune growth at Oysterville, WA during the 2016–2017 period was driven by swash processes, similar in magnitude to the modeled 10.1 m${}^3$/m/year. It is important to note that the choice of the dune toe definition used for calculating volume changes has important implications for determining the relative contributions of marine and aeolian processes to dune growth. Traditional concavity-based metrics for defining the dune toe do not work well for low-gradient beach systems such as Oysterville. Therefore in Cohn et al. [46] and this study the dune toe is taken to be at 4 m based on local geomorphic and ecological characteristics. An alternative definition of the beach-dune interface would correspondingly alter the contributions of aeolian and marine processes to total foredune growth reported in this study. However, as this work demonstrates, the boundary between the beach and dune is not as simple as typically characterized in geomorphic studies.

While a positive $\Delta V_{dune}$ is simulated most of the time when TWLs are in the collision regime, there are also some oceanographic conditions which instead produce wave-driven dune erosion (Figure 5c and Figure 8b,c). Therefore, the assertion that marine-related erosion can be masked by aeolian processes on dissipative coasts e.g., [26] is supported by the model under some conditions. This finding complicates attempts to decouple marine and aeolian contributions to dune sediment supply from field observations alone e.g., [46].

5.2. Sensitivity of Dune Growth to Environmental Perturbations

The ability to hindcast profile evolution at the field site provides some confidence that Windsurf is resolving the dominant morphodynamic processes at Oysterville. However, process-based models also provide the ability to explore environmental conditions outside of what was observed. Here additional simulations are completed to gain insight into the relative role of oceanographic and meteorological processes on dune evolution at the annual scale. To accomplish this, the year-long environmental time series of u${}_w$, SWL, and H${}_s$, were independently modified from the baseline (calibrated parameter) case, in which 17.7 m${}^3$/m of annual dune growth was simulated, to explore how the system would respond to environmental perturbations. First, these parameters are increased or decreased by fixed quantities. For example, one simulation was completed in which 0.5 m is added to every H${}_s$ value input to the XBeach/Windsurf time series. While a fixed change in the SWL can be viewed as a proxy for mean sea level changes, there is no expectation for a uniform/constant increase in H${}_s$ or u${}_w$ even under climate change scenarios. However, these hypothetical exploratory simulations are devised to understand how these different environmental forcing factors may alter subsequent geomorphic evolution. Here constant additions or subtractions covering the range of −1.5 to 1.5 m (or m/s in the case of u${}_w$), with increments of 0.5 m (m/s), were made to the environmental time series. Only one variable (either u${}_w$, SWL, or H${}_s$) is altered in each simulation; all other attributes of the baseline case are maintained.

The results of this exercise are shown in Figure 10 and, as expected, indicate that decreasing wind speeds result in decreased volume gains to the dune and increasing wind speed drives increased $\Delta V_{dune}$. Because aeolian sediment transport scales to a third order power (Equation (2)) the relationship between u${}_w$ and $\Delta V_{dune}$ is nonlinear and therefore increases in wind speed (including higher values than those presented in Figure 10) result in an exponential type dune growth curve (not shown).

Increases in the SWL result in decreased dune growth relative to the baseline case (Figure 10) due to active wave breaking closer to the dune face and higher TWLs. The model simulates that even small increases in the SWL will actively destroy the dune ($\Delta V_{dune}$ relative to the baseline <−17.7 m${}^3$/m in Figure 10). However, decreases in the SWL also decrease the frequency that TWLs reach the dune toe, effectively cutting off any wave-induced sediment supply to the dune. Therefore, all cases with SWL that were lower than the hindcast time series had decreased $\Delta V_{dune}$ relative to the baseline case. The observed SWL time series resulted in the largest dune growth of any simulated SWL cases.

Although increases in wave height are often associated with enhanced beach and dune erosion, the model simulations here suggest that increases in wave height may contribute positively to additional dune growth at Oysterville. This is consistent with model simulations in Cohn et al. [51] that found that as long as the dynamic still water level (DSWL), which is defined as the combination of the SWL and wave setup, falls below the dune toe that collisional wave impacts may contribute positively to growth of the lower portion of the dune at Oysterville. In this study, the exploratory Windsurf simulations suggest that uniform wave height increases of up to 1.5 m may encourage increased sediment delivery by marine processes to the base of the dune. However, it is expected that for more extreme wave height scenarios, increased frequency of the DSWL above the dune toe will instead promoted dune erosion. Assessing these extreme wave characteristic changes would require modifications to the Windsurf grid, which only extended to −9 m in this study, and therefore additional wave height simulations were not completed.

5.3. Windsurf Improvements and Future Applications

While Windsurf demonstrates the ability to simulate some realistic profile evolution in this application, the need for parameter calibration and the lack of a perfect match to data likely represents the effects of unknown physics, unincorporated physical processes, and/or oversimplified dynamics. Incorporating new processes or improving existing algorithms has the potential to increase the reliability and applicability of these tools in other settings. However, as both marine and aeolian sediment transport processes are highly complex, an improved physical understanding of coastal morphodynamic processes gained from field efforts is also essential to further improving these numerical capabilities e.g., [30]. Below some of the current key limitations of the Windsurf framework are highlighted.

5.3.1. Model Processes and Parameterizations

Of the three model cores in Windsurf, XBeach has the largest number of configurable parameters and is also generally the most sensitive to those model parameter choices. The fact that a separate high and low energy facAs and facSk were required for realistic model behavior implies that new formulations for these parameters that are more physically based may be necessary to make reliable forecasts of coastal progradation using XBeach in surfbeat mode. Recent improvements to XBeach add a non-hydrostatic correction term [98], which enables short waves to be modeled directly, and may provide a way forward. In non-hydrostatic mode, the nonlinearity of the wave shape is calculated directly, avoiding the need for facSk and facAs as these parameters are inherently included within the flow field. Therefore, using the wave resolving mode of XBeach could enable less site-specific model tuning. However, the non-hydrostatic model requires finer grid resolution which, for most applications, would substantially increase model run times. Furthermore, new sediment transport algorithms relevant for the non-hydrostatic model first need to be developed in the open-source code.

Although not explored in detail in this study, the algorithms parameterizing spatio-temporal vegetation density in CDM are currently simplified (a single, spatially uniform maximum ${\rho }_{veg}$ is currently implemented). However, the vegetation plays a major role in the sand trapping capacity of sediment and resulting geomorphic shapes e.g., [21,99]. Sensitivity to dune grass parameters (e.g., height and density) were not explored in this study other than altering the m parameter during model calibration. As these ecological properties vary in both space and time and not all coastal systems are characterized by a single grass species, more robust approaches could be developed for use within the model. Therefore, while the currently implemented formulations may be suitable for understanding first order behavior, the addition of more robust ecological algorithms will provide added capability for simulating these dune growth behaviors for more applied cases e.g., [80].

There are additional moisture-related limiters to aeolian sediment transport in the coastal zone that were not assessed for this study. Much of the functionality to incorporate both groundwater and precipitation related effects are already coded into XBeach and Aeolis. Including these additional supply limiting factors may be important for some locations, particularly those characterized by high rainfall such as at Oysterville, WA.

On coastal systems characterized by storm-induced scarping and sharp topographic changes, the presently incorporated Weng et al. [81] solution for wind flow will have limitations in its applicability. For more complex settings a steady state wind solver may not be suitable. In these cases, a Navier-Stokes-based numerical solver could be implemented, as computational fluid dynamics methods have been shown to be successful at modeling complex flow patterns over a range of dune shapes e.g., [100,101]. The framework of Windsurf is modular by nature which allows for new models to be implemented and tested with relative ease.

5.3.2. External Sediment Supply

While the assumption of local mass conservation for the time scales investigated for this study are appropriate, when extending simulations to longer time frames incorporating external sediment supply inputs is imperative e.g., [5,102]. The definition of the offshore boundary in Windsurf currently allows for new sediment to be added to the local sediment budget via the replenishment of eroded near-boundary sediments. However, a more robust procedure for introducing both cross-shore shoreface feeding and longshore transport gradients would add value. A recent implementation of longshore gradients in transport into XBeach (lsgrad) provides one mechanism to do this for 1D applications based on an imposed length-scale for longshore transport gradients. Although not tested for this application, using this feature in Windsurf would require no changes to the model framework. Including these 1D imposed longshore transport gradients, or extending Windsurf to 2D area applications, would better represent coastal systems, including Oysterville, WA, where longshore transport gradients are important.

5.3.3. Computational Efficiency

The present Windsurf application is run in an offline approach through MATLAB. A python version of Windsurf is also in development which implements the basic model interface [103] and enables models to be run in an on-line approach, with boundary conditions updated on the fly. The latter approach significantly decreases model run times by eliminating computational overhead associated with frequent writing of model input and outputs to the disk. For the present application, completing one-year simulations takes approximately 2.5 days using 4 processors on a Dell C6100 Intel Xeon server running the CentOS 7 operating system. Reducing run times, even fractionally, may make simulating longer term (>years, decades) simulations more practical.

5.3.4. Further Testing of Framework Limitations

These first proof of concept Windsurf simulations at Oysterville, WA suggest that a coupled approach for assessing the co-evolution of the nearshore-beach-dune system provides value over independent modeling of discrete morphological units. As this is the first application of the modeling framework, the limits of Windsurf at simulating realistic behavior in other morphodynamic settings and for other environmental cases have not yet been demonstrated. Ideally, Windsurf is flexible enough to explore dune erosion, recovery, and growth across the continuum of reflective to dissipative beaches. However, to understand model limitations and extend model capabilities, the tool needs to first be applied in other data rich settings. Similarly, additional future work should aim to characterize model sensitivity to input choices, particularly for spatially variable properties such as sediment grain size and ecological properties.

6. Conclusions

A new coupled numerical modeling framework, Windsurf, can simulate the co-evolution of the nearshore, beach, and dunes. Consistent with detailed morphologic measurements at a progradational field site, the model successfully simulates seasonal cycles of beach growth in summer, shoreline recession in winter, and net dune growth annually. The modeling results are also consistent with findings from field observations, both of which suggest that both marine and aeolian processes directly contribute to the growth of coastal foredunes on this high energy, dissipative coastline. The importance of, and challenges involved in, distinguishing between marine and aeolian controls on dune evolution revealed by field observations and model simulations highlight the need for further exploration and development of coupled tools bridging the land-sea interface for coastal landscape evolution. While significant efforts are required to further test the capabilities and limitations of Windsurf and its model cores, the Windsurf framework provides a new platform to explore complex interactions between the subaqueous and subaerial zones of the coastal profile for a variety of exploratory and applied applications.