2.1 Governing equations

Now, the mass conservation equation and Cauchy’s momentum equation in the fluid layer are given as

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\\ \displaystyle \unicode[STIX]{x1D70C}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}.\end{array}\right\}\end{eqnarray}$$

Equation (2.1) is simplified to obtain the $x$ -momentum equation as

(2.2) $$\begin{eqnarray}0=-\!\left(\frac{\text{d}p}{\text{d}x}\right)+\frac{\text{d}\unicode[STIX]{x1D70F}}{\text{d}z}.\end{eqnarray}$$

In the above equation, $\unicode[STIX]{x1D70F}$ is the shear stress of the fluid.

The constitutive relation for a Bingham fluid is expressed as

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}ll@{}}\displaystyle \unicode[STIX]{x1D70F}=\unicode[STIX]{x1D707}\frac{\text{d}u}{\text{d}z}+\unicode[STIX]{x1D70F}_{0}\,\text{sgn}\left(\frac{\text{d}u}{\text{d}z}\right), & |\unicode[STIX]{x1D70F}|\geqslant \unicode[STIX]{x1D70F}_{0},\\ \displaystyle \frac{\text{d}u}{\text{d}z}=0, & |\unicode[STIX]{x1D70F}|<\unicode[STIX]{x1D70F}_{0},\end{array}\right\}\end{eqnarray}$$

where $u=u(z)$ is the streamwise velocity component, $\unicode[STIX]{x1D70F}_{0}$ and $\unicode[STIX]{x1D707}$ represent the yield stress and the plastic viscosity of the Bingham fluid, and $\text{sgn}$ denotes the signum function.

For characterizing flow in the anisotropic porous medium, the governing equations are

(2.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{m}=0,\\ \displaystyle \left(\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D712}}\right)\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{m}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}p_{m}-\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D646}}\boldsymbol{u}_{m}.\end{array}\right\}\end{eqnarray}$$

In the above equation, $\boldsymbol{u}_{m}$ is the velocity vector in the porous layer, and $p_{m}$ is the interstitial pressure. Also, $\unicode[STIX]{x1D712}$ is the porosity, and $\unicode[STIX]{x1D646}$ is a diagonal tensor used to represent the anisotropic and inhomogeneous permeability of the porous medium. This tensor is expressed as $\unicode[STIX]{x1D646}=K_{x}\unicode[STIX]{x1D702}_{x}(z/d_{m})\boldsymbol{i}\boldsymbol{i}+K_{y}\unicode[STIX]{x1D702}_{y}(z/d_{m})\boldsymbol{j}\boldsymbol{j}+K_{z}\unicode[STIX]{x1D702}_{z}(z/d_{m})\boldsymbol{k}\boldsymbol{k}$ , where $\boldsymbol{i}$ , $\boldsymbol{j}$ and $\boldsymbol{k}$ are the unit vectors in the $x$ , $y$ and $z$ directions, respectively. In this expression, $K_{x}$ , $K_{y}$ and $K_{z}$ represent the permeabilities in the $x$ , $y$ and $z$ directions; and $\unicode[STIX]{x1D702}_{x}$ , $\unicode[STIX]{x1D702}_{y}$ and $\unicode[STIX]{x1D702}_{z}$ denote the respective inhomogeneity functions. A generalized formulation would allow each of these inhomogeneity functions to be multivariate. However, in the present study, $\unicode[STIX]{x1D702}_{x}$ , $\unicode[STIX]{x1D702}_{y}$ and $\unicode[STIX]{x1D702}_{z}$ have been assumed to be functions of $z$ only. This is done in order to ensure that the base solution is not multidimensional in nature. In this study, an exponential variation of the inhomogeneity function is assumed, $\unicode[STIX]{x1D702}_{x}=\unicode[STIX]{x1D702}_{y}=\unicode[STIX]{x1D702}_{z}=\exp (A_{inh}(1+\bar{z}_{m}))$ , where $A_{inh}$ is the inhomogeneity factor (Chen & Hsu Reference Chen and Hsu1991), and $\bar{z}_{m}$ is the dimensionless distance in the porous layer (as defined later in (2.10)). Continuity of pressure is assumed at the fluid–porous interface, i.e.  $p=p_{m}$ (Beavers & Joseph Reference Beavers and Joseph1967).

The velocity profile in the porous layer is given by the plane Buckingham–Reiner model (Rees Reference Rees and Vafai2015). In the present notation of variables, it is expressed as

(2.5) $$\begin{eqnarray}u_{m}=-\frac{K_{x}\unicode[STIX]{x1D702}_{x}}{\unicode[STIX]{x1D707}}\left[1-\frac{3}{2}\left(\frac{\unicode[STIX]{x1D6FC}_{0}}{\displaystyle \left|\frac{\text{d}p}{\text{d}x}\right|}\right)+\frac{1}{2}\left(\frac{\unicode[STIX]{x1D6FC}_{0}}{\displaystyle \left|\frac{\text{d}p}{\text{d}x}\right|}\right)^{\!3}\right]\left(\frac{\text{d}p}{\text{d}x}\right),\quad \text{for }\left|\frac{\text{d}p}{\text{d}x}\right|>\unicode[STIX]{x1D6FC}_{0},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{0}=(\unicode[STIX]{x1D6FD}_{0}\unicode[STIX]{x1D70F}_{0})/\sqrt{K_{x}}$ is essentially the threshold gradient, i.e. the limiting pressure gradient for flow to occur, and $\unicode[STIX]{x1D6FD}_{0}$  is a dimensionless parameter. When the applied pressure gradient is unable to exceed the threshold gradient, there is no flow in the porous layer. It may be observed that the velocity in the porous layer continuously approaches zero, as the limiting pressure gradient is attained. This is a characteristic feature of Bingham flow in porous media, as also confirmed by pore-scale studies (Balhoff & Thompson Reference Balhoff and Thompson2004; Nash & Rees Reference Nash and Rees2017).

The authors recognize the prevalence of the Darcy–Brinkman model for describing porous layer flow. Brinkman (Reference Brinkman1949) suggested the incorporation of a viscous stress component in the Darcy equation. The study describes an ‘effective viscosity’ to be used in the viscous stress term. However, it seems that there is a lack of consensus in the literature regarding the definition of this effective viscosity (Liu & Liu Reference Liu and Liu2009). In addition, it appears that most researchers have agreed upon the fact that the Brinkman model is applicable for high-porosity configurations (Nield & Bejan Reference Nield and Bejan2006; Auriault Reference Auriault2009). A detailed Stokesian dynamics study for flow in porous media was carried out by Durlofsky & Brady (Reference Durlofsky and Brady1987), and they concluded that Brinkman’s equation should be used for porosity greater than 0.95. On the contrary, in our analysis, we considered a porosity of only 0.1, for which the validity of the Brinkman model is undoubtedly debatable. Typical application of a viscoplastic fluid flow over a porous layer is in oil recovery, where drilling fluid flows over porous soil. The porosity of major varieties of soil lies in the range of 0.1–0.4 (Association of Swiss Road and Traffic Engineers 1999; Das Reference Das2013). Other researchers who have considered the porosity in a similar range (0.1–0.4) have also neglected the Brinkman model for the similar reason of a low-porosity configuration (Chang et al. Reference Chang, Chen and Straughan2006; Deepu et al. Reference Deepu, Anand and Basu2015, Reference Deepu, Kallurkar, Anand and Basu2016; Chang et al. Reference Chang, Chen and Chang2017). Only researchers who have carried out the analysis for a ‘highly porous material’ have adopted the Brinkman model (Hill & Straughan Reference Hill and Straughan2008, Reference Hill and Straughan2009). These are the major factors why the Brinkman model is not considered to describe the flow in porous medium in this study. Incorporation of the Brinkman model can be suitably done in a future study, exploring the stability of a flow configuration consisting of a highly porous medium.

2.2 Boundary conditions

At the upper wall, the no-slip boundary condition is assumed, i.e.

(2.6a,b ) $$\begin{eqnarray}\text{at }z=d,\quad u=0.\end{eqnarray}$$

The fluid–porous interface $(z=0)$ is described by the Beavers–Joseph boundary condition (Beavers & Joseph Reference Beavers and Joseph1967), i.e.

(2.7a,b ) $$\begin{eqnarray}\text{at }z=0,\quad \frac{\text{d}u}{\text{d}z}=\frac{\unicode[STIX]{x1D6FC}_{BJ}}{\sqrt{K_{x}\unicode[STIX]{x1D702}_{x}(0)}}(u-u_{m}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{BJ}$ is the Beavers–Joseph constant. For a Bingham fluid, the velocity remains invariant within the plug zone. In addition, both the velocity and its gradient are assumed to be continuous at both yield surfaces. Therefore,

(2.8) $$\begin{eqnarray}u|_{z=z_{01}}=u|_{z=z_{02}},\end{eqnarray}$$

and

(2.9) $$\begin{eqnarray}\left.\frac{\text{d}u}{\text{d}z}\right|_{z=z_{01}}=\left.\frac{\text{d}u}{\text{d}z}\right|_{z=z_{02}}=0.\end{eqnarray}$$

2.3 Non-dimensionalization

The variables are made non-dimensional as follows:

(2.10a-d ) $$\begin{eqnarray}\bar{z}=\frac{z}{d},\quad \bar{z}_{m}=\frac{z}{d_{m}},\quad \bar{u}=\frac{u}{U_{p}},\quad \bar{\unicode[STIX]{x1D70F}}=\frac{\unicode[STIX]{x1D70F}d}{\unicode[STIX]{x1D707}U_{p}},\end{eqnarray}$$

where $U_{p}=(d^{2}/\unicode[STIX]{x1D707})(-\text{d}p/\text{d}x)$ .

2.3.1 Non-dimensional governing equations

Equation (2.2) may be represented in terms of non-dimensional variables as

(2.11a ) $$\begin{eqnarray}1+\frac{\text{d}\bar{\unicode[STIX]{x1D70F}}}{\text{d}\bar{z}}=0,\end{eqnarray}$$

i.e.

(2.11b ) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D70F}}=-\bar{z}+l_{1},\end{eqnarray}$$

where $l_{1}$ is the integration constant.

