Holocene evolution of parabolic dunes, White River Badlands, South Dakota, USA, revealed by high-resolution mapping

Abstract

The White River Badlands (WRB) of South Dakota record eolian activity spanning the late Pleistocene through the latest Holocene (21 ka to modern), reflecting the effects of the last glacial period and Holocene climate fluctuations (Holocene Thermal Maximum, Medieval Climate Anomaly, and Little Ice Age). The WRB dune fields are important paleoclimate indicators in an area of the Great Plains with few climate proxies. The goal of this study is to use 1 m/pixel-resolution digital elevation models from drone imagery to distinguish Early to Middle Holocene parabolic dunes from Late Holocene parabolic dunes. Results indicate that relative ages of dunes are distinguished by slope and roughness (terrain ruggedness index). Morphological differences are attributed to postdepositional wind erosion, soil formation, and mass wasting. Early to Middle Holocene and Late Holocene paleowind directions, 324°± 13.1° (N = 7) and 323° ± 3.0° (N = 19), respectively, are similar to the modern wind regime. Results suggest significant landscape resilience to wind erosion, which resulted in preservation of a mosaic of Early and Late Holocene parabolic dunes. Quantification of dune characteristics will help refine the chronology of eolian activity in the WRB, provide insight into drought-driven landscape evolution, and integrate WRB eolian activity in a regional paleoenvironmental context.

INTRODUCTION

Eolian deposits—including sand dunes, sand sheets, and loess—preserve a proxy record of climate change, including drought, on the northern and central Great Plains (Fig. 1) spanning the late Pleistocene to the modern era. Eolian research on the Great Plains has focused primarily on the area's largest dune field, the Nebraska Sand Hills (NSH; Loope et al., 1995; Stokes and Swinehart, 1997; Goble et al., 2004; Mason et al., 2004, 2011; Forman et al., 2005; Miao et al., 2007b; Schmeisser McKean et al., 2015). However, studies of widely scattered dune fields throughout the Great Plains remain relevant as a way to distinguish the effects of local disturbances from regional shifts in hydroclimate related to global-scale climate events, including the Holocene Thermal Maximum (HTM), the Medieval Climate Anomaly (MCA), and the Little Ice Age (LIA). Eolian activity is driven largely by sand availability related to vegetation cover, which is influenced by hydroclimate. While other effects, including fire, anthropogenic activities, animal actions, and other local disturbances to the land surface can trigger eolian activity (Hesp, 2002; Arterburn et al., 2018; Barrineau et al., 2019), multiple lines of evidence indicate that droughts of greater than normal duration—referred to hereafter as “megadroughts” (Woodhouse and Overpeck, 1998; Cook et al., 2016)—have driven eolian activity on a regional scale in the Great Plains during the Late Holocene. Midcontinental megadroughts have been associated with two periods of global climate fluctuation, the MCA (1050 yr to 750 yr) and the LIA (approximately 600 yr to 300 yr), which are identified by climate proxies including eolian activity (Clarke and Rendell, 2003; Goble et al., 2004; Mason et al., 2004; Forman et al., 2005, 2008; Sridhar et al., 2006; Miao et al., 2007b; Halfen et al., 2010; Hanson et al., 2010; Schmeisser et al., 2010; Vimpere et al., 2021), tree rings (Cook et al., 2004, 2007, 2010), lake sediments (Laird et al., 1996; Fritz et al., 2000; Grimm et al., 2011; Hobbs et al., 2011; Schmieder et al., 2011, 2013), and alluvial incision (Burkhart et al., 2008).

Figure 1. North American Great Plains province shown in tan. Modified from U.S. Environmental Protection Agency (2006) and Muhs and Holliday (1995). The inset map shows the Great Plains province in stipple fill with dune fields in tan. Dune fields (light gray) and major streams (blue) of the northern and central Great Plains modified from Muhs and Holliday (1995).

Reviews of the state of understanding of the dune activity and climate change on the Great Plains (Halfen and Johnson, 2013; Halfen et al., 2016) indicate the need for more eolian chronologies using high-quality optically stimulated luminescence (OSL) ages, particularly where climate proxies are scarce, as is the case with the White River Badlands (WRB). Previous studies indicate the WRB dune fields experienced periods of eolian activity coincident with other dune fields on the central and northern Great Plains, supporting climate as the primary driver for eolian activity. WRB dune field experienced an early period of activity during the terminal Pleistocene (21 ka to 12 ka), a period of Early to Middle Holocene activity (~9 ka to 6 ka) during the HTM (~11 ka to 5 ka), and a post–700 yr period during the LIA that extends to the modern era (Rawling et al., 2003; Burkhart et al., 2008; Baldauf et al., 2019). These periods are similar, but not identical, to eolian activity determined for the nearby NSH. Understanding the variability of climate and eolian activity will require more dune field chronologies throughout the Great Plains, as well as a nuanced interpretation of the evolution of eolian landscapes.

