Wind-tunnel measurements of sensible turbulent heat fluxes over melting ice

Abstract

Accelerated glacier melt and the loss of perennial snowfields have been associated with increased warming in polar regions, at rates up to four times faster than the rest of the world, thereby reinforcing the critical need for improved models (and predictions) of glacier melt. An essential requirement for such models is an improved understanding of the sensible heat fluxes over glaciers. Since their complexity makes them difficult to model, and direct measurements of sensible turbulent heat fluxes over real glaciers are both rare and impractical, the present work involves simultaneous hot-wire anemometry and cold-wire thermometry measurements of two components of velocity and temperature above a melting glacier model in a series of wind-tunnel experiments. Both single- and multi-variable statistics were used to compare the turbulent velocity field measured over melting ice with that of a similar flow in the absence of ice. The results demonstrate that the ice's presence reduces the magnitude of the Reynolds stresses and vertical velocity variance, but also increases the streamwise velocity variance. The transient evolution of temperature statistics throughout the melt process was also investigated and found to be similar when suitably non-dimensionalized. The velocity and temperature fields were furthermore evaluated at an equivalent non-dimensional time during the melt process, in which statistics of the temperature field, and joint statistics of the vertical velocity and temperature, were studied. The present work lays the foundation for future laboratory-scale replications of the flow above melting glaciers, and provides additional insight into turbulent heat transfer over melting ice.

1. Introduction

As the effects of climate change become increasingly relevant to the lives of people across the globe, the need for accurate climate prediction models becomes even more critical. The loss of perennial snowfields and accelerated glacier melt are associated with increased warming in polar regions, at rates up to four times faster than the rest of the world (Rantanen et al. 2022). Melting ice caps, which drain to the oceans, hold a significant potential contribution to global sea level increases (Radić et al. 2014), and meltwater run-off from mountainous glaciers holds significant influence over water availability in downstream areas extending far beyond the source (Bash & Marshall 2014).

To reliably predict the melting of glaciers and perennial snowfields, accurate models of the heat fluxes occurring between the glacier surface and environment are required. However, these interactions are extremely complex, involving different energy sources and sinks that govern the melt process, making them difficult to model accurately. One of the most common models for predicting glacier melt involves the use of a surface energy balance, which is given by (Hock 2005)

(1.1)\begin{equation} 0 = Q_{N} + Q_{H} + Q_{L} + Q_{G} + Q_{R} + Q_{M}. \end{equation}

It accounts for the net heat fluxes from solar radiation ($Q_{N}$), sensible heat fluxes from convection ($Q_{H}$), latent heat fluxes from evaporation and condensation at the surface ($Q_{L}$), heat flux to the ground via conduction ($Q_{G}$) and heat fluxes caused by rainfall ($Q_{R}$). Since they are both transported by turbulent velocity fluctuations in the air above the surface, the sensible and latent heat fluxes are together referred to as ‘turbulent heat fluxes’. The sum of the aforementioned components in the surface energy balance determines the total amount of energy available for melting a glacier's ice ($Q_{M}$).

To understand and predict glacier melt, it is crucial to accurately quantify each of the different heat fluxes in the surface energy budget. In real-world glacier environments, the dominant term is the net radiation term ($Q_{N}$), followed by the sensible heat flux term ($Q_{H}$), which can contribute up to 20 %–40 % of the melt energy (Denby & Greuell 2000). Moreover, the importance of the turbulent fluxes in the surface energy budget increases during situations where radiative fluxes are reduced, such as in cloudy conditions. Although the radiative heat fluxes are generally well understood and have adequate models (Fitzpatrick, Radić & Menounos 2017), there exists a demonstrated need for improved modelling of turbulent heat fluxes, as current methods can significantly under- or overestimate overall glacier melt by poor estimations of sensible heat fluxes, which tend to be much larger than latent heat fluxes (Suter, Hoelzle & Ohmura 2004). It also bears noting that because the models for latent heat fluxes are similar to those for sensible heat fluxes, improvements to sensible heat flux modelling will likely also improve latent heat flux modelling.

During the melt season, air in the atmospheric boundary layer is warmer than the glacier surface, leading to stable stratification in glacier flows (Hock 2005). Stably stratified boundary layers have been studied in great depth, with some of the most important work in the area of stably stratified flows being performed by Monin and Obukhov, who developed a similarity theory for the profiles of the mean velocity and temperature in a statistically stationary, stably stratified, boundary layer flow (Monin & Obukhov 1954). In a stably stratified flow, turbulent fluctuations are suppressed by the mean temperature gradient. Here, these fluctuations are defined using Reynolds averaging (Tennekes & Lumley 1972):

(1.2a,b)\begin{equation} U = \langle U \rangle + u, \quad V = \langle V \rangle + v, \end{equation}

where upper-case letters denote an instantaneous velocity, angular brackets denote averaging and lower-case letters denote fluctuations about the average, with the $U$ and $V$ velocities being in the longitudinal and transverse (i.e. vertical) directions, respectively.

