基于地质随机场孔隙分析的地下水污染扩散特征模拟研究

杨晓丽

嘉兴市环境科学研究所有限公司 浙江嘉兴 314051

摘要：为分析地下水污染扩散规律，研究基于地质随机场孔隙分析的地下水污染扩散特征模拟方法。使用德劳内三角剖分结合约束关系获得约束劳内三角剖分三角网格，在此基础上依据地质随机场孔隙分析地下水流数值模拟参数的空间分布特征，划分空间参数分布类型。通过GIS空间分析功能结合点、面空间分布参数自动提取算法，获得构建水流模型所需参数，将这些参数输入到Visual MODFLOW软件中构建地下水污染扩散特征模型，确定该模型边界条件和含水层参数，通过该模型获得最终地下水污染扩散特征模拟结果。经实验验证，该方法模拟的污染物随水流流动出现的扩散特征具有较高准确性，并分析了获得近十年的中心点浓度与最大扩散距离。

关键词：地质研究；随机场；孔隙分析；地下水；污染扩散；特征模拟；Visual MODFLOW软件

中图分类号：X824;X52 文献标识码：A

1引言

我国对于新能源的研究逐渐取得客观成果，核能发电技术是一种安全、清洁、可靠的新能源技术，建设核电站能够有效使我国的电力结构得到优化，化解煤炭能源发电造成的环境污染问题，提升发电量的同时推动环境资源可持续发展，综合各种优势，核能发电已经成为我国能源发展的重点战略工程^[1-3]。但是伴随而来的问题是核能源的安全使用与储存。众所周知，核能具有较高的辐射性和污染性，人和动物长期生活在高辐射的环境中会出现严重身体疾病和遗传疾病^[4]。内陆建设的核电站一般会毗邻湖泊、河流等大型水体，一旦出现核泄漏事故会严重污染地表和地下水，其中地下水流速比较缓慢，因此研究地下水污染扩散特征具有重要意义^[5]。核能泄漏扩散实际上是一个集合化学、物理、生物多种学科共同作用的结果。针对目前水污染扩散特征模拟有学者提出针对硝态氮扩散特征模拟方法^[6]，该方法通过软件模拟硝态氮在地下水中的扩散、迁移、衰减规律特征，该方法能够分析地下水的污染扩散规律，但是部分准确性还有待改进；还有学者提出GMS作为基础的水污染扩散特征模拟方法^[7]，该方法使用GMS技术构建研究区数学模型和概念模型分析地下水污染的周期性迁移规律特征，但是该方法的模型计算过程过于复杂，在未来研究中还有很多有待改进的地方。

与人工建造的建筑存在区别，自然土体中的结构中的土性参数存在极大不确定性和离散性。随机场可以看做是一种多维的随机过程，随机内容主要是随机变量，对于二维随机场来说，随机变量分别为位置和地质参数^[8]。如果已经固定的地质参数随机场转变为位置随机函数，位置固定时随机场转变成岩土参数的随机函数。针对存在具体位置的地质参数来说，可以局部获得位置数据，而且地质中的岩石、土壤等元素参数分布存在固定规律。地下水蕴含在地质结构孔隙之中，包含饱水带与包气带孔隙中蕴含的水。利用地质随机场孔隙分析获得地下水污染扩散的特征参数，这些参数是模型构造的基础^[9]。地质随机场孔隙地下水数值模拟时一般使用三角形空间离散格网，提升模型赋值准确性和精确度，为后续污染扩散特征模拟奠定基础。

Visual MODFLOW软件是专门针对孔隙介质的地下水数值模拟软件，该软件具有强大的可视化特征和合理的菜单结构与软件支撑，能够良好模拟地下水污染扩散特征，在分析污染物扩散运动和扩散特征时能够准确模拟特征，研究出地下水的扩散趋势。

本文研究以地质随机场孔隙作为基础，提取地下水污染扩散模型的参数，构建水流模型，获取地下水污染扩散特征。

