Conceptualizing Aeolian Sediment Transport in a Cellular Automata Model to Simulate the Bio-Geomorphological Evolution of Beach–Dune Systems

: Understanding the dynamics of beach–dune systems is crucial for effective coastal management. The cellular automata model DuBeVeg provides a powerful tool for simulating and understanding the bio-geomorphological evolution of these systems, capturing key interactions of aeolian, hydro-, and vegetation dynamics in a simpliﬁed manner. In this study, we present an alternative representation of the aeolian transport component in DuBeVeg, aiming to better capture the saltation transport mode that prevails on beaches. This new representation is compared with the original aeolian transport representation in DuBeVeg, which is inspired by ripple migration. For three beach width scenarios, we considered the effects of the different aeolian transport representations on the predicted foredune morphology after 50 years, as well as the spatio-temporal evolution of the beach–dune system leading to that morphologic state. The saltation transport representation resulted in a more realistic simulation of the seaward expansion of the foredune compared with the original representation, particularly in scenarios with wide and prograding beaches. The new representation also more accurately captured the amplitude of aeolian bedforms emerging across the beach. These ﬁndings highlight the importance of selecting the representative transport mode when simulating the transient bio-geomorphological evolution of beach–dune systems.


Introduction
Sandy coasts provide many important ecosystem services such as coastal protection, habitat support, and opportunities for tourism and recreation [1][2][3][4][5]. Covering more than one-third of the world's coastlines [6], sandy coasts are highly vulnerable to coastal hazards such as erosion and sea-level rise [6][7][8][9][10][11][12][13]. The beach and the adjacent foredune (i.e., the beach-dune system) serve as the primary defense against coastal hazards, with their natural dynamics allowing the system to adapt to changing conditions, including sea-level rise [14][15][16]. The foredune is particularly crucial, acting as a natural barrier that prevents water from reaching inland areas while also serving as a natural source of sand for the beach [17,18]. Understanding the development of the beach-dune system is, therefore, essential for predicting its continued contribution to coastal safety and other ecosystem services under changing environmental conditions.
The beach-dune system is a sand-sharing system, where sediment is exchanged between the nearshore zone, the beach, and the foredune. This exchange is governed by natural processes, i.e., aeolian, hydrodynamic, and vegetation processes. Nearshore sediment can be transported onto the beach by the combined effects of waves and tides. Typically, these accretionary processes can reach as high up the beach as the spring high tide level, with an additional height for wave runup. When this lower part of the beach falls dry, the newly deposited sand can be transported higher up the beach and into the foredune by onshore winds. Vegetation can enhance dune formation by capturing sediment and stabilizing the topography. Under storm conditions, however, elevated water levels and energetic waves may erode significant volumes of sand from the foredune and transport this sand seaward [19][20][21]. Consequently, the beach-dune system evolves at different spatio-temporal scales ranging from the movement of sand grains over seconds to the advance/retreat of foredunes over decades, with each scale playing an important role in shaping the coastal landscape [22].
Self-organization is a useful concept for understanding the morphodynamic evolution of sandy systems. It refers to the emergence of structure or order in a complex system through its internal dynamics and feedback mechanisms, where feedback loops between the components of the system lead to the emergence of regular features or patterns [23][24][25][26][27][28][29][30][31][32].
In the context of beach-dune systems, there is ample opportunity for self-organization to emerge. Aeolian morphodynamics, wave-driven morphodynamics, and evolving vegetation all interact, often in a nonlinear way, making the dynamics of this sandy system, as a whole, highly complex [23,33,34].
Cellular automata (CA) modeling is a powerful tool for simulating the dynamics of complex nonlinear systems that exhibit self-organizing behavior and is widely used in geomorphology [35][36][37][38][39][40][41][42][43][44]. The Werner model is a CA model for aeolian dune formation [33], which was later extended by Nield and Baas [45] in the DECAL model by including vegetation dynamics and interactions to reproduce the dune patterns that emerge in coastal landscapes. The CA algorithm simulates sediment transport through the stochastic movement of slabs of sand, which are picked up by the wind and moved through the domain in the downwind direction by a certain jumping length (L). After being moved over a distance L, a slab can either be deposited or continue jumping downwind by another distance L depending on the probability of deposition. Positive feedback mechanisms allow for the self-organization of the slabs, resulting in the emergence of bedforms. These mechanisms include the deposition of slabs at the sheltered leeside of dunes and an avalanching mechanism that limits dune slopes to the angle of repose.
To simulate the behavior of coastal dunes in conjunction with the dynamics of the beach, the DuBeVeg (Dune-Beach-Vegetation) model was developed [46] by expanding DECAL to include the effects of waves, tides, and storm surges that bring sediments to the beach and can erode the beach profile. The DuBeVeg model includes the effects of aeolian transport, hydrodynamic erosion and accretion, groundwater levels, and vegetation growth. The DuBeVeg model considers the stochastic nature of the aeolian transport and wave erosion/accretion mechanisms, as well as the emergence of new vegetation in the beach-dune system and the limiting effect of groundwater on aeolian transport (i.e., wet sediments cannot be transported by the wind). The chance of the erosion and deposition of sand slabs caused by wind and waves is stochastically determined in each iteration. This allows the model to account for the temporal variability of these dynamics. The model output parameters (e.g., dune volume) should then be interpreted as one possible outcome of an uncertain reality, and several repetitions of the same model run should be considered to explore the bandwidth of the model's predictions [46,47]. The DuBeVeg model has been used successfully to study the effects of sea-level rise on dune evolution [46], the effects of beach width on dune morphology [47], and the effects of groundwater levels on dune morphology [48]. The DuBeVeg model was also used to study the effects of buildings on the development of aeolian bedforms at the beach [49]. While the DuBeVeg model has proven effective in simulating the emergent biogeomorphological patterns of sandy beach-dune systems, it was primarily developed to focus on the final morphological patterns that emerge on a decadal timescale rather than the rates of change or the temporal evolution toward those final states. Consequently, previous studies did not prioritize the accurate reproduction of the temporal evolution of sediment fluxes across different zones of the beach. Previous implementations of the model consider the migration of ripples the main process responsible for the movement of sand across the beach [45]. However, in reality, the saltation of sand grains is the predominant transport mode of sand at the beach [22,23]. Hence, it is hypothesized that a representation of aeolian transport that considers saltation as the main transport mode could better reproduce the temporal evolution of the sediment flux and the transient bedform dynamics across the beach. Accordingly, an adjustment to the representation of the aeolian transport mode in DuBeVeg is suggested to expand the model for application to a wider range of scenarios where saltation transport is dominant, particularly for anthropogenic sandy shores that typically include higher and wider beaches than would be naturally present and that typically show transient behavior with respect to foredune formation.
The aim of this study is to expand the cellular automata modeling of sandy beach-dune systems, with a particular focus on the aeolian dynamics occurring on the beach surface and the transient emergence of new foredunes. We propose a new representation of aeolian transport in the DuBeVeg model that is more in line with the saltation transport mode. We investigate how this new representation affects the emergence of bedforms at the beach and the temporal evolution of emerging foredune(s). We also examine which beach-dune configurations the new representation of aeolian transport in DuBeVeg is important to. Section 2 of this paper presents the DuBeVeg model and the adaptations that were made for representing the saltation transport mode, along with the model and parameter settings deployed in this study. Three scenarios are introduced to explore the predicted beach-dune development for the default and adapted aeolian transport representations. Section 3 then presents the simulation results of the default and adapted DuBeVeg representations for the three scenarios, focusing on the predicted foredune morphologies and the temporal beach-dune evolution. The impacts of the new representation of aeolian transport in the DuBeVeg model are discussed in Section 4, followed by conclusions in Section 5.

