基于WEPP模型的工程堆积体不同堆置方式的水土流失效应研究

刘 超，曹博召，王 健，任 艳

（1.西北农林科技大学水土保持研究所，陕西 杨凌 712100；2.西北农林科技大学资源环境学院，陕西 杨凌 712100）

引言

随着我国生产力水平不断发展，人民对美好生活提出了更高要求，生产建设项目在社会发展和进步下也逐渐增多［1-3］。然而，生产建设项目的增多造成地表裸露频现、土方堆置松散、人类机械活动频繁、施工工艺和运行方式多元，这些现象常诱发强烈的水土流失，该类水土流失具有侵蚀类型多样、发生规律复杂、时空分布不均和危害严重等特点［4］。针对生产建设项目水土流失，国内外学者开展了大量研究。Rovira等［5］研究表明，采矿区河流泥沙与生产建设项目工程堆积体水土流失密不可分；钟诚等［6］对广济水电站建设项目的研究表明，项目建设是造成区域水土流失的主要原因，且弃土场水土流失占新增水土流失总量的67.28%；蒲玉宏［7］对废弃堆积物坡面侵蚀特征研究发现，废弃堆积物不仅易造成砂砾化面蚀和泥石流等侵蚀，其堆放形式也易导致崩塌、滑塌等现象的发生；也有研究表明，工程堆积体坡面侵蚀量是裸露荒坡的10.76～12.33倍［8］；Riley［9］针对堆积体开展小尺度冲刷试验表明，弃渣场坡面可蚀性是原始坡面的10 倍～100 倍。可见，与裸露自然坡面相比，工程堆积体侵蚀特点改变，土壤侵蚀模数陡增。郭明明等［10］对神府地区弃土弃渣侵蚀特征的研究表明，由于工程堆积具有随意性、不均匀性，堆积体堆积形式不尽相同，极大程度加大了土壤侵蚀模数的预测难度。王继增等［11］对生产建设项目土壤侵蚀模数的研究表明，土壤侵蚀模数监测是生产建设项目水土保持工作的重点和难点。

为了有效计算土壤侵蚀模数，我国2018 年颁布的《生产建设项目水土保持技术规范》（GB50433-2018）对水土保持工作具有合理的指导性。但因生产建设项目水土保持监测时间序列短，水土流失量预测基础工作开展不及时，使得预测结果常与真实值偏差较大［12-13］。李宏伟等［14］针对黄土区工程堆积体可蚀性开展了较为深入的研究，对工程堆积体土壤侵蚀模数的确定产生推动作用。目前，关于工程堆积体水土流失规律的研究多集中于坡面侵蚀特征方面，针对工程堆积体土壤侵蚀模数预测的研究较少，对于堆积体堆置形式对土壤侵蚀模数的影响研究尚未发现。

在预测土壤侵蚀模数时，常运用水蚀预报模型［15］（Water Erosion Prediction Project，WEPP）。WEPP 模型是以天为步长预测径流量、泥沙量以及量化降雨事件中水土流失情况的土壤侵蚀预报模型，具有精确预估坡面尺度水土流失的优点。任柯蒙等［16］应用WEPP 模型提出了不同降雨强度下最优减流阻沙台面宽度；Rachman等［17］通过WEPP模型对流域降水和土壤沉积量模拟，所得结果与观测值基本一致；基于田间试验数据，熊勤学等［18］应用WEPP模型提出了5°坡度下4种植物篱的最佳种植间距；刘世梁等［19］模拟了不同坡度条件下梯田空间配置的水土保持效应，表明WEPP模型可以很好地定量模拟土壤在坡面的流失情况。然而，应用该模型进行生产建设项目工程堆积体堆置坡度优化的研究相对较少。因此，本研究通过分析生产建设项目的弃土堆置特征，运用WEPP模型预测生产建设项目工程堆积体土壤流失量，探索不同堆置方式下土壤侵蚀规律和不同坡型条件下最优堆置坡度，并以此为基础丰富生产建设项目水土保持理论，进而指导生产建设项目水土保持。

1 试验概况

1.1 研究区概况

