韓 忠，邵景力，崔亞莉，程湯培，李 玲，楊 程

1.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083 2.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所計(jì)算物理重點(diǎn)實(shí)驗(yàn)室，北京 100088 3.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所高性能計(jì)算中心，北京 100094

基于MODFLOW的地下水流模型前處理優(yōu)化

韓 忠1，邵景力1，崔亞莉1，程湯培2,3，李 玲1，楊 程1

1.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083 2.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所計(jì)算物理重點(diǎn)實(shí)驗(yàn)室，北京 100088 3.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所高性能計(jì)算中心，北京 100094

MODFLOW中現(xiàn)有的降水補(bǔ)給數(shù)據(jù)和井流數(shù)據(jù)是基于每個(gè)剖分網(wǎng)格輸入的，導(dǎo)致基于MODFLOW建立大區(qū)域水流數(shù)值模型前處理的降水補(bǔ)給文件和井流文件存儲(chǔ)開銷過大、讀取效率低。為此，通過改進(jìn)現(xiàn)有的降水補(bǔ)給子程序包RCH及開發(fā)新的子程序包RAW，給出了一種基于面狀補(bǔ)排項(xiàng)的數(shù)據(jù)輸入新方法。改進(jìn)后的MODFLOW程序通過讀取RCH文件或RAW文件中每個(gè)分區(qū)的補(bǔ)排量數(shù)據(jù)，以及每個(gè)網(wǎng)格對應(yīng)的分區(qū)編號(hào)，在程序內(nèi)部實(shí)現(xiàn)了補(bǔ)排量向每個(gè)模型網(wǎng)格的分配。RAW子程序包實(shí)現(xiàn)了多層面狀補(bǔ)排量的表達(dá)，可用于面狀的地下水開采、農(nóng)業(yè)回灌等源匯項(xiàng)的處理。相對于原始的源匯項(xiàng)數(shù)據(jù)存儲(chǔ)方式，基于新方法建立的華北平原地下水流模型，RCH及RAW文件大小分別減小為原來的1/145和1/255，整個(gè)模型數(shù)據(jù)的讀取時(shí)間的加速比為7.46。

MODFLOW；RCH；RAW；大區(qū)域地下水模擬；華北平原

0 前言

MODFLOW是由美國地質(zhì)調(diào)查局于20世紀(jì)80年代開發(fā)、用于模擬三維地下水流數(shù)值模擬模型的軟件，可以模擬井流、河流、排泄、蒸散和補(bǔ)給對非均質(zhì)和復(fù)雜邊界條件的水流系統(tǒng)的影響，是目前世界上應(yīng)用最廣泛的地下水模擬軟件[1-2]，并被大量應(yīng)用于孔隙介質(zhì)地下水流模擬中[3]。模塊化結(jié)構(gòu)是MODFLOW程序最顯著的特點(diǎn)，它包括一個(gè)主程序和一系列相對獨(dú)立的子程序包[4]，用戶可以按實(shí)際需要選用其中的某些子程序包對地下水流動(dòng)進(jìn)行數(shù)值模擬[5]。這種模塊化結(jié)構(gòu)使程序易于理解、修改，以及添加新的子程序包，使得MODFLOW的應(yīng)用范圍不斷擴(kuò)展[6]。MODFLOW所包含的子程序包可分為3種類型：點(diǎn)狀、線狀、面狀。其中：井流子程序包(well package，WEL)用來處理諸如開采井和注水井等點(diǎn)狀特征的地下水補(bǔ)給與排泄；降水補(bǔ)給子程序包(recharge package，RCH)可處理諸如降水補(bǔ)給、農(nóng)業(yè)灌溉入滲等面狀特征的地下水補(bǔ)給與排泄[7]。

