王玉林,謝康和,黃大中,李傳勛

(1-武夷學(xué)院土木工程與建筑學(xué)院，福建武夷山 354300;2-浙江大學(xué)軟弱土與環(huán)境土工教育部重點實驗室，杭州 310027;3-鐵道第三勘察設(shè)計院集團(tuán)有限公司，天津 300251;4-江蘇大學(xué)土木工程與力學(xué)學(xué)院，鎮(zhèn)江 212013)

1 引言

石油和地下水開采井在鉆進(jìn)和形成過程中井周土受到擠壓、涂抹和污染作用使井周擾動土的滲透性低于未擾動土的滲透性，而形成“正表皮層”(positive skin)；另一方面，在抽取地下流體的生產(chǎn)過程中，由于井周涂抹層剝落、開裂或土層中的細(xì)微顆粒被帶走使得井周土的滲透系數(shù)大于原有結(jié)構(gòu)層的滲透性，而形成“負(fù)表皮層”(negative skin)[1-4]．井周表皮層厚度往往從幾毫米到幾米不等，在井流理論研究中，根據(jù)表皮層厚度或者是否考慮表皮層的彈性釋水，建立了兩種具有很大差異的兩種數(shù)學(xué)模型進(jìn)行滲流問題求解[1-3]．當(dāng)表皮層厚度較薄或者不需要考慮表皮層彈性儲水時，可引入表皮效應(yīng)系數(shù)，建立徑向單層井流模型進(jìn)行地下滲流模擬[4,5]，而當(dāng)表皮層厚度較大，需要考慮表皮層彈性儲水時，應(yīng)采用徑向雙層的井流模型才能準(zhǔn)確的分析地下流體的滲流作用[1-3,6,7]．Yeh和Yang等學(xué)者對考慮表皮層彈性儲水的徑向雙層井流模型理論進(jìn)行較為深入的研究[1,6-8]．本文對已有的徑向雙層、各向同性承壓完整井二維井流模型作了改進(jìn)，考慮了表皮層各向異性及其徑向滲透系數(shù)的豎向非均質(zhì)性，考慮了未擾動承壓含水層的各向異性以及豎向越流補(bǔ)給作用，建立相應(yīng)的三維井流數(shù)學(xué)模型，通過Laplace變換和矩陣?yán)碚撉蟮帽砥雍臀磾_動承壓層內(nèi)水頭降分布和井壁流量的半解析解，應(yīng)用實例分析了表皮層和未擾動承壓層的滲透性對水頭降和井流量的影響．

2 數(shù)學(xué)模型

2.1 基本假設(shè)與滲流控制方程

對于表皮層為有限厚度，并且考慮其彈性儲水的承壓含水層完整井井流問題，Yeh和Yang等[1,6,7,9]引入了如下基本假設(shè)：

1)表皮層和未擾動承壓層為均質(zhì)、各向同性，含水層側(cè)向無限延伸，表皮層等厚分布于井周；

2)抽水井為有限直徑的完整井；

3)承壓含水層初始水頭為常數(shù)；

4)井周表皮層為有限厚度．

基于這些假設(shè)，他們建立了徑向雙層承壓含水層的二維井流模型，其中表皮層和未擾動承壓層的滲流控制方程為

其中sI,SsI和KIr分別為表皮層的水頭降、比彈性釋水系數(shù)(或貯水率)和徑向滲透系數(shù)；sII,SsII和KIIr分別為未擾動承壓層的水頭降、比彈性釋水系數(shù)(或貯水率)和徑向滲透系數(shù)．

實際上，表皮層和未擾動承壓含水層多為非均質(zhì)、各向異性，并且承壓地下水存在垂直方向的越流補(bǔ)給作用．據(jù)此，本文對Yeh和Yang的上述假設(shè)基礎(chǔ)上作了改進(jìn)，不僅考慮了表皮層各向異性及其徑向滲透系數(shù)的沿井深方向(豎向)的非均質(zhì)性，而且考慮了未擾動承壓含水層的各向異性，以及弱透水層的越流補(bǔ)給作用，建立徑向雙層、各向異性承壓含水層的三維井流數(shù)學(xué)模型，其控制方程為

