二元結構平原地區水庫浸沒預測方法研究

王 碧，杜興武，胡 成，羅 垚

(1.中國地質大學(武漢)環境學院，湖北 武漢 430074;2.湖南省水利水電勘測設計研究總院，湖南 長沙 410001)

水庫浸沒是指水庫蓄水后，庫區周圍地下水受庫水頂托作用，排泄受阻，導致地下水水位壅高，使其可能接近或者高于地表而導致該地區土壤次生鹽堿化、沼澤化以及建筑物地基條件惡化等危害[1]。一般來說，在平原地區修建水庫，正常蓄水位一般限制在一級階地附近，也有不少水庫利用兩岸堤防抬高水頭。由于平原地區庫岸地形平坦開闊，地面高程與正常蓄水位相差不大，地下水埋深較淺，易形成浸沒問題，且浸沒面積一般較大，因此對平原地區水庫浸沒范圍進行預測很有必要[2]。浸沒問題是平原地區水庫勘察設計和建設中最重要的水文地質問題，對于常見的沖積成因二元沉積結構的平原地區，我國現行的《水利水電工程地質勘察規范》(GB 50287—99)[3](以下簡稱《規范》)中采用的卡明斯基法在實際工程應用中存在諸多問題。為此，本文分析了利用卡明斯基法對二元結構平原地區水庫浸沒進行預測所存在的不足，并從黏性土的滲透性規律入手，在“啟動水力梯度”的基礎上，提出了適用于承壓含水層模型的折減系數法的計算公式，最后通過實例應用，對比分析了卡明斯基法、折減系數法和數值法3種預測方法對二元結構平原地區水庫浸沒范圍預測結果的適用性和準確性，以為類似工程水庫浸沒問題研究提供參考依據。

1 水庫浸沒的預測方法

1.1 卡明斯基法

對于雙層結構地層，《規范》中采用卡明斯基法進行水庫浸沒范圍預測。

其計算公式如下：

(1)

根據上式計算得出的地下水水頭壅高值H，結合水動力學法可推導出地下水壅高水頭線方程[4]為

(2)

上式中：M為下部含水層(砂礫石)的厚度(m)；K1為下部含水層(砂礫石)的滲透系數(m/d)；K2為上部含水層(黏性土)的滲透系數(m/d)；h1、H為起始斷面蓄水前、后地下水水位至含水層頂板的厚度(m)；h2、y為計算斷面壅水前、后地下水水位至含水層頂板的厚度(m)；Yx為壅水后任意斷面的地下水水位值(m)。

卡明斯基法以上層黏性土底板為計算基準面，水庫浸沒區計算示意圖見圖1。

圖1 水庫浸沒區計算示意圖Fig.1 Calculation diagram of the reservoir immersion zone

利用該方法對二元結構平原地區水庫浸沒進行預測存在以下不足：①卡明斯基法計算模式是地下水補給庫水，水庫蓄水后導致地下水水位整體抬高，但地下水水位仍大于水庫水位，整體的補給關系不變。而平原地區地形平坦、地面高程較小，平原地區修建的水庫水位一般在河流一級階地附近，甚至利用兩岸的堤壩抬高水頭，導致庫水位高于兩岸一級階地。水庫蓄水后庫水位往往高于地下水水位，補給關系發生了改變，由庫水補給地下水，導致地下水水位壅高；②分析卡明斯基法計算公式得出的地下水壅高水頭線(見圖1)可知，離水庫越遠地下水水位越高，其浸沒程度越大，在堤腳處浸沒程度最小。而工程實際情況是，在地形相差不大的情況下堤腳處浸沒程度最大，遠離水庫其浸沒程度逐漸減小；③平原地區二元結構地層中上層為滲透系數較小的黏土層，其厚度較大，而下層為滲透系數較大的砂礫層，兩者級差大于102，黏土層可視為弱透水層或相對隔水層，其阻水作用不能忽略，采用承壓含水層模型計算會更加符合實際。而基于潛水含水層模型的卡明斯基法計算公式顯然已不適用于這種情況，其忽略了上層黏性土的阻水性[5]。

以上三點原因是造成該規范提供的方法不適用于平原地區二元結構地層的主要原因，將會導致預測結果出現偏移，由于上層黏性土的特殊性，應當采取新的方法來預測水庫蓄水后的地下水水位壅高值，使得計算結果更加接近實際值。

1.2 折減系數法

1.2.1 考慮起始水力梯度(I0)

張忠胤[6]研究指出黏土中存在大量的結合水，其對水的流動起到阻礙作用，并提出了“起始水力梯度“的概念，即黏性土中當水力梯度很小時，不存在顯著的滲流，只有當水力梯度達到一定值后才會出現流動的水，這個值就是“起始水力梯度”(I0)。水在黏性土中的滲流運移具有較復雜的形態，其滲透性規律可用下式表示：

