廣州珠江感潮河網(wǎng)水流泥沙數(shù)值模擬*

李 彤,李適宇

(中山大學(xué)環(huán)境科學(xué)與工程學(xué)院，廣東 廣州 510275)

1 前 言

廣州珠江河網(wǎng)屬于典型的潮汐水道，上游徑流的水期差異以及下游河口的潮汐過程均能顯著影響河網(wǎng)水力特征；同時，由于河網(wǎng)汊口眾多、河道河涌交織成網(wǎng)，水流往往需幾經(jīng)分流、匯流，方能順利下泄，這也導(dǎo)致珠江河網(wǎng)廣州段水動力與物質(zhì)輸運過程具有較高的復(fù)雜性[1-3]。

在河網(wǎng)輸沙及河床演變研究方面，前人業(yè)已展開了較為豐富的研究工作。喬彭年[4]對廣州河道沖淤變化原因進行了分析，認為邊界條件的永久性改變及水沙條件的差異性是導(dǎo)致廣州珠江河網(wǎng)獨特的河床演變特征的主要原因；羅宗業(yè)[5]則主要從聯(lián)圍筑閘、河道淤積及河道束窄對珠江洪水位壅高度的貢獻進行了綜合分析；江沛霖[6]分析了1955-1975年上游淡水來量、河口斷面因素及河道縮窄對前后航道沖淤影響，并從中總結(jié)了徑流及河岸邊界變化條件下河床沖淤變化規(guī)律；何用等[3]采用二維潮流泥沙數(shù)學(xué)模型，對廣州珠江河網(wǎng)中的南航道在洪、枯季條件下，底泥沖淤過程以及不同水文條件下河床的極限沖刷量進行了數(shù)值模擬。但總的來看，過去的研究多以長尺度歷史地形資料對比分析為主[4-6]，借以數(shù)值模擬方法為預(yù)測手段的研究還相對較少，針對廣州河網(wǎng)整體水沙輸運過程及底泥沖淤過程的研究也不多見。

鑒于此，本文基于EFDC模型，構(gòu)建廣州珠江河網(wǎng)三維水動力與泥沙輸運模型，對洪枯季水動力過程、懸沙輸運過程以及河床沖淤過程進行模擬研究。探討不同水期條件下，懸沙在廣州河網(wǎng)水道中的輸移特征及沖淤分布特征，為航道管理部門預(yù)測航道整治范圍，同時也能為泥沙吸附態(tài)污染物的輸運、累積等環(huán)境問題提供有益的參考。

2 研究區(qū)域水文泥沙特征

本文研究的廣州珠江河網(wǎng)始于鴉崗斷面以下白坭河，經(jīng)白鵝潭汊口分別向東和向南分流至前航道和南航道；前者在魚珠碼頭附近經(jīng)黃埔斷面匯入獅子洋，后者在牙髻沙汊口與北江分流水道——平洲水道匯合流入瀝滘水道，并進一步分流，分別經(jīng)官洲河、新造水道經(jīng)黃埔匯入獅子洋，如圖1。

根據(jù)1999年7月與2001年2月兩次較大范圍的水文與泥沙聯(lián)合觀測結(jié)果，廣州珠江河網(wǎng)的水體懸移質(zhì)主要以粘粒和粉粒為主[8]。其中，典型枯水期時(2001年2月，后簡稱“0102”)，河網(wǎng)懸沙的中數(shù)粒徑一般在0.010～0.018 mm之間，而洪季時(1999年7月，后簡稱“9907”)在0.006～0.030 mm之間。河床質(zhì)上看，枯季時的河床泥沙較細，主要由粘粒和粉粒組成，中數(shù)粒徑與懸移質(zhì)相近，在0.016～0.021 mm之間；洪季時，老鴉崗和沙洛圍兩個主要輸沙斷面的底泥泥沙組成較細，中數(shù)粒徑分別為0.002和0.003 mm。但河網(wǎng)中的浮標廠粉粒比例有所增加，中數(shù)粒徑為0.018 mm；黃沙斷面沙粒組分較高，中數(shù)粒徑達1.080 mm。參考以往研究[9-12]，本文采用EFDC模型模擬珠江河網(wǎng)的泥沙輸運過程；此外，考慮到珠江河道推移質(zhì)輸沙率不及總輸沙率的15%[13]，研究中暫沒有考慮推移質(zhì)輸移過程，而主要以懸移質(zhì)模擬為主。

圖1 研究區(qū)域及模擬邊界及河網(wǎng)地形、網(wǎng)格及主要驗證站點分布

3 數(shù)學(xué)模型

