南黄海非线性内波的数值模拟研究(2)——内潮的非线性演变*

司宗尚,范植松,于万春,杜 凌

(中国海洋大学海洋环境学院,山东青岛266100)

在本工作的先行部分[1],在南黄海沿着假想的直壁渠道使用连续层化海洋的内Kelvin波模型,并依据若干历史观测结果,对半日和全日内潮波进行了初步的数值模拟。在本文中考虑内潮波从AB侧边界的反射和非线性演变,以此作为内潮波从朝鲜半岛的若干岛屿和边界的反射和非线性演变过程的近似。文中作者仅考虑第一模态半日内潮波分别从渠道侧边界的A站和B站向西、西南和西北方向的反射和非线性演变过程(参见文献[1]中的图1),所使用的数学物理模型为适用于连续层化海洋并且考虑背景正压落潮流与涨潮流不同作用的一般化的KdV模型,简称之为GKdV模型[2-3]①Fan Z S,Shi X G,Liu A K,et al.Effectsof ebb and flood background currentson nonlinear internal solitury waves in South China Sea[J].Chin JOceanol Lnmol,2011,in press.。该模型是作者在国外学者研究成果的基础上作出的进一步发展[4-5],是适用于宽阔陆架海域和边缘海域的非线性内波模型。类似于在南海北部的研究结果[2-3]①,下文中的数值模拟结果再次证明了在内潮波和内孤立波的传播过程中背景流场的作用是十分重要的,并且区分背景正压落潮流与涨潮流的不同作用是必须的。下文中的数值模拟结果表明,南黄海的内孤立波集中分布于南黄海的南部,而在其北部极少出现内孤立波的主要原因是在南黄海的南部存在较强的背景斜压环流和很强的背景正压潮流,而且在涨潮时段(相对于朝鲜半岛的西海岸而言)背景斜压环流和背景正压潮流的方向基本相同,而在南黄海的北部这两种背景流都很弱。背景流的这种分布状态导致在内潮波非线性演变的GKdV模型中的非线性系数α在南黄海南部的涨潮时段出现较大值,数倍于其在南黄海北部的值。因而在南黄海南部内潮波裂变生成内孤立波的事件频频发生。

1 GKdV模型

在非线性内波的GKdV模型中,实际上内波的等密度面的铅直位移为η(x,t)Φ(z),这里z是向上为正的铅直坐标,Φ(z)是某一模态内波铅直位移振幅的铅直结构函数,而η(x,t)满足如下形式的连续层化的GKdV方程[1-3]①:

上式中,t为时间,x为水平坐标,c为长内波的相速,α为非线性系数,β为频散系数,Q(x)反映深度的缓慢变化和背景密度及背景流场剪切的水平变化,κ为平方底摩擦中的经验系数,v是湍流水平涡旋黏性的经验系数,而代表内波模态的垂直尺度。

铅直结构函数Φ(z)和长内波相速度c由如下特征值问题的解确定:

上式中N(z)和U(z)分别为Brunt-V?is?l?频率和背景剪切流,H为局地水深,而且对Φ(z)的值使用其最大值进行正则化。

系数α,β和Q定义如下:

(5)式中的下标‘0’的值表示在任一固定点x0的值,为方便计算取为原点,具体在本文中即为A站位或B站位。

在取下列变换后

方程(1)简化为

方程(8)是本文中模拟内潮传播的基本方程。在方程(8)的数值解获得之后,依据方程(7)在任一站位(x)等密度面的铅直位移η(x,s,z)可以写为:

而内波场的铅直速度w(x,s,z)可以近似写为:

依据下列流函数Ψ的近似关系式[6]:

可以得到内波场的水平速度u(x,t,z)近似为:

由于底摩擦和水平扩散,内波的传播总伴随着能量的丧失。对于底部摩擦应力,底部湍流边界层被参数化,使用Chezy形式或平方底摩擦经验表达式如下[4]:

