Stability and dynamics of convection in dry salt lakes

Abstract Dry lakes covered with a salt crust organised into beautifully patterned networks of narrow ridges are common in arid regions. Here, we consider the initial instability and the ultimate fate of buoyancy-driven convection that could lead to such patterns. Specifically, we look at convection in a deep porous medium with a constant throughflow boundary condition on a horizontal surface, which resembles the situation found below an evaporating salt lake. The system is scaled to have only one free parameter, the Rayleigh number, which characterises the relative driving force for convection. We then solve the resulting linear stability problem for the onset of convection. Further exploring the nonlinear regime of this model with pseudo-spectral numerical methods, we demonstrate how the growth of small downwelling plumes is itself unstable to coarsening, as the system develops into a dynamic steady state. In this mature state we show how the typical speeds and length scales of the convective plumes scale with forcing conditions, and the Rayleigh number. Interestingly, a robust length scale emerges for the pattern wavelength, which is largely independent of the driving parameters. Finally, we introduce a spatially inhomogeneous boundary condition – a modulated evaporation rate – to mimic any feedback between a growing salt crust and the evaporation over the dry salt lake. We show how this boundary condition can introduce phase locking of the downwelling plumes below sites of low evaporation, such as at the ridges of salt polygons.


Introduction
This study of convection in a porous medium is motivated by the patterns shown in figure 1. The examples shown there are of dry salt lakes, or playa (Briere 2000), which are amongst the most inhospitable places on the surface of the earth. Dry lakes typically develop in arid environments, where evaporation outweighs precipitation and where mineral-rich groundwater is refreshed by inflow from surrounding regions of higher altitude (Lowenstein & Hardie 1985;Gill 1996;Briere 2000). The heat fluxes through the surface of such salt deserts are important to understanding the water and energy balances in arid regions (Bryant & Rainey 2002). Additionally, salt deserts are responsible for a significant part of the global emission of atmospheric dust (Gill 1996;Prospero 2002;Washington et al. 2003). However, despite the extreme conditions that prevail above ground, the water table of dry lakes often remains very near to the surface (Gill 1996;Briere 2000;Bryant 2003;Nield et al. 2015), allowing for active patterns of fluid flows within the pore spaces of the soil (e.g. Tyler et al. 1997;Wooding, Tyler & White 1997;Stevens et al. 2009;Van Dam et al. 2009). As evaporation rates are high DeMeo et al. 2003;Brunner et al. 2004;Groeneveld, Huntington & Barz 2010), soluble salts accumulate in such regions, and precipitate into a solid salt crust covering the desert floor. We have argued that these two processes, of subsurface flows and surface crust growth, are coupled (Lasser 2019;Lasser et al. 2019). In the main body of the present study we will focus on analysing the instabilities of the subsurface flow. Subsequently, we will discuss how this flow might support and interact with preferential precipitation of salt in certain areas on the surface.
Within the crust of a dry lake, captivating and beautiful patterns can emerge, developing into a network of polygons, as shown in figure 1. Around the world -from Owens Lake and Badwater Basin in California  or the Salt Lake Desert of Utah (Christiansen 1963) to the Great Salt Desert of Iran (Krinsley 1970), Salar de Uyuni in Bolivia and the Sua Pan of Botswana (Nield et al. 2015) -these salt polygons are usually expressed with a diameter of a couple metres, individually bounded by ridges a few centimetres high. These regular patterns immediately draw the eye of an observer and the question of their origin arises. The raised structures of polygonal ridge patterns in salt deserts also contribute to surface roughness (Nield et al. 2015). Under the effect of the often strong winds which blow over a desert's surface, dust is emitted from salt pans and is carried into the atmosphere. The surface roughness, alongside the salt chemistry and other crust characteristics, influences the uncertainty in modelling dust emission from desert landscapes (Raupach, Gillette & Leys 1993;Marticorena & Bergametti 1995). Christiansen (1963) and Krinsley (1970) attempted to explain the growth of salt polygons in dry lakes by the folding and cracking of the salt crust, respectively. However, neither of these explanations is sufficient to explain the emerging polygonal shapes and, in particular, the robust length scale observed in nature around the world. Specifically, both models would predict that the pattern wavelength is proportional to the thickness of the salt crust. However, this wavelength is consistently 1-3 m, in crusts ranging from sub-centimetre to several metres thick (Krinsley 1970;Lowenstein & Hardie 1985;Lokier 2012;Lasser et al. 2019). Recently, buoyancy-driven convection, taking place within the wet porous sand below the salt crusts, was brought forward as an alternative candidate for a driving mechanism for pattern formation ). Here, we will model the onset of this mechanism as well as the maturation and scaling of the dynamics of buoyancy-driven convection that can occur below the surface of a dry salt lake. Specifically, we will present a study of a model of solutal convection in a porous medium, as illustrated in figure 2. The porous medium can be interpreted as the sediment below the salt crust of a dry salt lake, and the solute as the salt dissolved in the groundwater that fills the porous medium. The salty groundwater is re-supplied by an influx of cleaner water from far below. Evaporation of water is enhanced both by wind and high temperatures, causing the precipitation of salt and growth of the crust at the surface. Due to the accumulation of salt near the surface the salinity, and thereby the density, of the water is higher there than it is further below. If the resulting density imbalance is large enough, this configuration is unstable, leading to buoyancy-driven convective motion.
We model this problem in an idealised and simple two-dimensional geometry, whose domain is deep compared with the dynamics that can arise near the surface. The domain has an upper boundary that is only permeable to fluid and which accounts for fluid loss through evaporation. The lower boundary provides for the recharge of fluid from a reservoir of fixed salt concentration. This model is investigated through a linear stability analysis, as well as with numerical simulations. In the latter case we use periodic boundary conditions in the horizontal direction, and a lower boundary with a constant flux condition. Both situations are given an initial condition in the form of a boundary layer of high salinity, formed just below the top boundary, in which salt diffusion is balanced by an upward advection of fluid. This set-up is designed to mirror the conditions found below an evaporating salt lake, such as Owens Lake or Badwater Basin. We investigate the onset of convective motion in the system at Rayleigh numbers Ra close to the critical value, Ra c , as well as the behaviour at higher Rayleigh numbers. Here, once the system is appropriately scaled, Ra is the only free parameter and can be interpreted as the dimensionless ratio between buoyancy forces and viscous dissipation. As the motivation for this work lies in the connection to pattern-forming processes in salt deserts, we also investigate the length scales and time scales of the resulting dynamics. We also consider the effects of a modulated, non-uniform top boundary condition, which serves as a connection to surface feedback processes in nature.
More broadly, we note that the situation of buoyancy-driven convection occurring below an evaporating salt lake is also closely related to the convective overturning of CO 2 , dissolved in a porous medium filled with brine. This mechanism has become known for its importance to the sequestration of CO 2 in underground aquifers (e.g. Metz et al. 2005;Neufeld et al. 2010;Slim & Ramakrishnan 2010;Slim et al. 2013;Slim 2014;Thomas, Dehaeck & De Wit 2018;Hewitt, Peng & Lister 2020) to help mitigate the anthropogenic impact of CO 2 in the atmosphere. The similarity of models is such that an exchange of methods, insights and results is possible, in both directions. In this context, our study relates to a variation on one-sided convection (Hewitt 2020). For example, analogously to the study by Slim (2014), we numerically investigate the time-dependence of the dynamics of solutal convection in a porous medium, driven from one side, and we find similar regimes of behaviour, ranging from initiation to coarsening and eventually a dynamic steady state.
The dynamics of thermally driven convection in porous media has also been extensively investigated in the past, for a variety of boundary conditions, and the equations are equivalent to the solutal-driven convection discussed here (see e.g. a recent review by Hewitt 2020). Of special interest has been the critical value of the Rayleigh number, Ra c , above which a system is unstable to convective motion, for a variety of situations. For example, a more well-studied case is where convection is driven from two sides, across a domain of fixed height but large width. Here, when two perfectly conducting boundary conditions are chosen, similar to the set-up of Rayleigh-Bénard convection, then transport of heat across the system is purely conductive for Ra < 4π 2 (Horton & Rogers 1945;Lapwood 1948). At higher Ra, first steady and then perturbed convective rolls occur (Busse & Joseph 1972;Graham & Steen 1994). For Ra 1300, the dynamics enters a chaotic regime (Otero et al. 2004;Hewitt, Neufeld & Lister 2012).
There are also certain details already known about the onset of convection in the set of equations and one-sided boundary conditions that we study, which we will briefly summarise. For a constant throughflow boundary condition -a situation that resembles a fluid-filled porous medium with surface evaporation -the onset of instability has been found to be at Ra c ≈ 14.35 (Wooding 1960;Homsy & Sherwood 1976;van Duijn et al. 2002) with a critical wavenumber of k c ≈ 0.76 (Wooding 1960). Finite amplitude perturbations in the range of 8.59 < Ra < 14.35 were also observed to grow in simulations and experiments (Wooding et al. 1997;van Duijn et al. 2002). In the slightly different case of a constant pressure boundary condition at the surface, a situation that resembles an evaporating salt lake covered with brine, the critical values of Ra c ≈ 6.95 and k c ≈ 0.43 are smaller and thus the system is unstable for a greater range of parameters (Wooding 1960). A linear stability analysis (Wooding 1960) as well as an energy minimisation method (van Duijn et al. 2002) were also used to calculate the neutral stability curves for both problems and to determine what range of wavenumbers are unstable at any given Ra.
In what follows we start by reproducing these theoretical results using an approach that develops particularly from the work of Wooding (1960). We are then able to complement the existing analysis by determining the most unstable mode for both constant pressure and constant throughflow boundary conditions, as well as solving for the growth rate of an arbitrary mode. This is followed by a numerical study of the ultimate fate of this instability, including a discussion of how the resulting convection coarsens and scales in the highly nonlinear regime.

