考慮非飽和帶作用及管道流的巖溶泉流量模擬

覃夏南,姜光輝,夏 源

(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模塊在暴雨條件下體現出泉流量對降雨的快速響應,使得巖溶泉流量模擬更加精確可信。

圖7 管道與非管道上水頭對比(2009年5月16—19日)Fig.7 Head comparison of pipe and non-conduit pipe from 16th to 19th on May of 2009