DuBeVeg Model
The Dune-Beach-Vegetation CA model (DuBeVeg) simulates the bio-morphological development of sandy beach-dune systems through the combined effects of aeolian transport, hydrodynamic erosion and accretion, groundwater levels, and vegetation development. A concise description of the original DuBeVeg model is provided below, complemented by a more elaborate description of the proposed model modifications. A full description of the original DuBeVeg model can be found in Keijsers [46].

Aeolian Module
Aeolian sand transport is the core process of the DuBeVeg model. Aeolian transport is represented by the stochastic reorganization of "sand slab" units, i.e., the model's representation of a standard volume of sand for a single grid cell, resulting in a default slab height for the uniform grid across the model domain. At each model iteration, a slab has a chance to be picked up and moved through the domain in the downwind direction by a certain jumping length, L. After being moved over a distance L, a slab can either be deposited or continue jumping over another distance L downwind until it is deposited (Figure 1a). This movement of the slabs over the grid cell domain simulates the transport of the sand by the wind.
A physical basis underlies the stochastic picking up and deposition of the slabs. Each sand slab has a basic probability of erosion (p e 0 ) and deposition (p d 0 ). These probabilities are a calibration parameter of the model and mainly relate to the mean grain size of the beach in combination with the wind climate. The basic probabilities of erosion and deposition are modified by the moisture content of the sediment and the presence of vegetation. The transport-reducing effect of moisture content is represented in the model by a groundwater surface (GWS), limiting the depth of erosion caused by wind in the dry sand above this surface. The effect of vegetation on the reduction in wind velocity and, therefore, the transport capacity of the wind is introduced in the model with the vegetation effectiveness parameter (ρ). This parameter represents the vegetation cover (as a percentage) for each individual cell. With a ρ equal to 0, the cell is a completely bare surface, while with ρ equal to 1, the cell is completely covered by vegetation. Thus, the total probability of erosion (P e ) and deposition (P d ) for each slab at each iteration are determined as follows: A physical basis underlies the stochastic picking up and deposition of the slabs. Each sand slab has a basic probability of erosion (p e 0 ) and deposition ( p d 0 ). These probabilities are a calibration parameter of the model and mainly relate to the mean grain size of the beach in combination with the wind climate. The basic probabilities of erosion and deposition are modified by the moisture content of the sediment and the presence of vegetation. The transport-reducing effect of moisture content is represented in the model by a groundwater surface (GWS), limiting the depth of erosion caused by wind in the dry sand above this surface. The effect of vegetation on the reduction in wind velocity and, therefore, the transport capacity of the wind is introduced in the model with the vegetation effectiveness parameter (ρ). This parameter represents the vegetation cover (as a percentage) for each individual cell. With a ρ equal to 0, the cell is a completely bare surface, while with ρ equal to 1, the cell is completely covered by vegetation. Thus, the total probability of erosion (P e ) and deposition (P d ) for each slab at each iteration are determined as follows: Two restrictions are considered for the movement of slabs. First, in the shadow zone, which is defined by a 15° angle at the leeward side of any protruding slab, slabs are forced to deposit and cannot be picked up. This mechanism reflects, in a very simple manner, the deceleration of wind behind any obstacle, such as dunes. Second, an avalanching mechanism is considered, only allowing slopes to be equal to or gentler than the angle of repose. This angle is set to 30° for bare slabs and 35° for vegetated slabs.
To apply the model for a specific wind climate, the transport rates obtained with the aeolian transport mechanism described above must agree with representative values from field studies. Nield and Baas [45] calculated the potential aeolian transport per meter alongshore, Q (m 3 /m/yr), in terms of the aeolian module parameters: where h represents the slab height and the number of iterations per year. The default Two restrictions are considered for the movement of slabs. First, in the shadow zone, which is defined by a 15 • angle at the leeward side of any protruding slab, slabs are forced to deposit and cannot be picked up. This mechanism reflects, in a very simple manner, the deceleration of wind behind any obstacle, such as dunes. Second, an avalanching mechanism is considered, only allowing slopes to be equal to or gentler than the angle of repose. This angle is set to 30 • for bare slabs and 35 • for vegetated slabs.
To apply the model for a specific wind climate, the transport rates obtained with the aeolian transport mechanism described above must agree with representative values from field studies. Nield and Baas [45] calculated the potential aeolian transport per meter alongshore, Q (m 3 /m/yr), in terms of the aeolian module parameters: where h represents the slab height and n the number of iterations per year. The default time-step of the aeolian module in DuBeVeg represents 7.3 days and, hence, n = 50 iterations per year. Equation (3) shows how the values of the aeolian model parameters impose the transport rates simulated with the model and, therefore, should represent the transport capacity of the wind climate of the beach-dune system under consideration. Equation (3) does not consider the limiting effect of groundwater or vegetation.

