叶 灿 成泽毅 高 宇 宋金宝 李 爽

地形和风速影响下的海气相互作用大涡模拟研究*

叶 灿 成泽毅 高 宇 宋金宝 李 爽①

(浙江大学海洋学院 浙江舟山 316021)

当水流经过海洋地形时, 水流的不稳定性会引起垂向混合并伴随大量湍流过程。针对传统海气耦合模式缺少在湍流尺度上讨论海洋地形与风速对海气相互作用影响的问题, 使用并行大涡模拟海气耦合模式(the parallelized large eddy simulation model, PALM)在5 m/s的背景风场下, 引入理想立方体地形, 对比有无地形的影响; 设置地形边长为, 高为3(其中大气部分高),与水深之比为/=1/2; 然后保持地形条件不变。设置5、10和15 m/s三种风速, 讨论风速对小尺度海气相互作用的影响。研究表明: 地形在大气部分减弱顺风向速度, 增强侧风向速度, 影响0～5的高度区域, 而对垂向作用较小; 无地形条件下湍流垂向涡黏系数m在-0.3时, 水深达到最大值0.024 m2/s, 有地形条件下m在-0.8时, 达到最大值为0.16 m2/s, 地形的存在使得上层海洋混合加强,m最大值增加1个数量级。随风速增大海洋和大气中的净热通量、淡水通量和浮力通量都相应增大, 在近海面处, 5 m/s和10 m/s风速下三个通量的数值接近, 当风速为15 m/s时净热通量和淡水通量相较于前者数值大小增加2倍, 浮力通量增加近3倍, 说明大风加剧了各通量在海表的交换; 海洋混合层中湍动能收支各项也响应风速的变化, 其中剪切项、Stokes剪切、耗散项随风速增大而增加, 且在区域-0.2～0变化明显, 在近海表面处剪切项、传输项、压力项和耗散项的值达到最大, 同时耗散项由传输项和剪切项平衡; 随风速增大,m达到最大值的深度基本一致为-0.8。

大涡模拟; 地形与风速; 海气通量; 湍流动能收支

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

围绕三维障碍物的湍流在自然界中很常见, 并出现在许多应用中。当水流经过三维障碍物时, 其不稳定性会导致在三维障碍物后方产生一系列中小尺度涡旋, 其中伴随着大量的湍流过程。舟山群岛是我国第一大群岛, 也是重要的渔场, 其岛屿众多, 水动力条件十分复杂(秦华伟等, 2014; 刘紫薇等, 2022), 影响渔场及其周围生态。舟山市嵊泗县海域目前拟投资建设中广核嵊泗5#、6#海上风电场工程, 工程的建设会大大影响附近海域的水动力过程(Menéndez- Vicente, 2023; Syed, 2023)。了解和分析此类地形条件对海洋以及大气的影响可以推动对海气相互作用的相关研究。由于海洋条件的复杂危险性, 相关海洋数据十分不易获得, 而且通常获取的数据不够详细。高分辨率卫星观测资料的大量增加为人们认识中纬度中小尺度的海气相互作用提供了条件, 但存在分辨率不够的问题, 研究表明分辨率低于1°的观测数据对海洋锋和海洋涡旋现象等描述不足, 严重低估了中小尺度上的感、潜热通量变化, 进而可能低估了海表热通量对热带外气旋活动的影响(Rouault, 2003)。虽然观测数据逐渐由单点向区域和全球尺度延展, 但是仍无法覆盖全部(McClain, 2009), 同时也无法完整体现海洋现象的年代际变化等特征(Capotondi, 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)以Risf实验室开展的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)～(4), 利用张量伸缩运算, 通过在时间差分之前将动量预算项乘以适当的扰动速度场, 推出空间水平平均的次网格湍动能方程(subgrid- scale tuebulence kinetic energy, SGS-TKE) (Skyllingstad, 2000):

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

表1 参数列表

Tab.1 List of parameters

1.2 海气耦合模式简介

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

1.3 模式设置

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

图1 初始温盐剖面

注: a: 盐度剖面; b: 温度剖面

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

表2 对比实验设置1

Tab.2 Setting of the control experiment 1

注: EXP_O: 无地形条件; EXP_T: 有地形条件

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

表3 对比实验设置2

2 有无地形对小尺度海气相互作用的影响

2.1 流场变化

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

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

图3 区域水平平均的速度剖面图

