海气相互作用影响着海洋生态要素的空间分布, 是影响海洋生态系统的重要过程, 是当前海洋科学的研究热点。而海洋上混合层作为一个中间层联系起大气底边界层和上层海洋, 是海气热量通量和动量交换的重要界面, 对海洋上混合层的热力学和动力学结构展开的研究已引起越来越多的关注(Phillips, 1980; Sullivan et al, 2010; Reichl et al, 2022; Gao et al, 2023; Jia et al, 2023)。湍流的生成及其垂向混合形成并维持了海洋上混合层, 作为中间层其同时受到多种时空尺度海洋和大气运动过程的影响, 因此具有复杂的热力学和动力学结构特征(吴凡, 2022; Pellichero et al, 2017; Von Storch et al, 2023)。同时海气界面动量通量也称为风应力, 是海流和表面海浪的主要驱动力, 是海洋从大气获得动量的重要途径。

围绕三维障碍物的湍流在自然界中很常见, 并出现在许多应用中。当水流经过三维障碍物时, 其不稳定性会导致在三维障碍物后方产生一系列中小尺度涡旋, 其中伴随着大量的湍流过程。舟山群岛是我国第一大群岛, 也是重要的渔场, 其岛屿众多, 水动力条件十分复杂(秦华伟等, 2014; 刘紫薇等, 2022), 影响渔场及其周围生态。舟山市嵊泗县海域目前拟投资建设中广核嵊泗5#、6#海上风电场工程, 工程的建设会大大影响附近海域的水动力过程(Menéndez- Vicente et al, 2023; Syed et al, 2023)。了解和分析此类地形条件对海洋以及大气的影响可以推动对海气相互作用的相关研究。由于海洋条件的复杂危险性, 相关海洋数据十分不易获得, 而且通常获取的数据不够详细。高分辨率卫星观测资料的大量增加为人们认识中纬度中小尺度的海气相互作用提供了条件, 但存在分辨率不够的问题, 研究表明分辨率低于1°的观测数据对海洋锋和海洋涡旋现象等描述不足, 严重低估了中小尺度上的感、潜热通量变化, 进而可能低估了海表热通量对热带外气旋活动的影响(Rouault et al, 2003)。虽然观测数据逐渐由单点向区域和全球尺度延展, 但是仍无法覆盖全部(McClain, 2009), 同时也无法完整体现海洋现象的年代际变化等特征(Capotondi et al, 2015)。随着超级计算机的出现, 利用数值模拟来研究这些相互作用已经成为可能。Walko等(1992)利用大涡模拟(large eddy simulation, LES)方法对比探究平坦地形和正弦丘陵地形之间可能存在的差异, 发现丘陵地形山顶附近上升运动的概率为70%, 山谷只有15%, 丘陵地形与平坦地形时变量的水平平均垂直廓线没有重要差别, 但该研究忽略了水平风的影响。He等(1997)用大涡模拟方法研究了单个建筑物周围的流场特征及建筑物顶部的涡旋结构。Shah等(1997)强调了大涡模拟方法用于立方体周围流场模拟的重要性。Wood(2000)在复杂地形流场模拟的回顾与LES前景展望中指出, LES在复杂地形流场模拟方面有独特的优越性, 但是也有很大的困难, 同时受计算量的限制。Gopalakrishnan等(2000)利用大涡模拟方法评估在哪些尺度上, 地形开始显著影响对流边界层(convective boundary layer, CBL)中湍流的平均特征和结构, 研究发现, 湍流与地形特征的尺度呈非线性关系。Lu等(2018)利用大涡模拟模式探究海底波状地形上湍动能收支的情况, 但其使用的是单一海洋模式, 未考虑海洋与大气之间的相互作用。李德顺等(2021)以Risφ实验室开展的Bolund岛实验结果为依据, 采用雷诺时均与大涡模拟方法对Bolund岛进行数值模拟, 将模拟区域方向坐标线为239°和270°的两条特征线命名为Line A和Line B, 研究结果表明: 针对Bolund岛的数值模拟, 在Line B特征线上5 m高度处模拟值比2 m高度处更接近实验值, 在岛前的位置模拟结果比岛后更准确, 同时说明LES方法能够有效揭示Bolund岛地形周围流动的绕流特征和尾流区涡特征。胡帅等(2021)考虑地形与风速的响应, 建立风速的空间相关性模型, 并将得到的风速数据转换为风电功率数据。

