A distinct ripple-formation regime on Mars revealed by the morphometrics of barchan dunes

Sand mobilized by wind forms decimeter-scale impact ripples and decameter-scale or larger dunes on Earth and Mars. In addition to those two bedform scales, orbital and in situ images revealed a third distinct class of larger meter-scale ripples on Mars. Since their discovery, two main hypotheses have been proposed to explain the formation of large martian ripples—that they originate from the growth in wavelength and height of decimeter-scale ripples or that they arise from the same hydrodynamic instability as windblown dunes or subaqueous bedforms instead. Here we provide evidence that large martian ripples form from the same hydrodynamic instability as windblown dunes and subaqueous ripples. Using an artificial neural network, we characterize the morphometrics of over a million isolated barchan dunes on Mars and analyze how their size and shape vary across Mars’ surface. We find that the size of Mars’ smallest dunes decreases with increasing atmospheric density with a power-law exponent predicted by hydrodynamic theory, similarly to meter-size ripples, tightly bounding a forbidden range in bedform sizes. Our results provide key evidence for a unifying model for the formation of subaqueous and windblown bedforms on planetary surfaces, offering a new quantitative tool to decipher Mars’ atmospheric evolution.


Measurements of barchan morphometrics from dune contours
After extracting the outlines of dunes using machine learning 1 , dunes morphometrics were measured by identifying six reference points along their outlines -the center of the slipface, the apexes of horns, the tail (edge of the dune in the upwind direction) and the sides. To measure the morphometrics of dunes, we developed the following pipeline: (a) The dune slipface center is found as the deepest convexity defect of the contour ( Figure  S1a,b). (b) The apexes of the horns are identified as the intersection points between the convex hull of the dune and the contour, closest to the slipface (red outline in Figure S1c,d). (c) After identification of the slipface center and horn apexes, the dune is rotated such that its migration direction faces "up". When the dune is symmetric (i.e., with long-to-short horn lengths ratio < 1.5), the migration direction is estimated as the bisector to the vectors connecting the slipface center to the horns' apexes. When the dune is asymmetric (i.e., with long-to-short horn lengths ratio > 1.5), its migration direction is approximated as the direction of the vector connecting the dune tail to the slipface center. The tail of the dune is identified as the point furthest away from the horns in the lower (downwind) half of the dune. The dune sides are calculated as the two extreme points in the direction perpendicular to the migration direction. (d) The height of the dune is measured automatically from the planview-projected length of the slipface (Figure S1e,f). (i) To measure the planview length of the slipface, we first automatically extract a profile of the pixel intensity of the CTX image along a line passing through the slipface center and perpendicular to the dune contour.
The profile is smoothed using a gaussian filter with window size of 2 pixel and calculate its absolute value (to only have positive smoothed grayscale values along the profile). (iii) We use a peak detection algorithm to find the width of the most prominent extremum with respect to the median brightness of pixels around it ( Figure S1e). (iv) The error in our method of measurement is at least 5 meters, the resolution of a CTX pixel. However, effects such as sand aprons and illumination may exacerbate the error -especially in high latitudes. See Figure S9.
Supplementary Figure S1: (a,b) Detecting dune slipfaces as the deepest convexity defect. Red contour and points shows the convex hull of the contour (c,d) Detecting the horns' apexes as the intersection points between the contour and its convex hull. Yellow arrows are for demonstrating the algorithm iteration along the dune contour (e,f) Inferring dune height from the planview length of the slipface (in map view). The width of the slipface is calculated as the width of the most prominent peak in pixel brightness along a profile passing through the slipface center perpendicular to the dune contour. Image derived from the global CTX catalog (-42.20°N, -31.89°E).

Uncertainty of derived scaling laws and mitigation
A full discussion of potential errors of our machine learning approach was previously published in an analysis of detection performance 1 . To compute the orthogonal distance regression parameters corresponding to the regression in Figure 2, we employed the Python Scipy ODR class, which is based on the ODRPACK Fortran implementation of the Levenberg-Marquardt algorithm 2 . As an estimate of the error on our scaling law parameters derived using orthogonal regression, we use the standard deviation of the fitted parameters (slope and intercept in log-log scale; see discussion in 3 ). For our orthogonal regression lines, we calculate ! = 1 − / , where RSS and TSS are the residual and total sum of squares. For the log-transformed relationship in Figure 3 and Supplementary Figures we use linear regression (assuming density is not a response variable, as it is measured from topography or by the GCM). We employ the Python Scipy linregress function and estimate the errors as the standard errors of the fitted parameters.
Additional examples of dunes detected using the machine learning model  Figure S1). North is up. Each panel describes a 800x800 pixel segment cropped out of the global CTX mosaic 4 , as described in the Methods. Data was subsequently filtered as described in the main text to remove spurious detections. Scale bars represent 1 km. North is up in all images.

