孙小龙,向 阳,2

(1. 中国地震局地壳应力研究所(地壳动力学重点实验室),北京 100085；2. 新疆维吾尔自治区地震局,新疆 乌鲁木齐 830002)

储水系数和渗透系数是含水层系统两个最重要的参数,前者反映含水层释放或储存水量的能力,后者表示水体通过孔隙骨架的难易程度。与储水系数相关的有储水率,表示单位体积的含水层在水头变化一个单位时释放(或存贮)的水量；与渗透系数相关的有导水系数,反映含水层全部厚度导水能力的参数,数值上等于渗透系数与含水层厚度的乘积。以上各参数通常基于室内实验、野外试验和模型反演获取：(1)通过实验室分析岩石或土壤样品的特定物理参数,如孔隙度、粒径、压缩系数等,确定其水文参数[1～2]；(2)通过野外抽水试验或锤击试验,基于Theis和 Hantush 模型,依据不同的含水层类型,通过试验中的水位变化估算以上各参数[3～5]；(3)通过井-含水层系统对外界刺激的响应特征,如固体潮变化[6]、气压响应[3,7]、同震响应[8～9]等,基于不同的含水层水流模型反演计算所需要的水文地质参数。以上三类方法中,方法(1)由于样品采集的随机性与实验的理想化,由此得到的结果与实际的参数可能存在一定的差别；方法(2)虽然最为普遍且实用,但其花费昂贵。因此,方法(3)常常是一个不错的选择,尤其是利用井水位对固体潮的响应特征反演含水层渗透率,近年来已在各类研究中广泛应用[11～12]。但是,并不是所有的观测井都对固体潮或气压有响应,因此该方法在使用时也会受到客观条件的限制。

Cooper 等研究了地表垂直运动引起的井-含水层系统响应特征,认为水位的响应是地表垂直运动和含水层孔隙压力波动共同引起,不同周期的信号引起的响应幅度也各不相同,而且响应幅度与含水层水文参数(导水系数、储水系数等)密切相关[4]。因此,可通过分析面波时段水位和地表运动的幅度比,确定各水文参数[9]。Shih 等通过频谱分析法,研究了井水位对地震面波的响应特征,并利用水位与地表垂直运动的频谱关系确定了井-含水层的储水系数[5]。本文通过分析新疆新10井水位对全球5次7级以上地震的同震响应特征,提出一种利用水震波信号(主要是瑞利波引起)反演获取承压含水层水文地质参数的方法。

1 新10井水文地质背景

新10井位于中国西部地区新疆乌鲁木齐市南部,海拔为1056 m。其位于柳树沟—红雁池逆冲断裂及其派生断裂的交汇部位(图1a),沿该断裂挤压破碎带广泛出露泉水,多数泉水流量常年稳定,涌水量达20～30 L/s。

新10井于1980年9月成井,其井孔结构如图1(b)所示。井深28 m,开孔孔径146 mm至9.2 m处,130 mm孔径9.2～13.04 m,13.04 m以下孔径为110 mm。井孔16～24 m岩芯破碎,为主要出水段,在该处设置滤水管。含水层为二叠系红雁池组泥质粉砂岩和砾岩,厚度约8 m,井水位埋深约1 m。2016年10月底,对其井水位观测仪器系统进行了改造,安置中国地震局地壳应力研究所研制生产的SWY-II型数字水位仪,分辨率为1 mm,监测时间间隔为1 s。

图1 新10井构造位置(a)和井孔柱状图(b)以及2016年11—12月全球7级地震分布Fig.1 Geological setting and location of the X10 well (black triangle) and epicenters of the five Mw≥7.0 earthquakes(beach balls) during November and December, 2016 (a); aquifer lithology and bore structure of the well (b)

新10井清晰记录了2016 年11—12月期间全球5次Mw≥7.0地震的水震波图像。为了对比分析研究该井水震波对地震波的响应特征,本文收集了新10井东北方向相距约10 km处的SMG台SLJ-100型强震计记录到的垂向速度波形数据。2016年12月9日1时38分(北京时间)所罗门群岛Mw7.8地震(S 10.665°,E 161.335°)引起的新10井水位同震响应和SMG台垂向速度波形对比,见图2。为了更加清楚地对比,利用S变换法[13]对二者进行了时频特征提取,见图2(c)和图2(d)。所罗门群岛Mw7.8地震引起的新10井水震波形态与SMG台地震波垂向分量的形态和频率基本相似,且可清晰地看到S波和面波。

