CN114357839A - 一种矿区地下水三维数值模型构建方法 - Google Patents

Abstract

本发明公开了一种矿区地下水三维数值模型构建方法，包括步骤S100，水流模型、溶质运移模型选取；S200，评价区地地下水数值模型采用Visual MODFLOW软件进行数值离散；S300，模型识别校正及拟合。本发明能够建立渣场区域的地下水三维数值模型，预测渣场区域的地下水影响程度。

Description

一种矿区地下水三维数值模型构建方法

技术领域

本发明涉及钻井钻探技术领域，具体来说，涉及一种矿区地下水三维数值模型构建方法。

背景技术

由于矿藏的不可移动性，导致矿山在长期开采运行过程中出现了严重的环境污染情况，不仅破坏了区域动植物系统，而且导致区域水系结构出现了严重的毁损。

现有的矿山环境治理方法与生态修复技术通常是基于渣场区域的地下水文勘查状况进行的，而目前的地下水文勘查不够直观，也无法用于进一步的预判，为更进一步的预测渣场区域的地下水影响程度，建立渣场区域的地下水三维数值模型势在必行。

发明内容

针对相关技术中的问题，本发明提出一种矿区地下水三维数值模型构建方法，以克服现有相关技术所存在的技术问题。

为此，本发明采用的具体技术方案如下：

一种矿区地下水三维数值模型构建方法，包括步骤S100，水流模型、溶质运移模型选取；S200，评价区地地下水数值模型采用Visual MODFLOW软件进行数值离散；S300，模型识别校正及拟合。

在可能的一个设计中，通过概化得到的非均值各向异性等效连续介质模型，地下水稳定运动数学模型为

式中：H—地下水水头(m)；Kx，Ky，Kz—各向异性主渗透系数(m/d)；Ss—含水层储水率(1/m)；Γ1—模拟区域第一类边界；Γ2—模拟区域第二类边界；H0(x，y，z)-含水层初始水头(m)；Hr(x，y，z)-第一类边界条件边界水头(m)，q0(x，y，z)-第二类边界单位面积过水断面补给流量(m2/d)；ε-源汇项强度(包括开采强度等)(1/d)；Ω-渗流区域。

在可能的一个设计中，溶质运移控制方程为

式中：R-阻滞系数；ρ_b-介质密度；θ-介质孔隙度；C-地下水中组分质量浓度；

介质骨架吸附的溶质质量浓度；t-时间；D_ij-水动力弥散系数张量；v_i-地下水渗流速度；W-水流的源和汇；C_s-源中组分的质量浓度；λ₁-溶解相一级反应速率；λ₂-吸附相反应速率。

在可能的一个设计中，步骤S100之前，还应进行渗透系数、降雨入渗系数以及弥散系数的确定。

在可能的一个设计中，渗透系数如表所示，

在可能的一个设计中，降雨入渗系数如表所示，

地貌+岩性	灰岩+峰丛	泥灰岩、粉砂岩	灰岩+沟谷+洼地
				入渗系数	0.2	0.1	0.6

。

在可能的一个设计中，根据区域弥散试验确定，纵向弥散系数DL为25.5m2/d，横向弥散系数DT为1.95m2/d。

在可能的一个设计中，步骤S200中采用网格剖分，区域的三维尺度在X方向上长度为为5000m，Y方向上长度为7500m，Z方向的长度为160-500m，模拟区域水平网格X为100个网格单元和Y方向网格为150个网格单元。

与现有技术相比，本发明的有益效果为：

本发明能够建立渣场区域的地下水三维数值模型，预测渣场区域的地下水影响程度，具有广阔的应用前景。

附图说明

为了更清楚地说明本申请实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1示出了本发明一种矿区地下水三维数值模型构建方法实施例的模拟区域剖分图。

具体实施方式

为进一步说明各实施例，本发明提供有附图，这些附图为本发明揭露内容的一部分，其主要用以说明实施例，并可配合说明书的相关描述来解释实施例的运作原理，配合参考这些内容，本领域普通技术人员应能理解其他可能的实施方式以及本发明的优点，图中的组件并未按比例绘制，而类似的组件符号通常用来表示类似的组件。

如图1所示，一种矿区地下水三维数值模型构建方法，包括以下步骤：

模型建立：

1.1概念模型

根据渣场区域的水文地质结构和相应的边界，工程区含水岩组与隔水岩组