Governing equations
The dynamics considered is that of fluid flow in the water-saturated porous soil of a dry salt lake and is described by the mass conservation of water and salt, as well as a detailed momentum balance. The governing equations are those of incompressible flow, the mass conservation of salt with advection and diffusion and Darcy's law in the presence of gravity and are, respectively, These equations describe a fluid of superficial velocity, or flux q and viscosity μ passing through a porous medium of porosity ϕ and permeability κ and carrying with it a dissolved salt of diffusivity D and relative salinity S. Flows are driven by a pressure p and by buoyancy effects due to a gravitational acceleration, −gẑ, and evolve over time t. This system of equations has been studied in a wide variety of contexts, including geysers (Wooding 1960), analogues of Rayleigh-Bénard convection (e.g. Horton & Rogers 1945;Lapwood 1948;Elder 1967;Hewitt et al. 2012), solutal convection (Wooding et al. 1997;Boufadel, Suidan & Venosa 1999) and carbon sequestration applications (e.g. Hewitt, Neufeld & Lister 2014;Loodts, Rongy & De Wit 2014). For example, Hewitt (2020) gives a recent review of their applications to vigorous convection in porous media. Our aim is to determine the main length scales and time scales that emerge from this model of an evaporating salt lake. The connection of the subsurface flow patterns to the surface crust can be seen in the form of (2.2), from which the salinity flux into the crust follows as E − ϕD∂ z S, evaluated at Z = 0. We will return to the interactions of crust and flows in § 5, after we describe how downwelling (low ∂ z S, so higher salinity flux into crust) and upwelling (high ∂ z S and lower salinity flux) structures arise and scale in this system.
The relative salinity S of the pore fluid depends on its density ρ and the boundary conditions. In our system, fluid enters from below (z → −∞) at a background density ρ 0 and evaporates as a saturated solution, of density ρ 1 , from the surface at z = 0. Between these limits where ρ = ρ 1 − ρ 0 . Thus, the salt-saturated fluid in contact with a solid salt crust has a relative salinity of S = 1, whereas the fluid feeding into the soil from some distant reservoir (which will still contain some dissolved salts) has S = 0. The formulation of (2.1)-(2.4) has a number of implicit assumptions that deserve mention. First, it assumes a Boussinesq flow, such that density variations only affect the buoyancy term of Darcy's law. It also neglects the effect of salt concentration on fluid viscosity: Wooding (1960) has discussed this approximation, and shows how a more realistic salinity-dependent viscosity will produce a small stabilising effect; a similar result was obtained by Boufadel et al. (1999), who showed that including a variable viscosity can slightly shift the balance of stability between competing modes of similar wavelengths. Furthermore, the model treats the diffusivity D as a constant, and thus ignores any velocity-induced dispersion (i.e. Taylor dispersion, Taylor 1953). Wooding et al. (1997) show how to consider such effects, but for the field conditions measured at the salt playa of Owens Lake  dispersive effects should be negligibly small. Additionally, the equations treat the porosity and permeability of the soil as constant in space and time. For exemplary research on the types of effects that might be expected in systems with heterogeneous permeability or porosity, see Chen & Meiburg (1998b), Sharp & Shi (2009), Li et al. (2019, Hewitt et al. (2020) and Harfash (2013). Finally, it neglects thermal contributions to density changes, as being small compared with solutal effects. This is justified by the relative magnitudes of these contributions: at Owens Lake, for example, we measured density changes due to salt content to be approximately 200 kg m −3 , whereas a 10 • C day-night temperature change would change the density of water by only approximately 1-2 kg m −3 . The buoyancy ratio, N, which gives the ratio of the solutal to thermal density differences, is thus of order N 100. Effectively, this means that we neglect phenomena like double-diffusive convection, since the driving forces and typical speeds of such flows will be reduced by a factor of 1/N, compared with solute-driven flows (see e.g. Mojtabi & Charrier-Mojtabi (2005) for a detailed review of such effects).
The non-dimensionalisation of our problem requires choosing a characteristic length L and velocity (i.e. flux) V. A natural time scale then follows as T = ϕL/V. Generally, for problems of convection in a porous medium, scaling results in one of two clear choices for dimensionless groups, as (2.2) and (2.3) become where U = q/V, τ = t/T, Z = z/L and where P is an appropriately rescaled pressure. For example, applications in carbon sequestration often use the limiting speed at which a dense parcel of fluid falls, V B = κ ρg/μ, as a natural velocity scale (e.g. Hewitt et al. 2014;Slim 2014;Hewitt 2020). Similarly, in many cases a slab geometry of fixed height is modelled, where a scaling based on the layer thickness can be made (Ruith & Meiburg 2000;Riaz & Meiburg 2003;Hewitt 2020). We will instead follow a scheme where ϕD/LV = 1 and where the group in (2.6) reduces to the Rayleigh number, Ra.
For this, we focus on the situation of a solid salt crust through which there is a constant and uniform evaporation; the more complex case of a modulated evaporation rate will be explored in § 5. As a boundary condition at the surface the upward fluid flux there must balance evaporation, such that the velocity component q z = E at z = 0. This will also give the average vertical flux of the pore fluid anywhere within the soil. Now, the governing equations allow for a simple stationary solution, which represents a dense boundary layer of fluid near the crust. For the one-dimensional problem, involving depth only, this is an attractive solution towards which transients will relax (see e.g. Wooding et al. 1997). Following Wooding (1960), we therefore use L = ϕD/E as the characteristic thickness of the heavy boundary layer which can potentially develop. This natural length scale is also the distance over which advective and diffusive effects will be of comparable magnitude. The natural time scale T = ϕ 2 D/E 2 is then the time fluid takes to cross the boundary layer, in this stationary state, and the characteristic speed V = E. Applying these transformations results in the following non-dimensional formulation of (2.1)-(2.3), Here, Ra = V B /E also represents the ratio of the characteristic speeds of a fluid parcel due to buoyancy and evaporation. The dimensionless salinity flux into the surface, affecting crust growth, is simply 1 − ∂ Z S. Finally, we note that the rescaled pressure, is now defined to include a contribution from the background fluid density. In the rescaled system the boundary conditions are U Z = S = 1 at Z = 0 (where U i are the components of U), with the salinity approaching the limit of S = 0 at large depths. These conditions are analogous to the thermal problem studied by Wooding (1960) and Homsy & Sherwood (1976). They are appropriate to a dry salt lake as long as there is no significant ponding of surface water (Wooding et al. 1997;van Duijn et al. 2002), although we will address the modifications to the model needed for ponding at the end of § 3.1. The stationary solution has U = (0, 0, 1) everywhere, a salinity boundary layer given by S = e Z and a corresponding pressure P = Ra(1 − e Z ) − Z. In the following section we will investigate the stability of this solution, along with a time-dependent solution representing the growth of this boundary layer from an initially homogeneous lake.

