初始裂隙对岩溶水紊流形成的影响

焦友军; 黄奇波; 于青春

doi:10.11932/karst20220401

初始裂隙对岩溶水紊流形成的影响

1.
中国地质大学(北京)水资源与环境学院, 北京 100083

2.
中国地质科学院岩溶地质研究所/自然资源部、广西岩溶动力学重点实验室, 广西桂林 541004

基金项目: 国家自然科学基金项目(41877196, U1612441, 41272387)；中国地质调查项目（DD20221758）

详细信息

作者简介: 焦友军（1990－），男，博士研究生，从事岩溶水资源研究。E-mail：jiaoyj@karst.ac.cn

通讯作者: 于青春（1963－），男，教授，博士研究生导师，从事岩溶水资源研究。E-mail：yuqch@cugb.edu.cn

中图分类号: P641.2

收稿日期: 2022-02-10

刊出日期: 2022-08-25

Influence of initial fractures on the occurrence of karst turbulent flow

1.
School of Water Resources and Environment, China University of Geosciences (Beijing) , Beijing 100083, China

2.
Institute of Karst Geology, CAGS/ Key Laboratory of Karst Dynamics,MNR&GZAR, Guilin, Guangxi 541004, China

More Information

Corresponding author: YU Qingchun, yuqch@cugb.edu.cn

Received Date: 10 February 2022

Publish Date: 25 August 2022

摘要

摘要:
岩溶地区地下发育着大量的溶洞和地下河管道，地下水流状态既有层流也有紊流，而紊流是溶洞管道形成的重要条件。紊流的形成受到岩石初始裂隙的影响，初始裂隙的张开度、分布、走向、迹长、密度等因素都影响着裂隙发育过程中水流状态的变化。通过对不同统计特征的初始裂隙网络进行水流和溶蚀的数值模拟发现，以张开度标准差反映的裂隙网络非均匀性越强，模拟紊流出现的时间就越早；主要裂隙的存在使裂隙网络的非均性增强，主要裂隙与水力梯度总方向的角度越小，紊流出现的时间就越早；当裂隙平均迹长过小时会导致裂隙连通性较差，影响裂隙水流和溶蚀作用；裂隙密度，尤其是主要裂隙密度，对岩溶发育的影响较大。相对于次要裂隙，如果主要裂隙密度偏小，紊流形成时间会大大增加，甚至很难形成紊流。当初始裂隙张开度小于0.001 cm，增大水力梯度仍没有紊流发生，岩溶几乎不发育。
- 初始裂隙 /
- 岩溶管道 /
- 紊流 /
- 溶蚀 /
- 数值模拟
Abstract:
In karst areas of China, there are many karst caves and underground river conduits in the underground aquifer. The development of large cave and conduit are controlled by many factors and it is not every karst aquifer would develop into conduit system, for example, the southern karst area with ample rainfall has many underground rivers, but in the northern arid and semi-arid karst area, underground rivers are relatively rare. The flow states in conduits may be laminar or turbulent. Turbulent flow is an important condition of forming large-scale caves and conduits. In turbulent state, the water flow begin to have the ability of mechanical transportation, which would carry the solid granule and cause the impact and erosion into surrounding rock. The powerful mechanical erosion of turbulence is very important for karst aquifer to develop into large scale conduits and caves. Furthermore, the dissolution rate of carbonate rock in turbulent flow is at least one order of magnitude faster than laminar flow condition. Therefore, the mechanical erosion and chemical dissolution in turbulent flow make the development of conduit and cave more possible. However, the occurrence of turbulent flow is rigidly affected by the hydraulic condition and the initial medium of the aquifer, such as the initial rock fracture, including the aperture, direction, length and density of initial fractures. So we designed different statistic features of fractures and different hydraulic gradients to study the flow state and dissolution widening rates of fractures by numerical simulation. The cubic law and Lomize equation were used to model the laminar and turbulent flow state in fracture. The Newton-Raphson iteration is high-performance to solve the nonlinear flow equation system of laminar and turbulent flow. Then the dissolution rate equation and the Ca²⁺ concentration equation system was employed to model the widening of the fractures. The simulation results and discussions were all under given outer environment with P_CO₂ 0.8% and hydraulic gradient 0.02. The modeling results of fracture aperture showed that when the standard deviation of aperture was 0.0005 cm and the mean of aperture is 0.006, there was no turbulence in the modeling period of 5,000 thousand years. But when the standard deviation was increased to 0.001 cm, the turbulent flow emerged on 189 thousand year and the time was greatly shorten. In the 8 modeling aperture statistic situations, as the mean and the standard deviation of aperture increased and the heterogeneity of fractures was more intensive, the turbulent flow began to appear and the time of turbulence became earlier. The existence of primary fractures led to much heterogeneous aquifer and earlier turbulent time. The results of fracture direction modeling scenes showed that when the angle between the direction of primary fractures and the direction of main hydraulic gradient was smaller, the turbulence time would be shorter. If the mean of fracture length is too little, the connectivity of fracture would become poor and the karst dissolution would be heavily restricted. The fracture density, especially the primary fracture density, had much influence on the karst development. Compared with the secondary fractures, if the density of primary fractures was too smaller, the turbulence time would largely increase, and even no turbulence in the whole simulation. We also discussed the influence of hydraulic gradient on the turbulent time. The hydraulic gradient varied from 0.001 to 1 and the mean of aperture was from 0.001 to 0.005. The results showed that for each aperture, it had the corresponding smallest hydraulic gradient to the occurrence of turbulence and the greater gradient the turbulent time would become earlier. Below the smallest hydraulic gradient the turbulence would never occur. If the mean of aperture is less than 0.001cm, no matter how the hydraulic gradient is increased, the flow state in fractures remained laminar and no turbulence occurred, in which karst is nearly not developed in the aquifer. In conclusion, the occurrence time of turbulent flow reflects the possibility of forming large karst conduits and caves in present aquifer. The shorter the time is, the greater the possibility will be.
- initial fractures /
- karst conduit /
- turbulent flow /
- dissolution /
- numerical modeling

