2. 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,

(2.1)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{q} = 0, \end{gather}

(2.2)\begin{gather}\varphi \frac{\partial S}{\partial t} + \boldsymbol{q} \boldsymbol{\cdot} \boldsymbol{\nabla} S = \varphi D\nabla^2 S, \end{gather}

(2.3)\begin{gather}\boldsymbol{q} ={-}\dfrac{\kappa}{\mu} (\boldsymbol{\nabla} p + \rho g \hat{\boldsymbol{z}}).\end{gather}

These equations describe a fluid of superficial velocity, or flux $\boldsymbol {q}$ and viscosity $\mu$ passing through a porous medium of porosity $\varphi$ and permeability $\kappa$ 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\hat {\boldsymbol {z}}$, and evolve over time $t$. This system of equations has been studied in a wide variety of contexts, including geysers (Wooding Reference Wooding1960), analogues of Rayleigh–Bénard convection (e.g. Horton & Rogers Reference Horton and Rogers1945; Lapwood Reference Lapwood1948; Elder Reference Elder1967; Hewitt et al. Reference Hewitt, Neufeld and Lister2012), solutal convection (Wooding et al. Reference Wooding, Tyler and White1997; Boufadel, Suidan & Venosa Reference Boufadel, Suidan and Venosa1999) and carbon sequestration applications (e.g. Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2014; Loodts, Rongy & De Wit Reference Loodts, Rongy and De Wit2014). For example, Hewitt (Reference Hewitt2020) 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 -\varphi D\partial _zS$, evaluated at $Z=0$. We will return to the interactions of crust and flows in § 5, after we describe how downwelling (low $\partial _zS$, so higher salinity flux into crust) and upwelling (high $\partial _zS$ and lower salinity flux) structures arise and scale in this system.

The relative salinity $S$ of the pore fluid depends on its density $\rho$ and the boundary conditions. In our system, fluid enters from below ($z\rightarrow -\infty$) at a background density $\rho _0$ and evaporates as a saturated solution, of density $\rho _1$, from the surface at $z = 0$. Between these limits

(2.4)\begin{equation} S = \frac{\rho - \rho_0}{\rho_1 - \rho_0} = \frac{\rho - \rho_0}{\Delta\rho}, \end{equation}

where $\Delta \rho = \rho _1 - \rho _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 (Reference Wooding1960) 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. (Reference Boufadel, Suidan and Venosa1999), 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 Reference Taylor1953). Wooding et al. (Reference Wooding, Tyler and White1997) show how to consider such effects, but for the field conditions measured at the salt playa of Owens Lake (Lasser et al. Reference Lasser, Nield, Ernst, Karius, Wiggs and Goehring2019, Reference Lasser, Nield and Goehring2020) 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 (Reference Chen and Meiburg1998b), Sharp & Shi (Reference Sharp and Shi2009), Li et al. (Reference Li, Cai, Li, Li and Chen2019), Hewitt et al. (Reference Hewitt, Peng and Lister2020) and Harfash (Reference Harfash2013). 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}$ (Lasser & Goehring Reference Lasser and Goehring2020; Lasser et al. Reference Lasser, Nield and Goehring2020), whereas a 10 $^\circ$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\simeq 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 (Reference Mojtabi and Charrier-Mojtabi2005) for a detailed review of such effects).

The non-dimensionalisation of our problem requires choosing a characteristic length $L$ and velocity (i.e. flux) $\mathcal {V}$. A natural time scale then follows as $T = \varphi L/\mathcal {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

(2.5)\begin{gather} \frac{\partial S}{\partial \tau} + \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} S = \frac{\varphi D}{L\mathcal{V}}\nabla^2 S, \end{gather}

(2.6)\begin{gather}\boldsymbol{U} ={-}\boldsymbol{\nabla} P - \frac{\kappa \Delta \rho g}{\mu \mathcal{V}} S \hat{\boldsymbol{Z}}, \end{gather}

where $\boldsymbol {U} = \boldsymbol {q}/\mathcal {V}$, $\tau = 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, $\mathcal {V}_{B} = \kappa \Delta \rho g / \mu$, as a natural velocity scale (e.g. Hewitt et al. Reference Hewitt, Neufeld and Lister2014; Slim Reference Slim2014; Hewitt Reference Hewitt2020). 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 Reference Ruith and Meiburg2000; Riaz & Meiburg Reference Riaz and Meiburg2003; Hewitt Reference Hewitt2020). We will instead follow a scheme where $\varphi D/L\mathcal {V} = 1$ and where the group in (2.6) reduces to the Rayleigh number, $\textit {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,

(2.7)\begin{equation} S(z) = \textrm{e}^{zE/\varphi D}, \end{equation}

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. Reference Wooding, Tyler and White1997). Following Wooding (Reference Wooding1960), we therefore use $L = \varphi 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 = \varphi ^2 D/E^2$ is then the time fluid takes to cross the boundary layer, in this stationary state, and the characteristic speed $\mathcal {V} = E$. Applying these transformations results in the following non-dimensional formulation of (2.1)–(2.3),

(2.8)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0, \end{gather}

