基于時(shí)空分?jǐn)?shù)階導(dǎo)數(shù)模型的潛流帶溶質(zhì)運(yùn)移模擬

魯程鵬，林雨竹，張勇，秦巍，吳成城，劉波，束龍倉(cāng)

(1.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098；2.Department of Geological Sciences of University of Alabama,Tuscaloosa USA AL 35487)

潛流帶是河道內(nèi)地表地下水相互作用的飽水沉積層，連接河道及地下水體，三者的生物、物理、化學(xué)特性彼此影響[1]。潛流帶中發(fā)生的物質(zhì)交換常被認(rèn)為是影響河道污染物遷移及能量交換的基本過(guò)程。河道中的溶質(zhì)粒子常隨水流進(jìn)入潛流帶停留一段時(shí)間，最終在下游某些點(diǎn)返回河流。潛流帶中的溶質(zhì)運(yùn)移很好地解釋了河道地表及地下水體間的物質(zhì)聯(lián)系，因此探究潛流帶中的溶質(zhì)遷移規(guī)律對(duì)于進(jìn)一步了解河道內(nèi)反應(yīng)性運(yùn)輸和生物吸收十分重要[2]。

潛流帶沉積物在空間上不均勻分布以及河道的物理化學(xué)性質(zhì)、生物地球化學(xué)過(guò)程、溫度變化、地表-地下水位波動(dòng)都會(huì)對(duì)河床介質(zhì)的非均質(zhì)性產(chǎn)生影響，滲透系數(shù)的空間差異性是這種非均質(zhì)性的重要表征：Tang等[3]設(shè)置了不同介質(zhì)條件模擬非均質(zhì)介質(zhì)中飽和流的水流通量交換情況，發(fā)現(xiàn)滲透系數(shù)的差異會(huì)對(duì)模擬結(jié)果影響較大；Wang等[4-6]發(fā)現(xiàn)天然河道中水平面內(nèi)順?biāo)鞣较虻暮哟矟B透系數(shù)沒(méi)有明顯變化，而縱深方向介質(zhì)具有較強(qiáng)非均質(zhì)性；Ren等[7]發(fā)現(xiàn)滲透系數(shù)對(duì)孔隙度變化十分敏感，具體表現(xiàn)在大孔隙是過(guò)流的主要通道，但較大的孔隙通道在壓實(shí)過(guò)程中總是被提前關(guān)閉；Leung 等[8]發(fā)現(xiàn)河床溫度可以影響水的黏度進(jìn)而改變介質(zhì)滲透系數(shù)；Zhu等[9]研究發(fā)現(xiàn)微生物空間中生長(zhǎng)產(chǎn)生阻塞作用會(huì)降低局部介質(zhì)的滲透性；活性礦物質(zhì)或離子的氧化和凝膠狀物質(zhì)的形成同樣也會(huì)通過(guò)影響河床電位梯度控制潛流交換的速率[10]。當(dāng)前針對(duì)模擬潛流帶溶質(zhì)過(guò)程的模型已經(jīng)十分豐富： Singh等[11]提出了一種降階模型用于模擬流量變化和河床地形相互作用影響下的潛流交換特征；Adu[12]運(yùn)用解析法構(gòu)造了一種非點(diǎn)源混合單元模型，模擬得到的穿透曲線與現(xiàn)場(chǎng)示蹤劑測(cè)試結(jié)果擬合程度較好；Sherman等[13]構(gòu)建了隨機(jī)游走-拉格朗日模型提高了溶質(zhì)粒子在介質(zhì)中的遷移運(yùn)動(dòng)模擬的精確度；對(duì)流-彌散方程作為地下水溶質(zhì)運(yùn)移的基礎(chǔ)理論方程也依然被廣泛應(yīng)用于示蹤模擬。但現(xiàn)有方法往往忽略了河床介質(zhì)的非均質(zhì)性，增大了對(duì)污染物遷移模擬的不準(zhǔn)確性[5]，近幾年出現(xiàn)的分?jǐn)?shù)階方法在數(shù)學(xué)層面上解釋了由介質(zhì)非均質(zhì)性引起的溶質(zhì)粒子隨機(jī)運(yùn)動(dòng)問(wèn)題，將影響粒子遷移的因素轉(zhuǎn)化為時(shí)間、空間分?jǐn)?shù)階導(dǎo)數(shù)綜合判斷，具有良好的模擬效果。

