线性稳定性分析(linear stability analysis, LSA)方法被广泛应用于大气与海洋研究中, 如航空动力学、机械工程力学及大气科学等。本文主要探讨线性稳定性分析方法在海洋, 尤其是小尺度不稳定过程研究中的应用。在层化海洋中, 上下层水体之间由于流速不同会产生强烈的剪切^[1]。同时, 当存在近惯性内波^[2-3]、内潮^[4-7]、高频内孤立波^[8]等波动时, 水体会发生不稳定过程, 从而产生湍流及垂直混合。影响海洋表层至次表层温度、盐度等水文变量及叶绿素等生化变量的分布。不稳定过程导致的跨等密度面湍流混合可能还会影响厄尔尼诺和南方涛动(El Niño-Southern Oscillation, ENSO)以及全球气候的发展^[9-10]。因此, 有关不稳定的研究是物理海洋研究中的热点问题之一。不稳定过程在海洋中普遍存在, 与其相关的不稳定特征可以由观测得到。但这对观测能力的要求极高, 尤其是当需要较大规模和较长时间的数据时, 直接观测总显得力有不逮。因此, 人们提出利用当前常规水文观测来研究不稳定过程, 线性稳定性分析方法发挥了显著的作用。

以Taylor-Goldstein方程(T-G方程)为理论基础, 将其转换为一个矩阵特征值问题后进行计算, 从而实现对观测流体的线性稳定性分析, 此举成功地将理论应用到了实际当中。通过LSA方法, Moum等^[11]假定不稳定模态与所测内孤立波的水平尺度相比较短, 对俄勒冈大陆架的内孤立波结构进行了详细观测, 发现这些波的界面会随着不稳定过程发生破裂, 并解释了湍流的结构与生成过程。在Moum等^[11]方法的基础上, Li等^[12]和Liu等^[13]考虑到小尺度湍流的影响, 对T-G方程进行了改进, 分别增加了水平与垂直涡流黏度和扩散系数项。刘志宇^[14]应用此方法研究了黄海与克莱德海两个典型高潮能驱动陆架海(强潮驱陆架海)中的湍流与混合特征并探讨了其发生机理。Sun^[15]通过LSA方法, 研究了1991年0°N, 140°W处21 d连续观测中, 层化剪切不稳定对赤道混合层下湍流循环(深循环)的影响。Smyth等^[16]用此法分析了2008年冬0°N, 140°W处剪切不稳定的昼夜变化及其导致的混合的昼夜变化(通常2~3 h)现象, 发现了混合昼夜变化的不对称性。Liu等^[17]也通过此法对赤道太平洋冷舌区的温跃层内的不稳定进行了研究, 发现不稳定现象不仅发生在温跃层上部, 还可以深达温跃层中心甚至温跃层以下, 在热带不稳定波(Tropical instability waves, TIWs)和拉尼娜(La Niña)存在的条件下因具有更强的剪切不稳定而发生得更为频繁, 深度也较平时更深, 产生强混合。最近, Liu等^[18]发现障碍层中也存在强的剪切与混合, 同样在拉尼娜年发生得更为频繁。本研究中, 我们将以热带太平洋赤道142°E进行的垂直剖面观测资料为例, 介绍LSA方法的具体应用并分析其结果。

1 LSA方法的应用 1.1 特征值问题的转换

已知T-G方程如下(详细推导见附录)。

$ {\hat w_{zz}} + \left[ {\frac{{{N^2}}}{{{{\left( {U - c} \right)}^2}}} - \frac{{{U_{zz}}}}{{U - c}} - {k^2}} \right]\hat w = 0 , $

(1)

其中, N为浮力频率, U为平均流, k为波数, w为垂向流速, c为相速度。

现将其转换为矩阵特征值问题。在实际观测下, N²和U是可以提供的, k可以假设, 那么则提供了一个以c或ω为特征值, 以w为特征函数的特征值问题。这个方程与典型的特征值问题方程不同, 我们可以进一步处理如下。

实际海洋中的流体在垂直方向上往往不是平行流。假设存在某一平行流, 且符合Boussinesq近似, 即低速流动中忽略压强, 只考虑密度变化受温度变化影响。

$ U=U(z), $

(2)

$ N=N(z). $

