基于VMOD模型的若爾蓋泥炭沼澤地下水數值模擬

魯瀚友，李志威，2①，胡旭躍，3（.長沙理工大學水利工程學院，湖南長沙 404；2.青海大學省部共建三江源生態和高原農牧業國家重點實驗室，青海西寧 8006；3.水沙科學與水災害防治湖南省重點實驗室，湖南長沙404）

泥炭濕地萎縮現象近年來在國內外普遍發生，造成濕地生態、水土資源面臨威脅［1-2］。若爾蓋泥炭沼澤曾達4 600 km2［3］，是全球最大的內陸高寒泥炭分布區。如今該地區自然溝道和人工溝渠密布，致使地下水位下降、地表水流失、泥炭層退化、有機碳氧化分解，宏觀上表現為泥炭沼澤面積大幅度萎縮，現存面積僅約2 200 km2［4］。影響泥炭濕地萎縮趨勢的重要指標是地形地貌［1］。2011—2017年若爾蓋沼澤濕地的野外調查表明，若爾蓋高原的泥炭地有2種典型的自然溝道：溝道切穿泥炭層和溝道未切穿泥炭層。因此，研究不同溝道對泥炭地地下水運動的影響對于揭示泥炭沼澤萎縮機制具有重要科學意義。

國內學者對若爾蓋高原的地下水變化也有一定關注，如王焱等［5］在朗州公路附近進行地下水的抽水實驗，通過有限差分法計算若爾蓋地下水的平均滲透系數與入滲率；李麗等［6］證實了地下水位降低是若爾蓋有機碳與全氮流失的重要原因；LI等［7］建立In?teractive Groundwater模型研究溝渠排水對濕地退化的影響。對于若爾蓋泥炭沼澤萎縮及地下水疏干，自然溝道的排水作用同樣較顯著［8］。泥炭地含水率高，在自重作用下沿溝道發生垮塌和崩岸，使泥炭地的溝道橫向侵蝕速率高，在一定程度上造成泥炭沼澤面積減少［9］。另一方面，溝道縱向的溯源下切，向泥炭沼澤內部延伸和發展，從而加速封閉泥炭沼澤的地表水及地下水流失［8］。自然溝道對于泥炭地的排水作用既是泥炭濕地萎縮的驅動力，也是當前若爾蓋流域侵蝕過程與物質輸運的集中體現。然而，受降雨頻率與大小的影響，泥炭地的地下水變化較快，波動幅度相對較大，加上缺少地下水的原位監測，使得泥炭地的地下水數據缺失，其地下水運動的定量分析較困難。為此，采用野外泥炭地的地下水原位觀測和地下水數值模擬相結合的方法，是研究泥炭地溝道對泥炭沼澤地下水運動的影響是一個可行的途徑。Visual MODFLOW（VMOD）軟件是一款用于模擬地下水運動的專業軟件，在國際上廣泛應用于地下水及溶質的運移研究［10-11］，如MELIKADZE等［12］用此軟件計算地下水在含水層中的滯留時間；CHETHAN等［13］通過建模模擬降雨對缺水地區地下水的改變。VMOD軟件在國內也有一些應用，如刁維杰等［14］運用此軟件對地下水資源進行了優化配置。該軟件也用于泥炭地及其溝道的數值模擬，如BUJAKOWSKI等［15］運用VMOD軟件研究了泥炭覆蓋對河谷水文地質條件的影響；CRUZ?FUENTES等［16］采用此模型研究了西班牙某干旱地區溝渠灌溉的方案。因此，采用VMOD軟件研究自然溝道對于泥炭地的地下水運動具有可行性。

筆者以若爾蓋高原黑河上游麥曲支流局部泥炭地的地下水野外觀測為基礎，通過野外實測數據，采用VMOD Flex 2015軟件建立切穿泥炭層和未切穿泥炭層2種情況下的泥炭地地下水模型。根據實際情況設計了不同地形坡度，共計20組工況，模擬溝道對泥炭地地下水運動的影響及水力梯度，有助于認識2種典型自然溝道對泥炭地的地下水及濕地萎縮的影響。

1 研究區域與研究方法

1.1 研究區域與野外觀測

