1 Introduction

The instability for a stationary solution of the local balance equations which govern the fluid flow is a cornerstone topic of the research in fluid mechanics over the last century. Its intrinsic importance is due to the close connection with the analysis of transitional flow and, ultimately, with turbulence. There is also a fundamental area of research where fluid mechanics is interrelated to heat transfer and convection. Within this area, the instability of a flow system pinpoints the conditions for the onset of convective cellular patterns as those identified in Bénard experiments [1] and modelled by Lord Rayleigh [2]. The typical setup, well-known as the Rayleigh-Bénard system, is a horizontal layer of fluid or fluid-saturated porous medium where the horizontal boundaries are kept isothermal with heating from below [3,4,5,6].

The classical strategy for the flow stability analysis is an implementation of Lyapunov’s idea of instability for a mechanical system. In fact, as is well-known, Lyapunov’s concept of stability or instability for an equilibrium state of a mechanical system is based on the analysis of the time evolution undergone by the system as a consequence of a slight initial perturbation. The instability occurs when the perturbation of the equilibrium state is amplified in the time evolution, while stability is found when the initial perturbation is damped at later times and eventually dies off. In fluid mechanics, such an approach is termed temporal analysis of the flow instability as opposed to the spatial analysis. The concept behind the spatial instability analysis is monitoring the evolution in space, typically in the streamwise direction, of a persistent perturbation signal localised at a given spatial position along the flow direction. Stated in these terms, the spatial instability interchanges the roles of space and time relatively to the temporal instability. The practical interest of the spatial instability analysis stems from the aerodynamics of jets in the pioneering papers by Betchov and Criminale Jr [7] and by Keller et al. [8]. In this field, the more recent papers by Alves et al. [9] and by Afzaal et al. [10] made significant contributions. An outlook into the research carried out in the area of spatially-developing instability can be found in textbooks such as Schmid and Henningson [11].

The aim of this paper is to survey the comparison between the temporal instability and the spatial instability by providing a simple example where such concepts are exploited: the one-dimensional Burgers’ equation with a driving linear force. The framework of spatial instability is applied to a real-world system made of a horizontal fluid-saturated porous layer with an anisotropic permeability bounded by impermeable horizontal walls kept at different uniform temperatures. This Rayleigh-Bénard system is studied along the steps illustrated in previous papers [12, 13] where the special case of an isotropic saturated porous layer has been analysed.

2 From the temporal analysis to the spatial analysis

The classical approach to the study of the linear instability in fluid mechanics, viz. the instability to small-amplitude perturbations of a given flow, consists in a direct application of Lyapunov’s idea: testing the instability of an equilibrium state, i.e., of a stationary solution of the governing equations of fluid flow means slightly altering the initial condition at time \(t=0\) and monitoring how this change modifies the system evolution at \(t>0\). In this framework, instability means a gradually amplifying discrepancy from the original stationary solution as time evolves. Furthermore, in this well-established approach, a perturbation is a disturbance of the initial condition at \(t=0\) with the boundary conditions for the system left unchanged.

Fig. 1
figure 1

The temporal analysis of instability versus the spatial analysis as a comparison between the effects of a time-evolving perturbation and the effects of a space-evolving perturbation (color figure online)

Full size image
Fig. 2
figure 2

The temporal analysis of instability through its two-fold implementation: the convective instability and the absolute instability (color figure online)

Full size image

Another approach to instability in fluid mechanics, the spatial analysis, marks a sharp difference with respect to the temporal analysis of time-evolving perturbations. In fact, the spatial analysis does not rely on Lyapunov’s idea of instability as a way to test the altered time evolution of a system as a consequence of a small change in the initial condition. If we imagine a flow system with a streamwise coordinate x, the spatial analysis is meant to test the development along the x direction of a perturbed boundary, or inlet, condition at \(x=0\). Then, the aim of the spatial analysis of instability is the determination of the effects produced downstream or upstream of a time-periodic perturbation signal, viz. a Fourier mode, set at a given position, conventionally at \(x=0\). A schematic illustration of the comparison between the analysis of time-evolving perturbations and the analysis of space-evolving perturbations is provided in Fig. 1.

Fig. 3
figure 3

The spatial analysis of instability (color figure online)

Full size image

There are two ways, traditionally employed, for the implementation of the temporal stability analysis and, hence, of the study of time-evolving perturbations: the modal analysis and the wavepacket analysis. The difference relies in the specification of the initial condition set at \(t=0\). The modal analysis arises from an initial perturbation expressed as a plane wave with wavenumber k along the x direction, namely a single Fourier mode. On the other hand, the wavepacket analysis is developed by assuming an initial perturbation given by a wavepacket, namely an envelope of infinite plane waves with all possible wavenumbers or, stated differently, a Fourier integral. A scheme of the two approaches to the temporal analysis is shown in Fig. 2. We mention that the wavepacket analysis leads to the definition of a parametric regime called absolute instability diverse from the convective instability defined via the modal analysis [14].

A sketch illustrating the spatial approach to instability and, hence, the study of space-evolving perturbations is displayed in Fig. 3. This figure shows the key point of the spatial instability where the response of the flow system is monitored downstream or upstream of a persistent time-periodic perturbation set up at a given spatial position, \(x=0\). The time-periodic signal at \(x=0\) is, in fact, a Fourier mode having an angular frequency \(\omega \). In Fig. 3, the already mentioned departure of the spatial instability concept from Lyapunov’s idea of instability is also highlighted.

Hereafter, we will focus on the comparison between the temporal analysis and the spatial analysis of instability in their modal formulations. In other terms, the focus will be on convective instability versus spatial instability. On the other hand, no further discussion of the wavepacket dynamics and the absolute instability will be provided here as this topic has been extensively presented in a recent book [14].

3 One-dimensional Burgers’ flow

Let us consider a flow system whose dynamics is governed by a one-dimensional Burgers’ equation with a driving linear force, namely

$$\begin{aligned} \frac{\partial {W}}{\partial {t}} + W \frac{\partial {W}}{\partial {x}} =\frac{\partial ^{2}W}{\partial {x}^{2}}+ \gamma \ ( W - a),\quad (x,t) \in {\mathbb {R}}^2, \end{aligned}$$
(1)

