海气界面的动量通量即界面间动量的垂向输送, 反映了海洋与大气之间的摩擦拖曳作用, 直接控制和影响着海气边界层过程和海气边界层结构, 在海气相互作用、海洋-大气耦合模式、表面波与流的预报中起着重要的作用。目前在很多海洋环流的基础模式中, 风的作用对于海水的运动非常重要, 准确估算风应力是确定海气动量通量以及海浪与海洋环流模式驱动力的基础。如果利用涡相关法计算风应力, 需要观测不同时间和不同空间的湍流脉动情况, 观测存在较大的困难, 因此在实际应用中, 通常用块体公式来计算海面风应力, 利用海面拖曳系数或者海面空气动力学粗糙度对海面风应力参数化(邹仲水等, 2014)。

关于风应力的研究, 早期, 人们认为海气界面的动量通量只受到湍流剪切应力影响, 如今, 大量研究表明, 发生在海气界面的表面波对海气动量交换也有重要影响(Janssen, 1989; Belcher et al, 1993; Edson et al, 1998; Hristov et al, 1998, 2003; Sullivan et al, 2004; Veron et al, 2008)。现场观测(Davidson et al, 1973; Benilov et al, 1974; Smedman et al, 1994, 1999; Donelan et al, 1997; Rutgersson et al, 2001; , Semedo et al, 2009)和大涡模拟(Sullivan et al, 2008)表明, 通常情况海气界面动量通量有两个主要分量, 一是湍流剪切应力, 另一个是波致应力。当波浪传播速度小于风速时, 风为波浪提供能量, 但当波浪传播速度大于风速时, 波浪为风提供能量。在弱风条件下, 快速运动的波浪将诱导一个向上的动量通量, 这时波致应力将起主导作用, 总动量从水面指向大气, 这个由海洋到大气的向上动量输送将直接影响和修正海上大气边界层的流体动力学和湍流结构:如:近海面射流的生成(Miller et al, 1999), 平均风速随高度增加的减少(Rutgersson et al, 2001), 表面应力与风向的不一致(Grachev et al, 2003)等。事实上, 由于近表面气流的加速, 垂直风廓线不再满足对数廓线, 一个近表面波驱动风速最大值(a wind jet)将会产生。这种现象最早由Harris(1966)在实验室发现, 并称其为波驱动风(wave-driven wind)。前人研究表明, 对于涌浪, 部分波动量将以波驱动风的形式从海洋输运到大气中, 这时Monin-Obukhov相似理论不再适用, 涌浪对整个大气边界层将产生一定的影响, 这与通常人们认为的波边界层深度很浅(如: Makin et al, 1996)相矛盾。为了揭示涌浪对近海面风速廓线的影响, 人们提出了MABL海洋大气边界层模型(Makin, 2008; Hanley et al, 2008; Semedo et al, 2009)。Kudryavtsev等(2004)基于两层波边界层近似, 提出了一个涌浪与气流相互作用模型。他们利用该模型再现了Harris(1966)的实验室实验结果, 并解释了波驱动风的形成机制, 但是该模型只在波边界层范围内有效。基于波修正的Ekman模型, Makin(2008)通过推广混合长度的定义, 并在动量方程中考虑了科氏力的影响, 将Kudryavtsev等(2004)的模型推广到整个大气边界层。最近, Song等人(2015)基于Makin(2008)以及Hanley等人(2008)给出的波修正Ekman模型, 研究了表面波对稳定的海洋大气边界层(ABL)近海表风廓线的影响。他们比较分析了考虑波致应力与否模型解析解之间的差异, 并将解析解与波罗的海岛屿 stergarnsholm观测塔的资料进行了比较, 研究表明, 表面波不仅对近海面平均风廓线, 而且对海洋大气边界层的湍流结构均有重要的影响, 因为它们影响和改变了海洋大气边界层的结构。

虽然以上研究从不同角度揭示了波浪对海气界面动量通量以及海洋上方大气边界层的影响, 但仍有许多重要的科学问题不够清晰, 值得进一步探讨。总而言之, 如何准确计算拖曳系数(即风应力)一直是一个比较复杂的问题。本文的研究目的是比较波致应力的不同计算方法, 同时探讨波致应力在太平洋的季节性变化。波致应力的准确计算对反演拖曳系数具有重要的意义, 进而对风应力即海气界面动量通量交换的研究具有重要意义。

