1 Introduction

Coastal sand dunes worldwide are the product of maritime winds that are capable of generating large amounts of surface shear stress. Shear stresses drive saltation, which changes the surface topography of dunes, and under certain conditions cause mineral dust emissions [1,2,3,4]. In addition, vegetation is common in coastal dune settings, which complicates the distribution of surface winds, shear stress, and sand transport. Ultimately, these interactions between wind, sand transport, vegetation, and variable topography give rise to dunes of varying morphology [5]. The Santa Maria Dune complex [6, 7] on California’s central coast (Fig. 1) contains the Oceano Dunes State Vehicular Recreation Area (ODSVRA), California State Park. The dunes within the ODSVRA are known to emit dust during high winds that can result in the exceedance of federal and state air quality particulate matter 24 h mean standard for particles less than or equal to \(10\,\upmu {\textrm{m}}\) aerodynamic diameter (\(\hbox {PM}_{10}\)) downwind on the Nipomo Mesa [8,9,10]. Approximately 500 ha of this dune environment allows off-highway recreational vehicle (OHV) activity. OHV activity has been demonstrated to increase mineral dust emissions [9] from impacted surfaces compared to non-impacted surfaces in the ODSVRA. A Stipulated Order of Abatement initiated in 2018 and implemented by the San Luis Obispo County Air Pollution Control District for infractions of state and federal air quality standards requires that dust originating from the ODSVRA be reduced to meet a mass emission reduction target and the state and Federal \(\hbox {PM}_{10}\) standards downwind of the park.

Fig. 1
figure 1

A map of the ODSVRA, surrounding areas, and location within California

Full size image

As part of a strategy to reduce dust emissions from the dunes, restoration techniques to enhance development of a foredune system are being evaluated [11]. Currently across much of the north–south extent of the beach area of the ODSVRA foredunes are absent, but they are established in the south of the park where OHV activity is not permitted. Historical aerial photography of the site indicates that hummocky foredunes did exist throughout the ODSVRA back into the 1930s, but have since been destroyed by human activities [12]. The expectation is that a well-developed foredune system will disrupt incoming boundary layer flow from the ocean, reduce surface shear stress, and decrease sand flux across its width and for some distance on the lee side, thereby reducing saltation-driven dust emissions and contributing to improved air quality downwind. Quantifying the benefit of the presence of a foredune on modulating sand transport and dust emissions is a challenge. Here we use a computational fluid dynamics (CFD) approach to characterize the value of foredunes typical of this part of California to assess its aerodynamic value to mitigate dust emissions.

CFD has gained prominence as a tool that can be applied to understand air flow in the natural [13,14,15] and built [16, 17] environments. Compared to point measurements made in the field, CFD simulations result in very dense datasets, offering a detailed view of the fluid flow of interest. However, simulations require boundary data and are subject to difficult computational and numerical issues. In many cases, available computational power may limit the size and scope of simulation. For flows of interest in micrometeorology, it is possible to combine sparse measurements to drive CFD simulations over landforms [15, 18]. This approach is advantageous, as it combines the strengths of both methods. Moreover, modification of simulation settings can be used to model conditions that may not be readily observable in the field.

Dune morphology studies have encompassed both experimental [19,20,21] and CFD approaches [14, 22]. Simulations of flow have been performed over idealized nebkha [23], reversing barchans [18], star dunes [20], and Martian sand ripples [19]. However, few CFD studies have studied the combined effects of two or more dune types in succession in the downwind direction on flow parameters. Multiple studies have validated the use of CFD in the simulation of flow over landforms [15, 18, 24]. An understanding of airflow over dunes offers insight into dune morphology and sediment transport. Flow separation behind dunes produces regions of shelter, and has been shown to vary with both dune type and wind speed [13, 25, 26]. Understanding of this effect may lead to a better understanding of dune height and periodicity for dunes produced under consistent wind fields. The coupling of saltation and flow field data can be used to predict the evolution of land forms, which is particularly useful for foredune restoration applications [27, 28].

In this paper, we focus on landscape-scale simulations of flow over a coastal dune system consisting of a flat sloping beach, an area of nebkha foredune followed by a flat upward sloping zone that grades into a series of large transverse dune forms. High resolution topographic data provide the digital surface model (DSM) of the dune system. Measurements of wind and turbulence characteristics in the modeling domain are used to provide a boundary condition and verify model-derived estimates of the flow parameters (e.g., wind speed, turbulence intensity, and shear stress). This work is motivated by the need to characterize the role of the dune morphology to modulate the sand transport and associated dust emissions and guide management decisions to improve air quality downwind of the ODSVRA. By digitally modifying the topography and performing paired simulations, we elucidate the underlying aerodynamic mechanisms that affect dust emissions and evaluate the contribution of the foredune, the vegetation on the foredune, and the transverse dunes to modulate the shear stress compared with a flat sloping sand surface representing the absence of the foredune in the landscape.

