基于數值模型的橡膠壩地下水回補作用研究

李 剛，付信偉

(山東省臨沂市蘭陵縣水利局，山東 蘭陵 277700)

0 引言

地下水與河水的相互補給是自然界水循環的重要環節。人造水源工程在防汛供水、蓄水發電的同時，在地下水-河水循環的環節中具有不容忽視的參與度和影響力，對地下水文和生態環境具有深刻的影響[1-2]。人與水的和諧可持續發展是我國現代水利工程建設的核心精神，意為在滿足發展需求的限度下最大程度地減少對生態環境的破壞[3]。為此需要對人造水源工程的生態影響效果進行評估。作為一種新型的高分子合成材料的人造水壩，橡膠壩的壩高調節更為靈活，理論上對地下水具有較大的回補作用，但具體的影響數值尚需進一步的研究和檢驗[4]。此次研究以蘭陵縣吳坦河橡膠壩水源工程為例，對該區地下水的水位變動數值進行模擬計算，以考察橡膠壩的建設對地下水環境的具體影響。

1 基于數值模型的地下水流場預測

1.1 工程概況及監測點設置

吳坦河屬武河行洪道支流的中運河水系。費縣和蘭陵縣交界處的鳳凰山是其濫觴[5]。現河道全長44.8km，流域面積達510.97km2，其東西支流匯集于吳坦村，由南至北縱貫蘭陵縣城，奔向邳蒼分洪道。吳坦河地處魯南郯蒼平原，地勢北高南低，北部上游為山丘區，出露巖層為寒武系張夏灰巖，以徐莊、毛莊組砂層地層，上部為第四系沖積物所覆蓋，土層厚2～10m。中游為河砂，巖層為奧陶系Q1.Q2地層，上部為第四系沖積物覆蓋。下游為砂礓盤摻雜，黑褐色黏土覆蓋，粘土層厚3m。

北外環橡膠壩工程壩址，位于淮河流域中運河水系吳坦河中游蘭陵縣城區段，吳坦河中泓樁號26+200處。上游260m處為G206國道吳坦河大橋。工程回水長度4.28km，回水至上游小嶺閘。壩址以上流域面積386.43km2，其中包括小馬莊水庫控制流域面積90km2，長新橋水庫控制流域面積36km2。

北外環橡膠壩設于主河槽段，壩高3.50m，凈寬70m，共1節壩袋。由壩室段及上下游連接段組成。上游連接段主要有鋪蓋、岸墻、30m長護坡組成，鋪蓋為順水流向15m長，500mm厚C30鋼筋砼結構，兩側為C30鋼筋砼岸墻；壩室段為順水流向14.0m長C30鋼筋砼底板，厚800mm；下游連接有順流向長22.8mC30鋼筋砼消力池，其中包括7.8m長斜坡段及15m長水平段，消力池深1.4m；消力池下游接順水流向長15m C25塊石砼海漫，兩側為C30鋼筋砼岸墻和混凝土聯鎖塊護坡50m；左岸灘地設有泵站，泵站與壩室段經設于壩前鋪蓋的充排水管路系統相連通。

為監測該區地下水水位的變動，在壩區各壩斷面布設6個地下水監測孔，分為6個區域觀測吳坦河水位變化、大氣降水、其他源匯因素等帶來的地下水位變動，以此作為數值模型的計算數據基礎。觀測孔的孔徑為60mm左右，配設自動檢測儀，每半小時自動記錄觀測數據。

1.2 數值模型構建

壩區淺層地下水主要的補給來源是降水和地表水的滲漏，深層地下水主要由上層徑流補給。自水源工程建立，地表水的滲漏發生改變，地下水的流場和水位不可避免地受其影響[6-7]。為研究工程對該區地下水的具體影響，可以借助構建數值模型預測地下水流場的狀態。這種構建數值模型模擬地下水數值、而后進行計算求解的方法即數值模擬法，是較為通用的地下水模擬方法之一。它有很好的流場分布模擬效果，可以更直觀地呈現地下水的水位變化[8-9]。通過模擬地下水數值，可以對地下水流的水量、水位和水質變化進行研究和分析，為評價和管理地下水資源提供幫助[10]。

在假定滲透系數恒定的條件下，采用有限單元法計算滲流數值，首先計算泛函(H)，計算過程如式(1)所示[11]。

