3.1. Characterization of the baseline flow in the empty tunnel without ice

Before presenting the results obtained in the experiments over melting ice, the baseline flow (i.e. in the absence of ice) was characterized. The principal quantities of interest were the mean longitudinal velocity ($\langle U \rangle$), RMS longitudinal velocity ($u_{rms}$) and the Reynolds stress ($\langle uv \rangle$). These quantities are presented non-dimensionally in figure 2 for free-stream velocities ($U_{\infty }$) of 1.1, 2.0 and 3.2 m s$^{-1}$, where they are scaled by the free-stream velocity, mean velocity at the given wall-normal position and product of $u_{rms}$ and $v_{rms}$, respectively. The latter two non-dimensional quantities form the turbulent intensity ($Ti$) and correlation coefficient ($\rho _{uv}$), respectively.

Figure 2. Non-dimensional profiles in the turbulent boundary layer. (a) Non-dimensional velocity, (b) turbulent intensity and (c) correlation coefficient.

The non-dimensional mean velocity profiles of figure 2(a) are similar for all three free-stream velocities investigated, but do show a Reynolds-number dependence as expected. They show good agreement with prior results for turbulent boundary layers (e.g. Ohya et al. Reference Ohya, Neff and Meroney1997). The shape factor ($H$) was calculated using the displacement thickness ($\delta ^{*}$) and the momentum thickness ($\theta ^{\delta }$) as follows (Pope Reference Pope2000):

(3.1)\begin{equation} H = \frac{\delta^{*}}{\theta^{\delta}}. \end{equation}

The shape factors were found to be 1.33, 1.25 and 1.21 for the 1.1, 2.0 and 3.2 m s$^{-1}$ free-stream velocities, respectively. From $y/\delta = 0.1$ to 0.5, the turbulent intensity profiles (figure 2b) are similar, but diverge outside of this range, where the turbulent intensity is greater for higher free-stream velocities when $y/\delta > 0.5$. Conversely, the turbulent intensity decreases for higher free-stream velocities when $y/\delta < 0.1$. The observations in figure 2(c) are consistent with the expected value of the correlation coefficient $\rho _{uv}$ (approximately $-$0.4) in the core region of a turbulent boundary layer (Tennekes & Lumley Reference Tennekes and Lumley1972).

The mean velocity profiles were then used to estimate the friction velocity ($u_{*}$) using two different methods – one using the velocity defect law (Tennekes & Lumley Reference Tennekes and Lumley1972) and the other using Prandtl's correlation for the skin friction in a turbulent boundary layer (Prandtl Reference Prandtl1925). Details of the methods for estimating the friction velocity can be found in Harrison (Reference Harrison2022).

The estimated friction velocity is presented in table 1, along with the inferred roughness Reynolds number ($R_{k}$) for the ice for a given free-stream velocity using the surface roughness measurements described in § 2. In all cases, $R_{k}$ is small enough (i.e. $<5$; see Tennekes & Lumley Reference Tennekes and Lumley1972) to justify the characterization of the glacier model as an aerodynamically smooth surface.

Table 1. Estimated friction velocities and corresponding values of $R_{k}$ ($k_{rms}=0.4$ mm).

To identify the non-dimensional wall-normal distances associated with the various measurement locations and free-stream velocities in this work, the non-dimensional wall-normal heights quantified using both wall units ($y^+ = u_{*}y / \nu$) and scaled by the boundary layer thickness ($y/\delta$) are tabulated in table 2.

Table 2. Non-dimensional wall-normal heights for each run shown in outer scaling ($y/\delta$) and wall units ($y^+$). The measurements at the lower two wall-normal positions ($y/\delta = 0.13$ and $0.33$) lie in the log-law region of the boundary layer (i.e. $y^+ \gtrsim 30, y/\delta \lesssim 0.3$; Pope Reference Pope2000). The measurements at the highest wall-normal position ($y/\delta = 0.53$) lie in the outer layer.

In addition to the statistics presented in figure 2, further characterization of the baseline flow was performed at the three heights for which the measurements over ice were made. This analysis consisted of calculating other relevant quantities (e.g. $\epsilon$, $\lambda$ and $\eta$), as well as spectra and PDFs of the velocity fluctuations. The dissipation rate ($\epsilon$) was calculated using Taylor's frozen-flow hypothesis (Taylor Reference Taylor1938) and assuming local isotropy such that

(3.2)\begin{equation} \epsilon = \frac{15\nu}{\langle U \rangle^{2}} \left\langle \left(\frac{\partial u}{\partial t}\right)^{2} \right\rangle. \end{equation}

The Taylor microscale ($\lambda$) was calculated using

(3.3)\begin{equation} \lambda^{2} \equiv \langle u^{2} \rangle \left\langle \left(\frac{\partial u}{\partial x}\right)^{2} \right\rangle^{-1} \end{equation}

to facilitate the calculation of the Taylor microscale Reynolds number, $R_{\lambda }$, defined as

(3.4)\begin{equation} R_{\lambda} \equiv \frac{u_{rms}\lambda}{\nu}. \end{equation}

The Kolmogorov microscale ($\eta$) is defined in the classical way:

(3.5)\begin{equation} \eta \equiv \left(\frac{\nu^{3}}{\epsilon}\right)^{1/4}. \end{equation}

The one-dimensional wavenumber spectra of the longitudinal and transverse velocity fluctuations ($E_{u}(\kappa _{1})$ and $E_{v}(\kappa _{1})$, respectively) are presented non-dimensionally in figure 3 for the 10 and 40 mm heights. The width of the inertial range of the spectra was found to increase with Reynolds number, and tended towards a $-5/3$ slope. The slope of the $u$ velocity spectra was typically steeper than that of the $v$ velocity, consistent with prior observations (e.g. Mydlarski & Warhaft Reference Mydlarski and Warhaft1996). Overall, the non-dimensional spectra show little sensitivity to the measurement height. The non-dimensional longitudinal velocity spectra collapse well for the two higher-speed runs, while the lowest-speed runs were smaller in magnitude. This difference is likely attributable to different Reynolds numbers of the three flows, given that $R_{\lambda }>80$ for the higher-speed cases and $R_{\lambda }<40$ for the lower-speed one.

Figure 3. Non-dimensional wavenumber spectra of longitudinal (a,c) and transverse (b,d) velocity fluctuations in the empty tunnel, where $\kappa _{1}=2{\rm \pi} f/\langle U \rangle$. The straight line is a $-$5/3 power law. (a) $E_{u}$, $y = 10$ mm, $y/\delta = 0.13$, (b) $E_{v}$, $y = 10$ mm, $y/\delta = 0.13$, (c) $E_{u}$, $y = 40$ mm, $y/\delta = 0.53$ and (d) $E_{v}$, $y = 40$ mm, $y/\delta = 0.53$.

The PDFs of the longitudinal ($u$) and vertical ($v$) velocity fluctuations were also computed, then non-dimensionalized using their RMS values. They were found to be very nearly Gaussian, with slight deviations from Gaussian behaviour in the tails. These deviations were more apparent in the vertical than the longitudinal velocity fluctuations. Additionally, PDFs of the instantaneous Reynolds stress ($uv$) were computed and non-dimensionalized similar to those of the single-component velocity PDFs. While general collapse was observed for all three speeds, it was found that the tails of these PDFs decreased more rapidly at lower speeds. There is an exponential decay in the tails of the PDFs of $uv$, an expected consequence of the nearly joint-Gaussian distribution (e.g. Mydlarski Reference Mydlarski2003). Characteristic PDFs of the longitudinal and transverse velocity fluctuations, as well as the instantaneous Reynolds stress, for the 25 mm height are presented in figure 4.

Figure 4. Non-dimensional PDFs of the longitudinal velocity fluctuations (a), transverse (vertical) velocity fluctuations (b) and instantaneous Reynolds stress (c) for the 25 mm height ($y/\delta = 0.33$) in the empty tunnel.

In addition to the spectra and PDFs of the two velocity components presented above, co-spectra and joint PDFs of the velocity fluctuations were also computed. The magnitude of the co-spectrum, $\varLambda _{uv}(\kappa _{1})$, was scaled by the corresponding velocity spectra to form the non-dimensional coherence spectrum:

(3.6)\begin{equation} H_{uv}(\kappa_{1}) \equiv \frac{\varLambda_{uv}(\kappa_{1})}{E_{u}(\kappa_{1})E_{v}(\kappa_{1})}. \end{equation}

Closer to the wall, there is higher coherence at larger Reynolds numbers (i.e. free-stream velocities) due to the stronger correlations of $u$ and $v$, as demonstrated in figure 5. Similarly, joint PDFs of $u$ and $v$ were also computed, and were found to be nearly joint-Gaussian. They did not show significant sensitivity to the height or free-stream velocity, such that only three of the nine cases in the test envelope are presented in figure 6 as characteristic representations of the joint PDFs of $u$ and $v$ in the empty tunnel without ice.

Figure 5. Non-dimensional coherence spectra of the longitudinal and transverse velocity fluctuations for the empty tunnel: (a) $y/\delta = 0.13$, (b) $y/\delta = 0.33$ and (c) $y/\delta = 0.53$.

Figure 6. Characteristic contour plots of joint PDFs of the longitudinal and transverse velocity fluctuations for the empty tunnel: (a) $y/\delta = 0.13$, $U_{\infty } = 1.1$ m s$^{-1}$, (b) $y/\delta = 0.33$, $U_{\infty } = 2.0$ m s$^{-1}$ and (c) $y/\delta = 0.53$, $U_{\infty } = 3.2$ m s$^{-1}$.

3.2. Stationary velocity field over melting ice

Having characterized the baseline flow in the absence of ice for the nine cases in the test envelope, the results could be compared with those of the experiments over melting ice to understand the effect of the ice surface on the velocity statistics. However, given that the warming and melting of an ice patch is a (thermally) unsteady problem, it is reasonable to expect that the statistics of both velocity and temperature measured in the turbulent boundary layer over melting ice would undergo a transient evolution. The transient evolution of the melting ice is discussed in detail in the next subsection; however, it was found that the velocity field statistics did not change significantly in time throughout the experiments. This is illustrated by plotting the local time averages of the mean longitudinal and the RMS vertical velocities as a function of time in figure 7 for the lowest measurement height (which will be the most influenced by changes in temperature). The local averages are non-dimensionalized by their final value (i.e. those of the last averaging period) to quantify the relative changes with time. It can be seen that the quantities do not change appreciably throughout the course of the experiments and, as such, the velocity field measured over melting ice for the experiments in the present work can be reasonably considered to be statistically stationary.

Figure 7. Evolution of local time averages of velocity quantities $y/\delta = 0.13$. Local averages were calculated in subsets of $8\times 10^{5}$ samples, lasting approximately 2.3 min each. (a) Mean longitudinal velocity and (b) RMS vertical velocity.

Having demonstrated the statistical stationarity of the velocity field over melting ice, the preceding analysis for the baseline cases can be repeated for the experiments over ice. Velocity statistics (averaged over the entire experimental run) are presented in table 3. Since it of interest to investigate the effects of the stratification on the velocity field over melting ice, quantification of the stratification (such as a Richardson number) is needed to make meaningful observations. While it has been demonstrated that the velocity statistics are statistically stationary, the same cannot be said for the temperature and combined velocity–temperature statistics. This is discussed in further detail in the following subsection. However for the purposes of analysing the stationary velocity statistics, the turbulent Richardson number (as defined by Yoon & Warhaft (Reference Yoon and Warhaft1990)) was found to vary only slightly in time after an initial decrease during rapid warming of the ice. For this reason, it has been selected as a parameter to quantify the relative effects of varying levels of stratification between the experiments performed in the present work. The bulk Richardson number was also found to be nearly stationary in time, but did not vary as significantly between runs as the turbulent Richardson number, as can be observed in table 3.

Table 3. Velocity statistics measured over melting ice. The integral length scale ($\ell$) was found by integrating the autocorrelation function of the longitudinal velocity.

The velocity statistics measured over melting ice and averaged over the entire time series of each experiment were then compared with those of the baseline cases. The quantities of particular interest are the second-order statistics (i.e. the variances and co-variance), the turbulent kinetic energy in the $x$–$y$ plane:

(3.7)\begin{equation} {\rm TKE}_{xy} \equiv \langle u^{2} \rangle + \langle v^{2} \rangle \end{equation}

and the anisotropy ratio:

(3.8)\begin{equation} \phi \equiv \frac{v_{rms}}{u_{rms}}. \end{equation}

To this end, the percentage changes in these values between the experiments over ice and the baseline cases were evaluated. The average percentage change for each of these quantities, as well as their correlation with the turbulent Richardson number, is presented in table 4. In this case, the correlation is measured using the Pearson correlation coefficient, calculated as

(3.9)\begin{equation} \rho_{(\,p, Ri_{t})} \equiv \frac{\displaystyle \sum_{n=1}^{N}[(\,p_{n}-\bar{p}) (Ri_{t_{n}}-\overline{Ri_{t}})]}{\sqrt{\displaystyle \sum_{n=1}^{N} (\,p_{n}-\bar{p})^{2}\sum_{n=1}^{N} (Ri_{t_{n}}-\overline{Ri_{t}})^{2}}}, \end{equation}

where $Ri_{t_{n}}$ is the average Richardson number for experiment $n$ in table 3, $p$ is the percentage change in a quantity between the experiments over ice versus the baseline cases, where an overbar represents the average across all $N = 9$ runs.

Table 4. Average percentage change relative to baseline cases and correlation with $Ri_{t}$ for velocity statistics of interest.

From the data in table 4, one observes a measurable increase in the magnitude of the streamwise velocity variance in the experiments over ice. This result is somewhat unexpected, since profiles of velocity statistics measured in stably stratified boundary layers over solid, non-melting surfaces show a decrease in all velocity components with increasing stratification (e.g. Ohya et al. Reference Ohya, Neff and Meroney1997). Furthermore, the observed increase in streamwise velocity variance was not found to be strongly correlated with the turbulent Richardson number, or any of the other quantities calculated in table 3. At the same time, there are clear decreases in both the vertical velocity variance and Reynolds stress in the experiments over ice when compared with the baseline cases, which agrees well with prior research relating to stably stratified boundary layer flows over solid surfaces (e.g. Arya Reference Arya1975; Ohya et al. Reference Ohya, Neff and Meroney1997; Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017), as buoyancy forces extract turbulent kinetic energy from the vertical fluctuations, thereby reducing the Reynolds stress (Arya Reference Arya1975). This is observed in the present work, wherein the suppression of vertical velocity variance is clearly negatively correlated with the turbulent Richardson number, and the decrease in the Reynolds stress is also negatively correlated with $Ri_{t}$.

However, this partial agreement in the second-order statistics with previous experiments remains intriguing, given that Arya (Reference Arya1975) and Williams et al. (Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017) both state that the interaction between the Reynolds stress and mean velocity gradient should also decrease the streamwise velocity variance. These results suggest that there may be effects (potentially related to the boundary between the air and melting ice, or perhaps related to latent fluxes of heat) which differentiate this flow from an otherwise identical stably stratified boundary layer flow over a non-melting surface. It was found that the difference in the turbulent kinetic energy in the $x$–$y$ plane between the baseline cases and experiments over ice was negligible. The magnitude of the decrease in vertical velocity variance was offset by the increase in streamwise velocity variance such that the total turbulent kinetic energy measured in the $x$–$y$ plane over melting ice was nearly unchanged relative to the baseline cases.

The anisotropy ratio has been previously demonstrated to be relatively insensitive to the effects of stratification below the critical Richardson number (Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017). The documented insensitivity of the anisotropy ratio in the literature is attributed to decreases in the streamwise velocity variance by the reduced values of the Reynolds stresses that are proportional to those of the vertical velocity variance due to buoyancy, such that $\phi$ is unchanged (Arya Reference Arya1975). However, in the present work, the anisotropy ratio was found to decrease in the experiments over ice compared with the baseline cases. This change is negatively correlated to the turbulent Richardson number, as shown in table 4. In other words, the anisotropy ratio does exhibit some sensitivity to changes in stratification, contrary to observations in stably stratified flows observed over solid surfaces. Conversely, studies of the downstream evolution of grid-generated stratified turbulence have found a relationship between the anisotropy and Richardson number (Yoon & Warhaft Reference Yoon and Warhaft1990; Piccirillo & van Atta Reference Piccirillo and van Atta1997). However, such an effect has not been observed in the boundary layer over solid (i.e. non-melting) surfaces, and may therefore be related to the change in the nature of the boundary condition at the surface when it melts.