對于模型剖分的每一個(gè)垂向柱體，MODFLOW中現(xiàn)有的RCH子程序包只能有一個(gè)單元可以設(shè)置補(bǔ)給量[8]，無法實(shí)現(xiàn)垂向上的多層面狀補(bǔ)排的處理；目前，多層補(bǔ)排多用點(diǎn)狀特征的WEL子程序包來實(shí)現(xiàn)。 然而，對于一些大區(qū)域地下水流數(shù)值模型的建立，用WEL子程序包處理多層面狀補(bǔ)排量需要逐層、逐網(wǎng)格定義補(bǔ)排量，形成巨大的WEL文件，這在建立精細(xì)網(wǎng)格的大區(qū)域地下水流數(shù)值模型中尤為突出[7-8]，已成為建立高精度大區(qū)域地下水流模型的瓶頸。董艷輝等[9]、徐海珍等[10]對MODFLOW主程序進(jìn)行了修改，通過多次調(diào)用RCH子程序包來實(shí)現(xiàn)垂向上多層補(bǔ)排的表達(dá)。然而，這種方法需要根據(jù)不同的模型添加RCH文件，模型垂向上剖分多少層即需添加多少個(gè)RCH文件，對于一些層數(shù)較多的大區(qū)域水流數(shù)值模型的建立仍存在缺陷。王仕琴等[11]在RCH程序包中增加了一個(gè)補(bǔ)給選項(xiàng)，用三維數(shù)組存儲(chǔ)每一層每個(gè)單元格的補(bǔ)排通量值，來實(shí)現(xiàn)垂向上的多層補(bǔ)排。上述2種方法在一定程度上優(yōu)化了MODFLOW程序?qū)τ诿鏍钛a(bǔ)排的表達(dá)，但是仍未解決大區(qū)域地下水流數(shù)值模型讀取數(shù)據(jù)量大的不足。

基于此，本研究對MODFLOW中已有的RCH子程序包進(jìn)行了改進(jìn)，并開發(fā)出新的子程序包 (recharge and well，RAW)；改進(jìn)后的MODFLOW程序不僅實(shí)現(xiàn)了垂向上對多層單元格補(bǔ)排量的表達(dá)，而且有效地壓縮了源匯項(xiàng)文件，減少了大區(qū)域地下水流數(shù)值模型所要讀取的數(shù)據(jù)量。

1 原理及程序?qū)崿F(xiàn)

1.1 RCH程序包的改進(jìn)

RCH程序包處理降水補(bǔ)給采用的是全網(wǎng)格賦值，即模型剖分的每個(gè)單元格都賦一個(gè)補(bǔ)給通量，并將每個(gè)單元格的補(bǔ)給通量以RCH文件的形式存儲(chǔ)在外部存儲(chǔ)器中。主程序調(diào)用RCH子程序包運(yùn)算時(shí)，每個(gè)應(yīng)力期循環(huán)都要從RCH文件里讀取每個(gè)單元格的補(bǔ)給通量值，再根據(jù)公式

Qi, j=Ii, j·DELR·DELC

計(jì)算每個(gè)單元格的補(bǔ)給量Qi, j。其中：Ii, j為補(bǔ)給通量；DELR為單元格(i,j)的寬度；DELC為單元格(i,j)的長度。這種數(shù)據(jù)讀取方式增加了內(nèi)存與外部存儲(chǔ)器的數(shù)據(jù)交換量，影響程序的運(yùn)行效率，而且建模過程中數(shù)據(jù)的前處理工作也相當(dāng)繁瑣。

實(shí)際建模過程中，一般可將整個(gè)模擬區(qū)域劃分為若干個(gè)子區(qū)域，每個(gè)子區(qū)域賦一個(gè)補(bǔ)排通量值，即同一子區(qū)域內(nèi)剖分的所有模型網(wǎng)格的補(bǔ)排通量值均相同。因此，本研究只將每個(gè)應(yīng)力期循環(huán)中每個(gè)子區(qū)域的補(bǔ)排通量值存儲(chǔ)在RCH文件中；同時(shí)，在RCH文件中建立了每個(gè)模型網(wǎng)格與每個(gè)子區(qū)域的對應(yīng)關(guān)系。由此，改進(jìn)后的RCH子程序包的輸入文件格式為：應(yīng)力期循環(huán)前為每個(gè)模型網(wǎng)格對應(yīng)的分區(qū)編號(hào)；進(jìn)入應(yīng)力期循環(huán)后為每個(gè)分區(qū)對應(yīng)的補(bǔ)排通量值。主程序調(diào)用改進(jìn)后的RCH子程序包進(jìn)行運(yùn)算時(shí)，首先確定每個(gè)模型網(wǎng)格所對應(yīng)的分區(qū)編號(hào)；進(jìn)入應(yīng)力期循環(huán)后，從RCH文件中讀取每個(gè)分區(qū)的補(bǔ)排通量值，并根據(jù)分區(qū)編號(hào)與每個(gè)模型網(wǎng)格的對應(yīng)關(guān)系，在程序內(nèi)部將讀入的補(bǔ)排通量值賦給所對應(yīng)的單元格。以此實(shí)現(xiàn)模型中每個(gè)模型網(wǎng)格補(bǔ)排通量值的輸入處理。RCH子程序包的其他部分如分配內(nèi)存、建立方程組、均衡計(jì)算等均未作改動(dòng)。

