Incorporating wind sheltering and sediment heat flux into 1-D models of small boreal lakes: a case study with the Canadian Small Lake Model V2.0

Lake models are increasingly being incorporated into global and regional climate and numerical weather prediction systems. Lakes interact with their surroundings through flux exchange at their bottom sediments and with the atmosphere at the surface, and these linkages must be well represented in fully coupled prognostic systems in order to completely elucidate the role of lakes in the climate system. In this study schemes for the inclusion of wind sheltering and sediment heat flux simple enough to be included in any 1-D lake model are presented. Example simulations with the Canadian Small Lake Model show improvements in surface-wind-driven mixing and temperature in summer and a reduction of the bias in the change in heat content under ice compared with a published simulation based on an earlier version of the model.

1 Introduction

The surface roughness of a small lake or reservoir nearly always contrasts sharply with that of its terrestrial surroundings and will thus be associated with a different atmospheric boundary layer than would be found on shore. The new boundary layer does not develop instantaneously over the entire lake: as the atmospheric flow encounters a sudden change (generally a decrease) in roughness at the shoreline, an internal boundary layer (IBL) develops with a transition region whose properties, including key lake mixing parameters like wind stress, vary with fetch. In fact, for sufficiently small lakes this transition region might occupy the entire surface area: wind stress varies with downstream distance and an equilibrium boundary layer never forms. In addition, for elongated lakes and reservoirs, the net response to wind forcing will vary depending on whether the wind is along or across the primary axis. The impact of fetch-varying surface wind speed due to the presence of IBLs on lake evaporation (e.g. Venäläinen et al., 1998) and gas flux (e.g. Kwan and Taylor, 1994) has been considered, but the effect of a sudden change in aerodynamic roughness on lake mixing does not appear to have been thoroughly examined in the literature.

On the other hand, the relationship between lake area and mixing has been investigated. Mazumder and Taylor (1994) related both lake size and water clarity to epilimnion depth and found linear relationships between epilimnion depth and log fetch (where fetch is taken as the square root of lake surface area) for different transparency classes. Fee et al. (1996) also investigated the effects of lake size and water clarity on mixed layer depths and found that lake size was the more important factor for determining mixed layer depth in a set of Canadian shield lakes, though transparency modulated this response for smaller lakes. The influence of fetch on epilimnion depth cannot be unambiguously attributed to variations in surface stress, since lake size influences seiching and thus shear production of turbulence at the thermocline, which can also contribute to mixed layer deepening (e.g. Gorham and Boyce, 1989), but these studies do hint at the possible importance of IBLs and fetch-varying surface stress for lake hydrodynamics. Most 1-D lake models assume that turbulent flux exchange with the atmosphere (heat and moisture as well as momentum) takes place under equilibrium conditions such as is assumed, for example, when employing Monin–Obukhov similarity theory.

Another physical process, particularly relevant for ice-covered lakes, is the sediment heat flux. This flux arises because lake sediment, especially in shallow littoral zones, warms under the influence of penetrating radiation (as well as thermal contact with warm water) during summer and releases this heat during the ice-cover season. Rizk et al. (2014) report results from numerous studies that found sediment heat fluxes of the order of a few Wm⁻², generally larger during early winter and tapering off as winter progresses. This is insignificant compared with other energy fluxes during the open water season and probably within the observational uncertainty of meteorological forcing (not to mention the uncertainty in process parameterization) in numerical modelling studies. During the ice-cover season, however, energy fluxes into the lake are small, and this source can become important for many applications. For example, sediment heat fluxes have been linked to basin-scale circulations in ice-covered lakes (Kirillin et al., 2012; Rizk et al., 2014), which can strongly influence the thermal structure and distribution of dissolved oxygen and nutrients in deep waters. The impact of sediment heat flux on surface conditions such as ice thickness and ice phenology is generally assumed small, except perhaps for very shallow lakes, though this has not been thoroughly examined in the literature.

The Canadian Small Lake Model (CSLM; MacKay, 2012; MacKay et al., 2017) is a 1-D thermodynamic lake scheme developed for coupling within global and regional climate and numerical weather prediction systems. This model computes the complete nonlinear surface energy balance in a thin layer (5 cm) and then solves the heat equation on a uniform finite-difference grid throughout the column. Shortwave radiation extinction follows Beer's law for both visible and near-infrared bands. A diurnal surface mixed layer is simulated based on an integrated turbulent kinetic energy (TKE) approach, with a variety of well-known sources and sinks of TKE parameterized. The seasonal thermocline arises naturally as a result of the daily excursions of the surface mixed layer. Congelation (i.e. black) ice forms when the energy balance in a layer is sufficiently negative to cool it below 0 ^∘C. Snow (i.e. white) ice forms when the weight of the overlying snowpack is sufficient to crack the ice and allow lake water to flood a layer of snow. The snow itself is represented with the complete snowpack parameterization component of the Canadian Land Surface Scheme (Verseghy and MacKay, 2017). Both fractional ice cover and fractional snow on ice are permitted. The current formulation of the model does not account for wind sheltering or sediment heat fluxes. To begin to address these shortcomings, this study proposes two new schemes which are simple and flexible enough to be easily incorporated into any 1-D lake model. Below these schemes are fully developed from theoretical considerations and some example CSLM simulations demonstrating their impacts are presented.