大涡模拟耦合模式作为一个好的方法已经开始被广泛应用于海洋和大气之间小尺度相互作用的研究中。Esau(2014)将大涡模拟的海洋模块与大气模块进行耦合, 探究了间接耦合机制在海气相互作用中的重要影响, 但未考虑风的作用。Noh等(2004)和Li等(2013)将Langmuir环流效应和波浪破碎效应加入了湍流大涡模拟中, 使得湍流大涡模拟在海洋上层的模拟更加接近实际海况; 孙丹译等(2020)使用大涡模拟海气耦合模式, 探究风速对小尺度海气相互作用的影响, 通过海气通量以及湍流动能收支进行表征; 但都缺少对有地形海域条件的分析。为此本文在前人的研究基础上使用并行大涡模拟海气耦合模式(the parallelized large-eddy simulation model, PALM)进行数值模拟, 通过流场分布、海气通量、湍流垂向涡黏系数和海洋湍流动能收支各项分布作为表征, 探究地形和风速对小尺度海气相互作用过程的影响。

1 模式简介及设置 1.1 PALM简介

本研究使用的是基于Boussinesq近似的非流体静力学数值模拟模式PALM, 对经过滤波的不可压缩Navier-Stokes方程进行求解。其基本方程如下公式(1)~(4)所示。

其中, 等式中尖括号表示水平域平均值, 上划线表示水平平均, 下标0表示海表面值, 双撇表示次网格尺度湍流模式变量。

在海洋模式: 将公式(1)中的 改为, 其中由海水状态方程得到; 考虑到盐度S对于密度的影响, 增加方程:

根据公式(1)~(4), 利用张量伸缩运算, 通过在时间差分之前将动量预算项乘以适当的扰动速度场, 推出空间水平平均的次网格湍动能方程(subgrid- scale tuebulence kinetic energy, SGS-TKE) (Skyllingstad et al, 2000):

公式(6)中等号右边各项依次是: 速度剪切项、浮力项、次网格耗散项、湍流动能传输项、压力项和耗散项。在表 1中列举出了上述方程中出现的所有参数及其说明。

1.2 海气耦合模式简介

PALM可以实现海气耦合模式, 讨论湍流过程在大气边界层(atmospheric boundary layer, ABL)和海洋混合层(oceanic mixed layer, OML)中的相互作用, 通过一个大气和一个海洋两个大涡模拟之间在海面上的各项变量实时在线交换实现耦合。根据水的密度来调整动量通量和从而保证海洋和大气之间的通量守恒, 如公式(8)和公式(9)所示:

海面的冷却由海水蒸发所导致, 海洋的热通量同时依赖于两个变量, 一个是大气的热通量, 另一个是大气在海表面的湿度通量, 如公式(10)所示:

公式(10)中表示水在恒压下的比热容。Steinhorn(1991)提出, 海洋部分盐分无法蒸发, 所以海水蒸发时会导致其含盐量增加, 这一过程可通过盐度通量表示, 在海面向下, 即为负值, 具体由公式(11)表示:

边界条件如式(12)所示:

其中, 表示海面位温, 和为水平速度分量。在PALM海气耦合模式中, 可以单独设置大气模式和海洋模式的时间步长, 不要求二者相等, 分别进行前驱运行后再以指定的频率执行耦合运行, 这样也可以减少计算量。在模式设置中, PALM要求大气部分和海洋部分的水平域大小相等, 因为耦合是通过海面数据的双向双线性插值实现。

