地球系统科学是目前地球科学发展的一个前缘，人们往往更关注水圈、大气圈及冻土圈等变化速率较明显的现象，却较少关注人类活动对地球表层和固体地球的影响.本文将提出华北地下水长期过量开采对地壳应力状态的影响问题，希望引起我国地球科学者的关注.

随着经济、人口的增长，水资源需求不断增加，地下水过度开采已成为全球性焦点问题.特别是近50年来，华北地区地下水持续过量开采引起区域地下水漏斗广泛分布.自20世纪70年代以来，华北地区地下水漏斗加速扩展(Shah et al.，2003;Zhang and Li，2013)，其中北京、石家庄等城市为主要的地下水沉降中心.地下水位埋深大于10 m的区域占华北平原总面积的40%，且潜水含水层最大埋深达65 m.目前，华北平原深层地下水水位下降速率高达3 m·a^-1，浅层水位下降速率高达0.9 m·a^-1(Zheng et al.，2010;费宇红等，2009).

地下水过度开采不但引起水质恶化、地面沉降等严重环境地质问题(Konikow and Kendy，2005;Sophocleous，2002)，同时，对地壳应力场也产生重要的影响.Holzer(1979)指出地下水开采会造成地球表层弹性扩张和地表抬升，长期的地下水开采则可能引起一些构造活动.已有研究表明：冰、雪等水文环境的变化会对地壳产生加/卸载作用，而这种作用会影响地震活动(Bettinelli et al.，2008;Heki，2003).González等(2012)通过对Lorca地震断层滑动的观测和模拟，得出地下水开采对地壳的卸载作用使断层浅部库仑应力增加，控制地震断层的滑动分布.Amos等(2014)据GPS观测数据认为San JoaquinValley南部地壳垂向上升速率主要归因于地下水开采，且卸载作用降低了正应力，促进微震的发生.

华北平原强震构造发育，是中国大陆强震多发区.近100年来华北地区已发生多次M>5.0地震，其中邢台、唐山等M>7.0地震对人民生活造成了严重影响，给社会经济带来了严重损失.那么，近50年来华北平原地下水过量开采对地壳的卸载作用有多大？是否可能对华北地震活动性造成影响？本研究建立了华北地区二维有限元模型，采用黏弹性Maxwell本构关系，定量计算地下水超采引起的地壳应力变化，进一步探讨地下水开采对华北平原5级以上地震分布的影响.

华北平原地下水沉降漏斗区域主要集中在北京、保定、石家庄和邯郸沿线.参考Stephen和Héctor(2004)给出的1960—2000年间华北平原浅层地下水水位下降分布信息，选取北东—南西向剖面AA′(总长500 km)和北西—南东向剖面BB′(总长200 km)来进行分析研究，见图 1.剖面AA′和BB′分别近似平行和垂直于北京—石家庄主要地下水沉降中心沿线，垂向深度为100 km.为降低固定边界的影响，将AA′和BB′两剖面两端分别延伸500 km和300 km，则AA′剖面计算模型总长1500 km，BB′剖面总长800 km.

图 1

Fig. 1

图 1 华北地区1960年以来水位下降图(Stephen and Héctor, 2004) 线AA′,BB′为计算模型剖面. Fig. 1 Distribution of groundwater decline since 1960 in North China (Stephen and Héctor, 2004) Line AA′, BB′ show the location of model profile.

Zhang等(2011)利用宽角反射剖面对华北地块研究，初步得出华北平原岩石圈厚度为70 km.Tian等(2014)通过对跨越华北克拉通的东西向地震反射/折射剖面数据分析，认为华北平原地壳平均厚度为30 km.Jia和Zhang(2005)据华北地区地震测深剖面数据，发现华北平原上地壳厚度约为20 km，下地壳厚度约为10 km.鉴于上述研究，建立了包括上、中、下地壳及地幔和软流层顶部总厚度为100 km的华北平原分层模型.模型中上、中、下地壳厚度均为10 km，岩石圈地幔40 km和软流层30 km，如图 2a.采用三角形网格对模型进行剖分，AA′剖面模型分辨率为1.5 km，节点总数为70193，单元总数为138250；BB′剖面模型分辨率为1.2 km，节点总数为70686，单元总数为139736.

图 2

Fig. 2

