王志勇，蔣夢(mèng)瑤

(1.黃河水利委員會(huì)河南水文水資源局，河南 鄭州 450000;2.蘭州大學(xué)資源環(huán)境學(xué)院，甘肅 蘭州 730000)

E-mail：2072802457@qq.com

20世紀(jì)以來(lái)隨著世界經(jīng)濟(jì)的快速發(fā)展，許多干旱地區(qū)地下水開(kāi)采量越來(lái)越大，以滿(mǎn)足日益增長(zhǎng)的用水需求。2012年12月12日，國(guó)務(wù)院發(fā)布并實(shí)施了最嚴(yán)格水資源管理制度[1]。數(shù)值模擬作為地下水資源研究領(lǐng)域中，定量分析水量、水質(zhì)不可或缺的手段，被廣泛應(yīng)用于最嚴(yán)格水資源管理制度的落實(shí)和預(yù)測(cè)[2-4]。本文運(yùn)用目前國(guó)際上應(yīng)用最廣泛的三維地下水模擬系統(tǒng)軟件GMS(Groundwater Modeling System)，建立了焉耆盆地地下水?dāng)?shù)值模型，并利用15個(gè)觀(guān)測(cè)井的地下水位監(jiān)測(cè)數(shù)據(jù)及遙感反演數(shù)據(jù)[5]進(jìn)行驗(yàn)證和校準(zhǔn)，分析了現(xiàn)狀地下水均衡狀態(tài)。最后，采用減少地下水開(kāi)采量、提高地表水引水量的處理方案，設(shè)計(jì)了5種不同的情景模式，并依據(jù)數(shù)值模擬結(jié)果提出焉耆盆地灌溉農(nóng)業(yè)的可持續(xù)開(kāi)采量。

1 研究區(qū)概況

焉耆盆地位于中國(guó)西北部的新疆，如圖1所示，地理坐標(biāo)為41°23′—43°31′N(xiāo)，86°39′—88°20′E，隸屬焉耆縣、和靜縣、和碩縣、博湖縣以及8個(gè)兵團(tuán)[9]，平均高程約為1050～1200m，占地總面積為1.32×106ha2，盆地最洼處為中國(guó)最大的內(nèi)陸淡水湖—博斯騰湖，面積約為1.0×105ha2[6]。

作為典型的中溫帶大陸性氣候區(qū)，年蒸發(fā)量約為2000～2450mm，而年降雨量只有50～80mm。開(kāi)都河是盆地內(nèi)最大的河流，也是最主要的地表水補(bǔ)給來(lái)源，多年平均徑流量為36.4×108m3[7]。焉耆盆地地下水開(kāi)采量的85.08%用于農(nóng)業(yè)灌溉[8]，且用水效益低，屬于典型的干旱綠洲灌區(qū)。隨著高效節(jié)水技術(shù)的深入推廣，焉耆盆地高效節(jié)水面積已從2000年的160ha2上升為2016年的1.47×105ha2，其中高效節(jié)水灌區(qū)每增加10000ha2，地下水位降低0.25m[9]，盆地內(nèi)地下水年開(kāi)采量達(dá)6.76×108m3，處于過(guò)度開(kāi)發(fā)狀態(tài)[10]。

焉耆盆地是一個(gè)完整的地下水系統(tǒng)，四周環(huán)山，西北、北部山地?cái)鄩K隆起強(qiáng)烈，形成高山區(qū)。盆地內(nèi)可大致劃分為三大地貌單元：洪積一沖積傾斜平原、開(kāi)都河三角洲以及博斯騰湖湖積平原[11]。

2 建立MODFLOW模型

本文基于GMS軟件建立干旱綠洲灌區(qū)地下水流數(shù)值模擬模型。模擬面積為10880km2，包括了盆地內(nèi)的5個(gè)灌區(qū)和博斯騰湖泊。

2.1 水文地質(zhì)概念模型

2.1.1含水層概化

