基于井水位潮汐响应探讨含水层的水力特征及其变化

赵可新，兰双双，谷洪彪，梁文宇，乔鹏

（1.北京工业大学，北京 100124；2.防灾科技学院，河北 三河 065201）

含水层的水力特征参数是表征地下水运动规律的重要指标，包括贮水系数、导水系数、渗透系数、越流系数等。当含水层所受应力状态发生变化时，可能会引起其水力特征发生变化，进一步导致地下水流发生相应的变化［1-2］，从而威胁水库坝基的稳定性以及地下水的安全性，甚至可能对深部地热、石油、天然气等重要地下资源的运移和储存［3-4］产生一定的影响，故研究应力荷载作用下含水层水力特征的变化具有非常重要的实际工程意义。

含水层水力特征参数的传统获取方法主要有抽水试验和注水试验，但这些方法无法推求含水层渗透系数的实时连续变化，因此无法用来对外加荷载作用下含水层渗透性的变化进行分析，而含水层水头对周期高频性自然力（气压、固体潮等）的响应模型为该方面的研究提供了新的思路［5-6］，井水位的气压效应可以用来推求含水层的骨架材料特性［7-8］，井水位的潮汐响应可以用来计算封闭性承压含水层的水平渗透性［9-11］。最新研究表明地震应力不仅可以改变地下水在井孔-含水层之间的水平流动状态［12-13］，也可能通过裂隙的张开或闭合引起相邻含水层垂向上的水力联系［14］，Wang等［14］引入越流含水层模型，并通过越流系数的解析解分析了美国Oklahoma深井的渗漏现象，由此地震引起含水层垂向水力特征参数的改变引起学者们的高度重视。人们［15］通过川滇陕甘渝地区井水位的潮汐响应发现，大部分含水层并不是完全承压的，在局部地区地下水存在着垂向运动，而以往的成果中对该方面的定量计算及讨论涉及较少，地震引起含水层垂向水力特征参数改变的可参考研究成果也相对较少。

综上，为了更加客观地、定量地反映地震引起含水层水平及垂向水力特征的改变，本文在系统分析总结井水位潮汐响应与含水层水力特征理论关系的基础上，以华蓥山断裂带附近大足井和北碚井为例，判断井孔附近含水层系统中地下水流运动状态，计算含水层水力特征参数，并以汶川地震为震例，讨论地震应力对含水层水平及垂直水力特征及地下水流运动状态的改变，探究地震引起两口井水位的变化原因。该研究成果可拓展传统地下水动力学的研究领域，进一步提高对外部荷载应力驱动含水层中地下水渗流运动的认识，为井孔-含水层系统对地震活动响应机制的研究提供参考依据。

1 理论基础

在日、月引潮力作用下地球会发生周期性的压缩与膨胀，造成含水层内部孔压随之增大与减小，引起观测井与含水层之间产生压力差并发生水流交换，最终导致井水位发生周期性波动，这种现象称为井水位的固体潮效应。具体可分解为两个过程：一是含水层内部水头对潮汐体应变和潮汐垂向流的响应，前者是不排水孔弹性响应过程，后者是含水层垂向上的不均匀导致含水层水头重新分布的过程，该过程只在垂向流存在的含水层中发生；二是井水位对含水层压力水头的响应，它主要受潮汐径向流的影响，是一个典型的水动力学过程，在垂向流和径向流中均存在此过程，井水位对压力水头的响应主要取决于含水层水动力学参数（如导水系数、贮水系数等）［16-17］。可见，井水位对固体潮的响应实质上反映了含水层的水力特征及地下水流运动状态，故通过分析井水位的潮汐特征可以达到推求含水层水力参数、判定地下水流运动特征的目的，从而为解释地震应力对井孔-含水层系统的影响提供依据。

1.1 计算贮水系数

贮水系数S和贮水率Ss是描述含水层释水能力的参数，贮水率乘以含水层厚度等于贮水系数。

1.1.1 地球对潮汐的响应

在均质、弹性条件下，地球的潮汐水平应变［18］可以表示为

