模型嵌套技术在济南趵突泉泉域数值模型中的应用

齐欢，秦品瑞，赵振华，丁冠涛

(山东省地矿工程勘察院，山东 济南 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.