采用美國EFDC(Environmental Fluid Dynamics Code)模型建立廣州珠江河網(wǎng)水動力模型、泥沙模型。EFDC是由美國弗吉尼亞海洋科學(xué)研究所(Virginia Institute of Marine Science)John Hamrick教授開發(fā)，并由數(shù)家科研單位(Tetre Tech Inc.，Dynamic Solutions等)后期維護下發(fā)展起來的，可用于近岸海洋及陸地地表水水動力-水質(zhì)模擬的三維水環(huán)境數(shù)學(xué)模型[14-15]。平面上EFDC采用曲線正交坐標，垂向采用Sigma無量綱化坐標，具有較好的邊界擬合能力。水動力模型采用EFDC-Hydro模塊進行模擬；泥沙模型采用EFDC-SED模塊，對懸浮泥沙的對流擴散、沉降再懸浮過程進行模擬[16-17]。

3.1 控制方程

EFDC基于靜水壓假定和Bossinesq近似，對不可壓縮、自由表面水體的連續(xù)、動量及狀態(tài)方程進行求解。各控制方程如下[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)

連續(xù)方程

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

(3)

其中，u，v和w分別是x,y,z方向流速；mx，my，mz為坐標轉(zhuǎn)換張量；ζ是水位，h為河床高程，H=h+ζ為總水深；p為壓強，g為重力加速度；Av和Ab為垂向紊動黏性系數(shù)和擴散系數(shù)，采用Mellor和Yamada模型求解。Qu和Qv為u,v方向動量源匯項；f為科氏力系數(shù)。

3.2 泥沙輸運方程

泥沙模型由懸浮泥沙的對流擴散、沉降再懸浮過程組成[15]

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

(4)

其中，KH和Kv分別為泥沙水平和垂向擴散系數(shù)；ws為泥沙沉降速度；S為水體含沙濃度。QE和Qt分別為懸沙源匯項。

“西湖山水還依舊……看到斷橋橋未斷，我寸腸斷，一片深情付東流！”白衣女子一揮水袖，哀怨的歌聲隱隱傳來。

3.3 模型設(shè)置

對廣東省航道勘察設(shè)計研究院測繪的1999年珠江河網(wǎng)地形圖(1∶5 000)進行數(shù)字化，提取河床底部高程后，對地形進行插值、平滑；采用GEFDC(EFDC模型自帶網(wǎng)格生成程序)對模擬區(qū)域進行網(wǎng)格劃分，生成貼體曲線網(wǎng)格(如圖1b所示)。計算域的水平網(wǎng)格為204行×250列，最小單元為46.7 m×22.6 m，最大為114.8 m×551.9 m；垂向水體按水深均勻分5個σ層。計算時間步長取值5 s。

基于珠江水利委員會1999年與2001年兩次水文泥沙聯(lián)測資料以及中山大學(xué)《廣州、佛山跨市水污染綜合整治方案》項目2001年1月和6月兩次水文聯(lián)測資料，選取2001年1月和2001年2月兩個枯季時段以及1999年7月和2001年6月兩個洪季時段的連續(xù)監(jiān)測數(shù)據(jù)，用于模型相關(guān)的率定與驗證。驗證的指標有水位、流量、流速、流向和泥沙濃度5項，其中水文要素為每小時監(jiān)測一次，泥沙濃度為兩小時監(jiān)測一次；流速、流向及水體懸沙數(shù)據(jù)為斷面平均值。率定與驗證案例的基礎(chǔ)信息如表1所示。水動力模型邊界采用“上游流量過程驅(qū)動，下游水位過程控制”的原則進行設(shè)置。模型在“冷啟”計算后，保留穩(wěn)定流場結(jié)果作為下次“熱啟”計算的初始場條件。

表1 水動力與泥沙模擬概況及驗證項目

3.4 參數(shù)率定與驗證結(jié)果

通過后 “0102”枯水期實測資料對模型糙率(EFDC模型中表現(xiàn)為摩擦高度)進行了率定，結(jié)果在在0.001～0.02 m之間，其中糙率值自上游河段向下游河段逐漸減?。簧嫌魏拥涝?.001 5～0.02 m之間，下游各河段在0.001 0～0.001 5 m之間；在水位模擬方面，驗證的 7個站點的水位誤差值都在10 cm以內(nèi)；流量結(jié)果的相對誤差較低，30%誤差以下的模擬時段占總模擬時段的67.7%～87.5%；斷面平均流速與流向與實測都基本吻合，說明水動力模型能夠客觀真實地反映河網(wǎng)內(nèi)的水流運動狀態(tài)，可以為后續(xù)的研究提供可靠的動力基礎(chǔ)。代表斷面的流量、水位及流向驗證如圖2-圖4。

