地震波作用下承压含水层地下水渗流运动规律的研究.pdf

DOI：10.16030/ki.issn.1000-3665.202212020梁文宇，兰双双，谷洪彪，等.地震波作用下承压含水层地下水渗流运动规律的研究 J.水文地质工程地质，2023，50(4):73-83.LIANG Wenyu,LAN Shuangshuang,GU Hongbiao,et al.A study of the law of groundwater seepage movement in a confined aquiferunder seismic wavesJ.Hydrogeology&Engineering Geology,2023,50(4):73-83.地震波作用下承压含水层地下水渗流运动规律的研究梁文宇1，兰双双1，谷洪彪2，乔鹏1，毛政潭1（1.北京工业大学城市建设学部，北京100124；2.南京工业大学交通工程学院，江苏 南京211816）摘要：经典的地下水渗流理论是基于水均衡原理建立的，无法解释天然地震活动等外部荷载引起的井水位变化现象，也不利于深入认识地下水渗流运动在各种环境地质灾害中的作用。针对这一问题，基于流固耦合的动力学理论构建了地震波应力驱动承压含水层孔隙压力变化的数学模型，借助软件 Comsol 实现了对模型的有限元数值求解及影响因素分析，利用Cooper 理论反演出井水位的变化特征，并将结果与强震引起川滇地区井水位的变化特征进行对比。结果表明：当地震波荷载作用于承压含水层时，孔隙压力会表现出与地震波同周期的振荡变化特征，且地震波振幅、频率以及水力坡度对孔隙压力的变化特征影响显著。地震波荷载施加初期，孔隙压力迅速增大，随后含水层中的部分水被缓慢排出，荷载压力逐渐转移到颗粒骨架上，孔隙压力的变化速率随之变缓并趋于新的平衡。井水位变化特征与孔隙压力密切相关，其振荡周期与变化形态与孔隙压力一致，但幅度不同，总体上表现出振荡上升且趋于稳定的变化特征，与川滇地区观测到的井水位变化形态基本相同。该成果对于应力作用下地下水渗流运动理论的建立和完善具有宝贵的探索意义，可以丰富和拓展传统地下水动力学和经典流固耦合理论的研究思路和应用领域。关键词：渗流；孔隙压力；井水位；流固耦合；动力学理论；数值模拟中图分类号：P641.2 文献标志码：A 文章编号：1000-3665（2023）04-0073-11A study of the law of groundwater seepage movement in a confinedaquifer under seismic wavesLIANG Wenyu1，LAN Shuangshuang1，GU Hongbiao2，QIAO Peng1，MAO Zhengtan1（1.Faculty of Architecture,Beijing University of Technology,Beijing100124,China；2.College ofTransportation Engineering,Nanjing Techbology University,Nanjing,Jiangsu211816,China）Abstract：The classical groundwater seepage theory was established on the basis of the principle of waterequilibrium,which cannot explain the phenomenon of well water level change caused by natural seismicactivities and other external loads,and is not conducive to the in-depth understanding of the role of groundwaterseepage movement in various environmental geological disasters.To solve this problem,a mathematical model ofpore pressure change of a confined aquifer driven by seismic wave stress is constructed based on the fluid-structure coupling dynamic theory.The numerical verification of the model is realized by using the software 收稿日期：2022-12-11；修订日期：2023-01-15投稿网址：基金项目：北京市自然基金面上项目（8222003）；国家自然科学基金青年项目（41807180）第一作者：梁文宇（1998-），男，硕士研究生，主要从事水动力学数值模拟研究。E-mail：通讯作者：兰双双（1982-），女，博士，副教授，主要从事水动力学数值模拟研究。E-mail： 第 50 卷 第 4 期水文地质工程地质Vol.50 No.42023 年 7 月HYDROGEOLOGY&ENGINEERING GEOLOGYJul.，2023Comsol.The change characteristics of well water level are inversely performed by using the Cooper theory,andthe results are compared with the change characteristics of well water level caused by strong earthquakes in theSichian-Yunnan region.The influencing factors of seepage movement of a confined aquifer under earthquake arestudied by changing the simulation parameters.The results show that when seismic wave loads act on the confinedaquifers,pore pressure oscillates in the same period as the seismic waves,and the amplitude,frequency andhydraulic slope of the seismic waves have significant effects on the pore pressure,while the coefficient ofpermeability and porosity have little effect on the change characteristics.At the initial stage of seismic waveloading,pore pressure increases rapidly,and then part of the water in the aquifer is slowly discharged.The loadpressure gradually transfers to the granular framework,and the rate of change of pore pressure slows down andtends to reach a new equilibrium.The variation characteristics of well water level are closely related to porepressure,and the oscillation period and variation pattern are consistent with pore pressure,but the amplitude isdifferent.Generally,the oscillation rises and tends to be stable,which is basically the same with the variationpattern of well water level observed in the Sichuan-Yunnan region.The results are of valuable explorationsignificance for the establishment and improvement of groundwater seepage theory under stress,and can enrichand expand the research ideas and application fields of traditional groundwater dynamics and classical fluid-structure coupling theory.Keywords：seepage；pore pressure；well water level；fluid-structure coupling；theory of dynamics；numericalsimulation 经典的地下水渗流理论假定含水层所受总应力不发生改变，其是建立在水均衡原理基础上的单纯水动力学过程，对地下水资源开发利用1 2、水污染预测与治理3 4等方面都具有重要的意义。随着交叉学科的发展，该理论无法解释地震活动、地铁荷载、地表荷载（建筑物、地表水体等）等外部荷载引起的水位变化现象。在这些问题中，含水层与外界一般不存在水量交换，多孔介质需视为可变形体，其具体过程为：应力荷载导致多孔介质发生变形，造成含水层内部孔隙压力发生变化；孔隙压力变化将驱动地下水在含水层与井孔之间发生瞬时水量交换，导致井水位发生变化，同时孔压的调整过程反过来也会引起有效应力的改变，影响多孔介质变形。可见，天然地震活动作用下的地下水渗流运动是水动力学和岩土力学相互作用的复杂过程，是流固耦合作用的结果。1923 年太沙基有效应力原理的建立奠定了流固耦合的理论基础，随后学者们不断对其理论进行完善和深化，并建立了描述小变形的饱和多孔介质动力学理论5 8。在解析解的发展过程中 Zienkiewics 等9建立了结合中心差和单边差的显-隐式交替格式算法；Simon 等10针对饱和土的一维问题，提出了线性内插的方法处理 2 节点单元；Huang 等11基于孔压修正法提出了一种无条件稳定的隐-隐式算法，为解决动荷载作用下的流固耦合问题提供了坚实的理论参考。随着计算机仿真技术的提升，数值算法也得到不断发展，赵成刚等12利用解耦技术，建立了流体饱和两相多孔介质动力反应分析显式有限元法；刘宝等13借助Comsol 有限元软件，模拟分析了饱和多孔介质的不同动力耦合形式；杨贝贝等14借助 Comsol有限元软件建立了饱和动力问题的有限元分析模型。值得注意的是，前人研究主要的关注点是岩土体的变形，而细致地刻画含水层内部地下水渗流运动规律的研究较少。Darcy 于 1856 年通过长期试验得出了水在饱和多孔介质中的渗透定律，奠定了地下水运动理论的基础。随后 Dupuit 在 Darcy 定律的基础上研究了一维稳定流动和向水井的二维稳定运动，Theis 提出了地下水向承压水井的非稳定流公式15。同时，水文地质学者Cooper 等16、Roeloffs 等17、Doan 等18、Hsieh 等19和Wang20和 Wang 等21将应力或应变变量引入经典的地下水运动控制方程，推导出简谐波孔压、气压和引潮力等外荷载作用下井水位变化的解析解，并用于解释天然地震活动引起井水位的变化现象，但这些模型未考虑含水层水头变化对介质变形的反作用，没有达到细致地刻画地下水渗流场与应力场完全耦合过程 74 水文地质工程地质第 4 期的目的。鉴于上述，本文拟以川滇地区封闭性良好的承压含水层系统为概化的模拟对象，基于流固耦合动力学理论和数值模拟的方法，重点研究地震波荷载作用下含水层内部孔隙压力变化控制方程的建立及求解，探讨地震波引起地下水在含水层内部及其在观测井与含水层之间的渗流运动规律及影响因素，揭示地震波对承压层孔隙压力及井水位变化的影响，旨在解释地震动态荷载引起井水位变化的多样性和复杂性，为进一步深化井孔-含水层系统对地震活动响应机制提供参考，同时拓展传统地下水动力学和经典流固耦合理论的研究思路和应用领域，深化认识渗流在流、固、热、化等多场耦合作用中的地位和意义，对预防和减少地面沉降、水电工程大坝失稳、矿井事故及边坡失稳等灾害的发生提供参考依据。1 流固耦合动力学方程的建立与求解 1.1 理论方程的建立地震波荷载驱动承压含水层渗流运动是动荷载作用下流固耦合共同作用的结果，其理论基础是多孔介质动力学理论。假设岩土颗粒是不可压缩的，孔隙水是可压缩的、有黏性的；岩土骨架是均质各向同性的弹性多孔介质，变形符合广义胡克定律；渗流符合Darcy 定律；应力连续条件符合太沙基有效应力原理；温度的影响忽略不计。规定应力方向以拉为正、压为负，含水层体积应变以膨胀为正、压缩为负，孔隙压力增加为正、减少为负。当不考虑重力场的影响时，基于质量守恒定律、达西定律、有效应力原理、胡克定律、动量守恒定律，推导出二维条件下以介质骨架位移 u 和孔隙压力 p 为变量的控制方程22：(2uxx2+2uxz2)+(+)(2uxx2+2uzxz)px=m2uxt2(2uzx2+2uzz2)+(+)(2uxxz+2uzz2)pz=m2uzt2t(uxx)+t(uzz)kff(2px2+2pz2)fkffx(2uxt2)+z(2uzt2)=nEfpt（1）=E(1+)(12)=E2(1+)m=(1n)s+nf（2）式中：u介质骨架位移/m；p孔隙压力/Pa；m多孔介质密度/（kgm3）；s固相密度/（kgm3）；f液相密度/（kgm3）；固相骨架材料的剪切弹性模量/MPa；Lame 常数；Ef流体压缩模量/MPa；f黏滞系数/（Pas）；kf 渗透率/m2；E杨氏模量/MPa；泊松比。1.2 数值求解方法由于饱和多孔介质动力响应问题的复杂性，很难用解析方法对描述该问题的二阶偏微分方程组进行求解，故本次利用模拟软件 Comsol Multiphysics 进行有限元数值求解，并通过经典的二维饱和多孔介质动力响应的解析解验证数值求解方法的合理性和可行性。但由于该软件自带物理模块无法解决流固耦合动力学方程的求解，需利用软件提供的偏微分方程（PDE）模块进行二次开发实现数值求解。Comsol 软件的 PDE 模块提供了 3 种形式定义偏微分方程：系数形式、一般形式以及弱形式。由于弱形式 PDE 应用范围最广，解决能力最强，尤其在求解空间和时间混合偏导时，结果最为精确，故本次采用 PDE 弱形式进行求解。软件中弱形式采用虚位移原理对偏微分方程的表达式进行描述，其内置的通用求解形式如下：wweakds=0（3）式中：weak求解区域中的弱表达式；求解区域；s积分变量。针对动荷载作用下流固耦合方程 u-p 形式的特征，首先推导出响应动力耦合控制方程的有限元弱形式，从而获取 u-p 方程组的数值解。对式（1）（3）乘以虚位移 u和虚孔隙压力 p，再利用 Green-Gauss 公式对方程进行分部积分，得到 u-p 方程组的弱形式如下：2023 年梁文宇，等：地震波作用下承压含水层地下水渗流运动规律的研究 75 w(uxx+uxz)+(+)(uxx+uzz)(uxx+uzz)ds+w(uxx+uxz)+(+)(uxx+uzz)uxdl+wp(uxx+uzz)dswpuxdlw2uxt2ds=0w(uzx+uzz)+(+)(uxz+uzz)(uxx+uzz)ds+w(uzx+uzz)+(+)(uxz+uzz)uzdl+wp(uxx+uzz)dswpuzdlw2uzt2ds=0wt(uxz)+t(uzz)pds+wkff(px+pz)(px+pz)dswkff(px+pz)pdlwfkffx(2uxt2)+z(2uzt2)pds+wnEfptpds=0（4）1.3 数值解的验证Xu 等23在 Huang 等24和 Zienkiewics 等25的基础上提出了双显式计算法，并给出了简单二维不可压缩饱和多孔介质动力响应的解析解。本文利用此解析解与数值解进行对比，从而达到验证数值方法的目的。算例模型如图1 所示，具体描述如下：建立一个30 mux=uz=0ux=0ux=uz=0p=030 m 的二维含水层模型如图 1（a）所示，含水层均质且各向同性，上部为自由透水边界且施加荷载 q(t)=1.0+0.5sin(2t)kPa 的循环荷载，底部为固定不透水边界（），两侧只能发生纵向位移（）且为不透水边界。孔隙位移和压力初始值均为 0（，），模型参数见表 1。015300.51.01.5F(t)/kPa 时间/s（b）循环荷载（a）几何模型及边界条件30 m30 m观测点C（0,27）B（0,15）zxA（0,3）图 1 循环荷载作用下二维算例的物理模型Fig.1 Physical model of a 2D example under cyclic load 表 1 二维算例的含水层参数Table 1 Aquifer model parameters for a 2D example 参数数值剪切模量/MPa65渗透系数/（ms1）0.001液体密度/（kgm3）1 000泊松比0.3孔隙度0.3Biot系数1固体密度/（kgm3）2 650重力加速度/（ms2）9.8 本次模拟采用自由三角网格进行剖分，计算时间为 30 s，计算步长为 0.1 s。分别取 A（0，3）、B（0，15）、C（0，27）3 点为计算点，将其孔隙压力及垂向位移的模拟值与文献 23 的解析解进行对比（图 2）。结果表明：A、B、C3 点的孔隙压力变化与颗粒的垂向位移皆呈现出初期迅速增大，随后缓慢下降且趋于新平衡的变化特征；且 u 和 p 的数值解与文献 23 的解析解在数值和趋势上基本一致，表明利用数值方法求解方程u-p 方程组具有可行性。2 地震波引起井水位变化的实例分析中国的川滇地区地处地中海喜马拉雅地震带，地震活动频率高，且地震地下流体观测网密度大，对地震响应敏感的观测井数量较多，故选取该区华蓥山 76 水文地质工程地质第 4 期断裂带附近北碚井（BB）、荣昌井（RC）、南溪井（NX）和巴南井（BN）为研究对象（图 3），以 2010 年 Ms8.0 级智利地震、2011 年 Ms9.0 级日本地震、2012 年 Ms7.4级日本地震和 2012 年 Ms8.6 级苏门答腊地震为震例，分析井水位对地震活动的响应特征。根据重庆市地震监测志26可知，观测井的主要观测层位岩性为侏罗系中统沙溪庙组长石砂岩，厚度几十米不等，为承压井。井水位观测数据来自于中国地震前兆监测台网，由 LN-3A 型水位仪观测，采样频率为 1 次/min。强震引起井水位的变化特征见图 4。2010 年智利地震导致井水位变化以振荡为主，幅度较小，见图 4（a）。2011 年日本地震发生时井水位变化以振荡为主，北碚、南溪、巴南井水位同时伴随有明显的阶升，阶升幅度较大，为 512 cm，见图 4（b）。2012 年日本地震导致北碚、荣昌和南溪井水位出现振荡现象，变化幅度较小，最大也不超过 2.7 cm，且几分钟内水位皆恢复到震前正常状态，见图 4（c）。2012 年苏门答腊地震引起的井水位变化以振荡为主，北碚、南溪井伴随有3.2，4.2 cm 的阶升，巴南井伴随有 6.2 cm 的阶降，阶变幅度较小，见图 4（d）。可见 4 次远场强震引起井水位变化均以振荡为主，2011 年日本地震和 2012 年苏门答腊地震发生时北碚井、南溪井和巴南井水位发生振荡的同时还伴随有阶变响应，其变化幅度为 3.212.5 cm。3 地震波驱动承压水渗流的数值模拟 3.1 计算模型的建立(pn=0)为了减少野外场地条件的复杂性及影响因素的多样性给计算结果带来的影响，本次假设含水层水平分布、均质、各向同性，固体颗粒变形为弹性小变形。研究区含水层主要岩性为长石砂岩，厚度几十米不等，且考虑实际地震流体观测井大部分处于承压含水层，故本次模拟对象为渗流条件下的承压含水层，含水层剖面可概化为 50 m100 m 的矩形区域（图 5），地下水流为二维稳定流运动。模型顶部和底部均为不透水边界，左右两侧为定水头边界（h=0.3 m），根据重庆市地震监测志26可知该区水力梯度范围为 1/1 0005/1 000，故将水力坡度概化为 3/1 000。为了更直观地反映孔隙压力的变化情况，本次模拟不考虑重力场的影响。在地震发生的过程中会产生横波、纵波和面波（主要为 Rayleigh 波）。Cooper 等16 0510152025300.500.51.01.52.02.5孔隙压力/kPa时间/s解析解（A）解析解（B）解析解（C）数值解（A）数值解（B）数值解（C）解析解（A）解析解（B）解析解（C）数值解（A）数值解（B）数值解（C）（a）孔隙压力051015202530时间/s（b）纵向位移0.140.120.100.080.060.040.020纵向位移/mm图 2 A、B、C 孔隙压力和纵向位移随时间变化模拟结果与解析解的对比Fig.2 Comparison between simulation results and analyticalexpression of pressure and longitudinal displacement changeswith time of point A,B and C 090180 km观测井N四川省重庆市贵州省南溪井荣昌井北碚井巴南井图 3 观测井的地理位置Fig.3 Geographic map of the water level observation well2023 年梁文宇，等：地震波作用下承压含水层地下水渗流运动规律的研究 77 f(t)=1 000+500sin(20.1t)Paux=0ux=uz=0认为远场地震作用下瑞利波是造成含水层发生体积变形的主要影响因素，Rayleigh 波是在地表传播的面波，故在顶部施加量化的动态荷载26，左右两侧边界受限只发生纵向位移（），初始位移为 0（），初始孔隙压力为 0 Pa，分别选取点 A（0，10），B（0，25），C（0，40）作为计算点。含水层水文地质参数及固体介质力学参数见表 1。f(t)=1 000+500sin(20.1t)Pa只可纵向位移只可纵向位移不透水边界不透水边界xzC（0,40）B（0,25）A（0,10）图 5 含水层物理模型的概化及空间网格剖分Fig.5 Generalization and spatial grid subdivision of the physicalmodel of an aquifer 3.2 计算结果与讨论采用 u-p 方程组描述地震波荷载驱动承压含水层渗流问题，基于方程的弱形式式（4），借助软件 Comsol的 PDE 模块，计算时间为 200 s，步长设置为 0.1 s，采用自由三角网络剖分方式进行有限元数值求解。由于本次模拟不考虑重力场的影响，故计算求得的孔隙压力变化值均为地震波荷载引起的。图 6 计算结果表明：A、B、C3 点处孔隙压力均表现为振荡上升，振动周期为 10 s，与地震波周期一致，孔隙压力最大值为 1.85 kPa，随后缓慢下降并最终趋于稳定；3 点中C 点孔隙压力振幅最大，为 0.76 kPa，其次是 B 点，A点值最小，为 0.18 kPa。A、B、C3 点所处含水层中的位置不同，3 点孔隙压力随时间的变化形态也有所不同，C 点距离荷载边界最近，受到地震波荷载影响最明显，在模拟初期其孔隙压力增量最大，同理，A 点处孔隙压力变幅最小。0501001502000.20.40.60.81.01.21.41.61.82.0孔隙压力/kPa时间/s A点 B点 C点1.85 kPa图 6 地震波作用下承压含水层中孔隙压力的变化Fig.6 Changes in pore pressure in pressurized aquifers underseismic wave 地震波荷载施加初期，由于孔隙水不能及时排出，荷载压力造成孔隙压力迅速增大，产生超静水孔隙压力。随后在水头差作用下地下水由孔隙中缓慢排出，荷载压力逐渐转移到颗粒骨架上，颗粒骨架受力压缩，孔隙体积随之变小，含水层仅发生体积变形，1.801.831.861.281.321.364.004.024.044.06北碚井水位埋深/mMs8.8（a）2010年智利荣昌井南溪井2月27日14:342月28日2:342月27日2:341.551.601.651.700.901.051.203.753.853.950.700.60北碚井Ms9.0荣昌井水位埋深/m水位埋深/m水位埋深/m水位/m水位/m南溪井（b）2011年日本（c）2012年日本（d）2012年苏门答腊巴南井3月11日13:463月12日1:463月11日1:460.920.940.961.041.061.083.913.933.95北碚井Ms7.4荣昌井南溪井12月7日4:1812月7日16:1812月8日4:181.141.171.201.231.081.201.324.024.054.083.403.363.323.28Ms8.2北碚井Ms8.6荣昌井南溪井巴南井4月11日4:384月11日16:384月11日4:38图 4 地震活动引起井水位的变化曲线Fig.4 Changes in well water level caused by seismic activities 78 水文地质工程地质第 4 期其结构没有或只有轻微的破坏，此时孔隙压力变化速率随之变缓并趋于稳定。由孔隙压力的变化可知，整个模拟期外加荷载的强度并未超过含水层临界动力强度，含水层的变形一直以弹性小变形为主。Ducan-Chang 模型27是基于 Kondner 三轴应力-应变关系曲线建立的非线性本构模型，可以用来解释单调加载和循环加载三轴试验所得到的数据(13)-的变化规律。如图 7 所示，当循环加载时，加载后卸载到原应力状态，介质变形一般不会恢复到原来的应变状态。可见每次应力循环都有可恢复的弹性应变和不可恢复的塑性应变。故当模拟含水层受到周期性地震波作用时，A、B、C3 点孔隙压力呈现规律性的周期波动变化特征，且含水层在发生弹性体积应变的同时也产生了不可恢复的变形，使含水层在地震波荷载作用下达到新的平衡。400（13）/kPav/%l/%20001246810循环加载单调加载pe图 7 单调加载与循环加载三轴试验应变-应力变化对比曲线27Fig.7 Comparison of strain-stress changes in triaxial tests undermonotonic and cyclic loading27注：e为弹性应变，p为塑性应变，(13)为主应力增量，1为轴应变，v为体应变。4 孔隙压力的影响因素分析地震波引起井水位变化形态复杂多样，会受到震源特征、含水层水文地质条件、井孔所处构造条件及井自身结构等多种因素的影响。由于影响因素多，地下水对地震波作用下渗流运动的响应机制十分复杂。故本次通过改变含水层特征参数和水力条件（如渗透系数、孔隙度、水力坡度）、地震波特征（如振幅与频率）讨论不同因素对 u-p 耦合结果的影响，分析不同因素对地震波引起承压含水层内部孔隙压力变化的影响。4.1 地震波振幅对孔隙压力的影响地震波振幅是影响井水位变化的重要因素，地震波振幅通常为 102106 Pa26，故将地震波振幅的工况分别设计为 P=250，500，1 000 Pa，地震波的最大峰值分别为 750，1 500，3 000 Pa。以点 B（0，25）为例，计算并对比 3 种工况下孔隙压力的变化特征。图 8 显示，孔隙压力皆呈现振荡增大并趋于稳定的变化趋势；当地震波振幅增大时，孔隙压力的振荡频率及最终稳定值几乎相同，但振荡振幅发生了明显的变化。当P=250 Pa 时，孔隙压力的振幅为 0.24 kPa，最大孔隙压力值为 1.60 kPa；当地震波振幅增大到 P=1 000 Pa 时，孔隙压力振幅相应地增加到 0.92 kPa，最大孔隙压力值为 1.96 kPa，可见孔隙压力振幅和孔隙压力的变幅都会随着地震波荷载振幅的增大而增大。地震波振幅增大含水层受力变形程度也相应增大，土骨架被挤压发生体积变形，导致孔隙压力振幅和最大变幅随之增大。01002000.51.01.52.0孔隙压力/kPa时间/s250 Pa500 Pa1 000 Pa1.96 kPa1.60 kPa图 8 不同地震波振幅作用下孔隙压力的变化曲线Fig.8 Changes in pore pressure under different seismicwave amplitudes 4.2 地震波频率对孔隙压力的影响Cooper 等16认为远场地震波的频率一般介于1/1001/10 Hz 之间，故本次地震波荷载频率的设计工况分别为 f=1/10，1/30，1/50，1/60 Hz，计算并对比4 种工况下 A（0，10）、B（0，25）和 C（0，40）3 点处孔隙压力变化特征。图 9 显示，孔隙压力皆呈现振荡增大并趋于稳定的变化趋势；当地震波频率从 1/60 Hz 增 0100200300400500600孔隙压力/kPa1.02.01.85 kPaf=1/10 Hz A点 B点 C点1.02.0f=1/30 Hz1.02.0f=1/50 Hz时间/s1.02.01.81 kPaf=1/60 Hz图 9 不同地震波频率作用下孔隙压力的变化曲线Fig.9 Changes in pore pressure under different seismic wavefrequencies2023 年梁文宇，等：地震波作用下承压含水层地下水渗流运动规律的研究 79 大到 1/10 Hz 时，孔隙压力的振动频率发生了明显的变化，始终与模拟地震波荷载一致；随着地震波频率的增大观测点 C 点的孔隙压力的最大值由 1.81 kPa 增加到 1.85 kPa，在孔隙压力值变化不明显，且孔隙压力达到平衡状态所需的时间也基本一致。随着振动频率的增大，含水层受力及岩土颗粒发生垂直位移的频率和幅度增大，导致岩土颗粒被压密的程度增大，孔隙压力不断累积无法消散，导致孔隙压力发生一定程度的增大。4.3 渗透系数对孔隙压力的影响u-p 形式偏微分方程适用于渗透系数较小的情况8，10，孔隙水和固相骨架之间的相对位移比较小，产生的惯性耦合效应比较弱，可以认为渗透系数（K）小于 103 m/s 时，惯性耦合效应可以忽略。遵循 u-p 方程的适用条件，并考虑研究区含水层岩性以砂岩为主，分别研究 K=1.0103，5.0104，1.0104，5.0105 m/s4 种工况28下 A（0，10）、B（0，25）、C（0，40）3 点处孔隙压力随时间的变化情况，如图 10 所示。结果显示，孔隙压力皆呈现振荡增大并趋于稳定的变化趋势；距离顶部荷载的 C 点孔隙压力最大，A 点的孔隙压力最小。在渗透系数由 5105 m/s 增加到 1.0103 m/s 时，当孔隙压力振荡达到稳定状态时，2 种工况下孔隙压力的最大值改变不大，分别为 1.88，1.86 kPa，所需时间分别约为 1 200，132 s（由于 K=1.0104，5.0105 m/s 2个模拟结果到达平衡所需的时间过长，考虑到图片的清晰程度，这 2 个工况只展示到 600 s）。可见渗透系数的变化对孔隙压力的平衡变幅影响不大。但渗透系数大时，孔隙压力在较短时间内骤增，且达到平衡所需时间短。根据达西定律可知，水头差为定值时，渗流速度会随着渗透系数的增大而变大，在相同条件下地下水在含水层中运动所受阻力减小，故孔隙压力骤升且达到平衡状态所需要的时间变短。4.4 孔隙度对孔隙压力的影响鉴于研究区含水层主要岩性为砂岩，本次分别探究孔隙度 n=0.1，0.2，0.3，0.4 工况28下 B 点孔隙压力随时间的变化情况。图 11 计算结果显示，孔隙压力呈现振荡增大并趋于稳定的变化趋势，当孔隙度由 0.1增大到 0.4 时，观测点 B 点的最大孔隙压力值由 1.72 kPa下降为 1.70 kPa；振幅由 0.464 kPa 下降为 0.446 kPa。表明孔隙度变化对孔隙压力的变化幅度存在一定影响，但是影响很小。不同孔隙度下孔隙压力振荡的振幅以及频率基本一致，孔隙压力的响应时间也基本一致，当 n 从 0.1 增大到 0.4 时，孔隙压力到达平衡所需的时间均为 142 s。可见，孔隙度对于孔隙压力的变幅、振荡频率以及到达平衡所需的时间影响很小。1.02.01.02.01.02.00501001502001.02.01.72 kPan=0.1142 sn=0.2n=0.3时间/s孔隙压力/kPa1.70 kPan=0.4142 s图 11 不同孔隙度条件下孔隙压力的变化曲线Fig.11 Changes in pore pressure under differentporosity conditions 4.5 水力坡度对孔隙压力的影响渗流条件中需要考虑一个重要的影响因素就是水力坡度（J）。根据重庆市地震监测志26及文献 29可知该区水力梯度范围为 1/1 0005/1 000，故本次水力坡度工况设计为 J=1/1 000，2/1 000，3/1 000，分别计算并分析 3 种工况下 A、B、C3 点的孔隙压力的变化特征。图 12 显示，当水力坡度为 3/1 000 时，3 点孔隙压力皆呈现振荡增大并趋于稳定的变化特征，最终稳定在 1.85 kPa 附近。随着水力坡度的减小，孔隙压力变化幅度也随之变小，当水力坡度为 1/1 000 时，孔隙压力最终稳定在 0.88 kPa 附近。这是由于当水力坡度增大，承压含水层中地下水的渗流速度和流量也会随之增大，最后导致孔隙压力的增大。值得注意的是，水力坡度的不同，A、B、C3 点孔隙压力随时间的变化形态也有所不同，尤其当水力坡度为 2/1 000，1/1 000 1.02.01.02.01.02.003006001.02.0孔隙压力/kPaK=1103 m/sK=5104 m/sK=5105 m/sK=1104 m/s1.86 kPa132 s时间/s A点 B点 C点图 10 不同渗透系数条件下孔隙压力的变化曲线Fig.10 Changes in pore pressure under different hydraulicconductivities 80 水文地质工程地质第 4 期时尤为明显。随着水力坡度减小为 J=1/1 000 时，离荷载边界最近的 C 点处孔隙压力的初始变幅最大，随后孔隙水压力整体表现为下降的变化趋势，直至达到新的平衡。该结果也表明，当含水层水力梯度很小时，孔隙压力的变化形态可能与观测点的埋深有关，不同深度的观测点可能表现为振荡上升或振荡下降或仅为振荡，这或许是井水位变化形态多样性的原因。0.51.01.50.51.01.50501001502000.51.01.5J=1/1 000J=2/1 0000.88 kPa孔隙压力/kPaJ=3/1 00032 s时间/s1.85 kPa162 s A点 B点 C点图 12 不同水力坡度条件下孔隙压力的变化曲线Fig.12 Changes in pore pressure under different