2 Fetch-varying wind stress and mixing in small lakes

In the analysis that follows we consider for simplicity the case of neutral atmospheric stability with no surface heat or vapour flux. Under steady-state conditions over horizontally homogeneous surfaces we expect the wind profile to be given by the familiar logarithmic form:

$$\stackrel{\mathrm{‾}}{u}\left(z\right)={\displaystyle \frac{u_{\ast }}{k}}\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right),$$

where $\stackrel{\mathrm{‾}}{u}$ is the mean wind profile, k is the von Karman constant, z₀ is the surface roughness, and u_∗ is the surface friction velocity of the air, which is related to surface stress by $\mathit{\tau }=\mathit{\rho }{u}_{\ast }^{\mathrm{2}}$, where ρ is the density of air. Thus upstream of a shoreline over a (likely vegetated or urbanized) terrestrial landscape we find (adopting the notation of Jensen, 1978)

$$\stackrel{\mathrm{‾}}{u}\left(z\right)={\displaystyle \frac{u_{-}}{k}}\ln \left({\displaystyle \frac{z}{z_{-}}}\right),$$

and far downstream of the shoreline under equilibrium conditions over a very large lake we expect

$$\stackrel{\mathrm{‾}}{u}\left(z\right)={\displaystyle \frac{u_{+}}{k}}\ln \left({\displaystyle \frac{z}{z_{+}}}\right),$$

where the − and + subscripts refer to upstream and far downstream quantities respectively. In what follows the ratio of the upstream-to-downstream roughness lengths, characterized by

$$M=\ln \left({\displaystyle \frac{z_{-}}{z_{+}}}\right),$$

is an important parameter governing the response of the surface stress due to the sudden change in surface roughness at the shoreline. Over lakes the roughness elements are provided by surface waves, and thus roughness is not static but rather a function of wind speed; nevertheless, we shall consider it O(10⁻³ m) in what follows. Common terrestrial surfaces are forest ($z_{-}=\mathrm{1.0}\mathrm{m}$ nominally), shrubland ($z_{-}={\mathrm{10}}^{-\mathrm{1}}\mathrm{m}$ nominally), and grassland ($z_{-}={\mathrm{10}}^{-\mathrm{2}}\mathrm{m}$ nominally), which yield M=6.9, 4.6, and 2.3 respectively.

The development of IBLs due to discontinuities in surface roughness over rigid surfaces has been studied for many years. This might not have been the case: as pointed out by Jensen (1978), in the planetary boundary layer the downstream equilibrium surface stress (${\mathit{\tau }}_{+}=\mathit{\rho }{u}_{+}^{\mathrm{2}})$ under realistic conditions is different from the upstream surface stress (${\mathit{\tau }}_{-}=\mathit{\rho }{u}_{-}^{\mathrm{2}})$ by at the very most a factor of 3. If the transition from upstream to downstream conditions were simply monotonic, then Jensen speculates that there would be little interest in the development of IBLs. For most purposes a simple interpolation rule could be developed to distribute wind stress over a surface, if one were to bother at all. However, experience shows that this transition is not so simple. For atmospheric flow from a rough to a smooth surface, both experimental and theoretical results indicate that surface stress drops suddenly across the roughness transition before asymptotically approaching its new equilibrium value (e.g. Garratt, 1990). On the other hand, for flow from a smooth surface to a rough surface, the surface stress initially “overshoots” the final equilibrium value before slowly asymptoting to the final value, though this situation seems less relevant for our purposes as most lakes are situated within environments that are usually rougher than the water surface.

Figure 1 Theoretical IBL depths (a) and surface stress (b) as a function of fetch for the rough-to-smooth transition case of Bradley (1968). Solid curves – Panofsky and Townsend (1964) approach; dashed curves – Jensen (1978) approach with A=1 (thin dash) and A=2 (thick dash) as indicated. Error bars in (b) represent the range of observed values from Bradley (1968), and surface stress values are normalized by the upstream value.

Panofsky and Townsend (1964) developed a simple, approximate analytical model describing this phenomenon, and Bradley (1968) showed that this model qualitatively captures the observed behaviour over the surfaces he examined, though he did find that the observed transition region was smaller than predicted. Panofsky and Townsend solve for a surface stress parameter given by (again adopting the notation of Jensen, 1978)

$$\begin{array}{}\text{(1)}& S=\left(u_{-}-u_{\mathrm{0}}\right)/u_{-},\end{array}$$

