巖溶塌陷水動力-力學耦合過程數值模擬
——以武漢市青菱鄉為例

陳冬琴，唐仲華

(1.中國地質大學(武漢)環境學院，武漢 430074；2.湖北城市建設職業技術學院，武漢 430000)

1 研究區巖溶塌陷區地質條件

烽火村喬木灣巖溶塌陷位于武漢市洪山區青菱鄉烽火村二、三組，上、下倒口湖之間，東距京廣鐵路750 m，西距107國道120 m，距長江直線距離1.8 km。1997-2005年先后共發生大小陷坑23處，且近幾年有新的塌陷坑形成，陷坑平面形態呈橢圓形，烽火村塌陷坑直徑一般<10 m，深跨比≥0.3，認為是相互連通的擾動帶和多個土洞發育而引起的塌陷。

垂向上分成三個發育段。

(1)上發育層。溶洞埋深20～40 m，標高段一般為(0～20 m)，鉆孔遇洞率15.71%，此埋深段為研究區溶洞垂向發育程度最強。幾乎存在于巖溶強發育區及較強發育地區，分布較廣泛，面狀多發育為溶孔、溶槽的形式，此層溶洞發育規模一般都比較大，最大溶洞高度為7.8 m。

上層發育帶分為四段：20～25，25～30，30～35，35～40 m；各埋深段的巖溶發育情況如圖1所示。從圖中可以看出20～25 m埋深段是一個弱發育段，25～30 m為巖溶最發育段。

圖1 巖溶發育程度與深度關系Fig.1 Degree of karst development vs. the depth

(2)中發育層。埋深在 40～65 m，標高段一般為-20～45 m， 該段巖層溶洞較發育，所有的石灰巖區普遍存在，遍布青菱鄉、毛坦港、阮家巷、烽火村等地區，溶洞規模不大，洞高一般在3 m以下，全部分布在塌陷區，表明地表塌陷較深層次的溶洞也起著至關重要的作用。

但50～55 m埋深段、60～65 m埋深段鉆孔遇洞率在3%以內，巖溶及溶洞發育相對較少，為巖溶弱發育段。

(3)下發育層。埋深約65 m以下(標高段一般為-45以下)，該段溶洞發育很不均勻，斷層帶溶洞發育規模較大，最高溶洞達7.08 m，這與地層的構造有密切的關系，此標高段鉆孔遇洞率為15.6%，分布在阮家巷及毛坦港附近，全部為重點塌陷區，區內深部巖溶應該較發育。

2 數學模型

2.1 地下水流動數學模型

根據Darcy定律和水均衡原理，忽略密度變化對地下水影響的三維流動數學模型為：

(1)

式中：h(x,y,z,t)為水位，L；Kxx、Kyy和Kzz分別為沿x、y、z主方向的滲透系數，L/T;w為單位體積單位時間內注入或抽出的水量，1/T；μs為單位儲水系數，1/L；t為時間，T；h0(x,y,z)為流場初始水位，L；B1為一類邊界條件；B2為二類邊界條件；q(x,y,z,t)為二類邊界上單位面積已知流量，L/T。

2.2 力學模型的簡化

假定地層在土壓力變化過程中引起的垂向總應力不變，土體側向變形忽略，只發生垂向變形。土體的有效應力變化及垂向應變為：

有效應力變化為：

(2)

Δσ′z=-αpΔP

(3)

平均有效應力變化為:

(4)

垂向應變為：

(5)

式中：ΔP為壓力增量；Δσ′x、Δσ′y和Δσ′z分別為x，y和z三個方向的有效應力變化；Δσ′M為平均有效應力，σ′M=(σ′x+σ′y+σ′z)/3；v為泊松比；K為體積模量；αp為Biot系數；εz為垂向應變。由此模型可以看出，孔隙水壓力的增加引起水平有效應力減小；由于地表的自由變形，垂向有效應力只受孔隙水壓力變化的影響；垂向應力應變滿足線彈性的Hooke定律。

2.3 耦合方法

在初始流場時有相對應初始應力值， 地下水位上升或下降都會引起土體有效應力的變化，利用GMS軟件模擬出不同時段地下水的流場值，計算出不同時段的應力，不同時段之間應力差值是壓力的增量為ΔP，不同時段之間有效應力差值為Δσ′x、Δσ′y和Δσ′z，再計算出不同時段土體的應力應變。飽和土體的有效應力表達式為：