(3)

假设不稳定是二维的, 即它们不在y方向上变化, 仅在x和z方向上变化。由基本方程组(附录A10-A12)结合不可压缩流体连续性方程$ \frac{{\partial u}}{{\partial x}} + \frac{{\partial w}}{{\partial z}} = 0 $推导得到平行层化剪切流的二维线性扰动方程组如下:

$ \frac{{\partial \xi }}{{\partial t}} = - U{\xi _x} + {U_{zz}}{\psi _x} + {b_x} , $

(4)

$ \frac{{\partial b}}{{\partial t}} = - U{b_x} - {N^2}{\psi _x} , $

(5)

$ {\nabla ^2}\psi = \xi . $

(6)

其中, u=$ - {\psi _z} $, w=$ {\psi _x} $, ξ(x, z, t)为涡度, ψ(x, z, t)为相应的流函数, b=–ρg/ρ₀(b为浮力; ρ₀为密度特征值, 是一个常数)。

因为只考虑了x和z方向的变化, 故可以引入水平波数k和相速度c将其扩展为常规模态, 即

$ \psi \left( {x, z, t} \right) = \hat \psi \left( z \right){e^{ik\left( {x - ct} \right)}}, $

(7)

$ b\left( {x, z, t} \right) = \hat b\left( z \right){e^{ik\left( {x - ct} \right)}}. $

(8)

其中, $\hat \psi $和$\hat b$是结构函数。通过如上变换可以将其转换为一个等价于T-G方程的矩阵特征值问题, 即

$ \left[ \begin{gathered} U - {U_{ZZ}}{\nabla ^{ - 2}}\;\;\;\; - 1 \hfill \\ {N^2}{\nabla ^{ - 2}}\;\;\;\;\;\;\;\;\;\;\;\;U \hfill \\ \end{gathered} \right]\left( {\begin{array}{*{20}{c}} {\hat \xi } \\ {\hat b} \end{array}} \right) = c\left( {\begin{array}{*{20}{c}} {\hat \xi } \\ {\hat b} \end{array}} \right) . $

(9)

如果给定了U和N剖面(也就是我们能够获得的某站点的流体垂向剖面), 那么此特征值问题可进行数值求解。其中, 扰动的相速度c为特征值, c=c_r+ic_i, 其可能是复数, $ \left( {\begin{array}{*{20}{c}} {\hat \xi } \\ {\hat b} \end{array}} \right) $是特征向量。

当c为复数, 且c_i > 0时, kc_i代表扰动模态的生长率(c_i < 0时, 模态衰减, 不构成不稳定模态), 此时ψ(x, z, t)可呈指数增长, 意味着扰动呈指数剧烈发展, 则会发生不稳定。不稳定发生的位置为临界层深度, 此时流速U(z)与相速度c_r相等。

当海洋背景状态同时满足多种不稳定条件时, 可能会产生多个不稳定模态。这些不稳定模态往往以其临界层深度归为显著不同的几类, 每一类成为一个模态族^[15-16]。一个观测得到的流剖面可能含有若干个模态族。在每个模态族中, 又可选出具有最大生长率的不稳定模态, 来代表该不稳定模态族的不稳定模态。确定潜在的不稳定模态需要采取一定的标准^[16], 包括“截止深度”标准、“截止不稳定生长率”标准、“截止波长”标准、理查森数(Richardson, Ri)标准、临界层深度标准等。

截止深度标准。该标准应排除距观测剖面的上下边界一定距离内的模态, 因为由于边界条件的设置, 及数据在边界附近的外推等因素会导致潜在虚假模态。有敏感性研究表明, 在边界附近获得到的不稳定模态的位置可能会偏离其真实位置^[19]。

截止不稳定生长率标准。应该保证不稳定的成长率大于观测数据的时间分辨率。

截止波长标准。该标准取决于数据的垂向采样间隔。无黏, 非扩散, 层化剪切层的不稳定性通常具有约为剪切层厚度2π倍的波长^{[16, 20-21]}。因此, 可以根据波长来剔除一些虚假模态。任何模态族都必须包括至少一个波长大于最快生长模态的解析模态和至少一个波长小于最快生长模态的解析模态, 否则该模态族可被排除掉。