为探寻工程堆积体最优堆放形式，试验以位于陕西省汉中市洋县某生产建设项目工程堆积体为例。工程位于陕西省南部，属洋县范围（33°02'N ～33°43'N，107°11'E ～108°33'E）。项目所在区地势东北高陡，南部低缓，中部低平；气候属大陆性季风气候，温和湿润，四季分明。该区多年平均气温14.5 ℃，年最高气温38.7 ℃，最低气温-10.1 ℃，降雨主要集中在7～10 月份，多年平均降雨量达839.7 mm。土壤以黄棕壤为主，土体强度遇水衰减快，抗风化能力差，土壤可蚀性因子参考值为0.003 9（t·hm2·h·hm-2·MJ-1·mm-1）［20］。

1.2 试验设计

项目设计土方量为5×105m3，堆土基底占地面积5×104m2，堆土形式为圆台状。堆土地基水平，以垂直于堆土地基的剖面确定坡型参数，设计坡型有直线型（图1（a），折线型（图1（b）），台阶型（图1（c））3种坡面。在保证气候文件，土壤文件，作物管理文件一致的情况下改变坡度文件，运行WEPP模型，计算得到不同坡型条件下土壤侵蚀模数。利用SPSS23.0 在显著性水平P＜0.05 的条件下多元拟合得到坡型参数与土壤侵蚀模数间最优拟合方程。

图1 工程堆积体坡型Fig.1 Slope type of engineering accumulation

1.3 WEPP模型文件的建立

WEPP 模型坡面版在科研、教学、生产单位的共同攻关下不断修正与完善，目前坡面版的应用较为成熟［21］。WEPP模型将土壤侵蚀因子概化为模型内部可以识别的语言，基于物理过程生成模型所需要的土壤、坡度、气候、作物管理文件［22-23］。

1.3.1 土壤文件的建立

文中以项目区典型土壤黄棕壤建立土壤文件。根据表1 实测土壤理化性质，计算［24］得到试验土壤反照率为0.596，初始饱和导水率为70%，土壤临界剪切力为1.80 Pa，细沟土壤可蚀性0.015 8 s/m，细沟间土壤可蚀性2.81×106kg·s/m4，有效水力传导系数为3.712 mm/h；同时输入砂粒含量（%）、黏粒含量（%）、有机质含量（%）、阳离子代换量（CEC）（meq/100 g）和砾石含量（%）五个分层土壤特性。

表1 试验用土壤物理化学性质Table 1 Physical and chemical properties of soil

1.3.2 坡度文件的建立

坡度坡长数据库主要包括坡长、坡宽、坡度、坡面形状4 类描述数据，利用WEPP 模型中的Slope Profile Editor编辑器编辑4类描述数据，根据所输入的坡型参数，通过单击页面Preview 按钮更新预览随着参数变化的坡面形状。

1.3.3 气候文件的建立

气候文件包括降雨量、降雨历时、最大雨强出现的历时与总历时的比率、最高温度、最低温度、太阳辐射、风向、风速、露点温度9个逐日参数。在预测土壤侵蚀量时，降雨、最高温度和最低温度是主要气象因子，而太阳辐射、风速、风向等气候因子为次要气候因子［25］。通过分析7 个气候参数（降雨量、最高温度、最低温度、太阳辐射量、风速、风向和露点温度）对土壤侵蚀量的敏感性发现降雨量对土壤侵蚀量最敏感，太阳辐射、风速、风向等次要气候因子对土壤侵蚀预测结果影响不大［26-28］。结合中国气象局气象数据中心公布的洋县1981年～2010年30年气候资料，建立洋县1981年～2010年的30个PAR气候文件。

1.3.4 作物管理文件的建立

作物管理文件可结合当地实际的耕作方式、耕作时间、作物种类等来确定，包括与作物生长和耕作等相关的详细资料。试验地初始条件为填土斜坡，坡面不布设排水设施，不采用耕作措施，按照裸露休闲地处理，采用模型自带的休闲地文件建立作物管理文件。

1.4 分析方法

1.4.1 模型验证与评价方法

模型通过相对误差（δ）进行验证：

式中：XP为模型模拟值；X0为实测值，一般认为当δ＜10%时，模型模拟效果较好。

