2018, Vol. 30 
(2: 鄱阳湖湿地与流域研究教育部重点实验室, 南昌 330022)
(2: Key Laboratory of Poyang Lake Wetland and Watershed Research, Nanchang 330022, P. R. China)
非点源污染与流域水文过程密切相关, 岩溶流域含水系统由于连通地表的落水洞等垂直管道将近水平的地下暗河联系起来，降水及其形成的地表径流可以通过这些垂直管道迅速灌入地下河系，从而改变了水及其所携带的非点源污染物质在垂直与水平方向的传输速度与数量，使岩溶流域内地表－地下之间的物质交换与传输过程变得比较复杂[1-2].应用分布式机理性模型对非点源污染进行数值模拟，是非点源研究中的重要方法[3-5]，但应用基于松散介质特性建立的分布式水文模型难以模拟岩溶流域的水文过程和与此相关的非点源污染过程[6-7]，应用广泛的SWAT模型在模拟岩溶地区的水文、水质时会存在一些不足与局限[8-9].因此，针对岩溶含水介质的特征，对现有的分布式非点源污染模型进行修正，建立适合于岩溶流域的分布式非点源污染模拟模型，对于开展非点源污染形成与传输的机理研究有着重要的作用.
SWAT模型广泛地应用于各流域的非点源污染研究，而岩溶地区的非点源污染的研究相对较少，由于SWAT模型是基于松散介质特性建立的分布式模型，缺乏必要的组件来刻画岩溶含水系统的特征，因而在模拟岩溶地区的水文问题时存在诸多局限.如Malagò等[10]应用SWAT模型对岩溶流域的基流进行模拟，结果表明模拟结果有较大的误差.在水质与非点源污染模拟方面，Amatya等[11]在应用SWAT模型于Chapel Branch Creek岩溶流域以评估河流的月流量动态变化时，也发现对岩溶泉流量的模拟存在高估或低估的现象，同时也表明SWAT在模拟岩溶特征流域的营养物质和沉积负荷时，应根据岩溶特征作适当的调整.就特殊的水文循环特点，国内有研究者针对这些局限对SWAT进行过一定的修正，如任启伟[12]在修正SWAT模型的基础上建立了双重尺度的岩溶流域分布式水文模型，其结构包括水循环的陆面部分和水面部分，考虑了岩溶流域的特殊结构，采用指数衰减方程刻画表层岩溶带的调蓄过程，应用线性水库刻画浅层岩溶裂隙网络的调蓄作用，并用马斯京根法概算地下河汇流过程；郑捷等[13]针对平原型灌区人工—自然复合的水文循环特点，考虑平原灌区灌溉渠道、排水沟和河道等人工干扰，在沟渠河网的提取方法、子流域与水文响应单元的划分以及作物耗水量计算模块等方面对SWAT模型进行了改进.因此，应用SWAT模型开展岩溶流域的非点源污染模拟，仍然是值得研究的问题.
我国的岩溶地区面积达3.44×106 km2，占国土面积的35.9 %，其中大部分集中在我国南方地区[14], 该地区从元古代至新生代地层发育齐全，隆起与拗陷复杂；构造运动、岩浆活动、沉积作用、变质作用具多旋回性、多阶段性和不平衡性，产生不同的碳酸盐可溶岩类构造；可溶岩地层与非可溶岩地层相间产出，隐伏岩溶发育并成带分布, 而亚热带湿润气候致使地表河湖水系发达，地下水网沟通，使不同时期的可溶岩体被内外地质营力不断改造，形成了典型的亚热带岩溶地貌.江西省岩溶分布具有“三带”和“三块”的特征，即瑞昌—彭泽发育带、萍乡—乐平发育带、崇义—宁都发育带和上饶发育块、吉安发育块及龙南发育块.横港河流域属于瑞昌—彭泽岩溶发育带，该带岩溶主要在奥陶系、石灰系、二迭系地层发育，石灰系、二迭系岩溶发育较好，峰丛、石芽等岩溶特征发育明显，地表呈现不同程度的植被覆盖；此外，天坑、暗河等岩溶构造众多，其中峨眉洞穴群含大小溶洞20余个，规模较大，是研究我国中东部岩溶地区水文特征的理想场所，以横港河流域岩溶为试验区，开展岩溶背景下非点源污染模拟的研究，具有基础意义和典型意义.
本文选择横港河流域的一个子流域作为典型试验区，以岩溶地区的非点源污染模拟为研究对象，针对岩溶系统的非均匀含水介质，引入落水洞、伏流/暗河、岩溶泉的水文过程及主要营养盐输移过程，修正SWAT模型原有的水文循环过程及相关算法，改变其只适用于松散均匀介质流域非点源污染模拟的单一特性，通过基于植被岩溶比重指数的岩溶流域覆盖信息提取、基于遥感影像/DEM的岩溶地表与地下水流向信息提取、岩溶地表径流水文曲线数等岩溶地貌特征参数的试验与研究，结合实地水文水质数据采集，建立了适合于岩溶流域的非点源污染模拟和相应的模拟方法，并应用修正模型，定量评估了落水洞、伏流/暗河等岩溶特征对氮、磷等主要非点源污染输移的影响及其带来的时空效应；探讨了落水洞、伏流/暗河等岩溶特征对地表—地下水文与营养盐交互作用及输移机理.
1 研究区概况及研究方法 1.1 研究区概况研究区(29°29′36″~29°32′59″N, 115°25′30″~115°30′25″E)位于横港河流域的子流域，横港河是长河上游的一条支流，也是长江在江西省境内的一条支流.研究区集水面积约27.63 km2，在行政上隶属于江西省瑞昌市和德安市面积，其中瑞昌市面积27.17 km2，德安市面积0.46 km2, 高程分布在54~562 m之间，平均高程为219 m，高程变差为106 m；该流域地处九岭东西向构造带－阳新东西向构造拗褶带东段，东西向构造甚为发育；地貌以低山丘陵为主，地势由西南向东北倾斜；岩溶发育强烈，拥有蛇狮洞、银丝洞、下畈洞、消水洞等大型溶洞，具有落水洞、地下暗河、涌泉、石芽、溶沟等岩溶形态特征，是峨嵋溶洞群的最主要部分，有关落水洞、地下暗河出入口的位置见图 1.
|
图 1 研究区位置及地形 Fig.1 Location and terrain of study area |
流域内气候温和，四季分明，属大陆温湿性气候带，年平均气温17.5℃，年降雨量1700 mm左右，年日照时数2000 h上下，年无霜期240~260 d.年均降雨量1614.3 mm，最大降雨量2180.3 mm(1998年), ≥100 mm暴雨日年平均1~3 d.每年4—8月降雨量占全年降雨量的63.4 %.按土壤发生原理，流域内土壤采用土类、亚类、土属、土种四级分别可归纳为6个土类，6个亚类，14个土属，19个土种，其中6个土类分别为水稻土、潮土、黄棕壤、红壤土、红色石灰土和石灰(岩)土；森林植被有常绿阔叶林、竹林、针阔叶混交林、常绿与落叶混交林、落叶阔叶林、暖性针叶林和灌丛七大类，森林覆盖率达68.3%.
流域内有2个小型村落，居民不足500人，以老人和小孩为主，青壮年外出务工；除了有少量中稻种植外，主要以玉米、大豆、油菜、棉花、西瓜和红薯等旱地作物为主；此外，流域内还少量的柑桔种植；流域内无点源污染.
1.2 数据准备及方法本文的数据及其处理包括两部分，其一为SWAT模型原始的基本数据需求，其二为针对落水洞、地下暗河、岩溶泉等主要岩溶特征及其水文过程，进行SWAT模拟的扩展数据需求.根据SWAT模型基本输入数据的说明文档，其输入数据主要有：地形DEM、河网、土地利用与覆盖、土壤分布与理化属性、水文与水质数据、气象数据，其中水文与水质数据主要用于模型的建立、模型参数的率定和验证.
1.2.1 岩溶流域岩溶覆盖信息提取及植被岩溶比重指数土地覆盖数据是SWAT模型的主要输入变量，同时也是针对落水洞这一岩溶特征及其水文过程进行SWAT模拟的重要参考数据.为了获得最新的土地利用/覆盖数据，本文采用了Landsat 8的数据作为遥感数据源，选用2013年8月9日、2014年10月15日和2015年1月3日OLI陆地成像仪(Operational Land Imager).结合通过Google Earth上的IKONOS 1 m分辨率遥感影像数据，通过目视解译方法对横港河流域土地利用/覆盖进行遥感解译，并对解译结果进行了多次调研并验证[15].
由于落水洞的水文过程和营养盐输移过程的定义和模拟均在水文响应单元HRU内进行，而SWAT模型HRU的划定是依据下垫面的土地利用/覆盖及土壤的理化属性来进行的.因此在流域的空间离散化过程中，若需将落水洞单独划分为一个HRU，必须对土地利用/覆盖类型和土壤的理化属性有新的定义[16]，其中土地利用/覆盖类型按计划将增加一种新的土地利用/覆盖类型，即落水洞，并在一级名称下增设5种二级名称，即林地落水洞、草地落水洞、城区落水洞、耕地落水洞和裸地落水洞，但实际上研究区的落水洞只有在清华村的蛇狮洞一个，其二级类型属于林地落水洞，为此，在解译的土地利用/覆盖类型中增加了单一的新类型，即林地落水洞.形成的最终类别为：水田、旱地、有林地、灌木林、疏林地、中覆盖度草地、河渠、农村居住用地和林地落水洞共9个类别(图 2).
|
图 2 研究区土地利用与覆盖解译的结果 Fig.2 Interpretation results of land use and land cover in the study area |
植被岩溶的比重可以反映出地表岩溶的发育程度，并进一步分析地表岩溶特征是如何影响营养盐产出的，对文章之后的非点源污染的机理研究有重要的意义.基于上述的遥感数据，获取NDVI、FV、LST以及卷云波段，基于FV、LST以及卷云波段提取第一主成分PFLC，并构建PFLC-NDVI特征空间[15].在此基础上，提出能够反映研究区岩溶覆盖信息的植被岩溶比重指数(VKPI)，利用分级阈值法对VKPI值进行等级分值生成VKPI分级图，结合野外实地调查获取的岩溶数据进行分析及验证.研究区的VKPI指数分布在0~0.53之间.该指数越大，表示地表岩溶特征越明显、植被覆盖越少(图 3).
|
图 3 计算的VKPI指数 Fig.3 Calculated VKPI index |
我国南方亚热带岩溶的突出特征是地表有不同程度的植被覆盖，对岩溶流域的水文过程带来重要影响[17]，对地表径流水文曲线数(CN)的影响就是其中之一，且SWAT模型对地表产流的计算是通过SCS-CN方程来计算的, 因此对研究区进行CN的研究具有重要的意义[18-19]，而对研究区进行岩溶地貌土地利用/覆盖的信息提取是进行CN试验研究的重要基础.
在研究区典型岩溶地貌下，选取1 km2范围内，完全岩溶、较多岩溶、较少岩溶和没有岩溶4种类型、4种坡度的试验样区，建立4个面积为1.5 m×1.5 m的地表径流样区，并在实验区内架设WatchDog自动气象站一台监测降雨量.通过地表径流监测试验，获取研究区的降水和地表径流资料；利用SCS-CN模型的基本原理得出CN值计算公式，利用反推的方法计算地表径流试验样区的CN值；利用遥感反演的研究区典型岩溶地貌下的土地利用/覆盖类型、坡度、土壤特性，估算整个研究区岩溶区域的CN值.由每个试验样区的实验数据，并根据SCS-CN方程，估算每个样区的CN值(表 1)，得到了岩溶区域不同类型的初步CN值.这为下一步估算研究区不同岩溶土地利用/覆盖、地形坡度和土壤属性的CN值提供了基础.
|
表 1 试验样区初步得出的CN值 Tab.1 The initial value of CN in the test sample area |
现有的地形数据(DEM)在SWAT模型内难于提取岩溶流域地表水流方向的细节[20]，且岩溶流域内地表复杂，裂隙、洞穴迂回曲折，纵横交错，将极大地影响到水及其所携带的非点源污染物质在垂直与水平方向的传输速度与数量.以DEM线性体为主，遥感影像线性体为辅，利用野外获取的数据对研究区的地下水流向进行验证，获取长河流域DEM和遥感影像的线性体[21].结合江西省地质部门2001—2002年完成的1:50000地质空间数据库资料，最终确定研究区地表及地下水流向的基本信息.
1.2.4 气象、水文水质数据的获取气象数据源自研究区周边的九江、瑞昌、星子、武宁和德安5个气象站的逐日观测数据，要素包括最高气温、最低气温、平均气温、相对湿度、降水、云量和风速，太阳辐射由SWAT自带的模型推求；雨量与流量数据来源于水文部门的横港、曾家山、下毕、铺头和下毕5个水文站的逐日实测数据.由于研究区空间范围较小，上述气象站和水文站位置都在研究区范围之外，所以这些站点的位置难于在相关图中标示.
水文水质数据主要用于SWAT模型的参数率定和模型验证，包括河流流量数据和河流主要营养盐浓度数据.为了对落水洞和地下暗河等岩溶流域含水系统等的水文过程和营养盐输移过程进行监测，设置了8个测点(图 4、表 2)，对流量与水质进行8次测定.采样要素包括水深、流速、总氮、总磷、氨氮、硝酸盐氮和亚硝酸盐氮等.将超声波水深探测仪所测得的水深数据转换为断面资料，结合每次采样所获得实测流速即可得到实时的流量数据.
|
图 4 流量和水质采样点的分布 Fig.4 Distribution of flow and water quality sampling points |
|
|
表 2 流量和水质采样点名称及其相关信息 Tab.2 Flow and water quality sampling points and their related information |
土壤基本数据主要包括土壤类型空间分布图、土壤理化属性数据，来源于中国科学院南京土壤研究所，其1: 50000矢量数据(shapefile格式)的主要字段包括土类、亚类、土属和土种，其中图斑名称备注了非土壤区域(居民区、湖泊、河流等)类型(图 5)；栅格数据文件包括有机碳含量、砂粒含量、粉粒含量、黏粒含量、全氮含量、全磷含量、速效磷含量等理化属性，均包含0~20、20~30和30~100 cm土层深度.
|
图 5 土壤类型空间分布 Fig.5 Spatial distribution of soil types |
SWAT模型使用的土壤数据包括物理属性和化学属性，物理属性存放在SWAT用户数据库(usersoil)中，一共有17个属性[22].我国的土壤质地标准与美国的土壤质地标准不同，中美间的土壤质地标准将根据朱秋潮等[23]的方法进行转换；土壤容重、有效田间持水量、饱和导水率等参数将根据数字化土壤资料的土壤空间分布及理化属性数据由软件SWAT来计算.
1.3 研究方法在上述相关数据、信息提取及获取的基础上，针对岩溶流域特殊的水文条件，引入落水洞和地下水的水文过程及氮磷营养盐的输移过程，修正SWAT模型原有的水文循环及相关算法，并进行模型的有效性分析.在修正模型建立以后，采用了汇流累积区面积为54 hm2的阈值进行了空间离散化，生成了29个子流域，并采用土地利用/覆盖、土壤、坡度比例阈值为20:10:20的水文响应单元(HRU)划分阈值，得到736个HRUs[24].模型参数率定先对水文参数进行率定，率定完水文参数以后，再对营养盐负荷参数进行率定，这期间不再对水文参数进行率定.使用的模型均为原始模型.所有参数率定完成之后，模型模拟结果作为模型对比和结果分析的对照(记为T1).然后采用控制性方法，控制所有变量保持不变，对修正模型所引入的新变量δgw_karst(岩溶地区地下水补给滞后系数)、sink(落水洞水量分配系数)和SS(地下暗河的营养盐分配系数)参数进行设置，对地下暗河所涉及到的HRU中的土壤水力传导系数调整到与河流的水力传导系数相同.设置了这些反映岩溶特征的参数之后，再对模型进行模拟, 所得结果记为T2.
前后的模拟过程中，输入变量除了岩溶特征变量之外，其余变量都相同，由此可以说明，T1和T2的结果差异完全是由于岩溶变量的引入所导致的.因此T1和T2的模拟结果对比，可以定量评估和分析岩溶流域对非点源污染的影响，并进行相关的形成机制分析.
2 模型的修正方法和有效性分析 2.1 SWAT模型修正方法 2.1.1 增加落水洞水文过程和营养物质输移过程落水洞的水文过程定义并计算在HRU内，为了表示落水洞的水文过程，引入新的变量wsink和sink，其中sink为落水洞水量分配系数，其值在0~1之间，通过sink可以判断HRU是否存在落水洞，如果存在落水洞，则sink > 0，否则sink=0；wsink为HRU内伏流与落水洞的渗透量，其计算公式为：
| $ {w_{{\rm{sink}}}} = \left( {{w_{{\rm{surf}}}} + {w_{{\rm{lat}}}}} \right) \cdot s\mathit{in}k $ | (1) |
式中，wsurf为地表径流流入落水洞的水量(mm)，wlat为通过侧流流入落水洞的水量(mm)，这2个变量的计算由SWAT模型原有的算法进行. sink值的大小通过落水洞汇水面积与所在子流域面积占比来确定.增加落水洞水文过程以后，子流域所有HRU的水文过程计算存在两种可能，一是该HRU为非落水洞(sink=0)，则HRU内的所有算法仍保留SWAT模型在HRU内的算法，见公式(1)；二是该HRU为落水洞(sink > 0)，此时对SWAT原有算法进行修正，其算法如下：
| $ w{\left( i \right)_{{\rm{rchrg\_karst}}}} = \left[ {1 - \exp \left( { - 1/{\delta _{{\rm{gw\_karst}}}}} \right)} \right] \cdot {w_{{\rm{sink}}}} + \exp \left( { - 1/{\delta _{{\rm{gw\_karst}}}}} \right) \cdot w{\left( {i - 1} \right)_{{\rm{rchrg\_karst}}}} $ | (2) |
式中，w(i)rchrg_karst是时间为i(日)岩溶地区非承压含水层的补给量(mm)；δgw_karst为岩溶落水洞的滞后时间(h)，其值可以通过率定或采用对岩溶泉的实测数据来估算，本文实际操作中，按δgw_karst=δgw·VKPI·cosθ来计算，其中VKPI为植被岩溶比重指数，θ为坡度. w(i-1)rchrg_kars是时间为i-1(日)的岩溶地区非承压含水层的补给量.而承压含水层HRU的日水量由下式计算：
| $ {w_{{\rm{deepst}}}} = {w_{{\rm{gwseep}}}} + \left( {1 - \mathit{sink}} \right)\left( {{w_{{\rm{surf}}}} + {w_{{\rm{lat}}}}} \right) $ | (3) |
式中，wdeepst为承压含水层HRU的日水量(mm)，wgwseep为HRU每日补给深含水层的水量(mm)，(wsurf+wlat)为HRU经由落水洞的每日水量(mm).此外，还要修正的算法是HRU流入主河道的日水量wdr计算方法，当sink=0时，该HRU为非落水洞HRU，wdr的计算公式如下：
| $ {w_{{\rm{dr}}}} = {w_{{\rm{day}}}} + {w_{{\rm{lat}}}} + {w_{{\rm{gw}}}} + {w_{{\rm{tile}}}} $ | (4) |
式中，wdr为从HRU流入主河道的每日水量(mm)，wday为从HRU流入主河道的地表径流每日水量(mm)，wlat为HRU的每日总侧流量(mm)，wgw为HRU对地下水贡献的每日水量(mm)，wtile为HRU内由土壤层人工排水管排放的日水量(mm).
当sink > 0时，该HRU为落水洞HRU，wdr不再计算wday和wlat，HRU的地表径流传输损失也不被模拟.
在岩溶环境下，经由落水洞补给到非承压含水层的营养盐量，应该概化等于该落水洞域被地表水和侧流输移的营养盐量的总和.为此，引入两类新变量Nchrg_sepbtm r和Nrchrg_karst，分别表示非岩溶HRU内通过土壤渗透补给含水层的营养盐日负荷和岩溶HRU补给含水层的营养盐日负荷(kg)，具体的营养盐(如氮(N)、磷(P)等元素)不再细述.当sink=0时，Nrchrg_sepbtm被计算；当sink > 0时，Nrchrg_karst被计算.它们的计算公式为：
| $ N{\left( i \right)_{{\rm{rchrg\_sepbtm}}}} = \left[ {1 - \exp \left( {\frac{{ - 1}}{{{\delta _{{\rm{gw}}}}}}} \right)} \right] \cdot {N_{{\rm{perc}}}} + \exp \left( {\frac{{ - 1}}{{{\delta _{{\rm{gw}}}}}}} \right) \cdot N{\left( {i - 1} \right)_{{\rm{rchrg\_sepbtm}}}}, 当\mathit{sink} = 0 $ | (5) |
| $ N{\left( i \right)_{{\rm{rchrg\_karst}}}} = \left[ {1 - \exp \left( {\frac{{ - 1}}{{{\delta _{{\rm{gw\_karst}}}}}}} \right)} \right] \cdot {N_{{\rm{sink}}}} + \exp \left( {\frac{{ - 1}}{{{\delta _{{\rm{gw\_karst}}}}}}} \right) \cdot N{\left( {i - 1} \right)_{{\rm{rchrg\_karst}}}}, 当\mathit{sink} > 0 $ | (6) |
式中，i为时间变量(天)，Nperc和Nsink分别为HRU内由渗透输移的营养盐日负荷量和由落水洞岩溶特征输移的营养盐日负荷量(kg). Nsink=(Nsurf+Nlat)·sink，分别表示由地表径流和侧流输移至落水洞岩溶特征的营养盐负荷. HRU内补给含水层的总营养盐日负荷由下式计算：
| $ {N_{{\rm{rchrg}}}} = {N_{{\rm{rchrg\_sepbtm}}}} + {N_{{\rm{rchrg\_karst}}}} $ | (7) |
由于伏流/暗河将河道里的水直接传输到非承压含水层，然后以基流的方式通过岩溶泉排放到主河道，因此伏流/暗河水文过程可以假定为河床具有高水力传导度的支流来表示，通过修正SWAT的“.rte”中的CH_K(2)参数来实现.其传输损失主要影响非承压含水层.在SWAT模型中支流的模拟是在子流域内进行，因此在模拟支流(伏流/暗河)时，非岩溶和岩溶的伏流/暗河是连通的.为此引入伏流/暗河的营养盐分配系数SS，由地表径流的传输损失来确定.其计算公式为：
| $ SS = {w_{{\rm{loss}}}}/{w_{{\rm{surf}}}} $ | (8) |
式中，wloss为伏流/暗河的传输损失量(mm)，wsurf的意义同前.
由伏流/暗河传输而进入到非承压含水层的营养盐日负荷量由下式计算：
| $ {N_{{\rm{sep\_direct}}}} = SS \cdot {N_{{\rm{surf}}}} $ | (9) |
式中，Nsurf为HRU内地表径流的日负荷量，伏流的营养盐传输损失用于含水层营养盐补给负荷的计算，计算式见公式(6).由HRU地表径流进入主河道的营养盐负荷计算按下式进行：
| $ {N_{{\rm{sub\_surf}}}} = \left( {1 - SS} \right) \cdot {N_{{\rm{sep\_direct}}}} \cdot {\xi _{{\rm{hru\_bafr}}}} $ | (10) |
式中，ξhru_bafr为HRU与流域面积的比例系数，其余同前.
2.2 模型的有效性分析 2.2.1 流量与营养盐参数的敏感性模拟及分析模型的参数率定采用SWAT-CUP软件进行，其参数的敏感性测试是采用t-Stat方法.流量的参数敏感性试验包括ALPHA_BF、GW_DELAY、GWQMN、GW_REVAP、ESCO、SOL_AWC等12个参数，各参数名称、所在数据库文件、物理意义和率定结果见表 3[25].
|
|
表 3 流量的率定参数及率定结果* Tab.3 Calibration parameter and calibration result of flow |
由t-Stat方法可知，P-Value < 0.05时，该参数认为是敏感的. 图 6给出了这12个参数的率定敏感程度，从P-Value值可以看出，在12个率定参数中，CH_K2、ALPHA_BF这2个参数比较敏感.
|
图 6 与流量相关的12个参数率定的敏感程度t-Stat和P-Value值 Fig.6 The sensitivity level of 12 calibration parameters related to flow of t-Stat values and P-Value |
氮(N)营养盐负荷的参数率定选择了ERORGN、SPCON、BIOMIX和SHALLST_N等8个参数，这8个参数的名称、所在数据库文件及物理意义见表 4.率定结果表明，在这8个参数中，ERORGN、CH_ONCO和NPERCO这3个参数对N营养盐的排放具有一定的敏感性，其中以NPERCO的P-Value为最小(图 7).各参数的率定结果见表 4.
|
图 7 与氮相关的8个参数率定的敏感程度t-Stat值和P-Value Fig.7 The sensitivity level of 8 calibration parameters related to nitrogen of t-Stat values and P-Value |
|
|
表 4 氮相关参数的率定结果 Tab.4 Calibration result of nitrogen related parameters |
磷(P)营养盐负荷的参数率定选择PSP、RSDCO、PHOSKD和PPERCO等10个参数，结果表明这10个参数中，只有ERORGP参数反应敏感，其P-Value小于0.05(图 8). 图 8给出了反映10个参数敏感程度的P-Value值和t-Stat值，各参数的名称、所在的数据库文件、物理意义和率定结果见表 5.
|
图 8 与磷相关的10个相关参数率定的敏感程度t-Stat值和P-Value Fig.8 The sensitivity level of 10 calibration parameters related to phosphorus of t-Stat values and P-Value |
|
|
表 5 磷相关参数的率定结果 Tab.5 Calibration result of phosphorus related parameters |
在N和P参数率定的基础上，对上述反应敏感的参数再利用SWAT-CUP进行进一步的微调率定，得到最优的参数，并利用这些率定好的参数进行模拟.在实测N和P负荷数据基础上对模拟结果进行验证，并分析模拟的有效性程度.
2.2.2 流量与营养盐模拟结果的验证及有效性分析使用SWAT-CUP参数率定的结果，用修正模型对2010—2013年进行模拟，并提取模型输出的2013年流量和营养盐负荷，利用瑞昌市水文站的实测流量和6次4个水质采样点的水质采样数据，进行模拟有效性分析. 图 9给出了修正模型模拟流量与实测流量及流域平均降雨量的对照，可以看出模拟与实测的流量基本趋势一致，但模拟的流量峰值与实测的流量峰值存在较大的出入，与实测值比较，模拟值具有峰值偏大的趋势.与流域平均降雨量比较可以看出，模拟值的变化趋势与流域降雨量的变化趋势有比较好的一致性.
|
图 9 修正模型模拟的流量与实测流量及流域平均降雨量对照 Fig.9 The comparison of watershed average rainfall between simulated flow and the measured flow |
修正模型模拟流量的有效性分析可以看出，模拟值与实测值的相关系数为0.94，达到极显著水平(表 6).此外模型的NS系数为0.89，说明修正模型在模拟岩溶流域的产流方面有效性得到了较好的改善[26].
|
|
表 6 修正模型模拟的流量有效性分析 Tab.6 The flow effectiveness analysis of correction model simulation |
在模拟营养盐负荷方面，图 10给出了修正模型模拟与实测TP、TN浓度对照，可以看出模型在模拟TP和TN负荷方面，时间流与空间流趋势都基本一致. 表 7给出了修正模型模拟营养盐方面的有效性分析，其中OrgN、OrgP、NO3、NH4、MinP分别为有机氮、有机磷、硝态氮、氨氮、矿物磷.从R2值和NS值来看，模拟效果略差于流量的模拟，TN和TP的模拟有效性指数R2和NS系数均在0.7以上，且TN的模拟效果比TP的模拟效果差.其模拟TN和TP的R2及NS系数分别为0.73、0.86和0.71、0.82，其余营养盐的有效性系数见表 7.
|
图 10 修正模型在TN和TP方面的实测值与模拟值对照 Fig.10 The comparison of measured values and simulated values of TN and TP about correction model |
|
|
表 7 修正模型模拟营养盐的有效性分析 Tab.7 The effectiveness analysis about correction model to simulation nutrient |
岩溶特征对非点源污染的影响最终体现在河道营养盐负荷的变化上面[27]，为了定量分析这些影响，在前述获取的水文水质数据中，选取2个地下暗河的出口河段作为观测点，对比T1和T2的模拟结果.这2个地下暗河的出口点分别为老社洞和洞下洞，其涉及的河段编码为18和15，其中18河段是2个落水洞的所在河段.对2个地下暗河出口河段洞下洞、老社洞的TN和TP变化百分比进行分析(图 11和表 8).
|
图 11 洞下洞、老社洞两个地下暗河出口河段的TN和TP变化百分比 Fig.11 The percentage change of TN and TP in two underground river outlets(Dongxiadong and Laoshedong) |
|
|
表 8 洞下洞、老社洞两个地下暗河出口河段不同月份TN和TP变化百分比(%) Tab.8 The percentage changes of TN and TP in different months in two underground river outlets |
18和15河段TN和TP均是增加的，其中18河段的增加比例要大于15河段的增加量，且TP的增加比例均大于TN. 15和18河段TN和TP的平均增加比例分别为0.41%、0.79%和0.86%、2.12%(表 8).由此可以看出，岩溶特征变量的引入导致的非点源负荷是增加的，且TP的增加比例大于TN的增加比例.
3.1 植被岩溶比重指数对水量及营养盐产出的影响分析植被岩溶比重指数VKPI反映了地表植被覆盖与地表岩溶出露程度的定量比例关系，强调的是地表岩溶发育程度，不同植被覆盖对营养物质数量以及氮磷形态的影响不同[28]. 图 3给出了研究区的VKPI指数，通过提取研究区各子流域的VPKI，并与模拟的各子流域水量(mm)进行相关分析，结果发现VKPI与子流域的地表产流量(SurQ)有极好的相关关系(图 12a)，与侧流(LATQ)有弱的负相关关系(通过α=0.1的检验)(图 12b)，而与地下水产流量(GWQ)则没有相关性(图 12c).由此可以说明，地表的岩溶特征会导致地表产流的增加和侧流的减少，而对地下水的产出则没有影响.这其中的原因可能与引入的岩溶特征变量只考虑了落水洞、地下暗河等大型的岩溶构造有关.
|
图 12 子流域产水量与VKPI的关系 Fig.12 The contact of sub basin water-output and VKPI |
VKPI指数与子流域营养盐产出量(kg/hm2)的关系分析，则表明VKPI与可溶性磷(SolP)、有机磷(OrgP)有相对较好的相关关系，其相关系数为0.42~0.43，与有机氮、地表产流中硝酸氮(NSurQ)和沉积磷(SedP)的相关系数的则较小，其值变化在0.20~0.39之间，而与地下水中的硝酸氮(GWNO3)关系极小(图 13).由此表明VKPI与子流域的营养盐产量存在一定的正相关关系，反映了地表岩溶特征对子流域可溶性磷及有机磷产出有一定的增加作用，对有机氮、地表产流中硝酸氮和沉积磷有微弱的增加作用，而对地下水中的硝酸氮产出则没有明显的作用.
|
图 13 子流域营养盐产量与VKPI的关系 Fig.13 The contact of sub basin water-output and nutrient production |
岩溶流域降水及其形成的地表径流可以通过岩溶地区垂直管道迅速灌入地下河系，从而改变了水及其所携带的非点源污染物质在垂直与水平方向的传输速度与数量.同时，这些大型构造也会带来相应的岩溶裂隙，这些岩溶裂隙同样会改变水及营养盐的输移方式[29-30].为此，将子流域的降雨量(mm)与典型岩溶子流域产水量(mm)及营养盐产出(kg/hm2)差异进行对比分析，结果表明落水洞、地下暗河等岩溶特征变量引入后，落水洞、地下暗河所在的子流域产水量及营养盐产量均有变化，而没有这些大型构造的子流域产水量及营养盐产量均无变化.这些存在大型构造的子流域其产水量及营养盐产量变化与降雨量有明显的关系. 图 14表明了含有大型岩溶构造的不同子流域地表产水量变化与降雨量的相关系数均值大于0.8，且是负相关关系；不同子流域之间降雨量的增加对地表产流量的减少量是不同的，其中以17子流域的减少量为最大，而17子流域正是落水洞所在地，其次是21子流域，是地下暗河的入口，入口处有小型的落水洞.
|
图 14 子流域地表产流变化与降雨量的关系 Fig.14 The contact of sub basin surface runoff and rainfall |
降雨量与子流域营养盐产出变化呈正相关关系，表明由于落水洞及地下暗河的存在导致降雨量越大，营养盐产出的增量就越大. 图 15表明落水洞所在的子流域17与21，无论OrgN或OrgP的增量都随降雨量的增大而增大，且落水洞越大型，其增量也越大，如17子流域的OrgN和OrgP增量分别在0~0.74和0~0.30 kg/hm2之间变化，而21子流域的OrgN和OrgP增量分别在0~0.28和0~0.13 kg/hm2之间变化. 28子流域有地下暗河但没有落水洞，其降雨量与OrgN和OrgP的增量相关系数大于17和21子流域，但其OrgN和OrgP的增量分别在0~0.11和0~0.05 kg/hm2之间变化.由此表明落水洞是造成子流域营养盐增加的最主要因素.
|
图 15 降雨强度与子流域营养盐变化的关系 Fig.15 The contact of rainfall and sub basin nutrition change |
从子流域的分析可以发现，由于降雨强度对落水洞、地下暗河等岩溶特征所引起的侧流及由此引起的侧流氮、磷输出均无明显变化，其他营养盐增量与降雨强度的关系也不明显.
4 结论与讨论通过引入岩溶特征变量，对岩溶流域的非点源污染进行控制性模拟，对比引入岩溶特征变量前后的河道流量与营养盐负荷、流域的产流及营养盐产出量变化，分析了植被岩溶比重指数VKPI对水量及营养盐产出的影响和降雨强度对典型岩溶子流域水文及营养盐产出的影响，得出了如下主要结论：
1) 岩溶特征对非点源污染的影响是增加的，且TP增加比例大于TN增加比例，包含落水洞的河流TN和TP增加比例分别为0.86%和2.12%，明显大于没有落水洞河流TN和TP增加比例(分别为0.14%和0.79%).
2) 岩溶特征对非点源污染影响的形成机制分析表明，地表的岩溶特征改变了流域的产流方式及营养盐的输移方式，地表的植被岩溶比重指数VKPI的增加会导致地表产流的增加和侧流的减少，而对地下水的产出则没有影响；VKPI的增加也导致了流域营养盐产出成分的变化，其中与可溶性磷(SolP)、有机磷(OrgP)有相对较好的相关关系，与有机氮(OrgN)、地表产流中硝酸氮(NSurQ)和沉积磷(SedP)的相关系数则居其次，而与地下水中的硝酸氮(GWNO3)关系极小.
3) 岩溶特征对非点源污染影响的形成机制分析还表明，落水洞改变了降雨产流的方式，降雨强度与流域营养盐增量关系分析反映了落水洞的降雨产流使OrgN和OrgP增量分别在0~0.7和0~0.3 kg/hm2之间变化，远大于其他没有落水洞的流域；落水洞与地下暗河的直接连通，缺少了土壤层对水中营养成分的物理及生化过程，使流域OrgN和OrgP产出增加，是导致落水洞所在河段TN和TP增加的重要因素，也是河段TN、TP增加的主要营养盐成分.
4) 虽然针对岩溶流域的特征引入了反映岩溶地貌的水文及营养盐过程，但岩溶系统是一个复杂的系统，多种岩溶特征难于量化也难于表达，如岩溶流域含水系统的裂隙，就很难量化，而该变量又极大地影响了岩溶系统含水层对地表及地下水的影响.因此，模拟结果有关地下水方面的内容，包括地下水延迟时间的率定、植被岩溶比重指数VKPI对侧流及地下水的影响可能存在一些不合理或不确定之处，这有待进一步的研究.
| [1] |
Pannoa SV, Kelly WR. Nitrate and herbicide loading in two groundwater basins of illinois, Sinkhole Plain. Journal of Hydrology, 2003, 290: 229-242. |
| [2] |
Keith ES, Matthew H. Tile drainage as karst:Conduit flow and diffuse flow in a tile-drained watershed. Journal of Hydrology, 2008, 349(3/4): 291-301. |
| [3] |
Lai GY, Qiu L, Zhang ZY et al. Modification and efficiency of SWAT Model based on Multi-plant Growth Model. J Lake Sci, 2018, 30(2): 472-487. [赖格英, 仇霖, 张智勇等. 基于多植物生长模式的SWAT模型的修正与有效性分析. 湖泊科学, 2018, 30(2): 472-487. DOI:10.18307/2018.0219] |
| [4] |
Hao FH, Cheng HG, Yang ST. Non-point source pollution model-Theoretical methods and applications. Beijing: China Environmental Science Press, 2006. [郝芳华, 程红光, 杨胜天. 非点源污染模型——理论方法与应用. 北京: 中国环境科学出版社, 2006.]
|
| [5] |
Gassman PW, Reyes MR, Green CH et al. The soil and water assessment tool:Historical development, applications, and future research directions. Transactions of the ASABE, 2007, 50(4): 1211-1250. DOI:10.13031/2013.23637 |
| [6] |
Scanlona BR, Robert EM, Michael EB et al. Can we simulate regional groundwater flow in a karst system using equivalent porous media models? Case study, barton spring edwards aquifer, USA. Journal of Hydrology, 2003, 276: 137-158. DOI:10.1016/S0022-1694(03)00064-7 |
| [7] |
Xue XW, Chen X, Zhang ZC et al. The effect of karst fissures for the saturated flow of the slope collects. Water Resources and Power, 2009, 27(6): 20-23. [薛显武, 陈喜, 张志才等. 岩溶裂隙对坡面饱和水流汇集的影响研究. 水电能源科学, 2009, 27(6): 20-23. DOI:10.3969/j.issn.1000-7709.2009.06.006] |
| [8] |
Amatya D, Manoj J, Amy E. Application of swat model on chapel branch creek watershed with a karst topography. The American Society of Agricultural and Biological Engineers, 2011. |
| [9] |
Fernandeza GP, Chescheira GM, Skaggsa RW et al. DRAINMOD-GIS:A lumped parameter watershed scale drainage and water quality model. Agricultural Water Management, 2005, 81(2): 77-97. |
| [10] |
Malagò A, Efstathiou D, Bouraoui F et al. Regional scale hydrologic modeling of a karst-dominant geomorphology:The case study of the Island of Crete. Journal of Hydrology, 2016, 540: 64-81. DOI:10.1016/j.jhydrol.2016.05.061 |
| [11] |
Amatya D, Manoj J, Amy E. Application of swat model on chapel branch creek watershed with a karst topography. The American Society of Agricultural and Biological Engineers, 2011. |
| [12] |
Ren QW. Water quantity evaluation methodology based on modified SWAT hydrological modeling in southwest karst area[Dissertation]. Beijing: China University of Geosciences, 2006. [任启伟. 基于改进SWAT模型的西南岩溶流域水量评价方法研究[学位论文].. 北京: 中国地质大学, 2006.]
|
| [13] |
Zhen J, Li GY, Han ZZ et al. Application of modified SWAT model in plain irrigation district. Journal of Hydraulic Engineering, 2011, 42(1): 88-97. [郑捷, 李光永, 韩振中等. 改进的SWAT模型在平原灌区的应用. 水利学报, 2011, 42(1): 88-97.] |
| [14] |
Meng HH, Wang LC. The progress of hydrological models in karst basin. Progress in Geography, 2010, 29(11): 1311-1318. [蒙海花, 王腊春. 岩溶流域水文模型研究进展. 地理科学进展, 2010, 29(11): 1311-1318. DOI:10.11820/dlkxjz.2010.11.008] |
| [15] |
Liu W, Lai GY, Peng XJ et al. Information extraction of karst in karst watershed based on vegetation-karst proportion index. Jiangxi Science, 2016, 34(2): 194-222. [刘维, 赖格英, 彭晓娟等. 基于植被岩溶比重指数的岩溶流域岩溶信息提取. 江西科学, 2016, 34(2): 194-222.] |
| [16] |
Pei JG, Liang MZ, Chen Z. Ground water system division and its main characteristic value statistics in the karst rock mountainous areas in southwest. China Karst, 2008, 27(1): 6-10. [裴建国, 梁茂珍, 陈阵. 西南岩溶石山地区岩溶地下水系统划分及其主要特征值统计. 中国岩溶, 2008, 27(1): 6-10. DOI:10.3969/j.issn.1001-4810.2008.01.002] |
| [17] |
Wang XZ, Gan CY, Liang ZX et al. The study of the dynamic changes of vegetation coverage at Lianjiang River basin in the north of Guangdong Karst mountain area. China Karst, 2010, 29(4): 425-433. [王兮之, 甘春英, 梁钊雄等. 粤北岩溶山区连江流域植被覆盖度动态变化研究. 中国岩溶, 2010, 29(4): 425-433. DOI:10.3969/j.issn.1001-4810.2010.04.013] |
| [18] |
Chen JN. Land use/cover type of karst surface in change river basin and impacts on CN value[Dissertation]. Nanchang: Jiangxi Normal University, 2014. [陈静妮. 长河流域岩溶地表土地利用/覆盖类型与构成及其对CN值的影响[学位论文]. 南昌: 江西师范大学, 2014.]
|
| [19] |
Xiao B, Wang QH, Fan J et al. Application of the SCS-CN model to runoff estimation in a small watershed with high spatial heterogeneity. Pedosphere, 2011, 21(6): 738-749. DOI:10.1016/S1002-0160(11)60177-X |
| [20] |
Wang YJ, Lv HJ, Jiang T. Sub basin division and DEM resolution to the influence of SWAT runoff simulation research. Hydrological, 2008, 28(3): 22-25. [王艳君, 吕宏军, 姜彤. 子流域划分和DEM分辨率对SWAT径流模拟的影响研究. 水文, 2008, 28(3): 22-25. DOI:10.3969/j.issn.1000-0852.2008.03.005] |
| [21] |
Liu W. Based on DEM and remote sensing images extraction of linear deriving groundwater flow in Karst Basin[Dissertation]. Nanchang: Jiangxi Normal University, 2016. [刘维. 基于DEM和遥感影像提取的线性体推求岩溶流域地下水流向的研究[学位论文]. 南昌: 江西师范大学, 2016.]
|
| [22] |
Wei HB, Zhang ZP, Yang JP. The method of establish soil database of SWAT model. Water Resources and Hydropower Engineering, 2007, 38(6): 15-18. [魏怀斌, 张占庞, 杨金鹏. SWAT模型土壤数据库建立方法. 水利水电技术, 2007, 38(6): 15-18.] |
| [23] |
Zhu QC, Fan HD. The conversion of classification standard for soil particle composition. Chinese Journal of Soil Science, 1999, 30(2): 53-54. [朱秋潮, 范浩定. 土壤颗粒组成分级标准的换算. 土壤通报, 1999, 30(2): 53-54. DOI:10.3321/j.issn:0564-3945.1999.02.003] |
| [24] |
Qin FL. SWAT-Based simulation on non-point source pollution in the north watershed of Miyun Reservoir[Dissertation]. Beijing: Capital Normal University, 2006. [秦福来. 基于SWAT模型的非点源污染模拟研究——以密云水库北部流域为例[学位论文]. 北京: 首都师范大学, 2006.]
|
| [25] |
SWAT-CUP4: SWAT Calibration and Uncertainty Programs-A User Manual. Swiss Federal Institute of Aquatic Science and Technology, 2009.
|
| [26] |
Moriasi DN, Arnold JG, Vanliew MW et al. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Transactions of the ASABE, 2007, 50(3): 885-900. DOI:10.13031/2013.23153 |
| [27] |
Yang YH, Yan BX. The transport of nitrogen and phosphorus in rainfall-soil-runoff system. Journal of Soil and Water Conservation, 2010, 24(5): 27-30. [杨育红, 阎百兴. 降雨-土壤-径流系统中氮磷的迁移. 水土保持学报, 2010, 24(5): 27-30.] |
| [28] |
Hu N, Yuan H, Lan JC et al. Factors influencing the distribution of inorganic phosphorus fraction in different vegetation restoration areas in Karst rocky desertification area. Acta Ecologica Sinica, 2014, 34(24): 7393-7402. [胡宁, 袁红, 蓝家程等. 岩溶石漠化区不同植被恢复模式土壤无机磷形态特征及影响因素. 生态学报, 2014, 34(24): 7393-7402.] |
| [29] |
Kang ZQ, Xiong ZB, Li QY et al. The precipitation effect for water cycle way in Karst underground river basin. Earth and Environment, 2011, 39(1): 26-36. [康志强, 熊志斌, 李清艳等. 岩溶地下河流域水循环方式的降水效应. 地球与环境, 2011, 39(1): 26-36.] |
| [30] |
Wang SF. The research progress precipitation recharge Karst Aquifer System. Hydrology, 2014, 34(6): 1-8. [王树芳. 岩溶含水系统降水入渗补给研究进展. 水文, 2014, 34(6): 1-8.] |