人工溶采作用下非均质盐湖储卤层渗透系数动态演化机制

Dynamic Evolution Mechanism of Hydraulic Conductivity in Heterogeneous Salt Lake Brine Aquifers during Artificial Solution Mining

摘要

人工溶采技术将储卤层蒸发盐矿物转化为卤水，对盐湖资源开发具有重要意义.然而补水溶矿引发含水层渗透性演化机制尚未得到充分阐述.研究基于Python开发耦合MODFLOW6和PhreeqcRM的数值模拟工具MF6PQC，系统研究卤水反应运移对储卤层渗透系数及溶采过程的影响.结果表明，含水层初始非均质结构决定了渗透系数的时空演变规律.溶矿初期在高渗区优先发生地球化学反应，光卤石等高活性矿物溶解使孔隙度与渗透系数显著增大，并在对流-弥散的正反馈作用下形成优势渗流通道.均质或高渗连通结构有利于溶浸剂均匀波及，而强连通或含低渗隔挡的弱连通地层则使固体矿物无法被有效接触，从而降低整体溶采效果.研究深化对储卤层渗透性变化的认识，为优化卤水溶采提供理论依据.

Abstract

Artificial solution mining technology, which converts evaporite minerals in brine aquifers into brine, is crucial for the sustainable development of salt lake resources. However, the dynamic evolution of aquifer hydraulic conductivity induced by mineral dissolution during water injection remains insufficiently understood, hindering accurate process prediction. In this study, a Python-based modeling tool, MF6PQC, coupling MODFLOW6 and PhreeqcRM, was developed to systematically investigate the effects of reactive transport on the hydraulic conductivity of brine aquifers and the overall solution mining process. Simulation results show that aquifer heterogeneity governs the spatiotemporal evolution of hydraulic conductivity. During the early stage of dissolution mining, hydrogeochemical reactions preferentially occur in high permeability zones. The dissolution of highly reactive minerals such as carnallite significantly enhances porosity and hydraulic conductivity, ultimately forming preferential flow paths driven by positive advection-dispersion feedback. Relatively homogeneous aquifers or those with extensive, well-connected high-permeability zones facilitate uniform lixiviant distribution and achieve higher solid-to-liquid conversion efficiency. In contrast, strongly preferential or poorly connected formations interrupted by low permeability barriers limit mineral contact and dissolution, thereby reducing overall solution mining efficiency. This study deepens the understanding of hydraulic conductivity evolution in brine aquifers during water injection and provides a theoretical basis for optimizing salt lake brine resource exploitation.

Graphical abstract

关键词

盐湖地下卤水 / 反应运移数值模拟 / 非均质随机场 / 储卤层渗透系数 / 优势渗流通道 / 环境科学.

Key words

salt-lake brine / reactive transport modeling / spatial heterogeneity / aquifer hydraulic conductivity / preferential flow / environmental science

0 引言

盐湖卤水中富含多种矿产资源，是生产氯化钾和硫酸钾等化工原料的主要来源（Li et al.，2022）.近年来，以柴达木盆地为代表的盐湖区因过度开发导致地下水位持续下降，卤水资源枯竭和贫化等现象突出.同时，不断扩张的水位降落漏斗使主采区潜水面以上赋存大量固体钾矿难以被利用（柴达木综合地质矿产勘查院，2023）.这些广泛分布的钾盐矿物虽品位较低（约<3%），但其资源总量（>2.94亿 t）仍然相当可观.

为保障钾资源持续开采，盐湖矿区广泛采用了“固转液”技术，通过回灌老卤或不饱和卤水以溶解目标固体矿物（宋高等，2024）.尽管一系列溶浸实验证明了人工溶采技术的可行性（王文祥，2013），然而大规模现场应用效果往往低于理论预期，溶矿效率不稳定且过程难以预测（李梦玲，2023）.有研究表明，溶浸开采过程是一个复杂的水-岩相互作用系统.注入溶剂与储卤层矿物发生一系列溶解与沉淀反应，不仅改变卤水的化学组分，且重塑了储卤层的多孔介质结构（Vásquez et al.，2013）.其中，固体矿物溶解使孔隙度与渗透系数提高，而次生矿物沉淀则引发孔隙堵塞并降低渗透性（常文静等，2024）.盐湖储层中的蒸发盐矿物（石盐、光卤石、钾石盐等）对卤水化学条件极为敏感，溶解速率快、反应活性高，进一步加剧了这一过程的复杂性.此外，盐湖地层本身具有极强的非均质特征，孔隙结构和水文地质参数在空间上变化剧烈（张岩等，2017；胡舒娅等，2022；赵宪福等，2023）.这种现象与水岩反应诱导的渗透性改变相互叠加，极易在储层中形成优势流通道（Weisbrod et al.，2012；Li et al.，2020），导致溶浸剂“逃逸”，使大面积固体矿产资源无法被有效接触和溶解，限制了疏干区固体矿资源开发利用.

现有研究从不同角度探讨了人工溶采过程中水岩反应与地下卤水的运移行为.聚焦于水文地球化学机制，研究者通过室内实验与热力学平衡/动力学建模，识别了关键盐类矿物的溶解结晶行为（洪显兰等，1994；王晓等，2015）、相变路径（李瑞琴等，2021）及反应速率参数（郝爱兵和李文鹏，2003），为理解盐湖固液转化规律提供了理论支撑.然而此类研究往往将流体流动过程简化处理或视为静态背景，未能充分计算实际渗流场对化学反应发生位置及空间分布格局的动态控制作用.另一方面，水文地质领域研究则更侧重于盐湖区卤水渗流及溶质运移模拟，采用有限差分、有限元等数值方法或相关模拟工具（黄鸿蓝等，2024），在区域或矿区尺度上求解不同开采情境下的水动力场与溶质浓度场的时空变化特征（周训等，2006；李辉，2010）.但这些模型大多忽略矿物溶解-沉淀反应过程，或将含水层渗透性视为固定不变的静态参数，未能反映反应活动对孔隙结构及渗流条件的反馈调控作用（Steefel et al.，2015）.此外，高盐度卤水体系中普遍存在的变密度现象常被忽视，导致模拟结果与实际物理过程产生偏差（任倩慧，2017）.

