Resolvent-based predictions for turbulent flow over anisotropic permeable substrates

Recent simulations indicate that streamwise-preferential porous materials have the potential to reduce drag in wall-bounded turbulent flows(Gomez-de-Segura&Garcia-Mayoral 2019). This paper extends the resolvent formulation to study the effect of such anisotropic permeable substrates on turbulent channel flow. Under the resolvent formulation, the Fourier-transformed Navier-Stokes equations are interpreted as a linear forcing-response system. The nonlinear terms are considered the endogenous forcing in the system that gives rise to a velocity and pressure response. A gain-based decomposition of the forcing-response transfer function---the resolvent operator---identifies response modes (resolvent modes) that are known to reproduce important structural and statistical features of wall-bounded turbulent flows. The effect of permeable substrates is introduced in this framework using the Volume-Averaged Navier-Stokes equations and a generalized form of Darcy's law. Substrates with high streamwise permeability and low spanwise permeability are found to suppress the forcing-response gain for the resolvent mode that serves as a surrogate for the energetic near-wall cycle. This reduction in mode gain is shown to be consistent with the drag reduction trends predicted by theory and observed in numerical simulations. Simulation results indicate that drag reduction is limited by the emergence of spanwise rollers resembling Kelvin-Helmholtz vortices beyond a threshold value of wall-normal permeability. The resolvent framework also predicts the conditions in which such energetic spanwise-coherent rollers emerge. These findings suggest that a limited set of resolvent modes can serve as the building blocks for computationally-efficient models that enable the design and optimization of permeable substrates for passive turbulence control.


Introduction
Many flows of engineering and scientific interest involve permeable substrates. The presence of such complex substrates can significantly alter the behavior 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 1991;Waleffe 1997;Jiménez & Pinelli 1999) and reduce drag in wall-bounded turbulent flows (Hahn et al. 2002;Itoh et al. 2006;Abderrahaman-Elena & García-Mayoral 2017;Gómez-de-Segura et al. 2018;Rosti et al. 2018;Gómez-de-Segura & García-Mayoral 2019). 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 2007). 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. 2002;Breugem et al. 2006;Manes et al. 2011;Zampogna & Bottaro 2016;Efstathiou & Luhar 2018;Rosti et al. 2018;Gómez-de-Segura & García-Mayoral 2019;Kim et al. 2020). 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 1999;Breugem et al. 2006;Efstathiou & Luhar 2018). 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 2019). 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 2010;McKeon 2017) 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.

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 2019). The physical mechanism underlying this drag reduction is similar to the mechanism that yields drag reduction over riblets (Walsh & Lindemann 1984;Luchini et al. 1991;Luchini 1996;Bechert et al. 1997;García-Mayoral & Jiménez 2011a). 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 to the spanwise cross-flows arising from turbulent fluctuations (Luchini et al. 1991). 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 1996;Gómezde-Segura & García-Mayoral 2019;García-Mayoral et al. 2019). 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 (2017) 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 ( + x , + 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 τ , and viscosity, ν. The analysis of Abderrahaman-Elena & García-Mayoral (2017) showed that + x ∝ K + x and + z ∝ K + z if the height of the substrate, H, is much larger than the permeability length scales K + x and K + z . The achievable drag reduction was estimated to be proportional to the difference between the streamwise and spanwise slip lengths, ∆D ∝ + x − + z , or equivalently, ∆D ∝ K + x − K + z . This is consistent with the findings of Busse & Sandham (2012), who studied the effect of anisotropic slip length boundary conditions in turbulent channel flow simulations.
Recent direct numerical simulation (DNS) results obtained by Gómez-de-Segura & García-Mayoral (2019) provide further support for the scaling developed by Abderrahaman-Elena & García-Mayoral (2017). 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. The permeable substrate was characterized by a permeability tensor of the form K = 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 φ xy = K + x / K + y ranged between 3.6 and 11.4. As predicted by Abderrahaman-Elena & García-Mayoral (2017), these simulations show that for surfaces for which 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, ∆D ∝ K + x − K + z . 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. 2006;Efstathiou & Luhar 2018), and they also contribute to performance degradation for riblets (García-Mayoral & Jiménez 2011b. 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 2017;Gómez-de-Segura et al. 2018;Gómez-de-Segura & García-Mayoral 2019). Specifically, simulation results show that spanwise rollers become increasingly energetic as the wall-normal permeability exceeds K + y ≈ 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 near-wall 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, ∆D ∝ K + x − K + z . However, it is unclear if this relationship also holds for more complex surfaces that are not characterized by diagonal permeability tensors of the form K = 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 2017;Lācis et al. 2020;Bottaro 2019). The effect of such slip-length models on the nearwall 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 K + y beyond which KH rollers become energetic (Abderrahaman-Elena & García-Mayoral 2017; Gómez-de-Segura et al. 2018). Further, the streamwise wavelengths predicted to be most unstable do not match the length scale of the spanwise rollers observed in simulations.

