利用示踪试验时间-浓度曲线分析岩溶管道结构特征

陈亚洲 ，董维红

（1.吉林大学建设工程学院，吉林 长春 130026；2.吉林大学水资源与环境研究所，吉林 长春 130021）

示踪试验常被用来探寻地下水补给源、污染源、水力滞留时间、岩溶地下河管道连通性及岩溶管道结构特征[1-6]，其中岩溶管道结构特征对水资源合理取用、保护及地下工程安全施工均具有重要影响。

由于示踪试验高效、准确的特点，相关学者依托实际工程进行了大量研究，并总结了大量经验，依据示踪剂浓度沿着地下水流运动方向逐渐降低的特点，根据浓度峰值与距离的关系可以得到含水层渗透系数、渗流速度、平均流速、地下河流量[7-9]。通过示踪试验示踪剂回收率，可确定投放点与出水点之间管道是否单一[10-11]。基于上述普遍认识，美国EPA开发了QTRACER软件用于获取水文地质参数[12-13]。张浪等[14]将示踪试验应用于岩溶区岩溶通道与泉水之间的水力联系研究，利用QTRACER软件获取研究区水文地质参数，并得出落水洞与JS02之间存在水力联系，认为地下水为典型的紊流流态。上述参数的获取为示踪试验地下水数值模拟提供了基础数据[15]。

在岩溶水文地质示踪试验中，利用示踪试验时间-浓度曲线可获取岩溶管道结构特征，时间-浓度曲线出峰个数即岩溶管道的条数，时间-浓度曲线的衰退过程出现钝峰（以下简称“梯度下降”）表征岩溶管道内部有地下湖[16-18]。梅正星[19-20]对喀斯特地区岩溶管道流示踪试验进行研究，结果表明单一管道曲线形态为单一的尖峰，略具对称特点；单一管道有地下湖，曲线形态为下降过程中存在钝锋；对于双管道，梅正星通过工程实例研究示踪试验时间-浓度双峰曲线形态，认为出现双尖峰曲线表明存在2条岩溶管道，无地下湖；先钝后尖，代表主管道存在地下湖；先尖后钝，代表支管道存在地下湖。

杨立铮等[21]通过分析地下水示踪剂运动的基本一维方程中时间-浓度曲线的基本形态，结合工程经验，提出了岩溶通道单管道、单管道有水池模型，多管道、多管道有水池模型。

由于岩溶管道天然发育的不均一性，投放点和接收点之间管道数量、地下湖、地下水承压状态等复杂形态导致示踪试验的时间-浓度曲线形态差异较大，在实际运用中，难以通过曲线叠加、钝锋、不规则上升和下降等形态准确解释多条岩溶管道连接关系、地下湖所处位置、地下水状态。针对上述问题，本文通过地下水溶质运移理论、地下水污染运移数值模拟和对有关参考文献的综合分析，利用示踪试验时间-浓度曲线分析岩溶管道的结构特征。

1 岩溶管道流溶质运移理论

地下水污染运移与地下水流速有关，地下水流速高，溶质运移以对流为主；地下水流速低，溶质运移以分子扩散为主[22]。岩溶地区地下水运移以岩溶管道流为主，径流方式属基岩裂隙流+岩溶管道流[23]。对岩溶地区岩溶管道流建立三维地下水溶质运移对流-弥散方程，见式（1）。

式中：θ——有效孔隙度/%；

C——示踪剂浓度/(mg·L-1)；

t——时间/s；

x——示踪剂运移距离/m；

Dij——弥散系数；

qi——单宽流量/(m2·s-1)；

qs·CS——源汇项；

i、j——三维坐标系任意方向。

式（1）右边有3项：第一项为基岩裂隙流弥散效应引起的溶质运动；第二项为对流引起的溶质运动，包括两种，即考虑渗透系数的基岩裂隙流溶质运动和不考虑渗透系数岩溶管道流溶质运动；第三项为有溶质参与的汇源项。

初始条件：示踪剂投放点浓度为：

边界条件：岩溶管道示踪试验，投放的示踪剂为背景中不存在或存在但浓度极低的溶质，因此其上、下、左、右、顶、底均不发生溶质交互，采用壁面边界，边界上浓度交换为零。

在上述边界条件的限定下，岩溶管道三维溶质运移模型接近于经典的一维对流-弥散方程，即式（3）-（6）[24]。

其中，M为投加示踪剂的量，单位为kg。其解析解为：

根据上式推导过程，可以发现，该解的影响因素为水流速度、弥散系数、运移时间；该方程未充分结合岩溶管道水流过程，在实际管道流过程中，其形态受多种因素控制，岩溶管道横截面、管道阻滞系数等最终导致管道内地下水在移动过程中并不是均匀流态，若管道以有压形式存在，则必然存在中心水流运动快，越靠近管道壁其水流速度逐渐减慢，该过程是现有数值解和解析解存在区别的主要原因。

