基于數(shù)值模型的橡膠壩地下水回補(bǔ)作用研究

李 剛，付信偉

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

0 引言

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

1 基于數(shù)值模型的地下水流場預(yù)測

1.1 工程概況及監(jiān)測點(diǎn)設(shè)置

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

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

北外環(huán)橡膠壩設(shè)于主河槽段，壩高3.50m，凈寬70m，共1節(jié)壩袋。由壩室段及上下游連接段組成。上游連接段主要有鋪蓋、岸墻、30m長護(hù)坡組成，鋪蓋為順?biāo)飨?5m長，500mm厚C30鋼筋砼結(jié)構(gòu)，兩側(cè)為C30鋼筋砼岸墻；壩室段為順?biāo)飨?4.0m長C30鋼筋砼底板，厚800mm；下游連接有順流向長22.8mC30鋼筋砼消力池，其中包括7.8m長斜坡段及15m長水平段，消力池深1.4m；消力池下游接順?biāo)飨蜷L15m C25塊石砼海漫，兩側(cè)為C30鋼筋砼岸墻和混凝土聯(lián)鎖塊護(hù)坡50m；左岸灘地設(shè)有泵站，泵站與壩室段經(jīng)設(shè)于壩前鋪蓋的充排水管路系統(tǒng)相連通。

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

1.2 數(shù)值模型構(gòu)建

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

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

(1)

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

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

(2)

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

然后使用二維數(shù)值模型模擬該區(qū)地下水流的時空關(guān)系，模型的表達(dá)式如式(3)所示。

(3)

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

1.3 模型主要參數(shù)設(shè)置

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

表1 各土層主要參數(shù)值

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

(4)

式中，ki—各區(qū)域大氣降水的滲入補(bǔ)給系數(shù)，m3/d；

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

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

表2 北外橡膠壩壩上水位-庫容-水面面積關(guān)系表

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

2 數(shù)值模型應(yīng)用結(jié)果分析

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

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

圖1 模型擬合結(jié)果圖

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)語

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