李 旭，蘇世林，文 章，許光泉

(1.安徽理工大學(xué)地球與環(huán)境學(xué)院，安徽 淮南 232001；2.中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北 武漢 430078)

裂隙介質(zhì)廣泛存在于自然界中，與人類的生產(chǎn)生活活動(dòng)密切相關(guān)。裂隙介質(zhì)中地下水滲流和溶質(zhì)運(yùn)移規(guī)律的研究一直是國(guó)內(nèi)外學(xué)者重點(diǎn)關(guān)注的問(wèn)題之一，但是由于裂隙介質(zhì)具有強(qiáng)烈的空間變異性，使得溶質(zhì)運(yùn)移特性研究變得十分復(fù)雜[1-4]。近幾十年來(lái)，國(guó)內(nèi)外學(xué)者對(duì)此開(kāi)展了大量的研究，并取得了一定的研究進(jìn)展，其主要研究成果可分為兩大類:第一類為裂隙網(wǎng)格中水流和溶質(zhì)運(yùn)移理論研究[5-6]；第二類為單裂隙介質(zhì)中水流和溶質(zhì)運(yùn)移理論研究[7-10]。單裂隙作為裂隙介質(zhì)水流和溶質(zhì)運(yùn)移研究的基礎(chǔ)單元，多數(shù)研究通常是從單裂隙開(kāi)始的。如鄭志成等[8]利用平行大理石板開(kāi)展了單個(gè)裂隙溶質(zhì)運(yùn)移實(shí)驗(yàn)研究，結(jié)果表明溶質(zhì)的穿透曲線表現(xiàn)出拖尾的現(xiàn)象。因此，開(kāi)展單裂隙溶質(zhì)運(yùn)移理論研究對(duì)于探索裂隙介質(zhì)中溶質(zhì)運(yùn)移規(guī)律具有重要的理論意義。

通常情況下，單裂隙溶質(zhì)運(yùn)移理論研究主要是在達(dá)西流的基礎(chǔ)上進(jìn)行的[11-12]，然而越來(lái)越多的研究表明達(dá)西定律只在一定的水力梯度范圍內(nèi)成立，裂隙介質(zhì)中地下水流速過(guò)快往往會(huì)導(dǎo)致地下水滲流速度與水力梯度呈非線性關(guān)系，即所謂的非達(dá)西流[13-16]。Qian等[17]通過(guò)單裂隙室內(nèi)實(shí)驗(yàn)研究發(fā)現(xiàn)，單裂隙介質(zhì)中地下水平均流速和水力梯度的關(guān)系能很好地與Izbash方程吻合，且水流速度與水力梯度的1/2次方成正比；李一鳴等[18]利用Izbash方程刻畫(huà)了單裂隙介質(zhì)中的非達(dá)西流，研究當(dāng)單裂隙的軸向與基巖水流方向斜交過(guò)程中，裂隙流的非達(dá)西程度對(duì)溶質(zhì)羽分布的影響。上述研究表明，Izbash方程能夠很好地描述單裂隙介質(zhì)中的非達(dá)西流，可以應(yīng)用于單裂隙介質(zhì)中溶質(zhì)運(yùn)移機(jī)理的研究。

目前國(guó)內(nèi)外學(xué)者對(duì)單裂隙的非達(dá)西滲流和溶質(zhì)運(yùn)移開(kāi)展了大量的理論研究，其研究通常假定地下水流是穩(wěn)定的[19-21]。然而，實(shí)際裂隙含水層由于受到補(bǔ)給/排泄模式、地表水體水位以及區(qū)域地下水抽采率的變化等因素的影響，導(dǎo)致裂隙介質(zhì)中地下水流速是隨時(shí)間發(fā)生變化的[22]。目前許多學(xué)者也開(kāi)展了變流速條件下的溶質(zhì)運(yùn)移理論研究[22-24]。如Singh等[24]利用拉普拉斯變換得到了地下水流速為正弦變化的二維溶質(zhì)運(yùn)移的解析解；Li等[22]以松散含水層為研究對(duì)象，提出了指數(shù)變化的地下水流速方程，并利用積分變換獲得變流速淋濾作用下的一維溶質(zhì)運(yùn)移的解析解。目前關(guān)于變流速條件下的溶質(zhì)運(yùn)移理論研究對(duì)象主要為松散含水層，而裂隙含水層中的變流速溶質(zhì)運(yùn)移研究幾乎還是空白。此外，裂隙介質(zhì)溶質(zhì)運(yùn)移過(guò)程中的地下水流時(shí)常不符合達(dá)西定律。因此，有必要開(kāi)展變流速條件下非達(dá)西裂隙流溶質(zhì)運(yùn)移規(guī)律的研究。

