基于HYCOM的西南黄海沿岸上升流分析

潘锡山，李春辉，王彬

（1.江苏省海涂研究中心，江苏 南京 210036；2.南京信息工程大学 海洋科学学院，江苏 南京 210044；3.江苏科技大学 船舶与海洋工程学院，江苏 镇江 210024）

海洋近岸区域最重要的海洋过程之一是上升流，其影响着海洋温盐垂直结构。上升流对于海水营养盐的分布有着重要影响，而营养盐的分布则影响着海洋生物的生存环境，因此上升流对海洋生态环境变化起着重要作用。黄海海域出现的上升流类型主要为沿岸上升流，其变化不仅需要考虑风力因素的影响，还应考虑地形对其调节作用，图1展示了地形影响上升流的一种方式[1]。

图1 沿岸流致上升流示意图，引自吕新刚[1]

通常情况下，上升流由海流侵蚀陆架、等深线辐散、海流局地加速这三种方式诱发。对于地形这一影响因素，有研究发现地形起伏对上升流沿岸分布的影响大于海岸形状的影响。此后有研究人员对于这一方面进行更进一步的探索，发现地形起伏可以引发风生上升流的不均匀分布，在跨陆架方向上存在的地形坡度对上升流的影响更大[1-2]。海岸形状对于上升流的影响更多体现诱导涡度上，部分区域会出现底摩擦诱发横向次级环流的情况增强上升流[1，3]。

近年来，针对黄海上升流，我国的研究人员获得诸多研究成果。苏纪兰和黄大吉发现黄海冷水团的垂向环流结构呈现双环结构[4]；吕新刚借助海浪-潮流-环流耦合模式提出在潮汐影响较强海域的夏季上升流具有斜压“潮生机制”[1]；苏育嵩等认为江苏沿岸低温带是潮混合、风混合和海底坡度共同作用的结果[3]；柴霞等在研究台风“布拉万”过境期间出现异常增温现象时，发现该台风强烈抑制该地区上升流，在风场对上升流产生影响的方面做出了有力的证明[5]。上升流是海洋温盐垂直结构和营养盐分布的重要影响因素，与海洋生物化学过程联系极为紧密，对海洋生态环境具有重大意义[6]。此外，上升流通过物理过程[7]、化学过程两方面影响着全球碳循环。

根据以往研究，我国沿岸海域属于沿岸上升流的多发区[8]。近年来，南黄海西南以及长江口外海域大规模的赤潮、绿潮等灾害频繁爆发，而海水富营养化是其背后的潜在原因[9-13]。径流（长江等）输入的陆源性营养盐逐年增加[14]，这曾被认为是造成该海区富营养化的主要原因。夏季上升流对该海域富营养化起着至关重要的作用，带来的影响不可忽视[15]。然而，迄今为止，该上升流区营养盐垂向输送通量特征尚不明确，上升流携带至上层的高营养盐水对南黄海西南海域的生态环境效应仍需要进一步研究。因此，本文针对西南黄海沿岸上升流的变化特征及成因机制开展相关研究与讨论。

1 数据资料

1.1 HYCOM 模式数据

本研究使用的HYCOM 模式数据为HYCOM 逐日数据，来自网站“http:∕∕apdrc.soest.hawaii.edu∕datadoc∕hycom.php”。HYCOM.org 提供了近实时的全球海洋预报系统输出，本研究的时段使用的数 据为1993 年1 月1 日至2012 年12 月31 日 的 全部数据集，数据集内包含了经度、纬度、水深、海水盐度、海水温度、水平流速等多类数据。实际使用部分的经纬度范围为20°N－40°N、117.04°E－144.96°E，水平分辨率为0.08°×0.08°。数据集内将水深划分为40 层，从海表（0 m） 到海底（5 000 m）不均匀分割，最终使用数据为前15 层（0～50 m）。

1.2 SST观测数据

研究使用的SST 观测数据来自哥白尼海洋服务公司所提供的数据产品——全球海洋OSTIA 海面温度和海冰再处理数据（Global Ocean OSTIA Sea Surface Temperature and Sea Ice Reprocessed），下载网址为：https:∕∕resources.marine.copernicus.eu∕produ ct-detail∕SST_GLO_SST_L4_REP_OBSERVATIONS_010_011∕INFORMATION。本研究实际使用的卫星观测SST 数据，覆盖区域为28.025°N－42.025°N、118.975°E－126.975°E，水平分辨率为0.05°×0.05°，时间范围为1993年1月1日至2012年12月31日。