根據(jù)地質(zhì)鉆孔資料、第四系巖性結(jié)構(gòu)特征以及含水層埋藏條件，將研究區(qū)內(nèi)模型模擬深度定為150～200m，模擬深度以上的含水介質(zhì)分為四層，區(qū)內(nèi)地下水主要開(kāi)采層在第四層。地下水流系統(tǒng)通過(guò)模型上邊界，即潛水含水層自由水面與外部發(fā)生垂向水量交換，模型的底部為第四系松散堆積物的下邊界。含水層巖性呈現(xiàn)非均勻的變化特征，地下水水力坡度較大，滲流速度隨著方向不同而發(fā)生變化，且具有垂直分量；地下水流的各運(yùn)動(dòng)要素隨時(shí)間而發(fā)生變化。綜上所述，研究區(qū)概化為四層非均質(zhì)各向異性的三維非穩(wěn)定流數(shù)值模型。

圖1 研究區(qū)地理位置圖

2.1.2邊界條件和源匯項(xiàng)處理

依據(jù)模擬區(qū)內(nèi)的水文地質(zhì)條件資料，模型底邊界為隔水邊界；上邊界參與入滲和蒸發(fā)過(guò)程，概化為開(kāi)放邊界。補(bǔ)給方式主要有山前側(cè)向補(bǔ)給、河道滲漏、渠系滲漏、灌溉入滲和少量降水入滲。主要排泄方式有農(nóng)業(yè)灌溉開(kāi)采、蒸發(fā)蒸騰、排水渠排水、河流排泄和向博斯騰湖排泄。區(qū)內(nèi)的北部和西部為接受山區(qū)側(cè)向徑流補(bǔ)給的山前地帶，其邊界概化為變流量邊界(GHB)；東部和南部以被看作等勢(shì)體的博斯騰湖為界，將其邊界概化為已知水頭邊界，其水頭值由博斯特湖實(shí)測(cè)水位指定，如圖2所示。

圖2 模型側(cè)邊界條件示意圖

2.2 數(shù)學(xué)模型

根據(jù)上述地下水流場(chǎng)、地下水運(yùn)動(dòng)定律及水文地質(zhì)條件等的系統(tǒng)分析，本文所建立的四層結(jié)構(gòu)的非均質(zhì)三維非穩(wěn)定流數(shù)值模型的數(shù)學(xué)表達(dá)式見(jiàn)式(1)。

(1)

式中，Kxx、Kyy、Kzz—X、Y、Z方向的滲透系數(shù)，m/d；Kxx=Kyy；H—水頭值，m；H0—初始水位，m；ε—源匯項(xiàng)，1/d；S—給水度，1/m；潛水含水層指重力給水度μd，第四層含水層在發(fā)生地下水開(kāi)采活動(dòng)時(shí)，會(huì)出現(xiàn)由承壓轉(zhuǎn)無(wú)壓的狀況，當(dāng)呈現(xiàn)承壓含水層時(shí)，模型計(jì)算使用彈性釋水系數(shù)，而含水層在無(wú)壓狀態(tài)下時(shí)，模型計(jì)算使用重力給水度，其它含水層則使用彈性釋水系數(shù)μe；Ω—模型模擬區(qū)域的范圍；n——邊界面的外法線(xiàn)方向；Г—側(cè)邊界；B—底部邊界。

2.3 地下水?dāng)?shù)值模型

根據(jù)研究區(qū)現(xiàn)有資料將模擬時(shí)間定為2000—2016年，劃分為17個(gè)應(yīng)力期，時(shí)間步長(zhǎng)為1天。在模型區(qū)域空間離散化的基礎(chǔ)上，采用GMS中1000m×1000m的三維網(wǎng)格將模擬區(qū)劃分為80個(gè)東西向和136個(gè)南北向的網(wǎng)格，如圖2所示，上下分共4層，共計(jì)有效計(jì)算單元數(shù)43520個(gè)。求解數(shù)值模型的過(guò)程中，根據(jù)研究區(qū)相關(guān)勘察報(bào)告、水文地質(zhì)條件和已有的抽水試驗(yàn)數(shù)據(jù)賦值模型中各分區(qū)參數(shù)的初值。在瞬態(tài)模型中，初始水頭等于2000年1月穩(wěn)定流模型生成的解。

