最近10年来，随着“人口剧增、粮食短缺、资源危机、环境恶化”等全球性问题的日益严峻，农林复合系统受到众多发展中国家和地区的普遍关注和高度重视，而且欧美一些发达国家也都非常重视农林复合系统的研究与开发(Buck, 1995; Mead, et al., 1980; Garrett, 1997; Zinkhan, et al., 1997; Herzog, 1998)。农林复合系统因农林争水现象的存在，故水分生态特征问题一直是水资源紧缺地区农林复合系统的重要研究内容之一，主要涉及蒸散耗水、土壤水分、根系吸水等水分生态因子的时空变化规律、影响机理及其与单作作物系统的差异特征等问题。由于蒸散耗水、根系吸水和土壤水分的运动变化是十分复杂、互为耦合的物理过程和生理过程，因此，要想阐明复合系统的水分生态特征问题，必须全面考虑植物、大气和土壤的综合作用。

SPAC水分运移模型的研究一直是农学、气象、水文学、土壤学、生态学、自然地理学等相关学科共同关注的重要课题，也是当前国际水文计划(IHP)、世界气候研究计划(WCRP)及全球水量与能量平衡计划(GEWEX)等国际性前沿研究项目的重要内容之一。目前，国内外相关学者做了大量的工作，为农作物生产模拟、水资源的利用与评价及节水农业的发展，提供了重要的理论指导作用，但在森林生态系统中，利用数学模型描述SPAC水分运移过程，尚处在半经验半理论阶段(马履一, 1997)。农林复合系统，因其结构的特殊性和各类物理、生物参数的处理复杂性，系统地模拟研究其水分运移过程至今未见详尽的文献报道。

本研究以果农复合系统为例，根据SPAC水分传输理论原理，充分考虑农林复合系统的特殊性和复杂性，建立太行山低山丘陵地区农林复合系统SPAC水分运移模拟模型，并利用实测数据对模型进行验证，旨在为该地区农林复合系统水分生态特征的分析提供重要的方法和途径，并对其它非均匀下垫面条件下SPAC水分运移模型的研究提供借鉴。

1 试验区概况及试验设计

试验地设在河南省济源市裴村“太行山低山丘陵区复合农林业综合研究”试验区内(35°11 ′N, 112°03′E)。试验区内各类农林复合模式总面积达160 hm²，主要为果农间作。试验区地处太行山南段南麓。温带大陆性季风气候，全年日照时数为2 367.7 h，历年平均降水量641.7 mm，主要集中在7—8月份，占全年的44.3%。土壤以石灰岩风化母质淋溶性褐土为主，土层厚度：50~80 cm；pH值：6.8~8.5；石砾含量为10%~18%，有机质含量在1.0%左右, 速效氮、速效磷、速效钾分别为21.4~80、5.4~16、50~103 mg·kg^－1。

试验于1998年冬小麦拔节期~腊熟期(4月1日~6月4日)在水平梯田条件下的苹果(Malus pumila)-小麦(Triticum aestivum)复合系统内进行。梯田南北宽36 m, 东西长200 m。土壤质地：沙壤－轻壤，土层厚度为80 cm，面积为0.72 hm²，有灌溉条件。苹果树栽植于1992年，品种为新红星(Starkrimson)。果树带行向为东西向，共9行，株行距为3 m×4 m。树盘直径为80 cm，株高2.1 m，南北冠幅1.4 m；冬小麦行距20 cm，品种：“温8”号，播种期：1997年9月25日，基本苗248.0万株·hm^-2。在复合系统的中部选择一代表性较好的间作通道，沿垂直于行带方向，距果树带南北两侧各0、80、160 cm及间作中心处(分别简称“S0”、“S80”、“S160”、“N0”、“N80”、“N160”及“SN200”, 下同)布置土壤水分观测点，于S0、S120、N0、S120及SN200处布置土壤水分运移参数土样采集点，于S80、N80及SN200处布置小气候观测点，于S60、S120、N120、N60及SN200处布置小麦吸水根根量观测点，于S0、S80、S160、N0、N80、N160及SN200处布置果树吸水根根量观测点。

1.1 小气候观测

除利用试验区所在地气象站日际资料外，还选择典型晴天日，于各观测点小麦活动面上0.5 m高度处及地面高度2 m处采用采用通风干湿球表和轻便风速表进行温、湿度和风速的观测，于小麦活动面上1 m高度处采用CN-11辐射平衡表测定冠层净辐射。采用热通量板测定地表土壤热通量。每发育时段内观测1~3 d。测定时刻：6:30、19:30，6:30~19:30期间每整点观测1次。

