模型嵌套技術在濟南趵突泉泉域數值模型中的應用

齊歡，秦品瑞，趙振華，丁冠濤

(山東省地礦工程勘察院，山東 濟南 250014)

?

模型嵌套技術在濟南趵突泉泉域數值模型中的應用

齊歡，秦品瑞，趙振華，丁冠濤

(山東省地礦工程勘察院，山東 濟南 250014)

在地下水數值模型中，經常需要對某一重點區進行詳細刻畫，模型嵌套技術僅需對原有模型進行較小的修改，并不需要重新進行網格剖分，就可以達到網格加密的目的。該文在系統分析濟南趵突泉泉域水文地質條件的基礎上，采用地下水數值模擬軟件GMS10.0，建立了覆蓋趵突泉泉域的地下水流數值模型，利用更加靈活、高效的共享節點的嵌套技術，對濟南西郊重點區域進行加密，并與未加密、全局加密的模擬結果進行對比。結果表明：LGR在提高加密區和非加密區模擬精度的同時，大大縮減了模型的運行時間，并減小了磁盤空間的占用。隨著濟南市地表水轉化地下水工程的建設，LGR可以應用于對斷裂、補源河流、回灌井等的精確刻畫，更好的為社會、工程服務。

模型嵌套(LGR)；GMS；數值模型；趵突泉；濟南市

0 引言

在地下水數值模型中，經常需要對某一重點區進行加密，來精確的刻畫該區域的水頭變化，污染運移等情況。通常需要精確刻畫的區域有：短距離內水力梯度突變的區域，如抽水井、注水井、河流、溝渠和集中強度的面狀補給處等；污染源附近污染羽的運動模擬；突變性地層的細化：如斷層、斷裂造成的巖性位移、夾層、透鏡體、尖滅等[1-2]。

國內外很多專家學者對加密技術進行了研究。1999年，孫雪梅提出了基于有限元計算的局部網格加密的理論，編制出了非常適合在大型有限元程序設計中應用的計算機軟件[3]。2005年，Steffen Mehl和Mary C.Hill開發出了采用迭代耦合的共享節點的局部網格加密技術，分析該技術的優缺點并對其加密的理論進行了詳細論述[1]。2007年，針對無法一步到位地確定剖分方案來準確分析混凝土裂縫的問題，韋未提出一種基于有限元的局部網格加密技術，該方法可以較大程度的減少網格尺寸對模擬精度的影響，并能大大減少運行時間[4]。2012年，Vilhelmsen Troels N和Christensen Steen對比了常用加密技術的加密原理，并分析了各自的優缺點[5-9]。2013年，楊楊對共享節點的局部網格技術進行理論分析和源程序的修改，使得該技術可以對任意形狀的不規則區域進行局部網格加密[10]。

該文在系統分析和對比幾種常用加密技術的優缺點的基礎上，采用GMS10.0軟件，建立了2012—2014年覆蓋趵突泉泉域的地下水流數值模型，并采用更加靈活、高效的共享節點的模型嵌套技術(LGR)，對濟南西郊進行模型嵌套處理，對粗網格、全局加密、模型嵌套3種網格剖分的模型運行結果進行了分析對比。

1 模型嵌套技術簡介

網格加密泛指能對已有的網格尺寸進行更細的剖分，從而達到更高模擬精度的過程。主要方法有全局網格加密和重建網格剖分重點區模型兩種。全局網格加密就是重新對整個模型區域進行更細的網格剖分。該方法簡單、有效，但是當模型相對復雜、點較多時，全局加密工作量較大，不僅延長了模型的運行時間，在非重點區可能出現某些網格的疏干或者溢出等問題，從而影響到整個模型的運行。如果邊界情況相對簡單，直接進行加密，模擬精度較高。但在實際應用中，地下水數值模型復雜，對邊界重新進行確定、識別和驗證較困難。

模型嵌套技術僅需對原有模型進行較小的修改，不需要重新進行網格剖分。因此只會在有限的區域內增加計算模型的單元節點數，不僅提高了加密區的模擬精度，而且減輕了對計算機內存和儲存磁盤空間的依賴，縮短了模型運行時間[3]。一般來說，局部加密對自身的邊界沒有太多的限制，它可以利用周圍的非加密區為加密區提供邊界條件，和重新建立模型來直接加密相比實用性更強，有時甚至能夠更好的反映加密區的實際情況。

