渤海作为我国北部典型的半封闭陆架海, 是我国最早进行悬沙分布和输运研究的海域。前人首先进行了大量的悬沙浓度现场观测, 得到了渤海悬沙浓度的空间分布特征及其季节变化规律(秦蕴珊等, 1982; 中国科学院海洋研究所海洋地质研究室, 1985; Jiang et al, 2004)。通过与渤海水动力环境相结合, 能够大体明确悬沙输运的模式和机理。早期应用的数值模型不论是二维模型(Fang et al, 1985; Zhang et al, 1986; 董礼先等, 1989)还是三维模型(Feng, 1987; Graber et al, 1989; Zhao et al, 1993; 王凯等, 1999), 大都仅考虑潮流或风浪的作用。随着数值模拟技术的发展, 对渤海悬沙输运的数值模拟开始综合考虑浪-潮-流等水动力耦合过程对悬沙输运的影响(江文胜等, 2000; 庞重光等, 2004)。

渤海水动力环境的研究和模拟是悬沙输运数值模拟的基础。渤海水动力环境的研究自20世纪60年代以来开展至今, 其中基于各类数值模式的研究较多(Fang et al, 1985, 2000; Zhang et al, 1986; Feng, 1987; Zhao et al, 1993; Guan, 1994; 乔方利, 2012), 数值模拟结果揭示了渤海水动力环境的基本特征。Xia等(2019)利用高分辨率波浪-潮汐-环流耦合模式模拟渤海夏季的环流, 发现辽东湾西部至黄河口近岸海域沿温度锋面存在明显的密度流, 在40°N以南海域环流呈双环流结构。马伟伟等(2016)利用ROMS (regional ocean modeling system)海洋模式对渤海冬季流场进行高分辨率的数值模拟, 发现冬季深度平均流场在辽东湾为顺时针环流, 且在渤海中央表现为顺时针的环流结构。相比于研究较为成熟的渤海水动力特征, 渤海悬沙的研究大都集中在黄河口, 渤海海峡等热点区域(Bi et al, 2011; Yang et al, 2011; Wang et al, 2020)。受现场观测范围的限制, 整个渤海海域的悬沙分布研究一般通过卫星遥感数据的反演得到(庞重光等, 2013, 2014; 周舟等, 2017; Zhao et al, 2022)。庞重光等(2014)通过分析长时间序列悬沙浓度和风场遥感反演资料, 指出渤海不同海域表层悬沙浓度差别较大, 高浓度集中在辽东湾、渤海湾和莱州湾, 低浓度主要位于秦皇岛海域以及渤海中部等。庞重光等(2013)基于SeaWiFs反演的渤海表层悬沙浓度资料, 通过经验正交函数分解, 展示了渤海表层悬浮泥沙的空间分布特征及其随时间的变化; 其中第一空间模态显示了渤海表层悬沙的季节变化特征。周舟等(2017)利用2003~2014年的MODIS (Moderate Resolution Imaging Spectroradiometer)卫星遥感数据, 对渤海表层悬浮物浓度的空间分布形态进行反演, 结果显示表层悬浮物浓度在时间尺度上有明显的季节变化, 且渤海表层悬浮物从2003~2014年呈现下降的趋势。近年来, 相对于卫星遥感资料反演悬沙浓度, 利用数值模式对整个渤海海域悬沙分布的研究相对较少(Lu et al, 2011; Liu et al, 2020; Wang et al, 2020)。

利用数值模拟方法不仅可以再现渤海悬沙浓度分布的时空特征, 同时数值模拟也是揭示悬沙浓度分布的动力成因的有效手段。本文基于ROMS海洋模式, 对渤海的水动力环境和悬沙分布进行气候态数值模拟, 重点展示渤海海域的冬夏季流场特征、悬沙分布, 并对悬沙分布及其影响机制进行讨论。

1 研究方法