1.3 模式设置

本文使用PALM海气耦合模式, 初始设置包括海洋和大气两个部分。考虑到在模式中加入理想立方体地形, 为了保证在模拟水域内可以充分讨论地形带来的影响, 设置足够大的模拟区域, 海洋部分网格数为640×320×100 (x, y, z)个, 网格大小是1.25 m× 1.25 m× 1 m, 即模拟区域大小为800×400×100 m³ (x, y, z); 大气部分的网格数是640×320×60 (x, y, z)个, 网格大小是1.25 m×1.25 m×5 m, 即模拟区域大小为800×400× 300 m³ (x, y, z)。PALM海气耦合可以分为前驱运行和耦合运行两个部分, 这两个部分的各项初始条件设置大致相同, 采用三阶Runge-Kutta方法为时间步长方案, 将侧边界条件设置为周期性边界条件, 同时海洋模式中的上边界和大气模式中的下边界均采用自由滑移边界条件。参照Esau (2014)的相关研究内容设置大气模式的表面温度是270.6 K, 温度梯度是0.015 K/m, 海洋模式温度是272.65 K, 盐度是34.8, 本文中纬度设置为45°N, 具体初始温盐剖面如图 1所示。

本文不考虑海洋分层和海洋背景流, 前驱运行时通过风速的输入使大气模式启动, 在海洋初始运行时引入一个微小的扰动使得海洋模式启动, 具体表现是输入模式中在海表以热通量形式加入, 模拟结果不会受该微小扰动所影响。前驱运行中设置运行时间大气模式和海洋模式均为28 800 s, 耦合运行的时间则为86 400 s, 耦合的时间步长均是30 s。本文首先通过在5 m/s的地转风背景下设置无立方体地形和有立方体地形两种实验来探究地形对流场分布、海气通量及湍流垂向涡黏系数的影响, 将无地形条件命名为EXP_O, 有地形条件命名为EXP_T。具体对比实验设置见表 2所示。

探究地形对海气相互作用的影响, 简单起见, 本文在耦合模式中加入理想立方体地形(图 2), 以立方体柱的形式模拟三维地形的存在。定义理想立方体柱长、宽为L, 高为3L (其中在海洋部分高占2L, 大气部分高占L), 模拟域总水深定义为H, H=2L。然后不改变地形条件, 设置三种地转风速度(5, 10和15 m/s)来讨论改变风速对小尺度海气相互作用的影响, 具体对比实验设置见表 3。

2 有无地形对小尺度海气相互作用的影响 2.1 流场变化

首先关注大气和海洋中速度分量在垂向上的变化, 将数据进行时间平均后得到大气风速和海洋流速三个分量u、v和w的垂向剖面图(图 3)。图 3中横坐标为速度大小, 黑线表示无地形条件EXP_O, 红线表示有地形条件EXP_T。为了便于展示, 本文将海洋部分的速度分量u、v数值放大100倍。

与EXP_O相比, EXP_T中在地形的作用下, 大气部分u方向风速明显减小, 模拟实验中设置立方体地形长宽以及在大气部分的高度相等, 与模拟域水深之比为L/H=1/2, 其对速度的影响高度可到5L, 随高度升高地形作用逐渐减弱, 风速有所回升, 并且在高度大于5L后u数值略大于无地形条件。海洋部分在海表面处EXP_T条件下速度最大为0.023 m/s, 相比于EXP_O条件u最大值为0.045 m/s减小一半; 随水深加深两种实验条件下速度分量u数值相接近。由于模式的地转风输入只在水平x方向, 所以在EXP_O中大气侧风向速度v随高度升高逐渐减小, 在高度大于5L时趋向于0, 表明在5 m/s的风速作用下对侧风向v的作用同样可达5L的高度。在EXP_T中v方向风速较之无地形情况下有所增大, 说明地形的存在可以降低顺风向速度, 增加侧风向速度。海洋部分海底v方向流速从0开始变化且都为正值, 同u一样两种实验条件下v数值也相近。速度分量w在大气部分两种实验条件下几乎相同且接近零, 表示低风速影响下大气垂向作用较弱, 值得注意的是, 在EXP_T中, 由于高1L立方体地形的存在, 在地形顶部w绝对值陡然反向增大然后恢复稳定变化趋势。在海洋中, 地形阻挡作用使w从总体上数值减小, 在海表由于地转风的作用, 随水深加深w逐渐增大, 当能量向下传递至逐渐耗散, w值开始减小直到遇到底边界数值为0。