研究區位于黃河源支流黑河流域的哈曲上游，屬于紅原縣色地鎮附近的泥炭沼澤濕地。該地區多年平均降雨量為765 mm，地勢低平，降雨匯流于此，排泄不暢形成濕地。1950年代以來，由于人工開挖溝渠和自然溝道的溯源侵蝕，大片濕地被疏干變成草原，沼澤面積與歷史時期相比萎縮近50%。黃河源干流的溯源下切起于180 Myr BP，至今仍保持1～2 mm·a-1的下切速率［17-18］，黃河河道以“U”型彎道切穿若爾蓋盆地，促進若爾蓋黑河、白河流域中自然溝道的橫向侵蝕和溯源下切，從而導致地表水及地下水橫向排水，形成溝道兩側的疏干帶［7］，這一點為實地考察所證實。研究區域選取自然溝道邊緣，觀測點1（32°59′11.130″N，102°59′17.922″E）和觀測點2（32°59′16.212″N，102°59′14.250″E）分別代表自然溝道切穿泥炭層和未切穿泥炭層（圖1）。研究區覆蓋約1.5 m厚的草本質泥炭層，是密集植被根系的淺變質炭化的多孔介質，干密度一般為0.2～0.4 g·cm-3，泥炭的質量含水率約為200%，泥炭下層為粉砂層，50%粒徑（d50）約為0.039 mm。不同的溝道深度導致2種溝道的排水能力及其對地下水運動的影響也不同。

2個觀測點的原位地下水觀測設計相似，均由3排共9個點組成。這些測量孔3排平行于溝道，3排垂直于溝道，且受局部地形影響，每個觀測區9個點間距離各不相同。測壓管由一個1.5 m長的鐵管組成，鐵管外徑2 cm，內徑1.8 cm，在管底部10 cm范圍內有4排穿孔以便地下水通過。每個測量管分別設置于地表以下100、70和40 cm位置，在每個深度插入鋼卷尺進測量管直至管內水面，以此測量地下水的水頭（water head），水頭測量誤差控制在3 mm以內。每測1次，測壓管垂直移動后要等8～17 h，以待水頭達到穩定。通過野外多次觀測，8 h后該地區水頭可達到平衡狀態。測量水頭之后，通過測量鉆孔內地表到地下水的距離，記錄地下水位值。測量方法與FISHER等［19］相似，水位值記錄3次，1次降雨前和2次降雨后。泥炭地表面地形變化較大，使用全站儀測量當地9個位置點的高程。每個觀測區共有27個地下水頭值、9個高程值，以及3個不同時刻共27個水位值。

圖1 研究區域原位測量位置示意Fig.1 Study area and location of two measuring points

除此之外，采用蠕動泵測量了35、72和100 cm的滲透系數K，滲透系數采樣位置在圖1標注的觀測點1與觀測點2，測量滲透系數的方法為CHASON等［20］于1986年基于靜水時差方法提出。先記錄起始水頭，使用蠕動泵抽水，抽水量約為87 cm3，抽水深度約10～15 cm，抽水時間取決于水頭高度，2～4 min不等。抽完水后，連續測量水頭變化直至回升到起始高度，通過水頭變化時間計算滲透系數［20］。

式（1）中，K為滲透系數，m·s-1；r為測壓管內半徑，m；T0為水頭恢復總時間，s。有時候水頭并不能完全恢復，所以測量的恢復時間T0可能會小于實際所需時間?？紤]到這一點，初始階段水位回升遵循線性趨勢，對后面的點進行擬合線性回歸，然后用線性模型預測總恢復時間，帶入式（1）計算滲透系數［21］。各深度滲透系數均測量3遍并取平均值以減少誤差（表1）。

表1 滲透系數測定Table 1 Determination of hydraulic conductivity m·s-1

根據2016—2017年野外實測結果，泥炭層的滲透系數約為 3.89×10-6～3.91×10-5m·s-1，而根據ZHENG等［22］的研究，粉砂滲透系數經驗值范圍為1×10-9～2×10-5m·s-1，故粉砂層滲透系數取值為10-8m·s-1。

1.2 模型建立與工況設計

1.2.1 模型建立

在VMOD模型中導入實測地形點，采用Natural Neighbors差值法生成地形面，通過定義模型結構體功能把地形面生成土體，最后定義土體的各項參數。對于該文討論的物理情景需要建立2種模型進行比較。以觀測點1的實測數據為基礎，建立溝道切穿泥炭層模型，以觀測點2的實測數據建立溝道未切穿泥炭層模型（圖2）。

