Exploring marine and aeolian controls on coastal foredune growth using a coupled numerical model

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.


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], 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 Sections 5 and 6, respectively.

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.

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) (Figures 1 and 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. Schematic representing the general Windsurf model framework which includes three standalone numerical models (XBeach, CDM, Aeolis) that are coupled offline through a back-end MATLAB interface. Major processes resolved within each model core and the outputs exchanged by the coupler are shown. Black arrows show the order of model processes and model core boundary condition exchanges. Morphology and total water level information is also exchanged between XBeach and Aeolis, as denoted by the dotted grey arrow.

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

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 (δτ) 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 (ρ veg ; [22]). The effective bed shear stress (τ s ) in the presence of vegetation, which is calculated by CDM, becomes: where τ o is the bed shear stress in the absence of variable topography or vegetation, m is the vegetation friction coefficient, σ is the ratio of the plant basal to frontal area, and β 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 δτ 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].

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] : where ρ 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 ) : where ρ 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].

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.

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

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

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 ∆X MHW over the full year. Volume changes in the beach and dune are also calculated at this monthly time scale (Figure 5b,c). ∆V beach has the same seasonal behavior as ∆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.

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. Locally measured (red) and transformed (grey) still water levels (a), significant wave heights (b), peak wave periods (c), wave direction (d), wind speed (e), and wind direction (f) for a 1 month period in summer 2016. The continuous, year-long time series of environmental parameters between summer 2016 and summer 2017, which includes the one-month period (shaded), for input to Windsurf is shown in panels (g-l) using the transformed variables, where applicable. Please note that, as per model core conventions, 270 • is shore-normal for D wave (nautical convention) and 0 • is shore-normal (onshore directed) for D wind .

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 (∆µ) and RMSEs of 0.12 m, 1.7 s, and 8.0 • for H s , T p , and D wave , respectively, relative to the AWAC data (Table 1; Figure 6b-d).

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 ∆µ (−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 • of the measurements, there is a fairly high overall RMSE value of 25.5 • 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.

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.

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 (ρ 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.

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 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 ∆V subaerial (∆V beach + ∆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 ∆V subaerial as calculated from all 10 ∆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].

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 ∆X MHW . A total of 17.7 m 3 /m/year of dune growth is simulated compared to the field measured ∆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 (Figures 4 and 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 (∆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 ∆V dune , small (<0.1 m 3 /m) negative ∆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).

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.

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 (Figures 5c and 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 (Figures 5c and 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 ∆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 ∆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.

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 (Figures 5, 7, and 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 ∆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 (Figures 5c and 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].

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 ∆V dune . Because aeolian sediment transport scales to a third order power (Equation (2)) the relationship between u w and ∆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 (∆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 ∆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.

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.

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

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.

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.

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.

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.