v=K(I-I0)

(3)

式中：v為滲透流速；K為黏性土的滲透系數；I為水力比降，即水力梯度；I0為起始水力梯度，即為v-I直線段延長在I軸線的截距，見圖2。

圖2 水在黏性土中的滲透規律曲線Fig.2 Curve of the water permeability in the cohesive soil

由公式(3)可知，當水力梯度I小于I0時，水在黏土中不產生顯著滲流現象；而當水力梯度I大于I0時，v-I曲線呈線性關系，因此I0是表征黏性土滲流特性的重要指標之一。

由于起始水力梯度的存在，承壓含水層中的地下水在頂托作用下逐漸向上覆的弱透水層滲透,在豎直方向上，弱透水層中的水頭損失與穿透厚度成正比例關系，該滲透規律已在部分工程實踐中得到了驗證[7]。如圖3所示，當鉆孔達到a點時，剛好能觀測到穩定的地下水水位，隨著鉆孔加深至b、c點，地下水水位分別上升了ΔH1和ΔH2，當鉆孔揭穿黏土層到達下部含水層時，地下水水位上升至距黏土層底板以上處H0，顯然：

圖3 弱透水層中的含水帶示意圖Fig.3 Diagram of the water-bearing zone in the aquitard

ΔH1/L1=ΔH2/L2=ΔH3/L3

=(ΔH1+ΔH2+ΔH3)/(L1+L2+L3)

=(H0-T)/T=I

(4)

由于各點地下水水位是穩定、靜止的，故v=0，由此得出I=I0。故可推導出任意斷面的黏土層中含水帶厚度T的表達式為

T=H0/(I0+1)

(5)

式中：H0為下層含水層的測壓水位高度；I0為黏土中的起始水力梯度；T為弱透水層中的含水帶厚度，其上界面即為水庫蓄水后地下水水位壅高位置。

1.2.2 承壓含水層模型

平原地區的二元結構地層中由于上層黏性土滲透系數很小且厚度較大，黏性土可視為弱透水層，下層可視為承壓含水層。因此，采用承壓含水層模型計算地下水水位壅高值更加符合實際情況[8]。

由于下層承壓的水頭不是一個定值，隨著距堤壩距離的增加，承壓的水頭呈線性降低，見圖4。由公式(5)可知，黏土層中的含水帶厚度與下層承壓水頭成正比例關系，因此該含水帶厚度也是一個變化值。為了計算出黏土層中含水帶厚度，只需求出下層承壓水頭的變化規律即可。

圖4 二元結構承壓含水層模型示意圖Fig.4 Diagram of the confined aquifer model with dual-formation structure

水庫蓄水前，在堤腳處通過鉆孔鉆穿上層黏土層，獲得砂礫層的承壓水頭H1；然后在計算剖面的最遠端距堤壩Lm處，同樣通過鉆孔獲得H2；利用水動力學中承壓水一維穩定運動模型，可得出單寬流量的計算公式為

(6)

(7)

以上層黏性土底板為計算基準面。將公式(7)代入公式(5)中，可得到黏土層中含水帶厚度的計算公式為

(8)

通過上式即可計算出平原地區二元結構地層中水庫蓄水后上層黏性土中地下水水位的壅高值，從而得出水庫浸沒范圍。該黏土層中含水帶厚度T是個變化值，隨著距堤壩距離x的增加，T減小，相應的水庫浸沒程度降低，符合實際情況。

1.3 數值法

在預測二元結構平原地區水庫浸沒范圍時，因研究區域內地下水水文地質條件較復雜，地下水往往呈空間三維運動狀態，滲流場的時空分布難以用解析方法計算。為了預測在庫水位抬升的影響下庫周地下水水位雍高及其所造成的水庫浸沒范圍，需要在分析研究區水文地質條件的基礎上，建立水文地質概念模型，并采用數值法求解。

利用數值法求解主要有以下4個步驟：一是通過對研究區的水文地質條件進行合理概化后，建立符合研究區實際情況的水文地質概念模型；二是將該水文地質概念模型進行轉換，建立與之對應的地下水滲流數學模型；三是將該數學模型在數值模擬軟件中轉為可視化的地下水滲流數值模型。經過上述轉換過程建立的地下水滲流數值模型是否能全面、客觀地表征研究區實際的水文地質條件和特征，還需要對該模型進行識別與驗證，也稱之為“反演”，即數學運算中的解逆問題，就是利用水頭函數求解地下水均衡方程。由于水頭函數受地下水均衡場內多個水文地質條件控制，因而它是一個多元函數，模型驗證時需對研究區內各種水文地質參數進行調整，具體做法是根據觀測孔的資料，反求水文地質參數；四是將“反演”算出的最優參數值代入地下水滲流數值模型中，計算得到預測結果。

