东中国海是宽广的陆架海, 东北通过对马海峡与日本海相连, 东部常年受黑潮水入侵影响, 南部通过台湾海峡与南海相连, 陆地有长江等径流的输入, 另外受季风影响显著, 导致东中国海环流结构很复杂, 有明显的季节变化。对于东中国海环流的数值模拟, 袁耀初等^[1-2]在较早时期主要运用二维正压诊断模式, 在计算机运算水平提高后, 朱建荣等^[3]使用三维正压诊断和三维斜压诊断模式, 苏纪兰等^[4]简述了中国近海的主要环流动力机制, 有学者^[5-7]利用数值模式模拟了东中国海环流及其季节变化, 分析了冬夏两季的环流特征。汤毓祥等^[8]根据实测资料分析了南黄海环流的若干特征, 赵瑞祥等^[9]利用ARGO数据分析了黑潮入侵东海陆架的季节特征, 张松等^[10]利用卫星数据分析了中国近海的海表面温度年际变化。王惠群等^[11]在利用诊断计算和半诊断计算对南海夏季环流结构进行模拟时认为半诊断计算的方法对密度场作出适当的动力调整, 使其与地形、风场等更加匹配。概括而言, 目前对于东中国海的温度、盐度的分布规律已经较为清楚, 也已经有前人借用数值模拟手段研究东中国海的环流结构, 其中对夏、冬两季模拟结果的描述居多, 而对秋季环流模拟特征的描述较少。目前观测船只大多安装了走航ADCP, 走航观测的海流资料包含了太多特定条件下的随机成分, 无法利用这些数据勾画研究海区总的流动状况^[12], 大面三维观测级海流数据现在仍主要依靠大面温、盐观测资料在其覆盖区域的诊断分析获得。2000年生态973航次观测覆盖了东中国海大部分区域, 因此本文基于该航次获得的2000年秋季实测温度、盐度数据, 应用FEOM有限元海洋模式, 进行诊断、强诊断和后报计算, 并对温、盐分布及模拟得到的环流结构进行分析。

1 CTD数据和风场

模拟区域的温度和盐度数据是由2000年10月19日~11月27日现场观测所得, 调查采用大面走航的方式, 得到每个站位(图 1)的剖面温盐数据, 温盐CTD数据垂向分辨率为1 m左右。

将调查得到的温度和盐度插值到模式网格上, 作为模式的初始温盐场。画出表层、25 m层和50 m层水深的温盐分布(图 2a~图 2c), 对温盐分布特征分析如下:东中国海秋季温度分布呈现东南高温并向西北递减的趋势。温度梯度密集区在东南海域和长江口入海口附近, 在东中国海中部附近有一温度锋面波动, 表现为黄海暖流暖舌入侵。在黄海北部表层存在一个低温中心, 黄海东南部海域底层有一个低温中心, 表层19℃等温线闭合, 底层温度低于10℃, 认为是黄海冷水团的残留部分。在长江口附近盐度梯度较大(图 2d~图 2f), 东南大部海域盐度都比较高, 在温度锋面波动海域同时出现盐度锋面的相似波动。25 m层在黄海西南部出现一个低盐中心, 而50 m层在黄海东南部出现一个盐度半封闭中心。北部呈现低温、低盐特征, 在中部存在一个温度和盐度的等值线密集区, 南部温度梯度较大, 盐度值较高, 但未发现明显分布规律。长江冲淡水的低温、低盐特征控制影响着长江口大部海域, 盐度极小值出现在这里。在长江口往南海域出现沿岸低温、低盐的沿岸流, 说明长江径流开始转向向南。在长江口北部, 长江冲淡水、台湾暖流和北部的黄海冷水团发生混合。黄海冷水团温盐表现为双核特征, 西部冷水团和东部冷水团温度都较低, 而东部冷水团平均盐度较高, 底部温度出现最低值。

