不同溶蚀速率对灰岩单裂隙渗流-溶解模型的影响

2019-11-15 02:11

人民长江 2019年10期

(1.长江水利委员会水文局长江下游水文水资源勘测局，江苏南京 210011； 2. 南京水利科学研究院水文水资源与水利工程科学国家重点实验室，江苏南京 210029)

我国岩溶分布面积广，引起的资源环境及工程问题也较为突出，如岩溶区隧洞涌突水、水库渗漏、边坡稳定等问题，而此类问题多与岩溶演化发育有关[1-4]。地下水在运移过程中与可溶性碳酸盐类发生反应，引起矿物溶解，在以裂隙为主的岩溶系统中使得裂隙开度发生扩张，并导致渗流场、化学场等发生不可逆的改变，此过程可通过渗流-溶解耦合模型加以反映刻画[5]。模型主要由水流模块、溶质运移模块、裂隙开度时变模块等组成，并通过各模块变量间的函数关系实现模型耦合[6-7]。其中，矿物溶蚀速率作为模型的重要组成，在以上3个模块中均有所体现。已有研究表明，在不考虑CO2对反应速率影响的封闭环境中，溶蚀速率主要受表面反应和扩散迁移双重因素的联合控制，该两项进程逐次进行，溶蚀速率值为两者间相对较小者[8]。但在实际模拟应用中，由于该速率方程表达较为复杂，常仅考虑受单因素控制，即只受表面反应或扩散迁移其中某一因素控制，从而使得不同研究中溶蚀速率方程有所差异[9-11]。可以看出，以往研究多侧重于研究裂隙系统的发育演变，以及由此带来的渗流动态、溶液组分等的改变，而未涉及在模拟过程中，当溶蚀速率表达不同时，对模拟结果存在何等影响。

针对目前研究中存在的不足，本文旨在量化解析不同溶蚀速率方程对模型结果的影响。首先根据溶蚀方程机理的不同，分为单因素控制(表面反应或扩散迁移)以及两者联合控制。然后假设裂隙形态皆为平行板裂隙，分别模拟获得不同溶蚀速率下裂隙开度、渗流性态以及溶液组分等的空间分布和时间演变规律，并进行比较获得其异同。研究成果可为今后实际应用中根据需要选择更为合适的模型提供指导。

1 溶蚀机理与速率模型

灰岩单裂隙渗流-溶解耦合模型主要由水流模块、溶质运移模块和裂隙开度时变模块3部分组成，具体可见文献[7，10]，在此不再赘述。由于灰岩溶蚀速率是渗流-溶解耦合模型的核心内容之一，在模型3个模块中皆有涉及。根据其受控机理的不同，一般包括以下3种：① 仅受表面反应单一因素控制；② 仅受扩散迁移单一因素控制；③ 受表面反应和扩散迁移双重因素联合控制。以上各机理及速率模型如下所述。

(1) 表面反应控制，即仅受可溶性矿物的溶解动力学控制，其速率大小与溶液中对应矿物的相对饱和度有关。对于碳酸盐类，在接近平衡态和远离平衡态时，其速率表达略有不同。当远离平衡态时反应相对较快，而接近饱和平衡态时，其速率陡降，该过程可表示为

(1)

式中，k1、k2为动力学反应常数；n1、n2为指数，通常n1为等于1或为接近于1的数，而n2在3～6之间；Cs为拐点浓度值，一般为0.7Ceq～0.9Ceq之间。

(2) 扩散迁移控制，即对应离子组分从固相表面进入溶液中的迁移速率，一般认为该速率不仅与溶液中离子浓度有关，还与该位置处裂隙开度b有关，表达如下：

(2)

式中，Dm为离子在水溶液中的扩散系数；ShD为Sherwood数，在层流情形下为8.24[12]；Ceq为离子饱和浓度。