式中：εhh为潮汐水平应变；hn和l n为勒夫数；Wn为引潮位；a为地球半径；n为含水层孔隙度。

垂向应变一般难以计算，利用泊松比v建立水平与垂向应变之间的关系，则体应变［19］表达式为

式中：εkk为潮汐体应变；εrr为潮汐垂向应变。

对于二阶引潮位的变化（ΔW2），利用式（2）得到潮汐体应变［20］为

1.1.2 地下水对固体潮的响应

承压含水层水位升降会引起含水层两部分体积变化［21］：一部分是水位上升后井中水体的体积变化量dV1＝Ah；另一部分是水位上升导致的静水压力增加使得含水层内水体体积变化，即dV2＝－gρwh·（Vw／Kw）＝－gnρwh（V／Kw），则固体潮体应变引起含水层体积的变化可表示为

式中：ρw为水的密度，kg／m3；Kw为水的体积模量，Pa；V为含水层体积，m3。因为KwA＋gρwV，则承压含水层水位变化为

假设含水层体积变化等于含水层孔隙度变化，由贮水率的定义可知含水层的贮水率［22］为

将式（6）代入式（5）可得井水位和贮水率的关系为

将式（3）代入上式可得到引潮位表示的贮水率［20］

式中：泊松比v＝0.25，h2＝0.606，l2＝0.084，地球半径a＝6 371.393 km；水头变化Δh是引潮位ΔW2引起的，W2＝g Kmbf（θ），Km是与地球质量、月球质量、地球与月球之间的距离、地球半径有关的常数，约等于0.537 m［23］；b是与潮汐波分量周期相关的常数，对于M2波，b＝0.908；f（θ）是纬度θ的函数，f（θ）＝0.5cos2θ。

1.2 计算导水系数

在引潮力的作用下，含水层岩石产生形变，压力水头改变，驱动井与井周含水层之间水流交换，地下水运动以水平流为主，简称潮汐径向流。值得注意的是，引潮力作用于含水层到井-含水层之间发生水流交换，存在着时间滞后，即相位差。Paul等［24］提出潮汐径向流条件下压力水头振幅比A和相位差ηr与含水层参数的解析解为

其中

式中：ω为潮汐分波频率，d－1；τ＝2π／ω是潮汐分波周期；Kei和Ker为零阶Kelvin的函数；T为径向导水系数，m2／d；rω为井滤水管半径，m；rc为井筒（水位变动段）半径，m；H为水井水位的变化量，m。

由式（9）、（10）可知，对于潮汐径向流而言，振幅比A和相位差ηr是关于导水系数T和贮水系数S的函数，图1显示了相位差ηr和无量纲参数之间的关系，对于径向流而言，相位差ηr＜0；当在一个可能的值域范围内，相位差随T或的增加而增大，故根据ηr的取值可估算参数，进而确定含水层的T值。

图1 不同条件下，径向流相位差（ηr）随的变化［9］Fig.1 Under different conditions，variation of phase shift in radial flow with

1.3 计算越流系数

自然界中存在的井-含水层系统上、下岩层往往并不是绝对隔水的，其中一个或两个可能是弱透水层：当这个含水层和相邻含水层之间存在水头差时，地下水便会发生越流运动；或当观测井本身为非完整井时，观测含水层与下伏含水层之间或许存在着水量交换。由此，Wang等［14］推导出越流含水层系统中潮汐力驱动地下水在含水层和井孔间运动的控制方程，并得出了振幅比和相位差的解析解

其中

在式（14）、（15）中，当K'＝0即没有越流存在时，该方程与Hsieh的模型解析解一致，如图2所示相位差是关于越流系数σ′＝K′／b′、贮水系数S以及导水系数T的函数：如果贮水系数S和导水系数T已知，当相位差为负值时，此时垂向上的水力联系较少，含水层以水平径向流为主；随着相位差增大为正值，此时σ′的值变大，含水层垂向渗漏增强。因此，在S和T已知的含水层系统中，由相位差可以得到σ′的值。

