周紅玉，劉 操，吳曉輝

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

0 引 言

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

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

1 研究區(qū)概況

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

2 密云水庫模型的搭建

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

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

2.1 MIKE21模型簡介

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

MIKE21模型的控制方程如下：

連續(xù)方程：

(1)

動量方程：

(3)

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

2.2 二維水動力模型計算條件設(shè)置

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

2.2.1 模型網(wǎng)格設(shè)置

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

圖1 研究區(qū)地形Fig.1 Topographic map of the study area

2.2.2 水動力模型計算條件的設(shè)置

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

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

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

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

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

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

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

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

其他沒有設(shè)置的參數(shù)采用默認(rèn)值。

2.3 水動力模型計算結(jié)果分析與率定

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

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

2.4 水動力模型的驗證

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

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

3 南水北調(diào)來水對密云水庫水文情勢的影響分析

3.1 工況設(shè)計

通過對實際調(diào)水情況的分析，本模型在綜合考慮當(dāng)?shù)貧庀髼l件，在不同的風(fēng)場變化的條件下以及不同的引排水情況分析不同情景對密云水庫水動力的影響。

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

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

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

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

3.2 水位及水面面積的影響

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

圖4 水位模擬結(jié)果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是在風(fēng)以及入流和出流的作用下形成的水體流動，而工況2是在無風(fēng)的條件下，水庫內(nèi)水體的流動在入流和出流的作用下，由入流和出流的慣性驅(qū)動產(chǎn)生[9]。由圖5和圖6工況1和工況2計算的流場分布圖可知，在調(diào)水同為10 m3/s的情況下，風(fēng)場對湖區(qū)整體流場形態(tài)產(chǎn)生了較大的影響，在SW風(fēng)的作用下整個湖區(qū)流場比無風(fēng)時明顯的多，而且在局部地區(qū)形成許多渦旋。工況2的流速大部分分布在0～0.006 m/s的范圍內(nèi)，工況1的流速大部分分布在0.0019～0.040 m/s的范圍內(nèi)，說明風(fēng)場對流速的增大發(fā)揮了明顯的作用。

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

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

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

4 結(jié) 語

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

圖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] 譚炳卿，張國平.淮河流域水質(zhì)管理模型[J].水資源保護(hù)，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在垃圾發(fā)電廠溫排水?dāng)?shù)值模擬中的應(yīng)用[J].安徽農(nóng)業(yè)科學(xué)，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]. 北京：科學(xué)出版社，2011:22.

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

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

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