研究区域的风场数据来源于欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts, ECMWF)发布的第三代再分析资料ERA- Interim, 采用该数据集提供耦合数值模式结果。本文选用的是2000年10和11月两个月平均风场数据, 水平分辨率为1/8°。2000年秋季东中国海南部风速大小为5 m/s左右, 以东北风为主, 北部转为北风, 风力略有减小。

2 三维数学模型 2.1 控制方程和数学方法

有限元法起源于20世纪50年代中期, 已发展成为一种有效的数值计算方法, 有限元法相比于有限差分法优势在于可以针对复杂的海域地形进行更有效的模拟, 网格可进行局部加密处理, 更好地贴合岸线。由德国阿尔弗雷德·魏格纳研究所开发的FEOM^[13]就是一种应用广泛的有限元模式, 它是基于流体不可压缩、静压近似和Boussinesq近似的三维原始方程组开发的, 由风、正压和斜压梯度力驱动的有限元海洋环流模型^[14], 模式控制方程如下:

水平动量方程为

$ \begin{array}{l} \frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} + \mathit{\boldsymbol{v}} \cdot {\nabla _3}\mathit{\boldsymbol{u}}{\rm{ + }}f\mathit{\boldsymbol{k}}{\rm{ \times }}\mathit{\boldsymbol{u}}{\rm{ = }} - \frac{{\rm{1}}}{{{\rho _0}}}\nabla p - g\nabla \eta + \nabla \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{A_{\rm{h}}}\nabla \mathit{\boldsymbol{u}} + \frac{\partial }{{\partial z}}\left( {{A_{\rm{v}}}\frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial z}}} \right) \end{array} $

(1)

连续方程为

$ \frac{{\partial \eta }}{{\partial t}} + \nabla \cdot \int_{z = - H}^{z = \eta } {\mathit{\boldsymbol{u}}\;{\rm{d}}z = P - E} $

(2)

静压方程为

$\frac{{\partial p}}{{\partial z}} = - g\rho $

(3)

温盐扩散方程为

$ \begin{array}{l} \frac{{\partial C}}{{\partial t}} + {\rm{}}\mathit{\boldsymbol{v}} \cdot {\nabla _3}{\rm{C}} - \nabla \cdot {K_{\rm{h}}}\nabla C - \frac{\partial }{{\partial z}}\left( {{K_{\rm{v}}}\frac{{\partial C}}{{\partial z}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{\tau }\left( {C - {C_0}} \right) = 0 \end{array} $

(4)

状态方程为

$\rho = \rho \left( {T, S, p} \right)$

(5)

式中, $\mathit{\boldsymbol{v}} \equiv \left( {\mathit{\boldsymbol{u}}, w} \right) \equiv \left( {u, v, w} \right)$表示球面坐标下的流体速度; $ρ_0$和$ρ$分别表示水体平均密度、密度距平; $p$为静压压强; $η$为水位; $H=h+η$表示水体总深度, $h$为平均海面深度; $f$表示科氏参数; $\mathit{\boldsymbol{k}}$表示垂直向量; $g$表示重力加速度; $∇_3$和$∇$分别为三维散度或二维梯度算子; ${A_{\rm{h}}}$和${A_{\rm{v}}}$表示水平和垂直的粘滞系数; $C $表示温度$T$或盐度$S$, $C_0$表示温盐的初始值; ${K_{\rm{h}}}$和${K_{\rm{v}}}$表示水平和垂直的扩散系数; $P、E$表示降水量、蒸发量; $τ$为松弛时间尺度, 通过调节$τ$的大小得到松弛尺度的强诊断方案, 当$τ→∞$时为后报计算。

海面边界条件为:

${A_{\rm{v}}}\frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial z}} = {\mathit{\boldsymbol{\tau}} _{\rm{w}}}, p = 0$

(6)

${K_{\rm{v}}}\frac{{\partial C}}{{\partial z}} = - q$

(7)

$w = \frac{{\partial \eta }}{{\partial t}} + \mathit{\boldsymbol{u}} \cdot \nabla \eta - \left( {P - E} \right)$

(8)

式中, ${\mathit{\boldsymbol{\tau }}_{\rm{w}}}$为风应力矢量; q为温盐的表面通量。

海底边界条件为:

${A_{\rm{v}}}\frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial z}} + {A_{\rm{h}}}\nabla H \cdot \nabla \mathit{\boldsymbol{u}} = {C_{\rm{d}}}\mathit{\boldsymbol{u}}\left| \mathit{\boldsymbol{u}} \right|$