依据某客运专线水土保持监测报告，位于洋县的某弃土弃渣区建设期实测土壤侵蚀模数为17 635 kg·m-2·a-1。将构建的参数文件导入模型，模拟结果换算后得到坡面土壤侵蚀模数为17 870 kg·m-2·a-1，相对误差δ=1.33%（δ＜5%），模拟效果较好，表明WEPP模型在洋县土壤流失量预测上具有适用性。

1.4.2 边坡稳定性分析

边坡的稳定性受坡高、坡角、边坡形式、土壤物理性质等较多因素的影响。本文基于洛巴索夫图解法，按照工程技术规范要求选取安全系数K=1.15［29］，运用极限平衡理论分析不同坡型边坡稳定性，以满足工程堆积体堆置要求的最小底坡坡角（14°）和试验区土壤休止角为最大坡角（圆锥法测得试验区土壤休止角为47°），在此条件下通过SPSS23.0 分析得到最大理论堆土高度与工程堆积体底坡坡角的最优拟合方程（式2）和关系曲线（图2）：

图2 最大堆积高度与底坡坡角关系Fig.2 The relationship between the maximum pile height and the slope angle of the bottom slope

式中：α1表示底坡坡角，（°）；Hm表示理论最大堆土高度，（m）。一般认为，当实际堆土高度H＜Hm时，满足工程稳定性要求。

1.4.3 堆积体尺寸求解圆台型堆积体可以看作是直角梯形绕纵轴一周而成的旋转体［30］，则该圆台体积为：

侧面积则按旋转体侧面积计算［31］：

式中：V表示堆积体体积，（m3，本试验选取V=5×105m3）；r1、r2分别表示堆积体的上下底面半径，（m）；S1、S2是与r1、r2于其相对应堆积体上下底面面积，满足Si= πri2，其中S2= 5× 104m2；h表示堆积体高度，（m）；f(x)是堆积体剖面形状的函数，即f(x)=φ(h)。

1.4.4 最优拟合方程分析

WEPP 模型涉及的参数众多，但在一定土方量和占地面积条件下工程堆积体坡面土壤侵蚀模数仅是底坡坡角α1、上部坡角α2和马道宽d的函数，采用SPSS23.0 软件多元拟合可得到不同坡型下最优拟合方程：

式中：MS表示土壤侵蚀模数，（kg·m-2·a-1）；α1表示底坡坡角，（°）；α2表示上部坡角，（°）；d表示马道宽度，（m）。

对于折线坡和台阶型坡，方程可写成形如Φ(α1)=f(α1,α2,d)的形式，有时，方程在α0处有极小值。

式中：α1ϵ[ 14°,47° ],α2ϵ[ 14°,47° ],d≥0。通过比较Φ( α0),Φ( 14 ),Φ( 47 )求得MS最小值，一般认为，MS取得最小值所对应的α1，α2，d为所需的不同坡型下最优断面尺寸。

2 结果与分析

2.1 直线型坡面对土壤侵蚀的影响

2.1.1 直线型坡面土壤侵蚀特征分析

直线型坡面是工程堆积体堆置过程中一种简单的形式，其土壤侵蚀模数模拟变化规律随坡度变化特征如图3所示。由图3模拟结果可见，随着设计底坡坡角α1的增大，直线型坡面土壤侵蚀模数总体表现出先减小后增大的趋势。当设计底坡坡角α1在14°～26°范围内变化时直线型坡面土壤侵蚀模数随着设计底坡坡角α1的增大而减小；在底坡坡角α1=26°时，土壤侵蚀模数取得最小值，此时最小土壤侵蚀模数MS=28.45 kg·m-2·a-1，工程堆积体表面积为5.2 hm2，坡面年平均土壤侵蚀量达W=1 478.93 t；当设计底坡坡角α1＞26°时，直线型坡面土壤侵蚀模数随着设计底坡坡角α1的增大而增大，在最小设计底坡坡角α1min=14°和最大设计底坡坡角α1max=47°时土壤侵蚀强度剧烈，此时工程堆积体坡面最大土壤侵蚀模数MS=35.6 kg·m-2·a-1，表面积达5.1 hm2，坡面年平均土壤侵蚀量达W=1669.42 t。

