林高松，廖國威，唐天均

(深圳市環(huán)境科學(xué)研究院，深圳 518001)

東深供水工程是為解決香港淡水供應(yīng)困難而興建的一項(xiàng)跨流域大型引水工程。目前，東深供水工程每年對香港、深圳和東莞沿線八鎮(zhèn)的供水量分別占其用水總量的75%、50%和80%以上，供水安全對香港、深圳和東莞等地社會經(jīng)濟(jì)發(fā)展和人民健康具有極其重要的意義。

東深供水工程取水口位于東莞橋頭鎮(zhèn)太園泵站。取水口水質(zhì)不僅受到東江上游來水影響，還受石馬河溢流的威脅。石馬河是東江的一級支流，發(fā)源于深圳市寶安區(qū)大腦殼山，河長約88 km，集水面積1 249 km2。石馬河受納了沿岸深圳、東莞的生活污水和生產(chǎn)廢水，水質(zhì)長期劣于地表水Ⅴ類標(biāo)準(zhǔn)。為防止石馬河河水排入東江而影響飲用水源安全，目前實(shí)施了石馬河調(diào)污工程，即在河口用橡膠壩阻攔將河水通過地下暗管排入小海河，流經(jīng)企石水閘入東莞運(yùn)河，從虎門出?？谂懦鯷1]。但橡膠壩截污能力有限，石馬河流量較大時(shí)河水會溢流過橡膠壩而進(jìn)入東江干流。石馬河河口雖處太園泵站取水口下游，但距離取水口僅300 m左右，如果發(fā)生突發(fā)性污染事故導(dǎo)致石馬河大量污染物溢流進(jìn)入東江，并遇到潮流上溯，則污染物將會影響太園泵站取水口，嚴(yán)重威脅東深供水的水質(zhì)安全。

石馬河影響東深供水工程取水口的前提是潮流上溯。東江水利樞紐工程位于取水口上游約28 km，必要時(shí)可通過調(diào)節(jié)下泄流量來控制潮流上溯強(qiáng)度。理論上下泄流量增加至某個(gè)“臨界值”后，取水口附近潮流不會上溯，從而避免石馬河對取水口的水質(zhì)污染。科學(xué)合理確定下泄流量的“臨界值”，對于制定水利應(yīng)急調(diào)度方案和保障東深供水水質(zhì)安全具有重要意義。

水利工程調(diào)度是處置突發(fā)性水污染事件的一個(gè)重要手段[2-5]。丁洪亮等[2]利用平面二維水動力水質(zhì)模型，分析了漢江丹襄段不同位置發(fā)生水污染事故時(shí)，丹江口水庫不同應(yīng)急調(diào)度方式對污染物的稀釋摻混作用；余真真等[3]研究了小浪底水庫運(yùn)行方式對污染事件的調(diào)控能力及實(shí)施效果。部分學(xué)者也針對突發(fā)性污染事故的東江水庫調(diào)度問題開展了研究[6,7]。例如，江濤等[6]模擬分析了不同下泄流量時(shí)石馬河對東深供水工程取水口的水質(zhì)影響，但其計(jì)算僅采用汛期一個(gè)典型潮位過程，未考慮下游各類復(fù)雜的潮汐條件，無法反映出最不利水文條件的影響，而且建模的水文驗(yàn)證數(shù)據(jù)較少，對結(jié)論的準(zhǔn)確性也產(chǎn)生一定影響。此外，現(xiàn)有研究均未進(jìn)一步探討調(diào)度水量的“臨界值”問題，即最不利條件下避免取水口附近潮流上溯所對應(yīng)的最小下泄流量。

本文重點(diǎn)研究東江水利樞紐工程下泄流量對取水口附近潮流的影響，采用實(shí)測數(shù)據(jù)統(tǒng)計(jì)分析并結(jié)合數(shù)值模擬方法，計(jì)算出最不利潮汐條件下取水口潮流不上溯所對應(yīng)的下泄流量臨界值，為東江水利樞紐工程應(yīng)急水量調(diào)度提供決策依據(jù)。