1.3 ERA5再分析风场数据

ERA5 是第五代ECMWF（欧洲中期天气预报中心）对全球气候的大气再分析数据[16]，属于综合性再分析数据。ERA5把来自模型的数据和全球各地区的观测数据组成一个覆盖全球的连贯数据集，可以提供从1979 年至今的实时数据，数据包括位势高度、温度、经向风、纬向风、垂直风、相对湿度、比湿、散度及涡度等参量的网格数据[17]。通常情况下，ERA5 数据集接收的观测值来自卫星观测与就地观测，分析数据采用数据同化方法[18]。ERA5 数据采用了相当先进的算法，整合了历史观测数据，这有效地提高了其准确性[19]。

本次研究采用ERA5再分析风场数据中的月平均数据，来自网站“https:∕∕cds.climate.copernicus.eu∕#!∕search?text=ERA5&type=dataset”。实 际 使 用部分的经纬度为27°N－41°N、117°E－132°E 区域，其水平分辨率为0.25°×0.25°，时间范围为1993年1月1日至2012年12 月31 日，为距离海表10 m处的水平风场数据。

表1 1998年夏季33.60°N断面垂向上升流流速和海表冷水斑面积的逐月相关系数

2 HYCOM 模式数据检验

2.1 夏季日均区域平均SST的检验

因本研究关注区域位于西南黄海，因此选取覆盖该区域范围（30.025°N－38.025°N、119.975°E－126.525°E）的数据进行检验。从观测的SST 中分别提取该区域范围内6 月、7 月、8 月三个月的日均值，从插值过的HYCOM 模式数据提取HYCOM模式数据在6 月、7 月、8 月三个月的日均值。将新获得的两组数据进行区域平均值的计算，并将结果逐日绘制在同一张图上，得到图2。图2 中蓝色部分为模式数据的SST 日均区域平均值，红色部分为观测数据的日均区域平均值。尽管图2 只绘制了夏季三个月的日均区域平均SST，以年为单位的周期性变化特征仍非常明显，6 月初通常为每年夏季SST区域平均的最低值，8月通常出现每年夏季区域平均的最高值。

图2 研究区域（30.025°N－38.025°N、119.975°E－126.525°E）内夏季平均SST观测数据和模式数据对比

为确定HYCOM 模式数据与观测数据在夏季的区域平均SST 关联度，绘制以模式数据为x 轴、观测数据为y轴的SST日均区域平均值散点图，得到图3 的散点部分。可看出两数据之间应大致为线性相关，于是进行线性最小二乘拟合。计算得到拟合曲线的k 值（拟合曲线公式中自变量的系数）约为0.957 5，可认为模式数据和观测数据之间相似度极高，也证明了模式数据的SST 区域平均值适用于本研究区域。

图3 研究区域（30.025°N－38.025°N、119.975°E－126.525°E）内平均SST的观测数据和模式数据

2.2 夏季西南黄海表层冷水斑变化

上升流可以影响海水的物理化学性质分布。上升流存在的海域，海表会因为上涌的深层或底层水产生能够被卫星遥感检测到的异常理化现象，例如低温、高盐度和高密度、高营养盐含量、高叶绿素、低氧等。本研究中主要关注海洋表面出现的异常低温这一指标，即夏季出现于黄海表层的冷水斑。在夏季西南黄海的江苏沿岸附近，出现的冷水斑中心位置一般位于33°N－34°N、122°E－123.5°E 之间，因此将对冷水斑面积进行统计的研究区域设置在32.575°N－35.075°N、121.125°E－123.525°E内。

图4（a）和 图4（b）分 别为SST 观测 数据 和HYCOM 模式数据夏季1993 年至2012 年多年平均SST 场。可以发现二者呈现出的温度分布均大致为南高北低的分布，而在33°N－34°N、122°E－123.5°E 区间内出现了较为明显的冷水斑；在模式数据多年平均SST 场中， 冷中心温度为22.107 4 ℃，比周围水体低4～5 ℃。

图4 1993－2012年夏季多年平均SST场

