Stratification and temporal evolution of mixing regimes in diurnally heated river flows

Direct numerical simulations of stratified open channel flows subject to a varying surface heat flux are performed. The influence of the diurnal heating time on the spatial and temporal variation of mixing in the flow and the characteristics of the mean flow state are examined. The control parameters are the bulk stability parameter λB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda_{B}$$\end{document}, defined through the ratio of the channel height δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta$$\end{document} and a bulk Obukhov length scale LB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathscr{L}_{B}$$\end{document}, and the diurnal time scale t^\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t}$$\end{document}, defined as the ratio of the heating time to an eddy turnover time. The Prandtl number Pr and Reynolds number Reτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Re_{\tau}$$\end{document} have values of 1 and 400. Simulations are performed over t^=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} = 1$$\end{document} to 24 and λB=0.6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda_{B} = 0.6$$\end{document} to 26. Two key flow features are used to classify the flow regimes observed, namely the laminar layer depth (LLD) and stratified layer depth (SLD) where the LLD is defined as the depth from the free surface when the buoyancy Reynolds number ReB≈7\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Re_{B} \approx 7$$\end{document} and the SLD is the depth from the free surface when the turbulent Froude number Fr≈1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Fr \approx 1$$\end{document}. This study attempts to characterise how these length scales vary across the diel cycle. The LLD is a viscous length scale and a regime map of a viscous parameter, the bulk Obhukov Reynolds number ReL\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Re_\mathscr{L}$$\end{document}, and t^\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t}$$\end{document} is presented to classify the LLD behaviour. A regime map of λB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{B}$$\end{document} and t^\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t}$$\end{document} is presented to classify the behaviour of the SLD. Three classifications for each layer depth behaviour within a diel cycle form the basis of the regime maps for this paper: a neutral flow where the LLD or SLD does not exist (denoted by NL and NS), a stratified flow where the LLD or SLD are diurnally varying (denoted as DL and DS) and a persistent layer of the LLD or SLD (denoted as PL and PS). The transition between the NL to DL is t^∝ReL4.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto Re_{\mathscr{L}}^{4.5}$$\end{document}, DL to PL is t^∝ReL-0.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto Re_{\mathscr{L}}^{- 0.5}$$\end{document}, NS to DS is t^∝λB0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto \lambda _{B}^{0}$$\end{document} and DS to PS is t^∝λB1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto \lambda _{B}^{1}$$\end{document}. The regime maps may be used as a predictive tool to determine when suppressed mixing regimes occur in rivers. At each flow depth, the flow sweeps though a range of mixing states across the diel cycle. The local mixing efficiency are briefly assessed and found to scale well with the instantaneous Fr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Fr$$\end{document} number according to the regimes proposed by Garanaik and Venayagamoorthy (J. Fluid Mech., vol. 867, 2019, pp. 323-333). This paper reports on direct numerical simulations of stratified open channel flows subject to a varying surface heat flux. The results have found that: Increasing the diurnal time scale allows the flow to sweep through a wider range of flow states from turbulent to strongly stratified, Simulation data of the mixing efficiency and turbulent Froude number from this temporally varying and spatially inhomogenous flow that undergoes strong temporal forcing collapses well onto the parameterisation scheme of Garanaik & Venayagamoorthy (J. Fluid Mech., vol. 867, 2019, pp. 323–333) found for homogenous stratified flows, and Three distinct classifications (a persistent layer, a diurnal layer and one where the layer does not exist) of the laminar layer depth (LLD) and stratified layer depth (SLD) behaviour persists throughout a diel cycle and form a regime map given a λB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{B}$$\end{document}, Reτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Re_{\tau }$$\end{document} and t^\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t}$$\end{document} value. The relationship between each transitions are: from no LLD (NL) to a diurnal LLD (DL) t^∝ReL4.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto Re_{\mathscr {L}}^{4.5}$$\end{document}, DL to a persistent LLD (PL) t^∝ReL-0.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto Re_{\mathscr {L}}^{-0.5}$$\end{document}, no SLD (NS) to a diurnal SLD (DS) t^∝λB0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto \lambda _{B}^{0}$$\end{document} and DS to a persistent SLD (PS) is t^∝λB1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto \lambda _{B}^{1}$$\end{document}. The regime maps may used as a predictive model to calculate when suppressed mixing transpires in rivers. Increasing the diurnal time scale allows the flow to sweep through a wider range of flow states from turbulent to strongly stratified, Simulation data of the mixing efficiency and turbulent Froude number from this temporally varying and spatially inhomogenous flow that undergoes strong temporal forcing collapses well onto the parameterisation scheme of Garanaik & Venayagamoorthy (J. Fluid Mech., vol. 867, 2019, pp. 323–333) found for homogenous stratified flows, and Three distinct classifications (a persistent layer, a diurnal layer and one where the layer does not exist) of the laminar layer depth (LLD) and stratified layer depth (SLD) behaviour persists throughout a diel cycle and form a regime map given a λB\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{B}$$\end{document}, Reτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Re_{\tau }$$\end{document} and t^\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t}$$\end{document} value. The relationship between each transitions are: from no LLD (NL) to a diurnal LLD (DL) t^∝ReL4.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto Re_{\mathscr {L}}^{4.5}$$\end{document}, DL to a persistent LLD (PL) t^∝ReL-0.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto Re_{\mathscr {L}}^{-0.5}$$\end{document}, no SLD (NS) to a diurnal SLD (DS) t^∝λB0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto \lambda _{B}^{0}$$\end{document} and DS to a persistent SLD (PS) is t^∝λB1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{t} \propto \lambda _{B}^{1}$$\end{document}. The regime maps may used as a predictive model to calculate when suppressed mixing transpires in rivers.


