寧波市地下水位動(dòng)態(tài)與地面沉降預(yù)測(cè)分析＊

付延玲

（河海大學(xué)地球科學(xué)與工程學(xué)院，南京 210098）

寧波市是我國(guó)著名的港口城市和歷史文化名城，位于浙江省東部，長(zhǎng)江三角洲東南角，浙江寧紹平原東端，東臨東海，北瀕杭州灣，西接紹興市，南接三門灣，東西寬175km，南北長(zhǎng)192km。區(qū)內(nèi)從上往下分布的第四紀(jì)松散孔隙含水層有全新統(tǒng)的潛水含水層、上更新統(tǒng)第Ⅰ承壓含水層（可劃分為上下兩個(gè)含水層I1、I2）和中更新統(tǒng)第II承壓含水層，各含水層之間均以弱含水的黏性土層相分隔，是我國(guó)第四紀(jì)孔隙承壓水開發(fā)利用最早的城市之一。上世紀(jì)六十年代以來，隨著工農(nóng)業(yè)的快速發(fā)展，深層孔隙承壓水開采量劇增，造成區(qū)域地下水位急劇下降，地下水資源衰竭、水質(zhì)惡化，引發(fā)了嚴(yán)重的地面沉降問題，并有進(jìn)一步擴(kuò)大的趨勢(shì)。為此，寧波市政府為了有效地控制地面沉降的進(jìn)一步發(fā)展，從2008年底起對(duì)第四紀(jì)孔隙承壓水實(shí)行全面禁采，孔隙承壓水水位得到了有效控制，并日趨回升。

為了查明禁采后地下水位的動(dòng)態(tài)變化趨勢(shì)及對(duì)地面沉降的控制效應(yīng)，本文將寧波市第四紀(jì)松散沉積層作為一個(gè)統(tǒng)一的水文地質(zhì)體，在概化出寧波市第四紀(jì)松散沉積層地下水系統(tǒng)水文地質(zhì)概念模型的基礎(chǔ)上，建立了寧波市第四紀(jì)松散沉積層地下水系統(tǒng)三維數(shù)值模擬模型和地面沉降與地下水位多元線性回歸模型［1－2］，在對(duì)各含水層及相應(yīng)黏性土弱含水層地下水位進(jìn)行模擬預(yù)測(cè)的同時(shí)，模擬預(yù)測(cè)了由各含水層及相應(yīng)黏性土弱含水層地下水位變化引起的地面沉降問題，為寧波市地面沉降防治規(guī)劃的確定提供了科學(xué)依據(jù)。

1 研究區(qū)地下水系統(tǒng)水文地質(zhì)概念模型

本次模擬計(jì)算將寧波市第四紀(jì)松散沉積層作為一個(gè)統(tǒng)一的地下水系統(tǒng)，平面上北部以江北區(qū)與鎮(zhèn)海區(qū)行政交接邊界為界，東部以鄞州區(qū)與北侖區(qū)行政交接邊界為界，東南部以第四系邊界為界，西南部以鄞州區(qū)與奉化市行政交接邊界為界，西部以第四系邊界為界，四周均概化為通用水頭邊界，計(jì)算面積共約835km2，包括了海曙區(qū)、江東區(qū)、江北區(qū)全區(qū)及鄞州區(qū)平原部分。剖面上自上往下，包括潛水含水層、淺部承壓含水層、第I1、第I2承壓含水層和第II承壓含水層，以及各含水層之間的黏性土弱含水層，各層均概化為非均質(zhì)各向異性（圖1）。系統(tǒng)頂部一方面接受大氣降水的補(bǔ)給，為補(bǔ)給邊界；另一方面地下水又通過其蒸發(fā)及植物的蒸騰，是一排泄邊界。系統(tǒng)的底部接受下伏基巖水的補(bǔ)給，是一補(bǔ)給邊界。地下水流態(tài)概化為三維非穩(wěn)定流。區(qū)內(nèi)地下水的開采量根據(jù)統(tǒng)計(jì)資料，分別按單井及以鄉(xiāng)鎮(zhèn)為單位的大井處理。

