Windblown dunes are common features in our solar system, forming on planetary surfaces that span wide ranges in gravity and both atmospheric and sediment properties. The patterns formed by their crests, which are readily visible from orbital images, can record information about recent changes in boundary conditions, such as shifts in wind regime or varying sediment availability. Here, we demonstrate that the density of dune interactions (where neighboring crestlines are close to each other) within a dune field is an indicator of such changes. Using orbiter-based images of 46 dune fields on Earth and Mars, we compiled a database of pattern parameters including dune spacing, crestline orientation, and interaction density. Combined with sediment fluxes derived from ERA5-Land data and a martian global circulation model, we also compiled dune turnover time scales (the time it takes for a dune to migrate one dune length) for each investigated dune field. First, we show that dune fields undergoing changes in boundary conditions display higher than expected dimensionless interaction indices. Second, dune fields with longer turnover times display a wider range in interaction indices on both Earth and Mars because they are more likely to be observed while still adjusting to recent changes in boundary conditions. Thus, a dune field’s interaction index offers a novel tool to detect and possibly quantify recent environmental change on planetary surfaces.

Windblown sediment self-organizes into bedforms, generating patterns that can be seen from orbit; these patterns have been reported on Venus, Earth, Mars, Titan, Io, Pluto, and the 67P/Churyumov-Gerasimenko comet (Fig. 1; Bourke et al., 2010; Thomas et al., 2015; Telfer et al., 2018; McDonald et al., 2022). Because dunes and their patterns form and evolve in concert with wind flow, their abundance in the solar system offers a unique opportunity to better understand interactions among planetary surfaces, gas flows (winds), and other near-surface volatiles over a wide range of boundary conditions, with important implications for interpretations of planetary surface processes, atmospheric circulation, and climate (Ewing et al., 2015; Lapôtre et al., 2020). Such boundary conditions control both dune morphology and pattern evolution (Kocurek et al., 2010; Ewing and Kocurek, 2010).

The crestline orientations of large dunes reflect their net drift integrated over relatively long time scales (often >1 k.y.; e.g., Gunn, 2023), such that dune-field pattern analyses may provide important insights into boundary conditions in the recent past, including changes in dune-forming winds over a range of time scales (individual storms to Milankovitch cycles). As dunes adjust to changing boundary conditions, e.g., by altering their orientation or wavelength or developing smaller superimposed dunes (e.g., Kocurek and Ewing, 2005), they interact with each other through exchanges of sediment, typically enhanced near dune defects (i.e., terminations and bifurcations; Werner and Kocurek, 1997; Ewing et al., 2015; Day, 2022).

The spatial densities of both dune defects (Ewing et al., 2006) and interactions (Day and Kocurek, 2018) have been shown to decrease with pattern coarsening (i.e., increase in dune size), and thus with dune-field maturity. Day and Kocurek (2018), in turn, showed that the interaction index (I*), which is a dimensionless form of interaction density (I), is a self-similar quantity; i.e., it remains statistically unchanged regardless of dune scale. However, the influences of boundary conditions and internal dynamics in controlling dune-pattern evolution remain elusive, challenging the usefulness of interaction or defect density as a proxy for time.

Here, we test the hypothesis that I reflects the maturity of dune-field patterns. Specifically, we demonstrate that a dune field’s adjustment to changes in boundary conditions results in transiently increased I (Fig. 2A). Furthermore, because dunes that take a longer time to adjust to boundary condition changes are more likely to be observed deviating from the expected pattern, dune fields with longer turnover time scales (Tt) would be expected to display a wider range in I on both Earth and Mars.

Site Selection

To test our hypothesis, we quantified dune-field patterns within 46 dune fields on Earth and Mars and constrained their time scales of pattern adjustment. Dunes were classified based on crestline morphology (linear or crescentic), as assessed visually from crestline sinuosity and lateral continuity. In our database, linear dunes were all longitudinal to oblique, whereas crescentic dunes encompassed transverse dune types such as barchans, laterally linked barchans, and crescentic ridges (Rubin and Hunter, 1987; Courrech du Pont et al., 2014; Gao et al., 2015). Terrestrial dune fields were predominantly selected from sites investigated by Day and Kocurek (2018) and spanned four continents (Fig. 1A). Martian dune fields were selected from within and outside of impact craters and spanned equatorial to polar latitudes (Fig. 1B). We also sought to include regions where dunes have been subjected to known perturbations like human-induced dune flattening (Fig. 3A; Ping et al., 2014; Lü et al., 2021, 2022), topographic confinement (Fig. 3C), or frost and local changes in wind conditions (Fig. 3E; Ward et al., 1985).

Dune-Field Digitization and Pattern Parameters