图 4和图 5分别展示了大气部分在0.2L高度处和海洋部分在−0.2L深度处的瞬时大气风速(海洋流速) u、v、w的变化情况, 时间取模式最后输出时刻, 用EXP_T减去EXP_O, 反映加入地形后流场的增减情况。图中横纵坐标将理想立方体地形存在的中间位置处数值设为0。在初始设置中纬度为45°N, 即科氏力, f=10⁻⁴ s⁻¹由于Ekman-Stokes边界层动力学的影响, 所以在图 4中可以看到大气中的风向沿水平x方向顺时针偏转。在地形斜后方水平距离4L的区域内两个实验下风速分量u和v的差值为负, 因为在EXP_T中受地形影响, 在地形后方一块区域的风速约为0, 所以速度分量差值为负。图 4c中在地形斜后方风场w分量数值变化较为混乱, 表明此区域垂向混合剧烈。

同样地, 受地形影响, 在其后方一块区域的流速u约为0, 被称为“死水区”, 所以图 5a显示在地形后方, EXP_T和EXP_O流速差值为负值, 图中死水区的上下宽度和地形直径大体一致, 且流速在纵向上(水平y方向)向两侧方向增大, 一段距离后恢复正常, 在横向上(水平x方向)上约6L处流速逐渐恢复正常, 但是受风和地形的持续影响, 之后的一段水流变得混乱。对速度分量v的变化进行讨论, 从图 5b中可以看出, 在地形前方产生两股相反的水流, 是由于地形存在的绕流效应导致的偶极子型v, 在地形前方两端对称。图 5c中, 速度分量w在地形后方剧烈变化, 表示由于地形作用其后方垂向混合作用加强。

2.2 海气通量

PALM中的海气耦合是ABL和OML之间的耦合, 耦合是通过一个大气-海洋界面来保存每个单元网格内的质量(湿度、盐度)、动量和热量的湍流通量。以ABL底部的湍流通量为OML的上边界条件, 以OML顶部的速度和海表温度作为ABL的下边界条件。ABL在上一个时间步长获得表面湍流通量(热、动量、质量)用于设置OML顶部边界, 在同一时间步长中, OML为下一个时间步长设置ABL底边界的表面温度和速度。

通过模式所得数据对其进行时间平均得到图 6, 图中横坐标为通量。由图 6可以直观地看到两个实验下海气通量的垂向变化。在海气界面, 加入地形条件后, 大气中的净热通量、淡水通量都相应增大, 在同一高度下有地形条件时净热通量、淡水通量要更大; 而浮力通量的数值在0~2.4L的区域内大于无地形条件且为正值, 在2.4L~5.0L的区域内小于无地形条件且为负值, 在5L~6L的高度区域净热通量、淡水通量和浮力通量在两个实验条件下数值相近且趋向于0, 说明地形对大气部分的影响高度在0~5L, 这与图 3中对速度剖面进行分析得出的结果一致。

