0 引言

苏锡常地区地面沉降始于20世纪70年代中期，由于规模化开采地下水，地下水位不断下降，地面沉降现象逐渐显现；20世纪80—90年代，地下水位漏斗面积不断增大，地面沉降范围由中心城区逐渐扩大至外围区域；20世纪90年代实行地下水限采措施以来，地下水环境改善效果显著，大部分地区地下水位回升幅度明显^[1-2]，地面沉降得到有效控制，地面沉降速率大幅减小。地面沉降的成因是多方面综合作用的结果，地下水资源超采是前期地面沉降发展的主因^[3]。但是苏州吴江盛泽地区由于近年来城镇化的快速发展，一方面建筑荷载短期内急剧增加，由建筑物基础桩的负摩擦力引发的地面沉降已成为一个不可忽视的重要因素，另一方面还存在较大的地下水位降落漏斗，地面沉降仍在发展之中。目前，地面沉降速率仍高于15 mm/a，成为苏锡常沉降速率最大的地区之一^[4]。因此，全面分析地面沉降的影响因素，对地面沉降发展趋势进行准确模拟预测，为科学防治地面沉降提供有效的解决方案已成为当务之急。

目前，国内外关于地面沉降模拟和预测研究的热点集中以地下水开采为主因^[5-9]，对于地面荷载引起的地面沉降研究处于初级阶段^[10-12]，针对此方面的研究较少。已有的研究都是将高层建筑荷载和地下水开采作为独立事件进行分析，未考虑两者间存在的耦合效应对地面沉降的影响^[13]。本文基于比奥固结理论^[14-15]，建立了吴江盛泽地区建筑荷载和地下水开采与地面沉降三维全耦合模型，模拟预测了建筑荷载单独作用、地下水开采单独作用及建筑荷载和地下水开采叠加作用下盛泽地区地面沉降发展趋势，分析二者对地面沉降造成的影响，以期为吴江盛泽地区地面沉降的科学防控提供决策依据。

1 研究区概况

地处江苏省苏州市吴江区南部的盛泽位于30°50′35″N—30°57′52″N，120°30′37″E—120°40′13″E，为江、浙、沪交界处，包括南麻、坛丘等，行政总面积150 km²。区内第四纪沉积层一般厚度为140~200 m，主要为河湖相、河口相沉积。主要含水层为全新统孔隙潜水含水层、上更新统第Ⅰ孔隙承压含水层、中更新统第Ⅱ孔隙承压含水层、下更新统第Ⅲ孔隙承压含水层，主采层为第Ⅰ承压含水层。含水层组分布详见表 1。

区内浅部地层中由高压缩性土层构成的软土相对发育，厚度一般为25~30 m，岩性、岩相变化不稳定，工程地质条件较差。

2 地质概念模型 2.1 边界条件

研究区积累了多年的钻孔勘探资料和地下水动态观测资料，根据地下水动态及渗流场特征，将研究区的侧向边界概化为流量边界。研究区的顶部主要接受大气降雨的补给，是补给边界，同时地下水又通过其蒸发，故顶部也是排泄边界。研究区的底部为一固结黏土层，概化为隔水边界。图 1为研究区水文地质剖面布置图。

2.2 地层结构特征

研究区第四纪地层厚度一般为140~200 m，垂向上为砂层和黏土层互层，平面上各土层的颗粒结构也极不均匀，水力学性质和土力学性质存在较大的差异，为非均质各向异性。图 2为Ⅰ-Ⅰ′南北向水文地质剖面略图。

2.3 地下水流态特征

由于早年间研究区不加节制地开采地下水，各含水层垂向产生较为强烈的水力联系，并形成降落漏斗。加之地下水开采和大气降雨季节性的影响，地下水位也发生强烈的波动，由此将地下水流态概化为三维非稳定流。

2.4 源汇项和建筑荷载

研究区的源汇项主要是地下水开采、大气降雨入渗补给与地下水蒸发排泄。地下水开采概化为大井，大气降雨入渗补给和地下水蒸发排泄分别概化为入渗补给强度和蒸发排泄强度，并结合在一起综合处理。将地面建筑荷载概化为点荷载处理。研究区井位分布情况见图 3，主要地面建筑荷载分布见图 4。

3 数学模型 3.1 比奥固结理论

本次研究以三维比奥固结方程为基础，建立了盛泽地区建筑荷载和地下水开采与地面沉降三维全耦合数学模型：

式中：G为剪切模量；▽为拉普拉斯算子；ω_x、ω_y、ω_z分别为x、y、z方向的位移；υ为泊松比；u为超静孔隙水压力；γ为土的重度；t为时间；K_x、K_y、K_z分别为x、y、z方向的渗透系数；γ_w为水的重度；W为源汇项。