需要指出的是，国内在模拟渗流-溶解耦合问题时主要基于“薄膜假设”，认为紧贴裂隙壁有一层极薄的“薄膜”，内部溶液浓度饱和，薄膜内组分可进入扩散区，导致薄膜内溶液浓度降低；此时，固壁表面发生溶解作用，离子进入“薄膜”补充其内物质损失，使之重新达到饱和；同时，由于固壁表面物质的溶失，造成裂隙开度的增大[13]。该假设本质上也为扩散迁移控制，但在该速率模型中不考虑裂隙开度变化对溶蚀速率的影响。

(3) 表面反应和扩散迁移两因素联合控制，即认为在灰岩溶解过程中，表面反应与扩散过程逐次依序进行，因此总的溶解速率取决于各时刻以上两过程中的最慢者，即以上式(1)与式(2)计算结果间相对较小者作为该时刻下灰岩的溶解速率。

2 不同溶蚀速率对模型影响分析

2.1 模型简介

岩溶地区裂隙相互交错形成裂隙网络，其真实形态极为复杂，因此常以对单裂隙渗流特性的研究作为认识岩体裂隙整体渗流特征的基础。将裂隙壁简化为光滑平板，同时隙宽相对较小时，可采用立方定律研究模型隙宽与单宽流量之间的关系，在此基础上加以修正可推广至粗糙裂隙及裂隙网络。因此，平行板裂隙下渗流特性的研究作为真实裂隙网络的基础具有较为重要的意义[14-15]。

考虑研究区为一长为L、宽为W的单裂隙，其初始裂隙隙宽为b0，裂隙两端水头分别为Hin和Hout，流入溶液的浓度为Cin，裂隙内初始水头和初始浓度分别为H0和C0，模型示意图如图1所示，具体取值见表1。

图1 模型示意Fig.1 Sketch of model

L/mW/mαLαTM/(g·mol-1)ρ/(kg·m-3)Cin/(mol·m-3)C0/(mol·m-3)b0/mmH0/mHin/mHout/m10.50.050.0210026600.0500.200.10

2.2 溶蚀速率参数及比较

表面反应和扩散迁移控制中涉及到的主要参数如表2所示。为研究不同溶蚀速率方程对模型结果的影响，分别建立以下模型，其溶蚀速率对应控制因素为：模型1为仅表面反应单一因素控制，模型2为仅扩散迁移单一因素控制，模型3为表面反应及扩散迁移两者联合控制。

表2 反应速率方程中相关参数Tab.2 Parameters in reaction rate equation

注：表中参数引自文献[8]。

在有限元软件COMSOL Multiphysics软件中建立尺寸为L×W(长×宽)的二维研究区域，并设定相应的模型方程和定解条件。采用瞬态分离式求解器( time dependent segregated) 对以上3个模型分别进行求解。在其它条件相同情形下，模拟时段(2 000 d)内选取500,1 000,1 500 d和2 000 d，在上述时刻时各模型裂隙开度变化如图2所示。由图2可以看出：

(1) 3个模型在上游进水口位置溶蚀皆较为严重，而向下游溶蚀逐渐减缓，从溶蚀锋面来看，在前期(500 d和1 000 d)皆较为平直，而在中后期(1 500 d和2 000 d)出现溶蚀差异不均衡现象。

(2) 3个模型的不同主要体现在溶蚀“程度”和“深度”两方面。从溶蚀“程度”来看，在研究区上游侧，尤其是靠近进水口部位，模型1结果较模型2和模型3明显偏大，而在研究区的中后部，模型1和模型2未发生明显溶蚀，其开度b较模型3明显偏小；从溶蚀“深度”看，模型3>模型2>模型1，其中模型3于2 000 d时在研究区靠近轴部的部位形成了一条较宽的集中渗流通道，且基本贯穿整条裂隙，模型2则形成了局部的渗流通道，达到了裂隙中部，但未贯穿，而模型1的溶蚀仅发生在靠近进水口的较短距离内(<0.2 m)。

图2 不同时刻下裂隙开度b的分布(由上至下依次为500,1 000 ,1 500 d和2 000 d)Fig.2 Distribution of fracture aperture at different time(from top to bottom are in 500,1 000,1 500 d and 2 000 d)

