孫小龍,向 陽,2

(1. 中國(guó)地震局地殼應(yīng)力研究所(地殼動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室),北京 100085；2. 新疆維吾爾自治區(qū)地震局,新疆 烏魯木齊 830002)

儲(chǔ)水系數(shù)和滲透系數(shù)是含水層系統(tǒng)兩個(gè)最重要的參數(shù),前者反映含水層釋放或儲(chǔ)存水量的能力,后者表示水體通過孔隙骨架的難易程度。與儲(chǔ)水系數(shù)相關(guān)的有儲(chǔ)水率,表示單位體積的含水層在水頭變化一個(gè)單位時(shí)釋放(或存貯)的水量；與滲透系數(shù)相關(guān)的有導(dǎo)水系數(shù),反映含水層全部厚度導(dǎo)水能力的參數(shù),數(shù)值上等于滲透系數(shù)與含水層厚度的乘積。以上各參數(shù)通?；谑覂?nèi)實(shí)驗(yàn)、野外試驗(yàn)和模型反演獲?。?1)通過實(shí)驗(yàn)室分析巖石或土壤樣品的特定物理參數(shù),如孔隙度、粒徑、壓縮系數(shù)等,確定其水文參數(shù)[1～2]；(2)通過野外抽水試驗(yàn)或錘擊試驗(yàn),基于Theis和 Hantush 模型,依據(jù)不同的含水層類型,通過試驗(yàn)中的水位變化估算以上各參數(shù)[3～5]；(3)通過井-含水層系統(tǒng)對(duì)外界刺激的響應(yīng)特征,如固體潮變化[6]、氣壓響應(yīng)[3,7]、同震響應(yīng)[8～9]等,基于不同的含水層水流模型反演計(jì)算所需要的水文地質(zhì)參數(shù)。以上三類方法中,方法(1)由于樣品采集的隨機(jī)性與實(shí)驗(yàn)的理想化,由此得到的結(jié)果與實(shí)際的參數(shù)可能存在一定的差別；方法(2)雖然最為普遍且實(shí)用,但其花費(fèi)昂貴。因此,方法(3)常常是一個(gè)不錯(cuò)的選擇,尤其是利用井水位對(duì)固體潮的響應(yīng)特征反演含水層滲透率,近年來已在各類研究中廣泛應(yīng)用[11～12]。但是,并不是所有的觀測(cè)井都對(duì)固體潮或氣壓有響應(yīng),因此該方法在使用時(shí)也會(huì)受到客觀條件的限制。

Cooper 等研究了地表垂直運(yùn)動(dòng)引起的井-含水層系統(tǒng)響應(yīng)特征,認(rèn)為水位的響應(yīng)是地表垂直運(yùn)動(dòng)和含水層孔隙壓力波動(dòng)共同引起,不同周期的信號(hào)引起的響應(yīng)幅度也各不相同,而且響應(yīng)幅度與含水層水文參數(shù)(導(dǎo)水系數(shù)、儲(chǔ)水系數(shù)等)密切相關(guān)[4]。因此,可通過分析面波時(shí)段水位和地表運(yùn)動(dòng)的幅度比,確定各水文參數(shù)[9]。Shih 等通過頻譜分析法,研究了井水位對(duì)地震面波的響應(yīng)特征,并利用水位與地表垂直運(yùn)動(dòng)的頻譜關(guān)系確定了井-含水層的儲(chǔ)水系數(shù)[5]。本文通過分析新疆新10井水位對(duì)全球5次7級(jí)以上地震的同震響應(yīng)特征,提出一種利用水震波信號(hào)(主要是瑞利波引起)反演獲取承壓含水層水文地質(zhì)參數(shù)的方法。

1 新10井水文地質(zhì)背景

新10井位于中國(guó)西部地區(qū)新疆烏魯木齊市南部,海拔為1056 m。其位于柳樹溝—紅雁池逆沖斷裂及其派生斷裂的交匯部位(圖1a),沿該斷裂擠壓破碎帶廣泛出露泉水,多數(shù)泉水流量常年穩(wěn)定,涌水量達(dá)20～30 L/s。