(9)

$ w = - \nabla H\cdot\mathit{\boldsymbol{u}} $

(10)

式中, ${C_{\rm{d}}}$为底摩擦系数, 取值为0.002 5。

水平粘性系数采用Biharmonic粘滞格式, 消除了Laplacian粘滞格式对于涡旋模拟时损耗比较严重的问题。同时引入新的变量:

$\mathit{\boldsymbol{a}} = \Delta \mathit{\boldsymbol{u}}$

(11)

$\mathit{\boldsymbol{b}} = - \nabla {A_{{\rm{hb}}}}\nabla \mathit{\boldsymbol{a}}$

(12)

由公式(11)首先计算得到a值, 接下来的计算和Laplacian粘滞格式的计算方式完全一致。A_h被A_hb替代, ${A_{{\rm{hb}}}}$为Biharmonic水平粘滞系数, 模式取值等于$2.0 \times {10^{11}}\;{{\rm{m}}^4}/{\rm{s}}$。

垂直粘滞和扩散系数参考基于Richardson数的Pacanowski & Philander^[15]混合层方案给定, 认为在每个体积单元为定常值。具体垂直粘滞系数和扩散系数的计算表达式为^[16]:

$ \nu = \frac{{{\nu _N}}}{{{{\left( {1 + \alpha Ri} \right)}^n}}} + {\nu _{\rm{b}}}, {\rm{}}k = \frac{\nu }{{\left( {1 + \alpha {R_i}} \right)}} + {k_b} $

(13)

式中, $ν_{\rm b}$、$k_{\rm b}$是对应的背景值; $α$和n是可调参数, ${\nu _{\rm{N}}}$表示中性混合值; ${R_i} = \frac{{{N^2}}}{{{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}^2} + {{\left( {\frac{{\partial v}}{{\partial z}}} \right)}^2}}}$是梯度Richardson数, 这里N为Brunt-Väisällä频率, ${N^2} = - \frac{g}{{{\rho _0}}} \times \frac{{\partial \rho }}{{\partial z}}$; 在本研究中, ${\nu _{\rm{b}}} = 1 \times {10^{ - 4}}{{\rm{m}}^2}/{\rm{s}}$, ${\rm{}}{k_{\rm{b}}} = 1 \times {10^{ - 5}}{{\rm{m}}^2}/{\rm{s}}$, n = 2, $α=10$, ${\nu _{\rm{N}}} = 0.1\;{{\rm{m}}^2}/{\rm{s}}$。

2.2 模式设置和计算方案

为了充分利用数据, 减少插值带来的误差, 采用观测站位最外缘连线为模型计算边界, 并全设为开边界, 利用改进逆方法^[1]计算得到的流量分布给定开边界的流量。水平网格使用SMS(Surface-water Modeling System)建立三角形网格, 共计1 360个节点, 2 619个单元, 最小、最大和平均网距分别为7.1、30.5和16.4 km左右(图 3)。垂向坐标采用z坐标, 分29层, 深度分别为0、2、4、6、8、10、15、20、25、30、35、40、45、50、55、60、65、70、75、80、85、90、95、100、110、120、130、140和150 m。数值模式计算所需要的温盐初始值由航次各层观测值应用客观分析法插值到模式网格上, 四周边界上的温盐初始值被保存下来, 作为温盐扩散方程求解时所需的开边界。不考虑热交换及降水和蒸发, 即q=0和P-E=0。计算所用的风场为气候态月平均数据, 由再分析风场数据插值到模式网格点上而得。