Since the stratification as quantified by the bulk Richardson number (given in table 3) is relatively weak when compared with the critical value of 0.25 (Ohya et al. Reference Ohya, Neff and Meroney1997), it is possible that the observed changes in second-order statistics measured over melting ice are not directly caused by the stratification, but rather are physical phenomena related to the phase change at the ice's surface. This physical phenomenon could be an aerodynamic effect related to the Reynolds number, a stratification effect related to the Richardson number, an effect caused by the latent fluxes of heat altering the flow, or a combination of these. In any case, it may be concluded that the nature of the turbulent velocity field over melting ice is different from that of a nearly identical field over a non-melting surface, such that prior observations and empirical relations derived from the latter may not be easily applicable to the former. This could have significant ramifications for the understanding and modelling of turbulent flows over melting glaciers, for which there is already a demonstrated need for improvement (Hock Reference Hock2005).

In an effort to explain the observed differences between the baseline cases and experiments measured over ice, the following hypothesis is offered. As the surface changes from solid ice to liquid water, a thin layer of water forms atop the ice. The interface between air and water (two fluids) compared with that between air and ice may cause additional streamwise velocity fluctuations in the turbulent boundary layer, thereby explaining the increased streamwise velocity variance which was not clearly related to Richardson number. The observed reduction in Reynolds stress could be a result of the lower velocity gradient arising from a non-zero velocity at the surface, which has been observed in studies of liquid-infused surfaces for the purpose of drag reduction (e.g. Fu et al. Reference Fu, Arenas, Leonardi and Hultmark2017). The decrease in vertical velocity variance may then be attributed either to the effects of stratification through the classical mechanism (Tennekes & Lumley Reference Tennekes and Lumley1972) or to the reduced wall shear stress.

The complex nature of this problem makes the evaluation of this hypothesis challenging, as well as the comparison to prior work, such that results found in adjacent areas of study may lend insight. In studies of large-scale atmospheric flows over sea ice, the buoyant motion of air can be attributed to differences in the surface heat conduction between solid ice and liquid water (Brümmer et al. Reference Brümmer, Busack, Hoeber and Kruspe1994), and that greater surface shear stress can be expected over ice than over water (Overland Reference Overland1985; Kantha & Mellor Reference Kantha and Mellor1989).

It could be hypothesized that these effects may be attributed to the decoupling of the flow due to the step change in surface temperature (Mott et al. Reference Mott, Schlögl, Dirks and Lehning2017). Boundary layer decoupling can be characterized by the suppression of turbulent mixing close to the surface, where shear effects are typically strongest (Mott et al. Reference Mott, Paterna, Horender, Crivelli and Lehning2016). However, the effects of boundary layer decoupling tend to amplify some of the expected effects of stratification, and do not adequately explain the clear increase in the longitudinal velocity variance which contradicts the decrease in the vertical velocity variance and momentum flux. As such, while the effects of boundary layer decoupling are important in the study of melting glaciers and are discussed later in this paper, they do not offer an explanation as to why the streamwise velocity variance measured over melting ice increased relative to the baseline cases when the vertical velocity variance and momentum flux both decreased, as expected.

In studies of small-scale turbulence over liquid films, surface waves will greatly influence the structure of the turbulence above a wavy surface (Kendall Reference Kendall1970), with impacts as high as 40 % (Sullivan et al. Reference Sullivan, McWilliams and Moeng2000). These effects have been observed for liquid films as thin as 1.8 mm (Cohen & Hanratty Reference Cohen and Hanratty1965). In the presence of stable stratification, the ‘effective’ surface wave speed is increased, further decreasing the skin friction coefficient (Sullivan & McWilliams Reference Sullivan and McWilliams2002) and surface shear stress. In both cases, there are explanations which partially support the present observations to some extent.

It is clear that the addition of the phase-changing ice surface adds further complexity to the problem, as there now may be effects that are attributable to neither the Reynolds nor Richardson numbers. Studies have already attempted to separate the effects of these two parameters (e.g. Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017); however, the addition of the phase-changing surface in the present work makes it especially difficult to attribute the changes between the velocity field over melting ice and that of a neutral boundary layer to either Reynolds number, Richardson number or surface phase effects. In the context of modelling the turbulent flow over melting glaciers, the present findings are significant, as the small-scale nature of the surface condition can have important consequences for the transport of momentum and heat in the outer region of the boundary layer over melting ice. Further work should be conducted to attempt to distinguish the change in the velocity statistics of a similarly stratified turbulent flow over both a dry (non-melting) surface and a melting one, as presently discussed.

Nevertheless, it is of interest to compare the (statistically stationary) results obtained over melting ice with those of the baseline cases using the PDFs and spectra. For all of the experiments herein, it was found that the non-dimensional wavenumber spectra measured over melting ice did not exhibit a significant change when compared with the baseline cases, i.e. changes were less than the experimental uncertainty. To illustrate this observation, a comparison of the longitudinal and transverse velocity spectra for the $U_{\infty }= 2.0$ and $U_{\infty }= 1.1$ cases, respectively, measured at the lowest height in the test envelope (where the effects of the ice should be strongest) is presented in figure 8. The (non-dimensional) spectra show little-to-no difference between the baseline cases and experiments over ice, implying that the turbulent cascade in these experiments is unchanged by the presence of the melting ice. To better understand the effects of stratification on the velocity spectra, they are plotted alongside the traditional Kolmogorov–Obukhov–Corrsin (KOC) $\kappa ^{-5/3}$ scaling (Tennekes & Lumley Reference Tennekes and Lumley1972) for passive scalar turbulence, as well as the $\kappa ^{-11/5}$ scaling for stably stratified turbulent flows (Bolgiano Reference Bolgiano1959; Obukhov Reference Obukhov1959). The spectra can exhibit both of these scalings in different regions of the turbulent cascade: $\kappa ^{-11/5}$ scaling for $\kappa < \kappa _{B}$, where $\kappa _{B}$ is the Bolgiano wavenumber, and $\kappa ^{-5/3}$ scaling when $\kappa > \kappa _{B}$ (Verma, Kumar & Pandey Reference Verma, Kumar and Pandey2017). This is because the effects of buoyancy are assumed to be negligible at smaller scales (Alam, Guha & Verma Reference Alam, Guha and Verma2019). Although it is difficult to compare the agreement of these two predicted scalings with accuracy, it seems that spectra of figure 8 are more consistent with KOC scaling. Given the strong spectral similarity between the baseline cases and the experiments performed over ice, it may therefore be concluded that the stratification in the present work is too weak to exert a significant influence on the velocity spectra. This is not unreasonable given the Reynolds numbers and Richardson numbers of the experiments of the present work, as they fall into a region of weakly stably stratified turbulence (Alam, Verma & Joshi Reference Alam, Verma and Joshi2023), wherein the temperature field over the ice may be effectively considered a passive scalar field. Moreover, given (i) the observed changes in the previously discussed single-point statistics, but (ii) the similarity in the power spectra, it may be appropriate to attribute the decrease in $\langle v^{2} \rangle$ and $\langle uv \rangle$ (at least partially) to the decoupling of the flow caused by the step change in surface temperature, since it does not appear that the effects of stratification are observable in the spectra measured in the present work.

Figure 8. Comparison of the baseline non-dimensional wavenumber spectra with those measured over melting ice for the (a) longitudinal ($y/\delta = 0.13$, $U_{\infty } = 2.0$ m s$^{-1}$) and (b) transverse ($y/\delta = 0.13$, $U_{\infty } = 1.1$ m s$^{-1}$) velocity fluctuations.

