## 2.6.1. Numerical Solution[¶](#numerical-solution "Permalink to this headline")
------------------------------------------------------------------------------

The soil column is discretized into 25 layers (section [2.2.2](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#vertical-discretization)) where \\(N\_{levgrnd} = 25\\) is the number of soil layers ([Table 2.2.3](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#table-soil-layer-structure)).

The overlying snow pack is modeled with up to twelve layers depending on the total snow depth. The layers from top to bottom are indexed in the Fortran code as \\(i=-4,-3,-2,-1,0\\), which permits the accumulation or ablation of snow at the top of the snow pack without renumbering the layers. Layer \\(i=0\\) is the snow layer next to the soil surface and layer \\(i=snl+1\\) is the top layer, where the variable \\(snl\\) is the negative of the number of snow layers. The number of snow layers and the thickness of each layer is a function of snow depth \\(z\_{sno}\\) (m) as follows.

\\\[\\begin{split}\\left\\{ \\begin{array}{l} snl=-1 \\\\ \\Delta z\_{0} = z\_{sno} \\end{array} \\right\\} & \\qquad {\\rm for\\; 0.01}\\le {\\rm z}\_{{\\rm sno}} \\le 0.03 \\\\\\end{split}\\\]

\\\[\\begin{split}\\left\\{ \\begin{array}{l} snl=-2 \\\\ \\Delta z\_{-1} ={z\_{sno} \\mathord{\\left/ {\\vphantom {z\_{sno} 2}} \\right.} 2} \\\\ \\Delta z\_{0} = \\Delta z\_{-1} \\end{array} \\right\\} & \\qquad {\\rm for\\; 0.03}\\, {\\rm <}\\, {\\rm z}\_{{\\rm sno}} \\le 0.04 \\\\\\end{split}\\\]

\\\[\\begin{split}\\left\\{ \\begin{array}{l} snl=-2 \\\\ \\Delta z\_{-1} = 0.02 \\\\ \\Delta z\_{0} = z\_{sno} -\\Delta z\_{-1} \\end{array} \\right\\} & \\qquad {\\rm for\\; 0.04}\\, {\\rm <}\\, {\\rm z}\_{{\\rm sno}} \\le 0.07 \\\\\\end{split}\\\]

\\\[\\begin{split}\\left\\{ \\begin{array}{l} snl=-3 \\\\ \\Delta z\_{-2} = 0.02 \\\\ \\Delta z\_{-1} = {\\left(z\_{sno} -0.02\\right)\\mathord{\\left/ {\\vphantom {\\left(z\_{sno} -0.02\\right) 2}} \\right.} 2} \\\\ \\Delta z\_{0} = \\Delta z\_{-1} \\end{array} \\right\\} & \\qquad {\\rm for\\; 0.07}\\, {\\rm <}\\, {\\rm z}\_{{\\rm sno}} \\le 0.12 \\\\\\end{split}\\\]

\\\[\\begin{split}\\left\\{ \\begin{array}{l} snl=-3 \\\\ \\Delta z\_{-2} = 0.02 \\\\ \\Delta z\_{-1} = 0.05 \\\\ \\Delta z\_{0} = z\_{sno} -\\Delta z\_{-2} -\\Delta z\_{-1} \\end{array} \\right\\} & \\qquad {\\rm for\\; 0.12}\\, {\\rm <}\\, {\\rm z}\_{{\\rm sno}} \\le 0.18 \\\\\\end{split}\\\]

\\\[\\begin{split}\\left\\{ \\begin{array}{l} snl=-4 \\\\ \\Delta z\_{-3} = 0.02 \\\\ \\Delta z\_{-2} = 0.05 \\\\ \\Delta z\_{-1} = {\\left(z\_{sno} -\\Delta z\_{-3} -\\Delta z\_{-2} \\right)\\mathord{\\left/ {\\vphantom {\\left(z\_{sno} -\\Delta z\_{-3} -\\Delta z\_{-2} \\right) 2}} \\right.} 2} \\\\ \\Delta z\_{0} =\\Delta z\_{-1} \\end{array} \\right\\} & \\qquad {\\rm for\\; 0.18}\\, {\\rm <}\\, {\\rm z}\_{{\\rm sno}} \\le 0.29 \\\\\\end{split}\\\]

\\\[\\begin{split}\\left\\{ \\begin{array}{l} snl=-4 \\\\ \\Delta z\_{-3} = 0.02 \\\\ \\Delta z\_{-2} = 0.05 \\\\ \\Delta z\_{-1} = 0.11 \\\\ \\Delta z\_{0} = z\_{sno} -\\Delta z\_{-3} -\\Delta z\_{-2} -\\Delta z\_{-1} \\end{array} \\right\\} & \\qquad {\\rm for\\; 0.29}\\, {\\rm <}\\, {\\rm z}\_{{\\rm sno}} \\le 0.41 \\\\\\end{split}\\\]

\\\[\\begin{split}\\left\\{ \\begin{array}{l} snl=-5 \\\\ \\Delta z\_{-4} = 0.02 \\\\ \\Delta z\_{-3} = 0.05 \\\\ \\Delta z\_{-2} = 0.11 \\\\ \\Delta z\_{-1} = {\\left(z\_{sno} -\\Delta z\_{-4} -\\Delta z\_{-3} -\\Delta z\_{-2} \\right)\\mathord{\\left/ {\\vphantom {\\left(z\_{sno} -\\Delta z\_{-4} -\\Delta z\_{-3} -\\Delta z\_{-2} \\right) 2}} \\right.} 2} \\\\ \\Delta z\_{0} = \\Delta z\_{-1} \\end{array} \\right\\} & \\qquad {\\rm for\\; 0.41}\\, {\\rm <}\\, {\\rm z}\_{{\\rm sno}} \\le 0.64 \\\\\\end{split}\\\]

\\\[\\begin{split}\\left\\{ \\begin{array}{l} snl=-5 \\\\ \\Delta z\_{-4} = 0.02 \\\\ \\Delta z\_{-3} = 0.05 \\\\ \\Delta z\_{-2} = 0.11 \\\\ \\Delta z\_{-1} = 0.23 \\\\ \\Delta z\_{0} = z\_{sno} -\\Delta z\_{-4} -\\Delta z\_{-3} -\\Delta z\_{-2} -\\Delta z\_{-1} \\end{array} \\right\\} & \\qquad {\\rm for\\; 0.64}\\, {\\rm <}\\, {\\rm z}\_{{\\rm sno}}\\end{split}\\\]

The node depths, which are located at the midpoint of the snow layers, and the layer interfaces are both referenced from the soil surface and are defined as negative values

(2.6.5)[¶](#equation-6-8 "Permalink to this equation")\\\[z\_{i} =z\_{h,\\, i} -0.5\\Delta z\_{i} \\qquad i=snl+1,\\ldots ,0\\\]

(2.6.6)[¶](#equation-6-9 "Permalink to this equation")\\\[z\_{h,\\, i} =z\_{h,\\, i+1} -\\Delta z\_{i+1} \\qquad i=snl,\\ldots ,-1.\\\]

Note that \\(z\_{h,\\, 0}\\), the interface between the bottom snow layer and the top soil layer, is zero. Thermal properties (i.e., temperature \\(T\_{i}\\) \[K\]; thermal conductivity \\(\\lambda \_{i}\\) \[W m\-1 K\-1\]; volumetric heat capacity \\(c\_{i}\\) \[J m\-3 K\-1\]) are defined for soil layers at the node depths ([Figure 2.6.1](#figure-soil-temperature-schematic)) and for snow layers at the layer midpoints. When present, snow occupies a fraction of a grid cell’s area, therefore snow depth represents the thickness of the snowpack averaged over only the snow covered area. The grid cell average snow depth is related to the depth of the snow covered area as \\(\\bar{z}\_{sno} =f\_{sno} z\_{sno}\\). By default, the grid cell average snow depth is written to the history file.

The heat flux \\(F\_{i}\\) (W m\-2) from layer \\(i\\) to layer \\(i+1\\) is

(2.6.7)[¶](#equation-6-10 "Permalink to this equation")\\\[F\_{i} =-\\lambda \\left\[z\_{h,\\, i} \\right\]\\left(\\frac{T\_{i} -T\_{i+1} }{z\_{i+1} -z\_{i} } \\right)\\\]

where the thermal conductivity at the interface \\(\\lambda \\left\[z\_{h,\\, i} \\right\]\\) is

(2.6.8)[¶](#equation-6-11 "Permalink to this equation")\\\[\\begin{split}\\lambda \\left\[z\_{h,\\, i} \\right\]=\\left\\{\\begin{array}{l} {\\frac{\\lambda \_{i} \\lambda \_{i+1} \\left(z\_{i+1} -z\_{i} \\right)}{\\lambda \_{i} \\left(z\_{i+1} -z\_{h,\\, i} \\right)+\\lambda \_{i+1} \\left(z\_{h,\\, i} -z\_{i} \\right)} \\qquad i=snl+1,\\ldots ,N\_{levgrnd} -1} \\\\ {0\\qquad i=N\_{levgrnd} } \\end{array}\\right\\}.\\end{split}\\\]

These equations are derived, with reference to [Figure 2.6.1](#figure-soil-temperature-schematic), assuming that the heat flux from \\(i\\) (depth \\(z\_{i}\\) ) to the interface between \\(i\\) and \\(i+1\\) (depth \\(z\_{h,\\, i}\\) ) equals the heat flux from the interface to \\(i+1\\) (depth \\(z\_{i+1}\\) ), i.e.,

(2.6.9)[¶](#equation-6-12 "Permalink to this equation")\\\[-\\lambda \_{i} \\frac{T\_{i} -T\_{m} }{z\_{h,\\, i} -z\_{i} } =-\\lambda \_{i+1} \\frac{T\_{m} -T\_{i+1} }{z\_{i+1} -z\_{h,\\, i} }\\\]

where \\(T\_{m}\\) is the temperature at the interface of layers \\(i\\) and \\(i+1\\).

Shown are three soil layers, \\(i-1\\), \\(i\\), and \\(i+1\\). The thermal conductivity \\(\\lambda\\), specific heat capacity \\(c\\), and temperature \\(T\\) are defined at the layer node depth \\(z\\). \\(T\_{m}\\) is the interface temperature. The thermal conductivity \\(\\lambda \\left\[z\_{h} \\right\]\\) is defined at the interface of two layers \\(z\_{h}\\). The layer thickness is \\(\\Delta z\\). The heat fluxes \\(F\_{i-1}\\) and \\(F\_{i}\\) are defined as positive upwards.

![Image 1: ../../_images/image17.png](https://escomp.github.io/ctsm-docs/versions/master/html/_images/image17.png)

Figure 2.6.1 Schematic diagram of numerical scheme used to solve for soil temperature.[¶](#id6 "Permalink to this image")

The energy balance for the \\(i^{th}\\) layer is

(2.6.10)[¶](#equation-6-13 "Permalink to this equation")\\\[\\frac{c\_{i} \\Delta z\_{i} }{\\Delta t} \\left(T\_{i}^{n+1} -T\_{i}^{n} \\right)=-F\_{i-1} +F\_{i}\\\]

where the superscripts \\(n\\) and \\(n+1\\) indicate values at the beginning and end of the time step, respectively, and \\(\\Delta t\\) is the time step (s). This equation is solved using the Crank-Nicholson method, which combines the explicit method with fluxes evaluated at \\(n\\) (\\(F\_{i-1}^{n},F\_{i}^{n}\\) ) and the implicit method with fluxes evaluated at \\(n+1\\) (\\(F\_{i-1}^{n+1},F\_{i}^{n+1}\\) )

(2.6.11)[¶](#equation-6-14 "Permalink to this equation")\\\[\\frac{c\_{i} \\Delta z\_{i} }{\\Delta t} \\left(T\_{i}^{n+1} -T\_{i}^{n} \\right)=\\alpha \\left(-F\_{i-1}^{n} +F\_{i}^{n} \\right)+\\left(1-\\alpha \\right)\\left(-F\_{i-1}^{n+1} +F\_{i}^{n+1} \\right)\\\]

where \\(\\alpha =0.5\\), resulting in a tridiagonal system of equations

(2.6.12)[¶](#equation-6-15 "Permalink to this equation")\\\[r\_{i} =a\_{i} T\_{i-1}^{n+1} +b\_{i} T\_{i}^{n+1} +c\_{i} T\_{i+1}^{n+1}\\\]

where \\(a\_{i}\\), \\(b\_{i}\\), and \\(c\_{i}\\) are the subdiagonal, diagonal, and superdiagonal elements in the tridiagonal matrix and \\(r\_{i}\\) is a column vector of constants. When surface water is present, the equation for the top soil layer has an additional term representing the surface water temperature; this results in a four element band-diagonal system of equations.

For the top soil layer \\(i=1\\), top snow layer \\(i=snl+1\\), or surface water layer, the heat flux from the overlying atmosphere \\(h\\) (W m\-2, defined as positive into the surface) is

(2.6.13)[¶](#equation-6-16 "Permalink to this equation")\\\[h^{n+1} =-\\alpha F\_{i-1}^{n} -\\left(1-\\alpha \\right)F\_{i-1}^{n+1} .\\\]

The energy balance for these layers is then

(2.6.14)[¶](#equation-6-17 "Permalink to this equation")\\\[\\frac{c\_{i} \\Delta z\_{i} }{\\Delta t} \\left(T\_{i}^{n+1} -T\_{i}^{n} \\right)=h^{n+1} +\\alpha F\_{i}^{n} +\\left(1-\\alpha \\right)F\_{i}^{n+1} .\\\]

The heat flux \\(h\\) at \\(n+1\\) may be approximated as follows

(2.6.15)[¶](#equation-6-18 "Permalink to this equation")\\\[h^{n+1} =h^{n} +\\frac{\\partial h}{\\partial T\_{i} } \\left(T\_{i}^{n+1} -T\_{i}^{n} \\right).\\\]

The resulting equations are then

(2.6.16)[¶](#equation-6-19 "Permalink to this equation")\\\[\\begin{split}\\begin{array}{rcl} {\\frac{c\_{i} \\Delta z\_{i} }{\\Delta t} \\left(T\_{i}^{n+1} -T\_{i}^{n} \\right)} & {=} & {h^{n} +\\frac{\\partial h}{\\partial T\_{i} } \\left(T\_{i}^{n+1} -T\_{i} \\right)} \\\\ {} & {} & {-\\alpha \\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]\\left(T\_{i}^{n} -T\_{i+1}^{n} \\right)}{z\_{i+1} -z\_{i} } -\\left(1-\\alpha \\right)\\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]\\left(T\_{i}^{n+1} -T\_{i+1}^{n+1} \\right)}{z\_{i+1} -z\_{i} } } \\end{array}\\end{split}\\\]

For the top snow layer, \\(i=snl+1\\), the coefficients are

(2.6.17)[¶](#equation-6-20 "Permalink to this equation")\\\[a\_{i} =0\\\]

(2.6.18)[¶](#equation-6-21 "Permalink to this equation")\\\[b\_{i} =1+\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left\[\\left(1-\\alpha \\right)\\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]}{z\_{i+1} -z\_{i} } -\\frac{\\partial h}{\\partial T\_{i} } \\right\]\\\]

(2.6.19)[¶](#equation-6-22 "Permalink to this equation")\\\[c\_{i} =-\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]}{z\_{i+1} -z\_{i} }\\\]

(2.6.20)[¶](#equation-6-23 "Permalink to this equation")\\\[r\_{i} =T\_{i}^{n} +\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left\[h\_{sno} ^{n} -\\frac{\\partial h}{\\partial T\_{i} } T\_{i}^{n} +\\alpha F\_{i} \\right\]\\\]

where

(2.6.21)[¶](#equation-6-24 "Permalink to this equation")\\\[F\_{i} =-\\lambda \\left\[z\_{h,\\, i} \\right\]\\left(\\frac{T\_{i}^{n} -T\_{i+1}^{n} }{z\_{i+1} -z\_{i} } \\right).\\\]

The heat flux into the snow surface from the overlying atmosphere \\(h\\) is

(2.6.22)[¶](#equation-6-25 "Permalink to this equation")\\\[h=\\overrightarrow{S}\_{sno} -\\overrightarrow{L}\_{sno} -H\_{sno} -\\lambda E\_{sno}\\\]

where \\(\\overrightarrow{S}\_{sno}\\) is the solar radiation absorbed by the top snow layer (section [2.3.2.1](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Surface_Albedos/CLM50_Tech_Note_Surface_Albedos.html#snow-albedo)), \\(\\overrightarrow{L}\_{sno}\\) is the longwave radiation absorbed by the snow (positive toward the atmosphere) (section [2.4.2](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Radiative_Fluxes/CLM50_Tech_Note_Radiative_Fluxes.html#longwave-fluxes)), \\(H\_{sno}\\) is the sensible heat flux from the snow (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)), and \\(\\lambda E\_{sno}\\) is the latent heat flux from the snow (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)). The partial derivative of the heat flux \\(h\\) with respect to temperature is

(2.6.23)[¶](#equation-6-26 "Permalink to this equation")\\\[\\frac{\\partial h}{\\partial T\_{} } =-\\frac{\\partial \\overrightarrow{L}\_{} }{\\partial T\_{} } -\\frac{\\partial H\_{} }{\\partial T\_{} } -\\frac{\\partial \\lambda E\_{} }{\\partial T\_{} }\\\]

where the partial derivative of the net longwave radiation is

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

and the partial derivatives of the sensible and latent heat fluxes are given by equations and for non-vegetated surfaces, and by equations and for vegetated surfaces. \\(\\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 \\(\\varepsilon \_{g}\\) is the ground emissivity (section [2.4.2](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Radiative_Fluxes/CLM50_Tech_Note_Radiative_Fluxes.html#longwave-fluxes)). For purposes of computing \\(h\\) and \\(\\frac{\\partial h}{\\partial T\_{g} }\\), the term \\(\\lambda\\) is arbitrarily assumed to be

(2.6.25)[¶](#equation-6-28 "Permalink to this equation")\\\[\\begin{split}\\lambda =\\left\\{\\begin{array}{l} {\\lambda \_{sub} \\qquad {\\rm if\\; }w\_{liq,\\, snl+1} =0{\\rm \\; and\\; }w\_{ice,\\, snl+1} >0} \\\\ {\\lambda \_{vap} \\qquad {\\rm otherwise}} \\end{array}\\right\\}\\end{split}\\\]

where \\(\\lambda \_{sub}\\) and \\(\\lambda \_{vap}\\) are the latent heat of sublimation and vaporization, respectively (J kg\-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)), and \\(w\_{liq,\\, snl+1}\\) and \\(w\_{ice,\\, snl+1}\\) are the liquid water and ice contents of the top snow/soil layer, respectively (kg m\-2) (Chapter [2.7](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Hydrology/CLM50_Tech_Note_Hydrology.html#rst-hydrology)).

For the top soil layer, \\(i=1\\), the coefficients are

(2.6.26)[¶](#equation-6-29 "Permalink to this equation")\\\[a\_{i} =-f\_{sno} \\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} }\\\]

(2.6.27)[¶](#equation-6-30 "Permalink to this equation")\\\[b\_{i} =1+\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left\[f\_{sno} \\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} } +\\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]}{z\_{i+1} -z\_{i} } \\right\]-\\left(1-f\_{sno} \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\frac{\\partial h}{\\partial T}\\\]

(2.6.28)[¶](#equation-6-31 "Permalink to this equation")\\\[c\_{i} =-\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]}{z\_{i+1} -z\_{i} }\\\]

(2.6.29)[¶](#equation-6-32 "Permalink to this equation")\\\[r\_{i} =T\_{i}^{n} +\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left\[\\left(1-f\_{sno} \\right)\\left(h\_{soil} ^{n} -\\frac{\\partial h}{\\partial T\_{} } T\_{i}^{n} \\right)+\\alpha \\left(F\_{i} -f\_{sno} F\_{i-1} \\right)\\right\]\\\]

The heat flux into the soil surface from the overlying atmosphere \\(h\\) is

(2.6.30)[¶](#equation-6-33 "Permalink to this equation")\\\[h=\\overrightarrow{S}\_{soil} -\\overrightarrow{L}\_{soil} -H\_{soil} -\\lambda E\_{soil}\\\]

It can be seen that when no snow is present (\\(f\_{sno} =0\\)), the expressions for the coefficients of the top soil layer have the same form as those for the top snow layer.

The surface snow/soil layer temperature computed in this way is the layer-averaged temperature and hence has somewhat reduced diurnal amplitude compared with surface temperature. An accurate surface temperature is provided that compensates for this effect and numerical error by tuning the heat capacity of the top layer (through adjustment of the layer thickness) to give an exact match to the analytic solution for diurnal heating. The top layer thickness for \\(i=snl+1\\) is given by

(2.6.31)[¶](#equation-6-34 "Permalink to this equation")\\\[\\Delta z\_{i\*} =0.5\\left\[z\_{i} -z\_{h,\\, i-1} +c\_{a} \\left(z\_{i+1} -z\_{h,\\, i-1} \\right)\\right\]\\\]

where \\(c\_{a}\\) is a tunable parameter, varying from 0 to 1, and is taken as 0.34 by comparing the numerical solution with the analytic solution ([Z.-L. Yang 1998, unpublished manuscript](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/References/CLM50_Tech_Note_References.html#yang1998)). \\(\\Delta z\_{i\*}\\) is used in place of \\(\\Delta z\_{i}\\) for \\(i=snl+1\\) in equations -. The top snow/soil layer temperature computed in this way is the ground surface temperature \\(T\_{g}^{n+1}\\).

The boundary condition at the bottom of the snow/soil column is zero heat flux, \\(F\_{i} =0\\), resulting in, for \\(i=N\_{levgrnd}\\),

(2.6.32)[¶](#equation-6-35 "Permalink to this equation")\\\[\\frac{c\_{i} \\Delta z\_{i} }{\\Delta t} \\left(T\_{i}^{n+1} -T\_{i}^{n} \\right)=\\alpha \\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]\\left(T\_{i-1}^{n} -T\_{i}^{n} \\right)}{z\_{i} -z\_{i-1} } +\\left(1-\\alpha \\right)\\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]\\left(T\_{i-1}^{n+1} -T\_{i}^{n+1} \\right)}{z\_{i} -z\_{i-1} }\\\]

(2.6.33)[¶](#equation-6-36 "Permalink to this equation")\\\[a\_{i} =-\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} }\\\]

(2.6.34)[¶](#equation-6-37 "Permalink to this equation")\\\[b\_{i} =1+\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} }\\\]

(2.6.35)[¶](#equation-6-38 "Permalink to this equation")\\\[c\_{i} =0\\\]

(2.6.36)[¶](#equation-6-39 "Permalink to this equation")\\\[r\_{i} =T\_{i}^{n} -\\alpha \\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } F\_{i-1}\\\]

where

(2.6.37)[¶](#equation-6-40 "Permalink to this equation")\\\[F\_{i-1} =-\\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} } \\left(T\_{i-1}^{n} -T\_{i}^{n} \\right).\\\]

For the interior snow/soil layers, \\(snl+1<i<N\_{levgrnd}\\), excluding the top soil layer,

(2.6.38)[¶](#equation-6-41 "Permalink to this equation")\\\[\\begin{split}\\begin{array}{rcl} {\\frac{c\_{i} \\Delta z\_{i} }{\\Delta t} \\left(T\_{i}^{n+1} -T\_{i}^{n} \\right)} & {=} & {-\\alpha \\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]\\left(T\_{i}^{n} -T\_{i+1}^{n} \\right)}{z\_{i+1} -z\_{i} } +\\alpha \\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]\\left(T\_{i-1}^{n} -T\_{i}^{n} \\right)}{z\_{i} -z\_{i-1} } } \\\\ {} \\end{array}\\end{split}\\\]

(2.6.39)[¶](#equation-6-42 "Permalink to this equation")\\\[a\_{i} =-\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} }\\\]

(2.6.40)[¶](#equation-6-43 "Permalink to this equation")\\\[b\_{i} =1+\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left\[\\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} } +\\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]}{z\_{i+1} -z\_{i} } \\right\]\\\]

(2.6.41)[¶](#equation-6-44 "Permalink to this equation")\\\[c\_{i} =-\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]}{z\_{i+1} -z\_{i} }\\\]

(2.6.42)[¶](#equation-6-45 "Permalink to this equation")\\\[r\_{i} =T\_{i}^{n} +\\alpha \\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left(F\_{i} -F\_{i-1} \\right)+\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\vec{S}\_{g,i} .\\\]

where \\(\\vec{S}\_{g,i}\\) is the absorbed solar flux in layer \\(i\\) (section [2.3.2.1](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Surface_Albedos/CLM50_Tech_Note_Surface_Albedos.html#snow-albedo)).

When surface water exists, the following top soil layer coefficients are modified

(2.6.43)[¶](#equation-6-46 "Permalink to this equation")\\\[\\begin{split}\\begin{array}{l} {b\_{i} =1+\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left\[f\_{h2osfc} \\frac{\\lambda \_{h2osfc} }{z\_{i} -z\_{h2osfc} } +f\_{sno} \\frac{\\lambda \\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} } +\\frac{\\lambda \\left\[z\_{h,\\, i} \\right\]}{z\_{i+1} -z\_{i} } \\right\]} \\\\ {\\quad \\quad -\\left(1-f\_{sno} -f\_{h2osfc} \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\frac{\\partial h}{\\partial T} } \\end{array}\\end{split}\\\]

(2.6.44)[¶](#equation-6-47 "Permalink to this equation")\\\[\\begin{split}r\_{i} =T\_{i}^{n} +\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left\[\\begin{array}{l} {\\left(1-f\_{sno} -f\_{h2osfc} \\right)\\left(h\_{soil} ^{n} -\\frac{\\partial h}{\\partial T\_{} } T\_{i}^{n} \\right)} \\\\ {+\\alpha \\left(F\_{i} -f\_{sno} F\_{i-1} +f\_{h2osfc} \\frac{\\lambda \_{h2osfc} }{z\_{1} -z\_{h2osfc} } \\left(T\_{1} -T\_{h2osfc} \\right)\\right)} \\end{array}\\right\]\\end{split}\\\]

(2.6.45)[¶](#equation-6-48 "Permalink to this equation")\\\[d\_{i} =-f\_{h2osfc} \\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left\[\\frac{\\lambda \_{h2osfc} }{z\_{i} -z\_{h2osfc} } \\right\]\\\]

where \\(d\_{i}\\) is an additional coefficient representing the heat flux from the surface water layer. The surface water layer coefficients are

(2.6.46)[¶](#equation-6-49 "Permalink to this equation")\\\[a\_{h2osfc} =0\\\]

(2.6.47)[¶](#equation-6-50 "Permalink to this equation")\\\[b\_{h2osfc} =1+\\frac{\\Delta t}{c\_{h2osfc} \\Delta z\_{h2osfc} } \\left\[\\left(1-\\alpha \\right)\\frac{\\lambda \_{h2osfc} }{z\_{1} -z\_{h2osfc} } -\\frac{\\partial h}{\\partial T} \\right\]\\\]

(2.6.48)[¶](#equation-6-51 "Permalink to this equation")\\\[c\_{h2osfc} =-\\left(1-\\alpha \\right)\\frac{\\Delta t}{c\_{h2osfc} \\Delta z\_{h2osfc} } \\frac{\\lambda \_{h2osfc} }{z\_{1} -z\_{h2osfc} }\\\]

(2.6.49)[¶](#equation-6-52 "Permalink to this equation")\\\[r\_{h2osfc} =T\_{h2osfc}^{n} +\\frac{\\Delta t}{c\_{i} \\Delta z\_{i} } \\left\[h\_{h2osfc} ^{n} -\\frac{\\partial h}{\\partial T\_{} } T\_{h2osfc}^{n} +\\alpha \\frac{\\lambda \_{h2osfc} }{z\_{1} -z\_{h2osfc} } \\left(T\_{1} -T\_{h2osfc} \\right)\\right\]\_{}\\\]