Hydrodynamic Module
The erosion and accretion of sand from the beach-dune system caused by the sea are modeled in the hydrodynamic module of DuBeVeg. Elevated seawater levels and energetic waves may erode sand from the beach and dunes, transporting this sand seaward. Part of the eroded sand will be lost to the open sea, and part will remain available on the nearshore to be transported back to the beach during calmer conditions, restoring the eroded beach profile. This loss or gain of sediment from the beach and foredune is modeled in the hydrodynamic module by the stochastic removal or appearance of slabs in the part of the beach-dune profile affected by the sea.
For each iteration of the hydrodynamic module, the maximum extent of the cross-shore profile affected by the hydrodynamic effects is determined from an input time-series of water levels and a statistically correlated wave height to estimate wave runup. The wave runup is estimated from the water levels in the same manner as presented in Galiforni-Silva et al. [47] and Keijsers et al. [46], based on the empirical relationship in Stockdon et al. [50] and is added to the input water levels. Each iteration of the hydrodynamic module is set to represent the total hydrodynamic effects during one full spring-neap cycle. Thus, the input water level is the maximum water level over 14.6 days, and 25 iterations of the hydrodynamic module are considered for a full year of model simulations.
For the part of the beach-dune profile affected by the sea, a probability (P hydro ) is calculated for the height of each grid cell to return to a predefined reference surface: where R V is a constant value that is multiplied by the vegetation cover ρ at each cell and represents the resistance of a vegetated surface to wave erosion; W E is a constant value that represents the wave erosive strength; W diss represents a wave dissipation factor that is multiplied by a function, D, which relates to the inverse water depth; and P inun represents a constant probability of erosion/accretion due to inundation without the effect of waves. For each grid cell below the input water level, P hydro then quantifies the probability of slabs above the reference surface being removed and new slabs appearing on top of those slabs below the reference surface. The reference surface can be interpreted as an equilibrium surface (or equilibrium profile), i.e., the surface shape that the beach will tend toward, in the long term, as a result of the balance of erosive versus accretive forces [51]. At the same time, beach groundwater has been found to be coherent with the erosive and accretionary trends of beach profiles [52]. Hence, the groundwater level in the model is defined as a surface proportional to the reference surface. Consequently, the reference surface has to agree with the observed or estimated equilibrium beach slopes for the location of the study and with the observed general trends (e.g., slopes) of the groundwater level. The reference surface used in previous implementations of DuBeVeg was a simplified representation of the initial beach surface, which was obtained by transforming the beach surface into a smoothly sloping plane excluding the dunes.
The groundwater surface (GWS) is defined relative to the reference surface for all cells. The relative elevation of the GWS with respect to the reference surface is defined by a constant groundwater depth factor (G) that can range between 0 and 1 (z GWS = G · z reference surface ). If G is equal to 0, the GWS is a horizontal plane at mean sea level (0 elevation is considered the reference mean sea level). If G is equal to 1, the GWS is equal to the reference surface. In cells where vegetation is present (ρ > 0), the growth and decay of the vegetation are determined by the annual sediment balance of each cell. Growth functions [53] determine a new value of ρ based on simulated annual elevation differences in each cell. The growth functions are segmented linear functions that establish the limits of erosion/accretion that vegetation can withstand, as well as the optimal erosion/accretion conditions for vegetation growth. Two types of vegetation are considered in the DuBeVeg model: a pioneer species and a stabilizer species. The pioneer species grow in a range of accretion rates of 0.2 m/yr up to 2.2 m/yr with optimal growth for 0.6 m/yr accretion. Higher or lower accretion/erosion rates cause a decrease in the pioneer species' vegetation. The stabilizer species' optimal growth occurs when there is no erosion/accretion in the cell, and its growth is defined by a lower range of accretion/erosion rates, ranging from −0.65 m/yr erosion to 0.2 m/yr accretion.
The vegetation module also stochastically accounts for the lateral expansion of vegetation. When a bare cell (ρ = 0) is adjacent to a vegetated cell, it has a probability of becoming vegetated and commencing vegetation growth based on the above growth functions. The probability of lateral expansion is uniform for all cells neighboring a vegetated cell, regardless of whether they are partially or fully vegetated. The model also considers the possibility of previously bare cells becoming vegetated even if no neighboring cell is vegetated. This dynamic accounts for the dispersion of seeds and rhizomes by the wind and hydrodynamics. Analogous to the lateral expansion, the settlement of vegetation is modeled probabilistically, with an equal probability assigned to each bare cell. Once a previously bare cell becomes vegetated, its vegetation starts to develop based on the annual sedimentation balance and the growth functions defined above.

DuBeVeg Adaptation to Represent Saltation
Existing DuBeVeg model applications consider a default jumping length, L, of 1 grid cell. The aeolian transport mechanism represented by this jumping length unit can be interpreted as the movement of ripples of sand over the beach/dune surface [45]. The analogy between the movement of slabs and the propagation of ripples is based on the relatively large slab size (in the order of 1 m 2 ), which is significantly bigger than individual grains, and the jumping length is restricted to one cell size. Accordingly, Nield and Baas [45] determined the size of the slab by relating the volume of sand transport to a realistic transport celerity of ripples and mentioned that a longer jumping distance of multiple slab sizes would be more appropriate for representing the saltation process.
Nield and Baas [45] related the aeolian transport mode in DECAL to ripple migration as they were aiming to apply the algorithm of aeolian transport to infinite and uniform coastal landscapes, i.e., initially flat domains lacking any hydrodynamic impacts and without the characteristic cross-shore profiles of the beach-dune system. In a beach-dune system, however, the saltation of sand grains is the predominant transport mode [22,23]. This suggests that, in order to simulate the dynamics of the beach-dune system, the movement of slabs should represent the saltation process. As suggested by Nield and Baas [45], this would require an increase in the unit jumping length. Accordingly, the original algorithm by Werner [33] applied a jumping length of L = 5 and related the simulated dynamics to saltation transport. Initial DECAL implementations [54] also applied a jumping distance of L = 5. This use of larger distances for the jumping length was later discouraged by Nield and Baas [45], as they suggested it only dilated the pattern emergence (i.e., increased the spacing between the emerging bedforms) without affecting their morphology. Moreover, with greater jumping lengths, jumping slabs would ignore the presence of vegetation in intermediate cells, as slabs only had a chance to deposit in vegetated cells at distances equal to (in multiples of) the jumping length, L, in the downwind direction.
In this study, the probability of deposition, P d , is updated to allow a picked-up slab to move across a jumping length of 5 cells (L = 5) in order to explore the effects of saltation transport on the model simulations. Considering the concerns of Nield and Baas [45], we adapted the transport mechanism to be able to move slabs over multiple grid cells without ignoring vegetated cells that are traversed. The vegetation in any intermediate cell is considered a catching element, causing a jumping slab to have a chance to be deposited at any location along its jumping length (Figure 1b). A vegetation effectiveness limit parameter (ρ * ) is introduced to define the minimum required vegetation effectiveness for a cell along the jumping path to be considered vegetated. Accordingly, the probability of deposition (P d ) is updated for each slab that is transported, distinguishing between the cells the slab jumps across (L < 5) and the destination cell if the slab does not become trapped along its jumping path (L = 5): This way, if there is not enough vegetation in the cells downwind of a picked-up slab to initiate deposition (ρ < ρ * ), the slab moves with a jumping length of 5 until it is potentially deposited in its destination cell (with a probability, P d ; see Equation (6)). However, if the vegetation effectiveness in downwind cells in between the origin of the slab and the jumping length is greater than the threshold value (ρ ≥ ρ * ), the slab has a chance of being deposited in those cells (also with a probability, P d ; see Equation (5)).
This new model representation allows slabs to jump longer distances than just one cell, similar to the first implementations of Baas [54] and Werner [33], but, at the same time, it takes into account the vegetation as a catching element, which effectively reduces the larger jumping length by promoting deposition in any vegetated cells along the transport path.

