李 莎,成建梅,宮輝力

(1. 中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院,湖北 武漢 430074；2.南京大學(xué)地球科學(xué)與工程學(xué)院,江蘇 南京 210023；3. 首都師范大學(xué)資源環(huán)境與旅游學(xué)院,北京 100048)

隨著人類(lèi)社會(huì)的快速發(fā)展,地下水資源被大量開(kāi)采,尤其是在經(jīng)濟(jì)較發(fā)達(dá)地區(qū),地下水的長(zhǎng)期超采導(dǎo)致地面沉降廣泛發(fā)育。作為首都的北京市是我國(guó)地面沉降發(fā)育較為嚴(yán)重的地區(qū)之一,目前仍處于地面沉降快速發(fā)展時(shí)期,地面沉降造成的經(jīng)濟(jì)損失和社會(huì)影響日趨明顯[1～2]。地下水的過(guò)量開(kāi)采是導(dǎo)致北京平原區(qū)地面沉降廣泛發(fā)育的主要原因[3～5]。地下水開(kāi)采引發(fā)地面沉降是一個(gè)十分復(fù)雜的過(guò)程,伴隨著地面沉降的發(fā)生發(fā)展,土體的孔隙率、滲透系數(shù)等參數(shù)也在不斷地發(fā)生變化,這種變化直接影響著土體的滲流和固結(jié)狀態(tài),進(jìn)而反饋在沉降特征中。但目前大部分地面沉降數(shù)值模擬的過(guò)程中為了減少計(jì)算量,通常不考慮滲透系數(shù)的變化[6～7],這種近似在小變形的假設(shè)下是合理的,但在現(xiàn)實(shí)發(fā)育的地面沉降中,累計(jì)沉降量通常達(dá)到了幾百甚至幾千毫米[8～11],此時(shí)忽略其變化所得到的模擬結(jié)果明顯和實(shí)際情況是有差距的。因而很有必要查明滲透固結(jié)過(guò)程中滲透系數(shù)的變化過(guò)程并明確忽略其變化對(duì)沉降過(guò)程及最終沉降量的具體影響,這對(duì)提高地面沉降數(shù)值模擬結(jié)果的精度具有重要意義。

本文在Biot三維固結(jié)模型的基礎(chǔ)上,利用Kozeny-Carman方程對(duì)滲透固結(jié)過(guò)程中滲透系數(shù)的變化規(guī)律進(jìn)行描述,在滲透固結(jié)的過(guò)程中增加考慮了滲透系數(shù)變化的影響因素。并以北京平原的潛水層、第一壓縮層及第一承壓含水層為模擬對(duì)象,在相同的抽水條件下分別進(jìn)行了變滲透系數(shù)和定滲透系數(shù)的模擬工作,主要研究?jī)煞N情況下地面沉降發(fā)展過(guò)程中的差異,對(duì)比二者水頭及沉降量的差距,討論忽略滲透系數(shù)變化引起的誤差以及抽水量、地層參數(shù)等條件對(duì)該誤差的影響。

1 控制方程

數(shù)值模擬是研究及預(yù)測(cè)地面沉降發(fā)生發(fā)展規(guī)律的重要手段[12～15]。1936年Terzaghi和Redulic[16]改進(jìn)了Terzaghi[17]提出的一維固結(jié)理論將其推廣到了二位和三維,但沒(méi)有考慮在土體固結(jié)過(guò)程中總應(yīng)力的變化,因而通常被稱(chēng)為“假三向固結(jié)理論”。Biot[18]提出了著名的Biot固結(jié)理論,這一理論考慮了地下水運(yùn)動(dòng)和土體變形的相互作用,將土體的變形模型和地下水流動(dòng)模型統(tǒng)一于相同的物理空間,真正的將一維固結(jié)理論推廣到了三維空間。在此之后,Biot固結(jié)理論在地面沉降的模擬計(jì)算方面得到了廣泛的應(yīng)用。

