Dune interactions record changes in boundary conditions

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 (T_t) would be expected to display a wider range in I on both Earth and Mars.

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

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 Material¹). 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 (L_s, used to determine dune height for crescentic dunes), and dune length (L_d) 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

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

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

A dune’s turnover time scale, T_t, 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, q̄. For a given sediment porosity (φ),

where dune height (H) and length (L_d) 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

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 (q_m). The flux at the crest, q_m,c, was estimated from q_m by including the acceleration of winds up dune slopes as q_m,c = q_m(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

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,

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 T_t for all 46 investigated dune fields; these values range over five orders of magnitude (T_t ≈ 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 T_t. Rather, the lowest I_* value for a given T_t decreases with T_t (inherently from geometry), whereas the largest I_* value for a given T_t increases with T_t (Fig. 4A), resulting in an increasing standard deviation in I_* with increasing T_t (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 T_t 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 T_t at any given time. Consistent with this expectation, we observe a systematically higher variance in I_* with increasing T_t (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, I₀ (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 T_t; Fig. 4A) ought to have mature patterns (i.e., I ≈ I₀, I_* ≈ I_*,0), some boundary conditions might, in principle, lead to I₀ values above that minimum. In turn, dune fields defining the upper bound in I_* for a given T_t (Fig. 4A) display the most immature patterns of our observed sites (e.g., time t_max 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 I₀ 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 T_t). 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.

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.