Pattern formation in anticonvective systems

The famous convection experiments performed by Bénard (1900, 1901) in 1900 initialized a great deal of effort to understand and describe convective instabilities in vertically heated fluid systems over the last century. Even though this convective instability is induced by thermocapillarity at the free surface, buoyancy as a driving force may also be seen in the framework of Bénard's experiments. The occurrence of the Bénard roles was explained theoretically with thermocapillary forces (Bénard–Marangoni) as well as with buoyancy forces (Rayleigh–Bénard) in early works. However, these investigations considered the upper gas layer as passive, and analyses were performed mainly in one-layer models. It is natural to extend the investigation of these systems to configurations with more than one layer.

Here we focus on systems heated from above. Experiments for different fluid combinations were performed by Juel et al (2000) and showed the rich variety of pattern formation in two-layer systems. Note that experimentally both thermocapillary and buoyancy effects are present. Theoretical investigations of the pure Bénard–Marangoni instabilities in two-layer systems in the frame of amplitude equations were performed by Colinet et al (1996) and showed also a rich variety of bifurcations. Note that these calculations were done with the Marangoni effect only and were done close to a codimension two point.

Multi-layer systems heated from above with buoyancy effects alone were not investigated until 1964. It was argued that multi-layer systems driven by buoyancy effects become only unstable on applying heating from below. Potential energy release was seen as being impossible for heating from above and consequently the system rests in its conductive state. However, Welander (1964) showed in 1964 that a two-layer system heated from above can overcome this stabilizing mechanism. He used the expression anticonvection for that paradoxical kind of instability. Linear investigations of vertical bounded (finite) systems were done by Gershuni and Zhukhovitskii (1980) in 1980. They found the range of material parameters in which anticonvection can occur. First results of nonlinear investigations were also done at this time by Simanovskii (1980) and Gershuni et al (1981). Since anticonvection arises due to various interactions at the interface, only special fluids allow for an instability. However, Ingel (1997) revealed that an (infinitely extended) water–air system is anticonvectively unstable if phase transitions (evaporation) or the salinity of water is taken into account. Nepomnyashchy et al (2000) and Nepomnyashchy and Simanovskii (2001) showed that every (bounded) two-layer system can become anticonvectively unstable for heat sources or sinks at the interface (e.g. evaporation, chemical reactions and so on). Furthermore, they proposed an experimental realization of interfacial heating (infrared light sources for a silicon oil–ethylenglycol system; (Nepomnyashchy and Simanovskii 2001)). Investigations on combined actions of buoyancy and thermocapillarity with heat sinks and sources at the interface performed by Nepomnyashchy and Simanovskii (2002) in 2002 showed that a distinction between thermocapillary and anticonvective domains is possible. Also Simanovskii et al (2009) dealt with the joint action of both effects and revealed some new types of nonlinear traveling waves. The effects of an imposed horizontal temperature gradient in anticonvective systems were reported by Simanovskii et al (2002) and may be considered to give an insight into more realistic scenarios. Note that experiments for two-layer Rayleigh–Bénard systems heated from above in the context of anticonvection have not be done so far, since surface tension effects can cover the anticonvective instabilities (see, e.g., Degen et al (1998) and references therein). However, all the above-mentioned theoretical and numerical investigations indicate possible setups and parameter regimes for future experiments.

The increased computer power available nowadays allows direct numerical three-dimensional (3D) investigations for two-layer systems. For instance, Boeck et al (2002) showed pattern formation in two-layer systems driven by surface tension effects. Therefore the numerical investigation of pattern formation in anticonvective systems is in the offing. Full 3D simulations of anticonvective systems can be found in Merkt (2005) and Boeck et al (2008).

In this paper, we perform analytical as well as numerical investigations of a two-layer system with buoyancy effects heated from above. First, we derive the governing equations for the perturbated field quantities from the basic hydrodynamic equations. Then we perform numerically a linear stability analysis of a vertically bounded system. This allows us to reveal the effect of the boundaries on the system. Our linear results are compared with the analytical expression for the marginal control parameter of a vertically infinitely extended system. The derivation of this analytical expression can be found in the appendix. Further on, we show 3D results of direct numerical simulations of an anticonvective system bounded by rigid plates at the top and bottom. We focus on a mercury–water system for several reasons. On the one hand, the large number of free parameters complicates general predictions for two-layer systems, even numerically. Furthermore, mercury–water setups are the only known systems which show anticonvection without heat sources at the interface. On the other hand, experiments can be performed more or less easily for a mercury–water system.

The mechanism of anticonvection has to be considered as a subtle combination of effects in the bulk and at the interface. Consider a system with large thermal expansion coefficient (α₁ > α₂) and large thermal diffusivity (κ₁ > κ₂) in the lower layer (figure 1). Due to buoyancy a positive temperature fluctuation close to the interface in the lower layer forces the fluid to move to the interface (point 1). At the interface, viscous coupling induces a (passive) convective motion in the upper layer (point 2). Since buoyancy and thermal diffusivity in the upper layer are small, hot fluid from above is transported to the interface (i.e. the upper layer serves as a heat supplier and is called a type I system; (Nepomnyashchy et al 2000)). Thermal coupling allows then for heat transfer to the lower layer (point 3) and, finally, thermal diffusion (point 4) can amplify the primary positive temperature fluctuation in the lower layer. Note that this is an essentially different instability mechanism than the usual Rayleigh–Bénard instability driven by buoyancy alone.

