王 喆，盧 麗，夏日元，石勝良

(1.中國地質(zhì)科學院a.巖溶地質(zhì)研究所;b.國土資源部巖溶動力學重點實驗室，廣西桂林 541004;2.廣東省地質(zhì)礦產(chǎn)公司，廣州 510083)

1 研究背景

我國巖溶地區(qū)發(fā)育面積較大，約占全國面積的1/3左右，在這些地區(qū)發(fā)育著各種泉，這些泉在我國經(jīng)濟發(fā)展的過程中，發(fā)揮著重要的作用。因此需要準確地模擬出泉流量，這樣才能為泉域的水資源地保護和開發(fā)利用提供必要的依據(jù)。目前，建立泉流量模型一般有4種方法:線性法、非線性法、數(shù)值法和耦合法。線性法一般包括多元回歸法［1］、偏最小二乘法等;非線性法包括灰色預測模型法［2］、神經(jīng)網(wǎng)絡法［3］、時間序列分析法［4］等;數(shù)值法主要是運用MODFLOW軟件進行預測分析［5］;耦合法則是目前較為流行的方法，有的是線性法和非線性法的耦合［6］，有的則是非線性法和數(shù)值法的耦合［7］。從上述論述可以看出，耦合法可以將2種或者多種方法進行融合，從而克服因算法本身缺陷所產(chǎn)生的問題，可以提高模型的精度，因此耦合法是泉流量模擬研究的發(fā)展趨勢。目前的耦合模型較少關注泉流量與巖溶含水系統(tǒng)演化的關系，只是追求泉流量的擬合效果，而且從未有過學者在分析巖溶含水系統(tǒng)的演化與泉流量關系中考慮過巖溶溶蝕過程，但溶蝕過程卻是一個重要階段，不能忽略。

因此本次研究利用基于VC++6.0的Fracture-ToKarst程序?qū)θ蛳到y(tǒng)的滲流模型和溶蝕模型進行耦合，分均質(zhì)和非均質(zhì)2種情況構建模型，并對泉流量進行分析。此程序適用于非連續(xù)裂隙網(wǎng)絡介質(zhì)模型，具體耦合過程是首先利用已知條件對滲流場進行計算，得到每一段裂隙水頭和流速，計算出鈣離子的濃度，依據(jù)公式算出每一段裂隙擴展速度及一個步長后裂隙的新寬度，在新裂隙寬度條件下重新計算滲流場，并依據(jù)上述步驟進行下一個步長的裂隙寬度計算，從而使?jié)B流模型和溶蝕模型達到耦合。

2 泉域系統(tǒng)演化數(shù)學模型構建

在構建耦合模型前必須選擇滲流模型。目前，巖溶地區(qū)滲流模型有3種:連續(xù)介質(zhì)模型、不連續(xù)介質(zhì)模型和耦合介質(zhì)模型。本文選用非連續(xù)裂隙網(wǎng)絡介質(zhì)模型。該模型認為裂隙與巖體間不存在水量交換，水流只在等效裂隙網(wǎng)絡中運動［8］。

模型的構建過程大致如下:首先將研究區(qū)的裂隙分為確定性裂隙和隨機裂隙2種，確定性裂隙直接輸入模型，隨機性裂隙通過Monte-Carlo產(chǎn)生，這2種裂隙構成原始裂隙網(wǎng)絡。將沒有導水能力的裂隙進行剔除，構建連通型裂隙網(wǎng)絡，接下來對裂隙網(wǎng)絡進行離散化處理。裂隙網(wǎng)絡建好后，利用滲流模型計算出裂隙節(jié)點水頭，構建初始的滲流場，然后利用初始滲流場和裂隙網(wǎng)絡來計算巖溶系統(tǒng)中裂隙的溶蝕過程，推演泉域系統(tǒng)演化過程，進而預測泉流量。

2.1 模型假設

(1)介質(zhì)場為等效裂隙網(wǎng)絡模型，由確定性裂隙和隨機裂隙構成，確定性裂隙直接輸入，隨機裂隙由Monte-Carlo法生成，裂隙兩兩正交，裂隙段假設為光滑平板裂隙。

