多孔介质固体颗粒粒径特征及膨胀作用对地下水孔隙尺度流场特征的影响

Influence of Particle Size Characteristics and Swelling of Solid Particles in Porous Media on Pore-Scale Flow Field Characteristics of Groundwater

摘要

地下水孔隙尺度流场特征的研究对于深入理解地下水渗流、溶质运移具有重要意义.目前固体颗粒粒径特征不同时介质微观结构对孔隙尺度流场的影响尚不清楚.基于迭代重排算法构建了固体颗粒粒径分布特征、膨胀程度不同的多孔介质，基于此研究了固体颗粒平均粒径、粒径方差及膨胀作用对地下水流场特征的影响.结果表明，在介质孔隙率相同的条件下，固体颗粒平均粒径、粒径方差对多孔介质中流速的非均质性、速度概率密度分布等流场特征影响较小.而在同一介质中，当固体颗粒膨胀引起孔隙率减小时，平均粒径较小幅度的变化就会对以上流场特征产生显著影响.例如当固体颗粒膨胀程度增强时，粒径小幅度的增大，就会导致流场中优势流区和不流动区比例同时大幅度增加，流场非均质性显著增强，流速概率密度分布更加发散.

关键词

地下水流场特征 / 多孔介质 / 固体颗粒粒径 / 孔隙尺度 / 数值模拟 / 水文地质学

Key words

groundwater flow field characteristics / porous media / solid particle size / pore scale / numerical simulation / hydrogeology

0 引言

饱和多孔介质中地下水流场特征的研究为相关领域所关注.从多角度实现多孔介质结构精细化分析，探究孔隙结构变化对地下水流场空间分布特征的影响，在研究地下水污染物迁移转化和进行地下水资源评价中尤为重要（惠威等，2020；柴波等，2022；周凯等，2022）.而目前多在达西尺度上采用数值模拟和观测统计等方法研究地下水流动，介质孔隙尺度结构差异对地下水流动影响的研究相对较少（Bijeljic et al.，2013a）.近期的研究表明介质孔隙尺度结构差异，如介质的空隙比和固体颗粒粗糙度、孔隙连通性、孔喉比等孔隙尺度结构特征对渗流场特征有显著影响（Dentz et al.，2004；Sánchez-Vila and Carrera，2004；赵贤正等，2023）.Edery et al. （2014）证明流场特征是介质统计参数与流动条件组合控制的结果，只基于关键路径分析和渗流理论的模型是不充分的.另外，孔隙尺度上地下水流速空间分布等流场特征也正是达西尺度上地下水渗流现象的本质（Lee et al.，2018）.魏赫鑫等（2020）利用低温液氮吸附和核磁共振实验证明了岩石介质的平均孔径及分形维数是影响致密气储层流体分布的主要因素之一；Aziz et al.（2020）发现在不同湿润条件下孔隙形态对两相流的相互驱替效果影响显著；姚兰兰（2021）基于对页岩油储层微观孔隙结构特征的定量表征，发现页岩样本中孔隙结构的连通性和大孔喉占比与渗透率呈正相关.多孔介质固体颗粒粒径大小及其概率分布决定了介质中孔隙的连通性和水力传导性能，Li et al.（2019， 2020）通过砂柱实验发现介质颗粒粒径大小和颗粒分选度对介质的导水率和流态有着不可忽视的影响.因此，深入研究多孔介质的粒径分布特征对孔隙尺度上地下水流场的影响，有助于在宏观尺度上深入理解地下水渗流现象及其内在机理（Li et al.，2017）.

另外，近期有多项研究表明多孔介质孔隙尺度结构差异对地下水流场的影响最终会导致地下水中溶质运移行为显著不同（刘咏等，2022）.Sharma et al.（2022）在研究多孔介质中的non-Fickian运移时，发现土壤充填物的不均匀性和土壤介质粒度的孔隙尺度不均匀性引起的流场变化是导致non-Fickian运移的重要因素；Bijeljic et al.（2013b）通过对比三组不同粒级多孔介质中的溶质运移现象，发现介质中有效不流动区大小对其运移的非菲克程度有较大影响.

诸如上述地下水渗流及溶质运移特征的变化，常与流体流速概率密度分布密切相关.其中优势流区和不流动区的占比和空间分布更是流场的主要特征，且与non-Fickian运移的提前穿透和拖尾现象有密切关联（Srzic et al.，2013；Rolle and Kitanidis，2014；张学羿和窦智，2018）.优势流区是地下水流场中流线最为聚集的位置，周围流体迅速汇集于此，产生显著的聚流效应（Werth et al.，2006；徐鹏等，2017）.在聚流效应下，流体中的溶质也会随之快速运移，这是导致溶质提前穿透的主要原因之一（Willingham et al.，2010）.与之相反，不流动区则是流场中流体近乎停滞的区域，该区域溶质运移主要以扩散为主，运移速度缓慢.不流动区的大幅增加，会使多孔介质对地下水污染物的滞留效果大幅加强，并导致溶质运移拖尾现象的显著增强（Wirner et al.，2014）.Di Palma et al.（2017）通过模拟存在渗透率较小的透镜体的含水介质中的溶质运移情况，发现在透镜体处形成的缓流区会导致溶质运移的重拖尾现象.因此，深入研究介质孔隙尺度结构差异对流场中优势流区和不流动区的影响，对探究地下水渗流和溶质运移规律具有重要意义.