图2 越流含水层系统中相位差与越流系数的关系Fig.2 Relationship between the coefficient of leakage and phase shift in leaky aquifer system

2 观测背景与数据

华蓥山断裂带位于中国西南地区，是一条右旋走滑型逆断裂带。该断裂带北起万源，向南西经达川、荣昌至宜宾西南，全长约460 km［25］，是川东隔挡式褶皱带和川中平缓构造的分界断层。地震流体观测井大足和北碚皆位于华蓥山断裂带附近10 km范围内［图3（a）］，隶属四川盆地的重庆市，钻孔揭露的地层主要为侏罗系中统沙溪庙组紫红色泥岩和浅灰白色长石砂岩互层，即“红层”。其中，长石砂岩为该区地下水赋存的主要含水层，单层厚度约10～50 m。含水砂岩上下均被相对隔水的泥岩所夹持，具有层间承压性，且该类型含水层广泛分布，在较大范围内存在着水力联系。多孔抽水试验观测数据［26］显示：位于同一含水砂岩层中相距数百米或千米之外的钻孔，一井抽水会引起另一井水位发生显著下降。该区域红层承压水［27］主要在露头区接受大气降水的补给，其次是地表水的垂直入渗和部分越流补给，降水和地表水通过含水层暴露于地表发育的裂隙系统下渗。在重力作用下顺层或沿砂岩裂隙发生径流或垂向越流补给，当含水层被切割时，地下水以泉的形式排泄于地表或地表水体，也可以通过上下相邻含水层进行越流排泄。

图3 观测井位置与柱状Fig.3 Location of the observation wells and stratigraphy of the borehole

大足井位于华蓥山断裂带中部，西山背斜北西翼［16］。如图3（b）所示，该井深108.7 m，46 m下设套管，直径127 mm，46 m以下为裸孔，井径150～120 mm。其中：地表出露地层主要是中侏罗统沙溪庙组泥岩夹砂岩；5.50～47.00 m为青灰色长石砂岩；层中32～35 m、38 m深度段集中发育倾角为30°～75°的裂隙；46.0～108.7 m岩性以泥岩为主。地下水类型属孔隙裂隙承压水，为非自流井。

北碚井位于华蓥山断裂带东侧，观音峡背斜东翼。如图3（c）所示，井深105.36 m，Φ127的套管下至42.1 m处，套管外用水泥固井止水。井孔岩性自地表往下：0～3.54 m为覆土；3.54～29.90 m为紫色泥岩，岩层倾角35°，层面风化裂隙发育；29.90～70.24 m为主要观测含水层，井孔岩性为中侏罗统沙溪庙组浅灰白色中-粗粒厚层状长石砂岩，该段裂隙较发育，裂隙倾角55°～80°，且30.40、43.30 m处的裂隙被方解石半充填，38.40 m处的裂隙被黄铁矿晶粒充填；70.24～105.36 m为紫红色砂质泥岩，夹紫灰色粉砂岩条带［16］。地下水类型属裂隙承压水，为非自流井。

北碚井和大足井水位观测方式均为静水位观测，观测仪器皆为LN-3A型水位仪，采样频率为每分钟1次，观测数据来自于中国地震前兆监测台网。如图4所示，北碚井和大足井年变化趋势基本一致，受水文因素影响明显，丰水期水位较高枯水期水位较低；井水位日动态变化稳定，正常变化幅度约6 cm，具有清晰的两峰两谷的潮汐变化形态。且北碚井和大足井水位的频谱分析结果（图5）表明：对两口井水位影响较大的潮汐成分是半日波M2和S2，周期分别为0.517 5 d和0.500 0 d，故可以利用井水位的潮汐效应探讨含水层的水力特征变化。

图4 大足井和北碚井水位的年变化曲线Fig.4 Annual variation of the water level in wells Dazu and Beibei

图5 大足井和北碚井水位频谱分析结果曲线Fig.5 Results of spectrum analysis of the water level in wells Dazu and Beibei

3 含水层水力特征参数的计算

