4.1. Overview

Figures 2 and 3 show the time evolution of the temperature and vorticity fields. Key features discussed below are marked with capital letters in the images. In addition, there are supplementary movies showing the time evolution of the flow field.

Figure 2.

Temperature $\theta$ in an $x$-$z$ plane during the evolution of the flow for Case 5. Slices were taken through the centre of the domain ($y = 0$). The colour scale varies in order to highlight the flow features. Some examples of these features are indicated: A – a weakly stratified turbulent region near the channel bottom; B – a strongly stratified region near the upper boundary in which the flow is close to laminar; C – Kelvin–Helmholtz-like instabilities in the stratified shear layer; D – downwelling thermal plumes; E – near-wall turbulence ejected into the central region by plumes as they reach the channel bottom. Supplementary movies showing the detailed time evolution of the temperature and vorticity fields for each of these time periods are available at https://doi.org/10.1017/jfm.2020.447.

Figure 3.

Vorticity magnitude $\omega$ in an $x$-$z$ plane during the evolution of the flow for Case 5. Slices were taken through the centre of the domain ($y = 0$). The colour scale is the same in all panels, with vorticity range 0 to 300. Some examples of important flow features are indicated: A – a weakly stratified turbulent region near the channel bottom containing hairpin-like vortical structures; B – a strongly stratified region near the upper boundary in which the flow is close to laminar; C – Kelvin–Helmholtz-like instabilities in the stratified shear layer; D – downwelling thermal plumes; E – near-wall turbulence ejected into the central region by plumes as they reach the channel bottom. See also supplementary movies showing the detailed time evolution of the temperature and vorticity fields for each of these time periods.

The initial temperature field in figure 2(a) shows that the channel is weakly stratified close to the bottom wall (A) and becomes progressively more strongly stratified towards the upper boundary (B). The initial state vorticity field shown in figure 3(a) exhibits characteristic features of turbulent channel flow such as hairpin-like vortical structures close to the bottom wall (A). Turbulence is noticeably damped with increasing height and the flow is essentially laminar in the strongly stratified region close to the upper boundary(B).

As the flow evolves the stratification breaks down as a result of the combination of shear from below and cooling from above. In KWAZ, we found that, for the case with no surface cooling, shear instabilities forming in the thermocline play an important role in the destratification process. These instabilities are similar to Kelvin–Helmholtz (K–H) waves, exhibiting the characteristic over-turns in the temperature field and braided cat's eyes in the vorticity field. Similar structures (C) are also clearly visible at a similar location ($z \approx 0.8$) in the visualisations of the flow with surface cooling at $t = 1$. At this time, thermal plumes (D) can also be seen starting to form close to the upper boundary. At $t = 2$ the shear instabilities are still apparent but are less prominent, while the thermal plumes have developed into larger structures. At $t = 3$ the thermal plumes are a dominant flow feature in the upper half of the channel. By $t = 4$ the stable stratification is essentially removed. Large plumes of cooler fluid are seen descending from the surface towards the channel bottom where they interact with the turbulence generated by shear at the channel bottom, ejecting it up into the central region of the channel (E). By this time the flow field is significantly more turbulent than the initial flow. By $t = 5$ the downwelling plumes (D) and ejections of wall turbulence (E) have expanded further as the flow approaches a new unstably stratified equilibrium state.

4.2. First- and second-order turbulence statistics

The time evolution of selected first- and second-order turbulence statistics for Case 5 are presented in figure 4. Statistics are shown for the initial equilibrium state, $t = 0$, and then at five subsequent times, $1 \le t \le 5$, during the flow evolution.

Figure 4.

Selected turbulence statistics during evolution of the flow for Case 5. Vertical profiles of (a) mean temperature $\langle \theta \rangle$ relative to the mean temperature at the channel bottom $\langle \theta _b \rangle$, (b) mean streamwise velocity $\langle u \rangle$, (c) a time series of $Ri_{\tau }$, (d–f) vertical profiles of turbulent shear stress $\langle u'w' \rangle$, root-mean-squared vertical velocity fluctuation $w^{rms} = \langle w'w' \rangle ^{1/2}$ and turbulent heat flux $\langle \theta 'w' \rangle$. Here, $\langle u'w' \rangle$ and $w^{rms}$ are renormalised in terms of the friction velocity $\tilde {u}_{\tau }(\tilde {t})$ measured at the time at which the statistic is calculated. $\langle \theta 'w' \rangle$ is normalised in terms of the surface cooling flux$q_s$.

As discussed in § 2, the friction velocity $\tilde {u}_{\tau }$ varies during the destratification process. For many statistics it is useful to renormalise in terms of the friction velocity measured at the time at which the statistic is calculated. Thus, in figure 4 for example, $\langle u'w' \rangle$, which is normalised in terms of $\tilde {u}_{\tau ,0}$, is presented in the form $\langle u'w' \rangle / u_{\tau }^2$. Since $u_{\tau }(t) = \tilde {u}_{\tau }(\tilde {t})/\tilde {u}_{\tau ,0}$, this has the effect of renormalising $\langle u'w' \rangle$ in terms of $\tilde {u}_{\tau }(\tilde {t})$. In some cases, such as $\langle u \rangle$ and $\langle \theta 'w' \rangle$, we chose not to renormalise in this way because the arguments we wish to make are made more clearly without renormalising. We refer the reader to KWAZ for a more detailed discussion of this approach.

The destratification of the channel is clearly apparent in the profiles of mean temperature in figure 4(a). Surface cooling reduces the temperature close to the upper boundary leading to a region of unstable stratification that gradually increases in depth over time. Shear driven mixing from below acts to reduce the temperature gradient in the stably stratified lower part of the channel. By $t \approx 4$ the initial stable stratification has been removed. This can also be seen from the time series of $Ri_{\tau }$ in figure 4(c), which reaches a value of $Ri_{\tau }=0$ at approximately $t \approx 4$. At this time the destratification process is considered complete. At $t = 5$ the flow is unstably stratified over the entire channel height and is evolving toward a new unstably stratified equilibrium state.

The reduction in turbulent mixing in the initial state due to strong stratification close to the upper boundary is apparent in the profile of $\langle u'w' \rangle$ in figure 4(d), which shows a reduction in magnitude in the near-surface region at $t = 0$. As a result the shear stress required to balance the streamwise pressure gradient in this region must be provided predominantly by viscous shear. This leads to the inflected velocity profile seen in figure 4(b), which shows a substantial increase in the mean streamwise velocity close to the upper boundary relative to the final unstably stratified flow.

