基于MODFLOW-SUB建立變滲透系數的地下水流-地面沉降模型

孟世豪，崔亞莉，田 芳，羅 勇，石鴻蕾

1.中國地質大學(北京)水資源與環境學院，北京 100083 2.北京市水文地質工程地質大隊(北京市地質環境監測總站)，北京 100195

0 引言

由地下水開采引起的地面沉降過程十分復雜，土體的壓縮變形將導致其孔隙比和滲透系數都產生相應的變化，進而影響沉降特征。數值模擬是研究及預測地面沉降演變規律的重要手段[1-7]，為了模擬實際固結沉降過程，不同學者嘗試采用不同的方法來解決變參數下的模擬問題，例如：Joseph[8]利用有限元法討論了變參數的一維固結問題；Chen等[9]采用積分差法刻畫了多層黏土一維非線性固結問題；葉淑君等[10]基于修正的麥欽特模型研發了變參數的地面沉降模擬軟件；劉波等[11]、李莎等[12]利用COMSOL建立了基于Kozeny-Carman方程的變滲透系數沉降模型；Xu等[13]、常嘯[14]基于修正的劍橋模型利用有限元方法對軟基固結進行了分析。從前人研究成果來看，眾多學者均是根據所選擇的方法及具體的地質、水文地質條件自主編程完成對地面沉降壓縮固結的模擬，廣泛推廣的可行性較弱；而COMSOL軟件對地下水補給徑流排泄刻畫較為粗糙，普適性不強。在廣泛應用的MODFLOW軟件中嵌入地面沉降SUB模塊，可用來模擬源于地下水開采造成的區域地下水流-地面沉降過程，但在SUB中設定黏土層的垂向滲透系數為常數，無法真實再現地面沉降隨時間延長的壓縮固結過程。

綜合以上因素，在充分分析前人有關沉降過程中滲透系數變化規律研究成果的基礎上[15-21]，本文結合呂衛清等[18]通過試驗結果擬合描述土體滲透性變化的關系式對MODFLOW-SUB進行改進，建立了變滲透系數的地下水流-地面沉降耦合模型。并以美國地質調查局建立的加州羚羊谷模型[22]為基礎，設計典型算例，在相同條件下，分別以常滲透系數和變滲透系數對該地進行模擬，對比兩者累計沉降量和壓縮釋水量的差異；同時通過數值模擬實驗，設定不同開采量和不同夾層厚度，分析不同條件下夾層內垂向滲透系數變化對沉降差值的影響，以期對改進前后模型的適用性作出評價。

1 基于MODFLOW-SUB建立的地下水流-地面沉降耦合模型

MODFLOW是一套廣泛用于孔隙介質中地下水流動數值模擬的軟件[23]。SUB程序包是MODFLOW中基于太沙基一維垂向有效應力原理的沉降模塊，可模擬由于夾層中地下水的排泄而造成的含水系統的壓縮，也可模擬弱透水層的排水和壓縮的滯后現象。

通過MODFLOW可建立地下水流模型，利用SUB可建立地面沉降模型[24]。在水流模型中，源匯項的改變會導致含水層水頭發生改變，水頭的改變在地面沉降模型中會產生相對應的沉降，在沉降壓縮過程中又會釋放一部分水量，這部分釋水量又會進入水流模型的源匯項中，最終完成地下水流模型與地面沉降模型的耦合。

由于構成夾層的沉積物(一般為黏土、亞黏土、粉土類)滲透系數較低，因此，在含水層系統中夾層的水頭變化通常滯后于周圍含水層的水頭變化。當含水層水頭下降時，夾層內的水頭無法立即與含水層水頭保持一致，此刻，含水層和夾層會產生一個水頭差，通過這個水頭差，從夾層的中間界面向兩邊含水層排水，完成雙面排水過程。

通過對消散夾層與含水層的水頭差過程進行求解[25]，可以得到表示滯后的時間常數τ0，其表示了在給定水頭變化下，所有沉降完成93%時所需要的時間[26]：

(1)

(2)

