基于gis-pmodflow的西部地下水潛水位數(shù)值模擬與水環(huán)境預(yù)警研究

基于gis-pmodflow的西部地下水潛水位數(shù)值模擬與水環(huán)境預(yù)警研究

應(yīng)急的概念最早應(yīng)用于軍事領(lǐng)域的雷達(dá)技術(shù)和彈頭防御系統(tǒng)。近年來，應(yīng)急理論在洪水報告、經(jīng)濟系統(tǒng)、氣象、機械工程等方面發(fā)揮了重要作用。然而，在資源和環(huán)境科學(xué)領(lǐng)域，對報警理論的研究和應(yīng)用相對較少。目前，關(guān)于水環(huán)境警報的研究文獻(xiàn)很少。松嫩平原西部地域遼闊,屬半干旱半濕潤的大陸性季風(fēng)氣候區(qū),多年平均降水量為400～500mm,其中6～9月份降水量占全年總降水量的74%～84%。多年平均蒸發(fā)量自東部向西部逐漸增大,為1400～2100mm。有較豐富的水、土、生物資源和石油資源,為我國重要的農(nóng)、牧業(yè)基地和能源基地。該區(qū)河流稀少,分布不均,閉流面積很大,地表徑流量少,地下水較為豐富。近年來,由于上游地區(qū)興建了水庫,河水被截流,西部和中部的洮兒河、霍林河及蛟流河下游河水水量減少,甚至河床干涸,地下水的補給量也呈減少趨勢。在全球氣候變暖及人類開發(fā)強度加大的雙重影響下,生態(tài)環(huán)境急劇惡化,水資源匱乏,土地鹽堿化、沙化和退化現(xiàn)象十分嚴(yán)重。干旱、土地鹽堿化和局部地下水位下降已成為該區(qū)主要的水環(huán)境問題。多年來,中國科學(xué)院、大專院校以及生產(chǎn)部門,在該區(qū)進(jìn)行了大量科學(xué)研究和生產(chǎn)實踐,為水資源的合理開發(fā)和生態(tài)環(huán)境的綜合治理做出了貢獻(xiàn)。但是從地下水水位變幅開展水資源環(huán)境預(yù)警尚屬首次。本文采用GIS-Modflow聯(lián)合系統(tǒng),通過對水文地質(zhì)條件的概化,建立起潛水水流數(shù)值模型,借助于PModflow軟件包,應(yīng)用有限差分法對潛水流場進(jìn)行了數(shù)值模擬并做出預(yù)測,以揭示潛水水流的運移機制,并對潛水水位進(jìn)行預(yù)報。在GIS平臺中通過現(xiàn)狀水位與警戒水位進(jìn)行比較運算,以此對現(xiàn)狀和未來的水環(huán)境進(jìn)行預(yù)警。1研究區(qū)地下水位與土壤鹽堿化關(guān)系的探討本區(qū)東、南、西三面地形較高,北部較低,為一個簸箕形的含水盆地。地下水匯水范圍廣闊,補給來源比較充沛,主要有大氣降水垂直補給、霍林河下游散流區(qū)垂直補給,山區(qū)地下水側(cè)向補給和河流側(cè)向補給。在西北山區(qū)和長嶺松遼分水嶺及高平原為地下水側(cè)向補給區(qū)。洮兒河、蛟流河中上游地表徑流在沿河兩側(cè)補給地下水,松花江、嫩江在汛期也有部分地表水補給地下水。本文以潛水水位變化幅度作為主要研究對象。采用以不至于發(fā)生土壤鹽堿化的潛水臨界水位作為上第四紀(jì)孔隙潛水(Q3-Q4),含水層巖性為粉砂、細(xì)砂、砂礫石。含水層厚度由南向北逐漸增厚,從1～5m到18～23m不等,最大為40m。潛水的埋深一般為1～3m、3～5m,崗地為5～10m,山前傾斜平原為10～30m。潛水是該區(qū)進(jìn)行農(nóng)業(yè)灌溉的一種重要水源。警戒線值,以潛水開采極限深度作為下警戒線值。經(jīng)過長期觀測和綜合研究,設(shè)定研究區(qū)潛水含水層的開采極限深度為由地表到潛水含水層厚度的1/2處的距離,超過該深度開發(fā)水資源將會產(chǎn)生枯竭。松嫩平原西部地域遼闊,地貌類型、土層結(jié)構(gòu)和水文地質(zhì)條件復(fù)雜多變,獲取全區(qū)的潛水臨界深度十分困難。為此,提出防止土壤積鹽的合理潛水深度的概念,即將土壤鹽漬化速度控制在較低值時地下潛水水位應(yīng)滿足的變幅。以1996年的潛水位為基準(zhǔn)值,利用地球動力學(xué)模型進(jìn)行反演,求出合理的潛水水位。根據(jù)地球動力學(xué)原理,若要確定自然因素和人為因素對地下水系統(tǒng)的影響,應(yīng)從以下方面考慮:地下水動態(tài)變化的產(chǎn)出要素的增長,是由投入要素的增加及人文作用對地下水系統(tǒng)的影響產(chǎn)生的(這里的“增長”可能是負(fù)的)。人文作用在地下水動態(tài)變化產(chǎn)出要素中的貢獻(xiàn)率可由下式計算得到:Κ1=|ΔΜ/ΜΔY/Y|×100％(1)K1=∣∣ΔM/MΔY/Y∣∣×100％(1)式中Y為潛水水位值,ΔY/Y為地下水潛水位變化;M為人類活動對地下水的影響程度,本文采用土壤鹽堿化程度,ΔM/M為人類活動的貢獻(xiàn)率;K1為影響強度系數(shù)。根據(jù)1983年、1996年潛水位資料和1983年、1995年土地鹽堿化遙感解譯數(shù)據(jù),利用GIS圖形疊加及分析功能,將潛水水位和土地鹽堿化圖疊加形成潛水水位對土地鹽堿化的影響強度柵格圖(略)和在人為作用下潛水水位對土壤鹽堿化的影響強度柵格圖(略)。通過對土壤類型圖、土地鹽堿化圖和地下水潛水位圖的疊加分析與計算發(fā)現(xiàn),潛水位對土壤鹽漬化的影響系數(shù)大于0.5的地區(qū)占計算區(qū)的10.28%,影響系數(shù)介于0.3～0.5之間的地區(qū)占12.77%,這些地區(qū)是鹽堿化發(fā)展嚴(yán)重的地區(qū)。假設(shè)在1996年至2015年將土壤鹽堿化的發(fā)展速度控制在1%,即該區(qū)每個單元土壤鹽堿化程度增加1%。其地下水潛水位和土壤鹽漬化關(guān)系如下式。Η2-Η1Η1×Μ1Μ2-Μ1＜50％或30％(2)H2?H1H1×M1M2?M1＜50％或30％(2)即表示在人類作用下潛水水位對土壤鹽堿化的影響強度系數(shù)大于0.5～0.3的單元達(dá)到50%或30%。式中H1為單元1996年的潛水水位;H2為單元2015年的潛水水位;M1為單元1996年的土壤鹽漬化強度;M2為單元2015年土壤鹽漬化的強度。若已知H1、M2、M1,通過式(2)可求出H2的值。通過對該區(qū)遙感解譯數(shù)據(jù)的分析,得到了本區(qū)1989-2001年12年間土壤鹽堿化的發(fā)展速度為2.61%/a,為保護(hù)生態(tài)環(huán)境,遏制土壤鹽堿化的發(fā)展,我們認(rèn)為該區(qū)土壤鹽漬化加重速度應(yīng)控制在1%,在此條件下2015年人類活動所控制的潛水水位值為H2。應(yīng)用GIS將各單元的H2值繪制成松嫩平原西部平原水環(huán)境預(yù)警潛水上警戒線圖,見圖1。2水位預(yù)測模型的構(gòu)建2.1系統(tǒng)模型的簡化(1)地質(zhì)概念模型研究區(qū)潛水含水層可劃分:①西部山前傾斜平原扇形地的孔隙潛水含水層,由中、上更新統(tǒng)沖洪積砂卵礫石、卵石層組成,顆粒粗大,滲透性極強(K=140～300m/d);②河谷平原、東部高平原及霍林河流域中游孔隙潛水含水層,由砂及砂礫石層組成,滲透性較強(K=10～100m/d);③松拉河間地塊及低平原西部的孔隙潛水含水層,由細(xì)砂、中細(xì)砂及粉細(xì)砂組成,滲透性較弱(K=10～20m/d);④低平原中部乾安、大安一帶的孔隙潛水含水層,多由黃土狀亞砂土、粉細(xì)砂及亞砂土組成,顆粒很細(xì),滲透性很差(K=3～10m/d)。研究區(qū)水文地質(zhì)概念模型見圖2所示。根據(jù)圖2可以了解研究區(qū)邊界的分布、性質(zhì)、水文地質(zhì)參數(shù)分區(qū)及觀測井分布等情況。(2)流量邊界及隔水邊界計算區(qū)西北角的洮兒河沖洪積扇形地,與其周邊以火山巖為主的基巖山區(qū)相鄰,滲透性較弱,可將西北部邊界概化為隔水邊界。對于常年有水的洮兒河和蛟流河河谷,則可作為流量邊界,利用上下游水文站的徑流資料和地下水等水位線圖,進(jìn)行綜合分析,利用水量均衡關(guān)系計算地表水向地下水的轉(zhuǎn)化量。南部邊界以地下分水嶺為界,概化為二類零流量邊界。計算區(qū)東北部以嫩江、松花江、拉林河和第二松花江常年有水的江河為界,均為潛水的排泄邊界,概化為一類水頭邊界,正東為高平原的地下分水嶺,概化為隔水邊界。除以上各邊界外,其它邊界上的水位動態(tài)觀測井較多,故可概化為已知水位的一類邊界。(3)隔水邊界為水量交換邊界,一個是下根據(jù)潛水位的動態(tài)變化資料,降水入滲、河道滲漏、人工開采、灌溉滲漏、潛水蒸發(fā)均會影響潛水的變化,因而,上部邊界為水量交換邊界;下部邊界則分兩種情況,當(dāng)在單一潛水含水層分布區(qū),如沖洪積扇形地孔隙潛水區(qū)或河谷平原區(qū),其下部均為相對隔水的底板,可概化為隔水邊界;上部為孔隙潛水,下部為孔隙承壓含水層的雙層含水層,則垂向邊界屬于越流邊界,研究區(qū)越流層為大青溝組的弱透水層。2.2含水層滲透系數(shù)計算基因根據(jù)上述的水文地質(zhì)概念模型,建立潛水水流數(shù)學(xué)模型如下:{??x[Κ(Η-B)?Η?x]+??y[Κ(Η-B)?Η?y]-ΚΜ(Η0-Η)+Qr-Qd-n∑i=1Qi(xi?yi?t)δ(x-xi?y-yi)=μ?Η?th(x?y?t)|t=0=h0(x?y?t)(x?y)∈Dh(x?y?t)|Γ1=h1(x?y?t)(x?y)∈D?t＞0Κ(h-B)?h?→n|Γ2=q(x?y?t)(x?y)∈D?t＞0(3)式中K為潛水含水層滲透系數(shù),m/d;μ為潛水含水層給水度;h、H0為地下水位、潛水水位,m;B為含水層底板標(biāo)高,m;Qr為入滲補給強度,m/d;Qd為排泄強度(包括開采),m/d;Qi為井開采量,m3/d;h0為初始水位,m;h1為類邊界點的水位,m;q為二類邊界單寬流量,m2/d;x,y為坐標(biāo),m;D為計算區(qū)范圍;Γ1,Γ2為一類及二類邊界;→n為邊界上的內(nèi)法線;n為開采井總數(shù)。2.3變采用單元計算單元本計算區(qū)域采用傳統(tǒng)的四邊形剖分,將含水層剖分為一層、38行、70列,共2660個單元。其中定水頭單元149個,無效計算單元753個,變水頭計算單元(即有效計算單元)1408個。單元尺寸與警戒線柵格單元尺寸保持一致。計算區(qū)面積為35201.17km2,占整個松嫩平原西部面積的74.77%,單元平均面積為25km2。行間距(ΔCi)和列間距(Δri)均為5000m,厚度(Δvk)視具體情況而定。2.4地下水位計算模型的驗證將計算單元的水頭與實測水頭進(jìn)行對比,進(jìn)而反求相關(guān)的水文地質(zhì)參數(shù)。模型校正時段選取1996年10月15日到1997年3月30日,時間步長為15d,共11個時段。在該時段內(nèi)的源匯項少,地下水位處于平水期,計算較簡單。潛水含水層計算模型參數(shù)分區(qū)基本上與區(qū)域水文地質(zhì)單元范圍一致。各參數(shù)區(qū)的初值根據(jù)已有勘查和研究階段所進(jìn)行的抽水試驗成果而定。模型識別時段為冬季,無大氣降水入滲補給和農(nóng)業(yè)灌溉,潛水的蒸發(fā)可以忽略不計。越流補給(排泄)與河流(湖泊)滲漏補給(排泄)等源匯項可直接加入相應(yīng)的計算單元,城鎮(zhèn)生活用水則按強度分配到單元上。最后將源匯項和初始地下水位輸入數(shù)學(xué)模型,通過水文地質(zhì)參數(shù)的匹配,直到調(diào)參時段末地下水水位計算值與實測值擬合誤差達(dá)到要求為止。檢驗結(jié)果表明,擬合誤差<0.5m者達(dá)到70%以上,說明所建立的水文地質(zhì)概念模型和數(shù)學(xué)模型是正確可靠的。為了進(jìn)一步驗證其可靠性,對所建立的區(qū)域數(shù)學(xué)模型進(jìn)行了檢驗。選擇地下水年內(nèi)低水位期(1997年3月11日至5月11日,共計61d)和地下水位上升期(1997年5月11日至8月11日,共計92d)作為檢驗時段,以低水位期的水位作為初始流場,將各時段的源匯項輸入已識別后的數(shù)學(xué)模型,并計算出兩個時段末刻的水位值。兩個檢驗時段末計算水位與實測水位的擬合誤差檢驗結(jié)果表明,各時段水位觀測點擬合誤差(0.5m的井點數(shù)達(dá)到總水位點總數(shù)的73%,所建的數(shù)學(xué)模型可用于地下水預(yù)報。2.5動滲流模型的建立及模型預(yù)測的依據(jù)根據(jù)研究區(qū)的遠(yuǎn)景發(fā)展規(guī)劃目標(biāo),將水位預(yù)報期限定為15年,即1999-2015年。鑒于1999年5月的地下水動態(tài)觀測資料較全,水位觀測井較多,故將此階段的潛水位作為預(yù)報基礎(chǔ)值。邊界及源匯項的處理降水量是地下水的主要補給來源,研究中采用Monte-Carlo方法進(jìn)行預(yù)報。對研究區(qū)一類邊界選取10個有長期水位觀測資料的邊界井,求出水位時間關(guān)系方程,相關(guān)檢驗結(jié)果均較顯著。其它邊界井則按插值法得出。潛水模型源匯項的預(yù)測主要依據(jù)該區(qū)各行業(yè)用水規(guī)劃及水利規(guī)劃,同時根據(jù)源匯項與時間或其它易于確定變量的相關(guān)關(guān)系進(jìn)行預(yù)測,如河流與地下水的關(guān)系,采用河水滲漏與降水及河流徑流量之間的相關(guān)模型進(jìn)行。水位預(yù)報結(jié)果及分析從預(yù)報結(jié)果(圖3)中可以看出,到2015年時,計算區(qū)內(nèi)潛水流場的整體形態(tài)及地下水流動方向均未發(fā)生明顯變化,水位降深值較小,最大降深值<5.5m。在預(yù)報時段內(nèi),區(qū)內(nèi)潛水位多數(shù)地段僅隨年降水量有微小的變動。在開采量較大的白城沖洪積扇區(qū),水位出現(xiàn)了持續(xù)的下降現(xiàn)象,但降幅較小,平均年降幅為0.25m。在前郭縣和鎮(zhèn)賚縣以地表水為水源的水田區(qū),由于灌溉水的回滲補給,潛水位有所上升,如143號和210號井潛水位在15年期間有12年水位上升。這一現(xiàn)象表明,水田的開發(fā)會增加潛水的補給量,但隨著潛水位的上升,將出現(xiàn)鹽堿化發(fā)展的危險。3超標(biāo)下東北部地區(qū)—潛水環(huán)境預(yù)警在水位預(yù)報的基礎(chǔ)上,通過計算現(xiàn)狀和預(yù)測期水位,并與警戒水位進(jìn)行對比,以達(dá)到預(yù)警的目的。利用GIS空間分析功能進(jìn)行柵格圖層疊加和屬性判斷,其中運用了多層?xùn)鸥駭?shù)據(jù)復(fù)合中的算術(shù)運算方法和矢量多邊形疊置分析方法。應(yīng)用MapInfo軟件,將警戒線圖層、潛水流線現(xiàn)狀圖層和預(yù)測圖層經(jīng)柵格化處理后,進(jìn)行預(yù)警判斷。分別以超出上、下警戒線一定的數(shù)據(jù)區(qū)間為不同警度值,對警情進(jìn)行分級。結(jié)果見圖4和表1。超出上警戒線警情主要分布在較大型的灌區(qū)及其周邊,如前郭灌區(qū)、鎮(zhèn)賚灌區(qū)等地,警情的存在和警情的等級說明灌區(qū)的排水工程未能充分發(fā)揮作用,加之灌區(qū)多年引地表水或地下水灌溉人為抬高了潛水水位所致。到2015年,若此狀況不能改變,將會引起警情的擴大。還有一部分警情發(fā)生在大安的西南與乾安、通榆一帶,該地區(qū)地勢較為低平,存在許多碟形封閉洼地、泡沼和低平地,農(nóng)田灌溉時排水不暢,引起潛水水位抬升,導(dǎo)致土壤次生鹽堿化。超出下警戒線的警情主要分布在松-拉河間地塊以及西部山前傾斜平原水資源相對較豐富、開采程度較大的地區(qū)。松-拉河間地塊大部分地區(qū)潛水水位高于承壓水位,在人工開采孔隙承壓水的影響下,潛水發(fā)生越流補給承壓含水層。在大量開采承壓含水層的情況,越流作用更加明顯,加之大量開采該層水資源,導(dǎo)致水位降至警戒線之下。西部山前傾斜平原的潛水也以越流的形式補給承壓含水層。近年來該區(qū)氣候漸趨干旱,降水量減少,地下水的補給來源明顯減小,同時,地下水的開采量逐年增加,使得部分地區(qū)水位持續(xù)下降以至超過了下警戒線。4松嫩平原西部水環(huán)境預(yù)警方法水資源環(huán)境預(yù)警是一個新的研究課題,它涉及到許多新理論和新方法。在地下水環(huán)境預(yù)警中,明確預(yù)警的目標(biāo)、確定警戒線十分重要。采用GIS-PModflow系統(tǒng)進(jìn)行地下水環(huán)境的預(yù)警是一個新探索。預(yù)警結(jié)果表明,該區(qū)