2 實例應用與分析

2.1 研究區概況

湘江長沙綜合樞紐工程位于湘江干流下游的長沙市望城縣境內，壩址選定在望城縣境內的湘江蔡家洲，位于長沙市城區下游約20 km。壩址至上游長沙市區瀏陽河庫段均屬于垸內型的庫區，堤壩內的平原區地勢低緩，地形平坦開闊，是水庫浸沒易發地帶，見圖5。

圖5 蘇托垸地區平面簡圖Fig.5 Plane diagram of the Sutuoyuan area

本文選擇水庫東岸蘇托垸地區為研究區，區內為典型的平原地區二元結構地層，上層為黏土層，下層為砂礫層，主要為第四系地層，其特點如下：

(4) 殘坡積地層(Qedl)：巖性為黃褐色、紅褐色黏土、黏土夾碎石，可硬塑狀，主要分布在低山、丘陵等坡麓地帶。

由鉆孔資料得出黏土層底板高程為25 m，天然情況下，研究區段撈刀河天然水位為25.45 m，達到正常蓄水位29.7 m高度時，按照撈刀河河面的水力坡降值得出河水位為29.76 m。當撈刀河河水位為27.5 m時，堤壩處承壓水位為27.6 m，與河水位基本一致，在遠離堤壩L=2 000 m處的觀測孔內承壓水位為27 m。由于承壓含水層與河水之間水力聯系非常緊密，且相距很近，故兩者水位近似相等。據現場壓水試驗結果，上層弱透水層的滲透系數K2為2.0×10-5cm/s，厚度為5 m，下層強透水層的滲透系數K1為1.5×10-3cm/s，厚度為7 m。

在研究區水庫浸沒評價時，本文根據蘇托垸地區的地形特點，并結合觀測孔布置情況，在遠離堤壩L=2 000 m處選取具有代表性的剖面(見圖5)，利用卡明斯基法、折減系數法和數值法3種方法對水庫浸沒范圍進行了預測。

2.2 卡明斯基法預測

本文利用卡明斯基法計算公式(1)，可計算得出距離河流L=2 000 m處的地下水壅高水頭值H，通過計算得出Lm處的地下水水位值為31.29 m。

將該參數代入公式(2)，化簡得到水庫蓄水后地下水壅高水頭線方程為

(9)

將利用卡明斯基法計算得出的地下水水位線繪制在研究區地形剖面上，并結合前述的水庫浸沒評價標準，繪制臨界地下水水位線，見圖6。

由圖6可見，利用卡明斯基法計算公式預測得到的水庫蓄水后蘇托垸地區地下水水位線接近地表，甚至在平坦及低洼地段溢出地表，而臨界地下水水位幾乎都高于地表，也就是說，水庫蓄水后整個蘇托垸地區均為浸沒區。

圖6 利用卡明斯基法預測得到的蘇托垸地區地下水水位剖面圖Fig.6 Cross section of groundwater level using Kamenski method in Sutuoyuan area

2.3 考慮起始水力梯度后的折減系數法預測

本文通過現場鉆孔試驗及長觀孔數據，得出起始水力梯度I0=0.36，將所有參數代入折減系數法計算公式(8)，可計算得到任意斷面的黏土層中含水帶厚度，從而推導出水庫蓄水后上層黏土層中地下水壅高水位線方程：

(10)

因此水庫蓄水后地下水壅高水頭

Yx=T+25=28.5-2.2×10-4x

(11)

將利用折減系數法計算得出的地下水水位線繪制在研究區地形剖面上，并結合前述的水庫浸沒評價標準，繪制臨界地下水水位線，見圖7。

圖7 利用折減系數法預測得到的蘇托垸地區地下水水位剖面圖Fig.7 Cross section of groundwater level using the reduction factor method in Sutuoyuan area

由圖7可見，距離撈刀河越近，地下水水位壅高值越大，遠離堤壩，地下水水位降低，符合實際情況。此外，距離堤壩越近，水庫浸沒程度也越大，符合常理。

2.4 數值法預測