其中KIr(z)為非均質(zhì)性表皮層的徑向滲透系數(shù)，它是關(guān)于z的函數(shù)；KIv和KIIv分別為表皮層和未擾動承壓含水層的豎向滲透系數(shù)．為便于問題求解，根據(jù)表皮層徑向滲透性沿井深方向的非均質(zhì)性，將表皮層和未擾動承壓層沿豎直方向(z方向)統(tǒng)一劃分N層，并自上而下依次編號，這樣可將原問題轉(zhuǎn)化為徑向為雙層、豎向為N層的承壓含水層系統(tǒng)井流概念模型，如圖1所示．

圖1: 考慮非均質(zhì)各向異性表皮層的徑向雙層承壓層的井流概念模型

2.2 徑向雙層各向異性承壓含水層非穩(wěn)態(tài)井流的數(shù)學(xué)模型

根據(jù)上述假設(shè)與圖1所示概念模型，可建立豎向N層越流系統(tǒng)中的表皮層和未擾動承壓層的地下水三維滲流方程，其中N層表皮層(rw≤r≤rs)內(nèi)的地下水三維滲流方程表達(dá)為

越流系統(tǒng)中N層未擾動承壓層(rs≤r＜∞)內(nèi)地下水三維滲流方程表達(dá)為

上式中sIi和sIIi分別為第i層表皮層和第i層未擾動承壓層的水頭降，假設(shè)上覆弱透水層無水頭降sI0=sII0=0，底部隔水層有sIN=sIIN+1=0；TIi=KIrihi和TIIi=KIIrihi分別為第i層表皮層和未擾動承壓層的水頭傳導(dǎo)系數(shù)；μeIi和μeIIi分別為第i層表皮層和未擾動承壓層的彈性釋水系數(shù)；

分別為表皮層和未擾動承壓層在第i層與第i?1層間之間的阻力系數(shù)，特別地，CI1=，CII1=和CIN→∞；KIri和KIvi分別為第i層表皮層徑向和豎向滲透系數(shù)，KIIri和KIIvi分別為第i層未擾動承壓層徑向和豎向滲透系數(shù)；hi為第i層土的厚度，h0為弱透水層厚度，i=1,2,3,···.

上述兩方程組系統(tǒng)的求解條件：

(a)抽水井定降升邊界條件：sIi(rw,t)=ssw；

(b)無限遠(yuǎn)處的自然邊界條件：sIIi(∞,t)=0；

(c)表皮層和未擾動承壓層之間的水頭連續(xù)條件：sIi(rs,t)=sIIi(rs,t),t＞0；

(d)表皮層和未擾動承壓層之間的流量連續(xù)條件：

(e)初始條件：sIi(r,0)=sIIi(r,0),r＞rw；

其中sw為豎井內(nèi)的降深，rw為抽水井半徑，rs為表皮層外徑．

3 問題求解

分別對方程組(5)和方程組(6)進(jìn)行Laplace變換可得到線性常微分方程組，并用微分算子進(jìn)一步簡化成矩陣形式的常微分方程組，形如

其中?2=為Laplace算子，I=[I1I2···IN]T為Laplace域內(nèi)各層表皮層水頭降組成的向量，II=[II1II2···IIN]T為Laplace域內(nèi)各層未擾動承壓層水頭降組成的向量，A和B分別為表皮層和未擾動承壓層的滲流控制方程對應(yīng)的矩陣，均為N×N的方陣，其中A矩陣可寫為

矩陣B的形式與矩陣A的形式一樣，將矩陣A各元素中的下標(biāo)I改為II即可得矩陣B．根據(jù)矩陣?yán)碚摚菀浊蟮镁仃嘇的N個特征值為：λI1,λI2,···,λIN，各特征值對應(yīng)向量組成的向量矩陣為YI，B的N個特征值為：λII1,λII2,···,λIIN，對應(yīng)向量組成的向量矩陣為YII．于是可利用A=YIλIYI?1和B=YIIλIIYII?1，將線性矩陣微分方程組(7)變?yōu)榻怦钚问降姆匠探M

將方程組(9)結(jié)合表皮層和未擾動承壓層中的滲流Laplace變換域求解條件，可得到表皮層和未擾動承壓層水頭降的Laplace域內(nèi)一般解

其中F,G和H為由邊界條件確定的待定常數(shù)向量，因為未擾動承壓層和表皮層沿著z方向劃分N層，因此每個常數(shù)向量均含有N個分量；為以零階修正Bessel函數(shù)為對角元素的對角方陣．