新10井于1980年9月成井,其井孔結(jié)構(gòu)如圖1(b)所示。井深28 m,開孔孔徑146 mm至9.2 m處,130 mm孔徑9.2～13.04 m,13.04 m以下孔徑為110 mm。井孔16～24 m巖芯破碎,為主要出水段,在該處設(shè)置濾水管。含水層為二疊系紅雁池組泥質(zhì)粉砂巖和礫巖,厚度約8 m,井水位埋深約1 m。2016年10月底,對(duì)其井水位觀測(cè)儀器系統(tǒng)進(jìn)行了改造,安置中國(guó)地震局地殼應(yīng)力研究所研制生產(chǎn)的SWY-II型數(shù)字水位儀,分辨率為1 mm,監(jiān)測(cè)時(shí)間間隔為1 s。

圖1 新10井構(gòu)造位置(a)和井孔柱狀圖(b)以及2016年11—12月全球7級(jí)地震分布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地震的水震波圖像。為了對(duì)比分析研究該井水震波對(duì)地震波的響應(yīng)特征,本文收集了新10井東北方向相距約10 km處的SMG臺(tái)SLJ-100型強(qiáng)震計(jì)記錄到的垂向速度波形數(shù)據(jù)。2016年12月9日1時(shí)38分(北京時(shí)間)所羅門群島Mw7.8地震(S 10.665°,E 161.335°)引起的新10井水位同震響應(yīng)和SMG臺(tái)垂向速度波形對(duì)比,見圖2。為了更加清楚地對(duì)比,利用S變換法[13]對(duì)二者進(jìn)行了時(shí)頻特征提取,見圖2(c)和圖2(d)。所羅門群島Mw7.8地震引起的新10井水震波形態(tài)與SMG臺(tái)地震波垂向分量的形態(tài)和頻率基本相似,且可清晰地看到S波和面波。

從圖2(c)和圖2(d)的時(shí)頻特征對(duì)比圖還可看出,水震波和地震波均存在一個(gè)顯著的頻率段,即0.03～0.06 Hz(15～35 s),持續(xù)時(shí)間約20 min。面波時(shí)段,頻率小、周期大的信號(hào)先出現(xiàn),頻率大、周期小的信號(hào)后出現(xiàn)。從能量值(紅色表示能量高、藍(lán)色表示能量小)的對(duì)比看,20 s左右周期信號(hào)水震波和地震波的能量最大,這也是面波(瑞利波)的主要周期成份。

