江 冰, 唐一格, 夏 強, 許 模

(1.地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護國家重點實驗室(成都理工大學(xué)), 成都 610059；2.重慶誠德巖土工程勘察有限公司,重慶 401147)

基于多孔非穩(wěn)定抽水試驗的緩傾層狀含水層各向異性滲透參數(shù)計算

江 冰1,2, 唐一格1, 夏 強1, 許 模1

(1.地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護國家重點實驗室(成都理工大學(xué)), 成都 610059；2.重慶誠德巖土工程勘察有限公司,重慶 401147)

探討通過多孔抽水試驗的降深(s)-時間(t)關(guān)系獲取水文地質(zhì)參數(shù)中的各向異性參數(shù)。利用重慶市渝北區(qū)某場地砂巖含水層的多孔非穩(wěn)定流抽水試驗數(shù)據(jù)，結(jié)合非軸向各向異性滲透系數(shù)張量的理論，通過運用Matlab編程迭代計算了滲透系數(shù)張量的比值及其主要方向。結(jié)果顯示，2組試驗結(jié)果所獲得的數(shù)據(jù)較為吻合，滲透系數(shù)(K)張量的方向θ值范圍為[8°, 18°]，Ky/Kx的取值范圍為[7.5, 13]。該結(jié)果證明了計算方法的適用性和合理性，運用該方法確定含水層各向異性滲透參數(shù)所需的核心參數(shù)為s-lnt直線在零降深上的截距。

抽水試驗；層狀含水層；各向異性；滲透系數(shù)張量

自20世紀60年代以來，通過抽水試驗來確定各向異性巖層滲透張量的研究在國外已有了一定的成果。I.S.Papadopulos[1]最先提出各向異性含水層中的井流公式；隨后M.S.Hantush[2]通過坐標(biāo)變換的方法將各向同性非穩(wěn)定井流公式轉(zhuǎn)變?yōu)楦飨虍愋詶l件下的井流公式。S.Way[3]則提出三維各向異性的滲透系數(shù)的確定方法。S.P.Neuman[4]則提出通過三孔兩次抽水的確定平面滲透張量的方法。在國內(nèi)，田開銘等[5]通過對裂隙巖體滲透性的學(xué)習(xí)并歸納總結(jié)，整合編著出版《各向異性裂隙介質(zhì)滲透性的研究與評價》。周志芳等則于1997年提出通過抽水試驗確定滲透系數(shù)張量的半解析法[6]，并于2015年通過單孔微水振蕩試驗現(xiàn)場確定裂隙導(dǎo)水系數(shù)及巖體滲透系數(shù)張量[7]。相對國外的研究進展，中國學(xué)者對滲透系數(shù)張量的研究多著眼于巖體裂隙特征，許模等[8]通過對巖體裂隙發(fā)育特征的研究確定巖體的各向異性滲透張量；劉青泉等[9]則通過引入裂隙的貫通系數(shù)概念，結(jié)合隨機裂隙網(wǎng)絡(luò)生成技術(shù)，得出估算等效二維裂隙等效滲透系數(shù)張量的疊加算法；而關(guān)于抽水試驗確定巖體各向異性參數(shù)的研究仍然較為缺乏。

本文通過非穩(wěn)定抽水試驗所得降深值，結(jié)合非軸向各向異性滲透系數(shù)張量的理論確定方法，直接計算出目標(biāo)含水層各向異性滲透系數(shù)張量比值及主方向。

1 多孔非穩(wěn)定抽水試驗求參數(shù)的理論與方法

本文研究內(nèi)容為通過定流量非穩(wěn)定流抽水試驗資料即已知地下水水位及流量反求目標(biāo)含水層參數(shù)。周志芳在參考I.S.Papadopulos, M.S.Hantush和S.Way的研究成果后，依據(jù)泰斯公式，通過坐標(biāo)軸的旋轉(zhuǎn)、壓縮變換提出，當(dāng)有3口觀測井在不同時刻的降深資料時，計算各向異性滲透系數(shù)張量主值和主方向的公式[10]

(1)

式中：s表示降深；Q表示流量；t表示時間；δ表示含水層厚度；Kx表示滲透系數(shù)在X軸方向的分量；Ky表示滲透系數(shù)在Y軸方向的分量；θ表示坐標(biāo)系XOY順時針旋轉(zhuǎn)變換的角度。

該公式表明，s與lnt呈線性關(guān)系。對于任一觀測井由直線斜率i，可以得

(2)

