### 2.12.5.7. Crank-Nicholson Solution[¶](#crank-nicholson-solution "Permalink to this headline")

The solution method for thermal diffusion is similar to that used for soil (Chapter [2.6](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Soil_Snow_Temperatures/CLM50_Tech_Note_Soil_Snow_Temperatures.html#rst-soil-and-snow-temperatures)), except that the lake body layers are sandwiched between the snow and soil layers (section [2.12.5.1](#introduction-lake)), 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:

(2.12.48)[¶](#equation-12-43 "Permalink to this equation")\\\[\\lambda \_{i} =\\frac{\\tau \_{i} \\tau \_{i+1} \\left(z\_{i+1} -z\_{i} \\right)}{\\tau \_{i} \\left(z\_{i+1} -\\hat{z}\_{i} \\right)+\\tau \_{i+1} \\left(\\hat{z}\_{i} -z\_{i} \\right)} ,\\\]

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

(2.12.49)[¶](#equation-12-44 "Permalink to this equation")\\\[\\frac{c\_{v,i} }{\\Delta t} \\left(T\_{i}^{n+1} -T\_{i}^{n} \\right)=F\_{i-1} -F\_{i} +\\phi \_{i}\\\]

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:

(2.12.50)[¶](#equation-12-45 "Permalink to this equation")\\\[\\begin{split}\\begin{array}{l} {r\_{i} =a\_{i} T\_{i-1}^{n+1} +b\_{i} T\_{i}^{n+1} +cT\_{i+1}^{n+1} ,} \\\\ {a\_{i} =-0.5\\frac{\\Delta t}{c\_{v,i} } \\frac{\\partial F\_{i-1} }{\\partial T\_{i-1}^{n} } ,} \\\\ {b\_{i} =1+0.5\\frac{\\Delta t}{c\_{v,i} } \\left(\\frac{\\partial F\_{i-1} }{\\partial T\_{i-1}^{n} } +\\frac{\\partial F\_{i} }{\\partial T\_{i}^{n} } \\right),} \\\\ {c\_{i} =-0.5\\frac{\\Delta t}{c\_{v,i} } \\frac{\\partial F\_{i} }{\\partial T\_{i}^{n} } ,} \\\\ {r\_{i} =T\_{i}^{n} +0.5\\frac{\\Delta t}{c\_{v,i} } \\left(F\_{i-1} -F\_{i} \\right)+\\frac{\\Delta t}{c\_{v,i} } \\phi \_{i} .} \\end{array}\\end{split}\\\]

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](#boundary-conditions-lake) (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.51)[¶](#equation-12-46 "Permalink to this equation")\\\[F\_{i} =\\lambda \_{i} \\frac{T\_{i} ^{n} -T\_{i+1}^{n} }{z\_{n+1} -z\_{n} } .\\\]

