长白山阔叶红松林的零平面位移和粗糙度

长白山阔叶红松林的零平面位移和粗糙度
3
赵晓松
1,233
　关德新1　吴家兵1　金昌杰1　韩士杰
1
(1中国科学院沈阳应用生态研究所,沈阳110016;2
中国科学院研究生院,北京100039)
摘　要　根据长白山阔叶红松林气象观测塔上16个月的风速、温度、湿度及气压的连续观测
资料,根据中性层结条件下风速随高度的对数变化规律,利用廓线法中的牛顿迭代法计算了该森林的零平面位移d 和粗糙度z 0,结果表明,d 和z 0均存在着较明显的季节变化,在生长季d 较大,z 0较小,而非生长季恰好相反,标准化的零平面位移d/h 和粗糙度z 0/h 在生长季和非生长季平均分别为01867,01764和4147×10-2,3159×10-2。

与叶面积指数对比分析发现,d/h 和z 0/h 与叶面积指数分别存在正相关和负相关的关系。

敏感性分析表明,牛顿迭代法求d 和z 0对风速的精度要求较高,需要精度高的仪器和长期的数据积累。

关键词　零平面位移,粗糙度,阔叶红松林,牛顿迭代法
中图分类号　S716 文献标识码　A 文章编号　1000-4890(2004)05-0084-05
Z ero 2plane displacement and roughness length of the mixed forest of broad 2leaved and K orean 2pine in Changb ai Mountain.ZHAO Xiaosong 1,2,GUAN Dexin 1,WU Jiabing 1,J IN Changjie 1,HAN Shi 2jie 1(1Institute of A pplied Ecology ,Chinese Academy of Sciences ,S henyang 110016,China ;2
Graduate School of Chinese Academy of Sciences ,Beijing 100039,China ).Chinese Journal of E 2cology ,2004,23(5):84～88.
Based on logarithm law of wind profile under neutral stratification ,zero 2plane displacement d and roughness length z 0of the mixed forest of broad 2leaved and K orean 2pine in Changbai Mountain were calculated with Newton iteration according to the 16months profile measurements of wind speed ,air temperature ,humidity and pressure.The results showed that d and z 0had obvious seasonal variation.In growing season was relatively high and z 0was lower.Normalized values ,d/h and z 0/h ,were in the magnitude of 01867,01764and 4147×10-2,3159×10-2in growing season and dormant sea 2son ,respectively.d/h and z 0/h were positively and negatively correlated with leaf area indexes of the forest.Sensitive analysis showed that the iteration approach of calculating d and z 0demanded high precision measurements of wind speed in the profiles.This work should be based on the long term measurements with high 2precision instruments.
K ey w ords zero 2plane displacement ,roughness length ,mixed forest of broad 2leaved and K orean 2pine ,Newton iteration approach.
3国家自然科学基金项目(30370293)和中国科学院知识创新工程重大资助项目(KZCX12SW 201201A1)。

33通讯作者
收稿日期:2004-02-27 改回日期:2004-07-01
1　引　言
零平面位移(d )和粗糙度(z 0)是植被动量、热量
和气体交换等参数化过程的重要因子,用来表达这些量的湍流交换系数[17]。

随着涡动相关技术越来越广泛地应用于森林等高大植被的CO 2、水汽等气体交换的测量,作为空气动力学分析中的重要因子,同时也是footprint 等多种模型的重要输入参数,d 和z 0得到了很大的关注。

求解d 和z 0的方法主要有廓线法[6,9,23]、解析分析[3,12,13,14,19]、数值模拟[20,22]和湍流综合法[1,5]等方法。

廓线法的优点是反映了流场实际观测结果,但需要高精度、运行稳定的自动观测仪器;数值模拟方法的流体力学机理明确,但需要野外实际观测的验证;湍流综合法建立
在湍流统计关系的基础上,并需要超声风温仪等精
密仪器,应用较少。

获取观测数据的方法主要有风洞模拟实验和野外实地观测,由于风洞模拟实验具有模型可随意改变的优点,所以很多学者进行了不同类型模型的风洞实验研究[2,4,7,8,16,24],但风洞实验中的模型很难实现与野外实际植被之间严格的几何相似(如枝叶尺度及其分布)和物理相似(如层结稳定度、枝杆的力学性质等),所以只能模拟特定条件的情况,与野外实际观测结果相互补充是十分必要的。

