大型跨流域調水工程對水源地水動力的影響研究

周紅玉，劉 操，吳曉輝

(1.首都師范大學 資源環境與旅游學院，北京 100048；2.北京市水科學技術研究院,北京 100048)

0 引 言

南水北調工程是中國有史以來最大的跨流域調水工程，其中線工程將漢江上丹江口水庫水由南向北調到北京市團城湖，再由團城湖通過京密引水渠輸送至密云水庫。密云水庫為北京市唯一的地表水源地，南水北調來水后密云水庫的水利聯系由原來的由北向南來水變為南北同時通水，必然會對其水動力、水質及其水生態產生一定的影響，而對水動力的影響是最基本、最重要的部分。為確保密云水庫水環境安全，本文就南水北調來水對密云水庫水動力的影響進行模擬研究，為研究南水北調來水對密云水庫水質及水生態的影響做準備。除了實際現場調研外，水文模型是研究該類問題較成熟的手段，而且方便成本低。水文模型對河流、湖泊、水庫等的模擬已發展成為對水環境模擬、研究、預測與治理的較為成熟的研究手段。MIKE21模型采用有限體積計算方法，在國內外水環境研究中得到廣泛的使用，其數值模擬的科學性得到世界公認[1-4]。

MIKE21模型利用二階精度的有限差分方法對連續方程和動量方程求解，可以模擬各種作用力下忽略分層的二維自由表面流的水位和水流變化[5]。本研究搭建2015年MIKE21水動力模型(HD)模擬密云水庫水位、面積、流速、流態等水文信息，通過與實際監測數據對比分析，經過不斷調試，找到適合密云水庫的模擬參數；通過2016年水文數據驗證分析檢驗水動力模型的可靠性；分析模擬了代表性引水方式、出入流量、不同風場等不同條件影響下密云水庫水動力的空間變化特征。本文通過構建密云水庫水動力模型，研究南水北調來水后密云水庫的水文變化情況，為下一步研究南水北調對密云水庫水質影響奠定基礎，為密云水庫的管理和決策提供技術支持。

1 研究區概況

密云水庫位于北京市密云縣境內，是京津唐地區第一大水庫，也是南水北調中線工程的重要調蓄水庫，原具有防洪、灌溉等重要功能，現作為北京市地表飲用水源地。密云水庫上游承接潮河和白河來水[6]，下游北京市第九水廠從中取水，流域面積為15 788 km2,占潮白河流域面積的88%[7]。密云水庫的經緯度范圍大致為北緯40°19′～41°31′，東經115°25′，～117°33′，氣候類型為溫帶大陸性季風氣候，冬夏季時間長，春秋季節短，降雨季節性比較明顯，夏季多雨，其他季節少雨，年降水量變幅較大。密云水庫最高蓄水面積188 km2，最大入庫容量43.75 億m3, 水深40～60 m，但由于近年來氣候干旱，水庫庫容銳減，水庫庫容在10 億m3左右，水位在135 m左右[7]。北京市缺水狀況一直以來是北京市水環境的重要問題，為緩解缺水嚴重問題，南水北調實施后將丹江口水庫水調到北京團城湖，再由京密引水渠輸送到密云水庫。

2 密云水庫模型的搭建

密云水庫有兩大支流，一條是白河，源于河北省沽源縣，流經赤城縣、延慶縣、懷柔區，流入密云水庫；潮河起源于河北省豐寧縣，自古北口入密云水庫。模擬區域的邊界設為進湖兩處(白河和潮河)，白河來水以張家墳水文站每日實測流量控制，潮河來水以下會水文站每日實測流量控制。模擬區域的出口有一處，密云水庫下游出水口是第九自來水廠(簡稱水九)，以水九水文站每日實測水位控制。

南水北調中線工程由丹江口水口引水至北京的團城湖，再由京密引水渠輸送至密云水庫，2014年年底正式通水，計劃在2014年后平均每年為北京帶來10～14 億m3的水資源，平均每年有2 億m3左右的水進入密云水庫。在模型驗證時將南水北調來水設置為進水口，以實際監測的每日流量控制。

