基于Visual MODFLOW的黑龍江省巴彥縣地下水流數(shù)值模擬

婁安宇，陳 末，臧珊珊，吳志剛

(黑龍江大學(xué) a.水利電力學(xué)院; b.寒區(qū)地下水研究所, 哈爾濱 150080)

0 引 言

隨著社會的發(fā)展，巴彥縣的城鎮(zhèn)化建設(shè)也迅速發(fā)展，地表水資源已經(jīng)供不應(yīng)求，地下水資源的利用在該地區(qū)已達(dá)到不可替代的地位。不過，地下水不是無窮無盡的，所以要對該地區(qū)的地下水進(jìn)行分析研究。地下空間錯(cuò)綜復(fù)雜，地下水含水層中有著各種形狀復(fù)雜的空隙[1]。因此，分析研究區(qū)地下水的補(bǔ)給來源、排泄途徑，進(jìn)行邊界條件的概化，構(gòu)建研究區(qū)地下水流概念模型，從而得到相關(guān)的水文地質(zhì)參數(shù)是一個(gè)極為簡單的方法[2]。

數(shù)值模擬方法在地下水研究方面的效率、靈活性和低成本管理等方面已成為一種重要的方法，并日益引起注意和廣泛的應(yīng)用[3]。數(shù)值模擬方法是一種數(shù)值模擬程序的微分方程，根據(jù)邊界和初始基本條件，來描述地下水水位隨著時(shí)間的推移而產(chǎn)生連續(xù)變動(dòng)得數(shù)學(xué)模型,并且通過電腦軟件中的離散化方法來解決這個(gè)數(shù)學(xué)問題[4]。

本文應(yīng)用數(shù)值模擬方法，通過使用Visual MODFLOW 軟件建立巴彥縣地下水流數(shù)值模型，得到該地區(qū)相應(yīng)的水文地質(zhì)參數(shù)，并對這一地區(qū)地下水水位進(jìn)行預(yù)測。

1 研究區(qū)概況

巴彥縣位于黑龍江省西部偏南，隸屬于哈爾濱市，巴彥縣境南北長85 km、東西寬72.7 km，幅員3 137.7 km2。境內(nèi)主要有“一江三河”，其中“一江”是松花江，“三河”分別為泥河、少陵河和五岳河。見圖1。

圖1 研究區(qū)示意圖

巴彥縣地區(qū)屬于中溫帶大陸性季風(fēng)氣候，每個(gè)季節(jié)的溫差很大，溫度最高可以達(dá)到36.5℃，最低可以到-40.2℃，年平均氣溫約為2.9℃[5]。地下水開采量達(dá)12 734.19 m3，主要為農(nóng)業(yè)灌溉用水，部分用于鄉(xiāng)鎮(zhèn)生活用水、工業(yè)生產(chǎn)和生態(tài)環(huán)境改善等方面。全縣地下水資源開發(fā)利用量約為7 400×104m3[6]。

地下水貯存條件及分布規(guī)律為，一級階地區(qū)含水層分布均勻，厚度穩(wěn)定，堆積物松散，水量較豐富，地下水類型為階地松散層孔隙承壓水。漫灘區(qū)孔隙潛水在2月下旬至3月初為枯水期，進(jìn)入雨季后水位逐漸上升，至8-9月初上升到最高值，水位年變幅1～4.5 m。階地區(qū)孔隙承壓水在4月下旬至5月上旬為枯水期，9月末至10月初出現(xiàn)最高值，地下水動(dòng)態(tài)類型為降水入滲-徑流型[7]。

綜合全區(qū)地貌成因和形態(tài)，將巴彥縣劃分為3種成因類型，3個(gè)形態(tài)單元見表1。

表1 地貌單元分區(qū)表

2 地下水流數(shù)值模型

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

