大型跨流域調(diào)水工程對(duì)水源地水動(dòng)力的影響研究

周紅玉，劉 操，吳曉輝

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

0 引 言

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

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

1 研究區(qū)概況

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

2 密云水庫(kù)模型的搭建

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

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

2.1 MIKE21模型簡(jiǎn)介

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

MIKE21模型的控制方程如下：

連續(xù)方程：

(1)

動(dòng)量方程：

(3)

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

2.2 二維水動(dòng)力模型計(jì)算條件設(shè)置

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

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

在密云水庫(kù)水力聯(lián)系中，潮河和白河從上游來水進(jìn)入到密云水庫(kù)，密云水庫(kù)位于兩河下游，系攔蓄兩河河水而成，模型將兩河入庫(kù)處設(shè)置為開邊界，第九水廠出水以點(diǎn)源的形式輸出，其余為閉邊界。模型采用非結(jié)構(gòu)三角網(wǎng)格劃分研究區(qū)，三角網(wǎng)格能較好的擬合研究區(qū)海陸邊界，且能隨意調(diào)整網(wǎng)格密度和網(wǎng)格大小。由密云水庫(kù)地形圖得到研究區(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個(gè)，三角網(wǎng)格的邊長(zhǎ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 水動(dòng)力模型計(jì)算條件的設(shè)置

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

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

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

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

克朗值(CFL number)：網(wǎng)格分辨率、水深和時(shí)間步長(zhǎng)決定了模型設(shè)置中的克朗值?？死手翟谛∮?的情況下模型才能正常運(yùn)行，所以為了保證模型運(yùn)行的穩(wěn)定性將CFL值設(shè)置為0.8。

干濕邊界(Flood and dry)：MIKE21 可以根據(jù)每個(gè)網(wǎng)格的水深情況調(diào)整計(jì)算條件，靈活的調(diào)用公式參與計(jì)算[8]。干濕水深是用來判斷單元網(wǎng)格是否參與到模型的計(jì)算中來。在該模型中干水深為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 水動(dòng)力模型計(jì)算結(jié)果分析與率定

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

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

2.4 水動(dòng)力模型的驗(yàn)證

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

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

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

3.1 工況設(shè)計(jì)

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

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

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

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

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

3.2 水位及水面面積的影響

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

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

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

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

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

4 結(jié) 語(yǔ)

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

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

圖6 工況2流場(chǎng)分布(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] 譚炳卿，張國(guó)平.淮河流域水質(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] 王曉燕.非點(diǎn)源污染過程機(jī)理與控制管理——以北京密云水庫(kù)流域?yàn)槔齕M]. 北京：科學(xué)出版社，2011:22.

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

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

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