2014年8月3日16时30分云南鲁甸龙头山镇发生M_S6.5强震,震中位置(27.11°N,103.35°E),震源深度12 km。本次地震发生在巴颜喀拉块体、川滇块体和华南块体交汇区南部,马边断裂南段、昭通断裂、莲峰断裂之间。由于本次地震造成的断层破裂长度有限且未达到地表,因此破裂面方位不易识别。震后许多学者及机构通过地震烈度分布和地震破裂过程反演对鲁甸地震的发震断层进行了研究,得出鲁甸地震发震构造为一条新发现的近北西走向的断层——包谷垴—小河断裂(房立华等,2014; 李西等,2014; 徐锡伟等,2014; 李艳娥等,2015)。但又有学者通过视震源时间函数分析、余震分布及震源机制解研究提出鲁甸地震可能是由近东西向和近南北向的共轭断层先后破裂导致(张勇等,2015)。

本文研究参照中国地震局对于鲁甸地震的初步分析报告(Cheng et al,2015),鲁甸地震为一次高倾角左旋走滑地震,发震断层为走向340°的包谷垴—小河断裂(图1)。

图1 鲁甸地震震中及其附近地区发震构造图
Fig.1 Seismogenic structures in the 2014 Ludian MS6.5 earthquake region

云南会泽井(滇01号井)为静水位观测井,地理坐标为(26.52°N,103.15°E)。井孔地处川滇地块东南部,位于小江、莲峰、则木河3大断裂的交汇部位,构造位置特殊,井孔周边水文地质条件如图2a所示。会泽井受白雾街不对称双曲弧形构造控制,张性裂隙发育,为地下水的富集提供了有利条件,井孔周边致密玄武岩柱状节理发育,易于形成地下水的运移通道,加上两侧岩溶水的侧向补给以及大气降水的天然补给,使得井孔所处含水层具有较好的补、径、排条件,能够形成一定规模的含水系统。会泽井深103.15 m,井孔半径54 mm,套管深度87.8 m,水位埋深约30 m,滤水管为34.06～87.80 m,观测段为34.06～103.15 m,观测层岩性主要为第四系风化玄武岩(图2b),地下水类型为裂隙承压水,水温约16 ℃(Sun et al,2019)。

会泽井自2012年开始数字化观测,水位对其周边几次地震均有响应,映震性较好(图3a)。多年来,水位整体呈缓慢降低型,水位下降速率约为0.3 m/a,大气降水对水位影响较小,气压效应不明显,有固体潮效应。在2014年鲁甸M_S6.5地震时,距离震中约70 km的会泽井水位有同震响应,水位在地震发生时瞬时上升0.33 m,之后开始下降,并在震后50 d内逐步恢复至震前水平(图3b)。

图2 会泽井水文地质简图(a)和井孔柱状图(b)
Fig.2 Hydrogeological map of the region in which Huize Well located(a)and the borehole histogram of Huize Well(b)

图3 2012—2014年会泽井水位变化(a)和2014年鲁甸MS6.5地震会泽井水位的同震响应(b)
Fig.3 Variation of the water level in Huize Well from 2012 to 2014(a)and the coseismic variation of the water level in Huize Well(b)

含水层中孔压受地震静态应力的影响,在压力梯度的驱动下扩散恢复,这个过程中含水层的弹性变形和孔隙水压的扩散是一个流固耦合的过程,因此,需要建立含水层应变与孔压之间的耦合数学模型来定量解释二者之间的关系。本文对含水层孔压扰动及扩散控制方程的建立和求解作如下假设:①含水层为均质各向同性、连续的承压含水层,顶、底板水平; ②多孔介质骨架为各向同性弹性体,骨架变形为小变形,符合线性孔弹性理论; ③地下水流符合达西定律; ④地下水流为等温渗流; ⑤渗透率各向同性; ⑥单元体内流体无外部源汇项。

本文只关注静态应力对含水层孔压的扰动,忽略加速度的影响,在含水层中选取典型单元体并进行受力分析。利用静力平衡方程、广义胡克定律、有效应力原理、达西定律与质量守恒定律,建立含水层中任一单元体在受到静态应力作用下应变与孔压的耦合控制方程组,具体推导过程见魏海滨(2021)的研究:

式中:λ、G为孔弹性介质的拉梅参数; u_x、v_y、w_z分别表示位移在x、y、z方向上的导数; α为有效应力系数; p为孔隙水压; k为含水层渗透率; μ为流体动力粘滞系数,Ss₃'为修改的三维贮水率;