(2)滲流場中的水流為穩(wěn)定流態(tài)。

(3)溶蝕場的碳酸鹽溶蝕速率，利用White的半經(jīng)驗公式進行計算。

2.2 裂隙網(wǎng)絡模擬

裂隙網(wǎng)絡分確定性裂隙和隨機裂隙。對于確定性裂隙，根據(jù)長度、寬度、產(chǎn)狀等參數(shù)，構建確定性裂隙網(wǎng)絡，對于數(shù)量眾多的隨機裂隙，采用 Monte-Carlo模擬技術［9］生成隨機裂隙網(wǎng)絡，將確定性裂隙網(wǎng)絡和隨機裂隙網(wǎng)絡耦合在一起，統(tǒng)一編號，然后對統(tǒng)一后的裂隙網(wǎng)絡進行離散化處理。

2.3 裂隙網(wǎng)絡中滲流場模擬

根據(jù)水流為穩(wěn)定流態(tài)和裂隙段為光滑平板裂隙這2個假設條件，可得出其水流方程滿足立方定律，即

式中:q是單寬流量;ρ是質(zhì)量密度;μ是流動黏滯系數(shù);b為裂隙張開度;H是裂隙兩端節(jié)點的水頭。

按水量平衡原理，節(jié)點i處的穩(wěn)定流水流方程為

式中:qj為裂隙j單元流入(為正)或流出(為負)節(jié)點i的流量;Qi為節(jié)點i處的源匯項。

將式(1)和式(2)聯(lián)立起來，采用迭代方法進行求解，進而得出滲流場。

2.4 巖溶裂隙網(wǎng)絡溶蝕模擬

由模型假設可知，本次研究裂隙形態(tài)是基于平行板模型的單裂隙，因此本文采用White［10］所提出的單裂隙溶蝕擴展數(shù)學方程，即

式中:F是溶蝕速率(mol·cm-2·s-1);Kc是溶蝕速度常數(shù)(cm4·mol-1·s-1);Ceq是Ca2+的飽和濃度;C是水中某時刻Ca2+濃度。

在一段裂隙中，溶液帶走的鈣離子的數(shù)量可以表示為

式中:Q是每段裂隙中的流量;Δt是溶蝕作用的時間步長;C0是溶液流入一段裂隙時的初始鈣離子濃度。

由質(zhì)量守恒定律和式(4)可以推算出，單裂隙張開度擴張速度為

式中:ΔB是每一個步長內(nèi)裂隙張開度變化量;M是在一段裂隙中溶液帶走的鈣離子的數(shù)量;ρ1是可溶巖密度;L是裂隙段長度。

如圖1所示，對于一個裂隙均衡單元，假如裂隙交叉點i0的水頭及Ca2+濃度為H和C，其相鄰交叉點 in(n=1，2，3，4)的水頭及 Ca2+濃度為 Hn，Cn，則裂隙交叉點i0與其鄰點in有

其中:

式中:Ci為鄰點in的鈣離子濃度;ΔCi為水流從節(jié)點in到節(jié)點i0的鈣離子濃度增加量;qi為節(jié)點in到節(jié)點i0的流量;Li為節(jié)點in到節(jié)點i0的裂隙跡長;Bi為節(jié)點in和節(jié)點i0的寬度。

上述公式的推導過程參見袁道先所寫的《中國巖溶動力系統(tǒng)》［11］。對式(5)采用迭代法進行求解。具體步驟如下:首先利用式(6)至式(9)迭代算出節(jié)點Ca2+濃度，然后在根據(jù)式(4)和式(5)計算出裂隙張開度擴張速度。

圖1 裂隙網(wǎng)絡離散單元Fig.1 A unit in discrete fracture networks

3 泉域系統(tǒng)演化模擬與分析

在已有數(shù)學模型的基礎上，結(jié)合補徑排、邊界條件和裂隙網(wǎng)絡等情況，對泉域的演化進行模擬和分析。