圖1 寧波市水文地質(zhì)剖面圖Fig．1 Hydrogeological section of Ningbo City

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

根據(jù)上述研究區(qū)地下水系統(tǒng)的水文地質(zhì)概念模型，取坐標(biāo)軸方向和含水層主滲透方向相一致，可建立下列與之相適應(yīng)的地下水運(yùn)動(dòng)三維非穩(wěn)定流數(shù)學(xué)模型［3］

式中：

kxx、kyy，、kzz分別為含水層各向異性主方向滲透系數(shù)（m／d）；

H為點(diǎn)（x，y，z）在t時(shí)刻的水頭值（m）；

W為源匯項(xiàng)（1／d）；SS為貯水率（1／m）；

t為時(shí)間（d）；Ω為立體計(jì)算域；

H0（x，y，z，t0）為點(diǎn)（x，y，z）處的初始水位（m）；

cos （n，x） 、cos （n，y） 、cos （n，z） 分別為各向流量邊界外法線方向與坐標(biāo)軸方向夾角的余弦；

q（x，y，z，t）為第二類邊界上單位面積的側(cè)向補(bǔ)給量（m／d）；

μ為飽和差（自由面上升）或給水度（自由面下降），它表示在自由面改變單位高度下，從含水層單位截面積上吸收或排出的水量；

qw為自由面單位面積上的大氣降雨入滲補(bǔ)給量與地下水蒸發(fā)量之和（m／d）；

Γ2、Γ3分別為第二類邊界和自由面邊界。

3 模型的識(shí)別、驗(yàn)證

上述模型采用有限差分法求解［4－5］，并采用強(qiáng)隱式（SIP）聯(lián)立迭代求解代數(shù)方程組［6］。將整個(gè)研究區(qū)在平面上剖分成250×250個(gè)矩形網(wǎng)格單元，單元邊長(zhǎng)為270×268m。垂向上按潛水含水層、淺部孔隙承壓含水層、第I1、第I2承壓含水層、第II承壓含水層及各含水層之間的黏性土弱含水層進(jìn)行剖分，共分9層。每層的有效計(jì)算單元為11977個(gè)，共計(jì)107793個(gè)，詳見圖2。以2008年12月31日至2009年12月31日作為模型識(shí)別、驗(yàn)證的時(shí)段，一個(gè)月為一個(gè)應(yīng)力期，共分為12個(gè)應(yīng)力期，每個(gè)應(yīng)力期又分為10個(gè)計(jì)算時(shí)間步長(zhǎng)。

地下水系統(tǒng)各含水層均有一定數(shù)量的地下水位監(jiān)測(cè)井用來進(jìn)行水位擬合，共計(jì)76個(gè)，基本上控制了各含水層地下水流場(chǎng)。各含水層的初始流場(chǎng)均由實(shí)測(cè)水位經(jīng)Kriging插值給出，各含水層之間的黏性土弱含水層初始流場(chǎng)由上下含水層初始流場(chǎng)插值獲得，圖3為第II承壓含水層初始流場(chǎng)圖。2009年各含水層開采量由實(shí)際統(tǒng)計(jì)給出。各含水層及之間的黏性土弱含水層通用水頭邊界上的水頭值由實(shí)測(cè)值經(jīng)插值給出。各含水層及之間的黏性土弱含水層參數(shù)分區(qū)的參數(shù)初值均按實(shí)驗(yàn)分析、前人資料并結(jié)合經(jīng)驗(yàn)給出［7－8］。圖4和圖5為第II承壓含水層各監(jiān)測(cè)井2009年6月30日和12月30日的地下水位擬合精度，水位擬合誤差均在1m 以下。

圖2 研究區(qū)剖面圖Fig．2 Section maps of studied area

圖3 第Ⅱ承壓含水層2008年12月31日初始流場(chǎng)（單位：m）Fig．3 Initial flow fields of confined aquiferⅡon December 31st，2008

圖4 第Ⅱ承壓含水層2009年6月30日地下水位擬合圖Fig．4 Fitting of groundwater levels of confined aquifer Ⅱon June 30st，2009