理查德森(Ri)数标准。Ri定义为浮力频率平方与剪切平方的比值, 即Ri=N²/S², 其中$ {S^2} = {\left( {\frac{{\partial u}}{{\partial z}}} \right)^2} + {\left( {\frac{{\partial v}}{{\partial z}}} \right)^2} $。作为发生不稳定的一个判据, 剖面的最小Ri应该小于Miles-Howard临界值0.25。此外, 在潜在不稳定模态的临界层深度附近(±1/14波长), 最小Ri也要求小于0.25, 这是基于典型的Kelvin-Helmholtz不稳定的性质确定的^[22]。

临界层深度标准。内插或外推得到的速度剖面仍旧是离散的, 所以它可能无法保证满足临界层深度的定义: U(z)−c_r =0。一般可以取∣U(z)−c_r∣≤某一给定值, 如0.01。在此标准下, 可能存在不止一个满足上述标准的临界层。这些临界层可以全部保留下来, 进一步分析。

1.2 其他有关LSA的影响因素 1.2.1 上下边界条件

为了尽可能合理地刻画边界层对研究的影响, 前人分析研究了多种上下边界条件。其中, “刚盖”近似是一种较为经典的海表边界条件, 即在海表处要求垂直速度W(0)=0。它包含了压力变化, 同时去除了海表惯性重力波, 排除了海表高度变化的动力学影响。在Liu等^[17]有关赤道太平洋冷舌区温跃层混合的研究中, 上边界和下边界均设定垂直速度为0。

除了“刚盖”近似外, 还有辐射边界条件应用也较为广泛。比如, Hazel^[21]曾分别在上边界和下边界处施加了如下辐射边界条件:

$ \frac{{\partial p}}{{\partial y}} = 0 . $

(10)

它们代表了在没有固壁的情况下的远场条件, 且在东西方向上呈周期变化。Sun等^[15]应用“刚盖”近似条件将海表扰动的垂向速度设为0, 只在下边界约200 m深度处施加过辐射边界条件。Smyth等^[23]用上下辐射边界条件对非对称层化急流进行了稳定性分析和数值模拟, 研究了该体系中失稳和重力波辐射的一些基本性质。

值得一提的是, 在不同边界条件中, 不稳定的临界值标准依旧适用。例如, Miles^[24]和Howard^[25]证明, 在固壁边界条件下, 如果Ri处处大于临界值0.25, 则流体对于小振幅扰动是稳定的。在辐射边界条件下, Miles-Howard理论依旧适用, 即某一深度上, Ri小于0.25是剪切不稳定发生的必要条件。

1.2.2 混合对LSA的影响

混合中的涡流黏度效应与扩散效应等也在LSA方法中也扮演着重要角色。比如, Smyth等^[16]在湍流昼夜变化研究中发现, 在夜晚、湍流最强的时刻, 涡流黏性项会阻碍不稳定的产生。在不同条件下, 涡流黏度的影响会很微妙, 它可能还会引起层化剪切流的不稳定。Li等^[12]考虑了涡流黏性与扩散项对于LSA方法的影响, 将涡流黏性系数A与扩散系数Z的比值定义为普朗特数Pr, 并分了3种情况进行讨论。结果表明, 只有在弱的湍流混合效应且Ri接近0.25时, 湍流黏性项与扩散项对于剪切不稳定的影响才会与经典的Thorpe^[26]理论基本相符。一旦Ri偏离0.25, 湍流的垂直结构会导致结果出现差异。将Ri降低到远低于其临界值0.25的水平, 可以改变甚至逆转涡流黏度和扩散项的影响。

2 小尺度不稳定分析 2.1 数据处理

本研究使用了赤道太平洋142°E处垂直剖面观测数据。所选剖面的观测时间为2017年12月16日17时。观测变量包括由下放式声学多普勒流速剖面仪(lowered acoustic Doppler current profiler, LADCP)观测的流速数据及由温盐深测量仪(conductivity- temperature-depth system, CTD)观测的温度(T)、盐度(S)和深度数据。其垂直分辨率分别为1 m和10 m。对上述数据进行样条插值, 统一到1 m间隔格点上。其中, CTD数据表层至5 m、LADCP数据表层至10 m为外推。