1.2 RAW程序包

MODFLOW中，使用WEL子程序包來處理開采井和注水井等。然而，對于一些地下水開采程度較高、大面積開采地下水的區(qū)域，模擬區(qū)內(nèi)開采井的數(shù)量龐大，不可能用WEL模塊對開采井逐個(gè)描述。因此，實(shí)際建模時(shí)通常將開采量處理成面狀補(bǔ)排項(xiàng)。這樣同樣會(huì)出現(xiàn)一個(gè)問題，即：模型剖分的每個(gè)有效單元格幾乎都存在開采井，用WEL子程序包處理時(shí)需要讀入十分龐大的開采量數(shù)據(jù)，嚴(yán)重影響模型的運(yùn)行效率，并大大增加了儲(chǔ)存空間。

實(shí)際上，在收集開采量數(shù)據(jù)時(shí)，區(qū)域地下水開采量多數(shù)情況下是以行政區(qū)(縣、市)給定的[11]，適宜于將井開采量處理成面狀補(bǔ)排的表達(dá)方式。為此，筆者基于MODFLOW開發(fā)了一個(gè)新的子程序包RAW；其程序?qū)崿F(xiàn)類似于原有的RCH程序包，只是RAW程序通過一個(gè)層循環(huán)，以面狀的形式逐層的讀入每個(gè)補(bǔ)排層上每個(gè)單元格的補(bǔ)排量值，實(shí)現(xiàn)了模型中多層面狀補(bǔ)排量的表達(dá)。如同對RCH子程序包中數(shù)據(jù)讀取部分的修改，RAW程序包在處理面狀的補(bǔ)給或開采時(shí)，將整個(gè)模擬區(qū)域內(nèi)每個(gè)補(bǔ)排層劃分為若干個(gè)子區(qū)域，每個(gè)子區(qū)域內(nèi)剖分的所有模型網(wǎng)格給定相同的補(bǔ)排量值。因此，RAW文件同樣只存儲(chǔ)了每個(gè)模型網(wǎng)格對應(yīng)的分區(qū)編號(hào)，以及每個(gè)應(yīng)力期循環(huán)中每個(gè)子區(qū)域的補(bǔ)排量數(shù)據(jù)。RAW程序包的輸入文件格式為：應(yīng)力期循環(huán)前為每個(gè)補(bǔ)排層上每個(gè)單元格對應(yīng)的分區(qū)編號(hào)；進(jìn)入應(yīng)力期循環(huán)后為每個(gè)補(bǔ)排層每個(gè)分區(qū)的補(bǔ)排量值。RAW子程序包的計(jì)算步驟(圖1)如下。

MODFLOW主程序的代碼本文未詳細(xì)給出，以省略號(hào)代替。圖1 程序流程框圖Fig. 1 Program flow chart

1)進(jìn)入應(yīng)力期循環(huán)前，MODFLOW程序從NAM文件中讀取RAW文件的代號(hào)IUNIT。若IUNIT>0，則由主程序調(diào)用GWF1RAW6ALP函數(shù)，為RAW子程序中的數(shù)組變量分配內(nèi)存空間，并在該函數(shù)中確定了每個(gè)模型網(wǎng)格對應(yīng)的分區(qū)編號(hào)。該對應(yīng)關(guān)系在整個(gè)模擬過程中保持不變。

2)進(jìn)入應(yīng)力期循環(huán)后，調(diào)用GWF1RAW6RPSS函數(shù)，讀取每個(gè)補(bǔ)排層位上每個(gè)分區(qū)的面狀補(bǔ)排量數(shù)據(jù)；并根據(jù)分區(qū)編號(hào)與每個(gè)模型網(wǎng)格的對應(yīng)關(guān)系，在程序內(nèi)部將補(bǔ)排量數(shù)據(jù)分配到所對應(yīng)的網(wǎng)格。該補(bǔ)排量數(shù)據(jù)在1個(gè)應(yīng)力期內(nèi)保持不變。