基于Fick定律得到的經(jīng)典的對(duì)流-彌散方程已被廣泛運(yùn)用于地下水溶質(zhì)運(yùn)移模擬，但針對(duì)潛流帶這一結(jié)構(gòu)特殊、介質(zhì)非均質(zhì)性強(qiáng)且水流方向復(fù)雜的多孔介質(zhì)，溶質(zhì)易呈現(xiàn)拖尾及非高斯分布的反常擴(kuò)散特征，傳統(tǒng)對(duì)流-彌散方程模型很難準(zhǔn)確描述這類非局域性特征[14]。分?jǐn)?shù)階模型將偏微分方程中求導(dǎo)的整數(shù)階數(shù)變?yōu)榉謹(jǐn)?shù)，已被證明在描述物質(zhì)的反常擴(kuò)散方面具有物理意義明確、可描述介質(zhì)非均質(zhì)性等優(yōu)勢(shì)[14-15]。當(dāng)前針對(duì)分?jǐn)?shù)階方法的研究多集中于算法的開(kāi)發(fā)(數(shù)值穩(wěn)定、快速求解)[16-20]，而在分?jǐn)?shù)階的應(yīng)用領(lǐng)域的探索主要集中于模擬飽和熱傳導(dǎo)問(wèn)題[21-22]、非飽和土壤水分?jǐn)U散[23-24]、土柱沙箱試驗(yàn)[25-26]、泥沙運(yùn)動(dòng)等方面[27-28]。本文將時(shí)空分?jǐn)?shù)階模型首次引入潛流帶溶質(zhì)運(yùn)移問(wèn)題，與傳統(tǒng)整數(shù)階模型結(jié)果進(jìn)行對(duì)比，試圖探究分?jǐn)?shù)階方法對(duì)應(yīng)實(shí)際問(wèn)題的物理意義以及對(duì)于水流介質(zhì)參數(shù)的靈敏性，最后結(jié)合野外現(xiàn)場(chǎng)示蹤試驗(yàn)數(shù)據(jù)，分別在一維和二維兩個(gè)維度探索分?jǐn)?shù)階方法在潛流帶溶質(zhì)運(yùn)移方面的應(yīng)用效果與前景。

1 時(shí)空分?jǐn)?shù)階模型的建立及數(shù)值解法

1.1 模型建立

在傳統(tǒng)對(duì)流-彌散方程基礎(chǔ)上分別引入時(shí)間、空間分?jǐn)?shù)階導(dǎo)數(shù)，建立時(shí)空分?jǐn)?shù)階二維對(duì)流-彌散方程

(1)

式中：C(x,y,t)為溶質(zhì)濃度；V和D分別為介質(zhì)中的平均流速和彌散系數(shù)；Rd表示阻滯因子；f(x,y,t)表示源匯項(xiàng)。α(0<α?1)為時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)，使用Caputo式求解：當(dāng)α=1時(shí)，方程化為空間分?jǐn)?shù)階對(duì)流-彌散方程；β(1<β?2)為空間分?jǐn)?shù)階導(dǎo)數(shù)，使用Riemman-Liouville式求解：當(dāng)β=2時(shí)，方程化為時(shí)間分?jǐn)?shù)階對(duì)流-彌散方程；當(dāng)α=1,β=2時(shí)，方程簡(jiǎn)化為傳統(tǒng)對(duì)流-彌散方程。設(shè)置初始條件和邊界條件分別為

C(x,y,0)=Ψ(x)

C(0,0,t)=0,C(x,0,t)=φ1(t),C(0,y,t)=φ2(t)

(2)

式中：x、y和t分別為x、y和t方向的邊界。

1.2 數(shù)值解法

采用Liu[29]開(kāi)發(fā)的已被證明穩(wěn)定收斂的一維時(shí)空分?jǐn)?shù)階數(shù)值解法，并采用穩(wěn)定的差分法將方程擴(kuò)展至二維條件[30-31]，實(shí)現(xiàn)二維時(shí)空分?jǐn)?shù)階數(shù)值解的運(yùn)算。

分?jǐn)?shù)階導(dǎo)數(shù)分別離散為

(3)

O(hx)

(4)

O(hy)

(5)

其中，bp=(p+1)1-α-p1-α,j=0,1,2,…，tn。

p=1,2,…，i+1。

(6)

q=1,2,…,j+1。

(7)

代入式(1)得到隱式差分格式為

(8)

其中，i=1,2,…,m;j=1,2…,n;k=0,1,2,…,tn-1。

得到時(shí)空分?jǐn)?shù)階方程的數(shù)值解為

(9)

2 分?jǐn)?shù)階階數(shù)對(duì)溶質(zhì)運(yùn)移過(guò)程的影響

2.1 分?jǐn)?shù)階階數(shù)物理意義

參考周璐瑩等[31]對(duì)理想算例的討論方法，設(shè)置理想條件(溶質(zhì)濃度為均一化濃度)：一維空間(0，10 m)，計(jì)算時(shí)間100 min。初始時(shí)刻，在x=60 m瞬時(shí)釋放100單位點(diǎn)源溶質(zhì)，取彌散系數(shù)D=0.15 m2/min，對(duì)流流速v=0.01 m/min。

