1 引言

青藏高原 (下称高原) 地处全球气候变化最活跃的地区, 是北半球气候变化的启动器和调节器 (冯松等, 1988)。高原的气候变化不仅直接驱动中国东部和西南部气候的变化 (吴国雄等, 2013), 而且对整个亚洲和北半球甚至对全球的气候变化也具有明显的敏感性、超前性和调节性 (肖志祥和段安民, 2015)。另外, 高原在亚洲季风中扮演重要角色 (Krishnamurti and Ramanathan, 1982), 相对大气而言冬季是冷源, 夏季是热源 (叶笃正等, 1957)。由于独特的地理环境, 青藏高原有特殊的下垫面, 它们的存在和变化对气候变化有明显的反馈调节作用, 而且对生态环境也有着重要影响。随着全球气候变化的深入研究, 高原地-气耦合系统变化及其全球气候效应的重要性越来越明显, 已成为重要的国际气候研究和地球系统科学研究的前沿。高原的地气物理过程对全球气候异常与东亚大气环流及中国的灾害性天气的发生、发展都有着重大影响 (温克刚等, 1999)。高原的热力强迫影响着亚洲季风的形成和时空变化。因此, 高原气候变化及其影响研究已成为地球系统圈层相互作用的重要前沿领域。另外, 高原的热力作用又是由地-气相互作用决定的, 因此深入研究土壤有关物理参数 (土壤热扩散率、土壤热传导率、土壤热通量) 有利于更好地分析能量分配 (Wu et al, 1999), 从而为利用模式模拟地-气相互作用和天气过程提供参数化方案 (Gao et al, 2000)。

土壤热通量是地表能量平衡中的重要组成部分, 土壤表面热量的吸收和损失都与气候变化息息相关 (徐兆生和马玉堂, 1984)。土壤的热力学特性影响着地表和大气间的物质和能量交换, 对地球大气边界层物理过程、大气环流和区域气候产生重要影响。另外, 土壤温度、热扩散率及热通量作为重要的热力指标, 是陆面模式中重要的输入量 (叶笃正等, 1991)。由于地理位置的不同, 地表接受的太阳辐射存在差异, 加之不同的生态系统、环境条件、土壤特性都会造成土壤热通量的时空差异 (马柱国等, 2000)。用模式或观测数据分析地球表面的能量平衡都需要正确估计土壤的热通量和土壤表面的温度, 这些都迫切需要深入认识土壤的热性质。而且对土壤热通量特征的研究, 可以使我们深入了解其物理过程, 为天气、气候及生态模型的建立和检验提供参数和基础理论支持 (韦志刚等, 2005)。

由于受观测资料的限制, 国内外对长时间尺度土壤热通量变化特征的研究工作并不多。过去对青藏高原热源的研究, 多是从地表能量平衡的角度来确定高原热源的强度。徐祥德等 (2001)在能量平衡中多未考虑土壤热流的贡献, 即使考虑土壤热流的贡献, 也仅仅取其气候平均值, 其年际变化甚至认为是零。Tang and Gao (1997)研究中发现, 土壤热流也有其时间变化, 且其最大变化幅度比气候平均值大2~3个数量级。高国栋和陆渝蓉 (1981)计算了我国地表面与下层土壤之间的热交换量和分布。杨振刚 (1985)在均质土壤常温带内和温度变化在余弦波形假定条件下计算了土壤热流值。翁笃鸣等 (1994)利用地表热量平衡各分量的气候计算方法, 集合全国热平衡站资料, 计算并讨论了地表热源在全国的时空分布特征, 分析了地面对大气加热的特点以及地表热源的参数化问题。杨红娟等 (2009)利用遥感地表温度模拟土壤热通量。张宏等 (2012)分析了我国土壤热通量的季节和空间变化规律。以上学者基本上是基于20 cm以上浅层土壤进行的研究, 对深层80 cm以下土壤热通量的研究并不多。本文拟采用热传导对流法计算了青藏高原深层土壤热通量, 并初步分析其时空分布特征及其与高原季风的关系。