Test Scenarios
This paper compares the newly adapted representation of the DuBeVeg model, which incorporates saltation transport, with the existing default representation of the model, which accounts for sand transport via ripple migration. We designed three simple scenarios to evaluate morphological beach-dune evolution predicted by the two transport representations. A simple initial morphology with a smooth beach and an existing foredune was used to simulate 50 years of bio-geomorphological development. The three idealized initial morphologies account for a narrow beach with a fixed shoreline position (100 m beach width), a wide beach with a fixed shoreline position (500 m beach width), and a prograding shoreline scenario with a beach that widens at a rate of 4 m/yr (from 100 m to 300 m) to represent a gradual increase in accommodation space for new dunes to form ( Figure 2). All three scenarios have simplified beach profiles, where the elevation increases linearly from mean sea level (MSL) at the shoreline to 3 m +MSL at the toe of the foredune. Below MSL, the profiles follow a Dean equilibrium profile [51]. The initial foredune is schematized as a simple dune ridge with a maximum elevation of 15 m +MSL and a width of 100 m (Figure 2). In the prograding scenario, the beach profile builds out from the elevation of a berm at 2 m +MSL down to a depth of closure at 8 m −MSL, i.e., the depth beyond which the profile is not affected by the longshore transport (Figure 2b), similar to the simplification of the cross-shore extension following longshore sediment transport in one-line models [55]. The dimensions of the narrow and wide beach scenarios are based All three scenarios have simplified beach profiles, where the elevation increases linearly from mean sea level (MSL) at the shoreline to 3 m +MSL at the toe of the foredune. Below MSL, the profiles follow a Dean equilibrium profile [51]. The initial foredune is schematized as a simple dune ridge with a maximum elevation of 15 m +MSL and a width of 100 m ( Figure 2). In the prograding scenario, the beach profile builds out from the elevation of a berm at 2 m +MSL down to a depth of closure at 8 m −MSL, i.e., the depth beyond which the profile is not affected by the longshore transport (Figure 2b), similar to the simplification of the cross-shore extension following longshore sediment transport in one-line models [55]. The dimensions of the narrow and wide beach scenarios are based on the observed beach-dune profiles on the Dutch coast. The prograding scenario starts with the narrow beach scenario but increases in width to a situation midway between the former two scenarios.
The applied model domain has a size of 400 × 1250 cells of 1 × 1 m 2 -corresponding to 400 m in the alongshore direction and 1250 m in the cross-shore direction-for all three scenarios. This domain size was sufficient to detect non-uniform alongshore morphological development.
To account for the internal stochasticity of the model, we conducted 5 simulations for each scenario (i.e., repetitions) for both model representations. This number of repetitions is consistent with the number used in previous studies and meets the precision requirements recommended by Galiforni-Silva et al. [17] for the dune parameter predictions used in this study. We then compared the results of both model representations of the three scenarios after 50 years of simulation. This was considered a sufficient timespan for the full development of aeolian bedforms in vegetated sandy environments, such as nebkhas and parabolic dunes [45]. Additionally, we examined the spatio-temporal evolution of the beach-foredune system in the 50-year morphology of each of the scenarios.

Model Parameters DuBeVeg Default Representation
For the default representation of DuBeVeg (representing sand transport through ripple migration), we implemented the model parameters used by Galiforni-Silva et al. [47,48] to study idealized scenarios and a real beach-dune system on the Dutch coast ( Table 1). Most of these parameters are based on the calibration executed by Keijsers [46].

DuBeVeg Adaptation to Saltation
For the newly adapted representation of the model that considers saltation the main transport mode, we considered slabs of 1 × 1 × 0.02 m 3 ( Table 1). The slab height (h) was decreased by a factor of 5 to compensate for the 5-fold increase in the jumping length (L) in order to maintain the same potential aeolian transport per meter alongshore (Q; see Equation (3)) that was used by Galiforni-Silva et al. [47,48].
A vegetation effectiveness limit (ρ * ) of 0.3 was set for this study. This value was based on the same limit that is used in the DuBeVeg model for the threshold of the vegetation effectiveness required to affect the angle of repose of sand. This implementation was previously used by Nield and Baas [45] in the DECAL model. The parameters of the hydrodynamic and vegetation module are equal to the default representation.

Boundary and Initial Conditions
For all scenarios, identical boundaries and initial conditions were applied. A constant onshore wind was considered, as in previous implementations of DuBeVeg [46][47][48]. This means that a slab of sand, once it is picked up by the wind, can only jump in the cross-shore landward direction (toward the dunes). As Q is independent of the instantaneous wind speed (see Equation (3)), no time-series of wind data was required. Observed water levels at Scheveningen Harbor for a 50-year period from 1970 to 2020 were used as inputs for the hydrodynamic module (Equation (4)).
At the start of each simulation, vegetation was only present at the initial foredune. Every foredune cell above 6 m +MSL was fully vegetated, and for cells between 4 m and 6 m +MSL, a random value of ρ ranging between zero and one was assigned. The reference surface was defined as the initial beach surface. Excluding the dunes, analogous to Galiforni-Silva et al. [47,48]. In the narrow and wide scenario, the reference surface remained constant throughout the whole simulation. In the prograding scenario, the reference surface moved seaward at the same rate as the seaward extension of the beach. A groundwater depth factor of G = 0.8 was chosen for all scenarios, using the same value as Keijsers et al. [46].

Foredune Morphology after 50 Years
With a default jumping of length L = 1, the model predicted a seaward change in the foredune location after 50 years for all three scenarios (Figure 3a,c,e). For the narrow scenario, a new foredune formed directly in front of the original foredune, where both the initial foredune and the new foredune retained a shore-parallel orientation with a high degree of longshore uniformity (Figure 3a). For the wide scenario, after 50 years, the new foredune had a barchanoid sinuous shape, with smaller discrete barchan-shaped bedforms appearing on the beach, seaward of the new foredune (Figure 3c). The prograding scenario also showed a new foredune in front of the original foredune, with a barchanoid sinuous shape (Figure 3e). In this scenario, two smaller foredunes with a similar shape also formed closer to the shoreline. The initial foredune showed more irregularities than for the narrow and wide beach scenarios, with the presence of blowouts (i.e., topographic depressions) with adjacent depositional lobes after 50 years.
The adapted model-with an increased jumping length of L = 5 that accounts for saltation transport-predicted a different foredune morphology after 50 years for all three scenarios (Figure 3b,d,f). Overall, the adapted model predicted the seaward expansion of the original foredune, with a less distinct separation between the consecutive foredune ridges than the predictions with the default jumping length setting. Moreover, compared with L = 1, sand passing the foredune front reached greater distances inland. In the narrow scenario, similar to L = 1, a new foredune formed in front of the initial foredune, but the initial foredune lost its longshore uniformity with the presence of blowouts and depositional-lobe-type features (Figure 3b). In the wide and prograding beach scenarios, the new foredune ridge seaward of the initial foredune was predicted to be much wider (Figure 3d,f). This new foredune ridge showed an irregular topography with multiple small blowouts. Moreover, the initial foredune showed an irregular elevation after 50 years, with blowouts and depositional features but a lower height than the final original foredune of the narrow scenario.
The width and volume of the dunes predicted for each model run were calculated by considering the dune toe the most seaward point of the cross-shore profile, with an elevation above 3 m +MSL. The width of the dune zone was then obtained by determining the distance between this dune toe and the landward end of the cross-shore profile. This was done for all profiles within the model domain (400 in total). The dune volume was obtained by integrating the volume of all cells above the 3 m +MSL plane per profile within the model. For both properties, the results of the 5 model repetitions were aggregated, resulting in 2000 predictions of dune width and volume per scenario (Figure 4). In Appendix A (Figure A1), the results for each of the five repetitions are displayed individually, demonstrating that the predicted distribution of both properties was similar across all five model repetitions and that observed differences between the scenarios and jumping length settings were consistent with the aggregated data presented in Figure 4.