因为d
生态学杂志Chinese Journal of Ecology 2004,23(5):84～88
和z0是植被与大气之间动力、热力作用研究的基本参数,分析其长期变化(如季节变化)有重要的意义,特别是落叶林的叶面积季节变化很大,d和z0也应有明显的改变,Shaw(1982)的数值模拟结果证明了这一点,但至今还没有实际观测的研究结果来完整地描述落叶森林d和z0的季节变化过程,绝大部分是生长季某一时段的研究[1,5,10,11,18,21,23]。

本文拟利用先进的自动观测仪器获得的较长时间的观测数据(16个月),用廓线方法中的牛顿迭代法计算长白山阔叶红松林的d和z0,分析其季节动态变化规律以及与叶面积指数(LAI)的关系,为森林-大气物质、能量交换的研究提供基本参数。

2　研究地区与方法
211　自然概况
本项研究在中国科学院长白山森林生态系统定位站一号标准地阔叶红松林内进行(42°24′N,128°6′E,海拔738m)。

年平均气温为316℃,年平均降雨量713mm。

土壤为山地暗棕色森林土。

观测场附近下垫面地势平坦,林型为成熟原始林,主乔木为红松(Pi nus koraiensis)、椴树(Tilia am urensis)、蒙古栎(Q uercus mongolica)、水曲柳(Fraxi nus m and2 shurica)和色木(A cer m i mo)。

林分为复层结构,下木覆盖度40%,平均株高26m,立木株数约560株・hm-2。

212　方法
21211　观测方法　观测林地建有高62m的微气象观测塔,在观测塔上安装了7层气象要素探头,高度分别为215,8,22,26,32,50和60m,测量要素有风速(A100R,Vector Instruments,Denbighshire,U K)、大气温度和湿度(HMP45C,Vaisala,Helsinki,Fin2 land)。

另外,215和60m分别测量气压(CS105, Vaisala,Helsinki,Finland)和风向(W200P,Vector Instruments,Denbighshire,U K),而40m的风速用三维超声风速仪(CAST3,Campbell,USA)测量。

原始采样频率为2Hz,通过数据采集器(CR23X and CR10X,Campbell,U T,USA)采集并按30min计算平均值进行存储。

观测时间从2003年8月24日～2004年12月31日。

同步进行叶面积指数LAI的测量(Li2000, LiCor Inc1,Lincoln,N E,USA),时间从2004年4月～10月末。

21212　数据处理　在中性条件下,风速廓线表达式为
U(z)=
u3
k
ln
z-d
z0
(1)式中,U(z)为z高度上的平均风速,u0为摩擦风速,k为von Karman常数,取014,z为测量高度,d 为零平面位移,z0为粗糙度。

用梯度理查孙数(R i)来判断近地面层的大气稳定度,
Ri=
g
θv
5 θv
5z
5 U
5z
2
(2)式中,Ri为梯度理查孙数,g为重力加速度,取10m ・s-2,θv为虚位温,根据公式θv=θ・(1+0161r)计算得到,其中r为未饱和空气混合比(g・g-1),θ为位温(θ=T(p0/p)01286,T为大气温度(K),p为气压(kPa),p0为基准气压,取100kPa)。

根据方程(2)的差分形式,选择26和50m两个高度差计算每一时刻的R i值,选取R i值区间(-0103,0103)作为中性层结[15]。

根据方程(1),应用牛顿迭代法计算出d。

具体算法如下:
设z1,z2,z3三个测量高度的风速分别为u1, u2,u3,则首先确定迭代公式为
f(d)=
u1-u2
u1-u3
=
ln(z1-d)-ln(z2-d)
ln(z1-d)-ln(z3-d)
(3)和迭代函数g(d)=d-f
(d)
f′(d)
(4)给d的初值,求得g(d),再令d=g(d),代入迭代函数求g(d),反复迭代,直到相邻两次g(d)之间的差值小于某个给定的误差(本文取01001),则g(d)即为所求的d。