Linear stability analysis
Here, we perform a linear stability analysis of our model for small perturbations around an initially unpatterned state. The methods used are inspired by those of Wooding (1960) and extend the range of his results to include an analytic series solution to the linear stability problem, along with predictions of the most unstable mode and first unstable mode. For a base state we will focus on the stationary solution of a well-developed boundary layer, but will also consider the instabilities of a more general time-dependent solution. In § 4 we will confirm these results with a numerical implementation of our model of convection below a dry salt lake and further explore how they are modified at larger amplitudes, in other words in the nonlinear regime of the dynamics.
As set up in § 2, we consider an infinite half-space (Z ≤ 0) of a three-dimensional porous medium saturated with water, which evaporates at the top boundary with a constant evaporation rate E and which is recharged from below by a reservoir of constant-salinity water. To investigate the stability of a horizontally homogeneous base state, we add perturbationsŨ,S,P to its velocity, salinity and pressure fields, respectively. The magnitudes of these perturbations are taken to be proportional to a small parameter, ε 0 . To leading order in ε 0 (specifically, ignoring theŨ · ∇S term, which is of order ε 2 0 ) the perturbations will grow or decay according to Here, the only remaining term related to the base state is the base salinity S 0 = S 0 (Z, τ ), which arises from the nonlinear aspects of the material derivative in (2.9). The pressure term can be eliminated by taking the curl of (3.3), as ∇ × (∇P) = 0. By applying the curl twice, simplifying with (3.1) and considering the Z-component of the result, one finds that (3.3) implies that Solving for the dynamics of the perturbations requires providing for the appropriate boundary conditions. In order to respect the original boundary conditions,Ũ Z =S = 0 at Z = 0. As a second velocity condition, we assume a steady flow far from the unstable surface layer, such that the vertical perturbation to the velocity must decay to zero as Z → −∞. Then, rearranging (3.2) gives Now, we can look at the evolution of the salinity perturbation,S. For this we follow the approach of Pellew & Southwell (1940) as well as Wooding (1960) and assume a separation of variables, such that where Φ is a harmonic function satisfying (∂ 2 X + ∂ 2 Y + k 2 )Φ = 0. Here, k is the characteristic wavenumber of the perturbation in the horizontal directions and α is its growth rate. For α > 0 the amplitude of the perturbation increases and the system is unstable, whereas for α < 0 the perturbation decays and the system is stable; the α = 0 case will give the neutral stability curve. Substituting (3.7) into (3.4) and using (3.5) to simplify the result leads to the following eigenvalue equation for the height-dependent function F(Z), At this point, following the structure of (3.5) and (3.8), we will also introduce (3.9) At Z = 0 the boundary conditions of F = G = 0 follow from the fact thatŨ Z =S = 0 there, and from (3.5). In § 3.1 we will solve this problem for the stationary base state of S 0 = e Z , whereas in § 3.2 we will explore instabilities of the time-dependent case of a boundary layer developing from an initially homogeneous salinity field.
3.1. Instabilities of the stationary base state Using DSolve in Mathematica we obtained an analytical solution of the differential equation (3.8) operating on F(Z), for the stationary base state S 0 = e Z . Wooding (1960) solved a similar problem for the special case of neutral stability (α = 0), whereas we show a more general solution. This solution is potentially a superposition of up to four independent infinite series. Specifically, we first factor (3.8) such that (3.10) where Ψ = 1 + 4(k 2 + α). The solutions then read (3.12) where H i (Z) is defined by (3.14) These are hypergeometric functions of the form r F s with r = 0 and s = 3 (Koekoek & Swarttouw 1998;Askey & Daalhuis 2010). They can be evaluated as series, for example, where (a) n = a(a + 1)(a + 2) · · · (a + n − 1) is the Pochhammer symbol for the rising factorial. Since r < s + 1 these series will always converge, except in the special cases where one of the terms in the denominator is 0 (Koekoek & Swarttouw 1998, p. 12); these exceptions occur when a term in the rising factorial sequence is 0. For example, H 0 is convergent except where 3/2 − Ψ/2 + k = 0, −1, −2 . . ., while H 2 will converge unless the same condition holds for 1/2 + Ψ/2 − k. Hence, the solutions form a dense set in the (α, k) space. From (3.6) it follows that F(Z) has to decay faster than e Z for Z → −∞ (Wooding 1960). Therefore, of the four possible series described above, only those with c i > 1 are allowed, since F i ∝ e Zc i . As k ≥ 0, this condition eliminates the c 1 = 1 − k case, whereas c 0 = 1 + k always leads to a valid solution. Similarly, as long as k 2 + α > 0, then c 3 < 1 and c 2 > 1. Thus, a general solution under these conditions (which include all unstable cases) can be given by F(Z) = C 0 F 0 (Z) + C 2 F 2 (Z), for some real constants C 0 and C 2 .
For our model, as mentioned earlier, the top boundary conditions are satisfied if and only if F(0) = G(0) = 0. Using (3.9) this condition can be written as Since only the non-trivial solution is physically relevant, applying this constraint gives a relationship between Ra, k and α and allows one of these parameters to be determined by fixing the other two. For example, as in figure 3(a), we can use Newton's method to find the roots of the determinant given in (3.16) for some particular values of k and α, and hence find the smallest Ra > 0 satisfying the boundary conditions (e.g. for the neutral stability curve, α = 0). Finally, we note that Wooding (1960) applied slightly different upper boundary conditions to his model, which here would correspond to constant salinity, S = 1, but with a fixed surface pressure instead of a constant surface throughflow. This could simulate a shallow layer of salt-saturated water at the surface, for example, or a water-logged crust. For the fixed pressure case there would be no horizontal flows along the surface (i.e. it would act as a no-slip boundary), and the incompressibility condition (3.1) then implies that ∂ ZŨZ = 0 there. For that scenario the arguments and solutions presented above remain valid, but require the modified boundary conditions of Some results from the solutions to our linear stability problem are given in figure 3 and table 1. In figure 3(a) we show the neutral stability curve and most unstable mode for both types of boundary conditions: results for constant evaporative flux are shown in blue and constant pressure in red. We also indicate the critical points for both cases and details of the critical Rayleigh number, Ra c , and its corresponding critical wavenumber, k c , are given numerically in table 1. For the constant pressure case the neutral stability curve and critical point are consistent with the results of Wooding (1960). For the constant flux case, the corresponding results are consistent with Homsy & Sherwood (1976) and van Duijn et al. (2002). The most unstable mode calculations provide additional predictions, which we will use to validate our numerical model of convection.
The eigenfunctions, F(Z), corresponding to the most unstable modes of scenarios of different Ra are shown in figure 3(b). These have been normalised to have a maximum value of 1. In all cases the general shape of F(Z) is similar to the solution sketched by Wooding (1960) for constant pressure conditions and at the critical point. Furthermore, the eigenmodes do not undergo any significant qualitative changes as Ra increases, but rather the peak gradually narrows and shifts to shallower depths, reflecting the k-dependence of (3.12).
In figure 3(c) we show the growth rate α for various modes and conditions just above the critical point (alongside corresponding results from our numerical model, for purposes of validation). These are consistent with a type-I (finite wavenumber, see Cross & Greenside 2009) instability. In the following section we will confirm these results with a numerical implementation of our model of convection below a dry salt lake and further explore 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 10 2   how they are modified at larger amplitudes, in other words in the nonlinear regime of the dynamics. Finally, and introducing = (Ra − Ra c )/Ra c to describe the proximity of the system to the critical point, figure 3(d) shows that the growth rate of the most unstable mode, α m , scales linearly with near the critical point, as expected for a type-I instability (Cross & Greenside 2009). Given this relationship, when we make quantitative comparisons of e.g. velocities at different Ra in what follows, we will often find it convenient to re-scale time, and thus defineτ = τ .
917 A14-11  3.2. Instabilities of a transitory base state Before proceeding to a full numerical simulation of the salt playa problem, we will consider how instabilities would arise for the case of a transient initial condition. This analysis is inspired by that of Slim & Ramakrishnan (2010), who presented a time-dependent linear stability analysis of the related problem of convection without throughflow, but including the case of a permeable upper surface. Specifically, we will consider the situation where the initial salinity everywhere is equal to that of the reservoir, S = 0, and where the S = 1 boundary condition is suddenly applied to the surface at τ = 0. This could describe the situation of a rapid change in conditions, such as the abrupt flooding of a surface by brine, for example, or the rise of a buried water table to the surface, which could reactivate crust growth.
For this initial condition, the system has the transient solution (Wooding et al. 1997) As shown in figure 4(a), this solution relaxes to the stationary base state, e Z , over a time scale of τ ∼ 1. The question now is whether the instabilities that can occur during such a transient phase will be consistent with those of the stationary base state, or whether they will allow for additional behaviours. We will show that the potential instabilities during the diffusive growth of the boundary layer are broadly similar to the instabilities of the well-developed boundary layer given in § 3.1. Equation (3.8) can be rearranged into an eigenvalue equation of the form AF = αBF, The largest eigenvalue here gives the growth rate of the most unstable perturbation for any particular mode k at any instant τ and Rayleigh number Ra. We solved this eigenvalue problem numerically, using a Chebyshev differentiation matrix (Trefethen 2000) for the differential operators, ∂ Z , acting on F and S 0 , the base state given in (3.17), and on a domain of a finite height h. The lower boundary has a constant flux U Z = 1 of water with relative salinity S = 0, as will be used throughout § 4. The late-time (τ = 100) solutions were used to validate the method against the results shown in figure 3(b,c), and agreed well for domains with a height of at least h 5. Since large h can introduce issues with numerical precision, due to the very rapid decay of (3.17) with depth, we therefore used h = 7.5 in what follows.
For the eigenvalue analysis of the transient base state we used Newton's method to search for the time that the first mode went unstable, for any given Ra. In figure 4(b) we show how the time of the first instability, and the value of the first unstable mode, depend on Ra. For comparison, we also include the most unstable mode of the stationary base state, from figure 3(a). The first unstable mode is similar to the most unstable mode, but consistently slightly higher. The instability sets in rapidly, suggesting that for these initial conditions there will be competition between the growth of the boundary layer, and the growth of the instability. In § 4.4 we will explore this competition further, using fully nonlinear numerical simulations.
To investigate how the spectrum of unstable modes evolves through time, as the boundary layer fills up, we also calculated the growth rate for a grid of different τ and k. Examples are shown in figure 4(c,d) for the cases of Ra = 20 and 40, respectively. We find that the range of unstable modes does not change dramatically over time. Rather, once a mode becomes unstable, it generally remains so. For example, for Ra = 40, modes between k = 0.212 and 3.412 are unstable at τ = 10, and this range is effectively coincident with the range of unstable modes as τ → ∞ (namely, from k = 0.211 to 3.406). The exception to this behaviour is a range of wavenumbers at the highest end of the unstable spectrum. Continuing our example for Ra = 40, the wavenumbers from k = 3.41 to 4.09 are stable to small perturbations at long times, but unstable for some period of the transient. Note that this behaviour is different from the case of the time-dependent diffusion into a finite porous layer without throughflow (for either an impermeable or permeable surface, Slim & Ramakrishnan 2010), where all modes eventually return to a stable situation.
Finally, we note that a more comprehensive analysis of the time-dependent stability problem could be made by non-modal stability theory, as has been done for the related problem of solutal convection without throughflow (see e.g. Rapaka et al. 2008;Slim & Ramakrishnan 2010). Alternatively, we will return to present numerical simulations with a time-dependent base state in § 4.4. What we can conclude here, however, is that the range of unstable modes does not change significantly throughout a transient phase, and that the first unstable mode for the time-dependent case is close to the most unstable mode found in § 3.1, for the stationary base state.

