基于Penman-Monteith方程模拟青海云杉生长季日蒸腾过程

左亚凡,贺康宁,*,柴世秀,俞国峰,李远航,林 莎,陈 琪,王琼琳

1 北京林业大学水土保持学院,水土保持国家林业局重点实验室；北京市水土保持工程技术研究中心；林业生态工程教育部工程研究中心,北京 100083 2 大通县气象局,青海 810100

森林作为地球上主要的陆地类生态系统,随着经济的发展,由于人类的不合理活动对环境造成一定破坏,导致温室效应的加剧,降水格局发生变化,尤其是黄土高原和青藏高原等干旱及半干旱地区,近五十年来的年均地表气温每10年就会上升0.2—0.3 ℃,增幅高于全国均值5—10倍[1-2]。林分蒸腾量的准确模拟,能够为林草植被优化与林分结构调整提供合理的理论依据,加强林分稳定性,提高水分利用效率,对于解决黄土高原地区因全球变暖所带来的干旱问题具有重要意义。

蒸腾作用是树木耗水的主要途径,它通过控制植物叶片气孔的开闭,将叶片内的水分散失到空气中的过程；植物的气孔成为了土壤-植物-大气连续体(SPAC)间进行水分传输的重要通道,气孔的开闭行为与气孔阻力密切相关,如何准确的确定气孔阻力,将成为精确模拟林分蒸腾量的前提。随着蒸散研究的不断发展,国内外学者陆续提出了不同的蒸散估算模型,其中 Penman-Moteith方程(以下简称PM方程)是目前公认的具有精度高、适用性强、可靠性高等特点的蒸散估算模型[3]。PM方程不仅能够以天为单位,还能以月为单位对蒸散进行模拟估算。然而PM方程的参数众多,其中气孔阻力是最为敏感且难以直接测得的参数[4],国内外对如何准确计算冠层气孔阻力已经做了大量的研究,Hatifield[5]、李召宝[6]通过实际观测得到的蒸腾量和各项环境因子,利用不同公式对小麦的冠层气孔阻力进行了反推估算；Szeicz等[7]则利用实测的叶片气孔阻力,结合叶面积的空间垂直分布结构,对大麦的冠层气孔阻力进行了估算；Tuner等[8]利用能量平衡原理结合冠层温度估算出水稻的冠层气孔阻力；叶片气孔阻力虽然可以通过Li- 1600等仪器直接测得,但是植物冠层内的叶片所接受的太阳辐射差异较大,并且叶片之间的生理特性也不同,无法精确得出总的冠层气孔阻力。

目前的研究还大多集中于田间农作物,对乔木的研究比较少,本文将利用反推法对冠层阻力进行估算,旨在建立能够准确估算蒸腾量的模型,以评价PM方程在林木冠层尺度的适用性。然而PM方程在实际运用中,常常把植被和土壤作为一个整体来计算,无法区分林木蒸腾和土壤蒸发,笔者利用冠层整体气孔阻力rsT代替原方程中的气孔阻力rc,对PM方程进行修正,从而把林木冠层蒸腾分离出来,更加准确的预测冠层蒸腾量。因此如何精确模拟rsT将成为本研究的重点。以青海省大通县塔尔沟小流域为例,采用反推法,利用修正后的PM方程反推求出rsT,并建立rsT与常规气象因子之间的回归方程,进一步模拟蒸腾量,辅以茎流计测定的蒸腾量实测值作为模型的检验指标,以期建立可以准确估算青海云杉生长季(5—10月)蒸腾量的模型。

青海云杉(Piceacrassifolia)是我国青藏高原东部祁连山地区的主要人工林树种,是该地区的顶级群落树种,多为纯林或与华北落叶松(Larixprincipis-rupprechtii)、白桦(Betulaplatyphylla)等乔木形成混交林,在维持区域生物多样性、促进水分循环、能量流动、碳增汇及生态系统维护中发挥着非常重要作用[9]。

1 研究区概况

