# Resolvent analysis of stratification effects on wall-bounded shear flows

M. A. Ahmed<EMAIL_ADDRESS>Graduate Aerospace Laboratories, California
Institute of Technology, Pasadena, CA 91125, USA  H. J. Bae<EMAIL_ADDRESS>Harvard University, Cambridge, MA 02139, USA  Graduate Aerospace Laboratories,
California Institute of Technology, Pasadena, CA 91125, USA  A. F. Thompson
<EMAIL_ADDRESS>Geophysical and Planetary Sciences, California Institute
of Technology, Pasadena, CA 91125, USA  B. J. McKeon<EMAIL_ADDRESS>Graduate Aerospace Laboratories, California Institute of Technology, Pasadena,
CA 91125, USA

###### Abstract

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 [1] 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 [2]. We use the mean velocity
profiles from a direct numerical simulation (DNS) of a stably-stratified
turbulent channel flow at varying friction Richardson number $Ri_{\tau}$. 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.

##  I Introduction

Stable boundary layers can be generated by the advection of warm air over a
colder surface. Stably-stratified atmospheric boundary layers are observed
during clear nights as a result of radiative cooling of the ground surface [3,
4]. Oceans, unlike the lower atmosphere, 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 stably-stratified 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 stably-stratified atmospheric boundary layer to study
turbulent energy budgets [19], heat and momentum transfer [20], regime
characterization [21, 14], 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 large-scale 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 non-turbulent 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
[26, 11, 13]. 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}\geq 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}\simeq 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 simulations 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 wall-bounded 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.1, 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.2, 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.3, 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.

##  II Modeling active scalar dynamics in the Navier-Stokes equations

