0 前言
随着计算机性能的增强，数值模拟技术可以有效模拟各类流体的流动机制和沉积过程。然而迄今为止数值模拟技术主要用于单一浊流模式[16-17]或者单一底流模式[18-19]，较少用来模拟研究重力流与底流交互作用对沉积物分布的控制效应及其沉积-侵蚀机制。Rebesco等[4]提出应用数值模拟研究重力流与底流交互作用的设想。本次研究借助商业有限元分析软件ANSYS，通过求解连续性方程和动量守恒方程，在建立重力流与底流交互作用三维数值模型的基础上，开展鲁伍马盆地上始新统重力流与底流交互作用过程中的水动力机制数值模拟研究，正演重力流与底流交互作用过程，探究底流对重力流沉积展布的影响，为重力流与底流交互作用研究提供数值模拟视角。研究成果对揭示深水重力流与底流交互的沉积动力学沉积-侵蚀机制具有重要的意义。
1 区域地质概况
鲁伍马盆地上始新统为一个三级层序，其顶界面H390、底界面H400表现为连续、稳定的强振幅地震反射特征，界面下削、上超特征明显。其中H400界面在水道发育位置为一水道侧向迁移所形成的侵蚀包络面（图2）。根据地震反射特征与沉积结构单元的组成及构型差异，在该三级层序内进一步识别出p1和p2两个四级层序界面，从而将上始新统三级层序划分为Ps1、Ps2及Ps3三个四级层序，Ps1和Ps2构成三级层序的低位体系域（LST），Ps3为三级层序的海侵+高位体系域（TST+HST），三个四级层序对应海底扇水道-朵叶复合体的3个沉积演化期次（图1b，图2）。Ps1时期为限制性水道-朵叶复合体发育阶段，该时期发育SW—NE走向的限制性水道-朵叶复合体（图1a），水道在剖面上呈“V”形，内部为中—强振幅、中连续地震反射，两侧与水道壁呈上超接触关系，与围岩弱振幅地震反射特征差异明显，水道两侧天然堤不甚发育。Ps2时期为半限制性迁移水道-朵叶体复合体发育阶段。相较于Ps1层序水道-朵叶复合体整体向南迁移，该时期体系末端的朵叶部分向南偏转（图1a）。水道内部呈中—强振幅、中连续地震反射，由3期依次向南迁移的水道构成（图2中的C1、C2、C3），水道北侧继承性发育宽度较大的溢岸天然堤/漂积体沉积；与北侧明显不同，水道南侧由于不断侵蚀，仅最后一期水道保留完整，且其南侧主要表现为侵蚀特征，天然堤体系不发育。Ps3时期为非限制性朵叶复合体发育阶段。朵叶复合体呈连续中—强振幅地震反射，外形略呈丘形，底部侵蚀不明显，南侧超覆于H400界面之上（图2）。
2 沉积动力学模拟与方法
2.1 沉积动力学模拟原理
本次研究基于雷诺平均Navier-Stokes方程在数值模拟中考虑湍流的影响[23-24]，分解流体瞬时湍流量，将湍流波动转换为雷诺兹应力，并应用动量方程进行求解[25]。ANSYS软件中的RNG k-ε模型（k为湍动能，ε为湍流耗散率）考虑低雷诺数效应和流体湍动所传递的二次相色散，设置相间湍动能传递和相位加权平均过程，使连续性方程中无需引入体积分数波动就能预测连续相的湍流，从而使涡流和快速应变流的模拟更为准确，能更好地揭示不同悬浮颗粒物质在环境流体中的流动过程。重力流与底流可以看成是环境流体和固体颗粒的多相流，而ANSYS软件提供的Mixture模型将不同相态在数学上视为相互穿透的连续体，且以不同的速度运动，能很好地揭示重力流与底流交互作用中混合沉积物与周围环境水体间复杂的相互关系，是本次研究的首选[24,26]。综上所述，本次研究借助ANSYS软件，基于RNG k-ε模型和Mixture模型对鲁伍马盆地上始新统沉积时期的重力流与底流交互作用过程进行沉积水动力机制数值模拟研究。
2.2 数据来源
将鲁伍马盆地上始新统的3个不同演化阶段视为数值模拟中的3个不同模型，分别对其进行模拟。各期几何模型数据见表1。
表1 沉积动力学模拟主要参数Table 1 Main parameters of the simulation of sedimentary hydrodynamics |
| 参数 | Ps1时期 | Ps2时期 | Ps3时期 |
|---|---|---|---|
| 模型长度/m | 400.00 | 400.28 | 400.96 |
| 模型宽度/m | 302.0 | 307.0 | 304.5 |
| 底流入流速度/(m·s-1) | 0.4 | ||
| 重力流入流速度/(m·s-1) | 1.5 | ||
| 砂质颗粒直径/mm | 0.08 | ||
| 泥质颗粒直径/mm | 0.005 | ||
| 砂泥比 | 4：1 | ||
| 海水密度/(kg·m-3) | 1 025 | ||
| 砂质颗粒密度/(kg·m-3) | 2 650 | ||
| 泥质颗粒密度/(kg·m-3) | 2 700 | ||
| 海水黏度/(kg·m-1·s-1） | 0.001 | ||
| 颗粒黏度/(kg·m-1·s-1） | 0.01 | ||
| 底面粗糙系数 | 0.8 | ||
| 底面粗糙高度/mm | 0.08 | ||
Heerema等[27]在蒙特利峡谷内观测到浊流的流速为0.5~8 m/s，本次研究应用现有的井震资料并通过重力流计算公式，结合模拟过程中扇体的展布变化特征，将重力流流速设置为1.5 m/s。海底实测数据显示，深层底流的流速平均0.2~0.4 m/s，最大可达1.2 m/s[28]。Miramontes等[29]通过水槽实验探究了不同底流流速对重力流与底流交互作用的影响，发现改变底流流速并不会对颗粒物造成不同的沉积响应，即共同表现出颗粒物顺底流溢出和水道逆底流迁移。因此，参考实测的底流流速范围，本次研究中给各模拟模型均设置了无底流流速和底流流速为0.4 m/s两种底流条件。环境流体密度设置为常规海水的密度1 025 kg/m3。参考Straub等[30]关于深水水道浊流的物理实验，将沉积颗粒体积浓度设置为2%。钻井资料揭示研究区上始新统Ps3时期砂泥岩含量之比约为4:1，因此在数值模拟中设置它们的占比分别为80%和20%。砂质颗粒密度为2 650 kg/m3，泥质颗粒密度为2 700 kg/m3，底面粗糙高度在数值上与砂质粒径一致（表1）。
2.3 模拟流程
本次模拟流程分为以下5步：①地质数据处理。通过井震资料建立层序地层格架，对关键层序界面数据进行时深转换、压实校正、负载均衡校正和古
水深校正，并将得到的地质数据在Surfer和Rhino软件中进行滤波分析和样条光滑处理，进而建立较光滑的研究区底界面地形（图3a）。②几何模型构建。在ANSYS的Geometry模块中建立不同时期的几何模型，并根据研究需要切割出研究区的范围（图3b）。③网格划分。在ANSYS的ICEM模块中，将几何模型进行结构网格划分，网格密度设置为1.5 m。为了更好地揭示重力流与底流交互作用底部的流体结构和沉积-侵蚀机制，将垂向网格进行向下递增加密，网格数变为原来的两倍，Ps1、Ps2和Ps3时期模型网格总数分别为290万、330万和320万个。网格质量越接近于1代表网格质量越好，模拟结果的收敛效果将越好。本次研究各期次的网格质量均大于0.9，表明模型的网格质量都很高，满足收敛条件（图3c）。④边界条件设定。在ANSYS的Fluent模块中设置了模型的边界条件，其中包括：入口为速度入口、顶部为对称边界、西侧为无滑移条件的固定壁面、南北两侧壁面设置为底流出入口、东侧为压力出口、底面为包含粗糙度的壁面。为了让流体以较稳定的初始状态进入模型，在模型西侧设置了长5 m的引流入口（图3b，3d）。⑤沉积过程模拟。在ANSYS的Fluent模块中按表1中的参数设置初始条件并进行模拟，适时更改时间步长，合理调控变量参数，以得到收敛性良好的模拟结果。
3 重力流与底流交互作用下的海底扇平面展布响应特征
3.1 模拟结果表征参数选取
在Walker(1978)提出的综合扇模式[31]中，细粒沉积物和粗粒沉积物在海底扇中具有相同的沉积展布趋势，但细粒沉积物（如粉砂岩和泥岩）在海底扇体系中具有更显著的远源沉积特征，而粗粒沉积物（如砂岩）并不会超出细粒沉积物的沉积范围。据此，为探究重力流与底流交互作用下海底扇整体的沉积展布，本次实验选择泥质颗粒浓度等值线的平面展布特征作为表征重力流海底扇空间分布形态的关键指标。鉴于细粒物质对底流作用的敏感性[5,32]，对比有底流与无底流条件下泥质颗粒浓度平面展布的变化，可以进一步分析重力流与底流交互作用的沉积响应过程。另外，在沉积动力学数值模拟中，流体携带的颗粒物质不断沉积和再悬浮，当迭代次数（流动时间）达到使悬浮沉积物浓度场呈现稳定状态时（即物质输运通量趋于平衡），认定系统进入沉积相态稳定阶段，并以此时刻的模拟结果作为各沉积期次的最终输出（图4）。
3.2 海底扇平面展布响应特征
对比3个时期的模拟结果：在Ps1层序限制性水道-朵叶复合体发育的地形地貌背景下，底流对重力流海底扇的平面展布特征影响较小，海底扇平面展布范围略微发生变化；在Ps2层序半限制性水道-朵叶复合体发育的地形地貌背景下，海底扇在有底流时沉积范围向北明显偏移、扩展；在Ps3层序非限制性朵叶复合体发育的地形地貌背景下，海底扇在有底流时沉积范围向北扩展更加明显。由此可见，底流作用会导致重力流海底扇平面展布特征发生变化，尽管海底地形地貌的限制性在一定程度上会削弱这种响应特征的表现，但是以泥质颗粒浓度表征的海底扇展布范围整体表现为顺底流方向偏移、扩展的特征。
4 重力流与底流交互作用沉积动力学机制分析
4.1 沉积机制
综上所述，相较于砂质颗粒，底流主要影响重力流中的泥质颗粒。泥质颗粒受到底流作用向北偏移，相比于无底流时的水道-朵叶复合体向北溢出并漂积，在水道北侧沉积充填，形成原水道北侧的泥质溢岸天然堤/漂积体。在重力流与底流持续地交互作用下，泥质颗粒不断向北偏移并沉积，从而影响了海底地形的演变。
4.2 侵蚀机制
在图6中，当无底流时，横向流速高速度核心位于水道轴部，侵蚀较强的区域位于水道轴部，表现为以垂向侵蚀为主（图6a）。当底流和重力流同时存在时，水道南侧出现矢量箭头较大的较高流速矢量，高速度核心靠近南侧水道壁，因此流体在垂向侵蚀的同时存在较明显的侧向侵蚀分量，对南侧水道壁造成向南的侧向侵蚀（图6b）。将无底流算例的剪切应力减去有底流算例的剪切应力，得出剪切应力差值图（图6c）。图6c中的正值区（黄色域—红色域）表示有底流作用情况下，该区域的剪切应力变小，侵蚀作用减弱；负值区（浅蓝色域—深蓝色域）则表示有底流作用情况下，该区域的剪切应力变大，侵蚀作用增强；而接近零值的区域（绿色域）表示有底流作用情况下，该区域的剪切应力及侵蚀作用大致不变。可以看出，在水道中部和北侧呈现一定面积的黄色域，说明在有底流时该区域剪切应力减小；在水道南侧呈现淡蓝色至蓝色域，说明在有底流时该区域剪切应力增加（图6c）。在模型最南端和最北端都呈现淡蓝色至蓝色域，这是来自模型南北两端设置的底流入水口提供的速度剪切力；模型东南角呈现红色域，则是由于该区域东南侧的地势较低并向南倾斜，底流受阻形成回流。总体表明：在重力流与底流交互作用下，重力流产生的高速度核心与自南向北的底流相结合，导致高速度核心和强剪切作用向南偏移，侵蚀区随之向南偏移。水道南侧受到比北侧更强的侧向侵蚀，削蚀逆底流一侧水道壁，水道向南迁移。
综上所述，底流对重力流沉积过程的影响具有显著的粒度选择性，主要表现为对细粒泥质颗粒沉积的显著控制作用，而对粗粒砂质颗粒的影响相对有限。底流作用导致泥质颗粒顺底流向北偏移、漂积，在水道北侧发育溢岸天然堤/漂积体，水道南北两侧不对称充填。无底流作用时，流体高速度核心位于水道轴部，以垂向侵蚀为主；有底流作用时，流体高速度核心和强剪切力向南偏移而靠近南侧水道壁，流体在垂向侵蚀的同时存在一定的侧向侵蚀，对南侧水道壁产生削蚀。在沉积-侵蚀机制的共同作用下，水道南侧持续被侵蚀而逆底流迁移，水道北侧持续充填顺底流漂积的细粒沉积物。在重力流与底流交互作用的持续影响下，侵蚀作用较强的高速度核心向南偏移，侵蚀南侧水道壁，北侧泥质颗粒溢出堆积，从而形成水道逆底流（向南）迁移现象。
5 结论
本文以鲁伍马盆地上始新统沉积为地质约束背景，基于Mixture模型和RNG k-ε模型，应用ANSYS软件模拟了鲁伍马盆地晚始新世重力流与底流交互作用的沉积动力学过程，揭示了逆底流侧向迁移水道的沉积-侵蚀机制。模拟结果表明：
（1）与粗粒砂质颗粒相比，重力流中的细粒泥质颗粒更容易受到底流作用的影响，从而顺底流方向偏移并沉积，在水道顺底流一侧沉积形成泥质溢岸天然堤/漂积体，顺底流偏移并漂积的细粒沉积物在水道南北两侧不对称充填的过程中形成不对称地貌。
（2）在重力流与底流交互作用下，流体高速度核心逆底流偏移，增强水道逆底流一侧的侵蚀作用，在泥质颗粒顺底流偏移、漂积形成不对称地貌的基础上，削蚀逆底流一侧的水道壁，促使水道逆底流方向迁移。