为了最大限度地减少模型边界对孔隙压力扩散的影响,将模型范围定为200 km(东西向)×200 km(南北向)×0.5 km(深度)(图5),此时模型边界处孔压扰动值很小可忽略不计。初始时刻,系统各项位移均为0,含水层孔压为地震静态应力导致的孔压扰动值,即初始条件为:u=v=w=0,Δp=-B/3Δσ_kk(其中Δσ_kk=2G(1+v)/(1-2v)Δε_kk)。研究区地表允许自由变形,侧边界及底边界固定位移; 含水层顶板为黏土层,底板为致密的玄武岩,故顶、底板为零渗透边界,侧边界处孔压不受扰动即为0。

图5 含水层孔压扰动及扩散模型范围及震中、井孔分布图
Fig.5 Disturbance of the pore pressure in the confined aquifer,the range of diffusion model,the epicenter and the borehole loacation

在井-含水层系统中,当地下水在压力梯度的驱动下由含水层向井中流动时,水流方向由纵向变为垂向,且井中为纯流体域。因此,地下水在井中的运动有别于水流在含水层中的运动,在井-含水层系统中孔隙水压扰动从压力扰动源到井中水位变幅,应分为含水层中压力扩散与井中水柱运动两部分。为了刻画在静态应力的影响下孔压在含水层中扩散并在井中以水柱的形式运动这一过程,本文建立了以井为中心的二维轴对称的小尺度模型(图6a)。

该模型包含承压含水层与井两部分,井孔半径54 mm,含水层概化为厚100 m、半径100 m。在侧边界施加地震静态应力导致的孔压变化值p,边界处孔压变化值的作用使得含水层内产生孔压梯度。在孔压梯度的驱动下,含水层中的水流向井或井水中流向含水层运动(井孔只在侧壁进水流速为u)(图6b)。当地下水在承压含水层中运动时,假设含水层均质各向同性,顶、底板水平,水流符合达西定律:

式中:S_s为贮水率; K_r为径向渗透系数。

当地下水进入井孔后,水流由径向流动变为垂向流动,且井孔为纯流体域,故采用Navier-Stokes方程描述井孔中水流的运动:

式中:V_z为井孔中水流速度; g为重力加速度; p为孔压; ρ和μ分别为水的密度和动力粘滞系数。

在初始时刻,含水层中各点孔压为静水压力,即p₀=ρgh; 边界条件为:p=ρgh+Δp。

图6 以井为中心的二维轴对称小尺度模型(a)及水位同震变化模拟过程(b)
Fig.6 Two-dimensional axisymmetric model of the water level response of the aquifer-well system(a) and simulation process of water level co-seismic response(b)

将鲁甸地震静态应变场转化为含水层孔压扰动值,以此作为初始条件,利用COMSOL Multiphysics中的PDE模块,在200 km×200 km×0.5 km的研究区范围内对建立的含水层孔压扰动及扩散数学模型进行求解,获得研究区0.1 km深度处含水层孔压演化规律。最后,在COMSOL中利用内置的达西定律、层流模块,建立以井孔为中心、半径100 m的小尺度井-含水层系统模型,以静水压强为初始条件,利用含水层孔压扰动及扩散模型中会泽井孔所处剖分网格的平均孔压扰动值作为边界条件,对会泽井水位的变化趋势进行模拟。

3.3.1 含水层孔压扰动及扩散数值模拟

基于研究区概念模型(图5)构建三维几何模型用于数值模拟,将式(1)表示为广义矩阵形式,输入COMSOL Multiphysics的PDE模块:

其中,

研究区浅层含水层岩性主要为二叠系峨眉山组风化玄武岩,本文利用等效连续介质理论,将裂隙介质等效为多孔介质,参照水文地质手册经验值对模型进行赋值,参数设置见表2。结合模型初始条件和边界条件,确定可能影响含水层孔压扩散的因素,设计不同情景模拟鲁甸地震后含水层孔压扩散的过程,模拟情景设置见表3。采用自由四面体网格剖分计算网格,考虑计算机的算力限制及计算效率,在保证模型计算精度的前提下共剖分5 417 817个网格,模拟时间步长为1 d,模拟时段为研究区0.1 km深度孔压扰动值消散99%以上所需时长。

表2 数值模拟参数
Tab.2 Parameters for numerical simulation

表3 不同模拟情景参数
Tab.3 Parameters under different operating conditions