上式中ρ为海水密度,u为边界层之外近底层的速度,而κ为经验系数(κ～0.001-0.002 6)。在GKdV模型中湍流水平涡旋黏性的经验系数v是1个较大的值。Liu等[7]在Sulu海孤立子模拟中取v=(10～30)m2·s-1;在New York湾的内波研究中,Liu[8]估计v为1 m2·s-1的量级;Sandstrom和Oakey[9]在Scotian陆架获得v的平均值为0.2 m2·s-1。

2 模型的初始化

考虑第一模态半日内潮波分别从渠道侧边界的A站和B站向西、西南和西北方向的反射和非线性演变过程(参见文献[1]中的图1),所考虑的站位名称列于表1和表2。自A站向西的4个站位分别为AW,AW 1,AW 2,AW 3;向西北的4个站位分别为ANW,ANW 1,ANW 2,ANW 3;向西南的4个站位分别为ASW,ASW 1,ASW 2,ASW 3。向西的每2个站位之间的距离约为27.5 km,而向西北和西南的每2个站位之间的距离约为38.9 km。实际上,AW,ANW和ASW站位的地理坐标均与A站位相同,即为(34°4′N,125°6′E)。为了区别内波不同的传播方向并且考虑到背景流的方向及量值不同,选用不同的名称以示区别。对于B站位(36°34′N,123°51′E),各站位名称的定义及每2个站位之间的距离与上述A站位的情况是类似的。在本文中使用GKdV模型进行模拟时,水平坐标的原点设置在A站位(或B站位),x轴的正向为向西(或向西北和西南)。

本文仅给出在8月份的数值模拟结果。本文中使用的温度、盐度资料及处理方法已在本工作的先行部分作了介绍。依据这些资料计算得到的A站位的密度的铅直分布和Brunt-V?is?l?频率的铅直分布分别示于文献[1]中的图2和图3。

图1 8月份在黄海海域斜压环流的表层流速分布图Fig.1 The distribution of surface current of baroclinic circulation in August in the Yellow Sea

本文使用的斜压环流数据为L ICOM 1.0模式[10-11]稳定积分950a后输出的线性插值到0.25(°)×0.25(°)网格的月平均值。图1为8月在黄海海域该斜压环流的表层流速分布图。图2和图3分别为8月该斜压环流的U分量在A站位和B站位的铅直剖面图。由于在A站位和B站位取x轴的正向为向西,所以该斜压环流的U分量为负值(见图1)。图1～3清楚地显示,斜压环流的强度在A站位和B站位之间相差悬殊,在A站位斜压环流的强度比其在B站位之值约大1个量级。

实际上,依据上述途径获得的资料显示,在6～9月份黄海的温度和盐度的气候分布状态,以及月平均斜压环流的流速分布状态均呈现类似的特性。本文以8月份的资料以及模拟的结果作为南黄海夏季的典型代表。

表1 自A站向西、西南和西北方向模拟内潮非线性演变经历的各站位名称和相关参数值Table 1 The station names and the parameter values used in the modeling of nonlinear evolution of the internal tidal wave from Station A in a westward,southwestward and northwestward direction respectively

表2 自B站向西、西南和西北方向模拟内潮非线性演变经历的各站位名称和相关参数值Table 2 The station names and the parameter values used in the modeling of nonlinear evolution of the internal tidal wave from Station B in a westward,southwestward and northwestward direction respectively

图2 8月在A站位斜压环流的U分量的铅直剖面图Fig.2 The vertical profile of the Ucomponent of baroclinic circulation in August at the station A

图3 8月在B站位斜压环流的U分量的铅直剖面图Fig.3 The vertical profile of the Ucomponent of baroclinic circulation in August at the station B

图4 8月在A站位表层潮流U分量的模拟结果Fig.4 The result of the U component of tidal current in the surface layer in August at the station A

