考虑变流速及尺度效应的一维溶质运移模型及半解析解

李艳芳

(鄂尔多斯市国能神东监理有限责任公司,内蒙古 鄂尔多斯 719315)

0 引言

地下水溶质运移理论一直受到水文地质学者的广泛关注,主要原因是溶质在多孔介质运移过程中涉及复杂的物理、化学和生物过程[1-3]。许多学者基于对流-弥散方程(Advection-Dispersion Equation,ADE)构建了大量的解析和数值的溶质运移理论模型,预测污染物在含水层的运移特征。例如,在均质含水层中建立了考虑常系数的一维、二维、三维ADE解析模型[4-5]。然而,传统的常系数ADE模型在室内和野外应用过程中受到了广泛的质疑,主要是由于受到变流速、尺度效应及非均质等因素的影响,无法有效地解译溶质的穿透曲线,尤其是很难刻画非费克溶质运移现象(如拖尾、提前穿透及多峰等现象)[6-7]。

国内外学者开展溶质运移解析模型研究时,为了避免地下水流速变化导致ADE模型难以求解,通常假定地下水流是稳定的,即为常数。而实际含水层由于受到补、径、排模式与地表水体水位及地下水开采量变化等因素的影响,含水层中地下水流速会随时间发生变化。例如,河流和湖泊水位季节性变化或洪水会导致沿岸附近地下水位随着地表水的变化而变化[8]。此外,物理、化学和生物堵塞也会导致室内土柱试验地下水流速不断衰减,进而导致实验的参数反演出现偏差。例如,Zaheer[9]等针对低渗透介质进行了一系列土柱溶质运移实验发现,地下水流速随时间呈指数下降。许多学者针对变流速溶质运移问题开展了变流速条件下的溶质运移理论研究。Jaiswal[10]等提出了线性、渐近及指数等形式来刻画随时间变化的流速,并利用时间-空间变换方法获得变流速条件下的一维ADE解析解。提出的线性及指数函数的流速将分别趋近于无穷和零,而实际地下水流速应当趋近于某个有限值。为此,Li[11]等提出了指数变化的地下水流速渐进方程,并利用积分变换获得变流速淋滤作用下的一维溶质运移的解析解。

除了地下水流速随着时间变化外,弥散度与迁移尺度有着密切联系[12-14]。Gelhar[15]等通过分析59个不同现场的观测数据得到,含水层的弥散度具有尺度效应。Wang[6]等进行土柱示踪试验,研究表明,尺度效应在溶质运移过程中发挥了重要的作用。大量野外与室内实验分析认为,多孔介质的非均性是导致弥散度尺度效应的主要原因。目前,许多学者综合多方面的实验数据,总结了迁移距离和弥散度之间的经验关系,包括线性、指数及抛物线等多种形式[16-17],构建了相应的尺度效应理论模型。例如,You和Zhan[18]利用Laplace变化推导了考虑线性渐近和指数变化的一维半解析解,拟合了Huang[19]等土柱试验数据,拟合效果优于传统的ADE模型。

图1 一维溶质运移示意图Fig.1 Schematic diagram of one-dimensional solute transport

在溶质运移理论中,变流速和弥散度的尺度效应是不容忽视的。通过文献研究发现,同时考虑变流速与尺度效应的研究较少,为此构建了考虑变流速及尺度效应的一维溶质运移模型,利用指数变化的地下水流速方程来刻画多孔介质中地下水流速的变化,并耦合弥散度随着距离指数变化的函数。利用积分变换及Laplace变换方法获得变流速及尺度效应溶质运移的半解析解,探究指数变化的地下水流速及弥散度对溶质运移的影响机理。

1 数学模型构建

关于多孔介质中溶质运移模型中的地下水流速通常设定为常数,但实际过程中地下水流速是随时间变化的,时常表现出指数增加或指数衰减趋势[11]。例如,当含水层受到生物、化学及物理堵塞时会导致渗透性减小,使得地下水流速表现出指数衰减[20-21]。河岸带及海岸带的松散含水层由于地表水位的波动(如洪水、潮汐等作用)会引起含水层地下水流速呈指数增或减的趋势[8]。因此考虑地下水流速可以表示为指数增或指数减的形式[11]:

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

(1)

式中:v0为多孔介质地下水初始的流速[L/T];v1为最终稳定的地下水流速[L/T];λ为流速变化指数[1/T]。当v0>v1时,式(1)表示地下水流速为指数减的情况;当v0

(2)

为了概化为溶质运移距离与弥散度之间的函数关系,许多学者提出多种经验方程,目前常见的经验方程主要有线性、指数、抛物线及渐进方程[16]。室内实验及野外试验研究表明,指数方程更符合实际情况[18],因此选取弥散度随溶质运移距离为指数变化方程,可表示为:

α(x)=α0(1-e-bx/L)

(3)

D(x,t)=α(x)v(t)

(4)

式中:α0(x)为随溶质运移距离为指数变化的弥散度[L];α0为迁移距离足够大时渐进的弥散度[L];b/L为弥散度变化指数[1/L];L为多孔介质的空间距离[L]。

在建立流速及尺度效应的一维溶质运移模型之前,为了简化数学模型,假定:①多孔介质地下水流动为非稳定流,流速为指数变化。②多孔介质弥散度迁移随迁移距离的变化而变化,弥散度为指数变化。③溶质为惰性溶质,仅考虑吸附作用。考虑变流速及尺度效应的对流-弥散方程(ADE),可表示为:

(5)

式中:C为溶质的浓度[M/L3];x为距离[L];t为时间[T];D(x,t)为随时间和空间变化的弥散系数[L2/T],且v(t)为随时间变化的地下水流速[L/T];R是迟滞系数[M/L3]。将式(3)、(4)带入式(5)可得:

(6)

初始及边界条件为:

C(t=0,x)=0

(7)

(8)

(9)

式中:C0为x=0处的给定浓度[M/L3];t0为脉冲注入的时间[T]。将该模型简称为EE模型(Advection-dispersion equation with exponentially time-dependent flow velocity and distance-dependent dispersivity)。

2 数学模型求解

为了化简上述的数学模型,式(6)两边同除以τ(t),可得:

(10)

为了获取模型的解,引入一个新的积分变换:

(11)

则式(6)可表示为:

(12)

边界条件及初始条件在相同的积分变化下可表示为:

C(T=0,x)=0

(13)

(14)

(15)

(16)

C(TD=0,xD)=0

(17)

(18)

(19)

对式(16)～(19)分别对TD作Laplace变换,得:

(20)

(21)

(22)

定义新的变量z=e-bxD,则式(20)可以表示为:

(23)

式(23)可表示为超几何函数的形式:

(24)

式中:Q= 0

(25)

(26)

(27)

求解式(25)～(26),可得:

(28)

(29)

式(24)的通解为:

(30)

若T0

(31)

X[m(e-bm)F(m+1,m+1;m-n+1;e-b)]+Y[n(e-bn)F(n+1,n+1;n-m+1;e-b)]=0

(32)

根据函数F(m,m+1;m-n+1; 1)的性质,当m-n+1-(m+m+1)>0(Gao et al., 2010),可得:

(33)

(34)

式中:Γ(·) 为Gama函数。解式(31)～(32),可获得系数X及Y(表1):

表1 式(30)的系数表达式Tab.1 Coefficient expression of Equation (30)

(35)

(36)

若T0>TD0,利用相同的计算方法可得式(30)的系数X及Y为:

(37)

(38)

上述推导的EE在Laplace空间的解析解包含Gama函数等特殊函数,很难利用解析逆变换的方法进行求解。因此采用数值逆变换的方法来获得实空间下的径向溶质迁移的解,利用Stehfest数值逆变换方法即可得到EE在实空间的解。

3 结果与讨论

根据以上模型的推导和计算,利用脉冲注入边界的解析模型来计算讨论表皮区域的弥散度、孔隙度及厚度对径向溶质迁移的影响。对于以下分析,默认的模型参数见表2。

表2 模型中参数取值Tab.2 Parameter value in the model

3.1 不同模型对比分析

为了研究变流速与尺度效应耦合作用对多孔介质中溶质运移的影响,分别对比分析了ADE、Li[11]等及You和Zhan[18]的解(图2)。Li等的模型考虑地下水流速是指数衰减的,弥散度是定值。You和Zhan的模型考虑指数变化的弥撒度,地下水流速是定值。ADE模型中地下水流速v及弥散度α为定值,分别设置为1 m/day、2 m,同时Li等与You和Zhan及EE模型相同参数设置是相同的,见表2。由图2可知,ADE的穿透曲线在早期溶质浓度高于其他模型,而后期的浓度低于其他模型。这是由于ADE模型中的地下水流速和弥散度为常数,而Li等的模型地下水流速为指数衰减,逐渐衰减到v1=0.5 m/day,地下水的流速小于ADE,而You和Zhan的模型弥散度随着迁移距离逐渐增大至α=2 m,其弥散度也小于ADE,故ADE相比其他模型表现出提前穿透现象。此外,EE解的穿透曲线表现出明显的滞后现象,主要是由于该模型的地下水流速指数减小趋近于1 m/day,且弥散度指数增加趋近于2 m,故其穿透曲线早期浓度最低,后期浓度高于其他模型。说明地下水流速的衰减与尺度效应耦合作用会对溶质运移产生较大影响。

