秦政，雷坤，黃國鮮,3，孫明東，代丹，程全國*

1.沈陽大學環(huán)境科學與工程學院 2.中國環(huán)境科學研究院 3.青海大學三江源生態(tài)與高原農(nóng)牧業(yè)國家重點實驗室

官廳水庫是永定河流域最重要的大型控制性水庫，水庫上游社會經(jīng)濟的迅速發(fā)展使需水量日益增加，導致官廳水庫入庫水量從20世紀50年代的20.40億m3a驟降到21世紀初的1.13億m3a，近年來略有回升，達2.00億m3a，入庫水量銳減導致官廳水庫持續(xù)多年低水位運行。官廳水庫曾是北京市的重要供水水源地之一，由于水環(huán)境污染嚴重等原因，1997年被迫退出北京市生活飲用水體系，隨著北京市水資源供需矛盾的日益突出，2007年其被重新啟用為北京市的備用飲用水源地，但只保留了對北京市局部地區(qū)工業(yè)和農(nóng)業(yè)的有限供水[1-2]。為緩解官廳水庫及其下游水資源緊缺問題，自2003年起通過永定河上游水庫陸續(xù)向官廳水庫調(diào)水，在一定程度上改善了官廳水庫蓄水量不足的狀況，同時提升了水庫水質(zhì)；但官廳水庫的蓄水水位仍較低，水庫下泄流量不能滿足連續(xù)下泄和下游河道最小生態(tài)流量的需求，使永定河下游的北京—天津段仍長期處于斷流狀態(tài)[3-4]。自2016年起從黃河萬家寨調(diào)水，使流域獲得一定的附加水資源[5]，永定河上游桑干河源頭區(qū)—三家店區(qū)域都得到一定程度的補水。

為了更好地解決官廳水庫及其下游水資源供需矛盾和可持續(xù)發(fā)展，同時維持官廳水庫及其下游河道生態(tài)環(huán)境健康，迫切需要定量求解官廳水庫水資源的收支過程和平衡，開展不同調(diào)水方案和生態(tài)流量下泄條件下的湖庫水量過程和平衡問題的研究[6]。目前在官廳水庫水量平衡方面的總量研究比較缺乏，已有研究主要基于統(tǒng)計學方法，圍繞實測水量和降水等水文特性，進行長序列和年內(nèi)變化的分析，或?qū)購d水庫流域的降水、蒸發(fā)和徑流等變量的變化特點進行研究[7-8]，缺乏利用平衡方程對官廳水庫不同水文條件下水量平衡的整體分析。筆者采用2002—2016年的氣象、水文和水庫特征等數(shù)據(jù)，基于水庫水量平衡方程建立模型，求解官廳水庫蓄水量和水位的日均過程，在模型中有效考慮水庫出入流、蒸發(fā)、降水和下滲等因素的具體過程和影響，分析官廳水庫的水量平衡狀況；定量分析官廳水庫在不同收水量和生態(tài)流量下泄條件下水庫蓄水和水位變化情況，并定量分析黃河萬家寨調(diào)水能否滿足官廳水庫最小生態(tài)流量下泄的要求，以期為未來官廳水庫水資源的合理調(diào)度和改善官廳水庫下游河道的生態(tài)環(huán)境狀況提供技術支撐和決策參考。

1 研究區(qū)概況與研究方法

1.1 研究區(qū)概況

官廳水庫位于華北沿海平原向中溫帶干旱的內(nèi)蒙古高原過渡的地帶、河北省懷來縣與北京市延慶區(qū)交界處(115°34′02"E～115°49′30"E，40°13′46"N～40°25′42"N)。官廳水庫建于1954年，水庫總庫容為41.6億m3，控制流域面積為43 402 km2，正常蓄水位為479 m，死水位為471.47 m。流入官廳水庫的河流主要有洋河、桑干河和媯水河，其中洋河與桑干河在懷來縣夾河村匯流進入永定河，永定河向東8 km注入官廳水庫，而媯水河自東向西流經(jīng)延慶區(qū)后注入官廳水庫(圖1)。永定河下泄水流經(jīng)官廳水庫的攔蓄和調(diào)節(jié)，沿永定河中下游進入北京和天津段。區(qū)域氣候特征為夏短冬長，年平均氣溫為6～8 ℃，多年平均降水量為405 mm，降水主要集中在6—9月，占全年的64%～76%；庫區(qū)多年平均蒸發(fā)量為1 200～1 400 mm，水面蒸發(fā)量大于降水量[9]。