研究区位于青海省西宁市大通回族土族自治县,该县地处青海省东部,祁连山东段的南麓,属于黄土高原与青藏高原的过渡区,位于东经100°51′—101°56′,北纬36°43′—37°23′,海拔大约2280—4622 m。大通深处内陆,远离海洋,地居中纬度地区, 海拔较高,属于大陆性高原气候。年辐射总量为547.03—615.93 kJ/cm2,全境年平均日照时数为2605 h。该地区气温垂直分布明显,昼夜温差较大,年平均气温为2.8℃,一年中气温最高的三个月为6、7、8月。大通县降水日数多、强度小,降雨量的季节分配极不均匀,多年平均降雨量为508 mm左右,年降雨量达到450—820 mm不等,其中植物生长季(5—10月)降雨量约占全年降雨量的87%。大通县境内土壤垂直分布较为明显,主要的土地类型有高山石质土、高山草甸土、山地棕褐土、黑钙土、栗钙土、潮砂土等。

大通县的森林植被主要分布在北川河及其支流的河谷两岸,主要的造林树种包括青海云杉(Piceacrassifolia)、青杨(Populuscathayana)、白桦(Betulaplatyphylla)、华北落叶松(Larixprincipis-rupprechtii)、祁连圆柏(Sabinaprzewalskii)、柠条(Caraganaintermedia)、沙棘(HippophaerhamnoidesLinn)、甘蒙柽柳(Tamarixaustromongolica)、宁夏枸杞(Lyciumbarbarum)、山杏(Armeniacavulgaris)等。

2 研究材料和方法

2.1 实验设计

实验地位于青海省大通县的青海云杉纯林固定样地(中坡,37°1′53.712″N,101°40′50.022″E),样地面积为25 m×25 m；样地平均树高为10.96 m；平均胸径为11.71 cm；平均冠幅度为7.05 m2,详见表1。在样地内选取三棵长势良好的青海云杉标准木作为实验对象,每木检尺数据如表2所示。

表1 样地基本信息

表2 样树信息

(1)茎流速率的测定

采用植物茎流计(Dynamax公司生产,Sapflow32,美国)进行全天候不间断观测。该茎流测量系统可以安装基于热扩散原理的TDP探针式热扩散传感器,安装方法为先在磨好的树干上钻孔,再将探针小心插入,最后用用橡皮泥和泡沫塑料密封空隙,防止水分渗入,最后用保温毯包裹。每间隔1 min测量一次数据,每30 min记录存储一次数据。所测的日蒸腾量为三棵标准木的平均值。

液流密度采用Granier于1987年提出的模型进行计算[10]:

(1)

(2)

式中,Js为液流密度(g cm-2min-1)；Tmax为探针之间的最大温差；T0为瞬时温差；EC为蒸腾速率(mm/min)；SA为边材面积(cm2)；A为冠幅投影面积(m2)。EC为30 min内的平均蒸腾速率,将时间尺度扩大到1 d,即为日蒸腾量。

(2)环境因子的测定

测量的环境因子主要为气象因素,测量的环境因子主要为气象参数,使用布设于样地林间空地的便携式气象站(DAVIS公司生产,Vantage Pro2 Plus,美国)对小气候环境进行观测。测量的主要参数有太阳辐射,大气温度,大气相对湿度,气压,降水量,风速等。

本文还考虑了温度(T,℃)和空气相对湿度(RH,kPa)的协同作用,采用水汽压亏缺(VPD,kPa)来表示:

(3)

式中,a、b、c分别为0.611 kPa、17.052、240.97℃[11]。

(3)消光系数、叶面积指数测定

叶面积指数(LAI)用冠层分析仪拍照定点观测,然后使用WinSCANOPY 2006a软件对照片进行分析处理计算。消光系数(k)的计算方法为:

(4)

式中,Rn为林冠上方接收到的太阳净辐射(J m-2d-1),Rns为穿透林冠到达林冠下的辐射(J m-2d-1)。

2.2 应用理论

(1)Penman-Monteith方程及其修正公式

Penman联合法是目前公认的适用性最强、计算结果精确可靠的估算方法。Penman于1948年首先提出无水汽水平运输情况下参考作物蒸发蒸腾量计算公式；Monteith在Penman公式的基础上,引入了气孔阻力(rc)的概念；随后,Bowen和Rosenberg提出在能量平衡模式中利用表面温度和近地面湍流运动的关系来改进Monteith提出的公式,并引入了中性层结下的空气动力学阻力(ra)。

