通用的土壤水热传输耦合模型的发展和改进研究

中国科学 D辑: 地球科学

《中国科学》杂志社 2007年 第37卷 第11期: 1522~1535

http://www.scichina.com SCIENCE IN CHINA PRESS

通用的土壤水热传输耦合模型的发展和改进研究

李 倩

①②

孙菽芬*

①

(① 中国科学院大气物理研究所大气科学和地球流体力学数值模拟国家重点实验室(LASG), 北京 100029;

② 中国科学院研究生院, 北京 100047)

摘要 一个既真实又简化、且适用于湿润与干旱、冻土与非冻土和均质与非均质多种情景下的通用土壤水热传输耦合模式对于陆面过程的模型发展研究至关重要. 研究首先通过量级估计和模型数值模拟结果的分析, 发展了简化且精度较好的土壤水热传输耦合统一模型. 为了克服该模型计算过程中由于需预估冰水相变速率项产生的误差造成的不确定性, 进一步对该简化的统一土壤模式进行变量变换, 以土壤总焓和土壤水总质量替代温度和体积含水量作为方程预报量, 建立了新的通用土壤模型统一体系, 并设计了一套行之有效、省时的数值计算方案. 此模型既可用于一般情况下的裸土, 也可用于较为难处理的非均质土、冻融土壤和干旱土壤等. 与观测结果相比, 改进后的统一土壤模型能很好地模拟出在湿润与干旱、冻土与非冻土和均质与非均质土壤中的水热传输过程. 且由于它的简化, 也适应当今陆面过程模式发展的需要. 关键词 通用、简化的统一土壤模型 水热传输耦合 冻土 干旱 非均质 大量的用陆面模式和GCM耦合所进行的敏感性实验[1~4]表明: 发展一个能真实刻画各种下垫面上的陆面过程模式至关重要. 土壤是陆面过程模型要研究的最基本、最重要的陆地下垫面, 它分布广泛、质地繁多且内部物理过程复杂, 在年际、年代际或更长的时间尺度上对各圈层的活动产生很大的影响, 其本身也可经历湿润与干旱、冻结与融化不断交替的变化. 在陆面过程模型的发展过程中, 这部份的模拟一直是重要的关注点, 已得到很大的提高和改善. 但由于认识所限, 目前常用的陆面过程模式对土壤下垫面的考虑仍有以下不足:

(1) 大部分模型所考虑的土壤状态主要为湿润的土壤, 土壤参数化方案使用的是忽略了水-热传输耦合作用的等温土壤模式[5,6], 仅分别考虑液态水流动对水分分布和热传导过程对温度分布影响, 忽略了水汽运动对土壤水-热平衡可能的影响, 不能很好

用于干旱沙漠和半干旱土壤. 对于冻土冻融相变过程, 也未有很好的参数化方案, 尽管一些有影响的模型都考虑了冻融过程的事实, 但定量描述该物理过程的参数化方案比较简单, 偏离真实的物理关系: 如SSiB和CLM模式[6,7]采用了冻(融)过程和相变热交换仅发生在0℃时(冰点)的简化假设; BATS模式对冻融过程的考虑相对复杂一些, 假设冻融过程发生在冰点附近某一温度(0~−4℃)范围内[8]. 而真实的冻融过程应是一个连续缓慢的变化过程, 不存在固定的冻融的临界温度(或范围), 必须根据平衡态情况下土壤水势和温度之间的热力学平衡关系以及固有的(但目前只能有实验确定的)土壤水力学特征本构关系来确定土壤含水量、含冰量和温度之间的定量关系.

(2) 大部分模型[9~11]基本上采用含水量梯度替代水势梯度的表达式描述水流运动, 使描述水分平衡的方程仅能用于均质土壤, 难于用于非均质的土壤、

收稿日期: 2007-05-11; 接受日期: 2007-06-25

国家自然科学基金(批准号: 40575043, 40233034)和中国科学院知识创新工程(批准号: KZCX3-SW-229)资助项目 * 联系人, E-mail: ssf@lasg.iap.ac.cn

第11期

李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1523

冻土物理过程的研究和与饱和地下水运动方程耦合的研究.

(3) 几乎现有常用的土壤模型中, 能量和质量平衡方程的描述基本上都以土壤温度、体积含水量作为预报量, 使与含冰量变化率有关的项显性地出现在控制方程中, 在计算迭代求解时对该变化率预估的误差会对冻土温度的计算产生较大偏差, 从而影响计算迭代收敛速度和模拟结果.

因此, 要想很好地模拟复杂的地气交换过程, 发展一个通用于冻土和非冻土, 均质和非均质土壤、干旱、半干旱和湿润状态下其内部水热输送过程的统一土壤模式无疑是非常重要的, 而且为了可用于陆面模型和气候模型的发展研究所需, 模型应该力求既简化而又不失真.

纵观已经发展的各种模型, 在处理土壤中水热传输过程的复杂程度上均有不同, 可应用的范围也不同[9~25]. 它们又大致可分为三类: 第一类, 完全的土壤水热传输耦合模型

[15]

耦合模型.

1 非饱和统一土壤水热传输耦合模型

1.1 完整复杂的土壤水热传输耦合模型

Zhao等[15]提出的完整、复杂的土壤水热传输耦合模型可作为本研究的出发点. 它共有10个方程, 其中4个预报方程:

体积含冰量θ i的控制方程:

󰀅+M󰀅∂θiMi,vi,l+=0, (1)

ρi∂t体积含水量θ1的控制方程:

󰀅−M󰀅∂θl∂ulMl,vi,l

++=0, (2)

ρl∂t∂Z水蒸汽质量体积含量 θvρv的控制方程:

∂θvρv

󰀅−M󰀅=∂⎡D∂ρv⎤, (3) −Mi,vl,veff

∂t∂Z⎢∂Z⎥⎣⎦

(10个方程), 完整考虑了水-

描述温度T*(℃)(=T(°K)−Tf(=273.15 K))变化的能量

冰-气三相的变化及三相之间的相变过程、热传导和水汽扩散相变对土壤能量和质量平衡的影响; 第二类, 相对简化的土壤水热传输耦合模型[16,17](5个方程), 忽略了水汽局地变化的影响; 第三类, 当前常用的简化的土壤水热传输耦合模型. 其水质量平衡方程未考虑水汽扩散对质量平衡的作用[9,10,18~24]. 能量方程中也不考虑水汽扩散相变和液态水运动携带的热量[9,22,23]. 此外, 还有一种介于第二类和第三类之间的土壤模型, 它们考虑了液态水运动携带的热量对土壤能量平衡的影响[18~21,24,25].

