基于FEFLOW的廣州金沙洲地區(qū)地下水流場數(shù)值模擬研究

陳小月，黃健民，盧 薇

(1.廣東省廣州市地質調查院，廣東 廣州 510440，2.廣州地理研究所，廣東 廣州 510070)

自2007年以來，廣州金沙洲陸續(xù)發(fā)生巖溶地面塌陷、地面沉降地質災害[1-3]，據(jù)已有的相關研究，上述地質災害的發(fā)生與地下水的異常波動密切相關，并主要發(fā)生在地下水位劇烈波動的下降及恢復期間。為了更好、更深入分析金沙洲地下水變化情況與地面塌陷、地面沉降地質災害的關系，本文根據(jù)金沙洲地帶的水文地質條件、鉆孔、地下水的水位動態(tài)及抽水試驗成果，應用 FEFLOW軟件建立金沙洲地下水系統(tǒng)的數(shù)值計算模擬模型。模擬區(qū)范圍以金沙洲界線為界，將整個金沙洲作為地下水系統(tǒng)的模擬計算范圍。

FEFLOW是由德國柏林水資源規(guī)劃與系統(tǒng)研究所(WASY)開發(fā)的一種對地下水水量及水質進行模擬的軟件系統(tǒng)，利用它不僅可以計算出水位、溶液濃度和溫度等標量數(shù)據(jù)，而且還可以模擬降水、地表水、地下水的流動與轉換，計算出流速、流線及流徑線等向量數(shù)據(jù)，很適合分析模擬區(qū)地質災害與地下水流場的關系[4-6]。

1 模擬區(qū)水文地質概況

金沙洲四面環(huán)水，東臨沙貝?！咨澈?，西側約1 km為水口涌，沙貝?！咨澈蛹八谟烤鶠榱飨又Я?，兩支流在金沙洲南部匯入珠江。區(qū)域水文主要受珠江干流控制，局部受流溪河影響。在金沙洲模擬區(qū)內(nèi)，其西側低丘陵的地層為石炭系下統(tǒng)大賽壩組的砂巖、頁巖、泥巖等，這些巖層隔水性好;東側石炭系壺天群、石磴子組灰?guī)r及白堊系大塱山組泥質粉砂巖與珠江白沙河及沙貝海相連，地下水與珠江水存在水力聯(lián)系;南面石磴子組灰?guī)r及大塱山組泥質粉砂巖與佛山市黃岐相接，與區(qū)外地下水存在水量交換。因此，金沙洲具有完整的補、徑、排條件，形成了一個相對獨立的一級水文地質單元。在自然狀態(tài)下，金沙洲地下水順勢總體自北西向南東流動，最終排泄于珠江水系。

2 地下水系統(tǒng)概念模型

為建立模擬區(qū)地下水系統(tǒng)的概念模型，首先要對實際的水文地質條件加以概化，簡化或忽略與系統(tǒng)無關的某些系統(tǒng)要素和狀態(tài)，以便于數(shù)學描述，并建立地下水系統(tǒng)模擬模型。

2.1 模擬區(qū)含水層結構的概化

根據(jù)金沙洲區(qū)內(nèi)的地下水系統(tǒng)含水介質的物質組成及水文地質特性，可將其內(nèi)部結構概化為四層，由上至下分別為:

(1)弱含水層(弱透水層):由人工填土、淤泥質粘土、淤泥組成。

(2)第一含水層(第四系含水層):主要由砂礫層夾砂，中粗砂、中砂、粗砂組成。

(3)隔(弱透)水層:由殘積礫質粘性土組成。

(4)第二含水層(基巖巖溶裂隙含水層):主要由壺天群和石磴子組灰?guī)r組成。

2.2 模擬區(qū)邊界條件概化

金沙洲模擬區(qū)東面緊鄰珠江水系，并與模擬區(qū)地下水有水量較緩，可視為定水頭邊界;西面丘陵地區(qū)為金沙洲地下水含水系統(tǒng)的分水嶺，其山脊線可視為零流量邊界，即為隔水邊界;在模擬區(qū)南側，第一層含水層和第二含水層與外圍地下水有水量交換，主要為地下徑流側向補給，可作為徑流邊界;人類工程活動抽排地下水導致的降水中心概化為模型的第四類邊界。潛水含水層自由水面為系統(tǒng)的上邊界，通過該邊界的潛水與系統(tǒng)外發(fā)生垂向水量交換;下部以基巖中風化層底板為底部邊界，按零通量邊界處理。

