辛安泉泉水流量數學模擬研究

牛二偉

（長治市水文水資源勘測分局）

辛安泉是山西省19個巖溶大泉之一，位于晉東南太行山中段西側。泉域總面積10 950 km2，地跨長治和晉中2市，其中長治市境內面積9 430 km2，占到86.12%，涉及全市除沁源縣以外的11個縣（區(qū)）。泉域裸露巖溶區(qū)面積2 200 km2，占總面積20.09%；隱伏巖溶區(qū)面積8 750 km2，占總面積79.91%。以泉域集中出露帶為主體的長治市城區(qū)辛安泉飲用水水源地，為當地最為主要的供水水源地，已被列入全國重要水源地名錄。2000年以來，辛安泉泉水流量與供水能力呈現下降趨勢，不僅直接威脅到當地居民的生活飲水安全，而且大大制約了社會經濟發(fā)展。因此，開展辛安泉泉水流量研究，分析泉水流量的主要影響因素，為切實保護辛安泉水資源提供堅實的科學依據，是非常必要的。

1 辛安泉域自然環(huán)境特征與泉水流量

辛安泉是一個泉群，為全排型巖溶大泉，泉水初露于潞城區(qū)辛安泉鎮(zhèn)濁漳河河谷地帶。泉域巖溶水受來自于北、西、南方向的水力梯度作用向排泄區(qū)匯流，形成王曲和實會兩大泉群。辛安泉泉域地形復雜，在不同的地質作用下形成了以山間盆地為中心、四周山區(qū)丘陵環(huán)抱的地貌景觀，其中山地丘陵面積占泉域面積90%以上。

1.1 泉域地質特征

辛安泉泉域構造為單斜構造，傾向北西，自西向東依次為沁水塊坳、長治新裂陷盆地、晉獲褶斷帶。泉域內斷層、地壘等地質構造為泉域巖溶水的主要控制因素，以長治大斷層和北耽車斷層最為明顯;其次為王山地壘、二崗山地壘等。區(qū)域出露地層，包括太古界古老變質巖，元古界礫巖、石英砂巖和頁巖，寒武—奧陶系碳酸鹽巖，石炭—二疊系礫巖、砂巖、頁巖和泥巖，三疊系砂巖，侏羅系石英砂巖，以及第三系至第四系松散巖層。泉域主要含水層，以奧陶系中統(tǒng)灰?guī)r、泥灰?guī)r為主構成；次要含水層，以寒武系中上統(tǒng)灰?guī)r、白云質灰?guī)r為主構成；相對隔水層，則由奧陶系下統(tǒng)（排泄區(qū)除外）為主構成。

1.2 泉域巖溶水系統(tǒng)特征

辛安泉泉域作為一個完整的水文地質單元，具有獨立的補、徑、排途徑。結合泉域內水文地質及富水性情況，可將泉域劃分為四個區(qū)域，即：強富水區(qū)、較強富水區(qū)、一般富水區(qū)和弱富水區(qū)。強富水區(qū)位于辛安泉排泄區(qū)；較強富水區(qū)分布于長治—潞城、襄垣、黎城三個徑流區(qū)；一般富水區(qū)分布于泉域強徑流帶外側和地下水補給區(qū)或補給徑流區(qū)；弱富水區(qū)主要位于泉域的灰?guī)r裸露區(qū)邊緣地帶和長子、屯留以西碳酸鹽巖類底層深埋區(qū)。辛安泉域泉點170余個，可分為王曲泉群和實會泉群兩大部分。泉域內各泉組分布于西流村至北耽車之間約16 km長的河谷中，出露標高600-615m。

1.3 泉域補、徑、排循環(huán)特征