2 Methods

2.1 Study site

Foredunes in this part of California, as typified by the Oso Flaco study area of the ODSVRA, are typically an assemblage of nebkhas, i.e., hummocky sand deposits topped with vegetation with elongated shadow dunes aligned with the resultant sand transport direction. Foredune formation, in general, depends on the presence of vegetation, sand supply, and other factors such as wind regime, surfzone beach type, plant species present, vegetation density and distribution, and the frequency and impacts of storms [21, 29,30,31]. Hesp et al. [32] suggest that nebkha foredune development on a coast is influenced strongly by the precipitation regime, with vegetation (density and morphology), sediment availability, and wind activity playing subordinate roles. The vegetation component of the nebkhas in the study area, identified as the north Oso Flaco dunes, are typically the following plants: Abronia laitfolia (yellow sand verbena) Abronia maritima (red sand verbena) Achillea millefolium (common yarrow), Ambrosia chamissonis (beach bur), Atriplex leucophylla (beach saltbush), Camissoniopsis cheiranthifolia (beach evening primrose) Erigonum parvifolium (coastal buckwheat), Eriophyllum staechadifolium (seaside golden yarrow), Fragaria chiloensis (beach strawberry), and Malacothrix incana (dunedelion).

2.2 Topography creation

Detailed aerial imagery of the Oso Flaco beach, nebkha foredune, and back dune regions was acquired using a fixed-wing, fully autonomous WingtraOne uncrewed aerial system (UAS) equipped with a Sony RX1RII 42 MP full-frame sensor and flown at approximately 110 m above ground level to collect high resolution (\(\le 2\,\hbox {cm}\)) orthographic RGB imagery. An on-board 5-band multispectral sensor provided data to identify vegetation. The details of the data collection and processing are provided in [11]. Photogrammetry was performed on the images to obtain a DSM of resolution better than \(10\,\hbox {cm}\).

To be used in a CFD context, additional processing was performed on the topographic dataset for the modeled domain. A geotiff of the region was exported to a point cloud, which was made into a 3D surface using 2D Delaunay triangulation. This results in a topographic surface that can be used to create a computational finite volume mesh (see Sect. 2.4). This surface is shown in Fig. 2.

Fig. 2
figure 2

A 3D representation of the unmodified topography of the mature Oso Flaco foredune. The shoreline is located on the left side of the surface

Full size image

Two additional topographies were created that test hypothetical scenarios. In the first, the foredune is effectively flattened to a gently sloping surface, which represents much of the current topography present within the ODSVRA along the extent of the north–south beach area. This is accomplished in the DSM by slicing a rectangular region covering the foredune out of the 3D pointcloud used to create the surface shown in Fig. 2. A data processing code was written in Julia to re-assign elevation values to the rectangular region, based on inverse distance weighting. The same 2D Delaunay triangulation process was then used to create a 3D surface. The result is a smooth topography, where the rectangular region containing the nebkha foredune has been replaced with a gently sloping surface. This sloping surface, shown in the upper panel of Fig. 3, is the result of spatially interpolating the edges of the clipped areas together. In the second scenario, a horizontal surface is desired behind the foredune. This is accomplished by clipping the pointcloud used to create Fig. 2, and interpolating a transition zone behind the foredune, which smoothly joins the foredune topography to a flat plane. The height of this plane is 4.2 m above sea level, which is based on the average height of the beach. The resulting 3D surface topography can be viewed in the lower panel of Fig. 3.

Fig. 3
figure 3

Two additional surface topographies created to test hypothetical scenarios. The upper panel shows the land surface with a flattened foredune; the lower panel shows the foredune topography joined to a flat plane on its leeward side

Full size image

2.3 Atmospheric measurements

A velocity measurement campaign was carried out with six sonic anemometers (R.M. Young, model 81000 three axis (U,V,W)) mounted on two towers to collect data for the purposes of boundary condition creation and simulation validation. Sonic anemometers were mounted to two \(3\,\hbox {m}\) mobile towers at heights of \(0.25\,\hbox {m}\), \(1.56\,\hbox {m}\) and \(3.26\,\hbox {m}\). The two towers were transported to regions of interest, with one tower installed in a static upwind location on the beach \({10}\,\hbox {m}\) in front of the foredune, and the other moved along a west to east transect of the foredune and left in place to collect data for approximately \({45}\,\hbox {min}\) per location. A photo of the upwind tower used to measure data for creating the simulation boundary condition is shown in Fig. 4. An aerial view of the tower locations within the Oso Flaco foredune area can be seen in Fig. 5.

