基于137Cs示踪和WaTEM／SEDEM模型的黑土区小流域侵蚀产沙模拟

李国强，方海燕

（1杭州市水文水资源监测总站，浙江 杭州310016；2中国科学院地理科学与资源研究所 陆地水循环及地表过程重点实验室，北京100101）

土壤侵蚀是世界性的环境问题，它不但引起土壤及其土壤养分的流失，侵蚀掉的泥沙进入河道后还会引起湖泊、水库的淤积以及环境污染等问题［1］。东北黑土区是我国重要的商品粮基地，黑土有机质含量丰富、肥力高、性状好，土地农作物产量高。自20世纪50年代以来，随着人口的增加和开垦强度的增大，黑土层厚度从开垦初期的60～70cm下降到目前的20～30cm，有些地方已出现破皮黄［2］，严重威胁到国家的粮食安全。

土壤侵蚀模型是土壤侵蚀和水土保持研究的重要手段。然而，由于东北黑土区土壤侵蚀研究相对较晚，小流域径流、产沙观测相对较少，该区土壤侵蚀模拟研究多采用 （R）USLE（（Revised）Universal Soil Loss Equation）模型［3－5］，但（R）USLE 模型计算侵蚀量时忽略了泥沙的沉积和空间的变异性［6］。其他的一些物理机制模型如 WEPP和LISEM等对数据要求十分严格，在黑土区流域尺度上尚无应用。比利时鲁汶大学研发的 WaTEM／SEDEM（Water and Tillage Erosion Model and Sediment Delivery Model）模型所需数据少，在土壤、地形、土地利用和降雨等基础数据具备的情况下，即可运行［7］。此外，该模型结构简单，所需校正参数少，同时又能模拟土壤侵蚀产沙及分布特征；除了产沙观测数据外，湖泊、水库等的泥沙淤积信息也可对该模型进行校正和验证。因此，自 WaTEM／SEDEM模型提出以来，在国外许多国家［8－10］得到了很好的应用，在我国的红壤区小流域［11－12］也利用该模型开展了工作。

截至目前，国内外多应用出口产沙校正和验证WaTEM／SEDEM模型，使得不同的模型参数组合都能很好地模拟流域出口产沙量，表明流域内侵蚀沉积的多种组合均会有相同的流域出口产沙［13］。然而，利用流域内部侵蚀沉积强度数据校正模型，将可能真实地再现流域侵蚀沉积特征，同时也能够模拟流域出口产沙量。

东北黑土区土壤侵蚀研究开展较晚，一些地区缺乏水沙观测数据，但该区有许多的湖泊和水库泥沙淤积信息数据，从而为该模型的应用提供了很好数据支持。此外，核示踪流域内侵蚀强度在黑土区已成功开展，利用核示踪得到的土壤侵蚀强度可校正模型［13］。因而，在黑土区联合应用 WaTEM／SEDEM模型及核示踪技术开展土壤侵蚀研究，将是弥补该区侵蚀产沙数据较少、模型模拟研究不足的重要举措。

1 材料与方法

1.1 研究区概况

拜泉县（E125°30′～126°31′，N47°20′～47°55′）位于黑龙江省齐齐哈尔市东北部，地处小兴安岭余脉与松嫩平原的过渡地带，属于东北典型黑土区，地形特征为漫川漫岗。拜泉县为中温带大陆性气候，多年平均气温为1.28℃，全年平均降水量500mm左右，多集中在7、8、9三个月，降雨历时短，强度大，易形成径流侵蚀。土壤以典型黑土和草甸黑土为主。流经县境内的河流有通肯河、双阳河、润津河。为了发展农业，拜泉县从20世纪60年代开始修建水库塘坝等蓄水工程，沉积了大量泥沙。本文选取该县境内的24个水库，水库控制流域土地利用较为单一，以耕地为主（表1）。

1.2 数据来源