2地下水污染扩散特征模拟

2.1基于地质随机场孔隙分析参数自动提取

2.1.1有限元三角网生成

由于德劳内三角剖分（Delaunay Trianglation，DT）存在最大最小特性，在二维情况下自动避免德劳内三角剖分生成小内角的长薄单元，在生成有限元三角网格时具有绝对优势。数值模拟实际过程中，约束关系普遍存在于模拟区域的离散点中，三角剖分该模拟区域时，应该仍然保留已经存在的约束关系，这种存在约束条件的德劳内三角剖分也被称作约束德劳内三角剖分^[10]。本文构造的约束劳内三角剖分过程如下：

（1）为获得最优三角形网格，需要先完成无约束条件下模拟区域离散点的德劳内三角剖分；

（2）考虑约束条件，获得优化后的三角剖分结果。

利用点结构存储德劳内三角剖分的离散点数据结构，利用面结构存储三角形的数据结构，为使后续处理简便，使用逆时针方向记录三角形顶点。

以下为自动生成有限元三角网格具体步骤：

（1）利用逐点插入法，完成德劳内三角剖分离散点，过程如下：

①生成离散点集的凸包，该凸包中包含全部点集；

②逆时针方向排序生成的凸包点集 ，把这些点集通过连接形成有向线段；

③依次插入凸包点集内的离散点，离散点外接圆包含的三角形需要在德劳内三角形链表中找寻，删除三角形公共边，连接插入点和三角形的全部顶点，构成一个全新的三角形；

④依据LOP优化准则对新构成的三角形实行优化，并且在德劳内三角形链表中保存优化结束后的三角形；

⑤循环步骤③，离散点插入结束后停止循环。

（2）在约束条件下实现优化德劳内三角剖分：

符合一般德劳内准则基础上开展约束德劳内三角剖分，在最终的三角划分中包括含约束线段，内边界多边形不属于内边界约束，本文研究中暂且将两种情况综合称为约束德劳内三角剖分。

约束德劳内三角剖分划分描述： 表示包含点集 、 的平面域；点在区域内边界上的集合使用 表示，假如点集 、 使用德劳内三角剖分作为三角划分，而且把点集 内全部点构成的内边界约束全部去除，此时的划分为 中的内边界约束德劳内三角剖分。约束德劳内三角剖分需要同时符合约束条件和空Circle准则。

实现约束德劳内三角剖分的主要算法思路如下：

标准约束德劳内三角剖分域 ，也就是说剖分时不将点的性质考虑进去， 作为剖分结果。由于内控不包含其它离散点，因此内孔边界点集合 内一定包含内控三角形的上个顶点，根据该原则实现约束德劳内三角剖分：

①获取任意区域中孔的约束边界，以 表示；

②判断 是否已经存在于 的三角形边界中，假如结果为存在，可以直接在三角形链表内把含有 公共边的左右两个三角形都取出来；

③查找点 和点 是否属于内孔边界上的点集 中的 ，如果结果为肯定，就把这个三角形从三角形链表中删除，如果结果是否定的就保留该三角形；

④重复步骤①至步骤④，检测全部内边界约束线，完成后可停止循环。

2.1.2参数自动提取

依据地质随机场孔隙分析地下水流数值模拟参数的空间分布特征，可以把空间划分成点与面空间分布参数类型。点、面空间分布参数自动提取算法通过GIS空间分析功能实现^[11]。获得模型需要的数据文件时需要结合有限元计算过程中的参数数据文件结构。

（1）点状参数自动提取方法

开展有限元数值模拟时主要点状参数是研究区域地下水污染扩散量。研究区域的地下水扩散分布从垂直面上来看污染位置是根据含水层的空间分布确定的；从平面上来看地下水污染空间分布属于离散情况，由此可以看出地下水污染扩散具有空间离散分布特征值^[12]。自动提取点状参数的具体路线：自动生成地下水污染扩散专题图、地下水污染点与三角网单元之间空间位置的拓扑关系、提取并统计污染扩散量、组织扩散量数据文件。