Fig. 4
figure 4

Tower configuration for data collected with the sonic anemometers

Full size image
Fig. 5
figure 5

Monitoring locations for the sonic anemometers. Sites P1–P4 were visited by one mobile tower, while another was held fixed at the upwind location

Full size image

Sufficient wind velocity data were collected to build an inlet boundary condition and validation points for numerical simulation. An average velocity profile was constructed over the range of instrument heights by averaging over the time periods of 1:00 PM to 2:00PM on 05-20-2021, which coincided with observed persistent saltation on the beach and along the measurement transect in the nebkha foredune. A log-law velocity profile was fit to values from the three sonic anemometers, for the purpose of providing velocity boundary values to a domain extending to a height of \({100}\,\hbox {m}\). A roughness length of \(z_0=.002\) was used based on ocean roughness [33], resulting in a fitted friction velocity value of \(U_\star =.45\). Turbulence intensity I is calculated as,

$$\begin{aligned} I = \frac{u^\prime }{U}, \end{aligned}$$
(1)

where \(u^\prime\) is the turbulent fluctuation of velocity, and U is the average velocity. A profile of I was created by linearly interpolating calculated values between the three anemometers, and extending the value at the top anemometer to the top of the computational domain. A profile of turbulent kinetic energy k is calculated from the average velocity U and turbulence intensity I as,

$$\begin{aligned} k= \frac{3}{2}\left( UI \right) ^{2}. \end{aligned}$$
(2)

The specific dissipation rate \(\omega\) can then be calculated,

$$\begin{aligned} \omega = C_{\mu }^{.75} \frac{k^{\frac{1}{2}}}{l}, \end{aligned}$$
(3)

where the constant \(C_{\mu }=.09\), and the turbulence length scale is taken to be \(l=5\). This value of turbulence length scale is an hypothetical choice based on the average elevation drop from sea level to the high tide, frequent presence of large ocean waves, and common values for a neutrally-stratified atmosphere [34]. The true value varies widely in both space and time, though this has little impact on the goals of this paper. The kinematic energy eddy viscosity can then be calculated,

$$\begin{aligned} \nu _{t}=.31 \frac{k}{\omega }. \end{aligned}$$
(4)

Profiles of U and I are given in Fig. 6. Boundary conditions for the other boundary patches resemble that of a typical wind tunnel simulation. Zero-gradient pressure was used at the inlet, ground, side, and top walls, while a fixed value of zero was applied at the exit. Slip velocity was applied to the top and side walls, no slip velocity was applied to the ground, and the outlet was given a zero gradient condition.

Fig. 6
figure 6

Profiles of U (left) and I (right) are constructed for the inlet boundary condition. Data points measured near the ground are shown as blue dots

Full size image

2.4 Computational methods

Simulations were performed by numerically solving the steady-state incompressible Navier–Stokes equations,

$$\begin{aligned} \rho {\textbf{u}}\cdot \nabla {\textbf{u}} = -\nabla p + \rho {\textbf{g}} + \mu \Delta {\textbf{u}}, \end{aligned}$$
(5)
$$\begin{aligned} \nabla \cdot {\textbf{u}}=0, \end{aligned}$$
(6)

where, \({\textbf{u}}\) is velocity, p is pressure, \({\textbf{g}}\) is acceleration due to gravity, and \(\mu\) fluid is viscosity. This isothermal model corresponds to the atmospheric case of a neutrally stratified atmosphere. This was chosen as it is the case with the most general applicability. Solutions were obtained with the computational finite volume toolbox openFoam [35] using the solver simpleFoam. Turbulence modeling was performed with a Menter’s shear stress transport (SST) turbulence model [36], chosen for its low sensitivity to normalized wall distance \(y+\), as the large and complex domain of interest presents a difficult meshing problem.

The computational domain is a \(900\,\hbox {m} \times 200\,\hbox {m} \times 100\,\hbox {m}\) prism, rotated about the vertical axis to align with the prevailing wind direction parallel to the \({900}\,\hbox {m}\) fetch. An angled plane was used to cut the inlet boundary, ensuring that the inlet boundary condition is applied across the beach at the most uniform elevation possible. The beach acts as an inlet fetch for the imperfect alignment of the boundary condition and beach elevation at any given point. Results generated from the simulations only include data from regions downwind of the start of the foredune. An illustration of the domain imposed over the topography can be seen in Fig. 7. A domain height of \({100}\,\hbox {m}\) was chosen to minimize flow acceleration that may occur due to changes in cross-sectional area caused by the topography. The chosen height exceeds five element heights of the largest topographic features in the domain.

