随着城市、人口和产业增长,城市面源对地表水环境的污染不断加剧^[1,2] ,在全球范围已成为水质问题的首要因素^[3,4,5] 。其首要污染物氮、磷可导致淤积、富营养化、藻类增殖、水体缺氧等^[6,7] ,威胁水资源和水环境容量^[8,9] 。水土流失(土地元素的淋蚀、输运、沉积等地表过程)通常是面源污染的重要支配条件^[6,9] 。在中国,水土流失面积占到了全国面积37%^[10] ,受人工干扰的山地丘陵地区尤为严重^[11] ;山坡易蚀的农用地常导致更多面源污染^[12,13] 。在中国南方山地和丘陵地区,水土流失现象正逐年加重^[6] 。在独特地貌作用下,丘陵盆地对水土流失及污染的更为敏感。

　　 水土流失等面源污染难以采用常规的点源防治措施^[2,3,14] ,对土地利用和土地覆盖(LULC)应采取综合性、战略性的环境对策。作为减轻开发建设对原始水文干扰的雨洪管理理念,海绵城市提倡自然渗透、自然净化的低影响开发(LID)途径,目前在中国正得到大力推广^[15] 。为了增加决策的可靠性,LULC改变后的水土流失对揭示规划效果和趋势至关重要^[5] 。海绵城市实践(或LID应用)已被很多研究提及,然而现阶段依然缺乏对土地利用和土地覆盖变化(LUCC)的情景分析及对措施效果和功能的预判^[15] 。

　　 诸多水土流失模型被开发以模拟面源污染的关键行为^[5] 。通用水土流失方程(USLE)^[16] 自1960年代起被广泛应用^[17] ;修正通用水土流失方程(RUSLE)^[18] 被用于估计长期事件的水土流失;为了对土壤产生、运移和沉积进行空间显式模拟,揭示水土流失的格局、速率和水体接纳总量,一些GIS技术和空间数据已整合到RUSLE方程,例如CRE-AMS,SWAT,ANSWERS,SWRRB,AGNPS,EP-IC,WEPP,EUROSEM,INCA,LASCAM,HSPF,WATEM-SEDEM等^[6,17,19,20] ,部分模型提供了基于LULC、地貌、土壤可蚀度的集水区污染物长期格局和负荷模拟^{[2,3,7,13,19,21]} 。另外,将简单过程(降雨径流、土壤侵蚀、产沙运移等)整合以模拟复杂过程也是常用方法^[22] ,其代表如IDRISI(Clark大学,1987)等。

　　 本研究将基于IDRISI对长期高频事件导致的丘陵盆地降雨、径流、侵蚀和沉积进行概化,得出水土流失格局。先对基地现状进行模拟,为特定类型和规模的海绵设施布局提供引导;然后对规划情景进行模拟,通过结果对比评估规划方案对现状的改进程度。

　　 黄石市位于长江中游南岸、湖北东南部(30°08′~30°16′N,114°56′~115°09′E),属亚热带季风气候,四季分明,雨量充沛。年均气温17℃,最热月和最冷月均温分别为29.2℃和4.3℃,年降雨量1 400mm,年降雨日132天,降水时间分布不均及全年湿润的气象条件为水力侵蚀提供重要驱动力。磁湖坐落于市区中部,三面群山环抱,环湖地区属于典型的丘陵盆地。上述气候及地理条件皆为基地的水土流失污染防治提出了迫切须求。

　　 磁湖水域面积约10km²,是黄石市最重要的城中湖,发挥着园林景观、洪水调蓄、生态调节等功能。由于平均水深1.75m,磁湖的水环境相当脆弱。随着城市发展和人口集聚,大量污水和径流污染导致了湖泊水质的恶化。在城市点源实现全面截污后,地表径流引起的面源污染将成为磁湖主要污染源。为此,在流域范围开展综合环境治理,引入海绵设施的保持水土、截污自净功能,对环湖区域可持续发展有重要意义。

