王 玥，徐 磊，趙建世

（1.中國(guó)灌溉排水發(fā)展中心，北京 100054；2.清華大學(xué)，北京 100084）

1 研究背景

黑河發(fā)源于祁連山南部向北流至下游的綠洲居延海，潛在蒸發(fā)量為1325 mm/a，而降水量?jī)H為190 mm/a，表明該地區(qū)氣候極為干旱。含水層的地下水位從上游的約200 m到下游的約5 m不等，地表水和地下水交換頻繁［1-3］。黑河流域人口約95%聚集在農(nóng)業(yè)發(fā)達(dá)的中游地區(qū)，中游灌溉行為與下游生態(tài)系統(tǒng)的關(guān)系對(duì)整個(gè)流域的可持續(xù)發(fā)展影響重大［4-6］。

為合理利用黑河流域水資源，1992年國(guó)務(wù)院批準(zhǔn)了黑河干流（含梨園河）分水方案。為提高分水方案在實(shí)際管理中的可操作性，水利部黃河水利委員會(huì)編制了《黑河流域水量分配方案》（簡(jiǎn)稱(chēng)“97方案”），規(guī)定了不同的鶯落峽徑流量條件下的正義峽下泄水量目標(biāo)［7］。在豐、平、枯和特枯年份（對(duì)于25%、50%、75%和99%頻率），上游鶯落峽對(duì)應(yīng)徑流量為17.1、15.8、14.2和12.9億m3，下游正義峽的分水控制目標(biāo)分別為10.9億、9.5億、7.6億和6.3億m3。其它來(lái)水條件下，分配正義峽下泄水量按以上保證率水量直線(xiàn)內(nèi)插求得。在黑河分水方案實(shí)施的過(guò)程中，由于分水曲線(xiàn)自身技術(shù)特點(diǎn)，出現(xiàn)了豐水年份的分水目標(biāo)實(shí)現(xiàn)壓力大、“越豐水越難完成任務(wù)”的現(xiàn)象：在枯水年一般比較容易完成水量下泄目標(biāo)，但在豐水年下泄水量經(jīng)常達(dá)不到目標(biāo)，尤其當(dāng)出現(xiàn)來(lái)水前枯后豐的年份調(diào)度困難更大，為調(diào)度管理工作帶來(lái)了技術(shù)挑戰(zhàn)。分析這一現(xiàn)象的原因及應(yīng)對(duì)措施是黑河流域管理中亟需回答和解決的問(wèn)題。

地表-地下水模擬是分析這一問(wèn)題的重要技術(shù)手段。地表-地下水系統(tǒng)模擬的數(shù)值方法從1980年代開(kāi)始建立并迅速發(fā)展［8］，目前已經(jīng)成為流域水循環(huán)的主要手段［9-11］。在眾多模型方法中，美國(guó)地質(zhì)調(diào)查局開(kāi)發(fā)的MODFLOW為地表-地下水模擬提供了重要的平臺(tái)［12-13］，將不同類(lèi)型的地表水模擬模型和MODFLOW進(jìn)行耦合是應(yīng)用最廣泛的方法［14-15］。另外，MODFLOW也開(kāi)發(fā)了地表水模塊Stream子程序，為地表-地下水系統(tǒng)模擬提供了便捷的工具。

本文基于MODFLOW平臺(tái)構(gòu)建了黑河中游地表-地下水模擬模型，通過(guò)情景模擬分析，對(duì)黑河流域“越豐水越難完成任務(wù)”這一問(wèn)題背后的物理機(jī)制進(jìn)行了模擬和解釋?zhuān)嚓P(guān)研究結(jié)論可以為分水方案的實(shí)施提供科學(xué)依據(jù)。

2 地表-地下水?dāng)?shù)值模擬模型

本研究采用Visual MODFLOW軟件平臺(tái)構(gòu)建了黑河中游地表-地下水聯(lián)合模擬模型，其中地表水模塊采用平臺(tái)自帶的Stream子程序。模型可以用于模擬黑河干流鶯落峽至正義峽區(qū)間的地表水部分以及河道水與地下水的轉(zhuǎn)化關(guān)系，該模塊還可模擬沿程各灌區(qū)從河道取水，并計(jì)算正義峽河道流量。

