## 1. 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 Reference Briere2000), 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 Reference Lowenstein and Hardie1985; Gill Reference Gill1996; Briere Reference Briere2000). The heat fluxes through the surface of such salt deserts are important to understanding the water and energy balances in arid regions (Bryant & Rainey Reference Bryant and Rainey2002). Additionally, salt deserts are responsible for a significant part of the global emission of atmospheric dust (Gill Reference Gill1996; Prospero Reference Prospero2002; Washington *et al.* Reference Washington, Todd, Middleton and Goudie2003). However, despite the extreme conditions that prevail above ground, the water table of dry lakes often remains very near to the surface (Gill Reference Gill1996; Briere Reference Briere2000; Bryant Reference Bryant2003; Nield *et al.* Reference Nield, Bryant, Wiggs, King, Thomas, Eckardt and Washington2015), allowing for active patterns of fluid flows within the pore spaces of the soil (e.g. Tyler *et al.* Reference Tyler, Kranz, Parlange, Albertson, Katul, Cochran, Lyles and Holder1997; Wooding, Tyler & White Reference Wooding, Tyler and White1997; Stevens *et al.* Reference Stevens, Sharp, Simmons and Fenstemaker2009; Van Dam *et al.* Reference Van Dam, Simmons, Hyndman and Wood2009). As evaporation rates are high (Tyler *et al.* Reference Tyler, Kranz, Parlange, Albertson, Katul, Cochran, Lyles and Holder1997; DeMeo *et al.* Reference DeMeo, Laczniak, Boyd, Smith and Nylund2003; Brunner *et al.* Reference Brunner, Bauer, Eugster and Kinzelbach2004; Groeneveld, Huntington & Barz Reference Groeneveld, Huntington and Barz2010), 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 Reference Lasser2019; Lasser *et al.* Reference Lasser, Nield, Ernst, Karius, Wiggs and Goehring2019). 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 (Lasser *et al.* Reference Lasser, Nield, Ernst, Karius, Wiggs and Goehring2019; Lasser, Nield & Goehring Reference Lasser, Nield and Goehring2020) or the Salt Lake Desert of Utah (Christiansen Reference Christiansen1963) to the Great Salt Desert of Iran (Krinsley Reference Krinsley1970), Salar de Uyuni in Bolivia and the Sua Pan of Botswana (Nield *et al.* Reference Nield, Bryant, Wiggs, King, Thomas, Eckardt and Washington2015) – 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.* Reference Nield, Bryant, Wiggs, King, Thomas, Eckardt and Washington2015). 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 Reference Raupach, Gillette and Leys1993; Marticorena & Bergametti Reference Marticorena and Bergametti1995). Christiansen (Reference Christiansen1963) and Krinsley (Reference Krinsley1970) 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 Reference Krinsley1970; Lowenstein & Hardie Reference Lowenstein and Hardie1985; Lokier Reference Lokier2012; Lasser *et al.* Reference Lasser, Nield, Ernst, Karius, Wiggs and Goehring2019). 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 (Lasser *et al.* Reference Lasser, Nield, Ernst, Karius, Wiggs and Goehring2019). 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 $\textit {Ra}$ close to the critical value, $\textit {Ra}_{c}$, as well as the behaviour at higher Rayleigh numbers. Here, once the system is appropriately scaled, $\textit {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.* Reference Metz, Davidson, de Coninck, Loos and Meyer2005; Neufeld *et al.* Reference Neufeld, Hesse, Riaz, Hallworth, Tchelepi and Huppert2010; Slim & Ramakrishnan Reference Slim and Ramakrishnan2010; Slim *et al.* Reference Slim, Bandi, Miller and Mahadevan2013; Slim Reference Slim2014; Thomas, Dehaeck & De Wit Reference Thomas, Dehaeck and De Wit2018; Hewitt, Peng & Lister Reference Hewitt, Peng and Lister2020) 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 Reference Hewitt2020). For example, analogously to the study by Slim (Reference Slim2014), 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 Reference Hewitt2020). Of special interest has been the critical value of the Rayleigh number, $\textit {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 $\textit {Ra} < 4{\rm \pi} ^2$ (Horton & Rogers Reference Horton and Rogers1945; Lapwood Reference Lapwood1948). At higher $\textit {Ra}$, first steady and then perturbed convective rolls occur (Busse & Joseph Reference Busse and Joseph1972; Graham & Steen Reference Graham and Steen1994). For $\textit {Ra} \gtrsim 1300$, the dynamics enters a chaotic regime (Otero *et al.* Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004; Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2012).

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 $\textit {Ra}_{c}\approx 14.35$ (Wooding Reference Wooding1960; Homsy & Sherwood Reference Homsy and Sherwood1976; van Duijn *et al.* Reference van Duijn, Pieters, Wooding and van der Ploeg2002) with a critical wavenumber of $k_{c} \approx 0.76$ (Wooding Reference Wooding1960). Finite amplitude perturbations in the range of $8.59 < \textit {Ra} < 14.35$ were also observed to grow in simulations and experiments (Wooding *et al.* Reference Wooding, Tyler and White1997; van Duijn *et al.* Reference van Duijn, Pieters, Wooding and van der Ploeg2002). 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 $\textit {Ra}_{c} \approx 6.95$ and $k_{c}\approx 0.43$ are smaller and thus the system is unstable for a greater range of parameters (Wooding Reference Wooding1960). A linear stability analysis (Wooding Reference Wooding1960) as well as an energy minimisation method (van Duijn *et al.* Reference van Duijn, Pieters, Wooding and van der Ploeg2002) were also used to calculate the neutral stability curves for both problems and to determine what range of wavenumbers are unstable at any given $\textit {Ra}$. In what follows we start by reproducing these theoretical results using an approach that develops particularly from the work of Wooding (Reference Wooding1960). 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.