補徑排條件:假設模型只接受降雨入滲補給和河流側(cè)向補給。這其中，降雨入滲補給的范圍是圖2(a)綠線部分，年降雨量為400 mm，降雨入滲系數(shù)為0.15，模擬期內(nèi)降雨入滲補給量不變(本次模擬期為5 000 a);河流側(cè)向補給的范圍是圖2(a)藍線部分，水位高程為30 m，排泄河流高程為10.5 m，河流側(cè)向補給量在模擬期內(nèi)不變，且小于降雨入滲補給量。

圖2 巖溶泉域系統(tǒng)概念模型Fig.2 Conceptual model of the karst spring system

邊界條件:補給河流設定為定水頭邊界，補給河流與排泄河流之間的地面線處的邊界設定為降雨入滲補給邊界，泉水點設定為自由排泄邊界，泉水點之間的節(jié)點也按自由排泄邊界進行處理，而系統(tǒng)左邊、底邊及右下邊界為隔水邊界(如圖2(a)中紅線所示)。

裂隙網(wǎng)絡:由確定性裂隙和隨機裂隙組成的(圖2(b))，圖2中所顯示裂隙網(wǎng)絡的只是形象化顯示，并非實際情況。真實的裂隙網(wǎng)絡分均質(zhì)系統(tǒng)和非均質(zhì)系統(tǒng)2種情況，均質(zhì)系統(tǒng)裂隙網(wǎng)絡無確定性裂隙，由隨機裂隙組成，且根據(jù)模型假設條件，隨機裂隙兩兩垂直，且水平裂隙和垂直裂隙跡長分別為400 cm和200 cm，間距分別為200 cm和100 cm，表示為Lh和Lv，橫縱向的張開度均為0.02 cm，裂隙條數(shù)分別為841條和735條;非均質(zhì)系統(tǒng)裂隙網(wǎng)絡由一條確定性裂隙和隨機裂隙組成，其中確定性裂隙跡長為100 m，在侵蝕基準面處，張開度為0.08 cm，用ap來表示，而隨機裂隙則與均質(zhì)系統(tǒng)的相同。

3.1 均質(zhì)性巖溶泉域系統(tǒng)演化

根據(jù)上節(jié)所介紹的裂隙網(wǎng)絡情況，均質(zhì)性巖溶泉域系統(tǒng)的裂隙網(wǎng)絡是通過如下步驟進行構建:首先根據(jù)模型假設，均質(zhì)系統(tǒng)的隨機裂隙分為水平裂隙和垂直裂隙，然后根據(jù)上節(jié)所介紹的裂隙參數(shù)情況將水平裂隙和垂直裂隙兩者的跡長、間距、傾角和張開度等參數(shù)輸入FractureToKarst程序中，然后根據(jù)Monte-Carlo法構造各參數(shù)的隨機數(shù)發(fā)生器，進而形成裂隙網(wǎng)絡。

均質(zhì)性巖溶泉域系統(tǒng)演化過程如圖3所示，圖內(nèi)的紅線為含水系統(tǒng)內(nèi)等水位線。下面詳述一下各個時間點上含水系統(tǒng)水位和張開度的變化情況。

圖3 均質(zhì)性巖溶泉域系統(tǒng)演化過程Fig.3 Evolution of a ho mogeneous karst spring system

100 a:巖溶泉域系統(tǒng)由于處于溶蝕的初期，其排泄區(qū)水頭值較高，在右側(cè)邊界處形成了6個泉點排泄口，發(fā)育高程分別為 11，12，13，14，15，16 m，排泄河流高程為10.5 m，泉流量為0.37 mL/s。而在張開度方面，只在河流補給處，表層裂隙部分被溶蝕，張開度變大，深層裂隙沒有大的變化。

1 000 a:排泄區(qū)水位較100 a時下降明顯，下降約4～8 m，泉點變成2個，發(fā)育高程分別為11 m和12 m，排泄河流高程為10.5 m，泉流量為1.96 ×104mL/s，中部水位下降1～2 m，其他地方水位沒有變化。而在張開度方面，排泄區(qū)和中部裂隙進一步被溶蝕，張開度變大，仍在表層范圍內(nèi)，沒有向深層發(fā)展。