經(jīng)識(shí)別、驗(yàn)證，得出了各層的滲透系數(shù)和儲(chǔ)水率等參數(shù)分區(qū)及其各分區(qū)的參數(shù)值。其中：潛水含水層共分6個(gè)參數(shù)區(qū)，水平滲透系數(shù)為3×10－6～1×102m／d，垂向滲透系數(shù)為7×10－7～40 m／d，給水度為1．54×10－5～3×10－1；淺部孔隙承壓含水層共分3個(gè)參數(shù)區(qū)，水平滲透系數(shù)為5×10－4～60m／d，垂向滲透系數(shù)為1×10－5～5m／d，彈性儲(chǔ)水率為1×10－4～2×10－3m－1；第I1、第I2承壓含水層共分19個(gè)參數(shù)區(qū)，水平滲透系數(shù)為5×10－1～25 m／d，垂向滲透系數(shù)為1×10－2～2m／d，彈性儲(chǔ)水率為1×10－8～2×10－5m－1；第II承壓含水層共分18個(gè)參數(shù)區(qū)，水平滲透系數(shù)為1×10－2～30m／d，垂向滲透系數(shù)為1×10－3～3m／d，彈性儲(chǔ)水率為5×10－8～3×10－4m－1；淺層黏性土弱含水層共分5個(gè)參數(shù)區(qū)，水平滲透系數(shù)為8×10－7～80m／d，垂向滲透系數(shù)為8×10－8～8m／d，彈性儲(chǔ)水率為1×10－8～2×10－3m－1；第I黏性土弱含水層（I1、I2間黏性土隔水層分布不穩(wěn)定）分為1個(gè)參數(shù)區(qū)，水平滲透系數(shù)為5．0×10－5m／d，垂向透系數(shù)為1．5×10－6m／d，彈性儲(chǔ)水率為1．0×10－6m－1；第II黏性土弱含水層也分為1個(gè)參數(shù)區(qū)，水平滲透系數(shù)為1．5×10－4m／d，垂向滲透系數(shù)為1．5×10－5m／d，彈性儲(chǔ)水率為1．0×10－8m－1。圖6為第II承壓含水層水文地質(zhì)參數(shù)分區(qū)特征，表1為第II承壓含水層的水文地質(zhì)參數(shù)值。

圖5 第Ⅱ承壓含水層2009年12月30日地下水位擬合圖Fig．5 Fitting of groundwater levels of confined aquifer Ⅱon December 30st，2009

圖6 第Ⅱ承壓含水層水文地質(zhì)參數(shù)分區(qū)圖Fig．6 Division map of hydrogeological parameters of confined aquiferⅡ

表1 第Ⅱ承壓含水層水文地質(zhì)參數(shù)分區(qū)值Table 1 Division values of hydrogeological parameters of confined aquiferⅡ

從識(shí)別結(jié)果來看，地下水系統(tǒng)各參數(shù)分區(qū)參數(shù)值的級(jí)別大小均符合常規(guī)，邊界上的水量交換強(qiáng)弱程度和實(shí)際基本一致，較好地反映了寧波市第四紀(jì)松散沉積層地下水系統(tǒng)的結(jié)構(gòu)與功能特征。各監(jiān)測(cè)井水位計(jì)算值和實(shí)測(cè)值擬合程度也較好，說明模型計(jì)算所得結(jié)果正確、可信，可以用來預(yù)測(cè)寧波市第四紀(jì)松散沉積層地下水位的動(dòng)態(tài)變化。

4 地下水位動(dòng)態(tài)預(yù)測(cè)

保持2009年的開采量和開采布局，以2009年12月31日作為預(yù)測(cè)初始時(shí)刻，預(yù)測(cè)2009年12月31日至2020年12月31日逐月的地下水流場(chǎng)變化特征。