2.4 模型識(shí)別與驗(yàn)證

以2001年1月1日到2007年12月31日作為模型參數(shù)的率定期，通過(guò)參數(shù)區(qū)劃和試點(diǎn)法相結(jié)合的方法，依據(jù)不同時(shí)段的水位實(shí)測(cè)值對(duì)模型含水層參數(shù)反復(fù)調(diào)整，見(jiàn)表1，得到含水層的滲透系數(shù)范圍為0～90m/d，給水度范圍為0.1～0.25。以2007年1月1日到2016年12月31日作為模型參數(shù)的驗(yàn)證期。本文選取圖2所示的15個(gè)觀(guān)測(cè)井進(jìn)行參數(shù)的率定及驗(yàn)證。選用3眼典型的長(zhǎng)觀(guān)井，對(duì)比分析其模擬水位與實(shí)測(cè)水位，如圖3所示，各個(gè)觀(guān)測(cè)點(diǎn)地下水位模擬值與實(shí)測(cè)值擬合程度較高，相關(guān)系數(shù)為0.89，水位動(dòng)態(tài)變化趨勢(shì)基本一致。且依據(jù)相關(guān)勘察報(bào)告可知，自2000年以后，地下水位不斷下降，降幅度在2～3m之間，部分地區(qū)水位下降超過(guò)5m，與模擬結(jié)果相符。可見(jiàn)所建模型能較為真實(shí)地反映焉耆盆地地下含水層的實(shí)際情況，可以用來(lái)進(jìn)行預(yù)測(cè)分析。

表1 含水層特性及水文地質(zhì)參數(shù)取值

圖3 部分觀(guān)測(cè)井水位動(dòng)態(tài)擬合圖

圖4 實(shí)測(cè)水位與模擬將結(jié)果擬合圖

3 地下水開(kāi)采方案及預(yù)測(cè)

3.1 地下水量均衡分析

灌區(qū)地下水的過(guò)度開(kāi)采將會(huì)對(duì)農(nóng)業(yè)生產(chǎn)和生態(tài)產(chǎn)生不利影響，破壞水資源平衡。在保持其他源匯項(xiàng)不變條件下，以現(xiàn)狀年(2016年)6.75×108m3/a的開(kāi)采量運(yùn)行MODFLOW模型至2050年，得到地下水變化均衡表，見(jiàn)表2。

3.1.1現(xiàn)狀水均衡

區(qū)內(nèi)模擬所得現(xiàn)狀年地下水總補(bǔ)給量為143658×104m3/a，入滲補(bǔ)給作為主要補(bǔ)給來(lái)源的補(bǔ)給量是77862.6×104m3/a，其他補(bǔ)給來(lái)源為河流滲漏補(bǔ)給、山前側(cè)向補(bǔ)給、和湖泊補(bǔ)給，補(bǔ)給量分別為16664.3、46114.2和3016.8×104m3/a；區(qū)內(nèi)地下水排泄量為144466×104m3/a，其中地下水開(kāi)采量為67566.0×104m3/a，占總排泄量的46.2%，排水渠排泄量為7926.0×104m3/a，河流排泄量為12535.0×104m3/a，蒸發(fā)蒸騰量為56439.0×104m3/a。該數(shù)值模擬結(jié)果與已有的焉耆盆相關(guān)工作報(bào)告中均衡法計(jì)算的水均衡結(jié)果和補(bǔ)排關(guān)系相一致。

3.1.2水均衡量變化

運(yùn)行到2050年，發(fā)現(xiàn)盆地內(nèi)的博斯騰湖北部和開(kāi)都河南部部分區(qū)域出現(xiàn)疏干現(xiàn)象，地下水位不斷下降，地下水儲(chǔ)量也減少。這也是為什么2025年的地下水開(kāi)采量還沒(méi)有發(fā)生變化，但到2050年已經(jīng)減少到了64258.0×104m3；蒸散發(fā)量增加至2025年后也開(kāi)始減小。

