考虑非饱和带作用及管道流的岩溶泉流量模拟

2019-10-25覃夏南姜光辉

桂林理工大学学报 2019年3期

覃夏南,姜光辉,夏源

(1.桂林理工大学 a.岩溶地区水污染控制与用水安全保障协同创新中心; b.广西环境污染控制理论与技术重点实验室, 广西桂林 541006; 2.中国地质科学院岩溶地质研究所/自然资源部、广西岩溶动力学重点实验室, 广西桂林 541004)

我国西南岩溶地区普遍存在干旱缺水、洪水内涝等资源环境问题,开展岩溶区地下水水质与水量数值模拟对岩溶水资源评价有重要意义。在表层岩溶带发育的地区, 非饱和带的模拟不可缺少, Fox等[1]利用传统的MODFLOW河流模块对其进行调整, 模拟河流与地下水系统的饱和以及非饱和带流动; 美国地质调查局开发联合降雨-径流系统(PRMS)和以MODFLOW-2005的地下水和地表水相互作用的GSFLOW; 凌敏华等[2]运用UZF模块与MODFLOW模拟地表水过程与地下水动力过程耦合; 程勤波等[3]探讨了饱和与非饱和带土壤水动力耦合的模拟及其入渗试验。岩溶区饱水带普遍分布着管道流[4], 管道水流模拟也是近年来岩溶水流研究的热点。对此，成建梅等[5]以广西环江北山岩溶含水系统为例, 建立岩溶管道-裂隙-孔隙-三重介质的地下水模型; 陈崇希等[6]阐述入渗-管道流耦合模型的原理和应用方法; 孙晨等[7]利用裂隙-管道介质物理模型模拟了岩溶管道裂隙的泉流量衰减过程; 常勇等[8]用水箱-CFP模型模拟丫吉试验场泉流量衰减; 赵良杰等[9]基于MODFLOW用Drain与River模块对岩溶管道水流的模拟进行对比讨论；袁道先等[10]以丫吉试验场为例，探讨岩溶峰丛山区岩溶水系统及其数学模型;章程等[11]运用SWMM模型模拟岩溶峰丛洼地系统降雨径流过程;夏源等[12]在岩溶流域用分数阶模型模拟泉流量。

目前,很少看到同时考虑岩溶区饱和-非饱和带的数值模拟。因此，本文尝试利用包含UZF模块以及CFP管道流模块的MODFLOW软件建立一个岩溶泉流域的数学模型并模拟岩溶泉流量,研究在岩溶峰丛洼地非饱和带内表层岩溶带及管道的水文过程及其水文功能。

1 研究区概况

丫吉试验场总面积约2 km2,位于桂林市东郊约8 km处峰丛洼地与峰林平原交接地带的丫吉村附近(图1)。试验场为峰丛洼地地貌,内部含有多个大小不等、高低不同的洼地,洼地底部标高最低250 m,最高400 m,最高峰625 m,与周围标高只有150 m左右的平原地区差别较大。

图1 桂林丫吉试验场水文地质图(据文献[11])Fig.1 Hydrogeological map of Yaji experiment site in Guilin

桂林市处于低纬度,属于亚热带季风气候,雨量充足,四季分明,并且高雨量与高热度基本同季,降雨在季节上分布明显不均衡。根据桂林市北部气象局气象站多年观测结果显示,多年平均降雨量为1 949 mm,年平均蒸发量为1 490～1 905 mm,降雨集中在每年4—8月,大约占全年的65%。

试验场出露的地层为上泥盆统融县组(D3r)灰岩(浅灰—灰白色致密质状中厚层泥亮晶颗粒灰岩),灰岩纯度较高,总厚度达387 m。

本文研究丫吉试验场流域面积最大的S31号泉流域,S31泉是一个较为独立的水文系统,主要受降雨补给，而无其他外源水流入, 流域面积约1.0 km2, 水文地质剖面图见图2。根据以往的示踪试验[10], 该泉主要受场地内1、 3、 4号洼地的补给, 其中3号与4号洼地高于1号洼地40 m左右。流域内表层岩溶带发育,分布5个表层岩溶泉。泉域内3个洼地底部均存在若干个落水洞,所以可以推测S31号泉域底部存在一条主要由管道、溶蚀裂隙与溶蚀空洞组成的通道连接泉域内各个洼地底部的落水洞，并排泄洼地内地下水。该泉对降雨响应灵敏,泉流量动态变化较大,枯季时一般为每秒数升,丰季暴雨时可达每秒数立方米,变化幅度很大。

