斑块状植被覆盖下坡面流水动力学特性

2020-05-19刘西乐张宽地王静雯

农业工程学报 2020年6期

刘西乐，张宽地,2※，杨婕，孟鹤，王静雯

（1.西北农林科技大学水利与建筑工程学院，杨凌712100；2.中国科学院水利部水土保持研究所黄土高原土壤侵蚀与旱地农业国家重点实验室，杨凌712100；3.河套学院，巴彦淖尔015000）

关键字：水动力学；坡面流；植被斑块；覆盖度；淹没度

0 引言

坡面水流受坡面地形地貌、土壤质地、植被特征以及降雨强度和历时等条件影响[1-2]，其阻力规律时空变化十分明显[3-4]，远比河道水流复杂。一般认为，植被具有壅高水深[5]，减缓流速[6-7]，增大阻力[8-9]，降低湍流强度[10]，促进沉积的作用[11]。目前，植被覆盖下坡面流水动力学特性研究大多从植被类型、覆盖度、刚度、径粗、淹没度及高度等方面入手，通过将一些植物学属性与水流阻力建立关系。Jarvela[12]采用天然水草、莎草以及柳枝研究不同植被类型对阻力系数的影响，指出柳树的叶子对阻力系数影响很大。赵璐等[13]认为坡面阻力系数并非随覆盖度增大而增大，而是存在一个理想的临界值。Zhang等[14]认为达西阻力系数随着植被茎干直径的增加而增加。Lawrence[15]基于淹没度的概念研究了水流阻力模型，得出部分淹没下水流阻力与淹没度和覆盖度呈正相关。孙菊英等[16]认为植被高度对阻力系数和雷诺数间的关系影响较大，并基于不同高度柔性植被试验得出裸坡或植被高15 cm时，阻力系数与雷诺数呈负相关，而植被高为4和8 cm时，阻力系数随雷诺数增大先减后增。近年来，植被空间位置[17]和植被的分布格局吸引了越来越多学者的关注[18]，在坡面尺度上表现尤为明显，植被空间格局的影响不可忽视[19]。不同于受人类活动影响的坡耕地、城市地表铺装绿地，由于自然条件下的坡地植被受人为因素干扰较小[20]，坡面植被大多处于非淹没状态[21-22]，且极易遭到破坏造成生长不好或分布不均，有植被处和无植被处受冲蚀的程度不同，呈高度破碎的斑块状分布[11,19]，这种簇状丛生、近似刚性的植被斑块对坡面流的流向和流量分配、流型流态及阻力规律有着深远影响[23]。因此，确定随机分布斑块状植被覆盖下的坡面流最佳水流条件对水土保持和水力侵蚀控制至关重要。

为此，本文对斑块状植被条件下的坡面流水动力特性进行了理论分析和试验研究，探讨了植被条件下坡面流流型流态的归属以及植被特征对阻力机制的影响，并基于水流阻力等效原则，建立水力坡度、覆盖度、雷诺数、淹没度4种影响因子作用下坡面流阻力计算公式，以期阐明斑块植被覆盖下坡面流阻力的内在规律，深入了解坡面水文过程和水力侵蚀过程的驱动因素。

1 理论分析

目前对于植被覆盖下坡面流阻力的计算大多基于传统明渠均匀流公式。故借鉴明渠流中研究植被阻力的推理方式，以蒋北寒等[24]非淹没工况刚性植被作用下阻力公式为基础，结合坡面流特征，推导植被条件下坡面阻力的计算模型。含植被坡面流中主要存在的力包括植被引起的绕流阻力FD、床面切应力Fb、坡面水体的有效重力FG。斑块状植被覆盖下坡面流植被段控制水体受力分析如图1所示。

图1 植被斑块覆盖下坡面流示意Fig.1 Diagram of overland flow under vegetation patch coverage

根据力学平衡原理，选取长度为2 m的植被区流段进行受力分析，其受力关系为

式中U*为摩阻流速，m/s；U为断面平均流速，m/s；ρ为水的密度，g/cm3；B1为有效过流宽度，m；L为植被铺设长度，本试验取2 m；f为达西—威斯巴哈阻力系数。

式中CD为拖曳系数；Ai为植被迎水面面积m2；N为草被株数，m；D为草被直径，m；h为断面平均水深，m。

式中B为水槽宽度，本试验取0.3 m；Cr为植被覆盖度。

式中R为水力半径，m；J为水力坡度，实验室条件下近似取sinθ；ν为水流运动黏性系数；Re为雷诺数；g为重力加速度，m/s2。

式中t为水温度，℃。

