滄州市地下水開采-地面沉降數值模擬

馬青山,駱祖江

(河海大學地球科學與工程學院,江蘇 南京 210098)

?

滄州市地下水開采-地面沉降數值模擬

馬青山,駱祖江

(河海大學地球科學與工程學院,江蘇 南京 210098)

以滄州市地下水開采引起的地面沉降為例,在分析研究區地質和水文地質條件的基礎上,得到研究區地下水開采與地面沉降地質概化模型,建立地下水三維滲流和土體垂向一維變形耦合數值模型,在完成模型識別驗證的基礎上,預測了現狀開采條件下從2010年末到2025年末地下水流場的變化特征和地面沉降的變化趨勢。結果表明,到2025年底,整個滄州市范圍內累計地面沉降量最大值為475 mm,最小值為270 mm。

地下水開采;地面沉降;預測模型;耦合數值模型

滄州市是我國第四系孔隙承壓水開發利用最早、開采量最大的城市之一,特別是20世紀80年代以來,隨著工農業的快速發展,地下水開采量大幅增加,造成了區域地下水位急劇下降,地下水資源衰竭,誘發了嚴重的地面沉降問題。因此建立區域性的地下水開采與地面沉降模擬模型,準確預測地面沉降發展趨勢已成為當務之急。目前,國內外關于地下水滲流與地面沉降問題已經提出了許多計算模型,主要有基于達西定律和太沙基一維固結理論的地下水滲流與地面沉降耦合模型及基于比奧固結理論的全耦合模型。由于全耦合模型存在諸多缺點,如，計算參數多、收斂速度慢、計算耗時長、模型難以穩定等,所以實用性較差[1-2]。本文嘗試采用地下水滲流與地面沉降耦合模型對滄州市地下水位和地面沉降的動態變化進行模擬和預測，旨在為滄州市地面沉降的科學防控提供決策依據。

1 模型基礎

1.1 地下水開采與地面沉降地質概念模型

1.1.1 水文地質概念模型

滄州市沉積有巨厚、復雜的第四系松散沉積層,厚度一般為350～550 m,厚者可達550～600 m。根據地層巖性和時代,第四系松散沉積層自上而下可以劃分為全新統第Ⅰ含水層組、上更新統第Ⅱ含水層組、中更新統第Ⅲ含水層組、下更新統第Ⅳ含水層組[2]。考慮到含水層中亞黏土、黏土透鏡體分布的不均勻性和不連續性,本文不單獨劃出黏土層作為含水層之間的弱透水層,而是將黏土透鏡體和含水層看作為一個整體的地下水系統。各含水層之間的垂向水力聯系復雜,地下水開采直接影響各層地下水位,故模型中將地下水流態概化為三維非穩定流。研究區的四周概化為通用水頭邊界。含水層的頂部一方面接受大氣降水入滲補給,是一補給邊界,另一方面地下水又通過其蒸發,是一排泄邊界；系統底部為隔水底板,與系統不存在水力聯系,可作隔水邊界[3],如圖1。

圖1 研究區三維模型示意圖

1.1.2 土體變形的概化

根據Terzaghi有效應力原理,在假定研究區土層總應力保持不變的情況下,孔隙水壓力的減小量等于骨架有效應力的增加量[4]。

嚴格地講,土體變形都是三維的,但考慮到研究區范圍廣,地面水平位移十分微小,所以本文忽略了土層的水平變形[3],認為土體的變形只發生在垂向上,即土體變形是垂向一維的。

1.2 地下水滲流與地面沉降耦合數學模型

1.2.1 水流數學模型

依據研究區水文地質概念模型,其數學模型可以描述為取多孔介質各向異性主滲透方向與坐標軸方向一致[5-6],即

(1)

(2)

(3)

(4)

1.2.2 土體垂向一維變形模型

由地下水位下降引起的含水層變形量的計算模型[7]如下式所示。

承壓含水層的彈性變形量:

(5)

承壓含水層的非彈性變形量:

(6)

潛水含水層彈性變形量:

(7)

潛水含水層非彈性變形量:

(8)

其中

(9)

(10)

在地面沉降模型中,含水層及包含在其中的弱透水層被視為一個整體,假設整個含水層的骨架都可發生壓縮并產生沉降作用,參數Sfe、Sfv是在整個含水層規模上的等效參數[9],其參數值經過模型調參反演得到。等效原理:假設n層含水層的儲水率分別為Ss1，Ss2，…，Ssn,厚度分別為b1，b2,…,bn,并使用了一個等效儲水因子Ssystem：

(11)

當各含水層的變形量確定后,總的地面沉降量便可由各含水層變形量求和而得出。將水流模型與土體垂向一維變形模型通過水頭項耦合起來,便構成了地下水開采與地面沉降耦合模擬模型。