现如今,Penman-Monteith的通用公式为:

(5)

式中,

(6)

(7)

(8)

在P-M公式中,假定大气为中性层结构,ra空气动力学阻力可采用下式计算[12]:

(9)

式中,λ为气化潜热(J/g),在气温为t时,λ=2498.9-2.33t；E为茎流速率(mm/d);α为时间尺度系数,时间尺度为小时时,α=3600,时间尺度为日时,α=86400；Δ为饱和水汽压关于温度的曲线的斜率(kPa/℃)；Rn为林冠接受的太阳净辐射(J m-2d-1)；ρ为大气密度(kg/m3)；Cp为定压比热(1012 J kg-1K-1)；ea为实际空气水汽压(kPa)；et为大气温度为t时的饱和水汽压；ra为空气动力学阻力(s/m)；rc为冠层阻力(s/m)t为大气温度(℃)；RH为空气相对湿度(%)；Z为蒸散层2 m以上高度(m)；Z0为粗糙度(m),此处取0.075H,H为树高；d为零平面位移(m)；u为Z高度处的实测风速(m/s)。

在假定动量、热量以及水汽运送的边界层阻力相差较小的情况下,即令rah(热传输阻力)≈rav(水汽传输阻力)≈ra也就是将林分的林冠层看作一个整体,假设是一片大叶子,而且在考虑气压订正后,可用冠层整体气孔阻力rsT替换冠层阻力rc[13],修正后的公式如下:

(10)

式中,P0/P为气压订正,P0/P=10L/(18400(1+t/273)),L为实验地的海拔高度(m)；k为消光系数；LAI为叶面积指数。

(2)潜在蒸散量计算

本研究采用基于PM方程的FOA-PM公式计算样地的潜在蒸散量,该公式在考虑了多种气象因素后,将物质传送方法与能量平衡相结合,是目前公认的误差最小的潜在蒸散量估算方法[14-15]。公式如下:

(11)

式中,ET0为潜在蒸散量(mm/d)；u2为2 m处的风速(m/s)；T为日平均温度(℃)；其余参数均与上文中介绍的参数意义相同

(3)冠层整体气孔阻力rsT的计算

冠层整体气孔阻力是修正后的PM公式中的关键参数,其大小由叶片气孔的开关程度所决定,而叶片气孔开闭程度主要受环境因子的影响[16]。上述公式中除冠层整体气孔阻力rsT外的各项参数,均可通过布设的仪器精确测量出来,目前,冠层整体气孔阻力rsT的主要推算方法为经验公式反推法,即通过测量出PM修正公式中的各项参数,包括实际的蒸腾速率,从而反推出rsT值；再构建出rsT和环境因子的回归模型,预测出生长季的冠层整体气孔阻力。

本文选取了日间(08:00—20:00)的三种常规气象因子(饱和水汽压差VPD,空气相对湿度RH,大气温度T)作为基本变量,探究其与冠层整体气孔阻力的响应关系。以往的研究发现,在实际的回归中,一个回归模型表达一整天的数据其模拟结果误差较大,因为rsT受到很多环境因子的影响,通常把一天的数据分时段回归,可达到较为理想的效果[17]。由于植被的蒸腾主要动力来源于太阳辐射,而夜间树木气孔几乎处于关闭状态,并受到很小的辐射通量,因此可假定青海云杉上方的净辐射量小于或等于土壤热通量时其蒸腾量为零[18],只对有太阳辐射的时间段(08:00—20:00)进行模拟。

为降低生长季内气候变化造成的误差进行分月模拟,利用以下日期内:5月10日—5月19日；6月1日—6月6日；7月1日—7月5日；8月1日—8月5日；9月1日—9月5日；10月1日—10月5日获取的数据进行建模。建立每个月份的回归关系式,最后带入PM方程对蒸腾量进行预测。

(4)误差分析

为进一步检验模型的准确性,将采用均方根误差(Root Mean Squared Error, RMSE)、平均绝对误差(Mean Absolute Deviation, MAE)和平均相对误差(Mean Relative Error, MRE)分析日蒸腾量的模拟值与实测值的误差,计算公式如下:

(12)

(13)

(14)

2.3 数据处理

使用EXCEL 2016对数据进行整理、误差分析、拟合系数计算,使用Origin 2018进行绘图,并用SPSS 20.0进行数据分析。

3 实验结果与分析

3.1 生长季日蒸腾量变化趋势

使用茎流计所测得的日蒸腾量如图1所示,在整个生长季中,青海云杉日蒸腾量的总体变化趋势不大,其峰值(4.55 mm)出现在7月10日,在8月27日达到最低值(0.05 mm),平均日蒸腾速率为1.20 mm,生长季总蒸腾量为220.21 mm,总潜在蒸散量为305.20 mm,日蒸腾量与日潜在蒸散量变化趋势基本相同。各月蒸腾量与潜在蒸散量的比值为7月(79.68%)>8月(72.71%)>6月(72.53%)>5月(67.08%)>9月(66.48%)>10月(64.29%),生长季日蒸腾量与潜在蒸散量变化曲线均为单峰型。

图1 青海云杉生长季内日蒸腾量、日潜在蒸散量、日均温及日降雨量变化Fig.1 Changes in daily evapotranspiration, daily potential evapotranspiration, average daily temperature and daily rainfall during the growing season of Qinghai Spruce

3.2 时滞效应

由于本研究采用反推法利用茎流速率实测值求得rsT,因此务必考虑冠层蒸腾与树干茎流是否存在时间滞后问题。本文将通过反推求出的rsT与饱和水汽压差(VPD)、空气相对湿度(RH)和温度(T)分别进行相关性分析,以0.5 h的间隔频率进行错位移动,探求滞前3 h到滞后3 h内的相关系数变化规律,具体如图2所示。经过对比发现,rsT与T和VPD的相关系数均在滞后0.5 h时达到最高值,分别为0.714和0.685；而rsT与RH的相关系数在无时滞0 h达到了峰值(R=0.734),与滞后0.5 h(R=0.722)结果十分接近,造成这种情况的原因可能是进行错位移动的间隔略大,其相关系数的最大值可能出现在0—0.5 h之间,为了简化研究过程,本文的后续研究均以滞后0.5 h作为校正。

图2 树干茎流对不同气象因子的时滞变化 Fig.2 Time-lapse changes of stem flow to different meteorological factors T:温度Tempeture；RH:空气相对湿度Relative humidity；VPD:饱和水汽压差Saturated water vapor pressure difference

3.3 冠层气孔阻力rsT与环境因子的响应

水分胁迫会影响树木冠层整体气孔阻力,进一步影响蒸腾速率。青海云杉生长季内,样地的土壤水分动态变化如图3所示,土壤深度50 cm以下的含水量变化基本控制在15%—25%之间,在整个生长季内变化不大。对土壤含水量与日蒸腾量做进一步的相关性分析,发现两者之间的相关性极低,因此样树在生长季内的蒸腾速率受到土壤水分胁迫的影响很小。本实验在研究冠层整体气孔阻力rsT时,可忽略土壤水分对其造成的影响,主要考虑常规气象因子。

图3 生长季样地土壤水分动态图Fig.3 Dynamic map of soil moisture in growing season plots

3.3.1rsT与气象因子的单因素回归

以冠层整体气孔阻力rsT为因变量,温度(T)、空气相对湿度(RH)和饱和水汽压差(VPD)为自变量,应用SPSS 20.0软件进行回归拟合。

冠层整体气孔阻力rsT总体上与温度(T)呈负相关关系,随着温度升高,叶片活性提高,气孔导度增大,气孔阻力减小,以进行蒸腾作用,如图4所示；而空气相对湿度(RH)与冠层气孔阻力rsT呈正相关关系,具体如图5所示,这是因为空气相对湿度的增加,rsT增大,会造成植物叶片内外的水汽压差减小,从而抑制蒸腾作用。如图6所示,rsT与饱和水汽压差(VPD)呈负相关,是因为温度的升高会导致叶片活性增加,气孔导度增大,从而使冠层整体气孔阻力减小；而林木冠层的饱和水汽压差是水汽扩散的主要动力,影响冠层蒸腾作用中水汽传输过程,所以当VPD增大时,叶片通过调节自身气孔导度来保持体内的水分平衡。

