張 琳,孫 娟,翟家佳,顏芬芬,張 奕

(1.長江科學院 河流研究所，武漢 430010；2.南京大學 金陵學院，南京 210089；3.河海大學 環(huán)境學院，南京 210098)

?

水庫突發(fā)性污染事故水質影響過程分析

張 琳1,孫 娟2,翟家佳3,顏芬芬3,張 奕3

(1.長江科學院 河流研究所，武漢 430010；2.南京大學 金陵學院，南京 210089；3.河海大學 環(huán)境學院，南京 210098)

為了有效地控制突發(fā)性水污染事故帶來的危害，以東北某大型水庫為典型研究案例，綜合考慮水庫水文、氣象條件和污染物遷移特征，擬定典型工況，基于已建立的三維水動力和水質耦合模型，模擬分析突發(fā)性污染事故發(fā)生后庫區(qū)污染物三維空間分布及隨時間變化的特征。結果表明：受入庫流量和全年主導風向東北風的影響，總體上污染物隨著水流沿西南岸線向壩址方向擴散，且污染物的超標面積逐漸增大；但隨著流程的增加，污染物的濃度逐漸減小。通過對水庫水質的數(shù)值模擬，定量掌握了水庫水質實時狀態(tài)以及污染物空間分布，可以為有關部門制定與實施應急預案提供相關的技術依據(jù)。

突發(fā)性水污染事故；水質；水庫；污染物三維空間分布；三維數(shù)值模擬

1 研究背景

近年來，突發(fā)性水污染事故在我國重點流域頻繁發(fā)生，由于其發(fā)生的不確定性、處理處置的艱巨性、水污染信息的不完整性與不可比性等特點，往往對水環(huán)境安全造成重大威脅[1]。為了預防和控制潛在的水環(huán)境風險，準確預測水體中污染物的遷移和擴散情況非常必要[2-3]。

隨著水質模型研究的深入，國際上已經(jīng)有很多成熟的水質模型軟件[4]，然而這些模型輸入?yún)?shù)多，分析工作量大，對于突發(fā)性水污染事故風險預測難度較大。本文針對大中型水庫水動力及污染物遷移擴散特征，采用以N-S方程為基礎的紊流運動方程和物質輸運方程建立三維水動力[5]和水質耦合模型[6-7]，能夠快速獲取數(shù)據(jù)，快速模擬污染團的遷移軌跡和空間分布[8]，實時預報校正，為有關部門制定與實施應急預案提供科學依據(jù)，對于制定切實可行的應急措施，有效控制事故的危害，具有重要的技術支撐作用。

2 三維水動力模型數(shù)學描述

2.1 三維水動力數(shù)學模型

笛卡爾坐標系下描述水庫三維水流運動的控制方程組[9-10]如下：

(1)

(2)

(3)

(4)

2.2 三維水質數(shù)學模型

笛卡爾坐標系下污染物三維輸運方程為

(5)

式中：C為污染物濃度；K為污染物衰減系數(shù)；Dv為垂向擴散系數(shù)；Dk為水平擴散系數(shù)。

3 三維水動力及水質模型驗證

3.1 水動力數(shù)值模擬驗證

3.1.1 明渠恒定流垂向流速分布

選用指數(shù)律流速分布經(jīng)驗公式[11]來驗證數(shù)值模擬的結果。設計試驗明渠長5 000 m，寬600 m，水深10 m，東、西邊界分別為上、下游邊界。將水體流動穩(wěn)定后的中心斷面模擬結果與根據(jù)垂向流速指數(shù)律分布公式計算的結果對比進行模型驗證，結果見圖1。經(jīng)統(tǒng)計，數(shù)值解與解析解的絕對誤差范圍為-0.015～0.003 1 m/s，相對誤差范圍為-11.86%～1.23%，模型較好地模擬了明渠水流的流速分布。

圖1 水動力數(shù)值模擬驗證結果Fig.1 Verification results of numericalhydrodynamics simulation

3.1.2 矩形水槽內風生環(huán)流分布

本文選用封閉矩形水槽內靜止水體在風應力的作用下引起的回流對模型進行驗證。設計試驗水槽長1 000 m，寬100 m，水深10 m，垂向分20層，垂向紊動黏性系數(shù)取定值0.012。