1.1 Biot固結(jié)模型

Biot固結(jié)理論是在力學(xué)平衡、有效應(yīng)力原理、土體本構(gòu)關(guān)系、幾何變形關(guān)系和地下水連續(xù)性基礎(chǔ)上建立起來(lái)的,耦合上述5種關(guān)系方程得到Biot固結(jié)模型,其方程組表達(dá)式[18]為：

(1)

式中:γ——土的重度/(N·m-3)；

γw——水的重度/(N·m-3)；

G——土體的剪切模量/Pa；

v——土體的泊松比；

us、vs、ws——x、y、z方向的位移/m；

K——滲透系數(shù)/(m·d-1)；

t——時(shí)間/d；

i——水力坡度；

f——水作用在土顆粒上的滲透力/N。

f=i×γw

其中的未知量是滲透力f和各方向的位移量,它們都是x、y、z坐標(biāo)、水頭和時(shí)間t的函數(shù),在給出一定的初始條件和邊界條件的情況下即可求解。由式(1)可知滲透系數(shù)K將會(huì)直接影響土體形變量。

1.2 滲透系數(shù)變化模型

在多孔介質(zhì)滲透系數(shù)的變化規(guī)律研究及其預(yù)測(cè)方面,1880年Seeheim[19]首次提出了多孔介質(zhì)滲透率與孔隙大小之間的相關(guān)關(guān)系。之后,隨著相關(guān)研究深入,大量的經(jīng)驗(yàn)方程、毛細(xì)管模型以及統(tǒng)計(jì)模型被提出。其中最著名的同時(shí)應(yīng)用最廣泛的是Kozeny-Carman方程,這個(gè)方程是由Kozeny在1927年第一次提出,1937年被Carman[20]修正。之后這個(gè)方程被許多學(xué)者推廣發(fā)展成了多種形式,其中最常見(jiàn)的形式為[21～29]：

(2)

式中：k——多孔介質(zhì)的滲透率/m2；

n——多孔介質(zhì)的孔隙度/%；

C——K-C常數(shù)；

S——多孔介質(zhì)的比表面積/(m-1)；

在式(2)基礎(chǔ)上結(jié)合土體變形前后的體應(yīng)變關(guān)系得[21]：

(3)

式中：k——變形后土體滲透率/m2；

k0——土體初始滲透率/m2；

n0——初始孔隙度；

εv——體應(yīng)變量。

式(3)反映了體應(yīng)變變化對(duì)滲透率的影響關(guān)系。

滲透系數(shù)與滲透率之間的關(guān)系為：

(4)

式中：ρw——水的密度/(kg·m-3)；

g——重力加速度/(m·s-2)；

μ——水的動(dòng)力粘滯系數(shù)/(Pa·s)；

K——滲透系數(shù)/(m·s-1)。

利用式(4)將式(3)加入到Biot固結(jié)模型的方程組中,即可耦合計(jì)算考慮滲透系數(shù)變化情況下的土體變形。

2 北京平原區(qū)模擬背景及模型構(gòu)建

北京平原區(qū)的地層結(jié)構(gòu)是境內(nèi)多條河流的沖洪積作用形成,沿西北—東南方向,由單一地層逐漸過(guò)渡為多層地層結(jié)構(gòu)。地下水主要賦存于河流沖洪積作用形成的砂及砂卵礫石中,一般情況下其成層性較好,易形成較連續(xù)的含水層?；趯?shí)際的地質(zhì)結(jié)構(gòu)和地下水的水力聯(lián)系,結(jié)合地下水的實(shí)際開(kāi)采情況,一般將北京平原第四系地層中的含水層和弱透水層劃分為4個(gè)含水層組。第一含水層組底板埋深25 m左右,屬于潛水,主要接受大氣降水、農(nóng)田灌溉入滲和河流入滲等補(bǔ)給；第二含水層組底板埋深為80～100 m,在開(kāi)采條件下可得到上部潛水的越流補(bǔ)給,地下水流動(dòng)主要以水平徑流為主,排泄主要以農(nóng)業(yè)開(kāi)采為主,與第一含水層組水力聯(lián)系密切；第三含水層組底板埋深150～180 m,主要以水平徑流為主,接受側(cè)向補(bǔ)給和越流補(bǔ)給,開(kāi)采主要用于生活及工業(yè)用水；第四含水層組底板為第四系的底面,埋深一般在200～300 m。第二、三、四含水層均屬于承壓含水層。

