On the Stability of Density Stratified Flow Below a Ponded Surface

Flooding of coastal areas with seawater often leads to density stratification. The stability of the density-depth profile in a porous medium initially saturated with a fluid of density ρf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _\mathrm{f}$$\end{document} after flooding with a salt solution of higher density ρs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _\mathrm{s}$$\end{document} is analyzed. The standard convection/diffusion equation subject to the so-called Boussinesq approximation is used. The depth of the porous medium is assumed to be infinite in the analytical approaches and finite in the numerical simulations. Two cases are distinguished: the laterally unbounded CASEA\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathbf {{\small {\uppercase {case~A}}}}}}}$$\end{document} and the laterally bounded CASEB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathbf {{\small {\uppercase {case~B}}}}}}}$$\end{document}. The ratio of the diffusivity and the density difference (ρs-ρf)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\rho _\mathrm{s} - \rho _\mathrm{f})$$\end{document} induced gravitational shear flow is an intrinsic length scale of the problem. In the unbounded CASEA\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathbf {{\small {\uppercase {case~A}}}}}}}$$\end{document}, this geometric length scale is the only length scale and using it to write the problem in dimensionless form results in a formulation with Rayleigh number R=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R = 1$$\end{document}. In the bounded CASEB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathbf {{\small {\uppercase {case~B}}}}}}}$$\end{document}, the lateral geometry provides another length scale. Using this geometrical length scale to write the problem in dimensionless form results in a formulation with a Rayleigh number R given by the ratio of the geometric and intrinsic length scales. For both CASEA\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathbf {{\small {\uppercase {case~A}}}}}}}$$\end{document} and CASEB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathbf {{\small {\uppercase {case~B}}}}}}}$$\end{document}, the well-known Boltzmann similarity solution provides the ground state. Three analytical approaches are used to study the stability of this ground state, the first two based on the linearized perturbation equation for the concentration and the third based on the full nonlinear equation. For the first linear approach, the surface spatial density gradient is used as an approximation of the entire background density profile. This results in a crude estimate of the L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}-norm of the concentration showing that the perturbation at first grows, but eventually decays in time. For the other two approaches, the full ground-state solution is used, although for the second linear approach subject to the restriction that the ground state slowly evolves in time (the so-called frozen profile approximation). Just like the ground state, the resulting eigenvalue problems can be written in terms of the Boltzmann variable. The linearized stability approach holds only for infinitesimal small perturbations, whereas the nonlinear, variational energy approach is not subject to such a restriction. The results for all three approaches can be expressed in terms of Boltzmann t\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{t}$$\end{document} transformed relationships between the system Rayleigh number and perturbation wave number. The results of the linear and nonlinear approaches are surprisingly close to each other. Based on the system Rayleigh number, this allows delineation of systems that are unconditionally stable, marginally stable, or transiently unstable. These analytical predictions are confirmed by direct two-dimensional numerical simulations, which also show the details of the transient instabilities as function of the wave number for CASEA\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathbf {{\small {\uppercase {case~A}}}}}}}$$\end{document} and the wave number and Rayleigh number for CASEB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathbf {{\small {\uppercase {case~B}}}}}}}$$\end{document}. A numerical example of the effect of a layer with low permeability is also shown. Using typical values of the physical parameters, the analytical and numerical results are interpreted in terms of dimensional length and time scales. In particular, an explicit stability criterion is given for vertical column experiments.