本文所采用的ROMS模式是基于三维非线性斜压原始方程的开源海洋模型, 模拟范围为36°~41°N, 117°~123°E, 如图 1所示。地形数据主要采用ETOPO1数据, 水平方向上采用正交经纬网格, 网格数为73×61, 水平分辨率为5′, 垂向网格采用σ坐标, 共分为20层, 最小水深设置为1 m。模型计算采用冷启动, 初始水位和流速均设为0, 初始场的温盐数据来自WOA13 (World Ocean Atlas 2013), 在初始场中将沉积物分为3类, 分别为黏土、粉砂和砂, 三种沉积物的占比分布来自Wang等(2020)的研究结果。其大气强迫场数据来自COADS (comprehensive ocean atmosphere data set)资料, 边界条件采用SODA (simple ocean data assimilation)气候态平均场, 潮强迫是在开边界上采用10个分潮(M₂, S₂, N₂, K₂, K₁, O₁, P₁, Q₁, M_f, M_m)的调和常数来确定水位变化, 其数据来源于全球潮波模式TPXO 7.0。波浪强迫是将SWAN (simulating waves nearshore)单独计算的有效波高、波向、波长以及波周期等波浪参数作为外在强迫添加到ROMS中。在边界强迫中, 开边界水位采用Chapman边界条件, 正压流速采用Flather边界条件, 三维斜压流速及温盐场均采用Rad边界条件。在本模型中, 由于在渤海入海的河流中, 辽河、滦河以及海河流量远小于黄河, 因此仅考虑黄河径流及悬浮物的输运作用。对于黄河径流量和悬浮物的数据来自于利津水文站近30 a的气候态数据, 且黄河中的黏土、粉砂和砂的占比分别为20%、70%和10% (Wang et al, 2020)。

本文用气候态数据将模式运行了2 a, 取第2年的12个月份的月平均值来表征每个月的特征。本文采用了2018年7月在渤海进行观测的大面站温盐深(conductivity temperature depth, CTD)数据与模拟的气候态的月平均值进行对比, 同时采用2017年2~3月与2018年8月在秦皇岛海域的实测数据, 对其求10 d的平均值, 与模拟10 d的平均值进行对比。

2 模式结果验证 2.1 潮位对比

对于潮汐的准确模拟是正确模拟渤海水动力环境的前提。取第二年1月份的水位数据, 采用调和分析方法分离出潮位的调和常数。由于渤海的主要半日分潮和全日分潮分别为M₂和K₁, 因此本文仅对M₂和K₁进行分析。图 2c和2d分别为M₂和K₁分潮的同潮图。

通过与Fang (1986)的渤海图集中同潮图进行对比, 可以发现两者之间分潮M₂与K₁的无潮点位置, 振幅以及相位的分布均较为一致。通过对模拟的调和常数与7个渤海验潮站实测的M₂和K₁分潮的调和常数进行定量对比, M₂分潮振幅误差的绝对平均为8.2 cm, 相对误差的绝对平均为11.8%; 迟角误差的绝对平均为14.3°, 相对误差的绝对平均为8.7%。K₁分潮振幅误差的绝对平均为0.7 cm, 相对误差的绝对平均为4.6%; 迟角误差的绝对平均为22.3°, 相对误差的绝对误差为19.4%。以上结果表明, 该模型对渤海潮汐的模拟是准确的。

2.2 流场对比

图 3a和图 3b分别是渤海冬季(1月)和夏季(7月)的气候态垂向平均流场, 从模拟结果来看, 在渤海海峡终年呈现“北进南出”的结构; 冬季沿岸线流动的海流最强, 辽东湾以及渤海中部均呈现顺时针环流结构, 渤海湾和莱州湾海流逆时针流动, 整体海流流速较小, 一般都在0.05 m/s以下, 这与已有的数值模拟结果较为类似(毕聪聪等, 2015; 马伟伟等, 2016; Jiang et al, 2020; Liu et al, 2021)。夏季环流明显强于冬季, 流速很多在0.05~0.1 m/s。夏季, 20 m等深线以深的渤海中部, 环流呈现“双环”结构, 即外部的逆时针环包围着内部的顺时针环, 这与Xia等(2019), Liu等(2021)的研究结果一致。本文的气候态模拟结果可以再现渤海环流的基本结构及其季节变化。