1.2.1137Cs数据采集与处理 开展研究的24个水库中，齐心水库控制流域地形适中，土地利用类型较为全面，具有一定的代表性。因此，2013年10月，在齐心水库控制流域，根据土地利用类型、地形特点，利用5cm直径的土钻采集土壤样品，采样深度视黑土层厚度而定，坡面侵蚀区采集深度40cm左右，沉积区100cm左右，以保证采样点所有的137Cs含量都能够采集到。采集的25个样点土壤经晾干、研磨后过2mm筛，剔除草根和石块等杂物，称300g左右供核元素测试。137Cs放射性活度用美国CANBERRA公司配备的高纯锗探头λ能谱仪（测定范围为3～3×103keV）测定。在662keV处求算137Cs放射性活度（Bq／kg），核素测定时间为 80 000s左右，测定误差 ＜8%。

表1 研究区24座水库的控制流域面积、耕地比例和年均产沙量Tab.1 The areas，percentages of cultivated lands and sediment yields of the 24reservoir catchments

其中：A（t）为单位面积累积137Cs活度（Bq／m2，由采样点核素放射性活度与土钻横截面积得到）；R为侵蚀强度（kg／（m2·a））；t为137Cs释放以来经历的时间；d为累积质量深度（kg／m2），λ为137Cs衰变常数；I（t）为137Cs年沉降通量（Bq／（m2·a））；Γ为土壤侵蚀新沉降的137Cs在耕作之前被侵蚀移动的比例，P为颗粒校正系数。

在沉积样点（土壤样品核素活度大于背景值区核素活度）：

1.2.2 土壤侵蚀强度计算 MBM2模型在黑土区核素含量和土壤侵蚀转换中已得到成功运用［14］。在侵蚀样点（土壤样品核素活度小于背景值区核素活度）：

其中：Aex为测量得到的137Cs活度与其背景值Aref之差；Cd（t′）为沉积采样点土壤的137Cs活度（Bq／kg）。

1.3 WaTEM／SEDEM 模型的运行

1.3.1 模型简介 WaTEM／SEDEM是在RUSLE基础上构建的分布式模型，它分为侵蚀模块、泥沙输移能力模块和泥沙运移模块［15］。模型年均输沙能力Ct（kg／（m·a））用来反映地表植被覆盖对输沙能力的影响，模型仅需要校正泥沙输移能力系数Ktc。地表植被覆盖较好，Ktc值较低；地表植被覆盖较差，则Ktc较高。

1.3.2 模型数据来源及处理

（1）数字高程模型（DEM，分辨率为20m）由国家基础地理信息中心1991年航摄的1∶5万地形图数字化得到；（2）土地利用数据利用2000年TM影像（30m分辨率）采用人工目视解译的方式获取，并进行了野外考察校正；（3）降水数据由拜泉县气象局及中国气象科学数据共享网获得；（4）拜泉县土壤类型空间分布图（1∶10万）由中国科学院东北地理与农业生态研究所获得；（5）水库控制流域年均产沙量（SY）由黑龙江水文图集中44个观测站36年实测产沙利用水文比拟法［16］获取，并经三道镇水库和齐心水库高精度GPS差分法和实地采样得到了验证。

WaTEM／SEDEM模型输入数据有DEM、降雨侵蚀力因子R、土壤可蚀性因子K、土地管理因子C和水土保持措施因子P，河流水系图等。河流水系由DEM通过ArcGIS水文模块提取，汇水面积大于1 km2即可视为河道［9］；研究区C和P值参考黑土区已有研究［5，17］及实地调查得到。土壤可蚀性因子K根据拜泉县1984年土壤普查资料及黑土区有关研究成果［17－18］。研究区及周边地区有些站点缺少逐日降水数据，利用年均降雨量方法获取降雨侵蚀力R 值［19］：

其中：R为多年平均降雨侵蚀力（MJ·mm／（h·hm2·a）），P为多年平均降雨量（mm）。在ArcGIS9.3中通过克里金插值方法得到各水库控制流域降雨侵蚀力值。

1.3.3 模型校正 利用137Cs示踪法，根据25个样点的土壤侵蚀强度对模型参数进行校正。先初步根据已发表的论文确定Ktc的范围，以减少模型运行次数，之后根据Ktc组合，得到流域侵蚀沉积特征，提取各采样点的侵蚀沉积强度，并与137Cs示踪得到的侵蚀沉积速率对比；根据Ktc高值和低值的多次组合运行模型。采用模型有效系数（cNS）确定Ktc最佳参数组合。cNS计算方法如下：