1 取水口斷面潮流上溯特征分析

東深供水工程取水口的水文狀況是下游潮汐和上游徑流動力綜合作用的結(jié)果。上游徑流量主要受東江水利樞紐工程控制；下游潮汐是南海潮波通過伶仃洋-獅子洋-東江河口傳導(dǎo)進(jìn)入，屬于不規(guī)則半日潮[8]。由于獅子洋等還受東江、北江、珠江的徑流影響，潮汐狀況較為復(fù)雜。為反映不同類型的下游潮汐和上游徑流作用，篩選出潮流上溯對應(yīng)的最不利水文組合，采用2017年3月22日-6月12日逐時(shí)實(shí)測流量數(shù)據(jù)進(jìn)行分析，站點(diǎn)包括代表取水口斷面的東岸站和代表東江水利樞紐工程的博羅站，剔除部分缺失數(shù)據(jù)后，共有樣本1 922 組。樣本的時(shí)間跨度82 d，包括了多種潮汐條件，博羅站流量范圍從129～2 130 m3/s，平均值661 m3/s，基本代表了不同上游流量的影響。盡管樣本無法反映出所有水文組合特征，但對上下游綜合動力作用也具有一定的代表性。

根據(jù)統(tǒng)計(jì)，東岸站潮流上溯在大中小潮均可能發(fā)生，取決于上游徑流和下游潮汐強(qiáng)弱對比。樣本中東岸潮流上溯的比例為9.8%，平均上溯流量為203 m3/s；最大連續(xù)上溯累計(jì)流量為1 703 m3/s(2017年4月10日0∶00-2∶00)，期間博羅站平均流量為541 m3/s；最長連續(xù)上溯時(shí)間持續(xù)達(dá)4 h。

實(shí)測的東岸流量與博羅下泄流量的對應(yīng)關(guān)系如圖1所示。由于東岸流量還受下游潮汐影響，因此兩者并沒有顯著的相關(guān)關(guān)系，但從圖1仍可以看出，總體上隨著博羅流量增加，東岸流量也相應(yīng)增大。值得注意的是，東岸流量為負(fù)值(潮流上溯)所對應(yīng)的博羅站流量往往較低，范圍為128～880 m3/s，即樣本顯示博羅流量超過880 m3/s時(shí)東岸流量均為正值，未發(fā)生潮流上溯現(xiàn)象?？紤]樣本的有限性，對最右下側(cè)兩點(diǎn)進(jìn)行連線，確保所有點(diǎn)均在連線左上方(見圖1)，粗略判斷博羅流量大于1 165 m3/s時(shí)可確保東岸流量為正。

圖1 實(shí)測東岸流量與博羅流量的對應(yīng)關(guān)系Fig.1 The relationship between measured Dongan's flow and Buoluo's flow

2 水動力模型

為盡可能全面反映所有水文組合特征和定量研究分析下泄流量對取水口潮流的影響，本文采用丹麥水力研究所(DHI)研發(fā)的MIKE21FM建立東江下游水動力模型。

2.1 模型基本原理

MIKE 21 FM水動力模型的控制方程為基于Boussinesq假定和流體靜壓假定的二維不可壓雷諾平均N-S方程，即淺水方程。

(1)

(2)

(3)

MIKE 21 FM對控制方程的空間離散采用基于網(wǎng)格中心的有限體積法。計(jì)算域采用非結(jié)構(gòu)化網(wǎng)格進(jìn)行概化，網(wǎng)格可以是三角形、四邊形的混合網(wǎng)格。非結(jié)構(gòu)化網(wǎng)格不僅可以對復(fù)雜幾何地形提供最優(yōu)程度的擬合，對邊界進(jìn)行光滑處理，而且可以在重點(diǎn)區(qū)域布置較小的網(wǎng)格單元，非重點(diǎn)區(qū)域布置較大的網(wǎng)格單元，在提高預(yù)測精度的同時(shí)也提高了計(jì)算速度。