(1)

式中，H—總水頭，m；σ—滲流的實際閾值；kij—以張量形式表示的滲透系數，m/d；i和j—系數矩陣中的行和列；m—流量增量；s1和s2—邊界面。然后如式(2)所示，通過求解等效結點虛流量法計算滲流場的水頭分布和河道的滲流量。

[I]{H}={Q0}+{ΔQ1}-{Q2}

(2)

式中，[I]—計算域中總的傳導矩陣；{H}—未知的結點水頭矩陣；{Q0}—總等效結點的流量矩陣；{Q2}—結點水頭已知的條件下，在虛域中對未知水頭的等效流量矩陣；{ΔQ1}—結點虛流量矩陣，作用于滲流虛區虛單元。

然后使用二維數值模型模擬該區地下水流的時空關系，模型的表達式如式(3)所示。

(3)

式中，k—滲透系數；(x，y)—坐標變量，m；H—地下水位；H0—初始水頭，m；B—含水層地板的標高，m；λ—給水度；i—時間變量，d；D—模擬區域范圍；T1、T2和T3—一類邊界、二類邊界和三類邊界；n—二類邊界的外法方向；q—二類邊界的單位寬度滲流量，m2/d；m′—河床堆積層厚度，m；e—源匯項；Q—河流滲漏補給量，m3/d m2；r—抽水井(注水井)的徑向距離，m。模型求解使用開放性和適用性兼備的Modflow軟件。

1.3 模型主要參數設置

模型應用之前需設置初始參數值，所需設置的水文地質參數包括基質的滲透系數、給水度、含水率孔隙比等。該區主要巖土層各參數的數值見表1。土層主要成分為壤土，混雜少量中粗砂，分布在河床兩岸。壤土的滲透系數為1.5E-05，并且含水率較高，中粗砂的滲透系數為6.8E-03。巖層以溶蝕風化白云質灰巖和裂隙性溶蝕風化白云質灰巖為主，滲透系數分別為1.2E-8～4.6E-9和1.2E-7～7E-10。

表1 各土層主要參數值

設置邊界條件時，如式(4)所示計算大氣降水的滲入補給量。

(4)

式中，ki—各區域大氣降水的滲入補給系數，m3/d；

Ji—各區域大氣降水量。m/d；Si—區域面積，m2。

根據實測壩上河道橫斷面資料量算，壩上水位-蓄水量-水面面積關系見表2。從表中可見，水庫蓄水位、庫容和水面面積三者正向相關，為模型模擬水庫的規模和效益和模擬計算水庫滲流提供依據。

表2 北外橡膠壩壩上水位-庫容-水面面積關系表

然后使用Modflow的River模塊模擬并計算吳坦河地表河水與地下水之間的交換量。運用達西定理分別計算地下水各段的徑流補給水量，該部分數據以注水井注水量或抽水井的抽水量的形式引入數值模型。對于潛水的蒸發量，使用軟件的Evaportranspiration模塊進行模擬計算。該區地下水蒸發的極限深度預估為6m，當超過此界限時，地下水的潛水蒸發量極小，因此忽略補給。蒸發度和蒸發量依據阿維揚諾夫潛水蒸發公式進行計算[12]。對地下水水源的人為開采量數據由歷年統計得出，并作為抽水井的抽水量引入所構建的模型。模型參數設定完成后，使用WHS迭代法進行求解。

2 數值模型應用結果分析

模型求解時，WHS算法的最大外部迭代次數為500，最大內部迭代次數為300，算法的收斂精度為0.001m，殘差收斂精度為0.001m。阻尼系數是1.0。固定各項滲透系數、邊界條件、給水度等參數條件，將模擬的降水、蒸發量、水庫水位、潛水開采量等數據輸入所建模型，并以2021年10月至2022年9月的水位作為均衡期對照，進行模型驗證。驗證后的模型利用所用軟件的Zone Budget模塊均衡分析壩區河水和地下水之間的補給關系，研究橡膠壩的蓄水對該區地下滲流場的作用。模型模擬分析了4種情景下的河水-地下水水位關系。情景一是建壩之后加入河道底部難透水層的情景；情景二是建壩之后不加入河道底部難透水層的情景；情景三是未建壩時加入河道底部難透水層的情景；情景四是未建壩時不加入河道底部難透水層的情景。難透水層雖然滲透性較差，但通過越流也可交換較大的水量，因此作為重要變量納入考量。

