Latitude and bathymetry modify lake warming under ice

In late winter, solar radiation is the main driver of water motion in ice-covered lakes. The resulting circulation and mixing determine the spatial distribution of heat within the lake and affect the heat budget of the ice cover. Although under10 ice lake warming is often modeled as a one-dimensional vertical process, lake bathymetry induces a relative excess heating of shallow waters, creating horizontal density gradients. This study shows that the dynamic response to these gradients depends sensitively on lake size and latitude —Earth rotation— and is controlled by the Rossby number. In the ageostrophic limit, horizontal density gradients drive cross-shore circulation that transports excess heat to the lake interior, accelerating the underice warming there. In the geostrophic regime, the circulation of the nearand off-shore waters decouple and excess heat is 15 retained in the shallows. The flow regime controls the fate of this excess heat and its contribution to water-induced ice-melt.


Introduction
Convection in rotating fluids is a ubiquitous process redistributing heat in planetary systems (Afanasyev and Zhang, 2018;Heimpel et al., 2016). In the case of the Earth, its spin controls the outer core convective flows responsible for the 'geodynamo' (e.g. Buffett, 2000), as well as the thermohaline oceanic circulation, and the large-scale atmospheric 'cell' circulation (Vallis, 20 2017) and its climate (e.g. Cabré et al., 2017). Yet, as we move from planetary to smaller-scale systems, the importance of rotation in affecting convective processes remains overlooked. This is the case in lakes. Although field observations have been able to characterize different convective processes in lakes , the influence of rotation on the distribution of heat is not well understood. The classical example of convection in lakes is the winter cooling that controls deep mixing (Schwefel et al., 2016). Yet, the presence of sloping lateral boundaries also sets up lateral convective flows that connect 25 rivers entering lakes (e.g. Wells, 2009). For horizontal convection, the scales U and L correspond to the cross-shore velocity of the gravity currents and the distance from the shore to the center of the lake basin, respectively. Gravity currents triggered by horizontal convection are usually "slow" (U of O(cm s -1 ), Fer et al., 2002;Monismith et al., 1990;Wells and Sherman, 2001), so their Rossby numbers are usually less than 1. Here we show that rotational effects on horizontal convection are 35 especially relevant for ice-covered lakes in late winter since they modulate their warming regime.
Lake ice largely insulates the water column from the wind and from most of the atmospheric heat fluxes (Kirillin et al., 2012).
Nevertheless, under thin snow cover conditions, solar radiation penetrates the ice and heats the water adjacent to it (Kirillin et al., 2012). In ice-covered freshwater lakes, water temperature is typically increasing from that of the freezing point at the icewater interface, to the temperature of maximum water density, Tmd (~ 3.98°C), near the bottom of the lake. The density of 40 water can be estimated from its temperature through an equation of state. For temperature values below Tmd, an increase in temperature results in an increase in density. Thus, by heating the uppermost coldest water, it becomes denser, causing it to plunge, initiating convectively-driven turbulence (Farmer, 1975). Radiatively-driven convection tends to both homogenize and heat the upper portion of the water column. The latter can have important implications for winter lake ecology (Farmer et al., 2015) and biogeochemical processes (e.g. Hampton et al., 2017;Karlsson et al., 2013). In addition, lake warming enhances ice 45 melting by increasing the diffusive heat flux from the lake water to the ice. Having a robust mechanistic understanding of the processes controlling the vertical heat transport and distribution is especially relevant to assess the fate of lakes under climate change scenarios. Long term predictions based on one-dimensional (1D, vertical) models indicate that up to 25% of the current seasonally ice-covered lakes are projected to be permanently ice-free by the end of the 21st century (Woolway and Merchant, 2019). The latter might have broad impacts. However, such models and predictions do not integrate the effect of lateral 50 boundaries in the vertical heat distribution.
Lake warming under ice is often modelled as a 1D vertical process (Farmer, 1975;Mironov et al., 2002). In the 1D framework, convective plumes driven by solar radiation lead to the formation of a convective mixed layer (CML, Fig.1a) above the relative quiescent stratified deep layers of the lake and below the diffusive layer at the ice-water interface (Farmer, 1975) (Fig. 1a). On a timescale of days to weeks during late winter, the cumulative insolation increases and the CML becomes warmer and deeper 55 (e.g., Bouffard et al., 2019a;Farmer, 1975), and, in shallow lakes, complete thermal mixing of the water column can occur before the end of the ice-on season (e.g., Huang et al., 2019;Salonen et al., 2014). However, even if the radiative flux through the ice remains spatially uniform over the lake area, the interaction of solar radiation and lake bathymetry can drive horizontal convection due to differential heating. Horizontal density gradients trigger gravity currents that propagate downslope and intrude into deep water at the base of the CML (Cortés and MacIntyre, 2020;Kirillin et al., 2015;Salonen et al., 2014;60 Stefanovic and Stefan, 2002) (Figs. 1b-c). Return flows towards the littoral regions develop at the top of the CML, leading to the formation of cross-shore circulation cells (Fig. 1b). Recently, it has been shown via 2D numerical experiments, that horizontal advective fluxes of heat due to differential heating speed up the warming and deepening of the CML (Ulloa et al., 2019). The degree of departure from the 1D formulation increases as the fraction of the basin volume occupied by the littoral region -the region whose maximum depths are shallower than the CML depth-increases. By assuming that all the excess heat in the littoral region is effectively flushed into the lake interior, the underestimation of the contribution of advection and, thus, of the warming of the CML, can be quantified by a geometrical factor G (Ulloa et al., 2019) where the vertical bars refer to absolute values, Atotal and Ashallow are the surface areas of the lake and the littoral region ( Fig.   1b and Fig. S1), respectively; and ℎ � and hcml are the average depth of the littoral region and that of the CML, respectively (Fig. 70 1d and Fig. S1). However, depending on the lake size -controlling the scale L in Ro-and location (latitude) -controlling f in Ro -, lake dynamics can be also strongly affected by Earth rotation, which affects the heat exchange between shallow and deep waters. There is field evidence on the influence of Coriolis in the circulation of lakes under ice (Forrest et al., 2013;Kirillin et al., 2015;Rizk et al., 2014), where the presence of horizontal cyclonic or anticyclonic gyres have been inferred or measured. These gyres can potentially modulate the advection of heat from the littoral region and modify lake warming. 75 Here we report, through 3D real-scale numerical simulations, that indeed for a fixed radiative heating and bathymetry, the resulting under-ice basin-scale heat distribution strongly depends on rotation, and, thus, potentially on the latitudinal location of a given lake. We found two distinctive dynamical regimes, controlled by the Rossby number, that result in a remarkably divergent distribution of incoming heat beneath the ice. In the ageostrophic regime, Ro O(10 -1 -10 0 ), advection of heat from the littoral region effectively increases the deepening and warming rates of the CML in the lake interior. In the geostrophic 80 regime, Ro ≤ O(10 -2 ), the dynamics of the littoral and the offshore region decouple: the excess heat is retained in the littoral region while the deepening and warming of the CML in the lake interior is well approximated using the 1D framework.