式中λ为与植被类型、配置相关的系数[25]。

将式(2)～式(5)代入(1)中，得：

将植被绕流阻力均匀的分布到床面上，得到等效过水断面E，如图2所示。等效过水断面一系列水力参数下标记为“e”。

图2 等效过水断面Fig.2 Equivalent water cross section

根据流体质量守恒定律可知，等效前后控制水体体积不变，植被段控制水体体积为

则等效平均水深he为

根据液流连续性方程可得：

等效床面剪应力可以写为：

对等效过水断面控制体受力分析，由力学平衡原理可得：

联立式(1)和(13)，可得：

从式(14)可以看出，等效过水断面床面剪应力Feb可以视为由原床面剪应力Fb和等效植被附加阻力FD两部分组成，符合摩阻作用的可加性[26]。

式(15)整理后可得非淹没条件下坡面流等效阻力系数：

式中β为淹没度，β=h/H，H为模拟草被平均高度，m。上式表明：f=F(Cr，Re，J，β，D，H，λ)，植被覆盖下等效阻力系数是覆盖度、水力坡度、淹没度、雷诺数、植株直径、植被高度、植被类型和植被空间格局的函数，基于此建立植被覆盖下等效阻力系数计算模型：

式中a为植株直径、植被高度、植被类型和植被空间的综合参数，在本试验对植被自身属性进行概化，模拟草被采用单一种类仿真水草，单颗植株直径、植被高度均取平均的定值，其中空间格局呈簇状随机分布，故a为常数值。b、c、d、e分别为覆盖度、坡度、雷诺数、淹没度的指数，其绝对值代表各自变化对综合阻力系数的影响大小。通过Nash-Sutcliffe系数来检验等效阻力系数计算模型的模拟效果，具体计算公式为

式中Pi为模拟值；Oi为实测值；-O为实测值的平均值；n为样本个数。

2 试验

试验在中国科学院水利部水土保持研究所黄土高原土壤侵蚀与旱地农业国家重点试验室进行，为消除床面形态等边界条件变化对试验结果的干扰，采用定床阻力冲刷试验。试验装置主要由供水装置、径流收集装置、试验水槽及流量和坡度调节装置四部分组成，试验装置结构如图3所示。

图3 试验装置结构示意Fig.3 Structure diagram of test device

供水装置为长2 m、宽0.5 m、高1 m的恒定水箱和扬程为12 m的潜水泵。试验水槽由长4.5 m、宽0.3 m、高0.1 m的坡度可调的有机玻璃板制成，坡度设置黄土高原地区常见坡度为0.034 9、0.069 8、0.104 5、0.139 2、0.173 6和0.207 9共6个坡度。黄土高原地区临界侵蚀雨强范围为30～36 mm/h，最大峰值雨强可达150 mm/h。单宽流量范围设置为0.278×10-3～2.500×10-3m2/s，相当于面积为8.33～90.00 m2产生的汇流。流量由流量调节阀控制，采用称重法测定。

本试验概化坡面植被类型[27]，采用人工矮草模拟黄土高原区旱生丛生的小禾草及旱生矮草组、半灌木组等优势种植物，采用人工加糙粒径为0.380 mm的水砂布模拟天然土壤，床面粗糙度由所粘贴砂粒粒径替代。人工矮草由60根平均直径1 mm的针状圆柱棒固定在直径20 mm薄片底座上组成一簇模拟矮草，草平均高20 mm。试验设置覆盖度Cr分别为0、1.26%、2.51%、3.77%和5.03%共5种工况。为模拟天然状态下植被分布形态，模拟草被采用斑块状随机布设，每个植被斑块由4～5簇人工矮草聚集而成，形状近似圆形，以消除因斑块形状的不同而对试验结果的影响。距水槽起始位置0.75 m处布设植被，沿水槽设纵向观测断面5个，沿坡面自上而下分别为0+1.00 m、0+1.50 m、0+2.00 m、0+2.50 m、0+3.00 m。每个断面横向设3个测点，均观测断面表层优势流速(KMnO4染色示踪法测流速)和水深(采用重庆水文仪器厂生产的SX402数显测针仪测定，精度为0.01 mm)。取断面0+1.00 m和断面0+3.00 m间水体作为控制水体，长度为2 m。

图4 不同植被覆盖度Cr下人工模拟植被斑块状分布俯视图Fig.4 Top view of patch distribution of artificially simulated vegetation under different vegetation coverage Cr

3 结果与分析

3.1 植被斑块覆盖下坡面水流机理分析