1 波致应力的计算方法

海气界面的动量通量即海面风应力一般分为两部分, 湍流剪切应力和波致应力(徐德伦等, 2001):

式中, τ为风应力; τ_t为湍流剪切应力, 根据普朗特混合长假设的理论可以参数化表示, 本文不做详细介绍; τ_w为波致应力, 正值表示风产生的进入波浪的动量, 负值表示波浪诱导的进入风的动量, 波致应力定义为(Markin et al, 1995; Semedo et al, 2009; Song et al, 2015):

式中, ρ_w是海水密度; z是距海平面的距离; ω是波浪的角频率; k是波数; K=(k_x, k_y)是水平方向的波数矢量; E(k, θ)是波浪波数方向谱; β_g是波浪增长率。

根据波致应力的计算公式, 波致应力与波浪增长率β_g和海浪谱E(k, θ)有关, 因此分别对波浪增长率公式和海浪谱计算公式进行敏感性分析, 并通过比较分析, 选择合适的波致应力计算方法。

1.1 波浪增长率公式的敏感性分析 1.1.1 不同的波浪增长率公式

已知波浪的波数方向谱E(k, θ), 波浪增长率β_g定义为:

不同的波浪增长率公式对波致应力的计算结果有重要影响, 这里给出4种常用的波浪增长率公式。

第一种, Komen等(1984)提到的波浪增长率公式:

式中, ρ_a、ρ_w分别是空气的密度和海水的密度; c是波浪的相速度; u_*是摩擦速度; θ是波向; θ_w是风向。u_*可以通过拖曳系数以及10m风速u₁₀进行计算:

C_d是在10米高度的拖曳系数, u₁₀是10米高度风速, 该文采用Zijlema等(2012)给出拖曳系数计算公式:

式中, ũ=u₁₀/u_ref，u_ref=31.5m/s。

第二种, Janssen(1991)提出的波浪增长率公式:

c_β为Miles常数, 是无因次临界高度的函数:

式中, κ=0.41是冯卡曼常数, g是重力加速度, J=1.2是常数, z₀是海表面空气动力粗糙度(Andrea, 2004):

第三种, Belcher等(1993)提到的波浪增长率公式:

Stewart(1974)提到波浪增长率系数c_β为:

式中, κ是冯卡曼常数, 取0.41; k是波浪的波数; z₀是海表面空气动力粗糙度; c是波浪的相速度; u_*是摩擦速度。

第四种, 波浪增长率β_g的公式形式与式(10)相同, 但是波浪增长率系数c_β的计算公式不同, 选用的是Polnikov(2011)根据实验数据得到的经验公式:

式中, θ_w是风向; c是波浪的相速度; u_*是摩擦速度。

1.1.2 敏感性分析

首先作两点假设:一、风向沿着x轴正向, 波向与风向一致, 即θ=θ_w=0;二、只考虑深水情况, 有深水频散关系。

实验观测数据使用的是Plant(1982)用到的经典实验数据, 其中三角形表示Shemdin等(1967)的实验数据; 实心圆表示Larson(1975)的实验数据; “x”表示Wu等(1971、1979)的实验数据; 空心圆表示Snyder等(1981)利用固定传感器得到的实验数据; 方形表示Snyder等(1981)利用波跟踪传感器得到的实验数据。

计算波浪增长率时, 取定10m高度处的风速值, 作无因次波浪增长率β_g/f随反波龄u_*/c的变化曲线图。

图 1表明, 对于Komen(1984)的计算公式, 当波数较大, 即反波龄较大的时候, 误差较大, 使得总体计算结果有一定的误差; 对于Janssen(1991)的计算公式, 波数较大时与实验数据吻合较好, 但当波数较小的时候, 误差很大, 使得总体计算结果有一定的误差; 对于Polnikov(2011)的计算公式, 计算结果普遍高于测量值, 因此相对误差较大; 相比之下, Belcher等(1993)的计算公式虽然在反波龄较大时与实验观测数据的误差稍微较大, 但总体的吻合度较好。

将4种公式的计算结果与经典实验观测数据进行对比, 考虑相对误差(RE)的平均值和最大值以及均方根误差(RMS):

式中, y_i为公式的计算值, x_i为实验的数据值。