As the flow evolves there are significant changes to the shear stress profile. At $t = 1$ there is an enhancement in $\langle u'w' \rangle$ from the upper boundary down to $z \approx 0.6$ with a peak at $z \approx 0.82$. By $t = 2$ the enhancement has increased substantially. The peak has descended to $z \approx 0.6$ and there is a noticeable increase in $\langle u'w' \rangle$ all the way down to $z \approx 0.05$, close to the channel bottom. This enhancement is much greater than that seen for the flow with no surface cooling (see KWAZ, figure 6d). From figure 4(e) it can be seen that the enhancement in turbulent shear stress for $t = 1$ and $t = 2$ corresponds to similar enhancements in the vertical velocity fluctuations. As plumes of cooled fluid descend they increase the vertical velocity fluctuation. Interaction of the downwelling plumes with the streamwise velocity fluctuations leads to augmentation of the turbulent shear stress. This enhancement extends downward as surface plumes penetrate deeper into the channel. The profiles for $3 \le t \le 5$ show a region of near constant $\langle u'w' \rangle$ for $0.05 < z < 0.4$ where the plumes interact with the near-wall shear-generated turbulence. In particular, as the plumes approach the bottom of the channel they displace the locally generated turbulence and eject it up into the flow, increasing the turbulent shear stress in the region above the channel bottom. This is seen in the flow visualisations in figure 3 (regions marked E) and more clearly in the flow visualisations showing transverse sections through the flow field in figure 5. Interestingly, at $t = 4$ the magnitude of the turbulent shear stress in the region $0.05 < z < 0.4$ exceeds the wall shear stress. This imbalance is due to the transient nature of the flow as it is evolving in response to the change in conditions. By $t=5$ the flow is evolving toward its new equilibrium state and the balance between wall stress and shear stresses in the near-wall region is restored.

Figure 5.

Temperature and vorticity on a transverse section ($y$-$z$ plane) at $t = 5$ for Case 5. Slices were taken through the centre of the domain ($x=0$) looking upstream. Key flow features are marked: D – downwelling thermal plumes; E – near-wall turbulence ejected into the central region by plumes as they reach the channel bottom.

The turbulent heat flux $\langle \theta 'w' \rangle$ shown in figure 4(f) also undergoes substantial changes. In the initial heated equilibrium state the net downwards turbulent heat flux (negative $\langle \theta 'w' \rangle$) balances the radiative heat source. When the flow switches from internal heating to surface cooling at $t = 0$ this balance is disturbed. As the flow destratifies a downwards heat flux continues in the lower, stably stratified region of the channel. Here, as with the turbulent shear stress, descending thermal plumes interact with temperature fluctuations, generating a descending region of strongly enhanced downwards turbulent heat flux at $t = 1$ and $t = 2$. The heat flux profile differs from the shear stress profile, however, because this region of downwards heat flux is located below a region of upwards heat flux (positive $\langle \theta 'w' \rangle$) imposed by the cooling flux applied at the upper boundary. The region of upwards heat flux gradually extends into the channel until by $t = 5$ there is an approximately linear profile from $\langle \theta 'w' \rangle \approx q_s$ close to the upper boundary to $\langle \theta 'w' \rangle = 0$ at the adiabatic channel bottom. This is similar to the heat flux profiles seen in other asymmetric free convection flows (for example, D'Asaro, Winters & Lien Reference D'Asaro, Winters and Lien2002; Zikanov, Slinn & Dhanak Reference Zikanov, Slinn and Dhanak2002).

4.3. Turbulent kinetic energy and temperature variance budgets

Figure 6 shows the time evolution for Case 5 of the dominant terms in the turbulent kinetic energy budget along with turbulent kinetic energy itself. For this flow, which is homogeneous on $x$-$y$ planes, the turbulent kinetic energy equation can be written as

(4.1)\begin{equation} \frac{\partial k}{\partial t} = \frac{\partial}{\partial z} \left[-\frac{1}{2} \langle w'u'_iu'_i \rangle - \langle w'p' \rangle + 2\nu\langle s_{i3}u'_i \rangle \right] - \langle u'w' \rangle \frac{\partial \langle u \rangle}{\partial z} + \gamma \langle \theta'w' \rangle - 2\nu\langle s_{ij} s_{ij} \rangle, \end{equation}

where turbulent kinetic energy $k = 1/2 \langle u'_iu'_i \rangle$ and $s_{ij}$ is the strain rate due to velocity fluctuations given by

(4.2)\begin{equation} s_{ij} = \frac{1}{2}\left(\frac{\partial u'_i}{\partial x_j} + \frac{\partial u'_j}{\partial x_i}\right). \end{equation}

The dominant terms in our flow are: shear production $P=-\langle u'w' \rangle {\partial \langle u \rangle }/{\partial z}$, downwards buoyancy flux $B = -\gamma \langle \theta 'w' \rangle$, dissipation rate $\varepsilon = 2\nu \langle s_{ij} s_{ij} \rangle$ and transport due to turbulent fluctuations $T = -1/2{\partial \langle w'u'_iu'_i \rangle }/{\partial z}$.

Figure 6.

Vertical profiles of turbulent kinetic energy budget terms for Case 5: (a) shear production, (b) dissipation rate, (c) downwards buoyancy flux, (d) turbulent transport and (e) turbulent kinetic energy. All terms are renormalised with respect to the friction velocity $\tilde {u}_{\tau }(\tilde {t})$ measured at the time at which the statistic is calculated.

During the first half of the destratification process (up to $t = 2$) there is substantial enhancement of $P$, $B$ and $\varepsilon$ in the upper two thirds of the channel, with peaks at $z \approx 0.82$ at $t = 1$ that descend to $z \approx 0.6$ to 0.7 at $t = 2$. The peaks for $B$ are at a slightly lower height than those for $P$ and $\varepsilon$ due to the region of negative $B$ (upwards buoyancy flux) above. In the initial equilibrium state, as can be seen in the mean velocity profile in figure 4(b), there is significant shear present between $z=0.5$ and 0.9. As described in KWAZ, removal of the heat source allows this region to destabilise, leading to the formation of Kelvin–Helmholtz-like instabilities. For the case without surface cooling, these instabilities were found to be associated with an increase in turbulent mixing, leading to increased $P$, $B$ and $\varepsilon$ in this region. For the current case of flow with surface cooling similar K–H-like instabilities are also present, as is evident from the flow field visualisations (regions marked C in figures 2 and 3). In this case, however, the increases in $P$, $B$ and $\varepsilon$ are much more substantial. This difference can be explained by the fact that, as seen above, when surface cooling is present the thermal plumes generate vertical velocity fluctuations that cause substantial enhancements in the turbulent shear stress and heat flux, which in turn contribute to the shear production and downwards buoyancy flux. Thus, for the flow with surface cooling, enhanced mixing is due to a combination of K–H-like shear instabilities and downwelling thermal plumes. The interactions between these flow features is most clearly seen in the supplementary movies.

Enhancement of shear production by thermal plumes through enhancement of the turbulent shear stress occurs only while the vertical mean velocity gradient is maintained. The reduction in $P$ for $z > 0.8$ at $t = 2$ relative to $t = 1$ is due to a reduction in the mean velocity gradient in this region (see figure 4b) as a result of turbulent mixing by the evolving convective instabilities. As the plumes penetrate deeper the region over which shear production is attenuated extends downwards. In the near-wall region ($z < 0.1$), however, (see inset in figure 6a) shear production increases in the later stages of the flow evolution ($t > 3$) as the plumes interact strongly with wall turbulence and shear.

The imposition of surface cooling to a system that is initially stably stratified leads to a profile for the buoyancy flux, $B$, that changes sign within the channel. In the upper part of the channel the buoyancy flux is upwards ($B < 0$) indicating net transfer of energy from potential energy to turbulent kinetic energy. In the lower part of the channel the buoyancy flux is downwards ($B > 0$), indicating net transfer of energy from turbulent kinetic energy to potential energy.

Profiles of dissipation $\varepsilon$ and $k$ follow similar trends to shear production. Globally, there is an imbalance between $P$, $B$ and $\varepsilon$ so that by the end of the flow evolution there is an increase in turbulent kinetic energy $k$ in the channel.

The profiles of turbulent transport $T$ are quite noisy; however, the profiles at $t=1$ and $t=2$ indicate regions of negative $T$ around $z \approx 0.8$ and $z \approx 0.7$ respectively at these times, which may represent downwards transport of turbulent kinetic energy by descending thermal plumes. At $t = 2$ there is also a region of positive $T$ around $z \approx 0.9$ which indicates transport of turbulence from the more energetic central region of the channel up into the now unstably stratified region close to the top of the channel.

Figure 7 shows the transient response of dominant terms in the transport equation for temperature variance $\langle \theta '^2 \rangle$. For our flow the temperature variance equation can be written as

(4.3)\begin{equation} \frac{\partial \langle \theta'^2 \rangle}{\partial t} = \frac{\partial}{\partial z} \left[ -\langle w'\theta'^2 \rangle + \kappa\frac{\partial \langle \theta'^2 \rangle}{\partial z} \right]- 2\langle \theta'w' \rangle \frac{\partial \langle \theta \rangle}{\partial z} - 2\kappa\left\langle\frac{\partial \theta'}{\partial x_j}\frac{\partial \theta'}{\partial x_j}\right\rangle. \end{equation}