2.2 模擬范圍

模擬范圍包括了東江水利樞紐工程以下的東江干流段、東江河網(wǎng)和獅子洋(圖2)。上游邊界包括東江干流的水利樞紐工程和曾江的麒麟嘴，下游邊界包括黃埔和虎門。

圖2 計(jì)算范圍Fig.2 Location of study area

本文的研究重點(diǎn)是水利樞紐工程至石龍的東江干流河段。但考慮到石龍受徑流影響較大，如果模擬范圍下邊界僅至石龍，則建模后作情景計(jì)算時(shí)上游邊界流量可能會影響下游邊界潮位，導(dǎo)致結(jié)果不準(zhǔn)確。因此，將計(jì)算范圍擴(kuò)大至獅子洋。需要指出，開邊界理論上還應(yīng)包括沙灣水道，但由于缺乏該邊界實(shí)測潮位數(shù)據(jù)，因此簡化作陸邊界處理。由于沙灣水道距離東江干流較遠(yuǎn)，相互作用不明顯，這種處理方式并不影響東江干流模擬結(jié)果的準(zhǔn)確性。

網(wǎng)格劃分采用了三角形、四邊形的混合網(wǎng)格，其中東江干流段尤其是取水口附近作加密處理，計(jì)算范圍內(nèi)共劃分網(wǎng)格9 620 個(gè)。

2.3 數(shù)據(jù)來源與參數(shù)取值

東江干流及河網(wǎng)區(qū)的地形采用2016年航道圖，獅子洋的地形采用2015海圖。上游(東江水利樞紐工程、麒麟嘴)邊界條件采用流量逐時(shí)實(shí)測值，下游(黃埔、虎門)邊界條件采用潮位逐時(shí)實(shí)測值。時(shí)間從2017年1月1日至6月30日。

模型的渦黏系數(shù)選取Smagorinsky公式，Smagorinsky系數(shù)取0.28。曼寧系數(shù)根據(jù)實(shí)測水文數(shù)據(jù)進(jìn)行校正，取值范圍30～70 m1/3/s。

2.4 模型驗(yàn)證

采用2017年3月16日-6月12日東岸的流量和水位實(shí)測數(shù)據(jù)進(jìn)行驗(yàn)證，部分驗(yàn)證結(jié)果如圖3和圖4所示(4月3日0∶00-4月11日0∶00)。東岸水位模擬結(jié)果與實(shí)測值吻合較好，平均誤差0.06 m；流量平均誤差112 m3/s?？傮w而言，模型誤差在可接受范圍內(nèi)，并能較準(zhǔn)確地反映出研究區(qū)域的水動力特征。

圖3 東岸斷面水位驗(yàn)證Fig.3 The verification result of water level of Dongan section

圖4 東岸斷面流量驗(yàn)證Fig.4 The verification result of flow of Dongan section

3 潮流上溯影響因素模擬分析

3.1 東江水利樞紐工程下泄流量的影響

為定量分析上游流量對取水口潮流的影響，采用模型模擬水利樞紐工程不同下泄流量對應(yīng)的潮流狀況。模擬時(shí)段從2017年1月1日至6月30日，時(shí)間上涵蓋了旱季、雨季，基本包括各種下游潮汐邊界條件。本研究重點(diǎn)關(guān)注石馬河污染取水口的潛在風(fēng)險(xiǎn)，為體現(xiàn)最不利潮汐條件，潮流上溯強(qiáng)度采用模擬時(shí)段內(nèi)取水口斷面逐時(shí)流量的最小值表示，流量越小則潮流上溯強(qiáng)度越大。東江水利樞紐工程下泄流量分別取800、1 000、1 200、1 300、1 400和1 500 m3/s幾種情況計(jì)算，結(jié)果如圖5所示。