Scenarios
For validation and comparison purposes, our scenarios extend the simulations presented in Ulloa et al. (2019). We specifically 85 kept the same lake bathymetry as in their 2D spectral large eddy simulations (LES), but expanded into three dimensions to obtain a basin with radial symmetry (Fig. 1b). The fluid was initially at rest and linearly stratified, with T = 0°C at the surface and T = Tmd at the maximum depth. We fix T at 0°C at the ice-water boundary and subject the under-ice water to a laterally uniform but time and depth dependent radiative flux, I(t,z) (°C m s -1 ). This is equivalent to a surface heat flux, Q(t,z) = cp ρ0 I(t,z) (W m -2 ), with cp and ρ0 being the heat capacity and reference density of fresh water near its freezing point, respectively. 90 The vertical distribution was modelled using Beer's law (Fig. 1a), whereas the time evolution was prescribed to roughly mimic the daily cycle of solar insolation ( Fig. 1b and see Sect. 2.3). The forcing parameters reflect under-ice conditions in late winter.
By modifying the value of f in our simulations, we investigate three scenarios ranging from weak (Ro O(10 -1 )) to stronger (Ro O(10 -3 )) rotational influence (Table 1 and see Sect. 2.4 and Sect. S1.1 in the supplementary material for the Ro calculations).
This range of Ro spans the expected range of values typical of the varying size and latitudinal distribution of ice-covered lakes on Earth (see Sect. 4). We focus on a rotational regime where convective-plume plunging and mixing occur on time scales much shorter than the inertial period (~ f -1 ), but the temperature gradients that result from differentially heating the shallow littoral and deep interior regions give rise to slower motions that are affected by rotation over longer time scales.