The dominant terms are production $P_{\theta }=-2\langle \theta 'w' \rangle {\partial \langle \theta \rangle }/{\partial z}$, turbulent transport $T_{\theta } = -{\partial \langle w'\theta '^2 \rangle }/{\partial z}$ and dissipation rate $\chi = 2\kappa \langle (\partial \theta '/\partial x_j)^2 \rangle$.

Figure 7.

Vertical profiles of temperature variance budget terms for Case 5: (a) production, (b) dissipation rate, (c) turbulent transport and (d) temperature variance. All terms are renormalised with respect to the friction velocity $\tilde {u}_{\tau }(\tilde {t})$ measured at the time at which the statistic is calculated.

In the initial flow, temperature variance is produced predominantly in a region centred around the thermocline at $z \approx 0.8$, where both the vertical temperature gradient and turbulent heat flux are high (see figure 4). As a result, $P_{\theta }$, $\chi$ and $\langle \theta '^2 \rangle$ have a similar distribution.

Imposition of surface cooling adds a new source of temperature variance at the upper boundary, leading to a rapid increase in $P_{\theta }$ just below the surface. This, combined with transition of the near-surface region to turbulent flow, leads to a substantial increase in $\langle \theta '^2 \rangle$ close to the upper boundary at $t = 1$. The increased intensity of turbulence around $z = 0.8$ also leads to an increase in $P_{\theta }$ in this region at $t = 1$; however, this is balanced by an increase in $\chi$ and transport $T_{\theta }$ out of this region ($0.7 < z < 0.85$) into the region above ($0.85 < z < 0.95$), so that $\langle \theta '^2 \rangle$ remains relatively constant for $z < 0.8$ at this time, while it increases for $z > 0.9$. In the region $0.9 < z < 0.95$ the temperature gradient is close to zero so $P_{\theta }$ becomes very small, and in fact becomes negative, indicating a small region of counter-gradient turbulent heat flux at $t=1$.

Similar dynamics is in play at $t = 2$; however, the thermal plumes have now penetrated deeper into the channel. For $t \ge 3$ the mean temperature gradient in regions below $z = 0.9$ is very small and production of temperature variance $P_{\theta }$ is confined predominantly to the region close to the upper boundary. The profiles of $T_{\theta }$ show that downwelling plumes transport temperature fluctuations out of this region into the lower regions leading to a temperature variance profile that gradually decreases with increasing depth in the late stages of the flow evolution.

4.4. Local flow parameters and length scales

In this section we present vertical profiles of gradient Richardson number $Ri$, shear Reynolds number $Re_S$, buoyancy Reynolds number ${\textit {Re}}_b$ and related turbulence length scales, the Kolmogorov $\eta$, Corrsin $l_C$ and Ozmidov $l_O$ length scales. These parameters and length scales can be used to characterise local flow conditions in stratified turbulence (see for example, Chung & Matheou Reference Chung and Matheou2012). All parameters and length scales are non-dimensionalised in terms of $\tilde {h}$, $\tilde {u}_{\tau ,0}$ and $\tilde {\Theta }_{N}$ (see (2.15)).

Figure 8(a) shows the time evolution of vertical profiles of the gradient Richardson number. The gradient Richardson number $Ri = N^2/S^2$ gives an indication of the stability of the flow with respect to density stratification characterised by buoyancy frequency $N = (\gamma \, \partial \langle \theta \rangle / \partial z)^{1/2}$, and velocity shear $S = \partial \langle u \rangle / \partial z$. Alternatively, $Ri$ can be written in terms of the ratio between the Corrsin $l_C$ and Ozmidov $l_O$ length scales, $Ri \equiv ({l_C}/{l_O})^{4/3}$, and hence can be interpreted as the degree of separation between scales above which turbulence is strongly affected by shear and buoyancy (Chung & Matheou Reference Chung and Matheou2012).

Figure 8.

Selected non-dimensional parameters and length scales characterising stratified turbulence for Case 5: (a) the gradient Richardson number $Ri = ({l_C}/{l_O})^{4/3}$, (b) the shear Reynolds number $Re_S = ({l_C}/{\eta })^{4/3}$, (c) the local buoyancy Reynolds number ${\textit {Re}}_b = ({l_O}/{\eta })^{4/3}$, (d) the Kolmogorov scale $\eta =({\nu ^3}/{\varepsilon })^{1/4}$, (e) the Corrsin scale $l_C = ({\varepsilon }/{S^3})^{1/2}$ and (f) the Ozmidov scale $l_O = ({\varepsilon }/{N^3})^{1/2}$.

The profile for $Ri$ in the initial equilibrium state (figure 8a) shows three distinct regimes spanning three regions across the channel. For $0 < z < 0.5$, $Ri$ increases from $Ri \approx 0$ at the channel bottom to a value of $Ri \approx 0.18$. In this region $l_C$ (figure 8e) is relatively constant except in the region very close to the channel bottom, while $l_O$ (figure 8f) decreases with height due to increasing stratification. For $0.5 < z < 0.75$, $l_C$ and $l_O$ both decrease at approximately the same rate with increasing height, so that $Ri$ is approximately constant at what appears to be a critical value of $Ri_c \approx 0.18$. Comparable values of $Ri_c$ have been reported for a range of stratified shear flows including the following: stratified channel flow (Garcia-Villalba & del Alamo Reference Garcia-Villalba and del Alamo2011), stationary homogeneous stratified sheared turbulence (Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000; Chung & Matheou Reference Chung and Matheou2012), stratified plane Couette flow (Zhou et al. Reference Zhou, Taylor and Caulfield2017) and the stably stratified atmospheric boundary (Basu & Porté-Agel Reference Basu and Porté-Agel2006). For $0.75 < z < 1$, $l_O$ decreases more rapidly than $l_C$, so that $Ri$ increases significantly. As discussed above with respect to the flow visualisations in figures 2(a) and 3(a), above this height turbulence is intermittent and the flow becomes essentially laminar close to the upper boundary.

As the flow evolves, the profiles of $Ri$ show an expanding region of negative $Ri$ close to the upper boundary, which represents unstable stratification resulting from the formation of convective plumes due to the surface cooling. At any time, the point at which $Ri$ changes sign delineates the boundary between the unstably and stably stratified regions of the flow. When no surface cooling is present, KWAZ found that the value of $Ri$ in the central region of the flow remains constant at $Ri \approx 0.18$ until the late stages of the flow evolution. In contrast, the profiles shown here for the case with surface cooling show a progressive reduction in $Ri$ in the central region of the channel, which can be attributed to the effect of thermal plumes penetrating into this region. The extra turbulent mixing reduces both temperature and velocity gradient; however, the fact that $Ri$ decreases indicates that the temperature gradient is reduced more rapidly.

Figure 8(b) shows the time evolution of vertical profiles of the shear Reynolds number $\textit {Re}_S = {\varepsilon }/{\nu S^2}$, which can also be written as a function of the ratio of Corrsin and Kolmogorov length scales, $\textit {Re}_S = ({l_C}/{\eta })^{4/3}$, and hence gives an indication of the extent of the inertial subrange (Chung & Matheou Reference Chung and Matheou2012).

In the initial equilibrium state, $\textit {Re}_S$ also shows multiple flow regimes across the height of the channel resulting from changes in $l_C$ and $\eta$. For $0 < z < 0.1$, shear at the channel bottom drives $l_C$ and $\textit {Re}_S$ toward zero and viscous effects dominate. In the weakly stratified turbulent boundary layer, $0.1 < z < 0.5$, $\textit {Re}_S$ is in the range 10–20, indicating approximately an order of magnitude separation between $\eta$ and $l_C$, both of which are approximately constant in this region. For $0.5 < z < 0.8$, $\eta$ remains approximately constant while both $l_C$ and $\textit {Re}_S$ decrease as a result of increasing shear (see velocity profile in figures 4a and 4b). In the region $0.8 < z < 0.95$, $\eta$ increases as significant damping of turbulence by stratification decreases viscous turbulence dissipation $\varepsilon$. Meanwhile $l_C$ continues to decrease, so that $Re_S$ drops below 1. Very close to the upper boundary, $0.95 < z < 1$, shear decreases again leading to an increase in $l_C$ and $Re_S$.

As the flow evolves, the main change seen is a significant increase in $Re_S$ due to a reduction in shear and hence an increase in $l_C$. This process starts in the region close to the upper surface as thermal plumes break down the stably stratified shear layer, replacing it with a strongly mixed, convectively unstable region in which shear is close to zero. The region of increased $l_C$ and $Re_S$ expands downwards over time as the thermal plumes penetrate further into the channel.

Figure 8(c) shows the time evolution of vertical profiles of the local buoyancy Reynolds number ${\textit {Re}}_b = {\varepsilon }/{\nu N^2}$, which can also be written as a function of the ratio of Ozmidov and Kolmogorov scales ${\textit {Re}}_b = ({l_O}/{\eta })^{4/3}$, giving a measure of the degree of separation between the scale above which fluid motion is influenced by stable buoyancy stratification and the smallest scales of motion (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007).

The buoyancy Reynolds number has received considerable attention in the literature due to its importance in geophysical measurements and modelling, where it is used in parameterisations for eddy viscosity (for example, Osborn Reference Osborn1980; Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005), and mixing efficiency (for example, Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014; Walter et al. Reference Walter, Squibb, Woodson, Koseff and Monismith2014; Scotti & White Reference Scotti and White2016; Ivey et al. Reference Ivey, Bluteau and Jones2018). In KWAZ we give an in-depth discussion and comparison of a number of these parameterisations with results obtained in our stratified channel flow simulations, and show that parameterisations developed for the ocean and stably stratified atmospheric boundary layer also accurately predict turbulence properties in the central region of stably stratified open-channel flow. In the following, we focus primarily on features resulting from the introduction of surface cooling.

Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) define three regimes for ${\textit {Re}}_b$ in stably stratified turbulent shear flows: a diffusive regime for ${\textit {Re}}_b < 7$ in which turbulence is strongly damped and viscous effects dominate, an intermediate regime $7 < {\textit {Re}}_b < 100$ in which turbulence is energetic but significantly affected by stratification, and an energetic regime, ${\textit {Re}}_b > 100$, in which the effects of stratification become progressively weaker as ${\textit {Re}}_b$ increases. In the initial equilibrium state for our simulations ${\textit {Re}}_b$ covers all three regions described above, with the energetic regime seen for $0 < z < 0.4$, the intermediate regime for $0.4 < z < 0.8$, and the diffusive regime for $0.8 < z < 1$. Comparison with $l_O$ in panel (f) shows that this sweep of regimes is due primarily to the three orders of magnitude decrease in $l_O$, which in turn is due primarily to the fact that stratification in our case is driven by progressive absorption of radiation from above and hence is strongest close to the upper boundary. Dissipation decreases with height but this decrease is relatively small for $z > 0.15$.