海洋部分在加入地形条件后, 净热通量和浮力通量都减小。而淡水通量(在海洋中用盐度通量代替)在−1L~0的深度区域内未加入地形条件时为负值, 加入地形条件后为正值, 表明地形存在使海表面受热蒸发减少; 在−2L~−1L的水深区域内则相反, 未加入地形条件下为正值且接近于0, 表明在水深到达−1L直到海底蒸发很弱, 在加入地形条件后为负值, 地形的存在给水深−1L到海底的区域内带来一定的(热量)盐度交换。同时加入地形后, 各通量在海洋向下传递的距离和在大气中向上传递的距离与在无地形条件下相差不大。

2.3 垂向涡黏系数

可通过湍流垂向涡黏系数来指示海洋的垂向混合过程。Boussinesq(1877)第一次提出关于湍流垂向涡黏度的假设, 其计算公式为

其中, 水平方向上速度二阶矩的平均值表示为, 平均速度沿垂向梯度的平均值为, 则表示水平方向上的湍流垂向涡黏系数。随后一段时间里, 湍流垂向混合参数化方案被不断完善。目前, 常用的参数化方案包括Large等(1994)提出的KPP (K-Profile Parameterization)模型、Wilcox(1998)提出并被Umlauf等(2003)完善的k-ε模型、Mellor等(1982)提出的MY2.5模型等等。

根据公式(13)计算得到海洋部分的垂向涡黏系数随深度变化如图 7所示, 图中横坐标为湍流垂向涡黏系数K_m数值。因为本节模式设置风速输入只沿水平x方向, 所以只研究沿该方向上的K_m变化。模式基本达到稳定状态的时间为运行开始后20 min, 因此取模式运行20 min后的水平平均值。

从图中可以看出, 在EXP_O垂向涡黏系数K_m的值在海表从0开始, 这是因为模式未考虑波浪破碎效应和海表通量的输入等条件; 然后随水深增大K_m的值逐渐增大, 当水深为−0.3L时达到最大值, 为0.024 m²/s, 是由于表面风场的作用使得垂向混合加剧; 接下来由于风作用向下传递的能量逐渐耗散, 所剩余能量不足以维持水体的混合, 所以随水深继续增加K_m值开始减小。与EXP_O相比较, 在EXP_T中加入立方体地形后整体K_m值的变化趋势保持不变, 即增大到最大值后再减小, 但在相同水深处K_m数值明显增大, 其达到最大值的深度也有所加深为−0.8L, 同时K_m最大值相比于EXP_O增加1个数量级, 为0.16 m²/s。说明在地形作用下, 海洋垂向混合增强。总体上K_m值均小于1, 根据公式(13)说明影响K_m变化的主要因素为垂向剪切。

3 有地形条件下风速对小尺度海气相互作用的影响 3.1 垂向动能时间序列

本节保持1.3节的初始条件不变, 改变风速, 设置5、10、15 m/s三种不同地转风速度, 在立方体地形存在条件下, 对不同风速下的小尺度海气相互作用情况进行分析。

图 8展示了在不同风速下大气和海洋部分的垂向动能时间序列, 图中横坐标为时间, 纵坐标为垂向动量通量, 为负值表示动量向下传递。从图 8中可以看出, 耦合模式下, 相较于海洋, 大气部分可以更快达到稳定状态, 在各风速下海洋模式约在3 h后均能达到稳定。大气模式中, 随风速增大动量在垂向上的传递逐渐加强, 在风速为5 m/s时, 对所有时间序列输出的垂向动量通量数值进行平均, 得到其平均值为−2.2×10⁻² m²/s²; 风速为10 m/s和15 m/s时, 垂向动能值上下波动。海洋模式中, 垂向动能的变化并不完全响应风速的增大, 在不同风速下达到稳定状态后其值上下波动变化。

3.2 流场分布