1.2 土壤水分运移参数 1.2.1 土壤水分特征曲线

采用压力板膜法测定土壤水分特征曲线(华孟，1992)，按 Van Genuchten(1980)方程形式进行拟合。

1.2.2 饱和导水率(K_s)和非饱和导水率(K)

田间采集原状土，采用室内定位水头法测定K_s(华孟，1992)；根据K=K_s×e^bh的方程形式拟合得到k(b为拟合参数)。

1.3 土壤水分观测

测定深度0~80 cm。S0、N0处各土层土壤水分采用水银张力计测定，以及其他各测点除70~80cm土层土壤水分采用水银张力计测定外，其余各层次均则采用土钻烘干法测定，每隔10cm取样测定。

测定时期：张力计法每隔2 d测定1次，土钻烘干法每隔10 d测定1次，雨后均加测。

1.4 树干径流和冠层截留的观测

利用胶泥，将半剖面的圆形塑料胶管粘贴在处理平滑的树干0.5~1.0 m高处，呈“S ”型绕树体3~4周后，将胶管引流至一个体积为1.0 m³的塑料桶内，用以测定树干径流量；采用MS-Ⅰ型雨量计分别测定对照地的降雨量和林冠下降雨量，利用二者的差值表示冠层截留量。

1.5 小麦与苹果叶面积、根长密度及株高的测定 1.5.1 叶面积测定

苹果与小麦叶面积均每隔10~15 d测定1次。选取4株代表性较好的苹果植株，采用CI-110冠层分析仪测定苹果叶面积指数。并从每株中抽取上、中、下层若干标准叶片，带回室内，先将样本叶片描写在坐标纸上，再估算总面积，最后求取单株叶面积。在试验小区内随机抽取15株小麦，带回室内，量取叶长和叶宽，然后采用下式(康绍忠等, 1994)计算某单株叶面积：

(1)

式中，A:为某单株总绿色叶面积，n:某植株总叶片数，L_i:叶长度，W_i:叶宽。

1.5.2 根长密度测定

假设苹果树根系在以树干为中心、在南北方向上的分布具有对称性，采用剖面挖掘法，分层取样。清洗土壤和杂质后，以根直径小于等于1mm为吸水根的界限(马秀玲等, 1997)，进行分类。于小麦拔节初期和乳熟末期各取样1次；小麦根样采用根钻法于0~80 cm土层内每隔10 cm分层取样，样品重复数3个，清洗土壤和杂质。测定时期：3月14日、4月1日、4月15日、4月29日、5月8日及5月26日，分别正值冬小麦起身期、拔节期、孕穗期、扬花期、灌浆期和乳熟期。

苹果与小麦根长(LR)均利用交叉法(薛德榕, 1985)或直接测定，再将LR除以土样体积可得根长密度。

1.5.3 株高测定

于3月10日—6月1日期间，采用直尺测定苹果、小麦株高，每隔10天测定1次。

2 概念模型构建

综合模拟模型包括3个互为耦合的子模型：蒸腾蒸发子模型、土壤水分运动子模型及根系吸水子模型。

2.1 蒸腾蒸发子模型

蒸腾蒸发模型包括小麦和果树蒸腾子模型、林地及小麦棵间土壤蒸发子模型。

2.1.1 蒸腾(TR_a)

(2)

式中，TR_a苹果或小麦实际蒸腾；TR_p为苹果或小麦潜在蒸腾量，其计算模式可分别见文献(张劲松等, 2001, 2004)；K_s为土壤水分胁迫系数，可采用Jensen(1971)计算公式确定：

(3)

式中：θ、θ_f、θ_w分别为小麦或苹果根层土壤实际含水量、田间持水量和凋萎湿度。

2.1.2 土壤蒸发(EV_a)

(4)

式中，K_su为表层0~10 cm土壤水分胁迫系数，计算公式同式(3)。EV_p为土壤潜在蒸发，本文根据康绍忠等(1994)、余新晓¹⁾提出的蒸腾与蒸发比例关系计算原理，将TR_p换算为EV_p。计算公式如下：

1) 余新晓.晋西黄土地区防护林生态系统土壤动力水文学理论及其应用研究.北京林业大学博士学位论文.1994

(5)

式中，K为小麦或苹果消光系数，LAI为小麦或苹果叶面积指数。

2.2 土壤水分运动模型

