华北平原地下水位驱动下的地面沉降现状与研究展望

郭海朋，李文鹏，王丽亚，陈 晔，臧西胜，王云龙，朱菊艳，卞跃跃

（1.京津冀平原地下水与地面沉降野外科学观测研究站，北京 100081；2.中国地质环境监测院，北京100081；3.北京市水文地质工程地质大队，北京 100195；4.中国地质大学（北京）水资源与环境学院，北京 100083；5.中国地质博物馆，北京 100034）

地面沉降是世界各地普遍发生的环境地质问题，而地下水动力与土层变形的互馈关系研究是地面沉降区地下水资源量评价的前提，是区域地面沉降预测预警的基础，也是环境地质学科发展的需求。目前国际上关于地下水动力与土层变形互馈机制的研究，主要集中在地下水位驱动下黏性土和砂土变形规律研究，时序沉降数据用于计算区域水文地质参数和校正模型参数，以及土-水耦合模型研究三个方面。华北平原是我国地面沉降影响面积最大的区域，地面沉降对京津冀协同发展、雄安新区、北京城市副中心等国家重大战略区规划建设和南水北调、高速铁路等国家重大工程区安全运营存在重要影响。地下水位变化与土层变形互馈关系研究，对华北平原地面沉降监测网的优化、监测数据的定量化分析、地面沉降模型仿真性和预测预警准确性的提高等方面具有科学意义，可为有效缓解该地区地下水利用引起的地面沉降灾害提供技术支撑。

1 国内外研究现状

1.1 地下水位驱动下黏土和砂土变形规律研究

近几十年来，含水层系统中地下水位变化影响下的土层变形规律受到国内外许多学者的高度关注。例如，He等[1]实施了物理模型试验，研究地下水补给和排泄引起的砂层和黏土层变形特征及孔隙水压力变化规律，提出了三种有利于减缓地面沉降发展的地下水开发模式。郭海朋等[2]初步总结了天津平原地下水水位变化模式，分析了不同地下水位变化模式下的土层变形特征。总体而言，该领域研究成果主要集中在土层尤其是黏性土层的流变变形特征研究和土层变形应力应变关系定量化研究等方面。

1.1.1 地下水位驱动下土层变形流变特征研究

随着地面沉降问题研究的深入，国内外越来越多研究者逐渐认识到抽水诱发下土层的压缩量中包含了一定比例的流变变形量，不仅黏性土可能发生蠕变[3-5]，砂砾含水层也可能发生蠕变[6-7]。过去在研究地面沉降时，通常将含水砂土层作为线弹性体处理，或者认为砂土层的变形对地面总沉降的贡献可以忽略不计[8]。试验结果表明在单调荷载下，砂土的蠕变变形可以达到总变形的10%[9]。分层标监测数据也表明，含水砂层的变形不仅包括弹性变形，而且包括塑性变形和蠕变变形[2,10]。研究人员根据长三角地区大量地面沉降和水位观测资料，结合室内土工试验，揭示了砂土层在不同的应力条件下，有的表现为弹性变形，有的则表现为以塑性变形为主并包含有蠕变的非线性变形特征[10-12]。李玉岐等[13]利用自行设计的砂土排灌水试验装置，分析地下水排灌引起的砂样宏观竖向变形及细观移动，结果表明：在砂样排灌水的初期阶段，砂样的结构发生了明显重组，砂样沉降不仅发生在排水时，而且在回灌时继续增加，砂样产生了较大的、不可恢复的塑性和黏性变形；砂样结构达到相对稳定后，排水时产生的竖向变形变小，而回灌时砂样发生较大的回弹。郭海朋等[2]通过土体高压固结试验表明冀中坳陷内武清凹陷地区及沧县隆起地区中晚更新世地层黏性土体在不同荷载条件下，蠕变特征明显。

1.1.2 地下水位驱动下土层变形本构关系研究

