山丘区河流滩地边缘内凹形态对溶质运移及水质恶化风险的影响

周之悦，夏继红，蔡旺炜，叶继兵，刘 瀚，窦传彬

(河海大学农业科学与工程学院，南京 210098)

0 引 言

滩地作为水域与陆域的过渡区域，是河道中水文、环境和生态过程相互作用的敏感区，对洪水调蓄、水质净化、生物栖息等具有重要价值[1,2]。我国广大的山丘区中分布着大量的中小型河流，其中不乏形态多样的滩地。然而人类活动、异常水文过程、地质活动等因素易造成滩地形态及格局的破坏，使得滩地功能退化，影响河流生态系统健康[3,4]。Mejia等[5]通过开发基于瞬时响应函数和流域地貌关系的分析框架，探讨了滩地格局塑造的河道断面形态对流域水文过程的影响；蔡建楠[6]等通过构建含河道断面形态、河岸带宽度等指标在内的城市河流形态评价体系，发现河流形态评分较高的河段水体自净能力也越强。伊紫函等[7]对滩地平面几何形态进行了分类以研究滩地整体形态演变；王金平等[8]应用灰色关联方法分析滩地景观格局时空演化特征，发现近年来滩地破碎化特征明显且常出现内凹向滩地腹地的边缘形态。然而目前关注滩地形态的研究多涉及流域尺度上的地貌格局探讨，对中、小尺度相关问题的定量研究较少[9]，尤其缺乏对滩地边缘形态的研究，如，滩地边缘内凹形态变化对水体溶质运移过程及水质恶化风险影响机理等问题尚待深入研究。本文以龙游县灵山港为典型，通过构建定量表征滩地边缘内凹的边缘形态指数(Ie)，探讨滩地内凹程度等级划分方法，深入研究不同内凹程度对溶质运移的影响及水质恶化的风险，以期为山丘区河流滩地修复提供数据支撑和理论依据。

1 材料与方法

1.1 研究区概况

研究区域(E119°07′～119°12′，N28°49′～29°03′)为浙江省西部的灵山港，是钱塘江上游衢江右岸的一条重要支流。流域总面积 726. 9 km2，属亚热带季风气候区，年平均降水量为1 815 mm，年内分布不均，主要集中在3-9月份。多年平均径流量20.8 m3/s，径流系数为0.62。流域内山地土壤以黄红壤亚类、侵蚀性红壤、黄壤为主，母质为凝灰岩、片麻岩、板岩等。河床及滩地以堆积和侵蚀作用为主，堆积物一般厚为3～8 m。灵山港主流长43.79 km，河道平均比降为0.245%，沿线设有步坑口、龙游水文站。近40年来，灵山港内的滩地呈现逐渐破碎和萎缩的变化特征，下游河段(姜席堰-兰石大桥 段)最为明显[8]。该河段(图1)总长6.1 km，平均比降为0.181%，河道宽度为50～160 m不等，共有18处滩地，其分布格局复杂多样，流速分布空间差异较大。

图1 研究区滩地分布格局及计算网格剖分示意图Fig.1 The patterns of floodplains in the study area and sketch map of computational mesh

1.2 数据采集

分别于2017年4月和10月、2018年4月和10月对灵山港下游段开展了实地调查，采集河道地形、滩地边缘形态、流量与水位、溶质(总磷TP)浓度等基础数据。方法：①河道地形数据由当地水利部门提供；②滩地边缘形态参数采用手持GPS记录仪(MAGELLAN，美国)实地踏勘而得；③流量与水位数据主要源于水文监测站，另于校核断面(图1(a))处设置了临时监测站点，记录逐日水位；④水质指标TP通过现场采集水样带回实验室参照《流动注射-钼酸铵分光光度法 HJ 671-2013》测定。

1.3 二维水动力-水质模型构建

基于Mike 21构建灵山港下游段水动力-水质耦合模型。水动力学控制方程为基于Boussinesq假定和静水压假设的二维不可压雷诺平均N-S方程组(即浅水方程)：

连续性方程：

(1)

动量方程：

(2)

(3)

基于对流-扩散方程耦合水动力模块模拟溶质运移过程，水质模块控制方程如下：

(4)

式中：c为标量浓度；Dx和Dy分别为x和y方向上的扩散系数。

通过对计算区域进行网格划分，采用基于单元中心的有限体积法对物理控制方程进行数值求解。初始条件、边界条件等参数通过具体模块设置输入，其中水动力模块为基础核心模块。由于山丘区河流河床地貌形态变化复杂，水流方向多变，滩地的出露形态随河流水位的动态变化而改变。遂依据2016年河道及滩地实测地形图，选择三角形非结构化网格对计算区域进行剖分，共生成网格330 377个，允许的最大网格面积为5 m2，网格最小角度为32°。网格分布特征如图1(b)所示。