圖2 1999年7月浮標廠斷面水位(a)、流量(b)驗證

圖3 1999年7月黃沙斷面水位(a)、流量(b)驗證

圖4 2001年1月4日至11日沙磚廠斷面平均流速(a)與流向驗證(b)

基于枯季實測泥沙過程，對泥沙模型中的沉降速率WSED0、臨界沉積應(yīng)力TAUD、臨界再懸浮速率WRSP0以及臨界沖刷應(yīng)力TAUR進行率定。泥沙顆粒的沉降速率在斯托克斯靜水沉速公式基礎(chǔ)上率定得到。臨界沖刷應(yīng)力與臨界沉積剪切應(yīng)力經(jīng)率得到的比值為1.6倍，與建議值1.2倍較為接近[18-19]，本文在各參數(shù)取值為如表2所示。泥沙模型對各站點模擬的斷面平均泥沙濃度與實際值的相對誤差在27.2%～43.5%之間，平均誤差為35.5%。代表性斷面懸沙驗證如圖5。

表2 泥沙模型主要參數(shù)率定結(jié)果

圖5 1999年7月黃沙站(a)與浮標廠站(b)斷面平均懸沙濃度驗證

4 結(jié)果與討論

4.1 廣州珠江河網(wǎng)水流與輸沙特征分析

廣州珠江河網(wǎng)縱橫，水道貫通，受徑流與潮汐雙重影響，水流運動十分復(fù)雜；加之歷史上人為河道窄縮工程，使得河網(wǎng)岸線發(fā)生一定程度的形變，進一步影響了不同河段水流運動特點。其中，前航道形成的獨特的“上窄下寬”的喇叭型結(jié)構(gòu)，使得水流在前航道上段流速相對較大，而行至下端，尤其是獵德大橋以下時，河道逐漸放寬，水流減緩。以“9907”典型洪水期為例，上段海珠橋附近斷面絕對流速平均值約為0.46 m/s，而在獵德大橋附近時，流速降至0.39 m/s；至下游“喇叭口”附近時，流速進一步降至約0.26 m/s。由于水流流速沿程逐漸降低，也勢必對水體懸沙的輸運與淤積過程產(chǎn)生影響。模擬發(fā)現(xiàn)，南航道漲潮流速大于落潮流速。以浮標廠斷面為例，7月16日至20日之間，模擬的斷面平均最大漲潮流速在0.54～0.74 m/s之間，落潮最大流速在0.22～0.34 m/s之間，與實測值基本一致[1]。模擬結(jié)果還發(fā)現(xiàn)，上溯水流還能長驅(qū)直入白坭水道，并通過落潮過程分流進至前航道。

在河網(wǎng)輸沙方面，由于白坭水道上游修建了蘆苞涌、西南涌水閘的修建，以往北江分流水沙通量均被大幅截斷，北江分流至平洲水道的水流攜沙是珠江河網(wǎng)廣州段水體懸沙主要來源。從河網(wǎng)水期輸沙量上看，“9907”洪季平洲水道經(jīng)沙洛圍斷面輸沙量達108.28萬噸，而“0102”枯季僅為2.96萬噸[3]。由此可見，北江向廣州珠江河網(wǎng)洪季的側(cè)向分流能夠深刻影響到河網(wǎng)的水量以及泥沙濃度，把握住洪季輸沙過程對于研究廣州珠江河網(wǎng)泥沙輸運與分布特征至關(guān)重要。洪季模擬結(jié)果表明，高濃度的懸浮泥沙主要輸沙路徑是經(jīng)由平洲水道進入廣州珠江河網(wǎng)，并沿瀝滘水道、三枝香水道再進入后航道，并分別經(jīng)新造水道、官洲水道向黃埔方向輸運(如圖6所示)。同時，平洲水道注入河網(wǎng)區(qū)的高濃度懸沙在潮汐頂托作用下，乘漲潮過程沿南航道上溯。大潮時，能夠穿過白鵝潭汊點進入西航道。當發(fā)生落潮時，高濃度泥沙水團又經(jīng)由白鵝潭分流進入前航道，形成間歇性的高濃度含沙水團，這與朱繼偉等[7]在對廣州內(nèi)港航道的水體懸沙分布分析中的發(fā)現(xiàn)的珠江南航道—前航道泥沙輸運過程是一致的。這一特點將使得北江泥沙能夠在潮汐過程中分布到前后航道的各個河段水體，并在落淤過程沉積在河段底泥，這也勢必影響到泥沙表面的各類吸附態(tài)污染物在河網(wǎng)中的輸移和分布。

