林高松，廖國威，唐天均

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

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

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

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

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

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

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

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

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

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

圖1 實測東岸流量與博羅流量的對應關(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)格中心的有限體積法。計算域采用非結(jié)構(gòu)化網(wǎng)格進行概化，網(wǎng)格可以是三角形、四邊形的混合網(wǎng)格。非結(jié)構(gòu)化網(wǎng)格不僅可以對復雜幾何地形提供最優(yōu)程度的擬合，對邊界進行光滑處理，而且可以在重點區(qū)域布置較小的網(wǎng)格單元，非重點區(qū)域布置較大的網(wǎng)格單元，在提高預測精度的同時也提高了計算速度。

2.2 模擬范圍

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

圖2 計算范圍Fig.2 Location of study area

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

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

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

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

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

2.4 模型驗證

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

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

圖4 東岸斷面流量驗證Fig.4 The verification result of flow of Dongan section

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

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

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

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

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

3.2 石馬河流量的影響

石馬河河口距離取水口僅300 m，而連續(xù)降雨情況下石馬河流量較大[1]。為反映汛期石馬流量對取水口斷面水文的影響，這里假定下泄流量為1 380 m3/s，分石馬河無溢流和石馬河溢流量為200 m3/s兩種情況，計算取水口的水文狀況，最不利下游邊界邊界條件下的流量過程如圖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ì)的影響進行了研究分析，結(jié)果顯示石馬河泄洪量越大則東江水利樞紐下泄流量越大才能使取水口氨氮達標，其中石馬河泄洪量為200 m3/s時，下泄流量需達到1 770 m3/s[6]。本文研究表明，在最不利潮汐條件下，當東江水利樞紐工程下泄流量達到1 380 m3/s時，即可抑制取水口附近潮流上溯，此時石馬河泄洪量并不會改變?nèi)∷诔绷髁飨?，也不會影響取水口水質(zhì)。本文研究結(jié)果與文獻[6]結(jié)論存在較大差異。

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

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

3.3 河床地形的影響

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

4 結(jié) 語

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

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

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

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

□

參考文獻：

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

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

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

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

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

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

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

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

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

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