以鲁甸地震同震静态应变场作为初始条件,对不同模拟情景下研究区含水层孔压扰动及扩散过程进行模拟。初始时刻,研究区孔压分布与同震静态应变场分布一致,但符号相反,最大正、负超孔隙水压变化为分别3.52×10⁵ Pa、-3.74×10⁵ Pa。研究区内距离断层两端较远的区域产生的同震孔压扰动值在±(10²～10³)Pa,相较断层两端附近的孔压扰动值小2～3个数量级。

情景1、2模拟不同渗透系数条件下含水层孔压演化规律,如图7所示。模拟结果表明,断层两端的孔压梯度决定了孔压的扩散模式。地下水沿着断层东西两侧的压力最大梯度流动。情景1中含水层孔压经过50 d的演化,孔压变化的最大值已降低了99%以上,孔压基本恢复至震前水平。情景2中鲁甸地震造成的孔压扰动值在低渗透系数模型中经过500 d消散99%以上,基本恢复至震前水平。

图7 在情景1(a)和情景2(b)下模拟不同时间研究区0.1 km深度孔压分布
Fig.7 Initial pore pressure in the beginning,and the simulated pore pressure at 0.1 km depth under Condition 1(a)and under Condition 2(b)at different times

通过比较情景1、2模拟结果可知,两个模型的孔压消散模式类似,但当渗透系数降低一个量级时,较低的渗透系数使孔压消散的时间增加了约10倍。因此,渗透系数是地震后孔压扩散时长的主控因素,研究区地质介质类型是决定地震后孔压扩散速率的主要变量。

情景3、4分别将孔隙度设置为0.2、0.6,其它条件与情景1保持一致,模拟不同孔隙度对研究区孔压扩散规律的影响。模拟结果显示(图8),不同孔隙度并不会影响孔压消散模式,在相同时间内不同的孔隙度只会对孔压的消散速率有所影响但影响有限,孔隙度越大孔压消散越慢,不同孔隙度下含水层孔压基本都在50 d内恢复至震前水平。

情景5、6分别将杨氏模量设置为8×10⁸和8×10¹⁰,其它条件与情景1保持一致,模拟不同风化程度的玄武岩含水层对研究区孔压扩散规律的影响。模拟结果显示(图9),杨氏模量会显著改变研究区孔压扩散的时间,在其余相同的条件下,杨氏模量越大,孔压扩散时间越短,风化严重的玄武岩孔压会在350 d内恢复至震前水平,而未风化的玄武岩孔压会在20 d内恢复至震前水平。

实际统计结果显示,研究区内多口监测井水位基本均在50 d内恢复至震前稳定值,因此,情景1、3、4模拟得到的孔压恢复时间较为吻合。由模拟结果可知,含水层渗透系数与杨氏模量会显著改变孔压扩散时长,而孔隙度对含水层孔压扩散时长影响较小。

图8 在情景3(a)和情景4(b)下模拟不同时间研究区0.1 km深度孔压分布
Fig.8 Simulated pore pressure at 0.1 km depth under Condition 3(a)and under Condition 4(b)at different times

图9 在情景5(a)和情景6(b)下模拟不同时间研究区0.1 km深度孔压分布
Fig.9 Simulated pore pressure at 0.1 km depth under Condition 5(a)and under Condition 6(b)at different times

3.3.2 压力扰动下含水层-井系统水位响应数值模拟

获得研究区含水层内孔压演化过程后,利用COMSOL Multiphysics中内置的达西定律及层流模块,以静水压强为初始条件,选择情景1中会泽井孔处网格单元含水层平均孔压演化曲线,作为压力扰动下含水层-井系统水位响应模型的边界条件,模拟会泽井水位对鲁甸地震同震响应过程及震后恢复过程,结果如图 10所示。

图 10 会泽井同震水位变化模拟及实测对比图
Fig.10 Simulated and the actual coseismic variations of the water level in Huize well

模拟结果显示,会泽井水位先呈阶跃上升形式,之后在50 d内恢复至震前水平,井水位最大振幅为0.45 m。模拟与实测水位变化趋势一致,但在水位最大振幅的幅值上有所差异,实测水位最大振幅为0.33 m,较模拟值偏小,可能由于水位同震响应实际上升过程中有井损的产生所导致。震后井水位实测值下降速率较缓,并且在下降期间还存在几次小幅度的上升,查询当地气象资料,发现研究区在8月和9月降雨较多,均匀的降水使井水位恢复曲线较为平缓,但几次暴雨可能使实测水位小幅上升,但也不排除地震或其它构造应力的影响。本文对井水位恢复过程的模拟,未考虑降雨及其他构造应力的影响,因此,模拟与实测水位恢复曲线有所偏差。