其中n是指观测的数目；Oi是指观测值；Omean是指观测平均值；Pi是指预测值。cNS的取值范围是（－∞，1），cNS值越接近于1，模拟效果越好。

2 结果分析

2.1 采样点侵蚀沉积强度

有研究［20］已经测量了研究区的137Cs背景值，因而，根据背景值和25个采样点核素面积活度及MBM2（公式1和2），得到了采样点侵蚀强度。侵蚀强度受土地利用类型和地形的影响，坡耕地侵蚀强度平均值为－24.5t／（hm2·a），最大侵蚀强度达到－45.7 t／（hm2·a），最大沉积速率为15.4t／（hm2·a）。林地平均侵蚀强度为－1.1t／（hm2·a），最大侵蚀强度为－4.5t／（hm2·a），最小值（沉积速率）为0.6t／（hm2·a）。灌木林地侵蚀强度介于二者之间，平均侵蚀强度为－5.6t／（hm2·a）。研究区草地平均侵蚀强度为－9.5t／（hm2·a）。整个齐心流域采样点侵蚀速率变异很大，平均侵蚀强度为－19.8t／（hm2·a）。

2.2 模型校正和验证

利用以上25个采样点的侵蚀沉积强度数据对模型进行校正。根据地表径流泥沙输移敏感性的不同，对两种截然不同的土地利用类型赋予不同的Ktc值。林草地泥沙输移能力系数赋予低值Ktc－low，耕地则赋予高值 Ktc－high，Ktc变化范围在0～12之间。Ktc－low和Ktc－high不同组合下得到采样点模拟侵蚀强度，通过137Cs示踪和模拟侵蚀强度公式（4）计算 NS变化值。图1显示最佳 Ktc－low和 Ktc－high值组合为4m和8m，此时cNS为0.62。

根据137Cs示踪法校正后的模型，利用表1中24个水库（不包括齐心水库）控制流域产沙量（YS）与产沙模数YSS（YS与流域面积的比值）进行模型验证。图2表明，模型模拟的YSS和YS值和观测值相近，其相关系数分别达到0.86和0.7，模型对于年产沙的模拟效果优于产沙模数。

表2 核素137Cs示踪得到的采样点土壤侵蚀（负数）沉积（正数）强度Tab.2 Activities of 137Cs and its derived soil erosion rates for the sampling points

3 讨论

在东北黑土区，已成功利用核示踪技术开展了土壤侵蚀研究。方海燕等［14］研究发现，在东北鹤山农场面积为28.5hm2的农耕地小流域上，土壤净侵蚀强度为2.0t／（hm2·a），远小于本研究中耕地平均土壤侵蚀强度；而流域内侵蚀区的侵蚀强度均值为－15.0t／（hm2·a），略小于本研究得到的土壤侵蚀强度。在拜泉县小流域坡面上，于寒青等［20］发现等高种植、梯田和顺坡坡面上土壤平均侵蚀强度为分别为－23、－25和－35t／（hm2·a），大于本研究平均值，这是因为他们的采样点位于坡度较大的坡面侵蚀区。从以上比较来看，研究中得到的侵蚀速 率数据可信。

图1 利用采样点的侵蚀速率对泥沙输移能力系数Ktc－low（a）和Ktc－high（b）进行校正Fig.1 Validations of the Ktc－low（a）and Ktc－high（b）values using 137Cs－derived erosion rates of the sampling points

图2 水库控制流域产沙观测值与WaTEM／SEDEM模拟值：（a）年均产沙量，（b）年均产沙模数Fig.2 Modeled sediment yield per year（a）and specific sediment yield（b）versus the counterparts from observed values for the 24reservoir catchments

研究发现，利用流域出口产沙数据校正模型时，虽然有时流域出口产沙的模拟效果很好，但流域内土壤侵蚀沉积模式可能有多个。如在西班牙Barasona水库流域，Alatorre等［13］发现有一系列的Ktc高值和低值组合，均能得到相同的cNS值。137Cs之所以能够校正WaTEM／SEDEM模型，除了能够反映空间上土壤侵蚀沉积模式外，另外一个重要特点即是二者的时间尺度相同或相似。137Cs示踪得到的是50年以来土壤的平均侵蚀强度；WaTEM／SEDEM也是一个模拟多年时间尺度上土壤平均侵蚀强度的模型，只要模型模拟时将二者的时间跨度调为一致，即可消除时间尺度的影响。目前，利用137Cs示踪得到的土壤侵蚀强度校正 WaTEM／SEDEM模型已成功应用［13］。

