地形和径潮条件变化对东江三角洲洪季连通性影响

2021-09-27余明辉王睿璞陈小齐刘长杰刘画眉

水科学进展 2021年5期

余明辉,王睿璞,,陈小齐,刘长杰,刘画眉

(1.武汉大学水资源与水电工程科学国家重点实验室,湖北武汉 430072；2.广东省水利水电科学研究院,广东广州 510635； 3.深圳市城市公共安全技术研究院,广东深圳 518000；4.哈尔滨工业大学(深圳),广东深圳 518000)

水系连通性反映河道干支流、湖泊及其他湿地等水流的连续性,具有自然属性和社会属性[1-2]。水系连通研究和连通性定量评价可为防洪、水生态保护、水环境改善及水资源开发与分配等提供基础支撑。近年来，国内外学者对连通性概念机制、定量评价及其在城市水系规划等方面应用开展了大量研究[3-6]。关于河流水系连通性的研究成果层出不穷,孟慧芳等[7]对比分析现有连通性评价方法,提出基于水流阻力和水文连通函数的评价方法可较好地评价平原河网区河流连通性；崔广柏等[8]认为水系结构连通性是水力连通性的基础,水系连通条件的提高可改善水环境；王延贵等[9]构建了包括河流边界、水流、泥沙及生态环境等在内的河道(湖库)连通评价体系。在筑闸建坝、河道采砂等人为干预下,部分地区的河网水系结构变化和功能退化,出现河流连通性衰退的局面[10]。此外,河流水系的水位和流量等存在明显的季节性差异和年际变化,河流水系连通性也会动态变化[2]。在河口三角洲水文连通性的研究中,Tejedor等[11]分析了三角洲河网脆弱性,量化了流量在河网区河段之间的分配过程与动态变化；Hiatt和Passalacqua[12]较为全面地考虑三角洲河道和岛屿,研究发现岛屿滩地在三角洲水文连通性中具有重要作用。

20世纪80年代以来,东江干流及三角洲采砂活动大规模兴起,导致河道演变特性由缓慢淤积为主转变为剧烈普遍下切,使采砂区河段河槽容积大幅度扩大,水深明显增加,河道横断面愈向窄深发展,而采砂区以上河段河床则表现出明显的溯源冲刷[13-14]。东江干流及三角洲地形剧烈变化引起河道水动力条件的改变:一方面,20世纪60年代前东江三角洲南、北水道枯水期分流比为45%和55%,但90年代后分流比逆转为60%和40%,且上游来流越小南支流分配越大[15]；另一方面,河床下切导致东江干流及东江三角洲同流量级下水位均有大幅度的下降,但1996年后樊屋站、石龙站、大盛站受上下游河道地形不均匀下切及河床纵比降减小的综合影响,3站水位不同程度的上升[16-17]。现有研究较多分析东江干流和东江三角洲河床演变以及水动力特征的改变,但聚焦到该区域水系形态和连通性研究的成果较少。

本文梳理1999—2009年间东江干流博罗—石龙河段及东江三角洲河网区地形变化,基于图论方法描述河流结构连通性和河流输运连通性,借助数值模拟手段研究径潮条件改变和地形变化对东江干流及三角洲河流连通性的综合影响,研究成果以期为改善东江水系河流连通性及东江防洪、水利工程规划和水生态环境的研究提供支撑。

1 东江干流及三角洲地形变化

东江发源于江西省寻乌县,干流博罗至石龙河段长约56 km,在广东省东莞市石龙汇入珠江三角洲,石龙以下河网水系称为东江三角洲(图1)。东江三角洲集水面积约1 380 km2,在石龙分出东江北干流和东江南支流,水系纵横交错,经大盛、麻涌、漳澎、泗盛注入狮子洋经虎门入海。

图1 东江干流及东江三角洲示意Fig.1 Map of the lower reaches of Dongjiang River and Dongjiang River Delta (DRD)