Zoom In Zoom Out Reset image size

To summarize, anticonvection can arise if at least four timescales match: buoyant timescale in the bulk, viscous timescale, thermal timescale and heat flux timescale at the interface.

Two-layer systems with layer thickness d₁ for the lower fluid and d₂ for the upper fluid (ratio d = d₂/d₁) are described by the Navier–Stokes equations, energy equation and incompressibility condition for each fluid layer:

$$\begin{equation} \rho_i\left(\partial_t{\bf v}_i + ({\bf v}_i\cdot\nabla)\cdot {\bf v}_i\right) = \rho_i\nu_i\triangle{\bf v}_i - \nabla P_i + g\rho_i{\bf e}_z, \end{equation} \tag{ 1a }$$

$$\begin{equation} \partial_t T_i + {\bf v}_i\cdot(\nabla T_i) = \kappa_i\triangle T_i, \end{equation} \tag{ 1b }$$

$$\begin{equation} \nabla\cdot{\bf v}_i = 0. \end{equation} \tag{ 1c }$$

The subscripts i = 1,2 denote the lower and the upper fluid, respectively. ρ_i are the densities, ν_i the kinematic viscosities, g the acceleration due to gravity, κ_i the thermal diffusivities and e_z the unit vector in the z-direction. v_i are the velocity vectors, P_i the pressures and T_i the temperatures.

Rigid boundaries at the bottom and top of the system provide (z = 0 and z = d₁ + d₂)

$$\begin{equation} \mbox{at } z=0: \quad {\bf v}_1 = 0,\quad \frac{\partial}{\partial z}v_{1z} = 0,\quad T_1 = T_{\mathrm{b}}, \end{equation} \tag{ 2a }$$

$$\begin{equation} \mbox{at } z=d_1+d_2: \quad {\bf v}_2 = 0,\quad \frac{\partial}{\partial z}v_{2z} = 0,\quad T_2 = T_{\mathrm{t}}, \end{equation} \tag{ 2b }$$

with T_b and T_t being the temperature of the lower and upper plates, respectively. The conditions ∂_zv_iz = 0 result from both v_i = 0 and the continuity equations (1c). We consider an undeformable interface. Therefore the interface conditions can be written as (z = d₁)

$$\begin{equation} v_{1x} = v_{2x}, \quad v_{1y}=v_{2y},\quad v_{1z} = v_{2z} =0, \end{equation} \tag{ 3a }$$

$$\begin{equation} \rho_1\nu_1\frac{\partial}{\partial z} v_{1x} - \rho_2\nu_2\frac{\partial}{\partial z}v_{2x} = 0, \quad \rho_1\nu_1\frac{\partial}{\partial z} v_{1y} - \rho_2\nu_2\frac{\partial}{\partial z}v_{2y} = 0, \end{equation} \tag{ 3b }$$

$$\begin{equation} T_1 = T_2, \quad \lambda_1\frac{\partial}{\partial z} T_1= \lambda_2\frac{\partial}{\partial z}T_2 \end{equation} \tag{ 3c }$$

with λ_i for thermal conductivities.

In the following, we use the Boussinesq approximation, i.e. all material parameters are independent of temperature except for densities in the gravity term (equation (1a)). A linear expansion provides

$$\begin{equation} \rho_1(T) = \rho_{1_0}\left(1- \alpha_1(T-T_{\mathrm{b}})\right), \end{equation} \tag{ 4a }$$

$$\begin{equation} \rho_2(T) = \rho_{2_0}\left(1- \alpha_1(T-T_i)\right), \end{equation} \tag{ 4b }$$

where T_i is a (stationary) reference temperature of the interface. ρ₁₀ and ρ₂₀ are reference densities in fluid 1 (at T_b) and fluid 2 (at T_i), respectively. α₁ = −1/ρ₁₀·∂ρ₁/∂T|_Tb and α₂ = −1/ρ₂₀·∂ρ₂/∂T|_Ti denote the thermal expansion coefficients in fluids 1 and 2, respectively.

We use the nondimensional variables

$$\begin{equation} x_i = d_1 {\tilde x_i}, \quad v_i = \frac{\kappa_1}{d_1} {\tilde v_i}, \quad t = \frac{d_1^2}{\kappa_1} {\tilde t},\quad T = \Delta T\; {\tilde T},\quad P = \frac{\nu_1\kappa_1\rho_{1_0}}{d_1^2}{\tilde P}, \end{equation} \tag{ 5 }$$

where ΔT = T_b − T_t denotes the temperature difference of the bottom and top plates. We find that, for the stationary nondimensional temperature profile (we dropped the tildes),

$$\begin{equation} T_1(z) = T_{\mathrm{b}} + \beta_1 z, \end{equation} \tag{ 6a }$$

$$\begin{equation} T_2(z) = T_{\mathrm{b}} + (\beta_1 - \beta_2) - \beta_2 z, \end{equation} \tag{ 6b }$$

where

$$\begin{equation} \beta_1 = \frac{1}{1+\frac{d}{\lambda}}, \end{equation} \tag{ 7a }$$

$$\begin{equation} \beta_2 = \frac{1}{\lambda + d} \end{equation} \tag{ 7b }$$