本次研究選擇北京平原地區(qū)沖洪積扇中下部的潛水層、第一壓縮層及第一承壓含水層組成的地下水系統(tǒng)為模擬對(duì)象,模型中的土體參數(shù)及含水層厚度均按照實(shí)際情況選取(表1)。在北京平原地區(qū)承壓含水層的開(kāi)采量較大[3, 26],因而選擇第一承壓含水層為開(kāi)采層。

圖1 模型結(jié)構(gòu)圖Fig.1 Model structure

模擬含水層基本結(jié)構(gòu),見(jiàn)圖1,上部第一層為潛水含水層,厚度為20 m,底部為承壓含水層,厚度為20 m,二者之間夾有一層弱透水層,厚度為15 m。模擬區(qū)域是以抽水井為中心,截面半徑為1 000 m,高為55 m的圓柱形地質(zhì)體(圖1),抽水井半徑為0.15 m,貫穿整個(gè)地質(zhì)結(jié)構(gòu),濾水段位于穿過(guò)兩層含水層的井壁處。由于整個(gè)模擬過(guò)程中只有水頭的變化量對(duì)于地面沉降具有實(shí)際意義,因而為了簡(jiǎn)化計(jì)算將地表作為基準(zhǔn)面,即地表高程為0 m。

針對(duì)上述地質(zhì)結(jié)構(gòu)和條件分別進(jìn)行了變滲透系數(shù)和定滲透系數(shù)的模擬。在正式抽水之前首先給定合理的邊界條件計(jì)算出合理的初始流場(chǎng)及應(yīng)力場(chǎng),在此基礎(chǔ)上開(kāi)始抽水模擬,模擬的具體條件見(jiàn)表2。由于第一階段僅僅是為了得到合理的初始條件,并不是本文的研究重點(diǎn),后期模擬結(jié)果分析中的沉降量及水頭變化值不包括第一階段的變化量。本文模型利用COMSOL Multiphysics進(jìn)行求解。

表1 含水層參數(shù)Table 1 Parameters of the aquifer

表2 模擬條件Table 2 Conditions of simulation

3 模擬結(jié)果分析

在抽水開(kāi)始后的一段時(shí)間內(nèi),承壓水水位不斷下降,同時(shí)產(chǎn)生超靜孔隙水壓力,并隨著排水過(guò)程逐漸消散。一段時(shí)間后補(bǔ)給量等于排泄量,各層地下水水位保持不變,此時(shí)超靜孔隙水壓力消散殆盡,在整個(gè)消散過(guò)程中伴隨著土體的固結(jié)變形,超靜孔隙水壓力為零時(shí),地下水系統(tǒng)達(dá)到另一個(gè)穩(wěn)定狀態(tài)。

通過(guò)多次模擬對(duì)比分析發(fā)現(xiàn),抽水作用下的滲透固結(jié)過(guò)程在開(kāi)始抽水后400 h時(shí)已基本完成,此時(shí)各參數(shù)及形變量已達(dá)到變化終點(diǎn),形成穩(wěn)定狀態(tài)。

3.1 滲透系數(shù)變化

由于假設(shè)在模擬過(guò)程中的流體性質(zhì)不發(fā)生變化,所以滲透系數(shù)的變化規(guī)律與滲透率一致。抽水的初始時(shí)刻同一水平面上滲透率相等,在變滲透系數(shù)的情況下抽水初期承壓含水層井壁附近滲透率開(kāi)始降低(圖2),隨著不斷抽水滲透率降低的范圍逐漸擴(kuò)大,但擴(kuò)大的速度逐漸減小,最終趨于穩(wěn)定。