式(1)可以简写成：

式中：ω=[ω_x, ω_y, ω_z]^T；D为应力应变关系矩阵；f=[f_x，f_y，f_z]^T，f_x，f_y，f_z分别为体积力在x, y, z方向上的分量；M=[1, 1, 1, 0, 0, 0]^T；Q为计算点的等效流量；C为弹性柔度矩阵；B为应变位移矩阵；K为固体刚度矩阵。

3.2 土体本构关系

根据土体的流变特性，将土体本构关系拓展到黏弹塑性，则该土体的任一点在任一时刻的应变增量^[16]可表示为

式中：dε为总应变增量；dε^ep为弹塑性应变增量；dε^ve为黏弹性应变增量；dε^vp为黏塑性应变增量。

3.3 定解条件 3.3.1 初始条件

1) 孔隙水压力初始条件

孔隙水压力初始条件为

式中：u₀(x, y, z)为初始孔压；h₀(x, y, z, t₀)为初始水头值；ρ_w为水的密度；g为重力加速度；Ω为整个研究区域。2014-09-01各含水层组的初始流场为该模型的初始水位值。

2) 地应力初始条件

已知土体的自重应力，可以估算土体的初始应力：

式中：σ_x、σ_z为土体的初始水平向和垂向应力；z′为各层计算点深度；k₁为静止侧压力系数。

式中，φ′为有效内摩擦角。

3) 位移初始条件

模型的位移初始条件为2014-09-01各含水层的初始位移值，初始位移值设为0。

3.3.2 边界条件

1) 流量边界条件Γ₂：

式中：K为渗透系数；H为水位；n为孔隙度；q_L为边界Γ₂上的已知单位面积流量。

2) 自由面边界条件Γ₃：

式中：μ为土体给水度；θ为自由面外法线方向与垂线的交角；q为通过自由面边界Γ₃的单位面积流量；Z为自由面所在的海拔高度。

3) 位移边界条件Γ₄：

式中：为位移边界Γ₄上3个方向的已知位移，本次模拟取3个方向位移为0。

3.4 参数动态变化 3.4.1 孔隙度与渗透系数的非线性

建筑荷载和地下水开采作用下产生孔隙水应力消散，土体骨架应力进而重新分布，土体固结变形为其宏观表现。基于比奥固结理论，结合孔隙度相关定义和Kozeny-Carman方程，得到渗透系数K和孔隙度n之间的关系式^[17]：

式中：n₀为初始孔隙度；K₀为初始渗透系数；ε_v为体应变,

3.4.2 弹性模量和泊松比的非线性

根据Duncan-Zhang提出的非线性模型，将土体本构关系推广到非线性，应力应变关系矩阵D的弹性常数E和υ为随应力状态改变的变量，则切线弹性模量E_t和切线泊松比υ_t的表达式为：

式中：p_a为大气压强；σ₁为第一主应力；σ₃为第三主应力；R_f为破坏比；φ为内摩擦角；c为黏聚力；v₀为围压为100 kPa时的初始泊松比；E₀为σ₃=100 kPa时的初始弹性模量；F为试验常数，本文F=0.04；m为弹性模量与固结压力曲线的斜率；D为土体试验参数，D=3。

4 模型求解与校正

根据前述数学模型，运用伽辽金有限元方法求解，采用Fortran95语言编制了建筑荷载和地下水开采与地面沉降三维全耦合数值分析程序^[18]，如图 5所示。

研究区平面上共计剖分成4 635个矩形网格，垂向上剖分成7个层位，平面尺度为18 481 m×13 578 m，垂向深度为140~200 m，面积约为150 km²。经过离散的模型中共有节点38 888个，单元32 445个。三维剖分网格见图 6。选取为期1 a(2014-09-01— 2015-08-31)的观测资料为研究基础，对模型参数进行反演计算，以每月作为1个应力期，1个应力期对应1个时间步长，共12个应力期。

根据盛泽地区的具体情况，降雨量采用年平均降雨量，经模型调试，扣除蒸发量以后的有效降雨入渗补给系数确定为0.014，有效入渗面积为72.27 km²。研究区共有41口井，地下水开采量在调查统计的基础上，结合现有地下水位降落漏斗形态特征，经模型反演计算确定。单井日开采量为50~200 m³/d，日均抽水量约为5 710 m³/d，其中镇区中南部茅塔、杨扇、圣塘、大谢处日均抽水量约为4 150 m³/d，盛虹社区以北思安桥处日均抽水量约为1 560 m³/d。本程序中流量的计算原理为，从外部水源流入计算单元的水流量Q与外部水源的水头h₂和该计算单元水头h₁之差成正比，公式为