Similarly, turbulent fluctuations of temperature ($\theta$) are defined relative to the instantaneous temperature ($T$) and average temperature ($\langle T \rangle$), such that

(1.3)\begin{equation} T = \langle T \rangle + \theta. \end{equation}

An important quantity in this theory of stably stratified flows is the Monin–Obukhov length scale, defined as (Tennekes & Lumley 1972; Nieuwstadt 1978)

(1.4)\begin{equation} L \equiv-\frac{\langle T \rangle u_{*}^{3}}{\kappa g\langle v\theta \rangle} = \frac{\langle T \rangle u_{*}^{2}}{\kappa g \theta_{*}}, \end{equation}

where $u_{*}$ is the friction velocity, $\kappa$ is the von Kármán constant, $g$ is the acceleration due to gravity, $\langle v\theta \rangle$ is the sensible turbulent heat flux in the vertical direction and $\theta _{*}$ is the temperature field analogue of the friction velocity, defined as $\theta _{*} \equiv -\langle v\theta \rangle / u_*$. In the case of stable stratification, in which $L>0$, the Monin–Obukhov length is used to scale the height to form the Monin–Obukhov stability parameter ($y/L$). In the surface layer, defined as the region in which $y/L<0.1$, $u_{*}$ and $\theta _{*}$ are constant, and as a result, the sensible turbulent heat flux is assumed to be constant. This represents a significant assumption (and limitation) of Monin–Obukhov theory, as it can be demonstrated that the constant-flux assumption is in the surface layer not always valid (Hock 2005).

This assumption has significant consequences for the modelling of atmospheric flows, since the majority of the turbulent activity occurs near the ground where the shear is highest (Pope 2000). This is especially important for atmospheric flows, wherein the surface layer can be several metres high due to the scale of the flow (Högström, Hunt & Smedman 2002). While Monin–Obukhov theory has been largely successful for predicting the mean velocity and temperature gradients in stably stratified flows (Salesky & Anderson 2018), it has been demonstrated that the constant-flux assumption in the surface layer is not always valid in the study of large-scale atmospheric flows (e.g. Mahrt 1998; Foken 2006).

There are numerous studies dedicated to the study of stably stratified turbulence in a wall-bounded flow, the most relevant of which to the present work are discussed below. One of the first of such studies was performed by Arya (1975), who measured turbulent velocity and temperature fluctuations in the presence of both stable and unstable stratification of varying intensities over a flat-plate boundary layer. The intensity, or strength, of the stable stratification in their work was quantified using the bulk Richardson number, defined as

(1.5)\begin{equation} Ri_{\delta} = \frac{g\delta}{T_{\infty}}\frac{T_{\infty}-T_{s}}{U_{\infty}^{2}}, \end{equation}

where $\delta$ is the boundary layer height, $T_{\infty }$ is the free-stream temperature, $T_{s}$ is the surface temperature and $U_{\infty }$ is the free-stream velocity. Larger bulk Richardson numbers indicate stronger, or more stable, stratification. Arya (1975) reported that turbulent quantities (such as intensities, fluxes and correlation coefficients) were sensitive to the effects of stratification, noting the suppression of turbulence in stronger stable stratification.

While not investigating a wall-bounded flow, Yoon & Warhaft (1990) performed wind-tunnel experiments to investigate the downstream evolution of grid-generated turbulence under stable stratification. Although these experiments were performed in the absence of mean shear, they found a similar destruction of the sensible turbulent heat fluxes with increasing Richardson number, and also observed counter-gradient fluxes under very high levels of stratification. Furthermore, they observed increased dissipation rates of scalar variance with increasing stratification. They defined a turbulent Richardson number as

(1.6)\begin{equation} Ri_{t} \equiv \frac{g\beta}{T_{\infty}}\left(\frac{\ell}{v_{rms}}\right)^{2}\end{equation}

to quantify the level of stratification, wherein $\beta$ is the temperature gradient ($\partial T/\partial y$) and $\ell$ is the integral length scale of the turbulence found using the autocorrelation of the longitudinal ($u$) velocity component.

Piccirillo & van Atta (1997) performed a similar study of the downstream evolution of grid-generated turbulence under stable stratification, in a uniformly sheared flow. They found that the critical Richardson number to distinguish between weak and strong stable stratification correlated well with the level of shear in the flow, such that increased shear levels corresponded to a larger critical Richardson number. Furthermore, they found that stably stratified turbulent flows were highly anisotropic at all scales, persisting down to the smallest scales of the turbulent cascade.