2 土壤资料和计算方法 2.1 土壤资料

利用青藏高原的39个观测站 (图 1) 实测的0. 8 m和3.2 m的月平均土壤温度资料计算土壤热通量, 资料长度为1980-2001年。由于青藏高原各站点的分布较少且不均, 而且资料长度不同, 本文选取高原中东部地区资料完整、空间分布较均匀的5个气象站作为代表站点进行分析 (图 1圆点), 分别为冷湖、安多、泽当, 刚察、西宁5个站。分别计算了各个站点的传导、对流及总热通量, 并分析各站土壤热通量的年变化和年际变化特征。

2.2 计算方法

计算土壤热通量首先需要确定热力学参数:土壤温度、热扩散率、土壤体积热容量, 其中土壤的体积热容是由土壤的成分决定 (Van Wijk et al, 1963)。计算土壤热通量的方法有振幅法、相位法、谐波法、热传导对流法。Gao et al (2005)验证了热传导对流法的可靠性, 认为热传导对流法能更好地描述土壤中的热传输过程。

2.2.1 热传导对流法

实际上, 土壤的热传输包含了传导过程和多孔介质中的对流过程 (de Silans et al, 1996)。范新岗和汤懋苍 (1994)及高志球等 (2002)指出土壤温度不仅受土壤热传导作用的影响, 还与土壤中水分的热对流有关 (Gao et al, 2003), 并在以上理论的基础上发展了一套包括热传导和热对流的土壤热通量计算方法, 该方法可以利用气象站地温资料反演获得地热流场信息。从热传导方程和热对流方程出发, 推导出了包括传导和对流的线性传热方程:

$ \frac{{\partial T}}{{\partial t}} = K\frac{{{\partial ^2}T}}{{\partial _z^2}} - W\frac{{\partial T}}{{{\partial _z}}} , $

(1)

式中: T是土壤温度 (单位: ℃), t为时间 (单位: s), K是土壤的热扩散率 (单位: m²· s^-1), $ - W\frac{{\partial T}}{{{\partial _z}}} $表示土壤的对流过程, W表示垂直对流速度 (向下为正)(单位: m·s^-1), 在半无界空间, 给定波动边界$T\left({{z_0}, t} \right) = \bar T{|_{z = 0}} + A\sin \left({{\omega}t + \beta } \right) $, 可以推导出式 (1) 的解析解和土壤热扩散率K与垂直对流速度W的计算公式 (Gao et al, 2005):

$ \begin{gathered} T\left({z, t} \right) = \bar T\left(z \right) + A\exp \left[ { - M\left({z - {z_0}} \right)} \right] \hfill \\ \times \sin \left[ {{\omega}t + \beta + N\left({{z_0} - z} \right)} \right] , \hfill \\ \end{gathered} $

(2)

$ K = \frac{{{{\left({{z_1} - {z_2}} \right)}^2}{\omega}\ln \left({{A_2}/{A_1}} \right)}}{{\left({{\phi _2} - {\phi _1}} \right)\left[ {{{\left({{\phi _2} - {\phi _1}} \right)}^2} + {{\left({\ln \left({{A_2}/{A_1}} \right)} \right)}^2}} \right]}} , $

(3)

$ W = \frac{{{\omega}\left({{z_1} - {z_2}} \right)}}{{{\phi _2} - {\phi _1}}}\left[ {1 - \frac{{2{{\left({\ln \left({{A_2}/{A_1}} \right)} \right)}^2}}}{{{{\left({{\phi _2} - {\phi _1}} \right)}^2} + {{\left({\ln \left({{A_2}/{A_1}} \right)} \right)}^2}}}} \right] , $

(4)

式中:

$ M = - \frac{W}{{2K}} + \frac{{\sqrt 2 }}{{4K}}\sqrt {{W^2} + \sqrt {{W^4} + 16{K^2}{{\omega}^2}} } , $

(5)

$ N = \frac{{\sqrt 2 w}}{{\sqrt {{W^2} + \sqrt {{W^4} + 16{K^2}{{\omega}^2}} } }} , $