根据中国地震台网可知2012年1月1日至5月1日川滇地区无Ms 5.0以上地震发生，且该时段属于枯水期，井水位受降雨因素干扰少，故本次选取该时间段内大足井和北碚井水位小时值观测数据进行潮汐分析，进而推求含水层的水力特征参数。潮汐波是由一系列不同周期的潮汐分量组成，其中M2波（半日波）水位数据振幅较大，受气压等因素干扰小，误差相对较小，因此本次借助国际通用潮汐分析软件Baytap-G程序进行井水位对M2波的潮汐分析。数据间隔设置为30 d，滑动窗长设置为15 d，分别提取两口井水位对M2波的响应特征参数（A和η），计算结果见表1。

表1 大足井和北碚井对M2波的潮汐分析结果Tab.1 Results of tidal response of M2 wave in wells Dazu and Beibei

据表1可知大足井水位对M2波响应的相位差η＜0，表明该井附近含水层地下水流类型以水平径向流为主。大足井管揭穿了整个含水层厚度，含水层分布稳定，地层倾角约为20°，地下水流运动方向受岩层倾角影响，接近水平流。利用大足井M2波的理论引潮高振幅和该井观测段水头响应振幅数据，据公式（8）推求该井的观测含水段的贮水能力，结果见表2，大足井观测含水层的贮水率为1.55×10－6m－1；由该井钻孔资料将观测含水层厚度取值为41.5 m，则含水层的贮水系数为6.45×10－5。由此利用公式（9）和（10）及相位差可计算观测井附近含水层的导水系数见表3，导水系数平均值为1.74 m2／d，渗透系数（表3）为0.041 9 m／d，与前人［28］抽水试验的结果0.053 4 m／d一致，表明利用井水位的潮汐响应计算含水层水力参数具有可行性。

表2 大足井和北碚井所处含水层贮水能力计算Tab.2 Calculation results of the storage capacity in the aquifer near wells Dazu and Beibei

表3 大足井所处含水层渗透系数Tab.3 The coefficient of permeability in the aquifer near well Dazu

北碚井的潮汐相位差η＞0（表1），说明北碚井附近含水层系统以垂向流为主、径向流为辅。根据前人钻孔资料，该井钻孔揭露本区29.9～70.24 m的长石砂岩含水层，井深30.40 m和36.32 m处涌水，水头高度分别为1.77 m和2.74 m，表明该含水层中地下水局部存在着自下而上的垂向流运动。原因为：一是据钻孔资料知井观测含水段裂隙较发育，裂隙倾角约55°～80°，判断局部发育的裂隙为含水层内地下水垂向运动提供了优势通道；二是井深30.40、43.30 m处的裂隙被方解石半充填，38.40 m处的裂隙被黄铁矿晶粒充填，形成局部相对弱透水层与含水层互层，观测段（42.1～70.24 m）可通过弱透水层与上覆含水层发生垂向水量交换。由于Hsieh模型适用于封闭性良好的承压含水层径向流（η＜0），Wang模型适用于存在垂向水量交换的越流含水层系统，故本次采用Wang模型计算北碚井附近含水层的水力特征参数，用以表征垂向上水量交换的能力。首先根据北碚井钻孔资料将观测含水层厚度M取值为40.34 m，利用潮汐M2波引潮高振幅和该井观测段的井水位数据，据公式（8）估算北碚井含水段的贮水能力，结果见表2，求得观测含水层的贮水率为1.41×10－6m－1，贮水系数为5.68×10－5。通过前人抽水试验［28］可知该井观测段含水层的水平渗透系数K＝0.014 9 m／d，则导水系数T＝KM＝0.60 m2／d，依据方程（14）和（15）可绘制出北碚井的越流系数和相位差的关系曲线（图6）。由表4可知北碚井水位潮汐分析的相位差变化范围约为25°～35°，则可估算越流系数的范围为1.5×10－8～2.5×10－8s－1。

表4 北碚井所处含水层系统越流系数计算Tab.4 Calculation results of the coefficient of leakage in the aquifer near well Beibei