Having observed spectral similarity between the baseline cases and experiments over ice, it is also useful to evaluate how the presence of the ice changes the PDFs of $u$, $v$ and $uv$. Here, the single-component PDFs are compared with a Gaussian distribution, given the previously observed Gaussian behaviour in the baseline cases. Such behaviour was observed once again for the experiments over melting ice for all of the nine test cases. However, the greatest deviation from Gaussian behaviour was observed for the runs nearest the ice surface ($y/\delta =0.13$), where the temperature gradient is strongest. Figure 9 depicts the normalized PDFs for the longitudinal and transverse velocity fluctuations. Compared with the baseline cases, the deviation from Gaussian behaviour occurs in the tails (i.e. $|{v}/{v_{rms}}| > 2$) of the transverse velocity PDFs and is generally stronger than that of the longitudinal velocity PDFs. Furthermore, while the baseline PDFs of the longitudinal velocity at $y/\delta = 0.13$ showed no significant deviation from a Gaussian PDF, there is an observable difference at the tails in figure 9(a). As the height is increased, this deviation disappears, indicating a decreased effect of the temperature field on the longitudinal velocity fluctuations. This effect of the ice on the PDFs of the transverse velocity fluctuations does not decrease with height as rapidly as for the longitudinal velocity fluctuations. This is illustrated by comparing the PDFs for the longitudinal and transverse velocity fluctuations at $y/\delta = 0.33$ (figures 9d and 9e). At this height, the tails of the longitudinal velocity PDFs retain their Gaussian behaviour as observed in the baseline cases. In contrast, the tails of the transverse velocity PDFs at this height are super-Gaussian. At $y/\delta = 0.53$, PDFs of both longitudinal and transverse velocity appear nearly Gaussian, indicating negligible influence of the ice surface on the velocity PDFs at this height.

Figure 9. Non-dimensional PDFs of the longitudinal velocity fluctuations (a,d,g), transverse velocity fluctuations (b,e,h) and instantaneous Reynolds stress (c, f,i) measured over melting ice: (a–c) $y/\delta = 0.13$, (d–f) $y/\delta = 0.33$ and (g–i) $y/\delta = 0.53$.

Lastly, the PDFs of $uv$ over melting ice are considered. They have similar shapes to the baseline cases shown in figure 4; however, the behaviour at the tails is opposite to that of the baseline cases – the PDFs of $uv$ over melting ice show better collapse at all speeds farther from the wall, and show the most sensitivity to free-stream velocity nearest to the wall. The differences with free-stream velocity appear to be restricted to the positive tail of the PDF, whereas the negative tail shows good collapse for all nine cases. Given that that deviations only exist in the region where $(uv)>5(uv)_{rms}$, the significance of these differences between the baseline cases and experiments over ice will be small.

Having compared the spectra and PDFs for both velocity components with the respective baseline cases, combined velocity statistics are now investigated. Since no significant difference between the experiments over melting ice and the baseline cases was observed in the non-dimensional wavenumber spectra for either of the two velocity components, a similar result was expected for the coherence spectra. However, this was not found to be the case. Nearest to the surface (at $y/\delta = 0.13$), the coherence spectra exhibit a decrease relative to the baseline cases that grows with Reynolds number, as shown in figure 10. This observation is presumably a result of the effects of boundary layer decoupling observed by Mott et al. (Reference Mott, Paterna, Horender, Crivelli and Lehning2016), as previously discussed.

Figure 10. Comparison of the baseline coherence spectra with those measured over melting ice at $y/\delta = 0.13$: (a) $U_{\infty }=1.1$ m s$^{-1}$, (b) $U_{\infty }=2.0$ m s$^{-1}$ and (c) $U_{\infty }=3.2$ m s$^{-1}$.

For the lowest free-stream velocity (i.e. when the Richardson number is largest), the coherence over melting ice shows the most similarity to the baseline case. The difference increases as the Richardson number decreases, such that it does not seem likely that the decrease in coherence is caused by the effects of stratification. For the other two measurement heights, there was no observed difference between the coherence of the experiments over ice and the baseline cases for any of the free-stream velocities investigated. The effects of the ice on the coherence appear to be limited to higher Reynolds numbers nearest to the surface, when the shear is highest. These results agree well with the literature surrounding boundary layer decoupling, wherein turbulent energy is most strongly suppressed closest to the ground (Derbyshire Reference Derbyshire1999).

It is then worth comparing the joint PDFs of $u$ and $v$ for two cases with the same free-stream velocity, but measured at different heights, to understand how the stronger reduction in Reynolds stress (i.e. coherence) observed near the wall vanishes as the height is increased. To this end, the joint PDFs are plotted for two heights at $U_{\infty }=3.2$ m s$^{-1}$, beside the baseline joint PDFs for comparison, in figure 11. Similar to the PDFs of $u$ and $v$, the joint PDFs appear largely similar with only marginal differences near the tails.

Figure 11. Comparison of the joint PDFs for the baseline cases (a,c) with those measured over ice (b,d) for $U_{\infty }=$3.2 m s$^{-1}$: (a,b) $y/\delta = 0.13$ and (c,d) $y/\delta = 0.33$.

3.3. Transient evolution of temperature-related quantities over melting ice surface

In this subsection, the transient evolutions of the melting ice surface and the turbulent flow over it are investigated to address one of the principal objectives of this work. To this end, the evolution of the glacier model, and the temperature-related statistics measured in the turbulent boundary layer above it, are discussed below.

During the experiments, the ice surface transitioned through three distinct ‘phases’ of the melt process. At the beginning of each experiment, the subsurface thermocouple reading was typically around $-10\,^\circ$C. When the wind tunnel was turned on, the ice surface would begin to warm due to the convective heat transfer from the room temperature air blowing over it, and a thin layer of frost would then form on the ice surface. The subsurface temperature would rapidly increase during this time. This phase is referred to as the warming phase, wherein the surface of the glacier model in contact with the airflow was solid ice.

Once the ice surface reached 0 $^\circ$C, it would begin to melt. This process would also begin to remove any minor surface topographical or roughness features. During this time, it can be assumed that the surface temperature of the ice was constant at 0 $^\circ$C, given that the temperature of a phase-changing substance is constant if its pressure does not change (Bergman et al. Reference Bergman, Bergman, Incropera, DeWitt and Lavine2011). The subsurface temperature would remain constant during this time, slightly below 0 $^\circ$C since it was embedded just below the ice's surface. This phase is referred to as the melting phase, wherein the surface layer of the glacier model in contact with the airflow was a thin layer of water sitting atop solid ice.

As the ice continued to melt, the water layer on top of the ice would grow deeper, removing all surface topographical features as the water layer thickness exceeded the surface roughness. Eventually, the water layer would become thick enough to have an internal temperature gradient of its own. This is referred to as the pooling phase. During the pooling phase, the surface was completely smooth due to the surface tension of the water, but the temperature at the air–water interface was unknown. Measurement of the true surface temperature of the liquid was not possible in this work, given the challenges associated with assessing the temperature jump at the air–water interface (Gatapova et al. Reference Gatapova, Graur, Kabov, Aniskin, Filipenko, Sharipov and Tadrist2017). Attempts were made to measure the surface temperature by use of infrared thermometry. However, the measurements were deemed to be of insufficient quality, given the complications arising from changes in the emissivity of the ice with time as its surface increased in temperature, followed by the even more significant changes in the emissivity as the surface melted and changed phase. Compounding these effects with directional dependencies and obstructed views of the surface arising from the hot-wire probes in the wind tunnel, accurate radiative estimates of the surface temperature were found to be unfeasible. For these reasons, we determined the most consistent and least error-prone approach to be measurement using a subsurface thermocouple, as employed herein.

As a consequence of this decision, the exact condition necessary for transition between the melting and pooling phases is somewhat arbitrary without precise knowledge of the surface temperature, but may be inferred from an increase in the subsurface temperature measured by the embedded thermocouple. By the end of the experiments, the thickness of the water layer was never more than 10 mm. Depending on a glacier's environmental conditions, the existence of a water layer on its surface could represent a key difference between the boundary condition that exists at the surface of an actual glacier and that of the present fundamental investigation. The influence of the water layer thickness on the turbulent quantities measured over its surface was discussed in detail in Harrison (Reference Harrison2022) and remains to be fully understood. But given the previously demonstrated stationarity of the velocity statistics over the course of the experiments (during which time the water layer thickness was continually increasing), it may be reasonably assumed that the thickness of the water layer above the ice did not significantly influence the results in the present work.

The point in time at which the ice would transition between these phases was dependent on the free-stream velocity, which controlled the melt rate of the glacier model. An appropriate time scale for the bulk melt rate ($\tau _{m}$) can be determined using an approximate energy balance between the melt heat flux ($Q_{M}$) and convective heat flux ($Q_{H}$). The melt heat flux can be expressed as

