辛安泉泉水流量数学模拟研究

牛二伟

（长治市水文水资源勘测分局）

辛安泉是山西省19个岩溶大泉之一，位于晋东南太行山中段西侧。泉域总面积10 950 km2，地跨长治和晋中2市，其中长治市境内面积9 430 km2，占到86.12%，涉及全市除沁源县以外的11个县（区）。泉域裸露岩溶区面积2 200 km2，占总面积20.09%；隐伏岩溶区面积8 750 km2，占总面积79.91%。以泉域集中出露带为主体的长治市城区辛安泉饮用水水源地，为当地最为主要的供水水源地，已被列入全国重要水源地名录。2000年以来，辛安泉泉水流量与供水能力呈现下降趋势，不仅直接威胁到当地居民的生活饮水安全，而且大大制约了社会经济发展。因此，开展辛安泉泉水流量研究，分析泉水流量的主要影响因素，为切实保护辛安泉水资源提供坚实的科学依据，是非常必要的。

1 辛安泉域自然环境特征与泉水流量

辛安泉是一个泉群，为全排型岩溶大泉，泉水初露于潞城区辛安泉镇浊漳河河谷地带。泉域岩溶水受来自于北、西、南方向的水力梯度作用向排泄区汇流，形成王曲和实会两大泉群。辛安泉泉域地形复杂，在不同的地质作用下形成了以山间盆地为中心、四周山区丘陵环抱的地貌景观，其中山地丘陵面积占泉域面积90%以上。

1.1 泉域地质特征

辛安泉泉域构造为单斜构造，倾向北西，自西向东依次为沁水块坳、长治新裂陷盆地、晋获褶断带。泉域内断层、地垒等地质构造为泉域岩溶水的主要控制因素，以长治大断层和北耽车断层最为明显;其次为王山地垒、二岗山地垒等。区域出露地层，包括太古界古老变质岩，元古界砾岩、石英砂岩和页岩，寒武—奥陶系碳酸盐岩，石炭—二叠系砾岩、砂岩、页岩和泥岩，三叠系砂岩，侏罗系石英砂岩，以及第三系至第四系松散岩层。泉域主要含水层，以奥陶系中统灰岩、泥灰岩为主构成；次要含水层，以寒武系中上统灰岩、白云质灰岩为主构成；相对隔水层，则由奥陶系下统（排泄区除外）为主构成。

1.2 泉域岩溶水系统特征

辛安泉泉域作为一个完整的水文地质单元，具有独立的补、径、排途径。结合泉域内水文地质及富水性情况，可将泉域划分为四个区域，即：强富水区、较强富水区、一般富水区和弱富水区。强富水区位于辛安泉排泄区；较强富水区分布于长治—潞城、襄垣、黎城三个径流区；一般富水区分布于泉域强径流带外侧和地下水补给区或补给径流区；弱富水区主要位于泉域的灰岩裸露区边缘地带和长子、屯留以西碳酸盐岩类底层深埋区。辛安泉域泉点170余个，可分为王曲泉群和实会泉群两大部分。泉域内各泉组分布于西流村至北耽车之间约16 km长的河谷中，出露标高600-615m。

1.3 泉域补、径、排循环特征

泉域水资源补给来源主要为大气降水和地表水渗漏补给。其中：大气降水为主要补给来源，在碳酸岩类地区尤为明显；其次，为地表水的渗漏补给。含水层岩性，以松散岩类、碎屑岩类、碳酸盐岩类为主。泉域碳酸盐岩类裂隙岩溶水的运动形式，主要是从不同位置及不同的方向，以各种不同的途径，向排泄区汇流，并呈泉群出露。由于浊漳河干流对太行山的强烈侵蚀作用，使潞城区西流至平顺县北耽车一带浊漳河河谷，成为辛安泉域内标高最低的排泄基准面。在辛安泉湾里一带出露下奥陶统白云岩弱透水层，使上游富水性较强的中奥陶统岩溶地下水受阻而汇集溢流，形成了王曲泉群。北耽车以下，沿河床出露的长城系、寒武系下统隔水层具有良好的阻水作用，形成实会泉群。

1.4 泉域泉水流量动态特征

通过收集辛安泉历年实测流量资料，可以准确地反映出泉水流量的变化趋势。历年流量监测情况见表1。根据1956-2017年62 a观测资料，辛安泉实测多年平均流量8.21m3/s，最大流量16.03m3/s（1964年），最小流量2.83m3/s（2012年）。根据资料分析，泉水流量总体上呈减少趋势，从20世纪80年代开始，泉水流量明显减少（见图1）。

