基于系统动力学的桉树生长用水模型构建

刘 挺, 苟 思, 刘 超, 贺宇欣, 庄文化, 刘铁刚

(四川大学 水力学与山区河流开发保护国家重点实验室 水利水电学院, 四川 成都 610065)

桉树(Eucalyptus)作为世界上主要的速生造林树种之一，有着生长速度快、环境适应力强、产量高、易于种植等优点[1]。但其种植和管理的合理性也引起不少争议[2-4]，桉树种植是否会对当地水资源产生不利影响以及和人类用水是否相冲突等问题都有待解决。解决相关争论需要对桉树的生长和用水过程进行研究。目前的研究方法主要有两种，一种是通过试验观测的方法，比如通过设置对照流域进行森林集雨区试验、通过人工种植或砍伐一定面积的桉树等手段，观测人工林面积变化对径流的影响，但这需要较长的时间，实施和管理起来有诸多不便[5-7]。另一种是基于试验数据的模型模拟研究。本文开发模型研究桉树生长和用水过程。模型开发需要准确描述桉树的生长特性以及用水特性。桉树用水特性复杂，其根系具有二态性，即在距离地表较近的浅层土壤中分布有大量细微的侧根，而较粗的主根则径直向下，可穿透较深的土壤层甚至直达地下水层[8]。在不同季节桉树对土壤水和地下水的使用方式可能会有所不同。在雨季，桉树主要使用浅层土壤水，而在旱季当浅层土壤变干后，则会逐步使用深层土壤水或地下水。桉树根系存在水力再分配现象。水力再分配是指水分在水势差的驱动下，通过根系在不同土壤层间的被动运移。Adams等[9]和Schymanski等[10]在澳大利亚对桉树进行研究后指出水力再分配现象会显著影响桉树的用水过程。模拟桉树根系对不同水源的吸收和水力再分配现象，需要合理的根系吸水模型。在已有模型研究中，传统的模拟方法是先计算潜在腾发量，然后根据植物根系在土壤中的分布比例，从不同土壤层中同时吸水。在这样的模拟方法中，植物会同时使用浅层土壤水、深层土壤水和地下水等不同水源，与实际观测结论不符。Li等[11]开发的根系吸水的补偿模型机制简单，容易过高估计植被对于高含水量土壤层的吸水量。在模拟桉树根系的水力再分配方面，已有研究中Ryel等[14]提出的模型比较简单，模拟机制有待改进。Lee等[15]开发的模型有所改进，但无法模拟植物休眠或植物死亡情况下的水力再分配。而Gou，Miuer[10]提出的根系吸水模型有所创新，模型通过势能差驱动水分运移，能够模拟植物在不同季节对于不同水源的利用情况的变化，同时模拟了水力再分配作用。桉树生长受到诸多因素的影响。气象条件、人工管理、立地条件都会不同程度的影响桉树的生长。合理的考虑这些影响因素是模拟桉树生长的关键。已有研究中，ProMod模型[16]和3-PG模型[17]是模拟植物生长中两个最常用的模型。其中3-PG模型被成功应用于国内外桉树生长模拟中。在国外，3-PG模型较好的模拟出了澳大利亚西部以及塔斯马尼亚地区7个站点的桉树生长[18]；在国内花利忠等利用3-PG模型成功模拟了雷州半岛尾叶桉的生长规律[19]。桉树生长迅速，生长过程中用水量变化明显，桉树用水也会反过来影响桉树生长速率，因此需要将这两个框架有机结合在一起。本文构建的模型由两个主要框架组成，即桉树生长过程模拟和桉树用水过程及相关生态水文过程模拟。本文根据3-PG模型原理模拟桉树生长过程[17]，根据Gou和Miller[10]提出的根系吸水模拟、水分胁迫模拟、水力再分配模拟等原理模拟桉树用水和相关生态水文过程，通过合理地设置耦合点将两大框架耦合起来。利用系统动力学软件Vensim搭建桉树生长用水模型。最后在国内外分别选择两个研究区，验证模型的合理性，以期为科学地种植桉树提供理论帮助。

1 材料与方法

1.1 模型原理

1.1.1 桉树生长过程模拟 描述桉树生长的公式，主要利用Landsberg等[17]开发的3-PG植物生长模型原理来建立。该模型模拟林分水平的植物生长过程主要由生产量的计算、生产量的分配、新陈代谢和死亡、林分特征几个部分组成。生产量计算是通过植被光合作用将太阳光照拦截吸收以后合成有机物质。模型将呼吸作用消耗的能量占总能量的比值视为一个常数 ，默认为0.47。净初级生产量的计算如公式(1)所示[20]：