After the application of surface cooling, the profiles of ${\textit {Re}}_b$ show three distinct regions: an upper region in which $\textit {Re}_b$ becomes very large and is then undefined, a central region in which ${\textit {Re}}_b$ is approximately constant with height and a lower region in which $Re_b$ remains close to its initial equilibrium state profile. The upper region corresponds to the unstably stratified convective region where $l_O$ and ${\textit {Re}}_b$ are undefined. The central region is stably stratified in a horizontally averaged sense, but convective plumes penetrate into this region providing extra mixing and leading to relatively uniform stratification conditions across this region. The lower region represents the region that is only weakly affected by surface convection.

Interestingly, at $t = 5$ there is a region close to the bottom boundary ($0 < z < 0.3$) in which $l_O > 1$. Since $l_O$ has been non-dimensionalised in terms of $\tilde {h}$, this indicates that the distance a parcel of fluid must travel in order to be significantly affected by stable stratification in this region is greater than the channel height. Since this is not possible, $l_O > 1$ indicates that, while the background stratification is still stable in this region, it is so weak that flow conditions are effectively neutral. By $t = 6$ this region has become unstable and $l_O$ is undefined throughout the channel.

4.5. Large scale flow structures in the destratified flow

In this section we describe large scale flow structures that are observed to form after all of the background stratification has been removed and the flow is evolving toward a new unstably stratified equilibrium state.