Ohya, Neff & Meroney (1997) performed a wind-tunnel investigation of the structure of turbulence in a stably stratified boundary layer for a wide range of bulk Richardson numbers (ranging from 0.12 to 1.33). They found that the magnitudes of turbulent intensities and fluxes of both momentum and heat were greatly reduced by the effects of stable stratification, and that the level of suppression correlated well with the bulk Richardson number. Furthermore, they found significant differences between weakly stable flows ($0< Ri_{\delta }<0.25$) and strongly stable flows ($Ri_{\delta }>0.25$), where 0.25 is the critical Richardson number found by a linearized theory for inviscid flow (Ohya et al. 1997). Under very stable stratification, the sensible turbulent heat fluxes tended to vanish, or even occurred against the buoyancy gradient. This work was then extended in Ohya (2001) to investigate the effects of roughness on the profiles of sensible turbulent heat fluxes in a stable boundary layer. They found that the presence of roughness elements tended to significantly increase the sensible turbulent heat fluxes for weakly stable cases when compared to smooth-wall experiments.

More recently, Williams et al. (2017) investigated the effects of stable stratification on turbulent boundary layer statistics, attempting to consider the effects of both the Reynolds number as well as the Richardson number. In their work, the best parameter for quantifying the stratification and differentiating between weak and strong stability cases was the gradient Richardson number, defined as

(1.7)\begin{equation} Ri_{g} \equiv \frac{g}{\varTheta_{\infty}}\frac{\partial\varTheta/\partial y}{(\partial U/\partial y)^{2}}. \end{equation}

They found that for weak to moderately stable flows, the sensible turbulent heat fluxes scaled with the wall shear stress. They suggested that this scaling was useful for differentiating between Reynolds and Richardson number effects. However, for higher levels of stratification, this shear stress scaling ceased to hold. Furthermore, they observed that the turbulence was preferentially damped in the outer region of the flow, maintaining its intensity near the wall, even for the case of strongly stable stratification. A consequence was that characteristics (such as the large-scale anisotropy) were unaffected by the effects of stable stratification up until a critical point. They argued that there was no change in the structure of turbulence under weak stratification, but under strong stratification (beyond the critical Richardson number) there was a ‘marked change in the flow response’.

Despite the well-developed literature surrounding stably stratified wall-bounded flows, there are nevertheless several important distinctions between real-world glacier flows and the stably stratified boundary layer flows described by Monin–Obukhov theory. One of these was investigated by Denby & Smeets (2000), who noted that the atmospheric boundary layer flow over glaciers is often dominated by katabatic (downslope) winds. They subsequently derived a momentum budget used to model katabatic glacier flows. Such flows typically originate from high-elevation air that cools enough to increase its density and travel downstream at high speeds. Simulated profiles of different quantities in a katabatic boundary layer were presented by Denby & Greuell (2000), using their derived energy and momentum budgets. In these katabatic flows, there is a distinct local velocity maximum that is lower in magnitude than the free-stream velocity found outside the boundary layer. The presence of this local wind speed maximum is inconsistent with the logarithmic profiles described by Monin–Obukhov theory, and is not captured in the previously discussed stably stratified studies, thereby limiting their applicability to glacier flows.

A detailed description of the katabatic flow over glaciers was then provided by Oerlemans & Grisogono (2002), who found that the katabatic layer (close to the surface) is characterized by a local wind speed maximum, typically within 10 m of the surface. The local wind speed maximum acts as a turbulent jet nested within the atmospheric boundary layer, resulting in wind speeds of the order of 6 m s$^{-1}$ for small-to-medium glaciers. The height and magnitude of the wind speed maximum served as characteristic length and velocity scales for the katabatic flow. Despite intense thermal gradients (up to 10 K m$^{-1}$) which act to stratify the flow, high levels of turbulence were observed in the katabatic boundary layer. These high levels of turbulence were found to dominate the exchange of energy and heat near the surface, easily overcoming the suppression from the stratification. Profiles of wind speed and temperature near the surface showed a nonlinear velocity profile and a linear temperature profile in the region between the surface and wind speed maximum. They explain that the temperature profiles are adequately modelled using the large-scale atmospheric boundary layer (since no temperature differential maximum occurs in the katabatic layer), but that velocity profiles require small-scale resolution of the order of the height of the wind speed maximum. Given their findings, a representative model of the katabatic boundary layer over a glacier (such as in a wind tunnel experiment) may be obtained using a highly turbulent background flow in the presence of a linear thermal gradient.

As a result of katabatic forcing, the flow over a glacier's surface can both be (weakly) stably stratified, while also containing high levels of background turbulence. Litt et al. (2015) studied the turbulence characteristics above a melting glacier using sonic anemometers and profile mast measurements made up to 6 m above the glacier's surface. They found that for low-speed wind conditions, the local wind speed maximum was approximately 2 m above the surface. However, for high-speed wind conditions, no wind speed maximum was observed. This would imply that the height of the local wind speed maximum was likely above the 6 m observational envelope of their study, and that its position scales with wind speed. They also reported exceptionally high turbulent intensities (as high as 50 % in some cases). For katabatic flows, they observed the typical turbulence intensity to be $32 \pm 12$ %, while the bulk Richardson number was around 0.09. In another related study, Litt et al. (2017) observed turbulence intensities to be around 10 %–15 %, although sometimes reaching as high as 50 %. The high turbulent intensities again highlight the dominance of the velocity field over the effects of stable stratification when considering heat transfer near the surface.