图4 青海云杉生长季各月份冠层整体气孔阻力与温度的关系Fig.4 Relationship between overall stomata resistance and temperature of spruce in Qinghai in each month of growing season

图5 青海云杉生长季各月份冠层整体气孔阻力与空气相对湿度的关系Fig.5 Relationship between the overall stomata resistance of the canopy and the relative humidity of the air of Qinghai spruce in each month of the growing season

图6 青海云杉生长季各月份冠层整体气孔阻力与饱和水汽压差的关系Fig.6 Relationship between the overall stomata resistance of the canopy and the saturated water vapor pressure of Qinghai spruce in each month of the growing season

3.3.2rsT与气象因子的多因素回归

冠层整体气孔阻力rsT与气象因子的多因素回归模型如表3所示,通过对比发现多因素回归模型的决定系数R2普遍高于单因素回归模型。

表3 冠层整体气孔阻力rsT与气象因子的多元回归方程

3.4 模型验证

为了进一步验证模型的模拟精度,以更加准确的评价PM方程在林木冠层尺度的适用性,将所建立的气象因子与rsT的多元回归模型带入PM方程中,采用5月20日—5月31日；6月7日—6月30日；7月6日—7月31日；8月6日—8月31日；9月6日—9月31日；10月6日—10月31日的连续蒸腾数据对模型预测结果进行验证,具体如图8所示,发现蒸腾量的变化趋势基本相同,与实测结果相近,但也有极个别天数的变化趋势与实测结果相反。

图7 青海云杉生长季内各月份的日蒸腾量实测值与模拟值的比较Fig.7 Comparison of Measured and Simulated Daily Transpiration of Qinghai Spruce in Each Month

误差分析结果如表4 所示,平均相对误差介于9.118%—19.254%之间,最低均方根误差为0.082,最低平均绝对误差为0.07 mm；累计均方根误差为0.200,累计平均绝对误差为0.160,累计平均相对误差为14.381%,误差范围比较合理。

表4 青海云杉生长季各月份日蒸腾量的模拟精度

4 讨论

4.1 生长季日间蒸腾特征

根据茎流计实测值,青海云杉在生长季内(5—10月)总日间(08:00—20:00)蒸腾量为220.21 mm,低于刘敏等[19]在青海省大通县所测得的青海云杉耗水量(306.45 mm),造成差异的原因可能是树龄不同,本研究的青海云杉为20—30年生,而刘敏等测量的云杉为50年生。万艳芳[20]研究发现祁连山地区的青海云杉平均液流密度随月份变化呈先增大后减小的趋势,即7月份最大,其次是6月和8月,9月和10月明显减小,这与本研究的结果一致。同时研究发现,8月份的日均温度略大于7月份,但是8月蒸腾量(35.58 mm)却显著低于7月份(61.01 mm),可能是由于8月份该地区阴雨天较多,日照时间小于7月份,导致空气湿度增加,造成植物叶片内外的水汽压差减小,从而抑制了植物的蒸腾作用。李世荣等在青海省大通县的研究发现,小麦、华北落叶松林和青海云杉天然林蒸腾量约占当地平均潜在蒸散量的80.5%—95.0%；幼龄人工混交林林相对较小,只占潜在蒸散的52.8%—66.8%[21]；与本文研究结果相近(64.29%—79.28%)。

4.2 时滞效应

本文采用错位对比分析法确定树干茎流滞后于饱和水汽压差、温度0.5 h,胡兴波等基于茎流计确定冠层蒸腾对气象因子变化的响应时滞为15 min[3]；陈琪等对华北落叶松的研究发现树干液流滞后于气象因子70 min[22]；赵平等研究表明马占思树具有较强的滞后效应,样树的液流变化滞后于光和有效辐射40—110 min[23]。

4.3 模型模拟

在青海云杉生长季内,饱和水汽压差、温度分别与冠层整体气孔阻力呈负相关关系,空气相对湿度与冠层整体气孔阻力则呈正相关关系。这与李召宝等[6]对冬小麦冠层气孔阻力的研究,孙林等[24]和陈琪等[22]对华北落叶松的研究结果相一致,但是环境因子对于不同树种的不同生长阶段有不同程度的响应,因此在相关性上会有一定差异。