2 水文地质概念模型和数学模型

2.1 研究区范围和边界条件概化

根据试验场水文地质条件及地下水的主要流向,本次研究范围为丫吉试验场S31号泉流域，包含1、 3、 4号洼地(1号洼地面积0.39 km2, 3号洼地面积0.27 km2, 4号洼地面积0.32 km2),泉域面积约1 km2。

图2 桂林丫吉试验场S31号泉岩溶水文地质剖面图Fig.2 Karst hydrogeological profile map of Well No.S31 in Yaji experiment site of Guilin1—土壤覆盖;2—表层的岩溶泉;3—饱水带泉;4—石灰岩;5—洼地编号

侧向边界: 因S31泉域由3个洼地组成，且试验场无其他外源水流入,S31号泉主要受1、 3、 4号洼地补给; S32号泉则主要受5号洼地远距离补给; S291号泉受到2、10号洼地补给; S29则主要受到11号洼地的补给, 所以按照与系统排泄端发生主要水力联系的洼地以及地形等高线图来划分地表以及地下分水岭,以此来圈定泉域的边界,并设置为隔水边界。

垂向边界:潜水含水层的潜水面为系统上边界,通过此边界,潜水与系统外产生垂直方向的水量交换,如大气降水入渗补给、蒸发排泄等,整个饱水带在垂向上剖分为一层,场地内的钻探钻孔岩心资料表明,在泉口附近65 m以下基本没有溶蚀裂隙,深部潜流带虽然存在,但排泄量与泉水相比可以忽略。边界上的侧向径流也是存在的,但是对于侧向径流量并未开展过多的研究,示踪试验的结果表明，泉域内绝大部分地下水还是通过管道向S31号泉排泄的,所以把泉域边界作为模型的隔水边界。由于模拟过程中,饱和带水位变化远小于地表标高,因此模型顶部采用场地地表的实际高程,并根据地形图的等高线对模型进行插值计算出地表标高。

2.2 数学模型

2.2.1 地下水数值模型根据上述水文地质模型,将S31泉域岩溶水系统概化为非均质、各向异性的三维非稳定地下水流模型,按照《地下水动力学》建立数学模型[13]:

(1)

式中:Ω—渗透区域;h—含水层的水位标高, m;Kxx、Kyy、Kzz—x、y、z方向上的渗透系数, m·d-1;Kn—边界面的法向方向上的渗透系数, m·d-1;μ—含水层的贮水率, m-1;ε—含水层的源汇项, d-1;h0—含水层的初始水头分布, m;Γ1—含水层渗流区域的第二类零流量边界;Γ2—含水层渗流区域的第二类已知流量边界;q—Γ2边界上的流量,m3·d-1。

模型空间剖分为30行×60列, 单元格为35 m×35 m, 有效单元格797个, 根据洼地形态及土壤覆盖分布分为两个区域：一区为S31号泉所在的西坡到1号洼地, 二区为3、4号洼地。模型模拟识别期为10天, 即从2009年7月22日0时—8月1日0时, 以天为单位划分为10个应力期, 用于模拟总的水均衡以及参数率定。模型拟合期选取2009年5月16日0时至19日0时, 以小时为单位分为72个应力期, 时间为3天, 用于实际泉流量与模拟泉流量对比。采用当月的岩溶地下水水位, 通过插值获得含水层的初始流场。

模型在垂直方向上分为一层,输入参数主要包括含水层顶底板高程、渗透系数、给水度、初始水头、入渗系数。模型的顶底板高程通过实际高程差值获得,模型四周为隔水边界。大气降水入渗补给量采用模拟期泉域实测资料,降雨入渗系数参照文献[11],有效降雨量按80%计算,计算研究区降雨入渗补给量。

2.2.2 Drain模块 MODFLOW中的Drain模块主要用泉口标高以及底层的传导系数来进行泉对含水层排泄的模拟,将传导系数设置足够大,以保证水都能排出来。当泉口标高低于含水层的水位时,此模块将体现出排泄地下水;当泉口标高高于含水层的水位时,此模块不起作用。

