波浪作用下考虑海床变形影响的溶质运移数值模拟

刘小丽,刘明珠

(1. 海洋环境与生态教育部重点实验室,山东青岛 266100；2. 中国海洋大学环境科学与工程学院,山东青岛 266100)

海底沉积物是海洋生态系统的重要组成部分,海水与沉积物之间的物质交换对海洋生态环境和海洋生物地球化学过程有重要影响[1-2]。自20世纪70年代起,较多学者以溶质运移理论为基础,研究了波浪水动力对海水中溶质进入海底沉积物运移过程的影响。

溶质在多孔介质中的运移过程常采用对流-弥散方程来描述[3-4]。基于一维对流-弥散方程,Harrison等[5]和Clark等[6]对行进波作用下沉积物中的溶质运移进行了计算分析,研究表明,对于平底海床,波浪循环荷载主要通过机械弥散作用加速溶质在海底沉积物中的运移,可使海水中溶质进入海底沉积物中的速率比单纯的分子扩散作用增加几个数量级。由于波浪荷载作用下流入和流出海底沉积物中的水体交换量基本相当,因而在波浪荷载作用较为显著的浅海区域,溶质运移过程中的对流作用基本可以忽略[5,7]。Qian等[1,8-9]基于二维对流-弥散方程,对波浪作用下粒径为0.001～0.1 m的粗砂和砾石沉积物中的溶质运移过程进行了一系列模拟计算,通过引入“竖向增强弥散系数”的概念,提出了将二维溶质运移问题转化为一维问题进行分析的方法；同时,研究了周期性波浪压力下将二维各向异性弥散张量简化为各向同性弥散张量的条件；结果表明,波浪水动力会加快沉积物中的溶质运移。虽然较多文献对波浪作用下沉积物中的溶质运移过程进行了研究,但其计算分析中均未考虑波浪导致的沉积物变形效应的影响。Biot固结理论可以考虑荷载作用下多孔介质变形场和孔隙水压力场的双向耦合效应[10-11]。Yamamoto[10]、Jeng和Hsu[11]均基于弹性介质Biot固结理论,提出了线性波浪荷载作用下海床瞬态响应的解析解,通过与波浪作用下砂质海床的试验数据对比,验证了解析解的可靠性和Biot固结理论的适用性,并进一步分析了不同计算参数下砂质海床的应力、孔隙水压力和变形响应特征;研究结果表明,波浪会导致砂质海床(沉积物)发生一定程度的变形,并影响沉积物中孔隙水的渗流过程。根据溶质运移理论,孔隙水的渗流与溶质在沉积物中的运移过程密切相关,因此,有必要研究波浪导致的海床土变形效应对溶质运移过程的影响,目前这方面的研究工作尚未开展。

为揭示波浪作用下海床变形对溶质运移过程的影响规律和作用机制,以砂质海床为研究对象,基于多孔弹性介质Biot固结理论[12]及多孔介质中溶质的对流-弥散方程[3-4],构建波浪作用下考虑海床土体变形影响的溶质运移数值计算模型,在验证数值模型正确性的基础上,研究线性行进波作用下海底沉积物变形对非反应性溶质运移过程的影响及其作用机制,并进一步分析海床土饱和度对溶质运移过程的影响,为海洋环境中的物质交换和生态健康评估提供参考。

1 数值模型

图1为一阶线性行进波作用下海底沉积物中孔隙水渗流和溶质运移计算的几何模型示意图。平底海床厚度为h,坐标轴原点位于海床的左上角,x轴位于海床上表面,y轴以向上为正。线性行进波的波长为Lw,波高为Hw,周期为Tw,水深为dw。

图1 海床孔隙水渗流和溶质运移计算的几何模型示意Fig.1 Geometric model for simulation of pore fluid seepage and solute transport in seabed

参照文献[1,8-9],设砂质海底沉积物为均匀各向同性体,沉积物-海水交界面处的溶质浓度c0=1 mg/L,且随时间保持不变,砂质海底沉积物中的初始溶质浓度为cs0=0。需要说明的是,c0的具体数值并不会影响对溶质运移规律的分析。

1.1 控制方程

1.1.1 海底沉积物中孔隙水渗流

假设海床土骨架为线弹性体,应力-应变关系服从广义胡克定律,根据弹性多孔介质Biot固结理论,海床土体的二维平衡方程为[12]：

