On the nonlinear dynamics of a saline boundary layer formed by throughflow near the surface of a porous medium

Physica D 237 (2008) 3075–3088

Contents lists available at ScienceDirect

Physica D journal homepage: www.elsevier.com/locate/physd

On the nonlinear dynamics of a saline boundary layer formed by throughflow near the surface of a porous medium G.J.M. Pieters a,∗ , H.M. Schuttelaars b a

Centre for Analysis, Scientific computing and Applications (casa), Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands

b

Delft Institute of Applied Mathematics (diam), Delft University of Technology, P.O. Box 5031, 2628 CD Delft, The Netherlands

article

info

Article history: Received 26 February 2007 Received in revised form 31 March 2008 Accepted 11 June 2008 Available online 26 June 2008 Communicated by H.A. Dijkstra Keywords: Model reduction Galerkin projection Numerical bifurcation analysis Non-self-adjointness Flow in porous media

a b s t r a c t We consider gravitational instability of saline boundary layers, observed at the subsurface of salt lakes. This boundary layer is the result of the convective transport induced by the evaporation at the horizontal surface of a confined porous medium. When this upward transport is balanced by salt dispersion, a steady state boundary layer is formed. However, this boundary layer can be unstable when perturbed. This results in complex groundwater motion and density fields. The aim of this paper is to investigate the existence of finite amplitude solutions describing these resulting patterns (both the number of solutions and their structure), their stability, and their dependency on the system Rayleigh and Péclet numbers. For this purpose we construct a low-dimensional dynamical system (a reduced model) by projecting the nonlinear model equations onto a relatively small set of eigenfunctions of the problem linearized at criticality. The Galerkin projection approach is complicated by the fact that the problem under consideration is non-self-adjoint due to the existing evaporation. This implies that the eigenfunctions do not form an orthogonal set and therefore the adjoint eigenfunctions are used for the projection. The reduced model is constructed in such a way that it is capable of providing solutions in the strongly nonlinear regime as well. Convergence of these solutions towards the fully nonlinear model results is shown by means of direct numerical simulations. Further, the reduced model seems to partly capture the complex nonlinear behavior as seen in Hele–Shaw experiments by Wooding et al. [R.A. Wooding, S.W. Tyler, I. White, P.A. Anderson, Convection in groundwater below an evaporating salt lake: 2. evolution of fingers or plumes, Water Resour. Res. 33 (6) (1997) 1219–1228]. The physical transition mechanism that explains the occurrence of some observed bifurcation types is presented as well. © 2008 Elsevier B.V. All rights reserved.

1. Introduction Salt lakes are observed in arid and semiarid environments, for example in Central Australia and the High Plains of Texas, USA. The accumulation of salt minerals at the subsurface of some of these lakes is a result of upward water flow induced through the lake bed to supply the evaporative demand at the lakes. Dispersive transport takes place in the opposite direction of the evaporation induced convection. As a result a saline boundary layer is formed at the surface. Disturbances of this boundary layer can result in a complex groundwater motion that is organized in convection cells. These cells are believed to result in regular hexagonal patterns that are often observed at the surface of the salt lake. Wooding et al. [1] were the first to mathematically formulate the above-sketched process. From this point on we will refer to it

∗

Corresponding author. E-mail addresses: [email protected] (G.J.M. Pieters), [email protected] (H.M. Schuttelaars). 0167-2789/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.physd.2008.06.004

as the Wooding problem. Using both experimental and numerical methods, they showed that the stability of the boundary layer can be inferred from the system Rayleigh number. It is defined as the ratio of the gravitational convective flow speed and the upward evaporation. The permeability of the water-bearing underground layer (for saline lake environments approximately 10−14 m2 ) and the evaporation rate (typically 10–100 mm/yr) at the lakes surface are the principal controlling factors. For small enough Rayleigh numbers, the evaporative flux is strong enough to stabilize the unstable density distribution of the boundary layer. For larger Rayleigh numbers, convective motion induced by the unstable density profile dominates and strongly affects the transport of salt and, whenever present, contaminants. The Wooding problem as described above is closely related to the Lapwood problem [2], which is the porous medium equivalent of the well-known Rayleigh–Bénard convection problem for free fluids. The Lapwood problem is the limiting case for the Wooding problem when the evaporation goes to zero. Although the Wooding problem and the Lapwood/Rayleigh–Bénard problem are closely related from a physical point of view, there is an important

3076

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

difference between the two problems: The Wooding problem is non-self-adjoint whereas the Rayleigh–Bénard problem is selfadjoint. The asymmetry of the Wooding problem results in a more complex and unexpected behavior at the onset of instability and in the strongly nonlinear regime. The onset of instabilities was studied by van Duijn et al. [3] who investigated the stability of the saline boundary layer with both infinitesimally small perturbations (method of linearized stability) and finite perturbations (variational method). The linear stability analysis results in a Rayleigh number above which the system is certainly unstable, whereas the variational approach gives an upper bound to the Rayleigh number below which the boundary layer is definitely stable. The calculated bounds are separated by a region in which subcritical instabilities may be encountered, depending upon the nature of finite perturbations present [4]. Note, however, that neither the solutions, the number of solutions, nor their stability can be inferred from the approach followed by van Duijn et al. [3]. To study finite-amplitude behavior for arbitrary parameter settings, both laboratory experiments were performed [1,5] and direct time-integrations of the discretized model equations were carried out [1,3,5–7]. These laboratory and numerical experiments show the existence of finite-amplitude convection cells and finger-like salt distributions with a well-defined horizontal periodicity. The direct time-integrations give interesting first results concerning convective solutions, but an extensive survey of the influence of the initial conditions and parameter settings on the model results is virtually impossible. The aim of this paper is to investigate the existence of finite amplitude solutions (both the number of solutions and their structure), their stability, and their dependency on the system parameters by means of bifurcation analysis and continuation techniques. For this purpose, we numerically construct a low-order dynamical system by approximating the original nonlinear problem by a low-dimensional system of ordinary differential equations using a numerical Galerkin projection onto an appropriately chosen set of basis functions [8–10]. The resulting nonlinear set of amplitude equations, the Landau equations, can then be analyzed in a systematic way with standard tools from dynamical systems theory. For this purpose we introduce, next to the Rayleigh number, the Péclet number Pe. It is defined as the ratio of the salt lake underground layer depth and the saline boundary layer intrinsic length scale. We investigate the dependence of the solutions on the parameters Ra and Pe, their stability, the number of solutions and their spatial structure and characteristic length scale. Further, the sensitivity of the underlying dynamics with respect to initial conditions is explored relatively fast by time-integration of the reduced model. The paper is organized as follows. We start in Section 2 with the formulation of the mathematical model for Wooding’s salt lake problem. In Section 3 we introduce a general framework for the construction of a low-dimensional model, the so-called reduced model. For the reduced model we discuss the weakly nonlinear limit. In Section 4 we apply a pseudo-arclength continuation technique to the reduced model and present some bifurcation diagrams for relevant Péclet and Rayleigh parameter ranges. Further, using these bifurcation diagrams, (pointwise) convergence of the Galerkin projection method is shown. The temporal behavior of both the full model and the reduced model is presented in Section 5. In Section 6 we explain the stability properties of the equilibria by discussing the underlying physical (in)stability mechanisms. We conclude with a general discussion in Section 7, where the obtained results are put in perspective with laboratory experiments of Wooding et al. [1].

2. The salt lake problem We consider a uniform isotropic porous medium occupying the domain of depth H and length L

  1 1 Ω ? := (x, z ) : − L < x < L, 0 < z < H , 2

(1)

2

