波流作用下黄河三角洲硬壳层液化渗流形成机制研究*

冷 浩 胡瑞庚 刘红军② 王兆耀

(①中国海洋大学环境科学与工程学院， 青岛 266100， 中国) (②山东省海洋环境地质工程重点实验室， 青岛 266100， 中国)

0 引 言

黄河三角洲地区普遍发育着表面硬壳层，又称“铁板砂”，该表面硬壳层具有密度大、颗粒粗、渗透性大等特点，厚度一般1～3 m左右，在床面分布不连续。由于其上述特点，黄河三角洲表面硬壳层常作为海底管线和海洋平台的持力层，所以对硬壳层的形成机制以及其在波浪作用下的性质变化研究具有重要意义(贾永刚等， 2007)。

贾永刚等(2007)通过在黄河三角洲海床表面进行人工振动模拟波浪循环荷载以及室内模型试验，每轮荷载施加期间以及停止作用后的一段时间内海床表层一定范围内的土体发生液化，孔压消散后进行下一轮荷载，长期循环荷载作用下海床表层土体呈现重复液化，研究发现海床重复液化的结果导致了海床强度增强并且一定范围内的细颗粒物质向海床表面发生移动，海床表层土体密度增大。Yang et al. (2021)在黄河三角洲的新滩和广利港地区进行了原位静力触探试验，海床表层土体中的锥尖阻力qt和侧摩阻力fs数值较大(图 1)，表明该地区海床表层土体具有较高的土体强度，常方强等(2009)取黄河三角洲硬壳层存在区域原位1 m深度的土体进行土工试验，研究了硬壳层土体的粒度分布、孔隙比、超固结比等，砂粒、粉粒含量较高，砂粒含量约为16%，粉粒含量在44%～87%之间，而且该层土体相较于下伏土层处于明显的超固结状态，由此表明黄河三角洲硬壳层土体强度高，密度大，颗粒粗等特点。王虎等(2019)通过从原生沉积、冲刷分选和液化渗流3种模式对国内外铁板砂形成机制研究进行了总结分析，认为黄河三角洲铁板砂分布不连续的特点以及现场试验研究都表明液化渗流模式导致了黄河三角洲硬壳层的形成与发展。

图 1 原位CPT测量结果Fig. 1 In-situ CPT measured results

波浪作用下会在海床内部产生附加应力，循环附加应力会引起土体单元产生塑性应变，这种应变使得土体产生不可恢复的变形体积压缩，由于土体内部孔隙水压力来不及排出而不断累积，形成超孔隙水压力，当超孔隙水压力超过该点处的有效应力时土体发生完全液化(王虎， 2012； Hu et al.，2021)。

根据超孔隙水压力发展模式，波致海床液化可分为两种机制： (1)瞬态液化，在波谷位置，土骨架受到向上的渗流作用力，当向上的渗流作用力大于竖向有效应力时，发生瞬态液化(Sakai et al.，1992; Tsai， 1995;Jeng et al.，2015； 李正辉， 2019)； (2)累积液化，波浪作用下会在海床内部产生附加应力，循环附加应力会引起土体单元产生塑性应变，这种应变使得土体产生不可恢复的变形体积压缩，由于土体内部孔隙水压力来不及排出而不断累积，形成超孔隙水压力，当超孔隙水压力超过该点处的有效应力时土体发生完全液化(刘占阁， 2008； Sumer et al.，2012； 王虎， 2012； 廖晨聪， 2016； 段钊等， 2020； 刘晓磊等， 2020; Hu et al.，2021)。在海洋环境中存在海流作用时，波浪与海流耦合作用会引起海床表面波压力的变化，进而影响海床内部累积孔隙水压力的大小，顺流条件下会导致孔隙水压力累积幅值增大； 逆流条件则相反(刘小丽等， 2018)。本文以数值模拟方法研究波流联合作用下黄河三角洲海床内部孔隙水压力的累积液化特征。