(1)

(2)

式中:μ为土体泊松比；G为土体剪切模量；p为波致海床孔隙水压力；wx、wy分别为海床土在x和y方向的位移；ε为土骨架体积应变。

孔隙水流动服从达西定律,根据流体质量守恒,孔隙水渗流连续方程为[12-13]

(3)

式中:γw为海水的容重；n0为海床土孔隙率；K为海床土渗透系数；t为时间；土骨架体积应变ε的表达式为

(4)

β为孔隙流体压缩系数,考虑溶解于或以微小气泡封闭于孔隙水中的气体,其表达式可取为[11,13]

(5)

式中:Kw0为孔隙水体积压缩模量,2×109N/m2；Pw0为绝对静水压力；Sr为海床土的饱和度。

根据达西定律,多孔介质中孔隙水的实际渗流速度为

(6)

式中:u和v分别为沉积物中孔隙水的水平向和竖向实际渗流速度。

式(1)—式(3)构成Biot固结理论的控制方程,可以看出,Biot固结方程考虑了土骨架变形与孔隙水渗流的双向耦合作用,基于Biot固结方程计算得到的海床孔隙水压力和渗流场反映了土体的变形效应。根据式(5)可知,土体饱和度直接影响孔隙流体的压缩性,孔隙水和气被视为同一种可压缩流体考虑,饱和度越小,孔隙中的气体越多,压缩性越大。

当不考虑沉积物变形效应时,渗流连续方程可以由式(3)退化得到,即此时的渗流连续方程中土骨架的体积应变为0,孔隙流体不可压缩,具体表达为[14]

(7)

式中:φ为孔隙水的压力水头；p=γwφ。

1.1.2 海底沉积物中溶质运移

沉积物中溶质的运移遵循质量守恒定律,非反应性溶质的二维对流-弥散方程表达式为[3-4]：

(8)

综上,联立Biot固结方程和溶质对流-弥散方程,结合适当的边界条件,可以得到线性行进波作用下考虑土体变形影响的溶质运移计算模型,该模型属于沉积物变形与溶质运移的单向耦合模型,即沉积物变形影响溶质运移,但溶质运移不影响沉积物变形。

1.2 边界条件

结合图1所示,在一阶线性行进波作用下,海床渗流场计算的边界条件为：

(2) 海床底部(y=-h)为刚性不透水边界,孔压竖向梯度和土体位移均为0。

(3) 海床两侧(x=0和Lw)为周期性边界条件：p|x=0=p|x=Lw；wx|x=0=wx|x=Lw；wy|x=0=wy|x=Lw。

溶质运移计算的边界条件为：海水与沉积物接触面的溶质浓度保持为c0=1 mg/L；沉积物底部不透水,溶质的竖向通量为0,∂c(x,-h,t)/∂y=0；在沉积物计算区域的两侧,溶质浓度保持为周期性对称边界,c|x=0=c|x=Lw,沉积物孔隙水中的初始溶质浓度为0。

1.3 模型的求解

利用上述边界条件,同时求解式(1)—式(3)和式(8),可以得到考虑沉积物变形影响的海床位移场、孔压场、流速场以及溶质浓度场；联立求解式(7)和式(8),可以得到未考虑沉积物变形影响的波致孔隙水渗流场和溶质浓度场。

上述控制微分方程组的求解,借助于Comsol Multiphysics®数值计算平台完成[17]。Comsol Multiphysics®是一个可以进行多物理场耦合分析的有限元数值计算平台,利用其提供的偏微分方程模式(PDE),用户可以将需要求解的偏微分方程组,输入进行求解。首先基于该平台,建立所求解问题的几何模型,进行网格剖分；然后在PDE模式下输入需要求解的控制微分方程组,设置适当的边界条件,提交至平台利用有限元法对方程组进行求解。数值计算中采用二阶拉格朗日单元进行网格划分,所划分单元的尺寸以对计算结果无影响为原则,具体的单元数量和大小详见各算例说明。

2 数值模型验证

分别将波浪作用下沉积物中溶质运移和孔隙水压力的数值计算结果与试验结果进行对比分析,以验证数值模型计算结果的正确性。

2.1 溶质运移计算结果验证

