胶体颗粒在潜流带中沉积过程及机制

2021-09-27张中天金光球陈鹤翔唐洪武

水科学进展 2021年5期

张中天,金光球,陈鹤翔,唐洪武,李凌

(1.河海大学水文水资源与水利工程科学国家重点实验室,江苏南京 210098；2.河海大学水利水电学院,江苏南京 210098；3.西湖大学,浙江杭州 310000)

潜流带是河流系统中地表水和地下水双向交换、相互作用的区域[1-2],同时含有地表水和地下水是潜流带的一个重要特征[3]。Triska等[4-5]对潜流带进行了量化,结果表明潜流带中地表水含量大于10%但小于98%。潜流带是重要的生物化学反应场所,物质和能量交换频繁[6],也是许多水生生物的栖息地和重要鱼类的产卵场[7-9],具有调节水温、涵养水源、促进氮、磷等物质循环和保护生物多样性等众多功能,其健康尤为重要[10-12]。然而,河流中的胶体颗粒会在潜流交换、颗粒沉降与河床截留等因素共同作用下沉积于潜流带[13]。受到水流、河床表面压差等因素作用,潜流带中沉积的胶体或形成一定的分布[13-14],对区域水动力特征和生态环境产生不同形式的影响[15]。研究胶体颗粒在潜流带中的沉积分布规律及其响应特征,对潜流生态修复和河流健康管理都有积极意义。

胶体颗粒在潜流带中的沉积会对区域水生态环境产生重要的影响。研究表明胶体颗粒在潜流带中迁移会被截留在多孔性河床中,改变区域河床孔隙率和渗透系数[14,16],从而减缓区域水流作用,为潜流带沉积物中的生物膜提供较为稳定的生境[17-18]。另外潜流带中的胶体颗粒也可成为污染物和生源物质的良好载体[19-21],对于颗粒粒径较小的胶体颗粒,较大表面积和特殊的表面物理化学性质使之对于重金离子、有机化合物、放射性物质等污染物质具有很强的吸附能力[22-23],并促进这些污染物在河床中富集[24-25]。一些胶体颗粒会因为水流和沉降作用,携带污染物进入河床更深的区域[26-27]。

有关颗粒物在河流系统中的迁移和沉积特征,前人已做过相关的研究,Drummond等[28-30]对河流中细颗粒物的输运进行了现场调查,发现大量颗粒被截留在河床中,甚至是低密度的微塑料颗粒。张维等[31]通过现场实验,发现水动力扰动作用可促进胶体颗粒在多孔介质中的释放与迁移。Zheng和Chen[32]及Zhang和Selim[33]发现地下水中的胶体颗粒迁移会受到高盐度的阻碍。Karwan和Saiers[34]对悬浮颗粒的潜流交换过程进行了数值模拟,但没有考虑胶体颗粒的沉降作用。Packman等[35]通过室内物理实验,研究了对流和沉降作用下的胶体颗粒在潜流带的地表-地下水界面通量规律,却没有探讨其在河床中的分布情况。Jin等[14,36]利用纳米颗粒为材料,通过室内实验探究了纳米颗粒的在河床中释放作用,并利用数值模拟手段,探讨了细颗粒在潜流带截留产生的负反馈机制。多孔介质中胶体颗粒的迁移能力强,能极大地促进其他污染物的迁移[37]。Karathanasis[38]通过野外调查发现在胶体颗粒存在的情况下铜和锌在多孔介质中的迁移速度是没有胶体存在情况下迁移速度的5～50倍,而铅的迁移速度更是达到10～300倍。Kan和Tomson[39]发现溶解的有机物质或者胶体颗粒可以促进DDT的迁移。O′Carroll等[40]发现,砷在潜流带的迁移深度和距离会因为铁氧化物胶体的存在而增加。然而这些研究多集中于细颗粒对污染物的作用以及细颗粒的地表-地下水界面通量上,尽管Jin等[13]对胶体颗粒在河床中的截留做出了初步的探究,但是仍有许多问题尚未详尽,如胶体颗粒在潜流带中沉积的范围及其分布呈什么特征,胶体颗粒在潜流带沉积分布特征对不同因素的响应机制是什么。

本文通过室内水槽实验和数值模拟,研究胶体颗粒在潜流带的截留面积及水平方向截留分布的变化规律,并对渗透系数、上覆水流速、沉降速度、截留系数、沙波高度对沉积分布的影响做量化分析。以探明潜流交换、颗粒沉降与河床截留作用对胶体沉积分布的作用机制,为胶体在潜流带中的生态环境效应研究提供科学依据。

1 循环水槽实验