Fig. 7
figure 7

An illustration of the computational domain imposed over the topography. The angled inlet is shown in blue to apply the inlet boundary conditions on the beach across a uniform distance from the start of the foredune

Full size image

Meshing was performed with the CFmesh utility. A maximum cell size of \({10}\,\hbox {m}\) was reached in the upper atmosphere, which is not a region of interest, to conserve computing resources. Cells were progressively refined with decreasing elevation to \({0.3}\,\hbox {m}\) within \({1}\,\hbox {m}\) of the ground. The addition of 7 layers against the ground surface brought the smallest wall-adjacent cells to heights of approximately \({0.009}\,\hbox {m}\). Total cell counts for the simulations performed vary by topography, but are on the order of 70 million for each scenario. For the unmodified geometry, a histogram of the \(y+\) values on the ground surface is shown in Fig. 8. The average value of \(y+=130.2\) is within the appropriate range for the \(k-\omega\) turbulence model when used with high Reynolds number wall functions in the log-law region [37, 38].

Fig. 8
figure 8

The histogram of \(y+\) values on the ground boundary for simulated flow over the unmodified Oso Flaco topography. The mean \(y+\) value is 130.2

Full size image

3 Results

3.1 Comparison with field data

The dataset measured in the field comes from two towers: one stationary tower positioned at the shoreline boundary of the computational domain, and a second that was moved throughout the day. Consequently, a direct validation cannot be carried out at all of the second tower’s locations for the time period that was used to construct the boundary condition. Moreover, the atmospheric conditions at the site evolved throughout the day, while the boundary condition only corresponds to an hour of data. Significant thermal effects likely existed at the site due to solar heating of the sand and cold onshore winds, though simulation was perfomed with isothermal flow equations for the purpose of generality. We compare the ratio of upwind velocity (i.e., measured at the beach) to the velocity measured at the mobile tower foredune locations from Fig. 5. A summary of these data are shown in Table 1. Over the course of data collection, significant changes were observed in both wind speed and direction. In particular, the difference in wind direction between the upwind and downwind tower increases with distance into the foredune.

Table 1 A summary of mean wind speed (WS) and direction (WD) for the validation times at p1–p4, and the upwind tower shown in Fig. 5, for the anemometer placed at \({3.26}\,\hbox {m}\). All data were collected on May 20, 2021
Full size table

A plot of wind speed ratio for the sonic anemometers at \({0.25}\,\hbox {m}\) and \({3.26}\,\hbox {m}\) is shown in Fig. 9. An electrical malfunction prevented data from being collected by the sonic anemometer placed at \({1.56}\,\hbox {m}\) on the downwind tower. The result in Fig. 9 should be interpreted in the context that the measured data were collected in a natural wind field, which varied in intensity, direction, and thermal gradient with field measurements taken in the lee of dune structures that often exceeded the height of the sonic anemometers. By contrast, the simulated results are modeled with a steady-state boundary condition and incompressible flow equations. Slight variations in GPS accuracy affecting the placement of the towers and measurement of the surface can result in the difference between measurements taken in the wake of surface features versus open area. Measured data were collected over the course of a day, which resulted in differing average velocities and Reynolds numbers being used to compute the wind speed ratio at each location. Thermal effects developed over the course of data collection, but these were not modeled. This is particularly evident for P3 and P4, where the land slope, distance from the ocean, and time of day increased the heating of the ground. At these points, the discrepancy is greater for the measurement height of \({0.25}\,\hbox {m}\), as a strong thermal gradient exists close to the sand surface.

Fig. 9
figure 9

Wind speed ratio for the four validation points shown in Fig. 5

Full size image

3.2 Flow and shear stress visualization

Wall shear stress was calculated from the fluid flow over the unmodified Oso Flaco surface topography. A visualization of the shear stress magnitude is shown in Fig. 10. As seen in the upper panel of Fig. 10, the foredune is characterized by spots of high shear stress. This is a result of the vegetated nebkahs being the first major obstacle encountered by the inlet wind profile. The low elevation area behind the nebkha foredune is a region of low shear, while the increasing elevation of the non-vegetated transverse dune system results in an elevated shear magnitude (see the lower panel of Fig. 10).

Fig. 10
figure 10

A visualization of the shear stress on the topography. A view of the foredune looking downwind is shown in the upper panel, while a view of the non-vegetated dune system looking upwind is shown in the lower panel

Full size image