在对林木蒸腾估算的研究中,大多是以小气候因子为自变量,与蒸腾建立起相关的经验预测模型。张岩等[25]通过建立rsT与气象因子的关系式,应用常规气象数据与PM方程对杨树生长季的蒸腾量进行连续模拟计算；戴剑锋等[26]通过建立气孔阻力与太阳辐射之间的定量关系,模拟了番茄的累计蒸腾量,但其仅建立出冠层气孔阻力和环境因子的单因素回归式来进行预测,而冠层气孔阻力往往会受到多种环境因子的影响,利用多因素回归方程来确定冠层整体气孔阻力,得到的模拟结果会更加可靠。李海光等[27]基于MATLAB对华北落叶松的蒸腾量进行仿真设计,构建了气象因子和蒸腾耗水的数学模型,但是该方法的局限性较大,普适性较低。葛亮等[28]以三种环境因子为自变量建立了估算柑桔树植株蒸腾量的经验回归模型,相较于PM方程更接近实际数据,但是其对气象数据的测量精度的要求较高,且目前只停留在盆栽试验。

本文通过建立冠层整体气孔阻力的多元回归方程,对青海云杉生长季的日间蒸腾量进行模拟,平均相对误差控制在了20%以内,最小为6月份9.862%,最大为7月份18.725%,母艳梅等[29]采用Noilhan提出的经验模型计算冠层阻力,通过PM方程模拟出太行山南麓栓皮栎-侧柏-刺槐人工混交林的蒸腾量比实际值相差21%；李仙岳等[30]利用Jarvis模型对樱桃冠层气孔阻力进行推导,蒸腾量的相对误差控制在20%以内,最小为12.12%,模拟误差与本研究较为相近。研究发现,7月份的rsT与蒸腾量的拟合精度均为生长季各月份中的最低值,赵华等[31]的研究结果显示,极端的高温天气以及高温日数可能会对模拟精度产生影响；另外,李仙岳等[30]对于樱桃的研究发现高辐射强度下的冠层气孔阻力的拟合精度小于低辐射强度；根据本研究所在地区林业站资料,7月份日照时数达到206.9 h,超过25℃高温的天数达到9天,为生长季中最高值,这可能是导致生长季中7月份模拟蒸腾量误差较高的原因之一。

本文并未对夜间时段进行研究,由于夜间不存在太阳辐射,植物叶片的气孔开闭行为变得较为复杂,不能单纯的建立回归方程来估算rsT,而且植物的夜间液流通量包括组织补水和夜间蒸腾,无论是从模拟还是实测的角度来看,都无法有效的将夜间蒸腾分离出来,因此还需进一步探讨。同时,随着科学技术的不断发展,涡度协方差技术越来越成熟,未来有望将该技术引入PM方程中,能够更准确的获取相关的气象因子从而提高模拟精度。

5 结论

(1) 青海云杉生长季(5—10月)日蒸腾量随月份而变化,呈先增高后降低的趋势,总体为单峰曲线,月蒸腾量占潜在蒸散量的比例分别为7月(79.68%)>8月(72.71%)>6月(72.53%)>5月(67.08%)>9月(66.48%)>10月(64.29%)。

(2)青海云杉树干茎流与气象因子表现出一定时滞效应,具体表现为:树干茎流与饱和水汽压差(VPD)和温度(T)均在滞后0.5 h相关系数达到最高值,分别为0.714和0.685,与空气相对湿度(RH)的相关系数在无时滞0 h时达到最高值(R=0.744),与滞后0.5 h时(R=0.722)非常接近。

(3)冠层整体气孔阻力(rsT)与常规气象因子建立的回归模型中,多因素回归模型的决定系数普遍高于单因素回归模型,rsT与温度(T)、饱和水汽压差(VPD)均为负相关关系,与空气相对湿度(RH)为正相关关系。

(4)蒸腾量的模拟验证表明,蒸腾量模拟值与实测值的变化趋势基本相同,累计平均相对误差为14.381%,累计平均绝对误差为0.16 mm,累计均方根误差为0.2；证明所建立的模型模拟效果基本可行。