圖2（a）模型寬9 m，長10 m，溝道寬2 m，深2m，沿x軸正方向坡度為0.5。表層為1.5 m的泥炭沼澤，其下為粉砂的弱透水層，弱透水層是含水層重要組成部分，具有低滲透性和高儲水性，其滲透系數一般小于 10-8m·s-1［22-23］，降雨入滲后水流主要通過泥炭層運動。圖2（b）溝道深度0.5 m且未切穿泥炭層，9個觀測點的位置受地形影響不同。

泥炭層的地下水水面一般在地表附近，該模型研究以潛水流動為主。二維潛水流方程是根據Du?puit假設和質量守恒原理推導的潛水非恒定流動Boussinesq方程。由此基礎上不采用Dupuit假設可得到三維流潛水流動的一般方程為

式（2）中，Kx、Ky和Kz為滲透系數分量，m·s－1；H為水頭，m；Ss為貯水率，m-1；w為源匯項，d-1；t為時間，s。

通過VMOD模型計算20種工況條件，這些工況都是根據實際觀測改變系數實現。實際工況有2種，即切穿泥炭層和未切穿泥炭層情況。這2種工況由于觀測地點不同，井的觀測位置改變（圖3），同樣常水頭與溝道水頭的值也隨實測值發生變化。式（2）給出了地下水潛水流運移方程，但并不能確定具體運動狀態，需要加上定解條件。

圖2 泥炭地的三維數值模型示意Fig.2 Sketch of 3D simulation model in peatland

圖3 邊界條件設定Fig.3 Definition of boundary conditions

1.2.2 邊界條件設定

邊界條件是定解條件的一種。側向常水頭邊界是模型的上邊界，溝道內有水流，屬于模型的下邊界，表達式為

式（3）中，φ1（x，y，z，t）為邊界Γ1上的已知函數，該模型中的邊界均屬于第1類邊界條件。

模型區域的邊界條件包含溝道與常水頭，溝道寬2 m，方向向右。常水頭邊界在坡頂與遠離溝道的側邊，代表地下水來水方向，取值為。坡度方向與溝道同向，2種情況下9個觀測井位置受地形影響各不相同（圖3）。坡腳位置為地下水去水方式之一，溝道切穿泥炭層情況下，地下水位取值在其地表以下，溝道切穿泥炭層情況下，地下水位取值在其地表以下。其余未設置邊界條件的邊緣，默認為0流量邊界［24］。下文研究以模型內9個觀察點數據展開，此小區域內有流量的出入，且距坡腳相對較遠，故能使用以上邊界條件分析瞬時局部水流趨勢。

非穩定流問題需要初始條件，對于地下水模型來水，初始條件指0時刻水頭值，模型中默認的初始水頭為0 m。

地形的影響主要指坡度i改變對地下水的影響，實測點的坡度i=0.5，這里選取以0.01為等距間隔，即從0～0.1的坡度共10組模擬地下水變化情況。

綜合有無溝道、坡度的不同組合，此次共模擬計算了20種工況。每種工況除自變量參數，其他參數與原始值相同，不考慮河岸形狀變化對河岸的影響。此次模擬計算的時間為穩態。有限差分網格劃分為50×45個，并對觀測井附近網格加密1倍。

2 結果與討論

2.1 實際工況下模型模擬結果

圖4為溝道切穿泥炭層觀測點1與未切穿泥炭層觀測點2的VMOD模擬運行結果。圖4（a）中地下水在岸上遞減相對較緩慢，觀測區域水頭每米約下降0.1 m?？拷吰滤^遞減迅速，溝道切穿泥炭層后地下水只能通過邊壁向下匯流或通過弱透水層入滲。這使得溝道邊壁附近地下水的水力梯度達到0.89，等水頭線密集程度是岸上9倍。圖4（b）水頭變化較均勻，觀測區水頭約每米下降0.047 m，溝道未切穿泥炭層時溝道影響相對較小。以圖4（a）中5號點為例，該點約66%的地下水偏向于垂直溝道方向，而圖4（b）中5號點這個數值約為59%，證明了溝道的存在促使泥炭層產生更大的橫向水力梯度，使溝道具有排水的作用，而且溝道切穿泥炭層后排水效果更加明顯。

