岩溶关键带年内水文动态变化规律

毛龙富， 刘 宏， 周恩民， 易 琦， 党翠翠， 杨丽瑞， 付 舒， 张净净

(1.云南大学 地球科学学院， 云南 昆明 650500； 2.云南大学 国际喀斯特联合研究中心， 云南 昆明 650223)

1 研究背景

岩溶关键带位于大气-水-岩石-土壤-生物等圈层交汇地带，受五大圈层共同的物理、化学、生物及其耦合作用[1]，垂向包括植被冠层的顶端至含水层的底部[2]，岩溶管道-洞穴-地下河-隔水层组成了巨大岩溶地下空间[3]。岩溶水是岩溶关键带物质、信息、能量等交流的纽带，也是岩溶地区的一种重要水资源[4]，世界人口的至少四分之一依赖岩溶水资源[5]，其水文功能对岩溶区生态水文的进一步研究具有深远意义[6]。裸露碳酸盐占据全球陆地表面积的15%，世界上最大的岩溶区位于中国，约达363.10×104km2，其中滇、贵、川、渝、桂、湘、鄂、粤等南方7 省1市碳酸盐岩出露总面积为114.26×104km2，占土地面积的27.36%[7-8]。由于特殊的岩溶地质环境、显著的季风气候、脆弱的生态环境，中国南方岩溶地区水资源具有脆弱性和敏感性，同时存在岩溶干旱、内涝、地下水易污染难治理等问题。受全球气候变化和人类活动等共同影响，岩溶地区水资源供需矛盾日益突出，并导致水土流失、石漠化、贫困等恶性循环。进一步解译岩溶水文过程是改善民生、保障岩溶地区生态优良和经济高质量发展的关键。

目前，国内外学者主要通过岩溶水动力学、水文地球化学、同位素水文学、数值模拟、灰箱模型、地球物理等方法来探究岩溶关键带水文过程。岩溶关键带是同时存在溶蚀裂隙、溶蚀孔、管道、溶井和溶洞等多重水流的系统[9]，岩溶含水层的非均匀性和各向异性给岩溶水文过程的研究带来了挑战。Mangiarotti等[10]运用混沌模型分析了两个不同水文特征、不同地质气候和人为作用条件下的岩溶泉，进一步认识到岩溶系统水文动态的低可预测性和非线性，揭示了泉水与气候条件之间的非线性关系。Jourde等[11]、Zhang等[12]发现非饱和带-饱和带子系统的岩溶水动力具有垂直变化规律，表层岩溶泉水流量对降雨可产生短中期、长期响应。于苗等[13]基于时间序列分形理论对岩溶大泉的动态变化特征进行了研究，并探讨了大气降水、人工开采和人工补给量对其的影响关系。δ18O是联系降水、洞穴滴水的重要线索[14]，通过氢氧同位素可预估每月降水滞留在岩溶包气带上的水量[15]。曾成等[16]采用水文水化学自动监测仪器对贵州南部的板寨岩溶地下河系统进行了野外监测，发现表层岩溶水文地球化学具有年内动态变化规律。

在岩溶地区，岩溶含水层复杂且特殊，给岩溶水文过程的研究增加了难度，洞穴滴水的存在为深入探究岩溶泉水的水文运移过程提供了便利条件。黄龙泉域地处典型的岩溶区，岩溶地貌发育成熟，水文地质构造复杂，老黄龙洞洞穴滴水-泉水近乎垂直分布，是进一步解析岩溶关键带水文动态变化规律的天然实验场。

对黄龙泉域老黄龙洞的降水—洞穴滴水-泉水(即降水—地下水—地表水)的水文差异特征进行综合研究，旨在剖析典型岩溶关键带水文动态变化规律及其原因，加深认识岩溶地下水运移模式和水循环机制，为广大岩溶地区水资源调控、管理、防旱抗涝等工作提供科学参考。

2 数据采集与研究方法

2.1 研究区概况

黄龙泉位于云南省昆明市嵩明县阿子营乡朵格村，是典型的暗河式岩溶泉，其泉口(25°25′49.552′′N，102°55′55.016′′E)海拔高程为2 296.70 m。老黄龙洞口海拔高程为2 331.64 m，处在垂直高度比黄龙泉口高约35 m的半山腰，洞穴内发育有落水洞、洞穴滴水和石笋、石柱、石花、石坝等岩溶景观。研究区老黄龙洞-黄龙泉区位概况见图1。