目前，关于固体颗粒粒径大小、固体颗粒分选程度不同时多孔介质孔隙尺度结构差异对地下水流场特征影响的研究还较少.因此，本文基于孔隙尺度模拟方法，深入研究了多孔介质固体颗粒平均粒径、粒径方差及固体颗粒膨胀作用强度对孔隙尺度上地下水流场特征的影响.

1 介质构建及介质特征表征

1.1　介质构建方法

本文基于改进的迭代重排算法构建了含水介质孔隙尺度结构模型.并利用颗粒半径概率密度分布曲线有效控制了介质颗粒粒径大小及分布范围.天然含水介质具有较强的非均质性，其颗粒形状具有高度的异质性，介质构成的控制条件较多，流道错综复杂.在研究多孔介质参数对流场变化的影响时，难以控制单一变量对流场产生的影响.为此，本文采用二维球型介质模型来研究孔隙尺度上的流场变化（Yang et al.，1996），可更直观地观察到不同介质中流场特征间的差异，且更容易控制介质的参数变化，为研究流场特征差异提供方便.

二维球型介质模型中用圆形粒子代表介质中的固体颗粒；由于天然土壤等含水介质中固体颗粒粒径概率密度多符合对数正态分布，本文中固体颗粒粒径由截断的对数正态分布决定（Molins et al.，2012）.具体构建过程如下：

（1）在研究区域内连续生成圆形固体颗粒，根据给定的对数正态分布的颗粒粒径概率密度函数确定每个颗粒的半径 r_i，随机确定固体颗粒初始位置，直至初始孔隙度 P_{\mathrm{0}}达到预设值，其中初始孔隙度 P_{\mathrm{0}}为：

P_{\mathrm{0}} = \mathrm{1} - {\sum }_{i=\mathrm{1}}^n \frac{\mathrm{\pi }{r}_i^{\mathrm{2}}}{S_{\mathrm{t}}}，（1）

其中， S_{\mathrm{t}}为介质区域的面积； r_i为固体颗粒i的半径；n为生成的固体颗粒个数.

（2）采用迭代重排算法消除颗粒重叠.在介质区域内沿y方向，坐标由小到大搜索重叠颗粒 G_i，然后，将该颗粒按向量 {\mathit{O}}_{\mathbf{t}}进行平移，并再次判断是否重叠，若重叠则继续重复上述执行过程，直至不再重叠.平移向量 {\mathit{O}}_{\mathbf{t}}为（Hochstetler et al.，2013；Hou et al.，2019）：

{\mathit{O}}_{\mathrm{t}} = {\mathit{O}}_x + {\mathit{O}}_y，（2）

{\mathit{O}}_x = {\mathit{i}}_x {\sum }_{k=\mathrm{1}}^n \left[r_i+r_k-d_{ik}\frac{x_i-x_k}{d_{ik}}\right]，（3）

{\mathit{O}}_y = {\mathit{i}}_y {\sum }_{k=\mathrm{1}}^n \left[r_i+r_k-d_{ik}\frac{y_i-y_k}{d_{ik}}\right]，（4）

d_{ik} = \sqrt[]{(x_i-x_k{)}^{\mathrm{2}}+(y_i-y_k{)}^{\mathrm{2}}}，（5）

其中 {\mathit{i}}_x、 {\mathit{i}}_y为x和y方向的单位矢量；n为与移动的固体颗粒 G_i重叠的颗粒的个数； r_i、 r_k分别为固体颗粒 G_i的半径和与其重叠的颗粒 G_k的半径；（ x_i， y_i）和（ x_k， y_k）为相互重叠固体颗粒的圆心位置.图1为运用迭代重排算法构建的孔隙率为0.45、平均粒径2 mm、粒径方差0.5 mm的多孔介质.

1.2　不同粒径特征的多孔介质构建

1.2.1　固体颗粒平均粒径不同的多孔介质构建