海床表面为排水界面，由波浪作用而产生的超孔隙水压力通过排水边界消散，从而形成稳定的垂直于床面的渗流场，海床内部土颗粒受到重力、水流拖拽力、上举力以及渗流力的作用(图 2)。渗流力会显著降低泥沙启动的临界切应力，当土颗粒受到的切应力大于临界切应力时泥沙开始启动(Cheng et al.，1999； 程永舟等， 2012； 曹志刚等， 2019)。细颗粒物质受上述作用力影响明显首先被启动上浮，通过粗颗粒之间的孔隙被带到海床表面并逐渐发展成渗流通道，导致海床表层一定范围内的土颗粒粗化(图 3)。与此同时海床受到波致循环剪应力作用产生的塑性变形使海床表层土体密化，密化和粗化作用两者联合作用于海床从而形成了铁板砂。王虎等(2014)通过将渗流力考虑到传统泥沙启动公式中：推导并验证了波浪作用下海床临界冲刷深度的计算方法，有效解释了波浪对海床的改造作用。

图 2 泥沙颗粒受力情况Fig. 2 Stress of sediment particles

图 3 液化渗流模式下硬壳层的形成(王虎等， 2019)Fig. 3 Formation of the upper seabed layer under liquefaction seepage model(Wang et al.，2019)

基于土力学与泥沙运动学相结合，本文以黄河三角洲粉土海床为研究对象，采用数值模型研究波浪与海流同一时间共同作用下海床内部孔隙水压力累积特征，分析不同流速以及不同硬壳层厚度对累积孔隙水压力的影响。根据数值计算求得的累积孔隙水压力沿深度的分布值，基于王虎推导建立的海床临界冲刷深度计算方法，研究考虑渗流力作用的海床泥沙启动以及临界冲刷深度，分析不同波流条件对铁板砂形成深度的影响。

1 数值模型建立与验证

1.1 控制方程

利用数值模拟软件COMSOL来对波浪与海流同时刻共同作用于海床时的累积孔隙水压力响应进行计算。累积孔隙孔压响应计算模型控制方程和边界条件如下。

本文基于biot固结理论和渗流连续方程，同时考虑土骨架和流体的压缩性质，海床土体假定为多孔弹性介质，海床土体在x、z方向上的二维控制方程为

(1)

(2)

(3)

(4)

(5)

式中：u，w分别为土体的水平和竖向位移;εV为土单元体应变；pt为波致孔隙水压力；β为流体压缩系数；G为土体的压缩模量；υ为泊松比；k为渗透系数；γ为海床土重度；n为孔隙度；τ为土单元剪应力；γW为海水重度。

(6)

式中：K为孔隙流体体积模量；pW为绝对静水压力；Sr为饱和度。

海床累积孔隙水压力发展控制方程为

(7)

式中：pres为累积孔隙水压力；CV为土体固结系数；f为孔隙水压力累积源项。

采用已建立的适用于黄河三角洲海床的孔隙水压力累积源项f(胡瑞庚等， 2021)：

(8)

(9)

(10)

式中：N为动荷载循环次数；Nl土体达到液化时的循环次数；a，b，c为经验系数，根据动三轴试验分别得到a=0.326，b=0.8，c=0.75(刘红军等， 2005)。ar，βr为土体类型与密度相关参数，根据试验分别取0.221， 0.312(Sumer et al.，2012)。σ′0为平均有效应力，当累积孔隙水压力超过该值时海床即发生液化。

(11)

式中：K0为土体侧压力系数；γ′为土体浮容重；z为海床深度。

1.2 边界条件

对于波流共同作用在海床表面生成的波压力，采用hsu et al.(2009)提出的解析解：

(12)

式中：ρW为海水密度；H为波高；d为水深；λ为波数；U0为海流流速；ω，ω0，ω2分别为总圆频率、第一阶圆频率和第三阶圆频率， 表达式为:

(13)

(14)