1.1 实验用沙与胶体颗粒

实验用沙为石英砂。实验前要利用酸洗和碱洗等方法,去除沙中杂质以及吸附在沙粒表面的金属、有机质等。将沙洗净后,测量其各项物理指标:粒径范围为0.25～0.6 mm,中值粒径d50=0.30 mm,孔隙率θ=0.33,渗透系数K=8.84×10-4m/s[41-42]。

实验选用高岭土胶体。为了研究胶体沉积的物理特性,防止在迁移实验中出现絮凝,本文利用NaCl对胶体进行预处理,方法详见文献[43]。胶体浓度通过紫外分光光度计来测量。事先标定了胶体颗粒吸光度和其浓度之间的关系,实验得知所测胶体颗粒浓度范围内(0.01～0.4 g/L)的吸光度和浓度有很好的相关性[13],整个实验温度不变(19℃)。

1.2 循环水槽实验操作方法

本文采用的水槽是可变坡循环水槽(图1(a))。参照Elliott和Brooks[44]水槽设计总长11.5 m,宽30 cm,高50 cm。水槽壁面使用的是钢化玻璃,方便实验观察。水槽进口段设置自制能量均匀扩散装置,以减小或消除在进水口部分水流的翻滚和波动,从而使进入河床的水流更加稳定。水槽的底板和进出口部分用的是耐腐蚀性的PVC板。循环回水管道用的是PVC管、可视化的电磁流量计(ORBLDBE-50)、抗腐蚀性水泵(IRG-50-100)。这些材料具有较强的抗酸碱和盐腐蚀性,同时对金属离子没有吸附性。在水槽的尾部安装了MPtroll(9500),能够自动监测上覆水的pH和盐度。

本次实验的有效沙床段长为10 m,厚度为14～16 cm。实验前,先在水槽中注满去离子水,然后铺设三角形沙床(图1(b)),以保证沙床处于全饱和的状态。之后进行水槽调试,在水箱中加入NaCl溶液,并运行水槽,通过测量上覆水和空隙水中NaCl溶液的质量浓度,判断河流上覆水和河床之间交换是否达到平衡,最终水槽中(包括孔隙水中)NaCl质量浓度为115.5 mg/L,这个过程一般用时12～15 h。

图1 循环水槽实验示意Fig.1 Diagram of circulating flume experiment

正式实验时,首先用NaCl处理过的胶体在去离子水中充分搅拌,使其均匀,在循环水1个周期内加入到水箱中。实验过程中,按预定的时间抽取上覆水,并测量其胶体颗粒质量浓度,每次质量浓度测量需要在12 h之内完成,以防止胶体沉淀。由于实验过程有蒸发等水量损失,需要每4 h向水箱中加入400～450 mL的去离子水。整个实验持续68 h,直到上覆水胶体颗粒的质量浓度达到基本稳定时停止实验。实验中所用的各项参数见表1。

表1 实验及数值模拟中参数取值Table 1 Values of parameters in experiment and simulation

水槽实验结束后,则需测量胶体颗粒在河床中的最终截留分布情况。此时对沙波进行分层分块提取(图1(b),图1(c))。沙波的迎水一侧水平方向平均分成4列,每块的宽度为2.5 cm,背水一侧水平方向平均分成2列,每列的宽度为3.0 cm。为了方便取样,本实验特制了高1 cm、宽4 cm的取沙器使得每次取样厚度为1 cm。将以上沙样加入试管,并加入15 mL去离子水。为了使胶体脱离沙粒进入溶液中,利用超声波振荡仪进行振荡,振荡时间为25 min[13]。最后测量上层溶液的吸光度,并换算成截留量。

2 数值模拟方法

本文利用胶体颗粒运移理论的控制方程,进行数值模拟。具体操作方法是利用数学模型将上覆水流动场、河床多孔介质中孔隙水流场和胶体颗粒运移场三方面进行耦合(表2)。

表2 3种物理场关系表Table 2 Relationships among three physical fields

2.1 上覆水流场

本文利用FLUENT 软件提供的有限体积法[46]对河流上覆水进行数值模拟,模型采用的是二维紊流模型。模拟中用的控制方程是雷诺时均化的纳维-斯托克斯方程(图2,Cc为孔隙水中胶体颗粒的质量浓度)。流体考虑为不可压缩的各向同性的连续体。由于k-ω模型在模拟有反向压力梯度的分离流和沙波引起的漩涡时效果最好,故本文采用Wilcox的k-ω模型[45,47]来进行紊流模拟。

图2 边界条件示意Fig.2 Diagram of boundary conditions

