Quantifying transport by strongly stratified turbulence in low Prandtl number ($Pr$) fluids is critically important for the development of better models for the structure and evolution of stellar and planetary interiors. Motivated by recent numerical simulations showing strongly anisotropic flows suggestive of a scale-separated dynamics, we perform a multiscale asymptotic analysis of the governing equations. We find that, in all cases, the resulting slow–fast systems take a quasilinear form. Our analysis also reveals the existence of several distinct dynamical regimes depending on the emergent buoyancy Reynolds and Péclet numbers, $Re_b = \alpha ^2 Re$ and $Pe_b = Pr Re_b$, respectively, where $\alpha$ is the aspect ratio of the large-scale turbulent flow structures, and $Re$ is the outer-scale Reynolds number. Scaling relationships relating the aspect ratio, the characteristic vertical velocity and the strength of the stratification (measured by the Froude number $Fr$) naturally emerge from the analysis. When $Pe_b \ll \alpha$, the dynamics at all scales is dominated by buoyancy diffusion, and our results recover the scaling laws empirically obtained from direct numerical simulations by Cope et al. (J. Fluid Mech., vol. 903, 2020, A1). For $Pe_b \ge O(1)$, diffusion is negligible (or at least subdominant) at all scales and our results are consistent with those of Chini et al. (J. Fluid Mech., vol. 933, 2022) for strongly stratified geophysical turbulence at $Pr =O(1)$. Finally, we have identified a new regime for $\alpha \ll Pe_b \ll 1$, in which slow, large scales are diffusive while fast, small scales are not. We conclude by presenting a map of parameter space that clearly indicates the transitions between isotropic turbulence, non-diffusive stratified turbulence, diffusive stratified turbulence and viscously dominated flows, and by proposing parameterisations of the buoyancy flux, mixing efficiency and turbulent diffusion coefficient for each regime.

In the generalized quasilinear approximation (GQLA) (Marston et al., Phys. Rev. Lett., vol. 116, 2016, 214501), a threshold wavenumber ($k_0$) in the direction of translational symmetry segregates the total into large- ($l$) and small-scale ($h$) fields. While the governing equation for the large-scale field is fully nonlinear, that for the small scales is linearized with respect to the large-scale field. In addition, some nonlinear triad interactions are omitted in the GQLA. Herein, the GQLA is applied to two-dimensional planar Rayleigh–Bénard convection (RBC). A scale separation between the large-scale convection rolls and small-scale turbulent fluctuations is typical in RBC. The present work explores the efficacy of GQLA in capturing the scale-by-scale energy transfer processes in RBC. The initial condition for the GQLA simulations was either the statistically stationary state obtained in direct numerical simulation (DNS) or random fluctuations superimposed on the linear conductive temperature profile and $\boldsymbol {u}=0$. The GQLA simulations can capture the convection rolls for $k_0$ larger than or equal to the dominant wavenumber for thermal driving of the flow ($k_0 \ge k_{\hat {Q}}$). Additionally, the GQLA emulates the fully nonlinear dynamics for $k_0$ larger than or equal to the first harmonic of the convection-roll wavenumber ($k_0 \ge 2k_{roll}$). In the intermediate regime with $2k_{roll} > k_0 \ge k_{\hat {Q}}$, the dynamics captured in GQLA simulations is different from the DNS. In DNS, two primary energy transfer processes dominate: (i) the energy transfer to/from the convection rolls and (ii) the scale-by-scale inverse kinetic energy and forward thermal energy cascades mediated by the convection rolls. The fully nonlinear dynamics is emulated by GQLA when these energy transfer processes are faithfully reproduced. Utilizing the framework of altering the triad interactions in GQLA, an additional intrusive calculation, including target triad interactions, is performed here to study their influence. This intrusive calculation shows that the convection rolls are not captured in GQLA for $k_0 < k_{\hat {Q}}$ because of the exclusion of the $h \to l \to l$ and $l \to h \to h$ triad interactions in GQLA. The inclusion of these triad interactions in the intrusive calculation yields the convection rolls, and the reproduced dynamics is similar to that of the intermediate GQLA regime with $2k_{roll} > k_0 \ge k_{\hat {Q}}$.