底部為無滑移邊界條件，水面風應力為恒定，忽略對流、哥氏力和水平擴散的作用情況下，流速垂向分布解析解為[12]

(6)

式中：u為水平流速(m/s)；σ為相對水深；KM為垂向紊動黏性系數(shù)；τ為均勻風應力；D為總水深。

均勻風應力計算公式為

(7)

式中：ρa為空氣密度；C為風的拖拽系數(shù)；U，V分別為風速在x，y方向的分量。

穩(wěn)定狀態(tài)下，中心斷面垂向流速分布及其與解析解的對比見圖2。經(jīng)統(tǒng)計，數(shù)值解與解析解的絕對誤差范圍為-0.005～0.015 m/s，相對誤差范圍為-9.59%～14.72%，模型較好地模擬了矩形水槽風生環(huán)流水平流速的垂向分布。

圖2 矩形水槽內水動力驗證結果Fig.2 Verification results of hydrodynamics simulation in rectangular flume

3.2 水質數(shù)值模擬驗證

通過理想水槽點源物質輸移擴散數(shù)值模擬與恒定狀態(tài)下非保守物質連續(xù)排放的三維對流擴散方程解析解的比較，驗證水質數(shù)學模型的準確度與精度。設計試驗水槽長10 000 m，寬5 000 m，水深20 m，垂向分10層。東、西邊界分別為上、下游邊界。合理設置邊界條件，保證水槽內為恒定均勻流，污染物排放為連續(xù)排放，排污口位于第5層。模擬足夠長時間，取y=2 500，第5層M水平線為縱向對比線；以排污口為原點，取過點 (1 000，0)的N垂線上的水質分布作對比垂線。對比結果見圖3，數(shù)值解與解析解基本吻合，可見，模型能較好地模擬非保守污染物的輸運過程。

圖3 水質數(shù)值模擬驗證結果Fig.3 Verification results of numerical water quality simulation

4 基于三維數(shù)值模型的庫區(qū)污染物濃度變化特征分析

大伙房水庫位于遼寧省東北部，是中國第一個五年計劃中的第一個大型水庫，是沈陽、撫順2市居民的重要飲用水源地。S202和G1212兩條公路從大伙房水庫渾河入庫區(qū)域北側經(jīng)過，最近點距大伙房水庫83 m，車輛泄漏的危險載運物品有可能流入大伙房水庫，因此選取附近水域為突發(fā)性事故風險源研究水域。具體位置如圖4所示。

圖4 水污染事故研究區(qū)域Fig.4 Research area of pollution accident

大伙房水庫最大水深達41.5 m，屬于深水庫，流速、污染物濃度沿垂向分布差異較大，因此綜合考慮水庫水文、氣象條件和污染物遷移特征[13]，選取水庫實測多年平均入庫流量作為水流邊界條件[14]，以常年平均風速2.3 m/s和主導風向(NE)作為自由表面條件[15]，垂向上采用等距離分層進行剖分，共分為7層，建立水庫三維水動力、水質模型模擬分析突發(fā)性污染事故后庫區(qū)污染物三維空間分布及隨時間變化的特征。

參照貨車載重量，源強取為20 t，污染因子為苯酚，假定污染物瞬時排入庫區(qū)。

圖5 水平流場分布Fig.5 Distribution of horizontal flow fields

4.1 參數(shù)取值

水平渦黏系數(shù)采用Smagorinsky亞網(wǎng)格尺度紊動模型計算，cs取值為0.28；垂向渦黏系數(shù)采用Kolmogorov-Prandl求解，其中k和ε通過求解k-ε方程可以得到；底部粗糙高度取決于底部局部粗糙程度，在缺乏相關資料的情況下，通常取為0.002～0.01 m，根據(jù)相關研究，本文取值為0.01 m。風速選取為春季平均風速3 m/s(NE)。

紊動擴散系數(shù)是影響污染物擴散的主要因素之一，它是由于水流的紊動造成的擴散，紊動擴散與渦黏系數(shù)的比例因子由1/σT估算，σT為施密特數(shù)，取1.0。參考相關文獻，苯酚的衰減系數(shù)取為0.05/d。

4.2 三維水動力特征分析

4.2.1 水平流場特征

在模型計算達到穩(wěn)定狀態(tài)時，取庫區(qū)自由水面層(表層)，第4層(中層)及床面層(底層)的水流流態(tài)進行分析。