本文使用的正压潮流数据取自文献[12]。杜凌[12]在其工作中采用正压POM 98模型模拟渤海,黄海和东海的潮波。温度和盐度设为常数(T=10℃,S=35),表面热通量,盐通量和风应力均设为0。底摩擦系数在渤海取为0.001,在黄海和东海取为0.002 5。考虑的区域范围为:24°N～41°N,117°E～131°E。在开边界上取水位边界条件。台湾海峡、对马海峡和琉球群岛岛链处的开边界数据是由附近的测站资料插值得到的。输出的潮流结果按照σ坐标分层为表层,中层和底层(σ值分别为0.0,-0.3,-1.0)。空间分辨率为10(′)×10(′)。考虑M2,S22个半日分潮和K1,O12个全日分潮的合成输出潮流结果。本文使用的正压潮流数据为8月1日1时～9月1日1时的模拟结果,数据的时间间隔为1 h。图4和图5分别为8月份在A站表层潮流U分量和V分量的模拟结果。图6和图7分别为8月份在B站表层潮流U分量和V分量的模拟结果。由于在中层和底层的输出潮流结果与表层的结果相差很小,故在以后的数值模拟工作中作者取表层的输出潮流结果作为正压潮流的代表值。图4～7表明,在黄海南部半日潮流强盛,并且在A站位潮流的强度比其在B站位之值大出2倍以上。

在铅直结构函数Φ(z)和长内波相速度c的确定过程中,即在特征值问题(2)的求解中,本文仅考虑月平均斜压环流的作用。求解特征值问题(2)的数值方法请见文献[13]。第一模态半日内潮波在AW,AW 1,AW 2,AW 3站位的铅直结构函数Φ(z)的铅直剖面图示于图8。而第一模态半日内潮波在BW,BW1,BW2,BW3站位的铅直结构函数Φ(z)的铅直剖面图示于图9。在图9中,在BW1,BW2,BW3站位的铅直结构函数明显不同于在BW站位的结果是由于水深变浅的缘故。

图5 8月在A站位表层潮流V分量的模拟结果Fig.5 The result of the V component of tidal current in the surface layer in August at the station A

图6 8月在B站位表层潮流U分量的模拟结果Fig.6 The result of the Ucomponent of tidal current in the surface layer in August at the station B

图7 8月在B站位表层潮流V分量的模拟结果Fig.7 The result of the Vcomponent of tidal current in the surface layer in August at the station B

图8 8月在AW(实线),AW 1(虚线),AW 2(点划线),AW 3(双划线)站位第一模态半日内潮波的铅直结构函数Φ(z)的铅直剖面图Fig.8 The vertical profiles ofΦ(z)of semi diurnal internalidal wave for the first mode at the stations AW(solid curve),AW 1(dotted curve),AW 2(dotted-and-dash curve)and AW 3(double dash curve)respectively in August

图9 同图8,但为BW(实线),BW 1(虚线),BW 2(点划线),BW 3(双划线)站位Fig.9 A s Fig.8,but fo r the stations BW(solid curve),BW 1(dotted curve),BW 2(dotted-and-dash curve)and BW 3(double dash curve)respectively