可由井周邊界條件(a)得到

由表皮層和未擾動承壓層邊界上的水頭連續(xù)條件(c)得到

再由表皮層和未擾動承壓層邊界上的流量連續(xù)條件(d)得到

為了便于計算機(jī)的編程計算，將式(11)—(13)合并成矩陣的形式

可用數(shù)學(xué)軟件程序求線性方程組(14)的解，即求得常數(shù)向量F,G和H：

將求解得到的常數(shù)向量F,G和H代入式(10)，即可分別得到表皮層和未擾動承壓層的水頭降 Laplace域上的解I和II．

同樣地，可對I關(guān)于變量r求導(dǎo)得到井周邊界處r=rw各分層的流量

由于上述解較復(fù)雜，直接用Laplace逆變換求得物理域解析解表達(dá)式比較困難，因此，可采用Stefest數(shù)值反演方法求得時間域上的物理量[10]．

4 解的討論

1)當(dāng)未擾動承壓層和表皮層均為均勻、各向同性時，即CIi,CIIi,μeIi,μeIIi,TIi以及TIIi各層相等，取為常數(shù)，則本文徑向雙層、豎向N層的承壓含水層系統(tǒng)抽水井流問題退化為徑向雙層、豎向單層的承壓含水層越流完整抽水井流問題[11]；進(jìn)一步，當(dāng)承壓含水層上層為隔水邊界，無越流補(bǔ)給時，即CIi→∞，則本文解可為考慮表皮層彈性釋水的定降深承壓井流問題的解[9]．

2)當(dāng)考慮承壓含水層的層狀特性，但不考慮表皮效應(yīng)時，即CI,CII和TI=TII，本文解退化為越流多層承壓含水層系統(tǒng)井流解[12]．

3)當(dāng)抽水時間較長時，即p→0，承壓含水層中的地下水達(dá)到穩(wěn)態(tài)流狀態(tài)，此時對應(yīng)的穩(wěn)態(tài)流水頭降解為

當(dāng)承壓含水層達(dá)到穩(wěn)態(tài)流狀態(tài)時，抽水井周各分層對應(yīng)的流量分量為

其中和分別為矩陣A和矩陣B中p→0時所求得的特征值．

5 表皮層非均質(zhì)的徑向雙層各向異性承壓層井流分析

應(yīng)用數(shù)學(xué)軟件編制計算程序，通過算例對徑向雙層各向異性承壓含水層水頭降和流量及影響因素進(jìn)行分析．由算例得到的分析曲線圖，坐標(biāo)系物理量采用無量綱形式：水位降、豎井的抽水流量、徑向坐標(biāo)、豎向坐標(biāo)．

在算例中，假設(shè)表皮層為層狀非均質(zhì)、各向異性的，徑向滲透系數(shù)KIr沿豎井方向(z方向)層狀隨機(jī)變化，可用分段連續(xù)函數(shù)表示為：KIr(z)=βj(z)(βj(z)為第j段函數(shù))，如圖2所示，分段函數(shù)共7個分段(層數(shù))，自上而下各層徑向滲透系數(shù)的取值為

在計算過程中，將每一分段函數(shù)KIr(z)對應(yīng)層位的層表皮層和非擾動承壓層(zj?1≤z≤zj)再劃分6層，因此，共劃分成42層，取計算層中點位置處的數(shù)據(jù)繪制各相應(yīng)圖形．

圖2: 表皮層徑向滲透系數(shù)KIr豎向?qū)訝铍S機(jī)變化示意圖

5.1 承壓含水層水頭降沿豎向分布規(guī)律及影響因素分析

根據(jù)圖2所示表皮層徑向滲透系數(shù)KIr沿z方向的層狀隨機(jī)變化，繪制得到不同井距處承壓含水層水頭降沿豎向變化的曲線s??z?，如圖3所示，圖中各計算點的位置分別為：r?=10和r?=100在表皮層內(nèi)，r?=300在表皮層和未擾動承壓含水層交界面處，r?=1000和r?=1500在未擾動承壓含水層內(nèi)．

