内孤立波是一种自然界中常见的非线性内波^[1], 是发生在层结海洋内的一种扰动^[2], 其最大振幅出现在海洋内部^[3]。内孤立波具有周期短、振幅大、流速高、剪切强的特性^[4-5]。其引起的等密度面起伏形成“水下断崖”威胁水下航行器, 流动引发的强剪切影响海上平台^[6]。此外, 内孤立波引起物质输运与混合, 沉积物再悬浮, 对海洋地质、海洋生物有一定影响^[7-9]。

内孤立波多发于全球各大洋的陆架边缘海^[10], 其中南海北部是非线性内波最为活跃的区域之一, 国内外学者针对此海区的内孤立波开展了大量现场观测与遥感观测, 如亚洲海联合声学实验(Asian Seas International Acoustics Experiment, 简称ASIAEX)、非线性内波计划(Nonlinear Internal Waves Initiative, 简称NLIWI)、海峡内波实验(Internal Waves in Straits Experiment, 简称IWISE)等观测计划^{[2, 3, 11-14]}, 揭示了内孤立波由深水区向浅水区的传播演变过程。受地形与层结的影响, 内孤立波在传播过程可分为不同阶段: 第一阶段, 在深水区内孤立波维持波形传播。第二阶段, 随水深的变浅, 内孤立波会发生变形和裂变, 波前变平, 波后变陡, 波形变宽, 能量辐射给尾波, 发展成为波列^[15]。第三阶段, 内孤立波到达临界深度时, 发生由下凹型到上凸型的极性转化^[16-18]。第四阶段, 内孤立波在陆架海中因不稳定性引发破碎, 引起海洋混合与能量耗散^[19-20]。

针对内孤立波向陆架海传播这一重要演化过程, 理论模型与数值模拟相比于传统观测, 可完整刻画内波的传播演化过程, 是研究内孤立波动力机制的有效手段。Korteweg等^[21]发展了经典弱非线性模型KdV方程, 并被广泛应用于描述内孤立波这类波动。在非线性与频散作用相平衡时, 内孤立波会维持波形传播。但实验室实验与现场观测结果均证实KdV方程局限于弱非线性假定, 只适用于小振幅内孤立波, 在刻画大振幅内孤立波波形时与实际波形差异大^[22-24]。在临界深度附近, 平方非线性系数趋近零无法有效刻画内孤立波动力演化过程。因此描述大振幅内孤立波时需要引入高阶非线性修正项, 如扩展KdV(extended KdV, 简称eKdV)方程, 修正KdV(modified KdV, 简称mKdV)方程与Michallet-Camassa and Choi(MCC)方程等^{[22, 25-26]}。基于弱非线性理论模型, 前人研究了内孤立波在陆架地形上的浅水演化过程, 模拟了极性转化过程^[17]、裂变过程^{[17, 27]}、探讨了演化过程中的动力机制, 如地转影响^[28]、波形绝热演化^[29-30]、特征参数分布等^[31-32], 指出地形与层结是影响内孤立波传播演化的关键因素。

前人研究多集中在内孤立波过陆架陆坡地形的波形演化, 对弱非线性模型中的动力机制需要进一步讨论。为后续开展实验室实验作验证, 本文使用实验室尺度的弱非线性规范化长波(regularized long wave, 简称RLW)模型^[31], 同时对比于全非线性的MIT General Circulation Model(MITgcm)模型^[33]模拟不同入射波振幅、地形高度条件下内孤立波的传播演化过程, 分析演化中的动力机制, 比较了两模型的演化过程差异, 评估弱非线性模型的适用性, 并应用于南海讨论实际海洋情况。

1 模型设置与数据分析 1.1 数值模拟设置

Cai等^[31]建立的RLW模型在eKdV基础上考虑内孤立波实际传播过程, 增加耗散项与浅滩效应项, 控制方程为:

$ {\eta _t} + {c_0}{\eta _x} + \alpha \eta {\eta _x} + {\alpha _1}{\eta ^2}{\eta _x} + \beta {\eta _{xxx}} - \frac{1}{2}A{\eta _{xx}} + \gamma \eta = 0, $