1972年以前东江干流及东江三角洲的采砂活动较少；20世纪80年代初改革开放以来,东江干流及三角洲地区需砂量大增,主要采砂区为东江干流、东江北干流和东莞水道等；20世纪90年代以后东江干流采砂区逐年上移到达合竹洲,这一阶段的采砂长期处于无序、无度状态,1997年东江流域主管机构介入严格控制采砂,实施河砂分段招标措施,其中东江干流博罗附近河段仍划为招标采砂的河段,加上偷采,采砂量仍较大。2004—2007年博罗以下河段采砂活动进一步严控,2009年以后采砂规模有序减小,直至停止[18](1972—2002年间主要采砂区及采砂量变化如表1)。

表1 东江干流及东江三角洲采砂统计[18-19]Table 1 Statistics of sand excavation in the lower reaches of Dongjiang River and DRD

1997—2008年间,东江典型河床断面博罗断面6.50 m高程以下河床平均下切3.62 m,年均下切0.33 m,深泓高程降低6.69 m,10 m水位下过水断面面积增加约59%。根据贾良文等[18]、Luo等[20]和胡晓张等[21]统计调查结论,1997—2009年间东江干流河床平均高程下切4.03 m,年均下切0.34 m；东江北干流石龙—新塘河段河床平均高程下切2.72 m,年均下切0.23 m；新塘—大盛河段河床平均高程下切0.07 m,但距大盛9 km的河段有冲有淤,河床平均高程抬高0.20 m；1997—2003间中堂水道河床平均高程增大1.35 m、厚街水道河床平均高程减少0.35 m。调查表明采砂主要采挖河床沙,多在河槽中进行，而近岸浅滩挖沙较少,且主要采砂区位于东江干流、东江北干流、东莞水道等河道[22]。以“1999年地形”为基础,由河床平均高程年均变化累加得到1999—2009年间主槽河床高程变化,结果如表2所示,另外，潢涌水道、大汾河、厚街水道河底平均高程变化按参考值选取,入海口处附近河段河床平均淤积0.20 m,以此推算得到2009年东江博罗以下河道断面地形(简称“2009年地形”),河道深泓和典型断面变化见图2、图3。

表2 1999—2009年东江河床平均高程变化推算值 mTable 2 Calculated values of average channel elevation change in the lower reaches of Dongjiang River and DRD from 1999 to 2009

图2 1999年和2009年东江河道深泓高程沿程变化Fig.2 Changes of thalweg elevation in the lower reaches of Dongjiang River and DRD from 1999 to 2009

图3 1999年和2009年东江干流及东江三角洲典型河床断面Fig.3 Typical cross sections in the lower reaches of Dongjiang River and DRD in 1999 and 2009

2 研究方法

2.1 河网区一维水动力模型

建立河网区一维河网水动力数学模型,模型的上游边界博罗站采用实测流量过程,增江麒麟咀站处的径流过程在新家铺作点源入汇处理,下游边界大虎站采用实测潮位过程(图4)。高程基准面采用珠江基面,建模地形采用1999年实测地形。

图4 模型边界及观察断面布置Fig.4 Model boundary and the location of observation sections

采用1998年6月21日至6月23日水文观测资料进行率定,得到东江河道糙率为0.023～0.031。采用同年6月25日至6月28日水文观测资料进行验证,图5为石龙站水位验证结果,均方根误差为0.09 m,未超过0.10 m[23],表明该模型可用来模拟东江干流及三角洲径潮动力过程。

图5 石龙站水位模拟值与观测值对比Fig.5 Comparison of simulated and measured water level at Shilong station

2.2 数值试验方案

相比于枯季,洪季河流径流较大,连通性较为显著,选取洪季作为本文研究时期。由于缺乏与建模采用的1999年地形年份和变化后2009年地形相匹配的洪水实测资料,选取2场与地形年份接近的实测洪水(发生于1998年6月、2008年6月,简称为“98·6”洪水、“08·6”洪水),洪潮过程如图6所示。表3为2场洪水来流及下游潮位特性,“08·6”洪水两站总来流量比“98·6”洪水增加约46%,但外海最大潮差减小0.04 m、平均潮差减小0.27 m。为探究不同径潮条件和地形条件下东江干流及三角洲河网区结构连通性和输运连通性的变化,计算过程中根据不同年份地形和不同径潮条件,组合成4种工况：工况1为“1999年地形”+“98·6”径潮过程实测资料；工况2为“1999年地形”+“08·6”径潮过程实测资料；工况3为“2009年地形”+“08·6”径潮过程实测资料；工况4为“2009年地形”+“98·6”径潮过程实测资料。