圖2 土體滲透率變化量空間變化圖Fig.2 Spatial variation in soil permeability

滲透系數(shù)的變化是由于抽水過(guò)程中靠近井壁處的水頭首先開(kāi)始降低,土體的有效應(yīng)力增加,靠近井壁的土體首先被壓縮,土體變密實(shí)且孔隙率減小,進(jìn)而導(dǎo)致靠近井壁的土體滲透率首先減小。隨著抽水的進(jìn)行,水頭降低范圍逐漸擴(kuò)大,土體壓縮范圍也增大,所以滲透率減小范圍也隨之延伸,滲透系數(shù)的變化規(guī)律在一定程度上和水頭變化的過(guò)程是一致的。上述考慮土體固結(jié)變形過(guò)程中滲透系數(shù)變化的模擬結(jié)果與實(shí)際情況相符[27]。

圖3分別選取了位于潛水含水層深度10.0 m和弱透水層中部深度27.5 m處過(guò)抽水井的截線,將其上500 h時(shí)滲透率的變化曲線與承壓含水層中深度為45.0 m處截線上的滲透率變化曲線進(jìn)行對(duì)比,發(fā)現(xiàn)潛水含水層和弱透水層在靠近井壁附近的范圍滲透率有所降低,但減少量較小。這是由于在承壓含水層抽水，各層之間的水頭差增大,導(dǎo)致越流補(bǔ)給,造成潛水層及弱透水層內(nèi)的水頭降低,進(jìn)而被壓縮且滲透率減小。綜上所述,隨著抽水的進(jìn)行,地層的滲透率體現(xiàn)出了不均勻變化的現(xiàn)象。

圖3 含水層及弱透水層中滲透率最終變化量對(duì)比圖Fig.3 Comparison of the final variation in permeability of the aquifer and aquitard

3.2 沉降變化

由于抽水作用在地表形成以抽水井為中心的沉降漏斗(圖4)。對(duì)比兩次模擬結(jié)果發(fā)現(xiàn),在沉降漏斗中心部位,考慮滲透系數(shù)變化時(shí)的最終沉降量大于不考慮滲透系數(shù)變化時(shí)的最終沉降量(圖5),前者的最終沉降量是24.126 mm,后者的最終沉降量是23.233 mm,二者相差0.893 mm,這一差值是由于模擬中不考慮滲透固結(jié)過(guò)程中滲透系數(shù)的變化造成的誤差。在沿沉降漏斗中心指向邊界的方向上,二者沉降量之差逐漸減小。這是由于隨著抽水的進(jìn)行,在承壓含水層井壁附近的土體滲透系數(shù)減小,導(dǎo)致在抽水量相同的情況下水頭下降速度變快且穩(wěn)定時(shí)水頭降低幅度較大,即孔隙水壓力減小幅度較大,導(dǎo)致有效應(yīng)力增加量較大,所以土體的最終壓縮量較不考慮滲透系數(shù)變化時(shí)大。在實(shí)際的模擬應(yīng)用中,變滲透系數(shù)的模擬結(jié)果更接近實(shí)際情況,忽略滲透固結(jié)過(guò)程中滲透系數(shù)的變化會(huì)造成得到的最終沉降量偏小,具體偏小的程度與抽水量、土體性質(zhì)以及地層結(jié)構(gòu)有密切的聯(lián)系,因而在精度要求較高的模擬工作中應(yīng)當(dāng)考慮滲透系數(shù)的變化,否則會(huì)造成最終沉降量的錯(cuò)誤估計(jì),可能引發(fā)嚴(yán)重后果。

