基于时空分数阶导数模型的潜流带溶质运移模拟

鲁程鹏，林雨竹，张勇，秦巍，吴成城，刘波，束龙仓

(1.河海大学水文水资源学院，江苏 南京 210098；2.Department of Geological Sciences of University of Alabama,Tuscaloosa USA AL 35487)

潜流带是河道内地表地下水相互作用的饱水沉积层，连接河道及地下水体，三者的生物、物理、化学特性彼此影响[1]。潜流带中发生的物质交换常被认为是影响河道污染物迁移及能量交换的基本过程。河道中的溶质粒子常随水流进入潜流带停留一段时间，最终在下游某些点返回河流。潜流带中的溶质运移很好地解释了河道地表及地下水体间的物质联系，因此探究潜流带中的溶质迁移规律对于进一步了解河道内反应性运输和生物吸收十分重要[2]。

潜流带沉积物在空间上不均匀分布以及河道的物理化学性质、生物地球化学过程、温度变化、地表-地下水位波动都会对河床介质的非均质性产生影响，渗透系数的空间差异性是这种非均质性的重要表征：Tang等[3]设置了不同介质条件模拟非均质介质中饱和流的水流通量交换情况，发现渗透系数的差异会对模拟结果影响较大；Wang等[4-6]发现天然河道中水平面内顺水流方向的河床渗透系数没有明显变化，而纵深方向介质具有较强非均质性；Ren等[7]发现渗透系数对孔隙度变化十分敏感，具体表现在大孔隙是过流的主要通道，但较大的孔隙通道在压实过程中总是被提前关闭；Leung 等[8]发现河床温度可以影响水的黏度进而改变介质渗透系数；Zhu等[9]研究发现微生物空间中生长产生阻塞作用会降低局部介质的渗透性；活性矿物质或离子的氧化和凝胶状物质的形成同样也会通过影响河床电位梯度控制潜流交换的速率[10]。当前针对模拟潜流带溶质过程的模型已经十分丰富： Singh等[11]提出了一种降阶模型用于模拟流量变化和河床地形相互作用影响下的潜流交换特征；Adu[12]运用解析法构造了一种非点源混合单元模型，模拟得到的穿透曲线与现场示踪剂测试结果拟合程度较好；Sherman等[13]构建了随机游走-拉格朗日模型提高了溶质粒子在介质中的迁移运动模拟的精确度；对流-弥散方程作为地下水溶质运移的基础理论方程也依然被广泛应用于示踪模拟。但现有方法往往忽略了河床介质的非均质性，增大了对污染物迁移模拟的不准确性[5]，近几年出现的分数阶方法在数学层面上解释了由介质非均质性引起的溶质粒子随机运动问题，将影响粒子迁移的因素转化为时间、空间分数阶导数综合判断，具有良好的模拟效果。

基于Fick定律得到的经典的对流-弥散方程已被广泛运用于地下水溶质运移模拟，但针对潜流带这一结构特殊、介质非均质性强且水流方向复杂的多孔介质，溶质易呈现拖尾及非高斯分布的反常扩散特征，传统对流-弥散方程模型很难准确描述这类非局域性特征[14]。分数阶模型将偏微分方程中求导的整数阶数变为分数，已被证明在描述物质的反常扩散方面具有物理意义明确、可描述介质非均质性等优势[14-15]。当前针对分数阶方法的研究多集中于算法的开发(数值稳定、快速求解)[16-20]，而在分数阶的应用领域的探索主要集中于模拟饱和热传导问题[21-22]、非饱和土壤水分扩散[23-24]、土柱沙箱试验[25-26]、泥沙运动等方面[27-28]。本文将时空分数阶模型首次引入潜流带溶质运移问题，与传统整数阶模型结果进行对比，试图探究分数阶方法对应实际问题的物理意义以及对于水流介质参数的灵敏性，最后结合野外现场示踪试验数据，分别在一维和二维两个维度探索分数阶方法在潜流带溶质运移方面的应用效果与前景。