are nondimensional temperature gradients in fluids 1 and 2, respectively. The conductive state is described by

$$\begin{equation} {\mathbf v}_{1,2} = 0. \end{equation} \tag{ 8 }$$

Finally, we obtain from equations (1a)–(1c) for the perturbated quantities (u for the velocity, θ for the temperature and p for the pressure) a set of equations

$$\begin{equation} \frac{1}{Pr}\left(\partial_t {\bf u}_1 + ({\bf u}_1\cdot\nabla)\cdot{\bf u}_1\right) + \nabla p_1 = R\;{\bf e}_z\; \theta_1 + \triangle {\bf u}_1, \end{equation} \tag{ 9a }$$

$$\begin{equation} \partial_t \theta_1 + {\bf u}_1\cdot(\nabla\theta_1) = \beta_1 u_{1z} + \triangle\theta_1, \end{equation} \tag{ 9b }$$

$$\begin{equation} \nabla\cdot {\bf u}_1 = 0, \end{equation} \tag{ 9c }$$

$$\begin{equation} \frac{1}{Pr}\left(\partial_t {\bf u}_2 + ({\bf u}_2\cdot\nabla)\cdot{\bf u}_2\right) + \nabla p_2 = \alpha\;R\;{\bf e}_z\; \theta_1 + \nu\triangle {\bf u}_1, \end{equation} \tag{ 9d }$$

$$\begin{equation} \partial_t \theta_2 + {\bf u}_2\cdot(\nabla\theta_2) = \beta_2 u_{2z} + \kappa\triangle\theta_2, \end{equation} \tag{ 9e }$$

$$\begin{equation} \nabla\cdot {\bf u}_2 = 0. \end{equation} \tag{ 9f }$$

Due to the Boussinesq approximation in each layer, the thermal expansion coefficients α_i appear in the bulk equations (hidden in R).

The boundary conditions (2a) and (2b) transform to

$$\begin{equation} \mbox{bottom}\quad (z=0):~ {\bf u}_1 = 0, ~ \frac{\partial}{\partial z}u_{1z} = 0, \quad \theta_1 = 0, \end{equation} \tag{ 10a }$$

$$\begin{equation} \mbox{top}\quad (z=1+d): {\bf u}_2 = 0, \quad \frac{\partial}{\partial z}u_{2z} = 0, \quad \theta_2 = 0, \end{equation} \tag{ 10b }$$

and from equations (3a)–(3c) the interface conditions at z = 1 read

$$\begin{equation} u_{1x} = u_{2x}, \quad u_{1y} = u_{2y}, \quad u_{1z} = u_{2z} = 0, \end{equation} \tag{ 11a }$$

$$\begin{equation} \frac{\partial}{\partial z} u_{1x} - \rho\nu \frac{\partial}{\partial z} u_{2x} = 0, \quad \frac{\partial}{\partial z} u_{1y} - \rho\nu \frac{\partial}{\partial z} u_{2y} = 0, \end{equation} \tag{ 11b }$$

$$\begin{equation} \theta_1 = \theta_2, \quad \frac{\partial}{\partial z}\theta_1 = \lambda\frac{\partial}{\partial z}\theta_2. \end{equation} \tag{ 11c }$$

The necessary material parameters of the system are the ratios of the layer thicknesses d = d₂/d₁, the densities ρ = ρ₂/ρ₁, thermal conductivities λ = λ₂/λ₁, thermal diffusivities κ = κ₂/κ₁, kinematic viscosities ν = ν₂/ν₁ and the thermal expansion coefficients α = α₂/α₁.

Further on, two dimensionless numbers result:

$$\begin{equation} Pr = \frac{\nu_1}{\kappa_1}\quad \mbox{Prandtl number,} \end{equation} \tag{ 12a }$$

$$\begin{equation} R = \frac{\alpha_1\,g \Delta T}{\nu_1\kappa_1}d_1^3\quad \mbox{Rayleigh number}, \end{equation} \tag{ 12b }$$

where the Rayleigh number R is the control parameter of the system.

2.1. Linearized equations

Equations (9a)–(9f) and the corresponding boundary and interface conditions (equations (10a) and (10b) and (11a)–(11c)) allow formulation as a general linearized eigenvalue problem. Therefore a normal mode ansatz for both fluids is used in the usual manner

$$\begin{equation} u_{z_i}(x,y,z,t) = \varphi_i(z) \exp{(\mathrm{i}\,\vec k\cdot\vec x+\chi t),} \end{equation} \tag{ 13a }$$

$$\begin{equation} {\theta}_i(x,y,z,t) = \vartheta_i(z) \exp{(\mathrm{i}\,\vec k\cdot\vec x - \mathrm{i}\pi/2 +\chi t)}, \end{equation} \tag{ 13b }$$

where $\vec x= (x,y)$, $\vec k = (k_x,k_y)$, $k=|\vec k|$ and χ is the growth rate. The velocities u_xi and u_yi are eliminated using the incompressibility conditions equations (9c) and (9f)). The phase shift −iπ/2 in equation (13b) ensures real equations in the following. Subsequently applying the curl and neglecting nonlinearities provides linear equations for fluid 1 (see also Chandrasekhar (1961) and Merkt and Bestehorn (2003))