The constitutive relation for a Bingham fluid (2.3) may be represented in dimensionless form as

(2.12) $$\begin{eqnarray}\left.\begin{array}{@{}ll@{}}\displaystyle \bar{\unicode[STIX]{x1D70F}}=\frac{\text{d}\bar{u}}{\text{d}\bar{z}}+Bn\,\text{sgn}\left(\frac{\text{d}\bar{u}}{\text{d}\bar{z}}\right), & |\bar{\unicode[STIX]{x1D70F}}|\geqslant Bn,\\ \displaystyle \frac{\text{d}\bar{u}}{\text{d}\bar{z}}=0, & |\bar{\unicode[STIX]{x1D70F}}|<Bn.\end{array}\right\}\end{eqnarray}$$

In the above equation, $Bn=(\unicode[STIX]{x1D70F}_{0}d)/(\unicode[STIX]{x1D707}U_{p})$ is the Bingham number of the fluid.

2.3.2 Non-dimensional boundary conditions

The boundary conditions (2.6)–(2.9) may be recast in terms of dimensionless variables:

(2.13a,b ) $$\begin{eqnarray}\text{at }\bar{z}=1,\quad \bar{u}=0,\end{eqnarray}$$

(2.14a,b ) $$\begin{eqnarray}\text{at }\bar{z}=0,\quad \frac{\text{d}\bar{u}}{\text{d}\bar{z}}=\frac{\unicode[STIX]{x1D6FC}_{BJ}\hat{d}}{\unicode[STIX]{x1D6FF}\sqrt{\unicode[STIX]{x1D702}_{x}(0)}}(\bar{u}-\bar{u}_{m}),\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{u}|_{\bar{z}=\bar{z}_{01}}=\bar{u}|_{\bar{z}=\bar{z}_{02}}, & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\text{d}\bar{u}}{\text{d}\bar{z}}\right|_{\bar{z}=\bar{z}_{01}}=\left.\frac{\text{d}\bar{u}}{\text{d}\bar{z}}\right|_{\bar{z}=\bar{z}_{02}}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}=\sqrt{K_{x}}/d_{m}$ represents the Darcy number, and $\hat{d}=d/d_{m}$ is the ratio of the thicknesses of the fluid and the porous layers; $\hat{d}$ is referred to as the depth ratio of the channel.

2.4 Base flow – analytical solution

The velocity profile may be derived by solving (2.11) subject to boundary conditions (2.13)–(2.16). For the region $0\leqslant \bar{z}\leqslant \bar{z}_{01}$ (i.e. the lower shear zone), $\text{d}\bar{u}/\text{d}\bar{z}>0$ , i.e. the velocity gradient is positive. The velocity may be given as

(2.17a ) $$\begin{eqnarray}\bar{u}=\left(l_{1}\bar{z}-\frac{\bar{z}^{2}}{2}-Bn\,\bar{z}\right)+l_{2}.\end{eqnarray}$$

In the upper shear zone (i.e.  $\bar{z}_{02}\leqslant \bar{z}\leqslant 1$ ), $\text{d}\bar{u}/\text{d}\bar{z}<0$ . The velocity is expressed as

(2.17b ) $$\begin{eqnarray}\bar{u}=\left(l_{1}\bar{z}-\frac{\bar{z}^{2}}{2}+Bn\,\bar{z}\right)+l_{3},\end{eqnarray}$$

where $l_{2}$ and $l_{3}$ are the integration constants.

For the intermediate unyielded zone $(\bar{z}_{01}\leqslant \bar{z}\leqslant \bar{z}_{02})$ , the velocity is constant, given by $\bar{u}|_{\bar{z}=\bar{z}_{01}}=\bar{u}|_{\bar{z}=\bar{z}_{02}}$ .

2.4.1 Velocity profile in the fluid layer

The dimensionless base velocity in the fluid layer may be given as ( $\bar{v}$ and $\bar{w}$ represent the dimensionless spanwise and normal velocity components)

(2.18) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \begin{array}{@{}rcl@{}}\displaystyle \bar{u}(\bar{z})\, & = & \displaystyle \left(l_{1}\bar{z}-\frac{\bar{z}^{2}}{2}-Bn\,\bar{z}\right)+l_{2},\quad 0\leqslant \bar{z}\leqslant \bar{z}_{01},\\ \displaystyle \, & = & \displaystyle \frac{\bar{z}_{01}^{2}}{2}+l_{2},\quad \bar{z}_{01}\leqslant \bar{z}\leqslant \bar{z}_{02},\\ \, & = & \displaystyle \left(l_{1}\bar{z}-\frac{\bar{z}^{2}}{2}+Bn\,\bar{z}\right)+l_{3},\quad \bar{z}_{02}\leqslant \bar{z}\leqslant 1,\end{array}\\ \bar{v}=0,\quad \bar{w}=0,\quad 0\leqslant \bar{z}\leqslant 1,\end{array}\right\}\end{eqnarray}$$

where

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle l_{1}=Bn+\frac{\unicode[STIX]{x1D6FC}_{BJ}\hat{d}}{\unicode[STIX]{x1D6FF}\sqrt{\unicode[STIX]{x1D702}_{x}(0)}}[l_{2}-\bar{u}_{m}(0)], & \displaystyle\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle & \displaystyle l_{2}=\left[\frac{\displaystyle 1-4Bn(1-Bn)+\frac{2\bar{u}_{m}(0)\unicode[STIX]{x1D6FC}_{BJ}\hat{d}(1-2Bn)}{\unicode[STIX]{x1D6FF}\sqrt{\unicode[STIX]{x1D702}_{x}(0)}}}{\displaystyle 2\left(1+\frac{\unicode[STIX]{x1D6FC}_{BJ}\hat{d}(1-2Bn)}{\unicode[STIX]{x1D6FF}\sqrt{\unicode[STIX]{x1D702}_{x}(0)}}\right)}\right], & \displaystyle\end{eqnarray}$$

(2.21) $$\begin{eqnarray}\displaystyle & \displaystyle l_{3}=\frac{1}{2}-2Bn-\frac{\unicode[STIX]{x1D6FC}_{BJ}\hat{d}}{\unicode[STIX]{x1D6FF}\sqrt{\unicode[STIX]{x1D702}_{x}(0)}}[l_{2}-\bar{u}_{m}(0)]. & \displaystyle\end{eqnarray}$$

The locations of the yield surfaces at $\bar{z}_{01}$ and $\bar{z}_{02}$ are expressed as

(2.22) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{z}_{01}=\frac{\unicode[STIX]{x1D6FC}_{BJ}\hat{d}}{\unicode[STIX]{x1D6FF}\sqrt{\unicode[STIX]{x1D702}_{x}(0)}}[l_{2}-\bar{u}_{m}(0)], & \displaystyle\end{eqnarray}$$

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{z}_{02}=2Bn+\frac{\unicode[STIX]{x1D6FC}_{BJ}\hat{d}}{\unicode[STIX]{x1D6FF}\sqrt{\unicode[STIX]{x1D702}_{x}(0)}}[l_{2}-\bar{u}_{m}(0)]. & \displaystyle\end{eqnarray}$$

Using (2.22) and (2.23), it is observed that

(2.24) $$\begin{eqnarray}\bar{z}_{02}=2Bn+\bar{z}_{01}.\end{eqnarray}$$

It is evident that the plug region should be within the fluid volume. Therefore, the following constraint is imposed on $\bar{z}_{01}$ and $\bar{z}_{02}$ :

(2.25a,b ) $$\begin{eqnarray}0<\bar{z}_{01}<1,\quad 0<\bar{z}_{02}<1.\end{eqnarray}$$

In addition, it may be inferred from (2.24) and (2.25) that

(2.26) $$\begin{eqnarray}Bn<0.5.\end{eqnarray}$$

2.4.2 Velocity profile in the porous layer

The non-dimensional velocity profile in the porous layer is obtained as

(2.27) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\bar{u}_{m}(\bar{z}_{m})=l_{4}\unicode[STIX]{x1D702}_{x}(\bar{z}_{m}),\quad -1\leqslant \bar{z}_{m}\leqslant 0,\\ \bar{v}_{m}=0,\quad \bar{w}_{m}=0,\quad -1\leqslant \bar{z}_{m}\leqslant 0,\end{array}\right\},\quad \text{for }\frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FD}_{0}\hat{d}}>Bn,\end{eqnarray}$$

where

(2.28) $$\begin{eqnarray}l_{4}=\frac{\unicode[STIX]{x1D6FF}^{2}}{\hat{d}^{2}}\left[1-\frac{3}{2}\left(\frac{\unicode[STIX]{x1D6FD}_{0}Bn\,\hat{d}}{\unicode[STIX]{x1D6FF}}\right)+\frac{1}{2}\left(\frac{\unicode[STIX]{x1D6FD}_{0}Bn\,\hat{d}}{\unicode[STIX]{x1D6FF}}\right)^{\!3}\right].\end{eqnarray}$$

When the threshold gradient is not reached (i.e. when $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FD}_{0}\hat{d}<Bn$ ), the velocity is zero in the porous layer.

2.5 Stability analysis

In earlier studies on stability analysis of viscoplastic fluids, there was no existence of the porous layer. As a result, the flow profile was symmetric about the middle of the plug zone. This eliminated the need to consider both the shear zones in the analysis. Stability analysis for the decoupled problem was carried out only for the upper shear zone. Here, it may be noted that the flow configuration consists of an upper shear zone, a plug zone, a lower shear zone and a porous layer. The complexity of the current study is enhanced by the fact that the velocity profile is asymmetric about the middle of the plug zone. There are two flow domains, separated by a plug zone discontinuity. One is the domain I comprising the porous layer and the lower shear zone, extending from $z=-d_{m}$ to $z_{01}$ . The other is the domain II, consisting of the upper shear zone, extending from $z=z_{02}$ to $d$ . Fortunately, it may be noted that the maximum velocity is the same for both domains, and is realized at the interface of the shear and plug zones. It is evident that the governing perturbation equations will be the same for both the shear zones. Since the maximum realizable velocity is the same in either of the two shear zones, for the current study, we carry out the stability analysis for the domain I. In other words, we study the zone extending from $z=-d_{m}$ to $z_{01}$ (comprising the porous layer and the lower shear zone). Henceforth, $z_{01}$ is represented as $z_{0}$ for the sake of simplicity.