Packman等[18]对上覆水中溶质进入砾石底床的过程进行了一系列水槽试验,利用一定流速的水流流经波浪状沉积底床的方式,研究不同波动压力条件下上覆水中溶质向底床的运移情况,给出了上覆水中溶质浓度随时间的变化曲线。基于试验砾石底床渗透性大和基本不变形的特点,同时求解式(7)和式(8),计算从上覆水进入到沉积底床中的溶质通量以及底床中的溶质浓度分布。

采用与试验中完全相同的参数,选取2个试验工况进行计算对比,试验砾石底床的参数为：颗粒平均粒径0.6 cm,孔隙率0.38,渗透系数0.15 m/s,沉积物厚度0.19 m。根据试验流速水头和波浪状底床参数,利用文献中公式计算得到床面上的波动压力。试验工况1和工况2中底床表面的波压力表达式分别为：pb1=6.67 sin(19.63x),pb2=2.91 sin(19.63x),波压力的单位Pa；x为水平坐标,单位m；波压力对应的波长均为0.32 m。有限元数值计算模型中共划分2 579个单元,最大单元边长0.005 m。

为了与试验结果对比,根据Packman等[18]研究成果,通过进入沉积物中的溶质量计算得到上覆水中溶质浓度cw随时间的变化曲线,结果如图2所示,为上覆水中溶质平均浓度随时间变化的数值计算结果和试验结果的对比图,可以看出,总体而言,数值计算结果与试验结果一致性较好,表明了溶质运移模拟计算结果的可靠性。

2.2 波致海床孔隙水压力计算结果验证

Qi等[19]通过水槽试验模拟了波浪和水流作用下砂质沉积物中的孔隙水压力响应。选取波浪作用下的试验工况进行计算,采用与试验中完全相同的计算参数：砂土底床的孔隙率n0=0.435,泊松比μ=0.3,饱和度Sr=0.995,平均粒径dg=0.38 mm,渗透系数K=1.88×10-4m/s,剪切模量G=1×107Pa,底床厚度h=0.5 m；试验水深dw=0.5 m,波浪Hw=0.095 m,周期Tw=1.2 s。

试验底床为中砂,波浪作用下存在一定的变形效应。为说明海床土变形效应对孔隙水压力的影响,利用数值模型分别计算了考虑和不考虑海床土变形时的孔压响应,并将本文数值计算结果与试验结果、线性波下的解析解进行对比,以验证数值计算结果的正确性。有限元计算模型中共划分1 622个单元,最大单元边长0.04 m。线性波浪作用下有限厚海床孔隙水压力解析解的计算公式来源于文献[11],该解析解是基于Biot固结方程推导而来,考虑了海床土变形的影响。计算对比结果如图3所示。

图3中纵坐标为量纲—深度ky,k为波数。从图中可以看出,考虑海床土变形影响的数值解与解析解吻合较好,但均与试验值在海床上部存在一定偏差,推测是试验中波浪的非线性以及部分海床土计算参数与实际参数之间存在些许偏差等原因所致。对比考虑海床土变形与不考虑海床土变形的孔压计算结果可知,考虑海床土的变形效应,计算得到的海床孔隙水压力更符合实际情况。

图2 数值计算结果与试验结果对比Fig.2Comparison of numerical and experimental results

图3 数值计算结果与解析解和试验结果的对比Fig.3Comparison of numerical,analytical and experimental results

3 结果与讨论

考虑砂质海床土体的变形效应,对波浪作用下上覆水中的溶质进入砂质海底沉积物中的运移过程进行数值模拟,分析海床土的变形效应对砂质沉积物中溶质运移过程的影响机制。波浪参数为：波高Hw=5 m,周期Tw=10 s,水深dw=20 m,波长Lw=121.2 m；海底沉积物的计算参数如表1中所示,沉积物-海水交界面处的溶质浓度为c0=1 mg/L。

根据文献[10-11],线弹性砂质海床瞬态波浪响应分析中,剪切模量取值范围为106～109Pa。有限元计算模型中共划分9 286个单元,海床面附近单元尺寸较小,边长约0.02 m,海床底部单元尺寸较大,边长约1 m。

表1 海底沉积物计算参数Table 1Parameters of marine sediments

3.1 波浪作用下海床变形对溶质运移的影响机制