模型中固定β=2，分別設(shè)置α=0.2、0.4、0.6、0.8、1.0，計(jì)算x=6 m處的穿透曲線以及t=15 min時(shí)的溶質(zhì)擴(kuò)散曲線，研究α的變化對(duì)溶質(zhì)運(yùn)移的影響，結(jié)果見(jiàn)圖1。由圖1可以看出：在穿透曲線圖1(a)中，隨α減小，穿透曲線的峰現(xiàn)時(shí)間滯后，峰值變小，同時(shí)濃度的下降速度也減緩，穿透曲線的偏態(tài)性增強(qiáng)，擴(kuò)散穩(wěn)定后溶質(zhì)濃度高，即α減小擴(kuò)散拖尾現(xiàn)象越明顯；在溶質(zhì)擴(kuò)散曲線圖1(b)中，隨α減小，溶質(zhì)的峰向上游移動(dòng)，擴(kuò)散減緩，擴(kuò)散曲線的偏態(tài)性增強(qiáng)，擴(kuò)散趨于穩(wěn)定后溶質(zhì)的濃度變高，這也與圖1(a)中反映的擴(kuò)散曲線特征吻合。

圖1 不同α和β在x=6 m處的歸一化穿透曲線及t=15 min時(shí)刻的歸一化溶質(zhì)擴(kuò)散曲線Fig.1 Breakthrough curves of different α and β at x=6 m and solute diffusion curves at t=15 min

結(jié)合擴(kuò)散行為的物理意義及分?jǐn)?shù)階階數(shù)的物理意義：時(shí)間分?jǐn)?shù)階α表征擴(kuò)散行為的歷史記憶性，即歷史上某點(diǎn)所有時(shí)刻的溶質(zhì)濃度均對(duì)當(dāng)前濃度產(chǎn)生影響，也即擴(kuò)散在非均質(zhì)介質(zhì)中的滯時(shí)效應(yīng)。α值越小，溶質(zhì)在介質(zhì)中的滯時(shí)效應(yīng)越大，擴(kuò)散減緩，穿透曲線趨于扁平，擴(kuò)散達(dá)到穩(wěn)定時(shí)的溶質(zhì)濃度也越大。歷史記憶性使得反常擴(kuò)散比正常的擴(kuò)散過(guò)程更慢。

在前述計(jì)算條件下，固定α=1，分別設(shè)置β=1.2、1.4、1.6、1.8、2.0，計(jì)算x=6 m處的穿透曲線以及t=15 min時(shí)的溶質(zhì)擴(kuò)散曲線，研究β的變化對(duì)溶質(zhì)運(yùn)移的影響，結(jié)果如見(jiàn)圖1(c)、1(d)。可觀察到：在圖1(c)中，隨β增大，穿透曲線的峰現(xiàn)時(shí)間及趨于穩(wěn)定的時(shí)間均延后，峰值變大，同時(shí)濃度的下降速度也增加，即β增大，溶質(zhì)的擴(kuò)散過(guò)程加快；在圖1(d)中，隨β減小，溶質(zhì)的峰向上游移動(dòng)，且峰值變大，即溶質(zhì)的擴(kuò)散過(guò)程減緩，擴(kuò)散趨于穩(wěn)定后溶質(zhì)的濃度增大，這也與圖1(c)中反映的擴(kuò)散曲線特征吻合。

結(jié)合擴(kuò)散行為的物理意義及分?jǐn)?shù)階階數(shù)的物理意義：空間分?jǐn)?shù)階β表征擴(kuò)散行為的空間非局域性，即同一時(shí)刻空間上的所有點(diǎn)(而非相鄰點(diǎn))都對(duì)某點(diǎn)的濃度產(chǎn)生影響，也即擴(kuò)散在非均質(zhì)介質(zhì)中的尺度效應(yīng)。這種非局域性使得擴(kuò)散過(guò)程加快，β值越小，溶質(zhì)擴(kuò)散的非局域性越強(qiáng)，非均質(zhì)介質(zhì)中形成的快速溶質(zhì)通道對(duì)于計(jì)算結(jié)果的影響越明顯，從而導(dǎo)致擴(kuò)散加快，峰值減小，峰現(xiàn)時(shí)間提前，溶質(zhì)濃度越早趨于穩(wěn)定，擴(kuò)散達(dá)到穩(wěn)定時(shí)的溶質(zhì)濃度也越小。非局域性使得反常擴(kuò)散比正常的擴(kuò)散過(guò)程更快。

2.2 分?jǐn)?shù)階對(duì)溶質(zhì)運(yùn)移物理參數(shù)的影響

在前述理想計(jì)算條件下設(shè)置不同的流速v、彌散系數(shù)D，對(duì)比分?jǐn)?shù)階與整數(shù)階方法對(duì)溶質(zhì)運(yùn)移過(guò)程中物理參數(shù)改變的影響。圖2(a)、2(b)設(shè)置彌散系數(shù)D=0.15 m2/min，分別計(jì)算引入分?jǐn)?shù)階前后(α=1.0、0.6、0.2，β=2.0、1.6、1.2)改變流速(v=0、0.01、0.02、0.03 m/min)對(duì)穿透曲線形態(tài)的影響；圖2(c)、2(d)設(shè)置流速v=0.01，分別計(jì)算引入分?jǐn)?shù)階前后(α=1.0、0.6、0.2，β=2.0、1.6、1.2)改變彌散系數(shù)(D=0.10、0.15、0.20、0.25 m2/min)對(duì)穿透曲線形態(tài)的影響。

圖2 α和β對(duì)對(duì)流彌散方程系數(shù)(v、D)敏感性分析Fig.2 Comparison of sensitivity of different α and β to velocity v and Dispersion coefficient D