由于植被占据坡面部分面积，且斑块内和斑块之间存在一定间隙，植被条件下坡面流被植被划分成部分连续的空间，使得水流与植被有效接触面积复杂化，坡面流内在结构、流量分配及流向的确定、能量的传递与耗散极为复杂[28]。

小流量下，坡面水流尚未铺满整个坡面，水层很薄。由于植被的存在，水流时分时合，在斑块间呈股状行进，形成辫状交织的水网[29]。该工况下，每一斑块独立发挥阻水作用，迎水面会产生壅水，根据最小能耗原理，水流大部分绕斑块行进，且流动方向不唯一(图5a)。随着流量的增加，平均水深与平均流速随之增大，该工况下，斑块间阻水效果相互影响，水流紊动强度大，流线极为复杂，在斑块下游产生互相干扰的尾迹和折冲水流，当流量和坡度较大时，甚至会产生涡旋[30](图5b)。植被覆盖和裸露地带阻力不同，分配的流量、水深和流速产生较大差异。由于植被的阻碍，斑块内水流分配的流量相对较少，流速小、水位高，因此在重力作用下斑块内水流会产生横向径流[31]，水流从斑块两侧及下游流出，汇入斑块间主流。斑块间水流由于缺少植被阻碍从而获得较多的流量，流速大，水位高低不一。沿程汇合斑块内的流量后向坡度较低处流动时往往会撞向下游斑块，作为下游斑块内水流的流量与能量来源(图5b)。裸坡和植被斑块的组合形成了一个由源和汇组成的空间镶嵌图，这与在大空间尺度上观察到裸地通常携带径流和沉积物(源)，而植被则捕获水和沉积物(汇)的趋势极为相似[31]。

图5 植被斑块作用下坡面水流Fig.5 Overland flow under action of vegetation patches

3.2 随机斑块状植被覆盖下坡面流流型流态

坡面水流流型流态的差异是由水流内部结构决定的，不同的流动型态有着不同的阻力机制[32]。不同的流态使同一流量具有不同的水深及流速，从而使得相应的坡面流阻力难以确定。坡面流属薄层水流范畴，受外部及边界条件影响显著，尤其植被存在扰乱坡面流流动，致使坡面流流态极为复杂。图6绘制了不同坡度、不同覆盖度条件下，各试验工况水流流区分布。

本试验雷诺数变化范围为205～2 622，弗劳德数变化范围为0.224～2.574。由图5可以看出，各试验工况下，坡面植被水流均分布在虚拟层流区[29]和过渡流区，未达到紊流区。水流流态发育受覆盖度和坡度共同制约。坡度增加，水流流态由缓流向急流延伸，覆盖度增大，流态逆向发育，由急流向缓流延伸。且随着单宽流量的增加，水流的平均流速与平均水深逐渐增大[9]。

由图6a可知，当坡度一定，单宽流量由0.278×10-3m2/s增至2.500×10-3m2/s时，平均流速由0.102 6 m/s增大到0.545 1 m/s，平均水深由0.222 cm增大到2.415 cm，增大速率受制于覆盖度。裸坡条件下（即Cr=0时），平均水深与平均流速线性相关。随着覆盖度增加，相同流速增幅条件下，水深的增幅随着植被覆盖度的增加而逐渐增大，且流量和植被覆盖度均较大时，水深增加极为显著，而流速增加非常有限。其原因为植被的阻水作用，覆盖度越大，植被壅水现象越明显，对流速增加有一定抑制作用[33]；随着单宽流量增加，水深逐渐增大，坡面水流绕植被流动，产生壅水和相互干扰的尾迹[34]，当流量较大时，水流紊动强度进一步增大，流线较为复杂，坡面水流由虚拟层流区逐渐发育为过渡流区；相同水流条件下，坡度增加，水流趋向急流，而覆盖度越大，植被阻水效果愈明显，流态朝着缓流区发育。

图6 不同坡度J下坡面植被水流流区分布Fig.6 Distribution of vegetation flow zones on slopes with different slopes J

3.3 随机斑块状植被条件下坡面流阻力机制

3.3.1 颗粒阻力占比与单宽流量和覆盖度的关系

坡面流阻力的划分实质是根据紊动产生的根源对阻力来源进行合并归类。一般认为植被覆盖下坡面流阻力来源于下垫面，由水土界面的摩擦力[35]和由于植被存在使得水流内部质点紊动而引起的绕流阻力组成[28]不同于明渠，坡面流由于水深较浅，颗粒阻力的影响不可忽略[29]，因此有必要从阻力占比的角度分析坡面阻力机制。颗粒阻力占比用α表示，α=f裸坡/f，其中f裸坡为裸坡条件下坡面达西阻力系数。