针对以上的讨论, 显而易见以下几个问题亟需解决: 什么形式的土壤水热输送耦合方程能既简单又不失真地描述在广泛的土壤温度和湿度变化范围内、不同质地的土壤内部其水热传输过程, 并能较好地避免土壤模式计算过程中由于冰-水相变速率的预估误差产生的计算不确定性？不同质地土壤冻融相变过程中, 含水量、含冰量和温度之间合理的相变热力学关系应该如何确定？本文将以完整复杂的土壤水热传输耦合模型为出发点, 考虑到非均质土、湿润土壤、干旱、半干旱土壤及冻土的特殊性, 对以上这几个问题进行仔细的分析研究, 对模型进行分层次的合理简化和数学变换, 得到一个真实合理且适用于陆面模型和气候模型研究的统一通用的土壤水热

控制方程:

∂CvT*∂ρθ∂ρθ−Li,lii+Ll,vvv∂t∂t∂t=−ρlcl

∂ulT∂⎛

+Ll,v⎜Deff∂Z∂Z⎝

*

∂ρv⎞∂⎛∂T*⎞

⎜Keff⎟⎟+

∂Z⎠∂Z⎜∂Z⎟⎝⎠

(4)

和5个诊断方程: 由推广的Darcy定律描述液态水流动速率ul的方程:

⎡∂ψ⎤

＋1⎥, (5) ul=Kl⎢−⎣∂Z⎦

冰点水势方程:

ψ=

Li,l(T−273.15)

gTf

, (6)

ψ为土壤基膜水势(这里不考虑土壤溶质势). 决定ψ与体积含冰量θi和体积含水量θl之间的推广的Clapp&Hornberger经验关系为[10]

⎛θ⎞

ψ=ψ0⎜l⎟

⎝θs⎠

−b

(1+ckθi)2, (7)

此外还有, 水汽密度ρv与当地水蒸汽压ev的关系; 体积水蒸汽含量θv与土壤空隙度θs、体积含冰量θi和体积含水量θl的关系[15].

在以上方程中, Z和t表示土壤的深度和时间; T

1524

中国科学 D辑 地球科学

第37卷

󰀅,M󰀅和 M󰀅分别表示冰到蒸汽、(K)为温度; M液i,vl,vi,l

∂CvT*∂ρiθi∂ulT*

−Li,l=−ρlcl ∂t∂t∂Z

∂⎛∂T*∂ψ⎞∂⎛∂T*⎞

++LlvDψV+K⎜Ll,vDTV⎟⎟∂Z⎜⎜eff∂Z⎟⎟, (10) ∂Z⎜∂∂ZZ⎝⎠⎝⎠

态水到蒸汽和冰到液态水的相变速率; ρl和ρi分别为液态水和冰的固有密度; c1为液态水的热容量; Cv为平均的体积热容量; Li,v, Ll,v和Li,l分别为冰到蒸汽、液态水到蒸汽和冰到液态水的相变潜热; Kl为水文导度; Keff为有效热传导系数; Deff为水蒸汽在有效空间θv中的有效扩散系数; g为重力加速度; ψ0为饱和水势; b为Clapp和Hornberger常数[26]; ck为一常数. 以上共有9个方程, 10个未知数(θi,θl,θv,ρv, T, 󰀅, M󰀅,M󰀅), 需补充一个条件以封闭方ul,ψ,Mi,vl,vi,l

其中DTV为温度梯度引起的水汽扩散系数, DψV为水势梯度引起的水汽扩散系数. 可以看出, 上述方程组已对复杂完整的方程组作了合理初步的简化, 并去掉了不必要的假定关系, 且仍是通用于均质和非均质土, 冻土和非冻土, 以及不同湿度状态的水热耦合模型.

以下将从已初步简化的耦合模型方程(9)和(10)出发, 在不影响精度的前提下, 考虑到不同土壤特性、不同温度变化及不同湿度变化的情况, 对其进行深入的量级分析, 尽可能的简化, 以得到一个普遍适用的土壤水热传输耦合模型.

首先, 根据量级分析的方法, 对控制方程(9), (10)中水热通量的各项进行比较, 再根据由该模式计算的水热通量垂直分布大小, 来分析方程中各项的大小关系.

可看出, 在控制方程中, 与方程(9)有关的水通量贡献项分别为

QL=−Kl

∂ψ+Kl, ∂Z

∂ψ, (11) ∂Z

程组. 为此, Zhao等[15]曾提出一个相变假设关系来封闭方程组, 而Jordan[27]则提出另一组不同的假设关系, 它们的正确性都没有得到证明, 从而降低了模式的准确性. 尽管以上方程体系完备通用, 它可通用于冻土和非冻土, 干旱、半干旱和湿润的状态. 但很显然, 这是十分复杂的方程组, 而且还包含一个正确性并没有得到证明的假设, 并不适用于当今陆面模型的发展.

1.2 通用、统一土壤水热传输耦合方程的合理简化

可以证明, 根据土壤温度剖面[28,29]变化的规律, 对水汽运动方程中水汽局地项和水汽扩散项的大小进行量级分析, 并利用数值模拟实验的结果, 发现在⎛∂θρ水汽方程(3)左端的水汽局地变化项⎜vv

⎝∂t⎛∂

右端的水汽扩散项⎜

⎝∂Z

⎞

⎟远小于⎠

QMV=−DψVQTV=−DTV

∂ρv⎤⎞⎡Deff⎢⎥⎟. 因此可忽略水∂Z⎣⎦⎠

∂T*

; ∂Z∂T*

, ∂Z

⎞

⎟, (12) ⎠

与方程(10)有关的热通量贡献项分别为

QHconduct=−Keff

汽的局地变化项对土壤水分的贡献1), 故水汽方程(3)退化为

󰀅−M󰀅Mi,vl,v

∂

=∂Z

∂ρv⎤⎡D⎢eff∂Z⎥, (8) ⎣⎦

⎛∂ψQHconvect=ρlclT*Kl⎜1−

⎝∂ZQHvapor

利用该方程, 上述方程组中四个预报方程(1)~(4)可初步简化为描述土壤中水分变化的质量方程和描述土壤中温度分布的能量方程:

∂θl∂θi∂ul∂⎛∂T*∂ψ⎞ρlDD,+ρi=−ρl++⎜TV⎟ψV⎟ ZZ∂t∂t∂Z∂Z⎜∂∂⎝⎠

(9)

1) Li Q, Sun S F. Development of a simplified frozen soil model

⎛∂ψ∂T*⎞=−⎜Ll,vDψV+Ll,vDTV,⎟⎜⎟ ∂∂ZZ⎝⎠