3)數(shù)據(jù)輸入部分完成后，調(diào)用GWF1RAW6FM函數(shù)，從有限差分方程組的右端項(xiàng)中減去補(bǔ)排量，方程組建立完成后,交給求解子程序包通過迭代法對水頭值求解。

4)迭代求解結(jié)束后，調(diào)用GWF1RAW6BD函數(shù)進(jìn)行均衡量的計(jì)算，并打印相關(guān)的計(jì)算結(jié)果。

RAW程序包的這種數(shù)據(jù)存儲(chǔ)方式不僅減少了人工處理數(shù)據(jù)的工作量、加快了模型數(shù)據(jù)的讀取速度，而且壓縮后的源匯項(xiàng)文件更有利于數(shù)據(jù)的修改及模型的調(diào)試。另外，RAW程序包的輸入文件格式是根據(jù)實(shí)際建模前處理過程中ArcView、GIS等軟件的輸出數(shù)據(jù)格式確定的，簡化了模型數(shù)據(jù)的前處理過程，使得大區(qū)域地下水流數(shù)值模型的建立更加高效、便捷。由于RAW子程序包實(shí)現(xiàn)了多層面狀補(bǔ)排量的表達(dá)，因此該程序包可用于農(nóng)業(yè)開采、地下水回灌等面狀源匯項(xiàng)的處理，拓展了MODFLOW的模擬功能。

2 程序測試

2.1 理想模型

圖2 數(shù)據(jù)讀取部分(a)及整個(gè)模擬過程(b)運(yùn)行時(shí)間Fig. 2 Running time for reading date (a) and the whole model (b)

擬建模型的研究區(qū)域是一個(gè)邊長為7 500 m的規(guī)則正方形。模型邊界西部設(shè)定為定水頭，初始水頭設(shè)定為20 m。模型網(wǎng)格剖分為3層，每層剖分500行、500列，外部源匯項(xiàng)只有降水補(bǔ)給和井流開采。模型假定了一種理想的補(bǔ)排情況，最上層每個(gè)單元格均存在降水補(bǔ)給，井流開采則發(fā)生在模型剖分的每個(gè)單元格。運(yùn)行基于改進(jìn)前MODFLOW程序所建的地下水流模型，測試程序各部分運(yùn)行所占的時(shí)間比例，其中，數(shù)據(jù)讀取部分約占整個(gè)模型運(yùn)行時(shí)間的53%?？梢姡瑢?shù)據(jù)讀取部分的優(yōu)化處理可以有效提高整個(gè)模型的運(yùn)行效率。

運(yùn)用改進(jìn)后RCH和新建的RAW模塊來處理降水入滲補(bǔ)給和多層面狀開采量，運(yùn)行改進(jìn)后的地下水流數(shù)值模型，得到的水頭值和地下水均衡項(xiàng)與改進(jìn)前完全相同，說明改進(jìn)后的程序未影響模型的運(yùn)行結(jié)果。本研究根據(jù)不同的分區(qū)數(shù)(4，100，200，2 000，20 000)對改進(jìn)后的程序進(jìn)行了5組測試。測試結(jié)果顯示：隨著分區(qū)數(shù)的增加，模型數(shù)據(jù)讀取部分耗時(shí)會(huì)略有增加；然而即使分區(qū)數(shù)達(dá)到20 000之多，數(shù)據(jù)讀取部分耗時(shí)仍不到1 s，較之原模型程序數(shù)據(jù)讀取速度提高了上百倍(表1)。從文件大小來看，改進(jìn)后程序的輸入文件明顯減?。耗Ｐ头謪^(qū)數(shù)為4時(shí)，存儲(chǔ)開采量數(shù)據(jù)的輸入文件由改進(jìn)前的288.77 MB減小為2.14 MB，存儲(chǔ)降水補(bǔ)給數(shù)據(jù)的RCH文件也由改進(jìn)前的26.90 MB減小為0.72 MB；分區(qū)數(shù)增大到20 000時(shí)，存儲(chǔ)開采量數(shù)據(jù)的輸入文件減小為11.40 MB，存儲(chǔ)降水補(bǔ)給數(shù)據(jù)的RCH文件減小為4.95 MB。