2.1 Visual MODFLOW地下水模擬原理 Visual MODFLOW由加拿大Waterloo水文地質(zhì)公司基于原MODFLOW軟件開(kāi)發(fā)，是目前國(guó)際上主流的三維地下水流和溶質(zhì)運(yùn)移模擬的標(biāo)準(zhǔn)可視化專(zhuān)業(yè)軟件。這個(gè)完整的集成軟件平臺(tái)將水流模型MODFLOW、平面和剖面流線(xiàn)示蹤模型MODPATH和溶質(zhì)運(yùn)移模型MT3D同最直觀強(qiáng)大的圖形用戶(hù)界面結(jié)合在一起［16］，分為輸入、運(yùn)行和輸出3大模塊。

MODFLOW的地下水流運(yùn)動(dòng)規(guī)律描述采用有限差分的方法，首先在空間和時(shí)間上對(duì)研究區(qū)進(jìn)行離散，然后建立研究區(qū)域網(wǎng)格點(diǎn)的水平衡關(guān)系方程，將相關(guān)方程聯(lián)立形成線(xiàn)性方程組，最后采用迭代的方式進(jìn)行求解，得到每個(gè)網(wǎng)格地下水狀態(tài)量。在不考慮水的密度變化的前提下，多孔介質(zhì)中地下水在三維空間中的流動(dòng)可用偏微分方程表示為：

式中：W為單位體積流量，代表來(lái)自源或流進(jìn)匯的水量，d-1；h為水頭，m；Kxx、Kyy、Kzz分別為沿三個(gè)坐標(biāo)軸方向的水力傳導(dǎo)系數(shù)，m/d；Ss為多孔介質(zhì)貯水系數(shù)，m-1；t為時(shí)間，d。

2.2 地表水模擬模塊及其與地下水的耦合 MODFLOW的Stream模塊可以為地表-地下水循環(huán)模擬、人類(lèi)活動(dòng)影響和地下水資源評(píng)價(jià)等提供技術(shù)支撐［17-19］。本研究采用Stream子程序模擬黑河中游段鶯落峽至正義峽的河道水流過(guò)程以及地表-地下水轉(zhuǎn)化關(guān)系，同時(shí)還可以模擬中游各灌區(qū)從河道取水帶來(lái)的河道流量和地下水改變，并最終計(jì)算下游正義峽的流量。需要注意的是，Stream模塊對(duì)地表水模擬的精度不如專(zhuān)業(yè)的地表水模擬模型，未來(lái)可以將更加精確的地表水模型與MODFLOW結(jié)合，改進(jìn)模擬的精度。

Stream子程序?qū)⒑拥浪髋c地下水的交換概化為一維流動(dòng)，其中的重要參數(shù)河床水力傳導(dǎo)系數(shù)可以表示為：

式中：K為河床淤積層滲透系數(shù)；L為計(jì)算單元河段長(zhǎng)；N為河床寬度；M為河床淤積層的厚度；CRIV為河床的水力傳導(dǎo)系數(shù)［20］。

假設(shè)河道與地下含水層間的水頭損失僅受到河床淤積層的影響，且河床淤積層處于完全飽和狀態(tài)，那么河道與地下含水層交換水量可以表述為：

式中：HRIV為河道水位；CRIV為地下含水層河道間的水力傳導(dǎo)系數(shù)；hi,j,k為河流所在單元的水頭值；QRIV為地下含水層與河道之間交換量（河水補(bǔ)給地下水時(shí)為正）。當(dāng)?shù)叵潞畬铀^下降至河床淤積層之下并產(chǎn)生非飽和帶的情況下，設(shè)河床底部高程為RBOT，則河道與地下含水層的水量交換可以表達(dá)為：

從式（5）可以看出，當(dāng)河道水頭HRIV低于RBOT時(shí)，河道和地下含水層的水量交換為負(fù)值（地下水補(bǔ)給河水）。