(6)

式中: ω是地球的自转角速度, $ {\omega} = 2{\pi /}p\left({{\text{rad}} \cdot {{\text{s}}^{ - 1}}} \right)$, p为z₀处土壤温度的变化周期。A₁、A₂、ϕ₁、ϕ₂分别是深度z₁、z₂处的土壤温度振幅和相位 (z₁ < z₂)。最后得到计算总热通量的公式为 (热流向上为正):

$ \left\{ \begin{gathered} {Q_d} = {C_m}K\frac{{\Delta T}}{{\Delta z}} \hfill \\ {Q_v} = - {C_w}W\Delta T \hfill \\ {Q_t} = {Q_d} + {Q_v} \hfill \\ \end{gathered} \right. , $

(7)

式中: Q_d, Q_v, Q_t分别为传导热通量, 对流热通量和总的热通量 (W·m^-2)。C_w为液态水的容积热容为4.186×10⁶ J·(m³·K)^-1。C_m为土壤的容积热容, 它的取值与季节无关, 主要由土壤类型决定 (董文杰和汤懋苍, 1992)。因青藏高原大部分地区为壤土, 所以孔隙度为40%~50%。根据魏东平等 (1992)的研究取C_m的值为1.256×10⁶ J·(m³·K)^-1。

利用0. 8 m和3.2 m两层地温资料, 由最小二乘法拟合求得振幅和相位, 再利用式 (3)、(4) 求出热扩散率和垂直对流速度, 最后算得土壤热通量。

2.2.2 振幅法和相位法

假设土壤的体积热容C_g(单位: J·cm^-3·K^-1) 和热传导率$\lambda $(单位: W·m^-1·K^-1) 不随深度变化, 则热传导方程为:

$ \frac{{\partial T}}{{\partial t}} = k\frac{{{\partial ^2}T}}{{\partial {z^2}}} , $

(8)

式中: k为土壤热扩散率 (单位: m²·s^-1), $ k = \lambda /{C_g}$。

给定土壤温度的上边界条件:

$ T\left({{z_0}, t} \right) = \bar T\left({{z_0}} \right) + A\sin \left({{\omega}t + \varphi } \right) , $

(9)

式中:$\bar T\left({{z_0}} \right) $为深度为z₀处的年平均土壤温度, A为z₀处的土壤温度振幅, φ为z₀处的土壤温度相位。根据边界条件式 (9) 可以得到式 (8) 的解析解:

$ T\left({z, t} \right) = \bar T\left(z \right) + A\exp \left({ - Bz} \right)\sin \left({{\omega}t + \varphi - Bz} \right) , $

(10)

$B = \sqrt {{\omega}/2k} $来表示深度对振幅和相位的影响, 其中$ \bar T\left(z \right)$为土壤的深度z处的年平均温度。由式 (9) 和式 (10) 可知, 任意深度的土壤温度都可以用一个正弦函数来表示, 则深度为z₁和z₂的土壤温度实测值可以表示为:

$ T\left({{z_1}, t} \right) = \bar T\left({{z_1}} \right) + {A_1}\sin \left({{\omega}t + {\varphi _1}} \right) , $

(11)

$ T\left({{z_2}, t} \right) = \bar T\left({{z_2}} \right) + {A_2}\sin \left({{\omega}t + {\varphi _2}} \right) , $

(12)

式中: $\bar T\left({{z_1}} \right) $和$ \bar T\left({{z_2}} \right)$分别是深度z₁、z₂处的年平均值土壤温度, z₁ < z₂, A₁、A₂和φ₁、φ₂分别对应土壤温度振幅和相位。振幅A可以通过 (13) 式计算:

$ A = \frac{{{T_{\max }} - {T_{\min }}}}{2} , $

(13)

式中: T_max为该深度土壤的年最高温度, T_min为该深度土壤的年最低温度。相位φ₁和φ₂是由最小二乘法拟合式 (11), (12) 得到。比较z₁和z₂的振幅和相位, 由式 (10)~(12) 可以得到热扩散率K的振幅法和相位法计算公式:

$ {k_{\text{A}}} = \frac{{{\omega}{{\left({{z_1} - {z_2}} \right)}^2}}}{{2\ln \left({{A_1}/{A_2}} \right)}} , $

(14)

$ {k_{\text{p}}} = \frac{{{\omega}{{\left({{z_1} - {z_2}} \right)}^2}}}{{2{{\left({{\varphi _2} - {\varphi _1}} \right)}^2}}} , $

(15)

式中: K_A和K_P分别为振幅法和相位法求得的热扩散率。

3 青藏高原土壤温度和土壤热流计算结果分析 3.1 土壤温度年变化

土壤温度是计算土壤热通量的输入项, 因此在进行土壤热通量的计算前, 先了解不同深度土壤温度的差异。通过39个站点22年间0. 8 m和3.2 m深度处土壤平均温度随时间的变化曲线 (图 2) 中可以看出, 土壤温度表现为定常的正弦波, 深层土壤温度存在显著的年变化趋势。土壤温度随着深度的增加振幅衰减、位相延迟。从1980-2001年青藏高原东部站点土壤温度的月平均值 (图 3) 可以看出, 0. 8 m深度一年当中土壤温度最大值出现在7月份, 3.2 m土壤温度最大值要比0. 8 m延迟两个月左右, 出现在每年的9月, 说明年波从0. 8 m传到3.2 m大约需要两个月的时间, 这与范新岗和汤懋苍 (1994)的研究结果一致。

3.2 土壤热扩散率年变化

由于热扩散率受土壤质地、温度和湿度的影响, 在冬季和夏季会有差别, 为使得K更接近实际, 本文采用分段求法, 将一年的土壤温度资料分成冬夏两段拟合年波获得振幅和相位。根据上述结果 (图 3), 土壤温度年波从0. 8 m传到3.2 m大约需要两个月的时间, 因此本文将每年的资料分段如下: 0. 8 m以5-10月为夏季 (青藏高原夏季风时段), 以11月到次年4月为冬季 (青藏高原冬季风时段), 3.2 m以7-12月为夏季, 以1-6月为冬季。用振幅法和相位法及热传导对流法分别计算了1980-2001年青藏高原月平均土壤热扩散率 (图 4)。从图 4中可以看出用振幅法求出的土壤热扩散率K_A、相位法求出的K_P和热传导对流法求出K的大小排序分别为:振幅法 < 热传导对流法 < 相位法。缪育聪等 (2012)的研究认为, 振幅法和相位法分别低估和高估了土壤热扩散率, 热传导对流法拟合效果要好于振幅法和相位法计算的土壤热扩散率, 因此本文以下的研究将用热传导对流法求得的热扩散率计算土壤热通量。从图 4中可看出K值有明显的季节变化, 8月、11月和3月为相对高值期, 9月、1月和4月量值较小。土壤热扩散率是指单位时间, 在单位温度梯度作用下, 单位体积土壤温度的变化, 热扩散率是热传导率与热容量的函数, 由于土壤热容量与热传导率都是含水量的函数, 并且都随含水量的增加而增加, 但土壤热容量与含水量为线性关系, 而导热率与含水量呈现幂函数增加趋势, 两者增加趋势不同步, 因此热扩散率与含水量的关系基本为抛物线变化趋势。因此热扩散率开始随着含水量的增加而增加, 当含水量达到一定程度后, 反而随含水量的增加而减小 (邵明安等, 2006), 因此土壤热扩散率与土壤湿度的关系很密切。受高原季风的影响, 高原地区的降水量大部分都集中在5-8月左右, 在此期间活动层土壤水分增大, 土壤热扩散率随之增大。而在11月和3月左右, 土壤热扩散率的增加应该与土壤冻融过程和土壤湿度的变化有关, 由于没有土壤湿度观测资料, 这里不再深入分析, 但土壤热扩散率这种年变化趋势的原因值得进一步研究。用传导对流法计算的土壤热扩散率K的年际变化 (图 5) 中可以看出, 1980-2001年K值总体为下降趋势, 在80年代前期波动较大, 在1991-2001年有小幅度波动。这与张雯等 (2015)研究的青藏高原东部土壤湿度的年际变化特征一致。K的年际变化范围为3.4×10^-7~8. 9×10^-7 m²·s^-1。

3.3 土壤热通量年变化和年际变化分析 3.3.1 土壤热流年变化

从1980年1月至2001年12月5个站0. 8~3.2 m各月土壤热通量的年变化曲线 (图 6) 可以看出, 青藏高原北部的冷湖站土壤热通量的负极值出现在每年冬季1月前后, 正极值出现在每年夏季6月前后 (图 6a), 青藏高原东南部安多和泽当站土壤热流的逐月变化为先增大后减小, 土壤热通量的负极值出现在每年的12月左右, 正极值出现在7月份左右 (图 6b, c)。青藏高原东部刚察站 (图 6d) 和西宁站 (图 6e) 土壤热通量的逐月变化与南部一致。土壤热通量的年变化负极值比冷湖站稍有超前, 出现在每年的12月左右, 正极值出现在6月前后。总之, 青藏高原5个站三个地区土壤热流的年变化趋势类似。在夏季的6或7月土壤热通量达到正极值, 根据公式 (7) 热流向上为正, 说明深层土壤热流是从下往上传输。冬季12或者1月土壤热流出现负极值, 说明深层土壤热流是从上往下传输。根据张立杰等 (2005, 2006) 的研究, 出现这种现象的原因应该与土壤的不同冻融过程有关, 冻土的年变化过程分为秋季冻结过程、冬季降温过程、春季升温过程、夏季融化过程。在冬季降温期, 土壤内部热量渐渐流逝, 中层0. 8~2 m土壤热流从上往下传输, 2 m以下土壤热流微弱几乎为零, 所以在冬季0. 8 m以下土壤热流总体向下传输。在夏季融化期, 土壤深层存在较强的热交换, 0. 8~2 m土壤热流从下往上传输, 2 m以下从上往下传输, 土壤热通量非常微弱, 所以在夏季0. 8 m以下整体向上传输 (张立杰等, 2005)。另外, 张立杰等 (2006)对2~10 cm土壤热通量的研究发现, 由于夏季太阳辐射强烈, 地面接受净辐射较多, 使得地表温度高于地下土壤温度, 所以热量在表层为自上往下传输。冬季地表温度低于地下土壤温度, 热量为自下往上传输。说明0. 8~3.2 m土壤深层的热流变化与表层热流变化具有相反的相位。

由于海拔和降水量的差异, 土壤热通量的空间分布有所不同。青藏高原海拔高, 空气稀薄, 太阳辐射通过的大气光学路径较短, 因此太阳辐射被削弱的少, 太阳总辐射量要高于全国其他地区, 年总辐射量分布趋势自东南向西北增多, 藏东南地区为低值区, 藏北高原、阿里地区、柴达木盆地年总辐射量为高值区。冷湖位于柴达木盆地, 根据陈芳等的研究, 由于柴达木盆地西北端云雨天气少, 而晴朗天气多 (陈芳等, 2005), 冷湖的年总辐射量高达7269 MJ·m^-2, 由此向东向南, 随着云雨天气增多, 年总辐射量随之减少。根据李韧等的研究, 位于青藏高原中部的安多、泽当地区, 年总辐射量约为6400 MJ·m^-2(李韧等, 2012)。刚察和西宁位于青藏高原北部的祁连山地区, 祁连山地区的年总辐射量较柴达木盆地少, 普遍少于6200 MJ·m^-2。根据图 6a~e可以看出以冷湖为代表的柴达木盆地地区平均土壤热通量的最大值为2. 5 W·m^-2左右, 以刚察和西宁为代表的祁连山地区的平均土壤热通量最大值为2 W·m^-2左右, 而以安多和泽当为代表的青藏高原南部地区的平均土壤热通量最大值为1. 9 W·m^-2左右。说明土壤热通量的大小与总辐射大小密切相关。通过1980-2001年青藏高原39个站点0. 8~3.2 m土壤热通量平均年变化 (图 7) 可以看出, 传导热通量和对流热通量变化的相位相反, 总热通量与对流热通量变化具有相同的相位。另外从图 7中还可看出青藏高原土壤热通量在3月左右由负值转为正值, 说明深层土壤热量由获得热量向下传输转变为失去热量向上传输, 夏季土壤热通量也为正值, 热量向上传输, 而且在7月左右达到全年最大, 8月、9月土壤热通量值开始下降。且在秋季土壤热通量的值由正值转为负值, 即由向上传输改为向下传输, 这与张宏等研究的土壤表层的土壤热通量的变化相反, 但与前文研究的0. 8~3.2 m土壤深层的热通量变化与土壤表层热通量变化具有相反的相位一致。在冬季土壤热通量都为负值, 且在12月、1月左右向下传输的土壤热通量值达到最大, 随后减小。

另外, 从中东部空间分布图 (图 8) 可以看出, 最大值出现在青藏高原东北部约为0. 2 W·m^-2, 祁连山地区土壤热通量值相对较小, 沿柴达木盆往东南地区为高值区, 沿唐古拉山脉往东南地区土壤热通量非常微弱, 昆仑山脉以南地区为相对低值区约为0. 02 W·m^-2。

3.3.2 土壤热通量年际变化特征及其与高原季风的关系

从1980-2001年青藏高原平均土壤热通量年际变化特征 (图 9) 中可以看出, 青藏高原地区传导热通量在0线附近轻微变化, 最大值的绝对值在0. 1 W·m^-2左右, 最小值的绝对值在0. 003 W·m^-2左右, 量值较小。对流热通量与总的热通量具有相同的年变化趋势。对流热通量最大值在0. 6 W·m^-2左右, 最小值在0. 001 W·m^-2左右。另外, 将青藏高原地区年平均土壤热通量分别与白虎志等 (2001)、马振峰和高文良 (2003)研究的高原夏季风指数和汤懋苍等 (1984)定义的高原冬季风指数进行相关分析, 发现土壤热通量的年际变化曲线与高原冬季风的年际变化趋势相反, 与高原夏季风指数的年际变化一致 (图 10), 相关系数分别为-0. 53和0. 58, 都通过了95%的显著性检验。说明深层土壤热通量的变化与高原季风的变化有密切的关系。高原季风和土壤热通量呈这种变化关系的原因还有待进一步研究。

4 结论

(1) 青藏高原深层土壤温度随深度的增加振幅减小、位相延迟, 0. 8 m土壤温度最大值出现在7月左右, 3.2 m土壤温度最大值出现在9月左右, 比0. 8 m的延迟两个月左右, 说明年波从0. 8 m传到3.2 m大约需要两个月的时间。

(2) 土壤热扩散率受温度和湿度的影响, 积雪的覆盖和降水量的增加都会导致土壤扩散率的增加。k的变化范围在3.4×10^-7 ~8. 9×10^-7 m²·s^-1。

(3) 传导热通量与对流热通量的年变化具有相反的相位, 总热通量与对流热通量的年变化具有相同的相位。对流热通量与传导热通量具有相同的量级, 均为10^-2~10⁰ W·m^-2。而且由于土壤不同冻融过程的影响, 土壤热通量在春季由负值转为正值, 夏季土壤热通量都为正值, 秋季由正值转为负值, 在冬季都为负值。土壤热通量有显著的空间分布特征:在祁连山地区土壤热通量值相对较小, 沿柴达木盆往东南地区为高值区, 沿唐古拉山脉往东南地区土壤热通量非常微弱, 昆仑山脉以南地区为相对低值区。

(4) 深层土壤热通量的年际变化非常小, 热通量的量级为10^-3~10^-1 W·m^-2。另外土壤热通量的年际变化与高原冬季风指数变化趋势相反, 与高原夏季风指数趋势一致, 相关系数分别为-0. 53和0. 58, 都通过了95%的显著性检验。