其中, QL, QMV和QTV分别为液态水流通量、由土 壤水势梯度和温度梯度引起的水汽扩散通量; QHconduct, QHconvect和QHvapor分别为热传导通量、 液态水流动携带的热通量和水汽扩散相变引起的热

第11期

李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1525

通量.

本文选取了6种不同的土壤: 砂土(sand), 砂壤土(sandy loam), 粉壤土(silt loam), 壤土(loam), 砂黏土(sandy clay), 黏土(clay), 并利用Clapp&Hornberger 提出的经验水势和水文导度关系(所用的各项水热物理参数见表1).

在不同的温度, 不同的含水量和含冰量情况下, 计算每种土壤的K1, DTV/ρl, DψV/ρl, Keff, Ll,vDψV, Ll,vDTV, ρlclKlT*及∂ψ/∂θl. 作为典型示

略. 由于在水热通量各项中还与温度梯度∂T*/∂Z和 水势梯度项∂ψ/∂Z有关, 而土壤表层附近的∂T*/∂Z和∂ψ/∂Z是比较显著和重要的, 所以可以用表层附近的温度和水势梯度的量级作为∂T*/∂Z和∂ψ/∂Z具有代表意义的值. 根据分析[28,29], 土壤中日温度波的传播决定了土壤表层温度随深度的变化, 而日温度波的振幅随深度呈指数衰减(exp(−z/d), d是日温度波衰减深度, 一般为O(10 cm) ). 那么在10 cm厚的土壤表层, 如果地表的日温度波振幅大约为∆T≈O(15~30℃), 则在深度d处的日温度波振幅约

范, 仅给出了黏土(表2)和砂土(表3)的结果, 其余从

a) 表1 不同土壤类型的物理参数

类型

砂土/%

黏土/%b) θ sc) ψSd) KSe) bf) 92 3 0.339 1.07E-06 2.79 砂土 −0.069 58 10 0.434 −0.141 5.23E-06 4.74 砂壤土

17 13 0.476 −0.759 2.81E-06 5.33 粉壤土

43 18 0.439 −0.355 3.38E-06 5.25 壤土

52 42 0.406 0.098 7.22E-06 10.73 砂黏土

22 58 0.468 −0.468 9.74E-07 11.55 黏土

a) 砂土比例; b) 黏土比例; c) 土壤孔隙度; d) 土壤饱和水势; e) 饱和导水度; f) Clapp&Hornberger常数

表2

黏土在非冻结时, 不同温度下(>0℃)各物理量大小

/m·m T=0.0℃

0.05

0.2 0.46 T=10.0℃

0.05

0.2 0.46 T=30.0℃

0.05 0.2 0.46

3

θ1

−3

K1 /m·s−1 4.30E-32 2.20E-16 6.20E-07 4.30E-32 2.20E-16 6.20E-07

4.30E-32 2.20E-16 6.20E-07

DTV/ρ1 /m2·(K·s)−1

0.00E+00 7.20E-13 4.00E-14 0.00E+00 1.40E-12 7.70E-14

0.00E+00 4.50E-12 2.40E-13

DψV/ρ1 /m·s−1 0.00E+00 7.80E-16 4.60E-17

0.00E+00 1.60E-15 9.10E-17

0.00E+00 5.70E-15 3.10E-16

Keff /J·(K·m·s)−1

2.40E-01 1.20E+00 1.70E+00 2.40E-01 1.20E+00 1.70E+00

2.40E-01 1.20E+00 1.70E+00

L1VDψV /J·m−2·s−1

0.00E+00 2.00E-06 1.10E-07

0.00E+00 4.00E-06 2.30E-07

0.00E+00 1.40E-05 7.70E-07

L1,VDTV /J·(m·K·s)−1

0.00E+00 1.80E-03 1.00E-04

0.00E+00 3.50E-03 1.90E-04

0.00E+00 1.10E-02 6.00E-04

ρ1c1K1T*

−2

/J·m·s

0.00E+00 0.00E+00 0.00E+00 1.80E-24 9.40E-09 2.60E+01

5.50E-24 2.80E-08 7.80E+01

−1

∂ψ ∂Z

1.80E+13 5.00E+05 1.40E+01

1.80E+13 5.00E+05 1.40E+01

1.80E+13 5.00E+05 1.40E+01

黏土在冻结时, 不同温度下(<0℃)各物理量大小

/m·m T=−2.044℃

0.015 0 2.60E-18 1.40E-12 1.50E-15 5.80E-01 3.90E-06 3.50E-03 −2.20E-11 0.027 0.1 1.30E-17 1.00E-12 0.00E+00 2.10E+00 0.00E+00 2.60E-03 −1.20E-10 0.039 0.25 1.80E-18 2.40E-13 0.00E+00 4.40E+00 0.00E+00 6.10E-04 −1.50E-11 T=−9.98℃

0.009 0 3.20E-20 7.30E-13 7.80E-16 5.80E-01 1.90E-06 1.80E-03 −1.30E-12 0.011 0.01 9.90E-20 8.20E-13 0.00E+00 5.50E-01 0.00E+00 2.10E-03 −4.20E-12 0.02 0.2 3.10E-20 3.10E-13 0.00E+00 3.50E+00 0.00E+00 7.70E-04 −1.30E-12 T=−20.423℃

0.007 0 3.70E-21 2.80E-13 2.90E-16 5.80E-01 7.10E-07 7.10E-04 −3.20E-13 0.015 0.185 4.30E-21 1.40E-13 0.00E+00 3.30E+00 0.00E+00 3.60E-04 −3.70E-13 0.017 0.25 1.50E-21 7.40E-14 0.00E+00 4.30E+00 0.00E+00 1.90E-04 −1.30E-13

7.70E+04 2.60E+04 1.80E+04 5.30E+05 3.20E+05 1.70E+05 1.40E+06 4.70E+05 4.10E+05

3

θ1

−3

/m·m

3

θi

−3

K1 /m·s−1 DTV/ρ1 /m·(K·s)−1

2DψV/ρ1 /m·s−1 Keff /J·(K·m·s)−1L1VDψV /J·m−2·s−1 L1,VDTV /J·(m·K·s)−1

/J·m−2·s−1

ρ1c1K1T*

∂ψ ∂Z

1526

中国科学 D辑 地球科学

第37卷

表3

砂土在非冻结时, 不同温度下(>0 ℃)各物理量的大小 /m·m T=0.0℃

0.05 0.2 0.33 T=10.0℃

0.05 0.2 0.33 T=30.0℃

0.05 0.2 0.33

3

θ1

−3

/m·m

