基于Visual Modflow的太子河流域地下水流數(shù)值模擬分析

徐世俊

(遼寧省遼陽(yáng)水文局，遼寧 遼陽(yáng) 111000)

水是人們賴以生存的必要條件，日常生活用水以及農(nóng)業(yè)灌溉用水更為重要。地下水在生活用水、工業(yè)用水以及農(nóng)田灌溉用水等方面發(fā)揮著至關(guān)重要的作用[1]。為了預(yù)防水資源污染、過(guò)度開(kāi)采水資源引起的嚴(yán)重缺水并合理化配置水資源，水流動(dòng)態(tài)模擬分析一直被國(guó)內(nèi)外學(xué)者認(rèn)為是一種非常重要的數(shù)學(xué)模擬方法。該數(shù)值模擬軟件是由美國(guó)著名學(xué)者采用三維有限差分法開(kāi)發(fā)研制而成的，其適用范圍為孔隙介質(zhì)地下水。Visual Modflow模型軟件是一種地下水模擬可視軟件，由多種地下水模型軟件的功能屬性如Modflow、Modpath、MT3D、RT3D、WIN PEST等共同形成[2]。

模擬分析范圍的確定受多個(gè)水文地質(zhì)因素影響，如地形地貌、地質(zhì)結(jié)構(gòu)、巖層分布、天氣環(huán)境、農(nóng)業(yè)實(shí)用、地下水流情況等，其中地質(zhì)條件屬性受均質(zhì)性、各向性、孔隙裂隙、流體密度等因素影響，地下水流情況又可分為穩(wěn)定流態(tài)和非穩(wěn)定流態(tài)水流、飽和流態(tài)和非飽和流態(tài)水流。軟件模型由一維的點(diǎn)式模型逐漸發(fā)展成二維平面模型，隨著科學(xué)計(jì)算機(jī)的進(jìn)步和發(fā)展逐漸形成了現(xiàn)在適用性最好的三維立體模型。邊界范圍和初始條件也是軟件進(jìn)行模擬分析范圍確定的主要參數(shù)。同時(shí)為了提高軟件模擬的精確性和可適用性，還應(yīng)考慮野外現(xiàn)場(chǎng)測(cè)試以及試驗(yàn)測(cè)量數(shù)據(jù)如彌散系數(shù)、滲透系數(shù)、分配系數(shù)、反應(yīng)速率常數(shù)等[3]。

文章以太子河流域地下水為例，考慮以上對(duì)模擬范圍確定的影響因素，分析該流域的水文地質(zhì)構(gòu)造條件、地下水流特性、農(nóng)業(yè)用水以及生產(chǎn)生活用水等情況。通過(guò)分析太子河流域水文地質(zhì)、巖層結(jié)構(gòu)、水流動(dòng)態(tài)方向、地理位置等條件，制定各參數(shù)取值范圍，建立了Visual Modflow軟件模型，然后對(duì)模型的計(jì)算結(jié)果與實(shí)際檢測(cè)結(jié)果進(jìn)行對(duì)比驗(yàn)證，并結(jié)合當(dāng)?shù)厮临Y源監(jiān)管部門對(duì)地下水位的控制要求，為太子河流域地下水資源的開(kāi)采和利用提出了合理化建議[4]。

1 太子河流域地下水狀況

1.1 水文地質(zhì)條件

處于我國(guó)遼寧省東部的太子河流域其大氣環(huán)境屬于溫暖帶濕潤(rùn)半濕潤(rùn)氣候區(qū)，氣候適宜，平均溫度為5～10℃，全年降水量充足，其中71.2%降雨量集中在6～8月。太子河全長(zhǎng)413km，流經(jīng)本溪、遼陽(yáng)以及鞍山等市，流域面積約為13880km2。太子河流域上游是以落葉闊葉林為主的山地森林富水區(qū)，植被所占面積超過(guò)50%，為山地森林多水區(qū)；低山丘陵區(qū)即為太子河中游地區(qū)，所占面積較少，為6.1%；其下游區(qū)為農(nóng)作物平原區(qū)，也稱為少水區(qū)，比較適合農(nóng)作物生長(zhǎng)及城鎮(zhèn)居住，該區(qū)域面積適中，占流域總面積的24.9%。

