基于改进的双作物系数法估算辽河三角洲芦苇湿地蒸散量

于文颖,纪瑞鹏,*,贾庆宇,冯 锐,武晋雯,张玉书

1 中国气象局沈阳大气环境研究所, 沈阳 110166 2 辽宁省农业气象灾害重点实验室, 沈阳 110166

蒸散过程是陆地水文循环的重要组成部分,估算蒸散量对理解土壤-植物-大气连续体(SPAC,Soil-Plant-Atmosphere Continuum)系统中的水分平衡和水分传输具有重要意义[1]。蒸散(ET,Evapotranspiration)的两个组分蒸发(E,Evaporation)和蒸腾(T,Transpiration)分别反映了水分从土壤和植物向大气传输时的非生物过程和生物过程[2],ET组分的拆分能更明晰不同水源对生态水文过程的贡献[3],目前对于不同生态系统ET每个组分的重要性仍未达成共识[4]。未来气候变化将显著影响和改变全球水循环过程,蒸散过程及其组分的拆分是理解碳水循环变化的关键,准确估算各组分可以深入了解植被对水文过程的影响程度,准确估算蒸散量对于我们理解陆地生态系统至关重要。

蒸散是湿地水分支出的主要途径,湿地蒸散量组分的拆分,对于了解湿地的水文过程和湿地蒸散过程对区域水循环的贡献等方面有着重要意义[5]。目前,湿地蒸散量的监测方法主要有蒸发皿法、同位素法、大孔径闪烁仪法、涡动相关法以及遥感法等,其中,涡动相关法被公认为精度较高的观测方法,获取的数据可用于模型估算法的校正。湿地蒸散量模型估算法主要有：Thornthwaite 公式、Penman模型、Penman-Monteith模型、FAO- 56模型、Hammer-Kadle 经验方程、Priestly-Taylor模型、波文比能量平衡法等,其中FAO Penman-Moteith公式被认为是目前计算参考作物蒸散量的标准方法[6]。湿地蒸散量通常利用参考作物蒸散量和作物系数来估算,该方法的关键是确定作物系数。作物系数的计算方法包括单作物系数法和双作物系数法,单作物系数把植被蒸腾和土壤蒸发作为一个整体来考虑,适合于均一完全覆盖植被,双作物系数将土壤蒸发系数和植被蒸腾系数区分开,反映了植被表面蒸散过程机理[7]。

国际粮农组织(FAO,Food and Agriculture Organization)规定了各时期典型作物系数的经验常数[8]。研究表明,作物系数受生物因子和环境因子的共同影响,主要包括植被的高度、覆盖度、叶面积指数、地表湿度、生长季节、地理位置、气候和环境因子等[9-11]。为了选择适用于不同类型湿地植被的作物系数以准确地评估湿地蒸散量,国内外学者对FAO推荐的作物系数进行了修正。贾志军等[12]利用线性方程修正了FAO推荐的3个发育期作物系数,并计算了三江平原毛苔草沼泽湿地蒸散量；周林飞等[13]根据环境因子的变化将FAO推荐作物系数调整为4个发育期的作物系数,并估算了石佛寺人工芦苇湿地的蒸散量。Zhou等[14]利用空气温度、相对湿度和净辐射作为输入参数建立了日尺度作物系数模型,并估算了芦苇湿地蒸散量。

芦苇是湿地生态系统的主要植被类型之一,在我国北方分布广泛[15]。湿地蒸散包括湿地水面-土壤-植被表面的蒸发过程、植被蒸腾过程和植被表面截留蒸发过程,对于芦苇湿地,表面截留蒸发量较小,湿地蒸散以水(土)表面蒸发过程和芦苇群落蒸腾过程为主。常用的双源蒸散模型包括Shuttleworth-Wallace模型(S-W)和FAO- 56双作物系数法,S-W模型计算参数多,相对复杂[16]。FAO- 56模型的双作物系数法[17]是一个简单易行的方法,但由于不同类型湿地植被双作物系数的缺乏而使得FAO- 56模型的应用受到限制[18-19]。因此,双作物系数的确定是模拟芦苇湿地蒸发与蒸腾过程需要解决的关键问题。