The purpose of this study was to leverage previously reported OSL ages (n = 37), together with a newly acquired set of very-high-resolution (1 m/pixel) digital elevation models (DEMs), to identify morphological criteria that could be used to distinguish dunes of Early to Middle Holocene from Late Holocene dunes. Results from previous studies of the WRB dunes (Baldauf et al., 2019) suggest that there is a correlation between dune age, determined from OSL, and dune morphology. Generally, field observations indicate that dune crests with older OSL ages have gentler slopes and rounded contours, whereas dune crests with younger OSL ages have sharp contours and steep slopes, similar to observations from other parabolic dune fields (Patton et al., 2019, 2022a, 2022b). No dune forms were found from the 21 ka to 12 ka eolian event, so these were not considered in this study. Geomorphological characteristics from the Early to Middle and Late Holocene dunes were analyzed, including length to width ratios, relief, slope, and terrain ruggedness index (roughness). Further, soil and vegetation maps (Von Loh et al., 1999; USDA-NRCS, 2022) were compared to dune distribution to determine whether these variables are correlated with dunes of different ages. The end goals of this study include refining the chronology of eolian activity in WRB, creating a model for landscape evolution, and putting eolian activity in the WRB into a regional paleoenvironmental perspective.

STUDY AREA

The WRB are located in the western section of the northern Great Plains, east of the Black Hills uplift in South Dakota, USA (Fig.1). The WRB are defined as the area of deeply eroded strata of late Eocene through Oligocene-age White River Group, bounded by three major E-W oriented streams—the Cheyenne, the Bad, and the White Rivers (from north to south, respectively; Fig. 1). Streams that bound the WRB originate east of the glacial limits of the Rocky Mountain Front and west of the limits of last glacial ice margins, delimited by the Missouri River. A well-developed system of tributaries to the major streams has eroded the prairie into a series of mesas and buttes, known locally as tables (Fig. 2). Large tables range in relief from tens of meters to approximately 100 m and are typically named features on 1:24,000-scale topographic maps of the area. Smaller tables and buttes, known locally as sod tables (Burkhart et al., 2008; Benton et al., 2015), are typically less than 10 m in height and unnamed on topographic maps. Both tables and sod tables are covered with alluvial and eolian deposits of varying thickness, deposited variously on Cretaceous marine strata of the Pierre Shale and Eocene–Oligocene fluvial overbank, lacustrine, and volcanic ash deposits of the White River Group. The presence of fluvial gravels, found at the sediment–bedrock interface of the tables, supports the interpretation that the tables are dissected fluvial terraces of the ancestral White River, Cheyenne River, or their tributaries (Harksen, 1974; Stamm et al., 2013).

Figure 2. Upland tables in the White River Badlands with locations of optically stimulated luminescence (OSL) samples from Baldauf et al. (2019) denoted by diamonds and Burkhart et al. (2008) denoted by circles. Rectangles mark the approximate locations of the study areas. Adapted from U.S. Geological Survey 1:5,000,000 scale surficial deposit map (Soller and Reheis, 2004).

Unconsolidated surficial eolian deposits are mapped on these tables throughout the area between the Cheyenne and White Rivers, from south of the Black Hills uplift to the Missouri River (Raymond et al., 1976; Fig. 1). However, most of the dune fields are located on tables adjacent to the White River from south of the Black Hills uplift to their eastern limit near the 102°W meridian (Fig. 2). Three sets of Quaternary eolian deposits have been identified in the WRB: (1) eolian sand of the sand sheets and dunes in the WRB, which is the focus of this paper; (2) eolian clifftop (ECT) deposits; and (3) the Red Dog Loess. Eolian sand occurs as sand sheets, parabolic dunes, and blowout lobes on upland tables throughout the WRB. Thickness of the eolian sand ranges from less than 1 m to 60 m, increasing in thickness in the downwind (southern) sections of the tables. Dune sand is typically composed of moderately well-sorted sand with less than 10% silt and clay, and CaCO₃ content ranging from less than 3% to 4.5% (Baldauf et al., 2019). Dune chronology published previously (Baldauf et al., 2019) indicates three periods of eolian activity at 21 to 12 ka, 9 to 6 ka, and post–700 yr. Chronology of ECT activity indicates three periods at 4.5 ka, 3.6 ka, and 1000 yr (Rawling et al., 2003). The complete range of ages for the Red Dog Loess remains largely unknown, but published OSL ages indicate loess deposition during the Late Pleistocene (18.5 ± 1.6 ka to 12.7 ± 1.3 ka) (Burkhart et al., 2008).