对于天然土壤等含水介质，一般将固体颗粒平均半径 的称为细砂、 的称为中砂、 \overline{R} > \mathrm{0} . \mathrm{25} \mathrm{m} \mathrm{m}的称为粗砂.由颗粒物堆积形成的土壤等多孔介质孔隙率多在0.44左右，且土地休耕和曝气后的孔隙率会进一步提高（Xue et al. 2023）.本文在80 mm×60 mm的介质区域范围内生成颗粒平均半径分别为 {\overline{R}}_{\mathrm{1}} =0.23 mm、 {\overline{R}}_{\mathrm{2}} =0.46 mm和 {\overline{R}}_{\mathrm{3}} = 0.74 mm的介质颗粒，并控制生成固体颗粒粒径的方差在0.22 mm左右.固体颗粒平均半径与天然含水介质中的中砂和粗砂的平均半径范围一致.多孔介质的初始孔隙率设置为0.38，最终孔隙率基本处于0.5±0.01范围内.在构建介质模型时为了避免出现固体颗粒过大或过小的不合理情况，特将颗粒半径分布范围分别控制在0.01~ 0.8 mm、0.2~1.0 mm和0.5~1.3 mm之间，最终构建固体颗粒平均粒径不同的三类多孔介质.对于每一粒级，构建三组介质模型作为参照.不同平均粒径的固体颗粒半径的概率密度分布如图2所示.

1.2.2　固体颗粒不同膨胀程度的多孔介质构建

本文在不改变固体颗粒分布位置且颗粒相互之间完全不重叠的前提下，将所有固体颗粒半径增加相同值，从而模拟多孔介质固体颗粒膨胀且不胶结的过程.随着固体颗粒粒径增大，孔隙率相应减小.在模拟固体颗粒膨胀过程中未改变其空间分布状态，能够更直观地反映颗粒膨胀作用对多孔介质流场分布特征的影响.由于固体颗粒相互之间未产生胶结，不会导致连通的孔隙转变为封闭状态，能够最大限度的保持孔隙原有的拓扑结构.

膨胀作用强度用变量 T_i来定量衡量：

T_i = \sqrt[]{\mathrm{1}-\frac{n_i}{n_{\mathrm{0}}}}，（6）

其中 n_i为固体颗粒膨胀后的多孔介质孔隙率； n_{\mathrm{0}}为介质的原孔隙率.本文通过这种方法生成了五组粒径膨胀强度分别为 T_{\mathrm{1}}=0、 T_{\mathrm{2}} =0.138、 T_{\mathrm{3}} =0.197、 T_{\mathrm{4}} =0.243、 T_{\mathrm{5}} =0.285的多孔介质，膨胀强度处于自由膨胀率小于10%的范围内（李先，2018；Wang et al. 2022）.膨胀强度与相应孔隙率如图3所示.

1.2.3　固体颗粒粒径分选度不同的多孔介质构建

多孔介质的分选性是指固体颗粒大小的均匀程度，本文主要以粒径的方差大小对分选性进行评价，方差越大，颗粒分布越分散、分选度越差，反之分选度越好（Qiao et al.，2020；Wang et al.，2020；Chen et al.，2021）.为了研究分选度对多孔介质流场分布的影响，本文构建了固体颗粒粒径方差分别为 {\sigma }_{d\mathrm{1}} =0.18 mm、 {\sigma }_{d\mathrm{2}} =0.24 mm、 {\sigma }_{d\mathrm{3}} =0.30 mm的三种多孔介质.在三种介质中固体颗粒半径分布范围均在0.2~1 mm之间，颗粒平均半径为0.46 mm，孔隙率为0.49，保证了粒径方差为唯一影响变量.同一方差下，生成三组多孔介质进行模拟对照，分选度不同的固体颗粒半径的概率密度分布如图4所示.

2 流场模拟

本文利用COMSOL Multiphysics 5.5软件中层流模块求解连续性方程及Navier-Stokes方程模拟多孔介质中孔隙部分的流场（Zhu et al.，2016；Hou et al.，2018）.连续性方程和Navier-Stokes方程如下所示：

\nabla \bullet \mathit{u} = \mathrm{0}，（7）

\rho \left(\mathit{u}\bullet \nabla \right) \mathit{u} = \nabla \bullet [ - p \mathit{l} + K ]，（8）

K = \mathit{\mu } [ \nabla \bullet \mathit{u} + {\left(\nabla \bullet \mathit{u}\right)}^{\mathrm{T}} ]，（9）

其中 \mathit{u}为速度矢量； \rho为水的密度； p为压力； \mathit{l}为单位向量； \mu为水的动力粘度.在模拟过程中运用稳态求解器进行求解，以保证模拟结束时流场达到稳定状态.模型的左边界、右边界为给定水头边界，水体在水头差的作用下，从左向右流动.上下边界和内部流固边界为无流动边界.通过调整压力水头的大小，获得孔隙部分流体平均流速分别为10 m/d、1 m/d和0.1 m/d的流场.

3 结果分析

3.1　固体颗粒平均粒径差异对流场的影响