本文在推导理论穿透曲线的过程中，选择恒定的雷诺系数、弥散系数和管道尺寸，其中管道雷诺系数为5 000，弥散系数为0.5 m2/s，管道长度为0.5 km，管道直径1 m。综合考虑管道阻滞系数及水流粘滞系数，管道内流速设定为变化值，即对管道流速进行切割，管道中心点流速设定为0.1 m/s，周边点流速沿距离递减，到管道壁流速为0.05 m/s。瞬时投加示踪剂20 kg，随着水流运动，溶质扩散呈现图1的三维溶质模型，接收点理论时间-浓度曲线见图2。

图1 岩溶管道流溶质运移三维概念图Fig.1 3D concept map of solute transport in karst pipeline flow

图2 接收点理论时间-浓度曲线图Fig.2 Theoretical time-concentration curve of receiving point

由于本文主要研究内容为岩溶管道流，模拟过程中不考虑基岩裂隙水的补给和储存作用；同时人为添加管道阻滞系数和水流运动黏滞系数，并将其刻画为沿中心至管道壁递减的地下水流速。该结果更符合野外管道特点，因为现实中，管道壁及管道形状不规则，导致管道内不同层位水流速度不同。

2 多岩溶管道连接关系分析

如果投放点到接收点之间存在多条岩溶管道，则岩溶管道存在3种连接关系：汇流、管道之间并联、分支。以2条岩溶管道为例，具体连接关系见图3。

图3 多管道连接方式Fig.3 Multiple pipeline connections

上述 3 种连接方式，接收点对应时间-浓度存在如下关系：（1）汇流岩溶管道，示踪剂接收点时间-浓度曲线最大浓度比预算浓度低；（2）多管道并联，示踪剂接收点的时间-浓度曲线为多峰曲线；（3）分支岩溶管道，接收点示踪剂回收率降低。

对第二种情况，管道并联的多峰曲线，采用式（1）依托Python软件进行数值模拟，模拟双管道不同长度下接收点浓度（管道尺寸恒定），并绘制时间-浓度变化曲线，其中弥散系数、雷诺数取定值，具体参数取值见表1。根据管道径流长度及流速将双峰曲线形态划分成如下3种模型（图4）。

表1 管道连接关系模拟参数Table 1 Simulation parameters of karst conduit connection relationship

图4（a），模型一，时间-浓度曲线中出现 2个孤立的单峰，岩溶主管道地下水径流路径短，溶质运移速度快；岩溶支管道，地下径流路径长，溶质运移速度慢（注:本文主管道指流量相对较大的岩溶管道，支管道指流量相对较小的岩溶管道）。主管道示踪剂浓度恢复至初始水平一段时间后，岩溶支管道示踪剂才运移至接收点，由于运移时间长，在弥散的作用下，后峰浓度最高值下降。

图4 两条岩溶管道并联接收点时间-浓度曲线图Fig.4 Time-concentration curve of two parallel receiving points of the karst pipelines

图4（b），模型二，存在双峰叠加现象，岩溶主管道径流路径短，岩溶支管道相对于图4（a）中径流路径短，在主管道主峰未完全衰退时，支峰开始出现。

图4（c），模型三，存在双峰叠加现象，支流管道率先到达，径流路径短，运移速度快，主峰滞后到达，溶质运移速度相对慢，径流路径长。支峰快速到达，浓度上升快，随后主峰到达，浓度上升速率减慢，最后形成2个连续上升的叠加峰。

3 地下湖分布位置分析

和地下湖位置具体分析。以双管道存在一个地下湖为例（图5），岩溶管道与地下湖位置存在下述 4种关系：（1）地下湖在未分支管道上（可以是分支管道的上部，或分支管道的下部）；（2）地下湖在主管道上；（3）地下湖在支管道上；（4）地下湖在主管道与分支管道上。

图5 地下湖位置与管道关系示意图Fig.5 Diagram of the relationship between the location of karst pool and the pipeline

如果投放点和接收点之间存在地下湖，按岩溶管道数量分析，则地下湖分布位置可划分 2种情况，即单管道存在地下湖和多管道存在地下湖。

对于第一种情况，利用Python软件进行数值模拟，在管道内部设置不同数量的地下湖，地下湖尺寸恒定，为储水200 m3的圆柱体，其余条件不变。单一岩溶管道之间存在地下湖，时间-浓度曲线为单峰曲线，曲线下降过程中存在梯度下降，梯度下降的台阶个数表征管道内地下湖数量。该结论与参考文献[21]结论一致。

对第二种情况，多管道存在地下湖，其梯度下降的台阶个数不完全表征地下湖数量，需结合管道条数

利用Python软件，对图5中的4种情况进行模拟，地下湖尺寸恒定，计算接收点浓度，并绘制时间-浓度变化曲线进行分析。