where \(\gamma \) and a are two real constants with \(a \ne 0\). The assumption \(a > 0\) is not restrictive as one can always recover the behaviour for a negative a by applying to Eq. (1) the transformation \(W \rightarrow - W\) with \(x \rightarrow -x\). Hence, in the following, we will implicitly consider \(a > 0\). Although Burgers’ equation is a partial differential equation originally formulated as a toy model for developing turbulence [15], we will not mind about its physical meaning and just use it as an arena for the implementation of the temporal and spatial stability analysis.

Equation (1) admits a simple constant solution,

$$\begin{aligned} W(x,t) = a, \end{aligned}$$
(2)

whose linear instability can be investigated by employing either the temporal analysis or the spatial analysis. In both cases, the starting point is the definition of a small-amplitude perturbation of the equilibrium solution given by Eq. (2), namely

$$\begin{aligned} W(x, t) = a + \varepsilon w(x,t), \end{aligned}$$
(3)

where \(\varepsilon \) is a perturbation amplitude parameter. Hence, one can substitute Eq. (3) into Eq. (1) and neglect terms \({\mathcal {O}}({\varepsilon ^2})\), so that we obtain

$$\begin{aligned} \frac{\partial {w}}{\partial {t}} + a \frac{\partial {w}}{\partial {x}} = \frac{\partial ^{2}w}{\partial {x}^{2}} + \gamma w. \end{aligned}$$
(4)

3.1 Temporal analysis

One may employ an x-based Fourier transform to solve Eq. (4), namely

$$\begin{aligned} {\tilde{w}}(k,t) = \frac{1}{\sqrt{2 \pi }} \int _{-\infty }^{\infty } w(x,t) e^{-i k x} \textrm{d}x,\quad w(x,t) = \frac{1}{\sqrt{2 \pi }} \int _{-\infty }^{\infty } {\tilde{w}}(k,t) e^{i k x} \textrm{d}k, \end{aligned}$$
(5)

where k has the meaning of a wavenumber. Due to the properties of the Fourier transform for derivatives with respect to x, Eqs. (4) and (5) yield

$$\begin{aligned} \frac{\partial {{\tilde{w}}}}{\partial {t}} = \lambda (k) {\tilde{w}}, \text { for } \lambda (k) = \gamma - k^2 - i k a. \end{aligned}$$
(6)

Here, the polynomial expression of \(\lambda (k)\) is the stability dispersion relation. Furthermore, from Eq. (6), we determine the Fourier transform of the perturbation w, namely

$$\begin{aligned} {\tilde{w}}(k,t) = {\tilde{w}}(k,0) e^{\lambda (k) t}. \end{aligned}$$
(7)

The analysis of convective instability is modal, which means that the growth or decay in time of the perturbation is assessed for the single Fourier mode which is periodic in the x coordinate with wavelength \(2\pi /k\), as sketched in Fig. 2. The exponential growth or decay of \({\tilde{w}}(k,t)\) is evidently regulated by the real part of \(\lambda (k)\), so that one predicts convective instability, on the basis of Eq. (6), for

$$\begin{aligned} \gamma > k^2. \end{aligned}$$
(8)

This condition is evidently independent of the constant a, while this constant determines the angular frequency \(\omega \) of the perturbation wave. In fact, by substituting the expression of \(\lambda (k)\) given by Eq. (6) into Eq. (7), each single Fourier mode in Eq. (5) contains the exponential with imaginary argument

$$\begin{aligned} e^{i(k x - \omega t)} \text { for } \omega = k a. \end{aligned}$$

Deciding whether the Fourier integral expressing w(x, t) unboundedly grows for large times t, at a fixed position x, is a different matter which does rely on the steepest-descent approximation of the integral [14]. The answer to this question leads to the condition of absolute instability. We refer the reader to Barletta [14] for details on this point. We just mention that satisfying Eq. (8) means having an unbounded growth of some Fourier modes, i.e. those with \(|k| < \sqrt{\gamma }\), but this does not imply in general an unbounded growth at large times of the Fourier integral, Eq. (5), expressing w(x, t).

3.2 Spatial analysis

One can use a t-based Fourier transform to solve Eq. (4), namely

$$\begin{aligned} {\hat{w}}(x, \omega ) = \frac{1}{\sqrt{2 \pi }} \int _{-\infty }^{\infty } w(x,t) e^{i \omega t} \textrm{d}t, \quad w(x,t) = \frac{1}{\sqrt{2 \pi }} \int _{-\infty }^{\infty } {\hat{w}}(x,\omega ) e^{-i \omega t} \textrm{d}\omega . \end{aligned}$$
(9)

On account of the properties of the Fourier transform for the derivative with respect to t, Eqs. (4) and (9) yield

$$\begin{aligned} \frac{\partial ^{2}{\hat{w}}}{\partial {x}^{2}} - a \frac{\partial {{\hat{w}}}}{\partial {x}} + \left( \gamma + i \omega \right) {\hat{w}} = 0. \end{aligned}$$
(10)

The solution of the linear ordinary differential equation at constant coefficients, Eq. (10), is given by

$$\begin{aligned} {\hat{w}}(x,\omega ) = c_{+}(\omega ) e^{\eta _{+}(\omega ) x} + c_{-}(\omega ) e^{\eta _{-}(\omega ) x}, \end{aligned}$$
(11)

where \(c_{+}(\omega )\) and \(c_{-}(\omega )\) are integration constants determined by the inlet conditions prescribed at \(x=0\), while \(\eta _{+}(\omega )\) and \(\eta _{-}(\omega )\) are the roots of the quadratic equation

$$\begin{aligned} \eta ^2 - a \eta + \gamma + i \omega = 0. \end{aligned}$$
(12)

Let us denote with s and k the real and the imaginary parts of \(\eta \), respectively,

$$\begin{aligned} \eta = s +i k. \end{aligned}$$
(13)