为了研究固体颗粒平均粒径对多孔介质中优势流区、不流动区占比等流场特征的影响，本文在中砂和粗砂粒级范围内模拟了固体颗粒平均半径分别为0.23 mm、0.46 mm和0.74 mm的多孔介质中主流方向平均流速分别为10 m/d、1 m/d、0.1 m/d的流场，其流速概率密度分布如图5所示.为了更好地表现流场速度分布的差异性，流速概率密度分布图绘制为双对数坐标图.从图5中可看出，不同颗粒半径下，流速概率密度分布差异性较小，只有在流速极其微小的区域流速分布有小幅度改变，但并不影响整体的流速概率分布状态.另外，当平均流速不同时，速度的概率密度分布函数基本不变.以上结果表明，当介质固体颗粒粒径概率分布、孔隙率相同时，介质颗粒平均粒径的变化对多孔介质中流体流速概率分布特征的影响较小.

如引言中所述，地下水流场中优势流区、不流动区对地下水渗流及溶质运移具有显著影响.通常将流速大于平均流速5倍的区域定义为优势流区，将流速小于100倍平均流速的区域定义为不流动区（Dou et al.，2018）.本文统计不同介质中优势流区、不流动区的占比，如表1所示.结果表明随着粒径的增大，随着固体颗粒平均半径增大，三种平均流速下流场中优势流区和不流动区的变化趋势、变化范围基本一致，均有小幅度的增加.当平均半径由0.23 mm增加至0.46 mm，优势流区平均占比由0.28%增加至0.39%，不流动区平均占比由3.07%增加至3.28%；当平均半径由0.46 mm增加至0.74 mm，优势流区平均占比由0.39%增加至0.45%，不流动区平均占比由3.28%增加至3.56%.整体来看，优势流区和不流动区的大小随平均粒径的增大存在增加的趋势，但是整体幅度较小，对整体上的流速概率密度分布特征影响有限.这是因为随着粒径的增加，固液边界接触长度减小，导致边界层对流动的影响减小，如表1所示.随着介质颗粒平均粒径的增加，固体颗粒表面流固边界层长度随之减小，进而导致流场内优势流区域的增加.但是因为粘滞性对流体流动的影响大多集中在边界层位置，其作用范围较小，所以优势流区域增加的幅度并不明显.而不流动区面积的微小增加，主要是由于多孔介质粒径的增加，使固体颗粒与介质研究区域上下不流动边界相切和相交位置处死端孔隙面积增加，从而导致流场中不流动区面积的微小增加.

同时，为评价粒径变化情况下流场分布的均匀程度，本文采用流速的变异系数CV_u 作为速度场空间变异性的整体度量参数，对不同粒径的介质中的流场均质性进行评价，CV_u 值越大流场非均质性越强，公式如下：

C V_u = \frac{{\sigma }_u}{\overline{u}} = \frac{\sqrt[]{{\int }_{V_{\mathrm{p}}}^{}{(u-\overline{u})}^{\mathrm{2}}\mathrm{d}{V}_{\mathrm{p}}/V_{\mathrm{p}}}}{{\int }_{V_{\mathrm{p}}}^{}u\mathrm{d}{V}_{\mathrm{p}}/V_{\mathrm{p}}}，（10）

其中， {\sigma }_u为流速在水平方向上分量的方差； \overline{u}为流速的平均值； V_{\mathrm{p}}为流体占据的孔隙空间体积.从表1可以看出随固体颗粒粒径的增加，CV_u 值从0.947增至0.994，总体变化幅度较小，这说明在固体颗粒粒径概率密度分布、介质孔隙率不变的情况下，不同平均粒径的多孔介质中流场的非均质性差异很小，这与上文优势流和不流动区小幅的增加情况一致.

3.2　固体颗粒膨胀作用对流场的影响

3.1节的模拟结果表明在介质孔隙率、固体颗粒粒径概率密度分布相同的条件下，平均粒径的变化对流场的非均质性、优势流区、不流动区占比影响较小.本节模拟固体膨胀作用，研究当固体颗粒粒径增大导致孔隙率减小时，固体颗粒平均粒径差异对流场的影响.本文选取固体颗粒平均半径为 0.74 mm的多孔介质，在不改变颗粒空间位置且保证固体颗粒完全不重叠的条件下，将所有固体颗粒粒径增大相同值，得到对应固体颗粒膨胀作用强度分别为 T_{\mathrm{1}} = \mathrm{0}、 T_{\mathrm{2}} =0.138、 T_{\mathrm{3}} =0.197、 T_{\mathrm{4}} =0.243、 T_{\mathrm{5}} =0.285的五组介质模型，并分别模拟其在不同平均流速条件下的流场，其中膨胀作用强度为 T_{\mathrm{1}} = \mathrm{0}、 T_{\mathrm{3}} =0.197、 T_{\mathrm{5}} =0.285的流场模拟图如图6所示.