ΔPn=0.552αcγ(1-e-kL)Q

(1)

式中：ΔPn——单位时间内所产生的净初级生产量；k——吸收系数；L——叶面积指数；Q——辐射数据；αc——冠层量子效率，用于度量植被将拦截的太阳辐射转化为生产量的效率。αc受到其他环境因子比如温度、土壤营养状况、霜冻等的影响。接下来植被将净初级生产量分配到植物的根、茎、叶，分配情况取决于相应的分配参数，分配参数又和土壤肥力、胸高直径等参数有关。除了呼吸作用导致的损耗以外，生产量也按照一定的规律有所损耗。这部分损耗分为植物正常新陈代谢和植株死亡所造成的损耗，植株死亡又可以分为植株按照正常速率的死亡和由于植株密度达到一定程度以后出现资源的争夺所导致的植株自然稀疏现象。在林分特征方面，模型可以模拟树高、胸高直径以及林分材积的变化情况。这几个参数对于评估桉树人工林生产力，合理规划具有重要意义。相应的计算公式[20]为：

H=aHBnHBNnHN

(2)

B=(1 000WS/asN)1/ns

(3)

VS=(1-ρBB)WS/ρ

(4)

式中：H——树高(m)；N——植株密度(株/hm2)；B——胸高直径(cm)；ns——净初级生产量分配到茎干的比例；WS——茎干生物量；ρBB,ρ和林龄有关，参数aH,nHB,nHN,aS无实际意义，均为常数。

1.1.2 生态水文过程模拟 描述生态水文过程的公式，主要利用Gou和Miller[10]开发的生态水文模型原理来建立。桉树是有可能使用地下水的，该模型可以模拟植被对土壤水和地下水的联合使用。该模型通过势能差驱动水分在土壤和植被中的运移，并引入了新的根系吸水和根系水力再分配方程对传统的根系吸水模拟做出了改进；模型还引入了新的植物水分胁迫模拟方法[10]。

(1) 根系吸水模拟。在水势差的驱动下，根系吸水部分从土壤层中吸收水分，然后水分经过根系输水部分再向上或向下运输。如公式(5)所示[10]：

(5)

式中：Qs,i——根系从第i层土壤中吸收(或释放)水分的速率；n——总的根系数量；froot,i——第i层土壤中根系数量占总根系数量的比例；Aa_root,i——第i层土壤根系吸水部分区域的面积；Γs,i——第i层土壤中单位长度根系吸水部分根围的阻力；Γr,i——第i层土壤中单位长度根系输水部分阻力；ψsoil,i——第i层土壤的水势；ψref——参考点水势(参考点被定义为植被和地表齐平的点)；zi——第i层土壤中心离地表的距离。将各土层包括地下水层的根系吸水速率相加便可以计算出根系吸水总速率。在湿润季节，表层土壤含水量相对较高，水势较高，ψsoil,i-ψref-zi>0，Qs,i>0,根系从表层土壤中吸收水分，而深层土壤含水量相对较低，ψsoil,i-ψref-zi<0，Qs,i<0，根系反向释放水分给深层土壤。在干旱季节，深层土壤含水量相对较高，水势较高，ψsoil,i-ψref-zi>0，Qs,i>0,根系从深层土壤中吸收水分，而浅层土壤含水量相对较低，ψsoil,i-ψref-zi<0，Qs,i<0，根系反向释放水分给浅层土壤。因此该计算式可以模拟根系吸水和根系的水力再分配。模型认为水分的运移是连续的，因此根系从各层土壤及地下水层中吸水的总速率应该等于植被的腾发速率，从而得出公式(6)[10]：

(6)

(2) 水分胁迫模拟。在处理水分胁迫的问题时，传统的方法是基于土壤含水量的情况来考虑的。认为当土壤含水量降低到某个临界值时，植物开始受到水分胁迫，随着土壤含水量下降，植被的腾发速率也随之下降；当含水量下降到凋萎点，植被完全停止腾发。但对于能够利用较深层土壤水或地下水的深根系植被，在旱季，当表层土壤中的根系由于水分胁迫而减少甚至停止水分的吸收时，在其他水分充足的含水层中的根系可能仍然可以继续吸收水分。因此，模型建立了新的水分胁迫模拟方法：从植物空穴化角度出发，考虑水势、水力传导度、气孔开度等多个因素的相互影响，最终通过植被腾发速率的变化得以体现，即公式(7)[10]：