表 1给出了风速分别为5m/s、10m/s、15m/s、20m/s情况下, 不同波浪增长率公式与实验观测数据的误差比较。可以看出, 改变给定风速的大小, 不同公式的计算结果与实验观测数据的相对误差和均方根误差都有所改变, 当风速较小时(< 15m/s左右)Belcher等的计算公式与实验观测数据吻合度较好, 而当风速较大时(> 20m/s左右)Janssen的计算公式与实验观测数据吻合度较好。

本文使用的数据是欧洲中期数值预报中心(ECMWF)的ERA-interim数据, 根据再分析数据的特点, 10米风速中大多数为中低风速, 同时综合来看, 无论是低风速情况还是高风速情况, Belcher等的计算公式与实验测量值都较为吻合, 相对误差的最大值明显比其他公式小, 因此下文中计算波致应力时, 我们选用Belcher等的波浪增长率计算公式。

1.2 海浪谱敏感性分析 1.2.1 单波近似

为使波致应力的求解更加快速有效, 首先只考虑单频率波, 单波情况不考虑方向的改变, 设波浪的能量:

式中, a是波浪的振幅; 波浪的能量; 有效波高, 则振幅a与有效波高H_s存在以下关系:

由公式(2)、(15)可以得到单波近似的波致应力:

式中, τ_m是单波情况[波谱为(15)式]下的波致应力, 是一般情形(2)式的特例。假设:风向沿着x轴; 只存在单波, 波向θ=0;只考虑深水情况。则有:

对于单波近似, 取波浪的波数k=k_p, k_p为谱峰波数, 则:

根据公式(16)可以将振幅用有效波高进行代换, 考虑海表面的波致应力:

关于谱峰频率与波龄的关系, Elfouhaily等(1997)提到:

式中, Ω是波龄, Zhao等(2003)提到波龄与有效波高、风速的关系:

那么, 当已知10米风速u₁₀、有效波高H_s, 根据公式(20)、(21)、(22)以及波浪增长率的计算公式(10)、(11)就可以计算得到近似求解的波致应力。

1.2.2 JONSWAP谱和Elfouhaily谱

海浪波面位移的均值一般取为0, 海浪谱(功率谱)定义为波面位移自相关函数的傅里叶变换。海浪谱表示海浪能量在不同频率、波数和方向上的分布, 对海浪谱积分得到海浪的总能量。

Lucas(2015)提到JONSWAP频率方向谱的表达式:

式中, ω_p为谱峰频率; 当ω≤ω_p时, σ=0.07, 当ω>ω_p时, σ=0.09;只考虑深水情况, 即有; γ称为峰升高因子, 观测值介于1.5—6之间, 取均值3.3。

常数α与有效波高有关:

对于风浪状态, Toba(1972)提到有效波高的计算表达式:

式中:谱峰相速度c_p=g/ω_p, 谱峰频率随无因次风区和风速增大而减小; u_*为摩擦速度, 可利用公式(5)、(6)计算。

但是一般的JONSWAP谱是风浪低频谱, 对于高频部分没有考虑, 使得计算结果与真实海浪谱可能存在较大的差异。因此, 需要考虑另外一种包含低频部分和高频部分的海浪谱, 该谱由Elfouhaily等人在1997年提出(后文统称作Elfouhaily谱):

式中l和h分别表示低频和高频; B表示曲率谱(B=k³E)。

B_l是长波曲率谱(低频部分):

式中波浪的相速度c=ω/k;

另外, B_h是短波曲率谱(高频部分):

式中波龄Ω=c_p/u₁₀。Elfouhaily谱只适用于反波龄大于0.84(波龄小于1.2)的情况, 即Elfouhaily谱也是风浪谱。

1.2.3 敏感性分析

利用单波近似(20)、JONSWAP谱积分公式(23)以及Elfoufaily谱积分公式(26)计算波致应力, 分别控制10m风速变量、波龄变量, 比较三种不同海浪谱公式计算的波致应力随波龄、风速的变化情况以及相互间的差异。

图 2和图 3表明, 在风浪占优的情况下, 单波公式的近似程度并不理想, 计算结果相比JONSWAP谱积分公式和Elfoufaily谱积分公式的计算结果都要小得多, 这与单波公式只考虑谱峰频率有关, JONSWAP谱积分公式没有考虑高频影响, 在波龄、风速增大时计算结果相对Elfoufaily谱积分公式偏小因此当波龄小于1.2时, 使用Elfoufaily谱积分公式计算波致应力较为合理。当波龄大于1.2时, 海域涌浪占优, JONSWAP谱和Elfoufaily谱都是风浪谱, 而目前没有普遍认可的涌浪谱的表达式, 考虑到涌浪频谱通常比较窄, 可以用一个主波近似代替, 为了简单起见, 在涌浪情况下, 我们使用单波公式进行近似求解。