河床的表面是自由表面无摩擦的边界,整个上覆水是基于刚盖假定[24],自由边界是0剪切应力状态(图2(a))。这种假定在小沙波情况下被广泛使用[36,45]。在水平尺度的求解域中只有1个沙波,因为认为沙波在空间上沿水流方向是周期性变化的。所以把周期性变化的压力和流速边界条件,加在求解域的两侧,包括河床[36,41]。水沙界面被视为无滑移壁边界,参数选取见表1。

2.2 孔隙水流场

本文利用COMSOL Multiphysics进行对沙波中的孔隙水流场的数值模拟,该软件是一个基于有限单元法进行数值计算的软件[48]。假设沙床为不可压缩的多孔介质,使用达西定律和不可压缩流体连续性方程作为控制方程,进行孔隙水流场的模拟。其方程如式(1)和式(2)所示。

(1)

(2)

式中:xi为位移,m,i=1表示水平方向,向右为正,i=2表示竖直方向,向上为正；ui为xi方向上的流速,m/s；P为孔隙水压力,Pa。

选取COMSOL Multiphysics软件中的达西定律计算模块。本文针对孔隙水流场数值模拟方法,考虑到θ和K不随时间、空间或其他环境条件的变化而变化,可被当作常量。对于孔隙水流场边界条件的设置,利用2.1节中上覆水流场计算结果,导出下边界的压力分布作为孔隙水流场计算的上边界条件。对于沙波的两侧,设置在空间上沿水流方向周期性变化的压力边界,在沙波的顶部,水沙交界面设置为压力边界条件。沙波的底端设置为无流动边界条件(图2(a))。孔隙水计算中的参数选取见表1,孔隙水流场见图3(a),胶体颗粒运移场和孔隙水流场和沉降场的叠加见图3(b)。

图3 孔隙水流场和胶体颗粒迁移速度Fig.3 Flow fields of porewater and velocity of colloidal particles

2.3 胶体颗粒运移场

本文考虑的胶体颗粒运移场作用区域在沙波中。传统的胶粒在多孔介质中运动方程考虑了对流、扩散和截留的影响,考虑到胶体的沉降作用,本文创新性地在竖直方向的流速上加上了沉降速度,如式(3):

(3)

(4)

式中:t为时间,s；Cc为孔隙水中胶体颗粒的质量浓度,kg/m3；S为单位体积胶体颗粒截留量,kg/m3；vs为胶体颗粒沉降速度,m/s,取向上为正方向,横向为0；Dij为水动力弥散系数张量,m2/s,i=1,2,j=1,2，计算公式如下:

(5)

(6)

式中:|u|为速度向量的大小,m/s。

对于胶体颗粒运移场的边界条件,水沙分界面设置为物质通量边界条件(图2(b)),表达形式如下:

(7)

式中:C为上覆水胶体颗粒质量浓度,kg/m3；n为垂直于界面的单位向量,方向指向沙床；u为速度向量。在t时刻COMSOL Multiphysics软件模拟中上覆水浓度C(t)的计算式如下式8所示。

(8)

式中:Vs为除沙床孔隙水外的水体总体积,m3；Lf为水槽有效长度,m；B为水槽宽度,m。参数取值见表1。

2.4 数值模拟工况设计与沉积特征选择

胶体颗粒在潜流带中迁移时,由于水流和自身沉降作用[13,35],使之向河床中和下游迁移(图4(a))。在不考虑胶体颗粒的情况下,孔隙水的流场如图3(a)所示,沙波中会形成1个地表-地下水循环区(潜流带),以及1个不与地表水有直接交换的区域地下水[41]。考虑到胶体颗粒的沉降作用,胶体颗粒运移场将是孔隙水的流场和沉降场的叠加(图3(b))。由于河床的截留作用,部分胶体颗粒会截留于河床中,最终在河床中形成一定的截留分布。

为了讨论截留和沉降对胶体颗粒沉积的影响,本文通过数值模拟参数调整,设计了4种工况:考虑胶体颗粒的截留和沉降(工况1)、不考虑截留和沉降(工况2)、只考虑截留不考虑沉降(工况3)和只考虑沉降不考虑截留(工况4)。若考虑截留和沉降,则将截留系数和沉降速度设置为表1中的数值；若不考虑截留和沉降则将截留系数和沉降速度设置为0。

图4 胶体颗粒在潜流带迁移截留示意Fig.4 Diagram of colloidal particles transport and retention in hyporheic zone

在考虑截留面积和水平方向截留分布对某种因素(X)的响应关系时,本文则选取了渗透系数、上覆水流速、沉降速度、截留系数及沙波高度进行敏感性分析,讨论其改变为原始值不同倍数(X/X0)时对截留面积和水平方向截留分布的影响。