(1)

其中, η为内波引起的跃层起伏, c₀为线性相速度, α为平方非线性系数, α₁为立方非线性系数, β为频散系数, A为耗散系数, γ为浅滩效应系数。各项动力系数由下式求出:

$ {c_0} = \sqrt {\frac{{g'{h_1}{h_2}}}{{{h_1} + {h_2}}}} , $

(2)

$ \alpha = \frac{{3{c_0}}}{2}\left( {\frac{1}{{{h_2}}} - \frac{1}{{{h_1}}}} \right), $

(3)

$ {\alpha _1} = - \frac{3}{8}\left( {\frac{1}{{h_1^2}} + \frac{1}{{h_2^2}} + \frac{6}{{{h_1}{h_2}}}} \right), $

(4)

$ \beta = \frac{{{c_0}{h_1}{h_2}}}{6}, $

(5)

$ \gamma = {\gamma _0}\frac{{{h_x}}}{{{{\left| {{h_x}} \right|}_{\max }}}}, $

(6)

其中, h₁、h₂分别为水体上层厚度与下层厚度, g'为约化重力加速度, h_x为地形水平梯度。

弱非线性模型设置参考实验室条件, 如图 1所示, 总长度L为10 m, 总水深D为0.3 m。模型水平分辨率为0.01 m, 对应网格数为10 000, 时间步长设置为0.001 s, 计算100 s。耗散系数A、浅滩效应系数γ设置为10^–6 m²·s^–1。网格设置满足Courant-Friedrichs-Lewy(CFL)条件: Δx/Δt远大于内孤立波传播速度。采用蛙跃格式离散控制方程迭代求解。

弱非线性模型为两层模型, 上下层水体密度参考物模实验的层结条件, 分别取为1 030 kg/m³与1 050 kg/m³, 上层厚度h₁=0.05 m, 下层厚度h₂=0.25 m。采用输入波形的方法, 在左侧x=1.5 m处输入eKdV方程中内孤立波的第一模态解析解:

$ \eta = \frac{{{\eta _0}}}{{b + (1 - b){{\cosh }^2}\nu (x - Vt)}}. $

(7)

解析解中, η₀表示内波引起的最大界面起伏, 非线性相速度为$V = {c_0} + \frac{{{\eta _0}}}{3}(\alpha + \frac{1}{2}{\alpha _1}{\eta _0})$, 特征半波宽为${\nu ^2} = \frac{{{\eta _0}(\alpha + \frac{1}{2}{\alpha _1}{\eta _0})}}{{12\beta }}$, 特征系数为$b = - \frac{{{\eta _0}{\alpha _1}}}{{2\alpha + {\alpha _1}{\eta _0}}}$。

本文设置数值水槽长度大于内孤立波传播距离, 内孤立波在传播中未碰撞边界, 因此采用固壁边界条件。

本文工况设置见表 1。设置了效应组E1—E4, 效应组中不存在地形, 在保留控制方程的时变项与线性项的基础上, 仅考虑特定效应项, 探究平方非线性效应、立方非线性效应与频散效应的影响。设置了地形组H1—H4, 采用不同的入射波振幅η₀与地形高度H₀, 以实现内孤立波过陆架陆坡地形的不同演化过程模拟。参考Michallet等^[22]使用水深D将入射波振幅η₀无量纲化, 划分为η₀/D < 0.1的小振幅情况与η₀/D > 0.1的大振幅情况, 每种地形高度下采用两种入射波振幅η₀=0.01 m, η₀=0.05 m分别作为小振幅情况与大振幅情况。

针对H1—H4实验组, 设置了相同条件的全非线性的MITgcm模型实验组, 用于对比弱非线性RLW模型, 评估模型差异。

1.2 数据使用情况

在南海动力系数分布的讨论中, 温盐数据使用了Global Ocean Forecasting System 3.1(GOFS 3.1)的2015年再分析资料, 该系统基于HYCOM数据, 结合Navy Coupled Ocean Data Assimilation (NCODA)数据进行数据同化。水平分辨率为1/12°, 垂向共40层, 采用1月与7月的月平均数据以代表夏季与冬季的不同情况。地形数据采用分辨率为1/60°的ETOPO1数据。

1.3 分析方法

内孤立波传播过程中的能量分析方法参照Sutherland等^[34], 基于流速场与密度场可分别计算内孤立波的动能E_k(KE)与有效势能E_a(APE)。总能量E=E_k+E_a, 由于弱非线性模型无法计算流速场分析动能, 本文采用了内孤立波在传播过程中动能与有效势能相等的假设以计算总能量^[35]。能量计算公式为

$ {E_{\text{k}}}(x, t) = \int\limits_{ - D}^0 {\frac{1}{2}{\rho _0}({u^2} + {w^2}){\text{d}}z} , $

(8)

$ {E_{\text{a}}}(x, t) = \int\limits_{ - D}^0 {g\left( {\rho - \overline \rho } \right)\left( {z + {h_1}} \right){\text{d}}z} , $

(9)

其中, ρ₀为Boussinesq近似下密度, $ \overline \rho $为背景密度, u、w分别为水平速度、垂向速度。

弱非线性模型中动力系数由式(2)—(6)计算。南海的上下层厚度通过指数前差梯度(ELG)法识别跃层计算^[36]。

2 结果与分析 2.1 非线性与频散效应

E1至E4模拟结果如图 2所示, 从图 2(a)中可得, 非线性效应使内孤立波在传播过程中发生变形, 波前变陡, 波宽变窄, 波后变化相对较小。内孤立波不断陡化并发生裂变, 形成振幅由大到小排列的波列。这一演化过程对应内潮生成内孤立波的分裂机制, 内潮在深水传播过程中受非线性效应的影响, 发生陡化, 生成振幅大, 波长相对较短的内孤立波波列。图 2(b)对应高阶非线性效应, 不同于非线性效应, 内孤立波波前变化相对较小, 波后上凸前涌, 呈翻转趋势, 随后发生裂变生成振幅由大到小排列的线性波列。这一过程对应内孤立波在浅水中波后上凸发生裂变的演化过程。图 2(c)为频散效应, 内孤立波在传播过程中发生频散, 生成振幅由大到小的波列, 波与波的间距不断变大。图 2(d)中, 内孤立波保持波形匀速向前传播, 说明内孤立波是非线性波陡化与线性波频散相平衡的结果。

2.2 不同模型结果

图 3为地形高度为0.16 m的H1组模拟结果, 其中虚线为跃层位置, 灰色区域为底地形。小振幅情况下, 两模型中内孤立波在深水区的传播速度均为0.095 m/s。在内孤立波传播至陆架陆坡区后, 由于局地水深变化, 内孤立波传播速度变慢, 全非线性模型中内孤立波波速由0.095 m/s减小至0.076 m/s, 弱非线性模型波速减小至0.08 m/s, 相对全非线性模型在陆架地形上传播速度快。内孤立波在陆架陆坡传播过程中受地形影响小, 波形无显著变化。大振幅情况下, 全非线性模型结果中内孤立波传播速度也相对弱非线性模型快。在陆坡前的稳定传播过程中, 两模型差异小。在15~25 s, 内孤立波传播至陆架坡折前, 弱非线性模型振幅减小, 波后发生上凸, 波前变平, 倾斜程度与陆坡坡度相当。演化过程与高阶非线性效应过程一致。25 s后内孤立波传播至陆架地形, 形成波列, 其中头波为波前平波后陡的非对称结构, 尾波由大到小依次排列。全非线性模型在15 s后的演化过程相对弱非线性模型存在明显差异。在15~25 s, 内孤立波发生陡化, 振幅显著增大, 波后形成上凸, 但上凸高度小于弱非线性模型, 波后持续翻转倾斜, 在30 s发生翻转破碎, 头波振幅随破碎减小, 并受地形影响发生变形, 随后形成波前平波后陡的非对称结构, 波长大于弱非线性模型结果。尾波在破碎后逐渐生成, 仅存在2个尾波, 明显少于弱非线性模型。

图 4—图 6分别展示了地形高度为0.18 m、0.20 m、0.22 m的H2—H4组模拟结果。小振幅情况下, 随地形高度增加, 弱非线性与全非线性模拟结果中内孤立波演化过程的转变规律一致: 内孤立波在陆架区的传播速度变慢, 由不受地形影响转变为发生变形, 后转变为发生裂变。立方非线性效应引起的内孤立波陡化与裂变现象增强。但弱非线性模型结果中陆架地形上的头波波长相对较大, 最后一个尾波振幅偏大。大振幅情况下, 随地形高度的增加, 全非线性模型中内孤立波在陆架坡折前的陡化作用加剧, 波后翻转强, 逐渐垂直于地形, 破碎混合过程更剧烈, 波后不稳定界面的水平尺度增加, 尾波列不明显, 头波随破碎混合加强, 振幅减小, 波前变平, 逐渐演变为上凸型内孤立波。弱非线性模型无法有效模拟内孤立波的破碎混合, 在内孤立波的陡化发生翻转破碎阶段, 弱非线性直接发生陡化与裂变, 不断在波后发展出新的尾波以频散能量。陆架坡折前尾波能量明显增强, 振幅持续增大。

地形高度为0.22 m的H4组中内孤立波与地形相互作用结果最明显, 弱非线性与全非线性模型模拟结果差异大, 因此分析H4组两模型模拟结果中的能量分布规律以探讨模型差异。图 7为小振幅情况下的能量时空分布。在0~30 s的陆坡坡折前的传播过程中, 内孤立波波形不变, 两模型能量分布无显著差异, 在陆架坡折处(x=4 m), 全非线性模型能量密度增强, 对应内孤立波发生陡化, 振幅增大, 波后发生翻转的过程, 弱非线性模型无能量密度增强现象, 直接发生裂变生成尾波, 表现为30~60 s中3.5 m至5 m处的多个弱能量的分布。在陆架坡折后, 受地形影响, 头波速度变缓, 能量分布中的头波斜率变大。图 8为大振幅情况下的能量时空分布。两模型模拟结果中大振幅情况存在显著差异, 弱非线性模型在20~60 s内不断发生裂变生成尾波列, 能量集中于陆坡区域, 透射波能量显著减弱。全非线性模型中内孤立波翻转破碎发生位置为陆架坡折前, 能量密度增强现象更明显, 陆架地形上的透射波能量也相对较弱。30~60 s内, 内孤立波发生破碎, 湍流混合加强, 导致能量耗散的加强。破碎后能量分布中存在两个反射波结构, 第一个反射波能量相对较强。透射波能量相对较弱, 头波后的界面不稳定, 但能量持续前传。两模型入射波与透射波波列中的头波能量结构一致。在陆架坡折附近的演化过程中, 弱非线性过程无法有效模拟破碎混合与反射过程, 以辐射能量的方式分裂尾波, 裂变产生的尾波列不符合实际。

图 9为H4组空间积分的总能量分布曲线, 使用初始入射波总能量E_i作能量归一化。小振幅情况下, 在0~25 s的传播过程中, 内孤立波能量持续衰减, 衰减速率恒定。在25~35 s的陆架坡折附近演化阶段, 弱非线性模型裂变产生波列, 能量迅速衰减, 5 s内衰减至初始能量的60%, 而全非线性模型能量无明显衰减速率的变化, 随内孤立波传播持续衰减。大振幅情况下, 弱非线性模型在15~27 s内能量衰减更剧烈, 最低仅为入射波能量的20%, 而全非线性模型发生剧烈破碎, 25 s后能量衰减速率加大, 35 s时总能量为入射波能量的62%。在能量分布规律中, 两模型差异大。当内孤立波发生破碎混合时, 全非线性模型可直接刻画破碎混合过程, 破碎加强能量耗散。弱非线性仅能通过波形演化, 发生裂变生成波列以频散能量, 能量耗散机制与全非线性不同。原因是弱非线性无法通过强非线性的破碎过程耗散能量。

2.3 动力系数分布

不同地形高度下的动力系数分布见图 10。图 10(a)为线性相速度的空间分布。受地形影响, 在陆架陆坡区, 线性相速度随地形高度增加而减小, 表现为内孤立波波速减小。图 10(b)为频散系数的空间分布, 在深水区频散系数大, 内孤立波受频散效应影响大, 但由于非线性作用的平衡, 内孤立波不发生频散。在陆架陆坡地形上频散系数量值迅速减小趋近于零, 不同地形高度下频散系数差异小, 说明在浅水区域内孤立波的频散效应相对较弱, 地形高度对频散效应强弱影响小。图 10(c)为平方非线性系数空间分布。在深水区平方非线性系数大, 非线性效应强, 易于内潮波陡化产生内孤立波。随地形高度增加, 平方非线性系数由负变正, 内孤立波极性发生转化, 由下凹型向上凸型转变。与深水区相比, 浅水区的非线性系数在量值上显著减小, 说明在浅水区域内孤立波的平方非线性效应相对较弱, 内孤立波不易发生波前变陡的演化过程。图 10(d)为立方非线性系数空间分布。在深水区立方非线性系数小, 随地形高度增加, 立方非线性系数量值增大, 绝对值增长了一个量级。说明立方非线性效应主导内孤立波的浅水演化过程。受立方非线性系数影响, 内孤立波波后变陡产生裂变, 这一动力机制与前文波形演变中的立方非线性效应组结果一致。

2.4 基于弱非线性模型的南海动力系数估计

弱非线性模型可计算南海的动力系数分布, 探讨内孤立波不同季节的演化规律与差异。图 11为南海1月与7月的跃层深度分布, 1月吕宋海峡东部与南海北部陆架陆坡区的跃层深度在100~200 m, 7月南海跃层在0~50 m, 即冬季跃层较深, 夏季跃层较浅。

图 12为1月南海的动力系数分布。图 12(a)为线性相速度分布, 受深跃层影响, 南海深水区的线性相速度平均在2~3 m/s, 这一范围与Grimshaw等^[37]分析的深水区线性相速度为2.5 m/s, Cai等^[38]分析所得线性相速度范围2.25~2.75 m/s结果相符。在水深变化大的陆坡区, 线性相速度显著下降至1 m/s左右, 随着水深进一步减小, 陆架区的波速减小到0.5 m/s左右。南海内孤立波典型观测结果中, 内孤立波的传播速度平均为3.23 m/s^[1], 大于本文结果, 其差异来源于背景流与内孤立波振幅引起的非线性相速度。图 12(b)为频散系数分布, 吕宋海峡东部频散系数大于2×10⁵ m⁶·s^–1, 内孤立波受频散效应影响大, 说明在吕宋海峡东部, 内孤立波难以维持波形, 与Liao等^[32]结论一致。南海中部与南部的深水区频散系数在(0.5~1.5)×10⁵ m⁶·s^–1范围内, 内孤立波受频散效应的影响, 南海北部的陆架陆坡区频散系数趋于零, 内孤立波在陆架陆坡区不易发生频散。平方非线性系数见图 12(c), 在吕宋海峡东部, 平方非线性系数在0~2×10^–2 s^–1范围内, 平方非线性效应弱。在南海中部与南部, 平方非线性系数大, 易于内孤立波的生成, 同时受强频散效应的影响, 非线性效应与频散效应平衡, 内孤立波可以保持波形传播。在陆架陆坡区, 平方非线性系数由负变正。图 12(d)为立方非线性系数分布, 与平方非线性系数相反, 在深水区立方非线性系数趋近于零。

图 13为7月南海的动力系数分布, 其中图 13(a)表示线性相速度分布。与1月分布规律一致, 南海中部与南部线性相速度平均在1~2 m/s的范围内, 但超过2 m/s的范围相对较小, 整体线性相速度小于1月。7月南海北部的陆架陆坡区波速相对1月较大, 为0.5 m/s左右, 内孤立波传播速度相对较快。7月频散系数分布与1月类似, 见图 13(b)。平方非线性系数分布见图 13(c), 明显不同于冬季, 7月南海的非线性系数整体在4×10^–2 s^–1以上, 在深水区平方非线性效应强。在南海北部陆架陆坡区, 平方非线性系数大于零的区域向近岸靠近, 说明夏季内孤立波易在近岸处发生极性转化, 陆架陆坡区的极性转化发生频率小。立方非线性系数分布见图 13(d), 与1月分布规律相似, 在深水区立方非线性系数小, 在浅水区系数大。与1月不同, 7月南海中部深水区的立方非线性系数分布范围为(–3 ~ –5)×10^–4 m^–1·s^–1, 这一范围显著大于1月–2×10^–4 m^–1·s^–1范围, 说明夏季内孤立波受立方非线性效应影响大。

南海21°N断面为典型的内孤立波生成演化传播路径^[2], 图 14为该断面陆架陆坡区的动力系数分布。图 14(a)为线性相速度分布, 在深水区1月线性相速度高于7月, 说明在冬季深水区内孤立波传播速度快, 陆架区传播速度1月与7月均在0~1 m/s左右。在内孤立波向陆架陆坡区的传播过程中, 随水深的减小, 内孤立波线性相速度减小, 传播速度变慢, 这一规律已在观测中证实^[3]。图 14(b)为频散分布, 深水区1月频散系数显著大于7月, 随地形增高趋于零, 说明冬季深水区内潮波频散强, 内孤立波能量频散强, 不易生成大振幅内孤立波。图 14(c)与图 14(d)分别为平方非线性系数与立方非线性系数, 两系数绝对值7月均显著大于1月, 说明7月内孤立波受非线性效应影响大。此外, 7月频散系数相对1月小, 说明在非线性效应与频散效应的共同调制下, 内潮在夏季更易于发生陡化产生内孤立波且振幅更大。夏季内孤立波的发生次数大于冬季, 这与Chen等^[39]观测结果一致。1月平方非线性系数由负变正, 表现为内孤立波易发生由下凹型向上凸型的极性转化。随地形高度增加, 平方非线性系数趋近零, 而立方非线性系数1月与7月分别保持在2×10^–4 m^–1·s^–1, 4×10^–4 m^–1·s^–1, 在浅水区域量值增大了1个量级, 说明在陆架区平方非线性效应不再显著影响内孤立波的演化过程, 立方非线性效应主导陆架陆坡区的演化过程, 内孤立波易发生波后变陡, 变形上凸, 易裂变为波列在陆架上传播。

3 结论

本文应用两层的弱非线性RLW模型与全非线性MITgcm模型, 模拟了不同振幅、地形高度条件下内孤立波在陆架地形上演化的过程, 对比了两模型的差异, 探讨了内孤立波非线性与频散效应对波形演变的内在动力机制, 并基于弱非线性模型探讨了南海的非线性内波动力特征。研究发现:

1) 内孤立波通过频散效应与非线性效应相平衡保持波形不变。平方非线性效应使内孤立波波前发生陡化, 演化为内孤立波波列, 立方非线性效应使内孤立波波后发生陡化上凸, 演化为内孤立波波列, 频散效应使内孤立波生成振幅由大到小的频散波列。

2) 内孤立波在过陆架陆坡地形的演化中波前变平, 波后变陡, 波尾前涌, 最终以波列形式传播。相比于全非线性模型, 弱非线性模型模拟头波结果一致。不发生破碎的情况下, 弱非线性模型模拟结果较好。但发生破碎过程时, 弱非线性无法通过强非线性的破碎混合过程耗散能量, 只以频散能量的方式分裂尾波耗散能量, 演化为不符实际的波列。

3) 内孤立波在陆架陆坡地形的浅化过程中, 随地形高度增加, 陆架区频散系数由(0.5~1.5)×10^–4 m⁶·s^–1减小到0, 平方非线性系数由负转正, 量值趋近于0, 而立方非线性系数在4×10^–4 m^–1·s^–1的基础上, 绝对值增长了一个量级, 说明传播过程中频散效应与平方非线性效应较弱, 立方非线性效应主导陆架地形上内孤立波的演化过程。

4) 通过对比南海夏季与冬季动力系数空间分布, 发现内孤立波在传播过程由于夏季平方非线性效应、立方非线性效应与频散效应较强, 其在夏季更易发生陡化与裂变, 波列发生频率高。