基于Modflow的虎林灌区地下水开采三维数值模拟

韩　雪，　马大男

（黑龙江科技大学 建筑工程学院，哈尔滨 150022）

基于Modflow的虎林灌区地下水开采三维数值模拟

韩雪，马大男

（黑龙江科技大学 建筑工程学院，哈尔滨 150022）

针对地下水开采引起的地面沉降问题，以鸡西市虎林灌区为例，采用国际通用的地下水三维渗流与地面沉降耦合计算软件——Processing Modflow，建立虎林灌区地下水-地面沉降耦合模型。以2014年年均开采量为依据，对2015年的地下水位及沉降进行数值模拟计算，并通过钻孔及室内土工实验分析其沉降机理。结果表明：2015年水位呈现北部及东南部水位高、西部较低，沉降量由中心区向四周逐渐增大的趋势；沉降严重区域土样的天然含水率、密度随深度增加呈现减小的趋势，相对密度随深度变化不明显；土样的塑性指数与其相对密度呈现正相关，变形量与其压力呈正相关，土体的压缩系数和渗透系数与压力呈负相关。

地下水开采；室内土工实验；沉降机理；Processing Modflow；耦合模型；相关性

虎林灌区是三江平原的粮食主产区，为保证黑龙江省粮食安全做出了显著贡献。长期以来，由于灌区地表水灌溉设施不完善、水利工程老化情况较为严重，灌区出现地下水超采现象，导致地面沉降，对区域农业生产和地区经济的可持续发展造成了不利影响。因此，建立地下水-地面沉降水土耦合模型，研究地下水超采引起的地面沉降问题及其沉降机理，为有关部门的科学决策提供依据，显得十分必要［1-4］。

为此，笔者以虎林灌区为研究对象，应用软件Processing Modflow建立地下水-地面沉降耦合模型，并对沉降严重区域利用现场钻孔及室内土工实验进行沉降机理分析。

1　水流模型与土力学模型

1.1水流模型

对于非均质、空间三维非稳定流系统，可用地下水流的三维运动偏微分方程来描述［5］：

式中：kxx、kyy、kzz——含水层坐标轴方向的渗透系数，md-1；

h——水头值，L；

W——汇源项，d-1；

Ss——储水率，m-1；

t——时间，d；

Ω——立体计算域。

1.2土力学模型

假定地层总应力不发生变化，且土层以竖向变形为主，水平方向的应变变化微小，可忽略不计。则土体有效应力变化值，可由式（2）求得［5］:

式中：ρw——水的密度，g·cm-3；

g——重力系数，N·kg-1；

Δh——水头，L。

土体有效应力的增加或减小，使得黏性土层在垂直方向上的压缩量或回弹量呈线性增长。土层弹性变形量可通过式（3）计算:

Δb——压缩量，m；

Δt——土层压缩经历时间，s；

dA——土层厚度，m；

Ss'k——弹性骨架释水系数，压缩为正，回弹为负。

将式（3）中的qi代入到方程（1）中，给定初始条件和边界条件就可以得到地下水-地面沉降耦合模型，其方程为

初始条件：

边界条件：

式中：n——外法线单位向量；

Г0——水流区域上的第二类边界；

q——边界上的流量。

2　数值模拟

2.1模拟区边界特性及地层岩性

灌区的边界范围：北以阿布沁河为界，南至穆棱河，西边与庆丰农场接壤，东至乌苏里江沿岸的珍宝岛湿地自然保护区。模拟区地层岩性见图1。

图1　模拟区地层岩性Fig.1 Simulation area’s formation lithology

灌区地层共为五层，第一层土层为0～5 m灰黄色的亚黏土，土层厚度为5.0 m；第二层土层为5.0～7.6 m灰绿黄色的淤泥质黏土，厚度为2.6 m，为不含水土层；第三层土层为7.6～17.4 m灰黄色的粉质黏土，粉质黏土夹亚黏土，土层厚度为9.8 m，土层含水饱和；第四层土层为17.4～40.2 m灰白及灰色的粉质黏土，土层厚度为22.8 m，土层含水饱和；第五层土层为40.2～51.5 m灰绿灰黑色的淤泥质黏土，土层致密，不含水。

2.2模型边界概化与网格划分

采用Processing Modflow进行孔隙介质中地下水运动的数值模拟［6-7］，首先建立离散化的三维模型（网格剖分、层数等）［7］。根据水文地质特征和模拟计算需求，将地层概化为三层，第一层为浅部弱潜水层，位于承压含水层顶部，0～7.6 m，为亚黏土和淤泥质黏土层，厚度为7.6 m，不含水，其渗透性差，概化为隔水边界；第二层为承压含水层，7.6～40.2 m，为粉质黏土层，厚度为32.6 m，含水饱和；第三层为不透水层，位于承压含水层底部，40.2～51.5 m，为淤泥质黏土层，厚度为11.3 m，致密不含水，其渗透性差，概化为隔水边界。