圖6 廣州珠江河網(wǎng)洪季8日平均懸沙濃度場

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

為了進一步分析廣州珠江河網(wǎng)河床演變規(guī)律與趨勢，對珠江河網(wǎng)在典型洪、枯季河床沖淤過程進行分段統(tǒng)計。參考朱繼偉等[7]對廣州內(nèi)河航道1964-2006年近50年的河床演變過程的河網(wǎng)統(tǒng)計中的劃分方法，本文同樣將研究河網(wǎng)區(qū)劃分為白坭河、白沙河、沙貝海、前航道、南河道、瀝滘水道、官洲河以及新造水道九河段(如圖7)，并以“9907”和“0102”典型洪、枯季河床沖淤特征及沖淤量進行模擬，分別計算各時段的沖淤量，并按等比倍數(shù)疊加得到全年河床形變量。

圖7 廣州珠江河網(wǎng)沖淤分析河道分區(qū)示意圖

圖8 1999年7月(a)與2001年2月(b)洪、枯季廣州珠江河床沖刷與淤積量,m

模擬的珠江河網(wǎng)沖淤結(jié)果如圖8所示。以水體懸沙模擬結(jié)果穩(wěn)定后的168 h(共模擬192 h)的河床演變量為基礎(chǔ)，統(tǒng)計7 d內(nèi)各網(wǎng)格的河床底泥容積變化量，通過加和各單元沖淤體積得到河段總沖淤量后，最后除以各河段總面積得到各河段平均沖淤厚度。從洪季模擬的河床沖淤態(tài)勢上看，廣州珠江河網(wǎng)主要以淤積為主。西航道的水口水道和白坭水道的鴉崗與硬頸海段、前航道“喇叭口”末段、南航道、瀝滘水道、三枝香以及官洲河下段、新造水道上段均以淤積趨勢為主。其中，以南航道淤積速率最大，達到18 cm/a。如前所述，受潮汐頂托作用，西北江經(jīng)平洲水道分流分沙與漲潮水沙疊加后涌入南航道，漲潮流速大，使得高濃度泥沙水流能長驅(qū)上溯，將南航道河段平均泥沙濃度維持在較高的水平(如圖6)；憩潮時水體懸沙大量落淤，形成較高的淤量。前航道上段及中段河道由于河段相對窄縮，水流較快，不易發(fā)生淤積而主要表現(xiàn)為沖刷，末端河口表現(xiàn)為淤積特點?？菁緯r，河網(wǎng)主要以沖刷過程為主，白坭水道、白沙河、沙貝海、南航道—新造水道以及前航道均表現(xiàn)出明顯的沖深特點。其中以新造水道、白沙河為最，年沖刷率可達6 cm/a和7 cm/a。

整體上看，廣州珠江河網(wǎng)(特指本節(jié)納入統(tǒng)計的九條河段)表現(xiàn)為洪季淤積枯季沖刷的特點，洪季淤積速率約為5.7 cm/a，枯季沖刷速率約為1.7 cm/a，全年呈緩慢淤積趨勢，速率約為2.0 cm/a(全河段淤積量除以河段總面積)。

圖9 珠江河網(wǎng)年沖淤厚度預(yù)測值與實測值比較