表1 基于不同分區(qū)數(shù)的模型運(yùn)行時(shí)間

Table 1 Running time of the model based on different partition number

分區(qū)數(shù)運(yùn)行時(shí)間/s數(shù)據(jù)讀取部分整個(gè)模型40.17108.341000.17108.952000.18107.9620000.28110.35200000.89110.46未分區(qū)124.60236.86

2.2 百萬網(wǎng)格模型

分別用網(wǎng)格數(shù)為100萬、200萬、400萬、800萬的理想模型，對新開發(fā)的RAW程序包的性能進(jìn)行進(jìn)一步的測試。結(jié)果見圖2、表2。模型中將研究區(qū)分為200個(gè)子區(qū)域，模型剖分的每個(gè)單元格都存在開采井。對程序改進(jìn)前后的建模效果進(jìn)行對比分析，改進(jìn)后程序在程序加速及對輸入文件的壓縮方面效果顯著，而且隨著網(wǎng)格數(shù)的增加，程序的優(yōu)勢越明顯。

表2 程序改進(jìn)前后開采量文件大小

Table 2 Size of exploitation files before and after the program improvement

模型網(wǎng)格數(shù)/萬原文件/MB改后文件/MB1007275.22001679.4010.14003379.2020.18006768.6039.9

3 實(shí)例應(yīng)用

3.1 地下水流模擬模型

華北平原是位于我國東部太行山以東、黃河以北、燕山以南、東至渤海的廣大平原，總面積13.92×104km2，是我國重要的糧食基地和工業(yè)基地[12]。由于人類活動(dòng)不斷加劇，該地區(qū)水資源緊缺日趨嚴(yán)峻，地下水超采嚴(yán)重[13]，是建立大區(qū)域地下水流數(shù)值模型的典型區(qū)域，其地下水合理開發(fā)利用問題是目前乃至今后一段時(shí)間本領(lǐng)域研究的重點(diǎn)。

本研究在前人工作的基礎(chǔ)上[14-16]，將模擬區(qū)細(xì)化為1 km×1 km網(wǎng)格，共剖分為656行、592列規(guī)則網(wǎng)格，在垂向上仍舊分為3層。模擬期由原來的2002--2003年加長至2001--2010年。模型根據(jù)研究區(qū)的行政區(qū)劃分為181個(gè)子區(qū)域，如圖3所示。研究區(qū)的補(bǔ)給項(xiàng)主要為降水入滲、農(nóng)業(yè)回灌、側(cè)向流入、河流入滲補(bǔ)給，排泄項(xiàng)主要為蒸發(fā)和開采。華北平原地區(qū)地下水開采強(qiáng)度大，到目前為止華北平原的開采井?dāng)?shù)已超過120萬[17]，且地下水開采同時(shí)發(fā)生在潛水含水層及下部的深層含水層。另外，所收集到的研究區(qū)域的降水補(bǔ)給、開采量以及農(nóng)業(yè)回灌資料都是以各個(gè)縣市的量給定的。因此，模型中處理降水補(bǔ)給、地下水開采和農(nóng)業(yè)回灌時(shí)，均采用了面狀補(bǔ)排的表達(dá)方式，即每個(gè)行政區(qū)域給定一個(gè)面狀補(bǔ)排數(shù)據(jù)?；诟倪M(jìn)后的MODFLOW程序重新建立RCH文件和RAW文件，具體是：用改進(jìn)后的RCH程序包實(shí)現(xiàn)降水補(bǔ)給；用新開發(fā)的RAW程序包進(jìn)行對農(nóng)業(yè)灌溉入滲和各層開采量的表達(dá)。

圖3 華北平原以區(qū)縣為單元的面狀補(bǔ)排量分區(qū)Fig. 3 Partition map of area recharge and discharge on county scale in North China Plain

3.2 模擬結(jié)果分析

