广州珠江感潮河网水流泥沙数值模拟*

李 彤,李适宇

(中山大学环境科学与工程学院，广东 广州 510275)

1 前 言

广州珠江河网属于典型的潮汐水道，上游径流的水期差异以及下游河口的潮汐过程均能显著影响河网水力特征；同时，由于河网汊口众多、河道河涌交织成网，水流往往需几经分流、汇流，方能顺利下泄，这也导致珠江河网广州段水动力与物质输运过程具有较高的复杂性[1-3]。

在河网输沙及河床演变研究方面，前人业已展开了较为丰富的研究工作。乔彭年[4]对广州河道冲淤变化原因进行了分析，认为边界条件的永久性改变及水沙条件的差异性是导致广州珠江河网独特的河床演变特征的主要原因；罗宗业[5]则主要从联围筑闸、河道淤积及河道束窄对珠江洪水位壅高度的贡献进行了综合分析；江沛霖[6]分析了1955-1975年上游淡水来量、河口断面因素及河道缩窄对前后航道冲淤影响，并从中总结了径流及河岸边界变化条件下河床冲淤变化规律；何用等[3]采用二维潮流泥沙数学模型，对广州珠江河网中的南航道在洪、枯季条件下，底泥冲淤过程以及不同水文条件下河床的极限冲刷量进行了数值模拟。但总的来看，过去的研究多以长尺度历史地形资料对比分析为主[4-6]，借以数值模拟方法为预测手段的研究还相对较少，针对广州河网整体水沙输运过程及底泥冲淤过程的研究也不多见。

鉴于此，本文基于EFDC模型，构建广州珠江河网三维水动力与泥沙输运模型，对洪枯季水动力过程、悬沙输运过程以及河床冲淤过程进行模拟研究。探讨不同水期条件下，悬沙在广州河网水道中的输移特征及冲淤分布特征，为航道管理部门预测航道整治范围，同时也能为泥沙吸附态污染物的输运、累积等环境问题提供有益的参考。

2 研究区域水文泥沙特征

本文研究的广州珠江河网始于鸦岗断面以下白坭河，经白鹅潭汊口分别向东和向南分流至前航道和南航道；前者在鱼珠码头附近经黄埔断面汇入狮子洋，后者在牙髻沙汊口与北江分流水道——平洲水道汇合流入沥滘水道，并进一步分流，分别经官洲河、新造水道经黄埔汇入狮子洋，如图1。

根据1999年7月与2001年2月两次较大范围的水文与泥沙联合观测结果，广州珠江河网的水体悬移质主要以粘粒和粉粒为主[8]。其中，典型枯水期时(2001年2月，后简称“0102”)，河网悬沙的中数粒径一般在0.010～0.018 mm之间，而洪季时(1999年7月，后简称“9907”)在0.006～0.030 mm之间。河床质上看，枯季时的河床泥沙较细，主要由粘粒和粉粒组成，中数粒径与悬移质相近，在0.016～0.021 mm之间；洪季时，老鸦岗和沙洛围两个主要输沙断面的底泥泥沙组成较细，中数粒径分别为0.002和0.003 mm。但河网中的浮标厂粉粒比例有所增加，中数粒径为0.018 mm；黄沙断面沙粒组分较高，中数粒径达1.080 mm。参考以往研究[9-12]，本文采用EFDC模型模拟珠江河网的泥沙输运过程；此外，考虑到珠江河道推移质输沙率不及总输沙率的15%[13]，研究中暂没有考虑推移质输移过程，而主要以悬移质模拟为主。

图1 研究区域及模拟边界及河网地形、网格及主要验证站点分布

3 数学模型

采用美国EFDC(Environmental Fluid Dynamics Code)模型建立广州珠江河网水动力模型、泥沙模型。EFDC是由美国弗吉尼亚海洋科学研究所(Virginia Institute of Marine Science)John Hamrick教授开发，并由数家科研单位(Tetre Tech Inc.，Dynamic Solutions等)后期维护下发展起来的，可用于近岸海洋及陆地地表水水动力-水质模拟的三维水环境数学模型[14-15]。平面上EFDC采用曲线正交坐标，垂向采用Sigma无量纲化坐标，具有较好的边界拟合能力。水动力模型采用EFDC-Hydro模块进行模拟；泥沙模型采用EFDC-SED模块，对悬浮泥沙的对流扩散、沉降再悬浮过程进行模拟[16-17]。

