大气科学  2015, Vol. 39 Issue (4): 827-838   PDF    
陆面过程模型中垂直非均匀土壤的水分传输及相变的模拟
李倩1, 孙菽芬2     
1 中国科学院大气物理研究所季风系统研究中心, 北京100190;
2 中国科学院大气物理研究所大气科学和地球流体力学数值模拟国家重点实验室, 北京100029
摘要:土壤湿度在陆气相互作用中的重要性体现在它既能影响陆地和大气之间水循环的速率, 又能改变地表的能量分配。本文针对陆面过程模型中描述土壤湿度变化的方程进行了理论分析, 指出在非均匀土壤和冻土中采用土壤水势梯度描述垂直非均匀土壤水分流动的合理性。基于描述土壤内部水热传输的统一土壤模型, 并利用推广的表征土壤水分特征的Clapp-Hornberger关系式, 研究了非冻结和冻结的土壤湿度对于垂直非均匀土壤的敏感性。结果表明, 由土壤质地决定的土壤水势和导水率对土壤湿度的模拟有重要的影响。具体地, 在决定土壤性质的Clapp-Hornberger关系式中, 与土壤质地有关的饱和水势、饱和导水率以及土壤孔隙大小分布指数B, 对土壤湿度的模拟起到了关键作用。参数B的重要性尤为突出, 它的增加会引起导水率的大大下降, 从而对水分在土壤中的垂直分布产生重要影响。饱和水势的绝对值和参数B的增加会使得土壤水势绝对值增加明显, 使土壤的结冰(融化)过程延迟, 土壤温度因为没有结冰(融化)释放(吸收)的潜热加热(冷却)而持续下降(上升), 因此在冻融时期土壤温度会比观测值振幅偏大。上述结果揭示了考虑土壤垂直非均匀性并采用有效的土壤特性参数对于陆面过程模型的重要性。
关键词陆气相互作用     土壤湿度     土壤垂直非均匀性     土壤冻融     Clapp-Hornberger     关系式     土壤孔隙大小分布参数    
The Simulation of Soil Water Flow and Phase Change in Vertically Inhomogeneous Soil in Land Surface Models
LI Qian1, SUN Shufen2     
1 Center for Monsoon System Research, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100190;
2 State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029
Abstract: Soil moisture plays an important role in land-atmosphere interaction because it not only has effects on the water cycle between the land and atmosphere, but also the surface energy distribution. This work theoretically analyzes the equations for soil moisture change in land surface models, and shows the rationality for the equations in which the soil water potential is used to describe the soil water flow in vertically inhomogeneous soil and frozen soil. Based on a Simple Unified Soil Model (SUSM) and the Clapp-Hornberger relationship, the sensitivity of soil moisture in frozen or unfrozen soil to the vertical heterogeneity is also investigated. The results show that soil water potential and hydraulic conductivity have crucial effects on the simulation of soil moisture. The parameter B (soil porosity distribution) in the Clapp-Hornberger relationship is the most important because, when it increases, it leads to decreasing hydraulic conductivity, which results in the soil water vertical distribution. Meanwhile, increases in both the saturated soil water potential and parameter B also lead to a higher absolute value of soil water potential. When soil begins to freeze, this effect would delay soil water from freezing (melting), meaning the soil temperature cannot persistently decrease (increase) due to a lack of release (absorption) of latent heat fluxes originating from the soil freezing (melting) process. Therefore, the amplitude of soil temperature when freezing or melting is larger than observed. The results of this study reveal the importance of soil heterogeneity and the use of effective parameters for soil characteristics in land surface models.
Key words: Land-atmosphere interaction     Soil moisture     Soil vertical heterogeneity     Soil freezing and melting     Clapp- Hornberger relationship     Soil porosity distribution parameter    
1 引言

土壤湿度是控制陆面与大气相互作用的一个重要变量(Shukla and Mintz,1982Koster et al.,2004)。它可以通过对地表蒸发及其他地表能量通量的影响来影响天气过程。土壤湿度的异常可以持续几个月,虽然缺乏足够的观测事实证明土壤湿度对于降水的影响,但是它对降水的作用在大气环流模型(Atmospheric General Circulation Model- AGCM)中是常常可见的(Rind,1982Yeh et al.,1984)。事实上,很多AGCM的研究已经表明,海洋对于夏季中纬度陆地上降水的影响要小于土壤湿度对降水的影响。且考虑了土壤湿度的异常对于降水的季节预报有很好的改进作用(郭维栋等,2007)。

不同复杂程度的陆面模式对土壤湿度的准确模拟关系到陆—气间水平衡和水交换的模拟结果,因此显得十分重要(Mahmood and Hubbard,2003)。在现有众多的陆面过程模型中,对于土壤湿度的定量描述是通过以下两个方面来实现的:(1)在垂向分布均匀的介质基础上,基于推广的达西定律(即Richard方程)来描述土壤水的流动;(2)同时赋以地表土壤不同的水文特征和热力学特性。但在现实中,土壤质地的垂直非均匀分布是很常见的。很多理论研究(Yeh et al.,1985a1985bMantoglou,1992)和实验室(Yeh and Harvey,1990Sassner et al.,1994Destouni et al.,1994)及野外观测(Jensen and Mantoglou,1992)已表明,这样的非均匀性对于非饱和土壤内部的水分输送有很大影响。且在模型中也显现出由于土壤特性的垂直非均匀性而造成土壤导水率、水势等土壤水力性质的非均匀性,进而影响土壤水分平衡和土壤湿度的模拟结果(Mahmood and Hubbard,2003Jhorar et al.,2004)。此外,当土壤发生冻融时,土壤中除了有液态水之外还有固态水—冰,即使含水量一样,冰含量垂直分布的不同也会使土壤的导水率垂直分布不均 匀,从而形成水势大小和梯度的垂直不均匀分布。因此,即使是垂直均质和含水量均匀分布的土壤中,当土壤伴有冻融过程发生时,其水分也会受到由冰含量垂直分布不同而形成的水势梯度力的驱动而发生流动,这相当于强化了土壤非均匀性的特点。

尽管目前部分陆面过程模型在水分的模拟上通过分层和土壤性质的参数化方案来刻画土壤的非均匀性,但是不同的模型对于同一土壤模拟的结果还是有差别的。这不仅与模型选用的土壤水性质的参数化方案有关,还与模型如何处理描述土壤水分变化的方案有关。