where u₀=u₀(x) is the surface friction velocity (of air) and is a function of fetch, x, with x= 0 representing the location of the roughness discontinuity (in our case the shoreline). Panofsky and Townsend argue that S can be approximated by

$$\begin{array}{}\text{(2)}& S\approx {\displaystyle \frac{M}{\ln \left(\frac{d}{z_{+}}\right)-\mathrm{1}}},\end{array}$$

where d is the depth of the IBL. They develop a relation between boundary layer depth d and fetch which can be solved iteratively and substituted into Eq. (2) to solve for S as a function of fetch. From Eq. (1) we can then evaluate the ratio of the over-lake wind stress to the upstream terrestrial surface wind stress as

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{0}}}{{\mathit{\tau }}_{-}}}={\left(\mathrm{1}-S\right)}^{\mathrm{2}}.\end{array}$$

Jensen (1978) suggests a different approach to this problem. He argues that the depth of the IBL can be solved from

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{d}{z_{+}}}\left(\ln \left({\displaystyle \frac{d}{z_{+}}}\right)-\mathrm{1}\right)=A{\displaystyle \frac{x}{z_{+}}}-\mathrm{1},\end{array}$$

for some constant A, and that the wind stress ratio becomes

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{0}}}{{\mathit{\tau }}_{-}}}={\left(\mathrm{1}-S^{\prime }\right)}^{\mathrm{2}},\end{array}$$

where

$$S^{\prime }={\displaystyle \frac{M}{\ln \left(\frac{d}{z_{+}}\right)}}.$$

Before examining the results of these models for forest, shrubland, and grassland environments, we first revisit the experiments of Bradley (1968) in order to gain some insight into the models. Bradley examined the flow transition from a rough surface (created with a mesh of wire spikes: $z_{-}=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$ m) to a smooth surface (runway tarmac embedded with sand and small pebbles: $z_{+}=\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m), and vice versa, in an experiment at an airfield in New South Wales, Australia, yielding M=4.8 – approximately the value for our shrubland environment. IBL depths and stress ratios τ₀/τ₋ as a function of fetch are indicated in Fig. 1 for the Panofsky–Townsend model (solid) and Jensen (A=1 thin dash; A=2 thick dash) models. Jensen (1978) reports that for the smooth-to-rough transition case, setting A=1 reproduces the observed Bradley data very well and yields results that are nearly identical to the theoretical approach of Rao et al. (1974) which is based on a second-order turbulence closure scheme. However, Fig. 1 suggests that for the rough-to-smooth case, setting A=2 may be equally appropriate based on the limited data from Bradley (the maximum fetch was only about 12 m). In either case the results from Jensen appear to better represent the observed data than does the Panofsky–Townsend approach. In addition, the free parameter A could be adjusted should more observed data become available, giving this approach some appeal.

The final downstream equilibrium surface stress τ₊ (or equivalently u₊) is not discussed by Panofsky and Townsend, but Jensen finds that

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{{\mathit{\tau }}_{+}}{{\mathit{\tau }}_{-}}}={\left[\mathrm{1}-{\displaystyle \frac{M}{\ln \left(Ro\right)}}\right]}^{\mathrm{2}},\end{array}$$

where Ro is the downstream surface Rossby number given by

$$Ro={\displaystyle \frac{G}{f{z}_{+}}}.$$