設(shè)第k(k=1,2,3)個觀測井的井位為(Xk，Yk)，且對應(yīng)s-lnt直線在零降深上的截距為[tk]，則由公式得

(3)

(4)

(5)

令Pyk=(xksinθ+kkcosθ)2，Pxk=(xkcosθ-yksinθ)2，帶入上式中可整理為

[t1]2(Py2Px3-Py3Px2)+[t1][t2](Py3Px1-Py1Px3)+[t1][t3](Py1Px2-Py2Px1)=0

(6)

設(shè)

f(θ)=[t1]2(Py2Px3-Py3Px2)+[t1][t2](Py3Px1-Py1Px3)+[t1][t3](Py1Px2-Py2Px1)

(7)

設(shè)

F(θ)=[f(θ)]2≥0

當(dāng)θ=θ*(精確解)時，F(xiàn)(θ)取得極小值(0)，即求目標(biāo)函數(shù)

E(θ)=minF(θ) (A0≤θ≤B0)

的最優(yōu)解。A0和B0為θ角變化的范圍，在此設(shè)為0°和90°。根據(jù)求得的θ角，可帶入公式(3)、(4)、(5)，從而計算出Ky與Kx的比值。

利用抽水試驗數(shù)據(jù)帶入理論公式具體步驟如下：

a.利用降深s與時間t半對數(shù)曲線獲取各觀測孔所對應(yīng)的截距[tk](k=1,2,3)。

b.將[tk](k=1,2,3)和各觀測孔坐標(biāo)(Xk，Yk)代入式(7)和(8)，得到滲透系數(shù)張量的主要方向θ值。

c.將θ值帶入式(3)、(4)、(5)聯(lián)立對比最終得到Ky/Kx。

通過使用Matlab 2014將上述理論公式程序化，可直接輸入各監(jiān)測孔坐標(biāo)和對應(yīng)截距而得出滲透系數(shù)張量方向θ值和Ky/Kx，這樣簡化了計算工作量，并且提高了計算效率和準確性，程序計算流程見圖1。

圖1 Matlab程序計算流程Fig.1 Calculating process of Matlab

2 計算實例

2.1 試驗場地及試驗過程

試驗場地位于重慶市渝北區(qū)，共布設(shè)4口鉆孔(圖2)，以抽水孔ZK1為原點，確定3個觀測孔坐標(biāo)(表1)。試驗場地覆蓋層厚1.5 m，下伏地層巖性為砂泥巖互層的單斜地層，傾角lt;10°，為緩傾角巖層，與上述公式的含水層產(chǎn)狀近似水平的基本條件一致。第一層泥巖厚約15 m，目標(biāo)含水層砂巖層厚約25 m(圖3)。

在此共進行9次抽水試驗，水位降深值由儀器Diver自動記錄，文中選取其中2次試驗數(shù)據(jù)進行計算，2次試驗s-t曲線見圖4。

圖2 試驗場地平面圖Fig.2 Plane map showing experiment site

圖3 試驗場地剖面圖Fig.3 Profile of experiment site

試驗一于2016年8月15日9時30分開始，8月16日11時00分停止抽水，持續(xù)時長25.5 h，平均抽水速率3.17 m3/h。抽水孔ZK1內(nèi)最大降深為2.68 m；各觀測孔ZK2、CK1、CK2內(nèi)的最大降深分別為2.42 m、2.43 m、2.55 m。

由s-t曲線可以看出，試驗一進行過程中，曾在500～900 min時間段內(nèi)水位發(fā)生震蕩，根據(jù)Q-t曲線基本可以判斷是由于水泵功率不穩(wěn)定，在這段試驗過程中，流量起伏過大而造成，但對整體試驗結(jié)果無實質(zhì)性影響。

圖4 試驗一與試驗二s-t及Q-t曲線Fig.4 Q-t and s-t curves of Test 1 and Test 2

試驗二于 2016年8月19日12時10分開始， 8月20日6時10分停止抽水，持續(xù)時長18 h。累計抽水量273.43 m3，平均抽水速率3.15 m3/h。抽水孔ZK1內(nèi)最大降深為4.1 m，各觀測孔ZK2、CK1、CK2內(nèi)的最大降深分別為3.51 m、3.54 m、3.51 m。

2.2 各向異性滲透參數(shù)求解

本文計算各向異性滲透參數(shù)的核心在于作出降深值(s)與時間(t)的半對數(shù)曲線，以獲取能夠直接帶入Matlab程序中進行迭代計算的s-lnt直線在零降深上的截距為[tk](k=1,2,3)，在此兩次試驗統(tǒng)一設(shè)監(jiān)測孔CK1截距為[t1]，ZK2為[t2],CK2為[t3]。