對比圖3中不同r?位置處的s??z?曲線可知：在表皮層內(nèi)以及在表皮層附近處，水頭降沿豎向的變化趨勢與圖2中表皮層徑向滲透系數(shù)KIr沿井深變化趨勢具有較好的一致性，即，各層承壓含水層的水頭降隨同層位表皮層徑向滲透系數(shù)值βj的大小變化而變化，水頭降沿豎向分布曲線隨層位變化較明顯，能夠從水頭降變化直觀的反映出不同層位表皮層滲透性．但是，隨著與表皮層距離增大，表皮效應(yīng)逐漸減弱，承壓含水層的水頭降s??z?曲線隨表皮層分層變化的特征逐漸消失．

圖3: 表皮層徑向滲透系數(shù)KIr沿井深分層變化時不同井距處的水頭降沿豎向分布s??z?曲線

為研究表皮層對不同滲透性的各向異性承壓含水層的影響，計算了四種滲透性的承壓含水層的水頭降，繪制得到在距井較近處(r?=400)和距井較遠(yuǎn)處(r?=1500)的s??z?曲線，分別如圖4和圖5所示．

對比圖4和圖5可以知道：表皮層對承壓含水層的影響范圍與承壓含水層的豎向滲透系數(shù)KIIv密切相關(guān)，即，承壓含水層的豎向滲透系數(shù)KIIv越小(大)，表皮層對承壓含水層滲流產(chǎn)生影響的范圍就越大(小)，而且表皮效應(yīng)越明顯(不明顯)．

圖4: KIr沿在較小井距r?=400處水頭降沿豎向分布s??z?曲線

圖5: KIr沿在較大井距r?=1500處水頭降沿豎向分布s??z?曲線

具體分析來看，承壓含水層豎向滲透系數(shù)較小時(KIIv=10?6m/s)，在距井較近處(r?=400)和較遠(yuǎn)處(r?=1500)水頭降沿豎向的分布形態(tài)均與圖2中的表皮層徑向滲透系數(shù)KIr沿井深層狀變化有較好的一致性，水頭降較明顯地隨表皮層的層位不同而變化，表明：KIIv較小時，在較遠(yuǎn)處承壓含水層地下水仍受到較強(qiáng)的表皮效應(yīng)作用；然而，當(dāng)豎向滲透系數(shù)較大時(KIIv=10?4m/s)，僅在距井較近處(r?=400)水頭降沿豎向的分布形態(tài)才與表皮層徑向滲透系數(shù)KIr的層狀變化有較好的一致性，而在距井較遠(yuǎn)處(r?=1500)水頭降沿豎向的分布形態(tài)不再具有與表皮層對應(yīng)的層狀性，表明：KIIv較大時，在較遠(yuǎn)處承壓含水層地下水受表皮效應(yīng)的作用不明顯．

產(chǎn)生上述現(xiàn)象的原因在于：承壓含水層豎向滲透系數(shù)KIIv較小時，上下相鄰層之間的越流較弱，各層中的地下水滲流受豎向水流干擾較小，而受同層位表皮層的影響較強(qiáng)，因此水頭降沿豎向的分布形態(tài)能較明顯地與表皮層的層狀性對應(yīng)；但是，當(dāng)承壓含水層豎向滲透系數(shù)KIIv較大時，上下相鄰層之間的越流較強(qiáng)，各層中的地下水滲流受豎向水流干擾較大，而受同層位表皮層的影響較弱，因而水頭降沿豎向的分布形態(tài)不再與表皮層的層狀性對應(yīng)．

5.2 井壁流量分布規(guī)律及影響因素分析

圖6為不同表皮層厚度對井壁處流量Q??z?分布曲線的影響．該圖表明：

1)表皮層的層狀性導(dǎo)致水流量沿抽水井井壁非均勻分布，流量沿井深變化趨勢與表皮層徑向滲透系數(shù)變化趨勢一致，表皮層徑向滲透系數(shù)值βj越大的層位，相應(yīng)位置處井壁的水流量就越大，反之則??；

2)對表皮層厚度的影響而言，有正、負(fù)表皮效應(yīng)之分，在“正表皮層”層位處，表皮層的厚度越大，則相應(yīng)層位處井壁的水流量越小，而“負(fù)表皮層”層位處，表皮層的厚度越大，則相應(yīng)層位處井壁的水流量也越大．

由此可以看出，在定降深條件下，可通過增大表皮層滲透系數(shù)或增大“負(fù)表皮層”厚度提高抽水井產(chǎn)量，因此改善井周巖土介質(zhì)滲透性或設(shè)置滲透性大于承壓含水層滲透性的圍填材料，對產(chǎn)量增加具有明顯的效果．