圖4 地下水水頭等值線圖Fig.4 Contour plot of water head of groundwater in peat layer

圖5是溝道附近某點的來水路徑。泥炭地地下水流運動受溝道與地形影響。地形影響下圖5（a）的跡線有向坡腳運動的趨勢，沿地形坡度方向水力梯度為0.043。在溝道影響下地下水偏向溝道方向運動，沿垂直溝道方向水力梯度為0.084。圖5（b）沿地形坡度方向水力梯度為0.031，沿垂直溝道方向水力梯度為0.044。圖5（a）溝道內水位較低，且并未與泥炭層相連，岸上水頭較高，在邊坡上有一個陡降的過程，地下水位在邊壁下降近1 m。圖5（b）溝道內水位相對較高且與泥炭層連接，跡線較均勻。

圖5 地下水流動軌跡追蹤Fig.5 Path line of groundwater movement

建立模型采用的地形和邊界條件都需簡化且形式不定。數值模型邊緣會選擇特征邊界，如河流、山脊、湖泊等以減少不確定因素，但邊界條件往往是根據需要和對水流運動的判斷建立的，具有一定的不確定性［23］。因此，需要把建立的泥炭地地下水模型的模擬結果與實測值進行對比，以檢驗其可靠性。圖6是觀測值與數值模擬的計算值的對比。圖6（a）的1、2、3、4、6號點觀測值與運算結果相對吻合，5、7、8、9號點相對離散，觀測值大于計算值。分析可能的原因是：（1）受地質、降雨、植物根系等影響，實際情況下溝道附近地下水運動更為復雜；（2）模型是在實際自然條件簡化后建立的，并不能完全模擬地下水的實際流態；（3）泥炭地的地下水水位及水頭是時變的，缺少在時間上高分辨率數據作為輸入條件。對于溝道未切穿泥炭層的情況，圖6（b）的結果類似，但誤差在0.1 m以內，小于圖6（a）。圖6（a）的平均殘差為-0.15 m，均方根誤差為0.21 m；圖6（b）的平均殘差為0.022 m，均方根誤差為0.045 m，說明模擬的計算值較為正確，可反映泥炭層地下水觀測的真實狀況。

圖6 觀測值與模擬值的比較Fig.6 Comparison between observed water head and calculated water head

2.2 地形對地下水的影響

在原始工況的條件下，按表1從0.01～0.10取10個坡度。對研究區9個位置地下水位進行計算，2、5、8號觀測點為垂直于溝道方向斷面（圖7），4、5、6號觀測點為地形坡度方向斷面（圖8）。

地形坡度的變化能引起地下水流動的改變，2、5、8號觀測點是溝道方向的斷面上3個點（圖7）。圖7（a）的3個點坡度每增加0.01，水頭平均下降0.4 cm，2號點與5號點的平均水頭差為7.49 cm，5號點與8號點平均水頭差11.39 cm。這表明沿坡度方向，越往坡腳流向，溝道的地下水水頭下降比例越大，且不受坡度改變的影響。對于溝道未切穿泥炭層，圖7（b）的線間距同步擴大，從0.01坡度時的2 cm，到0.1坡度時擴大至7 cm。說明坡頭到坡腳流向溝道的地下水比例變化不大，但隨坡度的增加地下水向溝道出流加快。

對比圖8（a）垂直溝道方向的3個點，2點間水頭差均為0.07 cm，正如前文描述的溝道切穿泥炭層時岸上遞減緩慢且均勻，在邊壁附近陡降。溝道未切穿泥炭層時，地下水沿垂直溝道方向先快速后緩慢下降，0.05坡度時5～6號點地下水下降速度比4～5號點慢35%左右。

圖7 垂直于溝道方向的坡度-水頭關系曲線Fig.7 Slope‐head relation curve in perpendicular to gully direction

2.3 水力梯度對比