HTML全文

0. 引　言

岩溶地区地下发育着众多裂隙、溶洞和管道，尤其是大型溶洞和地下河，一方面作为景观资源开发可带来经济、文化价值，另一方面作为强烈非均质性的含水空间，给岩溶水资源的研究和开发利用带来巨大挑战^[1-2]。然而，地下河管道和溶洞的形成受到许多因素的控制，并非每个岩溶水系统都可以发育大型的溶洞和管道。比如在南方岩溶区降雨充沛，地下河分布较多，而在北方干旱半干旱岩溶区地下河较少^[2]。即使在同一地区，地下河管道的分布也会受到地层岩性和构造等多种条件的影响^[1-4]。

地下河管道的形成与地下水流运动状态密切相关。当紊流出现时，地下水开始具有一定的机械搬运能力，紊流运动携带固体颗粒对围岩进行撞击和磨蚀^[2]，这种强大的机械侵蚀能力是岩溶管道形成的重要条件。另外在紊流条件下，从岩石表面溶解下来的钙离子进入水中的溶蚀速度相比层流至少高一个数量级^[5]。紊流的机械侵蚀和化学溶蚀能力共同加速了管道的发育^[6]，因此紊流是地下河管道形成的重要条件。紊流条件下模拟裂隙水流的求解从线性问题变成了非线性问题，需要通过迭代法进行求解。Newton-Raphson迭代是一种将非线性方程组线性化迭代求解的经典方法，很多岩溶模拟研究都应用了这种方法进行裂隙水流和溶蚀方程的非线性求解^[7-10]。

岩溶水系统能否形成紊流受到初始裂隙的影响。岩溶水系统的发育受到含水系统边界的水力条件和溶蚀条件控制^[11]，同时与含水介质本身的裂隙分布有关。紊流状态能否出现也受到岩体初始裂隙渗流场的影响。就裂隙介质本身而言，在外部边界条件不变的情况下，初始裂隙的分布越不均匀，裂隙张开度分布差异越大，越有利于优势裂隙的快速发育和紊流的出现^[2]。当岩体较为完整，初始裂隙张开度小于一定程度时，不利于紊流的形成，甚至在裂隙发育过程中不会出现紊流^[6,12]。因此，我们在控制其它影响因素不变的情况下，就初始裂隙在岩溶发育中对紊流形成的影响开展了数值模拟研究。

本研究有助于判断岩溶水系统的岩溶发育程度，对模拟岩溶管道和裂隙网络的内部结构推测起辅助作用，提高岩溶水模拟的精度。

1. 研究方法

1.1 裂隙层流和紊流

岩体裂隙网络是由节点和连通的裂隙段组成的，相邻的两个节点构成单条裂隙，裂隙可通过人工输入或蒙特卡罗方法随机产生。当一条裂隙段为内部裂隙网络的末端裂隙时，其中的水流不能发生流动，因此也不能发生溶蚀。有水流流动的裂隙段为导水裂隙段。因此在水流计算之前需要将末端裂隙删除，末端裂隙和末端节点不参与水流和溶蚀计算。删除不导水的末端裂隙后剩下的即为导水裂隙。裂隙生成和连通的算法在以往的研究中已较为成熟^[3,13]。

在裂隙网络中，任选一个节点作为中心节点i，它周围为相邻节点j，根据裂隙节点水头和张开度计算裂隙单宽流量^[14]。对于单条光滑平行板裂隙水流满足立方定律，用公式（1）来表示。当裂隙水流状态变为紊流时，采用Lomize经验公式^[15]，如公式（2）所示：

${q}_{i,j}=\dfrac{g{{b}_{i,j}}^{3}}{12v}\dfrac{({H}_{i}-{H}_{i,j})}{{L}_{i,j}},\mathrm{R}\mathrm{e} \lt 2300$

(1)

${q}_{i,j}=\dfrac{{H}_{i}-{H}_{i,j}}{\left|{H}_{i}-{H}_{i,j}\right|}4.7{b}_{i,j}{\left[\dfrac{{g}^{4}}{\upsilon }{{b}_{i,j}}^{5}{\left(\dfrac{{H}_{i}-{H}_{i,j}}{{L}_{i,j}}\right)}^{4}\right]}^{\frac{1}{7}},\mathrm{R}\mathrm{e} \gt 2300$

(2)

式中：q_{i, j}为裂隙单宽流量（cm²·s⁻¹）；H_i和H_{i, j}为中心节点i和相邻节点j的水头（m）；L_{i, j}为裂隙段的长度（m）； b_{i, j}为张开度（cm）； v为水的运动粘滞系数，在20 ℃标准大气压下取值为0.01 cm²·s⁻¹；g为重力加速度，取值为980 cm·s⁻²。临界雷诺数Re取值为2 300。

根据水均衡原理，裂隙网络中每个节点水头都可以作为中心节点与相邻节点建立水量平衡方程，并作为一个基本计算单元，如公式（3）所示：

$\displaystyle\sum \nolimits _{j=1}^{n}{q}_{i,j}=0$

(3)

