不同胶结度的多孔介质中溶质横向弥散的孔隙尺度模拟研究

Pore scale simulation study of transverse dispersion of solute in porous media with different cementation degrees

摘要

深入理解、精确刻画溶质孔隙尺度横向弥散行为是揭示溶质运移现象内在机制、提高溶质运移模拟精度的关键。基于孔隙尺度模拟方法研究了不同胶结度的多孔介质中溶质横向弥散特征。结果表明:当胶结度增大时,溶质横向弥散程度有所增强,但与纵向弥散相比增强幅度较小,且在模拟的胶结度范围内横向弥散未展现显著的尺度效应。研究还发现多孔介质中横向弥散系数D_T与纵向弥散系数D_L间的差异程度不是固定不变的,而是随着溶质运移距离逐渐增大最终达到稳态的过程,并且D_L/D_T的渐进值及达到渐进值所需时间均受胶结度的影响。因此,将D_T视为比D_L小一个数量级的值与溶质实际弥散行为存在较大偏差。同时,通过对比不同介质中流体流速横向、纵向分量概率密度及溶质对流速率概率密度的差异,揭示了以上溶质横向弥散特征变化及其与纵向弥散特征变化不同的内在机制。

关键词

横向弥散 / 多孔介质 / 胶结度 / 孔隙尺度 / 数值模拟

Key words

transverse dispersion / porous media / cementation degree / pore scale / numerical simulation

0 引言

深入研究、精确刻画饱水多孔介质中可溶性溶质运移在地下水污染防治、污染物地质封存、多相催化等领域至关重要^[1⇓-3]。而如何正确刻画溶质弥散是精确模拟饱水多孔介质中溶质运移的关键问题。空间中心二阶矩M₂及弥散系数D常被用来刻画溶质的弥散行为,其中空间中心二阶矩表示溶质羽在其质量中心位置周围的弥散程度。弥散系数则表示溶质的弥散速率,可由中心二阶矩对时间求导得到^[4]。根据溶质弥散方向与地下水主要流动方向的关系,又可将溶质弥散分为横向弥散和纵向弥散。垂直于水流主流方向上的溶质弥散称为横向弥散,沿水流主要流动方向上的溶质弥散为纵向弥散。

饱和多孔介质中的溶质运移多采用对流弥散模型(ADE)进行刻画。但是,近期研究表明饱水多孔介质中溶质弥散行为往往偏离对流弥散模型的预测结果。对流弥散模型应建立在溶质羽在其质量中心周围的弥散行为符合菲克定律的基础上。因此,溶质运移过程中弥散系数为常数,与运移时间和运移距离无关,溶质羽的空间中心一阶矩和二阶矩均随时间线性增大。但是,很多实验和模拟结果表明,溶质羽空间中心一阶矩、二阶矩常随时间非线性增大,溶质运移过程中弥散系数不是常数,而是随着运移时间的推移不断增大并最终达到一个渐近值^[5]。换而言之,弥散系数依赖于溶质运移距离,这种对流弥散模型无法精确刻画的溶质弥散现象被定义为尺度效应或非菲克运移^[6]。由于模拟溶质运移时弥散系数对模拟结果精度有重要作用,深入研究溶质弥散尺度效应有助于提高多孔介质中溶质运移模拟精度^[7]。

国内外学者针对溶质弥散开展了一系列卓有成效的研究工作。由于横向弥散远小于纵向弥散,早期研究多集中于溶质纵向弥散^[8]。对于横向弥散系数,常将其简单视为比纵向弥散系数小一个数量级的值^[9]。但是,随着研究的深入,学者们发现虽然横向弥散在强度上远小于纵向弥散,却对溶质的迁移有着不可忽略的影响^[10]。如横向弥散作用可中和溶质浓度波动,使得溶质羽浓度梯度趋于平滑^[11]。另外,横向弥散能够影响溶质羽的形状以及溶质羽与周围基质的接触面积,最终影响溶质的稀释作用强度及溶质与周围基质的混合程度,横向弥散的次线性行为是引起溶质非菲克弥散及不充分稀释的重要因素^[12-13]。近年来,随着对其重要性认识的深入,学者们在不同空间尺度上深入研究了溶质横向弥散,取得了较多进展^[14-15]。但是,目前关于介质孔隙尺度结构对溶质横向弥散影响的研究尚少。