2.1 MIKE21模型簡介

MIKE21模型是MIKE模型中的核心基礎模塊，是丹麥水力學研究所(DHI)開發的系列軟件之一。主要應用領域包括湖泊、河流、水庫、海洋、波浪及泥沙等的模擬研究，主要應用于與水有關的工程領域及環境問題方面，包括河流水體等的水動力、水質模擬，水環境功能模擬，海灣波浪、潮汐等的模擬，水環境規劃研究等。

MIKE21模型的控制方程如下：

連續方程：

(1)

動量方程：

(3)

式中：ε為水位高程；p、q分別是x、y方向上的流量；h為水深；t為時間；g為重力加速度；C為謝才系數；Ω為Coriol系數；f為風摩擦力；V為風速；Vx、Vy分別為風在x、y方向上的分量。

2.2 二維水動力模型計算條件設置

二維水動力模型需要輸入的數據包括：密云水庫庫區網格數據(庫區閉合水陸邊界)、庫區水深高程數據、開邊界流量和水位時間序列數據、模型模擬參數、模擬初始水位及流量、模擬時間和時間步長的設定。經過多次調試得出適合密云水庫的二維水動力模型。

2.2.1 模型網格設置

在密云水庫水力聯系中，潮河和白河從上游來水進入到密云水庫，密云水庫位于兩河下游，系攔蓄兩河河水而成，模型將兩河入庫處設置為開邊界，第九水廠出水以點源的形式輸出，其余為閉邊界。模型采用非結構三角網格劃分研究區，三角網格能較好的擬合研究區海陸邊界，且能隨意調整網格密度和網格大小。由密云水庫地形圖得到研究區海陸邊界線數據(Land.xyz)以及水深數據(Water.xyz),然后將海陸邊界線數據(Land.xyz)導入網格生成器Mesh Generator生成研究區網格并設置模型開邊界，模型網格數為41 267個，三角網格的邊長在100 m之內。然后將水深數據(Water.xyz)導入網格中并插值生成研究區地形圖如圖1。

圖1 研究區地形Fig.1 Topographic map of the study area

2.2.2 水動力模型計算條件的設置

MIKE21 水動力模型需要設置的計算條件包括：模擬時間與時間步長、地形、克朗值、干濕邊界、渦黏系數、降雨和蒸發數據、源和匯、邊界條件以及初始條件等。

初始條件：初始條件設為2015年年初水面高程133.18 m，水平速度和垂直速度均為0 m/s。

時間序列：模型中需要添加的時間序列數據有白河和潮河兩個水文站張家墳、下會每日實測流量數據，水九每日水位實測數據以及庫區實測降雨和蒸發數據，時間序列數據由密云水庫管理處提供的2015-01-01～2015-09-30的實測數據。兩個邊界的位置以及第九水廠的位置如圖1所示。

模擬時間與時間步長(time step)：水動力模型的模擬時間為2015年1月1日8∶00至2015年9月30日8:00。模擬時間步長(time step)為100 s，時間步數為235 008步。

克朗值(CFL number)：網格分辨率、水深和時間步長決定了模型設置中的克朗值。克朗值在小于1的情況下模型才能正常運行，所以為了保證模型運行的穩定性將CFL值設置為0.8。

干濕邊界(Flood and dry)：MIKE21 可以根據每個網格的水深情況調整計算條件，靈活的調用公式參與計算[8]。干濕水深是用來判斷單元網格是否參與到模型的計算中來。在該模型中干水深為0.001m,淹沒水深為0.05 m,濕水深為0.1 m。

渦黏系數(Eddy Viscosity)：選用Smagorinsky公式，取值0.28 m1/3/s。

底床摩擦力(Bed Resistance)：底床摩擦力可以選擇謝才數(Chezy Number)、曼寧數(Manning Number)或者不選，這里選用曼寧系數(Manning Number),取值32 m1/3/s。

其他沒有設置的參數采用默認值。

2.3 水動力模型計算結果分析與率定

