Study of unsteady natural convection induced by absorption of radiation based on a three-waveband attenuation model

The present study considers unsteady natural convection induced by the absorption of radiation for possible applications in the water quality management for the shallow regions of lakes and reservoirs. The direct absorption of the incoming radiation by the water body forms a stable thermal stratification, whilst residual radiation reaching the bottom bathymetry is re-emitted as a boundary flux, forming an unstable thermal stratification, which is a potential source for a Rayleigh-Benard type instability. The bottom boundary layer instability drives intermittent vertical convection in the form of rising plumes. The plume rise is, however, limited by the stable thermal stratification due to the direct absorption, which is controlled by the attenuation coefficient of water. The attenuation coefficient is therefore an important parameter in determining the plume rise and the associated vertical mixing. The wavelength dependency of the attenuation coefficient of water is accounted for by using a three-waveband model. A theoretical prediction is made for the plume rise distance, which represents the region of vigorous mixing. Two-dimensional numerical simulation provides verification for the accuracy of the theoretical prediction.


Introduction
Natural convection has a significant impact on transport and mixing processes in the shallow regions of lakes and reservoirs [1,2]. The radiative thermal forcing from the sun induces both horizontal exchange flow and vertical mixing flow in the water reservoirs. The direct absorption of the incoming radiation by a water body follows Beer's law and decays exponentially with depth, and therefore forms a stable thermal stratification, whilst residual radiation reaching the bottom bathymetry is absorbed and then re-emitted as a boundary flux, forming an unstable thermal stratification in that region [3], which is a potential source for a Rayleigh-Benard type instability. A steady horizontal exchange flow is induced by the differential heating of water columns with different local depths, whilst the bottom boundary layer instability drives intermittent vertical convection in the form of rising plumes. The present study focuses on the latter flow mechanism in isolation, and therefore uses a horizontal fluid layer with the spatial variation of the bathymetry being discounted. Further, the present study considers an idealised and controlled situation in laboratory-scale experiments, rather than a full-scale field situation, and assumes a constant radiation intensity at the water surface with a 90° incident angle, and the heat exchange between the water surface and the ambient being negligible.
As expected from the studies of laminar plume behaviour in a quiescent environment with a stable thermal stratification [4,5], the plume rise is limited by the stable thermal stratification which forms due to the direct absorption of the incoming radiation. The stable stratification is governed by the decay rate of the radiation intensity, which is controlled by the attenuation coefficient of water. The attenuation coefficient is therefore an important parameter in determining the plume rise and the associated vertical mixing. While the attenuation coefficient depends on the turbidity of water and the wavelength of the radiation source [6], most of the existing and closely related literature assumes a single bulk value for this parameter. However, this simple assumption does not correctly model the absorption of the long wave radiation, which contains a significant portion of the total radiative energy. In the present study, a three-waveband attenuation model is developed based on a blackbody radiation distribution and the spectrum of the attenuation coefficient for water available in the literature [7]. A theoretical prediction is then made for the plume rise distance, based on a onedimensional analytical solution for the temperature increase. The theoretical prediction is validated by two-dimensional numerical simulation of the flow.

Problem formulation
The flow under consideration is governed by the incompressible Navier-Stokes equations with the Boussinesq approximation and a volumetric heating source term in the energy equation, representing the direct absorption of the incoming radiation. The quantities in the governing equations are normalised by the following scales: Length scale: , ℎ ~ The present study uses 0 = 10 −1 , which is a typical bulk value for the PAR (photosynthetic active range) of the spectrum measured in lakes [6]. It is noted that 0 is only used for the normalisation, hence the choice of 0 does not affect the results and discussions in this paper. 0 is the total surface intensity. is the gravitational acceleration, is the thermal expansion coefficient, 0 is the density, is the specific heat and is the thermal diffusivity. The thermal properties of water are assumed to be constant. The governing equations are then given in non-dimensional forms as follows: With the spectrum of the attenuation coefficient of water divided into N wavebands, and represent the average attenuation coefficient and average surface intensity over the nth waveband, respectively, i.e. 0 = ∑ =1 . Ra is the Rayleigh number �= 0 / 0 2 0 4 � and Pr is the Prandtl number (= / ) fixed at = 7 for water, where is the kinematic viscosity. ,2 is the Kronecker delta.
The present study assumes the initial condition of the flow to be stationary and iso-thermal. The radiation is instantaneously applied at = 0, and maintained thereafter. The fluid layer is bounded by a stress-free, adiabatic top ( = 0) and a rigid bottom ( = −ℎ) with a fixed flux condition associated with the re-emission of the residual radiation as mentioned earlier, which is given as: (2.4)