辽河三角洲湿地是亚洲最大的暖温带滨海湿地,其中芦苇面积达到756km2,是亚洲第一大芦苇湿地[20],模拟辽河三角洲芦苇群落的蒸发与蒸腾过程,对于明确芦苇湿地的水面-土壤-植被表面的水分交换、传输过程和平衡关系有着十分重要的意义。本文以辽河三角洲湿地芦苇群落为研究对象,利用涡动相关水汽通量、小气候梯度要素、群落内水面蒸发和芦苇群落生长参数等数据和FAO- 56模型计算双作物系数,分析作物系数动态变化过程及其主导影响因子,基于生物因子和环境因子建立时间尺度为小时的双作物系数模型,为实现芦苇湿地蒸发与蒸腾过程的分离提供参数,同时可为湿地蒸散量的估算和水资源评估提供依据。

1 数据与方法

1.1 研究区概况

辽河三角洲位于40°45′—41°10′N, 121°30′—122°00′ E,湿地生物类型多样,自然湿地植被主要包括芦苇(Phragmitesaustralis)、碱蓬(Suaedaglauca)、怪柳(Tamarixchinensis)、獐毛(Aeluropussinensis)、羊草(Leymuschinensis)、罗布麻(Apocynumvenetum)、香蒲(Typhaorientalis)等[21-22],其中,芦苇沼泽湿地的面积占辽河三角洲总面积的20%以上[21]。该区域属暖温带大陆性半湿润季风气候,四季分明,雨热同期,年平均气温8.3℃,年降水量611.6mm[23]。

研究区(图1)选在中国气象局沈阳大气环境研究所盘锦湿地生态系统野外观测站(40°56′N, 121°57′E)。该站位于辽宁省盘锦市大洼县赵圈河芦苇自然保护区内,保护区内芦苇湿地保存完好,其下垫面代表性较好。芦苇生长季和淹水时间为4—10月,芦苇平均高度2—3m,10月以后芦苇进入枯萎期[24]。

图1 辽河三角洲芦苇湿地观测站位置及实景图Fig.1 Locations of the research station and actual scenery in the Liaohe River Delta wetland, China

1.2 试验观测及数据处理

芦苇观测场配备有小气候梯度观测系统(CAMS620-GS, Huatron, China)、涡动相关观测系统(Li-cor, Inc, USA)、固液态蒸发与降水记录仪(PD210, Huatron, China)等仪器设备。其中涡动相关观测系统传感器安装高度3m,提供30min频率的瞬时通量数据,小气候梯度观测系统提供了5个高度(1, 2, 8, 16, 30m)30min频率气象数据,包括风速、气温、相对湿度、水汽压等[25]。本文数据包括2016年芦苇生长季涡动相关系统的潜热通量、感热通量、小气候梯度的气象数据、降水量、水面蒸发量,土壤热通量等数据。在芦苇生长季不同生长阶段,选择天气晴朗日观测冠层叶面积指数和植株高度,并记录发育期。采用Li- 2000冠层分析仪观测叶面积指数,每次选择芦苇样方5个,重复10次,取其平均值；每次选择10株芦苇,利用量尺测量其高度,取其平均值。

涡动相关实测数据常用于模型计算或遥感反演蒸散值的精度检验,且广泛应用于各种类型湿地。利用Eddypro 5.0.1软件对涡动相关系统获取的通量数据进行坐标旋转、密度效应校正等处理,转化为小时尺度观测频率的潜热通量数据和蒸散量,并进行QA/AC奇异值剔除。

1.3 研究方法

1.3.1参考作物蒸散量的计算

利用FAO- 56 Penman-Monteith模型计算参考作物蒸散量,其表达式为：

(1)

式中,ET0为参考作物蒸散量,mm/d；Rn为作物表面净辐射,MJ m-2d-1；G为土壤热通量,MJ m-2d-1；T为平均气温,℃；es为饱和水汽压,KPa；ea为实际水汽压,KPa；Δ为饱和水汽压-温度曲线斜率,KPa/℃；γ为湿度计常数,KPa/℃；u2为距地面2m高处风速,m/s。

1.3.2双作物系数法

FAO- 56模型的双作物系数法将蒸散量分为土壤蒸发和植物蒸腾两部分计算,作物系数是实际蒸散量(ETc)和参考作物蒸散量(ET0)的比值,它综合反映了环境因素和植物对蒸散的影响,包括空气动力学阻力、表面阻力、植物品种和植被长势等。