0 0 0 0 0 0 0 0 0

3

θi

−3

K1 /m·s−1 7.90E-14 1.20E-08 8.50E-07

7.90E-14 1.20E-08 8.50E-07

7.90E-14 1.20E-08 8.50E-07

DTV/ρ1 /m·(K·s)−1

2DψV/ρ1 /m·s−1 1.60E-15 7.90E-16 5.10E-17

3.30E-15 1.60E-15 1.00E-16

1.10E-14 5.40E-15 3.50E-16

Keff /J·(K·m·s)−1

8.30E-01 2.70E+00 3.40E+00

8.30E-01 2.70E+00 3.40E+00

8.30E-01 2.70E+00 3.40E+00

L1VDψV /J·m−2·s−1

4.10E-06 2.00E-06 1.30E-07

8.20E-06 4.00E-06 2.60E-07

2.80E-05 1.30E-05 8.70E-07

L1,VDTV /J·(m·K·s)−1

3.70E-03 1.80E-03 1.10E-04

6.90E-03 3.30E-03 2.20E-04

2.20E-02 1.00E-02 6.70E-04

/J·m−2·s−1

0.00E+00 0.00E+00 0.00E+00

3.30E-06 4.80E-01 3.60E+01

9.90E-06 1.50E+00 1.10E+02

ρ1c1K1T*

∂ψ ∂Z

8.00E+02 4.20E+00 6.30E-01

8.00E+02 4.20E+00 6.30E-01

8.00E+02 4.20E+00 6.30E-01

1.50E-12 7.00E-13 4.60E-14

2.80E-12 1.30E-12 8.60E-14

8.60E-12 4.10E-12 2.70E-13

砂土在冻结时, 不同温度下(<0 ℃)各物理量大小

/m·m T=−2.044℃

0.015 0 2.60E-18 1.40E-12 1.50E-15 5.80E-01 3.90E-06 3.50E-03 −2.20E-11 0.025 0.075 1.50E-17 1.20E-12 0.00E+00 1.70E+00 0.00E+00 2.90E-03 −1.30E-10 0.039 0.25 1.80E-18 2.40E-13 0.00E+00 4.40E+00 0.00E+00 6.10E-04 −1.50E-11 T=−9.98℃

0.009 0 3.20E-20 7.30E-13 7.80E-16 5.80E-01 1.90E-06 1.80E-03 −1.30E-12 0.011 0.01 9.90E-20 8.20E-13 0.00E+00 5.50E-01 0.00E+00 2.10E-03 −4.20E-12 0.02 0.2 3.10E-20 3.10E-13 0.00E+00 3.50E+00 0.00E+00 7.70E-04 −1.30E-12 T=−20.423℃

0.007 0 3.70E-21 2.80E-13 2.90E-16 5.80E-01 7.10E-07 7.10E-04 −3.20E-13 0.01 0.05 1.40E-20 2.90E-13 0.00E+00 1.10E+00 0.00E+00 7.20E-04 −1.20E-12 0.017 0.25 1.50E-21 7.40E-14 0.00E+00 4.30E+00 0.00E+00 1.90E-04 −1.30E-13

1.40E+06 7.10E+05 4.10E+05 5.30E+05 3.20E+05 1.70E+05 7.70E+04 2.80E+04 1.80E+04

3

θ1

−3

/m·m

3

θi

−3

K1 /m·s−1 DTV/ρ1 /m·(K·s)−1

2DψV/ρ1 /m·s−1 Keff /J·(K·m·s)−1L1VDψV /J·m−2·s−1 L1,VDTV /J·(m·K·s)−1

/J·m−2·s−1

ρ1c1K1T*

∂ψ ∂Z

为O(e−1∆T)≈5~11℃, 所以在土壤表层附近温度梯度的量级大小约为O(10℃/10 cm)=O(100℃/m). 对⎛∂ψ∂θl⎞

⋅于表层的水势梯度来说, 可根据O⎜⎟来确定, ∂∂Zθ⎝l⎠⎛∂θ一般来说O⎜l

⎝∂Z⎛∂ψO⎜⎝∂Z

⎛∂ψ⎞

⎟≈O⎜⎠⎝∂θl

⎞

⎟ ≈O(0.1/10 cm)≈O(1/m), 所以⎠⎞⎛∂ψ⎟. O⎜⎠⎝∂θl

⎞

⎟可根据各种土壤的水力⎠

从以上的量级分析及数值模拟实验结果都可以看出: 对于与热通量有关的各项来说, 当土壤非常干 或者冻结比较厉害的时候, ρlclT*Kl(∂ψ/∂Z)总是比Keff(∂T*/∂Z)小2个量级以上. 而当土壤接近饱和

时, Keff(∂T*/∂Z)与ρlclT*Kl的量级差别不会很 大. 因此为了模型的完整性ρlcl(∂ulT*/∂Z)可保 留. 且对于土壤表层来说, 液态水运动携带的热 通量ρlcl(∂ulT*/∂Z)相当于地表入渗携带的热通量. 计算表明, 当日平均降水量大于10 mm时, 入渗携带的热通量其量级与Keff(∂T*/∂Z)相当. 因此在土壤的能量平衡方程中应考虑液态水运动携带的热通量对温度分布的影响. QHvapor无论在什么情况下都很小. 因此, 在能量方程中可以忽略QHvapor的作用. 观

学特征关系式来确定.

为了进一步分析比较水热通量的各项, 我们利用较为完整的模式(9)和(10)计算了土壤在冻融阶段和非冻融阶段水热通量的变化情况. 图1和2给出了模式模拟的D66站的土壤在未冻结-冻结-融化时期几个典型时刻的水热通量垂直分布.

第11期

李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1527

图1 模拟的D66站土壤在未冻结-冻结-融化阶段不同时刻的水通量(water flux)垂直分布. QTV, QMV和QL 分别代表QTV,

QMV和QL

未冻结阶段: (a), (b); 冻结阶段: (c), (d); 融化阶段: (e), (f)

察与水通量有关的各项, 发现不同的土壤在不同的DTV∂T*∂ψ和Kl相比, 两者含水(冰)量情况下,

∂Zρl∂Z

能小于、等于或大于后者. 所以在水质量平衡方 程中QTV应该保留. 同样对于

DψV∂ψ来说, 在冻

ρl∂Z

的量级大小并不确定, 由于土壤湿度不同, 前者有可

结情况下, 由于冻土空隙中的水汽密度只是温度

1528

中国科学 D辑 地球科学

第37卷

图2 模拟的D66站土壤在未冻结-冻结-融化阶段不同时刻的热通量(heat flux)垂直分布