相间分布，西侧、北侧以白水河为界、东部以翁项组(S1-2w)粉砂岩、钙质砂岩为界，南部以南侧分水岭及沟谷小溪为界；评价区面积约为18km²。

(1)第一类边界条件(已知水头边界条件)

西侧、北侧以白水河、南侧沟谷小溪为已知边界条件，其水位根据地表水水位确定。

(2)第二类边界条件(已知流量边界条件)

东部以翁项组(S1-2w)粉砂岩、钙质砂岩，南部以南侧分水岭为隔水边界，即垂直于该边界方向流量为0。

底部边界定为第二类边界条件，根据钻孔资料，压水试验孔深度50m时，其渗透率为16Lu-160Lu，渗透系数在10^-3-10^-2数量级，为中等透水-极强透水，而钻孔深度在100m左右的抽水试验显示，其渗透系数在10^-7-10^-8m/s，说明随深度其渗透性明显减弱。因此高程700m的底部边界定为隔水边界。

1.2数学模型：

(1)水流模型

通过概化得到的非均值各向异性等效连续介质模型，地下水稳定运动数学模型为

式中：H-地下水水头(m)；K_x，K_y，K_z-各向异性主渗透系数(m/d)；S_s-含水层储水率(1/m)；

-模拟区域第一类边界；

-模拟区域第二类边界；H₀(x，y，z)-含水层初始水头(m)；H_r(x，y，z)-第一类边界条件边界水头(m)q0(x，y，z)-第二类边界单位面积过水断面补给流量(m2/d)；ε-源汇项强度(包括开采强度等)(1/d)；Ω-渗流区域。

(2)溶质运移模型

溶质运移控制方程为

式中：R-阻滞系数；ρb-介质密度；θ-介质孔隙度；C-地下水中组分质量浓度；

介质骨架吸附的溶质质量浓度；t-时间；D_ij-水动力弥散系数张量；v_i-地下水渗流速度；W-水流的源和汇；C_s-源中组分的质量浓度；λ₁-溶解相一级反应速率；λ₂-吸附相反应速率。

1)初始条件

初始条件是指在初始时刻t＝0时评价区域Ω内各点上的浓度分布

C(x，y，z，t)＝C₀(x，y，z)t＝0，(x，y，)∈Ω

式中：C₀(x，y，z)-评价区内已知浓度分布。

2)边界条件通常是指在评价区域的边界线上溶质浓度或浓度通量的变化情况。通常以第一类边界条件为常见。在边界

处，溶质浓度已知为f(x，y，z，t)，则边界条件称为已知浓度边界或称第一类类边界，可表示为：

C(x，y，z，t)＝f(x，y，z，t)(x，y，z)∈Γ₁

对于边界流速比较大的已知浓度的入渗问题，经常可以表达为这类边界条件。边界

处，已知浓度梯度，称为第二类边界，即：

式中：q是已知函数，ni是方向余弦，当多孔介质的外界为隔水、隔溶质的不透水岩体时，通过边界的流量与溶质通量都为0。此时q＝0。

1.3数值模型

根据项目要求，评价区地地下水数值模型采用Visual MODFLOW软件进行数值离散。

(1)网格剖分

建立了地下水渗流的概念模型和数学模型之后，要对渗流区进行离散化(剖分)。将复杂的渗流问题处理成在剖分单元内简单的规则的渗流问题。无论是用有限元法或是用有限差分法进行数值计算。计算结果的精度和可靠性、收敛性及稳定性在很大程度上取决于单元的剖分方法及单元剖分程度，在离散化时遵循两条基本原则。1)几何相似。要求物理模拟模型从几何形状方面接近真实被模拟体。2)物理相似。要求离散单元的特性从物理性质方面(含水层结构、水流状态)近似于真实结构在这个区域的物理性质。网格剖分对计算的精度，及计算的效率有很重要的影响。区域的三维尺度在X方向上长度为为5000m，Y方向上长度为7500m，Z方向的长度为160-500m。模拟区域水平网格X为100个网格单元和Y方向网格为150个网格单元，由于区域内地层单一，地形高差较大，垂向上未进行分层。因此模型模拟区平面上分为15000个网格，节点80802个。模拟区域剖分图如图1。

(2)边界条件

1)侧向边界：

边界类型为第一和第二类边界。

2)上部边界：

上部边界为自由潜水面，作为潜水面边界，其边界条件由大气降水入渗、蒸发排泄等因素确定，这些参数将在下文描述。

2、模型识别校正