圖2 2016年12月9日所羅門Mw7.8級(jí)地震引起的新10井水震波(a)及其頻譜圖(c)、水磨溝地震臺(tái)記錄到的地震波垂直分量(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 水文地質(zhì)參數(shù)估算

水震波和地震波(垂向分量)在形態(tài)和頻率上均高度相似,尤其在面波時(shí)段。本文通過分析新10井對(duì)全球5次Mw≥7.0地震的水震波特征,確定其井-含水層水文參數(shù)。

2.1 儲(chǔ)水系數(shù)

儲(chǔ)水系數(shù),反映含水層水頭下降或上升單位高度時(shí),從單位水平面積和高度等于含水層厚度的柱體中釋放或儲(chǔ)存水體積能力的一個(gè)參數(shù)。水位對(duì)地表運(yùn)動(dòng)的響應(yīng)幅度與儲(chǔ)水系數(shù)、含水層厚度和面波波數(shù)等因素密切相關(guān)[8]：

(1)

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

w(t)——地表位移；

S——儲(chǔ)水系數(shù)；

b——含水層厚度；

kx——地震波波數(shù)。

通過傅立葉變換,可提取出不同周期段信號(hào)的響應(yīng)幅度[14]：

(2)

式中：Zh(ω) ——水位h(t)的傅立葉變換；

Zw(ω) ——地表位移w(t)的傅立葉變換；

ω——頻率；

t——時(shí)間；

圖3所示為2016年12月9日1時(shí)38分(北京時(shí)間)所羅門群島Mw7.8地震引起的水位、垂向速度的時(shí)序圖和頻譜圖,可以看出,相比其它周期的波,20～35 s周期波的響應(yīng)幅度較大。

圖3 面波時(shí)段的水位(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

利用水位和垂向位移的功率譜密度關(guān)系,可以求得井-含水層系統(tǒng)的儲(chǔ)水系數(shù)：

(3)

(4)

(5)

式中：λ——波長(zhǎng),可以由地震波波速和周期確定；

Shh和Sww——水位h(t)和地面運(yùn)動(dòng)垂向位移w(t)的自功率譜密度(PSD)；

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

實(shí)際上,對(duì)于一組有限長(zhǎng)度Tr的記錄數(shù)據(jù),可表示為[15]：

(6)

其中,數(shù)學(xué)期望E表示某個(gè)采樣周期段內(nèi)的平均值。

圖4所示為2016年12月9日1時(shí)38分(北京時(shí)間)所羅門群島Mw7.8地震引起的水位、垂向速度的自功率譜密度(a)和互功率譜密度(b)圖?？梢钥闯?在0.03～0.04 Hz頻段內(nèi)二者的功率譜密度最大,且一致性較高,可用這一顯著頻段的信號(hào)來估算儲(chǔ)水系數(shù)[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,計(jì)算波長(zhǎng)λ=v/ω所需的波速v值由震中距除以面波到時(shí)與發(fā)震時(shí)刻之差獲得,頻率ω取決于圖4中的功率譜密度的最高值。地表運(yùn)動(dòng)(地震波垂向分量)的質(zhì)點(diǎn)位移幅度由速度幅度轉(zhuǎn)化得到,即smax=vmax·τ/2π。依據(jù)式(3)～式(6),利用2016年12月9日1時(shí)38分(北京時(shí)間)所羅門群島Mw7.8地震引起的水位、垂向速度的功率譜密度關(guān)系,估算的新10井含水層平均儲(chǔ)水系數(shù)為1.89×10-5、儲(chǔ)水率為2.36×10-6/m。

2.2 滲透系數(shù)

地震波引起的水震波,其響應(yīng)幅度與井孔條件、含水層參數(shù)、地震波周期等密切相關(guān)。地震波作用除了會(huì)引起地表介質(zhì)的上下波動(dòng)之外,還會(huì)引起含水層內(nèi)孔隙壓力的波動(dòng)變化。井水位的振蕩由地表垂向運(yùn)動(dòng)和孔隙壓力變化共同造成,且井水位對(duì)二者均有一定的放大作用。因此,水震波響應(yīng)幅度既受到井-含水層系統(tǒng)內(nèi)的孔隙壓力波動(dòng)影響,也受到地表垂向運(yùn)動(dòng)的影響。在承壓性的井-含水層系統(tǒng)中,井水位對(duì)孔隙壓力和地表垂向運(yùn)動(dòng)的響應(yīng)可用放大因子(或幅度響應(yīng)比)表示[4]：

(7)

式中：A——水位對(duì)含水層孔隙壓力波動(dòng)的放大因子；

A′——水位對(duì)地表垂向位移的放大因子；

rw——井孔內(nèi)徑；

He——觀測(cè)井有效水柱高度；

T——導(dǎo)水系數(shù)；

Τ——地震波周期；

ω——地震波頻率；

g——重力加速度；

Kei和Ker——零階開爾文函數(shù)。

如果含水層孔隙壓力波動(dòng)引起的水位幅度為h、地表垂向運(yùn)動(dòng)引起的水位幅度為s,那么水位整體波動(dòng)幅度為x=h+s。令R為含水層孔隙壓力波動(dòng)引起的水位變化幅度與地表垂向運(yùn)動(dòng)引起的水位變化幅度之比、M為水位整體波動(dòng)幅度與地表垂向位移幅度之比：

(8)

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

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

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

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

n——含水層孔隙度,可由氣壓系數(shù)、儲(chǔ)水率求得[3]。

依據(jù)式(7)和式(8),利用實(shí)測(cè)的水震波幅度和地震波幅度,可計(jì)算得到相應(yīng)的R和M。

如果令放大因子比A/A′ =R′,那么,幅度比R與地震波周期τ成反比關(guān)系,而放大因子比R’則與周期τ成二次線性正比關(guān)系。R實(shí)測(cè)值計(jì)算所需的波速v值由震中距除以面波到時(shí)與發(fā)震時(shí)刻之差獲得,周期τ值可依據(jù)頻譜圖(圖3)讀?。籖′值計(jì)算需知道有效水柱高度He,可依據(jù)He=H+3d/8計(jì)算得到[4],H為隔水頂板以上的井內(nèi)水柱高度。隔水頂板埋深16 m,水位埋深1 m,井內(nèi)水柱高15 m,那么He=H+3d/8≈18 m。

依據(jù)頻譜圖可逐一統(tǒng)計(jì)得到各周期τ的M值,依據(jù)其R值可計(jì)算得到各自的A及A′,結(jié)果見圖5。圖5中實(shí)測(cè)值由圖3頻譜幅度比計(jì)算所得,理論曲線是利用頻譜幅度比擬合所得,依據(jù)式(7)計(jì)算出與實(shí)際統(tǒng)計(jì)值最接近的滲透系數(shù)K值。計(jì)算中井徑rw=54 mm,含水層厚度d=8 m,有效水柱高度He=18 m,儲(chǔ)水系數(shù)由上述方法獲得S=1.89×10-5。最終,滲透系數(shù)K值擬合結(jié)果為214 m/d。

圖5 新10井水震波放大因子實(shí)測(cè)值與估算值 (a) 水位對(duì)孔隙壓波動(dòng)的放大因子; (b)水位對(duì)地表垂直運(yùn)動(dòng)的放大因子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 結(jié)果與討論

與所羅門群島7.8級(jí)地震引起的同震響應(yīng)類似,新10井水位同樣清晰記錄到了圖1(a)中其它4次地震所引起的同震響應(yīng),它們分別是2016年11月13日新西蘭7.8級(jí)地震 (S 42.757°, E 173.077°)、2016年11月25日中美洲沿岸遠(yuǎn)海7.0級(jí)地震 (N 11.960°, W 88.836°)、2016年12月17日新愛爾蘭地區(qū)7.9級(jí)地震(S 4.509°, E 153.450°) 和2016年12月25日智利7.6級(jí)地震(S 43.416°, W 73.951°)。

基于以上方法,利用其它4次地震引起的新10井水震波和SMG臺(tái)垂直位移波形數(shù)據(jù),估算井-含水層系統(tǒng)參數(shù),見表1。雖然5次地震的震級(jí)、震中距均不相同,尤其是震中距,最小為8 451 km,最大為18 537 km(接近地球周長(zhǎng)的一半),除2016年12月9日估值偏高外,其余參數(shù)基本相當(dāng)。儲(chǔ)水系數(shù)S值為1.33×10-5～1.64×10-5(相對(duì)應(yīng)的儲(chǔ)水率Ss=S/d在1.67×10-6～2.05×10-6/m之間),滲透系數(shù)K值其余值在57～87 m/d之間(相對(duì)應(yīng)的導(dǎo)水系數(shù)T=K·d在456～696 m2/d之間)?；?016年12月9日所羅門群島地震的水震波估算的含水層參數(shù)值偏高的原因,主要是由于在此次地震發(fā)生前的12 h,距新10井約110 km處發(fā)生了呼圖壁6.2級(jí)地震,這種近場(chǎng)地震造成的區(qū)域應(yīng)力調(diào)整,可導(dǎo)致新10井含水層系統(tǒng)介質(zhì)參數(shù)發(fā)生短時(shí)間內(nèi)的明顯變化,如孔隙度、儲(chǔ)水系數(shù)、滲透系數(shù)等增大,進(jìn)而表現(xiàn)為水位的階變上升。隨著區(qū)域應(yīng)力調(diào)整的結(jié)束,各參數(shù)也會(huì)隨即恢復(fù)正常,如2016年12月17日的估值和其它三次地震的估值基本一致。

表1 新10井-含水層系統(tǒng)水文地質(zhì)參數(shù)估算結(jié)果Table 1 Estimated hydraulic parameters of the X10 well-aquifer system

3.1 方法的適用性

地下水位是地震監(jiān)測(cè)預(yù)報(bào)中重要的觀測(cè)項(xiàng)目之一,新10井即屬于這類觀測(cè),如果要獲取這類井的含水層參數(shù),只能通過模型反演實(shí)現(xiàn)。因?yàn)?考慮到這類觀測(cè)井的水位觀測(cè)數(shù)據(jù)必須要連續(xù),抽水試驗(yàn)不能實(shí)現(xiàn)。模型反演獲取水文參數(shù),常用的是水位固體潮調(diào)和分析法[6],即利用水位對(duì)潮汐的響應(yīng)幅度和相位滯后來估算其導(dǎo)水系數(shù)、儲(chǔ)水系數(shù)等。但是,由于井-含水層系統(tǒng)條件的差異,并不是所有觀測(cè)井的水位都能反映出明顯的固體潮汐變化,有些井水位根本就沒有潮汐形態(tài)。

如圖6所示,新10井水位由于氣壓對(duì)其影響較大,導(dǎo)致其水位固體潮汐響應(yīng)看上去很微弱。新10井水位,不論是長(zhǎng)期還是短期,均明顯受到氣壓的影響,水位和氣壓呈負(fù)相關(guān)關(guān)系(相關(guān)系數(shù)R=-0.85),其氣壓系數(shù)[17]為5.2 mm/hPa。但是,水位的固體潮汐形態(tài)卻很微弱,雖然在有些時(shí)候能看到(11月20日—12月1日),但相比氣壓的影響幅度,其變化幅度較小。因此,如果利用新10井水位的固體潮汐變化反演井-含水層水文參數(shù),由于數(shù)據(jù)誤差較大其結(jié)果的可信度會(huì)很低。另外,由于氣壓與固體潮在日波、半日波的某些頻率段其信號(hào)周期是重合的,這會(huì)導(dǎo)致在固體潮調(diào)和分析中必然存在某些誤差。

圖6 2016年11—12月新10井水位與氣壓分鐘值對(duì)比Fig.6 Time series comparisons of the groundwater level (blue) and barometric pressure (green) at the X10 well in November and December, 2016

因此,對(duì)于類似于新10井無清晰固體潮形態(tài)、無抽水試驗(yàn)條件的地震信息觀測(cè)井,利用水位對(duì)面波(瑞利波)的同震響應(yīng)來反演其井-含水層水文參數(shù),是一種不錯(cuò)的選擇。由于遠(yuǎn)場(chǎng)地震引起的水震波只能是地震波作用所引起,其井水位對(duì)地震波的響應(yīng)幅度與井孔條件、地震波特征、含水層參數(shù)等密切相關(guān)。在無其它作用源的情況下,井孔條件一般不會(huì)發(fā)生變化,與地震波特征相關(guān)的水震波形態(tài),其主要影響因素就是含水層參數(shù)。因此,在井孔條件明確的情況下,基于水震波和地震波聯(lián)合反演得到的含水層參數(shù),其結(jié)果的可信度也相對(duì)較高。

3.2 方法的局限性

從時(shí)頻特征分析看,水震波的周期與地震波(瑞利波)基本相當(dāng),多在20 s左右。因此,利用水震波和地震波聯(lián)合反演含水層參數(shù)時(shí),水位的監(jiān)測(cè)時(shí)間間隔應(yīng)小于水震波的周期,本文所用數(shù)據(jù)為秒值。目前,全國(guó)各地震水位觀測(cè)井的數(shù)據(jù)監(jiān)測(cè)時(shí)間間隔多為分鐘值或整點(diǎn)值,無法滿足本文所提出的參數(shù)反演方法需求。因此,推廣地震地下水位的秒值采樣,可更清晰完整地記錄到遠(yuǎn)場(chǎng)地震所引起的水震波,為同震響應(yīng)機(jī)理研究、含水層參數(shù)求取提供可靠的數(shù)據(jù)保障。

地震波作用除了引起地表介質(zhì)的上下波動(dòng)之外,還會(huì)引起含水層內(nèi)孔隙壓力的波動(dòng)變化。井水位的振蕩由地表垂向運(yùn)動(dòng)和孔隙壓力變化共同造成。在對(duì)比分析新10井水震波與地表垂向運(yùn)動(dòng)時(shí),地震波形數(shù)據(jù)選用了離新10井約10 km處的水磨溝臺(tái)垂向速度數(shù)據(jù)。雖然相距10 km波形數(shù)據(jù)的差異會(huì)很小,但畢竟由于地下介質(zhì)的細(xì)微差異,會(huì)導(dǎo)致與實(shí)際值之間存在微小的誤差。因此,理想情況下,利用水震波和地震波反演含水層參數(shù)時(shí),應(yīng)選用與水位觀測(cè)井同一位置處的波形數(shù)據(jù),這樣反演得到的參數(shù)值會(huì)更加準(zhǔn)確。

另外,在用本文探討的方法反演含水層參數(shù)時(shí),應(yīng)盡可能與其它方法(如室內(nèi)實(shí)驗(yàn)[18]、抽水試驗(yàn)[19]、水位固體潮反演等)得到的參數(shù)進(jìn)行對(duì)比分析,這樣得到的結(jié)果可信度將更高。筆者曾用該方法計(jì)算了北京昌平井的含水層參數(shù)[9],計(jì)算結(jié)果與利用水位固體潮反演得到的參數(shù)基本一致,表明利用水震波和地震波垂向速度分量聯(lián)合反演得到的含水層參數(shù)具有較好的可信度。

4 結(jié)論

本文基于新10井水位對(duì)遠(yuǎn)場(chǎng)地震的同震響應(yīng),通過頻譜法分析了該井水震波與地震波的相關(guān)性,提出了一種利用水震波信號(hào)(主要是瑞利波引起)反演獲取承壓含水層水文地質(zhì)參數(shù)的方法。得出以下結(jié)論：

(1)新10井水震波與地震波在變化形態(tài)、頻譜特征上均具有很好的一致性。

(2)通過頻譜關(guān)系分析,估算得到了新10井-含水層系統(tǒng)的水文地質(zhì)參數(shù),其儲(chǔ)水系數(shù)S值在1.33×10-5～1.64×10-5之間,滲透系數(shù)K值在57～87 m/d之間。

(3)基于水震波和地震波信號(hào)估算承壓井-含水層系統(tǒng)水文參數(shù),可為地下水管理、地震前兆異常分析等提供一種新的技術(shù)方法。

本文探討的方法僅應(yīng)用了少數(shù)幾次地震所記錄的數(shù)據(jù),雖然方法本身具有較好的應(yīng)用價(jià)值,但仍需積累更多震例、更多觀測(cè)井?dāng)?shù)據(jù)來加以完善或改進(jìn)。如果能收集到觀測(cè)井自身或附近區(qū)域的抽水試驗(yàn)資料,或者是觀測(cè)井水位具有較清晰的固體潮響應(yīng)形態(tài),就可以利用抽水試驗(yàn)計(jì)算或水位固體潮反演得到的含水層參數(shù)來印證本文探討方法得到的參數(shù)值的可靠性,這也是筆者下一步的研究方向。

致謝：新疆維吾爾自治區(qū)地震局監(jiān)測(cè)中心為本項(xiàng)研究提供了SMG臺(tái)地震波形數(shù)據(jù),在此表示衷心的感謝！

參考文獻(xiàn)：

[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] 葛勤, 梁杏, 龔緒龍,等. 低滲透巖土有效擴(kuò)散系數(shù)的室內(nèi)測(cè)定與分析[J]. 水文地質(zhì)工程地質(zhì), 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] 周志芳, 徐海洋. 一種實(shí)驗(yàn)確定弱透水層水文地質(zhì)參數(shù)的原理與方法[J]. 水文地質(zhì)工程地質(zhì), 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)]