在系数α,β和Q的计算中,为了确定背景流本文不仅考虑月平均斜压环流的作用同时考虑正压潮流的作用。依据上述正压潮流的模拟结果,本文在各个站位分别选取在涨潮时段和落潮时段的潮流的最大值(绝对值)近似作为该时段正压潮流的特征值。然后将该特征值与月平均斜压环流值的合成结果作为在各个站位分别在涨潮时段和落潮时段的背景流。比较特殊的情况出现在ANW,ANW 1,ANW 2,ANW 3站位(见表1),即出现正压潮流的值(绝对值)超过长内波相速度c的情况。例如在ANW站位,涨潮时段潮流的最大值(绝对值)为-0.75 m/s(在表1中用圆括号列出),而长内波相速度c为0.469 m/s。此时背景流的值(在表1中没有列出)为-0.88 m/s,其绝对值超过长内波相速度c的值,因而内潮波处于不稳定状态①Fan Z S,Shi X G,Liu A K,et al.Effects of ebb and flood background currents on nonlinear internal solitury waves in South China Sea[J].Chin JOceanol Lnmol,2011,in press.[5]。在本文中不考虑这种不稳定状态。因此,在以上4个站位选取较小的正压潮流的特征值与月平均斜压环流值的合成结果作为其背景流值。非线性系数α,频散系数β以及长内波相速度c在各站位的计算结果列于表1和表2中。值得注意的是非线性系数α的特性。在各站位非线性系数α在涨潮时段的值(绝对值)均大于其在落潮时段的值。在表1中,这2个时段的非线性系数α的值相差约3倍。另外,比较表1和表2中的结果可见,表1中多数站位(即南黄海的南部)的非线性系数α在涨潮时段的值(绝对值)大于表2中各相对应站位(即南黄海的北部)的值,二者相差约2～3倍。

3 数值模拟结果

当使用方程(1)进行模拟时,参照赵俊生等[15]的观测结果,并依据文献[5]中的式(7)和(8),作者假定η(x,t)在各初始站位(AW,ANW,ASW,BW,BNW,BSW)的初始波形均为具有振幅5.5 m的余弦波形。

图10 第一模态半日内潮波的波形均为对应于最大振幅所在深度等密度面(即Φ(z)=1)的铅直位移η(x,t)Fig.10 The displacementη(x,t)of isopycnic of semidiurnal internal tidal wave at the dep th of the maximum ofΦ(z)(Φ(z)=1)for the first mode

图11 同图10,但是对应的站位和背景流情况不同Fig.11 As Fig.10,but the station and the background current case are difference

图12 同图10,但是对应的站位和背景流情况不同Fig.12 As Fig.10,but the station and the background current case are difference

图13 同图10,但是对应的站位和背景流情况不同Fig.13 As Fig.10,but the station and the background current case are difference

本文的数值模拟工作为诊断性的数值模拟研究。假定内潮波自初始站位开始的历经4个站位的非线性演变过程分为如下4种不同的背景流情况:(1)各站位均为涨潮流状态(相对于朝鲜半岛的西海岸而言);(2)各站位均为落潮流状态;(3)第1和第3站位为落潮流状态,而第2和第4站位为涨潮流状态;以及(4)第1和第3站位为涨潮流状态,而第2和第4站位为落潮流状态。为简短起见,本文仅给出第一模态半日内潮波从A站位向西对应于以上4种不同背景流情况的非线性演变过程的模拟结果,分别示于图10～13。从A站位向西北和西南的非线性演变过程的模拟结果与之类似。此外,本文给出第一模态半日内潮波从B站位向西对应于上述(1)和(2)2种背景流情况的非线性演变过程的模拟结果,分别示于图14和图15,其它情况的模拟结果与之类似。这些图中给出的内波的波形均为对应于最大振幅所在深度等密度面(即Φ(z)=1)的铅直位移η(x,t)。

与表1和表2中显示的非线性系数α的特性相一致,图10～13表明,在南黄海的南部海域,自A站位向西传播的内潮波在各站位均为涨潮流状态和AW 3站位为涨潮流状态出现显著的非线性裂变,并且产生内孤立波,而在其它背景流情况基本上为线性内波状态,仅出现轻微的非线性变化。图14和图15表明,在南黄海的北部海域,自B站位向西传播的内潮波即便在各站位均为涨潮流状态也不出现显著的非线性裂变。

图14 同图10,但是对应的站位和背景流情况不同Fig.14 As Fig.10,but the station and the background current case are difference

图15 同图10,但是对应的站位和背景流情况不同Fig.15 As Fig.10,but the station and the background current case are difference

4 结论