圖5 模擬下泄流量與潮流上溯強(qiáng)度關(guān)系Fig.5 The relationship between the simulated discharge flow and the trend of upstream

由圖5可見，下泄流量與取水口斷面潮流上溯強(qiáng)度總體上呈線性負(fù)相關(guān)關(guān)系，即潮流上溯強(qiáng)度隨下泄流量增加而減小。當(dāng)下泄流量增加至1 400 m3/s左右時(shí)，取水口斷面最小流量為0，即水流不再上漲。經(jīng)試算，下泄流量臨界值為1 380 m3/s。

3.2 石馬河流量的影響

石馬河河口距離取水口僅300 m，而連續(xù)降雨情況下石馬河流量較大[1]。為反映汛期石馬流量對取水口斷面水文的影響，這里假定下泄流量為1 380 m3/s，分石馬河無溢流和石馬河溢流量為200 m3/s兩種情況，計(jì)算取水口的水文狀況，最不利下游邊界邊界條件下的流量過程如圖6所示。

圖6 石馬河溢流對取水口斷面流量影響Fig.6 Influence of Shima River overflow on flow of the water intake section

由圖6可見，石馬河溢流量為200 m3/s情況下取水口斷面的流量最小值略有上升，總體而言兩種情況下取水口斷面流量過程變化并不明顯。主要原因是取水口位于石馬河上游，而石馬河溢流主要影響其下游的水文過程。因此，東江水利樞紐工程下泄流量增加至1 380 m3/s后，石馬河溢流并不會改變?nèi)∷诔绷髁飨?，甚至對潮汐上溯略有抑制作用?/p>

江濤等[6]針對石馬河泄洪對東深供水工程取水口水質(zhì)的影響進(jìn)行了研究分析，結(jié)果顯示石馬河泄洪量越大則東江水利樞紐下泄流量越大才能使取水口氨氮達(dá)標(biāo)，其中石馬河泄洪量為200 m3/s時(shí)，下泄流量需達(dá)到1 770 m3/s[6]。本文研究表明，在最不利潮汐條件下，當(dāng)東江水利樞紐工程下泄流量達(dá)到1 380 m3/s時(shí)，即可抑制取水口附近潮流上溯，此時(shí)石馬河泄洪量并不會改變?nèi)∷诔绷髁飨?，也不會影響取水口水質(zhì)。本文研究結(jié)果與文獻(xiàn)[6]結(jié)論存在較大差異。

根據(jù)2017年3-6月的逐時(shí)實(shí)測數(shù)據(jù)，東岸流量為負(fù)值(潮流上溯)時(shí)的博羅站流量范圍為128～880 m3/s，即下泄流量超過880 m3/s時(shí)取水口附近均未發(fā)生潮流上溯現(xiàn)象?？紤]樣本的有限性，本文第1節(jié)對實(shí)測數(shù)據(jù)進(jìn)行了統(tǒng)計(jì)分析，初步判斷博羅流量大于1 165 m3/s時(shí)可確保取水口流量為正。監(jiān)測數(shù)據(jù)的下泄流量范圍從129～2 130 m3/s，基本包括了不同下泄流量的影響，雖然不能反映出所有水文組合特征，但仍具有較強(qiáng)的代表性。

本文模型模擬結(jié)果(1 380 m3/s)與實(shí)測數(shù)據(jù)的統(tǒng)計(jì)分析結(jié)果(1 165 m3/s)比較接近。由于模型模擬的周期更長，可以更方便和全面地設(shè)置各種水文條件進(jìn)行分析，計(jì)算結(jié)果會更偏保守(即臨界值更大)。因此，結(jié)合實(shí)測數(shù)據(jù)的分析結(jié)果與模型模擬結(jié)果，可認(rèn)為本文的研究結(jié)論是合理的。

3.3 河床地形的影響