為此，本文將建立變流速條件下考慮非達(dá)西流的單裂隙一維溶質(zhì)運(yùn)移模型。利用指數(shù)變化的地下水流速方程來(lái)刻畫(huà)單裂隙介質(zhì)中地下水流速的變化規(guī)律，并耦合Izbash方程，利用積分變換的方法獲得變流速條件下單裂隙非達(dá)西流溶質(zhì)運(yùn)移的解析解，以此來(lái)探究非達(dá)西流和指數(shù)變化的地下水流速對(duì)裂隙介質(zhì)中溶質(zhì)運(yùn)移的影響機(jī)理，以為裂隙地下水污染研究提供理論依據(jù)。

1 非達(dá)西滲流與變流速水流

1.1 裂隙達(dá)西流與非達(dá)西流

在裂隙介質(zhì)中，當(dāng)?shù)叵滤魉佥^小時(shí)，黏滯力占主導(dǎo)地位，慣性力的影響可以忽略，通常用線性的達(dá)西定律來(lái)表示，即：

vd=-KJ

(1)

式中：vd為達(dá)西滲流速度[L/T]；J為水力梯度[無(wú)量綱]；K為裂隙介質(zhì)的滲透系數(shù)[L/T]。

但當(dāng)?shù)叵滤魉佥^大、慣性力占主導(dǎo)地位時(shí)，地下水流速與水力梯度之間不再呈現(xiàn)線性的關(guān)系，發(fā)生非達(dá)西滲流。為此，國(guó)內(nèi)外學(xué)者提出了不同類型的非線性流方程，在方程中因參數(shù)太多，限制了其應(yīng)用[13]。本文采用Izbash方程來(lái)描述裂隙介質(zhì)中的非達(dá)西滲流[18]，其數(shù)學(xué)表達(dá)式為

vn=-kJ

(2)