图 2 计算模型及AA′和BB′ 地下水水位变化速率 (a) 岩石圈剖面分层模型及边界条件信息.W为剖面总宽度，AA′剖面W=1500 km，BB′剖面W=800 km; w为上边界加载宽度，AA′剖面w=500 km，BB′剖面w=200 km; (b)、(c) 加载区域地下水位年平均下降速率沿两剖面加载区分布图(Stephen and Héctor, 2004).Fig. 2 Numerical model and groundwater decline rate on lines AA′ and BB′ (a) Layered lithosphere and boundary conditions. W is width of the model, which is 1500 km on line AA′ and 800 km on line BB′, respectively. w is the width of upper boundary which is 500 km on line AA′ and 200 km on line BB′, respectively. (b) and (c) Annual groundwater decline rates on AA′ and BB′, respectively (Stephen and Héctor, 2004).

由于计算模型水平宽度选取了地下水开采区宽度3倍的范围，且垂向计算深度为100 km.考虑到地表卸载对于足够深度且足够远的范围影响较小，边界垂向位移可近似为0，所以模型中两侧边界法向位移约束为0，切向自由；底边界垂向位移约束为0，水平方向自由.考虑地下水开采对地壳产生卸载作用，参考AA′和BB′段的地下水位下降信息，对模型上边界施加垂直向上的法向力.由于本次研究采用Maxwell黏弹性模型，计算时间步长选取为1年，计算总时间为50年，来研究地下水累计开采50年卸载作用对地壳应力的影响.因此，上边界的初边值条件设定为：初始法向应力为0，每个时间步法向应力增量为地下水位年平均变化速率、含水层孔隙度和水的容重之积.两剖面地下水水位年平均下降速率分布见图 2b—2c(Stephen and Héctor, 2004).

由于地下水持续开采50年对地壳、岩石圈产生卸载作用，地球介质在几十年的黏性松弛效应需要考虑，本文采用黏弹性Maxwell本构方程，计算总时间为50年，时间步长为1年.

黏弹性Maxwell体满足以下方程形式：

平衡方程

几何方程

Maxwell体应力-应变(本构)关系

其中：

式中：σ_ij为应力张量，f_i为体力项，u_i为位移分量，ε_ij为应变张量， 为应变速率， 为应力速率，[D]为弹性矩阵，E为弹性模量，ν为泊松比，η为黏滞系数，[Y]为黏性项矩阵.

前人对中国大陆岩石圈流变性质已有大量的研究.例如，Zang等(2005)通过对鄂尔多斯及其邻近地区岩石圈地震波速度结构、热结构和岩石组成等研究，得到鄂尔多斯地块及其周围三维岩石圈流变结构；石耀霖和曹建玲(2008)基于岩石圈温度和应变速率研究成果，利用实验室流变试验结果、计算得出中国大陆岩石圈等效黏滞系数；孙玉军等(2013)利用地热和GPS得到的应变速率数据采用有限元方法计算得到中国大陆的三维流变结构.基于华北地块现有的地震反射/折射数据得到的岩石圈P波和S波速度结构(Huang et al.，2009;Jia and Zhang，2005;Tian et al.，2014)，以及Crust2.0模型给出的地壳密度分布(Christensen and Mooney，1995)，可计算得到岩石圈各分层结构的弹性模量和泊松比等地壳力学参数.表 1给出了各分层力学参数.华北平原含水层主要为细砂和粉砂岩，模型选取含水层岩石孔隙度为0.4(弗里泽和彻里，1987；王大纯等，2006).

表 1(Table 1)

表 1 计算模型采用的力学参数 Table 1 Mechanical parameters of calculation model 厚度 (km) 密度 (kg·m-3) 泊松比 弹性模量 (Pa) 黏滞系数 (Pa·s) 上地壳 10 2760 0.246 8.9×10-10 5.0×10-23 中地壳 10 2810 0.259 9.43×10-10 2.0×10-23 下地壳 10 2910 0.255 1.1×10-11 1.0×10-22 岩石圈 地幔 40 3320 0.296 1.52×10-11 3.0×10-22 软流圈 30 3320 0.296 1.52×10-11 2.0×10-21

表 1 计算模型采用的力学参数 Table 1 Mechanical parameters of calculation model