把d代入(1)中,则(1)变为l n(z-d)为自变量、U(z)为因变量的一元线性方程,用线性回归方法即可求出z0。

本文选择[z1,z2,z3]=[26,32,50m]的风速梯度计算d和z0。

为了减小计算结果的波动性,取5 d的平均值进行时间系列分析。

3　结果与分析
311　层结状态分布
理查孙数R i为负表示静力不稳定层结,多出现在白天,R i为正表示静力稳定层结,多出现在夜间, R i等于零为中性层结,一般出现在日出前或日落后一段时间内。

根据方程(2)计算出16个月梯度理查
58
赵晓松等:长白山阔叶红松林的零平面位移和粗糙度
孙数(R i )的分布如图1所示,负值的比例略大于正
值,R i 在(-100,0)范围占50171%,R i 在(0,100)占48121%,说明从16个月的累积情况来看,不稳定层结略大于稳定层结。

本文所选取的中性稳定层结(R i (-0103,0103))占20135%。

图1　观测期间的理查孙数分布Fig 11　Distribution of Rich ardson number during the measurement
312　零平面位移d 和粗糙度z 0的季节变化
在中性层结条件下进行迭代计算,得到5d 平均的标准化零平面位移(d/h )和粗糙度(z 0/h )时间系列(图2,3)(其中h =26m 为林分平均高度)。

从图2可以看出,相对于z 0/h ,d/h 值的波动较小,具有明显的季节变化,表现为生长季高,非生长季低的趋势。

按时间进程大致可分为4个阶段:①从观测开始的2002年8月24日至9月下旬的生长季,d/h 值较高,
平均为01858;②2002年10月到2003年5月中旬,为跨年度的非生长季,占林分较大比例的阔叶树落叶,d/h 随之减小,d/h 的平均值为01746;③2003年5月下旬到9月中旬,为2003年生长旺季,叶片密度迅速增加,阻碍了空气进入冠层,使得d/h 值增大保持较高值,平均为01884。

④2003年9月中旬到年末,d/h 又下降到01782。

图2　5d 平均的d/h 的季节变化Fig 12　Seasonal variation of 52d ay averaged d/h
z 0/h 值的波动相对较大(图3),从趋势线(周期为5的滑动平均)的变化看,在2003年能够观察到一定的季节变化,表现为生长季低,非生长季高的趋势,而在2002年观测的短暂的生长季中这种趋势并未体现出来。

2002年8月底到2003年5月底之间z 0/h 的平均值为4158×10-2,2003年6～9月的生长季,叶片密度大,森林郁闭度较高,从而使森林植被冠层的粗糙度降低,z 0/h 的平均值减小为3159×10-2,随着生长季的结束又上升到4136×10-2。

图3　5d 平均的z 0/h 的季节变化Fig 13　Seasonal variation of 52d ay averaged z 0/h
z 0/h 和1-d/h 的关系(图4),为线性关系,斜率为01093,截距为01025,r 2=01344。

如果强制使直线经过原点,即使其满足z 0/h =λ(1-d/h )方程,则λ=0121。

图4　z 0/h 和12d/h 的关系
Fig 14　Plot of norm alized roughness length z 0/h against 12d/h
313　d/h 和z 0/h 与叶面积指数(LAI )的相关关系
标准化的零平面位移d/h 和粗糙度z 0/h 与LAI 的关系(图5)。

d/h 与LAI 呈很好的正相关关
系,即d/h 随着LAI 的增加而升高,r 2=01628;而
z 0/h 与LAI 呈负相关关系,但相关系数较低,r 2
=
6
8 生态学杂志　第23卷　第5期　
01306,这与z 0/h 波动大有关。

图5　标准化的零平面位移d/h 和粗糙度z 0/h 与叶面积指数的相
关关系
Fig 15　R elationship bet w een norm alized displacement height (d/h ),
norm alized roughness length (z 0/h )and leaf area index
Shaw 等[20]认为零平面位移随LAI 单调递增,而粗糙度与LAI 呈单峰曲线,即当LAI 很小时,z 0/
h 随LAI 增大而增大,z 0/h 达到峰值后,又随LAI
增大而减小,z 0/h 取得峰值时LAI 的大小与冠层的垂直结构有关,Z max /h (Z max 为植被LAI 最大处的高度)越大,取得峰值时的LAI 越小。

