基于反向水文地球化学模拟的塔城盆地北区潜水水化学演化机制

崔小顺，郑昭贤，程中双，苏 晨，陈宗宇

（1.河北省煤田地质局环境地质调查院，河北石家庄 050085；2.中国地质科学院水文地质环境地质研究所，河北石家庄 050061；3.自然资源部地下水科学与工程重点实验室，河北正定 050803）

0 引言

塔城盆地位于我国西北端，与哈萨克斯坦国接壤，区内设有我国西部通往中亚及欧洲的交通要道巴克图口岸，是新疆丝绸之路经济带建设的核心区［1］。塔城盆地位于温带干旱和半干旱气候区，地下水是区内水资源的最重要组成部分，是区内社会经济发展的基础资源［2］。地下水化学演化机制是地下水化学组分及其含量形成和变化的主要驱动因素，是控制地下水质的内在成因［3］。科学的水功能区划依赖于对地下水化学形成及演化机制的深入了解，从而避免干旱内陆地区的缺水从资源性进一步向水质性缺水发展［4］。因此，正确认识地下水化学演化机制，可为科学管理和合理开发地下水资源提供科学依据，促进干旱区生态环境与社会经济的协调发展［5］。现阶段，只有少量的区域水文地质测绘和地下水质量评价工作涉及塔城盆地的地下水化学内容。1979年，李庆年［6］等在1∶20万塔城盆地水文地质测绘工作中，首次调查研究了盆地尺度潜水水化学特征；2019-2021 年，吕晓立［7-9］等在塔城盆地地下水质量评价中研究了三氮、铁锰和氟化物的分布特征和主控因素。2019年和2020年，崔小顺［10］等和吕晓立［11］等分别查明了塔城盆地承压地下水和潜水的水化学空间变化特征。上述成果为该区的地下水化学后续研究提供了良好的支撑，但当前仍缺乏从地下水含水和流动系统角度出发，以盆地地下水循环条件为基础，随径流路径进行的区域尺度潜水水化学演化研究。而且，上述研究主要采用图解法、离子组合比值法、多元统计分析法，定性分析和识别含水层中主要地下水化学作用，随典型水文地质剖面的地下水流路径，定量模拟地下水化学演化的研究更为鲜见。

塔城盆地北区是该盆地内一个完整的一级地下水系统，也是区内主要人口聚集区和农业生产区。区内潜水是当地农牧业的主要开采水源。本研究基于前人研究成果，根据塔城盆地北区的水文地质条件，计算潜水蒸发程度，之后采用反向水文地球化学模拟，定量评估地下水系统中发生的水文地球化学作用，从而揭示地下水化学演化机制。这对该区地下水资源的优化配置具有重要意义。

1 研究区概况

塔城盆地的地势整体北高南低，自北缘地形最高处的塔尔哈巴台山向南缘地形最低处的额敏河，地貌类型呈东西向带状分布，依次为冰水台地、强倾斜洪积平原、黄土岗地、弱倾斜洪积平原和冲积平原。研究区潜水主要为第四系松散岩类孔隙水，自北向南主要赋存于强倾斜洪积平原的上更新统洪积砂砾石和中更新统冰水堆积砾石层；黄土岗地的上更新统含砾亚砂土层和砂质亚黏土层；弱倾斜洪积平原的上更新统及全新统洪积砂卵砾石层；冲积平原的全新统冲积砾质砂及含土砾石层。潜水主要在山前冰水台地及洪积平原接收山区地表潜流及裂隙水的侧向径流补给和少量大气降水入渗补给，大致由北向南先后流经黄土岗地和冲积平原至额敏河阶地，随后排泄至该河。塔城盆地的新构造运动较为强烈。喀拉哈巴克乡以东，发育东西向延伸的断裂翘起，该区石炭—二叠系及新近系基岩埋藏不深，上覆第四系黄色砂质亚黏土，从而在区内形成黄土岗地地貌。由于该断裂翘起对地下水循环的影响，以黄土岗地西缘为界，可将研究区分为东区地下水子系统和西区地下水子系统（以下简称东区和西区）。东区和西区的第四系沉积特征有所差异，即整体上东区的第四系沉积物相比西区包含较多的细颗粒物质，西区的颗粒分选性相比东区较差，这使东区的地下水径流条件相比西区较差。西区地下水埋深自山前大于90 m逐渐向额敏河呈平稳减小趋势，基本从也门勒乡起，下游地下水埋深普遍小于5 m；东区地下水埋深也与西区呈现相似的特征，但受断裂翘起影响，东区地下水位在黄土岗地北缘明显上升，水位埋深小于5 m，地下水径流通过黄土岗地后，随导水断层快速下降，在与之相连的洪积弱倾斜平原形成“跌水”现象，水位埋深再一次大于30 m，之后向下游呈平稳上升趋势，自博孜达克农场南部起，地下水埋深小于5 m（图1）。虽然“跌水”现象使黄土岗地南北测的地下水埋深突变，但南北侧的第四系松散岩类孔隙地下水在导水断层的快速下跌过程中仍具有统一流面，保证了地下水流系统的连续性。