2.5.1 Perturbation equations

Now, we allow a perturbation to the system in the form

(2.29) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u=\overline{u}+u^{\prime },\quad v=\overline{v}+v^{\prime },\quad w=\overline{w}+w^{\prime },\\ u_{m}=\overline{u_{m}}+u_{m}^{\prime },\quad v_{m}=\overline{v_{m}}+v_{m}^{\prime },\quad w_{m}=\overline{w_{m}}+w_{m}^{\prime },\\ p=\overline{p}+p^{\prime },\quad z_{0}=\overline{z_{0}}+h^{\prime },\end{array}\right\}\end{eqnarray}$$

where $\bar{u}$ and $u^{\prime }$ respectively denote $x$ -direction steady-state velocity, and the velocity perturbation. The rest of the variables also have analogous meanings.

The perturbation equations for the fluid layer are given as ( $\text{D}=\text{d}/\text{d}z$ ; in addition, we write $\bar{u}=U$ for convenience in the following equations)

(2.30) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}v^{\prime }}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}w^{\prime }}{\unicode[STIX]{x2202}z}=0,\end{eqnarray}$$

(2.31) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}t}+U\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}x}+w^{\prime }\text{D}U & = & \displaystyle -\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}x}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}u^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{Bn}{Re}\frac{1}{|\text{D}U|}\left(\unicode[STIX]{x1D6FB}^{2}u^{\prime }-\frac{\unicode[STIX]{x2202}^{2}w^{\prime }}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}-\frac{\unicode[STIX]{x2202}^{2}u^{\prime }}{\unicode[STIX]{x2202}z^{2}}\right),\end{eqnarray}$$

(2.32) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}v^{\prime }}{\unicode[STIX]{x2202}t}+U\frac{\unicode[STIX]{x2202}v^{\prime }}{\unicode[STIX]{x2202}x} & = & \displaystyle -\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}y}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}v^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{Bn}{Re}\left[\frac{\unicode[STIX]{x1D6FB}^{2}v^{\prime }}{|\text{D}U|}+\left(\frac{\unicode[STIX]{x2202}v^{\prime }}{\unicode[STIX]{x2202}z}+\frac{\unicode[STIX]{x2202}w^{\prime }}{\unicode[STIX]{x2202}y}\right)\left\{\frac{\text{d}}{\text{d}z}\left(\frac{1}{|\text{D}U|}\right)\right\}\right],\end{eqnarray}$$

(2.33) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}w^{\prime }}{\unicode[STIX]{x2202}t}+U\frac{\unicode[STIX]{x2202}w^{\prime }}{\unicode[STIX]{x2202}x} & = & \displaystyle -\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}z}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}w^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{Bn}{Re}\left[2\left(\frac{\unicode[STIX]{x2202}w^{\prime }}{\unicode[STIX]{x2202}z}\right)\left\{\frac{\text{d}}{\text{d}z}\left(\frac{1}{|\text{D}U|}\right)\right\}+\frac{\displaystyle \unicode[STIX]{x1D6FB}^{2}w^{\prime }-\frac{\unicode[STIX]{x2202}^{2}w^{\prime }}{\unicode[STIX]{x2202}x^{2}}-\frac{\unicode[STIX]{x2202}^{2}u^{\prime }}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}}{|\text{D}U|}\right].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The perturbation equations for the porous layer are expressed as

(2.34) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{m}^{\prime }}{\unicode[STIX]{x2202}x_{m}}+\frac{\unicode[STIX]{x2202}v_{m}^{\prime }}{\unicode[STIX]{x2202}y_{m}}+\frac{\unicode[STIX]{x2202}w_{m}^{\prime }}{\unicode[STIX]{x2202}z_{m}}=0, & \displaystyle\end{eqnarray}$$

(2.35) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{m}^{\prime }}{\unicode[STIX]{x2202}t_{m}}=-\left(\frac{\unicode[STIX]{x1D712}}{Re_{m}}\right)\left[\frac{\unicode[STIX]{x2202}p_{m}^{\prime }}{\unicode[STIX]{x2202}x_{m}}+\frac{u_{m}^{\prime }}{\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D702}_{x}}\right], & \displaystyle\end{eqnarray}$$

(2.36) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}v_{m}^{\prime }}{\unicode[STIX]{x2202}t_{m}}=-\left(\frac{\unicode[STIX]{x1D712}}{Re_{m}}\right)\left[\frac{\unicode[STIX]{x2202}p_{m}^{\prime }}{\unicode[STIX]{x2202}y_{m}}+\frac{\unicode[STIX]{x1D709}_{1}v_{m}^{\prime }}{\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D702}_{y}}\right], & \displaystyle\end{eqnarray}$$

(2.37) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}w_{m}^{\prime }}{\unicode[STIX]{x2202}t_{m}}=-\left(\frac{\unicode[STIX]{x1D712}}{Re_{m}}\right)\left[\frac{\unicode[STIX]{x2202}p_{m}^{\prime }}{\unicode[STIX]{x2202}z_{m}}+\frac{\unicode[STIX]{x1D709}_{2}w_{m}^{\prime }}{\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D702}_{z}}\right]. & \displaystyle\end{eqnarray}$$

It is known that Squire’s transformation is not valid for Bingham fluids (Frigaard et al. Reference Frigaard, Howison and Sobey1994; Métivier & Nouar Reference Métivier and Nouar2011). Thus, to perform stability analysis, we assume three-dimensional disturbances of the form

(2.38) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u^{\prime }=\hat{u} (z,t)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y)},\quad v^{\prime }=\hat{v}(z,t)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y)},\quad w^{\prime }={\hat{w}}(z,t)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y)},\\ u_{m}^{\prime }=\hat{u} _{m}(z_{m},t_{m})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}_{m}x_{m}+\unicode[STIX]{x1D6FD}_{m}y_{m})},\quad v_{m}^{\prime }=\hat{v}_{m}(z_{m},t_{m})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}_{m}x_{m}+\unicode[STIX]{x1D6FD}_{m}y_{m})},\\ w_{m}^{\prime }={\hat{w}}_{m}(z_{m},t_{m})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}_{m}x_{m}+\unicode[STIX]{x1D6FD}_{m}y_{m})},\\ p^{\prime }=\hat{p}(z,t)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y)},\quad p_{m}^{\prime }=\hat{p}_{m}(z_{m},t_{m})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y)},\quad h^{\prime }={\hat{h}}(t)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y)},\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are the real streamwise and spanwise wavenumbers, respectively. It is to be noted that yield stress fluids are unique in the sense that the yield surface also undergoes perturbation. Hence, perturbation in the yield surface is also being considered. Incorporation of the above disturbance expressions (2.38) in the governing perturbation equations (2.30)–(2.37) results in an initial value problem (the hats are dropped from $u$ , $v$ , etc., for convenience). Thus, in the following equations, $u$ is essentially $\hat{u}$ . Also, $\mathbb{Q}=\text{D}^{2}-(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})$ . In addition, we write $\bar{u}=U$ for convenience in the following equations, wherever applicable.

For the fluid layer, the following set of equations is obtained:

(2.39) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\text{i}(\unicode[STIX]{x1D6FC}u+\unicode[STIX]{x1D6FD}v)+\text{D}w, & \displaystyle\end{eqnarray}$$

(2.40) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}t}=-\text{i}\unicode[STIX]{x1D6FC}Uu-w\text{D}U-\text{i}\unicode[STIX]{x1D6FC}p+\frac{1}{Re}\mathbb{Q}u+\frac{Bn}{Re}\left[\frac{-(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})u-\text{i}\unicode[STIX]{x1D6FC}\text{D}w}{|\text{D}U|}\right], & \displaystyle\end{eqnarray}$$

(2.41) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}t} & = & \displaystyle -\text{i}\unicode[STIX]{x1D6FC}Uv-\text{i}\unicode[STIX]{x1D6FD}p+\frac{1}{Re}\mathbb{Q}v\nonumber\\ \displaystyle & & \displaystyle +\,\frac{Bn}{Re}\left[\frac{1}{|\text{D}U|}\mathbb{Q}v+(\text{D}v+\text{i}\unicode[STIX]{x1D6FD}w)\left\{\text{D}\left(\frac{1}{|\text{D}U|}\right)\right\}\right],\end{eqnarray}$$

(2.42) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}t} & = & \displaystyle -\text{i}\unicode[STIX]{x1D6FC}Uw-\text{D}p+\frac{1}{Re}\mathbb{Q}w\nonumber\\ \displaystyle & & \displaystyle +\,\frac{Bn}{Re}\left[(2\text{D}w)\left\{\text{D}\left(\frac{1}{|\text{D}U|}\right)\right\}+\frac{(\text{D}^{2}-\unicode[STIX]{x1D6FD}^{2})w-\text{i}\unicode[STIX]{x1D6FC}\text{D}u}{|\text{D}U|}\right].\end{eqnarray}$$

For the porous layer, the following initial value problem is obtained:

(2.43) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}(\unicode[STIX]{x1D6FC}_{m}u_{m}+\unicode[STIX]{x1D6FD}_{m}v_{m})+\text{D}_{p}w_{m}=0, & \displaystyle\end{eqnarray}$$

(2.44) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{m}}{\unicode[STIX]{x2202}t_{m}}=-\left(\frac{\unicode[STIX]{x1D712}}{Re_{m}}\right)\left[\frac{1}{\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D702}_{x}}u_{m}+\text{i}\unicode[STIX]{x1D6FC}_{m}p_{m}\right], & \displaystyle\end{eqnarray}$$

(2.45) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}v_{m}}{\unicode[STIX]{x2202}t_{m}}=-\left(\frac{\unicode[STIX]{x1D712}}{Re_{m}}\right)\left[\frac{\unicode[STIX]{x1D709}_{1}}{\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D702}_{y}}v_{m}+\text{i}\unicode[STIX]{x1D6FD}_{m}p_{m}\right], & \displaystyle\end{eqnarray}$$