本研究中，利用核示踪法得到的最佳Ktc组合为4m和8m，该值小于Shi等人［11］在王家桥流域Ktc组合值（15m 和55m），也小于 Van等人［8］利用西班牙21个流域产沙值得到的Ktc值组合（75m和250m）。在西班牙小流域，Alatorre等人［13］利用137Cs方法得到的Ktc最佳组合值为2.84m和9.05 m，与本研究较为接近。在数据不充分时，有些研究甚至将 Ktc－low和 Ktc－high固定为一个比值［11］。有研究发现［21］，不同侵蚀类型区和不同数字高程模型分辨率情况下，校正得到的参数 Ktc－low和 Ktc－high值差别会很大。齐心水库控制面积为9.9km2，研究中25个土壤采样点偏少，但利用采样点获得的侵蚀沉积速率对模型校正后，流域出口产沙的模拟效果很好（图2）。值得注意的是，YS较YSS的模拟精度高，这是因为YS值变化幅度大，因而对于YS的模拟效果优于YSS模拟能力。

为进一步探讨模型模拟的精确性，利用校正好的模型得到了各个小流域总侵蚀量，进而根据小流域总产沙量与总侵蚀量的比值，得到了泥沙输移比RSD（图3），结果表明研究区小流域泥沙输移比变化于0.16和0.40之间，且随着流域面积的增加，泥沙输移比呈现减小的趋势。这符合传统的产沙模数－流域面积关系［22］。此外，在鹤山农场，Dong等［23］研究发现面积为27.29km2的小流域泥沙输移比为0.1；简金世［24］在东北漫川漫岗区三个水文站得到平均泥沙输移比为0.145；Fang等［14］发现黑土区28.5hm2的小流域RSD大于0.5，这可能与流域面积小有关。此外，野外考察发现，河道上游的沟道表现为侵蚀，而河道中下游多有泥沙沉积。这一方面说明图3的合理性，另一方面也表明与面积有关的YSS模拟效果要低于YS，这也从侧面证实了图2的准确性。

以上表明，利用137Cs示踪得到的土壤侵蚀沉积速率能够很好校正WaTEM／SEDEM模型。因此，在水沙监测资料稀少地区，核素示踪和模型 Wa－TEM／SEDEM的结合运用，对于该区及类似研究区的土壤侵蚀研究具有很好的应用前景。

图3 水库控制小流域泥沙输移比随流域面积的变化特征Fig.3 Changing pattern of sediment delivery ratio with catchment areas

4 结论

WaTEM／SEDEM是一个基于RUSLE的模型，模型参数少，数据需求量低，水文站实际观测数据和淤积信息及137Cs法都可对模型进行校正和验证。利用24个水库流域产沙量验证后，发现模型模拟的效果较好。取得的结论如下：

（1）采样点坡耕地侵蚀强度为－45.7～－15.4 t／（hm2·a），平均为－24.5t／（hm2·a），林地侵蚀强度较小，仅为－1.1t／（hm2·a）。灌木林地侵蚀强度介于二者之间，平均侵蚀强度为－5.6t／（hm2·a）。整个齐心流域采样点侵蚀速率变异很大，平均侵蚀强度为－19.8t／（hm2·a）。

（2）研究区模型泥沙输移系数Ktc最佳组合为4 m和8m，此时模型有效系数cNS为0.62。通过24个水库控制流域验证得到产沙量YS和产沙强度YSS的模拟值与实测结果相近，模拟效果很好。

（3）核示踪技术和 WaTEM／SEDEM模型联合运用，能够很好地开展黑土区土壤侵蚀和水土保持研究，也可为该区的水土保持工作提供科学依据和理论指导。

［1］Novotny V，Chesters G.Delivery of sediments and pollutants from nonpoint sources：A water quality perspective［J］.Journal of Soil Water Conservation，1989，44：568－576.