进行模拟计算时以常水位流量(Q=15 m3/s)为上游边界条件，以相应水位(h=47 m)为下游边界条件。因研究河段较短，不考虑降雨、蒸发及物质降解。参照实际情况，设置河道全水域TP背景值为0.023 8 mg/L。分析模拟结果时，当TP浓度<0.040 mg/L即认为恢复至背景值水平。待水动力条件稳定后，上游边界持续2 h输入TP浓度为8 mg/L。以溶质投放时刻作为模拟计算的初始时刻，计算步长为30 s，共计3 600步(30 h)。

1.4 滩地边缘内凹形态定量表征

借鉴河道蜿蜒度的概念[10]，本文将滩地边缘曲线长度(S)与近滩水域开口宽度(L，与主河道水流方向平行)的比值定义为弯曲度指数Ie，即Ie=S/L，用以定量表征滩地边缘内凹形态。Ie值越大表明滩地边缘越弯曲，内凹形态向滩地腹地延伸的程度越高。

1.5 水质恶化风险的表征方法

溶质运移是指溶解在水中的物质输运和转移的空间和时间模式[11]。当某处水域的营养盐或污染物的浓度在较长时间内偏高时，则该水域的水质恶化风险较高。在本研究中，①采用近滩水域TP与主河道TP均值的比值(即TP的相对浓度R)来反映溶质在近滩水域的聚积程度，R越高水质恶化风险越高；②采用近滩水域溶质浓度等值线图反映TP浓度的空间梯度，等值线越密集表明水质恶化风险的空间变异性越剧烈，而等值线图分布的空间范围越大则说明内凹形态影响的水域面积越大；③根据近滩水域形心所在断面的溶质平均浓度的时间变化曲线，及该浓度与相应主河道断面的溶质平均浓度差值ΔC的时间变化曲线，来反映溶质在近滩水域中的累积与削减过程。

2 结果与分析

2.1 滩地边缘内凹形态分布特征

野外调查结果显示，共有14处滩地边缘内凹形态(编号为F1～F14)分布于灵山港下游段的11个滩地上。由表1可知，这些内凹形态的边缘曲线长度S、近滩水域开口宽度L(平行于河道主流方向)、近滩水域面积A等特征参数的分布范围分别是43.1～191.0 m、16.8～59.6 m和206～1 460 m2，内凹形态边缘弯曲度Ie的分布范围则为1.51～6.70。根据下游段滩地Ie的分布频率，将滩地边缘内凹形态划分为低弯(Ie<2)、中弯(2 ≤Ie< 4)和高弯(Ie≥ 4)。则F1～F6为低弯，F7～F9为中弯，F10～F14为高弯。

表1 灵山港滩地边缘内凹形态特征参数Tab.1 Floodplain edge morphological characteristics in the downstream of Lingshan River

注：S为滩地内凹形态边缘曲线长度；L为内凹形态开口宽度；A为内凹形态所围近滩水域面积。

2.2 模型校核及可靠性分析

基于2016年的流量监测数据进行模型参数率定，确定主河道曼宁系数为28～30，滩地曼宁系数为17～25。基于水文监测数据得校核断面2017年逐日水位的模拟与实测值吻合良好，Nash-Sutcliffe效率系数(NSE)为0.98，接近于1，模型质量较好，总体相对误差(RE)为6.29%(见图2)。表明所构建的水动力数值模型具有良好的稳定性及可靠性，适用于模拟溶质运移过程。

图2 校核断面2017年水位验证曲线Fig. 2 Validation of daily surface water level in 2017

2.3 内凹形态对溶质浓度空间分布的影响

选取位于研究区中部的内凹形态F11作为典型分析滩地内凹形态对近滩水域溶质浓度空间分布的影响。通过调整局部边缘形态使F11的Ie值从4.86(高弯)逐渐下降为2.43(中弯)和1.82(低弯)，形态特征如图3所示。由近滩水域溶质浓度等值线图(图3)可知，在近滩水域中，TP浓度与主河道存在明显差异的水域范围随着Ie的降低而减小，同时等值线的数量也随之减少。这一结果表明，滩地边缘内凹形态的弯曲度等级越高，其近滩水域中存在水质恶化风险的水域面积将会越大。

图3 不同边缘形态下近滩水域溶质浓度等值线图Fig.3 Solute concentration map of nearshore water

2.4 内凹形态对溶质浓度运移时间的影响

模型输出不同形态条件下近滩水域溶质浓度及其与主河道浓度差值的时间变化曲线。由图4(a)知，当Ie等级从高弯(M0)变为低弯(M2)后，近滩水域溶质浓度从峰值下降至背景水平(<0.040 mg/L)所需的时间从超过22 h缩短至11.8 h。而由图4(b)知，溶质浓度在近滩水域高于主河道所持续的时间从超过23 h(M0-高弯)下降至9.8 h(M2-低弯)，且势差峰值明显降低。上述结果表明，滩地边缘形态的变化对近滩水域溶质变化历程具有显著的影响。表现为随着Ie下降，近滩水域溶质浓度的变化过程与主河道的差异缩小。