由此表明，焉耆盆地地下水處于嚴(yán)重超采的狀態(tài)。為防止研究區(qū)內(nèi)地下水位進(jìn)一步下降，減少地下水開(kāi)采量是最有效的手段。

3.2 開(kāi)采情景下的地下水模擬

考慮高效節(jié)水灌溉工程的實(shí)施、控制地下水超采和用水紅線(xiàn)指標(biāo)[12]3個(gè)方面，通過(guò)減少地下水開(kāi)采，結(jié)合地表水置換地下水的原理設(shè)置了5個(gè)開(kāi)采方案，各開(kāi)采方案的年開(kāi)采量分別為6.38×108、5.35×108、4.52×108、4.29×108、4.12×108和3.83×108m3，并以2016年為現(xiàn)狀年，將該方案及現(xiàn)狀年的開(kāi)采規(guī)模輸入識(shí)別后的模型，在保持其他源匯項(xiàng)不變的情況下，根據(jù)各方案開(kāi)采量?jī)H對(duì)Well和Recharge參數(shù)作相應(yīng)修改，運(yùn)行數(shù)值模擬至規(guī)劃年(2050年)，得到各開(kāi)采方案下的地下水流場(chǎng)圖，如圖5所示。

表2 焉耆盆地地下水均衡表 單位：104m3

圖5 2050年各開(kāi)采方案下地下水流場(chǎng)圖

本文將焉耆盆地合理地下生態(tài)水位埋深下限設(shè)為4.5m，上限埋深定為3.2m[13]。按照方案1的開(kāi)采規(guī)模運(yùn)行到2050年，有超過(guò)65%的區(qū)域地下水埋深大于生態(tài)水位下限；方案2和方案3在規(guī)劃年的地下水埋深有較為明顯的上升，但只有超過(guò)45%左右的區(qū)域地下水埋深介于合理的生態(tài)水位范圍內(nèi)。方案4和方案5在規(guī)劃年的地下水埋深60%以上都在合理生態(tài)水位之間，且與初始地下水流向及流態(tài)基本保持一致，地下水系統(tǒng)的流入流出量基本均衡。但考慮到國(guó)務(wù)院下達(dá)的最嚴(yán)格用水制度對(duì)地下水和地表水利用指標(biāo)的要求，認(rèn)為在有限的水資源情況下應(yīng)最大限度的利用地表和地下水資源，方案4開(kāi)采規(guī)模為4.12×108m3/a時(shí)，地表水引水量達(dá)12.63×108m3，下水位不會(huì)再持續(xù)降低，且逐漸恢復(fù)穩(wěn)定，滿(mǎn)足焉耆盆地用水紅線(xiàn)指標(biāo)。

4 結(jié)語(yǔ)

(1)本文所建立的焉耆盆地的地下水?dāng)?shù)值模擬模型可真實(shí)反映當(dāng)?shù)氐膶?shí)際情況，預(yù)測(cè)的結(jié)果具有較強(qiáng)的參照性。其中，滲透系數(shù)范圍為0～90m/d，給水度范圍為0.1～0.25，典型觀(guān)測(cè)孔的實(shí)測(cè)與模擬水位擬合的相關(guān)系數(shù)為0.89。

(2)方案4的地下水開(kāi)采情景運(yùn)行模型至2050年，研究區(qū)內(nèi)的地表水引水量達(dá)12.63×108m3，60%以上的地下水埋深都在合理的生態(tài)水位之間，確定焉耆盆地符合灌溉農(nóng)業(yè)可持續(xù)發(fā)展需求的地下水可持續(xù)開(kāi)采量為4.12×108m3/a。

(3)本文主要對(duì)盆地內(nèi)不同開(kāi)采規(guī)模下的地下水位進(jìn)行了預(yù)測(cè)，在往后的科研中應(yīng)綜合、系統(tǒng)的研究地下水開(kāi)采量和分布對(duì)各補(bǔ)排項(xiàng)及水位的影響。