罗布泊盐湖罗北凹地矿区卤水数值模拟研究

王凌芬，于咏梅，李博昀，张凡凯，邓宇飞，李文学，钦贺

(1.中化地质矿山总局地质研究院， 北京 100101； 2.国投新疆罗布泊钾盐有限责任公司， 新疆 哈密 839000)

0 引言

罗布泊位于塔里木盆地东部，罗北凹地是从罗布泊分割出来的一个次级盆地(刘成林等，2002)，行政区划隶属新疆维吾尔自治区巴音郭勒蒙古自治州若羌县，地理坐标：东经 90°40′00″～91°22′00″，北纬 40°33′00″～41°05′00″，南北长约60 km，东西宽约25 km，面积为1534 km2(李文学等，2018)，是一个赋存多层含钾卤水的大型液体钾盐矿床(陈永志等，2001；王弭力等，2006；李浩等，2008；刘成林等，2009)。目前，罗北凹地液体钾盐矿以开采潜卤水为主，潜卤水矿层平均厚度为17.54 m，矿层平均埋深1.89 m，矿层的平均孔隙度为18.58%，给水度11.36%，卤水平均比重1.22，氯化钾平均品位1.45%，氯化钾孔隙度储量8268.03万t、给水度储量3894.97万t，占罗北凹地内储量的46.5%(李文学等，2018；刘成林等，2020)。随着公司生产规模的扩大，采卤范围逐渐扩大到深部承压卤水，卤水资源不断消耗，地下卤水水位也随之下降(刘传福，2010；韩积斌等，2015)，由于液体钾盐矿与固体钾盐矿的不同，卤水的赋存特征和矿体的厚度及组分也随之发生较大的变化(王新民等，2003；李文学等，2017；王凯等，2020)。实践证明，含卤层为较为复杂的层组，地下卤水的非均质性及富水性的不均一性较为明显(焦鹏程等，2003)。近十年来，前人在矿区开展了一系列勘查、研究工作。焦鹏程等(2003，2006)采用深部承压含水层卤水抽水试验、14C测年技术、同位素和化学示踪技术，研究晶间卤水循环速率、卤水流向和流速及钾富集机理等卤水演化过程；李文学等(2021)利用大气联通法及辅助孔注入卤水法提高采卤井涌水量；作者所在团队于2019—2021年先后在研究区开展疏干开采抽水试验、承压水监测网建设、矿区资源储量核实工作、储卤层钾资源勘探和中深部卤水钾资源调查以及承压卤水开采方法研究等工作，积累了丰富的数据资料和研究基础，但目前对地下卤水资源的管理、规划和开采仍处于较初级的人工计算阶段，存在耗时长、计算量大、效率低、精度不够等问题，本文开展了矿区地下卤水数值模拟研究，通过数值模型的建立，实现对区内钻孔、采卤井及地质剖面、水文地质参数等原始数据的有效利用，细化对罗北凹地矿区地层结构、各含水层空间分布及相互间水力联系的认识，进行地下卤水流场的模拟预测，指导下一步采卤井的布置，实现地下卤水的合理化开采。

1 矿区地质和水文地质

矿区地表出露和钻孔已揭露地层均为第四系化学盐类沉积和碎屑沉积，并呈互层状循环交替产出，反映了本区咸化和淡化沉积环境的规律性变化，其中盐类沉积层是区内地下富钾卤水的主要赋存层位(顾新鲁等，2003，2009；袁文虎等，2020)。以往资料将罗北凹地含盐系划分为7个(编号为S1～S7)含盐组(图1，表1)，含盐组分别对应相应的盐层和碎屑层，富钾卤水主要赋存于各盐层中，形成相应的储卤层，其中罗北凹地浅部(一般小于90 m)划分出W1～W4四个储卤层(1个潜卤水层、3个承压卤水层)。

2 模型概化及建立

2.1 地下卤水渗流模型的建立

研究区地下水流整体上以水平运动为主、垂向运动为辅，地下水系统符合质量守恒定律和能量守恒定律；在常温常压下地下水运动符合达西定律；考虑含水层之间的流量交换，地下水运动可以概化为空间三维流；地下水系统的垂向运动是由层间水头差异引起的；地下水系统的输入输出随时间、空间变化，故地下水为非稳定流；参数随空间变化，体现了系统的非均质性，所以含水介质概化为非均质各向同性介质。综上所述，研究区可概化成非均质、水平方向各向同性、垂向存在变异、空间三维结构、非稳定地下水流系统，即地下水系统的概念模型。