含水介质的孔隙尺度结构特征对溶质弥散有显著影响^[16⇓-18]。而胶结度是自然界中土壤等含水多孔介质的基本属性之一,固体颗粒胶结程度的变化是自然界中最为常见的现象^[19]。已有研究表明,不同胶结度多孔介质的孔隙尺度结构具有显著差异,并对溶质纵向弥散有显著影响。但是,胶结度不同时,多孔介质孔隙尺度结构差异对溶质横向弥散有何影响尚不明确。近年来孔隙尺度模拟成为研究溶质弥散的重要手段^[20-21]。孔隙空间中的溶质迁移行为是宏观运移现象的本质。深入分析孔隙尺度溶质弥散是研究宏观运移现象机制、提出或改进运移模型的重要途径^[2,22]。本文针对以上问题,基于孔隙尺度模拟方法研究不同胶结度的多孔介质孔隙尺度结构差异对溶质横向弥散行为的影响及其机制。

本文的主要内容安排如下:在第一部分介绍主要的研究方法,包括不同胶结度的多孔介质构建,流场及溶质运移的模拟方法;第二部分内容主要针对模拟结果进行分析,包括不同胶结度的多孔介质中溶质横向弥散行为差异及其与纵向弥散行为间的关系;第三部分为总结得到的主要结论。

1 研究方法

1.1 孔隙介质构建

本文基于沉淀生成法构建不同胶结度的三维多孔介质,其中固体颗粒由球形粒子代表^[23]。介质构建过程主要包括三个步骤:(1)在模拟区域顶部随机的位置生成一个粒径为预先设定值的固体颗粒。(2)使生成的固体颗粒从其初始位置不断地沉降,直到有已经放置好的三个固体颗粒支撑或到达模拟区域的底部为止。通过该方法,当沉降过程结束时,沉降的固体颗粒到达其势能的最低点,从而保持多孔介质的拓扑结构。当只有三个或三个以上的固体颗粒支撑或达到底部时,固体颗粒才能达到稳定。循环执行以上两个步骤直到模拟区域内填满固体颗粒及粒间孔隙,生成所有固体颗粒相互不重叠即完全不胶结的多孔介质。(3)基于生成的完全不胶结的多孔介质,所有固体颗粒位置保持不变,粒径增加相同值,使固体颗粒轻微重叠,即固体颗粒开始轻微胶结,固体颗粒粒径持续增大,从而得到不同胶结度的多孔介质,如图1a所示。

本文生成的完全不胶结多孔介质中固体颗粒的粒径为2 mm,介质区域为12 mm×12 mm×40 mm(x×y×z),生成的介质如图1a所示,其孔隙率约为0.43。

基于该介质,增大固体粒径构建了三种固体颗粒不同程度胶结的多孔介质,并将介质按胶结度从小到大依次标记为PM43、PM35、PM30和PM25,分别代表孔隙率为0.43、0.35、0.30和0.25的多孔介质。统计固体颗粒胶结体积占固体总体积的比例P_c,并以此来定量衡量多孔介质胶结程度,结果如表1所示。从表中可以看出随着孔隙率从0.43降低到0.25,P_c从0.0%增大到了8.13%,表明从PM43到PM25,固体颗粒胶结度显著增大。

1.2 流场模拟

格子玻尔兹曼法(Lattice Boltzmann Method,LBM)能够处理复杂的边界条件,同时具有便于编程及并行计算等优点,已成为模拟多孔介质中地下水渗流的最常用方法之一^[24]。大量研究证明LBM是一种非常稳定、精确的流体流动模拟方法^[25]。因此,本文运用LBM中D3Q19模型(离散速度方向数为19个的三维模型)来模拟多孔介质内孔隙部分的流场^[26]。模拟中,将多孔介质的上边界以及下边界设置为定水头边界,水流方向为沿z轴正方向,平均流速为10 m/d。

1.3 溶质运移模拟

本文运用粒子追踪方法模拟溶质运移过程。该方法通过追踪大量代表溶质的粒子的迁移来模拟溶质运移过程。其中用粒子沿流线的漂移来模拟溶质对流作用,即在一个时间步长Δt内,溶质粒子以其所在位置处流体流速v沿流线方向的迁移距离为vΔt。用粒子在其位置附近的随机跳跃模拟溶质扩散作用,即溶质粒子向随机方向跳跃的随机距离为 \sqrt{2D_m\mathrm{\Delta }t}ξ。其中,D_m为溶质的扩散系数,ξ为符合标准正态分布的随机数。因此,溶质运移过程^[27]可概化为