Introduction
Water bodies such as rivers, estuaries and canals are subject to diurnal solar radiation.This unsteady heating causes these channels to experience time-varying flow states within their diel cycle.During the day, the short-wave radiation at the free surface is transmitted and progressively absorbed through the water column, leading to a non-uniform stable thermocline.Therefore, stable stratification is stronger near the free surface and gradually weakens towards the channel bed.This stratification then acts to dampen the turbulent mixing within these regions.Accounting for the diurnal variation of solar fluxes, solar radiation tends to strengthen stratification while the absence of radiation will cause the temperature field to relax, increasing turbulent mixing [1].This mixing determines the overall ecological health of the limnological ecosystems [2].As turbulence becomes suppressed as a result of stratification due to daytime heating, in extreme cases, this reduction in mixing limits the transportation of scalars such as dissolved oxygen, sediments and nutrients, causing the benthic region of the water column to experience conditions of eutrophication or hypoxia [3][4][5].
Among Australian inland rivers, during long periods of drought and high radiative forcings, flow rates are significantly reduced, causing a decrease in the turbulence production from the shear at the bottom of the channel.These low flows and persistently stratified flow states contribute to the extensive algal outbreaks resulting in mass fish mortality as seen in the Menindee Region of NSW over the summer of 2018 and 2019 [6].Cyanobacterial blooms contribute to oxygen depletion at the near-wall or bottom region of the river channel as their substantial biological matter sinks to the bottom of the channel and consumes the oxygen reserves for decomposition.If a rapid breakdown of stratification occurs, the sudden overturning of the oxygenated epilimnion layer with the anoxic hypolimnion layer will cause the water column to dilute and therefore suffocate the wildlife that concentrates in the previously oxygen-rich region [7].Understanding the processes that affect the mixing and turbulence of diurnally heated stratified flows can help predict and mitigate these detrimental events.
Using global energy balances, Simpson and Hunter [8] developed a criterion for the onset of stratified conditions in continental shelves, that included the mechanical work done by tidal stresses and the solar heat input.The mixing criteria is therefore just the ratio of the rates of stratifying thermal energy and destratifying turbulent kinetic energy.The work was latter extended by Holloway [9] for estuaries that incorporated a windmixed system and depth dependent heating function with no tidal forces and by Bormans and Webster [1] for turbid rivers.These works determined that the degree of stratification depends on the magnitude of tidal or river discharge, wind speed, and the attenuation of the radiative heat flux through the water column, which is dependent on the water body's turbidity and channel depth .In Borman and Condie [10], their study into the stratification dynamics of rivers modelled the diurnal variation of solar heating as a depth-varying volumetric heat source following the Beer-Lambert law, Q(Y, T) = I s (T) e (Y− ) , where T is the time, Y is the vertical position, and I s is the radiant heat flux through the surface.Williamson et al. [11] characterised stratification in radiatively heated open channels using DNS.The bulk stability parameter in this case was equal to the ratio of a confinement length scale related to the domain height and the Obukhov length scale L = o C p U 3 1  2 − 1 ∕ g I s , where o is the reference den- sity, C p is the specific heat of the fluid, U is the friction velocity, g is the gravitational acceleration and is the coefficient of thermal expansion [11, 12].Flows with = 1 is found to be strongly stratified with the flow in local energetic equilibrium over much of the channel and at Re L = LU ∕ ≃ 400 , where is the kinematic viscosity, laminarisa- tion occurs at the upper layer y ≈ 0.8 [11].Because descriptions of the state of turbu- lence in stratified flows can be characterised by the stratified outer layer, the definition of the Reynolds number can be defined using the bulk stability parameter B and , resulting in the bulk parameter Re L = Re ∕ B .Here, Re is the friction Reynolds number and is defined as Stably stratified open channel flows are used to investigate many of the environmental flows previously mentioned above [13,14].For open channel flows where stratification is spatially inhomogeneous, the vertical domain can acquire a wide range of gradient Richardson numbers Ri g = N 2 ∕S 2 , Froude numbers, specifically the turbulent Froude number Fr = ∕NK , and buoyancy Reynolds numbers Re B = ∕ N 2 [15][16][17][18][19]. Here, is the tur- bulent dissipation rate, N is the buoyancy frequency, S is the mean vertical shear, and K is the turbulent kinetic energy.These parameters are often considered important parameters in governing the effects of buoyancy of stratified flows.Though Ri g has been considered impractical for quantifying mixing of moderately to strongly stratified flows [13,20], Re B and Fr remain parameters of interest in parameterising the effects of stratification on mixing in turbulent flows.For Re B , the parameter indicates the range of length scales over which motions are considered to be largely unaffected by stratification [21,22] while Fr is a ratio of the characteristic timescale of turbulence and stratification and is suggested to be a good indicator of the local state of turbulence in a stably stratified flow [16,22].In terms of local states of turbulence, Re B < 7 is thought to be within a molecular regime [15] and Fr ≪ 1 in a strongly stratified regime.Understanding where these local mixing regimes exist within the water column and how they vary with bulk parameters is often useful when aiming to predict the behaviour of the flow.
While Williamson et al. [11] and Issaev et al. [23] studied the transition to statically stably stratified states and parameterization of mixing in stratified flows with a constant heating force, and Kirkpatrick et al. [24] investigated the removal of surface radiation from a thermally stratified open channel, only a handful of studies examined a flow that undergoes diurnal variation of a radiant heat source.Lei and Patterson [25] numerically investigated the natural convection induced by diurnal heating and cooling in varying topographical reservoirs.Their studies involve heating and cooling in shallow and deep reservoirs to observe their temperature profiles and flow rates.Their findings show a distinct time lag in response to the transient thermal forcing beyond one quarter of the heating-cooling cycle and dependence on their control parameter, the Grashof number.Lei and Patterson [25] further observed thermal instabilities to be the dominant driving mechanism responsible for vertical mixing in the reservoirs.While this paper does not study varying topographical reservoirs, it does simulate the diurnal solar radiation behaviour experienced in river flows.This current paper expands on the work of Williamson et al. [11] by introducing an unsteady radiant heat flux through the surface of an open channel flow and aims to uncover the local flow states, where these states exist within the water column, and the extent of their changes across the diel cycle as the flow goes through its cyclic heating (stratifying) and adiabatic -destratifying stages.
This paper is divided into seven sections.Section 2 gives the problem formulation, describing the mathematical set-up of the flow along with the governing equations and subsequent non-dimensionalisation of important parameters.Section 3 details the numerical simulations, stating all simulation parameters tested for each case.Sections 4 to 6 contain the analysis of the flow's transient response, the layer depths and their regimes and the relationship between the flow's local flux Richardson number to the turbulent Froude number, respectively.Finally, section 7 concludes this paper.