We present a new derivation of the kinetic equation for weak, non-hydrostatic internal gravity wave turbulence. The equation is equivalent to the one obtained by Caillol & Zeitlin (Dyn. Atmos. Oceans, vol. 32, issue 2, 2000, pp. 81–112), but it takes a canonical form. We show that it conserves the energy without involving the resonance condition in frequency, and look for the isotropic part of the steady, scale-invariant solutions. We provide a parametrization of the resonant manifold of non-hydrostatic internal gravity wave triadic interactions. This allows us to simplify the collision integral, and to evaluate the transfer coefficients of all triadic interactions. In the hydrostatic limit, our equation is equivalent to the Hamiltonian description of Lvov & Tabak (Phys. Rev. Lett., vol. 87, issue 16, 2001, 168501).

In this study, the resolvent-based estimation (RBE) is further generalised to cases with arbitrarily sampled measurements in time, where the generalised RBE is denoted as GRBE in this study. Different from the RBE that constructs the transfer function at each frequency, the GRBE minimises the estimation error energy in the physical temporal domain by considering the forcing and noise statistics. The GRBE is validated by estimating the complex Ginzburg–Landau equation and turbulent channel flows with the friction Reynolds number $Re_{\tau }=186$, 547 and 934, where the results from the RBE are also included. When the measurements are temporally resolved, the estimation results of the two approaches are equivalent to each other, and both match well with the reference numerical results. For the temporally unresolved cases, the estimation errors from the GRBE are obviously lower than those from the RBE. After validation, the GRBE is applied to investigate the impacts of the abundance of the measured information, including the temporal information and sensor types, on the estimation accuracy. Compared with the mean square error (MSE) in the estimation with temporally resolved measurements, that with measurements at only one snapshot, i.e. without any temporal information, increases by approximately $15\,\%$. On the other hand, it can effectively improve the estimation accuracy by increasing the number of sensor types. With temporally resolved measurements, the relative MSE decreases by $12.3\,\%$ when the sensor types increase from $\lbrace \tau _u \rbrace$ to $\lbrace \tau _u,\tau _w,p \rbrace$, where $\tau_u$, $\tau_w$ and p are the streamwise shear stress, spanwise shear stress and pressure at the wall. Finally, several existing forcing models are incorporated into the GRBE to investigate their performance in the linear estimation of flow state. The wall-distance-dependent model (W-model) results match well with the optimal linear estimations when the measurements are temporally unresolved. Meanwhile, with the increase of temporal information of the measurement, the estimation errors from the tested W-model and the scale-dependent model ($\lambda$-model) both increase, which contradicts the tendency observed in the optimal linear GRBE estimation results. Such a phenomenon highlights the importance of proper modelling of the forcing in the temporal domain for the accuracy of flow state estimation.

Time-resolved two-dimensional two-component particle image velocimetry measurements with high spatial resolution are carried out in a water tank agitated by four blades rotating at constant speed. Different blade geometries and rotation speeds are tested for the purpose of modifying turbulent flow conditions. In all cases where no baffles are used to break the rotation, the Zeman length is an order of magnitude smaller than the Taylor length. Compared with the cases with baffles which break the rotation, in the unbaffled cases turbulence production and/or mean advection are significant and the turbulence nonlinearity is dramatically reduced for the horizontal (i.e. normal to the axis of rotation) two-point turbulence fluctuating velocities. This nonlinearity reduction is manifest not only in the interscale turbulent energy transfer but also in the interspace turbulent energy transfer, which nearly vanishes. However, the nonlinearity is not reduced for the vertical two-point turbulence fluctuating velocities: the corresponding interscale turbulent transfer rate is in fact intensified, and its dependence on the two-point separation distance, as well as that of the corresponding interspace turbulent transfer rate, which does not vanish, is significantly modified. Even though non-homogeneities are very different for different blades and rotation speeds in the unbaffled cases, the horizontal fluctuating velocity's second-order structure function collapses with scalings which resemble predictions for homogeneous turbulence subject to strong rotation. The vertical fluctuating velocity's second-order structure function does not collapse for different blade geometries by neither these nor the Kolmogorov predictions.