时间平均后得到大气风速和海洋流速分量u、v、w的垂向剖面如图 9所示, 图中横坐标为速度。图 9a显示在大气模式中, 同一风速下随高度升高u逐渐增大且增大趋势随高度上升渐缓, 同时随风速增大u逐渐增大。在图 9a海洋部分, 风对流速的作用在海表处较为明显, 风速为5 m/s时, 海表u最大值为0.023 m/s; 风速为10 m/s时, 海表u最大值为0.044 m/s增大将近1倍; 风速为15 m/s时, 海表u最大值为0.06 m/s。随水深加深, 风给速度分量带来的变化逐渐减小, 在不同风速条件下u数值接近。图 9b大气模式中, 速度分量v在10和15 m/s的地转风速下随高度升高风速增大且增大趋势逐渐变缓, 而在5 m/s风速下一开始随高度升高v缓慢增大, 在高度大于3L后又缓慢减小; 总体上随初始地转风风速增大v逐渐增大, 说明在低风速下风的作用带来的侧风向速度增大效果没有高风速下好。在海洋部分, 同u一样, v在近海面处变化迅速, 不同的是在同一水深处随风速增大v数值呈较明显的增大趋势。

从图 9c中可以看出, 在大气模式中垂向流速w总体上变化不明显, 但在高度为1L处, w值反向突增又迅速减小, 是由于大气模式中1L高的地形作用, 在地形顶部形成反向涡旋, 同时初始设置的地转风速度越大, w值突增越明显, 说明在大风速作用下风和地形的相互影响更强。风应力可以使海洋表层产生水平流动, 从而破坏混合层的稳定性, 使混合层中的水体发生垂向运动, 因此海洋混合层中的垂向速度比海表处的垂向速度要大得多。在风速作用下随水深增大w逐渐增大, 在约−0.3L处达到最大然后逐渐减小, 海底减至0, 说明风的作用深度大概为水下0.3L。为了更清晰地展示流场变化, 取模式最后输出时刻作出大气模式0.2L高度处和海洋模式−0.2L深度处的瞬时水平流速图(图 10), 用有地形条件的结果减去无地形条件的结果, 显示加入立方体地形后流速的增减情况。随风速增大, 大气部分在地形右后方出现明显的下降气流, 且影响范围逐渐增大, 当风速为15 m/s时, 影响可达地形后方3L距离。海洋部分, 在风速为5 m/s时地形后方受到影响水流显得混乱的距离约为6L, 随着风速的增大, 影响范围也逐渐增大, 在整个模拟水域内水流变得混乱。

3.3 瞬时涡度

图 11为在不同风速下水下−0.2L处的瞬时(模式最后输出时刻)涡度水平截面图, 图中白色区域表示立方体地形, 灰色区域代表涡度值为0。从图中可以看出在立方体地形后面形成大小近似相等方向相反的对称涡结构, 其涡度数上半为负值, 下半为正值。当风速为5 m/s时, 在水平x距离约为2.3L处涡度上下分布慢慢汇合, 形成完整的“涡街”结构; 风速为10 m/s时, 上下分布的涡汇合距离变短, 约为2L, 而风速为15 m/s时则在1.8L处形成完整的“涡街”结构, 可见风的增加对尾涡有一定聚拢作用。不同风速下, 立方体地形后方尾部的涡结构相似, 涡在水平y方向的移动范围也很相近, 但当风速增大可以影响的水平x距离更远, 涡度分布也更为分散。

3.4 海气通量

图 12展示了不同风速下海气通量的垂向变化, 图中横坐标为通量大小。总体上来看, 大气部分随着高度升高各通量值逐渐减小; 在海气界面处, 随着风速的增大, 大气中的净热通量、淡水通量和浮力通量也都相应增大; 图 12c中, 由于地形的存在, 在高度为1L处浮力通量陡然减小再迅速恢复至正常变化趋势, 且风速越大突变程度越大。