All terrestrial dune fields except the Tengger Desert site (Fig. 3A; Google Earth Pro historical images) were mapped from ESRI satellite imagery. Context Camera (CTX; Malin et al., 2007) tiles from the global mosaic of Dickson et al. (2018) were used for martian dune fields (Supplemental Material1). A counting circle of surface area A was drawn within each dune field, in which all dune crestlines were traced, and mean dune spacing (λ), crestline azimuth (δc), slipface length (Ls, used to determine dune height for crescentic dunes), and dune length (Ld) were measured (Supplemental Material). Dune interactions, defined as the locations in A where two neighboring crestlines were traced within a distance <10% of the field-averaged λ from one another (Day and Kocurek, 2018), and their number (n) was counted within A. Dune interaction density, I, was then defined after Day and Kocurek (2018) as

formula

and the nondimensional interaction index, I*, was defined as

formula

A sensitivity analysis was conducted to investigate the influence of circle size, placement, and the threshold distance used to define an interaction and estimate typical uncertainty in I and I* (Supplemental Material).

Turnover Time Scales

A dune’s turnover time scale, Tt, was defined as the time it takes to migrate one dune length (Allen, 1974, 1976; Myrow et al., 2018) under a given average sand flux, . For a given sediment porosity (φ),

formula

where dune height (H) and length (Ld) were derived from elevation data and images, respectively (Supplemental Material).

Sediment flux was estimated from ERA5-Land (Muñoz Sabater, 2019) for Earth and a MarsWRF simulation (Richardson et al., 2007; Supplemental Material) for Mars. The 10 m instantaneous wind components were used to calculate wind shear velocity (u*; Supplemental Material) and wind azimuth (δw). The saturated sediment mass flux over a flat bed, q, was calculated as

formula

where γ = 5, ρf is atmospheric density, g is the acceleration of gravity, and u*,it is the impact threshold shear velocity (Martin and Kok, 2017). We employed a constant sediment density (quartz) and air density for Earth, and a constant sediment density (basalt) and atmospheric density from the MarsWRF simulation for Mars (Supplemental Material). Because Equation 4 was calibrated for terrestrial conditions (Martin and Kok, 2017), we empirically corrected Mars fluxes with those observed from repeat orbiter-based images (Bridges et al., 2012; Supplemental Material).

Saturated mass flux (Eq. 4) was then converted into a volumetric flux and projected in the crest-normal direction (qm). The flux at the crest, qm,c, was estimated from qm by including the acceleration of winds up dune slopes as qm,c = qm(1+βS), where β ~9.4 is a speed-up factor (Jackson and Hunt, 1975; Courrech du Pont et al., 2014; Gunn, 2023), and S is bed slope upwind of the crest as measured in the flux direction. Assuming the dune-toe flux to be ~0 (Supplemental Material; Gunn et al., 2022; Gunn, 2023), the average sediment flux over a dune could then be estimated as

formula

Despite large differences in boundary conditions, linear and crescentic dunes on Earth and Mars display the same overall decreasing trend in I with increasing λ (Fig. 2B). That trend roughly follows a power-law relationship,

formula

where a is a constant and b ≈ –2.0, consistent with the findings of Day and Kocurek (2018). However, in contrast with previous findings, linear and crescentic dunes do not appear to follow distinct trends (i.e., they do not have distinct a values). We attribute this difference from the results of Day and Kocurek (2018) to our refined methodology and, notably, to the independence of our counting surface area on λ and dune type (Supplemental Material).

First, we investigate how dune patterns respond to known perturbations in boundary conditions. In 2014, a part of the Tengger Desert in China was flattened (Fig. 3A) to test dune formation theory in a landscape-scale experiment (Ping et al., 2014; Lü et al., 2021, 2022). Applying our methodology to five images collected from 2016 to 2022, we observed a relatively high I* within emerging dunes, followed by a consistent decrease in I* as the pattern coarsened over time (Fig. 3B). Next, we investigated the case of dunes migrating through a wind gap in the Namib Desert (Fig. 3C). There, topographic confinement locally induces wind and sediment convergence and forces migrating dunes together, increasing I*. As dunes emerge from the valley, they readjust to their unconfined conditions, and I* decreases again (Fig. 3D). Finally, a set of barchan dunes in the north polar sand sea of Mars migrates east within a largely cyclonic ring of dunes driven by westerly winds, but local perturbation by northerly winds (likely driven by thermals off Chasma Boreale; Ward et al., 1985) and an apparent local increase in surface frost cause collision and merging into connected dunes (similar to dunes at White Sands National Park in response to local water table–controlled variations in sediment availability and internal boundary-layer effects; Kocurek et al., 2010; Ewing and Kocurek, 2010; Jerolmack et al., 2012). Those dunes separate into isolated barchan dunes farther downwind to the west (Fig. 3E). Here, too, I* increases before the pattern adjusts again to a lower I* downwind (Fig. 3F).