注:为水深,为理想立方体长度; EXP_O: 无地形条件; EXP_T: 有地形条件; 背景风向沿水平方向

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

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

图4 10 m高度瞬时风速图(大气部分)

注: a:; b:; c:; 空白方块表示理想立方体地形

图5 水下10 m瞬时流速图(海洋部分)

注: a:; b:; c:; 空白方块表示理想立方体地形

2.2 海气通量

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

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

图6 净热通量、淡水通量和浮力通量的垂向剖面图

注: a: 大气和海洋的净热通量; b: 大气的淡水通量, 海洋中用盐度通量代替; c: 大气和海洋的浮力通量

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

2.3 垂向涡黏系数

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

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

图7 区域水平平均的垂向涡黏系数

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

3 有地形条件下风速对小尺度海气相互作用的影响

3.1 垂向动能时间序列

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

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

图8 垂向动能时间序列

注: a: 大气的垂向动能; b: 海洋的垂向动能

3.2 流场分布

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

图9 速度剖面图

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

3.3 瞬时涡度

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

图10 大气和海洋部分速度w分量的水平截面

注: a～b: 5 m/s; c～d: 10 m/s; e～f: 15 m/s; 空白方块表示理想立方体地形

图11 不同风速下水下5 m瞬时涡度分布图

注: a～b: 5 m/s; c～d: 10 m/s; e～f: 15 m/s; 空白方块表示理想立方体地形

3.4 海气通量

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

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

图12 净热通量、淡水通量和浮力通量的垂向剖面图

注: a: 大气和海洋的净热通量; b: 大气的淡水通量, 海洋中用盐度通量代替; c: 大气和海洋的浮力通量

3.5 湍流动能收支

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

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

3.6 垂向涡黏系数

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

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

图13 不同风速下的湍流动能收支图

注: a: 5 m/s, b: a的部分放大图; c: 10 m/s, d: c的部分放大图; e: 15 m/s, f: d的部分放大图

图14 不同风速下垂向涡黏系数剖面

4 结论与展望

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

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

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

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

孙丹译, 李爽, 2020. 基于大涡模拟耦合模式的小尺度海气相互作用研究[J]. 海洋与湖沼, 51(6): 1310-1319.

刘紫薇, 赵帅康, 魏笑然, 等, 2022. 主流海表风场资料在舟山群岛近海的性能评估[J]. 海洋与湖沼, 53(3): 528-537.

李德顺, 董彦斌, 傅学刚, 等, 2021. Bolund岛对风场流动影响的数值模拟研究[J]. 液压气动与密封, 41(5): 6-9.

吴凡, 2022. 南海上混合层结构表征[D]. 连云港: 江苏海洋大学: 1-2.

陆宗泽, 2018. 风应力对海洋上层垂向混合影响的湍流大涡模拟[D]. 杭州: 浙江大学: 33-34.

胡帅, 向月, 沈晓东, 等, 2021. 计及气象因素和风速空间相关性的风电功率预测模型[J]. 电力系统自动化, 45(7): 28-36.

秦华伟, 蔡真, 周红伟, 等, 2014. 舟山特定海域三维水流数值模拟研究[J]. 机电工程, 31(4): 541-544.

BOUSSINESQ J, 1877. Essai sur la théorie des eaux courantes[M]. Paris: Imprimerie nationale, 663-666.

CAPOTONDI A, WITTENBERG A T, NEWMAN M,, 2015. Understanding ENSO diversity [J]. Bulletin of the American Meteorological Society, 96(6): 921-938, doi: 10.1175/ BAMS-D-13-00117.1.

ESAU I, 2014. Indirect air–sea interactions simulated with a coupled turbulence-resolving model [J]. Ocean Dynamics, 64(5): 689-705.

FURUICHI N, HIBIYA T, 2015. Assessment of the upper-ocean mixed layer parameterizations using a large eddy simulation model [J]. Journal of Geophysical Research: Oceans, 120(3): 2350-2369, doi: 10.1002/2014JC010665.

GAO Z, LONG S M, SHI J R,, 2023. Indian Ocean mixed layer depth changes under global warming [J]. Frontiers in Climate, 5: 1112713, doi: 10.3389/fclim.2023.1112713.