利用MIKE 21FM水動力學模型對密云水庫庫區進行模擬，將白河庫區的水位結果輸出，白河庫區位置如圖1所示。將白河庫區的實測水位與模擬水位進行比較率定，結果如圖2所示。由圖2可以得出白河庫區的模擬結果與實際結果相比，相對誤差在0.2%以內，絕對誤差在-0.015～0.268 m范圍內。從總體上看，模擬的效果比較好，模擬誤差控制在模型計算要求的范圍內。

圖2 白河庫區模擬水位與實測水位對比圖Fig.2Comparison of simulated and measured water level in Baihe reservoir area

2.4 水動力模型的驗證

為了驗證水動力模型的可靠性，在率定好的水動力模型的基礎上，用2016年1月1日至2016年11月8日水文數據對模型進行驗證，同時將南水北調來水以點源的形式通入密云水庫，采用每日實際監測數據。模型同樣將白河庫區水位輸出，通過與實際監測數據對比，驗證模型的可靠性。如圖3所示，白河庫區實測水位與模擬水位對比，絕對誤差98%在-0.206～0.30 m范圍內。通過驗證，模擬水位與實際水位吻合較好，模型比較可靠。

圖3 2016年白河庫區水位模擬值與實測值對比Fig.3 Comparison of simulated and measured water level in Baihe reservoir area in 2016

3 南水北調來水對密云水庫水文情勢的影響分析

3.1 工況設計

通過對實際調水情況的分析，本模型在綜合考慮當地氣象條件，在不同的風場變化的條件下以及不同的引排水情況分析不同情景對密云水庫水動力的影響。

密云水庫地區常年主導風向為SW風，多年平均風速為2.24 m/s。綜合考慮現有氣象水文資料，本研究設計兩種風場情況分別為2.24 m/s SW風和無風，來驗證密云水庫受不同風場影響下水動力的變化情況。

根據密云水庫水力聯系及調水情況可知，水體流動情況為北部上游白河和潮河來水與西南部南水北調來水匯兌在東南的第九水廠流出。根據實際調水情景可知，現有調水工程設計為調水量在10 m3/s左右,且冬季(12月-次年2月)不調水。為驗證加大調水量以及不調水時不同的流量組合對密云水庫水動力條件的影響，分別設計調水口進水量取0、10、20 m3/s。考慮到調水量增大后第九水廠的取水量相應的增加，水九的出水流量取1.5、3.0 m3/s。不同工況組合情況見表1。

表1 各計算工況方案Tab.2 Configuration of different water diversion schemes

表1中，工況1、2主要用于比較在相同引排水量及出入流組合方式下，密云水庫在受不同風場的影響下流場的變化情況。工況1、3、4主要用于比較在受相同風場影響的情況下，不同的流量組合對密云水庫水動力的影響變化情況。

3.2 水位及水面面積的影響

南水北調來水后對密云水庫水動力最直接的影響就是水量的增加引起了水位的增高同時水面面積增大。圖4為不同工況計算的水位隨時間變化情況，可得，工況1與工況2的水位變化線完全重合，說明風對水位變化的影響不大。表2和3分別列出了不同工況下水位、水面面積的模擬情況以及與不調水相比水位、水面面積的變化情況。與不調水(工況3)相比，調水量為10 m3/s(工況1)時水位升高了2.4 m,變幅為1.7%，水面面積增大了7.43 km2，變幅為8.73%；調水量為20 m3/s(工況4)時水位升高了4.24 m,變幅為3.04%，水面面積增大12.24 km2，變幅為14.39%。

圖4 水位模擬結果Fig.4Water level simulation result

工況水位/m最小值最大值差值水面面積/km2最小值最大值差值1135.76142.016.2568.7792.5023.732135.76142.026.2668.7792.5023.733135.76139.613.8568.7785.0716.34135.76143.858.0968.7797.3128.54

3.3 流場的影響