Large-scale coherent structures in incompressible turbulent pipe flow are studied for a wide range of Reynolds numbers ($Re_\tau =180, 550, 1000, 2000$ and $5200$). Employing the Karhunen–Loève decomposition and a novel approach based on the Voronoi diagram, we identify and classify statistically coherent structures based on their location, dimensions and $Re_{\tau }$. With increasing $Re_{\tau }$, two distinct classes of structures become more energetic, namely wall-attached and detached eddies. The Voronoi methodology is shown to delineate these two classes without the need for specific criteria or thresholds. At the highest $Re_{\tau }$, the attached eddies scale linearly with the wall-normal distance with a slope of approximately $l_y\sim 1.2y/R$, while the detached eddies remain constant at the size of $l_y \approx 0.26R$, with a progressive shift towards the pipe centre. We extract these two classes of structures and describe their spatial characteristics, including radial size, helix angle and azimuthal self-similarity. The spatial distribution could help explain the differences in mean velocity between pipe and channel flows, as well as in modelling large and very-large-scale motions (LSM and VLSM). In addition, a comprehensive description is provided for both wall-attached and detached structures in terms of LSM and VLSM.

The dynamics of turbulent kinetic energy (TKE), turbulence dissipation rate (TDR) and turbulence production rate (TPR) are explored in fully developed turbulent channel flow using direct numerical simulations up to $\textit {Re}_\tau \approx 2000$ with minimal computational box for large-scale structures. Time correlation analysis based on volume-averaged TKE and TDR shows a well-defined average time lag, as in periodic/homogeneous turbulence, which, unlike periodic/homogeneous turbulence, appears to be Reynolds-number-dependent. On the basis of a spatio-temporal correlation analysis, we show that plane-averaged TKE fluctuations in the near-equilibrium region are transported towards both the core and near-wall regions, and are positively correlated with plane-averaged TDR fluctuations there with combined wall-distance and time lags. In the path towards the core region, the wall-distance lag is very close to the time lag multiplied by the friction velocity. The path towards the near-wall region has a wide spread of time lags, which increases with Reynolds number. The spatio-temporal correlation paths both towards the core and towards the wall are reproduced when the reference plane TKE is conditionally averaged on either ejections or sweeps, and are in fact stronger in correlation values in the case of ejections, which are better organised than sweeps. While volume-averaged TPR evidently precedes volume-averaged TKE, a more complex picture of non-local space–time correlations between reference plane TKE and TPR is revealed. A mechanistic model is proposed to elucidate these correlations between TKE and TPR through the interaction between the mean shear and the Reynolds shear stress.

Two-dimensional turbulence transfers its energy towards the lowest mode in the domain, but domain geometry exerts a powerful control. On the sphere, with its three axes of rotational symmetry, angular momentum conservation prevents energy from entering the three lowest modes – those corresponding to the spherical harmonics $Y_1^0$ and $Y_1^{\pm 1}$ – because the amplitudes of these three modes are proportional to the three conserved components of the angular momentum vector. Non-spherical ellipsoids partly or completely break the rotational symmetry corresponding to angular momentum conservation. The flow on spheroids, which have only one axis of rotational symmetry, conserves only a single component of angular momentum. If the axis of symmetry is taken to be the $z$-axis, then only the $z$-component of angular momentum is conserved. Energy can flow into the other two lowest modes. The general triaxial ellipsoid breaks all rotational symmetries, thus angular momentum is not conserved, and energy can flow into any mode. We describe numerical experiments that confirm these predictions.