3.1 控制方程

EFDC基于静水压假定和Bossinesq近似，对不可压缩、自由表面水体的连续、动量及状态方程进行求解。各控制方程如下[14]

u,v方向动量方程

∂t(mHu)+∂x(myHuu)+∂y(mxhvu)+∂z(mwu)-

(mf+v∂xmy-u∂ymx)Hv=-myH∂x(gζ+p)-

my(∂xh-z∂xH)∂zp+∂z(mH-1Av∂zu)+Qu

(1)

∂t(mHv)+∂x(myHuv)+∂y(mxhvv)+∂z(mwv)-

(mf+v∂xmy-u∂ymx)Hu=-mxH∂y(gζ+p)-

mx(∂yh-z∂yH)∂zp+∂z(mH-1Av∂zu)+Qv

(2)

连续方程

∂t(mξ)+∂x(myHu)+∂y(mxHv)+∂z(mw)=0

(3)

其中，u，v和w分别是x,y,z方向流速；mx，my，mz为坐标转换张量；ζ是水位，h为河床高程，H=h+ζ为总水深；p为压强，g为重力加速度；Av和Ab为垂向紊动黏性系数和扩散系数，采用Mellor和Yamada模型求解。Qu和Qv为u,v方向动量源汇项；f为科氏力系数。

3.2 泥沙输运方程

泥沙模型由悬浮泥沙的对流扩散、沉降再悬浮过程组成[15]

∂t(mxmyHS)+∂x(myHuS)+∂y(mxHvS)+

(4)

其中，KH和Kv分别为泥沙水平和垂向扩散系数；ws为泥沙沉降速度；S为水体含沙浓度。QE和Qt分别为悬沙源汇项。

“西湖山水还依旧……看到断桥桥未断，我寸肠断，一片深情付东流！”白衣女子一挥水袖，哀怨的歌声隐隐传来。

3.3 模型设置

对广东省航道勘察设计研究院测绘的1999年珠江河网地形图(1∶5 000)进行数字化，提取河床底部高程后，对地形进行插值、平滑；采用GEFDC(EFDC模型自带网格生成程序)对模拟区域进行网格划分，生成贴体曲线网格(如图1b所示)。计算域的水平网格为204行×250列，最小单元为46.7 m×22.6 m，最大为114.8 m×551.9 m；垂向水体按水深均匀分5个σ层。计算时间步长取值5 s。

基于珠江水利委员会1999年与2001年两次水文泥沙联测资料以及中山大学《广州、佛山跨市水污染综合整治方案》项目2001年1月和6月两次水文联测资料，选取2001年1月和2001年2月两个枯季时段以及1999年7月和2001年6月两个洪季时段的连续监测数据，用于模型相关的率定与验证。验证的指标有水位、流量、流速、流向和泥沙浓度5项，其中水文要素为每小时监测一次，泥沙浓度为两小时监测一次；流速、流向及水体悬沙数据为断面平均值。率定与验证案例的基础信息如表1所示。水动力模型边界采用“上游流量过程驱动，下游水位过程控制”的原则进行设置。模型在“冷启”计算后，保留稳定流场结果作为下次“热启”计算的初始场条件。

表1 水动力与泥沙模拟概况及验证项目

3.4 参数率定与验证结果

通过后 “0102”枯水期实测资料对模型糙率(EFDC模型中表现为摩擦高度)进行了率定，结果在在0.001～0.02 m之间，其中糙率值自上游河段向下游河段逐渐减小；上游河道在0.001 5～0.02 m之间，下游各河段在0.001 0～0.001 5 m之间；在水位模拟方面，验证的 7个站点的水位误差值都在10 cm以内；流量结果的相对误差较低，30%误差以下的模拟时段占总模拟时段的67.7%～87.5%；断面平均流速与流向与实测都基本吻合，说明水动力模型能够客观真实地反映河网内的水流运动状态，可以为后续的研究提供可靠的动力基础。代表断面的流量、水位及流向验证如图2-图4。

图2 1999年7月浮标厂断面水位(a)、流量(b)验证

图3 1999年7月黄沙断面水位(a)、流量(b)验证

图4 2001年1月4日至11日沙砖厂断面平均流速(a)与流向验证(b)