式中：A为过水断面面积；L为两单元中心之间的距离。结合各应力期边界上地下水进出量，经模型反演计算确定边界条件。

建筑荷载所产生的地面沉降量包括基础上部荷载产生的沉降量和相邻基础荷载因应力扩散产生的沉降量。根据实地调查的盛泽地区建筑物荷载分布特征情况确定建筑荷载，建筑荷载按每层自重1.2 t/m²、每层承重0.3 t/m²计算。将实际建筑荷载区域概化为矩形区域，计算概化后的矩形区域建筑荷载影响下的地基附加应力，得出影响范围内模型中各节点的附加应力。则矩形基础受垂直分布荷载作用下地基土中任意节点的竖向附加应力系数a_sr为：

其中：

式中：l为矩形的长；b为矩形的宽；x、y、z为计算点坐标值。

则任意节点的附加应力为

式中, P_d为高层建筑荷载。

全区有9口观测井进行水位拟合，结果见图 3。经调参、反演后，模型共7层38个分区。以第Ⅰ承压含水层为例，图 7为分区情况，各分区的水文地质参数值详见表 2。

模型识别计算所得的各含水层第12应力期的水均衡要素见表 3，可见地下水流入量和流出量基本相当。研究区各水准点2014-09—2015-08地面沉降拟合情况见表 4。以盛泽镇盛溪河泵站水位为例，其水位拟合图见图 8。由表 4和图 8可见，地面沉降量计算值与实测值拟合精度较高。累计地面沉降量计算值与实测值拟合图见图 9，可见累计地面沉降量计算值与实测值吻合程度较高，水准观测点沉降计算值和实测值模拟误差均小于10%。综上，模型收敛性及稳定性较好，模型识别所得结果可信，精度较高。

图 10为2015-08-31第Ⅰ承压含水层计算流场图，根据研究区建筑荷载和抽水井分布情况，结合图 8，可见模拟计算得到的地下水流场变化符合研究区内实际地下水流场的宏观流动变化规律。

5 地面沉降发展趋势预测

模型计算周期为2015-09-01—2030-08-31，以1 a作为1个应力期，每个应力期为1个时间步长，共15个应力期。分别模拟建筑荷载单独作用、地下水开采单独作用以及建筑荷载和地下水开采叠加情况下地面沉降趋势，根据不同情况下2015—2030年预测地面沉降量分布图，进行对比分析。

5.1 建筑荷载对地面沉降的影响

在地下水禁采、地面建筑荷载单独作用下，2015-09-01—2030-08-31，盛泽地区从上至下的潜水含水层、3个承压含水层及其间的3个黏性土弱含水层共计7个层位中，除潜水含水层、第Ⅲ黏性土弱含水层部分地区水位下降外，其余含水层水位均呈现上升趋势。

停止地下水开采后，盛泽地区地下水位有一定程度的回升，但建筑荷载继续作用于各含水层，使其产生一定程度的压缩变形，从而引发地面沉降的发生。从2015-09-01—2030-08-31，盛泽地区从上至下7个层位的最大沉降量分别为27.99、190.28、87.97、74.10、27.28、24.63、9.86 mm，累计最大地面沉降量为442.11 mm。第Ⅰ黏性土弱含水层是建筑荷载引发土体变形的主压缩层，单层压缩量占比43.04%。

杨扇、龙桥、坛丘等建筑荷载较密集的地方为盛泽地区建筑荷载单独作用引发地面沉降的主要集中区域。远离建筑荷载及无建筑荷载分布的区域，地面沉降量逐渐减少，并趋近于0。此外，建筑荷载对土体变形的影响随着深度的增加也逐渐减弱。建筑荷载单独作用使盛泽地区出现漏斗状的地面沉降，以建筑物为中心向四周辐射。地面沉降量与影响范围与建筑荷载密集程度有关，建筑荷载分布越密集，地面沉降量越大，影响范围越广。建筑荷载单独作用下2015—2030年预测地面沉降量分布见图 11。

5.2 地下水开采对地面沉降的影响