通過式(1)和式(2)可以發現：夾層內垂向滲透系數和厚度的改變會影響夾層壓縮釋水的過程，影響沉降的滯后性，從而影響夾層的沉降量；且夾層厚度越大，影響程度越大。

2 SUB源代碼的改進

SUB程序包包含5個主要的子程序，分別為數據的讀取、水頭的調整、公式的迭代計算、均衡計算和結果輸出。改進思路是在水頭調整完成之后，讓夾層內垂向滲透系數根據夾層內水頭變化發生相應改變。

針對沉降過程中滲透系數的變化規律，很多學者[13, 28-30]通過試驗總結出了滲透系數與孔隙度、平均孔徑、塑性指數等參數的關系式，但在SUB程序包中，未對夾層的這些參數進行定義，其關系式無法很好地與模型進行耦合。而呂衛清等[18]從固結系數的定義出發，將壓縮系數視為變量，建立了滲透系數與固結應力的微分方程，并通過試驗數據建立了滲透系數與有效應力σ′的關系式：

(3)

本次研究利用有效應力變化與水頭變化的關系將式(3)進行進一步推導，得到夾層內滲透系數與水頭的變化規律(式(4))。

(4)

式中：ρw為水的密度(kg/m3)；g為重力加速度(m/s2);h′ 為當前夾層內水頭(m)；hc為夾層內預固結水頭(m)。

式(4)中以夾層內中心水頭為自變量，可以很好地耦合到SUB建立的地面沉降模型當中，完成對夾層內垂向滲透系數的調整。

對SUB程序包的改進主要包括以下3部分：

1)創建新的數組使得每個網格都獲得單獨的夾層內垂向滲透系數。在SUB中，夾層中的垂向滲透系數是以參數區的形式賦值，即多個網格共用一個垂向滲透系數值。想要使每個網格的夾層內垂向滲透系數能夠隨著該網格水頭的變化而變化，首先要設定一個新的數組使每個網格都獲得單獨的垂向滲透系數。

2)創建新的子程序GWF2SUB7YHK來完成夾層內變滲透系數的功能，并把該子程序插入到參數數據讀取之后、公式迭代計算之前(圖1)。在改進的子程序中讀入每個夾層網格的垂向滲透系數、初始預固結水頭和當前應力期水頭，并按照式(4)求得對應當前應力期水頭的新垂向滲透系數，使得夾層內的滲透系數可以隨著水頭的變化而變化，即完成變滲透系數的調整。

3)在運行改進后的軟件時，需要填寫新的輸入文件。輸入文件所需要內容包括：①改進式(4)中所需的常數參數C，其可以根據不同地區的黏土性質進行計算；②控制輸出某個網格某個應力期夾層內垂向滲透系數的參數。

圖1 SUB改進后的邏輯框圖

3 變滲透系數模型的驗證分析

3.1 典型算例建立

研究區位于美國加州羚羊谷流域(圖2)，羚羊谷是莫哈韋沙漠西部的一個地形封閉盆地，位于洛杉磯東北約80 km處。羚羊谷地下水盆地面積約為2 400 km2。地下水流系統由兩個含水層組組成：第一含水層組為潛水含水層，厚度為150～275 m，主要接受降雨和邊界流入補給；第二含水層組為承壓含水層，厚度為480～600 m，在開采條件下可得到上部潛水的越流補給，與第一含水層組水力聯系密切。模型被離散成一個由43行60列組成的網格(圖3)。模擬的時間跨度為80 a，從1915—1995年，第一年為穩定狀態。

圖2 研究區位置示意圖[22]

典型算例以美國地質調查局建立的加州羚羊谷模型為基礎，地層結構參數和邊界條件與原模型保持一致(表1、2)，對源匯項進行簡化，只在盆地中心地區的第二含水層(圖3b)設置固定開采量(8 500 m3/d)，設置夾層初始水頭和初始預固結水頭與含水層初始水頭一致。針對上述地質結構和條件分別進行了變滲透系數和常滲透系數的模擬。

3.2 模擬結果分析

在原MODFLOW2005(常滲透系數)的模擬下，經過80 a的定量抽水，在抽水區域形成了降落漏斗和沉降漏斗(圖4)。水位降深最大值為60.4 m，沉降量最大值為350.5 mm。