Another important factor differentiating the flow over glaciers from a typical stably stratified boundary layer is the topographical heterogeneity. Localized topographical variations can have significant effects on the near-surface flow (Mott, Vionnet & Grünewald 2018), including effects such as boundary layer decoupling (Mott et al. 2017) leading to the creation of stable internal boundary layers (Haugeneder et al. 2023) which influence the velocity and temperature fields above the surface. Nicholson & Stiperski (2020) compared the structures of turbulence over exposed and debris-covered glacier ice, finding strong similarity of the turbulent quantities for both surface types. Although there was strong aerodynamic similarity between the two cases, they found that the vertical temperature profiles were different, since debris covering the surface of glacier ice can influence its thermal and radiative properties, thereby influencing the temperature field.

Mott et al. (2016) performed a series of wind-tunnel experiments to investigate the effects of cold-air pooling observed above melting snow. In their experiments, they found that the effects of local topography features are significant for low to moderate wind speeds, where heat transfer at the snow's surface was strongly suppressed. However, at higher wind speeds (in this case, approximately 3 m s$^{-1}$), the difference between a flat and concave snow surface on the magnitude of sensible turbulent heat fluxes was negligible. Furthermore, they found that for more intensely stratified flows, the region above a flat snow patch in which turbulent mixing was suppressed was very shallow, perhaps suggesting that stratification effects may only be significant for a small region of the boundary layer, and only at lower Reynolds numbers.

A final important differentiation between the atmospheric flow over melting glaciers and a stably stratified boundary layer flow is related to the changing surface of a melting glacier. As the surface layer of ice melts, it changes phase from solid ice to liquid water. Most of the water is percolated through cracks in the glacier surface, potentially refreezing (Hock 2005) or runs off into rivers, streams, lakes, etc. (Bash & Marshall 2014). The stably stratified boundary layer over a thin film of liquid water may differ from that of an equivalent flow over a solid surface at the same temperature. Cohen & Hanratty (1965) performed experiments in which two- and three-dimensional surface waves were observed in a liquid film with a thickness of between 1.8 and 7.6 mm underneath a gaseous flow of sufficiently large velocity. Kendall (1970) found that the existence of surface waves will greatly influence the structure of the turbulence above a wavy (solid) surface. In the work of Sullivan, McWilliams & Moeng (2000), the momentum flux above a wavy surface was altered by up to 40 %, wherein an increase in surface wave speed was related to a decrease in the surface skin friction coefficient. In the presence of stable stratification, the ‘effective’ surface wave speed is increased, further decreasing the skin friction coefficient (Sullivan & McWilliams 2002). The depth of the liquid layer required to observe a significant effect on the turbulent statistics of the boundary layer remains to be seen. It has been observed that even in the case of very thin films, the liquid layer can influence the stability of the boundary layer (Náraigh et al. 2011), or potentially even lead to localized separation (Pelekasis & Tsamopoulos 2001).

The problem of understanding the nature of the turbulent velocity and temperature fields above melting glacier ice is unique in that the current (well-developed) areas of research do not adequately capture all of its distinct characteristics. Classical Monin–Obukhov theory and the results of studies of traditional stably stratified boundary layer flows have two key differences that hinder their application to the atmospheric flow over melting glaciers. First, Monin–Obukhov theory does not account for the presence of the local wind speed maximum commonly observed in real-world flows over melting glaciers. Second, typical studies of stably stratified boundary layer turbulence are performed over either solid or liquid surfaces, and therefore do not have the same surface conditions as a phase-changing boundary, such as melting ice. Moreover, field studies of glaciers, including those that use the most sensitive modern equipment, cannot adequately resolve the smallest (dissipative) scales of the turbulence (Bou-Zeid et al. 2010). As such, there is a demonstrated need for an improved understanding of weakly stratified turbulent flows over melting ice, in which the measurement equipment resolves the smallest scales of the turbulent cascade.

The present work is therefore motivated by the study of the atmospheric boundary layer in the region above a melting glacier, wherein there is a demonstrated need for an improved understanding of the velocity and temperature fields. This paper is intended to serve as a fundamental study into the nature of these fields above a melting ice surface by way of an investigation of the statistics of the longitudinal velocity, vertical velocity and temperature, simultaneously measured (with high spatial and temporal resolution) above such a surface in a series of wind-tunnel experiments to better understand the fluxes of momentum and sensible heat over melting ice. By way of these experiments, the present work attempts to answer the following questions:

1. (i) How do the statistics of the velocity field measured in the boundary layer above a melting ice compare with those measured in a similar flow in the absence of ice (i.e. over a solid surface)?

2. (ii) How do the statistics of velocity and temperature measured over a melting ice surface evolve in time as the surface warms and undergoes phase change?

3. (iii) How do the combined statistics of velocity and temperature measured over a melting ice surface change with surface height and free-stream velocity?

