基于GMS-SUB的地面沉降数值模拟研究——以太谷东站为例

惠 翔, 刘丽红

(安徽理工大学地球与环境学院，安徽 淮南 232001)

地面沉降在全世界范围内出现频繁且影响范围广，不仅会造成城市建筑物、道路交通、管道系统线路发生不均匀沉降而被破坏，还会给排水、防洪、城市内涝的防治带来诸多困难[1-2]。据统计，我国自从新中国成立以来至21世纪初由于地面沉降造成的损失高达4 500～5 000亿元，年均总体损失90～100亿元，严重影响经济社会发展[3]。我国已有超过17个省、96个城市出现明显的地面沉降，覆盖面积超过1.0×105km2[4-5]。随着近年来社会经济的快速发展，尤其是在第四纪沉积层广泛分布的地区，大规模地进行地下水抽采是引发大面积地面沉降的主要原因[6-7]。 目前， 国内外关于沉降的计算理论和预测方法，在构建形式上主要分为两类： 第一类是基于研究区详细统计资料， 通过沉降量与其相关因素的关系建立数理统计模型和随机统计模型； 第二类是确定性的地下水-地面沉降水土模型[8-11]。

数理统计模型一般包括回归分析模型和灰色理论模型。文献[12]利用对天津地区历史开采水量、历史沉降量数据的整合与分析，建立了地下水位与地面沉降量之间的线性回归方程，之后对天津地区未来9a地面沉降趋势作出预测。文献[13]依据太原市5个沉降中心观测点的沉降量建立了灰色Verhulst预测模型对未来沉降进行了预测。数理统计模型和随机统计模型重在考虑地下水开采量与地面沉降量之间的数量关系，并没有考虑土体内部水体与土体骨架间的力学性质和沉降内在机理，也没有考虑地下水流场的动态演化特征，模型不能很好地反映沉降内在动力学机制。

水土模型分为两步走计算模型、部分耦合模型、完全耦合模型。1989a上海地质勘察院与比利时合作使用两步走模型模拟了上海地区的沉降[14]。两步走模型是将水流模型与沉降模型分开考虑的，不能充分考虑土体不均一性以及沉降水土耦合情况。文献[15]在比奥固结理论的基础上建立水土全耦合模型模拟地面沉降，这类模型计算参数复杂，计算耗时长，模型难以稳定并不适合用于大区域沉降模型计算。

部分耦合模型可以很好地模拟含水层抽水引起的水位降深，并建立起含水层垂向渗透系数和贮水率与孔隙率之间的关系，水位变化时刻影响着沉降，沉降会导致含水层的压缩，同时孔隙水压力降低含水层水体膨胀从而含水层释水，释水量当作源汇项又添加到水流方程中，从而完成水土耦合。部分耦合模型能够模拟多种含水介质的变形并且模型涉及的参数简洁，计算明确。目前大多数抽水型地面沉降模型是基于太沙基有效应力原理建立的地下水三维渗流模型与一维垂向固结沉降的数学模型，本文在研究区的每个网格剖分点位置都能计算出随时间变化的水位降深，通过水位的降深和水量的变化得出每个网格点的垂向一维沉降量，每个点的沉降扩展至平面得到平面分布的三维沉降场，从而将地下水动力的三维模型与固结沉降的一维模型进行耦合。本文使用GMS-SUB模块建立地下水渗流与地面沉降部分耦合模型对太谷东站研究区地下水位和地面沉降的动态变化进行模拟，为定量化研究地面沉降问题提供一种可靠的手段。

1 研究区概况

太焦铁路太谷东站位于山西省晋中市太谷区胡村镇的桑梓村和朝阳村之间，自太谷区城东北约6km处由东南向西北穿过。车站地理坐标为37°40′～38°00′N，112°25′～112°45′E。研究区位于太原盆地中部偏东，地貌属冲洪积倾斜平原区，地形平坦开阔，地势总体东南高、西北低，地面高程在826.4～788.1m之间，最大相对高差38.8m(见图1)。

图1 研究区位置图

1.1 地层条件

太原盆地是山西省中部一个大型新生代断陷盆地，内部沉积有巨厚的新生界地层。太谷东站研究区有复杂的第四系松散沉积层，通过地质钻孔揭露，深度达230m。第四系地层在研究区内广泛分布，总厚160～175m，区内自上而下依次分布第四系下更新统河湖相沉积物、中更新统冲洪积物、上更新统冲洪积物、全新统冲洪积物，本次研究重点层位为第四系。

1.2 地质构造

研究区位于太原盆地东侧侯城-平遥陷隆内，基底总体上是一个复向斜构造，轴向北东。距离东阳-净化断裂约2km，该断裂是控制盆地东部的一条主要断层，其总体延伸方向为北东-南西向，属压性兼扭性断裂。

1.3 水文地质特征

