基于MODFLOW的撓力河流域地下水承載力綜合評價(jià)

郗鴻峰,劉春偉,杜新強(qiáng),劉小天,馬河寬,于令芹

(1.山東省地質(zhì)礦產(chǎn)勘查開發(fā)局八〇一水文地質(zhì)工程地質(zhì)大隊(duì),山東 濟(jì)南 250000;2.吉林大學(xué) 新能源與環(huán)境學(xué)院,吉林 長春 130000)

目前,地下水承載力綜合評價(jià)在全國各個(gè)地區(qū)得到了廣泛應(yīng)用,各學(xué)者也構(gòu)建了適合不同區(qū)域的評價(jià)指標(biāo)體系及評價(jià)等級;在東北農(nóng)灌區(qū)應(yīng)用較廣的8個(gè)指標(biāo)為地下水資源耕地灌溉率、地下水開發(fā)利用程度、地下水利用率、供水模數(shù)、需水模數(shù)、重復(fù)利用率、人均供水量、生態(tài)環(huán)境用水率[1-2]。這一套指標(biāo)體系已經(jīng)應(yīng)用到撓力河流域,并采用主成分分析法、集對分析法進(jìn)行計(jì)算評價(jià)[3]。但諸多實(shí)踐表明,承載力綜合評價(jià)所得到的結(jié)果不盡合理,且地下水承載力綜合評價(jià)結(jié)果的合理性尚未得到有效驗(yàn)證,由此就出現(xiàn)了綜合評價(jià)指標(biāo)體系及評價(jià)等級合理性的問題[4-5]。

地下水承載力綜合評價(jià)是否合理,目前主要從評價(jià)指標(biāo)的選取、指標(biāo)體系的構(gòu)建、評價(jià)方法的選擇等方面控制[5],但從地下水承載力綜合評價(jià)結(jié)果合理性角度,目前學(xué)者只是對評價(jià)結(jié)論進(jìn)行分析,發(fā)表的相關(guān)研究論文較少。評價(jià)結(jié)果的合理性,直接影響地下水資源承載能力真實(shí)性;因此,本文以水文實(shí)測數(shù)據(jù)為基礎(chǔ),利用較成熟的地下水?dāng)?shù)值模擬軟件MODFLOW,對撓力河流域農(nóng)灌區(qū)地下水進(jìn)行建模,評價(jià)各灌區(qū)開采潛力,并以所得結(jié)果作為參考,分析地下水承載力綜合評價(jià)結(jié)果的合理性,以期促進(jìn)本地區(qū)相關(guān)指標(biāo)體系的構(gòu)建與優(yōu)化。

1 研究區(qū)概況

撓力河流域由撓力河不斷遷移和泛濫形成,位于黑龍江省東北區(qū)域[3],流域總面積約2.5萬 km2,南部包括寶清縣行政區(qū)全部,西、北部包括友誼、富錦、集賢、雙鴨山等縣市行政區(qū)部分區(qū)域,東部包括饒河縣的部分區(qū)域[3];并建有紅興隆和建三江2個(gè)農(nóng)墾管理局,流域包括東泄總灌區(qū)、錦西灌區(qū)、八五三灌區(qū)、大興灌區(qū)、龍頭橋灌區(qū)、錦南灌區(qū)、友誼西部灌區(qū)等12個(gè)灌區(qū)及其他灌區(qū)的部分區(qū)域[3]。

該流域地處中緯度,為大陸季風(fēng)氣候,雨熱同期,年均氣溫1℃～4℃,年降水量500～600 mm,集中在6-9月。流域內(nèi)優(yōu)勢河流主要有撓力河和七星河;撓力河全長約596 km,流域總面積24 863 km2,其中撓力河干流流域面積為13 084 km2,支流七星河流域面積為10 815 km2,其他流域面積964 km2;撓力河流域整體為條狀,支流基本呈現(xiàn)葉狀。該流域因地形地勢和水文地質(zhì)影響,形成了較大面積的沼澤濕地。