2 000 a:較1 000 a的情況，排泄區(qū)水位下降1～2 m，系統(tǒng)中部水位大幅度下降，下降約10 m，泉點仍為2個，發(fā)育高程分別為11 m和12 m，排泄河流高程仍為10.5 m，泉流量為2.27 ×106mL/s。在張開度方面，溶蝕范圍已經(jīng)由表層發(fā)展到了深層，基準面以下的裂隙大部分被溶蝕。

5 000 a:與2 000 a相比，含水層中部水位稍有下降，泉點仍為2個，發(fā)育高程分別為11 m和12 m，排泄河流高程仍為10.5 m，泉流量為3.90 ×108mL/s。而在張開度方面，裂隙擴張的范圍已經(jīng)發(fā)展到河流補給區(qū)。

由上述過程可以看出，2 000 a以后含水系統(tǒng)侵蝕基準面附近的裂隙大部分被溶蝕，這個現(xiàn)象在5 000 a時更加明顯，溶蝕范圍已經(jīng)擴展到河流補給區(qū)。由于2 000 a后泉點始終維持在2個，因此隨著溶蝕過程的不斷發(fā)展，逐漸形成了侵蝕基準面附近的2個大型巖溶管道。由于巖溶管道的存在，進一步的加大了整個系統(tǒng)的溶蝕進程。

為了考察裂隙張開度溶蝕擴展情況，分別做橫剖面圖和縱剖面圖(如圖4)。如圖4(a)(高為11 m的橫剖面)所示，100 a時裂隙張開度幾乎沒有變化，仍為0.02 cm;到了1 000 a時，裂隙張開度略有變化，最大值為1.62 cm;2 000 a時，張開度變化較大，水平方向30～100 m內(nèi)，張開度增大明顯，最大值為9.69 cm，而0～30 m內(nèi)則沒有變化;4 000 a和5 000 a張開度增幅情況與2 000 a時大致相同，大約為 8 cm/ka，最大值分別為25.84 cm和33.91 cm，張開度變化范圍則進一步擴大，為14～100 m。從2 000～5 000 a，張開度的增幅近似恒定，對比圖3中后3幅圖可以發(fā)現(xiàn)，這段時間內(nèi)，裂隙被溶蝕的范圍很大，并且在侵蝕基準面處形成了2個巖溶管道，匯集了大部分的補給水量，由于補給量恒定，因此裂隙張開度的擴張速率也是相同的。

圖4 裂隙張開度演化情況Fig.4 Evolution of the fracture aperture

如圖4(b)(長60 m的縱剖面)所示，1 000 a前裂隙張開度的變化不是很大，而2 000 a后張開度的增幅趨于恒定，為8 cm/ka，但14～30 m范圍內(nèi)的裂隙張開度幾乎沒有變化，原因是2 000 a后，巖溶管道在侵蝕基準面處基本形成，使得這些裂隙都無法被進一步溶蝕變寬。

3.2 非均質(zhì)性巖溶泉域系統(tǒng)演化

本文在均質(zhì)性模型的基礎上，在侵蝕基準面處加上一條張開度為0.08 cm的大裂隙，以此構建非均質(zhì)性模型。

如圖5所示，跟均質(zhì)性系統(tǒng)相比，非均質(zhì)系統(tǒng)在演化初期水位下降較快，100 a時已相當均質(zhì)性系統(tǒng)2 000 a時的水平，但在演化中后期水位幾乎沒有變化，這是因為侵蝕基準面處有一條張開度為0.08 cm的大裂隙，加速了系統(tǒng)的演化，并在1 000 a時發(fā)育成巖溶管道，水流大部分從巖溶管道流走，形成了穩(wěn)定的排泄通道，水位得以相對穩(wěn)定。泉點則由初期的2個變?yōu)橹泻笃谇治g基準面處的1個，演化初期泉點的高程為11 m和12 m，中后期則為12 m，而排泄河流的高程始終為10.5 m。

圖5 非均質(zhì)性巖溶泉域系統(tǒng)演化過程Fig.5 Evolution of a heterogeneous karst spring system