圖1 研究區(qū)及各控制站分布Fig.1 Study area and location of control stations

1.2 數(shù)據(jù)來源

根據(jù)模型計算要求，需要收集庫區(qū)降水量與蒸發(fā)量，進、出流流量，調(diào)水量和壩前下滲過程的水文數(shù)據(jù)，水庫水位-面積-庫容曲線和水庫特征水位等參數(shù)，數(shù)據(jù)來源見表1。

表1 水量平衡計算數(shù)據(jù)來源

由于官廳水庫的底泥淤積較為嚴重，淤積量達6.51億m3，而水庫缺少近年的水位-面積曲線數(shù)據(jù)，因此基于1979—2014年的長序列遙感影像水體面積解譯分析結(jié)果[13-15]，結(jié)合當時水庫對應的水位，通過分析水庫水位與對應面積的擬合曲線，求解得到水庫的水位-面積曲線，結(jié)果如圖2所示。由于所選用的分析資料大部分與本研究計算的時段(2002—2016年)重疊，因此得到的水位-面積曲線具有較好的代表性，可用于模擬時段內(nèi)水庫水量計算。官廳水庫水位-庫容關系采用1980年水位-庫容曲線(圖2)。

圖2 官廳水庫的水位-面積-庫容曲線Fig.2 Water level-area-capacity curve of Guanting Reservoir

1.3 研究方法

根據(jù)質(zhì)量守恒原理，利用官廳水庫日時間尺度數(shù)據(jù)，建立水庫水量平衡方程。公式如下：

(1)

式中：QIn(i)、QOut(j)分別為第i條入庫河流、第j個水庫出口的日平均流量，m3s；NI為入庫河流條數(shù)；NO為水庫出口數(shù)量；TsecDay為時間換算常量，表示1 d內(nèi)包含的秒數(shù)，取值86 400；P為庫區(qū)日降水量，mmd；E為庫區(qū)日水面蒸發(fā)量，mmd，由直徑20 cm蒸發(fā)皿的實測數(shù)據(jù)(E20)及折算系數(shù)(k0，取值0.60)(E=k0×E20)計算得出；A為不同水位下水庫的水面面積，m2，可根據(jù)水庫的水位-面積曲線求解得出；ΔV為水庫蓄水量的日變化，ΔV=Vt+1-Vt(t為時間，d)；WMan為人工引水水量日均過程；Sp為庫區(qū)滲透率，主要包括水庫床面滲透率(Sp1)和水庫壩體有無截滲墻條件下的水庫滲漏損失率(Sp2+Sp3)。Sp的計算公式為：

(2)

式(2)第一項為根據(jù)土壤滲流達西公式計算得到的庫底滲漏量。其中K為地基的滲透系數(shù)，md，由于官廳水庫河道及庫區(qū)底部細顆粒泥沙沉積厚度平均約為20 m，因此K取1×10-6～1×10-7cms[16]；Zs為水庫水位，m；Zsa為庫區(qū)地下水水位，m；ΔL為河床平均高程與地下水位的垂直距離[17]，m。式(2)第二項為鋪設防滲墻的土壩的滲漏率。其中B1為防滲墻的寬度，m；K0為黏土地基滲透系數(shù)，md；T為透水地基厚度，m；Le為覆蓋面材料的等效長度，m；Δh1為覆蓋面材料的末端的水頭損失，m。式(2)第三項為未鋪設防滲墻的土壩的滲漏率。其中B2為未鋪設防滲墻的土壩寬度，m；ΔH為上下游水位差，m；M1為坡度系數(shù)[18]。