在地下水數值模擬中，目前技術比較成熟，應用較廣的加密方式主要有可變網格加密和伸縮式網格加密[6-8]。

1.1 可變網格加密

通過使用可變間距的有限差分網格，使需要加密區域的網格間距縮小，越遠離加密區，網格間距越大(圖1)。和全局加密相比，縮短了運行時間，但要求加密區每個方向的網格線都要延伸到母模型的邊界，導致加密區的外圍也被加密。不僅引入了額外的節點，導致更長的運行時間，而且新生成的存在較大縱橫比的網格也會導致數值上的誤差。

圖1 可變網格加密(黑色框內區域為加密區，黑色框外區域為非加密區)

1.2 伸縮式網格加密

首先對非加密進行模擬，然后將非加密的模擬水位或者流量作為加密區的邊界條件，對加密區進行模擬(圖2)。這種方法簡單、靈活、而且比較容易實現，但也存在較大的局限性。由于僅存在粗網格到加密區單方向信息交流，缺少從加密區到粗網格的信息反饋，所以模型運行后，在共享邊界處會產生顯著的水頭或者流量誤差。這種誤差會沿著共享邊界逐漸向加密區的內部傳遞，從而影響加密區的精度，而且誤差傳播的范圍也是未知和多變的。LGR(The shared node Local Grid Refinement)也稱通過迭代耦合共享節點的局部網格加密，該文中特指Steffen W. Mehl和Mary C. Hill以傳統的伸縮式網格加密技術為基礎，開發出的基于MODFLOW的程序包[1]。

圖2 伸縮式網格加密(黑色框內區域為加密區，黑色框外區域為非加密區)

LGR將模型分為子模型(加密區)和母模型(非加密區)，對母模型進行全區域求解，得到子模型邊界的水位(圖3)。采用cage-shell插值法對共享邊界進行插值，得到子模型的定水頭邊界。對子模型進行模擬，確定并計算邊界流量并反饋給母模型。當各個子模型同時滿足收斂標準時(迭代次數大于LGR控制文件設置的迭代最大值，或者水頭變化小于收斂標準，或者流量變化小于收斂標準時)，開始下一個應力期的模擬[1]。

圖3 LGR加密示意圖

由于LGR加密方法需要母模型提供給子模型水頭邊界，而母模型單元格的水頭值取決于該單元格中心節點處的水頭，為了將該節點的母模型的水頭值傳遞給子模型，就需要子模型共享邊界處的單元格的節點與母模型共享邊界處的單元格的節點在該處重合，這也就是“共享節點”名稱的由來。這就決定了對母模型的加密比只能是1∶1，3∶1，5∶1這樣的奇數，同時也決定了在水平方向上橫向和縱向的加密比必須相同。在垂直方向加密時，對于母模型的不同層，加密比卻可以不相同。

2 模型建立

2.1 水文地質概念模型

地下水流場模擬范圍的確定，應該以研究區水文地質條件為依據，同時還應充分考慮地下水系統的完整性和獨立性。該次工作的范圍為：東至東塢斷裂；西至馬山斷裂；北側以一系列近EW向隔水斷裂和奧陶系灰巖頂板埋深600m一線為界；南至泰山巖群變質巖[11-13]，面積1315km2(圖4)。

圖4 地下水流場模擬范圍

研究區含水系統主要包括第四系孔隙含水層和裂隙巖溶含水層。在該次模型中，將研究區分為3層，概化為非均質各向異性的承壓三維非穩定流[14-17]。

第一層為潛水含水層，該目的層主要為第四系全新統及上更新統地層，含水層巖性為中粗砂、砂礫石層，概化為二維流；第二層為越流層，為潛水含水層底板以下的粘土層，石炭、二疊系等弱透水層，只考慮垂向一維流；第三層為承壓水含水層，該目的層是寒武系和奧陶系地層，根據鉆孔資料，灰巖埋藏深度，確定模型研究深度為600m左右。

2.2 地下水流數值模型

對于上述非均質各向同性、空間三維結構、非穩定地下水流系統，可用如下方程的定解問題來描述：

