Resolvent analysis of stratification effects on wall-bounded shear flows

The interaction between shear driven turbulence and stratification is a key process in a wide array of geophysical flows with spatio-temporal scales that span many orders of magnitude. A quick numerical model prediction based on external parameters of stratified boundary layers could greatly benefit the understanding of the interaction between velocity and scalar flux at varying scales. For these reasons, here, we use the resolvent framework to investigate the effects of an active scalar on incompressible wall-bounded turbulence. We obtain the state of the flow system by applying the linear resolvent operator to the nonlinear terms in the governing Navier-Stokes equations with the Boussinesq approximation. This extends the formulation to include the scalar advection equation with the scalar component acting in the wall-normal direction in the momentum equations. We use the mean velocity profiles from a DNS of a stably-stratified turbulent channel flow at varying friction Richardson number. The results obtained from the resolvent analysis are compared to the premultiplied energy spectra, auto-correlation coefficient, and the energy budget terms obtained from the DNS. It is shown that despite using only a very limited range of representative scales, the resolvent model is able to reproduce the balance of energy budget terms as well as provide meaningful insight of coherent structures occurring in the flow. Computation of the leading resolvent models, despite considering a limited range of scales, reproduces the balance of energy budget terms, provides meaningful predictions of coherent structures in the flow, and is more cost-effective than performing full-scale simulations. This quick model can provide further understanding of stratified flows with only information about the mean profile and prior knowledge of energetic scales of motion in the neutrally-buoyant boundary layers.