式中：n为中心节点i具有的相邻节点数，n≤4，当其中的1条或2条裂隙为末端裂隙时，不参与计算，此时n取值为2或3。将公式（1）、（2）代入公式（3），得到关于H_i，H_{i, 1}，H_{i, 2}，H_{i, 3}，H_{i, 4}的方程，中心节点周围处于紊流状态的裂隙数N_turb变化从0到4。

所有节点根据公式（3）组成水量平衡方程方程组，当紊流存在时为非线性方程组，需进行迭代求解，本研究采用Newton-Raphson迭代求解：

${{H}}_{1} = {{H}}_{0} + {{f}}_{0} / {{F}}_{0 }$

(4)

${{H}}_{n+1} = {{H}}_{n} + {{f}}_{n} / {{F}}_{n}$

(5)

${{H}}_{{\rm{error}}} = {{H}}_{n+1} - {{H}}_{n}$

(6)

式中：H为所有节点水头向量，f为水头函数向量，按照公式（3）左端水头多项式计算得出。F为导数矩阵，对公式（3）左端多项式关于H的每个元素变量进行求导得到。H₀、f₀、F₀和H_n、f_n、 F_n分别为初始和第n次迭代计算的水头向量、水头函数向量和导数矩阵。H_error为第n+1次和第n次迭代的水头误差向量。假设迭代允许的最大水头误差为errormax，当H_error中所有元素绝对值的最大值小于errormax时，迭代过程达到收敛，H_n+1即为满足收敛标准的水头分布。

1.2 裂隙溶蚀扩宽

1.2.1 单裂隙溶蚀

裂隙中的溶蚀采用如下溶解速率公式^[4]：

$F\left(C\right)={k}_{1}\left(1-\dfrac{C}{{C}_{eq}}\right),\left(C \lt 0.9{C}_{eq}\right)$

(7)

$F\left(C\right)={{k}_{n}\left(1-\dfrac{C}{{C}_{eq}}\right)}^{n},\left(C \gt 0.9{C}_{eq},n \gt 1\right)$

(8)

式中：C为钙离子浓度，C_eq为钙离子平衡浓度，当C<0.9C_eq时采用一阶公式计算溶解量n=1，当C >0.9C_eq时，n=4。在层流条件下溶解速率常数k₁取值为4×10⁻¹¹（mol·cm⁻²·s⁻¹）；高阶溶解速率常数k₄取值为4×10⁻⁸（mol·cm⁻²·s⁻¹）^{[4, 7]}。紊流条件下溶解速率常数比层流条件高一个数量级^{[5, 16]}，即紊流条件下k₁取值为4×10⁻¹⁰（mol·cm⁻²·s⁻¹），k₄取值为4×10⁻⁷（mol·cm⁻²·s⁻¹）。

对于单条裂隙，溶液从左端流入裂隙之前的初始钙离子浓度为C₀（mol·L⁻¹）。由于模拟灰岩裂隙中一维层流运动的质量运输方程仅由岩石与反应溶液的接触时间决定^[3]。而且由于假设中光滑裂隙的两壁与水流之间没有摩擦力，同时钙离子溶解速率远大于钙离子分子扩散速率，在张开度较小的情况下，可认为在横向上钙离子瞬间达到溶解平衡，而在纵向上可忽略分子扩散项，仅考虑对流和溶解反应。

选取裂隙中一段微小体积溶液dV（cm³）作为研究对象，其长度为dx（m），单位宽度为1 cm，裂隙张开度b（cm）。该微小体积从裂隙左端向右移动，在移动的过程中发生动态溶解，溶液与裂隙隙壁一侧的接触面积为dA（cm²），于是该微小体积内钙离子浓度与时间的关系为：

$dV\dfrac{dC}{dt}=2\mathrm{ }\mathrm{d}AF\left(C\right)$

(9)

式中：C为微小体积内溶液的Ca²⁺浓度；F(C)为溶解速率。进一步可得：

$\dfrac{dC}{F\left(C\right)}=\dfrac{2}{b}dt$

(10)

当该微小体积经t时刻到位置x，浓度由C₀变为C，对t积分得C与t的关系。然后由单宽流量q和b得t与x的关系，当q > 0时，t = b×x÷q，可得到浓度C与位置x的关系。进一步可计算裂隙出口的浓度为C(L)，裂隙两端浓度差为ΔC = C(L)− C₀。假设中心节点的钙离子浓度为C_i，其相邻节点的钙离子浓度分别为C_{i, j}。从裂隙节点(i, j)流向节点i的浓度差ΔC_{i, j}为：

$\Delta {C}_{i,j}={C}_{eq}-{C}_{i,j}-\left({C}_{eq}-{C}_{i,j}\right){e}^{-\left(\tfrac{{2k}_{1}}{{{q}_{i,j}C}_{eq}}{L}_{i,j}\right)},\left({C}_{i} \lt 0.9{C}_{eq}\right)$

(11)

$\begin{split} \Delta {C}_{i,j}=&{C}_{eq}-{C}_{i,j}-{C}_{eq}\left({C}_{eq}-{C}_{i,j}\right)\cdot\\ & \sqrt[3]{\dfrac{{q}_{ij}{C}_{eq}}{{\left({C}_{eq}-{C}_{i,j}\right)}^{3}6{k}_{4}{L}_{i,j}+{q}_{i,j}{{C}_{eq}}^{4}}}, \end{split}$

$\left(0.9{C}_{eq} \lt {C}_{i} \lt {C}_{eq}\right)$

(12)

1.2.2 裂隙网络节点浓度计算

在裂隙网络节点钙离子浓度的计算中，只有上游流入的流量和浓度的裂隙段才会对中心节点浓度有影响，而流出的下游裂隙段和相邻节点不影响中心节点浓度的计算。当由相邻节点流入中心节点时，裂隙中流量为负值，设与中心节点相连且流量为负的裂隙数为n，依据钙离子质量守恒，中心节点的钙离子浓度为：

