鹏山水源地开采条件下的含水层固结模型.pdf

鹏山水源地开采条件下的含水层固结模型鹏山水源地开采条件下的含水层固结模型 摘摘 要：要：以鹏山水源地地下水系统作为研究对象，针对研究区地下水开采引发的地面沉降问题，应用土力学和地下水动力学原理，建立了该水源地开采条件下的含水层固结模型。为保证模型的准确性，从渗流空间场和地下水运动过程两个方面，对模型误差积累、误差分布和模型稳定性进行了较全面的评估，使模型具有较高的可信度。关关 键键 词：词：地下水资源；地面沉降；含水层固结；有限单元法；可持续发展 中图分类号：中图分类号：TU 433 文献标识码：文献标识码：A Consolidation model of aquifer under pumping condition at Pengshan water well field ZHANG Li-sheng1,ZHU Guo-rong2（1.Shandong Geotechnical Investigation Corporation,Jinan 250014,China;2.Department of Earth Sciences,Nanjing University,Nanjing 210093,China）Abstract：This study is based on the groundwater system in Pengshan water resources area,aiming at the problem of ground subsidence in study area,applying the theory of soil mechanics and dynamics of groundwater,the aquifer consolidation model under the condition of exploitation in water resources is established.In order to ensure the exactness of this model,making use of seepage spatial field and groundwater movement progress respectively,the author comprehensively evaluates the error cumulating and distribution and the model stability so as to lead to better reliability in the model.Key words：groundwater resources;ground subsidence;aquifer consolidation;finite element method;sustainable development 1 引 言 鹏山水源地位于莱芜市鹏山脚下蟠龙河与辛庄河交汇处，供水能力达莱芜市工业与居民用水的76%左右。自 1999 年底开始，水源地原有部分供水井开始出现水质混浊、涌砂、出水量降低等现象，于是在水源地西侧新凿 Z11水井，并于 2000 年 2 月18 日正式启用，同年 4 月即在百嘴红村的东部一带陆续产生地面沉降、地面塌陷、地裂缝及房屋墙体开裂等地质灾害。经调查，地面变形随水位的升降有所起伏。本文针对研究区地面沉降问题，从地下水资源合理开发利用角度，分析研究地面沉降的机理和固结模式，以探索缓解地面沉降继续恶化的对策，减轻或避免由此产生的重大经济损失，为莱芜市经济可持续发展战略提供决策依据。2 研究区地质构造概况 莱芜盆地在大地构造上受鲁中纬向构造和鲁西旋卷构造的控制，为一中生代形成的断陷盆地，北、东、南三面环山，中部为低缓起伏的平原，西部较为开阔。整个地势由东向西倾斜，北、东、南三面又向盆地中部倾斜(见图 1)。鹏山水源地位于莱芜盆地东南部，地貌上属构造侵蚀、中切割低缓丘陵至中低山区，因受莱芜弧形断裂的控制，呈单面山地形向西南倾斜。与本研究有关的主要构造有：(1)F1断裂 F1断裂是莱芜弧形断裂的东部延伸段，其东北盘为泰山群花岗片麻岩，西南盘为奥陶系、寒武系灰岩。具压扭性，走向325，倾向SW，倾角6070。破碎带具强烈的高岭土化和糜棱岩化，隔水性好。F1断裂从百嘴红小庄中间穿过。更多岩土工程资料 h t t p:/w w w.d o c i n.c o m/7 5 1 7 4 1 9 5 2 图图 1 莱芜盆地构造略图莱芜盆地构造略图 Fig.1 The structure outline of Laiwu basin (2)F3，F5断裂 为两条近乎平行的断裂，两断裂之间为石炭系地层。石炭系地层与断层上下地层呈不整合接触。该断裂走向o310，倾向NE，倾角约o60。F3的东北盘为奥陶系灰岩，F5西南盘为第三系砾岩。(3)F14断裂 位于百嘴红村西南，呈舒缓波状，总体走向o40，倾向SE，倾角约70。断裂东南盘为第四系，西北盘北部为奥陶系灰岩，南部为第三系砾岩，断裂带内岩溶裂隙相对发育，断裂接受西北和北部山区的地下水补给，再通过F14断裂的伴生断层间接补给研究区。(4)岩溶裂隙 由于本区断裂构造异常发育，断裂带内岩石破碎，致使研究区内上部岩溶裂隙极发育，成为地下水赋存和运移的主要场所。而在此范围以下，岩性一般致密完整，岩溶裂隙均不甚发育，连通性较差，可视为相对隔水底板处理。3 因抽水引起的土层固结模型 3.1 水文地质概念模型水文地质概念模型 这是建立数学模型的基础，它必须突出反映实际水文地质条件的基本特征，又充分保留原型的普遍规律。为查明鹏山水源地的水文地质条件，笔者在充分收集与研究前人成果的前提下，开展了大量的基础工作，基本查清了研究区在现状气象条件及水文地质条件下的地下水开采情况及水位变化情况，为建立研究区地下水模型和开展地面沉降预报提供了可靠的最基础的资料。3.2 定解问题定解问题 能否使数学模型合理而准确地反映客观的水文地质条件，将直接关系到模型的精度。因此，在对研究区的水文地质条件进行概化的基础上，深入了解概化模型中的每一个细节，是判明概化模型是否合理的前提。(1)模型范围 根据抽水试验结果和莱芜市对水源地的用水需求规划，将本次模型平面区间圈定在蟠龙河、辛庄河河谷地区 4.379 km2范围内，见图 2。图图 2 研究区地质构造略图研究区地质构造略图 Fig.2 The structure outline of study region 在垂直方向上，根据含水层空间分布特征，将含水层底板确定在绝对高程150 m处。(2)边界条件 研究区的东南边界外和西北边界外的南部为第三系致密状红色砂岩、砂页岩及第四系残坡积地层分布，渗透性能远弱于研究目的层，确定为隔水边界；研究区的西北边界外的北部及东北边界外的中部为太古界地层，渗透性能较差，作为水源地隔水边界处理；研究区的西北边界外的中部为寒武系、奥陶系碳酸盐岩分布所在，岩溶裂隙相对发育，该处在地表上表现为一山坡陡坎，接受西北和北部山区的地下水补给，再通过F14断裂的伴生断层间接补给研究区，本次研究将其处理成第二类补给边界；研究区的东北边界的北部、南部及西南部边界所在为本区的河谷地区，是本区地下水的径流通道。抽水试验的观测结果表明，分布在这两个边界附近的地下水动态观测孔的水位在抽水试验期间基本上没有受到影响，在不增加水源地开采强度的前提下，完全可以将它作为第一类水头边界对待。辛河庄河 蟠 龙 1142(3)含水层结构的概化 从水文地质角度分析，研究区含水系统主要包括第四系松散沉积物孔隙含水层和岩溶裂隙含水层，其主要补给来源为降雨入渗、河流渗漏及边界侧向补给。岩溶裂隙水多年动态变化主要受控于降水量的多寡，而年内动态变化则受降水及工农业开采的综合影响。但由于裂隙、断层的强烈切割，以及“天窗”的大量存在，二者之间的水力联系十分密切，水位动态亦基本呈现同步变化的趋势。因此，模拟时将两个含水层概化为一个统一的含水系统，即接受大气降雨补给的潜水系统。3.3 源汇项处理源汇项处理(1)地下水开采量pmw 现状水文地质条件下，研究区地下水的开采以水源地供水为主，农用水井开采量极小。本次在模拟计算时，农用水井开采量不予考虑。(2)降雨入渗补给量pew 研究区因地处河谷地带，地表全被第四系覆盖，河谷两侧的斜坡地带将降雨表流输送到河谷地区，参与河谷地区的实际降水的入渗补给。但由于本区第四系含水层之上普遍存在着一层稳定分布、厚度在4.05.0 m左右的粉质粘土层，无形中减弱了降雨的入渗补给强度，本次模型建造中直接引用以往的工作成果，取降水入渗补给系数值为0.11。(3)蒸发量ew 潜水蒸发是潜水在土水势作用下运移至包气带并蒸发为水汽的现象。对潜水蒸发影响最主要的因素是潜水的埋藏深度2,3。试验表明，当潜水埋深达到一定深度时，潜水蒸发微弱甚至趋于零，这一深度称为潜水蒸发的极限深度。研究区地下水位埋深一般都已接近或超过蒸发临界值，为简化模型，计算时不考虑地下水的蒸发量。(4)蟠龙河、辛庄河水系的河道渗漏补给rw 蟠龙河与辛庄河是横贯研究区的主要河流，河床深度一般在5.012.0 m左右，河床底部普遍沉积来自河谷两侧山区的粘土质物质，削弱了地表水通过河床渗漏的强度；另外，由于在河流上游地区业已兴修了一系列的中、小型水库（如乔店水库、杨家横水库），加之近几年上游水库的防渗堵漏及农田灌溉截流等，致使河流常年断流，仅在雨季才出现明流。在通常情况下，进入本区河道的水量极小，通过本区地表河流渗漏补给地下水的能力更是微不足道。本文对河道的渗漏补给未做专门研究，根据以往的工作经验4，取渗漏补给系数值为0.014参加模型建造。3.4 数学模型数学模型 根据以上讨论，鹏山水源地的地下水运动可用二维非均质各向异性定解问题来表达。形式如下：tHyHhKyxHhKxyx=+()0,tDyx (1)(),(),(|,00DyxyxHtyxHt=(2)0 ，),(),(11=ttyxHtyxH (3)0 ，),(),cos(),cos(2=+ttyxqynyHhKxnxHhKyx(4)rpepmwww+=(5)式中 H 为研究区域的水头函数(m)；Kx,Ky分别为坐标轴 x，y 方向的渗透系数(m/s)；h 为含水层厚度(m)；为给水度，无量纲；为源汇项，包括开采量 wpm、降雨入渗量 wpe及河道渗漏补给量 wr(m/s)；q 为第二类边界单宽流量(m2/s)；1为研究区第一类边界；2 为研究区第二类边界；D 为研究区的分布范围。研究区的东北边界的北部、南部及西南部边界为本区与区外的地下水通道，根据抽水试验的结果，可用式(3)描述；西北边界的中部则用式(4)定义；研究区的东南边界、西北边界的南部及北部、东北边界的中部可用式(4)表达，但需令等号右边的q(x,y,t)=0。需要说明的是，由于本区碳酸盐岩地下水的运动具有沿裂隙和溶缝运动的现象，导致地下水的运动具有明显的各向异性特征，根据本区构造和水文地质条件：本区的主干构造线走向为320。以F14为例，该断裂两侧的破碎带和断层影响带构成了一个与断裂走向一致的地下水优势流带，而在断裂垂直方向，地下水运动明显受阻，因此造成了在断裂上盘（西北盘）因接受北部山区地下水而在断裂带处沿裂隙上升成泉的现象。为此，将模型的坐标轴沿顺时针方向旋转50，使坐标系 x，y 方向与渗透的主向量方向保持一致。4 土层固结(地面沉降)计算原理 4.1 土层固结计算公式土层固结计算公式 本研究利用KTerzaghi有效应力原理和一维固结理论58，将地下水流动问题与土层因水头下降而1143产生的固结问题依据贮水率的物理含义合二为一，从而把固结问题联系到上述地下水流动问题之中。有必要说明，用来表征含水层储水能力的给水度、贮水系数、贮水率等水文地质参数，在潜水含水层和承压含水层中有着相同的作用实质，而不管含水层性质、以及排水机理与过程如何9。本研究引用“贮水率”这一术语表示两种含水层的储水能力。贮水率（Storativity）Ss定义为1015：含水层的水头在垂向上变化一个单位时，从单位面积的含水层上所释放或储存的水量。可用下式表示：()nSs+=w (6)式中 Ss为贮水率(m-1)；w为水的重度(kNm-3)；n 为孔隙率，een+=1，e为孔隙比，无量纲；为 水的体积压缩系数(水的压缩性很小，本文在计算地面沉降时，未考虑水体本身压缩的影响)，为水的体积压缩模量wE的倒数(Pa-1)；为土的体积压缩系数，为含水层固体骨架（土颗粒）的体积压缩模量sE的倒数(Pa-1)；根据土的压缩系数va的定义，有ea+=1v。基于上述，利用可把土的压缩系数va（Pa-1）与贮水率Ss联系起来，亦即把土层固结（沉降）与地下水流动问题联系起来。()e1vws+=aeS (7)另外，当土体垂向总应力保持不变时，含水层骨架有效应力的增加可用含水层中水头的下降量来表征。即 Hddw=(8)根据体积压缩系数的定义，对于单位面积地层来说，其单位厚度的垂直压缩量等于有效应力的增量与土的体积压缩系数的乘积。因此，地下水开采引发的地面沉降总量w为=2121ddwHHHHzszw (9)式中 H1，H2分别为可压缩土层底板和顶板的标高(m)；z为垂直坐标；s为抽水引起的水头下降值(m)。由此可见，对于因抽取地下水而引起的地面沉降，当水头下降值已知时，可以利用式（9）直接获得同一平面位置的土层固结总量（即地面沉降量）。计算中利用的是水文地质参数贮水率Ss，而不是土力学参数压缩系数av等。4.2 含水层固结公式的适用范围含水层固结公式的适用范围(1)本文推导的含水层固结公式，是根据 K.Terzaghi有效应力原理和含水层固结理论建立的。根据其假定条件，要求含水层具有如下特征：土体是均匀、饱和的理想弹性材料，只发生竖向压缩变形，且承受的总应力不随时间变化；土颗粒与孔隙水均不可压缩；地下水流服从Darcy定律。因此，本文的含水层固结公式只适用于土体的弹性变形阶段，实际含水层与假定条件的差别会造成一定的计算误差。另外，对因抽水造成土颗粒流失而引起的地面塌陷问题亦不可生搬硬套。(2)土体的变形是孔隙水流出、颗粒重新排列、粒间距离缩短及骨架发生错动的结果。无论是孔隙体积的变化，还是颗粒的重新排列都需要一个时间过程，而潜水含水层的重力排水又存在着明显的滞后现象，短时间内土层中的水不能随着地下水位的下降同时排出1013。因此，土层的固结变形也不会瞬时完成。这将使得抽水初期的实际地面沉降量有可能与模型的计算值存在一定的差异。5 模型识别与评价 5.1 观测孔水位拟合观测孔水位拟合 作者将研究区在平面上剖分成221个三角单元，内结点为91个，边界结点44个。研究区单元剖分图见图3。图图 3 研究区单元剖分图研究区单元剖分图 Fig.3 The element subdivision map of study region 利用有限单元法对所建立的数学模型进行求解、识别，得到表1、图4所示8个地下水动态观测孔的水位拟合情况。拟合结果表明：逆问题求解中存在的平均残差平方和E=0.061，与抽水井的最大水位1144 降深相比：平均拟合误差为0.61%。计算水位与实测水位的误差在0.060.39 m之间的拟合点占总数的95.3%；误差大于0.50 m者仅占总数的4.7%；最大误差绝对值为0.82 m（G7观测孔第6时段）。表表1 观测孔的实测与计算水位观测孔的实测与计算水位 (单位单位:m)Table 1 The observation boreholes computed and observed water levels(unit:m)时段 观测孔 水位值/m 1 2 3 4 5 6 7 8 实测 214.21 214.30 214.44 214.47 214.48 214.46 214.47 214.52 G1 计算 214.28 214.21 214.30 214.25 214.23 213.78 214.19 214.17 实测 214.10 213.92 213.75 213.85 213.85 211.91 212.92 213.15 G2 计算 214.18 214.00 213.65 213.78 213.68 212.27 212.78 213.06 实测 213.98 213.48 213.18 213.21 213.19 210.04 211.61 211.88 G3 计算 213.82 213.69 213.43 213.41 212.90 210.58 211.36 212.10 实测 213.59 212.97 212.57 212.46 212.41 210.13 211.24 211.38 G4 计算 213.69 213.08 212.69 212.33 212.26 210.52 211.34 211.53 实测 215.31 214.93 214.50 215.09 215.13 213.09 214.11 214.25 G5 计算 215.50 215.25 214.84 214.82 214.86 213.46 214.34 214.51 实测 216.51 216.50 216.02 216.54 216.42 215.30 215.71 215.79 G6 计算 216.30 216.32 216.15 216.32 216.29 215.62 215.55 215.55 实测 215.96 216.01 216.01 216.11 216.05 216.03 216.01 215.95 G7 计算 215.85 215.84 215.78 215.78 215.73 215.21 215.70 215.62 实测 216.27 216.16 215.71 216.42 216.16 215.15 215.61 215.72 G8 计算 216.35 216.24 215.86 216.33 216.22 215.39 215.68 215.60 图图4 水位水位-时间过程线拟合图时间过程线拟合图 Fig.4 The water table-time hydrograph fitting map 11455.2 模型解精度分析与评价模型解精度分析与评价 本研究采用实测值-计算值关系分析法（简称H-H分析法）对模型解精度的进行分析与评价。在一个以实测水位H为横坐标、以计算水位H为纵坐标的直角坐标系中，若各时刻的实测水位值与计算水位值完全一致，则所有时刻的实测值与计算值将在这个坐标系统中投影成一条通过原点、斜率为1的直线，作者称其为标准关系曲线。事实上，数值方法所得的解为一近似解，该解是模型解簇中的最优解之一，因此，所得的HH 的投影点不可能完全落在标准关系曲线上，但是，由这些点决定的、通过原点的回归曲线（本文称之为实际关系曲线）将与标准关系曲线形成一个夹角，称其为偏离角。当实际关系曲线与横坐标轴的夹角小于45时称为负偏离，反之称为正偏离。显然，当模型不失真时，这些投影点应该集中在标准曲线附近的一个带上，带的宽度越小，说明模型的拟合程度越高(图5)。图图 5 计算水位计算水位H与实测水位与实测水位H关系图关系图 Fig.5 The relation of computed water tables and observed water tables 采用线性回归分析方法对图5中显示的资料做进一步分析，得到计算水位H与实测水位H之间存在以下的关系：HH9999.0=,相关系数 r=0.987 1 (10)从式(10)看，计算水位整体上略低于实测水位，回归曲线与标准曲线的偏离角为20.6，最大偏离值为0.68 m。以上数据表明：若采用所获得的水流模型进行预测，所得到的水头预测值必将落在由式(11)决定的置信区间：()68.06.2045tan=HH (11)式中 H为计算水头；H为实测水头；20.6为拟合结果偏离角。它不仅可以用来检验水头预测结果的可信度，而且也可以用来检验预测模型的稳定性。5.3 地面沉降拟合结果地面沉降拟合结果 在对模型可信度的评价后，即可利用相点的计算水头下降值按含水层固结计算式(9)计算其地面沉降量。表2 列出了部分监测点地面沉降量观测值与计算值结果比较。表表2 各监测点地面沉降量观测值与计算值比较各监测点地面沉降量观测值与计算值比较 Table 2 Computed and observed values comparison for land subsidence of different monitoring points 监测点编号水位降深/m 压缩层厚度/m 实测沉降量/mm 计算沉降量/mm 绝对误差/mm C10 0.20 9.80 0.881 0.980 0.099 C10 0.27 9.30 0.971 1.256 0.285 C11 0.24 8.00 0.862 0.960 0.098 C11 0.27 8.50 0.973 1.148 0.175 C29 0.20 7.50 0.599 0.750 0.151 C29 0.26 6.20 0.615 0.806 0.191 6 地面沉降(反)预报 6.1 预报时段的选取预报时段的选取 本次地面沉降预报其实是个反预报，主要是考虑到如下几点：(1)本次抽水试验结束后，鹏山水源地即停止了对地面沉降的监测，也就是说，如果仿照以往的地下水模拟研究那样，利用所建立的模型预测未来时刻的流场，则所进行的预报结果将不可能得到检验；(2)所选择的预报时间段内，积累了一些地面沉降监测资料，也有相应的地下水开采资料和水位观测资料供利用；(3)预报时段为本区的枯水期末，区内基本无降雨，河流无明显表流，可简化预测模型中包括降雨预测、边界条件预测、河流流量预测、开采方案预测等太多的条件预测，无形中减少了“根据预测值进行预测”的不可预测误差。总之，本文的地面沉降反预报亦可当成根据历史资料对模型的又一次检验。6.2 让历史评判让历史评判 与以往的预测模型最大的不同之处在于，表3所列出的预测值其实是对历史的一种回顾，计算是否准确，可以利用历史资料进行检验。表4列出的是2000年6月30日地面沉降预报结果与实际监测资料的对比，表中同时给出了用土工计算方法得到的相应时刻的计算结果参与对比。1146 表表3 2000年各监测点计算沉降量年各监测点计算沉降量(单位单位:mm)Table 3 Land subsidence volume of different monitoring points on 2000(unit:mm)预测时间/月-日 监测点号 06-05 06-10 06-1506-20 06-2506-30C10 0.068 0.163 0.2790.405 0.599 0.882C10 0.107 0.182 0.2870.473 0.727 1.395C11 0.085 0.190 0.3280.517 0.694 1.001C11 0.088 0.225 0.4290.677 0.796 0.978C29 0.070 0.120 0.2680.405 0.465 0.675C29 0.071 0.135 0.2560.377 0.524 0.713 表表 4 地面沉降量计算对比表地面沉降量计算对比表 Table 4 Comparison of land subsidence volume computation 地面沉降/mm 编号 水位降深/m 压缩层厚度/m 实测值 模型计算值 预测误差 土工计算值C10 0.18 9.80 0.853 0.882 0.029 1.059 C10 0.30 9.30 1.173 1.395 0.222 1.573 C11 0.25 8.00 0.958 1.001 0.043 1.038 C11 0.23 8.50 1.005 0.978-0.028 1.241 C29 0.18 7.50 0.648 0.675 0.027 0.811 C29 0.23 6.20 0.620 0.713 0.093 0.871 从计算结果看，由模型计算的地面沉降量普遍略大于实际观测值，其比值约介于0.971.19之间，笔者分析这主要是地层固结滞后于水位变化的结果，亦即水位降低可以在短时间内完成，但土层固结需要经历漫长的时间过程；另外，本模型未考虑给水度、渗透系数K等参数因压密固结引起的变化，也会产生一定误差。从整体来说，经模型计算的地面沉降量基本与实际观测数据相吻合。7 结 论(1)迄今为止，对地面沉降的监测都是针对粘性土地层进行的，尽管人们亦认同由砂性土组成的地层同样存在着地面沉降这一现象，但一直缺少直接的证据。本文以鹏山水源地地表变形研究为契机，首次成功地完成了由砂性土组成的潜水含水层在开采条件下的沉降观测，获得了较系统的监测资料，为研究砂性土含水层因过量开采地下水而造成的地面沉降积累了宝贵的资料；(2)含水层固结是由于含水层中水头下降引起的，本文分析的地面沉降模式完全依赖于准确的地下水流模型。为评价鹏山水源地地下水流模型的可信度，本文建立了计算值与实测值的相关方程，提出了偏离角和最大偏离值概念，并用于评价预测结果的可信性；(3)本文利用鹏山水源地开采条件下的含水层固结模型，再现了一段发生在研究区的地面沉降 历史，预测结果经置信检验，证明所建立的模型是可靠的。参参 考考 文文 献献 1 张世杰.山东省莱芜市自来水公司鹏山水源地供水水文地质评价报告R.济南：山东岩土工程勘察总公司,2000.2 唐哲.原平水源地潜水蒸发规律分析与研究J.工程勘察,2003,(2):2629.TANG Zhe.Study on the law of groundwater evaporation for the water resource in YuanpingJ.Geotechnical Investingation&Surveying,2003(2):2629.3 朱学愚,钱孝星,刘新仁.地下水资源评价M.南京:南京大学出版社,1987.4 莱芜市水利水产局.莱芜市水资源评价与开发利用研究R.莱芜：莱芜市水利水产局,1992.5 罗国煜,陈新民,李晓昭,等.城市环境岩土工程M.南京:南京大学出版社,2002.6 罗国煜,李生林.工程地质学基础M.南京:南京大学出版社,1990.7 钱家欢,殷宗泽.土工计算与原理（第 2 版）M.北京:中国水利水电出版社,1996.8 殷宗泽.土体的沉降与固结M.北京:中国电力出版社,1998.9 美国内政部水力资源局.地下水手册M.北京:地质出版社,1989.10 薛禹群,朱学愚.地下水动力学M.北京:地质出版社,1979.11 陈崇希.地下水不稳定井流计算方法M.北京:地质出版社,1983.12 陈崇希,唐仲华.地下水流动问题数值方法M.武汉:中国地质大学出版社,1994.13 孙讷正.地下水流的数值模型和数值方法M.北京:地质出版社,1981.14 孙峰根,王心田,王晓明.水文地质计算的数值方法M.徐州:中国矿业大学出版社,1995.15 张蔚榛.地下水非稳定流计算和地下水资源评价M.北京:科学出版社,1983.1147