Near-wall streamlines were computed for passive particles released into the flow at a height of \({10}\,\hbox {cm}\) above the ground surface. As the simulations are steady-state, these can also be interpreted as average pathlines for air molecules. Visualizations of the streamlines for the foredune and non-vegetated transverse dune system behind the foredune can be seen in Fig. 11. This visualization provides insight into the aerodynamics of the two dune systems. The individual vegetation-topped hummocks within the foredune are characterized by very small regions of flow separation, as their shape is governed by the equilibrium between plant sheltering and saltation, resulting in elongated, teardrop-like forms in the downwind direction. In contrast, the unvegetated transverse dune system shows significant flow separation behind the slip faces.

Fig. 11
figure 11

Upper panel: Near-wall streamlines are shown at the foredune, viewed looking upwind. Lower panel: Near-wall streamlines are shown for the non-vegetated dune behind the foredune, looking upwind

Full size image

The aerodynamics of the two dune types motivates an investigation into the role of plants in the shear stress distribution across the foredune. It can be seen from Fig. 10 that the individual hummocks receive large amounts of shear, and the low areas between them less so, despite the lack of flow separation. A raster vegetation mask of the foredune was created from spectral data acquired with the UAS system. A map of this raster is shown superimposed over the magnitude of shear stress in Fig. 12. Areas of high shear are stabilized by plant cover, while leeward channels of bare sand surface characterized by low shear where the saltation flux will be lower or zero if the shear stress is below the threshold for saltation. The effect of plants on this system are quantified in the following sections and the implication for their effect on dust emissions are evaluated.

Fig. 12
figure 12

A raster vegetation mask of the foredune superimposed on the foredune shear stress magnitude mapping. Black denotes locations of vegetation at \({10}\,\hbox {cm}\) scale as viewed from above. The upwind boundary of the simulation is shown at the left of the image, with flow from the left to right

Full size image

3.3 The sheltering effect of the foredune

The presence of the foredune inevitably has downwind effects on the wind profile and resultant shear stress. However, these are difficult to quantify in field studies, as the foredune cannot be removed for a paired study and observations would inevitably be made under variable atmospheric conditions. In simulation, these obstacles are overcome without high costs or environmental damage. This motivates the creation of the geometry shown in the top panel of Fig. 3, where the foredune has effectively been removed, and replaced with a smooth sloping surface. Shear from simulation of flow over this surface can be compared with shear from the simulation shown in Fig. 10, to determine the effect of the presence of the foredune.

A comparison between the two topographies is achieved by performing a surface integral of shear \(\tau\) [\(\hbox {N}\,\hbox {m}^{-2}\)] over each surface. This integral is computed over the downwind (transverse dunes) section of the unmodified foredune and flattened foredune simulations, where the two topographies are the same. This corresponds to approximately the right hand sides of the topographies shown in Fig. 2 and the upper panel of Fig. 3. The result is \({27{,}009.0}\,\hbox {N}\) for simulation over the unmodified topography, and \({27{,}722.2}\,\hbox {N}\) for simulation over the flattened foredune topography. This difference in shear can be attributed to changes in the wind field as it passes over the foredune, with less energy extracted near the ground for the flattened surface. For the boundary condition tested in these simulations, removal of the foredune would result in a 2.6 % increase in shear on the non-vegetated transverse dune system.

For the entire domain of the two simulations, the total shear force on the flattened geometry of \({47{,}949.7}\,\hbox {N}\) exceeds the total shear of the unmodified geometry of \({46{,}769.1}\,\hbox {N}\). This is notable, as the unmodified geometry of the foredune has a higher surface area than the flattened surface due to its topographic complexity. The total surface area of the unmodified geometry is \({167{,}590.5}\,\hbox {m}^{2}\), compared to \({164{,}083.5}\,\hbox {m}^{2}\) for the geometry with the foredune removed (a difference of \({3507.0}\,\hbox {m}^{2}\)). Moreover, the foredune topography is characterized by spikes of high shear on the vegetated nebkha, as seen in the upper panel of Fig. 10. For the boundary condition tested in these simulations, removal of the foredune would result in a 2.5 % increase in shear on the entire dune system. In addition to significant shelter effects behind the foredune, the total effect of the foredune topography is a reduction in shear as compared to a flattened surface, regardless of the effects of vegetation on the wind field or saltation. Taking into account the area covered by vegetation induces a further decrease in shear. The area-integrated shear of the unmodified geometry excluding vegetated areas is \({42{,}425.3}\,\hbox {N}\). The total shear on the geometry with the foredune removed is a 13.0 % increase.