where z points vertically downwards. The medium is assumed to be fully saturated with a fluid of variable density ρ : i.e. water with dissolved salt. Following Wooding et al. [1,5], we prescribe along the surface {z = 0} the density and fluid flow that correspond to a ‘dry lake’ bed. That is, along the upper boundary we prescribe a constant outflow E > 0 of groundwater with dissolved salt and a constant density ρT , where ρT represents the maximum density of the salt-saturated solution. Along the lower boundary {z = H } we prescribe the density of fresh ground water, which is denoted by ρB . The fluid flow at the lower boundary is, by continuity, induced by the fluid flow at the surface boundary. Along the lateral boundaries of Ω ? we prescribe no-flow and no-flux boundary conditions. The flow equations in terms of the Boussinesq approximation are given by the continuity equation for incompressible fluids

∇ · q = 0,

(2)

the equation for fluid flow in porous media (Darcy’s law)

µ q + ∇ p − ρ gez = 0, κ

(3)

and the convection-diffusion equation for salt transport

φ ∂t ρ + q · ∇ρ = ∇ · (D∇ρ).

(4)

Here q denotes fluid discharge, µ (constant) fluid viscosity, κ (constant) medium permeability, p fluid pressure, g gravity constant, ρ density of the fluid, φ porosity and D an appropriately defined (constant) dispersivity or diffusivity. Further, ez denotes the unit vector in z-direction, pointing downwards. We recast the problem in dimensionless form by introducing

∆ρ g κ , ∆ρ = ρT − ρB , µ t {x, z } , t := . {x, z } := H τ

τ=

uc =

φH E

, (5)

The quantity uc in definition (5) denotes the scale for gravitational convective flow, i.e. the maximal shear flow between fluids with densities ρT and ρB [11]. The time scale τ represents, apart from a factor φ , the time needed to travel a distance H with evaporation speed E. Using (5), we introduce the scaled quantities S :=

ρ − ρB , ∆ρ

U :=

q uc

,

P :=

p − ρB gHz

∆ρ gH

,

(6)

where 0 6 S 6 1 represents the saturation. The dimensionless quantity P represents departures of the dimensionless pressure from hydrostatic conditions. Substitution yields

 ∇ · U = 0, (a) U + ∇ P − Sez = 0, (b) ∂ S + Ra U · ∇ S = Pe−1 ∇ 2 S , t

(7) (c)

in

  1 1 Ω := (x, z ) : − ` < x < `, 0 < z < 1 , 2

2

for t > 0, ` := L/H, subject to the horizontal boundary conditions U · ez |z =0,1 =

−1 Ra

,

S |z =0 = 1,

S |z =1 = 0,

(8a)

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

and lateral boundary conditions U · ex |x=± 1 ` = 0, 2

3.1. Choice of the Galerkin projection basis

∂x S |x=± 1 ` = 0.

(8b)

2

The initial condition is rescaled accordingly. The (Neumann) boundary conditions for P follow directly from (8), and (7)(b). The system Rayleigh number Ra and the Péclet number Pe which appear in (7)(c) are defined by Ra =

uc ∆ρ g κ = , µE E

Pe =

HE

D

3077

.

(9)

The Rayleigh number measures the ratio of the gravitational convective flow speed and the upward evaporation. Hence, keeping the evaporation rate fixed, an increase of the Rayleigh number implies that gravitational convection due to density differences becomes more important. Assuming the saline boundary layer intrinsic thickness D/E to be constant, an increase of the Péclet number implies an increase of the depth H. Therefore, the Péclet number can be seen as a ‘‘depth’’ parameter. Both the Rayleigh and Péclet number play a crucial role in the analysis of the model equations. Remark 1. Instead of using H, van Duijn et al. [3] used the ratio D/E as length scale, implying τ = φ D/E2 and Pe ≡ 1. However, using the transformations z = Pe−1 ζ , x = Pe−1 ξ , and t = Pe−1 τ , and letting Pe → ∞, it can be shown that the transformed problem is similar to the one discussed in [3]. Redefining the system Rayleigh number by Ra := Ra Pe and letting E & 0, H = 1 such that Pe & 0, we obtain the Lapwood problem.

A key feature in our approach is now to use the set of eigenfunctions corresponding to the ground state solution as the Galerkin projection basis. In this paper, we consider the ground state solution implied by the constant boundary data (8a) and (8b). This solution is uniform in the x-direction and is characterized by a uniform upflow U0 := −Ra−1 ez

in Ω ,

(14)

and a one-dimensional steady boundary layer, adjacent to the salt lake surface z = 0: S0 (z ) =

ePe(1−z ) − 1 ePe − 1

,

0 < z < 1.

(15)

We assume that all eigenvalues {σ (i) }, i > 1, calculated with respect to this ground state are simple and that the corresponding eigenfunctions {φ (i) } form a complete set. Hence

φ(x, z , t ) =

X

A(i) (t ) φ (i) (x, z ),

(16)

i>1

where i indexes the eigenfunctions in the lateral z-direction and where A(i) (t ) denote time-dependent amplitudes. Since Φ0 = Φ0 (z ), the linear operator J contains coefficients which do not depend on the x-coordinate. Therefore, taking into account the (Neumann) boundary conditions (8b), it is easily shown that φ (i) (x, z ) can be expanded in Fourier components in the x-direction with each horizontal Fourier component given by

φj(i) (x, z ) := φˆ j(i) (z ) cos(kj x), (17)  > (i) (i) (i) with φˆ j := pˆ j , sˆj and where the horizontal wavenumber kj (i)

3. Reduction to amplitude equations

is given by jπ /` for j > 0. The eigenfunction φˆ j (z ) is obtained by solving

By eliminating U in (7), we can recast the flow problem in the compact form

Jkj φˆ j

M ∂t Φ = L Φ + N (Φ, Φ),

numerically, with Jkj = J |∂ 2 →−k2 ,∂z →d/dz , cf. Eq. (13). Using (16)

(10)

where the solution vector Φ = (P , S ) , M denotes the (singular) mass-matrix, L a linear elliptic and N a nonlinear quadratic operator. Explicit expressions for M , L, and N are given in Appendix A. Let Φ0 := (P0 , S0 )> denote a steady-state solution of (10), i.e. >

(i)

= σj(i) M φˆ j(i)

(18) x

and (17), φ can be expanded as

φ(x, z , t ) =

XX

(i)

j

(i)

Aj (t )φˆ j (z ) cos(kj x),

(19)

i>1 j>0

where j indexes the modes in the horizontal plane. (i) In the theory of linearized stability, the amplitudes Aj (t ) in (19) (i)

L Φ0 + N (Φ0 , Φ0 ) = 0.

(i) (i) σ can be explicitly written in the form Aj (t ) = A¯ j e j

Next we write the solution vector Φ as a linear combination of the steady-state profile with respect to some (not necessarily small) perturbation φ = φ(x, z , t ), i.e.

constant amplitudes A¯ j are determined by the initial perturbation. For fixed Pe and `, the system is linearly unstable if and only (1) if σj > 0 for some Rayleigh number Ra and wavenumber kj .

Φ = Φ0 + φ,

The smallest value for Ra and corresponding kj for which σj = 0 are denoted by Rcrit and kcrit respectively and correspond to neutral stability. For a detailed description on how to compute Rcrit efficiently, we refer to Appendix B. In Fig. 1 we computed for Péclet numbers up to 40 the neutral stability curves R1 (kj ; Pe). The bold envelope curve represents Rcrit as function of Pe. For 0 6 Pe < 8.71, approximately, the critical wavenumbers are given by mode number j = 1, i.e. kcrit = π /`. For 8.71 < Pe < 15, approximately, we find that the critical wavenumbers now have mode number j = 2, i.e. the j = 1 mode is no longer the critical mode. The special case Pe = 8.71 is called a double-zero (DZ) codimensiontwo bifurcation [12]. At this point, both mode j = 1 and mode j = 2 become unstable. Further increasing Pe gives a cascade in which structures with larger wavenumbers take over the smaller ones, at least in the sense of linearized stability. The corresponding values

(12)

where the Jacobian J is defined by J := L + N (·, Φ0 ) + N (Φ0 , ·). Linearization of (12) gives the generalized eigenvalue problem J φ = σ M φ.

, where the

(i)

(1)

(11)