G is the geostrophic wind speed and f is the Coriolis parameter. Thus for a given value of Ro, Eq. (6) can be used to recast results of the Panofsky–Townsend and Jensen models in terms of the approach to equilibrium of surface stress (τ₀/τ₊) by simply dividing Eqs. (3) or (5) by Eq. (6). For midlatitude lakes, we take $f={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ and G=10 ms⁻¹. Jensen's primary interest was the planetary boundary, and he considered equilibrium achieved when the IBL filled the entire PBL. As indicated below, this is not an appropriate equilibrium value even for midlatitude lakes, and fails entirely in the tropics.

Results for lakes within forest (M=6.9; black curves), shrubland (M=4.6; blue curves), and grassland (M=2.3; red curves) environments are shown in Fig. 2. Figure 2a highlights the fact that IBL depth is not a function of M in the Jensen model (dashed curves), though it is obviously a strong function of A. Observed IBL depth data would help constrain the value of A and would be much easier to measure than surface stress – especially over open water for a sufficient range in fetch. Figure 2b shows that surface stress values are far from the equilibria proposed by Jensen (thin horizontal lines), but for all three surface types the surface stress approaches an asymptotic value at fetches around 5 km or less.

Figure 2 Theoretical IBL depths (a) and surface stress (b) as a function of fetch over water for a rough-to-smooth transition (offshore flow) for forested landscapes (black curves, M=6.9), shrubland (blue curves, M=4.6), and grassland (red curves, M=2.3). Solid curves – Panofsky and Townsend (1964) approach; dashed curves – Jensen (1978) approach with A=1 (thin dashed) and A=2 (thick dashed). Surface stress values are normalized by the upstream value. Thin horizontal lines represent equilibrium downstream values proposed by Jensen (1978) for downstream surface Rossby number Ro=10⁸.

The problem at hand is to estimate a mean surface wind stress value for small lakes, given that they are embedded within terrestrial environments of vastly different surface roughness and will thus be subject to an internal boundary layer that is almost certainly not in equilibrium with the lake surface. Figure 2 shows that the actual surface stress τ₀ is a function of fetch and always less than the equilibrium value τ₊. We have chosen the Jensen model with A=2 and a land cover of forest with M=6.9. A reasonable approach to estimate a lake-averaged surface wind stress is then the following.

1. Iteratively solve Eq. (4) for IBL depth d.

2. Solve Eq. (5) for the stress ratio ${\mathit{\tau }}_{\mathrm{0}}/{\mathit{\tau }}_{-}$.

3. Compute the average stress ratio over the mean fetch of the lake. The mean fetch can be taken as the square root of the surface area. This is exact for circular lakes but becomes progressively inaccurate as lakes become elongated and wind direction becomes important.

4. Divide this mean value by asymptotic value ${\mathit{\tau }}_{+}/{\mathit{\tau }}_{-}$. This is easily read from Fig. 2b or can be computed in the model by setting a threshold (e.g. 1 %) which the change in the stress ratio with fetch does not exceed. This yields the required stress reduction factor ${\mathit{\tau }}_{\mathrm{0}}/{\mathit{\tau }}_{+}$, i.e. the reduction in mean surface stress the lake experiences compared to the equilibrium surface layer estimate based on Monin–Obukhov similarity theory.

It is worth pointing out that Vickers and Mahrt (1997) found that surface drag decreased with fetch, because younger, growing waves are associated with larger surface drag than older waves (unless they are breaking). In reality, both processes may be occurring, but it is clear that for the smallest lakes the mechanism described here will dominate. Including the impact of wave state on the results reported here is certainly beyond the scope of the present study. Stepanenko et al. (2014) have proposed that excessive drag in the LAKE 1-D model for a simulation of Lake Valkea–Kotinen, a 4.1 ha boreal lake in southern Finland, can be compensated for by partitioning atmospheric momentum flux between wave development and surface currents. Our results provide an alternative (or perhaps additional) explanation. Fetch in this lake varies from about 100 to 400 m, and the lake is surrounded by forest. Examination of Fig. 2 suggests that surface stress (or equivalently the surface drag coefficient) is only a fraction of the equilibrium value: perhaps between 25 % and 50 %, depending on the model chosen.

We have applied this technique to a simulation of the CSLM for a small boreal lake (L239 at the Experimental Lakes Area, Canada) for 2013–2014 that has been described in MacKay et al. (2017). In that study the primary focus was on the winter season, though the simulation actually began in July. The first few days of this simulation (17–25 July 2013) are illustrated in Fig. 3, where observations are indicated with black curves, the original simulation with blue curves, and a modified simulation with red curves. Both 19 and 22 July were relatively windy days (Fig. 3a), with winds generally from the west or north-west (not shown) yielding a maximum fetch around 500 m. Based on the technique outlined above, we compute a drag reduction factor of 0.5 in the modified simulation. The impact has been an improved simulation of surface (i.e. 0.5 m depth) temperature (Fig. 3b) caused by a reduction in simulated depth of turbulent mixing (Fig. 3c), especially during or shortly following the wind events. A few observed and simulated temperature profiles during this period are shown in Fig. 3d, e. Observations suggest that between 19 and 20 July the epilimnion cooled from 24.7 to 23.4 ^∘C and deepened from 3.0 to 3.8 m. On 19 July both simulations and the observations show (Fig. 3d) identical epilimnion temperatures and depths (the simulations were initialized just 2 days earlier). However, by 20 July the original simulation (blue curves) indicates excessive deepening and cooling (4.5 m, 22.5 ^∘C) compared with the modified simulation (4.0 m, 23.2 ^∘C). Note that by 23 July (Fig. 3e), even though both simulations now show good agreement with the observed epilimnion temperature, an error of 0.5 m in the original simulated epilimnion depth persists.

Figure 3 Observed 2 m wind speed (a), observed and simulated 0.5 m water temperatures (b), and simulated diurnal turbulent mixed layer depths (c) for 17–25 July 2013. Temperature profiles (midnight) for (d) 19 July (dash) and 20 July (solid); and (e) 21 July (dash) and 23 July (solid). Standard simulation results (blue curves) taken from MacKay et al. (2017). Modified simulation results (red curves) have a surface drag coefficient reduced by 50 %. Observed curves are black.

With this approach we must keep in mind a number of caveats. Very large changes in roughness or very dense canopies are strictly speaking not consistently handled as no allowance has been made for changes in the zero plane displacement. Recent wind tunnel (Markfort et al., 2014), large-eddy simulation (Schlegel et al., 2015), and field (Detto et al., 2008) experiments all show complex vertical and horizontal structure in the Reynolds stress resulting from streamline displacement and the turbulent wake at a canopy edge. Markfort et al. (2014) showed that flow separation and a shallow recirculating zone can develop in the immediate lee of the forest edge, with flow reattachment some distance downstream – reminiscent of the well-known circulation pattern for a backward-facing step (e.g. Driver and Seegmiller, 1985). If one accounts for this reattachment distance (i.e. the origin in Fig. 2 is at the point of reattachment rather than the shoreline), then the empirical model for surface stress of Markfort et al. (2014) (their Fig. 16) appears similar to that proposed here. However, both studies consider only neutrally stratified boundary layer conditions. Unstable conditions such as would occur in autumn and early winter before ice-on would almost certainly have an impact on the turbulent canopy wake and the distribution of Reynolds stress. This is currently an active area of research. In addition, unlike for rigid surfaces, steady wind over open water will generally lead to an increase in surface roughness with time. Note that computing the surface drag coefficient based on lake surface roughness values of 10⁻² and 10⁻⁴ m changes the stress reduction factors to 0.59 and 0.45 respectively, similar to the factor of 0.5 proposed here. All of these issues could likely be incorporated into the theory at the expense of increasing complexity; nevertheless, this approach should at least qualitatively describe the importance of IBLs and fetch-varying wind stress over small lakes and appears a suitable first step.

3 Sediment heat flux

Currently the CSLM employs an adiabatic boundary condition at the lake bottom. In the recent simulation of L239 mentioned above, it was found that during a 100 d ice-covered period (2013–2014) the simulated change in lake heat content corresponded to a mean bias of about −2.5 W m² compared to observations (MacKay et al., 2017). The authors noted that this was of the same order of magnitude as the wintertime sediment heat flux found in a number of previous studies and suggested this warranted further investigation. Here we test whether the inclusion of a simple sediment heat flux scheme can ameliorate this bias.

A simple scheme for the storage and subsequent flux of heat from lake sediments can be constructed by considering a sediment slab of fixed thickness, and uniform thermal conductivity and volumetric heat capacity. Such a scheme is simpler than some existing multilayer sediment models (Stepanenko et al., 2016), but has the appeal of not requiring many additional levels to be carried in global weather and climate simulations. Mean temperature in the slab evolves due to thermal and radiative fluxes at the lake water–sediment interface in such a way as to conserve energy. The boundary condition at the base of the slab is isothermal, rather than adiabatic as is sometimes assumed (e.g. Stepanenko et al., 2016), which places a constraint on the minimum slab thickness. Note that for terrestrial (i.e. soil or rock) surfaces the diurnal temperature wave is believed to penetrate about 1 m and the annual temperature wave penetrates about 20 m (e.g. Carslaw and Jaeger, 1959), below which geothermal heating acts to maintain a constant temperature gradient. A number of studies have found that temperature oscillations in lake sediments are substantially damped after only a few meters, below which temperatures are essentially isothermal (e.g. Likens and Johnson, 1969; Tsay et al., 1992). Here we consider a slab thickness of 10 m, with thermal properties consistent with sand. A layer of pure sand has a thermal conductivity of about 2.5 W K⁻¹ m⁻¹ (about 4 times larger than that for liquid water) and a volumetric heat capacity of about 2.13×10⁶ J m⁻³ K⁻¹ (about half that of liquid water), and these values are adopted here. Values for clay would be similar, though sediments with significant organic matter would differ. The system is solved by asserting the continuity of both temperature and heat flux at the water–sediment interface (i.e. option 1 in Stepanenko et al., 2016) and a lower boundary condition temperature of 6.0 ^∘C.

The difficulty with applying such an approach here is that the CSLM does not include lake bathymetry information, and lake bottom is everywhere assumed to be at the mean lake depth. In fact, many lakes have extensive shallow littoral zones in which we would expect much larger shortwave (SW) insolation, and thus sediment heat flux, than in more “bathtub”-shaped lakes which nevertheless have the same mean depth. Here we propose a scheme to compute the net sediment SW insolation based on minimal bathymetric data – maximum and mean depth only.

In what follows, all lower-case variables are dimensionless – scaled by appropriate length and depth scales for the lake in question. Figure 4a shows a family of one-parameter theoretical “hypsographs” given by

$$y=x^s,$$

where x is the (dimensionless) radial distance from the lake center, y is the (dimensionless) height from the lake bottom at maximum depth, and s is a shape factor. Theoretical, axially symmetric lakes are formed by rotating these curves around the line x=0. Thus for s=1 the lake is conical, for values of s less than 1 the lake takes on a “birdbath” shape, while for values of s greater than 1 it is more “bathtub” shaped. In the limit of very large s the lake becomes a right circular cylinder (the default shape assumed in CSLM). Birdbath lakes have extensive shallow littoral zones whose sediment would absorb much more SW radiation than bathtub lakes of the same surface area.

Figure 4 (a) Theoretical “hypsographs” for idealized axially symmetric lakes with different shape factors s indicated; (b) schematic to aid in determining the volume of an axially symmetric lake; (c) normalized observed hypsograph (blue) and theoretical estimate (red) for L239. See text for details.

Reference to Fig. 4b shows that the volume of revolution of a thin disk at height y is

$$\mathrm{d}v=\mathit{\pi }{x}^{\mathrm{2}}\mathrm{d}y,$$

so that the total volume is simply

$$v=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathit{\pi }{x}^{\mathrm{2}}\mathrm{d}y=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathit{\pi }{y}^{\mathrm{2}/s}\mathrm{d}y=\mathit{\pi }{\displaystyle \frac{s}{s+\mathrm{2}}}.$$

The dimensional volume is thus

$$V=\mathit{\pi }{\displaystyle \frac{s}{s+\mathrm{2}}}{L}^{\mathrm{2}}{H}_{\max },$$

where H_max is the maximum depth and L is the radius of a circular lake with the same surface area as the lake in question. We can determine the value for s by equating this volume to the actual volume of the lake, or equivalently,

$$\mathit{\pi }{\displaystyle \frac{s}{s+\mathrm{2}}}{L}^{\mathrm{2}}{H}_{\max }=\mathit{\pi }{L}^{\mathrm{2}}\stackrel{\mathrm{‾}}{H},$$

where $\stackrel{\mathrm{‾}}{H}$ is the mean lake depth. This yields

$$s={\displaystyle \frac{\mathrm{2}\frac{\stackrel{\mathrm{‾}}{H}}{H_{\max }}}{\mathrm{1}-\frac{\stackrel{\mathrm{‾}}{H}}{H_{\max }}}}.$$

For L239 we have $\stackrel{\mathrm{‾}}{H}=\mathrm{11.0}$ m and H_max=30.0 m, so that s=1.16. Figure 4c compares the actual (normalized) hypsograph for L239 compared with our theoretical curve.

Figure 5 Schematic to aid in the estimation of sediment shortwave insolation for an axially symmetric lake with a given shape factor s<1.

The task now is to estimate the mean absorbed SW radiation at the sediment–water interface given this lake shape profile. It is clear from Fig. 5 that if i₀ is the intensity of SW radiation reaching the lake surface (per unit area), then the intensity reaching the sediment at radial distance x is given by

$$\begin{array}{}\text{(7)}& i_{\mathrm{sed}}\left(x\right)=i_{\mathrm{0}}\cos \mathit{\theta }{e}^{-\mathit{\beta }\left(\mathrm{1}-x^s\right)},\end{array}$$

where β is the (dimensionless) shortwave extinction (i.e. $\mathit{\beta }=\widehat{\mathit{\beta }}{H}_{\max }$ where $\widehat{\mathit{\beta }}$ is the dimensional extinction) and 1-x^s is the depth at x. Consider the surface area of a thin ring of radius x and width dl. Clearly

$$\mathrm{d}a=\mathrm{2}\mathit{\pi }x\mathrm{d}l,$$

so that the total radiation reaching this elemental surface is

$$\begin{array}{r}{i}_{\mathrm{sed}}\left(x\right)\mathrm{d}a=i_{\mathrm{0}}\cos \mathit{\theta }{e}^{-\mathit{\beta }\left(\mathrm{1}-x^s\right)}\mathrm{2}\mathit{\pi }x\mathrm{d}l\\ =\mathrm{2}\mathit{\pi }{i}_{\mathrm{0}}\mathrm{x}{e}^{-\mathit{\beta }\left(\mathrm{1}-x^s\right)}\mathrm{d}x.\end{array}$$

The total insolation of the lake sediment surface is found by integrating over all x, and the mean sediment insolation per unit lake surface area is found by dividing this integral by the lake surface area. Recalling that the dimensionless lake surface area is simply π, we get

$$\begin{array}{}\text{(8)}& \begin{array}{r}{\stackrel{\mathrm{‾}}{i}}_{\mathrm{sed}}={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\mathrm{2}\mathit{\pi }{i}_{\mathrm{0}}x{e}^{-\mathit{\beta }\left(\mathrm{1}-x^s\right)}\mathrm{d}x\\ =i_{\mathrm{0}}\left[\mathrm{2}{e}^{-\mathit{\beta }}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}x{e}^{\mathit{\beta }{x}^s}\mathrm{d}x\right].\end{array}\end{array}$$

This has analytic solutions in terms of the gamma and incomplete gamma functions:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{i}}_{\mathrm{sed}}& =\mathrm{2}{i}_{\mathrm{0}}{e}^{-\mathit{\beta }}\left\{{\displaystyle \frac{\mathrm{1}}{s}}{\left(-\mathit{\beta }\right)}^{-\frac{\mathrm{2}}{s}}\left[\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{2}}{s}}\right)-\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{2}}{s}},-\mathit{\beta }\right)\right]\right\};\\ s& >\mathrm{0}.\end{array}\end{array}$$