图 3给出了AA′和BB′剖面上地下水持续50年开采后引起的位移场分布.图 3a和3b分别显示了AA′剖面上垂向和水平向位移分布，可以看出：地下水开采区地表整体抬升，而外围区域的垂向位移较小；开采区内地表垂向位移空间分布变化与水位下降分布相对应，由于保定—邢台沿线存在大规模水位下降，该区地表垂向抬升剧烈，整体大于6 cm.在石家庄附近，地表抬升高达12.4 cm；受地下水开采影响北京地区地表上升达6 cm.垂向上升量随深度的增加而减小，在地下20 km深处，保定—邢台沿线地壳，垂向最大抬升量为7.8 cm.相对垂向位移，地下水开采引起的水平位移较小，地表最大水平位移主要集中在北京附近，约为2.3 cm；而保定—邢台沿线地表水平位移值均小于0.5 cm.图 3c和3d分别显示了BB′剖面上垂向和水平向位移分布，同样地，开采区地表出现明显的抬升，抬升最大值约12 cm，集中在石家庄东南附近的地下水最大沉降漏斗区域；地下20 km深处的最大垂向位移约为9.0 cm；开采区内水平位移大小在0～1.9 cm范围内，地下水超采对开采区外围区域的水平位移影响较小.

图 3

Fig. 3

图 3 地下水持续开采50年引起的位移变化 (a) AA′剖面垂向位移结果; (b) AA′剖面水平位移结果; (c) BB′剖面垂向位移结果; (d) BB′剖面水平位移结果. Fig. 3 Displacement due to groundwater mining in 50 years (a) Vertical displacement on line AA′; (b) Horizontal displacement on line AA′; (c) Vertical displacement on line BB′; (d) Horizontal displacement on line BB′.

岩石圈深部和软流圈流变性质较强，深部对地下水卸载的响应较明显，因此，地下深部水平位移整体大于地表位移.地表垂向位移空间分布与地下水水位下降分布曲线一致.水位下降越大，对应地壳垂向位移越明显.由于地下水赋存于松散沉积层中，地下水开采引起的沉积物颗粒压实、地面沉降现象更为明显.地表沉积物的压实量抵消了地下水卸载对应的地壳抬升作用，所以沉积层覆盖区域地表无法观测到垂向上升.从地壳、岩石圈规模考虑，地下水开采对地壳、Moho面垂向变形的影响因素不可忽视

地下水卸载会引起岩石圈应力变化，对区域构造应力场产生干扰.图 4为地下水持续开采50年后引起的地壳水平应力变化结果.图 4a为AA′剖面上水平应力变化，由图可知，地下水卸载使开采区产生水平向拉张变形，拉应力增量高达89 kPa，位于石家庄北东方向地下水沉降漏斗区.邢台市周围最大拉应力增量约为80 kPa；保定西南附近区域拉应力增量达62 kPa，其东北侧区域受地下水非均匀分布卸载影响，产生局部挤压变形，压应力增量高达20 kPa；北京地区拉应力增量最大为78 kPa.在开采区的外围约50 km范围内，地表产生挤压变形，最大压应力增量达38 kPa.开采区50 km以外，地壳应力受地下水影响甚微.由于华北地区地震震源深度范围为5~30 km，本文分别选取震源深度范围内5 km、15 km和25 km深度，分析不同深度上水平应力变化及其对构造应力的扰动程度，见图 4b.开采区地壳深部整体受拉张作用；北京附近区域，地下5 km深处最大水平拉应力达41 kPa，15 km深最大水平拉应力增量为21 kPa，地下25 km深处拉应力达15 kPa；石家庄附近水平最大拉应力增量在5 km深度高达49 kPa，15 km深度为25 kPa，25 km深度为18 kPa；邢台地区地下5 km、15 km和25 km处最大水平拉应力增量分别为36 kPa、25 kPa、18 kPa；保定地区相对其他地区，受卸载影响地壳应力变化较小.图 4c为BB′剖面水平应力分布，该剖面上开采区受地下水卸载作用整体呈现拉张变形，水平拉应力增量高达118 kPa，分布在石家庄东南方向地下水最大漏斗区.开采区外围一定范围同样产生压缩变形，压应力增量高达38 kPa.图 4d为BB′剖面不同深度水平应力变化，该剖面5 km、15 km和25 km深度对应最大水平拉应力增量分别为72 kPa、55 kPa和31 kPa.受地下水水位下降空间分布差异影响，BB′剖面上开采区应力变化高于AA′剖面.

图 4

Fig. 4

图 4 地下水持续开采50年引起的水平正应力分布 (a) AA′剖面； (b) AA′剖面不同深度； (c) BB′剖面； (d) BB′剖面不同深度. Fig. 4 Horizontal normal stress changes due to groundwater mining in 50 years (a) On line AA′; (b) At different depths on line AA′; (c) On line BB′; (d) At different depths on line BB′.