本文采用有限元數值模擬軟件FEFLOW對蘇托垸地區水庫浸沒區進行數值模擬預測。根據蘇托垸研究區實際條件，可將該地區地下水流概化為非均質各向同性3D水流模型，邊界條件設置如下：東南西三個方向均以撈刀河為界，設為一類邊界類型，邊界上各位置的河水位依據水位站的實測數據以及河水位的平均水力坡降進行估算；北部邊界設為二類流量邊界；上邊界取為自由潛水浸潤面；底邊界取為隔水邊界。大氣降雨和潛水蒸發作為源匯項處理。模擬區在平面上共劃分為44 328個單元、29 175個節點(見圖8)。水文地質參數按照巖性分層賦值。模型運行期為10 a。為了節省篇幅，詳細模擬過程在此省略，將另文闡述。

圖8 蘇托垸地區水庫浸沒區3D空間離散化示意圖Fig.8 3D space discretization diagram of the simulated reservior immersion zone in Sutuoyuan area

本文利用數值法通過數值模擬得到蘇托垸地區水庫蓄水1 a后、5 a后、10 a后的浸沒分區圖(見圖9至圖11)以及水庫蓄水10 a后地下水水位剖面圖(見圖12)。

圖9 蘇托垸地區水庫蓄水1 a后的浸沒分區Fig.9 Distribution of the immersed zone in Sutuoyuan area 1 a after the impoundment of the reservoir

圖10 蘇托垸地區水庫蓄水5 a后的浸沒分區Fig.10 Distribution of the immersed zone in Sutuoyuan area 5 a after the impoundment of the reservoir

圖11 蘇托垸地區水庫蓄水10 a后的浸沒分區Fig.11 Distribution of the immersed area in Sutuoyuan area 10 a after the impoundment of the reservoir

由圖9至圖12可見，水庫蓄水后，研究區內地下水水位明顯雍高，并且部分區域地下水水位略高于撈刀河水位,具體表現如下：水庫蓄水初期，由于河水位的抬升，地表水補給模擬區地下水，致使區域地下水水位整體抬升，逐漸接近地表，產生浸沒災害；在水庫正常蓄水條件下，短期內研究區內出現大范圍的輕微浸沒災害，在近岸段浸沒程度小，而在遠岸段浸沒程度反而更大，主要原因可能是水庫蓄水前期遠岸段地下水水位本身就高，加之區內地下水排泄受阻，因而更早出現浸沒現象；水庫蓄水5 a時，浸沒范圍大大增加，浸沒面積達70%，其中在研究區內西南部的近岸段浸沒程度非常嚴重，東北部浸沒程度最弱;水庫蓄水10a后浸沒范圍基本穩定，浸沒面積約占75%，與水庫蓄水5 a時浸沒情況相比，浸沒面積增加不顯著，主要集中在浸沒程度輕微區。因此，湘江長沙樞紐工程蓄水后，研究區內近岸段浸沒程度最嚴重，遠岸段浸沒程度最輕；西南部基本都處于浸沒狀態，而東北部浸沒程度最輕。

圖12 利用數值法預測得到的蘇托垸地區地下水水位剖面圖Fig.12 Cross section of groundwater level using the numerical method in Shutuoyuan area

3 結論與建議

通過對比3種方法預測得到的蘇托垸地區的地下水水位剖面圖，可以看出：卡明斯基法預測得到的地下水水位線幾乎全部高于地表面，其預測結果是水庫蓄水后蘇托垸地區均為浸沒區；而折減系數法與數值法預測得到的地下水水位線相近，距離撈刀河越近，地下水水位壅高值越大，浸沒程度也越大，遠離堤壩，地下水水位降低，浸沒程度減小，這兩種方法預測得到的結果是蘇托垸部分地區為浸沒區。由此得出以下結論：

(1) 通過上述3種預測方法在具體工程實例中的應用，結果顯示卡明斯基法的預測結果偏大，計算得到的地下水水位值偏高超過約1 m。

(2) 在二元結構平原地區進行水庫浸沒預測時，需要考慮起始水力梯度I0的影響，而考慮起始水力梯度后的折減系數法預測得到的地下水水位線明顯低于卡明斯基法的預測結果，且浸沒水位變化明顯，其預測結果更加合理。

(3)I0的原始定義是“起始水力梯度”，對于這一概念目前尚存爭議，但它是折減系數法的理論基礎。為了驗證I0的存在，應當通過對已建成的水庫工程所引發的浸沒災害情況進行重訪調查，將實際浸沒情況與預測值進行對比，才能得到更加有說服力的證據，并得出更加準確的計算公式或結論。

(4) 解析法(即折減系數法)主要適用于剖面狀地下水水位的計算，而數值法適用于浸沒災害的時空變化預測分析。由于二元結構平原地區地下水滲流規律較為復雜，因此采用數值法對該類地區水庫浸沒進行預測更加準確。

通訊作者：胡 成(1976—)，男，博士，副教授，主要從事工程地質、水文地質和環境水文地質方面的研究。E-mail:hu_cheng@cug.edu.cn