(1)r(t+Δt)=r(t)+v(r(t))Δt+\sqrt{2D_m\Delta t}ξ

模拟过程中,将时间步长Δt设置为10^-4 s,以保障溶质粒子在每个时间步长内的迁移距离小于一个格点边长,扩散系数D_m为2.2×10^-9 m²/s。在模拟开始时,将10⁵个溶质粒子均匀放置于介质上边界下方0.5 mm<x<11.5 mm、0.5 mm<y<11.5 mm和2 mm<z<22 mm范围内的孔隙部分,溶质粒子初始位置、t=1 000 s和t=7 000 s时的空间分布(剖面图)如图1b所示。孔隙与固体间的边界为反射边界,上边界、下边界以及其他四个侧面边界设置为周期性边界。

1.4 弥散系数计算

本文计算了所有介质中垂直于水流主要流动方向(x方向)上的溶质羽的空间二阶中心矩M_2,T。同时,为对比不同胶结度的多孔介质中溶质横向弥散与纵向弥散的差异,从而系统分析胶结度对溶质弥散行为的影响,本文同时计算了沿水流主要流动方向(z方向)上的溶质羽空间二阶中心矩M_2,L,横向、纵向上的二阶空间中心矩计算公式如式(2)、(3)所示。关于不同胶结度的多孔介质中溶质纵向的弥散行为及机理可参考文献[28],本文仅从与横向弥散对比的角度进行分析。

(2)

M_2,T=<[x_i(t)-x_i(0)-<x_i(t)-x_i(0)>]²>

(3)

M_2,L=<[z_i(t)-z_i(0)-<z_i(t)-z_i(0)>]²>

其中,x_i(t)、z_i(t)分别为t时刻溶质粒子i的x、z坐标,<·>则表示所有粒子相应变量的平均值。

为定量研究胶结度不同时介质孔隙尺度结构特征对溶质横向弥散速率的影响及其与对纵向弥散速率影响的差异,本文计算了不同介质中溶质横向、纵向弥散系数D_T、D_L:

(4)D_T= \frac{\Delta M_{2,T}}{2\mathrm{\Delta }t}; D_L= \frac{\Delta M_{2,L}}{2\mathrm{\Delta }t}

式中,ΔM_2,T、ΔM_2,L分别为时间Δt内横向、纵向二阶中心矩的差值。

2 结果及讨论

2.1 弥散程度分析

横向、纵向二阶空间中心矩的计算结果如图2和表2所示。从图2中可以看出,固体颗粒胶结度对横向二阶矩M_2,T随时间的演化有显著影响,但与其对纵向二阶矩M_2,L的影响又有所不同。首先,随着胶结度增大,同时刻的M_2,T的值显著增大。这表明当胶结度增大时,溶质横向弥散程度有所增强。但与纵向二阶矩M_2,L的计算结果对比可发现,横向弥散增强幅度要远小于纵向弥散。如当t=7 500 s时, PM43、PM35、PM30和PM25中M_2,T分别是2.23×10^-4、2.75×10^-4、3.46×10^-4和3.82×10^-4 m²,而M_2,L的值分别是3.27×10^-3、6.83×10^-3、1.04×10^-2和2.28×10^-2 m²。随着胶结度P_c从0%增大到8.13%,M_2,T的值增大了71.3%,与M_2,T增大了6.97倍相比,增大幅度较小。其次,从图2a中还可看出,在不同胶结度的多孔介质中,M_2,T均随时间线性增大。如在介质PM43、PM35、PM30和PM25中,M_2,T的幂律拟合指数分别为1.029、0.988、1.007和1.003,均近似为1。这表明在多孔介质中,溶质横向弥散均基本符合菲克定律。即便是在胶结度最大的介质PM25中,溶质横向弥散尺度效应亦较弱。这与M_2,L时间曲线随着胶结度增大的变化有所区别。图2b显示随着胶结度增大,M_2,L随时间的演化规律逐渐偏离线性,转变为随时间呈幂律增大的规律。如在胶结度最大的介质PM25中,M_2,L的幂律拟合指数为1.295。这表明随着胶结度的增大,溶质纵向弥散的尺度效应显著增强,与前期研究结论一致。

2.2 弥散速率分析