Three dune fields within the WRB were analyzed in this study, from west to east: Cuny Table, Conata Ranch Table, and Bouquet Table (Fig. 2). Currently, no active dune migration is taking place on any of these tables, and dunes are stabilized by a drought-resistant community of shrubs and grasses. There is no organized drainage on any of the tables, and maximum relief relative to surrounding drainages varies from approximately 100 m at Cuny Table to 30 m at Bouquet Table. Soils present on tables and buttes formed on dunes, sand sheets, and loess plains, and range from moderately well-developed Mollisol–Entisol complexes to poorly developed Entisols (USDA-NRCS, 2022). Vegetation maps for the area (Von Loh et al., 1999) show that the tables in the study area include two vegetation associations: (1) sand sagebrush shrubland and (2) blue grama grasslands. Both communities are associated with Mollisol–Entisol soil complexes (Valent–Valentine–Anselmo associations), but with sand sagebrush found more commonly on dunes, and blue grama grasslands on the sand sheets and loess plains (Von Loh et al., 1999). Lack of water and rugged relief of the table edges limits agricultural use, in most cases, to grazing of domestic livestock. There is no evidence that the tables surveyed in this study have been plowed (Larson and Whitman, 1942), with the exception of Cuny Table, which has been used for row crops since at least the 1940s (USGS).

METHODOLOGY

Photogrammetric methodology using a small unmanned aerial system (sUAS)

For the quantitative analyses of WRB dune morphology, we determined upwind and downwind relief; maximum and mean slope; maximum and mean terrain ruggedness index (TRI), or roughness; and inferred and dominant wind direction. We focused on areas of the WRB dune fields where we collected high-resolution imagery, including the northernmost section of the Conata Ranch Table and the central and southern sections of Bouquet Table. Photogrammetry data were acquired via sUAS (aka “drone”) at two of the three study sites (Bouquet Table and Conata Ranch Table), while data at Cuny Table were not acquired due to logistical constraints. Data were collected using a Dà-Jiāng Innovations (DJI) Phantom 4 Pro quadcopter system with a single gimbal-mounted camera in three separate campaigns (summers of 2017, 2018, and 2019), and for final orthorectified photomosaics and elevation models, final acquisition parameters and equipment were consistent. Multiple flights were necessary for complete coverage due to the large areas, 14.5 km² (5.6 mi²) and 5.1 km² (2.0 mi²) for Bouquet Table and Conata Ranch Table, respectively. Resolution of images ranged from 3.0 to 4.3 cm/pixel.

We collected several large data sets containing overlapping photographs and reconstructed a digital 3D model of the terrain using the structure-from-motion (SfM) process. The main processing algorithm uses the difference of Gaussians pyramid maxima through a scale-invariant feature transform (Lowe, 1999) to estimate the spatial positioning of each pixel in each photograph. Our resulting SfM data are readily converted to orthophoto mosaics (2D projection of the 3D feature), a digital 3D object, or a DEM. We used all three results from the SfM process. Raw elevation and RGB files were imported into ArcMap and projected into a WGS84 datum. For slope and roughness analysis, we downsampled DEMs to 1 m/pixel resolution to reduce the effects of high-frequency changes in elevation and the effects of plant cover. Projected raw elevation files were used to produce hillshade maps with elevations corrected using the World Elevation Terrain layer from ESRI. Projected raw elevation files were also used to create cross sections of selected eolian landforms using an ArcMap spatial analyst tool.

For each landform chosen for analysis, we drew a polygon that captured the approximate extent of the feature (Supplementary Fig. 1), then used the ArcMap 3-D analyst tool to measure length, width, and relief from high-resolution DEMs, using orthoimagery, hillshade maps, and contour maps to aid interpretation. We measured longitudinal distance (length) from the upwind edge or back ridge of the deflation hollow to the downwind crest (or nose) of the depositional lobe (Supplementary Fig. 1). Where curvature was evident in the landform polygon, we curved the longitudinal profile tracks to avoid measuring an apparent feature length. We measured transverse distance as the distance between the crests of each trailing arm of parabolic dunes or blowout, and relief from the upwind and downwind edges or crests of the dunes to the low point in the deflation hollow.

We used the ArcMap slope tool to calculate maximum and mean slope of the polygonal areas outlined for each landform, and ArcPro to determine TRI. ArcPro calculates TRI based on the sum of the elevation change of a cell in a DEM and its eight surrounding cells (Riley et al., 1999). Where dune forms extended outside the limits of the high-resolution DEMs, we calculated slope and TRI only on the part of the polygon within the high-resolution DEM. The resulting population of dune forms includes four Early to Middle Holocene dunes and seven Late Holocene dunes from Conata Ranch Table (N = 11), two Early to Middle Holocene dunes and 10 Late Holocene dunes from Bouquet Table (N = 12), and one Early to Middle Holocene and two Late Holocene dunes from Cuny Table (N = 3; Table 1).

Table 1. Geomorphological and wind data for the White River Badlands (WRB): Cuny Table, Conata Ranch Table, and Bouquet Table. The eolian event category indicates the period when the eolian landform stabilized based on optically stimulated luminescence age, cross-cutting relationships with landforms of known age (Relative), or direct observation (Historic). Mean values are given with 1σ error for slope, terrain ruggedness index (TRI), and wind. Calculations of the summary mean at the bottom of the table group values from each table by age.