图6 北碚井相位差与越流系数关系Fig.6 Relationship between the coefficient of leakage and phase shift in well Beibei

4 地震对含水层水力特征的影响

4.1 汶川地震引起的井水位变化

2008年5月12日14时28分，我国四川省汶川县映秀镇与卧龙镇之间发生了Ms 8.0级地震，地震发生在青藏高原东缘的松潘-甘孜地块与扬子地块交界的龙门山主中央断裂带上，该地震造成了我国大部分井水位发生同震响应［29］，其中包括大足和北碚井。如图7（a）所示，汶川地震发生时，大足井水位发生94 cm的阶降变化，2天后呈现正常潮汐波形，但数月内未恢复到震前水平。北碚井水位出现明显阶降现象，幅度为93.6 cm，3个小时后水位重新达到新的平稳状态，但数月内未恢复到震前状态，见图7（b）。

图7 汶川地震引起井水位变化Fig.7 Variation of well water level caused by the Wenchuan earthquake

4.2 汶川地震引起含水层水力特征的变化

汶川地震发生时大足和北碚井水位皆出现了明显的同震响应，利用井水位的潮汐效应定量分析地震前后含水层水力特征是否发生了变化。以大足井和北碚井2008年的水位小时值观测数据为基础，进行固体潮调和分析，为了准确地描述汶川地震前后含水层水力特征的变化，此次数据间隔设置为20 d，滑动窗长设置为10 d，计算地震前后大足和北碚井固体潮M2波的相位差。

如图8（a）所示，汶川地震前大足井的相位差在0°上下浮动，地震发生后相位差最高达到11°，根据相位差和地下水流动类型的关系判断，震后地下水发生了垂向流运动，表明含水层在垂向上存在着水量交换，由此引入越流系统模型来估计大足井附近含水层的垂向渗透能力，基于水平渗透系数以及含水层的贮水能力，根据公式（14）和（15），得到井水位固体潮响应相位差与含水层越流系数之间的关系（图9）以及越流系数的计算结果（图10），可见汶川地震发生前后大足井的水位相位差上升为正值，相应的越流系数增大，即垂向渗透性增强。

图8 汶川地震前后相位差变化Fig.8 Variation of phase shift before and after the Wenchuan earthquake

图9 大足井相位差与越流系数关系Fig.9 Relationship between the coefficient of leakage and phase shift in well Dazu

如图8（b）所示，汶川地震的发生造成北碚井的相位差由震前平均值26°上升到67°，上升幅度为41°，变化量超过3倍均方差，表明含水层内部地下水运动仍以垂向流为主，且垂向流运动增强。由此引入越流含水层系统模型推求北碚井附近含水层越流系数的变化，以表征含水层垂向水量交换的能力。计算结果见图10，汶川地震的发生导致含水层越流系数由震前平均值1.692×10－8s－1上升到2.541×10－6s－1，根据越流系数的定义可知，观测含水层与相邻含水层之间垂向水量交换增大。可见，北碚井和大足井的潮汐参数在汶川地震后均有所增大且大于0，表明两口井附近含水层垂向流运动加强。引入越流系统模型推求含水层水力特征参数的变化，结果显示两口井所在含水层越流系数皆增大，定量说明了观测含水层与相邻含水层的垂向水流交换增强。

图10 汶川地震前后含水层系统越流系数变化Fig.10 Variation of the coefficient of leakage in the aquifer system before and after the Wenchuan earthquake

4.3 含水层水力参数变化机制探讨

