基于GMS-SUB的地面沉降數(shù)值模擬研究
——以太谷東站為例

惠 翔, 劉麗紅

(安徽理工大學(xué)地球與環(huán)境學(xué)院，安徽 淮南 232001)

地面沉降在全世界范圍內(nèi)出現(xiàn)頻繁且影響范圍廣，不僅會造成城市建筑物、道路交通、管道系統(tǒng)線路發(fā)生不均勻沉降而被破壞，還會給排水、防洪、城市內(nèi)澇的防治帶來諸多困難[1-2]。據(jù)統(tǒng)計(jì)，我國自從新中國成立以來至21世紀(jì)初由于地面沉降造成的損失高達(dá)4 500～5 000億元，年均總體損失90～100億元，嚴(yán)重影響經(jīng)濟(jì)社會發(fā)展[3]。我國已有超過17個(gè)省、96個(gè)城市出現(xiàn)明顯的地面沉降，覆蓋面積超過1.0×105km2[4-5]。隨著近年來社會經(jīng)濟(jì)的快速發(fā)展，尤其是在第四紀(jì)沉積層廣泛分布的地區(qū)，大規(guī)模地進(jìn)行地下水抽采是引發(fā)大面積地面沉降的主要原因[6-7]。 目前， 國內(nèi)外關(guān)于沉降的計(jì)算理論和預(yù)測方法，在構(gòu)建形式上主要分為兩類： 第一類是基于研究區(qū)詳細(xì)統(tǒng)計(jì)資料， 通過沉降量與其相關(guān)因素的關(guān)系建立數(shù)理統(tǒng)計(jì)模型和隨機(jī)統(tǒng)計(jì)模型； 第二類是確定性的地下水-地面沉降水土模型[8-11]。

數(shù)理統(tǒng)計(jì)模型一般包括回歸分析模型和灰色理論模型。文獻(xiàn)[12]利用對天津地區(qū)歷史開采水量、歷史沉降量數(shù)據(jù)的整合與分析，建立了地下水位與地面沉降量之間的線性回歸方程，之后對天津地區(qū)未來9a地面沉降趨勢作出預(yù)測。文獻(xiàn)[13]依據(jù)太原市5個(gè)沉降中心觀測點(diǎn)的沉降量建立了灰色Verhulst預(yù)測模型對未來沉降進(jìn)行了預(yù)測。數(shù)理統(tǒng)計(jì)模型和隨機(jī)統(tǒng)計(jì)模型重在考慮地下水開采量與地面沉降量之間的數(shù)量關(guān)系，并沒有考慮土體內(nèi)部水體與土體骨架間的力學(xué)性質(zhì)和沉降內(nèi)在機(jī)理，也沒有考慮地下水流場的動(dòng)態(tài)演化特征，模型不能很好地反映沉降內(nèi)在動(dòng)力學(xué)機(jī)制。

水土模型分為兩步走計(jì)算模型、部分耦合模型、完全耦合模型。1989a上海地質(zhì)勘察院與比利時(shí)合作使用兩步走模型模擬了上海地區(qū)的沉降[14]。兩步走模型是將水流模型與沉降模型分開考慮的，不能充分考慮土體不均一性以及沉降水土耦合情況。文獻(xiàn)[15]在比奧固結(jié)理論的基礎(chǔ)上建立水土全耦合模型模擬地面沉降，這類模型計(jì)算參數(shù)復(fù)雜，計(jì)算耗時(shí)長，模型難以穩(wěn)定并不適合用于大區(qū)域沉降模型計(jì)算。