分?jǐn)?shù)階導(dǎo)數(shù)對(duì)彌散系數(shù)改變的影響見(jiàn)圖2(a)、2(b)。當(dāng)α、β確定時(shí)，增大流速v主要導(dǎo)致穿透曲線下降階段斜率增加，而峰現(xiàn)時(shí)間幾乎不變。這是由于流速為時(shí)間相關(guān)參數(shù)，主要影響穿透曲線的滯時(shí)效應(yīng)(即拖尾現(xiàn)象)，而對(duì)空間差異導(dǎo)致的非高斯分布影響較小，只有當(dāng)β為非整數(shù)階時(shí)，v的改變才使得峰現(xiàn)時(shí)間出現(xiàn)小幅變化，這應(yīng)與分?jǐn)?shù)階階數(shù)性質(zhì)相關(guān)，而與流速的改變無(wú)關(guān)。減小α和β都會(huì)使得穿透曲線的峰值及尾部相對(duì)變化率增加，并且隨著v的增大，階數(shù)不同帶來(lái)的差異性逐漸增大。這種現(xiàn)象可以被解釋為α的時(shí)間記憶性和β的空間非局域性使得方程對(duì)于流速參數(shù)v更加敏感，使得穿透曲線的拖尾現(xiàn)象更加明顯，且α和β越小這種特性的影響越強(qiáng)。

分?jǐn)?shù)階導(dǎo)數(shù)對(duì)流速改變的影響見(jiàn)圖2(c)、2(d)。當(dāng)α、β確定時(shí)，增大彌散系數(shù)D主要使得穿透曲線的峰現(xiàn)時(shí)間提前，而對(duì)峰值幾乎沒(méi)有影響。這是由于彌散系數(shù)為空間相關(guān)參數(shù)，主要影響穿透曲線的非高斯分布，而對(duì)冪律拖尾現(xiàn)象幾乎沒(méi)有影響，只有當(dāng)α為分?jǐn)?shù)階時(shí)，D的改變才使得峰值出現(xiàn)小幅變化，這應(yīng)與分?jǐn)?shù)階階數(shù)性質(zhì)相關(guān)，而與彌散系數(shù)的改變無(wú)關(guān)。減小α、β使得穿透曲線的峰現(xiàn)時(shí)間變化率增大，并且隨著D的增大，階數(shù)不同帶來(lái)的差異性逐漸增大。這種現(xiàn)象可以被解釋為分?jǐn)?shù)階的引入使得方程對(duì)參數(shù)D更加敏感，β的空間非局域性使得同一時(shí)間空間上的所有點(diǎn)都對(duì)觀測(cè)點(diǎn)的溶質(zhì)濃度產(chǎn)生影響，因而使得溶質(zhì)擴(kuò)散加快，β越小這種影響更加明顯，故減小β值會(huì)削弱彌散系數(shù)改變帶來(lái)的差異。

3 案例分析

3.1 野外試驗(yàn)情況

潛流帶場(chǎng)地位于新安江水文試驗(yàn)站的試驗(yàn)場(chǎng)(117°47′E，29°43′N)，區(qū)域內(nèi)地形切割強(qiáng)烈，河道表層介質(zhì)以礫石、砂石為主，松散巖類孔隙水和紅層孔隙裂隙水是流域內(nèi)最為主要的地下水，介質(zhì)非均質(zhì)性強(qiáng)、貯水性差，地表-地下水交換作用強(qiáng)烈。選擇河道近岸處5 m×7 m區(qū)域?yàn)樵囼?yàn)區(qū)，上游設(shè)置一注射井，觀測(cè)井下游1.5 m設(shè)置間隔1 m的3×5觀測(cè)井群(編號(hào)A1、A2、A3、B1、B2、…、E1、E2、E3)，觀測(cè)井及注射井均打入距離河道表面0.5 m深處，注射井底部0.1 m寬度側(cè)壁上均勻分布細(xì)孔，每個(gè)觀測(cè)點(diǎn)位設(shè)置3個(gè)觀測(cè)井，分別在地下距離河床表面0.1、0.3、0.5 m處設(shè)置一細(xì)孔(分別編號(hào)1、2、3)，3個(gè)縱深下的介質(zhì)從上至下由大塊碎石逐漸轉(zhuǎn)為砂礫，試驗(yàn)期間河道溫度為23.5～24.3 ℃。試驗(yàn)場(chǎng)地與點(diǎn)位分布見(jiàn)圖3。

圖3 試驗(yàn)場(chǎng)地與點(diǎn)位分布Fig.3 Schematic diagram of test site and the locations of injection and observation points