　　 本研究将包括:①在田间尺度对长期高频水土流失进行空间显式模拟,计算基地范围平均土壤侵蚀模数[MSEM,t/(hm²·y)]和沉积物入湖总量(NASLW,t/y),并识别现状水土流失关键地区;②诊断现状土壤侵蚀重点地区,针对性布局海绵城市措施优化现有土地格局;③对规划情景采用类似模拟得出规划后的MSEM和NASLW,并通过水土流失变化面积揭示方案对水土流失污染防治的影响。

　　 将主要采用IDRISI的两大模块:预测非渠化系统土壤流失的RUSLE模块和评估土壤流失运移的SEDIMENTATION模块。IDRISI的分布式模拟通过将计算区域抽象为许多微小HRU集合^[3,22,23] ,对每个集合采用RUSLE同质化地进行土壤流失量[t/(hm²·y)]计算:

　　 式中R———量化降雨强度的土壤侵蚀效应,MJ·mm/(hm²·h·y);

　　 L和S———分别为坡度和坡长因子,共同决定径流和沉积的运动行为;

　　 K———量化降雨对土壤的侵蚀度,t·h/(MJ·mm);

　　 C———表征不同LULC对水力侵蚀的敏感性;

　　 P———人工干预和管理措施下相对原始情况的侵蚀度^{[5,6,13,18,21]},将进一步用于海绵设施的建模中。

　　 泥沙输移比(SDR)将随后计量侵蚀土壤输运至水系部分的比例,评估净土壤流失的格局^[20,23] ;沉积物在HRU间的运动行为将通过比较运移路径前后相对高差和土壤流失量,通过水流算法确定^[23] 。

　　 由于模型数据的功能尺度决定集合流域计算结果^[3] ,本研究将水土流失过程尺度设为30m^[20,12]。另外,由于水文响应单元(HRU)影响水体溶质通量及其时空变化^[24] ,只有磁湖流域内的因子才会实质影响水环境,本研究首先基于排水方向划分磁湖集水盆地(面积60.0km²);考虑到研究尺度下规划不会改变水系和地形,该汇水区将作为所有研究的计算域。

　　 上述演绎研究将首先基于基地的现状条件。通过考察当前条件下的土壤流失量和沉积量的非空间总量数据,揭示现状急需弥补或改善的水土流失环节或地区,为海绵设施的规划布局提供指导。设施的种类和规模将遵循一般规律和导则,规划措施只会改变现状LULC的细节(导致P因子输入值变化)。由于规划后水土流失因素发生改变,为验证规划是否能缓解现状问题,规划情景将被以相同的方式模拟和评估。

　　 RUSLE和SEDIMENTION的数据准备如下:

　　 R将参考相关案例和RUSLE手册由年降雨量和年侵蚀量的关系得出^[18,20,23] ;参考当地气象局降雨资料,R将统一设为778 MJ·mm/(hm²·h·y)。

　　 L和S将基于RUSLE模块的内部函数(对中国相似地形的案例适用)确定^[13,23] 。L拟计算如下:

　　 式中l———坡长,m;

　　 β、m———通过坡度θ计算^[18] 。

　　 S拟计算如下:

　　 式中α———坡度,%。

　　 L和S的数据源为30m粒度的数字高程模型(DEM),由中国科学院计算机网络信息中心提供。本研究忽略规划对大尺度地形的改变,L和S在规划前后将保持一致。

　　 K由相关RUSLE模型的公式得到:

　　 式中,Sa、Si、Cl和C分别为土壤的沙粒量(%)、粉粒量(%)、粘粒量(%)和有机碳量(%)^[13] 。

　　 K的数据源为30”精度的HWSD栅格数据库,由FAO、IIASA、ISRIC、ISSCAS、JRC等机构提供。研究将忽略规划对土壤属性的扰动;K分布图由Krigin方法差值到30m精度。

　　 C/P由规划前、后的LULC特征得到^{[13,19,20,23]} ,原则上包括现状和规划对应的两种格局;赋值将基于RUSLE手册的取值范围^[18] 。P将为海绵设施功效提供定量方法。

