鄭玉萍，王 巍，韓 曄，王 淼，陸 陽(yáng)

(1.天津市控制地面沉降工作辦公室，天津300061;2.天津華北地質(zhì)勘查局地質(zhì)研究所，天津300170)

天津市屬于典型的因資源性缺水而引發(fā)地面沉降的城市。地面沉降的主要原因是由于水資源的短缺，被迫過量開采地下水，從而引起包括地面沉降在內(nèi)的一系列不良環(huán)境地質(zhì)效應(yīng)。第四系松散巖類含水層的孔隙水壓力(水位)是與第四系地層壓縮密切相關(guān)的一個(gè)重要指標(biāo)。本文選取西青區(qū)辛口鎮(zhèn)及其周邊地區(qū)，選用Processing Modflow軟件，利用2005年7水位統(tǒng)測(cè)數(shù)據(jù)和分層標(biāo)數(shù)據(jù)，開展了第四系孔隙水流場(chǎng)和地面沉降量的數(shù)值模擬，并對(duì)模擬結(jié)果進(jìn)行分析，以求對(duì)今后相關(guān)工作起到參考性意義。

1 研究區(qū)概況

1.1 模擬區(qū)簡(jiǎn)介

本次為研究地下水開采對(duì)地面沉降的影響，選取以西青鄭莊子分層標(biāo)為中心的數(shù)值模擬范圍，模擬區(qū)范圍10 km×10 km×566 m(長(zhǎng)×寬×高);其模擬范圍內(nèi)沉降有明顯的分帶性，模擬范圍內(nèi)東北部沉降在30～50 mm/a之間，西南部沉降大于50 mm/a;與本區(qū)水位等值線比較后，初步判斷該模擬范圍地面沉降主要受本區(qū)地下水開采影響。

1.2 地質(zhì)及水文地質(zhì)背景

西青區(qū)位于華北平原東北端，構(gòu)造單元屬新華夏構(gòu)造體系第二沉降帶滄縣隆起中部，發(fā)育有一系列次級(jí)凸起和凹陷。根據(jù)已有鉆孔資料顯示，本區(qū)566.27 m以淺范圍內(nèi)自下而上發(fā)育有上新統(tǒng)明化鎮(zhèn)組、下更新統(tǒng)楊柳青組、中更新統(tǒng)佟樓組、上更新統(tǒng)塘沽組以及全新統(tǒng)天津組地層，為一套350萬(wàn) a以來(lái)的沉積地層(見表1)。

表1 研究區(qū)沉積地層年代及巖性表

1.3 地面沉降背景

自上世紀(jì)六、七十年代以來(lái)，由于社會(huì)經(jīng)濟(jì)發(fā)展的需要，在水資源匱乏的情況下，全市長(zhǎng)期超采地下水，至1985年沉降最大速率達(dá)到了80～100 mm/a。1986年后進(jìn)入沉降治理階段，大部分地區(qū)沉降明顯減緩，市區(qū)沉降速率降低到10～15 mm/a。距市區(qū)較近的西青區(qū)是近年來(lái)沉降比較嚴(yán)重的地區(qū)之一，西青區(qū)2009年度平均沉降量43 mm，2010年度平均沉降量53 mm，其年沉降速率遠(yuǎn)大于天津市區(qū)。研究區(qū)所處地段1985年至今的累計(jì)沉降量達(dá)到了1 400 mm以上。2004年7月西青區(qū)鄭莊子分層標(biāo)建成，2006、2007、2009年沉降速率分別為 64.43、57.06、47.27 mm/a。同時(shí)，本次模擬研究區(qū)沉降速率較西青區(qū)全區(qū)平均沉降量要大，具有充分的代表性。

2 水文地質(zhì)概念模型及數(shù)學(xué)模型

2.1 水文地質(zhì)概念模型

本區(qū)處于平原區(qū)。通過收集大量的地質(zhì)資料及以往野外勘探情況，深入分析模擬區(qū)的含水層結(jié)構(gòu)、地下水流場(chǎng)和流動(dòng)特征、地下水均衡項(xiàng)，并建立水文地質(zhì)概念模型，為建立數(shù)值模型提供依據(jù)。

2.1.1 含水層劃分

根據(jù)地層形成的時(shí)代、成因、巖性和含水層的埋藏條件、水質(zhì)及地下水動(dòng)態(tài)特征等各方面因素，將埋深566m內(nèi)地層中的含水巖組概化為5層(見表2)。

表2 研究區(qū)概化含水層一覽表