基于枯季实测泥沙过程，对泥沙模型中的沉降速率WSED0、临界沉积应力TAUD、临界再悬浮速率WRSP0以及临界冲刷应力TAUR进行率定。泥沙颗粒的沉降速率在斯托克斯静水沉速公式基础上率定得到。临界冲刷应力与临界沉积剪切应力经率得到的比值为1.6倍，与建议值1.2倍较为接近[18-19]，本文在各参数取值为如表2所示。泥沙模型对各站点模拟的断面平均泥沙浓度与实际值的相对误差在27.2%～43.5%之间，平均误差为35.5%。代表性断面悬沙验证如图5。

表2 泥沙模型主要参数率定结果

图5 1999年7月黄沙站(a)与浮标厂站(b)断面平均悬沙浓度验证

4 结果与讨论

4.1 广州珠江河网水流与输沙特征分析

广州珠江河网纵横，水道贯通，受径流与潮汐双重影响，水流运动十分复杂；加之历史上人为河道窄缩工程，使得河网岸线发生一定程度的形变，进一步影响了不同河段水流运动特点。其中，前航道形成的独特的“上窄下宽”的喇叭型结构，使得水流在前航道上段流速相对较大，而行至下端，尤其是猎德大桥以下时，河道逐渐放宽，水流减缓。以“9907”典型洪水期为例，上段海珠桥附近断面绝对流速平均值约为0.46 m/s，而在猎德大桥附近时，流速降至0.39 m/s；至下游“喇叭口”附近时，流速进一步降至约0.26 m/s。由于水流流速沿程逐渐降低，也势必对水体悬沙的输运与淤积过程产生影响。模拟发现，南航道涨潮流速大于落潮流速。以浮标厂断面为例，7月16日至20日之间，模拟的断面平均最大涨潮流速在0.54～0.74 m/s之间，落潮最大流速在0.22～0.34 m/s之间，与实测值基本一致[1]。模拟结果还发现，上溯水流还能长驱直入白坭水道，并通过落潮过程分流进至前航道。

在河网输沙方面，由于白坭水道上游修建了芦苞涌、西南涌水闸的修建，以往北江分流水沙通量均被大幅截断，北江分流至平洲水道的水流携沙是珠江河网广州段水体悬沙主要来源。从河网水期输沙量上看，“9907”洪季平洲水道经沙洛围断面输沙量达108.28万吨，而“0102”枯季仅为2.96万吨[3]。由此可见，北江向广州珠江河网洪季的侧向分流能够深刻影响到河网的水量以及泥沙浓度，把握住洪季输沙过程对于研究广州珠江河网泥沙输运与分布特征至关重要。洪季模拟结果表明，高浓度的悬浮泥沙主要输沙路径是经由平洲水道进入广州珠江河网，并沿沥滘水道、三枝香水道再进入后航道，并分别经新造水道、官洲水道向黄埔方向输运(如图6所示)。同时，平洲水道注入河网区的高浓度悬沙在潮汐顶托作用下，乘涨潮过程沿南航道上溯。大潮时，能够穿过白鹅潭汊点进入西航道。当发生落潮时，高浓度泥沙水团又经由白鹅潭分流进入前航道，形成间歇性的高浓度含沙水团，这与朱继伟等[7]在对广州内港航道的水体悬沙分布分析中的发现的珠江南航道—前航道泥沙输运过程是一致的。这一特点将使得北江泥沙能够在潮汐过程中分布到前后航道的各个河段水体，并在落淤过程沉积在河段底泥，这也势必影响到泥沙表面的各类吸附态污染物在河网中的输移和分布。

图6 广州珠江河网洪季8日平均悬沙浓度场

4.2 不同河段冲淤分布初步分析

为了进一步分析广州珠江河网河床演变规律与趋势，对珠江河网在典型洪、枯季河床冲淤过程进行分段统计。参考朱继伟等[7]对广州内河航道1964-2006年近50年的河床演变过程的河网统计中的划分方法，本文同样将研究河网区划分为白坭河、白沙河、沙贝海、前航道、南河道、沥滘水道、官洲河以及新造水道九河段(如图7)，并以“9907”和“0102”典型洪、枯季河床冲淤特征及冲淤量进行模拟，分别计算各时段的冲淤量，并按等比倍数叠加得到全年河床形变量。

图7 广州珠江河网冲淤分析河道分区示意图