2.1识别方法拟合在地下水流数值模拟中占有很重要的地位，正确理解和进行拟合对于提高数值模拟的质量是至关重要的。地下水流数值模型是是人们对地下水流客观认识的一种表达形式。具体的地下水流是模型的原型。拟合则是要求数学模型与流场原型达到科学的相似，即达到与原型等效。进行拟合的过程就是不断改正错误认识，不断修正数学模型。对于地下水流问题，要使模型和原型拟合很好，必须对原型有足够的认识，即对含水层的结构、边界条件、地下水的形成机制以及运动规律等有全面了解，这样才容易拟合。模型识别是一个参数调节的过程，大多采用间接法反求水文地质参数，用尽可能多且分布比较均匀的观测孔进行水位拟合。一般在资料充足的条件下：(1)根据评价区内尽可能多的观测孔进行拟合，作为对模型的识别；(2)绘出模型运行过程中某些时刻的实测水头等值线图，将其与模拟得出的相应时刻的水头等值线图进行对比，即流场对比；(3)对区域内水均衡进行对比，从而确定模型参数。本区域的面积很大，钻孔水位资料缺乏，因此无法绘制区域的地下水等水位线图。因此，本区域模型的拟合将结合现场勘查所得的水位资料进行。

2.2识别资料

利用掌握的钻孔水位定为校正水位。通过数值模型计算出水头分布值并与观测水头值进行对比。如果计算水位与实际水位相差较大，则需先分析出现误差的原因，然后根据参数的变化范围进行调整，也可调整部分补排项(实测资料也可能存在误差)，直到水位拟合较好为止。从而最终确定地下水系统的空间参数分布及补排关系。

2.3拟合

拟合程度的优劣是检验模拟模型能否充分反映水文地质实体模型的重要依据，它是确定地下水数值模拟模型的关键。由于地下水水位的时空变化规律能充分反映地下水系统的结构及各种输入信息(排水量、降水量等)的特征，所以数值法是以水位观测数据为准，以使计算水位与其误差最小为目的，对模型及参数进行调整。通过适当地调整模型参数，使拟合孔实际水位与计算水位相吻合。判断拟合效果好坏，要计算水位与观测水位特征是否一致；拟合后反求出的模型参数是否与实际的水文地质条件匹配。

2.4地下水质分布结合其他赤泥堆场特点，渗漏液中钠的浓度约为2000mg/L。现状未堆存状态下，给定浓度2000mg/l，天然流场条件下污染物(钠)的分布迁移情况，由模拟结果可知，污染物的渗漏会对周围地下水体造成污染，渣场渗漏时污染物的排泄通道为XXX河。

根据项目设计要求，研究区的三维水文地质结构模型和三维地下水渗流模型采用Visual MODFLOW 10.2软件进行建模。Visual MODFLOW软件模块多，功能全，可以用来模拟与地下水相关的所有水流和溶质运移问题。它集成了MODFLOW、MT3D、MODPATH等程序包。同其它类软件相比，GMS软件除模块更多之外，各模块的功能也更趋完善。现将本研究需要用到的模块简介如下。(1)MODFLOW模块程序结构的模块化：MODFLOW包括一个主程序和若干个相对独立的子程序包(Package)。每个子程序中有数个模块，每个模块用以完成数值模拟的一部分。例如河流子程序包用来模拟河流与含水层之间水力联系；井流子程序包用来模拟抽水井和注水井对含水层的影响。离散方法的简单化：MODFLOW采用有限差分法对地下水流进行数值模拟。差分法易于程序的普及和数据文件的规范。其主要缺点是当对某些单元网格加密时，会增加许多额外不必要的计算单元，延长程序的运行时间，随着计算机速度的迅速提高，计算机受网格数量的限制越来越小，差分法的优势越来越大，MODFLOW解决地下水流运动问题已经将含水层剖分到多达360×360×18个网格单元。MODFLOW引进了应力期(StressPeriod)概念：它将整个模拟时间分为若干个应力期，每个应力期又可再分为若干个时间段。在同一应力期，各时间段既可以按等步长，也可以按一个规定的几何序列逐渐增长。而在每个应力期内，所有的外部源汇项的强度应保持不变。这样就简化、规范了数据文件的输入，而且使得物理概念更为明确。求解方法的多样化：迄今为止，MODFLOW已经含有强隐式法、逐次超松弛迭代法、预调共轭梯度法等子程序包。MODFLOW的求解子程序包多样化，应用范围也更为广泛。大量实际工作表明，只要恰当使用，MODFLOW也可以用来解决裂隙介质中的地下水流动问题。(2)MT3D模块，MT3D可以用来模拟可溶性污染物在地下水中的对流、弥散、扩散作用和一些基本的化学反应过程，能够有效处理各种边界条件和外部源汇项。模型中的化学反应主要是一些比较简单的单组分反应，包括平衡或非平衡状态的线性或非线性吸附作用、一阶不可逆反应(如生物降解等)和可逆的动态反应等。MT3D能够适用于各种水文地质条件，包括：①承压、无压或承压-无压含水层；②倾斜或厚度变化的含水层；③指定浓度或通量边界条件；④蒸发、抽水井、河流等多种外部源汇项。MT3D采用块中心的有限差分网格来对评价区进行离散，其节点位于每个网格的中心。在计算时，首先采用MODFLOW模型确定评价区各节点上的水头，然后将水头作为MT3D的输入数据，通过求解对流-弥散方程即可得到污染物浓度在这些节点上随时间的变化规律。(3)MODPATH模块MODPATH是确定给定时间内地下水稳定或非稳定流中质点运移路径的三维示踪模型。它和MODFLOW一起使用，根据MODFLOW计算出来的流场,MODPATH可以追踪一系列虚拟的粒子来模拟从用户指定地点溢出污染物的运动。这种追溯跟踪方法可以用来描述给定时间内井的截获区。(4)PEST模块是用于自动调参的两个模块。在自动进行参数估计时，交替运用PEST来调整选定的参数，并且重复用于MODFLOW等的计算，直到计算结果和野外观测值相吻合。