不同介质中溶质横向、纵向弥散系数D_T、D_L的计算结果如图3所示。从图中可看出,在初始阶段,不同胶结度的多孔介质中的D_T均有一个快速下降的阶段,其最大值并不是最终的渐进值,而是出现在运移初始阶段。经过较短时间的运移,所有介质中D_T值即趋于稳定。在不同胶结度的多孔介质中,D_T随时间的演化规律基本相同,区别仅在于D_T值的大小。介质PM43、PM35、PM30和PM25中,D_T的渐进值D_T,asym分别为1.50×10^-8、1.76×10^-8、2.23×10^-8和2.48×10^-8 m²/s。D_T时间曲线的以上特征表明,随着胶结度增大溶质横向弥散速率增大。但是,胶结度对溶质横向弥散尺度效应的影响较弱,在所有四种介质中溶质横向弥散系数均为常数,为菲克弥散。而从图3b中可看出,仅在完全不胶结的介质PM43中,D_L在较短时间内可达到渐进值并保持为常数。随着胶结度增大,D_L不仅表现为渐进值增大,达到渐进值所需时间也显著增长。在模拟结束t=7 500 s时,PM25中的D_L的值仍未达到渐进值。与前期关于纵向弥散的研究是一致的,即随着胶结度增大,在纵向上不仅溶质弥散速率增大,且将长期处于非菲克弥散状态,尺度效应显著增强。以上模拟结果表明不同胶结度的多孔介质孔隙尺度结构差异对溶质横向弥散速率有显著影响,且与对纵向弥散速率的影响不同。

如引言中所述,前期研究中常把横向弥散系数视为比纵向弥散系数小一个数量级的值。为进一步系统认识胶结度不同的介质中溶质的弥散速率变化规律,本文计算了不同胶结度介质中溶质纵向弥散系数与横向弥散系数的比值D_L/D_T,结果如图4所示。从图中可以看出,在初始阶段,D_L/D_T的值在固体颗粒完全不胶结的介质中最大。而在固体颗粒轻微胶结、严重胶结的介质中,D_L/D_T值会略小于PM43中的值。随着溶质的运移,所有介质中D_L/D_T的值快速增大,且胶结度越大,D_L/D_T值随时间增长得越快。固体颗粒胶结的介质中的D_L/D_T值在短时间内就增大到高于介质PM43中的水平。如在t=42 s时,介质PM43、PM35、PM30和PM25中,D_L/D_T的值分别为4.80、4.94、4.97和5.82,固体颗粒胶结的介质中D_L/D_T的值均大于PM43中的值。随着溶质运移,所有介质中D_L/D_T值会逐渐达到一个渐进值。从图4和表2中可看出,在固体颗粒完全不胶结的介质PM43中,D_L/D_T的渐进值为14.40,即D_T值基本为比D_L 小一个数量级的值。而随着胶结度的增大,D_L/D_T 的渐进值显著增大。如在固体颗粒轻微胶结的介质PM35、PM30中,D_L/D_T的渐进值分别为25.71、31.04。而在胶结度最大的介质PM25中, D_L/D_T的渐进值已高达77.59。此外,随着胶结度增大,D_L/D_T值达到渐进值所需的时间越来越长,如介质PM43、PM35中,在1 080 s时D_L/D_T的值已基本接近渐进值。介质PM30中D_L/D_T基本接近渐进值的时间为1 896 s,而在胶结度最大的介质PM25中,在模拟结束t=7 500 s时才逐渐接近渐进值。这说明在运移初期开始所有介质中溶质纵向弥散速率与横向弥散速率的差异程度迅速扩大,且胶结度越大差异程度增大越快,持续时间越长,最终达到稳态时两方向上弥散速率的差异程度也越大。从以上D_L/D_T值随时间的变化规律可看出,简单将横向弥散系数D_T视为比D_L小一个数量级的值是不恰当的,特别是在固体颗粒胶结的介质中。首先,即便是在相对均质的固体颗粒完全不胶结的介质PM43中,D_L/D_T值也不是一个常数,而是会随着溶质的运移不断增大并最后趋于稳定。随着胶结度增大,D_L/D_T值处于逐渐增大状态的时间显著增长。其次,即便经过较长时间运移D_L/D_T值达到渐进值后,将D_T视为一个比D_L小一个数量级的值,也与溶质实际弥散行为存在较大偏差。

2.3 弥散行为机理分析