Introduction
Stability of density stratified fluids has been studied for more than 100 years. In the book 'The self-made tapestry', Ball (2001) gives in chapter 7 a lucid description of the history of the analysis of pattern formation in such fluids. The extension to fluids in porous media is nearly as old and has been reviewed thoroughly by Straughan (2008). At first, stability was explored primarily by analytical methods applied to the linearized perturbation equations. Later it became feasible to use analytical and numerical methods to study the full nonlinear equations.
Employing those modern methods, van Duijn et al. (2002) analyzed gravitational stability of a saline boundary layer below an evaporating salt lake. They compared stability criteria obtained by three different methods: the method of linearized stability; the traditional energy method using as constraint the integrated Darcy equation; an alternative energy method using as constraint the point-wise Darcy equation. The two constraints are, respectively, referred to as integral constraint and differential constraint. They showed that the integral constraint gives a stability bound of the Rayleigh number equal to the square of the first root of the Bessel function J 0 , in agreement with previous numerical results of Homsy and Sherwood (1976). The differential constraint gave results that were in excellent agreement with experiments by Wooding et al. (1997a, b). The stability of the saline boundary layer formed by upward flow was treated in greater detail by Pieters (2004), Pieters and Schuttelaars (2008).
In this paper, we use the same approach as in van Duijn et al. (2002) and Pieters (2004) to analyze stability following flooding with salt water of a soil initially saturated with less saline water. In the natural environment, such flooding is quite common. Marine transgression is a dominant feature of Holocene geo-hydrological history, starting almost 12 millennia ago. For example, see Raats (2015) for a historical perspective, the Dutch coastal region is shaped by marine transgression and human interference in response to this. In the region where now the freshwater IJsselmeer and the IJsselmeerpolders are located, there was in Roman times the freshwater Flevomeer, later known as Almere. Due to marine transgression, the inland lake gradually evolved into the inland Zuiderzee, which stood in open connection with the North Sea. From about 1600 AD onward, the bottom of Zuiderzee became more and more saline. With the double purpose of protection from floods and the creation of new polders, in 1932 the tidal inlet to the Zuiderzee was closed off by a dam. In a few years, the brackish Zuiderzee changed into the freshwater IJsselmeer and as a result of that the bottom started to desalinize. That process is still continuing.
Earlier studies on gravitationally unstable systems having a nontrivial base steady-state condition showed that this base condition may be unstable for a given system configuration. Perturbations to some existing condition will then lead to other nontrivial and persistent convective steady-state conditions. However, in the studies of Nield et al. (2008) and Selim and Rees (2010), as well as in the present study, the base steady-state condition only temporar-ily evolves to such nontrivial convective conditions. Eventually, the system will return to a new base steady-state condition which can be considered as a (concentration or temperature) redistributed version of the original base steady state. The redistributed base steady-state condition may again be unstable.
In the 1930s and early 1940s, again following Raats (2015), a group of civil engineers inspired by the mathematical physicist Burgers did pioneering work on transport of salts across the water-sediment interface. They formulated the linear convection-dispersion equation and solved this equation for appropriate boundary conditions. The calculated salinity profiles in the top 10 m were mostly in good agreement with profiles measured in the 1950s by van der Molen (see Raats 2015) and the references therein). The spatial distribution of the measured salinity profiles over a large area, such as the North-East polder, was used as an indicator of the spatial distribution of the flow velocities. However, not all field measurements were in line with the computations. In a few scattered places in the North-East polder higher salinities were found at depths of 10-15 m, far deeper than predicted and measured in other places. In some cases, this could reflect seepage of water from the former Zuiderzee to lower lying polders along the coast. In others, there was a highly permeable Pleistocene deposit reaching the surface. Van der Molen speculated that for the latter 'this phenomenon is probably due to convective currents in the bottom of the Zuiderzee between 1600 and 1931 A.D'. Specifically, he noted that the small difference in density between the freshwater present in the soil and the supernatant seawater is sufficient to cause convection currents. This is in line with the theory and computations presented in this paper.
On the timescale of centuries and under certain circumstances, marine transgression may cause rapid salinization of entire aquifers. In northwest Germany, Holocene transgressions of a few thousands of years have brought salt water of corresponding age to a depth of over 200 m (Mull and Battermanna 1980). Earlier Geirnaert (1973) had reported similar observations along the Dutch coast. In view of such results, it may seem surprising that at many other places all around the world fresh and brackish waters have been found beneath saline waters on the continental shelves (Post et al. 2013). Wherever this occurs, the deeplying freshwater is invariably protected from invading ocean waters by sediments with a low hydraulic conductivity. Post and Simmons (2009) illustrate by means of sand tank laboratory experiments and numerical modeling how low-permeability lenses protect freshwater from mixing with overlying saline water. This is also in line with the simulations presented in this paper.
On a seasonal time scale, exchange of substances between freshwater in sediments and supernatant saline water may occur. Smetacek et al. (1976) reported that in the water of the Kiel Bight high nutrient concentrations and low oxygen concentrations were found following an influx of higher density water. They suggested that this could be related to density-induced natural convection in the sediments. Webster et al. (1996) used a computational model and laboratory experiments to demonstrate that gravitational convection can make an important contribution to the exchange of water and solutes between sediments and a supernatant water column in regions subject to significant temporal variations in salinity, such as estuaries. In effect, they verified the suggestion put forward by Smetacek et al. (1976).
In coastal regions, tsunamis often cause flooding by seawater and natural convection may then strongly contribute to the mixing of this seawater with the underlying freshwater. Illangasekare et al. (2006) described the impact of the 2004 tsunami in Sri Lanka on groundwater resources. They included results of an experiment in a 53 × 30 × 2.7 cm PlexiglasTM tank filled with glass beads. Vithanage et al. (2008) used a laboratory model to study natural convection in coastal aquifers following flooding by sea water.
The analysis of the stability of density stratified flow below a ponded surface in this paper has important connections with a pioneering paper by Wooding (1962b). In that paper, Wooding uses the method of linearized stability to analyze the behavior of an initially sharp, horizontal interface separating two miscible fluids in an infinitely long, vertical, porous column. In the upper half of the column, the fluid has a higher density than in the lower half. He studied the time dependence of the linearized equations by means of an expansion in terms of Hermite polynomials. In spite of the rather primitive numerical tools available in those days, the approach undertaken by Wooding in this and other papers-both theoretical and experimental: see, e.g., Wooding (1959Wooding ( , 1960Wooding ( , 1962aWooding ( , 1969; Wooding et al. (1997a, b)-is still a source of inspiration for those studying stability problems in porous media. For this reason, we dedicate this paper in his honor.
Long before we all were made energy conscious in the 1970s, already in the mid-1950s Robin Wooding (1926Wooding ( -2007 (Wooding 1969). In 1970/71, he took leave without pay to spend a year at The Johns Hopkins University and the University of Wisconsin to work on flow and transport in porous media, after which he decided to return not to CSIRO but to DSIR, where he stayed until he retired in 1987. He then returned again to Canberra as Honorary Fellow at the CSIRO Pye Laboratory. There for the last two decades of his life his main interest was the salt lake problem mentioned above.
After this introduction, the outline of this paper is as follows. The mathematical model for studying the stability of density stratified flow below a ponded surface is formulated in Sect. 2, where we distinguish between laterally unconfined and confined domains. A special solution of this model is presented in Sect. 3. This solution describes the formation of a onedimensional salt layer that penetrates into the subsurface by diffusion/dispersion. It is the main goal of this paper to investigate the gravitational (in)stability of this layer. In particular, we study the influence of the system's dimensionless Rayleigh number and the role of time. To this end, we perturb the special solution and derive the corresponding perturbation equations. This is done in Sect. 4. Next, in Sect. 5, we investigate the stability by two different methods, namely the standard method of linearized stability, where we employ two approaches, and the energy method with differential constraint. Stability results for homogeneous soils are discussed in Sect. 6.1 and for nonhomogeneous (layered) soils in Sect. 6.2. An interpretation of the results in terms of the physical parameters of the problem is given in Sect. 7, and conclusions are presented in Sect. 8.

Remark 1
Technical details related to the stability analyses and the numerical methods are omitted in this paper. They are given in a comprehensive version of this work that appeared as the technical report Pieters et al. (2018).

Problem Formulation
We consider an isotropic porous medium occupying either CASE A: the laterally unbounded half-space Ω = Ω A = {z > 0}, or CASE B: the three dimensional, laterally confined, semi-infinite region Here z is the vertical coordinate pointing downwards. We assume that the permeability κ is constant in the x, y direction, i.e., κ = κ(z) only, and that the porosity φ is constant throughout Ω. Adopting the Boussinesq approximation (Oberbeck 1879;Boussinesq 1903;Chandrasekhar 1981;Landman and Schotting 2007), we consider the equations in Ω and for t > 0. In these equations, ρ denotes the fluid density, U the Darcy velocity, and P the fluid pressure. These are the unknowns in the problem that need to be determined. The (given) physical parameters are: D (dispersion/diffusion coefficient), μ (fluid viscosity), g (gravity constant), and ρ f (freshwater density). Both D and μ are assumed constant throughout this paper.
The boundary conditions at the top of the porous medium are given by: where ρ s is the salt (sea) water density and h the height of the overlying salt water layer. Boundary condition (2a) is a simplification of what happens in practice. The laboratory experiments by Webster et al. (1996) show that growing instabilities induce local outflow from the porous medium into the ponded layer. It is assumed in this paper that outflow across the top of the porous medium {z = 0} has no effect on the density of the fluid in the overlying layer (which is ρ s for all t > 0). The initial condition is given by: where the initial density satisfies ρ f ρ i ρ s .
Remark 2 (i) The results for CASE B can be generalized to arbitrarily shaped bounded cross sections. In Sect. 7, we consider the case of a circular cross section as well. (ii) Since Darcy's law contains the gradient of the pressure, one may replace P by P + P 0 , with P 0 constant in space, and obtain the same discharge and salt distribution. Hence, the specific value of P at the top z = 0 does not play a role. In other words, the discharge and salt distribution do not depend on the height h of the ponded layer. In essence, this due to the assumed constraint of incompressibility of the fluid.
It is assumed that ρ i → ρ f as z → ∞, implying that ρ → ρ f as z → ∞ for all t > 0. Thus ρ f acts as a far field reference value. We further assume that as z →∞, the pressure gradient force is balanced by the gravitational force. Thus, we consider those solutions of (1) that satisfy: ρ → ρ f and ∂ P ∂z → ρ f g as z → ∞, for all t > 0.
For the horizontally bounded CASE B, it is assumed that on lateral boundary Γ v (see Fig. 2) the fluxes of water and solute vanish: where n is the unit vector on Γ v . Using the constant densities ρ s and ρ f , we define: Here, P is the pressure based on the far field reference density ρ f . Introducing (3) into (1) and (2) gives subject to: and for CASE B also Next, we introduce the characteristic scales for length (L c ), time (T c ), pressure (P c ), discharge (Q c ) and permeability (κ c = κ(0)). Using these scales, we substitute the dimensionless variables into Eqs. (4) and conditions (5) to obtain in the scaled domain Ω * : subject to and for CASE B also The transformed region Ω * is given by and in (7) and (8) the Rayleigh number R is given by Note that in (7b) we introduced Balancing the first and third term in transport equation (7a) gives for the timescale T c : and similarly, balancing all terms in Darcy's law (7a) gives the pressure scale P c : and discharge scale Q c : The discharge scale Q c can be interpreted as a density difference induced, gravitational shear flow (de Josselin de Jong 1960). The timescale T c and pressure scale P c are known once the length scale L c has been selected. In the laterally unbounded CASE A, there is no obvious length scale in terms of the spatial dimensions of the problem. Therefore, we use the intrinsic length scale L c = D/Q c . In the laterally bounded CASE B, we choose L c = L x /π or some other characteristic horizontal length scale. The factor π is for notational and technical convenience. All scales are summarized in Table 1.
For CASE B, the Rayleigh number is, apart from the factor π, the ratio of the characteristic geometric length scale L x and the intrinsic length scale D/Q c . Its value may vary considerably due to L x and, in particular, κ c (see Sect. 7).
The scaled region Ω * is given by Thus, after scaling, we obtain the dimensionless form of the equations, dropping the superscript *, in Ω and for t > 0, subject to and for CASE B also U · n = 0, (RUC − ∇C) · n = 0 along Γ v , for t > 0.
The transformed region Ω is given by Remark 3 Suppose one compares different values of R in CASE B. These differences are due to differences in κ c and L x . Note that if L x changes, then κ(z) changes accordingly, because of the scaling of z. This is expressed by (11). For CASE A, κ c only appears in the length scale.
subject to Applying (18b) in (17b) shows that U 0 = 0 for z → ∞. In view of (17c), this means that U 0 = 0 for all z. Using this in (17) gives: Solving (19a) subject to (18) yields This one-dimensional solution describes the formation of a salt layer that enters the subsurface by diffusion. Integrating (19b) subject to (18a) and (20a) results in Typical profiles for various t > 0 are shown in Fig. 1. Note that Equation (21) reflects the fact that ponding with a layer of fluid having constant density ρ s , results in the half-space Ω being filled with stagnant fluid of density ρ s as t → ∞.

The Perturbation Equations
To investigate the stability of the ground state {C 0 , U 0 , P 0 } we write the original variables as where u = (u, v, w) T . Here c, u and p are the perturbations whose behavior we need to study. For this, we have to first write the perturbation equations. Substituting (22) in Eq. (14), the initial condition (15a) and the boundary conditions (15b) give subject to and for CASE B also In view of the properties of the ground-state solution, there results for the perturbations the equations and the boundary conditions and for CASE B also Multiplying (25b) by κ(z) and applying the divergence gives Using this identity, the pressure boundary condition in (26a) implies for w the Neumann condition Following Lapwood (1948), we take the Laplacian of the vertical component of the perturbation Darcy equation (25b). This gives Further, taking the divergence of the perturbed Darcy equation (25b) and using (25c), yielding upon differentiation with respect to z Subtracting (30) from (29) results in The Neumann condition (28) and identity (31) will be used frequently in the next sections.

Stability Analyses
In the unbounded CASE A, we consider x, y-periodic solutions of the perturbation equations with respect to x, y-periodic cells of the form where T ⊂ R 2 is an element (tile) of the horizontal periodic structure. This implies a redefinition of Ω A as introduced in Sect. 2. As noted in Chandrasekhar (1981) or Straughan (2008), T can be a rectangle (or square), a strip or a hexagon. In this paper, the focus is on nonlinear stability and on the onset of infinitesimal small perturbations (linear approach). Details about the characterization of the periodic structure as for instance in Mielke (2002), who studied the spatial structure in the Rayleigh-Bénard problem, are beyond the scope of this paper. Some aspects, though, are discussed in Pieters et al. (2018).
In the linear approach, the horizontal periodicity is described by a function f = f (x, y) that satisfies where a is a characteristic wavenumber. If T is a hexagon with side L one has (see Chandrasekhar 1981) f (x, y) = cos 4π 3L and In this paper, we visualize the horizontal periodicity in the form of rectangles T = (0, x ) × (0, y ), where x and y are for the moment unspecified. This yields periodic cells of the form In (35), x , y are dimensionless lengths. A sketch of a two-dimensional cross section of the flow domain in the x, z-plane is shown in Fig. 2.
In the x, y-periodic assumption, periodic boundary conditions are applied along the vertical boundary of Ω A , which we denote by Γ v . This means that c, u, and p can be extended in a periodic way over the whole Ω, and belong to C ∞ (Ω).
Sketch of a two-dimensional cross section of the flow domain in the x, z-plane. In the laterally unbounded CASE A: Periodic boundary conditions for the concentration c, the discharge u, and pressure p apply at vertical boundary Γ v = Γ 2 ∪ Γ 4 . In the laterally bounded CASE B: No-flow and no-transport conditions apply at the vertical boundary Γ v . For numerical simulations, we consider a truncated column of depth d In the bounded CASE B, we consider the perturbation equations in the flow domain Ω = Ω B , which we write as where = L y L x is a dimensionless length ratio. Here, we use no-flow and no-transport conditions along the vertical boundary Γ v .
For computational purposes, see Sect. 6.1.3, we use the depth-truncated versions of Ω A and Ω B : and where d denotes the dimensionless depth. In Sect. 5.1, we consider the method of linearized stability, valid for small perturbations, and in Sect. 5.2, the nonlinear, variational energy method.

Linearized Stability
Assuming the perturbations to be small, one disregards the nonlinear term u · ∇c in (25a). The result is the perturbation system in and with Note that (39) is a linear system in terms of c, w (z-component of Darcy velocity), and p. Because of the linearity and the nature of the conditions (40a) and (40b), we may assume that solutions are of the form Then (40c), (40d) are satisfied provided the wavenumbers satisfy Substituting (41) in (39) gives for the amplitudes c(z, t), w(z, t), and p(z, t) the equations for z > 0 and t > 0, subject to conditions (40a), (40b). Here a 2 = a 2 x + a 2 y . We will study the linearized problem (43) for all a > 0 and R > 0. However, when interpreting the results in terms of stability we consider (CASE A) a unspecified and R ≡ 1, (CASE B) a = a n,m = n 2 + (m/ ) 2 and R = R S , where R S denotes the system Rayleigh number as specified in Table 1.
It suffices to consider Eqs. (43a) and (43b) for c and w only, since p follows directly from equation (43c).
In the following Sects. 5.1.1 and 5.1.2, we use two distinct methods to obtain linear stability results.

Direct L 2 -estimate
To obtain an estimate for the L 2 -norm of the concentration perturbation, we multiply equation (43a) by c and integrate the result in z. This gives the expression for each t > 0.
Since the ground-state saturation C 0 satisfies we can estimate (44) by, disregarding the second term in the left side as well, Introducing the L 2 -norm (i.e., ||c|| 2 2 = ∞ 0 c 2 dz) and dropping t from the notation, this inequality is written in the compact form At this point, we use the Cauchy-Schwarz inequality ∞ 0 |cw| dz ||c|| 2 ||w|| 2 , to estimate (46) further as Next, we multiply (43b) by w and integrate the result. Using the Neumann condition (28), this yields Now let κ = max z 0 κ(z). Then, we deduce from (48) where we used again Cauchy-Schwarz. Thus Applying this inequality in (47) gives Integrating this expression in t yields the estimate By (49), the same holds for ||w(t)|| 2 . Next, we apply (49) in (48) Using this inequality and (43c), we obtain where κ = min z 0 κ(z). Thus both the pressure p and the vertical flow w decay like c. Estimate (51) says that x, y-periodic perturbations, characterized by the wave number a, eventually decay in time if the initial amplitude c(0) ∈ L 2 ([0, ∞)) (i.e., ∞ 0 c(0) 2 dz < ∞). This is due to the term −a 2 t in the exponent. In essence, this is a diffusive term. This can be seen when going back to the dimensional form where now a and t are dimensional. A sketch of the upperbound in (51) is shown in Fig. 3. Note that the maximum is reached at The stability result expressed by (51) was derived using the rather crude estimate (45). In the next section, we do not use (45). Instead, we consider the growth or decay of small perturbations of C 0 (z, t) for each t > 0.
For future reference, we introduce Using this in (50) gives the form  (51) implying that if R < R M (a, t), then L 2 -stability of the linearized system is guaranteed.

Eigenvalue Formulation
Again we start from the linear system (43). In case of a steady, i.e., time independent, ground state, as in Wooding et al. (1997a, b), one looks for amplitudes c and w that have an exponential growth rate in time, i.e., {c, w}(z, t) = {c, w}(z)e σ t . If σ < 0, then all small perturbations decay in time and we call the system linearly asymptotic stable. If, on the other hand, σ > 0, then small perturbations grow exponentially in time and we call the system unstable.
In this paper, the ground state depends on time as well. However, since ∂C 0 ∂z (z, t + τ ) ≈ ∂C 0 ∂z (z, t) for |τ | being sufficiently small, we proceed similarly as in the steady case. This is called the 'frozen-profile' approach (van Duijn et al. 2002). Hence, for given t > 0, we consider instead of (43a) the approximate equation for τ > 0, but sufficiently small. Now time t appears as a parameter in the equation. We then assume that the perturbations have an exponential growth rate denoted by σ , i.e., Substituting (56) into (55) and (43b) yields the following eigenvalue problem: For given a > 0, t > 0 and σ , find the smallest R > 0 (denoted by R = R 1 (a, t, σ )) so that the problem for z > 0, with has a nontrivial solution.
The eigenvalues now depend on a, t, and on the growth rate σ . In the salt lake problem (Pieters and van Duijn 2006), it was shown that the smallest positive eigenvalue R 1 (a, t, σ ) satisfies This expresses exchange of stability. We assume here that (58) also holds for the eigenvalue problem (57). Under this assumption, which is supported by numerical evidence in Pieters et al. (2018), it suffices to analyze (57) for the case of neutral stability, i.e., σ ≡ 0. For completeness, we write the eigenvalue problem resulting from linearized stability.
with c and w satisfying has a nontrivial solution.
Note that here a > 0 and t > 0 act as parameters. The smallest positive eigenvalue is denoted by R L = R L (a, t). We return to this problem in Sect. 6.1.1.
The behavior of large perturbations is studied in the next section.

The Variational Energy Method
The method of linearized stability holds for infinitesimal small perturbations. What can be said about the system behavior for larger perturbations. To address this point, we apply the variational energy approach, see van Duijn et al. (2002) for the specific application to the salt lake problem and Straughan (1992Straughan ( , 2008 for a general overview.

Remark 4
The energy method can be applied for both domain Ω A (CASE A) and Ω B (CASE B). In the derivation that follows, we use Ω i to denote the flow domain. The index i can refer to A or B. Both geometrical cases would give the same equations, although the interpretation of the stability results is different.
In the energy method, one estimates the time derivative of the L 2 -norm of the concentration perturbation c. In particular, the aim is to find the largest R-interval for which Here and in the integrals below, we disregard the infinitesimal volume elements in the notation. The related maximum R-value clearly will depend on the wavenumber a and, because C 0 = C 0 (z, t), on time t. Starting point here is the nonlinear system (25).
To investigate (60), we multiply (25a) by c and integrate over Ω i . Using (25c) and the boundary conditions along Γ v , we find the identity: Thus, if R is chosen such that the right-hand side of (61) is negative for all perturbations satisfying a given constraint which is closest to the solution of the problem, then stability is guaranteed.
Roughly, there exist two approaches to incorporate the constraint on the admissible perturbations. One is based on (25c) (see Homsy and Sherwood 1976), and the integrated Darcy equation (25b): The other approach, and the one we focus on in this paper, uses the point-wise differential expression (31) as constraint: Applied to the salt lake problem, van Duijn et al. (2002) showed that this gives sharper results than constraint (62). It leads to the following class of admissible perturbations: In the energy method, the goal is to determine the largest R-interval for which To this end, we consider the maximum problem The boundedness of the functional in this expression is a direct consequence of the choice of space H and the behavior of ∂C 0 ∂z . This is detailed in Pieters et al. (2018), "Appendix E".
Suppose this maximum problem has a solution (i.e., there exists such a R E ∈ (0, ∞)). Then any solution of (25), (26), and (28) satisfies and thus for any R < R E one obtains from (61) Hence, for R < R E , the system is unconditionally stable. Maximum problem (63) is equivalent to the following minimum problem: Find the small- To solve (64), we introduce the functional where π acts as a Lagrange multiplier. We need to find points (c, w, π), with (c, w) ∈ H and π is Ω i -periodic, lim z→∞ π(z) = 0, and values of R where the first variation of (65) vanishes. This leads to the Euler-Lagrange equations in Ω i . Moreover, ∂π ∂z (0) = 0 appears as a natural boundary condition. Details of the derivation are given in Pieters et al. (2018).
Thus far Ω i has been left general. Explicitly considering periodicity of the form (41) for c, π, and w, which is allowed since the equations in (66) are linear, yields the eigenvalue problem with c, w, and π satisfying has a nontrivial solution. Here, the smallest positive eigenvalue is denoted by R E = R E (a, t). We return to this problem in Sect. 6.

Homogeneous Soils
When considering homogeneous soils, we set κ(z) = 1 in problems (P L ) and (P E ). We show below that this makes it possible to eliminate the time t from the eigenvalue problems. As a consequence, we obtain a clear interpretation of the role of a, R S and t in the stability behavior.

Eigenvalue Approach
The interpretation of the eigenvalues in terms of the wave number a and time t is not straightforward. Taking t as parameter, the eigenvalues R L (·, t) and R E (·, t) as function of a > 0, do not show a clear ordering. However, since we can remove t from (P L ) and (P E ) by introducing the variables Applying (68) to (P L ) and (P E ) yields transformed eigenvalue problems in terms of the same similarity variable as the ground state. The resulting problems read with c and w satisfying has a nontrivial solution, with c, w, and π satisfying has a nontrivial solution.
Clearly, now R * L and R * E depend on b only. This allows for a more transparent interpretation of the stability results.
To find the smallest possible eigenvalues for b > 0, we solve (P * L ) and (P * E ) by means of a Chebyshev-Petrov-Galerkin spectral method. Details and references are given in Pieters et al. (2018). The result is shown in Fig. 4. Note that R * L (b) > R * E (b) for all b > 0, but the difference is small for b > 0.5. This is different from the salt lake problem (van Duijn et al. 2002), where there is a significant gap between R L and R E . In this gap, small perturbations vanish, but large perturbations may grow in time. Since for all practical purposes, see Fig. 4, the ponded problem considered in this paper has a sharp stability bound: below R * E there is unconditional stability and above R * L instability. Remark 6 Using Green's functions, one can write w and π in (P * L ) and (P * E ) in terms of c. Substituting this into (69a) and (70a), and applying again the Green's function, yields eigenvalue problems in terms of integral equations for c. These transformations and some mathematical properties of the integral forms are detailed in Pieters et al. (2018). In particular, since c(0) = dw dη (0) = dπ dη (0) = 0, the corresponding integral equations are nonselfadjoint.
We show in Pieters et al. (2018) that both R * L (b) and R * E (b) have the lower bound: . This estimate is obtained by manipulating the equations in (P * L ) and (P * E ). Recalling that, see (53) with κ = 1, R M = R M (a, t) = a 2 √ πt, we find after applying (68), Hence, the stability bound obtained by the direct L 2 -estimate of the linearized equations is a lower bound for both R * L (b) and R * E (b).

Stability Interpretation
How to interpret and estimate the (in)stability from Fig. 4? For simplicity, we assume that (71) holds. The fact that R * L (b) > R * E (b) close to b = 0 requires only a minor modification of the argument. Let a > 0 be a given wavenumber and let R S be the system Rayleigh number. In view of (68), we have Hence, R * S (b) is a straight line passing through the origin in the (b, R * )-plane. With reference to Fig. 5, three cases can be distinguished, depending on whether the straight line lies (a) entirely below R * E , or (b) touches R * E at b * , or (c) intersects R * E at b * L and b * R . We successively consider cases (b), (a), and (c).
Writing R * (b) = αb, we have for α < α crit then the ground state is unconditional stable. How to interpret the result? CASE A : When considering a specific horizontal structure, for certain x and y , the wavenumber a has the discrete values a n,m = nπ then the ground state is L 2 -stable for perturbations that are x, y-periodic with respect to a rectangular structure with x and y satisfying (74). But perturbations with larger wavelengths (or smaller wavenumbers) may grow in time to form finger-like patterns. CASE B : Then a takes the discrete values a n,m = n 2 + m 2 n, m = 1, 2, . . . .

If the horizontal dimensions of the vertical column Ω B satisfy
then the ground state is L 2 -stable in Ω B implying that no fingers can occur as a result of initial perturbations. Fig. 5: elsewhere. Clearly, the intersection points b * L and b * R depend on a and R S : b * i = b * i (a, R S ) for i = L, R. For a given system, Rayleigh number R S , the boundary curves in the (a, t)-plane corresponding to a √ t = b * i (a, R S ) (i = L, R) are shown in Fig. 6 (Left: CASE A, Right: CASE B). These curves meet when R S a = 6.82, since then b * L = b * R = b * . Above the composite curve (red and blue curves in Fig. 5), one has L 2 -stability and below the curve small perturbations grow in time. Let a = a 1,1 , see (74) or (75), denote again the smallest wavenumber for a given Ω A or Ω B : i.e., the most unstable mode for a given geometry. Then, with reference to Fig. 6, perturbations decay in time when 0 < t < t L (a), they grow when t L (a) < t < t R (a), and finally again decay when t > t R (a). Here t i (a) = (b * i /a) 2 with i = L, R.
For t < t L (a) and t > t R (a) the system is stable (indicated by s), while for t L (a) < t < t R (a) the system has an episode of transient instability (indicated by u): perturbations having wave number a will first decay, then grow and eventually decay again in time In Sect. 7, we return to Figs. 5 and 6 for the interpretation in terms of the physical parameters.

Direct Numerical Simulations
To verify and interpret our theoretical findings we make a comparison with the nonlinear dynamics described by the full set of Eq. (14). To this end, we perform direct numerical simulations of (14) in two-dimensional cells: Having a two-dimensional flow domain in which ∇ · U = 0, allows us to formulate the convective flow part in (14) in terms of the stream function Ψ ; i.e., In the computations presented in this Section, we use κ(z) = 1, since the soil is considered homogeneous. Later, in Sect. 6.2, we allow κ = κ(z).

The salt transport equation in terms of Ψ becomes
At the top Γ 3 , we have P = 1. This means that q x = − ∂ P ∂ x = 0 at Γ 3 . Hence, from (76) At the bottom Γ 1 , we prescribe C = 0 and ∂ P ∂z = 0, see (14c). This means that 1 κ(d) q z = 0 at Γ 1 . Hence, again from (76), For CASE A, we impose periodic boundary conditions on C and Ψ along the vertical boundary Γ v . This means that for all t > 0 and For CASE B, we impose no-flow and no-transport conditions along Γ v , i.e., for t > 0 and We prescribe an initial concentration (77) and (78) and conditions (79)-(85) define the problem (N S P), the numerical solution problem. This problem has a unique solution {C, Ψ } for any reasonable initial concentration (85), e.g., see Clément et al. (1992).

Remark 7
Repeating Ω d A in horizontal direction yields a smooth 2 x -periodic solution in the strip Ω d := {0 < z < d}. If the initial concentration C i is symmetric with respect to x = 0, then the same holds for the evolving concentration C and anti-symmetry holds for the corresponding stream function Ψ : i.e., for (x, z) ∈ Ω d A and t > 0.
In view of the periodic boundary conditions, this implies for all t > 0. Thus for symmetric initial data, there is no transport across the vertical boundaries. Clearly, this does not hold for arbitrary initial data.
Problem (N S P) is solved by means of a characteristic-Galerkin finite element method. In the method, the convective operator is split from the diffusive one and is written in the form of the material derivative. A detailed description of this method is given in Pieters et al. (2018).
Let the flow domain Ω d i , with i = A, B, be discretized by a set of elements, forming a mesh M. The element size of the mesh M will be denoted by ΔM. Further, time is discretized by the time-grid t n := nΔt, n = 1, 2, . . . , N . We denote the numerical solution of problem (N S P) by { C n m , Ψ n m }, where m indexes over the elements in M, and n indexes over a discrete time grid.
To quantify the temporal behavior of the full system subject to an initial perturbation, we consider as a stability measure the discrete L 2 -norm in Ω i : In this definition, we do not use the analytical ground state C 0 evaluated per element and time-grid point, but rather its numerical approximation C n 0,m . We discuss (direct) numerical simulations for both CASE A and CASE B. For CASE A, an overview of the stability results is shown in Table 2. There we present for different wavenumbers a stability overview based on the eigenvalue approach (Fig. 6a) and on direct numerical simulations using the quantity E(t) (Fig. 7a). The stability behavior in Table 2 is derived from Figs. 6a and 7a. Figure 8 shows direct numerical simulations of the salt concentration for the cases 6, 7, and 8. In those simulations, we have used ΔM = π/5 ≈ 0.6, Δt = 0.05, and N = 50000 which corresponds to a (dimensionless) time horizon of 2500. Further, we perturb the initial concentration by a small symmetric (with respect to x = 0) perturbation of the form in which = 0.001. For CASE B, we have a fixed geometrical configuration (i.e., a 1,1 ≡ 1). An overview of the stability results is shown in Table 3. There we present for different Rayleigh numbers a stability overview based on the eigenvalue approach (Fig. 6b) and on direct numerical simulations using the quantity E(t) (Fig. 7b). This time the stability behavior in Table 3 is derived from Fig. 6b. Figure 9 shows direct numerical simulations of the salt concentration for the cases 1, 7, and 8. In those simulations, we have used ΔM = 0.15, Δt = 0.005, and N = 30,000 which corresponds to a (dimensionless) time horizon of 150. Figure 7a shows E over time for an x range as given in Table 2. For x = 5π, we have a monotonically decreasing E, which means that the solution remains close to the ground state C 0 (z, t) and the discrete L 2 (Ω A )-norm of the perturbation decays in time. Growth is observed for x > 7π: there exists a time t * and t * * for which Here, denotes the derivative with respect to time. The interval (t * , t * * ) is referred to as the unstable episode obtained by direct numerical simulations. Numerically, we find that it is consistent with the theoretical predicted range (t L (a 1,1 ), t R (a 1,1 )), as given by Fig. 6. t L (a 1,1 ) < t * < t * * < t R (a 1,1 ).
Also note that for those two system Rayleigh numbers, the used wavenumber of the initial perturbation (a = 1) is relatively far away from its critical value 2.2 and 2.93, respectively (see also Fig. 6, right). As a result, the numerical bounds are not sharp with respect to their corresponding theoretical bounds.   Table 2. b CASE B, with data from Table 3 Remark 8 As shown in Fig. 7, quantity E decays to zero for large times, indicating that the system will return to some one-dimensional profile. However, this profile is not the onedimensional ground-state solution due to the transient convective instability. The long-term one-dimensional profile can still be considered as the solution of (17), but with an adjusted initial condition.
The concentration profiles as function of time are presented in Fig. 8 for CASE A, and Fig. 9 for CASE B. In both cases, we observe for the selected cases transient instability. The instability introduces convective flow for a specific period of time. This convective dominated regime enhances the diffusive downward transport of salt into the porous column. The speed of the fingertips is higher for larger domains (CASE A) and larger Rayleigh numbers (CASE B). Along Γ 3 , near the corners, there is upward flow. Such upward flow was observed experimentally by Webster et al. (1996). The combination of the concentration boundary condition along Γ 3 and upward flow leads near the corners to a thin concentration boundary layer.

Layered Soils
Until now, we assumed for CASE A and CASE B homogeneous media only. We now consider, for CASE B only, a heterogeneous medium having three layers, in which one of the layers has a reduced permeability. The layered permeability structure is represented by Since the dimensional permeability is scaled by κ, we have κ 1.
In each vertical layer, the permeability is assumed to be constant. Hence, in each layer a 'local' system Rayleigh number, denoted by R κ S , can be defined. We consider two values for  The initial condition is again symmetrically perturbed by (86) with a 1,1 ≡ 1. The nature of the stream function ensures that continuity across the layer interface and yields continuity of the normal component of the fluid discharge. Hence, continuity of Ψ is the only condition imposed at those interfaces.
Numerical simulations for this three-layer setup are shown in Fig. 10 for a system Rayleigh number R S = 20. This Rayleigh number corresponds to a system with a transient instability, see also Table 3: The perturbation develops in layer I until the fingertip reaches z = d 1 , the start of the layer II. The finger pattern travels through layer II and as can be clearly seen for t = 4 and t = 5, the finger has become more narrow. For the special case in which layer II is almost impermeable, the convection cells in layer I develop and disappear over time again, whereas in layer III a reverse process takes place: No convection cells exist for small times, but for larger times they appear. However, for R S = 20 those convection cells in layer III are too weak to create a finger structure. Ultimately, the finger structure disappears and diffusion becomes the dominant process again.
The effect of κ is also reflected by the stability measure E(t), as shown in Fig. 11. For smaller values of κ, we observe a shorter episode of instability (t * , t * * ). Further, the magnitude of the growing instability is, as to be expected, less pronounced because the almost impermeable layer blocks the convective transport of the finger to deeper regions.

Interpretation in Terms of Dimensional Physical Parameters
Thus far the entire analytical and numerical analysis has been carried out in terms of dimensionless variables. In the concrete, dimensional world, the analysis covers a vast number of special cases. To interpret the results for specific cases, one has to return to dimensional variables by multiplying the dimensionless variables by the corresponding characteristic scales. We use subscript D to denote the dimensional variables. Specifically, multiplying the spatial dimensionless coordinates x, y, z by the characteristic length scale L c gives spatial dimensional coordinates x D , y D , z D and multiplying the dimensionless time t by the characteristic time T c gives the dimensional time t D : Also of interest are the dimensional wave number a D and wave length λ D given by To illustrate the conversion, we use the typical values of quantities listed in Table 4. Introducing the values of φ, ρ f , ρ s , g, μ, κ c and D from Tables 4 in 1 gives the characteristic    Table 1 the values for φ, ρ f , ρ s , g, μ and D from Table 4 Laterally unbounded CASE A Laterally bounded CASE B Table 6 Overview of parameters for CASE A and CASE B obtained by setting κ c = 10 −12 m 2 in Table 5 Laterally unbounded CASE A Laterally bounded CASE B scales for discharge (Q c ), length (L c ), time (T c ), pressure (P c ), and Rayleigh number R. If κ c is left to be chosen, then Table 5 results. Motivated by solute exchange by convection within estuarine sediments, Webster et al. (1996) made an experimental and numerical analysis of unstable penetration of ponded saline water into a layer of porous medium initially saturated with freshwater, using κ c = 5.4 × 10 −11 m 2 . Post and Kooi (2003), Post (2004) numerically modeled the rates of salinization in the high-permeability Pleistocene coastal aquifer of the Netherlands, using κ c = 10 −11 m 2 . More or less in line with those values of κ c , we consider in the first instance a porous medium with κ c = 10 −12 m 2 and subsequently infer results for coarser porous media with larger values of κ c and for finer porous media with lower values of κ c . Substitution of κ c = 10 −12 m 2 in Table 5 gives the values of Q c , L c , T c , P c and R in Table 6.

Laterally Unbounded CASE A
According to Table 5 (left column), for CASE A the Rayleigh number R = 1 and the characteristic scales Q c , L c , T c and P c are proportional to powers of the characteristic permeability κ c . To the wide range of values of κ c in Table 4 from 10 −10 m 2 for very coarse porous media to 10 −16 m 2 for very fine porous media correspond wide ranges of the values of the scales Q c , L c , T c and P c .
The approximate equality of the stability bounds from the linear and energy stability led for CASE A to the dimensionless condition (74) for L 2 -stability: In dimensional variables, this is equivalent to: where L c given by the left-hand side column in Table 5. Only perturbations with wave numbers that are smaller will grow and decay in a certain time interval t L,D = T c t L to t R,D = T c t R , where t L and t R can be read from Fig. 6a. The corresponding condition for the wave length is derived from (89): Only perturbations with wave lengths that are larger will grow and decay in a certain time interval Using the values in the left column of Tables 6 and 2 is translated to dimensional form by multiplying x and d by L c = 6 × 10 −3 m. Figure 7a is translated to dimensional form by multiplying the t-axis by T c = 8.4 × 10 3 s and x by L c = 6 × 10 −3 m. The dimensionless numerical results shown in Fig. 8 are be translated to dimensional form by multiplying the depth d = 250 and the column width x = 10, 11, and 12 by L c = 6 × 10 −3 m, giving d D = 1.5 m and x D = 6.0, 6.6, 7. For a coarser porous medium with a tenfold larger value κ c = 10 −11 m 2 , the values of x D = 6.0, 6.6, 7.2 mm are tenfold smaller and the values of t D = 0.2425, 0.4850, . . . , 3.3950 days are 100-fold smaller. For still coarser porous media, the width of the fingers will eventually be of the same order as the size of the pores, so that the limit of the Darcy model will be approached.

Laterally Bounded CASE B
According to Table 5 (right column), for CASE B, the characteristic length scale L c and pressure scale P c and the Rayleigh number R depend linearly on the horizontal length scale L x , while for the timescale T c this dependence is quadratic; further, the flux scale Q c and R are proportional to κ c . With the wide range of values of κ c in Table 4, the corresponding wide range of values of the Rayleigh number R is 53.05 × 10 2 L x for very coarse porous media to 53.05 × 10 −4 L x for very fine porous media.
For a finer porous medium with a tenfold smaller permeability κ c = 10 −13 m 2 , the values of L c are tenfold larger and the values of T c are 100-fold larger. For a coarser porous medium with a tenfold larger permeability κ c = 10 −11 m 2 , the values of L c are tenfold smaller and the values of T c are 100-fold smaller. The approximate equality of the stability bounds from the linear and energy methods led for CASE B to the dimensionless condition (75) for L 2 -stability: Here, represents the aspect ratio of a rectangular box (0, L x ) × (0, L y ), i.e., = L y /L x . In dimensional variables, this is equivalent to: where L c given by the right-hand side column in Table 5. Only perturbations with wave numbers that are smaller will grow and decay in a certain time interval t L,D = T c t L to t R,D = T c t R , where t L and t R , for a given R S , can be read from Fig. 6b. The corresponding condition for the wave length is derived from (93): Only perturbations with wave lengths that are larger will grow and decay in a certain time interval t L,D = T c t L to t R,D = T c t R . Criterion (93) means that in a rectangular column the ground state is L 2 -stable provided or, using Table 1, When the cross section is a square, i.e., = 1, L 2 -stability holds if In practice, experiments are often performed in cylindrical columns see, e.g., Wooding (1959). Suppose we consider in CASE B a semi-infinite cylinder of radius R. This radius will then be used, conform CASE B, in the characteristic length scale L c = R π . Then, instead of (36), we have Ω B = (x, y, z) : x 2 + y 2 < π, z > 0 .
As in (33), we need to consider now the eigenvalue problem In particular, for the given geometry, one needs to determine the smallest a > 0 for which the eigenvalue problem (97) has a solution. As explained in Sect. 6.1.2, this solution describes the most unstable mode.
Solutions of (97a) are found by separation of variables: where J m denotes the Bessel function of the first kind of order m. Boundary condition (97b) is satisfied when This yields a =â 1,1 = β 1,1 π as the smallest a > 0 for which (97) has a nontrivial solution. Usingâ 1,1 in (93) yields for cylinders the stability criterion Flow patterns corresponding to (98) are discussed in Pieters et al. (2018). The dimensionless time axis of Fig. 7 for the stability measures E is translated to a dimensional time axis by multiplying it by the characteristic time T c given by (92). Since T c is proportional to R 2 S , the tendency is that the transient instability occurs later for larger R S . The dimensionless numerical results shown in Fig. 9 are be translated to dimensional form by multiplying the dimensionless depth d = 75 and column width 2π by the linearly R dependent characteristic length L c given by (91) and multiplying the dimensionless times t = 1, 2, 3, 4, 5, 10, 15, 20, 30, 35 by the quadratically R dependent characteristic time T c is given by (92).

Conclusions
We analyzed the stability of the density-depth profiles in porous media initially saturated with a stagnant fluid of density ρ f after flooding with a salt solution of higher density ρ s . Such flooding of coastal areas with seawater regularly occurs, either due to climatological circumstances or as part of engineering projects. Historically, flooding was sometimes also imposed as a defensive measure.
We used the standard convection-diffusion equation subject to the so-called Boussinesq approximation. The variable density arising from the diffusive penetration of salt then affects the fluid flow only by its effect on the gravitational driving force in the Darcy equation. To explore the density difference-induced flow, we used both analytical and numerical methods. We assumed the depth of the porous medium to be infinite in the analytical approaches and finite in the numerical analyses. Geometrically, two cases were distinguished: the laterally unbounded CASE A and the laterally bounded CASE B. The ratio of the diffusivity D and the density difference-induced gravitational shear flow (κ c /μ)(ρ s − ρ f )g is an intrinsic length scale of the problem. In the unbounded CASE A, this is the only length scale and using it leads to a dimensionless formulation with Rayleigh number R = 1. In the bounded CASE B, the lateral geometry provides another length scale and using it leads to a dimensionless formulation with a Rayleigh number R given by the ratio of the geometric and intrinsic length scales.
The analytical approaches lead to general statements about stability or instability under given circumstances. The background diffusion of salt from the ponded layer into the fluid saturated porous medium, the so-called ground-state solution, is expressed in terms of the Boltzmann similarity variable z/ √ t. Three methods were used to study the stability of this ground state, the first two based on the linearized perturbation equation for the concentration and the third based on the full nonlinear perturbation equation. They yield eigenvalue problems that contain the horizontal wave number a and the Rayleigh number R as parameters. For homogeneous soils, using the Boltzmann variable z/ √ t, the time t can be eliminated from the eigenvalue problem by absorbing √ t into a and R: b = a √ t and R * = R √ t. This time-independence of the eigenvalue problems leads to the three stability curves relating the transformed Rayleigh number R * and wave number b shown in Fig. 4.
For the first linear approach, the surface spatial density gradient was used as an approximation of the entire background density profile. This results in a crude estimate of the L 2 -norm of the concentration shown in Fig. 3, with the perturbation at first growing, but eventually decaying in time.
For the second linear approach, the full ground-state solution subject to the restriction that the ground state slowly evolves in time was used: the so-called frozen profile approximation. This method of linearized stability only holds for infinitesimal perturbations. To get results for larger perturbations, a specific version of the variational energy method was used to estimate from the full nonlinear perturbation equation the L 2 -norm of the concentration perturbation c. The stability curve R * L (b) resulting from the linearized perturbation approach is surprisingly close to the stability curve R * E (b) resulting from the variational energy method. This means that there is only a small gap between stability curves R * L (b) and R * E (b) in which small perturbations vanish, but large perturbations may grow. In this respect, the current problem differs markedly from the problem with an imposed vertical background flow analyzed earlier (van Duijn et al. 2002), for which the gap was found to be large. The ponded problem considered in this paper has a sharp stability bound: below R * E there is unconditional stability and above R * L instability. The curve R * M = b 2 resulting from the first linear approach is shown to be a lower bound for both R * L and R * E . For both CASE A and CASE B, the √ t transformed system Rayleigh number R * S = R S √ t is a linear function of the √ t transformed system wave number b. To determine the stability for a particular wave number b -system Rayleigh number R S combination, three cases can be distinguished, depending on whether the straight line R * S = αb lies entirely below R * E , touches R * E at b = b * , or intersects R * E at b * L and b * R . Numerically, it was found that the line with α = α crit = 3.85 touches R * E at (R * , b) = (3.89, 1.01). For α < α crit the ground state is unconditionally stable. For α > α crit , the ground state is unconditionally stable for b < b * L and b > b * R and experiences an episode of transient instability for b * L < b < b * R . In Sect. 6.1.2, by reversal of the Boltzmann transformations, the analytical stability results were interpreted for both CASE A and CASE B in terms of system Rayleigh number R S , wave number a and time t (see in particular Fig. 6). In Sect. 6.1.3, these results were compared with direct numerical simulations based on the full, spatially 2-D flow and transport equations in x, z, t, with a range of lateral length scales for CASE A (see Table 2) and a range of system Rayleigh numbers for CASE B (see Table 3). Qualitatively, the numerical results turned out to agree with the eigenvalue results. They also clearly show the details of the transient instabilities as function of the wave number for CASE A and of the wave number and system Rayleigh number for CASE B. The transient instability clearly leads to deeper penetration of salt.
The case of a layered soil was considered in Sect. 6.2. There the Boltzmann transformation does not apply and one has to rely on numerical techniques. Results of the salt distribution in an aquifer having a horizontal and lower permeable layer are shown in Fig. 10. This figure clearly demonstrates the effect of the layer acting as a barrier for the flow and thus for the salt transport. This illustrates the protective role of low-permeability layers against invasion of saline water into underlying freshwater mentioned in the Introduction.
In Sect. 7, it was shown how the dimensionless analytical and numerical results can be translated to dimensional form and thereby interpreted in terms of dimensional length and time scales. This provided an opportunity to highlight the important role of the permeability of the porous medium: The higher the permeability the stronger is the destabilizing influence of the gravitational force.
The laterally bounded CASE B with its no-flow and no-solute transport conditions at the vertical boundary Γ ν is not restricted to rectangles and hexagons, but covers columns of any shape, such as the cylinders discussed at the end of Sect. 7. Thus, it provides a good basis for the interpretation of experiments in columns. In particular, column experiments are stable provided, see also (95) (iii) Circular cross section with radius R κ c (ρ s − ρ f )gR μD < 12.56.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.