图3 直线型坡面土壤侵蚀模数模拟结果Fig.3 Simulation results of soil erosion modulus on straight slope

结果表明，直线型坡面在底坡坡角α1为26°附近土壤侵蚀强度最小，其原因是：坡度、坡长的不同组合对坡面土壤侵蚀的影响分为2 个阶段：图4 可见，直线型底坡坡角从最小坡角14°增加到26°过程中，坡度增加百分比占总坡度变化的36.4%，坡长减少百分比总坡长变化的80.7%，该阶段坡长的变化对土壤侵蚀的影响起主导作用，随着坡长的减小，坡面径流不断分散，侵蚀能力不断降低，水流携带泥沙能力降低，当坡度减小到26°后，坡长变化对土壤侵蚀的影响微弱；当直线型底坡坡角从26°增加到最大坡角过程中，坡度增加百分比占总坡度变化的63.6%，坡长减少百分比总坡长变化的19.3%，该阶段坡度变化对土壤侵蚀的影响占主导作用，随着坡度增大，坡面径流量（入渗量）变化不大，坡面流速增大，土壤侵蚀力加强，坡面土壤侵蚀模数随之增大。

图4 直线型坡面坡长随坡角变化关系Fig.4 The relationship between slope length and slope angle on straight slope

2.1.2 直线型堆土方式优化

基于SPSS23.0对直线型工程堆积体底坡坡角α1和土壤侵蚀模数Ms分析得到最优拟合方程如下（式6）：

式中：MS表示土壤侵蚀模数，（kg·m-2·a-1）；α1表示底坡坡角，（°）。由式（7）和（8）可求得直线型坡面参数最优解。直线型最优拟合方程MS是α1的函数，对α1在区间［14°，47°］内求得MS最小值：

2.2 折线型坡面对土壤侵蚀的影响

2.2.1 折线型坡面土壤侵蚀特征分析

堆积体堆放采用折线型时，坡面形状取决于底坡和上部坡角。为了方便计算，先固定一个坡角，根据设定土方体积，再调整另一坡角大小，采用WEPP模型计算坡面土壤侵蚀模数。图5为上部坡角α2保持不变的条件下，土壤侵蚀模数随底坡坡角变化关系。可以看出，在上部坡角α2保持不变的条件下，工程堆积体土壤侵蚀模数随底坡坡角α1的增大而减小；在保证工程堆积体边坡稳定的前提下，研究取当地土壤休止角（47°）为上限，此时工程堆积体侵蚀程度微弱；在底坡坡角α1=14°时土壤侵蚀强度剧烈，设计条件下最大土壤侵蚀模数MS=36.29 kg·m-2·a-1。

图5 土壤侵蚀模数随底坡坡角变化关系Fig.5 The relationship between soil erosion modulus and bottom angle

在工程堆积体底坡坡角α1保持不变的条件下，土壤侵蚀模数随着上部坡角α2的增大整体呈先减小后增大的趋势，底坡坡角α1在［14°，22°］的范围内土壤侵蚀模数MS呈下降趋势，在α2=22°时设计条件下土壤侵蚀模数最小值MS=26.98 kg·m-2·a-1，此时工程堆积体表面积达5.3 hm2，坡面年平均土壤侵蚀量W=1 418.59 t；当α1＞22°时，工程堆积体土壤侵蚀量明显增大，且存在明显的差异（图6）；即上部坡角较大时所对应的土壤侵蚀模数大于上部边坡较小时所对应的土壤侵蚀模数，且在上部坡角趋于项目区土壤休止角（47°）时，工程堆积体土壤侵蚀程度越发剧烈。

图6 土壤侵蚀模数随上部坡角变化关系Fig.6 The relationship between soil erosion modulus and upper slope angle

2.2.2 折线型堆土方式优化

基于WEPP模型得到折线型坡面土壤侵蚀模数预测结果，通过多元拟合得到α1、α2与折线型坡面土壤侵蚀模数MS预测结果最优拟合方程（式9）：

