森林作为最大的陆地类生态系统，能够涵养水源、拦洪蓄洪、减少地表径流中的泥沙含量、改变局部地区小气候以及改善水质，并且在水文过程中发挥着重要的调节作用^[1-3]。林冠层位于森林的最顶层，是与降水直接接触的第一个部分，在对降水的分配过程中具有先导作用。降雨通过林冠层后将被分为穿透雨、树干茎流以及冠层截留3部分，降雨强度及分布状况均发生改变，最终体现在对降雨的"能"和"量"的调节上，从而对森林生态水文过程产生重要影响^[4]。

冠层截留是评价森林水源涵养功能的一个重要指标，但由于不同立地的林分结构和气候条件差异较大，影响因素众多，如何准确测量冠层截留成为近年来国内外学者的研究热点之一。自1919年Horton提出了冠层截留的概念后^[5]，除实际测量外，很多学者也竞相建立起不同的模型对冠层截留进行模拟和预测，主要分为理论模型、经验模型和半理论半经验模型3大类，比如Horton模型^[5]；Rutter模型^[6-8]；Gash模型^[9]；以及我国学者建立的刘家冈模型等^[10]。其中被应用最广泛的模型为半经验半理论的Gash模型，半经验半理论模型相较于理论模型和经验模型的优点在于简化理论模型中过于繁琐的推导和计算过程并增加经验公式中所没有涉及到的截留过程，同时又保留二者的优点，具有较好的适用性^[11-12]。Gash模型是在Rutter模型的基础上进一步发展而来的，Rutter模型在实际应用中由于测量和计算的气象要素多而杂，极不方便，Gash于1979年在Rutter模型的微气候原理的基础上对其进行简化，之后经过Gash及多位学者的修正，最终形成了现在流行的Gash修正模型。Gash修正模型在国内外已经有大量的应用，徐丽娜等^[13]应用修正Gash模型对长白山地区的白桦林和云冷杉针叶混交林的冠层截留进行模拟以评价模型的适用性；何常清等^[14]在岷江上游亚高山川滇对高山栎林的林冠截留量进行模拟，并定量化该种森林类型在应用Gash模型时的参数值；曹光秀等^[15]基于修正后的Gash模型模拟亚热带常绿阔叶林的截留过程，为当地森林生态水文研究提供科学依据。一些国外学者譬如Murakami^[16]、Fathizadeh等^[17]和Motahari等^[18]也应用Gash修正模型分别对日本柏树幼林、伊朗半干旱地区扎格罗斯森林中的布兰特栎以及半干旱气候区松柏人工林的截留过程进行模拟。

笔者应用修证后的Gash模型，对青海省高寒地区的2种典型针叶林(华北落叶松和青海云杉)的截留过程进行模拟，旨在了解同一立地条件不同林分的降水分配特征，同时又能够为森林的水源涵养综合能力评价提供科学理论依据，维持林分稳定性，以充分发挥该地区林分的水土保持功能和水源涵养功能。

1 研究区概况

研究区在青海省西宁市大通回族自治区(E 100°51′~101°56′，N 36°43′~37°23′)闇门滩流域，该地区位于青藏高原和黄土高原的过度带，地势自东南到西北逐渐升高，海拔2 280~4 600 m，属于大陆性高原气候；全年降雨量由东南至西北部呈增加趋势，全年降水量达450 mm以上且分布极不均匀，5—9月的降雨量占全年降雨量的85%左右，全年日照充足，年蒸发量在1 750 mm左右。该地区主要的造林树种为青海云杉(Picea crassifolia)和华北落叶松(Larix principis-rupprechtii)。其中青海云杉为自然更新的顶级群落，而华北落叶松多为人工栽植养护。

2 研究方法 2.1 样地设置

林分样方使用手持GPS、指南针和坡度仪测定海拔、坡位、坡度、坡向并记录。样方按照20 m×20 m布置，4个边角使用喷漆做出标记，并用玻璃绳围起。从样地左下角作为起点呈S型对样地内树木进行编号，使用测高仪、胸径尺分别测量样地内每棵编号树的树高和胸径并记录其平均值；每个样地内选取5棵标准木，使用生长锥由南北方向钻入树干，将获取的随心带回实验室泡水后数出其年轮数量，估算林龄。样地具体信息见表 1。

2.2 观测数据的获取