以相對完整的水文地質(zhì)單元為數(shù)值模擬計(jì)算區(qū)，盡量將計(jì)算區(qū)邊界設(shè)置在自然邊界處，或者設(shè)置在容易確定流量或地下水位的人為邊界處。本次研究區(qū)位于哈爾濱市巴彥縣，北部邊界為泥河，西部邊界為漂河和少陵河，南至松花江干流。位于松嫩平原東部，松花江中游左岸，依據(jù)地形地貌特征、地層結(jié)構(gòu)及水文地質(zhì)特征，將區(qū)內(nèi)地下水劃分為河漫灘松散層孔隙潛水、一級階地松散層孔隙承壓水和高平原(松花江二級階地)松散層孔隙承壓水。隨著系統(tǒng)的變化，三維空間參數(shù)向量的一般地下水流動(dòng)速度反映了輸入和輸出系統(tǒng)的異質(zhì)性，輸入和輸出系統(tǒng)隨時(shí)間和空間的變化而變化，來顯示地下水的流動(dòng)情況。

2.2 含水層系統(tǒng)的概化

該研究區(qū)的地下水主要是以水平運(yùn)動(dòng)為水流的運(yùn)動(dòng)方向，而垂直方向的運(yùn)動(dòng)則很少，水文地質(zhì)參數(shù)伴隨空間分布的變化而改變，這反映了地下水系統(tǒng)的均質(zhì)性。地下水系統(tǒng)輸入和輸出隨時(shí)間和空間的變化而變化，地下水流是不穩(wěn)定的[8]。綜上所述，將該模擬區(qū)的含水系統(tǒng)概化為均質(zhì)的各向同性的三維非穩(wěn)定流。

2.3 模型邊界條件的概化

地下水流向主要從泥河到松花江，泥河、大荒溝、少陵河以及漂河的流向切割了該區(qū)域的含水層并且地下水與它們有補(bǔ)排關(guān)系。當(dāng)河流或湖泊切割含水層，兩者有直接水力聯(lián)系時(shí)，可作為第一類邊界條件處理[9]，因此將泥河、大荒溝、少陵河以及漂河作為定水頭邊界。由于其他邊界和地下水流方向平行，所以把其定為地質(zhì)零通量邊界，沒有側(cè)向流量的補(bǔ)給，所以作為隔水邊界。

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

數(shù)學(xué)模型描述了一個(gè)研究領(lǐng)域即地下流動(dòng)的水流，并確定數(shù)學(xué)方程的解決方案，這些條件構(gòu)成一個(gè)真正的數(shù)學(xué)結(jié)構(gòu)問題。含水系統(tǒng)概化為均質(zhì)的各向同性的三維非穩(wěn)定流，其相應(yīng)的偏微分方程以及定解條件如下：

式中：h為水位，m；t為時(shí)間，d；W為源匯項(xiàng)，1/d；K為滲透系數(shù)，m/d；Ω為模型模擬區(qū)；Γ1為第一類邊界；Γ2為第二類邊界；n為邊界面的外法線方向；H0為地下水初始水位函數(shù)；q(x,y,z,t)為第二類邊界上已知流量函數(shù)。

3 模型的識別與驗(yàn)證

3.1 研究區(qū)的剖分

對于研究區(qū)的剖分采用有限差分方法求解，并且采用強(qiáng)隱式法聯(lián)立迭代求解代數(shù)方程組[10]，在水平方向上對潛水含水層用相互垂直的平行線對研究區(qū)進(jìn)行網(wǎng)格剖分，對研究區(qū)剖分成100×100的矩形單元格，其中有效單元計(jì)算點(diǎn)共3 164個(gè)。

3.2 計(jì)算時(shí)限與步長

識別期為2015年1月1日至2016年12月31日，驗(yàn)正期是2017年1月1日至2017年12月31日。由于地下水位受地下水開采與利用的影響，所以在對模型進(jìn)行輸入?yún)?shù)時(shí)將每個(gè)季度作為一個(gè)地下水開采期。

3.3 參數(shù)的選取

3.3.1 滲透系數(shù)K

根據(jù)巴彥縣當(dāng)?shù)氐乃牡刭|(zhì)條件以及含水層巖性，將其劃分為3個(gè)不同滲透系數(shù)值的區(qū)域，見圖2。通過查閱資料以及查看巴彥縣的水文地質(zhì)圖可以發(fā)現(xiàn)，巴彥縣的含水層巖性基本為亞黏土、粉土質(zhì)砂及粉砂等地質(zhì)；根據(jù)水文地質(zhì)手冊可以得到，所有分區(qū)的滲透系數(shù)經(jīng)驗(yàn)值的區(qū)間為0.5～2 m/d。不同的分區(qū)輸入相應(yīng)的滲透系數(shù)初始值，由于此次模型屬于各向同性，所以水平滲透系數(shù)KX、KY與垂向滲透系數(shù)KZ相同。