同時(shí)，Stream子程序可模擬灌區(qū)從河道引水的量，灌區(qū)引水量可以通過(guò)設(shè)定河段起點(diǎn)與終點(diǎn)流量來(lái)代表。另外，由于Stream是耦合在MODFLOW軟件包中的一個(gè)子模塊，因此可以很好地處理地表地下水之間數(shù)據(jù)交換和相互作用，計(jì)算方便、效率高。

3 黑河中游地表-地下水模型構(gòu)建

3.1 研究區(qū)域概況 黑河流域共發(fā)育有30多條河流，可分為東部、中部、西部三個(gè)相互獨(dú)立的子水系（見(jiàn)圖1）。本次研究區(qū)為東部的黑河干流中游（含梨園河及沿山20余條小支流，如圖1中黑框部分所示），研究區(qū)域總面積為8717km2。中部子水系包括馬營(yíng)河等小河流，都是淺山短流，流量小且出山后很快被消耗盡，對(duì)研究區(qū)水循環(huán)無(wú)影響；西部子水系為討賴(lài)河流域，包括討賴(lài)河和洪水河等，匯合后流入下游金塔盆地，對(duì)研究區(qū)無(wú)影響。

研究區(qū)地質(zhì)構(gòu)造復(fù)雜，邊緣被斷裂帶控制，覆蓋第四紀(jì)松散地層可存儲(chǔ)大量地下水，基底巖性是第三系泥巖；山前堆積物主要是河流沖洪積物，含水層以砂礫石為主，向北巖性逐漸過(guò)渡到細(xì)土平原區(qū)含亞砂土和亞黏土的含水層，含水層逐漸變淺、滲透性減弱。山前南半部的洪積扇最上層為大厚度單層潛水區(qū)，下游為細(xì)土平原區(qū)，逐步由單一潛水含水層區(qū)變?yōu)槎鄬拥臐撍?承壓水層混合區(qū)［21-22］。黑河流域地表地下水交換頻繁，出山口流出的地表徑流在南部山前大量滲漏補(bǔ)給地下，地下水到了張掖盆地中北部的細(xì)土平原帶，因含水層滲透性減弱并且地表河床下切，導(dǎo)致地下水大量溢出地表形成地表徑流和泉水，其中泉水溢出帶主要分布在洪積扇緣和細(xì)土平原地，而中游用水灌溉后又有水量補(bǔ)給到地下含水層，形成多次地表-地下交換［21-22］。山前扇群帶地下水埋深變化劇烈，扇頂為200～500 m，近山側(cè)為100～300 m。扇中張掖盆地地下水埋深約在50～100 m之間，酒泉東盆地地下水埋深在100～250 m之間。扇緣張掖盆地地下水埋深約為10～20 m，而酒泉東盆地則在80～200 m之間。細(xì)土平原區(qū)地下水埋深變化也較大，其中張掖盆地大多在5 m以?xún)?nèi)，而酒泉東盆地南半部則高達(dá)10～50 m，而北半部只有1～5 m。

黑河干流水系多年平均出山口徑流量為25.11億m3，其中鶯落峽站多年平均徑流量為16.19億m3。黑河流域礦化度小于2 g/L的多年平均地下水資源量為21.76億m3，其中山丘區(qū)為9.36億m3，平原區(qū)為20.32億m3，二者重復(fù)量為7.92億m3。全流域共有中小型水庫(kù)57座，總庫(kù)容為2.76億m3；干流灌區(qū)引水工程96處，設(shè)計(jì)引水能力約269.8 m3/s，其中直接從干流引水的口門(mén)有36處，設(shè)計(jì)引水能力為223.8 m3/s；機(jī)電井11 076眼，其中配套機(jī)井9771眼，年提水量5.81億m3。流域灌溉面積為465.67萬(wàn)畝，其中林草灌溉面積為95.91萬(wàn)畝，農(nóng)田灌溉面積為369.75萬(wàn)畝，30萬(wàn)畝以上的大灌區(qū)有8處共計(jì)301.44萬(wàn)畝。

3.2 黑河中游地表-地下水模型構(gòu)建