Second, we calculated Tt for all 46 investigated dune fields; these values range over five orders of magnitude (Tt ≈ 1.7–43,000 Earth years) and are consistent with independent values are available (Supplemental Material). We demonstrate that I* does not follow a unique trend with Tt. Rather, the lowest I* value for a given Tt decreases with Tt (inherently from geometry), whereas the largest I* value for a given Tt increases with Tt (Fig. 4A), resulting in an increasing standard deviation in I* with increasing Tt (Fig. 4B).

Dunes on Earth and Mars follow the same decreasing trend in I with increasing λ (Fig. 2B), regardless of dune type. This observation supports the hypothesis that universal dynamics, likely controlled by geometry rather than process, govern steady-state dune-pattern statistics as dunes grow (Eq. 6). Scatter around that overall trend can be explained by recent changes in boundary conditions among investigated dune fields. Support for this interpretation is twofold.

First, we show that dune fields undergoing changes such as human-made perturbations (Fig. 3A), topographic confinement (Fig. 3C), and local changes in wind conditions (Fig. 3E) record transient spatiotemporal variations in I*. Second, because dunes with long Tt take a longer time to adjust to new boundary conditions (Fig. 2A), there should be a wider range in I* values among dune fields with high Tt at any given time. Consistent with this expectation, we observe a systematically higher variance in I* with increasing Tt (Fig. 4). Altogether, our findings suggest that I* is an indicator of spatiotemporal changes in boundary conditions.

Several dune fields display low I* values, but no dune field was found to have a total absence of dune interactions (Table S3). This observation supports the existence of a nonzero background level of interactions, I0 (Fig. 2A), which arises from geometry, autogenic dune behavior, and local perturbations. Whereas dune fields with the lowest I for a given dune size (Fig. 2B; or lowest I* for a given Tt; Fig. 4A) ought to have mature patterns (i.e., I ≈ I0, I*I*,0), some boundary conditions might, in principle, lead to I0 values above that minimum. In turn, dune fields defining the upper bound in I* for a given Tt (Fig. 4A) display the most immature patterns of our observed sites (e.g., time tmax in Fig. 2A).

Most investigated dune fields display intermediate I* values (Fig. 4A), either indicating relatively high I*0 or, more likely, that they were still adjusting to recent perturbations (Myrow et al., 2018). Dune fields with * > 2.0λ–0.76 (upper limit of the 95% confidence interval for the power-law fit to the lower boundary; Fig. 2B; Table S3) deviate the most from expected pattern statistics and are thus most likely to indicate recent changes in boundary conditions.

Though not tested here, the same quantitative framework may apply to subaqueous bedforms. The interaction density of fluvial and eolian dunes follows the same quantitative relationship with dune spacing (Mason et al., 2020); furthermore, smaller dune sizes make fluvial dune patterns more likely to be observed in a mature state, as supported by their close proximity to the I0 line derived from eolian dunes (Figs. 2B and 4A). Therefore, we propose that I* can be used to detect changes in boundary conditions across environments if rudimentary constraints on sediment flux and bedform dimensions are available. In the case of eolian dunes, the latter can be readily gathered from satellite images and climate models, opening new avenues to search for hints of recent environmental change from dune-field patterns on planetary surfaces.

We show that, unlike previously suggested, dune interaction density follows a universal relationship with dune spacing that is independent of dune type on Earth and Mars. Through a combined analysis of spatiotemporal changes in pattern statistics where known perturbations have occurred and over 46 dune fields on two different planets, we attribute scatter around this universal trend to pattern adjustment in response to recent changes in boundary conditions within individual dune fields. Restated, under dynamic steady-state boundary conditions, dune-field patterns mature through coarsening and a decrease in I, but with sufficient perturbation of boundary conditions, the system responds by pattern regeneration, reflected by a transient increase in dune interaction density. Thus, as long as depositional conditions persist, dune fields recently subjected to changes in boundary conditions display an abnormally high I* for their size (or Tt). To date, knowledge of environmental conditions is limited to nonexistent on most planetary surfaces in the solar system. Therefore, dune interaction density provides a powerful new tool to identify the signatures of recent or local changes on other planetary bodies and, possibly, across transport environments.

1Supplemental Material. Supplemental text describing methods for dune-field digitization, digitization sensitivity analysis, sediment flux calculations, and turnover time scales; validation for turnover time scales; terrestrial wind data sensitivity analysis; three supplemental tables and five figures; and all compiled data. Please visit https://doi.org/10.1130/GEOL.S.23631873 to access the supplemental material, and contact editing@geosociety.org with any questions. ERA5-Land data are available from the Copernicus Climate Data Store (https://cds.climate.copernicus.eu/).

We thank editor Andrew Barth, Gary Kocurek, Paul Myrow, and two anonymous reviewers, whose comments improved our manuscript. A National Center for Airborne Laser Mapping (NCALM) seed grant was provided to Marvin. This work was partly supported by the National Aeronautics and Space Administration (NASA) under grant 80NSSC20K0145 to Lapôtre.

Gold Open Access: This paper is published under the terms of the CC-BY license.