图8 1999年7月(a)与2001年2月(b)洪、枯季广州珠江河床冲刷与淤积量,m

模拟的珠江河网冲淤结果如图8所示。以水体悬沙模拟结果稳定后的168 h(共模拟192 h)的河床演变量为基础，统计7 d内各网格的河床底泥容积变化量，通过加和各单元冲淤体积得到河段总冲淤量后，最后除以各河段总面积得到各河段平均冲淤厚度。从洪季模拟的河床冲淤态势上看，广州珠江河网主要以淤积为主。西航道的水口水道和白坭水道的鸦岗与硬颈海段、前航道“喇叭口”末段、南航道、沥滘水道、三枝香以及官洲河下段、新造水道上段均以淤积趋势为主。其中，以南航道淤积速率最大，达到18 cm/a。如前所述，受潮汐顶托作用，西北江经平洲水道分流分沙与涨潮水沙叠加后涌入南航道，涨潮流速大，使得高浓度泥沙水流能长驱上溯，将南航道河段平均泥沙浓度维持在较高的水平(如图6)；憩潮时水体悬沙大量落淤，形成较高的淤量。前航道上段及中段河道由于河段相对窄缩，水流较快，不易发生淤积而主要表现为冲刷，末端河口表现为淤积特点。枯季时，河网主要以冲刷过程为主，白坭水道、白沙河、沙贝海、南航道—新造水道以及前航道均表现出明显的冲深特点。其中以新造水道、白沙河为最，年冲刷率可达6 cm/a和7 cm/a。

整体上看，广州珠江河网(特指本节纳入统计的九条河段)表现为洪季淤积枯季冲刷的特点，洪季淤积速率约为5.7 cm/a，枯季冲刷速率约为1.7 cm/a，全年呈缓慢淤积趋势，速率约为2.0 cm/a(全河段淤积量除以河段总面积)。

图9 珠江河网年冲淤厚度预测值与实测值比较

模拟的9条河段河床的冲淤厚度(正数表示淤积，负数表示冲刷)及与前人研究的实际冲淤量对比如图9所示。可见，模型计算得到的珠江各河段冲淤厚度与前人的实测结果较为接近。其中，在河网主要淤积河段——南航道平均淤积速率模拟值约为8.1 cm/a，与实测淤积速率8.9 cm/a较为吻合，与前人模拟结果(8 cm/a)也十分接近[3]。沥滘水道淤积速率约为8.5 cm/a，与实测结果(9.3 cm/a)也比较接近。然而，白沙河和沙贝海两个河段模拟的冲淤状态与实测结果正好相反，并且白鹅潭和新造水道的模拟结果与实测结果也存在较大的误差。经分析，这可能有以下几点原因：首先限于研究资料的不足，本文尚不能对河网全年水流输沙过程进行模拟，而是采用两个不同年份洪、枯季节的冲淤过程加权得到，因此无法完整地反映河床长时间尺度上的形变过程，而且引用的实测值也是前人对1999-2006年(或2005年)河网平均冲淤值，因此一定程度上影响到计算结果的精确性。其次，沙贝海虽表现为淤积趋势，但1999-2006年，沙贝海河道的平均河宽减少了近10 m，河道过水面积减少约36.1 m2，因此从断面过水面积上看，该河段体现为淤积，但其水深却有所增大，平均水深从4.49 m增加至4.54 m，年均冲深约1 cm[7]，这与本文模拟的趋势是一致的。此外，新造河道的模拟冲刷量小于实测值，这可能历史上与人为的采砂活动和航道疏浚工程有关[3]。如1998年洪、枯季，珠江干流分别就有8艘采砂船从事采砂活动(其中有7条船是在新造水道作业)，年采砂量达161万吨[13]。这对新造水道的河道地形势必造成一定的影响，也可能是模拟值小于实际值的原因之一。从整体的冲淤模拟结果上看，本文使用的泥沙模型基本上把握住了珠江河网西航道、南航道及沥滘水道的淤积的主要特点，计算得到冲淤量与实际值总体上比较接近，能初步反映出珠江河网河床冲刷淤积的变化过程。

5 结 论

1)基于EFDC模型，构建了广州珠江河网三维水动力与泥沙数学模型。通过1999年7月和2001年2月典型洪、枯季水流泥沙实例验证，表明模型能较好地模拟和反映河网感潮水流及输沙特点，计算结果与实测资料较为吻合。