Problem formulation
The schematic of the open channel flow, illustrated in Fig. 1, has a free-slip, adiabatic upper surface and a no-slip, adiabatic lower wall.The model is periodic in both the horizontal planes, and the stream-wise direction is driven by a constant uniform pressure gradient.The water column experiences a progressive absorption of the surface solar radiation that gives rise to the thermal stratified state of the flow.This internal heat source is characterised by a volumetric depth-varying heat source Q(Y, T), conveyed by the Beer-Lambert law [1], where is the height of the channel, is the absorption coefficient and I s is the short-wave heat flux which varies with time.Here, the diurnality of I s is based on the piecewise function defined below [26], with I m signifying the maximum irradiance at noon and D signifying the diel cycle, a 24 h period where D = 86, 400s and therefore, can be used to determine the sunrise time T = 0.25D , length of light period T = 0.5D and sunset time, T = 0.75D.
The diurnal timescale t is non-dimensionalised by the channel depth, diel cycle, and friction velocity: with U defined as the friction velocity related to the shear stress at the lower surface of the channel and as the channel height.The dimensionless new day, sunrise time, length of the light period and sunset time can therefore be defined to be; The temperature field is decomposed as: where � (X, T) is the statically fluctuating factor and φ(T) is the domain-average tempera- ture at time T.
Non-dimensionalising the depth-dependent heat forcing q h and the temperature fluctuating field with the normalising components Q b and b for one diurnal period, respectively, can be shown as: where, and with Q(T) being the domain-average radiant heat flux for one diurnal cycle, ̃ (T) being the domain-average temperature field, C p is the specific heat of the fluid, and o is the reference density.Throughout this paper, capitalisation of a variable signifies a dimensional quantity, whereas a lower-case variable indicates non-dimensionality.Here, Q b is described as, where with I s as the average radiant heat flux through the surface for one diel cycle.Since strati- fied layer depths are predominately above y = 0.6 , the normalisation of q h and , along with this paper's modified definition of the bulk stability parameter B , is defined over the distance y = 0.6 to .As previous work of Williamson et al. [11], Issaev et al. [23] and Kirkpatrick et al. [24] have shown that local conditions are well aligned with mean flux profiles, the definition of the normalising components and the bulk stability parameter in this paper is more directly connected with local dynamics than for one based on the entire channel height as used in Williamson et al. [11], Issaev et al. [23] and Kirkpatrick et al. [24].
As mentioned above, stratification of the channel is characterised by the bulk stability parameter B , which is the ratio of the confinement scale to a modified bulk Obukhov length scale L B .These terms are defined as: where g is the gravitational acceleration and is the coefficient of thermal expansion which relates the fluid density to the temperature by d ∕ o = − d .b is defined as: where s is the surface shear from wind and w is the wall shear.
Direct numerical simulation (DNS) is applied to solve the Navier-Stokes equation for the stratified open channel flow.Using the Oberbeck-Boussinesq approximation for buoyancy [11], the governing equation for the non-dimensional conservation of mass, momentum and energy can be written as: and, with u i being the Cartesian components of the velocity vector u = (u, v, w) , x i the compo- nents of the position vector x = (x, y, z) , p is the pressure and, e x and e y are the unit vectors in the x and y directions.The non-dimensionalised form of the components above can be described as: with X, U, T, P, U and as the above equations dimensional counterpart.The friction Reynolds number is Re = 400 and Prandtl number Pr = ∕ = 1, where is the scalar diffusivity.The problem can therefore be fully defined by specifying Re , Pr, B , and the non-dimensional turbidity parameter found in Eq. 1.This paper only considers cases where = 8 so that the solar radiation does not fully penetrate the channel bed and stratification is confined only to the upper layer of the domain.
Referring back to the boundary conditions; the adiabatic no-slip bottom surface ( y = 0 ) and the adiabatic, stress-free top surface ( y = 1 ) can now be described as:

Numerical simulation
Simulations were solved using the three dimensional, Cartesian structured and fractionalstep finite-volume method described in Armfield et al. [27].Cell face velocities were calculated using a Rhie-Chow interpolation with the spatial discretisation using a fourth-order central differencing for advective approximation and a second-order central differencing for diffusion approximation for both velocity and scalar terms.A second-order accurate in time Adams-Bashforth scheme was employed for nonlinear terms, while a Crank-Nicolson scheme was used to advance the diffusive terms.A stabilised bi-conjugate gradient solver was used to solve the pressure correction equation, while a Jacobi solver was employed to solve the momentum and temperature equations.A Courant number between 0.16 and 0.17 was used to adjust the time step, ensuring that the Courant number stays within the set range.
Table 1 shows the simulation parameters and grid and domain sizes that were tested.The horizontal axes, x and z, are uniform in grid size with cell sizes in viscous wall units being △x + = 5 and △z + = 2.5 .The vertical axis is non-uniform, stretched, and set on a log mesh that grows from both ends.In other words, between y = 0 − 0.5 , the cell size in viscous wall units range from △y + = 0.5 − 3.3 and symmetrical over the rest of the domain height.Grid sizes 512 × 165 × 512 of the constant forcing case for a domain of 2 × 1 × were employed and compared to Williamson et al. [11] to which these results yielded indiscernible differences allowing this grid resolution to be used.Comparisons between a half-span domain of 2 × 1 × 0.5 and a full-span domain of 2 × 1 × indicate negligible dif- ferences in results when run with a constant heat force and identical parameters for each simulation.The vertical profiles of Fig. 2 are averaged over horizontal planes and time with perturbation from the mean denoted by a prime.Here, ⟨⋅⟩ indicate averaging over the horizontal planes and ⋅ symbolise averaging in time.Figure 2 are profiles of the mean tem- perature, velocity, scalar and turbulent shear fluxes, buoyancy Reynolds number Re B and turbulent Froude number Fr.The similarity in results shown in Fig. 2 for both flux and important bulk parameters such as Re B and Fr demonstrate that the half-span domain can be applied, and this in effect helps in lowering computational power and run time for each simulation.Simulations were initialised from a realisation of a stably stratified constant forcing flow.Preliminary tests showed the flow to be independent of the initial conditions.Once quasisteady-state is achieved, statistics are comparable to simulations with identical parameters and grid and domain size regardless of the initialisation states.where h is the simulation domain height, ⟨Δ ⟩ is the horizontally averaged temperature at the top and bottom of the channel, u b is the bulk velocity and, where K is turbulent kinetic energy and ⋅ indicate averaging over the entire domain.The turbulent dissipation rate and buoyancy frequency N are represented as follows: with S ij is the strain rate due to velocity fluctuations given by S ij = 1∕2( U � i ∕ X j + U � i ∕ X j ) and U ′ indicates the fluctuation stream-wise velocity.
Figure 3 shows the results of the initial condition independence tests.There is an observable difference between the start-up non-dimensional time t = 0 ; however, as time progresses, this difference decreases until the two simulation results collapse onto one another.Convergence here is defined when the average parameter of each diel cycle differs by less than five percent from the previous cycle.Once achieved, the oscillatory flow is fully developed, and quasi-steady-state is established.Therefore, for this study, quasisteady-state is defined when the cycle averages, as well as the maximum and minimum amplitudes of essential parameters such as the bulk Richardson number Ri B (Fig. 3a), bulk Froude number Fr B (Fig. 3b), and the friction Richardson number Ri = ⟨Δ ⟩ B h∕u 2 (Fig. 5) differs less than five percent from the previous cycle.Under these conditions, the time taken for the flow to converge and reach a quasi-steady-state is t ≈ 12 .Figure 4 illus- trates the horizontally averaged and time-averaged at each new day, rise time, mid-day, and set time for the scalar flux profile of Case INF and Case ISF.It is evident in all figures that the profiles are similar when statistics are taken after t ≈ 12 from the initialisation of the flow.
Provided that the mean Re and B of the initialised flow field are similar to the simula- tion conditions, the flow will reach full development when t ≈ 12 .That is, regardless of the ( 20) diurnal timescale t value of the simulation, when the flow is initialised with similar B and Re values to the current simulation, the development time is insensitive to t .This is shown in Fig. 5, a plot of the friction Richardson number Ri time-series for a constant forcing case and varying t .For all cases in Fig. 5, initial conditions are started from a flow with parameters B = 5.8 , Re = 400 , Pr = 1 and = 8 .As Fig. 5b, c and d have an unsteady thermal force applied to the surface of the open channel, the flow cycles through different flow states in space and time.Ri captures the buoyancy effects on the channel and for each case aside from Fig. 5a, buoyancy is at its peak between t = 0.5t and t = 0.75t while the lowest Ri lies between t = 0.25t and t = 0.5t .The flow throughout the diurnal cycle will experience a repeating cyclic trend where the flow becomes progressively stratified as the levels of thermal forcing increase and decrease with the reduction or absence of the thermal forcing.
Figure 5 shows a fully developed, quasi-steady-state flow that is established from t = 0 for different t values but initialised from the same flow field.This fully developed, quasi- steady-state flow developed right at t = 0 suggests that the flow fields from one t result can be used to initialise simulations at other t values efficiently if B and Re are kept constant.Shown in Section 6, the diurnal average flow quantities such as Fr are relatively insensitive to t which supports this conclusion.For the remainder of this paper's simulations, flow fields were initialised with parameters B = 5.8 , Re = 400 , Pr = 1 and = 8 for cases 2 to 5 with statistics being taken after one diurnal timescale cycle.The remaining cases are initialised from the aforementioned parameters though statistics are taken after four diurnal timescale cycles.