(3.10)\begin{equation} Q_{M} = \frac{m L_{f}}{\tau_{m} A} = \frac{\rho d L_{f}}{\tau_{m}},\end{equation}

where $m$ is the mass of the melting ice, $L_{f}$ is the specific latent heat of fusion and $A$ is the surface area in contact with the turbulent airflow. For the rectangular-prism-shaped block of ice used herein, the mass can be re-expressed as $\rho _{ice}\times A\times d$, where $\rho _{ice}$ is the density of the ice and $d$ is the depth of ice. The bulk convective heat flux from the air to the ice surface is given by

(3.11)\begin{equation} Q_{H} = h(T_{s}-T_{\infty}),\end{equation}

where $h$ is the convective heat transfer coefficient. For convective heat transfer to an isothermal flat plate (a reasonable approximation of the glacier model), the characteristic length scale in the local Nusselt number (defined here as $Nu_x \equiv hx/k_{f}$, where $k_{f}$ is the thermal conductivity of the fluid) is given by the streamwise position ($x$). In this series of wind-tunnel experiments, there is an unheated starting length ($\zeta$) over which the velocity boundary layer grows before reaching the ice surface. In this case, the semi-empirical correlation from the literature (Bergman et al. Reference Bergman, Bergman, Incropera, DeWitt and Lavine2011) enables the Nusselt number to be estimated as

(3.12)\begin{equation} Nu_{\zeta} = \frac{0.0296 (Re_{x}^{0.8}) Pr^{0.33}}{\left[1-\left(\dfrac{\zeta}{x}\right)^{0.9}\right]^{1/9}},\end{equation}

where $x$ is the total downstream distance from the grid to the measurement probe. Combining (3.10)–(3.12), the melt time scale can be found to be

(3.13)\begin{equation} \tau_{m} = \frac{\rho_{ice} L_{f} d}{T_{s}-T_{\infty}}\left(\frac{x}{k_{f}}\right) \left(\frac{\left[1-\left(\dfrac{\zeta}{x}\right)^{0.9}\right]^{1/9}}{0.0296 (Re_{x}^{0.8}) Pr^{0.33}}\right) .\end{equation}

The exact value of the characteristic depth ($d$) is somewhat arbitrary, since most of the terms in (3.13), apart from $Re_{x}$, are constant for all of the experiments performed herein. In this case, a characteristic ice depth of 1 mm was selected. Given this choice, the melt time scales ($\tau _{m}$) were found to be approximately $6.1\times 10^{4}$, $3.8\times 10^{4}$ and $2.6\times 10^{4}$ s for free-stream velocities of 1.1, 2.0 and 3.2 m s$^{-1}$, respectively. Using this melt time scale, the evolution of the subsurface temperature from the embedded thermocouple is plotted as a function of the non-dimensional time ($t/\tau _{m}$) in figure 12. The subsurface temperature is herein referred to as $T_{ss}$, serving as the closest possible estimate of the true surface temperature in this work. Figure 12 shows that $\tau _{m}$ is an appropriate time scale, since the plots of the subsurface temperature collapse well when plotted as a function of non-dimensional time. The warming phase is observed in the region where $T_{ss}$ is changing the most rapidly. The melting phase begins at approximately $0.06\tau _{m}$. The increase in temperature at approximately $0.2\tau _{m}$ corresponds the point in time at which the pool of water atop the ice surface reaches the depth of the thermocouple (approximately 5 mm). Since the surface temperature is a key parameter affecting the statistics of temperature and its related quantities in a thermal boundary layer, its transient evolution resulted in a transient evolution in the temperature and combined velocity–temperature statistics of the thermal boundary layer, which are discussed below.

Figure 12. Evolution of ice subsurface temperature plotted as a function of non-dimensional time.

The evolutions of both $\theta _{rms}$ and $\langle v\theta \rangle$ when locally averaged in time were found to be non-increasing functions of time. They were plotted non-dimensionally by normalizing the locally averaged values evaluated at time $t$ by their value at the initial time ($t_{0}$). The non-dimensional evolutions of the local averages are plotted as a function of the non-dimensional time (using the melt time scale determined above) in figure 13.

Figure 13. Evolution of locally averaged non-dimensional RMS temperature (a–c) and sensible turbulent heat flux (d–f) as a function of the non-dimensional time: (a,d) $y/\delta = 0.13$, (b,e) $y/\delta = 0.33$ and (c,f) $y/\delta = 0.53$.

The non-dimensional evolutions of the local RMS temperature and sensible turbulent heat flux show that the nature of the decay of $\theta _{rms}$ and $\langle v\theta \rangle$ in time generally does not change as a function of $U_{\infty }$. For the non-dimensional RMS temperature fluctuations measured at $y/\delta = 0.13$, the highest-speed run does not fully collapse with the lower-speed runs, although the observed deviation from $0.1< t/\tau _{m}<0.3$ was temporary. A notable exception to the generally observed collapse is the lowest-speed run in figure 13(f), which did not appear to achieve a quasi-steady state at all. However, as alluded to previously, the values of $\langle v\theta \rangle$ were smallest for this run out of all nine cases studied herein, and approached the uncertainty of the measurements. The continued decrease in figure 13(f) is most likely an artefact of the exceptionally low magnitudes of the measured velocity and temperature fluctuations for this height and free-stream velocity ($\langle v\theta \rangle \sim 0.0003$ K m s$^{-1}$). This evolution can be attributed to the fact that, for very small temperature fluctuations, a cold-wire sensor can behave like a hot-wire sensor, becoming sensitive to velocity fluctuations as the signal-to-noise ratio of the temperature fluctuations decreases (Bruun Reference Bruun1995; Lemay & Benaïssa Reference Lemay and Benaïssa2001). In this limit, the cold-wire thermometer predominantly measures $u$ instead of $\theta$, and measurements of $\langle v\theta \rangle$ will effectively become measurements of $\langle uv \rangle$, thus explaining the the tendency of $\langle v\theta \rangle / (v\theta )_0$ in this one case to negative values, given that $u$ and $v$ are inversely correlated in a boundary layer ($\rho _{uv} \sim -0.4$; Tennekes & Lumley Reference Tennekes and Lumley1972).

Interestingly, the RMS temperature and sensible turbulent heat flux continue to decrease after $0.1\tau _{m}$, which corresponds to the point in time after which changes in the subsurface temperature were small (see figure 12). The continued change in the temperature-related statistics measured in the turbulent boundary layer above melting ice may imply that the surface temperature is continually changing while the subsurface temperature remains relatively constant, otherwise a plateau in the plots of $\theta _{rms}$ and $\langle v\theta \rangle$ would be observed in the region where $t>0.1\tau _{m}$.

The magnitude of the decrease in $\theta _{rms}$ and $\langle v\theta \rangle$ when compared with their initial values was similar for the two larger heights ($y/\delta = 0.33$ and 0.53), which typically decreased by about 60 % (i.e. ending at approximately 40 % of the initial value). However, closer to the surface at $y/\delta = 0.13$, they decreased by about 40 % (i.e. ending at approximately 60 % of the initial value). Furthermore, the amount of time needed to reach the quasi-stationary state decreased with increasing height, from approximately $0.4\tau _{m}$ at $y/\delta = 0.13$ to $0.15\tau _{m}$ at $y/\delta = 0.53$. In all cases, the period of time needed to reach this quasi-steady state was far beyond the beginning of the melting period (approximately $0.1\tau _{m}$). The point in time at which the pooling phase begins may be estimated using the amount of time needed to reach quasi-stationarity, wherein the thickness of the surface water layer would be sufficient to act as a heat sink as it develops an internal temperature gradient, thereby reducing the rate of change of the surface temperature and corresponding temperature measurements in the boundary layer. While, strictly speaking, the continued addition of heat and the continued melting of ice render this problem consistently non-stationary, it is still possible to make an assumption of quasi-stationarity. More specifically, this assumption is that the changing thickness of the water layer will not have a significant effect on the statistics measured in the boundary layer after a certain point. Within the scope of this work, it was not possible to determine the thickness of the water layer required to fulfil the assumption of quasi-stationarity (since the meltwater thickness was not continuously measured in time), though this would be an interesting opportunity for a future study.

