Natural convection in a reservoir induced by sinusoidally varying temperature at the water surface

International Journal of Heat and Mass Transfer 134 (2019) 610–627

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Natural convection in a reservoir induced by sinusoidally varying temperature at the water surface Yadan Mao a,⇑, Chengwang Lei b, John C. Patterson b a b

Hubei Subsurface Multi-Scale Imaging Key Laboratory, Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China Centre for Wind, Waves and Water, School of Civil Engineering, The University of Sydney, NSW 2006, Australia

a r t i c l e

i n f o

Article history: Received 22 October 2018 Received in revised form 12 December 2018 Accepted 16 January 2019

Keywords: Scaling analysis Periodic thermal forcing Natural convection Nearshore water Thermal boundary layer

a b s t r a c t Subject to daytime heating and nighttime cooling, a horizontal temperature gradient is generated in nearshore waters owing to the gradually deepening topography offshore. The resulting natural convection has significant biological and environmental implications. The present investigation is concerned with natural convection in a reservoir model induced by a periodically varying surface temperature. A hybrid of semi-analytical approach coupled with scaling analysis and numerical simulation is adopted to characterize the periodic flow. The present scaling focuses on the long period scenario for which a quasisteady state is reached within a thermal forcing cycle. This is the case in natural water bodies under typical field conditions. Compared to the previously reported constant heating or cooling cases, the present scaling analysis is more complex and the developed scales depend on the period of the thermal forcing. Two critical functions of the Rayleigh number have been derived to identify the distinctness and stability of the thermal boundary layer. Distinct subregions are identified, and the scales of flow velocity and phase lag are derived for different subregions. These scaling results are verified by numerical simulations. Numerical results show that the flow response varies significantly with the length of period. For short periods, the large scale flow never reverses from cooling-induced flow. For the long period scenario, flow reversal occurs, with much more intense flow during the cooling phase than during the heating phase. The phase lag of flow response decreases with increasing period length. Finally, the present scaling results are applied to field situations, which agree well with field observations. Ó 2019 Published by Elsevier Ltd.

1. Introduction Adjoining the edge of landscape, the quality of nearshore waters tends to be severely affected by terrestrially derived solutes and particulates. Nevertheless, the water body is self-resilient to some extent due to exchange and mixing generated by various natural mechanisms, such as wind-generated currents and waves, tides, and circulation generated by salinity gradients in estuaries and coastal ocean. These physical processes reduce the residence time of the water body and have been reported comprehensively in reviews of physical limnology [1–3]. Under relatively calm conditions, another natural mechanism becomes prominent in refreshing the water and has attracted increasing recent interest: that is, natural convection resulting from temperature gradient. This mechanism is related to the intrinsic nature of the nearshore topography – an increasing water depth in the offshore direction. When subject to daytime heating ⇑ Corresponding author. E-mail address: [email protected] (Y. Mao). https://doi.org/10.1016/j.ijheatmasstransfer.2019.01.071 0017-9310/Ó 2019 Published by Elsevier Ltd.

or nighttime cooling, the shallow waters near the shore heat and cool more rapidly than the deeper waters offshore, leading to a cross-shore temperature gradient. This gradient generates a diurnally reversing cross-shore circulation, which is often referred to as a ‘thermal siphon’. 1.1. Field observations Field observations [4–11] have demonstrated that the thermal siphon plays a dominant role in driving cross-shore circulation in calm nearshore regions where wind-driven and tidal circulations are limited. It substantially reduces the residence time of nearshore water bodies and lowers the risk of algae blooms. Increasing interest has been paid to the biological and environmental impact of this cross-shore circulation. The nighttime convective circulation was identified to significantly enhance the phosphorus exchange between littoral and pelagic zones in a reservoir [7]. It was also observed to provide an important pathway for nearshore reef-bone material, such as phytoplankton, to the deeper strata of the Red Sea [10]. Moreover, the thermal siphon was sug-

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

611

Nomenclature A L f1(x) f2(x) f3(x)

g h k, n P p P0 Pr Q(x, t) Q(t) Ra Rac RaL T T0 t tc td

aspect ratio of the sidearm of the reservoir length of the reservoir (m) critical function of Rayleigh number for distinct thermal boundary layer critical function of Rayleigh number for unstable thermal boundary layer critical function of Rayleigh number for the long period scenario (a quasi-steady state can be reached within a thermal cycle) acceleration due to gravity (m2 s
1) the maximum water depth (m) non-negative integer non-dimensional period of the thermal forcing cycle non-dimensional pressure non-dimensional critical period for the long ramp scenario Prandtl number non-dimensional horizontal flow rate integrated over the local depth at offshore distance x and at time t non-dimensional horizontal flow rate integrated over the flow domain at time t global Rayleigh number critical Rayleigh number for instability local Rayleigh number of the thermal boundary layer temperature (K) reference temperature (K) non-dimensional time non-dimensional time for the thermal boundary layer to reach a steady state non-dimensional time for the thermal boundary layer to diffuse over the local depth

gested to be a general mechanism that enhances the connectivity between the near-reef waters and the open ocean [6,12]. In icecovered lakes where the effect of wind on flow is minimized, thermally driven flows are observed to play a significant ecologic role by transporting algae cells and dissolved oxygen [13–16]. A more comprehensive understanding of the thermal siphon requires more general approaches to quantify the flow and an isolation of this mechanism from other mechanisms concurrent in the field. Understanding of the flow mechanisms and quantification of flow properties have been achieved through various complementary approaches, including asymptotic solutions, scaling analysis, numerical simulations and laboratory experiments, and by assuming either a constant or a time-varying thermal forcing. 1.2. Theoretical and laboratory investigations 1.2.1. Constant thermal forcing For the daytime heating case, asymptotic solutions [17] and scaling analyses [18,19] have quantified natural convection induced by absorption of radiation. For the unstable regime, the stability properties at the initial and the early-transitional stages were revealed by linear stability analysis [17] and direct threedimensional numerical simulations [20], while those at the quasi-steady stage are revealed by spectral analysis [21]. In the absence of solar radiation such as on a cloudy day or at night, sensible heat transfer is still functioning to generate natural convection. This type of thermal forcing represented by an ideal constant isothermal surface heating was considered in Mao et al. [22] through scaling analysis and numerical simulations. For the nighttime surface cooling case, a two-dimensional wedge model was firstly proposed by Horsch and Stefan [23] to

u non-dimensional velocity component in the x direction umax, umin non-dimensional maximum and minimum u along the local water depth v non-dimensional velocity component in the y direction x non-dimensional horizontal coordinate x1 non-dimensional dividing position between distinct and indistinct thermal boundary layer x2 non-dimensional dividing position between stable and unstable subregions y non-dimensional vertical coordinate Greek symbols b thermal expansion coefficient (K
1) DT amplitude of the sinusoidally varying temperature at the water surface (K) Dt non-dimensional time delay of flow response dT non-dimensional thickness of the thermal boundary layer dV non-dimensional thickness of the viscous boundary layer j thermal diffusivity (m2 s
1) kk intermediate variable in the derivation of the semianalytical solution for s s non-dimensional temperature ss non-dimensional temperature variation across the thermal boundary layer q0 density at the reference temperature (kg m
3) h phase lag of flow response (rad) m kinematic viscosity (m2 s
1)

represent the nearshore geometry. Subsequently, a mixed-cellsin-series model was established and applied to field data to estimate the cooling-induced convective exchanges [24]. Later, the range of the Rayleigh numbers considered in Horch & Stefan [23] was extended numerically by Horsch et al. [25], in which different regimes of the developing flow were observed. The steady-state flow was reinvestigated by Sturman et al. [26] through scaling analysis, field and laboratory measurements. A different geometric configuration, with a shallow region gradually connected to a deep region by a sloping bottom, has also been adopted to study cooling induced stratification in lakes through laboratory experiments and field measurements [27]. Later, a detailed scaling analysis was provided for the cooling-induced flow by Lei and Patterson [28]. The proposed onset times of instability were subsequently verified by numerical simulations and laboratory experiments [29,30]. The scaling analysis of Lei and Patterson [28] was later extended by Mao et al. [31]. For both the daytime heating and nighttime cooling, the extended scaling analyses of Mao et al. [19,22,31] capture and quantify the variation of flow properties with offshore distance x. Two critical Rayleigh number functions with respect to x have been established to identify the distinctness and the stability of the local thermal boundary layer (TBL). Distinct flow subregions are identified and their extents are quantified. For the unstable regime, the dominant mode of heat transfer has been revealed to change from conduction to stable convection and finally to unstable convection as offshore distance increases. The scaling analyses of Mao et al. [19,31] were later extended by Dittko et al. [32] to suit the geometry of a shallow tetrahedron, and a three-dimensional numerical simulation was conducted to verify the results. Their simulation results show that the tetrahe-

612

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