在分析研究區(qū)域含水介質(zhì)以及水流特性的基礎(chǔ)上，將太子河流域地下水層合理地劃分為三類即第四系松散巖類孔隙潛水、白堅(jiān)系碎屑巖裂隙孔隙潛水及承壓水。含水層受地質(zhì)條件和沉積巖體變化影響，其分布規(guī)律和形成條件相差較大。綜合考慮，太子河上游區(qū)為第四系松散巖類孔隙水，該層地質(zhì)基巖以砂巖居多，其中含有少量泥質(zhì)或鈣質(zhì)膠結(jié)的泥巖和碎巖，厚度約為200～350m。地表具有很厚的松散堆積物，為地下水的形成和存儲(chǔ)提供了必要的先決條件。該區(qū)域地勢(shì)平坦，降雨量較多，有利于對(duì)降雨進(jìn)行存儲(chǔ)補(bǔ)給，是該流域富水地區(qū)，有地下水庫(kù)之稱。中游區(qū)為白堅(jiān)系碎屑巖裂隙，巖層厚度較小，約10～30m，最大不超過(guò)50m，上層含有少量泥沙，透水性較差。該區(qū)域分布狹窄，寬度較小，一般在800mm以內(nèi)，低山丘陵地面高低起伏，土質(zhì)結(jié)構(gòu)疏松，雖有一定降水量但受自然地貌環(huán)境因素制約，對(duì)水分的吸收和儲(chǔ)存能力較差，富水性差，不利于大氣環(huán)境的水氣循環(huán)。下游區(qū)域土地平緩，為承壓水層，水層厚度20～80m，水位埋深5～20m，適于農(nóng)作物生長(zhǎng)及人們生活居住。該層地下水主要被用于植被生長(zhǎng)灌溉、供應(yīng)生活飲用水、工廠生產(chǎn)化工用水。據(jù)統(tǒng)計(jì)，該流域各方面總供水量高達(dá)53.382×108m3/a，其中農(nóng)作物灌溉用水47.993×108m3/a，占用水量的90%，生活用水為4.382×108m3/a，工業(yè)用水2.074×108m3/a，林牧業(yè)用水5.558×108m3/a。

1.2 孔隙潛水運(yùn)動(dòng)特征及化學(xué)特性

由于含水層之間存在相對(duì)不穩(wěn)定的含水層，故第三層含水層與太子河潛水層在相同部位的結(jié)構(gòu)動(dòng)態(tài)特征大致一樣，且其地下水運(yùn)動(dòng)規(guī)律同樣受天氣氣候影響。中游山地丘陵地區(qū)，其降水期主要集中在3～5月，地下水位在5月份開(kāi)始逐漸上升至7月份達(dá)到極值點(diǎn)，然后開(kāi)始下降，到第二年3月份下降到最低。在少水區(qū)即平原地區(qū)，地下水不僅與降水量有關(guān)，而且與農(nóng)業(yè)灌溉以及生產(chǎn)生活用水量有關(guān)，同樣也受氣象氣候影響[5]。降雨補(bǔ)給、河流運(yùn)動(dòng)以及日照蒸發(fā)是研究地區(qū)地表水匯入和流出的重要途徑，占總徑流量的92%以上。根據(jù)大氣環(huán)境和年降水量，每年地表水流經(jīng)徑流量變化幅度不大，屬于均衡分配。結(jié)合每年水流量變化情況，當(dāng)年8月份至來(lái)年5月份是流量最大時(shí)期，流量在6月份會(huì)急劇下降，其原因可能是農(nóng)牧業(yè)使用或灌溉用水。6月中旬至7月中旬流量將至最低，8月份開(kāi)始逐漸上升至來(lái)年5月份。白堅(jiān)系碎屑巖裂隙潛水存于低山丘陵區(qū)，該地區(qū)以風(fēng)積碎石、膠泥砂以及吸粉砂為主，地表結(jié)構(gòu)質(zhì)地松散，對(duì)于降水和地下流經(jīng)方面產(chǎn)生不利影響，不利于地下水的交流循環(huán)以及鹽含量較低。因此在個(gè)別山丘低洼地區(qū)，地下水位較高，陽(yáng)光蒸發(fā)擴(kuò)散作用強(qiáng)烈，風(fēng)化作用嚴(yán)重。本地區(qū)礦物化學(xué)成分主要為HCO3- SO4- Na- Ca，礦化度大于1200mg/L，溶解性固體271.25～393.84mg/L，小于600mg/L，水質(zhì)情況一般?？傆捕仍?72.35～198.86mg/L范圍之內(nèi)。pH值在7.25～8.18區(qū)間，屬于弱堿性微硬型淡水。