工況1是在風以及入流和出流的作用下形成的水體流動，而工況2是在無風的條件下，水庫內水體的流動在入流和出流的作用下，由入流和出流的慣性驅動產生[9]。由圖5和圖6工況1和工況2計算的流場分布圖可知，在調水同為10 m3/s的情況下，風場對湖區整體流場形態產生了較大的影響，在SW風的作用下整個湖區流場比無風時明顯的多，而且在局部地區形成許多渦旋。工況2的流速大部分分布在0～0.006 m/s的范圍內，工況1的流速大部分分布在0.0019～0.040 m/s的范圍內，說明風場對流速的增大發揮了明顯的作用。

表3 調水對密云水庫水位、水面面積的影響Tab.3 Effects of water diversion on water level and surfacearea of Miyun reservoir

工況2、3、4為在不同引水流量的情況下流場的分布情況，工況2、4與工況3相比，除了在水九附近形成渦旋外，在調水口附近也形成一個渦旋，而且工況4形成的渦旋比工況2的大。從局部來看，調水口受調水的影響，流速明顯增大，由原來的0.1 cm/s增加到0.46 cm/s左右。

根據上述不同工況的計算結果可得，南水北調來水后帶動了水庫水體的流動，水體流動性增大，減小了水體富營養化發生的可能性。同時水位升高，水面面積增大，外來物種的進入，可能會對密云水庫原有水生植物、動物的生存環境帶來一定的影響，接下來將做進一步的研究。

4 結 語

本文基于MIKE21FM模型建立了密云水庫的二維水動力數值模型，對密云水庫的水位變化過程進行了驗證分析。模擬結果表明，水動力模擬結果較好，誤差在計算要求的范圍內。并對模型進行了檢驗，驗證結果較好。在驗證結果較好的基礎上,分4種工況進行模擬分析，在不同的風場變化的條件下以及不同的引排水情況分析不同情景對密云水庫水動力的影響，得出的結論如下：風場對水位及水面面積的變化無影響，不同的調水量引起水位升高和水面面積增大的幅度不同；風場的變化對流場的分布及流速的影響較大，在SW風的作用下整個湖區流場比無風時明顯的多，而且在局部地區形成許多渦旋，流速明顯增大；不同的調水量對流場的影響，與不調水相比隨著調水量增大，在調水口附近形成一個渦旋，調水量增大渦旋也隨著增大，而且調水口附近流速明顯增大。

圖5 工況1流場分布(2016年3月8日)Fig.5The distribution of a flow field for the No.1 scenario (March 8, 2016)

圖6 工況2流場分布(2016年3月8日)Fig.6The distribution of a flow field for the No.2 scenario(March 8, 2016)

□

[1] COX B A.A review of currently available in-stream water-quality models and their applicability for simulating dissolved oxygen in lowland rivers[J].Sci Total Environ,2003,314/315/316(1):335-337.

[2] 譚炳卿，張國平.淮河流域水質管理模型[J].水資源保護，2001,(3)：15-18.

[3] KELLER V.Risk assessment of “down-the-drain”chemicals:search for a suitable model[J]. Sci Total Environ,2006,360(1/2/3):305-318.

[4] 陳沐松，宋蘭蘭.Mike21在垃圾發電廠溫排水數值模擬中的應用[J].安徽農業科學，2011,39(15)：9 221-9 223.

[5] Danish Hydraulic Institute (DIH).MIKE21 Flow Model:Hydrodynamic Module Scientific Documentation[M]. Denmark :Danish Hydraulic Institute, 2007.

[6] 王曉燕.非點源污染過程機理與控制管理——以北京密云水庫流域為例[M]. 北京：科學出版社，2011:22.

[7] 吳曉輝，吳 鋼，潘軻旻，等.南水北調來水對密云水庫水質和水生生物影響的預測分析[J].北京水務,2015,(6)：4-6.

[8] 李 娜，葉 閔.基于MIKE21的三峽庫區涪陵段排污口COD擴散特征模擬及對下游水質的影響[J].華北水利水電學院學報,2011，32(1):128-131.

[9] 王 輝.大伙房水庫流場影響因素的數值模擬研究[J].中國水運,2015,15(6)：76-77，81.