如图6所示，当 T_{\mathrm{1}}=0时，流场中流速分布较为均匀，没有明显的优势渗流通道出现，随着膨胀效果的增加，逐渐出现较为显著的优势渗流通道.在主流方向平均流速为1 m/d的情况下，固体颗粒未膨胀的多孔介质流场最大流速为 u_{\mathrm{m}\mathrm{a}\mathrm{x}} =11 m/d，而随着膨胀效果增加至 T_{\mathrm{5}} =0.285时，流速上限增至 u_{\mathrm{m}\mathrm{a}\mathrm{x}} =22 m/d.表明随着膨胀程度的加强，多孔介质孔隙率降低，多孔介质中流速分布范围增大，最大流速显著增大.从图6a、6b中可看出，在固体颗粒膨胀程度最大的介质中（ T_{\mathrm{5}} =0.285），优势渗流区域显著增大，原先零散分布的优势渗流区相互连接形成较为连续的渗流通道.这表明随着膨胀程度的加强，优势渗流通道在整个地下水流动中起到更为显著的作用，即有更多的地下水通过优势渗流通道流动.

同时，从图6中可看出，随着膨胀作用的加强，多孔介质中流速较小的区域也显著地增加.为了定量地研究优势流区和不流动区变化的差异性，本文统计了不同膨胀作用程度介质中的流速概率密度函数，分别绘制了流速概率密度分布的常数坐标及半对数坐标图，如图7所示.并统计了不同膨胀程度的介质中流场的参数，包括流速最大值与流速平均值的比值、优势流区及不流动区占比、流场变异系数、速度概率密度峰值等，如表2所示.模拟结果表明在基于迭代重排算法构建的多孔介质中模拟得到的地下水流场与Bijeljic （2013b）通过CT扫描获取的天然含水介质结构中模拟得到的流场具有相似的流速概率分布.

结果表明，当固体颗粒膨胀程度从 T_{\mathrm{1}}=0增加到 T_{\mathrm{5}} =0.285时， u_x/{\overline{u}}_x<0.01的不流动区占比扩大了2倍多，从3.5%增加至8.1%； u_x/{\overline{u}}_x>5的优势流区占比扩大了6倍，从0.4%增大至2.59%，优势流区和不流动区占比增加明显.另外，流速接近平均流速的区域所对应概率密度峰值随膨胀程度的加强而逐步下降，从7.21%降低至5.84%.这表明流场的流速分布逐渐发散，流速接近平均流速的区域减少，局部流速偏离平均流速的区域增加，流场的非均质性增强.为了量化流场的非均质性程度，本文计算了不同膨胀程度的多孔介质流场中的流速变异系数CV_u，如表2所示，当膨胀程度从 T_{\mathrm{1}}=0逐步加强到 T_{\mathrm{5}} =0.285时，CV_u 值由0.998增至1.383，流场的非均质性显著增加.

同时本文研究了膨胀作用强度不同时，介质中的流场特征差异产生的内在机理.多孔介质中固体颗粒的膨胀作用会压缩当前多孔介质的流道空间，进而导致部分孔隙孔径变小，靠近固体颗粒表面的边界层对地下水流动的影响增强，极大地阻碍了流体的流动，使流体通过该区域的速度变慢甚至近乎不再流动.如图8A所示，因膨胀作用孔喉直径进一步缩小，导致在主流方向上，部分孔径较小的单体孔喉区域的流体流动十分缓慢，增大了不流动区的占比；而类似于图8B所示的情况则是因为膨胀后产生的小直径孔喉的随机组合，会在流场中分隔出部分与主要流动区水力联系较弱的区域，进一步增大不流动区的占比.同时，随着膨胀作用的增强，孔喉数量增加，易流通道减少，导致更多的流体涌向优势渗流通道以更大的速度流动，最终导致优势流区域和不流动区同时增长，流场非均质性显著增强.

3.3　固体颗粒分选性变化对流场的影响

本文在粗砂粒径范围内，模拟了平均粒径为0.73 mm、孔隙率为0.484，但粒径方差分别为 {\sigma }_{d\mathrm{1}} =0.18 mm、 {\sigma }_{d\mathrm{2}} =0.24 mm、 {\sigma }_{d\mathrm{3}} =0.30 mm的三种多孔介质中的流场，并控制其平均流速分别为 10 m/d、1 m/d、0.1 m/d.随着固体颗粒粒径方差的增大，多孔介质粒径差异性变大、分选度变差，但多孔介质中的流速概率密度分布特征基本不变，对应的优势流区域和不流动区的面积，平均占比分别为3.55%、0.396%，且其只因为不同含水介质模型之间孔隙构造的细微差异性发生小范围的上下波动变化.

从整体的流场分布来看，多孔介质粒径方差的变化并没有引起整体的介质流道宽度和数量发生变化，多孔介质中的孔隙依旧保持着较好的连通性，流场内不同位置处流速大小的差异较小，流场分布较为均一，没有明显的流速过大和过小的区域.如图9所示，粒径方差不同的流场流速概率密度分布状态吻合较好，流场内的流速分布并没有随之发生本质性改变.同时，本文统计了粒径方差由小到大的多孔介质流速变异系数 C V_u值，如表3所示.其中 C V_{u_x}为流速在主流方向分量的变异系数值， C V_{u_y}为流速在垂直于主流方向分量的变异系数值.从中可以看出，变异系数值并未随粒径方差的增大产生明显差异和趋势性变化，这说明粒径方差大小控制的分选度差异对多孔介质中流场的非均质性影响较小.