表1 辛安泉历年流量监测基本情况表

图1 辛安泉泉水历年天然流量变化图

2 泉水流量影响因素分析

根据辛安泉长系列泉水流量观测资料分析，认为流量呈现减少的原因主要有以下方面。

2.1 降水量因素

辛安泉排泄区多为碳酸盐裸露区，入渗补给速度快，年降水量与泉水流量显著相关。年降水量多时，地下水补给量就大，泉水流量也大。

2.2 开采量因素

随着城市扩大和工农业生产的发展，泉域内岩溶水开采量不断增大，对泉水流量影响亦逐渐增大。1980年开始开采岩溶水，开采量只有0.03m3/s；1987年长治市城市供水公司和天脊集团两个大型水源地的投入开采，开采量增大到了1.2m3/s；2017年开采量达到了3.12m3/s。

2.3 煤矿矿井排水因素

石炭、二叠系裂隙含水层的地下水，可沿裂隙向深部岩溶含水层产生越流补给，当地采煤排水使得煤系地层的地下水被疏干，从而影响到地下水的下渗补给。

3 辛安泉泉域水文资料分析

除长治市水文水资源勘测分局在石梁、天桥设立的2处水文测站对辛安泉泉水流量进行观测外，在1993年，辛安泉域管理中心在排泄区设立了王曲水文站，对苇蒲泉、南流二泉、南流三泉、西流泉、战备渠、浊漳河主河道泉水出露区上游南堡断面、泉水汇入后浊漳河下游安乐断面以及监测区间农业用水与天脊集团退水等辛安泉域的8个泉水流量观测断面进行了观测，观测断面分布见图2。

据对1956-2017年观测资料统计，年均降水量561.8mm，天然泉水流量9.46m3/s，折算径流深约27.2mm，泉水产流系数为0.048 4。

3.1 泉水对降水响应的滞后效应

一般发生降水后，需经过入渗补给过程，才能进入含水层中。因此，降水对泉水流量的影响存在着一定的时间差，即滞后性。采用错位时间不同相关度分析降水量、泉水流量逐月数据，可以大致推断出降水影响的滞后时间长度（见图3）。

从图3可以看出，互相关曲线第一个峰值对应的时间为λ=7个月，除7个月外，曲线峰值所对应的时间还包括19、31、43个月……。从峰值对应的时间序列可以得出，降水量对泉水流量影响滞后周期大致为12个月。该现象可归结于辛安泉岩溶水系统降水量的年内周期性和相似性。根据最小峰值时间可知，降水对泉水流量的影响滞后时间在12个月内。

天然降水量对泉水流量的影响，不仅表现为当年降水量，还包含上年降水量。为便于分析，将当年降水量值与相应的泉水径流深绘制过程线图，同时增加近2 a平均降水量、近3 a平均降水量、近4 a平均降水量以及近n年平均降水量。据分析，近2 a平均降水量起伏与泉水径流深基本一致，与大气降水对辛安泉流量的影响具有滞后效应结论一致。

图2 辛安泉域王曲水文站测验断面分布图

图3 泉水流量与降水量互相关分析曲线

3.2 降水量与径流深关联分析

泉水流量与降水量关系密切，据此进行相关分析。根据泉水流量和降水量观测值，建立二者之间的相关方程。将历年累计降水量与累计泉水径流深点绘在一起（图4），发现趋势基本吻合。二者间的相关方程为：y=0.049 5x+59.012，相关系数R2=0.991 6，说明高度相关。

结合泉水开采情况，进一步细分析发现，1956-1982年泉水未大量开采时，关系曲线的斜率较小，为低斜率区；1982年后受泉水开采影响，斜率明显增大，1999年以来斜率更大。

4 辛安泉泉水流量数学模拟研究

辛安泉属全排型泉水，泉水流量主要来自降水入渗补给。根据上述分析，泉水流量不但受当年降水量影响，而且受上年甚至几年降水量影响。除降水量影响之外，还受上游地下水开采量、煤矿矿井水排放影响。考虑大气降水输入的滞后及延迟时间，以年均降水量与泉水径流深建立的关系来推求泉水径流深，分时段分析结果见图5和表2。

由此可知：1956-1985年泉水未大量开采时的径流系数较大，其中：1976-1982年为大系数区，1961-1962年为小系数区；1982年后受降水与泉水开采影响径流系数逐年减小，且1999-2017年为系数最小区。

图4 历年累计降雨量与累计泉水径流深关系图

图5 辛安泉泉域降水径流系数变化曲线图