2.1.2 邊界條件

側(cè)向邊界:模擬區(qū)地勢(shì)平坦但不是完整的水文地質(zhì)單元，邊界的流入流出量難以計(jì)算。因此模擬區(qū)側(cè)向邊界均設(shè)置為通用水頭邊界。通用水頭邊界是模型特有的邊界處理方式，可通過給入每個(gè)應(yīng)力期邊界的水流參數(shù)、模型外的水位，模型自行計(jì)算流入或流出水量，從而調(diào)整邊界上的水位，來(lái)計(jì)算邊界的側(cè)向交換量。

垂向邊界:整個(gè)研究區(qū)在垂向上分為5個(gè)含水層組，根據(jù)含水層富水特征和開采條件來(lái)看，第二至第四含水層組研究意義較大，這幾個(gè)含水層組在上部主要接受上一含水層組的越流補(bǔ)給，第五含水層組下部邊界因下層沒有分層標(biāo)監(jiān)測(cè)數(shù)據(jù)，不考慮沉降，模型中處理成隔水邊界。第一含水巖組接受大氣降水入滲補(bǔ)給，蒸發(fā)和越流是其主要排泄方式，開采量小。從整個(gè)系統(tǒng)來(lái)看，垂向上只接受大氣降水補(bǔ)給，蒸發(fā)、開采是其主要的排泄方式，系統(tǒng)內(nèi)部各含水層組之間在垂向上通過越流產(chǎn)生部分交換量。

2.1.3 地下水流特征

從空間上看，模擬區(qū)地下水流整體上以水平運(yùn)動(dòng)為主、垂向運(yùn)動(dòng)為輔，地下水系統(tǒng)符合質(zhì)量守恒定律和能量守恒定律;在常溫常壓下地下水運(yùn)動(dòng)符合達(dá)西定律;考慮含水層之間的流量交換，地下水運(yùn)動(dòng)可以概化為空間三維流;地下水系統(tǒng)的垂向運(yùn)動(dòng)是由層間水頭差異引起的;地下水系統(tǒng)的輸入輸出隨時(shí)間、空間變化，故地下水為非穩(wěn)定流;參數(shù)隨空間變化，體現(xiàn)了系統(tǒng)的非均質(zhì)性，所以含水介質(zhì)概化為非均質(zhì)各向同性介質(zhì)。

綜上所述，模擬區(qū)可概化成非均質(zhì)、水平方向各向同性、垂向上存在變異、空間三維結(jié)構(gòu)、非穩(wěn)定地下水流系統(tǒng)，即地下水系統(tǒng)的概念模型。

2.2 數(shù)學(xué)模型

2.2.1 水流模型

對(duì)于非均質(zhì)、空間三維非穩(wěn)定流系統(tǒng)，依據(jù)滲流的連續(xù)性方程和達(dá)西定律，遵從質(zhì)量守恒和能量轉(zhuǎn)換定律，各參數(shù)符合地下水流連續(xù)性方程及其定解條件方程式(公式1):

式中:Ω為滲流區(qū)域;h為含水層的水位標(biāo)高(m);S為儲(chǔ)水率(1/m);K為含水層水平方向的滲透系數(shù)(m/d);Kz為垂向滲透系數(shù)(m/d);Kn為邊界面法向方向的滲透系數(shù)(m/d);ε為源匯項(xiàng)(1/d)，包括開采量和壓縮釋水或儲(chǔ)水量;h0為初始水位(m);Γ1為通用水頭邊界;hb為通用水頭邊界上的水位;σ'為通用水頭邊界的阻力系數(shù)，σ'=L/K，L為模型邊界到通用水頭邊界的水平距離(m)，K為模型邊界到通用水頭邊界之間的平均滲透系數(shù)(m/d)，Γ2為上邊界;hp為潛水含水層水位(m);Kz為邊界內(nèi)側(cè)含水層垂向滲透系數(shù)(m/d);σ″為阻力系數(shù)，σ″=M'/K'，M 為弱透水層的厚度(m)，K'為弱透水層的垂向滲透系數(shù)(m/d);Γ3為下邊界;

2.2.2 土力學(xué)特征模型

假定地層總應(yīng)力不發(fā)生變化，且土層為垂向變形，水平方向的應(yīng)變甚微，可忽略不計(jì)。則粘土性孔隙水壓力的任何變化，有效應(yīng)力均跟著等值反方向變化。各參數(shù)遵循有效應(yīng)力公式(2):