預(yù)測(cè)結(jié)果表明：山區(qū)溝谷孔隙潛水地下水位降落漏斗逐漸擴(kuò)大，至2020年12月底，如鳳岙水廠開采形成的地下水位降落漏斗中心水位標(biāo)高為－3．4m，鄞江水廠、鎮(zhèn)電集水廠開采形成的地下水位降落漏斗中心水位標(biāo)高為0．0m。平原區(qū)潛水水位年內(nèi)變化隨大氣降水量變化而變化，總體上水位基本處于穩(wěn)定；山前淺部孔隙承壓水動(dòng)態(tài)變化與潛水基本一致，平原區(qū)淺部孔隙承壓含水層水位年際變化不大，地下水流場(chǎng)年際變化不明顯；全區(qū)深層孔隙承壓水水位降落漏斗恢復(fù)較快，第I1與第I2承壓含水層地下水位降落漏斗在2010年4月之后即已恢復(fù)，2010年6月第II承壓含水層降落漏斗恢復(fù)。深層孔隙承壓水水位年際變化不大，年內(nèi)水位一般4月到8月份水位較低，9月到次年3月水位較高。圖7 和圖8為第II承壓含水層2012年12月31日和2020年12月31日的地下水預(yù)測(cè)流場(chǎng)，可見，2012年后深層孔隙承壓水的地下水流場(chǎng)基本穩(wěn)定，年際變化不大，第I、第II承壓含水層位于市區(qū)的地下水位降落漏斗也均已得到恢復(fù)。

圖7 第Ⅱ承壓含水層2012年12月31日預(yù)測(cè)流場(chǎng)（單位：m）Fig．7 Forecast flow fields of confined aquiferⅡon December 31st，2012（unit：m）

5 地面沉降動(dòng)態(tài)預(yù)測(cè)

（1）地面沉降量與地下水位變化的多元線性回歸模型

根據(jù)太沙基有效應(yīng)力原理，土體壓縮變形量與孔隙水壓力的變化量存在線性相關(guān)關(guān)系［9］，由此可以推斷含水層的壓縮沉降量與地下水位變幅存在線性相關(guān)關(guān)系［10］。根據(jù)寧波市江東北路監(jiān)測(cè)點(diǎn)1985年至2008年各含水層及之間的黏性土弱含水層的月壓縮沉降量與相應(yīng)的月水位變幅數(shù)據(jù)，采用多元線性回歸方法，建立地面沉降量與各含水層及之間的黏性土弱含水層地下水位變幅的多元線性回歸模型：

Δb＝a0＋a1ΔH1＋a2ΔH2＋a3ΔH3＋a4ΔH4＋a5ΔH5＋a6ΔH6

其中：

Δb為月地面沉降量（mm）；

ΔH1為潛水月水位變幅（m）；

ΔH2為淺部孔隙承壓含水層月水位變幅（m）；

ΔH3為淺部孔隙承壓含水層與第Ⅰ承壓含水層間黏性土弱含水層月水位變幅（m）；

ΔH4為第Ⅰ承壓含水層月水位變幅（m）；

ΔH5為第Ⅰ承壓含水層與第II承壓含水層間黏性土弱含水層月水位變幅（m）；

圖8 第Ⅱ承壓含水層2020年12月31日預(yù)測(cè)流場(chǎng)（單位：m）Fig．8 Forecast flow fields of confined aquiferⅡon December 31st，2020（unit：m）

ΔH6為第II承壓含水層月水位變幅（m）；

a0為常數(shù)項(xiàng)，與各含水層和黏性土弱含水層的監(jiān)測(cè)控制程度和監(jiān)測(cè)精度及非地下水開采引起的地面沉降有關(guān)；

a1、a2、a3、a4、a5、a6為系數(shù)項(xiàng)，與各含水層和黏性土弱含水層的骨架釋水系數(shù)有關(guān)。

根據(jù)已有266 組因變量（月沉降量）與自變量（各層月水位變幅）觀測(cè)統(tǒng)計(jì)數(shù)據(jù)，經(jīng)計(jì)算得到了下列多元線性回歸方程：

Δb＝－0．856＋1．222ΔH1＋0．457ΔH2＋2．605ΔH3＋0．319ΔH4＋0．626ΔH5＋0．832ΔH6

上述模型經(jīng)擬合度檢驗(yàn)分析、方差分析，共線性檢驗(yàn)分析和殘差分析［11］可知，模型正確可信，可以用來預(yù)測(cè)地面沉降的發(fā)展趨勢(shì)。