第四系松散沉积层自上而下可以分为四个含水岩组。第一层为全新统潜水-微承压含水层，由于垂向上含水层无较好的隔水性能，不具有承压水性质，地下水水位呈自由水面，在局部地带含水层与黏性土分布连续且稳定，微具有承压水性质，现已疏干，富水性极弱，为透水不含水的透漏水层。第二层为上更新统第一承压含水层，由于浅层水超采严重，导致浅层水位持续下降出现疏干含水层状况，形成区域性超采漏斗区。第三层为中更新统第二承压含水层，此层水为中层水，地下水以上部越流及降水入渗补给为主，径流条件好，交替循环性好，但随着地下水集中开采，人为改变了地下水流场。第四层为下更新统第三承压含水层，为深层水，地下水补给主要为其他含水层越流补给。承压含水层之间弱透水层相对于含水层厚度较薄，并且很多位置土层不连续存在缺失。中深层水是研究区混合集中开采层位，所以模拟的重点为第二、第三承压含水层中深层地下水。

2 太谷东站地面沉降数值模拟模型

2.1 概念模型的建立

根据提供的水文地质资料，潜水-微承压含水岩组和第一承压含水岩组均已经达到疏干状态，因此目前研究区主要开采中深层地下水，数值模拟主要模拟的含水层为第二承压含水层和第三承压含水层。模拟研究区面积约24.3km2，根据水文地质条件概化结果对研究区进行网格剖分，垂向上剖分为2层，平面上将单元格剖分为200行200列，共计80 000个单元格，平面和垂直剖分图如图2所示。研究区外侧为无效单元格，其他网格为模型计算范围单元格，有效单元格数为41 462个，无效单元格数为38 538个。

图2 研究区网格剖分图

2.2 地下水渗流与地面沉降耦合数学模型

1)地下水三维渗流数学模型 将研究区地下水流概化成非均质非稳定各向异性三维地下水流系统， 并建立相应的数学模型， 该运动微分方程如下

(1)

式中：kx、ky、kz为渗透系数在x、y、z方向的分量(假定渗透系数主轴方向与坐标轴的方向一致m·d-1)；w为单位体积流量，用以代表流进源或流出汇的水量；μs含水层的储水率，m-1；D为模拟区范围；Γ1为第一类给定水头边界条件；Γ2为第二类给定流量边界条件；q为边界流量。

2)地层垂向一维变形模型 在一维垂直方向上，地面沉降遵循有效应力原理。在某一地下水位水平面上，由于大量开采地下水导致含水层孔隙水压力减小，含水层有效应力增大，从而引起含水层骨架的压缩。有效应力的变化能够引起组成含水层的沉积层压缩及膨胀，沉积层的压缩主要是含水层弱透水夹层厚度的减少，一维可压缩系数α可定义为

(2)

式中：db为压缩层厚度的变化量，m3；dσ有效应力变化量，Pa。

如果有效应力的变化仅仅是由于孔隙水压力的变化，则可以

(3)

式中：Ssk=ρwgα为骨架储水率，m-1；Sk=Sskb为骨架储水系数，无量纲；dh为水头变化量，m。

则压缩量可进一步表达为

Δb=SkΔh

(4)

式中：Δb为夹层厚度变化量，m；Sk为骨架储水系数；Δh为水头变化量，m。

3)水流模型与土体模型耦合 三维地下水流的运动方程为

(5)

式中：Kxx、Kyy、Kzz分别为x、y、z方向上的渗透系数，m·d-1；h为水头，m；W为源汇项，m·d-1；Ss为储水率，m-1。可压缩夹层的水量变化可以添加到等式右边，同样也可以并入到源汇项W中。由于含水层抽水导致弱透水夹层水头降低，有效应力增大导致地面沉降，之后夹层压缩会将夹层内的水排出，进入地下水流系统中，同时沉降方程与水流方程通过水头项进行耦合，由水流方程求解出水头变化值再代入沉降方程计算变形值，进而完成部分耦合。

2.3 地下水开采

根据已有地下水开采资料，研究区对中深层含水层进行抽水，模拟识别验证期概化共计31个抽水井，单井抽水量为1 125m3/d。

2.4 各含水层参数的确定

本研究区范围面积不大，水平渗透系数不存在较大波动，根据对当地土质的分析，确定了大致含水层水平渗透系数范围，同时经过不断地调整参数来调整模型。根据调参过程不断拟合模型大致分为两个参数分区。通过水文资料的收集以及模型调参验证(见图3)，各含水层参数如表1所示。

表1 含水层参数取值表

图3 研究区水平渗透系数分区图

2.5 边界条件

将模型东部西部边界作为隔水边界，南北部边界概化为一般水头边界。

2.6 夹层参数的确定

利用研究区钻孔剖面柱状图将第二承压含水组划分为3个压缩沉降层，将第三承压含水组划分为3个沉降压缩层。夹层弹性储水系数Sfe由弹性储水率Sske乘以夹层厚度得到，夹层非弹性储水系数Sfv由非弹性储水率Sskv乘以夹层厚度得到。夹层的非弹性储水率Sskv由土工试验测得的原位土压缩系数和式确定。

(6)

式中：av为压缩系数，kPa-1；e0为初始孔隙比；rw为水的重度，kN·m-3。

Sske取比Sskv少一个数量级(见表2)。

表2 非滞后夹层计算参数取值表