海洋中随水深由浅至深净热通量、淡水通量和浮力通量逐渐增大, 且随风速增大相同水深处各项通量值也相应增大, 在近海面处, 5和10 m/s地转风条件下, 其净热通量、淡水通量和浮力通量数值接近, 而当风速为15 m/s时净热通量和淡水通量相较于前者数值大小增加2倍, 浮力通量增加近3倍, 说明大风加剧了各通量在海表处的交换。对不同条件下各通量向下传递的距离进行分析, 风速越大各通量值向下减小的趋势越缓慢, 说明风速的增加会减缓各通量在海洋上层中的消耗, 这与孙丹译等(2020)的研究结果相同。

3.5 湍流动能收支

图 13是海洋部分在三种不同风速条件下的湍流动能收支图, 根据公式(6)计算得出各项值, 图中横坐标为湍动能收支, 红线代表速度剪切项, 天蓝线代表Stokes剪切项, 绿线代表传输项, 蓝线代表压力项, 黑线代表耗散项, 为了便于展示将各项数值×10⁶。左侧图为不同风速条件下海洋湍流动能收支情况, 右侧均为左侧在水深−0.4L~0范围的局部放大图。

从图 13a、13c和13e中可以发现, 海洋湍流动能收支各项约在−0.2L处迅速减小。在低风速条件下(即风速为5 m/s), 模式中地形的存在抑制了剪切项的生成, 所以剪切项数值较小, 在距离海面约−0.1L~0范围内剪切生成又迅速消失, 剪切项与耗散项在近海面处的值达到最大值且两项大致持平, 最大值约为4.35×10⁻⁸ m²/s³。耗散项表示摩擦引起的动能损失, 已有研究表明动量从近海表传递到整个海洋混合层, 大约有40%的输入能量通过湍流耗散损失(Skyllingstad et al, 2000)。剪切和耗散之间在近海面处的近乎平衡表明, 大部分的湍流能量集中在直接受次网格耗散影响的小规模流动特征中。当风速增大至15 m/s时, 剪切项明显增大, 地形抑制作用减弱, 其最大值为8×10⁻⁷ m²/s³。海洋混合层中湍流动能收支也响应风速的变化, 其中剪切项、Stokes剪切、耗散项随着风速增大而增加, 且在区域−0.2L~0变化明显, 由于海洋模式受到大气风速的驱动, 所以在近海表面处各项值都达到最大, 之后随深度增加逐渐减弱。同时在上层海洋中耗散项通过传输项和剪切项来平衡。

3.6 垂向涡黏系数

根据公式(13)计算得到海洋部分的垂向涡黏系数随深度变化如图 14所示, 图中横坐标为湍流垂向涡黏系数。同2.3节一样, 本节模式风速输入只有水平x方向, 所以只研究沿水平x方向的K_m变化, 取模式运行20 min后的水平平均值。

从图 14中可以看出, K_m的垂向分布趋势受风速影响, 呈现先增大后减小的趋势。对比不同风速条件可以看出, 随着风速的增大, 总体上在相同水深K_m值也增大; 当风速为15 m/s时, K_m最大值可达0.33 m²/s, 相较于风速为5 m/s时数值增大1倍。但是, 风速不同K_m值达到最大值的深度基本一致, 这与通过MY2.5等模式所得出来的结果有出入, 之前的结果显示风速越大K_m最大值的出现位置越向水域更深处移动(Furuichi et al, 2015), 但是在实验室实验以及数值模拟中有被捕捉到类似的结论。陆宗泽(2018)将MY2.5垂向湍流参数化的结果与湍流大涡模拟的结果进行了对比, 结果发现, 大涡模拟的结果和MY2.5的结果吻合得很好, 并且大涡模拟在湍动方面体现得更好, 大涡模拟区别于传统的垂向混合参数化方案, 直接计算出垂向湍流涡黏系数。文章设置4组不同的海表面10 m风速(10、15、20、25 m/s), 结果表明随风速增大垂向涡黏系数K_m达到最大值几乎保持在一个深度。同时K_m值均小于1, 在不同风速条件下垂向剪切依旧是影响K_m变化的主要因素。

