2.1 Energetics of a radiatively forced ice-covered waterbody

Mass and momentum budgets in (2.1) allow a simple derivation of the mechanical energy budget, ${\dot{E}}_{M}$ , that involves the rate of change of $E_{k}$ and  $E_{p}$ . Here, we examine the energy quantities per unit mass, $\unicode[STIX]{x1D70C}_{0}V$ , adopting the standard $\unicode[STIX]{x1D6F7}$ -notation to denote energetics in ${\dot{E}}_{M}$ (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995, among others). To derive the kinetic energy budget, ${\dot{E}}_{k}$ , we multiply (2.1c ) by the velocity field $\boldsymbol{v}$ . Then, using continuity, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0$ , and the boundary conditions on  $S$ , the rate of change of $E_{k}\equiv (1/2)\langle \boldsymbol{v}^{2}\rangle _{V}$ in $V$ is reduced to

(2.2) $$\begin{eqnarray}{\dot{E}}_{k}=-\underbrace{\frac{g}{\unicode[STIX]{x1D70C}_{0}}\langle \unicode[STIX]{x1D70C}w\rangle _{V}}_{\unicode[STIX]{x1D6F7}_{z}}-\underbrace{\unicode[STIX]{x1D708}\langle |\unicode[STIX]{x1D735}u|^{2}+|\unicode[STIX]{x1D735}v|^{2}+|\unicode[STIX]{x1D735}w|^{2}\rangle _{V}}_{\unicode[STIX]{x1D716}},\end{eqnarray}$$

where the rate of work across $S$ is null. $\unicode[STIX]{x1D6F7}_{z}$ is the volume-averaged, $\langle \cdot \rangle _{V}$ , reversible buoyancy flux, and $\unicode[STIX]{x1D716}$ is the volume-averaged kinetic energy dissipation rate. Similarly, we multiply the mass transport equation (2.1a,b ) by ‘ $gz$ ’ and integrate over $V$ to derive the rate of change of potential energy $E_{p}\equiv (g/(\unicode[STIX]{x1D70C}_{0}V_{S}))\langle z\unicode[STIX]{x1D70C}\rangle _{V}$ in  $V$ ,

(2.3) $$\begin{eqnarray}{\dot{E}}_{p}=\unicode[STIX]{x1D6F7}_{z}+\underbrace{\frac{\unicode[STIX]{x1D705}g}{\unicode[STIX]{x1D70C}_{0}V}\oint _{S}z\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\,\text{d}S}_{\unicode[STIX]{x1D6F7}_{b}}+\underbrace{-\frac{\unicode[STIX]{x1D705}g}{\unicode[STIX]{x1D70C}_{0}}\left\langle \frac{\text{d}\unicode[STIX]{x1D70C}}{\text{d}z}\right\rangle _{V}}_{\unicode[STIX]{x1D6F7}_{i}}+\underbrace{\frac{g}{\unicode[STIX]{x1D70C}_{0}}\langle z\unicode[STIX]{x1D735}F_{r}\rangle _{V}}_{\unicode[STIX]{x1D6F7}_{r}},\end{eqnarray}$$

where $\hat{\boldsymbol{n}}$ is the outward normal vector to the surface  $S$ , $\unicode[STIX]{x1D6F7}_{b}$ is the rate of diffusive buoyancy flux across the surface  $S$ , $\unicode[STIX]{x1D6F7}_{i}$ is the rate at which $E_{p}$ increases due to vertical mass redistribution resulting from molecular diffusion (conversion rate from internal energy, $E_{i}$ , to  $E_{p}$ ), and $\unicode[STIX]{x1D6F7}_{r}$ is the conversion rate from the radiative heating to  $E_{p}$ . Understanding the role and contribution of $\unicode[STIX]{x1D6F7}_{r}$ in ${\dot{E}}_{M}$ is the main focus of § 2. The potential energy can be decomposed as a background and an available component, $E_{p}=E_{b}+E_{a}$ . The background potential energy, $E_{b}\equiv (g/(\unicode[STIX]{x1D70C}_{0}V_{S}))\langle z_{\ast }\unicode[STIX]{x1D70C}\rangle _{V}$ , is defined in terms of the global reference position $z_{\ast }(\boldsymbol{x},t)$ (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995) that characterises the height to which each fluid parcel is allocated to reach its stable position in the volume  $V$ . The evolution equation of $E_{b}$ in $V$ is given by