ω=(λH)2ω2+ω0

(15)

海床表面处波致孔隙压力为波流共同作用产生的波压力，海床表面为排水界面，累积孔隙水压力为0，并且平均有效应力和剪应力消失：

(16)

海床两侧为周期性边界条件，上下两层海床交界面处为连续边界条件，即：

pres左=pres右;pres上=pres下

(17)

海床底部为不可压缩不可排水边界，即

(18)

1.3 数值模拟验证

通过Sumer et al.(2012)的试验数据来验证本文孔压数值模型的有效性，将试验的土体和波浪参数输入到本文模型中，计算结果与试验数据对比如图4所示，可以看出本文孔压数值模型计算结果与Sumer et al.(2012)的试验数据基本吻合，表明了本模型的可靠性。

图 4 累积孔隙水压力随时间的变化(z=-0.085 m)Fig. 4 Variation of residual pure water pressure with time

2 累积孔隙水压力响应分析

基于建立的数值模型，根据试验以及相关资料获得了黄河三角洲地区原始海床以及表面硬壳层的物理力学指标，波浪条件为黄河三角洲地区一年一遇海况，并选取不同的海流流速，研究波浪共同作用下不同流速对累积孔隙水压力分布的影响。波流参数与海床土体参数见表 1。

表 1 波流参数与海床土体参数取值Table 1 Wave-current parameters and seabed soil parameters

2.1 原始单层海床累积孔压分布特征

考虑未形成表面硬壳层的原始单层海床，观察不同波浪作用时间累积孔隙水压力的变化以及液化情况。

由图 5 表明随着波浪作用时间的持续增长，海床累积孔隙水压力逐渐增加，但累积孔压每隔30 min的增量随着时间而减小，这是由于孔压的消散以及附加应力导致的海床孔隙压缩极限无法再产生塑性变形，最终稳定后的孔隙水压力称为稳态超孔压力。由于真实海况波浪持续时间以及数值模型计算量的原因，本文只考虑波浪作用120 min内的累积孔隙水压力结果。由图 5 可知，累积孔压沿深度分布存在极大值，该极值随波浪持续作用时间逐渐下移，并且由图 5 还可以看出累积孔隙水压力沿深度分布曲线与平均有效应力线的交点随时间逐渐下移，海床液化深度范围逐渐增大，海床的液化深度在波浪持续作用90 min的时间间隔内，液化深度由0.62 m增加至2.14 m。

图 5 波流作用下累积孔隙水压力沿深度分布随时间的变化Fig. 5 Variation of residual pore water distribution along depth with time under wave and current action

海流与波浪联合作用会显著影响海床累积孔隙水压力的分布，取表中的不同海流流速参数，分析不同流速对海床累积孔压分布特征的影响规律。图 6为波流持续作用120 min时，不同流速下海床孔隙水压力沿深度的分布图。

由图 6 可知，随着流速的增大，海床累积孔压相应增大，流速由0 m · s-1增大至2 m · s-1时，海床累积孔隙水压力极值由12.15 kPa增加到29.63 kPa。但流速为负值，即流速流向与波浪行进方向相反时，海床孔隙水压力则相应减小，流速为- 2 m · s-1时的累积孔压与流速为0 m · s-1相比仅下降了约40%。所以，当流速流向与波浪行进方向一致时，海流会对海床累积孔隙水压力的增长起促进作用，显著增加了海床累积孔隙水压力，而流向相反时，会对海床累积孔隙水压力起到抑制作用，这种抑制作用相对于流速一致时累积孔压的变化，影响较小。

图 6 波流联合作用不同流速下累积孔隙水压力沿深度分布Fig. 6 Vertical distribution of residual pore water pressure for different velocity under wave and current action

2.2 硬壳层厚度对海床累积孔压的影响