黄龙泉域在昆明市水源保护区内，属于滇池流域北部分水岭地带，是昆明市重要饮用水源松华坝水库的补给河流——牧羊河的主要源头之一(图1(a))，每年向水库注入700×104m3的优质水，也是当地村民的饮用水源。由黄龙泉域的水文地质状况(图1(b))来看，其补给区以碳酸盐岩地层为主，由上司段威宁组(C1d2)、石炭系下统大塘组(D2h1-2)和上泥盆统宰格群地层(D3zg)泥盆系华宁组海口段(O1t)构成，主要含水层为泥盆系、石炭系碳酸盐岩。洞穴主要发育在石炭系的下统大塘组上司段地层(C1d2)，上部主要为白云岩与灰岩，下部主要为灰岩质角砾岩。洞穴深部即尾端以威宁组地层为主，主要发育有灰白色厚层、浅灰和块状纯灰岩。洞穴局部岩层发育有角砾状灰岩及薄层砂页岩[17]。黄龙泉主要由降水入渗补给，是典型的岩溶包气带补给型岩溶泉(图1(c))，泉水平均流量为106 L/s，矿化度为137 mg/L。

图1 研究区老黄龙洞-黄龙泉区位概况

该区域同时受到东亚季风和印度季风的影响，属亚热带高原季风气候区，多年平均气温为15.07 ℃，多年平均降水量为962.94 mm，干湿季交替显著，雨季(5-10月)的降水量占全年的86.8%；干季(11-次年4月)的降水量则低于各月平均降水量，仅占年降水量的13.2%。该区碳酸盐岩裸露，以黄棕壤为主，土层发育厚度不一，由山顶向周边山麓递减，山顶土层厚度约100 cm，山麓土层厚度在0～30 cm之间。植被主要以云南松林、高山针叶林及灌丛草甸等耐旱植物交错分布，山顶处种有大量人工松林。

2.2 数据采集与监测

在黄龙泉口安置HOBO U24-001水文监测计，每30 min自动监测保存电导、水位、水温等数据。依据流量和水化学变化特征，可分为季节性滴水(A)、溢流性即敏感性滴水(B)、常年性即稳定性滴水(C靠近洞口、D在洞穴深部)[18-19](图1(d))。在老黄龙洞里选取上述4个不同类型的典型滴水点，借助安全帽灯和高精度秒表现场监测读取洞穴滴水的滴率，为减小或避免误差，其中滴率较快的滴水点取25～35滴的均值，滴率较慢的滴水点取3～5滴的均值，滴率中等的滴水点一般取10～15滴的均值，每滴水流量取0.05 mL。距泉口约100 m外安置雨量计记录降水、气温。选取的数据采集日期为2015-11-30-2017-10-05和2019-01-01-2020-04-17，洞穴滴水的采集日期为2016-03-06-2017-10-05，逐周监测和采集样品。1951-2020年的气象数据源自中国国家气象局，其中蒸发量指蒸发皿大型蒸发量。

2.3 分析方法

(1)自相关性分析。自相关(autocorrelation function)由Larocque于1998年提出[20]，指的是时间序列相关性，量化一段时间内连续序列值x=(x1,x2，…，xt，…，xn)的线性依赖关系。自相关系数r(k)可定义为：

(1)

(2)

(k=0,1,2,…,m)

从水文地质角度分析，r(k)揭示了岩溶含水层状况、储水结构以及水文过程，记忆时间长短说明了t与t+k系列之间的关联程度，反映出岩溶系统及其变化程度[24]。

(2)互相关性分析。互相关(cross correlation function)[25]指两个变量的时域上输入信号x=(x1,x2，…，xt，…，xn)与输出信号y=(y1,y2，…，yt,…，yn)的相关关系，当k≥0时，互相关系数rxy(k)可定义如下[26]：

(3)

(4)

互相关图最大振幅rxy(k)出现时，所对应的滞后时间(k)最能反映水文脉冲传输至含水层的延滞时间，当rxy(k)≤0.2时表示两组时间序列互相关性消失或不明显。根据显著性界限确定自相关性、互相关性的显著程度。

(3)Pearson相关系数。通过计算Pearson相关系数(rxy)分析探讨各要素之间的复杂关系[27]。

(5)

(4)年内分配不均匀系数与完全调节系数。降水、水位、滴率等要素的年内分配不均匀系数Cu与完全调节系数Cr定义如下[28]：

(6)

(7)

(8)

(9)

(10)