Thus, each Fourier mode in the integral representation of w(x, t), Eq. (9), contains the product of two exponential factors

$$\begin{aligned} e^{s x} e^{i(k x - \omega t)}. \end{aligned}$$
(14)

In fact, the product given by Eq. (14) involves an amplifying/damping factor, \(e^{s x}\). Along the positive x direction, we get amplification when \(s > 0\) and damping when \(s < 0\). The other factor in the product given by Eq. (14) is the same as that encountered on carrying out the temporal analysis. It defines a travelling wave along the x direction with phase velocity \(\omega /k\). Hence, the sign of \(\omega /k\) serves to determine whether the wave travels in the positive \((\omega /k > 0)\) or in the negative \((\omega /k < 0)\) direction of the x axis. The difference with respect to the convective instability is apparent. In the temporal framework, we have a one way evolution along the positive direction of t, due to the causality principle. In the spatial framework, the evolution can be in the positive x direction or in the negative x direction. Therefore, to assess the spatial instability of a Fourier mode, the sign of the spatial growth rate s is not sufficient, as we also need to know the direction of the x axis where the wave is heading and, hence, the sign of \(\omega /k\). One can define the spatial instability by checking if a given Fourier mode grows or decays in the direction of the x axis where this wave is travelling.

Definition 1

A Fourier mode with angular frequency \(\omega \) and \(k \ne 0\) is spatially unstable when

$$\begin{aligned} \frac{s \omega }{k} > 0, \end{aligned}$$
(15)

where s and k are real solutions of