表2 1956-2017年辛安泉泉水径流与区域降水量相关系数分析表

岩溶水系统的输入项为大气降水量，输出项为辛安泉流量。模型的输入值选取泉水流量（1956-2017）、降水量及岩溶水开采量。其中降水量采用流域内雨量站水资源公报选用站点的降水量均值。所建数学模型如下：

式中：Q实测——泉水实测流量，m3/s；

Q天——天然泉水流量，m3/s；

Q开——泉水开采量，m3/s；

R泉——泉水年径流深，mm；

A——泉域面积，采用10 950 km2；

T——时间，s；

a——降水径流系数；

P——年降水量，mm。

应用该模型，先采用近2 a降水量P通过（3）式求得辛安泉年径流深，再通过（2）式求得辛安泉泉水流量，用（1）式减去当年辛安泉开采利用量，即得到辛安泉的实测泉水出露流量。

经计算，预测值与实际径流深误差在-4.48-3.84mm之间，占实测径流深13%-14%。相对误差在10%以内的样本为48个，占总样本比例78.7%；相对误差在10%-15%之间的样本为13个，占总样本比例21.3%。具体预测结果见图6。

图6 辛安泉泉水径流深预测值与实测值比较图

采用所建的数学模型，对辛安泉2018年的泉水流量进行预测。经统计，2018年1-10月降水量为422.0mm，11-12月采用历年平均值22.3mm，则2018年降水量为444.3mm，2017年降水量614.8mm，预测采用2 a降水量的平均值。将2017年的岩溶水开采量（3.12m3/s）作为2018年开采量，取a=0.034 0。则有：

R泉=529.6×0.034 0=18.0mm

Q天=R泉×A/T×10-3=6.25m3/s

Q实测=Q天-Q开=6.25-3.12=3.13m3/s。

即：2018年实测泉水流量为3.13m3/s。

现应用该数学模型，对2019-2025年间的泉水流量变化趋势进行预测。在预测前，需要获取预测时段的降水量与开采量数据，而目前还难以提前多年做出准确的降水预报，拟利用频率分析方法来获取未来的降水量数据。即：首先以1956-2017年间62 a的降水量资料为原始数据，通过理论频率法计算出不同出现频率的降水量值（表3）；再选取出现频率为20%（丰水年）、50%（平水年）和90%（枯水年）的降水量，分别作为预测时段可能的降水量值，进而对不同降水保证率条件下的泉水流量进行预测。至于地下水开采量的处理，则采取如下方式：假设2019-2025年间的岩溶水开采强度不会进一步增加，将2017年的岩溶水开采量（3.12m3/s）作为今后的开采量。

当采用出现频率为20%的降水量值（664.9mm）作为今后几年的降水量数据时，预测结果可代表2019年后出现连续丰水年情况下辛安泉未来的流量；当采用出现频率为50%的降水量值（563.8mm）作为今后几年的降水量数据时，预测结果可代表2019年后出现连续平水年情况下辛安泉未来的流量；而当采用出现频率为90%的降水量值（425.5mm）时，预测结果则代表连年干旱气候条件下辛安泉流量的变化趋势。不同降水量保证率下的泉水流量预测结果见表4。

表3 降水量理论频率计算结果

表4 2019-2025年泉水流量预测结果

按照目前的地下水开采强度，未来的降水量直接影响岩溶泉水流量。若未来的降水量不出现极端情况，仅在563.8mm左右波动时，辛安泉实测流量将在4.69m3/s左右波动；当出现多年连续干旱的情况时，即使不再增加地下水开采强度，辛安泉流量也将呈下降趋势，实测流量将在2.77m3/s左右波动。泉水流量和开采量相关，开采量增加，泉水流量将减小；开采量保持平衡，则泉水流量波动较小。如果加大开采量，泉水流量下降趋势将更为明显。

5 结语

辛安泉作为长治市城乡最为重要的饮用水集中供水水源地，其泉水流量大小与水质变化，倍受当地关注。开展泉水流量模拟研究，预测发展动向，及时采取控制开采措施，同时加强泉域水土保持与水生态环境修复改善，对于向城乡持续供水及促进市域社会经济发展，均具有重大的现实意义与战略意义。作为水文部门进行的这项理论研究，是一次有益的大胆尝试，尽管考虑的因素还较少，所建模型还不够精细，预测的结果还有待时间验证。下一步，我们将从水量、水质、水环境等方面对辛安泉进行深入研究，以便为地方政府在供水与水环境保护宏观决策方面，提供更加确切的科学依据。