###  II.1 Navier-Stokes equation with active scalar

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 $\bm{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,

$\displaystyle\frac{\partial\widetilde{\bm{u}}}{\partial
t}+(\widetilde{\bm{u}}\cdot\nabla)\widetilde{\bm{u}}$
$\displaystyle=-\nabla\widetilde{p}+\frac{\nabla^{2}\widetilde{\bm{u}}}{Re_{\tau}}-Ri_{\tau}\widetilde{\rho}\bm{e}_{y},$
(1a) $\displaystyle\frac{\partial\widetilde{\rho}}{\partial
t}+(\widetilde{\bm{u}}\cdot\nabla)\widetilde{\rho}$
$\displaystyle=\frac{\nabla^{2}\widetilde{\rho}}{Re_{\tau}Pr},$ (1b)
$\displaystyle\nabla\cdot\widetilde{\bm{u}}$ $\displaystyle=0.$ (1c)

Here, $\widetilde{\bm{u}}=(\tilde{u},\tilde{v},\tilde{w})$ is the
instantaneous velocity vector in the reference system $(x,y,z)$, $t$ is time,
$\widetilde{p}$ is the kinematic pressure field that remains after removing
the part that is in hydrostatic balance with the mean density field,
$\widetilde{\rho}$ is the density deviation from the reference density
$\rho_{0}$ ($\widetilde{\rho}\ll\rho_{0}$), and $\bm{e}_{y}$ is the unit
vector acting in the $y$-direction. The velocity and length scales are non-
dimensionalized using the friction velocity $u_{\tau}$ and channel half-height
$\delta$, respectively, and the density is non-dimensionalized using
$\Delta\rho$, 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

$Re_{\tau}=\frac{u_{\tau}\delta}{\nu},\qquad Pr=\frac{\nu}{\gamma},\qquad
Ri_{\tau}=\frac{g\Delta\rho\delta}{\rho_{0}u_{\tau}^{2}},$ $None$

where $\nu$ is the kinematic viscosity, $\gamma$ is the molecular diffusivity
of density, and $g$ is the acceleration due to gravity.

###  II.2 Resolvent framework with an active scalar

The total fields $\widetilde{\bm{u}}$, $\widetilde{p}$ and $\widetilde{\rho}$
can be split into mean and fluctuating parts as

$\displaystyle\widetilde{\bm{u}}(\bm{x},t)$
$\displaystyle={\bm{\overline{u}}}(y)+\bm{u}(\bm{x},t),$ (3a)
$\displaystyle\widetilde{p}(\bm{x},t)$
$\displaystyle=\overline{p}(y)+p(\bm{x},t),$ (3b)
$\displaystyle\widetilde{\rho}(\bm{x},t)$
$\displaystyle=\overline{\rho}(y)+\rho(\bm{x},t),$ (3c)

where the mean is taken in the homogeneous directions, $x$ and $z$, and time.
Note that ${\bm{\overline{u}}}=(\bar{u},\bar{v},\bar{w})$ and
$\bar{v}=\bar{w}=0$. We substitute the decomposed variables into Eq. (1) to
obtain the fluctuation equations

$\displaystyle\partial_{t}\bm{u}+(\overline{u}\cdot\nabla)\bm{u}+(\bm{u}\cdot\nabla)\overline{u}$
$\displaystyle=-\nabla
p+\frac{\nabla^{2}\bm{u}}{Re_{\tau}}-Ri_{\tau}\rho\bm{e}_{y}+\bm{f}_{\bm{u}}$
(4a)
$\displaystyle\partial_{t}\rho+(\overline{u}\cdot\nabla)\rho+(\bm{u}\cdot\nabla)\overline{\rho}$
$\displaystyle=\frac{\nabla^{2}\rho}{Re_{\tau}Pr}+f_{\rho},$ (4b)
$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0,$ (4c)

where $\bm{f}_{\bm{u}}=-\bm{u}\cdot\nabla\bm{u}$ and
$f_{\rho}=-\bm{u}\cdot\nabla\rho$ are the nonlinear terms.

Taking the Fourier transform of the fluctuation equations above in homogeneous
directions and time, the variables can be expressed as

$\begin{bmatrix}\bm{u}(x,y,z,t)\\\ p(x,y,z,t)\\\
\rho(x,y,z,t)\end{bmatrix}=\iiint^{\infty}_{-\infty}\begin{bmatrix}\hat{\bm{u}}(y;k_{x},k_{z},\omega)\\\
\hat{p}(y;k_{x},k_{z},\omega)\\\
\hat{\rho}(y;k_{x},k_{z},\omega)\end{bmatrix}e^{\text{i}(k_{x}x+k_{z}z-\omega
t)}dk_{x}dk_{z}d\omega,$ (5)

for $\bm{k}=(k_{x},k_{z},\omega)\neq(0,0,0)$, where $(\hat{\cdot})$ denotes
the Fourier transformed variables. Here, the streamwise and spanwise
wavenumbers are $k_{x}$ and $k_{z}$, respectively, and $\omega$ is the
temporal frequency defined as $\omega=ck_{x}$, where $c$ is the wavespeed. The
streamwise and spanwise wavelengths are defined as $\lambda_{x}=2\pi/k_{x}$
and $\lambda_{z}=2\pi/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=\overline{u}(y_{c})$. Assuming the
mean velocity and density profiles are known, the fluctuations equations are
expressed compactly in a linear equation as

$-\text{i}\omega\hat{\bm{q}}-\mathcal{A}\hat{\bm{q}}=\hat{\bm{f}},$ (6)

where we define
$\hat{\bm{q}}=[\hat{u}\;\hat{v}\;\hat{w}\;\hat{p}\;\hat{\rho}]^{T}$ as the
state vector and
$\hat{\bm{f}}=[\hat{f}_{u}\;\hat{f}_{v}\;\hat{f}_{w}\;0\;\hat{f}_{\rho}]^{T}$
as the forcing vector. The linear operator is given by

$\mathcal{A}=\begin{pmatrix}A&-\partial\overline{u}/\partial
y&0&-\text{i}k_{x}&0\\\ 0&A&0&-D_{y}&-Ri_{\tau}\\\ 0&0&A&-\text{i}k_{z}&0\\\
-\text{i}k_{x}&-D_{y}&-\text{i}k_{z}&0&0\\\
0&-\partial\overline{\rho}/\partial y&0&0&A_{\rho}\end{pmatrix},$ (7)

where

$\displaystyle A$
$\displaystyle=-\text{i}k_{x}\overline{u}+\frac{\hat{\Delta}}{Re_{\tau}},$
(8a) $\displaystyle A_{\rho}$
$\displaystyle=-\text{i}k_{x}\overline{u}+\frac{\hat{\Delta}}{Re_{\tau}Pr},$
(8b)

$D_{y}$ is the wall-normal derivative operator and $\hat{\Delta}\equiv
D_{yy}-k_{\perp}^{2}$ is the Laplacian with
$k_{\perp}^{2}=k_{x}^{2}+k_{z}^{2}$. The block matrix $\mathcal{A}$ describes
the linear dynamics of the system. Equation (6) can be rearranged to yield

$\hat{\bm{q}}\;=\;\mathcal{H}(\bm{k})\;\hat{\bm{f}},$ (9)

where $\mathcal{H}(\bm{k})=(-\text{i}\omega I-\mathcal{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]

$\frac{1}{2}\|\bm{q}\|^{2}_{E}=\frac{1}{2}(\bm{q},\bm{q})_{E}=\frac{1}{2}\int_{0}^{2}\left(u^{*}u+v^{*}v+w^{*}w+Ri_{\tau}(\rho^{*}\rho)\right)dy,$
(10)

where $(\cdot)^{*}$ denotes the conjugate transpose.

We perform the Schmidt decomposition of the resolvent operator $\mathcal{H}$
to generate a basis based on the most highly amplified forcing and response
directions such that

$\mathcal{H}(\bm{k})=\sum_{j=1}^{\infty}\sigma_{j}(\bm{k})\bm{\hat{\psi}}_{j}(y;\bm{k})\bm{\hat{\phi}}^{*}_{j}(y;\bm{k}),$
(11)

where the right and left Schmidt bases (or singular vectors in the discrete
case) are given by $\bm{\hat{\phi}}_{j}$ and $\bm{\hat{\psi}}_{j}$ along with
their corresponding gains $\sigma_{j}$. The singular values are in descending
order such that $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq 0$. The forcing and
resolvent modes are orthonormal such that

$(\bm{\hat{\phi}}_{j},\bm{\hat{\phi}}_{k})_{E}=(\bm{\hat{\psi}}_{j},\bm{\hat{\psi}}_{k})_{E}=\delta_{jk},$
(12)

where $\delta_{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 as

$\displaystyle\hat{\bm{f}}(y;\bm{k})$
$\displaystyle=\sum_{j=1}^{\infty}\bm{\hat{\phi}}_{j}(y;\bm{k})\chi_{j}(\bm{k}),$
(13a) $\displaystyle\hat{\bm{q}}(y;\bm{k})$
$\displaystyle=\sum_{j=1}^{\infty}\chi_{j}(\bm{k})\sigma_{j}(\bm{k})\bm{\hat{\psi}}_{j}(y;\bm{k}).$
(13b)

Here, $\chi_{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 $\chi_{j}=\delta_{j1}$.

###  II.3 Computational approach

####  II.3.1 Mean velocity and density profiles

Table 1: Comparison of our DNS and the results of García-Villalba & del Álamo [41] denoted under columns titled GV11, both at $Re_{\tau}=180$. $Re_{B}$ is the bulk Reynolds number defined as $u_{B}\delta/\nu$ where the bulk velocity is $u_{B}=\int_{0}^{2}\overline{u}dy/2$. $Ri_{B}$ is the bulk Richardson number which is defined as $Ri_{B}=Ri_{\tau}(u_{\tau}/u_{B})^{2}/2$. $Nu$ is the Nusselt number defined as $Nu=2\delta q_{w}/(\gamma\Delta\rho)$. For laminar flow $Nu=1$. | $Re_{B}$ | $Ri_{B}$ | $Nu$
---|---|---|---
$Ri_{\tau}$ | GV11 | DNS | GV11 | DNS | GV11 | DNS
0 | 2820 | 2823 | 0.000 | 0.000 | 6.03 | 6.08
10 | - | 2970 | - | 0.018 | - | 4.78
18 | 3043 | 3060 | 0.031 | 0.031 | 4.02 | 4.15
60 | 3436 | 3473 | 0.082 | 0.081 | 2.80 | 2.82
100 | - | 3850 | - | 0.109 | - | 2.37

(a)  (b)

(c)  (d)

Figure 1: Mean (a) streamwise velocity and (b) density profiles and root-mean-
square (r.m.s.) (c) streamwise velocity and (d) density profiles from the
current DNS for $Ri_{\tau}=0,10,18,60,100$ (solid lines darker to lighter),
compared to the mean profiles of Ref. [41] for $Ri_{\tau}=0,18,120$ (dashed
lines darker to lighter). The friction density is defined as
$\rho_{\tau}=q_{w}/u_{\tau}$, where $q_{w}$ is the density flux at the wall.

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_{\tau}=180$ for a wide range of
$Ri_{\tau}$. 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].

Periodic boundary conditions are imposed in the streamwise and spanwise
directions, the no-slip and no-penetration condition with $\tilde{\rho}=0$ is
applied at the bottom boundary, and a no-slip and no-penetration condition
with $\tilde{\rho}=1$ is applied at the top boundary. The streamwise, wall-
normal, and spanwise domain sizes are $4\pi$, $2$, and $2\pi$ respectively.
The grid spacings in the streamwise and spanwise directions are uniform with
$\Delta x^{+}=8.8$ and $\Delta z^{+}=4.4$; non-uniform meshes are used in the
wall-normal direction, with the grid stretched toward the wall according to a
hyperbolic tangent distribution with $\min(\Delta y^{+})=0.31$ and
$\max(\Delta y^{+})=5.19$, where the superscript $+$ indicates length scales
in wall units normalized by $\nu/u_{\tau}$ rather than $\delta$. A constant
pressure gradient is applied to drive the flow. The simulation was run over
$100$ eddy-turnover times, defined as $\delta/u_{\tau}$, after transients.

The work of García-Villalba & del Álamo [41] at $Re_{\tau}=180$ is used to
validate the results. The comparison of a few key quantities is shown in Table
1, 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.

####  II.3.2 Resolvent mode computation

The Schmidt decomposition of the resolvent operator outlined in §II.2 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
$\hat{\bm{q}}\in\mathbb{C}^{5N_{y}\times 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_{y}^{3})$ (less if randomized
algorithms are employed [49, 58]), often only requiring a leading order
singular value decomposition (see §III.1 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\in[0,2]$ rather than
$y\in[-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 no-penetration
boundary conditions for the fluctuating velocities and density, i.e.
$u,v,w,\rho=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,
$(\sigma_{1}^{2}+\sigma_{2}^{2})/\sum_{j}\sigma_{j}^{2}$, 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.

###  III.1 Resolvent energy spectra

The resolvent norm is the principal singular value,
$\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}}$ in this case, of the resolvent operator
$\mathcal{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 $\mathcal{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.2, 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 the important component of the flow, e.g. Ref. [1, 44]. McKeon &
Sharma [1] demonstrated that the characteristics of the leading response modes
for a range of wavenumber-frequency combinations agree with experimental
observations in pipe flow and with scaling concepts in wall-bounded
turbulence. Moarref _et al._ [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
($\sigma_{1}\simeq\sigma_{2}\gg\sigma_{3}$) allows us to approximate the
operator as

$\mathcal{H}(\bm{k})\;\approx\;\sigma_{1}\;\bm{\hat{\psi}}_{1}\;\bm{\hat{\phi}}^{*}_{1}+\sigma_{2}\;\bm{\hat{\psi}}_{2}\;\bm{\hat{\phi}}^{*}_{2},$
(14)

for each $\bm{k}$ since most of the energy in the system is modelled by the
principal singular value. The low-rank behavior of $\mathcal{H}$ is typically
representative of there being a dynamically significant physical, spatio-
temporal structure at the scale dictated by $\bm{k}$.

Figure 2: Contour plots depicting the energy contained in the leading response
mode relative to the total response,
$(\sigma_{1}^{2}+\sigma_{2}^{2})/\Sigma_{j}\sigma_{j}^{2}$, for different
streamwise and spanwise wavelengths at (a) $c=\overline{u}(y^{+}=15)$, (b)
$c=\overline{u}(y^{+}=30)$ and (c) $c=\overline{u}(y^{+}=100)$ for
$Ri_{\tau}=$ 0, 10, 18, 60, 100 (top to bottom). Green dashed lines are (a)
$\lambda_{x}=15\lambda_{z}$, (b) $\lambda_{x}=10\lambda_{z}$ and (c)
$\lambda_{x}=5\lambda_{z}$.  Figure 3: Contour plots depicting the
premultiplied streamwise kinetic energy spectra as functions of the streamwise
and spanwise wavelengths obtained from DNS at (a) $y^{+}=15$, (b) $y^{+}=30$,
and (c) $y^{+}=100$ for $Ri_{\tau}=0$ (solid line), $Ri_{\tau}=60$ (dashed
line) and $Ri_{\tau}=100$ (dotted line). The shaded contours are from the
$Re_{\tau}=180$ neutral channel [60]. The levels plotted are $0.1,0.3,0.5$
times the maximum value of the corresponding spectrum.

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 $\bm{k}$ quantified by
$(\sigma_{1}^{2}+\sigma_{2}^{2})/\Sigma_{j}\sigma_{j}^{2}$ 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_{\tau}=0$), we see that $\mathcal{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_{\tau}$, 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_{\tau}$ and
$Ri_{\tau}$.

However, after $Ri_{\tau}=18$, the low-rank behavior of the principal
resolvent modes intensifies along a vertical band $\lambda_{x}/\delta\geq 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
($\lambda_{x}=15\lambda_{z}$, $10\lambda_{z}$ and $5\lambda_{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 $\lambda_{z}\gg\lambda_{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_{\tau}$ indicate high energy
in the vertical band as well [41].

###  III.2 Mode shapes

Table 2: Representative wavenumber combinations that we will explore in §III.2. Mode name | $k_{x}$ | $k_{z}$ | $c$
---|---|---|---
E1: most energetic mode for $y^{+}=15$ | $\pi/2$ | $4\pi$ | $\overline{u}(y^{+}=15)$
E2: most energetic mode for $y^{+}=30$ | $\pi/2$ | $3\pi$ | $\overline{u}(y^{+}=30)$
E3: most energetic mode for $y^{+}=100$ | $\pi/2$ | $2\pi$ | $\overline{u}(y^{+}=100)$

(a)  (b)

Figure 4: Amplitudes of the leading resolvent response modes for the (a)
streamwise velocity and (b) density for $Ri_{\tau}=0,10,18,60,100$ (darker to
lighter) at $c=\bar{u}(y^{+}=15)$ (dashed line), $\bar{u}(y^{+}=30)$ (dot-
dashed line) and $\bar{u}(y^{+}=100)$ (dotted line) for wave-parameters
corresponding to E1, E2 and E3, respectively. The subscripts $u$ and $\rho$
indicate the corresponding components of the resolvent response mode.

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_{\tau}$ 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 2. 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_{\tau}$
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.

(a)  (b)  (c)

Figure 5: Two-dimensional response mode shapes for
$(k_{x},k_{z})=(\pi/2,4\pi)$ at a critical-layer location of $y^{+}=15$ for
(a) $Ri_{\tau}=0$, (b) $18$, and (c) $100$. Red and blue contours represent
positive and negative fluctuations, respectively. The contour levels are
scaled by the maximum of each mode component. The dashed black line in each
sub-plot is the location of the critical-layer where
$c=\overline{u}(y^{+}=15)$.

(a)  (b)

Figure 6: Autocorrelation coefficients $C_{uu}$, $C_{vv}$, $C_{ww}$ and
$C_{\rho\rho}$ of the DNS at $y^{+}=15$ for (a) $Ri_{\tau}=0$ and (b) $100$.
Red and blue contours represent positive (0.4, 0.6, 0.8) and negative (-0.2)
correlation, respectively, with each contour level signifying 0.2 increments.
The horizontal dashed line is $y^{+}=30$ and the vertical dotted line is
$\Delta x=0$.

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 as

$\hat{R}_{qq}(k_{x},y,y^{\prime},k_{z})=\langle\hat{q}(k_{x},y,k_{z})\hat{q}^{*}(k_{x},y^{\prime},k_{z})\rangle,$
(15)

where $q$ is a generic variable of zero mean and $\langle\cdot\rangle$ is the
expected value. The auto-covariance in physical space, $R_{qq}(\Delta
x,y,y^{\prime},\Delta z)$, is obtained as the inverse Fourier transform of
$\hat{R}$, where $\Delta x=x-x^{\prime}$ and $\Delta z=z-z^{\prime}$ are the
distances between the two points in the homogeneous directions. The
autocorrelation coefficient,

$C_{qq}(\bm{x},\bm{x}^{\prime})=\frac{R_{qq}(\bm{x},\bm{x}^{\prime})}{\varsigma_{q}(\bm{x})\varsigma_{q}(\bm{x}^{\prime})},$
(16)

is obtained by normalising the covariance with the product of the standard
deviations, $\varsigma$, 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 Fig. 5. The autocorrelations of the
streamwise, wall-normal and spanwise velocity fields as well as the density
field are shown in Fig. 6, for a two-dimensional slice at $\Delta z=0$. The
reference location $y^{\prime+}=15$.

The LES of Armenio _et al._ [34] and the DNS of García-Villalba & del Álamo
[41] demonstrated that structures in the near-wall region are largely
unaffected by stable stratification. As expected, both the resolvent response
modes and the correlations do not change significantly for the range of
$Ri_{\tau}$ considered. For the velocities, the main difference is a reduction
in the autocorrelation coefficient in the stratified case. The largest
difference occurs for density properties as the phase in the wall-normal
direction along the resolvent response modes are shifted, creating structures
that are more detached from the wall. Similarly, the density correlations are
wall-attached for $Ri_{\tau}=0$ whereas they are more detached for
$Ri_{\tau}=100$.

(a)  (b)  (c)

Figure 7: Two-dimensional response mode shapes for
$(k_{x},k_{z})=(\pi/2,3\pi)$ at a critical-layer location of $y^{+}=30$ for
(a) $Ri_{\tau}=0$, (b) $18$, and (c) $100$. Red and blue contours represent
positive and negative fluctuations, respectively. The contour levels are
scaled by the maximum of each mode component. The dashed black line in each
sub-plot is the location of the critical-layer where
$c=\overline{u}(y^{+}=30)$.

(a)  (b)

Figure 8:  Autocorrelation coefficients $C_{uu}$, $C_{vv}$, $C_{ww}$ and
$C_{\rho\rho}$ of the DNS at $y^{+}=30$ for (a) $Ri_{\tau}=0$ and (b) $100$.
Red and blue contours represent positive (0.4, 0.6, 0.8) and negative (-0.2)
correlation, respectively, with each contour level signifying 0.2 increments.
The horizontal dashed line is $y^{+}=30$ and the vertical dotted line is
$\Delta x=0$.

We plot the resolvent response modes for the wavenumbers and wavespeed
corresponding to E2 (Fig. 7) and the correlations for $y^{\prime+}=30$ at
$\Delta z=0$ (Fig. 8). The results are similar to that of E1, since the
velocity response modes do not vary across $Ri_{\tau}$, 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_{\tau}=0$ and
$Ri_{\tau}=100$ case, although the centre of the density correlation for the
$Ri_{\tau}=100$ case lies farther away from the wall.

(a)  (b)  (c)

Figure 9: Two-dimensional response mode shapes for
$(k_{x},k_{z})=(\pi/2,2\pi)$ at a critical-layer location of $y^{+}=100$ for
(a) $Ri_{\tau}=0$, (b) $18$, and (c) $100$. Red and blue contours represent
positive and negative fluctuations, respectively. The contour levels are
scaled by the maximum of each mode component. The dashed black line in each
sub-plot is the location of the critical-layer where
$c=\overline{u}(y^{+}=100)$.

(a)  (b)

Figure 10: Autocorrelation coefficients $C_{uu}$, $C_{vv}$, $C_{ww}$ and
$C_{\rho\rho}$ of the DNS at $y^{+}=100$ for (a) $Ri_{\tau}=0$ and (b) $100$.
Red and blue contours represent positive (0.4, 0.6, 0.8) and negative (-0.2)
correlation, respectively, with each contour level signifying 0.2 increments.
The horizontal dashed line is $y^{+}=30$ and the vertical dotted line is
$\Delta x=0$.

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^{\prime+}=100$ at
$\Delta 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_{\rho\rho}$ term
for $Ri_{\tau}=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.

###  III.3 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 [70, 50] as

$\displaystyle\mathcal{P}_{\text{tot}}(y)$
$\displaystyle=\mathbb{R}\left[-\frac{\partial\overline{u}}{\partial
y}\sum_{j}\int_{-\infty}^{\infty}\sigma_{j}^{2}\chi_{j}^{2}\Big{(}\bm{\hat{\psi}}^{*}_{j,u}\bm{\hat{\psi}}_{j,v}\Big{)}d\bm{k}\right],$
(17a) $\displaystyle\mathcal{T}_{\text{tot}}(y)$
$\displaystyle=\mathbb{R}\left[\sum_{j}\sum_{i}\int_{-\infty}^{\infty}\sigma_{j}\chi_{j}\chi_{i}D_{y}\Big{(}\bm{\hat{\phi}}^{*}_{i,u}\bm{\hat{\psi}}_{j,v}+\bm{\hat{\phi}}^{*}_{i,v}\bm{\hat{\psi}}_{j,v}+\bm{\hat{\phi}}^{*}_{i,w}\bm{\hat{\psi}}_{j,v}\Big{)}d\bm{k}\right],$
(17b) $\displaystyle\mathcal{B}_{\text{tot}}(y)$
$\displaystyle=\mathbb{R}\left[-Ri_{\tau}\sum_{j}\int_{-\infty}^{\infty}\sigma_{j}^{2}\chi_{j}^{2}\Big{(}\bm{\hat{\psi}}^{*}_{j,v}\bm{\hat{\psi}}_{j,\rho}\Big{)}d\bm{k}\right],$
(17c) $\displaystyle\mathcal{V}_{\text{tot}}(y)$
$\displaystyle=\mathbb{R}\left[\frac{1}{Re_{\tau}}\sum_{j}\int_{-\infty}^{\infty}\sigma_{j}^{2}\chi_{j}^{2}\Big{(}\bm{\hat{\psi}}^{*}_{j,u}\hat{\Delta}\bm{\hat{\psi}}_{j,u}+\bm{\hat{\psi}}^{*}_{j,v}\hat{\Delta}\bm{\hat{\psi}}_{j,v}+\bm{\hat{\psi}}^{*}_{j,w}\hat{\Delta}\bm{\hat{\psi}}_{j,w}\Big{)}d\bm{k}\right],$
(17d)

where $\chi_{j}$, $\sigma_{j}$, $\bm{\hat{\psi}}_{j}$ and
$\bm{\hat{\phi}}_{j}$ are functions of $\bm{k}$ and the subscript $u,v,w,\rho$
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
$\bm{k}$, defined as

$\displaystyle\mathcal{P}(y,\bm{k})$
$\displaystyle=\mathbb{R}\left[-\frac{\partial\overline{u}}{\partial
y}\sigma_{1}^{2}\Big{(}\bm{\hat{\psi}}^{*}_{1,u}\bm{\hat{\psi}}_{1,v}\Big{)}\right],$
(18a) $\displaystyle\mathcal{T}(y,\bm{k})$
$\displaystyle=\mathbb{R}\left[\sigma_{1}D_{y}\Big{(}\bm{\hat{\phi}}^{*}_{1,u}\bm{\hat{\psi}}_{1,v}+\bm{\hat{\phi}}^{*}_{1,v}\bm{\hat{\psi}}_{1,v}+\bm{\hat{\phi}}^{*}_{1,w}\bm{\hat{\psi}}_{1,v}\Big{)}\right],$
(18b) $\displaystyle\mathcal{B}(y,\bm{k})$
$\displaystyle=\mathbb{R}\left[-Ri_{\tau}\sigma_{1}^{2}\Big{(}\bm{\hat{\psi}}^{*}_{1,v}\bm{\hat{\psi}}_{1,\rho}\Big{)}\right],$
(18c) $\displaystyle\mathcal{V}(y,\bm{k})$
$\displaystyle=\mathbb{R}\left[\frac{1}{Re_{\tau}}\sigma_{1}^{2}\Big{(}\bm{\hat{\psi}}^{*}_{1,u}\hat{\Delta}\bm{\hat{\psi}}_{1,u}+\bm{\hat{\psi}}^{*}_{1,v}\hat{\Delta}\bm{\hat{\psi}}_{1,v}+\bm{\hat{\psi}}^{*}_{1,w}\hat{\Delta}\bm{\hat{\psi}}_{1,w}\Big{)}\right].$
(18d)

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.

Figure 11: Energy budget terms computed from resolvent modes Eq. (18a–18d) for
wavenumbers given by (a) E1 at $c=\overline{u}(y^{+}=15)$, (b) E2 at
$c=\overline{u}(y^{+}=30)$, and (c) E3 at $c=\overline{u}(y^{+}=100)$ for
$Ri_{\tau}=0,10,18,60,100$ (darker to lighter). Figure 12: (Energy budget
terms computed from the DNS for $Ri_{\tau}=0,10,18,60,100$ (darker to
lighter).

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_{\tau}$ 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_{\tau}$ increases farther away from the wall, and at
$y^{+}=100$, the production (and viscous diffusion) of the $Ri_{\tau}=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
$\mathcal{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 DNS energy budget term as in Ref. [71], but
this would require a time-resolved dataset for a longer time domain. The
resolvent turbulent transport stays relatively similar among different
$Ri_{\tau}$, as does the turbulent transport term from DNS. The buoyancy flux
is much more dependent on $Ri_{\tau}$, with variations becoming greater
farther away from the wall in both the DNS and resolvent results.

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.

(a)  (b)

Figure 13: (a) Resolvent energy budget terms for wavenumber combinations E1,
E2 and E3 evaluated at $y^{+}=15,30,100$, respectively, for
$Ri_{\tau}=0,10,18,60,100$ (darker to lighter) normalized with
$\mathcal{P}(y^{+}=15)$ for $Ri_{\tau}=0$ and wavenumber combination E1. (b)
DNS energy budget terms at $y^{+}=15,30,100$ for $Ri_{\tau}=0,10,18,60,100$
(darker to lighter) normalized with $\mathcal{P}_{\text{DNS}}(y^{+}=15)$ for
$Ri_{\tau}=0$. Symbols are $\mathcal{P}$, circles; $\mathcal{T}$, triangles;
$\mathcal{V}$, crosses; and $\mathcal{B}$, asterisks.

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.

##  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_{\tau}$
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_{\tau}$. 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.

## Acknowledgements

The support of a Vannevar Bush Faculty Fellowship administered under the U.S.
Office of Naval Research, grant #N00014-17-1-3022, is gratefully acknowledged.
Additionally, the authors would like to thank Dr. Angeliki Laskari for
insightful discussions.

## References

  * McKeon and Sharma [2010] B. J. McKeon and A. S. Sharma, A critical-layer framework for turbulent pipe flow, J. Fluid Mech. 658, 336 (2010).
  * Dawson  _et al._ [2018] S. T. Dawson, T. Saxton-Fox, and B. J. McKeon, Modeling passive scalar dynamics in wall-bounded turbulence using resolvent analysis, in  _2018 AIAA Fluid Dynamics Conference_ (2018) p. 4042.
  * Nieuwstadt [1984] F. T. M. Nieuwstadt, The turbulent structure of the stable, nocturnal boundary layer, J. Atm. Sci. 41, 2202 (1984).
  * Stull [2000] R. B. Stull,  _Meteorology for scientists and engineers_ (Brooks/Cole, 2000).
  * Wunsch and Ferrari [2004] C. Wunsch and R. Ferrari, Vertical mixing, energy, and the general circulation of the oceans, Annu. Rev. Fluid Mech. 36, 281 (2004).
  * Thorpe [2005] S. A. Thorpe,  _The Turbulent Ocean_ (Cambridge University Press, 2005).
  * Panofsky and Dutton [1984] H. A. Panofsky and J. A. Dutton,  _Atmospheric Turbulence: Models and methods for engineering applications_ (Wiley, 1984).
  * Sorbjan [1989] Z. Sorbjan,  _Structure of the atmospheric boundary layer_ (Prentice Hal, 1989).
  * Stull [1988] R. B. Stull,  _An Introduction to Boundary Layer Meteorology_ , Vol. 13 (Springer Science & Business Media, 1988).
  * Wyngaard [2010] J. C. Wyngaard,  _Turbulence in the Atmosphere_ (Cambridge University Press, 2010).
  * Garratt [1994] J. R. Garratt, The atmospheric boundary layer, Earth-Sci. Rev. 37, 89 (1994).
  * Ivey  _et al._ [2008] G. N. Ivey, K. B. Winters, and J. R. Koseff, Density stratification, turbulence, but how much mixing?, Annu. Rev. Fluid Mech. 40, 169 (2008).
  * Mahrt [2014] L. Mahrt, Stably stratified atmospheric boundary layers, Annu. Rev. Fluid Mech. 46, 23 (2014).
  * Mahrt [1999] L. Mahrt, Stratified atmospheric boundary layers, Bound.-Layer Meteorol. 90, 375 (1999).
  * Williams  _et al._ [2017] O. Williams, T. Hohman, T. Van Buren, E. Bou-Zeid, and A. J. Smits, The effect of stable thermal stratification on turbulent boundary layer statistics, J. Fluid Mech. 812, 1039 (2017).
  * Chauhan  _et al._ [2013] K. Chauhan, N. Hutchins, J. Monty, and I. Marusic, Structure inclination angles in the convective atmospheric surface layer, Bound.-Layer Meteorol. 147, 41 (2013).
  * Salesky and Anderson [2018] S. T. Salesky and W. Anderson, Buoyancy effects on large-scale motions in convective atmospheric boundary layers: implications for modulation of near-wall processes, J. Fluid Mech. 856, 135 (2018).
  * Salesky and Anderson [2020] S. T. Salesky and W. Anderson, Revisiting inclination of large-scale motions in unstably stratified channel flow, J. Fluid Mech. 884 (2020).
  * Wyngaard and Coté [1971] J. C. Wyngaard and O. R. Coté, The budgets of turbulent kinetic energy and temperature variance in the atmospheric surface layer, J. Atm. Sci. 28, 190 (1971).
  * Kondo  _et al._ [1978] J. Kondo, O. Kanechika, and N. Yasuda, Heat and momentum transfers under strong stability in the atmospheric surface layer, J. Atm. Sci. 35, 1012 (1978).
  * Mahrt [1998] L. Mahrt, Nocturnal boundary-layer regimes, Bound.-Layer Meteorol. 88, 255 (1998).
  * Fernando and Weil [2010] H. J. S. Fernando and J. C. Weil, Whither the stable boundary layer? A shift in the research agenda, Bull. Am. Meteorol. Soc. 91, 1475 (2010).
  * Smedman  _et al._ [1994] A.-S. Smedman, M. Tjernström, and U. Högström, The near-neutral marine atmospheric boundary layer with no surface shearing stress: A case study, J. Atm. Sci. 51, 3399 (1994).
  * Stacey  _et al._ [1999] M. T. Stacey, S. G. Monismith, and J. R. Burau, Observations of turbulence in a partially stratified estuary, J. Phys. Oceanogr. 29, 1950 (1999).
  * Lu  _et al._ [2000] Y. Lu, R. G. Lueck, and D. Huang, Turbulence characteristics in a tidal channel, J. Phys. Oceanogr. 30, 855 (2000).
  * Large  _et al._ [1994] W. G. Large, J. C. McWilliams, and S. C. Doney, Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization, Rev. Geophys. 32, 363 (1994).
  * Arya [1975] S. P. S. Arya, Buoyancy effects in a horizontal flat-plate boundary layer, J. Fluid Mech. 68, 321 (1975).
  * Britter [1974] R. E. Britter,  _An experiment on turbulence in a density stratified fluid_ , Ph.D. thesis, Monash University, Australia (1974).
  * Piat and Hopfinger [1981] J. F. Piat and E. J. Hopfinger, A boundary layer topped by a density interface, J. Fluid Mech. 113, 411 (1981).
  * Komori  _et al._ [1983] S. Komori, H. Ueda, F. Ogino, and T. Mizushina, Turbulence structure in stably stratified open-channel flow, J. Fluid Mech. 130, 13 (1983).
  * Fukui  _et al._ [1983] K. Fukui, M. Nakajima, and H. Ueda, A laboratory experiment on momentum and heat transfer in the stratified surface layer, Q. J. R. Meteorol. Soc. 109, 661 (1983).
  * Miles [1961] J. W. Miles, On the stability of heterogeneous shear flows, J. Fluid Mech. 10, 496 (1961).
  * Garg  _et al._ [2000] R. P. Garg, J. H. Ferziger, S. G. Monismith, and J. R. Koseff, Stably stratified turbulent channel flows. I. Stratification regimes and turbulence suppression mechanism, Phys. Fluids 12, 2569 (2000).
  * Armenio and Sarkar [2002] V. Armenio and S. Sarkar, An investigation of stably stratified turbulent channel flow using large-eddy simulation, J. Fluid Mech. 459, 1 (2002).
  * Basu and Porté-Agel [2006] S. Basu and F. Porté-Agel, Large-eddy simulation of stably stratified atmospheric boundary layer turbulence: a scale-dependent dynamic modeling approach, J. Atmos. Sci. 63, 2074 (2006).
  * Stoll and Porté-Agel [2008] R. Stoll and F. Porté-Agel, Large-eddy simulation of the stable atmospheric boundary layer using dynamic models with different averaging schemes, Bound.-Layer Meteorol. 126, 1 (2008).
  * Iida  _et al._ [2002] O. Iida, N. Kasagi, and Y. Nagano, Direct numerical simulation of turbulent channel flow under stable density stratification, Int. J. Heat Mass Transf. 45, 1693 (2002).
  * Nieuwstadt [2005] F. T. M. Nieuwstadt, Direct numerical simulation of stable channel flow at large stability, Bound.-Layer Meteorol. 116, 277 (2005).
  * Brethouwer  _et al._ [2007] G. Brethouwer, P. Billant, E. Lindborg, and J.-M. Chomaz, Scaling analysis and simulation of strongly stratified turbulent flows, J. Fluid Mech. 585, 343–368 (2007).
  * Flores and Riley [2011] O. Flores and J. J. Riley, Analysis of turbulence collapse in the stably stratified surface layer using direct numerical simulation, Bound.-Layer Meteorol. 139, 241 (2011).
  * García-Villalba and Del Alamo [2011] M. García-Villalba and J. C. Del Alamo, Turbulence modification by stable stratification in channel flow, Phys. Fluids 23, 045104 (2011).
  * Yeh and Taira [2018] C.-A. Yeh and K. Taira, Resolvent-analysis-based design of airfoil separation control, J. Fluid Mech. 867, 572 (2018).
  * Towne  _et al._ [2018] A. Towne, O. T. Schmidt, and T. Colonius, Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis, J. Fluid Mech. 847, 821 (2018).
  * Bae  _et al._ [2020] H. J. Bae, S. T. Dawson, and B. J. McKeon, Resolvent-based study of compressibility effects on supersonic turbulent boundary layers, J. Fluid Mech. 883 (2020).
  * McMullen  _et al._ [2020] R. M. McMullen, K. Rosenberg, and B. J. McKeon, Interaction of forced Orr-Sommerfeld and Squire modes in a low-order representation of turbulent channel flow, Phys. Rev. Fluids 5, 084607 (2020).
  * Nogueira  _et al._ [2021] P. A. S. Nogueira, P. Morra, E. Martini, A. V. G. Cavalieri, and D. S. Henningson, Forcing statistics in resolvent analysis: application in minimal turbulent couette flow, J. Fluid Mech. 908, A32 (2021).
  * McKeon [2017] B. J. McKeon, The engine behind (wall) turbulence: perspectives on scale interactions, J. Fluid Mech. 817 (2017).
  * Jovanović [2020] M. R. Jovanović, From bypass transition to flow control and data-driven turbulence modeling: An input-output viewpoint, Annu. Rev. Fluid Mech. 53, null (2020).
  * Moarref  _et al._ [2013] R. Moarref, A. S. Sharma, J. A. Tropp, and B. J. McKeon, Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels, J. Fluid Mech. 734, 275 (2013).
  * Madhusudanan [2020] A. Madhusudanan,  _Coherent structures from the linearized Navier-Stokes equations for wall-bounded turbulent flows_ , Ph.D. thesis, University of Melbourne (2020).
  * Lorenz [1955] E. N. Lorenz, Available potential energy and the maintenance of the general circulation, Tellus 7, 157 (1955).
  * Turner [1979] J. S. Turner,  _Buoyancy Effects in Fluids_ (Cambridge University Press, 1979).
  * Orlandi [2000] P. Orlandi,  _Fluid Flow Phenomena: A Numerical Toolkit_ (Springer, 2000).
  * Wray [1990] A. A. Wray,  _Minimal-storage time advancement schemes for spectral methods_ , Tech. Rep. (NASA Ames Research Center, 1990).
  * Chorin [1968] A. J. Chorin, Numerical solution of the Navier-Stokes equations, Math. Comput. 22, 745 (1968).
  * Bae  _et al._ [2019] H. J. Bae, A. Lozano-Durán, S. T. Bose, and P. Moin, Dynamic slip wall model for large-eddy simulation, J. Fluid Mech. 859, 400 (2019).
  * Lozano-Durán and Bae [2019] A. Lozano-Durán and H. J. Bae, Characteristic scales of Townsend’s wall-attached eddies, J. Fluid Mech. 868, 698 (2019).
  * Ribeiro  _et al._ [2020] J. H. M. Ribeiro, C.-A. Yeh, and K. Taira, Randomized resolvent analysis, Phys. Rev. Fluids 5, 033902 (2020).
  * Symon  _et al._ [2018] S. Symon, K. Rosenberg, S. T. M. Dawson, and B. J. McKeon, Non-normality and classification of amplification mechanisms in stability and resolvent analysis, Phys. Rev. Fluids 3, 053902 (2018).
  * del Alamo  _et al._ [2004] J. C. del Alamo, J. Jiménez, P. Zandonade, and R. D. Moser, Scaling of the energy spectra of turbulent channels, J. Fluid Mech. 500, 135 (2004).
  * Hopfinger [1987] E. J. Hopfinger, Turbulence in stratified fluids: A review, J. Geophys. Research 92, 5287 (1987).
  * Tritton [1967] D. J. Tritton, Some new correlation measurements in a turbulent boundary layer, J. Fluid Mech. 28, 439 (1967).
  * Liu  _et al._ [2001] Z. Liu, R. J. Adrian, and T. J. Hanratty, Large-scale modes of turbulent channel flow: transport and structure, J. Fluid Mech. 448, 53 (2001).
  * Ganapathisubramani  _et al._ [2005] B. Ganapathisubramani, N. Hutchins, W. T. Hambleton, E. K. Longmire, and I. Marusic, Investigation of large-scale coherence in a turbulent boundary layer using two-point correlations, J. Fluid Mech. 524, 57 (2005).
  * Lee and Sung [2011] J. H. Lee and H. J. Sung, Very-large-scale motions in a turbulent boundary layer, J. Fluid Mech. 673, 80 (2011).
  * Pirozzoli and Bernardini [2011] S. Pirozzoli and M. Bernardini, Turbulence in supersonic boundary layers at moderate Reynolds number, J. Fluid Mech. 688, 120 (2011).
  * Sillero  _et al._ [2014] J. A. Sillero, J. Jiménez, and R. D. Moser, Two-point statistics for turbulent boundary layers and channels at Reynolds numbers up to $\delta^{+}\approx 2000$, Phys. Fluids 26, 105109 (2014).
  * Kline  _et al._ [1967] S. J. Kline, W. C. Reynolds, F. A. Schraub, and P. W. Runstadler, The structure of turbulent boundary layers, J. Fluid Mech. 30, 741 (1967).
  * Smith and Metzler [1983] C. R. Smith and S. P. Metzler, The characteristics of low-speed streaks in the near-wall region of a turbulent boundary layer, J. Fluid Mech. 129, 27 (1983).
  * Symon  _et al._ [2020] S. Symon, S. J. Illingworth, and I. Marusic, Energy transfer in turbulent channel flows and implications for resolvent modelling, J. Fluid Mech. (to appear), arXiv:2004.13266 (2020).
  * Mizuno [2016] Y. Mizuno, Spectra of energy transport in turbulent channel flows for moderate Reynolds numbers, J. Fluid Mech. 805, 171 (2016).
  * Moarref  _et al._ [2014] R. Moarref, M. R. Jovanović, J. A. Tropp, A. S. Sharma, and B. J. McKeon, A low-order decomposition of turbulent channel flow via resolvent analysis and convex optimization, Phys. Fluids 26, 051701 (2014).
  * Rosenberg and McKeon [2019] K. Rosenberg and B. J. McKeon, Efficient representation of exact coherent states of the Navier-Stokes equations using resolvent analysis, Fluid Dyn. Res. 51 (2019).