(2.4) $$\begin{eqnarray}{\dot{E}}_{b}=\underbrace{\frac{\unicode[STIX]{x1D705}g}{\unicode[STIX]{x1D70C}_{0}V_{S}}\oint _{S}z_{\ast }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\,\text{d}S}_{\unicode[STIX]{x1D6F7}_{b_{\ast }}}+\underbrace{-\frac{\unicode[STIX]{x1D705}g}{\unicode[STIX]{x1D70C}_{0}}\left\langle \frac{\text{d}z_{\ast }}{\text{d}\unicode[STIX]{x1D70C}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}|^{2}\right\rangle _{V}}_{\unicode[STIX]{x1D6F7}_{d}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6F7}_{b_{\ast }}$ is the rate of change of $E_{b}$ due to diffusive buoyancy fluxes across $A_{ice}$ , and $\unicode[STIX]{x1D6F7}_{d}$ is the rate of change of $E_{b}$ resulting from irreversible mixing of the density field in  $V$ . The evolution equation for $E_{a}$ results from subtracting (2.4) from (2.3),

(2.5) $$\begin{eqnarray}{\dot{E}}_{a}=\unicode[STIX]{x1D6F7}_{r}+\unicode[STIX]{x1D6F7}_{z}+(\unicode[STIX]{x1D6F7}_{b}-\unicode[STIX]{x1D6F7}_{b_{\ast }})-(\unicode[STIX]{x1D6F7}_{d}-\unicode[STIX]{x1D6F7}_{i}).\end{eqnarray}$$

We stress that the surface buoyancy flux $\unicode[STIX]{x1D6F7}_{b}$ is defined null across $S_{b}$ but might differ from zero across the ice-surface boundary $A_{ice}$ , at $z=z_{ice}$ . Thus

(2.6) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{b}=\frac{\unicode[STIX]{x1D705}g}{\unicode[STIX]{x1D70C}_{0}V}\oint _{A_{ice}}z\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}\,\text{d}S=\frac{\unicode[STIX]{x1D705}gA_{ice}z_{ice}}{\unicode[STIX]{x1D70C}_{0}V}\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}\right\rangle _{A_{ice}}.\end{eqnarray}$$

Similarly, the surface buoyancy flux associated with the evolution of $E_{b}$ at the level $z_{\ast }$ is determined by

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{b_{\ast }}=\frac{\unicode[STIX]{x1D705}g}{\unicode[STIX]{x1D70C}_{0}V}\oint _{A_{ice}}z_{\ast }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}\,\text{d}S=\frac{\unicode[STIX]{x1D705}gA_{ice}z_{\ast }(z_{ice})}{\unicode[STIX]{x1D70C}_{0}V}\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}\right\rangle _{A_{ice}}.\end{eqnarray}$$

Thus, we can express the net surface buoyancy flux, $(\unicode[STIX]{x1D6F7}_{b}-\unicode[STIX]{x1D6F7}_{b_{\ast }})$ , across the ice boundary working on the rate of change of $E_{a}$ , as

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{b}-\unicode[STIX]{x1D6F7}_{b_{\ast }}=[z_{ice}-z_{\ast }(z_{ice})]\frac{\unicode[STIX]{x1D705}gA_{ice}}{\unicode[STIX]{x1D70C}_{0}V}\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}\right\rangle _{A_{ice}}.\end{eqnarray}$$

However, since fluid parcels in the vicinity of the ice–water interface become lighter and strongly stable as $z\rightarrow z_{ice}$ , the reference height of a fluid parcel of density $\unicode[STIX]{x1D70C}(z\approx z_{ice})$ is $z_{\ast }(z\approx z_{ice})=z_{ice}$ . This fundamental thermodynamic characteristic that the ice–water boundary imposes to the fluid stability leads to the following equality:

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{b_{\ast }}|_{A_{ice}}=\unicode[STIX]{x1D6F7}_{b}|_{A_{ice}}.\end{eqnarray}$$

Hence, by (2.9), the net surface buoyancy flux (2.8) is null and therefore plays no role in ${\dot{E}}_{a}$ (2.5). Notice that $(\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}z)|_{z=z_{ice}}\leqslant 0$ , so $\unicode[STIX]{x1D6F7}_{b_{\ast }}|_{z=z_{ice}}=\unicode[STIX]{x1D6F7}_{b}|_{z=z_{ice}}\leqslant 0$ is a negative flux that decreases both $E_{p}$ and $E_{b}$ equally.

2.2 The pathway of $\unicode[STIX]{x1D6F7}_{r}$ in the mechanical energy budget

The evolution equation for the mechanical energy results by merging (2.2) and (2.3),