参数选择参数确定是建立数值模型的关键。一般来说，模型参数与实际参数有一定的区别，这种区别表现在：(1)由于模型结构与实际地下水系统结构不一致，只能在一定精度下建立等效模型，因此，用模型识别方法得到的参数，只能是等效参数；(2)用抽水试验、压水试验等方法确定的参数具有局部性，是不能够满足区域范围的参数要求。为此，确定模型参数时，以岩性分布为基础，以抽水试验和压水试验参数为约束，以水位动态观测资料为拟合要素，应用数值模拟模型识别方法来确定模型参数。

渗透系数岩体渗透性参数目前主要通过压水试验、注水试验等现场原位获得，此外还可通过取样进行室内试验的方式获得。前者尺度大，对于工程有直接应用的价值，后者尺度较小，可作为参考。本项目结合钻孔压水试验，并对比了同类型的工程实例，进行了渗透系数的选取。考虑岩体各向异性，其渗透系数按照三维取值，各岩性的渗透系数见表1。

降水对地下水的补给，受到多种因素影响，影响结果综合反映在地下水位的变化上，影响降水入渗的主要因素包括：①岩石和土壤的物理化学结构的影响土壤颗粒级配决定了土壤的渗透性，颗粒越粗，磨圆及分选性越好，渗透性越强，湿润峰面下移越快，入渗补给量就大，反之入渗补给量小。对于岩溶区域，发育的裂隙或岩溶越多，风化强度越强则入渗补给量就越大。②地面坡降的影响平缓的地面地表径流微弱，可形成良好的补给区，地面坡降越大，易形成地面径流，入渗补给受到一定的影响。在岩溶洼地，降水则全部入渗形成地下径流。③地下水埋深的影响降水的时候地下水位埋深较浅易形成地面径流，使入渗量减小，而埋深增大，下渗雨量大部分充填于非饱和带土壤，使补给量逐渐减小。对于岩溶地区，地下水多以岩溶管道形式向周边排泄，因此溶蚀洼地中降水全部补给地下水；当排泄量小于补给量时，地表水存于溶蚀洼地中，此时补给量会受到蒸发的影响。④降水量及降水形式的影响当地下水埋深大于最佳埋深时，入渗补给量随雨量的增大而增加，当埋深小于最佳埋深时，受地下水库容制约，易形成蓄满产流，入渗补给量与降水量成反比关系。各种降水形式中以雨对地下水的补给量最大，其它降水对地下水补给价值不大。在各种降雨类型中，霪雨的雨率不大，但雨时长，雨面广，对地下水的补给具有很大意义；细雨的雨量、雨率皆不大，易被蒸发所消耗，暴雨的雨率大，但雨时短，雨水大部分来不及渗入地下水，多呈地表径流流失。⑤雨强的影响降水强度小于入渗速度时，雨强增大，入渗量增加，反之，入渗量减小。当雨强与入渗速度相近时，补给地下水为最佳。⑥植被的影响植被可缓减地表径流，延长入渗补给时间，增加入渗补给量，但同时，植被使土壤水分蒸散发加快，相对减少了入渗补给。参考其他文献和区域的地形和岩性特征，区域的降入入渗分区及参数取值可见表2。