FAO- 56 双作物系数法数学表达式为：

ETc=KcET0

(2)

式中,ETc为实际蒸散量,Kc为作物系数,可分解为如下式：

Kc=KsKcb+Ke

(3)

式中,Kcb为反映植物蒸腾的基础作物系数；Ke为土壤表面蒸发的蒸发系数；Ks为水分胁迫系数,辽河三角洲湿地水分供应充足,故本研究中Ks=1。

Kc=Kcb+Ke

(4)

在芦苇生长季(芦苇生长季和淹水时间均为4—10月),芦苇根部完全被水淹没,因此,将公式中的土壤表面蒸发系数(Ke)代换为水面蒸发系数Kw。

Kc=Kcb+Kw

(5)

若已知作物系数的各分量,就可以利用FAO- 56模型(式2)计算芦苇湿地的蒸腾量、水(土)面蒸发量以及实际蒸散量。

1.3.3双作物系数模型改进

将辽河三角洲湿地芦苇发育期划分为4个阶段：生长初期、快速生长期、稳定生长期和生长末期。对芦苇不同发育阶段的Kc、Kcb和Kw与生物因子和环境因子进行拟合。首先分析Kc与冠层叶面积指数LAI、芦苇高度h等生物因子及大气温度Ta、大气相对湿度RH等环境因子的相关关系,找出其主导影响因子,利用主导影响因子重新构建双作物系数模型,对其参数进行拟合。拟构建方程如下：

Kc=a×f(LAI,h) ×f(Ta,RH,…) +b

(6)

安装在芦苇群落中的水面蒸发仪连续监测得到水面蒸发量E,利用下式计算实际水面蒸发系数Kw。

Kw=E/ET0

(7)

将涡动相关数据观测的蒸散量ETc作为实际值计算Kc。

Kc=ETc/ET0

(8)

已知Kc和Kw即可用式(5)计算出Kcb。

2 结果与分析

2.1 芦苇湿地ET0、ETc和E日变化

芦苇群落ETc日变化表现为单峰或多峰曲线。晴天条件下,凌晨0时—6时蒸散低,随着气温升高,蒸散量不断升高,到12时和14时左右出现峰值；而后,随着气温下降,蒸散不断下降,夜间蒸散量下降到最低(图2)。多云或阴雨天气下,蒸散将呈现多峰曲线变化。基于FAO- 56 Penman-Monteith模型,利用小气候梯度、通量等数据计算的ET0与利用涡动相关系统观测的ETc变化趋势基本一致,但芦苇湿地ETc在数值上高于ET0。

晴天条件下,E的日变化与ETc的日变化趋势基本一致,呈现早晚低、中午高的波动曲线变化(图2)；多云或阴雨天气条件下,E的日变化与ETc的日变化略有不同,但同样呈现多峰波动变化。从数值上看,E在芦苇生长初期高于快速生长期、稳定生长期和生长末期。

图2 实际蒸散量ETc、水面蒸发量E与参考作物蒸散量ET0日变化Fig.2 Daily dynamics of actual evapotranspiration (ETc), water evaporation (E) and reference evapotranspiration (ET0) ETc: 实际蒸散量 Actual evapotranspiration; E: 水面蒸发量Water evaporation; ET0: 参考作物蒸散量Reference evapotranspiration

2.2 作物系数日动态及其影响主导因子

2.2.1Kc、Kcb和Kw日动态

在芦苇生长初期,Kc和Kcb的日动态呈现多峰波动变化,全天波动幅度相对稳定,波动范围在0.6—1.2之间,其中白天(7时—17时,下同)波动范围较小,在0.9—1.1之间；在快速生长期和稳定生长期,Kc和Kcb夜晚波动幅度较大,而白天波动幅度较小,波动范围在0.9—1.3之间；生长末期,Kc和Kcb夜晚波动幅度较大,白天呈现多峰波动曲线,波动范围在0.8—1.2之间；Kw基本呈现白天低、夜晚高的波动曲线,在快速生长期、稳定生长期和生长末期白天的数值显著低于生长初期(图3)。

