引言

当前人类社会面临的诸多生态和环境问题, 吸引了众多科学家从事于地圈-生物圈-大气圈之间相互作用问题的研究, 其中气候系统与陆地生态系统之间的相互作用尤为世人关注.这一方面是由于气候变化直接关系到人类的生存环境, 另一方面是因为陆地生态系统是人类社会各种生活必需品最基本的来源地.本文依据我们近年在黄淮平原中心区大面积农林复合系统中进行的边界层低空探测资料, 设计了一个三维U中尺度原始方程模式, 模拟此类生态系统与低层大气间的相互作用, 以求在更高分辨率的基础上揭示大气与生态系统的相互作用机理.

农林复合是指在有限的土地资源上寻求农林作物的最佳配制, 最大限度的利用光、热、水、肥、土等自然资源, 使单位土地面积上的生产力达到最高.农林复合实践在我国已有上千年的历史, 有关农林复合的研究无论是深度或广度, 我国均处于世界领先水平.黄淮海平原绵延数百公里规格化的农田林网更是世界上独一无二的.研究这一独具中国特色的生态系统对低层大气环流的影响及其有关的反馈机制, 对于目前普遍关心的气候与生态系统的双向耦合问题有其特别意义.另一方面, 该区农林植物的规格化配制和平坦的下垫面特征也为数值模式的设计和植被参数化方法提供了便利条件.

1 实验区概况及边界层低空探测

实验区位于黄淮海平原南部中心区 (33.7°N, 116.6°E) 安徽省涡阳县境内.边界层内的低空探测于1992年4月21~24日进行.当时构成林网的主要树种杨树已放半叶, 林下小麦正值灌浆抽穗期.平均株高度65 cm, 肥水适中, 生长良好.冠层叶面积指数4.6, 林带在区域上平均叶面积指数为1.0.农林复合系统群体叶面积指数为5.6, 属全郁闭绿色冠层.

除了常规的微气象观测以外, 在林网中心设立简易观测塔架 (图 1), 以加大冠层内部的观测密度, 塔上梯度观测共设4层:即2.0 m, 4.0 m, 5.8 m和8.5 m.改进的低空探测仪可以探测不同高度的风向、风速、温度、湿度和气压分布.要素记录的垂直分辨率为15 m, 低探每3小时释放一次持续72小时 (从4月21日14:00到4月24日11:00, 地方时, 下同).探空曲线的初步分析结果及大气稳定性特征已有专门论述^[1].

2 模式简介 2.1 基本方程

包括地形起伏的三维热力-动力学方程组可以写成如下形式:

(1)

(2)

(3)

(4)

(5)

(6)

除了各种常规气象学符号外, z为地形追随坐标中的高度:

(7)

c=C_p (p/p₀) ^R/C_p为表示气压的Exner函数, p₀=1000 hPa, w为 (x, y, z) 坐标系中的垂直速度, 且有

(8)

z_g=z_g (x, y) 为二维地形起伏, H为模式顶高度 (H=1000 m), θ_v为虚位温.方程 (3) 中的Q^*为净辐射, L为蒸发潜热, E_t是大气中水汽的变率.在农林复合集区, 大气中水汽的增加包括地面蒸发和植被蒸腾两部分, 在无植被区, 则主要由土壤蒸发引起, 单位为[kg·m^-3·s^-1]; 而大气中水汽的减少主要由凝结所致, 在模式高度范围内凝结量采用饱和凝结法计算.使用不同下标的F分别代表湍流扩散作用于风速 (F_u, F_v)、位温 (F_θ) 和空气比湿 (F_q) 的贡献, 它可用下面的通式表示

(9)

森林植被的动力学摩擦对于风速的减弱由下式计算

(10)

式中C_d为拖曳系数, 它在农林复合植被区的特殊形式根据Van E K和Draaigers^[2]的方案求得, L_ad为植被叶面积密度垂直分布函数, 农林复合系统中各层的叶面积密度数值见表 1.方程 (1)~(10) 通过下面的地表能量平衡方程、植被蒸发散方程以及湍流扰动动能变化方程, 组成闭合方程组.

2.2 地表能量平衡方程和植物蒸发散

地表能量平衡过程同时考虑太阳辐射穿过植被冠层时的减弱, 作物冠层对地面和大气长波辐射的阻挡, 土壤热流入量及土壤热学特性等地表过程.具体计算方案与早期的生态边界层二维模式相同^{[3, 4]}.

在植被区上空, 植被的蒸腾作用对于大气中水汽含量的贡献是不能忽略的, 方程 (3)、(4) 中的E_t即是反映了这一项的作用, 根据经典的Penman-Monteith公式, 某一层单位面积内的蒸发、蒸腾速率为

(11)

式中S为饱和水汽压曲线在空气温度T_a处的斜率, γ=0.67 hPa/k, R_nl为某层植被所截获的净辐射, e_s、e_a分别为温度T_a所对应的饱和水汽压和实际水汽压, r_b为叶面边界层空气动力学阻力, 取为60 s/m, r_s为叶片气孔的阻力数^[4]由下式求算^[5]:

(12)