如圖5所示，水體流向基本為自渾河入庫口向壩址處流動，表現(xiàn)出良好的順岸流特征，庫區(qū)岸邊相對中部水深較深處流速明顯變大，平均流速約為0.045 m/s。來自上游渾河、蘇子河的水流，一部分流向壩址，一部分流向社河的入庫口，與從社河流入水庫的水流形成順時針的大環(huán)流。中層水體水流結構與表層水體呈現(xiàn)相似的水動力特征，水體流向基本為自渾河入庫口向壩址處流動，受上游來流的影響，渾河、蘇子河入庫口區(qū)域流速較大，在來自上游渾河、蘇子河的水流與從社河流入水庫的水流共同作用下，社河入庫口區(qū)域形成順時針的大環(huán)流。庫區(qū)中段卡口處，水流受彎道和地形的影響形成一個順時針方向的環(huán)流。水庫底部水流結果較為復雜，來自上游渾河、蘇子河的水流，一部分流向壩址，一部分流向社河的入庫口，與從社河流入水庫的水流形成順時針的大環(huán)流，來自社河的水流一部分流向壩址處，一部分在補償作用的影響下，向上游流動，與從渾河、蘇子河向下游流動的水流在水庫凹岸處形成一個順時針方向的環(huán)流。

4.2.2 垂向流場特征

選取斷面1、斷面2、斷面3分別對水庫入庫口、水庫庫區(qū)中段、水庫壩址處的垂向流場進行分析。斷面上流速矢量分布如圖6所示。

圖6 各斷面垂向流速分布Fig.6 Distribution of vertical velocity in different sections

由圖6可知，水庫入庫區(qū)域流速在垂向上存在較大梯度，表現(xiàn)為從表層向底層流速大小依次遞減；縱向流場分布形態(tài)與河床地形形態(tài)相似，水流流態(tài)受地形影響明顯，與地形變化趨勢吻合，在地形變化急劇的地方出現(xiàn)局部垂向分量較大的情況。庫區(qū)中段上游部分垂向水流結構表現(xiàn)為以吞吐流為主的水動力特征，水流方向總體上自上游向壩址處流動，流場分布形態(tài)與河床地形形態(tài)相似，水流流態(tài)受地形影響明顯，與地形變化趨勢吻合，流速大小呈現(xiàn)從表層向底層依次遞減的規(guī)律；庫區(qū)中段下游部分表層至中層水流結構表層為以吞吐流為主的水動力特征，流速大小從上層至下層逐漸減小，中層至底層水流結構表現(xiàn)為以補償流支配為主的水動力特征，流速大小從上層至下層逐漸增大。斷面3垂向流場空間分布較規(guī)則，總體上呈現(xiàn)以吞吐流支配為主的水動力特征，水流方向為自上游向壩址處流動，流速大小從表層向底層逐漸減小，在壩址處，受出水口位置的影響，水流出現(xiàn)向上的流速。

4.3 污染物濃度變化特征分析

4.3.1 平面分布特征分析

水平方向上，取自由水面層(表層)，第4層(中層)及床面層(底層)的苯酚分布狀態(tài)進行特征分析。苯酚最大超標范圍(濃度>0.02 mg/L)包絡線幾何參數(shù)見表1，圖7為事故發(fā)生后不同時刻苯酚表層濃度場分布。

表1 不同時刻污染物濃度超標范圍

如圖7所示，污染物為岸邊排放，在水流輸移作用下，隨著時間的推移，污染團沿東北岸隨流運動，至環(huán)流區(qū)受環(huán)流影響，污染團隨之環(huán)向擴散，但受水體剪切分散的影響，污染團整體仍向壩址方向移動。在此過程中污染團超標范圍(濃度>0.02 mg/L)的面積也逐漸擴大至最大范圍然后減小，在事故發(fā)生57 h后超標面積最大，表層最大值為3.07 km2，中層為3.09 km2，底層為3.11 km2。由表1可知，各層污染團超標范圍面積相差不大，接近垂向均勻混合狀態(tài)。

圖7 不同時刻污染物表層濃度場Fig.7 Surface concentration fields of contaminants at different moments

4.3.2 縱向立面分布特征分析

