2.1 Fundamentals of TS

Thermal siphons are a fascinating convective dynamics phenomenon. Their emergence hinges on the presence of a destabilising density gradient, which creates a pressure gradient that ultimately drives a convective flow. If we consider the fluid to be freshwater, its density and spatial gradients are determined by temperature. The action of gravity on horizontal density gradients creates buoyancy anomalies, which become a source of available potential energy that can be transformed into kinetic energy – i.e. motion (Reference Winters, Lombard, Riley and D'AsaroWinters et al. 1995; Reference Hughes, Gayen and GriffithsHughes, Gayen & Griffiths 2013). Although other fluid phenomena share the above driving mechanisms, including the canonical problem of horizontal convection (Reference RossbyRossby 1965; Reference Sturman, Ivey and TaylorSturman, Ivey & Taylor 1996; Reference Hughes and GriffithsHughes & Griffiths 2008; Reference Winters and YoungWinters & Young 2009; Reference Gayen, Griffiths, Hughes and SaenzGayen et al. 2013; Reference Griffiths, Hughes and GayenGriffiths, Hughes & Gayen 2013; Reference Amber and O'DonovanAmber & O'Donovan 2018) and more complex problems such as the stratified horizontal convection (Reference Noto, Ulloa, Yanagisawa and TasakaNoto et al. 2023), these phenomena result from a direct gradient of heating or cooling of the bottom or the top surface of the fluid. The case of TS is different: differential cooling or heating refers to the interplay between spatially uniform atmospheric forcing and varying bathymetry that leads to spatially varying rates of temperature change (cooling or heating).

As an analogue model for aquatic systems, we consider freshwater bodies with a horizontal surface with a normal vector parallel to the gravity field and a sloping bottom boundary, such that the systems have distinctive shallow and deep regions (figure 3). A uniform cooling of the surface, or internal heating of the body, induces a horizontal density gradient, eventually initiating an overturning flow between the shallow and deep waters. This horizontal density gradient occurs because the net heating or cooling rate in the fluid column depends on the local depth: shallower waters cool or warm faster than deep waters. Therefore, a critical element for TS development is the presence of a sloping bottom boundary, a fundamental geometric characteristic of natural aquatic systems. Table 1 provides a comprehensive overview of the key variables and parameters essential for investigating and characterising thermal siphons, along with their respective magnitudes.

Figure 3. (a,b) Conceptual model adopted by Reference Horsch and StefanHorsch & Stefan (1988) and Reference Sturman, Oldham and IveySturman et al. (1999) for the study of convective circulation due to a destabilising surface buoyancy flux $B_{0}$ in littoral waters of characteristic slope $s = h_{m}/\ell _{s}$, the ratio of the interior surface mixed layer $h_{m}$ to the horizontal length of the nearshore sloping bottom $\ell _{s}$. (c,d) Conceptual model adopted by Reference Wells and ShermanWells & Sherman (2001) and Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. (2022) for littoral regions of water basins. The system considers a shallow plateau (P) of characteristic depth $d$ and horizontal extent $\ell _{p}$ joined through a sloping (S) bottom region of horizontal extent $\ell _{s}$ and characteristic slope $\bar {s}$ to the interior (I) stratified basin. The reference system is positive upwards with the origin at the base of the water body. The horizontal distance $\ell _{ml}=\ell _{p}+\ell _{s}$ characterises the littoral region. The initial temperature distribution is characterised by a two-layer stratification with a surface well-mixed layer of thickness $h_{m}$, a smooth metalimnion region of length $\delta _{m}$ and a maximum depth $D$. The free surface has an inhomogeneous Neumann boundary condition that models a uniform surface heat loss that sets a uniform destabilising surface buoyancy flux $B_{0}$. For both models, the solid bottom is considered adiabatic and no-slip. Colours provide a conceptual temperature distribution, with colder temperatures depicted with darker blue colours.

Table 1. List of relevant variables and parameters to characterise TS.

2.2 Fluid properties

We focus our review on shallow freshwater systems, explicitly ignoring the effects of salt and its gradients throughout the water body. We also assume that variations in pressure are negligible in littoral regions where TS occur (e.g. compressibility effects are neglected). In this scenario, the equation of state of water, linking the fluid temperature $T$ to its density $\rho$, can be characterised by a polynomial function. A reference model is provided by (Reference Chen and MilleroChen & Millero 1986):

(2.1)\begin{align} \rho(T) & = 999.8395 + 6.7914\times 10^{{-}2}T-9.0894\times 10^{{-}3}T^{2} + 1.0171\times 10^{{-}4}T^{3} \nonumber\\ &\quad -1.2846\times 10^{{-}6}T^{4} + 1.1592\times 10^{{-}8}T^{5} -5.0125\times 10^{{-}11}T^{6}\,({\rm kg}\,{\rm m}^{{-}3}). \end{align}

An important fluid property for studying TS is the thermal expansivity of water:

(2.2)\begin{equation} \alpha = {-}\frac{1}{\rho}\left(\frac{\partial \rho}{\partial T}\right)(^{{\circ}}{\rm C}^{{-}1}). \end{equation}

This quantity can be readily determined from (2.1). The thermal expansivity $\alpha$ determines the magnitude of the density anomalies due to temperature differences. Lateral temperature gradients $\partial T/\partial x$ originating from differential heating or cooling create lateral density gradients $\partial \rho /\partial x=-\alpha \rho \partial T/\partial x$ that drive TS, where $x$ is the cross-shore coordinate increasing with distance from the shore (figure 3). The maximal water density ($\alpha =0$) is reached around $T_{md}\approx 3.98\,^{\circ }{\rm C}$. For systems with $T < T_{md}$ ($\alpha <0$), net heating leads to an increase in $\rho$ and thereby an increase in density in the littoral region compared with the pelagic region due to differential heating ($\partial \rho /\partial x<0$). The same process occurs for a net cooling in the case $T > T_{md}$ ($\alpha >0$), leading to $\partial \rho /\partial x<0$ by differential cooling. Reversed lateral density gradients $\partial \rho /\partial x>0$ are created by differential cooling if $T < T_{md}$ and by differential heating if $T > T_{md}$.

Additional relevant fluid properties include the molecular kinematic viscosity $\nu$ and thermal diffusivity $\kappa$. Despite being broadly considered constant, both parameters are temperature dependent. The relevant non-dimensional number that encapsulates the information of both molecular parameters is the Prandtl number, $Pr = \nu /\kappa$. In the case of surface water bodies with temperature around 20 $^{\circ }$C, $\nu \approx 10^{-6}$ ${\rm m}^{2}\,{\rm s}^{-1}$ and $\kappa \approx 1.43\times 10^{-7}$ ${\rm m}^{2}\,{\rm s}^{-1}$, $Pr \approx 7$. However, for freshwater bodies with surface temperatures ranging between $0\,^{\circ }$C and $30\,^{\circ }$C, $Pr$ varies in the range $5< Pr<14$.

2.3 Boundary conditions and external forcing

2.3.2 Boundary conditions

Conceptual models assume a no-slip bottom boundary and generally ignore any exchange at the sediment–water interface. In the case of the top boundary, we consider the open air–water interface free-to-slip and the ice–water interface of ice-covered waterbodies to be a no-slip surface. In this review, as in most of the literature on TS, we assume no heat flux at the sediment–water interface. Readers interested in the heat exchange at the sediment–water interface (and its impacts) are referred to the works by Reference Hondzo and StefanHondzo & Stefan (1993), Reference Fang and StefanFang & Stefan (1996), Reference Zdorovennova, Terzhevik, Palshin, Efremova, Bogdanov and ZdorovennovZdorovennova et al. (2021), Reference MacIntyre, Cortés and SadroMacIntyre, Cortés & Sadro (2018), Reference Cortés and MacIntyreCortés & MacIntyre (2020) and Reference PergaPerga et al. (2023). To resolve the evolution of heat in the aquatic system, we need to characterise the surface heat exchange and the internal heating caused by solar radiation throughout the water column. Details are provided in Appendix A and we only report the main results and final equations below.

The surface heat flux $H_{Q,0}$ (${\rm W}\,{\rm m}^{-2}$) results from the combined action of four processes exchanging heat between the atmosphere and the water body: sensible heat exchange, latent heat exchange, incoming and outgoing long-wave radiation (Reference Schmid and ReadSchmid & Read 2022).

According to the convention in field studies, the heat content of the water body increases if $H_{net,0}<0$ (net heating) and decreases if $H_{net,0}>0$ (net cooling). The surface heat flux changes the potential energy of the water column, which is expressed by a surface buoyancy flux per unit mass $B_{0} = \alpha g H_{Q,0}/(\rho _{0} C_p)$ (${\rm W}\,{\rm kg}^{-1}$) (Reference Soloviev and LukasSoloviev & Lukas 2014), where $g$ is the gravitational acceleration (${\rm m}\,{\rm s}^{-2}$), $C_p$ is the specific heat of water (${\rm J}\,^{\circ }{\rm C}^{-1}\,{\rm kg}^{-1}$) and $\rho _{0}$ is a reference density. Surface heat exchange stabilises (stratifies) the water column when $B_{0} < 0$ and destabilises it (promoting convection) when $B_{0} > 0$. For $T > T_{md}$ ($\alpha > 0$), surface heating stabilises and surface cooling destabilises the water column. The process reverses for $T < T_{md}$ and $\alpha < 0$. Considering the same approximation made in (A4), the net buoyancy flux is

(2.3)\begin{equation} B_{net,0} = B_{0} + B_{SW,0} = \frac{g \alpha}{\rho_{0} C_p}H_{Q,0} + \frac{g \alpha}{\rho_{0} C_p}H_{SW,0}, \end{equation}

where $H_{SW,0}$ is the solar radiation reaching the air–water interface and $B_{SW,0}$ is the solar radiative buoyancy flux at the surface.

3.1 Steady differential heating above $T_{md}$ (stabilising forcing)

Differential heating above $T_{md}$ corresponds to the case of stabilising surface buoyancy flux over a sloping boundary. A series of numerical investigations explored the influence of uniform surface heating and internal heating on freshwater bodies with varying depths, employing a basic geometry to model the littoral environment – a triangular domain connected to a flat interior basin, as shown in figure 3(a) (Reference Lei and PattersonLei & Patterson 2002, Reference Lei and Patterson2003; Reference Mao, Lei and PattersonMao et al. 2009; Reference Mao, Lei and PattersonMao, Lei & Patterson 2012; Reference Yu, Patterson and LeiYu, Patterson & Lei 2015). The equation of state for water was taken as a linear function between temperature and density, with temperatures above $T_{md}$. These numerical experiments consistently showed the emergence of TS propelled by differential heating. In this scenario, shallow waters, warming more rapidly than their deeper counterparts, become less dense, initiating a two-layer exchange flow, directing upslope deep waters toward the littoral and offshore shallow waters across the surface, as illustrated in figure 4(a). Thermal siphons exhibit three distinct regimes across the sloping region. From the very shallow to deeper waters, the system shows a ‘conductive’ region I with vertical isotherms, a ‘stable convective’ region II characterised by a two-layer exchange flow, and an ‘unstable convective’ region III where a surface flow coexists with a bottom convective layer. The spatial boundaries of these regions and regimes were determined by the slope of the wedge-like basin and the global Rayleigh number of the system (Reference Mao, Lei and PattersonMao et al. 2009). As mentioned in the introduction, this case is rarely observed in the field and is not further developed in this review. It is however important to remind the reader that the dynamics induced by stabilising surface buoyancy flux over a sloping boundary are not analogous to those induced by a destabilising surface buoyancy flux. For instance, the downslope flow resulting from the destabilising forcing is affected by the development of a convective mixed layer (see § 6.3) that is obviously not present in the stabilising forcing case.

Figure 4. Schematic of isotherm distributions for different convective flow sub-regions in nearshore waters with $T>T_{md}$ due to (a) differential heating and (b) differential cooling. Adapted from Reference Mao, Lei and PattersonMao, Lei & Patterson (2009, Reference Mao, Lei and Patterson2010). Note that Reference Mao, Lei and PattersonMao et al. (2009) reported rising thermal plumes at the bottom of the region (III). These instabilities can be attributed to the model's bottom boundary condition, where reflected heat leads to localised excess heating.

3.2 Steady differential cooling above $T_{md}$ (destabilising forcing)

Differential cooling above $T_{md}$ corresponds to the case of destabilising surface buoyancy flux over a sloping boundary. In a wedge-shaped water basin with a uniform slope, Reference Horsch and StefanHorsch & Stefan (1988) and Reference Horsch, Stefan and GavaliHorsch, Stefan & Gavali (1994) conducted pioneering experimental and numerical investigations, unravelling the physical processes behind the formation of TS by differential cooling (see schematic in figure 4b). As thermal plumes sinking from the surface interact with the sloping bottom, they contribute to the creation of a downslope gravity current in the shallowest region. By mass conservation, this thermally induced gravity current is compensated by a surface current toward the shore, thus establishing an overturning circulation characterised by a two-layer exchange flow: the bottom layer transports colder, denser waters downslope (offshore) and the surface layer conveys warmer, lighter waters towards the shore. Numerical studies exploring the consequences of surface cooling in wedge-like domains identified three regimes across the sloping region, as in the case of differential heating (Reference Lei and PattersonLei & Patterson 2005; Reference Mao, Lei and PattersonMao et al. 2010; Reference Papaioannou and PrinosPapaioannou & Prinos 2023).

Reference Wells and ShermanWells & Sherman (2001) explored TS via laboratory experiments using a distinct basin design, akin to the schematic in figure 3(c). The basin comprised a flat shelf connected to a deeper area via a sloping bottom. Before the TS formation, the authors noted an active phase of vertical mixing, resulting in a significant horizontal temperature difference between shallow and deep regions. They visualised the circulation using dye, as shown in figure 5, and noticed that the dye injected in the surface region was transported downward into the convective mixing layer by thermal plumes, while the dye injected on the shallow plateau was transported as a density current ultimately filling in the deep water. To estimate gravity current initiation, Reference Wells and ShermanWells & Sherman (2001) utilised the ‘transition time scale’ from Reference Finnigan and IveyFinnigan & Ivey (1999), initially derived for buoyancy-driven exchange between basins separated by a sill with a localised destabilising surface buoyancy flux. A fundamental result by Reference Wells and ShermanWells & Sherman (2001) was that a system can reach a quasi-steady state if the area of the shallow basin is equal to the area of the deep basin. As we discuss later, the theoretical and experimental study by Reference Finnigan and IveyFinnigan & Ivey (1999) on the buoyancy-driven interbasin exchange flow over a sill remains pivotal for understanding the transient dynamics of a TS. The latter authors identified three dynamic regimes in the exchange flow across the sill: localised convection, progressive growth of density difference and the development of a transient exchange flow until a quasi-steady state. The characterisation of the regime linked to horizontal convection, observed by Reference Finnigan and IveyFinnigan & Ivey (1999), relied on an inertia–buoyancy balance, similar to the analysis by Reference PhillipsPhillips (1966) describing convective circulation in the Red Sea.

Figure 5. Photograph showing the circulation of dye in a laboratory experiment. The red dye is placed into the convective mixing layer and the blue dye is placed into the gravity current. Adapted from Reference WellsWells (2001).

Utilising laboratory and numerical experiments, Reference Bednarz, Lei and PattersonBednarz, Lei & Patterson (2008a, Reference Bednarz, Lei and Patterson2009a) studied the onset of natural convection and the subsequent horizontal convective flow resulting from a destabilising surface heat flux on an initially isothermal water body of uniform slope connected to a flat interior basin. Although the authors focused on characterising the early convective stage theoretically, their experiments showed the existence of spatially varying convective regimes that motivated further numerical studies, including the work by Reference Mao, Lei and PattersonMao et al. (2010).

4.1 Conceptual model for TS in stratified water bodies

Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. (2022) proposed a two-dimensional (2-D) conceptual model to investigate the development of TS driven by surface cooling in water bodies of varying depth with temperatures $T\geq T_{md}$. The model integrates canonical topographic features observed in natural aquatic systems: an almost flat shallow plateau (or shelf) of length $\ell _{p}$, joined through a sloping (S) bottom region to the interior (I) stratified basin, which has a surface layer $h_{m}$, as shown in figure 3(c). The characteristic scales of the plateau are its horizontal length $\ell _{p}$ and its depth $d$. The depth-dependent sloping zone has a horizontal length $\ell _{s}$ and a characteristic slope $\bar {s}=(h_{m}-d)/\ell _{s} = h'_{m}/\ell _{s}$. The sloping bottom zone extends until the point where the base of the surface mixed layer intersects the bottom boundary. In the case studied here, the maximum depth of the water body $D$ is much smaller than the horizontal length scale $\mathcal {L}$ of the system, i.e. $D/\mathcal {L}\ll 1$, which allows treating these system environments as ‘shallow waters’.

The conceptual model in figure 3(d) considers an initial water temperature distribution defined by a stable and smooth two-layer stratification, i.e.

(4.1)\begin{equation} T(x,z,t = 0) = T_{b} + \frac{\Delta T}{2}\Bigl\{1 + \tanh\Bigl(\frac{z-z_{m}}{\delta_{m}}\Bigr)\Bigr\}, \end{equation}

with $T_{b}$ and $T_{s} = T_{b}+\Delta T$ representing the bottom layer temperature and surface layer temperature, respectively. The transition between the bottom layer and the surface layer occurs over a length scale $\delta _{m}$ and the position of the maximum temperature gradient, or thermocline, is initially at $z_{m}$. Both layers may exhibit weak temperature gradients or weak thermal stratification before a phase of destabilising surface cooling starts. The model (4.1) for the background temperature distribution represents well various thermally stratified lake systems featuring distinct epilimnion, metalimnion and hypolimnion regions. From the equation of state of water (2.1) and (4.1), we can determine the background density structure and the density difference between the top layer and bottom layer, $\Delta \rho$, resulting from the temperature difference between layers $\Delta T$. This density difference sets the actual reduced gravity that fluid parcels experience within the stratified body, $g' = (\Delta \rho /\rho _{0})g$.

Considering the conceptual model in figure 3(c), Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. (2022) inferred the existence of various convective regimes and phases of the development of TS. In this section we revisit the derivation of the theoretical time scales associated with the emergence of such convective regimes, starting from the time taken for thermal instabilities to grow from the air–water interface to the development of a cross-shore overturning circulation between the shallow and deep interior waters. The time scales involved in the lifespan of TS and discussed in this review are summarised in figure 6. Some of these time scales can only be determined from high-fidelity numerical experiments or highly controlled laboratory experiments, such as the time scale associated with the onset of thermal instabilities (Reference Bednarz, Lei and PattersonBednarz, Lei & Patterson 2008b, Reference Bednarz, Lei and Patterson2009b). In contrast, other time scales can be readily estimated from parameters directly measured in field experiments and comprehensive numerical studies. These time scales enable the interpretation of the different convective dynamics involved in the development of TS, as shown by Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. (2022) and Reference Ramón, Ulloa, Doda and BouffardRamón et al. (2022).

Figure 6. Summary of convective regimes associated with the development of TS driven by night time cooling.

4.2 Onset of thermal instabilities – $\tau _{B}$

We start by discussing the emergence of thermal instabilities that will originate primal motions. Let us assume that the fluid is initially at rest, with a uniform temperature near the surface. The air–water interface is subject to a ‘surface cooling’ rate $H_{0}$. In this scenario, water molecules at the surface interface cool faster than those beneath the surface. Prior to the onset of thermal instabilities, the time-dependent temperature drop at the surface is determined by the following expression (Reference Bednarz, Lei and PattersonBednarz et al. 2009a):

(4.2)\begin{equation} \Delta T(t) = {-}2\frac{H_{0}}{K}\sqrt{\frac{t}{\kappa {\rm \pi}}}. \end{equation}

Here $K$ is the water thermal conductivity. Such a temperature drop occurs over the conductive thermal boundary layer beneath the surface, which scales as $\delta _{T}\sim \sqrt {\kappa t/{\rm \pi} }$. The onset time scale of the instabilities can be estimated considering the critical Rayleigh number for free-slip boundaries, $Ra_{c}\approx 657.5$. Utilising expression (4.2), we can determine a time-dependent local Rayleigh number

(4.3)\begin{equation} Ra_{L}(t) \equiv \frac{g\alpha\Delta T(t) \delta^{3}_{T}(t)}{\nu\kappa} = 2g\alpha \frac{H_{0}}{K\nu\kappa^{2}}\left(\frac{\kappa t}{\rm \pi}\right)^2 \end{equation}

and evaluate the time needed for (4.3) to overcome the critical Rayleigh number. From this analysis, it yields a time scale for the onset of thermal instabilities (Reference Bednarz, Lei and PattersonBednarz et al. 2009a):

(4.4)\begin{equation} \tau_{B} \simeq \frac{Ra_{c}}{Ra_{L}} = \sqrt{657.5}\sqrt{\frac{\nu}{B_{0}}}. \end{equation}

For $B_{0}\sim 10^{-10}$–$10^{-8}$ ${\rm W}\,{\rm kg}^{-1}$, $\tau _{B}$ may vary from minutes to tens of minutes.

4.3 Characteristic time scale of free convection – $\tau _{RB}$

Upon the initiation of convection, the sinking speed of thermal plumes can be estimated using an energy scaling relationship between kinetic energy transport and buoyancy flux, known as the Deardorff convective velocity scale $w_{c}\simeq (B_{0}h_{m})^{1/3}$ (Reference DeardorffDeardorff 1970). The magnitude of $w_{c}$ is of the order of $10^{-3}\,{\rm m}\,{\rm s}^{-1}$, and $h_{m}$ represents the convective (vertical) length scale and is in the range of $O(1\unicode{x2013}10)\,{\rm m}$. These scales allow us to derive a characteristic time scale for free convection, analogous to Rayleigh–Bénard-type convection:

(4.5)\begin{equation} \tau_{RB} \simeq \frac{h_{m}}{(B_{0}h_{m})^{1/3}}. \end{equation}

This represents the time needed for the initial convective plumes to interact with either the sediment–water interface or the thermocline. In simpler terms, the time scale $\tau _{RB}$ characterises the lifespan of the very early free convective regime, where the system remains unaffected by the presence of the bottom boundary.

The early stage of free convection is brief, and this review does not delve deeply into this initial convective phase. Interested readers are referred to Reference Bednarz, Lei and PattersonBednarz et al. (2008a, Reference Bednarz, Lei and Patterson2009a) for more details. Here, our emphasis is on understanding the convective processes that unfold after thermals engage with the lower physical boundaries, leading to the development of a large-scale overturning circulation across the sloping region – the TS.

4.4 Transition towards a quasi-steady convective regime

Figure 6 illustrates three distinct convective stages followed by a baroclinic adjustment (see § 5). In the initial phase, fluid motion is dominated by local quasi-isotropic convective cells across the various zones, as shown in figure 6 (regime I). Rayleigh–Bénard-type convection facilitates the local vertical mixing of temperature in the surface layer. Simultaneously, a cross-shore temperature gradient develops due to differential cooling between shallow and deep waters (e.g. volume difference), leading to the gradual horizontal expansion of convective cells. This initial phase occurs within the time window $\tau _{RB}< t\lesssim \tau _{t}$.

Following the ‘transition time scale’, $\tau _{t}$, the fluid motion becomes progressively dominated by a horizontal density gradient. This sets the stage for initiating a large-scale overturning circulation from zone (S), which holds the most significant cross-shore density gradient. This second phase is termed regime II. Figure 6 illustrates the expansion of the horizontal convective cell toward the shallowest lateral boundary due to baroclinic adjustment. The lateral growth happens at a rate $u_{f}\equiv \textrm {d}\ell _{f}/\textrm {d} t$, with $\ell _{f}$ denoting the position of the shore-directed front from where the horizontal convection region originates. The large-scale overturning circulation takes the form of a two-layer exchange flow, with a surface layer $h_{s}$ moving toward the lateral boundary and the bottom layer flowing into the interior. Regime II unfolds over a time scale $\tau _{a}$, representing the time required by the large-scale overturning circulation to reach the (shallowest) lateral boundary, as depicted in figure 6. After a time scale $\tau _{qs} \simeq \tau _{t}+\tau _{a}$, we anticipate the large-scale overturning circulation to attain a quasi-steady state, referred to as regime III.

Next, we derive the dynamic regimes and associated time scales that govern the transition from one regime to the subsequent one. For conciseness, mathematical details are reported in Appendices A–D.

The evolution of the depth-averaged temperatures in the three zones (P), (S) and (I), as illustrated in figure 6, can be approximated through a simplification of the heat balance during the free convection regime (e.g. regime I). The details of the calculation are provided in Appendix B. We denote the vertical average of a function $\varphi (t,z)$ over a layer $h$ as $\langle \mathcal {\varphi }\rangle _{{(\cdot )}} = h^{-1}\int _{h}\varphi (t,z)\,{\rm {d}}z$. Here, the subscript $(\cdot )$ represents whether it pertains to zone (P), (S) or (I).

Starting from the initial time $t_{0}=0$, the average temperature within zone (P) changes as follows:

(4.6)\begin{equation} \langle T\rangle_{{(P)}} = \langle T(t_{0})\rangle_{{(P)}} -\frac{I_{0}}{d}t. \end{equation}

Here $I_{0}$ ($= H_{0}/(\rho _{0}C_p)$) is the kinematic heat flux. In zone (I) the average temperature can be approximated as

(4.7)\begin{equation} \langle T\rangle_{{(I)}} \approx \langle T(t_{0})\rangle_{{(I)}} -\frac{I_{0}}{h_{m}}t. \end{equation}

In the sloping region (S), the $x$-dependent depth-averaged temperature evolution, $\langle T\rangle _{{(S)}}$, follows a similar pattern to (4.6) and (4.7) but is normalised by the local depth, $D_{B}(x)$.

Given our assumption that $T>T_{md}$ and the small temperature drops during the cooling phase throughout the day relative to $T_{s}$ (approximately $0.1^{\circ }{\rm day}^{-1}$), we can rely on a local linearization of the equation of state (2.1) to accurately estimate the density, $\rho$, and buoyancy, $b$, for each zone:

(4.8a,b)\begin{equation} \frac{\langle\rho\rangle_{{({\cdot})}}}{\rho_{0}}-1 = \alpha(\langle T \rangle_{{({\cdot})}} - T_{0})\quad {\rm and} \quad \langle b \rangle_{{({\cdot})}} = {-}g\left(\frac{\langle\rho\rangle_{{({\cdot})}}}{\rho_{0}}-1\right) = {-} g \alpha(\langle T \rangle_{{({\cdot})}} - T_{0}). \end{equation}

The time-dependent, cross-shore pressure gradient between zones (P) and (I) can be estimated as (see Appendix D)

(4.9)\begin{equation} \left.\frac{1}{\rho_{0}}\frac{\partial p}{\partial x} \right|_{z = D-d} \approx \frac{d g \alpha}{\ell_{s}}(\langle T \rangle_{{(I)}} - \langle T \rangle_{{(P)}}) = {-}\frac{B_{0}}{\ell_{s}}\left(1 -\frac{d}{h_{m}}\right)t. \end{equation}

In general, the horizontal temperature difference, $\langle T\rangle _{{(I)}} - \langle T \rangle _{{(P)}}$, is substantially smaller than the temperature difference between the top and bottom layers in the stratified interior region (I), $\Delta T_z$.

4.5 Transition time scale – $\tau _{t}$

We aim at determining the time scale at which the cross-shore pressure gradient balances the inertial terms in (D2), i.e.

(4.10)\begin{equation} -\frac{\partial}{\partial x}\left(\frac{p}{\rho_{0}}\right) \simeq \frac{\partial u}{\partial t}\simeq \frac{\partial}{\partial x}\left(\frac{u^{2}}{2}\right). \end{equation}

From (4.10) and the cross-shore pressure gradient in (4.9), we derive two horizontal velocity scales:

(4.11)\begin{gather} \frac{\partial u}{\partial t}\simeq{-}\frac{\partial}{\partial x} \left(\frac{p}{\rho_{0}}\right) \quad\rightarrow\quad u\simeq \frac{B_{0}}{\ell_{s}}\left( 1 -\frac{d}{h_{m}}\right)\frac{t^{2}}{2}, \end{gather}

(4.12)\begin{gather}\frac{\partial}{\partial x}\left(\frac{u^{2}}{2}\right)\simeq{-}\frac{\partial}{\partial x} \left(\frac{p}{\rho_{0}}\right) \quad\rightarrow\quad u\simeq \sqrt{2 B_{0} \left(1 -\frac{d}{h_{m}}\right)t}. \end{gather}

The three-way balance (4.10) is achieved at a time scale $\tau _{t}$ at which the velocity scale $u$ allows balancing the inertial terms with the cross-shore pressure gradient. Therefore, equating expressions (4.11) and (4.12) yields

(4.13)\begin{equation} \tau_{t}\simeq \frac{2\ell_{s}}{(B_{0}\ell_{s})^{1/3}}\left(1-\frac{d}{h_{m}}\right)^{{-}1/3}. \end{equation}

The time scale $\tau _{t}$ defines the transition from regime I to regime II, schematised in figure 6, and at $\tau _{t}$, the characteristic horizontal velocity scale, $u_{t}$, is given by

(4.14)\begin{equation} u_{t}\simeq 2(B_{0}\ell_{s})^{1/3} \left(1-\frac{d}{h_{m}}\right)^{1/3}. \end{equation}

The exchange flow associated with the large-scale overturning circulation is expected to develop in the zone with the maximum cross-shore pressure gradient, i.e. zone (S). Considering that this zone has a uniform slope, we expect the large-scale overturning circulation to emerge at about the middle of zone (S), as shown later via numerical experiments.

4.6 Adjustment and quasi-steady time scales – $\tau _{a}$, $\tau _{qs}$

We now aim to determine the time scale required by the large-scale overturning circulation to reach the lateral boundary in zone (P) and be adjusted to a new equilibrium state. This adjustment time scale, $\tau _{a}$, depends on the horizontal length of the nearshore area between the point from where the TS starts and the lateral boundary. During regime II, this nearshore area can be split into two regions: a region of length $\ell _{f}(t)$ (measured from the starting point) that grows in time due to the lateral expansion of the large-scale overturning circulation, and its neighbouring region that shrinks in time, where quasi-isotropic convective cells dominate the local flow, as depicted in figure 6. In this shrinking well-mixed zone, denoted as (LP), the average buoyancy results from the expressions (B3) and (4.8a,b),

(4.15)\begin{equation} \langle b \rangle_{{(LP)}} = {-}\frac{B_{0}}{d} t. \end{equation}

In zone (LP) the temperature and buoyancy are continuously decreasing, and $\langle b \rangle _{{(LP)}}$ matches the buoyancy at the left front of the horizontal convective cell, $\ell _{f}(t)$. Following the arguments by Reference PhillipsPhillips (1966), we assume that the momentum of the surface layer current $u_{s}$ and buoyancy are found in an inertia–buoyancy balance. This assumption implies that the flow within the horizontal convective cell scales as

(4.16)\begin{equation} u_{s}\simeq (2\langle b\rangle_{{(LP)}} h_{s})^{1/2}, \end{equation}

where $h_{s}$ is the thickness of the surface exchange layer, shown in figure 6. Additionally, we may also consider that the surface buoyancy flux, $B_{0}$, across $\ell _{f}$ is balanced with the production of kinetic energy, which leads to the scaling velocity $u_{s}\simeq (2B_{0}\ell _{f})^{1/3}$. Hence, the above two velocity scales are equal when

(4.17)\begin{equation} \langle b \rangle_{{(LP)}} \simeq \frac{(2B_{0}\ell_{f})^{2/3}}{2h_{s}}. \end{equation}

By equating (4.15) and (4.17), we can derive the time scale required by the exchange flow to reach the lateral boundary. The latter is reached when $\ell _{f}(\tau _{a})\simeq \ell _{p}+\ell _{s}/2$,

(4.18)\begin{equation} \tau_{a} \simeq \frac{d}{h_{s}} \frac{\ell_{f}}{(2\ell_{f}B_{0})^{1/3}}, \end{equation}

where $\tau _{a}$ represents the time scale needed for the large-scale overturning circulation to fully adjust across zones (P) and (S). As a consequence, the time scale $\tau _{qs}\simeq \tau _{t}+\tau _{a}$ represents the time required to achieve the quasi-steady state that characterises regime III.

4.7 What if the surface layer is thermally stratified?

The surface layer, schematised in figure 3(c), may have an initial temperature gradient, particularly in the vicinity of the very surface, characterised by a buoyancy frequency $N_{s}$ (e.g. Reference ImbergerImberger 1985). In this situation, the destabilising surface cooling will erode and mix the surface layer until it encounters a physical barrier. This barrier could be either the sediment–water interface in the littoral region or the seasonal thermocline zone in the interior or pelagic zone. We can estimate the time required to mix the surface layer by using the evolution equation for a convective mixing layer (Reference ZilitinkevicZilitinkevic 1991) given as follows:

(4.19)\begin{equation} \frac{{\rm d}h}{{\rm d}t} = (1 + 2A)\frac{B_{0}}{hN^{2}_{s}}. \end{equation}

Here, the coefficient $A\approx 0.2$ is an empirical parameter. For a detailed derivation of (4.19) and more information about $A$, readers are referred to Reference ZilitinkevicZilitinkevic (1991). It is important to note that the evolution equation (4.19) accounts explicitly for the expansion of the convective mixing layer of depth $h$ due to a destabilising surface buoyancy flux $B_{0}$.

Given that the lakes’ surface layers exhibit weak stratification due to daytime solar radiation absorption in the epilimnion, the density's vertical distribution tends to be approximately linear. Consequently, the buoyancy frequency $N_{s}$ is almost constant and significantly weaker than the buoyancy frequency associated with the seasonal thermocline, as illustrated in figure 3(c). Assuming $N_{s}$ to be ‘uniform’ through the surface layer, the evolution equation for the surface convective mixing layer (4.19) yields an analytical solution:

(4.20)\begin{equation} h(t) = \sqrt{ h^{2}_{0} + 2(1 + 2A)\frac{B_{0}}{N^{2}_{s}}(t-t_{0})}. \end{equation}

Here, $h_{0}$ represents the initial thickness of the convective mixing layer at the initial time $t_{0}$. If both initial conditions are taken as ‘zero’, a straightforward expression for the time taken to mix the surface layer is obtained. This time, denoted as the initial mixing time scale by Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. (2022), is given by

(4.21)\begin{equation} \tau_{mix} = \frac{h^{2}N^{2}_s}{2(1 + 2A)B_{0}}. \end{equation}

We should note that the time required to mix shallow waters may differ from the time needed to mix the surface layer within the interior stratified zone. If $\tau _{mix}>0$, time scales linked to the various convective regimes experience a temporal shift. In particular, the time scale for the ‘onset’ or development of a TS due to differential cooling is more accurately predicted by ‘$\tau _{mix}+\tau _{t}$’, as shown by Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. (2022).

4.8 Numerical and field experiments

The theoretical time scales predicting the different stages at which a TS develops can be tested with numerical and field experiments (Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. 2022; Reference Ramón, Ulloa, Doda and BouffardRamón et al. 2022; Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. 2022). Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. (2022) confirmed the theoretical transition time scale $\tau _{t}$ and quasi-steady time scale $\tau _{qs}$ through 2-D high-fidelity spectral numerical experiments at high Rayleigh numbers, covering a range of nearshore slopes observed in aquatic systems, $1\,\%<\bar {s}<10\,\%$. In their assessment, the authors introduced the flow geometry parameter, $F_{G}$, denoted as

(4.22)\begin{equation} F_{G} \equiv \frac{\sqrt{\langle u^{2}\rangle} + c}{\sqrt{\langle w^{2}\rangle} + c}, \end{equation}

where $\langle u\rangle$ and $\langle w\rangle$ represent the space-averaged horizontal and vertical velocity components over zone (P), defined as the plateau nearshore region, figure 6, respectively. The constant $c\ll \sqrt {\langle u\rangle ^{2}}, \sqrt {\langle w\rangle ^{2}}$ is a regularisation term preventing a singularity when fluid is at rest before the onset of thermal instabilities around $\tau _{B}$. Figure 7 shows the time series of $F_{G}$ for different nearshore topographic conditions, highlighting that $F_{G}$ captures the three convective regimes associated with TS development. During regime I, $t/\tau _{t}<1$, $F_{G}$ hovers around 1, representing a quasi-isotropic flow geometry where the magnitudes of horizontal and vertical velocity components are statistically similar – as shown in figure 7(a). In regime II, $t/\tau _{t}\geq 1$, $F_{G}$ accurately represents the progressive growth of an overturning circulation that starts in zone (S) but that propagates in time through zone (P) (figure 6). Here, the velocity field evolves into a more elliptical geometry, becoming predominantly horizontal, causing $F_{G}$ to boldly exceed 1. To identify regime III, figure 7(b) shows that $F_{G}$ reaches a plateau for $t/\tau _{qs}\geq 1$, indicating the theoretically expected quasi-steady phase of regime III. The final value is in the range $1< F_{G}<10$ for a wide range of forcing and bathymetry, yet, there is no study that systematically investigates a possible parameterisation for the plausible values of $F_{G}$ can attain at quasi-steady state. Before this phase, i.e. $t/\tau _{qs}<1$, $F_{G}$ exhibits the transient phase associated with regimes I and II. Similar parameterisation of the quasi-steady time scale can be obtained from three-dimensional (3-D) numerical experiments using a hydrostatic Reynolds-averaged Navier–Stokes (RANS) model (Reference Ramón, Ulloa, Doda and BouffardRamón et al. 2022).

Figure 7. Different time scales associated with the TS development, extracted from numerical experiments. Panel (a) shows the flow geometry, $F_{G}$, as a function of the non-dimensional time $t/\tau _{t}$. Circle A denotes the onset of thermal instabilities; circle B highlights the time for the transition to horizontal circulation (i.e. TS development). Panel (b) illustrates $F_{G}$ as a function of the non-dimensional time $t/\tau _{qs}$. All simulations hold $Pr =7$, and the destabilising thermal forcing, quantified by the Rayleigh number in (E1a–c), is $Ra=10^{5}$. Shaded areas correspond to the different convective regimes described in figure 6. Here $A^{(h)}$ is the ratio between the horizontal length scale of the plateau and the slope regions. Adapted from Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. (2022).

Field experiments carried out at Lake Rotsee, Switzerland, by Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. (2022) explored the validity of the conceptual model presented in § 4 and figure 6 in a real-scale system. These experiments accurately determined the development of TS in the littoral region of natural aquatic environments by the theoretical transition time scale, $\tau _{t}$. In-situ observations shown in figure 8 showcase two days characterised by almost uniform night-time surface cooling, resulting in a consistent destabilising surface buoyancy flux $B_{0}$ that triggered TS after an initiation period of duration $\tau _{mix}+\tau _{t}$. Using the time scales presented in § 4, it is also possible to infer the seasonality of TS in a specific lake. Thermal siphons would form on a given day if the cooling period when $B_0>0$ lasts longer than the initiation period $\tau _{mix}+\tau _{t}$. For the case of Lake Rotsee (Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. 2022), the mixing time scale $\tau _{mix}>O(10)$ h was too long in early summer (e.g. $N_s^2$ is too strong) and TS could barely develop despite the high $B_0$. The optimal period for TS formation was observed in fall with still elevated $B_0$ but reduced $N_s^2$, which decreased $\tau _{mix}< O(1)$ h. In winter, $B_0$ was strongly reduced and $\tau _{t}$ increased from $O(1)$ h in summer to $O(10)$ h, which limited the development of TS. Overall, the theoretical model outlined in § 4 precisely captures the fluid dynamics observed in lakes (Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. 2022).

Figure 8. Panel (a) illustrates Lake Rotsee bathymetry. Panel (b) shows the cross-shore topographic profile of the northeast side basin, highlighting the average depth of zone (P), $d\approx 1$ m, the average slope $\overline {s}\approx 3\,\%$ of zone (S) and the depth of the surface mixed layer $h_{m}\approx 4.8$ m. (c) Time series of surface buoyancy fluxes defining the cooling and heating phases on 9–10 September 2019. (d) Vertical velocity profile. Positive (purple) values correspond to a flow moving upward. Strong vertical movements are the signature of convective plumes. (e) Cross-shore velocity profile. Positive (purple) values correspond to a flow moving offshore (southwestward flow, highlighted by the blue line). In (d,e) ‘Height’ refers to the distance from the bottom at the observation location. Black lines are 0.05 $^{\circ }$C isotherms linearly interpolated from a thermistor chain (vertical resolution indicated with the horizontal ticks on the right axis). The theoretical time scales $\tau _{mix}$, $\tau _{t}$ and $\tau _{qs}$ are indicated. They define the duration between the start of the cooling phase and (i) the end of the vertical mixing period ($\tau _{mix}$), (ii) the end of the initiation period ($\tau _{t}+\tau _{mix}$) and (iii) the quasi-steady-state regime ($\tau _{t}+\tau _{qs}$). Two dynamic phases, defined as the convective (C) and relaxation (R) phases, are identified after the TS formation (see § 6.3). Adapted from Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. (2022) and Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. (2022).

6.1 Studies considering periodic thermal forcing

Motivated by the diurnal cycle of TS outlined by Reference Monismith, Imberger and MorisonMonismith et al. (1990), Reference Farrow and PattersonFarrow & Patterson (1993) developed a linear model incorporating a simple diurnal heat flux cycle. This model aimed to capture both the cooling and heating phases within a fluid body characterised by uniform and small slopes. In this idealised representation, the heat source and sink were vertically distributed across the water column with varying depths, establishing a horizontal density gradient that drove phases of differential heating and cooling over a diurnal cycle. The authors derived analytical expressions for the velocity field and asymptotic solutions, delineating the diffusive regime observed in shallow waters and the inertial regime observed in deeper waters (Reference Mao, Lei and PattersonMao et al. 2009). Their findings were validated through fully nonlinear numerical simulations conducted under weak thermal forcing conditions.

The elegant asymptotic approach pioneered by Reference Farrow and PattersonFarrow & Patterson (1993) in the context of periodically forced sloping water bodies was generalised for systems with gradually varying topography by Reference FarrowFarrow (2004). These solutions were subsequently extended to encompass more complex environments and diverse periodic forcing mechanisms by various authors (Reference Bednarz, Lei and PattersonBednarz et al. 2009b; Reference FarrowFarrow 2013b; Reference Lin and WuLin & Wu 2014; Reference LinLin 2015; Reference FarrowFarrow 2016; Reference Ulloa, Davis, Monismith and PawlakUlloa et al. 2018; Reference MaoMao 2019; Reference Mao, Lei and PattersonMao, Lei & Patterson 2019; Reference Coenen, Sánchez, Félez, Davis and PawlakCoenen et al. 2021). While these models may not capture the complexness of the highly nonlinear dynamics observed in nonlinear simulations, laboratory experiments and field observations, they effectively encapsulate the fundamental features of TS, especially their characteristic diurnal rhythm.

Developing novel experiments, Reference Sturman and IveySturman & Ivey (1998) explored the convective exchange phenomena between shallow and deep basins within a water body due to periodic surface thermal forcing. The authors induced differential heating and cooling in a system where a uniform slope region interconnected a shallow reservoir and a deep reservoir. The results revealed that a destabilising heat flux in the shallow region initiated a vigorous exchange flow between the shallow and deep basins. As in the case of steady surface cooling, the resulting discharge per unit width was proportional to $q_{d}\sim (B_{0}\ell _{p})^{1/3}d$, where $B_{0}$ denotes the destabilising surface buoyancy flux across the horizontal extent $\ell _{p}$ of the shallow plateau region with a depth $d$. In contrast, experiments involving a stabilising surface buoyancy flux within the same shallow plateau region exhibited a discharge per unit width proportional to $q_{\delta }\sim (B_{0}\ell _{p})^{1/3}\delta$, with $\delta < d$ characterising the thermal boundary layer thickness of the exchange flow.

6.2 Diurnal cycle in the field

In temperate water bodies, atmospheric heat exchange causes a diurnal cycle of daytime heating and night-time cooling, with respective intensities modulated by seasons. In a simplified conceptual model we can consider that the water body uniformly loses heat at a rate $H_{net,0}(t) = H_{Q,0}(t) > 0$ over a given time scale $\tau _{cool}$ and gain heat at a varying rate $H_{net,0}(t) < 0$ over the rest of the day $\tau _{warm}$, as can be seen, for instance, with the hourly resolved heat flux measurements in Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. (2022). These two phases define the cooling–heating diurnal cycle. For a sloping topography, each phase is associated with the formation of lateral temperature gradients: the cooling–heating diurnal cycle generates a differential cooling–heating cycle with some inertia (Reference Monismith, Imberger and MorisonMonismith et al. 1990; Reference James, Barko and EakinJames et al. 1994).

Recent measurements in a wind-sheltered elongated lake (Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. 2022) highlighted how forcing conditions affect the formation and destruction of TS over a diurnal cycle during the summertime period (figure 8). The net surface buoyancy flux $B_{net,0}$ ranged from $-2\times 10^{-7}$ W kg$^{-1}$ during the daytime to $1\times 10^{-7}$ W kg$^{-1}$ during night time (figure 8c), with the change in sign mostly due to the solar radiation during daytime ($B_{SW,0}<0$). While not shown, the wind remained negligible during and before the period of measurements. The vertical velocity recorded with acoustic Doppler current profilers (ADCPs, all details in Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. 2022) confirmed that the net night-time cooling was strong enough to trigger convective mixing, with organised upward and downward plumes reaching vertical velocities of 5 mm s$^{-1}$ (figure 8d). Such convective mixing could erode and deepen the daily shallow stratification, as can be observed with the 0.05 $^\circ$C isotherm lines gradually becoming vertical (e.g. homogeneous profile) during the first part of the night (figure 8e). Interestingly, a cross-shore circulation started 3 h after complete vertical mixing, with a marked downslope flow near the bottom and a weaker opposite flow above (note that the upward-looking ADCP did not resolve the velocity near the surface). An interesting dynamic observed by Reference Doda, Ulloa, Ramón, Wüest and BouffardDoda et al. (2023) was the striking vertical oscillations of the interface between the downslope flow and the onshore, with fluid excursion spanning approximately 1–2 m, representing nearly 100 % of the thickness of the TS.

Results reported by Reference Doda, Ulloa, Ramón, Wüest and BouffardDoda et al. (2023) showed that, by the end of the cooling phase, the water column along the sloping zone exhibited thermally controlled stratification near the bottom. This increase in water density at the base of the water column was attributed to the persistent flow of colder water downslope. Another striking observation is the intensification of the TS in phase with the weakening of vertical convection at the onset of radiative heating. The cross-shore velocity increased above 2 cm s$^{-1}$ and the induced bottom stratification reached $\partial T/\partial z\sim 0.1\,^\circ$C m$^{-1}$ at the beginning of the heating phase. This flushing lasted several hours while radiative heating increased and re-stratified the water column. This example, supported by consistent observations, challenges the traditional assumption that most of the transport due to night cooling occurs at night; instead, it highlights that the primary flushing occurs in the early morning. This underscores the need for a quantitative understanding of physical processes to elucidate the timing, formation and fate of TS. In the next section we delve into a more detailed review of the effects of convective mixing on TS dynamics.

6.3 Two dynamic phases of transport modulated by penetrative convection

Time-varying forcing highlights two contrasted dynamical structures of the TS convective flow (Reference Doda, Ulloa, Ramón, Wüest and BouffardDoda et al. 2023). Qualitatively, the first phase is characterised by a weak density current of fluctuating thickness and active free convection taking place in the upper layer. This phase is defined as the convective phase, C phase. The second phase is characterised by an intensification of the density current, a more stable current interface over time and a decay of convection in the upper layer, denoted as the relaxation phase, R phase. These two distinctive regimes are schematised in figure 10(a).

Figure 10. Subsets of the convective and relaxation phases captured in the sloping region on 6–7 November 2019 in Rotsee with schematics of the processes in (a). (b) High frequency cross-shore velocity above the sediment. Solid black and blue lines depict the high frequencyand smoothed upper interface of the gravity current, respectively. The dashed horizontal line indicates the average thickness $\overline {h_{d}}$. Grey lines are $0.02\,^{\circ }$C isotherms and numbers refer to temperatures. (c) Same as (b) for the vertical velocity. Arrows highlight strong convective downdrafts and updrafts penetrating across the interface. (d) Burst-averaged time series of penetrative length scale $\delta _c$, shear interface thickness $\delta _S$ and density interface thickness $\delta _{\rho }$ with standard deviation indicated with the shaded colours. Adapted from Reference Doda, Ulloa, Ramón, Wüest and BouffardDoda et al. (2023).

The lower limb observed during the C phase is characteristic of a density current propagating in an active turbulent environment. Such a density current does not follow established scaling. The propagation of density currents under turbulent background conditions has been studied in the laboratory (Reference Linden and SimpsonLinden & Simpson 1986; Reference SimpsonSimpson 1986). The findings of these laboratory experiments are supported by recent field experiments (Reference Doda, Ulloa, Ramón, Wüest and BouffardDoda et al. 2023), showing that turbulence in the convective upper layer eroded the stratified downslope flow. This process is similar to the case of a uniform destabilising buoyancy flux competing with a local buoyancy flux that produces a stable density stratification (Reference Wells, Griffiths and TurnerWells, Griffiths & Turner 1999). The convective Richardson number, $Ri_c$, presented in (E1a–c) is a convenient non-dimensional number to evaluate how convective plumes impinging at density interfaces alter the stratified flow (Reference BainesBaines 1975; Reference Noh, Fernando and ChingNoh, Fernando & Ching 1992; Reference Cotel and KudoCotel & Kudo 2008). Here $Ri_c$ is a measure of the relative intensity of penetrative convection with respect to the stratified downslope flow. Below a critical value of ${Ri}^{{(crit)}}_{c}\approx 10$, convective plumes can penetrate across the stratified layer (Reference Noh, Fernando and ChingNoh et al. 1992; Reference Cotel and KudoCotel & Kudo 2008). The penetration depth of convective plumes can be defined as

(6.1)\begin{equation} \delta_c = \frac{|w_c|}{N_d}, \end{equation}

where $N_d$ is the depth-averaged buoyancy frequency within the gravity current. The dominant role of convection over shear-driven disturbance across the current interface can be assessed by comparing $\delta _c$ with the thickness of the shear ($\delta _S$) and density ($\delta _\rho$) interfaces (Reference Zhu and LawrenceZhu & Lawrence 2001): $\delta _S={(U_d-U_a)}/{(\partial u/\partial z)_{max}}$ and $\delta _{\rho }=(\rho _d-\rho _a)/(\partial \rho /\partial z)_{max}$. Indices $d$ and $a$ refer to the density current and the ambient fluid, respectively.

The driving processes responsible for the distinct dynamics of the two phases can be retrieved from a statistical description of $Ri_{c}$, $\delta _c$, $\delta _S$ and $\delta _{\rho }$. The C phase is characterised by $\delta _c>\delta _S,\delta _{\rho }$, suggesting that convective plumes penetrate beyond the current interface defined from shear and density and erode the stratified layer (figure 10). Furthermore, ${Ri}_{c}$ increased from values close to $Ri^{{(crit)}}_{c}$ during the C phase to orders of magnitude larger, confirming that the TS dynamics depend on the interaction between penetrative convection and the gravity currents. The enhancement of vertical fluxes between the downslope flow and the ambient layer by convective mixing diluted transported tracers and diminished the cross-shore exchange during the C phase. A similar reduction in the lateral transport of tracers due to turbulent diffusion has been reported for exchange flows (Reference HelfrichHelfrich 1995; Reference Winters and SeimWinters & Seim 2000; Reference Hogg, Ivey and WintersHogg, Ivey & Winters 2001). Reference Linden and SimpsonLinden & Simpson (1986) observed that the flow became vertically mixed after a distance $L_x$. In the case of penetrative convection, this distance represents how far the gravity current propagates before being entirely eroded by penetrative plumes and can be defined as $L_x\sim U_s h_d/w_e$, where $w_e=A w_c/{Ri}_c=A B_0 h_m/(g' h_d)$ is the convective entrainment velocity, with $A$ an empirical coefficient (Reference Deardorff, Willis and StocktonDeardorff, Willis & Stockton 1980). Considering the presented case of Lake Rotsee, a numerical example is given with $g'\sim g (\partial \rho /\partial x)\Delta x/\rho _0\approx 4\times 10^{-4}\,{\rm m}\,{\rm s}^{-2}$, $h_d\approx h_m\approx 2\,{\rm m}$, $B_0\approx 1\times 10^{-7}\,{\rm W}\,{\rm kg}^{-1}$ and $A\approx 0.2$ (Reference Sullivan, Moeng, Stevens, Lenschow and MayorSullivan et al. 1998), and we retrieve $L_x\approx 0.5{g'}^{3/2}h_{l}^{1/2}h_d^2/(AB_0 h_m)\approx 700\,{\rm m}$, which is consistent with observations of complete mixing (Reference Doda, Ulloa, Ramón, Wüest and BouffardDoda et al. 2023). Such a degeneration of the downslope stratified layer before reaching the lake interior confirms the numerical results by Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. (2022) during the C phase. In these numerical experiments and in Rotsee, the deep stratified region represents a larger portion of the surface area compared with the shallow littoral region such that the ratio of deep to shallow length scales defined in Reference Wells and ShermanWells & Sherman (2001) is $\mathcal {R}=\ell _{deep}/\ell _{shallow}=(\ell _{basin}-\ell _p-0.5\ell _s)/(\ell _p+0.5\ell _s)>1$, where $\ell _{basin}$ is the distance from the shore to the deepest point. From laboratory experiments reaching steady state, Reference Wells and ShermanWells & Sherman (2001) found that the stratified layer induced by the gravity current is reduced for $\mathcal {R}>1$, which is consistent with the observed erosion of the downslope flow in Reference Ulloa, Ramón, Doda, Wüest and BouffardUlloa et al. (2022) and Reference Doda, Ulloa, Ramón, Wüest and BouffardDoda et al. (2023). Yet, the role of the basin topography in the competition between TS and penetrative convection in the field requires further research.

The transition from the C phase to the R phase occurred at sunrise when the net cooling heat flux diminished and the intensity of convection weakened. This caused a sudden drop of $|(\partial \rho /\partial x)|$ and a baroclinic adjustment with an intensification of the downslope gravity current. The maximal cross-shore transport was observed during the R phase in Lake Rotsee (Reference Doda, Ramón, Ulloa, Wüest and BouffardDoda et al. 2022), in agreement with numerical simulations of TS (Reference Chubarenko, Esiukova, Stepanova, Chubarenko and BaudlerChubarenko et al. 2013; Reference Safaie, Pawlak and DavisSafaie, Pawlak & Davis 2022). This relaxation is comparable to the frontogenesis observed by Reference Linden and SimpsonLinden & Simpson (1986) once turbulence was turned off. The intensification of gravity currents with weaker turbulent mixing has also been observed in the atmosphere (Reference Simpson, Mansfield and MilfordSimpson, Mansfield & Milford 1977; Reference Sha, Kawamura and UedaSha, Kawamura & Ueda 1991; Reference Parker, Burton, Diongue-Niang, Ellis, Felton, Taylor, Thorncroft, Bessemoulin and TompkinsParker et al. 2005) and estuaries (Reference Hetzel, Pattiaratchi, Lowe and HofmeisterHetzel et al. 2015). In this R phase, the averaged downslope velocity $U_{d}\approx 1.5\pm 0.4\,{\rm cm}\,{\rm s}^{-1}\approx 0.9 U_{s}$ matched the velocity scale based on horizontal density gradients in a quiescent environment $U_{s}\sim 0.5\sqrt {gh_{l}\Delta \rho /\rho _0}$, where $\Delta \rho$ is the lateral density difference between shallow and deep regions and $\rho _0=1000$ kg m$^{-3}$ is the reference density. In the presence of convection during the C phase, however, it was only $U_{d}\approx 0.8\pm 0.3\,{\rm cm}\,{\rm s}^{-1}\approx 0.5U_{s}$. This observation is consistent with Reference Linden and SimpsonLinden & Simpson (1986) who demonstrated that the front velocity of lock-exchange gravity currents in turbulent conditions was less than the velocity scale $U_s$.

In addition to modifying the flow dynamics, penetrative convection also affects the growth of the gravity current by changing the amount of ambient water entrained into the downslope flow. The net shear entrainment coefficient is expressed as $E_{net} = ({\partial }/{\partial x})({U_dh_d}/{(U_d-U_a)}),$ where $U_d$ and $h_d$ are the depth-averaged velocity and thickness of the downslope flow, respectively, and $U_a$ is the depth-averaged velocity in the ambient layer. Here $E_{net}$ has been parameterised as a function of the bulk Richardson number $Ri_B=g'h_d\cos \theta /(U_d-U_a)^2$ in quiescent environments, where $\theta$ is the bottom slope angle (Reference Ellison and TurnerEllison & Turner 1959; Reference Cenedese and AdduceCenedese & Adduce 2010). The effect of convective plumes on $E_{net}$ is not well understood. Reference Fer, Lemmin and ThorpeFer et al. (2002b) reported higher $E_{net}$ in TS interacting with convective plumes compared with the parametrisation of Reference Ellison and TurnerEllison & Turner (1959), but their estimated $E_{net}$ was lower than expected from the large Reynolds number $Re=U_{d}h_{d}/\nu \approx 10^5$ (Reference Cenedese and AdduceCenedese & Adduce 2010). The contribution of convective plumes to $E_{net}$ could be examined in the field from simultaneous velocity measurements along a cross-shore transect. The effects of the return flow $U_a$ on the entrainment must be included since $h_d/h_a\approx 1$, as in counterflows (Reference Moore and LongMoore & Long 1971; Reference ChristodoulouChristodoulou 1986).

7.1 Intrusion in the thermocline

Thermal siphons contribute to transporting dissolved gases and compounds, modifying the lake ecosystem dynamics. These convective flows are often viewed as facilitators of deep oxygen ventilation and mixing with density currents flowing underneath the thermocline. Yet, for a TS to penetrate the thermocline, the initial horizontal temperature difference between littoral and surface pelagic waters ($\Delta T_x$) must be greater than the vertical temperature difference in the vicinity of the thermocline ($\Delta T_z$). Knowing the characteristics of the basin and the thermal stratification, and assuming that (i) the build-up of the horizontal temperature gradient (see § 4) is faster than the flushing time scale $\tau _{F}$ (5.4) and (ii) neglecting any entrainment, it is possible to estimate the buoyancy flux needed over the flushing time scale period to generate such underflow. From (4.9), the horizontal temperature difference that drives TS is

(7.1)\begin{equation} \Delta T_x(t) = \frac{B_0}{g\alpha d}\left(1-\frac{d}{h_{m}}\right) t. \end{equation}

The condition for the TS to generate an underflow is that $\Delta T_z < \Delta T_x$. At $t = \tau _F$, it yields

(7.2)\begin{equation} \Delta T_z < \Delta T_x = \frac{B_0}{g\alpha d} \left(1-\frac{d}{h_{m}}\right)\frac{\ell_{l}}{c_q(B_0 \ell_{ml})^{1/3}}. \end{equation}

For a deep thermocline, $d/h_{m} \ll 1$ and $\ell _{l} \approx \ell _{ml}$, the previous inequality reduces to $\Delta T_z < {(B_0\ell _{ml})^{2/3}}/{c_q g \alpha d}$.

Numerical applications with $B_0 = O(10^{-7})\,{\rm m}^2\,{\rm s}^{-3}$, $d = 10$ m, $\ell _{ml} = 1000$ m and $\alpha \sim 10^{-4}\,^{\circ }{\rm C}^{-1}$ suggest that $\Delta T_{z}$ should be smaller than $\Delta T_{x} \approx O(10^{-2})\,^{\circ }$C, meaning that the interior basin should be nearly well mixed. Therefore, TS are unlikely to penetrate beneath the thermocline or to lead to density current splitting at the thermocline as observed in Reference Cortés, Fleenor, Wells, De Vicente and RuedaCortés et al. (2014). From our perspective, other processes must coexist with TS to facilitate such deep ventilation from lateral boundaries (see, for instance, Reference Reiss, Lemmin and BarryReiss, Lemmin & Barry 2022; Reference Peng, Lemmin, Mettra, Reiss and BarryPeng et al. 2024). Future research should explore the interplay of complex bathymetric features, including sills, with wind and destabilising heat fluxes. This topic falls into the category of interbasin exchanges and is not discussed here. Another intriguing finding was reported by Reference Wells and ShermanWells & Sherman (2001) through a combination of laboratory experiments and theoretical considerations. The authors suggested that, under specific bathymetric conditions (see figure 5), TS can sustain a winter stratification in the deep basin from initially fully mixed conditions. This phenomenon may occur when (i) the ratio between the deeper and shallower areas of the lake is less than 1 and (ii) the ratio between shallow depth and deep depth is less than 0.5, resulting in the development of a surface mixing layer and a stratified deep region. These results challenge the long-standing view that deep lakes (e.g. here taken as water depths larger than the direct vertical extent of mixing by wind energy) can be fully mixed. Instead, as long as the lake has significant shallow areas (i.e. $h_{mean} \ll h_{max}$), TS prevent full vertical mixing induced by destabilising surface cooling by maintaining a weak bottom stratification. Detailed conductivity, temperature and oxygen profiles could reveal fine-scale structures in the deep water associated with such processes and facilitate progress in this research area.

While reaching the conditions for underflow from TS seems unlikely to be achieved in the field, the intrusion of the thermally driven density current into the thermocline can be identified from temperature and current profiles (figure 11). Reference Doda, Ramón, Ulloa, Brennwald, Kipfer, Perga, Wüest, Schubert and BouffardDoda et al. (2024) followed the path of natural dissolved oxygen and krypton injected at the shore with portable mass spectrometers. The authors showed that the peaks in krypton and dissolved oxygen matched the timing and depth of the observed peak in velocity in the thermocline region. Similar approaches for river intrusion can be used to infer the pathway of the density currents (Reference Cortés, Fleenor, Wells, De Vicente and RuedaCortés et al. 2014), the interaction with the water body (Reference Ramón, Priet-Mahéo, Rueda and AndradóttirRamón et al. 2020) and the length scale over which the intrusion diffuses into the background fluid (Reference Hogg, Marti, Huppert and ImbergerHogg et al. 2013).

Figure 11. Temporally averaged profiles of temperature (blue) and cross-shore velocity (pink) in the TS intrusion region in Rotsee ($\sim$600 m from the shore). The background averaged temperature profile collected at the lake centre ($\sim$1 km from the shore) is shown in black. Shaded areas depict the standard deviation. The intrusion creates a cold layer at the base of the mixed layer that is not observed at the lake centre. Adapted from Reference Doda, Ramón, Ulloa, Brennwald, Kipfer, Perga, Wüest, Schubert and BouffardDoda et al. (2024).

7.2 Large-scale circulation induced by TS: Observations from ice-covered lakes

Thermal siphons challenge the one-dimensional view of lakes’ dynamics with the interaction between bathymetry and surface heat fluxes driving dense water into the deeper basin. However, isolating this source of heat and momentum proves challenging amidst the multitude of processes acting on lakes. An efficient way to investigate the large-scale effect of TS on lakes is to look at ice-covered lakes where wind-driven kinetic energy is largely halted, providing an opportunity to isolate buoyancy-driven processes. Studying ice-covered freshwater lakes corresponds to transitioning to the opposite side of the $T_{md}$, where the equation of state now links an increase in temperature to a rise in density. In late winter, also known as the winter II phase (Reference KirillinKirillin et al. 2012), a diurnal cycle of penetrative radiative heating surpasses diffusive cooling at the ice–water interface, resulting in a net loss of buoyancy. This triggers both penetrative turbulent convection (Reference FarmerFarmer 1975; Reference KirillinKirillin et al. 2012; Reference Yang, Young, Brown and WellsYang et al. 2017; Reference BouffardBouffard et al. 2019) and downslope gravity currents, the latter stemming from differential heating in the shallows (Reference KenneyKenney 1996; Reference Stefanovic and StefanStefanovic & Stefan 2002; Reference Salonen, Pulkkanen, Salmi and GriffithsSalonen et al. 2014; Reference Kirillin, Forrest, Graves, Fischer, Engelhardt and LavalKirillin et al. 2015; Reference Ulloa, Winters, Wüest and BouffardUlloa et al. 2019; Reference Cortés and MacIntyreCortés & MacIntyre 2020).

Beneath the ice, the heating and deepening of the convective mixed layer are conventionally viewed as a one-dimensional vertical process (Reference FarmerFarmer 1975; Reference Mironov, Terzhevik, Kirillin, Jonas, Malm and FarmerMironov et al. 2002). Often, the role of lateral heat advection caused by gravity currents is disregarded. However, field evidence suggests the importance of such lateral flows. For example, Reference Kirillin, Forrest, Graves, Fischer, Engelhardt and LavalKirillin et al. (2015) observed an axisymmetric basin-scale circulation in the ice-covered arctic Lake Kilpisjärvian, driven by heat flux from the shallow region. They notably identified a dense underflow creating upwelling at the lake centre (figure 12).

Figure 12. (a) Buoyancy distribution based on a transect of conductivity–temperature–depth (CTD) profiles with the approximate position of the two current velocity recorders relative to the lakeshore and their vertical coverage marked by the thick dotted line and thin coloured rectangles, respectively. The red and black arrows depict the cross-shore circulation and the central upwelling, respectively. The surface anticyclonic gyre is shown by the pink circles. (b) Combined individual (grey lines) and time-averaged (red and pink lines) vertical profiles for the onshore and alongshore velocity components collected with standard (top) and high-resolution (bottom) ADCPs. The blue line marks the logarithmic fit of the onshore velocity component. Adapted from Reference Kirillin, Forrest, Graves, Fischer, Engelhardt and LavalKirillin et al. (2015).

Through numerical investigation of the energy and heat balance of radiatively driven under-ice convection, Reference Ulloa, Winters, Wüest and BouffardUlloa et al. (2019) quantified the role of lateral flows and their impact on the evolution of the diurnally active convective mixed layer. They demonstrated that a heat balance can be employed to assess the rate of mixed-layer warming and that neglecting differential heating and lateral fluxes would lead to an underestimation of this rate (figure 13). The heat balance in an interior volume $\tilde{{\rm V}}$, shaded as a blue area in figure 13(a), extending from the water–ice interface to the depth of the convective mixed layer $h_{m}$ is

(7.3)\begin{equation} \frac{{\rm d}T}{{\rm d}t} = \frac{1}{h_{m}} (F_{ext} + F_{adv}), \end{equation}

where $F_{ext}$ is the sum of the solar heating rate, depicted by the red arrow in figure 13(a), and the diffusive heat lost towards the ice, green arrow in figure 13(a), and $F_{adv}$ is the summation of advective heat fluxes across the lateral and bottom boundaries of the control volume, shown in yellow arrows in figure 13(a). Essentially, $F_{adv}$ accounts for the contribution of lateral heat transport from the shallows to the heating of the lake interior. Reference Ulloa, Winters, Wüest and BouffardUlloa et al. (2019) demonstrated that the relative importance of the advective term $F_{adv}$ (differential heating) in comparison to the one-dimensional effects of radiative heating and diffusive cooling at the ice–water interface is geometrically controlled, and can be estimated a priori for any lake morphology, incorporating a shallow zone correction to the one-dimensional mixed-layer heat balance. In relative terms, the extent of underestimation of the contribution of advection to the one-dimensional estimate of the warming of the convective mixed layer can be quantified by the geometrical factor $G$:

(7.4)\begin{equation} G = \gamma \left(\frac{\bar{h}}{h_{m}}-1\right). \end{equation}

Here $\bar {h}$ is the average depth of the shallow regions, referring to locations where the basin depth is shallower than the current mixed-layer depth $h_{m}$; and $\gamma$ is the shallow area fraction ($=A_{s} / A_{0}$), where $A_{0}$ is the total lake surface area and $A_{s}$ is the surface area of the shallow regions. Since $G$ consistently remains lower than 0 in the presence of shallows, a one-dimensional estimate of mixed-layer warming is therefore always an underestimate. This factor underscores the significance of both the shallow area fraction and its mean depth relative to the mixed-layer depth in the lake's interior. Lakes with small values of $\gamma$ deliver relatively little heat from the shallows to the interior, resulting in minimal impact on the overall mixed-layer evolution due to lateral transport. For a given shallow area fraction, a greater relative mean depth of the shallows results in a reduced lateral temperature contrast and a proportionally smaller effect on the heating of the lake's interior.

Figure 13. (a) Control volume (blue shading) for interior heat balance using (7.3). (b) Laterally averaged, midday temperature profiles for basins with a sloping bottom (purple) and vertical sidewalls (orange). (c) Temperature in the interior of the convective mixed layer (left axis, circles) and net daily averaged heating rate (right axis, squares). (d) Ratio of daily averaged heat fluxes for the days shown in (b). Adapted from Reference Ulloa, Winters, Wüest and BouffardUlloa et al. (2019).

The low speed of horizontal flows and the typically high latitude of ice-covered lakes call for low Rossby numbers $Ro = U_{h} f^{-1}L^{-1}$, where $f$ is the Coriolis frequency and $U_h$ and $L$ are the characteristic velocity and length scales of the flow. This characteristic allows us to assess how the lateral transport of heat initially constrained by lake bathymetry is also modulated by the Earth's rotation. Reference Ramón, Ulloa, Doda, Winters and BouffardRamón et al. (2021) conducted a set of high-resolution simulations of a bowl-shaped lake under different Rossby numbers to evaluate how the dynamic response depends on lake size and latitude and is ultimately controlled by the Rossby number (figure 14). The latter was defined as $Ro = U/(f R)$, where $R$ is the surface radius of the bowl-shaped lake. In the ageostrophic limit ($Ro=O(10^{-1}$) in figure 14), horizontal density gradients drive cross-shore circulation that transports excess heat to the lake interior, accelerating the under-ice warming there. In the geostrophic regime ($Ro=O(10^{-3}-10^{-2}$) in figure 14), the circulation of the nearshore and offshore waters decouples, and excess heat is retained in the shallows. The flow regime controls the fate of this excess heat and its contribution to water-induced ice melt. Reference Ramón, Ulloa, Doda, Winters and BouffardRamón et al. (2021) proposed a modification of the geometrical factor $G$ that accounts for the Earth's rotation, which for a lake with circular surface area and radius $R$ could be simplified to

(7.5)\begin{equation} G_{Ro} = \left(1- \frac{(R-\ell_{ml})^2}{(R-\ell_{ml} + Ro_R)^2}\right)\left(\frac{h_{Ro}}{h_{m}}-1\right), \end{equation}

where $Ro_R$ is the Rossby radius ($Ro\cdot R$) and $h_{Ro}$ is the average depth of the littoral region within a distance $Ro_R$ from the lake interior. Note that for $Ro_R \geq \ell _{ml}$, $G_{Ro} = G$. Differential heating and the Earth's rotation also set up a circulation characterised by the formation of gyres within the convective mixed layer (figure 14). Under-ice gyre circulation has been also reported from field studies (Reference Likens and RagotzkieLikens & Ragotzkie 1966; Reference Forrest, Laval, Pieters and LimForrest et al. 2013; Reference Kirillin, Forrest, Graves, Fischer, Engelhardt and LavalKirillin et al. 2015). For example, Reference Kirillin, Forrest, Graves, Fischer, Engelhardt and LavalKirillin et al. (2015) observed that the return flow of the axisymmetric basin-scale circulation in the ice-covered arctic Lake Kilpisjärvian turned into a lake-wide anticyclonic gyre, achieving velocities of up to 4 cm s$^{-1}$ (figure 12). Coriolis effects on TS, as exemplified for ice-covered lakes in this section, have also been reported during winter cooling in large lakes at times of $T > T_{md}$. Reference Fer, Lemmin and ThorpeFer et al. (2001), for example, demonstrated that the cross-shore downward current turned into a cyclonic current in Lake Geneva. If we define a $Ro$ number for the littoral region as $Ro = U/(f\ell _{ml})$, low $Ro$ numbers are expected in the littoral regions of large lakes such as the Great Lakes, given their latitudinal location and particularly the extent of these regions ($\ell _{ml} O(10)$ km). An example of the Earth's rotation interacting with differential cooling/heating in such lakes is the formation of the spring and autumn thermal bars.

Figure 14. Snapshots of the thermal and velocity structure for different $Ro$ showing the simulated cross-sectional (a,e,i) temperature (0.02 $^{\circ }$C isotherm spacing) and (b,f,j) radial and (c,g,k) azimuthal velocities (0.002 m s$^{-1}$ isovel spacing) and depth-averaged (d,h,l) azimuthal velocities and flow streamlines at ${\rm depths} < h_{m}$. Radial and azimuthal velocities are positive towards the lake interior and for cyclonic circulation, respectively. Black dashed lines in (d,h,l) show the location of the cross-section displayed in (a–c), (e–g) and (i–k). Results for runs with (a–d) $Ro$ $O(10^{-1})$, (e–h) $Ro$ $O(10^{-2})$ and (i–l) $Ro$ $O(10^{-3})$. Sketches on top conceptualise the under-ice circulation induced by differential heating for the ageostrophic and geostrophic regimes. Adapted from Reference Ramón, Ulloa, Doda, Winters and BouffardRamón et al. (2021).

7.3 Thermal bar

A particularly interesting convective phenomenon, mediated by the nonlinear equation of state of water, occurs in the sloping waters of lakes whose water temperature drops below $T_{md}$, such as those that freeze seasonally. As autumn progresses and nearshore waters cool faster than deeper waters (differential cooling), the very nearshore region reaches temperatures below $T_{md}$, creating what is known as the ‘thermal bar’ (Reference ForelForel 1880): a transition zone characterised by having a well-mixed water column at $T_{md}$, from where a downwelling convective flow fosters the formation of two cross-shore convective cells with opposing circulation – which inhibit the exchange between onshore and offshore waters, with potential impacts in the ecological functioning of aquatic systems, as discussed by Reference Holland and KayHolland & Kay (2003) and Reference TerzhevikTerzhevik (2012). As the surface water cooling progresses, the cross-shore position of the thermal bar moves offshore. The same phenomenon emerges during spring. This time, nearshore waters warm faster than deeper water (differential heating), pushing the thermal bar offshore (see schematic in figure 2). The autumn thermal bar has been reported to be difficult to detect given the small temperature gradient between onshore and offshore regions and atmospheric instability (e.g. Reference HuangHuang 1972). The spring thermal bar, however, could last for 1 month (e.g. Reference Boyce, Donelan, Hamblin, Murthy and SimonsBoyce et al. 1989) and the evolution can be tracked by satellite.

Reference Zilitinkevich, Kreiman and TerzhevikZilitinkevich, Kreiman & Terzhevik (1992) proposed a scaling for the horizontal displacement rate of the thermal bar following a heat balance approach where the propagation speed is estimated as $H_{net,0}/(\bar {s}\rho C_p (T_{md}-T_0))$, with $\bar {s}$ being the bottom slope and $T_0$ the initial lake temperature. Following Reference Malm, Grahn, Mironov and TerzhevikMalm et al. (1993) observations on Lake Ladoga, the typical propagation rate is about 2 km day$^{-1}$ (using $\bar {s} = 10^{-3}$, $H_{net,0} = 200\,{\rm W}\,{\rm m}^{-2}$ and $T_0 = 2\,^{\circ }$C). Various model approaches have been developed to investigate the circulation induced by the thermal bar, gradually increasing in complexity (Reference TsydenovTsydenov 2018, Reference Tsydenov2019). These range from a steady-state temperature and circulation model that includes the Coriolis effect (Reference HuangHuang 1972), to asymptotic unsteady models (Reference Brooks and LickBrooks & Lick 1972; Reference FarrowFarrow 1995a,Reference Farrowb; Reference Farrow and McDonaldFarrow & McDonald 2002; Reference FarrowFarrow 2013a) and 2-D hydrodynamic models that incorporate wind stress (Reference MalmMalm 1995). Moreover, we note that other local processes, such as rivers entering lakes, can also trigger a local riverine thermal bar (e.g. Reference Carmack, Gray, Pharo and DaleyCarmack et al. 1979). Lastly, we acknowledge that a significant number of studies in thermal bar dynamics were led by the scientist A.I. Tikhomirov between the 50's and 80's; a list of his papers written in Russian are listed in Reference TerzhevikTerzhevik (2012).