从图2(c)和图2(d)的时频特征对比图还可看出,水震波和地震波均存在一个显著的频率段,即0.03～0.06 Hz(15～35 s),持续时间约20 min。面波时段,频率小、周期大的信号先出现,频率大、周期小的信号后出现。从能量值(红色表示能量高、蓝色表示能量小)的对比看,20 s左右周期信号水震波和地震波的能量最大,这也是面波(瑞利波)的主要周期成份。

图2 2016年12月9日所罗门Mw7.8级地震引起的新10井水震波(a)及其频谱图(c)、水磨沟地震台记录到的地震波垂直分量(b)及其频谱图(d)Fig.2 (a) Groundwater level changes in X10, (b) seismic wave in SMG of M7.8 Solomon Islands earthquake on December 9, 2016, (c) the time-frequency images of the groundwater level and (d) seismic wave

2 水文地质参数估算

水震波和地震波(垂向分量)在形态和频率上均高度相似,尤其在面波时段。本文通过分析新10井对全球5次Mw≥7.0地震的水震波特征,确定其井-含水层水文参数。

2.1 储水系数

储水系数,反映含水层水头下降或上升单位高度时,从单位水平面积和高度等于含水层厚度的柱体中释放或储存水体积能力的一个参数。水位对地表运动的响应幅度与储水系数、含水层厚度和面波波数等因素密切相关[8]：

(1)

式中：h(t)——水位；

w(t)——地表位移；

S——储水系数；

b——含水层厚度；

kx——地震波波数。

通过傅立叶变换,可提取出不同周期段信号的响应幅度[14]：

(2)

式中：Zh(ω) ——水位h(t)的傅立叶变换；

Zw(ω) ——地表位移w(t)的傅立叶变换；

ω——频率；

t——时间；

图3所示为2016年12月9日1时38分(北京时间)所罗门群岛Mw7.8地震引起的水位、垂向速度的时序图和频谱图,可以看出,相比其它周期的波,20～35 s周期波的响应幅度较大。

图3 面波时段的水位(a)及其频谱(b)和垂向速度(c)及其频谱(d)Fig.3 (a) Groundwater level and (b) its frequency, and (c) the vertical velocity and (d) its frequency at the surface wave period

利用水位和垂向位移的功率谱密度关系,可以求得井-含水层系统的储水系数：

(3)

(4)

(5)

式中：λ——波长,可以由地震波波速和周期确定；

Shh和Sww——水位h(t)和地面运动垂向位移w(t)的自功率谱密度(PSD)；

Swh——h(t)和w(t) 的互功率谱密度(CSD)。

实际上,对于一组有限长度Tr的记录数据,可表示为[15]：

(6)

其中,数学期望E表示某个采样周期段内的平均值。

图4所示为2016年12月9日1时38分(北京时间)所罗门群岛Mw7.8地震引起的水位、垂向速度的自功率谱密度(a)和互功率谱密度(b)图。可以看出,在0.03～0.04 Hz频段内二者的功率谱密度最大,且一致性较高,可用这一显著频段的信号来估算储水系数[16]。

图4 水位、垂向速度的自功率谱密度(a)和互功率谱密度(b)Fig.4 Power spectral densities (a) and cross-power spectrum densities of the groundwater level and vertical velocity (b)

新10井含水层在埋深16～24 m处(图1b),则含水层厚度d为8 m,计算波长λ=v/ω所需的波速v值由震中距除以面波到时与发震时刻之差获得,频率ω取决于图4中的功率谱密度的最高值。地表运动(地震波垂向分量)的质点位移幅度由速度幅度转化得到,即smax=vmax·τ/2π。依据式(3)～式(6),利用2016年12月9日1时38分(北京时间)所罗门群岛Mw7.8地震引起的水位、垂向速度的功率谱密度关系,估算的新10井含水层平均储水系数为1.89×10-5、储水率为2.36×10-6/m。

2.2 渗透系数

地震波引起的水震波,其响应幅度与井孔条件、含水层参数、地震波周期等密切相关。地震波作用除了会引起地表介质的上下波动之外,还会引起含水层内孔隙压力的波动变化。井水位的振荡由地表垂向运动和孔隙压力变化共同造成,且井水位对二者均有一定的放大作用。因此,水震波响应幅度既受到井-含水层系统内的孔隙压力波动影响,也受到地表垂向运动的影响。在承压性的井-含水层系统中,井水位对孔隙压力和地表垂向运动的响应可用放大因子(或幅度响应比)表示[4]：