The time at which the flow becomes fully destratified, in the sense that all background stratification is removed, corresponds to a fundamental regime change. At this time the difference between the maximum horizontally averaged temperature in the channel and the horizontally averaged temperature at the channel bottom reaches zero, that is $\Delta {\theta } = \overline { {\theta }}_m - \overline { {\theta }}_b = 0$, and remains at this value since $\overline { {\theta }}_m = \overline { {\theta }}_b$. As a result, the two Richardson numbers, $Ri_{\tau }$ and $Ri_*$ also become fixed equal to $0$, and hence no longer act as governing parameters. The flow is now evolving toward a new equilibrium flow state, namely, unstably stratified turbulent open-channel flow. For Case 5 this occurs at a time of $t\approx 4$ as seen in figure 4(c).

This new equilibrium state has been studied by Walker et al. (Reference Walker, Tejada-Martínez, Martinat and Grosch2014), who suggest that it is governed by the parameters ${\textit {Re}}_{\tau }$, $Pr$ and a Rayleigh number $Ra_{\tau }$, which they define as

(4.4)\begin{equation} Ra_{\tau} = \frac{\tilde{\beta} \tilde{g} \tilde{h}^2 \tilde{q}_s}{\tilde{\rho}_b \tilde{c}_p \tilde{\kappa} \tilde{u}_{\tau}^2}. \end{equation}

For Case 5, at $t = 5$, $u_{\tau } \approx 1.2$ so $Ra_{\tau } \approx 700$. This is somewhat lower than the value used by Walker et al. (Reference Walker, Tejada-Martínez, Martinat and Grosch2014) ($Ra_{\tau } \approx 5000$). Walker et al. (Reference Walker, Tejada-Martínez, Martinat and Grosch2014) give a detailed discussion of the turbulent structure of the flow field for statistically steady open-channel flow subject to surface cooling. This material will not be repeated here. Instead we focus on the formation of large scale flow structures.

For our flow, large scale structures were observed to start forming after the regime change noted above, that is, between $t=4$ and $t=5$ for Case 5. Figure 9 shows visualisations of the temperature field at these times on a horizontal section at the channel mid-depth ($z=0.5$). At $t=4$ isolated intense downwelling thermal plumes (marked D) can be seen penetrating the central plane; however, the remainder of the temperature field is relatively homogeneous. By $t=5$, large scale structures, elongated primarily in the streamwise direction, are apparent. The low temperature regions correspond to downwelling (D) of cooler fluid from the upper boundary, while the higher temperature regions correspond to upwelling (U) of warmer fluid from below. This can be seen clearly by comparison with the transverse sections in figure 5. Each structure (an upwelling/downwelling pair) has a lateral dimension of approximately half the domain width.

Figure 9.

Temperature at $t = 4$ and 5 for Case 5 on a horizontal plane at the channel mid-depth ($z=0.5$). (a) At $t = 4$, isolated cooler downwelling thermal plumes (D) are visible penetrating this layer. (b) At $t = 5$, large coherent regions of upwelling warmer fluid (U) and downwelling cooler fluid (D) are apparent.

Inspection of flow visualisations at later times showed that these structures grow and become more coherent and elongated until, by approximately $t=6$, they eventually reach a width approximately equal to the domain width, after which they were restricted from further growth by the domain. Thus it appears that, as the flow starts evolving towards its final unstably stratified equilibrium state, the domain size also starts having an effect on the lateral dimension of these large scale structures.

This was confirmed through observation of the results for Case 9, which used the wider domain of $L_y = 2 \pi h$. On the larger domain, the structures eventually reach a size approximately equal to the domain width. This can be seen in the visualisations of horizontal and transverse sections through the temperature field for Case 9 at $t=10$ shown in figures 10(a) and 10(b). The flow is clearly divided in the lateral dimension into a single pair of upwelling and downwelling zones. This is made more clear by the transverse section of the streamwise-averaged streamfunction in figure 10(c), which shows clear clockwise and anticlockwise circulations centred on the vertical boundaries between these two zones.

Figure 10.

Temperature and streamfunction at $t = 10$ for Case 9. (a) Shows temperature on a horizontal $x\hbox{--}y$ plane at $z=0.5$. (b and c) Transverse sections at $x=0$ looking upstream. Temperature range: $-30$ to 5. Streamfunction range: $-0.6$ to 0.5. U indicates upwelling flow, while D indicates downwelling flow.

Similar large scale structures were observed by Walker et al. (Reference Walker, Tejada-Martínez, Martinat and Grosch2014), who refer to them as ‘convective supercells’. They found that the supercells in their simulations have a width four times the channel depth. Given that their domain width was $L_y = 4 \pi h/3$, this is also approximately equal to their domain width.

Thus it appears that, for the steady-state unstably stratified channel flow, domain widths up to $2\pi h$ may be placing a constraint on the size of these supercells.

This effect of domain size is confined to the stages of the flow evolution that occur after the flow has fully destratified – that is, after the time period for which $\Delta {\theta } > 0$, which is the subject of the discussions presented in the previous subsections, and in the scaling analysis presented below. For the time period over which $\Delta {\theta } > 0$, the results for Case 9, which uses the larger domain, were found to be very similar to the results for Case 6, which has the identical parameter settings on the standard size domain, indicating that domain size has negligible effect on these earlier stages of the flow evolution.

5.1. Time scales

As shown above, the destratification process involves an interplay between three dominant characteristics of the flow: turbulent fluid motion generated by bottom friction, turbulent fluid motion generated by surface cooling, and damping of these motions by the stable background stratification. As such, it is useful to define characteristic time scales associated with these three dominant processes, namely a friction time scale $\tilde {t}_{\tau }$, convection time scale $\tilde {t}_*$ and buoyancy time scale $\tilde {t}_N$:

(5.1a–c)\begin{equation} \tilde{t}_{\tau}(\tilde{t}) = \frac{\tilde{h}}{\tilde{u}_{\tau}(\tilde{t})},\quad \tilde{t}_* = \frac{\tilde{h}}{\tilde{w}_*}, \quad \tilde{t}_N(\tilde{t}) = \left(\frac{\tilde{\beta} \tilde{g} \Delta\tilde{\theta}(\tilde{t})}{\tilde{h}}\right)^{-1/2}. \end{equation}

Written in terms of our non-dimensionalisation scheme (where $t = \tilde {t}/\tilde {t}_{\tau ,0}$ with $\tilde {t}_{\tau ,0} = \tilde {h}/\tilde {u}_{\tau ,0}$), the equivalent dimensionless time scales are

(5.2a–c)\begin{equation} {t}_{\tau}({t}) = \frac{{h}}{{u}_{\tau}({t})},\quad {t}_* = \frac{{h}}{{w}_*}, \quad {t}_N({t}) = \left(\frac{\gamma \Delta{\theta}({t})}{{h}}\right)^{-1/2}. \end{equation}