不同坡度下单宽流量和覆盖度与阻力构成的变化规律基本一致，因此以J=0.034 9和J=0.139 2为例进行分析。由图7可以看出，相同单宽流量下，α随Cr增加而减小，减小程度随流量递增。同一覆盖度下，α随q增加而减小。究其原因，同一单宽流量下，Cr越大，植被占据空间越多，迫使水流具有相较于裸坡条件下较高水位，水流与植被有效接触面积增大，植被阻力占比增大。同时由于植被存在，水流与裸露坡面接触面积减少，颗粒阻力减小。植被阻力增加和颗粒阻力减小共同作用，呈现出α随覆盖度增加而减小的变化趋势；同一覆盖度下，随着单宽流量的增加，水位逐步上升，水流与植被有效接触面积增大，植被阻力增大，而颗粒阻力变化不大，故水流α随单宽流量增加整体上呈减小趋势。

图7 两个水力坡度下单宽流量和覆盖度对颗粒阻力占比的影响Fig.7 Influence of unit discharge and coverage on particle resistance ratio under two hydraulic gradients

3.3.2 平均水深对颗粒阻力占比的影响

如图8所示，植被覆盖下，α和h呈幂函数负相关，相关性随Cr增加而增大，与坡度基本无关。

由图8可以看出，同等水力条件下，Cr越大，α越小。不同覆盖度下，α和h均呈幂函数负相关，且相关性随Cr增大，R2由0.803增大到0.982。当Cr=1.26%时，颗粒阻力占比变化范围为30.30%～98.08%(图8a)，当Cr=5.03%时，颗粒阻力占比变化范围为2.47%～55.34%(图8d)。覆盖度较低时(Cr=1.26%和Cr=2.51%)，水流流动相对平缓，水流主要在斑块间流动，坡面裸露程度较高，水流与植被接触面积有限，颗粒阻力处于优势地位。当覆盖度较高时(Cr=3.77%和Cr=5.03%)，坡面裸露程度低，颗粒阻力较小，此时植被阻力占主导，随着水深增大，α仅为2.47%。

图8 不同植被覆盖度Cr下颗粒阻力占比α与平均水深的关系Fig.8 Relationships between particle resistance ratioαand average water depth under different vegetation coverage Cr

3.3.3 阻力系数与雷诺数的关系

由图9可知，当水流处于虚拟层流区时，综合阻力系数与雷诺数呈负相关变化；当水流处于过渡流区时，f－Re关系受覆盖度和坡度共同作用。随着覆盖度的增大，综合阻力系数与雷诺数的相关性由负变正，在2.51%和3.77%的范围内存在一个临界值，且覆盖度越大，不同坡度下的f－Re关系曲线逐渐分化，坡度影响凸显。

如图9所示，当Re＜580时，水流处于虚拟层流区，综合阻力系数随雷诺数增加而降低。当580＜Re＜6500时，水流处于过渡流区。裸坡条件下（即Cr=0.00%），综合阻力系数随雷诺数增大而减小，减小趋势比虚拟层流区要缓(图9a)。对裸坡条件下f－Re关系进行拟合，具体函数关系为：f=102.012Re-0.865(R2=0.853)，这与张光辉[36]试验得出的指数-0.72很接近，此工况下坡面阻力仅为粗糙元产生的颗粒阻力。大于临界覆盖度时，f－Re关系曲线与裸坡条件下呈现出明显差异，二者相关性由负相关逐渐变为正相关，坡度影响逐渐凸显，这与闫旭峰等[30]得出的结论较为一致，其原因与植被覆盖下坡面流阻力机制、内部紊动程度有关。坡面植被水流处于虚拟层流区时，其水层较薄，坡面阻力以颗粒阻力为主，流线较为平顺，紊动强度低。坡面植被水流处于过渡流区时，壅水现象明显，由于植被的阻滞作用，水流产生掺混消能效应，增大了水流阻力。该工况下，横向流动凸显，水流紊动强度大。相同覆盖度和上游来水情况下，坡度越小，植被阻水效果越明显，水层较大坡度下厚，植被阻力影响明显。

图9 不同植被覆盖度Cr下，达西阻力系数与雷诺数的关系Fig.9 Relationships between Darcy resistance coefficient and Reynolds number under different vegetation coverages Cr

3.4 淹没度对综合阻力系数的影响

坡面流为薄层水流，水深较浅，坡面植被处于非淹没状态[22]，植被在水流中的淹没深度决定着植被阻力大小。因此有必要从淹没度角度分析坡面阻力。