3 结果与讨论 3.1 渤海温跃层的季节变化

为了研究渤海温跃层的季节性变化, 将模式模拟的气候态结果与实测数据进行对比。实测数据分为两部分, 第一部分是利用在S1站点(图 1a黑色空心圆所示)于2017年2月15日~3月9日(冬季), 2019年8月15~29日(夏季)进行的海床基平台观测数据, 并对其温度取10 d平均, 与模拟的气候态温度进行对比, 对比结果如图 4a所示; 第二部分是利用2018年7月对渤海进行的大面站观测数据与模式模拟结果(月平均)进行对比, 如图 5所示, 站点B29~B30, B53~B59 (红色空心圆)位置如图 1a所示。

从图 4a可以看出, 模式对于S1站点季节性变化的模拟准确。图 4b中的S1站点CTD实测数据在夏季的垂向分布显示, 表层温度为27 ℃, 底层温度为24.5 ℃, 没有明显的温跃层存在, 这也与图 4a中的模拟结果相对应。图 4a的模拟温度显示底层温度(Layer 1, Layer 2)与表层温度(Layer 20)相差无几。

图 5显示, 在渤海水深较深的海域, 夏季存在明显的温跃层。由于大面站观测获得的温度是瞬时值, 且这些站点的温度大都是白天所测, 故测得的温度普遍高于模式模拟的温度, 且温跃层以上温度更高, 垂向变化更大。而模式计算的温度是7月的平均值, 温度略低且垂向变化范围较小。此外, 模式模拟的盐度在31~32, 而B29~B30, B53~B59站点所测的盐度大都也在31~32, 进一步验证了模式模拟的准确性, 故可以用模式模拟的结果来分析渤海温度的季节性变化。

由于温跃层的存在使得表层温度与底层温度差距较大, 因此采用表层温度与底层温度的差值来体现温跃层的强度大小与位置分布, 其1~12月份的空间分布如图 6所示。图 6显示, 温跃层在4月开始形成, 7月达到最强, 在9月之后, 开始慢慢消退, 10月已经基本无温跃层的存在, 这与屠金钊(1992)的研究结果较为一致, 且在温跃层强盛期(6~8月)表底层温差的分布与等深线的分布较为一致。

3.2 渤海悬沙分布的季节变化

本文参考Bian等(2013)的泥沙参数方案, 对泥沙各级粒径参数的相关设置见表 1。数值模拟的渤海表层悬沙浓度的季节性分布如图 7所示, 渤海表层的悬沙分布具有明显的季节变化。整体上讲, 渤海的悬沙浓度在冬季(12~2月)最大, 秋季(9~11月)次之, 春季(3~5月)再次之, 夏季(6~8月)最小, 且不论在哪个季节, 悬沙浓度在辽东湾、渤海湾以及莱州湾均相对较高, 大都在70 mg/L以上, 悬沙的时空分布特征与庞重光等(2013)利用遥感资料对渤海表层悬沙的反演结果基本一致。图 7还显示, 秦皇岛海域悬沙浓度显著低于其他海域, 常年偏低, 其表层悬沙浓度基本都维持在20 mg/L以下, 这一现象也与前人对渤海悬沙浓度的反演结果类似(Jiang et al, 2004; Wang et al, 2014; 王震等, 2016; Zhao et al, 2022)。通过对模式模拟的表层悬沙浓度和卫星遥感反演的海表悬沙浓度进行对比, 如图 8所示(位置点A~F如图 1a所示), 发现在渤海的不同位置, 如典型的高、低悬沙浓度点, 模式均能够模拟出表层悬沙浓度的季节变化特征, 且大致与卫星遥感反演的悬沙浓度值一致, 冬季浓度最高, 夏季浓度最低。