In our flow the friction time scale $ {t}_{\tau }$ varies with time due to the time variation of $u_{\tau }$, while the buoyancy time scale $ {t}_N$ varies due to the change in the mean temperature gradient. On the other hand, the surface cooling flux remains constant, so the convection time scale $ {t}_*$ is constant.

The two Richardson numbers defined in (2.19) can be reformulated as ratios of these three time scales,

(5.3a,b)\begin{equation} Ri_{\tau}(t) = \left(\frac{\tilde{t}_{\tau}}{\tilde{t}_N}\right)^2 = \left(\frac{{t}_{\tau}}{{t}_N}\right)^2\quad \text{and} \quad Ri_*(t) = \left(\frac{\tilde{t}_*}{\tilde{t}_N}\right)^2 = \left(\frac{{t}_*}{{t}_N}\right)^2, \end{equation}

and hence can be interpreted as representing ratios between the time scales of shear- and convection-generated fluid motions to the time scale associated with the deceleration of fluid parcels by the background stratification.

The convection velocity $ {w}_*$ and time scale $ {t}_*$ given above use the channel height $h$ as their length scale. This approach has the advantage that it allows the development of the relatively simple scaling formulas presented below, and, in particular the formula for destratification time presented in § 5.3. The assumption implicit here is that the effect of the plumes is felt throughout the depth of the channel. This assumption is reasonable during later stages of the destratification process, but over-estimates $ {w}_*$ during the initial stages when plumes are confined to the region close to the surface.

An alternative approach is to use a length scale based on some measure of the penetration depth of the plumes in a manner analogous to the original scaling of Deardorff (Reference Deardorff1970). In this case, $ {w}_*$ and $ {t}_*$ vary with time as the plumes penetrate progressively further into the channel. An approach similar to this was adopted by Ulloa et al. (Reference Ulloa, Winters, Wüest and Bouffard2019) in their recent analysis of penetrative convection in ice-covered lakes. The resultant scaling is, however, significantly more complex because it requires knowledge or prediction of the time evolution of the penetration depth of the plumes. The development of an alternative scaling based on this more complex approach is the subject of ongoing work by the authors.

5.2. Destratification rate

In order to determine a scaling for the destratification rate, we first define a destratification rate, ${\mathcal {D}}_{\tau }$, non-dimensionalised in terms of the friction time scale and the temperature difference at that time:

(5.4)\begin{equation} {\mathcal{D}}_{\tau} = -\frac{\textrm{d}\Delta \tilde{\theta}}{\textrm{d} \tilde{t}}\frac{\tilde{t}_{\tau}}{\Delta \tilde{\theta}} = -\frac{\textrm{d}\Delta {\theta}}{\textrm{d} {t}}\frac{{t}_{\tau}}{\Delta {\theta}}. \end{equation}

This is equivalent to the destratification rate ${\mathcal {D}}$ defined in KWAZ for the case without surface cooling. In the current paper we will non-dimensionalise the destratification rate in terms of the three time scales $ {t}_{\tau }$, $ {t}_*$ and $ {t}_N$ defined above and distinguish between these by using subscripts $_{\tau }$, $_*$ and $_N$ appended to ${\mathcal {D}}$. We will also divide destratification rate into components attributed to the two main forcings – bottom shear and surface cooling. These will be indicated by superscripts $^s$ and $^c$, respectively. No superscript indicates total destratification rate. These symbols are defined explicitly in the table of symbols contained in table 3. Figure 11 shows ${\mathcal {D}}_{\tau }$ plotted as a function of $Ri_{\tau }$ for Cases 1–4, for which the surface cooling flux $q_s$ varies from 0 to 2.67 and $ {w}_*$ from 0 to 1.75. $Ri_{\tau }$ decreases as the flow evolves, so flow evolution proceeds from left to right. As $Ri_{\tau }$ decreases, the stability of the flow decreases, and the destratification rate increases. As expected, destratification rate also increases with increasing $q_s$.

Figure 11.

Value of ${\mathcal {D}}_{\tau }$ as a function of $Ri_{\tau }$ for Cases 1–4. The solid line is the model for destratification rate attributed to shear, ${\mathcal {D}}_{\tau }^{s,{\textrm {mod}}} = C_1 Ri_{\tau }^{-1/2}$, with $C_1 = 2.1$. The dashed line indicates the lower limit of validity of this model, $Ri_{\tau } = 15$.

In KWAZ we found that, for the flow without surface cooling, the destratification rate can be modelled using the power law relationship

(5.5)\begin{equation} {\mathcal{D}}_{\tau}^{s,{\textrm{mod}}} = C_1 Ri_{\tau}^{-1/2} \quad \text{for}\ Ri_{\tau} > 15, \end{equation}

where $C_1 = 2.1$. In Case 1 (no surface cooling) destratification is driven only by turbulence generated by shear, so this case isolates the effect of bottom shear on destratification rate and ${\mathcal {D}}_{\tau } = {\mathcal {D}}_{\tau }^s$.

The relationship (5.5) is shown in figure 11. For Case 1 ($q_s = 0$) the data follow (5.5) for $Ri_{\tau } > 15$. For weak stratification, where $Ri_{\tau } < 15$, the exponent decreases and ${\mathcal {D}}_{\tau }^s$ approaches an asymptotic value of ${\mathcal {D}}_{\tau }^s \approx 1.1$ for neutral conditions.

In riverine flows in which thermal stratification is an issue, the friction Richardson number of the thermally stratified state is typically of order $Ri_{\tau } = \textit {O}(1000)$ (see field measurement data in Bormans & Webster Reference Bormans and Webster1998; Sherman et al. Reference Sherman, Webster, Jones and Oliver1998; Mitrovic et al. Reference Mitrovic, Oliver, Rees, Bowling and Buckney2003). In this context the weakly stratified flow regime, $Ri_{\tau } < 15$, corresponds to a degree of stratification that is less than $\textit {O}(1\,\%)$ of its initial value and hence is negligible when determining quantities such as the time taken for effective destratification of the flow. In the following analysis we focus on the flow regime in which $Ri_{\tau } > 15$ and (5.5) is valid.

When $q_s > 0$, surface cooling provides a second driver for destratification. In order to isolate this second process we subtract the destratification due to bottom shear modelled by (5.5) from the total destratification rate:

(5.6)\begin{equation} {\mathcal{D}}_{\tau}^c = {\mathcal{D}}_{\tau} - {\mathcal{D}}_{\tau}^{s,{\textrm{mod}}} = -\frac{\textrm{d}\Delta {\theta}}{\textrm{d} {t}}\frac{{t}_{\tau}}{\Delta {\theta}} - C_1 Ri_{\tau}^{-1/2}. \end{equation}

Here, ${\mathcal {D}}_{\tau }^c$ denotes destratification rate due to surface cooling non-dimensionalised in terms of the friction time scale $ {t}_{\tau }$. Implied here is an assumption that coupling between processes associated with destratification by shear and convection is linear.