在試驗(yàn)區(qū)內(nèi)進(jìn)行示蹤試驗(yàn)。向注射井中注入質(zhì)量濃度為250 mg/L的紅色溴離子溶液7.5 L，從溶質(zhì)注入開(kāi)始間隔5 min抽取觀測(cè)井中河水并監(jiān)測(cè)溴離子質(zhì)量濃度。試驗(yàn)期內(nèi)，測(cè)得試驗(yàn)區(qū)域內(nèi)的平均水深為0.3 m，地表水流速為0.4～1.7 m/s，試驗(yàn)開(kāi)始5 min后在觀測(cè)井A2中監(jiān)測(cè)到溴離子質(zhì)量濃度發(fā)生明顯變化，同時(shí)觀察到紅色溶液溢出觀測(cè)井外壁，證明試驗(yàn)區(qū)內(nèi)發(fā)生強(qiáng)烈的潛流交換。

3.2 模擬結(jié)果分析

根據(jù)監(jiān)測(cè)所得溶質(zhì)數(shù)據(jù)，首先分別運(yùn)用經(jīng)典一維、二維整數(shù)階對(duì)流彌散方程解析解進(jìn)行參數(shù)反演。在此基礎(chǔ)上，運(yùn)用分?jǐn)?shù)階方程對(duì)溶質(zhì)數(shù)據(jù)進(jìn)行擬合，計(jì)算采用的流速v、彌散系數(shù)D、阻滯因子Rd與對(duì)應(yīng)解析解反演參數(shù)相同，得到時(shí)間、空間分?jǐn)?shù)階階數(shù)α、β值見(jiàn)表1。

表1 模型參數(shù)Tab.1 Model Parameters

值得注意的是：所有溶質(zhì)數(shù)據(jù)均來(lái)自A2點(diǎn)位，A2-1、A2-2、A2-3這3個(gè)觀測(cè)點(diǎn)位在水平面內(nèi)位置一致，而縱向深度有所差異。在進(jìn)行一維模型設(shè)計(jì)時(shí)，溶質(zhì)的遷移路徑被簡(jiǎn)化為由注射井注射孔出發(fā)到觀測(cè)井A2不同深度3個(gè)觀測(cè)點(diǎn)A2-1、A2-2、A2-3的一維直線。觀測(cè)井到不同點(diǎn)位的直線距離均近似于觀測(cè)井與注射井間的水平直線距離(1.5 m)，而潛流帶介質(zhì)在縱深方向上的強(qiáng)非均質(zhì)性給溶質(zhì)遷移提供了快速通道，使得溶質(zhì)運(yùn)動(dòng)的變異性增強(qiáng)。因此，溶質(zhì)穿透曲線表現(xiàn)出的不同形態(tài)的原因?yàn)槿苜|(zhì)在強(qiáng)非均質(zhì)性潛流帶介質(zhì)中遷移的路徑差異，表現(xiàn)在模型參數(shù)設(shè)計(jì)上即為遷移路徑的流速v、介質(zhì)的彌散系數(shù)D及阻滯因子Rd存在差異。針對(duì)二維模型，溶質(zhì)在介質(zhì)中傳輸路徑的差異已通過(guò)設(shè)置縱深不同體現(xiàn)，故3個(gè)觀測(cè)點(diǎn)的參數(shù)保持一致。

從模擬結(jié)果來(lái)看，3個(gè)測(cè)點(diǎn)的時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)α均取1，這是由于試驗(yàn)場(chǎng)地河床的介質(zhì)粗大，潛流流速大，溶質(zhì)遷移伴隨的拖尾現(xiàn)象并不顯著，溶質(zhì)遷移的滯時(shí)效應(yīng)不明顯。對(duì)比之下，3個(gè)測(cè)點(diǎn)的物理參數(shù)及空間分?jǐn)?shù)階導(dǎo)數(shù)存在較大差異，這與試驗(yàn)場(chǎng)地介質(zhì)強(qiáng)烈的非均質(zhì)性有關(guān)。

結(jié)合試驗(yàn)現(xiàn)場(chǎng)情況：測(cè)點(diǎn)A2-1周圍介質(zhì)以大塊卵石、碎石為主；測(cè)點(diǎn)A2-3周圍介質(zhì)以礫石、砂石為主，雖分屬不同介質(zhì)帶，但整體位于淺層埋深且介質(zhì)的非均質(zhì)性相較測(cè)點(diǎn)A2-2不高；而測(cè)點(diǎn)A2-2位于卵石與礫石地交界帶，介質(zhì)的非均質(zhì)性極強(qiáng)，從微觀角度看，溶質(zhì)在其中的遷移路徑更加復(fù)雜，見(jiàn)圖4。溶質(zhì)在孔隙度大的上層(A2-1)遷移時(shí)通過(guò)速度快，其中：中間層(A2-2)遷移時(shí)部分溶質(zhì)優(yōu)先經(jīng)過(guò)快速通道到達(dá)觀測(cè)點(diǎn)，下層(A2-3)介質(zhì)層致密導(dǎo)致溶質(zhì)通過(guò)速度較上層、中層低；彌散系數(shù)反映溶質(zhì)在介質(zhì)中的彌散作用強(qiáng)度，介質(zhì)固體骨架導(dǎo)致的流速分布不均是產(chǎn)生彌散的根本原因，在均質(zhì)、大孔隙介質(zhì)中的彌散作用強(qiáng)，而中間層礫石顆粒填充卵石孔隙的非均質(zhì)介質(zhì)大大降低了溶質(zhì)的彌散作用，因此中間層的彌散系數(shù)小于上下兩層；Rd的大小是介質(zhì)層吸附能力、水力條件、溶質(zhì)質(zhì)量濃度共同作用的結(jié)果，彌散和擴(kuò)散作用相對(duì)對(duì)流作用影響越大時(shí)Rd越小[32]，隨著介質(zhì)顆粒粒徑增大有效孔隙度增大、總孔隙度減小，中間層介質(zhì)非均質(zhì)性強(qiáng)、粒徑變異性大，有效孔隙度和總孔隙度小于上下兩層，對(duì)流作用減弱，因此Rd小。一維模型物理參數(shù)的反演結(jié)果得到驗(yàn)證且具有物理意義，同時(shí)進(jìn)一步證明了在潛流帶中介質(zhì)的強(qiáng)非均質(zhì)性極大地影響了溶質(zhì)運(yùn)移的特征。