Scaling relationships in the nonlinear regime
For our model of subsurface convection in a dry salt lake the linear stability analysis predicts that a salt-rich boundary layer is unstable to convective rolls above a Rayleigh number of approximately 10, depending on the exact nature of the boundary conditions. Now we focus on the longer-time behaviour of the model and on the scaling of the convective dynamics as the system matures towards a dynamic steady state. To this end, we numerically solved the governing equations, (2.8) to (2.10), on rectangular domains with periodic boundary conditions in the horizontal direction. The deep aquifer is modelled as a lower boundary at a fixed background salt concentration, such that S = 0, and a constant recharge rate, U Z = 1. As an initial condition, we focus on the S = e Z time-independent solution of a well-developed diffusive boundary layer; different initial conditions will be explored in § 4.4. Details of the implementation of the numerical simulation, which follow the pseudo-spectral approach of Riaz & Meiburg (2003), Ruith & Meiburg (2000) and Chen & Meiburg (1998a), are given in appendix A and the code itself is available on a public repository (Lasser & Ernst 2020).
The simulations were first validated against the theoretical predictions of growth rates. Specifically, for Ra between 15 and 40 we added small perturbations of a single wavenumber to the initial conditions. The growth rate of this mode, α fit (k, Ra), was measured by a fit to half the peak-to-peak amplitude of the perturbation in the linear growth regime (see validation section and figure 13 in appendix A). As shown by the crosses in figure 3(c), the simulated α fit agree with the results of the linear stability analysis.
More generally, we added low levels of random noise, at all k, to the initial conditions. We then performed simulations for a selection of Rayleigh numbers ranging from 15 to 2000. The spatial resolution of the simulations increased with increasing Ra, so as to be able to resolve all important features of the dynamics. Consequently, the system sizes were adjusted (smaller width W and height H for higher Ra) to keep the computational cost of the simulations manageable. Spatial resolutions and domain dimensions are listed in Appendix A. Furthermore, in what follows, all uncertainty ranges given represent the standard deviations of properties measured in ensembles of 5 to 10 runs for each set of conditions: they are intended to showcase the variability seen between simulations. Snapshots of an example simulation at Ra = 100 and at different times τ are displayed in figure 5. Supplementary movies S1, S2 and S3 available at https://doi.org/10.1017/jfm. 2021.225, also give the results of three example simulations at Ra = 30, Ra = 100 and Ra = 1000, respectively.
As the simulations proceed they pass through several distinct regimes of dynamics, which can be related to those given by Slim (2014) for a similar problem motivated by one-sided convection beneath a CO 2 pool (in other words, without the evaporative flux of our model). At early times, as shown in figure 5(b,c), we observe a regime of the linear growth of high-salinity plumes, at a wavelength corresponding to the most unstable mode of the linear instability. This is followed by a flux-growth regime, where the downwelling plumes strip the boundary layer of its heavy burden of solute. Such a thinning of the boundary layer can be seen in figure 5(d,e). The next regime is the merging regime, during which time the plumes begin to influence each other via long-range interactions in the horizontal velocity field. As a result, nearby plumes are attracted to each other and merge together to form larger plumes, as is happening in figure 5( f,g). Once enough plumes have merged, the high-salinity boundary layer feeding the plumes begins to grow again, and it thickens until small proto-plumes start emerging at the top boundary and we enter a re-initiation regime, as shown in figure 5(h). These proto-plumes are typically attracted to and then swept into the larger pre-existing plumes. After this time the system settles into  a long dynamic steady state period, which lasts until the deepest plumes start interacting with the lower boundary, after which point we typically stop the simulation. It is worth noting that these simulations start with a well-developed boundary layer, in which upwards advection and diffusion back down the concentration gradient balance. As such, they do not show a clearly defined diffusive regime, which would correspond to the initial growth of this layer, up to the first moment of instability, starting from a homogeneous solute distribution as initial condition (i.e. the dynamics shown in figure 4a). Instead, at early times the boundary layer is deformed by the growth of the unstable modes. To illustrate this, in figure 5(i) we show the horizontally averaged salinity distributions for the simulation snapshots displayed in panels (a-d). A diffusive regime, analogous to the first regime reported by Slim (2014), would be expected for a homogeneous initial condition. Indeed, for such an initial condition, the time of the first unstable mode, reported in figure 4(b) shows how the duration of the diffusive regime would depend on Ra. Similarly, we see a diffusive regime in simulations (see § 4.4) started with a less-developed boundary layer.
We also do not observe a clear shut-down regime, since our boundary conditions allow for a throughflow of solute (rather than a constant build-up of salt that would be seen for impermeable boundary conditions). Indeed, we argue that the variations in salt flux to the surface, caused by the presence of the plumes, are important for the surface patterning seen in dry salt lakes ).