1) 林内降水数据获取。在所选典型植被类型的样方内林下，采取6×6点阵式布设简易雨量桶，每个雨量桶间隔3 m，以观测林内穿透雨量。在每个样地内选取3株代表不同径阶的标准木，在每个典型林木的树干上缠绕一圈中部被剖开的PVC塑料软管作为导水槽，用U型针固定，并将软管和树干之间的缝隙采用玻璃胶进行密封，避免树干流无法全部流入塑料软管导水槽，下接简易集水器，以收集树干茎流。穿透雨和树干茎流同步观测，每次观测结束均采用量筒测量降水量。要想直接准确观测到林冠层截留量是非常困难的，通常采用的方法是使用其他数据间接计算得出，可用水量平衡法^[19]进行计算：

$ I = {P_{\rm{G}}} - {S_{\rm{f}}} - {T_{\rm{f}}}. $

式中：I为林冠层的截留量；P_G为林外降雨量，mm；S_f为林内穿透雨量，mm；T_f为树干茎流量，mm。

2) 常规气象参数测定。在林外空阔场地处及护林房布设便携式气象站(DAVIS公司生产，Vantage Pro2 Plus，美国)，该仪器采集数据的时间间隔为60 s，每30 min记录1次平均数据，该气象站测量的主要气象数据有太阳辐射、风速、温度、大气湿度、降水量等。

2.3 理论模型

在以往的研究中，平均蒸发速率(E)和平均降雨强度(R)往往采用多场降雨的平均值作为参数带入Gash模型中进行计算，这样会对模拟结果产生一定的误差^[20]，因此，笔者在引入郁闭度后的Gash模型的基础上，进一步做出如下调整，采用下式进行模拟：

$ \begin{array}{l} \sum\limits_{j = 1}^{m + n} {{I_{{\rm{aj}}}}} = c\sum\limits_{j = 1}^m {{P_{{\rm{G}}j}}} + \sum\limits_{j = 1}^n {\left( {{P_{{\rm{G}}j}} - P_{\rm{G}}^\prime } \right)} \left( {\frac{{{{\bar E}_j}}}{{{{\bar R}_j}}}} \right) + \\ c\sum\limits_{j = 1}^n {P_{{\rm{G}}j}^\prime } + q{S_{\rm{t}}} + {P_{\rm{t}}}\sum\limits_{j = 1}^{n - q} {\left( {1 - \frac{{{{\bar E}_{{\rm{cj}}}}}}{{{{\bar R}_j}}}} \right)} \left( {{P_{{\rm{G}}j}} - P_{{\rm{G}}j}^\prime } \right) \end{array} $

式中：m为未使林冠层持水量达到饱和的降雨次数；n为使林冠层持水量达到饱和状态的降雨次数；I_aj为第j场降雨的冠层截留量，mm；c为林分郁闭度；P_Gj为第j场降雨的降雨量，mm；P′_G为使林冠层达到饱和状态所需的最小降雨量，mm；E_j为第j场降雨对应平均蒸发速率，mm/h；R_j为第j场降雨的降雨强度，mm/h；P′_Gj为第j场降雨所对应的使林冠层达到饱和状态所需的最小降雨量，mm；q为产生树干茎流的降雨次数；S_t为树干的持水能力，mm；P_t为树干茎流系数；E _cj为第j场降雨的林冠平均蒸发强度，mm；E_cj=E_j/c，林冠饱和状态下的平均蒸发速率E_j可用利用Penman-Monteith方程进行模拟计算^[21]。

使林冠层达到饱和状态的最小降雨量可用下式计算：

$ P_{\rm{G}}^\prime = - \frac{{\bar R}}{{{{\bar E}_{\rm{c}}}}}{S_{\rm{c}}}\;\ln \;\;\left( {1 - \frac{{{{\bar E}_{\rm{c}}}}}{{\bar R}}} \right). $

式中：S_c=S/c，S为林冠饱和时的持水量，mm；E_c为林冠层平均蒸发强度，mm/h；E_c=E/c，E为林冠饱和状态下的平均蒸发速率，mm/h。

树干持水能力达到饱和状态所需的最小降雨量可用下式计算：

$ P_{\rm{G}}^{\prime \prime } = \frac{{{S_{\rm{t}}}}}{{{P_{\rm{t}}}}}\left( {\frac{{\bar R}}{{\bar R - {{\bar E}_c}}}} \right) + P. $

树干茎流量(S_a)和穿透雨量(T_a)可用下式进行计算^[22]：