Bouquet Table (43°39.62′N 102°16′W) Conata Table (43°41.61′N 102°26.27′W) Cuny Table (43°32.18′N 102°38′W)

Descriptive dune analyses

In addition to the quantitative analyses, we conducted descriptive analyses to classify dune forms. Where available, we used our high-resolution imagery, but outside those areas we used Google Earth, National Agriculture Imagery Program (NAIP) imagery, and USGS 1/3 arcsecond DEMs. Descriptive terms used include categories classifying dunes and blowouts based on length to width ratios, where “lobate” refers to dunes with length to width ratios between 1 and 3, and “elongate” dunes are those with length to width ratios greater than 3.0 (Pye, 1982; Goudie, 2011). We used the term “degraded” in the sense of Pye (1982) to qualitatively describe dunes with subdued relief and rounded dune crests. Terms used for the morphology of parabolic dunes largely follow Wolfe and David (1997), with dune heads divided into a crest and a slipface, with associated trailing arms, a central deflation depression or hollow, and back ridge. Dune head shapes are described as V- or U-shaped (Pye, 1982), and we estimated dune fill as unfilled, partially filled, or filled (Wolfe and David, 1997). Additional terms include “windrift dunes,” parallel linear ridges formed when the head of a hairpin parabolic has been blown out (Pye, 1982; Kilibarda and Blockland, 2011), and “rake-like compound dunes,” composed of the merged heads of multiple dunes (Goudie, 2011). The term “blowout” was used to describe deflation hollows with poorly developed dune heads showing little sign of migration, in contrast to parabolic dunes with well-developed dune heads and signs of significant migration (Pye, 1982).

Wind direction determination

To determine paleowind directions, we measured orientations of dune forms and blowouts of known or inferred age using the high-resolution imagery, in the case of Conata Ranch Table and Bouquet Table, and we used orthoimagery from the NAIP for Cuny Table. A line of best fit for paleowind direction was estimated for selected landforms by bisecting the angle between dune arms or blowout walls. The azimuth of the bisecting line was measured using the ArcMap coordinate geometry (COGO) tool.

RESULTS

Quantitative dune analyses

Quantitative analyses of dune morphology of Early to Middle Holocene (N = 7) and Late Holocene dunes (N = 15) were performed on newly acquired 1 m/pixel DEMs using ESRI tools in ArcMap and ArcPro (Fig. 3a–d, Table 1). The location of each dune and associated polygon is shown in Supplementary Figures 2 and 3, and profiles of longitudinal and transverse elevations are shown for Conata Ranch Table (Supplementary Fig. 4a–k) and Bouquet Table (Supplementary Fig. 5a–j). Dunes used in this analysis were selected based on location within the high-resolution imagery and known OSL stabilization age or crosscutting relationships with dunes of known age. Early to Middle Holocene and Late Holocene dunes could be distinguished at the 1σ error level for maximum and mean TRI and maximum and mean slope. Dunes in this data set met the expected correlation of age with morphology, except for Late Holocene dune forms 8, 9, and 12 from Bouquet Table, which showed lower than expected values for slope and roughness, and Early to Middle Holocene dune form 15 at Conata Ranch Table, which showed higher than expected values. In addition, we constructed seven cross sections approximately orthogonal to inferred paleowind direction connecting points of OSL stabilization ages for Cuny Table (Supplementary Fig. 6), Conata Ranch Table (Supplementary Fig. 7a–c), and Bouquet Table (Supplementary Fig. 8a–e). Cross sections show Late Holocene dunes have sharper crests and greater relief relative to surrounding interdunal areas.

Figure 3. Early to Middle Holocene (EMH) and Late Holocene (LH) plots for mean ruggedness (a), maximum ruggedness (b), mean slope (c), and maximum slope (d). In each box, the horizontal line is the median value, the line cross is the mean, and the whiskers are 1σ from the mean. Boxes represent interquartile ranges (50% of values for each data set).

Wind direction

Results of paleowind direction analysis (Table 1) are grouped by known or relative ages. Locations of estimated wind direction for each dune form are shown in Figure 4a–c. The Early to Middle Holocene data set consisted of eight points (N = 8) of known age, ranging from 10.8 ± 0.80 ka to 6.1 ± 0.50 ka. Paleowind direction estimates ranged from 302° to 340° with a mean of 322° and 1σ error of 12°. The Late Holocene data set consisted of a total of 17 points (N = 17), including 9 points of known age ranging from 720 ± 70 yr to 0 ± 20 yr, 6 points with ages inferred by crosscutting relationships, and 2 points showing evidence of historic activity. Wind direction values ranged from 285° to 337° with a mean 322° and 1σ error of 14°. Results for the older and younger dune periods are indistinguishable at the 1σ error level and similar to the modern wind regime.