2)洪季北江经平洲水道分流分沙过程是广州珠江河网主要泥沙通量来源，主要的输沙路径集中在后航道沥滘水道—新造水道、官洲河—黄埔线。而受潮汐能量顶托作用，高浓度悬沙能经由南航道上溯至白坭水道，并在落潮时通过白鹅潭分流进入前航道，形成间歇性高浓度泥沙团，使得北江分沙得以在广州珠江河网中广泛分布。

3)通过预测主要河段的冲淤趋势及厚度，结果显示：21世纪初期，广州珠江河网总体呈缓慢淤积态势，年淤积速率约为2.0 cm/a；且在不同水期主要表现为“枯冲洪淤”的特点。其中，南航道、沥滘水道是主要的淤积河段，速率约在8～9 cm/a；而沙贝海、新造水道是主要冲刷河段，但人为窄缩河道、航道疏浚及人为采砂活动对河道形变趋势都能产生一定影响。

4)由于广州珠江河网资料相对较少，模拟时段较短，本文仅对河网底泥冲淤分布进行了初步的研究。在今后的研究中仍需进一步收集长时间尺度的水文泥沙资料，对河网河床演变过程进行更深入的模拟分析。

参考文献：

[1]包芸，姚冬静.珠江网河二维水动力数值计算研究[J].水动力学研究与进展，2005，20(增刊)：881-886.

[2]包芸，刘欢.风暴潮极值状态下珠江河网白坭水道内的不稳定流动[J].大连大学学报，2007，28(3)：22-27.

[3]何用，徐峰俊.珠江洲头咀河段过江隧道工程极限冲刷模拟研究[J].人民珠江，2007(5)：22-25.

[4]乔彭年，珠江三角洲河道冲淤特征的初步分析(下)[J].人民珠江，1984(5)：9-17.

[5]罗宗业，珠江洪水上涨——广州水侵街[J].海岸工程，1994，13(4)：36-43.

[6]江沛霖.广州附近河道的冲淤因素[J].人民珠江，1983(1)：36-41.

[7]朱继伟，武亚菊，刘俊勇,等.广州内港航道近50年来河床演变浅析[J].人民珠江,2008(3):15-17.

[8]王兴奎，邵学军，李丹勋.河流动力学基础[M].北京: 中国水利水电出版社,2002: 26.

[9]HydroQual Inc.A primer for ECOMSED,User's mannual[R].Mahwah,N.J.07430.

[10]WRIGHT S A,BERECIARTUA P.Case study: sediment transport in proposed geomorphic channel for napa river[J].Journal of Hydrulic Engineering,2001,901-910.

[11]ZIEGLER K,NISBET B S.Long-term simulation of fine-grained sediment transport in large reservoir[J].Journal of Hydraulic Engineering,1995,121(11):773-780.

[12]SEBNEM Seker-Elci.Modeling of hydrodynamic circulation and cohesive sediment transport and prediction of shoreline erosion in Hartwell lake,SC/GA[D].Georgia institute of Technology，2004.

[13]罗宪林,杨清书，贾良文,等.珠江三角洲网河河床演变[M].广州： 中山大学出版社,2002：102.

[14]HAMRICK J M.A three dimensional environmental fluid dynamics computer code: theoretical and computational aspects [R].Virginia Institute of Marine Science,1992.

[15]TETRA Tech Inc.Theoretical and computational aspects of sediment and contaminant transport in the EFDC model [Z].Fairfax,Virginia 22030，2002.

[16]齐珺， 杨志峰，熊明,等.长江水系武汉段水动力过程三维数值模拟[J].动力学研究与进展，2008，23(2)：212-219.

[17]陈景秋，赵万星，季振刚.重庆两江汇流水动力模型[J].水动力学研究与进展，2005，20(增刊)：829-835.

[18]JI Z G.Hydrodynamics and water quality: modeling rivers,lakes and estuaries [M].Hoboken,New Jersey: John Wiley & Sons Inc.,2007:197.

[19]JIN K R,JI Z G.Case study: modeling of sediment transport and wind-wave impact in lake Okeechobee [J].Journal of Hydraulic Engineering,2004,130(11):1055-1067.

[20]HWANG K N,MEHTTA A J.Fine sediment erodibility in Lake Okechobee[R].Rep.No.UFL/COEL-89/019,Univ of Florida,1989:140.