Terzaghi固结理论是简单的弹性变形模型，由于机理简单和控制方程的线性特征，被广泛应用于土层压缩分析。考虑到过于复杂的黏弹塑性应力-应变模型很难直接用于区域地面沉降模型中，Corapcioglu等[14]基于划分主、次固结的黏土固结理论，提出了一维线性黏弹本构关系（Merchant模型）。吴林高等[15]在分析抽灌水引起地面沉降的力学机制基础上，得出了广义Kelvin本构模型。以上两个模型的适用性都缺乏试验数据验证。Ye等[16]通过分析长三角地区土层变形特征，采用修正的Merchant模型刻画土层的变形规律，在此基础上建立了三维变参数水流和垂向一维沉降组成的部分耦合模型，对长三角区域的地面沉降进行了模拟。如果能获取土层变形的野外监测和室内试验数据，可能会建立更加贴近实际的土层变形本构模型并用于地面沉降模拟和预测[17]。Tsai等[18]研究提出了一种黏弹塑性模型，该模型由一个黏弹性模型和一个黏塑性模型组成，并由单刀双掷开关控制黏塑性模型的使用，经验证，该模型计算结果具有较小的相对误差。

1.2 时序沉降数据用于计算区域水文地质参数和校正模型参数

利用分层标时序土层变形和地下水位数据反演计算区域水文地质参数已经为大家所熟悉，应力-应变图解法是常用的手段之一。该方法最早由Riley[19]提出，根据应力应变关系曲线可求出弹性贮水率和非弹性贮水率，Cleveland等[20]提出的图解法可进一步确定垂向渗透系数。叶淑君等[21]结合上海沉降和水位观测资料，求出了上海地区弱透水层的弹性贮水率、非弹性贮水率和垂向渗透系数，为上海地面沉降的数值模拟提供了较好的初始参数值。Zhang等[22]利用北京平原某分层标监测数据计算了弹性贮水率和非弹性贮水率，并与美国加利福尼亚州圣荷昆谷（San Joaquin Valley）的对应参数值进行了相互比较和评价。

区域地面沉降监测技术在过去的几十年经历了较快的发展，其主要监测手段从经典的水准测量发展至日趋成熟的GPS技术与合成孔径雷达干涉测量技术（Interferometric Synthetic Aperture Radar，InSAR），获取了大量的区域地表形变数据[23-25]。水文地质参数是地面沉降数值模拟过程中必不可少的信息，但由于含水层的非均质等复杂性质，许多水文地质参数的空间分布在现有技术条件下无法快速准确获取，因而利用一些易获取且精度较好的监测数据（如地下水位等）推估难以获取的水文地质参数分布成为了地下水领域的研究热点[26]。当含水层系统相对简单，不存在地下水位变化和地面沉降之间明显滞后的情况下，区域地面沉降和地下水位监测数据可以用来直接反演计算弹性贮水率、非弹性贮水率等参数，进而可以补齐时间或空间上地下水位观测数据的缺失[27-29]。Zhuang等[30]探索了利用地质统计方法反演求取弱透水层渗透系数、弹性和非弹性贮水率的方法，并成功用于常州某地弱透水层水文地质参数的计算。Zhuang等[31]基于欧拉公式，获得了地下水位变化驱动下超固结弱透水层弹性变形的解析解，并用来评估上海某35.54 m厚弱透水层的垂向渗透系数和弹性储（释）水率，结果表明该解析方法可以定量解释地面沉降滞后特征并有效评估超固结弱透水层的水力参数。

高分辨率水文地质参数是地面沉降模型实现准确预测的基础，而InSAR、水准、GPS等多种手段获取的高时空分辨率地面沉降监测数据使有效校正模型的水文地质参数成为可能。UCODE[32]是一种基于梯度的自动反演算法，研究人员已经开展了利用InSAR解译计算的地面沉降数据和地下水位监测数据并基于UCODE进行模型水文地质参数校正的研究[27,33-34]。Zhang等[34]基于UCODE并利用InSAR获取的高时空分辨率地面沉降数据和地下水位数据校正了含水层导水系数、弹性和非弹性贮水率、断层导水系数等水文地质参数，结果表明联合利用地面沉降和地下水位监测数据校正的模型比仅用地下水位数据校正的模型更符合实际。但是，基于梯度的反演方法主要缺点是计算量大，不能快捷地评估参数，而且预测具有不确定性。

1.3 土-水耦合模型研究