部分耦合模型可以很好地模擬含水層抽水引起的水位降深，并建立起含水層垂向滲透系數(shù)和貯水率與孔隙率之間的關(guān)系，水位變化時(shí)刻影響著沉降，沉降會導(dǎo)致含水層的壓縮，同時(shí)孔隙水壓力降低含水層水體膨脹從而含水層釋水，釋水量當(dāng)作源匯項(xiàng)又添加到水流方程中，從而完成水土耦合。部分耦合模型能夠模擬多種含水介質(zhì)的變形并且模型涉及的參數(shù)簡潔，計(jì)算明確。目前大多數(shù)抽水型地面沉降模型是基于太沙基有效應(yīng)力原理建立的地下水三維滲流模型與一維垂向固結(jié)沉降的數(shù)學(xué)模型，本文在研究區(qū)的每個(gè)網(wǎng)格剖分點(diǎn)位置都能計(jì)算出隨時(shí)間變化的水位降深，通過水位的降深和水量的變化得出每個(gè)網(wǎng)格點(diǎn)的垂向一維沉降量，每個(gè)點(diǎn)的沉降擴(kuò)展至平面得到平面分布的三維沉降場，從而將地下水動(dòng)力的三維模型與固結(jié)沉降的一維模型進(jìn)行耦合。本文使用GMS-SUB模塊建立地下水滲流與地面沉降部分耦合模型對太谷東站研究區(qū)地下水位和地面沉降的動(dòng)態(tài)變化進(jìn)行模擬，為定量化研究地面沉降問題提供一種可靠的手段。

1 研究區(qū)概況

太焦鐵路太谷東站位于山西省晉中市太谷區(qū)胡村鎮(zhèn)的桑梓村和朝陽村之間，自太谷區(qū)城東北約6km處由東南向西北穿過。車站地理坐標(biāo)為37°40′～38°00′N，112°25′～112°45′E。研究區(qū)位于太原盆地中部偏東，地貌屬沖洪積傾斜平原區(qū)，地形平坦開闊，地勢總體東南高、西北低，地面高程在826.4～788.1m之間，最大相對高差38.8m(見圖1)。

圖1 研究區(qū)位置圖

1.1 地層條件

太原盆地是山西省中部一個(gè)大型新生代斷陷盆地，內(nèi)部沉積有巨厚的新生界地層。太谷東站研究區(qū)有復(fù)雜的第四系松散沉積層，通過地質(zhì)鉆孔揭露，深度達(dá)230m。第四系地層在研究區(qū)內(nèi)廣泛分布，總厚160～175m，區(qū)內(nèi)自上而下依次分布第四系下更新統(tǒng)河湖相沉積物、中更新統(tǒng)沖洪積物、上更新統(tǒng)沖洪積物、全新統(tǒng)沖洪積物，本次研究重點(diǎn)層位為第四系。

1.2 地質(zhì)構(gòu)造

研究區(qū)位于太原盆地東側(cè)侯城-平遙陷隆內(nèi)，基底總體上是一個(gè)復(fù)向斜構(gòu)造，軸向北東。距離東陽-凈化斷裂約2km，該斷裂是控制盆地東部的一條主要斷層，其總體延伸方向?yàn)楸睎|-南西向，屬壓性兼扭性斷裂。

1.3 水文地質(zhì)特征

第四系松散沉積層自上而下可以分為四個(gè)含水巖組。第一層為全新統(tǒng)潛水-微承壓含水層，由于垂向上含水層無較好的隔水性能，不具有承壓水性質(zhì)，地下水水位呈自由水面，在局部地帶含水層與黏性土分布連續(xù)且穩(wěn)定，微具有承壓水性質(zhì)，現(xiàn)已疏干，富水性極弱，為透水不含水的透漏水層。第二層為上更新統(tǒng)第一承壓含水層，由于淺層水超采嚴(yán)重，導(dǎo)致淺層水位持續(xù)下降出現(xiàn)疏干含水層狀況，形成區(qū)域性超采漏斗區(qū)。第三層為中更新統(tǒng)第二承壓含水層，此層水為中層水，地下水以上部越流及降水入滲補(bǔ)給為主，徑流條件好，交替循環(huán)性好，但隨著地下水集中開采，人為改變了地下水流場。第四層為下更新統(tǒng)第三承壓含水層，為深層水，地下水補(bǔ)給主要為其他含水層越流補(bǔ)給。承壓含水層之間弱透水層相對于含水層厚度較薄，并且很多位置土層不連續(xù)存在缺失。中深層水是研究區(qū)混合集中開采層位，所以模擬的重點(diǎn)為第二、第三承壓含水層中深層地下水。

2 太谷東站地面沉降數(shù)值模擬模型

2.1 概念模型的建立