Validation of morphometric database with manually measurements
To validate individual dune measurements, we compare manual measurements made in JMARS to automatic measurements made using our algorithm for 10 random dune fields on Mars. Figure  S3 illustrates one of these comparisons for the dune field shown in Figure 2 (-42.20°N, -31.89°E).
Automatically detected morphometrics are further compared with manual measurements of barchan dunes on Mars (Sherman et al., 2021). The best-fit power law describing the planview dimensions of automatically detected barchans is close to that describing the planview dimensions of manually detected barchans (blue and orange lines in Figure S4).

Validation of dune-height measurements with HiRISE stereographic DTM
The correlation between barchan width and height is weaker than between length and height (see R 2 values in Figure 2), which may partially result from the higher uncertainty associated with measuring width perpendicular to dune migration direction. In addition, our assumptions that the dune brink is the highest part of the dune and that the slipface is close to the angle of repose may not always be accurate, potentially increasing data scatter.
We validate our height-determination algorithm by manually measuring the heights of 22 barchan dunes in Oyama Crater using a HiRISE digital elevation model 5 (image ESP_042714_2035) of a dune field centered around 23.23°N, -20.4°E. Vertical precision is < 1 m ( Figure S5).
Supplementary Figure S5. Comparison between automatically estimated dune heights and dune heights as manually measured from a HiRISE digital elevation model (ESP_042714_2035). Dashed line highlights a 1:1 relationship. The 'x' symbol (bottom right corner) denotes a dune whose contours were not adequately detected by the machine learning algorithm.