Contribution and Outline
In this study, we develop a reduced-order modeling framework grounded in resolvent analysis (McKeon & Sharma 2010;McKeon 2017) 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 gainbased 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 dynamicallyimportant flow features such as the NW cycle (Moarref et al. 2013). 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 et al. 2014;Nakashima et al. 2017;Toedtli et al. 2019) and passive (Luhar et al. 2015;Chavarin & Luhar 2020) control techniques in wall-bounded turbulent flows. In particular, recent work by Chavarin & Luhar (2020) 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 (2020) 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 2011b). Motivated by these prior modeling 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 (2019), we consider a symmetric channel geometry at friction Reynolds number Re τ = 180 and substrates with identical wall-normal and spanwise permeabilities, i.e., substrates with φ yz = K + y / K + z = 1. We model the flow in the substrate using the volume-averaged Navier-Stokes equations, in which the effect of the permeable substrate is included via a permeability tensor. However, this modeling framework can be extended to include more sophisticated interfacial boundary conditions (Lācis & Bagheri 2017;Lācis et al. 2020), and to account for inertial effects via the so-called Forchheimer term (Breugem et al. 2006;Whitaker 1999).
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 spanwisecoherent resolvent modes resembling KH rollers. These predictions are compared against DNS results from Gómez-de-Segura & García-Mayoral (2019). 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 against predictions generated using the mean velocity profiles obtained in DNS by Gómez-de-Segura & García-Mayoral (2019). Section 4 concludes the paper.

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 conditiobluns 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).