撓力河流域整體呈現(xiàn)西南向東北傾的特點(diǎn),地形相對簡單;受構(gòu)造及河流的影響,根據(jù)地貌基本成因類型,本區(qū)可劃分為堆積地形(I1、I2)、剝蝕堆積(II1)、火山堆積(II2)及侵蝕剝蝕地形(III1、III2)四種。平原上除零星分布孤山和殘丘,高度多在500 m以下,主要由古生代、中生代頁巖,中酸性火山巖和和花崗巖所構(gòu)成;多數(shù)區(qū)域在復(fù)合的褶皺基底上堆積有千米以上的中、新生代沉積蓋層,第四系厚度一般為120～200 m,最厚可達(dá)220 m;大部分孔隙水含水層為砂和砂礫石層,單井出水量一般為150～300 m3/h,適合農(nóng)業(yè)開采。此次模擬該流域第四系孔隙含水層系統(tǒng),由于第四系較厚,存在4～5個(gè)含水層,但各含水層間無明顯隔水層,且弱透水層不連續(xù),形成統(tǒng)一大厚度含水層[6]。

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

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

2.1.1 水文地質(zhì)條件概化

研究區(qū)域以撓力河流域灌區(qū)為核心,向西、向北擴(kuò)到撓力河與松花江的分水嶺,向南、向東分別延伸到平原區(qū)與山丘區(qū)的交界線;研究區(qū)較流域面積小,主要為第四系覆蓋區(qū),總面積約為1.40萬 km2。

撓力河流域表層廣泛覆蓋有3～20 m的亞黏土,因此地下水具弱承壓性(水頭6～7 m)。含水層為第四系松散巖類砂及砂礫石層。撓力河流域水平方向上由山前臺地?cái)嗬m(xù)分布的第四系微孔隙裂隙水到撓力河流域平原地區(qū)的第四系孔隙弱承壓水逐步過渡,松散巖類含水層逐漸變厚,粒徑逐漸變細(xì)。

依據(jù)本區(qū)域地質(zhì)情況和實(shí)測水文資料,建立適合本研究區(qū)的水文地質(zhì)概念模型;該模型在垂向上將含水層分為2層,即上層為潛水含水層,下層為承壓含水層。上層含水層(潛水層)巖性主要是亞黏土,在撓力河流域分布廣泛;下層含水層(承壓層)主要為砂及砂礫石。含水介質(zhì)概化成非均質(zhì)、各向同性。

研究區(qū)AB界線為河流流出區(qū),側(cè)向流出研究區(qū),以側(cè)向徑流為主;BC和CD界線為撓力河與松花江的流域分水嶺;DE和EA界線為山區(qū)與平原區(qū)的分界線,邊緣山地側(cè)向徑流補(bǔ)給平原區(qū)(圖1);因此,整個(gè)區(qū)域的側(cè)向邊界均概化成第三類邊界。撓力河、七星河在區(qū)內(nèi)流量較大,與第四系松散巖類含水系統(tǒng)水量交換廣且大,可概化為河流邊界。研究區(qū)上部以潛水面為界,經(jīng)由該邊界淺層地下水與外部環(huán)境產(chǎn)生水量交換;底部以第四系松散層底板為界,概化為隔水邊界[7]。撓力河流域松散巖類含水層分布廣且厚度大,符合達(dá)西定律,概化為層流。潛水和承壓含水層之間有水量交換,概化為三維非穩(wěn)定流。研究區(qū)的地下水補(bǔ)給項(xiàng)可概化為4項(xiàng):大氣降水入滲、灌溉回滲(井灌和渠灌回灌、渠系回滲)、基巖區(qū)地下水徑向徑流補(bǔ)給和河流滲漏補(bǔ)給[8];地下水排泄項(xiàng)也可以概化為4項(xiàng):蒸發(fā)量、人工開采量(包括工業(yè)、農(nóng)業(yè)和生活用水開采量)、河流排泄量、側(cè)向徑流排泄量[8]。

圖1 研究區(qū)水文地質(zhì)概念模型