東江下游受上游建庫及人工采沙等人類活動的強(qiáng)烈影響，河床高程顯著降低，已成為影響潮流上溯的一個(gè)重要因素[9,10]。假定東江干流從水利樞紐工程至石龍段河床高程分別下降2和6 m并進(jìn)行計(jì)算，則潮流上溯強(qiáng)度明顯增加。試算結(jié)果顯示，東江水利樞紐工程下泄流量需分別增加至1 550和1 660 m3/s，才能確保取水口附近不發(fā)生潮流上溯現(xiàn)象。

4 結(jié) 語

本文利用現(xiàn)有水文監(jiān)測數(shù)據(jù)并結(jié)合數(shù)值模擬方法，針對東江水利樞紐工程下泄流量對東深供水工程取水口潮流上溯的影響進(jìn)行了研究分析，結(jié)果表明：

(1)東深供水工程取水口的水文狀況與下游潮汐、上游徑流、地形等有關(guān)。東江水利樞紐工程的下泄流量對取水口潮流上溯具有抑制作用。當(dāng)下泄流量達(dá)到1 380 m3/s時(shí)，取水口附近潮流不再上溯。如果石馬河發(fā)生嚴(yán)重溢流事故，建議東江水利樞紐工程下泄流量不低于1 380 m3/s，確保避免潮流上溯影響取水口水質(zhì)。

(2)下泄流量大于取水口潮流上溯的臨界值(1 380 m3/s)時(shí)，石馬河溢流流量增加并不會改變?nèi)∷诔绷髁飨?，因此溢流污染物也不會影響取水口水質(zhì)。

(3)潮流上溯強(qiáng)度與地形密切相關(guān)，如果東江干流河床分別下切2和6 m，則為防止取水口潮流上溯的東江下泄流量分別增加至1 550和1 660 m3/s。

□

參考文獻(xiàn)：

[1] 高 磊,陳建耀,朱愛萍,等. 基于SCS模型的跨界小流域物質(zhì)通量估算——以東莞石馬河流域?yàn)槔齕J]. 中國環(huán)境科學(xué),2015,(3):925-933.

[2] 丁洪亮,張洪剛. 漢江丹襄段水污染事故水庫應(yīng)急調(diào)度措施研究[J]. 人民長江,2014,(5):75-78.

[3] 余真真,張建軍,馬秀梅,等. 小浪底水庫應(yīng)急調(diào)度對下游水污染事件的調(diào)控[J]. 人民黃河,2014,(8):73-75.

[4] 馮民權(quán),范世平,楊建明,等. 基于非恒定流的污染物遷移擴(kuò)散隨機(jī)模擬[J]. 自然災(zāi)害學(xué)報(bào),2011,(5):11-17.

[5] 余昭輝,夏建新,任華堂. 水庫不同調(diào)度措施對突發(fā)污染帶遷移的影響[J]. 人民黃河,2015,(3):84-88.

[6] 江 濤,鐘 鳴,鄒隆建,等. 石馬河泄洪與東江水利樞紐調(diào)節(jié)不同情景下東江水質(zhì)的模擬與分析[J]. 中山大學(xué)學(xué)報(bào)(自然科學(xué)版),2016,(2):117-123.

[7] 黃國如,李 春. 東江干流突發(fā)性重金屬污染的水庫調(diào)度效應(yīng)研究[J]. 水電能源科學(xué),2013,(6):82-85.

[8] 方升豹. 東江三角洲潮流特性分析[J]. 人民珠江,1990,(2):28-31.

[9] 季榮耀,陸永軍,左利欽. 東江下游博羅河段人類活動影響下的河床演變[J]. 泥沙研究,2010,(5):48-54.

[10] 賈良文,羅章仁,楊清書,等. 大量采沙對東江下游及東江三角洲河床地形和潮汐動力的影響[J]. 地理學(xué)報(bào),2006,(9):985-994.