(2.46) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}w_{m}}{\unicode[STIX]{x2202}t_{m}}=-\left(\frac{\unicode[STIX]{x1D712}}{Re_{m}}\right)\left[\frac{\unicode[STIX]{x1D709}_{2}}{\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D702}_{z}}w_{m}+\text{D}_{p}p_{m}\right], & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D709}_{1}=K_{x}/K_{y}$ and $\unicode[STIX]{x1D709}_{2}=K_{x}/K_{z}$ respectively denote the spanwise and normal anisotropy parameters. Also, $\text{D}_{p}=\text{d}/\text{d}z_{m}$ . The Reynolds numbers in the two layers are related as

(2.47) $$\begin{eqnarray}Re=\hat{d}Re_{m}\left(\frac{\displaystyle l_{2}+\frac{\bar{z}_{01}^{2}}{2}}{l_{4}\unicode[STIX]{x1D702}_{x}(0)}\right).\end{eqnarray}$$

For modal analysis, the disturbance expressions in the governing perturbation equations (2.30)–(2.37) are replaced, using the following relations:

(2.48) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u^{\prime }=\hat{u} (z)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y-ct)},\quad v^{\prime }=\hat{v}(z)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y-ct)},\quad w^{\prime }={\hat{w}}(z)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y-ct)},\\ u_{m}^{\prime }=\hat{u} _{m}(z_{m})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}_{m}x_{m}+\unicode[STIX]{x1D6FD}_{m}y_{m}-c_{m}t_{m})},\quad v_{m}^{\prime }=\hat{v}_{m}(z_{m})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}_{m}x_{m}+\unicode[STIX]{x1D6FD}_{m}y_{m}-c_{m}t_{m})},\\ w_{m}^{\prime }={\hat{w}}_{m}(z_{m})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}_{m}x_{m}+\unicode[STIX]{x1D6FD}_{m}y_{m}-c_{m}t_{m})},\\ p^{\prime }=\hat{p}(z)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y-ct)},\quad z_{0}^{\prime }=\hat{z}_{0}(z)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y-ct)},\\ p_{m}^{\prime }=\hat{p}_{m}(z_{m})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}_{m}x_{m}+\unicode[STIX]{x1D6FD}_{m}y_{m}-c_{m}t_{m})},\quad h^{\prime }={\hat{h}}(t)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y)}.\end{array}\right\}\end{eqnarray}$$

Thus, the following eigenvalue problem is obtained in the fluid and the porous layers:

(2.49) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\text{i}(\unicode[STIX]{x1D6FC}u+\unicode[STIX]{x1D6FD}v)+\text{D}w, & \displaystyle\end{eqnarray}$$

(2.50) $$\begin{eqnarray}\displaystyle & \displaystyle -\text{i}cu=-\text{i}\unicode[STIX]{x1D6FC}Uu-w\text{D}U-\text{i}\unicode[STIX]{x1D6FC}p+\frac{1}{Re}\mathbb{Q}u+\frac{Bn}{Re}\left[\frac{-(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})u-\text{i}\unicode[STIX]{x1D6FC}\text{D}w}{|\text{D}U|}\right], & \displaystyle\end{eqnarray}$$

(2.51) $$\begin{eqnarray}\displaystyle -\text{i}cv & = & \displaystyle -\text{i}\unicode[STIX]{x1D6FC}Uv-\text{i}\unicode[STIX]{x1D6FD}p+\frac{1}{Re}\mathbb{Q}v\nonumber\\ \displaystyle & & \displaystyle +\,\frac{Bn}{Re}\left[\frac{1}{|\text{D}U|}\mathbb{Q}v+(\text{D}v+\text{i}\unicode[STIX]{x1D6FD}w)\left\{\text{D}\left(\frac{1}{|\text{D}U|}\right)\right\}\right],\end{eqnarray}$$

(2.52) $$\begin{eqnarray}\displaystyle -\text{i}cw & = & \displaystyle -\text{i}\unicode[STIX]{x1D6FC}Uw-\text{D}p+\frac{1}{Re}\mathbb{Q}w\nonumber\\ \displaystyle & & \displaystyle +\,\frac{Bn}{Re}\left[(2\text{D}w)\left\{\text{D}\left(\frac{1}{|\text{D}U|}\right)\right\}+\frac{(\text{D}^{2}-\unicode[STIX]{x1D6FD}^{2})w-\text{i}\unicode[STIX]{x1D6FC}\text{D}u}{|\text{D}U|}\right],\end{eqnarray}$$

(2.53) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}(\unicode[STIX]{x1D6FC}_{m}u_{m}+\unicode[STIX]{x1D6FD}_{m}v_{m})+\text{D}_{p}w_{m}=0, & \displaystyle\end{eqnarray}$$

(2.54) $$\begin{eqnarray}\displaystyle & \displaystyle -\text{i}c_{m}u_{m}=-\left(\frac{\unicode[STIX]{x1D712}}{Re_{m}}\right)\left[\frac{1}{\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D702}_{x}}u_{m}+\text{i}\unicode[STIX]{x1D6FC}_{m}p_{m}\right], & \displaystyle\end{eqnarray}$$

(2.55) $$\begin{eqnarray}\displaystyle & \displaystyle -\text{i}c_{m}v_{m}=-\left(\frac{\unicode[STIX]{x1D712}}{Re_{m}}\right)\left[\frac{\unicode[STIX]{x1D709}_{1}}{\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D702}_{y}}v_{m}+\text{i}\unicode[STIX]{x1D6FD}_{m}p_{m}\right], & \displaystyle\end{eqnarray}$$

(2.56) $$\begin{eqnarray}\displaystyle & \displaystyle -\text{i}c_{m}w_{m}=-\left(\frac{\unicode[STIX]{x1D712}}{Re_{m}}\right)\left[\frac{\unicode[STIX]{x1D709}_{2}}{\unicode[STIX]{x1D6FF}^{2}\unicode[STIX]{x1D702}_{z}}w_{m}+\text{D}_{p}p_{m}\right]. & \displaystyle\end{eqnarray}$$

At the yield surface (interface of yielded and unyielded zone), the following boundary conditions prevail:

(2.57) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u(z_{0})=0,\quad v(z_{0})=0,\quad w(z_{0})=0,\\ \text{D}u(z_{0})=-h\text{D}^{2}U(z_{0}),\quad \text{D}v(z_{0})=0,\quad \text{D}w(z_{0})=0.\end{array}\right\}\end{eqnarray}$$

Other boundary conditions are obtained from the Beavers–Joseph condition at the fluid–porous interface, along with the continuity of normal and spanwise velocity, and pressure.

Thus, at the fluid–porous interface $(\bar{z}=\bar{z}_{m}=0)$ , we have

(2.58) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \text{D}u(0)=\frac{\unicode[STIX]{x1D6FC}_{BJ}\hat{d}}{\unicode[STIX]{x1D6FF}\sqrt{\unicode[STIX]{x1D702}_{x}(0)}}\left[u(0)-\frac{\hat{d}Re_{m}}{Re}u_{m}(0)\right],\\ \displaystyle Re\,v(0)=\hat{d}Re_{m}v_{m}(0),\\ Re\,w(0)=\hat{d}Re_{m}w_{m}(0),\\ \displaystyle p(0)=\frac{\hat{d}^{2}}{Re}Re_{m}p_{m}(0).\end{array}\right\}\end{eqnarray}$$

At the bottom surface of the porous layer ( $z=-d_{m}$ , or $\bar{z}_{m}=-1$ ), we have

(2.59a,b ) $$\begin{eqnarray}v_{m}(-1)=0,\quad w_{m}(-1)=0.\end{eqnarray}$$

The following relation exists among the wavenumbers $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FC}_{m},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FD}_{m})$ in the fluid and porous layers:

(2.60a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}=\hat{d}\unicode[STIX]{x1D6FC}_{m},\quad \unicode[STIX]{x1D6FD}=\hat{d}\unicode[STIX]{x1D6FD}_{m}.\end{eqnarray}$$

Equations (2.49)–(2.56) are second order in $u$ , $v$ , $w$ , and first order in $w_{m}$ , $p$ , $p_{m}$ . For solving the same, we are seemingly using 12 boundary conditions (2.57)–(2.59). In other words, the system is apparently over-specified. However, a careful observation will reveal the reality that is typical of yield stress fluids. In (2.57), the relation for $\text{D}u$ essentially specifies the condition for the amplitude of the perturbation of the yield surface. Thus, the relation basically defines $h$ , and not $u$ . The relations for $\text{D}v$ and $\text{D}w$ are also required, as they take care of the singularity of $|\text{D}U|^{-1}$ as the yield surface is approached. The same has been elucidated by other researchers as well (Frigaard & Nouar Reference Frigaard and Nouar2003; Nouar et al. Reference Nouar, Kabouya, Dusek and Mamou2007).

It may be highlighted here that the disturbance waves arising in the fluid and porous layers are the same, i.e. they have the same wavelength. As a result, the dimensional wavenumbers, being the reciprocal of wavelength, are also the same in the two layers. However, in our analysis, we are considering non-dimensional wavenumbers, namely $\unicode[STIX]{x1D6FC}$ in the fluid layer and $\unicode[STIX]{x1D6FC}_{m}$ in the porous layer: $\unicode[STIX]{x1D6FC}$ is non-dimensionalized with respect to $d$ (thickness of the fluid layer), while $\unicode[STIX]{x1D6FC}_{m}$ is made non-dimensional with respect to $d_{m}$ (thickness of the porous layer). Thus, $\unicode[STIX]{x1D6FC}\neq \unicode[STIX]{x1D6FC}_{m}$ (in fact, $\unicode[STIX]{x1D6FC}=\hat{d}\unicode[STIX]{x1D6FC}_{m}$ , as already stated; in fact, $\unicode[STIX]{x1D6FD}=\hat{d}\unicode[STIX]{x1D6FD}_{m}$ as well). In a similar way, the non-dimensional distances in the two layers are related as $\hat{d}x=x_{m}$ and $\hat{d}y=y_{m}$ . Thus, the Fourier expressions are identical in the two layers, i.e.  $\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y-ct)}=\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}_{m}x_{m}+\unicode[STIX]{x1D6FD}_{m}y_{m}-c_{m}t_{m})}$ . The modal expressions $u$ and $u^{\prime }$ differ in the $z$ direction, which is logical, as the base profiles are also different in the two layers. These assumptions are in conformity with earlier studies (Chang et al. Reference Chang, Chen and Straughan2006; Hill & Straughan Reference Hill and Straughan2008; Deepu et al. Reference Deepu, Anand and Basu2015; Chang et al. Reference Chang, Chen and Chang2017).