Since destratification due to surface cooling is expected to scale with the convection time scale $ {t}_*$ rather than the friction time scale, we re-scale ${\mathcal {D}}_{\tau }^c$ in terms of $ {t}_*$, to give the destratification rate due to cooling non-dimensionalised in terms of convection time scale:

(5.7)\begin{equation} {\mathcal{D}}_*^c = ({\mathcal{D}}_{\tau} - {\mathcal{D}}_{\tau}^{s,{\textrm{mod}}}) \frac{{t}_*}{{t}_{\tau}} = -\frac{\textrm{d}\Delta {\theta}}{\textrm{d} {t}}\frac{{t}_*}{\Delta {\theta}} - C_1\frac{{t}_*}{{t}_{\tau}} Ri_{\tau}^{-1/2}. \end{equation}

This can be written more concisely as

(5.8)\begin{equation} {\mathcal{D}}_*^c = {\mathcal{D}}_* - {\mathcal{D}}_*^{s,{\textrm{mod}}}, \end{equation}

where ${\mathcal {D}}_*$ is the total destratification rate non-dimensionalised in terms of the convection time scale,

(5.9)\begin{equation} {\mathcal{D}}_* = -\frac{\textrm{d}\Delta {\theta}}{\textrm{d} {t}}\frac{{t}_*}{\Delta {\theta}}, \end{equation}

and ${\mathcal {D}}_*^{s,{\textrm {mod}}}$ is the modelled destratification rate due to bottom shear rescaled in terms of the convection time scale,

(5.10)\begin{equation} {\mathcal{D}}_*^{s,{\textrm{mod}}} = C_1\frac{{t}_*}{{t}_{\tau}} Ri_{\tau}^{-1/2}. \end{equation}

By analogy with ${\mathcal {D}}_{\tau } = f(Ri_{\tau })$ and the scaling relation (1.1) of Deardorff et al. (Reference Deardorff, Willis and Stockton1980), we expect the destratification rate associated with surface cooling ${\mathcal {D}}_*^c$ to be a function of the convection Richardson number $Ri_*$. Figure 12(a) shows ${\mathcal {D}}_*^c$ plotted as a function of $Ri_*$ for Cases 2–4, for which $q_s$ increases from 0.53 to 2.67 and $ {w}_*$ from 1.02 to 1.75. The data for different surface heat fluxes collapse to a relationship

(5.11)\begin{equation} {\mathcal{D}}_*^{c,{\textrm{mod}}} = C_2Ri_*^{-1/2}, \end{equation}

where $C_2 = 1.2$. There is some scatter in the data especially at high $Ri_*$. Inspection of the data for individual cases shows an initial oscillatory variation in ${\mathcal {D}}_*^c$ around (5.11). This occurs during the early relaxation period as the flow responds to the sudden change in conditions from radiative heating to surface cooling and the region close to the upper boundary undergoes a transition from laminar to turbulent flow. This can be clearly seen in figure 12(b) which shows the data for Case 3 only, and is also apparent in the data for the total destratification rate shown in figure 11. This deviation is also due to the fact that, in our scaling, we have used the channel height as the length scale for all processes, including convection, whereas, as discussed above, early in the process the influence of thermal plumes on turbulence within the channel extends to a smaller depth. We would expect a smoother transition to occur in environmental analogues of this flow such as riverine applications, where there are typically additional physical processes such as day-time surface cooling due to evaporation and long-wave radiative emission and a gradual decrease in solar heating (see Bormans et al. Reference Bormans, Maier, Burch and Baker1997, for example) which act to destabilise the region close to the upper surface before the night-time cooling period begins.

Figure 12.

Value of ${\mathcal {D}}_*^c$ as a function of $Ri_*$. The solid line is ${\mathcal {D}}_*^{c,{\textrm {mod}}} = C_2Ri_*^{-1/2}$, with $C_2 = 1.2$. (a) Cases 1–3 for which $q_s$ increases from 0.53 to 2.67 and $ {w}_*$ from 1.02 to 1.75. (b) Case 3 only, highlighting the response during the initial relaxation period.

In order to combine these two components of destratification into a single model we combine (5.11) with (5.7) to give

(5.12)\begin{equation} {\mathcal{D}}_* = -\frac{\textrm{d}\Delta {\theta}}{\textrm{d} {t}}\frac{{t}_*}{\Delta {\theta}} = {\mathcal{D}}_*^{s,{\textrm{mod}}} + {\mathcal{D}}_*^{c,{\textrm{mod}}} = C_1\frac{{t}_*}{{t}_{\tau}} Ri_{\tau}^{-1/2} + C_2Ri_*^{-1/2}, \end{equation}

and then rescale in terms of the time scale that is common to both processes, namely the buoyancy time scale, $ {t}_N$, to give

(5.13)\begin{equation} {\mathcal{D}}_N = -\frac{\textrm{d}\Delta {\theta}}{\textrm{d} {t}}\frac{{t}_N}{\Delta {\theta}} = C_1\frac{{t}_N}{{t}_{\tau}} Ri_{\tau}^{-1/2} + C_2\frac{{t}_N}{{t}_*}Ri_*^{-1/2}. \end{equation}

Rewriting $Ri_{\tau }$ and $Ri_*$ in terms of time scale ratios given in (4.3) leads to

(5.14)\begin{equation} {\mathcal{D}}_N = C_1 \frac{{t}_N}{{t}_{\tau}}\,\frac{{t}_N}{{t}_{\tau}} + C_2\frac{{t}_N}{{t}_*}\,\frac{{t}_N}{{t}_*}, \end{equation}

and, finally,

(5.15)\begin{equation} {\mathcal{D}}_N = C_1Ri_{\tau}^{-1} + C_2Ri_*^{-1}. \end{equation}

Thus, total destratification rate non-dimensionalised in terms of the buoyancy scale $ {t}_N$ can be modelled as a simple linear combination of power-law functions of $Ri_{\tau }$ and $Ri_*$ with exponent $n = -1$. This model is valid for $Ri_{\tau } > 15$.

Figure 13 shows measured values of ${\mathcal {D}}_N$ determined from the DNS data plotted against the model (5.15) for various subsets of simulation cases. Figure 13(a) shows data for Cases 1–4 for which $q_s$ varies from 0 to 2.67 and $ {w}_*$ from 0 to 1.75. These data span a large range of $Ri_{\tau }$ and $Ri_*$ and collapse convincingly to the model relationship. Figure 13(b) shows data for cases in which the convection velocity is kept constant at $ {w}_* = 1.39$ and the parameters $Ri_{\tau ,0}$, ${\textit {Re}}_{\tau }$, $\lambda$ and $\alpha$ vary. Again, the data are well predicted by the model. Figure 13(c) shows data for the three cases using the same flow parameters on different grids and domains. The data for all three cases are very similar, indicating that the standard grid resolution and domain size used are sufficient to capture the scales of motion that contribute to the destratification process.

Figure 13.