尽管图4 展示出的两数据温度分布大致相似，但为了确保冷水斑面积的适用性，需要对冷水斑面积这一要素进行单独的检验。

检验使用的观测数据为30.025°N－38.025°N、119.975°E－126.525°E 提取的夏季日均数据，模式数据使用在SST 检验中插值到观测数据网格上的表层夏季日均数据。

以往对于苏北沿岸夏季表层冷水斑的研究中，通常采用0.03 ℃∕km 的等值线作为冷水斑的边界[20]，但是这一等值线不能保证每年均存在，于是本研究试图采用更为适合的方式划分冷水斑。上升流引起的深层冷水上涌会导致对应区域SST骤降，这会使上升流区边界与周围海水形成较大的温差，即表层冷水斑边界处的温度梯度会比周围区域更大，也就意味着可以用温度梯度最大值对应SST作为冷水斑的SST边界值。考虑到近岸区域的温度梯度最大值对应的SST，存在SST等值线不能闭合的可能，因此采用SST 等值线较为稀疏的东北方向对边界值进行筛选。

分别计算SST 观测数据和模式数据在经向和纬向上的两个日均温度梯度值。对于同一天的同一经纬度，将经向梯度和纬向梯度中较大的一个作为后续处理中所使用的梯度值进行储存。随后考虑到冷水斑中心位置一般位于33°N－34°N、122°E－123.5°E之间，在其东北边界区域（34°N－35°N、122°E－124.5°E）寻找温度梯度最大值对应SST 作为冷水斑的SST 边界值。在32.575°N－35.075°N、121.125°E－123.525°E 范围内，SST 值低于冷水斑边界值的网格点均属于表层冷水斑。

在实际数据处理中，冷水斑东北边界区域（34°N－35°N、122°E－124.5°E）所求得的温度梯度最大值可能存在多个网格点与之对应。因此，先记录所求得的温度梯度最大值，后求东北边界区域内所有温度梯度值与之相同的网格点对应的SST平均值，把这一平均值作为冷水斑的SST边界值。遥感SST数据和HYCOM 模式数据均按照这一方式进行处理。

求得冷水斑的SST 边界值后，将32.575°N－35.075°N、121.125°E－123.525°E 范围内全部网格点上的SST值与冷水斑SST边界值进行对比，记录所有SST值低于冷水斑SST边界值的点的个数和位置，这些点均视为夏季西南黄海表层冷水斑范围内的点。对于多年平均的SST 场，夏季苏北外海海表冷水斑冷中心温度为22.1 ℃，与边界值温度22.4 ℃相差约0.3 ℃。

为确定这种冷水斑SST 边界值确定方法的使用效果，在地图上绘制出记录下的点进行确认。例如图5 所示，这种边界值确定方法对于1993 年至2012 年夏季西南黄海表层冷水斑的使用效果较为优秀，计算得出的冷水斑基本与图4 所示冷水斑范围基本一致。

图5 夏季西南黄海表层冷水斑位置示意图

经过上述一系列操作得到夏季西南黄海表层冷水斑内日均网格点个数，可以求得每年夏季（6月、7 月、8 月三个月）西南黄海表层冷水斑内平均网格点个数。在数据处理中使用的网格为0.05°×0.05°的经纬度网格，于是可将一个网格点换算为30.802 5 km2的面积。将每年夏季西南黄海表层冷水斑内平均网格点个数转换为每年夏季西南黄海表层冷水斑平均面积，将SST 观测数据和HYCOM 模式数据对应的面积变化曲线绘制在同一张图上，得到图6。根据图像信息，对于观测数据而言，夏季西南黄海表层冷水斑年平均面积在1997 年达最大值3.234×104km2，在1998 年达最小值1.519 8 ×104km2；对于HYCOM 模式数据而言，夏季西南黄海表层冷水斑年平均面积在2012 年达最 大 值3.413 3 ×104km2，在1998 年 达 最 小 值2.028 9×104km2。

图6 1993－2012年夏季西南黄海冷水斑年平均面积变化

为确定HYCOM 模式数据与SST观测数据在夏季西南黄海年平均冷水斑面积变化上的关联度，绘制以模式数据为x轴、观测数据为y轴的夏季西南黄海年平均冷水斑面积散点图，随后对离散数据进行线性拟合，得到拟合曲线“y=0.9468x-0.0977”。将拟合曲线绘制在散点图上，得到图7。拟合曲线k=0.946 8，可认为模式数据和观测数据之间相似度极高，也证明了模式数据在夏季西南黄海年平均冷水斑面积这一要素上能够适用于研究区域。