Note that for s=0 (i.e. a birdbath of zero depth everywhere except at the center) the solution Eq. (9) does not apply, but Eq. (8) can be integrated trivially to get

$${\stackrel{\mathrm{‾}}{i}}_{\mathrm{sed}}\left(s=\mathrm{0}\right)=\mathrm{2}{i}_{\mathrm{0}}{e}^{-\mathit{\beta }}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}x{e}^{\mathit{\beta }}\mathrm{d}x=i_{\mathrm{0}},$$

which is the correct limit. On the other hand, for very large s our lake becomes cylindrical and Eq. (8) can also be solved trivially. Expanding the exponential in a Taylor series, we find

$$\begin{array}{rl}{\stackrel{\mathrm{‾}}{i}}_{\mathrm{sed},\mathrm{\infty }}=& \stackrel{\lim }{s\to \mathrm{\infty }}\left\{\mathrm{2}{i}_{\mathrm{0}}{e}^{-\mathit{\beta }}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}x{e}^{\mathit{\beta }{x}^s}\mathrm{d}x\right\}\\ & =\mathrm{2}{i}_{\mathrm{0}}{e}^{-\mathit{\beta }}\stackrel{\lim }{s\to \mathrm{\infty }}\left\{\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}x\left(\mathrm{1}+\mathit{\beta }{x}^s+{\displaystyle \frac{(\mathit{\beta }{x}^s{)}^{\mathrm{2}}}{\mathrm{2}\mathrm{!}}}+\mathrm{\dots }\right)\mathrm{d}x\right\}\\ & =\mathrm{2}{i}_{\mathrm{0}}{e}^{-\mathit{\beta }}\stackrel{\lim }{s\to \mathrm{\infty }}\left\{\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{x}^{\mathrm{2}}+{\displaystyle \frac{\mathit{\beta }}{s+\mathrm{2}}}{x}^{s+\mathrm{2}}+{\displaystyle \frac{{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{2}(\mathrm{2}s+\mathrm{2})}}{x}^{\mathrm{2}s+\mathrm{2}}\right.\right.\\ & \left.{\left.+\mathrm{\dots }\right]}_{\mathrm{0}}^{\mathrm{1}}\right\}=i_{\mathrm{0}}{e}^{-\mathit{\beta }},\end{array}$$

