我国地面沉降模拟现状及需要解决的问题.doc

我国地面沉降模拟现状及需要解决的问题*薛禹群(南京大学 地球科学系 南京 210093)摘要：我国现有模型有下列特点：水流模型都是准三维模型，模型中的参数都是常数，沉降模型采用线弹性模型，水流模型和沉降模型的耦合分二步进行。因此，模型还存在难以精确刻画由弱透水层释水引起的地面沉降，不符合准三维模型的应用条件。水流方程应是一个变系数的动态水流方程，渗透系数、贮水率不应视为常数。土层压縮非线性，有蠕变、塑性变形存在，沉降模型采用线弹性模型不适合，对土层厚度变化大的地区沉降模型必须是三维的，迟后现象必须克服。如何描述引起地裂缝的水平位移，分两步进行的耦合不合适，要做到水流模型与沉降模型的真正耦合等八个方面的步骤。本文提出分三步进行的解决方案：首先采用三维水流模型，以便不仅给弱透水层建立方程，还能考虑它的各向异性，参数值能随着土层的变形沉降而不断改变。沉降模型为能考虑蠕变、塑性变形的非线性模型；水流模型和沉降模型要做到能真正耦合。其次，解决迟后问题。最后，建立三维沉降模型，考虑水平位移和地裂缝问题。主题词：地面沉降; 地裂缝; 地面沉降模型; 塑性变形; 蠕变中图分类号：P641.2; 文献标识码：A 文章编号：1000-3665（2003）05-0000-05 1 现状我国地面沉降模拟工作的开展已有10多年的历史。早期上海地质工程勘察院与比利时合作开展地面沉降研究，1989年完成我国第一个地面沉降模型1。1995年地矿部水文地质工程地质研究所等在此基础上，又提出“地下水水量水位沉降联合数学模型”2。上世纪末上海市地质调查院李勤奋等完成了第三个模型3，目前他们仍在继续研究如何考虑非线性变形的问题。1995年天津地质环境监测站与英国合作，完成了天津市地面沉降模型4。虽然我国发生不同程度地面沉降的城市已达四十多个，但进行过实际地面沉降模拟的城市迄今仅此二个。陈崇希等（1999）曾对苏州的地面沉降进行了有益的探讨，但由于该地区缺少地下水水位和地面沉降的实际观测资料，难以真正实施5,6。上海、天津的这些模型有下列共同特点：（1）水流模型都是准三维模型，只对含水层建立相应的水流方程，研究它的水位变化。对含水层之间的弱透水层或“隔水层”一般不建立方程，刻画它的水位变化，弱透水层的释水量作源汇项处理或由估算得。如上海的模型（1999）1为式中：各个含水层的水头(m)；越流量（m3/d）；土层变形释水量或回弹吸水量(m3/d)；-收稿日期：2003-04-08；修订日期：2003-05-09基金项目：国家自然科学基金（4017082），博士点基金（20010284002）资助作者简介： 薛禹群（1931-）， 男， 中国科学院院士、教授、从事水文地质科研与教学 抽水量(m3/d)；灌水量(m3/d)；S贮（释）水系数；T导水系数(m2/d)。近年来李勤奋对上述模型作了改进，为弱透水层建立了垂向一维水流模型，可以说这是迄今我国考虑因素最多的一个地面沉降水流模型3；（2）水流模型中的参数（渗透系数或导水系数、贮水系数等）都是常数；（3）沉降模型是垂向一维模型，采用根据Terzaghi理论建立的线弹性模型；（4）水流模型和沉降模型的耦合分二步进行，先由水流模型求出水位（水压），作为沉降模型的边界条件，再进行计算。2 问题这样的模型虽然在前一阶段的地面沉降预测预报中起了很大作用，但毕竟有它的很多不足，既无法满足描述复杂地质、水文地质条件的要求，又无法满足日益增长的精确预警预报的需要。主要原因和存在的问题有718：（1）不对弱透水层建立水流方程，模拟它的水位变化，显然难以精确刻画由它释水引起的地面沉降或者即使建立一个垂向一维水流模型，也无法满足下面将要提到的实际问题是各向异性的情况。（2）准三维模型要求弱透水层是均质、各向同性的，此时误差很小（5%。但上海、天津以及我国很多地区并不符合这种情况，冲积成因的弱透水层往往是各向异性的，一层粉砂，一层粘土或亚粘土，单层韵律厚度只有12或23mm，如此重复，总厚度可达数十米或更多。这种岩层水平方向的渗透系数远大于垂直方向的渗透系数，因此,不符合准三维模型各向同性的应用条件，也不满足垂向一维水流模型可以忽略水平方向流动的要求。勉强使用误差显然就会很大。诚如上世纪90年代初，在我国开展的一场有关“拟合”的大讨论中，笔者曾指出19，计算所得结果不好，不在公式或数值法本身，因为每个公式或每种数值方法都有它的适用条件，现在的问题出在没有正确应用公式和数值法，违反了它的使用条件。只有改用真正的三维模型才能解决问题。（3）随着固结沉降，土层的孔隙度必然不断减小，渗透系数也会随之减小。因此，水流方程中的参数如渗透系数应该不仅是x ,y ,z 的函数，还应是随时间在不断减小的孔隙度的函数。也就是说，这种情况下的水流方程应该是一个变系数的动态水流方程。渗透系数、贮水率都不应视为常数。（4）土层的压縮是非线性的，有蠕变、塑性变形存在。不仅粘土、亚粘土的变形表现为非线性，组成上海地区含水层的砂土取样后，经我们做试验（三轴试验条件下的应力应变关系；单向压缩试验条件下荷载与应变的关系；分级加载时每一级荷载作用下变形与时间的关系），均证实含水砂层的变形具有明显的非线性和蠕变性。分层标所测得的砂层实际变形与水位变化关系，也显示砂层变形表现出较大的非线性和与时间的相关性，水位上升导致的回弹主要表现为弹形变形，但变形量非常小（如上海劳动公园分层标1986年以来所测第四承压含水砂层变形的观测资料显示1991年以前水位上升时，砂层略有回弹，但1991年以后，特别是1993年以后水位上升时，砂层变形仍在继续增加，土体变形表现出较大的非线性和与时间的相关性）。因此，沉降模型采用线弹性模型显然是不适合的。在抽、灌水的反复作用下，土层内的孔隙水压力在上升下降的循环往复中有效应力也会相应地变化，而且变化速率较慢，和一般荷载情况下的变形不同。由于土层的塑性变形和蠕变变形的存在，使得每循环一次土层的变形特性都有变化，这些变化与循环次数之间究竟存在什么关系现在还研究得不够, 上述反复荷载作用和迟后有没有关系也有待研究。顺便说一句，各地岩性、地层结构、沉积相差别很大，导致土层的变形特性也有很大差别，必须分别研究。一些地区有很厚的粘土层，深埋于地面以下数百米、甚至近千米，为了做实验，不仅取样困难，更困难的是这些土样拿到地面后，压力改变了，失去了原有的应力状态，如何保持土样在近千米地下原有的状态也非易事。（5）土层变形模型采用一维模型对土层厚度变化小的区域（如上海）研究区域性沉降是合适的，但对土层厚度变化大的地区未必合适，因为那儿（如无锡的一些地区）的变形表现出明显的差异性沉降，变形应该是三维的。（6）水位下降与沉降变形并不同步，有明显的迟后存在，像天津即使水位恢复后，仍有明显沉降。迟后是一个比较复杂的问题，它主要与蠕变有关，但可能更复杂。如何用数学表达式来描述这种迟后现象是模拟这些地区地面沉降必须克服的一个难题。在此以前，首先得弄清楚迟后的机理。目前一些表明存在迟后现象的地区，实际上对机理并没有完全搞清楚，实测资料中的所谓沉降部分，事实上包含了一部分严格说来并不属于迟后的固结沉降。因为一般所谓的含水层水位恢复后仍继续沉降，这里面可能包含有两种情况：首先，含水砂层的蠕变导致含水层水位恢复后，继续发生沉降变形；其次，由于抽水，含水层水位下降，导致水位较高的相邻弱透水层因越流而水位下降，出现固结沉降。当含水层水位恢复时，只要水位仍低于相邻弱透水层的水位，有越流就会有弱透水层的固结沉降。有些迟后有很大一部分实际上指的是这种情况。两者水位差消失，没有越流后，在蠕变影响下，弱透水层的沉降仍将继续很长一段时间。这三种情况在地面上的表现是相似的，都是沉降，但从严格意义上说，弱透水层的固结沉降不宜算作迟后，只有与蠕变有关的变形现象才是迟后。另一方面，即使是含水砂层中也往往含有一些粘土、亚粘土凸镜体或夹层，一般不可能为这些凸镜体或夹层单独建立分层标，所以即使是含水层中的分层标所反映的沉降值中很可能包含了一部分由于这些粘土、亚粘土凸镜体或夹层失水而造成的固结沉降。至于没有建立分层标的地区，在地面测得的沉降总量所反映的所谓迟后就更复杂了。（7）地裂缝问题，如何描述引起地裂缝的水平位移。（8）目前分两步走的耦合，不仅由于沉降后计算水流运动的参数变了，再用那个水流模型计算得出的水头并不反映真正的水头；而且应该明白，一旦算出水头事实上就意味着越流量和该层的沉降量已经定了。两步走这种耦合办法中再用算得的水头来计算沉降量还有什么意义呢？如何做到水流模型与沉降模型的真正耦合，随着土层的压缩沉降、孔隙度不断减小，渗透系数、贮水率也不断自动改变，从而改变水流方程的模拟结果是下一步要研究解决的。3 解决方案这些问题的解决有一定难度，有些问题国外也没有很好解决。因此，笔者认为同步解决可能有困难。分步走，考虑轻重缓急，由易而难逐步解决可能更稳妥，也更现实些。由于我国发生地面沉降的城市和地区已有四十多个，累计沉降量相当可观，给国民经济造成了重大损失，国家自然科学基金会、中国地质调查局等有关部门在制定资助计划或工作计划时宜给予优先考虑，以使这方面的研究有较快的发展，适合国民经济发展的需要。下列研究顺序仅供讨论，并请国内同人指正。（1）首先采用三维水流模型，以便不仅给弱透水层建立方程，还能考虑它的各向异性。要求所有含水层和弱透水层的水流方程都是变系数方程，参数值能随着土层的变形沉降而不断改变。这里隐含着一个问题，即孔隙比或孔隙度和渗透系数的关系问题需要解决；其次，沉降模型仍先采用一维垂向模型，但不是线弹性模型，而是能考虑蠕变、塑性变形的非线性模型；最后，水流模型和沉降模型要做到能真正耦合。要实现这一步，需要有配套的弱透水层的孔隙水压力观测资料，以便检验模型。但目前由于多种原因，特别是观测装置很容易失效，国内至今还没有这方面完整的资料，需要地质调查研究部门设法加以解决。（2）解决迟后问题，首先要查明一个地区迟后的机理，在此基础上建立不仅考虑土层非线性变形，还考虑变形迟后的一维垂向沉降模型。还有在上述反复荷载作用下土的变形问题。显然，上述两类情况将给模拟带来很大困难。正是由于这种复杂性，我们把它单独列为一步，但和前一步有密不可分的关系。（3）建立三维沉降模型，同时考虑水平位移的计算和地裂缝问题。还有一个问题在进行三维模拟时必须考虑的。由于它不属于水文地质工程地质问题，而是数值方法问题，所以前面没有列入。以上海为例，面积达6400km2，目前使用的模型在平面上剖分了约1500个单元，还觉得单元太大，难以刻画各类变化。所以上海从事这方面工作的有一个想法觉得剖分成大约5000个单元为好。上海有5个含水层，显然至少还有5个弱透水层（事实上由于岩性的变化，还不止5个）需要考虑。如果采用每层5000个单元，则总单元数将达到50000个以上，如此众多单元的三维模型如果还沿用已有的有限元法在一般计算机上显然是有困难的。即使改为15000个单元（每层1500个），也是有不少困难的。出路就是突破有限元法原有的思维，研究新的算法。在同一单元内就能考虑岩性的变化，如能实现这点，就能大大减少单元总数了。就也是需要单独研究解决的。上述步骤只是一个初步想法，有兴趣的读者完全可以不拘泥于此步骤，可以合并，也可以多分几步做。由于问题的复杂性，单独一个部门、一个学科来解决可能都有困难，地质调查研究部门、教学单位、科研部门密切合作以及水文地质、工程地质、土力学、数学等多个学科的真正交叉是必要的。 “第七届国际地面沉降学术讨论会”将于2005年10月在中国上海举行，模拟是地面沉降研究的一个极其重要部分，在这方面取得一些重要成果是必要的。联合国教科文组织地面沉降工作组第一次把会议地点选在发展中国家也希望我国学者能提交一些出色的成果。参考文献1 上海地质环境监测站，葡萄牙国家科技部. 上海市区地面沉降的数学模型研究R. 2 朱锡冰,佴磊,奚建国,等.上海浦东新区“地下水水量水位沉降联合数学模型”应用R.1995.3 李勤奋,方正,王寒梅.上海市地下水可开采量模型计算及预测J.上海地质.2000，（2）:36-43.4 吴铁钧、金东锡，天津地面沉降防治措施及效果J.中国地质灾害与防治学报，1998,9(2):5 陈崇希,裴顺平. 地下水开采-地面沉降数值模拟及防治对策研究以江苏省苏州市为例M. 武汉:中国地质大学出版社,2001.6 陈崇希，地下水开采-地面沉降模型研究J. 水文地质工程地质，2001，28(2)：5-8. 7R Bravo, J R.Rogers, et al.休斯顿地区地下水流和地面沉降的新三维有限差分模型A.第四届地面沉降国际讨论会论文集C. 北京：地质出版社，1993.8 A Rivera, E.Ledovx, G de Marsily. 墨西哥城含水层弱透水层系统的地下水径流和地面沉降总量的非线性模拟A.第四届地面沉降国际讨论会论文集C北京：地质出版社，1993,11-21.9 G Gambolati, G Ricceri,et al.由水、气开采引起的拉温万纳地面沉降的数值分析A.第四届地面沉降国际讨论会论文集C北京：地质出版社，1993,34-40.10 K Daito, M Mizuno,et al. 日本大鳄平原为防止地面沉降而控制地下水抽水量A. 第四届地面沉降国际讨论会论文集C. 北京：地质出版社，1993.11 S A Leake. 区域地下水流模型中的土垂向压缩模拟A.第四届地面沉降国际讨论会论文集C. 北京：地质出版社，1993.12 X Y Gu, D N Xu And W Deng. A computing model based cyclic consolidation testsA. Land Subsidence: Proceedings of 5th International Symposium on Land SubsidenceC. 1995: 295-303 .13 A Rossi, C Calore and U Pizzi, Land subsidence of pisa plain, Italy: experimental results and preliminary modeling studyA. Land Subsidence: Proceedings of 6th International Symposium on Land Subsidence, ItalyC. 2000,1:91-104.14 G.de Lange, G B M Brand, A J Bruijne et al., Modeling regional and local surface subsidence due to compaction of unconsolidation sedimentsA. Land Subsidence: Proceedings of 6th International Symposium on Land Subsidence, ItalyC. 2000, 2:289-304.15 M Sneed, M T Pavelko and D L Galloway. Modeling residual aquifer-system compaction: constraining the vertical hydraulic diffusivity of thick aquifersA. Land Subsidence (Proceedings of 6th International Symposium on Land Subsidence, ItalyC. 2000,(2)：343-355.16 Gambolati , P Teatini, L Tomasi. Stress-strain analysis in productive gas/oil reservoirsJ. International Journal for Numerical and Analytical Methods in Geomechanics,1999,23(13)：1495-1519.17 A Ortega-Guerrero, D L Rudolph, J A Cherry. Analysis of long-term land subsidence near Mexico City: field investigations and predictive modelingJ. Water Resources Research. 1999, 35 (11)：3327-3341.18 K J Larson, H Basagaoglu, M A Marino. Prediction of optimal safe ground water yield and land subsidence in the Los Banos-Kettleman City area, California, using a calibrated numerical simulation modelJ. Journal of Hydrology.2001,242 (1-2):79-102.19薛禹群、吴吉春. 数值模拟是反映客观规律和定量评价的重要手段J.水文地质工程地质.1992，（2）：2-4.Present status of modeling land subsidence in China and problems to be solvedXUE Yu-qun (Department of Earth Sciences, Nanjing University, Nanjing, 210093 China)Abstract: The characteristics of present land subsidence models in China are as follows: (1)flow models are quasi-three dimensional models; (2) parameters in the models are constants; (3)subsidence models are linear elastic models; (4) couple of flow model and subsidence model has two steps. Hence, there are eight disadvantages for present models. (1) They can not precisely describe land subsidence caused by drainage of aquitards; (2) They do not agree with the conditions of quasi-three dimensional model; (3) Flow equation is a dynamic equation in which the coefficients vary with time. So hydraulic conductivity and specific storativity should not be constants; (4) It is improper to use linear elastic model for a subsidence model, because soil shows nonlinear stress-strain relation, and there exits plastic deformation and creep;(5) Three dimensional subsidence model is necessary for the area with large variation of thickness of the soil layer (Quaternary sediments); (6) The