計算模型主要采用Visual Studio和Inter Visual Fortran 2013組合平臺下Fortran 95語言編寫，計算時間步長為1 d，為了保證計算精度，在同一時間步長內(nèi)采用有限次(最大200次)迭代算法進行水庫水量平衡求解，水庫接收的有效降水量、庫面蒸發(fā)量和下滲量等主要與相鄰時間步長內(nèi)水庫水面面積平均值有關，因此需要迭代計算的間接變量包括有效降水量、蒸發(fā)量和下滲量，最終停止迭代計算的決策變量是相鄰迭代步內(nèi)計算得出的水位、水面面積和水體體積差足夠小。在迭代求解過程中，通過保持各出、入流流量不變條件下，迭代更新求解式(1)中每個時間步長內(nèi)水庫水面區(qū)域內(nèi)的有效降水量、蒸發(fā)量和下滲損失量，最大循環(huán)200次，在迭代過程中發(fā)現(xiàn)絕大多數(shù)條件下一般經(jīng)過3次迭代，相鄰時段內(nèi)計算的水位差絕對值就可以控制在1.0×10-8m以內(nèi)。具體計算結(jié)果表明：迭代后相鄰迭代步的水面面積差絕對值小于1.1×10-3m2，水體體積差絕對值在2.0×10-2m3以內(nèi)，每日計算產(chǎn)生的誤差較小。若不進行迭代計算，每步水位計算誤差絕對值在9×10-5m以內(nèi)，面積誤差絕對值在1.2×104m2以內(nèi)，體積誤差絕對值在7.7×104m3以內(nèi)，水位平均誤差絕對值為3.4×10-6m，面積誤差絕對值為34.7 m2，體積誤差絕對值為186.7 m3；而迭代后的水位平均誤差絕對值為5.2×10-9m，面積平均誤差絕對值為5.4×10-2m2，體積平均誤差絕對值為0.31 m3。說明從模型計算角度迭代計算是必要的，通過迭代計算能顯著降低相鄰時段內(nèi)水位、面積和體積的計算誤差，進而減小庫區(qū)水面蒸發(fā)量和下滲量的計算誤差。同時應該指出：庫區(qū)水面蒸發(fā)量和下滲量的計算受到物理過程、計算公式類型和主要參數(shù)取值等多個環(huán)節(jié)的影響[19-20]，迭代求解只是減少計算誤差的一個方面。

該模型計算所需輸入的數(shù)據(jù)包括庫區(qū)氣溫，水面蒸發(fā)量與降水量，水庫進、出口各水文站數(shù)據(jù)，水庫水位-面積-庫容曲線數(shù)據(jù)，水庫下泄的最小生態(tài)流量月過程數(shù)據(jù)，黃河萬家寨調(diào)水月流量過程，水庫各特征水位，水庫歷史實測水位和水庫滲漏流量等。

1.4 平衡因子分析

1.4.1入庫徑流量

由于永定河流入官廳水庫尾閭段的8號橋站沒有長時間觀測的水文站，所以選取位于洋河上的響水堡站和位于桑干河上的石匣里站的流量之和作為永定河入官廳水庫徑流量。洋河、桑干河、媯水河3條主要入?yún)R河流多年入庫徑流量如圖3所示。

圖3 官廳水庫各水文站2002—2016年入庫徑流量變化Fig.3 Variation of annual inflow water volume of each hydrological station of Guanting Reservoir in 2002-2016

1.4.2庫區(qū)降水量及蒸發(fā)量

官廳水庫庫區(qū)水面的降水和蒸發(fā)對水庫水量有一定程度的影響。利用延慶和沙城等氣象站的日降水量和日蒸發(fā)量數(shù)據(jù)，采用反距離加權空間插值法求得官廳水庫庫面中心點的日降水量和日蒸發(fā)量。其中上述參考氣象站的水面蒸發(fā)量采用直徑為20 cm的小型蒸發(fā)皿實測蒸發(fā)數(shù)據(jù)乘以蒸發(fā)折算系數(shù)確定[21-23]。通過統(tǒng)計得到2002—2016年官廳水庫多年水面降水量平均值為0.319億m3a，蒸發(fā)量平均值為0.535億m3a，降水量和蒸發(fā)量年變化如圖4所示。