Scaling of the plume velocity
The convective dynamics seen in the simulations tends to become more vigorous at higher Rayleigh numbers. This reflects the interpretation of Ra as the ratio of the natural speeds of flows driven by buoyancy and evaporation. To quantify this relationship we measured the maximum speed of the plumes relative to the background flow by calculating V(τ ) = max(|U Z − 1|) over the entire numerical domain at each time step of various simulations.
In figure 6(a) the time development of V is given for a range of Rayleigh numbers. There are initially small fluctuations in this speed, as the dominant unstable modes are selected from our broad-spectrum perturbation. After this brief initial transient the plume speed, characterising convection, increases until it reaches a peak and then plateaus at an essentially constant value. For Ra > 100 the initial peak can overshoot its plateau value significantly, as can be seen in figure 6(b) for a simulation at Ra = 1000. In all these simulations the initial variations in V correspond to the linear growth, flux-growth and merging regimes whereas the plateau corresponds to the re-initiation regime, which behaves as a dynamic steady state. Again, if we instead used homogeneous initial conditions we would also expect an initial diffusive regime to precede the linear growth phase (i.e. up to the time of first instability given in figure 4b). At very long times the plumes start to interact with the lower boundary and the convection starts to weaken, hence V decreases.
For every simulation we identified the peak speed, V p = max(V(τ )), as shown for example by the blue dot in figure 6(b). We also estimated the beginning of the subsequent dynamic steady state, τ stst , as the time after which ∂V/∂τ < 0.1. This is found to be, empirically, a good measure of the time (e.g. grey dashed line in figure 6b) after which the relative plume speed is no longer systematically changing, although it continues to fluctuate randomly after this point. As shown in figure 6(c), this time to develop a steady state approaches a value of τ stst ≈ 2.5/ (or, alternatively,τ stst ≈ 2.5), for large enough . This can be explained by the time the perturbations need to grow, since the initial amplitude of the perturbations to the salinity field is independent of Ra, whereas the growth rate of the most unstable mode is proportional to . Although the time needed for the disturbances to saturate will depend on the initial amplitude of these perturbations, the scaling of τ stst demonstrates that remains a useful characterisation of the relative vigour of the convective process well into the nonlinear regimes.
To emphasise the scaling of the system times and speeds with , we also calculated a steady state speed, V stst , as the average of V in the time window between τ stst and τ stst + 7 . This window was chosen to be as wide as possible, so as to provide a stable average value, while still avoiding the start of the shut-down regime for the largest Ra in our study, where the shut down of convection also approaches fastest. For example, for the Ra = 1000 case this time window and V stst are indicated as the grey shaded area and the red dashed line in figure 6(b), respectively. In figure 6(d) we show how both V p and V stst vary with , and hence Ra. This shows that the characteristic plume speed in the steady state, i.e. the re-initiation regime, increases linearly with . It also shows how the overshoot of the plume speeds becomes more significant at larger Ra, when the system starts out in a more unstable initial configuration.
Finally, we note that the use of a maximum speed to characterise the rate of convection has the potential to lead to overestimates. As such, we confirmed the scaling of the plume speed by considering the average speed of all downwelling plumes at a depth of Z = −1 and at the moment where the first plume tip reached a depth of Z = −2 (specifically, when S = 0.5 was first exceeded there). The results of this measurement are consistent with those otherwise presented here, in that the average plume speed scales linearly with (Ernst 2017). Similarly computing the average upwelling speed at this time and depth gives the same scaling. Thus, for a number of different metrics, the plume velocity in the system is an indicator of how fast the dynamics passes through the different regimes, and shows  a simple linear scaling with . This can be readily understood, given that Ra = V B /E represents the ratio of the characteristic speeds of buoyancy and evaporation.

Scaling of the plume wavelength
As the buoyancy-driven convection in our model evolves away from its initial instability, the number and spacing of the salt plumes can change. Already, in figure 5, we showed that the pattern has a tendency to coarsen during the merging regime, before a steady state develops where plume merging and re-initiation balance each other. Thus, the most unstable mode calculated by the linear stability analysis is unlikely to be representative of the long-term pattern that would be seen for examples of this convective process occurring in nature. In the following, we will quantify this coarsening behaviour with the aim of determining the spatial scale of convection expected in systems that are in a dynamic steady state, and long after their onset in time.
To characterise the spatial structure of the dynamics in a simulation at any given depth Z and timeτ = τ , we measured the effective wavelength Λ = W/N of the plumes, where W is the width of the simulated domain and N is the number of downwelling plumes of high salinity. For this, the positions of the plumes are identified as the maxima in the salinity along a horizontal profile of depth Z. We then take k = 2π/Λ to be the corresponding wavenumber of the plumes.
In figure 7(a) we show the time development of the plume wavenumber k, as measured at different depths for a simulation run at Ra = 100. The various regimes of the dynamics, described earlier, are reflected in the development of this wavenumber. The growth, merging and re-initiation of plumes are also apparent in the corresponding space-time diagram displayed in figure 7(b), measured at a depth of Z = −1. At first, both these panels show the linear and flux-growth regimes, with closely spaced plumes of large k and small wavelength. Choosing a shallow depth of Z = −1 and early time ofτ = 1 to characterise the initial response, we find that the wavenumber of k = 3.55 ± 0.54 seen at that time is similar to the most unstable mode of the linear stability analysis, k = 3.01.
After the first plumes have more fully formed, at intermediate times the wavenumber of the simulations declines as these plumes begin to strip the boundary layer of salinity and then merge into larger structures. For the Ra = 100 simulation shown in figure 7(a,b) this happens between approximatelyτ = 1 and 10. Fromτ ≈ 10 onward the plume wavenumber measured deeper into the simulated domain approaches a stable value of k ≈ 1 (see blue dashed curve in figure 7a for the example of Z = −10). The wavenumber closer to the top boundary (e.g. red dashed line in figure 7a, for Z = −1), fluctuates between that value and a higher value of k ≈ 3. This is indicative of the episodic re-initiation of proto-plumes that can also be seen by the herringbone pattern in the space-time diagram of figure 7(b). These proto-plumes are usually ephemeral, and merge into a larger plume before they can reach into, and be noticed at, the larger depths. A similar process of intermittent plume initiation was seen in the related model (i.e. without evaporation) studied by Slim (2014), who referred to them as proto-plume pulses.
In order to quantify the coarsening of the plume wavelength more generally, we measured Λ and k for a range of conditions, with results given in figure 7(c) for simulations with Rayleigh numbers between 14.9 and 1900. For a limiting case that is near the onset of the initial instability, we looked at the early (τ = 1) and near-surface (Z = −1) response in all simulations. As shown in figure 7(c) by the red filled circles, these values closely follow the theoretical prediction for the most unstable mode (grey line, from figure 3a), although they tend to slightly exceed it. As the system ages, however, the dependence of the plume spacing on Ra weakens. We tracked this behaviour by measuring k at progressively later times and lower depths (where values are less volatile), as shown by the data in figure 7(c), with the various measurement depths indicated by the figure inset. These results demonstrate how, for a wide range of Ra, the long-time limit of the plume wavelength appears to gradually approach the value that would occur at the critical point, namely k c 0.76. This is highlighted by the blue circles in figure 7(c), which show the wavenumbers as measured atτ = 30 and a depth of Z = −10. We will argue in the following section that the relative independence of k on Ra results from the depletion of the salt-rich upper boundary layer by plume formation, which effectively reduces it to the thickness of a system forced just beyond its critical point.

Dynamics of the high-salinity boundary layer
The emergence of plumes in our model is driven by the negative buoyancy of the salt in the diffusive boundary layer near the surface of the dry salt lake. Therefore a closer look  at the dynamics of the effective thickness of this layer is warranted, along with its link to how the mature plume wavelength emerges. We have already shown several instances of where plumes drain the boundary layer of solute, thereby limiting the convective drive: compare the salinity distributions near the upper boundaries in figure 5(a,h), for example. This trend can also be seen in figure 8(a),  where we show the horizontally averaged salinity distributions, S(Z) , as measured in the dynamic steady state regime of simulations at different Ra. Near the surface these results all demonstrate a rapid decay of the salinity with depth, which is stronger for higher Ra. Figure 8(a) also shows how, just below this boundary layer, and especially for higher Ra, the salinity may also pass through a small local maximum of around S ≈ 0.3. A similar salinity peak can be noticed in Slim (2014) (figure 3). At intermediate depths, below Z ≈ −2, the salinity then either approaches a constant value or gradually trails off. We note, however, that some of the apparent difference in internal structure seen at these lower depths may simply be due to the restricted height of the simulations at higher Ra (see appendix A). Similarly, for very high Ra the salinity starts building up below Z ≈ −9, since the bottom boundary of the simulated domain is only at Z = −10 in these cases. Given the shape of the horizontally averaged salinity distribution seen in both the initial and mature states of our simulations, we estimated the effective boundary layer thickness at various times and Ra. To this end, we fit an exponential decay function S(Z) = S exp(Z/L ) + S bulk to the rapidly decaying part of the salinity distribution found just below the top boundary, as demonstrated in the inset to figure 8(a). For reference, in the stationary solution to (2.8) to (2.10) the relative salinity S = exp(Z), so for our initial conditions L = 1. The results for the evolution of L with time are shown in figure 8(b).
As with the other metrics discussed here, the width of the salinity boundary layer evolves in different ways as the simulation passes through its various regimes. Initially, in the linear growth phase the perturbations are not large enough to affect the salinity field, and L 1 for short times (inset to figure 8b). For higher Ra the growth rate of these perturbations is faster, so this initial regime is shorter (again, this speeding up is well captured by the use ofτ for time). As the instability enters the flux growth regime, the boundary layer then shrinks as solute is removed by the growing plumes. Eventually this process stabilises and, as the plumes merge into fewer but larger structures, reverses itself while the boundary layer between the plumes is recharged (we note that this transition can also be related to the velocity overshoot shown in figure 6b). Then, when L becomes sufficiently large it allows for the re-initiation of proto-plumes, which repeat the cycle of growing in amplitude, depleting the boundary layer and merging into larger plumes. Thus, in the dynamic steady state of the simulations the boundary layer thickness fluctuates around a value of L < 1, corresponding to the episodic emergence of proto-plumes. The frequency of these proto-plume pulses, and the average value of L , depend on the Ra of the system: for higher Ra the boundary layer is thinner and the fluctuations are more rapid.
For convection driven by one side, Slim et al. (2013) suggested that (in our notation) an effective Rayleigh number of Ra(1 − S bulk )L would describe a boundary layer that had been disturbed by, for example, loss of material to downwelling plumes. Essentially, (1 − S bulk )L gives the relative buoyancy forces available in the boundary layer, whereas Ra characterises the system's general ability to act on those forces. In all our simulations, the thickness of the effective boundary layer, L , appears to scale approximately with 1/Ra in the steady state regime. More specifically, as shown in figure 8(c), the ratio Ra(1 − S bulk )L /Ra c approaches a constant value of about 1.8 for Ra above about 100. As further evidence of this scaling, figure 8(d) shows a data collapse of the near-surface salinity distributions, S(Z) , as rescaled by Ra/Ra c . These results all suggest that the proto-plume pulses continuously trim back the boundary layer and maintain it in a state that is close to, but just above, a critical condition. This is consistent with the tendency of the mature plume wavenumber to stabilise at a value near what is expected at these conditions, namely k c .