a. 第一含水層；b. 第二含水層及開采范圍。

表2 夾層參數

a. 最終水位降深(1995-12-31T20:15:00)；b. 最終地面沉降量(1995-12-31T20:15:00)。

在變滲透系數的模擬下，含水層水位變化與常滲透系數模擬基本相同，但其夾層內水位和沉降量均有減小。下面對第二含水層(承壓含水層)降落漏斗和沉降漏斗中心網格(23，46)不同結果進行對比分析。

1)夾層內垂向滲透系數

由常滲透系數和變滲透系數模擬下夾層內垂向滲透系數對比分析(圖5a)可以看出：模擬前期(0～20 a)由于滯后夾層中水頭變化不大，因此兩種模擬情況下垂向滲透系數變化都很??；隨著抽水的持續進行，含水層中水頭不斷下降，導致夾層內壓縮釋水不斷進行，夾層中的水頭不斷下降，在20 a后變滲透系數模擬下垂向滲透系數發生顯著變化，其中在20～35 a中夾層內變滲透系數模擬下垂向滲透系數下降速度很快，隨后下降速度減緩，至模擬期結束(80 a)，變滲透系數模擬下垂向滲透系數減小了2個數量級。

2)夾層內壓縮釋水量

在抽水初期，含水層水頭下降快，從第二含水層每個應力期的壓縮釋水量變化對比(圖5b)中可以看出：常滲透系數模擬下，夾層的壓縮釋水量逐漸增大，在第15 a達到最大壓縮釋水量19.5 m3；隨著降落漏斗的增大，含水層水頭下降速度減緩，夾層內壓縮釋水量也逐漸減小。變滲透系數模擬下，壓縮釋水量變化趨勢與改進前一致，隨著垂向滲透系數的減小，夾層釋水過程變得困難；在第20 a，夾層內垂向滲透系數發生明顯差別的同時，其壓縮釋水量也明顯減??；模擬80 a后，累積壓縮釋水量比常滲透系數模型少了393.6 m3，約減少45.1%。

3)夾層內中心水頭

常滲透系數模擬下，夾層內中心水頭由725 m下降至678 m(圖5c)。變滲透系數模擬下，由于夾層內垂向滲透系數的減小，壓縮釋水能力減弱，夾層內的水頭消散速度也相對延緩，夾層內最終水頭為692 m。

4)累計沉降量

隨著夾層內垂向滲透系數的逐漸減小，變滲透系數模型的壓縮釋水量逐漸減少。根據算例的模擬結果來看，變滲透系數模型的最終沉降量(圖5d)也會相對減少，依據本算例中模擬80 a地面沉降量來看，改進后的第二含水層累計沉降量比改進之前的總沉降量減少了大約15.6%。

通過以上對比分析可以發現，兩個模型在20 a內，計算結果具有較好的一致性；而20 a后，兩個模型計算結果差別逐漸顯現出來，說明在這種水文地質條件和開采條件下，由于黏性土夾層壓縮造成的垂向滲透系數的變化滯后20 a左右。

a. 夾層內垂向滲透系數；b. 逐應力期壓縮釋水量；c. 夾層內中心水頭；d. 累計沉降量。

4 數值模擬實驗

4.1 不同開采量對沉降的影響

為探究在不同開采量下夾層內垂向滲透系數變化所導致的沉降差值大小，設定開采量分別為0.85×104、1.42×104、1.98×104、2.55×104m3/d。

由80 a模擬結果顯示，隨著開采量的增加，水頭下降速度和下降幅度都會加大，夾層內垂向滲透系數(圖6)發生明顯差別的時間由第20 年逐步提前至第8 年，最終垂向滲透系數值依次減小。

圖6 研究區第二含水層夾層內垂向滲透系數變化對比

在第二含水層的累計沉降量增大的同時，沉降量累計變化幅度(沉降差值比)也由15.6%逐步增大到28.7%(圖7)。由此可見，隨著地下水開采量加大，兩個模型計算的地面沉降量的差異越來越大，即夾層內垂向滲透系數的變化對沉降差值的影響程度增大。

