基于MODFLOW的地下水流模型前处理优化

2014-07-01邵景力崔亚莉程汤培

吉林大学学报（地球科学版） 2014年4期

韩忠，邵景力，崔亚莉，程汤培，李玲，杨程

1.中国地质大学(北京)水资源与环境学院，北京 100083 2.北京应用物理与计算数学研究所计算物理重点实验室，北京 100088 3.北京应用物理与计算数学研究所高性能计算中心，北京 100094

韩忠1，邵景力1，崔亚莉1，程汤培2,3，李玲1，杨程1

MODFLOW中现有的降水补给数据和井流数据是基于每个剖分网格输入的，导致基于MODFLOW建立大区域水流数值模型前处理的降水补给文件和井流文件存储开销过大、读取效率低。为此，通过改进现有的降水补给子程序包RCH及开发新的子程序包RAW，给出了一种基于面状补排项的数据输入新方法。改进后的MODFLOW程序通过读取RCH文件或RAW文件中每个分区的补排量数据，以及每个网格对应的分区编号，在程序内部实现了补排量向每个模型网格的分配。RAW子程序包实现了多层面状补排量的表达，可用于面状的地下水开采、农业回灌等源汇项的处理。相对于原始的源汇项数据存储方式，基于新方法建立的华北平原地下水流模型，RCH及RAW文件大小分别减小为原来的1/145和1/255，整个模型数据的读取时间的加速比为7.46。

MODFLOW；RCH；RAW；大区域地下水模拟；华北平原

0 前言

MODFLOW是由美国地质调查局于20世纪80年代开发、用于模拟三维地下水流数值模拟模型的软件，可以模拟井流、河流、排泄、蒸散和补给对非均质和复杂边界条件的水流系统的影响，是目前世界上应用最广泛的地下水模拟软件[1-2]，并被大量应用于孔隙介质地下水流模拟中[3]。模块化结构是MODFLOW程序最显著的特点，它包括一个主程序和一系列相对独立的子程序包[4]，用户可以按实际需要选用其中的某些子程序包对地下水流动进行数值模拟[5]。这种模块化结构使程序易于理解、修改，以及添加新的子程序包，使得MODFLOW的应用范围不断扩展[6]。MODFLOW所包含的子程序包可分为3种类型：点状、线状、面状。其中：井流子程序包(well package，WEL)用来处理诸如开采井和注水井等点状特征的地下水补给与排泄；降水补给子程序包(recharge package，RCH)可处理诸如降水补给、农业灌溉入渗等面状特征的地下水补给与排泄[7]。

对于模型剖分的每一个垂向柱体，MODFLOW中现有的RCH子程序包只能有一个单元可以设置补给量[8]，无法实现垂向上的多层面状补排的处理；目前，多层补排多用点状特征的WEL子程序包来实现。然而，对于一些大区域地下水流数值模型的建立，用WEL子程序包处理多层面状补排量需要逐层、逐网格定义补排量，形成巨大的WEL文件，这在建立精细网格的大区域地下水流数值模型中尤为突出[7-8]，已成为建立高精度大区域地下水流模型的瓶颈。董艳辉等[9]、徐海珍等[10]对MODFLOW主程序进行了修改，通过多次调用RCH子程序包来实现垂向上多层补排的表达。然而，这种方法需要根据不同的模型添加RCH文件，模型垂向上剖分多少层即需添加多少个RCH文件，对于一些层数较多的大区域水流数值模型的建立仍存在缺陷。王仕琴等[11]在RCH程序包中增加了一个补给选项，用三维数组存储每一层每个单元格的补排通量值，来实现垂向上的多层补排。上述2种方法在一定程度上优化了MODFLOW程序对于面状补排的表达，但是仍未解决大区域地下水流数值模型读取数据量大的不足。

基于此，本研究对MODFLOW中已有的RCH子程序包进行了改进，并开发出新的子程序包 (recharge and well，RAW)；改进后的MODFLOW程序不仅实现了垂向上对多层单元格补排量的表达，而且有效地压缩了源汇项文件，减少了大区域地下水流数值模型所要读取的数据量。

1 原理及程序实现

1.1 RCH程序包的改进

RCH程序包处理降水补给采用的是全网格赋值，即模型剖分的每个单元格都赋一个补给通量，并将每个单元格的补给通量以RCH文件的形式存储在外部存储器中。主程序调用RCH子程序包运算时，每个应力期循环都要从RCH文件里读取每个单元格的补给通量值，再根据公式