表2是由插值得到的3個點在不同方向上的水力梯度。2號點遠離溝道，8號點靠近溝道，兩者地形坡度方向上的水力梯度之差在溝道切穿泥炭層情況下是0.023 2，溝道未切穿泥炭層情況下0.003 8，前者情況下溝道排水能力更強。垂直溝道上3個點水力梯度值變化不大，切穿時距溝道3 m的點水力梯度均在0.85左右，未切穿時距溝道2.5 m的點則在0.45附近，溝道對等距點的影響相同。5號觀測點位于觀測區域中心，對比5號點橫縱水力梯度得到圖9。切穿溝道泥炭層的5號點在垂直方向的水力梯度偏小，表明其位置地下水相比4號更偏向于地形坡度方向。未切穿溝道泥炭層的5號點在地形坡度方向的水力梯度偏小，也說明其位置地下水相比2號更偏向于垂直溝道方向。

圖8 地形坡度方向坡度水頭關系曲線Fig.8 Slope‐head relation curve along terrain slope direction

表2 典型方向上水力梯度值Table 2 Hydraulic gradient in typical directions

2種典型溝道情況下垂直溝道方向的水力梯度均大于地形坡度方向的水力梯度（圖9）。溝道切穿泥炭層時2個方向水力梯度差值平均為0.04，溝道未切穿泥炭層時這個差值為0.014，表明溝道可加速兩岸沼澤的疏干，而且溝道切穿泥炭層時其作用更強。坡度對水力梯度的影響在2種情況下是不同的，溝道切穿泥炭層時垂直溝道方向的水力梯度在0.01～0.1坡度范圍內只減少0.003，而溝道未切穿泥炭層在同等情況下增加0.052。溝道切穿泥炭層且未與兩岸泥炭層直接接觸，地下水在邊壁附近陡降，通過邊壁流入溝道，地形坡度對水力梯度影響小。溝道未切穿泥炭層時與泥炭層連接，溝道底部以下1 m依舊是泥炭層，地下水的變動受坡度影響較大。

圖9 垂直于溝道方向和沿地形坡度方向的水力梯度對比Fig.9 Comparison of hydraulic gradient in perpendicular to gully and along the terrain slope

結合實測資料，利用VMOD模擬結果的分析表明，溝道附近泥炭地地下水水力梯度偏向溝道方向，證明溝道確實具有疏干地下水的能力。這與LI等［7］的結論一致，認為排水溝渠對濕地退化有重要影響。研究定量分析了2種典型溝道附近垂直溝道方向水力梯度的大小，溝道切穿泥炭層時5號點位置水力梯度是溝道未切穿泥炭層的1.79倍，證明溝道從未切穿泥炭層情況逐漸溯源下切直至切穿泥炭層的過程會加速地下水的流失［8］。通過0.01～0.10地形坡度的對比證明溝道未切穿泥炭層情況受地形影響大，而溝道切穿泥炭層后地下水受地形影響小，這與溝道內水流與泥炭層是否相連有關。

該研究還存在一些不足之處，例如0流量邊界使計算存在誤差，常水頭位置僅代表瞬時水位，并非實際的動水頭邊界。但前人證明VMOD模擬泥炭地地下水有一定可行性［15，25］。這些問題有待增加實測野外數據和深入理論分析以進一步解決。

3 結論

基于若爾蓋局部泥炭沼澤地的地水下實測數據和VMOD模型，建立當地2種典型自然溝道影響下泥炭地的地下水數值模型，對比了若爾蓋2種典型溝道的排水能力。結果表明，濕地的保護可以針對排水作用更強的溝道切穿泥炭層區域進行。主要結論如下：

（1）觀測區超過50%的地下水流入溝道，溝道具有局部疏干排水的能力。若爾蓋泥炭地的自然溝道不僅排走部分地表水，還在非降雨期疏干溝道兩側的地下水，使其兩側的地下水潛水位降低，從而形成泥炭沼澤疏干帶。

（2）自然溝道切穿泥炭層比未切穿泥炭層的距溝道2 m位置垂直于溝道的水力梯度增大79%，即排水能力顯著增強。在溝道溯源下切作用下，未切穿溝道逐漸向切穿溝道發展，疏干帶在溝道兩側放射式擴張。溝道的存在及發育促使若爾蓋泥炭沼澤的地下水不斷流失。

（3）泥炭地的溝道排水能力始終大于小幅度地形坡度影響，即自然溝道的排水功能基本不受局部地形變化而改變。這使得在溝道密布的若爾蓋高原自然溝道排水現象普遍存在，其兩側的疏干帶連結成網絡，從而大范圍降低若爾蓋泥炭濕地的地下水位。