根据CFL判据, 并考虑计算效率和模式稳定性, 将时间步长定为90 s。共分三种计算方案: (1)诊断计算模式, 对观测温盐数据连续积分100 d, 过程中温盐场保持不变, 对第100 d的结果进行分析; (2)诊断计算60 d后, 再进行后报计算40 d, 确定强诊断松弛的时间尺度; (3)诊断计算60 d后, 进行40 d的强诊断计算, 对强诊断第40 d的结果进行分析。

为判断模式运行稳定性, 作出整个区块的体积平均动能变化过程曲线(图 4a), 发现体积平均动能在模式运行60 d时已经达到稳定, 可以进行后报计算和强诊断计算。另外, 可以根据密度位相速度$α$的变化过程判断模式调整是否达到稳定, $α$表示如下, 其中“—”代表对整个计算区域做体积平均。

$\alpha = \overline {\left| {\frac{{\partial \rho }}{{\partial t}}} \right|/\sqrt {{{\left( {\frac{{\partial \rho }}{{\partial x}}} \right)}^2} + {{\left( {\frac{{\partial \rho }}{{\partial y}}} \right)}^2}} } $

(14)

诊断计算60 d后开展后报计算, 根据后报计算的体积平均密度位相速度变化曲线(图 4b)可知, 在后报计算开始5 d左右后, 模式的体积平均密度位相曲线和体积平均动能曲线变化都趋于平缓, 由此判定温盐场匹配调整的时间尺度为5 d, 并最终采用松弛时间尺度为5 d的强诊断计算方案。

3 结果与讨论

本文将对诊断计算和强诊断模式下得到的稳定海表面高度和水平环流进行讨论, 并开展对比分析。先对海表面高度分布进行分析讨论, 并进一步根据表层、10、25和50 m层的流场分布, 分析和讨论东中国海秋季的环流的结构。

3.1 诊断计算 3.1.1 海表面高度分布及对比验证

