Ice needles weave patterns of stones in freezing landscapes.

Patterned ground, defined by the segregation of stones in soil according to size, is one of the most strikingly self-organized characteristics of polar and high-alpine landscapes. The presence of such patterns on Mars has been proposed as evidence for the past presence of surface liquid water. Despite their ubiquity, the dearth of quantitative field data on the patterns and their slow dynamics have hindered fundamental understanding of the pattern formation mechanisms. Here, we use laboratory experiments to show that stone transport is strongly dependent on local stone concentration and the height of ice needles, leading effectively to pattern formation driven by needle ice activity. Through numerical simulations, theory, and experiments, we show that the nonlinear amplification of long wavelength instabilities leads to self-similar dynamics that resemble phase separation patterns in binary alloys, characterized by scaling laws and spatial structure formation. Our results illustrate insights to be gained into patterns in landscapes by viewing the pattern formation through the lens of phase separation. Moreover, they may help interpret spatial structures that arise on diverse planetary landscapes, including ground patterns recently examined using the rover Curiosity on Mars.

Patterned ground is widespread in cold regions (1–8). Distinct sorted circles (1, 2, 5), stripes (1, 3), and labyrinths (5, 9), ranging in scale from a few centimeters to several meters, form during repeated freeze–thaw cycles as stones in the soil segregate according to size due to complex interactions among water and heat transport processes (10–13). However, stone patterns can evolve over thousands of years or more in the field (14, 15), making it difficult to determine their spatiotemporal behavior (see Fig. 1 and for more details). As a result, the fundamental physical principle underlying the phenomena and the quantitative characterization of the self-organization process remain unclear. Aside from creating stunning nonbiologic patterns in the natural world, the phenomenon merits close attention as it manifests soil processes important in carbon cycling in permafrost regions (5, 16, 17), which both affect and are affected by the rapid current warming in the Arctic, and it constitutes an extreme natural example of self-organization in excited granular materials (18–20).

Fig. 1.