图 5为受地下水持续50年开采影响AA′和BB′剖面上垂向应力变化结果.地下水开采区岩石圈垂向上整体表现为拉张变形，区域垂向拉应力最大值集中在地下水漏斗中心区，地表垂向拉应力增量与地下水水位下降呈正比.最大垂向拉应力增量达227 kPa，位于石家庄东南方向最大沉降漏斗中心区域，邢台、石家庄地区垂向拉应力增量高达170 kPa.北京地区垂向拉应力增量高达145 kPa.垂向应力变化随深度逐渐减小；开采区外围受地下水开采影响不明显，垂向应力增量小于5 kPa.

图 5

Fig. 5

图 5 地下水持续开采50年引起的垂向正应力变化分布 (a) AA′剖面； (b) BB′剖面. Fig. 5 Vertical normal stress changes due to groundwater mining in 50 years (a) On line AA′; (b) On line BB′.

受太平洋板块和菲律宾板块向欧亚大陆的俯冲联合作用，华北平原现代构造应力场的主体特征表现为北东东—南西西方向的挤压和北西—南东方向的拉张作用(谢富仁等，2004).主压应力方位为北东至北东东方向，中间主应力基本是垂直的，最大主应力和最小主应力接近水平，处于走滑应力状态(陈连旺等，2001;崔效锋和谢富仁，2001;邓起东等，1979)，华北地区构造压应力积累为0.5 kPa·a^-1(柳畅等，2012;朱守彪等，2010)，而本文计算得到华北地区地下水开采引起的水平拉应力增量在震源深度范围高达70 kPa.对比区域构造0.5 kPa·a^-1的压应力增量，地下水开采对华北地区区域应力场的扰动不可忽视.近50年来地下水的累积卸载有效降低了华北平原北东—南西方向的主压应力构造积累，同时又增加了北西—南东向区域性的拉应力.

地震的发生受现今区域构造应力状态控制(于慎谔等，2000).改革开放以来随着经济的发展，水资源需求迅速增加，地下水水位下降速率明显增加.即近50年来地下水位下降主要集中在1980年以来.参考研究区1980年前后35年两个时间段5.0级以上地震分布，见图 6，研究区内1980年之前地震事件数目明显多于1980年之后：河北宁晋—束鹿一代，1966—1968年间发生了近20个5级以上地震，震源深度大于10 km；1980之后，附近区域仅发生一次震源深度为20 km的地震.唐山地区，于1976年发生7.8级地震，其后4年伴随有一系列余震；1980年之后，仅在1991年和1995年发生两次5.0级地震.由于构造应力是控制地震发生的主要驱动力，构造应力积累造成地震发生的周期性.华北地区地震构造应力积累受多种因素干扰，2011日本东北M9.0地震对华北地区压应力释放达1.0 kPa，其应力降抵消了5～10年的构造应力积累(张贝等，2015).而华北地区地下水几十年内的累积超采引起震源深度范围压应力释放达70 kPa，有效降低了区域构造应力积累.由于地震复发周期具有不确定性，通过20世纪80年代前后地震空间分布的差异，本文不能断定地下水卸载减少了华北地震的发生.但考虑到卸载引起的地壳应力变化，本文认为地下水开采是干扰构造应力积累的重要因素之一，地下水开采对岩石圈的卸载过程可能在一定程度上减缓了地震的孕育过程.

图 6

Fig. 6

图 6 1980年前后35年时间段内5级以上地震分布 (a) 1945—1979年前地震； (b) 1980—2015年地震. Fig. 6 Earthquake (>MS5.0) epicenter distribution of two periods before and after 1980 (a) Earthquakes in 1945—1979； (b) Earthquakes in 1980—2015.

综上，华北地区持续50年的地下水开采对区域地壳抬升和应力产生重要的影响，同时在一定程度上影响到华北地区地震活动性.具体表现在以下几个方面：

(1)华北地区地下水过量开采对岩石圈的卸载作用会引起地表微小的水平位移；而地表垂向上升可达12.4 cm，垂向位移随深度的增加逐渐减小，卸载作用对开采区以外的区域垂向位移影响甚微.

(2)地下水超采引起开采区整体拉张变形：震源深度处，水平拉应力增加高达70 kPa；对局部区域产生挤压作用，压应力增量达17 kPa.开采外围出现局部小范围压缩变形：震源深度范围内，水平压应力增量达20 kPa.而华北平原构造压应力积累约为0.5 kPa·a^-1，可见，地下水过量开采整体上有效释放了华北地区构造主压应力积累，可能对华北地区地震孕震过程有减缓作用.

由于本次模拟采用二维模型模拟计算地下水开采对岩石圈变形的影响，不能精确地给出受地下水卸载影响华北地区三维地壳应力场及断层库仑应力的变化，该工作有待于三维数值模拟研究解决.