诊断计算结果得到的海表面高度分布如图 5a所示, 观测区域的西北部和东南部海表面高度较高, 在西南部出现一个低值区, 中心位置为(28.5°N, 123.5°E), 浙江近岸海域海平面高度较高, 长江口东北部也出现了海表面高度较大区域, 中心位置为(32°N, 122.5°E)。海表面高度异常会形成压强梯度力, 当压强梯度力和科氏力平衡时, 就形成了地转流。海表面高度一定程度上反映了上层环流结构, 将10 m层流速分布(图 6b)和海表面高度进行对比可知, 东南部的台湾暖流流速较大, 对应海表面高度梯度较大, 同样浙江近岸海域的海表面高度梯度较大。在(30°N, 126°E)出现一个相对平坦的海表面高度分布, 该海域平均流速较小, 同样在其他海域海表面高度和10 m层流速也有较好的对应关系。为了验证诊断计算的海表面高度分布(图 5a), 我们给出了2000年10~11月平均的卫星高度计数据(https://www.aviso.altimetry.fr/en/data/products/sea-surface-height-products/global/ssha.html.)、HYCOM(https://hycom.org/dataserver/gofs-3pt0/analysis)和SODA(http://www.atmos.umd.edu/~ocean/index_files/soda3.3.1_mn_download.htm)各自对应的海面高度结果(图 5b~图 5d)。从卫星高度计海面高度分布(图 5b)来看, 东中国近海潮汐作用强, 而卫星观测的时间分辨率较低(7天1次), 故2个月左右的平均不能完全滤掉潮的影响。一般来说, 50 m以浅潮的影响特别大, 其它区域潮的影响会相对小一些。计算域南部海域水深基本都大于50 m, 卫星高度计海面高度分布在观测区域南部, 南向的浙闽沿岸流和北向的台湾暖流海域两个高值区与诊断结果对应较好, 两个高值区间的低值中心的走向卫星高度计和诊断计算结果基本一致。计算域北部大部分区域水深小于50 m, 受潮的影响比较大, 但大体结构还是有几分相似的, 只是可能在位置的对应上有些出入。HYCOM的水平分辨率为1/12°, 应该算比较高的分辨率了, 但是它计算得到的海面高度分布(图 5c)却差强人意, 原因可能是它没有针对东中国海作特别调校。HYCOM海面高度分布结构单一, 除了计算域东南角与黑潮相关的高值区结构与我们的诊断计算结果一致之外, 没有其它可圈可点的地方。SODA海面高度分布(图 5d)水平分辨率为1/2°。SODA结果中, 计算域南部由于没有出现沿岸流系统, 所以近岸处没有海面高度高值区, 也就没有中间的低值分布区, 所以分布结构显得单一起来。计算域北部SODA的一个海面高度中心与诊断计算结果较为一致, 只是其等值线梯度太大, 甚至比南部强流区都大。

由此可见, 将观测点的外缘连接起来作为计算区域, 划分相连且各自闭合的计算盒子, 开展考虑风效应且满足每个盒子水量和热量平衡的改进逆方法计算, 因此得到模型所需的流边界条件, 并利用动力学完备、且满足水量和动量等守恒的有限元模型FEOM诊断计算得出的结果应该是可靠的。

3.1.2 流速场的分布

表层受到季风的影响较大, 东中国海东南部是广阔的大陆架, 在东北风的影响下, 表层呈现风生流的特征, 北部区域的流速较小, 大小在3.5 cm/s左右, 中部表现为西向流, 大小集中在13~17 cm/s左右, 南部流向为西南向, 在长江口附近出现流速较大区域。中部环流到达长江口附近转向南北近岸, 南部与沿岸流流向一致, 北部为北向沿岸流, 这应该是中部西南风的作用下水体在长江口附近堆积而向南北两岸扩展所致。表层最大流速分布在长江口以南近岸海域, 为浙闽沿岸流主轴, 流向为南向偏西, 在台湾海峡北部有一个气旋式涡旋。

在10 m层流速受季风影响较弱, 环流结构和表层差异大。观测区域中部海域流速较弱, 长江口附近流向和表层一致, 南北向流动依然存在, 流向和表层一致, 流速大小为12 cm/s左右。观测海域南部大部分是北向的台湾暖流, 大小为7 cm/s左右, 在向北流动过程中台湾暖流主轴发生弯曲, 并在29°N附近台湾暖流由北向流转而向东汇入黑潮继续向北, 在31°N附近有分支转向西北向形成黄海暖流。在(34°N, 124°E)有一个气旋式冷涡存在, 同时在10 m层观测温度同样有一个低温中心, 对应较好。在25 m层的环流结构和10 m层的流向大体一致, 中部流向和10 m层变化较大, 出现弱的南向流动。在台湾南部存在一个气旋式涡旋, 涡旋西侧浙闽沿岸流流向为西南向, 大小为20 cm/s左右, 东侧台湾暖流主轴依然发生明显弯曲, 在31°N附近的西北向黄海暖流减弱。

在50 m层浙闽沿岸流仍然明显, 大小为19 cm/s左右, 可以认为, 浙闽沿岸流从表至底存在, 影响可到50 m等深线处, 台湾暖流依旧北向弯曲前进, 最大流速在台湾岛东北部。在观测区域北部, 之前一直存在以上各层的南向的流动减弱了很多, 东北部的冷水团核心比25 m层向东北偏移, 中心在(35°N, 125°E)左右, 范围变大, 周围流动很弱, 核心温度只有7℃左右, 应该是黄海冷水团的残留部分。

3.2 强诊断计算 3.2.1 水位场分布

在进行60 d诊断计算之后进行松弛时间尺度为5 d强诊断计算会对温盐场进行适当的调整, 可使计算结果更加真实合理, 对强诊断结果(图 7)进行对比发现海表面高度场比诊断计算时(图 5a)更加平滑, 但总体趋势未改变。主要变化发生在台湾岛东北部, 强诊断计算的海表面高度略低。观测区域东北部的等高线也趋于平滑, 在长江口北侧海面高度梯度略有减小, 高水位面积略有增大。

3.2.2 流速场分布

由松弛时间尺度为5 d的第40 d强诊断结果速度分布(图 8)可知, 强诊断模拟计算结果和诊断计算结果大体一致, 这里只分析两者的差异。首先, 对比表层的速度场可知, 强诊断结果流速比诊断计算流速较小, 在观测区域北部, 流速大小为5 cm/s左右, 在观测区域中部流速大小为9~11 cm/s左右, 最大流速出现在浙闽沿岸, 大小为34 cm/s左右, 在台湾岛东北部, 最大流速为20 cm/s左右。由10 m层流场分布可知, 在观测区域中部海域, 强诊断计算结果黑潮在10 m层向西发生弯曲, 流速大小为6 cm/s左右, 而诊断计算该海域流速只有3 cm/s左右。浙闽沿岸流主轴最大流速为19 cm/s, 中心位置比诊断计算更靠南。强诊断结果显示, 台湾暖流速度为8~10 cm/s左右, 比诊断计算大3 cm/s左右, 向北运动到31°N附近与北边来的逆流混合后向西运动至长江口。在25 m层强诊断计算结果台湾暖流最大流速为8 cm/s左右, 诊断计算结果流速为9 cm/s, 略有变小。黄海冷水团区域的流速也比诊断计算低一点, 低温水核心面积也减小。在50 m层强诊断计算结果浙闽沿岸流最大流速26 cm/s左右, 台湾暖流最大流速7 cm/s左右, 但是在台湾岛西北部最大流速25 cm/s, 相比诊断计算位置没发生改变, 最大流速也基本不变。综上所述, 强诊断计算结果除了在10 m层流速变大, 其他层强诊断结果流速略有变小, 但是流向基本一致。

4 结论

本文利用在东中国海2000年10月19日~11月27日调查得到的温盐资料, 结合三维非结构网格FEOM模式, 对观测区域2000年秋季三维环流状况进行了诊断和强诊断计算, 得出如下主要结论:

(1) 由观测数据可知, 表层东南部温度比较高, 梯度也较大; 西北部温度较低, 分布相对均匀, 在济州岛西北有一个低温中心, 中心温度在18℃左右。低盐区域在长江口附近, 最大盐度梯度也在长江口附近, 长江口至济州岛连线有东北-西南向的盐度梯度线。观测区域北部25 m层的低盐中心在(33°N, 123°E)附近, 低温中心在(34°N, 124°E)附近, 台湾暖流的高盐水可以到达长江口附近, 但是温度却较低。50 m层在观测区域东部黑潮高温区域和高盐区域基本重合, 而北部黄海冷水团低温和低盐区域差异很大。

(2) 通过诊断计算, FEOM模式较好的再现了观测区域主要的环流特征:秋季受东北季风的影响, 表层流向以西向偏南为主, 表现为风生流的特征; 10 m层及以下流速减小, 流向基本一致, 台湾暖流表现为弯曲向北运动, 后又汇入黑潮, 在10 m层台湾暖流分支对长江口附近水体有影响; 在25 m层, (32°N, 125°E)附近出现一条向西南向的流动, 与南北两侧形成气旋与反气旋涡旋结构; 50 m层浙闽沿岸流依旧强劲, 黄海北部环流结构一直很弱, 气旋结构只在25 m层以上明显。台湾暖流西北部从表层至50 m层流向一致为东北向, 且流速并未明显减小。

(3) 比较诊断计算和强诊断计算结果, 流速、流向在定性上较为一致, 只是大小略有差异:强诊断结果流速在表层流速比诊断计算结果要小, 东南部黑潮流向更偏向北一点; 在10 m层仍有西向流的特征, 表现为仍受表层风的影响, 黑潮有向西的偏转, 这个特征诊断计算出现在25 m层; 强诊断计算温度模拟结果显示, 对于个别区域, 高温、低温的影响面积影响变大, 对于整个区域, 温度等温线的分布更加光滑。