抽排水引起的地面沉降问题实际上就是一个渗流场与应力场耦合的问题，土-水耦合模型按渗流场与应力场结合方式分为不耦合的两步计算模型、部分耦合模型和完全耦合模型。不耦合的两步计算模型分为完全独立的两步，先计算孔隙水压力，再计算变形，水流及沉降方程中各参数不随时间变化。两步计算模型最初在研究威尼斯的由多层含水层与弱透水层组成的地下水系统因抽水引起的地面沉降问题时提出[35]。Shearer[8]以太沙基一维固结理论为基础，在Modflow程序的基础上开发了考虑黏性土层孔隙水压力变化滞后于含水层水位变化的夹层排水模块（IDP），通过沉降计算获取弱透水层的释水，实现了大厚度含水层释水压缩的模拟。部分耦合模型中孔隙水压力和变形分步计算，但两者之间相互影响。含水层水头下降会导致相邻弱透水层中的地下水产生渗流，弱透水层的变形具有明显的非线性特征；同时，土层的孔隙比、压缩性和透水性等参数也随着土体变形量的变化而变化。国内外地面沉降模型大多基于土-水部分耦合模型[36-39]。完全耦合模型基于Biot固结理论，孔隙水压力和土层变形同时计算，在理论上这种模型最符合沉降物理机制，但是物理场控制方程更为复杂、计算量大、需要参数多，很难在实际问题中应用。近年来研究人员在完全耦合模型的数值计算方法方面进行了探索[40-41]，并将完全耦合模型成功应用于意大利威尼斯、美国拉斯维加斯和中国上海、沧州、德州等地的地面沉降机理研究及模拟预测[42-46]。骆勇等[47]分别采用基于渗流-沉降分步计算的土力学经验公式、渗流-沉降部分耦合的GMS软件中SUB模型和基于渗流-沉降完全耦合的COMSOL Multiphysics模型，对疏排水引起的地面沉降进行了模拟计算，研究结果表明，完全耦合模型计算结果更为合理、更符合实际沉降特征。Wang等[48]探索研究了随机非均质多孔介质中地下水流动和土层变形的动态相互作用。Pham等[49]基于多孔介质弹性力学理论，开发了程序模块SUB+，可以实现水-应力的完全耦合模拟。TOUGH2是一套功能强大、应用广泛的模拟孔隙或裂隙介质中多相流的系列程序，在水气两相流数值模拟中得到了广泛应用[50-51]。研究人员以TOUGH2为基础进行二次开发[52]或者将其与其它地质力学模拟程序结合[53]，实现了地下水流与地质力学的耦合模拟，显示出该程序在土-水耦合模拟方面有广泛的应用前景。

2 华北平原地面沉降现状

我国地面沉降主要发生在华北平原、长江三角洲、汾渭盆地、淮河平原、珠江三角洲、东北平原以及山区断陷盆地等地区，共有22个省（区、市）超过100个地级以上城市发生地面沉降。其中，华北平原是目前我国地面沉降影响面积最大、沉降速率最快的区域[54]。华北平原人均水资源量296 m3，为全国人均水资源的1/8，与干旱缺水国家以色列的人均水资源量相当，属于水资源极度短缺地区。因地表水严重短缺，生产生活供水主要依赖地下水[55]，2000年以来，地下水占总供水量的比例达到64%。南水北调工程一定程度上缓解了供水紧张局面，但依然不能满足社会经济快速发展需求。地下水长期超采形成多个地下水降落漏斗，并引发严重的地面沉降，累计沉降量较大地区与深层地下水漏斗区的分布基本吻合，天津、沧州、衡水、德州地面沉降已经连成一片（图1）。依据累计地面沉降量和地面沉降速率两个指标（表1），对京津冀平原地面沉降现状发育程度进行评价，结果表明地面沉降发育程度强的区域主要分布在北京平原东部、天津平原中南部、河北平原沧州、衡水、廊坊、邢台、邯郸等地，总面积达1.4×104km2（图2）。

图1 京津冀平原主要浅层和深层地下水漏斗以及累计沉降量较大地区（大于50 cm）分布图Fig.1 Distribution of shallow and deep groundwater cones and the area with the cumulative subsidence greater than 50 cm in Beijing—Tianjin—Hebei Plain

表1 地面沉降现状发育程度评价标准Table 1 Evaluation standard of the development of land subsidence

图2 京津冀平原地面沉降现状发育程度分区图Fig.2 Zoning map of land subsidence development in Beijing—Tianjin—Hebei Plain