泉域水資源補給來源主要為大氣降水和地表水滲漏補給。其中：大氣降水為主要補給來源，在碳酸巖類地區(qū)尤為明顯；其次，為地表水的滲漏補給。含水層巖性，以松散巖類、碎屑巖類、碳酸鹽巖類為主。泉域碳酸鹽巖類裂隙巖溶水的運動形式，主要是從不同位置及不同的方向，以各種不同的途徑，向排泄區(qū)匯流，并呈泉群出露。由于濁漳河干流對太行山的強烈侵蝕作用，使潞城區(qū)西流至平順縣北耽車一帶濁漳河河谷，成為辛安泉域內標高最低的排泄基準面。在辛安泉灣里一帶出露下奧陶統(tǒng)白云巖弱透水層，使上游富水性較強的中奧陶統(tǒng)巖溶地下水受阻而匯集溢流，形成了王曲泉群。北耽車以下，沿河床出露的長城系、寒武系下統(tǒng)隔水層具有良好的阻水作用，形成實會泉群。

1.4 泉域泉水流量動態(tài)特征

通過收集辛安泉歷年實測流量資料，可以準確地反映出泉水流量的變化趨勢。歷年流量監(jiān)測情況見表1。根據1956-2017年62 a觀測資料，辛安泉實測多年平均流量8.21m3/s，最大流量16.03m3/s（1964年），最小流量2.83m3/s（2012年）。根據資料分析，泉水流量總體上呈減少趨勢，從20世紀80年代開始，泉水流量明顯減少（見圖1）。

表1 辛安泉歷年流量監(jiān)測基本情況表

圖1 辛安泉泉水歷年天然流量變化圖

2 泉水流量影響因素分析

根據辛安泉長系列泉水流量觀測資料分析，認為流量呈現減少的原因主要有以下方面。

2.1 降水量因素

辛安泉排泄區(qū)多為碳酸鹽裸露區(qū)，入滲補給速度快，年降水量與泉水流量顯著相關。年降水量多時，地下水補給量就大，泉水流量也大。

2.2 開采量因素

隨著城市擴大和工農業(yè)生產的發(fā)展，泉域內巖溶水開采量不斷增大，對泉水流量影響亦逐漸增大。1980年開始開采巖溶水，開采量只有0.03m3/s；1987年長治市城市供水公司和天脊集團兩個大型水源地的投入開采，開采量增大到了1.2m3/s；2017年開采量達到了3.12m3/s。

2.3 煤礦礦井排水因素

石炭、二疊系裂隙含水層的地下水，可沿裂隙向深部巖溶含水層產生越流補給，當地采煤排水使得煤系地層的地下水被疏干，從而影響到地下水的下滲補給。

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

除長治市水文水資源勘測分局在石梁、天橋設立的2處水文測站對辛安泉泉水流量進行觀測外，在1993年，辛安泉域管理中心在排泄區(qū)設立了王曲水文站，對葦蒲泉、南流二泉、南流三泉、西流泉、戰(zhàn)備渠、濁漳河主河道泉水出露區(qū)上游南堡斷面、泉水匯入后濁漳河下游安樂斷面以及監(jiān)測區(qū)間農業(yè)用水與天脊集團退水等辛安泉域的8個泉水流量觀測斷面進行了觀測，觀測斷面分布見圖2。

據對1956-2017年觀測資料統(tǒng)計，年均降水量561.8mm，天然泉水流量9.46m3/s，折算徑流深約27.2mm，泉水產流系數為0.048 4。

3.1 泉水對降水響應的滯后效應

一般發(fā)生降水后，需經過入滲補給過程，才能進入含水層中。因此，降水對泉水流量的影響存在著一定的時間差，即滯后性。采用錯位時間不同相關度分析降水量、泉水流量逐月數據，可以大致推斷出降水影響的滯后時間長度（見圖3）。

從圖3可以看出，互相關曲線第一個峰值對應的時間為λ=7個月，除7個月外，曲線峰值所對應的時間還包括19、31、43個月……。從峰值對應的時間序列可以得出，降水量對泉水流量影響滯后周期大致為12個月。該現象可歸結于辛安泉巖溶水系統(tǒng)降水量的年內周期性和相似性。根據最小峰值時間可知，降水對泉水流量的影響滯后時間在12個月內。

天然降水量對泉水流量的影響，不僅表現為當年降水量，還包含上年降水量。為便于分析，將當年降水量值與相應的泉水徑流深繪制過程線圖，同時增加近2 a平均降水量、近3 a平均降水量、近4 a平均降水量以及近n年平均降水量。據分析，近2 a平均降水量起伏與泉水徑流深基本一致，與大氣降水對辛安泉流量的影響具有滯后效應結論一致。