(2.9)\begin{gather}\frac{\partial S}{\partial \tau} + \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} S = \nabla^2 S, \end{gather}

(2.10)\begin{gather}\boldsymbol{U} ={-}\boldsymbol{\nabla} P - \textit{Ra}\,S \hat{\boldsymbol{Z}}, \end{gather}

which is controlled by the Rayleigh number

(2.11)\begin{equation} \textit{Ra} = \dfrac{\kappa \Delta\rho g}{\mu E}. \end{equation}

Here, $\textit {Ra} = \mathcal {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-\partial _ZS$. Finally, we note that the rescaled pressure,

(2.12)\begin{equation} P = \dfrac{\kappa}{\varphi \mu D} (p+\rho_0 g z), \end{equation}

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 $\boldsymbol {U}$), with the salinity approaching the limit of $S = 0$ at large depths. These conditions are analogous to the thermal problem studied by Wooding (Reference Wooding1960) and Homsy & Sherwood (Reference Homsy and Sherwood1976). They are appropriate to a dry salt lake as long as there is no significant ponding of surface water (Wooding et al. Reference Wooding, Tyler and White1997; van Duijn et al. Reference van Duijn, Pieters, Wooding and van der Ploeg2002), although we will address the modifications to the model needed for ponding at the end of § 3.1. The stationary solution has $\boldsymbol {U} = (0,0,1)$ everywhere, a salinity boundary layer given by $S = \textrm {e}^Z$ and a corresponding pressure $P = \textit {Ra}(1-\textrm {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.

3. 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 (Reference Wooding1960) 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\leq 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 $\tilde {\boldsymbol {U}}$, $\tilde {S}$, $\tilde {P}$ to its velocity, salinity and pressure fields, respectively. The magnitudes of these perturbations are taken to be proportional to a small parameter, $\varepsilon _0$. To leading order in $\varepsilon _0$ (specifically, ignoring the $\tilde {\boldsymbol {U}}\boldsymbol {\cdot }\boldsymbol {\nabla }\tilde {S}$ term, which is of order $\varepsilon _0^2$) the perturbations will grow or decay according to

(3.1)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{U}} = 0, \end{gather}

(3.2)\begin{gather}\frac{\partial\tilde{S}}{\partial \tau} + \frac{\partial\tilde{S}}{\partial Z} + \tilde{U}_Z \frac{\partial S_0}{\partial Z} - \nabla^2 \tilde{S} = 0, \end{gather}

(3.3)\begin{gather}\tilde{\boldsymbol{U}} + \boldsymbol{\nabla} \tilde{P} + \textit{Ra}\,\tilde{S}\hat{\boldsymbol{Z}} = \boldsymbol{0}. \end{gather}

Here, the only remaining term related to the base state is the base salinity $S_0 = S_0(Z,\tau )$, 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 $\boldsymbol {\nabla }\times (\boldsymbol {\nabla }\tilde {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

(3.4)\begin{equation} {\nabla}^2\tilde{U}_Z + Ra(\partial_X^2 + \partial_Y^2) \tilde{S} = 0. \end{equation}

Solving for the dynamics of the perturbations requires providing for the appropriate boundary conditions. In order to respect the original boundary conditions, $\tilde {U}_Z = \tilde {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\rightarrow -\infty$. Then, rearranging (3.2) gives

(3.5)\begin{equation} \tilde{U}_Z = \partial S_0/\partial Z \left(\nabla^2 \tilde{S} - \frac{\partial\tilde{S}}{\partial Z} - \frac{\partial\tilde{S}}{\partial \tau}\right), \end{equation}

which implies that

(3.6)\begin{equation} \frac{1}{\partial S_0/\partial Z} \frac{\partial\tilde{S}}{\partial Z}\rightarrow 0 \quad \textrm{for} \ Z\rightarrow-\infty. \end{equation}

Now, we can look at the evolution of the salinity perturbation, $\tilde {S}$. For this we follow the approach of Pellew & Southwell (Reference Pellew and Southwell1940) as well as Wooding (Reference Wooding1960) and assume a separation of variables, such that

(3.7)\begin{equation} \tilde{S}(X,Y,Z,\tau) = F(Z) \varPhi(X,Y)\ \textrm{e}^{\alpha \tau}, \end{equation}

where $\varPhi$ is a harmonic function satisfying $(\partial _X^2 + \partial _Y^2 + k^2)\varPhi = 0$. Here, $k$ is the characteristic wavenumber of the perturbation in the horizontal directions and $\alpha$ is its growth rate. For $\alpha > 0$ the amplitude of the perturbation increases and the system is unstable, whereas for $\alpha < 0$ the perturbation decays and the system is stable; the $\alpha = 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)$,

(3.8)\begin{equation} (\partial_Z^2 - k^2)\left[\frac{1}{\partial S_0/\partial Z}(\partial_Z^2 - \partial_Z - k^2 - \alpha) F\right] = \textit{Ra}\,k^2F. \end{equation}

At this point, following the structure of (3.5) and (3.8), we will also introduce

(3.9)\begin{equation} G(Z,\tau) := \frac{1}{\partial S_0/\partial Z}(\partial_Z^2 - \partial_Z - k^2 - \alpha ) F(Z). \end{equation}

At $Z = 0$ the boundary conditions of $F = G = 0$ follow from the fact that $\tilde {U}_Z=\tilde {S}=0$ there, and from (3.5). In § 3.1 we will solve this problem for the stationary base state of $S_0 = \textrm {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 = \textrm {e}^Z$. Wooding (Reference Wooding1960) solved a similar problem for the special case of neutral stability ($\alpha = 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)\begin{gather} c_{0/1} = 1 \pm k, \end{gather}

(3.11)\begin{gather}c_{2/3} = \frac{1\pm \sqrt{1 + 4 (k^2 + \alpha)}}{2} = \frac{1}{2} ( 1 \pm \varPsi), \end{gather}

where $\varPsi = \sqrt {1 + 4 (k^2 + \alpha )}$. The solutions then read

(3.12)\begin{equation} F_{i}(Z) = ({-}k^2\, \textit{Ra}\ \textrm{e}^Z)^{c_i}\,H_i(Z)\quad \text{for }i \in \{1\ldots 4\}, \end{equation}

where $H_i(Z)$ is defined by

(3.13)\begin{gather} H_{0/1}(Z) = \,{}_0\mathrm{F}_3(;\{2c_{0/1}-1,c_{0/1}+c_2,c_{0/1}+c_3\};-k^2\,\textit{Ra}\ \textrm{e}^{Z}), \end{gather}

(3.14)\begin{gather}H_{2/3}(Z)=\,{}_0\mathrm{F}_3(;\{2c_{2/3},c_{2/3}+k,c_{2/3}-k\};-k^2\,\textit{Ra}\ \textrm{e}^{Z}). \end{gather}

These are hypergeometric functions of the form ${}_r \mathrm {F}_s$ with $r=0$ and $s=3$ (Koekoek & Swarttouw Reference Koekoek and Swarttouw1998; Askey & Daalhuis Reference Askey and Daalhuis2010). They can be evaluated as series, for example,

(3.15)\begin{equation} H_{0}(Z) =\sum _{n = 0}^{\infty}\left[ \frac{(- k^2\,\textit{Ra}\ \textrm{e}^{Z})^n}{n!(2c_0-1)_n(c_0+c_2)_n(c_0+c_3)_n} \right], \end{equation}

where $(a)_n = a(a+1)(a+2)\cdots (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 Reference Koekoek and Swarttouw1998, 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 - \varPsi /2+k=0,-1,-2\ldots$, while $H_2$ will converge unless the same condition holds for $1/2 + \varPsi /2-k$. Hence, the solutions form a dense set in the $(\alpha ,k)$ space.

From (3.6) it follows that $F(Z)$ has to decay faster than $\textrm {e}^Z$ for $Z\rightarrow -\infty$ (Wooding Reference Wooding1960). Therefore, of the four possible series described above, only those with $c_i>1$ are allowed, since $F_i \propto \textrm {e}^{Zc_i}$. As $k\geq 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+\alpha >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 $F(0) = C_0F_0(0) + C_2F_2(0) = 0$ and $G(0) = C_0G_0(0) + C_2G_2(0) = 0$. Consequently, we know that either $C_0 = C_2 = 0$ or

(3.16)\begin{equation} \left|\begin{matrix} F_0(Z) & F_2(Z) \\ G_0(Z) & G_2(Z) \end{matrix}\right|_{Z=0} = F_0(0) G_2(0) - F_2(0)G_0(0) = 0. \end{equation}

Since only the non-trivial solution is physically relevant, applying this constraint gives a relationship between $\textit {Ra}$, $k$ and $\alpha$ 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 $\alpha$, and hence find the smallest $\textit {Ra}>0$ satisfying the boundary conditions (e.g. for the neutral stability curve, $\alpha = 0$). Finally, we note that Wooding (Reference Wooding1960) 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 $\partial _Z \tilde {U}_Z=0$ there. For that scenario the arguments and solutions presented above remain valid, but require the modified boundary conditions of $F = \partial _Z G = 0$ at $Z = 0$. The determinant in (3.16) is similarly revised, to read $F_0(0) \partial _Z G_2(0) - F_2(0)\partial _ZG_0(0) = 0$.

Figure 3. Linear stability analysis and behaviour of numerical simulations near the critical point. (a) Neutral stability curves (solid lines), critical points at $\textit {Ra}_{c}$, $k_{c}$ (filled circles) and the most unstable mode (dashed lines), according to the linear stability analysis for a constant salinity at the upper surface and either uniform vertical flow (blue) or constant pressure (red) there. For all subsequent panels, only the uniform flow problem is considered. (b) Normalised eigenfunctions, $F(Z)$, for the most unstable mode, calculated at various $\textit {Ra}$. (c) Theoretical prediction of growth rates $\alpha$ as a function of wavenumber $k$ for a selection of $\textit {Ra}$ close to the onset of instability (solid lines). Their maximum values, $\alpha _{m}$, are indicated as red triangles. Measurements of the growth rates of single-mode perturbations ($\alpha _{fit}$) in the corresponding simulations are indicated as crosses, as discussed in appendix A (see figure 13). (d) Theoretical values for the growth rate of the most unstable mode, $\alpha _{m}$, for different conditions. Near the critical point, where $\epsilon = (\textit {Ra} - \textit {Ra}_{c})/\textit {Ra}_{c}$, the growth rate follows a power-law scaling, $\alpha _{\mathrm {m}} = \epsilon ^\gamma$, but begins to deviate from $\gamma = 1$ when $\epsilon \gg 0$ (see inset).

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, $\textit {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 (Reference Wooding1960). For the constant flux case, the corresponding results are consistent with Homsy & Sherwood (Reference Homsy and Sherwood1976) and van Duijn et al. (Reference van Duijn, Pieters, Wooding and van der Ploeg2002). The most unstable mode calculations provide additional predictions, which we will use to validate our numerical model of convection.

Table 1. Critical Rayleigh number, $\textit {Ra}_{c}$, and wavenumber, $k_{c}$, from linear stability analysis.

The eigenfunctions, $F(Z)$, corresponding to the most unstable modes of scenarios of different $\textit {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 (Reference Wooding1960) for constant pressure conditions and at the critical point. Furthermore, the eigenmodes do not undergo any significant qualitative changes as $\textit {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 $\alpha$ 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 Reference Cross and Greenside2009) 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 how they are modified at larger amplitudes, in other words in the nonlinear regime of the dynamics.

Finally, and introducing $\epsilon = (\textit {Ra} - \textit {Ra}_{c})/\textit {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, $\alpha _{\mathrm {m}}$, scales linearly with $\epsilon$ near the critical point, as expected for a type-I instability (Cross & Greenside Reference Cross and Greenside2009). Given this relationship, when we make quantitative comparisons of e.g. velocities at different $\textit {Ra}$ in what follows, we will often find it convenient to re-scale time, and thus define $\hat {\tau } = \tau \epsilon$.

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 (Reference Slim and Ramakrishnan2010), 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 $\tau =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. Reference Wooding, Tyler and White1997)

(3.17)\begin{equation} S_0(Z,\tau) = \textrm{e}^{Z/2}\left[\frac{1}{2}\ \textrm{e}^{Z/2}\, \mathrm{erfc}\left(\frac{-Z-\tau}{2\sqrt{\tau}}\right) + \frac{1}{2}\ \textrm{e}^{{-}Z/2}\, \mathrm{erfc}\left(\frac{-Z+\tau}{2\sqrt{\tau}}\right) \right]. \end{equation}

As shown in figure 4(a), this solution relaxes to the stationary base state, $\textrm {e}^Z$, over a time scale of $\tau \sim 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.

Figure 4. Linear stability analysis of time-dependent solution. (a) For the case of an evaporating salt lake with an initially constant salinity everywhere ($S=0$), and where the surface conditions are changed to $S=1$ at $\tau =0$, a time-dependent boundary layer will develop (dashed lines). This layer will relax to the stationary base state $\textrm {e}^Z$ over a time scale $\tau \sim 1$. (b) The linear stability of this boundary layer was evaluated to determine the first moment of instability (right $y$-axis), along with the first unstable mode (left $y$-axis). The most unstable mode of the stationary base state (from figure 3c, corresponding to late-time solution) is included for comparison. For (c) $\textit {Ra}=20$ and (d) $\textit {Ra}=40$ we show how the growth rates of the spectrum of unstable modes depend on time.

Equation (3.8) can be rearranged into an eigenvalue equation of the form ${\boldsymbol{\mathsf{A}}}F = \alpha {\boldsymbol{\mathsf{B}}}F$,

(3.18)\begin{equation} \left((\partial_Z^2 - k^2)\left[\frac{1}{\partial S_0/\partial Z}(\partial_Z^2 - \partial_Z - k^2) \right]-\textit{Ra} k^2\right)F = \alpha(\partial_Z^2 - k^2) \left[\frac{1}{\partial S_0/\partial Z}\right]F. \end{equation}

The largest eigenvalue here gives the growth rate of the most unstable perturbation for any particular mode $k$ at any instant $\tau$ and Rayleigh number $\textit {Ra}$. We solved this eigenvalue problem numerically, using a Chebyshev differentiation matrix (Trefethen Reference Trefethen2000) for the differential operators, $\partial _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 ($\tau = 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\simeq 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 $\textit {Ra}$. In figure 4(b) we show how the time of the first instability, and the value of the first unstable mode, depend on $\textit {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 $\tau$ and $k$. Examples are shown in figure 4(c,d) for the cases of $\textit {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 $\textit {Ra} = 40$, modes between $k=0.212$ and $3.412$ are unstable at $\tau = 10$, and this range is effectively coincident with the range of unstable modes as $\tau \to \infty$ (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 $\textit {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 Reference Slim and Ramakrishnan2010), 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. Reference Rapaka, Chen, Pawar, Stauffer and Zhang2008; Slim & Ramakrishnan Reference Slim and Ramakrishnan2010). 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.

4. 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 = \textrm {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 (Reference Riaz and Meiburg2003), Ruith & Meiburg (Reference Ruith and Meiburg2000) and Chen & Meiburg (Reference Chen and Meiburg1998a), are given in appendix A and the code itself is available on a public repository (Lasser & Ernst Reference Lasser and Ernst2020).

The simulations were first validated against the theoretical predictions of growth rates. Specifically, for $\textit {Ra}$ between 15 and 40 we added small perturbations of a single wavenumber to the initial conditions. The growth rate of this mode, $\alpha _{fit}(k,\textit {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 $\alpha _{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 $\textit {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 $\textit {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 $\textit {Ra} = 100$ and at different times $\tau$ 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 $\textit {Ra} = 30$, $\textit {Ra}=100$ and $\textit {Ra}=1000$, respectively.

Figure 5. Snapshots of the relative salinity $S$ at different times $\tau$ for a simulation with $\textit {Ra} = 100$. Panels show (a) the initial salinity distribution at $\tau = 0.01$, the linear growth regime at (b) $\tau = 0.1$ and (c) $\tau = 0.2$, the flux-growth regime at (d) $\tau = 0.3$ and (e) $\tau = 0.5$, the merging regime at (f) $\tau = 0.75$ and (g) $\tau = 1$ and the re-initiation regime at (h) $\tau = 3$. Panel (i) shows the horizontally averaged salinity distributions for the snapshots (a–d). The domain size of the simulation was $40\,L\times 40\,L$, and only part of the domain is shown in each case. The scale is the same for all snapshots, and given in panel (e).

As the simulations proceed they pass through several distinct regimes of dynamics, which can be related to those given by Slim (Reference Slim2014) 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 (Reference Slim2014), 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 $\textit {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 (Lasser et al. Reference Lasser, Nield, Ernst, Karius, Wiggs and Goehring2019).

4.1. 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 $\textit {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(\tau ) = \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 $\textit {Ra} > 100$ the initial peak can overshoot its plateau value significantly, as can be seen in figure 6(b) for a simulation at $\textit {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.

Figure 6. Development of the maximum plume speed, $V$, for different driving conditions and where $\epsilon = (\textit {Ra} - \textit {Ra}_{c}) / \textit {Ra}_{c}$. (a) Evolution of the plume speeds in simulations at different $\textit {Ra}$ show an initial transient, followed by a peak in velocity as the first plumes grow and merge and a subsequent steady state condition. (b) For the case of $\textit {Ra} = 1000$ we show results from a single simulation as a solid line, while the standard deviation within the corresponding ensemble is given by the shaded area. The blue dot, horizontal red dashed line, vertical grey dashed line and the grey rectangle indicate the measured peak speed $V_p$, the dynamic steady state speed $V_{stst}$, the beginning of the dynamic steady state $\hat {\tau }_{stst}$ and the time period used to calculate $V_{stst}$, respectively. (c) The time taken until the system reaches a steady state, $\hat \tau _{stst} = \tau _{stst}\epsilon$, settles down to a constant value of $\hat \tau _{stst}\approx 2.5$ for large enough $\textit {Ra}$ and $\epsilon$. (d) Development of $V_\mathrm {p}$ and $V_{stst}$ with $\epsilon$, showing how the plume velocity in the steady state scales linearly with the proximity to the critical conditions, $\epsilon$, whereas the initial overshoot of the plume speeds becomes more apparent at larger $\textit {Ra}$.

For every simulation we identified the peak speed, $V_\mathrm {p} = \max (V(\tau ))$, as shown for example by the blue dot in figure 6(b). We also estimated the beginning of the subsequent dynamic steady state, $\tau _{stst}$, as the time after which $\partial V/\partial \tau < 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 $\tau _{stst} \approx 2.5/\epsilon$ (or, alternatively, $\hat \tau _{stst} \approx 2.5$), for large enough $\epsilon$. 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 $\textit {Ra}$, whereas the growth rate of the most unstable mode is proportional to $\epsilon$. Although the time needed for the disturbances to saturate will depend on the initial amplitude of these perturbations, the scaling of $\tau _{stst}$ demonstrates that $\epsilon$ 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 $\epsilon$, we also calculated a steady state speed, $V_{stst}$, as the average of $V$ in the time window between $\tau _{stst}$ and $\tau _{stst} + 7 \epsilon$. 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 $\textit {Ra}$ in our study, where the shut down of convection also approaches fastest. For example, for the $\textit {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 $\epsilon$, and hence $\textit {Ra}$. This shows that the characteristic plume speed in the steady state, i.e. the re-initiation regime, increases linearly with $\epsilon$. It also shows how the overshoot of the plume speeds becomes more significant at larger $\textit {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 $\epsilon$ (Ernst Reference Ernst2017). 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 $\epsilon$. This can be readily understood, given that $\textit {Ra}= \mathcal {V}_{B}/E$ represents the ratio of the characteristic speeds of buoyancy and evaporation.

4.2. 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 $\hat \tau = \tau \epsilon$, we measured the effective wavelength $\varLambda =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{\rm \pi} /\varLambda$ 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 $\textit {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 $\hat \tau = 1$ to characterise the initial response, we find that the wavenumber of $k=3.55\pm 0.54$ seen at that time is similar to the most unstable mode of the linear stability analysis, $k=3.01$.

Figure 7. Time evolution of the plume wavelength, $\varLambda$, and wavenumber, $k$. The observed wavenumber depends on the depth $Z$ and time $\hat \tau$ where it is measured, as shown in (a) for a simulation at $\textit {Ra}=100$. Here, the red dashed curve gives $k$ as measured close to the top boundary, at $Z=-1$. This metric captures the initially high $k$ values before coarsening occurs, but also fluctuates at later times due to the re-initiation and merging of proto-plumes. Deeper, the blue dashed curve shows the wavenumber at $Z = -10$, where it has a much more stable value (at least, after the plumes have first descended to this level). (b) A space–time diagram is shown for the same simulation, where each horizontal line of data represents the state of the simulation for a moment at a depth of $Z=-1$. (c) Over time, all simulations show coarsening behaviour. Initially, the $k$ measured near the surface, at $Z=-1$, is similar to the most unstable mode of the linear stability analysis (dashed line). However, as the system ages $k$ decreases and ultimately shows little variation with $\textit {Ra}$ in the steady state regime, depicted here for the late-time case of $\hat \tau = 30$ and at a depth of $Z=-10$. Data for times and depths intermediate between these two limiting cases are shown as variously coloured circles. Note that wavenumbers for high $Ra$ and $\hat {\tau }$ are not available, since simulations at high $\textit {Ra}$ were usually terminated before the required duration was reached.

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 $\textit {Ra}=100$ simulation shown in figure 7(a,b) this happens between approximately $\hat {\tau } = 1$ and 10. From $\hat {\tau }\approx 10$ onward the plume wavenumber measured deeper into the simulated domain approaches a stable value of $k \approx 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 \approx 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 (Reference Slim2014), who referred to them as proto-plume pulses.

In order to quantify the coarsening of the plume wavelength more generally, we measured $\varLambda$ 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 ($\hat \tau =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 $\textit {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 $\textit {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} \simeq 0.76$. This is highlighted by the blue circles in figure 7(c), which show the wavenumbers as measured at $\hat \tau = 30$ and a depth of $Z=-10$. We will argue in the following section that the relative independence of $k$ on $\textit {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.

4.3. 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, $\langle S(Z)\rangle$, as measured in the dynamic steady state regime of simulations at different $\textit {Ra}$. Near the surface these results all demonstrate a rapid decay of the salinity with depth, which is stronger for higher $\textit {Ra}$. Figure 8(a) also shows how, just below this boundary layer, and especially for higher $\textit {Ra}$, the salinity may also pass through a small local maximum of around $S\approx 0.3$. A similar salinity peak can be noticed in Slim (Reference Slim2014) (figure 3). At intermediate depths, below $Z \approx -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 $\textit {Ra}$ (see appendix A). Similarly, for very high $\textit {Ra}$ the salinity starts building up below $Z\approx -9$, since the bottom boundary of the simulated domain is only at $Z=-10$ in these cases.

Figure 8. The dynamics of the salinity distribution and its effective boundary layer thickness, $L'$, depend on $\textit {Ra}$. (a) The horizontally averaged salinity, $\langle S(Z)\rangle$, is shown for a range of $\textit {Ra}$. Results are averaged over 5–10 ensembles for every $\textit {Ra}$ and taken from times within the dynamic steady state. The inset shows how an exponential fit is made to the near-surface results (filled dots), to determine $L^\prime$. (b) This thickness evolves over time and in the steady state regime the boundary layer fluctuates around an average value (dashed lines) that depends on $\textit {Ra}$. Here, each data set corresponds to a single simulation run. The inset highlights the initial evolution of the boundary layer during the linear growth and flux growth regimes. (c) The value of $L^\prime$ in the dynamic steady state is inversely proportional to $\textit {Ra}$, over a wide range of conditions. In other words, in this regime the buoyancy forces available in the boundary layer are always only slightly above what is required to trigger an instability. (d) Data collapse showing the salinity distribution in the boundary layers of simulations in their steady state, scaled by the Rayleigh number.

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 $\textit {Ra}$. To this end, we fit an exponential decay function $S(Z) = \Delta S \exp (Z/L^\prime ) + S_{\mathrm {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^\prime =1$. The results for the evolution of $L^\prime$ 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^\prime \simeq 1$ for short times (inset to figure 8b). For higher $\textit {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 $\hat \tau$ 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^\prime$ 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^\prime$, depend on the $\textit {Ra}$ of the system: for higher $\textit {Ra}$ the boundary layer is thinner and the fluctuations are more rapid.

For convection driven by one side, Slim et al. (Reference Slim, Bandi, Miller and Mahadevan2013) suggested that (in our notation) an effective Rayleigh number of $\textit {Ra} (1-S_{bulk} )L^\prime$ would describe a boundary layer that had been disturbed by, for example, loss of material to downwelling plumes. Essentially, $(1-S_{bulk})L^\prime$ gives the relative buoyancy forces available in the boundary layer, whereas $\textit {Ra}$ characterises the system's general ability to act on those forces. In all our simulations, the thickness of the effective boundary layer, $L^\prime$, appears to scale approximately with $1/\textit {Ra}$ in the steady state regime. More specifically, as shown in figure 8(c), the ratio $\textit {Ra} (1-S_{bulk})L^\prime /\textit {Ra}_{c}$ approaches a constant value of about 1.8 for $\textit {Ra}$ above about 100. As further evidence of this scaling, figure 8(d) shows a data collapse of the near-surface salinity distributions, $\langle S(Z)\rangle$, as rescaled by $\textit {Ra}/\textit {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$.

4.4. 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^\prime )$, where $L_0^\prime =1$, and with random perturbations characterised by an amplitude of $\eta = 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^\prime =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.

Figure 9. Development of the salinity boundary layer for different initial conditions. In all cases, salinity distributions were averaged horizontally and over ensembles of 9–11 simulations run at $\textit {Ra} = 100$. (a) Development of the near-surface salinity distribution in a simulation with our typical initial condition of $L^\prime _0 = 1.0$, showing its approach to a statistically steady state. (b) Early-time development of the boundary layer thickness $L'$ for different initial conditions $L^\prime _0$. (c) Development of $L^\prime$ starting with $L^\prime _0=0.4$ for different initial noise amplitudes $\eta$. (d) Development of $L^\prime$ for simulations started at the same initial conditions as in (b), but showing longer time scales. Values converge towards the corresponding steady state value (dashed line) seen in figure 8(b).

In figure 9(b) we show the development of the effective boundary layer $L^\prime$ for a selection of initial conditions $L_0^\prime$. 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^\prime =0.4$. When $L^\prime _0<1$ the boundary layer first grows towards the stationary solution of $L^\prime =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 $\hat \tau \approx 1$.

For all these cases, the simulations converge towards the same salinity distributions above times of $\hat \tau \approx 2$. This is emphasised in figure 9(d), which gives the late-time behaviour of $L^\prime$ for the simulations with different initial boundary layer thicknesses. Here, all simulations show the same value of $L^\prime$ 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.

5. 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-\varphi D\partial _z S$. This flux vanishes in the stationary solution corresponding to our initial condition, $S = \exp (ZE/\varphi 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 $\rho _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 $\partial _zS$ 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 $\textit {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}\approx 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.

Figure 10. Effect of a modulated evaporation rate on convection. (a) Space–time diagrams are shown for four simulations at $\textit {Ra}=100$ but with surface evaporation rates modulated at an amplitude of $A_{m} = 1$ and (from top to bottom) wavenumbers $k_{m} = 0.63$, 1.26, 1.88 and 2.51. These give the time evolution of the salinity profiles at a depth of $Z=-1$. In each case the evaporation rate at $Z=0$ is indicated by the sine wave on top of the diagram; the high-salinity downwellings show a preference for areas of lower evaporation. (b) Development of the order parameter $\xi$ with re-scaled time $\hat {\tau }$ for the same four $k_{m}$, also at a depth of $Z=-1$, and averaged over ensembles of five simulation runs.

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

(5.1)\begin{equation} \xi ={-} \frac{1}{N} \sum_{i=1}^N \cos(k_{m} X_i), \end{equation}

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 $\xi = 1$ then the plumes are perfectly aligned with minima in the surface evaporation, whereas when $\xi =-1$ the plumes all lie below maxima and if there is no preferred arrangement of plumes, then $\xi = 0$. A large $\xi$ 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 $\xi$ 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 $\xi$ up to about 0.9 at early times. Since $\textit {Ra} = \mathcal {V}_{B}/E$, lower-evaporation regions will appear, locally, as having a higher effective $\textit {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 $\xi$ 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 $\xi$ 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 $\textit {Ra}=100$ for homogeneous boundary conditions (see figure 7c). For the early-time ordering we define a phase-locking time, $\hat \tau _\phi$, as the time it takes for the order parameter to drop below $\xi = 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 $\textit {Ra}=100$ this is $k\approx 3.1$.

Figure 11. Effects of surface evaporation modulation at (a) early and (b) late times, for simulations at $\textit {Ra}=100$ and averaged over ensembles of 5 realisations for each set of conditions. Plume locations are determined by maxima in the salinity at a depth of $Z=-1$. (a) The duration of any initial phase locking, $\hat \tau _\phi$, depends on the modulation amplitude, $A_{m}$, and wavenumber, $k_{m}$. In particular, high-amplitude modulation at wavenumbers between 1.9 and 3.1 can noticeably delay plume merging. (b) Long-term behaviour of the average order parameter $\langle \xi \rangle$ for different modulation amplitudes and wavenumbers. This was calculated as the time-averaged value of $\xi$ in the range of $30 \leq \hat {\tau } \leq 60$. Surface evaporation modulation with $k_{m}$ between 0.3 and 0.9 is particularly effective at pinning downwelling plume locations, even down to relatively modest amplitudes of $A_{m}=0.3$.

To measure the long-time synchronisation of the plumes with the surface modulation we instead averaged the order parameter $\xi$ from $\hat \tau =30$ to 60. Figure 11(b) shows how this average order parameter, $\langle \xi \rangle$, 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, $\xi$ 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} \approx 0.76$ as well as the wavenumbers of polygonal salt crust patterns observed in nature, where $k_{crust} = 0.78\pm 0.43$ with no observable dependence on $\textit {Ra}$ (Lasser et al. Reference Lasser, Nield, Ernst, Karius, Wiggs and Goehring2019, Reference Lasser, Nield and Goehring2020).

5.1. 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. Reference Eloukabi, Sghaier, Nasrallah and Prat2013; Nield et al. Reference Nield, Neuman, O'Brien, Bryant and Wiggs2016; Bergstad et al. Reference Bergstad, Or, Withers and Shokri2017; Farhat Reference Farhat2018; Nachshon et al. Reference Nachshon, Weisbrod, Katzir and Nasser2018), 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 ${\pm }0.01\,^\circ$C and $\pm$0.1 %, respectively. The resulting data are deposited on a public repository (Nield, Lasser & Goehring Reference Nield, Lasser and Goehring2020) and the protocol for their collection is described in more detail by Lasser et al. (Reference Lasser, Nield and Goehring2020), 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\,^\circ$C at 2:00 PM and $0\,^\circ$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 (Lasser et al. Reference Lasser, Nield and Goehring2020). 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. Reference Nield, Bryant, Wiggs, King, Thomas, Eckardt and Washington2015). Additionally, the highest temperatures recorded during the day are similar between the centre of the crust polygon and the ridge.

Figure 12. Temperature (a) and relative humidity (b) measurements from salt polygons at Owens Lake, California, conducted in winter 2016. Sensors were embedded either within a ridge or in the centre of a nearby polygon. Displayed values are averaged over 20 minute intervals.

Temperatures inside the ridge, however, drop faster and to lower values at night; the temperature difference is approximately $2\,^\circ$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. Reference Nachshon, Weisbrod, Katzir and Nasser2018 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\,^\circ$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 Reference Farhat2018). 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} \geq 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.

6. 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 Reference Lasser2019; Lasser et al. Reference Lasser, Nield, Ernst, Karius, Wiggs and Goehring2019). When rescaled by the evaporation rate $E$ and a length scale that balances advection and diffusion, $L=\varphi D/E$, the model is controlled by a single dimensionless group, the Rayleigh number $\textit {Ra}$. In this context, $\textit {Ra}$ can be interpreted as the speed at which a large blob of salt-rich fluid would naturally descend, $\mathcal {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 Reference Wooding1960; Homsy & Sherwood Reference Homsy and Sherwood1976; Wooding et al. Reference Wooding, Tyler and White1997; van Duijn et al. Reference van Duijn, Pieters, Wooding and van der Ploeg2002). 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 $\textit {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 (Reference Slim2014) 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, $\epsilon = (\textit {Ra} - \textit {Ra}_{c})/\textit {Ra}_{c}$. In the linear growth regime the growth rate of the most unstable mode $\alpha _m\sim \epsilon$, as expected for a type-I instability. Similarly, for later regimes the plume velocities, measured in various ways, were found to be proportional to $\epsilon$. This follows naturally from the definition of $\textit {Ra}=\mathcal {V}_{B}/E$, since the velocity scaling of our system is based on the evaporation rate $E$. Given this, we found a time scale $\hat \tau = \tau \epsilon = \epsilon t E^2/\varphi ^2D$ to be a convenient rescaling. Specifically, the various regime transitions occur at times of order $\hat \tau =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 $\textit {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 $\textit {Ra}$, such that the wavenumber of the mature plumes always approached the critical wavenumber of $k_{c} \approx 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 $\textit {Ra} =\mathcal {V}_{B} H/\varphi D$ (e.g. Hewitt et al. (Reference Hewitt, Neufeld and Lister2014) who suggest $k\sim \sqrt {\textit {Ra}}$ or Fu, Cueto-Felgueroso & Juanes (Reference Fu, Cueto-Felgueroso and Juanes2014) who suggest $k\sim \textit {Ra}$). In a two-sided system, if a plume naturally moves at a speed $\sim \mathcal {V}_{B}$ across the domain height $H$, then a plume spacing of order $1/\sqrt {\textit {Ra}}$ reflects the distance over which concentration gradients would diffuse over the course of the plume's fall (Liang et al. Reference Liang, Wen, Hesse and DiCarlo2018). 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 $\textit {Ra}=\mathcal {V}_{B} L/\varphi 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 $\textit {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^\prime$, is proportional to $1/\textit {Ra}$ over a wide range of $\textit {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 $\textit {Ra}$, for Rayleigh numbers far enough above $\textit {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 Reference Christiansen1963; Krinsley Reference Krinsley1970; Nield et al. Reference Nield, Bryant, Wiggs, King, Thomas, Eckardt and Washington2015). A plume wavelength of $2{\rm \pi} /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. Reference Stevens, Sharp, Simmons and Fenstemaker2009; Van Dam et al. Reference Van Dam, Simmons, Hyndman and Wood2009). We intend to develop this argument further with detailed comparisons with field data elsewhere (Lasser et al. Reference Lasser, Nield, Ernst, Karius, Wiggs and Goehring2019, Reference Lasser, Nield and Goehring2020). 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.