${C}_{i}=\dfrac{\displaystyle\sum \nolimits _{j=1}^{n}{q}_{i,j}\left({C}_{i,j}+\Delta {C}_{i,j}\right)}{\displaystyle\sum \nolimits _{j=1}^{n}{q}_{i,j}},（{q}_{ij} \lt 0）$

(13)

式中：n的取值范围为3、2、1。

进一步整理得：

${C}_{i}×\displaystyle\sum \nolimits _{j=1}^{n}{q}_{i,j}-\displaystyle\sum \nolimits _{j=1}^{n}{C}_{i,j}×{q}_{i,j}-\displaystyle\sum \nolimits _{j=1}^{n}\Delta {C}_{i,j}×{q}_{i,j}=0\mathrm{ },$

(14)

对于所有内部节点i作为中心节点组成浓度方程组，当C >0.9C_eq时，方程组为非线性方程组，同水流方程一样可采用Newton-Raphson方法进行迭代求解。

1.2.3 裂隙网络扩宽

对于裂隙网络内部任一中心节点i连接的一条裂隙，相邻节点(i, j)作为裂隙进口浓度为C_{i, j}，经过裂隙长度L_{i, j}溶蚀到达中心节点i即裂隙出口，浓度增量为ΔC_{i, j}。由此可得裂隙内dt时间的溶蚀量dM，进一步可得单条裂隙的溶蚀速率R_{i, j}：

${R}_{i,j}=\dfrac{dM}{dt}=\Delta {C}_{i,j}×{q}_{i,j}$

(15)

若模拟的时间步长为Δt，经溶蚀后隙宽增量为Δb_{i, j}（cm），那么隙宽增量为：

${\Delta b}_{i,j}=\dfrac{{R}_{i,j}{M}_{{{\rm{CaCO}}}_{3}}}{\rho {L}_{i,j}}\Delta t$

(16)

式中：ρ为碳酸盐岩的密度（2.5 g·cm⁻³）； ${M}_{{{\rm{CaCO}}}_{3}}$ 为碳酸钙的摩尔质量（100 g·mol⁻¹）。

2. 裂隙岩溶发育模拟

本文以二维承压岩溶水系统作为模拟对象，含水层长度为800 m，深度为600 m。左、右两侧为河流定水头边界，水头分别为620 m和600 m。上下两侧为隔水边界。左侧边界河流底部由于水生植物和微生物呼吸作用释放CO₂，其CO₂分压设为P_CO₂ = 0.8%^{[1, 5]}，由开放系统中方解石溶解度随P_CO₂变化规律^[12]对应Ca²⁺的平衡浓度为Ca_eq=2.0×10⁻³ mol·L⁻¹。河流流入封闭含水层的Ca²⁺设为平衡浓度的0.8倍，即流入浓度为1.6×10⁻³ mol·L⁻¹。保持以上边界条件不变，对含水层内部裂隙进行刻画和模拟，保证除裂隙介质本身外其他控制因素都保持不变。

含水层包括两组随机生成不同走向的裂隙，裂隙统计参数如表1。裂隙中心点位置由均匀分布产生，走向服从正态分布，迹长服从对数正态分布，张开度由正态分布随机产生。

表 1. 随机裂隙网络统计参数

Table 1. Statistic parameters of the random fracture network

裂隙组	统计参数	服从分布	均值	标准差	最小值	最大值
	走向	正态分布	30	5	15	45
第一组	迹长/m	对数正态分布	130	10	100	160
	张开度/cm	正态分布	0.005	0.001	0.002	0.008
	走向	正态分布	120	5	105	135
第二组	迹长/m	对数正态分布	130	10	100	160
	张开度/cm	正态分布	0.005	0.001	0.002	0.008

| Show Table

DownLoad: CSV

在裂隙生成过程中去除末端裂隙后，含水层共包含3 200条随机裂隙。在裂隙水流模拟中，迭代计算的水头误差和水均衡误差均满足给定的误差要求。在裂隙溶蚀计算中，根据水流计算收敛情况确定和调整时间步长。本次模拟溶蚀时间步长设为1 000年，模拟时间为500万年，结果如图1，裂隙水流均为层流，没有产生紊流，裂隙张开度均小于0.01 cm。初始裂隙张开度分布图与图1对比均相同，没有发生明显的裂隙扩宽现象。

图 1. 裂隙含水层模拟500万年的水流状态和张开度分布（含水层长宽单位为m，张开度单位为cm）

绿色代表层流，初始时刻模拟图与500万年相同

Figure 1. Flow state and apertures of the fracture aquifer at 5 million year which is similar to the initial state, where the green lines represent laminar flow, the unit of the aquifer length is m, the aperture is shown with the line width, and the unit of aperture is cm

下载: 全尺寸图片幻灯片

尝试将上述岩溶含水层第1组初始裂隙张开度扩大，均值改为0.01 cm，标准差仍为0.001 cm，最小值为0.001 cm，最大值为0.02 cm。图2（a）为初始裂隙网络分布，第一组裂隙的张开度总体上比第二组大1倍。模拟结果显示，在83.9万年出现了紊流，相比原来的含水层岩溶发育速度大大增加。前50万年（图2（b））仅在进口处形成了明显扩宽的裂隙，到70万年（图2（c））含水层内部出现了明显的优势裂隙，但出口处的裂隙和流量仍然一直为微小的变化，直到81.0万年后含水层总流量才超过0.01 cm²·s⁻¹，在这之后出口总流量快速增加，裂隙扩宽也迅速增加，在83.9万年（图2（d））出现紊流，促进了地下河管道的形成。