2 波致应力求解在太平洋的应用 2.1 数据介绍

本文使用的数据是欧洲中期数值预报中心(ECMWF)的ERA-interim数据, ERA-Interim项目于2006年启动, 为ECMWF之前的再分析[ERA-40 (1957—2002)]和ECMWF设想的下一代扩展再分析提供了桥梁, 同时做了一定的改变:从地面到0.1hpa高度垂向分为60级; T255球面谐波表示为基本动力学领域; 对于表面和其他网格场, 具有大致均匀79km间距的降低高斯网格。时间跨度为1979年至今, 采用的同化方法为4D-VAR, 近地面风速的时间分辨率为6h, 模式近地面风速高度为10m。再分析数据对高风速通常存在一定的低估, 对低风速存在一定的高估, 但整体影响不大。本文中使用到的数据主要有10m风速(包括U分量和V分量)、包含风浪和涌浪的有效波高以及平均波向。所有数据选用的是观测数据, 数据时间为2009年1月、4月、7月以及10月, 时间分辨率为6h, 空间分辨率为85.5°S—65.5°N、90°E—60°W, 数据网格为0.125°×0.125°。

2.2 太平洋波致应力的求解

依据上文的分析, 由于Elfouhaily谱和JONSWAP谱是风浪谱, 而单波公式可以计算涌浪占优情况下波致应力的大小, 因此波龄大于等于1.2时使用单波公式来近似求解波致应力。当波龄小于1.2时, 根据分析得出单波公式的计算误差较大, Elfouhaily谱积分公式的计算结果较好, 为更准确的计算不同海况下波致应力的大小, 使用风速和有效波高数据构造Elfouhaily谱.。对不同海况条件, 我们假设:

式中, S_0i为i点构造的谱能; α_i为i点的反演系数; S_i为i点原来的Elfouhaily谱能。有效波高与波能的关系:

式中, H_si是i点的有效波高, 根据公式(38)、(39)以及不同海域的风速、有效波高数据得到反演的系数α:

对不同的海况, 波龄小于1.2时, 利用Elfouhaily谱积分公式得到反演系数α, 波致应力变为:

因此, 本文计算太平洋波致应力之前, 对每一个网格点计算其反演系数, 将计算得出的反演系数与理论Elfouhaily谱积分公式相乘用来修正海浪谱公式, 使得计算得到的海浪谱更加接近真实海况。同时在计算过程中, 首先计算每一个网格点的波龄, 当波龄大于等于1.2时使用单波公式计算波致应力, 当波龄小于1.2时使用Elfouhaily谱积分公式计算波致应力。

2.3 波致应力分布与风场

前人的大量研究表明, 浪与风存在着复杂的相互影响的关系, 风时、风区、风速等都会对波浪场产生较大的影响。为分析和解释波致应力的分布与风场存在一定相关性, 选择2009年1月1日00:00时刻的风速和有效波高数据计算波致应力的大小, 得到波致应力在太平洋的总体分布情况。

波致应力大于0, 即风产生的动量通量进入波浪, 此时为风浪占优的情况。波致应力小于0, 快速运动的波浪诱导了一个向上的动量通量, 总动量从水面指向大气, 此时为涌浪占优的情况。图 4表明, 总体而言, 此时刻赤道地区的波致应力的绝对值较小, 接近0;此时刻西北太平洋和澳大利亚以南的波致应力的绝对值较大, 且正值的绝对值要大于负值的绝对值, 即风浪占优情况下的波致应力绝对值大于涌浪占优情况下的波致应力绝对值。并且可以注意到波致应力的正值较大值与波致应力的负值较小值基本会同时出现在一块海域, 多数为风场产生的涡旋处。

根据风场图可以发现, 风速较大(接近或大于20m/s)的海域对应的波致应力大多数为正值且绝对值较大, 说明此海域为风浪占优的情况, 即风产生的进入波浪的动量通量较大; 风速为中等大小(接近10 m/s)的海域对应的波致应力大多数为负值, 说明此海域为涌浪占优的情况, 即波浪相速度大于风速, 在弱风条件下, 快速运动的波浪诱导了一个向上的动量通量, 总动量从水面指向大气, 即为负值; 风速很小的海域对应的波致应力接近0, 此处风与浪都较弱, 风产生的进入海浪的动量通量以及波浪诱导的向上的动量通量都很小。