在土壤水力性质的参数化方面,土壤水分特征关系式是一个重要的水力学关系,如何确定关系式中的各参数,会对模型的模拟结果产生影响。土壤水分特征关系式包括了土壤基质势与土壤水分的关系以及土壤导水率与土壤水分的关系。且当土壤冻结时,这些关系式中会增加相应的参数(CkE)分别将土壤含冰量对土壤水势和导水率的关系联系起来,因此相对于未冻结土壤其水分特征关系式会有所变化。Zhang et al.(2007)利用模型对CkE的存在对于模拟结果进行了敏感性试验,其结果表明,当Ck为零,即不考虑含冰量对土壤水势的影响时,模型就会高估土壤中冰含量而低估液态水含量,土壤温度的日变化也被低估,潜热通量偏小。同时,如果E为零,即不考虑含冰量对土壤导水率的影响,有时会出现土壤水分流动过大,造成与实际相悖的某层土壤水过饱和的情况。除此之外,其中与土壤本身性质有关的个别参数是以指数形式对土壤基质势及导水率产生影响的。而在目前大多数陆面过程模型中这些与土壤质地有关的参数都是根据经验或者部分实验室结果而得,具体大小对土壤水分的模拟会产生怎样的影响需要进一步探讨,且对发生冻融的土壤会产生怎样的影响也尚未有具体分析。

在处理描述土壤水分变化方案方面,目前常用的处理土壤湿度变化的推广达西定律,即Richard方程是一个高度非线性方程,一般很难有解析解,只有依靠数值解。而在这一求解过程中(包括基本方程的建立和推演,数值求解方法的发展),如果处理不慎、简化不当就会丢失土壤非均匀性的特征对水流的影响。例如,在Richard方程中,土壤的水分流动速率是由与土壤基质势梯度有关的项来驱动的。但在一些模型中为了求解方便,将土壤基质势的梯度项简单地转换成与土壤含水量梯度有关的项(张述文等,2009),正是这一转换就将土壤性质的空间非均匀性和土壤含冰量对土壤水势的影响丢失了。张述文等(2009)的研究曾明确指出,由于土壤质地在垂直方向上的不均匀性,很多陆面过程模型如果使用土壤湿度的梯度进行描述土壤水分变化,模拟出来的土壤水分是连续分布的,但是模拟的土壤基质势出现不连续的情况,这是与实际相悖的情形。而用土壤基质势的梯度进行描述时,模拟出的基质势是连续的,但是土壤水分是空间分布不均的,这与实际相符,因此在描述非均匀土壤中水流运动时不应该利用土壤含水量梯度来表征水流运动。

基于以上提到的不同模型对土壤水力学特性和土壤湿度变化方程处理的不同,也鉴于未有深入的研究指出不同的陆面过程模型对于表征土壤质地的参数选取方面的敏感性,本文中选用了Li and Sun(2008)发展的简化的通用土壤模型对土壤水分模拟的敏感性进行试验研究。研究中选取目前常用的表征土壤水分特征的关系式(Clapp-Hornberger 关系式)及不同土壤水分特性参数,并重点考虑了土壤冻结期间土壤湿度的模拟对各参数的敏感性。

2 非饱和土壤中的土壤湿度方程 2.1 陆面过程模型中常用的非饱和土壤湿度方程

基于土壤水分守恒建立的用于描述土壤非饱和流垂直运动的方程可写为

$\frac{{\partial {\theta _{\rm{l}}}}}{{\partial t}} = - \frac{{{\rho _{\rm{i}}}}}{{{\rho _{\rm{l}}}}}\frac{{\partial {\theta _{\rm{i}}}}}{{\partial t}} - \frac{{\partial ({q_{\rm{l}}})}}{{\partial Z}}$, (1)

如果土壤未冻结,方程写为

$\frac{{\partial {\theta _{\rm{l}}}}}{{\partial t}} = - \frac{{\partial ({q_{\rm{l}}})}}{{\partial Z}}$, (1a)

式中,${\theta _{\rm{l}}}$为土壤单位体积液态水含量(m3 m−3,文中以下简称土壤体积含水量或含水量);${\theta _{\rm{i}}}$为体积含冰量(m3 m−3);${\rho _{\rm{i}}}$和${\rho _{\rm{l}}}$分别为冰和液态水的密度(kg m−3);t为时间(s);z为土壤深度(m),也代表了竖向的坐标(向下为正)。${q_{\rm{l}}}$是液态水流通量(向下为正),遵循非饱和介质中推广的达西定律(即Richard方程),可写为

${q_{\rm{l}}} = {K_{\rm{l}}}\left[ { - \frac{{\partial \psi }}{{\partial z}} + 1} \right]$, (2)

其中,$\psi $为土壤基质势(m),${K_{\rm{l}}}$是土壤的导水率(m s−1);将上式带入(1)式,得到土壤体积含水量的变化方程:

$\frac{{\partial {\theta _{\rm{l}}}}}{{\partial t}} = - \frac{{{\rho _{\rm{i}}}}}{{{\rho _{\rm{l}}}}}\frac{{\partial {\theta _{\rm{i}}}}}{{\partial t}} - \frac{\partial }{{\partial Z}}\left( { - {K_{\rm{l}}}\frac{{\partial \psi }}{{\partial Z}} + {K_{\rm{l}}}} \right)$, (3)

2.2 土壤水分特征关系式

非饱和土壤中,忽略渗透势的影响,土壤水分特征曲线表征了土壤水势和土壤含水量及导水率之间的本构关系。这一本构关系与土壤中土壤颗粒的尺度和分布及孔隙空间的结构有关,随着土壤质地的变化而变化。根据大量的田间和实验室采样试验可以得出不同土壤质地的水分特征曲线关系式。目前常用的表征土壤水分特征曲线的经验关系式基本有三种,即Brooks-Corey 关系式(Brooks and Corey,1966)、van Genuchten关系式(van Genuchten,1980)和Clapp-Hornberger 关系式(Clapp and Hornberger,1978)。这三个关系式都是简洁、方便使用的数学关系式,能与很多的土壤特征曲线拟合得较好。

Brooks-Corey 关系式:

$\left\{ \begin{array}{l} \psi ({S_{\rm{e}}}) = {\psi _{\rm{s}}}S_{\rm{e}}^{ - B},\\ {K_l}({S_{\rm{e}}}) = {K_{{\rm{sat}}}}S_{\rm{e}}^{{\rm{3 + 2}}B}, \end{array} \right.$

van Genuchten关系式:

$\left\{ \begin{array}{l} \psi ({S_{\rm{e}}}) = \frac{1}{\alpha }{(S_{\rm{e}}^{ - {\rm{1/m}}} - 1)^{{\rm{1}} - {\rm{m}}}},\\ {K_{\rm{l}}}({S_e}) = {K_{{\rm{sat}}}}S_{\rm{e}}^{{\rm{1/2}}}{(1 - {(1 - S_{\rm{e}}^{{\rm{1/m}}}{\rm{)}}^{\rm{m}}})^2}, \end{array} \right.$

Clapp-Hornberger 关系式,也是Brooks-Corey 关系式的简化和近似:

$\left\{ \begin{array}{l} \psi ({\theta _{\rm{l}}}) = {\psi _{\rm{s}}}{\left( {\frac{{{\theta _{\rm{l}}}}}{{{\theta _{\rm{s}}}}}} \right)^{ - B}},\\ {K_{\rm{l}}}({\theta _{\rm{l}}}) = {K_{{\rm{sat}}}}{\left( {\frac{{{\theta _{\rm{l}}}}}{{{\theta _{\rm{s}}}}}} \right)^{{\rm{3 + 2}}B}}, \end{array} \right.$

这些关系式中,${K_{{\rm{sat}}}}$为土壤的饱和导水率(m s−1);${\psi _{\rm{s}}},{\theta _{\rm{s}}}$分别为土壤饱和水势(m;为一负值)和土壤孔隙度(m3 m−3);${S_{\rm{e}}}$为土壤有效孔隙度;参数$B,\alpha ,m$均为随着土壤质地而变化的常数(Burdine,1953)。参数B称为土壤孔隙大小分布指数,孔隙 的大小和分布决定了土壤吸附力和水分之间的关系。根据Clapp and Hornberger 的大量试验(1978)证明,砂土的参数B最小,粘土的最大。Clapp-Hornberger 关系式和Brooks-Corey关系式仅当$\psi $<${\psi _{\rm{s}}}$时成 立,且当$\psi $>${\psi _{\rm{s}}}$时${\theta _{\rm{l}}} = {\theta _{\rm{s}}}$。van Genuchten关系式能更好地描述土壤湿度和土壤含水量之间的关系,尤其在土壤接近饱和时,但是由于其复杂的非线性关系,大多数的陆面过程模型并不采用。而用较为简洁的Clapp-Hornberger 关系式。

当土壤发生冻结时,冰的存在对土壤水势和导水率均产生影响,因此在冻土情况下Clapp- Hornberger关系式可推广写为(Koren et al.,1999; Jame and Norum,1980):

$\psi = {\psi _{\rm{s}}}{\left( {\frac{{{\theta _{\rm{l}}}}}{{{\theta _{\rm{s}}}}}} \right)^{ - B}}{\left( {1 + {C_{\rm{k}}}{\theta _{\rm{i}}}} \right)^2}$, (4)

${K_{\rm{l}}}(\theta ) = {10^{ - E{\theta _{\rm{i}}}}}{K_{{\rm{sat}}}}{\left( {\frac{{{\theta _{\rm{l}}}}}{{{\theta _{\rm{s}}}}}} \right)^{ - B}}$, (5)

式中,${C_{\rm{k}}}$为经验常数(有一定取值范围),本文中取为8;E根据Shoop and Bigl(1997)提出的与饱和导水率有关的经验关系式而定。

因此可见,当土壤质地在垂直方向上有改变、呈现明显的非均匀性时,相应的土壤水分特征曲线中的与土壤质地有关的参数(如${K_{{\rm{sat}}}},{\rm{ }}{\psi _{\rm{s}}},{\rm{ }}{\theta _{\rm{s}}}$,BCkE等)都会有所不同。在这里需要指出的是,目前很多陆面过程模型(Engelmark and Svensson,1993Koren et al.,1999Mölders and Walsh,2004)中将(2)式改写为

${q_{\rm{l}}} = \left[ { - D({\theta _l})\frac{{\partial {\theta _{\rm{l}}}}}{{\partial z}} + {K_{\rm{l}}}} \right]$, (2a)

$D({\theta _{\rm{l}}})$为非饱和土壤水分扩散率(m s−2),且,即仅用土壤水分梯度来表征非饱和土壤液态水流动。应该说,在均匀、非冻结土中,土壤水势和土壤含水量均连续变化、且变化方向相同,方程(2a)能合理地描述液态水流动速率的变化情况。但是在非均匀土或冻土中,土壤水势应是空间分布连续的,而含水量在不同土壤质地的交界面处会不连续,发生跳跃,且水势增(减)的方向与含水量增(减)方向有可能不一致,因此由方程(2)决定的水流方向也许与土壤湿度降低的方向相反,所以由含水量梯度决定[即由方程(2a)决定]的水流动并不能正确地反应土壤液态水的流动。且从关系式(4)中可以看出,非均匀土壤中,土壤水势不仅仅与土壤液态水含量有关,当土壤冻结时还与含冰量有关,与影响土壤质地的土壤饱和水势及参数B等一起都会对土壤水势有很大影响。因此如果要用土壤水分扩散率来表征非饱和土壤中液态水流动时,关系式(2)应写为

${q_{\rm{l}}} = {K_{\rm{l}}}\left[ { - {\rm{(}}\frac{{\partial \psi }}{{\partial {\theta _{\rm{l}}}}} \cdot \frac{{\partial {\theta _{\rm{l}}}}}{{\partial z}}{\rm{ + }}\frac{{\partial \psi }}{{\partial {\theta _{\rm{i}}}}} \cdot \frac{{\partial {\theta _{\rm{i}}}}}{{\partial z}}{\rm{ + }}\frac{{\partial \psi }}{{\partial {\psi _{\rm{s}}}}} \cdot \frac{{\partial {\psi _{\rm{s}}}}}{{\partial z}}{\rm{ + }}\frac{{\partial \psi }}{{\partial B}} \cdot } \right.$$\left. {\frac{{\partial B}}{{\partial z}}{\rm{ + }}\frac{{\partial \psi }}{{\partial {C_{\rm{k}}}}} \cdot \frac{{\partial {C_{\rm{k}}}}}{{\partial z}}{\rm{)}} + 1} \right]$, (2b)

为了简化,可将方程(2b)写为

${q_{\rm{l}}} = {K_{\rm{l}}}\left[ { - {\rm{(}}\frac{{\partial \psi }}{{\partial {\theta _{\rm{l}}}}} \cdot \frac{{\partial {\theta _{\rm{l}}}}}{{\partial z}}{\rm{ + }}\frac{{\partial \psi }}{{\partial {\theta _{\rm{i}}}}} \cdot \frac{{\partial {\theta _{\rm{i}}}}}{{\partial z}}{\rm{ + }}{{\left. {\frac{{\partial \psi (B,{C_{\rm{k}}},{\psi _{\rm{s}}})}}{{\partial z}}} \right|}_{{\theta _{{\rm{l,}}}}{\theta _{\rm{i}}}}}{\rm{)}} + 1} \right]$ (2c)

式(2c)中的第一项表示了由于土壤液态水含量的梯度引起的水分流动;第二项表示由于含冰量在垂直方向上的不均匀对水势影响而造成的液态水流动;第三项表征在液态水含量和含冰量一定的情况下,由于土壤本身质地的不同($B$、${C_k}$、${\psi _s}$的不同)而引起的非均匀性对液态水流动的影响。因此推广的Darcy定律中描述土壤水分的垂直运动并不是简单的。所以,在陆面 过程模式中考虑到非均质土和冻土的情况,描述 的液态水流动只能或是直接采用土壤水势梯度()或是采用展开式(2b)才是合理的,但从模型简洁和应用的角度来说,采用土壤水势梯度的表达更合理。

2.3 Clapp-Hornberger 关系式中各参数对土壤水力特性的影响

表征非饱和土壤的水分特征曲线关系式(4)、(5)中,影响水分特征曲线的参数包括${K_{{\rm{sat}}}},{\psi _{\rm{s}}},{\theta _{\rm{s}}}$,BCkE等。鉴于Zhang et al.(2007)已经对CkE对冻融土壤水热输送影响有过相关研究,在这里就不再赘述。仔细观察后发现,参数B是通过指数形式对土壤水力特性进行影响的。因此这里主 要以参数B为例,定量地描述参数B对于土壤水 势和导水率的影响,我们探讨了在不同的土壤湿 度情况下,土壤水势和土壤导水率随参数B的变化情况。

图 1a中显示了当土壤总含水量(${\theta _{\rm{\tau }}}$)为0.4,孔隙度(${\theta _{\rm{s}}}$)为0.45时,不同的液态水含量情况下,土壤水势绝对值的对数($\lg ( - \psi )$)随着参数B的变化而变化的情况。可看出,随着B的不断增大,土壤水势绝对值均有所增加。但是,较干的或含冰量越多的土壤水势绝对值的增幅明显要高于较湿润或者含冰量较少的土壤。同时也意味着,在同一含水量情况下,粘土(B=11.4)比砂土(B=4.05)的土壤水势的绝对值要大。对于土壤导水率来说,随着B的增加,导水率是下降的(图 1b)。且对于较干的或者冻结程度较深的土壤,其随着B增加而下降的趋势比湿润或未冻结的土壤更明显。并且当B越大,导水率随含水量的变化也越大。也就意味着,粘土的导水率明显低于砂土,且当土壤质地偏向粘土或者是粘土时,土壤含水量的变化能引起导水率较大的变化。此外,导水率的变化除了与参数B有关外,饱和导水率也在一定程度上影响了土壤导水率。明显地,土壤导水率随着饱和导水率的增大而增大(图 1c)。即在同一含水量情况下,砂土的导水率高于壤土和粘土。

图 1 (a)土壤水势绝对值的对数()与参数B的关系(,,);(b)土壤导水率的对数()与参数B的关系(,,);(c)土壤导水率的对数()与饱和导水率的对数()的关系(,,)Fig.1 (a)The relationship between and parameter B(,,);(b)the relationship between and parameter B(,,);(c)the relationship between and (,,)

观察关系式(4)还发现,不同质地土壤的饱和水势也对土壤水势有明显影响。可以清楚地看出,当土壤饱和水势绝对值越大,土壤水势的绝对值也越大。即在同一个温度下,土壤颗粒对于土壤水的吸附力就会增大,需要更低的土壤温度才能使土壤水开始冻结。因此掌握了基本的土壤性质参数对土壤水分模拟的影响后,在利用陆面过程模型对垂直非均匀性的土壤进行模拟时,就能较为容易地判断出模型模拟结果的偏差原因,对模型的校验有一定的帮助。

3 模型及数据简介 3.1 简化的统一土壤模型(Simple Unified Soil Model-SUSM)

SUSM是一个基于物理基础建立的可用于气候研究的统一土壤模型(Li et al.,2010)。模型中,采用预报变量替换的方法改写了目前土壤模型中常用的预报方程,即用土壤焓和土壤水总质量分别代替土壤温度和土壤含水量作为控制方程的预报变量,使之能有效地处理目前土壤模型中有冻土时存在数值解法不稳定或不收敛的问题。且SUSM保留了用方程(2)(即)来刻画液态水流通量,能很好地描述由土壤质地本身的空间非均匀性,和由空间非均匀的冻结—融化过程引起的土壤水—热性质的非均匀性对液态水流通量的影响,从而能很好地应对由此引起的土壤中带有非均匀性特点的水热传输过程的模拟(Li and Sun,2008),模拟出寒冷地区有冻—融过程的土壤温度、含水量和含冰量变化。

SUSM中每层土壤的质量平衡和能量平衡方程分别为

$\frac{{\partial ({M_{{\rm{sl}}}})}}{{\partial t}} = - {\rho _{\rm{l}}}\frac{\partial }{{\partial Z}}\left( { - {K_{\rm{l}}}\frac{{\partial \psi }}{{\partial Z}} + {K_{\rm{l}}}} \right) + $$\frac{\partial }{{\partial Z}}\left( {{D_{{\rm{TV}}}}\frac{{\partial {T_{{\rm{sl}}}}}}{{\partial Z}} + {D_{{\rm{\psi V}}}}\frac{{\partial \psi }}{{\partial Z}}} \right)$, (6)

$\frac{{\partial {H_{{\rm{sl}}}}}}{{\partial t}} = \frac{\partial }{{\partial Z}}\left( {{K_{{\rm{eff\_sl}}}}\frac{{\partial {T_{{\rm{sl}}}}}}{{\partial Z}}} \right)$, (7)

其中,方程(6)为质量平衡方程,${M_{{\rm{sl}}}}$为土壤水总质量(为液态水和固态冰的质量之和)(kg m−2)。方程(7)为能量平衡方程,${H_{{\rm{sl}}}}$为土壤总焓(J m−2)(为温度和含冰量的函数),${K_{{\rm{eff\_sl}}}}$为土壤有效导热率(W m−1 K−1),${T_{{\rm{sl}}}}$为土壤温度。方程(6)右边的第一项代表土壤水分流动对土壤含水量变化的影响,第二项代表土壤中由于温度梯度和水势梯度引起的土壤水汽运动对土壤含水量变化的影响。

此外,方程(6)、(7)中有三个未知量(${T_{{\rm{sl}}}}$,${\theta _{\rm{l}}}$,${\theta _{\rm{i}}}$),要使方程体系完整封闭,还需一个关 系式或者限制条件才能决定这三个土壤中的未知量。其中,根据土壤中水—冰—汽三相平衡决定的土壤水势—冰点的函数关系(称为冰点水势方程):

$\psi = \frac{{{L_{{\rm{i,l}}}}{T_{{\rm{sl}}}}}}{{g{T_{\rm{f}}}}}$,

当其与采用的Clapp-Hornberger关系式(4)联立后,可得温度、未冻水含量和含冰量(${T_{{\rm{sl}}}}$,${\theta _{\rm{l}}}$,${\theta _{\rm{i}}}$)三者之间的函数关系:

, (8)

其中,${L_{{\rm{i,l}}}}$为融化潜热,${T_{{\rm{sl}}}}$为土壤温度(°C),$g$为重力加速度(m s−2),${T_{\rm{f}}}$为273.15 K。Li et al.(2010)分析了目前陆面过程模型中常用的几类关系式和限制条件后发现,这一关系式是在基于热力学平衡的基础上建立起来的,比其他的关系式或者限制条件合理,且采用这个关系式后模型的模拟结果比其他的更加与观测接近。

针对方程(6)、(7)、(8),Li et al.(2009)发展了相应的有效的数值解法对其进行求解。试验表明,采用了有效的数值解法的SUSM在运行时更加省时,也更适合于气候模拟的研究。

3.2 数据简介

为了研究季风系统中青藏高原陆面过程与大气的相互作用,从1996年开始,中日科学家进行了“全球能量水分平衡试验—青藏高原亚洲季风试验(GAME-Tibet)”的国际合作项目。在1997~1998年分别进行了第一阶段(预试验)和第二阶段(加强观测)的野外工作,在藏北高原的不同地点分别建立了自动气象站和埋设了土壤温度湿度观测系统,进行相关的观测。我们选择了位于青藏公路66道班附近D66站的观测资料对土壤非均匀性的试验进行验证。

D66站位于青藏高原北部(35°31' N,93°47' E),海拔4560 m,属于大陆性高原气候,年降水量较少。地表植被稀疏,土壤为非均质永冻土(杨梅学等,2000),最深温度观测表明至少有2.6 m的活动层,地下水位2.3 m。D66自动气象站提供了观测高度在1.5 m的每30分钟大气强迫场资料,包括入射短波辐射通量、气温、气压、相对湿度和风速。土壤温度的观测由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。土壤的观测数据每小时自动采集记录一次。

据分析,D66站的土壤类型主要是砂壤土。D66站土壤类型在垂直方向上还有很大的不均匀性。土壤上部导水率比下部大。此外,很多的研究已表明(Zhang et al.,2007Li et al.,2010),D66站数据的可靠性有助于研究季节性冻土的水热输送变化。因此利用这套资料能对模型针对非均匀性土壤水热传输过程的模拟能力进行验证。

4 试验和结果 4.1 SUSM对土壤质地的敏感性试验

1.3节中已表明Clapp-Hornberger关系式中与土壤质地有关的参数(B,${\psi _{\rm{s}}}$,s等)对土壤水势和土壤导水率有较大影响。因此在陆面过程模型中应根据土壤的质地来选取不同的参数才能较为准确地模拟出非饱和土壤中内部的水热输运过程。尤其当土壤呈现明显的垂直非均匀性时,合适的参数选取是模型模拟结果好坏的关键。本文中为了说明SUSM模型能较好地反应出土壤质地的垂直非均匀对于土壤水分传输的影响,且适用于土壤中垂直非均匀的冻融过程,我们设计了几个试验(见表1),利用SUSM在D66站的模拟情况来鉴别模型对土壤质地的敏感性。

表 1 不同试验中简化的统一土壤模型采用的土壤特性参数 Table 1 The parameters for soil texture used by SUSM (Simple Unified Soil Model) in different experiments

三个试验中用到的土壤参数均已列在表1中。其中控制试验能使SUSM模型能很好地模拟出土壤内部不同深度的土壤含水量和土壤温度的变化,包括土壤发生冻结时土壤温度和含水量的日变化特征。在其他参数不变的情况下,试验一中的B由控制试验中的3.86升高到8.86,其主要鉴别模型对于参数B的敏感性。而试验二用于鉴别模型对于饱和水势的敏感性,其饱和水势为-0.881,明显低于控制试验的-0.131。由于土壤孔隙度和土壤饱和导水率对于土壤水分的变化较为直观,且表现上不如参数B和饱和水势的影响明显,因此这里只显示模型对于参数B和饱和水势的敏感性试验结果。

从图2中可看出,控制试验(黑色实线)基本反映出了由于气温的日变化造成的D66站不同深度的土壤温度和湿度的日变化特征。4 cm的土壤体积液态水含量由于温度在0°C上下的日变化而相应地出现日循环波动(图 2b),这与观测结果非常相符。受到日温度波的影响,20 cm的温度在10月13日以后下降到0°C以下并缓慢降低,出现了明显的日变化(图 2d),同时土壤出现了缓慢增长的冰,液态水含量也相应减少(图 2c)。60 cm以下的土壤温度没有明显的日变化,只是缓慢降低(图略)。而试验一中,B的增大使得模拟出的含水量明显增大,并且未能体现出含水量的日变化特征(绿色实线)。因为B的增大使得土壤的导水率明显下降(图 1b),所以土壤水在层与层之间的流动明显减弱,相对维持一个较高的水平(图 2a,b,c)。即使10月9日的一次降水也只有较少的进入表层以下,其余的作为径流流走。同时,B的增大也使得土壤水势的绝对值明显增加(图 1a),土壤对液态水的吸附力增强,开始冻结的温度也降低,因此当土壤温度降到摄氏零度以下后,在4 cm和20 cm处,与控制试验比,即使温度较低土壤所持有的未冻结水含量也较高(图 2b,c),试验一中土壤并未发生土壤的冻结。此外,从图 2d中明显看出在夜间,有结冰过程的控制试验比未发生冻结的试验一模拟的土壤温度高,这主要是由于土壤结冰过程中释放的凝结潜热所致。因为这种情况下,此时控制试验的土壤已经结冰,但试验一中的土壤并未结冰,土壤温度因为没有结冰释放的潜热加热而持续下降,从而在观测的土壤冻结时期,模拟的土壤温度振幅偏大。同理,在白天温度升高后土壤融化需要的热量会使温度有所降低,而试验一中并未能体现出这一土壤内部冻融过程对土壤温度的影响。

图 2 简化的统一土壤模型模拟的D66站(a)表面、(b)4 cm、(c)20 cm深度的土壤液态水含量(m3 m−3)和(d)20 cm处土壤温度(°C)与观测的比较(红色实线:观测;黑色实线:控制试验;绿色实线:试验一;蓝色实线:试验二)Fig.2 Comparison of soil water content(m3 m−3) and soil temperature(°C)at different soil depths at the D66 site between the simulation results and observations:(a)Soil water content at the surface;(b)soil water content at 4 cm;(c)soil water content at 20 cm;(d)soil temperature at 20 cm

试验二中土壤水势绝对值的增加使得土壤开始冻结的温度也下降,结冰过程较之控制试验有所延迟(蓝色实线)。所以表现出在4 cm及20 cm处土壤液态水含量的日变化特征不如控制试验的明显(图 2b,c)。但是仅增加饱和水势绝对值(试验二)不如增加B的值(试验一)对水分和温度的影响明显。因为B是通过指数形式影响土壤水势和导水率的。

以上的试验可以看出,当土壤未冻结时,参数B越大导水率越小,使得土壤水分在上部的流动性减弱,因此试验一模拟的土壤水分相对于观测值及试验二的模拟值偏高。如图 2b中的绿色实线所示。同时,参数B越大饱和水势的绝对值也越大,当温度降到零摄氏度以下时,土壤冻结的临界温度较低,液态水含量明显也偏高。因此关系式(4)中的参数B的选取对于模型在模拟非冻结和冻结土壤时的表现起到很大作用。

同时也看出,SUSM模型在一定程度上能真实地反应出不同质地的土壤中水分和热量的传输,尤其对由冻融过程引起的土壤性质的空间非均匀性的模拟和分析具有很坚实的理论基础和实际意义。

4.2 SUSM对垂直非均匀土壤和垂直均匀土壤的模拟对比分析

以上的三个试验(表 1)中,SUSM模型中土壤上下各层采用的参数都是统一的,如表 1所示。而实际上,土壤的垂直分布不均匀特征在我国青藏高原等地比较明显。利用D66站的数据,我们选取了砂土和粘土这两个对比较明显的土壤质地类型在垂直方向上进行组合,同时为了对比垂直均匀与非均匀的土壤在水热传输方面的不同,我们设计了以下四组不同的试验(表 2所示):

表 2 SUSM模型模拟的四种不同的土壤垂直分层组合及用到的Clapp-Hornberger关系式中的各参数 Table 2 Parameters in the Clapp-Hornberger relationship used for four combinations of soil vertical stratification in SUSM

(一)土壤模型中垂直方向上表现为非均匀性。

组合1:上部(从地表到10 cm深处)为砂土,

下部(从10 cm深处到5 m)为粘土;

组合2:上部为粘土,下部为砂土。

(二)土壤模型中垂直方向表现为均一性。

组合3:上下均为粘土;组合4:上下均为砂土。

其中,砂土和粘土的Clapp-Hornberger关系式中的各参数均来自于根据1845个土壤样本分析得到的具体数值(Tarboton,2003)(见表 2)。与3.1节不同的是,这里不只是探讨模型对于某个参数的敏感性,而是多个参数的共同作用下的模型对于垂直均匀和非均匀土壤的敏感性。

从图中可以看出,组合1中(蓝色实线)上部大、下部小的饱和导水率导致土壤表面及上部的水分往下流,下部的含水量增大(图 3c)。同时,参数B随着深度的增加而增加,也会导致土壤导水率越往下越小,从而水分出现明显的上部偏小(图 3a,b),下部明显偏大的情形(图 3c)。而根据前面的试验一和试验二可看出,参数B对导水率的影响会明显大于饱和导水率对导水率的影响。组合1中饱和水势的绝对值也是随着深度的增加而增大,使得土壤开始结冰的临界温度越来越低,结冰过程延后,在20 cm处即使土壤温度到达零摄氏度以下依然没有结冰。也就意味着模式中的土壤在实际观测当中应该结冰(融化)时没有结冰(融化),所以土壤温度因为没有结冰(融化)释放(吸收)的潜热加热(冷却)而持续下降(上升)因此温度的振幅明显大于观测值。

图 3 SUSM模拟的D66站(a)表面、(b)4 cm、(c)20 cm深度的土壤液态水含量(m3 m−3)和(d)20 cm处土壤温度(ºC)与观测的比较。红色实线:观测;蓝色实线:组合1;绿色实线:组合2;黄色实线:组合3;紫色实线:组合4&Fig.3 Comparison of soil water content(m3 m−3) and soil temperature(°C)at different soil depths at the D66 site between the simulated results and observations:(a)Soil water content at the surface;(b)soil water content at 4 cm;(c)soil water content at 20 cm;(d)soil temperature at 20 cm&

组合2(绿色实线)与组合1在土壤湿度的模拟上完全相反,上部大下部小的参数B,以及上部小下部大的饱和导水率使得土壤导水率随着深度的增加而增大,因此10 cm以上的土壤含水量均比观测值偏大。同时上部饱和水势绝对值较大的粘土其结冰的临界温度也偏低,因此即使在4 cm处当温度小于零摄氏度,且含水量达到0.26都未结冰。而在20 cm处,砂土的饱和水势绝对值明显降低,即使含水量较低(约0.08),土壤依然随着温度的日变化而出现冻结和融化现象。这一日时间尺度上出现的上部未结冰而下部有结冰发生的现象可以说是由于土壤在垂直分布上的不均匀性造成的。在目前大部分模型的模拟结果中是很难见到的。也很少有研究关注。在较短的时间尺度上(日尺度),虽然无足够的观测事实来证明这一现象的普遍性,但是从理论上来说,由于土壤的垂直非均匀性造成的这一现象应该也是存在的,因为上下部土壤性质本身的差别造成了土壤冻结临界温度的差别,那么即使温度有到零摄氏度以下,且含水量较高,而未达到使之冻结的临界温度土壤仍然不会冻结。

所以从上面的分析可以看出,认识到土壤垂直分布的非均匀性以后如何在模型中加以考虑将直接影响到模型模拟的水热传输等过程。

以上分析了模型在模拟上下部不同土壤质地的含水量和土壤温度的变化,为了对比其与垂直均匀的土壤的情形,组合3(黄色实线)试验和组合4(紫色实线)试验分别展示了上下均为粘土和上下均为砂土的具体模拟结果。可看出,对于组合3来说,在上部(图 3a,b)与组合2(绿色实线)的结果接近,都是土壤含水量由于较低的导水率而偏大。而在下部(图 3c),组合3由于饱和水势的绝对值较组合2中的较大,土壤冻结的临界温度较组合2中的砂土偏低,因此土壤并未结冰,即含水量维持不变。相应地,20 cm处的温度振幅也明显偏大。

上下均为砂土的组合4较之组合1来说,明显的差别在4 cm和20 cm处的土壤湿度的模拟上,组合4中下部偏大的导水率使得4 cm和20 cm的含水量分配更为合理,并不像组合2中由于下部小的导水率使得20 cm处偏高。同时,较组合1偏小的饱和水势绝对值,也使得组合4中的土壤在20 cm处即使含水量不及组合1中的含水量高,也可以有冻融现象发生。此外,从与观测的对比上看,似乎上下均为砂土的组合4更接近于观测值,但这并不能完全认为D66站的土壤实质为上下均一的均匀性土壤。只可以说土壤类型应该是偏砂土类型的,因为这还与模型的分层等有关。

5 结论和讨论

本文主要是针对陆面过程模型中非均匀土壤的湿度模拟进行理论分析和敏感性试验研究。关于土壤非均匀性的概念,文中不仅考虑土壤质地的垂直非均匀分布,还包括了土壤发生冻结后有冰存在导致的土壤性质的不均匀性。对陆面过程模型中常用的描述非均匀土壤水分垂直流动的推广的达西定律进行理论分析,指出采用土壤水势梯度表征水分流动的合理性和优越性,也指出目前常用的部分陆面过程模型中采用土壤含水量梯度表征水分流动不合理性,即可能会出现水分流动沿土壤湿度减小、但水势增大(水势绝对值减小)方向的不协调情况。并利用了表征土壤水分特征的经典的 Clapp- Hornberger关系式对土壤水分特性进行了分析。结果表明,关系式中的参数B,饱和水势及饱和导水率对土壤湿度的模拟起到了关键作用。参数B的增加会导致导水率的大大下降。饱和水势的绝对值和参数B的增加会使得土壤水势绝对值增加明显,从而土壤的冻融过程发展的向低温方向偏移,即在同一温度、同一总含水当量下,具有较大值的参数B和饱和水势绝对值的土壤冻结的可能性小而所含的液态水含量较大。

尽管文中我们重点讨论了土壤的非均匀性对土壤水分特征和水分传输的影响。实际上,从图中可以看出,土壤中的温度和水分是相互作用的。例如,组合1和组合3中的土壤在20 cm处未发生冻结,使得此深度的土壤温度未受到土壤冻(融)所释放(吸收)热量的影响,从而出现了比观测值振幅偏大的情况。这一现象明显体现出土壤中的水热传输的相互作用。本文虽主要集中在用Clapp- Hornberger关系式中用到的各参数来讨论土壤非均匀性对水量平衡的影响,实际上,土壤中与热传输紧密相关的土壤导热率也受到土壤质地的影响。王愚等(2013)的工作曾指出过采用不同的土壤热传导方案对青藏高原站点土壤温度的模拟结果有较大影响。土壤的孔隙度,饱和度和晶粒尺寸分布及矿物质含量等都会对土壤的热传导率有影响,当有冰存在时,由于冰的导热系数是水的4倍,因此土壤中冰的存在又会大大增加土壤的热传导率。所以可以看出,土壤湿度和温度的模拟都应该合适准确地考虑到土壤的垂直非均匀性,并且考虑两者是相互影响和作用下进行的。

文中的工作仅仅是利用模式对于土壤垂直非均匀性的敏感性试验来说明土壤性质参数选取的一个重要性,它对于陆面过程模型起到了关键作用。但是在实际应用中,如果考虑到不同地区不同深度的土壤质地非均匀性,那么将是一个繁重的工作,对于计算空间和机时等的要求也会比较大。所以在认识到土壤垂直非均匀性的重要性之后,我们可以在模型中选用一个较为理想的替代方法,将某个尺度的非均匀性土壤看作是一个等效的均匀介质或等效简化的非均匀介质,用一组有效的土壤性质参数来表征这一均匀介质或等效简化的非均匀介质,以求对这一尺度的土壤平均流做出合适的预报。例如3.1节中所提到的控制试验,其中所采用的参数并非根据某一类型土壤的具体特性而定,而可以将这个非均匀土壤看作一个等效简化的均匀介质来对待,所以这一方法就回归到对某一组有效参数的鉴定上。但是这样做仅仅是为了计算和应用的方便,并不能真正地忽视土壤垂直非均匀性在水热传输模拟中的重要性。

此外,本研究仅仅利用了一个没有考虑植被覆盖的土壤模型,当有植被覆盖时,植被的根抽吸作用对于土壤水分的影响也应该考虑。

参考文献
[1] Brooks R H, Corey A T. 1966. Properties of porous media affecting fluid flow [J]. J. Irrig. and Drain. Div. ASCE, 92: 61-90.
[2] Burdine N T. 1953. Relative permeability calculations from pore size distribution data [J]. Journal of Petroleum Technology, 5 (3): 71-78.
[3] Clapp R B, Hornberger G M. 1978. Empirical equations for some soil hydraulic properties [J]. Water Resour. Res., 14: 601-604.
[4] Destouni G, Sassner M, Jensen K H. 1994. Chloride migration in heterogeneous soil: 2. Stochastic modeling [J]. Water Resour. Res., 30: 747-758.
[5] Engelmark H, Svensson U. 1993. Numerical modelling of phase change in freezing and thawing unsaturated soil [J]. Nord. Hydrol., 24: 95-110.
[6] 郭维栋, 马柱国, 王会军. 2007. 土壤湿度——一个跨季度降水预测中的重要因子及其应用探讨 [J]. 气候与环境研究, 12 (1): 20-28. Guo Weidong, Ma Zhuguo, Wang Huijun. 2007. Soil Moisture—An important factor of seasonal precipitation prediction and its application [J]. Climatic Enviromental Research (in Chinese), 12 (1): 20-28.
[7] Jame Y W, Norum D I. 1980. Heat and mass transfer in a freezing unsaturated porous medium [J]. Water Resour. Res., 16: 811-819.
[8] Jensen K H, Mantoglou A. 1992. Application of stochastic unsaturated flow theory, numerical simulations, and comparisons to field observations [J]. Water Resour. Res., 28: 269-284.
[9] Jhorar R K, Dam J C, Bastiaanssen W G M, et al. 2004. Calibration of effective soil hydraulic parameters of heterogeneous soil profiles [J]. J. Hydrol., 285: 233-247.
[10] Koren V, Schaake J, Mitchell K, et al. 1999. A parameterization of snowpack and frozen ground intended for NCEP weather and climate models [J]. J. Geophys. Res., 104: 19569-19585.
[11] Koster R D, Dirmeyer P A, Guo Z, et al. 2004. Regions of strong coupling between soil moisture and precipitation [J]. Science, 305: 1138-1140.
[12] Li Q, Sun S F. 2008. Development of the universal and simplified soil model coupling heat and water transport [J]. Sci. China Ser. D Earth Sci., 51: 88-102.
[13] Li Q, Sun S F, Dai Q D. 2009. The numerical scheme development of a simplified frozen soil model [J]. Adv. Atmos. Sci., 26: 940-950.
[14] Li Q, Sun S F, Xue Y K. 2010. Analyses and development of a hierarchy of frozen soil models for cold region study [J]. J. Geophys. Res., 115, doi:10.1029/2009JD012530.
[15] Mahmood R, Hubbard K G. 2003. Simulating sensitivity of soil moisture and evapotranspiration under heterogeneous soils and land uses [J]. J. Hydrol., 280: 72-90.
[16] Mantoglou A. 1992. A theoretical approach for modeling unsaturated flow in spatially variable soils: Effective flow models in finite domains and nonstationarity [J]. Water Resour. Res., 28: 257-267.
[17] Mölders N, Walsh J E. 2004. Atmospheric response to soil-frost and snow in Alaska in March [J]. Theor. Appl. Climatol., 77: 77-105.
[18] Rind D. 1982. The influence of ground moisture conditions in North America on summer climate as modeled in the GISS GCM [J]. Mon. Wea. Rev., 110: 1487-1494.
[19] Sassner M, Jensen K H, Destouni G. 1994. Chloride migration in heterogeneous soil: 1. Experimental methodology and results [J]. Water Resour. Res., 30: 735-746.
[20] Shoop S A, Bigl S R. 1997. Moisture migration during freeze and thaw of unsaturated soils: Modeling and large scale experiments [J]. Cold Reg. Sci. Technol., 25: 33-45.
[21] Shukla J, Mintz Y. 1982. Influence of land-surface evapotranspiration on the Earth's climate [J]. Science, 215: 1498-1501.
[22] Tarboton D G. 2003. Rainfall-Runoff Processes [R]. Workbook, Utah State University.
[23] van Genuchten M T. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils [J]. Soil Sci. Soc., 44: 892-898.
[24] 王愚, 胡泽勇, 荀学义, 等. 2013. 藏北高原土壤热传导率参数化方案的优化和检验 [J]. 高原气象, 32: 645-653. Wang Yu, Hu Zeyong, Xun Xueyi, et al. 2013. Optimization and test of soil thermal conductivity parameterization schemes in northern Qinghai-Xizang Plateau [J]. Plateau Meteor. (in Chinese), 32: 645-653.
[25] 杨梅学, 姚檀栋, 勾晓华. 2000. 青藏公路沿线土壤的冻融过程及水热分布特征 [J]. 自然科学进展, 10: 443-450. Yang Meixue, Yao Tandong, Gou Xiaohua. 2000. Soil freezing and thawing process and its hydrothermal distribution along the Qinghai-Tibet highway [J]. Progress in Natural Science (in Chinese), 10: 443-450.
[26] Yeh T C, Harvey D J. 1990. Effective unsaturated hydraulic conductivity of layered sands [J]. Water Resour. Res., 26: 1271-1279.
[27] Yeh T C, Wetherald R T, Manabe S. 1984. The effect of soil moisture on the short-term climate and hydrology change—A numerical experiment [J]. Mon. Wea. Rev., 112: 474-490.
[28] Yeh T C, Gelhar L W, Gutjahr A L. 1985a. Stochastic analysis of unsaturated flow in heterogeneous soils: 1. Statistically isotropic media [J]. Water Resour. Res., 21: 447-456.
[29] Yeh T C, Gelhar L W, Gutjahr A L. 1985b. Stochastic analysis of unsaturated flow in heterogeneous soils: 2. Statistically anisotropic media with variable α [J]. Water Resour. Res., 21: 457-464.
[30] 张述文, 李得勤, 邱崇践. 2009. 三类陆面模式模拟土壤湿度廓线的对比研究 [J]. 高原气象, 28: 988-996. Zhang Shuwen, Li Deqin, Qiu Chongjian. 2009. A comparative study of the three land surface models in simulating the soil moisture profile [J]. Plateau Meteor. (in Chinese), 28: 988-996.
[31] Zhang X, Sun S F, Xue Y K. 2007. Development and testing of a frozen soil parameterization for cold region studies [J]. J. Hydrometeor., 8 (4): 690- 701.