(6)

2.4 土體破壞判別準則

土體是否產生應力破壞，采用庫倫-莫爾準則：

τ=c+σ′ntanφ

(7)

首先需要根據應力狀態得到應力莫爾圓，然后確定主應力(有效應力)方向，再得到最有可能發生剪切破壞的方向，最后采用破壞準則判定巖土體是否破壞。

(8)

式中：σ1、σ3最大、最小主應力，上式表明巖土體處于極限平衡狀態，還有另外兩種情況存在。可以得到如下判別關系式：

(9)

每一個單元體相應的最大應力及最小主應力為：

(10)

3 數值模型

3.1 模型邊界條件及網格剖分

(1)側向邊界。研究區西臨長江，可設此邊界為河流邊界，根據資料分析，長江歷年的平均水位16.70 m，因長江切割了砂層及基巖，且與承壓含水層及巖溶水有水力聯系，在汛期側滲補給前緣孔隙水，通過越流形式補給巖溶水，所以砂層以下定為隔水邊界，砂層以上至長江底部高程以下定為定流量邊界，分段量化補給量。模擬區的東、南、北面邊界定為定水頭邊界，根據觀測孔的資料分析得到。

(2)底邊界。選在灰巖埋深165 m處，此邊界基本沒有產生水力交換，可以作為隔水層邊界處理。

(3)上邊界。潛水層的自由水面作為系統的上邊界，作為潛水面邊界，系統與外界發生垂向水流交換。如：接受河流滲漏補給、農田灌溉補給、降雨入滲補給、潛水面的蒸發排泄等。各含水層的頂、底板高程根據鉆孔資料確定。

通過對該研究區的地質結構模型、地下水水流特征及補、徑、排關系分析，將蓋層地下水流模型概化為非均質各向同性、三維非穩定流模型。基巖層定義為非均質各向異性、三維非穩定流模型。

把研究區共剖分網格單元6 315個(圖2)。垂向上由各含水層底板埋深及厚度決定，總模擬深度180 m。

圖2 網格系統Fig.2 Grid system

3.2 模型參數分區及參數值

(1)滲透系數、彈性儲水系數、給水度。滲透系數、彈性儲水系數、給水度是正確建立水流模型的關鍵，在空間上具有很強的變異性，本次建立水流模型采用等效厚度法，結合鉆孔、地質剖面圖、觀測孔資料、研究區前人的研究成果確定初始參數值，在模型識別階段進一步調整，圖3為各含水層參數分區圖。灰巖及其他基巖的參數根據灰巖的埋藏及覆蓋情況進行分區。

圖3 潛水、弱透水層、承壓含水層，其他基巖及灰巖參數分區Fig.3 Phreatic aquifer, aquitard, confined aquifer, other bedrock and limestone parameters partition

(2)力學參數初次選取。力學參數有壓縮體積模量、回彈體積模量、天然重度、飽和重度、泊松比、黏聚力、摩擦角等。

3.3 模型識別

模型識別的目的是為了提高模型的精度，模型識別的過程是反復調整模型參數(包括邊界條件)使計算結果盡可能與實際監測資料一致。模型識別資料選用2009年1月全區的統測水位作為模型識別的初始流場，2012年1月和7月的數據作為擬合檢驗。

(1)流場擬合。根據2009年1月為初始流場，計算出位于土洞上層即第三層的初始應變值及破壞分布如圖4、5所示。

圖4 2009年1月第三層垂向應變Fig.4 Vertical strain in lay 3 in January 2009

圖5 2009年1月第三層破壞區Fig.5 Destroyed zone in lay 3 in January 2009

1月是武漢地區的枯水期、7月是豐水期，長江水位的上漲、降雨集中在7月，2012年1月流場實測值與計算值比較(圖6)、7月土體應變及破壞區的分布情況如圖7、8所示。

圖6 2012年1月第三層地下水流場實測值與計算值比較Fig.6 Comparison of measured and calculated valuesfor groundwater flow field contours in lay 3 in January 2012

圖7 2012年7月第三層垂向應變Fig.7 Vertical strain in lay 3 in July 2012

圖8 2012年7月第三層破壞區Fig.8 Destroyed zone in lay 3 in July 2012

(2)應力及變形擬合。2012年監測點土壓力及土體變形擬合效果如圖9、10所示。

