1. Introduction
Wall turbulence regulates surface exchange rates of momentum, heat, mass and other quantities. These exchanges are vitally important to flows in geophysics, affecting emissions of dust, water vapour, greenhouse gases, etc. In addition, external and internal flows in engineering applications affect aerodynamic profiles of lifting bodies and thermal efficiency of vapour power systems. For the case of inertia-dominated, neutrally stratified wall flows, transport of Reynolds-averaged vorticity, 
$\langle {\boldsymbol \omega } \rangle _t = \boldsymbol {\nabla } \times \langle {\boldsymbol u} \rangle _t$, where 
${\boldsymbol u}$ is velocity and 
$\langle \cdots \rangle _t$ denotes time averaging, is sustained via the so-called stretching and tilting and ‘turbulent torque’ mechanisms (Bradshaw (Reference Bradshaw1987), and references therein):
\begin{equation} \langle {\boldsymbol u} \rangle_t \boldsymbol{\cdot} \boldsymbol{\nabla} \langle {\boldsymbol \omega} \rangle_t = \underbrace{{\langle {\boldsymbol \omega} \rangle_t \boldsymbol{\cdot}\boldsymbol{\nabla} \langle {\boldsymbol u} \rangle_t}}_{Stretching\ and\ tilting} - \quad \underbrace{\boldsymbol{\nabla} \times \left(\boldsymbol{\nabla}\boldsymbol{\cdot} \langle \boldsymbol{\mathsf{T}} \rangle_t \right)}_{Turbulent\ torque}, \end{equation}
where 
$\boldsymbol{\mathsf{T}} = {\boldsymbol u}^\prime \otimes {\boldsymbol u}^\prime$ is the Reynolds (turbulent) stress tensor. Note that in this article, the first, second and third components of any vector correspond to its constituent magnitudes in the streamwise (
$x$), spanwise (
$y$) and vertical (
$z$) directions, respectively, where velocity is denoted via 
${\boldsymbol u} = u \hat {i}_x + v \hat {i}_y + w \hat {i}_z$. The stretching and tilting and turbulent torque terms are mechanistically responsible for sustenance of Prandtl's secondary flow of the first and second kinds, respectively. The streamwise component of (1.1) for adiabatic channel flows reduces to
\begin{equation} \langle u_i \rangle_t \partial_i \langle {\omega}_x \rangle_t ={-} \epsilon_{xqi} \partial_q \partial_j \langle \textsf{{T}}_{ji} \rangle_t, \end{equation}
where the Reynolds-averaged stretching and tilting term vanishes for streamwise homogeneity; early contributions in this area were provided by those studying Reynolds-averaged secondary flows in ducts (Prandtl Reference Prandtl1952; Hoagland Reference Hoagland1960; Brundrett & Baines Reference Brundrett and Baines1964; Perkins Reference Perkins1970; Gessner Reference Gessner1973). More recently, there has been significant interest in adiabatic, isothermal turbulent wall flows responding to spanwise heterogeneity in surface texture, including surface hydrophobicity (Jelly, Jung & Zaki Reference Jelly, Jung and Zaki2014; Lee, Jelly & Zaki Reference Lee, Jelly and Zaki2015), hydraulic flows responding to spanwise variation in bed gravel roughness (Wang & Cheng Reference Wang and Cheng2005; Vermaas, Uijttewall & Hoitink Reference Vermaas, Uijttewall and Hoitink2011), complex roughness with predominant spanwise heterogeneity (Barros & Christensen Reference Barros and Christensen2014; Pathikonda & Christensen Reference Pathikonda and Christensen2017), and canonical arrangements (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Kevin et al. Reference Kevin, Monty, Bai, Pathikonda, Nugroho, Barros, Christensen and Hutchins2017; Yang & Anderson Reference Yang and Anderson2017; Chung, Monty & Hutchins Reference Chung, Monty and Hutchins2018; Hwang & Lee Reference Hwang and Lee2018; Medjnoun, Vanderwel & Ganapathisubramani Reference Medjnoun, Vanderwel and Ganapathisubramani2018; de Silva et al. Reference de Silva, Kevin, Baidya, Hutchins and Marusic2018; Anderson Reference Anderson2019; Kevin, Monty & Hutchins Reference Kevin, Monty and Hutchins2019; Vanderwel et al. Reference Vanderwel, Stroh, Kriegseis, Frohnapfel and Ganapathisubramani2019; Stroh et al. Reference Stroh, Schäfer, Frohnapfel and Forooghi2020; Wangsawijaya et al. Reference Wangsawijaya, Baidya, Chung, Marusic and Hutchins2020; Zampiron, Cameron & Nikora Reference Zampiron, Cameron and Nikora2020, Reference Zampiron, Cameron and Nikora2021; Schafer et al. Reference Schafer, Stroh, Forooghi and Frohnapfel2022). In the context of spanwise roughness heterogeneity, Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015) demonstrated that the resultant roughness-driven secondary flows are, indeed, a manifestation of Prandtl's secondary flow of the second kind. This was accomplished via rough surfaces composed of ‘strips’ of relatively high and low aerodynamic roughness, 
$z_0^H$ and 
$z_0^L$, respectively, as shown in figure 1(a).
Figure 1. Schematic illustrating spatial heterogeneity in surface conditions due to spanwise variation in aerodynamic roughness length (a) and heat flux (b), where the former has been considered previously (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015). Flux heterogeneity caused by strips of width 
$\ell _H$, composed of relatively higher roughness length 
$z_0^H$ (a), or heat flux 
$Q_0^H$ (b), adjacent to strips of width 
$\ell _L$, composed of relatively smaller roughness length 
$z_0^L$ (a), or heat flux 
$Q_0^L$ (b). For this article, the heat flux ratio 
$\beta = Q_0^H/Q_0^L$ was prescribed to exceed unity for all cases and was varied systematically; the area-averaged heat flux, 
$Q_0 = \int _{\mathrm {d}^2 {\boldsymbol x}} Q_0 ({\boldsymbol x})\,\mathrm {d}^2 {\boldsymbol x} / \int _{\mathrm {d}^2 {\boldsymbol x}}\mathrm {d}^2 {\boldsymbol x}$, is constant for all spanwise-heterogeneous cases and equivalent to the homogeneous case. Problem parameters addressed in this study are summarized in table 1.
Table 1. Properties of LES cases, including: surface heat fluxes of high and low heat flux strips, 
$Q_0^H$ and 
$Q_0^L$, respectively; ratio of heat fluxes, 
$\beta = Q_0^H/Q_0^L$; width of high and low heat flux strips, 
$\ell _H$ and 
$\ell _L$, respectively (see also figure 1b); friction velocity, 
$u_\tau$; Deardorff convective velocity scale for aggregate surface, 
$w_\star$; boundary layer depth (
$z_i$); Obukhov length magnitude (
$|L|$); global stability parameter, 
$-z_i/L$; and heterogeneity parameter, 
$\mathcal {H}$. The mean pressure gradient force 
$\rho _0^{-1} \partial _x P = 5.0 \times 10^{-4}$ m s
$^{-2}$ and the horizontally-averaged surface heat flux 
$\langle Q_0 \rangle _{xy} = 0.16$ K m s
$^{-1}$ were identical for all cases.
Beyond spatial heterogeneity in aerodynamic/hydrodynamic surface roughness, one can envision a multitude of practically-important turbulent wall flows affected by spatial heterogeneity in surface thermal conditions. Foremost engineering and geophysics examples include heat exchangers and landscape heterogeneities in sensible heat flux, respectively. Recent review articles have provided literature benchmarks in the area of land–atmospheric interactions and spatial heterogeneity in thermal conditions (Bou-Zeid et al. (Reference Bou-Zeid, Anderson, Katul and Mahrt2020), Stoll et al. (Reference Stoll, Gibbs, Salesky, Anderson and Calaf2020) and references therein), while yet more recent work based upon extensive field measurements has provided additional insights on the surface flux-driven mechanisms responsible for flow response (Margairaz, Pardyjak & Calaf Reference Margairaz, Pardyjak and Calaf2020a,Reference Margairaz, Pardyjak and Calafb; Morrison et al. Reference Morrison, Calaf, Higgins, Drake, Perelet and Pardyjak2021). Other recent work has provided new insight into the spatial and dynamical nature of unstably stratified inertia-dominated channel turbulence (Salesky & Anderson Reference Salesky and Anderson2018, Reference Salesky and Anderson2019, Reference Salesky and Anderson2020). However, no prior study has specifically addressed canonical spanwise variability in heat flux, illustrated in figure 1(b). This panel shows spanwise-aligned strips of widths 
$\ell _H$ and 
$\ell _L$, where relatively high and low heat fluxes, 
$Q_0^H$ and 
$Q_0^L$, are prescribed over the former and latter, respectively. In this sense, the present cases represent a ‘thermal analogue’ of the preceding spanwise roughness cases illustrated by figure 1(a), but the resulting unstable flow physics result in a distinctly different secondary flow response.
It is important to distinguish between turbulent secondary flows that form in response to thermally heterogeneous surface conditions, and horizontal convective rolls (HCRs) (Lemone Reference Lemone1973; Moeng & Sullivan Reference Moeng and Sullivan1994; Khanna & Brasseur Reference Khanna and Brasseur1998; Salesky, Chamecki & Bou-Zeid Reference Salesky, Chamecki and Bou-Zeid2017) that form in the convective atmospheric boundary layer over thermally homogeneous surfaces when the surface heat flux is small and mean wind shear is strong. Reviews of HCR structure and formation mechanisms can be found elsewhere (Young et al. Reference Young, Kristovich, Hjelmfelt and Foster2002; Salesky et al. Reference Salesky, Chamecki and Bou-Zeid2017). While HCRs also have horizontal vorticity that is generated by thermal torque, they are typically aligned 10
$^\circ$–20
$^\circ$ to the left of the geostrophic wind in the northern hemisphere, and advect due to the mean wind. While HCRs and other quasi-linear structures (Young et al. Reference Young, Kristovich, Hjelmfelt and Foster2002) are of interest for convective boundary layer dynamics, we focus herein on the mechanisms responsible for the formation of thermal secondary flows that arise in response to thermally heterogeneous surface conditions and therefore are ‘locked in’ to the heterogeneity. Furthermore, in this work we establish the link between the surface-heterogeneity-induced circulations and their corresponding mathematical representation as a Prandtl's secondary flow of the third kind.
The presence of unstable stratification within a channel flow is prescribed mathematically via a buoyancy force, 
${\boldsymbol f}_{B}$, in the momentum transport equation
\begin{equation} {\boldsymbol f}_{B} = {\boldsymbol e}_{x} 0 + {\boldsymbol e}_{y} 0 + {\boldsymbol e}_{z} \frac{g}{\theta_0}\,\delta \theta, \end{equation}
where 
${\boldsymbol e}_x = \{1, 0, 0 \}$, 
${\boldsymbol e}_y = \{0, 1, 0 \}$ and 
${\boldsymbol e}_z = \{0, 0, 1 \}$ are the standard basis functions, 
$g$ is acceleration due to gravity, 
$\theta _0$ is a reference potential temperature, and 
$\delta \theta$ is thermal variation, defined as 
$\delta \theta = (\theta - \theta _0)$. Application of the curl operation to the buoyancy force yields a thermal torque, 
${\boldsymbol f}_{B}^{\boldsymbol \omega }$, within (1.1) and (1.2):
\begin{equation} {\boldsymbol f}_{B}^{\boldsymbol \omega} = \boldsymbol{\nabla} \times {\boldsymbol f}_{B} = {\boldsymbol e}_{x} \frac{g}{\theta_0} \partial_y\,\delta \theta + {\boldsymbol e}_{y} 0 + {\boldsymbol e}_{z} 0. \end{equation}With this, the streamwise component of the thermal torque vector (1.4) appears within a revised equation (1.2):
\begin{equation} \underbrace{ \langle u_i \rangle_t \partial_i \langle {\omega}_x \rangle_t }_{\varOmega_{A}} = \underbrace{ - \epsilon_{xqi} \partial_q \partial_j \langle \textsf{{T}}_{ji} \rangle_t }_{\varOmega_{T}} + \underbrace{ \frac{g}{\theta_0} \partial_y\,\delta \theta }_{\varOmega_{B}}, \end{equation}
where underbraces are used to introduce symbols for advection 
$\varOmega _{A}$, turbulent torque 
$\varOmega _{T}$, and thermal torque 
$\varOmega _{B}$. In the absence of spanwise-variable thermal gradients – i.e. 
$\delta \theta (y,z) \rightarrow \delta \theta (z)$ – the second right-hand-side term vanishes, 
$\varOmega _{B} \rightarrow 0$. However, via prescribed spanwise-variable heat flux, the second right-hand-side term provides sustenance for the balance of Reynolds-averaged streamwise vorticity. And although this mechanism acts in conjunction with the turbulent torque term, its distinctly different origins render the term exclusive, herein defined as Prandtl's secondary flow of the third kind.
This third kind transport can, alternatively, be appreciated via transport of Reynolds-averaged turbulent kinetic energy, 
$\langle k \rangle _t = (1/2) \langle {\boldsymbol u}^\prime : {\boldsymbol u}^\prime \rangle _t$:
\begin{equation} 0 = \mathcal{A}_k + \mathcal{P}_{S} + \mathcal{P}_{B} - \epsilon, \end{equation}
where the first right-hand-side term is advection of 
$\langle k \rangle _t$, 
$\mathcal {A}_k = - \langle {\boldsymbol u} \rangle _t \boldsymbol {\cdot } \boldsymbol {\nabla } \langle k \rangle _t$, the second and third right-hand-side terms represent production via shear and buoyancy, 
$\mathcal {P}_{S} = - \langle {\boldsymbol u}^\prime \otimes {\boldsymbol u}^\prime \rangle _t : \boldsymbol {\nabla } \langle {\boldsymbol u} \rangle _t$ and 
$\mathcal {P}_{B} = - ({g}/{\theta _0}) \langle w^\prime \theta ^\prime \rangle _t$, respectively, and the fourth right-hand-side term is dissipation, 
$\epsilon = \langle - \boldsymbol{\mathsf{T}} : \boldsymbol{\mathsf{S}} \rangle _{xt}$. Note that under the Boussinesq (eddy viscosity) modelling framework, the deviatoric stresses are proportional to the Reynolds-averaged strain-rate tensor, 
$\boldsymbol{\mathsf{T}} = - \nu _t \boldsymbol{\mathsf{S}}$, where 
$\nu _t$ is eddy viscosity. With this, 
$\epsilon = \nu _t \langle \boldsymbol{\mathsf{S}} : \boldsymbol{\mathsf{S}} \rangle _{xt}$, illustrating that dissipation is positive definite, such that negative 
$\epsilon$ – as it appears in (1.6) – serves to deplete 
$\langle k \rangle _t$. For the ‘fully rough’ or inertial-dominated conditions assessed here, 
$\langle k \rangle _t$ transport due to pressure fluctuations, turbulent fluctuations and diffusion is vanishingly small and subsequently omitted; this is typical of high-Reynolds-number flows of practical importance in engineering and geophysics (Pope Reference Pope2000). Inspection of (1.6) highlights readily how any relative production–dissipation imbalance necessitates the existence of a finite advective term (even in the Reynolds-averaged limit). In subsequent developments, we will demonstrate how thermal torque (1.5) and buoyancy-driven production of turbulence (turbulent torque) affect resultant Reynolds-averaged secondary flows in flow over surfaces with spanwise-variable prescribed heat flux (i.e. figure 1b).
1.1. Dimensional analysis of the vorticity equation
Prior to investigating the formation and attributes of thermally-driven secondary circulations across a range of relevant parameters using large-eddy simulation, it is instructive to first leverage dimensional arguments of the reduced streamwise vorticity equation (1.5). As demonstrated in the following developments, these rudimentary arguments yield robust guidance on the role of thermal torque in the balance of streamwise vorticity in turbulent flows responding to surfaces with spanwise thermal variability. These scaling arguments offer preliminary evidence that thermal torque is a distinct driver of streamwise vorticity, which herein is defined as Prandtl's secondary flow of the third kind.
We first normalize (1.5) by the mean advection of vorticity, and rewrite (1.5) as a simplified balance between two non-dimensional ratios: production of vorticity by turbulent torque over advection (
$\varOmega _{T}/\varOmega _{A}$), and production of vorticity by thermal torque over advection (
$\varOmega _{B}/\varOmega _{A}$). This is
\begin{equation} \frac{\varOmega_{T}}{\varOmega_{A}} = 1- \frac{\varOmega_{B}}{\varOmega_{A}}. \end{equation}
By adopting this simplified non-dimensional representation of (1.5), the streamwise-averaged vorticity equation can be described advantageously by a single non-dimensional term, which we choose as 
$\varOmega _{B}/\varOmega _{A}$. In this case, when 
$\varOmega _{B}/\varOmega _{A} = 0$, turbulent torque is solely responsible for the existence of any advection. Alternatively, when 
$\varOmega _{B}/\varOmega _{A} = 1$, it is indicative that either the turbulent torque is zero, or it is very small in comparison to the advective transport since 
$\varOmega _{T}/\varOmega _{A} = 0$. For values of 
$\varOmega _{B}/\varOmega _{A} > 1$, 
$\varOmega _{T}/\varOmega _{A}$ becomes negative, indicating that the turbulent torque depletes Reynolds-averaged streamwise vorticity, or that it is balanced by thermal torque.
Because we have selected the ratio 
$\varOmega _{B}/\varOmega _{A}$ to be the relevant non-dimensional term, we proceed by developing a scaling analysis for this term only; the implications of variation in this term have been established in the preceding developments (equation (1.7) and accompanying text). The ratio 
$\varOmega _{B}/\varOmega _{A}$ is considered to be relatively more instructive, compared to 
$\varOmega _{T}/\varOmega _{A}$, since it precludes the need for a priori definition of characteristic scales for the range of different turbulent stress terms within the turbulent torque term (
$\epsilon _{xqi} \partial _q \partial _j \langle \textsf{{T}}_{ji} \rangle _t = (\partial _z^2 - \partial _y^2) \langle \textsf{{T}}_{yz} \rangle _t + \partial _{yz} (\langle \textsf{{T}}_{yy} \rangle _t - \langle \textsf{{T}}_{zz} \rangle _t$).
For the flows considered herein, mean spanwise velocity scales as surface shear velocity, 
$\langle v \rangle _t \sim u_\tau$, while mean vertical velocity scales with the Deardorff convective velocity scale, 
$\langle w \rangle _t \sim w_\star$, where 
$w_\star = (g Q_0 z_i / \theta _0 )^{1/3}$ (
$g$ is acceleration due to gravity, 
$Q_0$ is surface heat flux, 
$z_i$ is vertical extent of the flow, and 
$\theta _0$ is surface temperature). The ‘
$\sim$’ symbol is used here to mean ‘scales on’, and should not be interpreted as an exact equality. Similarly, spanwise position can be normalized by the spanwise heat flux heterogeneity extent, 
$y \sim \ell _H$. These arguments can be used within the Reynolds-averaged continuity equation to recover 
$\ell _d$, a scale for the vertical extent of thermally-driven secondary cells:
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \langle {\boldsymbol u} \rangle_t = \partial_y \langle v \rangle_t + \partial_z \langle w \rangle_t \rightarrow \ell_d \sim \ell_H \frac{w_\star}{u_\tau}, \end{equation}
where streamwise variability vanishes for the thermally-unstable channels considered herein. Note that 
$\ell _d \leq z_i$, that is, the depth of the secondary flow may not exceed the boundary layer depth. Using these scaling variables, and the (1.8) result, we recover a subsequent normalizing scale for streamwise vorticity:
\begin{equation} \langle \omega_x \rangle_t = \epsilon_{xjk} \partial_j \langle u_k \rangle_t \sim \left(\frac{w_\star}{\ell_H}, \frac{u_\tau}{\ell_d}\right) \sim \left(\frac{w_\star}{\ell_H}, \frac{w_\star}{\ell_H} \left[\frac{u_\tau}{w_\star} \right]^2 \right) \sim \frac{w_\star}{\ell_H}, \end{equation}
since 
$u_\tau / w_\star \sim {O}(1)$, which supports the final deduction leading to (1.9) (discussion of simulation attributes to follow). Note that 
$u_\tau /w_\star \sim {O}(1)$ for the parameter space that we consider in this work; however, the scaling arguments hold true more generally, as long as 
$w_\star$ remains finite. With these developments, we recover a scaling argument for the advective term:
\begin{equation} \varOmega_{A} = \langle v \rangle_t \partial_y \langle {\omega}_x \rangle_t + \langle w \rangle_t \partial_z \langle {\omega}_x \rangle_t \sim \left(\frac{u_\tau w_\star}{\ell_H^2}, \frac{w_\star^2}{\ell_H \ell_d} \right) \sim \frac{u_\tau w_\star}{\ell_H^2}. \end{equation}
Similarly, we recover a scaling argument for thermal torque, 
$\varOmega _{B} \sim (g / \ell _H) (\delta \varTheta / \theta _0)$, where 
$\delta \varTheta$ is a scale for the near-wall spanwise thermal difference while 
$\theta _0$ is a representative surface temperature. With this, we subsequently recover scaling arguments for the aforementioned ratio of thermal torque to vorticity advection term:
\begin{equation} \frac{\varOmega_{B}}{\varOmega_{A}} \sim \left(\frac{g \ell_H}{u_\tau^2}\, \frac{\delta \varTheta}{\theta_0}\right) \left(\frac{u_\tau}{w_\star}\right) \sim \mathcal{H} \frac{u_\tau}{w_\star}, \end{equation}
Equation (1.11) indicates that the normalized contribution of the thermal torque can be characterized by the product of two main non-dimensional terms. Readers will recognize 
$\mathcal {H} = ({g \ell _H}/{u_\tau ^2}) ( {\delta \varTheta }/{\theta _0})$ as a modified Richardson number (Stull Reference Stull1988), which quantifies the capacity for thermal gradients to sustain mixing in the spanwise–wall-normal plane. Such a deduction is thoroughly consistent with the arguments proposed herein: spanwise thermal gradients sustain mean secondary cells, and this mechanism for secondary flow sustenance differs distinctly from turbulent torque. The second non-dimensional term, 
${u_\tau }/{w_\star }$, quantifies the relative efficiency of shear- and thermally-driven mixing. Using the definitions of the Deardorff convective velocity scale and the Obukhov length
\begin{equation} L ={-} \frac{u_\tau^3 \varTheta_0}{\kappa g Q_0}, \end{equation}
where 
$\kappa$ is the von Kármán constant, one can write
\begin{equation} \frac{u_\tau}{w_\star} = \kappa^{1/3} \left({-}z_i/L \right)^{{-}1/3}; \end{equation}
thus one could express the velocity ratio in (1.11) in terms of the global stability parameter 
$-z_i/L$. Note that 
$-z_i/L \rightarrow 0$ for neutral stratification, and 
$-z_i/L \rightarrow +\infty$ in the free-convective limit. In the present work, we consider cases where 
$-z_i/L \approx 4$ and vary the properties of the surface thermal heterogeneity.
It is stressed, also, that Margairaz et al. (Reference Margairaz, Pardyjak and Calaf2020b) recently demonstrated how 
$\mathcal {H}$ serves to determine when mean horizontal transport of heat is large in comparison to turbulent transport when the flow is perturbed by randomly-distributed surface thermal heterogeneities. In their work this result was used to determine the conditions when horizontal transport of heat cannot be neglected in the study of atmospheric boundary layer flows. Similarly, the velocity ratio 
${u_\tau }/{w_\star }$ was recently identified in De Roo et al. (Reference De Roo, Zhang, Huq and Mauder2018) as a scaling variable determining the relative contribution of large secondary circulations in the imbalance for the closure of the surface energy. To this extent, the arguments leading to (1.11) exhibit underlying consistency with prior contributions in complementary areas.
1.2. This study
The scaling arguments outlined in § 1.1 provide direct guidance regarding the underlying importance of spanwise thermal heterogeneity (Margairaz et al. Reference Margairaz, Pardyjak and Calaf2020a,Reference Margairaz, Pardyjak and Calafb). Large-eddy simulation (LES) has been used to model flow over a series of surfaces with a priori prescribed spanwise-variable heat flux; a schematic of the surface configuration can be found in figure 1(b). We have systematically varied two problem parameters: heat flux ratio, 
$\beta = Q_0^H/Q_0^L$, and strip width ratio, 
$\lambda = \ell _H/\ell _L$, while for all cases the aggregate heat flux is held constant, 
$\langle Q_0 \rangle _{xy} = \int _{\mathrm {d}^2 {\boldsymbol x}} Q_0 ({\boldsymbol x}) \,\mathrm {d}^2 {\boldsymbol x} / \int _{\mathrm {d}^2 {\boldsymbol x}} \,\mathrm {d}^2 {\boldsymbol x}$. Specifically, we impose 
$\langle Q_0 \rangle _{xy} = 0.16\,{\rm K}\,{\rm m}\,{\rm s}^{-1}$. For the value of pressure gradient force selected (
$\rho_0^{-1} \partial_x P = 5.0 \times 10^{-4}\ \textrm{m s}^{-2}$), the global stability parameter is 
$-z_i/L \approx 4$, which corresponds to sheared convection. Our main focus is on how secondary flow properties vary with the heat flux ratio 
$\beta$, the strip width ratio 
$\lambda$, and the product of the modified Richardson number and the velocity scale ratio 
$\mathcal {H} u_\tau / w_\star$. While variability in 
$-z_i/L$ is also of interest, this would expand the parameter space significantly beyond what we are able to address in the present article. We show instantaneous and Reynolds-averaged first- and second-order turbulence statistics, in addition to constituent terms regulating transport of 
$\langle k \rangle _t$, all of which point consistently to the ability of spanwise heat flux variability to sustain Reynolds-averaged secondary flows (as per (1.4) and (1.5)). Turbulence statistics are used within the transport equation for Reynolds-averaged 
$k$, which provides support for the observed secondary flow direction. The results are cast against prior results for flow over spanwise-heterogeneous surface roughness (Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015), which provides context for the relative differences in secondary flow response due to thermal heterogeneity, and supports defining the flow response as Prandtl's secondary flow of the third kind. Additional turbulence statistics are used to demonstrate monotonic trends in secondary flow intensity with variation in salient forcing parameters.
In § 2, we outline briefly the LES code, accompanying numerical procedures, and details of the considered cases. Results are presented in § 3, which culminates in demonstrating that thermal torque is a distinctly different secondary flow driver and thus produces Prandtl's secondary flow of the third kind; this result is summarized via a concluding schematic. Concluding remarks and perspectives are presented in § 4.
2. Numerical procedure and cases
2.1. Large-eddy simulation
The LES code employed in the present work, described in Albertson & Parlange (Reference Albertson and Parlange1999) and Kumar et al. (Reference Kumar, Kleissl, Meneveau and Parlange2006), solves the three-dimensional filtered momentum and potential temperature equations written in rotational form. Horizontal spatial derivatives are calculated in Fourier space (pseudo-spectral projection), while vertical derivatives are calculated via second-order centred finite differencing. Time integration is performed through the fully explicit second-order Adams–Bashforth method. Nonlinear terms are fully dealiased, following the 3/2 rule (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Thomas2012). The subgrid-scale (SGS) model for momentum uses the Lagrangian scale-dependent (LASD) dynamic model (Bou-Zeid, Meneveau & Parlange Reference Bou-Zeid, Meneveau and Parlange2005), where the dynamic procedure (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991) is applied by averaging over Lagrangian trajectories of fluid parcels (Meneveau, Lund & Cabot Reference Meneveau, Lund and Cabot1996) to determine the optimal value of the Smagorinsky coefficient. The SGS heat flux is modelled using a constant SGS Prandtl number model: 
${\boldsymbol q}^{sgs} = -\nu _t^{sgs}/{Pr}^{sgs} \boldsymbol {\nabla } \tilde {\theta }$, where 
$\nu _t^{sgs} = (c_s \varDelta )^2 |\tilde {\boldsymbol{\mathsf{S}}}|$ is the SGS eddy viscosity, 
$\varDelta = (\varDelta _x \varDelta _y \varDelta _z)^{1/3}$ is the LES filter width, 
$|\tilde {\boldsymbol{\mathsf{S}}}| = (\tilde {\boldsymbol{\mathsf{S}}} : \tilde {\boldsymbol{\mathsf{S}}})^{1/2}$ is the magnitude of the resolved scale strain rate tensor, and 
$c_s$ is the dynamic Smagorinsky coefficient obtained from the LASD model for momentum. The SGS Prandtl number is set to 
$\mbox { {Pr}}^{sgs} = 0.4$ (Kang & Meneveau Reference Kang and Meneveau2002; Kleissl et al. Reference Kleissl, Kumar, Meneveau and Parlange2006).
In the wall model for momentum, Monin–Obukhov similarity is imposed in a local sense (Kumar et al. Reference Kumar, Kleissl, Meneveau and Parlange2006), with filtering at scale 
$2\varDelta$, which has been shown to better reproduce the mean surface stress (Bou-Zeid et al. Reference Bou-Zeid, Meneveau and Parlange2005). For potential temperature, the lower boundary condition corresponds to a spatially variable surface heat flux 
$Q_0(x,y)$ that is imposed at the first vertical grid level, located at height 
$\varDelta z / 2$. The upper boundary condition is stress-free with no flow through the upper boundary (
${\partial (\cdot ) }/{\partial x} = {\partial (\cdot )}/{\partial y} = \tilde {w} = 0 \rvert _{z=L_z}$, where 
$L_z$ is the domain height). In order to prevent the reflection of gravity waves from the upper boundary, a sponge layer is used in the upper 25 % of the domain following the method of Nieuwstadt et al. (Reference Nieuwstadt, Mason, Moeng and Schumann1993). The LES algorithm is parallelized in vertical slabs using message passing interface (MPI).
2.2. Cases
Simulations were conducted on an 
$\{L_x,L_y,L_z\} = \{6\,\mathrm {km}, 6\,\mathrm {km}, 2\,\mathrm {km}\}$ domain at resolution, 
$N_x N_y N_z = 160^{3}$. (Demonstration of resolution insensitivity has been provided in numerous prior articles, for example Salesky & Anderson Reference Salesky and Anderson2018, Reference Salesky and Anderson2019, Reference Salesky and Anderson2020.) Such insensitivity has been established via assessment of relative contribution from resolved and unresolved quantities, and via computation of spectral density of resolved velocity and potential temperature fluctuations. This resolution is less than the highest-resolution environmental LES completed in recent times (Wilczek, Stevens & Meneveau Reference Wilczek, Stevens and Meneveau2015), but is nonetheless adequate for the purposes of this study. The time step was set to 
$\varDelta _t = 0.05$ s. Grid convergence for the homogeneous heat flux convective boundary layer (CBL) simulations on this domain was examined by Salesky et al. (Reference Salesky, Chamecki and Bou-Zeid2017), who found that first- and second-order moments were well-converged on the 
$256^3$ grid (for 
$\varDelta = 25.8$ m, where 
$\varDelta = (\varDelta _x \varDelta _y \varDelta _z)^{1/3}$), and there was not a significant difference between mean vertical profiles on the 
$160^3$ grid (
$\varDelta = 41.3$ m) and the 
$256^3$ grid.
Table 1 provides salient attributes on the range of simulations considered for this article; case H is for a homogeneous heat flux, i.e. 
$\beta = Q_0^H/Q_0^L = 1$, which provides a benchmark against which the relative influences of 
$\beta$ and 
$\lambda = \ell _H/\ell _L$ can be quantified. Alphanumeric case numbering is used for the cases, where ‘A’, ‘B’ and ‘C’ correspond to 
$\lambda = 1/5$, 
$1/2$ and 2, respectively (see also figure 1(b) for graphical depiction), while numeric values are equivalent to 
$\beta$. The spacing between the centres of the high heat flux strips is equivalent, 
$s/z_i = (1/2) L_y/z_i = (\ell _H + \ell _L)/z_i = 2.5$. In complementary work on flow response to aerodynamic roughness heterogeneities (example cases shown in figure 1a), Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015) and Yang & Anderson (Reference Yang and Anderson2017) have shown that this spacing is optimal for sustenance of 
$z_i$-scale secondary flows, which is the focus of this work. The aerodynamic and thermal heterogeneous cases cannot be compared directly, but, as will be shown, the figure 1 cases both induce 
$z_i$-scale circulations.
For all cases, the aggregate heat flux is equivalent to case H – that is, 
$Q_0 = \int _{\mathrm {d}^2 {\boldsymbol x}} Q_0 ({\boldsymbol x}) \,\mathrm {d}^2 {\boldsymbol x} / \int _{\mathrm {d}^2 {\boldsymbol x}} \mathrm {d}^2 {\boldsymbol x}$ – yet we record differing bulk flow properties, as evidenced by the stability parameter, 
$-z_i/L$, and the heterogeneity parameter, 
$\mathcal {H}$ (discussion to follow). For perspective, figure 2 shows bulk properties for the table 1 cases. Figure 2(a) shows the bulk stability parameter and figure 2(b) the heterogeneity parameter, against 
$\beta$ for cases A (squares), B (circles) and C (asterisks). Shown in figure 2(a) is 
$-z_i/L$ for case H; it is emphasized once more that the aggregate heat flux is held constant for all cases, indicating that variable 
$\beta$, and to a far lesser extent 
$\lambda$, causes 
$-z_i/L$ to rise. This is a product of enhanced mixing associated with a priori prescribed spanwise variability in surface fluxes, which induces concurrent variability in Obukhov length that subsequently affects the bulk flow attributes (results to follow).
Figure 2. Graphical summary of table 1 case attributes shown against 
$\beta$. Panels (a,b) show the Reynolds-averaged bulk stability parameter and the heterogeneity parameter, respectively. In (a), 
$-z_i/L$ for the homogeneous case is superimposed, for context, while the direction of increasing 
$\lambda$ is shown in (b). In both panels, square, circle and asterisk datapoints correspond with cases A, B and C, respectively.
Figure 2(b) shows 
$\mathcal {H}$ (see (1.11)) against 
$\beta$ (note that, by definition, 
$\mathcal {H} = 0$ for case H). Here, 
$\mathcal {H}$ is computed via a characteristic thermal gradient
\begin{equation} \delta \theta = \underbrace{\mathrm{arg\ max}}_{y}( \langle \theta \rangle_{xt}(y,z_c) - \theta_0), \end{equation}
where 
$z_c$ is a predefined elevation at which the characteristic thermal gradient is recovered; for the purpose of this work, we will use the first computational level at which data are stored. Here, 
$\theta _0$ is a reference surface temperature, 
$\theta _0 = \langle \theta \rangle _{xyt} (z_c)$. The datapoints illustrate monotonic response to variation in 
$\beta$ and 
$\lambda$. Thus secondary flow intensity – a quantity predicated upon 
$\mathcal {H}$, representative of thermal torque through a modified Richardson number using (1.5) and (1.11) – is clearly evident. There is a temptation to conflate 
$\beta = Q_0^H/Q_0^L$ and 
$\mathcal {H}$, given that the latter is dependent intrinsically on the former since 
$\mathcal {H} \sim \delta \theta$ and 
$\beta = Q_0^H/Q_0^L = (Q_0^L + \delta Q_0)/Q_0^L = 1 + \delta Q_0/Q_0^L$. However, it is also stressed that ambient shear-driven mixing precludes any such direct correlation between a priori prescribed heat fluxes and resultant thermal distributions. It is for this reason that 
$\mathcal {H} \sim \beta ^n$ in figure 2(b), where 
$n \neq 1$; it is also for this reason that 
$\mathcal {H}$ and 
$\beta$ are separated throughout this narrative.
Simulations were forced by a mean pressure gradient force in the streamwise direction 
$\rho _0^{-1} \partial _x P = 5.0 \times 10^{-4}$ m s
$^{-2}$. For reference, this would correspond to a geostrophic wind magnitude 
$\lvert \boldsymbol {U_g} \rvert = 5$ m s
$^{-1}$ at latitude 
$\varphi = 43.3^\circ$ N, where the Coriolis parameter is 
$f = 1.0 \times 10^{-4}$ s
$^{-1}$. However, in order to simplify the analysis of secondary flow processes, we do not consider the Rossby number dependence of secondary flow processes in the present work, and therefore omit the Coriolis force from our present simulations (similar to Salesky & Anderson Reference Salesky and Anderson2019). The aerodynamic roughness length was set to 
$z_0 = 0.10$ m, and the initial CBL depth was set to 
$z_i = 1000\,\text {m}$. The CBL depth was calculated as the height corresponding to the minimum kinematic heat flux 
$\overline {w'\theta '}$. Note that an order of magnitude estimate of the ratio of 
$z_0$ to roughness element height 
$h$ can be taken as 
$z_0/h = 0.1$ (e.g. Raupach, Antonia & Rajagopalan Reference Raupach, Antonia and Rajagopalan1991). Thus 
$z_i/h = 0.1 z_i/z_0 \sim O(10^3)$, which is well within the range where outer similarity can be considered valid (Raupach et al. Reference Raupach, Antonia and Rajagopalan1991; Jimenez Reference Jimenez2004; Volino, Schultz & Flack Reference Volino, Schultz and Flack2007; Wu & Christensen Reference Wu and Christensen2007; Flack & Schultz Reference Flack and Schultz2010). The present results therefore are not sensitive to the value of 
$z_0$ employed. The initial temperature profile was imposed using the three-layer profile described in Sullivan & Patton (Reference Sullivan and Patton2011), i.e.
\begin{equation} \varTheta(z) = \begin{cases} 300\,\text{K}, & z \leq 1000\,\text{m},\\ 300\,\text{K} + (z - 1000\,\text{m}) \, \varGamma_1, & 1000\,\text{m} \leq z < 1100\,\text{m},\\ 308\,\text{K} +(z - 1100\,\text{m}) \varGamma_2, & z \geq 1100\,\text{m}, \end{cases} \end{equation}
where 
$\varGamma _1 = 0.08\,\text {K}\,\text {m}^{-1}$ and 
$\varGamma _2 = 0.003\,\text {K}\,\text {m}^{-1}$. For these simulations, the large-eddy turnover time is 
$T_l \approx z_i/w_\star = 637$ s. Simulations were run for 5 hours of dimensional time, or approximately 
$28.5 T_l$, which required 360 000 steps for LES on the 
$160^3$ grid. Turbulence statistics are derived from the final 
$\approx 11.4 T_l$ of the simulation.
It is emphasized that for the present work, the lower boundary condition for temperature is imposed through surface variable heat fluxes, rather than through spatially-variable surface temperature. In the case of the latter, sensible heat fluxes are a consequence of relative thermal gradients and momentum, as opposed to a prescribed surface heat flux. This modelling choice offers practical convenience, and it has been proven to be a useful approach in past studies of daytime conditions over land as discussed in Basu et al. (Reference Basu, Holtslag, Van De Wiel, Moene and Steeneveld2008).
3. Results
Herein we present first- and second-order turbulence statistics, in addition to constituent terms regulating transport of 
$\langle k \rangle _t$, all of which provide the foundation to demonstrate efficacy of the aforementioned third kind arguments. Where possible, the results are discussed in the context of prior studies on roughness-driven secondary flows, for which ‘second kind’ secondary flows occur due to turbulent torque (Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015). These comparisons are needed to highlight distinct mechanistic differences between the turbulent and thermal torque mechanisms, as reviewed in § 1.
3.1. First- and second-order turbulence statistics
In figure 3, we provide instantaneous flow visualization in the streamwise–spanwise plane (panels a–f) and spanwise–wall-normal plane (panels g–l) at 
$z/z_i = 0.225$ and 
$x/z_i = 2.6$, respectively. The panels include annotations for the width and centre of the elevated heat flux strips, which are defined in the caption. We have shown arbitrarily selected instantaneous realizations for cases A1.25 (figure 3a–c,g–i) and C6 (figure 3d–f ,j–l), which are the limiting cases addressed here with respect to 
$\beta = Q_0^H/Q_0^L$ and 
$\lambda = \ell _H/\ell _L$, as per table 1.
Figure 3. Instantaneous realizations of first-order flow quantities selected at an arbitrary computational time, 
$t_i$, during LES. Panels (a–c,g–i) and (d–f, j–l) correspond to cases A1.25 and C6, respectively. From left to right, columns (a,d,g,j), (b,e,h,k) and (c, f,i,l) correspond to 
$\tilde {u}({\boldsymbol x},t_i)$, 
$\tilde {w}({\boldsymbol x},t_i)$ and 
$\tilde {\theta }({\boldsymbol x},t_i)$, respectively. Panels (a–f) and (g–l) show instantaneous realizations in the streamwise–spanwise plane at elevation 
$z/z_i = 0.225$, and the spanwise–wall-normal plane at streamwise location 
$x/z_i = 2.6$, respectively (
$x/z_i = 2.6$ denoted with solid vertical lines on (a–f); 
$z/z_i = 0.225$ denoted with solid horizontal lines on (g–l)). Dashed lines indicate locations of spanwise heterogeneity; solid horizontal lines on (a–f) and solid vertical lines on (g–l) denote centre of warm strip (see also figure 1).
In the streamwise–spanwise plane, elevated heat flux manifests as an added drag, evident from the relative streamwise momentum deficit aloft the warm strip (figures 3a,d). This relative streamwise momentum deficit corresponds to upwelling (figures 3b,e). The flow response is a product of elevated heat fluxes, which manifest as elevated potential temperature above the strips (figures 3c, f). This flow response appears to be, at least based on instantaneous realizations, consistent with changes in 
$\beta$ and 
$\lambda$: variation in relative strip width and heat flux ratio does not alter the underlying flow physics. In complementary prior work wherein the roughness of adjacent strips was varied, downwelling occurs above the relatively more rough strips (figure 1(a) and Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015). Here, however, elevated heat injection above the warm strips manifests as a buoyancy force that drives upwelling.
In the spanwise–wall-normal plane, we observe that the relative streamwise momentum deficits in figures 3(a,d) occupy the extent of the flow (figures 3g,j). These streamwise momentum deficits exist as ‘plumes’ above the elevated heat flux regions. This flow response is concurrent with the instantaneous upwelling above 
$Q_0^H$, as per figures 3(h,k). Figures 3(i,l) show potential temperature, where we observe the greatest wall-normal thermal gradients above 
$Q_0^H$. Recall (1.3), (1.4) and (1.11), where rudimentary deductions related to the thermal torque illustrate how spanwise thermal gradients can sustain a streamwise secondary flow. Inspection of figures 3(j,l) shows how positive and negative thermal gradients across the warm strips yield a thermal torque responsible for sustenance of the secondary flows, which are evident from inspection of figures 3(g,h,j,k).
The instantaneous flow response is also exhibited in the Reynolds-averaged statistics, as reported in figure 4. Panels (a,c) and (b,d) show Reynolds-averaged streamwise velocity and potential temperature, respectively, with vector maps, 
$\{ \langle \tilde {v} \rangle _{xt}(y,z), \langle \tilde {w} \rangle _{xt}(y,z) \}$, superimposed to highlight secondary flow structure. As per figure 3, we show only results for the limiting parameters captured with cases A1.25 and C6; this structural flow response was consistent across the 
$\lambda$ and 
$\beta$ parameter space considered – although relative secondary flow intensity varied monotonically, as will be reviewed in subsequent results – and for brevity visualizations for intermediate cases are omitted.
Figure 4. Reynolds-averaged first-order flow quantities in the spanwise–wall-normal plane at streamwise location 
$x/z_i = 2.6$, where (a,b) and (c,d) correspond to cases A1.25 and C6, respectively. From left to right, columns (a,c) and (b,d) correspond to 
$\langle \tilde {u} \rangle _t({\boldsymbol x})$ and 
$\langle \tilde {\theta } \rangle _t ({\boldsymbol x})$, respectively, where corresponding vector maps, 
$\{ \langle \tilde {v} \rangle _{xt}(y,z), \langle \tilde {w} \rangle _{xt}(y,z) \}$, are superimposed to demonstrate vortical structure. Dashed lines indicate locations of spanwise heterogeneity; solid vertical lines denote centre of warm strip (see also figure 1). Grey curves denote spanwise profiles of outer-normalized Obukhov length, 
$-L/z_i$, recovered a posteriori.
The Reynolds-averaged contours are structurally consistent with the instantaneous observations: elevated wall-normal thermal gradients above the relatively warm strip impose buoyancy fluxes that induce upwelling over the depth of the flow, introducing the relative streamwise momentum deficit over the depth of the flow and sustaining the secondary flow. Alternatively, this can be interpreted via the thermal torque due to spanwise thermal gradients, as per (1.3) and (1.4). Note, too, that the spatial extent of the mean secondary cells is roughly equivalent between case A1.25 (figures 4a,b) and case C6 (figures 4c,d) – roughly 
$z_i$-scale counter-rotating secondary cells flanking the centreline of the warm strips – yet the width of the high heat flux strips varies significantly between the cases. This result is similar conceptually to prior work on variation of the width of high roughness strips (figure 1(a) and Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015).
In order to consolidate presentation of first-order Reynolds-averaged flow statistics into a single figure, we have elected to show vertical profiles of dispersive fluctuations of streamwise velocity, vertical velocity and temperature in figures 5(a–c), 5(d–f) and 5(g–i), respectively. The profiles are shown above the centre of the high (dashed lines) and low (solid lines) heat flux regions, respectively, as annotated on the figure ordinate labels. In each panel, the effect of varying 
$\beta$ is shown for 
$\lambda = 0.2$ (figures 5a,d,g), 
$\lambda = 0.5$ (figures 5b,e,h) and 
$\lambda = 2$ (figures 5c, f,i), where the outer-normalized Obukhov length is superimposed, for reference.
Figure 5. Vertical profiles of large-scale spatial variability relative to the underlying Reynolds-averaged distribution (dispersive fluctuation) shown above the high (dashed lines) and low (solid lines) heat flux regions, where the former and latter are defined with 
$\delta \tilde {\boldsymbol u}^H(z)/u_\tau = (\langle \tilde {\boldsymbol u} \rangle _{xt} (y^H,z) - \langle \tilde {\boldsymbol u} \rangle _{xyt} (z))/u_\tau$ and 
$\delta \tilde {\boldsymbol u}^L(z)/u_\tau = (\langle \tilde {\boldsymbol u} \rangle _{xt} (y^L,z) - \langle \tilde {\boldsymbol u} \rangle _{xyt} (z))/u_\tau$, respectively, where (a–c) and (d–f) show streamwise and vertical velocity, respectively. Panels (g–i) show dispersive potential temperature, where 
$\delta \tilde {\theta }^H(z)/\theta _\tau = (\langle \tilde {\theta } \rangle _{xt} (y^H,z) - \langle \tilde {\theta } \rangle _{xyt} (z))/u_\tau$ and 
$\delta \tilde {\theta }^L(z)/u_\tau = (\langle \tilde {\theta } \rangle _{xt} (y^L,z) - \langle \tilde {\theta } \rangle _{xyt} (z))/u_\tau$. Here, 
$y^H$ and 
$y^L$ are discrete spanwise locations coincident with the centre of the high and low heat flux regions, respectively (see also figure 1). Panels (a,d,g), (b,e,h) and (c, f,i) correspond to 
$\lambda = 1/5$, 
$1/2$ and 2, respectively. In all panels, the outer-normalized Obukhov length, 
$-L/z_i$, is denoted by a vertical grey band.
Elevated heat fluxes above the high flux region manifest as an ‘added drag’, inducing a relative streamwise momentum deficit (figures 5a–c) and a relative vertical momentum excess (figures 5d–f) associated with plumes of low streamwise momentum fluid throughout the domain. This corresponds with elevated temperatures throughout the depth of the flow (figures 5g–i). This flow construct provides underlying sustenance of the secondary cells, as conservation of mass necessitates concurrent streamwise momentum excess and vertical momentum deficit above the relatively low heat flux region (solid lines). This ‘downwelling’ provides an ongoing supply of relatively cool fluid, as observed in figures 5(g–i). The broad response in all quantities to variable 
$\beta$ is, by itself, preliminary evidence that spanwise heterogeneity in surface heat flux forcing and the resultant thermal torque provide distinctly different mechanistic sustenance for Reynolds-averaged secondary flows (i.e. Prandtl's secondary flow of the third kind).
In all panels of figure 5, it is clear that with increasing 
$\beta$ – i.e. increasingly dramatic spanwise gradient in thermal forcing – the relative dispersive fluctuations increase monotonically. The response to 
$\lambda$ is somewhat less evident, particularly for momentum (figures 5a–f). The exception is temperature, where we observe monotonically declining temperature above the high heat flux region with increasing 
$\lambda$. This result is a consequence of the underlying surface attributes: since cumulative heat flux is equivalent across all cases, even with variable 
$\beta$, and because 
$\lambda = 0.2$ and 
$\lambda = 2$ correspond to relatively narrow and wide high heat flux strips, respectively, it follows that for 
$\lambda = 0.2$ a relatively elevated heat flux is necessitated across the high heat flux area. Likewise, for 
$\lambda = 2$, it follows that a relatively lesser heat flux is needed over the high heat flux region to accomplish the same values of 
$\beta$. The inverse relationship is true, too, for temperature across the low heat flux region.
For 
$\lambda = 0.2$, elevated temperature over the high heat flux regions induces significant upwelling (figure 5d), which declines monotonically with decreasing temperature (figures 5e, f). As would be expected, maximum values in dispersive vertical velocity fluctuation occur above the Obukhov length (vertical grey bar), which highlights the transition between shear- and thermal-dominated production of turbulent kinetic energy. The kinks in the dispersive fluctuations in the temperature profiles near 
$z/z_i = 1$ are due to the increase in temperature in the entrainment zone and are evidence that the inversion strength over the high and low heat flux regions differs from that of the planar average.
To further the analysis, we now shift our attention to higher-order turbulence statistics. Figure 6 shows salient second-order turbulence statistics for cases A1.25 (panels a–c), B2 (panels d–f) and C6 (panels g–i), where the vertical dashed and solid lines denote heat flux transitions and high heat flux centreline, respectively; we show streamwise–wall-normal turbulent stresses, 
$-\langle u^\prime w^\prime \rangle _t$ (panels a,d,g), turbulent transport of potential temperature, 
$-\langle w^\prime \theta ^\prime \rangle _t$ (panels b,e,h), and turbulent kinetic energy (panels c, f,i). In the interests of brevity, we have not presented all second-order statistics. The horizontal lines in figure 6 show 
$-L/z_i$, where 
$L$ is the Obukhov length, which collapses independent parameters governing the ensemble mean wall-normal gradient of streamwise velocity (Obukhov Reference Obukhov1946). The Obukhov length physically represents the wall-normal elevation at which the dominant turbulent kinetic energy production mechanism transitions from mechanical shear to buoyancy:
\begin{equation} L(y) = \frac{-u_\tau(y)^3\,\varTheta_0}{\kappa g\,Q_0(y)},\quad \mathrm{where}\ u_\tau(y) = \left[\frac{\tau^w}{\rho}\,(\hat{\boldsymbol n}(y) : \hat{\boldsymbol n}(y)) \right]^{1/2}. \end{equation}
Here, 
$\hat {\boldsymbol n}(y,z) = \langle \tilde {\boldsymbol u} \rangle _{xt}(y,z)/\langle \tilde {\boldsymbol u} : \tilde {\boldsymbol u})^{1/2} \rangle _{xt}(y,z)$ is the unit vector based on the velocity vector, 
$\kappa$ is the von Kármán constant, and 
$Q_0(y) = Q_0^L$ or 
$Q_0^H$, depending on location across the surface (see also figure 1 for illustration of heat flux variation). For unstable flows such as those considered in this article (table 1), 
$-L/z_i$ provides a literal visualization of the mechanisms responsible for 
$k$ production. In this sense, superposition of 
$(-z_i/L)^{-1}$ is instructive for the present cases.
Figure 6. Contours of turbulent stresses for cases A1.25 (a–c), B2 (d–f) and C6 (g–i) in the spanwise–wall-normal plane at streamwise location 
$x/z_i = 2.6$. From left to right, columns show 
$-\langle u^\prime w^\prime \rangle _t$ (a,d,g), 
$\langle w^\prime \theta ^\prime \rangle _t$ (b,e,h), and 
$k = (1/2)\langle {\boldsymbol u}^\prime : {\boldsymbol u}^\prime \rangle _t$ (c, f,i). Dashed lines indicate locations of spanwise heterogeneity, solid vertical lines denote centre of warm strip (see also figure 1). Grey curves denote spanwise profiles of outer-normalized Obukhov length, 
$-L/z_i$, recovered a posteriori.
Recall that a homogeneous roughness is prescribed for the present cases, yet figures 6(a,d,g) demonstrate that 
$-\langle u^\prime w^\prime \rangle _t$ exhibits large-scale spanwise heterogeneity across the surface. This highlights how buoyancy-driven secondary flows redistribute 
$k$ throughout the domain, such that the high and low heat flux regions manifest physically as relatively more and less rough, respectively. Here, 
$(-z_i/L)^{-1}$ is lowest over the 
$Q_0^H$ region, indicating that production of 
$k$ by buoyancy is larger over this region, relative to mechanical shear, even given the larger 
$-\langle u^\prime w^\prime \rangle _t$ across the high heat flux strip. In contrast, across the low heat flux region, 
$(-z_i/L)^{-1}$ rises dramatically, illustrating the relative dominance of mechanical shear production across the 
$Q_0^L$ region.
Figures 6(b,e,h) show the corresponding vertical turbulent flux of potential temperature, 
$\langle w^\prime \theta ^\prime \rangle _t$. The spatial distributions of 
$\langle w^\prime \theta ^\prime \rangle _t$ resemble 
$-\langle u^\prime w^\prime \rangle _t$; note, however, that the colourbar limits for (a,d,g) and (b,e,h) are 
$-0.1 \lesssim z/z_i \lesssim 1.5$ and 
$-2 \lesssim z/z_i \lesssim 6$, respectively, which confounds visual inference on the value of 
$(-z_i/L)^{-1}$. Thus although the 
$-\langle u^\prime w^\prime \rangle _t$ and 
$\langle w^\prime \theta ^\prime \rangle _t$ contours are similar, in fact the relative difference in magnitude illustrates why 
$-z_i/L$ declines across the 
$Q_0^H$ strip.
Finally, figures 6(c, f,i) provide spatial contours of turbulent kinetic energy for the cases noted in the caption. The panels show plumes of elevated 
$k$ above the 
$Q_0^H$ strips. The plume structure is consistent with those for 
$-\langle u^\prime w^\prime \rangle _t$ and 
$\langle w^\prime \theta ^\prime \rangle _t$, which is consistent with the underlying production mechanisms as per (1.6).
The figures demonstrate that second-order statistics increase monotonically with 
$\beta$, evidenced by relative increases in the magnitudes of 
$-\langle u^\prime w^\prime \rangle _t$, 
$\langle w^\prime \theta ^\prime \rangle _t$ and 
$k$ over the high heat flux region. Note, however, that the spatial distributions in these quantities do not respond materially to 
$\lambda$. This is consistent with the Reynolds-averaged first-order statistics shown in figures 4 and 5, where spanwise heat flux variability yields counter-rotating 
$z_i$-scale secondary cells centred about the centreline of the 
$Q_0^H$ region.
3.2. Interpreting secondary flows from turbulent kinetic energy transport
In this subsection, the inertia-dominated (Reynolds-averaged) 
$k$ transport equation (1.6) is used to discern the underlying physics responsible for the secondary flow processes presented in § 3.1. This follows from Hinze (Reference Hinze1967), who used 
$k$ transport to characterize secondary flow processes in turbulent duct flow, arguing that secondary flow advection is necessitated by an imbalance in production and dissipation, such that ‘turbulence poor’ fluid is entrained by regions of ‘turbulence rich’ fluid. This argument has been used to infer secondary transport in roughness-driven secondary flows in hydraulic flows (Wang & Cheng Reference Wang and Cheng2005; Vermaas et al. Reference Vermaas, Uijttewall and Hoitink2011) and engineering roughness (Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015). Others have used the Reynolds-averaged transport equation for in-plane vorticity (Hoagland Reference Hoagland1960), but balancing this equation is problematic. In particular, large residuals usually are found when calculating terms in the mean streamwise vorticity budget, regardless of whether one calculates it from experimental or numerical simulation data. This may be due to the spontaneous reversals of the rotational direction of secondary cells that complicate averaging (Anderson Reference Anderson2019), although the underlying reason is not abundantly clear. We therefore declined this avenue in favour of the 
$\langle k \rangle _t$ balance equation (1.6) (Madabhushi & Vanka Reference Madabhushi and Vanka1991; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015).
As summarized in (1.6), transport of Reynolds-averaged 
$\langle k \rangle _t$ is regulated by the interplay between production and dissipation, where the former occurs via mechanical shear (
$\mathcal {P}_{S}$) and buoyancy (
$\mathcal {P}_{B}$), while the latter is quantified here for the inertia-dominated (Reynolds-averaged) form of (1.6), such that 
$\epsilon = \langle - \boldsymbol{\mathsf{T}} : \boldsymbol{\mathsf{S}} \rangle _{xt}$. Figure 7 shows contours for constituent terms within (1.6) for cases A1.25 (panels a,b,g,h), B2 (panels c,d,i,j) and C6 (panels e, f,k,l), where the colour flood contours and isocontour values are defined in the figure caption.
Figure 7. Constituent terms of the 
$\langle k \rangle _t$ transport equation for cases A1.25 (a,b,g,h), B2 (c,d,i,j) and C6 (e, f,k,l). Panels show shear production, 
$\mathcal {P}_{S}$ (a,c,e), buoyancy production, 
$\mathcal {P}_{B}$ (b,d, f), advection, 
$\langle \tilde {\boldsymbol u} \rangle _t \boldsymbol {\cdot } \boldsymbol {\nabla } \langle k\rangle _t$ (g,i,k), and dissipation (as it appears in (1.6)), 
$- \epsilon = - \langle - \boldsymbol{\mathsf{T}} : \boldsymbol{\mathsf{S}} \rangle _{xt}$ (h,j,l); see also (1.6). Included in the panels are isocontours of vertical velocity with values 
$\tilde {w}({\boldsymbol x},t)/u_\tau = 0.2$ (red) and 
$-0.2$ (blue). Dashed lines indicate locations of spanwise heterogeneity; solid vertical lines denote centre of warm strip (see also figure 1). Horizontal line denotes spanwise profile of outer-normalized Obukhov length, 
$-L/z_i$, recovered a posteriori (see also (3.1)).
In the results shown – which, we emphasize, are also representative of intermediate table 1 cases – we see a consistent pattern of upwelling and downwelling over the high and low heat flux surface regions, respectively. Comparison of shear-driven 
$\langle k \rangle _t$ production 
$\mathcal {P}_{S}$ (figures 7a,e,i) shows similar patterns, although 
$\lambda$ changes substantially across these cases. We see a consistent ‘plume’ of elevated 
$\mathcal {P}_{S}$ across the centre of the 
$Q_0^H$ region, and this pattern is ostensibly insensitive to 
$\lambda$. This result contrasts against the superimposed profiles of outer-normalized Obukhov length 
$-L/z_i$, which show that shear driven-production is relatively weakest across the warm strips. Noting that 
$(-z_i/L)^{-1} \rightarrow 0$ and 
$(-z_i/L)^{-1} \rightarrow \infty$ for free convection and canonical channel turbulence, respectively, large-scale spatial spanwise variability in 
$(-z_i/L)^{-1}$ highlights the distinct ways that spanwise heat flux heterogeneity affects the resultant flow processes. Recall, however, that turbulent transport of temperature is relatively greater across the warm strip, as per figures 6(b,e,h). With this, the dominant production mechanism above the 
$Q_0^H$ strip transitions from mechanical shear (figures 7a,c,e) to buoyancy (figures 7b,d, f).
Advection of Reynolds-averaged turbulent kinetic energy, 
$\mathcal {A}_k = - \langle \tilde {\boldsymbol u} \rangle _t \boldsymbol {\cdot }\boldsymbol {\nabla } \langle k \rangle _t$, is shown in figures 7(g,i,k) (see also (1.6) and accompanying text). The advection contours provide a unique avenue to interpret the resultant secondary flow: across the centre of the 
$Q_0^H$ strip, 
$- \langle \tilde {\boldsymbol u} \rangle _t \boldsymbol {\cdot } \boldsymbol {\nabla } \langle k \rangle _t \rightarrow - \langle \tilde {w} \rangle _t\partial _z \langle k \rangle _t$. Here, vertical velocity is positive, and the vertical gradient of turbulent kinetic energy is negative; thus advection is necessarily positive given the relative imbalance between the production and dissipation of turbulence (Hinze Reference Hinze1967; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015). This deduction illustrates that since 
$\partial _z \langle k \rangle _t > 0$ (as per figures 6c, f,i), 
$\langle \tilde {w} \rangle _t$ must also be positive, which is reported above the centre of the warm strip (figures 3b,e,h,k, 4 and 5d,e, f).
The negative of Reynolds-averaged dissipation, which serves as the principal sink term in the 
$\langle k \rangle _t$ transport equation, 
$-\epsilon = -\langle - \boldsymbol{\mathsf{T}} : \boldsymbol{\mathsf{S}} \rangle _{xt}$, is shown in figures 7(h,j,l). For the non-Hamiltonian conservative systems considered herein, 
$\epsilon$ is everywhere positive; thus 
$-\epsilon$ is negative definite and serves as a sink of 
$\langle k \rangle _t$. The spatial distribution of 
$-\epsilon$ reaches its maximum magnitude in the roughness sublayer and within the base of the shear-dominated 
$\langle k \rangle _t$ production region above 
$Q_0^H$.
When characterizing mechanisms responsible for turbulent secondary flow, Hinze (Reference Hinze1967) introduced the notion of ‘turbulence poor’ and ‘turbulence rich’ fluid, where the former is necessarily entrained by the latter, setting in motion the secondary flow direction. Later, Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015) used this argument to discern underlying mechanisms responsible for roughness-driven secondary flows in channel flows over surfaces depicted by figure 1(a). For such flows, figure 8(a) illustrates how a pronounced production–dissipation imbalance across the high-roughness strip, 
$z_0^H$, necessitates downwelling of low-
$\langle k \rangle _t$ fluid from directly above (dark red region in figure 8a). This resultant mean flow yields a large-scale undulation in 
$\langle \tilde {u} \rangle _t$, such that a mean excess and deficit in 
$\langle \tilde {u} \rangle _t$ occurs above 
$z_0^H$ and 
$z_0^L$, respectively. This has previously been defined as low- and high-momentum pathways (LMP, HMP) (Barros & Christensen Reference Barros and Christensen2014).
Figure 8. Schematic illustrating problem parameters addressed in this study (see also table 1 for quantitative details). Panel (a) shows the case of roughness-driven secondary flows, where the turbulent torque term and production–dissipation imbalance above the relatively high-roughness region, 
$z_0^H$, induce downwelling and formation of high- and low-momentum pathways as shown (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Barros & Christensen Reference Barros and Christensen2014; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015); upwelling occurs above the relatively low-roughness region, 
$z_0^L$. In contrast, for the cases considered herein, shear-driven 
$k$ production is also concentrated above the relatively elevated flux region, but this region is itself capped by a region of buoyancy production (b). This additional 
$k$ production necessitates lateral entrainment of low-
$k$ fluid, reversing the secondary flow direction such that an LMP forms over 
$Q_0^H$. The outer-normalized Obukhov lengths for the spanwise heterogeneous (grey) and homogeneous (black) heat flux cases are illustrated in (b); shown also in (b) are 
$\ell _H$ and representative 
$\ell _d$.
In contrast, the case of spanwise heterogeneous heat fluxes – as per figure 8(b) – yields a complete reversal in secondary flow direction as compared to the roughness heterogeneity case. This is quantified by spanwise variability of the Obukhov length, where 
$-L/z_i$ exhibits a relative minimum above 
$Q_0^H$ due to elevated shear (dark red region in figure 8b). Above 
$-L/z_i$, however, buoyancy continues active production of turbulent kinetic energy (yellow region in figure 8b). As such, a relative production–dissipation imbalance persists over the depth of the flow, and ‘turbulence poor’ fluid can not readily be entrained from aloft. In fact, the production–dissipation imbalance can be resolved only via lateral entrainment, thereby reversing the mean secondary flow polarity relative to the roughness-driven case and resulting in formation of an LMP above 
$Q_0^H$. We have provided annotations on figure 8 for the locations of HMPs and LMPs (Barros & Christensen Reference Barros and Christensen2014), which summarizes polarity reversal due to thermal torque.
3.3. Reynolds-averaged flow attributes
In this final subsection, we focus on aggregating bulk statistics from the suite of cases used in this study and summarized in table 1. Combining bulk results in this manner provides an avenue to identify underlying trends in turbulence statistics with variation in 
$\beta$ and 
$\lambda$. We elected to provide profiles showing spanwise distribution of the similarity parameter, 
$-L/z_i$, and compensated circulation (discussion to follow). These results complement quantities already presented in preceding sections.
Figures 9(a), 9(b) and 9(c) provide spanwise profiles of 
$-L/z_i$ for cases in groups A, B and C, respectively, where specific cases are summarized in the panel legend. In all panels, we have superimposed the reference value, 
$(-z_i/L)^{-1} = 0.26$, for case H, which is the homogeneous heat flux case. For all table 1 cases, the area-averaged heat flux, 
$Q_0 = \int _{\mathrm {d}^2 {\boldsymbol x}} Q_0 ({\boldsymbol x})\,\mathrm {d}^2 {\boldsymbol x} / \int _{\mathrm {d}^2 {\boldsymbol x}}\,\mathrm {d}^2 {\boldsymbol x}$, is equivalent (including case H). Figure 9 illustrates monotonic trends in the deviation in 
$(-z_i/L)^{-1}$ relative to the values for a homogeneous heat flux: as 
$\beta$ increases, shear-driven production of turbulence is relatively weaker and stronger above 
$Q_0^H$ and 
$Q_0^L$, respectively. As per (3.1), the Obukhov length is based on the prescribed heat flux and corresponding surface friction velocity; the ‘step-wise’ profile in 
$(-z_i/L)^{-1}$ is thus a product of the prescribed heat fluxes as per figure 1(b).
Figure 9. Spanwise profiles of Obukhov length for all cases, recovered a posteriori from LES (profiles defined in panel legends). Dashed lines indicate locations of spanwise heterogeneity; solid vertical lines denote centre of warm strip (see also figure 1). Horizontal dashed line denotes 
$-L/z_i$ for case H.
In addition to figure 9, figure 10 shows datapoints for compensated circulation:
\begin{equation} {\varGamma}^* = \frac{{\varGamma}}{w_\star z_i} = \int_{\mathrm{d}^2 {\boldsymbol x}} |\langle \tilde{\omega}_x \rangle_{xt} (y,z)| \,\mathrm{d}^2 {\boldsymbol x}, \end{equation}
where the absolute value of Reynolds-averaged streamwise vorticity, 
$|\langle \tilde {\omega }_x \rangle _{xt} (y,z)|$, is used to quantify secondary flow intensity (Yang & Anderson Reference Yang and Anderson2017). Without computation based on the absolute value of streamwise vorticity, 
${\varGamma }^* \rightarrow 0$ due to the presence of Reynolds-averaged counter-rotating secondary cells. Figure 10 shows 
${\varGamma }^*$ against 
$\mathcal {H}$ (1.11), where the symbols are defined in the figure. Figure 10 is shown with ordinate range 
$0 \leq \varGamma ^* \leq 15$, where the lower limit, 
$\varGamma ^* = 0$, is an important benchmark for assessing the influence of spanwise heterogeneity: for spanwise spacing, 
$s/z_i = 0$ and 
$s/z_i = \infty$, there is no differential heat flux to induce thermal torque, thus 
$\varGamma ^* \rightarrow 0$ is expected in the Reynolds-averaged limit. Between these limits, however, the preceding results have illustrated how thermal torque sustains 
$z_i$-scale secondary cells (figures 3, 4, 5, 6 and 7).
Figure 10. Datapoints for compensated circulation (see (3.2)) for the cases considered in this study (with the exception of case H, for which 
${\varGamma }^* = 0$; see also table 1). Values of 
${\varGamma }^*$ are shown against 
$\mathcal {H} ({u_\tau }/{w_\star }) = \varOmega _{B} / \varOmega _{A}$ (abscissa; see table 1 for case attributes) for the values of 
$\lambda$ denoted in the figure.
The figure 10 results can be understood readily via consideration of the scaling arguments presented in § 1.1. Since 
$\langle \tilde {\omega }_x \rangle _{xt} \sim w_\star /\ell _H$ (see (1.9)) and 
$\mathrm {d}^2 {\boldsymbol x} \sim \ell _H \ell _d \sim \ell _H^2 (w_\star /u_\tau )$ (see (1.8)), it follows that 
$\varGamma ^* \sim \ell _d/z_i$. (Recall that 
$\ell _H$ and 
$\ell _d$ correspond to the heterogeneity length scale and the depth of the secondary cells, respectively.) Therefore, given that 
$z_i$ presents only small variations between the cases analysed in this work (see table 1), potential changes in compensated circulation should then be mostly attributed to changes in 
$\ell _d$. However, from the results presented earlier, in figures 3 and 7, it is found that 
$\ell _d$ does not significantly respond to changes in 
$\beta$ and 
$\lambda$. This explains the almost-constant behaviour of compensated circulation as a function of 
$\varOmega _{B} / \varOmega _{A}$, as shown in figure 10. Furthermore, this result is clear evidence that in all cases considered in this work, the thermal torque provides sustenance for the balance of Reynolds-averaged streamwise vorticity, and hence they are examples of Prandtl's secondary flow of the third kind.
4. Conclusions
In recent years, a significant research effort has been dedicated towards assessing turbulent wall flow response to surfaces with spanwise aerodynamic/hydrodynamic roughness heterogeneity (Wang & Cheng Reference Wang and Cheng2005; Vermaas et al. Reference Vermaas, Uijttewall and Hoitink2011; Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Kevin et al. Reference Kevin, Monty, Bai, Pathikonda, Nugroho, Barros, Christensen and Hutchins2017; Yang & Anderson Reference Yang and Anderson2017; Chung et al. Reference Chung, Monty and Hutchins2018; Hwang & Lee Reference Hwang and Lee2018; Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2018; de Silva et al. Reference de Silva, Kevin, Baidya, Hutchins and Marusic2018; Anderson Reference Anderson2019; Kevin et al. Reference Kevin, Monty and Hutchins2019; Vanderwel et al. Reference Vanderwel, Stroh, Kriegseis, Frohnapfel and Ganapathisubramani2019; Stroh et al. Reference Stroh, Schäfer, Frohnapfel and Forooghi2020; Wangsawijaya et al. Reference Wangsawijaya, Baidya, Chung, Marusic and Hutchins2020; Zampiron et al. Reference Zampiron, Cameron and Nikora2020, Reference Zampiron, Cameron and Nikora2021). This body of work has focused on flow responses to a variety of spanwise-heterogeneous roughness arrangements, but in all cases the resultant outcome has been equivalent conceptually: such surfaces induce and sustain Reynolds-averaged secondary flows (that is, secondary circulations aligned with the primary transport direction). For canonical and complex arrangements, it has been shown that these secondary flows are a manifestation of Prandtl's secondary flow of the second kind (Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015), as per Bradshaw (Reference Bradshaw1987) and references therein.
In the present work, we have used scaling arguments and a comprehensive suite of LES cases to quantify flow response to surfaces with equivalent aerodynamic roughness, but spanwise-variable prescribed heat flux (i.e. a ‘thermal analogue’ of the spanwise-heterogeneous rough walls, as per figure 1). A priori assessment of the equations governing transport of Reynolds-averaged streamwise vorticity and turbulent kinetic energy highlight the potential for spanwise-variable heat fluxes to induce a thermal torque, and the aforementioned scaling arguments indicate that this thermal torque is capable of a non-trivial contribution to the sustenance of Reynolds-averaged streamwise vorticity. This result, by itself, is significant: although spatial heterogeneity in heat fluxes induces subsequent spatial heterogeneity in the turbulent stresses – sustaining Prandtl's secondary flows of the second kind – scaling arguments indicate that thermal torque provides an additional contribution to the balance of streamwise vorticity. It is for this reason that the resultant mechanism compels generalized nomenclature, where herein we demonstrate that thermal torque induces Prandtl's secondary flows of the third kind.
We assembled a suite of LES cases designed to capture flow response over a range of salient parameters. The results indeed confirm that spanwise heterogeneous heat flux and thermal torque alter the balance of streamwise vorticity. We showed a series of first- and second-order turbulence statistics to illustrate this result, and in all cases we recover monotonic trends with variation in forcing parameters; these quantities were also used to capture transport of 
$\langle k \rangle _t$. We also illustrate how geometric attributes of the surface serve to regulate the spatial extent of secondary cells and their intensity.
Funding
This work was supported by the US Department of Energy Atmospheric System Research Program, grant no. DE-SC0022124 (STS), and the National Science Foundation, grant no. AGS-1500224 (WA). M.C. acknowledges the support of National Science Foundation grant nos PDM-1649067 and PDM-1712538, as well as the support of the Alexander von Humboldt Stiftung/Foundation, Humboldt Research Fellowship for Experienced Researchers.
Declaration of interests
The authors report no conflict of interest.
Open access