我们将流矢量依次投影到0°~180°不同的α方向分别构建特征值问题, 从而得出不同波数方向α上的不稳定特征。这里方向0°代表正东, 90°代表正北; 因–180°~0°与0°~180°对称, 故不需重复计算。

为减少计算量, 可以取一定度数为间隔。当确定波数方向α后, 可得到该方向上的平行流垂向剖面, 结合密度剖面, 则可计算出该波数方向上的不稳定模态的基本特征, 如不稳定成长率大小、波数大小及其临界层深度。

2.2 结果分析 2.2.1 剖面性质

图 1(a)为该剖面的流速分布, 黑线和蓝线分别代表纬向(u)、经向(v)速度分量。u分量东向最大值为0.74 m·s^–1, 位于–14 m, 西向最大值为0.23 m·s^–1, 位于–141 m; v北向分量最大值为0.16 m·s^–1, 位于–101 m, 南向最大值为0.70 m·s^–1, 位于表层(插值得)。图 1(b)为位温(黑线)和盐度(蓝线)剖面。由图可见, 上300 m盐度变化不大, 位温随深度增加而逐渐降低。相应地, 图 1(c)为位势密度, 随着深度增加而增加。图 1(d)为浮力频率的平方, 表征流体的层结变化。上层层结较弱, 95~100 m和140 m左右层结最强, 利于水体趋于稳定。

图 2(a)和2(b)为流体剪切的平方, 蓝色阴影分别代表剪切的纬向分量U_z²和经向分量V_z², 黑色实线为总剪切的平方。可看出, 纬向剪切分量占据主导地位。剪切的极大值出现在46 m、90 m和146 m附近, 且主要由纬向剪切分量贡献。

图 2(c)为Ri^–1, 蓝色虚线为阈值0.25。由图可见, 5个区域存在Ri低于0.25的现象: 分布在10 m、19~23 m、42~49 m、85 m左右和175 m左右。在表层, Ri趋近于0, 剪切非常强, 意味着存在很强的不稳定过程。然而, 不稳定过程受边界临近效应调制^[21], 故我们不考虑边界处, 认为50 m、85 m、175 m这3个深度附近是不稳定存在的主要位置。

2.2.2 线性分析结果

在本例的线性分析中, 我们将截止深度设置为300 m, 截止不稳定生长率取1 h^–1, 波长取5: 5: 200 m、210: 10: 400 m、420: 25: 1 000 m、1 030: 50: 2 000 m、2 050: 100: 5 000 m, 临界层深度标准取∣U(z) –c_r∣≤ 0.015。实测流方向α取0°~180°(0°代表正东), 以15°为间隔, 依次计算。

剖面的生长率和临界深度随经向、纬向波长分布如图 3所示。为了更直观地展示不稳定生长率与临界深度的对应关系, 我们将经向波长截取–6 000 m至6 000 m。纬向波长截取0至220 m, 与–220 m至0 m关于(0 m, 0 m)呈中心对称结构, 因此–220 m至0 m未画出。选取第3组来展示生长率和临界深度计算结果。由图 3(a)可见, 该剖面存在多个不稳定模态族, 且每个族存在一个最快生长模态, 分别对应不稳定生长率的各个极值点, 主要位于纬、经向波长(10, –2 050)m、(45, –770)m、(50, –370)m、(50, –1 000)m、(75, –150)m处等。

图 3(b)为相应的临界深度。值得注意的是, 不同极值点或许对应同一临界深度。由图可见, 极值点对应的临界深度大致在50 m左右、80 m左右以及150 ~ 200 m范围内, 此时不稳定的生长率为最大值。最快生长模态在较长的时间内会占据主导地位, 而在现实海洋中, 扰动会在有限时间内增长, 所以除最快生长模态外的模态也会有显著的振幅。因此, 每个模态族都有可能生成一个波包(wave packet)。从图 3不稳定生长率分布图也可看出被俯视投影出的三维波包。

为了更直观地得到不稳定生长率和临界深度分布特征, 我们分别将最快生长模态和所有模态的生长率与临界深度用直方图进行统计, 如图 4所示。不考虑上边界, 我们可以得到3个明显的不稳定发生的临界深度, 分别为47 m、84 m和174 m。同时, 这些深度与Ri判定出的不稳定发生深度(~50 m、~85 m和~175 m)对应良好。