are heated from above and are usually stably stratified [5,6]. In both the atmosphere and oceans, stratification has a significant effect on turbulence production, propagation, and decay. The interaction between shear-driven turbulence and stratification is a key process in a wide array of relevant geophysical flows for which the spatio-temporal scales span many orders of magnitudes.
Classical understanding of stably-stratified boundary layers is well described in a number of textbooks [7][8][9][10] and reviews [11][12][13]. However, fundamental features of the stablystratified turbulent boundary layer still remain elusive from a modeling standpoint. The strong intermittency observed in stable boundary layers causes the upper portion of the boundary layer to decouple from the near-wall region due to the inhibition in vertical mixing [9,14,15]. Strong stable stratification also significantly changes the flow structures prevalent in a boundary layer with additional features becoming prominent such as large-scale intermittency, gravity waves and Kelvin-Helmholtz instabilities [14], and the near parallel downstream tilting of flow structures [16][17][18].
One way to study the stably-stratified turbulent boundary layer is through on-site experiments. Researchers in the past decades have conducted field experiments in the stablystratified atmospheric boundary layer to study turbulent energy budgets [19], heat and momentum transfer [20], regime characterization [14,21], flow structures [16], and the complexities of atmospheric stable boundary layers [22]. Measurements of turbulence quantities in the ocean near the bottom boundary are difficult to measure and as such the literature is sparse. Smedman et al. [23], using data from a marine coastal experiment over the Baltic sea, found that the near-wall turbulence was virtually independent of forcing from largescale structures embedded in the flow. Experiments performed in the northern bay of San Francisco [24] found that active turbulence is confined near the wall. Additionally, tidal channel experiments [25] demonstrated that the production of turbulent kinetic energy is generally greatest near the bottom boundary while the buoyancy flux is weakest in this region. Still, real-world atmospheric and oceanic boundary layers are complicated by nonturbulent motions occurring simultaneously on a variety of scales, the possible importance of radiative flux divergence of the air within the boundary layer, surface condensation, and variable cloudiness [11,13,26]. In order to isolate instances where the secondary effects are minimized, restrictions on nonstationarity or conditions on the minimum allowed value of turbulence energy may be applied to the data collected. Nonetheless, certain assumptions that are applied for analyses of these real-world stratified boundary layers are not always valid. As such, researchers supplement their work with laboratory experiments as well as simulations.
Laboratory experiments of stratified wall-bounded flows show that buoyancy effects play an important role in the transfer of heat and momentum in both the inner and outer layers of the boundary layer [27][28][29][30][31]. In general, the experiments show that with increasing stratification, the turbulence shear production rate is strongly affected by buoyancy and greatly reduced far from the wall. One measure of stratification strength is the local gradient Richardson number, Ri g . Since shear originates at the wall, the local gradient Richardson number, which is inversely proportional to the shear, is generally smaller in the near-wall region as the shear term overpowers the buoyancy term. The stabilizing effect of stratification has a greater impact farther from the wall. Indeed, works listed here demonstrated that velocity fluctuations become weaker from the wall and in some cases, turbulence intensity is reduced as the buoyancy frequency in the system is increased. Linear inviscid stability analysis [32] showed that there exists a critical value for the gradient Richardson number, Ri g ≥ 0.25, that serves as a sufficient condition for stability. Additionally, the experiments of Komori et al. [30] show that the correlation coefficients associated with the Reynolds shear stress approach zero at values of Ri g ≃ 0.2 − 0.3.
There have been many large-eddy simulations (LES) [33][34][35][36] and direct numerical simulations (DNS) [37][38][39][40][41] of density stratified channel flows. The results support the experimental observations: strengthening the stratification leads to the reduction (or even suppression) of turbulent velocity fluctuations further from the wall. Garg et al. [33] showed in their work that the mean velocity profiles of the stratified channel were similar in the near-wall region but differed in the logarithmic region. The difference is characterized by a reduction in the value of both the slope of the log-law of the mean velocity and the gradient of the mean velocity profile. It should be noted that the authors used the friction Richardson number to categorize the stratification strengths investigated in their simulations and concluded that the friction Richardson number is superior to the local gradient Richardson number in characterizing flow regimes as it is a global flow property.
Performing experiments (both on-site and in laboratories) of stratified wall-bounded turbulence can be challenging for reasons such as topography or secondary effects and simulations suffer from computational constraints. Moreover, laboratory experiments and simula-tions can attain only a limited range of Reynolds and Richardson numbers that are often orders of magnitude smaller than real-world geophysical phenomena. A quick numerical model prediction of key features of stratified boundary layers could greatly benefit the understanding of the interaction between velocity and scalar flux at varying scales. For these reasons, in this paper, we aim to explore the interaction between velocity and scalar fluctuations using the resolvent model [1].
The resolvent model provides an optimal basis, in an energy sense, that allows an in-depth comparison of the underlying mechanisms in the flow. Moreover, the model is computationally efficient with only a singular value decomposition of the largest singular value required to obtain the leading order model. Resolvent analysis has been widely applied to a range of flow configurations to identify dominant flow structures and the underlying forcing, e.g. Ref.
[1, [42][43][44][45][46], and has been reviewed in detail in Ref. [47] and Ref. [48]. We use the model to provide analysis of the flow using only mean quantities, which are easy to obtain even in field experiments, along with knowledge from the energetics of the unstratified case, which is better documented than the stably-stratified case. The predictions from the resolvent model are then compared to the flow statistics from a DNS of a stably-stratified turbulent channel flow. The Reynolds number under consideration in the current study is considerably lower than those observed in geophysical flows, which is dictated by the available DNS data for comparison, rather than by the resolvent model. Resolvent analysis of unstratified wallbounded flows shows that the results of the model are still relevant for moderate Reynolds numbers [49] with the resolvent modes in the logarithmic layer showing self-similar behavior. We expect the capability of the model in stably-stratified regimes to extend to higher Reynolds numbers as well.
The paper is organized as follows. In §II, we introduce the resolvent framework with the inclusion of the scalar advection-diffusion equation and discuss the relevant energy norm, boundary conditions, and computational methods. In §III A, we examine the sensitivity of the low-rank properties of the resolvent operator to the stable stratification strength and compare these properties with the most energetic scales in each flow. In §III B, we analyze the characteristics of the forcing and response modes of both velocity and scalar.
We compare the mode shapes with correlations obtained from DNS data. In §III C, we study the turbulent kinetic energy budget in the resolvent formulation and compare the results with the energy budget obtained from the DNS data. Finally, our conclusions on the application of the resolvent framework to a stably-stratified boundary layer are given in §IV. We consider a density-stratified turbulent channel flow where the density acts in the direction of gravitational acceleration. We use a Cartesian co-ordinate system x = (x, y, z) such that the force of gravity acts in the −y direction, with x, y and z being the streamwise, wall-normal and spanwise directions, respectively. The governing equations are given by the non-dimensional Navier-Stokes equation under the Boussinesq approximation, Here, u = (ũ,ṽ,w) is the instantaneous velocity vector in the reference system (x, y, z), t is time, p is the kinematic pressure field that remains after removing the part that is in hydrostatic balance with the mean density field, ρ is the density deviation from the reference density ρ 0 ( ρ ≪ ρ 0 ), and e y is the unit vector acting in the y-direction. The velocity and length scales are non-dimensionalized using the friction velocity u τ and channel half-height δ, respectively, and the density is non-dimensionalized using ∆ρ, the difference in density between the two channel walls. We define the walls to be located at y = 0 and y = 2. The non-dimensional quantities are given by the Reynolds, Prandtl and Richardson numbers, defined as where ν is the kinematic viscosity, γ is the molecular diffusivity of density, and g is the acceleration due to gravity.

B. Resolvent framework with an active scalar
The total fields u, p and ρ can be split into mean and fluctuating parts as p(x, t) = p(y) + p(x, t), where the mean is taken in the homogeneous directions, x and z, and time. Note that u = (ū,v,w) andv =w = 0. We substitute the decomposed variables into Eq.
(1) to obtain the fluctuation equations where f u = −u · ∇u and f ρ = −u · ∇ρ are the nonlinear terms.
Taking the Fourier transform of the fluctuation equations above in homogeneous directions and time, the variables can be expressed as  for k = (k x , k z , ω) = (0, 0, 0), where (·) denotes the Fourier transformed variables. Here, the streamwise and spanwise wavenumbers are k x and k z , respectively, and ω is the temporal frequency defined as ω = ck x , where c is the wavespeed. The streamwise and spanwise wavelengths are defined as λ x = 2π/k x and λ z = 2π/k z , respectively. Critical-layers can be identified when the wavespeed c is equivalent to the mean velocity, i.e. y c is the critical layer location for wavespeed c = u(y c ). Assuming the mean velocity and density profiles are known, the fluctuations equations are expressed compactly in a linear equation as where we defineq = [ûvŵpρ] T as the state vector andf = [f ufvfw 0f ρ ] T as the forcing vector. The linear operator is given by where The block matrix A describes the linear dynamics of the system. Equation (6) can be rearranged to yieldq where H(k) = (−iωI − A) −1 is the resolvent of the linear operator and I is the identity matrix. A related analysis has been performed in Ref. [50].
From Eq. (9), we wish to find a decomposition of the resolvent operator that enables us to identify high gain input and output modes with respect to the linear operator. For resolvent analysis, this is given by the Schmidt decomposition. However, this decomposition must be accompanied by a choice of inner product and the corresponding norm. The natural and physically meaningful norm is given by the non-dimensionalized energy norm, which is the sum of kinetic and potential energies [51,52] where (·) * denotes the conjugate transpose.
We perform the Schmidt decomposition of the resolvent operator H to generate a basis based on the most highly amplified forcing and response directions such that where the right and left Schmidt bases (or singular vectors in the discrete case) are given bŷ φ j andψ j along with their corresponding gains σ j . The singular values are in descending order such that σ 1 ≥ σ 2 ≥ · · · ≥ 0. The forcing and resolvent modes are orthonormal such where δ jk denotes the Kronecker delta. The basis pair defined above is used to decompose the nonlinear forcing and response field at a specified wavenumber triplet aŝ Here, χ j is a projection variable that is obtained by projecting the nonlinear forcing onto the forcing modes, and subsequently use to weight the response modes. Note that the largest energy is obtained when the forcing is aligned with the leading singular vector, i.e. when C. Computational approach

Mean velocity and density profiles
Mean velocity and density profiles are required to close the resolvent model. We obtain the one-dimensional mean velocity and density profiles from a DNS of a stratified turbulent channel at Re τ = 180 for a wide range of Ri τ . The simulations are performed by discretizing the incompressible Navier-Stokes equations with a staggered, second-order accurate, central finite-difference method in space [53], and an explicit third-order accurate Runge-Kutta method for time advancement [54]. The system of equations is solved via an operator splitting approach [55]. The code has been verified for neutrally-buoyant cases in Ref. [56,57].   toward the wall according to a hyperbolic tangent distribution with min(∆y + ) = 0.31 and max(∆y + ) = 5.19, where the superscript + indicates length scales in wall units normalized by ν/u τ rather than δ. A constant pressure gradient is applied to drive the flow. The simulation was run over 100 eddy-turnover times, defined as δ/u τ , after transients.
The work of García-Villalba & delÁlamo [41] at Re τ = 180 is used to validate the results. The comparison of a few key quantities is shown in Table I, which indicate a good agreement for all Richardson numbers. The mean and root-mean-squared streamwise velocity and density profiles are shown in Fig. 1 for all current cases and select cases from Ref. [41]. The profiles show good agreement among all statistics.

Resolvent mode computation
The Schmidt decomposition of the resolvent operator outlined in §II B is numerically implemented as the singular value decomposition (SVD). We solve the discrete equations using a spectral collocation method with the number of points in the wall-normal direction given by N y , thus limiting the number of singular values to 5N y because the state vector  [41] for Ri τ = 0, 18, 120 (dashed lines darker to lighter). The friction density is defined as ρ τ = q w /u τ , where q w is the density flux at the wall.
q ∈ C 5Ny×1 . In this study, after conducting a grid convergence study examining the singular values, we selected a wall-normal grid resolution of N y = 400. Thus, the computational cost of the resolvent mode computation is at most O(N 3 y ) (less if randomized algorithms are employed [49,58]), often only requiring a leading order singular value decomposition (see §III A for more information), and can be performed in seconds on a personal computer.
The discretized linear operator is constructed using Chebyshev differentiation matrices and is shifted to integrate between y ∈ [0, 2] rather than y ∈ [−1, 1]. The mean velocity and density profiles obtained from DNS as well as their wall-normal derivatives are interpolated to the Chebyshev grid points to form the resolvent operator as in Eq. (7). The no-slip and nopenetration boundary conditions for the fluctuating velocities and density, i.e. u, v, w, ρ = 0, are applied at the walls.
In the case of a turbulent channel, due to the symmetry in the geometry, the resolvent modes appear in pairs that can be linearly combined to produce symmetric and antisymmetric modes. Depending on the support of these modes, the singular values may be identical or similar in magnitude. For the results in the following sections, only results in the bottom half-channel will be shown, but the corresponding upper half-channel results are analogous in all cases.

III. RESULTS
In this section, we explore how the resolvent analysis provides insight to changes in flow characteristics with increasing stratification from only a limited range of representative scales. We compare (i) the resolvent energy spectra, obtained from the ratio of the energy in the leading resolvent response mode to the total response, (σ 2 1 + σ 2 2 )/ j σ 2 j , to the premultiplied energy spectra of the DNS, (ii) the structure identified by the leading resolvent mode to the correlation computed from DNS, and (iii) the energy budgets of the resolvent modes to that of the DNS.
In order for full representation of the system, a wide range of scales as well as information of all other subsequent modes in addition to the leading resolvent modes are necessary [1,47].
However, the goal here is to provide a quick model for characterising the flow. The simplest and quickest model can be provided via a rank-one approximation, where only the leading resolvent mode is computed. Thus, our focus will be on the representation given by the leading resolvent mode for a limited number of scales.

A. Resolvent energy spectra
The resolvent norm is the principal singular value, σ 2 1 + σ 2 2 in this case, of the resolvent operator H, and quantifies the system's sensitivity to temporal forcing [59]. The energetic contribution from broadband forcing is quantified as the square of the resolvent norm.
The resolvent operator H can be described as low-rank if the majority of its response to broadband forcing in the wall-normal direction is captured by the first few response modes.
Theoretically, there are an infinite number of singular values and corresponding modes because the wall-normal coordinate is continuous. However, not all of the singular vectors are energetically significant. As described in §II B, a self-sustaining representation of the flow will correspond to a weighted assembly of forcing modes rather than a broadband forcing [46]; however, past studies have showed that broadband forcing is successful in identifying  [49] showed that the first two resolvent modes account for more than 80% of the total response in a channel. Bae et al. [44] investigated the low-rank nature of a compressible turbulent boundary layer and highlighted the similarities in the region where the low-rank approximation is valid for the incompressible regime.
Assuming the resolvent operator is low-rank (σ 1 ≃ σ 2 ≫ σ 3 ) allows us to approximate the operator as for each k since most of the energy in the system is modelled by the principal singular value. The low-rank behavior of H is typically representative of there being a dynamically significant physical, spatio-temporal structure at the scale dictated by k.
To study the variation in the low-rank behavior for different magnitudes of stratification, we plot the energetic contribution of the principal response mode to the total response in the model for a given k quantified by (σ 2 1 + σ 2 2 )/Σ j σ 2 j for a range of wall-parallel wavelengths (Fig. 2). The leading response modes account for more than 80% of the total response over a large range of homogeneous wavelengths for the three wavespeeds selected.
The range of wavenumbers for which the resolvent operator is low-rank changes significantly with stratification. In the neutrally-buoyant case (Ri τ = 0), we see that H is low-rank in a range of moderate-to-large streamwise wavelengths. For the neutrally-buoyant case, it is known that the low-rank region coincides with the most energetic wavenumbers from the premultiplied energy spectra of a turbulent channel [49]. As the friction Richardson number first increases, the low-rank behavior shifts to only a small range of streamwise wavelengths. We see a similar phenomenon in the premultiplied streamwise energy spectra from the DNS (Fig. 3), where with increasing Ri τ , the larger streamwise wavelength content is suppressed. This was also observed in the premultiplied energy spectra of García-Villalba & delÁlamo [41] for a wider range of Re τ and Ri τ .
However, after Ri τ = 18, the low-rank behavior of the principal resolvent modes intensifies along a vertical band λ x /δ ≥ 1 until the system becomes low-rank at large spanwise wavelengths with almost no low-rank behavior below the green dashed line in Fig. 2 (λ x = 15λ z , 10λ z and 5λ z for y + = 15, 30 and 100, respectively). This seems to indicate a low-rank  behavior in structures that are descriptive of quasi-two-dimensional flow where λ z ≫ λ x .
Hopfinger [61] details the emergence of two-dimensional modes for a variety of flows with strong stratification. Moreover, Mahrt [13] alludes to the emergence of two-dimensional modes (often referred to as pancake modes) owing to the conversion of vertical kinetic energy to potential energy in the presence of strong stable stratification. The premultiplied energy spectra for higher Ri τ indicate high energy in the vertical band as well [41].

B. Mode shapes
In order to study the flow structures, we compute the resolvent response modes for a set of wave parameters. The most energetic scales for the various Ri τ under consideration for the different wall-normal heights still coincide with the neutrally-buoyant case (Fig. 3), falling in the low-rank region despite the fact that including the scalar advection-diffusion equation in the governing equations changes the wavelengths at which the resolvent operator is low-rank ( Fig. 2). In this section, we study the resolvent response mode shapes for these wavenumber and wavespeed combinations. The list of mode combinations under consideration is listed in Table II. In particular, mode E1 is the most energetic mode for y + = 15, E2 for y + = 30 and E3 for y + = 100.
The predictive capabilities of the resolvent modes are first shown through the amplitudes of the leading resolvent response modes (Fig. 4) of the streamwise velocity and density. The resolvent modes compare well to the streamwise and density turbulence intensities in Fig.   3(c,d). The streamwise root-mean-square (r.m.s.) quantities and resolvent amplitudes show no variation among different Richardson numbers closer to the wall and increase slightly with Ri τ farther away from the wall. On the other hand, the density r.m.s. and resolvent amplitudes decrease significantly with Richardson number at all wall-normal heights. Despite only using the leading resolvent mode, the relative magnitude at each corresponding wall-normal height is well captured for the range of Richardson numbers considered here.
Additionally, we examine the response mode shapes in two dimensions for the different regions and compare the structures observed in the resolvent modes with the autocorrelation coefficient from the DNS data. We first define the streamwise auto-covariance aŝ where q is a generic variable of zero mean and · is the expected value. The auto-covariance in physical space, R qq (∆x, y, y ′ , ∆z), is obtained as the inverse Fourier transform ofR, where ∆x = x − x ′ and ∆z = z − z ′ are the distances between the two points in the homogeneous directions. The autocorrelation coefficient, is obtained by normalising the covariance with the product of the standard deviations, ς, at the two points involved in the measurements, which is the normalization adopted by most researchers [62][63][64][65][66][67].
The two-dimensional structures of mode E1, which coincides with the size of the near the wall structures observed previously in experiments and simulations [68,69], are plotted in We plot the resolvent response modes for the wavenumbers and wavespeed corresponding to E2 (Fig. 7) and the correlations for y ′+ = 30 at ∆z = 0 (Fig. 8). The results are similar to that of E1, since the velocity response modes do not vary across Ri τ , but a difference is observed in the density modes as a phase change along y. The correlation for the density modes are both wall-detached in the Ri τ = 0 and Ri τ = 100 case, although the centre of the density correlation for the Ri τ = 100 case lies farther away from the wall.
The biggest difference in the resolvent response modes for the different Richardson numbers can be seen for the wavenumber and wavespeed corresponding to E3. We plot the resolvent response modes for the wavenumbers and wavespeed corresponding to E3 (Fig. 9) and the correlations for y ′+ = 100 at ∆z = 0 (Fig. 10).
Here, all resolvent modes show significant differences in the stratified case compared to the unstratified case. In particular, the backwards tilting of the wall-normal velocity modes, the forward tilting of the density modes, as well as the phase difference across y of the density mode are pronounced. These phenomena occur in the correlations as well. There is noticeable backwards tilting in the C vv term and forwards tilting in the C ρρ term for Ri τ = 100 compared to the neutrally stratified case. The biggest differences come in the form of the wall-normal and density models because they are coupled through the Richardson number in the stratified Navier-Stokes equations.

C. Energy balance at selected scales
Finally, we study the energy budget terms of the stratified channel. We define the production, transport, buoyancy flux, and viscous dissipation budget terms in the resolvent formulation [50,70] as where χ j , σ j ,ψ j andφ j are functions of k and the subscript u, v, w, ρ indicate the corresponding components of the response or forcing mode. To get a global sense of the energy balance, the equations above are integrated over all wavenumber triplets. Here, we will examine only the principal resolvent mode contribution to the local components of the total budgets for particular k, defined as T (y, k) = R σ 1 D y φ * 1,uψ 1,v +φ * 1,vψ 1,v +φ * V(y, k) = R 1 Re τ σ 2 1 ψ * 1,u∆ψ 1,u +ψ * 1,v∆ψ 1,v +ψ * 1,w∆ψ 1,w .
The results for wavenumber combinations E1, E2 and E3 are shown in Fig. 11. Since the wavenumber combinations E1, E2 and E3 are the most energetic at each wavespeed, we predict that the local components of the budget term should indicate the overall trend of the total budget term at the corresponding wall-normal height. These quantities are compared to the energy budget computed from the DNS, shown in Fig. 12.
The trends observed in the energy budget computed from the DNS are also recovered in the resolvent budgets. The production is mostly balanced by viscous dissipation and has larger magnitudes compared to the transport (approximately 10% of the production term) or buoyancy flux (approximately 0.1-1%, depending on Ri τ of the production term) terms.
Comparing the quantities at the wall-normal heights of interest, we see that at y + = 15, there is little variation in the production and viscous dissipation terms in both DNS and resolvent modes. The difference in relative magnitude over the various values of Ri τ increases farther away from the wall, and at y + = 100, the production (and viscous diffusion) of the Ri τ = 100 case is double the production (and viscous diffusion) of the neutrally-buoyant case in both the DNS and resolvent.
Direct comparison of the integrated magnitudes is more difficult for the transport and buoyancy flux terms as they are not uniformly positive or negative. However, this indicates that, locally, buoyancy flux acts as a energy transfer term, much like the turbulent transport, as the term adds energy in one wall-normal location and removes it from another. Because the DNS energy budget is integrated for all spatio-temporal scales, it is impossible to deduce that the buoyance flux term acts as a local energy transfer term from Fig. 12, which shows a net negative energy balance from B at all wall-normal locations. In contrast, the resolvent buoyancy flux term indicates a non-monotonic distribution of energy in the wall-normal direction. Similar results could be obtained through spatio-temporal deconstruction of the These results can be better quantified by plotting the values at each wall-normal location normalized by the peak production at y + = 15 for each case, as shown in Fig. 13. This shows that the overall trend of the budget terms are well captured by the resolvent budget terms, with the exception of the transport term close to the wall. This discrepancy may be attenuated by integrating over more wavespeeds.
Note that the results are not expected to match that of DNS for all scales as the the energy captured in the wall-parallel resolvent modes are known to be overpredicted and the energy captured in the Reynolds stress and wall-normal resolvent modes underpredicted. This is a known issue for the resolvent analysis in the primitive variables due to the competing mechanisms of the Squire modes with the Orr-Sommerfeld modes [72,73]. Additionally, the underprediction of energy captured in the Reynolds stress and wall-normal resolvent modes could explain the underprediction of the transport term close to the wall. Crucially, though, the most energetic scale can reproduce the integrated effect of all scales, which enables a quick predictive model of stratified boundary layers. Ri τ = 0. Symbols are P, circles; T , triangles; V, crosses; and B, asterisks.

IV. CONCLUSIONS
The resolvent framework for the Navier-Stokes equations with the Boussinesq approximation was applied to a stratified turbulent boundary layer. Computation of the leading resolvent modes is more cost-effective than performing a full-scale simulation or experiment, while being able to provide information on the flow. This quick model can provide meaningful insight into stratified flows with only information about the mean profile and prior knowledge of energetic scales of motion in the neutrally-buoyant boundary layers.
The results show that despite using only a very limited range of representative scales, the resolvent model was able to reproduce the relative magnitude of turbulence intensities and the balance of the energy budget as well as provide meaningful analysis of structures in the flow. We studied the amplitude of the resolvent response modes and their two-dimensional mode shapes of the rank-one approximation, which were then compared to the turbulence intensities and the two-dimensional auto-correlation of the velocity and density fields of the DNS, respectively. The resolvent response modes were able to predict the relative variation in turbulence intensities as a function of wall-normal distance and Richardson number for the Ri τ under consideration in this study. The two-dimensional mode shapes also provided insight into how the auto-correlation coefficient might shift as a function of Ri τ . Finally, the energy budget terms for the turbulent kinetic energy of the system were computed both using the rank-one approximation of the resolvent analysis and the DNS data. Again, the resolvent energy budget predicts well the relative distribution of energy between production, dissipation, transport, and buoyancy flux as a function of wall-normal distance and Richardson number.
In the current study, the resolvent model was closed using mean velocity and density profiles obtained from DNS. The computational cost of calculating the forcing and response modes at certain scales was on the order of seconds on a laptop. Therefore, by obtaining only mean velocity and scalar profiles we could generate the salient modal structure for a given stratified wall-bounded flow. The next steps involve using in-situ data to generate modes that are representative of flow phenomena observed in nature.