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

中国科学D辑:地球科学2007年第37卷第11期: 1522~1535收稿日期: 2007-05-11; 接受日期: 2007-06-25国家自然科学基金(批准号: 40575043, 40233034)和中国科学院知识创新工程(批准号: KZCX3-SW-229)资助项目* 联系人, E-mail: ssf@ 《中国科学》杂志社SCIENCE IN CHINA PRESS通用的土壤水热传输耦合模型的发展和改进研究李倩①②孙菽芬①*(①中国科学院大气物理研究所大气科学和地球流体力学数值模拟国家重点实验室(LASG), 北京 100029;② 中国科学院研究生院, 北京 100047)摘要一个既真实又简化、且适用于湿润与干旱、冻土与非冻土和均质与非均质多种情景下的通用土壤水热传输耦合模式对于陆面过程的模型发展研究至关重要. 研究首先通过量级估计和模型数值模拟结果的分析, 发展了简化且精度较好的土壤水热传输耦合统一模型. 为了克服该模型计算过程中由于需预估冰水相变速率项产生的误差造成的不确定性, 进一步对该简化的统一土壤模式进行变量变换, 以土壤总焓和土壤水总质量替代温度和体积含水量作为方程预报量, 建立了新的通用土壤模型统一体系, 并设计了一套行之有效、省时的数值计算方案. 此模型既可用于一般情况下的裸土, 也可用于较为难处理的非均质土、冻融土壤和干旱土壤等. 与观测结果相比, 改进后的统一土壤模型能很好地模拟出在湿润与干旱、冻土与非冻土和均质与非均质土壤中的水热传输过程. 且由于它的简化, 也适应当今陆面过程模式发展的需要.关键词通用、简化的统一土壤模型水热传输耦合冻土干旱非均质大量的用陆面模式和GCM耦合所进行的敏感性实验[1~4]表明: 发展一个能真实刻画各种下垫面上的陆面过程模式至关重要. 土壤是陆面过程模型要研究的最基本、最重要的陆地下垫面, 它分布广泛、质地繁多且内部物理过程复杂, 在年际、年代际或更长的时间尺度上对各圈层的活动产生很大的影响, 其本身也可经历湿润与干旱、冻结与融化不断交替的变化. 在陆面过程模型的发展过程中, 这部份的模拟一直是重要的关注点, 已得到很大的提高和改善. 但由于认识所限, 目前常用的陆面过程模式对土壤下垫面的考虑仍有以下不足:(1) 大部分模型所考虑的土壤状态主要为湿润的土壤, 土壤参数化方案使用的是忽略了水-热传输耦合作用的等温土壤模式[5,6], 仅分别考虑液态水流动对水分分布和热传导过程对温度分布影响, 忽略了水汽运动对土壤水-热平衡可能的影响, 不能很好用于干旱沙漠和半干旱土壤. 对于冻土冻融相变过程, 也未有很好的参数化方案, 尽管一些有影响的模型都考虑了冻融过程的事实, 但定量描述该物理过程的参数化方案比较简单, 偏离真实的物理关系: 如SSiB和CLM模式[6,7]采用了冻(融)过程和相变热交换仅发生在0℃时(冰点)的简化假设; BA TS模式对冻融过程的考虑相对复杂一些, 假设冻融过程发生在冰点附近某一温度(0~−4℃)范围内[8]. 而真实的冻融过程应是一个连续缓慢的变化过程, 不存在固定的冻融的临界温度(或范围), 必须根据平衡态情况下土壤水势和温度之间的热力学平衡关系以及固有的(但目前只能有实验确定的)土壤水力学特征本构关系来确定土壤含水量、含冰量和温度之间的定量关系.(2) 大部分模型[9~11]基本上采用含水量梯度替代水势梯度的表达式描述水流运动, 使描述水分平衡的方程仅能用于均质土壤, 难于用于非均质的土壤、第11期李倩等: 通用的土壤水热传输耦合模型的发展和改进研究1523冻土物理过程的研究和与饱和地下水运动方程耦合的研究.(3) 几乎现有常用的土壤模型中, 能量和质量平衡方程的描述基本上都以土壤温度、体积含水量作为预报量, 使与含冰量变化率有关的项显性地出现在控制方程中, 在计算迭代求解时对该变化率预估的误差会对冻土温度的计算产生较大偏差, 从而影响计算迭代收敛速度和模拟结果.因此, 要想很好地模拟复杂的地气交换过程, 发展一个通用于冻土和非冻土, 均质和非均质土壤、干旱、半干旱和湿润状态下其内部水热输送过程的统一土壤模式无疑是非常重要的, 而且为了可用于陆面模型和气候模型的发展研究所需, 模型应该力求既简化而又不失真.纵观已经发展的各种模型, 在处理土壤中水热传输过程的复杂程度上均有不同, 可应用的范围也不同[9~25]. 它们又大致可分为三类: 第一类, 完全的土壤水热传输耦合模型[15](10个方程), 完整考虑了水- 冰-气三相的变化及三相之间的相变过程、热传导和水汽扩散相变对土壤能量和质量平衡的影响; 第二类, 相对简化的土壤水热传输耦合模型[16,17](5个方程), 忽略了水汽局地变化的影响; 第三类, 当前常用的简化的土壤水热传输耦合模型. 其水质量平衡方程未考虑水汽扩散对质量平衡的作用[9,10,18~24]. 能量方程中也不考虑水汽扩散相变和液态水运动携带的热量[9,22,23]. 此外, 还有一种介于第二类和第三类之间的土壤模型, 它们考虑了液态水运动携带的热量对土壤能量平衡的影响[18~21,24,25].针对以上的讨论, 显而易见以下几个问题亟需解决: 什么形式的土壤水热输送耦合方程能既简单又不失真地描述在广泛的土壤温度和湿度变化范围内、不同质地的土壤内部其水热传输过程, 并能较好地避免土壤模式计算过程中由于冰-水相变速率的预估误差产生的计算不确定性？不同质地土壤冻融相变过程中, 含水量、含冰量和温度之间合理的相变热力学关系应该如何确定？本文将以完整复杂的土壤水热传输耦合模型为出发点, 考虑到非均质土、湿润土壤、干旱、半干旱土壤及冻土的特殊性, 对以上这几个问题进行仔细的分析研究, 对模型进行分层次的合理简化和数学变换, 得到一个真实合理且适用于陆面模型和气候模型研究的统一通用的土壤水热耦合模型.1非饱和统一土壤水热传输耦合模型1.1完整复杂的土壤水热传输耦合模型Zhao等[15]提出的完整、复杂的土壤水热传输耦合模型可作为本研究的出发点. 它共有10个方程,其中4个预报方程:体积含冰量θ i的控制方程:i,v i,lii0,M Mtθρ+∂+=∂(1) 体积含水量θ1的控制方程:l,v i,ll ll0,M Mut Zθρ−∂∂++=∂∂(2) 水蒸汽质量体积含量v vθρ的控制方程:v v vi,v l,v effM M Dt Z Zθρρ∂∂∂⎡⎤−−=⎢⎥∂∂∂⎣⎦, (3) 描述温度T*(℃)(=T(°K)f(273.15 K))T−=变化的能量控制方程:*v v vi ii,l l,v**vll l l,v eff effC TL Lt t tu T Tc L D KZ Z Z Z Zρθρθρρ∂∂∂−+∂∂∂⎛⎞∂∂∂∂∂⎛⎞=−++⎜⎟⎜⎟⎜⎟∂∂∂∂∂⎝⎠⎝⎠(4) 和5个诊断方程: 由推广的Darcy定律描述液态水流动速率lu的方程:l l1u KZψ∂⎡⎤=−⎢⎥∂⎣⎦＋, (5) 冰点水势方程:i,lf(273.15),L TgTψ−=(6) ψ为土壤基膜水势(这里不考虑土壤溶质势). 决定ψ与体积含冰量iθ和体积含水量lθ之间的推广的Clapp&Hornberger经验关系为[10]()2l0is1,bkcθψψθθ−⎛⎞=+⎜⎟⎝⎠(7)此外还有, 水汽密度vρ与当地水蒸汽压v e的关系;体积水蒸汽含量vθ与土壤空隙度sθ、体积含冰量iθ和体积含水量lθ的关系[15].在以上方程中, Z和t表示土壤的深度和时间; T1524中国科学 D 辑 地球科学第37卷(K)为温度; i,v M ,l,v M 和 i,lM 分别表示冰到蒸汽、液态水到蒸汽和冰到液态水的相变速率; l ρ和i ρ分别为液态水和冰的固有密度; c 1为液态水的热容量; C v 为平均的体积热容量; i,v ,L l,v L 和i,l L 分别为冰到蒸汽、液态水到蒸汽和冰到液态水的相变潜热; K l 为水文导度; K eff 为有效热传导系数; D eff 为水蒸汽在有效空间v θ中的有效扩散系数; g 为重力加速度; 0ψ为饱和水势; b 为Clapp 和Hornberger 常数[26]; c k 为一常数. 以上共有9个方程, 10个未知数(i ,θl ,θv ,θv ,ρ T , l u ,,ψi,v M , l,v M ,i,lM ), 需补充一个条件以封闭方程组. 为此, Zhao 等[15]曾提出一个相变假设关系来封闭方程组, 而Jordan [27]则提出另一组不同的假设关系, 它们的正确性都没有得到证明, 从而降低了模式的准确性. 尽管以上方程体系完备通用, 它可通用于冻土和非冻土, 干旱、半干旱和湿润的状态. 但很显然, 这是十分复杂的方程组, 而且还包含一个正确性并没有得到证明的假设, 并不适用于当今陆面模型的发展.1.2 通用、统一土壤水热传输耦合方程的合理简化可以证明, 根据土壤温度剖面[28,29]变化的规律, 对水汽运动方程中水汽局地项和水汽扩散项的大小进行量级分析, 并利用数值模拟实验的结果, 发现在水汽方程(3)左端的水汽局地变化项v vt θρ∂⎛⎞⎜⎟∂⎝⎠远小于右端的水汽扩散项v eff D ZZ ρ⎛∂⎞∂⎡⎤⎜⎟⎢⎥∂∂⎣⎦⎝⎠. 因此可忽略水汽的局地变化项对土壤水分的贡献1), 故水汽方程(3)退化为v i,v l,veff M M D ZZ ρ∂∂⎡⎤−=⎢⎥∂∂⎣⎦, (8) 利用该方程, 上述方程组中四个预报方程(1)~(4)可初步简化为描述土壤中水分变化的质量方程和描述土壤中温度分布的能量方程:*l i l l i l TV ψV ,u T D D t t Z Z Z Z θθψρρρ⎛⎞∂∂∂∂∂∂+=−++⎜⎟⎜⎟∂∂∂∂∂∂⎝⎠(9)**v i i l i,l l l C T u T L c t t Zρθρ∂∂∂−=−∂∂∂ **l,v TV lv ψV eff T T L D L D K Z Z Z Z Z ψ⎛⎞⎛⎞∂∂∂∂∂+++⎜⎟⎜⎟⎜⎟⎜⎟∂∂∂∂∂⎝⎠⎝⎠, (10) 其中D TV 为温度梯度引起的水汽扩散系数, ψV D 为水势梯度引起的水汽扩散系数. 可以看出, 上述方程组已对复杂完整的方程组作了合理初步的简化, 并去掉了不必要的假定关系, 且仍是通用于均质和非均质土, 冻土和非冻土, 以及不同湿度状态的水热耦合模型.以下将从已初步简化的耦合模型方程(9)和(10)出发, 在不影响精度的前提下, 考虑到不同土壤特性、不同温度变化及不同湿度变化的情况, 对其进行深入的量级分析, 尽可能的简化, 以得到一个普遍适用的土壤水热传输耦合模型.首先, 根据量级分析的方法, 对控制方程(9), (10)中水热通量的各项进行比较, 再根据由该模式计算的水热通量垂直分布大小, 来分析方程中各项的大小关系.可看出, 在控制方程中, 与方程(9)有关的水通量贡献项分别为L ll ,Q K K Zψ∂=−+∂ MV ψV ,Q D Zψ∂=−∂ (11) *TV TV;T Q D Z ∂=−∂ 与方程(10)有关的热通量贡献项分别为*conduct eff,T QH K Z∂=−∂ *convect l l l 1,QH c T K Z ψρ∂⎛⎞=−⎜⎟∂⎝⎠(12) *vaporl,v ψV l,v TV ,T QH L D L D Z Z ψ⎛⎞∂∂=−+⎜⎟⎜⎟∂∂⎝⎠其中, L ,Q MV Q 和TV Q 分别为液态水流通量、由土 壤水势梯度和温度梯度引起的水汽扩散通量; conduct QH , convect QH 和vapor QH 分别为热传导通量、 液态水流动携带的热通量和水汽扩散相变引起的热1) Li Q, Sun S F. Development of a simplified frozen soil model第11期李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1525通量.本文选取了6种不同的土壤: 砂土(sand), 砂壤土(sandy loam), 粉壤土(silt loam), 壤土(loam), 砂黏土(sandy clay), 黏土(clay), 并利用Clapp&Hornberger 提出的经验水势和水文导度关系(所用的各项水热物理参数见表1).在不同的温度, 不同的含水量和含冰量情况下, 计算每种土壤的K 1, TV l /,D ρ ψV l /,D ρ eff ,K l,v ψV ,L D l,v TV ,L D *l l l c K T ρ及l /.ψθ∂∂ 作为典型示范, 仅给出了黏土(表2)和砂土(表3)的结果, 其余从略. 由于在水热通量各项中还与温度梯度*/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 处的日温度波振幅约表1 不同土壤类型的物理参数类型砂土/%a)黏土/%b)θ s c) ψS d) K S e) b f)砂土 92 3 0.339 −0.069 1.07E-06 2.79 砂壤土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℃)各物理量大小θ1/m 3·m −3K 1 /m ·s −1 D TV /ρ1 /m 2·(K ·s)−1D ψV /ρ1 /m ·s −1 K eff /J ·(K ·m ·s)−1L 1V D ψV /J ·m −2·s −1L 1,V D TV /J ·(m ·K ·s)−1ρ1c 1K 1T */J ·m −2·s −1Zψ∂∂ T =0.0℃0.054.30E-32 0.00E+00 0.00E+00 2.40E-01 0.00E+00 0.00E+00 0.00E+00 1.80E+13 0.2 2.20E-16 7.20E-13 7.80E-16 1.20E+00 2.00E-06 1.80E-03 0.00E+005.00E+05 0.466.20E-07 4.00E-14 4.60E-17 1.70E+00 1.10E-07 1.00E-04 0.00E+00 1.40E+01 T =10.0℃0.054.30E-32 0.00E+00 0.00E+00 2.40E-01 0.00E+00 0.00E+00 1.80E-24 1.80E+13 0.2 2.20E-16 1.40E-12 1.60E-15 1.20E+00 4.00E-06 3.50E-03 9.40E-095.00E+05 0.466.20E-077.70E-149.10E-171.70E+002.30E-071.90E-042.60E+011.40E+01T =30.0℃0.05 4.30E-32 0.00E+00 0.00E+00 2.40E-01 0.00E+00 0.00E+00 5.50E-24 1.80E+13 0.2 2.20E-16 4.50E-12 5.70E-15 1.20E+00 1.40E-05 1.10E-02 2.80E-08 5.00E+05 0.466.20E-072.40E-133.10E-161.70E+007.70E-076.00E-047.80E+011.40E+01黏土在冻结时, 不同温度下(<0℃)各物理量大小θ1/m 3·m −3θi/m 3·m −3K 1 /m ·s −1 D TV /ρ1 /m 2·(K ·s)−1D ψV /ρ1 /m ·s −1 K eff /J ·(K ·m ·s)−1L 1V D ψV /J ·m −2·s −1 L 1,V D TV /J ·(m ·K ·s)−1ρ1c 1K 1T */J ·m −2·s −1Zψ∂∂ 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 7.70E+04 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 2.60E+04 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 1.80E+04 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 5.30E+05 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 3.20E+05 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 1.70E+05 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 1.40E+06 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 4.70E+05 0.017 0.25 1.50E-21 7.40E-14 0.00E+00 4.30E+00 0.00E+00 1.90E-04 −1.30E-134.10E+051526中国科学 D 辑 地球科学第37卷表3砂土在非冻结时, 不同温度下(>0 ℃)各物理量的大小θ1/m 3·m −3θi/m 3·m −3K 1 /m ·s −1 D TV /ρ1 /m 2·(K ·s)−1D ψV /ρ1 /m ·s −1 K eff /J ·(K ·m ·s)−1L 1V D ψV /J ·m −2·s −1L 1,V D TV /J ·(m ·K ·s)−1ρ1c 1K 1T */J ·m −2·s −1Zψ∂∂ T =0.0℃0.05 0 7.90E-14 1.50E-12 1.60E-15 8.30E-01 4.10E-06 3.70E-03 0.00E+00 8.00E+02 0.2 0 1.20E-08 7.00E-13 7.90E-16 2.70E+00 2.00E-06 1.80E-03 0.00E+00 4.20E+00 0.33 0 8.50E-074.60E-145.10E-173.40E+001.30E-071.10E-040.00E+006.30E-01T =10.0℃0.05 0 7.90E-14 2.80E-12 3.30E-15 8.30E-01 8.20E-06 6.90E-03 3.30E-06 8.00E+02 0.2 0 1.20E-08 1.30E-12 1.60E-15 2.70E+00 4.00E-06 3.30E-03 4.80E-01 4.20E+00 0.33 0 8.50E-078.60E-141.00E-163.40E+002.60E-072.20E-043.60E+016.30E-01T =30.0℃0.05 0 7.90E-14 8.60E-12 1.10E-14 8.30E-01 2.80E-05 2.20E-02 9.90E-06 8.00E+02 0.2 0 1.20E-08 4.10E-12 5.40E-15 2.70E+00 1.30E-05 1.00E-02 1.50E+00 4.20E+00 0.338.50E-072.70E-133.50E-163.40E+008.70E-076.70E-041.10E+026.30E-01砂土在冻结时, 不同温度下(<0 ℃)各物理量大小θ1/m 3·m −3θi/m 3·m −3K 1 /m ·s −1 D TV /ρ1 /m 2·(K ·s)−1D ψV /ρ1 /m ·s −1 K eff /J ·(K ·m ·s)−1L 1V D ψV /J ·m −2·s −1 L 1,V D TV /J ·(m ·K ·s)−1ρ1c 1K 1T */J ·m −2·s −1Zψ∂∂ 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 7.70E+04 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 2.80E+04 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 1.80E+04T =−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 5.30E+05 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 3.20E+05 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 1.70E+05 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 1.40E+06 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 7.10E+05 0.017 0.25 1.50E-21 7.40E-14 0.00E+00 4.30E+00 0.00E+00 1.90E-04 −1.30E-134.10E+05为1()O e T −∆≈5~11℃, 所以在土壤表层附近温度梯度的量级大小约为O (10/10 cm)=℃O (100/m)℃. 对于表层的水势梯度来说, 可根据l l O Z θψθ⎛⎞∂∂⋅⎜⎟∂∂⎝⎠来确定, 一般来说lO Z θ∂⎛⎞⎜⎟∂⎝⎠ ≈O (0.1/10 cm)≈O (1/m), 所以lO O Zψψθ⎛⎞∂∂⎛⎞≈⎜⎟⎜⎟∂∂⎝⎠⎝⎠. lO ψθ⎛⎞∂⎜⎟∂⎝⎠可根据各种土壤的水力学特征关系式来确定.为了进一步分析比较水热通量的各项, 我们利用较为完整的模式(9)和(10)计算了土壤在冻融阶段和非冻融阶段水热通量的变化情况. 图1和2给出了模式模拟的D66站的土壤在未冻结-冻结-融化时期几个典型时刻的水热通量垂直分布.从以上的量级分析及数值模拟实验结果都可以看出: 对于与热通量有关的各项来说, 当土壤非常干 或者冻结比较厉害的时候, *l l l (/)c T K Z ρψ∂∂总是比*eff (/)K T Z ∂∂小2个量级以上. 而当土壤接近饱和时, *eff (/)K T Z ∂∂与*l l l c T K ρ的量级差别不会很 大. 因此为了模型的完整性*l l l (/)c u T Z ρ∂∂可保 留. 且对于土壤表层来说, 液态水运动携带的热 通量*l l l (/)c u T Z ρ∂∂相当于地表入渗携带的热通量. 计算表明, 当日平均降水量大于10 mm 时, 入渗携带的热通量其量级与*eff (/)K T Z ∂∂相当. 因此在土壤的能量平衡方程中应考虑液态水运动携带的热通量对温度分布的影响. vapor QH 无论在什么情况下都很小. 因此, 在能量方程中可以忽略QH vapor 的作用. 观第11期李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1527图1 模拟的D66站土壤在未冻结-冻结-融化阶段不同时刻的水通量(water flux)垂直分布. QTV , QMV 和QL 分别代表Q TV ,Q MV 和Q L未冻结阶段: (a), (b); 冻结阶段: (c), (d); 融化阶段: (e), (f)察与水通量有关的各项, 发现不同的土壤在不同的含水(冰)量情况下,*TV l D T Zρ∂∂和l K Z ψ∂∂相比, 两者的量级大小并不确定, 由于土壤湿度不同, 前者有可能小于、等于或大于后者. 所以在水质量平衡方 程中TV Q 应该保留. 同样对于ψV l D Zψρ∂∂来说, 在冻结情况下, 由于冻土空隙中的水汽密度只是温度1528中国科学 D 辑 地球科学第37卷图2 模拟的D66站土壤在未冻结-冻结-融化阶段不同时刻的热通量(heat flux)垂直分布hflux_conv, hflux_vap 和 hflux_cond 分别代表了QH convect , QH vapor 和QH conduct . 未冻结阶段: (a), (b); 冻结阶段: (c), (d); 融化阶段: (e), (f)的函数, 所以ψV D 为零. 而土壤未冻结时,ψV l D Zψρ∂∂在土壤湿度不太干情况下, 相比于液态水通量 可以省略, 但在土壤变得很干时, 虽然它和lK Zψ∂∂ 项都很小, 但它的相对重要性逐渐显现出来. 所 以, 在质量平衡方程中该水汽通量的作用也可保留.因此可以得到合理简化的统一土壤水热传输耦合方程:第11期李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1529l i i l l l K K t t Z Z θρθψρ∂∂∂∂⎛⎞=−−−+⎜⎟∂∂∂∂⎝⎠*TV ψV l 1T D D Z Z Z ψρ⎛⎞∂∂∂++⎜⎟⎜⎟∂∂∂⎝⎠, (13) ***v i i l i,l eff l l .C T u T T L K c t t Z Z Z ρθρ⎛⎞∂∂∂∂∂−=⎜⎟⎜⎟∂∂∂∂∂⎝⎠－ (14) 除了以上2个预报方程外, 模型中还有3个诊断方程: (5)~(7)或(15)i l (,)ψψθθ=. (15)这样, 5个方程(13), (14), (5), (6)和(7)或(15), 5个未知变量(θi , l θ, T *, l u , ψ )构成了一个简化的通用土壤水热传输耦合模型.需要指出的是, 部分模式[9~11]用了含水量梯度来描述液态水流动:l .D Z Zθ∂∂⎛⎞⎜⎟∂∂⎝⎠其中D 是液态水的扩散系数(m 2/s), 且ll.D K ψθ∂=∂ 应该明确, 在均匀土中以上方程是可以合理地描述液态水的变化情况, 但是在非均匀土或土壤冻结时, 由于土壤分层的交界处含水量不连续, 水流方向也许与土壤湿度降低的方向完全相反, 含水量的梯度并不能真正地反应土壤液态水的流动. 且当土壤发生冻结时, 土壤水势不仅是液态水体积含量的函数, 还是含冰量的函数, 因此推广的Darcy 定律中l i l i z Z Z θθψψψθθ⎛⎞⎛⎞∂∂∂∂∂=+⎜⎟⎜⎟∂∂∂∂∂⎝⎠⎝⎠, 而不是l l z Z θψψθ⎛⎞∂∂∂=⎜⎟∂∂∂⎝⎠. 所以, 在统一土壤模式中考虑到非均质土和冻土的情况, 采用土壤水势梯度来描述的液态水流动(lK Z Zψ∂∂⎛⎞⎜⎟∂∂⎝⎠)是更合理的. 1.3 统一土壤水热传输耦合模型中的冻融参数化方案在土壤模型中, 方程(6)和(7)或(15)是决定冻结土壤中体积含水量、含冰量和土壤温度之间定量关系的两个重要关系式. 其中, 方程(15)是描述冻土和非冻土中土壤水势和土壤含水量(或含冰量)之间的水力学特征关系, 也称为非饱和冻土中的本构关系,虽然目前这一关系式只能通过实验得到. 而方程(6)是在热力学平衡态基础上, 根据土壤中冰-水-汽三相之间的相平衡关系严格推导得到的冻融土壤中水势和温度之间的关系式. 在冻土中, 如果接受系统处于热力学平衡态的假设 (这是必须遵循的基本物理关系), 联立(15)和(6)就能得到非饱和土壤中关联含水量、含冰量和温度的冻融相变平衡关系: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)便可得冻融相变平衡关系:()1*2i,l l s i o f 1,bk L T c g T θθθψ−−⎡⎤⋅=⋅+⎢⎥⎢⎥⎣⎦(17)式中的b , o ψ和s θ是依赖于土壤特性的参数.从式(17)可看出, (1)土壤的冻结和融化是缓慢变化的过程, 并不是在某一固定温度(如0℃)或某一温度区间内就能完成的. 且即使在温度很低的时候非饱和冻土中仍然有液态水存在. 这与很多室内外观测事实吻合. (2)在(17)式中令i θ=0.0, 可得到某一温度T *(零摄氏度以下)时土壤中存在的未结冰的最大液态水含量max l θ. 即在温度T *的情况下, 只有当土壤中的液态水大于max l θ时多于的水才会冻结. 此外, 变换(17)式可得()l max2*0f i i,l s1,bk g T T c L θψθθ−⎛⎞⎜⎟=+⎜⎟⎝⎠(18a)令i θ=0.0, 也能确定土壤开始冻融的临界温度*critT :1530中国科学 D 辑 地球科学第37卷*0f criti,l s ,bT g T T L ψθθ−⎛⎞=⎜⎟⎝⎠(18b)此时的l θ等于总的积水当量θT . 上式表明, 如果土壤温度低于临界温度*crit,T 土壤开始结冰, 相反则开始融化. 明显地, 土壤冻融的临界温度并不是在零摄氏度.(3) 由于各种土壤的0ψ和s θ都不相同, 因此由(17)式所确定的max l θ和*critT 也有差异. 图3展示了不同质地的土壤在不同温度下所能含有的未结冰的最大液态水含量.图3 不同类型的土壤在0℃以下时, 土壤所能含有不结冰的最大液态水含量纵坐标表示土壤在未结冰的情况下所能含有的最大液态水含量可看出, 对于所有的土壤类型, 液态水含量随温度的变化在冰点附近的某一温度区间内(0~−2℃)最显著. 在−2℃时, 对于砂土, 大部分的液态水会变成冰, 因此可以规定土壤的冻融过程发生在冰点附近的某一区间内, 但对于其它土壤而言, 液态水变成冰的过程还远未完成. 对于黏土, 即使到−20℃, 土壤所能含有不结冰的最大体积液态水含量θ1仍可达0.22. 这说明假定一个统一的结冰区间可能会产生误差. 在统一的非饱和土壤水热传输耦合模型中, 采用土壤冻融参数化方案(17)可能会较为真实普遍地描述土壤的冻融过程.1.4 统一土壤水热传输耦合模型中的蒸发方案对于干旱土来说, 由于地表附近土壤水和温度的梯度都很大, 地表面空隙中的水汽压与当地平衡饱和水汽压差别也会很大, 并不能达到局地平衡 态[32]. 所以在计算地表蒸发的时候应该放弃局地平衡的假设. 针对这一情况, 一些早期的研究在蒸发公式中引入了土壤表面阻抗的概念[14], 对蒸发公式进行了修正. 很多陆面过程模式(SIB, BA TS)都将其引入了表层蒸发的参数化方案中, 它有助于我们理解和有效地描述土壤在干旱时表层蒸发的机制. 但是表面阻抗一般都是土壤表层含水量的经验函数, 并不是在所有情况下均适用.研究表明[33], 在土壤比较干时, 表面会形成干表面层(Dry Surface Layer-DSL), 蒸发的水汽从DSL 的下边界经过DSL 被输送到土壤表面. 之后, 很多研究工作中[34,35]也利用DSL 或蒸发面下移这个概念来计算干旱土壤蒸发的水汽通量, 较为合理明确地解决了干旱土壤蒸发的问题.干旱土壤的蒸发主要有以下两个物理过程, 首先水汽通过分子扩散的形式从干表面层传输到土壤表面, 这一过程中的水汽通量可表示为*a effl ,e sehq q E D z ρρ−= (19)*,e q h 分别为蒸发面上土壤空隙中空气的饱和比湿和相对湿度, s q 为土壤表面上的空气比湿, e z 是干表面层的厚度. 第二个过程, 水汽从土壤表面进入大气, 通量为a s al aq q E r ρρ−=, (20)a q 为大气的比湿, a r 为空气动力学阻抗.Kondo 等[36]指出这两个过程中的水汽通量应该是相同的, 因此联立以上两式, 可以得到实际蒸发通量:*a al a eff/e e hq q E r z D ρρ−=+, (21) 其中, eff /e z D 等价于早期研究中[14,37]提出的表面阻抗因子. 即:s eff /e r z D =, (22)只是此时的s r 不再是土壤表层含水量的经验函数. 所以, 在模式中, 采用以上的蒸发面下移或DSL 的方法来计算表面蒸发, 避免了经验关系式, 不仅适和于不同质地的土壤, 并且其物理过程也很明确.2 简化通用的土壤水热传输耦合统一模型的改进可看出, 上述通用的简化耦合模型其控制方程是以土壤温度、体积含水量作为预报量的非线性方程,第11期李倩等: 通用的土壤水热传输耦合模型的发展和改进研究 1531含冰量变化显性地出现在控制方程中. 当土壤发生冻融时, 能量平衡方程(14)中包含冰水相变产生的潜热项i,l i i L t ρθ∂, 由于i,l L 值很大(i,l L =3.336×105 J/kg ), 故冰水相变产生的巨大潜热将强烈地影响土壤温度的改变, 而一般在求解以上的控制方程时都是采用迭代法, 在每次迭代开始需对含冰量值进行估计, 其估计的偏差将直接影响土壤温度的正确计算, 进而影响冰水相变速率自身的正确计算, 从而对模式的模拟精度和计算速度产生影响.为了避免相变速率的不确定性引起的模式计算不稳定性, 本工作对当前通用的土壤水热传输控制方程中的预报变量进行变量变换, 用总土壤水质量a m (a m =i i ρθ+l l ρθ)和总焓H (H =*v C T −i,l i i L ρθ)作为预报变量, 即将原来方程中与冰水相变有关的项均合并到总土壤水质量和土壤总焓中, 发展出新的预报方程组. 具体方程如下:l l l a m K K t Z Z ψρ∂∂∂⎛⎞=−−+⎜⎟∂∂∂⎝⎠*TV ψV l 1,T D D Z Z Z ψρ⎛⎞∂∂∂++⎜⎟⎜⎟∂∂∂⎝⎠(23)**l eff l l ,u T H T K c t Z Z Z ρ⎛⎞∂∂∂∂=⎜⎟⎜⎟∂∂∂∂⎝⎠－ (24) 这样, 每次对方程(23)和(24)迭代时无需对冰水相变速率进行估计就能得到土壤水总质量和土壤总焓, 然后再根据方程(23)和(24)所求得的土壤水总质量m a 和总焓H 的值以及土壤中的相变平衡关系式(17)就可求解出土壤液态水含量l θ, 含冰量i θ以及土壤温度*T 三个诊断量.3 统一土壤水热传输耦合模型的验证3.1 资料本研究用到的资料包括美国明尼苏达州的罗斯芒特(Minnesota/Rosemount)田间实验资料、青藏高原季风试验D66站资料、敦煌双墩子戈壁微气象站观测资料以及黑河试验站观测资料. 4套观测数据包括了不同性质的下垫面: Minnesota/Rosemount 田间实验站代表了季节性冻土的特点; D66站既代表永冻土又代表非均质土壤的特点; 敦煌和黑河实验站代表了干旱地区的特点. 以下仅对两套数据进行简要说明.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. 土壤的观测数据每小时自动采集记录一次.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].3.2 模式简化的合理性验证由方程组(23), (24)构成的简化改进型模式较(9), (10)已经省略了很多项, 明显地简化了计算, 节省了机时. 虽然简化方程是较严格地用量级分析和计算结果推导而得, 理论上讲不应降低精度, 但为了进一步用实际情况证明,本研究利用简化且变换后的模式 (23), (24)模拟结果与简化前模式(9), (10)模拟结果进行比较. 要指出的是D66站的土壤在垂向上本身就是非均质的,本模型适用于非均质的特点正可发挥优势. 比较发现, 对4个站来说, 简化前后的2个模式在模拟不同深度土壤温度和含水量时结果十分一致(图略), 这说明, 简化后的模式与原完全模式具有一样的精度.3.3 改进型简化模式性能的评价为了评价改进后的简化通用的统一土壤模型。