2 研究方法

2.1 样品采集与分析测试

根据地下水流场，本研究分别在东区和西区的典型地下水流路径的始末点取样。地下水样品的采集时间为2015 年8 月，主要采自农田灌溉机井和乡村供水民井。西区模拟路径的始末点地下水均采自上更新统洪积层；东区模拟路径始末点地下水分别采自上更新统洪积层和全新统冲积层，虽然始末点地下水源自不同的含水层，但从东区典型水文地质剖面（图1）可知，以上两个含水层中地下水属同一地下水流系统，始末点地下水源自同一地下水流。

图1 研究区水文地质简图、典型水文地质剖面图及反向地球化学模拟水流路径始末点位置图Fig.1 The hydrogeological map，typical hydrogeological sections and groundwater path simulated in inverse modeling at study area

地下水样品依照样品采集技术指南，经过洗井、现场理化指标监测后进行采集，并密封冷藏保存于聚乙烯瓶。重碳酸根和游离二氧化碳采用滴定法测定；硫酸根、氯离子采用离子色谱仪（IC）测定；钾、钠、钙、镁等离子采用电感耦合等离子体发射光谱仪（ICP-OES）测定；H2SiO3采用分光光度计测定，并检验水中主要离子的电荷平衡误差小于5%，以此保障测试精度。水中氢、氧稳定同位素采用激光同位素光谱仪测定，其测试误差均小于0.1‰。

2.2 反向水文地球化学模型

基于电荷平衡和质量平衡构建的反向水文地球化学模型，主要计算反应相和生成相之间的摩尔转移量，即气体组分或矿物组分离开或进入水溶液的毫摩尔质量，从而确定地下水在同一水流路径上发生的水文地球化学作用［12］。反向水文地球化学概念模型为：

路径起点地下水的水化学组分+反应相=路径终点的地下水的水化学组分+生成相

反应相和生成相主要为参与或形成自水文地球化学反应的矿物、可交换态离子和气体。概念模型可基于电荷平衡和质量平衡的数学模型表达［13］：

式中：m为水文地球化学反应中反应相和生成相的总数；MT，i为地下水中第i种元素的总摩尔浓度；MT，i为地下水中第i种元素的总摩尔浓度；an为第n种矿物或气体离开或进入水溶液的摩尔浓度；bn，i为第n种矿物或气体中第i种元素的化学计量系数；k为模型中包含的元素总数；ΔRS为地下水的氧化还原状态；un为第n个相中的有效价态。

3 结果与分析

3.1 模拟水流路径的选取

为提高模拟结果反映的地下水化学演化机制的代表性，反向水文地球化学建模首先要选取一条典型水流路径。该水流路径需从地下水的补给区横跨至排泄区，而且水流路径上的含水层结构、矿物组成具有区域代表性。按照以上原则，西区选取的路径为QS1→QS2；东区选取路径为QS3→QS4（图1）。模拟路径起始点的地下水化学分析测试结果见表1。

表1 水文地球化学模拟路径始末点的地下水化学分析测试结果Tab.1 The hydrochemical and isotopic results of the initial and final point for groundwater flow simulated by inverse modeling

3.2 反应/生成相的确定

