地面沉降是孔隙水承担的孔隙水压力和土骨架承担的有效应力发生变化的结果。产生地面沉降的原因包含两个部分自然因素和人为因素。开采地下水是产生地面沉降的人为因素，松散未固结土层的存在是地面沉降产生的自然因素。因为地面沉降机理的复杂性，目前学术界尚未有一个完整的定论。但是鉴于地面沉降范围广，含水层因开采地下水引起的压缩量相对于含水层厚度相比太小，可假定土颗粒间无胶结，即土颗粒间作用力忽略不计汇，，运用一维太沙基固结理论中的有效应力原理近似地解释开采地下水所引起的地面沉降演变过程，这是目前学术界普遍应用的一种理论。
土体的有效应力原理
根据有效应力原理，饱和土体中任何一点的法向总应力不随着时间变化。开采地下水之前，含水层上覆荷载由含水层骨架及水体共同承担达到平衡，随着地下水开采的增加，孔隙水压力不断减小，而上覆荷载总应力没有改变，含水层中有效应力必然会增加，即原先由孔隙水体承担的一部分荷载转而由土体骨架承担，引起土体骨架所受附加应力增大导致颗粒间孔隙缩小，含水层在宏观上表现为压缩土颗粒压缩忽略不计。当抽取承压含水层地下水引起水位下降时，含水层的相邻弱透水层中孔隙水压力也随之减小，弱透水层也因有效应力增大而压 缩，所以地面沉降总量包括含水层压缩量与弱透水层压缩量之和，。下面对含水层组开采地下水后的应力变化进行分析。
含水层应力变化的基本规律
假定含水层因开采地下水孔隙水压力下降了△H，根据有效应力原理，含水层有效应力就增加了△σ
研究区域的土层概化
在进行地面沉降计算之前，需要对地下各土层运用水文地质中有关含水层分类的概念进行简化。含水层一般是指位于地下水面之下的透水层，且经常为地下水所饱和的土层。土层按其透水性的强弱和受压情况可分为潜水含水层、弱透水层、承压含水层三种。 潜水含水层：直接与包气带相接具有自由表面的含水层。
弱透水层： 在自然条件下几乎不能透水或只能透过数量很少水的土层称为弱透水层。在某些假定理想条件下也可称为隔水层，但这种情形在自然界中是不存在的，所以此类土层称为弱透水层较为合适。
承压含水层：含水层上面存在稳定弱透水层，含水层中的水具有承压性。开采后承压水水头出露于其上的弱透水层标高之上。
地面沉降的数学模型
本文研究对象为大面积区域的地面沉降，暂不考虑开采地下水所引起的土体水平运动，所以采用三维土体渗流基本微分方程与一维固结方程相结合的完全耦合模型模拟开采地下水引起的地面沉降运动基本上是可行的。
（1）土体渗流基本微分方程
x(kxhx)
y(kyhy)
z(kzhz)WSsht
式中：x，y，z——为笛卡尔坐标系中的三个方向并与渗透系数主方向相一致，单位（L）
-1 kx ky kz——分别为与坐标轴方向一致的主渗透张量，单位（LT）
H——含水层某计算单元的总水头，单位（L）
W——为计算单元的源汇项。定义为单位体积含水层中，单位时间内产生或消耗的水量，W为正值时表示为源；W为负值时表示为汇，单位（L-1）
Ss——含水层的单位释水（储存）系数，单位（L-1）
t——时间计算长度，单位（T）
（2）土体一维固结方程
一维固结方程基于Terzaghi—Jacob理论形式，在弹性变形阶段一般假设含水层（弱透水层和含水层）固结的压缩变形与有效应力增量成正比：
b
wSskeb0
Δb—含水层的压缩量，正值代表压缩;负值代表膨胀，单位（L）
—有效应力增量，正值代表增大，负值代表减小，单位（FL-2）
w—水的容重（KN/L3）
b0——含水层的初始厚度（L）
Sske——含水层的弹性释水系数，单位（L）
应用有效应力原理对于沉降变形的饱和含水层其总应力是常数，有效应力增量,等于水的容重与水头降深的乘积
-1
所以式可以表示为含水层固结变形的压缩量与含水层水头降深成正比：
适用于计算含水层组的一维垂向弹性变形。大量固结实验数据和观测数据表明，当含水层压缩变形所受到的有效应力大于前期固结应力，含水层变形为非弹性变形且是不可恢复的。在非弹性阶段，有效应力的增量幅度大于弹性阶段，当有效应力减小时，含水层变形又会出现弹性性质。如果有效应力再次大于新的前期固结应力，含水层变形又再次表现为非弹性变形。因此，Helm（1975）在发展一维非弹性固结模型时认为含水层组变形量仍然与有效应力增量成正比关系，其方程式与一维弹性固结方程相类似：
其中 ∆b'—含水层在非弹性阶段中的压缩量，单位（L）
Sskv—含水层非弹性释水系数，单位（L-1）
运用有效应力原理相应地可以得到含水层的非弹性固结压缩量与含水层水头降深成正比：
以上各式表达的是弹性和非弹性沉降，可统一为如下一个表达式来表示：
其中 b（t）—单位时间内，含水层的压缩量
Φ（t）—单位时间内，含水层的有效应力增量与水的容重之间的比值
将地下水运动与土体变形方程的耦合
式中
X，y，z——为笛卡儿坐标系中的三个方向并与渗透系数主方向相一致，单位（L） kx，ky，kz——分别为与坐标轴方向一致的主渗透系数张量，单位（LT-1） H——含水层某计算单元的总水头，单位（L）
W——为计算单元的汇源项。定义为单位体积含水层中，单位时间内产生或消耗水量，为正值时表示为源为负值时表示为汇，单位（L-1）
Ss——含水层的单位释水储存系数，单位（L-1）
t——时间计算长度，单位（T）
b（t）——单位时间内，含水层的压缩量
φ（t）——单位时间内，含水层的有效应力增量与水的容重之间的比值
应用此模型对美国的波士顿、德克萨斯州、加利福尼亚地区的地面沉降进行计算与预测，取得较为精确的结果，证明了该模型对于处理大面积区域的一结沉降是可行的。而一维非弹性固结方程的可靠性在固结试验中得到近似地证明。
模型的应用：
（1）单位释水系数一般定义为测压水头上升下降一个单位时，从单位体积含水层柱体所储存释放的体积水量。从式可以看出地面沉降量的大小取决于含水层的释水储水系数、含水层厚度、水头降深。含水层厚度可以通过地质资料得到，水头降深由三维水流方程得到，所以释水系数的确定至关重要。在地面沉降计算中常常将释水系数看作常数或看作相关变量的函数。
如果认为释水系数为常数可以利用大量沉降量与水头降深的实测数据运用上述模型反推出含水层的弹性释水系数和非弹性释水系数。在leake《区域地下水流模型中的土层垂向压缩模拟》中提出单位释水系数是孔隙比及有效应力的函数。
其中 e0——初始孔隙比 C—e-logσ‘曲线上确定的压缩指数和回弹指数，当有效
应力小于前期固结应力时为压缩指数Cc，当有效应力大于前期固结应力时为回弹指数Cs
（2）简化方法求释水系数的方法
去上述方法反演比较繁琐，需要长期的水位观测资料和沉降观测资料，在很多地区都难以满足这一条件，且上述模型只适用于般应用于一维准弹性固结模型，针对大面积区域的地面沉降，能够模拟计算出抽水引起的土层主固结基本上是满足沉降预测要求的，因此运用一维弹性固结方程是合理的。
根据土力学知识结合地下水动力学有
Ssw(1-2)(1)
E(1)
其中 μ——土的泊松比
E——土的变形模量
在一般情况下，因弱透水层释水产生的压缩量往往大于成员含水层，因为土体压缩量大其压缩系数就小，承压含水层多为砂砾颗粒所组成， 其压缩模量往往大于由粘土颗粒组成的弱透水层的压缩模量。
（3）软件的使用
对以上建立的三维水流方程和一维固结方程相结合的完全耦合模型进行有限差分法求解即可用于地面沉降预测工作实际，本文利用地下水流动模拟软件ProeessingModflow (PMWIN)求解模型。地下水流动模拟软件proeessingModflow(PMWIN)是由美国地质调查局于二十世纪80年代开发出来的一套专门用于地下水三维流动数值模拟的软件。该软件由许多标准模块组成，每一个模块解决一个特定问题。由于各模块的代码完全公开，所以可以对某一模块进行改进或编写一个新的模块加入其中。对于开采地下水引起的地面沉降问题主要调用Intethed一stora模块包，它是由Leake和 prudie(1991)编写加入的，能够模拟计算含水层组的弹性变形和非弹性变形。
下面
自编算例研究区域假设为5km*5km，含水层组有三层分别为潜水含水层、弱透水层、承压含水层:水位埋深为10m:含水层组初始水头相一致为14m;区域内有一口抽水深井一直打到承压含水层底板，出水量初始为5000(单位:吨/天)，以后每年逐步加大开采量，总共计算5年最后一年采量为9000t/d边界条件为:
东西向只有第一层有水头补给其余各层没有，南北向三层均有山脉割断，与外界没有水头联系;观测连续抽水五年后的沉降。根据参考文献，各含水层的水地质参数见下表
地面沉降是孔隙水承担的孔隙水压力和土骨架承担的有效应力发生变化的结果。产生地面沉降的原因包含两个部分自然因素和人为因素。开采地下水是产生地面沉降的人为因素，松散未固结土层的存在是地面沉降产生的自然因素。因为地面沉降机理的复杂性，目前学术界尚未有一个完整的定论。但是鉴于地面沉降范围广，含水层因开采地下水引起的压缩量相对于含水层厚度相比太小，可假定土颗粒间无胶结，即土颗粒间作用力忽略不计汇，，运用一维太沙基固结理论中的有效应力原理近似地解释开采地下水所引起的地面沉降演变过程，这是目前学术界普遍应用的一种理论。
土体的有效应力原理
根据有效应力原理，饱和土体中任何一点的法向总应力不随着时间变化。开采地下水之前，含水层上覆荷载由含水层骨架及水体共同承担达到平衡，随着地下水开采的增加，孔隙水压力不断减小，而上覆荷载总应力没有改变，含水层中有效应力必然会增加，即原先由孔隙水体承担的一部分荷载转而由土体骨架承担，引起土体骨架所受附加应力增大导致颗粒间孔隙缩小，含水层在宏观上表现为压缩土颗粒压缩忽略不计。当抽取承压含水层地下水引起水位下降时，含水层的相邻弱透水层中孔隙水压力也随之减小，弱透水层也因有效应力增大而压 缩，所以地面沉降总量包括含水层压缩量与弱透水层压缩量之和，。下面对含水层组开采地下水后的应力变化进行分析。
含水层应力变化的基本规律
假定含水层因开采地下水孔隙水压力下降了△H，根据有效应力原理，含水层有效应力就增加了△σ
研究区域的土层概化
在进行地面沉降计算之前，需要对地下各土层运用水文地质中有关含水层分类的概念进行简化。含水层一般是指位于地下水面之下的透水层，且经常为地下水所饱和的土层。土层按其透水性的强弱和受压情况可分为潜水含水层、弱透水层、承压含水层三种。 潜水含水层：直接与包气带相接具有自由表面的含水层。
弱透水层： 在自然条件下几乎不能透水或只能透过数量很少水的土层称为弱透水层。在某些假定理想条件下也可称为隔水层，但这种情形在自然界中是不存在的，所以此类土层称为弱透水层较为合适。
承压含水层：含水层上面存在稳定弱透水层，含水层中的水具有承压性。开采后承压水水头出露于其上的弱透水层标高之上。
地面沉降的数学模型
本文研究对象为大面积区域的地面沉降，暂不考虑开采地下水所引起的土体水平运动，所以采用三维土体渗流基本微分方程与一维固结方程相结合的完全耦合模型模拟开采地下水引起的地面沉降运动基本上是可行的。
（1）土体渗流基本微分方程
x(kxhx)
y(kyhy)
z(kzhz)WSsht
式中：x，y，z——为笛卡尔坐标系中的三个方向并与渗透系数主方向相一致，单位（L）
-1 kx ky kz——分别为与坐标轴方向一致的主渗透张量，单位（LT）
H——含水层某计算单元的总水头，单位（L）
W——为计算单元的源汇项。定义为单位体积含水层中，单位时间内产生或消耗的水量，W为正值时表示为源；W为负值时表示为汇，单位（L-1）
Ss——含水层的单位释水（储存）系数，单位（L-1）
t——时间计算长度，单位（T）
（2）土体一维固结方程
一维固结方程基于Terzaghi—Jacob理论形式，在弹性变形阶段一般假设含水层（弱透水层和含水层）固结的压缩变形与有效应力增量成正比：
b
wSskeb0
Δb—含水层的压缩量，正值代表压缩;负值代表膨胀，单位（L）
—有效应力增量，正值代表增大，负值代表减小，单位（FL-2）
w—水的容重（KN/L3）
b0——含水层的初始厚度（L）
Sske——含水层的弹性释水系数，单位（L）
应用有效应力原理对于沉降变形的饱和含水层其总应力是常数，有效应力增量,等于水的容重与水头降深的乘积
-1
所以式可以表示为含水层固结变形的压缩量与含水层水头降深成正比：
适用于计算含水层组的一维垂向弹性变形。大量固结实验数据和观测数据表明，当含水层压缩变形所受到的有效应力大于前期固结应力，含水层变形为非弹性变形且是不可恢复的。在非弹性阶段，有效应力的增量幅度大于弹性阶段，当有效应力减小时，含水层变形又会出现弹性性质。如果有效应力再次大于新的前期固结应力，含水层变形又再次表现为非弹性变形。因此，Helm（1975）在发展一维非弹性固结模型时认为含水层组变形量仍然与有效应力增量成正比关系，其方程式与一维弹性固结方程相类似：
其中 ∆b'—含水层在非弹性阶段中的压缩量，单位（L）
Sskv—含水层非弹性释水系数，单位（L-1）
运用有效应力原理相应地可以得到含水层的非弹性固结压缩量与含水层水头降深成正比：
以上各式表达的是弹性和非弹性沉降，可统一为如下一个表达式来表示：
其中 b（t）—单位时间内，含水层的压缩量
Φ（t）—单位时间内，含水层的有效应力增量与水的容重之间的比值
将地下水运动与土体变形方程的耦合
式中
X，y，z——为笛卡儿坐标系中的三个方向并与渗透系数主方向相一致，单位（L） kx，ky，kz——分别为与坐标轴方向一致的主渗透系数张量，单位（LT-1） H——含水层某计算单元的总水头，单位（L）
W——为计算单元的汇源项。定义为单位体积含水层中，单位时间内产生或消耗水量，为正值时表示为源为负值时表示为汇，单位（L-1）
Ss——含水层的单位释水储存系数，单位（L-1）
t——时间计算长度，单位（T）
b（t）——单位时间内，含水层的压缩量
φ（t）——单位时间内，含水层的有效应力增量与水的容重之间的比值
应用此模型对美国的波士顿、德克萨斯州、加利福尼亚地区的地面沉降进行计算与预测，取得较为精确的结果，证明了该模型对于处理大面积区域的一结沉降是可行的。而一维非弹性固结方程的可靠性在固结试验中得到近似地证明。
模型的应用：
（1）单位释水系数一般定义为测压水头上升下降一个单位时，从单位体积含水层柱体所储存释放的体积水量。从式可以看出地面沉降量的大小取决于含水层的释水储水系数、含水层厚度、水头降深。含水层厚度可以通过地质资料得到，水头降深由三维水流方程得到，所以释水系数的确定至关重要。在地面沉降计算中常常将释水系数看作常数或看作相关变量的函数。
如果认为释水系数为常数可以利用大量沉降量与水头降深的实测数据运用上述模型反推出含水层的弹性释水系数和非弹性释水系数。在leake《区域地下水流模型中的土层垂向压缩模拟》中提出单位释水系数是孔隙比及有效应力的函数。
其中 e0——初始孔隙比 C—e-logσ‘曲线上确定的压缩指数和回弹指数，当有效
应力小于前期固结应力时为压缩指数Cc，当有效应力大于前期固结应力时为回弹指数Cs
（2）简化方法求释水系数的方法
去上述方法反演比较繁琐，需要长期的水位观测资料和沉降观测资料，在很多地区都难以满足这一条件，且上述模型只适用于般应用于一维准弹性固结模型，针对大面积区域的地面沉降，能够模拟计算出抽水引起的土层主固结基本上是满足沉降预测要求的，因此运用一维弹性固结方程是合理的。
根据土力学知识结合地下水动力学有
Ssw(1-2)(1)
E(1)
其中 μ——土的泊松比
E——土的变形模量
在一般情况下，因弱透水层释水产生的压缩量往往大于成员含水层，因为土体压缩量大其压缩系数就小，承压含水层多为砂砾颗粒所组成， 其压缩模量往往大于由粘土颗粒组成的弱透水层的压缩模量。
（3）软件的使用
对以上建立的三维水流方程和一维固结方程相结合的完全耦合模型进行有限差分法求解即可用于地面沉降预测工作实际，本文利用地下水流动模拟软件ProeessingModflow (PMWIN)求解模型。地下水流动模拟软件proeessingModflow(PMWIN)是由美国地质调查局于二十世纪80年代开发出来的一套专门用于地下水三维流动数值模拟的软件。该软件由许多标准模块组成，每一个模块解决一个特定问题。由于各模块的代码完全公开，所以可以对某一模块进行改进或编写一个新的模块加入其中。对于开采地下水引起的地面沉降问题主要调用Intethed一stora模块包，它是由Leake和 prudie(1991)编写加入的，能够模拟计算含水层组的弹性变形和非弹性变形。
下面
自编算例研究区域假设为5km*5km，含水层组有三层分别为潜水含水层、弱透水层、承压含水层:水位埋深为10m:含水层组初始水头相一致为14m;区域内有一口抽水深井一直打到承压含水层底板，出水量初始为5000(单位:吨/天)，以后每年逐步加大开采量，总共计算5年最后一年采量为9000t/d边界条件为:
东西向只有第一层有水头补给其余各层没有，南北向三层均有山脉割断，与外界没有水头联系;观测连续抽水五年后的沉降。根据参考文献，各含水层的水地质参数见下表