T=ETpotξ

(7)

式中：ξ——整个根系系统水力传导度在空穴化影响下剩余的比例； ETpot——潜在蒸散速率。

1.1.3 模型框架耦合 为将桉树生长过程的模拟和生态水文过程的模拟联系在一起，选取了3个关键点进行耦合。①耦合点1：考虑植物生长过程中叶生物量的变化对生态水文过程的影响。当植被未受到胁迫时(即植被的水分和养分供应充足、光线充沛、温度适宜)，根据生长模型得到的叶生物量为最大值。将该最大值作为基准值，计算各个时刻实际叶生物量和基准值的比值，用于修正生态水文模型中冠层覆盖度的变化，进而影响植物腾发速率。②耦合点2：考虑植物生长过程中根生物量的变化对生态水文过程的影响。通过生长模型中根生物量和叶生物量的比例变化来修正生态水文模型中根面积和叶面积的比值。③耦合点3：植物生长受土壤水分的影响。将生态水文模型当前的土壤含水量计算结果反馈给植物生长模型。植物生长模型据此考虑植物生长受到土壤水分的影响，进而计算生长模型中的土壤水参数。

1.2 系统动力学模型搭建

模型结构将在系统动力学软件Vensim中完成搭建。模型结构由各种变量组成。通过相应的联系符号和数学公式将各变量联系起来。Vensim中包括变量(variable)、存量(stock)、流量(flow)等。变量的类型取决于想要如何定义和使用某个要素。存量和流量是模型的核心变量。存量是系统中产生行为的基础，而流量则是让存量产生改变的原因所在。

在模拟桉树生长时，关键的存量包括根生物量、茎干生物量、叶生物量；关键的流量包括根、茎干、叶分配的净初级生产量、根周转率、叶凋零速率、植株死亡速率。各变量比如温度、树龄、水汽压差、辐射等对净初级生产量产生影响。净初级生产量通过输入流量被分配到根、茎干、叶相应的存量中。而在根周转率、叶凋零速率、植株死亡速率等输出流量的作用下，根、茎干、叶相应的存量又有所减少，处于动态变化的过程。在模拟桉树用水和生态水文过程时，各土壤层含水量被设置为存量，各土壤层之间的水分下渗、不同深度根系之间水分的运移、降雨入渗等均作为流量，使水分在各存量之间运移。同时，各土壤层、根系、叶片等都被赋予相应的能量变量，并和水分的运移系统紧密联系。水分运移，则能量同步发生变化，二者此消彼长，能量变量驱动水分运移，水分的运移又反馈造成能量的变化，最终带动整个模型运转。

1.3 研究区模型构建

为验证构建模型的合理性，选取位于澳大利亚的Tumbarumba站点和广东省雷州半岛的纪家林场作为相应的研究区域。

1.3.1 澳大利亚Tumbarumba研究区 Tumbarumba研究区位于澳大利亚新南威尔士(35°39′20.26″S，148°09′07.5″E)，海拔1 200 m[21]。年平均降雨量在986 mm左右，属于海洋型气候。冬季温度在4 ℃左右[22]，年平均温度在8 ℃左右[23]。站点的气象资料来源于Tumbarumba通量站。该通量站建立于2000年3月，为澳大利亚通量站系统之一，由隶属澳大利亚政府的联邦科学工业组织(CSIRO)和陆地生态系统研究网络(TERN)提供管理和支持。气温、风速、湿度、降雨、辐射等都是通过涡流技术测量(eddy flux)，土壤体积含水率是用时域反射仪(TDR)测量。目前可以获取的气象资料是从2001—2014年一共14 a的每1 h的气象观测数据。因此，选取2001—2014年作为Tumbarumba站点的模拟时长，选取每3 h作为一个模拟步长。降雨、气温、辐射等数据为模型输入数据。根据从澳大利亚土壤资源信息系统(ASRIS)获取的数据，耦合模型将土壤层分为了4层外加1层地下水层。残余土壤含水量、饱和土壤含水量、土壤饱和水力传导度等数据则是在参考Haverd等[24]采样获取的数据基础上设定(表1)。

表1 Tumbarumba站点主要土壤参数

该地区的桉树种植有着100 a上的历史，主要的桉树树种为德利格特桉(Eucalyptusdelegatensis)，林下植被较稀疏，以草的覆盖为主[21]。植被根系分布随土壤深度的变化的计算如公式(8)所示：

γ=1-βd

(8)