反应∕生成相是水文地球化学作用的直接参与者，设定准确的反应相对于模型求解十分重要。首先，选取含水介质的主要矿物组分（石盐、石膏、方解石、白云石、石英、钾长石、斜长石、云母、蒙脱石、伊利石、高岭石）为本次模型的第一类反应∕生成相［6，14，15］。其次，选取吸附于黏土矿物表面的可交换态阳离子（KX、NaX、MgX2、CaX2）为第二类反应∕生成相。最后，由于潜水处于开放体系，选取可参与水文地球化学反应的游离CO2为第三类反应∕生成相。

在模拟水流路径上，除了反应∕生成相的选取，准确设置水文地球化学反应方向（“溶解”或“沉淀”）可提高模型求解的精度，避免模型输出偏离实际情况的解。采用PHREEQC 软件中SOLUTION 模块，计算模拟水流路径始末点的地下水中矿物饱和指数（SI），以此判定路径上水岩反应方向。从地下水SI 计算表（表2）可知，所有模拟路径的起、终点地下水中钠长石、钙长石、盐岩和石膏的SI均小于0，反映上述矿物组分在地下水中处非饱和态，故在模型中设定上述矿物在模拟路径上始终处于“溶解”；所有模拟路径的起点地下水中钾长石、白云石和方解石的SI均小于0，但在模拟路径的终点地下水中上述矿物的SI均大于0，反映上述矿物在模拟路径的上游区处于非饱和溶解态，但在下游区逐渐饱和转变为沉淀态，故在模型中设定钾长石、白云石和方解石设为“溶解”或“沉淀”；所有模拟路径起始点高岭石、伊利石、钙蒙脱石和云母的SI均大于0，反映上述矿物在地下水流路径上始终处于饱和沉淀态，故在模型中设定上述黏土矿物和云母为“沉淀”。另外，根据崔小顺［16］对塔城盆地潜水阳离子交替吸附作用的研究结果，设定NaX、KX为“溶解”，将相应的阳离子交替吸附产物，即CaX2、MgX2设定为“沉淀”。

表2 模拟路径起始点地下水的矿物饱和指数（SI）Tab.2 Saturation indices data of the initial and final point for simulating groundwater flow

3.3 模拟路径上潜水蒸发程度定量评估

除了水岩作用，研究区潜水的水化学组分在径流路径上还受蒸发浓缩作用的影响，有必要在模型中设定水流路径上地下水的蒸发倍数。地下水中氢氧稳定同位素可表征蒸发作用和水岩作用引起的同位素交换［17］。郝爱兵等［18］指出在我国西北的含水介质和第四系孔隙水之间基本不产生由水岩作用引发的氢氧稳定同位素交换。因此，地下水中氢氧同位素只受蒸发作用影响。基于瑞利分馏推演的氘盈余定量蒸发模型［19，20］，定量评估模拟路径始末点的地下水蒸发度，即经蒸发作用后地下水与地下水补给来源—大气降水的体积比。地下水点蒸发度和模拟路径蒸发倍数的具体数学公式如下：

式中：f、ffinal、finitial为地下水蒸发度、模拟路径终点和起点地下水蒸发度；dL0、dL为蒸发后大气降水和地下水的氘盈余；α1、α2为开放系统中非平衡分流条件下的气相对液相的2H、18O 分馏系数；k为当地大气降水的斜率；e为模拟路径的地下水蒸发倍数。

参数α1、α2依赖于地下水温度和大气湿度，具体参数赋值见参考文献［20］。

根据研究区大气降水的算数平均氢氧同位素比值（-13.3‰、-92.5‰；n=131）、塔城盆地地下水年均温度（13.3 ℃）和年均大气湿度（50%）、当地大气降水线（y=7.2x+3.2；n=131、R2=0.98）［10］，经计算得出模拟路径始末点地下水的蒸发度分别为95.9%（QS1）、89.8%（QS2）、98.8%（QS3）和90.8%（QS4）；模拟路径1、2的地下水蒸发倍数分别为0.936和0.919。

3.4 模型参数设置

设定完模拟地下水流路径、模拟路径上反应矿物相的种类及反应方向和模拟路径上地下水的蒸发倍数之后，可用PHREEQC 软件建立反向水文地球化学模型。本次模型只选取模拟路径始末点之间浓度变化较大的Na+、Mg2+、Ca2+、HCO3-、Cl-和设为模型收敛的变量条件，并设置实测值和模拟值的误差小于5.0%。模拟水流路径上的地下水温均小于20℃，故该模型选择适用于该温度区间水岩反应体系的热力学数据库phree‐qc.dat。