圖2 研究區網格剖分圖

2 模型的識別驗證

研究區東西長約197.5 km,南北寬約179 km,總面積14 056 km2,垂向厚度300～500 m,厚者可達550～600 m。將整個研究區在平面上剖分為150×150的矩形網格單元,在垂向上,第Ⅰ含水層組、第Ⅱ含水層組、第Ⅲ含水層組、第Ⅳ含水層組4個含水層均作為獨立層位進行計算。網格數共計22 500個,其中有效單元9 365個,無效單元13 135個。其平面和垂向剖分網格見圖2。以2008年12月31日至2010年12月31日作為模型的識別、驗證階段。模型一共分為24個應力期,每個應力期分為10個時間步長。

本次模擬的過程中,各含水層組均設有一定數量的觀測井來觀測整個含水層組的水位變化并用來進行水位擬合,其中第Ⅰ含水層組28個觀測井和第Ⅱ含水層組3個觀測井,第Ⅲ含水層組19個觀測井,第Ⅳ含水層組2個觀測井,共52個觀測井,基本控制全區。采用2008年12月31日各含水層組的水位作為模型的初始水位值。第Ⅰ和第Ⅲ含水層組的初始水位為已知,第Ⅱ和第Ⅳ含水層組的初始水位由相鄰含水層的水位線性插值得到。考慮到初始沉降量對模擬計算影響很小以及地面沉降監測資料的缺乏,故將計算所需的初始沉降量統一取零。各縣市各含水層地下水的開采量經實際調查獲得。模型所需參數初值均由以往資料結合現場試驗和室內試驗給出。在此基礎上，根據抽水試驗資料,對觀測井的水位進行水位擬合,反演計算出各地的水文地質參數[11-12]。以第Ⅲ含水層組為例,部分觀測井分布如圖3所示,第Ⅲ含水層組參數分區如圖4所示,圖5給出了部分觀測井的計算水位與實測水位擬合情況。圖6給出了2008年12月31日第Ⅲ含水層組地下水位實測和計算等值線。圖7給出了2008年12月31日至2010年12月31日期間累計的地面沉降量實測和計算等值線圖。表1舉例說明了第Ⅲ含水層組的有關參數分區及各分區的參數值。

圖3 部分觀測井分布示意圖

圖4 第Ⅲ含水層組參數分區

圖5 第Ⅲ含水層組觀測井地下水位擬合

表1 第Ⅲ含水層組各參數分區參數

分區號主方向滲透系數/(m·d-1)KxxKyyKzz儲水率Ss/(1·m-1)骨架成分彈性釋水率Sfe213 53 54 00×10-76 50×10-54 63×10-5224 54 54 50×10-62 55×10-51 75×10-52330 030 08 40×10-65 00×10-54 76×10-52420 020 03 00×10-61 60×10-41 07×10-42515 015 02 50×10-51 00×10-58 23×10-6264 54 55 30×10-51 02×10-59 43×10-6277 57 54 00×10-53 50×10-52 96×10-5288 08 03 50×10-63 00×10-52 56×10-5293 53 51 02×10-52 30×10-52 10×10-5302 02 05 10×10-41 50×10-51 38×10-5314 04 08 00×10-63 10×10-62 52×10-6326 06 05 00×10-64 00×10-53 78×10-5335 55 54 50×10-53 00×10-62 74×10-6346 56 53 60×10-54 00×10-63 66×10-6

圖6 2008年12月31日第Ⅲ含水層組地下水位實測和計算等值線(單位:m)

圖7 2009—2010年累計地面沉降量等值線(單位:mm)

從圖5可以看出模型的地下水位擬合精度較好。從圖6可以看出,地下水位計算與實際等值線基本相吻合。從圖7可以看出,地面沉降整體擬合出了沉降趨勢。綜上所述,該模型可以用來模擬預測滄州市地面沉降。

3 地下水位與地面沉降預測

圖8 第Ⅲ含水層組預測地下水位等值線(單位:m)

保持2010年的地下水開采布局及開采量不變,利用已經識別和驗證的地下水開采與地面沉降耦合模型,以2010年12月31日為預測的初始時刻,模擬計算各含水層組2010年12月31日至2025年12月31日逐年的地下水水位變化特征以及地面沉降變化趨勢。第Ⅲ含水層組是滄州市主要開采地下水的含水層,因此以第Ⅲ含水層組為例來說明地下水水位變化特征以及地面沉降變化趨勢。

3.1 地下水水位預測