由于大裂隙的存在，差異性溶蝕現(xiàn)象更加明顯，相同時間2種系統(tǒng)裂隙被溶蝕范圍和張開度拓寬速度都有很大的不同。比如，2 000 a時非均質(zhì)系統(tǒng)的溶蝕范圍已延伸到補給區(qū)，比均質(zhì)系統(tǒng)寬近30 m，但在深度方面則相對較淺。從圖6中我們就可以更直觀的看出非均質(zhì)系統(tǒng)中裂隙張開度的變化情況，1 000 a后張開度的增幅速度，與均質(zhì)系統(tǒng)增幅一致，但相同時段內(nèi)張開度的最值卻比均質(zhì)系統(tǒng)要大4～5 cm。

圖6 x=60 m處裂隙張開度縱剖面圖Fig.6 Evolution of the fracture aperture(width vs.aperture)at the 60m length

4 泉流量預測研究

圖7(a)、圖7(b)2個子圖分別是均質(zhì)性及非均質(zhì)性巖溶泉域系統(tǒng)泉流量曲線擬合圖。從子圖7(a)中可以看出，50 a至800 a這段時間，泉流量一直＜9 mL/s，而在900 a時，泉流量有了較大的增長，為1 265.16 mL/s，1 000 a 時增加到了19 590.45 mL/s，1 000 a后泉流量更是以指數(shù)函數(shù)的形式在增加。

產(chǎn)生上述情況的原因是由于巖溶含水系統(tǒng)具有差異性溶蝕現(xiàn)象，而侵蝕基準面具有較強的勢匯，其逐漸開始匯聚較多的補給水流(包括降雨補給和河流補給)，這些水流攜帶有較多的CO2，侵蝕能力強，而且此處水循環(huán)深度淺，源匯距離短，徑流強烈，因此逐漸在侵蝕基準面處形成了大型巖溶管道，不斷地襲奪整個含水系統(tǒng)中的水量，這導致了1 000 a后泉流量以指數(shù)遞增。

而深層原因是在1 000 a時，含水系統(tǒng)發(fā)生穿透現(xiàn)象，裂隙溶蝕已不再局限在表層范圍內(nèi)，而是向深層發(fā)展，張開度迅速增大，有幾條裂隙都是由表層直接穿透到侵蝕面，而且侵蝕基準面處已經(jīng)形成了小型巖溶管道，降雨量及河流補給量沿著貫穿裂隙和巖溶管道迅速排出系統(tǒng)，此時的優(yōu)勢通道最小裂隙張開度為2 cm左右，前人的研究將此作為巖溶系統(tǒng)演化前期和后期的分界線［12］，這就是1 000 a前泉流量增長緩慢，而1 000 a后指數(shù)增長的原因。

圖7 巖溶泉域系統(tǒng)泉流量曲線擬合Fig.7 Fitted curves of spring discharge in the karst spring system

根據(jù)補徑排條件可知，降雨入滲量為泉域主要補給來源，隨著含水系統(tǒng)不斷的溶蝕演化，水位不斷下降，這使得降雨補給深度變短，當發(fā)生穿透現(xiàn)象后，降雨補給量比河流側(cè)向補給量更快沿著貫穿裂隙到達巖溶管道，然后迅速排出系統(tǒng)，而河流側(cè)向補給水量和范圍有限，且影響范圍內(nèi)貫穿裂隙相對較少，因此泉流量主要來自降雨入滲補給量。由于泉水流出后匯入排泄河流內(nèi)，因此排泄河流流量與泉流量成正比，即隨著系統(tǒng)不斷的溶蝕演化，排泄河流流量也隨之增大，但是并不是無限增大，由于補給量一定，其溶蝕能力也是有限的，當侵蝕基準面處的巖溶管道發(fā)育到一定程度必然會停止，排泄河流流量也會趨于穩(wěn)定。

從對泉流量演化機理的分析可以看出，影響泉流量的因素有以下幾點。①補給水量:排泄量的大小與補給量密切相關，在巖溶地區(qū)更是如此，而本文含水系統(tǒng)在1 000 a時，發(fā)生了穿透現(xiàn)象，補給水通過貫穿裂隙和巖溶管道迅速流動，以泉的形式進行排泄，因此其成為了影響泉流量大小的重要因素。②張開度:由立方定律可知，流量與張開度的三次方成正比，1 000 a后，受穿透現(xiàn)象的影響，巖溶管道開始襲奪含水系統(tǒng)中的水流，而本文的補給量是恒定的，因此張開度的大小直接影響泉流量的大小。