圖4 官廳水庫庫區(qū)2002—2016年降水量與蒸發(fā)量變化Fig.4 Varation of annual precipitation and evaporation in Guanting Reservoir area in 2002-2016

1.4.3水庫下滲量

官廳水庫下滲主要分為水庫庫底下墊面的滲透、鋪設防滲墻的大壩或兩岸繞壩壩體的下滲、水庫大壩與未鋪設防滲墻的兩岸繞壩壩體的下滲。經(jīng)過多年運行，官廳水庫庫底已形成了一定厚度的淤積泥層，其主要由細沙和中沙組成，且分布較為均勻，滲透系數(shù)取1×10-7～1×10-6cms[24]；水庫攔河壩為黏土心墻壩且壩前鋪設水平的防滲蓋，因此壩前的滲漏量較??；左右肩壩壩坡經(jīng)過回填防滲土并與黏土防滲心墻連接等工程措施形成厚達5 m的防滲體，防滲效果較好；水庫左右岸繞壩未鋪設防滲體壩體，因此按照未鋪設防滲墻條件下的滲漏公式計算下滲量。

2 結(jié)果與討論

2.1 模擬結(jié)果分析

2.1.1水庫水位模擬驗證

在原有來流和調(diào)度作用下，模型計算的2002—2016年官廳水庫水位變化如圖5所示。由圖5可知，水位的模擬值與實測值相關系數(shù)為0.933 3，二者吻合良好。但在水庫低水位運行期，模擬值和實測值出現(xiàn)一定差距，主要是由于水庫嚴重淤積，模型采用的水庫水位-面積-庫容曲線與實際曲線分布可能存在一定的差異。在水庫水位較低條件下(如小于470 m)，官廳水庫的永定河干道入?yún)R口附近長期淤積形成的攔門沙將官廳水庫分成2個相對隔離的子水庫，此時實際參與調(diào)節(jié)的子水庫(主要為永定河入?yún)R子水庫)的庫容較小，在調(diào)度過程中，媯水河入?yún)R的子庫容不能自由地參與調(diào)度調(diào)節(jié)，因此整個水庫的調(diào)節(jié)功能呈削弱狀態(tài)，從而使水位的模擬值與實測值有一定差異。

圖5 官廳水庫水位模型驗證結(jié)果(2002—2016年)Fig.5 Water level verification results of Guanting Reservoir in 2002-2016 by model

2.1.2下滲量模擬驗證

根據(jù)模型計算的水庫下滲量模擬值與1998年水庫下滲量實測值對比如圖6所示。由圖6可以看出，模型計算的下滲量與水庫水位相關曲線和1998年對應的實測水位與下滲量相關點群基本吻合，因此本研究采用的下滲量算法是合理的。1998年官廳水庫下滲量隨水位增加而增加，其變化幅度為0.050～0.065 m3s，量級遠小于庫區(qū)的來流量、庫區(qū)降水量和蒸發(fā)量，因此可以判斷下滲量對庫區(qū)的水量平衡影響較小。

圖6 官廳水庫下滲量驗證(2002—2016年)Fig.6 Verification of infiltration discharge of Guanting Reservoir in 2002-2016

2.1.3現(xiàn)狀調(diào)度條件下的水庫水量收支過程

現(xiàn)有調(diào)度條件下2002—2016年官廳水庫水量平衡如表2所示。由表2可知，2002—2016年官廳水庫多年平均入庫徑流量為1.58億m3a，其中從桑干河、洋河和媯水河控制站直接入庫的多年平均徑流量為1.26億m3a，占總?cè)霂炝康?0%；而庫區(qū)多年平均降水量為0.32億m3a，占總?cè)霂炝康?0%。多年平均出庫水量為1.45億m3a，其中多年平均蒸發(fā)量為0.54億m3a，占總出庫量的37.2%；水庫多年平均下泄流量為0.9億m3a，占總出庫量的62.2%；多年平均下滲量為0.009億m3a，僅占總出庫量的0.6%，可見下滲量對水庫整體運行水位影響較小。水庫多年平均蓄變量為0.133億m3a，總體上看，2002—2016年官廳水庫屬于正收支平衡狀態(tài)，但2006年之前蓄變量多為負收支，之后多為正收支，水位呈逐漸升高態(tài)勢。