介質(zhì)的強(qiáng)非均質(zhì)性同時(shí)會(huì)影響空間分?jǐn)?shù)階導(dǎo)數(shù)β的取值：上下層的空間分?jǐn)?shù)階導(dǎo)數(shù)明顯大于中間層。這是由于空間分?jǐn)?shù)階β表示介質(zhì)空間的非均質(zhì)性，β越小，介質(zhì)的非均質(zhì)性越強(qiáng)、對(duì)溶質(zhì)運(yùn)移的阻滯作用也越強(qiáng)。下層介質(zhì)以砂礫為主，顆粒介質(zhì)細(xì)密、總孔隙度大，對(duì)應(yīng)空間分?jǐn)?shù)階導(dǎo)數(shù)β最大；上層介質(zhì)為大塊卵石，介質(zhì)孔隙度大且對(duì)溶質(zhì)的阻滯效果小；中間層屬于混合介質(zhì)層，砂石顆粒填充卵石間孔隙，形態(tài)不一的孔隙形成多種滲透通道，介質(zhì)的非均質(zhì)性極強(qiáng)，因此對(duì)應(yīng)2的空間分?jǐn)?shù)階導(dǎo)數(shù)明顯小于其余兩點(diǎn)，故也可證明參數(shù)設(shè)置的合理性。

圖4 溶質(zhì)在不同介質(zhì)中的遷移路徑Fig.4 solute migration paths in different media

整數(shù)階和分?jǐn)?shù)階解的擬合曲線見(jiàn)圖5。在兩種維度的模擬中，分?jǐn)?shù)階的模擬效果在描述峰現(xiàn)時(shí)間、峰值質(zhì)量濃度、穿透曲線形態(tài)、質(zhì)量濃度下降階段溶質(zhì)的拖尾現(xiàn)象中均具有更好的效果，證明時(shí)間、空間分?jǐn)?shù)階導(dǎo)數(shù)的引入可以更加精確地刻畫非均質(zhì)介質(zhì)中的溶質(zhì)擴(kuò)散過(guò)程，對(duì)于溶質(zhì)消散階段發(fā)生的拖尾現(xiàn)象描述更加理想。同時(shí)注意到，在計(jì)算初始上升階段，兩種方法的模擬結(jié)果均不理想，尤其是在初始兩個(gè)觀測(cè)時(shí)刻，A2-1點(diǎn)位的分?jǐn)?shù)階方法模擬效果反而不如傳統(tǒng)解析解。這是由于初期濃度上升速度快，數(shù)值解振蕩計(jì)算的誤差被放大,而在后期溶質(zhì)擴(kuò)散過(guò)程趨緩、計(jì)算穩(wěn)定后，數(shù)值解的計(jì)算誤差減小，分?jǐn)?shù)階數(shù)值解對(duì)于穿透曲線尾部描述的優(yōu)越性得以體現(xiàn)。

圖5 解析解、分?jǐn)?shù)階數(shù)值解穿透曲線對(duì)比Fig.5 Comparison of breakthrough curves of analytical solution and fractional numerical solution

與一維模擬結(jié)果相比，二維方程模擬得到的穿透曲線的擬合程度并不理想。A2-2測(cè)點(diǎn)的峰現(xiàn)時(shí)間明顯比實(shí)測(cè)值提前；A2-3測(cè)點(diǎn)穿透曲線后段模擬的溶質(zhì)質(zhì)量濃度也明顯比實(shí)測(cè)值要高；相較而言，A2-1測(cè)點(diǎn)的模擬結(jié)果仍比較理想。盡管二維分?jǐn)?shù)階模擬結(jié)果較解析解而言在上述方面有了進(jìn)步，但二維的模擬結(jié)果仍不如一維，具體原因?qū)⒃?.2節(jié)進(jìn)一步討論分析。

4 討 論

4.1 分?jǐn)?shù)階與整數(shù)階模擬結(jié)果對(duì)比