由图10可知，低于临界覆盖度时，整体上综合阻力系数随淹没度增大而减小。高于临界覆盖度时，综合阻力系数随着淹没度增大先减后增，坡度对其也有一定的影响。当覆盖度小于临界值时，淹没度较小工况下，综合阻力系数随淹没度增加急剧降低，淹没度较大时，综合阻力系数随淹没度增加反而呈现逐渐增大的趋势(图10a和图10b)。当覆盖度大于临界值时，综合阻力系数随着淹没度增加呈先减后增趋势[37](图10c和图10d)。其原因可能与水流阻力的构成有关。小淹没度下，水层相对较薄，主流在斑块间流动，水流贴坡流动表面张力较大，颗粒阻力起主导作用。随着淹没度的增大，植被有效阻水面积增大，植被阻力逐渐凸显。而覆盖度大于临界值时，坡面覆盖程度高，相同水力条件下水深较大，植被与水流接触面积较大，水流在斑块间与斑块内穿梭流动，斑块内横向侧流现象明显，此时水流紊动强度较大，植被阻力起主导作用。由图10d可以看出，f－β关系曲线随坡度减小整体向右上延伸。相同流量和覆盖度条件下，坡度小时，淹没度大，植被有效阻水面积同比增加，植被拖曳力凸显，综合阻力系数增大。

图10 不同植被覆盖度Cr下淹没度对达西阻力系数的影响Fig.10 Influence of submergence on Darcy resistance coefficient under different vegetation coverages

3.5 覆盖度对坡面阻力的影响

如图11所示，综合阻力系数与覆盖度呈现幂指函数相关特征，这与曹颖等[33]的研究结果一致。由图11a和11b可知，低于临界覆盖度时(0和1.26%)，阻力系数范围为0.692～3.037，综合阻力系数随Cr增加而增大的趋势较缓；高于临界覆盖度时(3.77%和5.03%)，阻力系数变化范围为0.118～4.759，综合阻力系数随Cr增加而增大的趋势较陡。究其原因，综合阻力系数的大小主要取决于植被有效阻水面积。Cr较小时，相同水流条件下植被排水体积小，植被与水流有效接触面积较小，植被对水流的拖曳力较小；随着Cr的增加，植被排水体积增大，有效接触面积增大，植被阻水能力增强。其次，由于表面张力存在，Cr较小时坡面水流尚未平铺整个坡面，大部分水流从斑块间穿梭，此工况下α较大，植被阻力有限。Cr增大后，流量重新分配，斑块间阻水效果相互影响，尾迹和折冲水流混掺，增加了能量的消耗[35]，综合阻力系数进一步增大。

图11 不同水力坡度J下达西阻力系数与覆盖度及单宽流量的关系Fig.11 Relationships between Darcy"s resistance coefficient and coverage and single-width flow at different hydraulic gradients J

此外，综合阻力系数与覆盖度的关系受制于流量和坡度，在不同坡度流量组合下，阻力系数随覆盖度变化趋势的差异越来越显著[33]。由图11a所示，随着覆盖度的增大，综合阻力系数与流量由负相关变为正相关。其原因可能与坡面阻力构成有关，水深随流量增大而增大，由前文可知α和h呈幂函数负相关，随流量增加α逐渐减小，颗粒阻力主导位置逐渐由植被阻力所代替。不同流量条件下f－Cr关系曲线会发生交叉，交叉点区域随坡度增加整体向大覆盖度方向延伸(图11)。随着坡度增大，f－Cr关系曲线逐渐放缓。这是由于综合阻力系数f随着坡度的增加而减小，随着坡度增大，f的降低趋势明显变缓[34]。

3.6 模型的建立

由已建模型可知，坡面等效阻力系数与覆盖度、淹没度、坡度、雷诺数有关，基于此，通过SPSS软件对168组试验实测数据进行回归拟合，得到不同流态下坡面综合阻力系数f的半经验计算公式：

式(19)和(20)中各因子指数绝对值的大小反映了其变化对阻力系数影响的大小，由此可以看出，各因子对综合阻力系数的影响与流态密切相关，当水流处于虚拟层流区时，雷诺数是影响综合阻力系数的主导因素，淹没度次之；当水流处于过渡流区时，覆盖度成为影响综合阻力系数的主导因素，雷诺数次之。同时采用Nash-Sutcliffe系数（NSE）对公式模拟情况进行校验，NSE值越接近1，说明模型模拟效果越好，计算得出层流区NSE值为0.836，过渡流区NSE值为0.784，因此该模型可以较好地模拟计算斑块状植被覆盖下坡面流综合阻力系数。