黄河三角洲地区分布着厚度不一的硬壳层，选取表 1 中的硬壳层相关物理力学指标，并分别取硬壳层厚度为1 m， 2 m， 3 m，研究不同硬壳层厚度对海床累积孔压的影响规律。图 7为波流持续作用120 min时，不同硬壳层厚度下海床累积孔隙水压力沿深度的分布图。

图 7 波流联合作用不同硬壳层厚度对累积孔隙水压力的影响(t=120 min)Fig. 7 Effect of hard shell thickness on cumulative pore water pressure under wave and current actiona. 硬壳层厚度1 m； b. 硬壳层厚度2 m； c. 硬壳层厚度3 m

由图 7可以看出硬壳层存在时，海床累积孔隙水压力显著减小，随着硬壳层厚度的增加而相应减小，累积孔隙水压力沿深度分布出现转折，在原始海床中沿深度增加迅速。这表明密度大、孔隙大的表面硬壳层有利于海床累积孔压的消散。海床累积孔隙水压力沿深度分布的极值均出现在下层原始单层海床中，流速为0 m · s-1时硬壳层厚度由1 m增加到3 m，极值点深度下降了1.38 m。

由图 7还可以看出累积孔压沿深度分布极值所处深度均靠近两层海床界面处，硬壳层厚度越小越接近，而且该极值大小随着厚度的增加也约接近海床界面处累积孔隙水压力值。这是由于波流作用于海床影响显著的区域大约在1～3 m之间(图 6)，然而硬壳层的存在有利于海床孔隙水压力的消散，因此随着硬壳层厚度的增加，极值的增量也相应减小。

3 渗流模式下硬壳层形成机制

波浪及潮流共同作用于海床表面，会使海床表层附近的细粒泥沙颗粒发生启动，造成表层颗粒较粗、密度较大的情况。对于厚度较大的硬壳层是由形成于海床内部的渗流力导致的，由上文可知在波浪以及潮流作用下，海床内部会产生较大的累积孔隙水压力，这些累积孔隙水压力会逐渐向外排出，海床相邻点处的累积孔隙水压力差值形成的水力梯度产生渗流力，进而引起水流动，因为海床表面为排水界面，从而形成向上的渗流。海床内部的细颗粒泥沙受到渗流力的作用打破受力平衡状态产生位移，继而沿着粗颗粒之间的孔隙运移至海床表层。海床内部累积孔隙水压力随着时间逐渐增大，因此在海床深度较大的位置也可产生指向床面的渗流力，形成较厚的粗颗粒层。累积孔隙水压力的形成是由于波浪作用于海床，在海床内部产生循环附加应力，使海床内部颗粒之间的孔隙压缩产生超孔隙水压力，来不及排除以至逐渐累积，压缩孔隙体积使海床变得致密，细颗粒向上运移和孔隙压缩两方面同时进行形成了硬壳层海床。

根据王虎等(2014)的研究， 当波流引起的海床表面最大切应力超过临界启动切应力时，土体颗粒启动，进而形成冲刷。很多学者在传统泥沙颗粒力学平衡方程中考虑了渗流力的作用，并建立了新的力学平衡方程，并且窦国仁等(2001)给出了该力学平衡方程中的重力FG、水流拖拽力FD、上举力FL以及渗流力FS的表达式。根据该新的力学平衡方程以及各个受力的表达式，王虎等(2014)推导建立了海床临界冲刷深度计算公式，研究渗流力对海床冲刷的影响。推导过程如下：

泥沙颗粒受力平衡方程(图 2) 表达式为：

FD=(FG-FL-FS)tanφ

(19)

式中：FD为水流拖拽力；φ为饱和土静止内摩擦角。其中，FG，FL，FD，FS表达式为(窦国仁等， 2001)：

(20)

(21)

(22)

(23)

式中：τs为临界启动切应力；d50为中值粒径； 阻力系数CD=0.4； 上举力系数CL=0.1； ΔP为相距ΔL两点间的累积孔隙水压力差值。

根据式(17)～式(21)，临界启动切应力τs表达式为：

(24)