從模型運行得出的結果可以看出,第Ⅲ含水層組是滄州市主要開采地下水的含水層,地下水水位下降較快。2010年12月31日和2025年12月31日,各縣市第Ⅲ含水層最低地下水水位變化情況如表2所示。

表2 各地區第Ⅲ含水層組最低地下水水位 m

從表2可以得出,到2025年12月31日,第Ⅲ含水層組各地區地下水位年平均最大下降速率為:任丘1.22 m/a;河間2.06 m/a;肅寧1.03 m/a;獻縣2.13 m/a;泊頭2.33 m/a;青縣2.02 m/a;南皮1.21 m/a;東光1.95 m/a;吳橋1.21 m/a;孟村0.93 m/a;鹽山2.02 m/a;黃驊1.75 m/a;海興0.98 m/a;滄州市區1.45 m/a。第Ⅲ含水層組地下水位等值線如圖8所示。

3.2 地面沉降預測

由于地下水位持續下降,截止2025年12月31日,滄州地區肅寧縣、任丘市、河間市、獻縣、泊頭市、青縣、南皮縣、東光縣、吳橋縣、黃驊市、鹽山縣、孟村縣、海興縣、滄州市區最大地面沉降量分別為475 mm、310 mm、408 mm、380 mm、300 mm、259 mm、305 mm、337.5 mm、336 mm、235 mm、298 mm、300 mm、255 mm、270 mm,最大沉降速率分別為31.67mm/a、20.67mm/a、27.20 mm/a、25.33mm/a、20.00 mm/a、17.27 mm/a、20.33 mm/a、22.5 mm/a、22.40 mm/a、15.67 mm/a、19.87 mm/a、20.00 mm/a、17.00 mm/a、18.00 mm/a。地面沉降量預測等值線見圖9。

圖9 地面沉降預測等值線 (單位:mm)

通過圖9可以看出,在地下水現狀開采條件下,滄州形成了任丘市、肅寧縣、孟村縣、鹽山縣、吳橋縣、黃驊市6個較大的沉降漏斗分布區,其中以肅寧縣沉降漏斗最為明顯,最大沉降速率達到31.67 mm/a。隨著地下水的逐年開采,沉降漏斗范圍有不斷擴大的趨勢。

4 結 語

在充分研究所收集和整理的滄州市水文地質和工程地質資料的基礎上,依據地下水三維滲流和Terzaghi垂向一維固結理論,建立滄州市地下水三維滲流與地面沉降耦合數值模型,并利用該模型預測了滄州市地下水在現有開采條件下從2010年底至2025年底逐年的地下水位變化特征及地面沉降發展趨勢。主要有以下結論。

a.依據Terzaghi有效應力原理,結合地下水三維滲流理論,以此建立的地下水三維非穩定流和土體一維垂向變形耦合模型，是滄州市較為完善的能夠正確刻畫出研究區實際水文地質特征的模型。

b.根據收集的實際資料,選用部分耦合模型對滄州市地面沉降進行研究,通過含水層的儲水性和土體的壓縮性將地下水滲流和土體變形模型耦合,模型再現了地下水水位和地面沉降的動態特征,反映出在開采條件下各含水層水位動態及地面沉降變化規律,能為地下水資源管理以及地面沉降防治工作提供科學依據。

[1] 陳興賢,駱祖江,安曉宇,等.深基坑降水三維變參數非穩定滲流與地面沉降耦合模型[J].吉林大學學報:地球科學版,2013,43(5):1574-1578.(CHEN Xinxian,LUO Zujiang,AN Xiaoyu,et al.Coupling model of groundwater three dimensional variable parametric non-steady seepage and land-subsidence [J].Journal of Jilin University:Earth Science Edition,2013,43(5):1574-1578.(in Chinese))

[2]李培超.地面沉降變形非線性完全耦合數學模型[J].河海大學學報:自然科學版,2011,39(6):665-670.LI Peichao.A nonlinear fully coupled mathematical model for land subsidence[J].Journal of Hohai University:Natural Sciences,2011,39(6):665-670.(in Chinese))

[3] 駱祖江,王琰,田小偉,等.滄州市地下水開采與地面沉降地裂縫模擬預測[J].水利學報,2013,44(2):198-204.(LUO Zujiang,WANG Yan,TIAN Xiaowei,et al.Simulating and forecasting of ground water exploitation,land subsidence and ground fissure in Cangzhou City [J].Journal of Hydraulic Engineering,2013,44 (2):198-204.(in Chinese))

[4] 萬偉峰.西安市地下水開采-地面沉降數值模擬及防治方案研究[D].西安:長安大學,2008.