## 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,

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

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 Meiburg1998*b*), 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

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,

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),

which is controlled by the Rayleigh number

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,

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

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

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

which implies that

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

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)$,

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

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

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

where $H_i(Z)$ is defined by

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,

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

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$.

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.

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)

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.

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

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 Meiburg1998*a*), 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.

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 4*a*). 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 4*b*). 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_\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 6*b*) 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$.

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 7*a* for the example of $Z = -10$). The wavenumber closer to the top boundary (e.g. red dashed line in figure 7*a*, 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 3*a*), 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.

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 8*b*). 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 6*b*). 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.

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.

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 $\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 7*c*). 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$.

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.

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.

## Supplementary movies

Supplementary material and movies are available at https://doi.org/10.1017/jfm.2021.225.

## Acknowledgements

We thank C. Beaume for discussions and a close reading of the manuscript; M. Threadgold for sharing the eigenvalue solver code adapted for § 3.2; J. M. Nield for assistance with field work; and A. Fourrière for early discussions on pattern formation mechanisms.

## Declaration of interests

The authors report no conflict of interest.

## Author contributions

M.E. and J.L. contributed equally to this work.

## Appendix A. Numerical simulation

We implemented a two-dimensional finite difference model of an evaporating salt lake, based on the non-dimensional system of equations given by (2.8)–(2.10). This uses a pseudo-spectral approach and a streamfunction–vorticity representation, similar to Chen & Meiburg (Reference Chen and Meiburg1998*a*), Ruith & Meiburg (Reference Ruith and Meiburg2000) and Riaz & Meiburg (Reference Riaz and Meiburg2003), a sixth-order compact finite difference scheme to compute spatial derivatives and an explicit fourth-order Runge–Kutta scheme for time stepping. The code has been made available on GitHub (Lasser & Ernst Reference Lasser and Ernst2020).

#### A.1. Model set-up