4 结论

本文运用迭代重排算法构建了多孔介质孔隙尺度结构，结合数值模拟方法在孔隙尺度上分析了多孔介质固体颗粒粒径特征及固体颗粒膨胀作用对流场特征的影响.最终得到以下结论：

（1）在介质孔隙率不变的条件下，固体颗粒平均粒径的变化对流场特征影响较小.比如当孔隙率在0.5±0.01范围内时，随着固体颗粒平均粒径从 0.23 mm增大到0.74 mm，流场的非均质性以及优势流区、不流动区占比基本不变.但是，当固体颗粒粒径增大引起孔隙率减小时，即固体颗粒膨胀作用增强时，即便固体颗粒空间排列、拓扑结构保持不变，固体颗粒粒径小幅度的增加就会对流场的空间非均质性以及优势流区域、不流动区域占比产生显著影响.固体颗粒平均粒径增大、孔隙率减小时，流速概率密度分布更加发散，流场中接近平均流速的区域占比逐渐减小，而远大于平均流速和远小于平均流速的优势流区域及不流动区占比显著增加.比如随着固体颗粒膨胀程度从 T_{\mathrm{1}}=0增加至 T_{\mathrm{5}} =0.285，固体颗粒平均粒径仅增加了0.06 mm（从0.74 mm增大到0.80 mm），优势流区占比就扩大了6倍，从0.4%增大至2.59%，不流动区占比从3.5%增至8.1%.这是因为固体颗粒的膨胀作用增强，多孔介质的孔隙度减小，会进一步压缩流体的流动空间，导致流道变窄、孔喉增多，从而导致在部分孔喉区域形成不流动区.同时，随着流体流动难度的加大，流体会优先汇入连通性较好的优势渗流通道快速流动，从而导致优势流区的占比同时显著增大.

（2）粒径方差控制下的分选度差异不会显著影响多孔介质中地下水流场特征.随着固体颗粒粒径方差的增大，沿流体主要流动方向和垂直于流体主要流动方向上的流速分量概率密度分布函数没有明显的差异，且流速的变异系数并没有表现出与粒径方差趋势性变化的显著相关性.这主要是因为粒径方差的变化并未导致多孔介质内孔隙特征的显著变化，其变化前后介质内的孔隙连通性基本维持不变，流速分布较为均匀.

本文在孔隙尺度揭示了固体颗粒粒径大小、分选度及固体颗粒膨胀作用对地下水流场特征的影响，丰富了对流场变化微观机理的研究，可为深入研究多孔介质中地下水渗流现象、溶质运移行为及其内在机理提供参考.另外，含水介质膨胀强度受矿物组成等自身特性及温度等环境因素的影响，具有较大的变化范围.本文关于膨胀作用对渗流场影响的研究，主要关注于粒径增大但孔隙基本拓扑结构不变，即膨胀作用处于固体颗粒发生膨胀且尚未胶结的强度范围内.更强程度的膨胀作用对渗流场的影响尚需进一步研究.

参考文献

[1] Aziz, R., Niasar, V., Erfani, H., et al., 2020. Impact of Pore Morphology on Two-Phase Flow Dynamics under Wettability Alteration. Fuel, 268: 117315. https://doi.org/10.1016/j.fuel.2020.117315

[2] Bijeljic, B., Mostaghimi, P., Blunt, M. J., 2013a. Insights into Non-Fickian Solute Transport in Carbonates. Water Resources Research, 49(5): 2714-2728. https://doi.org/10.1002/wrcr.20238

[3] Bijeljic, B., Raeini, A., Mostaghimi, P., et al., 2013b. Predictions of Non-Fickian Solute Transport in Different Classes of Porous Media Using Direct Simulation on Pore-Scale Images. Physical Review E, 87: 013011. https://doi.org/10.1103/physreve.87.013011

[4] Chai, B., Shi, X.S., Du, J., et al., 2022. How to Realize Elaborated Analysis of Regional Rock Mass Structure? A Review and Idea. Earth Science, 47(12): 4629-4646 (in Chinese with English abstract).

[5] Chen, S. B., Gong, Z., Li, X. Y., et al., 2021. Pore Structure and Heterogeneity of Shale Gas Reservoirs and Its Effect on Gas Storage Capacity in the Qiongzhusi Formation. Geoscience Frontiers, 12(6): 101244. https://doi.org/10.1016/j.gsf.2021.101244

[6] Dentz, M., Cortis, A., Scher, H., et al., 2004. Time Behavior of Solute Transport in Heterogeneous Media: Transition from Anomalous to Normal Transport. Advances in Water Resources, 27(2): 155-173. https://doi.org/10.1016/j.advwatres.2003.11.002