dron configuration introduces some complex features such as helical flow. Using the real bathymetry data, large eddy simulation was conducted by Dittko et al. [33] for two real complex sidearms subject to constant radiation heating and constant cooling from a surface flux. They conducted separate simulations for the heating and the cooling phases and focused on the quasi-steady state flow. It was confirmed that many of the flow features revealed in the previous scaling and numerical studies based on simple geometry are still present in the complex geometry, and thus the previously adopted simple geometry is useful for revealing fundamental flow mechanisms. 1.2.2. Time-varying thermal forcing Different from the constant thermal forcing considered in the above investigations, the thermal forcing varies diurnally in the field. A periodic thermal forcing results in a phase lag between the flow response relative to the switch of the thermal forcing from heating to cooling or vice versa. To quantify this phase lag, Farrow & Patterson [34] and Farrow [35] simplified the diurnal model by introducing a uniformly distributed source term over the local depth to the energy equation. Based on a small bottom slope assumption, the zero-order asymptotic solutions of Farrow and Patterson [34] clearly revealed the transition from a viscousdominated flow in the shallow region to an inertia-dominated flow in the deeper part where the flow response significantly lags the forcing. It was found that the temperature response lags the thermal forcing by one quarter of a period, corresponding to 6 h in the diurnal cycle. Later, Farrow [35] generalized the problem with an arbitrary bathymetry and obtained second-order asymptotic solutions. Due to the complexity of the higher order solution, no explicit expression of the phase lag was obtained. These asymptotic solutions provide theoretical explanation for the phase lag and the transition of the dominant driving force with offshore distance. Nevertheless, the assumption of a vertically uniform distribution of thermal forcing is a considerable simplification. It means that the possible occurrence of a distinct TBL and the associated thermal instability are beyond the scope of those investigations. Subsequently, the numerical simulations of Lei & Patterson [36] adopted a thermal forcing that was composed of daytime radiation heating and nighttime surface cooling with the intensity varying sinusoidally with time. It is found that the estimated phase lag is significantly less than a quarter of the cycle in shallow water (compared to the penetration depth of radiation), whereas it is beyond a quarter of the cycle in deep water. Later, a different type of periodic thermal forcing, with temperature at the water surface varying sinusoidally with time, was considered in an experimental study [37] and concurrent numerical modelling [38]. In the above investigations, the period of the thermal forcing was set constant and thus the effect of the thermal forcing period on the flow response is yet to be revealed. Furthermore, the relevance of the fixed thermal forcing period considered in the previous investigations to the diurnal cycle in the field is still unclear. More recently, a periodic surface stress, representing the coastal sea-breeze/gully wind, was introduced into the thermally driven flow and its effect was investigated through asymptotic solutions [39] and numerical simulations [40]. The thermal source term was simplified to be uniformly distributed vertically. The effect of surface stress is revealed to be limited to a thin layer adjacent to the surface and the buoyancy effects dominate at depth for the parameters considered, consistent with the field measurements of Molina et al. [11]. As a bridge from a constant to a periodic thermal forcing, a simpler time-varying thermal forcing has been considered previously, in which the temperature at water surface increases linearly with time over a certain duration, i.e. the ramped duration, and then becomes constant [41]. A comparison between the time required

for the establishment of a quasi-steady state and the ramped duration results in two scenarios [41]. If the latter is larger, a quasisteady state will be reached within the ramped duration. For the present periodic thermal forcing case, there is a similar consideration, in which the period of the thermal forcing is the counterpart of the ramped duration. 1.3. The present thermal forcing Many heat transfer modes are functioning simultaneously between the water body and the atmosphere. Apart from solar radiation during the day, there are latent and sensible heat fluxes and the long-wave radiation into and out of the water body. All these forcings usually result in an approximately sinusoidal variation of the water surface temperature over diurnal cycles, modulated by fluctuations of the amplitude over a few days and superimposed on long-term temperature fluctuations (e.g. Fig. 8 in [5], Fig. 9 in [8], Fig. 3 in [6]). For the total surface heat flux, the diurnal variation is more complicated than a simple sinusoidal function and it varies strongly with field conditions (Figs. 3 and 4 in [6], Fig. 5 in [11]). For the present investigation, the thermal forcing is imposed by a prescribed water surface temperature that varies sinusoidally with time. Although in field situations the temperature at water surface also varies with offshore distance, the variation is assumed to be relatively small compared to the overall diurnal temperature variation. This simplification enables the present study to focus on the new features associated with periodic thermal forcing compared to constant and ramped thermal forcing. For future investigations, a more complex function of the surface heat flux tailored to field situations will be more appropriate for a closer representation of the thermal forcing. 1.4. Aims of the present investigation In summary, the previous investigations have shed some light on nearshore natural convection, especially for cases under constant thermal forcing, for which scaling analysis has identified different flow regimes and distinct subregions based on the characteristics of TBL [19,22,31]. Nevertheless, for the periodically varying thermal forcing, which is more realistic in the field, a comprehensive scaling analysis is still lacking owing to the complexity of the flow response. The previous quantification of the flow under periodic thermal forcing is achieved through asymptotic solutions [34,35], but the model adopts an oversimplified (vertically uniform) thermal forcing, which reveals no characteristics of the local TBL. Therefore, detailed features of the TBL under a periodic thermal forcing are still unclear. A periodically varying thermal forcing introduces two complexities to scaling analyses as compared with the constant thermal forcing: (1) While a steady or quasi-steady state will finally be reached under constant thermal forcing, that is not necessarily the case under periodic thermal forcing. Intuitively, it will be dependent on whether the period of the thermal forcing is sufficiently long for a quasi-steady state to establish. The effect of period length has not been covered by the previous numeral simulations [36,38] or experiments [37], since only a single period is adopted in those studies. (2) Owing to inertia, a phase lag is expected to be present between the flow response and the switch of the thermal forcing from heating to cooling or vice versa. A phase lag of a quarter of a cycle (equivalent to 6 h for the diurnal variation in field situations) has been predicted by asymptotic solutions [34], but field observations reveal that the phase lag of flow response can be more than 6 h [5]. This inconsistency requires further detailed quantification of the phase lag. The present study aims to provide insights into the above issues by integrating the semi-analytical solutions with scaling analysis

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

and numerical simulation. The period of thermal forcing will be classified as short or long depending on its comparison with the quasi-steady time of the TBL. The present numerical simulations will cover both the short and the long period scenarios. The hybrid approach aims to: (1) Identify the different flow regimes and distinct flow subregions based on the characteristics of the TBL; (2) Quantify the phase lag between the flow response and the thermal forcing. (3) Derive the relevant scales of natural convection under periodic forcing. 2. Problem statement and governing equations A two-dimensional (2D) reservoir model is adopted to capture the basic mechanism of the cross-shore exchange flow. It consists of two sections with an equal horizontal length of L/2: one bounded by a sloping bottom with a slope A (A = 2 h/L) and the other with a uniform depth of h (Fig. 1). With the Boussinesq assumption, the normalized 2D Navier-Stokes and energy equations governing the flow and temperature evolution are written as:

@u @ v þ ¼ 0; @x @y

ð1Þ

@u @u @u @p þu þv ¼
PrRa þ Prr2 u; @t @x @y @x

ð2Þ

@v @v @v @p þu þv ¼
PrRa þ Prr2 v þ PrRas; @y @t @x @y

ð3Þ

@s @s @s þu þv ¼ r 2 s: @t @x @y

ð4Þ

where Pr and Ra are the Prandtl number and the Rayleigh number respectively, defined as:

Pr ¼ m=j;

3

Ra ¼ gbDTh =ðmjÞ

ð5Þ

The respective scales for the normalization are: the length scale x, y  h; the time scale t  h =j; the velocity scale u, v  j=h; and the pressure gradient scale px, py  q0 gbDT. The non-dimensional temperature s is defined as ðT
T 0 Þ=DT. The temperature at the water surface T y¼0 is varying periodically with time, with a dimensionless period P, as given below: 2

T y¼0 ¼ T 0 þ DT sinð2pt=PÞ

ð6Þ