图5(a )d/h 与LAI 的关系满足上述的单调递增的关系,而图5(b )中z 0/h 与LAI 只表现为随LAI 单调递减,这是因为长白山阔叶红松林为针阔混交林,即使在非生长季由于有针叶树的存在,LAI 也会很大,在冬季LAI 最小时仍>119,远超过了
z 0/h 达到峰值时的LAI ,因此,只表现出递减趋势。

314　牛顿迭代法的敏感性
本文采用的牛顿迭代法(即公式(2)(3)(4))可以利用不同的三个观测高度[z 1,z 2,z 3]。

[z 1,z 2,z 3]取值不同,则上述方法求出合理d 值对应的a =
(u 1-u 2)/(u 1-u 3)的取值范围不同(图6)。

[z 1,z 2,z 3]=[26,32,50m ]时,a 的取值范围为(0132,
016),[z 1,z 2,z 3]=[32,50,60m ]和[32,40,50m ]
时,a 的取值范围分别为(0150,0165)和(0171,0180)。

后两种情况下a 的范围太小,计算出的合理
d 值的数量(<30个)不具有统计学意义,
因此,本
文选择[26,32,50m ]的风速梯度计算d 和z 0。

图6　观测高度[z 1,z 2,z 3]不同时的零平面位移d 随a =(u 1-u 2)/(u 1-u 3)的变化
Fig 16　V ariation of displacement height (d)with a =(u 1-u 2)/(u 1
-u 3)by different observation heights[z 1,z 2,z 3]
表1　廓线法26m 风速变化对d/h 和z 0/h 计算结果的敏感性(计算值及变化量百分比)T ab 11　Sensitive of calculation for d/h and z 0/h to wind velocity at 26m height
风速提高比例
(%)4个时间段
Ⅰ(200218～9)
Ⅱ(2002110～200315)Ⅲ(200315～9)Ⅳ(2003110～12)215d/h 018485(-1110%)017029(-5177%)
018690(-1169%)017492(-4118%)z 0/h 010606(32134%)010438(25112%)010507(16134%)5d/h 018379(-2134%)016426(-13186%)018607(-2163%)017072(-9156%)z 0/h 010736(60167%)010466(33114%)010588(34187%)715d/h 018250(-3185%)015587(-25111%)018522(-3158%)016489(-17102%)z 0/h 011037(12615%)010493(40190%)010756(73141%)10
d/h 018089(-5172%)014399(-41102%)
018427(-4166%)015669(-27150%)z 0/h
011211(16415%)
010531(51183%)011012(13212%)
按粗糙表面上的大气边界层理论[23],用廓线法计算d 和z 0时,26m 观测高度显得偏低,这可能使
d 和z 0的计算结果产生一定的误差。

由于26m 探
头距林冠很近,使测得的风速偏小,进而可能导致计算的d 值偏大而z 0值偏小。

为了估计其影响,将26m 的风速分别提高215%,5%,715%,10%,结果
7
8赵晓松等:长白山阔叶红松林的零平面位移和粗糙度
如表1所示,非生长季的d/h值有所降低,而生长季的d值降低不明显。

当26m风速提高215%时,上述4个时间段d/h值分别降低了1110%, 5177%,1169%和4118%,提高5%时,分别降低2134%,13186%,2163%和9156%,提高715%和10%时,d/h值变化较大,趋于不合理。

z0/h的变化也表现为非生长季变化幅度大,而生长季变化相对较小的特点。

当26m风速提高215%时,z0/h 值在上述3个时间段分别升高了32134%,25112%和16134%,提高5%时,分别为60167%,33114%和34187%,提高715%和10%时,数值变化较大,趋于不合理。

根据分析,认为风速提高5%的计算结果比较合理,即4个时间段d/h值分别为01837, 01642,01861和01707,而3个时间段z0/h值分别为7136×10-2,4166×10-2,5188×10-2。

4　讨　论
牛顿迭代法计算结果表明,长白山阔叶红松林标准化的零平面位移d/h和粗糙度z0/h在2002年8月～2003年12月内划分的时间段内分别为01858,01746,01884,01782和4158×10-2,3159×10-2,4136×10-2。