2 建立Visual Modflow模型

2.1 概念模型與數(shù)學(xué)模型

2.1.1 概念模型

文章通過(guò)以上分析并結(jié)合實(shí)際情況，建立Visual Modflow模型，模型模擬的范圍為太子河流域下游區(qū)域，區(qū)域?qū)儆诤畬悠皆瓟嗔研运牡刭|(zhì)，其邊界上沒(méi)有側(cè)向地下水流入，僅有側(cè)向流出。含水層埋藏于盆地及堅(jiān)硬的沉積巖下面，形成了第一承壓含水層和第二含水層，地下水位因含水層邊界變化所引起的流量變化量相對(duì)較小。下游區(qū)域地勢(shì)平坦，從而形成了可約束的地質(zhì)邊界條件[6]。研究區(qū)域包括本溪、遼陽(yáng)等所有的灌溉區(qū)域，其中包括少量個(gè)別的鞍山灌溉渠，覆蓋面積近8600km2。下游地區(qū)為居民生活區(qū)，其地下水由北向東南方向運(yùn)動(dòng)，東南部的本溪地下水自東南向西北方向運(yùn)動(dòng)，最終均匯入太子河。太子河地下水源匯流包括：氣候降雨滲流、生活用水排放、農(nóng)業(yè)灌溉匯入、冰雪融化流入、地下深水供給、日照蒸發(fā)、植被吸收、人工開(kāi)采利用等，其中氣候降雨、人工開(kāi)采利用以及農(nóng)業(yè)灌溉是地下水匯流的主要影響因素。降水變化引用系數(shù)及降水匯流系數(shù)計(jì)算結(jié)果見(jiàn)表1。

表1 降水匯流系數(shù)計(jì)算表

通過(guò)大量的查閱資料并對(duì)水文地質(zhì)進(jìn)行現(xiàn)場(chǎng)勘探觀測(cè)，將水文地質(zhì)邊界條件概化為以南部、西部以及北部作為流入邊界，因?yàn)榇诉吔缢鲄R入相對(duì)較小，結(jié)構(gòu)層相對(duì)致密[7]。東部為地下水流流出邊界，東北部可視為零流量邊界。不同含水層之間存在的水頭差相互作用，產(chǎn)生不同的影響是將其作為一個(gè)水系統(tǒng)的關(guān)鍵因素。研究區(qū)內(nèi)存在斷裂層，斷裂層與含水層之間相互貫通、相互作用，因此將兩者概化為一個(gè)含水層進(jìn)行模擬分析。由于Visual Modflow無(wú)法在一個(gè)界面上做一對(duì)多層編輯，故將南區(qū)含水層劃分為兩層，孔隙潛水層和承壓含水層之間是存在密切關(guān)系的含水層。

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