（2）地面沉降量預(yù)測(cè)

應(yīng)用上述地下水系統(tǒng)數(shù)值模型，預(yù)測(cè)得到了江東北路監(jiān)測(cè)點(diǎn)各含水層及之間的黏性土弱含水層2010至2020年的月水位變幅值，再利用上述地面沉降與地下水位變化的多元線性回歸模型，預(yù)測(cè)得到了逐月的地面沉降量，在此基礎(chǔ)上，經(jīng)累計(jì)獲得了2010年至2020年的年地面沉降量，詳見表2。

表2 2010年至2020年地面沉降量預(yù)測(cè)（mm／a）Table 2 Prediction of land subsidence rates from 2010to 2020（mm／a）

由上述計(jì)算結(jié)果可知，寧波市江東北路監(jiān)測(cè)點(diǎn)2012年后年地面沉降量逐年減小，主要是因?yàn)?012年后各含水層及之間的黏性土弱含水層地下水位趨于穩(wěn)定，由地下水位下降引起的地面沉降均基本得到控制。

6 結(jié)論

（1）本次模擬計(jì)算將寧波市第四紀(jì)松散沉積層作為一個(gè)統(tǒng)一的地下水系統(tǒng)，采用真三維數(shù)值模型，分別建立了各含水層和黏性土弱含水層的地下水流動(dòng)方程，并通過垂向滲透系數(shù)發(fā)生水力聯(lián)系，克服了以往二維和準(zhǔn)三維模型將相鄰含水層之間的黏性土弱含水層概化為越流層給模擬計(jì)算帶來重復(fù)或缺失的不足，提高了模擬計(jì)算的精度。

（2）寧波市2012年后各含水層及之間的黏性土弱含水層地下水位趨于穩(wěn)定，地面沉降量逐年減小，由地下水位下降引起的地面沉降均基本得到控制。

［1］薛禹群，謝春紅．地下水?dāng)?shù)值模擬［M］．北京：科學(xué)出版社，2007：1－50．

［2］黃丹，肖偉，李勇．地下水三維數(shù)值模擬及其優(yōu)化開采［J］．資源調(diào)查與環(huán)境，2005，26（2）：137－145．

［3］付延玲．基于地面沉降控制的區(qū)域性松散沉積層地下水可采資源規(guī)劃評(píng)價(jià)［J］．吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），2012．42（2）：476－484．

［4］孫從軍，韓振波，趙振，等．地下水?dāng)?shù)值模擬的研究與應(yīng)用進(jìn)展［J］．環(huán)境工程，2013，31（5）：9－17．

［5］Wang Kun，Jin Scheng，Liu Gang．Numerical modeling of free－surface flows with bottom and surface－layer pressure treatment［J］．Journal of Hydrodynamics，2009，21（3）：352－358．

［6］PATRICKL A．Natural disasters［J］．Surficial Geology，1996，21（3）：34－36．

［7］吳天壽，葉俊能．寧波軌道交通水文地質(zhì)參數(shù)試驗(yàn)及工程應(yīng)用［J］．寧波大學(xué)學(xué)報(bào)，2013，26（3）：127－132．

［8］許燁霜，余恕國(guó)，沈水龍，等．地下水開采引起地面沉降預(yù)測(cè)方法的現(xiàn)狀與未來［J］．防災(zāi)減災(zāi)工程學(xué)報(bào)，2006，26（3）：352－356．

［9］ZHANG Yun，XUE Yu－qun．WU Ji－chun，et al．Land subsidence and earth fissures due to groundwater with－drawl in the Southern Yangtze Delta，China［J］．Environmental Geology 2008，55（4）：751－762．

［10］賈瑩媛，黃張?jiān)?，張蒙，等．含水組地下水位變化對(duì)地面沉降影響的多元回歸分析與預(yù)測(cè)［J］．工程勘察，2013（1）：77－80．

［11］袁志發(fā)，周靜芋．多元統(tǒng)計(jì)分析［M］．北京：科學(xué)出版社，2002．