$$\begin{equation} \chi\frac{1}{Pr}\left(\frac{\partial^2}{\partial z^2} -k^2 \right)\varphi_1 = \left(\frac{\partial^2}{\partial z^2} -k^2 \right)^2\varphi_1 - R k^2\vartheta_1, \end{equation} \noindent \tag{ 14a }$$

$$\begin{equation} \chi\vartheta_1 = \left(\frac{\partial^2}{\partial z^2} -k^2 \right)\vartheta_1 + \beta_1\varphi_1, \end{equation} \tag{ 14b }$$

and for fluid 2

$$\begin{equation} \chi\frac{1}{Pr}\left(\frac{\partial^2}{\partial z^2} -k^2 \right)\varphi_2 = \nu\left(\frac{\partial^2}{\partial z^2} -k^2 \right)^2\varphi_2 - \alpha R k^2\vartheta_2, \end{equation} \tag{ 14c }$$

$$\begin{equation} \chi\vartheta_2 = \kappa \left(\frac{\partial^2}{\partial z^2} -k^2 \right)\vartheta_2 + \beta_2\varphi_2. \end{equation} \tag{ 14d }$$

The boundary conditions (10a) and (10b) turn out to be

$$\begin{equation} \varphi_i = 0,\quad \frac{\partial}{\partial z}\varphi_i =0,\quad \vartheta_i =0, \end{equation} \noindent \tag{ 15 }$$

and the interface conditions (11a)–(11c) read

$$\begin{equation} \varphi_1 = \varphi_2 = 0,\quad \frac{\partial}{\partial z}\varphi_1 = \frac{\partial}{\partial z}\varphi_2,\quad \frac{\partial^2}{\partial z^2}\varphi_1 = \rho\nu\frac{\partial^2}{\partial z^2}\varphi_2, \end{equation} \tag{ 16a }$$

$$\begin{equation} \vartheta_1 = \vartheta_2,\quad \frac{\partial}{\partial z}\vartheta_1 = \lambda\frac{\partial}{\partial z}\vartheta_2. \end{equation} \tag{ 16b }$$

3.1. Linear results

The general eigenvalue problem of the bounded system (equations (14)–(16)) for the growth rate χ is solved numerically.

First, we investigate the influence of the ratio of thermal expansion coefficients on the instability. Figure 2 shows the marginal stability line in the R–α-plane for d = 1 with parameters from table 1. Decreasing α leads to a smaller absolute value of the critical Rayleigh number R_c. Due to the anomaly of water the thermal expansion coefficient of water α₂ will be negative for temperatures below 4 °C. Therefore we suppose a monotonically decreasing function for α₂ and also for the ratio α = α₂/α₁ with decreasing temperature. In the following, we use α = 0.1144. This is reasonable since experiments can be performed with an upper plate temperature less than 20 °C.

Zoom In Zoom Out Reset image size

Table 1. Material parameters of a mercury–water system at 20 °C. Due to the anomaly of water, the ratio α decreases for decreasing temperatures. We take α = 0.1144 for our investigations and assume the other parameters to be constants.

	Mercury (lower layer)	Water (upper layer)	Ratio
Density ρ (gm−3)	13.546×106	0.9982×106	0.0737
Viscosity ν (m2 s−1)	1.15×10−7	1.0×10−6	8.696
Thermal diffusivity κ (m2 s−1)	4.35×10−6	1.43×10−7	0.0328
Thermal expansion α (K−1)	1.81×10−4	2.07×10−4	1.144
Thermal conductivity λ (W (mK)−1)	8.2	0.598	0.073
Prandtl number Pr	0.026	7.0

Figure 3 (solid thick line) shows the marginal Rayleigh number R versus k for different layer ratios d. We find that the critical Rayleigh number R_c ≈ −2950 and the critical wave number k_c ≈ 1.33 for d = 1 (thick black line). Surprisingly, this critical wave number k_c can be stabilized by stronger heating (increasing the control parameter −R). This is in contrast to the normal Rayleigh–Bénard instability (heating from below) where the marginal R–k diagram shows a parabolic shape (i.e. an unstable wave number stays unstable by increasing the control parameter). Further on, we show the analytical expression

$$\begin{equation} R_{\mathrm{c}} = R_{\mathrm{c}}(k,\alpha,\nu,\kappa,\lambda,\beta_1)\propto - k^4 \end{equation} \tag{ 17 }$$

(see equation (A.11) in the appendix) for infinitely extended systems in the z-direction (thin lines). Equation (17) is valid for small k values and large Rayleigh numbers. Note that β₁ is the dimensionless temperature gradient in the lower fluid which is a parameter for the infinitely extended system. To compare vertically bounded and infinitely extended systems we use the numerical value of β₁ from equation (7a) (with d from the vertically bounded system). A detailed derivation and discussion of expression (17) can be found in the appendix. The effect of the rigid boundaries for small and large wavelengths is clearly visible. Both small and large wavelengths force a bending of the marginal stability line ((a) and (b) branches).

Zoom In Zoom Out Reset image size

Small wavelengths (a) cause two convection roles on top of each other in the upper layer. Therefore the heat transport from the upper boundary to the interface is decreased. Large wavelengths (b) inhibit a passive convective motion in the upper layer since the horizontal motion along the interface is slowed down and finally stopped by viscosity. Note that the latter mechanism already occurs in the 'normal' Rayleigh–Bénard convection.