图 2. 第1组裂隙张开度增大后的含水层水流状态和张开度分布（图中含水层长度单位为m, 张开度单位为cm）

（a）为初始时刻，（b）为50万年，（c）为70万年，（d）为83.9万年出现紊流，红色代表紊流

Figure 2. Flow state and apertures of the aquifer with the first group fracture aperture increased

(a) initial time, (b) 500 thousand year, (c) 700 thousand year, (d) 839 thousand year, where the red lines represent turbulent flow

下载: 全尺寸图片幻灯片

3. 讨论

3.1 裂隙张开度

裂隙网络的非均匀性包括了裂隙张开度的非均匀性，也包括裂隙走向的不均匀性、迹长的不均匀性和分布的不均匀性，后三者构成了裂隙网络结构的非均匀性。为研究张开度的非均匀性对紊流形成的影响，首先要去除裂隙网络结构的影响。将裂隙网络设为均匀的正方形网格状，当张开度相同时即为均质各项同性的含水介质，如图3（a）为均匀裂隙网格。研究共设计了8种不同裂隙张开度统计特征的模拟情形，均值分别为0.005 cm、0.006 cm和0.008 cm，表2为有效对比不同模拟情形，裂隙张开度均在其正态分布99.7%置信水平的区间内随机生成。

图 3. (a)模拟情形A 在500万年的结果与初始裂隙分布相同，(b) 模拟情形C裂隙网络在18.9万年产生了紊流（图中含水层长度单位为m, 张开度单位为cm）

Figure 3. (a) Results of the simulation A at 5 million year which is the same with initial fractures, (b) Tturbulent time 189 thousand year of simulation C

下载: 全尺寸图片幻灯片

表 2. 不同模拟情形的初始裂隙张开度统计参数和紊流出现时间

Table 2. Statistic parameters of the initial aperture and the turbulent time in different simulations

模拟情形	均值/cm	标准差/cm	99.7%置信区间/cm		紊流出现时间/万年
A	0.005	0.001	0.002	0.008	>500
B	0.006	0.000 5	0.004 5	0.007 5	>500
C	0.006	0.001	0.003	0.009	18.9
D	0.006	0.001 5	0.001 5	0.010 5	12.0
E	0.008	−	−	−	42.7
F	0.008	0.000 5	0.006 5	0.009 5	20.3
G	0.008	0.001	0.005	0.011	14.3
H	0.008	0.001 5	0.003 5	0.012 5	9.0

| Show Table

DownLoad: CSV

图3（a）代表模拟情形A的初始裂隙张开度分布，经过500万年的模拟张开度几乎没有变化。图3（b）为模拟情形C在18.9万年出现了紊流，其比模拟情形B的标准差大了一倍，但是模拟B并没有出现紊流，初步说明裂隙张开度的非均匀性更容易引起紊流的产生。从A、B、C三种模拟情形可以看出，在初始阶段裂隙张开度的均值较小的情况下，标准差所代表的非均匀性对紊流的形成起着重要作用，标准差越大说明裂隙张开度变化越大，含水层非均质性就越强，紊流出现的时间也越早。综合对比表2的8种模拟情形，在均值不变的情况下，标准差越大，裂隙网络的非均匀性越强，紊流出现的时间就越早。另外当一组裂隙的张开度明显大于另一组裂隙时，即存在主要裂隙，裂隙网络的非均性会更强。

3.2 裂隙走向

为考虑裂隙走向对紊流形成的影响，首先将裂隙张开度、裂隙中心点位置和裂隙长度都设置为固定值，在模拟过程中仅改变裂隙的走向。以3.1节中模拟情形E为基础，将水平裂隙逆时针旋转45°，模拟结果如图4（a），产生紊流的时间为82.3万年，比模拟情形E的紊流出现时间大了将近1倍。这是因为水平裂隙原本是最短的水力途径，当走向改为斜交后增大了差异性溶蚀寻找的最小水力途径，进在而导致含水层岩溶发育减慢。两次模拟结果表明，当裂隙走向与总体水力梯度方向越接近时，越容易形成紊流。

图 4. （a）模拟情形E水平裂隙逆时针旋转45度的含水层在82.3万年出现紊流；（b）在图2中将第二组裂隙扩大为主要裂隙的含水层在21.7万年出现紊流, 图中含水层长度单位为m, 张开度单位为cm

Figure 4. (a) Turbulent time 823 thousand year of the aquifer based on model E in which the plane fractures were rotated 45° counterclockwise, (b) Turbulent time 217 thousand year with the second group of fractures widen into primary fractures

下载: 全尺寸图片幻灯片

另外对于存在主要裂隙的含水层，主要裂隙的走向决定了紊流出现时间的长短。在第2节中将所述的第一组裂隙扩大为主要裂隙，第二组为次要裂隙，模拟结果如图2中在83.9万年形成紊流。现在我们尝试将第二组裂隙扩大作为主要裂隙，第一组裂隙为次要裂隙，模拟紊流结果如图4（b），紊流出现的时间为21.7万年，比第一组作为主要裂隙时提前了60万年，说明主要裂隙的走向对岩溶发育有很大影响。当主要裂隙方向（120°）与总体水力梯度方向（90°）的夹角小于次要裂隙（30°）与水力梯度方向的夹角时，主要裂隙在溶蚀途径中的比例就会大于次要裂隙，溶蚀途径较平直，溶蚀就会进一步加快。反之如果主要裂隙与总水力梯度方向的夹角大于次要裂隙，溶蚀途径需要经过较多的次要裂隙而变得弯曲而长远，溶蚀速度也会大大减小。因此主要裂隙的方向决定着岩溶发育的快慢，当主要裂隙与水力梯度总方向的角度越小，由主要裂隙构成大部分水力和溶蚀途径，溶蚀形成的优势裂隙连通的次要裂隙也越少，紊流出现的时间就越早。

