### 2.12.4.2. Surface Flux Solution[¶](#surface-flux-solution "Permalink to this headline")

Conservation of energy at the lake surface requires

(2.12.12)[¶](#equation-12-7 "Permalink to this equation")\\\[\\beta \\vec{S}\_{g} -\\vec{L}\_{g} -H\_{g} -\\lambda E\_{g} -G=0\\\]

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

(2.12.13)[¶](#equation-12-8 "Permalink to this equation")\\\[\\begin{split}\\lambda =\\left\\{\\begin{array}{l} {\\lambda \_{sub} \\qquad T\_{g} \\le T\_{f} } \\\\ {\\lambda \_{vap} \\qquad T\_{g} >T\_{f} } \\end{array}\\right\\}.\\end{split}\\\]

The sensible heat flux (W m\-2) is

(2.12.14)[¶](#equation-12-9 "Permalink to this equation")\\\[H\_{g} =-\\rho \_{atm} C\_{p} \\frac{\\left(\\theta \_{atm} -T\_{g} \\right)}{r\_{ah} }\\\]

where \\(\\rho \_{atm}\\) is the density of moist air (kg m\-3) (Chapter [2.5](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#rst-momentum-sensible-heat-and-latent-heat-fluxes)), \\(C\_{p}\\) is the specific heat capacity of air (J kg\-1 K\-1) ([Table 2.2.7](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#table-physical-constants)), \\(\\theta \_{atm}\\) is the atmospheric potential temperature (K) (Chapter [2.5](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#rst-momentum-sensible-heat-and-latent-heat-fluxes)), \\(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](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#monin-obukhov-similarity-theory)).

The water vapor flux (kg m\-2 s\-1) is

(2.12.15)[¶](#equation-12-10 "Permalink to this equation")\\\[E\_{g} =-\\frac{\\rho \_{atm} \\left(q\_{atm} -q\_{sat}^{T\_{g} } \\right)}{r\_{aw} }\\\]

where \\(q\_{atm}\\) is the atmospheric specific humidity (kg kg\-1) (section [2.2.3.1](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#atmospheric-coupling)), \\(q\_{sat}^{T\_{g} }\\) is the saturated specific humidity (kg kg\-1) (section [2.5.5](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#saturation-vapor-pressure)) 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](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#monin-obukhov-similarity-theory)).

The zonal and meridional momentum fluxes are

(2.12.16)[¶](#equation-12-11 "Permalink to this equation")\\\[\\tau \_{x} =-\\rho \_{atm} \\frac{u\_{atm} }{r\_{atm} }\\\]

(2.12.17)[¶](#equation-12-12 "Permalink to this equation")\\\[\\tau \_{y} =-\\rho \_{atm} \\frac{v\_{atm} }{r\_{atm} }\\\]

where \\(u\_{atm}\\) and \\(v\_{atm}\\) are the zonal and meridional atmospheric winds (m s\-1) (section [2.2.3.1](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#atmospheric-coupling)), and \\(r\_{am}\\) is the aerodynamic resistance for momentum (s m\-1) (section [2.5.1](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#monin-obukhov-similarity-theory)).

The heat flux into the lake surface \\(G\\) (W m\-2) is

(2.12.18)[¶](#equation-12-13 "Permalink to this equation")\\\[G=\\frac{2\\lambda \_{T} }{\\Delta z\_{T} } \\left(T\_{g} -T\_{T} \\right)\\\]

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](#eddy-diffusivity-and-thermal-conductivities)), as evaluated in the top lake layer at the previous timestep, where \\(\\lambda \_{liq}\\) is the thermal conductivity of water ([Table 2.2.7](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#table-physical-constants)). 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

(2.12.19)[¶](#equation-12-14 "Permalink to this equation")\\\[\\vec{S}\_{g} =\\sum \_{\\Lambda }S\_{atm} \\, \\downarrow \_{\\Lambda }^{\\mu } \\left(1-\\alpha \_{g,\\, \\Lambda }^{\\mu } \\right) +S\_{atm} \\, \\downarrow \_{\\Lambda } \\left(1-\\alpha \_{g,\\, \\Lambda } \\right)\\\]

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](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#atmospheric-coupling)), and \\(\\alpha \_{g,\\, \\Lambda }^{\\mu }\\) and \\(\\alpha \_{g,\\, \\mu }\\) are the direct beam and diffuse lake albedos (section [2.12.3](#surface-albedo-lake)).

The net emitted longwave radiation is

(2.12.20)[¶](#equation-12-15 "Permalink to this equation")\\\[\\vec{L}\_{g} =L\_{g} \\, \\uparrow -L\_{atm} \\, \\downarrow\\\]

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](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#atmospheric-coupling)). The upward longwave radiation from the surface is

(2.12.21)[¶](#equation-12-16 "Permalink to this equation")\\\[L\\, \\uparrow =\\left(1-\\varepsilon \_{g} \\right)L\_{atm} \\, \\downarrow +\\varepsilon \_{g} \\sigma \\left(T\_{g}^{n} \\right)^{4} +4\\varepsilon \_{g} \\sigma \\left(T\_{g}^{n} \\right)^{3} \\left(T\_{g}^{n+1} -T\_{g}^{n} \\right)\\\]

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](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#table-physical-constants)), 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

(2.12.22)[¶](#equation-12-17 "Permalink to this equation")\\\[\\Delta T\_{g} =\\frac{\\beta \\overrightarrow{S}\_{g} -\\overrightarrow{L}\_{g} -H\_{g} -\\lambda E\_{g} -G}{\\frac{\\partial \\overrightarrow{L}\_{g} }{\\partial T\_{g} } +\\frac{\\partial H\_{g} }{\\partial T\_{g} } +\\frac{\\partial \\lambda E\_{g} }{\\partial T\_{g} } +\\frac{\\partial G}{\\partial T\_{g} } }\\\]

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

(2.12.23)[¶](#equation-12-18 "Permalink to this equation")\\\[T\_{g}^{n+1} =\\frac{\\beta \\overrightarrow{S}\_{g} -\\overrightarrow{L}\_{g} -H\_{g} -\\lambda E\_{g} -G+T\_{g}^{n} \\left(\\frac{\\partial \\overrightarrow{L}\_{g} }{\\partial T\_{g} } +\\frac{\\partial H\_{g} }{\\partial T\_{g} } +\\frac{\\partial \\lambda E\_{g} }{\\partial T\_{g} } +\\frac{\\partial G}{\\partial T\_{g} } \\right)}{\\frac{\\partial \\overrightarrow{L}\_{g} }{\\partial T\_{g} } +\\frac{\\partial H\_{g} }{\\partial T\_{g} } +\\frac{\\partial \\lambda E\_{g} }{\\partial T\_{g} } +\\frac{\\partial G}{\\partial T\_{g} } }\\\]

where the partial derivatives are

(2.12.24)[¶](#equation-12-19 "Permalink to this equation")\\\[\\frac{\\partial \\overrightarrow{L}\_{g} }{\\partial T\_{g} } =4\\varepsilon \_{g} \\sigma \\left(T\_{g}^{n} \\right)^{3} ,\\\]

(2.12.25)[¶](#equation-12-20 "Permalink to this equation")\\\[\\frac{\\partial H\_{g} }{\\partial T\_{g} } =\\frac{\\rho \_{atm} C\_{p} }{r\_{ah} } ,\\\]

(2.12.26)[¶](#equation-12-21 "Permalink to this equation")\\\[\\frac{\\partial \\lambda E\_{g} }{\\partial T\_{g} } =\\frac{\\lambda \\rho \_{atm} }{r\_{aw} } \\frac{dq\_{sat}^{T\_{g} } }{dT\_{g} } ,\\\]

(2.12.27)[¶](#equation-12-22 "Permalink to this equation")\\\[\\frac{\\partial G}{\\partial T\_{g} } =\\frac{2\\lambda \_{T} }{\\Delta z\_{T} } .\\\]

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](#surface-properties-lake).

1.  An initial guess for the wind speed \\(V\_{a}\\) including the convective velocity \\(U\_{c}\\) is obtained from [(2.5.24)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-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)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-50)) and \\(U\_{c} =0.5\\) for unstable conditions (\\(\\theta \_{v,\\, atm} -\\theta \_{v,\\, s} <0\\)).
    
2.  An initial guess for the Monin-Obukhov length \\(L\\) is obtained from the bulk Richardson number using [(2.5.46)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-46) and [(2.5.48)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-48).
    
3.  The following system of equations is iterated four times:
    
4.  Heat of vaporization / sublimation \\(\\lambda\\) ([(2.12.13)](#equation-12-8))
    
5.  Thermal conductivity \\(\\lambda \_{T}\\) (above)
    
6.  Friction velocity \\(u\_{\*}\\) ([(2.5.32)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-32), [(2.5.33)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-33), [(2.5.34)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-34), [(2.5.35)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-35))
    
7.  Potential temperature scale \\(\\theta \_{\*}\\) ([(2.5.37)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-37), [(2.5.38)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-38), [(2.5.39)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-39), [(2.5.40)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-40))
    
8.  Humidity scale \\(q\_{\*}\\) ([(2.5.41)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-41), [(2.5.42)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-42), [(2.5.43)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-43), [(2.5.44)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-44))
    
9.  Aerodynamic resistances \\(r\_{am}\\), \\(r\_{ah}\\), and \\(r\_{aw}\\) ([(2.5.55)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-55), [(2.5.56)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-56), [(2.5.57)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-57))
    
10.  Lake surface temperature \\(T\_{g}^{n+1}\\) ([(2.12.23)](#equation-12-18))
    
11.  Heat of vaporization / sublimation \\(\\lambda\\) ([(2.12.13)](#equation-12-8))
    
12.  Sensible heat flux \\(H\_{g}\\) is updated for \\(T\_{g}^{n+1}\\) ([(2.12.14)](#equation-12-9))
    
13.  Water vapor flux \\(E\_{g}\\) is updated for \\(T\_{g}^{n+1}\\) as
    
    (2.12.28)[¶](#equation-12-23 "Permalink to this equation")\\\[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)](#equation-12-23) is the change in saturated specific humidity due to the change in \\(T\_{g}\\) between iterations.

1.  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](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#monin-obukhov-similarity-theory)).
    
2.  Virtual potential temperature scale \\(\\theta \_{v\*}\\) ([(2.5.17)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-17))
    
3.  Wind speed including the convective velocity, \\(V\_{a}\\) ([(2.5.24)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-24))
    
4.  Monin-Obukhov length \\(L\\) ([(2.5.49)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Fluxes/CLM50_Tech_Note_Fluxes.html#equation-5-49))
    
5.  Roughness lengths (section [2.12.4.1](#surface-properties-lake)).
    

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)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/References/CLM50_Tech_Note_References.html#subinetal2012a)).

(2.12.29)[¶](#equation-12-24 "Permalink to this equation")\\\[\\begin{split}\\begin{array}{l} {{\\rm 1)\\; }T\_{T} \\le T\_{f} <T\_{g} ^{{'} } \\Rightarrow T\_{g} =T\_{f} ,} \\\\ {{\\rm 2)\\; }T\_{T} >T\_{g} ^{{'} } >T\_{m} \\Rightarrow T\_{g} =T\_{T} ,} \\\\ {{\\rm 3)\\; }T\_{m} >T\_{g} ^{{'} } >T\_{T} >T\_{f} \\Rightarrow T\_{g} =T\_{T} } \\end{array}\\end{split}\\\]

where \\(T\_{m}\\) is the temperature of maximum liquid water density, 3.85°C ([Hostetler and Bartlein (1990)](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/References/CLM50_Tech_Note_References.html#hostetlerbartlein1990)). 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)](#equation-12-24) 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)](#equation-12-7)) 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](#boundary-conditions-lake)). 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](#energy-conservation-lake)).