[7] Di Palma, P. R., Parmigiani, A., Huber, C., et al., 2017. Pore-Scale Simulations of Concentration Tails in Heterogeneous Porous Media. Journal of Contaminant Hydrology, 205: 47-56. https://doi.org/10.1016/j.jconhyd.2017.08.003

[8] Dou, Z., Chen, Z., Zhou, Z. F., et al., 2018. Influence of Eddies on Conservative Solute Transport through a 2D Single Self-Affine Fracture. International Journal of Heat and Mass Transfer, 121: 597-606. https://doi.org/10.1016/j.ijheatmasstransfer.2018.01.037

[9] Edery, Y., Guadagnini, A., Scher, H., et al., 2014. Origins of Anomalous Transport in Heterogeneous Media: Structural and Dynamic Controls. Water Resources Research, 50(2): 1490-1505. https://doi.org/10.1002/2013wr015111

[10] Hochstetler, D. L., Rolle, M., Chiogna, G., et al., 2013. Effects of Compound-Specific Transverse Mixing on Steady-State Reactive Plumes: Insights from Pore-Scale Simulations and Darcy-Scale Experiments. Advances in Water Resources, 54: 1-10. https://doi.org/10.1016/j.advwatres.2012.12.007

[11] Hou, Y. S., Jiang, J. G., Wu, J. C., 2018. Anomalous Solute Transport in Cemented Porous Media: Pore-Scale Simulations. Soil Science Society of America Journal, 82(1): 10-19. https://doi.org/10.2136/sssaj2017.04.0125

[12] Hou, Y. S., Wu, J. C., Jiang, J. G., 2019. Time Behavior of Anomalous Solute Transport in Three-Dimensional Cemented Porous Media. Soil Science Society of America Journal, 83(4): 1012-1023. https://doi.org/10.2136/sssaj2018.12.0476

[13] Hui, W., Xue, Y.Z., Bai, X.L., et al., 2020. Influence of Micro-Pore Structure on the Movable Fluid Occurrence in Tight Sandstone Reservoir. Special Oil & Gas Reservoirs, 27(2): 87-92 (in Chinese with English abstract).

[14] Lee, J., Rolle, M., Kitanidis, P. K., 2018. Longitudinal Dispersion Coefficients for Numerical Modeling of Groundwater Solute Transport in Heterogeneous Formations. Journal of Contaminant Hydrology, 212: 41-54. https://doi.org/10.1016/j.jconhyd.2017.09.004

[15] Li, X., 2018. Study on the Expansion Characteristics of Neogene System in Xining Basin (Dissertation). Southwest Jiaotong University, Chengdu (in Chinese with English abstract).

[16] Li, Z. X., Wan, J. W., Huang, K., et al., 2017. Effects of Particle Diameter on Flow Characteristics in Sand Columns. International Journal of Heat and Mass Transfer, 104: 533-536. https://doi.org/10.1016/j.ijheatmasstransfer.2016.08.085

[17] Li, Z. X., Wan, J. W., Zhan, H. B., et al., 2019. Particle Size Distribution on Forchheimer Flow and Transition of Flow Regimes in Porous Media. Journal of Hydrology, 574: 1-11. https://doi.org/10.1016/j.jhydrol.2019.04.026

[18] Li, Z. X., Wan, J. W., Zhan, H. B., et al., 2020. An Energy Perspective of Pore Scale Simulation and Experimental Evidence of Fluid Flow in a Rough Conduit. Journal of Hydrology, 587: 125010. https://doi.org/10.1016/j.jhydrol.2020.125010

[19] Liu, Y., Zhang, Q., Qian, J.Z., et al., 2022. Simulation of Bimolecular Reactive Solute Transport in Porous Media via Image Analysis. Earth Science Frontiers, 29(3): 248-255 (in Chinese with English abstract).

[20] Molins, S., Trebotich, D., Steefel, C. I., et al., 2012. An Investigation of the Effect of Pore Scale Flow on Average Geochemical Reaction Rates Using Direct Numerical Simulation. Water Resources Research, 48(3): W03527. https://doi.org/10.1029/2011wr011404

[21] Qiao, J. C., Zeng, J. H., Ma, Y., et al., 2020. Effects of Mineralogy on Pore Structure and Fluid Flow Capacity of Deeply Buried Sandstone Reservoirs with a Case Study in the Junggar Basin. Journal of Petroleum Science and Engineering, 189: 106986. https://doi.org/10.1016/j.petrol.2020.106986

[22] Rolle, M., Kitanidis, P. K., 2014. Effects of Compound- Specific Dilution on Transient Transport and Solute Breakthrough: A Pore-Scale Analysis. Advances in Water Resources, 71: 186-199. https://doi.org/10.1016/j.advwatres.2014.06.012

[23] Sánchez-Vila, X., Carrera, J., 2004. On the Striking Similarity between the Moments of Breakthrough Curves for a Heterogeneous Medium and a Homogeneous Medium with a Matrix Diffusion Term. Journal of Hydrology, 294(1-3): 164-175. https://doi.org/10.1016/j.jhydrol.2003.12.046