3 结果分析

3.1 胶体颗粒质量分布

整个实验过程中,胶体颗粒可能存在于上覆水、孔隙水和截留在沙床这3种状态中。本文通过水槽实验,监测上覆水中胶体颗粒的质量浓度,则可以换算为上覆水中胶体颗粒质量的占比(mi/m0,i=1,2,3,分别代表上覆水、孔隙水和截留3种状态的胶体质量；m0表示这3种状态的总质量),见图5(a)。对比图5(a)中上覆水中胶体颗粒质量占比的实验值和模拟值可以看出,数值模拟很好地反演了上覆水中胶体颗粒质量的变化。上覆水中胶体颗粒的质量一直处于下降趋势。在前10 h,上覆水中胶体颗粒的质量下降很快；在10～50 h,上覆水中胶体颗粒的质量下降的趋势变缓；50 h之后,上覆水中胶体颗粒的质量趋于平稳,且接近于0。由于在实验前期(10 h之前),上覆水中的胶体颗粒质量占比较大(0.6～1.0),且沙床中的胶体颗粒质量浓度较小(0～0.4),胶体颗粒能较快速地通过对流、扩散和沉降等作用,进入沙波。到实验后期(40 h之后),上覆水中的胶体颗粒质量占比变小(小于0.2),沙波中的胶体颗粒质量占比增大(大于0.8),此时胶体颗粒进入沙床的过程变缓。

图5 胶体颗粒3种状态中的相对质量变化Fig.5 Relative mass variation in three phases of colloidal particles

对于孔隙水中的胶体颗粒质量占比而言,其在0～1 h的时间快速增大到峰值,之后又缓慢减少,直至趋近于0。孔隙水中的胶体颗粒总量在整个胶体颗粒迁移、截留过程中,占比都非常小,即使是峰值占比,也没达到总量的1%。这是由于孔隙水中的胶体颗粒是上覆水胶体颗粒和截留胶体颗粒的一个中间过渡态,从上覆水进入孔隙水的胶体颗粒能迅速被沙床截留。对于截留态的胶体颗粒质量占比而言,其趋势几乎正好和上覆水中胶体颗粒质量占比变化趋势相反,其占比先快速增加,后趋于平缓,并接近于1。

为探求胶体颗粒从上覆水进入沙波的动力机制,本文通过数值模拟扩展,讨论了4种工况下上覆水中胶体质量浓度的变化,如图5(b)所示。对比工况1和工况2可知,胶体颗粒的截留作用和沉降作用能导致上覆水质量浓度趋于0。对比工况1和工况3可知,截留作用可以使得上覆水中胶体颗粒质量浓度趋于0,且沉降速度的存在加速了这一过程。对比工况1和工况4可知,截留对于上覆水中胶体质量变化影响较小。综上所述,胶体颗粒从上覆水进入河床的动力过程主要受到截留和沉降作用影响,且沉降作用的效果更为明显。

3.2 胶体颗粒在沙床中的截留分布和截留面积

水槽实验结束后,本文进行了分块取样(图1(b)、图1(c)),并测试了每一分块中的胶体颗粒截留量(图6(a),A0表示单个沙波总面积),同时利用数值模拟得出的胶体颗粒截留分布如图6(b)所示。可见实验和数值模拟呈现出比较相似的胶体颗粒截留分布,浅层的胶体颗粒截留量大,深层的较少；迎水面胶体颗粒截留量较大,背水面较少。实验和数值模拟呈现出相似的结果,截留前锋线的位置基本一致。

图6(c)是通过数值模拟得出的截留前锋线随时间的变化。总体来说,截留面积在x2>-0.05 m区域,这是由于河床截留作用阻碍了胶体颗粒向更深处迁移。截留面积分布在实验前期(0～9 h)独立成块,其中背水面有1个区域内没有截留,这是由于背水面胶体颗粒进入通量较小(图3(b)),背水面进入沙床的胶体颗粒还不够使得2个区域相连接。截留面积分布在实验中后期(以约9 h作为临界点)连成一体,且截留前锋线缓慢向下推进。通过图6(c)得出截留面积变化趋势为在前期(0～0.5 h)快速增加,之后其增加趋势变缓。

图6 胶体颗粒截留分布和截留面积Fig.6 Distribution of retention and retention area for colloidal particles

3.3 胶体颗粒水平方向截留分布特征