通过数值计算模型求得的海床土体累积孔隙水压力值，进而可得到海床任意深度处的渗流压力梯度，然后根据式(22)求解该深度处土体颗粒的临界启动切应力。

波流作用引起的海床表面最大切应力(钱宁等， 2003)为：

(25)

式中：ρ1为床面附近流体密度；um为海底最大轨迹质点速度；fW为琼森阻力系数，决定于水流雷诺数Re和海床相对粗糙度2am/d50，两参数取值表达式如下：

(26)

(27)

取得Re和2am/d50值，根据Sassa et al.(2001)文献可得fW值，最终，可由式(25)求得海床表面最大切应力值。

床面最大切应力τm大于临界启动切应力τs时，泥沙可发生启动。细颗粒泥沙首先发生启动，通过粗颗粒之间形成的孔隙运移至海床表面。当τm小于τs时，则不会发生冲刷。相关计算结果见表 2。

表 2 相关计算参数及结果Table 2 Relevant calculation parameters and results

由表 2 可知，当流速U0=1 m·s-1时海床深度为0 m处的临界启动切应力为负值，流速U0=2 m · s-1时临界启动切应力负值已达到海床深度0.5 m处，这说明该两种流速下引起的ΔP/ΔL值较大，以致(ρ-ρW)g-ΔP/ΔL值为负，τs为负值。流速流向与波浪行进方向一致时可形成较大的渗流力作用于泥沙颗粒上，渗流力越大可以将较粗的泥沙颗粒运移至海床表面，形成更粗的硬壳层。并且由表 2 还可以看出随着流速的增大，泥沙颗粒的临界启动切应力也相应减少，这表明流速越大，泥沙颗粒发生启动所需要的海床表面最大切应力也越小，泥沙越容易启动。

图 8 不同流速下切应力沿深度变化Fig. 8 Shear stress changes with depth at different velocitya. 流速U0=0 m·s-1； b. 流速U0=1 m·s-1

由图 8 可知，τs随着海床深度的增加急剧下降，τm随着深度增加而增加，两值随深度快速接近，两条线相交处的地方泥沙启动停止。从图 8 也可以看出两条线并没有相交，但是这两个不同流速下的泥沙启动在1.5 m深度处停止，这是由于波浪引起的累积孔隙水压力沿深度分布极值在1.5 m深度附近，并且流速存在对于该极值的深度分布贡献也极小(图 6)，所以海床泥沙启动只发展到1.5 m深度左右。该深度符合黄河三角洲地区硬壳层厚度普遍为1～3 m的情况。黄河三角洲地区在遭受5年一遇以及50年一遇的复杂海况条件时，累积孔隙水压力沿深度分布极值会分布在海床更深的位置，并且引起较大的ΔP/ΔL值，从而形成较厚以及颗粒较粗的表面硬壳层。

4 结 论

(1)波浪引起的累积孔隙水压力与波浪持续时间成正相关，累积孔隙水压力随时间逐渐增加，但增速逐渐放缓，最终会达到稳定状态，海床孔隙不可再压缩。

(2)海流对海床累积孔隙水压力有重要影响。海流流向与波浪行进方向一致时，对累积孔隙水压力起促进作用，流速越大累积孔隙水压力越大，反之对累积孔隙水压力有抑制作用。

(3)表面硬壳层的存在会显著促进累积孔压的消散，累积孔隙水压力沿深度分布的极值均出现在下层原始海床中，流速U0=0 m · s-1时硬壳层厚度由1 m增加到3 m，极值点深度下降了1.38 m。

(4)累积孔隙水压力引起的渗流力对于海床泥沙启动影响显著，在流速U0=0 m · s-1，U0=1 m · s-1时泥沙启动深度均为海床1.5 m深度处，并且海流流向与波浪行进方向一致时，会产生较大ΔP/ΔL值带动较粗的泥沙颗粒至海床表层，但对泥沙启动的最大深度影响不大。