(2.10) $$\begin{eqnarray}{\dot{E}}_{M}={\dot{E}}_{k}+{\dot{E}}_{p}=-\unicode[STIX]{x1D716}+\unicode[STIX]{x1D6F7}_{b}+\unicode[STIX]{x1D6F7}_{i}+\unicode[STIX]{x1D6F7}_{r}.\end{eqnarray}$$

Replacing $\unicode[STIX]{x1D6F7}_{b}=\unicode[STIX]{x1D6F7}_{b_{\ast }}={\dot{E}}_{b}-\unicode[STIX]{x1D6F7}_{d}$ into (2.10), we obtain the evolution equation for the available mechanical energy,  $E_{AM}$ ,

(2.11) $$\begin{eqnarray}{\dot{E}}_{AM}\equiv {\dot{E}}_{k}+{\dot{E}}_{a}={\dot{E}}_{M}-{\dot{E}}_{b}=-\unicode[STIX]{x1D716}-(\unicode[STIX]{x1D6F7}_{d}-\unicode[STIX]{x1D6F7}_{i})+\unicode[STIX]{x1D6F7}_{r}.\end{eqnarray}$$

The energy budget in (2.11) retains the elemental energetics governing the loss and gain of available mechanical energy in the system. From (2.11), we distinguish that $\unicode[STIX]{x1D716}$ represents the net decay rate of  $E_{k}$ , while $(\unicode[STIX]{x1D6F7}_{d}-\unicode[STIX]{x1D6F7}_{i})-\unicode[STIX]{x1D6F7}_{r}$ represents the net rate of change (either decay or gain) of $E_{a}$ in the system. Thus, the evolution of $E_{a}$ results from the competition between its decay rate, controlled by the action of irreversible mixing, $\unicode[STIX]{x1D6F7}_{d}-\unicode[STIX]{x1D6F7}_{i}$ , and its increment, resulting from the conversion of radiative to potential energy,  $\unicode[STIX]{x1D6F7}_{r}$ . If $\unicode[STIX]{x1D6F7}_{r}$ becomes null, then $E_{AM}$ will decay in time. Therefore, it is inferred that $\unicode[STIX]{x1D6F7}_{r}$ is the fundamental generator of $E_{a}$ and the implicit driver for  $E_{k}$ .

The pathway of radiative energy flux, $\unicode[STIX]{x1D6F7}_{r}$ , into the mechanical energetics in an ice-covered enclosed waterbody is illustrated in figure 1(b). The energy pathway is described by six main processes. (1) During the daytime, $\unicode[STIX]{x1D6F7}_{r}$ works over upper waters driving localised increments of temperature and density leading to unstable density distributions and a gain of  $E_{a}$ . (2) Part of $E_{a}$ is transformed to buoyancy flux $\unicode[STIX]{x1D6F7}_{z}$ via mechanical convection that tries to restore the system to a stable density distribution. In this process, $\unicode[STIX]{x1D6F7}_{z}$ is transformed into  $E_{k}$ , energising the vertical and horizontal transport. (3) Simultaneously, diapycnal mixing $\unicode[STIX]{x1D6F7}_{d}$ transfers $E_{a}$ towards  $E_{b}$ , and this process is enhanced once convection starts stirring the well-mixed layer and deep stratified regions. (4) Convection also enhances the rate of decay of kinetic energy,  $\unicode[STIX]{x1D716}$ , by friction which is transformed into internal energy. (5)  $E_{p}$ increases at the expense of the internal energy contained in the stratified fluid at a rate of  $\unicode[STIX]{x1D6F7}_{i}$ . (6) Finally, the system loses $E_{b}$ via the surface buoyancy flux at a rate of $\unicode[STIX]{x1D6F7}_{b_{\ast }}$ controlled by the density gradient in the TBL. Thus, the energy pathway in our system can be understood as a spiral cycle characterised by the gain of $E_{AM}$ (which, as a consequence, enhances the rate of increment of  $E_{b}$ ) as long as the radiative buoyancy flux, $\unicode[STIX]{x1D6F7}_{r}$ , works on the fluid. At night, $\unicode[STIX]{x1D6F7}_{r}=0$ and $E_{AM}$ starts decaying until the beginning of the next daytime.

We do not examine the internal energy reservoir nor its link with ${\dot{E}}_{M}$ . For an extensive analysis of the internal energy decomposition, we refer to Tailleux (Reference Tailleux2009). Additionally, although we do not develop a radiative energy balance, it is assumed that $\unicode[STIX]{x1D6F7}_{r}$ results from the net component of the radiative energy that reaches the volume  $V$ . The radiative balance (including downward/upward fluxes) on ice-covered lakes is addressed by Kirillin et al. (Reference Kirillin, Leppäranta, Terzhevik, Granin, Bernhardt, Engelhardt, Efremova, Golosov, Palshin, Sherstyankin, Zdorovennova and Zdorovennov2012) and Bouffard et al. (Reference Bouffard, Zdorovennov, Zdorovennova, Pasche, Wüest and Terzhevik2016).