Effects of varying the initial conditions
Finally, we look at how different initial conditions will modify the system's approach to a mature convection pattern and demonstrate that the dynamic steady state discussed above is robust. Up to this point our simulations have begun with an initial condition that can be defined as S = exp(Z/L 0 ), where L 0 = 1, and with random perturbations characterised by an amplitude of η = 0.05 (see appendix A, (A6)). Now, we varied both the depth of the boundary layer used as the initial salinity distribution and the perturbation amplitude. For this parameter study we used a somewhat restricted domain, extending down to only Z = −20. The exemplary case of L 0 = 1 is shown in figure 9(a), which has similar features to those discussed at length in regards to figure 8. We note that at late times here there is an additional slight upward drift of the salinity throughout the domain, which can be explained by a gradual saturation of the salinity in the system once the plumes start interacting with the lower boundary of the simulation.  In figure 9(b) we show the development of the effective boundary layer L for a selection of initial conditions L 0 . Similarly, in figure 9(c) we show the effect of varying the initial level of noise in our simulation, for a case where L 0 = 0.4. When L 0 < 1 the boundary layer first grows towards the stationary solution of L = 1 and then shrinks again, as the convective instability sets in. Additionally, for lower levels of noise the instability takes longer to manifest itself, allowing the boundary layer slightly more time to be established. As such, the maximum value of L depends on the initial conditions, in line with the competition between the growth of the boundary layer and the growth of the most unstable mode. These processes are, respectively, a roughly exponential relaxation of the boundary layer thickness (consistent with the fact that the stationary solution is stable to long-wavelength perturbations, see figure 3) and an exponential growth of the first plumes. Thus, the cross-over time between these two processes is relatively insensitive to the initial conditions, and occurs aroundτ ≈ 1.
For all these cases, the simulations converge towards the same salinity distributions above times ofτ ≈ 2. This is emphasised in figure 9(d), which gives the late-time behaviour of L for the simulations with different initial boundary layer thicknesses.
Here, all simulations show the same value of L in the dynamic steady state (although, the smaller domain size adds a slight saturation of S as compared with figure 8). As described earlier, in this regime the boundary layer thickness self-organises into a state that balances the initiation and growth of new plumes with their coarsening and merging into larger structures.

Effects of spatially varying evaporation rates
Our model of convection in dry salt lakes is inspired by the polygonal patterns that are often seen in the salt crusts at their surface and we will end this work with a discussion of the possible interaction of this crust with the convective plumes. We first recall that the model also predicts the salinity flux into such a surface crust, which follows from (2.2). In dimensional terms, it can be given by E − ϕD∂ z S. This flux vanishes in the stationary solution corresponding to our initial condition, S = exp(ZE/ϕD), although this scenario will still leave a constant upward flux of salt into the crust at whatever concentration is supplied by the reservoir (since S = 0 corresponds to the background density ρ 0 , see (2.4)). By similar argument, when there are convective plumes there should be less salt flux into the crust above an upwelling -where ∂ z S is high -than above a downwelling.
A fully dynamical crust, with a thickness varying in response to the salinity flux from the convection beneath it, is challenging to model. In particular, realistically simulating the evaporation rate can be difficult (see § 5.1 below), but an aspect of this problem that we intend to address in the future. In this section we will briefly consider, instead, the other side of how a feedback between crust patterns and convection patterns could work. In other words, we will look at whether periodic variations in the properties of a salt crust could influence any convection pattern occurring beneath it. Since evaporation is the ultimate driver of the dynamics, this allows us to consider the effect of a salt ridge through its control over the local evaporation rate at its location. In particular, we ask whether particular wavelengths of surface features could stabilise plume locations, leading to potential for long-term feedback between subsurface flows and the crust pattern.
To this end, we modified our model by modulating the surface evaporation rate, such that U Z = 1 − A m cos(k m X) at Z = 0, and where A m and k m are the amplitude and wavenumber of this modulation, respectively (see appendix A for details of how this is implemented). As the average evaporation rate remains unchanged by the modulation, this change in boundary conditions does not affect the system-averaged Rayleigh number, merely the local conditions. Space-time diagrams for simulations where Ra = 100, A m = 1 and for four different modulation wavenumbers k m are shown in figure 10(a) and movies of the corresponding simulations are given in supplementary movies S4 through S7. Several things can be noticed in these simulations. First, in all cases there is a preference for downwelling plumes to originate at the minima of the evaporation profile. Second, in most cases the plumes then remain locked at the positions of these minima for some time, before the order breaks down and the dynamics returns to a state which resembles a system with a uniform boundary condition. Third, however, when the wavelength of the evaporation modulation is close to k c ≈ 0.76 there is a tendency for the larger plumes to remain trapped near the spots of lower evaporation, even in the dynamic steady state.
To quantify these points further, we introduce an order parameter that characterises how co-aligned the plumes and evaporation patterns are. For N plumes this is given by where X i is the horizontal position of plume i, as determined by a minimum in the salinity at a depth of Z = −1. If ξ = 1 then the plumes are perfectly aligned with minima in the surface evaporation, whereas when ξ = −1 the plumes all lie below maxima and if there is no preferred arrangement of plumes, then ξ = 0. A large ξ would then suggest a route for selecting preferred wavelengths in the crust features, with downwelling plumes trapped by lower-evaporation ridges, and sustaining them with enhanced salinity flux. In figure 10(b), we show the development of ξ with time for simulations with the same parameters as in figure 10(a). This plot confirms that, especially for larger k m , the initial plumes start very well aligned with the surface modulation, with ξ up to about 0.9 at early times. Since Ra = V B /E, lower-evaporation regions will appear, locally, as having a higher effective Ra, and it makes sense to expect a higher growth rate of plumes there. The figure also shows that the duration of the initial pinning of plumes under the evaporative minima depends on the modulation wavenumber k m ; i.e. the surface modulation can delay the transition from the flux growth to the merging regime. Furthermore, in all cases the long-term limit continues to show at least weak ordering, with ξ fluctuating around positive values of about 0.1-0.2 or greater. The k m = 0.63 case shows a marked contrast, however. Here, the surface modulation is well matched to the wavelength seen in the dynamic steady state of the simulations and the rise of ξ to about 0.6 documents how initially disordered plumes arrange themselves to synchronise with the spacing and phase of the modulation.  Figure 11 shows how the above results hold true for a wide range of k m and A m . Here, the choice of modulation wavenumbers is limited such that the system width is a multiple of the modulation wavelength. Within this constraint, we chose to investigate a range of modulation wavenumbers in an area of the parameter space that promised to exhibit interesting behaviour, including the full range of wavenumbers seen in the dynamical simulations of Ra = 100 for homogeneous boundary conditions (see figure 7c). For the early-time ordering we define a phase-locking time,τ φ , as the time it takes for the order parameter to drop below ξ = 0.75. Note that this is an arbitrary value, but results are similar for other parameter choices of how to characterise the early-time synchronisation of plumes with flux patterns. Figure 11(a) shows how this time depends on the modulation details. The initial phase locking is strongest, i.e. the plumes stay aligned with the minima in the surface modulation the longest, for modulation wavenumbers between k m = 1.6 and k m = 3.1. (n.b. a rapid period-doubling instability can be seen in figure 10(a) for k m = 1.26, suggesting how this alignment breaks down at low k m ). As might be expected, the modulation wavenumbers for which the initial plumes stay aligned the longest roughly corresponds to the wavenumbers first seen to be unstable in the simulations with homogeneous evaporation rates: for Ra = 100 this is k ≈ 3.1.
To measure the long-time synchronisation of the plumes with the surface modulation we instead averaged the order parameter ξ fromτ = 30 to 60. Figure 11(b) shows how this average order parameter, ξ , depends on the modulation details. Similar to the phase-locking time, the strength of the ordering increases with the modulation amplitude A m . More interestingly, ξ exhibits a pronounced maximum in the wavenumber range of k m = 0.3 to 0.9, which is already present for moderate amplitudes of A m = 0.3. This range broadly matches the critical wavenumber k c ≈ 0.76 as well as the wavenumbers of polygonal salt crust patterns observed in nature, where k crust = 0.78 ± 0.43 with no observable dependence on Ra .