此外, 除了上述剖面性质, 我们也使用相同的LSA方法分析了其他剖面的不稳定特征。研究发现, 在其他剖面上, 不稳定发生的临界深度与Ri判定标准所得深度也具有良好的对应关系。例如, 在2017年12月15日21时, 如图 5(a)和图 5(b)所示, 通过剖面的Ri分布可以看出, 不稳定集中发生在30 m附近和60~70 m深度范围内, 并且这些深度与LSA方法计算所得不稳定发生的临界深度相一致。

图 5(c)和图 5(d)分别为2017年12月17日23时Ri及临界深度的分布, 通过LSA方法计算得到的不稳定过程主要分布在20~30 m、40~50 m以及80 m附近, 与Ri小于临界值0.25所对应深度相一致。

因此, 通过LSA方法, 我们可以得到流体某一剖面(也可进一步延伸分析多个剖面)的不稳定性质, 如层结、剪切、生长率、临界深度等。这为我们下一步深入研究不稳定过程发生机制等奠定了基础。

3 总结与展望

至此, 我们系统地介绍与总结了LSA方法的产生背景、具体实现过程、其可能存在的影响因素及相关应用, 并通过剖面举例, 阐述了LSA方法在实际海洋中的具体应用过程, 首次给出了赤道142°E处观测数据中流体的不稳定特征。作为经典的、较为常见的稳定性分析手段, LSA方法对今后海洋动力过程研究具有极高的参考价值, 可以得到广泛适用。比如, 通过LSA方法可以找到海洋中小尺度不稳定过程, 而这些不稳定过程可能使水体失稳, 进一步引起湍流混合, 实现物质、能量的再分配。同时, 不稳定过程带来的混合效应甚至可能影响更大尺度的海洋-大气动力循环。

在今后的研究中, 我们将通过对LSA方法中的经验参数和相关变量进行更为精确的参数化调整。同时, 适当考虑边界条件及混合效应等对LSA方法的影响, 对其进一步修正和完善, 得到更加完整、准确的LSA方法并将其投入到研究应用中, 从而确定哪些剖面可以更广泛且有效地表征被研究水体的具体特征及相关性质。

附录 Taylor-Goldstein方程的推导

描述海洋波动的基本方程组为:

$ \frac{{\partial {\boldsymbol u_{\text{h}}}}}{{\partial t}} + f \times {\boldsymbol u_{\text{h}}} = - \frac{1}{{{\rho _0}}}{\nabla _{\text{h}}}p , $

(A1)

$ \frac{{\partial w}}{{\partial t}} = - \frac{1}{{{\rho _0}}}\frac{{\partial p}}{{\partial z}} - g\frac{\rho }{{{\rho _0}}} , $

(A2)

$ \frac{{\partial \rho }}{{\partial t}} - \frac{{{\rho _0}{N^2}}}{g}w = 0 , $

(A3)

$ {\nabla _{\text{h}}} \cdot {\boldsymbol u_{\text{h}}} + \frac{{\partial w}}{{\partial z}} = 0 . $

(A4)

其中, (u_h, w)为波动流速(u_h=(u₁, u₂)为水平分量, w为铅直分量), p为压力。N(z)为Brunt-Väisälä频率, 即浮力频率, ${N^2} = - \frac{g}{{\bar \rho }} \cdot \frac{{{\text{d}}\bar \rho }}{{{\text{d}}z}}$。g为重力加速度, $ \bar \rho $为平均流的密度, $ \frac{{{\text{d}}\bar \rho }}{{{\text{d}}z}} $表示密度在垂直方向z上的导数。将科氏参数取近似, 即f≈2Ω。$\varPhi = gz$, 为重力势。

由方程组(A1—A3)展开得到线性扰动方程组如下:

$ \frac{{\partial {\boldsymbol u_{\text{h}}}}}{{\partial t}} + \boldsymbol U \cdot {\nabla _{\text{h}}}{\boldsymbol u_{\text{h}}} + w \cdot \frac{{\partial \boldsymbol U}}{{\partial z}} + f \times {\boldsymbol u_{\text{h}}} = - \frac{1}{{{\rho _0}}}{\nabla _{\text{h}}}p , $