$$\begin{aligned} {\left\{ \begin{array}{ll} s^2 - k^2 - a s + \gamma = 0 \\ 2 k s - a k + \omega = 0 \end{array}\right. }. \end{aligned}$$
(16)

It must be noted that Eq. (16) is obtained by substituting Eq. (13) into Eq. (12). Moreover, for every \(\omega \), there are two pairs (s, k) satisfying Eq. (16) as there are two complex roots of Eq. (12), i.e. \(\eta _{+}(\omega )\) and \(\eta _{-}(\omega )\). We also note that Eqs. (15) and (16) are left invariant by the transformation

$$\begin{aligned} (\omega , s, k) \rightarrow (-\omega , s, - k). \end{aligned}$$
(17)

Hence, the assumption \(\omega \ge 0\) is not restrictive as we can always use Eq. (17) to extend our considerations to negative angular frequencies. Thus, we can reformulate the condition for spatial instability expressed by Definition 1. In fact, one may focus on the case \(\omega \ge 0\) by recognising, from Eq. (13), that \(\textrm{Im }(\eta ^2) = 2 s k\), where \(\textrm{Im }\) denotes the imaginary part of a complex number.

Fig. 4
figure 4

Burgers’ flow: \(\textrm{Im }(\eta ^2)/a^2\) versus \(\omega /a^2\) for the roots \(\eta _{+}\) and \(\eta _{-}\) with different values of \(\gamma /a^2\): a ranging from 0 to 100 in steps of 10; b ranging from 0 to 0.5 in steps of 0.05. Both in a and in b, the plot colour continuously changes from cyan \((\gamma /a^2 = 0)\) to green. The black dots denote the zeros for \(\omega ^2 = \gamma a^2\) (color figure online)

Full size image

Remark 1

A Fourier mode with angular frequency \(\omega > 0\) is spatially unstable when

$$\begin{aligned} \textrm{Im }(\eta ^2) > 0, \end{aligned}$$
(18)

where \(\eta \) is a complex root of Eq. (12).

Figure 4 shows that only the \(\eta _{+}\) root of Eq. (12) may satisfy the condition for spatial instability given by Eq. (18). Indeed, zeros of \(\textrm{Im }(\eta ^2) = 2 s k\) are possible for \(k=0\) with \(s \ne 0\), which yields \(\omega =0\) on account of Eq. (16), or for \(s=0\) with \(k \ne 0\). By employing Eq. (16), one may infer that the latter case leads to \(\omega = a k\) and \(\omega ^2 = \gamma a^2\). Such values of \(\omega \) are denoted with black dots in Fig. 4. Frames (b) in Fig. 4 just evidence the small-\(\omega \) and small-\(\gamma \) parametric region, whereas frames (a) are relative to a significantly larger parametric region.

An important finding is that \(s=0\) defines the subset of the spatial Fourier modes involved in the integral transform given by Eq. (9) that is included in the set of temporal Fourier modes discussed in Sect. 3.1. More precisely, the spatial modes with \(s=0\) coincide with the neutrally stable temporal modes, i.e. those modes having both a zero spatial growth rate and a zero temporal growth rate.

Equation (16) is the starting point for determining the parametric condition for the onset of the spatial instability.

Theorem 1

Spatially unstable modes with \(\omega > 0\) can exist only if

$$\begin{aligned} \gamma > \frac{\omega ^2}{a^2}. \end{aligned}$$
(19)

Proof

For the sake of brevity, we use the notation

$$\begin{aligned} r = \textrm{Im }(\eta ^2) = 2 s k. \end{aligned}$$
(20)

Then, the second Eq. (16) yields

$$\begin{aligned} k = \frac{r + \omega }{a}. \end{aligned}$$
(21)

Substitution into the first Eq. (16) leads to

$$\begin{aligned} \gamma = \frac{a^4 r (r+2 \omega )+4 (r+\omega )^4}{4 a^2 (r+\omega )^2}. \end{aligned}$$
(22)

For every given a and \(\omega > 0\), the right hand side of Eq. (22) can be considered as a function of r. Then, we can define

$$\begin{aligned} Y(r) = \frac{a^4 r (r+2 \omega )+4 (r+\omega )^4}{4 a^2 (r+\omega )^2}, \end{aligned}$$
(23)

whose derivative is given by

$$\begin{aligned} Y'(r) = \frac{a^4 \omega ^2+4 (r+\omega )^4}{2 a^2 (r+\omega )^3}. \end{aligned}$$
(24)

Thus, Eq. (24) allows one to conclude that \(Y'(r) > 0\) for every \(r > - \omega \). As a consequence, Y(r) is a monotonic increasing function of r for \(r \ge 0\). Thus, on account of Eqs. (22) and (23), we conclude that

$$\begin{aligned} r> 0 \quad \Longrightarrow \quad \gamma = Y(r) > Y(0) = \frac{\omega ^2}{a^2}. \end{aligned}$$
(25)

\(\square \)

Theorem 1 reveals that, according to the linear analysis, the spatial instability region in the \((\omega , \gamma )\) plane, for \(\omega > 0\), is equivalent to the temporal instability region in the \((k, \gamma )\) plane. In fact, the instability region lies in every case above the neutral stability curve which is equivalently given by \(\gamma = \omega ^2/a^2\) or by \(\gamma = k^2\), as a consequence of the equality \(\omega = a k\) valid for \(s=0\). However, one must bear in mind that the temporal instability and the spatial instability manifest themselves in different manners as, in the former case, the perturbation modes grow exponentially in time and, in the latter case, the perturbation modes grow exponentially in space along their direction of propagation.

Frames (b) of Fig. 4 show that, when \(\omega \rightarrow 0\), \(\textrm{Im }(\eta ^2) \ne 0\) only if \(\gamma > a^2/4\). In particular, one has \(\textrm{Im }(\eta ^2) = \pm a (\gamma - a^2/4)^{1/2}\) with a spatial growth rate \(s = a/2\). A comment on the possibility to satisfy Eq. (16) with \(k = 0\), which yields \(\textrm{Im }(\eta ^2) = 0\), and \(\omega = 0\), is definitely important. Indeed, this case identifies a type of time-independent modes where the spatial growth rate s is a root of

$$\begin{aligned} s^2 - a s + \gamma = 0. \end{aligned}$$
(26)

Real roots of Eq. (26) exist only for \(\gamma \le a^2/4\),

$$\begin{aligned} s = \frac{a \pm \sqrt{a^2 - 4 \gamma }}{2}. \end{aligned}$$
(27)

With a positive \(\gamma \le a^2/4\), both roots of Eq. (26) are positive. They yield perturbation modes undergoing a purely exponential growth along the positive x direction. With a negative \(\gamma \), Eq. (27) yields a positive and a negative s meaning a growing perturbation mode along the positive x direction and a growing perturbation mode along the negative x direction. The exponential growth in |x| of such time-independent modes defines a growing departure from the basic equilibrium state as |x| increases. In this sense, the spatial modes with \(k=0\) and \(\omega =0\) can be classified as spatially unstable even if they do not satisfy Eq. (18), whereas the left hand side of Eq. (15) is actually undefined. If one accepts this conception of spatial instability, then spatially unstable modes may exist in a parametric domain \((\gamma < 0)\) where temporal instability is not possible according to a linear analysis.

4 Horizontal anisotropic porous layer

Let us consider a horizontal porous layer with thickness L and infinite horizontal width saturated by a Newtonian fluid. We assume a two-dimensional flow field in the (x, y) plane (see Fig. 5) with seepage velocity \({\textbf {u}} =(u_x, u_y)\) and temperature T. Heating is supplied from below through impermeable and isothermal boundaries kept at different uniform temperatures, \(T_1\) and \(T_2\).

The porous material has a uniform, but anisotropic, permeability with the permeability tensor having principal axes along the x and y directions,

$$\begin{aligned} {\textbf {K}} = \left( \begin{array}{cc} K_x &{} 0 \\ 0 &{} K_y \end{array} \right) . \end{aligned}$$
(28)
Fig. 5
figure 5

Fluid saturated porous layer: sketch of the coordinate system and boundary conditions (color figure online)

Full size image

4.1 The governing equations

The momentum transfer is modelled according to Darcy’s law [3, 6, 16] and the Boussinesq approximation is claimed in order to model the thermal buoyancy force induced by the non-uniform temperature field. In a dimensionless formulation, the system of local mass, momentum and energy balance equations is given by [3, 6, 16],

$$\begin{aligned}{} & {} \frac{\partial {u_{x}}}{\partial {x}} + \frac{\partial {u_{y}}}{\partial {y}} =0, \nonumber \\{} & {} \quad \frac{\partial {u_{x}}}{\partial {y}} - \tau \frac{\partial {u_{y}}}{\partial {x}} + R \frac{\partial {T}}{\partial {x}} = 0, \nonumber \\{} & {} \quad \frac{\partial {T}}{\partial {t}} + u_x \frac{\partial {T}}{\partial {x}} + u_y \frac{\partial {T}}{\partial {y}} = \frac{\partial ^{2}T}{\partial {x}^{2}} + \frac{\partial ^{2}T}{\partial {y}^{2}}, \end{aligned}$$
(29)

where the local momentum balance equation is written in its vorticity formulation, so that the dependence on the pressure field is encompassed. Furthermore, the permeability ratio \(\tau \) and the Rayleigh number R are defined as

$$\begin{aligned} \tau = \frac{K_x}{K_y}, \quad R = \frac{g \beta (T_1 - T_2) K_x L}{\nu \alpha }. \end{aligned}$$
(30)

Here, \(\alpha \), \(\beta \) and \(\nu \) are the average thermal diffusivity of the saturated porous medium, the thermal expansion coefficient of the fluid and the kinematic viscosity of the fluid, respectively. The gravitational acceleration \({\textbf {g}}\) has a modulus g. Due to their definitions, both parameters \(\tau \) and R are to be considered as positive.

In order to obtain Eq. (29), we denote with \(\sigma \) the ratio between the heat capacity of the saturated porous medium and that of the fluid, while the governing balance equations are made dimensionless by employing the constants

$$\begin{aligned} L, \quad \frac{L^2 \sigma }{\alpha }, \quad \frac{\alpha }{L}, \quad T_1 - T_2, \end{aligned}$$
(31)

to scale the coordinates, time, velocity and temperature, respectively. More precisely, the dimensionless temperature is defined as the ratio between \(T - T_2\) and the constant \(T_1 - T_2\).

The dimensionless governing equation (29) can be reformulated by employing the streamfunction \(\Psi \) defined as

$$\begin{aligned} u_x = \frac{\partial {\Psi }}{\partial {y}}, \quad u_y = - \frac{\partial {\Psi }}{\partial {x}}. \end{aligned}$$
(32)

In fact, we have

$$\begin{aligned}{} & {} \tau \frac{\partial ^{2}{\Psi }}{\partial {x}^{2}} + \frac{\partial ^{2}{\Psi }}{\partial {y}^{2}} + R \frac{\partial {T}}{\partial {x}} =0, \nonumber \\{} & {} \quad \frac{\partial {T}}{\partial {t}} + \frac{\partial {\Psi }}{\partial {y}} \frac{\partial {T}}{\partial {x}} - \frac{\partial {\Psi }}{\partial {x}} \frac{\partial {T}}{\partial {y}} = \frac{\partial ^{2}T}{\partial {x}^{2}} + \frac{\partial ^{2}T}{\partial {y}^{2}}. \end{aligned}$$
(33)

The boundary planes \(y=0\) and \(y=1\) are considered impermeable and isothermal, namely

$$\begin{aligned} \frac{\partial {\Psi }}{\partial {x}}= & {} 0, \quad T = 1 \quad \text {for} \quad y = 0, \nonumber \\ \frac{\partial {\Psi }}{\partial {x}}= & {} 0, \quad T= 0 \quad \text {for} \quad y = 1. \end{aligned}$$
(34)

4.2 Basic solution

A stationary solution of Eqs. (33) and (34) describing the basic equilibrium state is given by

$$\begin{aligned} \Psi = 0, \quad T = 1 - y. \end{aligned}$$
(35)

4.3 Linearised perturbation dynamics

We define the small-amplitude streamfunction and temperature perturbations of the stationary solution given by Eq. (35) as

$$\begin{aligned} \Psi (x,y,t) = \varepsilon \psi (x,y,t), \quad T(x,y,t) = 1 - y + \varepsilon \theta (x,y,t), \end{aligned}$$
(36)

where \(\varepsilon \) is the perturbation amplitude parameter. Linearisation is carried out by substituting Eq. (36) into Eqs. (33) and (34) and neglecting terms \({\mathcal {O}}({\varepsilon ^2})\). Thus, we can write

$$\begin{aligned}{} & {} \tau \frac{\partial ^{2}{\Psi }}{\partial {x}^{2}} + \frac{\partial ^{2}{\Psi }}{\partial {y}^{2}} + R \frac{\partial {\theta }}{\partial {x}} = 0, \nonumber \\{} & {} \quad \frac{\partial {\theta }}{\partial {t}} + \frac{\partial {\psi }}{\partial {x}} = \frac{\partial ^{2}{\theta }}{\partial {x}^{2}} + \frac{\partial ^{2}{\theta }}{\partial {y}^{2}}, \nonumber \\{} & {} \quad \frac{\partial {\psi }}{\partial {x}} = 0, \quad \theta = 0 \quad \text {for} \quad y = 0, 1. \end{aligned}$$
(37)

Fourier series in y can be employed, so that the solution of Eq. (37) is expressed as

$$\begin{aligned} \psi (x,y,t) \!\!=\!\! \sum _{n=1}^{\infty } \Phi (n,x,t) \sin \, (n \pi y), \quad \theta (x,y,t) \!\!=\!\! \sum _{n=1}^{\infty } \Theta (n,x,t) \sin \, (n \pi y). \end{aligned}$$
(38)

We can now rewrite Eq. (37) as

$$\begin{aligned}{} & {} \tau \frac{\partial ^{2}{\Phi }}{\partial {x}^{2}} - n^2\pi ^2 \Phi + R \frac{\partial {\Theta }}{\partial {x}} = 0, \nonumber \\{} & {} \quad \frac{\partial {\Theta }}{\partial {t}} + \frac{\partial {\Phi }}{\partial {x}} = \frac{\partial ^{2}{\Theta }}{\partial {x}^{2}} - n^2 \pi ^2 \Theta . \end{aligned}$$
(39)

4.4 Temporal analysis

We use the x-based Fourier transform defined by Eq. (5), so that Eq. (39) is transformed into

$$\begin{aligned}{} & {} (k^2 \tau + n^2\pi ^2) {\tilde{\Phi }} - i k R {\tilde{\Theta }} = 0, \nonumber \\{} & {} \quad \frac{\partial {{\tilde{\Theta }}}}{\partial {t}} + i k {\tilde{\Phi }} = -(k^2 + n^2 \pi ^2) {\tilde{\Theta }}. \end{aligned}$$
(40)

In particular, Eq. (40) yields

$$\begin{aligned} \frac{\partial {{\tilde{\Theta }}}}{\partial {t}} +\left( k^2 + n^2 \pi ^2 - \frac{k^2 R}{k^2 \tau + n^2\pi ^2} \right) {\tilde{\Theta }} = 0, \end{aligned}$$
(41)

whose solution is

$$\begin{aligned} {\tilde{\Theta }}(n,k,t) = {\tilde{\Theta }}(n,k,0) e^{\lambda (n, k) t} \quad \text {for} \quad \lambda (n, k) = \frac{k^2 R}{k^2 \tau + n^2\pi ^2} - k^2 - n^2 \pi ^2. \qquad \end{aligned}$$
(42)

The convective instability arises when the real part of \(\lambda (n, k)\) is positive, namely for

$$\begin{aligned} R > \frac{(k^2 \tau + n^2 \pi ^2)(k^2 + n^2 \pi ^2)}{k^2}. \end{aligned}$$
(43)

The absolute minimum of R for the onset of the convective instability is obtained with the \(n=1\) modes and with the critical values [17,18,19,20,21]

$$\begin{aligned} k = k_c = \frac{\pi }{\tau ^{{1}/{4}}}, \quad R = R_c = \pi ^2(1 + \sqrt{\tau }\,)^2. \end{aligned}$$
(44)

The neutral stability curve in the (k, R) plane is the plot of the function of k defined by the right hand side of Eq. (43) with \(n=1\). As such, it depends on the permeability ratio \(\tau \).

4.5 Spatial analysis

The alternative to the temporal analysis is the use of the t-based Fourier transform given by Eq. (9). By this method, Eq. (39) is transformed to

$$\begin{aligned}{} & {} \tau \frac{\partial ^{2}{\hat{\Phi }}}{\partial {x}^{2}} - n^2\pi ^2 \hat{\Phi } + R \frac{\partial {\hat{\Theta }}}{\partial {x}} = 0, \nonumber \\{} & {} \quad - i \omega \hat{\Theta } + \frac{\partial {\hat{\Phi }}}{\partial {x}} = \frac{\partial ^{2}{\hat{\Theta }}}{\partial {x}^{2}}- n^2 \pi ^2 \hat{\Theta }. \end{aligned}$$
(45)

Equation (45) is a system of ordinary differential equations in the independent variable x. Its solution can be expressed as

$$\begin{aligned}{} & {} \hat{\Phi }(n,x,\omega ) = R \sum _{j=1}^4 \frac{\eta _j(n,\omega ) C_j (n, \omega ) e^{\eta _j(n,\omega ) x}}{n^2 \pi ^2 - \tau [\eta _j(n,\omega )]^2}, \nonumber \\{} & {} \quad \hat{\Theta }(n,x,\omega ) = \sum _{j=1}^4 C_j (n, \omega ) e^{\eta _j(n,\omega ) x}, \end{aligned}$$
(46)

where \(\eta _j (n,\omega )\), with \(j = 1, \ldots , 4\), coincides with a root \(\eta \) of the fourth-degree equation

$$\begin{aligned} (n^2 \pi ^2 - \tau \eta ^2)(n^2 \pi ^2 - \eta ^2 - i \omega ) + R \eta ^2 = 0, \end{aligned}$$
(47)

while \(C_j (n, \omega )\) are coefficients to be determined on the basis of the conditions prescribed at \(x=0\) for the perturbations \(\psi \) and \(\theta \).

The modal stability analysis can be scaled down to the case \(n=1\) by recognising that the transformation defined by

$$\begin{aligned} \frac{\eta }{n} \rightarrow \eta , \quad \frac{\omega }{n^2} \rightarrow \omega , \quad \frac{R}{n^2} \rightarrow R, \end{aligned}$$
(48)

maps Eq. (47) into its \(n=1\) version, namely

$$\begin{aligned} (\pi ^2 - \tau \eta ^2)(\pi ^2 - \eta ^2 - i \omega ) + R \eta ^2 = 0. \end{aligned}$$
(49)

Hence, our forthcoming discussion will rely on Eq. (49).

Fig. 6
figure 6

Porous medium with \(\tau = 1\): \(\textrm{Im }(\eta ^2)\) versus \(\omega \) for the roots \(\eta _{+}^2\) and \(\eta _{-}^2\) with different values of R ranging from 0 to 50 in steps of 5 with the plot colour continuously changing from cyan \((R=0)\) to green (color figure online)

Full size image
Fig. 7
figure 7

Porous medium with \(\tau = 2/3\): \(\textrm{Im }(\eta ^2)\) versus \(\omega \) for the roots \(\eta _{+}^2\) and \(\eta _{-}^2\) with different values of R ranging from 0 to 50 in steps of 5 with the plot colour continuously changing from cyan \((R=0)\) to green (color figure online)

Full size image
Fig. 8
figure 8

Porous medium with \(\tau = 3/2\): \(\textrm{Im }(\eta ^2)\) versus \(\omega \) for the roots \(\eta _{+}^2\) and \(\eta _{-}^2\) with different values of R ranging from 0 to 50 in steps of 5 with the plot colour continuously changing from cyan \((R=0)\) to green (color figure online)

Full size image

We note that Eq. (49) is endowed with the same invariance defined, for Burger’s flow, by Eq. (17). This means that we can just focus on positive values of \(\omega \), exactly as for Burger’s flow.

Fig. 9
figure 9

Porous medium with \(\tau = 1\): \(\textrm{Im }(\eta ^2)\) versus \(\omega \) for the roots \(\eta _{+}^2\) and \(\eta _{-}^2\) with different values of R ranging from 0 to \(-50\) in steps of \(-5\) with the plot colour continuously changing from cyan \((R=0)\) to green (color figure online)

Full size image
Fig. 10
figure 10

Porous medium with \(\tau = 2/3\): \(\textrm{Im }(\eta ^2)\) versus \(\omega \) for the roots \(\eta _{+}^2\) and \(\eta _{-}^2\) with different values of R ranging from 0 to \(-50\) in steps of \(-5\) with the plot colour continuously changing from cyan \((R=0)\) to green (color figure online)

Full size image
Fig. 11
figure 11

Porous medium with \(\tau = 3/2\): \(\textrm{Im }(\eta ^2)\) versus \(\omega \) for the roots \(\eta _{+}^2\) and \(\eta _{-}^2\) with different values of R ranging from 0 to \(-50\) in steps of \(-5\) with the plot colour continuously changing from cyan \((R=0)\) to green (color figure online)

Full size image

Figures 6, 7 and 8 display plots of \(\textrm{Im }(\eta ^2)\) versus \(\omega \) for the two roots of Eq. (49) and several values of R. These figures are relative to different anisotropy ratios \(\tau \), with Fig. 6 describing the isotropy limiting case \((\tau = 1)\). The isotropic medium is a special case already investigated by Barletta [12]. Figure 7 is for \(\tau = 2/3\), a case where the porous medium is more permeable in the vertical direction than in the horizontal direction. Figure 8 is for \(\tau = 3/2\), i.e. the porous medium is more permeable in the horizontal direction than in the vertical direction. In Figs. 6, 7 and 8, the cyan curves denote the limiting case where \(R \rightarrow 0\). As it is easily gathered from Eq. (49), there are two possible roots for \(R \rightarrow 0\): one is \(\eta ^2 = \pi ^2/\tau \) leading to \(\textrm{Im }(\eta ^2) = 0\), the other root yields \(\textrm{Im }(\eta ^2) = - \omega \). The former root is not acceptable as it would lead to a vanishing denominator in Eq. (46). Then, the corresponding plot is reported just as a reference case in Figs. 6, 7 and 8.

The most important information disclosed by Figs. 6, 7 and 8 is that, in every case, spatially unstable modes exist for every \(R > 0\). In fact, as pointed out in Remark 1, spatially unstable modes are identified by the condition \(\textrm{Im }(\eta ^2) > 0\). Having in mind Eqs. (43) and (44), which provide a statement of the condition for temporal instability and reveal that such a condition is allowed only when \(R > R_c = \pi ^2 (1 + \sqrt{\tau }\,)^2\), Figs. 6, 7 and 8 show that one may have \(\textrm{Im }(\eta ^2) > 0\) also for \(0< R < R_c\). We have reached the conclusion that the spatial instability is not equivalent to the temporal instability.

The analysis carried out so far revealed that every setup having \(R > 0\), may display modes satisfying the condition \(\textrm{Im }(\eta ^2) > 0\) with \(\omega > 0\) and, hence, the condition of spatial instability is satisfied for every case with heating from below, no matter how small is R. We also observed that \(R =0\) yields either \(\textrm{Im }(\eta ^2) = 0\) or \(\textrm{Im }(\eta ^2) = - \omega \) and, thus, it does not show up any spatial instability. Heating from above, i.e. \(R < 0\), is a parametric condition explicitly excluded in our initial description as we claimed \(T_1 > T_2\) (see Fig. 5), but considering \(R < 0\) may offer some further insights anyway. In fact, Figs. 9, 10 and 11 display the plots of \(\textrm{Im }(\eta ^2)\) versus \(\omega \) for \(R \le 0\) ranging from 0 to \(-50\). These figures are relative to the isotropic case \(\tau = 1\) (Fig. 9), \(\tau = 2/3\) (Fig. 10) and \(\tau = 3/2\) (Fig. 11). Exactly as for Figs. 6-8, all Figs. 9-11 contain two frames: one for each root, either \(\eta _{+}^2\) or \(\eta _{-}^2\), of Eq. (49). For all the three cases, \(\tau = 1\), \(\tau = 2/3\) and \(\tau = 3/2\), devised in Figs. 9, 10 and 11, we have always \(\textrm{Im }(\eta ^2) \le 0\) suggesting that no spatial instability is possible. However, as pointed out in Sect. 3 for the analysis of Burgers’ flow, one has to be quite careful in the evaluation of the limiting case of stationary modes, \(\omega \rightarrow 0\).

It must be emphasised that, on account of Eq. (49), the limit for \(\omega \rightarrow 0\) of \(\textrm{Im }(\eta ^2)\) can be nonzero only if

$$\begin{aligned} \pi ^2(1 - \sqrt{\tau }\,)^2< R < \pi ^2(1 + \sqrt{\tau }\,)^2, \end{aligned}$$
(50)

where the upper bound of the interval is precisely the critical value \(R_c\) given by Eq. (44). Consistently with Eq. (50), Fig. 6 displays \(\textrm{Im }(\eta ^2) \rightarrow 0\) when \(\omega \rightarrow 0\) only for \(R = 0\) and \(R \ge R_c \approx 39.5\), namely for \(R = 0\), 40, 45 and 50. Similarly, Fig. 7 displays \(\textrm{Im }(\eta ^2) \rightarrow 0\) only for \(R = 0\) and \(R \ge R_c \approx 32.6\), namely for \(R = 0\), 35, 40, 45 and 50, while Fig. 8 displays \(\textrm{Im }(\eta ^2) = 0\) only for \(R = 0\) and \(R \ge R_c \approx 48.8\), that is for \(R = 0\) and 50. On the other hand, for all cases reported in Figs. 9, 10 and 11 relative to \(R \le 0\), one has \(\textrm{Im }(\eta ^2) \rightarrow 0\) when \(\omega \rightarrow 0\). This is expected as a negative R always lies outside the interval defined by Eq. (50). Following the same reasoning presented for the analysis of Burgers’ flow, we may consider the cases where the limit \(\omega \rightarrow 0\) leads to \(k=0\), being s and k the real part and the imaginary part of \(\eta \), as reported in Eq. (13). This situation undoubtedly means that the limit \(\omega \rightarrow 0\) yields \(\textrm{Im }(\eta ^2) = 2 k s \rightarrow 0\). Hence, spatial modes with both \(\omega =0\) and \(k=0\) exist only when R is outside the interval defined by Eq. (50). By employing Eq. (47), one can determine \(s^2\) for such modes as

$$\begin{aligned} s^2 = \frac{-[R - \pi ^2(1 + \tau )] \pm \sqrt{[R - \pi ^2(1 + \tau )]^2 - 4 \pi ^4 \tau }}{2 \tau }. \end{aligned}$$
(51)

The right hand side of Eq. (51) is real and non-negative, with both \(+\) and −, provided that

$$\begin{aligned} R \le \pi ^2(1 - \sqrt{\tau })^2. \end{aligned}$$
(52)

Therefore, we can have spatial modes with \(\omega =0\) and \(k=0\) with either a positive or a negative s whenever Eq. (52) is satisfied. Such modes are stationary, have \(\textrm{Im }(\eta ^2) = 0\), but they display an exponential growth in space either in the positive or in the negative x direction. Exactly as for Burgers’ flow, we can consider these modes as a manifestation of spatial instability. The important point is that they do exist also for conditions of heating from above \((R < 0)\). Figure 12 displays a scheme of the different instability types in the \((\tau , R)\) parametric plane. The parametric region where growing spatial modes with \(\omega =0\) and \(k=0\) exist alongside spatially unstable modes with \(\omega > 0\) and \(\textrm{Im }(\eta ^2) > 0\) is labelled as “spatial instability with \(\omega \ge 0\)” for brevity. The parametric region \(R < 0\) where growing spatial modes with \(\omega =0\) and \(k=0\) exist is labelled as “spatial instability with \(\omega = 0\)”.

Fig. 12
figure 12

Porous medium: instability diagram in the \((\tau , R)\) parametric plane (color figure online)

Full size image

5 The scope of the spatial instability analysis

The study of two test cases, a toy model based on Burgers’ equation and a real-world system modelled as a fluid saturated porous medium with anisotropic permeability, have shown that the linear analysis of the perturbations within either a temporal instability concept or a spatial instability concept may lead to diverse scenarios. The temporal instability, in its modal formulation, allows one to evaluate the convective instability threshold bounding the unstable parametric region through the neutral stability curve and the critical values of the governing parameters. The spatial analysis, in its modal formulation, revealed conditions of linear instability for cases where the temporal analysis predicts linear stability. In the examples discussed so far, we found that the spatial instability may always exist either with the Fourier spatial modes having \(\omega \ne 0\) or, at the very least, with the stationary spatial modes obtained in the limit \(\omega \rightarrow 0\). Then, a question arises: does the linear analysis of spatial perturbation modes always lead to the prediction of instability no matter which values are prescribed for the parameters governing the system? The answer is negative. In fact, a neat example is provided by the one-dimensional diffusion equation,

$$\begin{aligned} \frac{\partial {u}}{\partial {t}} = \frac{\partial ^{2}{u}}{\partial {x}^{2}}. \end{aligned}$$
(53)

An equilibrium state, \(u = u_0\), is any stationary solution of Eq. (53) and, hence, any linear function of x. If we introduce a perturbation of the equilibrium state, \(u = u_0 + \varepsilon U\), the linear nature of Eq. (53) leads to a governing equation for the perturbation function U which coincides with Eq. (53) where u is now replaced by U. The convective instability analysis, carried out following the steps and the notations introduced in Sect. 3.1, leads us to the temporal dispersion relation,

$$\begin{aligned} \lambda (k) = - k^2. \end{aligned}$$
(54)

The well-known conclusion is utterly evident: every equilibrium state \(u_0\) is stable as the real part of \(\lambda \), and hence the time growth-rate, is negative for every \(k \ne 0\). The spatial instability analysis can be developed according to the steps and notations introduced in Sect. 3.2, leading us to the spatial dispersion relation,

$$\begin{aligned}{}[\eta (\omega )]^2 = - i \omega . \end{aligned}$$
(55)

We have \(\textrm{Im }(\eta ^2) = - \omega < 0\) for \(\omega > 0\) so that, by adopting the criterion defined in Remark 1, again we reach the conclusion that every equilibrium state \(u_0\) is stable. Extending the spatial instability analysis to the limit of stationary modes, \(\omega \rightarrow 0\), does not lead to any further insights. In fact, in this limit, Eq. (55) just yields \(k=0\) and \(s=0\).

The diffusion equation example provides an illustration of the selective nature of the spatial instability analysis. Despite what we concluded for the two test cases developed in Sects. 3 and 4, the spatial instability condition is not satisfied in every case.

6 Conclusions and final remarks

The onset of the linear instability in a flow system has been outlined by focussing on two diverse schemes termed temporal instability and spatial instability. If the former scheme is an exploitation of Lyapunov’s idea of instability in a mechanical system, the latter scheme is a an utterly different approach where the roles of space and time are interchanged on defining the evolution variable for the development of the instability. Indeed, Lyapunov’s or temporal instability can be described as the amplification in time of an initial perturbation superposed at \(t = 0\) to the equilibrium state. Spatial instability can be described as the amplification in the x spatial direction of an inlet perturbation localised at \(x = 0\). Within a modal approach, the temporal instability analysis tests the reaction of the system to space-periodic perturbations with a given wavenumber k, while spatial instability analysis tests the reaction of the system to time-periodic perturbations with a given angular frequency \(\omega \).

The comparison between the temporal instability analysis and the spatial instability analysis has been introduced via a simple flow system governed by a one-dimensional Burgers’ equation with a driving linear force. The condition of spatial instability has been formulated through an inequality where the imaginary part of the square complex growth rate \(\eta \) must be positive for positive values of \(\omega \). Such an inequality reflects the condition of spatial instability as a spatial amplitude growth of the modal perturbation along the direction where the perturbation wave propagates. The study of Burger’s flow leads to the conclusion that the domains of temporal instability and spatial instability are parametrically equivalent.

It has been demonstrated that the equivalence between temporal instability and spatial instability is rather an exception than the rule, by developing the analysis of a Rayleigh-Bénard system where an anisotropic porous layer saturated by a fluid is heated from below. A permeability anisotropy ratio \(\tau \) and the Rayleigh number R are the governing parameters for this flow system. It has been shown that, while the temporal instability can be started up only if the Rayleigh number exceeds its \(\tau \)-dependent critical value \(R_c\), the spatial instability may exist not only when \(R > R_c\), but also when \(0 < R \le R_c\). This conclusion rules out, in this case, the equivalence between temporal instability and spatial instability proved for the case of Burgers’ flow.

The special case where the Fourier modes employed for the spatial stability analysis have both a zero angular frequency and a zero wavenumber has been considered. Such modes are time-independent with an undefined direction of propagation. For the simple example of Burgers’ flow, it has been shown that this class of modes may involve a nonzero spatial growth rate even in a parametric domain where the equilibrium solution is both temporally and spatially stable. However, these modes actually display an exponential growth in the positive or in the negative spatial direction. Such a growth can be intended as a form of spatial instability even if not associated to any direction of propagation. The same conclusion has been found also for the Rayleigh-Bénard system in a saturated anisotropic porous layer. In the study of the saturated porous layer, these stationary Fourier modes with zero wavenumber yet growing in space exist also when the Rayleigh number is negative, i.e. for conditions of heating from above.

The temporal instability may be studied in a non-modal formulation by testing the time evolution of a perturbation wavepacket expressed via a Fourier integral. This variant of the temporal analysis leads to the definition of a typically supercritical condition called absolute instability [14]. In principle, a similar non-modal conception could be conceived also for the spatial instability. Bringing the spatial instability beyond the purely modal analysis, as well as beyond the merely linear formulation, could offer interesting opportunities for future research in the fluid mechanics of saturated porous media.