图6 洪水流量过程及大虎潮位Fig.6 Flood process and tidal water level in Dahu station

表3 东江流域“98·6”洪水、“08·6”洪水水文特性Table 3 Characteristics of flood “98·6”and “08·6”in the Dongjiang River basin

2.3 径潮动力描述指标

选用余水位、余水位比降及潮差描述东江三角洲径潮动力特性。余水位指1个太阴日25 h内的平均水位,即采用傅里叶变换方法对逐时水位序列进行以25 h为周期滤波后的低频成分；余水位比降为单位距离下余水位的增量,反映潮汐传播时遭受的非线性摩擦阻力,余水位比降越大对潮动力的抑制作用越显著；潮差指1个潮周期内,去除余水位影响后(即滤波后的高频成分)高潮位与低潮位的差值[24-25]。

2.4 连通性描述指标

(1)结构连通性。图论方法描述水系的形态格局时把河流的分汊汇流点概化为节点,河道概化为连接,水系概化为节点和连接组成的网状结构。

可选路径数(Number of Alternative Paths,Nap)描述河道平面形态上的成环程度。子河网为从上游至下游出口的沿程河道组成的输水通道。Nap为子河网中从上游进口至下游出口的路径数目,其值愈小,反映子河网中分汊、入汇节点较少,成环发育程度较差[11]。

连接共享系数(Linking Sharing Index,ILS)描述连接在子河网间的共享程度,计算式为

(1)

式中:k为子河网的序号；连接(ij)表示水流从节点i向节点j；Vk为子河网Sk的节点数目；bij为子河网Sk中连接(ij)在所有子河网中的重复次数[11]。

节点连接率(β)描述河网中每个节点与其他节点连接难易程度[10,26],计算式为

(2)

式中:L为河道数,0

(2)输运连通性。输运连通性包含水体的流动性和河流的连续性2个基本要素[1]。河道流量越大反映流动性越好、连续性越优、水量传输能力越强[27]。连通度因子计算式为

ωij=1-exp(-qij)

(3)

式中:ωij为河道输运连通度,值愈小代表节点i和节点j之间连通程度愈差；qij为河道时均流量,m3/s。将“98·6”洪水、“08·6”洪水进行同倍比缩小,缩小比取1960—2009年博罗站洪季多年平均流量1 053 m3/s。式(3)改写为

ωij=1-exp(-qij/1 053)

(4)

将连通度因子ωij代入邻接矩阵A=(aij)n×n得到加权邻接矩阵W,则判断矩阵

(5)

(6)

(7)

式中:D为河网连通矩阵；dij为节点i和节点j之间的平均连通度；nij为节点i和节点j之间各不同长度的连接的数目,nij≥1,nij∈N+；IDC为河网输运连通度。

3 结果及讨论

3.1 径潮动力特性变化分析

径潮动力特性的变化是输运连通性变化的内因。4种计算工况观察断面处(图4)的余水位和潮差指标如表4和图7所示。分别对比工况1和工况2、工况3和工况4,均可以发现:在同一地形条件下,当上游来流增大46%且大虎潮差减小0.27 m时,东江干流及三角洲河道余水位均增加而潮差均减小,且径流控制作用较强的三角洲上游河道余水位增幅大于潮流控制作用较强的三角洲下游河道,而下游河道潮差的降幅则大于上游河道。

表4 计算工况余水位和潮差变化Table 4 Changes of residual water level and tidal range of model cases

分别对比工况1与工况4、工况2与工况3,当上下游遭遇的径潮过程相同时,河道地形变化后,东江干流及三角洲河道余水位均降低,变化幅度沿程分布表现为位于径流控制区愈上游的河道余水位增幅愈大,而潮差变化与河道地形变化相关,大盛潮差减小0.22 m、泗盛潮差仅增大0.07 m,这与当地主槽河床高程平均淤积0.2 m的统计结果一致,除此之外河道潮差整体上增大,且潮流控制区愈下游的河道潮差增幅愈大。