[24] Sharma, P. K., Agarwal, P., Mehdinejadiani, B., 2022. Study on Non-Fickian Behavior for Solute Transport through Porous Media. ISH Journal of Hydraulic Engineering, 28(sup1): 171-179. https://doi.org/10.1080/09715010.2020.1727783

[25] Srzic, V., Cvetkovic, V., Andricevic, R., et al., 2013. Impact of Aquifer Heterogeneity Structure and Local-Scale Dispersion on Solute Concentration Uncertainty. Water Resources Research, 49(6): 3712-3728. https://doi.org/10.1002/wrcr.20314

[26] Wang, L. G., Zhang, Y. Z., Zhang, N. Y., et al., 2020. Pore Structure Characterization and Permeability Estimation with a Modified Multimodal Thomeer Pore Size Distribution Function for Carbonate Reservoirs. Journal of Petroleum Science and Engineering, 193: 107426. https://doi.org/10.1016/j.petrol.2020.107426

[27] Wang, L. L., Wang, Z. T., Ding, Z. P., et al., 2022. Factors Influencing Accuracy of Free Swelling Ratio of Expansive Soil. Journal of Central South University, 29(5): 1653-1662. https://doi.org/10.1007/s11771-022-4986-9

[28] Wei, H.X., Lai, F.P., Jiang, Z.Y., et al., 2020. Micropore Structure and Fluid Distribution Characteristics of Yanchang Tight Gas Reservoir. Fault-Block Oil & Gas Field, 27(2): 182-187 (in Chinese with English abstract).

[29] Werth, C. J., Cirpka, O. A., Grathwohl, P., 2006. Enhanced Mixing and Reaction through Flow Focusing in Heterogeneous Porous Media. Water Resources Research, 42(12): W12414. https://doi.org/10.1029/2005wr004511

[30] Willingham, T., Zhang, C. Y., Werth, C. J., et al., 2010. Using Dispersivity Values to Quantify the Effects of Pore-Scale Flow Focusing on Enhanced Reaction along a Transverse Mixing Zone. Advances in Water Resources, 33(4): 525-535. https://doi.org/10.1016/j.advwatres.2010.02.004

[31] Wirner, F., Scholz, C., Bechinger, C., 2014. Geometrical Interpretation of Long-Time Tails of First-Passage Time Distributions in Porous Media with Stagnant Parts. Physical Review E, 90(1): 013025. https://doi.org/10.1103/physreve.90.013025

[32] Xu, P., Li, C.H., Liu, H.C., et al., 2017. Fractal Features of the Effective Gas Transport Coefficient for Multiscale Porous Media. Earth Science, 42(8): 1373-1378 (in Chinese with English abstract).

[33] Xue, J. F., Qi, Z. W., Chen, J. L., et al., 2023. Dynamic of Soil Porosity and Water Content under Tillage during Summer Fallow in the Dryland Wheat Fields of the Loess Plateau in China. Land, 12(1): 230. https://doi.org/10.3390/land12010230

[34] Yang, A., Miller, C. T., Turcoliver, L. D., 1996. Simulation of Correlated and Uncorrelated Packing of Random Size Spheres. Physical Review E, 53(2): 1516-1524. https://doi.org/10.1103/physreve.53.1516

[35] Yao,L.L.,2021.Evaluation of Microscopic Pore Structure Characteristics and Flow Mechanism of Shale Oil Reservoirs (Dissertation).Institute of Porous Flow & Fluid Mechanics, Chinese Academy of Sciences, Langfang (in Chinese with English abstract).

[36] Zhang, X.Y., Dou, Z., 2018. Influence of Microscopic Pore Structure of Clay on Soluble Contaminant Transport. Hydrogeology & Engineering Geology, 45(4): 157-164 (in Chinese with English abstract).

[37] Zhao, X.Z., Chen, C.W., Song, S.Y., et al., 2023. Shale Oil Reservoir Structure Characteristics of the Second Member of Kongdian Formation in Cangdong Sag, Bohai Bay Basin. Earth Science, 48(1): 63-76 (in Chinese with English abstract).

[38] Zhou, K., Chen, X.P., Qu, X.B., 2022. Determining Methods of Micro-Pore and Liquids and Adsorption- Desorption Experiment for Shale Reservoir: Taking Member He-8 Reservoir in Ordos Yulin Area as an Example. Petroleum Geology & Oilfield Development in Daqing, 41(2): 139-146 (in Chinese with English abstract).

[39] Zhu, Y. H., Zhan, H. B., Jin, M. G., 2016. Analytical Solutions of Solute Transport in a Fracture-Matrix System with Different Reaction Rates for Fracture and Matrix. Journal of Hydrology, 539: 447-456. https://doi.org/10.1016/j.jhydrol.2016.05.056

基金资助

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

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

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

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

济南大学博士基金项目(XBS1911)

济南大学科技计划基金项目(XKY1922)