图1 罗北凹地含盐系剖面图(据新疆地矿局第二水文工程地质大队，2006①修改)

表1 罗北凹地含盐系划分表

(1)

式(1)中：Ω—渗流区域；h—地下水系统的水位标高(m)；K—含水介质的水平渗透系数(m/d)；Kz—含水介质垂向渗透系数(m/d)；ε—含水层的源汇项(1/d)；h0—初始水位(m)；S—自由面以下含水层储水率(1/m)；Γ0—渗流区域的上边界，即地下水的自由表面；μ—潜水含水层重力给水度；p—潜水面的蒸发和降水入渗强度等(m/d)；Γ1—已知水位边界；h1—已知边界水位值(m)；Γ2—渗流区域的流量边界；Kn—边界面法线方向的渗透系数(m/d)，n—边界面的法线方向；q—Γ2边界的单位面积上的流量(m/d)，流入为正，流出为负，隔水边界为0。

采用基于有限差分法的GMS软件对上述数学方程进行求解，地下水模拟系统(Groundwater Modeling System)，简称GMS，是由美国Brigham young University的环境模型研究实验室在综合Modflow、Modpath等已有地下水模型基础上研发而成的，是一个具有综合性、用于地下水模拟的图形界面软件。基于GMS软件的Solid模块建立的水文地质结构模型(图2、图3)，为分析研究区的沉积特征提供了方法和思路，为建立地下水数值模拟模型奠定了基础。

整体上说，模型在垂向上分为7层，包含4层含水层和3层隔水层。将模拟区剖分为220行、200列规则网格，各层均采用300 m×300 m网格大小进行剖分，共剖分网格数目为118201个，其中，第一层有效单元(活动网格)17077个，第二至第七层的活动单元格数目均为16854个。

图2 罗北凹地水文地质结构横向剖面示意图W1—潜卤水层；W2—第一承压卤水层；W3—第二承压卤水层；W4—第三承压卤水层

图3 罗北凹地水文地质结构纵向剖面示意图W1—潜卤水层；W2—第一承压卤水层；W3—第二承压卤水层；W4—第三承压卤水层

2.2 边界条件及水文地质参数的确定

图4 研究区边界条件概化图

根据区内流场特征和对地层结构的分析，将研究区侧向边界确定为三种类型(图4)。其一是流量边界，为北部、南部和西部台地与凹地之间的自然边界。这类边界接受山区融雪和降雨入渗到出山口的侧向补给、西台地的侧向补给以及盐田渗漏的侧向补给。补给量的大小按照各块段水文地质条件和渗透系数大小有所不同，具体在计算过程中进行计算分析；其二是定水头边界，将模型的南部边界根据等水位线分布图处理成定水头边界，从等水位线图上看，南部边界处为地下水位的高值区，呈现出类似地下水分水岭的特征，将其处理为一类；其三是通用水头边界，即第三类边界条件，罗北凹地与腾龙台地的边界处，受附近采卤井的影响，该边界定义为通用水头边界。研究区上边界为潜水含水层自由水面，模型底界以钻孔揭穿最大深度处W4的底板作为模型的底部边界，处理为隔水边界。

在建模工作中，首先根据水文地质条件和野外工作获得的水文地质参数，按参数分区给定参数初值，通过水位拟合进行参数识别，最后确定各参数分区值，此处以潜水含水层渗透系数、孔隙度、给水度、第一承压含水层弹性释水系数为例(图5～图8)，其他各层不一一罗列。

图5 研究区潜水含水层渗透系数分区图

图6 研究区潜水含水层给水度分区图

图7 研究区潜水含水层孔隙度分区图

图8 研究区第一承压含水层释水系数分区图

2.3 模拟期确定及定解条件

根据现有各均衡要素资料和水位动态资料情况，确定2019年11月—2020年11月这一个完整的水文年为模拟期，以月为单位共划分为12个应力期。

初始条件：以2019年11月地下水水位资料，采用内插法和外推法获得潜水含水层和第一承压含水层的初始水位。根据现有的第二和第三承压含水层水位观测数据，各承压含水层间水力联系密切，故第二承压含水层和第三承压含水层初始流场采用第一承压含水层初始流场。

2.4 模型识别与检验