Foredune Morphology after 50 Years
With a default jumping of length L = 1, the model predicted a seaward change in the foredune location after 50 years for all three scenarios (Figure 3a,c,e). For the narrow scenario, a new foredune formed directly in front of the original foredune, where both the initial foredune and the new foredune retained a shore-parallel orientation with a high degree of longshore uniformity (Figure 3a). For the wide scenario, after 50 years, the new foredune had a barchanoid sinuous shape, with smaller discrete barchan-shaped bedforms appearing on the beach, seaward of the new foredune (Figure 3c). The prograding scenario also showed a new foredune in front of the original foredune, with a barchanoid sinuous shape (Figure 3e). In this scenario, two smaller foredunes with a similar shape also formed closer to the shoreline. The initial foredune showed more irregularities than for the narrow and wide beach scenarios, with the presence of blowouts (i.e., topographic depressions) with adjacent depositional lobes after 50 years. The adapted model-with an increased jumping length of L = 5 that accounts for saltation transport-predicted a different foredune morphology after 50 years for all three scenarios (Figure 3b,d,f). Overall, the adapted model predicted the seaward expansion of the original foredune, with a less distinct separation between the consecutive foredune ridges than the predictions with the default jumping length setting. Moreover, compared with L = 1, sand passing the foredune front reached greater distances inland. In the narrow scenario, similar to L = 1, a new foredune formed in front of the initial foredune, but the initial foredune lost its longshore uniformity with the presence of blowouts and depositional-lobe-type features (Figure 3b). In the wide and prograding beach scenarios, the new foredune ridge seaward of the initial foredune was predicted to be much wider (Figure 3d,f). This new foredune ridge showed an irregular topography with multiple small blowouts. Moreover, the initial foredune showed an irregular elevation after 50 years, with blowouts and depositional features but a lower height than the final original foredune of the narrow scenario.
The width and volume of the dunes predicted for each model run were calculated by considering the dune toe the most seaward point of the cross-shore profile, with an elevation above 3 m +MSL. The width of the dune zone was then obtained by determining the distance between this dune toe and the landward end of the cross-shore profile. This was done for all profiles within the model domain (400 in total). The dune volume was obtained by integrating the volume of all cells above the 3 m +MSL plane per profile within the model. For both properties, the results of the 5 model repetitions were aggregated, resulting in 2000 predictions of dune width and volume per scenario ( Figure 4). In Appendix A ( Figure A1), the results for each of the five repetitions are displayed individually, demonstrating that the predicted distribution of both properties was similar across all five model repetitions and that observed differences between the scenarios and jumping length settings were consistent with the aggregated data presented in Figure 4. In terms of dune volume after 50 years, the narrow beach scenario resulted in smaller predicted dune volumes compared with the wide and prograding beach scenarios for both jumping length settings (Figure 4a). With the default jumping length (L = 1), the wide beach scenario showed the greatest average dune volume, with a median value of 1563 m 3 /m, while the narrow beach scenario predicted the lowest median value of 1500 m 3 /m. With an increased jumping length (L = 5), the differences in dune volume between the scenarios became more pronounced. The wide beach scenario showed the greatest In terms of dune volume after 50 years, the narrow beach scenario resulted in smaller predicted dune volumes compared with the wide and prograding beach scenarios for both jumping length settings (Figure 4a). With the default jumping length (L = 1), the wide beach scenario showed the greatest average dune volume, with a median value of 1563 m 3 /m, while the narrow beach scenario predicted the lowest median value of 1500 m 3 /m. With an increased jumping length (L = 5), the differences in dune volume between the scenarios became more pronounced. The wide beach scenario showed the greatest average dune volume, with a median value of 2234 m 3 /m, while the narrow scenario predicted the lowest median value of 1611 m 3 /m. Comparing the dune volumes for the two different jumping length settings for the same scenario, the volumes obtained with L = 5 were significantly higher than those with L = 1 for all scenarios, including the narrow beach scenario (t-test, α = 0.05). Additionally, there was a greater spread in the predicted dune volumes for the narrow beach scenario, with L = 5. This greater spread can be attributed to the non-uniformity in the alongshore direction of the initial dune after 50 years (Figure 3b).
In terms of dune width after 50 years, the narrow beach scenario also resulted in significantly smaller values compared with the wide and prograding scenarios, regardless of the jumping length setting (Figure 4b). With the default unit jumping length (L = 1), the prograding scenario showed the greatest dune widths, with a median value of 304 m. This can be attributed to the new foredune ridges that formed farthest from the original foredune ( Figure 3e). With the increased jumping length (L= 5), the prograding beach scenario did not stand out and presented similar dune widths to the wide scenario. The wide and prograding beach scenarios predicted median widths of 230 m and 218 m, respectively, for L = 5, with the difference being significant (t-test, α = 0.05). Comparing the dune widths for the two different jumping length settings for the same scenario, the prograding beach scenario exhibited the largest difference. In contrast, the narrow beach scenario showed a smaller yet significant difference between the two jumping length settings, with median dune widths of 145 m for L = 1 and 150 m for L = 5 (t-test, α = 0.05). This difference was significant across all scenarios, even though an overlap in the boxplots was observed for the wide scenario. Notably, the wide beach scenario with L = 1 exhibited a much wider distribution, which can be attributed to the sinuosity of the most seaward foredune ridge after 50 years (Figure 3c).

Sediment Flux Evolution
To analyze and compare the beach-foredune evolution during the 50 years of model simulations, we quantified the sediment flux of the dune toe by counting the slabs that crossed a vertical plane located at the alongshore-averaged dune toe position. The dune toe position was updated after each model iteration, aligning it with the evolving beach-foredune morphology. For each iteration, the most seaward position with an elevation of 3 m +MSL was determined for every cross-shore profile within the domain and then averaged across the model domain (i.e., in the alongshore direction). The sand fluxes across the dune toe were then computed for each cross-shore profile and were also averaged in the alongshore direction.
For the default setting of the jumping length (L = 1), the sediment flux across the dune toe after the first two model iterations dropped to almost half its initial value, from 22 m 3 /m/yr to 10 m 3 /m/yr, for all three scenarios ( Figure 5-c). In the narrow beach scenario, the peaks in the sediment flux for the default setting of the jumping length (Figure 5a) occurred right after elevated water levels induced by storms (Figure 5d). These storms eroded the bedforms that emerged on the beach during periods with only tide-induced water level fluctuations and interrupted the duneward-directed aeolian sediment flux. For the wide beach scenario, the beach had a gentler slope, and the erosive/accretive hydrodynamic effects occurred at greater distances from the dune toe, reducing its impact on the sediment flux across the dune toe compared with the narrow beach scenario (Figure 5b). The prograding beach scenario shows the greater temporal variability of the sediment fluxes for L = 1 (Figure 5c), which is related to the increased variability of the dune toe location as new bedforms form closer to the shoreline (Figure 3e). With the adapted jumping length setting (L = 5), the sediment flux into the foredunes remained more stable over time and across the different scenarios, with predicted fluxes consistently closer to the potential aeolian transport defined by the model parameters (Q = 25 m 3 /m/yr; Figure 5). Consequently, sediment fluxes across the dune toe remained substantially larger for the adjusted model (L = 5) than for the default setting of the jumping length (L = 1) for the full duration of the 50-year simulation. The only exemption to this trend was the narrow scenario, where the default setting of the jumping length intermittently resulted in slightly higher sediment fluxes compared with the adapted setting ( Figure 5a). Furthermore, a gradual decline in the sediment flux was observed for the narrow beach scenario until it obtained a stable value after 20 years (Figure 5a). This sediment flux decrease followed a decrease in beach width because of the expansion of the foredune in the seaward direction ( Figure 3b). Conversely, a gradual increase in the sediment flux over time was observed for the prograding beach scenario, which corresponded with the increase in beach width, corroborating that the prograding scenario was transitioning from a narrow beach to a wider beach scenario (Figure 5c). The wide beach scenario maintained a relatively constant sediment flux over time with the adapted jumping length (Figure 5b).