3.3 渤海底部切应力的季节变化

泥沙模型分别计算了时间步长为3 min的流致切应力、波致切应力以及波流相互作用产生的最大切应力。通过对剪切应力求季节平均, 得到不同季节剪切应力的大小, 波流相互作用切应力、流致切应力和波致切应力分别如图 9~11所示。从整体来看, 由波流相互作用引起的最大剪切应力在冬季最大、秋季次之、春季和夏季相对较小。流致切应力与波流切应力几乎相当, 波致切应力仅仅在沿岸的浅水区域(20 m等深线以浅)发挥一定的作用, 其中在辽东湾和莱州湾较为显著, 而在深水区, 波致切应力与流致切应力相比, 几乎可以忽略不计, 这与Wang等(2020)的研究较为一致。从波流切应力的季节分布来看, 辽东湾东部的波流切应力在冬季明显高于夏季, 冬季波流切应力基本在0.2 N/m²以上, 夏季波流切应力主要在0.12~0.20 N/m²。但在辽东湾北部, 波流切应力的季节变化不大, 常年维持在0.3 N/m²以上, 这也是辽东湾海域常年悬沙浓度较高的原因。在渤海湾, 波流切应力在夏季甚至要大于冬季, 这主要是因为夏季海流相对冬季更强。渤海湾与辽东湾都是流致切应力占主导, 而莱州湾的流致剪切应力相对较小, 其悬沙浓度则是波流相互作用与黄河径流输入共同作用的结果。此外, 秦皇岛海域其底部的波流切应力常年保持在0.06 N/m²以下, 远远小于其他海域底部的波流切应力。因此, 之所以秦皇岛海域悬沙浓度常年低于其他海域, 最根本的原因在于秦皇岛海域位于M₂无潮点附近, 属于弱潮流区, 底层流速相对较小, 流致切应力相应较小, 底部沉积物发生再悬浮的概率同步较小。

3.4 渤海夏季表层悬沙分布的控制因子

渤海表层悬沙浓度(图 7)和波流底切应力(图 9)的季节分布显示, 在冬季、秋季以及春季, 波流底切应力的空间分布与表层的悬沙分布较为一致。据此推论, 冬季、秋季以及春季渤海表层悬沙大都来源于底部再悬浮泥沙的垂向扩散。然而, 在夏季, 其波流底切应力的空间分布与春季分布相当接近(图 9b和9c), 但其表层悬沙浓度的分布却明显低于春季(图 7b和7c)。本文数值模式的垂向分层为20层, 因此能够清楚地展示季节平均温度、悬沙浓度和流速大小的垂向分布, 分别如图 12~14所示。如图 12a所示, 冬季渤海海域海水温度在垂向几乎没有变化, 与其对应的悬沙浓度在垂向上也基本保持不变(图 13a), 流速大小自海表向下略有衰减(图 14a), 如流速最大的辽东湾北部位置, 流速大小从海表的0.41 m/s, 减小到海水下层第4层的0.35 m/s, 继续衰减到底层第1层的0.30 m/s。垂向基本均匀分布的水温和悬沙浓度表明在渤海, 冬季整个深度都是混合层, 垂向混合强度大、混合较为充分。然而, 在夏季, 辽东湾、渤海湾、莱州湾以及其他近岸海域由于水深较浅, 温度在垂向上依然没有明显变化, 因此浅水区的变化主要还是受动力条件的季节性波动控制, 但在渤海中部, 表层温度与底层温度差别显著, 存在明显的温跃层现象, 如图 12c所示。夏季, 渤海中部的表层温度一般在22~26 ℃, 而底部温度则大都在12~20 ℃。图 12c还显示, 从底层到中层Layer 1~Layer 10, 即温跃层之下, 温度的上升幅度较小, 但相对应的悬沙浓度在垂向上递减的幅度却较大(图 13c)。主要原因是温跃层的存在降低了风、海浪等外力对海底的扰动, 抑制了均温层水体在垂向的混合, 不利于底层再悬浮泥沙的向上扩散。夏季, 虽然水体在垂向上混合强度明显减弱, 但其平流项却更加发育, 使得中上层海水流速增大, 不论是季节平均的流速(图 14)还是冬夏季深度平均流场(图 3), 都是夏季流速大于冬季, 这与孙运佳等(2016)的研究一致。