For the sake of generalization, it is useful to present the surface shear normalized by surface area for the two types of dune systems present in the simulation. This is accomplished by integrating shear over the surface of the nebkha foredune and transverse dune from the simulation with the unmodified geometries, and dividing by the respective areas. For the foredune, a normalized shear of \({0.231}\,\hbox {N m}^{-2}\) was observed (inclusive of vegetated area), while the equivalent region from the simulation with the flattened geometry yielded a value of \({0.302}\,\hbox {N m}^{-2}\). The transverse dune area has a value of \({0.275}\,\hbox {N m}^{-2}\) in the pure simulation. This value increases slightly to \({0.285}\,\hbox {N m}^{-2}\) in the simulation with the foredune removed, as a result of less shelter and upward deflection of the wind field in the absence of the foredune.

The shear stress analysis can be extended to evaluate the effect on dust emissions by incorporating surface emission relations developed from emissivity data collected at the ODSVRA with the PI SWERL instrument [39]. This instrument uses a rotating blade in an enclosed chamber to generate multiple shear velocities on the sand surface. For robustness, we include analyses with average emissions relations for the riding and non-riding areas of the ODSVRA [9]. The emissions relation for the riding area is made from averaged data collected inside the red polygon shown in Fig. 1, while the non-riding emissions relation is made from averaged data collected outside the red polygon and inside the blue polygon. The average emissions relation E (\(\hbox {mg}\,\hbox {m}^{-2}\,\hbox {s}^{-1}\)) for the riding area as a function of shear velocity \(u^\star\) is,

$$\begin{aligned} E = 23.65 (u^\star )^{5.59}, \end{aligned}$$
(7)

while the average emissions relation for the non-riding area is,

$$\begin{aligned} E = 21.51 (u^\star )^{6.85}. \end{aligned}$$
(8)

These relations are created by fitting a power law to three \(u^\star\) levels applied by the PI SWERL instrument on the sand surface. After converting the simulated shear stress values to shear velocity, \(u^\star = \sqrt{\tau \rho ^{-1}}\), this expression can be integrated over the land surface area, resulting in total emissions units of \(\hbox {mg}\,\hbox {s}^{-1}\) and area-normalized units of \(\hbox {mg}\,\hbox {s}^{-1}\,\hbox {m}^{-2}\). Total emissions and the corresponding total shear force are shown in Table 2. The third scenario of an unmodified topography shows significantly lower emissions than the scenarios that neglect vegetation sheltering or have a flattened foredune.

Table 2 Shear and emissions totals for the surfaces with the foredune removed (Scenario 1), the unmodified surface neglecting the effects of plant cover (Scenario 2), and the unmodified surface with zero emissions from vegetated surfaces (Scenario 3)
Full size table

3.4 Foredune shadow effects

Being able to quantify the sheltering effects of the foredune is a useful tool for quantifying the shear stress reduction in the lee of foredune restoration projects. This is the motivation for the geometry created with the flat plane behind the foredune, as seen in the middle panel of Fig. 3. As the actual shear on the land surface may vary widely depending on the topography behind the foredune, the flat surface offers the best scenario for a generalized result that can be applied to account for shear stress reduction in the lee of restoration projects like those implemented at the Oceano Dunes. Shear stress was calculated from a simulation for this surface using the boundary conditions outlined in Sect. 2.4. A plot of the mean shear stress magnitude calculated in 1 m bins transverse to the freestream wind direction is shown in Fig. 13. The recovery of this curve illustrates the sheltering behind the foredune. Figure 14 shows this relation normalized against the maximum mean shear stress on the plane behind the foredune. This curve defines the foredune sheltering effect with shear reaching approximately 90 % of its potential value within \({50}\,\hbox {m}\) of fetch in the lee of the foredune.

Fig. 13
figure 13

Mean shear magnitude as a function of distance in the freestream direction. The leeward edge of the foredune is located at \(x=0\), with the shoreline in located at the left side of the plot

Full size image
Fig. 14
figure 14

Mean shear \(\tau\) normalized by maximum shear \(\tau ^\star\) as a function of distance

Full size image

3.5 Saltation threshold