GOPALAKRISHNAN S G, ROY S B, AVISSAR R, 2000. An evaluation of the scale at which topographical features affect the convective boundary layer using large eddy simulations [J]. Journal of the Atmospheric Sciences, 57(2): 334-351, doi: 10.1175/1520-0469(2000)057<0334:AEOTSA>2.0.CO;2.

HE J M, SONG C C S, 1997. A numerical study of wind flow around the TTU building and the roof corner vortex [J]. Journal of Wind Engineering and Industrial Aerodynamics, 67/68: 547-558, doi: 10.1016/S0167-6105(97)00099-8.

JIA W T, SUN J L, ZHANG W M,, 2023. The effect of boreal summer intraseasonal oscillation on mixed layer and upper ocean temperature over the South China Sea [J]. Journal of Ocean University of China, 22(2): 285-296.

LARGE W G, MCWILLIAMS J C, DONEY S C, 1994. Oceanic vertical mixing: a review and a model with a nonlocal boundary layer parameterization [J]. Reviews of Geophysics, 32(4): 363-403, doi: 10.1029/94RG01872.

LI S, LI M, GERBI G P,, 2013. Roles of breaking waves and Langmuir circulation in the surface boundary layer of a coastal ocean [J]. Journal of Geophysical Research: Oceans, 118(10): 5173-5187, doi: 10.1002/jgrc.20387.

LU Z Z, FAN W, LI S,, 2018. Large-eddy simulation of the influence of a wavy lower boundary on the turbulence kinetic energy budget redistribution [J]. Journal of Oceanology and Limnology, 36(4): 1178-1188.

MCCLAIN C R, 2009. A decade of satellite ocean color observations [J]. Annual Review of Marine Science, 1: 19-42, doi: 10.1146/annurev.marine.010908.163650.

MELLOR G L, YAMADA T, 1982. Development of a turbulence closure model for geophysical fluid problems [J]. Reviews of Geophysics, 20(4): 851-875, doi: 10.1029/RG020i004p00851.

MENÉNDEZ-VICENTE C, LÓPEZ-QUEROL S, BHATTACHARYA S,, 2023. Numerical study on the effects of scour on monopile foundations for offshore wind turbines: the case of Robin Rigg wind farm [J]. Soil Dynamics and Earthquake Engineering, 167: 107803, doi: 10.1016/j.soildyn.2023. 107803.

NOH Y, MIN H S, RAASCH S, 2004. Large eddy simulation of the ocean mixed layer: the effects of wave breaking and Langmuir circulation [J]. Journal of Physical Oceanography, 34(4): 720-735, doi: 10.1175/1520-0485(2004)034<0720: LESOTO>2.0.CO;2.

PELLICHERO V, SALLÉE J B, SCHMIDTKO S,, 2017. The ocean mixed layer under Southern Ocean sea-ice: seasonal cycle and forcing [J]. Journal of Geophysical Research: Oceans, 122(2): 1608-1633, doi: 10.1002/ 2016JC011970.

PHILLIPS O M, 1980. The Dynamics of the Upper Ocean [M]. 2nd ed. Cambridge: Cambridge University Press, 331-332.

REICHL B G, ADCROFT A, GRIFFIES S M,, 2022. A potential energy analysis of ocean surface mixed layers [J]. Journal of Geophysical Research: Oceans, 127(7): e2021JC018140, doi: 10.1029/2021JC018140.

ROUAULT M, REASON C J C, LUTJEHARMS J R E,, 2003. Underestimation of latent and sensible heat fluxes above the Agulhas Current in NCEP and ECMWF analyses [J]. Journal of Climate, 16(4): 776-782, doi: 10.1175/1520- 0442(2003)016<0776:UOLASH>2.0.CO;2.

SHAH K B, FERZIGER J H, 1997. A fluid mechanicians view of wind engineering: large eddy simulation of flow past a cubic obstacle [J]. Journal of Wind Engineering and Industrial Aerodynamics, 67/68: 211-224, doi: 10.1016/S0167-6105 (97)00074-3.

SKYLLINGSTAD E D, SMYTH W D, CRAWFORD G B, 2000. Resonant wind-driven mixing in the ocean boundary layer [J]. Journal of Physical Oceanography, 30(8): 1866-1890, doi: 10.1175/1520-0485(2000)030<1866:RWDMIT>2.0.CO;2.

STEINHORN I, 1991. Salt flux and evaporation [J]. Journal of Physical Oceanography, 21(11): 1681-1683, doi: 10.1175/ 1520-0485(1991)021<1681:SFAE>2.0.CO;2.