Figure 4 shows the eigenfunctions in the unstable domain for different wave numbers k using the same Rayleigh number (R ≈ −2980). For the critical wave number (solid line) the eigenfunctions of the velocity and the temperature perturbation are positively correlated in the upper layer (sign(u_z) = sign(θ)), i.e. potential energy is released in the upper layer. The eigenfunctions are changed essentially with increasing or decreasing the wave number, respectively. Firstly, the stabilization for small wavelengths (k > k_c, dashed line) for finite extensions in the z-direction is caused by the appearance of an anticorrelated region of u_z and θ at the boundary of the upper fluid, namely sign(u_z) = −sign(θ). This reduces the release of potential energy in the upper layer. A further increase of k spreads this anticorrelated region and destabilization is not possible for a large enough k. Note that an increase of the absolute value of the control parameter has the same effect for any fixed wave number k. Consequently, a large enough −R stabilizes every wave number k and explains branch (a) in figure 3.

Zoom In Zoom Out Reset image size

Secondly, for large wavelengths (k < k_c, dotted line, see also branch (b) in figure 3) the eigenfunction of the velocity u_z tends to zero. Therefore the control parameter of the system has to be increased for large wavelengths.

3.2. Nonlinear results

In this subsection, we solve the fully nonlinear equations (9a)–(9f) numerically. We use periodic boundary conditions laterally and a semi-implicit pseudospectral method. The discretization in the z-direction is done with finite differences. The stationary state is perturbed by small randomly distributed fluctuations.

We redefine the control parameter as

$$\begin{equation} \epsilon = \frac{R - R_{\mathrm{c}}}{R_{\mathrm{c}}} \end{equation} \tag{ 18 }$$

and we use a 128 × 128 × 20 grid for each layer.

Figure 5 shows patterns of stationary states at t = 3000 and t = 1000 for  = 0.15 and  = 0.4, respectively. We use parameters of a mercury–water system (see table 1) with d = 1, α = 0.1144 and an aspect ratio of Γ = 18.8 (i.e. four times the critical wavelength). A transition from hexagons to squares takes place by increasing the control parameter . We have to note that squares are not found in the simulations of Boeck et al (2008). However, these transitions from hexagons to squares (and even weak turbulent regimes) are well known for single-layer Bénard–Marangoni convection via increasing the control parameter (Merkt and Bestehorn 2003).

Zoom In Zoom Out Reset image size

As was already stated in the introduction, every two-layer system may become anticonvectively unstable incorporating contaminations or heat sources (sinks) at the interface (Nepomnyashchy et al 2000). Therefore we investigate now an artificial situation via changing the ratio of thermal conductivities λ = λ₂/λ₁ for the former simulation. This implies that a change of the temperature profiles equations (7a) and (7b) may reflect in some way the influence of contaminations or heat sources at the interface. Figure 6 shows the stationary linear temperature profiles for different values of the ratio of thermal conductivities. For λ = 0.0365 the gradient is mostly in the upper layer and for λ = 20 mostly in the lower layer.

Zoom In Zoom Out Reset image size

Figure 7 shows the marginal Rayleigh number R versus k (the same representation as in figure 3) for different ratios of conductivities λ. The unstable states vanish for λ → 0 or λ → ∞. Further on, the critical wave number k_c (and the critical Rayleigh number R_c) depends on λ. Accumulation of these critical points provides the critical line. The minimal R_c is found for λ = 1 with this parameter set. The shapes for λ = 0.0365 and λ = 20 are more or less similar. To show the different nonlinear behavior, we solved the full problem numerically for these two ratios of conductivities with  = 0.3.

Zoom In Zoom Out Reset image size

In figure 8, 3D plots of the stationary temperature perturbation fields are shown (with t = 7000). The resolution is 128 × 128 × 20 for each layer and we take an aspect ratio Γ ≈ 42 that corresponds to eight critical modes. For λ = 0.0365 (plot (a)) the stationary patterns are squares and the temperature perturbations are extended into the bulk of the upper fluid. Further on, the temperature perturbations in the upper fluid are 10 times larger than those in the lower fluid, whereas the velocity perturbations are the other way round (not shown). Therefore a strong convective motion occurs basically in the lower layer.

Zoom In Zoom Out Reset image size

However, for λ = 20 (figure 8(b)) we find hexagonal patterns. The temperature and vertical velocity perturbations are of the same order in both layers. The temperature perturbations are limited to the interface. Furthermore, the absolute value of the velocity perturbations correspond to the velocity perturbations of the upper layer in the former case, i.e. the convective motions in both layers are of the same order as the weak convective motion in the upper layer for λ = 0.0365.

In both situations, the upper layer serves as a heat supplier (α < 1, κ < 1, see the introduction), whereas the underlying convective energy transport fed in the system is fundamentally different and results in different patterns and realized flows.