2.2.3 UZF模块用MODFLOW-2005中UZF模块模拟岩溶区表层岩溶带。UZF模型利用一维Richard方程运动波形式近似求解一维均质非饱和土壤垂直方向水分运动,计算地表水、非饱和带及饱和带之间水量交换。

一维均质非饱和土壤水流运动的Richard方程

(2)

式中:θ为土壤含水量;D(θ)为水动力扩散系数;K(θ)为非饱和渗透系数;i为蒸散发率;z为垂直方向坐标(向下为正);t为时间。

2.2.4 CFP模块 CFP模块中包含3个模式,CFPM1适用于管道流、CFPM2适用于有强透水能力的含水层、CFPM3是前两者的综合。根据研究区的水文地质条件，选用CFPM1模式。在该模式下可以较好地反映石灰岩地区的岩溶管道的形态特点,此模式适用于层流与紊流的条件,所以用于描述岩溶管道系统较为合适。CFPM1模式将传统的地下水流动方程与离散的圆柱形管道网络进行耦合。只需要在CFP模块中为各个管道赋予属性，建立独立于多孔介质的管道流控制方程,提供管道的直径d、弯曲程度τ、粗糙程度kc、管道的导水系数或管道壁的渗透系数等就可以进行模拟。

CFP模块适用于模拟S31号泉域饱水带岩溶管道内地下水的流动。前人示踪试验表明[10]，洼地内的落水洞均通过管道与S31号泉直接相连,因此，根据各洼地内落水洞分布位置及区域构造方向，可大致推断S31号泉域内管道分布,且为管道的出口。

3 模型参数确定

模型运行所需的基本数据，如顶底板标高根据实际高程差值可得,根据野外地质调查：一区的1号洼地与二区的3、4号洼地水动力条件差别很大,1号洼地相比于3、4号洼地,前者为深洼地,后者为高洼地,表明了二者岩溶发育程度差异,其中1号洼地坡度大于3、4洼地,且1号洼地内的表层岩溶泉数量比3、4号洼地内的多,所以前者比后者表层岩溶带发育。一区、二区的渗透系数确定参考场地的抽水试验以及文献[8]，设置为0.5与0.8 m/d;其给水度参照文献[14]丫吉试验场内的给水度并结合两区水文地质条件通过调整参数获得,分别设置为1×10-2与8×10-4;初始水头可插值获得;贮水率是通过野外非稳定流抽水试验,用配线法、直线图解法进行推求统一给出，为2×10-5m-1(表1)。

表1 参数分区

CFP管道模块中管道平均直径可以利用示踪试验时获得的管道平均横截面积反算; 管道壁粗糙系数参见文献[11], 因石灰岩表面比较接近混凝土,其粗糙系数为0.015; 本文假设管道是直的, 所以其管道曲折率设置为1.0; 管道节点高度可以根据试验场钻探资料获知，标高设置为160 m; 管道内地下水温根据S31号泉口的长期监测资料设置为19 ℃; 模型中管道壁的渗透系数与周边水文地质单元的渗透系数一致。

UZF模块模拟非饱和带, S31号泉域包气带最厚可超过100 m, 表层岩溶带的厚度约有3～10 m, 此部分的模拟不可或缺,非饱和带模块中降雨量的数据按每小时累计降雨输入;土壤蒸散发率参见文献[15],用蒸发能力与极限蒸散发深度的比值乘以折算系数得出,设置为1×10-4m/h;土壤极限蒸散发深度参见文献[16]，主要考虑土壤毛细管吸水深度，设置为1.5 m。

从2009年7月22日0时—8月1日0时这10天里实测降雨总量为226 mm,用来识别模型、率定参数。从图3可以看出模拟结果与实测流量有较好的拟合程度,两者的变化过程基本相同,10日内模拟总流量为1.04×105m3,实测总流量为0.89×105m3,说明模型可以用来模拟S31号泉域的降雨径流过程。

图3 识别期泉流量对比(2009年7月22日—8月1日)Fig.3 Comparison of spring discharge during the calibration period from July 22th to August 1st of 2009

4 模型结果验证及分析

4.1 泉流量拟合