图7 1993－2012年模式数据和观测数据的夏季西南黄海年平均冷水斑面积

3 夏季西南黄海表层冷水斑与上升流

3.1 气候态变化

在对夏季西南黄海表层冷水斑进行检验时，已经获得夏季西南黄海表层冷水斑的分布范围。将表层冷水斑范围内温度最低点作为冷水斑中心，计算得到多年平均的冷水斑中心纬度在33.60°N附近。

在HYCOM 模式数据中，仅提供了每层的水平流速，因此上升流的垂向流速需要使用流体的连续性方程计算得到。连续性方程的一般形式为：

式中：ρ为流体的密度，u、v、w分别为流体在x、y、z方向上的流速。

通常将海水视为不可压缩流体，即海水密度为常数，公式1可化简为：

根据化简后的连续性方程（公式2）进行计算HYCOM模式下1993－2012年的日均垂向流速。

根据上一步得到的日均垂向流速，计算多年夏季6月、7月、8月三个月平均垂向流速，考虑到上升流的流速量级非常小，一般为×10-7～×10-5m∕s，为方便观察，绘图时使用的数值均乘105。沿33.60°N 绘制121°E－124°E 范围内的多年平均断面垂向流速分布图，得到图8，图8 中正值为流速向上方向。从图中可以看出122.5°E附近存在明显的上升流，这一上升的垂向流速随深度的增加逐渐减小，最大值位于海表附近，大于3.6×10-5m∕s；在123.3°E附近存在一个强烈的下降流，这一下降的垂向流速也随深度的增加逐渐减小，最大值同样位于海表附近，大于5.7×10-5m∕s；此外，在123.7°E 附近存在一个较弱的上升流速，这一经度附近最大流速约为3×10-6m∕s。上升流与下降流进行配合，保证了海水的连续性。

图8 冷水斑中心平均纬度多年平均33.60°N断面垂向流速

为确定垂向流速与垂向的温度场匹配情况，绘制33.60°N 这一断面在121°E－124°E 范围内的多年夏季平均垂向温度场，得到图9。可以看出在122.5°E附近为明显的低温区域，海水温度随深度的增加而降低。将图8和图9结合，可认为该断面122.5°E附近存在强烈的冷水上涌现象，上升流主轴位置与图4（b）所展示的冷水斑中心位置相符。

图9 冷水斑中心平均纬度多年平均断面海温

上升流最强烈的分支和最大值均出现在海表冷水斑范围内，这种吻合证明多年变化中上升流对于夏季西南黄海表层冷水斑现象的形成具有直接作用。

3.2 年际变化

为确定年际变化情况下，夏季年均垂向上升流流速与表层冷水斑面积的相关情况，绘制33.60°N 断面处夏季年均垂向上升流流速与表层冷水斑面积的时间序列图像，得到图10。可以发现年际尺度下，1998－2012 年间垂向上升流流速与表层冷水斑面的变化趋势有部分相似。为了确定年平均夏季西南黄海表层冷水斑面积和33.60°N断面处的年平均上升流垂向流速间是否存在相关性，计算年平均夏季西南黄海表层冷水斑面积和33.60°N 断面处的年平均上升流垂向流速的相关系数，为0.445 3，采用P 检验方法对其显著性进行检验，得到对应p 值矩阵非对角元素值为0.045 1，小于显著水平0.05，通过了95%置信度检验。因此可以证明年际变化中，夏季西南黄海的沿岸上升流对于夏季西南黄海表层冷水斑面积有一定程度的影响。

图10 1993－2012年间在夏季年均垂向上升流流速与表层冷水斑面积