滞后夹层是用于表示在模拟中水头差消散时间显著大于时间步长的夹层。由于任何含水层都可能包含大量不同厚度的夹层，为减少计算量，可以将同一模型层内具有相同垂直渗透系数和弹性与非弹性骨架储水率的滞后夹层合为一个滞后夹层系统。垂直水力扩散率相似的n个滞后夹层系统的等效厚度bequiv可计算为

(7)

为了保证相同的夹层总量，为滞后夹层系统提供正确的压缩量，需要将压缩量和与周围含水层交换的水量乘以该比例因子nequiv

(8)

可以计算出夹层系统的总压缩量的大小和时间历程。因此，方程只求解了一个等效厚度bequiv的夹层，计算总的夹层压缩量需要乘以nequiv。滞后夹层各项参数如表3所示，先期固结水头取模型识别验证期初始水头。

表3 滞后夹层计算参数取值表

2.7 时间离散

模拟的识别验证期为2015a 8月～2020a 8月，识别验证阶段是根据提供的水文地质条件，研究区地下水位持续下降，局部降落漏斗逐步形成阶段。模型一共分为60个应力期，每个应力期分为3个时间步长，总计180个时间步长。

2.8 模型的识别验证

本研究通过山西省第三勘察院多年来对晋中盆地地下水动态监测成果资料，获得1995a至今地下水位下降速率约为1.11m/a。由模拟结果可以看出，在2016a 8月第二承压含水层已经有部分地区水位到达含水层顶板了，中深层水降落漏斗已经初步形成，降落漏斗位于太谷东站向东1公里处的桑梓村，漏斗面积可达3.26km2，漏斗中心水位埋深48.7m。到2020a 8月研究区水位进一步下降，降落漏斗也逐步向外扩大，漏斗面积达到4.16km2，漏斗中心水位埋深53.3m。模拟区地下水流场与2020a 8月中深层等水位线基本吻合，水位变化趋势不断接近真实观测水平，参数率定达到最优，得到了较为理想的模型模拟识别验证结果，此模型可以用来模拟太谷东站研究区地下水流场分布，为之后与地面沉降模型的耦合提供了保障。

中深层地下水降落漏斗中心与地面沉降中心基本一致，且降落漏斗形状扩展范围与沉降中心范围也基本相似，两者之间存在高度相关性。根据研究区INSAR卫星对研究区2019a 8月～2020a 8月地面沉降监测结果来看，太谷东站位于一个沉降漏斗中心，沉降量范围在20～26mm， 沉降漏斗范围9.36km2。将INSAR监测沉降线与模型运行结果比较，模型整体拟合了沉降趋势[16-17]。综上所述，该模型可以用来预测研究区地面沉降(见图4～图5)。

图4 含水层地下水位图(单位：m)

图5 2019a 8月～2020a 8月地面沉降量等值线图(单位：m)

3 谷东站地下水位与地面沉降预测

利用验证后的GMS-SUB水土耦合模型，以2020a 8月作为预测的初始时刻，模拟2025a 8月的地下水水位变化和地面沉降趋势。现取6个沉降预测点：太谷东站点、东站东北方向500m处的中药原料加工厂点、东北方向1.8km处的董村点、西南方向1.3km处的朝阳中学点，正东方向1km处的桑梓村点、东南方向1.6km处的沙沟西阁点。 6个预测点水位降深分别为8.1m、7.2m、 4.2m、 2.6m、7.5m和6.9m(见图6)。

图6 各沉降预测点水位预测图

采用经验公式法计算预测点的沉降量。各地层压缩量公式为

ΔSi=MiHiΔp

(9)

式中：ΔSi为不同岩性地层的压缩量，m；Mi分别为砂土、粘土压缩率，m-1；Hi为砂土、粘土的层厚，m；Δp为水位降深，m。

2025a 8月地下水位预测等值线和地面沉降预测等值线如图7～图8所示。相同水位变幅下经验公式法与GMS-SUB计算的累积沉降量如表4所示，经验公式与GMS-SUB数值模拟预测的沉降量值较为接近，建立的模型可靠，同时随着地下水的逐年开采，沉降漏斗范围有继续扩大的趋势。

表4 相同水位降深下经验公式与GMS-SUB计算沉降量

图7 含水层预测水位线(单位：m)

图8 地面沉降预测等值线(单位：m)

4 结论

(1)依据太沙基有效应力原理和三维渗流理论，建立了太谷东站地下水三维流和土体一维垂向变形耦合模型，可以刻画出研究区实际水文地质特征。

(2)利用水土部分耦合模型对研究区地面沉降进行研究，模型在现了地下水水位与地面沉降的动态特征与演化趋势，并预测了2025a 8月研究区的沉降。太谷东站在预测期累积地面沉降量最高达到51.2mm，沉降漏斗面积可达10.64km2。GMS-SUB模块计算地面沉降参数简明，速度快，实用性高，计算结果能反映地面沉降动态过程，为地面沉降防治工作提供科学依据。

(3)基于GMS-SUB模块的太谷东站地面沉降模型主要考虑的是一维垂向沉降，而实际的地面沉降往往也伴有水平向的变形，水土模型数值模拟还需要进一步完善。