-
地面沉降是在自然和人为因素作用下，由于地表松散未固结土层压缩而导致的区域性地面高程降低的地质现象，严重时会形成灾害（Hu Ruilin et al.，2004; 薛禹群等，2006; Motagh et al.，2008; 郭海朋等，2017）。目前，世界范围内有34个国家，大约200个城市相继发生地面沉降（Herrera-García et al.，2021），如美国的圣华金流域和圣克拉拉流域（Galloway et al.，1998；Pavelko et al.，2006；Jeanne et al.，2019）、日本的东京和大阪（Cao et al.，2020）、墨西哥的墨西哥城（Strozzi et al.，1999; Chaussard et al.，2014）、意大利的威尼斯（Tosi et al.，2007; Teatini et al.，2010）、泰国的曼谷（Phien-Wej et al.，2006）、印度尼西亚的雅加达（Abidin et al.，2008，2011）、中国的上海、天津、北京等（Hu Beibei et al.，2009; Xu Yeshuang et al.，2012; Gong Huili et al.，2018）。据大量研究表明（Galloway et al.，2011；Mahmoudpour et al.，2016；Ye Shujun et al.，2016；Chen Beibei et al.，2019；狄胜同等，2020；Zhu Lin et al.，2020；郭海朋等，2021），地下资源的不合理开发利用，尤其是地下水开采不当是导致地面沉降及次生灾害的主要诱因。截至2019年，我国因超量开采地下水导致的地面沉降严重区（年沉降量大于50 mm）面积达1.14×104 km2，主要分布在华北平原、长江三角洲、汾渭盆地和淮北平原（李文鹏等，2020）。地面沉降的发生、发展造成大量基础设施变形破坏，增加城市洪涝风险，诱发地裂缝，并且永久性地损失地下含水系统的储水能力，已成为全球性、综合性、对人类生存环境产生严重影响的地质环境问题（Herrera-García et al.，2021）。
-
北京是国际上为数不多的以地下水作为主要供水水源的特大型城市之一。人均水资源年占有量不足300 m3，仅是全国人均水平的1/8，是世界人均水平的1/32，与干旱缺水国家以色列的人均水资源量相当，属于水资源极度短缺地区（贾三满等，2018；郭海朋等，2021）。因地表水严重短缺，北京市的生产生活供水主要依赖于地下水。自2000年以来，地下水开采量占全市总供水量的比例保持在65%以上（雷坤超等，2022a）。在此期间，北京市建成4处应急水源地，并出台了最严格水资源管理制度，不断加强用水管控和节水力度，虽然在一定程度上缓解了供水紧张的局面，但依然不能满足社会经济快速发展的需求。地下水长期超采形成多个地下水位降落漏斗（王树芳等，2019），并引发严重的地面沉降，累计沉降量较大的地区与深层地下水位降落漏斗区基本吻合，多个沉降中心连片发展，历史最大沉降速率达159 mm/a（雷坤超等，2016，2022b；Chen Beibei et al.，2019，2020）。
-
2014年12月，南水北调中线一期工程正式通水，标志着江水正式进京。南水北调工程是我国迄今为止建造的最大水利工程，旨在通过南水北调来缓解我国中、北方水资源短缺问题，实现水资源的优化再分配。南水北调中线工程全长约1300 km，起点为丹江口水库（N 32° 43′，E 111° 34′），横跨河南、河北、北京和天津4个省（市），终点至北京市颐和园团城湖（N 39°55′，E 116°24′）。2014~2020年，北京市累计接收南水超过60.0×108 m3，按照“喝、存、补”的原则，约有40.0×108 m3用于城市水厂供水，8.0×108 m3存入大中型水库，12.0×108 m3用于回补地下水。其中，回补地下水主要采用两种方式：第一，通过潮白河冲洪积扇中上部河道和砂石坑入渗回补地下水源地；第二，通过密云、怀柔等水库汛期放水进入河道，向下游补给地下水（雷坤超等，2020）。与此同时，北京市相继开展了地下水禁、限采，城区自备井置换等措施，全市平原区平均地下水位上升了3.6 m（北京市水资源公报，2020），主要沉降区面积和沉降中心速率均呈现出逐渐减小的趋势（雷坤超等，2022a）。
-
截至目前，已有较多的学者围绕南水北调工程对沿线地下水环境和生态变化（Yang Yong et al.，2012；杨庆等，2017; Zhao Yong et al.，2017），以及地下水位变化（Zhang Yuan et al.，2014；Ye Aizhong et al.，2014; Li Xue et al.，2017; Zhang Menglin et al.，2018; Chen Beibei et al.，2020）等方面的影响开展了较为丰富的研究。同时，部分学者利用合成孔径雷达干涉测量（InSAR）技术分析了南水北调工程运行前后地面沉降的变化状况（Chen Beibei et al.，2020；Zhu Lin et al.，2020; Dong Jie et al.，2021; 张双成等，2022; 雷坤超等，2022a）。随着南水正式进京，北京市水情发生了明显变化，许多地区的地下水位由持续下降转为逐渐上升，然而综合利用多种监测数据，系统分析南水北调工程运行前后北京平原区地下水流场、地面沉降、土层变形特征以及相关水文地质参数变化的研究鲜有报道。因此，本文在充分收集、整理多种监测数据基础上，分析南水北调中线工程运行前后北京平原区地下水位和地面沉降的变化；研究不同水位变化模式下不同岩性及深度土层的变形特征；计算土层不同变形阶段的弹性和非弹性储水率；并对黏性土层产生较大残余变形和滞后变形的原因进行了探讨。研究成果对查明南水北调工程对北京地质环境的影响，进一步厘定土层变形本构关系，实现南水进京后地面沉降的精确模拟和地下水的科学调控具有重要意义。
-
1 研究区概况
-
1.1 区域地质条件
-
北京位于华北平原的西北边缘，除东南局部地区与天津接壤外，其余地区均与河北省毗邻。地形西北高，东南低，西部为太行山脉，北部为燕山山脉，东南为缓倾斜平原。地理坐标：E 115.7°~117.4°，N 39.4°~41.6°。北京总面积16410.54 km2，平原区面积6400 km2，约占全市总面积的39%。北京属典型暖温带半湿润半干旱大陆性季风气候，平均气温11.7℃左右，夏季最高可达42.6℃。多年平均降水量570 mm/a（1961~2020年），其中80%的降水主要集中在六月中旬至九月。
-
北京平原区第四纪地层主要由永定河、潮白河、北运河、大清河和蓟运河5大河流的冲、洪（湖）积作用形成。在山前和冲洪积扇顶部，第四系厚度在20~40 m左右。岩性为单一的砂、砂砾石层或顶部覆盖较薄的黏性土层。在冲洪积扇中下部和冲积平原地区，沉积物厚度逐渐增大，层次增多，颗粒逐渐由粗变细，岩性过渡为砂、砂砾石与黏性土层交错出现，并以黏性土为主（图1）。在部分沉积凹陷中心，第四系厚度达1000 m左右（蔡向民等，2009）。
-
图1 北京平原区位置与第四系厚度
-
Fig.1 The location and Quaternary thickness of Beijing Plain
-
根据平原区第四纪地层的形成时代、沉积环境等，可将平原区划分为3个主要的压缩层组。第一压缩层组（Q4+Q3），底板埋深小于100 m，广泛分布于平原区。岩性主要为全新统和上更新统冲积、冲湖积的粉土和黏性土，可塑，压缩性中等。第二压缩层组（Q2），底板埋深150~180 m左右，主要分布在冲洪积扇中下部、冲积平原地区。岩性主要为中更新统冲洪积、冲湖积粉土和黏性土层，可塑，压缩性中-低。第三压缩层组（Q1），底板为第四系基底，部分地区在260~300 m左右，主要分布在几大沉积凹陷地区（图2）。岩性主要为下更新统河湖相沉积的黏性土层，可塑-硬塑，压缩性低，大部分呈坚硬状态（雷坤超等，2016）。
-
1.2 地下水开采与地面沉降历史
-
北京地面沉降的发生、发展与地下水开发利用历史具有明显的对应关系（图3和图4）。1955~1973年，为地下水初步开发和地面沉降形成阶段。这一时期，北京地下水开采量较小，年均开采量约10.79×108 m3。地下水采补基本平衡，仅在城近郊局部地区出现地下水超采现象。地面沉降也仅发生在地下水超采的东郊和东北郊地区，地面沉降速率每年几毫米到十几毫米。
-
1973~1983年，是地下水开采增加和地面沉降发展阶段。这一时期，北京年均地下水开采量维持在25×108~28×108 m3之间，地下水呈负均衡，城近郊区严重超采。形成东郊大郊亭和东北郊来广营两个沉降中心，沉降速率在18~30 mm/a。
-
1983~1999年，是地下水平稳开采和地面沉降扩展阶段。这一时期，地下水开采相对稳定，年均开采量约26×108 m3，地下水累计储变量持续亏损。东郊地区地下水开采量减少，其他区县开采量增加。超采区占平原区面积的70%左右。此阶段地面沉降区面积扩展明显，形成“北”和“南”两大沉降区，五个主要沉降中心，沉降速率19~24 mm/a。
-
1999~2014年底，是地下水严重超采和地面沉降快速发展阶段。这一时期，北京遭遇了史上最长的连续枯水年。为保障城市供水，北京市采取了有效的开源节水措施。此阶段，地下水开采量稳中有降，但开采深度向深层发展，加之连续干旱，北京地下水资源量持续亏损（累计亏损70.90×108m3），超采区占平原区面积的90%以上。该时期，北京地面沉降快速发展，在以往沉降中心基础上，逐渐形成七个主要沉降中心，最大沉降速率达159 mm/a（贾三满等，2018）。
-
2015年至今，是地下水开采减少和地面沉降减缓阶段。这一时期，随着南水正式进京，北京市相继开展了地下水禁、限采，自备井置换和地下水回补等一系列措施，地下水开采量逐年减小，由2014年的19.6×108m3下降至2020年的13.5×108m3。另外，2015~2018年，北京市经历了连续的丰水年，年均降水量达619.0 mm，较多年平均降水量多49.0 mm。因此，在上述多种因素综合影响下，北京平原区地下水水位平均回升了3.6 m（北京市水资源公报，2020），地下水资源量出现明显盈余（累计盈余约12.8×108m3）。在此阶段，北京市地面沉降逐渐减缓，部分沉降中心消失，最大沉降速率为86.0 mm/a （雷坤超等，2020）。
-
图2 北京平原区A—A’处水文地质剖面（位置见图1）
-
Fig.2 Hydrogeological cross-section A—A’ in Beijing Plain (the location is indicated in Fig.1)
-
图3 北京市降水量、地下水累计储变量与地面沉降发展关系
-
Fig.3 The relationship between the precipitation, the change of groundwater accumulated storage and the land subsidence in Beijing
-
图4 北京市地下水开采量与平均沉降速率关系
-
Fig.4 The relationship between groundwater exploitation and average subsidence rate in Beijing
-
2 数据选取与研究方法
-
2.1 数据选取
-
2.1.1 RADARSAT-2数据
-
RADARSAT-2由加拿大太空署与MDA公司合作运营，是一颗搭载C波段传感器的高分辨率商用雷达卫星。可以对地表进行全天时、全天候观测，具有较高的分辨率和重访周期，数据质量较高。本次利用永久散射体差分干涉测量（PS-InSAR）技术，对覆盖北京平原区72景RADARSAT-2雷达卫星 SAR影像（Wide strip模式，C波段，幅宽150 km 分辨率30 m×30 m，入射角为27.8°，2013-01-16~2020-12-22），进行时序差分干涉处理及PS点提取，剔除了基线距大于300 m的干涉像对，利用30 m分辨率SRTM数据作为参考DEM，进行地形相位去除，经地理编码及视线向（LOS）到垂直向投影变换后，获取了北京平原区2013~2020年时序地面沉降分布图。
-
2.1.2 分层标和水准测量数据
-
文中选取平原区7个地面沉降监测站全部55个分层标和37眼分层地下水位监测井数据，时间序列为2006~2020年或2009~2020年，分析不同深度土层变形与地下水位变化之间的响应关系，查明南水进京前后地下水位变动条件下土层应力-应变特性，计算相应水文地质参数。以天竺站为例，分层标和地下水位监测井观测层位见图5所示。同时，利用2017年20个区域一等水准测量数据对PS-InSAR结果进行精度验证（位置见图1）。
-
2.1.3 区域地下水位动态监测数据
-
文中利用平原区2015年和2020年地下水位等值线数据，绘制了平原区2015年和2020年不同含水层组地下水位同期对比图。选取平原区（北部、东部和东南部地区）15眼地下水井长时序监测数据，分析南水进京前后不同含水层组地下水位变化趋势。地下水井位置如图1所示。其中，井1、井2、井3位于平原区北部，井深分别为45 m、60 m和89 m。井4、井5、井6位于平原区北部潮白河河道附近，井深分别为50 m、100 m和60 m。平原区北部含水层组主要为单一的砂砾石层。D1-1、D1-2、D-3、D-4、D1-5井位于平原区东部的王四营地面沉降监测站，井深分别为185 m、155 m、102 m、54 m和19 m。井7、井8、井9、井10位于平原区东南部，井深分别为182 m、300 m、166 m和300 m。平原区东部和东南部含水层组主要为砂层与黏性土层组成的多层结构。
-
2.2 研究方法
-
2.2.1 PS-InSAR方法
-
PS-InSAR技术最早由Ferretti et al.（2000，2001）提出，其基本思想是在SAR 图像的长期序列中选择高相干像素点作为研究对象，通过分解每个相干点的相位信息，包括地形、轨道和大气等的相位变化，进而得到地表变形。该方法能有效克服DInSAR时空失相干和大气延迟的影响，提高地表变形监测精度，可达毫米级。主要处理步骤如下（雷坤超等，2016）：① SAR影像配准。根据SAR数据时空基线分布，确定主影像，对其他辅影像进行坐标变换和重采样，完成影像配准。② 差分干涉处理。根据基线分布，将所有满足条件的干涉像对进行差分干涉处理，生成差分干涉纹图序列。③ 相干目标选取。根据序列SAR影像的幅度稳定性和相干系数变化，选择相干目标候选点（PS点）。④ 地表形变速率提取。对解缠相位进行滤波处理，估算主、辅影像轨道误差、大气效应误差，并去除地形残余相位等信息，进而获得各PS点形变速率。同时，利用一等水准测量数据对PS-InSAR结果进行验证。以水准点附近100 m为缓冲区，提取缓冲区内全部PS点并求平均值。将水准测量结果与该平均值做差处理，其差值范围为2~8 mm，互差的均方误差为5.04 mm（表1）。
-
图5 天竺地面沉降监测站内分层标和地下水位监测井埋设位置
-
Fig.5 Observation horizon of extensometer and groundwater level monitoring well of Tianzhu land subsidence station
-
（a）—分层标和地下水位观测层位；（b）—监测站标房；（c）—分层标标底
-
（a）—observation horizon of extensometer and groundwater well; （b）—interior view of monitoring station; （c）—bottom of the boreholeextensometer
-
2.2.2 储水率和滞后时间估算方法
-
储水率是描述含水层储存或释出水能力大小的一个非常重要指标。其定义为地下水位抬升或降低一个单位，从单位体积土层中储存或释放的水量。由于上述储水率概念仅能定量描述线弹性变形土层，不能刻画土层普遍存在塑性变形，所以储水率的概念被扩展，把储水率分为弹性和非弹性（塑性）（Helm，1975）。如果含水层的有效应力小于先期固结压力，含水层呈现可恢复的弹性变形。相反，如果含水层的有效应力大于先期固结压力，含水层表现出不可逆的非弹性变形，即塑性变形。由于先期固结压力是土层历史上最大的有效应力，可以对应含水层的历史最低水位，可以通过监测数据来确定。因此，如果地下水位低于先期固结压力所对应的水位，储水率称为非弹性储水率；否则称为弹性储水率。值得注意的是，对于土层变形的不同阶段，骨架储水率的变化并非常数，其取值取决于有效应力与先期固结压力的关系。
-
式中，Sske为骨架弹性储水率（m-1），Sskv为骨架非弹性（塑性）储水率（m-1），为变形量（m），为水位变化量（m），b0为土层初始厚度（m），h为地下水水位（m），hmin为先期历史最低水位（m）。对于黏性土的释水滞后，常采用“滞后时间”τ0来定量描述滞后时间的长短，即表示93%的超静孔隙水消散所需的时间（Riley，1969）。计算公式如下：
-
式中，Ss为黏性土层储水率（m-1），b0为土层初始厚度（m），Kv为黏性土层垂向渗透系数（m/d）。根据前人研究结果，可以利用现场监测数据绘制应力（水位）与应变（位移）的关系曲线，估算含水层系统的弹性和非弹性储水率。在此基础上，计算黏性土层释水滞后时间，具体计算方法可参考以下文献（Riley，1969; Helm，1975; Burbey，2001; 叶淑君等，2005; Liu Yi et al.，2008a，2008b）。文中利用平原区7个地面沉降监测站内分层标和地下水位监测数据，结合南水进京前后土层变形特征，分阶段计算了典型监测层位土层弹性和非弹性储水率，并对黏性土层释水滞后时间进行估算。
-
2.2.3 GIS、数理统计和应力-应变分析
-
文中利用GIS空间分析和数理统计方法，分析南水进京前后，不同含水层组地下水位同期对比情况，以及时序地下水位和地面沉降动态变化，进而查明北京平原区地下水位和地面沉降时空演变规律。采用应力-应变分析方法，研究了南水进京前后，不同水位变化模式下含水砂层和第二、第三压缩层组的变形特征。
-
3 结果与讨论
-
3.1 地下水流场和地下水位变化特征
-
文中利用2015年和2020年平原区分层地下水位监测数据，绘制了2015年和2020年不同含水层组地下水位同期对比图（图6）。可以发现，南水进京以后，北京平原区不同含水层组地下水位均出现不同程度的上升。特别是平原区北部，部分地区地下水位上升显著。2015~2020年，潜水含水层地下水位出现上升的地区主要分布在平原区北部的密怀顺地区、平谷城区、昌平部分地区以及平原区西部和南部的丰台、大兴等地（图6a）。第一承压含水层组地下水位在平原区北部的顺义、昌平、平谷以及朝阳和大兴部分地区同比出现上升（图6b）。第二承压含水层组地下水位在平原区北部的顺义、昌平、平谷、海淀以及朝阳大部分地区、通州和大兴北部等地均出现不同程度的上升（图6c）。第三承压含水层组地下水位除了通州东南部和大兴南部等地同比下降外，其余地区地下水位同比均上升（图6d）。
-
根据平原区北部的1号至6号地下水位监测井资料显示，平原区北部地下水位从2015年开始出现上升。截至2020年底，井1、井2、井3地下水位分别上升了21.62 m、21.72 m和17.24 m（图7a），井4、井5、井6地下水位分别上升了2.11 m、12.59 m和23.15 m（图7b）。这主要是因为南水进京后，北京市及时对平原区北部的密怀顺应急水源地进行了地下水减采和热备涵养，开采量减少了2/3（王树芳等，2019）。同时，利用南水通过潮白河河道向密怀顺应急水源地进行生态补水。截至2020年底，该地区已累计补水超过5×108m3，由此造成平原区北部地下水位出现较大幅度上升（雷坤超等，2020）。平原区东部地下水位变化趋势如图7c所示。2017年前，各含水层组地下水位呈持续下降的趋势，年均降幅0.13~1.82 m。2017年后，各含水层组地下水位由降转升，年均升幅0.45~1.87 m。然而，在平原区东南部和南部的通州和大兴部分地区，地下水位仍呈持续下降的趋势。根据通州区东南部的7号至10号地下水位监测井资料显示，2010~2020年，中深层和深层承压水水头仍在持续下降，年均降幅1.18~1.75 m，最大累计降幅达到19.28 m（图7d）。由此说明这些地区的地下水仍在超量开采。
-
图6 北京平原区2015年和2020年不同含水层组地下水位同期对比
-
Fig.6 Comparison of groundwater level of different aquifer groups in 2015 and 2020
-
（a）—潜水含水层；（b）—第一承压含水层；（c）—第二承压含水层；（d）—第三承压含水层
-
（a）—unconfined aquifer; （b）—first confined aquifer; (c) —second confined aquifer; (d) —third confined aquifer
-
3.2 地面沉降变化特征
-
文中采用PS-InSAR技术获取了2013~2020年北京平原区地面沉降时间序列信息（图8）。首先，从地面沉降分布来看，北京地面沉降严重区（沉降速率大于50 mm/a）主要集中在平原区东部的朝阳和通州部分地区，以及平原区北部的昌平、顺义和海淀北部等地。几个主要沉降中心连片发展，最大沉降速率多年保持在100 mm/a以上。其次，从地面沉降变化趋势来看，自2014年底南水进京后，北京市地面沉降总体呈减缓的趋势。其中，平原区北部地面沉降于2016年开始逐渐减缓，平原区东部地面沉降于2018年开始明显减缓。2013~2020年，平原区区域沉降速率由2013年的21.7 mm/a减小至2020年的10.85 mm/a，减小了10.85 mm/a。沉降严重区面积由2013年的572 km2减小至2020年的45 km2，减小了527 km2。地面沉降体积由2013年的13210×104 m3减小至2020年的5338×104m3，减小了8120×104m3。
-
对平原区7个地面沉降站内分层标监测结果分析发现：2013~2020年，平原区北部的天竺、望京、八仙庄、平各庄和东部的王四营、张家湾6个监测站，地面沉降量总体呈逐渐减小的趋势。其中，八仙庄站年沉降量减小最多，减小了67.0 mm/a；其次是王四营站，年沉降量减小了63.0 mm/a。平原区南部的榆垡站，2020年较2013年，年沉降量反而增加了5.0 mm/a（图9）。从地面沉降垂向分布来看，目前北京平原区地面沉降主要集中在第二和第二压缩层组。2013~2020年，除平原区北部的平各庄和南部的榆垡站外，其余各站第三压缩层组沉降占比逐渐增大，且部分监测站第三压缩层组沉降占比超过90%。平各庄站第三压缩层组沉降占比呈现出先增大后减小的趋势。榆垡站第二和第三压缩层组沉降占比则较为平稳，维持在20%和62%左右（表2）。
-
图7 典型地区不同含水层组地下水位变化曲线
-
Fig.7 Groundwater level variation curves of different aquifer groups in typical areas
-
（a）—平原区北部1、2、3号井；（b）—平原区北部潮白河河道附近4、5、6号井；（c）—平原区东部王四营地面沉降站D1-1至D1-5号井；（d）—平原区东南部7、8、9、10号井
-
（a）—wells 1, 2 and 3 in the northern plain; (b) —wells 4, 5 and 6 near Chaobai River in the northern plain; （c）—well D1-1 to D1-5 of Wangsiying land subsidence station in the eastern of plain; （d）—wells 7, 8, 9 and 10 in the southeast of plain
-
3.3 土层变形特征
-
文中利用7个地面沉降监测站分层标和地下水位监测数据，分析南水进京前后不同岩性、不同深度土层的变形特征。选取有代表性的土层变形进行分析，将特征相似的土层归纳到一起。由于第一压缩层组对应的潜水含水层水位容易受到外界因素的影响，导致第一压缩层组具有较为复杂的变形特征。同时，地面沉降主要集中在第二和第三压缩层组，所以下面主要分析第二、第三压缩层组和含水砂层的变形特征。图10~16中，累计变形的负值代表土层压缩，正值代表回弹。
-
3.3.1 含水砂层
-
天竺站分层标F3-8监测层位49~65 m，总厚度16 m。岩性主要为细砂和中粗砂，其顶部覆盖约2 m厚的黏土层。从图10a可以看出：2006~2016年，地下水位总体呈下降的趋势，年均水位由2006年的-4.37 m下降至2016年的-8.73 m。分层标F3-8累计压缩了5.19 mm。从2017年开始，地下水位逐渐上升。至2020年底，年均水位上升至-3.00 m，分层标F3-8回弹了1.94 mm。2006~2016年，随着地下水位往复升降，含水砂层表现出几乎同步的压缩-回弹。2017年后，地下水位由降转升过程中，含水砂层变形略有滞后。从图10b可以看出：2006~2016年，含水砂层在缓慢的循环加卸荷的作用下，加载与卸载曲线非常接近，表现出与水位变化基本同步的弹性变形。2017年后，地下水位开始上升，但砂层并没有立刻回弹，存在变形滞后现象。这可能与覆盖在砂层顶部的黏土层发生蠕变有关。
-
图8 北京平原区2013~2020年（a~h）地面沉降空间分布特征
-
Fig.8 Spatial distribution characteristics of land subsidence in Beijing Plain from 2013 to 2020 （a~h）
-
3.3.2 第二压缩层组
-
（1）平原区北部。天竺站分层标F3-7监测层位65~82 m，总厚度17 m。岩性以黏性土层为主，间夹1 m厚的粉细砂层。八仙庄站分层标F4-6监测层位108~144 m，总厚度36 m。岩性以粉质黏土为主，间夹约5.5 m厚的细砂层。从图11a、c可以看出：2006~2016年，D3-4和D4-4地下水位总体呈下降的趋势，2017年开始，地下水位出现回升。分层标F3-7和F4-6在2017年之前，随着地下水位下降，土层持续快速压缩。2017年后，随着地下水位回升，沉降曲线由陡变缓，土层压缩速率有所减缓。从图11b、d可以看出：2006~2016年，分层标F3-7塑性变形量较大，存在滞后现象。这一方面与土层超静孔隙水压力消散滞后于含水层水位变化有关，另外还可能和蠕变有关。2017~2020年，土层在水位上升阶段仍在持续压缩，但压缩速率有所减小，土层包括塑性变形和蠕变变形。同时，在部分年度出现回滞环，说明还存在弹性变形。然而分层标F4-6随着地下水位的往复升降，土层始终持续压缩，存在滞后效应，说明土层不仅存在残余变形量较大的塑性变形，还可能存在随时间发展的蠕变变形。
-
图9 北京平原区2013~2020年7个监测站年沉降量变化
-
Fig.9 Annual subsidence of 7 monitoring stations in Beijing Plain from 2013 to 2020
-
图10 天竺站分层标F3-8处含水砂层累计变形与水位关系
-
Fig.10 Relationship between cumulative deformation of sand layer with groundwater level at F3-8 of Tianzhu station
-
（a）—水位与累计变形量历时曲线；（b）—水位与累计变形量关系曲线
-
（a）—time series curve of relationship between groundwater level and cumulative deformation； （b）—relationship curve between groundwater level and cumulative deformation
-
注：表中占比负值表示回弹。
-
图11 天竺站分层标F3-7和八仙庄站分层标F4-6处土层累计变形与水位关系
-
Fig.11 Relationship between cumulative deformation of soil layer with groundwater level at F3-7 of Tianzhu station and F4-6 of Baxianzhuang station
-
（a）—分层标F3-7处水位与累计变形量历时曲线；（b）—分层标F3-7处水位与累计变形量关系曲线；（c）—分层标F4-6处水位与累计变形量历时曲线；（d）—分层标F4-6处土层累计变形与水位关系
-
（a）—time series curve of relationship between groundwater level and cumulative deformation at F3-7； （b）—relationship curve between groundwater level and cumulative deformation of F3-7; （c）—time series curve of relationship between groundwater level and cumulative deformation at F4-6； （d）—relationship curve between groundwater level and cumulative deformation of F4-6
-
（2）平原区东部。王四营站分层标F1-3监测层位66~94 m，总厚度28 m。由黏土层和砂层组成，层次较多，其中砂层约占该段地层总厚度的74%，黏土层占比为26%。从图12a可以看出：2006~2016年，D1-3地下水位总体呈下降的趋势，2017年开始，地下水位出现回升。2017年之前，地下水位年度内呈周期性往复升降，在每个周期内，地下水位下降的幅度多大于上升幅度，土层经受的有效应力增加，土层持续快速压缩。2017年后，地下水位回升明显，土层压缩速率有所减缓，并在2019年和2020年出现小幅回弹。从图12b可以看出：虽然该监测层位以砂层为主，但2006~2016年随着地下水位下降，土层甚至没有出现回弹而仍在持续快速压缩，塑性变形量较大，并存在滞后现象。这一方面与土层超静孔隙水压力消散明显滞后于含水层水位变化有关，另外还可能和蠕变有关。2017年之后，当地下水位回升时，土层塑性变形量很小，存在变形滞后。2019年和2020年，土层又出现弹性变形。
-
（3）平原区南部。榆垡站分层标F7-4监测层位82~116 m，总厚度34 m。岩性以粉质黏土为主，间夹3.8 m厚的含砾中砂层。从图13a可以看出：2009~2020年，D7-3地下水位总体呈下降的趋势，土层经受的有效应力不断增加，土层持续快速压缩。从图13b可以看出：2009~2020年，随着地下水位的下降，土层持续快速压缩，以塑性变形为主。2018年以后土层出现小幅回弹，存在滞后变形。这一方面与土层超静孔隙水压力消散明显滞后于含水层水位变化有关，另外还可能和蠕变有关。
-
图12 王四营站分层标F1-3处土层累计变形与水位关系
-
Fig.12 Relationship between cumulative deformation of soil layer with groundwater level at F1-3 of Wangsiying station
-
（a）—分层标F1-3处水位与累计变形量历时曲线；（b）—分层标F1-3处水位与累计变形量关系曲线
-
（a）—time series curve of relationship between groundwater level and cumulative deformation at F1-3； （b）—relationship curve between groundwater level and cumulative deformation of F1-3
-
3.3.3 第三压缩层组
-
（1）平原区北部。平各庄站分层标F5-3监测层位209~234 m，总厚度25 m。岩性以黏性土为主，间夹约2.5 m厚粉细砂层。望京站分层标F2-1监测层位133 m以下。主要为黏性土与砂砾石互层结构，层次较多，并以黏性土层为主。从图14a可以看出：2009~2016年，D5-2地下水位总体呈下降的趋势，2017年开始，地下水位出现回升。分层标F5-3在2017年之前，随着地下水位下降，土层持续快速压缩。2017年后，随着地下水位回升，沉降曲线由陡变缓，土层压缩速率有所减缓，并在部分年度出现小幅回弹。从图14c可以看出：2009~2014年，D2-2地下水位呈下降趋势，2015年开始，地下水位出现回升。但分层标F2-1始终持续快速压缩，压缩速率未见减缓。从图14b、d可以看出：2006~2016年，分层标F5-3主要呈塑性变形，存在滞后现象。2017~2020年，土层在水位上升阶段仍在持续压缩，但压缩速率有所减小，土层包括塑性变形和蠕变变形。同时，在2019年和2020年出现回滞环，说明还存在弹性变形。然而分层标F2-1随着地下水位的往复升降，土层始终持续压缩，存在滞后效应，土层不仅存在残余变形量较大的塑性变形，还可能存在随时间发展的蠕变变形。
-
图13 榆垡站分层标F7-4处土层累计变形与水位关系
-
Fig.13 Relationship between cumulative deformation of soil layer with groundwater level at F7-4 of Yufa station
-
（a）—分层标F7-4处水位与累计变形量历时曲线；（b）—分层标F7-4处水位与累计变形量关系曲线
-
（a）—time series curve of relationship between groundwater level and cumulative deformation at F7-4； （b）—relationship curve between groundwater level and cumulative deformation of F7-4
-
（2）平原区东部。张家湾站分层标F6-3监测层位126~193 m，总厚度67 m。主要为黏性土层与砂层互层，其中黏性土层约占该段地层总厚度的52%，砂层占比为48%。王四营站分层标F1-2监测层位94~148 m，总厚度54 m。为黏性土层与砂层互层，其中黏性土层约占该段地层总厚度的77%。从图15a和图15c可以看出：2017年之前，D6-2和D1-2地下水位总体呈下降的趋势，土层持续压缩。2017年后，地下水位逐步回升，土层压缩速率有所减缓。其中，分层标F6-3在2019年和2020年出现小幅回弹。然而分层标F1-2则始终处于压缩状态，仅速率有所减小。从图15b、d可以看出：2006~2016年，分层标F6-3主要呈塑性变形，存在滞后现象。2017~2020年，土层在水位上升阶段仍在持续压缩，但压缩速率有所减小，土层包括塑性变形和蠕变变形。同时，在2019年和2020年出现回滞环，说明还存在弹性变形。然而分层标F1-2随着地下水位的往复升降，土层始终持续压缩，存在滞后效应，土层不仅存在残余变形量较大的塑性变形，还存在随时间发展的蠕变变形。
-
（3）平原区南部。榆垡站分层标F7-1监测层位205 m以下，主要为砂卵砾石和粉质黏土互层，并以砂卵砾石层为主。从图16a可以看出：2009~2020年，D7-1地下水位呈下降的趋势，土层始终持续快速压缩。从图16b可以看出：2009~2020年，随着地下水位的下降，土层塑性变形量较大，并存在滞后现象。这一方面与土层超静孔隙水压力消散明显滞后于含水层水位变化有关，另外还可能和蠕变有关。
-
以上是平原区北部、东部和南部地区，不同水位变化模式下不同深度及岩性土层的典型变形曲线。可以看出，自2014年底南水进京后，平原区北部和东部，大部分地区的中深层和深层承压水水头在2017年前持续下降，2017年后逐渐回升。而平原区南部，部分地区中深层和深层承压水水头始终呈下降的趋势。因此，在综合考虑不同地区、不同地下水位变化模式下，可将平原区第二、第三压缩层组和含水砂层的变形特征进行归纳总结（表3）。
-
3.4 弹性和非弹性储水率变化
-
根据Galloway et al.（2011）研究结果，Ssk的值约为5×10-4是典型的骨架弹性储水系数，而Ssk的值约为5×10-3很可能就是骨架非弹性储水系数。Hoffmann et al.（2003）指出，骨架非弹性储水系数通常是弹性储水系数的数十至数百倍。当然，在计算储水率时，需要将储水系数除以土层的初始厚度。文中根据上述不同岩性、不同深度土层变形与水位关系，计算出不同阶段含水砂层与压缩层组的弹性和非弹性储水率。结果的平均值见表4所示。需要注意的是，在地下水位持续下降阶段，分层标F5-3监测层位的非弹性储水率为8.0×10-5，量值较小，这可能是因为分层标F5-3监测层位在209~234 m，土层呈微超固结状态，压缩性较低的缘故。而分层标F6-3监测层位的非弹性储水率量值也较小，为4.8×10-5，这主要是该监测层位含水砂层占该段地层总厚度的48%，可能会对土层的非弹性变形产生一定影响。
-
对于土层变形的不同阶段，储水率的变化并不是恒定的，其值取决于有效应力与先期固结压力的关系。例如分层标F3-7、F7-4和F5-3，随着地下水位的下降，土层持续压缩。在地下水位由持续下降转为持续上升的过程中，土层压缩速率减缓，甚至出现弹性变形，但由于受到前期持续压缩的影响，土层孔隙比不断减小，非弹性储水率也呈现出逐渐减小的趋势，平均减小了59%左右，由此导致土层的释水和压缩能力也逐渐减弱。这与前人利用室内试验得出的结论相一致（李兆峰等，2017）。同时，像分层标F1-3在地下水位持续下降阶段，土层表现为塑性变形，当地下水位持续上升时，土层由塑性变形转变为弹性变形。储水率也由非弹性转变为弹性。其余分层标处的储水率，在地下水位由降转升的过程中，也都表现出类似的变化特征。分层标F3-8监测层位主要为细砂和中粗砂，土层几乎呈弹性变形特征。通过上述分析表明，南水进京以后，不同岩性和深度的土层变形发生了明显变化，相应的水文地质参数也发生了变化。因此，在模拟南水进京后地下水位持续上升条件下的地面沉降时，需要根据土层变形特征选取合适的土层变形本构关系和相应的水文地质参数。
-
图14 平各庄站分层标F5-3和望京站分层标F2-1处土层累计变形与水位关系
-
Fig.14 Relationship between cumulative deformation of soil layer with groundwater level at F5-3 of Pinggezhuang station and F2-1 of Wangjing station
-
（a）—分层标F5-3处水位与累计变形量历时曲线；（b）—分层标F5-3处水位与累计变形量关系曲线；（c）—分层标F2-1处水位与累计变形量历时曲线；（d）—分层标F2-1处土层累计变形与水位关系
-
（a）—time series curve of relationship between groundwater level and cumulative deformation at F5-3； （b）—relationship curve between groundwater level and cumulative deformation of F5-3; （c）—time series curve of relationship between groundwater level and cumulative deformation at F2-1； （d）—relationship curve between groundwater level and cumulative deformation of F2-1
-
3.5 土层残余变形和变形滞后
-
从表4可以看出，第二和第三压缩层组的非弹性储水率比弹性储水率大1~2个数量级，说明当地下水位恢复后，土层会存在较大的残余变形，会出现持续压缩而不回弹的现象。例如分层标F1-3监测层位总厚度28 m，利用表4中的参数，当地下水位下降5 m，不考虑土层的变形滞后，则土层将压缩43.4 mm；若地下水位上升5 m，则土层仅回弹 1.2 mm，仅占压缩量的2.8%。因此，可以利用非弹性储水率大于弹性储水率来解释土层存在较大残余变形量，但无法解释土层的变形滞后现象，特别是黏性土层变形明显滞后于相邻含水层水位变化（罗跃等，2015）。例如分层标F4-6、F5-3和F1-2，虽然其对应的含水层水位自2017年后持续上升，但土层仍在压缩，仅压缩速率有所减缓，土层变形存在明显的滞后性。
-
图15 张家湾站分层标F6-3和王四营站分层标F1-2处土层累计变形与水位关系
-
Fig.15 Relationship between cumulative deformation of soil layer with groundwater level at F6-3 of Zhangjiawan station and F1-2 of Wangsiying station
-
（a）—分层标F6-3处水位与累计变形量历时曲线；（b）—分层标F6-3处水位与累计变形量关系曲线；（c）—分层标F1-2处水位与累计变形量历时曲线；（d）—分层标F1-2处土层累计变形与水位关系）
-
（a）—time series curve of relationship between groundwater level and cumulative deformation at F6-3; （b）—relationship curve between groundwater level and cumulative deformation of F6-3; （c）—time series curve of relationship between groundwater level and cumulative deformation at F1-2; （d）—relationship curve between groundwater level and cumulative deformation of F1-2
-
图16 榆垡站分层标F7-1处土层累计变形与水位关系
-
Fig.16 Relationship between cumulative deformation of soil layer with groundwater level at F7-1 of Yufa station
-
（a）—分层标F7-1处水位与累计变形量历时曲线；（b）—分层标F7-1处水位与累计变形量关系曲线
-
（a）—time series curve of relationship between groundwater level and cumulative deformation at F7-1； （b）—relationship curve between groundwater level and cumulative deformation of F7-1
-
注：“—”表示未能计算。
-
黏性土层的滞后性主要体现在两个方面：① 黏性土层释水滞后。在承压含水层水头下降初期，黏性土层中靠近开采层位的部分孔隙水压力率先减小，远离开采层位的部分孔隙水压力变化滞后，因此黏性土层出现滞后释水现象。② 黏性土层变形滞后。黏性土层中靠近开采层位的土体率先压缩，远离开采层位的土体滞后压缩，有效应力的变化在黏性土层中出现滞后，导致了黏性土层变形出现滞后。黏性土层表现出的滞后性，主要是由于其弱渗透性造成的，黏性土层的渗透系数比含水层小3~5个数量级，造成黏性土层释水滞后于相邻含水层水位变化，从而导致其变形出现明显的滞后。需要注意的是，黏性土层所表现出的这种变形滞后与黏性土蠕变无关，仅是土层边界处水位引起的主固结需要一定时间来完成。主固结过程可能持续几个月、几十甚至数百年（叶淑君等，2005）。如分层标F3-7监测层位中黏性土层厚度b0=17 m，渗透系数Kv=5.56×10-5 m/d，非弹性储水率Ss=4.3×10-4m-1，根据公式（2）计算可得含水层水头下降时，黏性土层的变形滞后时间为559天。也就是说与F3-7对应的D3-4水位虽然在2017年后出现回升，但F3-7监测的黏性土层仍在响应D3-4在559天前的水位下降条件。因此部分土层在含水层水位普遍上升的条件下，仍在持续压缩。同时，随着黏性土层中超静孔隙水压力不断消散，以及土层不断压密，黏性土层的压缩速率将逐渐减缓。另外，黏性土层变形表现出的滞后性，除了土层渗透固结成因以外，土体蠕变是另外一个重要原因。如天津和沧州地区中、晚更新世土层蠕变可占压缩变形的14.3%~28.3%（郭海朋等，2017）。因此在开展地面沉降中长期预报时，忽略土体蠕变变形，会引起较大误差，应根据含水层组地下水位、孔隙水压力和先期固结压力，分析土层所处的变形阶段，考虑是否加入蠕变变形。
-
4 结论
-
本文利用多种手段及监测数据，分析了南水进京前后，北京平原区地下水位、地面沉降、土层变形和水文地质参数的变化，探讨了黏性土层产生较大残余变形和变形滞后的原因。研究结果对进一步掌握南水进京后，北京平原区地下水和地面沉降状况，确定土层变形本构关系，建立合适的地下水与地面沉降耦合模型具有重要意义。
-
（1）2015~2020年，北京平原区大部分地区第一至第四含水层组地下水位逐渐上升，地面沉降呈减缓的趋势。然而，在平原区东南部和南部的通州和大兴部分地区，地下水位仍持续下降，地面沉降继续发展，但沉降速率有所减缓。
-
（2）南水进京以来，除榆垡站外，其余各站沉降速率均减小。第二和第三压缩层组是平原区地面沉降主要贡献层。除平各庄和榆垡站外，其余各站第三压缩层组沉降占比逐渐增大，沉降主控层有向深部转移的规律。
-
（3）平原区北部和东部，第二和第三压缩层组对应的地下水位由降转升，在水位下降阶段，土层呈塑性和蠕变变形；水位上升阶段，土层以塑性变形为主，部分时间出现弹性变形，具有黏弹塑性变形特征。然而平原区南部，地下水位持续下降，土层变形始终呈塑性和蠕变变形。含水砂层则主要呈弹性变形特征。
-
（4）土层变形的不同阶段，弹性和非弹性储水率并不是恒定的，其值取决于有效应力与先期固结压力的关系。非弹性储水率比弹性储水率大1~2个数量级，且随着地下水位下降，储水率呈减小的趋势。
-
（5）黏性土层存在较大残余变形和变形滞后的原因，可以从两方面解释：① 非弹性储水率大于弹性储水率。② 黏性土层的弱渗透性，造成其释水滞后于相邻含水层水位变化，进而导致土层变形出现明显滞后现象。
-
参考文献
-
Abidin H Z, Andreas H, Djaja R, Darmawan D, Gmal M. 2008. Land subsidence characteristics of Jakarta between 1997 and 2005, as estimated using GPS surveys. GPS Solutions, 12(1): 23~32.
-
Abidin H Z, Andreas H, Gumilar I, Fukuda Y, Pohan Y E, Deguchi T. 2011. Land subsidence of Jakarta (Indonesia) and its relation with urban development. Natural Hazards, 59(3): 1753~1771.
-
Burbey T. 2001. Stress-strain analyses for aquifer-system characterization. Ground Water, 39(1): 128~136.
-
Cai Xiangmin, Luan Yingbo, Guo Gaoxuan, Liang Yanan. 2009. 3D Quaternary geological structure of Beijing Plain. Geology in China, 36(5): 1021~1029.
-
Cao A, Esteban M, Mino T. 2020. Adapting wastewater treatment plants to sea level rise: Learning from land subsidence in Tohoku, Japan. Natural Hazards, 103(7): 885~902.
-
Chaussard E, Wdowinski S, Cabral-Cano E, Amelung F. 2014. Land subsidence in central Mexico detected by ALOS InSAR time-series. Remote Sensing of Environment, 140: 94~106.
-
Chen Beibei, Gong Huili, Lei Kunchao, Li Xiaojuan, Zhou Chaofan, Gao Mingliang, Guan Hongliang, Lv Wei. 2019. Land subsidence lagging quantification in the main exploration aquifer layers in Beijing Plain, China. International Journal of Applied Earth Observation and Geoinformation, 75: 54~67.
-
Chen Beibei, Gong Huili, Chen Yun, Li Xiaojuan, Zhou Chaofan, Lei Kunchao, Zhu Lin, Duan Li, Zhao Xiaoxiao. 2020. Land subsidence and its relation with groundwater aquifers in Beijing Plain of China, Science of the Total Environment, 735: 139111.
-
Di Shengtong, Jia Chao, Zhang Shaopeng, Ding Pengpeng, Shao Ming, Zhang Yongwei. 2020. Regional characteristics and evolutionary trend prediction of land subsidence caused by groundwater over exploitation in North Shandong of the North China Plain. Acta Geologica Sinica, 94(5): 1638~1654 (in Chinese with English abstract).
-
Dong Jie, Lai Shangjing, Wang Nan, Wang Yan, Zhang Lu, Liao Mingsheng. 2021. Multi-scale deformation monitoring with Sentinel-1 InSAR analyses along the Middle Route of the South-North Water Diversion Project in China. International Journal of Applied Earth Observation and Geoinformation, 100: 102324.
-
Ferretti A, Prati C. 2000. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE Transactions on Geoscience & Remote Sensing, 38(5): 2202~2212.
-
Ferretti A, Prati C, Rocca F. 2001. Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience & Remote Sensing, 39(1): 8~20.
-
Galloway D L, Hudnut K W, Ingebritsen S E, Phillips S P, Peltzer G, Rogez F, Rosen P A. 1998. Detection of aquifer system compaction and land subsidence using interferometric synthetic aperture radar, Antelope Valley, Mojave Desert, California. Water Resources Research, 34(10): 2573~2585.
-
Galloway D L, Burbey T J. 2011. Review: Regional land subsidence accompanying groundwater extraction, Hydrogeology Journal, 19(8): 1459~1486.
-
Gong Huili, Pan Yun, Zheng Longqun, Li Xiaojuan, Zhu Lin, Zhang Chong, Huang Zhiyong, Li Zhiping, Wang Haigang, Zhou Chaofan. 2018. Long-term groundwater storage changes and land subsidence development in the North China Plain (1971-2015). Hydrogeology Journal, 26(5): 1417~1427.
-
Guo Haipeng, Bai Jinbin, Zhang Youquan, Wang Liya, Shi Jusong, Li Wenpeng, Zhang Zuochen, Wang Yunlong, Zhu Juyan, Wang Haigang. 2017. The evolution characteristics and mechanism of the land subsidence in typical areas of the North China Plain. Geology in China, 44(6): 1115~1127 (in Chinese with English abstract).
-
Guo Haipeng, Li Wenpeng, Wang Liya, Chen Ye, Zang Xisheng, Wang Yunlong, Zhu Juyan, Bian Yueyue. 2021. Present situation and research prospects of the land subsidence driven by groundwater levels in the North China Plain. Hydrogeology & Engineering Geology, 48(3): 162~171 (in Chinese with English abstract).
-
Helm D C. 1975. One-dimensional simulation of aquifer system compaction near Pixley, California: 1. Constant parameters. Water Resources Research, 11(3): 465~478.
-
Herrera-García G, Ezquerro P, Tomás R, Béjar-Pizarro M, López-Vinielles J, Rossi M, Mateos R M, Carreón-Freyre D, Lambert J, Teatini P, Cabral-Cano E, Erkens G, Galloway D, Hung W C, Kakar N, Sneed M, Tosi L, Wang H M, Ye S J. 2021. Mapping the global threat of land subsidence. Science, 371(6524): 34~36.
-
Hoffmann J, Galloway D L, Zebker H A. 2003. Inverse modeling of interbed storage parameters using land subsidence observations, Antelope Valley, California. Water Resources Research, 39: 1031.
-
Hu Beibei, Zhou Jun, Wang Jun, Chen Zhenlou, Wang Dongqi, Xu Shiyuan. 2009. Risk assessment of land subsidence at Tianjin coastal area in China. Environmental Earth Sciences, 59(2): 269~276.
-
Hu R L, Yue Z Q, Wang L C, Wang S J. 2004. Review on current status and challenging issues of land subsidence in China. Engineering Geology, 76(1-2): 65~77.
-
Jeanne P, Farr T G, Rutqvist J, Vasco D. 2019. Role of agricultural activity on land subsidence in the San Joaquin Valley, California. Journal of Hydrology, 569: 462~469.
-
Jia Sanman, Ye Chao, Luo Yong, et al. 2018. Beijing Land Subsidence. Beijing: Geological Publishing House (in Chinese with English abstract).
-
Lei Kunchao, Luo Yong, Chen Beibei, Guo Gaoxuan, Zhou Yi. 2016. Distribution characteristics and influence factors of land subsidence in Beijing area. Geology in China, 43(6): 2216~2225 (in Chinese with English abstract).
-
Lei Kunchao, Luo Yong, Liu He, Wang Xinhui, Cui Wenjun, Tian Miaozhuang, Qi Minghuan. 2020. Land subsidence monitoring report of Beijing in 2019. Beijing Institute of Hydrogeology and Engineering Geology.
-
Lei Kunchao, Ma Fengshan, Luo Yong, Chen Beibei, Cui Wenjun, Tian Fang, Sha Te. 2022a. Main subsidence layers and deformation characteristics in Beijing Plain at present. Journal of Engineering Geology, 30(2): 442~458 (in Chinese with English abstract).
-
Lei Kunchao, Ma Fengshan, Chen Beibei, Luo Yong, Cui Wenjun, Liu He, Tian Fang. 2022b. Three-Dimensional (3D) surface displacement in Beijing Plain based on time series InSAR and GPS technologies. Journal of Engineering Geology, 30(2): 417~431 (in Chinese with English abstract).
-
Li Wenpeng, Wang Longfeng, Yang Huifeng, Zheng Yuejun, Cao Wengeng, Liu Ke. 2020. The groundwater over exploitation status and countermeasure suggestions of the North China Plain. China Water Resources, 13: 26~30 (in Chinese with English abstract).
-
Li Xue, Ye Siyuan, Wei Aihua, Zhou Pengpeng, Wang Liheng. 2017. Modelling the response of shallow groundwater levels to combined climate and water-diversion scenarios in Beijing-Tianjin-Hebei Plain, China. Hydrogeology Journal, 25(6): 1733~1744.
-
Li Zhaofeng, Zhou Zhifang, Li Mingyuan, Zhou Cuiying. 2017. Variation of hydraulic parameters of aquitard during water release. Journal of Hehai University (Natural Sciences), 45(4): 340~344(in Chinese with English abstract).
-
Liu Yi, Helm D C. 2008a. Inverse procedure for calibrating parameters that control land subsidence caused by subsurface fluid withdrawal: 1. Methods. Water Resources Research, 44(7): W07423.
-
Liu Yi, Helm D C. 2008b. Inverse procedure for calibrating parameters that control land subsidence caused by subsurface fluid withdrawal: 2. Field application. Water Resources Research, 44(7): W07424.
-
Luo Yue, Ye Shujun, Wu Jichun, Jiao Xun, Wang Hanmei. 2015. Characterization of land subsidence during recovery of groundwater levels in Shanghai. Geological Journal of China Universities, 21(2): 243~254(in Chinese with English abstract).
-
Mahmoudpour M, Khamehchiyan M, Nikudel M R, Ghassemi M R. 2016. Numerical simulation and prediction of regional land subsidence caused by groundwater exploitation in the southwest plain of Tehran, Iran. Engineering Geology, 201: 6~28.
-
Motagh M, Walter T R, Sharifi M A, Fielding E, Schenk A, Anderssohn J, Zschau J. 2008. Land subsidence in Iran caused by widespread water reservoir overexploitation. Geophysical Research Letters, 35(16): 1029~1039.
-
Pavelko M T, Hoffmann J, Damar N A. 2006. Interferograms Showing Land Subsidence and Uplift in Las Vegas Valley, Nevada, 1992-99. U. S. Geological Survey, Reston, Virginia.
-
Phien-Wej N, Giao P H, Nutalaya P. 2006. Land subsidence in Bangkok, Thailand. Engineering Geology, 82(4): 187~201.
-
Riley F S. 1969. Analysis of borehole extensometer data from central California. Land Subsidence (Vol. 2), International Association Sciences Hydrological Publication, 89: 423~432.
-
Strozzi T, Wegmuller U. 1999. Land subsidence in Mexico City mapped by ERS differential SAR interferometry. Geoscience and Remote Sensing Symposium, IGARSS ‘99 Proceedings. IEEE 1999 International, 4: 1940~1942.
-
Teatini P, Tosi L, Strozzi T, Carbognin L, Cecconi G, Rosselli R, Libardo S. 2010. Resolving land subsidence within the Venice Lagoon by persistent scatterer SAR interferometry. Physics and Chemistry of the Earth Parts A/B/C, 40: 72~79.
-
Tosi L, Teatini P, Carbognin L, Frankenfield J. 2007. A new project to monitor land subsidence in the northern Venice coastland (Italy). Environmental Geology, 52(5): 889~898.
-
Wang Shufang, Li Jie, Liu Yuanzhang, Liu Jiurong, Wang Xu. 2019. Impact of South-to-North Water Diversion on groundwater recovery in Beijing. China Water Resources, 7: 26~30 (in Chinese with English abstract).
-
Xu Yeshuang, Ma Lei, Du Yanjun, Shen Shuilong. 2012. Analysis of urbanisation-induced land subsidence in Shanghai. Natural Hazards, 63(2): 1255~1267.
-
Xue Yuqun, Zhang Yun, Ye Shujun, Wu Jichun, Wei Zixin, Li Qinfen, Yu Jun. 2006. Research on the problems of land subsidence in China. Geological Journal of China Universities, 12(2): 153~160 (in Chinese with English abstract).
-
Yang Qing, Jiang Yuan, Lin Jian, Cui Wenjun. 2017. The environmental impacts on recharging groundwater from South-to-North Water Diversion Project. Urban Geology, 12(4): 30~34 (in Chinese with English abstract).
-
Yang Yong, Li Guomin, Dong Yanhui, Li Ming, Yang Jianqing, Zhou Dong, Yang Zhongshan, Zheng Fandong. 2012. Influence of South to North Water Transfer on groundwater dynamic change in Beijing Plain. Environmental Earth Sciences, 65(4): 1323~1331.
-
Ye Aizhong, Duan Qingyun, Chu Wei, Xu Jing, Mao Yuna. 2014. The impact of the South-North Water Transfer Project (CTP)'s central route on groundwater table in the Hai River basin, North China. Hydrological Processes, 28(23): 5755~5768.
-
Ye Shujun, Luo Yue, Wu Jichun, Yan Xuexin, Wang Hanmei, Jiao Xun, Teatini P. 2016. Three-dimensional numerical modeling of land subsidence in Shanghai, China. Hydrogeology Journal, 24(3): 695~709.
-
Ye Shujun, Xue Yuqun. 2005. Stress-strain analysis for storage coefficients and vertical hydraulic conductivities of aquitards in Shanghai area. Rock and Soil Mechanics, 26(2): 256~260.
-
Zhang Menglin, Hu Litang, Yao Lili, Yin Wenjie. 2018. Numerical studies on the influences of the South-to-North Water Transfer Project on groundwater level changes in the Beijing Plain, China. Hydrological Processes, 32(12): 1858~1873.
-
Zhang Suangcheng, Xu Qiang, Luo Yong, Lei Kunchao, Niu Yufen, Pang Xiaoguang. 2022. Temporal and spatial variation of land subsidence in Beijing from 2017 to 2020 interpreted by time series InSAR. Journal of Geodesy and Geodynamics, 42(1): 48~53 (in Chinese with English abstract).
-
Zhang Yuan, Li Guomin. 2014. Influence of South-to-North Water Diversion on major cones of depression in North China Plain. Environmental Earth Sciences, 71(9): 3845~3853.
-
Zhao Hongyi, Liu Wenguang, Sun Chunyuan, Jiao Zhongzhi, Li Minshi, Wang Sufen, Liu Cuizhu. 2020. Beijing Water Resources Bulletin. Beijing: Beijing Water Authority (in Chinese with English abstract).
-
Zhao Yong, Zhu Yongnan, Lin Zhaohui, Wang Jianhua, He Guohua, Li Haihong, Li Lei, Wang Hao, Jiang Shan, He Fan, Zhai Jiaqi, Wang Lizhen, Wang Qingming. 2017. Energy reduction effect of the South-to-North Water Diversion Project in China. Scientific Reports, 7(1): 1~9.
-
Zhu Lin, Gong Huili, Chen Yun, Wang Shufang, Ke Yinghai, Guo Gaoxuan, Li Xiaojuan, Chen Beibei, Wang Haigang, Teatini P. 2020. Effects of Water Diversion Project on groundwater system and land subsidence in Beijing, China. Engineering Geology, 276: 105763.
-
赵泓漪, 刘文光, 孙春媛, 焦忠志, 李民诗, 王素芬, 刘翠珠. 2020. 北京市水资源公报. 北京: 北京市水务局.
-
蔡向民, 栾英波, 郭高轩, 梁亚南. 2009. 北京平原第四系的三维结构. 中国地质, 36(5): 1021~1029.
-
狄胜同, 贾超, 张少鹏, 丁朋朋, 邵明, 张永伟．2020．华北平原鲁北地区地下水超采导致地面沉降区域特征及演化趋势预测. 地质学报, 94(5): 1638~1654.
-
郭海朋, 白晋斌, 张有全, 王丽亚, 石菊松, 李文鹏, 张作辰, 王云龙, 朱菊艳, 王海刚. 2017. 华北平原典型地段地面沉降演化特征与机理研究. 中国地质, 44(6): 1115~1127.
-
郭海朋, 李文鹏, 王丽亚, 陈晔, 臧西胜, 王云龙, 朱菊艳, 卞跃跃. 2021. 华北平原地下水位驱动下的地面沉降现状与研究展望. 水文地质工程地质, 48(3): 162~171.
-
贾三满, 叶超, 罗勇, 等. 2018. 北京地面沉降. 北京: 地质出版社.
-
雷坤超, 罗勇, 陈蓓蓓, 郭高轩, 周毅. 2016. 北京平原区地面沉降分布特征及影响因素. 中国地质, 43(6): 2216~2225.
-
雷坤超, 罗勇, 刘贺, 王新惠, 崔文君, 田苗壮, 齐鸣欢. 2020. 北京市地面沉降监测年度报告（2020年）. 北京市水文地质工程地质大队（北京市地质环境监测总站）.
-
雷坤超, 马凤山, 罗勇, 陈蓓蓓, 崔文君, 田芳, 沙特. 2022a. 北京平原区现阶段主要沉降层位与土层变形特征. 工程地质学报, 30(2): 442~458.
-
雷坤超, 马凤山, 陈蓓蓓, 罗勇, 崔文君, 刘贺, 田芳．2022b. 基于时序InSAR和GPS技术的北京平原区地表三维形变场特征. 工程地质学报, 30(2): 417~431.
-
李文鹏, 王龙凤, 杨会峰, 郑跃军, 曹文庚, 刘可. 2020. 华北平原地下水超采状况与治理对策建议. 中国水利, 13: 26~30.
-
李兆峰, 周志芳, 李明远, 周翠英. 2017. 弱透水层释水过程中水力参数响应规律. 河海大学学报（自然科学版）, 45(4): 340~344.
-
罗跃, 叶淑君, 吴吉春, 焦珣, 王寒梅. 2015. 上海市地下水位大幅抬升条件下土层变形特征分析. 高校地质学报, 21(2): 243~254.
-
王树芳, 李捷, 刘元章, 刘久荣, 王旭. 2019. 南水北调对北京地下水涵养的影响. 中国水利, 7: 26~30.
-
薛禹群, 张云, 叶淑君, 吴吉春, 魏子新, 李勤奋, 于军. 2006. 我国地面沉降若干问题研究. 高校地质学报, 12(2): 153~160.
-
杨庆, 姜媛, 林健, 崔文君. 2017. 南水北调水回灌对地下水环境的影响研究. 城市地质, 12(4): 30~34.
-
叶淑君, 薛禹群. 2005. 应用沉降和水位数据计算上海地区弱透水层的参数. 岩土力学, 26(2): 256~260.
-
张双成, 许强, 罗勇, 雷坤超, 牛玉芬, 庞校光. 2022. 时序InSAR解译2017~2020年北京地面沉降时空变化. 大地测量与地球动力学, 42(1): 48~53.
-
摘要
超量开采地下水引发的地面沉降已成为制约北京区域社会经济可持续发展的重要因素之一。2014年12月，南水北调中线工程正式通水，每年向北京输水超过10×108 m3，改变了北京供水格局，也为地下水压采、涵养及控制地面沉降创造了条件。本文利用多种监测数据，分析南水进京前后，北京平原区地下水和地面沉降的变化；研究不同水位变化模式下不同岩性及深度土层的变形特征；计算土层不同变形阶段的弹性和非弹性储水率；并对黏性土层产生较大残余变形和滞后变形的原因进行了探讨。结果表明：① 2015~2020年，平原区大部分地区第一至第四含水层组地下水位逐渐上升，地面沉降呈减缓的趋势。② 第二和第三压缩层组是沉降主要贡献层，除平各庄和榆垡站外，其余各站第三压缩层组沉降占比逐渐增大，沉降主控层有向深部转移的规律。③ 平原区北部和东部，第二和第三压缩层组对应的地下水位由降转升。在水位下降阶段，土层呈塑性和蠕变变形；水位上升阶段，土层以塑性变形为主，部分时间出现弹性变形，具有黏弹塑性。平原区南部，地下水位始终持续下降，土层变形始终呈塑性和蠕变变形。含水砂层则主要呈弹性变形。④ 土层变形的不同阶段，弹性和非弹性储水率并不是恒定的，随着地下水位下降，储水率呈减小的趋势。⑤ 黏性土层存在较大残余变形和变形滞后的原因，一是非弹性储水率大于弹性储水率，二是黏性土层的弱渗透性。
Abstract
Land subsidence caused by over-exploitation of groundwater has become one of the important factors restricting the sustainable development of the society and economy of Beijing. In December 2014, the middle route of the South-to-North Water Diversion Project was officially opened, and more than 1 billion m3 of water was delivered to Beijing every year. As a result, the pattern of water supply in Beijing has changed, and conditions have been created to reduce the exploitation of groundwater, conserve groundwater, and control land subsidence in Beijing. This paper uses a variety of monitoring data to analyze the changes in groundwater and land subsidence in the Beijing Plain before and after the South-to-North Water Diversion Project. The deformation characteristics of soil layers with different lithologies and depths under different groundwater level variation patterns are studied. The elastic and inelastic skeletal-specific storage rates of the soil layers at different deformation stages are calculated. The reasons for the large residual deformation and hysteresis deformation of the cohesive soil layer are discussed. The results show that: ① From 2015 to 2020, the groundwater level of the first to fourth aquifer groups in most areas of the plain gradually increased, and the land subsidence gradually slowed down. ② The second and third compression layer groups in the plain are the main contributing layers to the subsidence. Except for Pinggezhuang and Yufa station, the subsidence proportion of the third compression layer group in the other stations gradually increased. The main subsidence layer is gradually transferred to the deep formation. ③ In the north and east of the plain, the groundwater levels corresponding to the second and third compression layer groups changed from falling to rising. During the drop stage of the groundwater level, the soil layer exhibits plastic and creep deformation. During the rise stage of groundwater level, the soil layer exhibits plastic deformation and elastic deformation for part of the time. which is viscoelastic plastic. The soil layer exhibits the characteristics of viscoelastic plastic. In the southern part of the plain, the groundwater levels continued to decline, and the deformation of the soil layer was always plastic and creep deformation. The sand layer mainly exhibits elastic deformation. ④ For different stages of soil deformation, the change in elastic and inelastic skeletal-specific storage rates is not constant. As the groundwater level decreases, the soil skeletal specific storage rate shows a decreasing trend. ⑤ The reasons for the large residual deformation and deformation hysteresis of the cohesive soil layer mainly include two aspects. Firstly, the inelastic skeletal specific storage rate is greater than the elastic skeletal specific storage rate. Secondly, the weak permeability of the cohesive soil layer.