使用GKdV模型,通过对第一模态半日内潮波分别从南黄海的南部海域的A站(34°4′N,125°6′E)和南黄海的北部海域的B站(36°34′N,123°51′E)向西、西南和西北方向的反射和非线性演变过程的数值模拟研究,可以得到如下结论:

(1)在宽阔陆架海域,在内潮波和内孤立波的传播过程中背景流场的作用是十分重要的,并且区分背景正压落潮流与涨潮流的不同作用是必须的,这是GKdV模型的明显优势;

(2)南黄海的内孤立波集中分布于南黄海的南部,而在其北部极少出现内孤立波的主要原因是在南黄海的南部存在较强的背景斜压环流和很强的背景正压潮流,而且在涨潮时段(相对于朝鲜半岛的西海岸而言)背景斜压环流和背景正压潮流的方向基本相同,而在南黄海的北部这两种背景流都很弱;

(3)在内潮波和内孤立波的传播过程中,非线性系数α与水深变化,层化状态以及背景流场等多种因子相关,并非仅仅取决于内波振幅的大小;

(4)由于出现高达0.75 m/s的正压潮流,自A站位向西北方向传播的半日内潮波处于不稳定状态。

感谢:本文作者感谢中国科学院大气物理研究所刘海龙博士和国家海洋局第二海洋研究所张自历副研究员提供的帮助。

[1] 于万春,范植松,司宗尚.南黄海非线性内波的数值模拟研究.(1)内潮的内Kelvin波模型[J].中国海洋大学学报:自然科学版,2011,41(4):28-32.

[2] Fan Z S,Zhang Y L,Song M.A study of SAR remote sensing of internal solitary waves in the north of the South China Sea:I.Simulation of internal tide transformation[J].Acta Oceanologica Sinica,2008,27(4):39-56.

[3] 石新刚,范植松,李培良.正压潮流对南海东北部深水海域大振幅内孤立波影响的数值模拟[J].中国海洋大学学报:自然科学版,2009,39(sup.II):297-302.

[4] Holloway P E,Pelinovsky E,Talipova T,et al.A nonlinear model of internal tide transformation on the Australian North West Shelf[J].J Phys Oceanogr,1997,27:871-896.

[5] Grimshaw R,Pelinovsky E,Poloukhina O.Higher-order Korteweg-de Vriesmodels for internal solitary waves in a stratified shear flow with a free surface[J].Nonlinear Processes in Geophysics,2002,9:221-235.

[6] Holloway P E,Pelinovsky E,Talipova T.A generalized Korteweg-de Vries model of internal tide transformation in the coastal zone[J].J Geophys Res,1999,104:18333-18350.

[7] Liu A K,Holbrook J R,Apel JR.Nonlinear internal wave evolution in the Sulu Sea[J].J Phys Oceanogr,1985,15:1613-1624.

[8] Liu A K.Analysis of nonlinear internal waves in the New York Bight[J].J Geophys Res,1988,93(C10):12317-12329.

[9] Sandstrom H,Oakey N S.Dissipation in internal tidesand solitary waves[J].J Phys Oceanogr,1995,25:604-614.

[10] 刘海龙,俞永强,李薇,等.LASG/IAP气候系统海洋模式(L ICOM 1.0)参考手册[M].北京:科学出版社,2004:107.

[11] Liu H L,Li W,Zhang X H.Climatology and variability of the Indonesian Throughflow in the eddy-permitting oceanic GCM[J].Advances in Atmospheric Sciences,2005,22(4):496-508.

[12] 杜凌.全球海平面变化规律及中国海特定海域潮波研究[D].青岛:中国海洋大学海洋环境学院,2005.

[13] Shi X G,Fan Z S,Liu H L.A numerical calculation method of eigenvalue problem of nonlinear internal waves[J].Journal of Hydrodynamics,Ser.B,2009,21(3):373-378.

[14] 赵俊生等.浅海内波研究[M].中国海洋学文集(3),北京:海洋出版社,1992:102.