式中：γ——根系分布的比例；β——Jackson等[25]综合很多组根系研究的数据以后总结出的经验参数；d——土壤深度。

1.3.2 广东省雷州半岛的纪家林场研究区 广东省雷州半岛的纪家林场位于(20°54′N，109°52′E)，多年平均降雨量在1 500 mm以上[26]。海拔70 m，年平均气温在21～26 ℃之间，属于海洋性季风气候，土壤为玄武岩质黏土[26]。气象资料的获取来源于寒区旱区科学数据中心，其时间分辨率为3 h，包含近地面气温、近地面气压、近地面空气比湿、近地面全风速、地面向下短波辐射、地面向下长波辐射、地面降水率7个要素。该站点的桉树种植是从幼苗到轮伐期，而中国的桉树的轮伐期一般为6～7 a[27]，结合获得的实测数据，选取1996年6月至2002年6月一共6 a时间作为纪家站点的模拟时长，选取每3 h作为一个模拟步长。

根据张宁南[26]关于雷州半岛桉树人工林耗水量的研究，将土壤分为4层外加一层地下水层。残余土壤含水量、孔隙率、土壤饱和水力传导度是根据张宁南[26]研究中提供的土壤颗粒级配资料推算；地下水位是根据花利忠等[28]在雷州半岛桉树研究中所提供的资料，主要土壤参数如表2所示。该地区主要的桉树品种包括尾叶桉(Eucalyptusurophylla)、刚果12 W5(EucalyptusABL12)[28]。此次选取尾叶桉为研究树种，同样采取公式(8)的方法计算桉树和草的根系分布。

表2 纪家站点主要土壤参数

2 结果与分析

将2个研究区的参数数据分别代入耦合模型得出模拟结果，与已有的实测数据进行对比。选取蒸散量数据验证模型对桉树用水的模拟；选取林分材积、树高、胸高直径、植株密度验证模型对桉树生长情况的模拟。

2.1 Tumbarumba站点

将2001—2014年蒸散量模拟值和实测值进行对比(图1)。耦合模型模拟了蒸散量变化趋势。模型模拟的2001—2014年均蒸散量为819 mm，实测的年平均蒸散量为867 mm。将2001年年内蒸散速率模拟值与实测值对比(图2)，模拟蒸散速率从1—6月呈下降态势，然后直到12月都呈上升态势；实测蒸散速率从1—5月呈下降态势，之后除了8月略有下降以外，呈上升态势。Tumbarumba地区一般每年的5—10月为冬季，11月至次年4月为夏季，并且降雨大多发生在冬季。降雨增多，可供植被利用的水量增多，植被蒸腾和土壤蒸发都会增加，因此蒸散速率模拟值和实测值都是从5月左右开始上升。将胸高直径、林分材积、树高、植株密度实测值与模拟值对比，分析桉树生长模拟情况(表3)。除了树高的模拟精度相对较低，其余3个参数的模拟精度都较好。

图1 Tumbarumba站点2001-2014年

以2001年模拟结果为例进一步分析Tumbarumba站点桉树年内用水方式(图2)。列出2001年Tumbarumba站点月蒸散量模拟值、降雨入渗值、浅层土壤层(模型中设置的第1—4土壤层)含水量变化量之和、深层土壤水或地下水吸水量(表4)。在Tumbarumba站点，模拟结果显示，桉树除了使用土壤水，地下水也可能是桉树的重要水分来源，这和Keith等[29]在研究Tumbarumba桉树林的碳循环时发现桉树可能使用地下水的证据相符。耦合模型可以较为准确的模拟出Tumbarumba站点桉树的用水和生长情况。

表3 Tumbarumba站点树高、胸径、林分材积、植株密度模拟值和实测值对比

表4Tumbarumba站点桉树2001年年内用水方式mm

月份蒸散量入渗量浅层土壤水变化量深层土壤水或净地下水吸收量1140.3521.49-20.2898.572124.5369.4630.1685.23386.9929.0628.3086.23453.1327.7742.2967.65535.729.0437.9664.63622.0873.0695.5044.52723.7827.5849.9745.55828.7157.6668.3239.36956.4946.5928.0937.991067.9579.7048.1836.4211103.3932.56-27.5843.2312117.299.59-31.3476.35

图2 Tumbarumba站点2001年内每月蒸散速率模拟值与实测值对比

2.2 纪家站点