①自动生成地下水污染扩散专题图

通过ArcGIS专题图制作模块在空间数据库中保存研究区域的地下水污染扩散空间分布特征信息，实际使用时直接调用这些信息利用不同图层管理水污染扩散的不同层位。

②地下水污染点与三角网单元之间空间位置的拓扑关系

地下水污染扩散以无规律状态分布在空间上，因此先需要空间确定污染点，同时分析地下水污染点与三角网单元之间空间位置关系^[13]。为了便于剖分三角网图层的地下水污染扩散位置构建全新的属性，确定那个三角网格的网格形单元中存在地下水污染位置，分析剖分三角网单元格与污染位置之间的叠加拓扑分析，可以使用铅垂线算法实现该分析。

③污染扩散量提取算法

研究区域的地下含水层确定需利用地下水污染空间分布，分布具有离散分布特征。参考经有限元数值模拟的地下水污染扩散量提取要求，需把扩散量按照有关规则在三角形单元的哪个顶点分配，针对多个扩散点三角形计算，分配时需要针对单个点，然后再叠加顶点分配到的扩散量。使用空间统计计算方法在剖分三角网的三个节点上重新分配空间离散分布的污染扩散点，这样才能提取出模型所需的参数文件。提取算法的具体流程见图1。

图1点状参数自动提取流程

Fig. 1 Automatic extraction process of point-like parameters

提取参数时先要拓扑分析地下水污染源位置和剖分三角形单元，确定二者位置关系。假设地下水污染点坐标和三角形剖分单元顶点坐标分别为 和 ，其中 ，按照逆时针方向记录三角形剖分单元的顶点。遍历比较地下水污染点的坐标 和剖分三角形单元的顶点坐标就能判断剖分三角形单元的某个顶点是否与地下水污染点的坐标 相等；通过以下关系分析地下水污染点与三角形单元的空间位置关系：检查剖分三角形的全部边，如果符合 那么污染点在 边上，如果符合 那么污染点在 边右侧，如果符合 那么污染点在 边左侧，由此说明整个地下水污染源位置都在此剖分三角型单元内部。

（2）面状参数自动提取方法

提取面状分布参数的路线为：生成等值线图层与参数分区、分析参数分区与有限元三角网的空间位置拓扑关系、提取并调用参数和模型所需的数据文件格式，等待模型调用

^[14]。参考点状参数自动提取方法获得面状参数专题图以及拓扑关系，在此基础上提取出面状参数。

生成研究区域有限元三角网剖分，将初始水头值赋值到各个离散节点上。如果研究区域不存在专题数据，离散点的初始水位需要人为给定；如果已经水位值数据给定，那么可以使用插值法提取初始水头值。

初始水头值计算流程见图2。

图2初始水头值提取流程

Fig. 2 Extraction process of initial head value

计算初始水头值时先遍历全部初始流程图，将水位等值线单元提取出来，再进一步将水位等值线单元上的点提取出来，提取全部点时遍历初始流场图，使用IDW插值依据提取出的点生成水位栅格图。提取初始水位，遍历三角网点获取三角网点单元，叠加分析生成的水位栅格，提取出该点的初始水位值^[15]。遍历全部三角网点，将初始水位点提取出来，储存相应文件格式。

参数分区提取具体流程见图3。

图3参数分区提取流程图

Fig. 3 Parameter partition extraction flow chart

通过剖分遍历获取三角网图层和三角网单元，拓扑分析参数分区图和三角单元，提取含水层的水平渗透系数给水度与有效孔隙度参数。

2.2基于Visual MODFLOW的地下水污染扩散特征模型

将上面含水层的水平渗透系数给水度与有效孔隙度参数输入到Visual MODFLOW软件中，构建地下水污染扩散特征模型，把研究区域划分成两层结构模型。模型平面剖分大小为6.5*6.5m，地下水污染物扩散过程精细刻画使用污染附近加密1*1m的网格。