式中：MS表示土壤侵蚀模数，（kg·m-2·a-1）；α1表示底坡坡角，（°）；α2表示上部坡角，（°）。由式（10）～（13）可求得直线型坡面参数最优解。

折线型最优拟合方程Ms的取值与底坡坡角α1和上部坡角α2有关，Ms的最小值可转换为求关于α2的二次函数最小值，即求：

则函数MS(α2)在处取得最小值

令Φ(α1)=Ms(α2)min

由于α1ϵ[14°,47°]，故在区间[14°,47°]恒小于零，故Φ(α1)在Φ(47)处取得最小值，则Ms(α2)min=Φ(47)，此时α2= 21.8°。故折线型工程堆积体最优坡度组合为α1= 47°,α2=22°。

2.3 台阶型坡面对土壤侵蚀的影响

2.3.1 台阶型坡面土壤侵蚀特征分析

基于前期建立的基础数据库，为了简化计算过程，先固定底坡坡角α1和马道宽度d，得到台阶型工程堆积体土壤侵蚀模数随上部坡角α2的变化规律曲线。图7直观的反映了在工程堆积体底坡坡角α1和马道宽度d保持不变的情况下，土壤侵蚀模数MS随着上部坡角α2的增大总体表现出先减小后增大的趋势，α2在［14°，25°］内工程堆积体土壤侵蚀模数MS在设计条件下有最小值，此时土壤侵蚀模数最小值MS=25.41 kg·m-2·a-1，工程堆积体表面积为5.4 hm2，坡面年平均土壤侵蚀量W=1 373.02 t；当上部坡角α2接近项目区土壤休止角（47°）时，工程堆积体土壤侵蚀模数无明显差异，设计条件下土壤侵蚀模数最大值MS=34.46 kg·m-2·a-1，此时工程堆积体表面积达5.4 hm2，坡面土壤侵蚀量W=1 851.91 t。

图7 台阶型坡面土壤侵蚀模数随上部坡角变化关系Fig.7 Relationship between soil erosion modulus of stepped slopes and upper slope angle

图8 为底坡坡角α1和上部坡角α2保持不变的情况下，土壤侵蚀模数随马道宽度d变化关系图。可以看出，在工程堆积体底坡坡角α1和上部坡角α2保持不变的情况下，工程堆积体土壤侵蚀量随着马道宽d的增大而增加，且增长速率趋于稳定；在满足工程实际堆土要求和堆积体边坡稳定的前提下，当马道宽d保持不变时，底坡坡角α1越小，上部边坡坡角α2越大，工程堆积体土壤侵蚀模数MS越大，土壤侵蚀程度越剧烈。

图8 台阶型坡面土壤侵蚀模数随马道变化关系Fig.8 Relationship between the change of soil erosion modulus of stepped slopes and the step

2.3.2 台阶型堆土方式优化

马道作为台阶型堆土方式的重要组成组分，其与底坡坡角α1和上部坡角α2共同影响坡面土壤侵蚀模数。运用SPSS23.0多元拟合得到台阶型坡面最优拟合方程（式14）：

式中：MS表示土壤侵蚀模数，（kg·m-2·a-1）；α1表示底坡坡角，（°）；α2表示上部坡角，（°）；d表示马道宽度，（m）。由式（14）和（15）可求得台阶型坡面坡型参数最优解。

台阶型最优拟合方程可看作关于上部坡角α2的二次函数，即：

则函数Ms(α2)在α2= 25.4处取得最小值。

由于α1ϵ［14°，47°］，d≥0，故α1=47°且d=0 m 时，Ms(α2)取得最小值。综上，台阶型工程堆积体最优边坡组合为α1=47°，α2=24.5°，d=0 m。在堆土高度较大时，考虑到工程实际情况，工程堆积体中部往往需要设计马道，此时为了减少坡面水土流失，马道不宜太宽，且应布设排水、苫盖等水土保持措施。

3 讨论