(7)

式中：A——水位对含水层孔隙压力波动的放大因子；

A′——水位对地表垂向位移的放大因子；

rw——井孔内径；

He——观测井有效水柱高度；

T——导水系数；

Τ——地震波周期；

ω——地震波频率；

g——重力加速度；

Kei和Ker——零阶开尔文函数。

如果含水层孔隙压力波动引起的水位幅度为h、地表垂向运动引起的水位幅度为s,那么水位整体波动幅度为x=h+s。令R为含水层孔隙压力波动引起的水位变化幅度与地表垂向运动引起的水位变化幅度之比、M为水位整体波动幅度与地表垂向位移幅度之比：

(8)

式中：Ew——水的弹性模量/(22×108N·m-2)；

ρ——密度/(103kg·m-3)；

g——重力加速度，取9.81 m/s2；

v和τ——地震波的波速和周期；

n——含水层孔隙度,可由气压系数、储水率求得[3]。

依据式(7)和式(8),利用实测的水震波幅度和地震波幅度,可计算得到相应的R和M。

如果令放大因子比A/A′ =R′,那么,幅度比R与地震波周期τ成反比关系,而放大因子比R’则与周期τ成二次线性正比关系。R实测值计算所需的波速v值由震中距除以面波到时与发震时刻之差获得,周期τ值可依据频谱图(图3)读取；R′值计算需知道有效水柱高度He,可依据He=H+3d/8计算得到[4],H为隔水顶板以上的井内水柱高度。隔水顶板埋深16 m,水位埋深1 m,井内水柱高15 m,那么He=H+3d/8≈18 m。

依据频谱图可逐一统计得到各周期τ的M值,依据其R值可计算得到各自的A及A′,结果见图5。图5中实测值由图3频谱幅度比计算所得,理论曲线是利用频谱幅度比拟合所得,依据式(7)计算出与实际统计值最接近的渗透系数K值。计算中井径rw=54 mm,含水层厚度d=8 m,有效水柱高度He=18 m,储水系数由上述方法获得S=1.89×10-5。最终,渗透系数K值拟合结果为214 m/d。

图5 新10井水震波放大因子实测值与估算值 (a) 水位对孔隙压波动的放大因子; (b)水位对地表垂直运动的放大因子Fig.5 Observed (error circles) and estimated (solid line) amplification factors based on the groundwater waves in the X10 well (a) due to pore pressure fluctuations, and (b) due to vertical surface motions

3 结果与讨论

与所罗门群岛7.8级地震引起的同震响应类似,新10井水位同样清晰记录到了图1(a)中其它4次地震所引起的同震响应,它们分别是2016年11月13日新西兰7.8级地震 (S 42.757°, E 173.077°)、2016年11月25日中美洲沿岸远海7.0级地震 (N 11.960°, W 88.836°)、2016年12月17日新爱尔兰地区7.9级地震(S 4.509°, E 153.450°) 和2016年12月25日智利7.6级地震(S 43.416°, W 73.951°)。

基于以上方法,利用其它4次地震引起的新10井水震波和SMG台垂直位移波形数据,估算井-含水层系统参数,见表1。虽然5次地震的震级、震中距均不相同,尤其是震中距,最小为8 451 km,最大为18 537 km(接近地球周长的一半),除2016年12月9日估值偏高外,其余参数基本相当。储水系数S值为1.33×10-5～1.64×10-5(相对应的储水率Ss=S/d在1.67×10-6～2.05×10-6/m之间),渗透系数K值其余值在57～87 m/d之间(相对应的导水系数T=K·d在456～696 m2/d之间)。基于2016年12月9日所罗门群岛地震的水震波估算的含水层参数值偏高的原因,主要是由于在此次地震发生前的12 h,距新10井约110 km处发生了呼图壁6.2级地震,这种近场地震造成的区域应力调整,可导致新10井含水层系统介质参数发生短时间内的明显变化,如孔隙度、储水系数、渗透系数等增大,进而表现为水位的阶变上升。随着区域应力调整的结束,各参数也会随即恢复正常,如2016年12月17日的估值和其它三次地震的估值基本一致。

表1 新10井-含水层系统水文地质参数估算结果Table 1 Estimated hydraulic parameters of the X10 well-aquifer system

3.1 方法的适用性