由圖6可以看出在抽水初期兩種情況下沉降量的差距較小,之后隨著抽水的繼續(xù)二者的沉降量之差逐漸增大,沉降穩(wěn)定后二者之差也達(dá)到穩(wěn)定狀態(tài)。所以在沉降初期忽略滲透系數(shù)變化對(duì)沉降量計(jì)算結(jié)果造成的誤差較小,但有繼續(xù)增大的趨勢(shì),中后期造成的累計(jì)誤差較大并逐漸趨于穩(wěn)定。

圖4 最終形變量三維圖Fig.4 3D diagram of the final deformation

圖5 地表最終形變量對(duì)比圖Fig.5 Comparison of the final deformation

圖6 地表漏斗中心點(diǎn)沉降量對(duì)比圖Fig.6 Comparison of the subsidence in the center point

通過(guò)對(duì)比各地層的沉降量(圖7),發(fā)現(xiàn)沉降量在承壓含水層和潛水層中與埋深線性相關(guān),在弱透水層中呈非線性變化。承壓含水層的壓縮量最大,弱透水層次之,潛水含水層的壓縮量最小,在一定程度上說(shuō)明離開(kāi)采層越近的地層壓縮量越大,忽略滲透系數(shù)造成的沉降誤差也隨著沉降量的增大而增大。

圖7 各地層沉降量對(duì)比圖Fig.7 Comparison of the subsidence of each layer

3.3 各因素對(duì)沉降量及誤差的影響

為了進(jìn)一步研究忽略滲透系數(shù)變化導(dǎo)致的沉降誤差與抽水量、地層參數(shù)及結(jié)構(gòu)之間的具體關(guān)系,在上述模型的基礎(chǔ)上針對(duì)抽水量和承壓含水層的壓縮模量、滲透系數(shù)以及含水層厚度分別進(jìn)行了單因素的變化模擬。

模擬結(jié)果顯示(圖8),隨著抽水量的增加,在最終沉降量增加的同時(shí),忽略滲透系數(shù)變化造成的誤差也隨之增大。隨著承壓含水層壓縮模量和滲透系數(shù)的增加,最終沉降量均有所減小,與此同時(shí)忽略滲透系數(shù)變化造成的誤差也隨之減小。值得注意的是沉降量與承壓含水層厚度之間并不呈現(xiàn)簡(jiǎn)單的正相關(guān)或負(fù)相關(guān)關(guān)系,而是承壓含水層厚度適中時(shí)最終沉降量和誤差都最大,這是由于部分壓縮量較小的弱透水層被部分壓縮量較大的承壓含水層代替后,導(dǎo)致壓縮量增大；但同時(shí)承壓含水層厚度增大也直接導(dǎo)致過(guò)水?dāng)嗝孀兇?因而在外部定水頭條件下使得補(bǔ)給量變大,模擬區(qū)內(nèi)水頭降深減小,進(jìn)而導(dǎo)致沉降量減小,在上述兩種機(jī)理的作用下,最終沉降量隨著含水層厚度的增大呈現(xiàn)先增大后減小的趨勢(shì)。

由對(duì)比可以發(fā)現(xiàn)忽略滲透系數(shù)變化所造成的沉降誤差與最終沉降量呈正相關(guān)關(guān)系。

圖8 不同抽水量、壓縮模量、滲透系數(shù)和承壓含水層厚度條件下沉降量對(duì)比圖Fig.8 Comparison of the subsidence with different pumping conditions, compression modulus, permeability and confined aquifer thickness

4 結(jié)論

利用Biot模型和K-C模型耦合得到變滲透系數(shù)固結(jié)模型,針對(duì)北京平原區(qū)典型含水層結(jié)構(gòu)進(jìn)行抽水引發(fā)地面沉降的模擬,通過(guò)對(duì)比分析得到結(jié)論：

(1)在考慮滲透系數(shù)變化的情況下,在抽水井井壁附近,滲透系數(shù)有所降低,且在遠(yuǎn)離井壁的方向上滲透系數(shù)的下降量逐漸減小,隨著抽水時(shí)間的增加變化量逐漸減小,進(jìn)而趨于穩(wěn)定。隨著抽水的進(jìn)行,地層的滲透率體現(xiàn)了不均勻變化的現(xiàn)象。