3.1 Cumulative mixing efficiency

We define the ‘cumulative mixing efficiency’ (e.g., Manins & Turner Reference Manins and Turner1978; Peltier & Caulfield Reference Peltier and Caulfield2003; Davies Wykes & Dalziel Reference Davies Wykes and Dalziel2014) as

(3.5) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{c}\equiv \frac{\displaystyle \int _{t}^{t+\unicode[STIX]{x1D70F}}\dot{E_{b}}\,\text{d}t^{\prime }}{\displaystyle \int _{t}^{t+\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D6F7}_{r}(t^{\prime })\,\text{d}t^{\prime }}\equiv \frac{\unicode[STIX]{x0394}E_{b,\unicode[STIX]{x1D70F}}}{E_{r,\unicode[STIX]{x1D70F}}},\end{eqnarray}$$

where $\unicode[STIX]{x0394}E_{b,\unicode[STIX]{x1D70F}}=E_{b}(t+\unicode[STIX]{x1D70F})-E_{b}(t)$ corresponds to the change of $E_{b}$ over a time scale $\unicode[STIX]{x1D70F}$ , while $E_{r,\unicode[STIX]{x1D70F}}$ is the radiative energy injected into the bulk via radiative buoyancy flux, $\unicode[STIX]{x1D6F7}_{r}$ , during  $\unicode[STIX]{x1D70F}$ . Notice that (3.5) will not be equal to the average of (3.4) over the time scale $\unicode[STIX]{x1D70F}$ unless the system is in steady state.

The ratio $\unicode[STIX]{x1D702}_{c}$ is straightforward to compute if the whole density field and radiative flux are known. Numerical experiments provide a robust approach for examining both $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D702}_{c}$ in ice-covered flows forced by $\unicode[STIX]{x1D6F7}_{r}$ . However, considering the standard methodologies used to obtain in situ observations, the ‘cumulative’ formulation fits better to estimate the rate of mixing in natural systems.

3.2 Cumulative mixing efficiency based on vertically integrated observations

We formulate a simple, handy methodology to estimate $\unicode[STIX]{x1D702}_{c}$ based on in situ observations of density $\unicode[STIX]{x1D70C}$ (obtained from casts of temperature, $T$ , and conductivity, $C$ , as a function of depth), solar irradiance at the ice–water interface, $H_{0}$ , and light distribution across the water column in a basin of volume  $V$ . For simplicity, we consider that $H_{0}$ is uniformly distributed over the ice–water interface of area $A_{ice}$ and that the hypsometric curve is known, ${\mathcal{A}}(z)$ , which characterises the horizontal area of the waterbody as a function of depth  $z$ . Notice that ${\mathcal{A}}(z=z_{ice})=A_{ice}$ .

Let us assume we measure $N_{p}$ vertical density profiles, $\unicode[STIX]{x1D70C}^{j}$ ( $j=1,\ldots ,N_{p}$ ), throughout  $V$ . Furthermore, consider that each profile spans from $z_{ice}$ to $z_{bot}^{(j)}$ , where $z_{bot}^{(j)}$ might be different for each profile  $j$ . The cumulative mixing efficiency can be estimated using the following steps.

1. (i) Sort all the density profiles $\unicode[STIX]{x1D70C}^{j}$ to obtain a unique stable density distribution $\unicode[STIX]{x1D70C}_{\ast }$ increasing monotonically with depth from the top $z_{ice}$ to the deepest bottom  $D$ .

2. (ii) Interpolate $\unicode[STIX]{x1D70C}_{\ast }$ into an equally spaced reference state $z_{\ast }$ that allocates fluid parcels at a stable position. This procedure establishes a one-to-one function between $z_{\ast }$ and  $\unicode[STIX]{x1D70C}_{\ast }$ .

3. (iii) Estimate $E_{b}$ ( $\text{J}~\text{kg}^{-1}$ ) based on $\unicode[STIX]{x1D70C}_{\ast }(z_{\ast },t)$ and the hypsometric curve ${\mathcal{A}}(z_{\ast })$ ,

   (3.6) $$\begin{eqnarray}\langle E_{b}\rangle _{N_{p}}\approx \frac{g}{\unicode[STIX]{x1D70C}_{0}V}\int _{D}^{z_{ice}}z_{\ast }(t)\unicode[STIX]{x1D70C}_{\ast }(z_{\ast },t){\mathcal{A}}(z_{\ast })\,\text{d}z_{\ast }.\end{eqnarray}$$