where Φ satisfies (10) and the corresponding initial and boundary conditions (8a) and (8b). Substituting (11) into (10) yields the following equation for φ M ∂t φ = J φ + N (φ, φ),

t

(13)

3078

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

(i)

(i)

and hence the Galerkin projection in which ϕj = φj is nonorthogonal. As illustrated in [14], non-orthogonal projections may result in wrong dynamical behavior in the weakly nonlinear limit. This will be demonstrated numerically in Section 4. Pieters and van Duijn [4], however, show that (after normalization)

hM φj(k) , ωφj(l) i = δkl ,

(22)

where ω = ω(z ) is given by e . In the sequel of this paper we (i) (i) redefine ϕj := ωφj to be the test functions for the Galerkin projection, unless stated otherwise. Due to relation (22) we have {α, λ, ρ}(j i,k) = 0 for i 6= k. Furthermore, to simplify notation we define Pe z

{α, λ, ρ}(j i) := {α, λ, ρ}(j i,i) .

Fig. 1. Neutral stability curves as function of Pe (grey curves) for the Fourier modes jπ/`, j = 0, . . . , 6 and ` = 23 . The bold envelope curves depict the critical Rayleigh number Rcrit . The corresponding critical wavenumbers divided by Pe are depicted as dots (•).

for the critical wavenumber kcrit divided by the Péclet number are depicted as dots. This is done to make a smooth connection to earlier work of van Duijn et al. [3]: Let k¯ crit be the dimensional critical wavenumber. Then kcrit

=

k¯ crit H

D

= k¯ crit

,

Pe Pe E i.e. in their work the critical wavenumber is measured in terms of the lengthscale D/E, and not in terms of the length scale H as in this paper. In the limit, for Pe → ∞, we find that Rcrit → 14.35 and kcrit /Pe → 0.759, i.e. linear stability results become invariant for large Pe numbers or, alternatively, for layers of large vertical extent. These numbers were also found by van Duijn et al. [3]. Since we want to investigate the behavior of the nonlinear equations (10) at least for neutral conditions, it is natural to take the eigenfunctions corresponding to Rcrit and kcrit as a projection basis [13]: Substitution of Rcrit in eigenvalue problem (17) provides (i) us the projection basis {φj }. Equations describing the time evolution of the coefficients {A(j i) (t )} in expansion (19) are determined via the Galerkin projection method. With this approach, the nonlinear partial differential equations reduce to a set of nonlinear ordinary differential equations which can be solved quickly by standard time-integration techniques. Substitution of (19) into (12), taking (i) the inner product with an as yet unspecified test function ϕj , integrating the result over the entire flow domain Ω , and using the orthogonality of the Fourier-cosines we can write for (j, i) ∈ 0, n − 1 × 1, m: m X

αj(i,l)

l =1

(l)

dAj

dt

=

m  X

m X n −1 X l,k=1 r =0

(i,k,l) (l) Ra γj,r ,j−r A(rk) Aj−r ,

(20)

αj(i,l) = hM φj(l) , ϕ(j i) i,

(21a)

λ(j i,l) + Ra ρj(i,l) = hJkj φj(l) , ϕ(j i) i,

(21b)

(l)

(i)

Ra γj,r ,j−r = hN (φr(k) , φj−r ), ϕj i.

(21c)

In Eq. (21), h·, ·i denotes the L (Ω ; R ) scalar product. A simple computation shows that (after normalization) 2

hM φj(k) , φj(l) i 6= δkl ,

Exploiting the weakly nonlinear theory as discussed in Newell and Whitehead [16] and Fujimura [17], we can derive the Stuart–Landau amplitude equation from Eq. (20). We consider the problem for Ra = Rcrit (1 + r ε 2 ) with ε  1. The parameter r = O(1) allows for the analysis of subcritical (r < 0) and supercritical (r > 0) bifurcations. For fixed Péclet number, let mode j be defined as the first unstable mode. Thus j is a given and fixed number. By a simple Taylor expansion we find that

σj(1) (Ra) ' σj(1) (Rcrit ) +r ε 2 Rcrit | {z }

2

∂σj(1) (Rcrit ) ∂R

.

≡0

(1)

For the least stable mode, Aj , we write down the amplitude equation (see Eq. (20)) in which we only include the interaction (i) (i) terms with the 0 and 2j modes, i.e. with A0 and A2j , for i = 1, m. It will turn out that only mode numbers 0 and 2j are dynamically important. It must be emphasized that the Rayleigh numbers at which these two modes become linearly unstable are relatively far away from Rcrit . Thus in the weakly nonlinear regime, they are linearly stable. We obtain

αj(1)

(1)

dAj

dt

  = λ(j 1) + Ra ρj(1) A(j 1) ( m   X + Ra γj(,01,,ji,1) + γj(,j1,,01,i) A(j 1) A(0i) i=1

) m   X (1,i,1) (1,1,i) (1) (i) + γj,2j,j + γj,j,2j Aj A2j + · · · .

(23)

i=1

where

(i,k,l)

3.2. The weakly nonlinear limit of the reduced model

 λ(j i,l) + Ra ρj(i,l) A(j l)

l =1

+

The coefficients in (21) are computed numerically with a standard Gauss–Chebyshev quadrature rule, see Pieters [15] for more details.

The approach followed in (23) differs from the standard weakly nonlinear theory in the sense that we take into account m eigenfunctions in the z-direction. This has to be done since for the non-self-adjoint eigenfunctions and the nonlinearity we are considering, the behavior of the solution at the bifurcation point (e.g. magnitude, and even the sign of Λ1 in Eq. (25)) depends on the projection of the first eigenfunction on all other eigenfunctions in the z-direction. Rescaling the time and amplitudes to get the significant (1) (1) degeneration, it is found that τ = ε 2 t and Aj = ε A˜ j , and the (i)

(i)

amplitudes A0 and A2j must be of order ε 2 . Using these scalings in Eq. (20) and using the relation

  (i) (i) σj (Rcrit ) ≡ λj + Rcrit ρj αj(i) (i)

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

we obtain

ε2 ε

2

(i)

dA˜ 0 dτ

(i) dA˜ 2j

dτ

3.3. The reduced model for general Rayleigh numbers

= σ0(i) A˜ (0i) +

γ0(,ij,,1j,1) Rcrit (i)



(1)

2

+ O(ε2 ),

γ2j(i,,j1,j,1) 

(1)

2

+ O(ε2 ).

A˜ j

α0

(i) ˜ (i)

= σ2j A2j + Rcrit

A˜ j

α2j(i)

(i)

(i)

For ε & 0, we can compute the slaved amplitudes A˜ 0 and A˜ 2j from (1)

A˜ j :

γ{(0i,,12j,},1)j,j

(i)



A˜ {0,2j} ' −Rcrit (i)

α{0,2j} σ{(0i),2j}

(1)

A˜ j

2

.

(1)

dτ

= Rcrit

∂σj(1) (Rcrit )



∂R



(1)

r + Λ1 A˜ j

2 

(1)

A˜ j .

(24)

Λ1 = Rcrit

! −1

m X

∂R

(i)

Λ1 ,

(25)

i=1

with (i)

Λ1

  −γ (i,1,1) 0,j,j (1,i,1) (1,1,i) = γj,0,j + γj,j,0 (i) (i) (1) α0 σ0 αj  −γ (i,1,1)  2j,j,j 1,i,1) (1,1,i) . + γj(,2j + γ ,j j,j,2j α2j(i) σ2j(i) αj(1) (1)

Steady solutions of Eq. (24) are given by A˜ j



˜ (1)

solution) and r + Λ1 Aj

2

(1)

 Ra = Rcrit



(1)

1 − Λ1 Aj

2 

.

= 0 (trivial

= 0 (nontrivial solutions branching

from a pitchfork bifurcation). Since Aj Rcrit )/Rcrit , we find