2.1.2 水文地質(zhì)參數(shù)分區(qū)

主要的水文地質(zhì)參數(shù)為滲透系數(shù)與給水度/貯水率。根據(jù)撓力河流域地貌類型和水文地質(zhì)特點(diǎn),將模擬區(qū)潛水含水層滲透系數(shù)和給水度劃分為2個(gè)參數(shù)分區(qū),深層含水層滲透系數(shù)和貯水率分為8個(gè)參數(shù)分區(qū)(圖2)。本次研究主要依據(jù)包氣帶巖性對降水入滲系數(shù)以及極限蒸發(fā)深度進(jìn)行分區(qū)(圖3)。

圖2 潛水含水層(a)和承壓含水層(b)水文地質(zhì)參數(shù)分區(qū)圖

圖3 包氣帶巖性分區(qū)及入滲、蒸發(fā)參數(shù)分區(qū)圖

2.2 地下水流數(shù)學(xué)模型

(1)

2.3 模型求解

采用Visual Modflow軟件系統(tǒng),對撓力河流域地下水?dāng)?shù)學(xué)模型進(jìn)行求解,該軟件具有強(qiáng)大計(jì)算功能以及良好的圖形可視界面功能,得到國際同行的廣泛認(rèn)可[9]。

2.3.1 區(qū)域剖分和時(shí)空離散

根據(jù)地下水位動態(tài)監(jiān)測數(shù)據(jù)的時(shí)空分布規(guī)律,識別期設(shè)為2012年1月1日-2012年12月31日,驗(yàn)證期為2013年1月1日-2013年12月31日,以5 d為1個(gè)時(shí)段,共69個(gè)時(shí)段;該時(shí)段經(jīng)歷了枯、平、豐水期,能更好的反映含水層系統(tǒng)的特征。

計(jì)算區(qū)采用矩形網(wǎng)格剖分法,設(shè)定網(wǎng)格剖分尺寸為1 385 m×1 385 m,單個(gè)計(jì)算單元面積為1.9 km2。

2.3.2 含水層參數(shù)

根據(jù)撓力河流域?qū)嶋H水文地質(zhì)情況及既往研究成果,按照圖3所示含水層分區(qū),設(shè)定模型滲透系數(shù)與給水度/貯水率,根據(jù)各單元含水層巖性特征賦各參數(shù)分區(qū)的初始值(表1)。按照圖4所示降水入滲系數(shù)、灌溉回歸系數(shù)以及極限蒸發(fā)深度分區(qū)進(jìn)行賦初值(表2)。

表1 含水層水文地質(zhì)參數(shù)初值

表2 大氣降水入滲系數(shù)以及極限蒸發(fā)深度初始值

圖4 模擬值與觀測值擬合散點(diǎn)圖(a)和模擬誤差分布圖(b)

2.3.3 模型識別和驗(yàn)證

模型識別是將所建立模型的邊界條件、源匯項(xiàng)、水力特征、水文地質(zhì)參數(shù)等與實(shí)際情況相吻合,通過不斷調(diào)整模型參數(shù),各時(shí)段的地下水位與實(shí)測的地下水位擬合度較好(圖4a),兩者之間的相關(guān)系數(shù)可達(dá)0.99;誤差分布直方圖基本呈現(xiàn)正態(tài)分布,擬合誤差最小點(diǎn)出現(xiàn)的頻率位于正態(tài)分布曲線的中心部位(圖4b)。

經(jīng)識別后的含水層滲透系數(shù)、給水度的參數(shù)分區(qū)值如表3;大氣降水入滲系數(shù)分區(qū)、大氣降水入滲系數(shù)分區(qū)值如圖5、表4所示;極限蒸發(fā)深度分區(qū)及賦值與初始一致(表2)。

表3 含水層滲透系數(shù)和給水度初值

表4 大氣降水入滲系數(shù)識別值

圖5 大氣降水入滲系數(shù)識別分區(qū)圖