d/h和z0/h与叶面积指数分别存在正相关和反相关的关系。

如果考虑26m观测风速偏低而进行5%的订正,d/h和z0/h的取值在四个时间段内分别为01837,01642,01861, 01707,和7136×10-2,4166×10-2,5188×10-2是比较合理的。

我们的计算结果与国内外的同类研究进行了比较,其中文献引用率最高的是Shaw等[20]的数据模拟结果,根据他描述LAI在119～6范围内且Z max/ h>016时(由于长白山阔叶红松林中多为成熟乔木,因此,估算Z max/h>016),d/h大约在016～0186,z0/h约在0104～0109。

本文计算的d/h (01746～01884)略大于Shaw等[20]给出的范围,而z0/h(010359～010458)略小于其范围。

如果对26 m风速进行5%订正后,d/h(01642～01861)和z0/ h(010466～010736)均在Shaw等[20]所描述的范围内,这也说明了对26m风速进行5%的订正比较合理。

刘和平等[1]根据1995年8月25日～9月6日测量数据,利用湍流综合法估算的同一森林d= 1915m,z0=116m,标准化后d/h=0175,z0/h= 010615,d/h显著低于本文2003年生长季的01884,而z0/h显著高于010436,即使去除26m风速的影响,本文计算的d/h仍明显高于0175,这可能是由于不同的计算方法引起的。

与本文生长季的结果接近的有针叶林(当LAI =415时,d/h=0178)[10],玉米地(当LAI=415时,d/h=018)[18]。

到目前为止,作者尚未发现以实际观测资料分析d和z0长期变化的研究报道,主要原因是廓线法对观测数据敏感性太强的缘故,所以在分析季节变化时,都引用Shaw等[20]的数据模拟结果。

本文在方法上进行了尝试,希望突破以往有关d和z0取值的经验性约束,为森林2大气之间物质能量交换的研究提供实测的参数。

参考文献
[1]　刘和平,刘树华,朱廷曜,等.1997.森林冠层空气动力学参数
的确定[J].北京大学学报(自然科学版),33(4):522～528. [2]　关德新,朱廷曜,金昌杰.1999.中性温度层结时森林郁闭度与
林冠上风速廓线的关系[J].北京林业大学学报,19(1):26～
30.
[3]　Businger JA.1974.Aerodynamics of vegetated surfaces[M].
Scripta,Wash.DC:Proc.Sem.on Heat and Mass Transfer in the
Environment of Vegetation(Dubrovnik,Yugoslavia),139～165.
[4]　Counihan J.1971.Wind tunnel determination of the roughness
length as a function of the fetch and the roughness density of
three2dimensional roughness elements[J].A t mos.Envi ron.,5:
637～642.
[5]　De Bruin HAR,Verhoef A.1997.A new method to determine
the zero2plane displacement[J].Boundary2L ayer Meteorol.,82:
159～164.
[6]　Dolman AJ.1986.Estimates of roughness length and zero plane
displacement for a foliated and non2foliated oak canopy[J].A2
gric.For.Meteorol.,36:241～248.
[7]　Hussain M,Lee BE.1980.A wind tunnel study of the mean pres2
sure forces acting on large groups of low2rise buildings[J].J.
Wi nd Engi n.Indust r.Aerodyn.,6:207～225.
[8]　Iqbal M,Khatry A,Seguin B.1977.A study of the roughness ef2
fects of multiple windbreaks[J].Boundary2L ayer Meteorol.,
11:187～203.
[9]　Jackson PS.1981.On the displacement height in the logarithmic
velocity profile[J].J.Fl ui d.Mech.,111:15～25.
[10]　Jarvis PG,James G B,Landsberg JJ.1976.Coniferous forest[A].
In:Monteith JL,ed.Vegetation and the Atmosphere.Vol.2
[C].London:Academic Press,171～240.
[11]　K imura R,Otsuki K,Kamichika M.1999.Relationships between
the zero2plane displacement and the roughness length over
sorghum and alfalfa canopies[J].J.A gric.Meteorol.,55:15～
24.
[12]　K ondo J.1971.Relationship between the roughness coefficient and
other parameters[J].J.Meteorol.Soc.Jap.,49:121～124. [13]　Leonard RE,Federer CA.1973.Estimated and measured rough2
ness parameters for a pine forest[J].J.A ppl.Meteorol.,12:
302～307.
[14]　Lettau HH.1969.Note on aerodynamic roughness2parameter es2
timation on the basis of roughness2element description[J].J.
A ppl.Meteorol.,8:828～832.(下转第115页)
88
生态学杂志　第23卷　第5期　
port capacity of chloroplasts prepared from needles of unhardened
and hardened seedlings of Pi nus sylvest ris[J].Plant Sci.L ett,
7:359～369.
[22]　Rai LC,Tyagi B,Mallick N,et al.1995.Interative effects of
UV2B and copper on photosynthetic activity of the cyanobacteri2
um A nabaena doliol um[J].Envi ron.Ex per.Bot.,35:177～
186.
[23]　Rai LC,Tyagi B,Mallick N.1996.Alterlation in photosynthetic
characteristics of A nabaena doliol um following exposure to UV2
B and Pb[J].Photochem.Photobiol.,64:658～663.
[24]　Rai LC,Tyagi B,Rai PK,et al.1998.Inferactive effects of UV2
B and heavy metals(Cu and Pb)on nitrogen and phosphorus
metabolism of a N22fixing cyanobacterium A nadaena doliol um
[J].Envi ron.Ex per.Bot.,39:221～231.
[25]　Teramura AH,Sullivan J H.1994.Effect of UV2B radiation on
photosynth esis and growth of terrestrial plants[J].Photosynth.
Res.,39:463～473.
[26]　Towill L,Weiser E.1974.Studies on the reduction of2,3,5,2
triphenyltetrazolium chloride as viability assay for plant tissue
cultures[J].Can.J.Bot.,53:1097～1102.
[27]　Vu CV,Allen JL,G arrard LA.1982.Effects of supplemental
UV2B radiation on primary photosynthetic carboxylating enzymes
and soluble proteins in leaves of C3and C4crop plants[J].Physi2
ol.Plant,55:11～16.
作者简介　杜英君,男,1963年生,讲师,博士。