分別作出試驗一與試驗二的s-lnt曲線(圖5)。從圖中可以看出，s-lnt的完整曲線是隨時間的增加逐漸趨近為一條直線，而并非絕對的直線。從圖中所作虛線可以看出，s-lnt曲線的斜率實際上需要確定一個有效范圍，而這樣會導(dǎo)致其截距也產(chǎn)生了相應(yīng)的取值區(qū)間，并且2個s-lnt曲線都能看出CK1與ZK2曲線基本重合，所以在確定各觀測孔對應(yīng)截距取值區(qū)間時，可取[t1]= [t2]。

試驗一中，根據(jù)曲線確定[t1]和[t2]取值區(qū)間為[55, 60]，[t3]取值區(qū)間為[45, 50]。為提高計算效率，截距取值均取整數(shù)。將所有取值組合輸入程序，最終得到以下計算結(jié)果(表2)。

試驗二中，確定[t1]和[t2]取值區(qū)間為[30,35]，[t3]取值區(qū)間為[25, 30]。將數(shù)據(jù)帶入程序中運行得到結(jié)果見表3。

圖5 試驗一與試驗二s-lnt曲線Fig.5 s-lnt curves of Test 1 and Test 2

[t3][t1]454647484950θKy/KxθKy/KxθKy/KxθKy/KxθKy/KxθKy/Kx5518.013.215.010.012.08.59.07.76.07.44.07.25620.519.317.017.514.59.611.08.39.07.76.07.45723.043.020.017.217.011.614.09.311.08.28.67.65819.015.616.511.013.69.110.88.05821.726.219.014.416.010.613.08.96024.0112.021.022.018.013.415.010.0

2.3 數(shù)據(jù)討論

據(jù)表2結(jié)果，各向異性滲透張量的主要方向θ取值范圍為[4°, 24°]，Ky/Kx取值范圍為[7.2, 112]。而據(jù)表3結(jié)果，θ值對應(yīng)取值區(qū)間為[1°,24°]，Ky/Kx取值區(qū)間為[7.3, 254]。2次試驗計算結(jié)果取其交集為θ值取值區(qū)間為[4°, 24°]，Ky/Kx取值區(qū)間為[7.3, 112]。但考慮抽水試驗各觀測孔的實際降深情況及泰斯公式自身的局限性，區(qū)間內(nèi)仍存在部分取值范圍不合理。結(jié)合實際抽水試驗過程，最終確定θ值范圍為[8°, 18°]，Ky/Kx值范圍為[7.5, 13] 。

四川盆地紅層分布區(qū)構(gòu)造作用較弱的緩傾斜砂巖巖體，其滲透性各向異性三維特征主要受控于巖層產(chǎn)狀，呈現(xiàn)出以巖層傾向、走向和法向為三向滲透優(yōu)勢主方向，統(tǒng)計上顯示以傾向與走向為滲透性最強的2個主要優(yōu)勢方向，法向(一般多為垂向)一般為最小滲透優(yōu)勢方向。本文計算結(jié)果推測出該區(qū)域巖體滲透性存在較明顯的方向性，其優(yōu)勢方向為近南北向，這與研究區(qū)砂巖地層走向一致，顯示的是巖層走向的主控作用。

表3 試驗二計算結(jié)果Table 3 Calculating results of Test 2

3 結(jié) 論

a.運用Matlab將公式程序化是使用黃金分割法迭代計算參數(shù)最高效的方式之一。

b.通過多孔非穩(wěn)定抽水試驗，運用試驗參數(shù)成功反求各向異性滲透參數(shù)，將2次試驗數(shù)據(jù)的反求結(jié)果對比分析得到以下結(jié)果：滲透系數(shù)張量的方向θ值范圍為[8°, 18°]，Ky/Kx的取值范圍為[7.5, 13]。

c.在計算過程中，參數(shù)的取值尤為關(guān)鍵，由于s-lnt曲線并不完全是一條直線，而是越來越趨近于直線，這樣會導(dǎo)致曲線在零降深上的截距取值為一個區(qū)間而非唯一值。

d.s-lnt曲線在零降深上的截距是一個敏感度極高的參數(shù)，參數(shù)值的微小變化甚至可以導(dǎo)致Ky/Kx發(fā)生大于10倍的變化，因此確定其有效取值范圍是極為重要的。