Temperature evolution
The temperature evolution of an unsteady thermal forcing on a flow is discussed in this section.Figure 6 shows the visualisation of the transient temperature field for Case 3 with simulation parameters; B = 5.8 , t = 6 , Re = 400 , Pr = 1 and = 8 .The colour bar is scaled between 0 − 1 of Fig. 6 and normalised by ( − min ) ∕ ( max − min ) to highlight the turbulent traits for each image.Initially, as shown in Fig. 6a, the flow is unstratified in the turbulent region of the lower half of the channel and is progressively stratified in the upper half of the domain as shown through the gradual colour changes at the near surface of the channel y = 0.8 to y = 1 .As the flow evolves, turbulence is noticeably more ener- getic as the temperature field mixes though the channel and stratification is broken down as presented in Fig. 6c.Between the period, t = 0t to t = 0.25t (Fig. 6a and b), the flow experiences a zero surface heat force and turbulence begins to extend to the upper half of the channel, reducing the thickness of the near surface laminar layer and the strength of the stratification.This behaviour continues past the introduction of the heat forcing at t = 0.25t until re-stratification occurs from t = 0.5t and t = 0.75t (Fig. 6c and d), where the near surface region of the channel y = 0.8 to y = 1 settles into an almost-laminar state and exhibits distinct shear instabilities.
Figure 7 shows the vertical temperature profiles of cases 2 to 5 (refer to Table 1) to highlight the effects of varying diurnal timescales on the flow as well as comparing the diurnal flows with its constant forcing counterpart.As the vertical temperature profiles vary throughout a diel cycle for diurnal cases, temperature profiles at t = 0t , t = 0.25t , t = 0.5t and t = 0.75t are shown.The constant forcing case in Fig. 7 shows a well-mixed turbulent region near the bottom wall transitioning to a thermocline region extending from y = 0.7 − 1 .Unsurprisingly for all diurnal cases, the flow sweeps through a wider range of temperature profiles from an almost neutral flow state, to a weakly stratified flow, and to a strongly stratified flow.This range becomes greater as t increases where Fig. 7b with t = 6 experiences weakly to strongly stratified flows compared to Fig. 7d where t = 24 , which transitions through to a complete isothermal state as shown by the t = 0.25t line, to a weakly stratified flow shown by line t = 0.5t and finally to a strongly stratified flow indi- cated by the t = 0.75t line.This sweep is attributed to the length of the diurnal cycle where shorter periods do not fully allow the flow to respond completely to the external heat forcing as compared to a longer diurnal cycle.Observing the near-wall to mid-height of the channel y = 0 − 0.6 , this region is shown to remain relatively unstratified throughout the diel cycle however as t increases, the flow begins to experience changes within this area.This behaviour can be seen when comparing the stratified region of Case 3 where t = 6 (Fig. 7b) and Case 5 where t = 24 (Fig. 7d).For Case 3, as it progresses through the diur- nal cycle, each profile line in Fig. 7b collapses onto one another between the y = 0 − 0.6 region, compared to Fig. 7d, where the varying time profiles of Case 5 do not collapse in this region.
Upon comparing profiles when t remains constant and B varies, unlike temperature profiles of Fig. 7, Fig. 8 shows that temporal temperature profiles of each case decrease in differences as B increases.Figure 8 compares the temperature profiles of B = 1.5 and  8, the temperature difference at the thermo- cline rises with the thermocline extending deeper through the channel as B increases.For changes in B at a certain value where the flow can be considered strongly stratified throughout the diel cycle, the lower mixed region of the channel remains insensitive to changes in the two parameters as shown through Case 3 (Fig. 7b) and Case 7 (Fig. 8b).

Removal of diurnal heat source
Considering the period of the flow where the heat source is suspended, the destratification rate for a thermally stratified flow after the removal of the heat source can be described as [28]: with Ri described as: where Δ is the temperature difference at the top and bottom of the channel, t = h∕u , t N = ( Δ ∕h) −1∕2 and the bulk stability parameter : where, The destratification process can be thought of as the time when the initial state of a flow moves from one that is strongly affected by buoyancy to a final state in which buoyancy effects are minimal [28].The relationship between this destratification rate D s and Ri when Ri  > 15 follows the equation below [28]: Equations 22 to 26 are used to model a destratification time t d for a flow with no surface cooling [28].The destratification time t d is defined as [28]: where the subscripts i and f represent the initial and final parameter values, u 2 avg is equal to the average u 2 between the initial and final time when Δ f = 0 .As these simulations do not experience a time when Δ f = 0 , u 2 avg is taken between the time when the radiant heat flux is removed from the flow.The modified bulk stability parameter B (given in equation 21) and Kirkpatrick et al. [28] bulk stability is related by = 0.061 B for the constant forcing case and = 0.286 B for the diurnal cases.
In Table 2, the results for the destratification time defined through Eq. 27 are taken from when the flow experiences a zero radiant heat flux, between t = 0 − 0.25t and t = 0.75 − t .According to Table 2, it is expected for all diurnal cases to experience a period where Δ = 0 as the predicted t d from Eq. 27 is approximately less than or equal to the phase where the surface heat source is discontinued t = 0.5t of each case.This expected destratification however does not occur for cases 3, 4 and 18 as shown in Fig. 9, and while Case 5 does destratify during its diurnal cycle, the t sim or the time taken for Case 5 to reach Δ f = 0 from t = 0.75t (when heating is completely removed from the flow) is t ≈ 11 , a 6.4 difference from t d .These discrepancies may be due to the variation of the paper's flow compared to Kirkpatrick et al. [28] who investigated thermally stratified turbulent channel flows after the removal of a radiant heat source instead of a diurnally varying radiant heat source.As Eq. 27 is dependent on the temperature difference, having diurnal heating affects Δ as shown in Fig. 9 with the maximum Δ in the constant forcing simulation being much greater than for a flow subject to diurnal heating.Results from Table 2 show that t sim compares well to 2.1t∕t d .Applying this will give Case  5 Layer depths and flow regimes

