孫小龍,向 陽,2

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

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

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

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

新10井位于中國西部地區(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月底,對其井水位觀測儀器系統(tǒng)進(jìn)行了改造,安置中國地震局地殼應(yīng)力研究所研制生產(chǎn)的SWY-II型數(shù)字水位儀,分辨率為1 mm,監(jiān)測時間間隔為1 s。

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

從圖2(c)和圖2(d)的時頻特征對比圖還可看出,水震波和地震波均存在一個顯著的頻率段,即0.03～0.06 Hz(15～35 s),持續(xù)時間約20 min。面波時段,頻率小、周期大的信號先出現(xiàn),頻率大、周期小的信號后出現(xiàn)。從能量值(紅色表示能量高、藍(lán)色表示能量小)的對比看,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 水文地質(zhì)參數(shù)估算

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

2.1 儲水系數(shù)

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

(1)

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

w(t)——地表位移；

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

b——含水層厚度；

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

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

(2)

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

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

ω——頻率；

t——時間；

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

圖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

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

(3)

(4)

(5)

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

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

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

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

(6)

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

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

2.2 滲透系數(shù)

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

(7)

式中：A——水位對含水層孔隙壓力波動的放大因子；

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

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

He——觀測井有效水柱高度；

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

Τ——地震波周期；

ω——地震波頻率；

g——重力加速度；

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

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

(8)

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

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

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

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

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

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

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

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

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

基于以上方法,利用其它4次地震引起的新10井水震波和SMG臺垂直位移波形數(shù)據(jù),估算井-含水層系統(tǒng)參數(shù),見表1。雖然5次地震的震級、震中距均不相同,尤其是震中距,最小為8 451 km,最大為18 537 km(接近地球周長的一半),除2016年12月9日估值偏高外,其余參數(shù)基本相當(dāng)。儲水系數(shù)S值為1.33×10-5～1.64×10-5(相對應(yīng)的儲水率Ss=S/d在1.67×10-6～2.05×10-6/m之間),滲透系數(shù)K值其余值在57～87 m/d之間(相對應(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級地震,這種近場地震造成的區(qū)域應(yīng)力調(diào)整,可導(dǎo)致新10井含水層系統(tǒng)介質(zhì)參數(shù)發(fā)生短時間內(nèi)的明顯變化,如孔隙度、儲水系數(shù)、滲透系數(shù)等增大,進(jìn)而表現(xiàn)為水位的階變上升。隨著區(qū)域應(yīng)力調(diào)整的結(jié)束,各參數(shù)也會隨即恢復(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)測預(yù)報中重要的觀測項目之一,新10井即屬于這類觀測,如果要獲取這類井的含水層參數(shù),只能通過模型反演實現(xiàn)。因為,考慮到這類觀測井的水位觀測數(shù)據(jù)必須要連續(xù),抽水試驗不能實現(xiàn)。模型反演獲取水文參數(shù),常用的是水位固體潮調(diào)和分析法[6],即利用水位對潮汐的響應(yīng)幅度和相位滯后來估算其導(dǎo)水系數(shù)、儲水系數(shù)等。但是,由于井-含水層系統(tǒng)條件的差異,并不是所有觀測井的水位都能反映出明顯的固體潮汐變化,有些井水位根本就沒有潮汐形態(tài)。

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

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

3.2 方法的局限性

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

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

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

4 結(jié)論

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

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

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

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

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

致謝：新疆維吾爾自治區(qū)地震局監(jiān)測中心為本項研究提供了SMG臺地震波形數(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)測定與分析[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] 周志芳, 徐海洋. 一種實驗確定弱透水層水文地質(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)]