進(jìn)一步量化兩種計(jì)算方法對(duì)于實(shí)測(cè)數(shù)據(jù)的擬合效果。引入均方根誤差(RMSE)和擬合優(yōu)度(R2)兩種評(píng)價(jià)指標(biāo)判斷模擬值同觀測(cè)值之間的偏差以及模型的優(yōu)劣，計(jì)算結(jié)果見(jiàn)表2。

表2 解析解與分?jǐn)?shù)階數(shù)值解模擬結(jié)果評(píng)價(jià)Tab.2 Evaluation of simulation results of analytical solutions and fractional numerical solutions

從兩種評(píng)價(jià)指標(biāo)來(lái)看，分?jǐn)?shù)階數(shù)值解法計(jì)算得到的結(jié)果均方根誤差更小、擬合優(yōu)度更高；從穿透曲線的擬合效果來(lái)看，分?jǐn)?shù)階數(shù)值解法模擬的穿透曲線形態(tài)更加貼近實(shí)測(cè)數(shù)據(jù)，尤其是在峰現(xiàn)時(shí)間和溶質(zhì)質(zhì)量濃度下降階段的表現(xiàn)更加突出。這是由于潛流帶介質(zhì)非均質(zhì)性強(qiáng)，不同尺度的孔隙在介質(zhì)中形成復(fù)雜的強(qiáng)滲透通道，溶質(zhì)的傳輸路徑多樣，通常使得穿透曲線峰現(xiàn)時(shí)間推后、尾部產(chǎn)生拖尾現(xiàn)象。傳統(tǒng)的對(duì)流-彌散方程對(duì)于介質(zhì)的描述單一，適用于模擬均質(zhì)介質(zhì)中的溶質(zhì)遷移過(guò)程，難以很好地描述潛流帶的溶質(zhì)運(yùn)移特性,而時(shí)間分?jǐn)?shù)階α的引入體現(xiàn)了溶質(zhì)運(yùn)移過(guò)程中的滯時(shí)效應(yīng)，空間分?jǐn)?shù)階β的引入更加詳細(xì)地描述介質(zhì)的空間非均質(zhì)性。在本試驗(yàn)中，由于礫石河床質(zhì)高流速、強(qiáng)非均質(zhì)性的特點(diǎn)，空間分?jǐn)?shù)階導(dǎo)數(shù)β起到主要作用，而時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)α的滯時(shí)效應(yīng)并不顯著，但在低流速的潛流帶溶質(zhì)運(yùn)移過(guò)程中，α的影響將凸顯。因此，α和β兩者共同作用對(duì)于描述潛流帶中的溶質(zhì)運(yùn)移現(xiàn)象具有標(biāo)志性意義。

4.2 一維與二維模擬差異

從模擬結(jié)果(圖5和表2)看，二維模擬的效果整體不及一維好，但二維分?jǐn)?shù)階的模擬效果依然優(yōu)于二維解析解。這種現(xiàn)象可以由以下原因解釋：研究區(qū)潛流帶介質(zhì)顆粒呈現(xiàn)上大下小分布，復(fù)雜的介質(zhì)顆粒排布在水流及溶質(zhì)的遷移路徑中形成快速通道，這導(dǎo)致部分溶質(zhì)從注射井出發(fā)通過(guò)快速通道更快地到達(dá)觀測(cè)井；而另一部分溶質(zhì)則在小孔隙通道中更加緩慢地遷移，強(qiáng)烈的介質(zhì)非均質(zhì)性離散了溶質(zhì)的遷移特征，增強(qiáng)粒子運(yùn)動(dòng)的隨機(jī)性，延長(zhǎng)了尾部溶質(zhì)的滯留時(shí)間。

從分?jǐn)?shù)階參數(shù)設(shè)置來(lái)看：采用一維模型描述溶質(zhì)運(yùn)移時(shí)，可以通過(guò)設(shè)置不同的空間分?jǐn)?shù)階導(dǎo)數(shù)β描述溶質(zhì)粒子在不同運(yùn)動(dòng)軌跡下的離散化特征，而二維模型中的空間分?jǐn)?shù)階導(dǎo)數(shù)取值實(shí)際上是二維空間內(nèi)溶質(zhì)粒子運(yùn)動(dòng)的最大共性解；從物理參數(shù)的設(shè)置來(lái)看，由于研究區(qū)潛流帶存在多維水流，使用一維對(duì)流-彌散方程進(jìn)行模擬時(shí)設(shè)置了不同的流速v和彌散系數(shù)D以模擬溶質(zhì)到達(dá)不同縱深點(diǎn)位所經(jīng)過(guò)的路徑存在差異，而二維對(duì)流-彌散方程在計(jì)算時(shí)設(shè)置統(tǒng)一的水流、介質(zhì)參數(shù)，僅通過(guò)y方向的坐標(biāo)體現(xiàn)3個(gè)測(cè)點(diǎn)間差異，難以對(duì)潛流帶介質(zhì)的復(fù)雜特性進(jìn)行準(zhǔn)確描述。此外二維模擬中，靠近河床表面的A2-1和A2-3測(cè)點(diǎn)的模擬結(jié)果優(yōu)于A2-2測(cè)點(diǎn)。結(jié)合現(xiàn)場(chǎng)試驗(yàn)的情況(圖4、圖5)，A2-1測(cè)點(diǎn)所在埋深處介質(zhì)仍以大塊卵石、碎石居多，水流流速和介質(zhì)孔隙度大，對(duì)于溶質(zhì)運(yùn)輸?shù)淖铚饔眯。珹2-3測(cè)點(diǎn)所處埋深介質(zhì)以礫石、砂石為主，由介質(zhì)非均質(zhì)性造成的影響小；而A2-2測(cè)點(diǎn)所處埋深介質(zhì)由大塊卵石轉(zhuǎn)為砂石，介質(zhì)非均質(zhì)性增強(qiáng)，水流形態(tài)復(fù)雜，介質(zhì)強(qiáng)烈的非均質(zhì)性對(duì)溶質(zhì)運(yùn)移作用，使得二維的模擬產(chǎn)生較大偏差，同時(shí)也證明了描述介質(zhì)的非均質(zhì)性對(duì)精確模擬潛流帶溶質(zhì)運(yùn)移過(guò)程的重要性。