式中：n為經(jīng)驗(yàn)系數(shù)，n值反映了裂隙水流流態(tài)的變化，即從達(dá)西流(n=1)到非達(dá)西流(1

裂隙水流動(dòng)狀態(tài)依據(jù)雷諾數(shù)(Re)判定，其判定方程為

(3)

式中：2b為裂隙的平均開(kāi)度[L]；q為裂隙介質(zhì)中地下水平均流速[L/T]；μ為地下水動(dòng)力黏度[L2/T]。

通常情況下，裂隙中水流的Re大于1～10之間的某個(gè)值時(shí)，地下水滲流速度與水力梯度之間不再呈現(xiàn)線性的關(guān)系，而表現(xiàn)出非達(dá)西滲流的現(xiàn)象。

1.2 變流速水流

目前關(guān)于裂隙介質(zhì)溶質(zhì)運(yùn)移模型中地下水流速通常設(shè)定為常數(shù)，但實(shí)際過(guò)程中地下水流速是隨著時(shí)間的變化而發(fā)生變化的(見(jiàn)圖1)，針對(duì)這一問(wèn)題，提出了采用線性方程和指數(shù)變化的方程來(lái)刻畫(huà)地下水流速的變化規(guī)律。實(shí)際上裂隙介質(zhì)中地下水流速的變化應(yīng)趨于某一確定的值，表現(xiàn)出指數(shù)增加或指數(shù)衰減的趨勢(shì)[22]。例如，當(dāng)含水層受到生物、化學(xué)及物理的堵塞會(huì)導(dǎo)致其滲透性減小，使得含水層中地下水流速表現(xiàn)出指數(shù)衰減[25-26]。再如，河岸帶及海岸帶的裂隙含水層由于地表水位的波動(dòng)(如洪水、潮汐等作用)也會(huì)引起含水層中地下水流速呈現(xiàn)出指數(shù)增加或指數(shù)衰減的趨勢(shì)[23]。因此，在達(dá)西流條件下的裂隙含水層中地下水流速可以表示為指數(shù)增加或指數(shù)衰減的形式[22]：

圖1 變流速條件下非達(dá)西裂隙流溶質(zhì)運(yùn)移示意圖Fig.1 Schematic diagram of solute transport under Non-Darcian flow and variable flow velocity in the fracture

vd(t)=vd0φ(t)=vd1+(vd0-vd1)(e-λt)

(4)

式中：vd0為單裂隙介質(zhì)中初始的地下水流速[L/T]；vd1為單裂隙介質(zhì)中最終穩(wěn)定的地下水流速[L/T]；λ為地下水流速變化指數(shù)[1/T]。

當(dāng)vd0>vd1時(shí)，公式(4)表示地下水流速為指數(shù)衰減的情況；當(dāng)vd0

[v(t)]n=[v0τ(t)]n=(v1)n+[(v0)n-(v1)n](e-λt)

(5)

式中：v0為非達(dá)西流條件下單裂隙介質(zhì)中初始的地下水流速[L/T]；v1為非達(dá)西流條件下最終穩(wěn)定的地下水流速[L/T]。

若n=1時(shí)，公式(5)與公式(4)相同，則τ(t)可表示為

(6)

2 數(shù)學(xué)模型及求解

在建立變流速條件下非達(dá)西裂隙流一維溶質(zhì)運(yùn)移模型之前，為了簡(jiǎn)化數(shù)學(xué)模型，假設(shè)條件為：①裂隙為單裂隙，其水流為非達(dá)西一維流動(dòng)；②裂隙的開(kāi)度不隨距離變化而變化，不考慮裂隙的粗糙度及充填物等因素的影響；③地下水流為非穩(wěn)定流，其流速隨時(shí)間呈指數(shù)變化；④溶質(zhì)為惰性溶質(zhì)。則單裂隙介質(zhì)中對(duì)流-彌散方程(ADE)可以表示為

(7)

式中：C為單裂隙介質(zhì)中溶質(zhì)的濃度[M/L3]；x為距離[L]；t為時(shí)間[T]；D(t)為隨時(shí)間變化的彌散系數(shù)[L2/T]，v(t)為隨時(shí)間變化的地下水流速[L/T]，兩者分別表示如下：

v(t)=v0τ(t)

D(t)=αv(t)=αv0τ(t)=D0τ(t)

(8)

式中：α為彌散度[L]。

將公式(8)代入到公式(7)中，可得：

(9)

初始條件和邊界條件可表示如下：

C(r,0)=0,x<0

(10)

C(0,t)=C0,t>0

(11)

C(∞,t)=0,t>0

(12)

式中：C0為x=0處給定的溶質(zhì)濃度[M/L3]。

為了獲取模型的解，這里引入一個(gè)新的積分變換：

(13)

則公式(9)可以表示為

(14)

邊界條件和初始條件在相同的積分變化下可以表示如下：

C(x,0)=0,x<0

(15)

C(0,T)=C0,T>0

(16)

C(∞,T)=0,T>0

(17)

通過(guò)積分變換，公式(13)～(16)為經(jīng)典的ADE模型。因此，一維半無(wú)限長(zhǎng)單裂隙介質(zhì)在x-T域中定濃度邊界條件下ADE模型的解為

(18)

式中：erfc(·)為余誤差函數(shù)。

將公式(13)代入到公式(17)中，可獲得一維半無(wú)限長(zhǎng)單裂隙介質(zhì)在x-t域中定濃度邊界條件下ADE模型的解為

(19)

上述推導(dǎo)的公式(19)為變流速條件下非達(dá)西裂隙流溶質(zhì)運(yùn)移模型的解析解。為了進(jìn)一步驗(yàn)證該模型的準(zhǔn)確性，將不考慮非達(dá)西流的溶質(zhì)運(yùn)移結(jié)果,將本文模型解析解與Li等[22]的模型解析解進(jìn)行了對(duì)比。Li等[22]提出了指數(shù)變化的地下水流速方程，并獲得了變流速條件下一維溶質(zhì)運(yùn)移模型的解析解，但是該模型沒(méi)有考慮非達(dá)西流的影響。當(dāng)非達(dá)西參數(shù)n=1時(shí)，本文模型解析解結(jié)果應(yīng)該與Li等[22]模型解析解結(jié)果一致。圖2為本文與Li等[22]模型解析解在x=20 m處的溶質(zhì)穿透曲線(BTC)對(duì)比結(jié)果。模型主要參數(shù)設(shè)置如下：n=1，v0=1×10-3m/s，v1=6×10-3m/s，λ=0.000 01 1/s，α=0.1 m、0.5 m和2 m。

圖2 本文與Li等[22]模型解析解在x=20 m處的溶質(zhì) 穿透曲線對(duì)比Fig.2 Comparison of the breakthrough curves at x=20 m computed by analytical solution of this study and Li et al.[22]

由圖2可見(jiàn)，本文模型解析解結(jié)果與Li等[22]模型解析解結(jié)果完全吻合，進(jìn)一步說(shuō)明本文的解析模型是準(zhǔn)確的。

3 結(jié)果分析與討論

3.1 非達(dá)西裂隙流的溶質(zhì)運(yùn)移特征

為了研究非達(dá)西流對(duì)單裂隙介質(zhì)中溶質(zhì)運(yùn)移的影響，首先應(yīng)確定非達(dá)西裂隙流地下水流速的范圍。本文將發(fā)生非達(dá)西流與達(dá)西流流態(tài)轉(zhuǎn)變的臨界雷諾數(shù)設(shè)定為10，由公式(3)可得到單裂隙介質(zhì)中臨界地下水流速的計(jì)算公式為

(20)

將參數(shù)設(shè)置為：2b=0.02 m，μ=1.0×10-6m2/s[18]，根據(jù)公式(20)，可計(jì)算得到單裂隙介質(zhì)中臨界地下水流速為0.000 5 m/s。

當(dāng)裂隙介質(zhì)中地下水流速為指數(shù)衰減條件下，在距離原始溶質(zhì)注入點(diǎn)x=40 m處，當(dāng)非達(dá)西參數(shù)n取不同值時(shí)的溶質(zhì)穿透曲線，見(jiàn)圖3。模型主要參數(shù)設(shè)置如下：v0=0.001 5 m/s，v1=0.000 5 m/s，λ=0.000 05 1/s，α=0.1 m，n=1、1.2、1.5和2.0。

圖3 變流速條件下不同非達(dá)西參數(shù)n在x=40 m處 所對(duì)應(yīng)的溶質(zhì)穿透曲線Fig.3 Breakthrough curves for different Non-Darcian parameters (n) at x=40 m under variable flow velocity

由圖3可見(jiàn)，溶質(zhì)穿透曲線的溶質(zhì)濃度隨著非達(dá)西參數(shù)n的增大而不斷增大，表明非達(dá)西參數(shù)n越大，地下水紊流程度越高，溶質(zhì)運(yùn)移速度越快。

為了進(jìn)一步分析導(dǎo)致溶質(zhì)運(yùn)移速度變快的原因，本文研究了不同非達(dá)西參數(shù)n時(shí)單裂隙介質(zhì)中地下水流速v的變化，其變化曲線見(jiàn)圖4。

圖4 不同非達(dá)西參數(shù)n下單裂隙介質(zhì)中地下水流 速的變化曲線Fig.4 Velocity curves of the groundwater flow in a single fracture for different Non-Darcian parameters (n)

由圖4可見(jiàn)，當(dāng)裂隙介質(zhì)中地下水流速由0.001 5 m/s減小到0.000 5 m/s時(shí)，非達(dá)西參數(shù)n越大，地下水流速衰減得越慢，地下水流速越快?？傊谧兞魉贄l件下裂隙介質(zhì)中的溶質(zhì)運(yùn)移受到非達(dá)西流的影響明顯，非達(dá)西流的存在改變了裂隙介質(zhì)中地下水流速的變化規(guī)律。因此，在裂隙介質(zhì)中刻畫(huà)變流速溶質(zhì)運(yùn)移過(guò)程不能忽視非達(dá)西流因素的影響。

3.2 變流速條件下非達(dá)西流裂隙流的溶質(zhì)運(yùn)移特征

非達(dá)西流條件下地下水流速衰減指數(shù)λ取不同值時(shí)在x=40 m處所對(duì)應(yīng)的溶質(zhì)穿透曲線，見(jiàn)圖5。模型主要參數(shù)設(shè)置如下：v0=0.001 5 m/s，v1= 0.000 5 m/s，α=0.2 m，n=1.5，λ=0.000 025 1/s、0.000 050 1/s和0.000 075 1/s。對(duì)于指數(shù)衰減的地下水流速，λ值的大小反映了地下水流速衰減得快慢，主要是由地下水水位變化速率或者滲透性衰減得快慢所決定的。

圖5 非達(dá)西流條件下不同取值的地下水流速衰減指數(shù) λ在x=40 m處所對(duì)應(yīng)的溶質(zhì)穿透曲線Fig.5 Breakthrough curves for different attenuation indexes λ of the exponentially decreasing groundwater velocity at x=40 m under Non-Darcian flow

由圖5可見(jiàn)，地下水流速衰減指數(shù)λ越大，溶質(zhì)穿透曲線的溶質(zhì)濃度越小，說(shuō)明地下水流速衰減指數(shù)λ越大，溶質(zhì)運(yùn)移速度越慢。這主要是由于地下水流速衰減指數(shù)λ越大，地下水流速衰減得越快，導(dǎo)致溶質(zhì)運(yùn)移速度越慢。當(dāng)v1=0.001 5 m/s時(shí)，地下水為穩(wěn)定流，本文通過(guò)對(duì)比變流速條件下的溶質(zhì)穿透曲線與穩(wěn)定流的溶質(zhì)穿透曲線可以發(fā)現(xiàn)，變流速條件下溶質(zhì)穿透曲線的溶質(zhì)濃度值均小于地下水穩(wěn)定流的情況，說(shuō)明地下水流速的變化對(duì)裂隙介質(zhì)中溶質(zhì)運(yùn)移有較大的影響。

非達(dá)西流條件下地下水流速v1取不同值時(shí)在x=40 m處所對(duì)應(yīng)的溶質(zhì)穿透曲線，見(jiàn)圖6。模型主要參數(shù)設(shè)置如下：v0=0.001 5 m/s，λ=0.000 050 1/s，α=0.2 m，n=1.5，v1=0.000 7 m/s、0.001 1 m/s、0.001 5 m/s、0.001 9 m/s和0.002 3 m/s。當(dāng)?shù)叵滤魉賤1取不同值時(shí)，公式(2)可以表示為地下水流速指數(shù)增加或指數(shù)衰減的函數(shù)形式。因此，當(dāng)v1>0.001 5 m/s 時(shí)，地下水流速為指數(shù)增加的函數(shù)形式；當(dāng)v1<0.001 5 m/s 時(shí)，地下水流速為指數(shù)衰減的形式；當(dāng)v1=0.001 5 m/s 時(shí)，地下水為穩(wěn)定流。