In the absence of large-scale coherent structures, a widely used statistical theory of two-dimensional turbulence developed by Kraichnan, Leith, and Batchelor (KLB) predicts a power-law scaling for the energy, $E(k)\propto k^\alpha$ with an integral exponent $\alpha ={-3}$, in the inertial range associated with the direct cascade. A power-law scaling is also observed in the presence of coherent structures, but the scaling exponent becomes fractal and often differs substantially from the value predicted by the KLB theory. Here we present a dynamical theory that sheds new light on the relationship between the spatial and temporal structure of the large-scale flow and the scaling of small-scale structures representing filamentary vorticity. Specifically, we find hyperbolic regions of the large-scale flow to play a key role in the flux of enstrophy between scales. Small-scale vorticity in these regions can be described by dynamically self-similar solutions of the Euler equation, which explains the power-law scaling. Furthermore, we find that correlations between different hyperbolic regions are responsible for the emergence of fractal scaling exponents.

Recently, direct numerical simulations (DNS) of stably stratified turbulence have shown that as the Prandtl number ($Pr$) is increased from 1 to 7, the mean turbulent potential energy dissipation rate (TPE-DR) drops dramatically, while the mean turbulent kinetic energy dissipation rate (TKE-DR) increases significantly. Through an analysis of the equations governing the fluctuating velocity and density gradients we provide a mechanistic explanation for this surprising behaviour and test the predictions using DNS. We show that the mean density gradient gives rise to a mechanism that opposes the production of fluctuating density gradients, and this is connected to the emergence of ramp cliffs. The same term appears in the velocity gradient equation but with the opposite sign, and is the contribution from buoyancy. This term is ultimately the reason why the TPE-DR reduces while the TKE-DR increases with increasing $Pr$. Our analysis also predicts that the effects of buoyancy on the smallest scales of the flow become stronger as $Pr$ is increased, and this is confirmed by our DNS data. A consequence of this is that the standard buoyancy Reynolds number does not correctly estimate the impact of buoyancy at the smallest scales when $Pr$ deviates from 1, and we derive a suitable alternative parameter. Finally, an analysis of the filtered gradient equations reveals that the mean density gradient term changes sign at sufficiently large scales, such that buoyancy acts as a source for velocity gradients at small scales, but as a sink at large scales.

We systematically study the dissipative anomaly in compressible magnetohydrodynamic (MHD) turbulence using direct numerical simulations, and show that the total dissipation remains finite as viscosity diminishes. The dimensionless dissipation rate $\mathcal {C}_{\varepsilon }$ fits well with the model $\mathcal {C}_{\varepsilon } = \mathcal {C}_{\varepsilon,\infty } + \mathcal {D}/R_L^-$ for all levels of flow compressibility considered here, where $R_L^-$ is the generalized large-scale Reynolds number. The asymptotic value $\mathcal {C}_{\varepsilon,\infty }$ describes the total energy transfer flux, and decreases with increase of the flow compressibility, indicating non-universality of the dimensionless dissipation rate in compressible MHD turbulence. After introducing an empirically modified dissipation rate, the data from compressible cases collapse to a form similar to the incompressible MHD case depending only on the modified Reynolds number.