前期研究发现,流速概率密度分布等流场特征可较好揭示不同胶结度的介质中溶质纵向弥散行为差异的内在机制^[28-29]。为进一步系统探究不同胶结度的介质中溶质横向弥散行为及其与纵向弥散行为差异的内在机制,本文统计了介质中地下水横向(x方向)、纵向(z方向)归一化流速概率密度分布函数,结果如图5所示。

从图中可看出,随着胶结度增大,介质中横向、纵向上的归一化流速分量概率密度分布函数均越来越发散。局部流速与平均流速相近区域的比例显著减小,而远大于、远小于平均流速的区域的比例显著增大。但是与纵向上相比,流速横向分量的以上变化趋势相对较弱。如在PM43、PM35、PM30和PM25四种介质中,|v_x|/|v_x|_avg值为1时所对应的比例值分别是6.94%、6.53%、6.25%和5.96%,而|v_z|/|v_z|_avg为1时所对应的比例值分别为7.24%、6.52%、5.89%和5.19%。随着胶结固体颗粒占比从0增大到8.1%,横向上局部流速值位于其平均值附近的比例减少了9.8%,而纵向上局部流速值位于其平均值附近的比例减少了20.5%,纵向上减少的比例是横向上减少的2.09倍。这表明两个方向上流速分量的非均质性均有所增强,且横向上流速非均质性增强相对较弱。流场特征的以上变化可较好地解释胶结度对溶质横向弥散行为的影响及其与对纵向弥散影响差异的内在机制。首先,随着胶结度增大,纵向、横向上流速分量非均质性增强,流速概率密度分布函数越来越发散,导致不同位置处溶质对流速率差异显著增大。最终引起溶质纵向、横向弥散程度均显著增大。其次,纵向上流速分量非均质性增强显著而横向上流速非均质性增强相对较弱,导致随着胶结度增大,虽然溶质横向弥散系数增大,横向弥散程度增强,但是与纵向弥散相比,增强幅度相对较小,且未展现明显的尺度效应。为进一步验证以上分析,本文统计了不同时刻溶质粒子横向、纵向对流速度概率分布,结果如图6所示。

从图6可以看出,胶结度增大引起的流场流速概率密度发散对溶质对流速率的影响是显著的。在完全不胶结的、轻微胶结的介质PM43、PM35中,溶质在横向、纵向上的对流速率概率密度分布基本与地下水流速概率密度分布一致;而在胶结度较大的介质PM30、PM25中,溶质粒子在横向、纵向上的概率密度分布仅在运移初期与地下水流速概率密度分布基本一致。随着溶质运移,x方向、z方向上的对流速率接近流体流速分量平均值的溶质比例逐渐减小,且该比例要小于流体流速分量平均值附近区域的比例。而对流速率远小于平均流速(缓流区)、远大于平均流速(优势渗流区)的溶质的比例逐渐增大,且二者的比例高于流体缓流区、优势渗流区的比例。以上结果表明在胶结度较小的介质中所有的溶质粒子能够较快经历所有的流速值,从而溶质对流速率概率密度分布能够快速接近流体流速概率密度分布并保持一致。而在胶结度较大的介质中,由于流速概率分布越来越发散,部分溶质进入缓流区。进入缓流区的溶质粒子对流速率较小,需要经过较长时间才能运移出该区域,而随着溶质运移有越来越多的溶质进入缓流区域,最终导致处于缓流区域的溶质比例增加。同时,也有更多的溶质进入流速相对较大的优势渗流区域。由于对流速率较大,进入优势渗流区域的溶质多沿优势渗流通道快速运移,再次运移出优势渗流区的溶质比例较小,最终导致优势渗流区的溶质比例也逐渐增大。缓流区域、优势渗流区域的溶质比例同时增大,最终导致流速值位于平均流速区域的溶质比例减少。溶质对流速率差异变大最终导致溶质横向、纵向上的弥散程度增大。另外从图6可以看出,PM25介质中,虽然随着溶质运移,两个方向上的溶质弥散速率均逐渐发散,但是在两个方向上溶质对流速率有所不同。在横向上虽然弥散速率逐渐发散,但是在缓流区没有出现显著的峰值,位于平均流速区域的溶质粒子的比例仍然是最高的;而在纵向上,在缓流区出现了明显的峰值,甚至位于缓流区域的溶质粒子的比例要大于位于平均流速区域的溶质粒子的比例。这引起溶质粒子在纵向上的对流速率差异要显著大于横向上对流速率差异。这最终导致:随着胶结度增大,溶质在横向、纵向上的弥散度均有所增强,但横向弥散增强幅度相对较小,且未展现明显的尺度效应;而纵向弥散程度增强显著,且展现了非菲克弥散现象。