由于盐湖区补水溶矿同时涉及多组分溶质运移、变密度流动、水-岩作用及含水层渗透性的动态演化，其本质是一个高度耦合的多物理场过程.溶剂迁移路径和速度决定了其与孔隙介质的接触效率，而矿物溶解与沉淀又会改变孔隙度和渗透率，反过来影响溶剂的运移行为（图1）.因此，构建能够耦合上述多场的反应运移模型，是精确描述盐湖地下卤水动态行为的必要条件.然而，目前针对多孔介质盐湖地下卤水的专用仿真技术尚未成熟.现有反应运移软件，如PHT3D、Crunchflow和TOUGHREACT等在盐湖系统应用中面临诸多局限，或不支持计算高盐溶液的离子活度，或无法自由定义储层渗透性的动态演化（Xie et al.，2015）.因此，亟需发展一套灵活耦合水岩反应与渗流演化的计算工具，为刻画储卤层水文地质参数演化提供技术支撑.

针对上述问题，本研究基于Python环境结合MODFLOW6和PhreeqcRM平台，开发了一套基于顺序非迭代方法的盐湖卤水反应运移工具MF6PQC.MF6PQC能够显式耦合化学反应与溶质运移过程，同时根据水岩作用实时更新孔隙度和渗透系数.通过构建典型盐湖的二维剖面溶矿模型，系统模拟不同类型的非均质条件对储卤层渗透系数演化路径与特征的影响，并分析了渗透系数变化与矿物相变、流体密度变化的空间耦合关系.研究旨在揭示非均质储卤层渗透系数的动态响应机制，为盐湖卤水资源的优化溶采和环境保护提供理论支持.

1 建模原理与方法

1.1　地下卤水渗流及溶质运移方程

在没有人为干扰的条件下，盐湖地下卤水几乎处于静止状态.人工溶采活动破坏了这一平衡，使卤水在孔隙和晶间介质中产生流动.在渗流模拟中，通常采用连续介质假设，将含水层介质的复杂微观结构抽象为宏观体系，并以孔隙度、渗透系数等水文地质参数加以表征.考虑到卤水具有高浓度和变密度的流体特性，需在渗流控制方程中引入密度项耦合（Langevin et al.，2022）：