Figure 4. Wind direction analysis for (a) Cuny Table, (b) Conata Ranch Table, and (c) Bouquet Table. Dashed arrows are wind direction from Early to Middle Holocene parabolic dunes, and solid arrows are wind direction from the Late Holocene. In each figure, the dashed line is the approximate boundary between Early to Middle Holocene and Late Holocene dunes.

Qualitative dune analyses

Dune morphology shows considerable variation from the upwind to the downwind sections of the tables. The dune field at Bouquet Table is located on the central and southern section of the table. The northern end of the table is a nearly planar area of loess soils (USDA-NRCS, 2022) underlain by as much as 30 m of coarse to fine, sandy silt where exposed on the table edge. From the upwind section of the dune field, high-relief parabolic dunes are unfilled, compound, and nested. Downwind, blowouts and high-relief parabolic dunes are filled, lobate to elongate, and closely spaced. Inliers of degraded parabolic dunes in the upwind section are U-shaped and lobate parabolic, but no degraded dunes are present downwind. Notably, landslides at the southern end of table have transported sand onto the White River floodplain, and where exposed in a landslide scar, sand thickness is as much as 60 m, the maximum thickness of eolian sand observed in the field.

Conata Ranch Table also shows considerable variability in dune style from the northern upwind to the downwind section. The northern end of the table is dominated by partially filled degraded lobate to elongate parabolic dunes and high-relief blowouts and elongate V-shaped parabolic dunes. Downwind, the high-relief parabolic dunes are more complex with more compound, nested dunes and rake-like dunes, and areas of degraded parabolic dunes are interspersed within areas of high-relief dunes. Low-relief windrift parabolic dunes up 1.6 km in length are common.

At Cuny Table there are three dune fields separated by agricultural plains covered with loess or alluvial soils (USDA-NRCS, 2022). We surveyed the easternmost of these dune fields, which consists of a northern section of degraded lobate to elongate parabolic dunes. To the south, dune morphology transitions to closely spaced blowouts and poorly developed parabolic dunes.

DISCUSSION

Drone data acquisition

During three field seasons in the study area, we collected approximately 19.6 km² of drone imagery with a resolution of 3.0 to 4.3 cm/pixel. While this resolution is comparable to LIDAR, the SfM is affected by plant cover, whereas bare earth LIDAR is not. Given that the vegetation communities are low, primarily grasses and shrubs, we downsampled the raw elevation data to 1 m/pixel to approximate bare earth. While methodologies exist for high-accuracy vertical SfM using ground control points (Hugenholtz et al., 2013; Cook, 2017), we did not calibrate our raw elevation files for vertical accuracy. Instead, we corrected our elevations using elevations from USGS 1/3 arcsecond DEMs accessed through ArcMap. Because this study focused on changes in elevation along measured horizontal distances, we believe this did not affect the outcome of calculations.

Quantitative analyses

High-resolution imagery acquired in this study provides an opportunity to refine the existing dune chronology, correlate WRB eolian activity with nearby dune fields (e.g., NSH), and improve correlation with regional and global climate events. Most dunes and blowouts in the WRB can be distinguished as Early to Middle Holocene or Late Holocene based on criteria of mean and maximum slope and mean and maximum roughness (Fig. 3a–d; Table 1), variables that decrease with reduced dune relief related to wind erosion and colluviation (Patton et al., 2019, 2022a, 2022b). While values overlap for size, relief, and length to width ratios of the landforms at Bouquet Table and Conata Ranch Table, slope and roughness factors distinguished the different ages of activity within 1σ error of the mean for each data set with minimal overlap. These criteria correlate the dunes with appropriate OSL age for 6 of the 7 Early to Middle Holocene dune forms and 12 of 15 of the Late Holocene dune forms.

Nonetheless, there were overlapping values within each age set that appeared anomalous. Two Late Holocene dunes from the eastern edge of Bouquet Table, dune forms 8 and 9 (Fig. 3a–d, Table 1, Supplementary Fig. 2), showed lower than predicted mean slope and mean roughness. Late Holocene dune form 12 from Bouquet Table also showed lower than predicted values for all four variables (Fig. 3a–d, Table 1, Supplementary Fig. 3). In the case of dune forms 8 and 9, we interpret this to be the effect of reduced sand availability due to the thinner section of sand that was observed on the eastern side of the table. In the case of dune form 12, it is possible that sand from active dunes upwind, or possibly anomalously rapid colluviation, may have lowered the slope and decreased the roughness. Similarly, Early to Middle Holocene dune form 15 from Conata Ranch Table showed greater than predicted values for all four variables. In this case, we cannot rule out the possibility of reactivation and overprinting of older dune sands with younger sand, where younger sand mantles the morphology of the older feature. Cross sections through points of known age on the tables (Supplementary Figs. 6–8) also support the correlation between steepness of slope and relative age.