2.6 Numerical solution

We have assumed porosity $\unicode[STIX]{x1D712}=0.1$ throughout. Unless otherwise explicitly specified, the parameters assume the following values: $Bn=0.2$ , $\unicode[STIX]{x1D6FD}_{0}=0.016$ , $\hat{d}=0.3$ , $\unicode[STIX]{x1D6FC}_{BJ}=0.2$ , $\unicode[STIX]{x1D6FF}=0.001$ , $\unicode[STIX]{x1D709}_{1}=2$ , $\unicode[STIX]{x1D709}_{2}=2$ and $A_{inh}=2$ . For modal analysis, the resultant generalized eigenvalue problem of the form $\unicode[STIX]{x1D63C}\unicode[STIX]{x1D653}=c\unicode[STIX]{x1D63D}\unicode[STIX]{x1D653}$ is solved via the Chebyshev spectral collocation method using QZ decomposition. Here, $\unicode[STIX]{x1D653}=(u,v,w,p)$ , and $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}$ are the matrices comprising Chebyshev derivatives. The discretization of the equations is carried out on the Gauss–Lobatto grid, consisting of $(N+1)$ collocation points. The details of the spectral collocation method using Chebyshev polynomials have been discussed by Schmid & Henningson (Reference Schmid and Henningson2001). MATLAB 2014b software was used for the computation. The calculations were checked with increasing number of Chebyshev polynomials ( $N=32$ , 64, 96, 128, 192, 256, 384, 512). It was observed that for $N=192$ , the first 50 eigenvalues (sorted in terms of increasing imaginary part) gave accuracy up to five digits. In other words, the first five digits remained invariant with further increase in $N$ . Hence, all the computations reported in this paper were carried out with $N=192$ .

2.7 Non-modal analysis and energy growth

For non-modal analysis, the following linear initial value problem needs to be solved:

(2.61) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\boldsymbol{q}=-\text{i}\unicode[STIX]{x1D647}\boldsymbol{q},\end{eqnarray}$$

where $\boldsymbol{q}$ is expressed as $\boldsymbol{q}=(\hat{u} ,\hat{v},{\hat{w}},\hat{p})^{\text{T}}$ , and $\unicode[STIX]{x1D647}$ is defined as $\unicode[STIX]{x1D647}=-\text{i}\unicode[STIX]{x1D63D}^{-1}\unicode[STIX]{x1D63C}$ . The solution is expressed as $\boldsymbol{q}(t)=\boldsymbol{q}(0)\exp (-\text{i}\unicode[STIX]{x1D647}t)$ , where $\boldsymbol{q}(0)$ is the initial value of the perturbation variables.

2.7.1 Response to initial conditions: growth function

In the above context, the growth function defines the maximum possible amplification of the initial condition, and is mathematically expressed as

(2.62) $$\begin{eqnarray}G(t)=\max _{q(0)\neq 0}\frac{\Vert \boldsymbol{q}(t)\Vert ^{2}}{\Vert \boldsymbol{q}(0)\Vert ^{2}}=\Vert \exp (-\text{i}\unicode[STIX]{x1D647}t)\Vert ^{2},\end{eqnarray}$$

where $\Vert \cdot \Vert$ is the norm.

In general, for a flow stability problem, if the Reynolds number $Re$ is less than the critical energy Reynolds number $Re_{1}$ , then the growth function $G(t)$ is never greater than unity, i.e.  $G_{max}=1$ . In such a case, optimal time $t_{opt}$ is zero. On the other hand, if the Reynolds number is greater than the critical Reynolds number $Re_{2}$ , then one of the eigenvalues of the linear Orr–Sommerfeld operator has a positive imaginary part, and the flow becomes unstable. Mathematically, the growth function becomes unbounded at infinite time. In the case of $Re_{1}<Re<Re_{2}$ , modal analysis yields no unstable eigenvalue and thus the flow is linearly stable. However, a transient growth is observed which decays with time. For the current work, stability analysis is performed for a vast range of Reynolds numbers, but no unstable eigenvalues are obtained. Thus, for analysing non-modal behaviour, a Reynolds number higher than the critical energy Reynolds number $Re_{1}$ is taken into consideration.

2.7.2 Response to external excitations: resolvent

Let us assume that a flow configuration is subjected to an input signal $V$ , represented as

(2.63) $$\begin{eqnarray}V(x,y,z,t)=\exp (-\text{i}\unicode[STIX]{x1D714}t)v(x,y,z),\end{eqnarray}$$

where $\unicode[STIX]{x1D714}$ is the frequency of $V$ , and $\unicode[STIX]{x1D714}$ is assumed to be complex. Then, the relationship among the input signal $V$ and the response $U$ is given by the following differential equation:

(2.64) $$\begin{eqnarray}\frac{\text{d}U}{\text{d}t}=-\text{i}\unicode[STIX]{x1D647}U+\exp (-\text{i}\unicode[STIX]{x1D714}t)v.\end{eqnarray}$$

Solution of the above equation leads to the following expression for $U$ :

(2.65) $$\begin{eqnarray}U(x,y,z,t)=\text{i}\exp (-\text{i}\unicode[STIX]{x1D714}t)u=\text{i}\exp (-\text{i}\unicode[STIX]{x1D714}t)(\unicode[STIX]{x1D714}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D647})^{-1}v.\end{eqnarray}$$

In the above relation, $\unicode[STIX]{x1D644}$ is the identity matrix, and $(\unicode[STIX]{x1D714}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D647})^{-1}$ is referred to as the resolvent; $\unicode[STIX]{x1D647}$ has already been defined earlier. For an input signal of frequency $\unicode[STIX]{x1D714}$ , the maximum amplification of the external excitation, $R(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})$ , is given by the norm of the resolvent operator. Thus,

(2.66) $$\begin{eqnarray}R(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D714})=\max _{v\neq 0}\frac{\Vert u\Vert }{\Vert v\Vert }=\Vert (\unicode[STIX]{x1D714}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D647})^{-1}\Vert ,\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the identity matrix, and $\unicode[STIX]{x1D714}$ is the eigenvalue of $\unicode[STIX]{x1D647}$ , such that $\Vert (\unicode[STIX]{x1D714}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D647})^{-1}\Vert \rightarrow \infty$ . In mathematical terms, the resolvent norm is equivalent to amplification of the response under external forcing. In this context, it is essential to discuss the definition of the ‘ $\unicode[STIX]{x1D700}$ -pseudospectrum’. A basic understanding of the pseudospectrum is critical to non-modal stability analysis. For any $\unicode[STIX]{x1D700}\geqslant 0$ , the pseudospectrum of $\unicode[STIX]{x1D647}$ is defined as (Trefethen & Embree Reference Trefethen and Embree2005)

(2.67) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D647})=\{\unicode[STIX]{x1D714}\in \mathbb{C}:\Vert (\unicode[STIX]{x1D714}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D647})^{-1}\Vert \geqslant \unicode[STIX]{x1D700}^{-1}\}.\end{eqnarray}$$

With a variation in $\unicode[STIX]{x1D700}$ , contours of pseudospectra are obtained. When $\unicode[STIX]{x1D700}=0$ , the pseudospectrum reduces to the spectrum, i.e. the set of eigenvalues obtained from normal mode analysis. Thus, with increasing $\unicode[STIX]{x1D700}$ , there is deviation from the traditional eigenvalue analysis. In a non-modal stability analysis problem, there are two major components that need to be explored. One is the response of the system to external excitations, the other being the transient growth of the initial conditions. The resolvent norm describes the maximum possible amplification in response to external excitations. On the other hand, transient growth of the initial conditions is expressed by the growth function $G(t)$ . In the current study, we have considered only the response to initial conditions, i.e.  $G(t)$ . The pseudospectrum has been studied only in the limiting case of Newtonian fluids, in order to validate our simulation methodology with that of the existing literature.

2.8 Energy stability

For carrying out energy stability analysis, the perturbation kinetic energy is quantified as

(2.68) $$\begin{eqnarray}\displaystyle E & = & \displaystyle \frac{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}{8\unicode[STIX]{x03C0}^{2}}\left[\frac{1}{z_{0}}\int _{0}^{z_{0}}\int _{0}^{2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}}\int _{0}^{2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FD}}(u_{r}^{\prime 2}+v_{r}^{\prime 2}+w_{r}^{\prime 2})\,\text{d}y\,\text{d}x\,\text{d}z\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x1D6FC}_{m}\unicode[STIX]{x1D6FD}_{m}}{8\unicode[STIX]{x03C0}^{2}}\left[\int _{-1}^{0}\int _{0}^{2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}_{m}}\int _{0}^{2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FD}_{m}}(u_{mr}^{\prime 2}+v_{mr}^{\prime 2}+w_{mr}^{\prime 2})\,\text{d}y_{m}\,\text{d}x_{m}\,\text{d}z_{m}\right],\end{eqnarray}$$

where $u_{r}^{\prime }=\text{Real}(\hat{u} (z,t)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x+\unicode[STIX]{x1D6FD}y)})$ , $u_{mr}^{\prime }=\text{Real}(\hat{u} _{m}(z_{m},t_{m})\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}_{m}x_{m}+\unicode[STIX]{x1D6FD}_{m}y_{m})})$ , and so on. No energy growth is possible if $\text{d}E/\text{d}t<0$ as $t\rightarrow \infty$ (Butler & Farrell Reference Butler and Farrell1992). Thus, the condition for no energy growth may be expressed as

(2.69) $$\begin{eqnarray}\frac{1}{Re_{1}}=\sup _{\boldsymbol{u}}\frac{\unicode[STIX]{x1D644}_{1}(\boldsymbol{u})}{\unicode[STIX]{x1D644}_{2}(\boldsymbol{u})+\unicode[STIX]{x1D644}_{3}(\boldsymbol{u})},\end{eqnarray}$$

with