垂向方向上，取事故點垂線斷面不同時刻濃度分布進行分析。苯酚泄漏進水體，在10 min后即垂向混合均勻，濃度場垂向分層明顯，5 h后事故點垂線位置苯酚濃度屬于達標狀態(tài)(<0.02 mg/L)。事故點苯酚濃度垂向分布見圖8。

圖8 事故點不同時刻縱向立面濃度場Fig.8 Facades of vertical concentration fields at the accident site at different moments

在事故點下游距事故點直線距離分別為1 200，1 500，1 800 m的位置設定監(jiān)測點S1，S2，S3。監(jiān)測點位置見圖4。圖9為事故發(fā)生5 d內監(jiān)測點S1，S2，S3的濃度變化曲線。

圖9 風生流表層監(jiān)測點濃度曲線Fig.9 Curves of concentration at monitoring points on the surface of wind-driven current

由圖9可知，事故發(fā)生后0.5 h左右，苯酚污染團即擴散到S1點，S1點的濃度隨即上升，在4 h左右達到濃度最大值0.84 mg/L，隨后濃度值降低，由于受到環(huán)流的影響，上游受污染的水體運動至S1點，因此在19 h時濃度達到最小極值后出現(xiàn)小幅度的上升，在50 h時左右濃度達到一定極值后再下降至對水質無影響；事故發(fā)生后2 h左右，苯酚污染團擴散至S2點，S2濃度隨即上升，在7h時濃度達到最大值0.42 mg/L，由于該監(jiān)測點水流結構呈現(xiàn)以吞吐流為支配的水動力特征，因此在對流輸移與水體紊動剪切分散作用下，污染物的濃度值逐步降低至對水質幾乎沒有影響；S3點水流結構與S1相似，因此該監(jiān)測點濃度隨時間變化規(guī)律大體同于S1，受地理原因影響，苯酚污染時間發(fā)生相應延遲。因污染物的對流輸移和水體紊動剪切分散作用，污染團中心濃度沿程減小，S2點的最大污染濃度小于S1，S3點小于S2。 圖10為S2不同時刻橫向立面濃度場分布，圖11為S2和S3兩點連線垂向斷面不同時刻垂向濃度場分布。

圖10 S2不同時刻橫向立面濃度場Fig.10 Facades of transverse concentration fields at monitoring point S2 at different moments

圖11 S2-S3不同時刻縱向立面濃度場Fig.11 FacadesoflongitudinalconcentrationfieldsatmonitoringpointS2andS3atdifferentmoments

如圖10、圖11所示，在事故發(fā)生9 h時，污染團中心移動至S2點，S2-S3斷面苯酚濃度開始急速上升，至50 h時，S2苯酚濃度達標，在此過程中，因斷面垂向流速梯度較大，苯酚濃度垂向均勻混合，垂向分層明顯。在120 h后，全斷面均達標(濃度<0.02 mg/L)。

模擬結果顯示，隨著污水團不斷向下游遷移，污水團中心可溶性污染物濃度沿流程逐漸減小。水流結構是影響斷面縱向立面濃度場隨時間變化特征的重要因素，當斷面水流結構呈現(xiàn)以風生環(huán)流為支配的水流特征時，濃度變化呈現(xiàn)從無影響→受到影響→影響增加→達到峰值→逐步衰減→再次影響(影響較小)→影響消失的過程；當水流結構呈現(xiàn)以風生環(huán)流為支配的水流特征時，濃度變化呈現(xiàn)從無影響→受到影響→影響增加→達到峰值→逐步衰減→影響消失的過程。

5 結 語

針對大中型水庫水動力及污染物遷移擴散特征，采用以N-S方程為基礎的紊流運動方程和物質輸運方程建立三維水動力和水質耦合模型。以東北某大型水庫為典型研究案例，該模型有效地模擬了大伙房水庫突發(fā)性污染事故發(fā)生后庫區(qū)污染物三維空間分布及隨時間變化的特征。模擬結果表明，受入庫流量和全年主導風向東北風的影響，總體上污染物隨著水流沿西南岸線向壩址方向擴散，且污染物的超標面積逐漸增大，但隨著流程的增加，污染物的濃度逐漸減小。水流結構是影響斷面縱向立面濃度場隨時間變化特征的重要因素，當斷面水流結構呈現(xiàn)以風生環(huán)流為支配的水流特征時，濃度變化呈現(xiàn)從無影響→受到影響→影響增加→達到峰值→逐步衰減→再次影響(影響較小)→影響消失的過程。該模型的建立，完善了水庫水污染應急反應體系，在一定程度上間接提高了水污染事故應急反應的水平。