SULLIVAN P P, MCWILLIAMS J C, 2010. Dynamics of winds and currents coupled to surface waves [J]. Annual Review of Fluid Mechanics, 42(1): 19-42, doi: 10.1146/annurev- fluid-121108-145541.

SYED A H, MANN J, PLATIS A,, 2023. Turbulence structures and entrainment length scales in large offshore wind farms [J]. Wind Energy Science, 8: 125-139, doi: 10. 5194/wes-8-125-2023.

UMLAUF L, BURCHARD H, 2003. A generic length-scale equation for geophysical turbulence models [J]. Journal of MarineResearch, 61(2): 235-265, doi: 10.1357/002224003322005087.

VON STORCH J S, LÜSCHOW V, 2023. Wind power input to ocean near-inertial waves diagnosed from a 5-km global coupled atmosphere-ocean general circulation model [J]. Journal of Geophysical Research: Oceans, 128(2): e2022JC019111, doi: 10.1029/2022JC019111.

WALKO R L, COTTON W R, PIELKE R A, 1992. Large-eddy simulations of the effects of hilly terrain on the convective boundary layer [J]. Boundary-Layer Meteorology, 58(1/2): 133-150.

WILCOX D C, 1988. Reassessment of the scale-determining equation for advanced turbulence models [J]. AIAA Journal, 26(11): 1299-1310, doi: 10.2514/3.10041.

WOOD N, 2000. Wind flow over complex terrain: a historical perspective and the prospect for large-eddy modelling [J]. Boundary-Layer Meteorology, 96(1/2): 11-32.

LARGE EDDY SIMULATION OF AIR-SEA INTERACTION UNDER THE INFLUENCE OF TOPOGRAPHY AND WIND SPEED

YE Can, CHENG Ze-Yi, GAO Yu, SONG Jin-Bao, LI Shuang

(Ocean College, Zhejiang University, Zhoushan 316021, China)

When water flows through ocean terrain, vertical mixing and a large amount of turbulent may occur due to water flow instability. Aiming at the issue that the traditional air-sea coupling model lacks of discussing the impact of ocean topography and wind speed on air-sea interaction on turbulent scale, the Parallelized Large Eddy Simulation Model (PALM) was used to introduce an ideal cube topography under background wind field of 5 m/s. The impact of topography was simulated under three wind speeds of 5, 10, and 15 m/s. In the simulation, the side length of the cube was 1L and the height was 3of which the height above water surface was 1, and the ratio ofand water depthwas 1/2. The impact of wind speed on small-scale air-sea interactions was discussed. Results show that topography could weaken the downwind velocity in the atmosphere, enhances the crosswind velocity, and affect the altitude range of 0～5, while the impact on the vertical direction is small. Under no-topography conditions, the turbulent vertical eddy viscosity coefficientmreached the maximum value of 0.024 m2/s at-0.3water depth without topography, while the maximum value of 0.16 m2/s at-0.8with topography. The presence of topography enhanced the mixing of the upper ocean and increased the maximum value ofmby one order of magnitude. As the wind speed increases, all the net heat flux, freshwater flux, and buoyancy flux in the ocean and atmosphere increased accordingly. At near the sea surface level, the values of the three fluxes under 5 m/s and 10 m/s wind speeds are similar, while at 15 m/s, the net heat flux and freshwater flux increased by two times compared to the former, and the buoyancy flux increased by nearly three times, indicating that strong winds could intensify the exchange of various fluxes at the sea surface. In addition, the turbulent kinetic energy budget in the mixed layer of the ocean responded to the changes of wind speed. The shear term, the Stokes shear term, and the dissipation term increased with the wind speed increase. The variation was significant in the region of-0.2～0, and the values of shear term, transmission term, pressure term and dissipation term reached the maximum near the sea surface. Meanwhile, the dissipation term was balanced by the transfer term and shear term. As the wind speed increased, the depth at whichmreached its maximum value was-0.8in overall.

large eddy simulation; topography and wind speed; air-sea flux; turbulent kinetic energy budget

* 国家自然科学基金项目，41830533号, 41876003号。叶 灿，硕士研究生, E-mail: yecan@zju.edu.cn

李 爽, 博士生导师, 副教授, E-mail: lshuang@zju.edu.cn

2023-04-23,

2023-06-20

P731.26

10.11693/hyhz20230400090