在五舟SR4803S工作站上分別運(yùn)行基于改進(jìn)前后程序所建立的地下水流數(shù)值模型，得到的水頭值完全相同，說明改進(jìn)后的MODFLOW程序可以用來建立大區(qū)域地下水流數(shù)值模型。模型中需要輸入的源匯項(xiàng)文件主要包括三大部分：降水入滲補(bǔ)給文件(RCH)、蒸發(fā)文件(EVT)、面狀補(bǔ)排文件(RAW)。其中：RAW文件由改進(jìn)前的1.40 GB減小為5.64 MB，數(shù)據(jù)占用儲(chǔ)存空間減小為原來的1/255；RCH文件也由改進(jìn)前的278.33 MB減小為1.90 MB，數(shù)據(jù)占用空間減小為原來的1/145。運(yùn)行結(jié)果(表3)顯示：改進(jìn)后程序的數(shù)據(jù)讀取部分時(shí)間開銷明顯減少，加速倍數(shù)約為7倍。其中，存儲(chǔ)降水補(bǔ)給文件(RCH)和面狀開采文件(RAW)的運(yùn)算時(shí)間分別僅為0.97 s和1.76 s。

表3 程序改進(jìn)前后模型的運(yùn)算時(shí)間及加速比

Table 3 Running time of the model before and after the program improvement and the time ratio for reading data

運(yùn)算時(shí)間/s程序改進(jìn)前程序改進(jìn)后加速比補(bǔ)給程序包90.400.9793.20井流程序包606.001.76344.32數(shù)據(jù)讀取部分779.10104.507.46整個(gè)模型1797.701105.601.62

基于改進(jìn)后MODFLOW程序建立的華北平原地下水流數(shù)值模型，改進(jìn)了傳統(tǒng)的源匯項(xiàng)數(shù)據(jù)讀取方式。新建的RCH及RAW文件只存儲(chǔ)每個(gè)縣市的補(bǔ)排量數(shù)據(jù)，將每個(gè)縣市的補(bǔ)排量值直接輸入到模型中，避免了補(bǔ)排量向每個(gè)模型網(wǎng)格的人為分配，簡化了模型數(shù)據(jù)的前處理過程。此外，降水補(bǔ)給、井流開采以及農(nóng)業(yè)回灌均采用面狀補(bǔ)排的方式處理，這樣更符合實(shí)際的地下水補(bǔ)排狀況，提高了模擬結(jié)果的準(zhǔn)確性。

4 結(jié)論和建議

1)對MODFLOW中現(xiàn)有的RCH子程序包進(jìn)行了改進(jìn)，并開發(fā)出新的子程序包RAW。改進(jìn)后的MODFLOW程序直接輸入每個(gè)子區(qū)域的補(bǔ)排量數(shù)據(jù)，簡化了源匯項(xiàng)文件的處理過程，提高了模型的調(diào)試效率。模型讀入數(shù)據(jù)量的大幅減少，解決了大區(qū)域地下水流數(shù)值模型數(shù)據(jù)讀取部分開銷大的問題，提高了程序的運(yùn)行速度。相對于原始的源匯項(xiàng)數(shù)據(jù)存儲(chǔ)方式，基于新方法建立的華北平原地下水流模型，RCH及RAW文件大小分別減小為原來的1/145和1/255，整個(gè)模型數(shù)據(jù)的讀取時(shí)間的加速比為7.46。

2)新開發(fā)的RAW程序包實(shí)現(xiàn)了模型中多個(gè)層位上面狀補(bǔ)排量的表達(dá)，可用于面狀的地下水開采、農(nóng)業(yè)回灌等源匯項(xiàng)的處理，是對MODFLOW處理面狀補(bǔ)排量程序的重要補(bǔ)充。

3)蒸發(fā)文件(EVT)也是全網(wǎng)格賦值格式，因而數(shù)據(jù)占空間很大。如果能根據(jù)包氣帶巖性特征將蒸發(fā)極限深度分區(qū)，利用上述原理亦可對蒸發(fā)子程序包進(jìn)行改進(jìn)，則模型整個(gè)數(shù)據(jù)讀取時(shí)間還將進(jìn)一步縮短。

[1] Michael M G, Harbaugh A W. A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model[M]. Washington: United States Government Printing Office,1988.