根據(jù)泉流量的演化機理、影響因素和泉流量隨時間變化曲線，推算出均質(zhì)性系統(tǒng)和非均質(zhì)性系統(tǒng)泉流量公式:

均質(zhì)

非均質(zhì)

非均質(zhì)性系統(tǒng)泉流量變化曲線形狀與均質(zhì)系統(tǒng)的大體一致，只是其泉流量增幅起始時間提前(由1 000 a提前到了600 a)。相比于均質(zhì)性系統(tǒng)，影響非均質(zhì)系統(tǒng)泉流量大小除了補給量和張開度兩個原因外，還有一個重要原因就是非均質(zhì)性，由于演化初期存在著大裂隙，差異性溶蝕現(xiàn)象明顯，大型巖溶管道形成早，穿透現(xiàn)象也出現(xiàn)得早(600 a)，張開度拓寬速度加快，導致泉流量增幅時間提前，同時段非均質(zhì)系統(tǒng)泉流量比均質(zhì)性系統(tǒng)的要大。

上述式(10)和式(11)只適用于裂隙巖溶介質(zhì)系統(tǒng)，由于本文是在假設條件下對泉流量進行擬合，因此上述公式只是在假設條件下才會出現(xiàn)。本文假設裂隙兩兩正交，相同條件下，裂隙段內(nèi)水流可以更快流動，使泉流量達到最大值，如果裂隙斜交，會使水流流動變慢，泉流量也會減小，因此裂隙相交的角度對泉流量會有一定影響，圖7所示的泉流量在2個系統(tǒng)內(nèi)均是逐漸增大的，但是并不會一直增大下去，本次模擬期為5 000 a，由于補給量是一定的，當巖溶溶蝕演化到一定階段時，裂隙的張開度就會不在變寬，泉流量也會趨于穩(wěn)定不在增長。

5 結(jié)論

本文在耦合滲流模型和溶蝕模型的基礎上，對巖溶泉域系統(tǒng)進行了概化，分別模擬了均質(zhì)性系統(tǒng)和非均質(zhì)性系統(tǒng)的演化過程。均質(zhì)性系統(tǒng)隨著時間的推移，其排泄區(qū)和中部水位逐漸下降，裂隙溶蝕范圍由表層發(fā)展到了深層，在2 000 a后在侵蝕基準面附近形成了2條大型巖溶管道;非均質(zhì)系統(tǒng)與其發(fā)展模式類似，只是差異性溶蝕更加明顯，在1 000 a左右的時間在侵蝕基準面處形成了大型巖溶管道。張開度的增幅速度2個系統(tǒng)大致相同，為8 cm/ka，但是起始時間不同，非均質(zhì)系統(tǒng)比均質(zhì)系統(tǒng)早了1 000 a，而且相同時段的張開度的最大值也要大4～5 cm。

隨著時間的推移，2個含水系統(tǒng)不斷被溶蝕，泉流量也在不斷地增大。由于穿透現(xiàn)象和巖溶管道的出現(xiàn)，使得均質(zhì)性系統(tǒng)1 000 a后泉流量有了大幅度的增加，而非均質(zhì)性系統(tǒng)泉流量拐點出現(xiàn)在600 a。根據(jù)泉流量的數(shù)據(jù)對2個系統(tǒng)的泉流量曲線進行了相應的預測，均質(zhì)系統(tǒng)為指數(shù)函數(shù)，非均質(zhì)系統(tǒng)為多項式函數(shù)。

［1］姜寶良，付北鋒，趙延濤.泉水動態(tài)分析預測和資源評價-以輝縣百泉為例［J］.水文地質(zhì)工程地質(zhì)，2002，(3):45-48.(JIANG Bao-liang，F(xiàn)U Bei-feng，ZHAO Yan-tao.Analysis Prediction and Resources Evaluation of Spring Dynamic:A Case Study of Bai Spring in Hui County［J］.Hydrogeology and Engineering Geology，2002，(3):45-48.(in Chinese))