\nabla \cdot \left[K\cdot \left(\nabla h+\frac{\rho -{\rho }_f}{{\rho }_f}\nabla z\right)\right]+W^{\text{'}}=S_s\frac{\partial h}{\partial t},(1)

\rho ={\rho }_f+{\displaystyle \sum _{i=\mathrm{1}}^N}{D}_i\left(C_i-C_i^{\text{'}}\right),(2)

（1）~（2）中，K是渗透系数（m/d）；h为等效淡水水头（m）；\frac{\rho -{\rho }_f}{{\rho }_f}\nabla z为浮力项；体积源汇项W'表示单位体积含水层中单位时间的注入或抽出体积比（d^-1）；S_s 为比储水系数（m^-1）.ρ和ρ_f 为卤水密度和参考密度（kg/m³），N为决定流体密度的溶质种类数量，C_i 和C_i'分别为溶质i的质量浓度和参考浓度（kg/m³）；D_i 为描述密度随化学物质浓度变化的斜率.在本研究中，密度由PhreeqcRM模块计算，因此D_i 为1.0，以确保模型接收热力学计算得到的真实密度.

溶质运移过程则通过多组分对流-弥散方程描述.对于任意组分i，其浓度变化满足以下控制方程（Patrick and Franklin，1998）：

\varphi \frac{\partial C_i}{\partial t}+\nabla \cdot J_i+u\cdot \nabla c_i=R_i+S_i,(3)

式（3）中，C_i 为组分i的摩尔浓度（mol/m³）；J_i 为组分i的扩散通量（m²/d），包含机械弥散和分子扩散；u·∇c_i ​为对流项；R_i 表示组分的i水-岩作用（mol/（m³·d）），包括溶解、沉淀和离子吸附等反应；S_i 为外部源汇项（mol/（m³·d））.本研究中，地下卤水渗流和溶质运移均采用FloPy（Hughes et al.，2024）调用MODFLOW6（U.S. Geological Survey，2017）数值求解.

1.2　卤水-矿物地球化学反应模型

干盐滩以其蒸发盐矿物为主导的地质特性，影响了矿区水文地球化学过程的动态特征.已有研究通过X射线粉晶衍射及扫描电子显微镜等技术手段，明确盐湖区主要蒸发盐矿物包括石盐、石膏、钾石盐、光卤石及杂卤石；粉砂碎屑层中的主要矿物则为石英、伊利石和钠长石（柴达木综合地质矿产勘查院，2023）.为构建一个代表性的地球化学系统，本研究选取5种主要蒸发盐矿物组合，其物化参数见表1.

蒸发盐矿物表现出高溶解速率，能迅速达到溶解-析出平衡（郝爱兵和李文鹏，2003）.因此，本文采用热力学平衡假设，将主要矿物组分的反应速率视为单位时间步长内完成.为保证高盐度卤水体系的模拟精度，采用PhreeqcRM进行水文地球化学模拟计算（Parkhurst and Wissmeier，2015），并选择Pitzer活度模型作为溶液化学计算数据库.Pitzer数据库尤其适用于高离子强度条件下的水化学系统，能够准确描述离子间相互作用（Lu et al.，2022）.表2总结了本文主要涉及的反应式及其平衡常数.

1.3　储卤层孔隙度-渗透系数关系模型

地下卤水与蒸发盐矿物间存在强烈水岩作用，易使多孔介质物理结构发生显著变化.这种动态过程可视为孔隙度ϕ与时间t的函数，定义为孔隙体积 V_p （t）与总有效体积V的比值（Kaufmann，2016）：

\varphi \left(t\right)=\frac{V_p\left(t\right)}{V}.(4)

然而，准确计算任意时刻t的孔隙度需全面考虑矿物反应和溶质运移机制，直接采用上式难以实现.在实际研究中，孔隙度变化常被简化为与矿物含量直接关联的形式，即t+1时刻孔隙度ϕ_t+1取自ϕ_t 和m种矿物在单位时间内的变化情况（Sun et al.，2023）：

{\varphi }_{t+\mathrm{1}}={\varphi }_t+{\varphi }_t\cdot {\displaystyle \sum _i^m}{V}_{molar,i}\cdot q_i\cdot M\left(q_i\right),(5)

M\left(q_i\right)=\left\{\begin{array}{c}+\mathrm{1},q_i\le \mathrm{0}\\ -\mathrm{1},q_i>\mathrm{0}\end{array}\right.,(6)

式（5）~（6）中，q_i 为矿物i在单位时间消耗或生成量（mol）；V_molar，i 为矿物i的摩尔体积（m³/mol），具体取值见表1；M（q_i ）为q_i 的阶跃函数，当矿物析出时取负值，溶解时取正值.

孔隙度-渗透率关系通常由Kozeny⁃Carman（K⁃C）方程表示（骆祖江等，2024；Rehman et al.，2024）：

\frac{k_{t+\mathrm{1}}}{k_t}=\frac{{\varphi }_{t+\mathrm{1}}^{\mathrm{3}}(\mathrm{1}-{\varphi }_t{)}^{\mathrm{2}}}{{\varphi }_t^{\mathrm{3}}(\mathrm{1}-{\varphi }_{t+\mathrm{1}}{)}^{\mathrm{2}}},(7)

式（7）中，k为渗透率（m²）.在本研究中，由于采用渗透系数K（m/d）表征含水层过水能力，需将渗透率转换为渗透系数.二者关系为：

K=\frac{k\rho g}{\mu },(8)

式（8）中，ρ为流体密度（kg/m³），g为重力加速度（m/s²），μ为动力黏度（Pa·s）.结合式（7）和（8），可得渗透系数的更新公式为：

\frac{K_{t+\mathrm{1}}}{K_t}=\frac{{\rho }_{t+\mathrm{1}}}{{\rho }_t}\cdot \frac{{\mu }_t}{{\mu }_{t+\mathrm{1}}}\cdot \frac{{\varphi }_{t+\mathrm{1}}^{\mathrm{3}}(\mathrm{1}-{\varphi }_t{)}^{\mathrm{2}}}{{\varphi }_t^{\mathrm{3}}(\mathrm{1}-{\varphi }_{t+\mathrm{1}}{)}^{\mathrm{2}}}.(9)

由式（9）可知，渗透系数变化不仅依赖孔隙度演化，还受流体密度与黏度的影响.鉴于目前缺乏适用于高盐度多组分卤水体系的成熟黏度-浓度公式模型，且研究区溶浸剂与原卤黏度差异有限，本研究将μ项近似视为常量处理，主要考虑孔隙度与流体密度的贡献，以避免引入难以验证的参数不确定性.

综上，式（5）~（9）构建了一个能动态刻画孔隙度演化及其对渗透性能影响的耦合模型框架，为后续数值模拟中同时更新含水层各网格孔隙度与渗透系数提供计算依据.

1.4　顺序非迭代耦合方法

在耦合溶质运移与化学反应的数值模拟中，复杂矿物体系的水‑岩作用引入了强非线性与多尺度耦合，使传统方法难以直接求解反应运移方程组.为解决这一难题，本研究引入算子分裂框架下的顺序非迭代方法（sequential non⁃iterative approach，SNIA）（Ahmadi et al.，2022），基于Python开发了MODFLOW6与PhreeqcRM的耦合工具MF6PQC，以描述多孔介质卤水中物理化学过程的动态演化.MF6PQC的执行流程（图2）如下.

（1）设定研究区水文地质边界条件、初始水头与溶质浓度，并赋予各网格单元初始矿物组成、孔隙度及渗透系数等参数.

（2）在时刻t内，通过MODFLOW6求解变密度地下水流方程，获得更新后的水头场和单元间流量.随后利用这些流动信息求解对流-弥散项，计算因运移引起的溶质浓度变化.该步骤通常采用隐式时间离散格式（Langevin et al.，2022）.

（3）将各单元溶质浓度与孔隙水体积等输送至PhreeqcRM接口，后者调用相应反应数据库，并行计算矿物溶解沉淀、离子交换及表面络合过程，输出离子浓度、矿物体积变动和流体密度ρ_t+1.在此框架下，密度不再依赖式2溶质浓度线性近似，而由PhreeqcRM直接返回反应后的密度.

（4）依据反应导致的矿物相体积变化q_i，运用式（5）~（6）更新孔隙度ϕ_t+1，并基于ϕ_t+1与流体密度ρ_t+1，按式（7）更新各网格渗透系数K_t+1.

（5）重复步骤（2）~（4），直至模拟达到预设的总时长，或关键参数（如特定位置的浓度、总溶解矿物量）的变化小于收敛容差.

在时间离散方面，MF6PQC采用固定步长的算子分裂策略.由于SNIA方法为一阶时间精度，步长选择对数值稳定性和模拟精度至关重要.为此，在模型应用前需对不同耦合步长进行了收敛性与可行性检验，以在精度与计算开销之间取得平衡.

相比传统全耦合方法，SNIA无需在每个时刻反复求解非线性方程组，降低了计算复杂性和收敛风为验证MF6PQC框架的可靠性，本研究在场景模拟前测试了多项对流-弥散与平衡/动力学反应相耦合的基准案例.模拟结果与PHT3D等模型的参考解和文献报道高度一致，证实耦合接口与核心模块的准确性.MF6PQC的开发代码及验证案例均已上传至代码托管平台（见文末说明），供后续研究参考.

2 数值模拟与情景设计

2.1　水文地质概念模型

为探究人工溶采作用下储卤层渗透系数演化过程，本研究构建了一个二维垂直剖面的数值模型.该模型旨在概化典型盐湖近地表环境中，通过补水渠和采卤渠进行固液转化溶采的代表性情景.模型计算域为水平方向400 m、垂直方向60 m的矩形区域.为进行数值求解，采用结构化网格对该区域进行离散化处理，水平方向划分为400个单元，垂直方向划分为100个单元，共包含40 000个网格单元（图3）.模拟总时长为5年（1 825天），时间上离散为180个步长，大约每10天进行一次运移-反应耦合计算.为验证时间离散的合理性，在模型应用前对多种步长方案进行了对比分析，并据此选取10天作为最终设定.

定解条件包括水力边界、溶质运移边界及初始条件.水力边界方面，模型剖面左上侧（0 m≤x≤2 m，55 m≤z≤60 m）被指定为60 m的恒定水头边界，代表补水渠持续注水作用；剖面右上侧（0 m≤x≤2 m，55 m≤z≤60 m）则设定为5 m³/d的固定流量抽采边界，模拟稳定的采卤过程.模型底部、顶部（除渠口影响段外）及两侧垂直边界均设为零水力通量边界，表明除人工干预外系统与外界无其他水力交换.针对溶质运移过程，补水渠内溶浸剂被设定为近饱和NaCl溶液，同时所有零水力通量边界亦为零溶质弥散通量边界.

初始条件方面，初始水头场被均匀设置为60 m.初始孔隙水化学成分参考了察尔汗盐湖达布逊段观测孔的实测水样数据（李梦玲，2023），并经电荷平衡计算及根据矿物组合调整至近饱和状态（表3），以表示相对稳定的地球化学背景.模型中各计算单元初始矿物含量依据表1提供的参考值设定（表4），每个单元格的总矿物含量设定为100 mol，其中阳离子交换作用由单一类型的交换位点X^-驱动.此外，区域内初始孔隙度为0.3，水平与垂向方向弥散度分别为0.005 m和0.000 5 m.

2.2　非均质渗透系数随机场情景设计

为有效模拟盐湖区储卤层固有的渗透系数空间非均质性，本研究采用地统计学方法生成随机对数渗透系数场.随机场生成依赖协方差模型描述空间相关性，其核心由半变异函数γ（r）定义（Müller et al.，2022）：

\gamma \left(h\right)={\sigma }^{\mathrm{2}}\cdot \left(\mathrm{1}-cor\left(s\cdot \frac{h}{l}\right)\right)+n,(10)

式（10）中，h为任意两点间的空间距离；σ²为结构化方差；ℓ为相关长度；s为尺度因子（本研究中设为1）；n为块金效应（本研究中设为0）；cor（·）为归一化的相关函数.考虑到地质介质渗透性变化的复杂性，本文选择Matern相关函数（Chaudhuri et al.，2025）定义，其表达式为：

cor\left(h\right)=\frac{\mathrm{1}}{{\mathrm{2}}^{\nu -\mathrm{1}}\Gamma \left(\nu \right)}{\left(\frac{h}{l}\right)}^{\nu }{K}_{\nu }\left(\frac{h}{l}\right),(11)

式（11）中，Γ（⋅）是伽马函数，K_ν （⋅）是第二类修正贝塞尔函数，ν是控制随机场平滑度的正参数.Matern函数生成的随机场具有ν-1阶可微性，相较于高斯函数（ν→∞）与指数函数（ν=0.5），能够在随机场连续性和地质特征粗糙度之间取得较好的平衡（Pardo-Iguzquiza and Chica⁃Olmo，2008）.本研究中，平滑度参数ν设为1，以模拟具有一定连续性但非过度平滑的渗透场.

在此基础上，共设计了5种具有不同非均质特征的场景（图4）.所有随机场均作用于对数渗透系数lnK.其中，场景A为均质场，作为基准对照；场景B是由式（8）~（9）生成的标准对数高斯场；场景C和D则是对场景B的Zinn⁃Harvey变换（Zinn and Harvey，2003），表示具有优先高渗通道或低渗阻隔特征的连通性结构（Li et al.，2022）；场景E则将场景B转换为反正弦分布.为体现渗透系数的各向异性，各场景水平均设相关长度ℓ_x =15 m，垂直方向ℓ_z =5 m.方差σ²则针对不同场景手动调整，确保所有场景下的K值均控制在0.001~85.000 m/d区间内，以便于横向对比分析.各场景的定义与参数见表5.

3 结果与讨论

3.1　矿物组分与渗透系数时空变化特征

研究表明含水层渗透性普遍存在对数高斯分布的连续性随机变化特征（Wang et al.，2024），因此本节聚焦于场景B的模拟结果.图5展示了采补剖面内溶矿引发的矿物动态演化过程.可以看出，溶浸剂运移优先沿初始高渗透区域运移，并打破其原有的地球化学平衡，诱发一系列水岩反应.在溶浸剂推进前缘，光卤石因高溶解度及不饱和状态（图5b）使其快速溶解，并向孔隙水中释放大量K⁺、Mg²⁺及Cl^-离子，进而产生次生钾石盐沉淀（图5d）；同时这一过程也造成微量（<3%）石盐过饱和析出（图5a）.

随着补水过程持续，在原生光卤石完全溶解区域，次生钾石盐受溶浸剂冲刷转为不饱和并发生再溶解.这一过程使石盐沉淀进一步增加（>4%）.当钾石盐也完全溶解后，石盐少量溶解并趋于稳定.此时溶液化学环境的变化促使低溶解度的杂卤石缓慢溶解，释放的Ca²⁺和SO₄^2-也引起了次生石膏沉淀.

上述溶解、沉淀与再溶解过程随溶剂前缘推进不断上演，并最终塑造了矿物蚀变区域的复杂几何形态.然而，矿物反应空间分布形态严格受初始非均质渗透系数场的控制，在高渗透区域反应活跃，而低渗透区域则几乎未发生明显的地球化学变化.这一模拟结果解释了溶矿区岩芯观测中，部分地区光卤石保存完好而附近区域却溶解严重的现象.

图6为上述矿物相变过程引起的孔隙度和渗透系数演化情况.其中，水文地质参数的显著变化区域与地球化学反应影响范围高度吻合，主要集中在溶剂波及区域，尤其是在靠近补水渠边界并沿初始高渗透路径延伸的地区.在5年的模拟期内，剖面整体平均孔隙度从初始0.3增长至0.33，平均渗透系数也从2.06 m/d提升至3.48 m/d，反映矿物净溶解对含水层导水性能的整体贡献.然而，参数变化在空间上分布极不均匀，局部区域的变化远超平均水平.在强烈局部不均衡的溶蚀作用下，孔隙度最大值可增至0.51，渗透系数最大值由85.50 m/d增长至峰值180.85 m/d，又回落至166.82 m/d.这种非单调变化趋势反映了高渗透区域光卤石溶解带来的渗透性急剧增大，并与次生矿物沉淀对孔隙结构产生复杂叠加的动态平衡.

从空间分布形态上看，孔隙度和渗透系数演化路径强化了初始非均质渗透场的结构特征，变化最剧烈区域恰好对应于初始渗透系数较高的地区.这说明初始渗透性越高的区域，渗透系数增幅也越大.这种正反馈作用使原有高渗透区域被进一步溶蚀、贯通和扩展，最终在含水层中形成了高度连通的优先渗流通道.这些通道主导了后期溶浸剂的运移方向和分布，也决定了化学反应发生的主要场所，使得溶浸剂更容易到达至采卤渠，而绕过了大部分低渗透性区域.

总体而言，孔隙结构与渗透性能的演化是光卤石等易溶矿物溶解造成的孔隙体积增大，并与钾石盐和石盐等次生矿物沉淀造成的孔隙体积减小之间动态竞争的结果.同时，竞争性的水岩反应过程受到含水层初始渗透性空间分布的控制.

3.2　储卤层渗透系数演化影响因素分析

为深入解析储卤层渗透系数动态演化控制机制，本研究结合运用了Pearson相关性分析（图7）与随机森林（RF）模型的特征重要性评估（图8）.其中，前者旨在揭示初始渗透系数（K₀）、流体密度以及各矿物含量变化与渗透系数变化量（ΔK）间的线性相关性，后者量化了模拟第1、3、5年各因素对ΔK空间变异性的相对贡献程度.

RF模型的构建基于单次模拟所生成的4万（400×100）个网格点时空演化数据.每个网格点被视为独立样本，其输入特征包括初始渗透率、矿物含量变化及流体密度变化，对应的输出目标为ΔK.经由RF模型训练后，通过统计各特征在决策树分裂节点中所引起的预测误差下降量，并在全部树模型中进行累积与归一化，得到特征重要性指标.该指标量化了各因素对ΔK预测的贡献权重，从而揭示其在渗透系数空间变异性中的相对影响程度.

在模拟初期，ΔK主要局限于补水渠附近，局部水文地球化学反应是控制渗透性改变的首要因素.此时光卤石含量变化与ΔK呈现显著的负相关性（r <-0.75），表明光卤石作为主要溶解矿物，其含量减少直接对应渗透系数增加.同期次生沉淀的石盐与ΔK表现为显著正相关（r >0.75），是因为石盐骨架有着较大的初始含量基数，使得小幅沉淀产生的总体积变化对孔隙结构的影响依然显著.RF评估结果亦证实该阶段光卤石（λ=0.456）和石盐（λ=0.382）是影响ΔK空间分布的最关键因素.

随着对流-弥散作用促使溶质运移和化学反应前缘不断扩展，水动力对化学反应空间格局的塑造作用开始凸显.K₀与ΔK的正相关性和权重逐渐增强，并在模拟接近第4年时达到峰值（r=0.660，λ=0.586）.与此同时，补水渠近端光卤石消耗殆尽，贡献率相应下降（λ=0.171），溶解度相对较低的杂卤石（λ=0.102）和石膏（λ=0.125）的溶解-沉淀作用在局部调整区域渗透性.石盐的影响降至最低（λ=0.004），表明已形成的优势流动通道内，石盐已接近局部的溶解-沉淀平衡，不再是ΔK的主要因素.这一阶段地球化学反应过程虽仍在进行，但溶质运移路径已基本成型.

此后，系统逐渐趋向动态平衡状态，各因素与ΔK相关性（图7）趋于稳定.K₀成为最主要的正相关因素（r =0.60，λ=0.524），流体密度变化则成为最显著的负相关因素（r =-0.50，λ=0.213）.这表明在长期溶采作用下，高渗透性通道持续受到密度相对较低的溶浸剂冲刷和强烈的溶蚀作用，其渗透性增幅最为显著；而密度较大的原生卤水在流动性较差的区域运移受阻，渗透系数增幅较小甚至可能因沉淀而降低.

由此可知，储卤层渗透系数演化是一个受多因素、多阶段复杂调控的过程.从时间演进角度看，模拟早期渗透性的剧烈变化主要由水文地球化学反应，特别是易溶矿物（如光卤石）的快速溶解所启动和主导；随着溶采过程的持续，优势流动通道的形成、发展与强化，使初始渗透系数场的结构性特征以及流体密度差异所驱动的流动分异逐渐成为主导因素，最终决定了渗透性演化的宏观空间格局和长期演变趋势.

3.3　非均质结构对渗透性空间演化格局及溶矿效率的影响

为进一步评估储卤层初始非均质空间结构对渗透系数演化路径及最终溶矿成效的影响，对比分析了5种非均质场景下ΔK模拟结果（图9）.基于各场景K₀的中位数，将研究区域划分为高渗透区与低渗透区，并统计二者随模拟时间的相对变化率（图10）.

在完全均质分布的场景A中，K₀空间一致性使溶浸剂宏观运移路径受控于水力梯度及流体间密度差异.由于溶浸剂密度（1.199 g/L）略低于原生卤水（1.278 g/L），模拟初期可见其溶剂沿剖面顶部向采卤渠运移.后续受弥散作用影响逐渐与原生卤水逐渐混合并向深部渗透.又因缺乏低渗透区障碍，溶剂能够相对均匀地波及含水层，导致整体渗透系数增加最为显著，模拟第1、3、5年时平均渗透系数分别较初始值增加了14.0%、47.4%和66.2%（图10a）.

相比之下，非均质场景下的渗透系数演化展现出明显的路径依赖性.场景C代表具有强连通性的高渗透网络，ΔK集中在少数几条溶蚀贯通的通道内，形成了清晰的树枝状优先流强化结构.这种高度集中的溶蚀效应使渗透性演化呈现极端的不均衡性，大量低渗透基质区域未能与溶剂有效接触.尽管通道内ΔK极高，但仅占区域总体积的一小部分，导致场景C平均渗透系数增幅有限（图10c）.

场景D模拟了低连通性的渗透场特征，渗透系数增加的区域形态更为曲折.低渗透障碍迫使流体在被分割、连通性差的高渗透块区运移和聚集，增加了溶浸剂在局部区域的滞留时间，深化了水岩反应程度，使渗透系数增幅较大.然而，整体的连通性差阻碍了溶液的推进范围.在此场景下，高渗透区平均渗透系数5年累计增幅为59.4%（图10d）.场景E通过反正弦变换生成具有双峰分布特征的K₀随机场，旨在模拟存在明确的、导水性能差异悬殊的两类介质的极端情况.此时溶浸剂更易沿贯穿性良好的高渗透区快速流动，但相对更宽的通道使得溶剂与矿物接触面积较大，溶蚀作用更强.因此高渗区平均渗透系数在第1、3、5年分别增加了17.5%、60.8%和85.9%（图10e），其幅度在所有非均质场景中最高.

为最终评价不同非均质结构对溶矿效果的影响，采用固液转化效率R和溶矿波及效率S两个指标进行定量评估.在任意时刻t，区域内矿物j的固液转化效率可由下式获得（汪子涛，2024）：

R=\frac{\mathrm{1}}{n}{\displaystyle \sum _{j=\mathrm{0}}^n} \frac{M_j^{\mathrm{0}}-M_j^t}{M_j^{\mathrm{0}}}\times \mathrm{100}\%,(12)

式（12）中，n为单元格总数；M_j ⁰和M_j ^t 分别表示矿物初始与时刻t的含量.波及效率S定义为矿物含量相对变化超过阈值的网格占比，本文分别取5%、10%、15% 3种阈值进行计算与对比.

表6汇总了模拟结束时（第5年）光卤石的R值及不同阈值下的S值.结果表明，初始非均质结构显著调控着最终的溶矿成效.均质分布（场景A）与反正弦分布（场景E）在各阈值下均取得较高的R值和S值，显示含水层初始良好连通性对提高资源采收率的优势.强连通（场景C）与低连通（场景D）两类结构的渗透系数演化路径差异明显，但在3种阈值下的R与S均处于较低水平，表明优先流贯通与低渗阻隔两端都会限制溶采；高斯随机场（场景B）是更普遍的非均质性代表，其溶矿效率则介于上述良好连通与极端结构之间.

综上所述，初始渗透场的非均质性不仅直接塑造了盐湖区人工溶采过程中渗透系数演化的空间形态和路径，更通过控制溶剂的运移效率、分布范围以及与目标矿物的有效接触面积和接触时间，决定了溶采过程的整体潜力与卤水溶采效率.本研究基于耦合模型的精细化模拟与多场景定量评估，为深入理解不同含水层结构下的溶采动态过程与效果提供了理论依据与技术支撑.需要强调的是，不同盐湖的水化学类型、离子浓度及储卤层介质差异可能会影响含水层渗透性演化的速率和幅度，但并不改变非均质性对演化过程的控制作用.

4 结论

本文构建了基于MODFLOW6⁃PhreeqcRM耦合的盐湖地下卤水反应运移数值工具MF6PQC，深入探讨了人工补水溶矿作用下的溶质迁移、水岩作用对储卤层渗透系数演化的控制机制.主要结论如下：

（1）储卤层渗透系数演化呈现明显的路径依赖性.数值模拟结果表明，在人工补水溶矿过程中，渗透系数变化并非均匀发生，而是高度受控于初始渗透场的非均质性.高渗透区域优先受到溶蚀作用，导致渗透系数显著增大，并形成优势渗流通道.这种正反馈机制使得溶浸剂更倾向于在这些高渗区快速运移，限制其在低渗区的扩散和水岩反应，最终导致渗透系数演化呈现出与初始非均质性高度相关的空间格局.

（2）储卤层渗透系数的动态演化受到多因素、多阶段复杂调控.相关性分析与随机森林特征评估发现，在溶采初期，渗透系数剧烈变化主要由水文地球化学反应，特别是光卤石等高活性矿物的快速溶解所主导.随着对流-弥散作用扩大溶浸剂波及范围，初始渗透系数场的空间结构特征及流体密度差异对溶剂运移路径的控制作用凸显，最终决定了渗透系数的长期演变趋势.

（3）含水层非均质空间结构特征对渗透系数演化模式及最终溶矿效率具有决定性影响.不同类型含水层非均质结构导致了迥异的渗透性演化路径.相对均质或具有分布广泛、连通性良好的高渗透区域有利于溶剂均匀波及，能获得较高的矿物固液转化效率.强连通性优先通道或含低渗隔挡的弱连通性结构均会使大量固体矿产资源无法被有效接触和溶解，从而显著降低整体溶采效果.

尽管本研究通过耦合数值模拟能够揭示补水溶矿过程中盐湖卤水反应运移机制与渗透系数的演化规律，但仍存在一定局限性.所采用的假设二维剖面模型难以完全捕捉真实环境的高度复杂性.后续研究将在进一步完善耦合模型的基础上，将其应用于具体盐湖矿区，并进行针对性的现场预测与优化溶采，从而实现对盐湖区固液转换过程的精准预测与高效调控，为卤水资源可持续开发提供更可靠的科学支撑.

说明：本文所开发的建模工具MF6PQC代码及原始数据，均由作者本人上传至Github网站（链接：https：//github.com/wangzitao21/mf6pqc），可供其他研究人员参考与复现.

参考文献

[1] Ahmadi,N.,Muniruzzaman,M.,Sprocati,R.,et al.,2022. Coupling Soil/Atmosphere Interactions and Geochemical Processes:A Multiphase and Multicomponent Reactive Transport Approach. Advances in Water Resources,169:104303. https://doi.org/10.1016/j.advwatres.2022.104303

[2] Chang,W.J.,Yuan,X.L.,Liu,J.B.,et al.,2024.Research on the Changes in Physical Properties during the Seepage⁃Dissolution Process of Potassium Salt⁃Bearing Reservoir in Qarhan Salt Lake.Journal of Salt Lake Research,32(1):99-106 (in Chinese with English abstract).

[3] Chaudhuri,S.,Barceló,M.A.,Juan,P.,et al.,2025.Enhanced Spatial Modeling on Linear Networks Using Gaussian Whittle⁃Matérn Fields.Stochastic Environmental Research and Risk Assessment,39(3):1143-1158.https://doi.org/10.1007/s00477⁃025⁃02912⁃6

[4] Hao,A.B.,Li,W.P.,2003.The Application of Pitzer Theory in the Geochemical Equilibrium Study of High Concentration Brine System with Variable Temperature.Journal of Salt Lake Research,11(3):24-30 (in Chinese with English abstract).

[5] Hong,X.L.,Xia,S.P.,Gao,S.Y.,1994.Dissolution Kinetics of Carnallite.Chinese Journal of Applied Chemistry,11(3):26-31 (in Chinese with English abstract).

[6] Hu,S.Y.,Ren,J.,Li,J.Q.,et al.,2022.Permeability and Brine Enrichment Mechanism of Brine Reservoir in the Mahai Salt Lake,Qaidam Basin.Scientia Geographica Sinica,42(11):2039-2046 (in Chinese with English abstract).

[7] Huang,H.L.,Song,J.,Yang,Y.,et al.,2024.Reactive Transport Numerical Modeling of Typical Heavy Metal Pollutants in Three⁃Dimensional Fracture Networks.Earth Science,49(8):2879-2890 (in Chinese with English abstract).

[8] Hughes,J.D.,Langevin,C.D.,Paulinski,S.R.,et al.,2024.FloPy Workflows for Creating Structured and Unstructured MODFLOW Models.Ground Water,62(1):124-139.https://doi.org/10.1111/gwat.13327

[9] Kaufmann,G.,2016.Modelling Karst Aquifer Evolution in Fractured,Porous Rocks.Journal of Hydrology,543:796-807.https://doi.org/10.1016/j.jhydrol.2016.10.049

[10] Langevin,C.D.,Provost,A.M.,Panday,S.,et al.,2022 Documentation for the Modflow 6 Groundwater Transport Model:6⁃A61.U.S.Geological Survey,U.S.A..

[11] Li,H.,2010.Analysis of Environmental Impact Due to Exploitation of K⁃Li⁃B Mine and Optimization of Exploitation Scheme in Qinghai West Taijinar Salt Lake (Dissertation).Changan University,Xi’an(in Chinese with English abstract).

[12] Li,H.K.,Lu,C.H.,Werner,A.D.,et al.,2022.Impacts of Heterogeneity on Aquifer Storage and Recovery in Saline Aquifers.Water Resources Research,58(5):e2021WR031306.https://doi.org/10.1029/2021wr031306

[13] Li,J.S.,Li,T.W.,Ma,Y.Q.,et al.,2022.Distribution and Origin of Brine⁃Type Li⁃Rb Mineralization in the Qaidam Basin,NW China.Science China Earth Sciences,65(3):477-489.https://doi.org/10.1007/s11430⁃021⁃9855⁃6

[14] Li,M.L.,2023.Dynamic Characteristics of Brine and Its Indicative Effect on Solid⁃Liquid Transformation Process in the Qarhan Salt Lake (Dissertation).Qinghai Institute of Salt Lakes,Chinese Academy of Sciences,Xining(in Chinese with English abstract).

[15] Li,R.Q.,Liu,C.L.,Jiao,P.C.,et al.,2020.The Present Situation,Existing Problems,and Countermeasures for Exploitation and Utilization of Low⁃Grade Potash Minerals in Qarhan Salt Lake,Qinghai Province,China.Carbonates and Evaporites,35(2):34.https://doi.org/10.1007/s13146⁃020⁃00562⁃z

[16] Li,R.Q.,Liu,C.L.,Jiao,P.C.,et al.,2021.Numerical Simulation Analysis on Potash Dissolution Extraction from Low-Grade Solid Potash Ore in Modern Salt Lake:A Case Study from Bieletan Mining Area in the Qarhan Salt Lake.Acta Geologica Sinica,95(7):2150-2159 (in Chinese with English abstract).

[17] Lu,P.,Zhang,G.R.,Apps,J.,et al.,2022.Comparison of Thermodynamic Data Files for PHREEQC.Earth⁃Science Reviews,225:103888.https://doi.org/10.1016/j.earscirev.2021.103888

[18] Luo,Z.J.,Wang,X.,Dai,J.,et al.,2024.Influence of Land Subsidence on Minable Groundwater Resources.Earth Science,49(1):238-252 (in Chinese with English abstract).

[19] Müller,S.,Schüler,L.,Zech,A.,et al.,2022.GSTools V1.3:A Toolbox for Geostatistical Modelling in Python.Geoscientific Model Development,15(7):3161-3182.https://doi.org/10.5194/gmd⁃15⁃3161⁃2022

[20] Pardo⁃Iguzquiza,E.,Chica⁃Olmo,M.,2008.Geostatistics with the Matern Semivariogram Model:A Library of Computer Programs for Inference,Kriging and Simulation.Computers & Geosciences,34(9):1073-1079.https://doi.org/10.1016/j.cageo.2007.09.020

[21] Parkhurst,D.L.,Wissmeier,L.,2015.PhreeqcRM:A Reaction Module for Transport Simulators Based on the Geochemical Model PHREEQC.Advances in Water Resources,83:176-189.https://doi.org/10.1016/j.advwatres.2015.06.001

[22] Patrick, A. D.,Franklin, W. S.,1998.Physical and Chemical Hydrogeology.Wiley,New York.

[23] Qaidam Comprehensive Geological and Mineral Exploration Institute,2023.Research Report on Recoverable Reserves of Brine Potassium Ore in Qaidam Basin.Salt Lake Administration Bureau of Haixi Prefecture,Qinghai Province,Golmud(in Chinese ).

[24] Rehman,M.,Hafeez,M.B.,Krawczuk,M.,2024.A Comprehensive Review:Applications of the Kozeny-Carman Model in Engineering with Permeability Dynamics.Archives of Computational Methods in Engineering,31(7):3843-3855.https://doi.org/10.1007/s11831⁃024⁃10094⁃7

[25] Ren,Q.H.,2017.The Numerical Simulation of the Evolution Mechanism of Brine under the Supplyment of Humans (Dissertation).Qinghai Institute of Salt Lakes,Chinese Academy of Sciences,Xining(in Chinese with English abstract).

[26] Song,G.,Zhao,Y.Y.,Ma,H.T.,et al.,2024.Research on the Sedimentary Characteristics of Solid Sylvite after Liquefaction in the Bieletan Area of Qaidam Basin.Acta Geologica Sinica,98(10):2873-2882 (in Chinese with English abstract).

[27] Steefel,C.I.,Appelo,C.A.J.,Arora,B.,et al.,2015.Reactive Transport Codes for Subsurface Environmental Simulation.Computational Geosciences,19(3):445-478.https://doi.org/10.1007/s10596⁃014⁃9443⁃x

[28] Sun,Q.M.,Gao,M.S.,Wen,Z.,et al.,2023.Reactive Transport Modeling for the Effect of Pumping Activities on the Groundwater Environment in Muddy Coasts.Journal of Hydrology,621:129614.https://doi.org/10.1016/j.jhydrol.2023.129614

[29] U.S.Geological Survey,2017.Modflow 6 - Description of Input and Output.U.S.Geological Survey,U.S.A..

[30] Vásquez,C.,Ortiz,C.,Suárez,F.,et al.,2013.Modeling Flow and Reactive Transport to Explain Mineral Zoning in the Atacama Salt Flat Aquifer,Chile.Journal of Hydrology,490:114-125.https://doi.org/10.1016/j.jhydrol.2013.03.028

[31] Wang,W.X.,2013.Driving and Dissolving Exploitation of Low⁃Grade Solid Potassium Research in Qarhan Salt Lake(Dissertation).China University of Geosciences,Beijing(in Chinese with English abstract).

[32] Wang,X.,Jiang,H.,Wang,A.L.,2015.Study on the Dissolution Kinetics of K⁺,Na⁺ in Low⁃Grade Salt Mines.Chemistry,78(11):1049-1052 (in Chinese with English abstract).

[33] Wang,Z.T.,2024.Groundwater Modeling and Inverse Problem for Brine in Salt Lakes.Qinghai Institute of Salt Lakes,Chinese Academy of Sciences,Xining(in Chinese with English abstract).

[34] Wang,Z.T.,Yue,C.,Wang,J.P.,2024.Evaluating Parameter Inversion Efficiency in Heterogeneous Groundwater Models Using Karhunen⁃Loève Expansion:A Comparative Study of Genetic Algorithm,Ensemble Smoother,and MCMC.Earth Science Informatics,17(4):3475-3491.https://doi.org/10.1007/s12145⁃024⁃01361⁃z

[35] Weisbrod,N.,Alon⁃Mordish,C.,Konen,E.,et al.,2012.Dynamic Dissolution of Halite Rock during Flow of Diluted Saline Solutions.Geophysical Research Letters,39(9):2012GL051306.https://doi.org/10.1029/2012GL051306

[36] Xie,M.L.,Mayer,K.U.,Claret,F.,et al.,2015.Implementation and Evaluation of Permeability⁃Porosity and Tortuosity⁃Porosity Relationships Linked to Mineral Dissolution⁃Precipitation.Computational Geosciences,19(3):655-671.https://doi.org/10.1007/s10596⁃014⁃9458⁃3

[37] Zhang,Y.,Gao,D.L.,Ren,Q.H.,et al.,2017.Permeability of Unconfined Aquifer under the Dayantan Mine Area near the Kunty Salt Lake in the Qaidam Basin.Arid Zone Research,34(1):36-42 (in Chinese with English abstract).

[38] Zhao,X.F.,Zhao,Y.J.,Jiao,P.C.,et al.,2023.The Reservoir Heterogeneity of Reservoirs and Its Impact on the Brine Recovery in the Yiliping Area,Qaidam Basin.Journal of Salt Lake Research,31(3):29-34 (in Chinese with English abstract).

[39] Zhou,X.,Fang,B.,Chen,M.Y.,et al.,2006.Numerical Simulation of Brine in Crystalline Halite in the Bieletan Area of the Charham Salt Lake Region,Qinghai Province.Arid Zone Research,23(2):258-263 (in Chinese with English abstract).

[40] Zinn,B.,Harvey,C.F.,2003.When Good Statistical Models of Aquifer Heterogeneity Go Bad:A Comparison of Flow,Dispersion,and Mass Transfer in Connected and Multivariate Gaussian Hydraulic Conductivity Fields.Water Resources Research,39(3):2001WR001146.https://doi.org/10.1029/2001wr001146

基金资助

国家重点研发计划青年科学家项目(2023YFC2908600)

青海省“昆仑英才·高端创新创业人才”计划项目(QHKLYC⁃GDCXCY⁃2024⁃049)

国家自然科学基金项目(42502241)

中国博士后科学基金面上资助项目(2024M760016)