地下湖位置在未分支的管道上，对应的时间-浓度曲线见图6。图6 (a)中，两个孤峰均存在一个梯度下降，两个下降梯度段示踪剂浓度存在一定的比例关系。图6 (b)中，混合峰下降段均存在一个下降梯度，右侧单峰存在一个下降梯度，两个下降梯度存在一定的比例关系，与图4 (b)相比，混合峰结合部位最低点有所提高。图6 (c)存在两个下降梯度，其中第一个下降梯度为支管道对应的水流，第二个下降梯度为主管道对应的水流，两个下降梯度存在一定的比例关系。

图6 地下湖位置在未分支管道上时间-浓度曲线图Fig.6 Time-concentration curve of karst pool on the unbranched pipeline

地下湖位置在主管道上，对应的时间-浓度曲线见 图7。图7 (a)中，左侧主管道单峰存在一个下降梯度，导致主管道下降速度缓慢。图7(b)中，混合峰中间主管道下降曲线存在一个梯度下降，与图4(b)相比，混合峰最低点有所提高。图7 (c)中，右侧下降曲线中存在一个下降梯度。

图7 地下湖位置在主管道上时间-浓度曲线图Fig.7 Time-concentration curve of solution pool location on the main pipeline

地下湖位置在支管道上，对应的时间-浓度曲线见图8。图8 (a)中，右侧支管道单峰存在一个梯度下降，导致支管道下降速度缓慢。图8 (b)中，混合峰右侧支管道单峰存在一个梯度下降，与图7 (b)相比，混合峰最低点降低；图8 (c)中，右侧下降曲线中存在一个梯度下降，与图7 (c)相比，右侧下降曲线中梯度下降时间要超前。

图8 地下湖位置在支管道上时间-浓度曲线图Fig.8 Time-concentration curve of solution pool location on branch pipeline

地下湖位置在主管道和支管道上均存在时，时间-浓度曲线图与图6地下湖未分支管道上形态类似，不同点在于由于主管道及支管道地下湖体积不同，主管道及支管道下降梯度示踪剂浓度不成比例。

4 岩溶管道流状态分析

岩溶管道内部地下水状态存在多种形式，地下水在岩溶管道内部径流时，多数情况下会发生多次表流、承压流之间的相互转化，如图9（a）所示。在潜水、承压水转化过程中，管道过流断面会发生相应变化，中间发生三段变径，如图9（b）所示。利用Python软件进行数值模拟，其中物理模型为单一管道，不同区段截面积不同，弥散系数不同（根据文献[23]，弥散系数与流速存在正相关，本文未做模拟试验，弥散系数取经验值），参数取值见表2。模拟结果见图10。

图9 岩溶管道水流状态变化示意图Fig.9 Schematic map of water flow state change in karst pipeline

图10 岩溶管道水流状态变化时间-浓度曲线图Fig.10 Time-concentration curve of flow state change in the karst pipeline

表2 岩溶管道流状态分析模拟参数表Table 2 Simulation parameters of karst conduit flow state analysis

图10中，当岩溶管道水流状态发生变化时，曲线上升段存在明显的形态变化。岩溶水以潜水状态存在时，流速相对缓慢，弥散系数小，曲线上升幅度较慢；岩溶水以承压水状态存在时，曲线上升幅度较快，即截面积不同，在流量恒定的条件下，潜水状态和承压水状态对应的流速不同，弥散系数不同，会影响示踪剂的浓度变化及接收端的接收时间。根据图中浓度变化过程及对应时间可分别获得潜水段及承压段直线距离及水流速度。由于岩溶地下水状态变化，导致地下水流速和弥散系数发生变化，示踪试验时间-浓度曲线极速下降段同极速上升段在数量和形式上存在相关性。但该极速下降不同于前节地下湖梯度下降，地下湖内由于示踪剂浓度以弥散形式扩散，对应时间-浓度曲线呈现缓慢下降；管道内地下水形态变化，导致示踪剂不均匀下降，且下降速度较快。

5 结论

（1）根据曲线出锋个数，可判断岩溶管道条数。对多岩溶管道，根据多峰曲线形态可判断多管道中的管道长度及水流速度。

（2）根据曲线下降梯度个数可判断地下湖个数。进一步的，若2条及以上岩溶管道中存在地下湖，则在分析野外示踪试验曲线时，应分析多曲线下降梯度之间的相关性，若多曲线下降梯度之间存在相关性，则地下湖存在于未分支的管道上，若其不存在相关性，则其分别存在于多条岩溶管道上。若仅其中1条单峰存在下降梯度，则地下湖仅存在于单峰所在的岩溶管道上。

（3）岩溶管道地下水状态变化，管道承压段对应的时间-浓度曲线，在上升段和下降段均存在快速变化；管道表流段对应的时间-浓度曲线，在上升段和下降段变化幅度相对较缓慢；不同水流状态的下降段同上升段在数量和形式上存在相关性。