Computational model
Simulations were conducted with the RANS model MITgcm (MIT General Circulation Model, Marshall et al., 1997aMarshall et al., , 1997b and details in http://mitgcm.org). MITgcm is a three-dimensional, non-hydrostatic, z-coordinate, finite volume model that solves the Boussinesq form of the Navier-Stokes equations for incompressible fluids. The momentum equations are discretized on an Arakawa-C grid. The variables are advanced in time using a quasi-second-order Adams-Bashforth time-stepping scheme and preconditioned conjugate-gradient methods are used in the 2D and 3D inversion of the hydrostatic and non-hydrostatic pressure, respectively. The advection terms in the transport equation for temperature were discretized with the non-linear 3rd 105 order DST (direct space-time) with flux limiter. A modified nonlinear UNESCO equation of state (Jackett and Mcdougall, 1995) and the non-hydrostatic capabilities of MITgcm were used. Horizontal and vertical viscosities were parameterized with the isotropic 3D Smagorinsky approach (0.0005 coefficient). Background grid-dependent lateral viscosities were set to 0.002 (equivalent to a lateral eddy viscosity, υh, of ~ 8 × 10 -4 m 2 s -1 for grid resolution in the x-direction, Δx, of 0.9 m and a time step, Δt, = 0.5 s). Background vertical viscosities, υz, were set to 10 -6 m 2 s -1 . Background horizontal and vertical diffusivities 110 for heat, Kh and Kz, were set to 10 -5 m 2 s -1 and 1.4 × 10 -7 m 2 s -1 , respectively. No-slip conditions were applied at the lateral walls and the bottom. Though we could potentially have done this initial study using 3D LES, our interest extends to larger lakes over seasonal time scales and this setup provides and validates a starting point for such studies. MITgcm was originally built for ocean modelling, but has been applied satisfactorily to study lake hydrodynamics (e.g. Cimatoribus et al., 2018;Dorostkar et al., 2017). We further validate MITgcm against spectral LES (Ulloa et al., 2019) of radiatively heated ice-covered 115 lakes (see Sect. S1.5 in the supplementary material).

Under ice hydrodynamic and transport model
The lake bathymetry is conical and axisymmetric ( Fig. 1), with a radius R = 114 m and depths, D, increasing towards the lake interior in the radial direction, r, as in Ulloa et al. (2019): where H = 27.1 m is the maximum height of the basin, dmin (= 0.56 m) is the minimum depth in the littoral region, and r0 (= 57.2 m) and β (= 22.9 m) are topographic (slope) parameters. The physical domain was discretized using a horizontally uniform Cartesian grid (Δx = Δy = 0.9 m) with vertically variable thickness (Δz). Δz increases with depth from Δz = 0.05 m within the first 2 m to a bottom cell of 0.3 m. In our simulations, temperature and pressure are the only scalars affecting density. The lake was initially quiescent with horizontal isotherms following a linear thermal stratification, with T = 0 ºC at the ice-water interface and T = Tmd at the bottom. Boundaries are prescribed as no slip and adiabatic, except at the ice-water interface. The ice-water boundary is modeled as a rigid lid with a constant temperature of 0 ºC. The radiative heat flux through the ice-water interface is specified as a sinusoidal daily cycle that attenuates with depth as in Ulloa et al. (2019): (3) Here the time t is expressed in days, I0 (= 1×10 -5 ºC m s -1 ) is the water surface radiative forcing, F(t) = sin(2πt) during the day 130 (t < 0.5) and zero otherwise, and 1/λ (= 2.5 m -1 ) is the depth scale for light attenuation. The order of magnitude of I0 and the value for λ selected, are representative of late-winter conditions in moderately clear lakes (Ulloa et al., 2019). For visualization purposes only, our results are shifted 0.25 days so the peak in the radiative heat flux matches midday (Fig. 1).
The model was first used to simulate rotational effects as in Ulloa et al. (2019), with a Coriolis frequency f = 1.26×10 -4 s -1 135 (~60° N). This corresponds to run 1 in Table1. For this value of f, gravity currents in our simulations reach velocities of ~5 mm s -1 consistent with the reported O(10 -3 -10 -2 ) m s -1 radial velocities under ice (Forrest et al., 2013;Kirillin et al., 2015;Rizk et al., 2014). To analyze the effect of rotation in the lake circulation and in the warming of the CML, two additional simulations were conducted where we increased rotational effects by increasing f up to two orders of magnitude (runs 2 and 3 in Table 1). To analyze bathymetric effects (differential heating), an additional simulation was conducted (reference 140 simulation) where forcing was kept as in run 1, but the bathymetry was modified to obtain a cylinder of depth D = H. Each run spans 12 radiative cycles (12 days). This number of cycles was long enough to expose and analyze the effect of rotation and bathymetry on lake warming under ice.

Rossby number
Since we are interested in evaluating the advection of heat from the littoral to the lake interior, the surface radius, R, was 145 selected as the characteristic length scale in the calculations of Ro. Ro was calculated using the maximum radial velocity in the littoral region, Urs-max, as the characteristic velocity scale (see details in Sect. S1.1 in supplementary material).
(4) R /Urs-max is the nominal time required for a gravity current at speed Urs-max to reach the center of the lake. When R /Urs-max > f -1 , gravity currents are affected by Earth rotation. 150

Differential heating, warming of the CML and heat balance
The magnitude of differential heating was quantified by comparing CML temperatures, Tcml, in the littoral region and at the lake center. We compared temperatures at 1 m depth (inside the CML in all runs and days > 1) at profile P1 in the littoral region (at a distance β from the outer edge and ~91 m from the lake center) and at profile P2 at the lake center (Fig. 1c). The average Tcml at the lake interior was calculated as the average mixed-layer temperature inside an interior control volume ( Fig.1d) of radius ~25 m (as in Ulloa et al., 2019). This control volume was also used to calculate the importance of the advection of heat from the littoral region in the mixed-layer heat balance at the lake interior (see Sect. 3).

Mixed layer warming and differential heating for different Rossby numbers
Radiatively-driven convection produces a CML with uniform temperature in each of the Ro scenarios simulated (see the 160 temperature profiles in Figs. 2a-c). This CML warms and deepens during the 12 daily cycles simulated, with higher temperatures reached as Ro increases (red line in Figs. 2a-c). Fig. 2D shows the time history of temperature in the CML, Tcml, in the lake interior (average value in the interior control volume CV, Fig. 1d) for the three Ro tested. The rate of warming of the interior mixing layer, dT/dt, increases with increasing Ro, with the net warming rate for Ro O(10 -1 ) 23% higher on average than that for Ro O(10 -3 ) (compare the red and blue lines in Fig. 2d). 165 The development and maintenance of a temperature contrast between the littoral and the lake interior is also affected by Ro.
To quantify the differential heating, we compared the temperature difference at 1m depth (inside the CML), ΔTh,z=1m, between two profiles: one in the littoral region (P1, Fig. 1d) and one at the lake center (P2, Fig.1d). The temporal evolution of ΔTh,z=1m is remarkably sensitive to Ro (  Fig. 2e), the rate of change of ΔTh,z=1m in time also decreases, which may also suggest the establishment of a quasi-steady state.
The simulations illustrate that horizontal temperature and density gradients increase by about an order of magnitude as Ro decreases from O(10 -1 ) to O(10 -3 -10 -2 ) and that the increasing cross-shore contrast over successive days indicates that excess 180 heat is effectively retained in the littoral region. As a consequence, the horizontal lake circulation becomes modified by rotational effects.

Radial and azimuthal lake circulation
The characteristic signal of differential heating under ice, with temperature at a given depth increasing towards the littoral region and the formation of cross-shore circulation cells, is clearly observed for Ro O(10 -1 ) (Figs. 1b and 3a-b). Because of 185 rotation, however, the gravity currents and the return flows are deflected towards the right in the northern hemisphere. Driven https://doi.org/10.5194/hess-2020-471 Preprint. Discussion started: 14 October 2020 c Author(s) 2020. CC BY 4.0 License.
by the offshore-directed gravity currents, a cyclonic gyre develops in the lake interior (positive azimuthal velocities, uθ, in Fig.   3c); while the onshore return flow drives an anticyclonic gyre in the littoral region (negative uθ in Fig. 3b). The gyres are governed by a three-way balance of horizontal pressure gradient, Coriolis and centripetal forces (cyclogeostrophic balance, Sect. S1.2 in the supplementary material). The signature of the double-gyre circulation is distinguishable from the CML-190 averaged azimuthal velocities shown in Fig. 3d. The strength of the gyres is maximal immediately after the peak in solar radiation, when gravity currents flush the littoral region. Although weaker in strength, the double gyre circulation does not fully dissipate at night (see video S1 in the supplementary material). The weakly-dissipative nature of the gyres is in agreement with estimates of the Ekman damping time scale, tEk = (fEk 1/2 ) -1 , where Ek ( = υz f -1 hgyre -2 ) is the Ekman number, υz is the vertical kinematic viscosity and hgyre is the gyre thickness. For υz = 10 -6 m 2 s -1 (see Sect. 2) and hgyre between 3 and 7 m (Fig.  195 3c), tEk ranges between 3 and 7 days, which is longer than the period of the radiative forcing (less than 1 day). Flow convergence in the lake interior in the radial direction and rotation upwell isotherms at the bottom of the CML. This is illustrated by the dome shape of the isotherms at the bottom of the CML in Fig. 1b and Fig. 3a and by the positive horizontal density anomaly at the center of the lake in Fig. 1c. Isotherm upwelling explains the presence of a weaker stratification immediately below the CML at the center of the lake (red profiles in Fig.1a and Figs. 2a-c). 200 Decreasing Ro from O(10 -1 ) to O(10 -2 ) alters the circulation described above. The bowl shape of isotherms above 10m depth in Fig. 3e shows that temperature increases towards the littoral region. However, the radial velocity in Fig. 3f shows that the downslope cross-shore gravity currents are strongly suppressed. There, a balance is set between the Coriolis force and the pressure gradient (geostrophic balance, Eq. (S1) in Sect. S1.2 in the supplementary material). This geostrophic balance results in the formation of an anticyclonic gyre that dominates the circulation within the CML (Figs. 3g and 3h) and restricts the cross-205 shore transport. For Ro O(10 -3 ), isotherms in the littoral region get closer together (Fig. 3i) and the cross-shore transport stops (see radial velocities in Fig. 3j). As the area where cross-shore density gradients develop becomes smaller, the spatial extension of the anticyclonic gyre becomes confined to the shallow littoral region (Figs 3k and 3l). Thus, water circulation under ice is strongly affected by Earth rotation. Once a geostrophic balance is set, the circulation in the CML evolves from a combination of radial (counter-rotating cells) and azimuthal (gyres) motions to a circulation that occurs preferentially in the azimuthal 210 direction. This change in the lake circulation could have consequences for the heat balance of the lake interior.

Earth rotation modifies the interior mixed-layer heat balance
We examine the impact of rotation on the heating rate of the interior mixed layer. For this, we consider an arbitrary control volume away from the littoral region (gray cylinder in Fig. 1d) with a surface area ACV and extending vertically from the icewater interface to a depth hcml. Neglecting interior diffusive fluxes, the heat balance in this control volume is (Ulloa et al.,215 2019): where Fext is the sum of the solar heating rate (yellow arrow in Fig. 1d) and the diffusive heat lost towards the ice (blue arrow in Fig. 1d), and Fadv is the sum of advective heat fluxes across the lateral and bottom boundaries of the control volume (red and orange arrows in Fig. 1d). The daily-averaged Fext and Fadv are calculated as: 220 Here, the angle brackets denote both a volume and a daily averaging in the CML, zcml is the height of the base of the CML, H is the total height of the basin (Figs. 1d and S1) and Kz is the vertical thermal diffusivity. Equation (7) integrates both vertical 225 and radial advective fluxes, Scv being the surface of the control volume with outward facing unit normal vector � and u being the vector velocity.
The relative contribution of advection to the total heat fluxes is, thus, δadv = Fadv/(Fext + Fadv). We compared this ratio, for the three Ro examined, as deviations from the contribution calculated in a reference numerical experiment where the 1D approach for deepening and warming of the CML applies, that is Δδadv = δadv -δadv-ref. This reference run is identical to run 1 in Table 1, 230 except for the lake bathymetry: the conical shape in the reference run was replaced by a cylinder of height H and of surface area equal to Atotal. By choosing a cylinder, we eliminated the effect of differential heating. Fig. 4a shows Δδadv for the three  (Figs. 1d and S1). Adapted to a circular surface area, this corresponds to: Dots in Fig. 4b show that only the calculated Δδadv values for Ro O(10 -1 ) are consistent with predictions using the G factor 245 (red dots in Fig. 4b), while the contribution of advection to the heat balance at the lake interior would be strongly over-predicted for Ro O(10 -3 -10 -2 ) (blue and yellow dots in Fig. 4b). Given that rotation controls the advection of heat to the lake interior, the https://doi.org/10.5194/hess-2020-471 Preprint. Discussion started: 14 October 2020 c Author(s) 2020. CC BY 4.0 License. characteristic length scale to define the littoral region is, however, better characterized by the Rossby radius (RoR = Ro×R, Sect. S1.1 in supplementary material) whenever RoR < Lshallow. This leads to a new definition of G, here called GRo, expressed as a function of the Rossby radius and, specifically for our bathymetry, as (see details in Fig. S1 and Sect. S1.3 in the 250 supplementary material): where ℎ ����� is the average depth of the littoral region within a distance RoR from the lake interior (Fig. S1). Note that for RoR ≥ Lshallow, GRo = G. By using GRo when RoR < Lshallow, the trends in the contribution of advection of heat from the littoral region to the heat balance at the lake interior are better represented, especially for the scenario with the lowest Ro (blue squares in 255 Fig. 4b).

Advection of heat from the littoral regions decreases with decreasing Ro (big lakes and/or high latitude lakes)
In this study, we show that, subjected to the same radiative forcing, the under-ice warming of lakes is strongly modulated by the combination of lake bathymetry and latitude (i.e., intensity of Earth rotation). Once solar radiation is able to penetrate the 260 ice and heat the water below it, under-ice convection becomes the main driver of lake circulation and mixing (Kirillin et al., 2012). Convective cells impinge on the stratified layer below, entraining water from below and deepening the CML. This process is often viewed as a one-dimensional process (Farmer, 1975). However, the presence of a bathymetry with sloping boundaries implies that as the CML becomes deeper in the lake interior, convective deepening is constrained by the bottom in the shallower regions. With a smaller volume to heat up, these shallower areas warm faster than the lake interior and a 265 horizontal density gradient develops. If this excess heat is advected towards the lake interior, the warming and deepening of the CML in the lake interior speeds up (Ulloa et al., 2019;and Fig. 2a-d). The advection of heat from the littoral region to the lake interior is, however, constrained by Earth rotation. Our findings show that the rotation regime plays a fundamental role on the heating rate of the CML and in the cross-shore and vertical heat distribution (Figs 2-3). As the rotation increases, the vertical and radial heat exchanges are dramatically inhibited (Figs. 3-4a). In fact, the warming and deepening rates of the 270 interior CML for the scenario with the smallest Rossby number (O(10 -3 )) examined, match the results of the reference simulation, which lacks differential heating and lateral neat heat transport. For this scenario, the density difference between the littoral region and the lake interior increases over time until the end of the simulation (Figs. 2e). This indicates that heat is effectively retained in the littoral region (see also Fig. 3k). Consistent with the results of our simulations, the new GRo (Eq. (9)) geometrical parameter that adopts RoR as the representative length scale when RoR < Lshallow, predicts a decrease in the lateral 275 advection of heat as Ro decreases.

Conceptual model for lake circulation
Differential heating and Earth rotation set up a circulation characterized by the formation of gyres within the CML (Fig. 3 and   Fig. 5). The double-gyre circulation (Fig. 5a) for Ro O(10 -1 -10 0 ) is consistent with past field-based work on ice-covered lakes. Forrest et al. (2013), conducted CTD measurements mounted on an autonomous underwater vehicle and observed a density 280 distribution in the interior of Pavilion Lake (50º N, R~0.4 km and Atotal = 5 km 2 ) that suggested the existence of a cyclonic gyre for Ro O(10 -1 ) (recalculated from Forrest et al. (2013) using L = R). Kirillin et al. (2015) conducted CTD profiles across lake Kilpisjärvi (69º N, R~1.5 km and Atotal = 37 km 2 ) and ADCP measurements at the littoral region when Ro was O(10 -1 ). With the former data, they observed warm "upwelling" at the lake center, which could be indicative of a strong cyclonic gyre in the lake interior; and with the latter, they measured radial velocities of ~3-5 cm s -1 and the presence of anticyclonic gyre in the 285 littoral region with azimuthal velocities of 2-4 cm s -1 . A double-gyre circulation was also reproduced in numerical simulations (Huttula et al., 2010) of early winter conditions in Lake Pääjäarvi (61° N), when under-ice circulation was dominated by the input of heat from the sediment instead of by radiatively-driven convection. Note, however, that for the ageostrophic regime to display a double-gyre circulation, the cross-shore cell circulation should exist over a significant fraction of an inertial period (~ f -1 ). In our simulations, I0 is close to its daily maximum for ca. 6 h, but cloudiness or mountain shading could decrease this 290 duration in real settings. When Ro ≥ 10 -1 , the horizontal heat transport is then accomplished by the ageostrophic components of the flow (downslope gravity currents). This cross-shore circulation might be considered analogue to the atmospheric Hadley cells (Fultz et al., 1959).
As Ro decreases, the circulation in the lake becomes geostrophic (Fig. S3), the daily cross-shore cells are rotationally suppressed and, because of the horizontal density gradient, an anticyclonic gyre (in the northern hemisphere) with velocities 295 decreasing with depth must develop (Eq. (S1) and Figs 5b). With decreasing Ro, the spatial extent of this anticyclonic gyre is progressively confined to the littoral region (Fig. 5c). To our knowledge, this circulation has not been reported in the field at times of under-ice radiatively-driven convection. The lake-wide anticyclonic circulation in Fig. 5b would, however, be consistent with the inferred lake-wide anticyclonic gyre reported by Rizk et al (2014) at a time when circulation in lake Pääjäarvi was dominated by a lateral gradient in the heat flux from the sediment and Ro was O(10 -3 -10 -2 ). The anticyclonic 300 gyre circulation is also consistent with the Rossby wave regime reported in laboratory studies (rotating cylinder and annulus) mimicking the mid-to-high-latitude atmospheric circulation (Fultz et al., 1959). Within the Rossby regime, vortices and waves traveling cyclonically develop (Fig. 3h) and as Ro decreases, the wave lengths decrease and the gyre circulation is concentrated into jets that meander in the radial direction and could finally break. The presence of waves and/or vortices as in the scenario with Ro O(10 -2 ) (Fig. 3h) is typical of transitional regimes (Fultz et al., 1959) and when they develop, the center of the 305 anticyclonic gyre is not static in time but fluctuates laterally (video S1). This lateral displacement favors the existence of pulses of gravity currents from the littoral region to the lake interior, leading to a higher contribution of advection in the heat balance of the lake interior. The lateral fluctuations of the gyre cannot be accounted for in the geometrical factor GRO and explains the https://doi.org/10.5194/hess-2020-471 Preprint. Discussion started: 14 October 2020 c Author(s) 2020. CC BY 4.0 License. under prediction of the contribution of advection for Ro O (10 -2 ) in Fig. 4b (note that the yellow squares in Fig. 4b are above the 1:1 line). 310

Implications
This study reveals that for a given radiative forcing and latitude, as a result of differential heating, a temperature mooring deployed at the deepest point of an ice-covered lake will record a higher warming rate if shallow regions contribute substantially to the lake volume. Thus, for a given lake surface area, a fjord-type lake (more vertical walls) would tend to warm at a slower rate than a bowl-shape lake with gentle slopes. However, if we could virtually move these two same lakes to higher 315 latitudes, that mooring could potentially record the same warming rate in the two lakes. That is because the contribution of the littoral region to the warming of the lake interior under ice is strongly dependent on its latitudinal location and size, both of which condition the importance of Earth rotation (parameterized by Ro) in driving lake water circulation. The range of Ro values investigated here lies within the natural range of variability of Ro that characterize lake circulation driven by under-ice differential heating in lakes on . This means we are covering ice-covered lakes from the Peruvian Andes (~9-320 13°S) and the Tibet Plateau (~30° N) to high-latitude polar regions. The bias in the potential for ice cover formation (e.g. Sharma et al., 2019) and the distribution of lakes on Earth towards mid to high latitudes (> 45° N, Fig. 5b) indicates that the ageostrophic regime with Ro O(10 -1 ) should be common under ice as shown by the mean (red lines) values in Fig. 5d. The geostrophic regime, Ro < O(10 -1 ), lies outside the interquartile range (gray bars in Fig. 5d), but is also naturally occurring in lakes in almost all latitudes (Fig. 5d). The retention of heat in the littoral region in the geostrophic regime could potentially 325 reinforce ice melting near the shoreline. This is suggested in our simulations where with Ro O(10 -3 ) the diffusive flux towards the ice is on average 23 % larger in the littoral region than at the center of the lake (Profiles P1 and P2 in Fig. 1c). However, heat retention in the littoral region does not necessarily imply an increase in the diffusive flux towards the ice. The diffusive flux, Kz ∂T/∂z, depends not only on the temperature in the shallow region, but also on the thickness of the diffusive boundary layer, which determines the value of ∂T/∂z. Due to the retention of heat in the littoral region, water temperature there could 330 potentially reach values ≥ Tmd. This would lead to the development of a stable stratification in the littoral region and to the suppression of convection that, in contrast, would continue in the lake interior.
Depending on the magnitude of the forcing, and the thickness of the convective mixing layer hcml, the Rayleigh number, which quantifies the strength of buoyancy over diffusion, Ra = gαI0 hcml 4 /(υz Kz 2 ), can easily reach magnitudes of ∼O(10 15 ) or even higher if we consider the horizontal convective length as the relevant scale. Here α is the thermal expansion coefficient. These 335 magnitudes are difficult to achieve in laboratory and numerical experiments (King et al., 2009). There have been significant advances in in-situ measurements over the last years that facilitate more robust measurement of turbulence and mixing properties (Bouffard et al., 2019b). Thus, ice-covered lakes might provide unique observations to test theoretical advances on the heat and mass transfer in rapidly rotating convection in fluid systems of planetary scale (e.g. Julien et al., 2012;King et al., 2009). Large and high-latitude ice-covered lakes provide, therefore, unique opportunities to investigate strong convective 340 regimes at low (≪ 1) Rossby and Ekman numbers and with Prandtl numbers ≈ 10.

Conclusions
Penetrative radiation affects the under-ice melting rate by regulating the water-to-ice heat transfer and is often modeled as a one-dimensional process. Our results show, however, that lake bathymetry and latitudinal location also affect the rate of warming of ice-covered lakes. Lake bathymetry induces a relative excess heating of shallow waters. The transport of this 345 excess heat to the lake interior depends on the intensity of Earth rotation and determines lake warming rates and the horizontal distribution of heat. This study stresses that accounting for the shape and size of the lake basin and its latitudinal location is essential for global estimations of lake ice cover that take into account the warming rates and the distribution of heat in the water column.

Code availability 350
The MITgcm code is publicly available at http://mitgcm.org

Data availability
The data displayed in the figures, temperature profiles at P1 and P2, time series of cross-sectional data and the input files to run the simulations are available for download at Zenodo (http://doi.org/10.5281/zenodo.4027393). The complete time series of the 3D model outputs used in this study can be requested from the corresponding author. 355

Author contribution
C.L.R has run the simulations, conducted the bulk analysis and led the writing of the manuscript. The other co-authors have equally contributed to the study design, interpretation of the data and editing of the manuscript.   Run 2 1.26 ˣ 10 -3 0.019 ± 0.011 Conical (Fig. 1B) Run 3 1.26 ˣ 10 -2 0.001 ± 0.001 Conical (Fig. 1B)