(2)考慮滲透系數(shù)變化下的地表最終形變量大于忽略滲透系數(shù)變化的情況,具體差距的大小與抽水量、土體的力學(xué)參數(shù)及含水層結(jié)構(gòu)有密切的關(guān)系,所以在精度要求較高的模擬中不可以輕易忽略滲透系數(shù)在滲透固結(jié)過(guò)程中的變化,否則會(huì)造成較大誤差。

(3)忽略滲透系數(shù)變化導(dǎo)致的模擬誤差在沉降初期較小且具有繼續(xù)發(fā)展的趨勢(shì),中后期較大并逐漸趨于平穩(wěn)。

(4)忽略滲透系數(shù)變化造成的沉降誤差隨著沉降量的增大而增大，且與抽水量呈正相關(guān),與承壓含水層壓縮模量和滲透系數(shù)負(fù)相關(guān),隨著承壓含水層厚度的增加呈現(xiàn)先增大后減小的趨勢(shì)。

參考文獻(xiàn)：

[1] 賈三滿(mǎn),王海剛,羅勇,等. 北京市地面沉降發(fā)展及對(duì)城市建設(shè)的影響[J]. 城市地質(zhì),2007, 1(4): 19-23.[JIA S M, WANG H G, LUO Y,etal. The impacts of land subsidence on city build of Beijing[J]. City Geology, 2007,1(4):19-23. (in Chinese)]

[2] 楊建民,紀(jì)森林. 抽水導(dǎo)致區(qū)域性地面沉降中的s-lnr線性關(guān)系[J]. 巖土工程學(xué)報(bào), 2016, 38(9): 1606-1614. [YANG J M, JI S L. Linear s-lnr relation of regional land subsidence induced by groundwater withdrawal[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(9): 1606-1614. (in Chinese)]

[3] 王祎萍. 北京市超量開(kāi)采地下水引起的地面沉降研究[J]. 勘察科學(xué)技術(shù),2004(5): 46-49. [WANG Y P. Study on land subsidence caused by overexploiting groundwater in Beijing[J]. Site Investigation Science and Technology, 2004(5): 46-49. (in Chinese)]

[4] 楊勇,鄭凡東,劉立才,等. 北京平原區(qū)地下水水位與地面沉降關(guān)系研究[J]. 工程勘察,2013, 41(8): 44-48.[YANG Y, ZHENG F D, LIU L C,etal. Study on the correlation between groundwater level and ground subsidence in Beijing plain areas[J]. Geotechnical Investigation & Surveying, 2013, 41(8): 44-48. (in Chinese)]

[5] 周毅,羅鄖,郭高軒,等. 沖洪積平原地面沉降特征及主控因素——以北京平原為例[J]. 地質(zhì)通報(bào),2016,35(12): 2100-2110. [ZHOU Y, LUO Y, GUO G X,etal. A study of the characteristics of land subsidence and the main control factors in the alluvial plain: a case study of Beijing plain[J]. Geological Bulletin of China, 2016,35(12): 2100-2110. (in Chinese)]

[6] 吳意謙,朱彥鵬. 考慮側(cè)向變形下抽降潛水引起地面沉降的計(jì)算[J]. 華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版),2016,44(4): 116-120.[WU Y Q, ZHU Y P. Calculation of settlement considering lateral deformation due to pumping of phreatic aquifer[J]. J Huazhong Univ of Sci & Tech (Natural Science Edition),2016,44(4): 116-120. (in Chinese)]

[7] 吳意謙,朱彥鵬. 考慮疏干帶非飽和土影響下基坑降水引起地面沉降的計(jì)算[J]. 工程力學(xué), 2016,33(3): 179-187. [WU Y Q, ZHU Y P. Calculation of settlement considering unsaturated soil influence on the dewatering of foundation pits[J]. Engineering Mechanics, 2016,33(3): 179-187. (in Chinese)]