［2］郝永紅，黃登宇，張文忠，等.山西神頭泉流量的灰色預測模型研究［J］.水利學報，2004，(4):111-114.(HAO Yong-hong，HUANG Deng-yu，ZHANG Wenzhong，et al.Gray Prediction Model for Forecasting Spring Discharge［J］.Journal of Hydraulic Engineering，2004，(4):111-114.(in Chinese))

［3］陳宏峰，朱明秋，夏日元，等.湖南洛塔干河豬場表層巖溶泉BP人工神經(jīng)網(wǎng)絡分析［J］.中國巖溶，2005，24(4):300-304.(CHEN Hong-feng，ZHU Ming-qiu，XIA Ri-yuan，et al.Analysis on Epikarst Spring with BP ANN in Luota，Hunan Province［J］.Carsologica Sinica，2005，24(4):300-304.(in Chinese))

［4］李國敏，楊 福.神頭泉流量衰減的時間序列分析［J］.中國巖溶，1993，12(2):123-130.(LI Guo-min，YANG Fu.Time Series Analysis for the Declining Discharge of Shentou Spring［J］.Carsologica Sinica，1993，12(2):123-130.(in Chinese))

［5］王慶兵，段秀銘，高贊東，等.濟南巖溶泉域地下水流模擬［J］.水文地質(zhì)工程地質(zhì)，2009，(5):53-60.(WANG Qing-bing，DUAN Xiu-ming，GAO Zan-dong，et al.Groundwater Flow Modelling in the Ji'nan Karst Spring Area［J］.Hydrogeology and Engineering Geology，2009，(5):53-60.(in Chinese))

［6］陳南祥，黃 強，曹連海.基于偏最小二乘回歸與神經(jīng)網(wǎng)絡耦合的巖溶泉預報模型［J］.水利學報，2004，(5):68-72.(CHEN Nan-xiang，HUANG Qiang，CAO Lian-hai.Model for Prediction of Karst Spring Flow Based on the Coupling of Neural Network Model with Partial Least Square Method［J］.Journal of Hydraulic Engineering，2004，(5):68-72.(in Chinese))

［7］鞠 琴，郝振純，余鐘波，等.泉流量模擬研究——以小南海泉域為例［J］.四川大學學報(工程科學版)，2011，43(1):77-82.(JU Qin，HAO Zhen-chun，YU Zhongbo，et al.Study of the Model of Spring Flow Forecast:A Case Study of the Xiaonanhai Spring Catchment［J］.Journal of Sichuan University(Engineering Science Edition)，2011，43(1):77-82.(in Chinese))

［8］張有天.裂隙巖體滲流數(shù)學模型研究現(xiàn)況［J］.人民長江，1991，22(3):1-10.(ZHANG You-tian.Current Research in Modeling Seepage Field Fractured Rock Mass［J］.Yangtze River，1991，22(3):1-10.(in Chinese))

［9］王環(huán)玲，余宏明.結(jié)構面模擬分析在巖體質(zhì)量研究中的應用［J］.勘察科學技術，2002，(4):25-29.(WANG Huan-ling，YU Hong-ming.Discontinuity Modeling Analysis in the Application of Rock Mass Quality Research［J］.Site Investigation Science and Technology，2002，(4):25-29.(in Chinese))

［10］WHITE W B.Role of Solution Kinetics in the Development of Karst Aquifers［C］∥IAH.Proceedings of the 12th IAH Congress，Huntsville，Alabama，September 21-27，1977:503-517.

［11］袁道先，劉再華，林玉石，等.中國巖溶動力系統(tǒng)［M］.北京:地質(zhì)出版社，2002.(YUAN Dao-xian，LIU Zaihua，LIN Yu-shi，et al.Karst Dynamic System in China［M］.Beijing:Geological Publishing House，2002.(in Chinese))

［12］GABROV?EK F，DREYBRODTB W.A Model of the Early Evolution of Karst Aquifers in Limestone in the Dimensions of Length and Depth［J］.Journal of Hydrology，2001，240(3/4):206-224.