3.5 模型最优解筛选

由于模拟选取的潜在反应∕生成相大于实际参与水岩反应的相，基于PHREEQC 的反向水文地球化学模型，经过求解模拟路径上反应∕生成相的质量平衡方程和电荷平衡方程［公式（1）、（2）］，输出所有符合平衡方程及反应热力学规律的可行解，但其中包含不符合实际水文地质条件的无效解，有必要筛选最优解。模拟路径1 的反向水文地球化学模型共输出了5 个可行解，见表3。可行解4和5中表明在水流路径上发生了大量的钠长石的不全等溶解（钠长石的高岭石化），其反应量甚至大于易溶的石膏。另外，可行解4 和5 中未发生正向阳离子交换（Na-X+Ca2+=Ca-X+2Na+），这在产生大量石膏溶解而释放钙离子的地下水中不合理。另外，作者在前期研究中对135 组均匀分布于研究区潜水的氯碱指数（CAI）、溶解性总固体（TDS）和地下水化学类型的相关性进行了定性分析，结果表明潜水自补给区至排泄区均发生正向阳离子交换［16］。因此，判定可行解4、5 不符合模拟路径上实际水文地球化学过程。可行解1、2、3 中石膏、石盐、白云石、方解石、长石的溶解量，阳离子交换量和高岭石的沉淀量基本相同，只区别于发生不全等溶解的长石种类。从反应热力学角度，在0.1 MPa、25 ℃的环境条件下，钠长石高岭石化反应的吉布斯自由能量增量（ΔG）（-77.84 kJ∕mol）小于钾长石高岭石化反应的相应值（-43.48 kJ∕mol）［21］。这反映地下水系统中钠长石相比钾长石更易发生不全等溶解，但可行解3 中钾长石的高岭石化反应优先于钠长石相应反应发生，这在钠长石和钾长石同时存在的含水介质中不合理。综上，可行解1 和2 的结果符合野外实际情况，可将其确定为模拟路径1 模型的最优解。

从模拟路径1 的最优解可知，方解石在整体径流路径上发生了沉淀，这是由于下游径流区和排泄区含水介质中石膏发生了大量溶解，使该段地下水中Ca2+浓度显著升高，并引发了同离子效应，导致地下水中CaCO3达到过饱和态，并沉淀析出方解石。这结果也反映在整个模拟路径的水岩相互相互作用体系中方解石的沉淀量大于溶解量。正向阳离子交换作用发生于研究区潜水的整个径流路径，并降低地下水中Mg2+浓度［16］，因此白云石溶解是地下水中Mg2+浓度升高的主控作用，该作用对潜水相应水化学组分的贡献量为1.45 mmol∕L。潜水的蒸发浓缩在整个模拟路径上对水化学组分的贡献有限，潜水中Cl-和的上升主要来自盐岩和石膏溶解，对地下水化学组分的贡献量分别为0.86 mmol∕L 和2.48 mmol∕L。模拟路径上还发生Na+与Ca2+之间的正向阳离子交换，该做用对相应地下水化学组分的贡献量为2.54 mmol∕L（以Na 为计）；另外，由于钠长石的不全等溶解量有限，其对地下水化学演化的贡献也较小。

模拟路径2 的反向水文地球化学模型共输出了4 个可行解，见表3。依照模拟路径1 的最优解筛选原则，排除只发生钾长石不全等溶解反应的可行解3、4，确定可行解1和2为模拟路径2 的最优解。模拟路径2 上对水化学组分贡献量较大的水文地球化学作用与模拟路径1基本一致，主要为石膏、石盐和白云石溶解、方解石沉淀和Na+-Ca2+之间的正向阳离子交换，而且以上作用对水化学组分的贡献量均大于模拟路径1。

表3 反向水文地球化学模型输出可行解统计表 mmol∕LTab.3 The models produced from the inverse modeling

4 讨论