3 地下水系統(tǒng)數(shù)值模擬模型

3.1 地下水系統(tǒng)數(shù)學模型

從金沙洲模擬區(qū)的空間上看，地下水流整體上以水平運動為主、垂向運動為輔，方向性明顯，且地下水系統(tǒng)的輸入、輸出隨時間和空間變化，為了準確模擬各含水層的相互影響，將模擬區(qū)內(nèi)的地下水流作為非均質、各向異性介質中的三維非穩(wěn)定流處理，建立地下水系統(tǒng)數(shù)學模型如下:

式中:Ω為滲流區(qū)域;h為含水層的水位標高(m);Kx、Ky、Kz為分別為 x、y、z方向的滲透系數(shù)(m/d);Kn為邊界面法向方向的滲透系數(shù)(m/d);S為自由面以下含水層儲水系數(shù)(1/m);μ為潛水含水層在潛水面上的重力給水度;ε為含水層的源匯項(1/d);p為人工開采和降水等(1/d);h0為含水層的初始水位分布(m);φ1為第一類邊界上的已知函數(shù);Γ1為滲流區(qū)域的上邊界，即地下水的自由表面;Γ2為滲流區(qū)域的側向邊界;Γ4為滲流區(qū)域的下邊界，即承壓含水層底部的隔水邊界;n→為邊界面的法線方向;q(x，y，z，t)為定義為二類邊界的單寬流量(m2/d·m)，流入為正，流出為負，隔水邊界為0。

在所建立地下水系統(tǒng)概念模型中，所有邊界均按照第二類邊界條件來處理。

3.2 模擬區(qū)網(wǎng)格劃分

將模擬區(qū)范圍內(nèi)平面(即每一個模擬層)剖分為9 350個結點，14 092個三角單元。模型垂向上共分4個模擬層，4層模擬層的總結點數(shù)為9 350×4=37 400個，總單元數(shù)為14 092×4=56 368個。金沙洲模擬區(qū)的三維網(wǎng)格剖分圖見圖1。

圖1 金沙洲地下水數(shù)值模擬區(qū)網(wǎng)格剖分圖

3.3 模擬區(qū)水文地質參數(shù)劃分

金沙洲地下水數(shù)值模擬模型主要的水文地質參數(shù)為滲透系數(shù)、給水度、儲水系數(shù)及彈性釋水系數(shù)等。

為保證對計算精度不產(chǎn)生影響，并使水文地質參數(shù)分區(qū)簡化，將每一個滲透系數(shù)分區(qū)同時視為給水度、儲水系數(shù)分區(qū)。滲透系數(shù)的初始分區(qū)按地形地貌、沉積類型和地層巖性等進行劃分，如表1、表2及圖2所示。各區(qū)的參數(shù)初始值根據(jù)抽水試驗、沉積類型和地層巖性、地層時代等特征進行估值，垂向滲透系數(shù)根據(jù)類似地區(qū)的計算經(jīng)驗估計給值。

表1 含水層給水度和釋水系數(shù)分區(qū)表

表2 第四系砂層滲透系數(shù)分區(qū)表 m/d

圖2 金沙洲地下水數(shù)值模擬計算參數(shù)分區(qū)圖

基巖裂隙含水層考慮到斷層對含水層導水性的影響，沿斷層的方向生成緩沖區(qū)，加大緩沖區(qū)的滲透系數(shù)，見表3。

3.4 源匯項處理

源匯項主要考慮地下水開采量、垂向入滲補給量和蒸發(fā)量?，F(xiàn)狀條件下，模擬區(qū)地下水開采主要為工程活動抽排地下水，概化為井單元處理。垂向入滲主要為大氣降水的入滲補給和地表水的入滲補給，分別取多年加權平均降雨量和河流補給系數(shù)法確定補給量，模擬過程中不考慮蒸發(fā)量。

表3 基巖巖溶裂隙含水層滲透系數(shù)分區(qū)表 m/d

3.5 模型檢驗與校正

為了驗證所建立的數(shù)學模型和模型識別后確定的水文地質參數(shù)的可靠性，采用試估—校正法，取 K36、K37和 K723三個鉆孔觀測水位數(shù)據(jù)與模擬水位數(shù)據(jù)進行對比，對數(shù)值模擬模型進行檢驗。根據(jù)收集到的地下水觀測資料，取模型識別的時段為2009年7月13日至2009年12月31日，取時間步長為1天，共計1 171個時間步長，以實際水位和計算水位的擬合情況作為調整參數(shù)的依據(jù)。從圖3中可以看出，經(jīng)過識別的參數(shù)計算的水位線與實測水位線趨勢相同，誤差都在容許范圍之內(nèi)。