运行模型后, 将模拟泉流量和实测值进行拟合, 得到泉流量拟合图4。可以看出, 对于泉流量模拟值和实际测量值的总体趋势一致, 模拟值可较好地反映泉流量的实际情况, 模型可以对该试验场岩溶峰丛洼地地区降雨径流过程进行较好的模拟。模拟时段3日内实测总泉流量为1.24×105m3，模拟总泉流量为1.21×105m3, S31号泉总流量相对误差不大,响应时间以及峰值都与实际泉流量曲线接近,满足模拟精度要求。

图4 泉流量拟合图(2009年5月16—19日)Fig.4 Fitting plot of the spring discharge from 16th to 19th on May of 2009

4.2 讨论及分析

(1)在UZF模块中,降雨入渗量是影响流量曲线最重要的因素,由于每个时间段的降雨强度不一样, 导致降雨入渗量随之变化, 表现为降雨强度越大入渗量越大, 反之越小。

图5中，降雨入渗量增大时流量曲线上升；降雨入渗量减小时, 流量曲线下降。由于输入量的不断变化, 使得非饱和带的输出是一个动态变化的过程,且降雨量与非饱和带输出量表现为正相关, 即降雨量越大, 通过非饱和带补给饱水带的量就越大。

图5 非饱和带输出量时间序列(2009年5月16—19日)Fig.5 Time series of unsaturated zone output from 16th to 19th on May of 2009

(2)从泉流量拟合图4可以看出,降雨开始于5月16日晚上,但到了17日凌晨泉流量才开始响应,说明降雨输入量先把表层岩溶带水量充填满后再由饱水带运动至泉口,需要经过这一个过程才能使得泉流量上升。通过时间轴可以看出,非饱和带输出是滞后于降雨输入的,在暴雨条件下滞后时间较短,约2小时左右,在中雨或小雨条件下,滞后时间更长一些,体现出表层岩溶带的调蓄作用。

(3)从非饱和带输出量与饱水带储存量的关系(图6)可以了解到,通过非饱和带输出的量随着降雨增大而增大,其输出量一部分主要转化为泉流量,另一部分主要转化为饱水带储存量,饱水带储存量随之增大,但是饱水带储存量最终也会通过泉排泄,此时饱水带储存量与泉流量表现为增大趋势。当降雨减少或者停止时,饱水带储存量增量随之减少,同时消耗饱水带存储量,不断补给泉流量,此时储存量表现为减少,泉流量表现出平缓下降。

图6 非饱和带输出量与饱水带储存量(2009年5月16—19日)Fig.6 Unsaturated zone output and saturated zone storage from 16th to 19th on May of 2009

(4)泉域水流方向为由北东东向南西西方向,在泉口附近垂直方向上取两个点对比水头值,其中一个位于管道模块上,另一个则位于含水层基质内。由图7可见，在暴雨条件下,模型中CFP管道模块响应较快,大部分降雨表现为直接由消水洞通过管道向泉口排泄,曲线上升比无管道地方快且早。在小雨或无降雨条件下,非饱和带内储存的水量缓慢补给饱水带,曲线下降趋势较慢, 其主要受裂隙系统影响,管道系统对其影响很小。

5 结论

研究结果表明,同时利用UZF与CFP模块对丫吉试验场岩溶峰丛洼地地区S31号泉流量的模拟是较可信的。本文采用UZF软件包模拟非饱和带水分运移,更为准确地反映降雨入渗补给过程,也体现出非饱和带对降雨的调蓄作用,在此基础上添加CFP管道流模块,对岩溶区泉流量模拟,更好地刻画了岩溶泉的水文过程。由于岩溶区含水层的复杂性与实测资料缺少,所以假设管道与实际管道存在一定差别,无论是管道大小，还是管道弯曲程度及有无分支,都不同程度地影响着泉流量的变化过程。在岩溶区的裂隙与管道的不规律分布以及强烈的非均质性,也让模型结构存在不确定性,使模型模拟结果具有一定的误差。但总体来说,模型较好体现了UZF模块对降雨输入的调蓄作用,对泉流量峰值的滞后作用,以及CFP模块在暴雨条件下体现出泉流量对降雨的快速响应,使得岩溶泉流量模拟更加精确可信。