We note that values determined for individual dune forms are sensitive to the shape of the polygon created for the form, as investigators may draw boundaries differently. As others have noted, the position of lines and polygons around dune forms is often subjective (Hugenholtz et al., 2012). Further, the calculated values are sensitive to the resolution of the DEM. In this case, we determined that calculated values from the submeter-resolution DEM was influenced strongly by high-frequency changes in elevation with position, due mostly to plant cover, but also likely affected by bioturbation and colluviation. Resolution of 1 m/pixel reduced the noise associated with these variables. Possibly other proxies for roughness, such as dune curvature (Patton et al., 2019, 2022a, 2022b), could be used effectively with this data set for distinguishing dunes by age. Curvature is less sensitive to user subjectivity, as high and low points on dune profiles are more objectively determined than polygon outlines. However, a comparison of these methodologies was not a part of this study.

Morphological analysis and OSL age

Results of this study support the hypothesis that changes with time in dune morphology (e.g., slope and roughness) can be quantified, which can then be used to develop relative age chronologies of eolian activity. In turn, relative age chronologies can be constrained with absolute age data from OSL or ¹⁴C dating. Spatial analysis of eolian landforms on the WRB tables in this study resulted in the identification of two statistically significant categories of dunes based on slope and roughness. These morphological groups correspond with two statically significant clusters of OSL ages from previous work (Baldauf et al., 2019). Similarly, work conducted by Patton and others (2019) found a comparable relationship between absolute age and morphology on the Great Sandy Coast of Queensland, Australia, but on a much larger spatial and temporal scale. While their spatial analysis metrics were different in some respects, their results also showed that changes in eolian landform morphology could be correlated to ages over great distances; in their case, several islands and hundreds of kilometers. Thus, this technique appears to work in this study, and in the case of Patton and others (2019), due to the episodic nature of the activation of eolian landforms.

Nonetheless, while chronologies were determined from spatial analysis in this study, these relative age chronologies cannot replace the need for absolute age information from OSL or ¹⁴C dating. Instead, spatial analysis provides a context for preliminary grouping of landforms by relative age, which can be used for targeted OSL sampling or for the efficient application of other costly or time-consuming field-based techniques. Several key points must be considered during the application of a morphological statistic (such as roughness) and linking it with temporal data from OSL. Most importantly, the ages acquired from dune crests represent the final phase of dune field activity and thus do not constrain the timing of initiation or the duration of the activity. This was recognized by Patton and others (2022b), and they used the OSL dates to represent terminal ages for dune fields, thus linking the temporal data with morphological statistics.

Wind direction

Wind direction determination from Early to Middle Holocene and Late Holocene dune forms (Fig. 4a–c, Table 1) shows dominant wind directions (ā = 322° ± 12° and ā = 322° ± 14°, respectively) that are indistinguishable at the 1σ error level. This direction is consistent with the modern winter dominant wind direction and modern potential sand drift direction calculated in other studies (Muhs et al., 1997; Schmeisser et al., 2010). WRB dune forms do not preserve evidence of Late Holocene SW winds, coincident with the MCA, responsible for linear dune development in central NSH (Sridhar et al., 2006). However, many surrounding dune fields lack a SW wind component in sand transport, including dune fields in northern Colorado, the western section of the NSH, and dune fields in North Dakota (Schmeisser et al., 2010). This suggests that either the WRB dune fields were outside the influence of SW winds during the Late Holocene or the SW winds during this period were not strong enough to transport sand in the WRB. Alternatively, if dunes developed in the WRB during the MCA, they were reworked completely by subsequent Late Holocene eolian activity, which cannot be ruled out using a surficial data set. Detailed mapping of dune field subsurface architecture and associated paleowind indicators would be required.

Topographic profiles through WRB dune fields, constructed using high-resolution imagery collected in this study, were used to determine crosscutting relationship and support paleowind analysis (Supplementary Figs. 6–8). Neither the depth profiles of the cored and augured samples nor the cross-stratified sediments found at exposure outcrops were sufficient to support a 3D analysis of Holocene dune forms. Nonetheless, geomorphological analysis of cross sections through points of known age in the WRB dune fields supports the distinction between the degraded Early to Middle Holocene dunes and the higher-relief Late Holocene dunes.

Dune field dynamics

Based on quantitative and qualitative analyses, we believe that upwind sections of the dune fields have remained stable, preserving a degraded Early to Middle Holocene eolian landscape, while downwind areas have been broadly reactivated in the Late Holocene. Dunes with shallow slopes and low roughness are found in the upwind sections of the three dune fields considered in this study. Three variables likely to be responsible for this preservation are (1) variations in sand availability, (2) soil development, and (3) stabilizing vegetation.