Beach-Dune Dynamics
Analyzing the 50-year topographic evolution of the beach-dune system revealed two distinct modes of spatio-temporal foredune development in the model results. The first mode will be referred to as the landward propagation mode and showed aeolian bedforms developing on the beach surface and growing in volume as they propagated landward. This mode was observed in the default setting of the jumping length (L = 1) for the wide and prograding scenarios (Figure 6c,e). The second mode will be referred to as the seaward expansion mode and showed a steady seaward build-out of the initial foredune front in the seaward direction. This seaward expansion mode was observed in the narrow beach scenario with the default jumping length setting (Figure 6a) and in all three scenarios with the adapted jumping length (L = 5) (Figure 6b,d,f). While Figure 6 shows the evolution of the cross-shore profiles in the center of the model domain for only a single repetition of the With the adapted jumping length setting (L = 5), the sediment flux into the foredunes remained more stable over time and across the different scenarios, with predicted fluxes consistently closer to the potential aeolian transport defined by the model parameters (Q = 25 m 3 /m/yr; Figure 5). Consequently, sediment fluxes across the dune toe remained substantially larger for the adjusted model (L = 5) than for the default setting of the jumping length (L = 1) for the full duration of the 50-year simulation. The only exemption to this trend was the narrow scenario, where the default setting of the jumping length intermittently resulted in slightly higher sediment fluxes compared with the adapted setting ( Figure 5a). Furthermore, a gradual decline in the sediment flux was observed for the narrow beach scenario until it obtained a stable value after 20 years (Figure 5a). This sediment flux decrease followed a decrease in beach width because of the expansion of the foredune in the seaward direction ( Figure 3b). Conversely, a gradual increase in the sediment flux over time was observed for the prograding beach scenario, which corresponded with the increase in beach width, corroborating that the prograding scenario was transitioning from a narrow beach to a wider beach scenario (Figure 5c). The wide beach scenario maintained a relatively constant sediment flux over time with the adapted jumping length (Figure 5b).

Beach-Dune Dynamics
Analyzing the 50-year topographic evolution of the beach-dune system revealed two distinct modes of spatio-temporal foredune development in the model results. The first mode will be referred to as the landward propagation mode and showed aeolian bedforms developing on the beach surface and growing in volume as they propagated landward. This mode was observed in the default setting of the jumping length (L = 1) for the wide and prograding scenarios (Figure 6c,e). The second mode will be referred to as the seaward expansion mode and showed a steady seaward build-out of the initial foredune front in the seaward direction. This seaward expansion mode was observed in the narrow beach scenario with the default jumping length setting (Figure 6a) and in all three scenarios with the adapted jumping length (L = 5) (Figure 6b,d,f). While Figure 6 shows the evolution of the cross-shore profiles in the center of the model domain for only a single repetition of the model for each scenario, the discussed development modes of the beach-dune dynamics were observed in all cross-shore profiles and all repetitions. model for each scenario, the discussed development modes of the beach-dune dynamics were observed in all cross-shore profiles and all repetitions.
In the landward propagation mode, some of the aeolian bedforms that were propagating landward merged with the initial foredune, while others stabilized before reaching the initial foredune and formed a new foredune ridge in front of the initial foredune ( Figure  6c,e). At the same time, the bedforms also grew in the alongshore direction and connected with each other, forming a continuous ridge, with some alongshore variability coming from the barchan shape of the aeolian bedforms (Figure 3c,e). Zooming in on the beach, distinctively higher aeolian bedforms emerged on the beach for the default jumping length than for the increased jumping length (Figure 7). These emerging aeolian bedforms on the beach intercepted the sand transport to the foredune and simultaneously created non-erodible zones on their leeward sides. Combined, these effects reduced the sediment flux across the beach and into the foredune for the default jumping length compared with the increased jumping length ( Figure 5). For the default jumping length setting, the sediment flux dropped from a value of 22 m 3 /m/yr to a value of 10 m 3 /m/yr across most of the beach profile in the wide beach scenario (Figure 8a), except for parts of the beach profile closer to the shoreline. These higher fluxes close to the shoreline were accommodated by the erosive effect of the sea, which removed the bedforms shown in Figure 7a and allowed for a greater sediment flux on the part of the profile affected by seawater levels. Thus, the width of the beach profile with greater sediment fluxes varied with time, in accordance with the water level fluctuations. With the adjusted jumping length setting (L = 5), the lower bedforms that emerged on the beach (Figure 7b) appeared to be not tall enough to disrupt the sediment flux, resulting in a relatively constant and high sediment flux value of 24 m 3 /m/yr across the beach domain (Figure 8b). In the landward propagation mode, some of the aeolian bedforms that were propagating landward merged with the initial foredune, while others stabilized before reaching the initial foredune and formed a new foredune ridge in front of the initial foredune (Figure 6c,e). At the same time, the bedforms also grew in the alongshore direction and connected with each other, forming a continuous ridge, with some alongshore variability coming from the barchan shape of the aeolian bedforms (Figure 3c,e).
Zooming in on the beach, distinctively higher aeolian bedforms emerged on the beach for the default jumping length than for the increased jumping length (Figure 7). These emerging aeolian bedforms on the beach intercepted the sand transport to the foredune and simultaneously created non-erodible zones on their leeward sides. Combined, these effects reduced the sediment flux across the beach and into the foredune for the default jumping length compared with the increased jumping length ( Figure 5). For the default jumping length setting, the sediment flux dropped from a value of 22 m 3 /m/yr to a value of 10 m 3 /m/yr across most of the beach profile in the wide beach scenario (Figure 8a), except for parts of the beach profile closer to the shoreline. These higher fluxes close to the shoreline were accommodated by the erosive effect of the sea, which removed the bedforms shown in Figure 7a and allowed for a greater sediment flux on the part of the profile affected by seawater levels. Thus, the width of the beach profile with greater sediment fluxes varied with time, in accordance with the water level fluctuations. With the adjusted jumping length setting (L = 5), the lower bedforms that emerged on the beach (Figure 7b) appeared to be not tall enough to disrupt the sediment flux, resulting in a relatively constant and high sediment flux value of 24 m 3 /m/yr across the beach domain (Figure 8b).
which removed the bedforms shown in Figure 7a and allowed for a greater sediment flux on the part of the profile affected by seawater levels. Thus, the width of the beach profile with greater sediment fluxes varied with time, in accordance with the water level fluctuations. With the adjusted jumping length setting (L = 5), the lower bedforms that emerged on the beach (Figure 7b) appeared to be not tall enough to disrupt the sediment flux, resulting in a relatively constant and high sediment flux value of 24 m 3 /m/yr across the beach domain (Figure 8b).