图3 作物系数Kc、基础作物系数Kcb和水面蒸发系数Kw日动态Fig.3 Daily dynamics of crop coefficient (Kc), basal crop coefficient (Kcb), and water evaporation coefficient (Kw ) Kc: 作物系数 Crop coefficient; Kcb: 基础作物系数 Basal crop coefficient; Kw: 水面蒸发系数 Water evaporation coefficient

2.2.2Kc、Kcb和Kw与生物和气象因子的相关性

基于2016年的涡动水汽通量、小气候梯度要素、水面蒸发量等数据,计算芦苇群落小时尺度和日尺度Kc、Kcb和Kw,并将其与叶面积指数LAI、株高h、气温Ta、相对湿度RH和风速u等因子进行相关性分析,结果表明,小时尺度Kc与Ta、LAI、h显著相关,Kw与RH显著相关,Kcb与Ta、RH、u、LAI、h均显著相关(表1)；日尺度3个系数与LAI、h显著相关,与环境因子相关性不显著(表2)。

为了避免估算小时尺度蒸散量时忽略掉生物因子的作用,在小时尺度数据相关分析中同样加入了生物因子。基于整个生长季的各要素相关性分析结果,发现日尺度3个系数主要受生物因子的影响,小时尺度3个系数在短期内主要受环境因子影响,而在较长时间范围(如整个生长季)受生物因子和环境因子的共同影响。

表1 小时尺度作物系数与生物因子、气象因子的相关关系

*P<0.05 水平显著相关, **P<0.01 水平显著相关

2.3 双作物系数模型构建及检验

基于观测站小时尺度的连续观测数据,参照FAO推荐的单、双作物系数方程,结合非线性回归法,建立基础作物系数Kcb和水面蒸发系数Kw的模型,方程如下：

Kcb =[a1(Ta-20)+b1(u-2)+c1(RH-45)]LAId+e (9)

*P<0.05 水平显著相关,**P<0.01 水平显著相关

利用非线性方程回归迭代法估算方程中的参数,参数拟合结果为：

a1=-0.702,b1=-0.651,c1=-0.066,d=-2.872,e=0.979,a2=-0.002,b2=-0.277,c2=0.138

2.4 改进的双作物系数法模拟芦苇湿地蒸腾与蒸发过程

图4 实测与模拟蒸散量比较 Fig.4 Comparison of the calculated evapotranspiration and measured values

利用模拟的作物系数Kc计算小时尺度的芦苇实际蒸散量,并与涡动相关系统观测的实际蒸散量进行比较,决定系数R2达0.894,模拟效果很好(图4)。

利用改进的作物系数法不仅可以计算芦苇实际蒸散量,同时也能拆分芦苇蒸腾过程和蒸发过程。分别以生长初期、快速生长期、稳定生长期和生长末期为例,模拟芦苇群落蒸腾和水面蒸发量的日动态变化过程(图5),其中,水面蒸发量分别占芦苇总蒸散量的9.0%—16.3%、7.7%—13.1%、 7.4%—12.3%和11.6%—15.9%。随着芦苇的快速生长,群落对水面的覆盖度越来越高,芦苇蒸散以蒸腾为主,棵间水面蒸发量占总蒸散量的比重越来越小。当芦苇叶面积指数、覆盖度达到最大值后,芦苇开始呈现衰落趋势,总蒸散量也将减少,水面蒸发量所占比重增加。

图5 芦苇湿地蒸腾和蒸发过程模拟Fig.5 The transpiration and evaporation process in reed wetland

辽河三角洲湿地芦苇4月20日左右萌芽,在生长初期蒸散量较低,随着芦苇的快速增长,一般叶面积指数在7月份达到最大值,湿地蒸散量也通常在7月份达到最大值,9月之后进入凋萎期,蒸散量随之减少[26](图6)。研究认为,芦苇覆盖度的增加能够增加湿地下垫面的蒸散量,同时,芦苇群落蒸腾量的增加量将远大于由于芦苇遮荫引起的棵间水(土)面蒸发的减少量[27],这与本研究的结果一致。

图6 芦苇湿地蒸散量日变化Fig.6 Diurnal variation of evapotranspiration in reed wetland

3 讨论

3.1 基于改进的双作物系数与FAO推荐作物系数模拟蒸散过程的结果比较