Qi, j=Ii, j·DELR·DELC

计算每个单元格的补给量Qi, j。其中：Ii, j为补给通量；DELR为单元格(i,j)的宽度；DELC为单元格(i,j)的长度。这种数据读取方式增加了内存与外部存储器的数据交换量，影响程序的运行效率，而且建模过程中数据的前处理工作也相当繁琐。

实际建模过程中，一般可将整个模拟区域划分为若干个子区域，每个子区域赋一个补排通量值，即同一子区域内剖分的所有模型网格的补排通量值均相同。因此，本研究只将每个应力期循环中每个子区域的补排通量值存储在RCH文件中；同时，在RCH文件中建立了每个模型网格与每个子区域的对应关系。由此，改进后的RCH子程序包的输入文件格式为：应力期循环前为每个模型网格对应的分区编号；进入应力期循环后为每个分区对应的补排通量值。主程序调用改进后的RCH子程序包进行运算时，首先确定每个模型网格所对应的分区编号；进入应力期循环后，从RCH文件中读取每个分区的补排通量值，并根据分区编号与每个模型网格的对应关系，在程序内部将读入的补排通量值赋给所对应的单元格。以此实现模型中每个模型网格补排通量值的输入处理。RCH子程序包的其他部分如分配内存、建立方程组、均衡计算等均未作改动。

1.2 RAW程序包

MODFLOW中，使用WEL子程序包来处理开采井和注水井等。然而，对于一些地下水开采程度较高、大面积开采地下水的区域，模拟区内开采井的数量庞大，不可能用WEL模块对开采井逐个描述。因此，实际建模时通常将开采量处理成面状补排项。这样同样会出现一个问题，即：模型剖分的每个有效单元格几乎都存在开采井，用WEL子程序包处理时需要读入十分庞大的开采量数据，严重影响模型的运行效率，并大大增加了储存空间。

实际上，在收集开采量数据时，区域地下水开采量多数情况下是以行政区(县、市)给定的[11]，适宜于将井开采量处理成面状补排的表达方式。为此，笔者基于MODFLOW开发了一个新的子程序包RAW；其程序实现类似于原有的RCH程序包，只是RAW程序通过一个层循环，以面状的形式逐层的读入每个补排层上每个单元格的补排量值，实现了模型中多层面状补排量的表达。如同对RCH子程序包中数据读取部分的修改，RAW程序包在处理面状的补给或开采时，将整个模拟区域内每个补排层划分为若干个子区域，每个子区域内剖分的所有模型网格给定相同的补排量值。因此，RAW文件同样只存储了每个模型网格对应的分区编号，以及每个应力期循环中每个子区域的补排量数据。RAW程序包的输入文件格式为：应力期循环前为每个补排层上每个单元格对应的分区编号；进入应力期循环后为每个补排层每个分区的补排量值。RAW子程序包的计算步骤(图1)如下。

MODFLOW主程序的代码本文未详细给出，以省略号代替。图1 程序流程框图Fig. 1 Program flow chart

1)进入应力期循环前，MODFLOW程序从NAM文件中读取RAW文件的代号IUNIT。若IUNIT>0，则由主程序调用GWF1RAW6ALP函数，为RAW子程序中的数组变量分配内存空间，并在该函数中确定了每个模型网格对应的分区编号。该对应关系在整个模拟过程中保持不变。

2)进入应力期循环后，调用GWF1RAW6RPSS函数，读取每个补排层位上每个分区的面状补排量数据；并根据分区编号与每个模型网格的对应关系，在程序内部将补排量数据分配到所对应的网格。该补排量数据在1个应力期内保持不变。

3)数据输入部分完成后，调用GWF1RAW6FM函数，从有限差分方程组的右端项中减去补排量，方程组建立完成后,交给求解子程序包通过迭代法对水头值求解。

4)迭代求解结束后，调用GWF1RAW6BD函数进行均衡量的计算，并打印相关的计算结果。

RAW程序包的这种数据存储方式不仅减少了人工处理数据的工作量、加快了模型数据的读取速度，而且压缩后的源汇项文件更有利于数据的修改及模型的调试。另外，RAW程序包的输入文件格式是根据实际建模前处理过程中ArcView、GIS等软件的输出数据格式确定的，简化了模型数据的前处理过程，使得大区域地下水流数值模型的建立更加高效、便捷。由于RAW子程序包实现了多层面状补排量的表达，因此该程序包可用于农业开采、地下水回灌等面状源汇项的处理，拓展了MODFLOW的模拟功能。