Beach-Dune Dynamics and Sediment Transport Modes
The model results obtained with the two different jumping length settings revealed distinct spatio-temporal dynamics of the beach-dune system. Two different modes of spatio-temporal foredune evolution were observed: seaward expansion and landward propagation ( Figure 9). The landward propagation mode was only observed for simulations with the default jumping length (L = 1) and was characterized by the emergence of aeolian bedforms near the shoreline, which grew in volume as they propagated landward toward the foredune. These bedforms either formed a new foredune ridge in front of the initial foredune or merged with the existing foredune. Interestingly, the aeolian bedforms in the landward propagation mode could reach elevations higher than the designated limit of the dune toe (3 m +MSL) as they traveled landward (Figure 6c,e). This mechanism of landward propagation does not align with the traditional schematizations of dune progression proposed by Arens and Wiersma [56] and Hesp [57]. Typically, foredune evolution is characterized by vertical growth in existing foredunes or by the formation of new foredunes in the seaward direction (Figure 9a), without a disconnection between the foredune ridges as observed in the landward propagation mode (Figure 9b). Thus, existing schematizations of foredune growth are in agreement with the seaward expansion mode that was obtained with the increased jumping length for all beach width scenarios.

Beach-Dune Dynamics and Sediment Transport Modes
The model results obtained with the two different jumping length settings revealed distinct spatio-temporal dynamics of the beach-dune system. Two different modes of spatio-temporal foredune evolution were observed: seaward expansion and landward propagation ( Figure 9). The landward propagation mode was only observed for simulations with the default jumping length (L = 1) and was characterized by the emergence of aeolian bedforms near the shoreline, which grew in volume as they propagated landward toward the foredune. These bedforms either formed a new foredune ridge in front of the initial foredune or merged with the existing foredune. Interestingly, the aeolian bedforms in the landward propagation mode could reach elevations higher than the designated limit of the dune toe (3 m +MSL) as they traveled landward (Figure 6c,e). This mechanism of landward propagation does not align with the traditional schematizations of dune progression proposed by Arens and Wiersma [56] and Hesp [57]. Typically, foredune evolution is characterized by vertical growth in existing foredunes or by the formation of new foredunes in the seaward direction (Figure 9a), without a disconnection between the foredune ridges as observed in the landward propagation mode (Figure 9b). Thus, existing schematizations of foredune growth are in agreement with the seaward expansion mode that was obtained with the increased jumping length for all beach width scenarios. Furthermore, the size of the non-vegetated aeolian bedforms on th the landward propagation mode (Figure 7a) did not correspond with the b in real beach environments. Non-vegetated bedforms on the beach c from small ripples of a few centimeters [23] to larger sand patches wit centimeters [58][59][60][61]. For the default jumping length, which produ propagation mode in the wide and prograding beach scenario, bedfo surface were predicted to have heights in the order of a meter (Figure than the observed dimensions in the field. The seaward expansion mode, c increased jumping length, produced bedforms with a much smaller order of tens of centimeters), in agreement with observed bedform h The reduced height of these bedforms on the beach also reduced th sediment and inhibited their stabilization, allowing for more direct tra toward distant zones of the beach. Accordingly, the sediment flux over with the increased jumping length demonstrated a trend where increased with cross-shore distance until it reached a maximum flux trend was corroborated by the expected evolution of the sediment flux its maximum, as described by Davidson-Arnott et al. [62]. With the seaw the maximum obtained sediment flux was close to the potential aeo while the sediment flux observed with the landward propagation mode, jumping length of = 1, was closer to half the potential aeolian transp Field observations further support the alignment of natural beac with the seaward expansion mode (Figure 9a) that was obtained with the length but not necessarily with the default jumping length (only fo scenario). To illustrate this agreement, we compared our model  [56] and Hesp [57]. Numbers indicate the order of appearance of new aeolian bedforms. New aeolian bedforms are indicated with a lighter shade.
Furthermore, the size of the non-vegetated aeolian bedforms on the beach surface for the landward propagation mode (Figure 7a) did not correspond with the bedforms observed in real beach environments. Non-vegetated bedforms on the beach can range in height from small ripples of a few centimeters [23] to larger sand patches with heights of tens of centimeters [58][59][60][61]. For the default jumping length, which produced the landward propagation mode in the wide and prograding beach scenario, bedforms on the beach surface were predicted to have heights in the order of a meter (Figure 7a), much greater than the observed dimensions in the field. The seaward expansion mode, characteristic of the increased jumping length, produced bedforms with a much smaller amplitude (in an order of tens of centimeters), in agreement with observed bedform heights in the field. The reduced height of these bedforms on the beach also reduced their ability to catch sediment and inhibited their stabilization, allowing for more direct transport of sediment toward distant zones of the beach. Accordingly, the sediment flux over the beach obtained with the increased jumping length demonstrated a trend where the sediment flux increased with cross-shore distance until it reached a maximum flux (Figure 8b). This trend was corroborated by the expected evolution of the sediment flux on a beach toward its maximum, as described by Davidson-Arnott et al. [62]. With the seaward expansion mode, the maximum obtained sediment flux was close to the potential aeolian transport (Q), while the sediment flux observed with the landward propagation mode, as obtained with a jumping length of L = 1, was closer to half the potential aeolian transport (Figure 8a).
Field observations further support the alignment of natural beach-dune dynamics with the seaward expansion mode (Figure 9a) that was obtained with the increased jumping length but not necessarily with the default jumping length (only for a narrow beach scenario). To illustrate this agreement, we compared our model results with the morphological evolution of a beach-dune profile located 100 m south of Eierlandse Dam, which has experienced accretion and a prograding shoreline over the past 25 years. Eierlandse Dam is a jetty on the northwest beach of Texel Island in the Netherlands built in 1995 (Figure 10a). Since its construction, no nourishments have been administered around the location of the jetty. Annual surveys of the Dutch coast, with a resolution of 5 m on the beach and in the dune area [63,64], showed the widening of the cross-shore beach profile over a 25-year period from 1997 to 2022. In this period, the initial foredune widened in the seaward direction (Figure 10b), in line with the seaward expansion mode (Figure 9a). Similar behavior was also observed in field data from measurements of different beach-dune systems in, e.g., Australia and the USA [65][66][67]. The integration of an increased jumping length ( = 5) in the DuBeVeg results that aligned better with the beach-dune evolution observed in re systems. Specifically, the use of the increased jumping length of =5 cons the seaward expansion mode of foredune evolution, whereas = 1 tended t landward propagation mode of foredune evolution, except in the narrow b These findings highlight the importance of selecting an appropriate value f length in DuBeVeg, particularly for situations with wide or progradin previously noted by Werner [33], a greater value of L better mimics the sal sand transport on beaches. As the saltation mode is the prevailing transpor beaches [22,23], this explains the improved alignment of the model r increased jumping length (L = 5) with observed beach-dune dynamics.
By better representing the prevailing sediment transport mode, representation of the jumping length in DuBeVeg allowed us to capture th states that occurred during the beach-dune system's evolution, provid understanding of the interactions that occurred at these stages. In contrast of the default jumping length of L = 1 was mainly aimed at investiga morphology of beach-dune systems after a period of several decades. The the temporal dynamics of a beach-dune system is especially important beach-dune systems that experience significant human impacts. Such a m used to investigate the effects of beach nourishment on foredune develop formation of new dunes, which are not yet well understood [68]. T reproduction of the beach-dune system's transient behavior, obtained w The integration of an increased jumping length (L = 5) in the DuBeVeg model yielded results that aligned better with the beach-dune evolution observed in real beach-dune systems. Specifically, the use of the increased jumping length of L = 5 consistently led to the seaward expansion mode of foredune evolution, whereas L = 1 tended to result in the landward propagation mode of foredune evolution, except in the narrow beach scenario. These findings highlight the importance of selecting an appropriate value for the jumping length in DuBeVeg, particularly for situations with wide or prograding beaches. As previously noted by Werner [33], a greater value of L better mimics the saltation mode of sand transport on beaches. As the saltation mode is the prevailing transport mode on real beaches [22,23], this explains the improved alignment of the model results for the increased jumping length (L = 5) with observed beach-dune dynamics.
By better representing the prevailing sediment transport mode, the updated representation of the jumping length in DuBeVeg allowed us to capture the intermediate states that occurred during the beach-dune system's evolution, providing a deeper understanding of the interactions that occurred at these stages. In contrast, the prior use of the default jumping length of L = 1 was mainly aimed at investigating the final morphology of beach-dune systems after a period of several decades. The simulation of the temporal dynamics of a beach-dune system is especially important for managed beach-dune systems that experience significant human impacts. Such a model could be used to investigate the effects of beach nourishment on foredune development and the formation of new dunes, which are not yet well understood [68]. The enhanced reproduction of the beach-dune system's transient behavior, obtained with a greater jumping length, makes the cellular automata approach an attractive tool for studying the evolution of these sand-nourished beach-dune systems [69] and, eventually, support coastal management.