(A5)

$ \frac{{\partial w}}{{\partial t}} +\boldsymbol U \cdot {\nabla _{\text{h}}}w = - \frac{1}{{{\rho _0}}}\frac{{\partial p}}{{\partial z}} - g\frac{\rho }{{{\rho _0}}} , $

(A6)

$ \frac{{\partial \rho }}{{\partial t}} + \boldsymbol U \cdot {\nabla _{\text{h}}}\rho - \frac{{{\rho _0}{N^2}}}{g}w = 0 , $

(A7)

$ {\nabla _{\text{h}}} \cdot {\boldsymbol u_{\text{h}}} + \frac{{\partial w}}{{\partial z}} = 0 . $

(A8)

为了便于介绍平均流U平行于波数方向的部分, 有$ {U_{||}}(z) = \frac{{\boldsymbol U(z) \cdot {\boldsymbol k_h}}}{{{k_h}}} $。后续为方便, 用U代替U_||。

地球自转的影响可以忽略, 即科氏参数f=0。拟设一水平波数k_h, 将式(A5—A8)代入$ \left\{ {\boldsymbol u, w, p, \rho } \right\} = $ $ \left\{ {\hat {\boldsymbol u}(z), \hat w(z), \hat p(z), \hat \rho (z)} \right\}\exp i({\boldsymbol k_h} \cdot \boldsymbol x - \omega t) $解得

$ - {\text{i}}(\omega - \boldsymbol U \cdot {\boldsymbol k_{\text{h}}})\hat {\boldsymbol u} + \hat w{\boldsymbol U_{\text{z}}} = - \frac{{\text{i}}}{{{\rho _0}}}{\boldsymbol k_{\text{h}}}\hat p , $

(A9)

$ - {\text{i}}(\omega - \boldsymbol U \cdot {\boldsymbol k_{\text{h}}})\hat w = - \frac{{{{\hat p}_{\text{z}}}}}{{{\rho _0}}} - g\frac{{\hat \rho }}{{{\rho _0}}} , $

(A10)

$ - {\text{i}}(\omega - \boldsymbol U \cdot {\boldsymbol k_{\text{h}}})\hat \rho - \frac{{{\rho _0}{N^2}}}{g}\hat w = 0 , $

(A11)

$ {\text{i}}{\boldsymbol k_{\text{h}}} \cdot \hat{\boldsymbol u} + {\hat w_z} = 0 . $

(A12)

此处$ {\hat w_z} = \frac{{\partial \hat w}}{{\partial z}} $, $ {\boldsymbol U_Z} = \frac{{\partial \boldsymbol U}}{{\partial z}} $等是导数简化符号。ω=iω_i+ω_r, 当ω_i > 0时, 振幅随时间呈指数增长, 不稳定发生。

由(A9, A12)乘k并在垂向上求导得

$ \left( {\omega - U\boldsymbol k} \right){\hat w_{zz}} + \hat w\boldsymbol k{U_{zz}} = - \frac{{\text{i}}}{{{\rho _0}}}{\boldsymbol k^2}{\hat p_z} . $

(A13)

由(A10, A11)联立消去$ \hat \rho $得

$ - {\left( {\omega - U\boldsymbol k} \right)^2}\hat w = {\text{i}}\left( {\omega - U\boldsymbol k} \right)\frac{{{{\hat p}_z}}}{{{\rho _0}}} - {N^2}\hat w . $

(A14)

将(A14)代入(A13)中消去$ {\hat p_z} $得

$ {\hat w_{zz}} + \left[ {{\boldsymbol k^2}\frac{{{N^2} - {{\left( {\omega - U\boldsymbol k} \right)}^2}}}{{{{\left( {\omega - U\boldsymbol k} \right)}^2}}} + \frac{{\boldsymbol k{U_{zz}}}}{{\omega - U\boldsymbol k}}} \right]\hat w = 0 . $

(A15)

相速度c=ω/k, 上述等式可写为

$ {\hat w_{zz}} + \left[ {\frac{{{N^2}}}{{{{\left( {U - c} \right)}^2}}} - \frac{{{U_{zz}}}}{{U - c}} - {\boldsymbol k^2}} \right]\hat w = 0 . $

(A16)

即Taylor-Goldstein方程。