表2 2002—2016年官廳水庫水量平衡分析統(tǒng)計

雖然2002—2007年永定河流域的平均降水量為424.20 mm，屬于多年偏枯年份，但由于上游對官廳水庫的調(diào)水作用，入庫徑流量仍為水庫來水的主要部分。為保證一定的下泄流量，官廳水庫年均水位從2002年的473.48 m降至2007年的470.29 m，為多年來的最低水位。2007年后水庫開始減小下泄流量，并且上游的調(diào)水水量減小直至2010年停止調(diào)水，但由于2008—2016年永定河流域多年平均降水量增至535.26 mm，流域降水對徑流的影響具有一定的滯后且上游各水庫的蓄水延長了滯后的時間，官廳水庫入庫徑流量逐年增加，加之下泄流量的減少，官廳水庫水位逐漸抬升。

將官廳水庫各出入庫項變化作為一個整體平衡過程(來水項為正，耗水項為負，下滲量因占比較小未列入)，得出水庫平衡分析及各出入庫項在平衡中的占比，結(jié)果如圖7所示。由圖7可知，官廳水庫來水量構(gòu)成比例并不穩(wěn)定，年際變化較大。以洋河響水堡站為例，該站來水量在平衡中的占比最小，僅為7%，在水量充沛時可達31%。入庫徑流對水庫蓄水量影響較大，各站徑流量總和約占水庫來水的38%～80%。在水庫耗水項中，蒸發(fā)量的占比較為穩(wěn)定，為12%～27%；水庫下泄流量年際波動十分明顯，為7%～51%，且整體呈遞減趨勢，下泄流量在水庫耗水項中占主導地位。入庫水量的不穩(wěn)定導致水庫下泄流量年際波動較大，這給下游水生態(tài)環(huán)境帶來了一定風險，也給水庫水資源調(diào)度管理帶來一定的困難與挑戰(zhàn)[25]。

圖7 2002—2016年官廳水庫水量平衡 Fig.7 Water balances of Guanting Reservoir in 2002-2016

2.2 水庫下游生態(tài)流量計算

由于官廳水庫上游來水量存在明顯不足，加上水庫水面蒸發(fā)量較大，水庫下泄流量逐年減少，下泄流量不足將會直接影響下游的生態(tài)功能。目前官廳水庫下游河段存在嚴重的水資源不足、河道斷流和河道沙化等生態(tài)退化等問題，保障官廳水庫下游生態(tài)流量對于維護永定河生態(tài)功能健康有著重要意義。

依據(jù)Tennant法計算官廳水庫下游永定河河道生態(tài)需水量[26]。Tennant法以下泄流量與多年平均流量的之比為基準，將河流流量分為幾個等級：最大允許極限、最佳范圍、極好狀態(tài)、良好狀態(tài)、一般或較差狀態(tài)、最小狀態(tài)和極差狀態(tài)[27-28](表3)。

表3 河道不同生態(tài)環(huán)境狀況與流量對應情況