3 结论

本文运用孔隙尺度模拟方法研究了不同胶结度的多孔介质中溶质横向弥散特征。得到如下结论:

(1)随着胶结度增大,溶质横向弥散程度有所增强,但与纵向弥散相比增强幅度较小。如当t=7 500 s时,随着胶结度P_c从0%增大到8.13%,横向上溶质空间二阶矩M_2,T仅增大了71.3%,而M_2,L增大了6.97倍。溶质空间二阶矩时间曲线、溶质横向弥散系数的模拟结果均表明与纵向弥散随胶结度增大尺度效应增强不同,当胶结度增大时,溶质横向弥散未展现显著的尺度效应。溶质空间二阶矩时间曲线的幂律拟合指数均在1左右,即多孔介质中溶质空间二阶矩均随时间线性增大,且多孔介质中溶质横向弥散系数均在较短的时间内达到稳定值并保持为常数。

(2)纵向弥散系数D_L与横向弥散系数D_T的比值则表明,将D_T视为比D_L小一个数量级的量与实际的溶质弥散行为存在较大偏差。这是因为多孔介质中D_L/D_T值不是常数,而会随着溶质运移逐渐增大并最终达到渐进值,且胶结度越大,D_L/D_T值处于不稳定状态的时间越长。此外,即便是达到稳态以后,横向弥散系数与纵向弥散系数之间的差异程度也因受介质结构的影响而具有较大差别:在介质PM43中D_L/D_T的最终渐进值为14.40, D_T值基本为比D_L 小一个数量级的值;而在PM25中, D_L/D_T的渐进值已高达77.59。

(3)流体流速概率密度模拟结果表明,胶结度增大时,流体横向流速分量概率密度分布越来越发散,非均质性增强。但是,与纵向流速分量相比,流速横向分量的变化程度相对较小。通过统计不同胶结度介质中溶质对流速率概率密度发现这正是不同胶结度的介质中溶质横向弥散特征与纵向弥散特征不同的内在机制。

衷心感谢各位专家及编辑在审稿过程中对本文提出的宝贵修改意见。

参考文献

[1] ROLLE M, HOCHSTETLER D, CHIOGNA G, et al. Experimental investigation and pore-scale modeling interpretation of compound-specific transverse dispersion in porous media[J]. Transport in Porous Media, 2012, 93(3): 347-362.

[2] SCHEVEN U M. Pore-scale mixing and transverse dispersivity of randomly packed monodisperse spheres[J]. Physical Review Letters, 2013, 110(21): 214504.

[3] BIJELJIC B, BLUNT M J. Pore-scale modeling of transverse dispersion in porous media[J]. Water Resources Research, 2007, 43(12): W12S11.

[4] KANG P K, DE ANNA P, NUNES J P, et al. Pore-scale intermittent velocity structure underpinning anomalous transport through 3-D porous media[J]. Geophysical Research Letters, 2014, 41(17): 6184-6190.

[5] 王景瑞, 赵建世, 胡诗若. 地下水溶质反常运移的分数阶对流扩散模型研究进展[J]. 中国环境科学, 2022, 42(12): 5845-5855.

[6] 郭芷琳, 马瑞, 张勇, 地下水污染物在高度非均质介质中的迁移过程: 机理与数值模拟综述[J]. 中国科学: 地球科学, 2021, 51(11): 1817-1836.

[7] DENTZ M, BORGNE T L, ENGLERT A, et al. Mixing, spreading and reaction in heterogeneous media: a brief review[J]. Journal of Contaminant Hydrology, 2011, 120/121(Spec): 1-17.

[8] CHIOGNA G, CIRPKA O A, GRATHWOHL P, et al. Transverse mixing of conservative and reactive tracers in porous media: quantification through the concepts of flux-related and critical dilution indices[J]. Water Resources Research, 2011, 47(2): W02505.

[9] DELGADO J M P Q. Longitudinal and transverse dispersion in porous media[J]. Chemical Engineering Research and Design, 2007, 85(9): 1245-1252.