圖2 辛安泉域王曲水文站測驗斷面分布圖

圖3 泉水流量與降水量互相關分析曲線

3.2 降水量與徑流深關聯分析

泉水流量與降水量關系密切，據此進行相關分析。根據泉水流量和降水量觀測值，建立二者之間的相關方程。將歷年累計降水量與累計泉水徑流深點繪在一起（圖4），發(fā)現趨勢基本吻合。二者間的相關方程為：y=0.049 5x+59.012，相關系數R2=0.991 6，說明高度相關。

結合泉水開采情況，進一步細分析發(fā)現，1956-1982年泉水未大量開采時，關系曲線的斜率較小，為低斜率區(qū)；1982年后受泉水開采影響，斜率明顯增大，1999年以來斜率更大。

4 辛安泉泉水流量數學模擬研究

辛安泉屬全排型泉水，泉水流量主要來自降水入滲補給。根據上述分析，泉水流量不但受當年降水量影響，而且受上年甚至幾年降水量影響。除降水量影響之外，還受上游地下水開采量、煤礦礦井水排放影響?？紤]大氣降水輸入的滯后及延遲時間，以年均降水量與泉水徑流深建立的關系來推求泉水徑流深，分時段分析結果見圖5和表2。

由此可知：1956-1985年泉水未大量開采時的徑流系數較大，其中：1976-1982年為大系數區(qū)，1961-1962年為小系數區(qū)；1982年后受降水與泉水開采影響徑流系數逐年減小，且1999-2017年為系數最小區(qū)。

圖4 歷年累計降雨量與累計泉水徑流深關系圖

圖5 辛安泉泉域降水徑流系數變化曲線圖

表2 1956-2017年辛安泉泉水徑流與區(qū)域降水量相關系數分析表

巖溶水系統(tǒng)的輸入項為大氣降水量，輸出項為辛安泉流量。模型的輸入值選取泉水流量（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年的泉水流量進行預測。經統(tǒng)計，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年后出現連續(xù)豐水年情況下辛安泉未來的流量；當采用出現頻率為50%的降水量值（563.8mm）作為今后幾年的降水量數據時，預測結果可代表2019年后出現連續(xù)平水年情況下辛安泉未來的流量；而當采用出現頻率為90%的降水量值（425.5mm）時，預測結果則代表連年干旱氣候條件下辛安泉流量的變化趨勢。不同降水量保證率下的泉水流量預測結果見表4。

表3 降水量理論頻率計算結果

表4 2019-2025年泉水流量預測結果

按照目前的地下水開采強度，未來的降水量直接影響巖溶泉水流量。若未來的降水量不出現極端情況，僅在563.8mm左右波動時，辛安泉實測流量將在4.69m3/s左右波動；當出現多年連續(xù)干旱的情況時，即使不再增加地下水開采強度，辛安泉流量也將呈下降趨勢，實測流量將在2.77m3/s左右波動。泉水流量和開采量相關，開采量增加，泉水流量將減??；開采量保持平衡，則泉水流量波動較小。如果加大開采量，泉水流量下降趨勢將更為明顯。

5 結語

辛安泉作為長治市城鄉(xiāng)最為重要的飲用水集中供水水源地，其泉水流量大小與水質變化，倍受當地關注。開展泉水流量模擬研究，預測發(fā)展動向，及時采取控制開采措施，同時加強泉域水土保持與水生態(tài)環(huán)境修復改善，對于向城鄉(xiāng)持續(xù)供水及促進市域社會經濟發(fā)展，均具有重大的現實意義與戰(zhàn)略意義。作為水文部門進行的這項理論研究，是一次有益的大膽嘗試，盡管考慮的因素還較少，所建模型還不夠精細，預測的結果還有待時間驗證。下一步，我們將從水量、水質、水環(huán)境等方面對辛安泉進行深入研究，以便為地方政府在供水與水環(huán)境保護宏觀決策方面，提供更加確切的科學依據。