The present work represents a novel experimental investigation of the turbulent velocity and temperature fields above melting ice. The aforementioned questions are posed in the context of a flow that shares similarities to the atmosphere above a melting glacier: weakly stably stratified flow with high levels of background turbulence above a melting ice surface. However, it goes without saying that we are not attempting in this work to replicate a true glacier surface and environment in a laboratory-scale wind tunnel.

The remainder of this paper is organized as follows. In § 2, the wind-tunnel apparatus is described, along with the measurement techniques and experimental procedures. In § 3, the results are discussed in the following order: first, the characterization of the baseline flow in the absence of ice is presented; next the velocity statistics measured over melting ice are discussed, followed by an analysis of the transient evolution of the measured statistics over melting ice; lastly, the combined velocity–temperature statistics over melting ice are discussed. The experiment and results are summarized in § 4, in which conclusions are given and opportunities for future work are identified.

2. Apparatus and procedures

2.1. Wind-tunnel apparatus

The experiments in this work were performed in an open-circuit wind tunnel in the Aerodynamics Laboratory at McGill University. It is described in detail in Cohen (2019); however, a schematic is also presented herein in figure 1. To best emulate highly turbulent atmospheric glacier flow – which has turbulent intensities of the order of 10 % (Litt et al. 2015) – a wind tunnel with an active grid (Makita & Sassa 1991) was used to induce higher levels of background turbulence in these experiments. Both the wind tunnel and active grid used for this work were the same as those of Mydlarski & Warhaft (1996), which consisted of a $7 \times 7$ grid of rotating aluminium bars that were 6.4 mm in diameter. The mesh spacing, $M$, was 50.8 mm, or $1/8$ of the tunnel width.

Figure 1. Schematic of wind tunnel.

The existing test section was 2.8 m long, with a cross-sectional area of $0.407~{\rm m} \times 0.407$ m, and had a slightly diverging ceiling to offset the effects of boundary layer growth on the tunnel walls, thereby maintaining a constant centreline velocity and zero pressure gradient in the wind tunnel. It was followed by a section containing the glacier model, above which the experiments in this work were performed. To emulate a melting glacier inside the wind tunnel, an ice patch was constructed to sit flush with the tunnel floor. The glacier model was formed inside of a fibreglass tray with an inner length of 78 cm, width of 36 cm and depth of 20 cm. The bottom half of the tray was filled with styrofoam to provide insulation, and restrict most of the melting to the top surface in contact with the airflow. The insulated portion was partitioned from the water basin by placing a thin (0.5 mm thick) aluminium sheet on top of the styrofoam, and then sealing the edges using a marine-grade sealant. Halfway along the length of the glacier model, a type E thermocouple wire was embedded along the outer edge, approximately 5 mm below the surface, to estimate the surface temperature. The upper portion of the tray was then filled with water and frozen solid. The result was an ice block that was approximately 78 cm long, 36 cm wide and 10 cm deep.

Although glacier ice is typically characterized by surface roughness (Hock 2005), for the experiments in the present work (which represent an initial attempt to replicate glacial ice in a laboratory-scale wind tunnel), the ice surface was kept as smooth as possible to avoid the introduction of an additional parameter (surface roughness) into our investigation. A procedure describing the preparation of the ice surface is given in Harrison (2022). Since the existence of roughness elements in a turbulent boundary layer can change the nature of the turbulence (Tennekes & Lumley 1972), the objective of this work did not include the additional effects of surface roughness, and so the present investigations considered turbulent heat fluxes only over a smooth ice surface. Using a Vernier caliper and a bar of known thickness, measurements of surface roughness were taken at 20 evenly spaced points along the centreline of the ice surface. Surface roughness measurements were taken for three different glacier models, totalling 60 data points. The root-mean-square (RMS) roughness value ($k_{rms}$) was found to be 0.4 mm, and was subsequently used to determine whether the ice surface could be reasonably considered as an aerodynamically smooth surface (see § 3.1).

Another important consideration was the management of meltwater accumulated during the experiment. In a real glacier environment, meltwater can percolate through cracks in the glacier surface (potentially refreezing), or run off into streams and rivers (Bash & Marshall 2014), as noted earlier. In the present experiments the meltwater was not drained from the tray given that removing meltwater during the experiment would change the height of the ice surface, changing the height of the probe relative to the surface during an experiment. In the atmosphere, small changes in height associated with surface ablation take place over long time scales (several months). However, in a laboratory, the smaller length scales mean that small changes in height could change the results significantly. Furthermore, if the surface moves sufficiently below the top of the tray, there would be flow separation beyond the leading edge of the ice tray, changing the nature of the turbulent flow above the ice.

2.2. Measurement techniques