式中：Kx，Ky，Kz為x，y，z方向滲透系數(m/d)；h—水位標高(m)；W—含水層的源匯項(1/d)；S—含水介質的貯水率(1/m)；h0—初始水位(m)；Γ2—二類邊界；Kn—邊界面法線方向的滲透系數(m/d)；n—Γ2邊界的外法線方向；q(x,y,z,t)—二類邊界上已知流量函數，流入為正，流出為負；隔水邊界為0(m/d)；Ω為滲流區域。

通過分析研究區的地質條件和水文地質條件，對各種源匯項資料進行數據分類整理，利用GMS 10軟件建立趵突泉泉域地下水數值模型。研究區剖分為400m×400m的網格，3層，125行，125列。

模擬期為2012年10月1日至2014年9月30日，其中2012年10月1日至2013年9月30日為模型的識別期，2013年10月1日至2014年9月30日為模型的驗證期，以5天作為一個應力期，每個應力期內包括若干時間步長，時間步長根據模擬計算量設定，嚴格控制每次迭代的誤差。

2.3 模型的識別和驗證

以2012年10月1日流場作為模型的初始流場，模擬期間源匯項根據實際搜集的資料，按照GMS要求的格式進行輸入，經過參數的調整，得到趵突泉、黑虎泉以及濟南西郊的幾個典型承壓水觀測井的水位過程擬合線(圖5)。可以看出，該模型可以很好的反映趵突泉泉域的泉水以及西郊的地下水流場變化趨勢，模型與實際情況具有良好的吻合性，符合研究區水文地質條件，所建模型合理可行[18-20]。

圖5 典型水位觀測孔擬合圖

3 模型嵌套的應用

該次模型嵌套范圍選擇濟南市西郊，該區域地下水長期觀測數據充足，研究程度高，便于對加密區內、外的實施效果進行評估，而且隨著濟南市“地表水轉換地下水工程”的實施，未來通過地下水流數值模型來評估補源工程對西郊水位的影響也顯得尤為重要。

通常認為全局網格加密后各個節點的水位為“真值”，為了對比的需要，在原有粗網格剖分的基礎上，對模型進行全局加密和局部加密，其中全局加密和局部加密的加密比為5∶1，局部網格加密范圍始于第1層、第41行、第9列，止于第3層、第52行、第20列(表1)。

表1 不同網格剖分對比

對于地下水流數值模型，400m×400m的網格剖分，很難精確地刻畫西郊附近的水位變化；而80m×80m的網格剖分，受計算機處理能力和GMS軟件自身的限制，390625個單元格的剖分很難實行，而后期的數據輸入和調參將更加困難。LGR網格單元數即為活動單元格數，不會造成單元格的“浪費”。加密區范圍如圖6所示，圖7～圖9為子模型、母模型以及整體運行結果圖。

圖6 加密區范圍

圖7 子模型運行結果

圖8 母模型運行結果

圖9 整體運行結果及長觀孔位置圖

從圖7、圖8可以看出，LGR運行時，子模型、母模型各自獨自運行，兩者通過共享邊界進行水位和流量的數據交換。圖9可以看出，子模型、母模型在共享邊界處的水位連續性較好。

通過加密區內3個長觀孔J65，J4，J3，非加密區3個長觀孔CX59，DY1，KD1粗網格和局部加密各個應力的水位與“真值”的對比，求取水位平均絕對誤差，對比情況如圖10、圖11所示。可以看出，同粗網格相比，采用LGR后，加密區內、外的水位模型精度都得到了提高，但無法達到全局加密的精度。加密區內J65，J4，J3的平均絕對誤差分別降低了20.77%，19.54%，15.79%；非加密區3個長觀孔CX59，DY1，KD1的平均絕對誤差分別降低了14.95%，14.86%，14.72%。

圖10 加密區內長觀孔誤差對比

圖11 非加密區內長觀孔誤差對比

利用4核3.2GHz，RAM為8G的計算機分別運行粗網格、全局加密、局部網格加密的模型，得出各自的運行時間以及運算前、后所占磁盤空間，對比結果如表2所示。可以看出：同粗網格相比，在運行時間上，全局加密時間增加了30.02倍，而LGR僅增加了2.32倍；運算前所占的磁盤空間，全局加密時間增加了5.18倍，而LGR僅增加了1.11倍；運算后所占的磁盤空間，全局加密時間增加了21.98倍，而LGR僅增加了1.36倍。