Other non-dimensionalizations for the RMS temperature and sensible turbulent heat flux exist, and can be used to investigate how the changing quantities relate to the physics of the melting surface. The evolutions were mostly similar for different heights, such that only plots for the height $y/\delta = 0.33$ are presented. In figure 14(a), the locally averaged RMS temperature was scaled by the temperature difference $\langle T \rangle -T_{ss}$, where $T_{ss}$ is used as a proxy of the true surface temperature. The non-dimensional local RMS temperature shown in figure 14(a) increases until reaching a distinct peak at approximately $t=0.05\tau _{m}$. After this point, the local non-dimensional RMS temperature decays until the aforementioned quasi-steady state is achieved. The physical explanation of this peak is likely related to the use of $T_{ss}$ as a proxy for $T_{s}$. Once $T_{s}$ reaches 0 $^\circ$C, it will remain constant while the surface layer changes phase. However, the subsurface temperature will continue increasing, since it is not melting. After a brief period where $T_{s}$ is not changing in time, a thin water layer will form and $T_{s}$ will once again begin to increase. Since $\theta _{rms}$ is mostly dependent on $T_{s}$ (and not $T_{ss}$), this peak may be interpreted as a ‘phase lag’ between the subsurface and true surface temperatures. The location of this peak can be indicative of the beginning of the melt phase, consistent with the observations of $T_{ss}$ in figure 12, which place the beginning of the melt phase near $0.1\tau _{m}$. The difference between the $U_{\infty } = 3.2$ m s$^{-1}$ run and the other two lower-speed runs may be attributed to slight differences in the actual depth of the subsurface thermocouple inside the ice, as the behaviour at this speed nearly matches identically that of the other two speeds and could be easily corrected using a constant offset term. From figure 14(a), it is clear that all three speeds show similar evolutions, slopes and asymptotic behaviour, so it is not likely that there is an underlying physical difference in the nature of the RMS temperature at $U_{\infty } = 3.2$ m s$^{-1}$ compared with that at 1.1 and 2.0 m s$^{-1}$.

Figure 14. Non-dimensional evolutions of turbulent quantities at $y/\delta = 0.33$. (a) The RMS temperature and (b) sensible turbulent heat flux.

In figure 14(b), the sensible turbulent heat flux was scaled by the molecular heat flux (following the method of Yoon & Warhaft (Reference Yoon and Warhaft1990)), resulting in an effective, turbulent Nusselt number. In this case, if the temperature gradient $\beta ={\partial T}/{\partial y}$ is assumed to be linear, ${\partial T}/{\partial y}$ may be approximated by $({\langle T \rangle -T_{s}})/{y}$ (Fitzpatrick et al. Reference Fitzpatrick, Radić and Menounos2017). Since $T_{s}$ was not known precisely in the present work, it was estimated using $T_{ss}$, such that $\beta$ was approximated as $({\langle T \rangle -T_{ss}})/{y}$. The deviations from this linear assumption are strongest near the surface of the ice. In a constant flux layer, the friction temperature fluctuation ($\theta _{*}$) could be used to estimate the magnitude of the error using, for example, a Monin–Obukhov profile (see e.g. Denby & Smeets Reference Denby and Smeets2000). However, $\langle v\theta \rangle$ clearly changes with height in the present work. As such, there is no constant flux layer, and $\theta _{*}$ would not be constant. Therefore, a Monin–Obukhov profile cannot be used to estimate ${\partial T}/{\partial y}$ for this work. Furthermore, due to the transient nature of these experiments, profile measurements of the temperature field were not possible and we therefore resorted to a linear approximation to the vertical temperature gradient. Figure 14(b) depicts the evolution of the ratio of the turbulent to molecular heat flux. The results show a rapid increase until approximately $0.1\tau _{m}$, after which the evolution in time is minimal, again supporting the notion that the melt phase starts at approximately $0.1\tau _{m}$. Interestingly, the effects of the phase lag between $T_{ss}$ and $T_{s}$ are not observed in figure 14(b). This may be because the temperature gradient ($\beta$) will be determined by the velocity field and boundary conditions; however, the sensitivity of the non-dimensionalization of $\langle v\theta \rangle$ to the difference $T_{ss}-T_{s}$ is not as high as for the RMS temperature.

A final non-dimensionalization to consider for the sensible turbulent heat flux is the correlation coefficient ($\rho _{v\theta }={\langle v\theta \rangle }/{v_{rms}\theta _{rms}}$). Since it depends upon both $\langle v\theta \rangle$ and $\theta _{rms}$, its evolution serves to show how $\langle v\theta \rangle$ evolves relative to $\theta _{rms}$ (since $v_{rms}$ is effectively constant in time). In these experiments it was found that, in general, the correlation coefficient ($\rho _{v\theta }$) did not vary significantly in time. A representative plot of the correlation coefficient as a function of the non-dimensional time is presented in figure 15 for $y/\delta = 0.33$. Its stationarity is notable, because it shows that the evolution of the sensible turbulent heat flux will be proportional to that of the RMS temperature. The values of $\rho _{v\theta }$ observed in the present work generally compare well with those of Arya (Reference Arya1975) for similar levels of stratification ($\rho _{v\theta }\approx -0.3$).

Figure 15. Time dependence of the correlation coefficient ($\rho _{v\theta }$) for $y/\delta = 0.33$.

Lastly, we also note that the effects on the sensible heat fluxes of any humidity-induced buoyancy changes will be negligibly small. If one conservatively assumes that the free-stream air is dry (i.e. 0 % relative humidity) and that the air over the ice is saturated (100 % relative humidity), the difference in density of the air for the most sensitive case ($y=10$ mm and $U_{\infty }=1.1$ m s$^{-1}$; see table 5) is 0.7 %. This difference is an order of magnitude smaller than the density difference associated with a 20 $^\circ$C temperature difference between that of the ambient air (${\sim }20\,^\circ$C) and that of the ice surface (0 $^\circ$C), which is 7.3 %. This conclusion is further validated by estimating the convective Richardson number using the method of Yano (Reference Yano2023):

(3.14)\begin{equation} Ri_{c} = \frac{g}{\rho}\Delta\rho\frac{y}{\langle U \rangle^{2}}, \end{equation}

which results in $Ri_{c} \approx -1 \times 10^{-3}$, such that the effects of changing humidity on the convective processes (i.e. the sensible heat fluxes) are negligible, even for an assumed relative humidity difference of 100 %.

Table 5. Statistics of velocity and temperature over melting ice, averaged over a 10-minute period centred at $0.1\tau _{m}$.

Overall, the temperature and combined temperature–velocity statistics exhibited a transient evolution throughout the course of the experiments as the ice melted. The relative decrease of the locally averaged RMS temperature and sensible turbulent heat flux was found to be similar for all heights and speeds when non-dimensionalized by their initial values and plotted as a function of the non-dimensional time. Other non-dimensionalizations of the locally averaged RMS temperature and sensible turbulent heat flux showed that the melt period can be reasonably expected to begin near $0.1\tau _{m}$ for all speeds investigated in this work. In the following section, the temperature and combined velocity–temperature statistics will be evaluated over an interval of time in the melting phase in which they are approximately stationary, and analysed in a similar manner to the results of the previously discussed stationary velocity field.

3.4. Combined velocity-temperature statistics over melting ice

The previous subsection discussed the transient nature of the temperature and combined velocity–temperature statistics over melting ice over the course of the experiments. These statistics were measured over solid, warming ice ($T_{s}<0\,^\circ$C), phase-changing melting ice ($T_{s} = 0\,^\circ$C) and liquid water atop melting ice ($T_{s}>0\,^\circ$C). However, the period of principal interest for the present work is the melting phase of the experiments, when $T_{s} = 0\,^\circ$C (by definition) and there is only a thin film of water atop the ice.