which is again the correct limit.

For L239 we have found s=1.16 and the (dimensionless) extinction is β=27, so that Eq. (8) or Eq. (9) becomes

$${\stackrel{\mathrm{‾}}{i}}_{\mathrm{sed}}=\mathrm{0.06}{i}_{\mathrm{0}}.$$

In other words, of the net surface SW radiation on L239, about 6 % reaches the sediments, leaving 94 % to be absorbed by the lake water on a lake-wide average. In the standard simulation assuming the entire lake is at the mean depth (11 m), the sediment insolation is only

$${\stackrel{\mathrm{‾}}{i}}_{\mathrm{sed}}=\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{5}}{i}_{\mathrm{0}},$$

3 orders of magnitude less. The approach taken here is to simply reduce the incoming net SW insolation at the lake surface by 6 % and apply this energy flux directly at the lake water–sediment interface. In this way we estimate the mean sediment insolation for any lake whose surface insolation is known, along with extinction, mean and maximum depths. For example, Table 1 lists sediment SW insolation fractions for a variety of shape factors s for lakes with the same (dimensionless) extinction as L239.

Table 1 Sediment shortwave insolation fractions for lakes with dimensionless shortwave extinction β=27 and various shape factors s.

The above approach was used in two new simulations of L239 for 2013–2014 (Fig. 6). In the standard simulation for this period (MacKay et al., 2017) the lake bottom was adiabatic. In order to isolate the impact of geothermal heating alone, we relaxed the adiabatic boundary condition as described above, but set the sediment SW insolation fraction to 0.0 (experiment X1; Fig. 6 blue curves). The second experiment (experiment X2) is identical to the first, except we set the sediment SW insolation fraction to the proposed value of 0.06 (Fig. 6, red curves). The simulated ice thicknesses (Fig. 6a) are virtually identical in the two experiments, though both the sediment temperatures (dotted curves) and lowest water layer temperatures (solid curves) differ even well before ice-on (Fig. 6b). Note that the sediment heat flux is close to zero or negative (i.e. into the water column) throughout these simulations, with values never exceeding a few W m⁻² (Fig. 6b, dashed curves, right-hand scale). During the ice-cover season the geothermal component is significant, generally about half of the total. Finally, the impact of the new sediment heat flux scheme on lake thermal structure during the ice-cover period is illustrated in Fig. 7. Ice phenology and thickness for the original and modified simulations are virtually identical (Fig. 7a). Temperature profiles in the original experiment show no sign of water column warming (Fig. 7b) – in fact, a general cooling trend is evident in the top 5 m until ablation begins (after snow cover is gone – not shown). The observations on the other hand show clear signs of warming throughout most of the column under ice (Fig. 7c). Results from the sediment heat flux experiment (X2) more closely resemble this pattern (Fig. 7d), though there are obvious deficiencies in the simulation, including a more stratified region for several meters below the ice as well as a general warm bias in deeper waters and a too strongly stratified surface layer near and following ice-off.

