付延玲

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

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

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

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

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

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

2 數(shù)學模型

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

式中：

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

H為點（x，y，z）在t時刻的水頭值（m）；

W為源匯項（1／d）；SS為貯水率（1／m）；

t為時間（d）；Ω為立體計算域；

H0（x，y，z，t0）為點（x，y，z）處的初始水位（m）；

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

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

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

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

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

3 模型的識別、驗證

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

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

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

圖3 第Ⅱ承壓含水層2008年12月31日初始流場（單位：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ù)和儲水率等參數(shù)分區(qū)及其各分區(qū)的參數(shù)值。其中：潛水含水層共分6個參數(shù)區(qū)，水平滲透系數(shù)為3×10－6～1×102m／d，垂向滲透系數(shù)為7×10－7～40 m／d，給水度為1．54×10－5～3×10－1；淺部孔隙承壓含水層共分3個參數(shù)區(qū)，水平滲透系數(shù)為5×10－4～60m／d，垂向滲透系數(shù)為1×10－5～5m／d，彈性儲水率為1×10－4～2×10－3m－1；第I1、第I2承壓含水層共分19個參數(shù)區(qū)，水平滲透系數(shù)為5×10－1～25 m／d，垂向滲透系數(shù)為1×10－2～2m／d，彈性儲水率為1×10－8～2×10－5m－1；第II承壓含水層共分18個參數(shù)區(qū)，水平滲透系數(shù)為1×10－2～30m／d，垂向滲透系數(shù)為1×10－3～3m／d，彈性儲水率為5×10－8～3×10－4m－1；淺層黏性土弱含水層共分5個參數(shù)區(qū)，水平滲透系數(shù)為8×10－7～80m／d，垂向滲透系數(shù)為8×10－8～8m／d，彈性儲水率為1×10－8～2×10－3m－1；第I黏性土弱含水層（I1、I2間黏性土隔水層分布不穩(wěn)定）分為1個參數(shù)區(qū)，水平滲透系數(shù)為5．0×10－5m／d，垂向透系數(shù)為1．5×10－6m／d，彈性儲水率為1．0×10－6m－1；第II黏性土弱含水層也分為1個參數(shù)區(qū)，水平滲透系數(shù)為1．5×10－4m／d，垂向滲透系數(shù)為1．5×10－5m／d，彈性儲水率為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Ⅱ

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

4 地下水位動態(tài)預測

保持2009年的開采量和開采布局，以2009年12月31日作為預測初始時刻，預測2009年12月31日至2020年12月31日逐月的地下水流場變化特征。

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

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

5 地面沉降動態(tài)預測

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

根據(jù)太沙基有效應力原理，土體壓縮變形量與孔隙水壓力的變化量存在線性相關關系［9］，由此可以推斷含水層的壓縮沉降量與地下水位變幅存在線性相關關系［10］。根據(jù)寧波市江東北路監(jiān)測點1985年至2008年各含水層及之間的黏性土弱含水層的月壓縮沉降量與相應的月水位變幅數(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日預測流場（單位：m）Fig．8 Forecast flow fields of confined aquiferⅡon December 31st，2020（unit：m）

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

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

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

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

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

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

（2）地面沉降量預測

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

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

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

6 結(jié)論

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

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

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

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

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

［4］孫從軍，韓振波，趙振，等．地下水數(shù)值模擬的研究與應用進展［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ù)試驗及工程應用［J］．寧波大學學報，2013，26（3）：127－132．

［8］許燁霜，余恕國，沈水龍，等．地下水開采引起地面沉降預測方法的現(xiàn)狀與未來［J］．防災減災工程學報，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］賈瑩媛，黃張裕，張蒙，等．含水組地下水位變化對地面沉降影響的多元回歸分析與預測［J］．工程勘察，2013（1）：77－80．

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