Compilation of terrestrial barchan dune morphometrics
Our morphometric dataset of terrestrial barchan dunes (Figure 2) was compiled from previously published studies with attention to providing a diversity of dune-field settings and geography (  Peru (N = 45) 10 Length was measured from the tail to slipface Brazil (N = 47) 11 Length was measured as the length of the windward side California (N = 42) 12,13 In 10 , length was measured as the length of the stoss side Kuwait (N = 10) 14 Length was derived by adding the slipface and the stoss lengths Multiple locations, including Algeria, Yemen, Saudi Arabia, Western Sahara, Chad, Peru, China, Namibia, Sudan, Iran, and Niger (N = 2686) 15 Length was measured from the tail to slipface Table S3. Locations and references for compiled terrestrial barchan-dune morphometrics.

Correlation between atmospheric density and dune width
To further validate the result presented in Figure 3, we additionally compute the power law relationship between dune width and atmospheric density. As expected, since the relationship between width and length is close to linear, the power law exponent relating dune width and atmospheric density is, too, close to -2/3 ( Figure S6). Figure S6: Decrease in dune width (blue-to-yellow circles) and large-ripple wavelength (light purple 16 and dark green 17 circles) with increasing atmospheric density. Ripple data follows a previously derived scaling relationships for drag ripples 18 ; best-fit power law to an updated global compilation of large-ripple wavelength has an exponent of -0.63 (green dashed line) 17 . The width of the hydrodynamic anomaly (gray shade) was calculated for a 125-μm grain size 19 and assuming that wind shear velocity is at the transport threshold 20 (Methods). The statistical significance of the regression slope was computed using a Wald Test, and margins (smaller than line width( indicate a 95% confidence interval.

Sensitivity of the derived bedform gap to transport threshold model
To test whether the width of the bedform gap we calculate for Mars is affected by the chosen transport threshold model, we tested three different formulations derived from theoretical and empirical considerations. We find that the bedform gap is bounded by the empirically derived scaling laws predicting ripple wavelength and dune length as a function of atmospheric density for all tested models; thus, our conclusions do not depend on the choice of the threshold model ( Figure S7). We also calculated the width of the bedform gap equating wind shear velocity to several impact threshold models [21][22][23] . We find that the predicted gaps overlap with some observed dune data for models that predict stronger transport hysteresis ( " * " " * # <~0.3., consistent with the fact that true formative shear velocities ought to be larger than the impact threshold for sand transport. We note that the model proposed by Andreotti et al. (2021) 20 fits dunes and ripples dimensions slightly better than the other two models we employed; our measurements could potentially be employed to calibrate threshold models.

Seasonal variations in atmospheric density
To test if the trends we detect are sensitive to seasonal variations in density, we compared results obtained using the hypsometric equation (Equation M3) with outputs from a global climate model (GCM) 26 . Because we are only interested in the value of atmospheric density at times of active sand transport, we calculate the average annual atmospheric density, weighted by the magnitude of sand flux as calculated from GCM outputs, as a characteristic value of atmospheric density during the dune-forming season. We use the GCM to extract atmospheric density and wind stress globally on Mars in intervals of 10° solar longitude, and estimate the sand flux using the formulation of 27 , where # = 6.1, * & is the threshold wind shear velocity (estimates using the value used in 28 as is appropriate for GCM modeling), * is the surface wind shear velocity, ' is the density and = 3.71 m s (! is the acceleration due to gravity. We verified that our result does not significantly change when other flux laws are employed 29 . The observed trend between dune size and atmospheric density does not significantly change when density is computed using the GCM ( Figure S8). The GCM predicts a slightly higher atmospheric density in latitudes > 70°N relative to when assuming an isothermal atmosphere (Figures 1 and S6-S7), which results in a slightly lower fitted power-law exponent. When only fitting data in latitudes < 70°N, we again obtain an exponent that is close to the theoretically predicted -2/3. We note that scatter in dune size is greater near the pole, implying other mechanisms may be involved in controlling dune size there, as discussed in the next section ( Figure S9). Furthermore, GCM predictions are likely to be less accurate near the poles, where significant filtering is applied due to the nature of traditional latitude/longitude GCM grids with converging meridians; as a result, the GCM may fail to capture local phenomena such as slope or thermal-contrast winds near the polar cap.
Supplementary Figure S8: Decrease in dune length with mean atmospheric density weighted by sand flux magnitude, computed using a global climate model. The slightly higher mean density in the polar region of Mars decreases the exponent of the fitted power law. However, when removing data > 70°N, we obtain a similar power-law exponent as in Figure 3 (-0.58). The resolution of the data is affected by the GCM resolution near the pole.

Potential impact of near-surface volatiles at high latitudes
Overall, Mars' smallest barchan dunes (2 nd percentile, Table 1) have widths and lengths of ~80 m, i.e., they are ~5-6 times larger than their terrestrial counterpart. Dune size and shape vary between equatorial and polar regions, which captures an 11 km latitudinal average elevation range and a three-fold increase in density. In the north polar region (>70°N), which has the greatest abundance of dunes at the lowest elevations, barchan dunes are on average smaller than barchan dunes at lower latitudes (70°S-70°N; Table 1) and are more elongated (see discussion in 30 ), with the longest and narrowest dunes concentrated along the sand sea's northernmost margin ( Figure S9). Compared with lower-latitude dunes, which are mainly located within impact craters ( Figure S10), north-polar barchans have a more contiguous distribution, providing an opportunity to investigate spatial patterns in sizes and shapes. Much like length-to-width ratios, dune volumes vary dramatically across the north pole ( Figure S9). In Olympia Undae, barchans are 2-3 orders of magnitude more voluminous relative to the global average.
Lower volumes and high length-to-width ratios identified near the north polar ice cap ( Figure  S9) might hint at a link between dune morphodynamics and volatile distribution. Areas with maximum modeled 26 annual surface CO2 frost >40 cm (e.g., ±30°E in Figure S9a) tend to host fewer dunes, supporting the hypothesis that volatiles exert a significant control on sand mobility. By limiting sand availability, surface frost could hinder dune growth, accounting for locally less voluminous dunes 31 . Furthermore, more efficient rebound of saltating grains on a frosted bed likely increases the saltation saturation length (the spatial scale over which sand flux equilibrates to its saturated value), which in turn could result in more elongated barchan morphologies 32,33 . However, a combination of other factors could play a role in creating the observed peculiarities of barchans along the northern edge of the north polar sand sea, including complex interactions between dune fields, the atmospheric boundary layer, and polar-cap winds [34][35][36] .

Comparison of dunes within and outside of impact craters
The effect of crater topography on the wind field and resulting eolian landforms was investigated locally, e.g., in Gale and Herschel craters 37,38 , and differences in dune morphology within and outside craters have been previously noted 39 . These studies found that the concave shape of craters affects dune morphometrics in three main ways -by deflecting regional winds, by inducing local katabatic winds 40 , and by trapping sediments 41 , locally increasing sediment supply. Our analysis provides quantitative evidence consistent with these effects on a global scale. We find, as predicted by models (e.g. 30 ) that intracrater dunes have greater volumes and are more elongated (i.e., they have higher L/W; Figure S10). The median volume of an intracrater barchan is about 2.28×10 5 m 3 , almost twice that of dunes outside of craters, ~1.45×10 5 m 3 (p≪0.01, 1-tailed t-test).