(2.70) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D644}_{1}(\boldsymbol{u})=-(\langle (u_{r}w_{r}+u_{i}w_{i})\text{D}U\rangle +\unicode[STIX]{x27E6}(u_{mr}w_{mr}+u_{mi}w_{mi})\text{D}U\unicode[STIX]{x27E7}), & \displaystyle\end{eqnarray}$$

(2.71) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D644}_{2}(\boldsymbol{u})=\langle |\text{D}\boldsymbol{u}|^{2}+(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})|\boldsymbol{u}|^{2}\rangle +\unicode[STIX]{x27E6}|\text{D}\boldsymbol{u}_{m}|^{2}+(\unicode[STIX]{x1D6FC}_{m}^{2}+\unicode[STIX]{x1D6FD}_{m}^{2})|\boldsymbol{u}_{m}|^{2}\unicode[STIX]{x27E7}, & \displaystyle\end{eqnarray}$$

(2.72) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D644}_{3}(\boldsymbol{u}) & = & \displaystyle Bn\left[\left\langle \frac{3|\text{D}w|^{2}+(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})|u|^{2}}{|\text{D}U|}\right\rangle \right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left\langle \frac{(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})|\unicode[STIX]{x1D6FC}u-\text{i}\text{D}w|^{2}+|\unicode[STIX]{x1D6FC}\text{D}u-\text{i}(\text{D}^{2}w+\unicode[STIX]{x1D6FD}w)|^{2}}{\unicode[STIX]{x1D6FD}^{2}|\text{D}U|}\right\rangle \right].\end{eqnarray}$$

In the above expressions, the following terminology is applicable:

(2.73) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle |u|^{2}=u_{r}^{2}+u_{i}^{2},\quad |\boldsymbol{u}|^{2}=|u|^{2}+\left|\frac{-\unicode[STIX]{x1D6FC}u+\text{i}\text{D}w}{\unicode[STIX]{x1D6FD}}\right|^{2}+|w|^{2},\\ \displaystyle \langle \cdot \rangle =\int _{0}^{z_{0}}(\cdot )\,\text{d}z,\quad \unicode[STIX]{x27E6}\cdot \unicode[STIX]{x27E7}=\int _{-1}^{0}(\cdot )\,\text{d}z_{m},\end{array}\right\}\end{eqnarray}$$

and $u_{r}$ and $u_{i}$ are the real and imaginary components of the perturbation. The corresponding variational problem needs to be solved. The Euler equations for the same may be given as

(2.74) $$\begin{eqnarray}\displaystyle & & \displaystyle (\text{D}^{2}-\unicode[STIX]{x1D6FC}^{2}-\unicode[STIX]{x1D6FD}^{2})^{2}w\nonumber\\ \displaystyle & & \displaystyle \qquad +\,Bn\left[-4(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})\text{D}\left(\frac{\text{D}w}{|\text{D}U|}\right)-\text{i}(\text{D}^{2}+\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})\left\{\frac{\text{D}(-\unicode[STIX]{x1D6FC}u+\text{i}\text{D}w)+\text{i}\unicode[STIX]{x1D6FD}^{2}w}{|\text{D}U|}\right\}\right]\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{\unicode[STIX]{x1D706}}{2}[(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})u\text{D}U+\text{i}\unicode[STIX]{x1D6FC}\text{D}(w\text{D}U)],\end{eqnarray}$$

(2.75) $$\begin{eqnarray}\displaystyle & & \displaystyle (\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})(\text{D}^{2}-\unicode[STIX]{x1D6FC}^{2}-\unicode[STIX]{x1D6FD}^{2})u-\text{i}\unicode[STIX]{x1D6FC}\text{D}(\text{D}^{2}-\unicode[STIX]{x1D6FC}^{2}-\unicode[STIX]{x1D6FD}^{2})w-Bn\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,\left[\frac{(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D6FD}^{2})^{2}u-\text{i}\unicode[STIX]{x1D6FC}^{3}\text{D}w}{|\text{D}U|}+\unicode[STIX]{x1D6FC}\text{D}\left(\frac{\text{D}(\text{i}\text{D}w-\unicode[STIX]{x1D6FC}u)}{|\text{D}U|}\right)\right]\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FD}^{2}}{2}w\text{D}U,\end{eqnarray}$$

(2.76) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{\unicode[STIX]{x1D6FF}^{2}}\left[\frac{\text{D}_{p}^{2}u_{m}}{\unicode[STIX]{x1D702}_{x}}-\left(\frac{\unicode[STIX]{x1D709}_{1}\unicode[STIX]{x1D6FC}_{m}^{2}}{\unicode[STIX]{x1D702}_{y}}+\frac{\unicode[STIX]{x1D709}_{2}\unicode[STIX]{x1D6FD}_{m}^{2}}{\unicode[STIX]{x1D702}_{z}}\right)(u_{m}+\text{i}\text{D}_{p}w_{m})-\frac{\text{D}_{p}\unicode[STIX]{x1D702}_{x}\text{D}_{p}(u_{m}+\text{i}\text{D}_{p}w_{m})}{\unicode[STIX]{x1D702}_{x}^{2}}\right]\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{\text{i}\unicode[STIX]{x1D706}}{\hat{d}\unicode[STIX]{x1D712}}\left(\frac{l_{4}\unicode[STIX]{x1D702}_{x}(0)}{\displaystyle l_{2}+\frac{\bar{z}_{01}^{2}}{2}}\right)(\text{D}_{p}^{2}-\unicode[STIX]{x1D6FC}_{m}^{2}-\unicode[STIX]{x1D6FD}_{m}^{2})(u_{m}+\text{i}\text{D}_{p}w_{m}).\end{eqnarray}$$

The above eigenvalue problem in $\unicode[STIX]{x1D706}$ is solved subject to the following boundary conditions:

(2.77) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u(z_{0})=0,\quad w(z_{0})=0,\quad \text{D}w(z_{0})=0,\\ \displaystyle \text{D}u(0)=\frac{\unicode[STIX]{x1D6FC}_{BJ}\hat{d}}{\unicode[STIX]{x1D6FF}\sqrt{\unicode[STIX]{x1D702}_{x}(0)}}\left[u(0)-\frac{\hat{d}Re_{m}}{Re}u_{m}(0)\right],\quad Re\,w(0)=\hat{d}Re_{m}w_{m}(0),\quad \text{D}w_{m}(0)=0,\\ u_{m}(-1)=0,\quad w_{m}(-1)=0,\quad \text{D}w_{m}(-1)=0.\end{array}\right\}\end{eqnarray}$$

As already emphasized earlier, the solution gives the conditions for no energy growth. The solution methodology for this eigenvalue problem set is also Chebyshev spectral collocation (presented in § 2.6). The results are discussed subsequently in § 3.4.

3.1 Validation

For the purpose of validation, we compare our findings with two different studies (both modal and non-modal) reported in the literature, in the limiting cases. First, we compare our non-modal results with that obtained for Poiseuille flow of a Newtonian fluid. As discussed in the preceding section, the pseudospectrum for Poiseuille flow of a Newtonian fluid $(Bn=0)$ is constructed in order to validate our findings with the literature. The same is depicted in figure 2(a,b). Both streamwise $(\unicode[STIX]{x1D6FC}\neq 0,~\unicode[STIX]{x1D6FD}=0)$ as well as spanwise $(\unicode[STIX]{x1D6FC}=0,~\unicode[STIX]{x1D6FD}\neq 0)$ perturbations are considered at $Re=3000$ . It is found that the response to streamwise perturbation is a spectrum consisting of the A, P and S branches. On the other hand, the spectrum for the spanwise perturbation is a single vertical branch. Our results coincide with those of Reddy & Henningson (Reference Reddy and Henningson1993), thereby validating the current numerical approach. In addition, we have also validated our modal stability results with Poiseuille flow of a Newtonian fluid overlying a porous layer, as reported by Chang et al. (Reference Chang, Chen and Straughan2006). Figure 2(c) shows the neutral stability curve in the case of streamwise perturbation for the set ( $Bn=0$ , $A_{inh}=0$ , $K_{y}=1$ , $K_{z}=1$ ). The neutral stability results reported in the literature by Chang et al. (Reference Chang, Chen and Straughan2006) are also displayed on the same plot (figure 2a in Chang et al. Reference Chang, Chen and Straughan2006). As evident from the figure, they show a satisfactory match.

Figure 2. Validation with literature. (a,b) Comparison with Reddy & Henningson (Reference Reddy and Henningson1993) for (a) streamwise ( $\unicode[STIX]{x1D6FC}=1,~\unicode[STIX]{x1D6FD}=0$ ) and (b) spanwise $(\unicode[STIX]{x1D6FC}=0,~\unicode[STIX]{x1D6FD}=2)$ perturbations. (c) Comparison with Chang et al. (Reference Chang, Chen and Straughan2006) for $Bn=0$ . For (a), the boundaries of the $\unicode[STIX]{x1D700}$ -pseudospectra from outer to inner are $\unicode[STIX]{x1D700}=10^{-1},~10^{-2},~10^{-3},~10^{-4}$ . For (b), they are $\unicode[STIX]{x1D700}=10^{-1},~10^{-2}$ . In (a,b), the legend denotes the following: RH-EV, eigenvalues reported by Reddy & Henningson (Reference Reddy and Henningson1993); CS-EV, eigenvalues obtained in the current study; RH-PS, pseudospectra reported by Reddy & Henningson (Reference Reddy and Henningson1993); CS-PS, pseudospectra obtained in the current study.

3.2 Variation of $\bar{z}_{01}$ and $\bar{z}_{02}$