模擬的9條河段河床的沖淤厚度(正數(shù)表示淤積，負數(shù)表示沖刷)及與前人研究的實際沖淤量對比如圖9所示?？梢姡Ｐ陀嬎愕玫降闹榻骱佣螞_淤厚度與前人的實測結(jié)果較為接近。其中，在河網(wǎng)主要淤積河段——南航道平均淤積速率模擬值約為8.1 cm/a，與實測淤積速率8.9 cm/a較為吻合，與前人模擬結(jié)果(8 cm/a)也十分接近[3]。瀝滘水道淤積速率約為8.5 cm/a，與實測結(jié)果(9.3 cm/a)也比較接近。然而，白沙河和沙貝海兩個河段模擬的沖淤狀態(tài)與實測結(jié)果正好相反，并且白鵝潭和新造水道的模擬結(jié)果與實測結(jié)果也存在較大的誤差。經(jīng)分析，這可能有以下幾點原因：首先限于研究資料的不足，本文尚不能對河網(wǎng)全年水流輸沙過程進行模擬，而是采用兩個不同年份洪、枯季節(jié)的沖淤過程加權(quán)得到，因此無法完整地反映河床長時間尺度上的形變過程，而且引用的實測值也是前人對1999-2006年(或2005年)河網(wǎng)平均沖淤值，因此一定程度上影響到計算結(jié)果的精確性。其次，沙貝海雖表現(xiàn)為淤積趨勢，但1999-2006年，沙貝海河道的平均河寬減少了近10 m，河道過水面積減少約36.1 m2，因此從斷面過水面積上看，該河段體現(xiàn)為淤積，但其水深卻有所增大，平均水深從4.49 m增加至4.54 m，年均沖深約1 cm[7]，這與本文模擬的趨勢是一致的。此外，新造河道的模擬沖刷量小于實測值，這可能歷史上與人為的采砂活動和航道疏浚工程有關(guān)[3]。如1998年洪、枯季，珠江干流分別就有8艘采砂船從事采砂活動(其中有7條船是在新造水道作業(yè))，年采砂量達161萬噸[13]。這對新造水道的河道地形勢必造成一定的影響，也可能是模擬值小于實際值的原因之一。從整體的沖淤模擬結(jié)果上看，本文使用的泥沙模型基本上把握住了珠江河網(wǎng)西航道、南航道及瀝滘水道的淤積的主要特點，計算得到?jīng)_淤量與實際值總體上比較接近，能初步反映出珠江河網(wǎng)河床沖刷淤積的變化過程。

5 結(jié) 論

1)基于EFDC模型，構(gòu)建了廣州珠江河網(wǎng)三維水動力與泥沙數(shù)學(xué)模型。通過1999年7月和2001年2月典型洪、枯季水流泥沙實例驗證，表明模型能較好地模擬和反映河網(wǎng)感潮水流及輸沙特點，計算結(jié)果與實測資料較為吻合。

2)洪季北江經(jīng)平洲水道分流分沙過程是廣州珠江河網(wǎng)主要泥沙通量來源，主要的輸沙路徑集中在后航道瀝滘水道—新造水道、官洲河—黃埔線。而受潮汐能量頂托作用，高濃度懸沙能經(jīng)由南航道上溯至白坭水道，并在落潮時通過白鵝潭分流進入前航道，形成間歇性高濃度泥沙團，使得北江分沙得以在廣州珠江河網(wǎng)中廣泛分布。

3)通過預(yù)測主要河段的沖淤趨勢及厚度，結(jié)果顯示：21世紀初期，廣州珠江河網(wǎng)總體呈緩慢淤積態(tài)勢，年淤積速率約為2.0 cm/a；且在不同水期主要表現(xiàn)為“枯沖洪淤”的特點。其中，南航道、瀝滘水道是主要的淤積河段，速率約在8～9 cm/a；而沙貝海、新造水道是主要沖刷河段，但人為窄縮河道、航道疏浚及人為采砂活動對河道形變趨勢都能產(chǎn)生一定影響。

4)由于廣州珠江河網(wǎng)資料相對較少，模擬時段較短，本文僅對河網(wǎng)底泥沖淤分布進行了初步的研究。在今后的研究中仍需進一步收集長時間尺度的水文泥沙資料，對河網(wǎng)河床演變過程進行更深入的模擬分析。

參考文獻：

[1]包蕓，姚冬靜.珠江網(wǎng)河二維水動力數(shù)值計算研究[J].水動力學(xué)研究與進展，2005，20(增刊)：881-886.

[2]包蕓，劉歡.風(fēng)暴潮極值狀態(tài)下珠江河網(wǎng)白坭水道內(nèi)的不穩(wěn)定流動[J].大連大學(xué)學(xué)報，2007，28(3)：22-27.

[3]何用，徐峰俊.珠江洲頭咀河段過江隧道工程極限沖刷模擬研究[J].人民珠江，2007(5)：22-25.

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

[5]羅宗業(yè)，珠江洪水上漲——廣州水侵街[J].海岸工程，1994，13(4)：36-43.

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

[7]朱繼偉，武亞菊，劉俊勇,等.廣州內(nèi)港航道近50年來河床演變淺析[J].人民珠江,2008(3):15-17.

[8]王興奎，邵學(xué)軍，李丹勛.河流動力學(xué)基礎(chǔ)[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]羅憲林,楊清書，賈良文,等.珠江三角洲網(wǎng)河河床演變[M].廣州： 中山大學(xué)出版社,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]齊珺， 楊志峰，熊明,等.長江水系武漢段水動力過程三維數(shù)值模擬[J].動力學(xué)研究與進展，2008，23(2)：212-219.

[17]陳景秋，趙萬星，季振剛.重慶兩江匯流水動力模型[J].水動力學(xué)研究與進展，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.