2.12. Lake Model¶
The lake model, denoted the Lake, Ice, Snow, and Sediment Simulator (LISSS), is from Subin et al. (2012a). It includes extensive modifications to the lake code of Zeng et al. (2002) used in CLM versions 2 through 4, which utilized concepts from the lake models of Bonan (1996), Henderson-Sellers (1985), Henderson-Sellers (1986), Hostetler and Bartlein (1990), and the coupled lake-atmosphere model of Hostetler et al. (1993), Hostetler et al. (1993). Lakes have spatially variable depth prescribed in the surface data (section External Data); the surface data optionally includes lake optical extinction coeffient and horizontal fetch, currently only used for site simulations. Lake physics includes freezing and thawing in the lake body, resolved snow layers, and “soil” and bedrock layers below the lake body. Temperatures and ice fractions are simulated for \(N_{levlak} =10\) layers (for global simulations) or \(N_{levlak} =25\) (for site simulations) with discretization described in section 2.12.1. Lake albedo is described in section 2.12.3. Lake surface fluxes (section 2.12.4) generally follow the formulations for non-vegetated surfaces, including the calculations of aerodynamic resistances (section 2.5.2); however, the lake surface temperature \(T_{g}\) (representing an infinitesimal interface layer between the top resolved lake layer and the atmosphere) is solved for simultaneously with the surface fluxes. After surface fluxes are evaluated, temperatures are solved simultaneously in the resolved snow layers (if present), the lake body, and the soil and bedrock, using the ground heat flux G as a top boundary condition. Snow, soil, and bedrock models generally follow the formulations for non-vegetated surfaces (Chapter 2.6), with modifications described below.
2.12.1. Vertical Discretization¶
Currently, there is one lake modeled in each grid cell (with prescribed or assumed depth d, extinction coefficient \(\eta\), and fetch f), although this could be modified with changes to the CLM subgrid decomposition algorithm in future model versions. As currently implemented, the lake consists of 0-5 snow layers; water and ice layers (10 for global simulations and 25 for site simulations) comprising the “lake body;” 10 “soil” layers; and 5 bedrock layers. Each lake body layer has a fixed water mass (set by the nominal layer thickness and the liquid density), with frozen mass-fraction I a state variable. Resolved snow layers are present if the snow thickness \(z_{sno} \ge s_{\min }\), where smin = 4 cm by default, and is adjusted for model timesteps other than 1800 s in order to maintain numerical stability (section 2.12.6.5). For global simulations with 10 body layers, the default (50 m lake) body layer thicknesses are given by: \(\Delta z_{i}\) of 0.1, 1, 2, 3, 4, 5, 7, 7, 10.45, and 10.45 m, with node depths \(z_{i}\) located at the center of each layer (i.e., 0.05, 0.6, 2.1, 4.6, 8.1, 12.6, 18.6, 25.6, 34.325, 44.775 m). For site simulations with 25 layers, the default thicknesses are (m): 0.1 for layer 1; 0.25 for layers 2-5; 0.5 for layers 6-9; 0.75 for layers 10-13; 2 for layers 14-15; 2.5 for layers 16-17; 3.5 for layers 18-21; and 5.225 for layers 22-25. For lakes with depth d \(\neq\) 50 m and d \(\ge\) 1 m, the top layer is kept at 10 cm and the other 9 layer thicknesses are adjusted to maintain fixed proportions. For lakes with d \(<\) 1 m, all layers have equal thickness. Thicknesses of snow, soil, and bedrock layers follow the scheme used over non-vegetated surfaces (Chapter 2.6), with modifications to the snow layer thickness rules to keep snow layers at least as thick as smin (section 2.12.6.5).
2.12.2. External Data¶
As discussed in Subin et al. (2012a, b), the Global Lake and Wetland Database (Lehner and Doll 2004) is currently used to prescribe lake fraction in each land model grid cell, for a total of 2.3 million km-2. As in Subin et al. (2012a, b), the Kourzeneva et al. (2012) global gridded dataset is currently used to estimate a mean lake depth in each grid cell, based on interpolated compilations of geographic information.
2.12.3. Surface Albedo¶
For direct radiation, the albedo a for lakes with ground temperature \({T}_{g}\) (K) above freezing is given by (Pivovarov, 1972)
where z is the zenith angle. For diffuse radiation, the expression in eq. is integrated over the full sky to yield a = 0.10.
For frozen lakes without resolved snow layers, the albedo at cold temperatures a0 is 0.60 for visible and 0.40 for near infrared radiation. As the temperature at the ice surface, \({T}_{g}\), approaches freezing [ \({T}_{f}\) (K) (Table 2.2.7)], the albedo is relaxed towards 0.10 based on Mironov et al. (2010):
where a is restricted to be no less than that given in (2.12.1).
For frozen lakes with resolved snow layers, the reflectance of the ice surface is fixed at a0, and the snow reflectance is calculated as over non-vegetated surfaces (Chapter 2.3). These two reflectances are combined to obtain the snow-fraction-weighted albedo as in over non-vegetated surfaces (Chapter 2.3).
2.12.4. Surface Fluxes and Surface Temperature¶
2.12.4.1. Surface Properties¶
The fraction of shortwave radiation absorbed at the surface, \(\beta\), depends on the lake state. If resolved snow layers are present, then \(\beta\) is set equal to the absorption fraction predicted by the snow-optics submodel (Chapter 2.3) for the top snow layer. Otherwise, \(\beta\) is set equal to the near infrared fraction of the shortwave radiation reaching the surface simulated by the atmospheric model or atmospheric data model used for offline simulations (Chapter 2.32). The remainder of the shortwave radiation fraction (1 \({-}\) \(\beta\)) is absorbed in the lake body or soil as described in section 2.12.5.5.
The surface roughnesses are functions of the lake state and atmospheric forcing.
For unfrozen lakes (\(T_{g} > T_{f}\)), \(z_{0m}\) is given by (Subin et al. (2012a))
where \(\alpha\) = 0.1, \(\nu\) is the kinematic viscosity of air given below, C is the effective Charnock coefficient given below, \(u_{*}\) is the friction velocity (m/s), and g is the acceleration of gravity (Table 2.2.7). The kinematic viscosity is given by
where \(\nu _{0} =1.51\times 10^{-5} {\textstyle\frac{{\rm m}^{{\rm 2}} }{{\rm s}}}\) , \(T_{0} ={\rm 293.15\; K}\), \(P_{0} =1.013\times 10^{5} {\rm \; Pa}\) , and \(P_{ref}\) is the pressure at the atmospheric reference height. The Charnock coefficient C is a function of the lake fetch F (m), given in the surface data or set to 25 times the lake depth d by default:
where A and B define the fetch- and depth-limitation, respectively; \(C_{\min } =0.01\) , \(C_{\max } =0.01\), \(\varepsilon =1\) , \(f_{c} =100\) , and u (m s-1) is the atmospheric forcing wind.
The scalar roughness lengths (\(z_{0q}\) for latent heat and \(z_{0h}\) for sensible heat) are given by (Subin et al. 2012a)
where \(R_{0}\) is the near-surface atmospheric roughness Reynolds number, \(k\) is the von Karman constant (Table 2.2.7), \(Pr = 0.713\) is the molecular Prandt number for air at neutral stability, \(Sc = 0.66\) is the Schmidt number for water in air at neutral stability. \(z_{0q}\) and \(z_{0h}\) are restricted to be no smaller than \(1 \times 10^{-10}\).
For frozen lakes ( \(T_{g} \le T_{f}\) ) without resolved snow layers ( \(snl = 0\) ), \(z_{0m} =z_{0m_{ice}} =2.3\times 10^{-3} {\rm m}\) (Meier et al. (2022)).
For frozen lakes with resolved snow layers ( \(snl > 0\) ), the momentum roughness length is evaluated based on accumulated snow melt \(M_{a} {\rm }\) (Meier et al. (2022)). For \(M_{a} >=1\times 10^{-5}\)
where \(M_{a}\) is accumulated snow melt (meters water equivalent), \(b_{1} =1.4\) and \(b_{4} =-0.31\). For \(M_{a} <1\times 10^{-5}\)
Accumulated snow melt \(M_{a}\) at the current time step \(t\) is defined as
where \(M ^{t}_{a}\) and \(M ^{t-1}_{a}\) are the accumulated snowmelt at the current time step and previous time step, respectively (m), \(q ^{t}_{sno} \Delta t\) is the freshly fallen snow (mm), and \(q ^{t}_{snowmelt} \Delta t\) is the melted snow (mm).
For frozen lakes without and with resolved snow layers, an initial guess for the scalar roughness lengths is derived by assuming \(\theta_{*} = 0 {\rm }\) (Meier et al. (2022))
where \(\nu=1.5 \times 10^{-5}\) is the kinematic viscosity of air (m2 s-1), and \(u_{*}\) is the friction velocity in the atmospheric surface layer (m s-1). Thereafter, the scalar roughness lengths are updated within the stability iteration described in section 2.12.4.2 as
where \(\beta\) = 7.2, and \(\theta_{*}\) is the potential temperature scale (section 2.12.4.2).
2.12.4.2. Surface Flux Solution¶
Conservation of energy at the lake surface requires
where \(\vec{S}_{g}\) is the absorbed solar radiation in the lake, \(\beta\) is the fraction absorbed at the surface, \(\vec{L}_{g}\) is the net emitted longwave radiation (+ upwards), \(H_{g}\) is the sensible heat flux (+ upwards), \(E_{g}\) is the water vapor flux (+ upwards), and G is the ground heat flux (+ downwards). All of these fluxes depend implicitly on the temperature at the lake surface \({T}_{g}\). \(\lambda\) converts \(E_{g}\) to an energy flux based on
The sensible heat flux (W m-2) is
where \(\rho _{atm}\) is the density of moist air (kg m-3) (Chapter 2.5), \(C_{p}\) is the specific heat capacity of air (J kg-1 K-1) (Table 2.2.7), \(\theta _{atm}\) is the atmospheric potential temperature (K) (Chapter 2.5), \(T_{g}\) is the lake surface temperature (K) (at an infinitesimal interface just above the top resolved model layer: snow, ice, or water), and \(r_{ah}\) is the aerodynamic resistance to sensible heat transfer (s m-1) (section 2.5.1).
The water vapor flux (kg m-2 s-1) is
where \(q_{atm}\) is the atmospheric specific humidity (kg kg-1) (section 2.2.3.1), \(q_{sat}^{T_{g} }\) is the saturated specific humidity (kg kg-1) (section 2.5.5) at the lake surface temperature \(T_{g}\), and \(r_{aw}\) is the aerodynamic resistance to water vapor transfer (s m-1) (section 2.5.1).
The zonal and meridional momentum fluxes are
where \(u_{atm}\) and \(v_{atm}\) are the zonal and meridional atmospheric winds (m s-1) (section 2.2.3.1), and \(r_{am}\) is the aerodynamic resistance for momentum (s m-1) (section 2.5.1).
The heat flux into the lake surface \(G\) (W m-2) is
where \(\lambda _{T}\) is the thermal conductivity (W m-1 K-1), \(\Delta z_{T}\) is the thickness (m), and \(T_{T}\) is the temperature (K) of the top resolved lake layer (snow, ice, or water). The top thermal conductivity \(\lambda _{T}\) of unfrozen lakes ( \(T_{g} >T_{f}\) ) includes conductivities due to molecular ( \(\lambda _{liq}\) ) and eddy (\(\lambda _{K}\) ) diffusivities (section 2.12.5.4), as evaluated in the top lake layer at the previous timestep, where \(\lambda _{liq}\) is the thermal conductivity of water (Table 2.2.7). For frozen lakes without resolved snow layers, \(\lambda _{T} =\lambda _{ice}\). When resolved snow layers are present, \(\lambda _{T}\) is calculated based on the water content, ice content, and thickness of the top snow layer, as for non-vegetated surfaces.
The absorbed solar radiation \(\vec{S}_{g}\) is
where \(S_{atm} \, \downarrow _{\Lambda }^{\mu }\) and \(S_{atm} \, \downarrow _{\Lambda }\) are the incident direct beam and diffuse solar fluxes (W m-2) and \(\Lambda\) denotes the visible (\(<\) 0.7\(\mu {\rm m}\)) and near-infrared (\(\ge\) 0.7\(\mu {\rm m}\)) wavebands (section 2.2.3.1), and \(\alpha _{g,\, \Lambda }^{\mu }\) and \(\alpha _{g,\, \mu }\) are the direct beam and diffuse lake albedos (section 2.12.3).
The net emitted longwave radiation is
where \(L_{g} \, \uparrow\) is the upward longwave radiation from the surface, \(L_{atm} \, \downarrow\) is the downward atmospheric longwave radiation (section 2.2.3.1). The upward longwave radiation from the surface is
where \(\varepsilon _{g} =0.97\) is the lake surface emissivity, \(\sigma\) is the Stefan-Boltzmann constant (W m-2 K-4) (Table 2.2.7), and \(T_{g}^{n+1} -T_{g}^{n}\) is the difference in lake surface temperature between Newton-Raphson iterations (see below).
The sensible heat \(H_{g}\), the water vapor flux \(E_{g}\) through its dependence on the saturated specific humidity, the net longwave radiation \(\vec{L}_{g}\), and the ground heat flux \(G\), all depend on the lake surface temperature \(T_{g}\). Newton-Raphson iteration is applied to solve for \(T_{g}\) and the surface fluxes as
where \(\Delta T_{g} =T_{g}^{n+1} -T_{g}^{n}\) and the subscript “n” indicates the iteration. Therefore, the surface temperature \(T_{g}^{n+1}\) can be written as
where the partial derivatives are
The fluxes of momentum, sensible heat, and water vapor are solved for simultaneously with lake surface temperature as follows. To begin, \(z_{0m}\) and the scalar roughness lengths are set as described in section 2.12.4.1.
An initial guess for the wind speed \(V_{a}\) including the convective velocity \(U_{c}\) is obtained from (2.5.24) assuming an initial convective velocity \(U_{c} =0\) m s-1 for stable conditions (\(\theta _{v,\, atm} -\theta _{v,\, s} \ge 0\) as evaluated from (2.5.50)) and \(U_{c} =0.5\) for unstable conditions (\(\theta _{v,\, atm} -\theta _{v,\, s} <0\)).
An initial guess for the Monin-Obukhov length \(L\) is obtained from the bulk Richardson number using (2.5.46) and (2.5.48).
The following system of equations is iterated four times:
Heat of vaporization / sublimation \(\lambda\) ((2.12.13))
Thermal conductivity \(\lambda _{T}\) (above)
Friction velocity \(u_{*}\) ((2.5.32), (2.5.33), (2.5.34), (2.5.35))
Potential temperature scale \(\theta _{*}\) ((2.5.37), (2.5.38), (2.5.39), (2.5.40))
Humidity scale \(q_{*}\) ((2.5.41), (2.5.42), (2.5.43), (2.5.44))
Aerodynamic resistances \(r_{am}\), \(r_{ah}\), and \(r_{aw}\) ((2.5.55), (2.5.56), (2.5.57))
Lake surface temperature \(T_{g}^{n+1}\) ((2.12.23))
Heat of vaporization / sublimation \(\lambda\) ((2.12.13))
Sensible heat flux \(H_{g}\) is updated for \(T_{g}^{n+1}\) ((2.12.14))
Water vapor flux \(E_{g}\) is updated for \(T_{g}^{n+1}\) as
(2.12.28)¶\[E_{g} =-\frac{\rho _{atm} }{r_{aw} } \left[q_{atm} -q_{sat}^{T_{g} } -\frac{\partial q_{sat}^{T_{g} } }{\partial T_{g} } \left(T_{g}^{n+1} -T_{g}^{n} \right)\right]\]
where the last term on the right side of equation (2.12.28) is the change in saturated specific humidity due to the change in \(T_{g}\) between iterations.
Saturated specific humidity \(q_{sat}^{T_{g} }\) and its derivative \(\frac{dq_{sat}^{T_{g} } }{dT_{g} }\) are updated for \(T_{g}^{n+1}\) (section 2.5.1).
Virtual potential temperature scale \(\theta _{v*}\) ((2.5.17))
Wind speed including the convective velocity, \(V_{a}\) ((2.5.24))
Monin-Obukhov length \(L\) ((2.5.49))
Roughness lengths (section 2.12.4.1).
Once the four iterations for lake surface temperature have been yielded a tentative solution \(T_{g} ^{{'} }\), several restrictions are imposed in order to maintain consistency with the top lake model layer temperature \(T_{T}\) (Subin et al. (2012a)).
where \(T_{m}\) is the temperature of maximum liquid water density, 3.85°C (Hostetler and Bartlein (1990)). The first condition requires that, if there is any snow or ice present, the surface temperature is restricted to be less than or equal to freezing. The second and third conditions maintain convective stability in the top lake layer.
If equation (2.12.29) is applied, the turbulent fluxes \(H_{g}\) and \(E_{g}\) are re-evaluated. The emitted longwave radiation and the momentum fluxes are re-evaluated in any case. The final ground heat flux \(G\) is calculated from the residual of the energy balance (equation (2.12.12)) in order to precisely conserve energy. This ground heat flux is taken as a prescribed flux boundary condition for the lake temperature solution (section 2.12.5.3). A check is included at each timestep to insure that energy balance is obeyed to within 0.1 W m-2 (see 2.12.5.10).
2.12.5. Lake Temperature¶
2.12.5.1. Introduction¶
The (optional-) snow, lake body (water and/or ice), soil, and bedrock system is unified for the lake temperature solution. The governing equation, similar to that for the snow-soil-bedrock system for vegetated land units (Chapter 2.6), is
where \(\tilde{c}_{v}\) is the volumetric heat capacity (J m-3 K-1), \(t\) is time (s), T is the temperature (K), \(\tau\) is the thermal conductivity (W m-1 K-1), and \(\phi\) is the solar radiation (W m-2) penetrating to depth z (m). The system is discretized into N layers, where
\(n_{sno}\) is the number of actively modeled snow layers at the current timestep (Chapter 2.8), and \(N_{levgrnd}\) is as for vegetated land units (Chapter 2.6). Energy is conserved as
where \(\tilde{c}_{v,j} (t)\)is the volumetric heat capacity of the jth layer (section 2.12.5.5), \(L_{j} (t)\)is the latent heat of fusion per unit volume of the jth layer (proportional to the mass of liquid water present), and the right-hand side represents the net influx of energy to the lake system. Note that \(\tilde{c}_{v,j} (t)\) can only change due to phase change (except for changing snow layer mass, which, apart from energy required to melt snow, represents an untracked energy flux in the land model, along with advected energy associated with water flows in general), and this is restricted to occur at \(T_{j} =T_{f}\) in the snow-lake-soil system, allowing eq. to be precisely enforced and justifying the exclusion of \(c_{v,j}\) from the time derivative in eq..
2.12.5.2. Overview of Changes from CLM4¶
Thermal conductivities include additional eddy diffusivity, beyond the Hostetler and Bartlein (1990) formulation, due to unresolved processes (Fang and Stefan 1996; Subin et al. (2012a)). Lake water is now allowed to freeze by an arbitrary fraction for each layer, which releases latent heat and changes thermal properties. Convective mixing occurs for all lakes, even if frozen. Soil and bedrock are included beneath the lake. The full snow model is used if the snow thickness exceeds a threshold; if there are resolved snow layers, radiation transfer is predicted by the snow-optics submodel (Chapter 2.3), and the remaining radiation penetrating the bottom snow layer is absorbed in the top layer of lake ice; conversely, if there are no snow layers, the solar radiation penetrating the bottom lake layer is absorbed in the top soil layer. The lakes have variable depth, and all physics is assumed valid for arbitrary depth, except for a depth-dependent enhanced mixing (section 2.12.5.4). Finally, a previous sign error in the calculation of eddy diffusivity (specifically, the Brunt-Väisälä frequency term; eq. ) was corrected.
2.12.5.3. Boundary Conditions¶
The top boundary condition, imposed at the top modeled layer \(i=j_{top}\), where \(j_{top} =-n_{sno} +1\), is the downwards surface flux G defined by the energy flux residual during the surface temperature solution (section 2.12.5.3). The bottom boundary condition, imposed at \(i=N_{levlak} +N_{levgrnd}\), is zero flux. The 2-m windspeed \(u_{2}\) (m s-1) is used in the calculation of eddy diffusivity:
where \(u_{*}\) is the friction velocity calculated in section 2.12.5.3 and k is the von Karman constant (Table 2.2.7).
2.12.5.4. Eddy Diffusivity and Thermal Conductivities¶
The total eddy diffusivity \(K_{W}\) (m2 s-1) for liquid water in the lake body is given by (Subin et al. (2012a))
where \(\kappa _{e}\) is due to wind-driven eddies (Hostetler and Bartlein (1990)), \(K_{ed}\) is a modest enhanced diffusivity intended to represent unresolved mixing processes (Fang and Stefan 1996), \(\kappa _{m} =\frac{\lambda _{liq} }{c_{liq} \rho _{liq} }\) is the molecular diffusivity of water (given by the ratio of its thermal conductivity (W m-1 K-1) to the product of its heat capacity (J kg-1 K-1) and density (kg m-3), values given in Table 2.2.7), and \(m_{d}\) (unitless) is a factor which increases the overall diffusivity for large lakes, intended to represent 3-dimensional mixing processes such as caused by horizontal temperature gradients. As currently implemented,
where d is the lake depth.
The wind-driven eddy diffusion coefficient \(\kappa _{e,\, i}\) (m2 s-1) for layers \(1\le i\le N_{levlak}\) is
where \(P_{0} =1\) is the neutral value of the turbulent Prandtl number, \(z_{i}\) is the node depth (m), the surface friction velocity (m s-1) is \(w^{*} =0.0012u_{2}\), and \(k^{*}\) varies with latitude \(\phi\) as \(k^{*} =6.6u_{2}^{-1.84} \sqrt{\left|\sin \phi \right|}\). For the bottom layer, \(\kappa _{e,\, N_{levlak} } =\kappa _{e,N_{levlak} -1\, }\). As in Hostetler and Bartlein (1990), the 2-m wind speed \(u_{2}\) (m s-1) (eq. ) is used to evaluate \(w^{*}\) and \(k^{*}\) rather than the 10-m wind used by Henderson-Sellers (1985).
The Richardson number is
where
and \(g\) is the acceleration due to gravity (m s-2) (Table 2.2.7), \(\rho _{i}\) is the density of water (kg m-3), and \(\frac{\partial \rho }{\partial z}\) is approximated as \(\frac{\rho _{i+1} -\rho _{i} }{z_{i+1} -z_{i} }\). Note that because here, z is increasing downwards (unlike in Hostetler and Bartlein (1990)), eq. contains no negative sign; this is a correction from CLM4. The density of water is (Hostetler and Bartlein (1990))
The enhanced diffusivity \(K_{ed}\) is given by (Fang and Stefan 1996)
where \(N^{2}\) is calculated as in eq. except for the minimum value imposed in.
The thermal conductivity for the liquid water portion of lake body layer i, \(\tau _{liq,i}\) (W m-1 K-1) is given by
The thermal conductivity of the ice portion of lake body layer i, \(\tau _{ice,eff}\) (W m-1 K-1), is constant among layers, and is given by
where \(\tau _{ice}\) (Table 2.2.7) is the nominal thermal conductivity of ice: \(\tau _{ice,eff}\) is adjusted for the fact that the nominal model layer thicknesses remain constant even while the physical ice thickness exceeds the water thickness.
The overall thermal conductivity \(\tau _{i}\) for layer i with ice mass-fraction \(I_{i}\) is the harmonic mean of the liquid and water fractions, assuming that they will be physically vertically stacked, and is given by
The thermal conductivity of snow, soil, and bedrock layers above and below the lake, respectively, are computed identically to those for vegetated land units (Chapter 2.6), except for the adjustment of thermal conductivity for frost heave or excess ice (Subin et al., 2012a, Supporting Information).
2.12.5.5. Radiation Penetration¶
If there are no resolved snow layers, the surface absorption fraction \(\beta\) is set according to the near-infrared fraction simulated by the atmospheric model. This is apportioned to the surface energy budget (section 2.12.4.1), and thus no additional radiation is absorbed in the top \(z_{a}\) (currently 0.6 m) of unfrozen lakes, for which the light extinction coefficient \(\eta\) (m-1) varies between lake columns (eq. ). For frozen lakes (\(T_{g} \le T_{f}\) ), the remaining \(\left(1-\beta \right)\vec{S}_{g}\) fraction of surface absorbed radiation that is not apportioned to the surface energy budget is absorbed in the top lake body layer. This is a simplification, as lake ice is partially transparent. If there are resolved snow layers, then the snow optics submodel (Chapter 2.3) is used to calculate the snow layer absorption (except for the absorption predicted for the top layer by the snow optics submodel, which is assigned to the surface energy budget), with the remainder penetrating snow layers absorbed in the top lake body ice layer.
For unfrozen lakes, the solar radiation remaining at depth \(z>z_{a}\) in the lake body is given by
For all lake body layers, the flux absorbed by the layer i, \(\phi _{i}\) , is
The argument of each exponent is constrained to be non-negative (so \(\phi _{i}\) = 0 for layers contained within \({z}_{a}\)). The remaining flux exiting the bottom of layer \(i=N_{levlak}\) is absorbed in the top soil layer.
The light extinction coefficient \(\eta\) (m-1), if not provided as external data, is a function of depth d (m) (Subin et al. (2012a)):
2.12.5.6. Heat Capacities¶
The vertically-integrated heat capacity for each lake layer, \(\text{c}_{v,i}\) (J m-2) is determined by the mass-weighted average over the heat capacities for the water and ice fractions:
Note that the density of water is used for both ice and water fractions, as the thickness of the layer is fixed.
The total heat capacity \(c_{v,i}\) for each soil, snow, and bedrock layer (J m-2) is determined as for vegetated land units (Chapter 2.6), as the sum of the heat capacities for the water, ice, and mineral constituents.
2.12.5.7. Crank-Nicholson Solution¶
The solution method for thermal diffusion is similar to that used for soil (Chapter 2.6), except that the lake body layers are sandwiched between the snow and soil layers (section 2.12.5.1), and radiation flux is absorbed throughout the lake layers. Before solution, layer temperatures \(T_{i}\) (K), thermal conductivities \(\tau _{i}\) (W m-1 K-1), heat capacities \(c_{v,i}\) (J m-2), and layer and interface depths from all components are transformed into a uniform set of vectors with length \(N=n_{sno} +N_{levlak} +N_{levgrnd}\) and consistent units to simplify the solution. Thermal conductivities at layer interfaces are calculated as the harmonic mean of the conductivities of the neighboring layers:
where \(\lambda _{i}\) is the conductivity at the interface between layer i and layer i + 1, \(z_{i}\) is the depth of the node of layer i, and \(\hat{z}_{i}\) is the depth of the interface below layer i. Care is taken at the boundaries between snow and lake and between lake and soil. The governing equation is discretized for each layer as
where superscripts n + 1 and n denote values at the end and beginning of the timestep \(\Delta t\), respectively, \(F_{i}\) (W m-2) is the downward heat flux at the bottom of layer i, and \(\phi _{i}\) is the solar radiation absorbed in layer i.
Eq. is solved using the semi-implicit Crank-Nicholson Method, resulting in a tridiagonal system of equations:
The fluxes \(F_{i}\) are defined as follows: for the top layer, \(F_{j_{top} -1} =2G;a_{j_{top} } =0\), where G is defined as in section 2.12.5.3 (the factor of 2 merely cancels out the Crank-Nicholson 0.5 in the equation for \(r_{j_{top} }\) ). For the bottom layer, \(F_{N_{levlak} +N_{levgrnd} } =0\). For all other layers:
2.12.5.8. Phase Change¶
Phase change in the lake, snow, and soil is done similarly to that done for the soil and snow for vegetated land units (Chapter 2.6), except without the allowance for freezing point depression in soil underlying lakes. After the heat diffusion is calculated, phase change occurs in a given layer if the temperature is below freezing and liquid water remains, or if the temperature is above freezing and ice remains.
If melting occurs, the available energy for melting, \(Q_{avail}\) (J m-2), is computed as
where \(T_{i}\) is the temperature of the layer after thermal diffusion (section 2.12.5.7), and \(c_{v,i}\) is as calculated in section 2.12.5.6. The mass of melt in the layer M (kg m-2) is given by
where \(H_{fus}\) (J kg-1) is the latent heat of fusion of water (Table 2.2.7), and \(M_{ice}\) is the mass of ice in the layer: \(I_{i} \rho _{liq} \Delta z_{i}\) for a lake body layer, or simply the soil / snow ice content state variable (\(w_{ice}\) ) for a soil / snow layer. The heat remainder, \(Q_{rem}\) is given by
Finally, the mass of ice in the layer \(M_{ice}\) is adjusted downwards by \(M\), and the temperature \(T_{i}\) of the layer is adjusted to
where \(c'_{v,i} =c_{v,i} +M\left(c_{liq} -c_{ice} \right)\).
If freezing occurs, \(Q_{avail}\) is again given by but will be negative. The melt \(M\), also negative, is given by
where \(M_{liq}\) is the mass of water in the layer: \(\left(1-I_{i} \right)\rho _{liq} \Delta z_{i}\) for a lake body layer, or the soil / snow water content state variable (\(w_{liq}\) ). The heat remainder \(Q_{rem}\) is given by eq. and will be negative or zero. Finally, \(M_{liq}\) is adjusted downwards by \(-M\) and the temperature is reset according to eq..
In the presence of nonzero snow water \(W_{sno}\) without resolved snow layers over an unfrozen top lake layer, the available energy in the top lake layer \(\left(T_{1} -T_{f} \right)c_{v,1}\) is used to melt the snow. Similar to above, \(W_{sno}\) is either completely melted and the remainder of heat returned to the top lake layer, or the available heat is exhausted and the top lake layer is set to freezing. The snow thickness is adjusted downwards in proportion to the amount of melt, maintaining constant density.
2.12.5.9. Convection¶
Convective mixing is based on Hostetler et al.’s (1993, 1994) coupled lake-atmosphere model, adjusting the lake temperature after diffusion and phase change to maintain a stable density profile. Unfrozen lakes overturn when \(\rho _{i} >\rho _{i+1}\), in which case the layer thickness weighted average temperature for layers 1 to \(i+1\) is applied to layers 1 to \(i+1\) and the densities are updated. This scheme is applied iteratively to layers \(1\le i<N_{levlak} -1\). Unstable profiles occurring at the bottom of the lake (i.e., between layers \(i=N_{levlak} -1\) and \(i=N_{levlak}\) ) are treated separately (Subin et al. (2012a)), as occasionally these can be induced by heat expelled from the sediments (not present in the original Hostetler et al. (1994) model). Mixing proceeds from the bottom upward in this case (i.e., first mixing layers \(i=N_{levlak} -1\) and \(i=N_{levlak}\), then checking \(i=N_{levlak} -2\) and \(i=N_{levlak} -1\) and mixing down to \(i=N_{levlak}\) if needed, and on to the top), so as not to mix in with warmer over-lying layers.
For frozen lakes, this algorithm is generalized to conserve total enthalpy and ice content, and to maintain ice contiguous at the top of the lake. Thus, an additional mixing criterion is added: the presence of ice in a layer that is below a layer which is not completely frozen. When this occurs, these two lake layers and all those above mix. Total enthalpy Q is conserved as
Once the average ice fraction \(I_{av}\) is calculated from
the temperatures are calculated. A separate temperature is calculated for the frozen (\(T_{froz}\) ) and unfrozen (\(T_{unfr}\) ) fractions of the mixed layers. If the total heat content Q is positive (e.g. some layers will be above freezing), then the extra heat is all assigned to the unfrozen layers, while the fully frozen layers are kept at freezing. Conversely, if \(Q < 0\), the heat deficit will all be assigned to the ice, and the liquid layers will be kept at freezing. For the layer that contains both ice and liquid (if present), a weighted average temperature will have to be calculated.
If \(Q > 0\), then \(T_{froz} =T_{f}\), and \(T_{unfr}\) is given by
If \(Q < 0\), then \(T_{unfr} =T_{f}\), and \(T_{froz}\) is given by
The ice is lumped together at the top. For each lake layer j from 1 to i + 1, the ice fraction and temperature are set as follows, where \(Z_{j} =\sum _{m=1}^{j}\Delta z_{m}\) :
If \(Z_{j} \le Z_{i+1} I_{av}\), then \(I_{j} =1\) and \(T_{j} =T_{froz}\).
Otherwise, if \(Z_{j-1} <Z_{i+1} I_{av}\), then the layer will contain both ice and water. The ice fraction is given by \(I_{j} =\frac{Z_{i+1} I_{av} -Z_{j-1} }{\Delta z_{j} }\). The temperature is set to conserve the desired heat content that would be present if the layer could have two temperatures, and then dividing by the heat capacity of the layer to yield
(2.12.61)¶\[T_{j} =\frac{T_{froz} I_{j} c_{ice} +T_{unfr} \left(1-I_{j} \right)c_{liq} }{I_{j} c_{ice} +\left(1-I_{j} \right)c_{liq} } .\]Otherwise, \(I_{j} =0\) and \(T_{j} =T_{unfr}\) .
2.12.5.10. Energy Conservation¶
To check energy conservation, the left-hand side of equation (2.12.32) is re-written to yield the total enthalpy of the lake system (J m-2) \(H_{tot}\) :
where \(M_{liq,i}\) is the water mass of the ith layer (similar to section 2.12.5.8), and \(W_{sno,bulk}\) is the mass of snow-ice not present in resolved snow layers. This expression is evaluated once at the beginning and once at the end of the timestep (re-evaluating each \(c_{v,i}\) ), and the change is compared with the net surface energy flux to yield the error flux \(E_{soi}\) (W m-2):
If \(\left|E_{soi} \right|<0.1\)W m-2, it is subtracted from the sensible heat flux and added to G. Otherwise, the model is aborted.
2.12.6. Lake Hydrology¶
2.12.6.1. Overview¶
Hydrology is done similarly to other impervious non-vegetated columns (e.g., glaciers) where snow layers may be resolved but infiltration into the permanent ground is not allowed. The water mass of lake columns is currently maintained constant, aside from overlying snow. The water budget is balanced with \(q_{rgwl}\) (eq.; kg m-2 s-1), a generalized runoff term for impervious land units that may be negative.
There are some modifications to the soil and snow parameterizations as compared with the soil in vegetated land units, or the snow overlying other impervious columns. The soil can freeze or thaw, with the allowance for frost heave (or the initialization of excess ice) (sections 2.12.5.4 and 2.12.5.8), but no air-filled pore space is allowed in the soil. To preserve numerical stability in the lake model (which uses a slightly different surface flux algorithm than over other non-vegetated land units), two changes are made to the snow model. First, dew or frost is not allowed to be absorbed by a top snow layer which has become completely melted during the timestep. Second, because occasional instabilities occurred during model testing when the Courant–Friedrichs–Lewy (CFL) condition was violated, due to the explicit time-stepping integration of the surface flux solution, resolved snow layers must be a minimum of \(s_{\min }\) = 4 cm thick rather than 1 cm when the default timestep of 1800 s is used.
2.12.6.2. Water Balance¶
The total water balance of the system is given by
where \(W_{sno}\) (kg m-2) is the total mass of snow (both liquid and ice, in resolved snow layers or bulk snow), \(w_{liq,i}\) and \(w_{ice,i}\) are the masses of water phases (kg m-2) in soil layer i, \(q_{rain}\) and \(q_{sno}\) are the precipitation forcing from the atmosphere (kg m-2 s-1), \(q_{snwcp,\, ice}\) is the ice runoff associated with snow-capping (below), \(E_{g}\) is the ground evaporation (section 2.12.4.2), and \(n_{levsoi}\) is the number of hydrologically active soil layers (as opposed to dry bedrock layers).
2.12.6.3. Precipitation, Evaporation, and Runoff¶
All precipitation reaches the ground, as there is no vegetated fraction. As for other land types, incident snowfall accumulates (with ice mass \(W_{sno}\) and thickness \(z_{sno}\) ) until its thickness exceeds a minimum thickness \(s_{\min }\), at which point a resolved snow layer is initiated, with water, ice, dissolved aerosol, snow-grain radius, etc., state variables tracked by the Snow Hydrology submodel (Chapter 2.8). The density of fresh snow is assigned as for other land types (Chapter 2.8). Solid precipitation is added immediately to the snow, while liquid precipitation is added to snow layers, if they exist, after accounting for dew, frost, and sublimation (below). If \(z_{sno}\) exceeds \(s_{\min }\) after solid precipitation is added but no snow layers are present, a new snow layer is initiated immediately, and then dew, frost, and sublimation are accounted for. Snow-capping is invoked if the snow depth \(z_{sno} >1000{\rm m}\), in which case additional precipitation and frost deposition is added to \(q_{snwcp,\, ice}\).
If there are resolved snow layers, the generalized “evaporation” \(E_{g}\) (i.e., evaporation, dew, frost, and sublimation) is treated as over other land units, except that the allowed evaporation from the ground is unlimited (though the top snow layer cannot lose more water mass than it contains). If there are no resolved snow layers but \(W_{sno} >0\) and \(E_{g} >0\), sublimation \(q_{sub,sno}\) (kg m-2 s-1) will be given by
If \(E_{g} <0,T_{g} \le T_{f}\), and there are no resolved snow layers or the top snow layer is not unfrozen, then the rate of frost production \(q_{frost} =\left|E_{g} \right|\). If \(E_{g} <0\) but the top snow layer has completely thawed during the Phase Change step of the Lake Temperature solution (section 2.12.5.8), then frost (or dew) is not allowed to accumulate (\(q_{frost} =0\)), to insure that the layer is eliminated by the Snow Hydrology (Chapter 2.8) code. (If \(T_{g} >T_{f}\), then no snow is present (section 2.12.4.2), and evaporation or dew deposition is balanced by \(q_{rgwl}\).) The snowpack is updated for frost and sublimation:
If there are resolved snow layers, then this update occurs using the Snow Hydrology submodel (Chapter 2.8). Otherwise, the snow ice mass is updated directly, and \(z_{sno}\) is adjusted by the same proportion as the snow ice (i.e., maintaining the same density), unless there was no snow before adding the frost, in which case the density is assumed to be 250 kg m-3.
2.12.6.4. Soil Hydrology¶
The combined water and ice soil volume fraction in a soil layer \(\theta _{i}\) is given by
If \(\theta _{i} <\theta _{sat,i}\), the pore volume fraction at saturation (as may occur when ice melts), then the liquid water mass is adjusted to
Otherwise, if excess ice is melting and \(w_{liq,i} >\theta _{sat,i} \rho _{liq} \Delta z_{i}\), then the water in the layer is reset to
This allows excess ice to be initialized (and begin to be lost only after the pore ice is melted, which is realistic if the excess ice is found in heterogeneous chunks) but irreversibly lost when melt occurs.
2.12.6.5. Modifications to Snow Layer Logic¶
A thickness difference \(z_{lsa} =s_{\min } -\tilde{s}_{\min }\) adjusts the minimum resolved snow layer thickness for lake columns as compared to non-lake columns. The value of \(z_{lsa}\) is chosen to satisfy the CFL condition for the model timestep. By default, \(\tilde{s}_{\min }\) = 1 cm and \(s_{\min }\) = 4 cm. See Subin et al. (2012a; including Supporting Information) for further discussion.
The rules for combining and sub-dividing snow layers (section 2.8.7) are adjusted for lakes to maintain minimum thicknesses of \(s_{\min }\) and to increase all target layer thicknesses by \(z_{lsa}\). The rules for combining layers are modified by simply increasing layer thickness thresholds by \(z_{lsa}\). The rules for dividing snow layers are contained in a separate subroutine that is modified for lakes, and is a function of the number of layers and the layer thicknesses. There are two types of operations: (a) subdividing layers in half, and (b) shifting some volume from higher layers to lower layers (without increasing the layer number). For subdivisions of type (a), the thickness thresholds triggering subdivision are increased by \(2z_{lsa}\) for lakes. For shifts of type (b), the thickness thresholds triggering the shifts are increased by \(z_{lsa}\). At the end of the modified subroutine, a snow ice and liquid balance check are performed.
In rare instances, resolved snow layers may be present over an unfrozen top lake body layer. In this case, the snow layers may be eliminated if enough heat is present in the top layer to melt the snow: see Subin et al. (2012a, Supporting Information).