A theory of incompressible turbulent plane jets (TPJs) is proposed by advancing an improved boundary layer approximation over the limiting classical – retaining more terms in the momentum balance equations. A pressure deficit inside the jet (with respect to the ambient) must exist due to transverse turbulence (Miller & Comings, J. Fluid Mech., vol. 3, 1957, pp. 1–16; Hussain & Clarke, Phys. Fluids, vol. 20, 1977, pp. 1416–1426). Contrary to the universally accepted invariance of the total momentum flux $J_T(x)$ (non-dimensionalized by its inlet value) as a function of the streamwise distance $x$, we prove that $J_T(x) >1$ – a condition that all TPJs must satisfy; surprisingly, prior theories and most experiments do not satisfy this condition. This motivated us to apply Lie symmetry analysis with translational and dilatational transformations of the modified equations (incorporating $J_T>1$), which yields scaling laws for key jet measures: the mean streamwise and transverse velocities $U(x,y)$ and $V(x,y)$, the turbulence intensities, the Reynolds shear stress $-\rho \,\overline {u'v'}(x,y)$, the mean pressure $P(x,y)$, etc. Experiments satisfying $J_T(x)>1$ validate our predictions for all jet measures, including, among others, the profiles of $U$, $V$ and $-\rho \,\overline {u'v'}$. We further predict $U \sim x^{-0.24}$, $V \sim x^{-0.45}$, $-\rho \,\overline {u'v'}\sim x^{-0.69}$, the mass flux $Q_m \sim x^{0.55}$, and $J_T$ increases to approximately 1.5. Contrary to the classical linear jet spread, we find sublinear spread, with the jet half-width growing like $b(x)\sim x^{0.79}$, indicating a narrower jet. Our predictions differ notably from most results reported in the literature. These contradictions demand revisiting jet studies involving carefully designed facilities and boundary conditions, and highly resolved simulations.

The scale-dependent variability of convective velocities and structure inclination angles in wall-bounded turbulence was studied experimentally via space–time energy spectrum measurements. We found that the variability of convection velocities for large-scale motions (LSMs) decreased inversely with streamwise wavenumbers, and that the variability trend was not fully explained by earlier applications of Kraichnan's ‘random-sweeping’ model of turbulence that assumed perfect scale separation. By analytically extending the random-sweeping model to allow for a dominant large scale in the random-sweeping signal that can interact with other LSMs, we showed how scale interactions can explain the variability trend in convection velocities for LSMs. The variability in convection velocities was also shown to correlate with the scale-dependent inclination angles of coherent structures that were obtained via cross-spectral analysis. Large-scale motions tended to exhibit shallower inclination to the wall with increasing convection velocity, while small-scale motions and those far from the wall exhibited the reverse behaviour. We proposed that these two opposite relationships between inclination angle and convection velocity can be explained in terms of a balance between opposing effects of the mean shear and the coherent structure geometry. Descriptions and models of convection velocity variability effects are useful both for modelling turbulence spectra and explaining the geometry of coherent structures.

Analysing the dynamics of velocity gradients is useful for understanding various nonlinear turbulence processes. This work focuses on how the vibrational non-equilibrium of the constituent molecules in a gaseous medium affects the dynamics of the velocity gradient and the pressure Hessian tensors. We first derive the exact evolution equation of the pressure-Hessian tensor in the presence of the vibrational non-equilibrium process. Subsequently, we perform several direct numerical simulations of compressible isotropic turbulence, including the vibrational relaxation process therein. Using flow fields extracted from these simulations, we conduct several parametric studies over a range of the Damköhler number (ratio of the relevant fluid time scale to that of the mean vibrational relaxation process) and the initial ratio of the vibrational temperature to the mean local temperature. We find that a variation in the initial Damköhler number does influence the evolution of the pressure-Hessian and the velocity gradient tensors. As the vibrational relaxation process becomes more rapid (an increase in the value of the initial Damköhler number), it causes a decrease in the strength of the pressure-Hessian tensor and simultaneous suppression of dilatational fluctuations in the flow field. On the other hand, a variation in the initial value of the ratio of the vibrational temperature to the local temperature does not seem to affect the pressure-Hessian or the velocity gradient tensor. These findings are expected to aid in the development of closure models for the pressure-Hessian tensor in compressible flows under vibrational non-equilibrium conditions.