Sorted patterned ground on Earth and Mars. Note range of length scales. ( (ESP_030222_1220) image of clusters of boulders southeast of the giant Hellas impact basin, Mars. Image credit: NASA/JPL/UArizona. (B) Labyrinths of stones in Svalbard, Norway. Image credit: B.H. (C) Polygons in the Swiss Alps. Image credit: N.M. (D) Circles in Svalbard, Norway. Image credit: B.H. (E) Stripes in Hawaii, USA. Image credit: B.H.

Phase separation has in the past decades become a central physical principle for self-organized patterning in cell structure (21, 22), gravitational fluid (23, 24), active matter (18, 25–28), and ecological systems (29–31). Herein, we suggest that the principle leads to insights into pattern formation in geomorphic systems. It is important to note that it is distinct from the principle at the core of sophisticated phase-field models, which arise from regularized partial differential equations designed to solve moving boundary problems (32, 33), such as solidification fronts (34, 35), fracture (36, 37), and so on. The concentration-dependent movement feedback that leads to the spatial aggregation and dispersion of particles resulting in diverse patterns is central to the phase separation principle underlying our models. It contrasts with classical models with scale-dependent feedbacks; these include models using the foundational and widely applied Turing principle (see refs. 38 and 39 for details) and diverse other models, including those of systems in which organisms attach themselves providing positive feedback to aggregate at a small scale while creating negative feedback at a large scale, and of spatial patterns that form in near solidification fronts (35, 40–42). When a binary mixture, such as stones and soils, is subjected to periodic external forcing, for example, freeze–thaw cycles or cyclic vibrations, particles autogenically separate and commonly form distinct spatial patterns. This self-organization is similar to the well-known phase separation of a mixed fluid into two phases, in which separation driven by minimum free energy results in distinct spatial regions (43–45). Here, we investigate patterned ground self-organization by examining the dynamics of stone movements in laboratory investigations. We show that phase separation is expected during repeated freeze–thaw cycles, and using two models, we demonstrate how the concentration-dependent movement of stones produces diverse spatial patterns.

Results

Laboratory Experiment.

Using a laboratory experiment with wet soil subjected to freeze–thaw cycles (46, 47), we demonstrate small-scale (1 to 10 cm) pattern formation resulting from cyclic growth and decay of needle ice on the ground surface developing either on level or sloping surfaces of a 100 by 50 cm soil area (see Materials and Methods, and ). Our experiments addressed two key factors controlling the formation of stone patterns: the fraction of the ground surface covered by stones (hereafter termed stone concentration) and the slope of the surface (). Stone patterns developed from stones that were initially evenly distributed. Within the first five freeze–thaw cycles, small stone clusters form and progressively merge with adjacent ones (Movies S1–S3). Coarsening of stone clusters occurs (48, 49), and well-defined patterns are eventually generated (see Fig. 2 and for more information). On level ground, the patterns range from dispersed stone clusters to connected irregular labyrinth forms as the stone cover increases from 20 to 60%; isolated stone-free islands form where the initial stone concentration is high (reference Movies S1–S3). On inclined ground, as the slope increases from 5° to 7°, the stone patterns transition to stripe-like patterns aligned downslope; clear patterns did not form on steeper slopes (46).

Fig. 2.

Self-organization of stone patterns and stone motion in the laboratory. (A) Spatial pattern development starting from a uniform 30% stone cover through 30 freeze–thaw cycles. The panels cover an area ∼0.4 m on a side. (B) Box violin plots of the speed of individual stones for three configurations, homogeneous state at 80% cover (no patterns), isolated stones, and arrested stones (within patterns) at 40% cover (patterned), respectively. The boxplot spans the 25th to 75th percentiles. Red dots indicate mean values. The statistical analysis was implemented with Welch’s one-way ANOVA, , , and the significant levels were adjusted with the Benjamini–Hochberg method for comparisons among groups. (C) Relationship between stone speed and local stone concentration (within ∼3-cm radius, see F) during ∼30 freeze–thaw cycles. The solid line is the function with a decay rate of , , and (reference for models’ selection and statistical parameters). (D) Similar to C, the residual speed but removal of the effect of concentration beyond twofold radius. (E) Inverse relationship between the needle ice height and number of stones in clusters ( per treatment). The significance of changes among treatment clusters: *, ,and “NS” for . Error bars represent one SD in C–E. (F) Schematic diagram shows the definition of the local cover and concentration for the center stone in the analyses of C and D. Reference Movies S1–S4 for sorted patterned ground and stone movement.

Physical Mechanisms Underlying Patterned Formation.

The stone patterns in our experimental system result from the interactions between stone concentration and the amount of needle ice that forms under the stones. Needle ice growth and subsequent collapse due to thawing are the direct drivers of stone motion (). The horizontal transport of individual stones is stochastic and dependent on the length and curvature of the underlying ice needles and local stone concentration. Contrary to common assumption, much of this transport occurs during the freezing phase (Movies S4 and S9). The amount of needle ice that forms and the resulting stone motion in areas with low-stone concentration (sparse stone areas) exceed those in areas of high-stone concentration (dense stone areas) (). This drives net lateral transport of stones toward areas with high-stone concentration and further increases the stone concentration there ( and Movie S4). This positive feedback is further enhanced by the recently recognized tendency for ice needles to curve and to move stones they lift toward stone-rich areas (46, 47). On the other hand, high-stone concentrations reduce needle ice growth, needle ice curvature, and stone motion. This effect results in the stacking of stones in a direction perpendicular to the soil surface (Movie S4), rather than moving and dispersing them laterally, and hence creates a negative feedback to stone aggregation. On sloping ground, in addition to lateral sorting processes, downslope frost creep and toppling or sliding failures during thaw phases play important roles in net downslope displacement of stones, leading to elongated patterns oriented downslope (46). Together, these physical processes that drive stone movement and pattern formation are critically dependent on the spatial variation of stone concentration, as well as the cumulative stone displacements caused by needle ice growth and collapse.

The fact that the stochastic transport of stones depends on the local stone concentration and surface slope implies that their dynamics are nonlinear, potentially leading to pattern forming instabilities. The transport of stones is strictly mass conserving, and pattern formation can arise either through nonlinearities arising from thermodynamics, as in phase separation in binary alloys (50), or from motility-induced phase separation (27, 28, 51). This latter possibility is clearly supported by our laboratory experiments that document how stones move as if self-propelled through their interaction with ice needles and the local stone concentration field around individual stones. In general, isolated stones move significantly faster than stones with close neighbors [Fig. 2; , ]. Our data analysis yields a strong negative relationship between stone speed and local stone concentration at multiple scales (Fig. 2 and ); stones slow down when the surrounding stone concentration increases. This negative relationship is robust in terms of residual analysis, even when excluding the larger-scale interactions ().

Phase Separation Implications for Patterned Ground Formation.

We now consider two classes of models that quantify the relation between the stone speed v(S(), H()), local stone concentration S(), and height of ice needles H(), where represents position. As a first approximation, we assume that these effects can be separated so that . This parametric velocity dependence on space implies that the dynamics of the stones will show important deviations from Fick’s Law, as is known to occur in other systems, such as bacterial motion, where the effective mobility is concentration dependent (52). We represent the dependence of stone motion on concentration (Fig. 2) as [this constitutes the core of Model 1 with ; see Fig. 3 and for more information]. Experimental data also revealed a relationship between local stone concentration and H, the height of ice needles. H declines with increasing number of stones in clusters (groups of overlapped and constrained stones) from 1 up to 64 (Fig. 2), showing that stone speed decreases with increasing concentration S and decreasing H. Thus, we assume that is an increasing function of H and explore consequences of specific assumed forms for in Model 2 (Fig. 3 , Top).

Fig. 3.

Phase separation and the diversity of self-organized stone patterns. (A and B) Phase diagrams illustrate the dependence of pattern formation for both models (Top) on stone concentration and speed decay rate () and define the regions in which the phase separation switches from the one- to two-phase regimes (Bottom): stone-poor and -rich phases. Each point () represents a parameter set (, ) for the model equations. The theoretical predictions (white lines) coincide with numerical results. (C) Simulated 2D stone patterns with increasing stone concentrations from 20 to 80% for both models. Numerical simulations were implemented using periodic boundary conditions with parameters , , , and (reference for additional details). Islands transform to labyrinths and polygons with increasing stone concentration, and no pattern develops at high-stone concentration. Color bar represents stone concentration in units of g/cm2. Spatial scale with 5.1 m. (D) Self-organized patterns formed in the laboratory within 240 to 360 h due to repeated freeze–thaw cycling with stones initially laid out evenly on the soil surface: islands (20% cover), labyrinths (30% cover), polygons (40% cover), homogeneous state (80% cover), and incipient stripes (20% cover on a sloping surface of 7°; white arrows show the general downslope direction) [reprinted with permission from ref. 46]. Reference Movies S5 and S6 for numerical simulations of Model 1 and 2, respectively.

For both models, the change of spatial stone concentration with time obeys the law of conservation of mass, , with a local stone concentration that is composed of gradients of S and H, along with rotationally invariant scalars of which the lowest order is . This results in the following:Here, , where is the collision rate per stone that depends upon position and direction of motion (52), is the spatial dimension (i.e., ref. 2), is the coefficient of potential energy forces (equivalent to the dispersal coefficient by a unit stone at a nonlocal scale) (52), and is the characteristic width of the mobile convergent front of the low- and high-concentration phases. In principle, there could be a higher-order term in H, but this is neglected here as the characteristic front width is likely to be much smaller than for S. The mobility of individual stones can be expressed by an effective diffusivity , where in the models. To complete our description of the coupled needle ice stone system, we represent the negative feedback between the stone concentration and needle ice height by the following equation (Fig. 3 ):On the right-hand side, the first term describes the water input in the soil. The second term describes the inhibitory effect of increasing stone concentration on needle ice growth; a and m are positive empirical coefficients, where m represents the sensitivity of ice growth to stone concentration ( is assumed here, because a linear coupling between S and H is ruled out by the upward curvature visible in Fig. 2 and ). is the specific rate of loss of the H2O (both water and ice) due to evaporation, and is the diffusion coefficient representing the H2O (both water and ice) transport processes during a freeze–thaw cycle. Eqs. and constitute a complete description of the system dynamics. The form of Eq. is reminiscent of the equations governing phase separation in alloys and is expected to lead to pattern formation. To check this, we compared our experimental results with the two models of patterning resulting from ice-induced stone displacements: one dependent on stone concentration only (Model 1, Fig. 3 and Movie S5), and the other also including explicitly the dependence on the height of ice needles (Model 2, Fig. 3 and Movie S6). Our models yield diverse spatial patterns, similar to those emerging under a range of experimental conditions (Fig. 3 ) and at various field sites (Fig. 1 ). They include well-defined stone islands, as well as stone circles, and labyrinths on level ground. There are two end members in our theoretical phase separation model: bare soil with no stones (S → 0) and soil with the ground surface fully covered by stones (S → 100%). According to theory, these single-phase regions are bound by curves, , so-called spinodal lines as shown in Fig. 3. In the two-phase regions, both islands and labyrinthine patterns emerge in a finite concentration range. This theoretical phase separation boundary coincides with laboratory experiments, where regions of intermediate stone concentration separate distinctly into homogenous regions with either sparse or dense stones (Fig. 3 ). These results illustrate the important role of needle ice height in shaping the sorted patterns with different stone concentrations (Fig. 3 and ). In the experiments, the self-organization was pronounced, forming distinct stone patterns, including stone islands, labyrinths, or circles, which emerged after only 30 freeze–thaw cycles in experiments with initial stone concentration below S (ca. 60%). Above this critical value, only a few large clusters of stones formed (Fig. 3 and Movie S7). Note that the phase separation models also reproduced stripe-like patterns on sloping ground by incorporating downslope soil creep effects into the models (Fig. 3 ).

As is well-known, phase separation initially arises from a linear instability in which particles diffuse from low concentration to high concentration, that is, opposite to Fick’s Law, until nonlinearities overwhelm the aggregation process, leading to pattern restabilization through scale-invariant dynamics (53). To test these predictions and the expected relationship between stone velocity and concentration in phase separation models, we conducted an experiment with three domains with different initial stone concentrations as shown in Fig. 4 (see Materials and Methods and ). The velocity fields in both the model and experimental data coincide with one another: velocities shifting from radially outward to inward with increasing distance from the center (Fig. 4 ). The net radial flux is outward for lower concentration, but it is inward for high concentration. It is noteworthy that two clear dips appear at the interfaces between stone-poor and -rich domains. Interfacial phenomena are strongly implicated in coarsening dynamics for phase separation described by the quartic gradient term of Eq. . Furthermore, as the stone concentration increases from 10 to 80%, needle ice crystals decrease in both height and curvature (Fig. 4 , respectively), confirming the negative feedbacks between stone concentration and motion and providing further impetus for extending the phase separation model with a single variable (Model 1) to one (Model 2) with two variables. Overall, our results indicate that the phase separation models capture the key elements of self-organization in wet soils undergoing diurnal freeze and thaw cycles.

Fig. 4.

Comparing experimental and theoretical results for well-defined, triple concentration patterns. (A) Vertical view of initial arrangement of stones in concentric rings of increasing stone concentration away from the center. (B) Observed net radial flux of stones, positive for outward and negative for inward. (C) Stone displacement field from phase separation Model 1 with . (D) Stones diverge from the stone-poor central domain and move inward from the inner side of stone-rich peripheral domain. (E) Radial decrease in needle ice height with increasing stone concentration during first freezing period. Inset shows the area photograph. (F) Box plots showing needle ice height in the three concentric 10, 40, and 80% stone domains shown in A. Boxes extend from the lower to upper quartile values of the data. Horizontal lines mark the median heights. Comparisons are as follows: * and **. Reference Movies S8 and S9 for experimental details.

Scale Invariance and Self-Similar Dynamics of Patterns.

To determine whether our data displayed nonlinear restabilization and dynamic self-similarity expected from the models, we analyzed all images in each experiment (e.g., ∼2,160 images for 30 freeze–thaw cycles in an experiment with 40% stone cover) and characterized the length scale of the observed patterns using two-dimensional (2D) Fourier analysis (54) (power spectra were derived using a square, moving window method). We calculated the temporal evolution of the dominant wavelengths for level and sloping ground surfaces (). The wavelengths in the early phase of the experiments are in quantitative agreement with the simulations. They increase considerably and nonlinearly in the first 180 h within 15 freeze–thaw cycles; subsequently, the patterns stabilize (Fig. 5 ). Importantly, the models qualitatively account for this early coarsening behavior, which is robust and independent of initial stone concentration and surface inclination. For sloping ground, experimental results coincide with the t scaling law that is standard in phase separation models (22, 29) during the first 180 h (∼15 freeze–thaw cycles; Fig. 5 and Table 1). In contrast, for level ground, experiments reveal a t scaling law (Fig. 5 and Table 1). This type of scaling law agrees with the results of a quite different grain-scale model developed by Kessler et al. (2, 10) who focused on seasonal frost heave and subsurface ice growth as formative mechanisms for the larger scale (∼2- to 4-m diameter) patterned ground.

Fig. 5.

Scaling behavior of the pattern coarsening from experiments and simulations. (A and B) Temporal evolution of dominant pattern wavelength on sloping and level ground. Note, the datasets are offset from one another for graphical clarity. Experimental data (colored solid lines with symbols) and numerical simulations with both phase separation models (solid lines). The dashed lines fit the experimental data with a power law at early stages. (C and D) Scaling behavior of the self-organized patterned ground from experiments and models. The rescaled structure function , as a function of wavenumber q/q, , versus the number of freeze–thaw cycles for experiments and simulations (see Materials and Methods for details). (Insets) The scaling of dominant wavenumbers plotted against the freeze–thaw cycles.

Table 1.

Statistical properties of coarsening dynamics on experimental evolution at early stages

Experiment	Exponent*	SD	P value	R 2	n
Sloping ground	0.34	±0.010	<0.0001	0.78	3
Level ground	0.23	±0.012	<0.0001	0.87	4

Note that the nonlinear regression model was used to obtain the exponent with “fitnlm” function in MATLAB 2020b. n indicates the independent experimental replicates.

It is noteworthy that the relatively fast coarsening behavior on sloping ground can be attributed to the more directional movement with a bias downslope than that random stone movement in level scenarios (), as well as the progressive elimination of complementary stripe defects, so-called stripe termination pairs (3). A similar phenomenon occurs in other systems in which there are long-range interactions, such as in block copolymers that only exhibit microphase separation at long times. However, at short times, even these systems exhibit phase separations, although the exponent describing the growing length scale is frequently found to be of order 1/4, at least when the interface width is not very small with respect to the domain size (45). Here, the oriented movement can result in overlapping stones on sloping ground, whereas level ground display quasiconservation (). The deviation from the scaling laws after 15 freeze–thaw cycles presumably occurs because the models do not consider other geomorphic processes (55). This slow down and 0.25 power-law behaviors could in principle be captured with an additional term that describes interfacial dynamics, such as the celebrated Kardar–Parisi–Zhang dynamics (32, 45, 56, 57).

To quantify the emergent length scales (58) more precisely, we calculated the dynamic structure factor , where is the modulus of the wave vector, that is, the Fourier variable conjugate to . At large , the functional form of is determined by the boundaries between phases and in 2D is expected to lead to a variant of Porod’s Law, where q-3 (59). Both model and experiment follow this universal trend. Furthermore, the rescaled structure factors collapse to form a master line when is plotted versus , where is the time-dependent wavenumber corresponding to the width of spatial patterns from , indicating that numerical simulations and experiments are undergoing the same coarsening process (Fig. 5 ). The reason that does not vanish as q tends to zero is that the stones overlap so that effectively in 2D the global conservation of mass does not hold. The time evolution of is shown in the insets of Fig. 5 . Whereas the scaling collapse reported above is expected in any phase separation process with a single, emergent length scale, the time dependence of the length scale reflects finite-time and finite-size effects and need not be power law. In theoretical studies of motility-induced phase separation, the exponents have also been found to be in the range 0.2 to 0.3, possibly indicative of slow crossover to a value 1/3. The observation of these dynamic scaling laws in our experiments and models is an example of motility-induced, sorted patterned ground formation going beyond phase separation in the active matter (24, 28).

Discussion

Two models currently exist for the origin of patterned ground; both are based on grain-scale numerical simulations. In the first, Werner and Hallet (3) suggested that differential growth of needle ice (between soil domains having different elevations and different abundances of stones) produced stone stripes with a typical spacing of ∼0.1 m. In the second, Kessler and Werner (2) proposed that a diversity of forms of meter-scale patterned ground emerge spontaneously from two feedback mechanisms associated with subsurface ice growth (lateral sorting and stone domain squeezing during ground ice freeze–thaw). These seminal models have been very instructive, but they have yet to be thoroughly tested or validated. Here, we took quite a different approach. Starting with unprecedented results from laboratory experiments, we developed a theoretical foundation, building on recent phase separation studies but with a modification that arises from the active or nonequilibrium nature of the needle ice–driving force. Our phase separation model provides a coherent explanation of the early stages of sorted patterned ground but does not aim to be a numerically predictive model of its long-term evolution (2, 3). By modeling the feedback between the amount of needle ice and the stone concentration, we reproduced the field-observed diversity of stone patterns, including stone circles, labyrinthine patterns, and islands on level ground. Incipient stripe-like patterns developed on sloping ground, but the formation of very well-defined regular stripes in laboratory may require a much larger experimental apparatus and larger numbers of freeze–thaw cycles.

Theoretically, we remark that the pattern formation mechanisms discussed here involve large-scale instabilities throughout the region of interest and are of a somewhat different character than other geophysical pattern formation phenomena that arise from front propagation. For these phenomena, such as crack formation in basalt columns (36, 37), advanced moving boundary methods can be used, for example, using so-called phase-field models (33, 34).

The pattern forming mechanisms driving the phase separation process in Arctic sorted patterns results from cyclic freezing and thawing. Because of that, changes in rock pattern formations may be indicative of changes in freezing conditions, pointing to permanent thawing in areas where patterns are observed to be lost. Our results support these predictable arguments that sorted patterned ground are not only by themselves sole local movement processes (Fig. 3) but also dependent on the feedback effect of needle ice on movement speed (Fig. 3 and ). Hence, changes in the characteristics of rock patterns may be a crucial first sign of local climate change impacts, pointing at potential future soil loss. This may be observed on Arctic regions on earth, but changes in rock pattern formations may equally be indicative of changing conditions on Mars or other planets and could in the future be used to study temperature changes as well as the dynamics of soil conditions with planetary studies.

Aside from its intriguing geometric regularity, patterned ground is important because it may provide valuable information about surface processes and conditions in remote or hostile regions where detailed observations or monitoring are difficult or impossible, both on Earth and beyond. For instance, changes in patterned ground may signal subsurface changes in the vast permafrost regions of the warming Arctic where instrumentation is extremely sparse at best (16). Visible changes in patterned ground could provide important clues about the release of greenhouse gases from the permafrost to the atmosphere (60). Novel types of patterned ground on Mars revealed in detail by the Curiosity rover are currently being studied to further understanding of the underlying processes and the clues they contain about energy and mass exchange between the atmosphere and lithosphere on our neighboring planet (61). Our mathematical description of the universal principles that govern pattern formation contributes to understanding the processes that shape ground surfaces found in cold regions and may help in identifying the impacts of global climate change on our own planet.

Materials and Methods

Sorted Patterned Ground Implementation in Laboratory.

The laboratory experiments were implemented with the same procedures and environmental conditions but on level ground and sloping ground, respectively. We first conducted a series of experiments in the controlled cold room to simulate patterned ground formation on level ground conditions with 20, 30, 40, 60, and 80% stone concentration, being ∼648, 912, 1,226, 1,916, and 2,597 stone individuals, respectively. The second series of experiments targets patterns on sloping ground by changing gradient from 5 to 11% (ref. 46) (). The height of needle ice, soil temperature, and moisture were recorded at 10-min intervals with two multichannel data acquisition systems in each experiment (46, 47). Once each experiment had finished, all images of the surface patterns from top and side view were connected to produce a video (Movies S1–S3 and S7).

Feedback Between Stone Clusters and the Height of Ice Needles.

In addition to concentration-dependent movement determining the dynamics of self-organized patterns, the height of ice needles contributes a crucial role to stone movement and pattern evolution. To quantify the relationship between local stone concentration and height of ice needles, we manually measured the height of ice needles below stone clusters of 1, 2, 4, 8, 16, 32, and 64 individual stones of 4 to 8 mm in diameter, which were stacked into seven circles naturally corresponding to diameters of 0.5, 0.8, 1.5, 2.0, 3.0, 4.0, and 6.0 cm, respectively. In general, we found that the needle ice height is inversely proportional to the number stones of the clusters from 1 up to 64 and is consistent with the negative relationship between stone movement and the height of ice needles. We see the same general relationship between stone concentration and stone speed.

Real-Time Trajectory and Local, Concentration-Dependent Movement Analysis.

To quantify the relationship between gradient in stone concentration and pattern formation, and how interactions between individual stones (clusters) determine pattern formation, we manually traced the movement of individual stones in successive images of our experiments and analyzed their characteristics in relation to gradient in concentration. All the movement trajectories of individual stones were recorded manually using the free software Fiji (https://imagej.net/Fiji/Cite, developed by the NIH) with the Track-Mate package (62).

We determined both the speed of and extent of stone concentration surrounding individual stones by extracting the specific stones as they moved from an isolated location toward the patterned inner and open areas (soil). The concentration of stones in the neighborhood of a tracked stone was estimated by measuring the fractional cover of the stones within the distances of 1.2-, 2.0-, 3.0-, 4.0-, 5.0-, 6.0-, 7.0-, 8.0-, and 9.0-fold of the diameter (see Fig. 2 for an example). These specific setups correspond to the spatial scales of about 12, 18, 24, 30, 42, 48, and 54 mm for stones with an ∼6-mm stone diameter. Following the methods proposed by van de Koppel and coauthors (63), images were converted to binary bitmaps indicating the presence or absence of stones using a custom-made MATLAB program. All circles were extracted from these bitmaps, with the tracked stones set at the center. The central circle of the onefold radius was excluded as it contains the tracked stone itself. To covert the cover estimates to concentration (as it is convenient to use continuum variables in the theoretical model), we precisely measured the local concentration of stones with 100% cover and their weight. Finally, we obtained the conversion factor of 4.0 in our laboratory experiments; that is, local stone concentration is about 4.0 g/cm2 when the stone cover is 100%.

Statistical Analyses.

A one-way ANOVA was used in R (64) to test whether speeds differ among stones in homogenous, isolated, and patterned scenarios. A Bonferroni correction was employed, and all P values below 0.05 were considered significant. Statistical details can be found in the main text and figures. A data point was considered an outlier if it was greater than the 75th quartile + (1.5 × interquartile range [IQR]) or lower than the 25th quartile − (1.5 × IQR). Furthermore, we analyzed the relation between local stone concentration and speed with a generalized linear model using hyperbolic and exponential functions, respectively. The best single-scale model as well as the two-scale multiple model were selected from all possible sets using Akaike’s information criterion. The summary results of the movement relationships are listed in . The two-scale multiple model shows that the negative feedback always occurs at the largest radius based on the Z-score test. This implies that a two-scale feedback relation is unrealistic for our self-organized patterned ground experiments. This was further confirmed by the correlation coefficient of changed spatial scales as shown in , in which the experimental data revealed positive feedback even at the scale of ninefold the diameter.

Velocity Field and Radial Flux Analyses.

To quantify feedback between stone concentration and displacement field, a triple stone concentration sorted annulus experiment had been designed. Inside diameters of outer, intermediate, and inner circles are 40, 30, and 10 cm, respectively. Out annulus is full with 80% stones, while intermediate annulus and inner circles are full with 40 and 10% stones. All the stones were placed on the surface of soil in the container and subjected to 20 freeze–thaw cycles with air temperature oscillating between −5 and 10 °C for 12 h in a cold room. We obtained the velocity fields with Particle Image by Velocimetry PIVLAB package (version 2.02) in MATLAB 2019a (65). Contrast limited adaptive histogram equalization was used to enhance contrast, and a high-pass method was used to filtrate out the low-frequency signal during preprocessing of the images. For the laboratory experimental images and simulated data, we defined the point with the minimum velocity within the circular patch as the center of radial direction. The displacement velocity was obtained with two consecutive snapshots. Furthermore, we calculated the amount of net radial flux at distance with formula .