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

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

0 引言

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

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

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

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

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

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

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

(1)

(2)

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

2 SUB源代碼的改進(jìn)

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

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

(3)

本次研究利用有效應(yīng)力變化與水頭變化的關(guān)系將式(3)進(jìn)行進(jìn)一步推導(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)中以?shī)A層內(nèi)中心水頭為自變量，可以很好地耦合到SUB建立的地面沉降模型當(dāng)中，完成對(duì)夾層內(nèi)垂向滲透系數(shù)的調(diào)整。

對(duì)SUB程序包的改進(jìn)主要包括以下3部分：

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

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

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

圖1 SUB改進(jìn)后的邏輯框圖

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

3.1 典型算例建立

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

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

典型算例以美國(guó)地質(zhì)調(diào)查局建立的加州羚羊谷模型為基礎(chǔ)，地層結(jié)構(gòu)參數(shù)和邊界條件與原模型保持一致(表1、2)，對(duì)源匯項(xiàng)進(jìn)行簡(jiǎn)化，只在盆地中心地區(qū)的第二含水層(圖3b)設(shè)置固定開采量(8 500 m3/d)，設(shè)置夾層初始水頭和初始預(yù)固結(jié)水頭與含水層初始水頭一致。針對(duì)上述地質(zhì)結(jié)構(gòu)和條件分別進(jìn)行了變滲透系數(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)水位和沉降量均有減小。下面對(duì)第二含水層(承壓含水層)降落漏斗和沉降漏斗中心網(wǎng)格(23，46)不同結(jié)果進(jìn)行對(duì)比分析。

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

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

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

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

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

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

4)累計(jì)沉降量

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

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

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

4 數(shù)值模擬實(shí)驗(yàn)

4.1 不同開采量對(duì)沉降的影響

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

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

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

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

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

4.2 不同夾層厚度對(duì)沉降的影響

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

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

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

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

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

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

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

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

5 結(jié)論與建議

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

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

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

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