3.2.1 模擬范圍及含水層條件 本文研究的空間范圍東起山丹縣祁家店水庫(kù)、西至酒泉東盆地清水車(chē)站，南北以山前為界。由于承壓區(qū)沒(méi)有穩(wěn)定的隔水層，研究區(qū)各含水層有不可分割的水力聯(lián)系，是一個(gè)較為完整的地下水盆地（如圖1所示）。近幾十年來(lái)，由于地下水的開(kāi)采量不斷增加，使得承壓含水層的連通性增加，不同層之間的水力聯(lián)系較歷史上更為緊密，導(dǎo)致上下層水位相差很小甚至基本相同?；诖?，本文將研究區(qū)含水層概化為非均質(zhì)的各項(xiàng)同性單層潛水層，同時(shí)由于區(qū)內(nèi)潛水含水層分布廣、厚度大，將研究區(qū)地下水運(yùn)移視為二維非穩(wěn)定流。

圖1 黑河流域水系及中游研究區(qū)

研究區(qū)模型采用正交網(wǎng)格剖分，剖分格距為ΔX=ΔY=1 km，利用GIS中大地坐標(biāo)的底圖導(dǎo)入MODFLOW圈定模擬區(qū)域范圍。東西向剖分為172（km）格，南北向剖分為155（km）格，垂向上為單層。河流寬度和深度由模型根據(jù)地形和流量計(jì)算。

3.2.2 模型邊界條件 模型邊界條件包括側(cè)向邊界和垂向邊界，具體設(shè)置如下。

側(cè)向邊界條件：研究區(qū)域周邊二類(lèi)流量邊界，其中南、東、北邊界均為弱透水邊界，有較弱的側(cè)向補(bǔ)給［22］：南邊界山前斷裂帶為基巖裂隙水和溝谷潛流補(bǔ)給，而東部民樂(lè)和山丹斷面以及北部山前斷面均為區(qū)外側(cè)向流入邊界。西邊界由于是酒泉西盆地與酒泉東盆地分水嶺，可視為零通量邊界。

垂向邊界條件：上邊界是天然潛水面，下邊界是侏羅系和第三系泥巖以及砂質(zhì)泥巖等隔水邊界，因此可忽略頂托越流補(bǔ)給。

3.2.3 模擬時(shí)段與參數(shù) 模型參數(shù)率定和驗(yàn)證的時(shí)段為1995年1月—2000年12月，計(jì)算步長(zhǎng)為月，共計(jì)72個(gè)時(shí)段。本文采用ΔX=ΔY=1 km格距的正交網(wǎng)格進(jìn)行研究區(qū)的網(wǎng)格剖分，河床寬度N數(shù)據(jù)采用代表性斷面的實(shí)測(cè)平均寬度。模擬區(qū)域總面積為8717 km2，東西方向可剖分為172格，南北方向可剖分為155格，垂直方向?yàn)閱螌?。跟?jù)《黑河流域水資源評(píng)價(jià)利用及分配方案研究》中提供的不同巖性土的給水度，以及“九五”和“十五”期間黑河中游地區(qū)的研究資料提供的水文地質(zhì)參數(shù)，確定研究區(qū)導(dǎo)水系數(shù)大約在500～14000 m3/d之間，而給水度的范圍大約在0.08～0.2之間。本文利用Stream子程序包模擬黑河干流中游段河水與地下水的頻繁轉(zhuǎn)化關(guān)系，河床底板高程和河床厚度及水力傳導(dǎo)系數(shù)等相關(guān)參數(shù)參考文獻(xiàn)［23-25］，并通過(guò)模型率定與驗(yàn)證進(jìn)行了調(diào)整，鶯落峽的實(shí)測(cè)流量和各灌區(qū)實(shí)測(cè)引水流量來(lái)自黑河流域水資源管理局統(tǒng)計(jì)資料。1995年初的地下水初始水位根據(jù)地下水觀測(cè)井?dāng)?shù)據(jù)空間插值得到，經(jīng)過(guò)模型率定過(guò)程將模型模擬的地下水位與觀測(cè)井實(shí)測(cè)值對(duì)比，2000年末研究區(qū)地下水埋深與觀測(cè)值擬合良好，以此作為模型計(jì)算分析時(shí)段2001年初的初始水位值和初始流場(chǎng)。