［2］范昊明，蔡强国，王红闪.中国东北黑土区土壤侵蚀环境［J］.水土保持学报，2004，18（2）：66－70.

［3］张宪奎，许靖华，卢秀琴，等.黑龙江省土壤流失方程的研究［J］.水土保持通报，1992，12（4）：1－9.

［4］张树文，王文娟，李颖，等.近50年来三江平原土壤侵蚀动态分析［J］.资源科学，2008，30（6）：843－849.

［5］王文娟，张树文，方海燕.东北典型黑土区坡沟侵蚀耦合关系［J］.自然资源学报，2012，27（12）：2113－2122.

［6］Verstraeten G.Regional scale modeling of hillslope sediment delivery with SRTM elevation data［J］.Geomorphology，2006，81（1／2）：128－140.

［7］WaTEM／SEDEM manual ［EB／OL］.http：／／www.kuleuven.be／geography／frg／modelling／，2014－11－30.

［8］Van R A，Paolo B，Robert J，et al.Modeling sediment yields in Italian catchments［J］.Geomorphology，2005，65（1／2）：157－169.

［9］Verstraeten G，Prosser lan P，Fogarty P.Predicting the spatial patterns of hillslope sediment delivery to river channels in the Murrumbidgee catchment，Australia［J］.Journal of Hydrology，2007，334（3／4）：440－454.

［10］Boix－Fayos C，De V J，María M M，et al.The impact of land use change and check－dams on catchment sediment yield［J］.Hydrological Process，2008，22（25）：4922－4935.

［11］Shi Z H，Ai L，Fang N F，et al.Modeling the impacts of integrated small watershed management on soil erosion and sediment delivery：A case study in the Three Gorges Area，China［J］.Journal of Hydrology，2012，438／439：156－167.

［12］方海燕，孙莉英，聂彬彬，等.基于 WaTEM／SEDEM模型的双枫潭流域侵蚀产沙模拟 ［J］.陕西师范大学学报：自然科学版，2014，42（1）：92－97.

［13］Alatorre L C，Begueria S，Lana－Renault N，et al.Soil erosion and sediment delivery in a mountain catchment under scenarios of land use change using a spatially distributed numerical model［J］.Hydrology and Earth System Sciences，2012，16（5）：1321－1334.

［14］Fang Haiyan，Sun Liying，Qi Deli，et al.Using137Cs technique to quantify soil erosion and deposition rates in an agricultural catchment in the black soil region，Northeast China［J］.Geomorphology，2012，169／170：142－150.

［15］盛美玲，方海燕.WaTEM／SEDEM模型及其应用研究进展与展望［J］.地理科学进展，2014，33（1）：85－91.

［16］孙和平.小型水库泥沙淤积成因分析及淤积量的初步估算［J］.地下水，2005，27（3）：221－222.

［17］魏建兵，肖笃宁，李秀珍，等.东北黑土区小流域农业景观结果与土壤侵蚀的关系［J］.生态学报，2006，26（8）：2608－2615.

［18］翟伟峰，许林书.东北典型黑土区土壤可蚀性K值研究［J］.土壤通报，2011，45（5）：1209－1213.

［19］章文波，付金生.不同类型雨量资料估算降雨侵蚀力［J］.资源科学，2003，25（1）：35－41.

［20］于寒青，李勇，Nguyen M L，等.基于FRN技术的我国不同地区典型土壤保持措施的有效性评价［J］.核农学报，2012，2（2）：340－347.

［21］Verstraeten G，Poesen J，Gillijins K，et al.The use of riparian vegetated filter strips to reduce river sediment loads：An overestimated control measure？［J］Hydrological Processes，2006，20：4295－4267.

［22］方海燕，蔡强国，李秋艳.产沙模数与流域面积关系研究进展［J］.地理科学进展，2008，27（6）：63－69.

［23］Dong Y F，Wu Y Q，Zhang T Y，et al.The sediment delivery ratio in a small catchment in the black soil region of Northeast China［J］.International Journal of Sediment Research，2013，28（1）：111－117.

［24］简金世.松花江流域不同侵蚀类型区泥沙输移比的估算［D］.西安：西北农林科技大学资源环境学院，2011.