As already discussed in the beginning, the existence of a plug zone is a characteristic feature of yield stress fluids. The thickness of the plug zone can be determined from the location of yield surfaces, $\bar{z}_{01}$ and $\bar{z}_{02}$ . Thus, studying the variations of $\bar{z}_{01}$ and $\bar{z}_{02}$ with various parameters can provide useful insight into flow characteristics. Figure 3 demonstrates the effect of Bingham number, as well as various porous layer parameters, on $\bar{z}_{01}$ and $\bar{z}_{02}$ . Figure 3(a) shows that, although the variations of both $\bar{z}_{01}$ and $\bar{z}_{02}$ are monotonic with $Bn$ , the trends are opposite to each other. The lower yield surface $\bar{z}_{01}$ reduces with $Bn$ , unlike $\bar{z}_{02}$ , which increases monotonically. Physically, this represents widening of the plug zone with $Bn$ . The values of $\bar{z}_{01}$ and $\bar{z}_{02}$ coincide at $Bn=0$ , indicating no plug zone (Newtonian fluid). Essentially, as yield stress increases, there is a gradual transition from single shear (Newtonian) to double shear flow with a plug zone in the middle. Figure 3(b,c) reveals that both $\bar{z}_{01}$ and $\bar{z}_{02}$ increase with depth ratio $\hat{d}$ and slip parameter $\unicode[STIX]{x1D6FC}_{BJ}$ . However, there is saturation at higher values of $\hat{d}$ and $\unicode[STIX]{x1D6FC}_{BJ}$ . Increase in $\unicode[STIX]{x1D6FC}_{BJ}$ translates to rise in velocity gradient at the fluid–porous interface. This results in higher momentum transfer in the transverse direction, leading to plug behaviour in the upper shear zone and a more pronounced lower shear zone, represented by an upward shift of both the yield surfaces (Sengupta & De Reference Sengupta and De2019). Of course, there cannot be prolonged monotonic rise of the interface velocity gradient, leading to saturation at higher values of $\unicode[STIX]{x1D6FC}_{BJ}$ . As per the Beavers–Joseph condition, depth ratio is also proportional to the interface velocity gradient, thereby exhibiting similar trend as that of $\unicode[STIX]{x1D6FC}_{BJ}$ . As shown in figure 3(d), the trend is reverse for $\unicode[STIX]{x1D6FF}$ . Darcy number represents the permeability of the porous layer. As the porous layer becomes more permeable, the porous layer velocity increases. Thus, there is a reduction in the velocity gradient at the fluid–porous interface, resulting in a downward shift of the yield surfaces with $\unicode[STIX]{x1D6FF}$ . In figure 3(e), the yield surfaces exhibit an asymmetric variation with inhomogeneity factor $A_{inh}$ with respect to $A_{inh}=0$ , in addition to being non-monotonic. For $A_{inh}>0$ , there is a reduction in the values of $\bar{z}_{01}$ and $\bar{z}_{02}$ . Increase in inhomogeneity amounts to increase in overall permeability of the porous medium. Owing to this, the trend of a downward shift in the yield surfaces with $A_{inh}$ somewhat resembles the behaviour of $\unicode[STIX]{x1D6FF}$ . High negative values of $A_{inh}$ indicate drastic (exponential) reduction in the permeability, due to which the fluid–porous interface almost behaves as an impermeable wall. Thus, the porous layer no longer affects the location of the yield surfaces, due to which both $\bar{z}_{01}$ and $\bar{z}_{02}$ become almost invariant at large negative  $A_{inh}$ .

Figure 3. Variation of $\bar{z}_{01}$ and $\bar{z}_{02}$ with various parameters: (a) Bingham number; (b) depth ratio; (c) slip coefficient; (d) Darcy number; and (e) inhomogeneity factor.

3.3 Eigenvalue spectra

The dependence of the normal modes on $Bn$ is shown in figures 4–6, in addition to the effect of the various parameters of the porous layer. It is observed that increase in $Bn$ from 0.01 to 0.1 leads to reduction in $c_{i}$ . However, there is no value of $c_{i}>0$ in each of the cases. Unlike non-porous channel shear flow, there is no vertical S branch here. Instead, there is an oblique branch at approximately $c_{r}=0.48$ (the oblique line is represented by dotted line in figure 4 a). The complex interplay of the fluid and porous modes makes the branches oblique. Also, two distinct branches can be located on either side of this oblique line, similar to A and P branches of non-porous Newtonian flow. Reduction in the value of $c_{i}$ is due to the increase in the viscoplasticity of the fluid with $Bn$ . Thus, the flow becomes increasingly stable with increased $Bn$ (even though it is never linearly unstable).

Figures 4–6 explore the effects of the porous layer parameters on the spectra. The P and S branches are not very sharp, unlike the case of non-porous flow configuration. Although these A and P branches are diffused, there is an interesting observation regarding the imaginary part of the eigenvalues. To explain the same, let us observe the two modes marked 1 and 2 in figure 4(a). Their coordinates in figure 4(a) are as follows: mode 1 $(0.53624,-0.06912)$ and mode 2 $(0.72222,-0.06976)$ . Thus, for $\hat{d}=0.1$ , mode 1 would dominate the transition behaviour. But as $\hat{d}$ is increased to 1 in figure 4(b), the coordinates for the two modes are: mode 1 $(0.53624,-0.11480)$ and mode 2 $(0.72222,-0.11411)$ . So, mode 2 becomes dominant in this case. Thus, with increase in depth ratio from 0.1 to 1, there is a ‘switching’ of the modes carrying the maximum value of $c_{i}$ , from mode 1 to mode 2. Similar behaviour of mode switching is also observed in the case of slip parameter and permeability, depicted in figures 5 and 6. Thus, slip parameter and permeability also induce mode switching in the eigenvalue spectra. However, the mode switching is not explicitly demonstrated in figures 5 and 6, to avoid repetition. The mode switching phenomenon leads to the development of unique behaviour in the neutral curves, as discussed later in § 3.4.

Figure 4. Eigenvalue spectra studying the effect of depth ratio $\hat{d}$ : (a)  $\hat{d}=0.1$ and (b)  $\hat{d}=1$ . Dotted lines represent the oblique branches.

Figure 5. Eigenvalue spectra studying the effect of slip parameter $\unicode[STIX]{x1D6FC}_{BJ}$ : (a)  $\unicode[STIX]{x1D6FC}_{BJ}=0.1$ and (b)  $\unicode[STIX]{x1D6FC}_{BJ}=0.4$ . Dotted lines represent the oblique branches.

Figure 6. Eigenvalue spectra studying the effect of Darcy number $\unicode[STIX]{x1D6FF}$ : (a)  $\unicode[STIX]{x1D6FF}=0.0001$ and (b)  $\unicode[STIX]{x1D6FF}=0.001$ . Dotted lines represent the oblique branches.

3.3.1 Effect of the porous modes

In order to have a better understanding of the effect of viscoplasticity on porous modes, the effect of $Bn$ on the eigenfunction is studied. As demonstrated in figure 7, for low $Bn$ ( $Bn=0.01$ ), there is a flow reversal near the interface. Also, there is significant protrusion of the momentum into the porous layer. With gradual increase in $Bn$ to $Bn=0.05$ , the degree of protrusion reduces. The extent of flow reversal is also less. But both these effects of flow reversal and protrusion are exhibited at this stage. However, with further increase in $Bn$ to 0.1, there is no more protrusion. Only flow reversal near the interface is observed. For higher $Bn$ ( $Bn=0.2$ ), even the flow reversal ceases to exist. All these effects reveal the effect of the porous layer in deciding the mode of stability. At lower $Bn$ , Newtonian flow characteristics are more dominant. Thus, there is significant perturbation in the porous layer as well, in addition to the fluid layer. This is the cause of the flow reversal and protrusion in figure 7(a,b). At this stage, the effects of the fluid layer and the porous layer on flow transition are comparable. As $Bn$ further increases to 0.1, the protrusion of the perturbation becomes negligible, owing to the onset of yield stress effects. At even higher $Bn$ , the viscoplastic effects dominate the porous layer completely. As a result, the porous layer is unable to influence the transition behaviour. The dominant role is now played by the fluid layer alone.

Figure 7. Variation of the eigenfunctions with $Bn$ : (a)  $Bn=0.01$ ; (b)  $Bn=0.05$ ; (c)  $Bn=0.1$ ; and (d)  $Bn=0.2$ . Here $u_{m}$ is the eigenfunction corresponding to the domain $-1\leqslant z_{m}\leqslant 0,$ whereas $u$ is the eigenfunction corresponding to the domain $0\leqslant z\leqslant 1.$ Solid lines refer to the real components, while the dashed lines refer to the imaginary components.

Figure 8. Plots of energy Reynolds number versus spanwise wavenumber ( $\unicode[STIX]{x1D6FC}=0$ ) to depict the effects of: (a) depth ratio $\hat{d}$ ; (b) slip parameter $\unicode[STIX]{x1D6FC}_{BJ}$ ; (c) Darcy number $\unicode[STIX]{x1D6FF}$ ; (d) inhomogeneity factor $A_{inh}$ ; (e) spanwise anisotropy parameter $\unicode[STIX]{x1D709}_{1}$ ; and (f) normal anisotropy parameter $\unicode[STIX]{x1D709}_{2}$ . The numbers in the panels denote the respective values of the parameters. Here $Bn=0.2$ in all cases.

3.4 Non-modal analysis

3.4.2 Growth rate

Figure 9 shows the effects of the porous layer parameters on the growth rate curves for a spanwise perturbation. It is observed that, with increase in depth ratio in figure 9(a), the maximum possible growth as well as the time corresponding to maximum growth is reduced. However, beyond $\hat{d}=$ 1, the reduction is not substantial. The slip parameter also depicts a similar tendency in figure 9(b), but the reduction in $G(t)$ is significant even at high slip coefficient, i.e. at high velocity gradient at the fluid–porous interface. As far as permeability is concerned in figure 9(c), increase in Darcy number leads to monotonic increase in the maximum growth, as well as the time of maximum growth. These effects illustrate the relative roles played by each of the porous layer parameters in deciding the non-modal growth. With increase in depth ratio, the fluid layer gradually grows in thickness. Thus, both the porous and the fluid layers simultaneously contribute to intermediate growth. This explains the growth characteristics with $\hat{d}$ . For the slip parameter, increase signifies increased momentum transfer across the interface, thus establishing the dominance of the fluid layer in influencing the criticality. Increase in permeability (characterized by the Darcy number) assists relatively greater flow in the porous layer. Thus, with increase in permeability, the porous layer is able to influence the transient amplifications. Thus, the complex interactions among these parameters decide the short-time growth.