圖2 研究區(qū)滲透系數(shù)分區(qū)圖

3.3.2 降雨入滲系數(shù)

通過查找資料可知，該地區(qū)是含水層研究的重要組成部分，橫向滲透、邊境供應(yīng)和其他形式的補(bǔ)給作用相對較弱，且觀測井的埋深都在1～5 m之間，所以判斷巴彥縣的地下水補(bǔ)給主要為大氣降水補(bǔ)給。根據(jù)巴彥縣水文地質(zhì)條件、流域水系、行政分區(qū)現(xiàn)狀、地貌和巖性等特征，把研究區(qū)的大氣降水滲透性分為6個(gè)區(qū)。見圖3。

圖3 研究區(qū)大氣降水滲透性分區(qū)圖

3.3.3 初始水頭

此次模型中初始水頭為各個(gè)觀測井中2015年1月1日的水頭值。觀測井位置見圖4。

圖4 觀測井位置圖

3.4 模型的識別

模型采用間接求參法達(dá)到識別目的, 通過不斷調(diào)整參數(shù)誤差的方法對測井計(jì)算的初值，使其計(jì)算水頭和測量水頭符合擬合要求, 即計(jì)算水頭和實(shí)測水頭之間相差很小。若計(jì)算水頭和實(shí)測水頭相差較大，說明設(shè)定參數(shù)或邊界條件不合適，需調(diào)整，只能通過模型識別、調(diào)整、再識別、再調(diào)整的多次反復(fù)過程，一直到計(jì)算水位與實(shí)測水位誤差滿足要求為止，此時(shí)所選用的參數(shù)即為校正后的參數(shù)。觀測井校準(zhǔn)水位圖見圖5，誤差分析見表2，觀測井水頭時(shí)間圖見圖6，單個(gè)觀測井水位擬合圖見圖7(以5號觀測井為例)。

圖5 觀測井水位校準(zhǔn)圖

圖6 觀測井水頭時(shí)間曲線

圖7 5號觀測井水位擬合曲線

表2 誤差分析表

通過以上圖表可以發(fā)現(xiàn)，各個(gè)觀測井的計(jì)算水位與觀測水位變化已基本一致，除個(gè)別點(diǎn)處的水位變化步調(diào)有偏差外，大部分點(diǎn)的變化基調(diào)基本上是同步的。通過該圖可知，該模型具有良好的擬合效果，說明建立的數(shù)學(xué)模型、各類邊界條件概化、水文地質(zhì)參數(shù)的選擇、源匯條件的處理是合理的，基本反映了研究區(qū)內(nèi)地下水真實(shí)的變化特征[11]。

3.5 模型的驗(yàn)證

為使建立的數(shù)值模型更為準(zhǔn)確地描述出地下水流系統(tǒng)的真實(shí)情況，需要利用現(xiàn)有的另一個(gè)時(shí)段的地下水動(dòng)態(tài)觀測資料對已經(jīng)建立的模型進(jìn)行檢驗(yàn)，以此來進(jìn)一步證明所建立的數(shù)值模型和所得模型參數(shù)的可靠性，本次驗(yàn)證的初始時(shí)間為2017年1月1日。觀測井校準(zhǔn)水位圖見圖8，誤差分析見表3，觀測井水頭時(shí)間圖見圖9，單個(gè)觀測井水位擬合圖見圖10(以5號觀測井為例)。

圖8 觀測井水位校準(zhǔn)圖

表3 誤差分析表

圖9 觀測井水頭時(shí)間曲線

圖10 5號觀測井水位擬合曲線

從模型識別結(jié)果和擬合效果可以看出，地下水位擬合程度非常高，建立的模型可以真實(shí)地反映該地區(qū)實(shí)際的地下水情況，同時(shí)也論證其數(shù)學(xué)模型、邊界條件、水文地質(zhì)參數(shù)和源匯條件建立的正確性。

3.6 水文地質(zhì)參數(shù)

經(jīng)過模型的識別與驗(yàn)證可以得到，研究區(qū)各個(gè)分區(qū)的滲透系數(shù)和降雨入滲系數(shù)值見表4、表5。

表4 分區(qū)滲透系數(shù)表