二維模型和一維模型本質(zhì)上存在理論差異，因此各適應(yīng)的模擬條件有所不同。其中：一維模型的優(yōu)勢(shì)在于可以通過(guò)參數(shù)設(shè)置描述縱深方向介質(zhì)非均質(zhì)性導(dǎo)致的溶質(zhì)粒子遷移路徑差異，對(duì)不同深度的監(jiān)測(cè)點(diǎn)位的溶質(zhì)質(zhì)量濃度變化模擬有良好表現(xiàn)；二維模型的優(yōu)勢(shì)在于模擬強(qiáng)非均質(zhì)性介質(zhì)條件下溶質(zhì)受到滯時(shí)效應(yīng)和超擴(kuò)散現(xiàn)象的影響，在水平面內(nèi)到達(dá)不同點(diǎn)的遷移特征差異。盡管二維分?jǐn)?shù)階模擬結(jié)果在本試驗(yàn)中并不理想，但可以預(yù)見(jiàn)，若能在等深水平面上觀測(cè)到更多點(diǎn)位的溶質(zhì)數(shù)據(jù)(如A1、A3、B1…)，將二維分?jǐn)?shù)階的坐標(biāo)系放置于水平面，在沒(méi)有縱深方向介質(zhì)差異導(dǎo)致的水流、介質(zhì)參數(shù)差異前提下，二維分?jǐn)?shù)階模型將具有更加適用的場(chǎng)景。此外，若想在縱深方向上提高二維分?jǐn)?shù)階模擬的準(zhǔn)確性，可以考慮將流速、介質(zhì)參數(shù)在縱深方向進(jìn)行分層設(shè)置以提高模擬的準(zhǔn)確性。

5 結(jié) 論

從分?jǐn)?shù)階方程的物理意義上看：α體現(xiàn)溶質(zhì)運(yùn)移過(guò)程的歷史記憶性，α減小，穿透曲線峰現(xiàn)時(shí)間滯后、溶質(zhì)質(zhì)量濃度下降速度減緩、曲線趨于扁平且拖尾性增強(qiáng)；β體現(xiàn)擴(kuò)散行為的空間非局域性，β減小，穿透曲線峰現(xiàn)時(shí)間提前、溶質(zhì)擴(kuò)散速度加快。

分?jǐn)?shù)階的引入使得方程對(duì)流速和彌散系數(shù)的敏感性增強(qiáng)，對(duì)流速的敏感性主要體現(xiàn)在峰值即尾部的濃度值變化，對(duì)彌散系數(shù)的敏感性主要體現(xiàn)在峰現(xiàn)時(shí)間的改變，同時(shí)發(fā)現(xiàn)，α和β越小，這種對(duì)于時(shí)間、空間參數(shù)的敏感性越高。

結(jié)合野外示蹤試驗(yàn)對(duì)比整數(shù)階和分?jǐn)?shù)階分別在一維和二維的模擬結(jié)果證明：分?jǐn)?shù)階方法對(duì)于潛流帶溶質(zhì)運(yùn)移的模擬效果較傳統(tǒng)方法更好，尤其是對(duì)于穿透曲線中溶質(zhì)的峰現(xiàn)時(shí)間及拖尾現(xiàn)象的描述準(zhǔn)確；二維方程的模擬結(jié)果由于在參數(shù)設(shè)置上將各測(cè)點(diǎn)概化為相同的參數(shù)，難以體現(xiàn)潛流帶非均質(zhì)性強(qiáng)的特點(diǎn)，故使用一維方程的模擬結(jié)果更優(yōu)。因此，一維分?jǐn)?shù)階方法最可以彌補(bǔ)傳統(tǒng)對(duì)流-擴(kuò)散方程無(wú)法描述潛流帶介質(zhì)非均質(zhì)性程度高的缺陷，二維分?jǐn)?shù)階方法在無(wú)須考慮縱深介質(zhì)均質(zhì)性導(dǎo)致的流速、介質(zhì)參數(shù)設(shè)置差異的水平面溶質(zhì)運(yùn)移模擬中也具有積極意義。