沧州市地下水开采-地面沉降数值模拟

马青山,骆祖江

(河海大学地球科学与工程学院,江苏 南京 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 编辑：高渭文)