地下水位是地震监测预报中重要的观测项目之一,新10井即属于这类观测,如果要获取这类井的含水层参数,只能通过模型反演实现。因为,考虑到这类观测井的水位观测数据必须要连续,抽水试验不能实现。模型反演获取水文参数,常用的是水位固体潮调和分析法[6],即利用水位对潮汐的响应幅度和相位滞后来估算其导水系数、储水系数等。但是,由于井-含水层系统条件的差异,并不是所有观测井的水位都能反映出明显的固体潮汐变化,有些井水位根本就没有潮汐形态。

如图6所示,新10井水位由于气压对其影响较大,导致其水位固体潮汐响应看上去很微弱。新10井水位,不论是长期还是短期,均明显受到气压的影响,水位和气压呈负相关关系(相关系数R=-0.85),其气压系数[17]为5.2 mm/hPa。但是,水位的固体潮汐形态却很微弱,虽然在有些时候能看到(11月20日—12月1日),但相比气压的影响幅度,其变化幅度较小。因此,如果利用新10井水位的固体潮汐变化反演井-含水层水文参数,由于数据误差较大其结果的可信度会很低。另外,由于气压与固体潮在日波、半日波的某些频率段其信号周期是重合的,这会导致在固体潮调和分析中必然存在某些误差。

图6 2016年11—12月新10井水位与气压分钟值对比Fig.6 Time series comparisons of the groundwater level (blue) and barometric pressure (green) at the X10 well in November and December, 2016

因此,对于类似于新10井无清晰固体潮形态、无抽水试验条件的地震信息观测井,利用水位对面波(瑞利波)的同震响应来反演其井-含水层水文参数,是一种不错的选择。由于远场地震引起的水震波只能是地震波作用所引起,其井水位对地震波的响应幅度与井孔条件、地震波特征、含水层参数等密切相关。在无其它作用源的情况下,井孔条件一般不会发生变化,与地震波特征相关的水震波形态,其主要影响因素就是含水层参数。因此,在井孔条件明确的情况下,基于水震波和地震波联合反演得到的含水层参数,其结果的可信度也相对较高。

3.2 方法的局限性

从时频特征分析看,水震波的周期与地震波(瑞利波)基本相当,多在20 s左右。因此,利用水震波和地震波联合反演含水层参数时,水位的监测时间间隔应小于水震波的周期,本文所用数据为秒值。目前,全国各地震水位观测井的数据监测时间间隔多为分钟值或整点值,无法满足本文所提出的参数反演方法需求。因此,推广地震地下水位的秒值采样,可更清晰完整地记录到远场地震所引起的水震波,为同震响应机理研究、含水层参数求取提供可靠的数据保障。

地震波作用除了引起地表介质的上下波动之外,还会引起含水层内孔隙压力的波动变化。井水位的振荡由地表垂向运动和孔隙压力变化共同造成。在对比分析新10井水震波与地表垂向运动时,地震波形数据选用了离新10井约10 km处的水磨沟台垂向速度数据。虽然相距10 km波形数据的差异会很小,但毕竟由于地下介质的细微差异,会导致与实际值之间存在微小的误差。因此,理想情况下,利用水震波和地震波反演含水层参数时,应选用与水位观测井同一位置处的波形数据,这样反演得到的参数值会更加准确。

另外,在用本文探讨的方法反演含水层参数时,应尽可能与其它方法(如室内实验[18]、抽水试验[19]、水位固体潮反演等)得到的参数进行对比分析,这样得到的结果可信度将更高。笔者曾用该方法计算了北京昌平井的含水层参数[9],计算结果与利用水位固体潮反演得到的参数基本一致,表明利用水震波和地震波垂向速度分量联合反演得到的含水层参数具有较好的可信度。

4 结论

本文基于新10井水位对远场地震的同震响应,通过频谱法分析了该井水震波与地震波的相关性,提出了一种利用水震波信号(主要是瑞利波引起)反演获取承压含水层水文地质参数的方法。得出以下结论：

(1)新10井水震波与地震波在变化形态、频谱特征上均具有很好的一致性。

(2)通过频谱关系分析,估算得到了新10井-含水层系统的水文地质参数,其储水系数S值在1.33×10-5～1.64×10-5之间,渗透系数K值在57～87 m/d之间。

(3)基于水震波和地震波信号估算承压井-含水层系统水文参数,可为地下水管理、地震前兆异常分析等提供一种新的技术方法。