文章在充分考慮了研究區(qū)域邊界概念模型后，將太子河流域地下水流按照一定的順序逐級(jí)劃分為潛水層、弱透水層以及深水層。為了簡(jiǎn)化計(jì)算，減少計(jì)算工作量，將水文地質(zhì)條件進(jìn)行了簡(jiǎn)化和概化，并對(duì)不同含水層的參數(shù)進(jìn)行了分區(qū)。根據(jù)軟件的計(jì)算原理和方法，并采用數(shù)學(xué)計(jì)算模型，則簡(jiǎn)化后的地下水流數(shù)學(xué)方程的建立和求解過(guò)程如式(1)所示。

(1)

式中，h—水頭；K—滲透系數(shù)；*μ—含水層的單位儲(chǔ)水系數(shù)；G—輔助計(jì)算項(xiàng)；W—地下水流排泄項(xiàng)；h0—初始水頭；h1(x，y，z，t)—第一類邊界上的水頭；f1—第二類邊界上的水流通量；n—邊界外法線方向；Γ1，Γ2—第一、第二類邊界。

(2)

式中，n—有效孔隙度；C(x，y，z，t)—污染濃度；Ce—污染源的濃度；D—彌散系數(shù)張量。

上述公式可用如下矩陣表示：

(3)

Visual Modflow軟件是通過(guò)將含水層模擬成像形成一個(gè)三維的網(wǎng)格圖形，將研究區(qū)域的地下含水層按照一定的參數(shù)設(shè)置并結(jié)合地質(zhì)條件劃分為若干個(gè)小的含水層，每個(gè)小的含水層又根據(jù)坐標(biāo)和數(shù)值大小劃分為行和列。因此含水層通過(guò)參數(shù)設(shè)置計(jì)算被劃分為多個(gè)小長(zhǎng)方體，其中心位置被稱為格點(diǎn)，格點(diǎn)的多少代表被劃分的多少。然后將劃分后的格點(diǎn)通過(guò)有限元差值分析法求解地下水流方程，求解方程并建立該區(qū)域的地下水模型。為更好地標(biāo)記和計(jì)算地下水流方程，用字母i代表格點(diǎn)所在位置的行號(hào)，字母j代表所在位置的列號(hào)，長(zhǎng)方體的的厚度用k表示。研究規(guī)定自上而下第一層為模型的頂層即k=1，隨深度的降低k值逐漸增大，并用x和y方向代表行和列的方向[8]。概念模型如圖1所示。

圖1 Visual Modflow含水層三維單元格劃分

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

2.2.1 研究區(qū)域的剖分

根據(jù)太子河流域的真實(shí)情況，并考慮了該地區(qū)的復(fù)雜性，文章建立的太子河流域下游地下水水文地質(zhì)概念模型將地下水流視為三維非穩(wěn)定流態(tài)，并將其地下含水層概化為單層非均質(zhì)各向異性潛水含水層[9]。采用有限差分求解，并采用強(qiáng)隱式法聯(lián)立迭代求解代數(shù)方程組，在垂直空間上，采用X、Y正交網(wǎng)格對(duì)模擬區(qū)域進(jìn)行平面上的等間距劃分，△X=△Y=1km。將太子河流域下游平原區(qū)劃分為210×180的矩形網(wǎng)格單元，因此位于計(jì)算區(qū)內(nèi)的單元有3614個(gè)在垂直方向上，考慮含水層結(jié)構(gòu)特性和概化措施，將含水層劃分為2層，每層有效計(jì)算單元為3076個(gè)，所以該研究區(qū)域內(nèi)有效單元計(jì)算格點(diǎn)總共為7228個(gè)。

2.2.2 計(jì)算的時(shí)限與步長(zhǎng)

受研究區(qū)域內(nèi)水文地質(zhì)歷史資料完整性、人工開(kāi)采數(shù)據(jù)真實(shí)性以及統(tǒng)計(jì)科學(xué)的發(fā)展局限性，并考慮本模型研究的使用性范圍和精確性，模型模擬期選擇2000年1月～2005年12月，共72個(gè)月，采用簡(jiǎn)介識(shí)別法驗(yàn)證模型中滲透參數(shù)和彈性儲(chǔ)水率。為降低模型模擬與實(shí)測(cè)水頭之間的誤差，文章選擇地下水流變化和降水較小的時(shí)間段作為模型識(shí)別期和驗(yàn)證期，故將模型識(shí)別期定為2005～2007年，模型驗(yàn)證區(qū)間為2008年。因地下水開(kāi)采利用對(duì)水位影響較大，故考慮了地下水開(kāi)采對(duì)水位變化影響，并以每季度視為一個(gè)地下水開(kāi)采期進(jìn)行模型參數(shù)輸入。