本文模拟的地下水流向是从西南方向向东北方向流动，符合地形地貌趋势，依据流畅处理，把西南和东北方向定义为水头边界，与等水线方向垂直的其它两侧被处理成隔水边界，降水入渗补给和蒸发边界被定义为上边界，研究区域的降水入渗系数为0.1，平均年降水量约为1350mm。

根据经验参数和研究区域的实际抽水试验，拟合实际水位和计算后得出的水位数据，不断调整含水层的水平渗透系数给水度与有效孔隙度参数，获得最终的含水层采纳数，由此确定出研究区域含水层的水平渗透系数在0.03m/s至0.24m/s之间，给水度与有效孔隙度分别为0.2和0.15。

3实验结果

以我国某核电站为研究对象，以MATLAB仿真平台为仿真实验平台。实验所选用的该核电站建成于1988年，位于我国东南沿海附近，是中国内陆地区第一座大型核电站截止到2018年年底总发电量达到7000亿度以上。该核电站面朝大海，受海洋气候影响多雨雾，夏季常伴随强降水，而且植被覆盖率较高，受海洋气候调节年平均气温在22℃上下，平均风速2.5m/s，舒适度较高；受地形影响地表径流坡陡流短，地表径流量与地下径流量受降水影响严重，水系属于海洋水系，河流网道密集，均为入海河流，水流方向为自西南向东北方向流动。在MATLAB仿真平台输入以上该核电站的地址信息和水文信息，设置泄漏点。

在实验中设置一个抽水井模拟地下水涌水，实验时间分别设置为核泄漏6个月以及核泄漏1年，分别模拟地下水流场中污染物扩散和迁移的情况（图4）。

（a）污染物扩散情况

（b）污染物迁移模拟

图4 核泄漏6个月

Fig. 4 Six months of nuclear leakage

（a）污染物扩散情况

（b）污染物迁移模拟

图5 核泄漏一年

Fig. 5 Nuclear leakage for one year

从图4（a）中能够看出，泄漏发生初期，污染物扩散方向呈现一种无序扩散情况，并没有随水流方向形成扩散路径，随着时间的推移，污染物的扩散方向与水流场方向一致。从图4（b）能够看出，污染物随着水流场在抽水井的位置形成一个稳定的漏斗，地下水向抽水井方向流动，污染物也随着水球方向流向抽水井。污染物溶解在水中，按照一定扩散范围随时间变化而变化扩散方向。

该核电站中的放射性物质包含锶和铯，分析不同时间放射性物质的污染羽扩散范围，结果见表1。

表1污染羽情况

放射性物质	分析时间	污染羽范围/m2
锶	3个月	181
	6个月	204
	9个月	214
	12个月	236
	15个月	258
	18个月	315
	21个月	458
	24个月	582
	5年	967
	8年	1287
	10年	2061
铯	3个月	95
	6个月	174
	9个月	198
	12个月	216
	15个月	237
	18个月	258
	21个月	296
	24个月	304
	5年	1066
	8年	1542
	10年	2684

从表1能够看出，两种放射性物质随着时间的推移，污染羽的范围不断扩大，但是随着时间增加，污染羽扩散速度变慢，增长速度不如污染发生初期快，说明随着地下水流动，对污染物存在一定代谢作用。

分析不同扩散距离与中心浓度在不同时间下时间变化情况，结果见表2。

表2污染物变化情况

放射性物质	分析时间	中心浓度/（mg/L)	最大扩散距离/m
锶	6个月	20.45	4.54
	1年	19.62	8.72
	2年	15.24	16.74
	3年	11.65	32.94
	4年	9.54	64.85
	5年	7.14	86.41
	6年	5.43	127.54
	7年	3.12	159.66
	8年	1.54	195.64
	9年	0.95	204.82
	10年	0.25	227.6