表2 不同網格剖分運行結果對比

全局加密時，單層的單元格數達到近40000個，受計算機運算能力和GMS軟件自身的限制，已經很難對網格的剖分、數據的輸入、文件的保存等進行操作，在實際工作中也多次出現軟件的卡死、電腦重啟等現象。隨著模型模擬范圍的擴大、層數的增加、網格尺寸的減小，這種情況將更加明顯；相比較而言，LGR僅需要較少的運行時間、較低的磁盤空間占用就可以提高加密區和非加密區的模擬結果，并達到接近全局加密的精度。

4 結論

(1)在地下水數值模型中，經常需要對某一重點區進行更細的網格剖分，來精確的刻畫該區域的水頭變化，污染運移等情況。全局網格加密工作量較大，不僅延長了模型的運行時間，更細的網格剖分可能出現某些網格的疏干或者溢出等問題；對重點區進行重新建模，需要重新對模型進行識別驗證較困難。模型嵌套技術僅需對原有模型進行較小的修改，而并不需要重新進行網格剖分，就可以達到網格加密的目的。

(2)LGR比伸縮式網格加密的模擬精度更高，同可變間距的網格加密相比計算方式更加靈活、效率更高。

(3)LGR提高了加密區、非加密區的水位模擬精度，加密區內平均絕對誤差降低了約19%，非加密區平均絕對誤差降低了約14%。

(4)同全局加密相比，LGR運行時間上縮短了92.27%，運算前所占磁盤空間減少了78.63%，運算后所占磁盤空間減少了93.80%。LGR在提高加密區和非加密區模擬精度的同時，可以大大縮減模型運行時間和較小磁盤空間占用。

(5)由于LGR主要用于詳細刻畫短距離內水力梯度突變區域和地層突變區域。濟南西郊加密區水力梯度較小，因此對于提高加密區的模擬精度上效果不太明顯；但其運行時間短、磁盤空間占用低的特點卻得到很好的體現。趵突泉泉域巖溶發育，地質構造復雜，隨著濟南市地表水轉化地下水工程的建設，LGR可以應用于對斷裂、補源河流、回灌井等的精確刻畫，更好的為社會、工程服務。

致謝：該項目得到濟南市市政公用事業局《海綿城市地下水補給水模型研究》課題、山東地礦重大科技攻關項目(2012-045)的經費資助，在此一并表示感謝。

[1] S. W. Mehl and M. C. Hill. MODFLOW-2005, the US Geological Survey Modular Ground-Water Model: Documentation of shared node local grid refinement (LGR) and the boundary flow and head (BFH) package [M]. 2005.

[2] 邢玉權.基于Modflow的局部網格加密技術的應用[D].北京：中國地質大學(北京)，2012.

[3] 孫雪梅.有限元單元網格加密方法的研究[J].西安工業學院學報.1999.19(1):54-57.

[4] 韋未，李同春，牛志偉，等.局部網格加密技術在混凝土裂縫擴展模擬中的應用[J].華南農業大學學報，2007，28(4):112-116.

[5] T. N. Vilhelmsen, S. Christensen and S. W. Mehl. Evaluation of MODFLOW‐LGR in Connection with a Synthetic Regional-Scale Model[J]. GROUND WATER，2012，50(1): 118-132.

[6] J. E. Dickinson,S. C. James, S. Mehl, etc. Eddebbarh. A new ghost-node method for linking different models and initial investigations of heterogeneity and nonmatching grids [J]. Advances in Water Resources，2007，30(8):1722-1736.

[7] M. Mansour and A. Spink. Grid Refinement in Cartesian Coordinates for Groundwater Flow Models Using the Divergence Theorem and Taylor's Series[J]. GROUND WATER，2013，51(1):66-75.

[8] S. Mehl, M. C. Hill and S. A. Leake. Comparison of local grid refinement methods for MODFLOW [J]. GROUND WATER，2006，44(6):792-796.