The velocity measurements were obtained in the present work using a Dantec Streamline Pro constant temperature anemometer (CTA) connected to a TSI 1241 X-wire probe, made of 5 $\mathrm {\mu }$m tungsten wire. The temperature measurements were obtained using a custom constant current anemometer made by Université Laval (Lemay & Benaïssa 2001), which was connected to a 1.25 $\mathrm {\mu }$m Wollaston wire mounted on a TSI 1210 single-wire probe. The analogue CTA signal was low-pass-filtered at 3 kHz using the onboard signal conditioner from the Dantec Streamline Pro CTA to remove electronic noise not caused by turbulent velocity fluctuations, then subsequently high-pass-filtered at 0.05 Hz using a Krohn-Hite 3384 band-pass filter. A gain was applied to the output signal of the high-pass-filtered hot-wire output, to optimize the analogue-to-digital conversion range. Both the low-pass- and band-pass-filtered signals from the hot-wires were recorded using a National Instruments BNC-2110 16-bit A/D board. The analogue constant current anemometer signal was also low-pass- and high-pass-filtered at 3 kHz and 0.05 Hz (respectively) using a Krohn-Hite 3384 filter. No gain was required for the filtered cold-wire signal. A Dantec Streamline thermistor temperature probe was placed upstream of the hot- and cold-wire assemblies (in a position so as to not disturb the flow) to record the ambient temperature throughout the course of the experiment. The hot- and cold-wire probes were calibrated using a TSI 1128 Air Velocity Calibrator, in which a controllable jet of compressed air was was cooled to sub-ambient temperatures by passing the compressed air line through an ice–water bath.

For the experiments over melting ice in the present work, data were acquired in blocks of 16 384 samples, measured at a frequency of 6000 Hz (i.e. twice the low-pass filter frequency) to satisfy the Nyquist criterion. The overall uncertainty associated with the velocity measurements, including the effects of the uncertainty in the temperature measurements, was found to be 4.1 % to achieve a 95 % confidence interval. For the temperature measurements, the uncertainty from the cold-wire was found to be approximately 0.97 % for a 95 % confidence interval (Harrison 2022).

2.3. Experimental procedures

Before making simultaneous measurements of velocity and temperature over melting ice, it was important to characterize the background flow in the tunnel in the absence of ice. Although no temperature fluctuations were expected to occur, simultaneous velocity and temperature measurements were nevertheless recorded. A substitute for the ice tray, composed of a geometrically identical tray with a sealed medium-density fibreboard lid, was inserted into the tunnel floor. The measurement probe was placed 3.6 m downstream of the active grid, wherein the boundary layer height, $\delta$ (defined as the height of the measurement probe above the tunnel floor at which the mean velocity was equal to 99 % of the free-stream velocity), was found to be approximately 7.5 cm. Vertical profiles of the local longitudinal and transverse velocities were measured for three different free-stream velocities (1.1, 2.0 and 3.2 m s$^{-1}$). These free-stream velocities were chosen to maximize the stratification (quantified using the Richardson number), while also being representative of the speeds typically observed in glacier flows (Mott et al. 2016). In this case, 1.1 m s$^{-1}$ was the slowest speed attainable in the wind tunnel, corresponding to a bulk Richardson number ($Ri_{\delta }$) of 0.04 which was similar to that observed by Oerlemans & Grisogono (2002) in real-world glacier flows.

For the main experiments in this research, statistics of velocity and temperature were measured over the surface of the glacier model as it melted. For each experiment, the glacier model was removed from the freezer and placed inside the glacier model test section. For these experiments, measurements were made at approximately 60 % of the total ice patch length. The unheated starting length (i.e. the distance between the grid and leading edge of the ice) was 3 m. Once it was aligned properly inside the tunnel, the measurement probe was positioned at either 10, 25 or 40 mm above the ice surface which correspond to $y/\delta =0.13$, $0.33$ and $0.53$. These heights were selected to cover a range of $y/\delta$ in the region below the wind speed maximum, but still above the ‘constant-flux’ surface layer predicted by Monin–Obukhov theory, and chosen such that the equivalent non-dimensional height could be reasonably attainable by those performing field observations in the atmospheric boundary layer over a glacier using commercially available weather stations, facilitating comparison of the values of longitudinal velocity and temperature measured in the present work with those measured by field scientists. For example, given a wind speed maximum height of 5 m over a glacier, $y/\delta =0.33$ corresponds to a height of 1.65 m. After successfully positioning the probe above the ice, the wind tunnel fan was turned on. Data collection was initiated once the thermocouple embedded in the ice reached $-10\,^\circ$C.

Due to the transient nature of the experiments, statistical moments and cross-moments of velocity and temperature were calculated by averaging over 50-block (approximately 2 minute) subintervals. This interval size was selected to ensure ‘instantaneous’ statistical convergence within any given time interval, while simultaneously avoiding transient changes caused by actual physical phenomena over the (4.5 hour) total duration of individual experiments. The subinterval duration was determined following an analysis of the convergence of the results in which 50 blocks was found to be the minimum sampling time required for convergence of the velocity statistics. Shortening this averaging window would result in unconverged statistics, while averaging over a longer period would unnecessarily smooth out transient effects.

3. Results and discussion

Given the complexity of the problem of understanding weakly stratified turbulent flows over melting ice, this section has been organized to highlight the findings from four different aspects of the present work.