4 结论与展望

本文通过大涡模拟海气耦合模式来探究地形和风速对小尺度海气相互作用的影响, 设置理想立方体地形, 地形长宽为L, 高为3L (其中大气部分高L), 与水深H之比为L/H=1/2, 在5 m/s背景风下, 对比有无地形的区别, 然后保持地形条件不变, 改变风速为5、10和15 m/s, 展示不同风速下大气风场和海洋流场、海气通量以及海洋混合层湍动能收支各项的分布。主要结论如下:

(1) 地形在大气部分减弱顺风向风速, 增强侧风向风速, 作用高度在0~5L的区域, 对垂向影响很小。在地形斜后方一小块区域风速分量u和v减小为0, w分量数值变化较为混乱, 表明垂向混合剧烈。在海洋部分水平流速变化不大, 在地形后面的一块区域出现“死水区”, 在横向上流速影响距离约6L。在海气界面, 加入地形条件后, 大气中的净热通量、淡水通量都相应增大, 海洋中则是净热通量和浮力通量减小, 淡水通量(用盐度通量表示)为正表明海面蒸发减少。同时相比于无地形条件下, 在有地形时K_m最大值增加1个数量级, 说明地形的存在使上层海洋混合加强。

(2) 随着风速的增加, 在大气部分, 垂向动能的分布呈现纵向加深的趋势, 水体混合程度加深; 在海气界面处, 随风速增大大气中的净热通量、淡水通量和浮力通量相应增大, 由于地形存在, 在高度为1L处浮力通量陡然减小后恢复正常变化趋势; 海洋中随水深由浅至深净热通量、淡水通量和浮力通量逐渐增大, 且随风速增大各项值也相应增大, 但在近海面处, 5和10 m/s地转风条件下, 其三个通量的数值接近, 当风速为15 m/s时净热通量和淡水通量相较于前者数值大小增加2倍, 浮力通量增加近3倍, 大风条件下加剧了各通量在海表的交换; 海洋混合层中湍流动能收支也响应风速的变化, 其中剪切项、Stokes剪切、耗散项随着风速增大而增加, 在区域−0.2L~0变化明显, 由于海洋模式受到大气风速的驱动, 所以在近海表面处几项值达到最大, 之后随深度增加逐渐减弱。同时在上层海洋传输项和剪切项来平衡耗散项。随着风速的增大, 总体上K_m值也增大, 风速为15 m/s时K_m最大值为0.33 m²/s, 是风速为5 m/s条件下K_m最大值的两倍。风速不同K_m值达到最大值的深度基本一致, 为−0.8L。

本文实现了基于大涡模拟的海气耦合模式, 对于海洋和大气之间通量交换等相互作用过程的研究有重要的参考意义。使用耦合模式充分考虑海洋部分与大气部分的协同作用, 探究地形与风速的影响, 模拟结果与真实情况更接近, 从而为真实海洋中岛屿存在处的海洋动力研究以及风电场等工程建设过程提供指导。但海洋混合层受风速影响存在两个重要湍流过程: 波浪破碎和Langmuir环流。波浪破碎在海面附近产生大量的小尺度湍流, 而由风驱动的表面流切变与Stokes漂移的相互作用产生Langmuir环流, 这两者对于海洋混合层有较大影响。Noh等(2004)对海洋混合层进行了波浪破碎和Langmuir环流的大涡模拟, 揭示了波浪破碎和Langmuir环流在海洋混合层垂向混合过程中的重要作用。Li等(2013)使用大涡模拟研究两者在海表边界层中产生湍流的作用, 模式结果与从CBLAST站点湍流测量分析得出的结论一致。波浪破碎与Langmuir环流均发生在海气边界上, 影响海气通量交换, 所以未来的工作重点便是在耦合模式中加入波浪破碎和Langmuir环流, 考虑二者对海气相互作用的影响。