圖9 監測點土壓力擬合Fig.9 Earth pressure fitting of monitoring points

圖10 監測點土體變形擬合Fig.10 Soil deformation fitting of monitoring points

(3)模擬結果分析。土體壓力自2009年開始模擬計算，模擬數據表明， 2009年1月第三層砂層處于破壞狀態，這與2009年發生在烽火村的巖溶塌陷時間一致，2012年7月模擬的土體破壞區與地下水流場變化有著密切的關系，與實際的塌陷點位置相符。充分說明利用流場值模擬土體的變形及巖溶塌陷是可行的。

土體受力狀況的變化值受地下水的影響較大，每年的7-8月份相對于全年的土壓力值及有效應力值偏低，是由于地下水位明顯升高，產生浮托力作用的結果。在模擬的年度內土壓力值的震蕩性較大，隨著地下水位的變化而壓實回彈，經過幾次回合后，土體出現塑性變形，之后發展到破壞引起巖溶塌陷。

土體垂向變形與地下水的升幅有直接的聯系，土體的沉降與回彈與每一層地下的流場走向相符，在沒有考慮施加荷載時，土體沉降的變形值一般小于4 mm。

土體破壞規律：地下水的活動和運移，將對土層產生潛蝕作用，從而形成土洞。此外，土洞形成后，其洞壁周圍將產生應力集中現象，當地下水位發生變化時(上升或下降時)，將進一步改變洞壁周圍土體的應力狀態，并致使洞體周邊產生破壞最終導致塌陷。

(4)識別后的參數，通過模型識別，得到的模型參數如表1-5所示。

表1 潛水層參數分區表Tab.1 Phreatic aquifer parameters partition table

表2 弱透水層參數分區表Tab.2 Aquitard parameters partition table

表3 承壓含水層參數分區表Tab.3 Confined aquifer parameters Partition table

表4 其他基巖及灰巖參數分區表Tab. 4 Other bedrock and limestone parameters Partition table

表5 壓縮回彈模型相關參數表Tab.5 The parameters of the compression and springback model

4 巖溶塌陷影響因素模擬評價

4.1 長江水位動態對巖溶塌陷的影響預測

實際調查監測表明，研究區內地下水水位動態受長江水位波動影響明顯，地下水位動態與長江水位具有顯著的相關性。因此，以長江水位變幅較大的2010年水位動態數據，在其他各項參數及源匯項不變的情況下，模擬預測了長江水位動態對巖溶塌陷的影響。

經過兩次交替長江水位上升、下降，土體也經歷了回彈、壓縮作用，預測的第三層(砂層)應變及土體破壞情況如圖11、12所示。結果表明，研究區地下水與長江水位具有很好的連通性，當長江水位在最高水位(年)與最底水位(低)交替出現時，地面塌陷點向長江方向偏移，這與地下水的排泄方向一致，巖溶塌陷很明顯是由地下水滲流作用引起。

圖12 長江水位交替時第三層破壞區Fig.12 Destroyed zone in lay 3 when water level changes in the Yangtze River water

4.2 工程建設排水施工對巖溶塌陷的影響預測

研究區近幾年拆遷還建、市政建設等工程大量增加，很多工程涉及到深基坑開挖，為保證施工安全，需要降低地下水。根據施工規范規定，地下水位必須降至基礎以下1～1.5 m，一般情況下基礎的開挖深度10 m左右，研究區地下水位需降至高程11 m位置。為此，利用數值模型模擬了當某基坑開挖排水時得到的垂向應變及土體破壞情況如圖13、14所示。模擬結果顯示，土體的垂向應變發生明顯變化，從降落漏斗最深處開始破壞。

圖13 基坑排水第三層垂向應變Fig.13 Vertical strain in lay 3 when drainage of foundation pit

圖14 基坑排水第三層破壞區 Fig.14 Destroyed zone in lay 3 when drainage of foundation pit

4.3 大氣降水對巖溶塌陷的影響預測

降雨是研究區地面塌陷發生的重要影響因素之一，實際調查表明，研究區內多起地面塌陷的觸發因素都是降雨。降雨對地面塌陷的直接作用主要體現在兩個方面：使上部土體飽水自重增加，物理力學性質降低；降雨入滲的滲流作用破壞土洞的穩定性。