由于灌区北部、南部是阿布沁河和穆棱河，西邻庆丰农场，东至乌苏里江，因此，北部、南部均为虎林灌区的分水岭，可以概化为定水头边界；西部为庆丰农场平原区，浅层地下水不与外界进行水流交换，可概化为隔水边界；东部以乌苏里江为界，常年地下水位变化不大，可视为定水头边界；中部七虎林河河流穿越，系统分析表明，七虎林河为模拟区的优势河流，可视为河流边界。虎林灌区概化图见图2。

图2　虎林灌区边界概化Fig.2 Hulin irrigation district’s generalized boundary

灌区共14口井，自上而下依次编号为1～14号。在模拟范围为479 km2平面内，将模拟区剖分为不同尺寸大小的矩形网格，一般区域的单元网格划分成500 m×500 m，井群重点区域网格加密，单元网格划分为85 m×85 m。垂直方向上剖分三层，加密后网格共分154行和156列，共剖分24 024个单元，模拟区网络剖分见图3。

2.3模型识别与验证

2.3.1模型基本参数设定

三维数值模型采用分区赋值方法，其主要参数包括水平向渗透系数、垂直渗透系数、贮水系数、有效孔隙度、降水入渗系数及灌溉入渗系数［8］。

图3　模拟区网格剖分Fig.3 Simulation area mesh

水平方向上，模拟区潜水含水层渗透系数和给水度可以分为两个区，七虎林河以北为Ⅰ区，七虎林河以南为Ⅱ区（见图4）。各分区初值参数如表1所示。

表1　虎林灌区含水层参数初值Table 1　Initial parameters of aquifer parameters in Hulin irrigation district

图4　虎林灌区含水层参数分区Fig.4　Parameter zoning map of aquifer in Hulin irrigation district

垂直方向上，浅部弱潜水层为亚黏土，渗透系数（kxx、kyy）为6 m/d，kzz为0.6 m/d，给水度0.06，弹性释水系数为0.000 27；中层承压水层为粉质黏土，渗透系数（kxx、kyy）为20.76 m/d，kzz为2.1 m/d，给水度为0.21，弹性释水系数为0.000 40；底部不透水层为淤泥质黏土，渗透系数（kxx、kyy）为5 m/d，kzz为0.5 m/d，给水度0.01，弹性释水系数为0.000 21。

2.3.2源汇项处理

灌区地下水补给来源，主要接受大气降水入渗补给、农田灌溉回渗补给及阿布沁河、穆棱河、七虎林河和乌苏里江渗漏补给。工农业及生活用水开采、浅层含水层的蒸发，越流排泄及植物蒸腾排泄是本区地下水的主要排泄方式。其中，入渗和蒸发利用recharge程序包和蒸发包处理，人工开采在模型中利用程序包well处理。

2.3.3模型验证

模拟的时间范围为2014年1月1日运行到2015年1月1日，共365个时段，时间歩长为365 d。将初始水位和其他要素输入模型后，进行水位输出，对比模拟流场和初始流场，反复调试参数，进行模型识别和校正。初始流场与实际末流场拟合图如图5所示。总体来看，模型的总体拟合程度较好。

图5　初始流场与实际末流场拟合Fig.5 Fitting of initial and actual flow field

2.4结果分析

基于Processing Modflow建立虎林灌区2015年全年的地下水流场与地面沉降耦合三维模型，如图6～8所示。

图6　含水层水位线等值线云图（单位：m）Fig.6 Contour lines of water level line（unit:m）

从图6、7可以看出，含水层在靠近乌苏里江，穆棱河和阿布沁河的水位变化比较明显，水位变化量为20 cm，中心区水位变化较初始水位有明显回升趋势。一方面，研究区北部、东部以及南部受乌苏里江和阿布沁河的入渗补给明显；另一方面含水层受到降水等面状补给，中部水位出现回升现象。平原区地下水位变化由北向南逐渐变浅，地下水位变化明显。

图7　365 d含水层水位线三维立体图（单位：m）Fig.7 365 d three-dimensional of water level（unit:m）

图8　365 d含水层沉降量等值线平面云图（单位：m）Fig.8 365 d equivalent line of water bearing stratum（unit:m）