(5)集中度与集中期。集中度(Cd)、集中期(Dm)是将一年内泉水、滴水等水文气象序列{R(i)，i= 1，2，…，12；i为年内某个月}作为向量分析，可定义为[32]：

(11)

(12)

Dm=arctan(Ry/Rx)

(13)

式中：Cd为集中度；Rx、Ry分别为R在x、y方向上的分解矢量；θ(i)为第i月的矢量角度，θ(i) = 2πi/12；Dm为集中期[33-34]。

3 结果与分析

3.1 水文要素动态变化特征

根据相关气象资料计算1951-2020年昆明市嵩明县降水距平百分数，结果如图2所示。昆明市嵩明县1951-2020年年均降水量为960 mm，其中雨季降水量占86.8%，干季降水量占13.2%，近70年该区域降水量呈减少趋势(平均变化率为-1.328 mm/a)(图2)、气温呈上升趋势(平均变化率为0.016 ℃/a)，总体呈趋暖干化，旱涝率占45.7%。该区域季节性特征显著，各月间降水量差距大。2016年属于平水年，降水量为870.1 mm，雨季降水量为728.2 mm，占年降水量的83.7%，其中7-9月降水量为477.1 mm，占年降水量的54.8%；干季降水量明显比雨季少，仅占年降水量的16.3%。

图2 1951-2020年昆明市嵩明县降水距平百分数

根据数据采集期记录资料得出研究区大气降水-洞穴滴水-泉水水文要素动态变化特征，如图3所示。

图3 数据采集期研究区大气降水-洞穴滴水-泉水水文要素动态变化特征

洞穴滴水数据采集期老黄龙洞不同类型(A、B、C、D)洞穴滴水的滴率和电导率年内变化特征以及电导率变化曲线分别见图4、5。

图4 数据采集期老黄龙洞不同类型洞穴滴水年内滴率及电导率变化特征

图5 数据采集期老黄龙洞不同类型洞穴滴水电导率变化曲线

分析图3中滴水类型A、B、C、D的滴率变化曲线可知，不同类型的洞穴滴水产生了不同的水文响应，季节性滴水(A)的年均滴率为176.47 mL/h，有显著的季节性特征，雨季滴水流量远大于干季，滴率水文动态变化曲线滞后于降水量峰值出现的时间，呈多峰曲线响应，峰值出现在雨季中、后期且时间较短，而干季降水量较少，未产生显著响应。敏感性滴水(B)的年均滴率为82.56 mL/h，其滴率变化曲线对降水的响应比较迅速，滴率峰值与降水量峰值相近，对降水事件均有明显响应。稳定性滴水(C、D)的年均滴率分别为2.71、1.08 mL/h，年内滴水滴率变化曲线变化较小，对降水量峰值响应微弱。

分析图4、5可知，老黄龙洞各类型洞穴滴水滴率及其变化程度由大到小排序为A>B>C>D(图4)，洞穴滴水A、B、C、D的平均电导率分别为592.61、620.17、648.10、1 560.53 μS/cm(图5)，总体上电导率大小排序为A

表1列出了老黄龙洞各类型洞穴滴水滴率与电导率之间的相关系数。由表1可看出，洞穴滴水的电导率与滴水滴率呈显著正相关，其相关程度排序为B>A>C>D，表明滴率大的洞穴滴水对降水响应明显。

表1 各类型洞穴滴水滴率与电导率之间的相关系数

黄龙泉水文气象各参数的相关系数见表2。由表2可见，黄龙泉水位、水温、电导率均对降水产生了明显的动态响应，其中泉水水位与电导率的相关系数达到0.608，呈显著正相关，两者变化具有良好的一致性。泉水水位的峰值与降水量的变化呈相似趋势，但峰值之间没有呈递减趋势，而电导率的峰值出现在雨季初期且向雨季末期递减，表明干季高度矿化的岩溶裂隙水被淋滤挤出，产生“活塞效应”，促使洞穴滴水和岩溶泉水的电导率偏高。雨季包气带水受到新进雨水的混合作用，洞穴滴水和岩溶泉水的电导率偏低，降水补给后表层岩溶带以弥散流形式渗透补给[35]。泉水水温与水位、电导率均呈显著负相关，相关性由冬季向夏季递增，至雨季时突然急剧减小，主要原因是受管道水快速补给，泉水中高度矿化的岩溶水比例减小。水文气象各参数的相关性分析表明，黄龙泉水受土壤、表层岩溶带、前期裂隙管道储存水、饱水带水流的共同补给。

表2 黄龙泉水文气象各参数的相关系数