本文探讨的方法仅应用了少数几次地震所记录的数据,虽然方法本身具有较好的应用价值,但仍需积累更多震例、更多观测井数据来加以完善或改进。如果能收集到观测井自身或附近区域的抽水试验资料,或者是观测井水位具有较清晰的固体潮响应形态,就可以利用抽水试验计算或水位固体潮反演得到的含水层参数来印证本文探讨方法得到的参数值的可靠性,这也是笔者下一步的研究方向。

致谢：新疆维吾尔自治区地震局监测中心为本项研究提供了SMG台地震波形数据,在此表示衷心的感谢！

参考文献：

[1] Chapuis R P. Similarity of internal stability criteria for granular soils[J]. Canadian Geotechnical Journal, 1992,29(4): 711-713.

[2] Ma L, Xu Y S, Shen S L,etal. Evaluation of the hydraulic conductivity of aquifers with piles[J]. Hydrogeology Journal, 2014,22(2):371-382.

[3] Jacob C E. On the flow of water in an elastic artesian aquifer [J]. Eos, Transactions American Geophysical Union, 1940,21(2): 574-586.

[4] Cooper H H, Jacob C E. A generalized graphical method for evaluating formation constants and summarizing well‐field history[J]. Eos, Transactions American Geophysical Union, 1946,27(4):526-534.

[5] Sethi R. A dual-well step drawdown method for the estimation of linear and non-linear flow parameters and wellbore skin factor in confined aquifer systems[J]. Journal of Hydrology, 2011,400(1/2):187-194.

[6] Hsieh P A, Bredehoeft J D, Farr J M. Determination of aquifer transmissivity from earth tide analysis[J]. Water resources research, 1987,23(10):1824-1832.

[7] Ritzi R W, Sorooshian S, Hsieh P A. The estimation of fluid flow properties from the response of water levels in wells to the combined atmospheric and earth tide forces[J]. Water Resources Research, 1991, 27(5): 883-893.

[8] Shih D C F. Storage in confined aquifer: Spectral analysis of groundwater responses to seismic Rayleigh waves[J]. Journal of Hydrology, 2009,374(1): 83-91.

[9] Sun X, Wang G, Yang X. Coseismic response of water level in Changping well, China, to the Mw 9.0 Tohoku earthquake[J]. Journal of Hydrology, 2015,531: 1028-1039.

[10] Xue L, Li H B, Brodsky E E,etal. Continuous permeability measurements record healing inside the Wenchuan earthquake fault zone[J]. Science, 2013,340(6140): 1555-1559.

[11] Wang C Y, Liao X, Wang L P,etal. Large earthquakes create vertical permeability by breaching aquitards[J]. Water Resources Research, 2016,52(8): 5923-5937.

[12] Shi Z, Wang G. Aquifers switched from confined to semiconfined by earthquakes[J]. Geophysical Research Letters, 2016,43:1-7.

[13] Stockwell R G, Mansinha L, Lowe R P. Localization of the complex spectrum: the S transform[J]. IEEE Transactions on Signal Processing, 1996,44(4): 998-1001.

[14] Duhamel P, Vetterli M. Fast Fourier transforms: a tutorial review and a state of the art[J]. Signal Processing, 1990,19(4): 259-299.

[15] Bendat J S, Piersol AG. Random Data, Analysis and Measurement Procedures[M]. New York: John Wiley and Sons, 2000.

[16] Shih D C F, Wu Y M, Chang C H. Significant coherence for groundwater and Rayleigh waves: Evidence in spectral response of groundwater level in Taiwan using 2011 Tohoku earthquake, Japan[J]. Journal of Hydrology, 2013,486(1): 57-70.

[17] Clark W E. Computing the barometric efficiency of a well[J]. Journal of the Hydraulics Division, 1967,93(4): 93-98.

[18] 葛勤, 梁杏, 龚绪龙,等. 低渗透岩土有效扩散系数的室内测定与分析[J]. 水文地质工程地质, 2017, 44(3):93-99. [GE Q, LIANG X, GONG X L,etal. Laboratory determination and analysis of dffective diffusion coefficients for low-permeability rock and clay[J]. Hydrogeology & Engineering Geology, 2014, 44(3):93-99.(in Chinese)]

[19] 周志芳, 徐海洋. 一种实验确定弱透水层水文地质参数的原理与方法[J]. 水文地质工程地质, 2014, 41(5):1-4.[ZHOU Z F, XU H Y. Theory and methods for determining hydrogeological parameters of an aquitard based on experimental data[J]. Hydrogeology & Engineering Geology, 2014, 41(5):1-4.(in Chinese)]