为更直观地比较3个模型各时刻下开度在沿裂隙延伸方向上的变化，以y=0.25 m作为典型断面，以上时刻各模型开度b在该典型断面上的分布见图3。

由图3可以看出：相同时刻下，在裂隙前端靠近进水口部位(<0.1 m)，模型2和模型3的开度分布较为近似，且都明显小于相同位置处模型1的。这表明模拟时段较长时，该部位溶解主要受扩散迁移控制，而采用表面反应控制的动力学表达，会导致该部位灰岩溶蚀的过量估计。在研究区中后段(>0.4 m)，尤其是模拟时刻的后期(如2 000 d)，模型1和模型2的结果明显小于模型3的，表明单一因素控制的模型在研究区范围相对较大时，其位于研究区中后部的结果误差相对较大。

可通过计算获得模拟时段内研究区平均隙宽随时间的变化。平均隙宽一方面可反映裂隙形态的变化，另一方面也可反映灰岩裂隙整个溶失量的变化。同时，为更直观地比较各模型间的差异，以模型3的结果作为基准，计算获得模型1、模型2与模型3间的比例。各模型平均隙宽及对应比例随时间变化见图4。

由图4可以看出，3个模型计算获得的实际平均隙宽在整个模拟时段内皆呈增大的趋势，但彼此间的大小关系在不同时刻有所不同。在该算例中，各模型间的差异性变化可分为以下4个阶段：① 300 d之前，各模型的平均隙宽表现为模型2>模型1>模型3；② 300～510 d之间，表现为模型1>模型2>模型3；③ 在510～1 370 d之间，则为模型1>模型3>模型2；④ 在1 370 d之后，模型3>模型1>模型2，且这段期间模型3与其他两个模型间的差异越来越大，在模拟结束时模型1和模型2仅为模型3值的0.739和0.593。以上4个阶段可简单概括为：前期单一因素控制下模型平均隙宽变化大于两因素联合控制的；而后模型3(两因素联合控制)的隙宽变化依次超过模型2(扩散迁移控制)和模型1(表面反应控制)，且与后两者间的差异越来越显著。

由以上可知，不同溶蚀速率方程对裂隙开度的变化有所影响，因此通过整个裂隙的流量也会有所差异。同样以模型3的结果作为基准，分别计算获得模型1、模型2与模型3间的比例，则出口位置处各模型通过裂隙的流量Q及对应比例随时间的变化如图5所示。

由图5可以看出：① 3种模型中通过裂隙的流量，总体上随时间的推移皆呈逐渐增长的趋势，而相同时刻下，三者比较，总体上呈“模型3>模型1>模型2”。该结果表明，由单一因素控制溶解速率的模型计算获得的流量结果普遍低于实际通过裂隙的流量，其中尤以模型2的差别相对较大。在290d时，模型2与模型3的差别即超过了10%，而模型1在1080d后才达到该数值。②模型3中流量增长出现“临界点”，在该时刻前增长速率较为平缓，而该时刻后，其速率陡增，并最终达到初始时的96.57倍，此与模型3中裂隙形成贯穿型渗流通道有关；而模型1和模型2在整个模拟时段内流量皆呈缓慢增长的趋势，无明显的流量激增时刻，而终了时刻时的流量也仅为初始时的3.22倍和1.83倍。以上结果表明，在对通过裂隙的渗流量进行预测时，采用单因素控制的模型会低估通过裂隙的渗流量，尤以扩散迁移控制的模型表现明显。

图3 不同时刻下y=0.25 m处裂隙开度变化Fig.3 Fracture aperture along y=0.25 m

图6 不同时刻下研究区内Ca2+浓度的分布(由上至下依次为500，1 000，1 500 d和2 000 d)Fig.6 Distribution of Ca2+ concentration in 500,1 000,1 500 d and 2 000 d