根據(jù)提供的水文地質(zhì)資料，潛水-微承壓含水巖組和第一承壓含水巖組均已經(jīng)達(dá)到疏干狀態(tài)，因此目前研究區(qū)主要開采中深層地下水，數(shù)值模擬主要模擬的含水層為第二承壓含水層和第三承壓含水層。模擬研究區(qū)面積約24.3km2，根據(jù)水文地質(zhì)條件概化結(jié)果對研究區(qū)進(jìn)行網(wǎng)格剖分，垂向上剖分為2層，平面上將單元格剖分為200行200列，共計(jì)80 000個(gè)單元格，平面和垂直剖分圖如圖2所示。研究區(qū)外側(cè)為無效單元格，其他網(wǎng)格為模型計(jì)算范圍單元格，有效單元格數(shù)為41 462個(gè)，無效單元格數(shù)為38 538個(gè)。

圖2 研究區(qū)網(wǎng)格剖分圖

2.2 地下水滲流與地面沉降耦合數(shù)學(xué)模型

1)地下水三維滲流數(shù)學(xué)模型 將研究區(qū)地下水流概化成非均質(zhì)非穩(wěn)定各向異性三維地下水流系統(tǒng)， 并建立相應(yīng)的數(shù)學(xué)模型， 該運(yùn)動(dòng)微分方程如下

(1)

式中：kx、ky、kz為滲透系數(shù)在x、y、z方向的分量(假定滲透系數(shù)主軸方向與坐標(biāo)軸的方向一致m·d-1)；w為單位體積流量，用以代表流進(jìn)源或流出匯的水量；μs含水層的儲水率，m-1；D為模擬區(qū)范圍；Γ1為第一類給定水頭邊界條件；Γ2為第二類給定流量邊界條件；q為邊界流量。

2)地層垂向一維變形模型 在一維垂直方向上，地面沉降遵循有效應(yīng)力原理。在某一地下水位水平面上，由于大量開采地下水導(dǎo)致含水層孔隙水壓力減小，含水層有效應(yīng)力增大，從而引起含水層骨架的壓縮。有效應(yīng)力的變化能夠引起組成含水層的沉積層壓縮及膨脹，沉積層的壓縮主要是含水層弱透水夾層厚度的減少，一維可壓縮系數(shù)α可定義為

(2)

式中：db為壓縮層厚度的變化量，m3；dσ有效應(yīng)力變化量，Pa。

如果有效應(yīng)力的變化僅僅是由于孔隙水壓力的變化，則可以

(3)

式中：Ssk=ρwgα為骨架儲水率，m-1；Sk=Sskb為骨架儲水系數(shù)，無量綱；dh為水頭變化量，m。

則壓縮量可進(jìn)一步表達(dá)為

Δb=SkΔh

(4)

式中：Δb為夾層厚度變化量，m；Sk為骨架儲水系數(shù)；Δh為水頭變化量，m。

3)水流模型與土體模型耦合 三維地下水流的運(yùn)動(dòng)方程為

(5)

式中：Kxx、Kyy、Kzz分別為x、y、z方向上的滲透系數(shù)，m·d-1；h為水頭，m；W為源匯項(xiàng)，m·d-1；Ss為儲水率，m-1?？蓧嚎s夾層的水量變化可以添加到等式右邊，同樣也可以并入到源匯項(xiàng)W中。由于含水層抽水導(dǎo)致弱透水夾層水頭降低，有效應(yīng)力增大導(dǎo)致地面沉降，之后夾層壓縮會將夾層內(nèi)的水排出，進(jìn)入地下水流系統(tǒng)中，同時(shí)沉降方程與水流方程通過水頭項(xiàng)進(jìn)行耦合，由水流方程求解出水頭變化值再代入沉降方程計(jì)算變形值，進(jìn)而完成部分耦合。

2.3 地下水開采

根據(jù)已有地下水開采資料，研究區(qū)對中深層含水層進(jìn)行抽水，模擬識別驗(yàn)證期概化共計(jì)31個(gè)抽水井，單井抽水量為1 125m3/d。

2.4 各含水層參數(shù)的確定