图2 不同模型的溶质穿透曲线Fig.2 Solute penetration curves for different models

通过对比不同模型结果发现,EE和Li等的解相比于ADE表现出明显的拖尾现象,说明地下水流速的减小会导致穿透曲线表现出拖尾现象。而拖尾现象(非费克运移)通常被认为是由于多孔介质的非均质性导致的,忽视地下水流速的衰减也会导致非费克运移现象,这为非费克运移研究提供了新思路。尤其是低渗介质溶质运移时常表现出拖尾现象,其中地下水流速由于受到物理、生物、化学堵塞及固结作用,导致地下水流速衰减,这可能是引起非费克运移的影响因素之一。

3.2 尺度效应对溶质运移的影响

图3为变流速条件下不同b/L值时在x=10 m的溶质穿透曲线。弥散度变化指数b/L分别为0.1、0.2、0.4,模型其他参数设置为L=50 m,v0=1 m/day,v1=0.5 m/day,α0=2 m 且λ=0.1。将EE模型的计算结果与Li等的结果进行对比分析发现,Li等的解计算溶质的穿透曲线早期浓度高于EE解,说明尺度效应会减缓溶质迁移。这是由于Li等解的弥散度是定值α=α0=2 m,而EE解的弥散度是随着迁移距离逐渐增大到α0= 2 m。因此Li等解的弥散度大于EE解。此外,随着b/L值的增大,EE解的穿透曲线早期浓度增大,后期浓度减小,且峰值不断减小。表明变流速条件下弥散度尺度效应会进一步导致产生穿透曲线削峰及拖尾现象。

图3 变流速条件下不同取值的弥散度变化指数b/L 在x=10 m处所对应的溶质穿透曲线Fig.3 Dispersion change index b/L with different values corresponding to the solute penetration curve at x=10 m under the condition of variable flow velocity

3.3 变流速对溶质运移的影响

图4为弥散度尺度效应下不同λ值时在x=10 m的溶质穿透曲线。地下水流速衰减指数λ分为0.1、0.2、0.5,模型其他参数设置为L=50 m, b/L=0.1/m,v0=1 m/day,v1=0.5 m/day,α0=2 m。将EE解的计算结果与You和Zhan的解进行对比分析发现,You和Zhan解的计算的穿透曲线表现出提前穿透现象,说明地下水流速的衰减会减缓溶质迁移。这是由于You和Zhan解的地下水流速是定值v=v0=1 m/day,而EE解的地下水流速是随着迁移时间逐渐由v0=1 m/day衰减到v1=0.5 m/day。因此You和Zhan解的地下水流速大于EE解。随着λ值的增大,EE解的穿透曲线早期浓度减小,后期浓度增大,溶质迁移的速率不断减小,表明弥散度尺度效应下地下水流速的变化对溶质运移迁移规律产生了较大的影响,其影响结果不容忽视。

图4 弥散度尺度效应下不同取值的地下水流速衰减指数λ在x=10 m处所对应的溶质穿透曲线Fig.4 Attenuation index λ of groundwater velocity with different values corresponding to the solute penetration curve at x=10 m under the diffusivity scale effect

图5为尺度效应下不同稳定流速v1在x=10 m的溶质穿透曲线。根据式(1)可知,不同的v1值表示地下水流速指数降低或增加,若v1>v0表示指数增加,若v1

图5 尺度效应下不同稳定流速v1在x=10 m的溶质穿透曲线Fig.5 Solute penetration curve at x=10 m at different stable velocity v1 under scale effect

4 结论

基于对流弥散方程,构建了考虑随时间和空间指数变化的地下水流速及弥散度的一维溶质运移模型,通过积分变换及Laplace变换获得了模型的半解析解,分析了指数变化的流速及弥散度对溶质运移的影响,得到的结论如下:

变流速条件下,随着弥散度增长速率的增大,溶质穿透曲线的拖尾现象明显且峰值逐渐减小。因此多孔介质溶质运过程中的尺度效应因素不能忽视。

在尺度效应影响下,多孔介质中指数变化的地下水流速对溶质运移有较大影响,流速的衰减指数越大,地下水流速衰减越快,穿透曲线早期浓度减小,后期浓度增大,导致穿透曲线拖尾现象的产生。

稳定地下水流速在溶质运移过程中起到了关键作用,其大小反映了地下水流速增加或减小的幅度,变化幅度越大,溶质运移快慢表现得越明显。稳定速度的增加会增大对流作用,导致穿透曲线峰值高,拖尾现象减弱。