式中:ρw為水的密度;g為重力加速度;Δh為水頭增量。

有效應(yīng)力增加或減小，使得粘性土層垂直方向的壓縮或回彈量呈線性增長(zhǎng)。土層彈性變形量遵從公式(3):

式中:Sskb為彈性骨架釋水系數(shù)，壓縮為正，回彈為負(fù);b0為壓縮單元的初始厚度。

當(dāng)粘土層的應(yīng)力水平小于前期最大應(yīng)力值，采用彈性骨架釋水系數(shù)(Sskb);當(dāng)粘土層所受的應(yīng)力水平超過前期最大應(yīng)力時(shí)，公式中Sskb則用一個(gè)非彈性骨架釋水系數(shù)代替(Sskv)。

3 西青區(qū)地面沉降數(shù)值模擬

3.1 數(shù)值模型建立

3.1.1 網(wǎng)格剖分

計(jì)算區(qū)內(nèi)層的劃分按含水組劃分5個(gè)層，含水組的厚度即為各層厚度，各層剖分采用相同的形式，主要是便于計(jì)算垂向各層之間的越流量。根據(jù)模擬區(qū)范圍，將模擬區(qū)平面上分為50行，50列，網(wǎng)格大小為200 m×200 m，總共2 500個(gè)單元格。

3.1.2 模擬期及初始條件

根據(jù)收集開采量和流場(chǎng)等數(shù)據(jù)，結(jié)合西青區(qū)分層標(biāo)的施工及監(jiān)測(cè)情況(2004年2月13號(hào)開鉆，2004年7月13日竣工，該標(biāo)于2005年1月開始監(jiān)測(cè))，擬選模型的模擬期為2005年7月—2010年6月。以月為單位，共劃分為60個(gè)應(yīng)力期，每個(gè)應(yīng)力期一個(gè)步長(zhǎng)。以2005年7月初的流場(chǎng)作為初始流場(chǎng)，以2010年8月的流場(chǎng)作為末擬合流場(chǎng)。

以2005年7月地下水水位統(tǒng)測(cè)數(shù)據(jù)，采用kriging插值法獲得各含水巖組的初始水位。

3.1.3 邊界條件

模型采用通用水頭邊界，邊界所在的網(wǎng)格的初始水位賦為通用水頭邊界的水位，參數(shù)C由經(jīng)驗(yàn)和流場(chǎng)給定，基本的計(jì)算方法是:

式中:L為單元格內(nèi)通用水頭邊界的長(zhǎng)度，K為滲透系數(shù)。

式中:Qb為通用水頭邊界上的流入或流出量。

3.1.4 源匯項(xiàng)的處理

模型外部源匯項(xiàng)主要是大氣降水入滲補(bǔ)給，蒸發(fā)，人工開采，其中入滲和蒸發(fā)通過recharge和蒸發(fā)包處理，人工開采在模型中都是以開采井的形式給入。

3.1.5 參數(shù)的選擇

參數(shù)選用是依據(jù)三維模型的計(jì)算需要，主要采取分區(qū)賦值的方法。通過收集的鉆孔資料，整理后的參數(shù)表如下(表3)。

表3 研究區(qū)水文地質(zhì)參數(shù)一覽表

3.1.6 力學(xué)參數(shù)的處理

在模擬地面沉降的模塊Interbed Storage中，需要輸入初始固結(jié)水頭，彈性儲(chǔ)水系數(shù)和非彈性儲(chǔ)水系數(shù)三個(gè)參數(shù)。

初始固結(jié)水頭:含水層的歷年最小水位值，模型中用初始水位代入。

彈性儲(chǔ)水系數(shù):經(jīng)驗(yàn)值，反應(yīng)土體骨架應(yīng)力小于前期最大應(yīng)力時(shí)，土體壓縮釋出的水量的參數(shù)，一般比非彈性儲(chǔ)水系數(shù)小，且從淺層到深層系數(shù)逐漸減小。根據(jù)其他地區(qū)的資料，模擬區(qū)各含水層的彈性儲(chǔ)水系數(shù)為5×10－6～2×10－5。

非彈性儲(chǔ)水系數(shù):經(jīng)驗(yàn)值，反應(yīng)由于土體骨架應(yīng)力小于前期最大應(yīng)力時(shí)，土體壓縮釋出的水量的參數(shù)，特點(diǎn)是越往深層，參數(shù)越小。根據(jù)前人資料和壓縮實(shí)驗(yàn)，模擬區(qū)各含水層的非彈性儲(chǔ)水系數(shù)為0.000 2～0.015。