当芦苇生长茂密,芦苇对下垫面覆盖度较高,此时水面蒸发很小,芦苇蒸散主要以冠层蒸腾为主,但在芦苇生长初期,水面裸露较多,水面蒸发无法忽略时,芦苇群落的蒸散过程包括了芦苇冠层蒸腾以及水面蒸发过程。辽河三角洲芦苇湿地的蒸散过程以芦苇蒸腾为主,5月份的水面蒸发量大于6—9月份,而6—9月份由于芦苇覆盖度较高,水面蒸发几乎可以忽略不计。在芦苇生长初期,作物系数随着环境因子的变化,日变化也呈现波动过程,因此,单一的作物系数经验常数不适于小时尺度的蒸散模拟计算。在芦苇生长快速和稳定生长期,作物系数受环境因子影响减小,白天的变化呈现较为稳定的过程,在计算白天的蒸散量时,可以采用简化的经验常数。

FAO推荐的芦苇沼泽湿地(有水层)的作物系数Kc在生长初期、中期和末期分别为1.0,1.2和1.0。利用FAO推荐的作物系数计算辽河三角洲芦苇湿地蒸散量与涡动相关系统观测的实际蒸散量进行比较,相关系数R2为0.837(图7)。以芦苇生长中后期为例,分别利用FAO推荐的作物系数和利用改进的双作物系数法计算芦苇湿地蒸散量(图8),基于改进双作物系数法的蒸散量、基于FAO推荐作物系数的蒸散量和涡动相关实测值相比,发现三者的趋势基本一致,但是利用FAO推荐的作物系数计算蒸散量模拟值偏高,尤其是在芦苇生长中期和后期；而利用改进的双作物系数法计算的结果与实测值比较接近。利用FAO推荐的作物系数模拟值偏高的原因可能是辽河三角洲湿地的芦苇群落与FAO推荐的芦苇沼泽生长环境和生态类型不同所引起的差异,另外单一的作物系数无法完全代表芦苇长势与环境的变化。

由于地区和品种等差异,直接采用FAO推荐的各类作物系数估算蒸散量可能产生较大误差[28],研究普遍认为,植被群落的作物系数需要根据当地天气、植被特点、土壤特性以及其他因素进行修正[29]。因此,一些学者考虑了地理因素和气候特点等因素对FAO推荐的双作物系数法进行了修正。Wright[30]根据当地气候环境特点修正FAO推荐的作物系数；Rosa等[10,31]基于双作物系数法开发了SIMdualKc模型；Paredes等[32]利用修正后的作物系数校准SIMdualKc模型用于豌豆蒸散过程模拟；彭世彰等[33]修订了冬小麦、夏玉米、棉花及水稻作物的作物系数和土壤蒸发系数公式；张强等[34]针对春小麦利用风速和相对湿度修正了Kumar的作物系数方程；李谦等[35]等参照当地气候条件和作物叶面积指数订正稻田作物系数。冯禹等[36]基于叶面积指数修正了双作物系数中的土壤蒸发系数,用于估算春玉米蒸散。

图7 实测与FAO法模拟蒸散量比较 Fig.7 Comparison of the calculated evapotranspiration used FAO method and measured values

图8 实测与模拟蒸散量比较Fig.8 Comparison of the calculated evapotranspiration and measured values

基于FAO公式改进双作物系数法和利用双作物系数法估算农田、森林等生态系统蒸腾和蒸发量的研究较多,但以湿地植被为研究对象的较少。芦苇群落在生长季的覆盖度较高,水(土)面蒸发系数虽然在一定范围内可以忽略,但是通过对作物系数法的改进,发现利用改进后的双作物系数法估算湿地蒸散量比FAO推荐系数法更为准确,因此,所构建的双作物系数模型更适用于辽河三角洲湿地的芦苇群落,同时还可以拆分芦苇湿地的蒸发与蒸腾量。

3.2 调整FAO推荐的单作物系数经验值简化芦苇湿地蒸散模拟过程

改进的双作物系数法是基于环境因子和生物因子的小时尺度动态作物系数,能够更有效地模拟芦苇群落蒸发与蒸腾过程。但在实际应用中,如果仅计算群落蒸散量,可以利用改进的单作物系数简化计算。将芦苇生长季划分为生长初期、快速生长期、稳定生长期和生长末期4个阶段。利用公式(7)计算Kc实际值,将各个时期每日的动态系数求取平均值,4个阶段的作物系数平均值分别为0.9,1.1,1.0,0.9(表3)。