Foredune Morphology and Vegetation Effectiveness Limit
All simulated scenarios in our study resulted in a seaward expansion of the foredune(s) (Figure 3), which is attributed to a positive sediment input from the sea. The imposed onshore wind direction tended to erode the beach profile near the shoreline, and the hydrodynamic module strived to restore this shoreline erosion to the reference surface. Therefore, new sediment was introduced to the system by the hydrodynamic module. The expansion of the foredune in the seaward direction agreed with the conceptual schematizations in previous studies where the expansion of the foredunes was predicted for beaches with positive sediment input from the sea [19,56,57,70].
While both jumping length settings demonstrated an expansion of the foredunes, they yielded different morphologies after 50 years. These divergent outcomes can be attributed to the spatio-temporal differences in sediment fluxes that occurred within the two models with the different jumping lengths. The greater sediment fluxes for the increased jumping length of L = 5 on the beach (Figure 8) resulted in a greater flux of sediment into the foredune (Figure 5), causing a greater foredune volume in all scenarios (Figure 4). The disparity in the final foredune volume highlights the critical role of selecting the appropriate value of L, as the foredune volume serves as an important parameter for flood protection provided by dunes with a coastal defense function [71]. Furthermore, the increased sediment flux with the L = 5 setting led to a more discontinuous morphology landward of the dune toe. The increased sediment flux caused more sediment to cross the dune crest, burying the initial vegetation of the original foredune and leading to more erosive/accretive activity from the wind, causing the emergence of blowouts and depositional lobes.
The updated jumping length setting, in combination with the incorporation of interception caused by the vegetation (Equation (5)), introduced a variability to the effective jumping length (L) across the simulated domain, taking into account the vegetation density of the cells. In the context of beach-dune systems, it is important to consider the trapping capacity of vegetation, as even a vegetation cover of 30% has been found to significantly disrupt aeolian transport [72,73]. In the adapted model, the vegetation effectiveness limit (ρ * ) determines the density of vegetation that is required for it to behave similarly to any terminal cell with the default jumping length setting. A value of ρ * = 0.3 was selected based on the angle of the repose threshold for bare cells to transition into vegetated cells in DuBeVeg. When ρ * = 0, indicating that any cell with a vegetation density greater than zero is considered vegetated, the effective jumping length becomes L = 1 for all cells with any vegetation. The predicted topography after 50 years obtained with ρ * = 0, therefore, did not deviate significantly from the initial topography (Figure 11a), as the potential aeolian transport for a slab height of h = 0.02 m (for the increased jumping length) and an effective jumping length of L = 1 resulted in a potential aeolian transport rate of Q = 5 m 3 /m/yr (Equation (3)). Increasing the slab height to 0.1 m yielded a potential aeolian transport rate similar to the default jumping length setting and predicted a similar final morphology (Figure 3c). Conversely, when ρ * = 1 was applied, the model behaved similarly to the previous simulations for ρ * = 0.3 (Figure 11b,c), as most of the beach remained mostly unvegetated during the simulation, causing little difference between the two vegetation effectiveness limits. Consequently, this setting resulted in an effective jumping length of L = 5 for the majority of the transported slabs. Nevertheless, a distinct difference could be observed for the foredune area, where the vegetation density was significantly greater than at the beach. With a unit vegetation effectiveness limit of ρ * = 1, here, the sand encountered less resistance from the vegetated part of the foredune near the toe, leading to a narrower foredune after 50 years than with the original vegetation effectiveness limit (Figure 11b,c).

Conclusions
In order to represent saltation transport, the updated representation of the aeolian transport mode in the DuBeVeg model provides valuable insights into the emergence of bedforms on the beach and the temporal evolution of foredunes. Our study underscores the significance of selecting the appropriate transport mode when simulating the temporal evolution of beach-foredune systems. A fivefold increase in the jumping length in the DuBeVeg model better captured the saltation transport mode, rather than the ripple migration mode represented by the default unit jumping length. By increasing the jumping length, the model more accurately predicted the transient dynamics of beachdune systems, particularly those with wide or prograding beaches. The adapted model more accurately reproduced the emergence of bedforms on the beach, predicting much smaller amplitudes than the default representation of the model. Accordingly, the adapted model predicted increased sediment transport across various zones of the beachdune system, resulting in a more accurate representation of the typical seaward foredune development. With the adapted model better predicting the temporal evolution of beachdune systems, it will also be better suited for model studies into the effects of human activities on the evolution of these systems. For such studies, the transient evolution of the beach-dune system matters as much as its final state, given the interconnected temporal variability of human activities and the morphological development of the beach-dune system.

Conclusions
In order to represent saltation transport, the updated representation of the aeolian transport mode in the DuBeVeg model provides valuable insights into the emergence of bedforms on the beach and the temporal evolution of foredunes. Our study underscores the significance of selecting the appropriate transport mode when simulating the temporal evolution of beach-foredune systems. A fivefold increase in the jumping length in the DuBeVeg model better captured the saltation transport mode, rather than the ripple migration mode represented by the default unit jumping length. By increasing the jumping length, the model more accurately predicted the transient dynamics of beach-dune systems, particularly those with wide or prograding beaches. The adapted model more accurately reproduced the emergence of bedforms on the beach, predicting much smaller amplitudes than the default representation of the model. Accordingly, the adapted model predicted increased sediment transport across various zones of the beach-dune system, resulting in a more accurate representation of the typical seaward foredune development. With the adapted model better predicting the temporal evolution of beach-dune systems, it will also be better suited for model studies into the effects of human activities on the evolution of these systems. For such studies, the transient evolution of the beach-dune system matters as much as its final state, given the interconnected temporal variability of human activities and the morphological development of the beach-dune system. Figure A1. Boxplots of the alongshore variability of the dune volume and width per scenario predicted with both model settings for each of the 5 model repetitions: (a,c,e,g,i) show the dune volume for the narrow, wide, and prograding beach scenarios; (b,d,f,h,j) show dune width for the narrow, wide, and prograding beach scenarios.