The water surface is assumed to be stress free (@u=@y ¼ 0, v = 0 at y = 0), whereas the sloping and flat bottoms are assumed adia^ ¼ 0), rigid and no-slip (u = 0, v = 0). An open boundary batic (@ s=@ n condition is adopted for the far end boundary (@u=@x ¼ 0, v = 0, @ s=@x ¼ 0 at x = 2/A). Any backflow through the far end boundary is at the reference temperature T0. This open boundary condition essentially caters for the exchange flow between the modeled domain and the much larger main body of the reservoir. It is worth noting that hereafter all of the mathematical expressions and quantities are presented in non-dimensional form unless mentioned otherwise.

613

3. Theoretical analysis As the first attempt to conduct scaling analysis for the thermal flow under periodic thermal forcing, the present analysis focuses on the flow development from the moment when the convective flow is about to reverse its direction to the time when convection becomes dominant. Over this period, convection is relatively weak and thus is negligible. In the absence of convection, we will first derive an asymptotic solution of one-dimensional conduction over the local water depth, and then carry out a scaling analysis based on the obtained asymptotic solution. 3.1. Semi-analytical solution A sinusoidally varying temperature at the water surface induces a warm/cold surface TBL when the temperature is increasing/ decreasing. Within this TBL, the temperature difference is greater vertically than horizontally and therefore vertical conduction dominates over horizontal conduction. Convection is negligible for conduction-dominated region near shore, and also when the convective flow is about to reverse. If convection and horizontal conduction are neglected, the problem can be simplified as a one-dimensional conduction problem with a variable local water depth of Ax:

@s @2s ¼ @t @y2

ð7Þ

s ¼ sin ð2pt=PÞ ðy ¼ 0Þ;

ð8Þ

@s ¼ 0 ðy ¼
AxÞ; @y

ð9Þ

s ¼ 0 ðt ¼ 0Þ:

ð10Þ

An isothermal water surface and an adiabatic bottom are specified by (8) and (9), respectively. Eq. (9) is the leading order ^ ¼ 0 (n ^ is normal to the sloping bottom) approximation of @ s=@ n under the assumption of small bottom slope A. An isothermal initial condition at the reference temperature T0 is embodied in (10). The Eqs. (7)–(10) can be solved through separation of variables, and the solution is given by


 1 2p 4p X 1 sin ðkk yÞ t
P AxP k¼0 kk "    # 2 k2k e
kk t
k2k cos 2Ppt
2Pp sin 2Ppt ;   2 k4 þ 2p  k 1P
kþ2 p where kk ¼ ; Ax

s ¼ sin

ð11Þ

where k is a non-negative integer. The exponential term in (11) approaches zero after the flow has experienced many cycles of peri2

odic forcing, so that e
kk t ? 0. If we assume that 2p=P is much smaller than the minimum of k2k , i.e. P  A2x2/p (In fact, we can always find a nearshore region where this is satisfied), then Eq. (11) can be further simplified as:



 1 2p 4p X 1 2p t ; sin ð k y Þ cos t þ k P P AxP k¼0 k3k  
k þ 12 p where kk ¼ Ax

s ¼ sin

Fig. 1. Sketch of the flow domain and the coordinate system.

ð12Þ

From this simplified form, the vertical temperature gradient @ s=@y at the water surface (y = 0) at the time when the temperature s there switches sign can be derived and approximated by the dominant term of k = 0.

614

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

@ s 16Ax cosð2pt=PÞ:  Pp @y

ð13Þ

Therefore, owing to the sinusoidally varying temperature s at the water surface, a cosinusoidally varying heat flux is imposed. Scale (13) is derived under the assumption of negligible convection. In this case, the time t in (13) applies only to the duration over which convection is not dominant. This is the context of the following scaling analysis and the later numerical verification, which is based on (13), 3.2. Scaling analysis 3.2.1. Scaling for the development of the thermal boundary layer With the periodic thermal forcing applied at the water surface, the growth of the TBL below the water surface is cyclic. When convection is negligible, heat conducted vertically into the TBL through the surface is larger than that convected away horizontally, and thus the TBL grows downwards. A balance between the unsteady term and the diffusion term in the energy equation (4) gives the scale for the thickness of the TBL:

dT  t 01=2

ð14Þ

The time t’ in (14) starts from when the TBL starts to grow. Owing to the cosinusoidal variation of the surface heat flux (refer to (13)), the continuous growth of the TBL can last for half a cycle at most, i.e. P/2, before which a quasi-steady state may be reached. A comparison between P/2 and the time it takes for the flow to reach a quasi-steady state classifies the flow into short and long period scenarios. In the following scaling analysis, we focus on the long period scenario for which a quasi-steady state is reached. This scenario is more relevant to the field situations as will be discussed in Section 5. The temperature within the TBL varies according to the surface heat flux in (13),

ss

16 Ax d cosð2pt=PÞ  p P T

ð15Þ

Since the bottom slope is assumed to be small to represent the field situations, the vertical velocity v within the TBL is negligible compared to the horizontal velocity u and so are its related terms in the vertical momentum equation (3). Therefore, the dominant balance in the vertical momentum equation (3) is between the pressure gradient term and the buoyancy term, which yields the scale for the pressure p,

16 Ax 2 d cosð2pt=PÞ: p p P T

ð16Þ

Within the horizontal momentum equation (2), a comparison among the unsteady inertia term Oðu=t 0 Þ, the advection term     O u2 =x , and the viscous term O Pru=d2T yields that the viscous 0 term dominates for Pr > 1 and ut =x < Pr. For the case with water (Pr = 7) at the early stage of the TBL growth, both conditions are satisfied. Therefore, in Eq. (2), the balance between the viscous term and the horizontal pressure gradient results in the scale for horizontal velocity:

u

16

p

A Ra d4T cosð2pt=PÞ: P

ð17Þ

3.2.2. Balance within the thermal boundary layer When the heat convected horizontally away from the TBL becomes sufficiently large to balance that conducted downward from the surface, a quasi-steady state is reached. For the long period scenario, this balance is achieved before the thermal forcing reverses and therefore a quasi-steady state is established. In the

energy equation (4), the balance between convection and conduction can be expressed as:

s

s

u  2: x dT

ð18Þ

After the balance is reached, the TBL will adjust its thickness to maintain the balance between conduction and convection within the TBL as shown in Mao et al. [41] and Patterson et al. [43]. Both studies have shown that, in order to maintain the balance, the TBL gradually becomes thinner after a quasi-steady state is reached. To simplify the discussion, we focus on the order of magnitude of flow properties and their dependency on the controlling parameters, such as Ra, P and x, but not their instantaneous variation with time. In this case, substituting the amplitude of the velocity scale (17) into (18), the thickness of the TBL can be derived as,

 p 16

dT 

16

1

x6 Ra
6 1


16 P : A

ð19Þ

Substituting (19) into (18) or (17), this thickness corresponds to a velocity scale of

u


13 16

p

2

1

x3 Ra3


13 A : P

ð20Þ

Therefore, the velocity at quasi-steady state increases with increasing Ra and x, but decreases with increasing P. Assuming that the TBL has been growing according to (14), the thickness of the TBL corresponds to a quasi-steady state time of

tc 

 p 13 16

1

x3 Ra
3 1


13 P : A

ð21Þ

According to the long period assumption, the quasi-steady state is reached within half of a cycle before the thermal forcing switches sign, i.e. tc < P/2, which transfers into

P > P0 ; where P0 

p12 2

1

x2 Ra
2 A
2 : 1

1

ð22Þ

It is clear that P0 decreases with increasing Ra. Therefore, the prerequisite for ‘the long period scenario’ is more relaxed as Ra increases. For the convection-dominated region, the TBL is distinct at the quasi-steady state. In other words, the thickness of the TBL is smaller than the local depth, i.e. dT < Ax, therefore,

Ra >

p 16

A
7 Px
5 :

ð23Þ

The right side of (23) is a critical function for Ra, and it is termed as f1(x)

f 1 ð xÞ 

p 16

A
7 Px
5 :

ð24Þ

For a given period P, f1(x) reaches the minimum at the maximum x over the sloping region (i.e. x = 1/A), which is A
2Pp/16. If Ra > A
2Pp/16, then Ra crosses with f1(x) at x1, the scale of which can be derived from (23) as

x1 

 p 15 16

7

Ra
5 A
5 P5 1

1

ð25Þ

This x1 marks the switch from an indistinct TBL subregion to a distinct TBL subregion as x increases. Its value decreases with increasing Ra, but increases with increasing P. For the nearshore conductive region, the TBL reaches the local bottom before convection becomes sufficiently strong to balance conduction. A quasi-steady state is reached locally when the TBL reaches the bottom [41]. Therefore, the thickness of the TBL at the quasi-steady state is

615

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

dT  Ax:

ð26Þ

And thus the time it takes the TBL to reach the local depth is

t d  A2 x2

ð27Þ

For the long period scenario, this is achieved within half a cycle, i.e. td < P/2. Using (27), this becomes,

P > 2A2 x2

ð28Þ

Substituting (26) into (17), the scale for the quasi-steady state velocity within the conductive region is:

u

16

p

Ra

A5 x4 cosð2pt=PÞ P

ð29Þ

3.2.3. The phase lag of flow response Intuitively it is expected that the inertia of the water body will result in a phase lag between the flow reversal and the switch of the thermal forcing. When the thermal forcing is about to change sign, it approaches zero and thus the pressure gradient term in the momentum equation generated by the buoyancy term is negligible. The previous buoyancy-viscous balance is replaced by inertia-viscous balance within the viscous boundary layer (VBL). Therefore, the balance in the horizontal momentum equation (2) is between the unsteady inertia term and the viscous term

u u  Pr 2 : Dt dV

ð30Þ

where Dt is the time scale over which the balance occurs and dV is the thickness of the viscous boundary layer, scaled as [18]. For the conduction-dominated region under the scenario, the VBL has reached the local depth before of thermal forcing, and therefore dV  Ax. Substituting the time delay of the flow response is derived as,

Dt 

A2 x2 : Pr

dV  Pr 1=2 dT long period the reverse it into (30),

ð31Þ

Multiplying (31) by 2p/P yields the corresponding phase lag h given by

h

2pA2 x2 : PPr

ð32Þ

It is clear from (32) that the phase lag h in the conductiondominated region increases with offshore distance x, but decreases with increasing thermal forcing period P. It is independent of the Rayleigh number Ra. This makes sense as convection is negligible in the conduction-dominated region and Ra is a measure of the strength of convection. For convection-dominated region, the thickness of the TBL, i.e., dT is specified by scale (19). Using the relation of dV  Pr 1=2 dT [18] and substituting (19) into (30) yields the time scale for the delay of the flow response as


13  p 13 1 1 d2 P Dt  V  d2T  x3 Ra
3 : A Pr 16

ð33Þ

The corresponding phase lag h is derived by multiplying scale (33) by 2p/P 1

h  2
3 p3 x3 Ra
3 A
3 P
3 : 1

4

1

1

2

ð34Þ

Therefore, the phase lag of the flow response decreases with increasing Ra and P, but increases with offshore distance x. For the long period scenario, the duration of the heating phase is sufficiently long to counteract the flow induced by the preceding cooling phase, so that the large-scale flow reverses well before the

heating phase ends. In this case, the time delay Dt is much less than half of the thermal forcing period, i.e.

Dt < P=2 or h < p:

ð35Þ

Substituting (33) or (34) into (35) leads to the same criterion for ‘the long period scenario’ as that specified by the condition of (22). Therefore, the requirement for a quasi-steady state to occur coincides with the requirement for flow reversal. 3.2.4. Instability for the convective region If the local Rayleigh number RaL exceeds a certain critical value Rac, that is, RaL > Rac, where Rac  657.5 based on a free-free boundary configuration [42], the Rayleigh-Bénard type instability will be present in the local flow. For the long period scenario, a quasi-steady state is reached, and therefore the local Rayleigh number RaL can be calculated using the scale of ss and dT specified in (15) and (19) respectively: 3

3

gbðDT ss ÞðhdT Þ gbDTh ss d3T   Rass d3T mk mk
13
13 1 16 A 5  x3 Ra3 : P p

RaL 

ð36Þ

It is worth noting that ss and dT are non-dimensional parameters. Multiplying them by DT and h respectively recovers the dimensional values of the temperature difference across the TBL and the thickness of the TBL. Local instability will be present if RaL > Rac. Therefore, the criterion for the occurrence of instability can be derived from (36) as

Ra >

p P
5 3 x Rac 16 A

ð37Þ

Similar to the criterion for the presence of a distinct TBL specified in (23), the right side of (37) is also a critical function of Ra, which varies with the offshore distance x, and is denoted by f2(x) as:

f 2 ð xÞ 

p P
5 3 x Rac : 16 A

ð38Þ

Similar to f1(x), f2(x) also increases with P but decreases with increasing x. For a fixed period, the minimum of f2(x) is reached at x = 1/A, which is PA4Rac3p/16. If Ra > PA4Rac3p/16, then Ra crosses with f2(x) at x2, the scale of which can be derived from (37) as follows:

x2 

 p 15 16

1

Ra
5


15 3 P Ra5c : A

ð39Þ

This x2 marks the switch from a stable nearshore region to an unstable offshore region. Similar to x1, x2 also decreases as Ra increases or as P decreases. 3.3. Prerequisites and possible flow regimes So far, two criteria have been derived to characterize the local boundary layer flow: (1) The TBL is distinct if Ra > f1(x). (2) The flow is unstable if Ra > f2(x). Both f1(x) and f2(x) decrease as x increases or as P decreases. A comparison between f1(x) and f2(x) yields two possible scenarios depending on the bottom slope A (Fig. 2): (a) For A > Rac
1/2, f1(x) < f2(x) for all x (Fig. 2a). (b) For A < Rac
1/2, f1(x) > f2(x) for all x (Fig. 2b).

616

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

Fig. 2. Profiles of f1(x), f2(x) and f3(x) for (a) A > Ra-1/2 , the result of A = 0.1 is plotted as an example (b) A < Ra-1/2 , the result of A = 0.01 is plotted as an example. c c

For the assumed value of Rac  657.5, the threshold of the bottom slope, i.e. Rac-1/2, is approximately 0.04. Although the prediction of two scenarios are similar to the cases of constant heating and cooling, the periodicity of the thermal forcing introduces more complexity; for instance, an additional parameter P appears in the two critical functions. Furthermore, the periodicity brings a precondition for the critical functions, the derivation of f1(x) and f2(x) is based on the assumption of long period during which a quasi-steady state is reached, as specified by (22). From (22), it can be derived that for a certain P, the long period scenario requires that,

Ra >

p 2

P
2 A
1 x:

ð40Þ

The right side of (40) is a function of x, which is termed as f3(x),

f 3 ð xÞ 

p 2

P
2 A
1 x:

ð41Þ

Contrary to f1(x) and f2(x), f3(x) increases with x and reaches the maximum value A
2P
2p/2 at x = 1/A. A quasi-steady state can be reached within the convective region if Ra > f3(x), and this is the precondition for the derivation of f1(x) and f2(x) and the associated flow properties at the quasi-steady state. Therefore, for f1(x) and f2(x) to be meaningful, Ra must be larger than f3(x). The discussion of flow regimes will focus on scenario (a) described above (i.e. A > Rac
1=2 ). Discussion of scenario (b) can be made following the same procedures outlined in Mao et al. [19,31]. For A > Rac-1/2, f1(x) < f2(x) for all x. Both f1(x) and f2(x) reach their minimum at x = 1/A, which are A
2Pp/16 and PA4Rac3p/16 respectively. f1(x) and f2(x) are derived based on the assumption of a quasi-steady state, which requires that Ra > f3(x). Comparing the maximum of f3(x), i.e. A
2P
2p/2, with the minimum of f1(x) and f2(x) results in three possible relations as P varies (Fig. 3). (1) For P > 2, f3(x) is smaller than both f1(x) and f2(x) for all x (Fig. 3a). (2) For 2A
2 Ra
1 c < P < 2, f3(x) crosses with f1(x) (Fig. 3b). (3) For P < 2A
2Rac-1, f3(x) crosses with both f1(x) and f2(x) (Fig. 3c). The two minimum values of f1(x) and f2(x), which are A
2Pp/16 and A4PRac3p/16 respectively, determine three possible flow regimes depending on the global Rayleigh number. The low (Ra < A
2 P p=16), medium (A
2Pp/16 < Ra < PA4Rac3p/16) and high Ra regimes (Ra > PA4Rac3p/16) are represented by lines L1, L2 and L3 respectively in Fig. 3.

(i) In the low-Ra regime of Ra < A
2Pp/16, represented by L1 in Fig. 3, Ra is smaller than both f1(x) and f2(x), and therefore, the TBL is indistinct and the flow is stable. Conduction dominates over convection over the entire domain. The value of f3(x) is irrelevant here since it determines whether a quasisteady state is reached within the convective region, and thus its value does not affect the conductive region. The quasi-steady time is described by scale (27), independent of Ra. (ii) In the medium-Ra regime of A
2Pp/16 < Ra < PA4Rac3p/16, represented by L2 in Fig. 3, Ra < f2(x) for all x and therefore the flow remains stable. The intersection x1(L2) between Ra and f1(x) divides two subregions: (1) for x < x1(L2), Ra < f1(x), the TBL is indistinct and conduction dominates. (2) For x > x1(L2), Ra > f1(x), the TBL becomes distinct and stable convection dominates. Nevertheless, f3(x) crosses with f1(x) if P < 2, and therefore it should be noted that any dividing position below f3(x) should not be taken as a real one, since the prerequisite for a quasi-steady state flow (Ra > f3(x)) is not satisfied. A typical example is shown in Fig. 3c, where the dividing position x1(L2) should not be taken as a meaningful one. (iii) In the high-Ra regime of Ra > PA4Rac3p/16, represented by L3 in Fig. 3, Ra is larger than the minimum of f2(x) and so it intersects with both f1(x) and f2(x) at x1(L3) and x2(L3) respectively. If these two dividing positions are over f3(x), then they divide the flow domain into three subregions: for x < x1(L3), Ra < f1(x) < f2(x), conduction dominates with an indistinct TBL; for x1(L3) < x < x2(L3), f1(x) < Ra < f2(x), stable convection dominates with a distinct TBL; for x > x2(L3), Ra > f2(x), the flow becomes unstable. Similarly, any dividing position below f3(x) should not be taken as a relevant one since the prerequisite for f1(x) and f2(x) to be meaningful is not satisfied. In the following numerical verification of the quasi-steady state scaling, these irrelevant dividing positions and regions where a quasi-steady state cannot be reached have been avoided. 4. Numerical verification of the scaling analysis 4.1. Numerical schemes In the numerical simulations, the water body is initially set to be stationary and isothermal at the reference temperature T0. The normalized governing Eqs. (1)–(4) along with the initial and boundary conditions specified in Section 2 are solved numerically using a finite-volume method. The SIMPLE (semi-implicit method

617

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

Fig. 3. Profiles of f1(x), f2(x) and f3(x) for scenario (a) (A > Rac-1/2, A = 0.1 is taken as an example) at three representative periods (a) P > 2 (P = 3 is taken as an example), f3(x)
2 < f1(x) < f2(x) for all x (b) 2A
2 Ra
1 Rac-1 (P = 0.1 is taken as an example), f3(x) crosses with both f1(x) c < P < 2 (P = 0.5 is taken as an example), f3(x) crosses with f1(x). (c) P < 2A and f2(x).

for pressure-linked equations) scheme is adopted for pressurevelocity coupling; and the QUICK (quadratic upstream interpolation for convective kinematics) scheme is applied for spatial derivatives. A second-order implicit scheme is applied for time discretization in calculating the transient flow. The simulation is conducted with a non-dimensional depth of 1, a fixed bottom slope of A = 0.1, a horizontal length of L = 20, and a fixed Prandtl number of Pr = 7. The section with a uniform depth is set to be of the same length as the sloping section. A mesh and time-step dependency test has been conducted using four different meshes, 300  50, 600  100, 900  150, 1200  200 for the case with the highest Ra, i.e., Ra = 2  107, and P = 0.5. The time-steps for different meshes are adjusted accordingly so that the Courant–Friedrichs–Lewy (CFL) number remains the same for all meshes. The maximum velocities umax over the local water depth at x = 7.0 obtained with different meshes are listed in Table 1 with the value at a specific time (t = 0.25P) for the stable heating phase, and with the averaged value (over t = 0.7P
0.8P) for the unstable cooling phase. The errors of the first three meshes are estimated based on the results of the finest mesh. As mesh size increases, the numerical error reduces and the value of umax converges. The effect of mesh size on the solution is expected to be even less significant for other cases with lower Ra. It is worth noting that the numerical errors during the cooling phase are much larger than during the heating phase (Table 1). This shows that the mesh size has a stronger effect on simulation accuracy for the cooling phase owing to the occurrence of instability. Based on these mesh tests, a mesh of 900  150 and a dimensionless time-step of 3.2  10
6 are adopted for the simulation, giving a maximum CFL number of approximately 0.89. To avoid singularity at the tip, a very small tip region (x < 0.02) was cut off and an extra vertical adiabatic wall was assumed. Similar treatment of the tip region was carried out in previous investigations [19,21,22,31,41], and no significant modification to the flow was observed. The flow domain is meshed with a non-uniform grid with an increasing grid density toward all the boundaries. For the chosen mesh, the minimum face area is 0.0044; the maximum face area is 0.5668; and the maximum grid-stretching factor is 1.04. Since the numerical simulations start with an initially isothermal and stationary flow domain, a start-up effect is expected in the simulations. As will be shown later, the start-up effect is pre-

sent over the first few cycles only, and the longer the thermal forcing period is, the fewer cycles the start-up effect will affect. The numerical results for verification are selected after the start-up effect is minimized. 4.2. Definition of flow properties To characterize the variation of thermal features with x, the horizontal heat transfer rate H(x, t) is calculated as:

Hðx; t Þ ¼

1 Ax

Z


 @s us
dy @x
Ax 0

ð42Þ

where H(x, t) is normalized by jDT=h. The first and the second terms in the bracket of Eq. (42) represent the contributions from convection and conduction respectively. To quantify the intensity of the horizontal exchange flow at a specified location x, the horizontal exchange rate Q(x, t) at position x is calculated as:

Qðx; tÞ ¼

1 2

Z

0


hx

jujdy;

ð43Þ

where hx is the local water depth at the horizontal position x. The average volumetric horizontal flow rate Q(t) is obtained by integrating Q(x, t) horizontally over the length L:

QðtÞ ¼

1 L

Z

L

Q ðx; t Þdx ¼ 0

1 2L

Z

L

0

Z

0


hx

jujdydx

ð44Þ

It is clear that the volumetric flow rate Q(t) describes the intensity of the horizontal exchange flow. 4.3. Different flow regimes 4.3.1. Ra < A
2Pp/16 In this low Ra regime, scaling has predicted that conduction dominates over the entire domain. The critical Ra value, A
2Pp/16, is around 196, and thus a simulation at Ra = 20 is conducted. For a quasi-steady state to be reached, the long period criterion specified by (28) (P > 2A2x2) has to be met, which transfers into P > 2 for the entire domain to be quasi-steady. This criterion is met by choosing a period of P = 10. The numerical results confirm the dominance of

Table 1 Effect of mesh size and time step on the maximum velocity umax at x = 7.0 over the local water depth for Ra = 2  107, P = 0.5. Mesh

300  50

600  100

900  150

1200  200

Time step umax (t = 0.25P) Error (%) (t = 0.25P) umax averaged over t = 0.7P
0.8P Error (%) (t = 0.7P
0.8P)

9.6  10
6 1946.20 0.070 28355.51 4.31

4.8  10
6 1945.40 0.029 28875.82 2.56

3.2  10
6 1944.94 0.0057 29384.29 0.84

2.4  10
6 1944.83 – 29631.97 –

618

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

conduction over the entire flow domain at different stages of the flow development (Fig. 4). Scaling (29) predicts that in the conductive region, the absolute flow speed reaches the maximum at t = nP and t = nP + 0.5P (n is a non-negative integer), and the minimum at t = nP + 0.25P and t = nP + 0.75P. Consistent with this prediction, the streamlines are much denser at t = nP + 0.5P and n + P than at t = nP + 0.25P and nP + 0.75P (Fig. 4). Since the results are from cycles when the start-up effect is minimized, we add nP in time t to indicate generalization of the results. 4.3.2. A
2Pp/16 < Ra < PA4Rac3p/16 In the medium Ra regime, scaling has predicted two flow subregions: a conduction-dominated subregion nearshore and a convection-dominated subregion offshore. The flow is expected to be stable over the entire domain. These predictions agree with the numerical results shown in Fig. 5. Scale (25) (x1  ðp=16Þ1=5 Ra
1=5 A
7=5 P 1=5 ) predicates that the dividing position x1 moves toward shore as Ra increases or as P decreases, which agrees with the trend shown by the simulation results (Fig. 5a–b, b–c). The dividing position x1 is derived under the assumption of the long period scenario, for which a quasi-steady state is reached before the thermal forcing reverses. Therefore, numerical data for verification have to meet the long period criterion of P > P0 (P0  ðp=2Þ1=2 x1=2 Ra
1=2 A
1=2 ), as specified in (22). Since P0 decreases with increasing Ra, a wider range of P can be selected as Ra increases. The values of x1 obtained from numerical simulations under different Ra and P are plotted against the scaling predictions in Fig. 6. The clear linear correlation demonstrates that x1 is well predicated by the scaling analysis. For the unstable regime of Ra > PA4Rac3p/16, such a dividing position x1 is also present (Fig. 7). Therefore, data in Fig. 6 include those obtained in the unstable regime. 4.3.3. Ra > p/16PA4Rac3 In the high Ra flow regime, scaling has predicted the coexistence of three subregions, with the dominant mode of heat transfer switching from conduction to stable convection, and then to unstable convection as x increases. The numerical results of the horizontal heat transfer rate H(x, t) at the start of 30 consecutive thermal forcing cycles are plotted in Fig. 7. The lower panel of Fig. 7 zooms in the nearshore region of the upper panel. It is clear that the dominant mode of heat transfer changes from conduction to convection as x increases. The dividing positions x1 predicated by scaling have been verified by simulation results (Fig. 6). Moreover, scaling has predicted that for the unstable flow regime, unstable convection

will eventually dominates as x increases to the unstable subregion (x > x2). The predicted presence of instability in the form of thermal plumes corresponds well with the observed oscillation of the horizontal heat transfer rate with x in Fig. 7. The oscillation is shown to increase with increasing x for all the parameters considered (Fig. 7a–c). This is expected since scaling (36) (RaL  ð16=pÞ1=3 x5=3 Ra1=3 ðA=P Þ1=3 ) has predicted that the local Rayleigh number RaL increases with increasing x, and therefore the local flow becomes more unstable as x increases and thus resulting in larger oscillations. A comparison between Fig. 7a and b suggests that the oscillation of the heat transfer rate increases with increasing Ra, while a comparison between Fig. 7b and c shows that the oscillation increases with decreasing P. These results are consistent with the scaling predication, since the local Rayleigh number RaL, specified by scaling (36), increases with increasing Ra or decreasing P, resulting in a more unstable flow and therefore larger oscillations. Consistent with this oscillation with x, the heat transfer rates at the same phase of different cycles begin to diverge due to the initiation of instability as x increases. The horizontal location where they begin to diverge marks the start of the unstable subregion, and thus represents the dividing position x2 between the stable and unstable region. To quantify the degree of divergence, the standard deviation of the calculated convective heat transfer rates at the start of the 30 consecutive cycles is calculated at different x and shown in Fig. 8a. The positions where the standard deviation exceeds a certain threshold (0.01) approximately indicate the dividing position x2 between the stable and unstable subregions. These positions obtained from the simulation results are plotted against the scaling prediction in Fig. 8b and are projected in Fig. 7 (lower panel). The good linear correlation in Fig. 8b verifies the scaling prediction of x2. Scaling (39) predicts that x2 moves toward shore as Ra increases or as P decreases, which agrees with the trend shown in Fig. 7. It is found that a variation of the threshold over an order of magnitude, e.g. from 0.001 to 0.01, does not affect the linear correlation. Since the scaling is catered to the long period scenario, the criterion of P > P0 (P0  (p/2)1/2x1/2Ra
1/2A
1/2), as specified in (22), has to be met. As Ra increases, the value of P0 decreases and thus the results of smaller P can be selected. Meanwhile, as Ra increases, larger P can be selected to meet the unstable criterion of Ra > p=16PA4 Ra3c . Therefore, a wider range of P is selected for verification as Ra increases (Fig. 8b). The horizontal heat transfer rate in Fig. 7 changes from positive to negative as x increases. This is because that as specified by scale (34), the phase lag increases with x, and therefore it takes longer time for the flow to reverse as x increases. Before the start of the

Fig. 4. Flow features at the low Ra regime (Ra < A
2Pp/16) with Ra = 20, P = 10. The horizontal heat transfer rate averaged over the local depth (the upper panel), and the isotherms with an interval of 0.02 overlaid by streamlines with an interval of 2  10
4 (the lower panel). The solid and dashed streamlines represent anti-clockwise and clockwise flow, respectively. The vertical scale on the contour plots is elongated four times as compared with the horizontal length.

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

619

Fig. 5. Flow features at the medium Ra regime (A
2Pp/16 < Ra < PA4Rac3p/16) at the start of a new thermal forcing cycle at the quasi-steady state. The horizontal heat transfer rate averaged over the local depth (the first row), the isotherms with an interval of 0.05 overlaid by streamlines with an interval of 0.05 (the second row) for (a) Ra = 2000, P = 5; (b) Ra = 104, P = 5; and (c) Ra = 104, P = 2.5.

Fig. 6. The dividing position x1 between the conduction-dominated subregion and the convection-dominated subregion derived from numerical simulations versus that predicated by scale (25). The data cover both the medium and the high Ra regimes.

current thermal cycle, i.e. when s at y = 0 switches to positive, the temperature at the water surface has already been increasing for a quarter of a cycle, and thus the nearshore region has already adjusted to this positive heat flux and reversed, whereas the further offshore region is still flowing anti-clockwise owing to inertia. 4.4. The effect of thermal forcing period on the flow response Scaling analysis has predicted two scenarios depending on the thermal forcing period P. For short periods (P < P0, P0  (p/2)1/2x1/2Ra
1/2A
1/2), a quasi-steady state cannot be reached, and the phase lag is longer than half of the cycle, meaning that a counter-rotating flow cannot be generated after the temperature at the water surface switches sign. The reverse is true for long periods.

Fig. 9 shows the isotherms and streamlines over a typical cycle for three different periods at Ra = 3  106. For the very short period of P = 0.0125, instability is present almost over the entire cycle as revealed by the wavy streamlines at all stages (Fig. 9, the left column). Over the entire heating phase (from t = nP to nP + 0.5P), the solid streamlines show that the large-scale circulation is still down the sloping bottom and toward the shore at the surface, which is typical of the flow induced by surface cooling. Therefore, the cooling-induced gravity currents over the previous half cycle persist over the entire heating phase. This persistence is also illustrated by the horizontal velocity u monitored near the far end of the sloping bottom at x = 9 (Fig. 10), which is irregular and almost always positive (in the offshore direction). Therefore, for the short period cases, cooling-induced flow dominates over the entire ther-

620

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

Fig. 7. Horizontal heat transfer rate averaged over the local depth in the high Ra regime (Ra > p=16PA4 Ra3c ) at the start of 30 consecutive thermal forcing cycles. Dashed and solid lines represent the heat transfer rates by conduction and convection respectively. The second row zooms in the near shore region. x1 is where conduction equals convection. Position of x2 is determined from Fig. 8, where the standard deviation exceeds the threshold of 0.01.

Fig. 8. (a) The standard deviation of horizontal convection averaged over the local depth at the start of 30 consecutive cycles. (b) The dividing positions x2 predicated by scaling (39) versus that from numerical simulations.

mal cycle over most of the flow domain. During the heating phase (from t = nP to nP + 0.5P), the streamlines show a gradual disappearance of small flow structures and a formation of large structures, indicating a decreasing intensity of flow instability. Subsequently, at t ¼ nP þ 0:75P, when the water surface temperature reaches the minimum, cold plumes form beneath the water surface. Strong thermal plumes are still observed near the end of the cooling phase at t = nP + 0.9P. As the period P increases to P = 0.125, the cooling-induced flow becomes less persistent in the heating phase, since the streamlines indicate that cooling-induced flow dominates only during the initial stage of the heating phase (from t = nP to nP + 0.1P, Fig. 9, the middle column). At t = nP + 0.25P, the dashed streamlines demonstrate that the flow within the sidearm has already reversed to

heating-induced flow structure, with an up-slope flow near the bottom and an offshore flow in the surface layer. In this case, the thermal forcing period P is sufficiently long to allow for a flow reversal. The reverse of the flow is also recorded by the nearbottom horizontal velocity u at x = 9 (Fig. 10), which switches sign at about t ¼ nP þ 0:2P. After t = nP + 0.25P, temperature at the water surface starts to decrease, and therefore a cold surface layer is generated beneath the water surface, which results in instability represented by thermal plumes at t = nP + 0.5P (Fig. 9 middle column). After t = nP + 0.5P, the cooling phase starts. The streamlines at t = nP + 0.75P and t = nP + 0.9P for P = 0.125 (Fig. 9 middle column) are smoother and consist of less small flow structures than those for the case of P = 0.0125 (Fig. 9 left column), suggesting a reduced flow instability as P increases. This is because that for suf-

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

621

Fig. 9. Isotherms and streamlines at different stages of a typical thermal forcing cycle for Ra ¼ 3  106 under different forcing periods: P = 0.0125 (left column), P = 0.125 (middle column), P = 1.25 (right column). The solid and dashed streamlines denote anti-clockwise and clockwise flow respectively.

Fig. 10. The horizontal velocity u near the end of the sloping region and close to the sloping bottom, at x = 9, y =
0.867 for cases with different forcing periods, Ra = 3  106.

ficiently large P, the cold plunging plumes that induce the small flow structures have sufficient time to merge into intense largescale flow down the bottom slope.

As P increases further to P = 1.25, the cooling-induced flow can hardly be observed during the heating phase. At an very early stage of the heating phase, i.e. t = nP + 0.1P, the flow has already reversed, as shown by both the dashed streamlines (Fig. 9, the right column) and the near-bottom flow velocity u (Fig. 10), in contrast to the clear presence of cooling-induced flow for the case of P = 0.125 at the same stage. Therefore, as P increases, the phase lag decreases, consistent with the scaling prediction specified in (34). Furthermore, at the start of a thermal cycle (t ¼ nP), the flow has already become stable with no wavy feature in the streamlines (Fig. 9). This is owing to the increased critical Ra required for the occurrence of instability as P increases, agreeing with the scaling relation specified by (37). From t = nP + 0.1P to nP + 0.25P, the density of the streamlines increases, suggesting an increased intensity of heating-induced flow with increasing water surface temperature. At the end of the heating phase, i.e. at t = nP + 0.5P, the presence of instability represented by the thermal plumes is also observed, similar to the case of P = 0.125, Subsequently, at t = nP + 0.75P, the consecutive streamlines reveal that a large scale flow is already formed across the entire reservoir, much earlier than the shorter period case of P = 0.125, for which the large scale flow across the reservoir is observed at a later stage, i.e. at t = nP + 0.9P. The horizontal exchange rate Qx(x, t) over a typical thermal forcing cycle is illustrated in Fig. 11 for different periods. The occurrence of plunging cold-water plumes is reflected by the wavy

622

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

undulations of Qx(x, t) with x. Smoother variation of Qx(x, t) with x suggests reduced intensity of plumes and the formation of large-scale continuous flow. For cases with short P (P = 0.0125), the cooling phase (t = 9.5P
10P) is dominated by plunging thermal plumes which do not have sufficient time to merge into a large-scale flow, as evidenced by the strong undulations of Qx with x (Fig. 11a). As P increases to P = 0.125, the cold plumes during the cooling phase eventually develop into a large scale flow, with Qx(x, t) varying rather smoothly with x before the end of the cooling phase (Fig. 11b). As P increases further to P = 1.25 (Fig. 11c), the large-scale down-slope flow is formed at a much earlier stage of the cooling phase. The time when Qx starts to increase dramatically can be identified from Fig. 11, and the isotherms and streamlines at the corresponding time suggest that this is when cold thermal plumes start to descend from the surface. 4.4.1. The phase lag of flow response The effect of period length on the flow response and phase lag is also reflected by the average volumetric horizontal flow rate Q in Fig. 12. The time series of the calculated Q starts from the very beginning of the simulation when the flow domain is stationary and isothermal, and thus the start-up effect can be observed. It is clear that the start-up effect persists in fewer cycles as P increases. The general trend is that Q becomes more in phase with the surface temperature s as P increases. For the short period cases (Fig. 12b– d), the value of Q never approaches zero. This is consistent with the fact that cooling-induced flow dominates over the entire cycle with no reversal of the large-scale flow (the left column of Fig. 9). In these cases, Q starts to decrease at a certain stage of the heating phase (Fig. 12b–e) since heating acts mainly to weaken the cooling-induced flow, and the reduction of Q proceeds into the next cooling phase. In contrast, for the very long period cases (Fig. 12i–k), the variation of Q is much more in phase with the sinusoidally varying surface temperature. The value of Q is much higher during the cooling phase than the heating phase, owing to the large scale down-slope flow formed by the cold plunging ther-

mal plumes during the cooling phase. If there were no inertia and the associated phase lag, the flow would respond instantaneously to the surface temperature s, reaching the minimum at the start of a cycle. Nevertheless, the minimum Q occurs at the early stage of the heating phase (Fig. 12i–k). The phase of this minimum Q represents the phase lag of the flow response, and it decreases with increasing P, consistent with the scaling prediction (34). Since in the convective flow regime, the value of Q is dominated by the flow in the convective region, scale (34) which describes the phase lag of u in the convective region applies for the overall Q as well. Substituting x = 1/A into (34) gives an estimation of the phase lag of Q. For the long period cases, scale (34) predicts a phase lag proportional to P
2/3 and Ra
1/3. This scale is confirmed by the numerical results as shown by the clear linear correlation, and a scaling factor (constant) of 2.9 is found by the least square curve fitting (Fig. 13). It should be noted that scale (34) applies only to the long period cases. According to the long period criterion (22) (P > P0, P0  (p/2)1/2x1/2Ra
1/2A
1/2), as Ra increases, the critical period P0 decreases, and thus a wider range of P are selected (Fig. 13). Similarly, according to (40) (Ra > p=2P
2 A
1 x), as P increases, the critical Ra decreases, and thus a wider range of Ra are selected for verification (Fig. 13). 4.5. Validation of quasi-steady-state scales for distinct flow regions 4.5.1. Conduction-dominated region For the conduction-dominated region, scales (29) and (32) provide the scales for the quasi-steady flow velocity and the phase lag, respectively. Being of the opposite sign, the maximum velocities in the upper and lower layers are averaged as (umax
umin)/2, where umax is positive and umin is negative. This average value are obtained from simulation results for different horizontal positions x and different periods P (Fig. 14a). As indicated by the cosine variation in scale (29) (u  16p
1 P
1 RaA5 x4 cosð2pt=PÞ), the velocity would reach a maximum at the start of a cycle if there were no phase lag. However, this maximum is evidently delayed and the

Fig. 11. The horizontal flow rate Qx(x) over a typical thermal forcing cycle (the 9th cycle) under different forcing periods for Ra = 3  106.

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

623

Fig. 12. Time series of average flow rate Q for different periods P at Ra = 3  106. Plotted are (a) time series of surface temperature. (b–k) time series of Q for different periods.

Fig. 13. Verification of the phase delay h specified by scale (34) by the simulation results of the average flow rate Q.

phase lag increases as x increases or as P decreases (Fig. 14a), consistent with the scaling prediction (32). The phase lag obtained from simulation is plotted against the scaling result in Fig. 14b, where a clear linear correlation is identified, verifying the dependency of the phase lag on x and P. For the verification of velocity scale in the conductive region, the time series of the maximum horizontal velocity (umax
umin)/2 are obtained over the local depth for different P, Ra and x (Fig. 15a). After normalizing by scale (29), the various time series almost collapse onto a single curve (Fig. 15b), verifying the dependency of the

magnitude of umax on P, Ra and x. Since the focus here is on the magnitude of the velocity, the effect of phase lag is minimized by selecting relatively large P and small x, as scale (32) (h  2pA2 x2 P
1 Pr
1 ) suggests. 4.5.2. Convective region The present scaling is based on the simplified heat flux specified in (13). This simplification is based on the assumption of negligible convection, which is satisfied when the flow is about to reverse and the flow velocity is low. For the long period cases, at the end of the

624

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

Fig. 14. Verification of the scale of phase lag h in the conductive region, Ra ¼ 104 . (a) Time series of 0.5(umax
umin), the vertical dashed lines mark the time when maximum occurs (b) h at different x and P obtained from simulation versus scaling, the x positions are extracted at equal intervals within the specified range.

heating phase (t = nP + 0.5P), the flow velocity is low, as shown in Fig. 16 (the upper panel), since the flow is about to reverse from heating-induced flow to cooling-induced flow. The flow remains stable at relatively large offshore distances (Fig. 9, right column). After normalizing the time series of the local maximum velocity umax by scale (20), the values of umax for different Ra, P, and x converge at the end of the heating phase (t = nP + 0.5P) (Fig. 16, the lower panel), verifying the dependency of u on Ra, P, and x. These dependencies are further verified by a clear linear correlation between the simulation results of flow velocity at t = nP + 0.5P and the those given by the velocity scale in (20) (Fig. 17). A clear phase lag is observed as the flow velocity continues to decrease at the start of a new cycle (t = nP), reaching the minimum at the initial stage of the heating phase (Fig. 16). The phases corresponding to this minimum value are obtained from the simulation results and are plotted against those predicted by the phase lag scale specified in (34) (Fig. 17b), and the linear correlation confirms the dependency of the phase lag on Ra, P and x. 5. Discussion In the laboratory experiment of Bednarz et al. [37], a phase lag of the flow response was observed as the thermal forcing switches from cooling to heating. The estimated time delay of the flow response was 11% of the forcing period. The parameter setting in that experiment corresponds to a Rayleigh number of 1.2  105

based on the present Ra definition, and a non-dimensional period of P = 0.534, which falls into the long period scenario of P > P0 ðP 0  ðp=2Þ1=2 x1=2 Ra
1=2 A
1=2 Þ as P0  0.036. Although the scaling factor of 2.9 obtained in the present study may not be relevant to the experiment due to the different boundary conditions at the far end (an open boundary here versus a rigid endwall in the experiment), the order of magnitude estimated by scaling can be compared with the experimental result. Substituting the experimental parameters into scale (33) results in a time delay of 8.3% of the forcing period, which is of the same order of magnitude as the reported experimental result. In contrast to the current numerical simulations which adopt many different periods, covering both the short and the long period scenarios, the previous experiment [36] and numerical simulations [37] for the periodic thermal forcing adopt only a single period which belongs to the long period scenario. Since the boundary conditions and the exact value of the controlling parameters of the present study are different from the previous ones, no direct comparison can be made. Nevertheless, the present results of the horizontal exchange flow rate Q for the long period scenario agree well with the previous experiments and simulations [37,38], showing a much higher flow rate during the cooling phase than the heating phase. Based on the present scaling, some estimation of the field situations can be made. The field parameters are adopted as follows: a diurnal temperature difference DT of 5 °C at the water surface, a diurnal period P of 24 h, a bottom slope A of 0.05 which falls into

Fig. 15. Verification of the dependency of umax on x, P and Ra, (a) Time series of 0:5ðumax
umin Þ at various P, Ra and x obtained from simulation. (b) Simulation results normalized by the amplitude of scaling (29).

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

625

Fig. 16. Time series of the horizontal velocity u over two typical cycles (the upper panel). Time series of u normalized by scale (20) (the lower panel). The normalized results of different Ra, P and x converged at the time of t = nP + 0.5P.

Fig. 17. Verification of the flow velocity and phase lag in the convective region. (a) The local maximum flow velocities umax at the end of the heating phase (t = nP + 0.5P) obtained from simulation results versus those given by the velocity scale (20). (b) The phase lag of flow reversal obtained from simulation results versus the scale of phase lag specified in (34). The interval of x is 0.5.

scenario (a) (A > Ra-1/2  0.04), and a maximum water depth h of c 10 m. The above parameters and the molecular properties of water give a Ra value of approximately 6.8  1013. For this set of parameters, the Rayleigh number crosses with both f1(x) and f2(x) and is larger than the maximum of f3(x) (4  1010), and therefore in real field cases, the flow falls into the unstable convective regime. Whether a quasi-steady state is reached or not depends on the comparison of the quasi-steady time with the period of the diurnal cycle. For the conductive subregion, according to scale (27) (t d  A2 x2 ), the maximum quasi-steady state time td is reached at the maximum x, which is scaled by x1 (scale (25)). Substituting the scale of x1 into scale (27), the scale of the maximum quasisteady state time within the conductive region becomes,

t d  ðp=16Þ2=5 Ra
2=5 A
4=5 P2=5

ð45Þ

Substituting the assumed field parameters into (45), the maximum dimensional value of t d is 0.09 h, much shorter than the diurnal cycle. Therefore, a quasi-steady state is reached for the nearshore conductive region. Similarly, for the convective subregion, the quasi-steady time specified by scale (21) reaches its maximum at the maximum offshore distance of x = 1/A, resulting in a maximum dimensional value of 1.0 h, also much shorter than the diurnal cycle, and thus a quasi-steady state is also expected to be reached for the convective subregion. Therefore, the diurnal period under the field conditions is classified as long period, for which the present scaling for the flow velocity and phase lag can be applied to. The flow velocity, specified by scale (20), reaches the maximum at x = 1/A. Substituting x = 1/A into scale (20) gives an estimation of the maximum flow velocity of 5.5 cm/s. This magnitude of veloc-

626

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

ity is significant under calm weather conditions and is of the same order of magnitude as those measured in the field [5]. Therefore, natural convection generated by the diurnal temperature variation at the water surface cannot be neglected in near-shore waters under calm weather conditions when wind speed is small. The phase lag of the flow response specified by scale (34) has been verified by the numerical results. Converting the phase lag to time delay by multiplying P/2p and using the scaling factor of 2.9 obtained from the previous curve fitting (Fig. 13) result in a time delay of about 3 h. It is worth noting that before the temperature s at the water surface switches sign, the surface heat flux (@ s=@y) has already switched sign a quarter of a cycle (i.e. 6 h) earlier. Therefore, the total phase lag between the reversal of the flow and the reversal of the surface heat flux is more than 6 h, which is consistent with field observations [5]. As mentioned in Mao et al. [19,31], many other geophysical mechanisms are generating flows simultaneously with the periodic thermal siphon, including tides, wind shear stress, stokes drift of surface waves, geostrophic pressure gradients, and salinity gradients, etc. Further, for thermally driven flows, there are other forms of thermal forcing, such as the absorption of solar radiation [19], heat flux from the sediment [13] and the shorelines [14], and a spatially non-uniform flux along one horizontal boundary which generates horizontal convection that drives the large scale ocean circulation [44,45]. A study of each of these mechanisms in isolation allows an estimation of the significance of different mechanisms and a more comprehensive understanding of the overall response of the natural water body. In most cases, wind exerts a strong influence on the surface flow through shear stress and waves. For lakes and reservoirs, wind is an important factor that should not be neglected for estimating surface current. The response of surface current to wind varies with latitude, fetch and duration conditions [46]. Wind generated surface current is composed of two ingredients: the Ekman-type current induced by the tangential stress of wind, and the Stokes drift generated by the nonlinearity of surface waves. The former typically generates a surface current 1–2.5% of the wind speed depending on the latitude [47], while the latter generates a current 0.8–1.6% of the wind speed depending on the fetch if waves are not duration limited, i.e. given sufficient time to develop [46]. The former is 45° to right/left of the wind vector at the water surface in the north/south hemisphere, while the latter is usually in the same direction as the wind. Given the angle difference, a relative broad and therefore inclusive range of 1–3.8% can be estimated for the ratio of the resultant surface current to the wind speed, with the lower end for flows at high latitudes and no development of surface waves. For a moderate wind speed of 3–6 m/s, a surface current of 3–23 cm/s will be generated. The present scaling analysis has estimated a thermally driven current of the order of 5 cm/s as discussed above. Therefore, the thermally driven flow is of a similar order to the lower end of the current generated by moderate wind, and thus is ineligible under calm weather conditions with a smooth water surface. For the long period scenario, the flow rate in the heating phase is shown to be much weaker than the cooling phase (Figs. 11 and 12). This agrees with field observations that the thermal flow during the daytime is easily overwhelmed by wind while the nighttime cooling induced flow is more persistent and stronger [10].

6. Conclusions A hybrid of semi-analytical solutions coupled with scaling analysis and numerical simulations is adopted to investigate the nearshore natural convection induced by a periodic surface temperature. By focusing on the specific time when the convective

flow is about to reverse, a set of scales for flow properties, the phase lag and the dividing positions between various flow subregions are derived. The thermal forcing period can be classified as long if P > P0, where P0  ðp=2Þ1=2 x1=2 Ra
1=2 A
1=2 is the quasisteady time. Correspondingly, a critical Rayleigh number function f3(x) has been derived: f3(x)  p/2P
2A
1x. If Ra > f3(x), the flow falls into the long period scenario for which a quasi-steady state is reached before the thermal forcing switches sign. The present scaling and numerical verification focus on the long period scenario, which is relevant to field situations. Two other critical functions f1(x) and f2(x) have been derived to identify the distinctness and stability of the TBL: f1(x)  p/16A
7Px
5 and f2(x)  p/16A
1Px
5Rac3. The TBL is distinct if Ra > f1(x), and it is unstable if Ra > f2(x). A comparison between f1(x) and f2(x) reveals two possible scenarios: (a) A > Rac-1/2, for which f1(x) < f2(x) for all x. (b) A < Rac-1/2, for which f1(x) > f2(x) for all x. The present study focuses on the scenario of relatively large slope, i.e. A > Ra
1=2 , which can be further classified into three flow regimes c by comparing Ra with the minimum of f1(x) and f2(x). (i) Ra < p/16A
2P, the entire flow domain is dominated by conduction (ii) For p/16A
2P < Ra < p/16PA4Rac3, the dominant heat transfer mode switches from conduction to stable convection as the offshore distance x increases. (iii) For Ra > p/16PA4Rac3, the dominant heat transfer mode changes from conduction to stable convection and finally to unstable convection as x increases The dividing position between conductive and convective subregions is at x1  ðp=16Þ1=5 Ra
1=5 A
7=5 P1=5 , and that between stable and unstable subregions is at x2  ðp=16Þ1=5 Ra
1=5 ðP=AÞ1=5 Ra3=5 c : Both x1 and x2 move toward the shore as Ra increases or as P decreases. In the conductive subregion, the scale of phase lag h of the flow response is h  2pA2 x2 =ðPPr Þ, and the flow velocity is proportional to 16p
1 RaA5 x4 P
1 . For the convection-dominated subre-

gion, the phase lag h is scaled as h  2
1=3 p4=3 x1=3 Ra
1=3 A
1=3 P
2=3 ,

and the flow velocity is scaled as u  ð16=pÞ1=3 x2=3 Ra1=3 A1=3 P
1=3 . The above scaling results have been verified by the results of numerical simulations. Results of numerical simulation show that the flow response varies significantly with the thermal forcing period P. For short P, the large-scale flow never reverses from the cooling-induced down-slope flow, with flow instability present over the entire cycle. As P increases, a large-scale flow reversal begins to occur during the heating phase, and the phase lag deceases with increasing P. The cooling-induced flow is much more intense than the heating-induced flow owing to the occurrence of flow instability during the cooling phase, with the cold thermal plumes plunging down from the surface and merging together to form a cool gravity current down the bottom slope. While the previous scaling analyses on nearshore natural convection are mainly concerned with constant thermal forcing, a periodic thermal forcing is more appropriate to represent field situations. It is the periodic variation of the thermal forcing that brings in many new aspects to the scaling analysis, making it more interesting and challenging. Some new features pertinent to the periodic thermal forcing conditions are: (1) The phase lag between the flow response and the thermal forcing; (2) The short and the long period scenarios depending on whether the thermal forcing period is sufficiently long for a quasi-steady state to be reached; and (3) The dependency of all the scaling results, including the critical functions, the dividing positions of different subregions, the thermal and momentum quantities of the flow, on the thermal forcing period. All these new aspects are not relevant to the cases

Y. Mao et al. / International Journal of Heat and Mass Transfer 134 (2019) 610–627

with constant or ramped heating or cooling that have been studied previously. Furthermore, for the periodic thermal forcing cases, much richer flow phenomena and possibilities may exist. For example, the flow responses in the short-period and the longperiod scenarios are very different (Fig. 9). The present scaling analysis reveals that the real field situation belongs to the longperiod scenario for which a quasi-steady state can be reached before the thermal forcing switches sign. The hybrid analytical approach adopted here may be applied to analyzing other similar periodic problems. Declarations of interest None. Acknowledgements This research was supported by National Natural Science Foundation of China (grant number 11472250, 41630317). References [1] H.B. Fisher, J.E. List, C.R. Koh, J. Imberger, N.H. Brooks, Mixing in Inland and Coastal Waters, Academic, 1979. [2] J. Imberger, J.C. Patterson, Physical limnology, Adv. Appl. Mech. 27 (1990) 303– 475. [3] J. Imberger, Flux paths in a stratified lake: a review, in: J. Imberger (Ed.), Physical Processes in Lakes and Oceans, Washington, D.C., American Geophysical Union, vol. 54, 1998, pp. 1–18. [4] E.E. Adams, S.A. Wells, Field measurements on side arms of lake, J. Hydraul. Eng. ASLE 110 (1984) 773–793. [5] S.G. Monismith, J. Imberger, M.L. Morison, Convective motions in the sidearm of a small reservoir, Limnol. Oceanogr. 35 (1990) 1676–1702. [6] S.G. Monismith, A. Genin, M.A. Reidenbach, G. Yahel, J.R. Koseff, Thermally driven exchanges between a coral reef and the adjoining ocean, J. Phys. Oceanogr. 36 (2006) 1332–1347. [7] W.F. James, J.W. Barko, Estimation of phosphorus exchange between littoral and pelagic zones during nighttime convective circulation, Limnol. Oceanogr. 36 (1991) 179–187. [8] W.F. James, J.W. Barko, H.L. Eakin, Convective water exchanges during differential cooling and heating: implications for dissolved constituent transport, Hydrobiologia 294 (1994) 167–176. [9] G. Symonds, R. Gardiner-Garden, Costal density currents forced by cooling events, Cont. Shelf Res. 14 (1994) 143–157. [10] H. Niemann, C. Richter, H.M. Jonkers, M.I. Badran, Red sea gravity currents cascade near-reef phytoplankton to the twilight zone, Mar. Ecol.-Prog. Ser. 269 (2004) 91–99. [11] L. Molina, G. Pawlak, J.R. Wells, S.G. Monismith, M.A. Merrifield, Diurnal crossshore thermal exchange on a tropical forereef, J. Geophys. Res. Oceans 119 (2014) 6101–6120. [12] S.G. Monismith, Hydrodynamics of coral reefs, Annu. Rev. Fluid Mech. 39 (2007) 37–55. [13] G.B. Kirillin, M. Leppäranta, A. Terzhevik, N. Granin, J. Bernhardt, C. Engelhardt, T. Efremova, S. Golosov, N. Palshin, P. Sherstyankin, G. Zdorovennova, R. Zdorovennov, Physics of seasonally ice-covered lakes: A review, Aquat. Sci. 74 (2012) 659–682. [14] G.B. Kirillin, A.L. Forrest, K.E. Graves, A. Fischer, C. Engelhardt, B.E. Laval, Axisymmetric circulation driven by marginal heating in ice-covered lakes, Geophys. Res. Lett. 42 (2015) 2893–2900. [15] K. Salonen, M. Pulkkanen, P. Salmi, R.W. Griffiths, Interannual variability of circulation under spring ice in a boreal lake, Limnol. Oceanogr. 59 (2014) 2121–2132. [16] B. Yang, J. Young, L. Brown, M. Wells, High-frequency observations of temperature and dissolved oxygen reveal under-ice convection in a large lake, Geophys. Res. Lett. 44 (2017) 12218–12226. [17] D.E. Farrow, J.C. Patterson, On the stability of the near shore waters of a lake when subject to solar heating, Int. J. Heat Mass Transf. 36 (1993) 89–100.

627

[18] C. Lei, J.C. Patterson, Unsteady natural convection in a triangular enclosure induced by absorption of radiation, J. Fluid Mech. 460 (2002) 181–209. [19] Y. Mao, C. Lei, J.C. Patterson, Unsteady natural convection in a triangular enclosure induced by absorption of radiation – a revisit by improved scaling analysis, J. Fluid Mech. 622 (2009) 75–102. [20] C. Lei, J.C. Patterson, A direct stability analysis of a radiation-induced natural convection boundary layer in a shallow wedge, J. Fluid Mech. 480 (2003) 161– 184. [21] Y. Mao, C. Lei, J.C. Patterson, Characteristics of instability of radiation-induced natural convection in shallow littoral waters, Intl. J. Therm. Sci. 49 (2010) 170– 181. [22] Y. Mao, C. Lei, J.C. Patterson, Unsteady nearshore natural convection induced by constant isothermal surface heating, J. Fluid Mech. 707 (2012) 342–368. [23] G.M. Horsch, H.G. Stefan, Convective circulation in littoral water due to surface cooling, Limnol. Oceanogr. 33 (1988) 1068–1083. [24] H.G. Stefan, G.M. Horsch, J.W. Barko, A model for the estimation of convective exchange in the littoral region of a shallow lake during cooling, Hyrobiologia 174 (1989) 225–234. [25] G.M. Horsch, H.G. Stefan, S. Gavali, Numerical simulation of cooling-induced convective currents on a littoral slope, Int. J. Numer. Methods Fluids 19 (1994) 105–134. [26] J.J. Sturman, C.E. Oldham, G.N. Ivey, Steady convective exchange flows down slopes, Aquatic Sci. 61 (1999) 260–278. [27] M.G. Wells, B. Sherman, Stratification produced by surface cooling in lakes with significant shallow regions, Limnol. Oceanogr. 46 (2001) 1747–1759. [28] C. Lei, J.C. Patterson, Unsteady natural convection in a triangular enclosure induced by surface cooling, Int. J. Heat Fluid Fl. 26 (2005) 307–321. [29] T.P. Bednarz, C. Lei, J.C. Patterson, An experimental study of unsteady natural convection in a reservoir model cooled from the water surface, Exp. Therm Fluid Sci. 32 (2008) 844–856. [30] T.P. Bednarz, C. Lei, J.C. Patterson, A numerical study of unsteady natural convection induced by iso-flux surface cooling in a reservoir model, Int. J. Heat Mass Transf. 52 (2009) 56–66. [31] Y. Mao, C. Lei, J.C. Patterson, Unsteady near-shore natural convection induced by surface cooling, J. Fluid Mech. 642 (2010) 213–233. [32] K.A. Dittko, M.P. Kirkpatrick, S.W. Armfield, Three-dimensional simulation of natural convection in a reservoir sidearm, Phys. Fluids 25 (2013) 025105. [33] K.A. Dittko, M.P. Kirkpatrick, S.W. Armfield, Large eddy simulation of complex sidearms subject to solar radiation and surface cooling, Water Res. 47 (2013) 4918–4927. [34] D.E. Farrow, J.C. Patterson, On the response of a reservoir sidearm to diurnal heating and cooling, J. Fluid Mech. 246 (1993) 143–161. [35] D.E. Farrow, Periodically forced natural convection over slowly varying topography, J. Fluid Mech. 508 (2004) 1–21. [36] C. Lei, J.C. Patterson, Natural convection induced by diurnal heating and cooling in a reservoir with slowly varying topography, JSME Int. J. B-Fluid T. 49 (2006) 605–615. [37] T.P. Bednarz, C. Lei, J.C. Patterson, An experimental study of unsteady natural convection in a reservoir model subject to periodic thermal forcing using combined PIV and PIT techniques, Exp. Fluids 47 (2009) 107–117. [38] S. Chen, C. Lei, J.C. Patterson, Numerical modelling of a concurrent PIT/PIV experiment with TLC particles in a reservoir model subject to periodic thermal forcing, Numer. Heat Trans. A-Appl. 66 (2014) 64–88. [39] D.E. Farrow, Periodically driven circulation near the shore of a lake, Environ. Fluid Mech. 13 (2013) 243–255. [40] D.E. Farrow, A numerical model of periodically forced circulation near the shore of a lake, Environ. Fluid Mech. 16 (2016) 983–995. [41] Y. Mao, C. Lei, J.C. Patterson, Unsteady natural convection in a reservoir sidearm induced by time-varying isothermal surface heating, Int. J. Therm. Sci. 71 (2013) 61–73. [42] P.G. Drazin, W.H. Reid, Hydrodynamic Stability, Cambridge University Press, 1981. [43] J.C. Patterson, C. Lei, S.W. Armfield, W. Lin, Scaling of unsteady natural convection boundary layers with a non-instantaneous initiation, Int. J. Therm. Sci. 48 (2008) 1843–1852. [44] G.O. Hughes, R.W. Griffiths, Horizontal convection, Annu. Rev. Fluid Mech. 40 (2008) 185–208. [45] O. Shishkina, S. Grossmann, D. Lohse, Heat and momentum transport scalings in horizontal convection, Geophys. Res. Lett. 43 (2016) 1219–1225. [46] Y. Mao, M.L. Heron, The influence of fetch on response of surface currents to wind studied by HF ocean surface radar, J. Phys. Oceanogr. 38 (2008) 1107– 1121. [47] R.H. Stewart, Response of the upper ocean to winds, Introduction to Physical Oceanography, Texas A&M University, 2005, pp. 133–138.