　　 作为影响土壤净流失的因素,面积较大、坡度较缓的HRU会截留更多土壤及面源污染^[20,22] 。对未与水体相连的HRU,SDR一般通过面积估算:

　　 式中A———集水区面积,km²;

　　 α、β———经验参数^[13,19,21] 。

　　 对与水体相连的HRU,SDR为1.00^[23] 。无论规划前后,SDR统一由上述原则赋值。坡度等地形因素由DEM确定^[23] 。

　　 模拟结果包括土壤流失(ASL,t/ha/yr)和沉积去除后的净流失(NASL,t/ha/yr)格局。ASL越高表示源头侵蚀强度越大;NASL为正值(负值)处表示流失过程的“源”(“汇”),正值越大代表土壤沉积的输出流量越大(见图1)。结果显示了磁湖的东南侧裸地、南侧道路沿线和西北侧道路沿线的ASL和NASL较高,这将作为海绵设施规划的现状依据。

　　 为探究土壤流失总量和入湖总量,现对ASL和NASL进行统计。MSEM即计算域内HRU的土壤流失平均值,描述盆地范围内年平均土壤流失量;NASLW由磁湖盆地系统内所有HRU的土壤净流失加和得到,反映直接沉积于湖体的年土壤总量。根据统计,现状MSEM和NASLW分别为315.29t/(hm²·y)和12543.52t/y。这些数据将作为评价规划效果的基准值。

　　 海绵设施拟布局在现状土壤流失较高处,首要目标是缓解现状MSEM和NASLW。参考常用的水土保持手段^[25] ,本研究选取了被广泛应用且有充分指标参数的LID设施,包括滞留池、湿地盆地、湿地调蓄渠、生物滞留系统、植草沟盆地、植被渗滤带、透水铺装^[15,26] ,作为这些设施对水土流失污染物(有实际意义的土壤流失)削减率的归纳,将引入D指数修正P的输入参数(见表1)。

　　 根据RUSLE定义,修正后的P(P’,水土保持措施下的土壤流失与原始状态下的比值^[13,18] )可量化海绵设施对水土流失污染物削减效果:

　　 式中S_i———设施i的面积比例;

　　 D_i———设施i的校正因子(在设施布局区域∈(0,1.00],非布局区域=0),设施交叠区域将采用累乘算法。

　　 一般地,S_i通过地方规范、侵蚀压力、布局条件、适建性等综合确定,D_i将基于设施对各类污染物的消减量,同时结合不同污染物本底含量和实际影响的相对重要性确定。P′的值将被控制在RUSLE手册规定的取值范围内^[18] 。

　　 此外,海绵设施将以合理规模(和现状土壤流失量正相关)被布局在合理区位(接近土壤流失源)(见图2)。海绵设施规划的情景将以相同技术路线被模拟,得出MSEM和NASLW等指标(见图3)。

　　 根据统计,规划后MSEM和NASLW分别为281.38t/ha/yr和12138.47t/yr,分别较现状下降了10.76%和3.337%(见表2)。在同质性指标上,海绵设施规划将缓解现状的土壤平均侵蚀情况和入湖污染物总负荷。

　　 为揭示规划影响的区域,本研究将比较规划前后ASL和NASL的空间变化。对ASL/NASL采用栅格计算“规划后的值-规划前的值”得到规划对土壤侵蚀强度和沉积物输出流量的格局影响图(见图4);其中,正/负值处代表规划对现状产生积极/消极影响的地区。根据结果,规划后ASL减小处占基地总面积15.91%,显著高于增大处(面积比例2.91%),表明规划所导致土壤侵蚀源头减弱的地区大于增强的地区;规划后NASL减小处占基地总面积14.93%,略微低于增大处(面积比例19.29%),反映规划所导致沉积物输出增加的区域将略大于减小的区域(见表3)。