Laminar and stratified layer depths
This section investigates the changing laminar and stratified layer depths as the flow evolves through time and is subject to a transient radiant heat source.In this analysis a laminar layer indicates minimal to no turbulent mixing within the region and certain flows, given the right conditions, can exhibit a persistent, diurnal or absent laminar layer close to the free surface throughout its diel cycle.These three behaviours are also exhibited for stratified layer depths which are depths where the variations of temperature are highest and tend to act as a thermal barrier for vertical mixing [29].
The laminar layer depth (LLD) for this paper is defined as the depth from the free surface where the buoyancy Reynolds number Re B , a buoyancy parameter used to indicate the separation of the smallest eddy affected by buoyancy to the smallest scale of turbulence, is equal to 7 [15].Within this diffusive range Re B < 7 for stably stratified shear flows, it was found that turbulence is strongly damped and minimal lateral mixing occurs [15].
The stratified layer depth (SLD) for this paper is defined as the distance from the free surface to where the turbulent Froude number Fr, defined as the ratio of buoyancy to turbulence times scales, is equal to 1 as Fr ≪ 1 lies within the transition into a strongly stratified regime [16].Since Fr is defined through local turbulent quantities, Fr can be considered a reasonable measure for the local state of turbulence in a stably stratified flow with Fr ∼ 1 shown to be a good definition for the onset of strong stratification in many studies of turbulent flow [16,30].Furthermore, the Froude number is a significant parameter in inferring the state of turbulence and parameterises well with mixing efficiency in stably stratified turbulent flows [16,[31][32][33].
Figure 10 shows the time-series of the LLD and SLD for simulation cases 2 to 5. For the constant forcing case, Fig. 10a shows an almost-constant layer depth with the average LLD = 0.05 and the SLD = 0.57 while oscillating layer depths occurs for diurnal cases, cases 3 to 5. All diurnal simulations upon observing the time-averaged LLD have an average LLD = 0.10 .The average SLD for diurnal cases vary as t increases.Case 3 ( t = 6 ) has an average SLD = 0.56 (much like the constant forcing case where SLD = 0.57 ) whereas Case 5 ( t = 24 ) has an average SLD = 0.46 .On the minimum depths of the two layers, as t increases, Fig. 10b to d shows the lowering of the minimum depth of the layers, and in some instances, these layer depths will go to zero.The maximum LLD and SLD, in all cases, increase with t.
Figure 10d reveals a rapid decline of the SLD at the end of each diurnal cycle, indicating that the channel is weakly stratified over its entire depth and is the only case that exhibits this complete breakdown before recovering to an almost constant SLD.This relationship also highlights that in extreme cases where t is very small, its behaviour, when averaged, mirrors that of its constant forcing counterpart.
When changing B and keeping t constant, the flow will progressively increase in strati- fication and its LLD and SLD will increase in size as shown in Fig. 11a where B = 1.5 will have no LLD throughout its diel cycle, while Fig. 10b where B = 5.8 will have a diurnal LLD and finally to Fig. 11b where the LLD is persistent throughout the diel .Furthermore, the length of the LLD will also grow and penetrate further down the water column as B increases.A similar behaviour is observed with the SLD.These behaviours exhibited in both the LLD and SLD are much like when t increases for the same B as shown above.
For all simulations, there is a distinct lag between the thermal forcing and the flow response in all diel simulations for both LLD and SLD shown in Figs. 10 and 11 where the LLD and SLD are slightly out of phase with the black dashed lines that represent the unsteady heat forcing.Furthermore, for all cases, the SLD is greater than LLD, with the maximum SLD ranging between 0.55 and 0.75.
The results in Fig. 10 demonstrate that for increasing t the flow transitions from persis- tently laminar in the near surface region to diurnally laminar.The SLD also demonstrates this transition although a flow which is diurnally laminar is not necessarily diurnally stratified as Fig. 10c shows.In certain cases (refer to Table 1), the flow may also attain a regime where the flow will never acquire an LLD nor an SLD (Case 8), a regime where the flow won't obtain an LLD but will have a diurnal SLD (Case 6, Fig. 11a) and a regime where LLD and SLD are persistent (Case 7, Fig. 11b).In the following section, a regime map of these behaviours for the LLD and SLD is presented where each simulation will fall into one region of the LLD regime map and another for the SLD regime map.The LLD map denotes all three flow regimes: never laminar (NL), diurnally laminar (DL) or persistently laminar (PL) while the SLD map will never stratified (NS), diurnally stratified (DS) or persistently stratified (PS).It can therefore be said that when increasing B , the flow will transition into higher stratification regimes (NL > DL > PL or NS > DS > PS).

Flow regimes
The location of the LLD and SLD depends on the governing parameters, B or Re L = L B U ∕ = Re ∕ B , and t .Simulations of varying B and t have been used to iden- tify the flow regimes defined through these governing parameters and locate these regimes on a regime map.The regime map in Fig. 12 reveals a flow's response against varying Re L for LLD and B for SLD against t.Fig. 12 Regime map for an absent, diurnal or persistent a LLD and b SLD for all cases given in Table 1.As LLD is a viscous length scale defined by Re B ≈ 7 , Re L is used as the relative vis- cous parameter to map the regime of the LLD.Re B is employed to defined the LLD as it is often used as an indicator of the collapse of turbulence or the transient relaminarisation in stratified flows.It is a parameter defined through the ratio of the characteristic length scales, the Ozmidov length scale l O to the Kolmogorov length scale in stratified turbulent flow [18,21,34].The Ozmidov length scale l O is a length scale defined as the smallest (vertical) eddy influenced by buoyancy whereas the Kolmogorov length scale characterises the smallest scales of motion [21] therefore, Re b suggests the dynamic range of scales over which motion remain unaffected by the buoyancy force that dampens the larger scales and the viscous dissipation that affects the small [21].Because Re L is a bulk scale analogous to Re B in equilibrium flows [11], it is on this basis that it is used to construct a regime map for the laminarisation the t and Re L space.
While Re L is used to characterise the regime map for LLD, B is the characterising parameter used for regime map based on the stratified layer depth.The SLD is defined as the location from the free surface where Fr ≈ 1 as a strongly stratified regime is found to be Fr ≪ O(1) [16,22,32].In the limit of strong stratification, buoyancy effects are dominant and hence it is reasonable to assume that the bulk stability parameter in this paper B is appropriate to map the SLD regime map.
The LLD and SLD may break down and re-establish throughout the diel cycle; however, at sufficiently low B or higher Re L , these layers may be completely non-existent, and the flow is fully turbulent over the entire channel depth.At very high B or low Re L , the layer depths may persist through the diel cycle as shown in Fig. 10.The first classification; a non-existent LLD or SLD throughout the diurnal cycle, is denoted by NL and NS, respectively.The second classification; a break down and re-established LLD or SLD throughout the diurnal cycle, is denoted by DL and DS, and the last classification; a persistent LLD or SLD, is denoted by PL and PS.
The transition lines between these regimes are found by identifying the region of the transition before fitting a trend-line onto the regime map.For Fig. 12a, the transition from NL to DL is given by the equation t = 2.7 × 10 −10 Re 4.5 L − 0.1 and the transition from DL to PL is dictated by t = 1.56 × 10 2 Re −0.5 L − 18 .For a constant zero stratified layer regime, NS to DS, for a neutral flow exists when B = 1.01 with the transition from DS to PS governed by t = 3 B .NS is not labelled in Fig. 12b but exists within the region when B ≤ 1.01.
The regimes maps indicate a direct relationship between t and B (or for Re L , which is indirectly related to t ) as the regimes begin to transition from NL / NS, DL / DS and PL / PS.This behaviour is as expected.B is a determining value for stratification levels.As B increases the flow will move from NL, DL and PL.The same principle applies to the SLD behaviour.

Regime map application
To apply the regime maps, values of B , Re and t are required and Eqs. 9 and 12 must be calculated.Here we give an example for the application of the regimes maps.Analysing data taken from a site at the Bourke Weir pool located on the Darling River in 1997, finds flow rates during late September, early October to be around 3.5 m 3 s −1 [35].Though surface heat flux over the Bourke Weir is difficult to find, average surface heat flux over Maude Weir can be approximated to be around ≈ 170 Wm −2 [36].Given a cross-sectional area of around 97 m 2 for the Bourke Weir and a depth of = 4 m , the flow velocity can be 1 3 calculated to be around U b = 0.04 ms −1 [35].The shear stress at the wall of the river can be determined by, where r p is the roughness of the bed channel which is equal to 0.05 [35].Here, the Reyn- olds number is equal to Re = U b ∕ = 1.4 × 10 5 .Analysis of average wind velocities at recording station in Bourke indicate minimal variations through the year with average speeds varying between 1.2 to 2.1 ms −1 between the period of 1991-1995 [37] taken at altitude 107 m.Using a logarithmic wind profile [38], U 10 can be interpolated to be around 1.7 ms −1 .From here, a surface shear stress can be found.
Turbidity at Bourke can range between 9-1740 NTU [39] though there is a strong correlation between the decline of turbidity and decreased flow rates [37].With turbidity at Bourke in August of 1995 dropping to around 20 NTU at flow rates 3.5 m 3 s −1 [37].20 NTU has been taken to calculate the attenuation coefficient .
With all the parameters listed above and summarised into Table 3, B , Re and t can be solved once Eqs. 9 and 12 are calculated.With these parameters, B = 17.6 , Re = 8.4 × 10 3 , t = 45 and Re L = 478 for the example of this flow.From Mitrovic et al.  [35], the period where U b was calculated (September-October of 1997), the flow was persistently stratified with persistent stratification defined when the temperature difference between the top and bottom of the water column is greater 0.5 • C for more than five days.Though temperatures at the surface of the water column will drop during the night, the difference between surface and bottom temperatures do not fall below 0.5 • C during this time.In the case of this flow, from the regime maps of Fig. 12, the flow will be in the NL regime for the LLD regime map and a PS regime for the SLD regime map.Overall, the flow for this instance will have no laminar layer depth within its diurnal timescale of t = 45 but will have a stratified layer depth that is persistent throughout the entire period of the diurnal timescale.Since the data indicates a persistent SLD, it is in agreement with the data from Mitrovic et al. [35] where they observe persistent stratification within this period.

Parameterisation of R * f as a function of Fr
This section aims to determine if previous parameterisation of the mixing efficiency can be applied to a diurnally heated, stratified channel flow and is distinct from previous sections with regards to the layer depths and their regimes.For many studies into stratified flows, there is a lot of emphasis placed on parametrisation of mixing with one such notable relation, the Froude number and flux Richardson number Fr − R * f framework.Previous data simulation have found the turbulent Froude number Fr to be a strong   parameter to parameterise the mixing efficiency and infer on the localised state of turbulence in stably stratified flow [16,23,32] though its relation has yet to be tested on diurnally heated open channel flows.For the paper, the irreversible flux Richardson number R * f is used as the mixing efficiency parameter.This definition of the flux Richardson number allows the formulation to quantify the mixing efficiency in stably stratified flows with the consideration of removing stirring, a large scale non-diffusive reversible flux from the equation [40].The R * f number can be defined as, where pe is the dissipation rate of the turbulent potential energy and is approximated by the density scalar variance dissipation rate p [41].Both parameters are defined respectively as: Here, is the density (temperature) and ′ i is the fluctuating density.For strongly stratified flows Fr ≪ 1 , the irreversible flux Richardson is approximately constant R * f ∝ Fr 0 , while moderately stratified flows Fr ∼ O(1) exhibit the relation R * f ∝ Fr −1 and for weakly stratified conditions Fr ≫ 1 ; R * f ∝ Fr −2 [16, 23, 32].The Fr-based framework defined by Garanaik & Venayagamoorthy [16] differs from our definition of SLD as the SLD is based on the length from the free surface when Fr ≈ 1 , making the location of the SLD within the transition region between weakly stratified and strongly stratified.Due to the nature of open channel flows and depth and time dependent heating function, it is expected for the paper's flow to exhibit a wide range of Fr and R * f spatially (in the vertical direction) and temporally.The temporal relevance is most notable for long diurnal timescales where the change in regime in one diurnal timescale cycle at a specific location is most evident and visible (Fig. 13d).
Figure 13 illustrates the irreversible flux Richardson number as a function of the turbulent Froude number for cases 2 to 5. Each point represents a specific time, while each shape indicates where the data lies along the channel height.The vertical dotted black line denotes where Fr = 1 while the solid black lines take the form of the function Fr −1 , the closely spaced dashed line shows the function Fr −2 and the horizontal and widely spaced dashed line when Fr 0 .The dot-dashed horizontal line indicate when R * f = 0.17 .Between Fig. 13a and d, the simulation data collapses reasonably well onto the parameterization made by Garanaik & Venayagamoorthy [16] with all cases supporting the theoretical critical R * f = 0.15 − 0.21 value above which turbulence is incapable to be sustained at steady-state [15,42,43].
With regards to the constant forcing case, the flow exists mostly within a moderately to strongly stratified regime as shown in Fig. 13a.For most of the constant forcing, the flow stays constant in its stratification regime along the vertical positions of the channel column except for when y = 0.3 .This can be noted in Fig. 13a where the y-position points stay rela- tively in the same Fr and R * f values as oppose to the scattering observed in Fig. 13d where the y = 0.3 red-circle points demonstrate differing Fr and R * f values.This scattering behaviour is seen in all diurnal cases within Fig. 13 where each region of the water column experiences a range of Fr numbers throughout its diel cycle, and this range broadens as t increases.Furthermore, the free surface is seen to undergo greater shifts in stratification levels than the near-wall regions between the diurnal periods.Figure 13d highlights all stratification regimes with the strongly stratified regime following the critical flux Richardson number value R * f ∝ Fr 0 .Where the case is moderately stratified, to the left of the vertical Fr = 1 line on the plot, the points follow the R * f ∝ Fr −1 equation well, and this additionally applies to the weakly stratified regime on the right of the Fr = 1 line which follows R * f ∝ Fr −2 adequately.This behaviour is similar to the time graphs and vertical profiles of t = 24 , revealing that at large t , the flow exhibits periods of stratification and an abrupt collapse of it for specific locations within the vertical water column.

Conclusion
This paper seeks to determine the effects of varying diurnal timescales on stably stratified channels as a model for radiatively heated river flow.This surface heating irradiance acts as potential energy and suppresses the turbulence in competition with the turbulent kinetic .Each marker represents a vertical position on the water column at a point in time with red circle marker at y = 0.3, blue diamond at y = 0.5, green squares at y = 0.6, purple crosses at y = 0.7 and orange triangles at y = 0.9 energy production through shearing from the bottom of the channel.In this flow, the near wall region is turbulent while the mid-channel and near surface region is stratified.As the thermal source is removed and re-introduced, the stratification alternatively increases and then decreases.
This paper provides a description of the local flow states of the open channel and show where they exist within the water column and their changes across the diel cycle as the flow goes through its cyclic heating period.The diurnal timescale between t = 6 to 24 with The flow is shown to exhibit either a fully neutral state, or a diurnal or persistently stratified state.It is essential to understand the conditions for the transition between these regimes.These conditions highlight whether the flow, within its diel cycle, will experience a breakdown of stratification or persist in its stratified state with minimal mixing.
A relationship was found between the diurnal timescale t and the vertical extent of the LLD and SLD.A regime map was produced with the transition from NL to DL given by t = 2.7 × 10 −10 Re 4.5 L − 0.1 , DL to PL by t = 1.56 × 10 2 Re −0.5 L − 18 , NS to DS by B = 1.01 and DS to PS by t = 3 B .This regime map may be of direct use in identifying flow condi- tions that may lead to adverse phenomena such as cyanobacterial blooms, caused by the rapid breakdown of a persistently stratified and stagnant flow [6,35].
For these flow conditions, the irreversible flux Richardson number R * f and its relationship with the turbulent Froude number Fr collapses well onto previous parameterizations of turbulent mixing in stratified flows [16].Increasing the diurnal timescale for this flow will broaden the range of stratification levels throughout a flow's diurnal cycle, and its thermal buoyancy influences extend further down the depth of the channel.At very high t the flow will move between a fully neutral state flow to a steady state equilibrium within its diurnal cycle.

Fig. 1
Fig. 1 Schematic of radiatively heated open channel flow

2 bFig. 2
Fig. 2 Comparison profiles of a full-and half-span channel for a constant forcing case with B = 5.8 , Re = 400 , Pr = 1 and = 8 : a stream-wise velocity, b instantaneous temperature difference to wall temperature w , c turbulent shear stress, d scalar flux, e buoyancy Reynolds number and f turbulent Froude number profile.Red solid lines = full-span and blue dashed lines = half-span domain

Fig. 3 Fig. 4 Fig. 5
Fig. 3 Time-series of a Ri B and b Fr B for Case INF and Case ISF with simulation parameters B = 2.9 , t = 3 , Re = 400 , Pr = 1 and = 8 .Red solid lines = neutral initial condition and blue dashed lines = stratified initial condition

Fig. 9
Fig. 9 Temperature difference Δ plotted against time t: a Case 2 -constant forcing case, b Case 4 -diurnal with t = 12 , c Case 18 -diurnal with t = 18 and d Case 5 -diurnal with t = 24

10
Fig. 11 Laminar and stratified layer depths for Case 6 and 7 with simulation parameters t = 6 , Re = 400 , Pr = 1 and = 8 : a Case 6 -B = 1.5 , and b Case 7 -B = 17.45 .Vertical axis is flipped to clearly show where layer depths lie within the water column.The dashed black lines signifies the surface radiant heat flux value I s that is scaled to a secondary y-axis defined by Eq. 2 located on the right of the plots.Red dotted lines indicate the LLD and blue solid lines indicate the SLD Fig.12 Regime map for an absent, diurnal or persistent a LLD and b SLD for all cases given in Table1.Figure a is plotted against Re L with the horizontal axis reversed and b is plotted against B .The blue squares indicate an absence of either the LLD or SLD, red circles indicate diurnal and black crosses are persistent laminarisation or stratification

Fig. 13 R
Fig. 13 R * f as a function of Fr for cases 2 to 5 with parameters B = 5.8 , Re = 400 , Pr = 1 and = 8.Each marker represents a vertical position on the water column at a point in time with red circle marker at y = 0.3, blue diamond at y = 0.5, green squares at y = 0.6, purple crosses at y = 0.7 and orange triangles at y = 0.9

B = 5 . 8
along with t = 6 with B = 1.5 and B = 17.4 for parameter values Re = 400 , Pr = 1 and = 8 are examined.Results of the temperature profile, laminar and stratified layers are shown and flow regimes are presented.

Table 1
Simulation cases and parametersCase

Table 2
Destratification time for each simulation where the superscript * represents the average parameter value at t = 0.75t for a given diurnal cycle while t d represents the estimate for destratification time defined by the scaling relation in Eq. 27 and t sim is the destratification time taken from the simulation results

Table 3
Summary of variables used to determine B , t and Re at Bourke Weir pool