1 时空分数阶模型的建立及数值解法

1.1 模型建立

在传统对流-弥散方程基础上分别引入时间、空间分数阶导数，建立时空分数阶二维对流-弥散方程

(1)

式中：C(x,y,t)为溶质浓度；V和D分别为介质中的平均流速和弥散系数；Rd表示阻滞因子；f(x,y,t)表示源汇项。α(0<α≪1)为时间分数阶导数，使用Caputo式求解：当α=1时，方程化为空间分数阶对流-弥散方程；β(1<β≪2)为空间分数阶导数，使用Riemman-Liouville式求解：当β=2时，方程化为时间分数阶对流-弥散方程；当α=1,β=2时，方程简化为传统对流-弥散方程。设置初始条件和边界条件分别为

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 数值解法

采用Liu[29]开发的已被证明稳定收敛的一维时空分数阶数值解法，并采用稳定的差分法将方程扩展至二维条件[30-31]，实现二维时空分数阶数值解的运算。

分数阶导数分别离散为

(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。

得到时空分数阶方程的数值解为

(9)

2 分数阶阶数对溶质运移过程的影响

2.1 分数阶阶数物理意义

参考周璐莹等[31]对理想算例的讨论方法，设置理想条件(溶质浓度为均一化浓度)：一维空间(0，10 m)，计算时间100 min。初始时刻，在x=60 m瞬时释放100单位点源溶质，取弥散系数D=0.15 m2/min，对流流速v=0.01 m/min。

模型中固定β=2，分别设置α=0.2、0.4、0.6、0.8、1.0，计算x=6 m处的穿透曲线以及t=15 min时的溶质扩散曲线，研究α的变化对溶质运移的影响，结果见图1。由图1可以看出：在穿透曲线图1(a)中，随α减小，穿透曲线的峰现时间滞后，峰值变小，同时浓度的下降速度也减缓，穿透曲线的偏态性增强，扩散稳定后溶质浓度高，即α减小扩散拖尾现象越明显；在溶质扩散曲线图1(b)中，随α减小，溶质的峰向上游移动，扩散减缓，扩散曲线的偏态性增强，扩散趋于稳定后溶质的浓度变高，这也与图1(a)中反映的扩散曲线特征吻合。

图1 不同α和β在x=6 m处的归一化穿透曲线及t=15 min时刻的归一化溶质扩散曲线Fig.1 Breakthrough curves of different α and β at x=6 m and solute diffusion curves at t=15 min

结合扩散行为的物理意义及分数阶阶数的物理意义：时间分数阶α表征扩散行为的历史记忆性，即历史上某点所有时刻的溶质浓度均对当前浓度产生影响，也即扩散在非均质介质中的滞时效应。α值越小，溶质在介质中的滞时效应越大，扩散减缓，穿透曲线趋于扁平，扩散达到稳定时的溶质浓度也越大。历史记忆性使得反常扩散比正常的扩散过程更慢。

在前述计算条件下，固定α=1，分别设置β=1.2、1.4、1.6、1.8、2.0，计算x=6 m处的穿透曲线以及t=15 min时的溶质扩散曲线，研究β的变化对溶质运移的影响，结果如见图1(c)、1(d)。可观察到：在图1(c)中，随β增大，穿透曲线的峰现时间及趋于稳定的时间均延后，峰值变大，同时浓度的下降速度也增加，即β增大，溶质的扩散过程加快；在图1(d)中，随β减小，溶质的峰向上游移动，且峰值变大，即溶质的扩散过程减缓，扩散趋于稳定后溶质的浓度增大，这也与图1(c)中反映的扩散曲线特征吻合。

结合扩散行为的物理意义及分数阶阶数的物理意义：空间分数阶β表征扩散行为的空间非局域性，即同一时刻空间上的所有点(而非相邻点)都对某点的浓度产生影响，也即扩散在非均质介质中的尺度效应。这种非局域性使得扩散过程加快，β值越小，溶质扩散的非局域性越强，非均质介质中形成的快速溶质通道对于计算结果的影响越明显，从而导致扩散加快，峰值减小，峰现时间提前，溶质浓度越早趋于稳定，扩散达到稳定时的溶质浓度也越小。非局域性使得反常扩散比正常的扩散过程更快。

2.2 分数阶对溶质运移物理参数的影响

在前述理想计算条件下设置不同的流速v、弥散系数D，对比分数阶与整数阶方法对溶质运移过程中物理参数改变的影响。图2(a)、2(b)设置弥散系数D=0.15 m2/min，分别计算引入分数阶前后(α=1.0、0.6、0.2，β=2.0、1.6、1.2)改变流速(v=0、0.01、0.02、0.03 m/min)对穿透曲线形态的影响；图2(c)、2(d)设置流速v=0.01，分别计算引入分数阶前后(α=1.0、0.6、0.2，β=2.0、1.6、1.2)改变弥散系数(D=0.10、0.15、0.20、0.25 m2/min)对穿透曲线形态的影响。

图2 α和β对对流弥散方程系数(v、D)敏感性分析Fig.2 Comparison of sensitivity of different α and β to velocity v and Dispersion coefficient D

分数阶导数对弥散系数改变的影响见图2(a)、2(b)。当α、β确定时，增大流速v主要导致穿透曲线下降阶段斜率增加，而峰现时间几乎不变。这是由于流速为时间相关参数，主要影响穿透曲线的滞时效应(即拖尾现象)，而对空间差异导致的非高斯分布影响较小，只有当β为非整数阶时，v的改变才使得峰现时间出现小幅变化，这应与分数阶阶数性质相关，而与流速的改变无关。减小α和β都会使得穿透曲线的峰值及尾部相对变化率增加，并且随着v的增大，阶数不同带来的差异性逐渐增大。这种现象可以被解释为α的时间记忆性和β的空间非局域性使得方程对于流速参数v更加敏感，使得穿透曲线的拖尾现象更加明显，且α和β越小这种特性的影响越强。

分数阶导数对流速改变的影响见图2(c)、2(d)。当α、β确定时，增大弥散系数D主要使得穿透曲线的峰现时间提前，而对峰值几乎没有影响。这是由于弥散系数为空间相关参数，主要影响穿透曲线的非高斯分布，而对幂律拖尾现象几乎没有影响，只有当α为分数阶时，D的改变才使得峰值出现小幅变化，这应与分数阶阶数性质相关，而与弥散系数的改变无关。减小α、β使得穿透曲线的峰现时间变化率增大，并且随着D的增大，阶数不同带来的差异性逐渐增大。这种现象可以被解释为分数阶的引入使得方程对参数D更加敏感，β的空间非局域性使得同一时间空间上的所有点都对观测点的溶质浓度产生影响，因而使得溶质扩散加快，β越小这种影响更加明显，故减小β值会削弱弥散系数改变带来的差异。

3 案例分析

3.1 野外试验情况

潜流带场地位于新安江水文试验站的试验场(117°47′E，29°43′N)，区域内地形切割强烈，河道表层介质以砾石、砂石为主，松散岩类孔隙水和红层孔隙裂隙水是流域内最为主要的地下水，介质非均质性强、贮水性差，地表-地下水交换作用强烈。选择河道近岸处5 m×7 m区域为试验区，上游设置一注射井，观测井下游1.5 m设置间隔1 m的3×5观测井群(编号A1、A2、A3、B1、B2、…、E1、E2、E3)，观测井及注射井均打入距离河道表面0.5 m深处，注射井底部0.1 m宽度侧壁上均匀分布细孔，每个观测点位设置3个观测井，分别在地下距离河床表面0.1、0.3、0.5 m处设置一细孔(分别编号1、2、3)，3个纵深下的介质从上至下由大块碎石逐渐转为砂砾，试验期间河道温度为23.5～24.3 ℃。试验场地与点位分布见图3。

图3 试验场地与点位分布Fig.3 Schematic diagram of test site and the locations of injection and observation points

在试验区内进行示踪试验。向注射井中注入质量浓度为250 mg/L的红色溴离子溶液7.5 L，从溶质注入开始间隔5 min抽取观测井中河水并监测溴离子质量浓度。试验期内，测得试验区域内的平均水深为0.3 m，地表水流速为0.4～1.7 m/s，试验开始5 min后在观测井A2中监测到溴离子质量浓度发生明显变化，同时观察到红色溶液溢出观测井外壁，证明试验区内发生强烈的潜流交换。

3.2 模拟结果分析

根据监测所得溶质数据，首先分别运用经典一维、二维整数阶对流弥散方程解析解进行参数反演。在此基础上，运用分数阶方程对溶质数据进行拟合，计算采用的流速v、弥散系数D、阻滞因子Rd与对应解析解反演参数相同，得到时间、空间分数阶阶数α、β值见表1。

表1 模型参数Tab.1 Model Parameters

值得注意的是：所有溶质数据均来自A2点位，A2-1、A2-2、A2-3这3个观测点位在水平面内位置一致，而纵向深度有所差异。在进行一维模型设计时，溶质的迁移路径被简化为由注射井注射孔出发到观测井A2不同深度3个观测点A2-1、A2-2、A2-3的一维直线。观测井到不同点位的直线距离均近似于观测井与注射井间的水平直线距离(1.5 m)，而潜流带介质在纵深方向上的强非均质性给溶质迁移提供了快速通道，使得溶质运动的变异性增强。因此，溶质穿透曲线表现出的不同形态的原因为溶质在强非均质性潜流带介质中迁移的路径差异，表现在模型参数设计上即为迁移路径的流速v、介质的弥散系数D及阻滞因子Rd存在差异。针对二维模型，溶质在介质中传输路径的差异已通过设置纵深不同体现，故3个观测点的参数保持一致。

从模拟结果来看，3个测点的时间分数阶导数α均取1，这是由于试验场地河床的介质粗大，潜流流速大，溶质迁移伴随的拖尾现象并不显著，溶质迁移的滞时效应不明显。对比之下，3个测点的物理参数及空间分数阶导数存在较大差异，这与试验场地介质强烈的非均质性有关。

结合试验现场情况：测点A2-1周围介质以大块卵石、碎石为主；测点A2-3周围介质以砾石、砂石为主，虽分属不同介质带，但整体位于浅层埋深且介质的非均质性相较测点A2-2不高；而测点A2-2位于卵石与砾石地交界带，介质的非均质性极强，从微观角度看，溶质在其中的迁移路径更加复杂，见图4。溶质在孔隙度大的上层(A2-1)迁移时通过速度快，其中：中间层(A2-2)迁移时部分溶质优先经过快速通道到达观测点，下层(A2-3)介质层致密导致溶质通过速度较上层、中层低；弥散系数反映溶质在介质中的弥散作用强度，介质固体骨架导致的流速分布不均是产生弥散的根本原因，在均质、大孔隙介质中的弥散作用强，而中间层砾石颗粒填充卵石孔隙的非均质介质大大降低了溶质的弥散作用，因此中间层的弥散系数小于上下两层；Rd的大小是介质层吸附能力、水力条件、溶质质量浓度共同作用的结果，弥散和扩散作用相对对流作用影响越大时Rd越小[32]，随着介质颗粒粒径增大有效孔隙度增大、总孔隙度减小，中间层介质非均质性强、粒径变异性大，有效孔隙度和总孔隙度小于上下两层，对流作用减弱，因此Rd小。一维模型物理参数的反演结果得到验证且具有物理意义，同时进一步证明了在潜流带中介质的强非均质性极大地影响了溶质运移的特征。

介质的强非均质性同时会影响空间分数阶导数β的取值：上下层的空间分数阶导数明显大于中间层。这是由于空间分数阶β表示介质空间的非均质性，β越小，介质的非均质性越强、对溶质运移的阻滞作用也越强。下层介质以砂砾为主，颗粒介质细密、总孔隙度大，对应空间分数阶导数β最大；上层介质为大块卵石，介质孔隙度大且对溶质的阻滞效果小；中间层属于混合介质层，砂石颗粒填充卵石间孔隙，形态不一的孔隙形成多种渗透通道，介质的非均质性极强，因此对应2的空间分数阶导数明显小于其余两点，故也可证明参数设置的合理性。

图4 溶质在不同介质中的迁移路径Fig.4 solute migration paths in different media

整数阶和分数阶解的拟合曲线见图5。在两种维度的模拟中，分数阶的模拟效果在描述峰现时间、峰值质量浓度、穿透曲线形态、质量浓度下降阶段溶质的拖尾现象中均具有更好的效果，证明时间、空间分数阶导数的引入可以更加精确地刻画非均质介质中的溶质扩散过程，对于溶质消散阶段发生的拖尾现象描述更加理想。同时注意到，在计算初始上升阶段，两种方法的模拟结果均不理想，尤其是在初始两个观测时刻，A2-1点位的分数阶方法模拟效果反而不如传统解析解。这是由于初期浓度上升速度快，数值解振荡计算的误差被放大,而在后期溶质扩散过程趋缓、计算稳定后，数值解的计算误差减小，分数阶数值解对于穿透曲线尾部描述的优越性得以体现。

图5 解析解、分数阶数值解穿透曲线对比Fig.5 Comparison of breakthrough curves of analytical solution and fractional numerical solution

与一维模拟结果相比，二维方程模拟得到的穿透曲线的拟合程度并不理想。A2-2测点的峰现时间明显比实测值提前；A2-3测点穿透曲线后段模拟的溶质质量浓度也明显比实测值要高；相较而言，A2-1测点的模拟结果仍比较理想。尽管二维分数阶模拟结果较解析解而言在上述方面有了进步，但二维的模拟结果仍不如一维，具体原因将在4.2节进一步讨论分析。

4 讨 论

4.1 分数阶与整数阶模拟结果对比

进一步量化两种计算方法对于实测数据的拟合效果。引入均方根误差(RMSE)和拟合优度(R2)两种评价指标判断模拟值同观测值之间的偏差以及模型的优劣，计算结果见表2。

表2 解析解与分数阶数值解模拟结果评价Tab.2 Evaluation of simulation results of analytical solutions and fractional numerical solutions

从两种评价指标来看，分数阶数值解法计算得到的结果均方根误差更小、拟合优度更高；从穿透曲线的拟合效果来看，分数阶数值解法模拟的穿透曲线形态更加贴近实测数据，尤其是在峰现时间和溶质质量浓度下降阶段的表现更加突出。这是由于潜流带介质非均质性强，不同尺度的孔隙在介质中形成复杂的强渗透通道，溶质的传输路径多样，通常使得穿透曲线峰现时间推后、尾部产生拖尾现象。传统的对流-弥散方程对于介质的描述单一，适用于模拟均质介质中的溶质迁移过程，难以很好地描述潜流带的溶质运移特性,而时间分数阶α的引入体现了溶质运移过程中的滞时效应，空间分数阶β的引入更加详细地描述介质的空间非均质性。在本试验中，由于砾石河床质高流速、强非均质性的特点，空间分数阶导数β起到主要作用，而时间分数阶导数α的滞时效应并不显著，但在低流速的潜流带溶质运移过程中，α的影响将凸显。因此，α和β两者共同作用对于描述潜流带中的溶质运移现象具有标志性意义。

4.2 一维与二维模拟差异

从模拟结果(图5和表2)看，二维模拟的效果整体不及一维好，但二维分数阶的模拟效果依然优于二维解析解。这种现象可以由以下原因解释：研究区潜流带介质颗粒呈现上大下小分布，复杂的介质颗粒排布在水流及溶质的迁移路径中形成快速通道，这导致部分溶质从注射井出发通过快速通道更快地到达观测井；而另一部分溶质则在小孔隙通道中更加缓慢地迁移，强烈的介质非均质性离散了溶质的迁移特征，增强粒子运动的随机性，延长了尾部溶质的滞留时间。

从分数阶参数设置来看：采用一维模型描述溶质运移时，可以通过设置不同的空间分数阶导数β描述溶质粒子在不同运动轨迹下的离散化特征，而二维模型中的空间分数阶导数取值实际上是二维空间内溶质粒子运动的最大共性解；从物理参数的设置来看，由于研究区潜流带存在多维水流，使用一维对流-弥散方程进行模拟时设置了不同的流速v和弥散系数D以模拟溶质到达不同纵深点位所经过的路径存在差异，而二维对流-弥散方程在计算时设置统一的水流、介质参数，仅通过y方向的坐标体现3个测点间差异，难以对潜流带介质的复杂特性进行准确描述。此外二维模拟中，靠近河床表面的A2-1和A2-3测点的模拟结果优于A2-2测点。结合现场试验的情况(图4、图5)，A2-1测点所在埋深处介质仍以大块卵石、碎石居多，水流流速和介质孔隙度大，对于溶质运输的阻滞作用小，A2-3测点所处埋深介质以砾石、砂石为主，由介质非均质性造成的影响小；而A2-2测点所处埋深介质由大块卵石转为砂石，介质非均质性增强，水流形态复杂，介质强烈的非均质性对溶质运移作用，使得二维的模拟产生较大偏差，同时也证明了描述介质的非均质性对精确模拟潜流带溶质运移过程的重要性。

二维模型和一维模型本质上存在理论差异，因此各适应的模拟条件有所不同。其中：一维模型的优势在于可以通过参数设置描述纵深方向介质非均质性导致的溶质粒子迁移路径差异，对不同深度的监测点位的溶质质量浓度变化模拟有良好表现；二维模型的优势在于模拟强非均质性介质条件下溶质受到滞时效应和超扩散现象的影响，在水平面内到达不同点的迁移特征差异。尽管二维分数阶模拟结果在本试验中并不理想，但可以预见，若能在等深水平面上观测到更多点位的溶质数据(如A1、A3、B1…)，将二维分数阶的坐标系放置于水平面，在没有纵深方向介质差异导致的水流、介质参数差异前提下，二维分数阶模型将具有更加适用的场景。此外，若想在纵深方向上提高二维分数阶模拟的准确性，可以考虑将流速、介质参数在纵深方向进行分层设置以提高模拟的准确性。

5 结 论

从分数阶方程的物理意义上看：α体现溶质运移过程的历史记忆性，α减小，穿透曲线峰现时间滞后、溶质质量浓度下降速度减缓、曲线趋于扁平且拖尾性增强；β体现扩散行为的空间非局域性，β减小，穿透曲线峰现时间提前、溶质扩散速度加快。

分数阶的引入使得方程对流速和弥散系数的敏感性增强，对流速的敏感性主要体现在峰值即尾部的浓度值变化，对弥散系数的敏感性主要体现在峰现时间的改变，同时发现，α和β越小，这种对于时间、空间参数的敏感性越高。

结合野外示踪试验对比整数阶和分数阶分别在一维和二维的模拟结果证明：分数阶方法对于潜流带溶质运移的模拟效果较传统方法更好，尤其是对于穿透曲线中溶质的峰现时间及拖尾现象的描述准确；二维方程的模拟结果由于在参数设置上将各测点概化为相同的参数，难以体现潜流带非均质性强的特点，故使用一维方程的模拟结果更优。因此，一维分数阶方法最可以弥补传统对流-扩散方程无法描述潜流带介质非均质性程度高的缺陷，二维分数阶方法在无须考虑纵深介质均质性导致的流速、介质参数设置差异的水平面溶质运移模拟中也具有积极意义。