主要从事生理生态学方面研究,发表论文数篇。

责任编辑　魏丽萍
(上接第88页)
[15]　Malek parison of the Bowen ration2energy balance
and stability2corrected aerodynamic methods for measurement of
evapotranspiration[J].Theor.A ppl.Cli matol.,48:167～178.
[16]　Marshall J K.1971.Drag measurements in roughness arrays of
varying density and distribution[J].A gric.Meteorol.8:269～
292.
[17]　Monteith JL,Unsworth MH.1990.Principles of Environmental
Physics[M].London:Edward Arnold.
[18]　Saito T,Nagai Y,Isobe S,et al.1970.An investigation of turbu2
lence within a crop canopy[J].J.A gric.Meteorol.,25:205～
214.
[19]　Seginer I.1974.Aerodynamic roughness of vegetated surfaces
[J].Boum dary2L ayer Meteorol.,5:383～393.
[20]　Shaw RH,Pereira AR.1982.Aerodynamic roughness of a plant
canopy:a numerical experiment[J].A gric.Meteorol.,26:51～
65.
[21]　Takagi K,Miyata A,Harazono Y,et al.2003.An alternative ap2
proach to determining zero2plane displacement,and its application
to a lotus paddy field[J].A gric.For.Meteorol.,115:173～
181.
[22]　Watanabe T,K ondo J.1990.The influence of canopy structure
and density upon the mixing length within and above vegetation
[J].J.Meteorol.Soc.Jap.,68:227～235.
[23]　Wieringa J.1993.Representative roughness parameters for ho2
mogeneous terrain[J].Boundary2L ayer Meteorol.,63:325～
363.
[24]　Wooding RA,Bradley EF,Marshal J K.1973.Drag due to regular
arrays of roughness elements of varying geometry[J].Boudary2
L ayer Meteorol.,5:285～308.
作者简介　赵晓松,女,1980年出生,在读硕士。

研究方向为生态气候学。

发表论文数篇。

责任编辑　王　伟
511
杜英君等:UV2B辐射和Hg2+复合处理对黑小麦生理代谢和生长的影响。