Despite the transient nature of the experiments in which the temperature-related statistics are known to evolve in time, a ‘snapshot’ of the temperature and combined velocity–temperature statistics can be used to assess the nature of the temperature field over melting ice (and not warming ice or liquid water). To this end, a subset of the data containing $3.3 \times 10^{6}$ samples (measured over approximately 10 minutes) was used to provide insight on the temperature field over melting ice. In the previous section, it was demonstrated that the melt phase begins at around $0.1\tau _{m}$, and so the subsets were taken at this point in time for each of the nine runs in the test envelope. During this period, an assumption of statistical quasi-stationarity is necessary to facilitate a similar analysis to that of the velocity field (i.e. PDFs and spectra). To this end, the validity of this assumption was verified by calculating the change in $\theta _{rms}$ over the length of a subset, since it underwent a more significant change in time than $\langle v\theta \rangle$. It was found that the average change in $\theta _{rms}$ was only 4 %, approximately equal to the uncertainty in the velocity measurements, making the assumption of statistical stationarity over this time interval a sensible one. Given this, the statistics of velocity and temperature quantities measured during this period are presented in table 5.

To investigate the spectra and PDFs of the temperature field in a manner similar to that of the velocity field, the non-dimensional wavenumber spectra and PDFs of the temperature fluctuations are presented in figure 16. They are the most similar nearest the wall at $y/\delta = 0.13$, where the temperature gradient is the strongest. As the height is increased, there is more dissimilarity between the spectra for different speeds at the same height. The slope of the spectra also becomes less steep with increasing height, and it is difficult to observe an inertial range with $\kappa _{1}^{-5/3}$ scaling, indicating that the Reynolds (Mydlarski & Warhaft Reference Mydlarski and Warhaft1998) and Péclet (Lepore & Mydlarski Reference Lepore and Mydlarski2009, Reference Lepore and Mydlarski2012) numbers are not sufficiently high to satisfy the second hypothesis of KOC theory.

Figure 16. Non-dimensional wavenumber spectra (a–c) and PDFs (d–f) of the temperature fluctuations: (a,d) $y/\delta = 0.13$, (b,e) $y/\delta = 0.33$ and (c,f) $y/\delta = 0.53$.

Unlike the PDFs of the velocity fluctuations over melting ice, the PDFs of the temperature fluctuations over melting ice are clearly non-Gaussian. The PDFs are asymmetric due to the cooled nature of the flow such that temperatures must be bounded between $T_{\infty }$ and $T_{s}$. The PDFs of the temperature fluctuations are better mixed near the wall, unlike the velocity PDFs, which were less Gaussian near the wall. As the wall-normal distance is increased, the negative tail of the PDF increases. There is thus an increasing likelihood of measuring a ‘large’ temperature fluctuation relative to $\theta _{rms}$ with increasing wall-normal distance resulting from rare, cold fluid particles reaching those wall-normal locations.

Again mirroring the analysis performed for the velocity field, the coherence spectra between the vertical velocity fluctuations and the RMS temperature fluctuations ($H_{v\theta }$) were evaluated to quantify the spectral composition of the sensible turbulent heat flux over melting ice, and are plotted in figure 17. In all cases, the coherence between $v$ and $\theta$ increases with free-stream velocity (i.e. Reynolds number), as increased turbulent activity increases the turbulent heat transfer. There is a peak coherence at approximately $\kappa _{1}\eta = 10^{-3}$ for the highest free-stream velocity at all heights, which was not observed for the coherence spectra of $uv$. The existence of this peak implies that there is a preferential scale for the exchange of heat in the vertical direction that is not the largest scale of the flow. A similar peak was observed by Mestayer (Reference Mestayer1982) in his investigation of turbulent velocity and temperature statistics in a near-neutral turbulent boundary layer over cooled water, whose bulk Richardson number was approximately 0.0013. Compared to Mestayer's results, the coherence spectra are much lower in magnitude, potentially due to the comparatively lower Reynolds numbers in the present work. In the presence of strongly stable stratification, this peak might be attributable to internal gravity waves at the Brunt–Väisälä frequency ($N_{b}=\sqrt {{g\beta }/{\langle T \rangle }}$; Tennekes & Lumley Reference Tennekes and Lumley1972). In the present work, the ratio of the peak frequency to the Brunt–Väisälä frequency was always found to be less than 0.01, making this unlikely. This likelihood is further reduced when considering the weak stratification in the experiments of the present work.

Figure 17. Non-dimensional coherence spectra of vertical velocity and temperature: (a) $y/\delta = 0.13$, (b) $y/\delta = 0.33$ and (c) $y/\delta = 0.53$.

The overall magnitude of the coherence (and therefore the correlation coefficient $\rho _{v\theta }$ shown in table 5) was not found to scale predictably with height. This result is consistent with those of Arya (Reference Arya1975) and Ohya et al. (Reference Ohya, Neff and Meroney1997), who found that profiles of $\rho _{v\theta }$ exhibited a distinct peak at varying wall-normal distances. For similar levels of stratification to the experiments performed herein, they found this peak to occur around $y/\delta = 0.33$. This corresponds to that of the present work, in which similar values of $\rho _{v\theta }$ were found at this height to those of Arya (Reference Arya1975) and Ohya et al. (Reference Ohya, Neff and Meroney1997) (approximately $-$0.3).

The joint PDFs of the vertical velocity and temperature fluctuations were also computed. Similar to those of $u$ and $v$, joint PDFs of $v$ and $\theta$ were found to vary primarily with height and showed little dependence on the free-stream velocity. As such, the joint PDFs of $v$ and $\theta$ are only shown for the $U_{\infty }= 2.0$ m s$^{-1}$ case in figure 18. As expected given the PDFs of $v$ (which deviated from Gaussian behaviour in the tails) and $\theta$ (which were increasingly non-Gaussian with increasing height), none of the joint PDFs are joint-Gaussian, although the tendency to joint-Gaussian PDFs increases with height apart from a long tail associated with rare temperature fluctuations. Near the surface, the likelihood of observing a positive (warm) temperature fluctuation where $\theta >\theta _{rms}$ is rare, which explains the flat region in figure 18(a). As the height is increased, the likelihood of measuring a positive temperature fluctuation relative to $\theta _{rms}$ increases, as demonstrated by the rounding of the region in the top half of the plots in figure 18 as the joint PDFs approach a joint-Gaussian shape in the region where $\theta >0$.

Figure 18. Joint PDFs of the vertical velocity and temperature fluctuations for $U_{\infty }= 2.0$ m s$^{-1}$: (a) $y/\delta = 0.13$, (b) $y/\delta = 0.33$ and (c) $y/\delta = 0.53$.

However, it remains important to make the distinction between a ‘large’ fluctuation relative to $\theta _{rms}$ and a ‘large’ fluctuation as quantified by its magnitude. Since the RMS temperature decreases with height, and since the bulk temperature difference ($T_{\infty }-T_{s}$) is constant, the magnitude of a certain temperature fluctuation relative to the RMS temperature will also increase with height. For example, consider a warm fluctuation at the top of the outermost contour in figure 18(a), which will have a magnitude of approximately $1.5\times \theta _{rms}=1.4\,^{\circ }$C (corresponding to an absolute temperature of $T_{\infty }-1.1\,^{\circ }$C, which is the mean temperature deficit from table 5). The same fluctuation measured at $y/\delta = 0.53$ corresponds to approximately $11.8\times \theta _{rms}$, an exceptionally rare occurrence as shown by figure 18(c). This is a consequence of the bounded nature of the temperature field – the minimum temperature measurable is 0 $^\circ$C, whereas the maximum is $T_{\infty }$. Since $\langle T \rangle$ is closer to $T_{\infty }$ than $T_{s}$, there will be an asymmetry in the joint PDFs of $v$ and $\theta$, which skews towards $\theta <0$.

The joint PDFs also show that the largest temperature fluctuations are typically associated with smaller velocity fluctuations. Closer to the surface, larger velocity fluctuations (both positive and negative) are associated with smaller temperature fluctuations, implying that the vertical movement of air is mechanically driven by the turbulence as opposed to fluctuations in buoyancy. In other words, despite the larger temperature gradient closer to the surface, large vertical velocity fluctuations are still negatively correlated with the temperature fluctuations. The physical interpretation of this is that fast-moving parcels of air are bringing warm air downwards, as opposed to being brought towards (or away from) the surface due to their local temperature. The result is a joint PDF which appears horizontally elongated. As the height is increased, larger negative velocity fluctuations (i.e. movement towards the surface) are still associated with small temperature fluctuations. However, larger positive velocity fluctuations become associated with larger positive temperature fluctuations. Physically, this may be explained by parcels of air being quickly brought away from the surface by buoyant effects as opposed to mechanical effects, since the strength of the turbulence is known to decrease with increasing height (see table 5). Near the surface, fast upward motion appears to be more strongly associated with mechanical effects, but as the height increases and vertical motion generated by mechanical shear is weaker, there may be more warm parcels rising rapidly due to buoyancy effects as indicated by the growing region where $\theta >0$ and $v>0$.