对研究区内长观孔的实测水位和模拟水位拟合误差进行了分析。其中潜水含水层拟合的28个钻孔中，平均误差绝对值小于0.5 m的钻孔有20个，平均误差绝对值在0.5～1 m之间的钻孔有7个，两者占总数的96.43%；平均误差绝对值大于 1 m 的钻孔有1个，占总数的3.57%。承压含水层长观孔点位较少，仅有11个，平均误差绝对值小于0.5 m的钻孔有8个，平均误差绝对值在0.5～1 m 之间的钻孔有3个，两者占总数的100%。拟合误差偏大的原因主要有以下三点：(1)误差较大的地区集中在开采降落漏斗附近，资料误差比较大；(2)部分控制点受生产用井影响，观测水位变动过大，造成误差偏大；(3)实测数据使用的是月平均值，计算导致误差。

图9 模拟期末潜水含水层流场拟合图

潜水含水层和承压含水层的流场拟合见图9和图10，典型长观孔水位过程线拟合情况见图11和图12。计算流场与实际流场的拟合程度较好，计算流场基本上反映了地下水的流动与趋势，含水层识别模型所对应的水文地质模型与实际水文地质模型是基本吻合的，边界概化较为合理，源汇项处理较正确，参数分区较为合理，参数基本正确。

图10 模拟期末承压含水层流场拟合图

图11 潜水观测孔水位拟合曲线图

图12 承压水观测孔水位拟合曲线图

3 卤水开采预测

利用所建立的地下水数值模型，结合现有开采情况，预测了地下水流场变化趋势，对地下水开采降深进行了评价。本次预测将2020年11月的统测水位确定的流场作为预测模型的初始流场，预测年限为3年和5年，即从2020年11月到2025年11月。图13和图14为3年末潜水位和承压水位等值线图；图15和图16分别为5年后潜水位和承压水位等值线图；图17和图18为3年后潜水水位降深和承压水水位降深等值线图；图19和图20为5年后潜水水位降深和承压水水位降深等值线图。

图13 预测2023年11月潜水等水位线分布

图14 预测2023年11月承压水等水位线分布

从图13与图14可见，研究区潜水地下水水流方向总体上和初始流场形态基本一致，而承压含水层在开采井附近流场形态发生较大变化，出现明显的降落漏斗。

经3年开采后，研究区潜水总体降深2.00～4.00 m左右(正为水位下降)，其中个别开采区中心的最大降深可达到6.00 m以上。降深面积分布较为广泛，开采区降深2.00 m以上区域面积达376.77 km2。而开采5年后降深2.00m以上区域面积达718.60 km2。

图16 预测2025年11月承压水等水位线分布

图15 预测2025年11月潜水等水位线分布

图17 预测2023年11月潜水水位降深分布图

图18 预测2023年11承压水水位降深分布图

图19 预测2025年11月潜水水位降深分布图

图20 预测2025年11月承压水水位降深分布图

经3年开采后，研究区承压水总体降深0.00～4.00 m左右，其中个别开采区中心的最大降深可达到10.00 m以上。降深面积分布较为广泛，开采区降深4.00 m以上区域面积达194.11 km2。而开采5年后降深4.00 m以上区域面积达338.00 km2。

4 结论

(1)使用GMS软件建立水文地质结构模型，将模型在垂向上概化为7层，给出了参数分区及参数初值，确定的水文地质概念模型为非均质、水平方向各向同性垂向存在变异、空间三维结构、非稳定地下水流系统。

(2)从模拟期末等值线拟合、典型长观孔水位过程线拟合、水文地质参数以及研究区地下水系统均衡分析来看，研究区地下水流数值模拟模型反映了罗北凹地地下水流动规律和特征，符合研究区实际的水文地质条件，反映了研究区内地下水的运动特点和动态变化趋势，基本达到模型精度要求，模型识别和均衡计算结果与实际水文地质条件相符。

(3)通过模型预测，得到开采3年和5年后的地下水流场与降深图。研究区潜水地下水水流方向总体上和初始流场形态基本一致，而承压含水层在开采井附近流场形态发生较大变化，出现明显的降落漏斗。后期在开采井部署时，可考虑远离降落漏斗分布区，在水位降深幅度较小地区加密采卤井设计，这样才有利于卤水的合理采出，减小生产成本。

注 释

① 新疆地矿局第二水文工程地质大队.2006.新疆若羌县罗北凹地钾盐矿详查报告[R].