图 12c展示了夏季温跃层中心所在层位(Layer 13), 该层温度位于上层(Layer 16)和下层(Layer 10)之间, 该层的悬沙浓度相对于下层(Layer 10)也有明显的降低, 除三湾外, 悬沙浓度基本上都低于10 mg/L, 如图 13c所示。由于温跃层具有层结强, 垂向扩散系数低的特征, 温跃层的存在对于底部悬浮颗粒物向表层的扩散起到明显的抑制作用(乔璐璐等, 2010; 雷栋等, 2019), 故温跃层上层的悬沙浓度(Layer 16和Layer 19)显著低于下层的悬沙浓度(Layer 1, Layer 4, Layer 7和Layer 10)。

水体层化现象对流速的垂向梯度具有明显的增强作用, 如图 14所示。在水体下层Layer 1~Layer 10, 夏季垂向流速梯度明显大于冬季, 体现在水体中层(Layer 10), 夏季流速明显大于冬季, 但在水体底层(Layer 1), 夏季流速要小于冬季。这一现象在辽东湾东部最为显著, 夏季中层水体流速约为0.48 m/s, 明显大于冬季中层水体的0.39 m/s, 而夏季的底层流速约为0.28 m/s, 小于冬季对应的底层流速0.32 m/s。水体的层化在一定程度上减弱了底层流速, 从而间接影响了底部沉积物的再悬浮, 这与齐昱恺等(2021) 在三沙湾的观测结果一致。与夏季相比, 春季在一定程度上也有温跃层现象(图 12b), 但强度较低, 这也使得春季的悬沙浓度(图 13b)比秋季(图 13d)更低, 因为秋季水体混合均匀, 几乎不存在温跃层的现象。

4 结论

本文基于ROMS模式对渤海水动力环境与悬沙分布开展数值模拟, 得出如下结论:

在渤海海峡环流终年“北进南出”; 夏季环流明显强于冬季, 并呈现外围逆时针环、内部顺时针环的“双环”结构。在平均深度仅为18 m的渤海, 夏季存在明显的温跃层现象, 温跃层在4月开始形成, 7月达到最强, 9月开始慢慢消退, 温跃层强度分布与等深线较为一致。

渤海悬沙分布具有明显的季节变化, 冬季悬沙浓度最大, 秋季次之, 春季再次之, 夏季最小。控制悬沙浓度的波流底切应力, 秋冬季明显大于春夏季, 且不论哪个季节, 流致切应力在波流剪切应力中占主导地位, 波致切应力仅仅在浅水海域有一定的作用。秦皇岛海域由于位于M₂无潮点附近, 属于弱潮流区, 底层流速相对较小, 相应的流致切应力较小, 故底部沉积物发生再悬浮的程度较小, 这是秦皇岛海域常年悬沙浓度较低的主要原因。

春夏季温跃层的存在对悬沙的垂向交换具有重要影响。温跃层在一定程度上减小了底边界层流速, 增大了流速的垂向梯度; 温跃层的形成对底部悬浮颗粒物向上层的扩散起到明显的抑制作用。夏季, 温跃层的存在是中上层悬沙浓度维持较低的重要原因。但在浅水区, 无论哪个季节, 温度在垂向上都没有明显的变化, 悬沙浓度的变化主要还是受动力条件的季节性波动控制。