[10] BOON M, BIJELJIC B, NIU B, et al. Observations of 3-D transverse dispersion and dilution in natural consolidated rock by X-ray tomography[J]. Advances in Water Resources, 2016, 96: 266-281.

[11] SONNENWALD F, HART J R, WEST P, et al. Transverse and longitudinal mixing in real emergent vegetation at low velocities[J]. Water Resources Research, 2017, 53(1): 961-978.

[12] HOU Y, HU B X, LIU S, et al. Using a pore-scale modelling approach to study solute dilution process through cemented porous media[J]. Geofluids, 2022, 2022: 8944803.

[13] ROLLE M, CHIOGNA G, HOCHSTETLER D L, et al. On the importance of diffusion and compound-specific mixing for groundwater transport: an investigation from pore to field scale[J]. Journal of Contaminant Hydrology, 2013, 153(2013): 51-68.

[14] DU Z P, CHEN J J, KE S N. Transverse mixing zone under dispersion in porous media: effects of medium heterogeneity and fluid rheology[J]. Physics of Fluids, 2023, 35(4): 043105.

[15] ZHANG D, LIU T, PRIGIOBBE V. Enhanced solute transport in porous media due to pH-dependent adsorption and transverse dispersion[J]. Advances in Water Resources, 2022, 164: 104195.

[16] BOON M, BIJELJIC B, KREVOR S. Observations of the impact of rock heterogeneity on solute spreading and mixing[J]. Water Resources Research, 2017, 53(6): 4624-4642.

[17] DOLAMORE F, FEE C, DIMARTINO S. Modelling ordered packed beds of spheres: the importance of bed orientation and the influence of tortuosity on dispersion[J]. Journal of Chromatography A, 2018, 1532: 150-160.

[18] 叶逾, 张宇, 蔡芳敏, 宏观各向异性多孔介质中流体变形及溶质运移[J]. 水科学进展, 2021, 32(6): 903-910.

[19] 宋依依, 曹阳, 段鑫盈, 秸秆还田深度对土壤团聚体组成及有机碳含量的影响[J]. 土壤, 2022, 54(2): 344-350.

[20] BLUNT M J, BIJELJIC B, DONG H, et al. Pore-scale imaging and modelling[J]. Advances in Water Resources, 2013, 51: 197-216.

[21] ZHOU J Q, GUO L G, JIAO J J, et al. Geometry-based prediction of solute transport process in single 3D rock fractures under Laminar flow regime[J]. Journal of Geophysical Research: Solid Earth, 2023, 128(3): e2022JB025542.

[22] 周鸿翔, 郑延丰, 吴劳生, 孔隙尺度多孔介质流体流动与溶质运移高性能模拟[J]. 水科学进展, 2020, 31(3): 422-432.

[23] PILOTTI M. Generation of realistic porous media by grains sedimentation[J]. Transport in Porous Media, 1998, 33(3): 257-278.

[24] KANG Q, LICHTNER P C, ZHANG D. Lattice Boltzmann pore-scale model for multicomponent reactive transport in porous media[J]. Journal of Geophysical Research: Solid Earth, 2006, 111(B5): B05203.

[25] ZHOU H, YU X, CHEN C, et al. Evaluating hydraulic properties of biochar-amended soil aggregates by high-performance pore-scale simulations[J]. Soil Science Society of America Journal, 2018, 82(1): 1-9.

[26] WOLF-GLADROW D A. Lattice-gas cellular automata and lattice Boltzmann models: an introduction[M]. Berlin: Springer, 2004.

[27] SRINIVASAN G, TARTAKOVSKY D M, DENTZ M, et al. Random walk particle tracking simulations of non-Fickian transport in heterogeneous media[J]. Journal of Computational Physics, 2010, 229(11): 4304-4314.

[28] HOU Y S, WU J, JIANG J. Time behavior of anomalous solute transport in three-dimensional cemented porous media[J]. Soil Science Society of America Journal, 2019, 83(4): 1012-1023.

[29] LU C, WANG Z, ZHAO Y, et al. A mobile-mobile transport model for simulating reactive transport in connected heterogeneous fields[J]. Journal of Hydrology, 2018, 560: 97-108.

基金资助

国家自然科学基金项目(42002257)

山东省自然科学基金项目(ZR2020QD123)

山东省新高校院所创新团队项目(2021GXR070)

山东省高校院所创新团队项目(2018GXRC012)