First, the results of baseline experiments are used to characterize the non-stratified velocity field over a smooth surface in the absence of ice. The nature of the baseline velocity field is investigated to facilitate an improved comparison between the results obtained over melting ice.

Second, the velocity statistics measured over melting ice are then compared to the baseline cases. The purpose of this comparison is to enable a subsequent understanding of the combined effects of stratification and the melting ice surface (i.e. a different surface condition) on the velocity field. It was found that despite the continued melting of the ice surface throughout the experiments, the velocity statistics did not show a significant time dependence. This allows for a time-averaged comparison of the results for the two flows where the second-order statistics, single-component spectra and probability density functions (PDFs), coherence spectra and joint PDFs are considered.

Third, after investigating the stationary velocity field over melting ice, the evolution of the temperature statistics is considered. Unlike the velocity field, there are clear time dependencies of the temperature and combined velocity–temperature statistics because the ice surface continues to melt and both its temperature and phase evolve with time, as is discussed in detail in § 3.3. A ‘melt’ time scale characterized using the bulk (i.e. free-stream) quantities of the flow is used to collapse the evolution of the temperature-related statistics in time with considerable success. Evaluation of the transient evolution of the temperature field provides insight into the nature of the thermal boundary layer of the present experiments.

Lastly, using the results of the transient evolution of the temperature-related quantities, a representative ‘snapshot’ of the combined velocity–temperature field measured over melting ice is considered. Throughout the course of the experiments in the present work, the surface condition transitions between three stages. Initially the surface is solid ice with a subzero temperature. It then transitions to melting ice at 0 $^\circ$C. Finally, the surface consists of a thin layer of water warmer than 0 $^\circ$C atop solid ice. However, the primary region of concern for the study and understanding of melting glaciers is the second stage, wherein the solid ice is isothermally melting at 0 $^\circ$C. The aforementioned melt time scale is used to compare the results of the combined velocity–temperature statistics over melting ice at an equivalent point in (non-dimensional) time for all three free-stream velocities investigated herein. During this period, the combined velocity–temperature field can be reasonably assumed to be stationary, such that the calculation of coherence spectra and joint PDFs is feasible.

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. 1997). The shape factor ($H$) was calculated using the displacement thickness ($\delta ^{*}$) and the momentum thickness ($\theta ^{\delta }$) as follows (Pope 2000):

(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 1972).

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 1972) and the other using Prandtl's correlation for the skin friction in a turbulent boundary layer (Prandtl 1925). Details of the methods for estimating the friction velocity can be found in Harrison (2022).

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 1972) 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 2000). 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 1938) 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 1996). 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 2003). 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 (1990)) 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. 1997). 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 1975; Ohya et al. 1997; Williams et al. 2017), as buoyancy forces extract turbulent kinetic energy from the vertical fluctuations, thereby reducing the Reynolds stress (Arya 1975). 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 (1975) and Williams et al. (2017) 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. 2017). 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 1975). 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 1990; Piccirillo & van Atta 1997). 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. 1997), 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 2005).

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. 2017). The decrease in vertical velocity variance may then be attributed either to the effects of stratification through the classical mechanism (Tennekes & Lumley 1972) 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. 1994), and that greater surface shear stress can be expected over ice than over water (Overland 1985; Kantha & Mellor 1989).

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. 2017). 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. 2016). 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 1970), with impacts as high as 40 % (Sullivan et al. 2000). These effects have been observed for liquid films as thin as 1.8 mm (Cohen & Hanratty 1965). In the presence of stable stratification, the ‘effective’ surface wave speed is increased, further decreasing the skin friction coefficient (Sullivan & McWilliams 2002) 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. 2017); 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 1972) for passive scalar turbulence, as well as the $\kappa ^{-11/5}$ scaling for stably stratified turbulent flows (Bolgiano 1959; Obukhov 1959). 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 2017). This is because the effects of buoyancy are assumed to be negligible at smaller scales (Alam, Guha & Verma 2019). 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 2023), 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. (2016), 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 1999).

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. 2011). 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. 2017). 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 (2022) 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. 2011) 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 1995; Lemay & Benaïssa 2001). 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 1972).

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 (1990)), 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. 2017). 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 2000). 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 (1975) 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 (2023):

(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 1998) and Péclet (Lepore & Mydlarski 2009, 2012) 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 (1982) 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 1972). 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 (1975) and Ohya et al. (1997), 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 (1975) and Ohya et al. (1997) (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 1990; Piccirillo & van Atta 1997), whereas the anisotropy ratio has been found to be relatively constant for different levels of stratification in boundary layer flows (Williams et al. 2017). 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 1972), 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 1958; Kays, Crawford & Weigand 2005) 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.

4. Conclusions

The present work aimed to improve the understanding of the turbulent velocity and temperature fields over melting ice, while parametrizing the sensible turbulent heat flux. To this end, a series of turbulent flows representative of those typically seen over glaciers were reproduced in a laboratory-scale wind tunnel using an active grid to generate homogeneous isotropic turbulence with similar turbulence intensities and Richardson numbers to those observed in katabatic glacier flows. A novel glacier model was constructed to fit inside the wind tunnel. Hot-wire anemometry and cold-wire thermometry were used to make simultaneous measurements of two components of velocity and temperature, respectively, at three different heights in the turbulent boundary layer above the glacier model for three different free-stream velocities.