2.2.3 初始?jí)毫λ^的確定

對(duì)格點(diǎn)初始?jí)毫λ^的確定可通過(guò)計(jì)算含水層平均水頭減去相應(yīng)格點(diǎn)的上下標(biāo)高求得，不同格點(diǎn)相應(yīng)的平均總水頭值是通過(guò)模擬含水層與隔水層在格點(diǎn)處厚度的方法求得。計(jì)算含水層地下水2005年的平均水頭可參考不同含水層所在水位觀測(cè)孔的數(shù)值[10]。

2.2.4 Visual Modflow模型識(shí)別結(jié)果

為了達(dá)到模型的識(shí)別目的，文章采用間接法求各參數(shù)值。通過(guò)試錯(cuò)法給參數(shù)設(shè)定變化范圍，并不斷調(diào)整，將擬合計(jì)算的水頭和實(shí)測(cè)水頭進(jìn)行對(duì)比分析，最終以擬合結(jié)果誤差最小值作為參數(shù)初始值。分區(qū)所計(jì)算的第Ⅰ、Ⅱ承壓水層共含有25個(gè)參數(shù)區(qū)間，弱水層含有10個(gè)參數(shù)區(qū)間。通過(guò)識(shí)別、驗(yàn)證以及計(jì)算，本次模擬潛水層共有7個(gè)參數(shù)，其中水平方向的滲透系數(shù)區(qū)間為5～135m/d，垂直方向滲透系數(shù)區(qū)間為0.03～15m/d，模擬計(jì)算彈性儲(chǔ)水系數(shù)區(qū)間為0.21～0.28。弱水層橫向滲透系數(shù)區(qū)間為0.03～1.8m/d，弱水層豎向滲透系數(shù)區(qū)間為0.01～0.12m/d，彈性區(qū)間系數(shù)約為0.18。Visual Modflow模型對(duì)實(shí)測(cè)點(diǎn)進(jìn)行模擬并擬合，模擬計(jì)算水位與實(shí)測(cè)水位擬合曲線如圖2、圖3所示。

圖2 模型第Ⅰ承壓水層觀測(cè)孔水位擬合圖

圖3 模型第Ⅱ承壓水層觀測(cè)孔水位擬合圖

為更好地分析擬合情況和實(shí)際情況之間的誤差，先將統(tǒng)誤差統(tǒng)計(jì)結(jié)果進(jìn)行匯總分析，見(jiàn)表2。

表2 模擬水頭與實(shí)測(cè)水頭誤差絕對(duì)值統(tǒng)計(jì)結(jié)果

3 地下水資源預(yù)測(cè)

采取不同抽水量測(cè)試模擬方法是對(duì)太子河流域地下水進(jìn)行數(shù)值模擬時(shí)減少模擬結(jié)果偏差的重要方法。因此文章建立的Visual Modflow模型是通過(guò)對(duì)不同抽水量的不同測(cè)試孔結(jié)果進(jìn)行分析計(jì)算該含水層的含水量變化，以避免過(guò)度抽取地下水造成抽水量超出實(shí)際資源量。同時(shí)采用線性回歸方程分解法對(duì)太子河流域邊界進(jìn)行預(yù)測(cè)，考慮到該河流特征基本保持不變，所以對(duì)其他參數(shù)不進(jìn)行設(shè)置，僅考慮河流水位變化情況。一般情況下，地下水位僅受降雨量、流經(jīng)方向等自然環(huán)境影響，故不做更加深入復(fù)雜的參數(shù)設(shè)定。綜合考慮以上因素，結(jié)合當(dāng)前對(duì)太子河開(kāi)采利用現(xiàn)狀對(duì)研究區(qū)域地下水位進(jìn)行預(yù)報(bào)，其各承壓水層地下水位變化情況見(jiàn)表3。