4. (iv) Estimate the solar irradiance flux at the ice–water interface, $H_{0}(t)$ ( $\text{W}~\text{m}^{-2}$ ), and the light distribution over the water column. From the light distribution, estimate a bulk e-folding decay depth scale,  $\unicode[STIX]{x1D706}$ . Then, the radiative-induced density flux, $F_{r}(z,t)$ ( $\text{kg}~\text{m}^{-2}~\text{s}^{-1}$ ), can be modelled as

   (3.7) $$\begin{eqnarray}F_{r}(z,t)\approx F_{0}(z,t)\text{e}^{z/\unicode[STIX]{x1D706}},\end{eqnarray}$$
   where $F_{0}=(\unicode[STIX]{x1D6FC}_{T}/c_{p})H_{0}(t)$ , while $\unicode[STIX]{x1D6FC}_{T}$ and $c_{p}$ are the thermal expansion coefficient and the heat capacity of water, respectively. Note that the procedure is not restricted to an exponential decrease of light but rather applies to any vertical light absorption distribution.

5. (v) Estimate the radiative buoyancy flux ( $\text{W}~\text{kg}^{-1}$ ):

   (3.8) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{r}(t)\approx \frac{g}{\unicode[STIX]{x1D70C}_{0}}\frac{1}{D}\int _{D}^{z_{ice}}z\frac{\text{d}F_{r}}{\text{d}z}\,\text{d}z.\end{eqnarray}$$

6. (vi) Estimate the time-integrated radiative energy over a period $\unicode[STIX]{x1D70F}$ , $E_{r,\unicode[STIX]{x1D70F}}$ ( $\text{J}~\text{kg}^{-1}$ ) defined in (3.5).

7. (vii) The cumulative mixing efficiency over a period $\unicode[STIX]{x1D70F}$ is estimated using (iii) and (vi):

   (3.9) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{c}\approx \frac{\langle E_{b}(t+\unicode[STIX]{x1D70F})\rangle _{N_{p}}-\langle E_{b}(t)\rangle _{N_{p}}}{E_{r,\unicode[STIX]{x1D70F}}}.\end{eqnarray}$$

It must be stressed that $\unicode[STIX]{x1D702}_{c}$ relies on a volumetric integration. In effect, $\unicode[STIX]{x1D702}_{c}$ will depend on the radiative flux strength and its vertical distribution, the background stratification (density/temperature structure), but also on the basin topography. Basin topography can boost the formation of cross-shore buoyancy gradients and consequently induce horizontal flows. Still, if there is little horizontal variability in the density field and the density at the lowest measured point used to estimate $\langle E_{b}\rangle _{N_{p}}$ is larger than the highest density in the unstable layer near the ice interface, then the introduced field-based method and the theoretical formulation in § 3.1 would be equivalent. In contrast, undersampling in a system with strong horizontal variability can lead to sub- or overestimation of the global mixing efficiency rate (Davies Wykes et al. Reference Davies Wykes, Hughes and Dalziel2015). Thus, the better the spatial resolution of the observations, the more robust the estimation of $\unicode[STIX]{x1D702}_{c}$ is, assuming that the system is closed and/or neglects other forcing mechanisms.

This paper does not aim to give practical information about techniques required for surveying under ice-covered waterbodies. A detailed analysis of under-ice measurements and the errors resulting from sampling is presented by Bouffard et al. (Reference Bouffard, Zdorovennova, Bogdanov, Efremova, Lavanchy, Palshin, Terzhevik, Råman Vinnå, Volkov, Wüest, Zdorovennov and Ulloa2018).

Figure 2.

Observations in the ice-covered Lake Onego (Russia) from 12th to 20th of March 2017. (a) Solar irradiance $H_{0}$ in the water column, (b) density field, $\unicode[STIX]{x1D70C}(z,t)$ , relative to the density at the ice–water interface, $\unicode[STIX]{x1D70C}_{0}$ . (c) Background potential energy $\langle E_{b}\rangle _{Np}$ (blue line, left vertical axis), linear fitting (black dashed line, left vertical axis), diurnally averaged $\langle E_{b}\rangle _{N_{p}}$ (green circles, left vertical axis) and radiative energy, $E_{r,\unicode[STIX]{x1D70F}}$ (red squares, right vertical axis).