本研究區(qū)范圍面積不大，水平滲透系數(shù)不存在較大波動(dòng)，根據(jù)對當(dāng)?shù)赝临|(zhì)的分析，確定了大致含水層水平滲透系數(shù)范圍，同時(shí)經(jīng)過不斷地調(diào)整參數(shù)來調(diào)整模型。根據(jù)調(diào)參過程不斷擬合模型大致分為兩個(gè)參數(shù)分區(qū)。通過水文資料的收集以及模型調(diào)參驗(yàn)證(見圖3)，各含水層參數(shù)如表1所示。

表1 含水層參數(shù)取值表

圖3 研究區(qū)水平滲透系數(shù)分區(qū)圖

2.5 邊界條件

將模型東部西部邊界作為隔水邊界，南北部邊界概化為一般水頭邊界。

2.6 夾層參數(shù)的確定

利用研究區(qū)鉆孔剖面柱狀圖將第二承壓含水組劃分為3個(gè)壓縮沉降層，將第三承壓含水組劃分為3個(gè)沉降壓縮層。夾層彈性儲水系數(shù)Sfe由彈性儲水率Sske乘以夾層厚度得到，夾層非彈性儲水系數(shù)Sfv由非彈性儲水率Sskv乘以夾層厚度得到。夾層的非彈性儲水率Sskv由土工試驗(yàn)測得的原位土壓縮系數(shù)和式確定。

(6)

式中：av為壓縮系數(shù)，kPa-1；e0為初始孔隙比；rw為水的重度，kN·m-3。

Sske取比Sskv少一個(gè)數(shù)量級(見表2)。

表2 非滯后夾層計(jì)算參數(shù)取值表

滯后夾層是用于表示在模擬中水頭差消散時(shí)間顯著大于時(shí)間步長的夾層。由于任何含水層都可能包含大量不同厚度的夾層，為減少計(jì)算量，可以將同一模型層內(nèi)具有相同垂直滲透系數(shù)和彈性與非彈性骨架儲水率的滯后夾層合為一個(gè)滯后夾層系統(tǒng)。垂直水力擴(kuò)散率相似的n個(gè)滯后夾層系統(tǒng)的等效厚度bequiv可計(jì)算為

(7)

為了保證相同的夾層總量，為滯后夾層系統(tǒng)提供正確的壓縮量，需要將壓縮量和與周圍含水層交換的水量乘以該比例因子nequiv

(8)

可以計(jì)算出夾層系統(tǒng)的總壓縮量的大小和時(shí)間歷程。因此，方程只求解了一個(gè)等效厚度bequiv的夾層，計(jì)算總的夾層壓縮量需要乘以nequiv。滯后夾層各項(xiàng)參數(shù)如表3所示，先期固結(jié)水頭取模型識別驗(yàn)證期初始水頭。

表3 滯后夾層計(jì)算參數(shù)取值表

2.7 時(shí)間離散

模擬的識別驗(yàn)證期為2015a 8月～2020a 8月，識別驗(yàn)證階段是根據(jù)提供的水文地質(zhì)條件，研究區(qū)地下水位持續(xù)下降，局部降落漏斗逐步形成階段。模型一共分為60個(gè)應(yīng)力期，每個(gè)應(yīng)力期分為3個(gè)時(shí)間步長，總計(jì)180個(gè)時(shí)間步長。

2.8 模型的識別驗(yàn)證

本研究通過山西省第三勘察院多年來對晉中盆地地下水動(dòng)態(tài)監(jiān)測成果資料，獲得1995a至今地下水位下降速率約為1.11m/a。由模擬結(jié)果可以看出，在2016a 8月第二承壓含水層已經(jīng)有部分地區(qū)水位到達(dá)含水層頂板了，中深層水降落漏斗已經(jīng)初步形成，降落漏斗位于太谷東站向東1公里處的桑梓村，漏斗面積可達(dá)3.26km2，漏斗中心水位埋深48.7m。到2020a 8月研究區(qū)水位進(jìn)一步下降，降落漏斗也逐步向外擴(kuò)大，漏斗面積達(dá)到4.16km2，漏斗中心水位埋深53.3m。模擬區(qū)地下水流場與2020a 8月中深層等水位線基本吻合，水位變化趨勢不斷接近真實(shí)觀測水平，參數(shù)率定達(dá)到最優(yōu)，得到了較為理想的模型模擬識別驗(yàn)證結(jié)果，此模型可以用來模擬太谷東站研究區(qū)地下水流場分布，為之后與地面沉降模型的耦合提供了保障。

