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

With reference to [Figure 2.7.2](#figure-water-flux-schematic), the equation for conservation of mass (equation [(2.7.41)](#equation-7-79)) can be integrated over each layer as

(2.7.64)[¶](#equation-7-101 "Permalink to this equation")\\\[\\int \_{-z\_{h,\\, i} }^{-z\_{h,\\, i-1} }\\frac{\\partial \\theta }{\\partial t} \\, dz=-\\int \_{-z\_{h,\\, i} }^{-z\_{h,\\, i-1} }\\frac{\\partial q}{\\partial z} \\, dz-\\int \_{-z\_{h,\\, i} }^{-z\_{h,\\, i-1} } e\\, dz .\\\]

Note that the integration limits are negative since \\(z\\) is defined as positive upward from the soil surface. This equation can be written as

(2.7.65)[¶](#equation-7-102 "Permalink to this equation")\\\[\\Delta z\_{i} \\frac{\\partial \\theta\_{liq,\\, i} }{\\partial t} =-q\_{i-1} +q\_{i} -e\_{i}\\\]

where \\(q\_{i}\\) is the flux of water across interface \\(z\_{h,\\, i}\\), \\(q\_{i-1}\\) is the flux of water across interface \\(z\_{h,\\, i-1}\\), and \\(e\_{i}\\) is a layer-averaged soil moisture sink term (ET loss) defined as positive for flow out of the layer (mm s\-1). Taking the finite difference with time and evaluating the fluxes implicitly at time \\(n+1\\) yields

(2.7.66)[¶](#equation-7-103 "Permalink to this equation")\\\[\\frac{\\Delta z\_{i} \\Delta \\theta\_{liq,\\, i} }{\\Delta t} =-q\_{i-1}^{n+1} +q\_{i}^{n+1} -e\_{i}\\\]

where \\(\\Delta \\theta\_{liq,\\, i} =\\theta\_{liq,\\, i}^{n+1} -\\theta\_{liq,\\, i}^{n}\\) is the change in volumetric soil liquid water of layer \\(i\\) in time \\(\\Delta t\\)and \\(\\Delta z\_{i}\\) is the thickness of layer \\(i\\) (mm).

The water removed by transpiration in each layer \\(e\_{i}\\) is a function of the total transpiration \\(E\_{v}^{t}\\) (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 the effective root fraction \\(r\_{e,\\, i}\\)

(2.7.67)[¶](#equation-7-104 "Permalink to this equation")\\\[e\_{i} =r\_{e,\\, i} E\_{v}^{t} .\\\]

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

Figure 2.7.2 Schematic diagram of numerical scheme used to solve for soil water fluxes.[¶](#id12 "Permalink to this image")

Shown are three soil layers, \\(i-1\\), \\(i\\), and \\(i+1\\). The soil matric potential \\(\\psi\\) and volumetric soil water \\(\\theta\_{liq}\\) are defined at the layer node depth \\(z\\). The hydraulic conductivity \\(k\\left\[z\_{h} \\right\]\\) is defined at the interface of two layers \\(z\_{h}\\). The layer thickness is \\(\\Delta z\\). The soil water fluxes \\(q\_{i-1}\\) and \\(q\_{i}\\) are defined as positive upwards. The soil moisture sink term \\(e\\) (ET loss) is defined as positive for flow out of the layer.

Note that because more than one plant functional type (PFT) may share a soil column, the transpiration \\(E\_{v}^{t}\\) is a weighted sum of transpiration from all PFTs whose weighting depends on PFT area as

(2.7.68)[¶](#equation-7-105 "Permalink to this equation")\\\[E\_{v}^{t} =\\sum \_{j=1}^{npft}\\left(E\_{v}^{t} \\right)\_{j} \\left(wt\\right)\_{j}\\\]

where \\(npft\\) is the number of PFTs sharing a soil column, \\(\\left(E\_{v}^{t} \\right)\_{j}\\) is the transpiration from the \\(j^{th}\\) PFT on the column, and \\(\\left(wt\\right)\_{j}\\) is the relative area of the \\(j^{th}\\) PFT with respect to the column. The effective root fraction \\(r\_{e,\\, i}\\) is also a column-level quantity that is a weighted sum over all PFTs. The weighting depends on the per unit area transpiration of each PFT and its relative area as

(2.7.69)[¶](#equation-7-106 "Permalink to this equation")\\\[r\_{e,\\, i} =\\frac{\\sum \_{j=1}^{npft}\\left(r\_{e,\\, i} \\right)\_{j} \\left(E\_{v}^{t} \\right)\_{j} \\left(wt\\right)\_{j} }{\\sum \_{j=1}^{npft}\\left(E\_{v}^{t} \\right)\_{j} \\left(wt\\right)\_{j} }\\\]

where \\(\\left(r\_{e,\\, i} \\right)\_{j}\\) is the effective root fraction for the \\(j^{th}\\) PFT

(2.7.70)[¶](#equation-7-107 "Permalink to this equation")\\\[\\begin{split}\\begin{array}{lr} \\left(r\_{e,\\, i} \\right)\_{j} =\\frac{\\left(r\_{i} \\right)\_{j} \\left(w\_{i} \\right)\_{j} }{\\left(\\beta \_{t} \\right)\_{j} } & \\qquad \\left(\\beta \_{t} \\right)\_{j} >0 \\\\ \\left(r\_{e,\\, i} \\right)\_{j} =0 & \\qquad \\left(\\beta \_{t} \\right)\_{j} =0 \\end{array}\\end{split}\\\]

and \\(\\left(r\_{i} \\right)\_{j}\\) is the fraction of roots in layer \\(i\\) (Chapter [2.9](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Photosynthesis/CLM50_Tech_Note_Photosynthesis.html#rst-stomatal-resistance-and-photosynthesis)), \\(\\left(w\_{i} \\right)\_{j}\\) is a soil dryness or plant wilting factor for layer \\(i\\) (Chapter [2.9](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Photosynthesis/CLM50_Tech_Note_Photosynthesis.html#rst-stomatal-resistance-and-photosynthesis)), and \\(\\left(\\beta\_{t} \\right)\_{j}\\) is a wetness factor for the total soil column for the \\(j^{th}\\) PFT (Chapter [2.9](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/Photosynthesis/CLM50_Tech_Note_Photosynthesis.html#rst-stomatal-resistance-and-photosynthesis)).

The soil water fluxes in [(2.7.66)](#equation-7-103),, which are a function of \\(\\theta\_{liq,\\, i}\\) and \\(\\theta\_{liq,\\, i+1}\\) because of their dependence on hydraulic conductivity and soil matric potential, can be linearized about \\(\\theta\\) using a Taylor series expansion as

(2.7.71)[¶](#equation-7-108 "Permalink to this equation")\\\[q\_{i}^{n+1} =q\_{i}^{n} +\\frac{\\partial q\_{i} }{\\partial \\theta\_{liq,\\, i} } \\Delta \\theta\_{liq,\\, i} +\\frac{\\partial q\_{i} }{\\partial \\theta\_{liq,\\, i+1} } \\Delta \\theta\_{liq,\\, i+1}\\\]

(2.7.72)[¶](#equation-7-109 "Permalink to this equation")\\\[q\_{i-1}^{n+1} =q\_{i-1}^{n} +\\frac{\\partial q\_{i-1} }{\\partial \\theta\_{liq,\\, i-1} } \\Delta \\theta\_{liq,\\, i-1} +\\frac{\\partial q\_{i-1} }{\\partial \\theta\_{liq,\\, i} } \\Delta \\theta\_{liq,\\, i} .\\\]

Substitution of these expressions for \\(q\_{i}^{n+1}\\) and \\(q\_{i-1}^{n+1}\\) into [(2.7.66)](#equation-7-103) results in a general tridiagonal equation set of the form

(2.7.73)[¶](#equation-7-110 "Permalink to this equation")\\\[r\_{i} =a\_{i} \\Delta \\theta\_{liq,\\, i-1} +b\_{i} \\Delta \\theta\_{liq,\\, i} +c\_{i} \\Delta \\theta\_{liq,\\, i+1}\\\]

where

(2.7.74)[¶](#equation-7-111 "Permalink to this equation")\\\[a\_{i} =-\\frac{\\partial q\_{i-1} }{\\partial \\theta\_{liq,\\, i-1} }\\\]

(2.7.75)[¶](#equation-7-112 "Permalink to this equation")\\\[b\_{i} =\\frac{\\partial q\_{i} }{\\partial \\theta\_{liq,\\, i} } -\\frac{\\partial q\_{i-1} }{\\partial \\theta\_{liq,\\, i} } -\\frac{\\Delta z\_{i} }{\\Delta t}\\\]

(2.7.76)[¶](#equation-7-113 "Permalink to this equation")\\\[c\_{i} =\\frac{\\partial q\_{i} }{\\partial \\theta\_{liq,\\, i+1} }\\\]

(2.7.77)[¶](#equation-7-114 "Permalink to this equation")\\\[r\_{i} =q\_{i-1}^{n} -q\_{i}^{n} +e\_{i} .\\\]

The tridiagonal equation set is solved over \\(i=1,\\ldots,N\_{levsoi}\\).

The finite-difference forms of the fluxes and partial derivatives in equations [(2.7.74)](#equation-7-111) - [(2.7.77)](#equation-7-114) can be obtained from equation [(2.7.44)](#equation-7-82) as

(2.7.78)[¶](#equation-7-115 "Permalink to this equation")\\\[q\_{i-1}^{n} =-k\\left\[z\_{h,\\, i-1} \\right\]\\left\[\\frac{\\left(\\psi \_{i-1} -\\psi \_{i} \\right)+\\left(z\_{i} - z\_{i-1} \\right)}{z\_{i} -z\_{i-1} } \\right\]\\\]

(2.7.79)[¶](#equation-7-116 "Permalink to this equation")\\\[q\_{i}^{n} =-k\\left\[z\_{h,\\, i} \\right\]\\left\[\\frac{\\left(\\psi \_{i} -\\psi \_{i+1} \\right)+\\left(z\_{i+1} - z\_{i} \\right)}{z\_{i+1} -z\_{i} } \\right\]\\\]

(2.7.80)[¶](#equation-7-117 "Permalink to this equation")\\\[\\frac{\\partial q\_{i-1} }{\\partial \\theta \_{liq,\\, i-1} } =-\\left\[\\frac{k\\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} } \\frac{\\partial \\psi \_{i-1} }{\\partial \\theta \_{liq,\\, i-1} } \\right\]-\\frac{\\partial k\\left\[z\_{h,\\, i-1} \\right\]}{\\partial \\theta \_{liq,\\, i-1} } \\left\[\\frac{\\left(\\psi \_{i-1} -\\psi \_{i} \\right)+\\left(z\_{i} - z\_{i-1} \\right)}{z\_{i} - z\_{i-1} } \\right\]\\\]

(2.7.81)[¶](#equation-7-118 "Permalink to this equation")\\\[\\frac{\\partial q\_{i-1} }{\\partial \\theta \_{liq,\\, i} } =\\left\[\\frac{k\\left\[z\_{h,\\, i-1} \\right\]}{z\_{i} -z\_{i-1} } \\frac{\\partial \\psi \_{i} }{\\partial \\theta \_{liq,\\, i} } \\right\]-\\frac{\\partial k\\left\[z\_{h,\\, i-1} \\right\]}{\\partial \\theta \_{liq,\\, i} } \\left\[\\frac{\\left(\\psi \_{i-1} -\\psi \_{i} \\right)+\\left(z\_{i} - z\_{i-1} \\right)}{z\_{i} - z\_{i-1} } \\right\]\\\]

(2.7.82)[¶](#equation-7-119 "Permalink to this equation")\\\[\\frac{\\partial q\_{i} }{\\partial \\theta \_{liq,\\, i} } =-\\left\[\\frac{k\\left\[z\_{h,\\, i} \\right\]}{z\_{i+1} -z\_{i} } \\frac{\\partial \\psi \_{i} }{\\partial \\theta \_{liq,\\, i} } \\right\]-\\frac{\\partial k\\left\[z\_{h,\\, i} \\right\]}{\\partial \\theta \_{liq,\\, i} } \\left\[\\frac{\\left(\\psi \_{i} -\\psi \_{i+1} \\right)+\\left(z\_{i+1} - z\_{i} \\right)}{z\_{i+1} - z\_{i} } \\right\]\\\]

(2.7.83)[¶](#equation-7-120 "Permalink to this equation")\\\[\\frac{\\partial q\_{i} }{\\partial \\theta \_{liq,\\, i+1} } =\\left\[\\frac{k\\left\[z\_{h,\\, i} \\right\]}{z\_{i+1} -z\_{i} } \\frac{\\partial \\psi \_{i+1} }{\\partial \\theta \_{liq,\\, i+1} } \\right\]-\\frac{\\partial k\\left\[z\_{h,\\, i} \\right\]}{\\partial \\theta \_{liq,\\, i+1} } \\left\[\\frac{\\left(\\psi \_{i} -\\psi \_{i+1} \\right)+\\left(z\_{i+1} - z\_{i} \\right)}{z\_{i+1} - z\_{i} } \\right\].\\\]

The derivatives of the soil matric potential at the node depth are derived from [(2.7.53)](#equation-7-94)

(2.7.84)[¶](#equation-7-121 "Permalink to this equation")\\\[\\frac{\\partial \\psi \_{i-1} }{\\partial \\theta\_{liq,\\, \\, i-1} } =-B\_{i-1} \\frac{\\psi \_{i-1} }{\\theta\_{\\, \\, i-1} }\\\]

(2.7.85)[¶](#equation-7-122 "Permalink to this equation")\\\[\\frac{\\partial \\psi \_{i} }{\\partial \\theta\_{\\, liq,\\, i} } =-B\_{i} \\frac{\\psi \_{i} }{\\theta\_{i} }\\\]

(2.7.86)[¶](#equation-7-123 "Permalink to this equation")\\\[\\frac{\\partial \\psi \_{i+1} }{\\partial \\theta\_{liq,\\, i+1} } =-B\_{i+1} \\frac{\\psi \_{i+1} }{\\theta\_{\\, i+1} }\\\]

with the constraint \\(0.01\\, \\theta\_{sat,\\, i} \\le \\theta\_{\\, i} \\le \\theta\_{sat,\\, i}\\).

The derivatives of the hydraulic conductivity at the layer interface are derived from [(2.7.47)](#equation-7-85)

(2.7.87)[¶](#equation-7-124 "Permalink to this equation")\\\[\\begin{array}{l} {\\frac{\\partial k\\left\[z\_{h,\\, i-1} \\right\]}{\\partial \\theta \_{liq,\\, i-1} } = \\frac{\\partial k\\left\[z\_{h,\\, i-1} \\right\]}{\\partial \\theta \_{liq,\\, i} } = \\left(2B\_{i-1} +3\\right) \\ \\overline{\\Theta}\_{ice} \\ k\_{sat} \\left\[z\_{h,\\, i-1} \\right\] \\ \\left\[\\frac{\\overline{\\theta}\_{liq}}{\\overline{\\theta}\_{sat}} \\right\]^{2B\_{i-1} +2} \\left(\\frac{0.5}{\\overline{\\theta}\_{sat}} \\right)} \\end{array}\\\]

where \\(\\overline{\\Theta}\_{ice} = \\Theta(\\overline{\\theta}\_{ice})\\) [(2.7.48)](#equation-7-86), \\(\\overline{\\theta}\_{ice} = 0.5\\left(\\theta\_{ice\\, i-1} +\\theta\_{ice\\, i} \\right)\\), \\(\\overline{\\theta}\_{liq} = 0.5\\left(\\theta\_{liq\\, i-1} +\\theta\_{liq\\, i} \\right)\\), and \\(\\overline{\\theta}\_{sat} = 0.5\\left(\\theta\_{sat,\\, i-1} +\\theta\_{sat,\\, i} \\right)\\)

and

(2.7.88)[¶](#equation-7-125 "Permalink to this equation")\\\[\\begin{array}{l} {\\frac{\\partial k\\left\[z\_{h,\\, i} \\right\]}{\\partial \\theta \_{liq,\\, i} } = \\frac{\\partial k\\left\[z\_{h,\\, i} \\right\]}{\\partial \\theta \_{liq,\\, i+1} } = \\left(2B\_{i} +3\\right) \\ \\overline{\\Theta}\_{ice} \\ k\_{sat} \\left\[z\_{h,\\, i} \\right\] \\ \\left\[\\frac{\\overline{\\theta}\_{liq}}{\\overline{\\theta}\_{sat}} \\right\]^{2B\_{i} +2} \\left(\\frac{0.5}{\\overline{\\theta}\_{sat}} \\right)} \\end{array}.\\\]

where \\(\\overline{\\theta}\_{liq} = 0.5\\left(\\theta\_{\\, i} +\\theta\_{\\, i+1} \\right)\\), \\(\\overline{\\theta}\_{sat} = 0.5\\left(\\theta\_{sat,\\, i} +\\theta\_{sat,\\, i+1} \\right)\\).

#### 2.7.3.2.1. Equation set for layer \\(i=1\\)[¶](#equation-set-for-layer-i-1 "Permalink to this headline")

For the top soil layer (\\(i=1\\)), the boundary condition is the infiltration rate (section [2.7.2.1](#surface-runoff)), \\(q\_{i-1}^{n+1} =-q\_{infl}^{n+1}\\), and the water balance equation is

(2.7.89)[¶](#equation-7-135 "Permalink to this equation")\\\[\\frac{\\Delta z\_{i} \\Delta \\theta\_{liq,\\, i} }{\\Delta t} =q\_{infl}^{n+1} +q\_{i}^{n+1} -e\_{i} .\\\]

After grouping like terms, the coefficients of the tridiagonal set of equations for \\(i=1\\) are

(2.7.90)[¶](#equation-7-136 "Permalink to this equation")\\\[a\_{i} =0\\\]

(2.7.91)[¶](#equation-7-137 "Permalink to this equation")\\\[b\_{i} =\\frac{\\partial q\_{i} }{\\partial \\theta\_{liq,\\, i} } -\\frac{\\Delta z\_{i} }{\\Delta t}\\\]

(2.7.92)[¶](#equation-7-138 "Permalink to this equation")\\\[c\_{i} =\\frac{\\partial q\_{i} }{\\partial \\theta\_{liq,\\, i+1} }\\\]

(2.7.93)[¶](#equation-7-139 "Permalink to this equation")\\\[r\_{i} =q\_{infl}^{n+1} -q\_{i}^{n} +e\_{i} .\\\]

#### 2.7.3.2.2. Equation set for layers \\(i=2,\\ldots ,N\_{levsoi} -1\\)[¶](#equation-set-for-layers-i-2-ldots-n-levsoi-1 "Permalink to this headline")

The coefficients of the tridiagonal set of equations for \\(i=2,\\ldots,N\_{levsoi} -1\\) are

(2.7.94)[¶](#equation-7-140 "Permalink to this equation")\\\[a\_{i} =-\\frac{\\partial q\_{i-1} }{\\partial \\theta\_{liq,\\, i-1} }\\\]

(2.7.95)[¶](#equation-7-141 "Permalink to this equation")\\\[b\_{i} =\\frac{\\partial q\_{i} }{\\partial \\theta\_{liq,\\, i} } -\\frac{\\partial q\_{i-1} }{\\partial \\theta\_{liq,\\, i} } -\\frac{\\Delta z\_{i} }{\\Delta t}\\\]

(2.7.96)[¶](#equation-7-142 "Permalink to this equation")\\\[c\_{i} =\\frac{\\partial q\_{i} }{\\partial \\theta\_{liq,\\, i+1} }\\\]

(2.7.97)[¶](#equation-7-143 "Permalink to this equation")\\\[r\_{i} =q\_{i-1}^{n} -q\_{i}^{n} +e\_{i} .\\\]

#### 2.7.3.2.3. Equation set for layer \\(i=N\_{levsoi}\\)[¶](#equation-set-for-layer-i-n-levsoi "Permalink to this headline")

For the lowest soil layer (\\(i=N\_{levsoi}\\) ), a zero-flux bottom boundary condition is applied (\\(q\_{i}^{n} =0\\)) and the coefficients of the tridiagonal set of equations for \\(i=N\_{levsoi}\\) are

(2.7.98)[¶](#equation-7-148 "Permalink to this equation")\\\[a\_{i} =-\\frac{\\partial q\_{i-1} }{\\partial \\theta\_{liq,\\, i-1} }\\\]

(2.7.99)[¶](#equation-7-149 "Permalink to this equation")\\\[b\_{i} =\\frac{\\partial q\_{i} }{\\partial \\theta\_{liq,\\, i} } -\\frac{\\partial q\_{i-1} }{\\partial \\theta\_{liq,\\, i} } -\\frac{\\Delta z\_{i} }{\\Delta t}\\\]

(2.7.100)[¶](#equation-7-150 "Permalink to this equation")\\\[c\_{i} =0\\\]

(2.7.101)[¶](#equation-7-151 "Permalink to this equation")\\\[r\_{i} =q\_{i-1}^{n} +e\_{i} .\\\]

#### 2.7.3.2.4. Adaptive Time Stepping[¶](#adaptive-time-stepping "Permalink to this headline")

The length of the time step is adjusted in order to improve the accuracy and stability of the numerical solutions. The difference between two numerical approximations is used to estimate the temporal truncation error, and then the step size \\(\\Delta t\_{sub}\\) is adjusted to meet a user-prescribed error tolerance [\[Kavetski et al., 2002\]](https://escomp.github.io/ctsm-docs/versions/master/html/tech_note/References/CLM50_Tech_Note_References.html#kavetskietal2002). The temporal truncation error is estimated by comparing the flux obtained from the first-order Taylor series expansion (\\(q\_{i-1}^{n+1}\\) and \\(q\_{i}^{n+1}\\), equations [(2.7.71)](#equation-7-108) and [(2.7.72)](#equation-7-109)) against the flux at the start of the time step (\\(q\_{i-1}^{n}\\) and \\(q\_{i}^{n}\\)). Since the tridiagonal solution already provides an estimate of \\(\\Delta \\theta\_{liq,i}\\), it is convenient to compute the error for each of the \\(i\\) layers from equation [(2.7.66)](#equation-7-103) as

(2.7.102)[¶](#equation-7-152 "Permalink to this equation")\\\[\\epsilon\_{i} = \\left\[ \\frac{\\Delta \\theta\_{liq,\\, i} \\Delta z\_{i}}{\\Delta t\_{sub}} - \\left( q\_{i-1}^{n} - q\_{i}^{n} + e\_{i}\\right) \\right\] \\ \\frac{\\Delta t\_{sub}}{2}\\\]

and the maximum absolute error across all layers as

(2.7.103)[¶](#equation-7-153 "Permalink to this equation")\\\[\\begin{array}{lr} \\epsilon\_{crit} = {\\rm max} \\left( \\left| \\epsilon\_{i} \\right| \\right) & \\qquad 1 \\le i \\le nlevsoi \\end{array} \\ .\\\]

The adaptive step size selection is based on specified upper and lower error tolerances, \\(\\tau\_{U}\\) and \\(\\tau\_{L}\\). The solution is accepted if \\(\\epsilon\_{crit} \\le \\tau\_{U}\\) and the procedure repeats until the adaptive sub-stepping spans the full model time step (the sub-steps are doubled if \\(\\epsilon\_{crit} \\le \\tau\_{L}\\), i.e., if the solution is very accurate). Conversely, the solution is rejected if \\(\\epsilon\_{crit} > \\tau\_{U}\\). In this case the length of the sub-steps is halved and a new solution is obtained. The halving of substeps continues until either \\(\\epsilon\_{crit} \\le \\tau\_{U}\\) or the specified minimum time step length is reached.

Upon solution of the tridiagonal equation set, the liquid water contents are updated as follows

(2.7.104)[¶](#equation-7-164 "Permalink to this equation")\\\[w\_{liq,\\, i}^{n+1} =w\_{liq,\\, i}^{n} +\\Delta \\theta\_{liq,\\, i} \\Delta z\_{i} \\qquad i=1,\\ldots ,N\_{levsoi} .\\\]

The volumetric water content is

(2.7.105)[¶](#equation-7-165 "Permalink to this equation")\\\[\\theta\_{i} =\\frac{w\_{liq,\\, i} }{\\Delta z\_{i} \\rho \_{liq} } +\\frac{w\_{ice,\\, i} }{\\Delta z\_{i} \\rho \_{ice} } .\\\]