从图8中可以看出，由于含水层补给、排泄方式不同，该层各区域的沉降量也不同，模拟区沉降量呈现由中心区向周边江河逐步增大的总体趋势，最大沉降值2.5 cm。

3　模拟区沉降机理分析

3.1综合钻孔的选取

研究模拟区浅层地下水开采，分别取不同深度含水层的原状土进行研究。利用DC-1水文地质双层抽水钻孔在沉降严重区钻孔取样。该孔位于黑龙江省农垦总局牡丹江分局八五四农场的大王家村西北，揭穿整个模拟区的浅层含水层，具有典型性和代表性。文中通过研究该孔土样的物理性指标分析虎林灌区地面沉降区的沉降机理。

3.2黏性土物理特性

此次DC-1孔内取样深度范围0～52 m，其中，5.5～5.8 m，37.5～37.8 m，49.1～49.4 m深度的土样编号分别为原001、原005、原009，实验项目均为土的基本物理参数测试、排水固结及渗透实验。为充分了解土体的物理性质，测定抽水钻孔原状土样密度、天然含水率、相对密度、液塑限四项指标，并总结实验所得的土体物理参数和压缩系数之间存在的规律。

3.2.1基本物理参数测试结果

土样基本物理参数与压缩系数随深度变化曲线如图9所示。从图9可以看出，测试的原状土土样的天然含水率在24.86%～25.71%；密度值是在1.73～1.95 g/cm3；相对密度在2.71～2.79。天然含水率值随深度增加呈现减小的趋势，同时呈现上下波动的变化特点，埋深在27.5 m以上的土样天然含水率随深度的增加而逐渐减小，埋深27.5～41.5 m之间土样的天然含水率有逐渐增大的趋势：密度随深度的增加也呈现减小的趋势，其波动并不明显，在埋深11.3 m处，密度突然变大，其值为1.95 g/cm3；其他深度值趋于平稳。相对密度随埋深增加的变化不明显，总体呈现上下波动的趋势，其值总体上在2.71～2.79之间变化。当然，随深度变化，土样的土体结构和所处应力环境对土体基本物理参数产生重要的影响。

图9　土样基本物理参数与压缩系数随深度变化曲线Fig.9 Soil sample’s basic physical parameters and compression coefficient change with depth

3.2.2液塑限实验结果

由实验结果可知，实验土样的液限（ωL）在29%左右，塑限（ωp）在20%左右，其值均随着深度的增加而逐渐增大，从而说明土体中黏粒的含量在逐步增加。从图9中可以看出，土样的塑性指数与其相对密度呈现正相关性。相对密度的大小与其矿物组成成分有关，随着土样黏粒含量的增加而增大；土样的塑性指数和相对密度又与土体的压缩系数成一定的负相关性，也就是说，塑性指数越大，相对密度越大，土体越不容易压缩。

3.3土样压缩变形规律

由实验可知，土体的初始孔隙比（e0）随着埋深的增加而呈现逐渐减小的趋势。压缩系数在深度范围11.3～21.6 m，27.0～36.0 m两处有峰值，在22.0～26.0 m，36.0～46.0 m处有两个谷值，其他深处的值变化不大，这与土样所处的地质环境和土体结构密切相关。

3.3.1变形量随压力的变化特征

不同深度土样的变形量随压力变化曲线如图10所示。由图10可以看出，土样的变形量随压力的增大而增大，在压力100～400 kPa荷载作用下，土样的变形量变化较为剧烈，当压力达到1 000 kPa时，土样的变形量趋于稳定，尤其对于原009土样，这种变化趋势更为明显。综合分析，相同结构的土体，其压缩变形量随土体埋深的增加而呈现逐渐减小的趋势。

图10　不同深度土样变形量随压力变化曲线Fig.10 Varying depth of soil sample deformation with pressure

3.3.2土体压缩系数随压力变化特征

不同深度土样的压缩系数随压力变化曲线如图11所示。由图11可以看出，土样的压缩系数随着压力的增加而呈现逐渐减小的趋势。原001和原005土样在压力范围100～250 kPa之间，波动较大，出现谷峰值。三种土样在前期压应力较小的情况下，压缩系数变化剧烈，其增大值和减小值变化幅度非常大。埋深在49.1～49.4 m的原009土样为低压缩性土，埋深在37.5～37.8 m的原005土样为中压缩性土，埋深在5.5～5.8 m的原001土样为中压缩性土。

3.3.3土体渗透系数随埋深变化特征

不同深度土样的渗透系数变化曲线如图12所示。由图12可知，0～50 m土样的渗透系数数量级为101～10-2，且随着土样埋深的增大，渗透系数逐渐减小。