塔城盆地作为我国西北干旱-半干旱区的典型内陆盆地，盆地内降水稀少，缺少自山前向盆地腹地的常年性河流，盆地内地下水在空间上较少受到降水和河水的入渗稀释作用，使地下水化学组分在空间上呈明显的水平分带规律，即潜水的TDS和主要离子自山前向额敏河呈递增特征，潜水水化学类型则从山前HCO3-Ca 型水逐渐过渡为SO4-Na 型水和SO4∙Cl-Na 型水［16］。因此，通过本研究中两个地下水子系统典型水流路径上发生的地下水化学作用，可定量揭示塔城盆地北区一级地下水系统的水化学演化机制。

由于研究区地下水几乎不受降水和河水的稀释作用，主要地下水化学作用为蒸发作用和水岩作用。模拟路径1 和2 的蒸发度分别为0.936 和0.919，表明地下水在径流过程中受到的蒸发作用有限。这是由于潜水在以洪积平原为主的补给区和上游径流区埋深较大，通常大于10 m，从而几乎不发生地下水蒸发作用。另外，虽然潜水在以冲积平原为主的下游径流区和排泄区的埋深小于5 m，满足发生地下水蒸发作用的极限埋深条件，但该区的包气带岩性主要为砂质亚黏土和亚黏土，其垂向渗透系数较低，从而降低了蒸发量。以综合体现地下水中主要离子的水化学指标TDS 为例，仅在以上蒸发程度下模拟路径1和2 的TDS 分别增加11.86 mg∕L 和10.30 mg∕L，只占实际地下水TDS 增量的1.78%和0.60%。因此，研究区蒸发作用对地下水化学组分的影响较小，也反映水岩作用对研究区潜水的水化学演化起到了主控作用。

在地下水流路径上，地下水补给区和上游径流区的潜水水力坡度较大，含水层渗透性较好，径流条件相比下游径流区和排泄区较好（图1）。即使补给区和上游径流区的含水层中赋存石盐、石膏等易溶盐，经过长期充分溶滤，易迁移的离子淋洗比较充分，地下水化学组分主要受碳酸盐和硅酸盐等难溶盐淋滤作用和阳离子交换的影响。由于潜水径流至下游径流区和排泄区时其中碳酸盐处于溶解平衡或饱和状态，该区的地下水化学组分主要受易溶盐溶解的影响。另外，由于下游径流区和排泄区的地下水径流滞缓，水岩反应时间较长，从而使易溶盐溶解对地下水化学组分的贡献量也明显比补给区和上游径流区难溶盐的溶滤作用大。虽然补给区和上游径流区的含水介质普遍比下游径流区和排泄区颗粒粗大，但受碎屑颗粒的分选影响，在整个径流路径上普遍存在黏性土颗粒，具有可发生阳离子交替吸附作用的物质基础。作者的前期研究结果表明潜水的CAI 在整个径流路径上均为负值，反映自补给区至排泄区的地下水中Ca2+、Mg2+持续置换岩土所吸附的Na+、K+，即发生正向阳离子交替吸附作用，而且CAI的绝对值与TDS呈负相关，表明随着地下水流路径阳离子交替吸附作用的驱动力逐渐减弱，反映相比于下游径流区和排泄区，补给区和上游径流区发生较为显著的正向阳离子交替吸附作用。这种规律也符合下游径流区和排泄区石盐溶解引发的地下水中Na+含量快速上升下的水文地球化学环境对阳离子交换作用的约束。在空间上，模拟路径2 代表的东区地下水子系统相比模拟路径1 代表的西区地下水子系统，受断裂翘起形成的黄土岗地影响，其含水介质颗粒整体较细，径流相对滞缓，水岩反应时间的拉长使相同矿物的溶滤对地下水化学组分的贡献量较大。

5 结论

基于定量蒸发程度计算及反向水文地球化学模拟结果的综合分析及讨论，可获得以下结论：

（1）研究区地下水化学演化主要受水岩作用控制，仅受少量蒸发作用影响，几乎不受降水和河水的稀释作用。

（2）从补给区至排泄区的径流路径上，地下水系统中主要发生方解石沉淀，白云石、石膏、盐岩溶解和Na+-Ca2+间的正向阳离子交换，并发生少量斜长石的不全等溶解。

（3）由于西区地下水子系统的径流条件相比东区较好，地下水循环速率较快，溶滤作用、阳离子交替吸附作用对地下水化学演化的贡献量较小，而且在西区地下水子系统基本不发生盐岩溶解作用，这使西区地下水子系统的水质相比东区较好。