3.2.4 模型率定與驗(yàn)證 根據(jù)上述資料給定研究區(qū)的初始參數(shù)，經(jīng)過(guò)多次調(diào)試和優(yōu)選，以各觀測(cè)孔處的模擬水頭值與實(shí)測(cè)值最接近為原則來(lái)率定模型參數(shù)。模型采用1995—2000年間的實(shí)測(cè)資料（包括氣象水文、地下水位觀測(cè)資料和水資源利用數(shù)據(jù)等）對(duì)地表-地下水模型進(jìn)行參數(shù)率定和模型驗(yàn)證。注意的對(duì)比數(shù)據(jù)包括：（1）研究區(qū)17個(gè)觀測(cè)孔水位動(dòng)態(tài)過(guò)程；（2）下游正義峽斷面的流量過(guò)程。

統(tǒng)計(jì)結(jié)果表明，模型率定期（1995—1997年），觀測(cè)井水位值模擬誤差在1.0m以?xún)?nèi)的數(shù)據(jù)占98%（其中誤差小于0.5m的占61%，誤差在0.5～1.0m之間的占37%），Nash-Sutcliffe效率系數(shù)為0.87；在模型驗(yàn)證期（1998-2000年），觀測(cè)井水位值模擬誤差在1.0m以?xún)?nèi)的數(shù)據(jù)占96%（其中誤差小于0.5m的占51%，誤差在0.5～1.0m之間的占45%），Nash-Sutcliffe效率系數(shù)為0.83，表明模型參數(shù)率定和驗(yàn)證效果良好。選取其中兩個(gè)具有代表性的觀測(cè)井率定和驗(yàn)證期內(nèi)的擬合過(guò)程如圖2所示，率定和驗(yàn)證期末（2000年12月末）的地下水位等高線(xiàn)以及所有觀測(cè)井的對(duì)比情況如圖3所示。

圖2 47號(hào)觀測(cè)井和65號(hào)觀測(cè)井的擬合結(jié)果

黑河下游正義峽斷面的流量過(guò)程模擬值與實(shí)測(cè)值對(duì)比結(jié)果如圖3所示。由于實(shí)測(cè)數(shù)據(jù)只有觀測(cè)井的水位值，圖3（a）只給出了模擬流場(chǎng)，無(wú)法給出模擬流場(chǎng)和實(shí)測(cè)流場(chǎng)的對(duì)比；圖3（b）給出了率定和驗(yàn)證期末觀測(cè)井模擬和實(shí)測(cè)水位的對(duì)比，表明觀測(cè)井模擬效果良好。圖4給出了下游正義峽斷面的模擬流量和實(shí)測(cè)流量對(duì)比，其中模擬流域?yàn)镾tream模塊輸出值，模型率定期（1995—1997年）的Nash-Sutcliffe效率系數(shù)為0.82，相關(guān)系數(shù)為0.88；模型驗(yàn)證期（1998—2000）的Nash-Sutcliffe效率系數(shù)為0.84，相關(guān)系數(shù)為0.89，表明模型對(duì)地表和地下水交換作用的模型效果整體上較好。

圖3 率定和驗(yàn)證期末地下水位等高線(xiàn)與觀測(cè)井水位對(duì)比

圖4 正義峽斷面流量過(guò)程擬合結(jié)果對(duì)比

4 地表-地下水交互作用及其影響分析

為了分析分水方案調(diào)度的情況，利用上述模型模擬了黑河統(tǒng)一調(diào)度的實(shí)際情景，時(shí)段為2001年1月至2006年12月。由于上述模型率定驗(yàn)證的結(jié)束時(shí)段為2000年12月，此時(shí)黑河中游的地下水埋深與觀測(cè)值擬合良好，可以作為模型模擬的初始水位值。