3.3 裂隙迹长

裂隙迹长的分布影响着裂隙网络的连通性。我们将3.2节随机裂隙含水层中裂隙的迹长平均值减小为100 m，其他裂隙参数保持不变，形成的裂隙网络及模拟结果如图5（a），裂隙连通后有效裂隙段的数量减少至1 882条，相比于原来的3 200条裂隙段，含水层的连通性明显减小，出现紊流的时间为47.1万年，比原来晚了25万年。因此在岩溶含水层中，同一组裂隙和不同组裂隙的迹长如果不能有效的连通形成水流运动，就等同于无效的末端裂隙，岩溶和紊流就不会发育。初始裂隙只有延伸到一定程度能够与其它裂隙连通，才被视为有效裂隙。如果裂隙平均迹长过小，或者裂隙迹长分布差异变化过大，都会导致裂隙连通性较差，影响裂隙水流和溶蚀作用。

图 5. （a）裂隙迹长减小后含水层在47.1万年出现紊流；（b）主要裂隙密度减小后模拟500万年仍没有出现紊流, 图中含水层长度单位为m, 张开度单位为cm

Figure 5. (a) Turbulent time 471 thousand year of the aquifer with the length of fractures decreased, (b) No turbulent flow occurred in 5 million years with the density of primary fracture decreased

下载: 全尺寸图片幻灯片

3.4 裂隙密度

裂隙中心点位置的分布影响着含水层中裂隙密度的变化，将3.2节中随机裂隙含水层中第二组主要裂隙密度减小一倍，第一组裂隙密度不变，其他裂隙特征参数保持不变。生成的裂隙中主要裂隙由240 条减小为120 条，与次要裂隙共形成连通裂隙段1 771条，结果模拟500万年仍然没有出现紊流（图5b），裂隙出口流量没有明显增加，含水层中裂隙张开度与初始时刻几乎相同。这说明裂隙密度尤其是主要裂隙密度对岩溶发育的影响较大。相对于次要裂隙，如果主要裂隙密度偏小，裂隙发育形成紊流的时间就会大大增加，甚至很难形成紊流。

3.5 水力梯度

水力梯度作为外部条件影响着裂隙含水层紊流形成的时间。3.1-3.4节的讨论都是在设定的水力条件下进行的，现在以第二节中的裂隙网络结构为基础，将含水层所有裂隙的张开度都设为一个值，然后通过改变含水层水流边界的水力梯度（0.001~1.0），模拟不同张开度（0.001~0.005 cm）含水层紊流出现的时间。结果如图6，总体上水力梯度越大，紊流出现的时间越早，当水力梯度较小时，模拟500万年仍没有紊流形成，这表示含水层有一定程度的外界水力梯度时才可能形成紊流。另外当裂隙张开度小于0.001 cm时，不管水力梯度有多大，模拟都没有出现紊流。这也说明当张开度小于一定程度时，裂隙岩溶发育很难产生紊流，岩溶几乎不发育。根据相关研究^{[6, 12]}，中大规模岩溶发育所需最小张开度在10⁻³ cm数量级，本次张开度为0.001 cm的模拟结果也与其相吻合。研究中对没出现紊流的情形继续延长模拟时间结果仍然没有紊流出现。

图 6. 不同水力梯度和张开度条件下紊流形成时间的三维柱状图，其中每个柱体上的数字为紊流出现的时间（万年），500万年表示没有出现紊流

Figure 6. 3D histogram of the turbulent time with different hydraulic gradients and apertures

下载: 全尺寸图片幻灯片

上述讨论可以辅助我们推测岩溶含水层的内部结构，根据岩溶水系统的调查资料进行初始裂隙的恢复，通过对含水层初始裂隙在不同裂隙特征条件下和不同水力条件下的岩溶发育模拟，可了解岩溶水系统在过去所处的水力条件和初始内部裂隙的非均一程度，与现今岩溶水系统进行对比分析，在实际岩溶水系统模拟中对裂隙管道内部结构和水流状态的刻画更有针对性。

4. 结　论

在外界流入水中CO₂分压为0.8%、平均水力梯度为0.02的条件下，对不同的初始灰岩裂隙进行渗流和溶蚀模拟，结论如下：裂隙张开度的标准差越大，初始裂隙网络的非均匀性越强，紊流出现的时间就越早；主要裂隙的存在使裂隙网络的非均性会更强，有利于紊流形成；裂隙走向与总水力梯度方向越接近，越容易形成紊流；裂隙平均迹长过小，或者裂隙迹长分布差异过大，都会导致裂隙连通性较差，影响裂隙水流和溶蚀作用；裂隙密度尤其是主要裂隙密度对岩溶发育的影响较大，相对于次要裂隙，如果主要裂隙密度偏小，裂隙发育形成紊流的时间就会大大增加，甚至很难形成紊流。含水层外部水力条件与初始裂隙张开度共同影响着紊流出现的时间，对于给定的含水层初始裂隙，紊流能够形成需要达到最小的水力梯度值才能形成紊流，水力梯度低于该值时含水层将不会有紊流产生，而水力梯度越大紊流形成所需时间也越短；当裂隙张开度小于0.001 cm后，不管如何增大水力梯度仍然没有紊流形成。研究含水层初始裂隙的紊流形成有助于我们了解岩溶水系统内部的裂隙管道结构，提高岩溶水模拟的精度。