Evaporation rate modulation in nature
In the following we will briefly show how a modulation of the evaporation rate could be realised in the setting of a salt desert with salt ridges in a surface crust. Although the effect of a salt crust on evaporation can be hard to evaluate exactly (see e.g. Eloukabi et al. 2013;Bergstad et al. 2017;Farhat 2018;Nachshon et al. 2018), the temperature and relative humidity are important parameters, next to salt concentration and air movement, controlling evaporation from a salt pan. We attempted to estimate the influence of ridges on the micro-climate at the crust level by embedding sensors within salt polygons found in Owens Lake (CA) and tracking the temperature and relative humidity between 26 November and 2 December 2016. We used HiTemp140 and RHTemp1000IS data loggers, which recorded temperatures and relative humidity every two minutes with a precision of ±0.01 • C and ±0.1 %, respectively. The resulting data are deposited on a public repository  and the protocol for their collection is described in more detail by , along with further details of the field site.
In figure 12 we show temperature and relative humidity measurements from sensors placed inside a salt ridge and within the crust at the centre of a polygon. In figure 12(a) the diurnal fluctuations of temperature between about 10 • C at 2:00 PM and 0 • C at 6:00 AM are clearly visible. We note that these temperature changes are unlikely to directly affect the density-driven flows by thermal expansion, as they would induce a density change of no more than 1 kg m −3 . For comparison, we measured the water immediately below the crust to have a density of at least 200 kg m −3 higher than that at depths of approximately 1 m . This difference in magnitude is what justified our original assumption (see § 2) of ignoring thermal contributions to fluid density, and double-diffusive effects. These periodic daily fluctuations are also fast compared with the growth of the crust, which occurs over weeks to months (Nield et al. 2015). Additionally, the highest temperatures recorded during the day are similar between the centre of the crust polygon and the ridge.
Temperatures inside the ridge, however, drop faster and to lower values at night; the temperature difference is approximately 2 • C on average. Figure 12(b) shows the development of the relative humidity over the same period. The relative humidity measured inside the ridge is approximately equal (on the first, third and sixth days) or up to 15 % higher (on the second, fourth and fifth days) to that in the polygon centre. Humidity differences are most pronounced during nights and mornings but are preserved to some extent over the course of the day. The difference in relative humidity can be explained by the trapping of moist air below the ridges (see Nachshon et al. 2018 for a discussion of how trapped, stagnant air can also reduce evaporation through salt crusts). Both reduced temperature and increased relative humidity inhibit evaporation from the surface of a salt ridge and could serve as part of the feedback mechanism proposed above, which modulates the evaporation rate.
In a similar geographic setting a temperature reduction of 1 • C and relative humidity increase of approximately 20 % resulted in a decrease of the evaporation rate from 3 mm day −1 to 2 mm day −1 , or approximately 30 % (Farhat 2018). These differences are comparable to the temperature and relative humidity differences we measured between a salt ridge and polygon centre. In figure 11(b) we have illustrated that modulation amplitudes of A m ≥ 0.3 are sufficient to cause a significant spatial ordering of the downwelling plumes. Our observations therefore suggest that the modulation of evaporation rates that could reasonably occur in a field setting would be sufficient to influence the plume positions.

