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

江 冰, 唐一格, 夏 強(qiáng), 許 模

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

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

江 冰1,2, 唐一格1, 夏 強(qiáng)1, 許 模1

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

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

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

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

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

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

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

(1)

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

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

(2)

設(shè)第k(k=1,2,3)個(gè)觀測(cè)井的井位為(Xk，Yk)，且對(duì)應(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)θ=θ*(精確解)時(shí)，F(xiàn)(θ)取得極小值(0)，即求目標(biāo)函數(shù)

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

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

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

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

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

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

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

圖1 Matlab程序計(jì)算流程Fig.1 Calculating process of Matlab

2 計(jì)算實(shí)例

2.1 試驗(yàn)場(chǎng)地及試驗(yàn)過(guò)程

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

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

圖2 試驗(yàn)場(chǎng)地平面圖Fig.2 Plane map showing experiment site

圖3 試驗(yàn)場(chǎng)地剖面圖Fig.3 Profile of experiment site

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

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

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

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

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

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

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

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

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

圖5 試驗(yàn)一與試驗(yàn)二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é)果，θ值對(duì)應(yīng)取值區(qū)間為[1°,24°]，Ky/Kx取值區(qū)間為[7.3, 254]。2次試驗(yàn)計(jì)算結(jié)果取其交集為θ值取值區(qū)間為[4°, 24°]，Ky/Kx取值區(qū)間為[7.3, 112]。但考慮抽水試驗(yàn)各觀測(cè)孔的實(shí)際降深情況及泰斯公式自身的局限性，區(qū)間內(nèi)仍存在部分取值范圍不合理。結(jié)合實(shí)際抽水試驗(yàn)過(guò)程，最終確定θ值范圍為[8°, 18°]，Ky/Kx值范圍為[7.5, 13] 。

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

表3 試驗(yàn)二計(jì)算結(jié)果Table 3 Calculating results of Test 2

3 結(jié) 論

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

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

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

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

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

[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] 田開銘,萬(wàn)力.各向異性裂隙介質(zhì)滲透性的研究與評(píng)價(jià)[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ì)算[J].水利學(xué)報(bào),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] 周志芳,莊超,戴云峰,等.單孔振蕩式微水試驗(yàn)確定裂隙巖體各向異性滲透參數(shù)[J].巖石力學(xué)與工程學(xué)報(bào),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] 許模,黃潤(rùn)秋.巖體滲透特性的滲透張量分析在某水電工程中的應(yīng)用[J].成都理工學(xué)院學(xué)報(bào), 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é)與實(shí)踐,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] 周志芳,王錦國(guó).地下水動(dòng)力學(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。

國(guó)家自然科學(xué)基金青年基金項(xiàng)目(41502237)。

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