Figure 9. Growth rate curves depicting the effects of: (a) depth ratio $\hat{d}$ ; (b) slip parameter $\unicode[STIX]{x1D6FC}_{BJ}$ ; and (c) Darcy number $\unicode[STIX]{x1D6FF}$ . The numbers in the panels denote the respective values of the parameters. Here $Re=600$ and $\unicode[STIX]{x1D6FD}=4$ in all cases.

In further discussions, we consider oblique perturbations ( $\unicode[STIX]{x1D6FC}\neq 0,~\unicode[STIX]{x1D6FD}\neq 0$ ), and the effect of the parameters of the porous layer on the same. Figure 10 describes the effect of viscoplasticity, in terms of contours of maximum growth in the plane of streamwise and spanwise perturbations. At $Bn=0.01$ , the maximum growth is achieved by a spanwise perturbation. As the perturbations become oblique, the maximum growth is reduced. It is observed from figure 10(a) that the disturbance is streamwise at lower values of optimal, whereas it turns into a spanwise disturbance at the highest value of $G_{max}$ . This is in conformity with earlier studies highlighting the role of spanwise disturbances in amplifying growth (Liu & Liu Reference Liu and Liu2014; Liu et al. Reference Liu, Ding and Hu2018). Nouar et al. (Reference Nouar, Kabouya, Dusek and Mamou2007) also constructed contours of $G_{max}$ for plane Bingham–Poiseuille flow. Their results are also presented in figure 10(a), for the sake of comparison. It is observed from figure 10(a) that, for a particular combination of $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ , the $G_{max}$ values obtained in the current study are higher than those of Nouar et al. (Reference Nouar, Kabouya, Dusek and Mamou2007). This is because of the dominant role played by the porous layer in the present study for low $Bn$ , as demonstrated in figure 7. Figure 10(b) depicts the contours of $G_{max}$ for $Bn=0.1$ . A decrease in $G_{max}$ with $Bn$ is observed, for the same combination of $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ . The decrease in maximum growth with increased viscoplasticity is due to the enhanced viscous dissipation. From figure 10(a,b), it is found that the magnitude of optimal growth reduces with $Bn$ . Thus, although there is asymptotic long-time decay, viscoplasticity plays a stabilizing role in amplifying the short-time growth.

Figure 10. Effect of Bingham number: contours of $G_{max}$ at (a)  $Bn=0.02$ and (b)  $Bn=0.1$ . In (a), dashed lines refer to the studies of Nouar et al. (Reference Nouar, Kabouya, Dusek and Mamou2007).

Figure 11. Effect of Darcy number: contours of $G_{max}$ at (a) $\unicode[STIX]{x1D6FF}=5\times 10^{-4}$ and (b)  $\unicode[STIX]{x1D6FF}=5\times 10^{-3}$ .

Figures 11–13 depict the effects of Darcy number, depth ratio and slip parameter with regards to the oblique perturbations. In figure 11, with increase in permeability, the optimal perturbation is increased. It increases from 150 to 290 as the permeability increases from $\unicode[STIX]{x1D6FF}=5\times 10^{-4}$ to $5\times 10^{-3}$ . This again demonstrates the destabilizing role played by permeability in creating short-time amplifications. The reason for this is the increased momentum in the porous layer with permeability. Also, as the perturbations turn from spanwise to oblique to streamwise, the maximum growth is lowered. However, at higher permeability, the perturbation producing maximum growth is oblique, although the streamwise component is negligible compared to the spanwise component.

Figure 12 expresses the effect of depth ratio. It is observed that the maximum growth is reduced with the increase in fluid layer thickness. Thus, as the relative influence of the fluid layer compared to the porous layer is increased in the flow system, there is a reduction in transient amplifications. The fluid modes restrict the growth of amplifications, leading to the reduction in the magnitudes of $G_{max}$ . At very low depth ratio, the fluid layer has a small thickness, and a large part of the momentum is influenced by the porous layer. This explains the trend of the maximum growth with depth ratio. A similar behaviour is depicted by slip parameter in figure 13. At lower slip coefficient, the maximum growth perturbation is spanwise in nature. However, as the slip coefficient is increased, the perturbation becomes oblique.

Figure 12. Effect of depth ratio: contours of $G_{max}$ at (a)  $\hat{d}=0.3$ and (b)  $\hat{d}=3$ .

Figure 13. Effect of slip parameter: contours of $G_{max}$ at (a) $\unicode[STIX]{x1D6FC}_{BJ}=0.05$ and (b)  $\unicode[STIX]{x1D6FC}_{BJ}=0.2$ .

3.5 Mechanism of transient growth

The shape of the optimal perturbation is explored in detail in figure 14. Figure 14(a) shows the shape of the optimal perturbation for $Bn=0.02$ and $\hat{d}=0.3$ . It is observed that the perturbation is spanwise. Moreover, the optimal perturbations also protrude into the porous layer. This explains the significance of the porous layer in deciding the flow transition characteristics, as also revealed in figure 7. Figure 14(b,c) displays the optimal perturbations for a high value of depth ratio ( $\hat{d}=3$ ) for two different Bingham numbers, $Bn=0.02$ and 0.3. From the plots, it is observed that, at lower $Bn$ , the optimal perturbation consists of spanwise streaks. However, it becomes oblique with larger $Bn$ . This obliqueness of the optimal perturbation is possibly because of the fact that the effective viscosity of the Bingham fluid varies in a nonlinear fashion across the width of the channel. Moreover, at higher depth ratio, there is no protrusion of the optimal perturbation into the porous layer. The perturbations are confined to the fluid layer alone. This explains the relative significance of the fluid and the porous layers at high depth ratio.

It is worth while to discuss the mechanism of transient growth occurring here. In the case of inviscid shear flow, the streamwise component of the velocity perturbation has been found to grow linearly with time (Ellingsen & Palm Reference Ellingsen and Palm1975). This mechanism is commonly referred to as ‘lift-up’ (Reddy & Henningson Reference Reddy and Henningson1993; Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993). However, for viscous shear flows, viscous dissipation plays an important role. In this case, there is initial linear growth of the streamwise velocity with time due to lift-up, and consequent damping due to viscosity (Butler & Farrell Reference Butler and Farrell1992; Farrell & Ioannou Reference Farrell and Ioannou1993; Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993). The interplay of lift-up and Orr mechanisms contributes to transient growth of the oblique disturbances (Farrell & Ioannou Reference Farrell and Ioannou1993). They both act simultaneously, and the interaction between the two leads to the development of the transient growth. In our analysis, in the case of lower $Bn$ , the shape of the optimal disturbance is found to be similar to that reported by Reddy & Henningson (Reference Reddy and Henningson1993) in the case of plane Newtonian–Poiseuille flow. It is known that the lift-up mechanism dominates transient growth in the case of such perturbations. However, at higher $Bn$ , the perturbation becomes oblique in shape. As discussed, when the perturbation is oblique, both lift-up and Orr mechanisms become comparable. Thus, it may be said that transient growth is dominated by lift-up at lower $Bn$ , while both lift-up and Orr mechanisms contribute to transient response at higher  $Bn$ .

Figure 14. Optimal perturbation at $x=0$ , $Re=600$ for: (a) $\hat{d}=0.3$ , $Bn=0.02$ ; (b)  $\hat{d}=3$ , $Bn=0.02$ ; and (c) $\hat{d}=3$ , $Bn=0.3$ . Contour levels from outer to inner denote the values of $u$ in the following order $[0.1,0.3,0.5,0.7,0.9]$ .

3.6 Description of the fluid–porous interface and the porous layer

As may be observed from the expressions of the velocity profiles, there is a discontinuity in the axial velocity profile at the interface between the fluid and the porous layers. This discontinuity arises because of the Beavers–Joseph boundary condition, which is generally obtained from experimental study, and thus backed by experimental verification. The authors acknowledge the existence of more refined interface conditions that have been proposed in recent times, for example, the ones by Ochoa-Tapia & Whitaker (Reference Ochoa-Tapia and Whitaker1995a ,Reference Ochoa-Tapia and Whitaker b ) and Bars & Worster (Reference Bars and Worster2006). However, each of these interface conditions also carries its own set of limitations. For example, the validity of the Bars–Worster boundary condition is questionable in the case of a sharp fluid–porous interface (Bars & Worster Reference Bars and Worster2006). In addition, some of these conditions, like the boundary condition by Ochoa-Tapia & Whitaker, are generally used together with the Brinkman model. The limitations of the Brinkman model with regard to the current study have already been discussed in § 2.1. Taking all these factors into consideration, the current study uses the Beavers–Joseph condition to model the fluid–porous interface. The same has been adopted by other researchers as well (Chang et al. Reference Chang, Chen and Straughan2006; Deepu et al. Reference Deepu, Anand and Basu2015, Reference Deepu, Kallurkar, Anand and Basu2016; Chang et al. Reference Chang, Chen and Chang2017).

It may be mentioned that the Forchheimer model is not considered in the current analysis. Incorporation of the Forchheimer model ensures that the inertial effects are taken care of (Lyubimova et al. Reference Lyubimova, Lyubimov, Baydina, Kolchanova and Tsiberkin2016). The Forchheimer model is generally considered when the Reynolds number in the porous layer $(Re_{m})$ is high. However, in the present study, $Re_{m}$ in the porous layer is less than $Re$ by approximately four orders of magnitude. This can be easily checked from the relation (2.47) between $Re$ and $Re_{m}$ . By taking representative values of the parameters ( $Bn=0.2$ , $\unicode[STIX]{x1D6FD}_{0}=0.016$ , $\hat{d}=0.3$ , $\unicode[STIX]{x1D6FC}_{BJ}=0.2$ , $\unicode[STIX]{x1D6FF}=0.001$ , $A_{inh}=2$ ), $Re/Re_{m}=1.29\times 10^{4}$ . Therefore, $Re_{m}$ is significantly small even for a seemingly high $Re$ . Thus, there is negligible inertial effect in the porous layer. This is the primary reason behind neglecting inertial effects by a majority of researchers (Chang et al. Reference Chang, Chen and Straughan2006; Hill & Straughan Reference Hill and Straughan2008, Reference Hill and Straughan2009; Deepu et al. Reference Deepu, Anand and Basu2015, Reference Deepu, Kallurkar, Anand and Basu2016; Chang et al. Reference Chang, Chen and Chang2017).