The numerical model simulates a two-dimensional $(X,Z)$ area of width $W$ in the $X$-direction and height $H$ in the $Z$-direction. For the salinity $S$ and fluid flux $\boldsymbol {U} = (U_X,U_Z)$ we assume periodic boundary conditions in the $X$-direction. The relative salinity $S = 1$ at the top boundary ($Z=0$) and $S = 0$ at the lower boundary ($Z=-H)$. At $Z=0$ surface evaporation is modelled as a boundary condition on $U_Z$. To allow for a modulation in 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 (Reference Ruith and Meiburg2000) and Riaz & Meiburg (Reference Riaz and Meiburg2003). For a two-dimensional and incompressible flow the flux $\boldsymbol {U}$ is related to the Lagrange streamfunction $\psi$ by $\boldsymbol {U} = (\partial _Z \psi ,-\partial _X\psi )$. The vorticity, $\omega = \partial _XU_Z-\partial _ZU_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 $\nabla \times (\nabla 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 $\psi$ and vorticity $\omega$. The constraints on flow give $\partial _X\psi = -1+A_{m}\cos (k_{m}X)$ at $Z=0$ and $\partial _X\psi = -1$ at $Z = -H$. These inhomogeneous boundary conditions can be accounted for by taking $\psi = \psi ^\prime + \varPsi$, where

such that $\partial _X\varPsi = 0$ at both $Z=0$ and $Z=-H$. The constant salinity conditions correspond to a vanishing vorticity, $\omega = 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 (Reference Riaz and Meiburg2003) 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 $\sigma = 3$ grid cells, such that

where $\eta$ gives the magnitude of the perturbation. The initial salinity field is then $S = S_0 + S^\prime$, at $\tau = 0$. Unless otherwise stated we used a default value of $\eta =0.05$.

#### A.2. Implementation

To numerically solve the governing equations of the simulation we largely follow the implementation described by Riaz & Meiburg (Reference Riaz and Meiburg2003) and Ruith & Meiburg (Reference Ruith and Meiburg2000). 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 Reference Lele1992); 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\times N$ points at positions $(x_m,z_n)$. The grid spacings $\Delta X$ and $\Delta Z$ are varied with $\textit {Ra}$, since the spatial resolution has to be increased at higher $\textit {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 $\textit {Ra}$. The time step $\Delta \tau$ follows the Courant–Friedrichs–Lewy (CFL) condition (Courant, Friedrichs & Lewy Reference Courant, Friedrichs and Lewy1928),

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 coefficients

and analogous expressions for $\hat \psi ^\prime _{\xi }$ and $\hat \omega _{\xi }$, where $\xi$ can take integer values between $\pm M/2$. In terms of these Fourier coefficients (A1) may be written as

where $\hat \omega _{\xi }$ are calculated at each time step by a similar transformation of (A2), and where $\hat {\omega }_{\xi 0} = \nabla ^2\hat \psi ^\prime _\xi$ 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 $\partial _{Z}^2\hat {\varPsi }_\xi (z_n)$ following the implicit sixth-order compact finite difference scheme given by Carpenter, Gottlieb & Abarbanel (Reference Carpenter, Gottlieb and Abarbanel1993). 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 $\alpha _i$, $\beta _i$, $a_i$ and $b_i$ listed in table 3, as calculated by Tyler (Reference Tyler2007). At each time step we then solve for the Fourier components $\hat \psi$, which are inverted to give the streamfunction $\psi (x_m,z_n)$. The flux $\boldsymbol {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 (Reference Süli2003, 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 (Reference Ruith and Meiburg2000), Riaz & Meiburg (Reference Riaz and Meiburg2003), Chen & Meiburg (Reference Chen and Meiburg1998*a*) and Tan & Homsy (Reference Tan and Homsy1988). For this, instead of (A6) we perturbed our initial conditions with

using a perturbation of initial magnitude $\eta = 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 $\alpha _{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 $\textit {Ra}$ and $k$ agree with the theoretical values to within a relative error of order $1\,\%$, as was shown in figure 3(*b*).