2020, Vol. 41
地下水是一种含有不同浓度离子成分的溶液, 地下水在岩体中运移过程中, 各种离子运移受到对流、扩散以及水岩化学反应三种机制控制.其中, 对流是由压力差引起的流体流动；扩散是由浓度差导致的溶质运移；水岩化学反应则是水溶液与矿物岩石间物质成分的相互交换作用的化学反应.裂隙中溶质运移使得溶质浓度分布不断变化, 影响岩体工程所处地下水环境, 同时, 裂隙开度和壁面形状不断发生变化, 影响岩体工程安全性[1].
多孔介质中反应性溶质运移研究较为充分.Moore等通过试验的方式得到了在流固热化学耦合条件下花岗岩裂隙渗透性的变化情况[2]；Wu等建立了动态水力传导率、动态饱和度及两者的动态模型[3]；Li等研究表明裂缝中的对流对溶质浓度的分布有很大的影响, 说明离散裂缝对加速物质的运移起着重要的作用[4]；Hosseini等推导了裂隙与周围基质间质量交换的积分控制方程, 采用适当的富化函数, 建立了裂缝性孔隙域的扩展有限元模型[5].Li等采用拉普拉斯变换和数值拉普拉斯逆变换方法求解分析了蒙皮效应对空间浓度分布和穿透曲线的影响[6].Yang等研究了基于变分原理, 提出了一种新的溶质输运理论, 并发现了溶质输运的两种模式:移动及升降模式[7]；Caceres-Jensen等通过对10种典型火山灰土的吸附试验, 研究了草甘膦的吸附动力学和吸附-解吸过程[8]；Bottacin-Busolin提出了一种基于物理的一维溶质运移模型, 该模型将弱对流混合描述为随深度指数衰减的扩散过程[9]；Kuva等提出了一种时间域-随机游走方案来数值逼近扩散方程的解[10]；Morales-Hernández等将二维溶质输运方程加入到二维浅水方程中, 求解耦合方程组中的流动和溶质相互作用[11].Zhao等的研究表明应力不仅改变了溶质在裂隙网络中的停留时间, 而且改变了溶质的运移路径, 当水力梯度较小时, 基质扩散在溶质运移中起主导作用[12]；Qiao等研究了地下水与岩体间水物理和水化学作用对矿物成分和岩石结构变化的影响[13].随着研究的不断深入, 裂隙中反应性溶质运移得到了较多关注.Grisak等建立了单裂隙岩体中溶质运移模型，考虑对流、扩散、化学反应方面的数学模型[14]；Tsang等研究了裂隙延展方向的扩散作用、弥散作用以及壁面的吸附作用和基质内部分子扩散作用, 同时建立了单裂隙岩体溶质运移数学模型, 推导出方程的解析解[15]；严小三等利用自行研制的粗糙裂隙和光滑平行裂隙试验仪器, 对裂隙内非反应性溶质和反应性溶质的运移规律进行研究[16-17]；Dijk等研究了光滑单裂隙条件下反应性溶质运移问题, 并总结了沉淀过程中反应性溶质浓度、裂隙开度等的变化情况[18].
上述研究分析了特定条件下裂隙中反应性溶质运移特征, 未系统地研究不同机制主导作用下裂隙溶解过程中反应性溶质运移特征和粗糙度对运移特征的影响.为此, 本文开展了不同机制主导作用下裂隙溶解过程中反应性溶质运移特征, 并分析粗糙度的影响.
1 研究方法 1.1 几何模型为分析粗糙度对裂隙溶解过程影响, 选取了4种不同粗糙度几何模型, 其粗糙度系数[19](JRC)分别为0, 4.30, 14.25, 19.98, 如图 1所示, 裂隙长度均为1 m, 平均几何开度设为5 cm.
|
图 1 裂隙初始形状 Fig.1 Initial fracture shape. (a)—JRC=0；(b)—JRC=4.30；(c)—JRC=14.25；(d)—JRC=19.98. |
本文所研究的流体是等温且易混合的, 假设溶液的黏滞度和密度不变, 模型的运移机制包括对流、扩散以及溶质与裂隙壁间的水岩化学反应.为了便于研究, 本文采用一阶不可逆的动力学表面反应模型.在这里反应速率可以用一阶动力学反应速率常数k1表示.裂隙反应性溶质运移基本物理方程由对流扩散、水岩化学反应以及流速变化方程组成.该系统溶液溶质质量守恒方程(1), 裂隙壁溶解质量平衡方程(2), 以及流速变化方程(3)为
|
(1) |
|
(2) |
|
(3) |
式中:Dm表示溶质扩散系数；umax表示局部最大横向速度；dφ/dx表示沿裂隙方向上的局部水头梯度；μ是动态流体黏滞度系数；ρs是固体壁面的密度；ρ为溶液密度；csat是饱和流体的溶质浓度；b是裂隙开度；k1是一阶动力学反应常数.
为表示不同机制在反应性溶质运移中所起到的作用, 引入贝克来数Pe和达姆科勒数Da:
|
(4) |
|
(5) |
其中:Pe表示溶质运移过程中对流所占的比重；Da表示运移过程中水岩化学反应所占的比重.
1.3 模型参数根据Dijk等[18]的研究, 溶液密度ρ和黏滞度系数μ分别为103 kg/m3和10-3 Pa·s, 水力梯度的大小一般为10-4~10, 溶质的自由水分子扩散系数Dm的大小是10-9 m2/s, 相应的反应动力学系数k1一般范围是10-12~10 s-1, 饱和浓度csat设为1 mol/L, 入口浓度cin设为0.
1.4 分析工况Dijk等[18]根据自然界中裂隙流动速度、开度范围和主要矿物反应常数, 将Pe和Da分为低、典型和高三种基本情况, 通过分析进一步得出4种典型工况:典型Pe、低Da；高Pe、典型Da；典型Pe、典型Da；典型Pe、高Da.前两种工况对流占主导, 后两种工况化学反应占主导.依据上述结果本文研究工况如表 1所示.
|
表 1 粗糙裂隙不同工况中模型参数取值 Table 1 Parameter values of rough fracture in different cases |
为验证分析方法与模型参数, 文献[20]首先模拟了光滑单裂隙的沉淀特征, 通过反分析发现这种研究方法与Dijk等[18]所得结论相同, 沉淀模拟结果吻合, 所以验证了本文分析方法的可靠性.
2 裂隙溶解过程特征由于4种粗糙度变化规律相似, 这里以初始JRC为14.25的裂隙模型为例进行分析.图 2为不同工况下裂隙溶解厚度与浓度分布图, 不同颜色反映了不同浓度分布, 黑线代表初始裂隙形状.图 3、图 4、图 5分别为各种工况下裂隙中溶质浓度、裂隙开度和溶液流速变化曲线.
|
图 2 不同工况下裂隙溶解厚度与浓度分布图 Fig.2 Dissolution depths and ion concentrations for fractures in different cases (a)—典型Pe、低Da；(b)—典型Pe、典型Da；(c)—典型Pe、高Da；(d)—高Pe、典型Da. |
|
图 3 粗糙裂隙浓度分布图 Fig.3 Concentration distribution of rough fracture (a)—典型Pe、低Da；(b)—典型Pe、典型Da；(c)—典型Pe、高Da；(d)—高Pe、典型Da. |
|
图 4 粗糙裂隙开度变化图 Fig.4 Aperture distribution of rough fracture (a)—典型Pe、低Da；(b)—典型Pe、典型Da；(c)—典型Pe、高Da；(d)—高Pe、典型Da. |
|
图 5 粗糙裂隙流速变化图 Fig.5 Velocity distribution of rough fracture (a)—典型Pe、低Da；(b)—典型Pe、典型Da；(c)—典型Pe、高Da；(d)—高Pe、典型Da. |
由图 2a可以看出典型Pe、低Da工况溶解时上、下游(图 2中左侧为上游，右侧为下游)溶解量相当, 基本保持均匀变化，由图 3a, 4a, 5a可以看到这种工况中的溶质浓度、裂隙开度、裂隙中流速的分布和变化情况.裂隙内的溶质浓度随时间和位置基本保持为零不变, 其原因是由于Da较低, 化学反应速率缓慢, 这时对流占主导作用, 由于流速相对较快, 溶解生成的溶质立刻被带走, 故浓度一直保持为零.开度基本呈均匀变化, 其原因与浓度的均匀分布有关.速度随开度的变化而变化, 故产生图示分布状态.这种情况溶解速率缓慢, 是4种工况中对岩体工程影响最小的情况.
2.2 典型Pe、典型Da由图 2b可以看出，典型Pe、典型Da工况溶解时上游溶解速度明显高于下游, 由图 3b可以看出，裂隙内的溶质浓度随着时间变化逐渐上升, 因为这时的化学反应速率与对流都发挥一定的作用, 所以溶质沿着延展方向运移并不断积累, 导致产生图示浓度分布.由图 4b、图 5b可以明显看出，上游溶解速率高于下游, 主要原因是浓度的不均匀分布导致裂隙开度不均匀变化, 上游溶解速率不断加快, 到后期入口处溶解速率急剧上升, 最终可能导致裂隙在入口处被完全溶解.
2.3 典型Pe、高Da图 2c为典型Pe、高Da工况溶解浓度云图, 入口处迅速溶解, 从图 3c看出在入口处迅速升高并维持饱和浓度, 原因为高Da下的化学反应速率极高, 这时溶质在进入裂隙的一瞬间迅速溶解，并产生大量溶质而充满裂隙, 直接使溶液饱和.图 4c、图 5c中开度也在入口处迅速增大, 其余基本保持初始开度不变, 原因同样是入口处反应剧烈, 而裂隙内溶质饱和不再发生反应.由于短时间内裂隙入口就可能被完全溶解, 所以这种工况对岩体工程的威胁最大, 实际工程中应加以防范.
2.4 高Pe、典型Da由图 2d可以看出，这种工况溶解时上下游溶解量相当, 但从时间上看它超前于典型Pe、低Da工况.从图 3d、图 4d、图 5d可以看出，这种工况与第一种工况相似, 浓度维持零不变, 但其开度在短时间内呈均匀变化, 但到反应后期, 呈现出开度大的地方变化速率更快的不均匀变化现象, 而典型Pe、低Da在更长的时间后才会出现这种现象, 因此这种情况下其可能更快地从开度最大的位置贯通, 所以也会对岩体工程造成影响.
3 裂隙粗糙度对溶解过程影响将JRC为0(光滑单裂隙)作为参照模型, 将另外三种粗糙裂隙模型的单位时间溶解量、裂隙内溶液平均流速与其对比, 研究不同粗糙度对溶解量变化及流速的影响, 如图 6、图 7所示.
|
图 6 粗糙度对溶解量的影响关系图 Fig.6 Effect of roughness on mass of dissolution |
|
图 7 粗糙度对流速的影响关系图 Fig.7 Effect of roughness on velocity |
这里由于典型Pe、高Da工况下溶解形式极度不均匀, 并且溶解也仅仅发生在入口处, 故其无法反映裂隙粗糙度的影响, 这里只讨论其他三种情况.
3.1 粗糙度对溶解量的影响图 6所示典型Pe、低Da工况中, 同一时间段光滑单裂隙模型溶解量为0.043 61 m2, JRC为4.30，14.25，19.98的溶解量分别为0.044 44，0.046 22，0.051 24 m2, 溶解速度分别是参照模型的101.91 %，104.00 %，110.86 %, 因而可以看出这种工况下溶解速率随粗糙度的增大而加快.主要原因为各个位置裂隙开度不同、流速不同而影响化学反应速率, 同时由于裂隙粗糙度越大, 裂隙反应面面积越大, 单位时间消耗溶质越多, 溶解量变大.所以粗糙度的增大实际上使得岩体安全性降低, 因此实际工程中应尽量避免粗糙度较大的位置, 或采取一定的措施对岩体工程加以防护.同样现象也出现在典型Pe、典型Da和高Pe、典型Da工况中.
图 6所示典型Pe、典型Da工况中, 同一时间段光滑单裂隙模型溶解量为0.005 52 m2, JRC为4.30, 14.25, 19.98的溶解量分别为0.005 63, 0.006 92, 0.007 24 m2, 溶解速度分别是参照模型的101.93 %, 125.35 %, 131.22 %.
图 6所示高Pe、典型Da工况中, 同一时间段光滑单裂隙模型溶解量为0.010 18 m2, JRC为4.30, 14.25, 19.98的溶解量分别为0.010 21, 0.010 30, 0.010 34 m2, 溶解速度分别是参照模型的100.29 %, 101.14 %, 101.53 %.
3.2 粗糙度对流速的影响图 7典型Pe、低Da工况中, 同一时刻光滑单裂隙模型溶液平均流速为4.29×10-15 m/s, JRC为4.30, 14.25, 19.98的平均流速分别为6.29×10-15, 6.43×10-15, 6.63×10-15 m/s, 溶液流动速度分别是参照模型的146.73 %, 150.08 %, 154.72 %, 溶质运移速率随粗糙度的增大而增大.
图 7典型Pe、典型Da工况中, 同一时刻光滑单裂隙模型溶液平均流速为9.13×10-16 m/s, JRC为4.30, 14.25, 19.98的平均流速为1.33×10-15, 1.39×10-15, 1.43×10-15 m/s, 溶液流动速度分别是参照模型的145.69 %, 152.47 %, 156.51 %, 溶质运移速率随粗糙度的增大而增大.
图 7高Pe、典型Da工况中, 同一时刻光滑单裂隙模型溶液平均流速为1.14×10-7 m/s, JRC为4.30, 14.25, 19.98的平均流速分别为1.69×10-7, 2.05×10-7, 2.05×10-7 m/s, 溶液流动速度分别是参照模型的148.49 %, 179.70 %, 179.92 %, 溶质运移速率随粗糙度的增大而增大.
从这三种工况所得结果来看, 流速与粗糙度呈正相关, 也即粗糙度增大后等效水力开度也随之增大.
3.3 裂隙溶液浓度分布分析同一时刻4种不同工况粗糙度裂隙内溶质浓度在反应壁面分布如图 8所示.从图 8a中可以看出典型Pe、低Da类型中裂隙粗糙度越大, 浓度升高幅度越大, 光滑单裂隙的溶质浓度上升幅度在4个工况中属于最慢的一组.原因是粗糙裂隙的粗糙度越大, 溶质与裂隙壁间的反应面积就越大, 这种情况下它们能够更多地与溶质发生反应，导致裂隙内溶质浓度相对更大.
|
图 8 4种工况同一时刻浓度分布 Fig.8 Ion concentration distributions in the four cases (a)—典型Pe、低Da；(b)—典型Pe、典型Da；(c)—典型Pe、高Da；(d)—高Pe、典型Da. |
从图 8b中可以看出, 粗糙度只会对浓度分布产生微弱的影响, 因为这时溶解不均匀, 反应更快, 裂隙变化与初始裂隙形状相关, 与粗糙度相关性小.图 8c中, 粗糙度对浓度分布影响更小, 这时的反应速率最快, 由于裂隙只有在入口处溶解, 所以无法体现粗糙度的影响.图 8d中, 粗糙度越大, 浓度升高幅度越大, 原因同图 8a.
3.4 JRC随时间变化规律图 9所示为4种工况中裂隙粗糙度系数JRC随时间变化规律, 可以看出图 9a中粗糙度随时间呈上升趋势, 因为随着时间的推移, 裂隙在溶解过程中开度变化差异随初始开度的增大而增大, 开度大的位置平均流速相对较低, 导致溶质停留时间增加, 加快了反应进程, 这时这个位置的溶解速率更快, 最终不同位置的开度差异越来越大, 所以出现JRC在反应性溶质运移过程中增大的现象.同样的现象也出现在另外3种工况中, 具体见图 9b，9c和9d.
|
图 9 4种工况JRC随时间变化 Fig.9 JRC change with time in the four cases (a)—典型Pe、低Da；(b)—典型Pe、典型Da；(c)—典型Pe、高Da；(d)—高Pe、典型Da. |
图 10所示为全部9种组合工况溶解形式.第2节分析得出, 典型Pe、低Da工况下, 裂隙开度随时间基本保持不变, 只有在更加漫长的地质年代后才会发生其他变化, 所以这种变化形式被称为均匀变化类型.典型Pe、典型Da工况下, 裂隙开度变化速率在上游快于下游, 产生不均匀溶解, 所以这种变化形式被称为非均匀变化类型.典型Pe、高Da工况下, 裂隙在入口处迅速发生溶解, 而由于反应物的消耗导致整个裂隙其他位置不发生溶解, 故把这种变化形式称之为极度非均匀变化类型.高Pe、典型Da工况下, 裂隙开度初始时刻随时间基本保持不变, 但与均匀变化类型相比, 它在更短的时间后就发生了不均匀溶解, 所以这种变化形式被称为均匀-非均匀变化类型.
|
图 10 不同工况下裂隙溶解形式 Fig.10 Fracture dissolution forms in different cases |
1) 粗糙裂隙溶解过程分为4种类型:①均匀溶解, 对流占主导作用(典型Pe、低Da)；②非均匀溶解, 对流与化学反应相当(典型Pe、典型Da)；③极度非均匀溶解, 化学反应占主导作用(典型Pe、高Da)；④均匀-非均匀溶解, 对流占主导作用(高Pe、典型Da).
2) 裂隙中溶液流动速度较快时, 对流占主导作用, 裂隙内溶质浓度基本保持饱和浓度不变, 前期裂隙均匀溶解, 后期产生不均匀溶解, 溶解速度随时间加快, 同时裂隙粗糙度增大.水岩化学反应常数较大时, 化学反应占主导地位, 裂隙前部溶解量高于后部.
3) 随初始粗糙度增大, 裂隙壁面溶解量增大, 反应性溶质浓度增大, 等效水力开度增大, 裂隙溶解趋势加快.
| [1] |
Wang Z, Glais Y, Qiao L, et al. Hydro-geochemical analysis of the interplay between the groundwater, host rock and water curtain system for an underground oil storage facility[J]. Tunnelling and Underground Space Technology, 2017, 71: 466-477. |
| [2] |
Moore D E, Lockner D A, Byerlee J D. Reduction of permeability in granite at elevated temperatures[J]. Science, 1994, 265(5178): 1558-1561. |
| [3] |
Wu S, Jeng D S, Seymour B R. Numerical modelling of consolidation-induced solute transport in unsaturated soil with dynamic hydraulic conductivity and degree of saturation[J]. Advances in Water Resources, 2020, 135: 103466. |
| [4] |
Li X, Li D, Xu Y, et al. A DFN based 3D numerical approach for modeling coupled groundwater flow and solute transport in fractured rock mass[J]. International Journal of Heat and Mass Transfer, 2020, 149: 119179. |
| [5] |
Hosseini N, Bajalan Z, Khoei A R. Numerical modeling of density-driven solute transport in fractured porous media with the extended finite element method[J]. Advances in Water Resources, 2020, 136: 103453. |
| [6] |
Li X, Wen Z, Zhu Q, et al. A mobile-immobile model for reactive solute transport in a radial two-zone confined aquifer[J]. Journal of Hydrology, 2020, 580: 124347. |
| [7] |
Yang J, Yuan Q, Zhao Y. Solute transport and interface evolution in dissolutive wetting[J]. Science China(Physics, Mechanics & Astronomy), 2019, 62(12): 52-60. |
| [8] |
Caceres-Jensen L, Rodríguez-Becerra J, Sierra-Rosales P, et al. Electrochemical method to study the environmental behavior of glyphosate on volcanic soils:proposal of adsorption-desorption and transport mechanisms[J]. Journal of Hazardous Materials, 2019, 379: 120746. |
| [9] |
Bottacin-Busolin A. Modeling the effect of hyporheic mixing on stream solute transport[J]. Water Resources Research, 2019, 55(11): 9995-10011. |
| [10] |
Kuva J, Voutilainen M, Mattila K. Modeling mass transfer in fracture flows with the time domain-random walk method[J]. Computational Geosciences, 2019, 23(5): 953-967. |
| [11] |
Morales-Hernández M, Murillo J, García-Navarro P. Diffusion-dispersion numerical discretization for solute transport in 2D transient shallow flows[J]. Environmental Fluid Mechanics, 2019, 19(5): 1217-1234. |
| [12] |
Zhao Z, Jing L, Neretnieks I, et al. Numerical modeling of stress effects on solute transport in fractured rocks[J]. Computers and Geotechnics, 2010, 38(2): 113-126. |
| [13] |
Qiao L, Wang Z, Huang A. Alteration of mesoscopic properties and mechanical behavior of sandstone due to hydro-physical and hydro-chemical effects[J]. Rock Mechanics and Rock Engineering, 2016, 50(2): 255-267. |
| [14] |
Grisak G E, Pickens J F. Solute transport through fracture media:1.the effect of matrix diffusion[J]. Water Resources Research, 1980, 16(4): 719-730. |
| [15] |
Tsang D H, Frind E O, Sudicky E A. Contaminant transport in fractured porous media:analytical solution for a single fracture[J]. Water Resources Research, 1981, 17(3): 555-564. |
| [16] |
严小三, 钱家忠, 陈冰宇. 大理石平板裂隙水运移试验与模拟研究[J]. 水动力学研究与进展, 2014, 29(5): 524-529. (Yan Xiao-san, Qian Jia-zhong, Chen Bing-yu, et al. Experimental and simulated study of water flow and transport in a single fracture of marble parallel plates[J]. Hydrodynamics Research and Progress, 2014, 29(5): 524-529.) |
| [17] |
高光耀.非均匀介质中溶质和水分运移模拟及弥散尺度效应研究[D].北京.中国农业大学, 2009. (Gao Guang-yao.Simulation of solute and water transport and dispersion scale effect in heterogeneous media[D].Beijing: China Agricultural University, 2009. ) |
| [18] |
Dijk P, Berkowitz B. Precipitation and dissolution of reactive solutes in fractures[J]. Water Resources Research, 1998, 34(3): 457-470. |
| [19] |
Barton N, Choubey V. The shear strength of rock strength of rock joints in theory and practice[J]. Rock Mechanics, 1977, 10(1/2): 1-54. |
| [20] |
毕竞超.典型条件下裂隙反应性溶质运移特征研究[D].沈阳: 东北大学, 2019. (Bi Jing-chao.Modeling reactive transport in fractures under typical conditions[D].Shenyang: Northeastern University, 2019. ) |