hflux_conv, hflux_vap 和 hflux_cond 分别代表了QHconvect, QHvapor和QHconduct. 未冻结阶段: (a), (b); 冻结阶段: (c), (d); 融化阶段: (e), (f)

的函数, 所以DψV为零. 而土壤未冻结时,

DψV∂ψρl∂Z

项都很小, 但它的相对重要性逐渐显现出来. 所 以, 在质量平衡方程中该水汽通量的作用也可保

在土壤湿度不太干情况下, 相比于液态水通量 留.

因此可以得到合理简化的统一土壤水热传输耦∂ψ 可以省略, 但在土壤变得很干时, 虽然它和Kl

合方程: ∂Z

第11期

李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1529

∂θlρ∂θ∂⎛∂ψ⎞

=−ii−+Kl⎟ ⎜−Kl∂tρl∂t∂Z⎝∂Z⎠

1∂⎛∂T*∂ψ⎞

DD++⎜TV⎟ψV⎟, (13) ρl∂Z⎜ZZ∂∂⎝⎠

∂CvT*∂ρiθi∂ulT*∂⎛∂T*⎞

−Li,l=－ρlcl. (14) ⎜Keff⎟⎟∂t∂t∂Z⎜∂∂ZZ⎝⎠

热力学平衡态基础上, 根据土壤中冰-水-汽三相之间的相平衡关系严格推导得到的冻融土壤中水势和温度之间的关系式. 在冻土中, 如果接受系统处于热力学平衡态的假设 (这是必须遵循的基本物理关系), 联立(15)和(6)就能得到非饱和土壤中关联含水量、含冰量和温度的冻融相变平衡关系:

除了以上2个预报方程外, 模型中还有3个诊断方程: (5)~(7)或(15)

θi=θi(θl,T) 或 θl=θl(θi,T), (16)

上式能定量地刻画不同质地的土壤在冻融状态下, 其含水量、含冰量和温度之间的关系, 因此也称为非饱和土壤冻融参数化方案. 在系统处于热力学平衡态的假设下, (15)和(16)式中只有一个是独立的, 另一个可根据(6)推出. 很多研究中, 在接受本构关系(15)的同时, 又根据实验直接得出一个经验关系式(16)来作为土壤冻融参数化方案, 使有严格物理基础的关系式(6)变得多余. 合理的途经是在冻土中从基本物理关系式(6)入手, 尽可能地得到一个合理且较为普遍适用的冻融参数化方案.

观察后发现, 在大多数研究中, 非饱和土壤水力学特征关系式(7)或(15)已得到广泛的应用和认 可[26,30]. 然而, Faruoki[31]指出冰的存在对土壤中水势有一定影响. 因此目前在冻土和非冻土的研究中, 大都采用包含了冰含量影响的经验水力学特征关系修正式(7). 由(7)和(6)便可得冻融相变平衡关系:

⎡Li,l⋅T*⎤

θl=θs⋅⎢(1+ckθi)−2⎥

⎢⎥⎣gψoTf⎦

−1bψ=ψ(θi,θl). (15)

这样, 5个方程(13), (14), (5), (6)和(7)或(15), 5个未知变量(θi, θl, T, ul, ψ )构成了一个简化的通用土壤水热传输耦合模型.

需要指出的是, 部分模式描述液态水流动:

∂

∂Z⎛∂θl⎜D⎝∂Z

[9~11]

*

用了含水量梯度来

⎞

⎟. 其中D是液态水的扩⎠

散系数(m2/s), 且D=Kl

∂ψ. 应该明确, 在均匀土中∂θl

以上方程是可以合理地描述液态水的变化情况, 但是在非均匀土或土壤冻结时, 由于土壤分层的交界处含水量不连续, 水流方向也许与土壤湿度降低的方向完全相反, 含水量的梯度并不能真正地反应土壤液态水的流动. 且当土壤发生冻结时, 土壤水势不仅是液态水体积含量的函数, 还是含冰量的函数, 因∂ψ⎛∂ψ∂θl⎞⎛∂ψ∂θi⎞

=⎜此推广的Darcy定律中⎟+⎜⎟, ∂z⎝∂θl∂Z⎠⎝∂θi∂Z⎠

, (17)

而不是

∂ψ⎛∂ψ∂θl⎞

=⎜⎟. 所以, 在统一土壤模式中考∂z⎝∂θl∂Z⎠

∂∂Z⎛∂ψ⎜Kl

⎝∂Z

⎞

⎟)是更合理的. ⎠

式中的b, ψo和θs是依赖于土壤特性的参数.

从式(17)可看出, (1)土壤的冻结和融化是缓慢变化的过程, 并不是在某一固定温度(如0℃)或某一温度区间内就能完成的. 且即使在温度很低的时候非饱和冻土中仍然有液态水存在. 这与很多室内外观测事实吻合. (2)在(17)式中令θi=0.0, 可得到某一温度T*(零摄氏度以下)时土壤中存在的未结冰的最大液态水含量θlmax. 即在温度T*的情况下, 只有当土壤中的液态水大于θlmax时多于的水才会冻结. 此外, 变换(17)式可得

max

gψ0Tf⎛θl*

⎜T=

Li,l⎜θs

⎝