如前所述,周期性波浪荷载作用下,可以忽略平底海床溶质运移过程中的对流作用,溶质主要通过机械弥散作用运移扩散。根据文献[1],结合计算结果分析可知,行进波作用下,沉积物沿水平方向均承受相同的周期性波浪压力,同一深度水平方向上各点的孔隙水压力响应相同,其溶质浓度沿水平方向是均匀分布的,溶质主要沿沉积物的深度方向(竖向)运移。因此,主要对行进波作用下砂质沉积物中孔隙水压力的竖向分布、孔隙水的竖向(y方向)流速、溶质在沉积物中的竖向(纵向)水动力弥散系数以及溶质沿沉积物深度方向的运移特征进行分析。

土体的剪切模量是衡量土体变形能力的一个重要指标[10]。G越大,相同荷载下土体体积应变越小,土体抵抗变形的能力就越强。对于弹性多孔海床,在周期性波浪荷载作用下,海床土的位移和变形随着波浪荷载的变化而呈现周期性变化,波峰和波谷的交替作用,伴随着孔隙水的流入和流出。海床土的这种周期性瞬态变形,会影响沉积物中孔隙水的渗流[10-11]。不同剪切模量条件下海床土骨架的体积应变幅值、孔隙水压力幅值和孔隙水竖向流速幅值沿深度的分布如图4所示。

图4 不同剪切模量条件下土骨架体积应变幅值、孔隙水压力幅值、孔隙水竖向流速幅值随深度的分布Fig.4 Vertical distributions of amplitude of volumetric strain of soil skeleton,pore water pressure and vertical velocity of pore water for various shear modulus

图4(a)为海床土骨架体积应变幅值随海床深度的分布曲线,不考虑土体变形效应时的土体体积应变为0。从图中可以看出,随着土体剪切模量的降低,在距海床表面约1.5 m深度范围内,土骨架的体积应变呈明显的非线性增加趋势。土骨架体积应变越大,单位时间内流经土体单元的水量变化就越大,孔隙水的流速就越大,根据达西定律可知其对应的孔隙水压力竖向梯度也越大,即表现为孔压随深度的衰减越快。如图4(b) 中所示,孔压在深度方向衰减速率的增加,对应于该深度内孔隙水竖向流速的增大。如图4(c)中所示,当剪切模量为106Pa时,海床表面孔隙水竖向流速幅值约为不考虑土体变形时的6.5倍。

从图4中还可以看出,随着土体剪切模量的增加,土体体积应变量减小,孔隙水压力和孔隙水流速均逐渐趋近于不考虑土体变形时的结果；当土体的剪切模量增加至109Pa时,海床土体的变形效应对孔隙水压力和孔隙水流速的影响已经很小。

图5为溶质纵向水动力弥散系数幅值以及波浪作用1 800个周期时溶质浓度和溶质通量随深度的变化曲线。从图5(a)可知,溶质纵向弥散系数在海床表面处的数值最大,随深度逐渐减小。土体变形效应对溶质纵向弥散系数的影响显著,随着海床土骨架体积应变量的增加(即随着剪切模量的降低),孔隙水的流速变大,导致距海床表面1.5 m深度内的溶质纵向弥散系数明显增大,如当G=106Pa时,海床表面纵向弥散系数的幅值约为不考虑土体变形效应时的8.5倍,约为分子扩散系数的545倍。

图5 沉积物中溶质纵向水动力弥散系数、溶质浓度和溶质通量随深度的分布(t=1 800Tw)Fig.5Vertical distribution of vertical hydrodynamic dispersion coefficient of solute,solute concentration and solute flux (t=1 800Tw)

溶质纵向水动力弥散系数的增大,加速了沉积物中溶质的运移。从图5(b)中可以明显看出,随着土体变形效应的增大,溶质运移的速率增大,运移深度增加。不考虑土体变形时,溶质运移深度约为12 cm；考虑土体变形,当G=106Pa时,溶质运移深度约30 cm,是不考虑变形效应时的2.5倍,且随着计算时间的增加,这种差异会不断加大。同时还可以看出,当G=109Pa时,由于土体的变形已经很小,其对溶质运移过程基本无影响。

图5(c)所示为沉积物上部35 cm深度内的溶质通量分布曲线,由于周期性波浪荷载作用下平底海床中的对流作用基本可以忽略[5,7],因此其溶质通量即为溶质的弥散通量,不同深度处的平均溶质通量可通过以下公式进行计算：