分层标监测数据显示，北京市100 m以下深部地层沉降贡献率为82%，其中天竺、张家湾等地深部地层贡献率超过90%，浅部部分地层出现回弹。天津市300 m以下深部地层贡献率占65%，其中汉沽营城、西青郑庄子等地深部地层贡献率超过70%，浅部部分地层出现回弹。河北省实施深层地下水禁限采的衡水、廊坊等地区，深部地层对地面沉降的贡献率逐渐减少，其中衡水市区150 m以深地层的沉降贡献率由2009—2015年的66%降至2019年的34%。必须高度重视地下水尤其是深层地下水超采问题，开展地面沉降防控，减少经济损失。

近年来，南水北调工程运行为地下水超采治理提供了置换水源，减少了地下水开采量，工程沿线大中型城市区地下水位下降趋势得到遏制，地面沉降防控取得了初步成效，天津、沧州、衡水等重点城市主城区地面沉降得到有效控制，北京、天津地面沉降严重区面积（本文地面沉降严重区指年沉降量大于50 mm的地区）总体呈现减少趋势（图3）。南水北调进京5年以来，北京地下水位整体处于恢复上升状态，尤其是100 m以浅的含水层，水位回升幅度较大，与2014年相比，平原区地下水位平均上升了3 m，地下水降落漏斗面积减少了38%（图4）。地下水位恢复导致地面沉降减缓，地面沉降严重区面积减少了71%。沧州市自2005年开始实施关停单位自备井的禁采措施，市区年沉降量已由60～80 mm降至目前的10～15 mm。沧州市区分层标组于2010年6月开始监测，共由五座分层标组成，分别监测第一压缩层（5～69 m）、第二压缩层（69～196 m）、第三压缩层（196～253 m）、第四压缩层（253～375 m）。2014年第一至四压缩层的变形量分别为−2.4，−5.4，−4.6，−8.9 mm（负号表示压缩，正号表示回弹，下同），2019年地面沉降进一步减缓（图5），四个压缩层的变形量分别为+0.6，−2.7，−2.2，+4.2 mm。2013—2018年，河北衡水市主城区年平均沉降量由15 mm减至10 mm。衡水分层标主要监测第一压缩层（41～150 m）、第二压缩层（150～267 m）、第三压缩层（267～401 m），其中第一压缩层对应浅层含水层组，第二、三压缩层对应深层含水层组。2009年10月—2015年8月，第一至三压缩层的平均变形速率为−19.3，−31.3，−6.1 mm/a，对地面沉降量的贡献率分别为34%、55%、11%。衡水地区从2014年开始进行地下水超采区综合治理，大幅压减地下水开采量，衡水市城区禁采深层地下水，深层地下水水位下降减缓甚至回升，有效遏制了地面沉降的快速发展。2019年第一至三压缩层的变形量减至−10.9，−3.4，−2.1 mm，对地面沉降量的贡献率为66.2%、20.8%、13%，主要压缩层位从之前的以第三、四含水层组为主逐步转变为以一、二含水层组为主。

图3 北京、天津和河北平原地面沉降严重区面积占各省（市）平原区面积比例（2012—2019年）Fig.3 The ratio between the severe subsidence area and the total area for Beijing,Tianjin and Hebei plain（from 2012 to 2019）

图4 北京平原地面沉降严重区、地下水漏斗面积占平原区总面积比例与平均水位埋深（2012—2019年）Fig.4 The ratio of area of severe subsidence and groundwater depression cones to the total plain area and the average waterlevel depth for Beijing plain（from 2012 to 2019）

图5 沧州市区各压缩层累计变形随时间变化曲线（2019年）Fig.5 The temporal variation of the deformation of different compression layers in the downtown of Cangzhou（2019）

虽然华北平原地面沉降出现减缓态势，北京、天津、河北沧州等治理区地面沉降速率大幅下降，但华北平原尤其是河北平原地面沉降总体上仍然处于较快发展阶段，地面沉降防控形势依然严峻。为满足农业用水需求，山前平原超采浅层地下水，中东部平原的衡水、沧州等地超采深层地下水，农业区超采是当前华北平原地面沉降严重的主要影响因素之一。2019年河北全省平均降水量437.6 mm，比2018年减少71.6 mm，为确保农业生产，加大了灌溉期地下水开采力度，使农灌区地下水位普遍下降。据中国地质调查局地下水位统测数据，2019—2020年，华北平原浅层和深层地下水位平均分别下降0.6 m和1.0 m，水位下降区主要分布在农业灌区，加剧了地面沉降发展。目前，华北平原地面沉降严重区面积超过全国沉降严重区总面积的80%，沉降速率较快地区主要分布在北京朝阳与通州的交界处，天津武清、静海、西青、北辰，河北廊坊、沧州、保定、衡水、邯郸等地。河北雄安新区、北京城市副中心等地存在地面沉降严重区，京津城际、京沪高铁等高速铁路穿越地面沉降严重区，应予以高度重视。