$ \begin{array}{l} \sum\limits_{j = 1}^q {{S_{{\rm{a}}j}}} = c{P_{\rm{t}}}\sum\limits_{j = 1}^q {\left( {1 - \frac{{{{\bar E}_j}}}{{{{\bar R}_j}}}} \right)} \left( {{P_{{\rm{G}}j}} - P_{{\rm{G}}j}^\prime } \right) - q{S_t}\\ \;\;\;\;\;\sum\limits_{j = 1}^{n + m} {{T_{{\rm{a}}j}}} = \sum\limits_{j = 1}^{m + n} {{P_{{\rm{G}}j}}} - \sum\limits_{j = 1}^{m + n} {{I_{{\rm{a}}j}}} - \sum\limits_{j = 1}^q {{S_{{\rm{a}}j}}} \end{array} $

式中：S_aj为第j场降雨产生的树干茎流量，mm；T_aj为第j场降雨产生的穿透雨量，mm。

在应用Gash修正模型时，需要满足以下几个基本假设条件：1)林冠没有达到饱和状态时，不会产生穿透雨；2)只有林冠处于饱和状态才会产生树干茎流；3)只有在1次降雨结束的时候，才会产生蒸散量。2次降雨之间要有足够的干燥期，以使林冠完全干燥，本研究采用8 h为1次降雨间隔，即2次降雨间隔时间＜8 h时，视为一场降雨。

2.4 数据处理

使用Excel 2016进行数据整理、误差分析和敏感性分析，通过计算模拟值(predicted)和实测值(observed)的相对误差(MRE)来确定模型模拟的准确性，相对误差计算公式如下：

$ {\rm{MRE}} = \frac{{\left| {{\rm{observed}} - {\rm{predicted}}} \right|}}{{{\rm{observed }}}} \times 100\% $

使用Excel 2016进行敏感性分析，对Gash模型中的6个主要参数进行分析：E、R、S、c、S_t、P_t。每个参数的调整范围为-50%~+50%，调整步长为10%，并记录调整后的模拟截留量对比初始模拟截留量(即参数调整幅度为0时的截留量模拟值)的变化幅度。使用Origin 2018进行绘图。

3 结果与分析 3.1 林外降雨及蒸发特征

在观测期内(5月15日—8月15日)共检测到29次降雨，单次最大的降雨量为30.38 mm，出现于7月25日，单次最小降雨量为0.4 mm发生于6月11日；观测期降雨量共计255.83 mm；降雨总历时大约为155 h；平均降雨强度为1.95 mm/h，降雨强度最大为4.8 mm/h，出现于6月16日；最小降雨强度为0.52 mm/h，发生在6月8日。单次降雨量较少，强度较高，是该地区的降水特征。利用Penman-Monteith方程计算的平均蒸发速率如图 1所示。

3.2 林内穿透雨特征

华北落叶松与青海云杉的穿透雨与降雨量保持着良好的正相关关系，穿透雨量随着降雨量的增大而不断增大，分别占总降雨量的74.06%和71.03%，其关系式如下：

$ \begin{array}{l} 华北落叶松:{S_{\rm{f}}} = 0.79651{P_{\rm{G}}} - 0.49334\\ \;\;青海云杉:{S_{\rm{f}}} = 0.7433{P_{\rm{G}}} - 0.29704 \end{array} $

由图 2可见，华北落叶松在降雨量为0.4 mm的降雨场次中，穿透雨量为0，青海云杉在降雨量为0.4 mm和0.73 mm的降雨场次中穿透雨为0。在相同的降雨条件下，华北落叶松样地的穿透雨比例要略大于青海云杉样地。

3.3 林内树干茎流特征

通过建立树干茎流和降雨量的回归方程可发现，两者基本上呈线性正相关关系，可用下式表示：

$ \begin{array}{l} 华北落叶松:{T_{\rm{f}}} = 0.0028{P_{\rm{G}}} - 0.011\\ \;\;青海云杉:{T_{\rm{f}}} = 0.0175{P_{\rm{G}}} - 0.049 \end{array} $

华北落叶松和青海云杉的树干茎流占总降雨量的比例分别为0.15%和1.195%，相较于穿透雨量来说要小很多。树干茎流并不是一开始降雨就会产生的，而是具有一定时滞性，由图 3可以看出，落叶松在降雨量为8.62 mm左右开始产生树干茎流，而青海云杉在降雨量5.3 mm时产生树干茎流。树干茎流量和降雨量的拟合关系相比于穿透雨量较为分散。