[8] 田芳,郭萌,羅勇,等. 北京地面沉降區(qū)土體變形特征[J]. 中國(guó)地質(zhì), 2012, 39(1): 236-242.[TIAN F, GUO M, LUO Y,etal. The deformation behavior of soil mass in the subsidence area of Beijing[J]. Geology in China, 2012, 39(1): 236-242. (in Chinese)]

[9] 張?chǎng)?宮輝力,陳蓓蓓,等. 北京典型區(qū)地面沉降演化特征與成因分析[J]. 地球信息科學(xué)學(xué)報(bào), 2015, 17(8): 909-916. [ZHANG W, GONG H L, CHEN B B,etal. Evolution and genetic analysis of land subsidence in Beijing typical area[J]. Geo Information Science, 2015, 17(8): 909-916. (in Chinese)]

[10] 姜洪濤. 蘇錫常地區(qū)地面沉降及其若干問(wèn)題探討[J]. 第四紀(jì)研究, 2005, 25(1): 29-33. [JIANG H T. Problems and discussions in the study of land subsidence in the suzhou-wuxi-changzhou area[J]. Quaternary Sciences, 2005, 25(1): 29-33. (in Chinese)]

[11] 胡蓓蓓,姜衍祥,周俊,等. 天津市濱海地區(qū)地面沉降災(zāi)害風(fēng)險(xiǎn)評(píng)估與區(qū)劃[J]. 地理科學(xué), 2008, 28(5): 693-697. [HU B B, JANG Y X, ZHOU J,etal. Assessment and zonation of land subsidence disaster risk of Tianjin Binhai area[J]. Scientia Geographica Sinica, 2008, 28(5): 693-697. (in Chinese)]

[12] 陳崇希,裴順平. 地下水開(kāi)采—地面沉降模型研究[J]. 水文地質(zhì)工程地質(zhì), 2001, 28(2): 5-8.[CHEN C X, PEI S P. Study on model of groundwater exploitation and ground subsidence[J]. Hydrogeology & Engineering Geology, 2001, 28(2): 5-8. (in Chinese)]

[13] 沈孝宇,孫愫文,周?chē)?guó)云,等. 寧波城市地面沉降物理數(shù)學(xué)模型及沉降預(yù)測(cè)[J]. 地球科學(xué), 1989,14(2): 135-144. [SHEN X Y, SUN S W, ZHOU G Y,etal. Physical mathematics model and prediction on subsidence of Ningbo city[J]. Earth Science, 1989,14(2): 135-144. (in Chinese)]

[14] 陳武,裴萬(wàn)勝,李雙洋,等. 飽和多孔介質(zhì)流-固耦合的數(shù)值模型及工程應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報(bào), 2013,32(增刊2): 3346-3354. [CHEN W, PEI W S, LI S Y,etal. Numerical simulation and engineering application of solid-fluid coupled model for saturated porous media[J]. Chinese Journal of Rock Mechanics and Engineering, 2013,32(Sup2): 3346-3354. (in Chinese)]

[15] 金瑋澤,駱祖江,陳興賢,等. 地下水滲流與地面沉降耦合模擬[J]. 地球科學(xué)-中國(guó)地質(zhì)大學(xué)學(xué)報(bào), 2014, 39(5): 611-619. [JIN W Z, LUO Z J, CHEN X X,etal. Coupling simulation of groundwater seepage and land subsidence[J]. Earth Science Journal of China University of Geosciences, 2014, 39(5): 611-619. (in Chinese)]

[16] 安慶軍. 太沙基固結(jié)理論的若干問(wèn)題研究[D]. 南京: 南京工業(yè)大學(xué), 2007.[AN Q J. Study on some problems of Terzaghi consolidation theory[D]. Nanjing: Nanjing University of Technology, 2007. (in Chinese)]

