Salt polygons and porous media convection

From fairy circles to patterned ground and columnar joints, natural patterns spontaneously appear in many complex geophysical settings. Here, we investigate the origins of polygonally patterned crusts of salt playa and salt pans. These beautifully regular features, approximately a meter in diameter, are found worldwide and are fundamentally important to the transport of salt and dust in arid regions. We show that they are consistent with the surface expression of buoyancy-driven convection in the porous soil beneath a salt crust. By combining quantitative results from direct field observations, analogue experiments and numerical simulations, we further determine the conditions under which salt polygons should form, as well as how their characteristic size emerges.

A challenge to any explanation is that the polygonal patterns in salt crusts are remarkably similar wherever they occur, despite local differences in geology, salt chemistry, and environmental conditions. For example, while crusts can vary from sub-centimeter to meters in thickness [2,12,18] they predominantly express polygons about 1-2 m across [2,11,12,[15][16][17][18][19][20]. The same patterns also appear for pure halite crusts such as those of Badwater Basin [11] and pans in Iran [12] as well as the sulfate-rich crusts found for example at Owens Lake (California) [11], Sua Pan [21], the coastal sabkhas of Abu Dhabi [14] and Salar de Uyuni [9].
The salt crusts of dry lakes are known to be dynamic over months to years [2,15,17,18,22,23], and couple to other environmental processes. They form in areas with a negative water balance (i.e. where evaporation significantly exceeds precipitation) and the dominant water balance is between a perennial groundwater flow seeping in from distant sources, and the evaporation of this water from the surface of the dry lake [2,4,24]. The groundwater brings dissolved salts with it, which accumulate at the surface in a crust of evaporites. Wind blowing over the crust entrains dust, which forms a significant proportion of the global atmospheric dust production [3,6,7] and of mineral transport to the oceans [25]. As one example, dust from the dry Owens Lake has posed health problems for people living nearby [26,27], and the site is the center of a decades-long, intense remediation effort [28]; a preferred method is the use of shallow flooding to encourage formation of a continuous, active surface crust [29,30].
While playas are inherently associated with dust sources, these sources are temporally and spatially variable [24,[31][32][33], and there is a lack of understanding of the relevant processes of how crusts form, are sustained, and degrade [7,34]. The height of polygon ridges affects the threshold wind speed for dust entrainment [7,35] and complete coverage of a crust reduces the sediment available for dust emission compared to a patchy crust. In order to better account for these complexities we need to understand the relationships between surface and subsurface drivers of playa dust sources [7,31].
Salt crusts also modify evaporation and heat flux from the playa surface [5], and hence the water and energy balances of these fragile environments. In particular, experimental and field research consistently show that continuous crusts reduce evaporation rates on playa surfaces by an order of magnitude compared to climate-derived estimates [8,29,[35][36][37]. As a result, the observed groundwater evaporation rates at sites with active and continuous salt crusts tend to lie in the relatively low range of 0.2-0.7 mm/day, including seasonal variations [29,[37][38][39][40][41].
Research on salt pans has typically focused on either the dynamics of their complex subsurface flows [42][43][44][45][46] or their crusts [15,17,18,35], without considering how these features interact. Here we show that by treating the surface of a salt playa together with the fluid in the porous media of the soil near its surface, the origins, dynamics, length-scale and shape of the polygonal patterns in salt crusts can be apprehended. As we develop and test such a model of salt crust formation, we will show how it makes explicit predictions of crust growth rates, the conditions necessary to form ridges and continuous crusts, and other inputs into dust or groundwater models of these extreme environments. In particular, by combining theoretical analysis, numerical simulations, experiments and field observations, we will demonstrate how the density-driven convective dynamics of groundwater can lead to variations in the salt flux into the crust, with faster precipitation along the boundaries between convective cells facilitating the growth of a polygonal network of salt ridges (see Fig. 2). Under measured field conditions the size of the expected convective cells matches the scale observed for the surface crust patterns and makes accurate predictions about the crust growth rates. Finally, we will present direct field-based evidence of convective dynamics that are correlated to the surface crust patterns.

II. PHENOMENOLOGY OF SALT POLYGONS
Any model for the emergence of salt crust patterns should be able to convincingly explain their commonly observed features. Polygonal networks are one of the more common emergent patterns in geophysics, seen for columnar jointing [47], ice-wedge polygons [48,49], polygonal terrain [50] some types of mud-cracks [51] as well as convection [52,53] for example. As such, a suite of more detailed or quantitative predictions are important to evaluate any model that displays polygonal features. Based on the known phenomenology of salt playa, the most significant tests for a predictive model include: (a) The driving mechanism should be specific to the geophysical conditions encountered in salt lakes [9][10][11][12][14][15][16] and (b) should spatially modulate salt transport, leading to increased salt precipitation at ridges in the crust [2,12,54]; (c) The response should be fast enough, and stable enough, to account for ridge growth on a scale of months [2,15,18,22]; (d) The model dynamics should generate closed polygonal patterns of relatively narrow ridges spaced a few meters apart [11,12,[15][16][17][18]; and (e) the characteristics and scale of the patterns should be robust to environmental variations, including different types of soil ranging from silt to sand [11], salt compositions ranging from pure halite crusts to those with significant amounts of sulphates and carbonates (e.g. trona) [9,11,12,14,37] and any local differences in climate and variables like the evaporation rate [29,[37][38][39].
Previous models of salt polygons have mainly considered mechanical effects. The earliest argument was made by Christiansen in 1963 [17]. He suggested that the shapes were the result of compressive stress, due to salt precipitation and temperature changes, which caused the crust to first buckle and then crack as the growing ridges thrust upwards. He noted that the spacing of the buckling features should be proportional to the thickness of the stressed crust and also depend on its strength; these predictions are consistent with the modern understanding of wrinkling phenomena [55][56][57]. However, wrinkles due to isotropic compression typically form parallel features, such as stripes, rather than closed shapes [55][56][57]; polygons can form by compression, but only by folds into an elastic layer, rather than ridges projecting out of a layer [57,58]. A more fracture-focused model was proposed by Krinsley [12] in his 1970 survey of Iranian salt pans. He argued that desiccation places the crust into tension, rather than compression, which is relieved by fracture. The cracks then allow for groundwater to more easily reach the surface to evaporate, so that ridges grow up along them. As the crust develops faster at the ridges, these can buckle or thrust upwards (as in Fig. 1(b)). Like wrinkling, for cracks the expected feature size is a few times the thickness of the cracking layer [59,60].
Both the wrinkling [2,54] and fracture [18][19][20]61] models can explain the qualitative appearance of cracks and thrusting at the ridges in salt crusts. The main difficulty with any purely mechanical model, however, is lengthscale selection. The natural spacing between buckles or cracks is proportional to the thickness of the layer under stress [55][56][57][58][59][60]. As some geophysical examples, the scale of ice-wedge polygons in permafrost is determined by the depth at which annual temperature changes diffuse into the ground [48,62], whereas the size of columnar joints in lava records how it cooled, and in particular the thickness of a cooling front invading into an initially molten lava formation [47]. For salt playa, crust thicknesses vary widely, so cannot explain the robust polygonal scale seen in nature, and no other convincing explanation for the feature spacing has been proposed.
As an important mechanism one may also consider the crust composition itself, which can include minerals in different hydration states. For example, thenardite (Na 2 SO 4 ) transforms into mirabilite (Na 2 SO 4 · 10H 2 O), and expands its volume by ∼320%, at temperatures easily seen in salt deserts [63]. While phase changes can undoubtably generate stresses, the difficulty of length-scale selection remains, given similar patterns in crusts of different thickness. In addition, salt polygons are routinely seen in crusts with almost pure halite (NaCl) composition [11,12], where such effects are absent.

III. BUOYANCY-DRIVEN CONVECTION AS A MECHANISM FOR PATTERN FORMATION
As discussed, explanations for the formation of salt polygons have so far focused on mechanical effects. Crust growth is implied in these cases, but the salt source feeding this growth is not modeled explicitly. In other contexts there is a well-developed literature regarding buoyancy-driven convection in porous media, including salt lakes [44,45,64], sabkhas [46,65], carbon sequestration applications [52,66,67] and thermally-driven cases [42,53,68]; the model we develop here builds on a recent study of the linear stability and approach to a statistically steady state of convection in a dry salt lake [64]. However, the effects of such flows on a salt crust have not yet been explored. Here, we will examine the implications of convection in the soil beneath a salt crust and its predictions for the emergence of salt polygons.
We consider a coupling of surface salt patterns to subsurface flows, as visualized in Fig. 2. Briefly, surface evaporation will leave the near-surface groundwater enriched in salt, and heavier than the fluid beneath it. This can lead to precipitation of salt at the surface as well as convection in the soil, with narrow, regularly spaced downwellings of high salinity [64]. We will show that convection should result in a higher salt flux into the surface above the downwellings, where salinity gradients in the groundwater are weaker. We argue that this preferential precipitation of salt templates the growth of ridges at the locations of the downwellings. At upwelling plumes the larger salinity gradients strengthen the diffusive transport away from the surface, rather than advective transport towards it, and crust growth is slower.
Our main field site, Owens Lake ( Fig. 1(b)), has a typical crust pattern in a well-studied and controlled landscape, and can introduce our modeling assumptions. This dry, terminal saline lake has an aquifer extending from the near-surface to over 150 m depth [69]. The groundwater carries dissolved salts [27,37,70], which collect in an evaporite pan of about 300 km 2 [3,69] and which are particularly concentrated in the fluid within a meter or so of the surface [11,37]. Average annual groundwater evaporation rates of E = 0.4 ± 0.1 mm/day have been measured in areas with active salt crusts [37]. The near-surface soil is a silt to fine sand with high porosity of ϕ = 0.7 [37]. A tortuosity-corrected diffusion constant, D = 1.00 ± 0.24 × 10 −9 m 2 s −1 is representative of the mix of dissolved salts present at the site (mainly NaCl and Na 2 SO 4 ). Efforts to control dust emission from the lake involve shallow flooding [30], vegetation [71], gravel cover and encouraging crust growth [29]. As shown in supplemental movie S1, after a shallow-flooding event the soil is saturated with water, which evaporates and leaves behind salts that crystallize into a continuous crust, covered by a network of slowly growing ridges. Given the high quality of the observational record at Owens Lake, and in order to make quantitative predictions relevant to the remediation of this important site, we will use representative values from it throughout this section.

A. Governing equations
The transport of fluid beneath a dry salt lake can be treated as a Darcy flow in a porous medium, with evaporation occurring at the surface and being fed from below by water with some background salinity (see e.g. [42,44,64,72]). Such a system will naturally develop a salinity gradient below the surface, which can become unstable to convective overturning [42,44,64].
As a prototype case we consider a large, flat salt lake with an average surface evaporation rate E, and where the groundwater is recharged from some deep, distant reservoir. For the subsurface flows we model the volumetric flux, or superficial fluid velocity, q, of a fluid of pressure p moving through a porous medium of constant uniform porosity ϕ and permeability κ. By mass conservation E is also the average of the upward component of the fluid flux, q z , at any depth as well as the recharge rate from the reservoir. The fluid has a viscosity µ and carries dissolved salt, whose diffusivity D can be corrected for the presence of different ions as well as tortuosity [73,74].
The dissolved salt contributes to the density of the fluid, ρ, and hence to buoyancy forces. In salt lakes, thermal buoyancy is several orders of magnitude weaker than solutal buoyancy, and can be neglected [44,64]. Using the Boussinesq approximation, ρ = ρ 0 +S∆ρ, where ρ 0 is the density of the reservoir fluid (including its dissolved salts, see Section III C), and ρ 1 = ρ 0 + ∆ρ is the density of a saturated solution. The relative salinity S mediates between these limits: S = 0 corresponds to the background salinity of fluid entering from a distant source (z → −∞), whereas S = 1 represents a salt-saturated brine in contact with the crust (at z = 0).
The governing equations are then a continuity equation for incompressible fluid flow, an advection-diffusion equation for the relative salinity, and Darcy's law: where g is the acceleration due to gravity, andẑ is an upward-pointing unit vector. These mass and momentum balances can describe a variety of porous media flows, including solutal flows in playas and sabkhas [44,45,64,72,75], CO 2 -rich flows in carbon sequestration applications [67,76] and, with an appropriate transformation of variables, flows driven by thermal buoyancy [42,68,77]. A review of these phenomena was recently given in Ref. [78]. As boundary conditions for the salt lake scenario we take: In order to non-dimensionalize the governing equations, one needs to specify a characteristic length scale L and velocity scale V, with a natural time scale following as T = ϕL/V. Taking advantage of how the evaporation rate sets the average vertical fluid flux everywhere, we use V = E. Then, in the absence of any geometric length scale in the problem (e.g. the layer thickness often used for two-sided convection [78]), we set L = ϕD/E as the distance over which advective and diffusive effects will balance, for fluids moving at the characteristic speed E.

Rescaling Eqs. (1)-(3) then leads to
with a rescaled velocity U = q/E, depth Z = zE/ϕD, time τ = tE 2 /ϕ 2 D and pressure P = κ ϕµD (p + ρ 0 gz). This system of equations is governed by a single dimensionless group, the Rayleigh number: Here, the Rayleigh number also gives the ratio of the speed at which a large heavy plume will fall under its own weight, V B = κ∆ρg/µ, to the upward drift due to surface evaporation: Ra = V B /E. Thus, it is frequently used to characterize the vigor of porous medium convection with a through-flow [42-44, 64, 72]. A more detailed derivation of this model, along with discussion of the assumptions behind it, a linear stability analysis of the onset of convection and a two-dimensional numerical investigation into its statistically steady state, was recently published by some of us [64]. Our aim here is to show how solutal convection beneath a salt lake will couple to the growth of patterns in the crust and to make predictions which can be compared to field measurements. As such, after a discussion of how the emergent dynamics of this model will appear in geophysically relevant conditions, we will develop predictions of the salt flux into the crust, which can give rise to crust patterns.

B. Convective dynamics and scale selection
The model above is designed to match the expected dynamics of groundwater in a dry salt lake with at least a thin salt crust developing on top of a fluid-saturated soil. For the one-dimensional case, where salinity depends only on depth, there is a unique stationary solution to Eqs. (4)-(6), namely S = e Z = exp(zE/ϕD), which represents a salt-rich layer of groundwater lying just below the surface and a balance between advective and diffusive transport of salt [42]. In the absence of any horizontal flows other plausible initial conditions, such as an initially uniform lake (i.e. S = 0 everywhere), will relax towards the stationary solution on the characteristic timescale T [44,64]. Similarly, an initially crust-free lake can be described at first by an unsaturated boundary at Z = 0, with no flux into the surface until the salinity there reaches saturation: this scenario will also develop a boundary layer that reaches saturation, and starts to precipitate, over the same timescale [44]. For representative parameter values from Owens Lake, T is about 6 months, consistent with the known timescale of crust development [15], and L = 15 ± 8 cm.
For a wide variety of initial conditions, therefore, a dry salt lake will naturally develop a heavy salt-rich boundary layer of fluid near its surface. Above some critical Rayleigh number Ra c this layer will be unstable to convection [42-44, 64, 72]. Linear stability analysis shows that for constant, uniform evaporation at the surface, the stationary solution has a finite-wavelength instability above Ra c 14.35 that leads to the growth of convective plumes [42,43,64]. At Ra c the critical wavenumber a c 0.76 for conditions at Owens Lake, corresponding to a wavelength of 1.3 m. Intriguingly, this characteristic scale of the convective instability is of the same order of magnitude as the patterns typically seen in salt crusts.
This instability is also robust to different assumptions about the initial and boundary conditions of the lake. For example, an initially uniform lake (where S = 0 at t = 0) has the same Ra c as the stationary state, along with a very similar spectrum of unstable modes [64]. Alternatively, as would be appropriate for a salt lake with a thin ponding of surface brine, the case of a constant pressure boundary condition at Z = 0 is unstable to convection above Ra c 6.95, with a c 0.43.
To predict the long-time dynamics of convection pattern in a dry salt lake, we simulated Eqs. (4)-(6) on rectangular domains (2D) or boxes (3D), where S = U Z = 1 at the top boundary, S = 0 and U Z = 1 at the bottom, and with periodic boundary conditions along the vertical walls. We conducted simulations from Ra = 20 to 1000. As initial conditions, we started from the stationary state for a bounded domain, S = (e Z − e −h )/(1 − e −h ), with small random perturbations applied to the salinity distribution at all grid points. Most simulation domains had a horizontal extent of 12π (for 2D) or 24π (for 3D) and depth of h = 10. The domain depth introduces an additional length scale, and thus dimensionless group, into the numerical problem. To confirm that the results are in a limit where they are essentially independent of h, we repeated selected simulations for domains with h = 20. Our numerical methods are summarized in Appendix A.
A finite-wavelength instability sets in for all simulations and rapidly grows beyond the linear regime. As shown for the 2D case in Fig. 3(a-c), perturbations around the most unstable mode develop into a set of narrow plumes of high salinity, which fall from the upper surface. In 3D the situation is similar, as in Fig. 3(d,e), with the convection near the surface instead organizing into a closed polygonal network of downwelling sheets, which surround patches of upwelling fluid. At lower depths the sheets break up into isolated plumes. Both the growth rate of the linear instability and the speed of the mature downwelling features scale with = (Ra−Ra c )/Ra c , which can be understood as Ra ∼ V B . As such, we make use of a rescaled time, τ = τ , in order to compare features at a similar stage of development.
As they grow, the initial downwelling features strip the boundary layer of dissolved salt, leaving it depleted of the buoyancy forces that initiated the convection. The plumes also tend to merge together into larger structures. This coarsening process causes the pattern length-scale to increase over time; for further details, see Ref. [64]. However, as the distance between downwelling areas grows, they become less efficient at draining salinity from the boundary layer, which can start to thicken again. This allows for the growth of additional instabilities. In 2D these take the form of small plumes that appear between, and are attracted to, the larger established downwellings. In 3D, although new point-like plumes do also occur, the initiation of linear features is more typical. As shown in Fig. 3(d,e), these are weaker downwelling sheets that connect two sides of polygonal convection cells. As they grow, they also move towards the boundaries of the more established, larger convection cells, eventually merging with an existing edge. These transient 'microplume' features become more prevalent at higher Ra, overlaid on the more stable 'megaplumes' of convection cells [53], which extend much deeper than the boundary layer. Over long enough times a statistically steady state develops [64] where the rates of plume formation and merging balance, and a well-defined length-scale emerges.
These dynamics can be compared with analogous results from two-sided convection without throughflow [52,53,68,79], which has become an important model of carbon geosequestration. For this, some caution is required, as the characteristic scales and Rayleigh number are defined differently in the two cases. Comparing the two systems, the patterns near the walls are remarkably similar, with transient ribs or sheets appearing between larger, more stable and regular polygonal cells [52,53,68]. Interestingly, however, the wavelength selection process is quite different. For two-sided convection the dominant circulation cells, sampled at the midplane, scale as a ∼ Ra 0.5 . Near the walls the density of the smaller, transient features increases more rapidly with Ra, with a dependence thought to approach a ∼ Ra at high enough Ra [52,53,78]. However, averaging out these fast-moving features over time, or through spatial filtering, shows a near-surface pattern that closely matches the internal convection cells [53].
For our one-sided salt lake model, quantitative details of the scale selection process are presented in Fig. 4, which also includes the neutral stability curve and most unstable mode of convection, as derived in Ref. [64]. Simulations with Ra > Ra c are unstable to convective overturning, which becomes more vigorous with increasing Ra. To begin with, pattern wavenumbers were calculated as in [68], as the first moments of the radially-averaged power spectra of the salinity distributions at a depth of Z ≈ 1. The initial instability was characterized at the rescaled time τ = 1. At these early times the spacing of downwellings in both the 2D and 3D simulations (red markers) agrees well with the most unstable mode of the linear stability analysis. When measured at τ = 30, however, many downwelling plumes have merged, leading to smaller wavenumbers (blue markers). Confirming that the system is in a statistically steady state by this time, simulations run to τ = 60 show no measurable systematic drift in k (values match with τ = 30 to ∆k = −0.004 ± 0.01, mean ± standard error). Similarly, additional simulations with h = 20 gave comparable results (values match with h = 10 to ∆k = −0.04 ± 0.04).
In the steady state, there is a gradual increase in the average wavenumber with Ra, as is the case also for 2-sided convection, although this trend is systematically weaker than the scalings evident in that system [52,53,68]. This increase appears to be driven by the proliferation of small, transient downwelling microplumes as Ra increases. These intermittently appear, then are swept into the more established large-scale convective features (see supplemental movies S2, S3). As the surface crust is built up over time, we would expect it to reflect the more stable features, i.e. megaplumes. To isolate the size of these features, we calculated the dominant wavenumber,ā, associated with the maximum of the power spectrum. As shown in Fig. 4 (turquoise markers), this predicts a pattern of polygonal cells, bounded by salt-rich downwellings, and a characteristic wavelength that consistently settles to values near to Λ = 2π/a c , or about 1.3 m for the conditions typical of Owens Lake.

C. Salt flux and ridge growth
At the surface of the soil water will evaporate and leave behind its dissolved burden of salt. In our model of this process the rate at which salt is added to the surface depends on its advection along with the water, but is also influenced by diffusion down concentration gradients away from the crust. In other words, the mass balance of Eq. (5) can be written in terms of a salinity flux where ∂S/∂τ = −∇ · J S . At the surface, Z = 0, the upwards water flux matches the evaporation rate, so U Z = 1. The presence of a salt crust also sets S = 1 there. As such, the salinity flux into the crust is given by This flux vanishes, J Z = 0, in the stationary solution S = exp(Z), although this scenario will still leave a constant upward flux of salt into the crust at whatever concentration is supplied by the reservoir, as we will describe

Rayleigh number Ra
Pattern wavenumbers a, as measured in twodimensional (circles) and three-dimensional (squares) simulations for different Ra. Average wavenumbers a are measured at two times: once early into the growth of the initial instability, at τ = 1 (red), and later at τ = 30 (blue), after the pattern has had time to develop towards a statistically steady state. We also report the wavenumber of the dominant modeā (turquoise), at late times. For the 2D simulations, the marker positions and error bars represent means and standard deviations as measured in ensembles of 6-10 realizations each.
shortly. Variations in the near-surface salinity gradient will then modulate crust growth around this reference case. Where vertical gradients are strong crust growth will be suppressed (J Z < 0), whereas over weaker gradients crust growth will be enhanced (J Z > 0).
To demonstrate this modulation quantitatively, we calculated the surface salinity flux according to Eq. (9) in simulations at Ra = 100, with results for 2D shown in Fig. 5, and for 3D shown in Fig. 3(d, e). Supplemental movies S2 and S3 show the full dynamics of the two-and three-dimensional systems, respectively. Above downwellings plumes (2D) or sheets (3D), the salinity flux into the surface is positive and shows marked peaks, whereas between downwellings the salinity flux into the surface is negative, and roughly constant. As sketched in Fig. 2, this predicts the development of narrow regions of faster salt precipitation above any convective downwellings, which we argue gives rise to ridges there.
As visible in the movies S2 and S3, the positions of downwelling plumes slowly fluctuate as the simulations progress, as does the pattern of flux into the crust. The motion of the more well-established plumes is gradual, especially in 3D, but the presence of large ridges at some field sites suggests that feedback mechanisms may be needed to help stabilize the patterns over long periods of time. Such feedback could act through a spatial modulation of the evaporation rate E based on salt coverage, which then acts to pin the position of the downwellings in space. Reference [64] gives some discussion on how such a feedback mechanism could work.
We also note that the dynamics in a sufficiently deep lake should approach the case where the average surface salinity flux into the crust vanishes. This follows from the divergence theorem. For a deep lake the descending plumes will mix diffusively with upwelling fluid over a timescale ∼ Λ 2 , for plume spacing Λ. Since the plumes start their fall from the surface at a characteristic velocity V B ∼ Ra [64], this mixing will occur over depths of order RaΛ 2 . As they mix and lose buoyancy the plumes will slow, setting a maximum depth that they can penetrate, of the same order. Below this depth the average upwards salinity flux of the flows from the reservoir must approach zero. Since in a statistically steady state the total salinity in the convecting part of the domain is constant, or fluctuating around a well-defined average value, the salinity flux into the surface must also approach zero. Probing this limit accurately, however, requires domain sizes well beyond the scope of our current simulations.
To make predictions of the crust dynamics at any particular field site, the salinity flux J Z needs to be converted into estimates of the real salt flux into the crust. Assuming minimal volume change due to mixing, the density of the groundwater can be written as ρ = ρ w +βc, where ρ w is the density of pure water, c is the mass concentration of dissolved salts in the solution, and β is a coefficient of expansion. Analogous to the salinity flux in Eq. (8) the mass flux of salt, in dimensional variables, is where the second equation follows by changing variables using ρ = ρ 0 + S∆ρ. As can be seen, there is an offset between the two fluxes-representing the fact that the groundwater feeding into the system from a distant source (S = 0) contains some initial level of dissolved salt, c 0 . In dimensional terms, the predicted salt mass flux into the crust is now We will apply this to Owens Lake in Section IV C.

D. Model summary & predictions
We have developed here a model of porous medium convection beneath a salt lake, which can accurately predict much of the known phenomena of the crust patterns in dry salt lakes, as summarized in Section II. The driving mechanism of salinity-driven convection and its associated governing equations are drawn from existing models of salt lake dynamics [44,45,64] and observations of convection beneath tidal sabkhas [46,65], but we extend these to consider interactions with the surface crust. Evaporation at the surface causes a vertical salinity gradient to emerge, which leads to salt precipitation and crust growth. If a critical Rayleigh number is surpassed, convective cells emerge and modulate the salinity flux into the surface, with faster precipitation expected over downwellings and slower crust growth over upwellings. We will provide field observations in Section IV to quantify these predictions, and show that they are consistent with average growth rates on the order of millimeters per month, and about a two-fold faster deposition of salt at downwellings, relative to upwelling regions. Observing the pattern structure in 3D, we found that over sufficient time the downwelling features arrange into narrow sheets that delimit a closed polygonal pattern of convection cells. The scale of these polygons, in a steady state, does not have any strong dependence on the Rayleigh number. In particular, from Ra = 20 to 1000 we showed that the average wavenumber of the near-surface convection pattern changes by at most a factor of four; the dominant wavenumber, which reflects the structure of the more stable, deeper convection cells, does not evolve significantly over this range and remains consistent with the critical wavenumber, a c = 0.76. This predicts that the pattern wavelength will be relatively insensitive to soil composition (which determines e.g. permeability), groundwater chemistry (e.g. sulphate rich, or pure halite, affecting ∆ρ and to a lesser extent D) and robust to fluctuations in environmental driving parameters.

IV. EVIDENCE FOR CONVECTION AS DRIVING MECHANISM
In Section II we outlined a set of criteria against which to test mechanisms of pattern formation in salt playa. Briefly, any mechanism should be active in salt lake settings, modulate salt transport by directing more salt towards ridges, account for the timescales of crust growth, explain the consistent wavelength of closed polygonal shapes in the crust and be robust enough to explain all this across the variety of conditions associated with patterned crusts. Here, we will explore these criteria in turn, focusing on empirical evidence from field work, supplemented by analogue experiments. Details of field data collection are given in Appendix B, data processing in Appendices C-H and experimental methods in Appendix I. The results point towards convective overturning as a suitable driving mechanism of salt crust patterns.

A. Convective instability
A driving mechanism should be specific to the geophysical conditions encountered in salt lakes. As our model is based on boundary conditions and governing equations of porous media transport that are appropriate for playa and sabkhas (e.g. [44,45,64,72,75]), we aim to prove here that typical dry salt lakes have conditions that would make them unstable to the onset of convective motion. To this end, we visited and characterized field sites at Owens Lake (CA, USA), Badwater Basin (CA, USA) and Sua Pan (Botswana), to show that they can be described by a Rayleigh number above the critical value of Ra c 14.35 [43,64,72]. All sites show well-developed polygonal patterns, and details of the field methodology are given in Appendix B. These results are supplemented by experiments that visualize buoyancydriven convective dynamics in conditions similar to those of real salt lakes (see Supplemental movie S4).
To determine the Rayleigh number appropriate for each field site, we measured the relevant parameters of Eq. (7). From grain-size distributions (data deposited at [80], analysis in Appendix C) we calculated d s , the Sauter diameter [81], of near-surface soil samples. The results, from 4.3 ± 0.6 µm to 138 ± 23 µm, represent a silt to fine sand. A high soil porosity, ϕ = 0.70 ± 0.02, was previously measured at Owens Lake [37]. From these, the relative permeability was estimated using the empirical relationship κ = 0.11 ϕ 5.6 d 2 s , which fits a broad set of experimental and simulation data [82]. Across all sites κ = (3 ± 2) × 10 −13 m 2 to (2.7 ± 1.2) × 10 −10 m 2 . At Owens Lake we also directly measured density differences of ∆ρ = 210 ± 10 kg m −3 between pore water samples taken from close to the surface and at approximately 1 m depth (data at [83]). Annual average evaporation rates of groundwater are taken as E = 0.4 ± 0.1 mm/day [29,37] for Owens Lake, 0.3 ± 0.1 mm/day [38] for Badwater Basin, and E = 0.7 ± 0.5 mm/day for Sua pan [35,39], with further details reviewed in Appendix D. These low levels of groundwater evaporation are characteristic of salt pan environments, and similar rates are seen in active playa of the Atacama desert (E = 0.5 ± 0.1 mm/day [40]) and sabkha near Abu Dhabi (E = 0.2 ± 0.05 mm/day [41]). Finally, we assume the dynamic viscosity of the groundwater to be a constant µ = 10 −3 Pa s. A detailed description of the data sets from Owens Lake and Badwater Basin is in Ref. [11].
From these observations we calculated Ra at twentyone sites around Owens Lake, five in Badwater Basin and seven at Sua Pan. The median values at these three regions were Ra = 3700, 36000 and 420, respectively. Values for all 33 sample locations were between Ra = 120 ± 50 and (1.2 ± 0.5) × 10 5 , well above Ra c . The Sauter diameters, permeabilities and Rayleigh numbers for all sites are given in the Appendix , Table IV. The conditions throughout these patterned salt playa are, therefore, suitable to expect a convective overturning of their groundwater, with plumes of high salinity sinking downwards from the surface.
It is interesting to note that convective plumes of saltrich water have also been observed by electrical resistivity measurements after a heavy rainfall on salt crusts near Abu Dhabi [46] and beneath wind-tidal flats in Texas [65]. While these are a slightly different phenomenon, involving the sudden addition of salt-rich brine to the surface as the rain dissolves salt, they demonstrate convection in similar geophysical conditions to salt lakes. The Rayleigh numbers calculated for these two cases are up to about 40,000 and 90, respectively [46,65].
We complemented these observations with analogue experiments in Hele-Shaw cells to demonstrate buoyancydriven convection in porous media under conditions comparable to the field. These experiments are inspired by Refs. [84][85][86] but instead of using a narrow and empty gap between the cell walls to create the porous medium, we used a 0.8 cm thick cell filled with glass beads. We applied a strong surface evaporation, driven by heat lamps and ventilation, and the base of the cells were connected to fluid reservoirs with a fixed salt concentration of 50 kg m −3 . Further details of the experimental methods are given in Appendix I and an example time-lapse movie of an experiment in progress is given as supplemental movie S4. To observe whether a convective instability arises for different conditions, we systematically varied the grain sizes of the glass beads. This modifies the relative permeability of the system and allowed us to change the Rayleigh number of the experiments. Convection was seen in all experiments above Ra c , but not in the finest-grained case, where Ra Ra c . These results are consistent with those of Wooding et al. [84], who experimentally confirmed theoretical predictions of the onset of convection for the slightly different case of a salt lake with surface ponding, or constant pressure boundary conditions. Finally, we note that the range over which the experiments show convection, between Ra = 20 ± 10 and (1.3 ± 0.3) × 10 3 , also coincides with much of the range of Rayleigh numbers we estimated from the field sites.

B. Mapping salt heterogeneity
If salt polygon growth is driven by convective dynamics happening beneath the surface, then horizontal differences in salt concentration should be detectable in soil and pore fluid under typical field conditions, and also in laboratory-based analogue experiments. Salt variations of this nature have been seen in sabkhas after rain [46,65]. We provide evidence here of salt-rich plumes beneath the crust of Owens Lake, and show that the plumes are correlated with the positions of the ridges in the crust.
First, to demonstrate the coupling of salt concentration and convective motion under controlled conditions, we sampled one of the Hele-Shaw experiments introduced in Section IV A. For this we dissected an experiment that was undergoing convection, and which had been evolving under constant conditions for two months. Samples were extracted from locations along the downwelling and upwelling plumes, as indicated by the motion of dye added into the reservoir fluid a few days before sampling. As shown in Fig. 6, the fluid flow in the analogue experiments is clearly driven by, and coupled to, variations in salinity. Pore water samples taken from the upwelling regions (dyed yellow in the figure) show systematically lower salinity than those taken from downwellings.
Next, from the field we collected samples of wet soil from fresh trenches dug at two unmanaged sites at Owens Lake. Surveys of the surface relief were made before sampling, using a terrestrial laser scanner (TLS, see Appendix E and Ref. [22] for methodology; data deposited at [87]). These show the presence of salt polygons of about 2 m in size delimited by high ridges, as in Fig. 7(a).
In each case we sampled soil along a grid pattern in a cross-section below a polygon, including beneath two bounding ridges. Analysis of the salt concentration of the samples with respect to pore water content shows direct evidence of plumes of high-salinity fluid below the salt ridges ( Fig. 7(b), methods in Appendix F and data deposited at [88]). Specifically, we tested whether the distribution of salt concentrations in an area below each ridge was different to that below the flat pan of the polygon; testing this hypothesis (two-sample KS test), shows that the distributions below ridges and flat crust are statistically distinct (p < 0.02), at both sites.
The measurements of the salt concentration in the pore water also shows an exponential decay in salinity with depth (Fig. 7(c)), consistent with a salt-rich boundary layer that is heavy enough to drive convection in a porous medium. As evidence for this, we recovered representative values for the boundary layer thickness from exponential fits to the horizontally-averaged salt concentrations at both trench sites (see Appendix H). Interestingly, the observed values of 13.5 ± 5.3 cm and 17.7 ± 1.5 cm are comparable to the natural length scale of L = ϕD/E = 15.1 ± 8.0 cm estimated for Owens Lake.
Thus, not only does direct field sampling of groundwater beneath a patterned salt crust show both horizontal and vertical variations in salt concentration, which support the hypothesis that the system is unstable and convecting, but it also demonstrates that the plumes of high salinity are co-localized with the surface ridges.

C. Surface salt flux in the field
A driving mechanism for pattern formation should be able to spatially modulate salt transport to the soil surface, leading to salt ridge growth on a scale of months. We demonstrate this growth at Owens Lake through a time-lapse video (Supplemental movie S1) which shows the surface crust dynamics over a period of approximately four months in spring 2018. At Sua Pan, crust dynamics have been measured by TLS scanning of several sites over the course of several months [15]. Active crusts there showed average growth rates of 0.7 to 1.5 mm/month. The fastest short-term growth rates of ridge features, briefly reaching up to 10 mm/month, were associated with ridge thrusting. In Section III C we explored the corresponding predictions of our model. Specifically, Eq. (12) describes the expected salt flux into the surface, J c,z , depending on the salt concentration gradient there. Here, we apply the framework of this model to the maps of the subsurface salt concentration measured at Owens Lake (Fig. 7). In particular, we will show that the observed variations are sufficient to drive average crust growth rates on the order of millimeters per month along with a preferential growth of the ridges, over the rest of the crust.
To calculate the expected salt flux into the surface, we estimated the near-surface gradient in salt concentration from the field data from our trench sites. Dividing this data between regions beneath the ridges and centers of each polygon, we used an exponential fit to recover the effective thickness of the salt-rich boundary layer L , the background salt concentration c bkg and the salt concentration at the surface c sat , for each such region. For this, the measured weight fractions of salt were re-scaled by the densities of salt solutions, to give concentrations c in terms of the mass of dissolved salt per unit volume of fluid. Further details of the fitting are given in Appendix H, including Figs. 11 and 12.
Substituting the exponential form of the fit into Eq. (12) yields a prediction for the mass flux of salt, in terms of measurable quantities, of Values for L and c bkg below each ridge and polygon center are extracted from the fits. For c sat , the fitted values are all consistent with each other, so a common value of a saturated solution, c sat = 316 kg m −3 , was used to minimize uncertainties. As in Section IV A we use values for the evaporation rate, E = 4.6 × 10 −9 m s −1 , and porosity, ϕ = 0.7, taken from previous observations at Owens Lake [29,37]. For an effective diffusion constant, we use the tortuosity-corrected diffusivity D = 1.00 × 10 −9 m 2 s −1 (see Appendix G). Finally, we measured the density of salt crust samples collected from Owens Lake to be ρ crust = (960 ± 50) kg m −3 . This value is used to convert mass flux rates into crust growth rates. The resulting salt flux and growth rates for the ridges and centers of the two trench sites from Fig. 7 are shown in Table I. The derived salt flux into the surface is about twice as high at the ridges, as compared to the centers of polygons. In order to allow a clear statistical comparison between these rates, only independent sources of errors were propagated (i.e. excluding uncertainties in E, D and ϕ, as these will affect all rates in the same way) to give the main uncertainties reported in Table I. Within these uncertainties, we see a consistent difference between site area Jc,z 10 −7 [kg m −2 s −1 ] r [mm/month] 1 ridge 1 5.6 ± 2.0 (3.9) 1.5 ± 0.6 (1. TABLE I. Salt flux into the surface Jc,z and crust growth rate r inferred at the polygon ridges and centers for the two sites at Owens Lake depicted in Fig. 7. The error bounds consider independent sources of uncertainty, while the values in brackets give ranges taking into account all sources of uncertainty (i.e. including lake-wide uncertainties in E, D and ϕ).
the crust growth rates at polygon ridges when compared to polygon centers. As salt is being added faster to the ridge areas, we would expect differential growth, focused on the ridges. Although the crust morphology can become quite complex [12,19,89], the faster expansion of the crust in these locations can plausibly contribute to differential strains, and features like ridge thrusting or local buckling. This is consistent with the faster dynamics seen at ridges in Sua Pan [15] and that we report at Owens Lake (see supplemental movie S1). Over long times the average crust growth rate will be essentially constrained by the evaporation rate and the amount of salt in the groundwater feeding into the lake. However, it would not be surprising if these rates varied with time, e.g seasonally, or via extreme events like surface flooding. As we have not accounted for the excess mass of the hydrated states of any of the crust materials (especially mirabilite, present as a significant minority species of salt at Owens Lake and Sua Pan, although not noticeably present at other pans, such as Badwater Basin), the rates we report in Table I also represent a lower bound on the anticipated crust growth rates.
We have shown here how the heterogeneous salinity distribution measured below salt polygons in the field infers variations in the salt flux to the surface, if the subsurface fluids are undergoing porous medium convection. Within the framework of our convection model, the measured pattern of salt concentration near the surface of the crust imply about a twofold increase in crust growth rate at the ridges, rather than the centers of salt polygons. The crust growth rates on the order of millimeters per month are also consistent with the time scale of crust growth observed in the field, which show that noticeable changes can occur over months.

D. Pattern length scale
If groundwater convection leads to preferential locations for salt precipitation, and from thence to salt crust patterning, then the convective cells and crust polygons should have similar length-scales and patterns. We presented model predictions for the length-scale selection of convection in Section III B. To test these predictions in the field, we measured the surface relief of the crusts at multiple sites from all three dry lakes using a terrestrial laser scanner (see Appendix B). The crusts display the typical polygonal patterns of salt lakes, as shown in Fig. 8(a-d), with further data in [11,15]. From each scan we extracted a characteristic wavelength for the pattern, based on the dominant frequency of the power spectrum, using the methods reported in [15]. Full results are given in the Appendix, Table IV. In all cases the pattern wavelengths ranged from λ = 0.53 ± 0.20 m to 3.02 ± 1.40 m. These were converted into dimensionless wavenumbers a = 2πL/λ, using the characteristic lengthscale L = ϕD/E as calculated for each lake (E from Section IV A, ϕ and D from Section IV C). The results are summarized in Fig. 8(e), which shows the Rayleigh number and wavenumber for each site studied.
We found that the wavenumbers characterizing the surface crust patterns in the field are clustered near the critical wavenumber for convection, a c = 0.76, and largely independent of Ra. In Section III B, we predicted how the statistically steady state of subsurface convection should look. In this state the average wavenumber increased slowly with Ra, by about a factor of four from Ra = 20 to 1000. This was attributed to the proliferation of small, transient features, superimposed on more stable convection cells. The dominant wavenumber, given by the peak of the power spectrum, showed little-to-no dependence on Ra, and was argued to better represent the long-term behavior of the convection cells. As shown in Fig. 8(e), although there is large scatter in the data, the wavelengths of the crust pattern broadly match the wavelengths of the dominant features of the sub-surface convection expected in a statistically steady state.
Furthermore, for the field sites there should be ample time for the convective dynamics to develop towards the steady-state condition. In the model, the timescale most appropriate to characterize the dynamics of convection is T / , where T = ϕ 2 D/E 2 describes how the boundary layer would develop in the absence of convection, and = (Ra − Ra c )/Ra c accounts for the increased vigor of convective motion as the driving forces increase (see also [64]). For the relevant parameter values of Owens Lake, T / is just over one day. Alternatively put, one year of evolution at Owens Lake, for a site with a representative value of Ra = 3700, would represent a rescaled time of τ = 350. This is long compared to the timescale of τ = 30 that we used to characterize the steady-state conditions in Figs. 4 and 8.
Finally, as shown in Fig. 8(a-d), polygon ridges are relatively narrow features. To characterize their size, we thresholded our surface relief maps using Otsu's method [90]. The relative area occupied by ridges, measured in this way, was 24 ± 8% (mean ± standard deviation); this may be a slight overestimate, due to occlusion of the TLS behind larger ridges. In the 3D simulations the comparable features are downwellings, which we characterized as areas of positive surface salinity flux at τ = 30. These covered 11-18% of the surface; narrower features were seen at higher Ra, but the range of Ra explored is small. Although the salt ridges can be relatively complex, involving breaking and thrusting for example, the convective model qualitatively captures the essential asymmetry between the narrow ridge features and the flat pans between them.
In summary, the pattern wavelength we measured for salt polygons at a range of sites across three different salt lakes is thus found to be consistent with the length scale expected for buoyancy driven convection going on in a statistically steady state. Simulations of the dynamics in a three-dimensional domain (see Fig. 3(d)) also show that the convective mechanism can produce patterns of closed polygonal shapes at this lengthscale. These patterns are qualitatively similar to the salt polygons observed in the field (Fig. 8(a-d)) with relatively narrow features delimiting the larger polygonal shapes.

E. Environmental variation & remaining challenges
To this point, we have focused on the effects of varying environmental conditions as they contribute to the Rayleigh number and natural lengthscale of the convection problem. We showed that a wide range of measured field conditions should lead to convection, and to convection cells with a size and pattern that is consistent with the observed surface crust patterns; faster crust growth, and thicker accumulations in ridges, was shown to correlate with the downwellings of the subsurface convection. In this section, we will briefly look at and discuss other factors, including the salt chemistry, crust thickness and features like cracks and crust wrinkling. We will also discuss some of the remaining assumptions and challenges of the convective model, in explaining crust patterning.
In addition to the soil, we characterized the composition of the salt crust at Owens Lake and Badwater Basin (see Appendix G for methods; data deposited at [83]). Both sites are predominantly halite, although Owens Lake additionally has significant concentrations of carbonates and sulfates [37,83]. Between all the field sites we also observed salt crusts with thicknesses from below a centimeter to approximately 30 cm. Across these differences in soil and salt crust composition there is a similar appearance of the salt polygons. This agrees with a mechanism like subsurface convection, which is relatively insensitive to salt chemistry. It would, conversely, present difficulties to any mechanism that was reliant on specific salts, such as mirabilite, or which varied strongly with crust thickness.
Other surface features, including cracks and wrinkling or buckling, are also associated with salt crusts [12,17]. While studying Owens Lake, we observed well-ordered hexagonal desiccation cracks in crust-free mud with approximately 10 cm spacing, as shown in Fig. 9(a). These appeared at the crust-free edge of a pan that was otherwise covered in salt ridges (with a spacing of 1-2 m). The hexagonal crack patterns can be explained as the result of the intermittent and episodic cracking of the soil, as water levels fluctuate [51]. However, they form at a very different scale to the salt ridges, and the cracks develop as localized depressions in the ground, rather than upthrust ridges. Further, despite their presence near the salt crust, there were no signs of preferential precipitation of salt in the cracks. We also observed surface buckling of the salt crust at Owens Lake, as shown in Fig. 9(b) (see also supplementaal video S1 on dates 05/07 through 05/11). These features appeared when the crust was solidifying, and may result from the stresses induced by phase changes or the addition of salt mass into the crust. However, the scale of their features is, again, much smaller than those of the dominant polygonal ridges in the crust, and the patterns formed are also different. Finally, we note that throughout our core samples and trench sites, as in Fig. 9(c), the soil horizons remained normal and undisturbed by the presence of surface ridges. This would not be the case if deep fractures were occurring along with the ridge formation (as in e.g. ice-wedge formation in permafrost terrain [48]).
More generally, we have treated the salt crust as twodimensional, of nominal thickness. This is a necessary simplification, given the novelty of the hypothesis advanced here. However, there will be some internal dynamics within the salt crust itself, beyond the scope of our model. For example, dissolution will be favored at the base of the crust, in contact with groundwater, and precipitation at the top, where evaporation occurs. Similar observations have been made in playa [89] and in lab [91,92] studies of efflorescence with cultural heritage applications. Dissolution-precipitation dynamics raises porosity in crusts, leading to a 'fluffy' texture, trapped air pockets, self-lifting effects [89,91,92], and a lower threshold for dust entrainment [24]. We accounted for the high crust porosity while calculating growth rates in section IV C, but a more nuanced approach may give further insight, including into dust models.
There is also likely to be an element of positive feedback between the growth of the salt crust, and the dynamics of the convective plumes. Our simulations showed that, without such a feedback, the positions of downwelling features slowly fluctuate in space, although this tendency is weaker in 3D than in the 2D simulations. A feedback cycle, where evaporation rates are affected locally by the presence of the ridges, can act to stabilize the locations of the plumes. This has been demonstrated [64] by simulations in which downwellings tend to move towards, and stay at, locations with lower evaporation rates. We have also measured temperature and relative humidity records from Owens Lake [64], which show that evaporation rates should be slower at the ridges, compared to the centers of polygons, in a way consistent with such a feedback mechanism. Further model development along these lines would be promising.
Incorporating feedback into crust evaporation may also help refine predictions of wavelength selection. The scale of salt polygons matches the dominant wavelength of our simulations, which we argued to be representative of a long-time integration of the more stable convection pattern. The effect of the smaller-scale, transient 'microplumes' remains uncertain. As captured by the average wavenumber, these are more prevalent at higher Ra, but their faster dynamics (changing visibly over simulation times of ∼ 0.1, or about two weeks in real terms) mean that they should only weakly contribute to the surface pattern, which develops over months. This prediction could be checked in simulations that modulate surface evaporation, based on tracking of the total amount of salt accumulating at the surface over time.
As a final point, one may also consider thermal contributions to buoyancy, and phenomena like doublediffusive convection. A typical 10 • C diurnal variation [64] will change water density by about ∆ρ T = 1 kg m −3 , compared to the ∆ρ S 200 kg m −3 density contrast due to solutes at Owens Lake. The buoyancy ratio, N = ∆ρ S /∆ρ T 200, means that the driving forces and speeds of thermally-driven flows will be reduced by a factor of 1/N , compared to solute-driven flows (see e.g. [93] for review of double-diffusive flow in porous media). Such effects can, therefore, be assumed to be negligible.

V. DISCUSSION AND CONCLUSION
Salt deserts, playa and pans are a common landform important for climate balances such as dust, energy and water, and express a rich repertoire of patterns and dynamics. Here we have shown that, in order to model and understand the surface expression of such deserts, insight is gained by considering the crust together with the subsurface dynamics. In particular, we have shown how the emergence of regular salt polygons, which are a common salt crust pattern, can result from their coupling to a convection process in the soil beneath them. The existence of salinity-driven convection in salt pans is itself already an anticipated finding (see e.g. [44]). Furthermore, a salinity gradient below the surface, strong enough to drive such convection, has been measured at Owen's Lake [37] and strong horizontal salinity gradients indicative of convective motion have been shown to exist in sabkhas [46,65]. What we have demonstrated is that this convection process develops into a statistically steady state characterized by large-scale convection cells with a wavelength of a few meters, after a few weeks of evolution under typical field conditions. We further showed how this can explain the remarkably consistent appearance, growth rates, polygonal shape, narrow ridges, and representative wavelengths of the salt crust patterns observed in different parts of the world, including our field sites at Sua Pan, Badwater Basin and Owens Lake. Our findings help to elucidate the relationships between subsurface drivers and the development of surface crust patterns. Better knowledge of crust development can improve dust emission models [7,24,31], including quantifying the connection of salt pan hydrology to controls on dust emission [33] such as those used at Owens Lake [29,30]. As further context, dusts and aerosols remain one of the larger uncertainties in modelling climate sensitivity [94], with e.g. feedback from soil biocrust degredation recently predicted to contribute half as much as direct anthropogenic aerosol emissions by 2070 [95]. While the surface variability of crusted surfaces can be probed using remote sensing [23], our model is a first step to improving predictions of how subsurface hydrologic variability links to surface crust variability. This knowledge is increasingly important as saline lakes like the Dead Sea or Great Salt Lake are shrinking [96] and their newly exposed beaches are candidate areas for salt crust formation and dust emission.
To establish the connection between surface features and subsurface flows, we demonstrated consistent results from theoretical and numerical modeling, analogue experiments and field studies. In contrast to previous explanations of crust patterns [12,17], this model is able to explain the robustness of the pattern length scale by considering the statistically steady state of porous media convection, based only on measured environmental parameters. In fully three-dimensional simulations, the convective dynamics were also shown to give rise to closed-form polygonal shapes. More importantly, it is able to predict a suite of quantitative details, such as the timescales required for convection to set in, the rates of crust and ridge growth, the positions of the ridges above salt-rich plumes, and the relative narrowness of the ridge features. At the downwellings the salinity is higher and therefore the salinity gradient between the crust and the underlying fluid is weaker (compare sketch of Fig. 2 to measurements in Fig. 7). As salt transport is a balance of advective and diffusive processes, this will lead to an increased rate of salt precipitation above downwelling features, contributing to the growth of ridges at the boundaries of convection cells. After the initial emergence of ridges, the growth process might be bolstered by feedback mechanisms such as a modulation of the evaporation rate by the presence of ridges, cracks or surface wicking phenomena.
As such, our results show how salt polygons are part of a growing list of geophysical phenomena, such as fairy circles [97], ice wedges [49], polygonal terrain [50] and columnar joints [51], which can be successfully explained as the result of the instability of a dynamical process.

VI. DATA AND CODE AVAILABILITY
Field site labels and locations as well as site-by-site data for Fig. 8(e) are given in the appendix Table IV. Full data sets for the soil particle sizes [80], surface height profiles [87], salt chemistry [83], subsurface salinity distribution and pore water density [88] as well as images describing the US field sites [98] are available under a CC-BY-SA 4.0 license at PANGAEA.
Code for the 2D simulation is available under an MIT license on GitHub, doi:10.5281/zenodo.3969492.
All scripts, data sets and images used to produce the figures in this work are available at https://doi.org/ 10.17605/OSF.IO/KJDWF.

VII. AUTHOR CONTRIBUTIONS
JL, ME and LG wrote the original draft of the manuscript; all authors reviewed the manuscript.
JL participated in fieldwork at Owens Lake and Badwater basin, analyzed the soil salinity and grain size distributions, performed numerical experiments of the 2D simulation, performed the Hele-Shaw experiments, analyzed all data except for the 3D simulations and TLS measurements and created the data visualizations.
JMN participated in fieldwork at Owens Lake, Badwater Basin and Sua Pan and performed and analyzed the TLS measurements. ME developed the theory of convective dynamics and programmed the 2D simulation.
VK performed the X-ray analysis, supervised JL for the grain size measurements and provided feedback for the soil composition and salt species analysis.
GFSW participated in fieldwork at Sua Pan. MRT performed numerical experiments and analyzed the data from the 3D simulations.
CB supervised MRT, programmed the 3D simulation and provided feedback for the theoretical model.
LG conceptualized the study, developed the theory of convective dynamics, supervised JL and ME, performed the water density measurements and participated in fieldwork at Owens Lake and Badwater Basin. The three-dimensional (3D) simulations use a pseudospectral element method where the horizontal directions, over which the numerical domain is periodic, are treated using fast Fourier transforms, and the vertical (bounded) direction is decomposed into elements, each of which is treated using Gauss-Lobatto-Legendre quadratures (see e.g. [99]). The time-derivative is approximated using a second-order backward differentiation formula and the nonlinear term extrapolated using a Taylor expansion of second order. The coefficients of the resulting schemes can be found in [100]. All simulations were carried out on a domain of depth 10 and horizontal area of 24π×24π.
To compute wavenumbers using two-dimensional simulations (cf. Fig. 4), we used the same numerical solver and imposed the solution invariance in one of the horizontal directions. For the two-dimensional (2D) simulations of wavenumber selection (Fig. 4), we used the same numerical solver and imposed the solution invariance in one of the horizontal directions. These simulations have horizontal domain size 12π and depth 10 or 20. Pattern wavenumbers were calculated from the salinity at the gridpoint closest to Z = 1 as in 3D, but without radiallyaveraging the power spectra. The dominant wavenumber was identified as the mode with the highest power spectral density. All other two-dimensional (2D) simulations were carried out using a stream-function-vorticity approach based on Refs. [101,102] with a detailed implementation described in Ref. [64].
Appendix B: Field data collection Field work was conducted at Owens Lake and Badwater Basin in November 2016 and January 2018; see e.g. [103,104] for geological descriptions of the area. The Owens Lake sampling sites are indicated in Fig. 10. At Badwater Basin five sites were visited ∼500 m south of the main tourist entrance to the playa. For both locations field methodology and data collection are fully documented in a separate data description publication [11].
At Owens Lake we use field site labels referring to surface management cells of the dust control project there [105]. These labels link either to managed cells or to unmanaged areas in the direct vicinity of a managed cell. Labels typically start with "TX-Y", where X is a number and Y is a number or letter. The first number refers to the water outlet taps along the main water pipeline that crosses the lakebed south to north and is used to irrigate management cells. Numbers/letters after the hyphen refer to subregions branching from the same tap. For some sites we investigated more than one polygon. This is indicated with brackets, e.g. T27-A (3) is the third polygon investigated at site T27-A, which corresponds to the Addition (A) region of the management cell next to the 27 th tap. Bracketed numbers are used for the sites visited at Badwater Basin.
To evaluate the profile of salt concentration below polygons, at most sites soil cores (4 cm Dutch gouge auger) were taken to a depth of up to 1 m. The soil showed normal bedding (see also Fig. 9(c)), indicative of sedimentation following flooding. Samples were collected from each visible soil horizon, or with a vertical resolution of ∆z =10-15 cm.
Trenches were also dug at sites T27-S (Site 1 in Figures 7 and 11) and T32-1-L1 (3) (Site 2), in order to take samples along cross-sections below salt polygons. The trenches were dug about 200 cm in length, 40 cm in width and down to a water table of ∼70 cm. Soil samples were taken from a freshly cleaned trench wall in a grid pattern with spacings of ∆x = 15 cm and ∆z = 10 cm. An example of such a trench along with sampling locations is shown in Fig. 9(c). The samples had an average volume of approximately 10 ml and were taken using a metal spatula, which was cleaned with distilled water and dried before each use. The samples were a mixture of soil with a grain size of medium sand to clay, water and salt (both dissolved and precipitated). After collection, samples were immediately stored in air-tight containers, which were sealed with parafilm.
To evaluate the density difference, ∆ρ, we collected pore water samples at Owens Lake from eight sites, including liquid taken from directly below the surface (in cases where the water table was at or within 10 cm of the surface), and at a depth of about 1 m.
Samples from Sua Pan were collected during a different field campaign, reported in Ref. [15]; site labelling follows that publication. Sediment samples at Sua Pan were collected 2 cm below the crust in August 2012. These were double bagged and subjected to grain size analysis only.
At all three dry lakes, surface height maps were collected with a Leica Terrestrial Laser Scanner (TLS); a P20 model was used at Owens Lake and Badwater, and a Scanstation at Sua Pan. The scanner head was positioned at least 2 m above the playa surface and scans were performed before the surface was disturbed by sampling.
GPS locations of all sites are provided in Tab. IV.

Appendix C: Soil characterization
Soil samples from Owens Lake, Badwater Basin and Sua Pan were analyzed to determine their distribution of grain sizes. This analysis was performed using the soil remaining after a sample's salinity was determined (Appendix F). Soil grain size distributions were measured by laser particle sizer (Coulter LS 13 320), from which the Sauter diameter (the mean diameter, respecting the soil's specific surface area [81]) was calculated. For each site a representative d S is calculated as the average Sauter diameter of all soil samples (for trenches, one sample per depth) from that site. For Sua Pan, only samples from sites B7 and L5 were available; d S for the other five sites is estimated as the mean of the measured values at these two sites. All grain size data is deposited at Ref. [80].
Soil porosity has previously been measured to be around ϕ ≈ 0.70 ± 0.02 [37] at Owens Lake. Because of lack of similar measurements at Badwater Basin and Sua Pan, we used the value measured at Owens for calculations of κ and Ra at these sites. For each site a permeability is then calculated based on the Sauter diameter and porosity, as κ = 0.11ϕ 5.6 d 2 S [82]. The Sauter diameter d S , permeability κ and Rayleigh number Ra for all sites investigated are provided in Table IV. Uncertainty ranges for Ra were calculated as systematic errors based on standard errors of all the input environmental parameters.

Appendix D: Evaporation rate
Groundwater evaporation rates are taken from the literature. At Owens Lake, Tyler et al. [37] combined lysimeter, Bowen ratio and eddy correlation methods, with repeat seasonal measurements conducted well after any rainfall. The value of E = 0.4 ± 0.1 mm/day captures the range of measurements at their NFIP site, near our northern sites at the lake. At Badwater Basin, De-Meo et al. [38] measured evaporation by Bowen ratio and eddy correlation methods, with continuous monitoring over several years. The value of E = 0.3 ± 0.1 mm/day is taken from their two playa sites south of Furnace Creek, and the 2001 conditions most consistent with the wetter season preceding our field work. At Sua Pan we use E = 0.7 ± 0.5 mm/day, estimated by remote sensing and energy balances [39], and wind tunnel experiments [35].
Appendix E: TLS data processing TLS Scan data was processed following Ref. [22]. Data were first gridded into a digital elevation map (DEM) with a lateral resolution of 1 cm and a vertical resolution of 0.3 mm. Dominant frequencies of surface roughness, which we use here for the pattern wavenumber a (and wavelength λ = 2π/a as reported in Table IV), were then quantified using the 90 th percentiles determined with the zero-upcrossing method from the DEMs [15]. TLS data from the US field sites are deposited at Ref. [87].
Appendix F: Salinity and density measurements Owens Lake and Badwater Basin soil samples were analyzed to determine the amount of salt in their pore fluid. These samples had been sealed immediately after collection, to preserve their water content. After unsealing, each sample was first transferred to a crystallizing dish and weighed, to give a combined mass of sand, salt and water. It was then dried at 80 • C until all moisture had visibly vanished, or for at least 24 h, and re-weighed to determine the mass of the (evaporated) water. Next, it was washed with 50 ml of deionized water to dissolve any salt, allowed to sediment for 24 hours, and the supernatant liquid was collected in another crystallizing dish. After two such washings the remaining soil and the recovered salt solution were dried and weighed. Measurement uncertainty is based on the difference between the initial sample mass and the sum of the separated water (m w ), salt (m s ) and soil masses. This gave a direct measure of the mass fraction of salt in solution as C = m s /(m w + m s ), reported as a weight % (wt.%).
Owens Lake pore water samples were analyzed to determine their density using a vibrating-tube densitometer (Anton Paar DMA4500). The near-surface (0-10 cm depth) groundwater density was consistently 1255 ± 8 kg m −3 , while water from 70-100 cm depth had a density of 1050 ± 2 kg m −3 . These values are broadly consistent with chloride concentration profiles previously measured at Owens Lake [37]. We note that thermal effects on the groundwater density will be comparatively negligible, as the mean annual variation in temperature at Owens Lake will allow for a density change of, at most, 5 kg m −3 . Similarly, the solubility of halite in water would change by less than 5 kg m −3 , seasonally. Data from salt concentration and pore water density measurements are deposited at Ref. [88].

Appendix G: Salt chemistry and diffusivity
Pore water samples from selected sites at Owens Lake and Badwater Basin were analyzed to determine the dominant salt species present. This analysis was performed on the dried salts remaining after the salinity measurements (Appendix F). Mineral identification was performed by quantitative X-ray powder diffraction analysis (Philips X'Pert MPD PW 3040). Samples from Owens Lake showed a mixture of salts with 53 ± 7 wt.% sodium chloride and 30 ± 5 wt.% hydrated sodium sulphate (mirabilite). Other minerals, such as natrite, sylvite and burkeite, were variously present at less than 10 wt.%, each. As such, in solution Na + and Cl − ions predominate (>70% by mass), with SO 2− 4 , CO 2− 3 and K + ions present in descending order of significance. All data are deposited at Ref. [83].
Based on the ionic species found in the pore water, we estimate an average aqueous diffusivity of D * = 1.37 × 10 −9 m 2 s −1 from measurements of ternary  TABLE II. Boundary layer thickness L , along with the salt content (wt.%) of background (C bkg ) and saturated (Csat) pore fluid, from exponential fits to the salt concentration distributions below the center and ridge areas of two trench sites.