1. Introduction
Many flows of engineering and scientific interest involve permeable substrates. The presence of such complex substrates can significantly alter the behaviour of the near-wall turbulence. A growing body of work suggests that appropriately designed anisotropic permeable substrates have the potential to suppress the dynamically important near-wall (NW) cycle (Robinson Reference Robinson1991; Waleffe Reference Waleffe1997; Jiménez & Pinelli Reference Jiménez and Pinelli1999) and reduce drag in wall-bounded turbulent flows (Hahn, Je & Choi Reference Hahn, Je and Choi2002; Itoh et al. Reference Itoh, Tamano, Iguchi, Yokota, Akino, Hino and Kubo2006; Abderrahaman-Elena & García-Mayoral Reference Abderrahaman-Elena and García-Mayoral2017; Gómez-de-Segura, Sharma & García-Mayoral Reference Gómez-de-Segura, Sharma and García-Mayoral2018; Rosti, Brandt & Pinelli Reference Rosti, Brandt and Pinelli2018; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). Permeable materials can also be used to enhance turbulent mixing for applications in the development of high-efficiency thermal management systems and chemical reactors (Gad-el Hak Reference Gad-el Hak2007). Previous laboratory experiments and numerical simulations have provided significant insight into the effect of both isotropic and anisotropic permeable substrates on turbulent boundary layer and channel flows (e.g. Hahn et al. Reference Hahn, Je and Choi2002; Breugem, Boersma & Uittenbogaar Reference Breugem, Boersma and Uittenbogaar2006; Manes, Poggi & Ridolfi Reference Manes, Poggi and Ridolfi2011; Zampogna & Bottaro Reference Zampogna and Bottaro2016; Efstathiou & Luhar Reference Efstathiou and Luhar2018; Rosti et al. Reference Rosti, Brandt and Pinelli2018; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019; Kim et al. Reference Kim, Blois, Best and Christensen2020). For instance, it is well known that flows over porous materials are susceptible to a Kelvin–Helmholtz (KH) instability that gives rise to spanwise-coherent energetic rollers (Jiménez & Pinelli Reference Jiménez and Pinelli1999; Breugem et al. Reference Breugem, Boersma and Uittenbogaar2006; Efstathiou & Luhar Reference Efstathiou and Luhar2018). The mechanism that could lead to drag reduction in turbulent flows over anisotropic materials is also reasonably well understood (Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). Despite these advances, there are few reduced-complexity models that can be used to predict how a given porous substrate will affect the turbulent flow, i.e. whether it will suppress the NW cycle or give rise to KH rollers. Given the vast parameter space available in the development of porous materials for passive flow control, such models can be useful tools for design and optimization. In this study, we extend the resolvent analysis framework (McKeon & Sharma Reference McKeon and Sharma2010; McKeon Reference McKeon2017) to develop reduced-complexity models for turbulent flows over porous substrates. We focus on evaluating the effect of anisotropic permeable materials that can give rise to drag reduction. However, these models can also be used to evaluate the effect of porous materials for other applications or to provide insight into environmental flows over granular beds and vegetation canopies.
1.1. Previous work
Recent numerical simulations show that streamwise-preferential permeable materials have the potential to yield as much as 25 % drag reduction in turbulent flows (Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). The physical mechanism underlying this drag reduction is similar to the mechanism that yields drag reduction over riblets (Walsh & Lindemann Reference Walsh and Lindemann1984; Luchini, Manzo & Pozzi Reference Luchini, Manzo and Pozzi1991; Luchini Reference Luchini1996; Bechert et al. Reference Bechert, Bruse, Hage, Van Der Hoeveb and Hoppe1997; García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011a). Specifically, for anisotropic materials that have larger streamwise permeability than spanwise permeability, the streamwise mean flow penetrates into the substrate to a larger extent compared with the spanwise cross-flows arising from turbulent fluctuations (Luchini et al. Reference Luchini, Manzo and Pozzi1991). In other words, there is an offset between the virtual origins perceived by the mean flow and the turbulent fluctuations. The virtual origin for the turbulent cross-flow can also be interpreted as the location for which the quasi-streamwise vortices associated with the NW cycle perceive a no-slip wall (Luchini Reference Luchini1996; García-Mayoral, Gómez de Segura & Fairhall Reference García-Mayoral, Gómez de Segura and Fairhall2019; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). Importantly, the offset in virtual origins for the mean streamwise flow and the turbulent cross-flows weakens the quasi-streamwise vortices associated with the NW cycle. This reduces turbulent momentum transfer towards the substrate, which leads to a decrease in skin friction.
Building on this concept, Abderrahaman-Elena & García-Mayoral (Reference Abderrahaman-Elena and García-Mayoral2017) used the Brinkman equations to establish a relationship between the streamwise and spanwise permeabilities (
$K^+_x,K^+_z$) and the streamwise and spanwise slip lengths (
$\ell ^+_x,\ell ^+_z$), i.e. the lengths from the interface where the virtual wall would be perceived. As shown in figure 1, 
$x$, 
$y$ and 
$z$ represent the streamwise, wall-normal and spanwise coordinates, respectively. A superscript 
$+$ denotes normalization with respect to the friction velocity, 
$u_\tau$, and viscosity, 
$\nu$. The analysis of Abderrahaman-Elena & García-Mayoral (Reference Abderrahaman-Elena and García-Mayoral2017) showed that 
$\ell ^+_x\propto \sqrt {\smash {K_x^+}\vphantom {K_x^+}}$ and 
$\ell ^+_z\propto \sqrt {\smash {K_z^+}\vphantom {K_x^+}}$ if the height of the substrate, 
$H$, is much larger than the permeability length scales 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}}$ and 
$\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$. The achievable drag reduction was estimated to be proportional to the difference between the streamwise and spanwise slip lengths, 
$\Delta D \propto \ell ^+_x-\ell ^+_z$, or equivalently, 
$\Delta D \propto \sqrt {\smash {K_x^+}\vphantom {K_x^+}}-\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$. This is consistent with the findings of Busse & Sandham (Reference Busse and Sandham2012), who studied the effect of anisotropic slip length boundary conditions in turbulent channel flow simulations.
Figure 1. Schematic showing the symmetric channel flow configuration considered in this paper.
Recent direct numerical simulation (DNS) results obtained by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) provide further support for the scaling developed by Abderrahaman-Elena & García-Mayoral (Reference Abderrahaman-Elena and García-Mayoral2017). For these simulations, the Brinkman equations were used to model flow inside the permeable substrate and the Navier–Stokes equations were used in the fluid domain. Stresses and velocities were matched at the interface between the permeable substrate and the unobstructed flow. Simulations carried out using this Brinkman model led to results that were very similar to the simulation results obtained by Breugem et al. (Reference Breugem, Boersma and Uittenbogaar2006) using the full volume-averaged Navier–Stokes (VANS) equations for an isotropic permeable medium with permeability 
$\sqrt {K^+} \approx 1$. Moreover, Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) also showed that the simplified Brinkman equations captured the dominant dynamics emerging from the full VANS equations. The anisotropic permeable substrates considered by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) were characterized by a permeability tensor of the form 
${\boldsymbol{\mathsf{K}}} =\text {diag}({K}_{x},{K}_{y},{K}_{z})$, and the wall-normal permeability was set to be equal to the spanwise permeability, 
$K_y = K_z$. The effect of substrate anisotropy was evaluated by systematically varying the streamwise and spanwise permeability, such that the anisotropy ratio 
$\phi _{xy}={\sqrt {\smash {K_x^+}\vphantom {K_x^+}}}/{\sqrt {\smash {K_y^+}\vphantom {K_x^+}}}$ ranged between 
$3.6$ and 
$11.4$. As predicted by Abderrahaman-Elena & García-Mayoral (Reference Abderrahaman-Elena and García-Mayoral2017), these simulations show that for surfaces for which the permeability length scale is smaller than the size of the near wall turbulent structures, the initial decrease in drag is proportional to the difference in the slip length scales, 
$\Delta D \propto \sqrt {\smash {K_x^+}\vphantom {K_x^+}}-\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$. However, the simulations also show that the achievable drag reduction is limited by the appearance of energetic spanwise rollers resembling KH vortices. Such rollers have been observed over isotropic permeable substrates (Breugem et al. Reference Breugem, Boersma and Uittenbogaar2006; Efstathiou & Luhar Reference Efstathiou and Luhar2018), and they also contribute to performance degradation for riblets (García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011b, Reference García-Mayoral and Jiménez2012). Linear stability analyses and simulation results show that the appearance of these spanwise rollers is controlled primarily by the wall-normal permeability (Abderrahaman-Elena & García-Mayoral Reference Abderrahaman-Elena and García-Mayoral2017; Gómez-de-Segura et al. Reference Gómez-de-Segura, Sharma and García-Mayoral2018; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). Specifically, simulation results show that spanwise rollers become increasingly energetic as the wall-normal permeability exceeds 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \approx 0.4$. The additional Reynolds shear stress produced by these rollers causes performance to deteriorate and ultimately leads to an increase in drag.
These prior efforts show that the drag-reducing performance of anisotropic permeable substrates is dictated by two key factors: the suppression of NW cycle and the emergence of energetic spanwise-coherent rollers. The physically motivated slip length arguments presented earlier provide useful insight into the first effect. These arguments predict that the initial decrease in drag is proportional to the difference between the streamwise and spanwise permeability length scales, 
$\Delta D \propto \sqrt {\smash {K_x^+}\vphantom {K_x^+}}-\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$. However, it is unclear if this relationship also holds for more complex surfaces that are not characterized by diagonal permeability tensors of the form 
${\boldsymbol{\mathsf{K}}} =\text {diag}({K}_{x},{K}_{y},{K}_{z})$. Moreover, these slip length arguments are based on solutions to the Brinkman equations in the porous medium, which are coupled to the fluid domain via velocity and stress-matching boundary conditions at the fluid–substrate interface. Recent studies show that the interfacial boundary conditions may be better characterized by a slip-length tensor (Lācis & Bagheri Reference LĀcis and Bagheri2017; Bottaro Reference Bottaro2019; Lācis et al. Reference LĀcis, Sudhakar, Pasche and Bagheri2020). The effect of such slip-length models on the near-wall turbulence remains to be studied. Similarly, linear stability analyses are able to predict the emergence of spanwise-coherent KH rollers over permeable substrates as the wall-normal permeability increases. However, such analyses fail to accurately predict the exact threshold for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$ beyond which KH rollers become energetic (Abderrahaman-Elena & García-Mayoral Reference Abderrahaman-Elena and García-Mayoral2017; Gómez-de-Segura et al. Reference Gómez-de-Segura, Sharma and García-Mayoral2018). Furthermore, the streamwise wavelengths predicted to be most unstable do not match the length scale of the spanwise rollers observed in simulations.
1.2. Contribution and outline
In this study, we develop a reduced-order modelling framework grounded in resolvent analysis (McKeon & Sharma Reference McKeon and Sharma2010; McKeon Reference McKeon2017) that can be used to predict the effect of substrates with known permeability on the NW cycle and test for the emergence of KH rollers. Under the resolvent formulation, the Navier–Stokes equations are interpreted as a forcing–response system in which the nonlinear convective terms are treated as the forcing to the linear system that generates a velocity and pressure response. A gain-based decomposition of the resolvent operator, which is the linear transfer function that maps the nonlinear forcing to the velocity and pressure response, is used to identify high-gain forcing and response modes across spectral space. Specific high-gain response modes (resolvent modes) have been shown to serve as useful models for dynamically important flow features such as the NW cycle (Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013). This means that, as a starting point, the effect of any control can be evaluated on these individual resolvent modes instead of the full turbulent flow field. Indeed, previous studies show that the gain for the NW resolvent mode is a useful predictor of drag reduction performance for both active (Luhar, Sharma & McKeon Reference Luhar, Sharma and McKeon2014; Nakashima, Fukagata & Luhar Reference Nakashima, Fukagata and Luhar2017; Toedtli, Luhar & McKeon Reference Toedtli, Luhar and McKeon2019) and passive (Luhar, Sharma & McKeon Reference Luhar, Sharma and McKeon2015; Chavarin & Luhar Reference Chavarin and Luhar2020) control techniques in wall-bounded turbulent flows. In particular, recent work by Chavarin & Luhar (Reference Chavarin and Luhar2020) shows that the gain for the NW resolvent mode is a useful surrogate for total drag reduction in turbulent flows over riblets. Riblet geometries that lead to drag reduction in experiments and high-fidelity simulations are found to reduce the forcing–response gain for the NW resolvent mode relative to its smooth wall value. In addition, Chavarin & Luhar (Reference Chavarin and Luhar2020) show that the resolvent framework is also able to predict the emergence of spanwise rollers resembling KH rollers over certain riblet geometries, which is consistent with previous DNS results (García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011b). Motivated by these prior modelling successes, we consider the effect of anisotropic porous substrates on resolvent modes resembling the NW cycle and spanwise-coherent structures resembling KH rollers. To enable a direct comparison with the simulation results of Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019), we consider a symmetric channel geometry at friction Reynolds number 
$\textit {Re}_{\tau } = 180$ and substrates with identical wall-normal and spanwise permeabilities, i.e. substrates with 
$\phi _{yz}=\sqrt {\smash {K_y^+}\vphantom {K_x^+}}/\sqrt {\smash {K_z^+}\vphantom {K_x^+}}=1$. We model the flow in the substrate using the VANS equations, in which the effect of the permeable substrate is included via a permeability tensor. However, this modelling framework can be extended to include more sophisticated interfacial boundary conditions (Lācis & Bagheri Reference LĀcis and Bagheri2017; Lācis et al. Reference LĀcis, Sudhakar, Pasche and Bagheri2020), and to account for inertial effects via the so-called Forchheimer term (Whitaker Reference Whitaker1999; Breugem et al. Reference Breugem, Boersma and Uittenbogaar2006).
The remainder of this paper is structured as follows. Section 2 describes the permeable substrate model used here, the extended resolvent formulation, as well as the numerical methods used for resolvent analysis. Section 3 presents model predictions for the effect of anisotropic permeable substrates on the NW resolvent mode as well as spanwise-coherent resolvent modes resembling KH rollers. These predictions are compared with DNS results from Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019). We also pursue a limited sensitivity analysis of model predictions to the exact form of the mean profile used to construct the resolvent operator. Specifically, we compare predictions made using a synthetic mean profile that is computed using an eddy viscosity model with predictions generated using the mean velocity profiles obtained in DNS by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019). Section 4 concludes the paper.
2. Methods
In this section, we present the equations used to model flow in the porous medium (§ 2.1), briefly describe the resolvent formulation and discuss its extension to account for permeable substrates (§ 2.2), present the boundary conditions imposed at the fluid–substrate interface (§ 2.3), discuss the mean velocity profiles used to construct the resolvent operator (§ 2.4) and describe the numerical method used to implement the analysis (§ 2.5).
2.1. Accounting for permeable substrates
The resolvent framework is reformulated using the VANS equations. Volume-averaging gives rise to two additional terms: a term representing the subfilter scale stresses and a surface filter term that accounts for the force exerted by the solid phase of the permeable medium on the fluid phase (Whitaker Reference Whitaker1969, Reference Whitaker1996, Reference Whitaker1999). A typical closure model for the surface filter term involves parameterizing the flow resistance using the Darcy permeability tensor and the so-called Forchheimer correction term that accounts for inertial effects (Whitaker Reference Whitaker1996). This model has been used in previous numerical simulations of flow over porous substrates (Breugem et al. Reference Breugem, Boersma and Uittenbogaar2006; Rosti, Cortelezzi & Quadrio Reference Rosti, Cortelezzi and Quadrio2015; Rosti et al. Reference Rosti, Brandt and Pinelli2018) as well as in linear stability analyses (Tilton & Cortelezzi Reference Tilton and Cortelezzi2006, Reference Tilton and Cortelezzi2008). The VANS equations and continuity constraint for flow through a porous medium with porosity 
$\varepsilon$, dimensionless permeability 
$\boldsymbol{\mathsf{K}}$ and dimensionless Forchheimer resistance 
$\boldsymbol{\mathsf{F}}$ can be expressed as
\begin{align} \frac{\partial \langle\boldsymbol{u}\rangle}{\partial t} + \frac{1}{\varepsilon}\boldsymbol{\nabla}\boldsymbol{\cdot} \left(\varepsilon\langle\boldsymbol{u}\rangle\langle\boldsymbol{u}\rangle+\varepsilon\boldsymbol{\tau} \right) ={-} \frac{1}{\varepsilon}\boldsymbol{\nabla}(\varepsilon \langle p\rangle) + \frac{1}{\varepsilon\textit{Re}_{\tau}}\boldsymbol{\nabla}^2 (\varepsilon\langle\boldsymbol{u}\rangle)-\frac{\varepsilon}{\textit{Re}_{\tau}}\boldsymbol{\mathsf{K}}^{{-}1}(\boldsymbol{\mathsf{I}}+\boldsymbol{\mathsf{F}})\langle\boldsymbol{u}\rangle \end{align}and
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}(\varepsilon\langle\boldsymbol{u}\rangle) = 0.\end{equation}
In the equations above, 
$\langle \cdot \rangle$ denotes an intrinsic volume average, 
$\langle \boldsymbol {u}\rangle$ is the dimensionless velocity, 
$\langle p \rangle$ is the dimensionless pressure and 
$t$ is dimensionless time. The equations presented above have been normalized using the channel half-height 
$h$ and the friction velocity 
$u_\tau$. The friction Reynolds number is given by 
$\textit {Re}_{\tau }=u_{\tau }h/\nu$ and the dimensionless permeability defined as 
$\boldsymbol{\mathsf{K}} = \boldsymbol{\mathsf{K}}^{{\dagger}ger}ger /h^2$, where 
$\boldsymbol{\mathsf{K}}^{{\dagger}ger}ger$ is the dimensional permeability. As discussed in Breugem et al. (Reference Breugem, Boersma and Uittenbogaar2006), the volume averaging operation serves to filter the flow field and only passes on information on the large-scale structure. The quantity 
$\boldsymbol {\tau }=\langle \boldsymbol {u}\boldsymbol {u}\rangle -\langle \boldsymbol {u}\rangle \langle \boldsymbol {u}\rangle$ represents the subfilter scale stresses which arise from volume averaging the Navier–Stokes equations. Note that the averaging volume must be chosen to ensure that the resulting flow field is continuous, i.e. defined in both the solid and the fluid phase. Moreover, the resulting flow field must contain negligible variations at scales smaller than the dimensions of the averaging volume. Thus, the size of the averaging volume must be larger than the characteristic length scales associated with the pore-scale geometry of the material being considered.
For our analysis we consider the following simplifications. Consistent with prior numerical simulations (Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019), we focus on substrates characterized by a permeability tensor of the form 
${\boldsymbol{\mathsf{K}}} =\text {diag}({K}_{x},{K}_{y},{K}_{z})$ with the ratio of the wall-normal and spanwise permeabilities set to unity, i.e. 
$K_y=K_z$. Note that the permeability tensor is symmetric, and so an eigenvalue decomposition can be used to identify its principal values and directions (or axes). The assumed form of the permeability tensor implies that its principal directions align with the streamwise, wall-normal and spanwise directions of the flow. Furthermore, we omit the nonlinear Forchheimer correction term, 
$\boldsymbol{\mathsf{F}}$, that is used to account for inertial effects in flows through porous media (Ochoa-Tapia & Whitaker Reference Ochoa-Tapia and Whitaker1995a,Reference Ochoa-Tapia and Whitakerb; Whitaker Reference Whitaker1999; Breugem et al. Reference Breugem, Boersma and Uittenbogaar2006). This assumption is made for several reasons. First, the Forchheimer term is expected to be small for the low values of 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}}$, 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$ and 
$\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$ considered below (
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} \lesssim 11$ and 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = \sqrt {\smash {K_z^+}\vphantom {K_x^+}} \lesssim 1$ in all cases). Second, neglecting the Forchheimer term ensures consistency with the Brinkman model employed by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) in the numerical simulations that serve as a basis for comparison for our model predictions. Third, including the Forchheimer term would introduce several additional parameters that are dependent on the pore-scale geometry and Reynolds number. Finally, resolvent analysis probes the linear forcing–response characteristics of the governing equations and therefore, only a linearized version of the Forchheimer term would be retained explicitly in the analysis. This is equivalent to simply considering a substrate with a lower apparent permeability (see e.g. Zampogna & Bottaro Reference Zampogna and Bottaro2016). We also assume that the porosity of permeable substrates is 
$\varepsilon \approx 1$ to maximize any potential drag reduction (Abderrahaman-Elena & García-Mayoral Reference Abderrahaman-Elena and García-Mayoral2017). Finally, since we are primarily interested in structures that are much larger than the characteristic length scale of the porous medium (i.e. NW cycle and KH rollers), we neglect the subfilter scale stresses (Breugem et al. Reference Breugem, Boersma and Uittenbogaar2006). With these assumptions, (2.1a) and (2.1b) can be expressed as
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\boldsymbol{u}\boldsymbol{u} \right) ={-}\boldsymbol{\nabla}(p) + \frac{1}{\textit{Re}_{\tau}}\boldsymbol{\nabla}^2 \boldsymbol{u} -\frac{1}{\textit{Re}_{\tau}}\boldsymbol{\mathsf{K}}^{{-}1}\boldsymbol{u}\end{equation}and
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0,\end{equation}
where the 
$\langle \cdot \rangle$ notation has been eliminated for simplicity. The unobstructed fluid domain is characterized by infinite permeability. In this region, the permeability term goes to zero and (2.2a) reduces to the standard Navier–Stokes momentum equation.
Figure 1 shows the symmetric channel flow configuration considered in this study. The unobstructed region corresponds to 
$y\in [-1,1]$. The regions occupied by the permeable substrates correspond to 
$y\in [-(1+H),-1)$ and 
$y\in (1,1+H]$. The height of the permeable substrate is 
$H$. Note that all lengths are normalized by the channel half-height.
2.2. Modified resolvent analysis
The resolvent formulation for wall-bounded turbulent flows proposed by McKeon & Sharma (Reference McKeon and Sharma2010) – and employed in several recent flow control studies (Luhar et al. Reference Luhar, Sharma and McKeon2014, Reference Luhar, Sharma and McKeon2015; Nakashima et al. Reference Nakashima, Fukagata and Luhar2017; Toedtli et al. Reference Toedtli, Luhar and McKeon2019; Chavarin & Luhar Reference Chavarin and Luhar2020) – is extended to account for the presence of permeable substrates as follows. For an extended discussion of resolvent analysis and its applications, the reader is referred to McKeon (Reference McKeon2017).
Construction of the modified resolvent operators begins with a standard Reynolds decomposition of the simplified VANS equations in (2.2). The velocity and pressure fields are decomposed into a time-averaged component (denoted by 
$\overline {(\cdot )}$) and a fluctuating component about this average (denoted by 
${(\cdot )^\prime }$). Under this decomposition, the velocity field is expressed as 
$\boldsymbol {u}(t,{\boldsymbol {x}})=\bar {\boldsymbol {U}}({\boldsymbol {x}})+\boldsymbol {u}^\prime (t,{\boldsymbol {x}})$ and the pressure field is expressed as 
$p(t,{\boldsymbol {x}}) = \bar {P}({\boldsymbol {x}}) + p^\prime (t,{\boldsymbol {x}})$. Note that 
$\bar {\boldsymbol {U}}({\boldsymbol {x}}) = [\bar {U}(y),0,0]^\textrm {T}$ represents the turbulent mean profile. Next, the velocity and pressure fluctuation are Fourier-transformed in the homogeneous streamwise and spanwise directions as well as in time as follows:
\begin{equation} \begin{bmatrix} \boldsymbol{u}^\prime(t,\boldsymbol{x}) \\ p^\prime(t,\boldsymbol{x}) \end{bmatrix} = {\iiint}\begin{bmatrix} \boldsymbol{u}_{\boldsymbol{{\kappa}}}(y) \\ p_{\boldsymbol{{\kappa}}}(y) \end{bmatrix}\exp(-\textrm{i}\omega t+\textrm{i}\kappa_x x +\textrm{i}\kappa_z z)\,\textrm{d}{\omega}\,\textrm{d} {\kappa_x}\,\textrm{d}{\kappa_z}. \end{equation}
In the expression above, 
$\kappa _x$ is the streamwise wavenumber, 
$\kappa _z$ is the spanwise wavenumber and 
$\omega$ is the frequency. The Fourier coefficients for the velocity and pressure field at a given wavenumber–frequency combination, 
$\boldsymbol {{\kappa }} = (\kappa _x,\kappa _z,\omega )$, are denoted 
$\boldsymbol {u}_{\boldsymbol {{\kappa }}}$ and 
$p_{\boldsymbol {{\kappa }}}$. Under the Fourier transform, the equations in (2.2) can be expressed compactly as
\begin{equation} \begin{bmatrix} {\boldsymbol{u}}_{\boldsymbol{{\kappa}}} \\ {p}_{\boldsymbol{{\kappa}}} \end{bmatrix} = \left(-\textrm{i}\omega \left[\begin{array}{@{}cc@{}} {\boldsymbol{\mathsf{I}}} & \\ & 0\end{array}\right] -\left[\begin{array}{@{}cc@{}} \boldsymbol{\mathcal{L}}_{\boldsymbol{{\kappa}}} & -\tilde{\boldsymbol{\nabla}} \\ -\tilde{\boldsymbol{\nabla}}^T & 0\end{array}\right]\right)^{{-}1} \left[\begin{array}{@{}c@{}} {\boldsymbol{\mathsf{I}}}\\ 0\end{array}\right]\boldsymbol{f}_{\boldsymbol{{\kappa}}} = \boldsymbol{\mathcal{H}}_{\boldsymbol{{\kappa}}}\boldsymbol{f}_{\boldsymbol{{\kappa}}}.\end{equation}In the expression above, the first row represents the momentum equations and the second row represents the continuity constraint. The operator
\begin{align}
\boldsymbol{\mathcal{L}}_{\boldsymbol{{\kappa}}} = -
\textit{Re}_{\tau}^{-1} \boldsymbol{\mathsf{K}}^{-1} +
\begin{bmatrix} -\textrm{i}\kappa_x \bar{U} +
\textit{Re}_{\tau}^{-1}\tilde{\boldsymbol{\nabla}}^2 & -
\dfrac{\textrm{d}{\bar{U}}}{\textrm{d}y} & 0
\\ 0 & -\textrm{i}\kappa_x \bar{U} +
\textit{Re}_{\tau}^{-1}\tilde{\boldsymbol{\nabla}}^2
& 0 \\ 0 & 0 & -\textrm{i}\kappa_x \bar{U} +
\textit{Re}_{\tau}^{-1}\tilde{\boldsymbol{\nabla}}^2 \end{bmatrix}
\end{align}
represents the linear dynamics in (2.2) and 
$\boldsymbol {f}_{\!\boldsymbol {{\kappa }}}$ is the Fourier coefficient for the nonlinear terms. The differential operator 
$\tilde {\boldsymbol {\nabla }}$ is defined as 
$\tilde {\boldsymbol {\nabla }}=(\textrm{i}\kappa _x,\partial /\partial y,\textrm{i}\kappa _z)$, and so the symbols 
$\tilde {\boldsymbol {\nabla }}^T$ and 
$\tilde {\boldsymbol {\nabla }}$ essentially represent the Fourier-transformed divergence and gradient operators, respectively. Similarly, the Laplacian is defined as 
$\tilde {\boldsymbol {\nabla }}^2=(-\kappa _x^2-\kappa _z^2+\partial ^2/\partial y^{2})$. Note that the velocity and pressure response at a given wavenumber–frequency combination constitute a travelling wave flow field with streamwise wavelength 
$\lambda _x=2{\rm \pi} /\kappa _x$ and spanwise wavelength 
$\lambda _z=2{\rm \pi} /\kappa _z$ that is moving downstream at speed 
$c=\omega /\kappa _x$. The transfer function that maps the nonlinear forcing 
$\boldsymbol {f}_{\!\boldsymbol {{\kappa }}}$ to the velocity and pressure response 
$[\boldsymbol {u}_{\boldsymbol {{\kappa }}},{p}_{\boldsymbol {{\kappa }}}]^\textrm {T}$ in (2.4) is the resolvent operator, 
$\boldsymbol{\mathcal{H}}_{\boldsymbol {{\kappa }}}$. The central difference between this resolvent operator for channel flow over permeable substrates and its smooth wall counterpart is the inclusion of the Darcy permeability tensor 
$\boldsymbol{\mathsf{K}}$ in (2.5). Note that the resistance exerted by the Darcy term is linear with respect to 
$\boldsymbol {u}$.
As detailed in prior studies (McKeon & Sharma Reference McKeon and Sharma2010; Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013; Luhar et al. Reference Luhar, Sharma and McKeon2014, Reference Luhar, Sharma and McKeon2015), a singular value decomposition (SVD) of the discretized resolvent operator (see § 2.5) is used to identify a set of orthonormal forcing and response modes that are ordered based on their forcing–response gain under an 
$L^2$ energy norm. To enforce this norm, the discretized resolvent operator in (2.4) is scaled as follows:
\begin{equation} \left[\boldsymbol{\mathsf{W}}_{\boldsymbol{u}}\enspace 0\right] \begin{bmatrix} \boldsymbol{u}_{\boldsymbol{{\kappa}}} \\ {p}_{\boldsymbol{{\kappa}}} \end{bmatrix} = \left(\left[\boldsymbol{\mathsf{W}}_{\boldsymbol{u}} \enspace 0\right]\ \boldsymbol{\mathcal{H}}_{\boldsymbol{{\kappa}}} \boldsymbol{\mathsf{W}}_{\boldsymbol{f}}^{-1}\right) \boldsymbol{\mathsf{W}}_{\boldsymbol{f}}\boldsymbol{f}_{\boldsymbol{{\kappa}}} \end{equation}or
\begin{equation} \boldsymbol{\mathsf{W}}_{\boldsymbol{u}}\boldsymbol{u}_{\boldsymbol{{\kappa}}} = \boldsymbol{\mathcal{H}}^w_{\boldsymbol{{\kappa}}} \boldsymbol{\mathsf{W}}_{\boldsymbol{f}}\boldsymbol{f}_{\boldsymbol{{\kappa}}}. \end{equation}
Here, 
$\boldsymbol{\mathsf{W}}_{\boldsymbol {u}}$ and 
$\boldsymbol{\mathsf{W}}_{\boldsymbol {u}}$ incorporate numerical quadrature weights for the entire domain spanning 
$y\in [-(1+H),(1+H)]$. With this weighting, the SVD of the scaled resolvent,
\begin{equation} \boldsymbol{\mathcal{H}}^w_{\boldsymbol{{\kappa}}}= \sum_{m} \psi_{\boldsymbol{{\kappa}},m} \sigma_{\boldsymbol{{\kappa}},m} \phi^{{\ast}}_{\boldsymbol{{\kappa}},m}, \end{equation}where
\begin{equation} \sigma_{\boldsymbol{{\kappa}},1}>\sigma_{\boldsymbol{{\kappa}},2}>\cdots>0, \quad \psi_{\boldsymbol{{\kappa}},l}^*\psi_{\boldsymbol{{\kappa}},m} = \delta_{lm}, \quad \phi_{\boldsymbol{{\kappa}},l}^*\phi_{\boldsymbol{{\kappa}},m} = \delta_{lm} \end{equation}
yields forcing modes 
$\boldsymbol {f}_{\!\boldsymbol {{\kappa }},m}=\boldsymbol{\mathsf{W}}_{\boldsymbol {f}}^{-1}\phi _{\boldsymbol {{\kappa }},m}$ and velocity responses 
$\boldsymbol {u}_{\boldsymbol {{\kappa }},m}={\boldsymbol{\mathsf{W}}_{\boldsymbol {u}}^{-1}}\psi _{\boldsymbol {{\kappa }},m}$ that have unit energy when integrated across the entire domain spanning 
$y\in [-(1+H),(1+H)]$. In other words, this scaling ensures that
\begin{equation} \int_{-(1+H)}^{(1+H)} \boldsymbol{u}^*_{\boldsymbol{{\kappa}},l} \boldsymbol{u}^{}_{\boldsymbol{{\kappa}},m} \,{\textrm{d} y} = \delta_{lm},\quad \int_{-(1+H)}^{(1+H)} \boldsymbol{f}^*_{\boldsymbol{{\kappa}},l} \,\boldsymbol{f}^{}_{\boldsymbol{{\kappa}},m} \,{\textrm{d}y} = \delta_{lm}. \end{equation}
In (2.7) and (2.8a,b), the superscript 
$(\cdot )^*$ denotes a conjugate transpose.
A major contribution of the resolvent framework lies in the finding that the forcing–response transfer function tends to low rank at 
$\boldsymbol {{\kappa }}$ combinations that are energetic in wall-bounded turbulent flows. Often, the first singular value is an order of magnitude larger than subsequent singular values, 
$\sigma _{\boldsymbol {{\kappa }},1} \gg \sigma _{\boldsymbol {{\kappa }},2}>\cdots$, and so the resolvent operator can be well approximated using a rank-1 truncation after the SVD (McKeon & Sharma Reference McKeon and Sharma2010; Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013),
\begin{equation} \boldsymbol{\mathcal{H}}^w_{\boldsymbol{{\kappa}}} \approx \psi_{\boldsymbol{{\kappa}},1} \sigma_{\boldsymbol{{\kappa}},1} \phi^{{\ast}}_{\boldsymbol{{\kappa}},1}. \end{equation}
The expressions in (2.6)–(2.9) show that forcing in the direction of the first forcing mode 
$\boldsymbol {f}_{\!\boldsymbol {{\kappa }},1}=\boldsymbol{\mathsf{W}}_{\boldsymbol {f}}^{-1}\phi _{\boldsymbol {{\kappa }},1}$ generates a velocity response 
$\boldsymbol {u}_{\boldsymbol {{\kappa }},1}={\boldsymbol{\mathsf{W}}_{\boldsymbol {u}}^{-1}}\psi _{\boldsymbol {{\kappa }},1}$ that is amplified by factor 
$\sigma _{\boldsymbol {{\kappa }},1}$. Put another way, a forcing of the form 
$\boldsymbol {f}_{\!\boldsymbol {{\kappa }},1}$ to the unscaled resolvent operator in (2.4) generates a velocity and pressure response 
$\sigma _{\boldsymbol {{\kappa }},1}[\boldsymbol {u}_{\boldsymbol {{\kappa }},1},p_{\boldsymbol {{\kappa }},1}]^\textrm {T}$. Under the 
$L^2$ scaling used here, 
$\sigma _{\boldsymbol {{\kappa }},1}^2$ is a measure of energy amplification.
Recent modelling efforts for active and passive flow control techniques show that specific rank-1 modes serve as useful surrogates for the dynamically important NW cycle. Specifically, the ability of a control technique to suppress the forcing–response gain for modes with wavenumber–frequency combinations corresponding to 
$\lambda _x^+ \approx 10^3$, 
$\lambda _z^+ \approx 10^2$ and 
$c^+ \approx 10$ (i.e. similar to the length and velocity scales associated with NW streaks) has been shown to be a useful predictor of drag reduction performance (Luhar et al. Reference Luhar, Sharma and McKeon2014; Nakashima et al. Reference Nakashima, Fukagata and Luhar2017; Chavarin & Luhar Reference Chavarin and Luhar2020). Building on these prior efforts, in this study we evaluate the effect of anisotropic permeable substrates on the rank-1 resolvent mode that serves as a surrogate for the NW cycle. A reduction in gain for this mode relative to the smooth-wall value is interpreted as mode suppression, which is indicative of drag reduction. In addition, we also test whether the permeable substrates lead to an increase in principal singular values for spanwise-coherent modes (e.g. with 
$\kappa _z = 0$) that resemble KH rollers. Since we only consider the rank-1 truncation shown in (2.9) for the remainder of this paper, we drop the additional subscript 
$1$ to simplify notation.
2.3. Boundary and interface conditions
As will be discussed in § 2.5 below, the resolvent operator is discretized using spectral discretization and rectangular block matrices as described by Aurentz & Trefethen (Reference Aurentz and Trefethen2017). This approach enables us to use two different sets of equations in the unobstructed region and porous domain (i.e. without and with the permeability term), and couple the two via appropriate interfacial conditions. As shown in figure 1, the unobstructed channel corresponds to the region corresponding to 
$y\in [-1,1]$ and the upper and lower permeable regions correspond to 
$y\in (1,1+H]$ and 
$y\in [-(1+H),-1)$, respectively. At the lower and upper substrate walls, 
$y=\pm (1+H)$, we apply no-slip boundary conditions. At the interfaces between the porous medium and the unobstructed flow, 
$y = \pm 1$, we impose continuity in all three components of velocity and pressure. We also impose continuity in the streamwise and spanwise shear at the interface. These boundary conditions can be summarized as follows:
\begin{gather} \boldsymbol{u} = 0 \quad\text{at}\ y ={\pm} (1+H), \end{gather}
\begin{gather}\boldsymbol{u}|_{y_+} =\boldsymbol{u}|_{y_-}\quad\text{and}\quad p |_{y_+} = p |_{y_-} \quad\text{at}\ y ={\pm} 1, \end{gather}
\begin{gather}\left.\frac{\partial u}{\partial y}\right\vert_{y_+} = \frac{1}{\varepsilon}\left.\frac{\partial u}{\partial y} \right\vert_{y_-} \quad\text{and}\quad \left.\frac{\partial w}{\partial y}\right\vert_{y_+} = \frac{1}{\varepsilon}\left.\frac{\partial w}{\partial y}\right\vert_{y_-} \quad \text{at} \ y ={-}1, \end{gather}
\begin{gather}\left.\frac{\partial u}{\partial y}\right\vert_{y_-} = \frac{1}{\varepsilon}\left. \frac{\partial u}{\partial y}\right\vert_{y_+} \quad\text{and}\quad \left.\frac{\partial w}{\partial y}\right\vert_{y_-} =\frac{1}{\varepsilon}\left. \frac{\partial w}{\partial y}\right\vert_{y_+} \quad \text{at} \ y = 1. \end{gather}
In the expressions above, 
$y_+$ and 
$y_-$ denote values taken on either side of the porous substrate-unobstructed flow interface. Note that the boundary conditions shown in (2.10c,d) do not penalize momentum transfer into the permeable substrate. These conditions also assume that the effective viscosity in the porous medium is 
$\tilde {\nu } = \nu /\varepsilon$. As per the discussion presented in Abderrahaman-Elena & García-Mayoral (Reference Abderrahaman-Elena and García-Mayoral2017) and Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019), these assumptions are reasonable for a highly connected medium with high porosity, 
$\varepsilon \approx 1$. For a poorly connected medium, previous studies show that a stress jump boundary condition may be more appropriate, and that more complex effective viscosity models may be needed (see e.g. Ochoa-Tapia & Whitaker Reference Ochoa-Tapia and Whitaker1995a; Minale Reference Minale2014). Also keep in mind that the equations above imply a sharp transition between the porous medium and the unobstructed fluid. Previous studies employing the VANS equations have typically assumed the existence of a finite transition zone between the homogeneous porous and homogeneous fluid regions (Ochoa-Tapia & Whitaker Reference Ochoa-Tapia and Whitaker1995a; Breugem et al. Reference Breugem, Boersma and Uittenbogaar2006). This continuous approach was also considered for the present study, but it led to poorer numerical convergence in singular values. Moreover, the converged values showed a significant dependence on the size of the transition zone.
Finally, since we are interested in permeable substrates that have the highest potential for drag reduction, we assume a porosity of 
$\varepsilon =1$. For this value of porosity, the boundary conditions shown in (2.10) are similar to those used in previous high-fidelity simulations (Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). Emulation of the channel configuration and boundary conditions used by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) allows for a direct comparison between model predictions and simulation results.
2.4. Mean velocity profile
As shown in (2.5), construction of the resolvent operator requires knowledge of the mean velocity profile 
$\bar {U}(y)$. Here, we use two different sets of mean velocity profiles to generate model predictions: (i) mean profiles obtained in DNS by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019); and (ii) mean profiles generated using a synthetic eddy viscosity profile. The DNS mean profile dataset consists of 22 different cases: eight profiles for 
$\phi _{xy} = 3.6$; seven profiles for 
$\phi _{xy} = 5.5$; and seven profiles for 
$\phi _{xy} = 11.4$. For all 22 cases tested in DNS, the ratio between the wall-normal and spanwise permeability length scales was 
$\phi _{yz}= \sqrt {\smash {K_y^+}\vphantom {K_x^+}}/\sqrt {\smash {K_z^+}\vphantom {K_x^+}} =1$. Each of the DNS profiles contained 153 points spanning the unobstructed region, 
$y\in [-1,1]$. Within the permeable medium, the analytical solutions developed by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) were used. These mean profiles were interpolated onto our Chebyshev grid using the modified Akima interpolation algorithm.
For each of the 22 different DNS cases, the synthetic profiles were generated as follows. The modified governing equation (2.2) was Reynolds-averaged to yield the following equation for the streamwise mean flow:
\begin{equation} \frac{1}{\textit{Re}_{\tau}}\frac{\textrm{d}^2\bar{U}}{\textrm{d} y^{2}}- \frac{\bar{U}}{\textit{Re}_{\tau} K_x} -\frac{\textrm{d}(\overline{u^\prime v^\prime})}{\textrm{d} y} = \frac{\textrm{d}\bar{p}}{\textrm{d} x}, \end{equation}
where 
$K_x$ is the streamwise component of the permeability tensor. The Reynolds shear stress term in (2.11) was estimated using a modified version of the eddy viscosity formulation proposed by Reynolds & Tiederman (Reference Reynolds and Tiederman1967) as follows:
\begin{equation} \nu_e =\left\lbrace \begin{array}{@{}ll} \dfrac{1}{2}\left\lbrack 1+\left\{\dfrac{c_2\textit{Re}_{\tau}} {3}(2y-y^2)(3-4y+2y^2)\right.\right. & \\ \qquad\times \left.\left.\left(1-\exp\left((|y-1|-1)\dfrac{\textit{Re}_{\tau}}{c_1}\right)\right)\right\}^2 \right\rbrack^{1/2}- \dfrac{1}{2}, & |y| \leqslant 1, \\ 0, & |y| > 1. \end{array}\right.\end{equation}
Here, the eddy viscosity has been normalized by the fluid viscosity 
$\nu$. To generate the synthetic profiles, we used the values 
$c_1 = 46.2$ and 
$c_2 = 0.61$, which were identified by Moarref & Jovanović (Reference Moarref and Jovanović2012) as yielding the best fit to the mean profiles obtained in smooth wall DNS at 
$Re_{\tau } \approx 180$. Thus, the Reynolds shear stress term was modelled using a standard smooth-wall eddy viscosity profile in the unobstructed region of the channel and assumed to be zero in the permeable substrate. In other words, this eddy viscosity model assumes that there is no turbulence penetration into the porous medium. In reality, the penetration of turbulent cross-flows into the permeable substrate will depend on the spanwise and wall-normal components of permeability, and so the model above is only valid for small 
$K_y$ and 
$K_z$. Finally, (2.11) and (2.12) were combined to yield the following equation for the mean velocity profile:
\begin{equation} \frac{1}{\textit{Re}_{\tau}}\left((1+\nu_e)\frac{\textrm{d}^{2}}{\textrm{d}y^{2}} + \frac{\textrm{d}\nu_e}{\textrm{d}y}\frac{\textrm{d}}{\textrm{d}y}- \frac{1}{K_x}\right)\bar{U} =\frac{\textrm{d}\bar{p}}{\textrm{d} x}, \end{equation}
which was solved numerically to yield 
$\bar {U}(y)$.
We recognize that the procedure outlined above to estimate the synthetic mean profile involves significant assumptions. However, as we show below, the interfacial slip velocities generated using this model are in good agreement with the DNS mean profiles in the initial drag reduction regime over anisotropic permeable substrates, i.e. until the point of performance degradation. Moreover, the resolvent-based predictions obtained using these profiles are also similar to those computed using the DNS mean profiles. Note that the resolvent operator shown in (2.4) is formulated using the full VANS equations inside the permeable substrate, even though the eddy viscosity used to compute the synthetic mean profiles is set to zero in this zone. In other words, the fluctuating velocity field associated with resolvent modes is allowed to penetrate into the porous substrate.
2.5. Numerical methods
The resolvent operator and the equations used to synthesize the mean velocity profile are discretized in the wall-normal direction using spectral discretization methods involving rectangular block operators, as described by Aurentz & Trefethen (Reference Aurentz and Trefethen2017). Each differential operator is discretized using Chebyshev polynomials and the resulting matrices are rectangular. The size of these matrices is 
$[N \times N+n]$ where 
$n$ represents the dimension of the operator null space. The overall block operator is made square by appending 
$n$ boundary conditions. By using these block operators, we can deal with a system of boundary value problems coupled through boundary conditions. For a more thorough discussion of these methods readers are directed to Aurentz & Trefethen (Reference Aurentz and Trefethen2017).
As shown in figure 1, for our problem configuration the channel is separated into three regions: the lower permeable substrate spanning 
$y \in [-(H+1),-1)$; the free channel spanning 
$y \in [-1,1]$; and the upper permeable substrate spanning 
$y \in (1,H+1]$. Although the configuration is symmetric across the channel centreline, we consider the entire channel to retain resolvent modes that are both symmetric and anti-symmetric across the channel centreline (Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013; Luhar et al. Reference Luhar, Sharma and McKeon2015). Each of the three regions is discretized using 
$N$ Chebyshev points and so the total number of Chebyshev points for the channel is 
$3N$. A grid convergence study showed that the normalized change in singular values is of 
$O(10^{-6})$ for grid sizes 
$N=80$, 
$N=112$, and 
$N=200$. Each resolvent mode computation takes approximately 0.7 s for 
$N=80$, approximately 2.0 s for 
$N=112$ and 7.0 s for 
$N=200$. The results presented below were generated using 
$N=112$, which corresponds to a total of 
$3N=336$ grid points across the entire channel.
3. Results and discussion
In this section, we compare resolvent-based predictions for the NW mode (§ 3.2) and the emergence of energetic spanwise rollers (§ 3.3) with DNS observations. Before that, we briefly compare the mean velocity profiles obtained from DNS with the synthetic profiles generated using the eddy viscosity model (§ 3.1).
3.1. Mean velocity profiles
Figure 2 compares the interfacial slip velocity 
$\bar {U}_s^+$ predicted using the synthetic eddy viscosity profile (2.12) with those obtained in DNS by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) for each of the 22 different cases being considered here, with anisotropy ratios 
$\phi _{xy} = 3.6$, 
$5.5$ and 
$11.4$. The slip velocities for the synthetic profiles all collapse together onto a straight line corresponding to 
$\bar {U}_s^+ \approx \sqrt {\smash {K_x^+}\vphantom {K_x^+}}$. In other words, the slip velocity only depends on the streamwise component of permeability, which is consistent with the assumptions outlined in § 2.4. The synthetic mean profiles do not account for the effect of wall-normal or spanwise permeability, which are likely to determine the extent to which turbulence penetrates into the permeable substrate.
Figure 2. Predicted slip velocity at the porous interface as a function of streamwise permeability 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}}$ for anisotropy ratios (a) 
$\phi _{xy} = 3.6$, (b) 
$\phi _{xy} = 5.5$ and (c) 
$\phi _{xy} = 11.4$. Black symbols show DNS results from Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019). White symbols show predictions made using (2.13). The vertical dashed lines show where the wall-normal permeability is 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \approx 0.4$. This is roughly the threshold above which drag reduction performance deteriorates in DNS.
The synthetic slip velocities (white symbols) show close agreement with DNS results (black symbols) for low streamwise permeabilities, but begin to deviate from DNS results at higher 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}}$. Moreover, the threshold values of 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}}$ above which the synthetic predictions begin to deviate from DNS results depend on the anisotropy ratio, 
$\phi _{xy}$. As an example, for 
$\phi _{xy} = 3.6$, synthetic 
$\bar {U}_s^+$ begin deviating from DNS results for 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} \gtrsim 2$ (see figure 2a). For 
$\phi _{xy} = 11.4$, the predicted slip velocities deviate significantly from DNS results for 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} \gtrsim 6$ (see figure 2c). Closer inspection of this trend shows that the synthetic slip velocities deviate from DNS data above a constant threshold value for the wall-normal permeability for all anisotropy ratios, 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$ (vertical dashed lines in figure 2). This value of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$ corresponds closely to the conditions in which KH-type rollers appear over the permeable substrates in the simulations of Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019). The appearance of these rollers is likely to generate significant interfacial turbulence that penetrates the porous medium. The eddy viscosity model used to generate the synthetic mean profiles does not account for these effects. The eddy viscosity model shown in (2.12) assumes that no turbulence penetrates into the porous medium, and that the turbulence in the unobstructed region remains similar to that over a smooth wall regardless of the value of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$ and 
$\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$.
For cases in which the interfacial slip velocities agree, the synthetic mean profiles are similar to the DNS profiles across the entire channel (data not shown here for brevity). For these cases, resolvent-based predictions are not very sensitive to the choice of mean profile. Moreover, as we show in § 3.3, resolvent analysis is able to predict the emergence of high-gain KH rollers over permeable substrates as the wall-normal permeability increases beyond 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$. These observations indicate that resolvent analysis carried out using the synthetic mean profile can be a useful design tool for the design of passive flow control using anisotropic permeable substrates.
3.2. Near-wall resolvent mode
Previous studies have shown that the resolvent mode with streamwise wavelength 
$\lambda _x^+ = 10^3$, spanwise wavelength 
$\lambda _z^+ = 10^2$ and phase speed 
$c^+ = (\omega /\kappa _x)^+ = 10$ serves as a useful surrogate for the dynamically important NW cycle characterized by the presence of quasi-streamwise vortices and alternating high- and low-speed streaks (McKeon & Sharma Reference McKeon and Sharma2010; Sharma & McKeon Reference Sharma and McKeon2013). The forcing–response gain (i.e. singular value) for this resolvent mode has also proven to be a useful predictor of control performance for both active (Luhar et al. Reference Luhar, Sharma and McKeon2014; Nakashima et al. Reference Nakashima, Fukagata and Luhar2017; Toedtli et al. Reference Toedtli, Luhar and McKeon2019) and passive (Luhar et al. Reference Luhar, Sharma and McKeon2015; Luhar, Sharma & McKeon Reference Luhar, Sharma and McKeon2016; Chavarin & Luhar Reference Chavarin and Luhar2020) techniques. Here, we evaluate whether it can serve as a useful reduced-complexity tool for the evaluation of anisotropic permeable substrates for passive drag reduction. Since the resolvent-based predictions are evaluated against DNS results from Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) carried out at 
$\textit {Re}_{\tau } = 180$, the resolvent mode with streamwise wavenumber 
$\kappa _x = 2{\rm \pi} \textit {Re}_{\tau }/10^3 \approx 1.1$, spanwise wavenumber 
$\kappa _z = 2{\rm \pi} \textit {Re}_{\tau }/10^2 \approx 11$ and frequency 
$\omega = 10 \kappa _x \approx 11$ is used as a surrogate for the NW cycle.
Figure 3(a) shows the forcing-gain for the NW resolvent mode over the porous substrate 
$\sigma _{\boldsymbol {{\kappa }},p}$ normalized by the smooth-wall value 
$\sigma _{\boldsymbol {{\kappa }},s}$ at 
$\textit {Re}_{\tau } = 180$ for all 22 substrates tested by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019). Resolvent-based gain predictions (black symbols) are compared with the drag reduction measured in DNS (grey symbols). Following Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019), the outward shift in the logarithmic region of the mean profile is used to quantify the change in drag: 
$\Delta \bar {U}^+ > 0$ denotes a decrease in drag and 
$\Delta \bar {U}^+ < 0$ denotes an increase in drag. For all three anisotropy ratios, there is qualitative agreement between NW mode suppression and drag reduction. For a given value of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$, materials with higher anisotropy ratios are predicted to yield greater suppression of the NW resolvent mode, which is consistent with the drag reduction trends observed in DNS. Furthermore, mode suppression and drag reduction both increase monotonically up to a value of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \approx 0.4$. Agreement between drag reduction and mode suppression is particularly good for substrates with anisotropy ratios 
$\phi _{xy} = 3.6$ and 
$5.5$. For 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$, the normalized gain begins to increase again for the substrates with 
$\phi _{xy} = 3.6$ and 
$5.5$, which is consistent with DNS drag reduction trends. The DNS results show an increase in drag relative to smooth-wall values for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.55$ for the substrates with 
$\phi _{xy} = 3.6$ and 
$5.5$. Resolvent-based predictions show that the gain of the NW mode over the permeable substrate exceeds the smooth-wall value beyond 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.7$ for these substrates. For substrates with the higher anisotropy ratio, 
$\phi _{xy}= 11.4$, there are visible discrepancies between the predicted change in NW mode gain and the measured drag reduction, particularly for higher values of the wall normal permeability, 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. For instance, the increase in normalized gain beyond 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$ is less pronounced for the material with 
$\phi _{xy}= 11.4$. For this substrate, DNS results show an increase in drag for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.7$. However, the resolvent predictions do not show an increase in gain beyond the smooth-wall value for the conditions tested (up to 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \approx 1$).
Figure 3. Predictions for the NW resolvent mode. (a) Comparison between normalized mode gain (
$\sigma _{\boldsymbol {{\kappa }},p}/\sigma _{\boldsymbol {{\kappa }},s}$, black symbols) and drag reduction observed in DNS (
$\Delta \bar {U}^+$, light grey symbols) by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019), plotted as a function of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. The resolvent-based predictions shown in this panel are obtained using the velocity profiles from DNS. (b) Comparison between normalized gains predicted using the DNS mean profiles (black symbols) and the synthetic mean profiles (white symbols) plotted as a function of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. (c) Comparison between normalized gains predicted using the DNS mean profiles (black symbols) and the synthetic mean profiles (white symbols) plotted as a function of 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. The black and red dashed lines show linear fits to the initial decrease in normalized gain for the DNS mean and synthetic mean predictions, respectively. For all panels, the 
$\square$ symbols represent substrates with 
$\phi _{xy} = 3.6$, 
$\triangle$ symbols represent substrates with 
$\phi _{xy}= 5.5$ and 
symbols represent substrates with 
$\phi _{xy}=11.4$.
Note that the observed disagreement between NW mode gain and drag reduction for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$ is not surprising. As discussed in prior studies (Abderrahaman-Elena & García-Mayoral Reference Abderrahaman-Elena and García-Mayoral2017; Gómez-de-Segura et al. Reference Gómez-de-Segura, Sharma and García-Mayoral2018; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019), the initial decrease in drag over anisotropic permeable substrates depends on the shift in the virtual origins for the streamwise mean flow and the turbulent cross-flow fluctuations produced by the quasi-streamwise vortices associated with the NW cycle. As we show below, this effect is reproduced by the NW resolvent mode. However, the increase in drag beyond 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \approx 0.4$ is associated with the emergence of spanwise rollers resembling KH vortices. This phenomenon cannot be captured by resolvent modes that serve as surrogates for the streamwise-elongated structures constituting the NW cycle. We use the resolvent framework to test for the emergence of energetic spanwise rollers in § 3.3.
Figure 3(b) shows normalized gain predictions made using the DNS mean profile (black symbols) and the synthetic mean profile (white symbols) for the NW resolvent mode, plotted as a function of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. For values of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \lesssim 0.4$, the normalized gain predictions are in close agreement. For this range of wall-normal permeabilities, the normalized gain predictions made using the synthetic mean profile are within 6 % of those made using the DNS mean profile. However, for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$, the synthetic profile predictions show a continued decrease in NW mode gain, while the DNS profile predictions show an increase in gain. Discrepancy between the predictions for higher values of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$ is consistent with the slip velocities shown in figure 2. The synthetic profile agrees with the DNS mean profile for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \lesssim 0.4$, but not beyond this value. As discussed earlier, this is because higher wall-normal permeabilities give rise to spanwise rollers resembling KH vortices in the simulations. The effect of the additional interfacial turbulence generated by these rollers is not captured in the eddy viscosity model (2.12). Instead, the eddy viscosity model assumes smooth-wall-like turbulence regardless of the value of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$ and 
$\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$.
Figure 3(c) shows normalized gain predictions made using the DNS mean profile (black symbols) and the synthetic mean profile (white symbols) for the NW resolvent mode, plotted as a function of 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. Since the permeable substrate model used here has identical wall-normal and spanwise permeabilities, the difference between streamwise and wall-normal permeability is also equal to the difference between streamwise and spanwise permeabilities, i.e. 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_y^+}\vphantom {K_x^+}} = \sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_z^+}\vphantom {K_x^+}}$. Slip length arguments (Luchini et al. Reference Luchini, Manzo and Pozzi1991; Abderrahaman-Elena & García-Mayoral Reference Abderrahaman-Elena and García-Mayoral2017) suggest that the initial decrease in drag over anisotropic permeable substrates is proportional to the difference between the streamwise and spanwise permeability length scales, which determine the virtual origins felt by the mean streamwise flow and turbulent cross-flow fluctuations, respectively. In other words, we expect 
$\Delta \bar {U}^+ \propto \sqrt {\smash {K_x^+}\vphantom {K_x^+}}-\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$. These arguments are supported by the simulation results of Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019), which show a linear decrease in drag that is proportional to 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}}-\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$ for substrates with small permeabilities. The NW mode gains shown in figure 3(c) agree well with this model. Specifically, the normalized NW-mode gains initially decrease linearly with 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_y^+}\vphantom {K_x^+}} ({=}\sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_z^+}\vphantom {K_x^+}})$ across all anisotropy ratios. Moreover, there is little difference between the predictions made using the synthetic and DNS mean profiles. Linear relationships fitted to the initial decrease in normalized gain are nearly identical for the predictions based on the DNS profiles (dashed black lines) and those based on the synthetic mean profiles (dashed red lines). Note that the predictions based on the synthetic mean profile all collapse together. This suggests that NW mode behaviour is self-similar when appropriately normalized. Predictions made using the DNS profiles diverge from this self-similarity when 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$.
Figure 4 provides physical insight into the observed changes in mode gain. The structure of the NW resolvent mode for the smooth-wall case is shown in figure 4(a). By construction, the mode structure shows alternating regions of positive and negative velocity fluctuations with spanwise wavelength 
$\lambda _z^+ = 10^2$. The wall-normal and spanwise velocity fluctuations show the presence of counter-rotating streamwise vortices. Regions of downwelling (
$v^\prime$ towards wall) coincide with regions of high streamwise velocity (red shading) and regions of upwelling (
$v^\prime$ away from wall) coincide with regions of low streamwise velocity (blue shading). These observations are consistent with known features of the NW cycle (Robinson Reference Robinson1991; Waleffe Reference Waleffe1997; Jiménez & Pinelli Reference Jiménez and Pinelli1999). When the solid wall is replaced by anisotropic permeable substrates with 
$\phi _{xy} = 5.5$, there are some subtle but important changes in mode structure. For the permeable substrates that yield mode suppression – and drag reduction – there is a downward shift in the velocity fluctuations and the mode is compressed in the wall-normal direction (see figure 4c,d for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 0.18$ and figure 4e,f for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 0.45$). To explain why these changes lead to mode suppression, we briefly review two mechanisms responsible for high gain in the resolvent framework (McKeon & Sharma Reference McKeon and Sharma2010; McKeon Reference McKeon2017). First, resolvent modes tend to localize around the critical layer, 
$y_c$, where the mode speed matches the local mean velocity, 
$\bar {U}(y_c) = c = \omega /\kappa _x$. This minimizes the term 
$-\textrm {i}\omega + \textrm {i}\kappa _x \bar {U}$ inside the inverted operator shown on the right-hand side of (2.4). Second, energy is extracted from the mean flow via the interaction between the wall-normal velocity fluctuations and the mean shear, i.e. the so-called lift-up mechanism that leads to energy transfer from the mean flow to the fluctuations via the term 
$-\overline {u^\prime v^\prime } (\textrm {d}\bar {U}/{\textrm {d} y})$ in the turbulent kinetic energy equation. As the streamwise permeability of the substrate increases, there is a shift in the virtual origin of the mean profile and so the critical layer location moves closer to the permeable interface, from 
$\hat {y}_c^+ \approx 13$ for the smooth-wall case to 
$\hat {y}_c^+ \approx 9.7$ for the case with 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 0.45$. Here, 
$\hat {y} = 1+y$ such that 
$\hat {y}=0$ represents the porous interface at 
$y=-1$ (see figure 1). As noted in the introduction, this shift in the mean profile depends on 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}}$. A downward shift in the critical layer location causes the NW resolvent mode to localize closer to the porous interface (see 
$\rightarrow -\leftarrow$ in figure 4). However, the downward shift for the NW mode is constrained by the low wall-normal and spanwise permeabilities, 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = \sqrt {\smash {K_z^+}\vphantom {K_x^+}}$, which leads to the observed wall-normal compression in mode structure. This wall-normal compression in mode structure reduces the energy extracted from the mean flow, which depends on the wall-normal integral of the product of the Reynolds shear stress and mean shear, i.e. 
$\propto \int Re(-u_{\boldsymbol {{\kappa }}}^* v_{\boldsymbol {{\kappa }}}) (\textrm {d}\bar {U}/{\textrm {d} y})\,{\textrm {d} y}$. Note that 
$Re(\cdot )$ represents the real component and so the Reynolds shear stress contribution from the NW resolvent mode is proportional to 
$Re(-u_{\boldsymbol {{\kappa }}}^* v_{\boldsymbol {{\kappa }}})$ (Luhar et al. Reference Luhar, Sharma and McKeon2014). A reduction in energy extraction from the mean flow leads to a reduction in mode gain. Importantly, this observation is consistent with the interpretation that drag reduction is only expected if the slip length for the mean streamwise flow 
$\ell _x^+ \propto \sqrt {\smash {K_x^+}\vphantom {K_x^+}}$ is larger than that for the turbulent cross-flows, 
$\ell _z^+ \propto \sqrt {\smash {K_z^+}\vphantom {K_x^+}}$, i.e. if the virtual origin for the mean flow is offset from the virtual origin for the turbulence. Additional resolvent-based predictions for 
$\sqrt {\smash {K_z^+}\vphantom {K_x^+}} > \sqrt {\smash {K_x^+}\vphantom {K_x^+}}$ (
$\ell _z^+ > \ell _x^+$) lead to an increase in mode gain (data not shown).
Figure 4. The NW mode structure predictions made using the DNS velocity profiles for (a,b) the smooth-wall case; (c,d) the substrate with 
$\phi _{xy} = 5.5$ and 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 0.18$; (e,f) the substrate with 
$\phi _{xy} = 5.5$ and 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 0.45$; (g,h) the substrate with 
$\phi _{xy} = 5.5$ and 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 1.0$. In panels (a,c,e,g), the shading shows normalized contours of positive (red) and negative (blue) streamwise velocity normalized by the maximum value. The vectors show the wall-normal and spanwise velocity fluctuations. In panels (b,d,f,h), the solid lines show the magnitude of the streamwise velocity for this resolvent mode, 
$|u_{\boldsymbol {{\kappa }}}|$ and the dashed lines show the wall-normal velocity magnitude multiplied by a factor of ten, 
$10 |v_{\boldsymbol {{\kappa }}}|$. The black lines represent the smooth-wall case while the grey lines represent the permeable substrates. These plots make use of the shifted coordinate 
$\hat {y}=1+y$, such that 
$\hat {y}=0$ represents the location of the smooth wall or the porous interface.
Figure 4(g,h) show predictions for NW mode structure for the anisotropic permeable substrate with 
$\phi _{xy} = 5.5$ and 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 1.0$ that leads to an increase in drag and forcing–response gain (see figure 3). Compared with the cases that lead to mode suppression and drag reduction, the wall-normal footprint for the NW mode increases substantially for this case. This can be attributed to the change in the DNS mean profiles for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$, for which spanwise-coherent structures resembling KH vortices emerge over the permeable substrate. The emergence of these rollers leads to a substantial increase in turbulent mixing which causes the mean profile to flatten in the interfacial region (Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). For the substrate with 
$\phi _{xy} = 5.5$ and 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 1.0$, this pushes the critical layer for the NW mode outwards to 
$\hat {y}_c^+ \approx 26$. This means that the NW mode is again able to localize around the critical layer, and there is no vertical compression in mode structure. This allows for greater energy extraction from the mean flow and leads to the observed increase in NW mode gain over the material with 
$\phi _{xy}= 5.5$ and 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 1.0$ in figure 3(a).
Note that the NW mode gain decreases monotonically with increasing 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_y^+}\vphantom {K_x^+}}$ for predictions made using the synthetic mean profiles (see figure 3c). This is because the synthetic mean profiles assume no change in turbulence characteristics near the interface relative to smooth-wall conditions, i.e. they do not account for the emergence of spanwise rollers and turbulence penetration into the porous medium. As shown in figure 3(a), NW mode gain begins to increase again for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$ when the DNS mean profiles are used in the resolvent operator. In the DNS, spanwise rollers resembling KH vortices emerge beyond this value of wall-normal permeability and alter the form of the mean profile.
The preceding discussion assumes that the phase speed of the near-wall structures (
$c^+ \approx 10$) remains unchanged relative to smooth wall conditions. It could also be argued that the phase speed changes over the permeable medium to reflect the interfacial slip velocity, which depends on the streamwise slip length (
$\bar {U}_s^+ \approx \ell _x^+$ for small 
$\ell _x^+$), and turbulence penetration into the permeable medium, which depends on the spanwise slip length (
$\ell _z^+$). To account for these effects, the phase speed can be modified to 
$c^+ = 10 + \ell _x^+ - \ell _z^+ \approx 10 + \sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_z^+}\vphantom {K_x^+}}$. Figure 8 in the Appendix shows additional predictions for NW mode gain made using this modified phase speed. These additional predictions are broadly similar to the predictions shown in figure 3 for 
$c^+ = 10$. Specific observations are discussed further in the Appendix.
Together, the predictions shown in figures 3 and 4 confirm that suppression of the dynamically important NW cycle requires high streamwise permeability and low spanwise/wall-normal permeabilities. In other words, materials with high anisotropy ratios (
$\phi _{xy} \gg 1$) are expected to yield the largest suppression of NW turbulence. However, the absolute value of the wall-normal permeability also plays an important role in dictating drag reduction performance. Specifically, the emergence of energetic spanwise coherent rollers is dictated primarily by 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$.
3.3. Spanwise-coherent resolvent modes
In this section, we use the resolvent analysis framework to evaluate the conditions in which energetic spanwise-coherent rollers resembling KH vortices emerge over anisotropic permeable substrates. Snapshots of the flow field from the simulations of Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) show that for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} < 0.4$, the flow is dominated by streamwise-elongated structures associated with the NW cycle. As the wall-normal permeability increases to 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$ the flow becomes dominated by spanwise coherent rollers with an approximate streamwise spacing of 
$\lambda _x^+ \approx 100-300$. Furthermore, using momentum balance arguments, Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) showed that the Reynolds stress contribution from structures with 
$\lambda _x^+ \approx 70 - 320$ and 
$\lambda _z^+ \geqslant 120$ is responsible for the increase in drag observed for substrates with 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$. We note that spanwise-coherent rollers have also been observed in previous simulations over isotropic porous materials (Breugem et al. Reference Breugem, Boersma and Uittenbogaar2006; Rosti et al. Reference Rosti, Cortelezzi and Quadrio2015).
Linear stability analyses successfully predict the emergence of KH vortices in channel flows over anisotropic permeable substrates (Abderrahaman-Elena & García-Mayoral Reference Abderrahaman-Elena and García-Mayoral2017; Gómez-de-Segura et al. Reference Gómez-de-Segura, Sharma and García-Mayoral2018). Gómez-de-Segura et al. (Reference Gómez-de-Segura, Sharma and García-Mayoral2018) also show that the wall-normal permeability becomes the driving parameter for the onset of KH rollers for streamwise-preferential substrates of large depth, i.e. substrates with 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} \gg \sqrt {\smash {K_y^+}\vphantom {K_x^+}}$ and 
$H^+ \gg \sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. However, models based on linear stability theory do not accurately predict the threshold value of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$ for the onset of the rollers or the wavelength of the rollers. Linear stability analyses predict structures with streamwise wavelength 
$\lambda _x^+ \approx 50\text {--}70$ to be most unstable, while DNS results suggest that structures with 
$\lambda _x^+ \approx 150$ are most energetic initially.
Figure 5 shows resolvent-based predictions for the normalized amplification of spanwise-coherent structures over the anisotropic permeable substrate with 
$\phi _{xy}= 5.5$ and varying 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. Based on the DNS results of Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019), we limit ourselves to the evaluation of structures with streamwise wavelength 
$\lambda _x^+\in [50,500]$ and mode speed 
$c^+\in [4,10]$. In addition, we consider structures that are infinitely long in the spanwise direction (
$\kappa _z = 0$) as well as structures with a finite, but large, spanwise wavelength (
$\kappa _z \approx 2.3$; 
$\lambda _z^+ = 500$).
Figure 5. Normalized gain for spanwise-coherent structures plotted as a function of streamwise length 
$\lambda _x^+$ and wave speed 
$c^+$ for substrates with an anisotropy 
$\phi _{xy}= 5.5$. These predictions make use of the synthetic mean profile. Panels (a–d) represent structures with spanwise wavenumber 
$\kappa _z\approx 2.3$ (
$\lambda _z^+=500$). Panels (e–h) represent structures with spanwise wavenumber 
$\kappa _z=0$. Red and blue shading represent mode amplification and suppression relative to the smooth-wall case, respectively. Wall-normal permeability increases from left to right. Structures with the largest amplification are labelled with a (
$\bullet$) marker.
Figure 5(a,e) shows that for low wall-normal permeability, 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}=0.18$, resolvent analysis predicts minimal amplification of spanwise-coherent structures over the anisotropic permeable substrate. The maximum increase in gain relative to smooth-wall values is less than 
$25\,\%$, i.e. 
$\sigma _{\boldsymbol {{\kappa }},p}/\sigma _{\boldsymbol {{\kappa }},s}<1.25$ for both 
$\kappa _z = 0$ and 
$\kappa _z \approx 2.3$. Moreover, the gain for these spanwise-coherent structures is also small in absolute terms. This is consistent with the simulation results of Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019), which show that the Reynolds stress contributions from structures with 
$\lambda _x^+\approx 70\text {--}300$ and 
$\lambda _z^+\gtrsim 120$ are negligible for materials with 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}\lesssim 0.31$. In other words, spanwise-coherent modes do not have a significant effect on drag for these conditions.
However, the amplification of these spanwise-coherent structures increases significantly for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \geqslant 0.39$. Figure 5(b,f) shows a clear increase in the size and intensity of the region showing mode amplification. For 
$\lambda _z^+=500$, the region of highly amplified structures is localized around 
$(\lambda _x^+,c^+)\approx (130,7)$ and the normalized gain for the most amplified structure is 
$\sigma _{\boldsymbol {{\kappa }},p}/\sigma _{\boldsymbol {{\kappa }},s} \approx 2$. For spanwise constant structures, the region of high amplification is localized around 
$(\lambda _x^+,c^+)\approx (180,6.7)$ and extends to structures with 
$\lambda _x^+ \approx 300$. The normalized gain for structures in this region is as high as 
$\sigma _{\boldsymbol {{\kappa }},p}/\sigma _{\boldsymbol {{\kappa }},s} \approx 3$, indicating a 200 % increase in gain relative to smooth-wall values. The extent and intensity of these high-gain regions continue to increase as wall-normal permeability increases further, as shown in figures 5(c,d) and 5(g,h). Indeed, for the case with 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 0.66$ shown in figure 5(d), the region of high amplification extends from 
$\lambda _x^+ \approx 100\text {--}500$ and 
$c^+ \approx 4\text {--}9$. The highest-amplification structures in this region show normalized gains 
$\sigma _{\boldsymbol {{\kappa }},p}/\sigma _{\boldsymbol {{\kappa }},s} \geqslant 100$.
The substantial increase in amplification of spanwise-coherent structures for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \geqslant 0.39$ is again consistent with the DNS results of Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019), which show that the drag-reducing performance of anisotropic substrates begins to degrade for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \geqslant 0.39$ due to the additional Reynolds shear stresses contributed by spanwise coherent structures with 
$\lambda _x^+ \approx 70\text {--}300$ and 
$\lambda _z^+ > 120$. The DNS results also show the emergence of a new peak in the premultiplied spectra for wall-normal velocity near 
$\lambda _x^+ \approx 150$ for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} = 0.39$. This is very close to the streamwise wavelength of the most amplified modes in figure 5(b,f). Moreover, the streamwise extent of the high-intensity region in the DNS premultiplied spectra increases with increasing 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$, which is similar to the expansion of the high-amplification region evident in figure 5(c,d,g,h).
Figure 6 shows a snapshot of the flow field associated with the highest-gain resolvent mode in figure 5(f) with 
$\lambda _x^+ \approx 180$ and 
$c^+ \approx 6.7$. As expected for KH type vortices, the structure shows the presence of counter-rotating rollers in the 
$x$–
$y$ plane, localized near the porous interface. Regions of prograde rotation (i.e. in the direction of the mean shear) are associated with negative pressure fluctuations and regions of retrograde rotation are associated with positive pressure fluctuations. Moreover, the velocity and pressure fields associated with these rollers span a significant portion of the buffer region of the flow, which is consistent with DNS observations.
Figure 6. Flow structure associated with the most amplified spanwise-constant resolvent mode over the permeable substrate with 
$\phi _{xy}=5.5$ and 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}=0.39$; see 
$\bullet$ in figure 5(f). The red and blue shaded contours show regions of positive and negative pressure, normalized by the maximum value.
Importantly, the model predictions shown in figure 5 for 
$\phi _{xy} = 5.5$ are representative of what is observed over the permeable substrates with 
$\phi _{xy} = 3.6$ and 
$\phi _{xy} = 11.4$ as well. This is illustrated well in figure 7, which shows the normalized gain for the most amplified resolvent mode in the window spanning 
$\lambda _z^+\in [50,500]$ and 
$c^+\in [4,10]$, for both 
$\kappa _z \approx 2.3$ and 
$\kappa _z = 0$. For all three anisotropy ratios, the normalized gain show a sharp increase for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$. For structures with 
$\lambda _z^+ = 500$ (
$\kappa _z \approx 2.3$), model predictions for all three anisotropy ratios collapse together nicely (figure 7a). For spanwise-constant structures (
$\kappa _z = 0$), there is a bit more scatter in the normalized gain values but the overall trend remains similar.
Figure 7. Comparison of resolvent-based gain predictions for spanwise-coherent structures with anisotropy ratios 
$\phi _{xy}=3.6$ (
$\square$), 
$5.5$ (
$\triangle$) and 
$11.4$ (
). The maximum normalized gain obtained for resolvent modes with 
$\lambda _x^+\in [50,500]$ and 
$c^+ \in [4,10]$ is shown as a function of wall-normal permeability for (a) 
$\kappa _z \approx 2.3$ (
$\lambda _z^+ = 500$) and for (b) 
$\kappa _z = 0$ (
$\lambda _z^+ = \infty$). All of these predictions were obtained using the synthetic mean profiles.
The DNS results in figure 3(a) clearly show drag reduction performance deteriorating for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$ for all three anisotropy ratios. The flow field snapshots, velocity spectra, and momentum-based analyses pursed by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) attribute this deterioration in performance to the emergence of high-gain spanwise-coherent structures resembling KH vortices. The results presented in this section confirm that resolvent analysis – based on synthetic mean profiles – is able to successfully predict the emergence of such spanwise-coherent structures over anisotropic permeable substrates. Since the synthetic mean profiles are generated using an eddy viscosity model, these results also show that the resolvent framework can generate a priori predictions for whether a porous substrate with known 
$\boldsymbol{\mathsf{K}}$ will give rise to KH rollers. Moreover, the resolvent-based predictions are in better quantitative agreement with the DNS results than linear stability predictions. Specifically, resolvent analysis yields better predictions for the wall-normal permeability threshold beyond which such structures become important (
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \approx 0.4$), as well as the streamwise wavelength of the structures that emerge first (
$\lambda _x^+ \approx 150$).
4. Conclusion
Recent theoretical efforts and numerical simulations show that anisotropic permeable substrates have the potential to yield significant drag reduction in wall-bounded turbulent flows. Specifically, numerical simulations by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) indicate that a drag reduction of up to 25 % can be achieved by appropriately tuning the streamwise, wall-normal and spanwise permeabilities of such substrates. The mechanism through which these substrates reduce drag is similar to that for riblets, and can be rationalized in terms of the offset in the virtual origin felt by the mean flow, which depends on streamwise permeability, and the virtual origin felt by the turbulent fluctuations, which depends primarily on spanwise permeability. Maximum drag reduction is limited by the appearance of spanwise rollers resembling KH vortices. The emergence of these vortices is linked to a relaxation in the wall-normal permeability. The extended resolvent analysis framework described in this paper reproduces all of these features with minimal computation.
As shown in § 3.2, the gain for a single resolvent mode that serves as a surrogate for the NW cycle reproduces the initial drag reduction trends observed in numerical simulations. Model predictions show that the reduction in gain for this mode depends on the difference between the streamwise and spanwise permeability length scales, 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_z^+}\vphantom {K_x^+}}$, such that substrates with higher anisotropy ratios lead to greater mode suppression. This is consistent with the trends observed in numerical simulations. Moreover, the link between resolvent mode gain and drag reduction also provides a complementary explanation to the slip-length based arguments used in previous studies. Specifically, resolvent-based predictions suggest that the difference between 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}}$ and 
$\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$ leads to a wall-normal compression of the quasi-streamwise vortices associated with the NW cycle, and this compression limits energy extraction from the mean flow. Furthermore, as we show in § 3.3, resolvent analysis also predicts the emergence of energetic spanwise rollers with streamwise wavelength 
$\lambda _x^+ \approx 150$ as the wall-normal permeability increases beyond 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \approx 0.4$. These features are quantitatively consistent with simulation results.
One weakness of the extended resolvent framework developed here is the requirement of a mean velocity profile over the permeable substrate. For turbulent flows over smooth walls, mean profile estimates can be generated using a variety of simplified models or obtained from previous simulations and experiments. However, such models or data are not readily available yet for flows over anisotropic permeable substrates. Since the form of the mean profile can have a significant effect on resolvent-based predictions, it would be interesting to revisit the predictions made in this study as more accurate eddy viscosity models become available for flows over (and within) anisotropic porous materials. Indeed, the comparison between predictions made using the mean profiles obtained in DNS and those generated using an eddy viscosity model presented in § 3.2 shows important differences in the gain for the near-wall mode for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.4$, i.e. when drag reduction performance begins to deteriorate. This is because the emergence of the spanwise rollers leads to a substantial increase in interfacial turbulence. This effect is not captured by the eddy viscosity model used to generate the synthetic profiles, which assumes that the turbulence remains smooth-wall-like and does not penetrate the permeable substrate (see § 2.4). Yet, resolvent analysis with the synthetic mean profile does show the emergence of these spanwise rollers beyond 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \approx 0.4$. In other words, resolvent analysis can generate useful a priori predictions even with a synthetic mean profile. Near-wall mode gain can be used as a measure of initial drag reduction performance, and the emergence of spanwise-coherent rollers can be used to estimate the point at which performance is likely to deteriorate. Since evaluation of these features does not require significant computation, resolvent analysis can serve as a useful design tool for more complex permeable substrates. For instance, resolvent analysis could be used to pursue formal optimization efforts for the full permeability tensor.
We also recognize that the VANS equations used in this study involved significant simplifying assumptions (see § 2.1). These simplifying assumptions were made primarily to allow for a direct comparison with previous simulation results. Moving forward, resolvent analysis could be used to consider how the inclusion of inertial effects (e.g. via the Forchheimer term) or more complex interfacial boundary conditions (e.g. involving stress jumps, or more complex slip and transpiration models) is likely to affect control performance. Recent work by Morra et al. (Reference Morra, Semeraro, Henningson and Cossu2019) also shows that the statistical predictions made using resolvent analysis improve significantly when the effect of the Reynolds stresses is included explicitly in the resolvent operator via an eddy viscosity term, i.e. if the linear operator in (2.5) is formulated using both molecular and eddy viscosities. In the context of the present study, the inclusion of an appropriate eddy viscosity model in the resolvent operator itself would enable an evaluation of how turbulent momentum and energy exchanges across the permeable interface are likely to modify NW cycle and influence the emergence of spanwise rollers.
Funding
This material is based on work supported by the Air Force Office of Scientific Research under awards FA9550-17-1-0142 (program manager Dr G. Abate) and FA9550-19-1-7027 (program manager: Dr D. Smith) and the National Science Foundation under grant no. 1943105.
Declaration of interests
The authors report no conflict of interest.
Appendix
As noted earlier, the predictions shown in § 3.2 assume that the phase speed of the resolvent mode that serves as a surrogate for the NW cycle (
$c^+ \approx 10$) remains unchanged over the permeable substrate. Figure 8 shows additional predictions for mode gain assuming the phase speed change with substrate permeability as 
$c^+ = 10 + \ell _x^+ - \ell _z^+ \approx 10 + \sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_z^+}\vphantom {K_x^+}}$. In effect, this modification accounts for the interfacial slip velocity, which depends on the streamwise slip length 
$\ell _x^+ \approx \sqrt {\smash {K_x^+}\vphantom {K_x^+}}$, and the virtual origin for the turbulence inside the permeable substrate, which depends on the spanwise slip length 
$\ell _z^+ \approx \sqrt {\smash {K_z^+}\vphantom {K_x^+}}$.
Figure 8. Predictions for the NW resolvent mode with a modified phase speed of 
$c^+ \approx 10 + \sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_z^+}\vphantom {K_x^+}}$. (a) Comparison between normalized mode gain (
$\sigma _{\boldsymbol {{\kappa }},p}/\sigma _{\boldsymbol {{\kappa }},s}$, black symbols) and drag reduction observed in DNS (
$\Delta \bar {U}^+$, light grey symbols) by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019), plotted as a function of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. The resolvent-based predictions shown in this panel are obtained using the velocity profiles from DNS. (b) Comparison between normalized gains predicted using the DNS mean profiles (black symbols) and the synthetic mean profiles (white symbols) plotted as a function of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. (c) Comparison between normalized gains predicted using the DNS mean profiles (black symbols) and the synthetic mean profiles (white symbols) plotted as a function of 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}} - \sqrt {\smash {K_y^+}\vphantom {K_x^+}}$. The black and red dashed lines show linear fits to the initial decrease in normalized gain for the DNS mean and synthetic mean predictions, respectively. For all panels, the 
$\square$ symbols represent substrates with 
$\phi _{xy} = 3.6$, 
$\triangle$ symbols represent substrates with 
$\phi _{xy}= 5.5$ and 
symbols represent substrates with 
$\phi _{xy}= 11.4$.
The modified phase speed does not lead to significant changes in mode gain compared with the predictions shown in figure 3 for 
$c^+ \approx 10$. Figure 8(a) shows that mode gain predictions made using the DNS mean profile are consistent with drag reduction trends obtained in DNS. Compared with the predictions shown in figure 3(a), singular value ratios for the substrates with 
$\phi _{xy}=3.6$ and 
$\phi _{xy}=5.5$ show slightly better agreement with the drag reduction trends from DNS for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} < 0.6$. However, for 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \gtrsim 0.6$, the predictions are more scattered with the modified phase speed. For the substrate with 
$\phi _{xy}= 11.4$, the modified phase speed leads to reduced mode suppression. Furthermore, the deterioration in performance – in terms of mode suppression – commences at a lower value of wall-normal permeability, 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}} \approx 0.2$, compared with the threshold observed in DNS.
Figure 8(b) shows that the predictions made using the synthetic mean profile are consistent with those made using the DNS mean profiles for small values of 
$\sqrt {\smash {K_y^+}\vphantom {K_x^+}}$, i.e. until the emergence of spanwise rollers resembling KH vortices in the DNS. Figure 8(c) shows that the initial reduction in mode gain depends primarily on 
$\sqrt {\smash {K_x^+}\vphantom {K_x^+}}-\sqrt {\smash {K_z^+}\vphantom {K_x^+}}$, even with the modified phase speed. These observations are again consistent with the predictions shown in figure 3(b,c). However, linear fits to the initial decrease in mode gain show that the modified phase speed leads to less pronounced mode suppression; see red and black dashed lines in figures 3(c) and 8(c).
Open access