工况1与工况3对比可探究地形变化和上下游不同径潮遭遇对东江干流及三角洲水动力的综合影响,由表4和图7可知,东江干流径潮特性受地形变化影响较大,而东江三角洲河网区受径潮过程的影响较大,表现为东江干流河道沿程余水位降低和潮差增大,东江北干流、东江南支流河道沿程余水位增大和潮差减小。

图7 余水位和潮差变化Fig.7 Change of residual water level and tidal range

3.2 结构连通性对地形和径潮条件变化的响应分析

东江干流及东江三角洲可划分为4个子河网,分别用符号S1、S2、S3、S4表示(图8)，顶点均为博罗,出口分别为大盛、麻涌、漳澎、泗盛。基于图论将东江干流及三角洲河网概化为有向图(见图9),其中节点23个,连接31条。东江干流及三角洲各河道连接重复次数bij值如图10所示。表5为结构连通性描述指标Nap、ILS和β计算值,子河网S4的Nap值最大,表明河网内有诸多河道分汊、入汇点,河网环形输水通道较多、成环发育程度较好,因而增加了水流入海的路径；子河网S1、S2的Nap均为1,表明其流路唯一,不存在环形输水通道,说明该子河网成环发育程度不高。由图10可见，bij=4的河道仅3条,占总河道数的9.4%；不存在bij=3的河道；bij=2的河道共10条,占总河道数的31.3%；bij=1的河道共19条,占59.3%,多为子河网S4的输水通道。由式(1)计算子河网连接共享系数值为0.28～0.55,其中，子河网S4的ILS值最小,说明S4子河网结构上共享性较低；子河网S1、S2、S3的ILS均高于0.50,表明共享性较高。由式(2)计算子河网节点连接率值变化范围为0.83～1.42,东江干流及三角洲的平均节点连接率β=1.39,表示平均每个节点有1～2个连接线,说明河网有较好的连接途径。计算结果表明，整体上东江干流及三角洲河网成环发育程度不高,但河道共享性较高,有较好的连通途径,总体结构连通性较好；不同径潮条件和地形条件下,结构连通性未受影响。

图8 东江干流及三角洲子河网示意Fig.8 Schematic of sub-network of the lower reaches of Dongjiang River and DRD

图9 东江干流及东江三角洲概化有向图Fig.9 Directed map of the lower reaches of Dongjiang River and DRD图10 连接重复次数Fig.10 Value of linking sharing frequency bij

表5 东江干流及三角洲结构连通性指标Table 5 Topologic connectivity metrics of the lower reaches of Dongjiang River and DRD

3.3 输运连通性对地形和径潮条件变化的响应分析

4种工况下输运连通度计算结果如图11和表6所示。分别对比工况1和工况2、工况3和工况4,其结果表明:在同一地形条件下,相比于上下游遭遇“98·6”径潮过程,在“08·6”径潮过程下河道输运连通度和河网输运连通度均增大,表明来流增大及外海潮动力减弱使东江干流博罗—石龙河段及东江三角洲河网内河道流量普遍增大,余水位比降增大反映潮波上溯的阻力增大,综合作用下东江干流及三角洲河道潮动力减弱,因此，来流较大、外海潮动力较弱的“08·6”径潮过程条件下河流输运连通性增强。分别对比工况1和工况4、工况2和工况3,其结果表明:在东江三角洲遭遇同一径潮过程时,地形变化对子河网输运连通度的影响无显著特征,但河网整体输运连通度降低,这表明东江干流博罗—石龙河段及三角洲地形不均匀下切引起河道流量重新分配,且余水位比降减小,反映潮波上溯的阻力减小,河道内潮动力增强、水量传输能力下降,说明地形下切使输运连通度减小、输运连通性减弱。对比工况1和工况3,工况3河道输运连通度变化范围为-0.231 3～0.272 9,其中输运连通度减小的河道共4条,占全部河道数的12.9%,子河网输运连通度增大范围为0.014 3～0.391 7,整体上河网东江干流及三角洲输运连通度增大0.017 6；且相比于地形下切使河网输运连通度减小了0.004 2～0.008 1,来流及外海潮动力带来的增大幅度在0.021 8～0.025 7,说明洪季径潮条件对输运连通性的影响程度较大,洪水流量越大则输运连通性越优。