图4 不同模型裂隙平均隙宽及相应比例Fig.4 Average aperture width and corresponding proportion of different models

图5 不同模型通过裂隙流量及比例过程线Fig.5 Flow through aperture of different models and corresponding proportion curve

除裂隙开度、渗流场外，模型间化学场分布亦有所不同，其中Ca2+分布如图6所示。由图6可以看出各模型计算获得的Ca2+分布具有以下异同点。

(1)对于每个模型而言，Ca2+浓度沿裂隙延伸方向逐渐增大，在进水口位置处等于流入溶液的Ca2+浓度，而在一段距离后接近饱和；且与裂隙开度分布相似，锋面在前期较为平直，后在中后期出现溶蚀差异。

图7 不同时刻y=0.25 m处Ca2+浓度分布Fig.7 Concentration of Ca2+ along y=0.25 m

(2)模型3由于在模拟时段后期在近裂隙轴部部位形成较为集中的渗流通道，且贯穿整个裂隙，因此在该段时期裂隙溶液中Ca2+浓度相对较小，而相同时段内模型1和模型2该部位的Ca2+浓度仍相对较大。

(3)在模拟时段前期(<1 500 d)，尽管3个模型在裂隙中后部的溶液Ca2+浓度皆接近饱和，但其值大小仍有区别，其中模型2的浓度值已基本等于灰岩溶解的饱和浓度值，即0.1768 mol/m3；而模型1和模型3的最大值仍小于0.16 mol/m3，这与表面反应控制在溶液达到拐点浓度后，其反应速率大大减慢有关。

为更清晰地比较各模型中Ca2+浓度沿裂隙延伸方向上的变化，以y=0.25 m作为典型断面，则各模型在该断面上Ca2+浓度分布如图7所示。由图7可以看出：

(1)无论是单一因素控制的模型1和模型2，还是双重因素控制的模型3，Ca2+浓度在沿裂隙延伸方向上皆呈逐渐增大的趋势，且在裂隙前部增长较快，而在裂隙中后部增长缓慢。

(2)模型1与模型3在裂隙中后部Ca2+的浓度值呈缓慢增长的趋势，但仍未达到饱和状态，而模型2则达到饱和平衡态。

3 结论

在建立灰岩单裂隙渗流-溶解耦合模型的基础上，研究了不同溶蚀速率对模型结果的影响，并以一平行板裂隙中渗流溶解为例，得出如下结论。

(1) 对于裂隙开度而言，从溶蚀深度来看，模拟前期3种模型基本一致，而中后期双重因素下模型明显更深，从而形成贯穿裂隙的渗流通道；就溶蚀程度来看，3种模型下进水口位置溶蚀更为严重，其中尤以表面反应单因素控制的更为明显，其在上游位置的溶蚀明显大于另外两种模型。

(2) 就总的溶失量而言，前期单一因素控制下模型平均隙宽变化大于两因素联合控制的，但相差不大；而后两因素联合控制依次超过扩散迁移控制和表面反应控制，且与后两者间的差异越来越显著。

(3) 从通过裂隙的渗流量来看，采用单因素控制的模型会低估通过裂隙的渗流量，尤以扩散迁移控制的模型表现明显。同时，溶蚀使得裂隙中可形成贯穿型集中渗流通道，从而引起流量陡增，单因素控制的模型预测的该临界点远远晚于实际。

(4) 从化学场分布来看，模拟前期3种模型差别并不显著，但随着两因素联合控制的模型形成贯穿型渗流通道，其裂隙内Ca2+浓度分布明显小于另外另种模型。此外，对于扩散迁移控制的模型，其在裂隙下游侧溶液中Ca2+浓度达到饱和，其值明显大于另两种模型。

(5) 当采用单因素控制对溶解速率进行简化时，会对模型结果造成一定影响，但可根据应用中关注的内容(如裂隙开度、溶蚀量、渗流量、Ca2+浓度等)以及模拟时段的长度，从表面反应和扩散迁移中选择合适的模型，亦可满足一定的精度要求。