Given the availability of data at the site and shear stress results from the simulation, regions above the saltation threshold can be mapped for the boundary condition used. A mean grain diameter of \({380}\,\upmu \hbox {m}\) was measured from sand samples taken at the study site. Application of the Bagnold threshold equation [40] results in a threshold shear velocity of \({0.28}\,\hbox {m s}^{-1}\), equivalent to a shear stress magnitude of \({0.097}\,\hbox {N m}^{-2}\). For the unmodified and removed nebkha foredune topographies, shear stress magnitude above this value is shown in Fig. 15. Surface area was integrated over the unmodified topography in the same transverse \({1}\,\hbox {m}\) bins used in Fig. 13. Figure 16 shows the surface area above threshold as a function of freestream direction. The blue line denotes when vegetated area is included as an erosive surface (EV), while the red line denotes vegetated area not included as a erosive surface (NEV). The effect of vegetation is primarily seen for negative distance values, with the value of zero being the arbitrarily chosen location of the downwind extent of the foredune, used in the alternate land surfaces seen in Fig. 3. The area above threshold is shown to rapidly decrease as the foredune is entered from the upwind side. Discrepancy between the two curves seen after a distance of \({450}\,\hbox {m}\) is a result of the transverse dune unevenly ending and giving way to shrub and tree cover in the upper right corner of the domain (see Fig. 12).

Fig. 15
figure 15

Locations above the Bagnold saltation threshold, colored by shear stress magnitude. Removed (white) areas represent surface regions below the threshold for sand transport. Above: the unmodified topography, below: the topography with the foredune removed

Full size image
Fig. 16
figure 16

Surface area above threshold, summed in \({1}\,\hbox {m}\) bins, transverse to the stream-wise wind direction. Vegetation effects are seen in the foredune (distance < \({0}\,\hbox {m}\)), and at the end on the end of the transverse dune (distance \(>{400}\,\hbox {m}\)). EV denotes vegetated area above threshold counted as an erodible surface, NEV denotes vegetated area counted as a non-erodible surface

Full size image

4 Discussion

Qualitatively, the vegetated foredune and non-vegetated transverse dune show distinctly different flow patterns. Extreme shear values in the foredune are highly localized, while shear across the transverse dune primarily changes as a function of streamwise distance. However, both dune types feature extended regions of shelter in the average flow field. While this paper makes efforts to separate the role of vegetation from topography, it is important to note that vegetation was not explicitly modeled in the simulations that were performed. This is justified by the very low profile of the vegetation found in the foredune (on the order of \({10}\,\hbox {cm}\)), and the generally coarse (by traditional CFD standards) mesh resolution used to achieve simulations over this large spatial extent. While the effects of the vegetation on the wind field are assumed to be small in this setting, it is clear that they play a primary role in shaping the topography of the foredune and affect the zones of saltation in their lee. For this reason, the modeling of flow-roughness interactions over the foredune merits further research, where an appropriate momentum extraction can be made from the Navier–Stokes equations (5) to model the effects of vegetation canopies on the wind field.

The visualization of flow separation occurring on the transverse dunes, shown in Fig. 11, shows significant flow features in the direction transverse to the freestream. This underscores the importance of 3D simulation in understanding the flow dynamics, sand transport pathways, and resulting morphodynamics of dune fields. Mechanisms for the complex flow observed are likely a combination of variations in both the height of individual transverse dunes and their wavelength. In many places, streamlines seeded at \({10}\,\hbox {cm}\) above the ground are seen leaving the crest of a transverse dune and quickly reconnecting to the surface below. In other places transverse motions combine to form areas of low velocity that persist over the crest of the following dune. This complicates the concept of flow separation, which is found to be highly variable over this topography. Baddock et al. [26] reported field-measured separation, finding lengths of 2.3h to 5h (where h is crest height with respect to the interdune zone) for crest-brink type barchan dunes. Performing measurements on the average flow lines of Fig. 11, we find similar values for the transverse dune in the range of 3.2h to 4.7h for this boundary condition. However, this is a crude approximation, given the complexity of flow observed transverse to the free-stream wind direction. For example, the helical flow observed in Fig. 17 emanates from a point and has a continuously evolving length, making it difficult to define a measurement for the individual form. The simulations performed here suggest that the true separation length in nature likely diverges from idealized cases most of the time. Flow separation behind the transverse dunes correlates closely with the saltation map of Fig. 15. Comparison of these two visualizations indicate that separation length is not solely responsible for shelter behind the transverse dune: highly irregular patterns of below-threshold surface are observed behind the transverse dunes due to low velocity attached flow. These features can be seen in an overlay of the near-surface streamlines and the actively saltating surface, shown in Fig. 17. Streamlines show that shelter zones extend beyond points of flow reattachment, and that flow separation is not always sufficient to stop saltation.

Fig. 17
figure 17

Streamlines displayed over a map of actively saltating surface. A helical flow (red outline) is observed to produce little shelter. Regions of extended shelter form under reattached flow streamlines (yellow outline) in the lower part of the image

Full size image