圖6: 不同厚度表皮層對井壁流量分布Q??z?影響的比較

圖7 為四種各向異性承壓含水層的井壁水流量Q??z?分布曲線．對比表皮層不同層位處對應(yīng)的井壁流量(產(chǎn)量)可看出：表皮層徑向滲透系數(shù)越大，不同滲透性的承壓含水層之間井流量(或產(chǎn)量)的差別越大，承壓含水層的滲透性對流量(或產(chǎn)量)影響就越大，因此承壓含水層的滲透性(包括徑向滲透系數(shù)KIIr和豎向滲透系數(shù)KIIv)對抽水井流量(或產(chǎn)量)的影響與表皮層徑向滲透系數(shù)有關(guān)．就具體分析來看，第1、4和5層表皮層的徑向滲透系數(shù)值較小，這些層位處各種承壓含水層之間的流量差異較小，相應(yīng)的Q??z?曲線接近重合，而其他層位處表皮層的徑向滲透系數(shù)值較大，各種承壓含水層之間的流量差異也較大，相應(yīng)的Q??z?曲線明顯分離．從總體上來看，七層表皮層的徑向滲透系數(shù)值從小到大排列為β1,β2,β3,β4,β5,β6,β7，各層位處四種不同承壓含水層的Q??z?曲線的分離程度隨βj值增大而增大，說明承壓含水層的滲透性對井壁流量Q?(產(chǎn)量)的影響隨表皮層的徑向滲透系數(shù)βj值增大而增大．

圖7: 不同類型各向異性承壓層的井壁流量分布s??z?影響的比較

5.3 承壓含水層和表皮層的滲透性對水頭降徑向分布s??r?的協(xié)同影響

根據(jù)圖2所示表皮層徑向滲透系數(shù)KIr沿井深的分布模式，對五種不同承壓含水層的水頭降沿徑向分布進(jìn)行計算，取其中表皮層徑向滲透系數(shù)KIr為較小值的層位(即頂層)和較大值的層位(即底層)的水頭降曲線進(jìn)行比較分析，如圖8和圖9所示．對比圖8和圖9水頭降曲線可知：

1)不管表皮層的徑向滲透系數(shù)值較小(β1)還是較大(β7)，水頭降幅值均隨著承壓含水層滲透系數(shù)(包括徑向滲透系數(shù)KIIr和豎向滲透系數(shù)KIIv)減小而增大；

2)表皮層徑向滲透系數(shù)KIr較小值時(如：β1=10?5m/s)，承壓含水層的水頭降幅度主要取決于承壓含水層徑向滲透系數(shù)KIIr的大小，承壓含水層豎向滲透系數(shù)KIIv對承壓含水層水頭降影響較?。@容易從圖8看出，對于承壓含水層豎向滲透系數(shù)KIIv相同，徑向滲透系數(shù)KIIr不同的任意兩條曲線，水頭降差異均較大，但是對于承壓含水層徑向滲透系數(shù)KIIr相同，豎向滲透系數(shù)KIIv不同的任意兩條曲線，水頭降差異均較小；

3)表皮層徑向滲透系數(shù)KIr的值為較大時(如：β7=50×10?5m/s)，承壓含水層徑向滲透系數(shù)KIIr和豎向滲透系數(shù)KIIv均對承壓含水層的水頭降有較大影響．這容易從圖9看出，對于承壓含水層豎向滲透系數(shù)KIIv相同，而徑向滲透系數(shù)KIIr不同的任意兩條曲線，水頭降差異較大，若對比承壓含水層徑向滲透系數(shù)KIIr相同，而豎向滲透系數(shù)KIIv不同的兩條曲線，水頭降差異均較大；

4)同一承壓含水層水頭降曲線s??r?在交界面處是呈向上折趨勢，還是呈下折趨勢，由所在層位的表皮層的負(fù)、正表皮效應(yīng)性質(zhì)決定．

圖8: 表皮層KIr取值較小為β1=10?5m/s時五種承壓含水層s??r?曲線比較

圖9: 表皮層KIr取值較大為β7=50×10?5m/s時五種承壓含水層s??r?曲線比較

6 結(jié)論