Anticonvection may appear in two-layer systems via heating from above with buoyancy forces only. This instability mechanism needs a sensitive interaction of bulk and interface forces. At least four timescales must match to yield an anticonvective instability. We have derived the governing equations for the perturbed field quantities. First, we performed a linear stability analysis of the bounded mercury–water system. The linear stability analysis has shown that one unstable k mode is stabilized by further increasing the control parameter −R. Subsequently, we compared vertically bounded with vertically infinitely extended systems. Therefore we derived an analytical expression for the critical Rayleigh number R_c = R_c(k,α,ρ,ν,λ,κ,β₁)∝ − k⁴ in vertically infinitely extended systems. We showed that rigid boundaries lead to an increase of the conductive regime. Subsequently, we presented the fully nonlinear dynamics of an anticonvective system. Hexagons arise for a mercury–water system close to threshold. Increasing the applied temperature gradient leads to a transition to squares. Moreover, two different stationary temperature profiles have revealed how convective heat and mass transport is realized and consequently may influence the resulting stationary patterns. These latter investigations may represent in some way the influence of contaminations or heat sources (sinks) at the interface and provide expectations for experimental setups. Finally, we note that we have not performed an exhaustive investigation of the parameter range for anticonvection due to the amount of parameters (namely at least four) that exceeds our computational facilities. However, experimental results and further theoretical insights may allow for a more restrictive parameter space.

Here we derive an analytical expression for the marginal Rayleigh number R_c of a vertically infinitely extended system. For the sake of simplicity we shift the interface to z = 0 in this section.

The linear equations for the marginal state (equations (14a)–(14d) with growth rate χ = 0) can be written as (e.g. (∂²_z − k²) ×  equation (14a) + R_ck² ×  equation (14b), see Welander (1964) and Chandrasekhar (1961)):

$$\begin{equation} (\partial_z^2-k^2)^3 \varphi_1 = -R_{\mathrm{c}}\beta_1 k^2 \varphi_1, \end{equation} \tag{ A.1a }$$

$$\begin{equation} (\partial_z^2-k^2)^3 \varphi_2 = -R_{\mathrm{c}}\beta_1 k^2 \varepsilon_4^3 \varphi_2. \end{equation} \tag{ A.1b }$$

The six conditions at the interface (see equations (16a) and (16b)) turn out to be

$$\begin{equation} \varphi_1 = 0, \end{equation} \tag{ A.2a }$$

$$\begin{equation} \varphi_2 = 0, \end{equation} \tag{ A.2b }$$

$$\begin{equation} \partial_z \varphi_1 = \partial_z \varphi_2, \end{equation} \tag{ A.2c }$$

$$\begin{equation} \partial_z^2 \varphi_1 = \varepsilon_1 \partial_z^2 \varphi_2, \end{equation} \tag{ A.2d }$$

$$\begin{equation} (\partial_z^2 - k^2)^2 \varphi_1 = \varepsilon_2 (\partial_z^2 -k^2)^2 \varphi_2, \end{equation} \tag{ A.2e }$$

$$\begin{equation} \partial_z(\partial_z^2 - k^2)^2 \varphi_1 = \varepsilon_3 \partial_z(\partial_z^2 -k^2)^2 \varphi_2, \end{equation} \noindent \tag{ A.2f }$$

with ε₁ = ρν, ε₂ = ν/α, ε₃ = νλ/α and ε₄ = (α/(νκλ))^(1/3). Note that β₁ is a parameter of the infinitely extended systems and cannot be expressed with equation (7a).

The six roots for fluid 1 (substitute $\varphi _1\propto \exp {(rz)}$ in equation (A.1a)) are

$$\begin{equation} r_{1,4} = \pm\underbrace{\sqrt{k^2 + G}}_{t_1}, \end{equation} \tag{ A.3a }$$

$$\begin{equation} r_{2,3} = \underbrace{\frac{1}{2}\sqrt{\sqrt{z_1^2 + 3G^2} + z_1}}_{t_2}\; \mp\; \mathrm{i}\underbrace{\frac{1}{2}\sqrt{\sqrt{z_1^2 + 3 G^2} - z_1}}_{t_3}, \end{equation} \tag{ A.3b }$$

$$\begin{equation} r_{5,6} = - r_{2,3}, \end{equation} \tag{ A.3c }$$

and for fluid 2 (substitute $\varphi _2\propto \exp {(qz)}$ in equation (A.1b))

$$\begin{equation} q_{1,4} = \pm\underbrace{\sqrt{k^2 + \varepsilon_4G}}_{s_1}, \end{equation} \tag{ A.4a }$$

$$\begin{equation} q_{2,3} = \underbrace{\frac{1}{2}\sqrt{\sqrt{z_2^2 + 3\varepsilon_4 G^2} + z_2}}_{s_2}\; \mp\; \mathrm{i} \underbrace{\frac{1}{2}\sqrt{\sqrt{z_2^2 + 3\varepsilon_4 G^2} - z_2}}_{s_3}, \end{equation} \tag{ A.4b }$$

$$\begin{equation} q_{5,6} = - q_{2,3}, \end{equation} \tag{ A.4c }$$

with

$$\begin{equation} G = \sqrt[3]{-R_{\mathrm{c}}\,k^2\beta_1},\quad z_1 = (2 k^2 - G), \quad z_2 = (2 k^2 - \varepsilon_4G) \end{equation} \tag{ A.5 }$$

and the assumption

$$\begin{eqnarray} &&G > 0 \quad \Rightarrow\quad R_{\mathrm{c}}<0 \end{eqnarray} \tag{ A.6a }$$

that ensures positive arguments of the square roots in equations (A.3) and (A.4).

Since perturbations decay at large distances from the interface, we take into account only the 'active' modes at the interface, i.e. three modes with the positive real part of the roots for fluid 1 and 3 modes with the negative real part of the roots for fluid 2:

$$\begin{equation} \varphi_1 = \sum_{n=1}^{3}{C_{n}\mathrm{e}^{r_{n} z}},\quad\quad \varphi_2 = \sum_{n=1}^{3}{D_{n}\mathrm{e}^{-q_{ n} z}}. \end{equation} \tag{ A.7 }$$

This ansatz provides an algebraic system of equations for the amplitudes C_n and D_n with nontrivial solutions det(M) = 0:

$$\begin{eqnarray*} \fl \left|\underbrace{ \begin{array}{@{}c@{\hskip6pt}c@{\hskip7pt}c@{\hskip8pt}c@{\hskip8pt}c@{\hskip7pt}c@{\hskip7pt}} 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 \\ r_1 & r_2 & r_3 & q_1 & q_2 & q_3 \\ r^2_1 & r^2_2 & r^2_3 & \!-\varepsilon_1 q^2_1 & \!-\varepsilon_1 q^2_2 & \!-\varepsilon_1 q^2_3\\ (r^2_1 \!-\! k^2)^2 & (r^2_2 \!-\! k^2)^2 & (r^2_3 \!-\! k^2)^2 & \!-\varepsilon_2 (q^2_1 \!-\! k^2)^2 & \!-\varepsilon_2 (q^2_2 \!-\! k^2)^2 & -\varepsilon_2 (q^2_3 \!-\! k^2)^2 \\ r_1(r^2_1 \!-\! k^2)^2 & r_2(r^2_2\!-\!k^2)^2 & r_3(r^2_3 \!-\! k^2)^2 & \varepsilon_3 q_1(q^2_1 \!-\! k^2)^2 & \varepsilon_3 q_2(q^2_2\!-\!k^2)^2 & \varepsilon_3 q_3(q^2_3 \!-\! k^2)^2 \end{array}}_{{\bf M}} \right| &\! = & 0. \end{eqnarray*}$$

After some manipulations we find the second-order surface equation of Welander (1964) in the ε₁, ε₂ and ε₃ space. It can also be written as a quadratic form in the real and imaginary parts of the roots in equations (A.3) and (A.4), namely $\vec V = (\vec t^{\mathrm {T}},\vec s^{\mathrm {T}})=(t_1,t_2,t_3,s_1,s_2,s_3)^{\mathrm {T}}$.

$$\begin{equation} \vec V\mathrm{^T} {\bf H} \vec V = 0,\quad\quad {\bf H}= \left( \begin{array}{@{}cc@{}} \bf H_{11} & \bf H_{12} \\ \bf H_{12}^T & \bf H_{22} \end{array} \right). \end{equation} \tag{ A.8 }$$

Note that H is symmetric and depends on material parameters alone. The 3 × 3 submatrices can be written as:

$$\begin{eqnarray*} {\bf H}_{11} = \left( \begin{array}{@{}ccc@{}} 0 & -h_1 & -\sqrt{3}h_1 \\ -h_1 & \quad 2 h_1\quad & 0 \\ -\sqrt{3}h_1 & 0 & 2 h_1 \end{array} \right), &\quad & {\bf H}_{22} = \frac{h_4}{h_1}\cdot {\bf H}_{11}, \end{eqnarray*}$$

$$\begin{eqnarray*} \fl {\bf H}_{12} & = & \left(\!\begin{array}{@{}ccc@{}} \frac{3}{2}\left(h_2-h_3-\frac{4}{81}\frac{h_1 h_4}{h_3}\right) & \frac{3}{2}\left(h_3-h_2-\frac{8}{81}\frac{h_1 h_4}{h_3}\right) & -\frac{3}{2}\sqrt{3}(h_2+h_3) \\ -\frac{3}{2}\left(h_2+2h_3-\frac{4}{81}\frac{h_1 h_4}{h_3}\right) & \quad\frac{3}{2}\left(h_2+2h_3-\frac{8}{81}\frac{h_1 h_4}{h_3}\right)\quad & \frac{3}{2}\sqrt{3}(h_2-2h_3) \\[8 pt] \frac{3}{2}\sqrt{3}\left(-h_2 - \frac{4}{81}\frac{h_1 h_4}{h_3}\right) & \frac{3}{2}\sqrt{3}\left(h_2 - \frac{8}{81}\frac{h_1 h_4}{h_3}\right) & \frac{9}{2}h_2 \end{array}\!\right) \end{eqnarray*}$$

with $h_1 = \frac {9}{2}\varepsilon _4^2\varepsilon _1\varepsilon _2$, h₂ = ε³₄ε₂ε₃ + ε₁, h₃ = ε₄ε₂ and $h_4 = \frac {9}{2}\varepsilon _4\varepsilon _3$.

A.1. Conclusions for infinitely extended systems

Firstly, a principal axis transformation of equation (A.8) provides (after back substitution of material parameters)

$$\begin{eqnarray} \fl \frac{12}{5}\nu^2\rho^2\kappa^2\left(\frac{\alpha}{\nu\lambda\kappa}\right)^{4/3} \left(\frac{1}{2}(1+\sqrt{5})\cdot\tilde t_1^2 + \frac{1}{2}(1-\sqrt{5})\cdot\tilde t_2^2 + \tilde t_3^2\right) & = & (\rho\kappa\alpha -1)^2\cdot\tilde s_1^2. \end{eqnarray} \tag{ A.9 }$$