铯	6个月	25.64	5.64
	1年	22.92	10.85
	2年	19.57	14.62
	3年	16.85	28.64
	4年	13.71	56.88
	5年	10.54	67.43
	6年	7.42	72.61
	7年	4.88	83.84
	8年	1.54	169.47
	9年	0.94	194.38
	10年	0.14	224.67

分析表2可知，随着时间推移，污染物中心浓度逐渐降低，最大扩散距离逐渐增大，说明污染物随着水流流动，逐渐被稀释，但是也伴随着水流的流动，流动到更远的地方，对其它区域造成污染。

4结论

本文研究基于地质随机场孔隙分析的地下水污染扩散特征模拟，通过地质随机场孔隙分析自动提取构建地下水污染扩散特征模型所需的参数，在此基础上构建水污染扩散特征模型，利用德劳内三角剖分实现研究区域的空间离散，使用GIS空间分析，提取点和面的空间分布参数，通过模型模拟出污染物扩散特征，通过仿真平台获取核电站的相关数据模拟出发生核泄漏的地下水域区域核污染扩散特征。经分析，随着水流场流动，污染物随之发生扩散，污染中心的污染物浓度逐渐降低，污染物向更远的方向扩散，而且污染初期污染物呈现无规律扩散，随着时间增长，污染物向水流涌动中心流动，由此分析出地下水污染物扩散规律与地下水水流方向一致。

参考文献

[1]周雨泽,董少刚,李铱,等.尾矿库渗滤液对地下水污染的数值模拟研究[J].湖南师范大学自然科学学报,2019,42(05):10-16.

[2]张雁,秦秋寒.基于地质成因的砂岩储层微观孔隙结构分形特征分析[J].湘潭大学自然科学学报,2018,40(03):95-99.

[3]李培熙,杨桂莲,李伟,等.基于含水层非均质性随机特征的地下水脆弱性评价[J].水文,2019,39(01):56-59+6.

[4]李莎,成建梅,宫辉力.基于变渗透系数的地下水开采-地面沉降三维模拟研究[J].水文地质工程地质,2018,45(03):14-21.

[5]彭小玉,周理程,吴文晖,等.水污染指数法在湘江长沙段支流水质评价中的应用分析[J].四川环境,2021,40(02):172-177.

[6]郑富新,尹芝华,杜青青,等.某污染场地地下水硝态氮迁移过程的数值模拟[J].环境工程,2018,36(12):103-107.

[7]陈言菲,李翠梅,齐国远,等.基于GMS的江南某地区浅层地下水溶质迁移规律分析[J].水电能源科学,2018,36(08):33-38.

[8]王平,韩占涛,张海领,等.某氨氮污染地下水体抽出-处理系统优化模拟研究[J].水文地质工程地质,2020,47(03):34-43.

[9]卢玉秋,高何凤,狄瑜,等.广西地下水污染现状与防治对策研究[J].四川环境,2021,40(01):100-103.

[10]骆成杰,刘建,刘丹,宋凯,等.某红层地区临河非正规生活垃圾填埋场地下水污染与修复模拟[J].环境工程,2018,36(10):24-29.

[11]潘紫东,卢文喜,范越,等.基于模拟-优化方法的地下水污染源溯源辨识[J].中国环境科学,2020,40(04):1698-1705.

[12]张小茅,周俊,熊小锋,等.地下水环境影响评价中污染物运移模拟软件的适宜性评估[J].环境科学研究,2019,32(01):10-16.

[13]杨悦锁,朱一丹,张文卿,等.地下水系统中镍污染和天然胶体共迁移特征[J].吉林大学学报(地球科学版),2020,50(01):226-233.

[14]李泽岩,黄福杨,刘丹丹,等.海河流域滹沱河冲洪积扇地下水中农药污染及分布特征[J].岩矿测试,2019,38(02):186-194.

[15]高月香,沈欢,张毅敏,等.基于FEFLOW的高尔夫球场地下水污染风险预测研究与效果评估[J].水利水电技术,2018,49(11):144-150.