In 1945–1949, Lars Onsager made an exact analysis of the high-Reynolds-number limit for individual turbulent flow realisations modelled by incompressible Navier–Stokes equations, motivated by experimental observations that dissipation of kinetic energy does not vanish. I review here developments spurred by his key idea that such flows are well described by distributional or ‘weak’ solutions of ideal Euler equations. 1/3 Hölder singularities of the velocity field were predicted by Onsager and since observed. His theory describes turbulent energy cascade without probabilistic assumptions and yields a local, deterministic version of the Kolmogorov $4/5$th law. The approach is closely related to renormalisation group methods in physics and envisages ‘conservation-law anomalies’, as discovered later in quantum field theory. There are also deep connections with large-eddy simulation modelling. More recently, dissipative Euler solutions of the type conjectured by Onsager have been constructed and his $1/3$ Hölder singularity proved to be the sharp threshold for anomalous dissipation. This progress has been achieved by an unexpected connection with work of John Nash on isometric embeddings of low regularity or ‘convex integration’ techniques. The dissipative Euler solutions yielded by this method are wildly non-unique for fixed initial data, suggesting ‘spontaneously stochastic’ behaviour of high-Reynolds-number solutions. I focus in particular on applications to wall-bounded turbulence, leading to novel concepts of spatial cascades of momentum, energy and vorticity to or from the wall as deterministic, space–time local phenomena. This theory thus makes testable predictions and offers new perspectives on large-eddy simulation in the presence of solid walls.

Two simulations of turbulent Couette flows were performed at friction Reynolds numbers of 1000 and 2000 in a large box of dimensions $L_x=16{\rm \pi} h$, $L_y=2h$ and $L_z=6{\rm \pi} h$, where h is the semi-height of the channel. The study focuses on the differences in the intensity and scaling of turbulence at these two Reynolds numbers. The 2000 case showed a lack of a clear log layer with a higher value of the Von Kármán constant $\kappa$ than Poiseuille channels. The intensities were well-scaled in the buffer layer and below, with a second maximum of the streamwise intensity at approximately 350 wall units. Contrary to Poiseuille channels, the dissipation scales close to the wall in wall units. This fact can be attributed to the constant value of the derivative of the streamwise intensity in wall units. The intensities of the 2000 case showed remarkable differences compared with those at Reynolds number 1000 at the channel centre, likely due to the organization of large scales of the streamwise fluctuactions, $u$. These large scales were thought to be considered ‘infinite’. However, for the 2000 case, while all the structures have a width of $\ell _z \approx 6/8{\rm \pi} h$, their length varies from $\ell _x \approx 6{\rm \pi} h$ to $\ell _x \approx 16{\rm \pi} h$, which clearly contradicts the trends obtained in the past. This is a new effect that has not been reported for turbulent Couette flow and points to the uncertainty and sensitivity that is observed for certain statistical quantities.

Two sets of measurements utilizing hot-wire anemometry and oil-film interferometry for flat-plate turbulent boundary layers, exposed to various controlled adverse and favourable pressure gradients, are used to evaluate history effects of the imposed and varying free-stream gradients. The results are from the NDF wind tunnel at Illinois Tech (IIT) and the MTL wind tunnel at KTH, over the range $800 < Re_\tau < 22\,000$ (where $Re_{\tau }$ is the friction Reynolds number). The streamwise pressure-gradient parameter $\beta \equiv (-\ell /\tau _{w}) \cdot (\partial P_{e}/\partial x)$ varied between $-2 < \beta < 7$, where $\ell$ is an outer length scale for boundary layers equivalent to the half-height of channel flow and the radius of pipe flow, and is estimated for each boundary-layer profile; note that $\tau_w$ is the wall-shear stress and $P_e$ is the free-stream static pressure. Extracting from each profile the three parameters of the overlap region, following the recent work of Monkewitz & Nagib (J. Fluid Mech., vol. 967, 2023, p. A15) that led to an overlap region of combined logarithmic and linear parts, we find minimum history effects in the overlap region. Thus, the overlap region in this range of pressure-gradient boundary layers appears to be in ‘quasiequilibrium’.

