巖溶塌陷水動力-力學耦合過程數(shù)值模擬
——以武漢市青菱鄉(xiāng)為例

陳冬琴，唐仲華

(1.中國地質(zhì)大學(武漢)環(huán)境學院，武漢 430074；2.湖北城市建設(shè)職業(yè)技術(shù)學院，武漢 430000)

1 研究區(qū)巖溶塌陷區(qū)地質(zhì)條件

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

垂向上分成三個發(fā)育段。

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

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

圖1 巖溶發(fā)育程度與深度關(guān)系Fig.1 Degree of karst development vs. the depth

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

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

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

2 數(shù)學模型

2.1 地下水流動數(shù)學模型

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

(1)

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

2.2 力學模型的簡化

假定地層在土壓力變化過程中引起的垂向總應(yīng)力不變，土體側(cè)向變形忽略，只發(fā)生垂向變形。土體的有效應(yīng)力變化及垂向應(yīng)變?yōu)椋?/p>

有效應(yīng)力變化為：

(2)

Δσ′z=-αpΔP

(3)

平均有效應(yīng)力變化為:

(4)

垂向應(yīng)變?yōu)椋?/p>

(5)

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

2.3 耦合方法

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

(6)

2.4 土體破壞判別準則

土體是否產(chǎn)生應(yīng)力破壞，采用庫倫-莫爾準則：

τ=c+σ′ntanφ

(7)

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

(8)

式中：σ1、σ3最大、最小主應(yīng)力，上式表明巖土體處于極限平衡狀態(tài)，還有另外兩種情況存在?？梢缘玫饺缦屡袆e關(guān)系式：

(9)

每一個單元體相應(yīng)的最大應(yīng)力及最小主應(yīng)力為：

(10)

3 數(shù)值模型

3.1 模型邊界條件及網(wǎng)格剖分

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

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

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

通過對該研究區(qū)的地質(zhì)結(jié)構(gòu)模型、地下水水流特征及補、徑、排關(guān)系分析，將蓋層地下水流模型概化為非均質(zhì)各向同性、三維非穩(wěn)定流模型?；鶐r層定義為非均質(zhì)各向異性、三維非穩(wěn)定流模型。

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

圖2 網(wǎng)格系統(tǒng)Fig.2 Grid system

3.2 模型參數(shù)分區(qū)及參數(shù)值

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

圖3 潛水、弱透水層、承壓含水層，其他基巖及灰?guī)r參數(shù)分區(qū)Fig.3 Phreatic aquifer, aquitard, confined aquifer, other bedrock and limestone parameters partition

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

3.3 模型識別

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

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

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

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

1月是武漢地區(qū)的枯水期、7月是豐水期，長江水位的上漲、降雨集中在7月，2012年1月流場實測值與計算值比較(圖6)、7月土體應(yīng)變及破壞區(qū)的分布情況如圖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月第三層垂向應(yīng)變Fig.7 Vertical strain in lay 3 in July 2012

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

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

圖9 監(jiān)測點土壓力擬合Fig.9 Earth pressure fitting of monitoring points

圖10 監(jiān)測點土體變形擬合Fig.10 Soil deformation fitting of monitoring points

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

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

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

土體破壞規(guī)律：地下水的活動和運移，將對土層產(chǎn)生潛蝕作用，從而形成土洞。此外，土洞形成后，其洞壁周圍將產(chǎn)生應(yīng)力集中現(xiàn)象，當?shù)叵滤话l(fā)生變化時(上升或下降時)，將進一步改變洞壁周圍土體的應(yīng)力狀態(tài)，并致使洞體周邊產(chǎn)生破壞最終導致塌陷。

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

表1 潛水層參數(shù)分區(qū)表Tab.1 Phreatic aquifer parameters partition table

表2 弱透水層參數(shù)分區(qū)表Tab.2 Aquitard parameters partition table

表3 承壓含水層參數(shù)分區(qū)表Tab.3 Confined aquifer parameters Partition table

表4 其他基巖及灰?guī)r參數(shù)分區(qū)表Tab. 4 Other bedrock and limestone parameters Partition table

表5 壓縮回彈模型相關(guān)參數(shù)表Tab.5 The parameters of the compression and springback model

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

4.1 長江水位動態(tài)對巖溶塌陷的影響預測

實際調(diào)查監(jiān)測表明，研究區(qū)內(nèi)地下水水位動態(tài)受長江水位波動影響明顯，地下水位動態(tài)與長江水位具有顯著的相關(guān)性。因此，以長江水位變幅較大的2010年水位動態(tài)數(shù)據(jù)，在其他各項參數(shù)及源匯項不變的情況下，模擬預測了長江水位動態(tài)對巖溶塌陷的影響。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

(2)利用該模型模擬了長江水位動態(tài)、基坑排水、降雨等因素對巖溶塌陷的影響，結(jié)果表明：長江水位上升、下降交替出現(xiàn)時，土體發(fā)生破壞且塌陷點偏移長江方向；周邊工程建設(shè)排水施工時，降落漏斗最深處沉降量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] 郭殿權(quán)，劉道喜. 武昌陸家街地面塌陷成因機理分析[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].地質(zhì)科技情報，2007，26(1):96-100.

[5] 張人權(quán)，梁 杏，靳孟貴，等.水文地質(zhì)學基礎(chǔ)[M].6版.北京，地質(zhì)出版社，2011.

[6] 方 云，林 彤，譚松林. 土力學[M]. 武漢: 中國地質(zhì)大學出版社，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.