纪家站点并不具备Tumbarumba站点类似的长时间序列观测资料，因此选取1999年10月至2000年9月月蒸散速率的实测值和模拟值进行对比，如图3所示：模拟的月蒸散速率从1999年10月至2000年2月呈下降，然后2—8月呈上升趋势，之后略有下降。实测的月蒸散速率从1999年10月至2000年2月之间除了11月有提升，整体呈下降趋势，然后直到2000年7月都呈上升，之后两个月下降。模拟和实测值变化情况基本一致。和Tumbarumba站点处一

样，纪家站点的蒸散速率变化模式和降雨密切相关。纪家站点大致4—9月为雨季，在此期间降雨增多，这和蒸散速率模拟值以及实测值变化情况大致吻合。选取林分材积、树高、胸径、植株密度几个参数作为桉树生长模拟情况对比依据(如图4所示)，从拟合结果来看，除个别点有一定偏差外，总体拟合较好。以2000年模拟结果为例进一步分析纪家站点桉树年内用水方式。2000年雷州半岛纪家站点每月蒸散量模拟值、降雨入渗值、浅层土壤层(模型中设置的第1—4层土壤层)含水量变化量之和、深层土壤水或地下水吸水量如表5所示。

图3 纪家站点1999年10月至2000年9月每月蒸散速率模拟值与实测值对比

图4 纪家站点材积、树高、胸径、植株密度模拟值与实测值对比

不同于Tumbarumba站点，纪家站点模拟的桉树2000年全年几乎不使用深层土壤水或地下水，依靠降雨入渗量和浅层土壤水储备的前期含水量即可维持其蒸散需要。这和张宁南等[26]提出该地降雨量超过1 200 mm时桉树将不会使用地下水的结论一致。可见，耦合模型可以较为准确地模拟出纪家桉树林的用水和生长情况。

表5 纪家站点桉树2000年年内用水方式 mm

3 讨论与结论

原有的2个模型即3-PG模型[17]和Gou与Miller[12]提出的生态水文模型都是较为成熟的模型，但前者侧重于模拟植被的生长过程，对植被水分利用过程的模拟则相对粗略；后者侧重于模拟植被的水分利用过程以及相关生态水文过程，忽略了植被生长的影响。本研究将这2个模型进行耦合，考虑桉树生长与桉树用水之间的相互影响，耦合模型能较好地模拟Tumbarumba站点和纪家站点的生长和用水情况。

在Tumbarumba站点，如果单独使用Gou和Miller[12]提出的生态水文模型模拟桉树用水，模拟的2001—2014年年均蒸散量为661 mm；耦合模型模拟的年均蒸散量为841 mm，而实测资料的年均蒸散量为866 mm，耦合模型能更好地模拟蒸散量。如果单独使用3-PG模型模拟桉树的生长，则胸高直径的平均模拟精度由之前的71%提高到95%；林分材积的平均模拟精度由93%降低为88%；树高的平均模拟精度由88%提高到90%；植株密度的模拟精度不变。耦合模型对胸高直径的模拟精度降低了24%，其他参数模拟精度变化不大。

在纪家站点，如果单独使用Gou和Miller[12]提出的生态水文模型模拟桉树的用水，模拟的1999年10月至2000年10月期间蒸散量为1 236 mm，耦合模型模拟蒸散量为1 067 mm，而实测资料的蒸散量为1 069 mm，耦合模型能更好地模拟蒸散量。如果单独使用3-PG模型模拟桉树的生长，则胸高直径的平均模拟精度由之前的94%提高到95%；林分材积的平均模拟精度不变；树高的平均模拟精度由之前的95%降低到89%；植株密度的模拟精度不变。耦合模型对于桉树生长状况的模拟相比已有生长模型模拟的精度大体不变。各参数模拟精度都较好，因此耦合模型在纪家站点模拟表现较好。

总体而言，耦合模型在2个站点模拟桉树用水方面均有提高，但在模拟桉树生长方面没有明显提升，在Tumbarumaba站点胸高直径参数模拟精度还出现了下降。原因可能是因为Tumbarumba站点所在的巴戈国家森林(Bago State Forest)面积较大，不同区域获得的生长指标实测资料存在着一定的差异。而纪家站点桉树林研究区范围相对小得多，而且不同于Tumbarumba处桉树天然林，此处均是人工控制管理的人工林，林分组成更加一致。因此在有关生长指标的实测资料获取上Tumbarumaba站点误差可能会更大。耦合模型能较为合理地模拟桉树用水和生长情况，可以为后期桉树的科学种植于管理提供一定的理论帮助。但耦合模型模型参数众多，获取模型参数较为困难，后续模型的改进可以考虑从参数的简化方面入手。