3.2 模型模擬與驗(yàn)證

3.2.1 擬合方法

模型的識(shí)別和驗(yàn)證應(yīng)遵循以下原則:(1)模擬的地下水流場(chǎng)要與實(shí)際地下水流場(chǎng)基本一致，即要求地下水模擬等值線與實(shí)測(cè)地下水位等值線形狀相似;(2)模擬的地面沉降總趨勢(shì)與實(shí)際總地面沉降基本一致;(3)模擬地下水的動(dòng)態(tài)過程要與實(shí)測(cè)的動(dòng)態(tài)過程基本相似，即要求模擬與實(shí)際地下水位過程線形狀相似;(4)識(shí)別的水文地質(zhì)參數(shù)要符合實(shí)際水文地質(zhì)條件。

在此前提下，通過調(diào)整參數(shù)，保證模型輸出結(jié)果與實(shí)際觀測(cè)值接近;保證調(diào)整參數(shù)結(jié)果與實(shí)際水文地質(zhì)條件相符。

3.2.2 水位及流場(chǎng)擬合

在模擬時(shí)間范圍內(nèi)，將初始水位及其它要素輸入模型后運(yùn)行程序。觀測(cè)孔水位動(dòng)態(tài)曲線與模型模擬曲線擬合對(duì)比見圖1—圖4;第Ⅲ含水巖組擬合中期實(shí)際流場(chǎng)與模擬流場(chǎng)對(duì)比如圖5、圖6。

圖1 第Ⅰ含水組F11G1孔擬合對(duì)比

圖2 第Ⅱ含水組XI204孔擬合對(duì)比

圖3 第Ⅳ含水組XI401孔擬合對(duì)比

圖4 第Ⅴ含水組XI514孔擬合對(duì)比

圖5 第Ⅲ含水組2005年7月流場(chǎng)擬合情況(虛線為觀測(cè)，實(shí)線為模擬值)

圖6 第Ⅲ含水組2006年7月流場(chǎng)擬合情況(虛線為觀測(cè)，實(shí)線為模擬值)

3.2.3 地面沉降模擬

在模擬期內(nèi)，兼顧地下水滲流場(chǎng)變化，通過地下水與地面沉降的耦合關(guān)系，計(jì)算出的沉降量與實(shí)際沉降量進(jìn)行比較。以分層標(biāo)點(diǎn)為例，擬合分層標(biāo)監(jiān)測(cè)深度566 m內(nèi)單月沉降及五年內(nèi)累計(jì)沉降變化對(duì)比曲線，其模擬期內(nèi)累計(jì)沉降實(shí)際觀測(cè)值與模型計(jì)算值在一定程度上相符(見圖7、圖8)。

圖7 分層標(biāo)點(diǎn)單月沉降擬合情況

圖8 分層標(biāo)點(diǎn)累計(jì)沉降擬合情況

4 模擬結(jié)果分析

4.1 誤差分析

由于影響地下水動(dòng)態(tài)變化的因素非常復(fù)雜，具體模型已經(jīng)過高度概化，會(huì)與實(shí)際情況有一定誤差。產(chǎn)生差異的影響因素可能有:

1)某些觀測(cè)孔水位不能準(zhǔn)確代表所含水組水位，如在同一含水組內(nèi)，淺部水位與深部水位是有差異的。井孔濾水管如貫穿整個(gè)含水組，則代表性強(qiáng)，如僅安裝在局部含水層，則代表性減弱。

2)盡管十分注重常觀孔的觀測(cè)質(zhì)量，但可能存在觀測(cè)誤差。

3)長(zhǎng)觀及統(tǒng)測(cè)在觀測(cè)過程中受開采、降水、及人為因素影響，其數(shù)據(jù)存在偶然性。

4)因剖分的關(guān)系，計(jì)算水位代表的是整個(gè)單元的平均水位，而觀測(cè)水位是某個(gè)點(diǎn)的水位。

5)實(shí)測(cè)等水位線圖是人工繪制的流場(chǎng)，在繪制過程中存在誤差。

6)水文地質(zhì)概化產(chǎn)生的誤差及模擬過程的誤差。

4.2 三維水流模型與地面沉降對(duì)應(yīng)關(guān)系分析

通過本次數(shù)值模型方法計(jì)算沉降量與實(shí)際觀測(cè)值存在誤差，但整體能反映其沉降變化過程。分析模擬過程總結(jié)沉降量與地下水開采關(guān)系如下:

1)模擬計(jì)算沉降主要源于水位的變化量、地層彈性儲(chǔ)釋水系數(shù)、非彈性儲(chǔ)釋水系數(shù);水位變化是主導(dǎo)因素，即有水位變化便能計(jì)算出沉降量;數(shù)值模型計(jì)算不能反映開采引起沉降的滯后性。

2)除模擬區(qū)內(nèi)開采影響，邊界交換水量對(duì)區(qū)域內(nèi)水位變化影響也很明顯，從而對(duì)區(qū)域沉降影響也較明顯。

3)模擬計(jì)算沉降調(diào)參過程顯示，水位變化對(duì)沉降量影響較儲(chǔ)釋水系數(shù)影響大的多，儲(chǔ)釋水系數(shù)起到微調(diào)作用。

［1］陳崇希，裴順平.地下水開采—地面沉降模型研究［J］.水文地質(zhì)工程地質(zhì).2001(02).

［2］周載陽(yáng).地下水開采引起地面沉降的機(jī)理研究［J］.工程勘察.2012(03).

［3］楊勇，李國(guó)敏，等.抽取地下水引起地面沉降的研究現(xiàn)狀與進(jìn)展［J］.工程勘察.2010(11).

［4］張偉，段保旭，等.考慮地層結(jié)構(gòu)和地下水位變化時(shí)地面沉降的模擬和預(yù)測(cè)［J］.工程勘察.2002(05).

［5］姜規(guī)模，韓鳳霞.西安市地面沉降與地下水位動(dòng)態(tài)監(jiān)測(cè)信息管理系統(tǒng)應(yīng)用研究［J］.工程勘察.2010(06).

［6］楊勇，鄭凡東，等.北京平原區(qū)地下水水位與地面沉降關(guān)系研究［J］.工程勘察.2013(08).

［7］張本平.地面沉降與地下水位動(dòng)態(tài)監(jiān)測(cè)數(shù)據(jù)模型的研究［J］.山西建筑.2006(08).

［8］董克剛，王威，等.天津市地面沉降防治歷史的調(diào)查研究及啟示［J］.中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào).2008(03).

［9］王福，鐘新寶，等.天津市及其沿海地區(qū)地表高程變化現(xiàn)狀及趨勢(shì)［J］.地質(zhì)通報(bào).2005(01).

［10］王家兵，李平.天津平原地面沉降條件下的深層地下水資源組成［J］.水文地質(zhì)工程地質(zhì).2004(05).

［11］吳鐵鈞，金東錫.天津地面沉降防治措施及效果［J］.中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào).1998(02).

［12］Rao，N.Subba，Need for geological studies to probe land subsidence［J］.CURRENT SCIENCE，2012，102(3):382 － 382.

［13］張紅超，李琳，宋霞.回歸分析和灰色理論在地面沉降預(yù)測(cè)中精度的對(duì)比［J］.中小企業(yè)管理與科技(上旬刊).2011(05).

［14］朱恒銀，王幼鳳.地面沉降監(jiān)測(cè)標(biāo)結(jié)構(gòu)設(shè)計(jì)及施工技術(shù)要點(diǎn)探討［J］.探礦工程(巖土鉆掘工程).2003(S1).

［15］張阿根，顧為棟.上海市地面沉降監(jiān)測(cè)標(biāo)的設(shè)計(jì)原理與施工技術(shù)［J］.探礦工程(巖土鉆掘工程).2000(05).

［16］朱恒銀，王幼鳳.地面沉降監(jiān)測(cè)標(biāo)結(jié)構(gòu)設(shè)計(jì)及施工技術(shù)要點(diǎn)探討［J］.探礦工程(巖土鉆掘工程).2003(S1).

［17］龔士良，顧為棟.上海大口徑豎井式基巖水準(zhǔn)標(biāo)設(shè)計(jì)與施工［J］.探礦工程(巖土鉆掘工程).2002(S1).

［18］黃自培.天津市控沉鉆探現(xiàn)狀［J］.探礦工程(巖土鉆掘工程).1994(02).

［19］王達(dá).天津市將埋設(shè)基巖標(biāo)［J］.探礦工程(巖土鉆掘工程).1987(05).

［20］趙建康，吳孟杰，劉思秀.浙江省杭嘉湖平原地面沉降及監(jiān)測(cè)網(wǎng)絡(luò)研究［J］.地質(zhì)災(zāi)害與環(huán)境保護(hù).2004(01).