圖3 實測水位與模擬水位檢驗圖

4 地下水流場數(shù)值模擬計算

模擬區(qū)人類工程活動抽排地下水造成地下水位波動大致分為如下五個階段(見表 4)，可概化為 1 200 m3/d、2 300 m3/d和3 200 m3/d等三種典型情況，現(xiàn)以抽排水量3 200 m3/d(第四階段)為例進行金沙洲地下水流場數(shù)值模擬計算，模擬結果見圖4和圖5。

表4 模擬區(qū)抽排地下水階段劃分及其抽排水量

圖4 抽排水量為3 200 m3/d時第四系砂層含水層模擬區(qū)流場分布圖

圖5 抽排水量為3 200 m3/d時基巖巖溶裂隙含水層模擬區(qū)流場分布圖

從地下水流場分布圖中可以看出，模擬區(qū)地下水向施工降水中心(豎井)流動，形成以人類工程施工降水點為中心的降落漏斗地段。地面塌陷的分布范圍位于橫沙村—沙貝村—鳳岡隔水條帶以西，說明隔水條帶阻擋了工程降水并接受來自于東側珠江河水的補給，使金沙中學附近的地下水沿斷層強徑流帶由東北側向工程施工降水點方向流動，導致這些地段的地下水流場動態(tài)變化加劇。

野外調查發(fā)現(xiàn)的T9-T14六宗地面塌陷地質災害均分布于工程施工降水點附近及其所控制地面沉降范圍內(nèi)(BX1)(圖4、圖5);地面塌陷根據(jù)地下水監(jiān)測資料，在2008年8月15日至2009年5月8日時間段內(nèi)，金沙洲地下水出現(xiàn)異常波動，總體處于快速下降過程，部分民井逐漸干涸，在此期間，金沙洲共發(fā)生地面塌陷地質災害14宗，占塌陷總數(shù)的58.3%(圖6)。因此，不論從空間還是時間角度來看，模擬區(qū)地面塌陷地質災害的發(fā)生都和模型分析結果吻合，說明了該模型的建立是合理的。

圖6 巖溶地面塌陷與地下水位動態(tài)變化關系圖

5 數(shù)值模擬的結論與認識

金沙洲模擬區(qū)內(nèi)斷裂構造發(fā)育，斷層及不同規(guī)模的節(jié)理裂隙密集，造成可溶巖發(fā)育多層裂隙和溶洞，洞隙之間連通性好，構成區(qū)域地下水的主要徑流場。人類工程施工降水是引起金沙洲地下水流場變化的主要原因:抽排水期間，地下水的流動方向由四周向降水中心流動，地下水流場分布圖內(nèi)斷層發(fā)育地段的地下水流線加長，并形成以施工抽水點為中心的降落漏斗，地下水位變化明顯，易引發(fā)地面塌陷和地面沉降地質災害。因此，基于上述模擬結果，應嚴格控制金沙洲區(qū)域的工程活動施工降水強度，以有效防止地面塌陷及地面沉降地質災害的發(fā)生。

[1]黃健民，鄭小戰(zhàn)，陳小月，等.廣州金沙洲廣州大學附屬實驗學校地面沉降形成機理研究[J].熱帶地理.2012，32(04):338-343.

[2]郭宇，黃健民，陳建新，等.廣州市白云區(qū)金沙洲地區(qū)地質災害風險區(qū)劃[J].熱帶地理.2013，33(6):659-665.

[3]郭宇，黃健民，周志遠，等.廣州市白云區(qū)金沙洲地區(qū)地質災害現(xiàn)狀及防治對策[J].中國地質災害防治學報.2013，24(3):100-104.

[4]趙國紅，寧立波，王現(xiàn)國.新鄭市淺層地下水流數(shù)值模擬及評價[J].地下水.2007，29(6):43-46.

[5]郭曉東，田輝，張梅桂，等.我國地下水數(shù)值模擬軟件應用進展[J].地下水.2010，32(4):5-7.

[6]張立志.基于 FEFLOW的大興區(qū)地下水動態(tài)模擬研究[D].地質大學(北京).1-66，2009.