根据相关气象资料，2019年黄龙泉域降水量为913.1 mm，最大降水量出现在7月，泉水水位最高值在9月；2016年降水量最大月为8月，泉水水位最高值在9月，不同年份的降水量差异使泉水水位具有年际、年内动态变化特征，表明雨水是黄龙泉的主要补给来源。计算得出的2016与2019年黄龙泉水文气象各参数的相关系数，见表3。

表3 2016与2019年黄龙泉水文气象各参数的相关系数

由表3可见，2016与2019年的降水量呈弱相关，而该两个年份的泉水水位、水温、电导率和气温均呈显著正相关，泉水水位过程线比降水稳定。由此表明黄龙泉水受基岩裂隙流和管道快速流的共同补给，泉水流量比较稳定，表层岩溶系统具有一定的调蓄功能。

3.2 年内Cu、Cr、Cd、Dm等系数特征

在季风气候区，降水和气温等要素具有明显的季节性变化特征，从而导致洞穴滴水、泉水等岩溶水呈季节性水文波动。以2016年为例，计算老黄龙洞-黄龙泉降水、洞穴滴水及泉水水位等参数的年内分配不均匀系数Cu、完全调节系数Cr、集中度Cd和集中期Dm，还计算了2019年降水和泉水水位的集中期Dm作为对比，结果见表4。

表4 老黄龙洞-黄龙泉降水、洞穴滴水及泉水水位的年内Cu、Cr、Cd和Dm值

由表4可知，2016年各类型洞穴滴水年内分配不均匀系数Cu最大为1.02、最小为0.10，平均值为0.42，Cu值排序为A>C>B>D，降水和泉水水位的Cu值分别为0.91和0.21。2016年各类型洞穴滴水年内分配完全调节系数Cr的最大值为0.47、最小为0.04，平均值为0.19，Cr值排序为A>C>B>D，降水和泉水水位的Cr值分别为0.38和0.09。由此可见，在岩溶关键带水文动态变化过程中，降水-洞穴滴水-泉水水文过程的Cu、Cr值均呈垂直递减趋势。2016年各类型洞穴滴水的集中度Cd也有显著差异，季节性滴水A的集中度最大，为0.66，远大于敏感性滴水B和稳定性滴水C、D的集中度，且大于降水和泉水水位的集中度。由于敏感性和稳定性滴水年内变化较小，故其集中度也较小。黄龙泉域常年降水的集中期Dm在6月初左右，2016年降水的集中期Dm在7月下旬至8月上旬之间，而洞穴滴水、泉水水位的集中期在8月中下旬，表明降水入渗形成洞穴滴水的过程中受上覆土壤厚度、表层岩溶带、包气带、滴水和补给管道运移方式等因素的影响[36]。

3.3 时滞分析

老黄龙洞各类型洞穴滴水及黄龙泉水位的自相关系数r(k)变化曲线见图6。

图6 老黄龙洞各类型洞穴滴水及黄龙泉水位的自相关系数r(k)变化曲线

由图6可看出， 作为泉域补给区的洞穴滴水在第1～2周内的自相关系数r(k)快速减小，之后随时间变化趋于缓和，这种自相关系数曲线变化解译了岩溶含水层的三维二元性。稳定性滴水(C)自相关系数减小缓慢，在第8～9周期间自相关系数下降至0.2，故滞后时间可达56～63 d；稳定性滴水(D)自相关系数曲线变化复杂，在第1、5～6、10周自相关系数值均大于0.2，即滞留时间分别为14、35～42、70 d，降水对该类滴水的影响持续时间为14～70 d，该类滴水对降水产生多重响应，滴水由历经土壤、表层岩溶带和管道、裂隙的雨水入渗补给。敏感性滴水(B)在开始时自相关系数急剧减小，对降水响应迅速，在第1～1.5周时，自相关系数降至0.2，降水滞留时间约为7～10.5 d。季节性滴水(A)的自相关系数曲线变化相对和缓，介于稳定性滴水C、D的自相关系数曲线之间，在第7～8周自相关系数降至0.2，表明滞后时间为49～56 d。黄龙泉水位初期的自相关系数变化也相对和缓，介于稳定性滴水C、D的自相关系数曲线之间，与季节性滴水(A)的自相关系数变化趋势相似，反映了泉水既受裂隙基质流补给也受管道流补给。在第7周时，水位自相关系数降至0.2，表明黄龙泉水位滞后时间为49 d。由于岩溶含水层透水介质的不均匀性，加之岩溶管道、裂隙、洞穴发育复杂程度的差异[37]以及地下水运移路径的差异，研究区域植被覆盖良好致使降水对洞穴滴水、泉水的持续影响具有分带运移性。