將識別所得的水文地質(zhì)參數(shù)輸入到模型中,驗(yàn)證期為2013年1月1日-2013年12月31日,以5天為1個(gè)時(shí)段,檢驗(yàn)下水位計(jì)算值與觀測值的擬合情況,結(jié)果顯示地下水位擬合效果良好,計(jì)算值與實(shí)測值的相關(guān)系數(shù)達(dá)到0.993(圖6a);誤差分布直方圖呈現(xiàn)正態(tài)分布,擬合誤差最小點(diǎn)出現(xiàn)的頻率位于正態(tài)分布曲線的中心部位(圖6b)。

圖6 模擬值與觀測值擬合散點(diǎn)圖(a)和模擬誤差分布圖(b)

經(jīng)過模型的識別與驗(yàn)證,所建模型的邊界條件、源匯項(xiàng)、水力特征、水文地質(zhì)參數(shù)等與實(shí)際情況擬合效果良好,能夠正確反映區(qū)域水文地質(zhì)條件以及地下水運(yùn)動規(guī)律,可以采用此模型對研究區(qū)地下水動態(tài)規(guī)律開展模擬研究。

3 地下水資源開采潛力評價(jià)

Modflow軟件中的Zone Budget模塊可對整個(gè)區(qū)域進(jìn)行水量均衡分析,也可對局部區(qū)域進(jìn)行水量均衡分析;利用Zone Budget模塊可統(tǒng)計(jì)出模型中的區(qū)域入滲補(bǔ)給量、邊界補(bǔ)給量、河流補(bǔ)給量以及排泄量等數(shù)據(jù)。利用Zone Budget模塊,將研究區(qū)內(nèi)的各個(gè)灌區(qū)設(shè)為具體的水量均衡區(qū)進(jìn)行統(tǒng)計(jì)分析(圖7)。

圖7 水量均衡統(tǒng)計(jì)分區(qū)示意圖

模型中地下水總補(bǔ)給量包括入滲補(bǔ)給量、河道滲漏補(bǔ)給量、邊界側(cè)向徑流補(bǔ)給量、區(qū)間側(cè)向徑流補(bǔ)給量等。即:

Q總補(bǔ)= Q滲補(bǔ)+ Q河補(bǔ)+ Q邊補(bǔ)+ Q側(cè)補(bǔ)

(2)

式中:Q總補(bǔ)為地下水總補(bǔ)給量,104m3/a;Q滲補(bǔ)為包括:降水入滲補(bǔ)給量、井灌回歸補(bǔ)給量、渠系滲漏補(bǔ)給量、渠灌田間入滲補(bǔ)給量,104m3/a;Q河補(bǔ)為河道滲漏補(bǔ)給量,104m3/a;Q邊補(bǔ)為由通用水頭邊界(GHB)計(jì)算出的邊界側(cè)向徑流補(bǔ)給量,104m3/a;Q側(cè)補(bǔ)為相鄰灌區(qū)之間側(cè)向徑流補(bǔ)給量,104m3/a。

灌區(qū)現(xiàn)狀年地下水可開采量的計(jì)算采用可開采系數(shù)法,即:

Q可采= ρ· Q總補(bǔ)

(3)

式中:Q可采為地下水可開采量,104m3/a;Q總補(bǔ)為地下水總補(bǔ)給量,104m3/a。ρ為可開采系數(shù)(ρ≤1),由于本流域地下水富水性較好,因此對于單井單位降深出水量>20 m3/h·m,且地下水埋深大、水位連年下降的超采區(qū),ρ值的取值范圍為0.8～1.0;對單井單位降深出水量在5～10 m3/h·m,地下水埋深大、實(shí)際開采程度較高地區(qū)或地下水埋深較小、實(shí)際開采程度較低地區(qū),ρ值的取值范圍為0.6～0.8[10-12]。

經(jīng)計(jì)算可得撓力河流域灌區(qū)現(xiàn)狀年地下水總補(bǔ)給量、地下水可開采量、資源量及可開采量,詳細(xì)計(jì)算見表5。

表5 灌區(qū)現(xiàn)狀年地下水總補(bǔ)給量、可開采量及資源量計(jì)算表

