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

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

1.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083 2.北京市水文地質(zhì)工程地質(zhì)大隊(北京市地質(zhì)環(huán)境監(jiān)測總站)，北京 100195

0 引言

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

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

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

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

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

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

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

(1)

(2)

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

2 SUB源代碼的改進

SUB程序包包含5個主要的子程序，分別為數(shù)據(jù)的讀取、水頭的調(diào)整、公式的迭代計算、均衡計算和結(jié)果輸出。改進思路是在水頭調(diào)整完成之后，讓夾層內(nèi)垂向滲透系數(shù)根據(jù)夾層內(nèi)水頭變化發(fā)生相應(yīng)改變。

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

(3)

本次研究利用有效應(yīng)力變化與水頭變化的關(guān)系將式(3)進行進一步推導(dǎo)，得到夾層內(nèi)滲透系數(shù)與水頭的變化規(guī)律(式(4))。

(4)

式中：ρw為水的密度(kg/m3)；g為重力加速度(m/s2);h′ 為當(dāng)前夾層內(nèi)水頭(m)；hc為夾層內(nèi)預(yù)固結(jié)水頭(m)。

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

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

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

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

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

圖1 SUB改進后的邏輯框圖

3 變滲透系數(shù)模型的驗證分析

3.1 典型算例建立

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

圖2 研究區(qū)位置示意圖[22]

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

3.2 模擬結(jié)果分析

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

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

表2 夾層參數(shù)

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

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

1)夾層內(nèi)垂向滲透系數(shù)

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

2)夾層內(nèi)壓縮釋水量

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

3)夾層內(nèi)中心水頭

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

4)累計沉降量

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

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

a. 夾層內(nèi)垂向滲透系數(shù)；b. 逐應(yīng)力期壓縮釋水量；c. 夾層內(nèi)中心水頭；d. 累計沉降量。

4 數(shù)值模擬實驗

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

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

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

圖6 研究區(qū)第二含水層夾層內(nèi)垂向滲透系數(shù)變化對比

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

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

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

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

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

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

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

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

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

圖8 研究區(qū)第80 年不同夾層厚度下第二含水層夾層內(nèi)垂向滲透系數(shù)對比

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

5 結(jié)論與建議

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

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

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

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