[9] S. W. Mehl and M. C. Hill. MODFLOW-2005, The U.S. Geological Survey Modular Ground-Water Model—Documentation of the Multiple-Refined-Areas Capability of Local Grid Refinement (LGR) and the Boundary Flow and Head (BFH) Package [J]，2007.

[10] 楊楊.基于LGR的不規則模型嵌套技術[D].北京：中國地質大學(北京)，2013.

[11] 唐永光,趙輝,蕫妍,等.濟南泉域重點地段可視化三維地質模型建立[J].山東國土資源，2012,28(5):29-32.

[12] 孫斌,彭玉明,李常鎖,等.濟南巖溶水系統劃分及典型泉域水力聯系[J].山東國土資源,2016,32(10):31-34+38.

[13] 孫斌,魏月,趙振華,等.濟南泉域巖溶地下水動態開采模型研究[J].山東國土資源,2013,29(9):65-68+84.

[14] 王慶兵,段秀銘,高贊東,等.濟南巖溶泉域地下水流模擬[J].水文地質工程地質,2009,(5):53-60.

[15] 徐軍祥,邢立亭.濟南泉域巖溶水數值預報與供水保泉對策[J].地質調查與研究,2008,(3):209-213.

[16] 秦品瑞,張慧,彭凱,荊浩森.濟南泉域巖溶水系統管理模型研究[J].山東國土資源,2012,28(12):32-36.

[17] 張杰,茅樵,宋玉琴.濟南市玉符河回灌補源保泉研究[J].水利水電科技進展,2002,(3):19-20+70.

[18] 鄭紅梅,劉明柱.Visual MODFLOW在天津市地下水數值模擬中的應用[J].華北水利水電學院學報,2007,28(2):8-11.

[19] 孫納正.地下水流的數學模型和數值方法[M].北京:地質出版社,1981:42-56.

[20] 李琴,葉永根.基于GMS的北京市房山平原區地下水數值模擬研究[J].中國農村水利水電,2012,(5):1-5.

Application of Local Grid Refinement Technique in Numerical Model of Baotu Spring Area in Jinan City

QI Huan，QIN Pinrui，ZHAO Zhenhua， DING Guantao

(Shandong Geo-engineering Exploration Institute，Shandong Jinan 250014，China)

In the numerical model of groundwater, it is often necessary to refine areas of some important areast. The model nesting technique only needs to modify the original model with little modification, and does not need to re-meshing. Previous studies have tended to focus on the theory of cryptography, but lack of practical application in groundwater numerical simulation. Based on the analysis of hydrogeological conditions of the Baotu Spring area in Jinan city, a numerical model of groundwater flow is established by using the groundwater numerical simulation software GMS10.0. By using more flexible and efficient sharing nodes of local grid refinement (LGR) , the key areas of western Jinan city have been refined. Comparing with the results of refined and global refinement, it is showed that LGR improves the accuracy of water level simulation in the refined and unrefined areas. LGR reduced the average absolute error of simulated water level in refined and unrefined areas by 19% and 14%, respectively. Compared with the global refinement, the running time of LGR has been reduced by 92.27% before and after operation. The disk space is reduced by 78.63% and 93.80% respectively. LGR can greatly reduce the running time of the model and disk space occupancy, and improve the simulation precision of refinement area and unrefinement area. With the construction of projects of groundwater transformation into groundwater in Jinan city, LGR can be applied to portray the fault, artificial recharge river and recharge wells precisely and server for the community, engineering better.

Local grid refinement(LGR); GMS; numerical model; Baotu Spring; Jinan city

2016-11-28；

2017-01-05；編輯：曹麗麗 基金項目：山東地礦重大科技攻關項目2012-045，項目名稱：濟南市市政公用事業局《海綿城市地下水補給水模型研究》課題 作者簡介：齊歡(1986—)，男，河北邢臺人，助理工程師，主要從事水工環地質工作；E-mail：943314446@qq.com

P208

B

齊歡，秦品瑞，趙振華，等.模型嵌套技術在濟南趵突泉泉域數值模型中的應用[J].山東國土資源，2017,33(4)：52-58.QI Huan,QIN Pinrui，ZHAO Zhenhua，etc.Application of Local Grid Refinement Technique in Numerical Model of Baotu Spring Area in Jinan City[J].Shandong Land and Resources, 2017,33(4)：52-58.