3.4 响应分析

将降水作为岩溶系统的输入，各类型洞穴滴水、泉水水位作为输出序列，计算降水与各类型洞穴滴水、泉水水位的互相关系数rxy(k)，其互相关系数rxy(k)变化曲线如图7所示。

图7 降水与老黄龙洞各类型洞穴滴水及黄龙泉水位的互相关系数rxy (k)变化曲线

由图7(a)可知，季节性滴水(A)的rxy(k)值在第1～11、14周均大于0.2，在第6周时rxy(k)出现峰值0.603，表明该滴水点对降水响应的滞后时间为42 d，且存在7～77 d的最明显降水响应。敏感性滴水(B)的rxy(k) 值在第1周就出现了峰值0.330，表明该滴水点对降水响应迅速且滞后时间为7 d，其他时段rxy(k)均未出现大于0.2的情况，说明其他时间未对降水产生显著响应。稳定性滴水(C、D)的互相关系数变化曲线走势大致相同，表明其对降水的响应基本一致，但滴水D的变化幅度小于滴水C。滴水C的rxy(k)在第30周时出现峰值0.220，说明该点滴水对降水响应的滞后时间为210 d，其他时期rxy(k)均小于0.2，表明响应不显著。滴水D的rxy(k)值在第0～45周均小于0.2，表明滴水D在该水文年无显著降水响应。上述分析表明，各类型洞穴滴水对降水的响应程度排序为季节性滴水>敏感性滴水>稳定性滴水。

由图7(b)可知，黄龙泉水位的rxy(k)值在第1～51、58～61 d均大于0.2，说明泉水水位对降水的响应滞后时间平均约为41 d，其中出现了4个依次递减的峰值(0.370、0.363、0.358和0.220)，分别为4、9、22和59 d，第4、9 d可能为竖井、落水洞等点状快速流补给泉水；第22、59 d可能受线状裂隙流、网状渗透流等中慢速弥散流补给泉水，各补给通道及洞穴滴水类型对泉水的补给影响程度排序分别为点状>线状>面状、季节性滴水≥敏感性滴水>稳定性滴水，体现了岩溶介质空隙水流和储存的双重性，表明岩溶含水介质的孔隙度、渗透率和导水率空间分布极不均匀[38]，泉水受裂隙介质控制的慢速弥散裂隙流与集中补给的管道流的共同补给，其中表层岩溶带对岩溶水起到了入渗补给和调蓄的双重作用，延长了岩溶系统对雨水入渗的滞时。

4 结 论

(1)1951-2020年昆明市嵩明县呈暖干化趋势，黄龙泉域洞穴滴水、泉水水文动态变化具有季节性和垂直分异性，黄龙泉水位动态变化特征与季节性滴水最为相似。年内分配不均匀系数Cu、年内分配完全调节系数Cr和集中度Cd在降水-洞穴滴水-泉水过程中总体呈垂直递减趋势，其中敏感性滴水≤稳定性滴水<泉水水位<季节性滴水。黄龙泉域常年降水集中期Dm在6月初左右，2016年洞穴滴水、泉水水位的集中期在8月中下旬。

(2)黄龙泉域洞穴滴水和泉水水位均对降水产生滞后响应，其滞后时间排序为敏感性滴水<泉水水位≤季节性滴水<稳定性滴水；对降水的响应程度由弱到强排序为稳定性滴水、敏感性滴水、季节性滴水。稳定性滴水的补给为弥散渗透流，对降水响应的滞后时间约为210 d；季节性滴水的补给为管道流，对降水响应的滞后时间为42 d；敏感性滴水为包气带底部补给，对降水响应迅速，滞后时间约为7 d。泉水水位响应滞后时间与季节性滴水最为接近，约为41 d。泉水水位对洞穴滴水类型的响应程度排序为稳定性<敏感性<季节性。

(3)黄龙泉域岩溶含水层中裂隙流与管道流并存、层流与紊流并行，泉水由流经岩溶区的风化裂隙水和降雨入渗补给，补给方式包括点状集中、线状慢速、面状入渗等，加之土壤水、表层岩溶带水、包气带水的调蓄作用，形成常流性暗河泉。

本研究进一步揭示了岩溶关键带地下水水文动态变化规律，以期为岩溶地区水资源管理、调控提供一定的科学参考。