Figure 6 (a) Simulated ice thickness; (b) simulated sediment temperature (dotted), lowest-level water temperature (solid), and sediment heat flux (dashed – right-hand scale). Simulated results are from experiments X1 (blue) and X2 (red). Zero sediment heat flux (dashed black curve) is also indicated in (b). The period shown covers 19 July 2013–25 May 2014.

Figure 7 Simulated ice thickness (a), and simulated (b) (d) and observed (c) temperature profiles for 1 November 2013–20 May 2014. The standard simulation (MacKay et al., 2017) is shown in (b) and the blue curve in (a). Simulation X2 is shown in (d) and the red curve in (a).

The change in total heat content under ice is, however, improved. MacKay et al. (2017) found that between 26 November 2013 and 6 March 2014 (100 d) the simulated 1–10 m water column lost 7.26×10⁶ (J m⁻²) corresponding to a mean heat flux of −0.84 (W m⁻²), whereas observations suggested warming corresponding to +1.66 (W m⁻²). Including sediment heat flux improves this. Geothermal heating alone (X1) brings the mean heat flux into the column up to 0.08 (W m⁻²), while including the radiative forcing (X2) yields 1.05 (W m⁻²) for a net bias of −0.61 (W m⁻²). These results are summarized in Table 2.

Table 2 Observed and simulated change in 1–10 m heat content from 26 November 2013 to 6 March 2014.