表3 辽河三角洲芦苇群落各时期调整作物系数

表4给出了FAO推荐的不同时期芦苇湿地作物系数经验值,以及其他文献修正后的芦苇湿地作物系数计算值。从表中可以发现,4个阶段芦苇湿地作物系数变化趋势基本一致,随着芦苇的生长,Kc值不断增加,并在成熟后下降。芦苇群落植被的高度、密度、覆盖度以及环境因子均影响着作物系数的大小。

作物系数受气候条件和下垫面特征等因素的影响[39],不同品种湿地植被的作物系数不尽相同。Drexler等[19]计算华金三角洲沼泽湿地芦苇的作物系数日均值为0.7—1.2,生长季平均值为0.95；Zhou等[14]计算的芦苇湿地作物系数范围为0.2—1.4,生长季平均值为0.71；Borin等[40]计算威尼斯和西西里岛芦苇湿地的作物系数范围分别为4.7—8.4和2.3—8.5；Anda等[40]计算巴拉顿湖芦苇湿地4—9月份作物系数分别为1.03,1.23,1.4,1.51,0.99,0.77,作物系数的季节平均值为0.73—1.37。Tuttolomondo等[41]计算西西里岛人工湿地植物伞竹在生长初期、快速生长期、稳定生长期和生长末期的作物系数分别为1.05,3.39,5.71,2.55,香蒲的作物系数分别为1.2,3.84,6.51,2.95。Queluz等[42]计算人工湿地香蒲的作物系数月平均范围为2.03—3.68。王昊等[38]利用FAO- 56模型和单作物系数法计算了扎龙芦苇湿地的蒸散量,在计算中调整了扎龙湿地芦苇各生长时期的作物系数,其生长初期、快速生长期、稳定生长期和生长末期的作物系数分别为1.2,1.191,1.178和1.083,与本文相比,各时期作物系数均略高于辽河三角洲芦苇湿地。贾志军等[6]研究认为三江平原沼泽湿地(主要植被毛苔草)的作物系数主要受叶面积指数LAI的影响,利用LAI建立了作物系数线性回归方程。本文的芦苇群落除了受自身株高、叶面积的影响外,由于研究区地处滨海地区,风速较大,风速是影响作物系数日变化的主导影响因子之一,因此作物系数方程的输入因子包括风速。周林飞等[13]根据FAO推荐的作物系数调整了石佛寺人工芦苇湿地的作物系数,调整后的生长初期作物系数为1.2,快速生长期为1.204,稳定生长期和生长末期加入风速、湿度和株高的影响因子,作物系数为变化值。邓雯等[43]认为采用FAO- 56推荐的单作物系数法计算滇池外海环湖湿地植被年蒸散量是合理的。李志威等[44]对若尔盖高原泥炭湿地植被初期、中期和末期修正后的作物系数分别为0.68,1.1,0.62,作物系数小于辽河三角洲芦苇湿地。

表4 FAO和其他研究的芦苇湿地作物系数(Kc)比较

FAO: 国际粮农组织Food and Agriculture Organization

4 结论

基于涡动通量、小气候梯度要素、水面蒸发等数据和作物系数的影响主导因子(气温、相对湿度、风速、叶面积指数),对辽河三角洲湿地芦苇群落的作物系数Kc、基础作物系数Kcb和水面蒸发系数Kw进行模拟,利用构建的双作物系数模型和FAO- 56模型,对辽河三角洲芦苇群落冠层尺度的蒸发与蒸腾过程进行模拟。基于生物因子和环境因子建立的动态双作物系数模型,适用于小时尺度的芦苇群落蒸散模拟,同时,实现了芦苇群落蒸发过程与蒸腾过程的分离,解决了无法通过实际观测直接获取芦苇群落蒸腾量的问题,芦苇群落蒸散量的模拟精度提高了5.7%,可为开展芦苇群落蒸发与蒸腾过程模拟提供参数。同时,调整了FAO推荐的芦苇单作物系数经验常数值,新调整的单作物系数更适用于辽河三角洲湿地芦苇群落的蒸散模拟。