表3 模型預(yù)測(cè)各承壓含水層地下水位變化情況 單位：m

由表3可知，在目前地下水開(kāi)采利用條件下，各承壓水層最低水位逐漸降低，地下儲(chǔ)水量逐漸減少，如果不能進(jìn)一步加強(qiáng)用水控制管理并采取合理規(guī)劃措施，未來(lái)太子河流域地下水將面臨嚴(yán)重缺水現(xiàn)象。

4 結(jié)論

文章通過(guò)對(duì)Visual Modflow軟件的研究應(yīng)用，建立了適用于太子河流域的地下水流數(shù)值預(yù)測(cè)三維模型，結(jié)果表明該模型能夠準(zhǔn)確反映太子河流域的水文地質(zhì)條件和特征；能夠客觀反映本地區(qū)實(shí)際情況，可以對(duì)該流域地下水流變化情況進(jìn)行有效模擬。該軟件操作簡(jiǎn)便、計(jì)算簡(jiǎn)單、功能性強(qiáng)、立體感好、適用性強(qiáng)，可對(duì)其他流域進(jìn)行推廣應(yīng)用，能夠及時(shí)準(zhǔn)確地對(duì)地下水流進(jìn)行預(yù)測(cè)。

文章通過(guò)Visual Modflow軟件對(duì)太子河流域地下水流的發(fā)展變化情況進(jìn)行預(yù)測(cè)，結(jié)果顯示在各承壓含水層水位逐漸降低，提出應(yīng)及時(shí)采取有效措施加強(qiáng)對(duì)地下水的開(kāi)采和管控。

[1] 涂亮, 宋漢周. 基于Visual Modflow的黑河中游地下水流數(shù)值模擬[J] . 勘察科學(xué)技術(shù), 2009(02): 19- 23.

[2] 黨素珍, 劉昌明, 王中根, 等. 近10年黑河流域上游積雪時(shí)空分布特征及變化趨勢(shì)[J]. 資源科學(xué), 2012(08): 1574- 1580.

[3] 邵景力, 崔雅莉, 趙云章, 等. 黃河下游影響帶(河南段)三維地下水流數(shù)值模擬模型及其應(yīng)用[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 2003, 33(01): 51- 55.

[4] 王巍巍. 錦州市地下水水質(zhì)現(xiàn)狀與變化趨勢(shì)研究[J]. 水利技術(shù)監(jiān)督, 2016(05): 59- 62.

[5] 王志民, 于福蘭, 劉志芬, 等. 阜新市地表水資源質(zhì)量評(píng)價(jià)[J]. 東北水利水電, 2003(09): 38- 39.

[6] 高素麗. 遼陽(yáng)市水資源開(kāi)發(fā)利用和管理保護(hù)對(duì)策[J]. 水土保持應(yīng)用技術(shù), 2011(04): 45- 47.

[7] 楊錦華. 佛寺水庫(kù)水源飲水安全與水環(huán)境保護(hù)[J]. 水利技術(shù)監(jiān)督, 2006(03): 50- 52.

[8] 溫樹(shù)影. 大伙房水庫(kù)水體富營(yíng)養(yǎng)化現(xiàn)狀分析及對(duì)策研究[J]. 水土保持應(yīng)用技術(shù), 2015(02): 31- 32.

[9] 馬宇, 王淑偉. 遼寧省水資源現(xiàn)狀分析及保護(hù)措施研究[J]. 水利規(guī)劃與設(shè)計(jì), 2015(11): 42- 44.

[10] 徐斌, 劉克巖, 時(shí)曉飛, 等. 地表水資源水質(zhì)水量結(jié)合評(píng)價(jià)方法研究[J]. 水利規(guī)劃與設(shè)計(jì), 2003(02): 25- 30.