圖7 研究區不同開采量下第二含水層沉降量對比

4.2 不同夾層厚度對沉降的影響

為了探究在不同夾層厚度下，夾層內垂向滲透系數變化所導致的沉降差值大小，分別設定第二含水層中夾層厚度為0.7、1.4、2.1、2.8、3.5、4.2 m(圖8，9)。其余參數不變，將模擬時間延長至400 a，80～400 a間水位基本保持穩定。

在不同夾層厚度的模擬過程中，其含水層水頭變化相同，但隨著夾層厚度的增加，壓縮釋水的滯后性不斷增大。前80 a水位變化所造成的沉降，0.7 m厚度的夾層需要到第100 年可以全部完成，而4.2 m厚度的夾層需要到400 a左右才全部完成。

通過對比400 a內不同夾層厚度對應的沉降差值比(圖9c)可以發現，當設定夾層厚度為1.4 m時，常滲透系數模型和變滲透系數模型的沉降差值在前25 a內基本保持一致，25 a以后，沉降差值比逐漸增大，在第100 年時完成全部沉降，最終沉降差值比為19.4%；而當夾層厚度擴大3倍，達到4.2 m時，兩模型在第75 年產生差異，最終在第400 年完成全部沉降時沉降差值比達32.4%。

隨著夾層厚度由0.7 m增加至4.2 m，兩模型計算的沉降量產生差異的時間從第17 年逐漸滯后到第75 年。當模擬時間達到400 a時(圖9b)，最終的累計沉降量隨著夾層厚度呈線性增加，同時兩模型的沉降差值比由11.6%上升到32.4%。

由此表明，隨著夾層厚度的增大，兩個模型計算的地面沉降量的差異越來越大，即夾層內垂向滲透系數的變化對沉降差值比的最終影響也會增大。

在眾多的實際應用中，由于監測資料的限制，大多數模型的模擬時間往往在幾十a，不足以完整體現沉降的滯后性。在本算例中，由80 a內不同夾層厚度對應的垂向滲透系數變化對比(圖8)可以看出，夾層厚度的增加導致夾層內垂向滲透系數發生顯著差別的時間逐步延遲。由于夾層厚度增加導致的滯后性，由垂向滲透系數變化所導致的沉降變化幅度(圖9a)先由9.7%上升至17.3%，再下降至1.8%。雖然在第80年夾層厚度為4.2 m的沉降差值比僅1.8%，但在隨后的50 a中，其沉降差值比迅速增長，忽略垂向滲透系數的變化會導致模型的預測和應用產生較大誤差。

圖8 研究區第80 年不同夾層厚度下第二含水層夾層內垂向滲透系數對比

a. 第80年累計沉降量及差值對比；b. 第400年累計沉降量及差值對比；c. 400 a內不同夾層厚度沉降差值比變化對比。

5 結論與建議

1)通過對SUB源代碼的改進，建立起變滲透系數模型。通過典型算例的模擬分析，改進后的變滲透系數模型的最終結果(夾層垂向滲透系數、壓縮釋水量、夾層內中心水頭下降幅度及累計沉降量)都比常滲透系數模型的結果有所減小。

2)在數值模擬實驗中，常滲透系數模型和變滲透系數模型的最終沉降差值比與開采量大小和夾層厚度大小均呈正相關；但隨著開采量的增加，兩模型產生沉降差值的時間逐漸提前，而隨著夾層厚度的增大，其時間逐漸滯后。

3)在典型算例中，地下水開采初期，常滲透系數模型和變滲透系數模型計算的沉降量基本一致，若僅對開采初期進行模擬，選用常滲透系數模型即可。從長時間序列的模擬和預測來看，變滲透系數模型更加符合實際沉降過程，隨著夾層厚度的增大，忽略垂向滲透系數變化所造成的最終沉降差值也逐漸增大，雖然在較短時間內較大厚度的夾層導致沉降差值不明顯，但其會對后期的預測和應用產生較大的影響。

4)本次研究僅對羚羊谷典型算例進行模擬，初步驗證了對MODFLOW-SUB模型改進的合理性，為應用到更復雜的實際場地奠定了基礎。