3.4 林冠层截留特征

笔者利用水量平衡的原理，通过上文中所提到的穿透雨、树干茎流及降雨量的观测值，间接推算出林冠截留量，并建立其与林外降雨量的回归方程，如下：

$ \begin{array}{l} 华北落叶松:I = 0.20074{P_{\rm{G}}} + 0.50438\\ \;\;\;青海云杉:I = 0.23921{P_{\rm{G}}} + 0.33977 \end{array} $

如图 4所示，两树种林分的林冠截留量与降雨量均呈正相关关系，随着降雨量的增大，截留量也在不断增大；但其离散程度较穿透雨量—降雨量的回归方程较大。在整个观测期内，华北落叶松的总截留量为65.98 mm，单场次降雨的截留率在13.12%~100%，观测期内总的截留率为25.79%；青海云杉的截留量为71.05 mm，略大于华北落叶松，单次降雨的截留率为15.95%-100%，总截留率为27.78%。

3.5 模型参数确定

1) 使用Leyton^[23]提出的残差回归法进行计算。最终得出华北落叶松的冠层饱和持水量为1.342 6 mm，青海云杉为1.576 3 mm(表 2)，大于使用浸水法测得的华北落叶松(1.28 mm)和青海云杉(1.37 mm)的冠层饱和持水量。该现象与芦新建等^[20]和王军帅^[24]的研究结果相同，即利用回归法模拟出的林冠枝叶部分的持水能力会被高估。

2) 树干最大持水量和树干茎流系数一般采用回归法得出^[25-26]。华北落叶松和青海云杉的树干最大持水量分别为0.011 mm和0.049 mm；树干茎流系数分别为0.002 8和0.017 5。

3.6 模型模拟

笔者应用修正后的Gash模型分别对华北落叶松和青海云杉2个林分进行了截留模拟，具体如表 3所示。在整个观测期内，两种林分的模拟值均为穿透雨量＞冠层截留量＞树干流，与实测值有相同的趋势。其中，华北落叶松的模拟冠层截留量为59.02 mm，低于实测值6.96 mm，相对误差为10.55%；模拟树干茎流量为0.27 mm，低于实测值0.11 mm，相对误差为28.95%；而穿透雨量的模拟值则略大于实测值，相对误差为3.73%。青海云杉的冠层截流量为61.36 mm，比实测值小9.69 mm，相对误差为13.64%；树干茎流量的模拟值为2.41 mm，比实测值小0.64 mm；穿透雨模拟值(192.05 mm)比实测值(181.73 mm)高10.32 mm，相对误差为5.68%。Gash修正模型对于华北落叶松的林冠截留量的模拟效果要优于青海云杉。

3.7 敏感性分析

敏感性分析结果如图 5所示，对于华北落叶松，对截留量影响最大的参数是S，其次是E，然后是R，当S、E和R增加10%时，截留量的模拟值变动幅度分别为4.56%、3.84%和-3.5%；对于青海云杉，对截留量模拟值影响最大的3个参数也为S、E和R，当其分别增加10%时，模拟值变化幅度分别为6.07%、2.78%和-2.53%。对于2种林分而言，影响最小的参数均为c、S_t和P_t，当其分别增加10%时，模拟值的增加幅度均小于2%，各项参数的敏感性排序为S＞E＞R＞c＞S_t＞P_t，除R对结果产生负影响外，其余参数对截留量模拟结果均产生正影响。

4 讨论 4.1 林分的降雨分配特征

相关文献表明，我国主要的林分类型的截留率大约在11.4%~36.5%之间^[27]，芦新建等^[20]对高寒地区25次降雨中的华北落叶松截留率的测量结果为23.34%；高婵婵等^[12]对祁连山东部地区28次降雨中的青海云杉截留率的实测值为25.30%，与本研究实验结果接近。笔者发现，青海云杉的截留量略高于华北落叶松，但是前者产生的树干流(3.05 mm)却远大于后者(0.38 mm)，魏曦等^[26]在河北围场县林场的研究中测得华北落叶松的树干茎流为0.17~0.51 mm远小于油松的树干流0.64~2.81 mm，这说明华北落叶松树皮的形态特征不易产生树干茎流。笔者发现，华北落叶松的林分密度和郁闭度均小于青海云杉样地，这可能是造成前者树干茎流小于后者的原因。