虑到非均质土和冻土的情况, 采用土壤水势梯度来描述的液态水流动(

1.3 统一土壤水热传输耦合模型中的冻融参数化方案

在土壤模型中, 方程(6)和(7)或(15)是决定冻结土壤中体积含水量、含冰量和土壤温度之间定量关系的两个重要关系式. 其中, 方程(15)是描述冻土和非冻土中土壤水势和土壤含水量(或含冰量)之间的水力学特征关系, 也称为非饱和冻土中的本构关系,虽然目前这一关系式只能通过实验得到. 而方程(6)是在

⎞

⎟⎟⎠

−b

(1+ckθi)2, (18a)

*

: 令θi=0.0, 也能确定土壤开始冻融的临界温度Tcrit

1530

中国科学 D辑 地球科学

第37卷

*Tcrit

gψ0Tf=

Li,l⎛θT⎞⎜⎟⎝θs⎠

−b

, (18b)

式中引入了土壤表面阻抗的概念[14], 对蒸发公式进行了修正. 很多陆面过程模式(SIB, BATS)都将其引入了表层蒸发的参数化方案中, 它有助于我们理解和有效地描述土壤在干旱时表层蒸发的机制. 但是表面阻抗一般都是土壤表层含水量的经验函数, 并不是在所有情况下均适用.

研究表明[33], 在土壤比较干时, 表面会形成干表面层(Dry Surface Layer-DSL), 蒸发的水汽从DSL的下边界经过DSL被输送到土壤表面. 之后, 很多研究工作中[34,35]也利用DSL或蒸发面下移这个概念来计算干旱土壤蒸发的水汽通量, 较为合理明确地解决了干旱土壤蒸发的问题.

干旱土壤的蒸发主要有以下两个物理过程, 首先水汽通过分子扩散的形式从干表面层传输到土壤表面, 这一过程中的水汽通量可表示为

E=Deff

*

ρahqe−qs

, (19)

ρlze

此时的θl等于总的积水当量θT. 上式表明, 如果土壤

*

, 温度低于临界温度Tcrit

土壤开始结冰, 相反则开始

融化. 明显地, 土壤冻融的临界温度并不是在零摄氏度.

(3) 由于各种土壤的ψ0和θs都不相同, 因此由

*

也有差异. 图3展示了不(17)式所确定的θlmax和Tcrit

同质地的土壤在不同温度下所能含有的未结冰的最大液态水含量.

*

qe,h分别为蒸发面上土壤空隙中空气的饱和比湿和

相对湿度, qs为土壤表面上的空气比湿, ze是干表

图3 不同类型的土壤在0℃以下时, 土壤所能含有不结冰

的最大液态水含量

纵坐标表示土壤在未结冰的情况下所能含有的最大液态水含量

面层的厚度. 第二个过程, 水汽从土壤表面进入大气, 通量为

E=

可看出, 对于所有的土壤类型, 液态水含量随温度的变化在冰点附近的某一温度区间内(0~−2℃)最显著. 在−2℃时, 对于砂土, 大部分的液态水会变成冰, 因此可以规定土壤的冻融过程发生在冰点附近的某一区间内, 但对于其它土壤而言, 液态水变成冰的过程还远未完成. 对于黏土, 即使到−20℃, 土壤所能含有不结冰的最大体积液态水含量θ1仍可达0.22. 这说明假定一个统一的结冰区间可能会产生误差. 在统一的非饱和土壤水热传输耦合模型中, 采用土壤冻融参数化方案(17)可能会较为真实普遍地描述土壤的冻融过程.

ρaqs−qa

, (20)

ρlra

qa为大气的比湿, ra为空气动力学阻抗.

Kondo等[36]指出这两个过程中的水汽通量应该是相同的, 因此联立以上两式, 可以得到实际蒸发通量:

*

ρahqe−qa

, (21) E=

ρlra+ze/Deff

其中, ze/Deff等价于早期研究中[14,37]提出的表面阻抗因子. 即:

rs=ze/Deff, (22)

只是此时的rs不再是土壤表层含水量的经验函数. 所以, 在模式中, 采用以上的蒸发面下移或DSL的方法来计算表面蒸发, 避免了经验关系式, 不仅适和于不同质地的土壤, 并且其物理过程也很明确.

1.4 统一土壤水热传输耦合模型中的蒸发方案

对于干旱土来说, 由于地表附近土壤水和温度的梯度都很大, 地表面空隙中的水汽压与当地平衡饱和水汽压差别也会很大, 并不能达到局地平衡 态[32]. 所以在计算地表蒸发的时候应该放弃局地平衡的假设. 针对这一情况, 一些早期的研究在蒸发公

2 简化通用的土壤水热传输耦合统一模型的改进

可看出, 上述通用的简化耦合模型其控制方程是以土壤温度、体积含水量作为预报量的非线性方程,

第11期

李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1531

含冰量变化显性地出现在控制方程中. 当土壤发生冻融时, 能量平衡方程(14)中包含冰水相变产生的潜热项Li,l∂ρiθi∂t, 由于Li,l值很大(Li,l=3.336×105 J/kg), 故冰水相变产生的巨大潜热将强烈地影响土壤温度的改变, 而一般在求解以上的控制方程时都是采用迭代法, 在每次迭代开始需对含冰量值进行估计, 其估计的偏差将直接影响土壤温度的正确计算, 进而影响冰水相变速率自身的正确计算, 从而对模式的模拟精度和计算速度产生影响.

为了避免相变速率的不确定性引起的模式计算不稳定性, 本工作对当前通用的土壤水热传输控制方程中的预报变量进行变量变换, 用总土壤水质量

3.1.1 青藏高原季风试验D66站资料

D66站位于青藏高原北部(35°31′N, 93°47′E), 海拔4560 m, 属于大陆性高原气候, 年降水量较少. 地表植被稀疏, 土壤为非均质永冻土[38], 质地以砂壤土为主. D66自动气象站提供了观测高度在1.5 m的每30 min大气强迫场资料, 包括入射短波辐射通量、气温、气压、相对湿度和风速. 土壤温度的观测由10个白金地温探头(Pt)和数采仪获得, 地温探头埋设的深度分别为4, 20, 40, 60, 80, 100, 130, 160, 200, 263 cm. 土壤湿度(含水量)的测量由6个时域反射仪(Time-Domain Reflectometer, 简称TDR)探头和数采仪完成, TDR探头的埋设深度分别为: 4, 20, 60, 100, 160, 225 cm. 土壤的观测数据每小时自动采集记录一次.

ma(ma=ρiθi+ρlθl)和总焓H(H=CvT*−Li,lρiθi)作

为预报变量, 即将原来方程中与冰水相变有关的项均合并到总土壤水质量和土壤总焓中, 发展出新的预报方程组. 具体方程如下:

∂ma∂⎛∂ψ⎞

=−ρl+Kl⎟ ⎜−Kl∂t∂Z⎝∂Z⎠

1∂⎛∂T*∂ψ⎞

DTV,++DψV⎜⎟⎟ ZZρl∂Z⎜∂∂⎝⎠

3.1.2 敦煌双墩子戈壁微气象站观测资料

敦煌试验站位于40°10′N, 94°31′E, 海拔1150 m, 属典型的温带大陆性气候, 年平均降水量为40 mm[39]. 地表为平坦的砂石戈壁, 土壤质地是砂土. 该试验站的大气观测资料包括分别位于1, 2, 8, 18 m

铁塔的风、温、湿观测, 以及1.5 m高的架子上观测的直接辐射, 总辐射, 反射辐射, 大气向下长波辐射和地表向上长波辐射. 地表有3个传感器, 互为120度角; 土壤温度共6层, 其深度分别为5, 10, 20, 40, 80, 180 cm; 土壤湿度传感器的深度分别为5, 10, 20, 80 cm. 以上的观测每小时记录一次. 相关仪器的名称和观测精度可参见文献[40].

(23)

∂ulT*∂H∂⎛∂T*⎞

－ρlcl, (24) = ⎜Keff⎟⎟ZZ∂t∂Z⎜∂∂⎝⎠

这样, 每次对方程(23)和(24)迭代时无需对冰水相变速率进行估计就能得到土壤水总质量和土壤总焓, 然后再根据方程(23)和(24)所求得的土壤水总质量ma和总焓H的值以及土壤中的相变平衡关系式(17)就可求解出土壤液态水含量θl, 含冰量θi以及土壤温度

T*三个诊断量.

3.2 模式简化的合理性验证

由方程组(23), (24)构成的简化改进型模式较(9), (10)已经省略了很多项, 明显地简化了计算, 节省了机时. 虽然简化方程是较严格地用量级分析和计算结果推导而得, 理论上讲不应降低精度, 但为了进一步用实际情况证明,本研究利用简化且变换后的模式 (23), (24)模拟结果与简化前模式(9), (10)模拟结果进行比较. 要指出的是D66站的土壤在垂向上本身就是非均质的,本模型适用于非均质的特点正可发挥优势. 比较发现, 对4个站来说, 简化前后的2个模式在模拟不同深度土壤温度和含水量时结果十分一致(图略), 这说明, 简化后的模式与原完全模式具有一样的精度.

3 统一土壤水热传输耦合模型的验证

3.1 资料

本研究用到的资料包括美国明尼苏达州的罗斯芒特(Minnesota/Rosemount)田间实验资料、青藏高原季风试验D66站资料、敦煌双墩子戈壁微气象站观测资料以及黑河试验站观测资料. 4套观测数据包括了不同性质的下垫面: Minnesota/Rosemount田间实验站代表了季节性冻土的特点; D66站既代表永冻土又代表非均质土壤的特点; 敦煌和黑河实验站代表

了干旱地区的特点. 以下仅对两套数据进行简要说明.

3.3 改进型简化模式性能的评价

为了评价改进后的简化通用的统一土壤模型

1532

中国科学 D辑 地球科学

第37卷

((23)和(24)式)的性能, 比较了四个观测站数据和模型的模拟结果. 仅以D66站和敦煌试验站为例, 来具体说明改进型简化模式应用于不同下垫面时的性能.

图4显示了在D66站土壤冻融过程中, 改进型简

化模式模拟的不同深度的土壤温度和液态水体积含量随时间变化的结果与观测的比较. 而图5显示了敦煌试验站2000年6月1日01时~6月9日16时模拟的不同深度土壤温度与观测值的比较情况. 由于敦

图4 改进后的统一土壤模式模拟的1997年10月1日到11月1日青藏高原D66站不同深度的土壤温度((a), (b))与土壤液

态水体积含量((c), (d))同观测值的比较

实线为模拟值, 虚线为观测值

第11期

李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1533

图5 统一土壤模式模拟的2000年6月1日到6月10日敦煌实验站不同深度的土壤温度与观测值的比较

实线为模拟值, 虚线为观测值

煌地区气候十分干旱, 年降水量仅几十毫米, 土壤湿度很低, 感热远大于潜热, 土壤内主要的过程为热传导过程, 土壤表层温度日较差很大, 可达40~50℃.

由图看出, 简化通用的统一土壤模型不仅在变化趋势上而且在数值上都能较好地重现观测结果.

以上几个模拟例子包括了季节性冻土下垫面、永冻土下垫面和干旱下垫面, 土壤也历经了湿度较湿的时间(D66站在10月8日有一次降雨过程), 还牵涉到非均质的土壤(D66站), 这说明模型的通用性很强. 此外, 实验表明, 较之原简化模型(方程(13), (14)), 改

1534

中国科学 D辑 地球科学

第37卷

进后的数值求解方案能使求解(方程(23), (24))过程中的迭代次数减少约三分之一, 大大地节省了运算时间. 所以改进后的统一土壤模式不仅更合理真实地刻画了土壤中的水热传输过程, 还满足了陆面模式和气候模式研究中省时的需要.

参 考 文 献

1 Kurbatkin G P, Manabe S, Hahn D G. The moistening of continents

and the intensity of monsoon circulation. Sov Met Hydrol, 1979, 11(6): 5—11

2 Sud Y C, Fennessy M J. A study of the influence of surface albedo

on July circulation in semi-arid regions using the GLASG CM. J Clim, 1982, 2(2): 105—125

3 Sud Y C, Smith W E. The influence of surface roughness of deserts

on the July circulation—A numerical study. Bound-Layer Me-teor, 1985, 33(1): 15—49

4 Cunningtion W M, Rowntree P R. Simulations of the Saharan at-mosphere: Dependence on moisture and albedo. Q J Roy Meteorol Soc, 1986, 112: 971—999

5 Dickinson R E, Henderson-Sellers A, Kennedy P J, et al. Bio-sphere-atmosphere transfer scheme (BATS) for the NCAR Com-munity Climate Model. NCAR/TN-275+STR. 1986

6 Sellers P J, Mintz Y, Sad Y C, et al. A simple biosphere model (SIB)

for use within general circulation models. J Atmos Sci, 1986, 43(6): 505—531

7 Dai Y, Zeng X, Dickinson R E, et al. The Common Land Model

(CLM). Bull Am Met Soc, 2003, 84(8): 1013—1023

8 Dickinson R E, Henderson-Sellers A, Kennedy P J. Bio-sphere-Atmosphere Transfer Scheme (BATS) Version 1e as Cou-pled to the NCAR Community Climate Model. NCAR/TN- 387+STR. 1993

9 Engelmark H, Svensson U. Numerical model of phase change in

freezing and thawing unsaturated soil. Nord Hydrol, 1993, 24: 95—110

10 Koren V, Schaake J, Mitchell K, et al. A parameterization of

snowpack and frozen ground intended for NCEP weather and cli-mate models. J Geophys Res, 1999, 104(D16): 19569—19585 11 Mölders N, Haferkorn U, Döring J, et al. Long-term numerical in-vestigations on the water budget quantities predicted by the hy-dro-thermodynamic soil vegetation scheme (HTSVS)-Part I: De-scription of the model and impact of long-wave radiation, roots, snow, and soil frost. Meteorol Atmos Phys, 2003, 84:115—135 12 Philip J R. Evaporation and moisture and heat fields in the soil. J

Meteor, 1957, 14: 354—366

13 Rosema A. A mathematical model for simulation of the thermal

behavior of the bare soil, based on the heat and moisture transfer. 11th ed. Delft: Niwars, 1975

14 Sun S F. Moisture and Heat Transport in a Soil Layer Forced by

Atmospheric Conditions. M.S.Thesis. Connecticut: University of Connecticut, 1982

15 Zhao L, Cray D M, Male D H. Numerical analysis of simultaneous

heat and water transfer during infiltration into frozen ground. J Hydrol, 1997, 200: 345—363

16 Flerchinger G N, Saxton K E. Simultaneous heat and water model

of a freezing snow-residue-soil system: I. Theory and development. Trans Amer Soc Agric Engr, 1989, 32(2): 565—571

17 Hansson K, Simunek J, Mizoguchi M, et al. Water flow and heat

transport in frozen soil: Numerical solution and freeze/thaw applications. Vadose Zone J, 2004, 3: 693—704

4 结论和讨论

建立了一个合理统一的土壤水热传输耦合通用模型. 该模型对当前陆面模型中较为不完善的干旱、半干旱区以及冻土区的土壤水热传输耦合过程有真实合理的描述.

(1) 首先, 利用量级分析对完整复杂的土壤水热传输耦合模型进行合理简化, 得到简化的统一土壤水热传输耦合通用模型. 并针对模型在处理土壤冰水相变过程中的不确定性进行了改进, 用土壤总焓和土壤水总质量替代温度和体积含水量作为方程预报量, 设计了一套合理省时的数值计算方案, 建立了统一土壤模型的新体系.

(2) 用土壤水势梯度来描述土壤中液态水的流动, 使得模型对于垂直质地的均匀土和非均匀土都适用.

(3) 在处理干旱区地表蒸发过程时, 放弃了地表局地平衡的假设, 用蒸发面下移(或DSL)的方法来计算表面蒸发, 使得物理过程更加合理明确.

(4) 采用了建立在热力学平衡态基础上的土壤冻融参数化方案, 能真实地再现土壤的冻融变化过程.

尽管本文未考虑植被在土壤水热平衡中的作用, 但是通用的土壤水热传输耦合模型仍然可作为模拟有植被覆盖区土壤水热传输过程的基础. 由于有植被覆盖时, 水从地表输送到大气很大一部分是由植被控制, 因此考虑到植被就需要引入植物根系对土壤水的抽吸作用、植被的蒸腾作用等过程. 除此之外, 土壤内部基本的水热传输过程同无植被覆盖时是完全一样的. 所以发展一个既真实又简化、且适用于多种情景下的通用土壤水热传输耦合模式对于进一步发展土壤-植被-大气相互作用模式至关重要, 也是陆面过程模型发展研究中必不可少的环节之一. 致谢 感谢中国科学院寒区旱区环境与工程研究所杨梅学、胡泽勇、卫国安以及美国明尼苏达大学的John Baker为模式验证提供了可靠的观测资料和相关的重要信息.

第11期

李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1535

18 Harlan R L. Analysis of coupled heat-fluid transport in partially

frozen soil. Water Resour Res, 1973, 9: 1314—1323

19 Guymon G L, Luthin J N. A Coupled heat and moisture transport

model for arctic soils. Water Resour Res, 1974, 10: 995—1001 20 Fuchs M, Campbell G S, Papendick R I. An analysis of sensible

and latent heat flow in a partially frozen unsaturated soil. Soil Sci Soc Am J, 1978, 42(3): 379—385

21 Shoop S A, Bigl S R. Moisture migration during freeze and thaw

of unsaturated soils: modeling and large scale experiments. Cold Regions Sci Tech, 1997, 25 (1): 33—45

22 Slater A G, Pitman A J, Desborough C E. Simulation of

freeze-thaw cycles in a general circulation model land surface scheme. J Geophys Res, 1998, 103: 11303—11312

23 Cherkauer K A, Lettenmaier D P. Hydrologic effects of frozen

soils in the upper Mississippi River basin. J Geophys Res, 1999, 104: 19599—19610

24 Pauwels V R N, Wood E F. A soil-vegetation-atmosphere transfer

scheme for the modeling of water and energy balance processes in high latitudes.1.Model improvements. J Geophys Res, 1999, 104: 27811—27822

25 Cary J W, Mayland H F. Salt and water movement in unsaturated

frozen soil. Soil Sci Soc Am J, 1972, 36: 549—555

26 Clapp R B, Hornberger G M. Empirical equations for some soil

hydraulic properties. Water Resour Res, 1978, 14(4): 601—604 27 Jordan R. A one dimensional temperature model for a snow cover:

Technical Docementation for SNTHERM.89. US. Army Corps of Enigineers, Cold Regions Research and Engineering Laboratory, Special Report 91-16. 1991

28 Deardorff J W. Efficient prediction of ground surface temperature

and moisture with inclusion of a layer of vegetation. J Geophys Res, 1978, 83: 1889—1903

29 Lin J D. On the force-restore method for prediction of ground sur-face temperature. J Geophys Res, 1980, 85: 3251—3254

30 Miller R D. Phase equilibria and soil freezing. In: Permafrost:

Proceedings of the Second International Conference. Washington DC: National Academy of Science-National Research Council, 1965. 287: 193—197

31 Farouki O T. Thermal Properties of Soils. In: Series on Rock and

Soil Mechanics. Vol 11. Germany: Trans Tech Publications, 1986. 102—119

32 Alvenas G, Jansson P E. Model for evaporation, moisture and

temperature of bare soil: calibration and sensitivity analysis. Agric Forest Meteorol, 1997, 88: 47—56

33 Hillel D. Soil and Water: Physical Principles and Processes. New

York: Academic Press, 1971

34 Choudhury B J, Monteith J L. A four-layer model for the heat

budget of homogeneous land surfaces. Q J Roy Meteorol Soc, 1988, 144: 373—398

35 Yamanaka T, Takeda A, Sugita F. A modified surface-resistance

approach for representing bare-soil evaporation: Wind tunnel ex-periments under various atmospheric conditions. Water Resour Res, 1997, 33: 2117—2128

36 Kondo J, Saigusa N, Sato T. A Parameterization of evaporation

from bare soil surfaces. J Appl Meteoro, 1990, 29: 385—389 37 Camillo P J, Gurney R J. A resistance parameter for bare soil

evaporation models. Soil Sci, 1986, 141(2): 95—105

38 杨梅学, 姚檀栋, 勾晓华. 青藏公路沿线土壤的冻融过程及水

热分布特征. 自然科学进展, 2000, 10(5): 443—450

39 张强, 卫国安, 黄荣辉. 西北干旱区荒漠戈壁动量和感热总体

输送系数. 中国科学D辑: 地球科学, 2001, 31(9): 783—792 40 张强, 黄荣辉, 王胜, 等. 西北干旱区陆-气相互作用试验

(NWC-ALIEX)及其研究进展. 地球科学进展, 2005, 20(4): 427—441