图11　不同深度土样的压缩系数随压力变化曲线Fig.11 Compression coefficient of different depth soil samples with pressure change curve

图12　不同深度土样的渗透系变化曲线Fig.12 Change curve of soil samples with different depth

综合上述分析，在研究区内伴随着地下水的开采，地下水位降低，土体孔隙水压力减小，而有效应力增大，土体的变形量与其压力呈正相关性，土体的压缩系数和渗透系数与压力呈负相关性。这些结果显示，由于土样中黏性土和淤泥质黏土占的比例比较大，所以容易形成永久性地面沉降。该研究结果对揭示研究区地面沉降机理有重要意义。

4　结　论

（1）基于Processing Modflow软件建立三维地下水-地面沉降耦合模型，模拟结果显示：2015年水位呈现北部及东南部水位高、西部较低、沉降量由中心区向四周逐渐增大的趋势。

（2）沉降严重区的现场钻孔与室内实验结果表明，沉降重灾区综合钻孔土样的天然含水率、密度随深度增加呈现减小的趋势 ，相对密度随深度变化不明显。

（3）土样的压缩变形规律为：土样的塑性指数与其相对密度呈现正相关，变形量与其压力呈正相关，土体的压缩系数和渗透系数与压力呈负相关。这一规律对揭示研究区地面沉降机理有重要意义。

［1］吴昌友，李中才，王福林.三江平原挠河流域地下水流数学模型的建立及仿真［J］.水动力学研究与进展，2011，26（3）:384-389.

［2］LIN LIN，YANG JIN ZHONG，ZHANG BIN.A simplified numerical model of 3-d groundwater and solute transport at large scale area［J］.Journal of Hydrodynamics，2010，22（3）:319-328.

［3］罗美芳，陈远法，杜丽坷，等.浙江温州市永强平原地面沉降分析与趋势预测［J］.中国地质灾害与防治学报，2012，23（3）:51-56.

［4］孙香泰，艾晓燕，张力春.三江平原地下水生态水位初步研究［J］.黑龙江水利技，2012，40（8）:166-168.

［5］付延玲，郭正法.Processing Modflow在地下水渗流与地面沉降研究中的应用［J］.勘察科学技术，2006，24（4）:19-23.

［6］陶月赞，姚 梅.地下水渗流力学的发展进程与动向［J］.吉林大学学报:地球科学版，2007，37（2）:221-230.

［7］CHIANG WENHSING，WOLFGANG KINZELBACH.3D-Groundwatermodeling with PMWIN-A Simulation S-ystem f or Modeling Ground water Flow and Pollution［M］.［S.l.］：Springer，2001.

［8］路瑞利，方树星，王红雨.基于Modflow的某水源区地下水开采三维数值模拟［J］.武汉大学学报：工学版，2011，44（5）:618-623.

［9］孙袖正.地下水流的数学模型和数值方法［M］.北京:地质出版社，1981.

（编辑荀海鑫）

MODFLOW -based 3D numerical simulation of groundwater exploitation in Hulin irrigation area

HAN Xue，MA Da’nan
（School of Civil Engineering，Heilongjiang University of Science＆Technology，Harbin 150022，China）

This paper is focused on addressing ground settlement resulting from groundwater exploitation.The research exemplified by a case study of irrigation areas in Hulin in Jixi city involves using Processing Modflow-an internationally accepted coupling calculation software for 3D groundwater seepage and surface subsidence and thereby developing the coupling-model between groundwater and ground subsidence typical of Hulin irrigation area.The study describes a simulation calculation of the groundwater level and subsidence in 2015，according to the annual mean production in 2014 and an analysis of the mechanism behind the settlement by drilling and laboratory soil test.Results show that，in 2015 there is a higher water level in the north and southeast，but a lower water level in the west，with the settlement tending to have a gradual increase from the central area to surrounding areas；the natural moisture content and density tend to decrease due to an increasing depth in the soil samples taken from the area with severe subsidence，but there is no obvious change in relative density；a positive correlation exists between the plastic index and relative density of the soil sample，and between the deformation amount and the pressure，but a negative correlation occurs between the compression coefficient and the osmotic coefficient of the soil on the one hand and the pressure on the other hand.

groundwater mining；laboratory soiltestin；mechanism of settlement；Processing Modflow；coupling-model；relevance

10.3969/j.issn.2095-7262.2015.06.016

P641.8

2095-7262（2015）06-0654-06

A

2015－10－10

黑龙江省自然科学基金项目（E201462）

韩雪（1969－），男，吉林省榆树人，教授，博士，研究方向:边坡稳定技术及其理论，E-mail:898531033@qq.com。