地貌+岩性	灰岩+峰丛	泥灰岩、粉砂岩	灰岩+沟谷+洼地
				入渗系数	0.2	0.1	0.6

水动力弥散系数是表征在一定流速下，多孔介质对某种污染物质弥散能力的参数，它在宏观上反映了多孔介质中地下水流动过程和空隙结构特征对溶质运移过程的影响，单位常用m²/s。水动力弥散系数是一个与流速及多孔介质有关的张量，即使几何上均质，且有均匀的水力传导系数的多孔介质，就弥散而论，仍然是有方向性的，即使在各向同性介质中，沿水流方向的纵向弥散和与水流方向垂直的横向弥散不同。一般地说，水动力弥散系数包括机械弥散系数与分子扩散系数。当地下水流速较大以致于可以忽略分子扩散系数，同时假设弥散系数与孔隙平均流速呈线性关系，这样可先求出弥散系数再除以孔隙平均流速便可获取弥散度。弥散度只与介质本身特征有关。根据区域弥散试验确定，纵向弥散系数DL为25.5m²/d，横向弥散系数DT为1.95m²/d。

尽管为使解释简单化将上述方法图示并描述为一系列动作，但是应理解并领会，这些方法不受动作的次序所限，因为根据一个或多个实施例，一些动作可按不同次序发生和/或与来自本文中图示和描述或本文中未图示和描述但本领域技术人员可以理解的其他动作并发地发生。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种矿区地下水三维数值模型构建方法，其特征在于，包括步骤S100，水流模型、溶质运移模型选取；S200，评价区地地下水数值模型采用Visual MODFLOW软件进行数值离散；S300，模型识别校正及拟合。

2.根据权利要求1所述的一种矿区地下水三维数值模型构建方法，其特征在于，通过概化得到的非均值各向异性等效连续介质模型，地下水稳定运动数学模型为

式中：H-地下水水头(m)；Kx，Ky，Kz-各向异性主渗透系数(m/d)；Ss-含水层储水率(1/m)；Γ1-模拟区域第一类边界；Γ2-模拟区域第二类边界；H0(x，y，z)-含水层初始水头(m)；Hr(x，y，z)-第一类边界条件边界水头(m)，q0(x，y，z)-第二类边界单位面积过水断面补给流量(m2/d)；ε-源汇项强度(包括开采强度等)(1/d)；Ω-渗流区域。

3.根据权利要求2所述的一种矿区地下水三维数值模型构建方法，其特征在于，溶质运移控制方程为

式中：R-阻滞系数；ρ_b-介质密度；θ-介质孔隙度；C-地下水中组分质量浓度；

-介质骨架吸附的溶质质量浓度；t-时间；D_ij-水动力弥散系数张量；v_i-地下水渗流速度；W-水流的源和汇；C_s-源中组分的质量浓度；λ₁-溶解相一级反应速率；λ₂-吸附相反应速率。

4.根据权利要求1-3任一项所述的一种矿区地下水三维数值模型构建方法，其特征在于，步骤S100之前，还应进行渗透系数、降雨入渗系数以及弥散系数的确定。

5.根据权利要求4所述的一种矿区地下水三维数值模型构建方法，其特征在于，渗透系数如表所示，

6.根据权利要求4所述的一种矿区地下水三维数值模型构建方法，其特征在于，降雨入渗系数如表所示，

地貌+岩性 灰岩+峰丛 泥灰岩、粉砂岩 灰岩+沟谷+洼地 入渗系数 0.2 0.1 0.6

。

7.根据权利要求4所述的一种矿区地下水三维数值模型构建方法，其特征在于，根据区域弥散试验确定，纵向弥散系数DL为25.5m²/d，横向弥散系数DT为1.95m²/d。

8.根据权利要求7所述的一种矿区地下水三维数值模型构建方法，其特征在于，步骤S200中采用网格剖分，区域的三维尺度在X方向上长度为为5000m，Y方向上长度为7500m，Z方向的长度为160-500m，模拟区域水平网格X为100个网格单元和Y方向网格为150个网格单元。

Images