假定第Ⅰ承压含水层作为开采层位，维持2014-09-01—2015-08-31单井开采量及开采布局不变。2015-09-01—2030-08-31，在不考虑建筑荷载，仅有地下水开采的情况下，盛泽地区从上至下7个层位中仅有潜水含水层、第Ⅲ承压含水层水位上升，其余层位水位均呈现下降趋势。地面沉降量最大的杨扇地区，2015—2030年地下水水位变化情况见图 12。由图 12可见, 杨扇地区2016年起地下水位进入平稳阶段，下降极其缓慢，说明地下水开采对地面沉降影响小于建筑荷载的影响。茅塔、杨扇、圣塘等纺织工厂分布较密集的地方为盛泽地区地下水开采单独作用引发的地面沉降主要集中区域，地下水开采单独作用引发的地面沉降漏斗分布范围与地下水位降落漏斗分布范围基本吻合，表现出了时间和空间的一致性。

从2015-09-01—2030-08-31，茅塔社区从上至下7个层位的最大沉降量分别为9.74、32.97、101.74、28.32、8.11、4.47、2.84 mm，茅塔社区为各层最大沉降量产生区域，累计最大地面沉降量为188.19 mm。第Ⅰ承压含水层是地下水开采引发土体变形的主压缩层，单层压缩量占比54.06%。地下水开采单独作用下2015—2030年预测地面沉降量分布见图 13。

5.3 建筑荷载和地下水开采叠加作用对地面沉降的影响

将现有建筑荷载和地下水开采的作用同时考虑，第Ⅰ承压含水层仍作为假定开采层位，保持2014-09-01—2015-08-31的假定开采布局及开采量不变。盛泽地区从上至下7个层位中仅有潜水含水层、第Ⅲ承压含水层水位上升，其余层位水位均呈现下降趋势。

从2015-09-01—2030-08-31，盛泽地区从上至下7个层位的最大沉降量分别为40.46、276.06、189.26、82.88、27.41、21.54、15.01 mm，累计最大地面沉降量为652.62 mm。第Ⅰ承压含水层及其上覆第Ⅰ黏性土弱含水层是建筑荷载和地下水开采叠加作用引发土体变形的主压缩层，二者压缩量之和占总压缩量的71.30%。

盛泽地区建筑荷载和地下水开采叠加作用下引发的地面沉降主要集中在杨扇、茅塔、圣塘等镇中南部区域，且形成了以杨扇社区为中心的沉降漏斗区，与建筑调查密集区和地下水位降落漏斗区分布一致。建筑荷载和地下水开采叠加作用下2015—2030年预测地面沉降量分布见图 14。

5.4 结果分析

盛泽部分地区不同条件下累计最大地面沉降量和地面沉降速率如表 5所示。根据表 5可知，除茅塔和圣塘外，盛泽地区建筑荷载引发的地面沉降量要大于地下水开采引发的地面沉降量，2015-09-01—2030-08-31，大部分区域建筑荷载和地下水开采对地面沉降的贡献比约为3.3:1，局部地段约为2.5:1，可知建筑荷载是造成盛泽地区地面沉降的主因，其次为地下水开采。此外，由表 5还可得出，盛泽地区建筑荷载和地下水开采单独作用下引发的地面沉降量的线性叠加之和大于二者叠加作用下引发的地面沉降量，说明建筑荷载和地下水开采叠加作用时产生地面沉降耦合效应。

以最大地面沉降速率为分析指标，由表 5可知，盛泽地区除茅塔和圣塘外，建筑荷载单独作用引发的地面沉降速率均大于地下水开采单独作用引发的地面沉降速率。建筑荷载和地下水开采单独作用下引发的地面沉降速率的线性叠加之和大于二者叠加作用下引发的地面沉降速率，这反映了建筑荷载和地下水开采叠加作用引发地面沉降时具有耦合效应。

6 结论

1) 以比奥固结理论为研究基础，建立了盛泽地区建筑荷载和地下水开采与地面沉降三维全耦合模型，该模型对盛泽地区不同情况下地面沉降发展趋势的模拟结果较为准确，与实际情况相吻合。

2) 第Ⅰ黏性土弱含水层和第Ⅰ承压含水层分别为建筑荷载和地下水开采单独作用下的主压缩层，单层压缩量占比为43.04%和54.06%；第Ⅰ承压含水层及其上覆第Ⅰ黏性土弱含水层是建筑荷载和地下水开采叠加作用引发土体变形的主压缩层，二者的压缩量之和占总压缩量的71.30%。

3) 建筑荷载对盛泽地区地面沉降的影响大于地下水开采。从2015-09-01—2030-08-31，大部分区域建筑荷载和地下水开采对地面沉降的贡献比约为3.3:1，局部地段约为2.5:1。

4) 盛泽地区建筑荷载和地下水开采单独作用下引发的地面沉降量及地面沉降速率的线性叠加之和均分别大于二者叠加作用下引发的地面沉降量及地面沉降速率，叠加作用引发的地面沉降具有耦合效应。