由于立地条件为水平梯田，故研究时忽略了地表迳流，并假设果树行带方向土壤物理属性、生物属性及小气候属性具有均一性，并忽略土壤温度对土壤水分运动的影响。取西果树行带东方向为y轴方向、北南方向为x轴方向、土壤垂直向下方向为z轴方向，取北侧果树带中心(x=0)至南侧果树带中心(x=L_x)为计算区域(L为带距，即林带行距)。根据Richard非饱和流方程原理，建立二维土壤水运动方程：

(6)

式中，C：比水容重；h：土壤水势; K：非饱和导水率; SWCR、SWTR分别为作物、果树根系吸水量；Z：土层深度；X：果树带距，t：时间值。

2.3 根系吸水模型

根系吸水模型的数学表达式，可视为土壤水分运动方程的变换式。即将式(6)中的根系吸水项移至方程的左项中，将水势的时间变化率移至方程的右项中，则可得到根系吸水模型的数学表达式：

(7)

根系吸水模型的关键在于根系吸水模式的确立，借鉴一维宏观半经验半理论吸水模式，根据根系的分布特征，建立二维吸水模式。

2.3.1 作物吸水模式

采用Silem吸水函数形式²⁾，表达式如下：

2) 龚元石.以土壤水为中心的农田水量平衡模型及其应用研究.北京农业大学博士后出站报告.1994

(8)

(9)

式(8)、(9)中，TRCR_q为作物潜在蒸腾量；SWCR_p、SWCR分别为作物根系潜在吸水量、实际吸水量；RDCR为吸收根根长密度，LC_z为作物根区深度，其他符号同式(6)。

2.3.2 果树吸水模式

梯田土层较浅，80 cm以下为石砾，而果树根系生长相对较深，部分根系可以延伸至石砾层。由于石砾层土壤水分和土壤非饱和导水率的测定较为困难，所以，本研究中果树根系吸水模式采用Delong吸水函数形式¹⁾，表达式如下：

1) 龚元石.以土壤水为中心的农田水量平衡模型及其应用研究.北京农业大学博士后出站报告.1994

(10)

(11)

式(10)、(11)中：TRTR_p为果树潜在蒸腾量，RDTR为果树吸收根根长密度；SWTR_p、SWTR分别为果树根系潜在吸水量、实际吸水量；L_z为果树根区深度; L_x为间作带距(林带行距)。其他符号同式(6)。

3 模型的求解 3.1 定解条件

将以上各式联立成方程组，以土壤水分运移为中心进行求解。求解土壤水分运移方程(式6)时，还需列立初始条件和边界层条件。

3.1.1 初始条件

(12)

3.1.2 边界条件

水平方向上、下边界条件：假设果树带中心为不透水层，即当x=0或x=L时，

(13)

垂直方向上、下边界条件：上边界条件：

式(14)中，R、JL和IR分别表示降雨量、树干径流量和灌溉量，可实测得到; IN_ap为果树冠层降雨截留量，可实测得到；IN_w为小麦冠层降雨截留量，利用下式(刘昌明，1992)计算得到：

(15)

EV_w、EV_ap分别为小麦棵间土壤蒸发量和果树棵间土壤蒸发量，由棵间蒸发子模型(式5)计算得到。

下边界条件：

(16)

3.2 求解方法与求解过程

由于式(6)是非线性方程，所以对它的求解需采用数值计算方法。本研究中的土壤水分流动具有二维性质，为此，差分格式采用交替隐式差分(ADI)，并采用迭代法进行线性化处理，前后两次迭代相对误差最大值取0.01，采用追赶法求解所形成的线性化方程组(雷志栋，1988)。具体求解过程如下：首先，将计算区域(0≤x≥L_x, 0≤z≥L_z)按矩形网格离散化(见图 1)，沿x方向的节点号为i=0, 1, 2, …m，步长Δx=40 cm；沿z 方向的节点号为j=0, 1, 2, …，n, 步长Δz=10 cm；将时间离散化，节点号为k=1, 2, …t，步长Δt=10。然后，对控制方程和定解条件进行交替隐式差分。

3.2.1 z方向上隐式差分、x方向上显式差分

对于计算区域内任一节点(i, j, k+1)处相应偏微分方程(式6)的差分方程可表示为：

(17)

令:

将r₁、r₂、r₃和A_{i, j}、B_{i, j}、E_{i, j}、F_{i, j}代入式(17)，经整理，并将式(17)简化成：

(18)

根据z方向的边界条件，补充以下差分方程：

(1) 当j=0时, 对式(14)进行前向差分，并经整理可得：

(19)

令:

将B_{i, 0}, E_{i, 0}和F_{i, 0}代入式(19), 可得：

(20)

式(20)即为方程组(式18)中第1个方程。

(2) 当j=n-1时, 式(18)为：

(21)

其中:

(22)

式(22)是方程组(式18)中的最后1个方程式。

由于方程组左段含有时段末(k+1)的三个相临节点的水势，因此，根据时段初(k)已知的水势求出时段末(k+1)的各节点的水势，必须联立求解代数方程组。按差分方程(式18、20、21)表述的方程组如下：

(23)

将式(23)简化写成：

(24)

其中，[A]为系数矩阵，[F]为常数项阵，[h]^k+1为求解未知量距阵。[A] 仅在主对角及相邻两侧对角线有非零元素，所以，求解方程组为三对角型方程，可用追赶法求解，由消元和回代两个过程组成(雷志栋，1988)。

求解代数方程组式(23)或式(24)时，系数矩阵[A]和常数项阵[F]中均含有隐含值C_{i, j}^k+1和K_{i, j}^k+1，而C_{i, j}^k+1、K _{i, j}^{k+ 1}本身就是未知量h_{i, j}^k+1的函数，所以，对求解方程组需要进行线性化处理。常用的方法有显式线性化、预报校正法和迭代法等(雷志栋，1988)。本研究采用迭代法，并将允许误差设置为0.01。另外，求解过程还存在节点间水分运动参数(如K_j+0.5, K_j-0.5 )如何处理问题。本研究在空间节点上取上下节点间的几何平均值，在时间节点上取前后节点间的算术平均值。

消元和回代由i=0开始，解出h_{0, 0}, h_{0, 1} h_{0, n-1}, 依次解出i=1, 2, ……沿z方向各节点的时段末h值，直至i=m。一个时段计算后，改用z方向上显式差分、x方向上隐式差分的差分方程进行下一时段的计算。

3.2.2 z方向上显式差分、x方向上隐式差分

当z方向上显式差分、x方向上隐式差分时，同理，可将式(6)差分方程写成以下形式：

(25)

其中:

与上一时段相似，从j=1开始，由边界(到i=0和i=m)条件与式(25)组成代数方程组，采用追赶法依次求出j=1, 2, …n-2, n-1的x方向上各节点时段末水势。交替使用上述两种差分格式，便可求出所要求的各个时段的水势值。研究中采用visual Basic(6.0)和Access(7.0)混合编程。

4 模型验证

本文利用实测土壤水分资料对所建立的模型进行验证，结果表明：土壤容积含水量模拟值与实测值吻合效果较好(见图 2，表 1)，168对模拟值与实测值相对误差的平均值为6.89%，最大值为16.19%，最小值为0.09%，线性相关系数r可达0.836 7(n＝168), t＝0.197<t_0.01³³⁴＝2.568, 说明模拟值与实测值无显著性差异，线性拟合线方程为：

(26)

式(26)中，X和Y分别表示土壤容积含水量的模拟值和实测值(%)。

关于根系吸水量和实际蒸腾量模拟结果的验证，因试验条件的限制，实测值的获取比较困难，但由于土壤水分运移、植物蒸腾耗水与根系吸水本就是互为耦合的水分传输过程，所以，土壤水分的验证结果间接地验证了根系吸水模型和蒸腾模型的合理性。因此，可以利用所建立的模型，模拟计算复合系统中作物与果树的日蒸腾耗水、根系吸水量和土壤水分的时空变化量。

5 结论与讨论

根据土壤-植物-大气连续体(SPAC)水分传输理论原理，考虑农林复合系统的特殊性和复杂性，以果树株行距为3 m×4 m、东西行向的苹果-小麦复合系统为具体研究对象，建立太行山低山丘陵地区农林复合系统SPAC水分运移模拟模型, 该模型具有二维空间属性。

采用数值计算法求解模型方程，利用田间土壤水分实测数据对模型进行验证，结果表明：模拟值与实测值吻合效果较好，说明模型具有可行的计算精度。可以利用所建立模型对苹果-小麦复合系统土壤水分、蒸腾耗水、根系吸水等水分因子进行时空动态模拟计算，能为农林复合系统水分生态特征的定量研究提供可行的方法，并对其他非均匀下垫面条件下SPAC水分运移模型的研究具有一定的借鉴意义。

本模型是针对太行山低山丘陵地区特定的自然条件和特定结构的农林复合模式而建立的，对于其它结构(如不同株行距、不同行带方向)条件下水分因子的模拟计算，有关参数则需进一步调整。