[5] 駱祖江,武永霞.鹽城市地下水資源規劃評價三維數值模型[J].水資源保護,2005,21(5):37-41.(LUO Zujiang,WU Yongxia.Three dimensional numerical model for groundwater resource evaluation and plan in Yancheng [J].Water Resources Protection,2005,21 (5):34-41.(in Chinese))

[6] 付延玲.鹽城市地下水資源預測評價[J].資源調查與環境,2005,26(1):54-59.(FU Yanling.Evaluation of groundwater resource in Yancheng[J].Resource Survey and Environment,2005,26(1):54-59.(in Chinese))

[7] 駱祖江,張月萍,劉金寶,等,江蘇沿江開發帶地下水開采與地面沉降三維數值模擬[J].地球科學與環境學報,2007,29(3):281-284.(LUO Zujiang,ZHANG Yueping,LIU Jinbao,et al.Three-dimensional numerical model for groundwater withdrawal and land-subsidence in development zone along Yangtze River in Jiangsu Province[J].Journal of Jilin University:Earth Science Edition,2007,29(3):281-284.(in Chinese))

[8] 郭華,金浩波,駱祖江.基于地面控制的江陰市淺層地下水資源評價[J].地質學刊,2009,33(1):64-69.(GUO Hua,JIN Haobo,LUO Zujiang.Shallow groundwater resource evaluation based on ground control in Jiangyin [J].Geological Journals,2009,33(1):64-69.(in Chinese))

[9] WEN H C.3D-groundwater modeling with PMWIN a simulation system for modeling groundwater flow and transport processes [M].New York:Springer,2001:120-132.

[10] 曾峰.第四系松散沉積層地下水開采與地面沉降耦合數值模擬[D].南京:河海大學,2008.

[11] 駱祖江,張英英,施春華.南通市第四紀松散沉積層地下水資源評價規劃三維數值模型[J].水資源保護,2009,25(5):19-23.(LUO Zujiang,ZHANG Yingying,SHI Chunhua.Three-dimensional numerical model for groundwater resource evaluation and plan in loose quaternary sediments in Nantong [J].Water Resources Protection,2009,25(5):19-23.(in Chinese))

[12] 駱祖江,姚炳魁,王曉梅.復合含水層系統地下水資源評價三維數值模型[J].吉林大學學報:地球科學版,2005,35(2):188-194.(LUO Zujiang,YAO Bingkui,WANG Xiaomei.Three dimension numerical model of groundwater resource evaluation of complex aquifer system [J].Journal of Jilin University:Earth Science Edition,2005,35(2):188-194.(in Chinese))

水利部批準《河湖生態修復與保護規劃編制導則》為水利行業標準

中華人民共和國水利部2015年6月2號發布2015年第44號水利部公告，批準SL709—2015《河湖生態修復與保護規劃編制導則》為水利行業標準，該標準將于2015年9月2號實施。《河湖生態修復與保護規劃編制導則》共分12 章36 節162 條和3 個附錄，主要包括以下內容：①標準的編制目的和適用范圍，規劃編制原則和要求、水平年設定、主要內容及執行的規范和標準等；②河湖生態現狀調查與評價的內容與要求；③規劃目標、指標、分區、總體布局的內容和要求；④生態需水、水質維護與改善、河湖地貌形態保護與修復、重要生物棲息地與生物多樣性保護、重要區域生態保護與修復等規劃措施的技術內容和要求；⑤河湖生態監控與管理、投資估算、規劃實施意見與綜合效益分析的內容和要求。

(本刊編輯部供稿)

Numerical simulation of groundwater exploitation and land subsidence in Cangzhou City

MA Qingshan,LUO Zujiang

(SchoolofEarthScienceandEngineering,HohaiUniversity,Nanjing210098,China)

Taking the land subsidence caused by groundwater exploitation in Cangzhou City as an example, based on the analysis of geology and hydrogeological conditions of the study area, a geological conceptual model of groundwater exploitation and land subsidence was generalized, setting up a corresponding coupled numerical model of three dimensional groundwater flow and one dimensional soil deformation.On the basis of the completion of identification and verification of model, we forecasted the change characteristics of groundwater flow field and the change trend of land subsidence from the end of 2025 to the end of 2010 under the current conditions of groundwater exploitation.The results show that the maximum cumulative land subsidence within the scope of Cangzhou City will be 475 mm at the end of 2025, the minimum value being 270 mm.

groundwater exploitation; land subsidence; prediction model; coupled numerical model

10.3880/j.issn.1004-6933.2015.04.004

河北省國土資源廳科技項目(CZCG2008008)

馬青山(1988—),男,博士研究生,研究方向為水文地質。E-mail:ssqm2007@163.com

P641.6

A

1004-6933(2015)04-0020-07

2014-11-12 編輯：高渭文)