本文改進(jìn)了Yeh和Yang的徑向雙層、各向同性承壓完整井二維井流模型，考慮表皮層和未擾動承壓層的非均質(zhì)和各向異性以及層間的豎向越流補(bǔ)給作用，通過將表皮層和未擾動承壓層沿豎向劃分為有限層，建立了徑向雙層豎向有限層的三維井流數(shù)學(xué)模型，采用Laplace變換和矩陣?yán)碚撉蟮贸袎汉畬铀^降和井壁流量的半解析解，通過算例分析得到如下結(jié)論：

1)在表皮層內(nèi)以及在表皮層附近處，水頭降沿豎向的變化趨勢與表皮層徑向滲透系數(shù)KIr變化趨勢具有較好的一致性，而在遠(yuǎn)處表皮效應(yīng)逐漸減弱，承壓層水頭降隨表皮層分層變化的特征逐漸消失；

2)承壓含水層的豎向滲透系數(shù)KIIv越小(大)，表皮層對承壓含水層滲流產(chǎn)生影響的范圍就越大(小)．原因在于：承壓含水層豎向滲透系數(shù)KIIv大小決定豎向各分層之間越流強(qiáng)弱，進(jìn)而影響表皮效應(yīng)的強(qiáng)弱；

3)流量沿井深變化趨勢與表皮層徑向滲透系數(shù)變化趨勢一致，通過改善井周巖土介質(zhì)滲透性或增大“負(fù)表皮層”厚度可提高抽水井產(chǎn)量；

4)表皮層徑向滲透系數(shù)KIr越大，不同滲透性承壓含水層之間的井流量(或產(chǎn)量)差別越大，承壓含水層的滲透性對流量(或產(chǎn)量)影響就越大；

5)表皮層徑向滲透系數(shù)KIr為較小值時，承壓含水層的水頭降幅度主要取決于承壓含水層徑向滲透系數(shù)KIIr的大小，而承壓含水層豎向滲透系數(shù)KIIv對承壓含水層水頭降影響較??；但是當(dāng)表皮層徑向滲透系數(shù)KIr較大時，承壓含水層徑向滲透系數(shù)KIIr和豎向滲透系數(shù)KIIv均對承壓含水層的水頭降有較大影響．

參考文獻(xiàn)：

[1]Zhang W,Zhan H B,Huang G H,et al.Constant-head test in a leaky aquifer with a finite-thickness skin[J].Journal of Hydrology,2011,399(3):326-334

[2]Chen Y J,Yeh H D.Parameter estimation/sensitivity analysis for an aquifer test with skin effect[J].Ground Water,2009,47(2):287-299

[3]Yan S Y,Yeh H D.Laplace-domain solutions for radial two-zone flow equations under the conditions of constant-head and partially penetrating well[J].Journal of Hydraulic Engineering,2005,131(3):209-216

[4]Chen C S,Chang C C.Well hydraulics theory and data analysis of the constant head test in an unconfined aquifer with the skin effect[J].Water Resources Research,2003,39(5):1121-1135

[5]Chen C S,Chang C C.Theoretical evaluation of non-uniform skin effect on aquifer response under constant rate pumping[J].Journal of Hydrology,2006,317(3):190-201

[6]Yeh H D,Yang S Y,Peng H Y.A new closed-form solution for a radial two-layer drawdown equation for groundwater under constant-flux pumping in a finite-radius well[J].Advances in Water Resources,2003,26(7):747-757

[7]Yeh H D,Yang S Y.A novel analytical solution for a slug test conducted in a well with a finite-thickness skin[J].Advances in Water Resources,2006,29(10):1479-1489

[8]Yeh H D,Chen Y J,Yang S Y.Semi-analytical solution for a slug test in partially penetrating wells including the effect of finite-thickness skin[J].Hydrological Processes,2008,22(18):3741-3748

[9]Yan S Y,Yeh H D.Solution for flow rates across the wellbore in a two-zone confined aquifer[J].Journal of Hydraulic Engineering,2002,128(2):175-183

[10]Stehfest H.Remark on Algorithm 368:numerical inversion of Laplace transforms[J].Communications of the Association for Computing Machinery,1970,13(10):624-625

[11]Perina T,Lee T C.General well function for pumping from a confined,leaky,or unconfined aquifer[J].Journal of Hydrology,2006,317(3):239-260

[12]Hemker C J.Transient well flow in leaky multiple-aquifer systems[J].Journal of Hydrology,1985,81(1):111-126