圖6 非達(dá)西流條件下不同v1取值在x=40 m處所 對(duì)應(yīng)的溶質(zhì)穿透曲線Fig.6 Breakthrough curves for different final steady velocities (v1) at x=40 m under Non-Darcian flow

由圖6可見(jiàn)，地下水流速v1越小，地下水流速衰減的幅度越大，溶質(zhì)穿透曲線的溶質(zhì)濃度越小，溶質(zhì)運(yùn)移速度越慢。例如，當(dāng)v1=0.000 7 m/s時(shí)，溶質(zhì)運(yùn)移速度最慢，其溶質(zhì)穿透曲線位于最右端(見(jiàn)圖6)。另外，地下水流速v1越大，地下水流速增加的幅度越大，溶質(zhì)穿透曲線的溶質(zhì)濃度越大，溶質(zhì)運(yùn)移速度越快。因此，利用公式(2)可以刻畫(huà)裂隙介質(zhì)中地下水流速的指數(shù)增加或指數(shù)衰減的情況，為研究地下水流速變化對(duì)溶質(zhì)運(yùn)移的影響提供了有效途徑。

4 結(jié) 論

本文依據(jù)Izbash和指數(shù)變化的地下水流速方程建立了單裂隙介質(zhì)非達(dá)西流溶質(zhì)運(yùn)移新模型，通過(guò)積分變換的方法得到了模型的解析解，并模擬分析了非達(dá)西流和指數(shù)變化的地下水流速對(duì)裂隙介質(zhì)中溶質(zhì)運(yùn)移的影響，得到如下結(jié)論：

(1) 變流速條件下，單裂隙介質(zhì)中非達(dá)西參數(shù)n越大，地下水紊流越強(qiáng)，地下水流速越大，導(dǎo)致溶質(zhì)運(yùn)移速度越快，溶質(zhì)穿透曲線的溶質(zhì)濃度越高。因此，單裂隙介質(zhì)中刻畫(huà)變流速溶質(zhì)運(yùn)移過(guò)程不能忽視非達(dá)西流因素的影響。

(2) 單裂隙介質(zhì)中指數(shù)變化的地下水流速對(duì)溶質(zhì)運(yùn)移有較大的影響，地下水流速的衰減指數(shù)λ越大，地下水流速衰減得越快，溶質(zhì)穿透曲線的溶質(zhì)濃度越低，溶質(zhì)運(yùn)移速度越慢。

(3) 地下水流速(或漸進(jìn)地下水流速)在單裂隙介質(zhì)溶質(zhì)運(yùn)移過(guò)程中起到了關(guān)鍵的作用，其大小反映了地下水流速增加或減小的幅度，其變化的幅度越大，溶質(zhì)運(yùn)移速度快慢表現(xiàn)得越明顯。