Direct numerical simulations of channel flow and temporal boundary layer at a Reynolds number $Re_{\tau } = 1500$ are used to assess the scale-by-scale mechanisms of wall turbulence. From the peak of turbulence production embedded at the small scales of the near-wall region, spatially ascending reverse cascades are generated that move through self-similar eddies growing in size with the wall distance. These fluxes are followed by spatially ascending forward cascades through detached eddies thus reaching sufficiently small scales where eventually scale energy is dissipated. This phenomenology is shared by both boundary layer and channel flow and is recognized as a robust physical feature characterizing wall turbulence in general. Specific features related to the flow configuration are indeed identified in the outer region. In particular, the central region of channels is characterized by a generalized Richardson energy cascade where large scales are in equilibrium with small scales at different wall distances through a combined forward cascade and spatial flux. On the contrary, the interface region of boundary layers is characterized by an almost two-dimensional physics where spatially ascending reverse cascades sustain long and wide interface structures with a forward cascade that survives only in the wall-normal scales. The overall scenario consists in a variety of scale motions that while protruding from the turbulent core towards the external region, squeeze at the interface thus sustaining vertical shear in a thin layer. The observed multidimensional physics sheds light on the complex interactions between outer entrainment and near-wall self-sustaining mechanisms with possible repercussions for theories.

In this work, we employ well-established relations for compressible turbulent mean flows, including the velocity transformation and algebraic temperature–velocity (TV) relation, to systematically improve the algebraic Baldwin–Lomax (BL) wall model for high-speed zero-pressure-gradient air boundary layers. Any new functions or coefficients fitted by ourselves are avoided. Twelve published direct numerical simulation (DNS) datasets are employed for a priori inspiration and a posteriori examination, with Mach numbers up to 14 under adiabatic, cold and heated walls. The baseline BL model is the widely used one with semilocal scalings. Three targeted modifications are made. First, we employ a total-stress-based transformation (Griffin et al., Proc. Natl Acad. Sci. USA, vol. 118, issue 34, 2021, e2111144118) to the inner-layer eddy viscosity for improved scaling up to the logarithmic region. Second, we utilize the van Driest transformation in the outer layer based on the compressible defect velocity scaling. Third, considering the difficulty in modelling the rapidly varying and singular turbulent Prandtl number near the temperature peak in cold-wall cases, we design a two-layer strategy and use the TV relation to formulate the inner-layer temperature. Numerical results prove that the modifications take effect as designed. The prediction accuracy for mean streamwise velocity is notably improved for diabatic cases, especially in the logarithmic region. Moreover, a significant improvement in mean temperature is realized for both adiabatic and diabatic cases. The mean relative errors of temperature to DNS for all cases are down to 0.4 % in the logarithmic wall-normal coordinate and 3.4 % in the outer coordinate, around one-third of those in the baseline model.

Stirring a fluid through a Gaussian forcing at a vanishingly small Reynolds number produces a Gaussian random field, while flows at higher Reynolds numbers exhibit non-Gaussianity, cascades, anomalous scaling and preferential alignments. Recent works (Yakhot & Donzis, Phys. Rev. Lett., vol. 119, 2017, 044501; Khurshid et al., Phys. Rev. E, vol. 107, 2023, 045102) investigated the onset of these turbulent hallmarks in low-Reynolds-number flows by focusing on the scaling of the velocity gradients and velocity increments. They showed that the onset of power-law scalings in the velocity gradient statistics occurs at low Reynolds numbers, with the scaling exponents being surprisingly similar to those in the inertial range of fully developed turbulence. In this work we address the onset of turbulent signatures in low-Reynolds-number flows from the viewpoint of the velocity gradient dynamics, giving insights into its rich statistical geometry. We combine a perturbation theory of the full Navier–Stokes equations with velocity gradient modelling. This procedure results in a stochastic model for the velocity gradient in which the model coefficients follow directly from the Navier–Stokes equations and statistical homogeneity constraints. The Fokker–Planck equation associated with our stochastic model admits an analytic solution that shows the onset of turbulent hallmarks at low Reynolds numbers: skewness, intermittency and preferential alignments arise in the velocity gradient statistics through a smooth transition as the Reynolds number increases. The model predictions are in excellent agreement with direct numerical simulations of low-Reynolds-number flows.