Value of ${\mathcal {D}}_N$ plotted against the model (5.15). The solid line is ${\mathcal {D}}_N = C_1Ri_{\tau }^{-1} + C_2Ri_*^{-1}$, with $C_1 = 2.1$ and $C_2 = 1.2$. (a) Data for Cases 1–4 for which $q_s$ varies from 0 to 2.67 and $ {w}_*$ from 0 to 1.75, while other parameters are constant. (b) Data for cases in which $ {w}_* = 1.39$ while $Ri_{\tau ,0}$, ${\textit {Re}}_{\tau ,0}$, $\lambda$ and $\alpha$ vary. (c) Data for cases using the same flow parameters on different grids and domains. (d) Data for all cases.

Finally, figure 13(d) shows the data for all DNS cases. The collapse of the data across all cases provides strong evidence for the validity of the model (5.15) over the range of parameters included in this study. The main deviations seen are due to the response of the system to the sudden initial change in conditions as discussed above.

5.3. Destratification time

As a further step, the model for destratification rate (5.15) can be integrated to give a model for destratification time. Expanding all the terms in (5.15) gives

(5.16)\begin{equation} -\frac{1}{\Delta{\theta}}\frac{\textrm{d}\Delta{\theta}}{\textrm{d} {t}}\left(\frac{\gamma \Delta{\theta}}{{h}}\right)^{-1/2} = C_1 \left(\frac{\gamma \Delta{\theta} {h}}{{u}_{\tau}^2}\right)^{-1} + C_2 \left(\frac{\gamma \Delta{\theta} {h}}{{w}_*^2}\right)^{-1}, \end{equation}

which simplifies to

(5.17)\begin{equation} -{\Delta{\theta}}^{-1/2}\frac{\textrm{d}\Delta{\theta}}{\textrm{d} {t}} = \gamma^{-1/2} {h}^{-3/2}(C_1 {u}_{\tau}^2 + C_2 {w}_*^2). \end{equation}

Separating and integrating from an initial state $i$ to final state $f$,

(5.18)\begin{equation} -\int_{\Delta{\theta}_i}^{\Delta{\theta}_f}{\Delta{\theta}}^{-1/2}{\textrm{d}\Delta{\theta}} = \gamma^{-1/2} {h}^{-3/2} \int_{t_i}^{t_f} \left(C_1 {u}_{\tau}^2 + C_2 {w}_*^2\right) \textrm{d}t, \end{equation}

leads to

(5.19)\begin{equation} {t}_d = \frac{2 \gamma^{1/2} {h}^{3/2} \left({\Delta{\theta}}_i^{1/2} - {\Delta{\theta}}_f^{1/2}\right)}{C_1 [{u}_{\tau}^2]_{av} + C_2 [{w}_*^2]_{av}}, \end{equation}

where $[.]_{av}$ indicates the average taken over the time period $t_i - t_f$ (for example $[u_{\tau }^2]_{av} = 1/(t_f - t_i)\int _{t_i}^{t_f} u_{\tau }^2\,\textrm {d}t$). This can be recast in terms of dimensional variables as

(5.20)\begin{equation} \tilde{t}_d = \frac{2 \left(\tilde{\beta} \tilde{g}\right)^{1/2} \tilde{h}^{3/2} \left({\Delta\tilde{\theta}}_i^{1/2} - {\Delta\tilde{\theta}}_f^{1/2}\right)}{C_1 [\tilde{u}_{\tau}^2]_{av} + C_2 [\tilde{w}_*^2]_{av}}. \end{equation}

In our simulations, $ {w}_*$ is constant; however, $u_{\tau }$ varies during the destratification process. The values of $[ {u}_{\tau }^2]_{av}$ for each of the simulations are shown in table 1. Figure 14(a) shows a comparison of the temperature difference $\Delta {\theta }$ plotted against time $t$ for all cases, while figure 14(b) shows $\Delta {\theta }$ scaled in terms of the difference between the initial and final values $\Delta {\theta }_i-\Delta {\theta }_f$, and $t$ scaled by the destratification time $t_d$ determined from (5.19). Here, we set $\Delta {\theta }_f = 0$ so that $t_d$ represents time taken to destratify completely. (The formula can also be used to determine the time required to reach a certain degree of destratification relative to the initial stratification.)

Figure 14.

Temperature difference $\Delta {\theta }$ plotted against time $t$. Panel (a) shows $\Delta {\theta }$ and $t$ as measured in the simulations. Panel (b) shows $\Delta {\theta }$ scaled in terms of the difference between the initial and final values $\Delta {\theta }_i-\Delta {\theta }_f$, and $t$ scaled by the destratification time $t_d$ determined using (5.19). Here we have used $\Delta {\theta }_f = 0$.

Scaling $t$ in this manner gives a good collapse across the various cases with respect to the time at which $\Delta {\theta }$ reaches a value of zero, which is what $t_d$ predicts. In this regard, all cases reach zero at approximately $t/t_d = 1$, indicating that the value of $t_d$ given by (5.19) gives an accurate prediction of the destratification time. There is, however, variation in the trajectories taken in different cases. This is due to the fact that the time evolution of transfers between kinetic and potential energy within the channel varies across the different cases due to differences in the initial stratification and the rate of surface cooling.

5.4. Application to predicting destratification in rivers

In principle our equations for destratification rate (5.15) and destratification time (5.20) could form the basis for methods to predict destratification in rivers. (5.15) could be incorporated into a one-dimensional river model and integrated in time alongside other equations in order to predict changes in stratification over time. Alternatively, (5.20) could be used directly to predict destratification time from a given initial state for a given surface cooling flux. In either approach $\tilde {w}_*$ can be estimated based on an estimate of heat fluxes determined from predicted or recorded weather conditions, while $\tilde {u}_{\tau }$ can be estimated from river flow rate using formulas such as the Manning or Chézy equations (Bjerklie, Dingman & Bolster Reference Bjerklie, Dingman and Bolster2005). These formulas do not account for the effect of stratification on the friction coefficient $C_f$ and hence $\tilde {u}_{\tau }$. In order to include this effect, a correlation for $C_f$ at values of ${\textit {Re}}_{\tau }$ and $\lambda$ typical of stratification conditions in rivers of interest (${\textit {Re}}_{\tau } \approx 25\,000$ and $\lambda \approx 5$) would be required.

The scaling relations above are based on simulations at values of ${\textit {Re}}_{\tau }$ and $Pr$ that are significantly lower than those seen in real rivers. As discussed in § 3, results from KWAZ and other authors indicate that destratification rate is likely to show some dependence on $Pr$. Given that our relations have been developed based on physical arguments, the hope here is that, as with similar scaling relations such as (1.1) of Deardorff et al. (Reference Deardorff, Willis and Stockton1980), the dependence of our scaling relationships on $Pr$ may also be restricted to the magnitude of the coefficients. While a dependence on Reynolds number is not apparent in our results, only a small range of ${\textit {Re}}_{\tau }$ has been tested here. Given that ${\textit {Re}}_{\tau }$ of typical stratified river flows is significantly higher, even a small ${\textit {Re}}_{\tau }$ dependence would have an effect. Thus, before using the relationships proposed here for predicting destratification rates and times in rivers, it is necessary to test these relationships through comparison with field or experimental measurements.