Sand availability at the northern end of the tables in the study area is limited by the relief of the tables relative to potential upwind sand sources (Fig. 2). We infer that the sand source for the dunes was the floodplain of the ancestral White River, and that the sand was deposited on the table before terrace incision by the White River and its tributaries, no later than 12 ka, the end of the first eolian event detected in the study area. No dune forms have been detected from this terminal Pleistocene event, but sand of this age, some of it cross-bedded, is exposed near the base of the eolian sequence on cliff edges and road cuts (Rawling et al., 2003; Burkhart et al., 2008; Baldauf et al., 2019). To the north, the Cheyenne River incised its channel during this same period, in response to stream piracy, a shift from glacial to interglacial climate, and isostatic rebound from the retreating Laurentide Ice Sheet (Stamm et al., 2013; Stetler and Hazelwood, 2020). The White River, like the Cheyenne, is a tributary to the Missouri River (Fig. 1), and likely incised its channel during this same period, although the timing of incision has yet to be resolved. Based on the elevation of gravels exposed on the cliff edges of the tables, the incision of the White River and its tributaries has been 20 to 30 m, comparable to incision seen on the Cheyenne River.

Dune activity that followed the late Pleistocene occurred as the tables were increasingly isolated from upwind sand sources and from adjacent dune fields. We found little evidence of sand migration onto the tables at the upwind end, with the exception of the formation of ECT deposits (Rawling et al., 2003), which represent a relatively small volume of eolian material. On the other hand, we found abundant evidence of thicker downwind eolian deposits and sand migration onto the White River floodplain to the south. Furthermore, increasingly complex and filled parabolic dune heads downwind suggest an increase in sand availability.

Soil development is strongly correlated with dune age, supporting its role in stabilization and preservation of eolian surfaces in the study area. Soil maps of the WRB (USDA-NRCS, 2022) show thicker, better developed soils in areas with older dunes. Middle Holocene dune forms on Cuny Table (Fig. 4a), for example, are found within the zone mapped with a Jeyem–Valentine soil complex, including soils with moderately developed B horizons, while downwind areas of Late Holocene dunes are mapped as Valentine sands, thin soils with no B horizon (USDA-NRCS, 2022). This relationship holds true for Bouquet Table as well, where an Anselmo–Valentine complex with moderate B-horizon development is associated with Early to Middle Holocene dunes, while downwind Valentine sands are mapped in areas dominated by Late Holocene dunes.

Similarly, grassland community distribution and dune age are correlated in the study area. Early to Middle Holocene dunes are mapped with a blue grama grasslands plant community, which is found in sandy soils on lower slopes (Von Loh et al., 1999). Higher-relief, Late Holocene dunes are mapped with sand sagebrush and prairie sandreed communities, which are associated with sandy soils and steeper slopes. Of the two communities, blue grama offers the higher vegetative cover, which may have acted to stabilize older dunes during Late Holocene droughts. Both communities experienced significant periods of drought during the 1930s Dust Bowl without reactivation of dunes in the WRB (Baldauf et al., 2019). Taken together, soil and vegetation may operate as reinforcing feedback for stabilization of the eolian surfaces; as slopes decrease with colluviation and deflation, soils develop and vegetation cover increases, increasing the landscape resistance to eolian reactivation and allowing further soil development. Reduced sand availability upwind reduces the chance of disturbance of the vegetative cover by migrating sand.

Paleoenvironmental interpretation

Early to Middle Holocene

Periods of eolian activity associated with dune sand in the WRB at 21 ka to 12 ka, 9 ka to 6 ka, and post–700 yr coincide broadly with periods of activity in other dune fields on the northern and central Great Plains, as well as with periods of drought identified by other climate proxies. The oldest eolian period at 21 ka to 12 ka coincides with the terminal Pleistocene and the end of the Younger Dryas in the Rocky Mountains (Gosse et al., 1995; Pierce, 2003; Briner, 2009). In the study area, dune forms from this period are not preserved, but sand from this age is preserved in cliff exposures and road cuts throughout the WRB (Rawling et al., 2003; Burkhart et al., 2008; Baldauf et al., 2019). During this period, climate transition from glacial to interglacial triggered stream reorganization, incision, and terrace formation (Stamm et al., 2013; Stetler and Hazelwood, 2020). Significant incision along the Cheyenne River is documented from this period (Stamm et al., 2013; Stetler and Hazelwood, 2020), and it is likely that significant incision of the White River and the formation of the tables occurred during this period as well, lowering water tables and cutting the tables off from upwind sources of sand. During this period, the Red Dog Loess was deposited on many tables throughout the study area, beginning no later than 18 ka (Burkhart et al., 2008), which falls within the period of Peoria Loess deposition in central Nebraska and western Iowa (Aleinikoff et al., 2008; Muhs et al., 2008; Yang et al., 2017).