致谢：感谢中国地质科学院岩溶地质研究所岩溶水资源团队对本研究的支持！

图 1 裂隙含水层模拟500万年的水流状态和张开度分布（含水层长宽单位为m，张开度单位为cm）

Figure 1.

下载: 全尺寸图片幻灯片

图 2 第1组裂隙张开度增大后的含水层水流状态和张开度分布（图中含水层长度单位为m, 张开度单位为cm）

Figure 2.

下载: 全尺寸图片幻灯片

图 3 (a)模拟情形A 在500万年的结果与初始裂隙分布相同，(b) 模拟情形C裂隙网络在18.9万年产生了紊流（图中含水层长度单位为m, 张开度单位为cm）

Figure 3.

下载: 全尺寸图片幻灯片

图 4 （a）模拟情形E水平裂隙逆时针旋转45度的含水层在82.3万年出现紊流；（b）在图2中将第二组裂隙扩大为主要裂隙的含水层在21.7万年出现紊流, 图中含水层长度单位为m, 张开度单位为cm

Figure 4.

下载: 全尺寸图片幻灯片

图 5 （a）裂隙迹长减小后含水层在47.1万年出现紊流；（b）主要裂隙密度减小后模拟500万年仍没有出现紊流, 图中含水层长度单位为m, 张开度单位为cm

Figure 5.

下载: 全尺寸图片幻灯片

图 6 不同水力梯度和张开度条件下紊流形成时间的三维柱状图，其中每个柱体上的数字为紊流出现的时间（万年），500万年表示没有出现紊流

Figure 6.

下载: 全尺寸图片幻灯片

表 1 随机裂隙网络统计参数

Table 1. Statistic parameters of the random fracture network

裂隙组	统计参数	服从分布	均值	标准差	最小值	最大值
	走向	正态分布	30	5	15	45
第一组	迹长/m	对数正态分布	130	10	100	160
	张开度/cm	正态分布	0.005	0.001	0.002	0.008
	走向	正态分布	120	5	105	135
第二组	迹长/m	对数正态分布	130	10	100	160
	张开度/cm	正态分布	0.005	0.001	0.002	0.008

下载: 导出CSV

表 2 不同模拟情形的初始裂隙张开度统计参数和紊流出现时间

Table 2. Statistic parameters of the initial aperture and the turbulent time in different simulations

模拟情形	均值/cm	标准差/cm	99.7%置信区间/cm		紊流出现时间/万年
A	0.005	0.001	0.002	0.008	>500
B	0.006	0.000 5	0.004 5	0.007 5	>500
C	0.006	0.001	0.003	0.009	18.9
D	0.006	0.001 5	0.001 5	0.010 5	12.0
E	0.008	−	−	−	42.7
F	0.008	0.000 5	0.006 5	0.009 5	20.3
G	0.008	0.001	0.005	0.011	14.3
H	0.008	0.001 5	0.003 5	0.012 5	9.0

下载: 导出CSV

参考文献(16)