中深層地下水降落漏斗中心與地面沉降中心基本一致，且降落漏斗形狀擴(kuò)展范圍與沉降中心范圍也基本相似，兩者之間存在高度相關(guān)性。根據(jù)研究區(qū)INSAR衛(wèi)星對研究區(qū)2019a 8月～2020a 8月地面沉降監(jiān)測結(jié)果來看，太谷東站位于一個(gè)沉降漏斗中心，沉降量范圍在20～26mm， 沉降漏斗范圍9.36km2。將INSAR監(jiān)測沉降線與模型運(yùn)行結(jié)果比較，模型整體擬合了沉降趨勢[16-17]。綜上所述，該模型可以用來預(yù)測研究區(qū)地面沉降(見圖4～圖5)。

圖4 含水層地下水位圖(單位：m)

圖5 2019a 8月～2020a 8月地面沉降量等值線圖(單位：m)

3 谷東站地下水位與地面沉降預(yù)測

利用驗(yàn)證后的GMS-SUB水土耦合模型，以2020a 8月作為預(yù)測的初始時(shí)刻，模擬2025a 8月的地下水水位變化和地面沉降趨勢?，F(xiàn)取6個(gè)沉降預(yù)測點(diǎn)：太谷東站點(diǎn)、東站東北方向500m處的中藥原料加工廠點(diǎn)、東北方向1.8km處的董村點(diǎn)、西南方向1.3km處的朝陽中學(xué)點(diǎn)，正東方向1km處的桑梓村點(diǎn)、東南方向1.6km處的沙溝西閣點(diǎn)。 6個(gè)預(yù)測點(diǎn)水位降深分別為8.1m、7.2m、 4.2m、 2.6m、7.5m和6.9m(見圖6)。

圖6 各沉降預(yù)測點(diǎn)水位預(yù)測圖

采用經(jīng)驗(yàn)公式法計(jì)算預(yù)測點(diǎn)的沉降量。各地層壓縮量公式為

ΔSi=MiHiΔp

(9)

式中：ΔSi為不同巖性地層的壓縮量，m；Mi分別為砂土、粘土壓縮率，m-1；Hi為砂土、粘土的層厚，m；Δp為水位降深，m。

2025a 8月地下水位預(yù)測等值線和地面沉降預(yù)測等值線如圖7～圖8所示。相同水位變幅下經(jīng)驗(yàn)公式法與GMS-SUB計(jì)算的累積沉降量如表4所示，經(jīng)驗(yàn)公式與GMS-SUB數(shù)值模擬預(yù)測的沉降量值較為接近，建立的模型可靠，同時(shí)隨著地下水的逐年開采，沉降漏斗范圍有繼續(xù)擴(kuò)大的趨勢。

表4 相同水位降深下經(jīng)驗(yàn)公式與GMS-SUB計(jì)算沉降量

圖7 含水層預(yù)測水位線(單位：m)

圖8 地面沉降預(yù)測等值線(單位：m)

4 結(jié)論

(1)依據(jù)太沙基有效應(yīng)力原理和三維滲流理論，建立了太谷東站地下水三維流和土體一維垂向變形耦合模型，可以刻畫出研究區(qū)實(shí)際水文地質(zhì)特征。

(2)利用水土部分耦合模型對研究區(qū)地面沉降進(jìn)行研究，模型在現(xiàn)了地下水水位與地面沉降的動(dòng)態(tài)特征與演化趨勢，并預(yù)測了2025a 8月研究區(qū)的沉降。太谷東站在預(yù)測期累積地面沉降量最高達(dá)到51.2mm，沉降漏斗面積可達(dá)10.64km2。GMS-SUB模塊計(jì)算地面沉降參數(shù)簡明，速度快，實(shí)用性高，計(jì)算結(jié)果能反映地面沉降動(dòng)態(tài)過程，為地面沉降防治工作提供科學(xué)依據(jù)。

(3)基于GMS-SUB模塊的太谷東站地面沉降模型主要考慮的是一維垂向沉降，而實(shí)際的地面沉降往往也伴有水平向的變形，水土模型數(shù)值模擬還需要進(jìn)一步完善。