根據武漢地區的降雨特征，在7月份流場的基礎上，通過改變降雨入滲補給量，模擬3種降雨模式(表6)，計算土壓力的變化值分析降雨對土體變形及巖溶塌陷的影響。

在補給條件不變的情況下，以上幾種降雨模式的土體變形及破壞情況如圖15-20。

表6 降雨預測模式有以下3種情況Tab.6 Rainfall prediction model has the following three cases

圖15 模式1第三層垂向應變Fig.15 Strain contour in Layer 3 for case 1

圖16 模式1第三土層破壞區Fig.16 Destroy zone in Layer 3 for case 1

圖17 模式2第三層垂向應變Fig.17 Strain contour in Layer 3 for case 2

圖18 模式2第三土層破壞區Fig.18 Destroy zone in Layer 3 for case 2

圖19 模式3第三層垂向應變Fig.19 Strain contour in Layer 3 for case 3

圖20 模式3第三層破壞區Fig.20 Destroy zone in Layer 3 for case 3

模擬結果表明，研究區降雨對于土洞穩定性有明顯影響。三種模式下應力分布規律大致相同，由于區內地下水與長江有著較直接的排泄途徑，短期的降雨沒有造成水位大幅度上升，但使水壓力增大導致土洞破壞。由以上模擬分析的變形位移圖可知，相同時間下，降雨強度越大，土洞變形越大，同時降雨強度較小但持續時間較長時，對土洞的影響幾乎和降雨強度較大的暴雨相同，也可成為塌陷的主要誘發因素。

5 結 論

(1)基于地下水動力學理論和Terzaghi固結沉降理論，采用MODFLOW軟件和自編Terzaghi模型計算程序，建立了青菱鄉烽火村巖溶塌陷區地下水動力學-巖土力學耦合數值模型。

(2)利用該模型模擬了長江水位動態、基坑排水、降雨等因素對巖溶塌陷的影響，結果表明：長江水位上升、下降交替出現時，土體發生破壞且塌陷點偏移長江方向；周邊工程建設排水施工時，降落漏斗最深處沉降量4 mm，對土洞的破壞性影響最大；降雨強度不大，但降雨時間長造成的土體變形大，最大處引起土體上浮1～2 mm，大大降低了蓋層的強度。

□

[1] Lu Y R，Liu Q，Zhang F E，et al. Environmental characteristics of karst in China and their effect on engineering[J]. Carbonates and Evaporites，2013，28(1-2)：251-258.

[2] 郭殿權，劉道喜. 武昌陸家街地面塌陷成因機理分析[J].工程勘察，1990,(6):11-13.

[3] He K Q，Zhang S Q，Wang F，et al.The karst collapses induced by environmental changes of the groundwater and their distribution rules in North China [J]. Environmental Earth Sciences，2010，61(5):1 075-1 084.

[4] 胡亞波，劉廣潤，肖尚德,等.一種復合型巖溶地面塌陷的形成機理——以武漢市烽火村塌陷為例[J].地質科技情報，2007，26(1):96-100.

[5] 張人權，梁 杏，靳孟貴，等.水文地質學基礎[M].6版.北京，地質出版社，2011.

[6] 方 云，林 彤，譚松林. 土力學[M]. 武漢: 中國地質大學出版社，2003.

[7] He K Q，Jia Y Y，Zhao M，et al. Comprehensive analysis and quantitative evaluation of the influencing factors of karst collapse in groundwater exploitation area of Shiliquan of Zaozhuang，China[J]. Environmental Earth Sciences，2012，66(8):2 531-2 541.

[8] Zhao H J，Ma F S，Guo J，et al. Regularity and formation mechanism of large-scale abrupt karst collapse in southern China in the first half of 2010[J]. Natural Hazards, 2012，6(3)：1 037-1 054.

[9] Liu X M，Chen C X. Prediction and evaluation of space collapse for typical covered karst[J]. Disaster Advances，2010，3(4): 120-126.

[10] He K Q，Jia Y Y，Wang B. et al.Comprehensive fuzzy evaluation model and evaluation of the karst collapse Susceptibility in Zaozhuang Region，China [J]. Natural Hazards，2013，68(2):613-629.

[11] Haijun Zhao o Fengshan Ma o Jie Guo. Regularity and formation mechanism of large-scale abrupt karst collapse in southern China in the first half of 2010[J].Natural Hazards，2012，60:1 037-1 054.