[1] 沈永明,鄭永紅,吳修廣.近岸海域污染物遷移轉化的三維水質水動力學模型[J].自然科學進展,2004,14(6):694-699.

[2] 侯國祥，翁立達，葉 閩,等.三峽水庫重慶庫區(qū)水質預測[J].長江科學院院報,2002,19(1):13-16.

[3] 唐建華，劉瑋祎，陶 靜,等.長江口北支倒灌對南支水域水質影響初探[J].長江科學院院報,2012,29(9):14-17.

[4] 黎育紅，賀石磊.淺水湖泊群聯(lián)通與調水的二維水動力-水質耦合模型研究[J].長江科學院院報,2015,32(1):21-27.

[5] DAVYDOV B I. On the Statistical Dynamics of an Incompressible Turbulent Fluid[J]. Soviet Physics Doklady,1961, 4: 769.

[6] 韓龍喜，陸東燕，李洪晶，等.高鹽度湖泊艾比湖風生流三維數(shù)值模擬[J].水科學進展,2011,22(1)：97-103.

[7] HARLOW F H, NAKAYMA P I. Transport of Turbulence Energy Decay Rate[R]. US: Los Alamos Scientific Lab., 1968.

[8] 張防修,王艷平,劉興盛,等.黃河下游突發(fā)性污染事件數(shù)值模擬[J].水利學報,2007,(增刊):613-618.

[9]SPALDING D B. The Prediction of Two-dimensional Steady Turbulent Flows[R]. London: Heat Transfer Section, Imperial College, 1969.

[10]UNESCO. Tenth Report of the Joint Panel on Oceanographic Tables and Standards[J]. UNESCO Technical Papers in Marine Science,1981, (36): 17-21.

[11]黃才安. 明渠流速分布指數(shù)律與對數(shù)律的統(tǒng)一及轉換[J]. 人民長江,1994,25(1)：42-44.

[12]KOUTITAS C, O’CONNOR B. Modeling Three-dimensional Wind-induced Flows[J]. Journal of Hydraulic Division, ASCE, 1980，106(11)：1843-1865.

[13]黃 瑞,張防修,孫 娟,等. 大伙房水庫監(jiān)測斷面水質特征分析[J].水科學與工程技術, 2012,(1):78-81.

[14]張 靜，何俊仕. 大伙房水庫汛期氣象水文特征分析[J]. 人民黃河, 2012, 34(11):45-47.

[15]白玉良，賀國靜，吳玉禾. 大伙房水庫環(huán)境水文氣象效應[J]. 水利管理技術,1994, (3):35-37.

(編輯：姜小蘭)

Numerical Simulation on the Process of Water Quality Influenced bySudden Pollution Accident of Reservoir

ZHANG Lin1, SUN Juan2, ZHAI Jia-jia3, YAN Fen-fen3, ZHANG Yi3

(1.River Department, Yangtze River Scientific Research Institute, Wuhan 430010, China; 2.Jinling College, Nanjing University, Nanjing 210089, China; 3.College of Environment, Hohai University, Nanjing 210098, China)

To prevent from and control the effect of emergent water pollution accident in reservoir area, the temporal and spatial variations of pollutants in the reservoir after the accident were simulated by a 3-D model coupling hydrodynamics and water quality at typical working conditions. Factors of hydrology, meteorology and contaminant move-ment were taken into account. A large reservoir in the northeastern area of China was taken as a case study. Results revealed that affected by inflow and dominant wind in northeast direction, the contaminants spread toward the direction of the dam site along the southwest coastline, and the area of standard-exceeding contaminants increased gradually; but the concentration of contaminants decreased gradually with the flow process. Through simulation of the water quality, the spatial distribution of contaminants is quantitatively obtained, which offers technical basis for emergency treatment.

emergent water pollution accident; water quality; reservoir; spatial distribution of contaminants; 3-D numerical simulation

2014-08-26；

2015-10-20

張 琳(1988-)，女，福建三明人，工程師，碩士，主要從事環(huán)境水力學研究，(電話)15342788110(電子信箱)zlincandy@163.com。

10.11988/ckyyb.20140746

2016,33(10):12-17,23

X24

A

1001-5485(2016)10-0012-06