图7中红色点划线和蓝色实线分别表示实验结束时实测和模拟的胶体颗粒水平方向截留分布。总体规律呈现出迎水面(x1<0.1 m)胶体颗粒截留量较大,背水面(x1>0.1 m)较少。这是由于迎水面沙波一般受到的水压要大于背水面,胶体颗粒在迎水面向下运动的速度高于背水面的速度(图3(b)),胶体颗粒更容易从迎水面进入并截留于此。

图7 胶体颗粒水平方向截留分布Fig.7 Distribution of horizontal retention for colloidal particles

图7同时显示随着时间推移胶体颗粒如何逐步截留累计的数值模拟结果。可以看出,无论何时,迎水面胶体颗粒截留量大于背水面。胶体颗粒截留量在前期明显快于后期。1～5 h和20～50 h的截留增量相当,而最后的50～68 h的截留量几乎不变。这也印证了图5中上覆水中和截留胶体颗粒质量的变化趋势。

图7中,胶体颗粒水平方向截留分布在迎水面x1=0.05 m处附近出现1个特殊的小谷值,这是由于胶体颗粒运移场在此处呈辐射状向下输运(图3(b)),此处进入的胶体颗粒除了有往下的趋势,还有往左右扩散的趋势,因此胶体颗粒在x1=0.05 m截留量反而变小。

4 讨论

4.1 胶体颗粒截留面积的敏感性分析

为了探究胶体颗粒截留面积对各种因素的响应关系,以实际渗透系数、上覆水流速、胶体颗粒沉降速度、沙波截留系数和沙波高度为初始值(K0、V0、Vs0、Kr0、Hb0),探讨实际截留面积(Ai0)为初始值的敏感性。

从图8中可以看出,K、U、vs和Ai呈正相关。这是由于K、U的增大,增加了孔隙水流速(式(2)),而孔隙水将胶体颗粒带入了更深的区域；vs速度向下,较大的沉降速度能直接将胶体颗粒带入更深的区域。vs对Ai的影响在0～0.5倍区间比较剧烈,在0.5～5倍区间比较小。这是由于沉降速度较小时,背水面以出流为主(图3(a)),使得上覆水中的胶体颗粒难以从背水面进入和截留[13]。随着沉降速度增大,胶体颗粒从背水面进入,同时会出现胶体截留(图9(d))。较大的截留系数能将更多的胶体颗粒截留在浅层,沉积截留更集中于浅层,kr和Ai呈负相关(图8)。沙波高度同样对截留面积有所影响,从上图8中可以看出随着沙波高度上升,截留面积有所下降。这是由于在其他条件不变时,较高的沙波高度导致了更多的胶体颗粒趋于水平运动,向下运动的趋势被弱化,故而截留面积减小。

图8 截留面积敏感性分析Fig.8 Sensitivity analysis of retention area

4.2 胶体颗粒水平方向截留分布的响应因子

图9 表示K、U、vs、kr、Hb变化对应水平方向(x1方向)截留量的分布情况。可以得出,对于所有的情况,迎水面(x1<0.1 m)胶体颗粒截留量较大,背水面(x1>0.1 m)较少。对于图9(a)和图9(b)而言,随着K和U增大,胶体颗粒水平方向截留分布的峰谷值差距也越大。根据公式2可知,孔隙水流速是随K、U的增大而增加的；再根据图3,孔隙水流速的增大,导致迎水面流速增大,进而胶体颗粒进入通量增加,背水面孔隙水出流流速同样增大,故背水面的胶体颗粒进入通量减小。因此,在截留量相同的情况下,K、U增大,扩大了其截留分布的峰谷值的差距。

对于图9(c)而言,随着kr增大,胶体颗粒水平方向截留的峰谷值差距也越大。这是由于在迎水面,较大的截留系数能导致胶体颗粒更迅速地截留于此；而在背水面,尽管截留系数较大,但是进入的胶体颗粒通量较小,所以其截留量较小。

对于图9(d)而言,随着vs增大,水平方向截留峰谷值差距也越小。这是由于沉降速度越大,胶体颗粒运移场将由沉降速度主导(图3),而本文假设的沉降速度在时空上是个定值。故较大的沉降速度将导致水平方向截留量趋于平稳。值得注意的是,当沉降速度为0和0.2倍时,沙波背水面只有出流速度(图3(a)),导致背水面部分区域没有胶体颗粒进入,截留量为0。

对于图9(e)而言,随着Hb增大,水平方向截留峰谷值差距也越大。这是由于沙波高度增大会导致迎水面和水流夹角更大,随水流运动的胶体颗粒直接冲击迎水面概率更大,胶体更容易进入迎水面,并且更难在背水面沉积,故水平方向截留峰谷值差距增大。对于沙波高度小的工况,水平方向截留分布将趋于均匀。