为确定年际变化情况下，33.60°N 断面上升流区域的垂向流速与对应区域的表层冷水斑平均温度的相关情况，绘制年均断面上升流区域的垂向流速与对应区域的年均海水温度之间的时间序列，得到图11。可发现年均垂向断面流速与上升流区的年平均海水温度呈现大致相似的波动上升趋势；33.60°N 断面上升流流速最大值出现在2010 年，最小值出现在1993 年；而上升流区SST 的最大值出现在1999年，最小值出现在1993年。因为图像未完全重合，所以需要计算年均断面上升流区域的垂向流速与对应的年均SST 之间的相关系数来确定具体相关情况。计算相关系数后，得到年均断面上升流区域的垂向流速与对应的年均SST 之间的相关系数为0.382 9，对应p 值矩阵非对角元素值为0.045 7，通过95%置信度检验。因此可以证明年际变化中，夏季西南黄海的沿岸上升流对于夏季西南黄海表层冷水斑的SST 有一定程度的影响。

图11 1993－2012年间在33.60°N断面每年夏季年均垂向上升流流速和SST

3.3 季节内变化

根据图6 所示信息，在HYCOM 模式数据中，1993－2012 年间夏季西南黄海表层冷水斑面积最小值年在1998年，为2.028 9万平方千米；夏季西南黄海表层冷水斑面积最大值年在2012 年，为3.413 3 万平方千米。为研究上升流与夏季西南黄海表层冷水斑之间的季节内变化关系，分别对1998 年和2012 年在33.60°N 断面的垂向上升流流速和表层冷水斑面积与平均温度进行相关性计算。

计算两个最值年（最小值和最大值所在年份，分别为1998年和2012年）在33.60°N断面121°E－124°E 范围内的夏季的垂向上升流流速和海表冷水斑面积的逐月相关系数，1998 年的逐月相关系数在表1 中展示，2012 年的逐月相关系数展示在表2 中，表中全部相关系数均通过95%置信度检验。可以看出两最值年月均相关系数主要分布在0.39～0.5 之间；最高可达0.608 0，在2012 年6 月份达到；最低值为0.389 1，在1998 年8 月份达到（全部相关系数均通过95%置信度检验）。

表2 2012年夏季33.60°N断面垂向上升流流速和海表冷水斑面积的逐月相关系数

因此，可以认为33.60°N断面夏季垂向上升流流速和表层冷水斑面积之间的相关性属于中度相关。

计算1998 年在33.60°N 断面121°E－124°E 范围内的6 月、7 月、8 月的垂向上升流流速和冷水斑平均温度的逐月相关系数，可以得到表3。计算2012 年在33.60°N 断面121°E－124°E 范围内的6 月、7 月、8 月的垂向上升流流速和冷水斑平均温度的逐月相关系数，得到表4。能够看出两最值年月均相关系数主要分布在0.49～0.6 之间；最高可达0.710 8，最低值为0.496 5，两个极值均在最大值年达到（全部相关系数均通过置信度检验）。

表3 1998年夏季33.60°N断面垂向上升流流速和冷水斑平均温度的逐月相关系数

表4 2012年夏季33.60°N断面垂向上升流流速和冷水斑平均温度的逐月相关系数

因此，可以认为33.60°N断面夏季垂向上升流流速和表层冷水斑平均海水温度之间的相关性属于中度相关。

4 西南黄海上升流的动力机制研究

已有研究表明，黄海上升流除了风生机制（Ekman），还与潮汐过程有关。由于潮流与风生流之间复杂、强烈的非线性作用，大多研究将两者分开简化讨论。本文基于的HYCOM数据，未引入潮汐过程，因此重点仅讨论上升流的风生机制。

Ekman 输运、Ekman 抽吸可促使上升流产生，而Ekman 输运与Ekman 抽吸均是因为风力因素产生的现象，本质作用是漂流使上层水体产生输运。因此研究中使用了ERA5 再分析风场数据，对多年风应力和风应力旋度分别对表层海水带来的影响进行探讨。

风应力计算的经验公式为：

相应的垂向风应力旋度为：

式中：curl(τ)为垂向风应力旋度，τ→为风应力矢量，U为距海表10 米处风应力矢量，x、y分别为向东和向北方向，Cd为拖曳系数（研究中取Cd=(0.8 + 0.065 |U|) × 10-3），ρa为空气密度（研究中使用值为1.29 kg∕m3）[15]。

4.1 风应力对上升流的影响