[2] 魏林宏, 束龍倉, 郝振純. 地下水流數(shù)值模擬的研究現(xiàn)狀和發(fā)展趨勢[J]. 重慶大學(xué)學(xué)報(bào): 自然科學(xué)版, 2000, 23(增刊1): 50-52. Wei Linhong, Shu Longcang, Hao Zhenchun. The Present Situation and Development Tendency of Groundwater Flow Numerical Simulation[J]. Journal of Chongqing University: Natural Science Edition,2000,23(Sup.1):50-52.

[3] Perkins S P,Sophocleous M.Development of a Compre-Hensive Watershed Model Applied to Study Stream Yield Under Drought Conditions[J].Ground Water, 1998, 37(3): 418-426.

[4] 沈媛媛, 蔣云鐘, 雷曉輝,等. 地下水?dāng)?shù)值模型在中國的應(yīng)用現(xiàn)狀及發(fā)展趨勢[J]. 中國水利水電科學(xué)研究院學(xué)報(bào), 2009, 7(1): 57-61. Shen Yuanyuan, Jiang Yunzhong, Lei Xiaohui, et al. Current Practice and Development Trend in Numerical Modeling of Groundwater in China[J]. Journal of China Institute of Water Resources and Hydropower Research, 2009, 7(1): 57-61.

[5] 蔣亞萍, 陳余道. MODFLOW：一套水文地質(zhì)學(xué)實(shí)用計(jì)算軟件[J]. 廣西地質(zhì), 1999, 12(3): 75-78. Jiang Yaping, Chen Yudao. MODFLOW：A Practical Calculating Software in the Field of Hydrogeology[J]. Guangxi Geology, 1999, 12(3): 75-78.

[6] 吳劍鋒, 朱學(xué)愚. 由MODFLOW淺談地下水流數(shù)值模擬軟件的發(fā)展趨勢[J]. 工程勘察, 1999 (2): 12-15. Wu Jianfeng, Zhu Xueyu.Study on Development Trend of Software of Subsurface Flow Numerical Simulation from MODFLOW[J]. Geotechnical Investigation & Surveying, 1999(2): 12-15.

[7] 王浩, 陸垂裕, 秦大庸,等. 地下水?dāng)?shù)值計(jì)算與應(yīng)用研究進(jìn)展綜述[J]. 地學(xué)前緣, 2010,17(6): 1-12. Wang Hao, Lu Chuiyu, Qin Dayong, et al. Advances in Method and Application of Groundwater Numerical Simulation[J]. Earth Science Frontiers, 2010, 17(6): 1-12.

[8] 姬亞東, 柴學(xué)周, 劉其聲,等. 大區(qū)域地下水流數(shù)值模擬研究現(xiàn)狀及存在問題[J]. 煤田地質(zhì)與勘探, 2009, 37(5): 32-36. Ji Yadong, Chai Xuezhou, Liu Qisheng, et al. Current Situation and Existing Problems of Numerical Simulation of Groundwater Flow in a Large-Scale Area[J]. Coal Geology & Exploration，2009, 37(5): 32-36.

[9] Dong Yanhui, Lin Guomin, Xu Haizhen. An Areal Recharge and Discharge Simulating Method for MODFLOW[J]. Computers & Geosciences, 2012, 42: 203-205.

[10] 徐海珍, 李國敏, 董艷輝,等. 基于MODFLOW的面狀通量處理方法研究[J]. 工程勘察, 2012 (2): 37-46. Xu Haizhen, Li Guomin, Dong Yanhui, et al. An Approach for Simulating Planar Flux Based on MODFLOW[J]. Geotechnical Investigation & Surveying,2012(2):37-46.

[11] 王仕琴.地下水模型MODFLOW與GIS的整合研究：以華北平原為例[D].北京:中國地質(zhì)大學(xué),2006. Wang Shiqin. Study of the Coupling of Groundwater Modeling Software and GIS for North China Plain[D].Beijing: China University of Geosciences,2006.

[12] 張兆吉, 費(fèi)宇紅, 郭春艷,等.華北平原區(qū)域地下水污染評價(jià)[J].吉林大學(xué)學(xué)報(bào)：地球科學(xué)版,2012,42(5):1456-1461. Zhang Zhaoji, Fei Yuhong, Guo Chunyan, et al. Regional Groundwater Contamination Assessment in the North China Plain[J].Jouranl of Jilin University：Earth Science Edition ,2012,42(5):1456-1461.