地震的发生会造成区域构造应力（静应力）的变化和地震波（动应力）的传播，导致含水层形变甚至引起含水层水力特征（渗透性）的改变。对于静应力引起的井水位变化的解释，孔隙线弹性理论［30］和松散含水层固结理论［31］比较成熟。根据Biot［30］孔隙线弹性理论，震源破裂产生的同震静应力引起地壳发生变形，导致含水层内部孔压发生变化，故表现为压缩区域井水位上升，膨胀区域井水位下降。固结理论认为，对于松散沉积物，剪应变引起颗粒运动进入到原有的孔隙中，导致孔隙体积减小，孔压增大，井水位发生变化。对于动应力对含水层造成的影响，一般认为地震波作用于井-含水层系统，导致裂隙的疏通［32］（裂隙中胶体颗粒、液滴和气泡的活化迁移）或堵塞，裂隙的张开或闭合［33］，甚至是新裂隙的形成［34］，引起含水层渗透性或层间水力联系的改变，从而造成含水层内部孔隙压力的重新调整与平衡，引起井水位的响应。对于地震引起井水位以及含水层水力特征变化的原因，可能并不是一种机制可以解释的，还需要结合井孔周边的水文地质条件以及井自身的状况进行具体分析。

汶川地震后，北碚井和大足井的水位发生了明显的阶降变化且地下水垂向流增强。一方面，据Shi等［35］得到汶川地震同震破裂体应变空间分布与井水位同震响应形态对比（图11），可知北碚井和大足井位于含水层膨胀区域，即汶川地震发生后，震源破裂产生的构造应力使北碚和大足井附近含水层发生了膨胀变形，孔压降低，井水位明显下降。另一方面，地震发生时，应力传播至含水层，加快了地下水流动速度，可能致使堵塞在含水层局部裂隙中的钙质、铁质胶体或细小颗粒被清除，从而造成大足井和北碚井附近含水层中裂隙疏通，连通性增大，且因为裂隙倾角近垂直，所以大足井表现出垂向流的特征，北碚井垂向交换量增强。值得注意的是，该区域主要观测含水层在较大范围内存在着统一的水力联系［26］，地震应力使得裂隙通道张开，各个井孔附近含水层间的水量交换加强，势必造成在岩层露头部位孔隙压力的聚集和释放，井水位迅速上升。如巴南井为该区一口自流井，两口井距离约80 km，该区地势较低，地下水沿砂岩裂隙通道出露地表，如图12所示，汶川地震的发生导致该井水位发生幅度约为120 cm的阶升现象，这也是区域地下水径流运动状态改变的结果。综上，该区域上沙溪庙组砂岩含水层分布面积广，且具有统一的水力联系，大足井和北碚井皆地处径流区，汶川地震的发生导致含水层发生膨胀变形和裂隙疏通，地震造成井水位下降及垂向水量交换增强，井水位变化形态基本一致，巴南井地势较低且为自流井，孔隙压力的聚集和释放导致井水位上升。

图11 汶川地震同震孔隙体应变［28］（膨胀为正）与井水位变化Fig.11 Co-seismic pore volume strain and well water level changes in the Wenchuan earthquake

图12 巴南井水位变化曲线Fig.12 Variation of well water level of well Banan

5 结 论

本文系统分析了井水位潮汐响应与含水层水力特征的理论关系，以华蓥山断裂带上大足井和北碚井为例，基于井水位的潮汐分析计算了观测段含水层水力特征参数（包括贮水系数、渗透系数、导水系数、越流系数），并探讨了汶川地震对含水层水力特征及地下水流运动状态的影响及机制。主要结论如下：

由潮汐分析知大足井附近含水层中地下水运动以水平径向流为主，北碚井孔-含水层系统中裂隙较发育，地下水垂向流和径向流并存；利用潮汐径向流模型计算出大足井含水层渗透系数为0.041 9 m／d，与以往抽水试验的结果（0.053 4 m／d）基本一致，表明利用井水位的潮汐响应计算含水层水力参数具有可行性，引入越流含水层系统模型计算得到北碚井越流系数范围为1.5×10－8～2.5×10－8s－1。

汶川地震引起北碚井和大足井的水位均发生了阶降变化，且含水层中地下水垂向流运动加强，原因是该区域侏罗系砂岩含水层在较大范围内存在着统一的水力联系，且局部裂隙发育，倾角较大，汶川地震的发生导致含水层发生膨胀变形和裂隙疏通，造成了径流区北碚井和大足井水位下降及垂向水量交换增强，排泄区井水位上升。