永定河流域降水年內(nèi)主要集中于6—9月，占全年降水量的76%，10月—翌年5月降水量占全年的24%。根據(jù)1956—2009年系列月天然徑流，計算出官廳水庫下游(壩下)河道的生態(tài)需水量為1.84億m3a，其中汛期為0.81億m3a，非汛期為1.03億m3a?；谠撋鷳B(tài)需水量數(shù)據(jù)，結(jié)合2002—2016年下游河道實測水量，計算官廳水庫汛期與非汛期實際生態(tài)需水量虧缺情況，結(jié)果如表4所示。由表4可知，官廳水庫下游河道生態(tài)水量虧缺較為嚴重，汛期多年平均虧缺量為0.50億m3a，虧缺率為61.73%；非汛期多年平均虧缺量為0.60億m3a，虧缺率為58.25%；官廳水庫下游河道多年平均生態(tài)水量虧缺率為60.43%。可見，造成下游河道生態(tài)水量虧缺嚴重的原因之一是官廳水庫下泄流量逐年減小。雖然多年來通過上游水庫調(diào)水以及減小下泄流量，官廳水庫蓄水水位有所升高，但下游河道的生態(tài)需水量沒有得到相應的補給。因此在水庫現(xiàn)有調(diào)度條件下難以滿足下游河道的生態(tài)需水量，需要流域外的水資源調(diào)度來補充[29-30]。通過山西省萬家寨引黃北線工程，連通桑干河與黃北干線向永定河補水(簡稱引黃調(diào)水)，是解決永定河生態(tài)水量虧缺的重要措施。引黃北線線路全長164 km，規(guī)劃引水量為5.6億m3a，可利用供水能力為2.64億m3a[31]。

表4 2002—2016年官廳水庫下游河道實測水量和生態(tài)需水量虧缺情況

2.3 不同調(diào)度方案下模型對水位的預測

利用官廳水量平衡模型，設置不同調(diào)水和生態(tài)流量下泄控制情景，模擬官廳水庫水位和下游生態(tài)流量的關系。為了檢測調(diào)水對官廳水庫正常運行水位和水庫下游生態(tài)需水量的影響，設置2種情景進行討論：1)情景1，不進行調(diào)水，保持原有水庫調(diào)度條件，只考慮滿足下游最小生態(tài)流量限制；2)情景2，考慮每年3—5月從萬家寨引黃調(diào)水1億m3a，官廳水庫收水0.5億m3a，同時滿足下游最小生態(tài)流量限制。針對以上2種情景進行模擬，同時與基于現(xiàn)狀真實水量(不考慮補水調(diào)水，也不考慮下游生態(tài)需水量)的模擬結(jié)果進行比較，如圖8所示。

圖8 2種情景下官廳水庫水位模擬結(jié)果Fig.8 Water level simulation results of two scenarios of Guanting Reservoir

由圖8可知，情景1時，官廳水庫為了滿足下游最小生態(tài)流量限制，水位將會持續(xù)下降，年均水位從2002年的473.48 m跌至2014年的466.11 m，遠小于官廳水庫正常運行水位。盡管在降水量較大，且來流充沛的情況下還能保證一定的下泄流量，但在降水量減小條件下，官廳水庫并不能滿足下游最小生態(tài)流量的需求。情景2時，官廳水庫水位保持穩(wěn)定并在降水量較大的年份略有升高，即使在沒有其他調(diào)水且在降水較少的年份也能保持一定的水位運行，既保證了官廳水庫的水位維持，也保證了下游的最小生態(tài)流量。

3 結(jié)論

(1)結(jié)合多源數(shù)據(jù)，建立官廳水庫水量日過程計算模型，編程實現(xiàn)相應算法，模型中詳細考慮了水庫河流流量入?yún)R、庫區(qū)水流輸出、降水、蒸發(fā)和下滲等過程，模型計算的2002—2016年官廳水庫庫區(qū)水位和相應實測值吻合良好。

(2)基于官廳水庫來水構(gòu)成及水庫水量平衡，分析得出官廳水庫水位受地表徑流量、庫區(qū)降水量及水庫下泄流量影響較大。庫區(qū)水面蒸發(fā)量、下滲量多年平均值分別為0.54億、0.009億m3a，且變化幅度較小，較為穩(wěn)定。水庫在沒有調(diào)水補充的情況下，需要減少下泄流量來保證水庫合適的運行水位。

(3)官廳水庫下游河道生態(tài)需水量為1.84億m3a，目前官廳水庫下游生態(tài)需水量虧缺情況較為嚴重，對比2002—2016年河道實測水量，平均生態(tài)水量虧欠率為60.43%。

(4)模型模擬結(jié)果表明，在引黃調(diào)水1億m3a，水庫收水0.5億m3a條件下，官廳水庫在降水較少的年份也能保持一定的水位運行，且保證水庫下泄流量基本滿足下游河道最小生態(tài)流量的需求。