[13] 張光輝, 連英立, 劉春華,等.華北平原水資源緊缺情勢與因源[J].地球科學(xué)與環(huán)境學(xué)報(bào),2011,33(2):172-176. Zhang Guanghui, Lian Yingli, Liu Chunhua, et al. Situation and Origin of Water Resources in Short Supply in North China Plain[J].Journal of Earth Sciences and Environment,2011,33(2):172-176.

[14] 邵景力, 趙宗壯, 崔亞莉,等.華北平原地下水流模擬及地下水資源評價(jià)[J]. 資源科學(xué), 2009, 31(3): 361-367. Shao Jingli, Zhao Zongzhuang, Cui Yali, et al. Application of Groundwater Modeling System to the Evaluation of Groundwater Resources in North China Plain[J]. Resources Science, 2009, 31(3): 361-367.

[15] 王仕琴, 邵景力, 宋獻(xiàn)方,等. 地下水模型MO-DFLOW和GIS在華北平原地下水資源評價(jià)中的應(yīng)用[J]. 地理研究, 2007, 26(5): 975-983. Wang Shiqin, Shao Jingli, Song Xianfang, et al. The Application of Groundwater Model, MODFLOW, and GIS Technology in the Dynamic Evaluation of Groundwater Resource in North China Plain[J]. Geographical Research, 2007, 26(5): 975-983.

[16] Shao Jingli, Li Ling, Cui Yali, et al. Ground-Water Flow Simulation and Its Application in Groundwater Resource Evaluation in the North China Plain, China[J]. Acta Geologica Sinica, 2013, 87(1):243-253.

[17] Liu Jie, Zheng Chunmiao, Zheng Li, et al. Ground Water Sustainability: Methodology and Application to the North China Plain[J]. Ground Water, 2008, 46(6): 897-909.

Preprocessing Optimization of Groundwater Flow Model Based on MODFLOW

Han Zhong1, Shao Jingli1, Cui Yali1, Cheng Tangpei2,3, Li Ling1,Yang Cheng1

1.School of Water Resources and Environment, China University of Geosciences, Beijing 100083, China 2.Key Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing 100088, China 3.High Performance Computing Center, Institute of Applied Physics and Computational Mathematics, Beijing 100094, China

Reading and writing precipitation and well data are all based on each grid in the existing MODFLOW. So building regional groundwater flow numerical models with refined spatial and temporal discretization often involves storing and reading large quantities of well and precipitation data. However, the storage requirement can be prohibitive and the program efficiency can be quite low. To address this problem, we have modified the existing package RCH and also developed a new package named RAW (recharge and well), which is based on the existing packages in MODFLOW. The improved MODFLOW program can only read values of every subarea and the subarea number of each grid which is stored in the RCH or RAW file, and then allocates the values to corresponding grids. With the RAW package, the area recharge and discharge could occur in all layers. Therefore, it could be used to simulate the source and sink term, like groundwater exploitation and agricultural irrigation return. The correctness has been verified and the efficiency has been tested through simulating a theoretical groundwater flow model and a real North China Plain model respectively. Compared with the traditional method of storing data, the size of RCH and RAW file decreases 1/145 and 1/255 respectively, and the time ratio for reading data is 7.46.

MODFLOW；RCH；RAW；large-scale groundwater simulation；North China Plain

2013-10-05

國家“973”計(jì)劃項(xiàng)目(2010CB428804)

韓忠(1988--)，男，碩士，主要從事地下水?dāng)?shù)值模擬方面的研究，E-mail:hanzhonggl@163.com

邵景力(1959--)，男，教授，主要從事水文學(xué)及水資源專業(yè)的教學(xué)和科研工作，E-mail:jshao@cugb.edu.cn。

10.13278/j.cnki.jjuese.201404207.

10.13278/j.cnki.jjuese.201404207

P641.1

A

韓忠，邵景力，崔亞莉，等. 基于MODFLOW的地下水流模型前處理優(yōu)化.吉林大學(xué)學(xué)報(bào):地球科學(xué)版,2014,44(4):1290-1296.

Han Zhong, Shao Jingli, Cui Yali, et al. Preprocessing Optimization of Groundwater Flow Model Based on MODFLOW.Journal of Jilin University:Earth Science Edition,2014,44(4):1290-1296.doi:10.13278/j.cnki.jjuese.201404207.