　　 综上,规划对沉积物输出产生的正面影响存在空间分布局限的问题,但在土壤侵蚀平均强度、入湖沉积物负荷、土壤侵蚀源头削减等层面显著优于现状。

　　 上述模拟研究验证了规划情景下海绵设施布局的实效,揭示了未来LULC变化将如何改善现状水土流失情况。随着城市的无序增长,LULC的变化将对与其关联的功能格局产生剧烈而不可预期的影响,这种不确定性是单向、线性传统规划方法无法把控的,会最终影响到对现状问题响应决策的实效性^[27] 。本文引入了情景规划概念,将现状和规划视作两种不同阶段的模型输入;这种方法论有利于提升海绵城市规划对满足既定期望和目标的可靠性。

　　 研究已证实了土地干扰类型和可蚀性、运移容量^[28] 、营养输出^[29] 、污染负荷^[30] 和水文系统衰退^[9,25,31] 间的互动关联。为此,规划前期就应基于整合过程探究不同阶段土地变化的协同效应;规划策略应遵循系统机制,至少应考虑LUCC对功能要素的最终影响^[5,9,31] 。由于子流域系统的级联性和耦合性^[20,27,29] ,LULC对水土流失过程的复合作用依然是模糊逻辑^[5,20] ,这将增加保护和治理决策的难度;如果海绵城市方案没有经历进一步评估和验证,仅根据现状问题制定的规划响应将难以确保实现优化措施的期望、目标和意义^[15,17] 。

　　 和其他决策法不同,情景规划具有考虑所有可能过程和情况的优点^[27] 。这种理念简化了存在于系统间和系统成分间多级、动态的交互作用,优先考虑未来的预景结果而非不可预知的突变和剧变。在空间决策早期阶段引入这一方法论能确保海绵设施在水土流失总量、沉积总量和改进地区面积方面的最低功效^[14,26] 。对任何模糊系统的优化决策,情景规划都宜成为相关量化程序的重要理念。

　　 在水文模型半个世纪的发展历程中出现过诸多对特定物理过程的概化^[17] 。本研究的首要目的是量化水土流失对水环境的影响,因此引入污染物削减率作为不同类型和规模设施的功能指标,以此校正模型输入参数。该方法将污染物的削减率等同于“功能性”水土流失的削减率,从而和水土防治措施的效用挂钩;然而,这也可能简化了海绵设施的关键机理。

　　 首先,污染物削减率和“功能性”水土流失削减率的等量关系是基于污染物(SS、P、N等)的土壤比例在研究范围内空间同质的假设。由于研究区域缺乏更详细的土壤参数,这些本底值在每个HRU内都假设一致。但在二维RUSLE模拟中,污染物行为往往基于其在土壤中的本底数据^[5,20] 。另外,海绵设施的功能被限定在土壤或污染物的源头控制,污染物运移过程的变化被抽象为源头的下游反应。但在一些污染行为建模中,地表径流基本的行为变化仍是影响总污染负荷的重要因素^[14,26] ;一些LULC模型的参数量化防治措施对土壤运移的作用(例如运输承载量)^[17,19,20] 。为此,需进一步细化污染物本底值、源头和中间过程的建模。

　　 由于受时空变化等复合因素影响,水土流失及其污染的模拟是极端复杂的^[6] ,使在不同土壤和LULC作用下的水土流失、径流总量、营养滞留和输运难以精准预测^[32] 。海绵设施对污染物的消减量还取决于其他因素,例如:过滤和吸纳设施所在的地貌、水文、植被-土壤类型、季节^[32,33] ;海绵设施的构型、结构、宽度、厚度^[26,33] 和植被种类^[15] ;地上设施的水力粗糙度和地下设施在多重过程影响下的复合流行为^[3,20,22] ;和雨滴动量、化学扩散、渗滤过程、温度、土壤湿度等有关的水文-化学过程^[2,34] 。另外,土壤的粒径、粗糙度、颗粒成分有时也会影响到治理措施的消减量^[35] 。对于这些相关参数的建模都值得进一步研究。