式中Q_s为达于植被冠层上的太阳短波净辐射, 最小气孔阻力系数r_{s, min}和最大气孔阻力系数r_{s, max}的数值见表 1.空旷区土壤的蒸发假定为同时刻植被区蒸腾量的15%^{[4, 5]}.

2.3 湍流扰动动能

湍流扰动动能 (TKE) 是度量边界层湍流特征的综合性指标, 在农林复合系统上空, 植被的动力学摩擦消耗一定量的空气运动动能, 并通过改变湍流运动的混合长影响边界层湍流结构.湍流扰动动能一般由下式求得

(13)

式中, 方程右端前三项的动力学意义是清楚的, 第四项是浮力生成项, 最后一项是湍流运动的分子耗散项, 其中l是混合长, 它在农林复合系统中特殊形式已作过详细讨论^[6].本文仍沿用文献[6]中提供的算法.其中湍流交换系数K_H、K_Z、K_E和K_θ与桑建国等介绍的计算方法类同^[7].

2.4 薄转换层

中午前后强烈的太阳辐射可使无植被地区的地面温度迅速升高, 造成超绝热现象, 这是根据地面能量平衡方程计算地面温度时经常遇到的难题, 也是保持数值积分稳定的主要障碍.Salerno等建议引进一薄转换层 (TTL), 用TTL高度上的温度作为模式的下边界温度, 它与地面温度T_G的关系由下式给出^[6]:

(14)

式中L₀为薄转换层的高度, T_d为土壤中深度为d处的土壤温度, 即温度日变化近似为0处的温度.w/k为地面附近大气平均的垂直速度和涡度垂直扩散系数之比.根据文献[4, 8]和我们的试验, 取w/k=1 m^-1, d=0.85 m, 在无植被区白天取L₀=0.25 m, 晚上取L₀=0.16 m, 而农林复合区取L₀=0.01 m.

2.5 模式网格设计

为了反映各要素的垂直方向上梯度分布的不均匀性, 垂直方向采用可变步长, 垂直格距随高度的变化见文献[3]和[4].水平方向则采用均匀步长, △x=△y=2000 m.模式区域分为植被区和无植被区两部分, 在x方向从第11个格点开始向右无限延伸为植被区, 之前为无植被区, 模式中所用的不同地表特征参数如表 1所列.这里的无植被区假设完全是从模式设计的通用性来考虑的, 据当地林业部门统计, 实验区所在县 (市) 的林网覆被率已达96%以上, 加上该区为我国十大商品粮生产基地, 农作物覆盖率更高, 实际的完全无植被区是不存在的.相应的差分方案和初边值条件设定与早期的二维模式相同^{[3, 4]}.

3 结果讨论 3.1 地面温度日变化

为了集中反映农林复合系统热力和动力学效应, 观测和模拟都是对晴好天气进行的.模式以1992年4月22日05:00的观测场为初始场开始积分.数值模拟显示, 该模式保持积分稳定的时间步长不大于30 s.图 2a是模式输出的地面温度变化曲线.为了在小麦田中放置地面温度表, 观测时拔除了2 m×2 m的小麦植株而形成了小片田间空地, 因此地面温度的观测值更接近无植被区的情形.比较合理的做法是保持植被结构不受破坏, 用纤细的热敏电阻测定植被区的地面温度^[11].模式输出的农林植被区的温度曲线振幅明显小于观测值和无植被区, 无植被区温度日较差最大可达28 ℃, 而农林植被区为14.5 ℃.考虑植被冠层对于辐射的阻挡减弱, 上述结论是合理的.由于植被对于太阳辐射的减弱作用, 植被区地面温度白天升温较慢.最大值出现较无植被区略迟.夜间植被冠层阻挡了大量地面长波辐射, 使得植被区地面温度高于无植被区.这也是农林复合系统中冬小麦冬季少受害的主要原因^[12].

至冠层高度以上15 m处, 模式输出值与观测结果的一致性有明显改善 (图 2b).植被区与无植被区的差异也明显减小, 但二者变化的位相仍显示差异.日出后无植被区温度升高较快, 最大值略高于农林复合区.夜间, 由于植被冠层的强烈辐射冷却, 植被区最低温度低于空旷区.在无植被空旷地区, 15 m高度已处于夜间辐射逆温层内.这说明植被的存在使辐射冷却层有所抬高.本文研究的农林复合系统在观测时的平均高度为12 m左右.辐射冷却层的抬升在20 m左右 (见图 3a和图 3b).某些高大的森林植被可达30~40 m, 它们对于辐射冷却层的抬升应引起模式设计者的重视.

3.2 大气温度的时空变化