Model reduction is based on the assumption that the series in Eq. (20) can be truncated in such a way that the essential nonlinear dynamics is adequately reproduced in given relevant Ra and Pe ranges. Unfortunately, there exists no universal truncation rule and to test the actual convergence is in general a tedious task. We distinguish two convergence types: (i) Convergence of solutions within the reduced model framework; (ii) Convergence of the reduced model solutions to solutions of the original nonlinear problem.

In Eq. (24), the Landau coefficient Λ1 is given by

∂σj(1) (Rcrit )

For moderate and large Rayleigh numbers the weakly nonlinear approach cannot be used. To explore the steady-state solutions for these Rayleigh numbers, a pseudo-arclength continuation method ([18]) is applied to the coupled system of equations (20) with (j, i) ∈ 0, n − 1 × 1, m. Hence a system of n × m coupled ordinary differential equations is investigated to find nontrivial steady state solutions of the Wooding problem in a specific part of the (Ra, Pe)parameter space. The linear stability of the steady state solutions is investigated by calculating the eigenvalues of the corresponding Jacobian. 3.4. Convergence issues

Applying the scalings and using these expressions in Eq. (23), the (1) amplitude equation for A˜ j then reads dA˜ j

3079

= ε A˜ (j 1) and r ε 2 = (Ra −

(26)

The sign of Λ1 in Eq. (26) determines the type of bifurcation at threshold: Negative Λ1 corresponds to supercritical bifurcations, positive Λ1 with subcritical bifurcations. For Λ1 = 0, the weakly nonlinear approach with only three horizontal modes (n = 3) degenerates and one needs to incorporate higher order horizontal modes.

Convergence of type (i) heavily depends on the projection method used or, alternatively, depends on the choice of the test (i) functions ϕj . This will be explored in the model results presented in Section 4. Non-self-adjointness of the linearized problem affects the convergence of type (ii). As the Wooding problem is not selfadjoint, the parameter range for Pe and Ra in which the reduced model captures the nonlinear dynamics of the original problem may be very small. To be more precisely, nonlinear interactions between the non-orthogonal basis functions that are essential in the low-order model to reproduce the dynamic behavior observed in the full model, may not be necessarily restricted to the basis functions used to construct the low-order model. As a consequence, convergence tests by successively increasing the number of basis functions can be very misleading. In van der Vaart et al. [9] for example, Table II shows that the accuracy of the approximation can decrease when one incorporates additional basis functions. Therefore, in view of convergence type (ii), it is essential to compare the results from the reduced model with those of the full set of equations. In this way it will be shown that both the steady state solutions and the (temporal) dynamics of the full model is captured for a large range of parameter values by the low-order model. This will be further discussed in Section 5, where a direct numerical simulation is performed to the fully nonlinear equations (10) by means of a finite element method (fem). 4. Model results

Remark 2. At double-zero bifurcation points we have not a single mode that becomes unstable, but instead we have two modes that become unstable simultaneously. In particular, in the neighborhood of these DZ points, the two modes will strongly interact and hence the assumed scaling of the amplitudes is not valid anymore. One can show that the magnitude of the two modes must be of order ε 2 , whereas all other modes must have amplitudes of order ε 3 or smaller. The scalings for time and the Rayleigh number remain unchanged. A detailed analysis of the Landau coefficient for Rayleigh numbers in these DZ neighborhoods is beyond the scope of this paper.

In their pioneering papers, Wooding et al. [1,5] considered both numerical and Hele–Shaw cell experiments. In this paper, we focus our attention on the cell experiments. For this purpose we distilled the parameters of their experimental setup in Table 1. From the table we compute ` = 23 ; this value will be fixed for all experiments in this paper. Although Pe = 166.67 in the experiments, it turns out that our approach is not capable to deal with such a large Péclet number, and therefore we limit our investigations to Péclet numbers up to 10 and show some results for Pe = 40. As will be discussed in Section 7, laboratory Hele–Shaw experiments [5] can be partly explained from the

3080

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

that in the x-direction only 3 modes are considered. We consider here two distinct cases. First we compute the Landau coefficient by (i) (i) means of the Galerkin method in which ϕj = ωφj . These curves are depicted in black. In the second case we compute Λ1 by means (i) (i) of the Galerkin method in which ϕj = φj . Since for this case we lack diagonality of the coefficients (21a) and (21b), we need (albeit not mentioned in this paper) a much more complicated version of expression (24). These Λ1 curves are depicted in grey. For Pe numbers up to ∼20 both approaches converge to the same Λ1 , but for larger Pe numbers the grey curves diverge. We may conclude that the weighting function ω critically stabilizes the convergence of the method. In the limiting case Pe & 0, both approaches are similar and self-adjoint. In particular, we have (1)

(1,1,1)

Λ1 = 0, (2)

Λ1 = − (i)

Λ1 = 0,

Fig. 2. The Landau coefficient Λ1 as function of m and n = 3, for several Pe values. Comparison between computation with (black curves) and without (grey curves) weighting function ω in the Galerkin projection. Table 1 Hele–Shaw cell parameters for Wooding’s experimental setup Cell depth H, m Cell evaporation length L, m Cell intrinsic permeability κ , m2 Density difference 1ρ at 20 ◦ C, kg/m3 Dynamic viscosity µ at 20 ◦ C, kg/m s Diffusion coefficient D at 20 ◦ C, m2 /s Hele–Shaw Péclet number Pe Rayleigh number Ra

0.75 × 10−1 0.5 × 10−1 3.33 × 10−9 16.8 1.103 × 10−3 0.9 × 10−9 166.67 2.3 – 40.1

Table 2 Numerical continuation experiments Experiment

Pe

Ra

I II III

{3, 5, 8.72, 10}

[0, 20] [0, 65] [0, 15]

10 40

nonlinear behavior at Pe = 40. A possible remedy to deal with larger Péclet numbers will be presented there as well. Three experiments will be performed in which Pe is fixed and Ra varied to investigate the solutions of system (20); see Table 2 for a summary. The various steady-state solutions are presented in bifurcation diagrams of the amplitude of the first unstable mode (solid branches) and its slaved amplitudes (dashed branches). The stability of these branches is determined by evaluating the eigenvalues of the Jacobian of system (20). We use the convention that stable branches are depicted as black curves, and unstable ones as grey curves. The results of the first two experiments are discussed in this section, the third experiment will be discussed in Section 6. 4.1. Weakly nonlinear results In Fig. 2 the convergence of the Landau coefficient Λ1 versus the number of eigenfunctions in the z-direction, m, is depicted. Recall

since γj,r ,j−r ≡ 0 for j, r ∈ 0, n − 1,

π2 4

 1+

1

`2

 =−

13 16

π 2,

for i > 2,

which implies that Λ1 ' −8.02. This value is recovered numerically in Fig. 2. (i) For Pe > 0, however, Λ1 6= 0 for i > 1. Hence the weakly nonlinear analysis cannot be limited to the incorporation of just (i) one depth mode (i = 2). Higher order contributions of Λ1 to Λ1 indicates the non-self-adjointness of the Wooding problem. More (i,1,1) (i,1,1) (i) generally, the coupling coefficients γ0,j,j and γ2j,j,j in Λ1 could serve as a measure for the non-self-adjointness of the underlying problem: The larger Pe the more non-self-adjoint the problem becomes, and hence the more depth eigenmodes m we have to incorporate for convergence. 4.2. Bifurcating structures appearing in Experiment I Convergence tests show that for the Pe- and Ra-ranges in Experiment I we have sufficient convergence for n = 10 and m = 10. These settings will be used throughout this subsection. The results for Pe ∈ {3, 5} are presented in Fig. 3 and for Pe ∈ {8.72, 10} in Fig. 4. According to Fig. 1, mode j = 1 is the first unstable mode for Pe = 3. Its critical Rayleigh number is RaP1 ' 18.75. The transition at criticality follows a stable supercritical pitchfork bifurcation. From this point, two stable steady-state solution branches appear which have a 12 -finger structure. Observe that the lower branch solutions are mirrored versions of the upper branch solutions. This behavior is well-known in the weakly nonlinear theory of Rayleigh–Bénard convection (i.e. Pe & 0). For Pe = 5 we find a similar behavior, but the transition point at criticality now follows an unstable subcritical pitchfork bifurcation at RaP2 ' 14.82. This is in agreement with Fig. 2 since for this Péclet number Λ1 has a positive sign, indicating a subcritical pitchfork bifurcation. At RaF1,2 ' 14.64 both the upper and lower branch become stable by a fold bifurcation. The patterns still have a 1 -finger structure but the fingers are ‘‘deeper’’. 2 In Fig. 3 the weakly-nonlinear approximation Eq. (26) is plotted as well for 0 6 ε < 1. They are depicted as slim dashed lines. For Pe = 8.72, the bifurcation diagram in Fig. 4 is more complicated since now a double-zero bifurcation occurs: both modes j = 1 and j = 2 become unstable at RaDZ = 15.35, see Fig. 1. This type of bifurcation degenerates the bifurcation diagram in the sense that both mode j = 1 and j = 2 form primary branches. To get a clear picture of the bifurcation structure, the amplitudes of both the j = 2 and j = 1 modes are shown in Fig. 4(a, b), respectively. They will be discussed separately:

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

3081

Fig. 3. Bifurcation diagram for Pe = 3 and 5, with n = 10 and m = 10. Here F denotes a fold bifurcation and P a pitchfork bifurcation. Weakly nonlinear amplitudes are depicted by slim dashed branches. The patterns give an impression of all other solutions on that branch. ‘‘s’’ denotes stable solutions.

(1)

(1)

– On the A2 > 0 branch for which A1 ≡ 0 we observe an unstable equilibrium with two real positive eigenvalues. At the fold bifurcation RaF3 ' 15.33, the equilibrium solution is still unstable with only one positive real eigenvalue. On this branch no further stable solutions exist for Ra 6 20. The unstable patterns all exhibit a symmetric structure which can be seen as two 12 -finger structures glued together at x = 0. (1)

(1)

– For A2 < 0 and A1 ≡ 0, the equilibrium changes stability at the fold bifurcation RaF4 ' 15.33, i.e. the number of eigenvalues with positive real part changes from one to zero. The solutions on this branch lose stability at RaP4 ' 15.70 through a pitchfork bifurcation. Following this branch even further, the solutions become stable at RaP5 ' 19.11, which is again a pitchfork bifurcation. The steady-solutions on the lower branch all have a symmetrical 1-finger structure. – At DZ two branches with finite amplitudes for both modes j = 1 and j = 2 exist with Ra > Rcrit . Their corresponding steadystate solutions have one unstable eigenvalues. They become more unstable at the fold bifurcations RaF7,8 ' 16.22. At RaH1,2 ' 15.93, the unstable equilibria become stable through a Hopf bifurcation. Following these branches, the amplitude of (1) A1 (mode j = 1) gradually decreases to zero and it connects to the solid branch at the pitchfork bifurcation P4, see also Fig. 4(a). The Hopf bifurcation is supercritical, i.e. stable timeperiodic solutions exist for Rayleigh numbers larger than RaH1,2 . The temporal behavior in the vicinity of these Hopf bifurcations will be investigated in Section 5. The steady-state solutions on these two branches are asymmetric since they can be seen as a mixture of mode j = 1 and mode j = 2. – Two other branches exist for Ra < Rcrit . At RaF5,6 ' 12.42 they become stable by two fold bifurcations, resulting in steadysolutions with a 21 -finger structure. Observe again that there exist stable solutions for Ra < Rcrit .

Concerning the number of stable solutions as a function of Ra, we can conclude the following: Up to RaF5,6 we have one stable solution, the trivial state. For Rayleigh numbers between RaF5,6 and RaP4 there exist three stable solutions, one trivial solution, and two solutions with a 12 -finger structure. For Rayleigh numbers between RaF4 and RaDZ three stable solutions exist, two 12 -finger structures and one 1-finger structures. From RaP4 to RaH1,2 we have four solutions (two 12 -finger structures and two asymmetric 1-finger

structures). From RaH1,2 to RaP5 we only have two stable 12 -finger structures. Furthermore, (at least) close to RaH1,2 stable periodic solutions exist. Finally from RaP5 up to Ra = 20 we have three stable steady-state solutions, two 21 -finger structures and one 1finger structure. Further increasing the Péclet number to Pe = 10 gives a similar picture, but now the mode j = 1 branch, bifurcating at RaP9 ' 14.7, is detached from the branch bifurcating at RaP7 ' 16.34, see Fig. 4(c, d). The pitchfork P7 is connected via two fold bifurcations F13,14 and two Hopf bifurcations at RaH3,4 ' 15.35 to pitchfork P10, as in the Pe = 8.72 case, see Fig. 4(d). 4.3. Bifurcation structures appearing in Experiment II The results presented in Fig. 4 for Pe = 10 are for ‘‘moderate’’ Rayleigh numbers. It is expected that for large Rayleigh numbers, one needs more horizontal Fourier modes n and more eigenmodes m for convergence. Unfortunately, these numbers can only be found by means of convergence tests. Convergence up to Ra ≈ 62 is reached for n = 100 and m = 20. The converged bifurcation diagram is presented in Fig. 5. It is the continuation of Fig. 4(c) for larger Rayleigh numbers. Concerning (1) the solid upper branch for which A1 ≡ 0, a pitchfork bifurcation P12 occurs at Ra ' 22.1, after which the steady-state solutions are stable.

3082

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

Fig. 4. Bifurcation diagram for Pe = 8.7 and 10, with n = 10 and m = 10. Here F denotes a fold bifurcation, P a pitchfork bifurcation, and H a Hopf bifurcation. Top figures correspond to the critical mode. The bottom figures correspond to the slaved mode that was critical in Fig. 3. The patterns give an impression of all other solutions on that branch. ‘‘s’’ denotes stable and ‘‘u’’ unstable solutions.

The mode j = 1 branch (dashed curve) bifurcating at P9 gives two new stable Hopf bifurcations at RaH5,6 ' 58.5. Following the mode j = 1 branch at P12 (dashed curve) results in two new unstable Hopf bifurcations at RaH7,8 ' 60.2. The two branches starting from P9 and P12 converge to each other, but they do not connect for Ra 6 62. Repeating the continuation for the branch

bifurcating at P11 (dashed curve) gives four new unstable Hopf bifurcations at H9,10 and H11,12 for Ra up to 62. To conclude, we note here that the number of stable solutions from Ra = 20 up to RaH5,6 is four: two of them have a 12 -finger structure, and the other two a 1-finger structure. In addition, for the subcritical solutions we can conclude the following: the larger

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

Fig. 5. Bifurcation diagram for Pe = 10 with n = 100 and m = 20.

the Péclet number, the smaller the Rayleigh number for which there exists stable subcritical solutions (see also Fig. 9 for Pe = 40). This corresponds with results obtained for Pe → ∞: a lower bound for the existence of stable subcritical steady-solutions is given by Ra = 8.59 [3]. This lower bound is found using variational arguments in which an optimum is found over the space of L2 -integrable solutions satisfying a pointwise differential constraint.

3083

For the direct numerical simulations we use the same parameters and initial conditions as for the reduced model approach. The reduced model (20) is integrated in time by means of a Runge–Kutta method. We consider three experiments for Pe = 10: One initial condition is given by an unstable 1-finger steadystate at Ra = 20 which evolves to a stable 12 -finger solution which is located on a different branch. For the other two experiments the initial conditions are given by the amplitudes at RaH3 = 15.35 and RaH5 = 62.5 respectively. These amplitudes correspond to the unstable equilibria which will evolve to stable Hopf bifurcations. Fig. 6 shows a comparison between the two methods in terms (i) (i) of Aˆ j (t ) (solid curves) and Aj (t ) (dashed curves) for some i and j. As can be seen from this figure, the amplitudes match qualitatively well. However, for large times the finite element (i) projected amplitudes Aˆ j in Fig. 6(a) deviate from the reduced model ones in the sense that mode j = 0 dissipates faster than higher order modes. For the Hopf bifurcations in Fig. 6(b, c) we observe a tiny shift in the frequency of both Hopf bifurcations and we observe in Fig. 6(c) a small secondary modulation, which indicates the fact that the Hopf bifurcation itself becomes unstable. (i) (i) The difference between Aj and Aˆ j might be explained from the numerical methods used. In the finite element approach used for this experiment, we implicitly introduce two numerical artifacts to the nonlinear dynamics of the underlying problem. Firstly, due to the upwind discretization we introduce artificial numerical diffusion to the problem. Hence, the actual Rayleigh and Péclet numbers in the discretization do not correspond to the physical ones, i.e. the ones used in the reduced model. This seriously affects the comparison with the reduced model solutions. Secondly, the numerical scheme decouples the nonlinear parabolic part from the elliptic part and uses Picard iterations to capture the nonlinear dynamics. In the elliptic part, a derivative from the previous timestep has to be taken in the right-hand side. Even when one takes small timesteps, such an iteration adds up numerical errors during the time-integration.

5. Temporal dynamic behavior

6. Transition mechanism

The reduced model approach is capable of revealing the dynamics of nontrivial solutions which are (in)directly connected to the ground state, i.e. the zero branch. The next question that one could ask is: What will be the temporal dynamics of the system when an (arbitrary) initial state is prescribed. It is very likely that the ‘‘path’’ in the phase-plane, which will eventually lead to some steady-state solution, will not be necessarily spanned by the modes that build up the steady-state as in the reduced model. In view of the convergence issues discussed in Section 3, we will compare some time-integration results of the reduced model with direct numerical simulations of the full model. Generally, direct numerical simulations may not be sufficiently accurate to capture subtle nonlinear dynamics due to e.g. the existence of additional non-physical diffusion and round-off errors. This implies that unstable solutions are hard to find with this approach, whereas knowledge of these solutions is often necessary to understand transitions, as well as the dynamical mechanisms underlying the physics of transition. Let Sfem (x, z , t ) be the solution of the fem direct numerical simulation as discussed in [3]. To compare this solution to the reduced model time-integration results, we project the deviation from the ground state solution, (Sfem − S0 ), onto the eigenfunctions given by Eq. (17), i.e.

In the previous sections the stability of the equilibrium solutions was determined using the information contained in the eigenvalues. Now we place the mathematically obtained result in a physical context. To that end the equilibrium solution

(i) S (x, z , t ) − S0 (z )) sj (x, z )ω(z ) dxdz (i) Ω ( fem Aˆ j (t ) := ,   2 R (i) s ( x , z )ω( z ) dxdz j Ω

R

(i)

(i)

where sj denotes the second component of the eigenfunction φj .

e Φ0 := Φ0 + φ eq ,

e Φ0 = (e P0 , e S0 )> ,

is perturbed by the least stable eigenmode e φ := (e p,e s)> of e Φ0 . Here Φ0 denotes the one-dimensional ground state solution and φ eq is some steady solution of Eq. (12). The perturbed salinity flux is given by F := −Ra (−∇e P0 + e S0 ez )e s + (−∇e p +e sez )e S0 + Pe−1 ∇e s,




(27)

where e p can be expressed in terms of e s using the linear relation given in Eq. (A.1). The (de-)stabilizing mechanisms can be identified by studying the various components of the salinity flux: pe S0 , FcX = −Ra −∂xe P0e s + ∂xe





FdX = Pe−1 ∂xe s, FcZ = −Ra (−∂ze P0 + e S0 )e s + (−∂ze p +e s)e S0 ,





(28)

FdZ = Pe−1 ∂ze s,

where the superscripts X and Z denote the horizontal and vertical direction respectively, and c (d) denotes the convective (diffusive) part. We refer to Appendix C for technical details on the numerical computation of these flux components.

3084

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

(i )

(i )

Fig. 6. Comparison of three direct numerical simulations of the full model Aˆ j (t ) (black curves) and time-integrations of the reduced model Aj (t ) (grey curves). The pair (j, i) denotes the amplitude corresponding to horizontal mode j and depth mode i. The reduced model time-integrations for (a) and (b) are performed with n = 10 and m = 10. For (c) we used n = 100 and m = 20.

Firstly we consider the case in which the pattern e Φ0 is just given by Φ0 , i.e. φ eq ≡ 0. This equilibrium is perturbed by the salinity pattern shown in Fig. 7(a). Here darker colors indicate a positive salinity disturbance, a negative salinity disturbance is indicated by lighter colors. The convective flux of salt is depicted by the arrows in Fig. 7(a; top). In Fig. 7(a; middle,bottom) the diffusive and total flux are shown respectively. We define a control volume V related to this perturbation as

 V :=

1

1

2

2

(x, z ) :e s (x, z ) > 0; ε < z < 1, − ` < x <



` ,

with ε  1. This control volume is indicated by the black lines in Fig. 7(a). Next we calculate the total diffusive and convective salt fluxes through the boundary of V , i.e. FdX

Z :=

FcX :=

Γ2 ∪Γ4

Z Γ2 ∪Γ4

FdX

,

FcX ,

FdZ

Z :=

FcZ :=

Γ1 ∪Γ3

Z Γ1 ∪Γ3

FdZ , FcZ ,

where Γ1 denotes the top horizontal, Γ2 the left vertical, Γ3 the bottom horizontal, and Γ4 the right vertical boundary of V , such that Γ1 ∪ Γ2 ∪ Γ3 ∪ Γ4 = ∂ V . If the total flux is positive, there is a net transport of salt from V , indicating stability. If the total flux is negative, salt is accumulating in the control volume and the perturbation starts to grow. In Fig. 7(b; top) the integrated horizontal (FcX ) and vertical (FcZ ) convective flux components through the boundaries are shown as a function of the Rayleigh number Ra. Fig. 7(b; middle, bottom) show these quantities for the diffusive and total flux respectively. From these figures we can conclude the following. For Ra < Rcrit the diffusive outflux is stronger than the convective influx,

i.e. the initial salt distribution of the least stable perturbation is stabilized by diffusion. At critical Rayleigh number, the convective influx balances with the diffusive outflux. Further, for Ra > Rcrit the convective influx becomes stronger than the diffusive outflux. This leads to an increase of salt in V and hence the equilibrium is unstable. When changing the Rayleigh number, it is mainly the vertical diffusive flux that increases and tries to balance the increased horizontal convective flux. Next we consider the Hopf bifurcation. This type of transition is more complicated since the perturbation exhibits time-periodic behavior. The approach is to consider the time-evolution of the pattern, in combination with the flux vector fields. In Fig. 8(b–e) four snapshots of the difference between the actual timedependent salt distribution and the time-independent equilibrium distribution are shown. In Fig. 8(b) a negative salt anomaly is present in the center of the cell, whereas the two upper corners have a weak positive salt anomaly. This pattern corresponds with the real part of the most unstable eigenfunction. This pattern is not stable and starts to develop: salt is transported from the righthand side of the cell towards the lefthand side of the cell, resulting in a negative salt anomaly at the righthand side and a positive one at the lefthand side, see Fig. 8(c). This pattern corresponds with the imaginary part of the eigenfunction. Next, salt fluxes result in the salt distribution as shown in Fig. 8(d). Compared to Fig. 8(b), it is interesting to see that the amplitude of the salt anomaly has decreased with approximately a factor 2, indicating that the underlying equilibrium is stable with respect to this perturbation. Next, the diffusive and convective processes result in the pattern shown in Fig. 8(e), after which the initial salt distribution is recovered, although the amplitude is appreciably smaller.

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

3085

˜ 0 , where darker colors indicate a positive salinity disturbance. The salinity fluxes are denoted by the arrows. Note that only the Fig. 7. (a) Least stable perturbation of Φ upper half of the cell is depicted. (b) Flux balance in a control volume V . Positive flux indicates outflow of salt, and negative flux inflow. See text for a detailed discussion.

To understand the cyclic behaviour of this perturbation, we will study the salt fluxes related to the various stages shown in Fig. 8(b–e). The salt perturbation shown in Fig. 8(b) results in a diffusive flux that is stabilizing and that the convective flux is destabilizing. The latter flux results in a divergence of salt in the center of the cell and a convergence at the two top corners. Due to the asymmetry in the perturbed salt field the two convection cells that result in the convective transport of salt are not of the same strength. The righthand side cell is slightly stronger, which can be easily inferred from the asymmetric part of the salt perturbation. This (slight) asymmetry results in two effects: Firstly, the diffusive flux results in a stronger flux of salt directed from the righthand upper corner to the center of the cell than from the lefthand side corner. Secondly, the destabilizing convective flux towards the upper righthand corner is slightly stronger than the flux towards the upper lefthand corner. It turns out that, due to the asymmetry, the total stabilizing flux from the upper righthand corner is stronger than the stabilizing flux from the upper lefthand corner. Hence the salt is more rapidly removed from the upper righthand corner than from the upper lefthand corner. As a result, the negative anomaly of salt is shifted to the right and a region with a positive anomaly is formed just left of the center as can be seen in

Fig. 8(c). Note that the salt concentration at the new maximum and minimum is much lower than in the situation in the left top corner figure. This can be understood as the formation of the resulting pattern takes place when the stabilizing diffusive fluxes dominate over the destabilizing convective fluxes. In this new situation, only a single convection cell exists. Again, the diffusive flux tries to damp the instability. Apart from a convergence–divergence pattern which almost balances the diffusive one, the convective flux results in a strong convergence of salt in the middle of the cell, approximately 20% below the surface, where a strong convergence of salt is observed. Hence the original perturbation will be damped by the diffusive fluxes, but due to the destabilizing convective fluxes a new pattern of positive and negative salt anomalies is formed, see Fig. 8(d). This pattern has a maximum and a minimum salinity larger than the one observed in Fig. 8(c), but smaller than the original perturbation (Fig. 8(b)). Hence the convective fluxes are destabilizing, but integrated over a complete cycle, the stabilizing diffusive fluxes are dominant. The development of this pattern into the one depicted in Fig. 8(e) can be understood using the same type of reasoning as for Fig. 8(b). Clearly, the above described transition mechanism involves the interaction of a salinity perturbation where the perturbed velocity

3086

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

Fig. 8. (a) Amplitudes as function of time, relative to the stable equilibrium for Ra = 15.139 and Pe = 10. The ⊗ denote snapshots at t = 0, 0.25, 0.45, and 0.8 respectively, and T denotes the length of one cycle. (b–e) Time evolution of the difference between the actual time-dependent and time-independent salt distribution. (f) Amplitudes as function of time, relative to the unstable equilibrium for Ra = 15.41 and Pe = 10.

consists of a 2-cell structure (the real part of the least stable eigenmode) and one with a perturbed velocity consisting of a 1-cell structure (the imaginary part of the least stable eigenmode), see Fig. 8(b–e) which visualize these cells by means of closed contours. The time-evolution of the amplitudes corresponding to these modes is shown in Fig. 8(a). In this figure the solid line shows the amplitude of the two-cell mode, A2−cell , the dashed line the amplitude of the one-cell mode, q A1−cell , and the thick solid line the total amplitude Atotal =

A22−cell + A21−cell . We start with a

2-cell mode in (b). Due to diffusive damping the curve decreases towards (c). From this point, the convective destabilization of the 1-cell mode takes over and moves the curve up to (d). Then the cycle repeats and the diffusive stabilization becomes stronger than the convective destabilization and the curve ends up in (e). The convective destabilization turns out to be weaker over the complete cycle than the diffusive stabilization and hence the amplitudes eventually decay to zero. If the Rayleigh number is increased, the destabilization due to the convective fluxes becomes dominant and the perturbation starts to grow. The asymmetry of the perturbation is just the other way around and hence the net convergences and divergences of the fluxes are now mirrored with respect to the {x = 0} axis, resulting in a periodic solution that ‘‘moves’’ in the opposite direction. In this case the damping due to diffusion during the first quarter of the cycle is not strong enough to prevent the amplitude of the perturbation to grow during the second part of the cycle. Hence convection dominates and the perturbation starts to grow. For the unstable case, the amplitudes of the 1- and 2-cell are shown in Fig. 8(f). Now it is seen that, compared to the diffusive stabilization, the convective destabilization turns out to

be stronger over the cycle, and hence the amplitudes of both modes grow in time, i.e. the time-periodic solutions becomes unstable. 7. Conclusions and future directions In this paper we have constructed a low-dimensional model to describe convective instabilities induced by a saline boundary layer adjacent to the horizontal surface of a porous medium. For this purpose a Galerkin projection of the nonlinear model equations onto a set of basis functions is applied. The projection basis is given by the eigenfunctions at neutral stability. Since this set of functions do not form an orthogonal basis due to the non-self-adjointness of the problem, we use the adjoint eigenfunctions in the Galerkin projection. The obtained Landau equations are then truncated in such a way that the principal nonlinear aspects of the unstable boundary layer is still preserved. This approach allows for a drastic reduction of computational cost. We have demonstrated that both the equilibrium and time-dependent behavior of the reduced model solutions are in good agreement with direct numerical simulations of the fully nonlinear model. Predictions of the length scale of the observed patterns that are based on the least stable wavenumber found by linear stability analysis, can be misleading. Hence it is not clear which of these patterns will eventually show up when doing (numerical) experiments. This is illustrated in Fig. 9, where a part of the bifurcation diagram for Pe = 40 is shown. The grey unstable curve corresponds to the primary bifurcation. According to linear stability analysis, 3-finger solutions appear at the onset. In the fully nonlinear regime, however, none of these 3-finger structures on this branch are stable. Further, stable subcritical solutions do exist but have a 12 -finger structure. A similar behavior is observed in

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

3087

To conclude, a careful analysis of the reduced model reveals a subtle competition between the vertical diffusive and horizontal convective fluxes. This competition is shown to be responsible for the observed bifurcation types. Appendix A. Definition of the differential operators We first combine (7)(b), and (7)(c) to obtain

∇ P = ∂z S . 2

(A.1)

Substitution of (7)(c) in (7)(a) gives

∂t S = Pe−1 ∇ 2 S + Ra ∇ P · ∇ S − Ra S ∂z S .

(A.2)

Eqs. (A.1) and (A.2) define the operators M , L, N and J :

 M :=



0 0



1 , 0

N (Φ1 , Φ2 ) := Ra

L :=



0 Pe−1 ∇ 2 + ∂z − Ra C ω

2



−∂z  2

∇ Φ11 − Φ21 ez · ∇ Φ2

Ra C ω−1 ∂z

∇

∇

2




 J :=

Pe−1 ∇ 2

0

,  −1

−∂z

,

,

Pe

where C := −Pe ePee −1 and with ω as defined in Section 3.1. Appendix B. Determining the critical Rayleigh number (i)

Consider again eigenvalue problem (18). The eigenvalues {σj } Fig. 9. Bifurcation diagram for Pe = 40 with n = 30 and m = 35.

the laboratory experiments of Wooding et al. [1] in which Pe = 166.67: we observe many fingers at onset, and during the course of the experiment these small structures disappear in favor of larger ones and eventually we observe a pattern that resembles the stable 21 -finger solution in Fig. 9, see also Fig. 4(a) and Fig. 5(a). To answer the question which of the patterns is most likely to pop up in an experiment, one has to determine for each stable solution its domain of attraction. These attraction domains can be determined numerically by means of stochastic methods, see for example Gardiner [19]. This will be a topic for future research. The approach followed in this paper has its limitations. For large Rayleigh numbers one has to incorporate many modes to get sufficient convergence. This is also to be expected: the modes for projection are based on the one-dimensional ground state Φ0 , and the patterns that may pop up in this large Rayleigh number regime do not necessarily have very much in common with the modes that form the projection basis. For example, for large Péclet numbers one tries to project relatively deep fingers onto the eigenfunctions in z-direction which support is very localized just below the surface. Therefore, it would be better to change the projection basis when one continues a specific branch for large Péclet and/or Rayleigh numbers. There are two approaches to obtain the new projection basis. The easiest way is to take the amplitudes of a specific pattern on the followed converged branch and to determine the eigenvectors of its Jacobian. These eigenvectors containing amplitudes are then used to recompute the coefficients (21). A major drawback of this approach is the fact that one sticks with the reduced model approximation. Furthermore, the set of eigenvectors do not necessarily span the correct pattern. To circumvent this problem, one could also plug the specific pattern into the fully nonlinear problem and perform a time-integration until the solution reaches the attractor. This direct numerical simulation approximation can then be used to derive a twodimensional eigenvalue problem as in the ground state case. In this way the number of projection modes can be drastically reduced to obtain convergence for large Rayleigh numbers.

(i)

clearly depend on Ra and Pe. Further, it can be shown that σj ∈ R for all j > 0 and i > 0. Then the first step in our approach (i) is to set σj = 0 in (18) (neutral stability condition) and solve for the set of wavenumbers {kj } a generalized eigenvalue problem for Ra = R. Here R is introduced to distinguish between R as eigenvalue and the parameter Ra for the actual physical system. By using a modified spectral Chebyshev–Galerkin method we obtain the set of eigenvalues 0 < R1 (kj ) < R2 (kj ) < R3 (kj ) < · · ·. From the neutral stability curve R1 (kj ) we determine the critical Rayleigh number beyond which the system just loses its stability, i.e. we define kcrit = arg min R1 (kj ) kj

and Rcrit := R1 (kcrit ).

Note that Rcrit depends indirectly on ` via the wavenumbers kj and on the Péclet number Pe. Appendix C. Computation of the flux components In calculating the flux components in Eqs. (28) one should realize that only their divergence n o can be directly expressed in (i)

terms of the eigenfunctions sj

. Hence to be consistent with the

Galerkin results, one should use the identity

D E D E ∇ · F , ωs(j i) = − FcX + FdX , ω∂x (s(j i) ) | {z } (i),c e(i),d e H +H j

− |

D

FcZ

+

j

E , ∂z (ωs(j i) ) . {z }

FdZ

(i),c e(i),d e V +V j

e(i),c

Using the amplitudes Hj the flux components as FcX , FdX (x, z ) =





j

(i),d

,e Hj

(i),c

,e Vj

(i),d

, and e Vj

m X n −1 n o X (i),c e(i),d (i) e H ,H H (x, z ), j

j

j

i =1 j =0

FcZ , FdZ (x, z ) =





m X n−1 n o X (i),c e(i),d (i) e V ,V V (x, z ), j

i=1 j=0

j

j

we reconstruct

3088

G.J.M. Pieters, H.M. Schuttelaars / Physica D 237 (2008) 3075–3088

with (i)

Hj (x, z ) =

x

Z

(i)

sj (ξ , z ) dξ , 0

(i)

Vj (x, z ) = −

1

Z

(i)

sj (x, ζ ) dζ . z

(i)

(i)

(i)

A simple computation shows that {Hj , ∂x (ωsj )} and {Vj , ∂z (i)

(ωsj )} are bi-orthogonal bases in L (Ω ). However, since the perturbations e s and −∂ze p +e s have zero boundary conditions at z = 0, 1, it can be immediately seen that FcZ also has zero boundary (i) conditions, whereas the projection basis {Vj } has a non-zero boundary condition at z = 0. Apparently the projection basis is 2

not very suitable for this flux component and it leads to the wellknown Gibbs phenomenon at z = 0 when approximating FcZ . A similar problem occurs for FcX since the boundary condition for

(i) ∂xe pe S0 at z = 0 does not match the ones of {Hj }. For this reason we exclude the boundary {z = 0} from the definition of a control

volume V as in Section 6. References [1] R.A. Wooding, S.W. Tyler, I. White, Convection in groundwater below an evaporating salt lake: 1. onset of instability, Water Resour. Res. 33 (1997) 1199–1217. [2] D.A. Nield, A. Bejan, Convection in Porous Media, 3rd ed., Springer-Verlag, New York, 2006. [3] C.J. van Duijn, G.J.M. Pieters, R.A. Wooding, A. van der Ploeg, Stability criteria for the vertical boundary layer formed by throughflow near the surface of a porous medium, in: P.A.C. Raats, D. Smiles, A.W. Warrick (Eds.), Environmental Mechanics — Water, Mass and Energy Transfer in the Biosphere, in: Geophysical Monographs, vol. 129, American Geophysical Union, 2002, pp. 155–169.

[4] G.J.M. Pieters, C.J. van Duijn, Transient growth in linearly stable gravitydriven flow in porous media, Eur. J. Mech. B/Fluids 25 (2006) 83–94. doi:10.1016/j.euromechflu.2005.04.008. [5] R.A. Wooding, S.W. Tyler, I. White, P.A. Anderson, Convection in groundwater below an evaporating salt lake: 2. evolution of fingers or plumes, Water Resour. Res. 33 (6) (1997) 1219–1228. [6] R.N. Horne, M.J. O’Sullivan, Oscillatory convection in a porous medium heated from below, J. Fluid Mech. 66 (1974) 339–352. [7] R.N. Horne, M.J. O’Sullivan, Origin of oscillatory convection in a porous medium heated from below, Phys. Fluids 21 (1978) 1260–1264. [8] H.M. Schuttelaars, Nonlinear long term equilibrium profiles in a short tidal embayment, in: J. Donkers, M.B.A.M. Scheffers (Eds.), Physics of Estuaries and Coastal Seas, Balkema, Rotterdam, 1998, pp. 337–343. [9] P.C.F. van der Vaart, H.M. Schuttelaars, D. Calvete, H.A. Dijkstra, Instability of time-dependent wind-driven ocean gyres, Phys. Fluids 14 (2002) 3601–3615. doi:10.1063/1.1503804. [10] D. Calvete, H. de Swart, A nonlinear model study on the long-term behavior of shoreface-connected ridges, J. Geophys. Res. 108 (C5) (2003) 3169, doi:10.1029/2001JC001091. [11] G. de Josselin de Jong, The simultaneous flow of fresh water in aquifers of large horizontal extension determined by shear flow and vortex theory, in: A. Verruijt, F.B.J. Barends (Eds.), Proceedings Euromech, vol. 143, Balkema, Rotterdam, 1981, pp. 132–149. [12] Y.A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer-Verlag, New York, 2004. [13] A. Mielke, The Ginzburg–Landau equation in its role as a modulation equation, in: B. Fiedler (Ed.), Handbook for Dynamical Systems, Elsevier, Holland, 2002, pp. 759–834. [14] W. Eckhaus, Studies in Non-Linear Stability Theory, in: Tracts in Natural Philosophy, vol. 6, Springer-Verlag, Berlin, 1965. [15] G.J.M. Pieters, Stability and evolution of gravity-driven flow in porous media, Ph.D. Thesis, Eindhoven University of Technology, 2004. [16] A.C. Newell, J.A. Whitehead, Finite bandwidth, finite amplitude convection, J. Fluid Mech. 28 (1969) 279–303. [17] K. Fujimura, Centre manifold reduction and the Stuart–Landau equation for fluid motions, Proc. R. Soc. Lond. A 453 (1997) 181–203. [18] H.B. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems, in: P.H. Rabinowitz (Ed.), Applications of Bifurcation Theory, Academic Press, New York, 1977, pp. 359–384. [19] C.W. Gardiner, Handbook of Stochastic Methods: For Physics, Chemistry and the Natural Sciences, 2nd ed., Springer Verlag, Berlin, 2004.