In general one finds two zero, two negative and two positive eigenvalues. This provides a single sheet hyperboloid in $(\tilde t_1,\tilde t_2,\tilde t_3)$-space ( $\tilde t_i$ depend on material parameters, whereas $\tilde s_1$ is independent). For ρ = 1/(κα) one negative eigenvalue becomes zero (rhs of equation (A.9)) and therefore the surface reduces to a cone. We note that a hyperboloid also results directly from equation (A.8) calculating the signature of H (using Sylvester's theorem (Golub and Van Loan 1996)).

Secondly, Welander (1964) has shown that R_c is a monotonic function of k. Therefore we try to extract a simple analytical expression for the marginal R_c from equation (A.8) for large-temperature gradients (i.e. large |R_c|). Equation (A.8) has 17 terms with products of the roots (equations (A.3) and (A.4)). We expand these 17 products of the roots at k = 0 to the order O(k^14/3) (i.e. three terms). For the sake of clarity, we present examples of only two of the 17 expanded products:

$$\begin{equation} t_1\cdot s_1 = \frac{1}{2}{\cal G}k^{2/3} + \frac{1}{2} k^{2} + \frac{3}{16{\cal G}}k^{10/3} +O(k^{14/3}), \end{equation} \tag{ A.10a }$$

$$\begin{equation} t_1\cdot s_3 = \frac{\sqrt{3\varepsilon_4}}{2}{\cal G}k^{2/3} + \frac{\sqrt{3}(\varepsilon_4-1)}{4\sqrt{\varepsilon_4}}k^{2} - \frac{\sqrt{3}(\varepsilon_4+2)}{16\sqrt{\varepsilon_4}{\cal G}}k^{10/3} + O(k^{14/3}), \end{equation} \tag{ A.10b }$$

with ${\cal G} = G\cdot k^{3/2}=\sqrt {-R_{\mathrm {c}}\beta _1}$ (see equation (A.5)). All expansions underlie the assumptions that, firstly, |R_c| is large (i.e. high orders in k vanish) and, secondly, moderate k values (k < 10) and ε₄ > 1 (i.e. expansion at k = 0 is valid). We have checked the validity of these assumptions numerically. The expansion for moderate k values and ε₄ > 1 provides a good approximation for the products of the roots (e.g. t₁·s₁,t₁·s₂, ...). To exemplify this, we plotted the products and the expansions of equations (A.10a) and (b)) in figure 9 for R_c = −1000. Even this small R_c shows good accordance for k < 4.

Zoom In Zoom Out Reset image size

Substituting the 17 expanded products into equation (A.8) and back substitution of ${\cal G}$ and h_i provides, finally, an analytical expression for the marginal Rayleigh number

$$\begin{equation} R_{\mathrm{c}} = \underbrace{-\left(\frac{A_2\pm \sqrt{A_2^2-4A_1A_3}}{2 A_1}\right)^3\frac{1}{\beta_1}}_{C}\,\cdot\, k^4 \end{equation} \tag{ A.11 }$$

with

$$\begin{eqnarray*} &&\fl A_1 = 4\varepsilon_4^2 \left[ \varepsilon_2\varepsilon_3\varepsilon_4^3 - 4\varepsilon_1\varepsilon_3\varepsilon_4^2 - 3(\varepsilon_1\varepsilon_2 + \varepsilon_3)\varepsilon_4^{3/2} - 4\varepsilon_2\varepsilon_4 + \varepsilon_1 \right], \\ &&\fl A_2 = 4\varepsilon_4\left[\varepsilon_2\varepsilon_3\varepsilon_4^4 +\varepsilon_3(\varepsilon_2-\varepsilon_1)\varepsilon_4^3 +3\varepsilon_1\varepsilon_2\varepsilon_4^{5/2} +2(\varepsilon_2 + \varepsilon_1\varepsilon_3)\varepsilon_4^2 +3\varepsilon_3\varepsilon_4^{3/2} +(\varepsilon_1-\varepsilon_2)\varepsilon_4 +\varepsilon_1 \right], \\ &&\fl A_3 = \varepsilon_2\varepsilon_3\varepsilon_4^5 + 2\varepsilon_3(\varepsilon_1 + 2\varepsilon_2)\varepsilon_4^4 +9 \varepsilon_1\varepsilon_2\varepsilon_4^{7/2} + (\varepsilon_2\varepsilon_3-\varepsilon_2 + 2\varepsilon_1\varepsilon_3)\varepsilon_4^3 + (\varepsilon_1-\varepsilon_1\varepsilon_3+2\varepsilon_2)\varepsilon_4^2 \\ &&+ 9\varepsilon_3\varepsilon_4^{3/2} +2(\varepsilon_2+2\varepsilon_1)\varepsilon_4 + \varepsilon_1, \end{eqnarray*}$$

where the proportionality factor C depends on material parameters alone. Note that two solutions exist for R_c in equation (A.11). However, only negative C has a physical meaning, since we restrict ourselves to negative Rayleigh numbers R. Therefore we take the first solution (i.e. ' + ' sign).

Finally, we have compared equation (A.11) with the results of the full 'Welander' equation (A.8). Surprisingly, the marginal stability line is perfectly described by a k⁴ dependence (even for violating the above-noted assumptions). We mention that the k⁴ proportionality results already from the expansion to order O(k^10/3). However, then the approximation of the coefficient of proportionality C' is no longer valid.