在西北太平洋和澳大利亚以南的东西两侧都存在较强的涡旋, 同时发现存在明显涡旋的这三处海域都有较大的波致应力正值和较小的波致应力负值, 即此处的风浪和涌浪现象都存在且较为明显。风浪离开风直接作用的水域于风区以外形成涌浪, 于风区以内, 由于风速骤然降低或风向骤然改变, 原来的风浪于当地也会变成涌浪(文圣常, 1984)。由此得出, 此处由于风向的骤然改变加上风速的骤然降低, 原来的风浪在当地变成涌浪, 即涡旋中风速大的地方对应波致应力为正值且绝对值大, 风浪占优; 风速小的地方对应波致应力为负值且绝对值大, 涌浪占优。

2.4 太平洋波致应力季节分布特征

上节中通过风场的特征解释和分析了波致应力分布的合理性, 此节中继续讨论波致应力分布的季节性变化特征。以2009年1月、4月、7月、10月代表春夏秋冬四个季节, 用它们的风速、有效波高数据来计算月平均波致应力, 并绘制月平均波致应力在太平洋的分布图, 以及使用四个月的风速数据绘制月平均风场图。

图 6表明, 2009年1月份太平洋波致应力的大小在不同海区有较明显的差异, 尤其是北太平洋和南太平洋的西风带海区, 波致应力明显强于其他海区, 西风带的风速较为强盛, 对应的波致应力也较强盛。

图 7表明, 2009年4月份的波致应力相比于1月份整体较弱, 这与4月整体的风场强度整体比较弱有关, 波致应力在南太平洋的西风带明显强于其他海区, 在北太平洋的西风带稍微强于其他海区, 与风场强度的分布比较一致。

图 8表明, 2009年7月份太平洋的波致应力大小分布存在较明显的差异, 南太平洋西风带的波致应力相比其他海域明显强盛, 与对应的风场强弱分布一致。

图 9表明, 2009年10月份太平洋的波致应力整体较强且分布差异明显, 同样南太平洋和北太平洋西风带的波致应力明显强盛。

分析比较2009年四个月份的波致应力分布得到:从空间分布上看, 太平洋西风带的波致应力都较为强盛, 在4月和7月, 南太平洋西风带的波致应力比北太平洋西风带的波致应力要强盛的多; 在时间变化上, 太平洋的波致应力分布整体在1月和10月的时候较为强盛, 波致应力大于0.25的海区面积最多, 同时波致应力的最大值也较大。

3 结论

波浪增长率公式对波致应力计算的敏感性分析表明, 相比Komen、Janssen和Polnikov的计算公式, Belcher and Hunt的计算公式虽然在反波龄较大时与实验数据的误差稍微增大, 但总体的吻合度最好; 海浪谱对波致应力计算的敏感性分析表明, 在低风速的情况下, 单波近似的计算结果与Elfouhaily谱积分公式的计算结果较为接近, 但是当高风速时, 单波近似和JONSWAP谱积分公式的计算结果明显小于Elfouhaily谱积分公式的计算结果, 当波龄小于1.2时, 单波公式的近似误差比较大, 当波龄大于等于1.2时, 单波公式可以用来计算涌浪占优情况下海域的波致应力大小。

波致应力的分布与风场存在较大的联系, 风速较大(接近或大于20m/s)的海域对应的波致应力大多数为正值且较大, 此海域为风浪占优, 即风产生的进入波浪的动量通量较大; 风速为中等大小(接近10 m/s)的海域对应的波致应力大多数为负值, 此时涌浪占优, 即波浪相速度大于风速, 在弱风条件下, 快速运动的波浪诱导了一个向上的动量通量, 总动量从水面指向大气; 风速很小的海域对应的波致应力接近0, 此处风与浪都较弱, 风产生的进入海浪的动量通量以及波浪诱导的向上的动量通量都很小。对于太平洋波致应力分布的季节性特征, 从空间分布上看, 西风带的波致应力较其他地区更为强盛, 而在时间变化上1月和10月为波致应力整体较为强盛的时期, 波致应力大于0.25的海区面积最多, 同时波致应力的最大值也较大。