Three-waveband attenuation model and plume-driven mixing
The spectrum of the attenuation coefficient of water obtained from the literature [7,8] is shown in figure 1(a). The smallest attenuation coefficient is found in the visible range (250 − 700 ) of the spectrum, whilst longer waves are increasingly attenuated. Otanicar et al. [8] argues that the accuracy of their measurement is reduced for the visible range where the largest variation from the data of Hale and Querry [7] is found. The present study therefore uses the data of Hale and Querry [7]. In figure  1(b), the blackbody radiation distribution is shown for the theatre spot light (3200 ) used as a radiation source in the laboratory-scale experiment [9] and the solar radiation (5800 ). For the colour temperature of 5800 , a large portion of the intensity is carried by the visible, less attenuated waves, whilst for the colour temperature of 3200 , a large portion of the intensity is carried by the infra-red waves with increased attenuation. To determine the effect of the wavelength dependence on the attenuation coefficient, for each colour temperature, the spectrum of the attenuation coefficient is divided into = 50 wavebands of equal radiation intensity, therefore the variation of the radiation intensity in a water body follows (   During the initial stage of the flow development, the flow is quiescent and the temperature variation is therefore purely diffusive with equation (2.2) reduced to:  After the bottom boundary layer grows to a sufficient thickness, it becomes unstable to a Rayleigh-Benard type instability. Any small disturbance then grows exponentially with time, eventually leading to the formation of rising plumes [10]. As discussed earlier, the stable thermal stratification in the surface boundary layer limits the plume rise. In the present study, the term 'mixing depth' is used to quantify the plume rise distance from the bottom boundary since this length scale represents the region of vigorous mixing. The plume rise is considered to be determined by the location of neutral buoyancy ( ), which is obtained as follows: where ( , ) is the analytical solution (3.2) and ≡ + ℎ is the mixing depth. Equation (3.3) is solved numerically to obtain and hence ( = ( ) for a given h).

Numerical validation
The validation of the above theoretical prediction of the mixing depth is conducted via a twodimensional numerical simulation of an idealised laboratory-scale experiment using a horizontal fluid layer subject to constant and uniform radiation at the water surface supplied by a theatre spot light (3200 ) with the Rayleigh number in the range 10 7 ≤ ≤ 10 9 . The numerical simulation uses the SnS code [11,12], which is a non-staggered, Cartesian mesh, finite volume code based on a fractional step method [13,14], with the Adams-Bashforth and Crank-Nicolson time discretisation schemes being used for the advection and diffusion terms, respectively. The spatial discretisation of the diffusion terms uses the second-order central differencing, whilst the spatial discretisation of the advection terms uses the second-order central differencing with the ULTRA (Universal Limiter for Tight Resolution and Accuracy) flux limiter [15]. A variable time step is used to maintain the CFL (Courant-Friedrichs-Lewy) number within the range 0.15~0.3; the time step is updated at every time step based on the CFL number.
The horizontal extent of the computational domain is 0 ≤ ≤ 10 , and the vertical extent is −ℎ ≤ ≤ 0 with h varied in the range 0.5 ≤ ℎ ≤ 2.0 at a 0.5 interval. The boundary conditions imposed at the horizontal boundaries ( = −ℎ and = 0) are as aforementioned, while the periodic condition with =0 ≡ =10 , where f is any flow quantity, is applied at the vertical boundaries. 2000 grid points are uniformly distributed in the horizontal (x) direction. In the vertical (y) direction, the grid size near the horizontal boundaries is 4 × 10 −4 < ∆ < 5 × 10 −4 with 1%~3% linear stretching applied, giving at least 30 grid points within the bottom boundary layer after the flow becomes unsteady; the largest grid size is of 8 × 10 −3 < ∆ < 9 × 10 −3 in the centre of the domain; and the total grid points is 200 , 360 , 540 and 700 for ℎ = 0.5 , 1.0 , 1.5 and 2.0 , respectively. A grid dependence test has shown that further increasing the number of grid points has only marginal effect on the solution accuracy. A typical flow structure obtained from the numerical simulation is shown in figure 3. The figure contains instantaneous temperature and vertical velocity fields for the case of = 10 8 and ℎ = 1.0 at = 0.0286. The temperature field shows the formation of rising plumes from the bottom boundary layer and the stably stratified layer below the surface limiting the plume rise. The vertical velocity field also shows that the plume-driven convection is limited in the lower part of the flow domain and the upper part remains quiescent.
In figure 4(a), the theoretical prediction for the mixing depth, ≡ + ℎ, obtained by solving equation (3.3), is plotted over time for h in the range, 0.5 ≤ ℎ ≤ 2.0. It is shown that the mixing depth undergoes a rapid reduction with time and approaches a constant value for a given h, i.e. quasi-steady state. In figure 4(b), the vertical distribution of , where is the standard deviation of the vertical velocity over a horizontal line, = (−ℎ ≤ ≤ 0), is shown for Ra in the range, 10 7 ≤ ≤ 10 9 , and h in the range, 0.5 ≤ ℎ ≤ 2.0, at different times. The abscissa is normalised by the maximum value over the depth. The ordinate is translated by h to fix the location of the bottom boundary for different cases and is then normalised by the theoretical plotted in figure 4(a). It is shown in this figure that rapidly reduces as ( + ℎ)/ → 1 and is negligible for ( + ℎ)/ > 1 in all the cases. The accuracy of the theoretical prediction based on the location of neutral buoyancy in equation (3.3) is therefore verified.

Conclusions
Unsteady natural convection induced by the absorption of radiation has been investigated in the present study for possible applications in the water quality management for the shallow regions of lakes and reservoirs. An unstable thermal stratification which forms due to the absorption and reemission of residual radiation by the bottom bathymetry triggers convective instability, which then drives intermittent vertical convection in the form of rising plumes. The plume rise is, however, limited by a stable thermal stratification which forms due to the direct absorption of the incoming radiation and is controlled by the attenuation coefficient of water. The attenuation coefficient is therefore an important parameter in determining the plume rise and the associated vertical mixing. The wavelength dependency of the attenuation coefficient of water has been accounted for by using a three-waveband model for a given colour temperature of the radiation source. A theoretical prediction has been made for the plume rise distance, i.e. 'mixing depth', based on a one-dimensional analytical solution for the temperature increase. The corresponding two-dimensional numerical simulation has provided verification for the accuracy of the theoretical prediction.