Accounting for Permeable Substrates
The resolvent framework is reformulated using the volume-averaged Navier-Stokes equations. Volume-averaging gives rise to two additional terms: a term representing the sub-filter 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 1969(Whitaker , 1996(Whitaker , 1999. 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 1996). This model has been used in previous numerical simulations of flow over porous substrates (Breugem et al. 2006;Rosti et al. 2015Rosti et al. , 2018 as well as in linear stability analyses (Tilton & Cortelezzi 2006, 2008. The volume-averaged Navier-Stokes equations and continuity constraint for flow through a porous medium with porosity ε, dimensionless permeability K , and dimensionless Forchheimer resistance F can be expressed as: and ∇ · (ε u ) = 0. (2.1b) In the equations above, · represents a volume-averaged quantity, u is the dimensionless velocity, p 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 τ . The friction Reynolds number is given by Re τ = u τ h/ν and the dimensionless permeability defined as K = K † /h 2 , where K † is the dimensional permeability. The quantity τ = uu − u u represents the sub-filter scale stresses which arise from volume averaging the Navier-Stokes equations. For our analysis we consider the following simplifications. Consistent with prior numerical simulations (Rosti et al. 2018; Gómez-de-Segura & García-Mayoral 2019), we focus on substrates characterized by a permeability tensor of the form K = 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, wallnormal, and spanwise directions of the flow. Second, we omit the nonlinear Forchheimer correction term, F, that is used to account for inertial effects in flows through porous media (Ochoa-Tapia & Whitaker 1995a,b;Breugem et al. 2006;Whitaker 1999). This assumption is again consistent with the numerical simulations we will use to test model predictions. Third, we assume that the porosity of permeable substrates is ε ≈ 1 to maximize any potential drag reduction (Abderrahaman-Elena & García-Mayoral 2017). 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 sub-filter scale stresses (Breugem et al. 2006). With these assumptions, (2.1a) and (2.1b) can be expressed as: and ∇ · u = 0, (2.2b) where the · 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 ∈ [−1, 1]. The regions occupied by the permeable substrates correspond to y ∈ [−(1 + H), −1) and y ∈ (1, 1 + H]. The height of the permeable substrate is H. Note that all lengths are normalized by the channel halfheight.

Modified Resolvent Analysis
The resolvent formulation for wall-bounded turbulent flows proposed by McKeon & Sharma (2010)-and employed in several recent flow control studies (Luhar et al. 2014(Luhar et al. , 2015Nakashima et al. 2017;Toedtli et al. 2019;Chavarin & Luhar 2020)-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 (2017).
Construction of the modified resolvent operators begins with a standard Reynoldsdecomposition of the simplified volume-averaged Navier-Stokes equations in (2.2). The velocity and pressure fields are decomposed into a time-averaged component (denoted by (·)) and a fluctuating component about this average (denoted by (·) ). Under this decomposition, the velocity field is expressed as u(t, x) = U (x) + u (t, x) and the pressure field is expressed as p(t, x) = P (x) + p (t, x). Note that U (x) = [U (y), 0, 0] 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: In the expression above, κ x is the streamwise wavenumber, κ z is the spanwise wavenumber, and ω is the frequency. The Fourier coefficients for the velocity and pressure field at a given wavenumber-frequency combination, κ = (κ x , κ z , ω), are denoted u κ and p κ .
Under the Fourier transform, the equations in (2.2) can be expressed compactly as: In the expression above, the first row represents the momentum equations and the second row represents the continuity constraint. The operator (2.5) represents the linear dynamics in (2.2), and f κ is the Fourier coefficient for the nonlinear terms. The differential operator∇ is defined as∇ = (iκ x , ∂/∂y , iκ z ), and so the symbols∇ T and∇ essentially represent the Fourier-transformed divergence and gradient operators, respectively. Similarly, the Laplacian is defined as∇ 2 = (−κ 2 x − κ 2 z + ∂ 2 ∂y 2 ). Note that the velocity and pressure response at a given wavenumberfrequency combination constitute a traveling wave flow field with streamwise wavelength λ x = 2π/κ x and spanwise wavelength λ z = 2π/κ z that is moving downstream at speed c = ω/κ x . The transfer function that maps the nonlinear forcing f κ to the velocity and pressure response [u κ , p κ ] T in (2.4) is the resolvent operator, H κ . 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 K in (2.5). Note that the resistance exerted by the Darcy term is linear with respect to u.
As detailed in prior studies (McKeon & Sharma 2010;Moarref et al. 2013;Luhar et al. 2014Luhar et al. , 2015, a singular value decomposition (SVD) of the discretized resolvent operator (see Section 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: Here, W u and W u incorporate numerical quadrature weights for the entire domain spanning y ∈ [−(1 + H), (1 + H)]. With this weighting, the SVD of the scaled resolvent: yields forcing modes f κ,m = W −1 f φ κ,m and velocity responses u κ,m = W −1 u ψ κ,m that have unit energy when integrated across the entire domain spanning y ∈ [−(1 + H), (1 + H)]. In other words, this scaling ensures that In (2.7)-(2.8), the superscript (·) * denotes a conjugate transpose.
A major contribution of the resolvent framework lies in the finding that the forcingresponse transfer function tends to low-rank at κ combinations that are energetic in wall-bounded turbulent flows. Often, the first singular value is an order of magnitude larger than subsequent singular values, σ κ,1 σ κ,2 > ..., and so the resolvent operator can be well approximated using a rank-1 truncation after the SVD (McKeon & Sharma 2010; Moarref et al. 2013): The expressions in (2.6)-(2.9) show that forcing in the direction of the first forcing mode f κ,1 = W −1 f φ κ,1 generates a velocity response u κ,1 = W −1 u ψ κ,1 that is amplified by factor σ κ,1 . Put another way, a forcing of the form f κ,1 to the unscaled resolvent operator in (2.4) generates a velocity and pressure response σ κ,1 [u κ,1 , p κ,1 ] T . Under the L 2 scaling used here, σ 2 κ,1 is a measure of energy amplification. Recent modeling 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 λ + x ≈ 10 3 , λ + z ≈ 10 2 , and c + ≈ 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. 2014;Nakashima et al. 2017;Chavarin & Luhar 2020). 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 κ 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.

Boundary and Interface Conditions
As discussed in Section 2.5 below, the resolvent operator is discretized using spectral discretization and rectangular block matrices as described by Aurentz & Trefethen (2017). 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 ∈ [−1, 1] and the upper and lower permeable regions correspond to y ∈ (1, 1 + H] and y ∈ [−(1 + H), −1), respectively. At the lower and upper substrate walls, y = ±(1 + H), we apply no-slip boundary conditions. At the interfaces between the porous medium and the unobstructed flow, y = ±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: 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ν = ν/ε. Per the discussion presented in Abderrahaman-Elena & García-Mayoral (2017) and Gómez-de-Segura & García-Mayoral (2019), these assumptions are reasonable for a highly-connected medium with high porosity, ε ≈ 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 1995a;Minale 2014). Also keep in mind that the equations above imply a sharp transition between the porous medium and the unobstructed fluid. Previous studies employing the volume-averaged Navier-Stokes equations have typically assumed the existence of a finite transition zone between the homogeneous porous and homogeneous fluid regions (Ochoa-Tapia & Whitaker 1995a;Breugem et al. 2006). Finally, since we are interested in permeable substrates that have the highest potential for drag reduction, we assume a porosity of ε = 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 2019). Emulation of the channel configuration and boundary conditions used by Gómez-de-Segura & García-Mayoral (2019) allows for a direct comparison between model predictions and simulation results.

Mean Velocity Profile
As shown in (2.5), construction of the resolvent operator requires knowledge of the mean velocity profile 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 (2019), and (ii) mean profiles generated using a synthetic eddy viscosity profile. The DNS mean profile dataset consists of 22 different cases: 8 profiles for φ xy = 3.6; 7 profiles for φ xy = 5.5; and 7 profiles for φ xy = 11.4. For all 22 cases tested in DNS, the ratio between the wall-normal and spanwise permeability length scales was φ yz = K + y / K + z = 1. Each of the DNS profiles contained 153 points spanning the unobstructed region, y ∈ [−1, 1]. Within the permeable medium, the analytical solutions developed by Gómez-de-Segura & García-Mayoral (2019) 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: 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 (1967): (2.12) Here, the eddy viscosity has been normalized by the fluid viscosity ν. To generate the synthetic profiles, we used the values c 1 = 46.2 and c 2 = 0.61, which were identified by Moarref & Jovanović (2012) as yielding the best fit to the mean profiles obtained in smooth wall DNS at Re τ ≈ 180. Thus, the Reynolds shear stress term was modeled 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 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, equations (2.11) and (2.12) were combined to yield the following equation for the mean velocity profile: which was solved numerically to yield 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.

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 (2017). Each differential operator is discretized using Chebyshev polynomials and the resulting matrices are rectangular. The size of these matrices is [N ×N +n] and 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 (2017).
As shown in figure 1, for our problem configuration the channel is separated into 3 regions: the lower permeable substrate spanning y ∈ [−H, −1), the free channel spanning y ∈ [−1, 1], and the upper permeable substrate spanning y ∈ (1, H]. Although the configuration is symmetric across the channel centerline, we consider the entire channel to retain resolvent modes that are both symmetric and anti-symmetric across the channel centerline (Moarref et al. 2013;Luhar et al. 2015). 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.7s for N = 80, approximately 2.0s for N = 112 and 7.0s for N = 200. The results presented below were generated using N = 112, which corresponds to a total of 3N = 336 grid points across the the entire channel.

Results and Discussion
In this section, we compare resolvent-based predictions for the NW mode (Section 3.2) and the emergence of energetic spanwise rollers (Section 3.3) against DNS observations. Before that, we briefly compare the mean velocity profiles obtained from DNS against the synthetic profiles generated using the eddy viscosity model (Section 3.1). Figure 2 compares the interfacial slip velocity U + s predicted using the synthetic eddy viscosity profile (2.12) against those obtained in DNS by Gómez-de-Segura & García-Mayoral (2019) for each of the 22 different cases being considered here, with anisotropy ratios φ xy = 3.6, 5.5, and 11.4. The slip velocities for the synthetic profiles all collapse together onto a straight line corresponding to U + s ≈ K + x . In other words, the slip velocity only depends on the streamwise component of permeability, which is consistent with the assumptions outlined in Section 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.

Mean Velocity Profiles
The synthetic slip velocities (open symbols) show close agreement with DNS results (closed symbols) for low streamwise permeabilities but begin to deviate from DNS results at higher K + x . Moreover, the threshold values of K + x above which the synthetic predictions begin to deviate from DNS results depend on the anisotropy ratio, φ xy . As an example, for φ xy = 3.6, synthetic U + s begin deviating from DNS results for K + x 2. For φ xy = 11.4, the predicted slip velocities deviate significantly from DNS results for K + x 6. 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, K + y 0.4. This value of K + y corresponds closely to the conditions in which Kelvin-Helmholtz-type rollers appear over the permeable substrates in the simulations of Gómez-de-Segura & García-Mayoral (2019). 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 x − K + y . The black and red dashed lines show linear fits to the initial decrease in normalized gain for the DNS-and synthetic-mean predictions, respectively. For all panels, the symbols represent substrates with φxy = 3.6, symbols represent substrates with φxy = 5.5, and symbols represent substrates with φxy = 11.4. unobstructed region remains similar to that over a smooth wall regardless of the value of K + y and K + z . 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 Section 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 K + y 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.

Near-Wall Resolvent Mode
Previous studies have shown that the resolvent mode with streamwise wavelength λ + x = 10 3 , spanwise wavelength λ + z = 10 2 , and phase speed c + = (ω/κ 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 2010;. 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. 2014;Nakashima et al. 2017;Toedtli et al. 2019) and passive (Luhar et al. 2015(Luhar et al. , 2016Chavarin & Luhar 2020) 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 (2019) carried out at Re τ = 180, the resolvent mode with streamwise wavenumber κ x = 2πRe τ /10 3 ≈ 1.1, spanwise wavenumber κ z = 2πRe τ /10 2 ≈ 11 and frequency ω = 10κ x ≈ 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 σ κ,p normalized by the smooth-wall value σ κ,s at Re τ = 180 for all 22 substrates tested by Gómez-de-Segura & García-Mayoral (2019). Resolvent-based gain predictions (filled black symbols) are compared against the drag reduction measured in DNS (filled gray symbols). Following Gómez-de-Segura & García-Mayoral (2019), the outward shift in the logarithmic region of the mean profile is used to quantify the change in drag: ∆U + > 0 denotes a decrease in drag and ∆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 K + y , 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. Further, mode suppression and drag reduction both increase monotonically up to a value of K + y ≈ 0.4. Agreement between drag reduction and mode suppression is particularly good for substrates with anisotropy ratios φ xy = 3.6 and 5.5. For K + y 0.4, the normalized gain begins to increase again for the substrates with φ xy = 3.6 and 5.5, which is consistent with DNS drag reduction trend. DNS results show an increase in drag relative to smooth-wall values for K + y 0.55 for the substrates with φ 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 K + y 0.7 for these substrates. For substrates with the higher anisotropy ratio, φ 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, K + y . For instance, the increase in normalized gain beyond K + y 0.4 is less pronounced for the material with φ xy = 11.4. For this substrate, DNS results show an increase in drag for K + y 0.7. However, the resolvent predictions do not show an increase in gain beyond the smooth-wall value for the conditions tested (up to K + y ≈ 1).
Note that the observed disagreement between NW mode gain and drag reduction for , 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 K + y ≈ 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 Section 3.3. Figure 3(b) shows normalized gain predictions made using the DNS mean profile (filled symbols) and the synthetic mean profile (open symbols) for the NW resolvent mode, plotted as a function of K + y . For values of K + y 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 K + y 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 K + y is consistent with the slip velocities shown in figure 2. The synthetic profile agrees with the DNS mean profile for K + y 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 K + y and K + z . Figure 3(c) shows normalized gain predictions made using the DNS mean profile (filled symbols) and the synthetic mean profile (open symbols) for the NW resolvent mode, plotted as a function of K + x − K + y . 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., K + x − K + y = K + x − K + z . Slip length arguments (Luchini et al. 1991;Abderrahaman-Elena & García-Mayoral 2017) 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 crossflow fluctuations, respectively. In other words, we expect ∆U + ∝ K + x − K + z . These arguments are supported by the simulation results of Gómez-de-Segura & García-Mayoral (2019), which show a linear decrease in drag that is proportional to K + x − K + z 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 K + x − K + y (= K + x − K + z ) 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 behavior is selfsimilar when appropriately normalized. Predictions made using the DNS profiles diverge from this self similarity when K + y 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 λ + z = 10 2 . The wall-normal and spanwise velocity fluctuations show the presence of counter-rotating streamwise vortices. Regions of downwelling (v towards wall) coincide with regions of high streamwise velocity (red shading) and regions of upwelling (v away from wall) coincide with regions of low streamwise velocity (blue shading). These observations are consistent with known features of the NW cycle (Robinson 1991;Waleffe 1997;Jiménez & Pinelli 1999). When the solid wall is replaced by anisotropic permeable substrates with φ 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 Figs. 4(c,d) for K + y = 0.18 and Figs. 4(e,f) for K + y = 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 2010;McKeon 2017). First, resolvent modes tend to localize around the critical layer, y c , where the mode speed matches the local mean velocity, U (y c ) = c = ω/κ x . This minimizes the term −iω + iκ x 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 −u v (dU /dy) 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 . NW mode structure predictions made using the DNS velocity profiles for (a,b) the smooth-wall case, (c,d) the substrate with φxy = 5.5 and K + y = 0.18, and (e,f) the substrate with φxy = 5.5 and K + y = 0.45, and (g,h) the substrate with φxy = 5.5 and K + y = 1.0. In (a,c,e,f), 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 (b,d,f,g), the solid lines show the magnitude of the streamwise velocity for this resolvent mode, |uκ| and the dashed lines show the wall-normal velocity magnitude multiplied by a factor of ten, 10|vκ|. The black lines represent the smooth-wall case while the gray lines represent the permeable substrates. These plots make use of the shifted coordinatê y = 1 + y, such thatŷ = 0 represents the location of the smooth wall or the porous interface.
permeable interface, fromŷ + c ≈ 13 for the smooth-wall case toŷ + c ≈ 9.7 for the case with K + y = 0.45. Here,ŷ = 1 + y such thatŷ = 0 represents the porous interface at y = −1 (see figure 1). As noted in the introduction, this shift in the mean profile depends on K + x . A downward shift in the critical layer location causes the NW resolvent mode to localize closer to the porous interface (see → − ← in figure 4). However, the downward shift for the NW mode is constrained by the low wall-normal and spanwise permeabilities, K + y = K + z , 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., ∝ Re(−u * κ v κ )(dU /dy)dy. Note that Re(·) represents the real component and so the Reynolds shear stress contribution from the NW resolvent mode is proportional to Re(−u * κ v κ ) (Luhar et al. 2014). 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 + x ∝ K + x is larger than that for the turbulent cross-flows, , if the virtual origin for the mean flow is offset from the virtual origin for the turbulence. Additional resolvent-based predictions for K + z > K + x ( + z > + x ) lead to an increase in mode gain (data not shown). Figures 4(g,h) show predictions for NW mode structure for the anisotropic permeable substrate with φ xy = 5.5 and K + y = 1.0 that leads to an increase in drag and forcingresponse gain (see figure 3). Compared to 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 K + y 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 2019). For the substrate with φ xy = 5.5 and K + y = 1.0, this pushes the critical layer for the NW mode outwards toŷ + c ≈ 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 φ xy = 5.5 and K + y = 1.0 in figure 3(a). Note that the NW mode gain decreases monotonically with increasing K + x − K + y for predictions made using the synthetic mean profiles (see figure 3(c)). 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 K + y 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 + ≈ 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 (U + s ≈ + x for small + x ), and turbulence penetration into the permeable medium, which depends on the spanwise slip length ( + z ). To account for these effects, the phase speed can be modified to c + = 10 + + x − + z ≈ 10 + K + x − K + z . Figure 8 in Appendix A 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 Appendix A.
Together, the predictions shown in figure 3 and figure 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 (φ xy 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 K + y . Figure 5. Normalized gain for spanwise-coherent structures plotted as a function of streamwise length λ + x and wave speed c + for substrates with an anisotropy φxy = 5.5. These predictions make use of the synthetic mean profile. Panels (a)-(d) represent structures with spanwise wavenumber κz ≈ 2.3 (λ + z = 500). Panels (e)-(h) represents structures with spanwise wavenumber κz = 0. Red and blue shading respectively represent mode amplification and suppression relative to the smooth-wall case. Wall-normal permeability increases from left to right. Structures with the largest amplification are labeled with a (•) marker.

Spanwise-Coherent Resolvent Modes
In this section, we use the resolvent analysis framework to evaluate the conditions in which energetic spanwise-coherent rollers resembling Kelvin-Helmholtz vortices emerge over anisotropic permeable substrates. Snapshots of the flow field from the simulations of Gómez-de-Segura & García-Mayoral (2019) show that for K + y < 0.4, the flow is dominated by streamwise-elongated structures associated with the NW cycle. As the wallnormal permeability increases to K + y 0.4 the flow becomes dominated by spanwise coherent rollers with an approximate streamwise spacing of λ + x ≈ 100 − 300. Further, using momentum balance arguments, Gómez-de-Segura & García-Mayoral (2019) showed that the Reynolds stress contribution from structures with λ + x ≈ 70 − 320 and λ + z 120 is responsible for the increase in drag observed for substrates with K + y 0.4. We note that spanwise-coherent rollers have also been observed in previous simulations over isotropic porous materials (Breugem et al. 2006;Rosti et al. 2015).
Linear stability analyses successfully predict the emergence of Kelvin-Helmholtz vortices in channel flows over anisotropic permeable substrates (Abderrahaman-Elena & García-Mayoral 2017;Gómez-de-Segura et al. 2018). Gómez-de-Segura et al. (2018) 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 K + x K + y and h + K + y . However, models based on linear stability theory do not accurately predict the threshold value of K + y for the onset of the rollers or the wavelength of the rollers. Linear stability analyses predict structures with streamwise wavelength λ + x ≈ 50 − 70 to be most unstable, while DNS results suggest that structures with λ + x ≈ 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 φ xy = 5.5 and varying K + y . Based on the DNS results of Gómez-de-Segura & García-Mayoral (2019), we limit ourselves to the evaluation of structures with streamwise wavelength λ + x ∈ [50, 500] and mode speed c + ∈ [4, 10]. In addition, we consider structures that are infinitely long in the spanwise direction (κ z = 0) as well as structures with a finite, but large, spanwise wavelength (κ z ≈ 2.3; λ + z = 500). Figures 5(a,e) show that for low wall-normal permeability, K + y = 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., σ κ,p /σ κ,s < 1.25 for both κ z = 0 and κ z ≈ 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 (2019), which show that the Reynolds stress contributions from structures with λ + x ≈ 70 − 300 and λ + z 120 are negligible for materials with K + y 0.31. In other words, spanwisecoherent modes do not have a significant effect on drag for these conditions. However, the amplification of these spanwise-coherent structures increases significantly for K + y 0.39. Figures 5(b,f) show a clear increase in the size and intensity of the region showing mode amplification. For λ + z = 500, the region of highly-amplified structures is localized around (λ + x , c + ) ≈ (130, 7) and the normalized gain for the most amplified structure is σ κ,p /σ κ,s ≈ 2. For spanwise constant structures, the region of high amplification is localized around (λ + x , c + ) ≈ (180, 6.7) and extends to structures with λ + x ≈ 300. The normalized gain for structures in this region is as high as σ κ,p /σ κ,s ≈ 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 Figs. 5(c,d) and (g,h). Indeed, for the case with K + y = 0.66 shown in figure 5(d), the region of high amplification extends from λ + x ≈ 100 − 500 and c + ≈ 4 − 9. The highest-amplification structures in this region show normalized gains σ κ,p /σ κ,s 100.
The substantial increase in amplification of spanwise-coherent structures for K + y 0.39 is again consistent with the DNS results of Gómez-de-Segura & García-Mayoral (2019), which show that the drag-reducing performance of anisotropic substrates begins to degrade for K + y 0.39 due to the additional Reynolds shear stresses contributed by spanwise coherent structures with λ + x ≈ 70 − 300 and λ + z > 120. The DNS results also show the emergence of a new peak in the premultiplied spectra for wall-normal velocity near λ + x ≈ 150 for K + y = 0.39. This is very close to the streamwise wavelength of the most amplified modes in Figs. 5(b,f). Moreover, the streamwise extent of the highintensity region in the DNS premultiplied spectra increases with increasing K + y , which is similar to the expansion of the high-amplification region evident in Figs. 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 λ + x ≈ 180 and c + ≈ 6.7. As expected for Kelvin-Helmholtz 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. Importantly, the model predictions shown in figure 5 for φ xy = 5.5 are representative of what is observed over the permeable substrates with φ xy = 3.6 and φ xy = 11.4 as well. This is illustrated well in figure 7, which shows the normalized gain for the most amplified  Comparison of resolvent-based gain predictions for spanwise-coherent structures with anisotropy ratios φxy = 3.6 ( ), 5.5 ( ), and 11.4 ( ). The maximum normalized gain obtained for resolvent modes with λ + x ∈ [50, 500] and c + ∈ [4, 10] is shown as a function of wall-normal permeability for (a) κz ≈ 2.3 (λ + z = 500) and for (b) κz = 0 (λ + z = ∞). All of these predictions were obtained using the synthetic mean profiles. resolvent mode in the window spanning λ + z ∈ [50, 500] and c + ∈ [4, 10], for both κ z ≈ 2.3 and κ z = 0. For all three anisotropy ratios, the normalized gain show a sharp increase for K + y 0.4. For structures with λ + z = 500 (κ z ≈ 2.3), model predictions for all three anisotropy ratios collapse together nicely ( figure 7(a)). For spanwise-constant structures (κ z = 0), there is a bit more scatter in the normalized gain values but the overall trend remains similar.
The DNS results in figure 3(a) clearly show drag reduction performance deteriorating for K + y 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 (2019) attribute this deterioration in performance to the emergence of high-gain spanwisecoherent structures resembling Kelvin-Helmholtz 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 used an eddy viscosity model, these results also show that the resolvent framework can generate a priori predictions for whether a porous substrate with known 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 ( K + y ≈ 0.4) as well as the streamwise wavelength of the structures that emerge first (λ + x ≈ 150).

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 (2019) 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 Kelvin-Helmholtz 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 near-wall 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, K + x − K + z , 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 K + x and K + z leads to a wall-normal compression of the quasi-streamwise vortices associated with the near-wall cycle, and this compression limits energy extraction from the mean flow. Further, as we show in §3.3, resolvent analysis also predicts the emergence of energetic spanwise rollers with streamwise wavelength λ + x ≈ 150 as the wall-normal permeability increases beyond K + y ≈ 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 for flows over anisotropic permeable substrates. The form of the mean profile can have a significant effect on resolvent-based predictions. 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 K + y 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 K + y ≈ 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 volume-averaged Navier-Stokes 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 also 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.
This material is based on work supported by the Air Force Office of Scientific Research under awards FA9550-17-1-0142 (Program Manager: Dr. Gregg Abate) and FA9550-19-1-7027 (Program Manager: Dr. Douglas Smith). Declaration of interests: the authors report no conflict of interest.

Appendix A.
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 + ≈ 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 + + x − + z ≈ 10 + K + x − K + z . In effect, this modification accounts for the interfacial slip velocity, which depends on the streamwise slip length + x ≈ K + x , and the virtual origin for the turbulence inside the permeable substrate, which depends on the spanwise slip length + z ≈ K + z . The modified phase speed does not lead to significant changes in mode gain compared to the predictions shown in figure 3 for c + ≈ 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 to the predictions shown in figure 3(a), singular value ratios for the substrates with φ xy = 3.6 and φ xy = 5.5 show slightly better agreement with the drag reduction trends from DNS for K + y < 0.6. However, for K + y 0.6, the predictions are more scattered with the modified phase speed. For the substrate with φ xy = 11.4, the modified phase speed leads to reduced mode suppression. Further, the deterioration in performance-in terms of mode suppression-commences at a lower value of wall-normal permeability, K + y ≈ 0.2, compared to 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 K + y , i.e., until the emergence of spanwise rollers resembling Kelvin-Helmholtz vortices in the DNS. Figure 8(c) shows that the initial reduction in mode gain depends primarily on K + x − K + z , 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). Figure 8. Predictions for the NW resolvent mode with a modified phase speed of c + ≈ 10 + K + x − K + z . (a) Comparison between normalized mode gain (σκ,p/σκ,s, black symbols) and drag reduction observed in DNS (∆U + , light gray symbols) by Gómez-de-Segura & García-Mayoral (2019), plotted as a function of K + y . 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 (filled symbols) and the synthetic mean profiles (open symbols) plotted as a function of K + y . (c) Comparison between normalized gains predicted using the DNS mean profiles (filled symbols) and the synthetic mean profiles (open symbols) plotted as a function of K + x − K + y . The black and red dashed lines show linear fits to the initial decrease in normalized gain for the DNS-and synthetic-mean predictions, respectively. For all panels, the symbols represent substrates with φxy = 3.6, symbols represent substrates with φxy = 5.5, and symbols represent substrates with φxy = 11.4.