(9)

从图5(c)中可以看出,海床剪切模量越小,溶质弥散通量相对越大,总体而言,溶质弥散通量随深度的增加而减小。当土体剪切模量G=106Pa时,从沉积物表面至5 cm深度的范围内,弥散通量存在一个先迅速降低然后再增加到极大值的现象,这与该深度范围内竖向溶质浓度梯度沿深度的分布有关。较大的土体变形会促进溶质的运移,在距离沉积物表面一定深度内(此处为0～5 cm)的溶质浓度均较大,由此导致竖向浓度梯度最大值的位置向更深处移动;由于纵向弥散系数随深度逐渐减小,根据二者乘积计算得到的溶质通量,在距离沉积物表面一定深度处会出现极大值,且该极大值的位置与竖向浓度梯度最大值的位置相对应。相比较而言,当土体变形效应较小时,竖向溶质浓度梯度在距沉积物表面非常小的深度处即达到最大值,对应的溶质弥散通量表现为在沉积物表面处最大、随深度不断减小的变化趋势。

根据上述分析可知,海床剪切模量越小,土体变形效应越明显；较大的土体变形会增大孔隙水的流速,促进溶质的运移,增加溶质运移深度,并导致沉积物表面以下一定深度处的溶质弥散通量出现极大值,该溶质通量极大值出现的深度会随计算时间和土体变形量的增大而有所增加。

3.2 海床土饱和度对溶质运移的影响

海底土体中常常因存在少量的气体而处于非饱和状态[20]。气体的存在会影响孔隙流体的压缩性,土体饱和度越低,孔隙流体的压缩性越大,对孔隙水渗流和溶质运移的影响也越大。对不同饱和度条件下波致海床渗流场和溶质浓度场进行计算,分析海底土体的饱和度对溶质运移的影响。海床的剪切模量取为G=106Pa,饱和度Sr的值分别取为0.94,0.96,0.98和1。

如图6所示,分别为不同饱和度下海床孔隙水压力幅值、孔隙水竖向流速幅值和溶质浓度随深度的变化曲线。从图6(a)可知,随着饱和度的降低,孔隙流体压缩性增大,导致海床上部约2 m深度内孔压随深度的衰减速率随饱和度降低而增大,孔隙水流速也相应增加(如图6(b)所示),孔隙水的竖向流速随饱和度的降低而增加。孔隙水流速的增大促进了溶质水动力弥散系数的增大,加速了沉积物中的溶质运移(如图6(c)

图6 不同饱和度下沉积物孔隙水压力幅值、孔隙水竖向流速幅值和溶质浓度随深度的分布(G=106 Pa)Fig.6 Vertical distributions of amplitude of pore pressure,velocity of pore water and solute concentration for various degree of saturation (G=106 Pa)

所示),饱和度越小,溶质的运移深度越大。在此算例中,Sr=0.94时,溶质运移深度约为42 cm,是饱和海床溶质运移深度的1.4倍。

4 结 论

以平底海床为例,构建了波浪作用下考虑海底沉积物变形效应的溶质运移数值计算模型,在模型可靠性验证的基础上,对波浪作用下溶质进入砂质海底沉积物的运移过程进行了数值模拟,主要得到以下结论：

(1) 波浪作用下砂质海床土体的变形效应,会加速非反应性溶质在其中的运移速率,增大其运移深度。其作用机制在于,土体的变形会使沉积物上部一定深度内的孔隙水流速增加,进而增大溶质纵向水动力弥散系数,增强溶质扩散过程中的机械弥散作用。海床的剪切模量越小,变形效应的影响越显著。算例分析表明,海床土剪切模量G=106Pa时,溶质最大纵向水动力弥散系数约为不考虑变形效应时的8.5倍,约为分子扩散系数的545倍。

(2) 波浪作用下随着砂质沉积物中溶质浓度的增大,会导致溶质通量在距海床面一定深度的位置出现极大值,该极大值出现的深度随溶质运移时间和土体变形效应的增大而有所增加,且与溶质竖向浓度梯度最大值的位置相同。

(3) 海床土饱和度降低,孔隙流体的压缩性增大,从而导致孔隙水流速的增加,并进一步加速波浪作用下沉积物中溶质的运移。