基于表5的結(jié)果,參考楊國強(qiáng)(2012)在地下水超采指標(biāo)體系建立中,利用開采系數(shù)法對研究區(qū)地下水開發(fā)利用情況劃分等級,即:小于95%為潛力區(qū),95%～105%為補(bǔ)給平衡區(qū),105%～120%為一般超采區(qū),大于120%為嚴(yán)重超采區(qū)[13]。對撓力河地下水開發(fā)利用情況劃分為5個(gè)等級,即:Ⅰ級開采率低于50%(潛力較大);Ⅱ級50%～80%(潛力中等);Ⅲ級80%～100%(潛力較小);Ⅳ級100%～120%(一般度超采);Ⅴ級大于120%(重度超采)。

由表6可知,撓力河流域灌區(qū)現(xiàn)狀年總的可開采量為1.44×109m3/a,總的開采量為1.63×109m3/a,地下水開采率為113.1%,灌區(qū)總體呈超采狀態(tài)。但各灌區(qū)地下水資源的承載力不同:(1)蛤蟆通灌區(qū)、七里沁灌區(qū)、尖山子灌區(qū)、三環(huán)泡灌區(qū)、龍頭橋灌區(qū)開采率低于50%,屬開采潛力較大的灌區(qū),尚具有較大的農(nóng)業(yè)灌溉潛力;(2)集賢小黃河灌區(qū)、幸福灌區(qū)開采率在50%～80%之間,屬開采潛力中等,還可支持一定規(guī)模的農(nóng)業(yè)灌溉面積發(fā)展;(3)富錦灌區(qū)地下水開采率為83.5%,開采率在80%～100%之間,具有較小的開采潛力,如果繼續(xù)大規(guī)模開采,將有超采的可能。(4)二九一南部灌區(qū)、友誼西部灌區(qū)開采率在100%～120%之間,屬一般超采區(qū)。灌區(qū)地下水開采量已經(jīng)超過該灌區(qū)地下水可開采量,如果灌區(qū)地下水利用量不加調(diào)節(jié)則地下水位將持續(xù)下降。(5)錦西灌區(qū)、五九七灌區(qū)、東泄總灌區(qū)、錦南灌區(qū)、大興灌區(qū)、八五三灌區(qū)地下水開采率大于120%,屬重度超采,部分灌區(qū)已出現(xiàn)地下水降落漏斗,若不加調(diào)節(jié)將可能引發(fā)嚴(yán)重的環(huán)境地質(zhì)問題。

表6 灌區(qū)地下水開發(fā)利用程度分析

4 承載力綜合評價(jià)結(jié)果合理性分析

地下水資源承載力強(qiáng)弱程度是區(qū)域內(nèi)地下水資源系統(tǒng)與社會、環(huán)境綜合作用的結(jié)果,直觀地表現(xiàn)在地下水開采潛力上,區(qū)域地下水承載力強(qiáng)則表現(xiàn)為地下水開采潛力大,反之則表現(xiàn)為地下水開采潛力小;根據(jù)這一特性,本文采用已有的撓力河流域灌區(qū)地下水承載力綜合評價(jià)結(jié)果[3],其結(jié)果見表8;并以模擬所得出的結(jié)果作為參考,評價(jià)其合理性。

因兩個(gè)體系的分級不同,因此為了便于兩者的比較,根據(jù)兩者等級劃分標(biāo)準(zhǔn),給出了兩者的對比標(biāo)準(zhǔn)(表7),當(dāng)兩者一致度達(dá)到80%以上時(shí)即可接受,當(dāng)?shù)陀?0%時(shí)需要重新修訂指標(biāo)體系[14]。參考已有的地下水承載力綜合評價(jià)結(jié)果,將地下水承載力綜合評價(jià)結(jié)果和地下水開采潛力分析結(jié)果匯總于表8。

表7 一致性判別表

表8 地下水承載力綜合評價(jià)結(jié)果合理性分析