2 程序测试

2.1 理想模型

图2 数据读取部分(a)及整个模拟过程(b)运行时间Fig. 2 Running time for reading date (a) and the whole model (b)

拟建模型的研究区域是一个边长为7 500 m的规则正方形。模型边界西部设定为定水头，初始水头设定为20 m。模型网格剖分为3层，每层剖分500行、500列，外部源汇项只有降水补给和井流开采。模型假定了一种理想的补排情况，最上层每个单元格均存在降水补给，井流开采则发生在模型剖分的每个单元格。运行基于改进前MODFLOW程序所建的地下水流模型，测试程序各部分运行所占的时间比例，其中，数据读取部分约占整个模型运行时间的53%。可见，对数据读取部分的优化处理可以有效提高整个模型的运行效率。

运用改进后RCH和新建的RAW模块来处理降水入渗补给和多层面状开采量，运行改进后的地下水流数值模型，得到的水头值和地下水均衡项与改进前完全相同，说明改进后的程序未影响模型的运行结果。本研究根据不同的分区数(4，100，200，2 000，20 000)对改进后的程序进行了5组测试。测试结果显示：随着分区数的增加，模型数据读取部分耗时会略有增加；然而即使分区数达到20 000之多，数据读取部分耗时仍不到1 s，较之原模型程序数据读取速度提高了上百倍(表1)。从文件大小来看，改进后程序的输入文件明显减小：模型分区数为4时，存储开采量数据的输入文件由改进前的288.77 MB减小为2.14 MB，存储降水补给数据的RCH文件也由改进前的26.90 MB减小为0.72 MB；分区数增大到20 000时，存储开采量数据的输入文件减小为11.40 MB，存储降水补给数据的RCH文件减小为4.95 MB。

表1 基于不同分区数的模型运行时间

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

分区数运行时间/s数据读取部分整个模型40.17108.341000.17108.952000.18107.9620000.28110.35200000.89110.46未分区124.60236.86

2.2 百万网格模型

分别用网格数为100万、200万、400万、800万的理想模型，对新开发的RAW程序包的性能进行进一步的测试。结果见图2、表2。模型中将研究区分为200个子区域，模型剖分的每个单元格都存在开采井。对程序改进前后的建模效果进行对比分析，改进后程序在程序加速及对输入文件的压缩方面效果显著，而且随着网格数的增加，程序的优势越明显。

表2 程序改进前后开采量文件大小

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

模型网格数/万原文件/MB改后文件/MB1007275.22001679.4010.14003379.2020.18006768.6039.9

3 实例应用

3.1 地下水流模拟模型

华北平原是位于我国东部太行山以东、黄河以北、燕山以南、东至渤海的广大平原，总面积13.92×104km2，是我国重要的粮食基地和工业基地[12]。由于人类活动不断加剧，该地区水资源紧缺日趋严峻，地下水超采严重[13]，是建立大区域地下水流数值模型的典型区域，其地下水合理开发利用问题是目前乃至今后一段时间本领域研究的重点。

本研究在前人工作的基础上[14-16]，将模拟区细化为1 km×1 km网格，共剖分为656行、592列规则网格，在垂向上仍旧分为3层。模拟期由原来的2002--2003年加长至2001--2010年。模型根据研究区的行政区划分为181个子区域，如图3所示。研究区的补给项主要为降水入渗、农业回灌、侧向流入、河流入渗补给，排泄项主要为蒸发和开采。华北平原地区地下水开采强度大，到目前为止华北平原的开采井数已超过120万[17]，且地下水开采同时发生在潜水含水层及下部的深层含水层。另外，所收集到的研究区域的降水补给、开采量以及农业回灌资料都是以各个县市的量给定的。因此，模型中处理降水补给、地下水开采和农业回灌时，均采用了面状补排的表达方式，即每个行政区域给定一个面状补排数据。基于改进后的MODFLOW程序重新建立RCH文件和RAW文件，具体是：用改进后的RCH程序包实现降水补给；用新开发的RAW程序包进行对农业灌溉入渗和各层开采量的表达。

图3 华北平原以区县为单元的面状补排量分区Fig. 3 Partition map of area recharge and discharge on county scale in North China Plain

3.2 模拟结果分析

