1. Introduction
The study of highly compressible turbulent boundary layers (TBLs) is of major importance for high-speed turbulence research. Compressibility acts upon the flow by influencing the mean and fluctuating fields of thermodynamic quantities, which are, in turn, coupled to the momentum, promoting the energy exchange between kinetic and thermal fields. This poses several difficulties in the prediction of drag and wall heat transfer, which makes engineering design increasingly difficult as higher speeds are attained.
In recent decades, supersonic TBLs have been studied extensively and compared with their incompressible counterparts, mainly focusing on the prediction of drag assuming adiabatic walls; see, e.g., Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011a), Duan, Beekman & Martín (Reference Duan, Beekman and Martín2011) and Wenzel et al. (Reference Wenzel, Selent, Kloker and Rist2018). In fact, at supersonic speeds the wall temperature can be considered for practical purposes very close to the recovery temperature of the flow, implying a very low heat exchange at the wall. The recovery temperature indicates the temperature that is attained by the flow when it is brought to rest in a non-isentropic manner, defined as
\begin{equation} T_r=T_{\infty}\left(1+r \frac{\gamma-1}{2}M_{\infty}^2\right), \end{equation}
where 
$r = Pr^{1/3}$ is the recovery factor (Zhang, Duan & Choudhari Reference Zhang, Duan and Choudhari2018) and 
$Pr$ is the Prandtl number. However, in hypersonic boundary layers, the recovery temperatures are so high that the wall temperature is usually lower (Fernholz & Finley Reference Fernholz and Finley1980; Urzay & Di Renzo Reference Urzay and Di Renzo2021), generating large heat fluxes to the wall. This affects the flow dynamics in concurrency with the Mach number, enriching the physical effects that have to be accounted for when developing theoretical relations and reduced order models.
A renewed interest in hypersonic flight, along with the computational advancements that render direct numerical simulations (DNS) more feasible, sparked the attention on these problems, e.g. Zhang et al. (Reference Zhang, Duan and Choudhari2018) and Wenzel, Gibis & Kloker (Reference Wenzel, Gibis and Kloker2022), but there is still a lack of understanding of the individual effect of different flow parameters.
The framework of theoretical relations applied to compressible flow for mean velocity and fluctuating fields aims at mapping compressible profiles onto incompressible reference by taking into account variations of mean properties such as density and viscosity. When applied to the mean velocity field, these relations are called compressibility transformations, first introduced by Van Driest (Reference Van Driest1956) by accounting for mean density variations in the wall-normal velocity profile. Among the plethora of relations proposed in recent years, Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and Griffin, Fu & Moin (Reference Griffin, Fu and Moin2021) stand out as capable of efficiently collapsing velocity profiles even at high Mach numbers. Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) used a mixed physical and data-driven approach to determine the optimal parameters that define the weight of density and viscosity, whereas Griffin et al. (Reference Griffin, Fu and Moin2021) based their arguments on the total-stress equation, allowing for separate assumptions for the viscous sublayer and the log layer.
Theoretical relations have also been derived to describe the interaction between kinetic and thermal fields, classically referred to as strong Reynolds analogy (SRA). First proposed by Morkovin (Reference Morkovin1962), SRA establishes a framework based on the similarity between the momentum and total enthalpy equations, from which a direct proportionality between velocity and total enthalpy can be inferred. Under the more restrictive condition of wall adiabaticity, a set of relations coupling velocity and temperature can be derived for both mean and fluctuating fields, in which the temperature resembles a passive scalar field.
These relations have been validated extensively for adiabatic TBLs at different Mach numbers, e.g. Guarini et al. (Reference Guarini, Moser, Shariff and Wray2000), Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011a) and Wenzel et al. (Reference Wenzel, Selent, Kloker and Rist2018), although at hypersonic speeds (
$M_{\infty }>5$) discrepancies start to arise (Zhang et al. Reference Zhang, Duan and Choudhari2018). Subsequent extensions of the SRA accounting for diabatic walls have been proposed recently (e.g. Zhang et al. Reference Zhang, Bi, Hussain and She2014), which obtained promising results for different flow conditions, even when thermochemical effects are present (Di Renzo & Urzay Reference Di Renzo and Urzay2021; Passiatore et al. Reference Passiatore, Sciacovelli, Cinnella and Pascazio2021, Reference Passiatore, Sciacovelli, Cinnella and Pascazio2022).
A cold wall imposes a change in the sign of the mean temperature gradient near the wall, affecting the production of temperature fluctuations, which may result in a severe loss of similarity between velocity and temperature fields, a building block of SRA, clearly visible in instantaneous snapshots of turbulent structures (Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022; Zhang et al. Reference Zhang, Wan, Liu, Sun and Lu2022). However, these studies also noted that comparing cases with different Mach numbers at a fixed wall-to-recovery temperature ratio 
$T_w/T_r$ (
$<1$) resulted in vastly different near-wall dynamics for temperature fluctuations, in a way that cold cases at high 
$M_{\infty }$ seemed ‘more adiabatic’ than their low 
$M_{\infty }$ counterparts. Recently, other definitions of the wall temperature condition have been proposed, such as the diabatic parameter 
${\varTheta }=(T_w-T_{\infty })/(T_r-T_{\infty })$ (Zhang et al. Reference Zhang, Bi, Hussain and She2014) or the Eckert number 
$Ec=(\gamma -1)M_{\infty }^2 T_{\infty }/(T_r-T_w)$ (Wenzel et al. Reference Wenzel, Gibis and Kloker2022), which are capable of accounting for the indirect effect of Mach number on the wall temperature condition. Although progress has been made to incorporate the effects of compressibility and heat transfer on these relations, their individual influence is still not well understood. While compressibility effects induced by the increase in Mach number can be similar to a change in wall temperature condition (and vice versa) for certain mechanisms, such as redistribution of turbulent kinetic energy (TKE) (Duan & Martin Reference Duan and Martin2011), their relative role is still unclear in other aspects, such as separations of turbulent scales (Huang, Duan & Choudhari Reference Huang, Duan and Choudhari2022). In this regard, wall-cooling has been shown to reduce the separation between the large and small turbulence scales in hypersonic flows (Fan, Li & Pirozzoli Reference Fan, Li and Pirozzoli2022; Huang et al. Reference Huang, Duan and Choudhari2022), but the specific role of the Mach number is still debated. Furthermore, while wall-cooling has been shown to be dominant in regulating energy exchanges in the near-wall region (Fan et al. Reference Fan, Li and Pirozzoli2022), the effect of the Mach number is still not clear. These and other authors called upon the need for additional computations to determine their individual effects.
The aim of this study is to unveil the physical mechanisms that yield similarities and differences between the effect of compressibility and wall-cooling. To pursue this objective, an extensive DNS database consisting of 12 simulations of zero-pressure-gradient TBLs has been computed fixing the friction Reynolds number (
$Re_\tau \approx 450$), while spanning three Mach numbers 
$M_{\infty } = [2,4,6]$ and four diabatic parameters 
${\varTheta }=[0.25,0.5,0.75,1.0]$, going from extremely cold walls, 
${\varTheta }=0.25$, to adiabatic cases, 
${\varTheta }=1$. The database is discussed in the present paper and made available to the scientific community to be used for the development of simplified models for high-speed wall-bounded flows with strong heat flux.
The remainder of the paper is organised as follows. The numerical method and details on the simulation set-up are outlined in § 2. In § 3, a general visualisation of instantaneous velocity and temperature fields is given, describing the individual effect of Mach and wall temperature conditions on the flow dynamics and turbulent structures. Then, first-order statistics for mean velocity and temperature are presented in § 4, which also compares different SRA formulations. Finally, second-order statistics are presented in § 5, focusing on the effect of wall-cooling on thermal production, and its implications on velocity–temperature correlations.
2. Simulation parameters and computational set-up
The three-dimensional compressible Navier–Stokes equations are numerically solved for a viscous, heat-conducting gas
\begin{equation} \left. \begin{gathered} \frac{\partial \rho}{\partial t}+\frac{\partial (\rho u_j)}{\partial x_j}=0,\\ \frac{\partial (\rho u_i)}{\partial t}+\frac{\partial (\rho u_i u_j)}{\partial x_j}+\frac{\partial p}{\partial x_i}-\frac{\partial \sigma_{ij}}{\partial x_j}=0,\\ \frac{\partial (\rho E)}{\partial t}+\frac{\partial (\rho E u_j+pu_j)}{\partial x_j}-\frac{\partial (\sigma_{ij}u_i-q_j)}{\partial x_j}=0, \end{gathered} \right\} \end{equation}
where 
$\rho$ is the density, 
$u_i$ denotes the velocity component in the ith Cartesian direction (
$i=1,2,3$), 
$p$ is the thermodynamic pressure, 
$E=c_v T+u_i u_i/2$ is the total energy per unit mass and
\begin{equation} \sigma_{ij}=\mu \left( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}\delta_{ij}\right), \quad q_j={-}k\frac{\partial T}{\partial x_j} \end{equation}
represents the viscous stress tensor and the heat flux vector, respectively. The molecular viscosity 
$\mu$ is assumed to follow Sutherland's law
\begin{equation} \frac{\mu}{\mu_{\infty}}=\left(\frac{T}{T_{\infty}}\right)^{1/2} \frac{1+C/T_{\infty}}{1+C/T}, \end{equation}
where 
$C=110.4$ K and 
$T_{\infty }=220.0$ K, representing the typical conditions that are met in the stratosphere. The thermal conductivity 
$k$ is related to the viscosity by the expression 
$k=c_p \mu /Pr$, where 
$c_p$ is the specific heat at constant pressure and the Prandtl number is 
$Pr=0.72$. The thermodynamic variables are correlated to each other by means of the equation of state for a calorically perfect gas. This choice was also assumed for cases at 
$M_{\infty }=6$, after having verified that by introducing a dependence of specific heat with temperature 
$c_p=f(T)$ differences in all statistics were negligible. Moreover, gas dissociation effects are also not expected in the present database, according to the observations of Passiatore et al. (Reference Passiatore, Sciacovelli, Cinnella and Pascazio2022) who observed negligible effects with 
$T_w=1800$ K (our highest imposed value is 
$T_w=1640$ K for 
$M_{\infty }=6$). The system of equations is solved on a Cartesian grid using the in-house code STREAmS (Bernardini et al. Reference Bernardini, Modesti, Salvadore and Pirozzoli2021, Reference Bernardini, Modesti, Salvadore, Sathyanarayana, Della Posta and Pirozzoli2023), which has been validated extensively in numerous canonical configurations (Bernardini, Pirozzoli & Grasso Reference Bernardini, Pirozzoli and Grasso2011; Bernardini & Pirozzoli Reference Bernardini and Pirozzoli2011b; Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022). Convective terms are discretised using high-order, energy-preserving schemes applied in shock-free regions, while a fifth-order shock capturing scheme (WENO) is applied when strong compressions are identified by the Ducros sensor (Ducros et al. Reference Ducros, Ferrand, Nicoud, Weber, Darracq, Gacherieu and Poinsot1999). Diffusive terms are discretised using a locally conservative formulation (De Vanna et al. Reference De Vanna, Benato, Picano and Benini2021), expanded to Laplacian form to ensure finite molecular dissipation at all resolved wavelengths. The solver takes advantage of a multi-GPU architecture by means of the CUDA Fortran paradigm. The domain is a rectangular box of length 
$L_x=100 \delta _{in}$, 
$L_y=15 \delta _{in}$, 
$L_z=9 \delta _{in}$, where 
$\delta _{in}$ is the boundary layer thickness at the inflow station, based on the 
$99\,\%$ of the freestream velocity 
$u_{\infty }$ (which is referred for other stations as 
$\delta _{99}$). For all cases, the domain size in terms of the boundary-layer thickness measured in the proximity of the outflow is larger than 
$L_x / \delta _{99}= 39$, 
$L_y/ \delta _{99} =5.8$, 
$L_z / \delta _{99}=3.5$, respectively. For each spatial direction, the number of computational points employed for all cases is 
$N_x=5120$, 
$N_y=320$, and 
$N_z=512$. We verified that selected the grid satisfies the resolution requirement of 
$\Delta x_i / \eta < 5$ throughout the entirety of the boundary layer thickness for all cases, being 
$\eta$ the Kolmogorov length scale. Periodic boundary conditions are enforced in the spanwise direction, purely non-reflecting boundary conditions are employed for the outflow and the top boundary, and unsteady characteristic boundary conditions are used at the bottom wall (Poinsot & Lele Reference Poinsot and Lele1992), where an isothermal wall temperature condition is enforced. The recycling–rescaling procedure (Pirozzoli, Bernardini & Grasso Reference Pirozzoli, Bernardini and Grasso2010) is applied at the inflow to reach a fully developed state, the recycling length being placed at a distance of 
$80 \delta _{in}$ from the inlet, ensuring a complete decorrelation of the fluctuations between the recycling station and the inflow plane (Morgan et al. Reference Morgan, Larsson, Kawai and Lele2011).
Table 1 summarises the flow conditions and grid resolutions for each run, where 
$M_{\infty }$ is the freestream Mach number and 
$Re_{\tau }$ is the friction Reynolds number, defined as the ratio between the boundary layer thickness 
$\delta _{99}$ and the viscous length scale 
$\delta _{\nu }=\nu _w/u_{\tau }$, where 
$u_{\tau }=\sqrt {\tau _w/\rho _w}$ is the friction velocity, 
$\tau _w$ is the mean wall shear stress and 
$\nu _w$ is the kinematic viscosity at the wall. Here 
$\Delta x^+=\Delta x / \delta _{\nu }$ and 
$\Delta z^+=\Delta z / \delta _{\nu }$ are the uniform grid spacings in the streamwise and spanwise directions and 
$\Delta y^+=\Delta y / \delta _{\nu }$ represents the non-uniform wall-normal grid spacing (the wall and edge values are reported). In the wall-normal direction, the stretching function of Pirozzoli & Orlandi (Reference Pirozzoli and Orlandi2021) is employed, which provides a more favourable scaling of the number of grid points with the Reynolds number. Furthermore, this function has the natural property of yielding precisely constant resolution in terms of the local Kolmogorov length scale 
$\eta$ in the outer part of the wall layer while maintaining a uniform near-wall spacing.
Summary of parameters for DNS study. Grid spacings are given in wall-units according to the stations selected in table 2. The values of 
$\Delta y^+_{w}$ and 
$\Delta y^+_{edge}$ refer to the wall-normal spacing at the wall and at the boundary layer edge, respectively. The range of 
$Re_{\tau }$ is representative of the statistical growth of the boundary layer's thickness along 
$x$. Ranges reported for 
$Re_{\tau }$ refer to the second half of the domain, while all cases share the same inflow friction Reynolds number of 
$Re_{\tau }=250$.
The present database is composed of a total of 12 simulations, spanning three Mach numbers 
$M_{\infty } = [2,4,6]$ and four diabatic parameters 
${\varTheta }=[0.25,0.5,0.75,1.0]$ (see table 1). We stress that the choice of using the diabatic parameter throughout this study is made in order to recover a similar degree of wall-cooling across different Mach numbers. In fact, the goal of a suitable parameter is not to be independent of 
$M_\infty$, but to incorporate it in order to have ‘the same integral behaviour between different cases, regardless of whether its variation is caused by the change of the Mach number or of the wall temperature’ (Wenzel et al. Reference Wenzel, Gibis and Kloker2022). This is the rationale with which Wenzel et al. (Reference Wenzel, Gibis and Kloker2022) argued that the Eckert number 
$Ec=(\gamma -1)M_{\infty }^2 T_{\infty }/(T_r-T_w)$ represents a more suitable option than the conventional 
$T_w/T_r$ ratio. The Eckert number happens to be directly related to the diabatic parameter 
${\varTheta }=(T_w-T_{\infty })/(T_r-T_{\infty })$ proposed by Zhang et al. (Reference Zhang, Bi, Hussain and She2014), since it can be shown that 
$Ec = 2/[r(1-{\varTheta })]$. This parameter shows more clearly the improvement over the conventional ratio 
$T_w/T_r$, showing that 
$T_{\infty }$ needs to be subtracted from both 
$T_w$ and 
$T_r$ to compare only the 
$\Delta T$ that is recovered when the flow is brought at rest, being the only one responsible for kinetic–internal energy exchange. In this study, we choose to use the diabatic parameter 
${\varTheta }$ over 
$Ec$ given its simplicity, but we also report the latter in table 1. Table 2 summarises the boundary layer parameters at selected locations where turbulence statistics are gathered. All the simulations were initially run for about 600 
$\delta _{in}/ u_{\infty }$, a time sufficient to achieve a fully developed turbulent condition before starting to collect statistics. The time averaging period for all cases exceeds 
$1000 \delta _{in}/ u_{\infty }$ (
$28 \delta _{99}/ u_{\tau }$) and approximately 1000 flow samples, equally spaced in time, were considered.
Boundary layer properties averaged at the selected station. Here 
$Re_{\tau }=\rho _w u_{\tau } \delta /\mu _w$; 
$Re_{\theta }=\rho _{\infty }u_{\infty }\theta /\mu _{\infty }$; 
$Re_{\delta _2}=\rho _{\infty }u_{\infty }\theta /\mu _{w}$; 
$Re_{\tau }^*=\sqrt {\rho _{\infty } \tau _w}\delta /\mu _{\infty }$; 
$H=\delta ^*/\theta$ (
$\delta ^*$ and 
$\theta$ are the boundary layer displacement and momentum thickness, respectively); 
$B_q=q_w/(\rho _w C_p u_{\tau }T_w)$ and 
$C_f=\tau _w/(1/2 \rho _{\infty } u_{\infty }^2)$ are the non-dimensional wall heat transfer 
$q_w=-\bar {k} \ \partial \tilde {T}/\partial y$ and skin friction coefficient 
$\tau _w = \bar {\mu } \partial \tilde {u} / \partial y$, respectively.
Throughout this study, we use the symbols 
$u$, 
$v$ and 
$w$ to denote the streamwise, wall-normal and spanwise velocity components and the decomposition of any variable is conducted using either the standard Reynolds decomposition (
$f=\bar {f}+f'$) or the density-weighted (Favre) representation (
$f=\tilde {f}+f^{\prime \prime }$), being 
$\tilde {f}=\bar {\rho f}/\bar {\rho }$. Here, the averaging is conducted using multiple samples and along the periodic direction 
$z$.
3. Instantaneous visualisation
To highlight the emerging features of the flow in a qualitative way, we selected the two extreme cases with regard to the wall-cooling condition, 
${\varTheta }=0.25$ and 
${\varTheta }=1.0$, for each Mach number in our database. Figure 1 shows a portion of wall-normal 
$x$–
$y$ planes coloured with the instantaneous density, whose variability is a clear sign of the degree of compressibility. The effect of Mach number is clearly apparent for all cases moving from top to bottom with a decrease of the minimum value of density and an increase of the general level of acoustic disturbances, generated in the boundary layer and emanated towards the far field. However, a stronger wall-cooling (lower 
${\varTheta }$, figure 1a,c,e) attenuates this effect, since lower wall temperatures generate higher density fields in the near-wall region.
Instantaneous density in wall-normal slices (
$x$–
$y$ plane), with a window size of 
$\Delta x =20\delta$–
$30 \delta$ and 
$\Delta y =0\delta$–
$2 \delta$. Here, all Mach numbers are shown while the two extremes are chosen with regard to wall-cooling (
${\varTheta }=0.25$ and 
${\varTheta }=1.0$): (a) 
$M_{\infty }=2,\ {\varTheta }=0.25$; (b) 
$M_{\infty }=2,\ {\varTheta }=1.0$; (c) 
$M_{\infty }=4,\ {\varTheta }=0.25$; (d) 
$M_{\infty }=4,\ {\varTheta }=1.0$; (e) 
$M_{\infty }=6,\ {\varTheta }=0.25$; ( f) 
$M_{\infty }=6,\ {\varTheta }=1.0$.
The intensity of wall-cooling strongly affects the coupling between velocity and temperature fluctuations, especially in the near-wall region. This is apparent in figure 2, which compares these quantities in wall-parallel slices located at approximately 
$y^*\approx 10$, representing the onset turbulence activity after the viscous sublayer. Here, 
$y^*=y/ \delta _{\nu,SL}$ is the semilocal scaled wall-normal coordinate, with 
$\delta _{\nu,SL}= \bar {\nu }/\sqrt {\tau _w/\bar {\rho }}$. The chosen 
$x$–
$z$ planes are centred at the selected stations of table 1 spanning a window of 
$\Delta x^*=4000$ and 
$\Delta z^*=600$. Velocity fluctuations 
$\sqrt {\bar {\rho }} u^{\prime }/\sqrt {\tau _w}$ are scaled according to the Morkovin's transformation (Morkovin Reference Morkovin1962) (also used in § 5.1), which enables comparison across different Mach numbers and wall temperature values by accounting for the variation of the mean properties of the flow. In other words, velocity fluctuations are scaled by the semilocal friction velocity 
$u_{\tau,SL}=\sqrt {\tau _w/\bar {\rho }}$, which differs from the conventional one by employing the mean density 
$\bar {\rho }$ instead of the wall density 
$\rho _w$. Temperature fluctuations 
$\bar {\rho } T^{\prime }/ (R \tau _w)$ are scaled in a similar fashion, assuming 
$\tau _w$ as proper parameter to scale pressure fluctuations, then 
$\tau _w/(R \bar {\rho })$ can be used to scale temperature (for further details refer to § 5.2). A general look at the velocity fluctuations shows the presence of near-wall streaks for all cases, representative of the near-wall self-sustaining cycle of turbulence. Similar values of intensities appear across all cases. This result is not observed for temperature fluctuations, where cold cases (
${\varTheta }<1$) show reduced intensity and a clear breakdown of elongated streaks, appearing more isotropic when compared with their adiabatic counterpart. Adiabatic cases maintain a streaky pattern, which shows a clear coupling with the velocity field. It is worth noting that the similarity between cold cases across different Mach numbers (figure 2a,c,e,g) is attained using the same value of 
${\varTheta }$, while effectively the wall-to-recovery temperature ratio 
$T_w/T_r$ varies. Clearly different patterns emerge matching the latter, as visible in Cogo et al. (Reference Cogo, Salvadore, Picano and Bernardini2022). Although this behaviour is further discussed in the following sections by analysing temperature fluctuations and thermal production profiles, these qualitative results are consistent with the discussion of Wenzel et al. (Reference Wenzel, Gibis and Kloker2022), which states that the same general behaviour due to the effect of wall-cooling is expected when comparing flows with the same Eckert number (or diabatic parameter 
${\varTheta }$).
Temperature fluctuations 
$\bar {\rho } T'/ \tau _w$ (top) and streamwise velocity fluctuations 
$\sqrt {\bar {\rho } } u'/ \sqrt {\tau _w}$ (bottom) in wall-parallel slices (
$x$–
$z$ plane) selected at 
$y^*\approx 10$. Here, all Mach numbers are shown while the two extremes are chosen with regard to wall-cooling (
${\varTheta }=0.25$ and 
${\varTheta }=1.0$): (c) 
$M_{\infty }=2,\ {\varTheta }=0.25$; (d) 
$M_{\infty }=2,\ {\varTheta }=1.0$; (e) 
$M_{\infty }=4,\ {\varTheta }=0.25$; ( f) 
$M_{\infty }=4,\ {\varTheta }=1.0$; (g) 
$M_{\infty }=6,\ {\varTheta }=0.25$; (h) 
$M_{\infty }=6,\ {\varTheta }=1.0$. Here, 
$x_{ref}$ is the streamwise location of the selected station.
4. Mean flow statistics
In this section, we present the wall-normal profiles of averaged quantities such as velocity and temperature, selected at stations listed in table 2. In the framework of compressibility transformations, we consider the classical relation of Van Driest (Reference Van Driest1951) (which has been the standard for several decades and widely employed in wall modelling) and the recent transformations of Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and Griffin et al. (Reference Griffin, Fu and Moin2021). In general, we find that the latter developments have a wide range of applicability over our database, while Van Driest (Reference Van Driest1951) transformation shows poor accuracy as Mach number and wall-cooling increase. The results confirm the overall behaviour noted by Cogo et al. (Reference Cogo, Salvadore, Picano and Bernardini2022), and the interested reader can find them in Appendix A.
Figures 3(a)–3(e) show the mean temperature profiles throughout the height of the boundary layer and in the near-wall region, respectively. In particular, all profiles are scaled considering their incremental variation and their relative intensity with respect to 
$T_w$. This is done in order to highlight the effects of wall-cooling on the peak intensity and location. As expected, the adiabatic wall temperature greatly increases with the Mach number, while enhanced wall-cooling (lower 
${\varTheta }$) forces the mean temperature profiles to slant towards lower wall temperatures (
$T_w< T_r$). The combination of these two conditions imposes a change in the sign of temperature gradient near the wall, which is necessary to adjust to a wall temperature lower than the recovery value. Thus, a local peak arises, whose prominence and location are directly connected to the phenomenon of aerodynamic heating, generating a net heat flux from the flow to the solid boundary. Local temperature peaks are marked in figures 3(b)–3(e) with dots. An increase in the Mach number generates more intense gradients and higher peak temperatures for non-adiabatic cases, enhancing aerodynamic heating. However, the wall-normal position of the peaks seems to be mainly affected by the change 
${\varTheta }$, and weakly dependent on the Mach number. This is apparent in figure 3( f), which shows a progressive departure from the wall of the peak location as the wall-cooling increases, with a mild downward shift at high Mach numbers. As anticipated in § 3, the position of the local maximum of the temperature profile has major implications in the generation of temperature fluctuations, which affect both their overall intensity and their spatial organisation (breakdown of near-wall streaks). Actually, the departure from a monotonic adiabatic profile to increasingly prominent local peaks of mean temperature profiles prevents the formation of organised temperature streaks that are generated by thermal production (see § 5.2).
(a) Mean temperature profiles for all cases of table 1 as a function of the wall-normal coordinate 
$y/ \delta _{99}$. (b–e) Mean temperature profiles and relative peaks as a function of the wall-normal coordinate 
$y^*$ scaled with 
$T_w$: (b) 
${\varTheta }=0.25$; (c) 
${\varTheta }=0.5$; (d) 
${\varTheta }=0.75$; (e) 
${\varTheta }=1$. ( f) Wall-normal position of mean temperature peaks as a function of the wall-cooling 
${\varTheta }$ parameter.
4.1. Reynolds analogy
In the framework of the Reynolds analogy, we discuss the coupling between velocity and temperature with a particular focus on the validity of theoretical relations across the present database. The relation between mean temperature and velocity can be approximated by a quadratic law, as apparent by the classical relation of Walz (Reference Walz1969)
\begin{equation} \frac{\bar{T}}{T_{\infty}}=\frac{T_w}{T_{\infty}}+\frac{T_r-T_w}{T_{\infty}}\frac{\bar{u}}{U_{\infty}}+ \frac{T_{\infty}-T_r}{T_{\infty}}\left(\frac{\bar{u}}{U_{\infty}}\right)^2. \end{equation}More recently, Zhang et al. (Reference Zhang, Bi, Hussain and She2014) improved this relation in order to account for high Mach number and large heat fluxes:
\begin{equation} \frac{\bar{T}}{T_{\infty}}=\frac{T_w}{T_{\infty}}+\frac{T_{rg}-T_w}{T_{\infty}}\frac{\bar{u}}{U_{\infty}}+ \frac{T_{\infty}-T_{rg}}{T_{\infty}}\left(\frac{\bar{u}}{U_{\infty}}\right)^2, \end{equation}
where 
$T_{rg}=T_{\infty }+r_g U_{\infty }^2/(2 c_p)$ and 
$r_g=2 c_p (T_w-T_{\infty })/U_{\infty }^2-2 \, Pr \, q_w/(U_{\infty } \tau _w)$. As reported in Appendix B, our database confirms the aptness of the latter expression, with minor deviations for the most challenging case (M6T025). For engineering design purposes, the value of 
$r_g$ can be difficult to evaluate given its dependence on the wall temperature 
$T_w$ and the ratio of the wall heat flux 
$q_w$ and the wall shear stress 
$\tau _w$. Following the discussion of Zhang et al. (Reference Zhang, Bi, Hussain and She2014), the Reynolds analogy factor 
$s$ comes into play to greatly simplify the calculation, since 
$r_g$ can be rewritten in terms of 
$s$
\begin{equation} r_g = r[s Pr + (1-s Pr) {\varTheta}] \end{equation}
with 
$s$ being defined as
\begin{equation} s = \frac{2 C_h}{C_f} = \frac{q_w u_{\infty}}{\tau_w c_p (T_w-T_r)}, \end{equation}
where 
$C_f=\tau _w/(1/2 \rho _{\infty } u_{\infty }^2)$ is the skin friction coefficient and 
$C_h=q_w/(\rho _{\infty } u_{\infty } c_p (T_w-T_r))$ the Stanton number. The simplification consists of the fact that several authors (Duan, Beekman & Martín Reference Duan, Beekman and Martín2010; Zhang et al. Reference Zhang, Bi, Hussain and She2014; Wenzel et al. Reference Wenzel, Gibis, Kloker and Rist2021) identified the term 
$s Pr$ to be an empirical constant around the value of 
$0.8\pm 0.03$ (data fitting of Zhang et al. Reference Zhang, Bi, Hussain and She2014) over several different flow cases, meaning that only 
$T_w$ would be needed to be evaluated to compute 
$r_g$.
Figure 4 reports the computed values of 
$s Pr$ in our database showing a good agreement with the fit of Zhang et al. (Reference Zhang, Bi, Hussain and She2014). A slight decreasing trend with 
${\varTheta }$ can be observed, and it is interesting to note that at a given 
${\varTheta }$ the values appear to be independent of 
$M_{\infty }$.
Reynolds analogy factor 
$s=2 C_h/C_f Pr$ as a function of the diabatic parameter 
${\varTheta }$ for different Mach numbers. The grey band refers to the data fitting of 
$0.8\pm 0.03$ of Zhang et al. (Reference Zhang, Bi, Hussain and She2014).
The data reported in figure 4 have a mean value and standard deviation of 
$0.78\pm 0.03$, which is close to the value reported by (Zhang et al. Reference Zhang, Bi, Hussain and She2014). By approximating 
$r_g$ in (4.3) with the mean value of 
$sPr$ and comparing it with DNS data, we obtain a maximum error of 
$5\,\%$ (for the case M6T025), which can be considered acceptable for engineering purposes.
Another important set of theoretical relations that couple the thermodynamic and kinetic fluctuating fields is given by the SRA (Morkovin Reference Morkovin1962). This set of relations consists have been modified over the years to account for finite heat flux at the wall and remove wall temperature dependence (Gaviglio Reference Gaviglio1987; Huang, Coleman & Bradshaw Reference Huang, Coleman and Bradshaw1995; Zhang et al. Reference Zhang, Bi, Hussain and She2014) (HSRA). Our database shows that the most recent developments clearly improve its accuracy and the insensitivity to the freestream Mach number and wall temperature condition, with only slight deviations at the edge of the boundary layer (see Appendix B). Moreover, results also highlight the ability of the diabatic parameter 
${\varTheta }$ in recovering the same behaviour in terms of wall-cooling across different Mach numbers.
5. Fluctuation statistics
5.1. Velocity fluctuations and length scales
The distribution of velocity fluctuation intensities and Reynolds shear stress is reported in figures 5(a,c,e) and 6(a,c,e,g), using the classical transformation of Morkovin (Reference Morkovin1962):
\begin{equation} \left(u_i^*\right)^2=\frac{\widetilde{u_i^{\prime \prime 2}}}{u_\tau^2} \frac{\bar{\rho}}{\bar{\rho}_w}, \quad(u v)^*=\frac{\widetilde{u^{\prime \prime} v^{\prime \prime}}}{u_\tau^2} \frac{\bar{\rho}}{\bar{\rho}_w}. \end{equation}
The profiles are shown as a function of the wall-normal distance in semilocal scaling 
$y^*$ (Huang et al. Reference Huang, Coleman and Bradshaw1995), considering its ability to collapse compressible profiles of different Mach numbers and wall temperature conditions, in particular with respect to the peak positions (Zhang et al. Reference Zhang, Duan and Choudhari2018, Reference Zhang, Wan, Liu, Sun and Lu2022). This choice facilitates the comparison across different cases of our database in terms of wall-normal location, even though differences in their intensity are preserved. Figures 5(b,d, f) and 6(b,d, f,h) show the corresponding TKE budget terms (with 
$k=\widetilde {u_i^{\prime \prime } u_i^{\prime \prime }}/2$ being the TKE) according to the derivation of Zhang et al. (Reference Zhang, Duan and Choudhari2018), where 
$P$ is the production term, 
$TT$ represents the turbulent transport, 
$\varPi$ includes the pressure diffusion and dilatation, 
$-\phi$ is the viscous dissipation and 
$D$ is the viscous diffusion. For these results, semilocal scaling is also employed in the normalisation of budget terms (refer to Zhang et al. Reference Zhang, Duan and Choudhari2018) and for the wall-normal distance 
$y^*$, enabling a good collapse between different profiles (Zhang et al. Reference Zhang, Duan and Choudhari2018; Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022).
Semilocal-scaled turbulent velocity fluctuations (a,c,e) and turbulent kinetic budget (b,d, f) as functions of the wall-normal distance 
$y^*$. Here, different diabatic parameters 
${\varTheta }$ are compared at a given Mach number 
$M_{\infty }$: (a,b) 
$M_{\infty }=2$; (c,d) 
$M_{\infty }=4$; (e, f) 
$M_{\infty }=6$.
Semilocal-scaled turbulent velocity fluctuations (a,c,e,g) and turbulent kinetic budget (b,d, f,h) as function of the wall-normal distance 
$y^*$. Here, different Mach number 
$M_{\infty }$ are compared at a given diabatic parameter 
${\varTheta }$: (a,b) 
${\varTheta }=0.25$; (c,d) 
${\varTheta }=0.5$; (e, f) 
${\varTheta }=0.75$; (g,h) 
${\varTheta }=1.0$.
The effect of wall-cooling on velocity fluctuations, shown in figure 5(a,c,e), is apparent as an increase in the peak of the streamwise component located at 
$y^*\approx 15$ that is more prominent at high Mach numbers. In contrast, the spanwise component of highly cooled cases shows the opposite behaviour, being reduced in intensity compared with the adiabatic reference. This implies an increase in the anisotropy of normal components of Reynolds stresses in the near-wall region, which is discussed in more detail at the end of this section. The semilocal scaling provides an excellent collapse of the peak positions for all cases, preventing the outward shift that is present for cold cases when plotted in wall units (not shown). This is also true for the position of the turbulent production peak (figure 5b,d, f), which would move farther from the wall if displayed in wall units. In general, the effect of wall-cooling on the TKE budget is marked in the very near-wall region, especially at high Mach numbers, while all profiles progressively collapse in the outer layer. The effect of the Mach number on velocity fluctuations is reported in figure 6(a,c,e,g), where an increase of the streamwise component peak with the Mach number is apparent, while the other normal components intensities are observed to weakly decrease until 
$y^*<40$. Unlike the wall-cooling effect, all normal components increase in the log layer as 
$M_{\infty }$ increases. We note that this effect could be reduced at the boundary layer edge by matching the semilocal friction Reynolds number 
$Re_{\tau }^*$ in place of the conventional definition (see table 1), which would allow all profiles to collapse when 
$y^* \approx Re_{\tau }^*$. In fact, 
$Re_{\tau }^*$ has been shown by several authors to better incorporate compressibility and wall-cooling effects on the separation of scales in highly compressible flows (e.g. Griffin et al. Reference Griffin, Fu and Moin2021; Hirai, Pecnik & Kawai Reference Hirai, Pecnik and Kawai2021). However, in the present study we preferred to match wall-related quantities (i.e. 
$Re_{\tau }$) and discuss the influence of 
$M_{\infty }$ and 
${\varTheta }$ on the separation of turbulent scales through the boundary layer thickness.
This suggests that compressibility acts in the direction of increasing the scale separation in the outer layer, while wall-cooling has the opposite effect (Fan et al. Reference Fan, Li and Pirozzoli2022). For all values of 
${\varTheta }$, the TKE budget (figure 6b,d, f,h), shows an increase of the production term 
$P$ in the buffer and log layers as the Mach number increases and a corresponding decrease of diffusion 
$D$ and turbulent transport 
$TT$ in the same regions, consistently with Cogo et al. (Reference Cogo, Salvadore, Picano and Bernardini2022), who noted the presence of this effect also in the outer region at higher 
$Re_{\tau }$. While the effect of wall-cooling on the TKE budget seems confined in the near-wall region, the influence of Mach number is more prominent after the peak of production and throughout the log layer. Further insights on the mechanism of redistribution of TKE in the near-wall region can be gained by looking at the ratio between the streamwise component of the pressure–strain term and the streamwise component of turbulent production (Duan et al. Reference Duan, Beekman and Martín2010):
\begin{equation} \mathcal{R} = \left.\left(\overline{p^{\prime} \frac{\partial u^{\prime\prime}}{\partial x}}\right) \right/\left(\overline{\rho u^{\prime\prime} v^{\prime\prime}}\frac{\partial \tilde{u} }{\partial y}\right), \end{equation}which is a measure of the energy transfer from the streamwise velocity fluctuations to the others.
The role of the pressure–strain term in increasing turbulence anisotropy was also noted for other flows (e.g. Foysi, Sarkar & Friedrich Reference Foysi, Sarkar and Friedrich2004). To gauge the respective effects of the Mach number and the wall temperature condition, figure 7 compares 
$-\mathcal {R}$ for cases at 
${\varTheta }=1.0$ (figure 7a) and at 
$M_{\infty }=6$ (figure 7b). In figure 7(a), profiles of 
$-\mathcal {R}$ are reduced in magnitude as compressibility increases, with greater intensity farther from the wall.
Ratio of streamwise components of pressure–strain and turbulent production terms for cases at (a) 
${\varTheta }=1.0$ and (b) 
$M_{\infty }=6$, as functions of the wall-normal distance in semilocal units.
This is consistent with the less efficient redistribution of TKE discussed previously, and is attributed to the absence of a solenoidal condition for the velocity field for highly compressible cases preventing an efficient energy transfer between velocity components. Looking at panel 7(b), we observe that wall-cooling acts similarly to an increase of compressibility, strongly decreasing the profiles of 
$-\mathcal {R}$, but its effect is localised in the near-wall region and strongly reduced after the buffer layer.
We attribute this effect to a localised stratification of flow properties in the near-wall region. As wall-cooling is increased and the location of the mean temperature peak approaches the buffer layer, the flow above and below the peak location is relatively colder and denser. This is true for all Mach numbers in our database (although with different intensities), since the temperature peak location remains unaffected (see figure 3f). However, the effect is marked when the peak intensity is higher (i.e. 
$M_{\infty }=6$). This is visible in figure 8, where case M6T025 of panel (a) exhibits a high-temperature substrate in the near-wall region surrounded by a colder flow above and below, which is in direct contrast with the adiabatic case shown in panel (b).
Instantaneous temperature contours in wall-normal slices (
$x$–
$y$ plane) at 
$M_{\infty }=6$: (a) case M6T025 (cold wall), 
${\varTheta }=0.25$; (b) case M6T100 (adiabatic), 
${\varTheta }=1.0$.
This localised stratification forces turbulent fluctuations to be active almost only in the streamwise direction, while the other components tend to be suppressed. This effect is quantified in figure 9 by showing the barycentric map of Banerjee et al. (Reference Banerjee, Krahl, Durst and Zenger2007), which shows principal components of turbulence anisotropy. The invariant map is composed of three limiting states, one-component (
$x_{1c}$), axisymmetric two-component (
$x_{2c}$) and isotropic (
$x_{3c}$), which are representative of the relative strengths of the fluctuating velocity components. Looking at figure 9(a), we note that the cusp point, which coincides approximately with the peak of velocity fluctuations in the buffer layer, shifts towards a one-component behaviour (
$x_{1c}$) as 
$M_{\infty }$ increases. This effect is strongly enhanced by wall-cooling, figure 9(b), which further promotes the one-dimensional state of the flow.
Barycentric map of Banerjee et al. (Reference Banerjee, Krahl, Durst and Zenger2007) for cases at (a) 
${\varTheta }=1.0$ and (b) 
$M_{\infty }=6$. Each point in the trajectories represents a different wall-normal location within the boundary layer up to the edge.
We note that although this effect resembles a promotion of compressibility, as noted by several authors (Duan & Martin Reference Duan and Martin2011; Chu, Zhuang & Lu Reference Chu, Zhuang and Lu2013), the underlying mechanism is strongly different and relevant only when 
$M_{\infty }$ is high. In fact, different wall-cooling and compressibility signatures are clearly noted for other effects, such as their effect on scale separation and their region of influence through the boundary layer.
To provide further insights on these differences, we analyse the characteristic turbulent lengths. We consider the length scale characterising large eddies as 
$L=\bar {\rho } k^{3/2}/ \phi$ (Pope Reference Pope2000), and the Kolmogorov length scale 
$\eta =[(\bar {\mu } / \bar {\rho })^3 /(\phi / \bar {\rho })]^{0.25}$ for the smallest ones, with 
$\phi$ being the local dissipation rate of TKE. The ratio of these two scales, reported in figure 10, measures the separation between large and small scales, which in our discussion can be ascribed to the effect of 
$M_{\infty }$ and 
${\varTheta }$ numbers (since 
$Re_{\tau }$ is fixed). In agreement with previous observations, figure 10(a) shows that the separation of scales in the outer layer increases with the Mach number, while the opposite behaviour is found reducing the diabatic parameter 
${\varTheta }$, see figure 10(b). The insets in figure 10(a,b) show the individual change of 
$L^+=L/ \delta _{\nu }$ and 
$\eta ^+=\eta / \delta _{\nu }$, revealing that 
$M_{\infty }$ and 
${\varTheta }$ strongly affect the Kolmogorov length 
$\eta ^+$, with a minor impact on large scales 
$L^+$, influencing the separation of large to small scales 
$L/ \eta$ in the outer layer.
Ratio of integral length scale 
$L$ and Kolmorov scale 
$\eta$ for cases at (a) 
${\varTheta }=1.0$ and (b) 
$M_{\infty }=6$, as function of the wall-normal distance 
$y^+$. The inset shows separately 
$L^+$ and 
$\eta ^+$, normalised with the viscous length 
$\delta _{\nu }$.
On this aspect, we remark that while an increase in compressibility, i.e. 
$M_{\infty }$, reduces the Kolmogorov length, the opposite holds decreasing the wall temperature, i.e. 
${\varTheta }$.
The variation of 
$L/\eta$ in the outer layer is effectively captured by the change of 
$Re_{\tau }^* = \mu _w / \mu _{\infty } \sqrt {\rho _{\infty } / \rho _w} Re_{\tau }$ (see table 2), which better account for density and viscosity variations in the outer layer. It should be noted, however, that the definition of a single similarity parameter among different flow cases concerning the scale separation is prevented by the strong change of flow properties across the boundary layer. In particular, while 
$Re_{\tau }$ essentially regulates the outer–inner scale separation, i.e. 
$L^+$, 
$Re_{\tau }^*$ controls the large-small scale separation in the outer layer, i.e. 
$L/ \eta$. These two variables are strongly related in incompressible flows and both growing functions of 
$y^+$ in the log-layer (Pope Reference Pope2000), while they appear to be decoupled for highly compressible flows due to the influence of 
$M_{\infty }$ and 
${\varTheta }$ numbers. For this reason, specific flow features associated with the outer–inner scale separation, such as the enhancement of outer layer motions at high 
$Re_{\tau }$ (Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022), are not visible here, even though 
$L/ \eta$ actually increases in the outer layer.
5.2. Thermodynamic quantities
Important insights into the respective roles of Mach number and wall-cooling can also be attained by looking at root-mean-square profiles of temperature and pressure fluctuations shown in figures 11 and 12. The semilocal scaling is used to better account for fluid property variations across the boundary layer and root-mean-square quantities are scaled accordingly. In particular, root-mean-square profiles of pressure are scaled with the wall shear stress 
$\tau _w$, while the resulting scaling for temperature is obtained using the ideal gas law 
$P=R \rho T$:
\begin{equation} \frac{\tau_w}{R \bar{\rho}}=\frac{\bar{\rho} u_{\tau, S L}^2}{R \bar{\rho}}=\frac{u_{\tau, SL}^2}{R}=\gamma \bar{T} \frac{u_{\tau, S L}^2}{R \gamma \bar{T}}=\gamma \bar{T} M_{\tau, S L}^2 \end{equation}
where 
$u_{\tau, S L}= \sqrt {\tau _w/\bar {\rho }}$ is the semilocal friction velocity and 
$M_{\tau, SL}=u_{\tau, SL} / \sqrt {\gamma R \bar {T}}$ is the semilocal friction Mach number. First, the effect of 
${\varTheta }$ at a given Mach number is presented in figure 11. Considering the region starting from 
$y^*>10$, both temperature and pressure fluctuations show a reduction in intensity as 
${\varTheta }$ decreases, although more intense for the temperature. In particular, the suppression of temperature fluctuations by wall-cooling forms a plateau for the coldest case that is due to the great attenuation of the turbulent heat flux in the log-layer, consistently with Fan et al. (Reference Fan, Li and Pirozzoli2022). Around 
$y^* \approx 10$, the aforementioned attenuation of temperature fluctuations reaches its maximum for highly cooled cases (
${\varTheta }=0.25$), which is the point where mean temperature gradients are close to zero. In the near-wall region (
$y^*<10$), strongly cooled cases exhibit a peculiar behaviour, which goes in direct contrast to the monotonic attenuation of adiabatic profiles. In fact, in this region, there is an increase in the intensity of the temperature fluctuations that forms a local peak. We attribute this phenomenon to the large increase of conductive heat flux close to the wall, which is able to overcome the expected attenuation of turbulent heat flux that concurs with the generation of thermal production (see figure 13d). This is due to the large increase of near-wall temperature gradients that generate steeper mean profiles and for a wider region of 
$y^*$ values (before reaching the temperature peak), as visible in figure 3(a). The increase in pressure fluctuations in this region is shared only by the high-Mach-number case, showing that additional physical interpretations are needed on the distinct role of vorticity and acoustic modes, for which we remind to the recent study of Zhang et al. (Reference Zhang, Wan, Liu, Sun and Lu2022). Figure 12 shows the effect of Mach number at a given 
${\varTheta }$. Here, temperature fluctuation profiles are very similar up to roughly 
$y^*<15$, while the main differences are present in the outer layer, where at higher Mach numbers a peak is formed. This result indicates that 
${\varTheta }$ is an adequate parameter to recover the same general behaviour with respect to wall-cooling at different Mach numbers, as noted by Wenzel et al. (Reference Wenzel, Gibis and Kloker2022) (we remark the similarity between 
${\varTheta }$ and Eckert number). Building on the choice of the diabatic parameter 
${\varTheta }$, we can clearly see how the Mach number exerts its influence mainly beyond the buffer layer, while wall-cooling dominates the near-wall region. This is in agreement with the trends of TKE budget previously reported in § 5.1. Moreover, as discussed for velocity fluctuations § 5.1, we note the tendency of compressibility to increase the separation of scales (figure 11), while the opposite is true for enhanced wall-cooling (figure 12). This effect would be greatly reduced if profiles were compared at the same 
$Re_{\tau }^*$, which incorporates these effects (not shown). Pressure fluctuations exhibit a good collapse at the peak location around 
$y^*\approx 30$, in accordance with Zhang et al. (Reference Zhang, Wan, Liu, Sun and Lu2022), but do not share the collapse between profiles in the near-wall region noted for temperatures.
Profiles of root-mean-square temperature (a,c,e) and pressure (b,d, f) in semilocal scaling. Here, different diabatic parameters 
${\varTheta }$ are compared at a given Mach number 
$M_{\infty }$: (a,b) 
$M_{\infty }=2$; (c,d) 
$M_{\infty }=4$; (e, f) 
$M_{\infty }=6$.
Profiles of root-mean-square temperature (a,c,e,g) and pressure (b,d, f,h) in semilocal scaling. Here, different Mach number 
$M_{\infty }$ are compared at a given diabatic parameter 
${\varTheta }$: (a,b) 
${\varTheta }=0.25$; (c,d) 
${\varTheta }=0.5$; (e, f) 
${\varTheta }=0.75$; (g,h) 
${\varTheta }=1.0$.
(a–c) Production of temperature variance 
$\mathcal {P}_T$ as function of 
$y^*$ and scaled by 
$\bar {\rho } u_{\tau, SL} \tilde {T}^2 / \delta _{\nu, SL}$. Here, different wall-cooling conditions are compared for each Mach number: (a) 
$M_{\infty }=2$; (b) 
$M_{\infty }=4$; (c) 
$M_{\infty }=6$. (d) Turbulent 
$\bar {\rho } \widetilde {v^{\prime \prime } T^{\prime \prime }}$ and conductive 
$\partial \tilde {T} / \partial y$ heat transfer terms in the thermal production. Here, different wall-cooling conditions are compared for case 
$M_{\infty }=6$.
Further insights on the sources of production of temperature fluctuations, which are highly influenced by wall-cooling, can be gained by considering the temperature variance budget 
$K_T=\widetilde {T^{\prime \prime 2}}$, which can be written as (Gatski & Bonnet Reference Gatski and Bonnet2013)
$$\begin{gather} \bar{\rho} \frac{D K_T}{D t}={-}\bar{\rho} \widetilde{u_k^{\prime \prime} T^{\prime \prime}} \frac{\partial \tilde{T}}{\partial x_k}-\frac{\partial}{\partial x_k}\left(\frac{\bar{\rho} \widetilde{u_k^{\prime \prime} T^{\prime \prime2}}}{2}\right)+\gamma \overline{T^{\prime \prime}} \frac{\partial}{\partial x_k}\left(\frac{\bar{k}_T}{c_p}\frac{\partial \bar{T}}{\partial x_k}\right)\nonumber\\ +\bar{\rho} D_T-\bar{\rho} \varepsilon_T+\bar{\rho} C_T, \end{gather}$$
where the terms on the right-hand side are in order of appearance: thermal production, turbulent velocity transport, mean thermal conduction, thermal diffusion, thermal dissipation rate and contributions due to pressure–dilatation and viscous dissipation, respectively. Details on the composition of each term can be found in Gatski & Bonnet (Reference Gatski and Bonnet2013). Here, we analyse the thermal production term, which acts in a similar way to turbulent production, transferring internal energy from the mean field to the fluctuating one Fan et al. (Reference Fan, Li and Pirozzoli2022). For TBLs, its wall-normal component is the main contributor, especially as we approach the wall, which we refer to as 
$\mathcal {P}_T$:
\begin{equation} \mathcal{P}_T ={-}\bar{\rho} \widetilde{v^{\prime\prime} T^{\prime\prime}} \frac{\partial \tilde{T}}{\partial y}. \end{equation}
Here, two terms concur to the heat exchange between different flow regions by two distinct processes: the first part 
$\bar {\rho } \widetilde {v^{\prime \prime } T^{\prime \prime }}$ is dominated by turbulence with the velocity–temperature fluctuations correlation, while 
$\partial \tilde {T} / \partial y$ represent the conductive part and is related to the mean temperature gradient. Profiles of 
$\mathcal {P}_T$ are reported in figure 13, showing the effect of wall-cooling at different Mach numbers. Similarly to temperature fluctuations, cold profiles behave differently before and after 
$y^*\approx 10$, where mean temperature gradients change after the peaks. While adiabatic profiles monotonically rise from zero to a clear peak at around 
$y^*\approx 15$, proving their coupling with velocity fluctuations, cold cases progressively exhibit a reduction and outward shift of the main peak, with the creation of another peak in the viscous sublayer. An insight to understand this process, which is more apparent at 
$M_{\infty }=6$, can be gained by analysing the individual behaviour of turbulent and convective heat exchange terms in thermal production (Fan et al. Reference Fan, Li and Pirozzoli2022), which are shown in figure 13(d). Here, it can be seen that while the turbulent term is significantly far from the wall with reduced intensity for cold cases, convective heat exchange dominates the near-wall region as wall-cooling increases. In this region, even though 
$\bar {\rho }\widetilde {v^{\prime \prime } T^{\prime \prime }}$ is close to zero for all cases, the temperature gradient raises considerably for 
${\varTheta }=0.25$ which result in a non-zero product that is visible in plots of thermal production. Thus, the formation of a peak of thermal fluctuation production in the viscous sublayer is promoted. The vanishing mean temperature gradient in the buffer layer reduces the production of temperature fluctuations and promotes a decorrelation with velocity fluctuations, as discussed in the previous sections. This mechanism is clearly visible for all Mach numbers in our database, although with progressively lower intensities. This fact corroborates the importance of the role of 
${\varTheta }$ in isolating the effect of wall-cooling form compressibility, which greatly helps in the identification of individual trends.
At this point, it is possible to reconsider the qualitative results presented in figure 2 in a more quantitative way. Wall parallel slices of velocity and temperature fluctuations were taken approximately at 
$y^*=10$, where 
$\partial \tilde {T}/\partial y \approx 0$ for extremely cold cases. It is now apparent that the decorrelation between 
$\bar {\rho } T^{\prime \prime }/ \tau _w$ and 
$\sqrt {\bar {\rho }} u^{\prime \prime }/\sqrt {\tau _w}$ can be explained with the interplay of the mean temperature gradient and 
$\bar {\rho } \widetilde {v^{\prime \prime } T^{\prime \prime }}$, which entirely damps the production of temperature fluctuations. This is also visible in figure 14, which shows the joint probability density function between velocity and temperature fluctuations only for extreme cases at 
$M_{\infty }=6$ (other cases are similar). Here, a direct contrast is present between figures 14(a) and 14(b). While the latter (M6T100) shows a good correlation between the two fields, supporting their similarity, the former (M6T025) shows a strong decorrelation, especially when velocity fluctuations are negative, which is due to the influence of wall-cooling.
Scatter plot of 
$\bar {\rho } T^{\prime \prime }/ \tau _w$ vs 
$\sqrt {\bar {\rho }} u^{\prime \prime }/\sqrt {\tau _w}$: (a) 
$M_{\infty }=6,\ {\varTheta }=0.25$; (b) 
$M_{\infty }=6,\ {\varTheta }=1.0$. Here, only the cases M6T025 and M6T100 are shown. Data were collected in the same plane shown in figure 2 (
$y^*\approx 10$).
6. Conclusions
In this paper, we have presented a systematic study on the effect of the Mach number and wall-cooling on zero-pressure-gradient TBLs using DNS. A total of 12 computations have been carried out spanning 3 Mach numbers and 4 values of the diabatic parameter 
${\varTheta }$, while the friction Reynolds number has been kept constant. In this parameter space, we put emphasis on the choice of the wall-cooling parameter 
${\varTheta }$, first proposed by Zhang et al. (Reference Zhang, Bi, Hussain and She2014), which can better incorporate the indirect effects of the Mach number on wall-cooling, yielding the same integral behaviour between different cases. It is worth noting that 
${\varTheta }$ can be directly related to the Eckert number 
$Ec$, whose relevance has recently been discussed by Wenzel et al. (Reference Wenzel, Gibis and Kloker2022). These parameters show an improved ability to account for the wall-cooling effect at different Mach numbers with respect to the more classically used wall-to-recovery temperature ratio 
$T_w/T_r$, which has been shown to produce vastly different effects of wall temperature on the flow dynamics in the near-wall region across different Mach numbers (e.g. Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022; Zhang et al. Reference Zhang, Wan, Liu, Sun and Lu2022).
A summary of the most important remarks is as follows.
(i) The instantaneous flow organisation of temperature fluctuations near the wall, which for adiabatic cases is clearly discernible with the presence of near-wall elongated streaks highly correlated to streamwise velocity, breaks down as the wall temperature is progressively lowered, showing an isotropic behaviour for extremely cold cases without organised patterns. Nevertheless, a similar flow organisation is attained at different Mach numbers when
${\varTheta }$ is fixed, a first sign of the aptness of this parameter to yield the same wall-cooling effects across different $M_{\infty }$.(ii) The recent compressibility transformations of Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and Griffin et al. (Reference Griffin, Fu and Moin2021) correctly collapse all mean velocity profiles of our database to the incompressible laws of the wall. Similarly, Zhang et al. (Reference Zhang, Bi, Hussain and She2014) mean velocity–temperature relations are able to capture non-adiabatic and compressibility effects in an excellent manner. When this relation is approximated with the computed mean value of the Reynolds analogy factor
$s=0.78\pm 0.03$ (which is close to the fit of Zhang et al. Reference Zhang, Bi, Hussain and She2014), an excellent estimate is recovered, with maximum errors of $5\,\%$ from the DNS data.(iii) As the Mach number increases, we observe an increased separation between large and small scales in the outer layer measured by the ratio
$L/ \eta$, which is mainly regulated by the strong reduction of the Kolmogorov length $\eta$, and only weakly by a growth of the largest scale $L$. This effect can be effectively described by the growth of the semilocal friction Reynolds number $Re_{\tau }^*$, even though the resulting flow dynamics is different from a pure increase of the friction Reynolds number $Re_{\tau }$, the latter also leading to an increase of the inner–outer scale separation $L^+$, feeding outer layer motions. In the near-wall region, compressibility enables a less efficient redistribution of TKE, which results in a promotion of the peak of the streamwise velocity component and a decrease of the others.(iv) The enhancement of wall-cooling appears as a reduction of the large–small scale separation in the outer layer (as opposed to the effect of Mach number), which is mainly due to an increase of the Kolmogorov length scale
$\eta$ that occurs throughout the whole boundary layer thickness. Lower wall temperatures force the rise of the mean temperature peak, inducing a stratification of flow properties localised around the buffer layer. This effect is visible as an apparent promotion of compressibility, since velocity fluctuations are enhanced in the streamwise direction, while the other components are damped.(v) In the near-wall region, a dominant effect of wall-cooling is present in the root-mean-square temperature profiles and TKE budget, while the Mach number exerts its influence mainly through the buffer and log layers. When the diabatic parameter
${\varTheta }$ is kept constant, the root-mean-square temperature profiles at different $M_{\infty }$ collapse into each other near the wall, displaying a similar wall-cooling effect.(vi) For extremely cold cases (in our database
${\varTheta }=0.25$), the effect of wall-cooling is so marked that temperature fluctuations are massively damped at the point where the mean temperature gradient is zero (thus, thermal production is also zero), and a second (minor) peak arises in the viscous sublayer. This phenomenon completely decorrelates velocity and temperature fields in the near-wall region, and is more pronounced at high Mach numbers. The different behaviour of thermal production for cold cases can be explained by looking at the mean temperature gradient, which persists with a positive value for a wider wall-normal region (before reaching the mean temperature peak at $y^*_{peak}$), and with a much stronger intensity than for adiabatic cases.
Funding
We acknowledge that the results reported in this paper have been achieved using the PRACE Research Infrastructure resource MARCONI 100 based at CINECA, Casalecchio di Reno, Italy, under project PRACE call 23 n. 0058. M.B., U.B. and M.C. have also been supported by the grant agreement number A0375-2020-36614, project GREEN H2 CFD, POR FESR LAZIO 2014-2020.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available at http://newton.dma.uniroma1.it/highspeedbl/.
Appendix A
In this appendix, we consider mean velocity profiles in the framework of compressibility transformations, which aim at incorporating compressibility effects in wall-bounded flow statistics in order to recover the incompressible behaviour. Since the pioneering work of Van Driest (Reference Van Driest1951), several relations have been proposed to account for the variations of mean fluid properties, such as density and viscosity. These relations can be cast in terms of mapping functions 
$f_I$ and 
$g_I$ for wall distance 
$y_I$ and mean velocity 
$u_I$, which denote the equivalent incompressible distributions obtained from the transformation 
$I$:
\begin{equation} y_{I}=\int_0^yf_I \,{{\rm d} y} ,\quad u_I=\int_0^{\tilde{u}} g_I \,{\rm d}\tilde{u}. \end{equation}
Table 3 lists the relative values of 
$f_I$ and 
$g_I$ for Van Driest (Reference Van Driest1951) and the recent transformation of Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020), which employs a partially data-driven approach to derive the mapping exponents.
Compressibility transformations for the wall distance and the mean velocity according to (A1a,b), where 
$R=\bar {\rho }/\bar {\rho }_w$ and 
$M=\bar {\mu }/\bar {\mu }_w$.
Griffin et al. (Reference Griffin, Fu and Moin2021) transformation, instead, is based on the total stress equation, which reads
\begin{equation} \tau^+=S_t^+ \left( \frac{\tau_v^+}{S_{TL}^+}+ \frac{\tau_R^+}{S_{eq}^+}\right), \end{equation}
where 
$\tau _v^+$ and 
$\tau _R^+$ are the scaled viscous and Reynolds shear stresses (whose sum is equal to 
$\tau ^+$), while 
$S_{TL}^+=\partial U_{TL}^+ / \partial y^*$ and 
$S_{eq}^+=\partial U_{eq}^+ / \partial y^*$ are the generalised non-dimensional mean shear stresses derived for the viscous region (the subscript 
$TL$ indicated the accordance with the Trettel & Larsson (Reference Trettel and Larsson2016) velocity transformation) and for the log-layer (the subscript 
$eq$ indicates the assumption of turbulence quasi-equilibrium). The generalised non-dimensional mean shear 
$S_t^{+}=\partial U_t^+ / \partial y^*$ is the unknown and once computed it can be integrated with respect to the semilocal wall-normal coordinate 
$y^*$, leading to the transformed velocity 
$u^+_{GR} = \int S_t^+ \,{\textrm {d} y}^*$.
We report in figure 15 the scaled profiles according to the classical law of Van Driest (Reference Van Driest1951) and the latest transformations of Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and Griffin et al. (Reference Griffin, Fu and Moin2021). Figure 15(a) reveals the main weaknesses of the Van Driest (Reference Van Driest1951) scaling, whose accuracy is affected both by the increase of the Mach number and wall-cooling. In particular, non-adiabatic cases at 
$M_{\infty }=4,6$ show a clear departure from the linear law of the wall, while even adiabatic cases show a positive shift from the log-law as compressibility increases. Figures 15(b) and 15(c) show a great improvement in collapsing all profiles to the laws of the wall, the only minor discrepancy being present in the log-layer for extremely cold cases at high Mach numbers. Overall, our database supports Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and Griffin et al. (Reference Griffin, Fu and Moin2021) transformations, proving their wide range of applicability.
Mean velocity profiles at stations listed in table 2 scaled according to (a) Van Driest (Reference Van Driest1951), (b) Volpiani et al. (Reference Volpiani, Iyer, Pirozzoli and Larsson2020) and (c) Griffin et al. (Reference Griffin, Fu and Moin2021) compressibility transformations. Profiles have been translated along the 
$x$ axis according to the law 
$10^{M_\infty /2-1}$ to enable better comparison.
Appendix B
In this appendix, we report DNS results regarding theoretical relations developed to predict the coupling between velocity and temperature for both mean and fluctuating fields. First, we consider the classical relations of Walz (Reference Walz1969) and the modified relation of Zhang et al. (Reference Zhang, Bi, Hussain and She2014), as described in § 4.1. Figure 16 compares the relations (4.1) and (4.2) with the present database. As expected, Walz (Reference Walz1969) relation greatly degrades its accuracy when the wall-cooling is increased, while Zhang et al. (Reference Zhang, Bi, Hussain and She2014) is able to better perform under these conditions, the only minor deviations being present for the case M6T025. However, we note that Walz (Reference Walz1969) law still excellently holds for adiabatic cases at high Mach numbers, while being incapable of correctly capturing wall-cooling effects.
Mean temperature profiles against mean velocity compared with the classical law of Walz (Reference Walz1969) (4.1) and the modified relation of Zhang et al. (Reference Zhang, Bi, Hussain and She2014) (4.2): (a) 
$M_{\infty }=2$; (b) 
$M_{\infty }=4$; (c) 
$M_{\infty }=6$.
We then consider the SRA (Morkovin Reference Morkovin1962), which consists of a set of theoretical relations that couple the thermodynamic and kinetic fluctuating fields. Originally derived for an adiabatic case, the three main relations can be expressed as
\begin{equation} \left. \begin{gathered} \frac{(\widetilde{T^{\prime\prime2}})^{1/2}/\tilde{T}}{(\gamma-1)\tilde{M}^2(\widetilde{u^{\prime\prime2}})^{1/2}/\tilde{u}} \approx 1,\\ R_{u^{\prime\prime} T^{\prime\prime}} = \frac{\widetilde{u^{\prime\prime} T^{\prime\prime}}}{\sqrt{\widetilde{u^{\prime\prime 2}}} \sqrt{\widetilde{T^{\prime\prime 2}}}} \approx 1,\\ Pr_t=\frac{\overline{\rho u^{\prime} v^{\prime}}(\partial \tilde{T} / \partial y)}{\overline{\rho T^{\prime} v^{\prime}}(\partial \tilde{u} / \partial y)} \approx 1, \end{gathered} \right\} \end{equation}
where we recall the Favre average definition 
$\tilde {f}=\bar {\rho f}/ \bar {\rho }$ and that 
$f^{\prime \prime }=f-\tilde {f}$.
Figure 17 shows the profiles of 
$R_{u^{\prime \prime } T^{\prime \prime }}$ that clearly deviate from unity, which is expected since it was derived assuming zero total temperature fluctuation (Morkovin Reference Morkovin1962). All profiles collapse around the value 
$-R_{u^{\prime \prime } T^{\prime \prime }}=0.6$, except in the near-wall region, which is marked with an inset (Duan et al. Reference Duan, Beekman and Martín2010). The inset of figure 17 shows that the crossover location, where 
$R_{u^{\prime \prime } T^{\prime \prime }}=0$, corresponds approximately to the location of the maximum mean temperature. Here, we observe that as the wall progressively cools, the crossover location moves at higher 
$y^*$ values, indicating a temperature–velocity decorrelation that is progressively moved farther from the wall. Our database also shows that this location is almost independent of the Mach number when 
${\varTheta }$ is fixed, whereas distinct Mach and wall-cooling effects are visible on the near-wall peak intensity and position of 
$R_{u^{\prime \prime } T^{\prime \prime }}$.
Velocity and temperature correlation 
$R_{u^{\prime \prime } T^{\prime \prime }}$ as function of 
$y/\delta _{99}$. Full lines indicate 
$M_{\infty }=2$, dashed lines indicate 
$M_{\infty }=4$ and dotted lines indicate 
$M_{\infty }=6$.
The remaining two relations of (B1) have been modified over the years to account for finite heat flux at the wall and remove wall temperature dependence (Huang et al. Reference Huang, Coleman and Bradshaw1995) (HSRA). The most recent improvement has been made by Zhang et al. (Reference Zhang, Bi, Hussain and She2014), who proposed another definition of the turbulent Prandtl number 
$\overline {Pr_t}$ which should perform better at high Mach numbers, yielding the following expression of the modified SRA (modified HRSA):
\begin{equation}
\frac{(\widetilde{T^{\prime\prime2}})^{1/2}/\tilde{T}}{(\gamma-1)\tilde{M}^2(\widetilde{u^{\prime\prime2}})^{1/2}/\tilde{u}}
\overline{Pr_t} (1-(\partial\tilde{T}_t / \partial\tilde{T}))\approx 1,
\end{equation}
where the proposed definition of 
$\overline {Pr_t}$ is
\begin{equation} \overline{Pr_t}=\frac{\overline{(\rho v)' u^{\prime}}\partial \tilde{T}/\partial y}{\overline{ (\rho v^{\prime})T^{\prime}}\partial \tilde{u}/\partial y}=Pr_t \frac{1+\bar{v}\overline{\rho' u'}/\overline{\rho v' u'}}{1+\bar{v}\overline{\rho' T'}/\overline{\rho v' T'}} \end{equation}
in which the difference from the classical definition is notable when both 
$\bar {v}$ and 
$\rho '$ are non-zero.
Figure 18 compares the wall-normal profiles obtained with the original SRA and the modified version of (B2), as well as the profiles of 
$Pr_t$ and the modified 
$\overline {Pr_t}$ of (B3). Figure 18(c) shows that the modified version of Zhang et al. (Reference Zhang, Bi, Hussain and She2014) clearly improves the insensitivity to the freestream Mach number and wall temperature condition, with only slight deviations at the edge of the boundary layer. It is also interesting to note the excellent collapse that the original SRA of panel 18(a) exhibits for profiles at fixed 
${\varTheta }$, independently of the Mach number, highlighting the relevance of the diabatic parameter 
${\varTheta }$ in accounting for the effects of different wall temperatures independently of the Mach number.
Comparison of the SRAs in the (a,b) original form (B1), (c,d) modified HSRA (Zhang et al. Reference Zhang, Bi, Hussain and She2014) ((B2), (B3)). Full lines indicate 
$M_{\infty }=2$, dashed lines indicate 
$M_{\infty }=4$ and dotted lines indicate 
$M_{\infty }=6$.
Open access