由表8可知,從排序上看,除二九一南部灌區(qū)、富錦灌區(qū)、五九七灌區(qū)和三環(huán)泡灌區(qū)外,兩種方法所得的評價(jià)結(jié)果都比較相近;但從評價(jià)等級上看兩者相差較大,二九一南部灌區(qū)、錦西灌區(qū)、錦南灌區(qū)、友誼西部灌區(qū)、五九七灌區(qū)地下水承載力為中或高等水平,但運(yùn)用數(shù)值模擬分析的實(shí)際情況為嚴(yán)重超采區(qū),開采潛力很弱。此外,地下水承載力綜合評價(jià)結(jié)果在等級劃分上不夠細(xì)化,如:集賢小黃河灌區(qū)、幸福灌區(qū)地下水承載力為高等,但該灌區(qū)實(shí)際開采潛力為中等水平,這就導(dǎo)致評價(jià)結(jié)果出現(xiàn)偏差,影響評價(jià)結(jié)果的合理性。從一致性判斷上分析,兩者的一致性為68.7%,一致性低,需要對指標(biāo)體系重新修訂。

分析其中的原因主要包括4個(gè)方面:(1)撓力河流域乃至整個(gè)三江平原區(qū)域內(nèi)耕地面積多,人口少,以人口量作為區(qū)域衡量指標(biāo)不合理,所選人均供水量指標(biāo)欠妥。(2)區(qū)域地下水資源可利用量計(jì)算較復(fù)雜,此外,如果已知區(qū)域地下水開發(fā)利用程度,也就已知區(qū)域地下水整體開采情況,再以地下水水量為基礎(chǔ)選取指標(biāo)建立指標(biāo)體系,進(jìn)而評價(jià)地下水資源承載力是不合理的。(3)撓力河流域灌區(qū)內(nèi)第二、三產(chǎn)業(yè)幾乎不存在,而農(nóng)業(yè)較發(fā)達(dá),在后續(xù)選取指標(biāo)時(shí),應(yīng)以“地下水—農(nóng)業(yè)”二者之間關(guān)系為主。(4)關(guān)于地下水資源量,其計(jì)算也比較復(fù)雜,對于撓力河流域各灌區(qū)灌溉期集中,且灌區(qū)地下水主要補(bǔ)給來源為降水入滲補(bǔ)給,可以降水入滲量作為地下水資源量。

綜合以上分析結(jié)果來看本次地下水資源承載力評價(jià)所采用的指標(biāo)體系尚存在一定的不合理性。因此,指標(biāo)選取需要仔細(xì)甄別,指標(biāo)分級標(biāo)準(zhǔn)也需要細(xì)化,以便構(gòu)建一套合理的地下水資源承載力評價(jià)指標(biāo)體系。

5 結(jié)語

(1)利用地下水?dāng)?shù)值法,建立撓力河流域灌區(qū)地下水?dāng)?shù)值模型,計(jì)算各個(gè)灌區(qū)的總補(bǔ)給量、地下水可開采量、地下水開采率等,進(jìn)而分析各灌區(qū)地下水資源開采潛力,可對基于指標(biāo)體系的地下水承載力評價(jià)結(jié)果的合理性和可靠性起到評估作用。

(2)以數(shù)值模擬所得出的結(jié)果作為參考,評價(jià)已有的撓力河流域灌區(qū)地下水承載力綜合評價(jià)結(jié)果的合理性;結(jié)果表明:所選指標(biāo)體系具有一定的參考價(jià)值,但存在選取不夠合理和分級不夠細(xì)化問題。

(3)指標(biāo)體系的構(gòu)建對地下水承載力評價(jià)結(jié)果的合理性和可靠性有較為顯著的影響,在實(shí)際工作中應(yīng)慎重選擇,并在可能的條件下采用多種方法進(jìn)行綜合分析。

(4)本文采用的評價(jià)方法較為成熟,但評價(jià)思路尚未有學(xué)者深入研究,結(jié)論仍需在今后實(shí)際工作研究中進(jìn)行驗(yàn)證。