The period following 12 ka coincides with a period of higher effective moisture on the Great Plains and soil formation (Brady Soil) following the deposition of the Peoria Loess (Miao et al., 2007a; Muhs et al., 2008; Tecsa et al., 2020). There is little evidence of dune activity in the WRB between 12 ka and 9 ka, however this period is marked by loess accumulation and soil formation on surrounding tables (Rawling et al., 2003; Burkhart et al., 2008; Baldauf et al., 2021). WRB dune activity resumes in the Early to Middle Holocene at ~9 ka and ends by 6 ka, approximately coincident with the HTM (11 ka to 5 ka; Fischer et al., 2018). The HTM coincides broadly with ecological changes indicative of hemispheric climate change (Fischer et al., 2018) and lower effective moisture on the Great Plains (Laird et al., 1996; Meltzer, 1999; Dean and Schwalb, 2000; Williams et al., 2009; Grimm et al., 2011; Davies et al., 2019; Wendt et al., 2022). During this eolian period, soils that formed on tables were buried by loess and ECT, and dune forms from the terminal Pleistocene were reworked. Early to Middle Holocene parabolic dunes survive as degraded dunes in upwind sections of the tables and in isolated areas within areas of Late Holocene dunes (Fig. 4b and c).

Late Holocene

While climate generally returned to a period of higher effective moisture, the Middle to Late Holocene was punctuated by periods of extreme drought based on eolian activity in the WRB and other dune fields on the Great Plains. Periods of eolian activity at 3.8 ka, 2.5 ka, and 1000 yr are correlated to drought according to lake sediment studies (Schmieder et al., 2011, 2013) and tree ring data (Cook et al., 2011), which indicate that these Late Holocene droughts were short-duration events in a period of otherwise higher effective moisture. Deposition of ECT on table edges throughout the WRB during this period (Rawling et al., 2003) suggests that these droughts affected the study area, but there is no evidence of dune reactivation until approximately 700 yr, at the end of the MCA. The LIA is a period of dune reactivation in other dune fields in western sections of the northern and central Great Plains (Arbogast and Johnson, 1998; Clarke and Rendell, 2003; Forman et al., 2008; Halfen et al., 2010, 2012; Werner et al., 2011), coinciding with periods of drought determined from tree ring data (Cook et al., 2010), lake sediments (Fritz et al., 1994), and alluvial incision (Burkhart et al., 2008). Thus, it is possible that the resilient grasslands communities in the study area resisted drought reactivation through the end of the MCA, then collapsed during drought conditions in the western section of the northern Great Plains during the LIA. Delayed dune reactivation during drought periods and delayed recovery to stabilization have been described in other dune fields (Werner et al., 2011; Xu et al., 2020).

CONCLUSION

In an area of the Great Plains with few available climate proxies, the WRB dune field activity spans regional, hemispheric, and global climate events, including the Pleistocene–Holocene transition, HTM, MCA, and LIA. Geomorphological analysis of a new database of high-resolution imagery has increased our understanding of eolian processes and dune field evolution. From our results, we draw four conclusions.

First, drone imagery and SfM DEMs have proven to be effective tools for studying dune evolution, especially the determination of relative age. We distinguished Early to Middle Holocene parabolic dunes from Late Holocene dunes by their shallow slopes and their low roughness (TRI), as determined from 1 m/pixel very-high-resolution imagery.

Second, WRB dune fields preserve a mosaic of older and younger eolian surfaces. Early to Middle Holocene parabolic dunes are preserved in upwind areas of the dune fields and in isolated areas of Late Holocene dunes.

Third, enhanced preservation of the Early to Middle Holocene landscape results from reduced sand availability, better soil development, and denser vegetative cover.

Fourth, wind directions determined from Early to Middle Holocene and Late Holocene dune forms are indistinguishable at 1σ error level and are similar to modern seasonal NW dominant wind.

We propose a general evolutionary sequence, from stability to reactivation, to explain the mosaic of older and younger eolian landscapes present on the tables in the WRB. During periods of moderate mesic climate, soil formation stabilizes eolian surfaces on tables and reduces wind erosion potential. Stability and better developed soils allow grama grass communities to succeed sage communities, increasing the density of vegetative cover, promoting soil development, and further increasing resistance to erosion. With a shift to xeric conditions during megadrought events, collapsing grass communities allow a strong seasonal NW wind to erode blowouts in table surfaces. This sand migration further destabilizes downwind grass communities, leading to broader eolian reactivation of table surfaces and formation of parabolic dunes. Over time, upwind areas become increasingly resistant to wind erosion through continued soil development and high-density vegetative cover. Downwind areas remain vulnerable to reactivation due to the steeper slopes of the parabolic dunes, thinner soils, and less dense vegetative cover. Nonetheless, isolated areas of older dune forms remain in downwind areas due to incomplete reworking of the table surface.

We envision two future directions of research in the WRB dune fields. First, we plan to refine the Late Holocene eolian chronology. To do this, we will target sampling of elongate parabolic dunes, from back ridge to dune head, to determine ages from initiation to stabilization. We hope to discriminate whether activity began as early as the MCA or was confined to the LIA. Second, we have begun an investigation of the Red Dog Loess, which is a last glacial period loess. The stratigraphy of this loess, the timing of its deposition, and its relationship to the Peoria Loess to the south are poorly understood. Putting the Red Dog Loess in an accurate age and stratigraphic context through high-resolution mapping and OSL sampling promises to shed light upon the paleoenvironment in this section of the northern Great Plains during the Pleistocene–Holocene transition.