[1]	袁道先, 朱德浩, 翁金桃, 朱学稳, 韩行瑞, 汪训一, 蔡桂鸿, 朱远峰, 崔光中, 邓自强. 中国岩溶学[M]. 北京: 地质出版社, 1993 YUAN Daoxian, ZHU Dehao, WENG Jintao, ZHU Xuewen, HAN Xingrui, WANG Xunyi, CAI Guihong, ZHU Yuanfeng, CUI Guangzhong, DENG Ziqiang. Karst of China[M]. Beijing: Geological Publishing House, 1993.
[2]	王大纯, 张人权, 史毅虹, 许绍倬, 于青春, 梁杏. 水文地质学基础[M]. 北京: 地质出版社, 1995 WANG Dachun, ZHANG Renquan, SHI Yihong, XU Shaozhuo, YU Qingchun, LIANG Xing. Fundamentals of hydrogeology[M]. Beijing: Geological Publishing House, 1995.
[3]	Yu Q, Shen J, Wan J, Ohnishi Y. Some investigation on early organization of karst system[J]. Journal of China University of Geosciences, 1999, 10:314-321.
[4]	Dreybrodt W. The role of dissolution kinetics in the development of karst aquifers in limestone: A model simulation of karst evolution[J]. The Journal of Geology, 1990, 98(5):639-655. doi: 10.1086/629431
[5]	Liu Z, Dreybrodt W. Dissolution kinetics of calcium carbonate minerals in H₂O CO₂ solutions in turbulent flow: The role of the diffusion boundary layer and the slow reaction H₂O+ CO₂→ H⁺+ HCO₃⁻. Geochimica et Cosmochimica Acta, 1997, 61(14): 2879-2889.
[6]	Groves C G, Howard A D. Minimum hydrochemical conditions allowing limestone cave development[J]. Water Resources Research, 1994, 30(3):607-615. doi: 10.1029/93WR02945
[7]	Gabrovšek F, Romanov D, Dreybrodt W. Early karstification in a dual-fracture aquifer: The role of exchange flow between prominent fractures and a dense net of fissures[J]. Journal of Hydrology, 2004, 299(1-2):45-66. doi: 10.1016/j.jhydrol.2004.02.005
[8]	Kaufmann G. Modelling karst geomorphology on different time scales[J]. Geomorphology, 2009, 106(1):62-77.
[9]	Reimann T, Rehrl C, Shoemaker W B, Geyer T, Birk S. The significance of turbulent flow representation in single‐continuum models[J]. Water Resources Research, 2011, 47(9):1-15.
[10]	Howard A D, Groves C G. Early development of karst systems: 2. turbulent flow[J]. Water Resources Research, 1995, 31(1):19-26. doi: 10.1029/94WR01964
[11]	Gabrovšek F, Peric B, Kaufmann G. Hydraulics of epiphreatic flow of a karst aquifer[J]. Journal of Hydrology, 2018, 560:56-74. doi: 10.1016/j.jhydrol.2018.03.019
[12]	Dreybrodt W. Principles of early development of karst conduits under natural and man‐made conditions revealed by mathematical analysis of numerical models[J]. Water Resources Research, 1996, 32(9):2923-2935. doi: 10.1029/96WR01332
[13]	于青春, 武雄, 大西有三. 非连续裂隙网络管状渗流模型及其校正[J]. 岩石力学与工程学报, 2006, 25(7):1469-1474. YU Qingchun, WU Xiong, Ohnishi Yuzo. Channel model for fluid flow in discrete fracture network and its modification[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(7):1469-1474.
[14]	王云, 于青春, 薛亮, 马浩. 裂隙岩溶含水系统溢流泉演化过程的数值模拟[J]. 中国岩溶, 2010, 29(4):378-384. doi: 10.3969/j.issn.1001-4810.2010.04.005 WANG Yun, YU Qingchun, XUE Liang, MA Hao. Numerical simulation for the evolution of the overflow spring in fracture-karst aquifer system[J]. Carsologica Sinica, 2010, 29(4):378-384. doi: 10.3969/j.issn.1001-4810.2010.04.005
[15]	高阳, 邱振忠, 于青春. 层流—紊流共存流场中岩溶裂隙网络演化过程的数值模拟方法[J]. 中国岩溶, 2019, 38(6):831-838. GAO Yang, QIU Zhenzhong, YU Qingchun. Numerical simulating method for the karst development of carbonate fracture networks with both laminar and turbulent flow[J]. Carsologica Sinica, 2019, 38(6):831-838.
[16]	刘再华, Dreybrodt W. DBL理论模型及方解石溶解沉积速率预报[J]. 中国岩溶, 1998, 17(1):1-7. LIU Zaihua, DREYBRODT Wolfgang. The DBL model and prediction of calcite dissolution / precipitation rates[J]. Carsologica Sinica, 1998, 17(1):1-7.

施引文献

期刊类型引用(2)

1.	徐子凡，陈喜，刘维翰，刘皓，张志才. 表层岩溶带土岩结构对降雨-径流响应特征的影响. 中国岩溶. 2024(04): 863-875 . 本站查看
2.	刘渝港，贺秋芳，沈立成，范佳鑫. 洞穴溶解有机质组分和循环过程的季节变化特征. 中国岩溶. 2023(03): 456-471 . 本站查看

其他类型引用(0)

资源附件(0)

访问统计

图(6)

表(2)

计量

文章访问数: 1484
PDF下载数: 57
施引文献: 2

0. 引　言
1. 研究方法
1.1 裂隙层流和紊流
1.2 裂隙溶蚀扩宽
2. 裂隙岩溶发育模拟
3. 讨论
3.1 裂隙张开度
3.2 裂隙走向
3.3 裂隙迹长
3.4 裂隙密度
3.5 水力梯度
4. 结　论

初始裂隙对岩溶水紊流形成的影响

1. 中国地质大学(北京)水资源与环境学院, 北京 100083 2. 中国地质科学院岩溶地质研究所/自然资源部、广西岩溶动力学重点实验室, 广西 桂林 541004

作者简介: 焦友军（1990－），男，博士研究生，从事岩溶水资源研究。E-mail：jiaoyj@karst.ac.cn

通讯作者: 于青春（1963－），男，教授，博士研究生导师，从事岩溶水资源研究。E-mail：yuqch@cugb.edu.cn

Influence of initial fractures on the occurrence of karst turbulent flow

1. School of Water Resources and Environment, China University of Geosciences (Beijing) , Beijing 100083, China 2. Institute of Karst Geology, CAGS/ Key Laboratory of Karst Dynamics,MNR&GZAR, Guilin, Guangxi 541004, China

Corresponding author: YU Qingchun, yuqch@cugb.edu.cn

0. 引 言

1. 研究方法

1.1 裂隙层流和紊流

1.2 裂隙溶蚀扩宽

1.2.1 单裂隙溶蚀

1.2.2 裂隙网络节点浓度计算

1.2.3 裂隙网络扩宽

2. 裂隙岩溶发育模拟

3. 讨论

3.1 裂隙张开度

3.2 裂隙走向

3.3 裂隙迹长

3.4 裂隙密度

3.5 水力梯度

4. 结 论

期刊类型引用(2)

其他类型引用(0)

计量

出版历程

目录

0. 引 言

1. 研究方法

1.1 裂隙层流和紊流

1.2 裂隙溶蚀扩宽

1.2.1 单裂隙溶蚀

1.2.2 裂隙网络节点浓度计算

1.2.3 裂隙网络扩宽

2. 裂隙岩溶发育模拟

3. 讨论

3.1 裂隙张开度

3.2 裂隙走向

3.3 裂隙迹长

3.4 裂隙密度

3.5 水力梯度

4. 结 论

1.
中国地质大学(北京)水资源与环境学院, 北京 100083

2.
中国地质科学院岩溶地质研究所/自然资源部、广西岩溶动力学重点实验室, 广西桂林 541004

1.
School of Water Resources and Environment, China University of Geosciences (Beijing) , Beijing 100083, China

2.
Institute of Karst Geology, CAGS/ Key Laboratory of Karst Dynamics,MNR&GZAR, Guilin, Guangxi 541004, China

0. 引　言

4. 结　论

0. 引　言

4. 结　论