3 研究展望

地下水位变化驱动下的土层变形特征及其机制研究仍属国内外地面沉降研究的前沿课题。如前所述，土层变形规律、水文地质参数反演校正和土水耦合模型应用等三个研究方向联成一个统一的整体，近年来在基础理论上有了较大的发展，在技术方法上更加注重利用InSAR等新技术、新方法获取的数据开展研究，为华北平原地面沉降防治研究提供了很好的经验借鉴。同时，华北平原具有独特的水文地质条件，是我国含水层结构最为复杂的大型平原，受长期地下水超采等人类活动影响，形成了世界上面积最大的地下水降落漏斗和最严重的地面沉降区。不同地区地质条件的复杂性和时空上的差异性导致许多研究结论可能无法在该地区适用，需要不断提高地下水与地面沉降监测技术水平与监测精度，因地制宜地开展相关研究工作。在总结地下水位变化模式和土层变形规律的基础上，研究地面沉降差异性特征及变形成因机理，将大幅提升华北平原地面沉降机理的认识水平，极大促进环境地质学科发展，为地面沉降防灾减灾提供重要技术支撑。

（1）深化地面沉降机理和预测预警研究。华北平原构建了地面沉降监测网，监测面积达9×104km2，覆盖了京津冀整个平原区；其中，建成基岩标14个，分层标31组，水准监测点3 677个，GPS监测点272个，地下水监测井5 000余眼，InSAR实现了京津冀平原区全覆盖连续动态监测。由于沉积过程、沉积相、地层岩性及应力加载过程条件（地下水位变化模式）等不同，华北平原不同埋深土层在地下水位变化影响下的变形特征存在较大差异。多年积累的地下水位、土层变形等信息对于地面沉降预测预警来说无疑十分宝贵，可是迄今对这些数据缺乏科学化、规律化的分析和提炼，没有对不同地下水位变化模式下土层变形规律这一关键问题进行系统的研究。如何联合利用室内土工试验和基于现场调查监测数据的数理分析揭示土层变形应力应变规律？如何量化这一规律从而形成简单实用的土层变形本构关系，甄别诱发地面沉降快速发展的地下水位及其下降速率等地面沉降控制指标？如何将土层变形本构模型与地下水流模型耦合用于地面沉降监测预警？对这些科学问题的解决，直接关系到对华北平原地面沉降形成机理的认识，也是区域地面沉降预测预警的迫切需求，有助于因地制宜地提出有针对性的地面沉降防控措施。图6概括了地下水位变化驱动下土层变形规律的研究思路。

（3）提升对华北平原地下水资源属性的认识。地下水动力与土层变形的互馈关系研究是地面沉降区地下水资源量评价的前提条件。以华北中东部平原为例，深层地下水补给条件较差，可更新能力低，在相同地下水开采强度下，更容易产生地面沉降。近40年以来，深层地下水位最大累计降幅已经超过100 m，目前天津、沧州、衡水、德州深层地下水漏斗已经连成一片，面积超过2.0×104km2。华北平原深层地下水开采量中含水层和弱透水土层压密释水量所占比例达40.1%～43.8%，沧州地区高达57.6%[56]。正确评价地面沉降区地下水资源量，需要充分利用地面沉降观测数据，科学确定沉降层位、沉降量与地下水开采量的定量关系。

图6 地下水位变化驱动下土层变形规律研究思路Fig.6 Research ideas on groundwater-derived deformation law of soil layers

（4）加强地热开发与地面沉降关系研究。京津冀地区地热资源开发始于20世纪70年代初，21世纪以来开发利用规模不断扩大，目前是我国地热资源开发利用程度最高的地区，主要开发新近系明化镇组、馆陶组、古近系东营组3个孔隙热储层和寒武-奥陶系、蓟县系2个岩溶热储层。地热水开发利用对地面沉降的影响目前尚无定论，主要原因是缺少观测数据和研究成果支持。深入研究地热水开发引起的水位变化与土层变形的关系，提出合理的水位控制指标，可以为实现地热资源持续利用和地面沉降有效控制提供理论依据和数据支撑。