在五舟SR4803S工作站上分别运行基于改进前后程序所建立的地下水流数值模型，得到的水头值完全相同，说明改进后的MODFLOW程序可以用来建立大区域地下水流数值模型。模型中需要输入的源汇项文件主要包括三大部分：降水入渗补给文件(RCH)、蒸发文件(EVT)、面状补排文件(RAW)。其中：RAW文件由改进前的1.40 GB减小为5.64 MB，数据占用储存空间减小为原来的1/255；RCH文件也由改进前的278.33 MB减小为1.90 MB，数据占用空间减小为原来的1/145。运行结果(表3)显示：改进后程序的数据读取部分时间开销明显减少，加速倍数约为7倍。其中，存储降水补给文件(RCH)和面状开采文件(RAW)的运算时间分别仅为0.97 s和1.76 s。

表3 程序改进前后模型的运算时间及加速比

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

运算时间/s程序改进前程序改进后加速比补给程序包90.400.9793.20井流程序包606.001.76344.32数据读取部分779.10104.507.46整个模型1797.701105.601.62

基于改进后MODFLOW程序建立的华北平原地下水流数值模型，改进了传统的源汇项数据读取方式。新建的RCH及RAW文件只存储每个县市的补排量数据，将每个县市的补排量值直接输入到模型中，避免了补排量向每个模型网格的人为分配，简化了模型数据的前处理过程。此外，降水补给、井流开采以及农业回灌均采用面状补排的方式处理，这样更符合实际的地下水补排状况，提高了模拟结果的准确性。

4 结论和建议

1)对MODFLOW中现有的RCH子程序包进行了改进，并开发出新的子程序包RAW。改进后的MODFLOW程序直接输入每个子区域的补排量数据，简化了源汇项文件的处理过程，提高了模型的调试效率。模型读入数据量的大幅减少，解决了大区域地下水流数值模型数据读取部分开销大的问题，提高了程序的运行速度。相对于原始的源汇项数据存储方式，基于新方法建立的华北平原地下水流模型，RCH及RAW文件大小分别减小为原来的1/145和1/255，整个模型数据的读取时间的加速比为7.46。

2)新开发的RAW程序包实现了模型中多个层位上面状补排量的表达，可用于面状的地下水开采、农业回灌等源汇项的处理，是对MODFLOW处理面状补排量程序的重要补充。

3)蒸发文件(EVT)也是全网格赋值格式，因而数据占空间很大。如果能根据包气带岩性特征将蒸发极限深度分区，利用上述原理亦可对蒸发子程序包进行改进，则模型整个数据读取时间还将进一步缩短。

[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] 魏林宏, 束龙仓, 郝振纯. 地下水流数值模拟的研究现状和发展趋势[J]. 重庆大学学报: 自然科学版, 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] 沈媛媛, 蒋云钟, 雷晓辉,等. 地下水数值模型在中国的应用现状及发展趋势[J]. 中国水利水电科学研究院学报, 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：一套水文地质学实用计算软件[J]. 广西地质, 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] 吴剑锋, 朱学愚. 由MODFLOW浅谈地下水流数值模拟软件的发展趋势[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] 王浩, 陆垂裕, 秦大庸,等. 地下水数值计算与应用研究进展综述[J]. 地学前缘, 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] 姬亚东, 柴学周, 刘其声,等. 大区域地下水流数值模拟研究现状及存在问题[J]. 煤田地质与勘探, 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].北京:中国地质大学,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] 张兆吉, 费宇红, 郭春艳,等.华北平原区域地下水污染评价[J].吉林大学学报：地球科学版,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].地球科学与环境学报,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] 邵景力, 赵宗壮, 崔亚莉,等.华北平原地下水流模拟及地下水资源评价[J]. 资源科学, 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] 王仕琴, 邵景力, 宋献方,等. 地下水模型MO-DFLOW和GIS在华北平原地下水资源评价中的应用[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”计划项目(2010CB428804)

韩忠(1988--)，男，硕士，主要从事地下水数值模拟方面的研究，E-mail:hanzhonggl@163.com

邵景力(1959--)，男，教授，主要从事水文学及水资源专业的教学和科研工作，E-mail:jshao@cugb.edu.cn。

10.13278/j.cnki.jjuese.201404207.

10.13278/j.cnki.jjuese.201404207

P641.1

韩忠，邵景力，崔亚莉，等. 基于MODFLOW的地下水流模型前处理优化.吉林大学学报:地球科学版,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.