本研究中使用的ERA5 再分析风场数据与HYCOM 模式数据使用的经纬度网格不同，为保证后续分析正常进行，将ERA5 再分析风场数据插值至HYCOM 模式数据所用经纬度网格。依据插值后的ERA5 再分析风场数据绘制28°N－36°N、120°E－128°E 范围内的多年夏季平均风场分布，具体如图12 所示。可以看出在夏季西南黄海表层冷水斑范围（32.575°N－35.075°N、121.125°E－123.525°E）内，盛行东南风，将会产生东北方向的Ekman输运。

图12 在28°N－36°N、120°E－128°E范围内多年夏季平均风应力矢量图

根据图12 的陆地分布，可以发现此时风场产生的Ekman 输运是一个离岸的方向，因此需要得到深层或底层水的补充来维持水体稳定。

根据ERA5 再分析风场数据插值结果计算风应力，单位为N∕m2，随后计算关键断面年平均上升流流速大小与风应力大小，将两者的时间序列绘制在同一张图上，得到图13。能够看出二者的变化趋势存在一定的相似性，为确定二者间关系，随后计算关键断面年平均上升流流速大小与风应力大小的相关系数，其值为0.487 2，对应p 值矩阵非对角元素值为0.040 3，小于显著水平0.05，通过95%置信度检验，证明断面平均上升流速大小与风应力大小之间的具有一定的相关性。

图13 1993－2012年间在33.60°N关键断面年平均上升流速大小与风应力的时间序列图像

因此可以认为因风应力而产生的Ekman 输运对于海洋上升流存在一定的影响。

4.2 风应力旋度对上升流的影响

根据垂向风应力旋度计算公式（公式4），计算得到垂向风应力旋度大小，单位为N∕m3，最大值为3.139 3×10-7N∕m3，最小值为-1.335 4×10-7N∕m3。将垂向风应力旋度分布绘制在地图上，可以得到图14。可以发现在28°N－34°N区域及附近的沿岸地区存在为正的风应力旋度，此外的大部分海域风应力旋度为负值。

风应力旋度会引起Ekman 抽吸，理论上可以影响上升流流速，但实际分析过程中对比图14 与图4（b），发现风应力旋度与上升流引起的SST 分布吻合情况并不明显，因此需要计算上升流流速和风应力旋度间的相关系数确认二者之间的关系。

依然采用33.60°N 这一关键断面进行相关分析，分别计算33.60°N断面上每月的上升流平均垂向流速和冷水斑附近的风应力旋度平均值，将两者的时间序列绘制在同一张图上，得到图15。能够看出二者的变化趋势存在一定的相似性，为确定二者间关系，随后计算两组平均值之间的相关系数。可以得出33.60°N断面处，平均上升流流速与风应力旋度的相关系数为0.420 9，对应p 值矩阵非对角元素值为0.037 6，小于显著水平0.05，通过95%置信度检验，表明风应力旋度与平均上升流流速之间具有一定的相关性。

图15 1993－2012年33.60°N关键断面年平均上升流速大小与风应力旋度的时间序列图像

可以认为，这一区域风应力旋度引起的Ekman 抽吸对关键断面平均上升流流速的存在一定的影响。

5 结论

本文应用HYCOM 模式数据、高分辨率卫星遥感SST 观测数据和ERA5 再分析风场数据，对HYCOM 模式数据在西南黄海江苏沿海地区的适用性进行了检验，对于1993年至2012年间夏季西南黄海表层冷水斑进行了统计，研究了长期情况下上升流与夏季西南黄海表层冷水斑间的关系、风应力与风应力旋度与上升流之间的相关性，计算所得相关系数均通过95%置信度检验。得到的主要结论如下：

（1）夏季苏北外海海表冷水斑多年平均冷中心位于33.60°N 断面，冷中心温度为22.1℃，比表层冷水斑范围外的周围水体温度低4～5℃；上升流主轴位置与冷水斑中心位置一致。

（2）夏季苏北外海海表冷水斑面积与上升流流速的相关系数在年际尺度上为0.445 3，季节内二者相关系数在0.39～0.5之间。

（3）夏季苏北外海海表冷水斑平均温度与上升流流速的相关系数在年际尺度上为0.382 9，季节内二者相关系数在0.49～0.6之间。

（4）对于上升流流速，风应力引起的Ekman输运比风应力旋度引起的Ekman 抽吸影响更大。在33.60°N断面处，计算得到风应力与平均上升流流速间相关系数为0.487 2，风应力旋度与平均上升流流速间相关系数为0.420 9。