图4 近滩水域溶质浓度变化曲线及 近滩水域-主河道浓度势差变化曲线Fig.4 The temporal curve of solute concentration in nearshore water and the concentration difference between nearshore water and the main channel

2.5 内凹形态对近滩水域水质恶化风险影响

野外调查结果显示14处滩地边缘内凹形态的近滩水域TP浓度(0.152～1.825 mg/L)与对应主河道断面的TP浓度(0.023～0.190 mg/L)差异较大，TP相对浓度R最大可达10倍。由图5可知，模拟结果(当主河道TP浓度恢复至背景水平时)与实测结果均显示，随着Ie从低弯升至高弯，R呈极显著(P<0.01，DUNCAN)升高。这一结果表明随着滩地边缘内凹形态的弯曲度等级升高，近滩水域水质恶化风险也随之增高。另外，相同弯曲度等级下，R的实测值均略低于模拟值。如高弯等级下，R的模拟值和实测值的均值分别为7.01±1.18和5.53±2.71。产生这一差别的主要原因是模拟计算时不考虑溶质的降解和植物吸收。

图5 不同弯曲度等级下近滩水域TP相对浓度RFig.5 Relative TP concentration under different Ie levels

3 讨 论

上述结果显示滩地边缘内凹形态的弯曲度指数Ie越大则近滩水域内溶质的相对浓度越高，溶质浓度下降速度越慢。根据余根听等[12]所开展的室内物理模型试验结果，正弦型蜿蜒岸的背水坡振幅越大(等价于Ie越大)，溶质在潜流带中的停留时间越长。因而，在近滩水域内，地表水与潜流可能会对溶质运移过程发生协同作用，进一步延长溶质在水体中的滞留时间，从而引发水质恶化风险，影响滩地生态系统的稳定性与功能性。

溶质动力学与水流的物理运动紧密相关。对于滩地错落分布的山丘区河流来说，对流作用对溶质运移的驱动影响要大大强于溶质势差主导的扩散作用。在暂态存储理论中，溶质被认为能够暂时存储在河流中的缓流或静水区，进而能够延长溶质向下游水体的释放时间[13]。而增强主流区与溶质暂态存储区之间的水流交换作用对河流系统健康至关重要[14]。从这个意义上说，对F11的形态改造本质上是通过控制边缘形态来控制主流区与溶质暂态存储区的对流关系，进而影响了暂态存储区对溶质的存储能力。

自然河流系统中，滩地周边广泛存在的缓水区、回水区及静水区能够为水生生物提供多样化的生境条件，尤其为鱼类在山丘区河流中提供了重要的休憩、产卵及索饵的场所[15]。在对滩地边缘形态进行优化改造时，一方面应当避免营养物质或污染物等溶质在近滩水域中停留时间过长而发生水质灾害，另一方面还应满足生境多样化的需求。根据本文研究结果，Ie为中弯(2 ≤Ie< 4)时，近滩水域不会出现过于明显的溶质积聚现象，水质恶化风险不高。实际上中弯条件还能提供一定面积比例的缓流区，能为水生生物提供良好的生境。对于高弯形态，通过改造局部形态就有可能使得其弯曲度从“高弯”下降为“中弯”。这类基于Ie指数进行设计的“微改造”既能避免对原有滩地整体形态及生态功能的强烈扰动，又能减少工程量。

此外，即便Ie处于相同等级，其近滩水域面积的大小可能差异较大(如F13和F14，见表1)。在实际情况中，当近滩水域的面积较大时，其TP相对浓度较高，反之，TP相对浓度较低(如F13的TP相对浓度R的实测值仅为1.84)。在后续研究中可将近滩水域面积大小与弯曲度指数Ie结合起来对滩地边缘内凹形态进行更精细的划分。

4 结 论

论文在野外调查的基础上，建立了滩地边缘内凹形态指数计算方法及等级划分方法，运用数值模拟探究不同弯曲等级对溶质运移及水质恶化的影响，为滩地治理和水环境改善提供参考。

(1)山丘区滩地边缘内凹形态可根据边缘弯曲度指数Ie划分为低弯(Ie< 2)、中弯(2 ≤Ie< 4)和高弯(Ie≥ 4)。

(2)滩地边缘内凹形态的变化能显著影响回水区的溶质分布和运移过程。Ie的变化对近滩水域溶质浓度的空间分布与时间变化过程均具有显著的影响。随着Ie降低，近滩水域溶质浓度的聚积程度、空间异质性、和削减耗时均明显降低。

(3)当Ie为低弯时，内凹形态所围近滩水域几乎不存在水质恶化风险；当Ie为中弯时，近滩水域水质恶化风险较低；当Ie为高弯时，水质恶化风险则较高。但是根据维持生境多样性的需求，不宜将内凹形态降至低弯。

□