2021年10月模擬水位與實測水位擬合如圖1所示，模型模擬的地下水位與實測的水位在南部和西北部部分地區的誤差較為明顯，模型整體擬合誤差較小，可用于進一步的水位預測分析。

圖1 模型擬合結果圖

根據對地下水水位的模擬計算，得到的四種情景下地下水的均衡分析結果。2021年10月—2022年9月，在情景一下的結果見表3，地下水的總補給量為1582.718萬m3，總排泄量為1573.812萬m3，均衡差值為8.906萬m3，總變化量為8.048萬m3。這說明該區地下水流為正均衡狀態，滿足均衡條件。

表3 2021年10月—2022年9月情景一下地下水模擬均衡分析結果 單位：萬m3

在情景二下，地下水模擬均衡分析結果見表4，總地下水得到的總補給量為7856.427萬m3，總排泄量為7852.759萬m3，均衡差值為3.668萬m3，總變化量為2.860萬m3。可見，在建壩之后不加入河道底部難透水層時，地下水一定程度上受到了難透水層的影響。

表4 2021年10月—2022年9月情景二下地下水模擬均衡分析結果 單位：萬m3

在情景三下，地下水模擬均衡分析結果見表5，總地下水得到的總補給量為1026.821萬m3，總排泄量為1020萬m3，均衡差值為6.821萬m3，總變化量為6.140萬m3。從此可知，未建壩時，在模型中加入河道底部難透水層時地下水狀態也處于均衡狀態。

表5 2021年10月—2022年9月情景三下地下水模擬均衡分析結果 單位：萬m3

在情景四下，地下水模擬均衡分析結果見表6，總地下水得到的總補給量為6896.796萬m3，總排泄量為6894.451萬m3，均衡差值為2.345萬m3，總變化量為2.120萬m3。未建壩時若忽略河道底部難透水層，壩區河水與地下的水量補給強度較高，二者發生的補給關系較為密切。

表6 2021年10月—2022年9月情景四下地下水模擬均衡分析結果 單位：萬m3

總的來看，無論在何種情景下，河流補給是該區地下水補給的主要方式，為其貢獻了約過半的補給量。在考慮難透水層的情景下，大壩的建設為該區地下水的水位變動有較大影響。潛水蒸發量與地下水的埋深密切相關，地下水位升高時，蒸發量也會隨之增加。在忽略難透水層時，大壩的建設導致了蒸發的加劇。

在4種情景下，壩區河水與該區地下水的水量交換如圖2所示。由圖可見，除了未建壩時加入河道底部難透水層的情景下秋冬部分月份交換量為負，即發生河流水位較低導致地下水補給河流水的情況，其他情景下交換量均為正。說明在建壩之后無論是否考慮難透水層，河流均對地下水有補給作用，未建壩時不考慮難透水層時河水也對地下水有補給作用。并且，建壩蓄水對該區地下水有較大的水量補給。在情景三下，模擬期的河流為地下水累計補給了36.02萬m3的水量。在情景一下，水壩庫區預計為地下水累計補給水量為177.21萬m3，增加量可達141.19萬m3，增幅預計為392.00%。可見，橡膠壩的建成對增加該區地下水水位、回補地下水資源和維護地下水生態具有積極作用。

圖2 各情景下河水-地下水逐月交換量

3 結語

橡膠壩工程的建設在改變地表河流面貌的同時，可能對地下水生態造成影響問題。為研究工程對地下水水位變動的影響，研究采用數值模擬法構建了數值模型，模擬了蘭陵縣吳坦河橡膠壩水源工程對地下水水位變化的影響。模擬分析的結果顯示，在諸情景下地下水均處于均衡狀態，河流回補了地下水。在考慮河道難透水層的條件下，壩區河流對地下水的預計補給量約為177.21萬m3。可見，該模型擬合和預測性能較好，可以為庫區地下水水資源開發和水生態保護提供數值分析和理論指導，有利于后續庫區水生態評估工作的展開。但研究缺少更充足的觀測數據，還需要增加數據來源以提高模型的可靠性。