A disadvantage of the method described here is that all sediment heat flux into the lake water takes place at the mean lake depth (which is the lowest model level). Circulation processes that redistribute heat are clearly active in the observations (Fig. 7c), but are not represented in the model. In fact, the only reason the simulated heat is redistributed vertically is that the bottom simulated temperatures are near the temperature of maximum density (4 ^∘C), so that warming of these waters produced convection. But temperatures in the bottom half of the simulated lake are clearly biased warm during ice cover, and free convection is likely not taking place in reality. So while the change in total heat content under ice has improved with the new scheme, work remains for parameterizing mixing processes to appropriately redistribute heat emanating from the sediments.

4 Conclusions

Lakes are increasingly being recognized as important components of the land surface, and 1-D modelling schemes are currently being developed for inclusion in a number of climate and numerical weather prediction systems around the world. Lakes interact with their environment through flux exchange with their bottom sediments and at the surface with the atmosphere, and even a perfect lake model will perform poorly in a prognostic sense if these linkages are represented poorly. However, these are challenging areas of research: high-quality data in lake sediments (e.g. temperature, thermal properties) are difficult to achieve at regional or global scales, and the process understanding itself of turbulent exchange under non-equilibrium conditions, ubiquitous for all but the largest of boreal lakes, is lacking. This study contributes to this discussion with two simple schemes that begin to represent these linkages in simple 1-D “bathtub”-like models like the CSLM.

Terrestrial landscapes play a role in lake hydrodynamics through the sudden drop in aerodynamic roughness that the atmosphere encounters at the shoreline of virtually any boreal lake. While some models do make adjustments to surface drag to account for various processes (or for “tuning”), none to our knowledge factors in the terrestrial roughness in a quantitative way. Here we propose a straightforward scheme based on earlier IBL research that begins to address this. When applied to a simulation of the CSLM we have found an improvement in the near-surface temperature following two wind events a few days apart that resulted from a more realistically simulated diurnal mixed layer depth.

It is well known that shallow lakes or lakes with extensive shallow littoral zones may be subject to sediment heating that is significant for many applications, and several models have attempted to account for this process. While it is trivial to add a sediment layer beneath a lake in order to absorb and subsequently release heat, problems arise when the lake model itself does not represent bathymetry. Since a significant fraction of the sediment heat content can arise from SW insolation at the sediment–water interface, lake models that assume uniform depth (i.e. constant depth set equal to the lake mean depth) would in most cases receive almost no radiative forcing at the sediment unless the lake was exceptionally shallow or clear. To account for this a scheme is proposed to estimate the actual sediment insolation by approximating any given lake as an axially symmetric lake with a shape factor determined from mean and maximum depth data only. When applied in a simulation of the CSLM we have found an improvement in the change in lake heat content under ice cover over a 100 d period. Issues remain regarding the distribution of the sediment heating throughout the water column, as currently all heating takes place at the lake mean depth.