The elongated tail of the joint PDFs in the region where $v>0$ and $\theta <0$ grows with increasing height. This region is associated with upward motion of cold air parcels, against the buoyancy gradient. The peak of this tail is close to $v = 0$, so the most intense temperature fluctuations are related to only small vertical velocity fluctuations. Nearest the surface, this may be interpreted as strong mechanical generation of turbulent motion overcoming the suppressing effects of stratification. However, with increasing height and decreasing strength of turbulent motions, it is difficult to make this same interpretation. Instead, the rare, intense cold fluctuations associated with small velocity fluctuations may indicate that at the largest height, the measurement probe sits at the edge of the thermal boundary layer and experiences only intermittent temperature fluctuations.

The PDFs of $v\theta$ were plotted non-dimensionally to evaluate the distribution of the sensible turbulent heat flux itself. They are presented in figure 19. The PDFs of $v\theta$ are more symmetric near the surface, but show an increasing negative skewness with increasing height.

Figure 19. Non-dimensional PDFs of $v\theta$ measured over melting ice: (a) $y/\delta = 0.13$, (b) $y/\delta = 0.33$ and (c) $y/\delta = 0.53$.

The computation of velocity and temperature statistics at equivalent non-dimensional times for all of the experiments facilitated further analysis of the relationship between the velocity and temperature fields to complement the analysis of the velocity statistics, as previously discussed. To understand how the intensity of the temperature fluctuations relates to increased anisotropy from the ice surface, the anisotropy ratio ($\phi =v_{rms}/u_{rms}$) is plotted as a function of the non-dimensional RMS temperature (herein referred to as the turbulent intensity of temperature, $Ti_{\theta } \equiv {\theta _{rms}}/{\langle T \rangle -T_{s}}$, where $T_{s}=0\,^{\circ }$C) in figure 20(a).

Figure 20. Anisotropy ratio plotted as a function of the turbulent intensity of temperature (a) and turbulent Richardson number (b), along with best-fit power laws ($\phi =0.224Ti_{\theta }^{-0.19}$ and $\phi =0.856Ri_{t}^{-0.20}$, respectively). Dashed lines represent a fit using all nine experiments, whereas the solid line in (b) represents a fit using only data measured at the $y/\delta = 0.13$ and 0.33 heights.

Figure 20(a) shows that the anisotropy ratio is closer to unity when the turbulent intensity of temperature is lower. This is not an unexpected result, given that the most dramatic anisotropy should be reasonably expected to be related to higher-intensity temperature field, which is observed at lower speeds and lower heights. However, figure 20(a) also shows that for the largest height ($y/\delta = 0.53$, triangular markers), the turbulent intensity of temperature is less than 1 %. At this height, the temperature fluctuations are exceptionally weak when compared with the mean temperature deficit ($\langle T \rangle -T_{s}$), such that they may be within the uncertainty of the measurements. This warranted another investigation of the relationship between the anisotropy ratio and the stratification (as quantified by the turbulent Richardson number).

To this end, the anisotropy ratio measured over melting ice is plotted as a function of $Ri_{t}$ in figure 20(b). Given the apparent collapse of the data at the two lower heights, it would appear that the (outlying) data at $y/\delta = 0.53$ are due to the fact that these data are in the outer, intermittent fringes of the thermal boundary layer. This conclusion is consistent with the values of $\rho _{v\theta }$ in table 5, which fall in the range of $-$0.05 to $-$0.018 for the $y/\delta = 0.53$ data set, and which are substantially lower than the values for $y/\delta = 0.13$ and 0.33, which are all approximately $-$0.3.

A dependence of the relationship between the anisotropy ratio on the turbulent Richardson number has only been previously observed in studies investigating the downstream evolution of stratified homogeneous turbulent flows (Yoon & Warhaft Reference Yoon and Warhaft1990; Piccirillo & van Atta Reference Piccirillo and van Atta1997), whereas the anisotropy ratio has been found to be relatively constant for different levels of stratification in boundary layer flows (Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017). Figure 20(b) validates the observed relationship between these two quantities previously noted in table 4. In this case, there is a clear relationship between the anisotropy ratio and turbulent Richardson number measured over melting ice for the cases where the turbulent intensity of the temperature field is greater than 1 %.

In addition to the observed Richardson number scaling of the anisotropy ratio for the experiments at $y/\delta = 0.13$ and 0.33, it was also of interest to investigate possible scalings of the sensible turbulent heat flux evaluated over melting ice. To this end, the sensible turbulent heat flux was non-dimensionalized by $\alpha \beta$ to form the turbulent Nusselt number ($Nu_{t}\equiv -\langle v\theta \rangle /\alpha \beta$) and plotted as a function of the Taylor-microscale Reynolds number. In the present work, where the flow is weakly stably stratified, the turbulent Nusselt number was also plotted as a function of the turbulent Richardson number. The plots are presented in figure 21, shown with corresponding best-fit power laws. Due to the low magnitude of the sensible turbulent heat flux at the highest measurement height ($y/\delta = 0.53$), two power laws were fitted to the data. One fit was applied to all nine cases in the test envelope (dashed lines), and the other was applied to only the results from the $y/\delta = 0.13$ and 0.33 heights (solid lines). The data for the $y/\delta = 0.53$ heights are significantly lower in magnitude than for the other two heights. The weak sensible turbulent heat flux at this height is presumably caused by external intermittency at the edge of the thermal boundary layer (Tennekes & Lumley Reference Tennekes and Lumley1972), as noted earlier, which will be ever more intermittent than that in a turbulent boundary layer with a laminar free stream, on account of the turbulent free stream generated by the active grid. The powers of the fits excluding the $y/\delta = 0.53$ data (solid lines) in figures 21(a) and 21(b) are 0.60 and $-$0.49, respectively. The exceptionally low values of the sensible turbulent heat flux at $y/\delta = 0.53$ observed in these plots are indicative of the measurement location's proximity to the edge of the thermal boundary layer. Due to the transient nature of these experiments, it was not possible to accurately measure the thermal boundary layer thickness ($\delta _\theta$) to confirm this conclusion. Moreover, classical analytical approaches to estimating the thermal boundary layer thickness by way of integral analysis for boundary layers with step changes in surface temperature (e.g. Reynolds, Kays & Kline Reference Reynolds, Kays and Kline1958; Kays, Crawford & Weigand Reference Kays, Crawford and Weigand2005) are not valid, given that these do not account for the (non-negligible) effect of the turbulent free stream. However, we approximated the thermal boundary layer thickness at $t = 0.1 \tau _m$ for the three free-stream velocities using the data in table 5 by assuming (i) a value of $\theta _{rms} = 0.05$ K in the isothermal free stream (which is typical for measurements in our laboratory, and includes contributions from electronic noise in the cold-wire thermometer and data acquisition system) and (ii) an exponential increase in $\theta _{rms}$ as the wall is approached (i.e. $\theta _{rms} = 0.05 \exp (-C(y/\delta _\theta - 1)))$. This is clearly invalid very close to the wall, but suitable for our range of measurements in which $0.13 \le y/\delta \le 0.53$. In this approximation, $C=C(U_\infty )$ is a constant that depends on the free-stream velocity. When $\delta _\theta$ (and $C$) are determined by least-squares regression, one obtains values of the thermal boundary layer thickness in the range $51.5 \le \delta _\theta \le 56.8$ mm (i.e. $0.69 \le \delta _\theta / \delta \le 0.76$). These estimates are consistent with measurements at the position farthest from the wall being near the edge of the thermal boundary layer, as hypothesized above, and thus subject to the effects of external intermittency. Moreover, these estimates give an indication of the development of the thermal boundary layer in the downstream direction.

Figure 21. Turbulent Nusselt number plotted as a function of the Taylor-microscale Reynolds number (a) and turbulent Richardson number (b). Dashed lines represent the best-fit power law using all nine experiments, whereas solid lines represent the best-fit power law using only the $y/\delta = 0.13$ and 0.33 heights.