地下水开采诱发地面沉降研究及其工程应用：以山东德州地区为例
贾超++张国荣++王嘉斌++张永伟
摘要：地下水开采诱发的不良环境及地质问题日益凸显，其中地面沉降致灾问题已广受关注。华北平原山东区是若干地面沉降显著区之一，地下水开采是其主要诱发因素。针对地面沉降问题，以华北平原德州区为典型工程背景，采用水文地质工程地质与岩土力学相结合的方法，通过现场及室内试验，在获得研究区岩土力学参数的基础上，建立流固耦合数值模型，对研究区地面沉降的历史进行重演，对地面沉降的发展趋势进行预测。在获得研究区流场及地面沉降基本规律的基础上，提出减少地下水开采量、定期回灌地下水、调整开采井开采层位及开采井平面布局等工程措施进行地面沉降控制，并采用建立的数值模型从定量分析角度评价这些措施的效果。
关键词：地下水；地面沉降；流固耦合；水文地质；工程地质；岩土力学；华北平原
中图分类号：P641；P642.26文献标志码：A
0引言
地面沉降是一种自然因素或人类工程活动引发的地下松散地层固结压缩，从而导致一定区域范围内地面高程降低的缓变地质灾害现象。在众多诱发因素中，过量抽取地下水是造成地面沉降最主要的原因之一[1]。根据相关资料统计报道，中国目前已有超过50个城市发生了明显的地面沉降现象，其中最典型区域有长江三角洲、华北平原和汾渭盆地[2]。针对这些地区的地面沉降灾害问题，已有薛禹群、吴吉春、骆祖江、彭建兵、殷跃平、邵景力等众多专家学者进行了长期持续的工作，并取得了显著的研究成果[118]。相比于其他地区，华北平原山东地区地面沉降的机理研究工作尚显不足。在华北平原，山东地面沉降尤以鲁北地区最为严重，德州是鲁北地区地面沉降的典型代表。德州地区地面沉降已对当地环境及生产造成严重影响，制约当地社会经济的可持续发展[3]，并且根据该区地面沉降长期监测资料显示，地面沉降尚无明显停止的迹象[48]。
已有研究成果大多针对德州地区地下水开采与地面沉降的关系，应用统计学方法进行数理统计回归分析，但仅有这些成果并不能准确认识到该区地面沉降的深层机理及沉降本质。地下水开采诱发地面沉降的本质原因是：地下水开采导致各地层应力调整，引发各地层压缩总和在地面上显现。因此，采用水文地质工程地质与岩土力学相结合的方法更能揭示德州地区地面沉降的本质规律。该方法已被广泛应用于其他地区地面沉降研究中，但在德州地区尚无这方面的研究成果。
基于此，本文通过现场工作获取进尺达800 m的全取芯钻孔试样，进行试样的实验室测试，取得翔实的物理力学参数；分析流场分布规律，建立地面沉降流固耦合分析模型，重演历史沉降情景，预测沉降发展趋势；针对沉降发展趋势及规律，提出若干控制地面沉降的工程措施，并对提出的措施进行综合对比分析[1920] ，为山东德州地区防治地面沉降理论及分析计算奠定基础。
1研究区概况
山东省德州市德城区沉降资料及前期工作相对较好。该区东与德州市陵城区接壤，南与山东省平原县、武城县毗邻，西、西北、东北分别与河北省故城县、景县、吴桥县连接。其经度范围为116°10′15″E～116°30′40″E，纬度范围为 37°10′15″N～37°30′30″N，面积为539 km2。
研究区地下水多年处于超量开采状态，地下水位呈持续下降趋势，并伴随严重的地面沉降问题。该区自1989年开始进行地面精密测量，2006年400 mm沉降等值线圈定沉降漏斗面积为344 km2，历年沉降发展过程见图1。
图1山东省德州市德城区地面沉降发展过程
Fig.1Land Subsidence Development Process in Decheng District of Dezhou City， Shandong Province
2现场及室内试验
山东省德州市德城区在工程地质深层钻孔资料方面比较欠缺。山东省国土资源部门在该区进行了2处800 m深层全取芯钻孔钻探工作，山东大学承担了京台高速公路德州服务区钻位所得的岩芯参数测试工作。该处钻孔共取173个原状土样、15个平行样和37个扰动砂样（图2）。
根据现场地质钻探和室内试验结果，对德城区地层以工程地质及水文地质相近的原则进行地层概化。概化后德城区工程地质特征层分布见图3；实验室测试所得各地层物理力学参数取值范围见表1。
3数值模型
3.1数值模型的建立
以图3划分的工程地质特征层为所考虑的地层分布，以德州市德城区行政边界为模型外边界，建立研究区数值模型，总面积为539 km2。
数值模型共纳入开采深井263口，按各井实测流量赋值，部分井点信息见表2。井位分布见图4。
德州地面沉降观测井2号孔396～400 m深度处
图2现场岩芯钻孔地质编录
Fig.2Geological Record of Field Core Hole
开采井直径相比于整个区域面积非常微小，不利于数值模型的建立，因此，借鉴Oh等的研究成果[15]，将开采井概化为线，开采井边界条件采用质量流量进行表示。该方法既不会影响计算结果的精确度，又可使大量开采井实际位置及实际开采量
图3工程地质特征层划分
Fig.3Division of Engineering Geological Feature Layers
加入大区域模型得以顺利实现；经过验证，该方法正确有效。
表1地层物理力学参数
Tab.1Physical Mechanics Parameters of Strata
土层厚度/m密度/（kg·m-3）弹性模量/MPa泊松比渗透系数/（m·d-1）贮水系数/m-1孔隙度
第1层2001 900～2 190100～1500.30～0.451.800 0～2.400 00.010 00.40～0.50
第2层1002 000～2 20040～600.25～0.040.005 0～0.008 00.005 00.24～0.35
第3层1502 000～2 300110～1600.25～0.301.500 0～1.900 00.008 00.30～0.40
第4层1002 050～2 25050～650.20～0.330.004 0～0.005 50.000 10.24～0.27
第5层2502 100～2 300190～2100.17～0.221.400 0～1.600 00.007 00.30～0.40
表2开采井开采流量
Tab.2Flows of Groundwater Withdrawal from Extraction Wells
开采井编号开采井经纬度实测流量/（m3·d-1）质量流量/（kg·（m·s）-1）
K1（37°28′10.86″N，116°17′55.53″E）-82.19-0.006 341 821
K2（37°27′3.56″N，116°18′20.69″E）-82.19-0.006 341 821
K3（37°26′43.02″N，116°21′37.20″E）-109.59-0.008 456 019
K4（37°26′43.02″N，116°21′37.21″E）-136.99-0.010 570 216
K5（37°27′58.89″N，116°19′12.09″E）-82.19-0.006 341 821
K7（37°28′8.18″N，116°18′34.75″E）-493.15-0.038 051 698
K8（37°28′3.26″N，116°18′44.72″E）-301.37-0.023 253 858
K9（37°29′39.51″N，116°19′1.54″E）-104.11-0.008 033 179
注：开采井K3已出边界。
图4开采井分布
Fig.4Distribution of Extraction Wells
建立的数值模型见图5。三维网格共划分为202 964个单元，计算平台采用Comsol Multiphisices多物理场仿真分析平台。
图5数值模型三维网格剖分
Fig.53D Mesh Subdivision of Numerical Model
3.2流固耦合的实现及其支配方程
德州市德城区中深层承压水的大量抽取超过承压水的补给能力，引起承压水位降低，也就是土层的孔隙水压力降低导致了承压含水层土体的变形。这是一个包含渗流与固体力学相互作用的过程，其压缩变形具有一定的时间效应。
流固耦合支配方程主要由流体控制方程和多孔介质控制方程组成。本文采用比奥固结理论，其基本控制方程有流体控制方程和多孔介质控制方程。
（1）流体控制方程为
SymbolQC@ σ=ρg（2）
式中：σ为总应力张量；ρ为总密度；g为重力加速度。
流固耦合的实现过程就是在初始条件和边界条件的基础上对流体控制方程和多孔介质控制方程的耦合求解过程。
3.3初始条件和边界条件
德州市德城区地面沉降监测资料自1991年开始，取为计算初始时间点。整个模型外边界有水头监控点，采用水头边界；模型内所含263口开采井取为流量边界，按实际开采量施加，考虑数值模型建立及计算简便，均转化为质量流量（图6）。
3.4模型识别与校验
采用实测岩土力学参数、初始条件及边界条件对德州市德城区地面沉降进行分析计算。图7为1991～2007年德州市德城区计算地面沉降量分布。将图7与图1的实测地面沉降进行对比，计算所得整体沉降趋势及沉降量与实测资料较为吻合。
图6三维开采井点位置
Fig.6Location of 3D Extraction Well Points
图71999～2007年累计地面沉降（单位：m）
Fig.7Accumulative Land Subsidence from 1999 to 2007 （Unit：m）
为进一步说明建立的数值模型合理性及正确性，从图7所得流固耦合三维数值模拟计算结果中提取德州市德城区8个点位（平均覆盖整个计算区域），与该8个点位的实测沉降量进行对比分析。
提取的8个点位逐年沉降模拟值及实测值对比曲线见图8。从图8可知，沉降量模拟值与实测值基本吻合，满足精度需求。建立的数值模型可以较好地重现德州市德城区地下水开采诱发地面沉降的历史过程。
图8沉降监测点实测值与模拟值对比曲线
Fig.8Comparison Curves Between Measured and Simulated Land Subsidences in the Monitoring Points
4地面沉降发展趋势
应用建立的数值模型预测德州市德城区地面沉降在现状地下水开采量情况下的地面沉降发展趋势。
以2010年作为地面沉降预测初始时刻，地下水开采量采用同期各开采井统计量。图9展示了计算得到的德州市德城区1 m沉降范围演化。由图9可知，研究区若保持现状的地下水开采条件，其沉降范围具有不断扩展的趋势。沉降中心位置基本保持稳定，但对同一位置沉降量随时间的推移不断加大，由此形成的沉降漏斗范围不断扩大，漏斗中心高程不断下降，因此，沉降灾害具有加剧趋势，应采取切实有效的控制措施加以处置。
图9不同年份1 m沉降范围演化预测（单位：m）
Fig.9Evolution Prediction of 1 m Subsidence Range in Different Years （Unit：m）
表3列出了德州市德城区代表性点位的沉降量预测值。根据数值模型预测结果，在现状条件下开
表3现状地下水开采条件下不同沉降监测点地面沉降预测结果
Tab.3Predictions of Land Subsidence in Different Monitoring Points Under the Condition of Current Groundwater Exploitation
年份监测点D62地面沉降量/mm监测点Q1地面沉降量/mm监测点Q6地面沉降量/mm监测点D101地面沉降量/mm监测点D14地面沉降量/mm监测点DJ9地面沉降量/mm监测点D110地面沉降量/mm监测点D111地面沉降量/mm
20151 316.21 269.5910.8962.8503.3801.6687.5910.8
20161 345.01 297.3931.8986.5514.3823.7705.5931.8
20171 373.71 324.9952.31 009.2525.7844.6723.0952.3
20181 402.51 352.5972.91 031.7537.2865.4740.4972.9
20191 431.41 380.1993.51 053.7548.8885.6757.7993.5
20201 460.61 407.91 014.21 075.6560.5905.7775.01 014.2
20211 490.41 436.11 035.01 097.2572.3925.6792.51 035.0
20221 520.21 464.21 055.71 118.9584.2945.5809.91 055.7
20231 550.01 492.41 076.51 140.5596.0965.3827.31 076.5
20241 580.21 521.01 097.41 162.3608.0985.3845.01 097.4
20251 611.41 550.31 118.51 184.3620.21 005.5863.11 118.5
采地下水10年，总体地面沉降速率相对于2007年之前有所减缓，地面沉降中心位置基本稳定，但累积沉降量不断加大。尤其是在开采井分布较为集中的地段，累计沉降量增大的趋势尤为严重。从开采井分布位置分析，沉降漏斗基本处于开采井集中分布区，因此，集中开采地下水将会加剧地面沉降量，导致漏斗高程加速下降，若不及时采取有效措施，地下水开采诱发的地面沉降将越发难以控制。
5地面沉降控制措施
德州市德城区非常迫切地需要进行地面沉降控制，本文提出减小地下水开采量、定期回灌地下水、调整开采井开采层位及开采井平面布局等工程措施，并采用建立的数值模型从定量分析角度评价这些措施的效果。
5.1减小地下水开采量
根据《山东省地面沉降防治规划（2012～2020年）》，从2015年起减少德州市德城区第一和第二承压含水层的地下水开采量。根据建立的数值模型预测未来20年地面沉降的变化情况，计算方案见表4。根据数值模拟结果，随着地下水开采量的减少，德州市德城区地面沉降速率明显降低（图10）。当地下水开采量减少30%及以上时，部分地面沉降监测点在减少地下水开采量的最初几年内出现了地面沉降趋于停止或少量回弹的现象，但是随着时间推移，回弹量减少，地面沉降继续发展。总体来说，地面沉降速率较未减少地下水开采量前明显降低，因此，减少地下水开采量在一定程度上可以有效控制德州市德城区地面沉降。
表4减少地下水开采量计算方案
Tab.4Schemes for Reducing the Groundwater Exploitation
方案1方案2方案3方案4方案5
开采量减少10%开采量减少20%开采量减少30%开采量减少40%开采量减少50%
图10减少地下水开采量不同计算方案下沉降监测点累计地面沉降量对比
Fig.10Comparisons of Cumulative Land Subsidence Under Different Schemes for Reducing the Groundwater Exploitation
5.2定期回灌地下水
人工回灌是指利用工程设施将地表水注入地下含水层，以增加地下水储量的措施。根据德州市德城区地面沉降中心分布状况，在其中心布设15口回灌井，这些井年回灌量相同，根据总回灌量的不同，设置4种回灌方案（表5）。假设回灌的地下水均达到承压含水层，根据计算结果，当地下水回灌量较少时，对地面沉降速率的影响较小；回灌井附近地面沉降量受地下水回灌影响较大，远离回灌井的地面效果不明显（图11）。
图11不同回灌方案下沉降监测点累计地面沉降量对比
Fig.11Comparisons of Cumulative Land Subsidence Under Various Recharge Schemes
表5地下水回灌方案
Tab.5Schemes for Groundwater Recharge
方案方案1方案2方案3方案4
总回灌量/（104m3·年-1）200400600800
5.3调整开采井开采层位
根据《德州市地面沉降地质调查与评价》，德州市德城区浅层微咸水、咸水盈余量每年为395.3×104 m3，尚处于未开发阶段。由此设计了调整开采井开采层位的方案：从2015年起，潜水层每年多开采390×104 m3的水量，将第一或第二承压含水层的地下水开采量每年减少390×104 m3。由于潜水层水质较差，可将开采的潜水层地下水作为城市绿化、消防和对水质要求不高的工业用水水源。
通过计算，调整地下水开采层位后的累计沉降量与未调整开采层位时相比有大幅度降低；同时，减少第二承压含水层的地下水开采量比减少等量的第一承压含水层的开采量对地面沉降的控制效果要好，而且沉降中心有所减缓。因此，在条件允许的情况下，地下水尽量从潜水层或上部第一承压含水层开采。
5.4调整开采井平面布局
在开采量一定的情况下，为缓解地面沉降梯度大的问题，可通过调整开采井及开采量的分布来实现对地面沉降的控制。结合德州市的实际情况，重点针对地面沉降严重的沉降中心区进行开采井平面布局研究。调整开采井平面布局后，沉降中心区的开采井间距至少达到1 000 m（图12）。假设沉降中心区开采井的开采流量重新调整，即各开采井的年开采量相同。
图12开采井平面布局的调整
Fig.12Adjustment of Plane Layout for Extraction Wells
根据布局调整进行计算，结果显示在总开采量不变的条件下，重新调整开采井和开采流量的布局，可以在很大程度上减小或减缓德州市德城区地面沉降速率，使沉降中心沉降梯度减小。因此，在开采量一定的情况下，为最大程度减少地面沉降，可以根据德州市德城区的实际情况，适当调整开采井的分布和各开采井的年开采量。
5.5小结
根据各工程措施所进行的计算分析结果表明，各工程措施均对控制研究区地面沉降产生效果，但也存在差异性。因此，可以考虑将多种方案进行综合利用：根据当地社会经济发展趋势，合理规划城市发展布局，这可以与地下水开采井平面布局相结合；在不同阶段可以设置不同的地下水开采量上限，在枯水季和丰水期可以设置不同的量值；在科学合理评估各开采层位不同地下水开采量情况下，可以随时间的推移适当调整不同开采层位，尽量保持各层位具有回灌回弹的特性，但在回灌过程中一定要注意水环境污染控制问题，因为地下水一旦遭到污染，治理将会非常困难。
6结语
地下水开采已导致华北平原德州区出现了大范围的地面沉降，且有加剧的趋势。已有研究成果已不能从根本上解决地面沉降控制的深层次问题。本文从地面沉降的本质出发，采用水文地质工程地质及岩土力学相结合的方法，建立山东省德州市德城区地面沉降流固耦合数值模型，成功反演该区地面沉降历史发展情景，并对该区未来地面沉降发展趋势进行了定量化预测；在此基础上，提出研究区地面沉降若干控制措施，并对提出的工程措施进行评价。本次研究成果是山东省德州市地面沉降研究总结，为该区在地面沉降机理方面的进一步研究奠定了基础。
参考文献：
References：
[1]薛禹群，吴吉春，张云，等.长江三角洲（南部）区域地面沉降模拟研究[J].中国科学：D辑，地球科学，2008，38（4）：477492.
XUE Yuqun，WU Jichun，ZHANG Yun，et al.Simulation of Land Subsidence in the Southern Yangtze River Delta Plain[J].Science in China：Series D，Earth Sciences，2008，38（4）：477492.
[2]殷跃平，张作辰，张开军.我国地面沉降现状及防治对策研究[J].中国地质灾害与防治学报，2005，16（2）：18.
YIN Yueping，ZHANG Zuochen，ZHANG Kaijun.Land Subsidence and Countermeasures for Its Prevention in China[J].The Chinese Journal of Geological Hazard and Control，2005，16（2）：18.
[3]何国荣，冯在敏，冯克印.水准测量网监测在德州市地面沉降中的应用[J].山东国土资源，2012，28（6）：2830.
HE Guorong，FENG Zaimin，FENG Keyin，et al.The Monitoring of Land Subsidence Based on the Leveling Surveying Network in Dezhou City[J].Land and Resources in Shandong Province，2012，28（6）：2830.
[4]王贞国，王彦俊，冯守涛，等.德州市深层地下水人工回灌试验浅探[J].山东国土资源，2005，21（6/7）：8285.
WANG Zhenguo，WANG Yanjun，FENG Shoutao，et al.Study on Artificial Recharge Experiment of Deep Underground Water in Dezhou City[J].Land and Resources in Shandong Province，2005，21（6/7）：8285.
[5]孟庆峰，戴鲁旗.德州市地面沉降现状及防治对策[J].山东国土资源，2008，24（11）：89.
MENG Qingfeng，DAI Luqi.Land Subsidence and Countermeasures for Its Prevention in Dezhou City[J].Land and Resources in Shandong Province，2008，24（11）：89.
[6]王彦俊，刘桂仪.德州深层地下水降落漏斗相关模型建立及其应用[J].山东地质，1998，14（2）：3845.
WANG Yanjun，LIU Guiyi.Correlation Models Establishment and Its Application in Dezhou Funnel of Internal Underground Water[J].Geology of Shandong，1998，14（2）：3845.
[7]杨丽芝.德州深层地下水位降落漏斗演变机制与可调控性研究[D].北京：中国地质科学院，2009.
YANG Lizhi.From Princple and Controllingadjusting Research About Deep Groundwater Depression Cone in Dezhou[D].Beijing：Chinese Academy of Geological Sciences，2009.
[8]蔡文晓.德州市深层地下水开采与地面沉降关系研究[D].长春：吉林大学，2009.
CAI Wenxiao.Research on Relationship Between Deep Groundwater Exploitation and Land Subsidence in Dezhou[D].Changchun：Jilin University，2009.
[9]姬亚东，柴学周，刘其声，等.大区域地下水流数值模拟研究现状及存在问题[J].煤田地质与勘探，2009，37（5）：3236.
JI Yadong，CHAI Xuezhou，LIU Qisheng，et al.Current Situation and Existing Problems of Numerical Simulation of Groundwater Flow in a Largescale Area[J].Coal Geology and Exploration，2009，37（5）：3236.
[10] ZHAO C Y，ZHANG Q，DING X L，et al.Monitoring of Land Subsidence and Ground Fissures in Xian，China 20052006：Mapped by SAR Interferometry[J].Environmental Geology，2009，58（7）：15331540.
[11]骆祖江，曾峰，李颖.地下水开采与地面沉降控制三维全耦合模型研究[J].吉林大学学报：地球科学版，2009，39（6）：10801086.
LUO Zujiang，ZENG Feng，LI Ying.Study on Three Dimensional Full Coupling Model of Groundwater Exploitation and Landsubsidence Control[J].Journal of Jilin University：Earth Science Edition，2009，39（6）：10801086.
[12]刘聪，彭建兵，陈立伟，等.西安地裂缝场地地震效应分析[J].工程地质学报，2012，20（6）：979985.
LIU Cong，PENG Jianbing，CHEN Liwei，et al.Finite Element Analysis of Seismic Effects on Ground Fissures in Xian[J].Journal of Engineering Geology，2012，20（6）：979985.
[13]魏加华，崔亚莉，邵景力，等.济宁市地下水与地面沉降三维有限元模拟[J].长春科技大学学报，2000，30（4）：376380.
WEI Jiahua，CUI Yali，SHAO Jingli，et al.3Dfinite Element Method Simulation of Groundwater and Land Subsidence in Jining City[J].Journal of Changchun University of Science and Technology，2000，30（4）：376380.
[14]万伟峰.西安市地下水开采地面沉降数值模拟及防治方案研究[D].西安：长安大学，2008.
WAN Weifeng.Numerical Simulation of Land Subsidence Due to Groundwater Withdraw in Xian and Its Prevention Schemes[D].Xian：Changan University，2008.
[15]OH H J，AHN S C，CHOI J K，et al.Sensitivity Analysis for the GISbased Mapping of the Ground Subsidence Hazard near Abandoned Underground Coal Mines[J].Environmental Earth Sciences，2011，64（2）：347358.
[16]SUN Y H，ZHANG X D，MAO W F，et al.Mechanism and Stability Evaluation of Goaf Ground Subsidence in the Third Mining Area in Gong Changling District，China[J].Arabian Journal of Geosciences，2015，8（2）：639646.
[17]KIM K D，LEE S，OH H J.Prediction of Ground Subsidence in Samcheok City，Korea Using Artificial Neural Networks and GIS[J].Environmental Geology，2009，58（1）：7172.
[18]OUYANG Z H.The Role of Potential Soil Cavity on Ground Subsidence and Collapse in Coal Mining Area[J].Journal of Coal Science and Engineering，2010，16（3）：240245.
[19]张永伟.华北平原德州地面沉降成生机理、监测预警与可控性研究[D].济南：山东大学，2014.
ZHANG Yongwei.Formation Mechanism，Monitoring and Warning，Controlling Research of Subsidence of Dezhou in North China Plain[D].Jinan：Shandong University，2014.
[20]贾超，张国荣，何奔，等.华北平原德州地面沉降监测预警系统建设报告[R].济南：山东大学，2014.
JIA Chao，ZHANG Guorong，HE Ben，et al.Land Subsidence Monitoring and Warning System Construction Report in Dezhou District of North China Plain[R].Jinan：Shandong University，2014.
摘要：地下水开采诱发的不良环境及地质问题日益凸显，其中地面沉降致灾问题已广受关注。华北平原山东区是若干地面沉降显著区之一，地下水开采是其主要诱发因素。针对地面沉降问题，以华北平原德州区为典型工程背景，采用水文地质工程地质与岩土力学相结合的方法，通过现场及室内试验，在获得研究区岩土力学参数的基础上，建立流固耦合数值模型，对研究区地面沉降的历史进行重演，对地面沉降的发展趋势进行预测。在获得研究区流场及地面沉降基本规律的基础上，提出减少地下水开采量、定期回灌地下水、调整开采井开采层位及开采井平面布局等工程措施进行地面沉降控制，并采用建立的数值模型从定量分析角度评价这些措施的效果。
关键词：地下水；地面沉降；流固耦合；水文地质；工程地质；岩土力学；华北平原
中图分类号：P641；P642.26文献标志码：A
0引言
地面沉降是一种自然因素或人类工程活动引发的地下松散地层固结压缩，从而导致一定区域范围内地面高程降低的缓变地质灾害现象。在众多诱发因素中，过量抽取地下水是造成地面沉降最主要的原因之一[1]。根据相关资料统计报道，中国目前已有超过50个城市发生了明显的地面沉降现象，其中最典型区域有长江三角洲、华北平原和汾渭盆地[2]。针对这些地区的地面沉降灾害问题，已有薛禹群、吴吉春、骆祖江、彭建兵、殷跃平、邵景力等众多专家学者进行了长期持续的工作，并取得了显著的研究成果[118]。相比于其他地区，华北平原山东地区地面沉降的机理研究工作尚显不足。在华北平原，山东地面沉降尤以鲁北地区最为严重，德州是鲁北地区地面沉降的典型代表。德州地区地面沉降已对当地环境及生产造成严重影响，制约当地社会经济的可持续发展[3]，并且根据该区地面沉降长期监测资料显示，地面沉降尚无明显停止的迹象[48]。
已有研究成果大多针对德州地区地下水开采与地面沉降的关系，应用统计学方法进行数理统计回归分析，但仅有这些成果并不能准确认识到该区地面沉降的深层机理及沉降本质。地下水开采诱发地面沉降的本质原因是：地下水开采导致各地层应力调整，引发各地层压缩总和在地面上显现。因此，采用水文地质工程地质与岩土力学相结合的方法更能揭示德州地区地面沉降的本质规律。该方法已被广泛应用于其他地区地面沉降研究中，但在德州地区尚无这方面的研究成果。
基于此，本文通过现场工作获取进尺达800 m的全取芯钻孔试样，进行试样的实验室测试，取得翔实的物理力学参数；分析流场分布规律，建立地面沉降流固耦合分析模型，重演历史沉降情景，预测沉降发展趋势；针对沉降发展趋势及规律，提出若干控制地面沉降的工程措施，并对提出的措施进行综合对比分析[1920] ，为山东德州地区防治地面沉降理论及分析计算奠定基础。
1研究区概况
山东省德州市德城区沉降资料及前期工作相对较好。该区东与德州市陵城区接壤，南与山东省平原县、武城县毗邻，西、西北、东北分别与河北省故城县、景县、吴桥县连接。其经度范围为116°10′15″E～116°30′40″E，纬度范围为 37°10′15″N～37°30′30″N，面积为539 km2。
研究区地下水多年处于超量开采状态，地下水位呈持续下降趋势，并伴随严重的地面沉降问题。该区自1989年开始进行地面精密测量，2006年400 mm沉降等值线圈定沉降漏斗面积为344 km2，历年沉降发展过程见图1。
图1山东省德州市德城区地面沉降发展过程
Fig.1Land Subsidence Development Process in Decheng District of Dezhou City， Shandong Province
2现场及室内试验
山东省德州市德城区在工程地质深层钻孔资料方面比较欠缺。山东省国土资源部门在该区进行了2处800 m深层全取芯钻孔钻探工作，山东大学承担了京台高速公路德州服务区钻位所得的岩芯参数测试工作。该处钻孔共取173个原状土样、15个平行样和37个扰动砂样（图2）。
根据现场地质钻探和室内试验结果，对德城区地层以工程地质及水文地质相近的原则进行地层概化。概化后德城区工程地质特征层分布见图3；实验室测试所得各地层物理力学参数取值范围见表1。
3数值模型
3.1数值模型的建立
以图3划分的工程地质特征层为所考虑的地层分布，以德州市德城区行政边界为模型外边界，建立研究区数值模型，总面积为539 km2。
数值模型共纳入开采深井263口，按各井实测流量赋值，部分井点信息见表2。井位分布见图4。
德州地面沉降观测井2号孔396～400 m深度处
图2现场岩芯钻孔地质编录
Fig.2Geological Record of Field Core Hole
开采井直径相比于整个区域面积非常微小，不利于数值模型的建立，因此，借鉴Oh等的研究成果[15]，将开采井概化为线，开采井边界条件采用质量流量进行表示。该方法既不会影响计算结果的精确度，又可使大量开采井实际位置及实际开采量
图3工程地质特征层划分
Fig.3Division of Engineering Geological Feature Layers
加入大区域模型得以顺利实现；经过验证，该方法正确有效。
表1地层物理力学参数
Tab.1Physical Mechanics Parameters of Strata
土层厚度/m密度/（kg·m-3）弹性模量/MPa泊松比渗透系数/（m·d-1）贮水系数/m-1孔隙度
第1层2001 900～2 190100～1500.30～0.451.800 0～2.400 00.010 00.40～0.50
第2层1002 000～2 20040～600.25～0.040.005 0～0.008 00.005 00.24～0.35
第3层1502 000～2 300110～1600.25～0.301.500 0～1.900 00.008 00.30～0.40
第4层1002 050～2 25050～650.20～0.330.004 0～0.005 50.000 10.24～0.27
第5层2502 100～2 300190～2100.17～0.221.400 0～1.600 00.007 00.30～0.40
表2开采井开采流量
Tab.2Flows of Groundwater Withdrawal from Extraction Wells
开采井编号开采井经纬度实测流量/（m3·d-1）质量流量/（kg·（m·s）-1）
K1（37°28′10.86″N，116°17′55.53″E）-82.19-0.006 341 821
K2（37°27′3.56″N，116°18′20.69″E）-82.19-0.006 341 821
K3（37°26′43.02″N，116°21′37.20″E）-109.59-0.008 456 019
K4（37°26′43.02″N，116°21′37.21″E）-136.99-0.010 570 216
K5（37°27′58.89″N，116°19′12.09″E）-82.19-0.006 341 821
K7（37°28′8.18″N，116°18′34.75″E）-493.15-0.038 051 698
K8（37°28′3.26″N，116°18′44.72″E）-301.37-0.023 253 858
K9（37°29′39.51″N，116°19′1.54″E）-104.11-0.008 033 179
注：开采井K3已出边界。
图4开采井分布
Fig.4Distribution of Extraction Wells
建立的数值模型见图5。三维网格共划分为202 964个单元，计算平台采用Comsol Multiphisices多物理场仿真分析平台。
图5数值模型三维网格剖分
Fig.53D Mesh Subdivision of Numerical Model
3.2流固耦合的实现及其支配方程
德州市德城区中深层承压水的大量抽取超过承压水的补给能力，引起承压水位降低，也就是土层的孔隙水压力降低导致了承压含水层土体的变形。这是一个包含渗流与固体力学相互作用的过程，其压缩变形具有一定的时间效应。
流固耦合支配方程主要由流体控制方程和多孔介质控制方程组成。本文采用比奥固结理论，其基本控制方程有流体控制方程和多孔介质控制方程。
（1）流体控制方程为
SymbolQC@ σ=ρg（2）
式中：σ为总应力张量；ρ为总密度；g为重力加速度。
流固耦合的实现过程就是在初始条件和边界条件的基础上对流体控制方程和多孔介质控制方程的耦合求解过程。
3.3初始条件和边界条件
德州市德城区地面沉降监测资料自1991年开始，取为计算初始时间点。整个模型外边界有水头监控点，采用水头边界；模型内所含263口开采井取为流量边界，按实际开采量施加，考虑数值模型建立及计算简便，均转化为质量流量（图6）。
3.4模型识别与校验
采用实测岩土力学参数、初始条件及边界条件对德州市德城区地面沉降进行分析计算。图7为1991～2007年德州市德城区计算地面沉降量分布。将图7与图1的实测地面沉降进行对比，计算所得整体沉降趋势及沉降量与实测资料较为吻合。
图6三维开采井点位置
Fig.6Location of 3D Extraction Well Points
图71999～2007年累计地面沉降（单位：m）
Fig.7Accumulative Land Subsidence from 1999 to 2007 （Unit：m）
为进一步说明建立的数值模型合理性及正确性，从图7所得流固耦合三维数值模拟计算结果中提取德州市德城区8个点位（平均覆盖整个计算区域），与该8个点位的实测沉降量进行对比分析。
提取的8个点位逐年沉降模拟值及实测值对比曲线见图8。从图8可知，沉降量模拟值与实测值基本吻合，满足精度需求。建立的数值模型可以较好地重现德州市德城区地下水开采诱发地面沉降的历史过程。
图8沉降监测点实测值与模拟值对比曲线
Fig.8Comparison Curves Between Measured and Simulated Land Subsidences in the Monitoring Points
4地面沉降发展趋势
应用建立的数值模型预测德州市德城区地面沉降在现状地下水开采量情况下的地面沉降发展趋势。
以2010年作为地面沉降预测初始时刻，地下水开采量采用同期各开采井统计量。图9展示了计算得到的德州市德城区1 m沉降范围演化。由图9可知，研究区若保持现状的地下水开采条件，其沉降范围具有不断扩展的趋势。沉降中心位置基本保持稳定，但对同一位置沉降量随时间的推移不断加大，由此形成的沉降漏斗范围不断扩大，漏斗中心高程不断下降，因此，沉降灾害具有加剧趋势，应采取切实有效的控制措施加以处置。
图9不同年份1 m沉降范围演化预测（单位：m）
Fig.9Evolution Prediction of 1 m Subsidence Range in Different Years （Unit：m）
表3列出了德州市德城区代表性点位的沉降量预测值。根据数值模型预测结果，在现状条件下开
表3现状地下水开采条件下不同沉降监测点地面沉降预测结果
Tab.3Predictions of Land Subsidence in Different Monitoring Points Under the Condition of Current Groundwater Exploitation
年份监测点D62地面沉降量/mm监测点Q1地面沉降量/mm监测点Q6地面沉降量/mm监测点D101地面沉降量/mm监测点D14地面沉降量/mm监测点DJ9地面沉降量/mm监测点D110地面沉降量/mm监测点D111地面沉降量/mm
20151 316.21 269.5910.8962.8503.3801.6687.5910.8
20161 345.01 297.3931.8986.5514.3823.7705.5931.8
20171 373.71 324.9952.31 009.2525.7844.6723.0952.3
20181 402.51 352.5972.91 031.7537.2865.4740.4972.9
20191 431.41 380.1993.51 053.7548.8885.6757.7993.5
20201 460.61 407.91 014.21 075.6560.5905.7775.01 014.2
20211 490.41 436.11 035.01 097.2572.3925.6792.51 035.0
20221 520.21 464.21 055.71 118.9584.2945.5809.91 055.7
20231 550.01 492.41 076.51 140.5596.0965.3827.31 076.5
20241 580.21 521.01 097.41 162.3608.0985.3845.01 097.4
20251 611.41 550.31 118.51 184.3620.21 005.5863.11 118.5
采地下水10年，总体地面沉降速率相对于2007年之前有所减缓，地面沉降中心位置基本稳定，但累积沉降量不断加大。尤其是在开采井分布较为集中的地段，累计沉降量增大的趋势尤为严重。从开采井分布位置分析，沉降漏斗基本处于开采井集中分布区，因此，集中开采地下水将会加剧地面沉降量，导致漏斗高程加速下降，若不及时采取有效措施，地下水开采诱发的地面沉降将越发难以控制。
5地面沉降控制措施
德州市德城区非常迫切地需要进行地面沉降控制，本文提出减小地下水开采量、定期回灌地下水、调整开采井开采层位及开采井平面布局等工程措施，并采用建立的数值模型从定量分析角度评价这些措施的效果。
5.1减小地下水开采量
根据《山东省地面沉降防治规划（2012～2020年）》，从2015年起减少德州市德城区第一和第二承压含水层的地下水开采量。根据建立的数值模型预测未来20年地面沉降的变化情况，计算方案见表4。根据数值模拟结果，随着地下水开采量的减少，德州市德城区地面沉降速率明显降低（图10）。当地下水开采量减少30%及以上时，部分地面沉降监测点在减少地下水开采量的最初几年内出现了地面沉降趋于停止或少量回弹的现象，但是随着时间推移，回弹量减少，地面沉降继续发展。总体来说，地面沉降速率较未减少地下水开采量前明显降低，因此，减少地下水开采量在一定程度上可以有效控制德州市德城区地面沉降。
表4减少地下水开采量计算方案
Tab.4Schemes for Reducing the Groundwater Exploitation
方案1方案2方案3方案4方案5
开采量减少10%开采量减少20%开采量减少30%开采量减少40%开采量减少50%
图10减少地下水开采量不同计算方案下沉降监测点累计地面沉降量对比
Fig.10Comparisons of Cumulative Land Subsidence Under Different Schemes for Reducing the Groundwater Exploitation
5.2定期回灌地下水
人工回灌是指利用工程设施将地表水注入地下含水层，以增加地下水储量的措施。根据德州市德城区地面沉降中心分布状况，在其中心布设15口回灌井，这些井年回灌量相同，根据总回灌量的不同，设置4种回灌方案（表5）。假设回灌的地下水均达到承压含水层，根据计算结果，当地下水回灌量较少时，对地面沉降速率的影响较小；回灌井附近地面沉降量受地下水回灌影响较大，远离回灌井的地面效果不明显（图11）。
图11不同回灌方案下沉降监测点累计地面沉降量对比
Fig.11Comparisons of Cumulative Land Subsidence Under Various Recharge Schemes
表5地下水回灌方案
Tab.5Schemes for Groundwater Recharge
方案方案1方案2方案3方案4
总回灌量/（104m3·年-1）200400600800
5.3调整开采井开采层位
根据《德州市地面沉降地质调查与评价》，德州市德城区浅层微咸水、咸水盈余量每年为395.3×104 m3，尚处于未开发阶段。由此设计了调整开采井开采层位的方案：从2015年起，潜水层每年多开采390×104 m3的水量，将第一或第二承压含水层的地下水开采量每年减少390×104 m3。由于潜水层水质较差，可将开采的潜水层地下水作为城市绿化、消防和对水质要求不高的工业用水水源。
通过计算，调整地下水开采层位后的累计沉降量与未调整开采层位时相比有大幅度降低；同时，减少第二承压含水层的地下水开采量比减少等量的第一承压含水层的开采量对地面沉降的控制效果要好，而且沉降中心有所减缓。因此，在条件允许的情况下，地下水尽量从潜水层或上部第一承压含水层开采。
5.4调整开采井平面布局
在开采量一定的情况下，为缓解地面沉降梯度大的问题，可通过调整开采井及开采量的分布来实现对地面沉降的控制。结合德州市的实际情况，重点针对地面沉降严重的沉降中心区进行开采井平面布局研究。调整开采井平面布局后，沉降中心区的开采井间距至少达到1 000 m（图12）。假设沉降中心区开采井的开采流量重新调整，即各开采井的年开采量相同。
图12开采井平面布局的调整
Fig.12Adjustment of Plane Layout for Extraction Wells
根据布局调整进行计算，结果显示在总开采量不变的条件下，重新调整开采井和开采流量的布局，可以在很大程度上减小或减缓德州市德城区地面沉降速率，使沉降中心沉降梯度减小。因此，在开采量一定的情况下，为最大程度减少地面沉降，可以根据德州市德城区的实际情况，适当调整开采井的分布和各开采井的年开采量。
5.5小结
根据各工程措施所进行的计算分析结果表明，各工程措施均对控制研究区地面沉降产生效果，但也存在差异性。因此，可以考虑将多种方案进行综合利用：根据当地社会经济发展趋势，合理规划城市发展布局，这可以与地下水开采井平面布局相结合；在不同阶段可以设置不同的地下水开采量上限，在枯水季和丰水期可以设置不同的量值；在科学合理评估各开采层位不同地下水开采量情况下，可以随时间的推移适当调整不同开采层位，尽量保持各层位具有回灌回弹的特性，但在回灌过程中一定要注意水环境污染控制问题，因为地下水一旦遭到污染，治理将会非常困难。
6结语
地下水开采已导致华北平原德州区出现了大范围的地面沉降，且有加剧的趋势。已有研究成果已不能从根本上解决地面沉降控制的深层次问题。本文从地面沉降的本质出发，采用水文地质工程地质及岩土力学相结合的方法，建立山东省德州市德城区地面沉降流固耦合数值模型，成功反演该区地面沉降历史发展情景，并对该区未来地面沉降发展趋势进行了定量化预测；在此基础上，提出研究区地面沉降若干控制措施，并对提出的工程措施进行评价。本次研究成果是山东省德州市地面沉降研究总结，为该区在地面沉降机理方面的进一步研究奠定了基础。
参考文献：
References：
[1]薛禹群，吴吉春，张云，等.长江三角洲（南部）区域地面沉降模拟研究[J].中国科学：D辑，地球科学，2008，38（4）：477492.
XUE Yuqun，WU Jichun，ZHANG Yun，et al.Simulation of Land Subsidence in the Southern Yangtze River Delta Plain[J].Science in China：Series D，Earth Sciences，2008，38（4）：477492.
[2]殷跃平，张作辰，张开军.我国地面沉降现状及防治对策研究[J].中国地质灾害与防治学报，2005，16（2）：18.
YIN Yueping，ZHANG Zuochen，ZHANG Kaijun.Land Subsidence and Countermeasures for Its Prevention in China[J].The Chinese Journal of Geological Hazard and Control，2005，16（2）：18.
[3]何国荣，冯在敏，冯克印.水准测量网监测在德州市地面沉降中的应用[J].山东国土资源，2012，28（6）：2830.
HE Guorong，FENG Zaimin，FENG Keyin，et al.The Monitoring of Land Subsidence Based on the Leveling Surveying Network in Dezhou City[J].Land and Resources in Shandong Province，2012，28（6）：2830.
[4]王贞国，王彦俊，冯守涛，等.德州市深层地下水人工回灌试验浅探[J].山东国土资源，2005，21（6/7）：8285.
WANG Zhenguo，WANG Yanjun，FENG Shoutao，et al.Study on Artificial Recharge Experiment of Deep Underground Water in Dezhou City[J].Land and Resources in Shandong Province，2005，21（6/7）：8285.
[5]孟庆峰，戴鲁旗.德州市地面沉降现状及防治对策[J].山东国土资源，2008，24（11）：89.
MENG Qingfeng，DAI Luqi.Land Subsidence and Countermeasures for Its Prevention in Dezhou City[J].Land and Resources in Shandong Province，2008，24（11）：89.
[6]王彦俊，刘桂仪.德州深层地下水降落漏斗相关模型建立及其应用[J].山东地质，1998，14（2）：3845.
WANG Yanjun，LIU Guiyi.Correlation Models Establishment and Its Application in Dezhou Funnel of Internal Underground Water[J].Geology of Shandong，1998，14（2）：3845.
[7]杨丽芝.德州深层地下水位降落漏斗演变机制与可调控性研究[D].北京：中国地质科学院，2009.
YANG Lizhi.From Princple and Controllingadjusting Research About Deep Groundwater Depression Cone in Dezhou[D].Beijing：Chinese Academy of Geological Sciences，2009.
[8]蔡文晓.德州市深层地下水开采与地面沉降关系研究[D].长春：吉林大学，2009.
CAI Wenxiao.Research on Relationship Between Deep Groundwater Exploitation and Land Subsidence in Dezhou[D].Changchun：Jilin University，2009.
[9]姬亚东，柴学周，刘其声，等.大区域地下水流数值模拟研究现状及存在问题[J].煤田地质与勘探，2009，37（5）：3236.
JI Yadong，CHAI Xuezhou，LIU Qisheng，et al.Current Situation and Existing Problems of Numerical Simulation of Groundwater Flow in a Largescale Area[J].Coal Geology and Exploration，2009，37（5）：3236.
[10] ZHAO C Y，ZHANG Q，DING X L，et al.Monitoring of Land Subsidence and Ground Fissures in Xian，China 20052006：Mapped by SAR Interferometry[J].Environmental Geology，2009，58（7）：15331540.
[11]骆祖江，曾峰，李颖.地下水开采与地面沉降控制三维全耦合模型研究[J].吉林大学学报：地球科学版，2009，39（6）：10801086.
LUO Zujiang，ZENG Feng，LI Ying.Study on Three Dimensional Full Coupling Model of Groundwater Exploitation and Landsubsidence Control[J].Journal of Jilin University：Earth Science Edition，2009，39（6）：10801086.
[12]刘聪，彭建兵，陈立伟，等.西安地裂缝场地地震效应分析[J].工程地质学报，2012，20（6）：979985.
LIU Cong，PENG Jianbing，CHEN Liwei，et al.Finite Element Analysis of Seismic Effects on Ground Fissures in Xian[J].Journal of Engineering Geology，2012，20（6）：979985.
[13]魏加华，崔亚莉，邵景力，等.济宁市地下水与地面沉降三维有限元模拟[J].长春科技大学学报，2000，30（4）：376380.
WEI Jiahua，CUI Yali，SHAO Jingli，et al.3Dfinite Element Method Simulation of Groundwater and Land Subsidence in Jining City[J].Journal of Changchun University of Science and Technology，2000，30（4）：376380.
[14]万伟峰.西安市地下水开采地面沉降数值模拟及防治方案研究[D].西安：长安大学，2008.
WAN Weifeng.Numerical Simulation of Land Subsidence Due to Groundwater Withdraw in Xian and Its Prevention Schemes[D].Xian：Changan University，2008.
[15]OH H J，AHN S C，CHOI J K，et al.Sensitivity Analysis for the GISbased Mapping of the Ground Subsidence Hazard near Abandoned Underground Coal Mines[J].Environmental Earth Sciences，2011，64（2）：347358.
[16]SUN Y H，ZHANG X D，MAO W F，et al.Mechanism and Stability Evaluation of Goaf Ground Subsidence in the Third Mining Area in Gong Changling District，China[J].Arabian Journal of Geosciences，2015，8（2）：639646.
[17]KIM K D，LEE S，OH H J.Prediction of Ground Subsidence in Samcheok City，Korea Using Artificial Neural Networks and GIS[J].Environmental Geology，2009，58（1）：7172.
[18]OUYANG Z H.The Role of Potential Soil Cavity on Ground Subsidence and Collapse in Coal Mining Area[J].Journal of Coal Science and Engineering，2010，16（3）：240245.
[19]张永伟.华北平原德州地面沉降成生机理、监测预警与可控性研究[D].济南：山东大学，2014.
ZHANG Yongwei.Formation Mechanism，Monitoring and Warning，Controlling Research of Subsidence of Dezhou in North China Plain[D].Jinan：Shandong University，2014.
[20]贾超，张国荣，何奔，等.华北平原德州地面沉降监测预警系统建设报告[R].济南：山东大学，2014.
JIA Chao，ZHANG Guorong，HE Ben，et al.Land Subsidence Monitoring and Warning System Construction Report in Dezhou District of North China Plain[R].Jinan：Shandong University，2014.