[17] Terzaghi K. Principle of soil mechanics IV: settlement and consolidation of Clay[J]. Engineering News Record, 1925, 5(8): 874-878.

[18] Biot M A. General theory of three-dimensional consolidation [J]. Journal of Applied Physics, 1941, 38(12): 155-164.

[19] Seelheim F. Methoden zur bestimmung der durchl?ssigkeit des bodens[J]. Zeitschrift für Analytische Chemie, 1880, 19(1): 387-402.

[20] Carman P C. Fluid flow through granular beds[J]. Transactions-Institution of Chemical Engineeres, 1937, 15: 150-166.

[21] 劉波,張功,江永華,等. 基于變滲透系數(shù)的深基坑單井抽水沉降研究[J]. 工程地質(zhì)學(xué)報(bào), 2014,22(6): 1123-1127.[LIU B, ZHANG G, JIANG Y H,etal. Settlement research on single pumping well of foundation pit using changeable permeability coefficient model[J]. Journal of Engineering Geology, 2014,22(6): 1123-1127. (in Chinese)]

[22] 劉江濤,廖東良,葛新民. 基于Kozeny-Carman方程的水相相對(duì)滲透率計(jì)算方法[J]. 科學(xué)技術(shù)與工程,2012,12(29): 7500-7504. [LIU J T, LIAO D L, GE X M. Water Phase Relative Permeability Calculation Based on Kozeny-Carman Equation[J]. Science Technology and Engineering, 2012,12(29): 7500-7504. (in Chinese)]

[23] 徐鵬,邱淑霞,姜舟婷,等. 各向同性多孔介質(zhì)中Kozeny-Carman常數(shù)的分形分析[J].重慶大學(xué)學(xué)報(bào), 2011,34(4):78-82. [XU P, QIU S X, JIANG Z T,etal. Fractal analysis of Kozeny-Carman constant in the homogenous porous media[J]. Journal of Chongqing University, 2011,34(4):78-82. (in Chinese)]

[24] 鄭斌,李菊花. 基于 Kozeny-Carman方程的滲透率分形模型[J].天然氣地球科學(xué),2015,26(1):193-198. [ZHENG B, LI J H. A New fractal permeability model for porous media based on Kozeny-Carman equation[J]. Natural Gas Geoscience, 2015, 26(1):193-198. (in Chinese)]

[25] Dvorkin J, Gvirtzman H, Nur A. Kozeny-Carman Relation for a Medium with Tapered Cracks[J]. Geophysical Research Letters, 1991, 18(5):877-880.

[26] Zhu L, Gong H, Li X,etal. Land subsidence due to groundwater withdrawal in the northern Beijing plain, China[J]. Engineering Geology, 2015, 193: 243-255.

[27] 王非,繆林昌,黎春林. 抽水地面沉降中含水層滲流變形耦合模型研究[J]. 鐵道工程學(xué)報(bào),2011(7): 1-5. [WANG F, MIAO L C, LI C L. Study on coupled model for aquifer seepage deformation during pumped ground subsidence[J]. Journal of Railway Engineering Society, 2011(7): 1-5. (in Chinese)]

[28] 熊小鋒,施小清,吳劍鋒,等.彈塑性變形條件下抽水引起的地面沉降三維數(shù)值模擬[J].水文地質(zhì)工程地質(zhì),2017,44(2):151-159. [XIONG X F, SHI X Q, WU J F,etal.3D numerical simulation of elas to-plastic land subsidence induced by groundwater pumping[J].Hydrogeology & Engineering Geology,2017,44(2):151-159. (in Chinese)]

[29] 蘇晨,崔亞莉,邵景力,等.基于理論法的地面沉降數(shù)值模型進(jìn)展[J].水文地質(zhì)工程地質(zhì),2014,41(6):147-152.[SU C, CUI Y L, SHAO J L,etal.Advances in numerical models of land subsidence based on comprehensive theories[J].Hydrogeology & Engineering Geology,2014,41(6):147-152. (in Chinese)]