Summary and discussion
We have presented a linear stability analysis and subsequent numerical study of buoyancy-driven convection in a fluid-saturated porous medium with surface evaporation and where fluid is replenished from a distant reservoir. This model is inspired by the situation below a dry lake or salt pan and by the possible connection of the convective dynamics below the ground to the emergence of regular polygonal patterns at the surface (Lasser 2019;Lasser et al. 2019). When rescaled by the evaporation rate E and a length scale that balances advection and diffusion, L = ϕD/E, the model is controlled by a single dimensionless group, the Rayleigh number Ra. In this context, Ra can be interpreted as the speed at which a large blob of salt-rich fluid would naturally descend, V B , relative to the upward flux of fluid required to balance evaporation.
There is a stationary solution for the resulting system of equations, corresponding to a salt-rich boundary layer of fluid, of thickness L, lying just below the evaporating surface. Our linear stability analysis considers whether this solution is stable or not, and complements earlier work concerning the onset of convection in equivalent models (Wooding 1960;Homsy & Sherwood 1976;Wooding et al. 1997;van Duijn et al. 2002). In particular, we confirmed the critical conditions and neutral stability curves given in those works for the surface boundary conditions of either constant evaporation or constant fluid pressure. Additionally, our analysis extends on previous approaches by solving for the growth rate of an arbitrary small-amplitude perturbation at any Rayleigh number. Through these methods the initial growth of convective plumes near the evaporating surface was shown to generally be a type-I/finite-wavelength instability, where the most unstable mode increases with increasing Ra.
In order to follow the evolution of the convective instability past its initial stages, we then performed a range of numerical simulations. The dynamics of the convection in these simulations passes through several regimes: the sequence is similar to one identified by Slim (2014) for the related case of porous-medium convection also driven from one side, but without any evaporation or throughflow. The plumes initially follow a linear growth regime, where their wavelength closely matches the predictions of the linear stability analysis. As their amplitude grows, however, they begin to deplete the salt-rich fluid near the surface, and the dynamics passes into a flux-growth regime. In the subsequent merging regime the system coarsens as the growing plumes begin to interact and join together into larger structures. This process leaves gaps between plumes and in the re-initiation regime new proto-plumes appear in these gaps, via similar instabilities in the boundary layer. The emerging proto-plumes are quickly drawn into larger and more stable downwelling plumes, allowing for new proto-plume pulses to occur episodically. We characterised this re-initiation regime as a dynamic steady state of the system.
Throughout this study we focused on exploring how the dynamical properties of the convection scaled with its driving parameters, here summarised by the Rayleigh number or by the related proximity of the system to its critical point, = (Ra − Ra c )/Ra c . In the linear growth regime the growth rate of the most unstable mode α m ∼ , as expected for a type-I instability. Similarly, for later regimes the plume velocities, measured in various ways, were found to be proportional to . This follows naturally from the definition of Ra = V B /E, since the velocity scaling of our system is based on the evaporation rate E. Given this, we found a time scaleτ = τ = tE 2 /ϕ 2 D to be a convenient rescaling. Specifically, the various regime transitions occur at times of orderτ = 1 for a wide range of initial conditions and Rayleigh numbers. In terms of length scales, the most unstable mode of the linear instability increases monotonically with Ra. However, due to the tendency of plumes to merge, this response is only seen at very short times -we would not expect it to affect the growth of salt crusts in realistic settings. Instead, in the dynamic steady state we found that the constant interplay between proto-plume initiation and merging resulted in a plume spacing that was largely independent of Ra, such that the wavenumber of the mature plumes always approached the critical wavenumber of k c ≈ 0.76.
We note that this last result is in contrast to the more well-studied case of two-sided convection, where the thickness of the entire convecting domain, H, imposes a scaling of the mature wavenumber with a Rayleigh number defined alternatively as Ra = V B H/ϕD (e.g. Hewitt et al. (2014) who suggest k ∼ √ Ra or Fu, Cueto-Felgueroso & Juanes (2014) who suggest k ∼ Ra). In a two-sided system, if a plume naturally moves at a speed ∼ V B across the domain height H, then a plume spacing of order 1/ √ Ra reflects the distance over which concentration gradients would diffuse over the course of the plume's fall (Liang et al. 2018). For our one-sided case of convection the system height is irrelevant. Instead, the boundary conditions provide the length scale, L, over which the typical contributions of advection and diffusion balance, and our Rayleigh number can be written as Ra = V B L/ϕD. The length L is also the equilibrium thickness of a heavy boundary layer of fluid that would otherwise naturally develop below the surface and so it characterises the potential driving force available for convection. We argued that the relative independence of the plume spacing with Ra, observed in out simulations, relates to how convection depletes this salt-rich boundary to leave a layer just barely thick enough to allow for convection. This conclusion was supported by a demonstration that in the steady state regime the effective thickness of the boundary layer, L , is proportional to 1/Ra over a wide range of Ra and by a data collapse of the shape of the horizontally averaged salinity field with the same scaling. Put simply, if the boundary layer grows much thicker than this, then it will favour the rapid growth of new proto-plumes, which will strip the layer back down to close to a critical thickness before they disappear through plume mergers. We argued that this balance also controls the plume spacing and is why the steady state wavenumber of the convection plumes is always maintained at a value near k c , regardless of the real Rayleigh number of the system.
Finally, we modified our model to allow for inhomogeneous boundary conditions, namely a sinusoidal modulation of the evaporation rate in space. This modulation is a first step for exploring how feedback between subsurface convection and surface crust growth could work. The influence of ridges on the evaporation rate is supported by data measured in the field, which shows a difference in temperature and relative humidity below ridges, as compared with salt polygon centres. We found that for a range of modulation wavenumbers and amplitudes, regions of locally suppressed evaporation could pin downwelling plumes in place for long periods of time. This pinning was particularly apparent for modulations around the critical wavenumber, k c .
Put together, our results indicate that the vigorous convection patterns that emerge in our model system are robust to differences in Ra, for Rayleigh numbers far enough above Ra c . Such a robustness to fluctuations in the environment is an important feature of any mechanism driving salt polygon emergence in nature, as these patterns occur in areas with vastly different conditions but nonetheless display remarkably consistent length scales (e.g. Christiansen 1963;Krinsley 1970;Nield et al. 2015). A plume wavelength of 2π/k c corresponds to features with a spacing of a few metres, assuming diffusion constants of order 10 −9 m s −2 and evaporation rates of order 1 mm day −1 (or, 10 −8 m s −1 ). This agrees both with the observed sizes of salt polygons and with the scales of convective plumes seen under similar conditions in tidal flats and sabkhas (Stevens et al. 2009;Van Dam et al. 2009). We intend to develop this argument further with detailed comparisons with field data elsewhere . Furthermore, although motivated by the problem of a dry salt lake, the model system developed here could also be applied to other cases of convection in a porous medium, where there is some background throughflow of fluid across the convecting domain. evaporation rate, U Z = 1 − A m cos (k m X) there, where A m represents the strength of the modulation and k m its wavenumber; for constant evaporation, A m = 0. In all cases the average value of U Z at the surface is 1. At the bottom boundary, Z = −H, fluid recharge is assumed to be uniform, such that U Z = 1 there.
We use potential functions to reformulate the equations of porous-medium flow, following Ruith & Meiburg (2000) and Riaz & Meiburg (2003). For a two-dimensional and incompressible flow the flux U is related to the Lagrange streamfunction ψ by U = (∂ Z ψ, −∂ X ψ). The vorticity, ω = ∂ X U Z − ∂ Z U X , is then given by the Poisson equation, Taking the curl of Darcy's law, (2.10), allows us to eliminate the pressure term, as ∇ × (∇P) = 0, and solve for the vorticity as The salt mass balance of (2.9) can also be written in terms of the streamfunction To solve (A1)-(A3) we need boundary conditions for the streamfunction ψ and vorticity ω. The constraints on flow give ∂ X ψ = −1 + A m cos(k m X) at Z = 0 and ∂ X ψ = −1 at Z = −H. These inhomogeneous boundary conditions can be accounted for by taking ψ = ψ + Ψ , where such that ∂ X Ψ = 0 at both Z = 0 and Z = −H. The constant salinity conditions correspond to a vanishing vorticity, ω = 0, at both these surfaces. As an initial condition we make use of which is the stationary solution of the salinity when A m = 0. Following Riaz & Meiburg (2003) we introduce perturbations by adding random fluctuations into the initial salinity.
To do so we generate random numbers f (X, Z), uniformly distributed in [−1, 1], for all grid points. In order to avoid artefacts in derivatives, these are then convolved with a Gaussian function of width σ = 3 grid cells, such that where η gives the magnitude of the perturbation. The initial salinity field is then S = S 0 + S , at τ = 0. Unless otherwise stated we used a default value of η = 0.05.
A.2. Implementation To numerically solve the governing equations of the simulation we largely follow the implementation described by Riaz & Meiburg (2003) and Ruith & Meiburg (2000). Specifically, at each time step we: (i) compute derivatives in the X-direction by first making use of a Fourier transform; (ii) compute derivatives in the Z-direction by using a compact finite difference scheme (Lele 1992); and (iii) use an explicit fourth-order Runge-Kutta scheme for the time integration of equation (A3).  The simulation is performed on a grid of M × N points at positions (x m , z n ). The grid spacings X and Z are varied with Ra, since the spatial resolution has to be increased at higher Ra to resolve all relevant features. Domain sizes were also adjusted, allowing for a similar number of grid points in most simulations. The simulation grid spacings and domain width W and height H are given in table 2 for all Ra. The time step τ follows the Courant-Friedrichs-Lewy (CFL) condition (Courant, Friedrichs & Lewy 1928), where U max is the maximum speed occurring at that time, and where we set C 0 = 0.1. Derivatives in X-direction: we employ Fourier expansions, with coefficientŝ and analogous expressions forψ ξ andω ξ , where ξ can take integer values between ±M/2. In terms of these Fourier coefficients (A1) may be written as α 0 α 1 a 1 a 2 a 3 a 4 a 5 a 6 interior -1/3 1 4 /9 1 /9 ---node 0 5 -197/60 −5/12 5 −5/3 5/12 −1/20 node 1 3/4 1 /8 −43/96 −5/6 9/8 1 /6 −1/96 - interior 12/97 −1/194 120/7 ----boundary 11/12 −131/4 177/16 −507/8 783/8 −201/4 81/16 −3/8 Table 3. Coefficients for sixth-order compact finite difference schemes used in the numerical simulations, following Lele (1992) and Tyler (2007). The α and a terms are used to calculate first derivatives in Z, while β and b terms are used for second derivatives. whereω ξ are calculated at each time step by a similar transformation of (A2), and wherê ω ξ 0 = ∇ 2ψ ξ are constant in time and arise from the boundary conditions. Derivatives in Z-direction: to solve the system of equations given by (A9) we computed ∂ 2 ZΨ ξ (z n ) following the implicit sixth-order compact finite difference scheme given by Carpenter, Gottlieb & Abarbanel (1993). The N linear differential equations can be described by two matrices such that To construct the matrices A left and A right , we use the coefficients α i , β i , a i and b i listed in table 3, as calculated by Tyler (2007). At each time step we then solve for the Fourier componentsψ, which are inverted to give the streamfunction ψ(x m , z n ). The flux U is then calculated from the streamfunction by again using an implicit sixth-order compact finite difference scheme for the spatial derivatives. Coefficients for the first and second-order derivatives of interior and boundary points are listed in table 3. Time integration: to update the salinity field via (A3) we used a fourth-order Runge-Kutta scheme and the coefficients given for example by Süli (2003, p. 352). For spatial derivatives we again use an implicit sixth-order compact finite difference scheme.

A.3. Validation
The model was validated by comparison with the results of the linear stability analysis. A similar approach is used in Ruith & Meiburg (2000), Riaz & Meiburg (2003), Chen & Meiburg (1998a) and Tan & Homsy (1988). For this, instead of (A6) we perturbed our initial conditions with using a perturbation of initial magnitude η = 10 −6 . This perturbation is consistent with the constant-salinity boundary condition and affects only a single wavenumber k. Its growth over time was measured from (half) the peak-to-peak amplitude of this mode evaluated at a depth of Z = −1. The simulated growth rate α fit was determined by fitting an exponential to the amplitude measurements. For this fit we focus on the linear growth phase, manually excluding any initial transient or later nonlinear saturation. Some example measurements, and fits, are shown in figure 13. Results for various Ra and k agree with the theoretical values to within a relative error of order 1 %, as was shown in figure 3(b).