表5 分區(qū)降雨入滲系數(shù)表

4 地下水位預(yù)測

本文建立的數(shù)值模型通過不同監(jiān)測井的結(jié)果經(jīng)過分析與計(jì)算得到含水層的含水量變化情況，并根據(jù)線性回歸方法對研究區(qū)地下水位進(jìn)行預(yù)測。大多數(shù)情況下，僅有流經(jīng)方向和降雨量等一些自然環(huán)境影響地下水位。將模型初始時(shí)間設(shè)定為2015年1月1日，不改變已經(jīng)建立的模型的任何條件與參數(shù)，對未來10年的地下水位進(jìn)行預(yù)測。研究區(qū)未來10年地下水位預(yù)計(jì)變化情況見表6,地下水位折線圖見圖11，單個(gè)觀測井地下水位折線圖見圖12-圖14。

表6 地下水位預(yù)計(jì)變化情況

圖11 地下水位預(yù)測折線圖

圖12 4號觀測井水位預(yù)計(jì)變化圖

圖13 5號觀測井水位預(yù)計(jì)變化圖

圖14 3號觀測井水位預(yù)計(jì)變化圖

根據(jù)表6、圖11總體來看，巴彥縣地下水位將會逐年下降，造成這一現(xiàn)象的根本原因是地下水開采量大于補(bǔ)給量。巴彥縣經(jīng)濟(jì)主要以農(nóng)業(yè)生產(chǎn)為主，是黑龍江省的糧食主產(chǎn)區(qū)，由于地表水的供水能力有限，需要不斷開采地下水來維持農(nóng)業(yè)的發(fā)展。地下水開采主要用途是農(nóng)業(yè)灌溉，部分用于鄉(xiāng)鎮(zhèn)生活用水、工業(yè)生產(chǎn)和生態(tài)環(huán)境改善等方面。

根據(jù)圖12-圖14可以看出，巴彥縣不同地區(qū)地下水位的變化情況并不相同，大部分地區(qū)地下水位有下降趨勢。4號井位于洼興鎮(zhèn)附近，地下水不僅用于農(nóng)業(yè)灌溉，附近的城鎮(zhèn)居民生活用水以及工業(yè)用水也需要通過地下水進(jìn)行補(bǔ)給，所以位于城鎮(zhèn)附近的地下水位下降很明顯。5號井位于村莊附近，地下水基本只用于農(nóng)業(yè)灌溉，地下水位只是略有下降。3號井位于少陵水庫和漂河附近，能夠得到足夠的地表水補(bǔ)給，這一地區(qū)的地下水位會有些上升。

針對地下水過度開采的情況，應(yīng)該嚴(yán)格控制地下水的開采量，可以通過增加工廠企業(yè)的中水回用，取締新水，逐漸降低對地下水的依賴。還可以考慮通過建設(shè)各種大中小型蓄水工程，在汛期時(shí)，對過境洪水進(jìn)行攔蓄使用。同時(shí)，相關(guān)部門應(yīng)加強(qiáng)監(jiān)測站網(wǎng)的建設(shè)，進(jìn)一步積累和完善監(jiān)測資料，提高地下水資源評價(jià)成果的可靠性和精度。

5 結(jié) 論

本文以巴彥縣為研究區(qū)域，根據(jù)前期水文地質(zhì)資料的搜集、水文地質(zhì)勘察與試驗(yàn)結(jié)果，建立了水文地質(zhì)概念模型， 并運(yùn)用Visual MODFLOW軟件進(jìn)行三維地下水流模擬，經(jīng)過模型識別與校正，建立了符合目標(biāo)的地下水流模型。得到以下結(jié)論：

1) 農(nóng)業(yè)用水主要來源是地下水，同時(shí)隨著城市規(guī)模的不斷擴(kuò)大，人口增長和經(jīng)濟(jì)建設(shè)的迅猛發(fā)展使得工業(yè)和生活用水量日益增大，區(qū)域地下水資源供需矛盾日益突出。

2) 通過建立巴彥縣地下水流數(shù)值模型，經(jīng)過模型的識別與驗(yàn)證，可以確定巴彥縣地區(qū)的滲透系數(shù)為：河漫灘地區(qū)0.8 m/d；山區(qū)1.4 m/d；山前地區(qū)1.7 m/d。