e.本文研究的目標(biāo)含水層實際為潛水，雖不完全滿足泰斯公式的承壓含水層條件，但計算結(jié)果與實際各鉆孔水位降深情況基本吻合，這樣的計算理論是適用的。

[1] Papadopulos I S. Nonsteady flow to a well in an infinite anisotropy aquifer[J]. Symposium of Dubrovnik, 1965, 1: 21-31.

[2] Huntash M S. A method for analyzing a drawdown test in anisotropic aquifer[J]. Water Resources Research, 1966, 2(2): 281-285.

[3] Way S C. In-situ determination of three-dimensional aquifer permeabilities[J]. Ground Water, 1982, 20(5): 595-603.

[4] Neuman S P, Walter G R, Bentley H W,etal. Determination of horizontal aquifer anisotropy with three wells[J]. Ground Water, 1984, 22(1): 66-72.

[5] 田開銘,萬力.各向異性裂隙介質(zhì)滲透性的研究與評價[M].北京:學(xué)苑出版社, 1989.

Tian K M, Wan L. Research and Evaluating of Anisotropic Fissured Medium Permeability[M]. Beijing: Xueyuan Publishing House, 1989. (in Chinese)

[6] 周志芳,朱學(xué)愚,李艷.巖體滲透系數(shù)張量的半解析計算[J].水利學(xué)報,1997(9):6-11.

Zhou Z F, Zhu X Y, Li Y. Semi-analytic method for conductivity tensor in rock mass[J]. Journal of Hydraulic Engineering, 1997(9): 6-11. (in Chinese)

[7] 周志芳,莊超,戴云峰,等.單孔振蕩式微水試驗確定裂隙巖體各向異性滲透參數(shù)[J].巖石力學(xué)與工程學(xué)報,2015,34(2):271-278.

Zhou Z F, Zhuang C, Dai Y F,etal. Determining anisotropic hydraulic conductivity in fractured rocks based on single-borehore slug tests[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(2): 271-278. (in Chinese)

[8] 許模,黃潤秋.巖體滲透特性的滲透張量分析在某水電工程中的應(yīng)用[J].成都理工學(xué)院學(xué)報, 1997,24(1):56-63.

Xu M, Huang R Q. Application of permeability tensor analysis to study rock mass permearbility of a certain hydroelectric project in southwest of China[J]. Journal of Chengdu University of Technology, 1997, 24(1): 56-63. (in Chinese)

[9] 劉青泉,樊紅光.二維裂隙滲流等效滲透系數(shù)張量的一種近似算法[J].力學(xué)與實踐,2012,34(5):16-20.

Liu Q Q, Fan H G. An approximate method for calculating the equivalent permeability tensor of seepage in complicated fractures[J]. Mechanics in Engrneering, 2012, 34(5):16-20. (in Chinese)

[10] 周志芳,王錦國.地下水動力學(xué)[M].北京：科學(xué)出版社,2013:135-137.

Zhou Z F, Wang J G. Groundwater Dynamics[M]. Beijing: Science Press, 2013: 135-137. (in Chinese)

Calculationofanisotropicpermeabilityparametersofflatpitchaquiferbasedonmultiple-wellunsteadypumpingtest

JIANG Bing1,2, TANG Yige1, XIA Qiang1, XU Mo1

1.StateKeyLaboratoryofGeohazardPreventionandGeoenvironmentProtection,ChengduUniversityofTechnology,Chengdu610059,China;2.ChongqingChengdeGeotechnicalEngineeringInvestigationCo.Ltd,Chongqing401147,China

The parameters of anisotropy of hydrogeology can be obtained through drawdown-time relationship of multiple-well pumping test. In this paper, the multiple-well unsteady pumping test data are used to calculate the main direction and specific value of permeability tensor by Matlab incorporating with non-axial anisotropic permeability tensor theory. The final result which is consistent with the actual drawdown condition is obtained after contrastive analysis from two separate test results, with the main direction range of [8°, 18°], and theKy/Kxrange of [7.5, 13]. The calculation method proves to be feasible and rational for the non-axial anisotropic permeability tensor theory. It points out that the key parameter to calculate anisotropic permeability parameters is the intercept ofs-lntcurve.

pumping test; layered aquifer; anisotropy; permeability tensor

P641.2

A

10.3969/j.issn.1671-9727.2017.06.13

1671-9727(2017)06-0756-06

2017-01-16。

國家自然科學(xué)基金青年基金項目(41502237)。

江冰(1963-)，男，高級工程師，研究方向：水文地質(zhì)與工程地質(zhì), E-mail:448124138@qq.com。