2001—2006年的模擬結(jié)果如表1和圖5所示。從結(jié)果可以看出，“鶯落峽來(lái)水越豐時(shí)越難完成分水指標(biāo)”這一現(xiàn)象確實(shí)存在。2001—2006年的模擬結(jié)果和1992年國(guó)務(wù)院批準(zhǔn)了黑河干流（含梨園河）分水方案中列出的泄水指標(biāo)分別相差-0.22億、-0.90億、-1.76億、-0.91億、-1.59億、-1.17億m3，在來(lái)水較豐的2003、2005和2006年，差異都超過(guò)了1.0億m3。而產(chǎn)生這一現(xiàn)象的主要原因，是由于豐水年河道水位相對(duì)較高、地下水向的補(bǔ)給量顯著減少，造成“越豐水越難完成分水目標(biāo)”現(xiàn)象。對(duì)比2003年和2001年，地下水向河道的凈補(bǔ)給量分別為1.05億和2.98億m3，豐水的2003年相對(duì)于枯水的2001年減少了1.93億m3，在這種情況下即使壓縮中游用水也無(wú)法完成正義峽分水指標(biāo)。這一現(xiàn)象體現(xiàn)了地下水作為一個(gè)天然調(diào)蓄池對(duì)水資源的調(diào)豐補(bǔ)枯作用，對(duì)于地表水的調(diào)度和管理具有不可忽視的影響。

表1 分水情景模擬 （單位：108m3）

圖5 分水情景對(duì)比分析

基于上述分析，建議對(duì)現(xiàn)行的1992年分水方案進(jìn)行修正，因?yàn)檫@一方案基本上按照線(xiàn)性關(guān)系來(lái)考慮不同豐枯年份水量分配，沒(méi)有考慮的地下水調(diào)蓄作用對(duì)實(shí)際調(diào)度過(guò)程的影響。如果將這一因素考慮在內(nèi)，建議將豐水年向正義峽的下泄比例降低5%～10%、平水年的下泄比例保持不變、枯水年的下泄比例增加5% ～10%，這樣即有利于方案的實(shí)際執(zhí)行，也有利于中游的用水保障和下游的生態(tài)保護(hù)。具體的調(diào)整方案，建議綜合考慮各方面因素進(jìn)行詳細(xì)的論證。另外，如果可以在上游修建控制性水庫(kù)，進(jìn)行多年豐枯調(diào)節(jié)，將更有利于整個(gè)分水方案的執(zhí)行和流域的可持續(xù)發(fā)展。

5 結(jié)語(yǔ)

本研究構(gòu)建了黑河中游地表-地下水模擬模型，利用上述模型模擬了2001年1月—2006年12月黑河統(tǒng)一調(diào)度的實(shí)際情景，分析了地表-地下水交互作用對(duì)水量調(diào)度的影響。結(jié)果表明，向正義峽下泄的實(shí)際水量與1992年分水方案指標(biāo)相比，在來(lái)水較豐的2003、2005和2006年，差異均超過(guò)了1.0億m3，產(chǎn)生了“越豐水越難完成分水目標(biāo)”的現(xiàn)象。豐水的2003年和枯水的2001年對(duì)比表明，由于豐水年份河道水位相對(duì)較高，地下水向河道補(bǔ)給量會(huì)有所減少了1.93億m3，進(jìn)而導(dǎo)致上游向正義峽的下泄比例減少，難以完成分水指標(biāo)。因此，建議對(duì)分水方案進(jìn)行調(diào)整，將豐水年向正義峽的下泄比例降低5%～10%、平水年的下泄比例保持不變、枯水年的下泄比例增加5%～10%，同時(shí)在有可能的情況下建設(shè)上游多年調(diào)節(jié)水庫(kù)，促進(jìn)分水方案的科學(xué)實(shí)施和流域水資源的高效利用。

鑒于調(diào)度問(wèn)題的復(fù)雜性，后續(xù)研究可以在本文地表-地下水模型的基礎(chǔ)上，對(duì)不同水文條件和經(jīng)濟(jì)社會(huì)發(fā)展情景下的分水方案進(jìn)行綜合對(duì)比分析，進(jìn)一步提出合理分水指標(biāo)調(diào)整方案和調(diào)度管理實(shí)施方案。