specific heat of water, cp,w (4217.7 J/kgK) or ice, cp,i ( 13.3 + 7.80Ta J/kgK) depending
on Tp, Up the mass precipitation flux (kg/m2s), and γ p the precipitation density (kg/m3). It
is assumed that the precipitation temperature is the same as the air temperature.
4.1.6 Final Top Boundary Equation
Combining Equations (4.2), (4.4), (4.5) (4.8) and (4.11), the top boundary condition is
reconfigured as
F (T ) = (1  α ) I s ↓ +ε Iι↓  εσ T 4 + (e0 + ρac p,aCDW )(Ta  T )
h
ρ
∂T
 ρaCDW ′l (qg  qa ) + U pc pTp + κ
e
@ z = 0m.
(4.12)
∂z
ρi ∂θi
+l fus
∆z  vc pT = 0
ρw ∂t
This is solved at each time step for the surface temperature.
4.2 Numerical Solution
Equation (4.1) is solved using a modified secondorder CrankNicholson approach.
Following the technique presented in Hornbeck (1975), Equation (4.1) is rewritten as
Tj+1,i  Tj,i ⎡ kthj ,i+1  kthj ,i ⎤ ⎡ Tj,i+1  Tj,i ⎤
=⎢
⎥+
⎥⎢
∆t
∆z
∆z
⎢
⎥⎣
⎦
⎣
⎦
kthj ,i ⎡ Tj+1,i+1  2Tj+1,i + Tj+1,i1 Tj,i+1  2T j,i + Tj,i1 ⎤
+
⎢
⎥
(4.13)
( ∆z )
( ∆z )
2
2
2 ⎢
⎥
⎣
⎦
⎡ Tj,i+1  Tj,i ⎤ l fus ρi ( ∆θi ) j,i
c p,w

⎥+
.
v j,i ⎢
⎦ c p ρw ∆t
∆z
c p j ,i
⎣
where the subscripts j and i represent time and depth, respectively. Combining like terms
and rearranging so that all terms involving Tj+1 are on the lefthand side of the equation,
Equation (4.13) becomes
T j+1,i1 + γ j,iT j+1,i + T j+1,i+1 =
(4.14)
T j,i1 + T j,i ( β j,i + α j,i  δ j,i ) + T j,i+1(1  β j,i + δ j,i ) + j,i ≡ φ j,i
where
33