图 3a和3b分别表示无植被空旷区和农林复合区温度的时空剖面图, 二者在200 m以上差异很小, 为了突出其间差异, 作图高度限制在300 m以下.图中示出二者的时空变化基本一致, 主要差别在于:无植被区 (图 3a) 日较差最大值出现在地面, 并随高度增加而减小.植被区 (图 3 b) 日较差最大值出现在地面以上10 m左右, 这一高度位于小麦冠层以上, 树冠高度以下, 是湍流活动最微弱的高度.农林作物冠层夜间的辐射冷却使得植被区最低温度出现于较高的层次, 针对本文研究的黄淮海平原地区的农业复合系统, 辐射冷却层在20 m左右.应该指出, 贴地层中热力特性的这些差异, 除非象本文使用的细密网格是很难捕捉的, 它对于大尺度的天气系统模拟可能是微不足道的, 但在研究生物生长的微环境问题时至关重要.近年来颇为流行的作物产量生态学数值模拟显示, 生长环境温度的微小变化, 可导致作物生育期、作物终极生产力, 乃致整个生态系统的CO₂同化功能的较大变异^[13].这类工作是生态边界层研究对经典气象学模式的延伸.

3.3 植物蒸发散与低层大气湿度

生态边界层模式与一般大气边界层模式的根本区别在于模式需重点考虑生态系统中农林作物的生命活动, 蒸腾作用是作物最重要的生命活动之一, 也是影响冠层上方空气湿度的主要因子.图 4a是根据方程 (11)、(12) 计算的农林复合系统上方植物的蒸发散量和可能发生的凝结量.由图可见, 植被区每天的蒸发量 (曲线1) 约4 mm, 这一数值与同季节在同类植被区的农业气象观测记录^{[11, 14]}是吻合的.相应的凝结 (露) 量 (曲线4) 约为每天0.8 mm.这与当时的田间露量观测十分一致.曲线2是裸地蒸发量计算值, 模式中假定它为植被区的15%.每天约为0.6 mm左右, 相应的凝结 (露) 量 (曲线3) 几乎为零, 这是裸地经常观测到的事实.本文对于裸地蒸发量的计算方法比较粗糙, 较精确的裸地蒸发量计算需要土壤水势或土壤水含量的观测数据^[12].但如前所述, 模式中的无植被区在黄淮海平原是很难找到的, 为了模式设计的通用性, 假设气流从无植被空旷区流向农林复合区.当模式区域范围内包含真实的裸地时, 需用比较精确的土壤蒸发量计算方案.图 4b是地面以上TTL高度处的空气湿度的模拟结果, 植被区和无植被区均显示出比较明确的日变化特点, 与观测结果的一致性也是令人满意的.图中示出, 农林复合区在白天的空气湿度一般大于无植被的空旷区.夜间的数值比较接近, 清晨以后, 空旷区的湍流活动发展较快, 土壤蒸发的水汽迅速在近地层扩散, 造成空旷区湿度较大, 植被区由于射达地面的短波辐射较少, 地面蒸发很小, 而植被顶层蒸腾较大, 且湍流发展滞后, 故TTL处空气湿度反而较空旷区低.田间观测证实, 至上午10:00麦田露水仍不能完全蒸发.之后, 随着蒸发蒸腾加剧, 植被区空气湿度迅速增加, 明显大于空旷区.至傍晚时分, 植被区湍流交换再次减弱, 地面蒸发很小, 使得此时的TTL处湿度再次小于空旷区的, 但在冠层顶附近及以上 (见图 5a和5b), 植被区空气湿度总大于空旷区, 这也是林网保护区小麦生长免受干热风危害, 产量相应提高的主要原因之一^{[12, 13]}.

模式输出的空气湿度时空变化特征如图 5所示.植被区湿度分别在上午10:00和傍晚出现两次极值, 这是模拟植被蒸腾的气孔模型经常揭示的事实.其原因是上午10:00前后及傍晚时分各种生态环境因子如光照、温度、湿度等的组合为最佳状态.此时气孔开张度最大, 而中午前后, 光照太强, 温、湿条件也不是最适, 作物生理活动出现“午休”现象.故此时蒸腾作用趋缓, 近地层湿度出现相对低值, 这一现象已为多数植物蒸腾的田间观测和各种生理生态学模式所证实^[10]

3.4 近地层湍流结构特征

湍流脉动动能 (TKE) 是综合反映边界层湍流发展强度的量, 图 6a和6b是模式输出的TKE时空剖面图.与无植被空旷区 (图 6a) 相比, 由于植被冠层摩擦造成的动能耗散, 农林复合区 (图 6b) TKE的绝对数值略小.但二者的日变化特征和最大值出现的高度都比较一致, 这说明在平坦的黄淮海平原上, 农林复合植被的动力学效应不足以改变大气流场的基本特征, 但它对于近地层 (300 m以下) 温、湿特性的影响是不容忽视的.它改变了农林作物生长的微环境, 进而影响作物的生长发育.在研究气候与生态系统的相互作用问题时, 已开始注意到所谓双向耦合过程.即:既考虑地面植被对大气流场的影响, 也考虑变化了的近地层环境对陆地生态系统反馈作用.这正是生态边界层模式希望解决的问题.

模式模拟的TKE日变化特点与边界层内的各种观测结果^[10]比较一致, 由于没有外来天气系统的影响, 近地层湍流发展主要为辐射加热所驱动, 日出后, 湍流开始发展加强, 至午后达最大值.日落后, 地面长波辐射增强, 逆温层逐渐形成, 大气层结趋于稳定, 湍流难以发展, TKE数值接近零值.