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

魯程鵬，林雨竹，張勇，秦巍，吳成城，劉波，束龍倉

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

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

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

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

1 時空分數(shù)階模型的建立及數(shù)值解法

1.1 模型建立

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

(1)

式中：C(x,y,t)為溶質(zhì)濃度；V和D分別為介質(zhì)中的平均流速和彌散系數(shù)；Rd表示阻滯因子；f(x,y,t)表示源匯項。α(0<α?1)為時間分數(shù)階導(dǎo)數(shù)，使用Caputo式求解：當(dāng)α=1時，方程化為空間分數(shù)階對流-彌散方程；β(1<β?2)為空間分數(shù)階導(dǎo)數(shù)，使用Riemman-Liouville式求解：當(dāng)β=2時，方程化為時間分數(shù)階對流-彌散方程；當(dāng)α=1,β=2時，方程簡化為傳統(tǒng)對流-彌散方程。設(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]開發(fā)的已被證明穩(wěn)定收斂的一維時空分數(shù)階數(shù)值解法，并采用穩(wěn)定的差分法將方程擴展至二維條件[30-31]，實現(xiàn)二維時空分數(shù)階數(shù)值解的運算。

分數(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ù)階方程的數(shù)值解為

(9)

2 分數(shù)階階數(shù)對溶質(zhì)運移過程的影響

2.1 分數(shù)階階數(shù)物理意義

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

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

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

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

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

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

2.2 分數(shù)階對溶質(zhì)運移物理參數(shù)的影響

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

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

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

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

3 案例分析

3.1 野外試驗情況

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

圖3 試驗場地與點位分布Fig.3 Schematic diagram of test site and the locations of injection and observation points

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

3.2 模擬結(jié)果分析

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

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

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

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

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

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

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

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

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

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

4 討 論

4.1 分數(shù)階與整數(shù)階模擬結(jié)果對比

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

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

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

4.2 一維與二維模擬差異

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

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

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

5 結(jié) 論

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

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

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