While this study has quantified characteristics of the foredune and transverse dune separately, the two regions are strongly coupled in that the atmospheric boundary layer develops over the foredune before it encounters the transverse dune. The 2.6 % increase in shear observed over the transverse dune (see Sect. 3.3) between the flattened and unmodified foredunes is more than the 2.5 % increase observed over the entire domain. However, these numbers are potentially misleading in that the distribution of shear in the foredune is characterized by extreme spikes of shear on the nebkhas, and relatively low shear between them. This is evidenced by visualization of the shear in Fig. 10, and the prevalence of zones below salatation threshold seen in Fig. 15. Additionally, the foredune topography contains more surface area than the flattened surface, and the nebkhas responsible for this increase in surface area, compared with the plane area, are covered by vegetation at what would be the position most exposed to the highest shear. This is illustrated by Fig. 16, which shows that total area above saltation threshold drops consistently as the foredune removes momentum from the boundary layer. However, when the surface area covered by vegetation is considered, this curve has a sudden drop and then a gradual increase with downwind distance through the foredune. Thus in the context of wind shelter, shear stress reduction from the presence of the foredune is carried to the transverse dune. When considering the potential to drive dust emissions, small regions of extreme shear may result in higher emissions than large regions of moderate shear, due to the power law nature of the emissions relations (Eqs. 7 and 8). It is clear that the area-normalized shear on the unmodified surface is significantly lower than if the foredune was flattened, as seen in Sect. 3.3.

CFD studies are commonly limited in scope by the size and number of simulations that can be performed with the given resources, and this investigation is no exception. It is important to not over-generalize the results, given that only one boundary condition has been simulated. Previous work shows that the separation length behind dunes changes as a function of free stream velocity [25], and unpredictable changes in the flow field with changing velocity is a general property of turbulent flow [37]. However, to examine the role a foredune plays in mediating mineral dust emissions, the results presented may be combined with experimental data to address a variety of atmospheric scenarios. The experimental work of Gillies et al. [9] demonstrated a linear relationship between wind power density (WPD, \(\hbox {W}\,\hbox {m}^{-2}\)) [41] and \(\hbox {PM}_{{10}}\) at time scales ranging from months to years over a variety of atmospheric conditions at the ODSVRA. This result may be applied to allow a generalized prediction of \(\hbox {PM}_{{10}}\) generated by emissions from a surface using a point measurement. From the simulation data for the unmodified topography, two data points that define the relation between WPD and \(\hbox {PM}_{{10}}\) can be collected from the boundary condition: one for the threshold value predicted using the Bagnold approximation [40] that we assume respresents the initiation of \(\hbox {PM}_{{10}}\) emissions, and one from the boundary condition used to run the simulation. For the simulation over the unmodified topography, the resulting affine function is \(E =.000197 WPD -.630\). This linear relationship may be integrated over observed meteorological data collected at a height of \({10}{400}\,\hbox {m}\) to obtain generalized predictions of \(\hbox {PM}_{10}\) concentration by converting wind speed to wind power density. The flexibility of this method allows for the partitioning of the land surface (e.g., foredune or transverse dune areas), and conversion of \(\hbox {PM}_{10}\) concentration to emissions attributed to different parts of the surface. A full reanalysis to validate this approach is outside the scope of this study, but could be performed at any site where air quality sensors are operational to provide measured values of \(\hbox {PM}_{10}\) and wind speed.

5 Conclusions

Simulations of flow over a coastal nebkha foredune to back dune system were conducted, using high resolution topographic data and measured flow conditions. Additional simulations were performed on modified topographies to isolate the downwind effects of the foredune on shear stress and dust emissions. It was determined that nebkha hummocks in the foredune are effective at providing shelter in their lee where shear is reduced and modulates saltation flux, despite the forms being highly aerodynamic. Higher shear was observed on the transverse dunes as compared to the foredune, although the foredune removes momentum from the wind field before it reaches the transverse dunes. The removal of the foredune was found to increase the area-normalized level of shear on the surface, and reduced downwind sheltering effects. Vegetation was found to play an out-sized role in providing shelter to the foredune, in terms of the reduction of total shear force, dust emissions, and surface area where saltation could occur within the domain studied. The data presented suggest that in terms of a management strategy to mitigate dust emissions from the ODSVRA that impact air quality downwind, the establishment of a foredune is justified for that purpose as it reduces the integrated shear across the domain. As \(\hbox {PM}_{{10}}\) mineral dust emissions scale non-linearly with increasing shear, even relatively modest reductions can have significant effects on the mass emissions and concentrations downwind. Foredunes also likely create additional benefits beyond contributing to improvements in air quality as they provide the opportunity to improve the environment by enhancing diversity of habitats for native plant and animal species.