The statistics of the turbulent velocity field were measured over melting ice, and subsequently compared with those of the baseline flow (i.e. the same free-stream velocity, but in the absence of ice). It was found that the spectra of the streamwise (longitudinal) and vertical (transverse) velocity components did not change significantly between the experiments over ice and the baseline cases; however, there were small differences observable in the tails of the PDFs. Investigations of the joint PDFs also showed a strong similarity between the experiments over ice and the baseline cases, whereas the coherence spectra showed a distinct reduction in magnitude for the experiments over ice at the lowest height, with larger reductions in magnitude being observed with increasing Reynolds number. Furthermore, decreases in the vertical velocity variance and Reynolds stress were observed over melting ice, which may likely be attributed to boundary layer decoupling. However, an increase in the streamwise velocity variance was also observed, which has not been observed in studies of stably stratified turbulence over solid, non-melting surfaces. This increase may be attributed to the different nature of the boundary condition at the surface of the melting ice when compared to a solid surface. A future comparison is required between stably stratified flows over both solid (dry) and melting (wet) surfaces to better understand the implications of these phenomena on the turbulent velocity field over melting glaciers.

The transient evolutions of the temperature and combined velocity–temperature statistics measured throughout the melting process were also investigated. Both types of statistics were found to evolve similarly in time for different heights and speeds when non-dimensionalized using a time scale associated with the bulk convective heat transfer to the surface, decreasing in magnitude as the ice surface melted. After experiencing an initial decrease, it was observed that the level of stratification (as measured by the bulk and turbulent Richardson numbers) did not change significantly throughout the course of the experiments, even though the RMS temperature fluctuation and the sensible turbulent heat flux continued to evolve in time. This was due in part to the formulation of these Richardson numbers using the subsurface temperature as an approximation of the true surface temperature, but was also attributed to the mean temperature measured at the measurement height, which was found to be less sensitive to continued change in the surface temperature than the RMS temperature and sensible turbulent heat flux.

The statistics of the temperature and combined velocity–temperature fields over melting ice were then assessed during a shortened interval occurring at an equivalent non-dimensional time for all nine runs, in which an assumption of quasi-stationarity was made. PDFs of the temperature fluctuations exhibited non-Gaussian behaviour, unlike the PDFs of velocity which were found to be nearly Gaussian. Coherence spectra of the vertical velocity and temperature fluctuations were evaluated and found to have a distinct peak for the measurement height nearest the surface, comparing well with prior measurements over cooled water. Joint PDFs of the vertical velocity and temperature fields showed a general trend towards joint-Gaussian behaviour at larger heights, apart from an elongated tail of large, negative temperature fluctuations which were transported upward by positive velocity fluctuations. Investigations of the correlation between the anisotropy ratio and turbulent intensity of the temperature field, as well as correlations between the turbulent Nusselt number and Reynolds and Richardson numbers, showed that the largest measurement height in the present work was on the externally intermittent fringe of the thermal boundary layer.

While the present work contributes to the understanding of the turbulent velocity and temperature fields over melting ice and highlights important phenomena for the modelling of turbulent flows over melting glaciers, it also provides many opportunities for future work. With respect to the experimental apparatus, measuring and/or controlling the internal temperature gradient within the ice could be considered for a future study that focuses on the surface energy balance for a melting block of ice. With respect to the parameters investigated in a future work, the effects of stable stratification on the velocity field over melting ice can be further studied. In the present work, the velocity field was investigated for bulk Richardson numbers in the range of 0.001 to 0.05 and turbulent intensities in the range of 10 %–15 %, which are consistent with those observed above melting glaciers. However, further studies should be performed at both greater Reynolds numbers and Richardson numbers, to validate the observations of the second-order statistics over melting ice, as well as to determine the extent of the difference between the effect of stable stratification over solid versus melting surfaces. Furthermore, the addition of measurements at different downstream locations would provide valuable insight into the development of the thermal boundary layer, from non-equilibrium to equilibrium. The effects of surface roughness and heterogeneity could also be investigated to further link such research to actual glaciers. Moreover, an analysis of the contributions of humidity fluctuations/latent turbulent heat fluxes could also provide insight. Lastly, the present work can be extended to better model katabatic glacier flows in wind-tunnel experiments, insofar as it considers only the region below the wind speed maximum, which is the dominant region for convective heat transfer to the surface. In real-world flows, the temperature profile in this region is approximately linear. However, the temperature profile in the present work was closer to the traditional logarithmic shape found in turbulent boundary layers. A more accurate model of katabatic glacier flow could be achieved by use of a turbulent jet mounted inside the wind tunnel, at varying heights above the floor. The height of the jet would represent the wind-speed maximum, and may allow for a more accurate representation of katabatic glacier flow.