生产建设造成的大量松散工程堆积体由于其高、陡等特点影响着堆积体的侵蚀特性，堆积体堆置坡度、坡长、坡度坡长组合等不同，不仅影响着堆积体的不稳定性和径流的重力顺坡分力，同时改变了堆积体侵蚀过程［32-33］。本试验结果表明直线型、折线形、台阶型堆积体中底坡坡角α1对坡面土壤侵蚀模数的影响有所差异，直线型坡面底坡坡角α1趋于土壤休止角47°时易发生土壤侵蚀，折线形和台阶型坡面在底坡坡角α1在47°时土壤侵蚀相对最小，主要是底部坡度最大，可有效增加土方堆置，减少上部堆置土方量。陈俊杰等［34］针对不同雨强与坡长条件下坡度对细沟侵蚀的影响认为，降雨条件相同时，由于坡长变短，水流对坡面的冲刷时间降低，水深逐渐减小，侵蚀强度相应减弱，与本文直线型堆积体坡面土壤侵蚀特征结论相近。陈利顶等［35］研究黄土丘陵区坡面形态的土壤侵蚀效应表明，上部为凹型坡面，下部为线型坡面，凹型坡土壤侵蚀的截留功能难以有效发挥，土壤侵蚀程度剧烈，使得该类坡型在整体上表现出较强“源”的作用，其研究结果与文中折线型（α2＞α1）和台阶型堆积体坡面土壤侵蚀特征相一致。

通过WEPP 模型确定不同堆积形式下坡面土壤侵蚀量，采用SPSS23.0 得到不同坡型条件下最优拟合方程，台阶型堆积体在马道宽度d=0 m时就转变成折线形堆积体，但最优拟合方程所求解的最优上部坡角α2存在差异，这是由于在建立拟合方程时折线型最优拟合方程使用了50 组数据点，而台阶型最优拟合方程使用了100组数据点，但二者所得到的最优拟合方程相关系数R2均趋近1（R32=0.956，R42=0.973），模拟得到令人满意的结果［36］。WEPP模型以建立的文件为基础，模型所需要的土壤、气候文件等直接影响着预测结果的可靠性。研究表明，坡面侵蚀产沙受土壤类型的影响较大［37］，魏晓燕等［38］针对不同土壤坡面产流产沙特征对比分析认为，土壤质地的不同，直接影响区域产沙能力；李建明等［39］针对不同土质工程堆积体径流产沙差异研究认为，砂土堆积体平均侵蚀速率分别是壤土和黏土的3.0 和2.3 倍。由此可见，不同土质土壤在侵蚀产沙能力上差异性显著，而笔者主要针对项目区典型土壤黄棕壤开展研究，因此其它类型土壤侵蚀特征有待进一步分析。本文主要针对堆积体土壤侵蚀特点进行初步探讨，将工程堆积体概化为简易圆台型，计算一定土方量和占地面积相对固定下的土壤侵蚀量，针对不同下垫面处理、土方量、占地面积条件下堆积体堆积形式的优化有待深入研究。

4 结论

利用气象站1981～2010 年30 年间的气象观测资料以及项目区实测土壤性质，建立气候文件和土壤文件，基于WEPP模型预测生产建设项目建设过程中不同堆置方式下工程堆积体土壤侵蚀模数，在一定土方量和占地面积条件下寻求优化不同坡型堆置坡度最优组合方式，得到以下结论：

（1）直线型堆土方式下工程堆积体土壤侵蚀模数与底坡坡角α1有关，设计条件下直线型工程堆积体最优底坡坡角α1=26°，此时坡面土壤侵蚀程度最弱；当工程堆积体α1接近14°或者趋于当地土壤休止角47°时，土壤侵蚀剧烈，易诱发严重的水土流失。

（2）折线型堆土方式下工程堆积体底坡坡角α1和上部坡角α2的不同组合直接影响工程堆积体土壤侵蚀模数，底坡坡角α1趋于项目区土壤休止角47°，上部坡角α2设计为22°时，折线型工程堆积体土壤侵蚀程度较弱，设计条件下折线型工程堆积体最优边坡组合为α1=47°，α2=22°。

（3）台阶型堆土方式下工程堆积体底坡坡角α1、上部坡角α2和马道宽度d的不同组合共同作用工程堆积体土壤侵蚀模数。随马道宽度增大，土壤侵蚀模数不断增大，在满足设计要求的前提下台阶型工程堆积体断面最优组合为α1=47°，α2=24.5°，d=0 m。