4.2 模型模拟

在应用Gash模型的模拟中，2种林分的模拟值均小于实测值，这与曹光秀等^[15]、柴汝杉等^[28]分别对亚热带常绿阔叶林和红松林的研究结果一致，造成模拟值偏小的原因可能是：Gash模型应用的前提是假定每次降雨前林冠都达到了完全干燥的状态，本研究划分两次降雨的时间间隔是>8 h算1次降雨，但是在实际应用时，尤其是在青海省高寒地区，该地区昼夜温差大，平均气温较低，在8 h内可能无法完全达到干燥状态，这就导致在下次降雨时，使林冠达到饱和状态的降雨量会一定程度上减小，同时林冠层内的枝叶和其他生物会吸收一定量的水分，但是Gash模型无法模拟这部分水分，从而导致模拟值偏小；另一方面，林分中一些较高灌木对雨水的截留可能会影响穿透雨的测量，最后导致截留流量实测值偏高。赵洋毅等^[29]在对缙云山毛竹林的研究中，截留量的模拟值比实测值略高；王艳萍等^[30]在黄土区人工刺槐林的研究中也得出截留量模拟值比实测值高的结论，笔者分析，在Gash模型的模拟过程中，其认为当降雨量大于使林冠层达到饱和状态的最小降雨量时(P_G≥P′_G)才会产生穿透雨，然而在实际情况中，一些长历时、雨量小的降雨，在降雨开始时就会产生少量的穿透雨并伴随整个降雨过程，使得截留量的模拟值偏大。

4.3 敏感性分析

在敏感性分析中得出，对于两种林分的模拟结果影响最大的3个参数为E、R和S，其中降雨强度对模拟结果产生负影响，其余参数均产生正影响，这与李亦然等^[25]的研究结果一致。另外，高婵婵等人在对青海云杉的研究中，将平均蒸发速率和降雨强度的比值作为一项参数进行分析，并得出结论该项参数对模拟结果的影响最大。而曹光秀等^[15]、Limousin^[31]和时忠杰等^[32]的实验结果认为影响最大的参数为林分郁闭度c。造成这差异的原因可能是树种和立地条件不同所导致的。

综合来看，Gash修正模型在对华北落叶松和青海云杉的冠层截留模拟中，取得了较为合理的模拟结果，是否适合该立地条件下的其他树种有待进一步研究。由于本研究是利用Penman-Monteith方程对平均蒸发速率进行估算，模型本身会受到多种环境因子的影响，势必会产生一定误差，该方程的模拟精度在以后的研究中有待进一步评估；另外由于实验中采用的气象站测量频率是30 min记录一次数据，对降雨历时的测量也有一定偏差，使降雨强度的计算产生误差，导致最终的模拟结果偏离实测值。随着科学技术的发展，涡度协方差技术越来越成熟，它能够直接将原始的观测数据直接处理成科学研究所需的通量结果，从而避免模型本身具有的误差性，未来有望将其应用到Gash模型的模拟以及森林水文方面的研究中，使得模拟结果更加精确。

5 结论

1) 在观测期内的总降雨量为255.83 mm；该地区2种典型针叶树种的降雨分配特征为：华北落叶松的穿透雨率为74.06%，冠层截留率为25.79%，树干茎流占总降雨量0.15%；青海云杉的穿透雨率为71.03%，冠层截留率为27.78%，树干茎流占降雨比例为1.195%。华北落叶松林相比青海云杉林更容易汇集雨水。

2) 树干茎流量的产生相对于降雨量具有一定滞后性，华北落叶松和青海云杉分别在降雨量为8.62 mm和5.3 mm时产生树干茎流。

3) 穿透雨量、冠层截留量和树干茎流量均随林外降雨量增大而增大，并保持良好线性正相关关系；穿透雨率和降雨量呈正相关对数关系，林冠截留率则和降雨量呈负相关对数关系，

4) 修正后的Gash模型能够较好的模拟2种林分的冠层截留量。华北落叶松(林分密度2 700株/hm²，郁闭度0.62)和青海云杉林(林分密度3 125株/hm²，郁闭度0.92)的冠层截留量模拟值相对误差分别为10.55%和13.64%，华北落叶松的冠层截留模拟效果要优于青海云杉。

5) Gash修正模型中各项参数的敏感性排序为S＞E＞R＞c＞S_t＞P_t，除R对结果产生负影响外，其余参数对截留量模拟结果均产生正影响。