长江口咸潮入侵预报系统的设计与应用
打开文本图片集
摘要：基于长江口高分辨率FVCOM数值模型，在率定验证可靠的基础上，建立了一个驱动数据下载、模型计算、可视化处理、预报成果发布的集成自动化操作的长江口咸潮入侵预报系统，采用WRF天气模式进行长江口及临近海域的风场预测，从而能够进行未来3 d（72 h）长江口盐度的三维空间分布以及时间变化预报.并通过该系统后报检验了2011年1月份长江口南支堡镇和崇西的盐度变化过程，后报检验良好.操作便捷、自动化运行等优点使得该系统具有一定实用价值，能够方便及时地对长江口咸潮入侵进行预报.
关键词：长江口；咸潮入侵；自动化预报；FVCOM
中图分类号：P751文献标识码：ADOI：10.3969/j.issn.1000-5641.2016.02.015
Abstract：Based on a high-resolution and well-validated FVCOM model， a saltwater intrusion forecasting system in the Changjiang Estuary has been set up.A WRF model is also applied to predict the wind speed and direction in the Changjiang Estuary and its adjacent sea. The forecasting system can forecast three-dimensional spatial distribution of salinity in the Changjiang Estuary at each moment during the following three days （72 hours）.The forecasting system downloads the forcing data， runs the WRF and FVCOM model， visualizes prediction results and publishes the results. All the above operations are automatically processed by computer.The system was validated by hindcasting salinity of two stations， which are Baozhen and Chongxi at South Branch， during January 2011. The comparison of salinity between system hindcasted value and observation show the system is reliable. The system， with advantages of convenient operation and automatic procession， forecasts the saltwater intrusion in the Changjiang Estuary easily and timely.
Key words：Changjiang Estuary；saltwater intrusion；automatic forecast；FVCOM
0引言
上海地处长江入海口，是中国的经济和金融中心，繁荣的国际大都市.淡水资源的取用影响着上海市的经济生产，与上海市民的日常生活息息相关.受长江径流季节性变化的影响，枯水季节长江口的咸潮入侵现象显著，对上海市的淡水资源取用有着严重的影响[1-2]，现今供给上海市用水的三大水库（陈行水库、青草沙水库和东风西沙水库）均位于长江口南支内，因此研究长江口咸潮入侵具有重要的生活和生产意义.以往的研究表明，长江口咸潮入侵受到诸如潮汐、径流、风、地形、水平垂向混合等因子的共同影响[3].近年来，在气候变化、河槽演变等自然因素和诸如三峡工程、青草沙水库等人为因素的作用下，径流、潮汐等对长江口咸潮入侵有着重要作用的因子随之发生了改变.因此，快捷地对长江口盐度的时间变化、空间分布变化进行预报，并且能够及时地向社会实时公布，对相关部门的工作以及市民的日常生活有着重大的作用.
对咸潮入侵的研究主要有实测数据分析和数值模拟两种方法.在大量观测数据的基础上，统计预测模型衍生而出，郑晓琴等通过建立盐度、径流和潮差三者之间的多元回归关系，得到青龙港盐度统计预测模型，较好地预测南支受到咸潮入侵影响的时间和强度[4].但是，基于实测数据的统计预测只能依据定点的数据推测该点盐度的时序变化，不能够很好地刻画其他区域以及整个盐度场的三维空间变化.随着计算机技术的发展，数值模型的预报应用越来越多.朱建荣、吴辉等率先将改进的ECOM-si数值模式应用于长江口数值模拟，较准确地模拟了长江口咸潮入侵，并在此基础上对盐水入侵的动力机制进行了分析[5-8]，还深入探讨了气候变化、重大工程等因素对盐水入侵的影响[9].林唐宇、朱建荣等基于ECOM-si数值模式，利用三维变分数据同化方法，将长江口水文观测资料同化在该模式中，建立了长江口盐水入侵预报系统，较准确地对陈行水库、青草沙水库提供了盐度季度预报[10]；但其预报主要基于季度平均的定常风场，对显著影响咸潮入侵的风场未进行动力学预报，且预报系统未将模型输入文件处理、模型计算、可视化处理等相关过程实现标准自动化，使用有所不便.
本文主要介绍基于FVCOM数值模型的长江口三维盐度场预报系统.建立长江口三维FVCOM咸潮入侵数值模型，在对2013年7月洪季和2014年1月枯季的流场、盐度率定验证较好的基础上，以该模型为核心，集成WRF天气预报模式，建立了长江口咸潮入侵预报系统.本系统实时更新径流数据，自动下载天气模式驱动数据，实现了输入数据的制作、模式计算、输出结果可视化、平台发布等一系列操作的自动化.最后通过野外台站OBS盐度数据进行后报检验，证明系统的可信性.
1数值模型的设置和率定检验
1.1FVCOM模型简介
FVCOM是无结构三角形网格架构、有限体积、自由表面、三维原始方程海洋数值模型，该模型水平上通过无结构三角形网格进行空间离散，对复杂的岸线地形有很好的拟合能力，垂向上采用σ坐标系或通用垂向s坐标系拟合底部不规则地形.数值计算上采用有限体积积分方式对控制方程进行离散求解，该方式能够较好地保证数值计算中质量、动量、温度、盐度等物理量的守恒性[11-13]，由于该模型在数值处理方法和岸线地形拟合上的优势使其在河口海岸地带有着广泛的应用.
1.2模型计算区域网格和设置
计算区域如图1所示，向西沿着长江河道延伸至大通，向东延伸至外海，覆盖了整个长江口、杭州湾、舟山群岛等区域.总计有10 6117个三角形元素，56 180个三角形节点，开边界处分辨率最低，约为14.7 km.如图2所示，徐六泾至口门区域具有较高的空间分辨率，最高可达150 m左右，能够很好地刻画该区域的岸线水深特征，该区域水深数据来源于2013年海图数字化资料，崇明岛的岸线来自2013年的遥感数据.东海陆架及大洋水深地形资料来源为美国国家海洋大气管理局（NOAA）国家地理数据中心（NGDC）的ETOPO1 Global Relief Model数据.
垂向上采用随地形变化的sigma坐标，垂向分20层.采用GOTM湍流闭合模型，湍流模型计算方程的选取、边界条件类型的设定以及湍流参数设定均参考GOTM使用说明书自带咸淡水混合的案例（Liverpool Bay）[14]，较模型原提供的Mellor-Yamada 2.5阶湍流闭合模型能更好地模拟近岸浅水区域的垂向混合.底摩擦系数根据离岸远近分区给定不同数值，从浅水向深水逐渐减小[15].
模型后报检验所使用的历史风场数据由欧洲中期天气预报中心（简称ECMWF，http：//apps.ecmwf.int/）提供，该数据按月份提供下载，空间分辨率达到0.125°×0.125°，时间分辨率为3—6 h，能够满足模型的计算需求.
外海开边界通过给出13个主要天文分潮（M2、S2、K1、O1、N2、K2、P1、Q1、M4、MS4、MF、MM和MN4）的调和常数计算的水位进行驱动，水位计算公式为ζ=ζ0+∑13i=1fiHicos（σit+（Vi+ui）-gi）.（1）式中，ζ为潮位，ζ0为余水位（此处设为零），fi为交点因子，Hi为振幅，σi为角频率，ui是交点订正角，Vi+ui则表示初相角，gi是迟角.调和常数，即Hi和gi的数值来源于葛建忠等建立的东中国海大区域的模型计算的结果，经过实测数据的率定检验，比较可靠[15].上游大通站的径流资料来源于长江水文网（http：//.cn/pages/sssq.html）.模式采用冷启动方式，流速和水位的初始场设为零，温盐初始场来自《渤海、黄海、东海海洋图集》数字化资料[16].
1.3模型的率定检验
华东师范大学河口海岸学国家重点实验室于2013年7月和2014年1月在长江口进行了两次大规模的定点水文调查，分别代表了长江口的洪季和枯季的水文特征，观测站位分布如图3所示，洪季和枯季观测位置一致，覆盖了长江口内主要汊道、口门及长江口外等区域.
由于南支区域是咸潮入侵的主要影响区域，本文主要给出南支至口外的主要站点的验证图.图4、图5是洪季验证图，图6、图7是枯季验证图，从左到右表示各站点的表层、中层、底层的验证结果，黑线表示模型计算值，红点表示实测值.
模型对于流速流向的验证精度比较高，能够较好地模拟实际的流场变化.流速验证统计参数见表1，CC表示相关性系数，RMSE表示均方根误差.各站点的表中底的流速验证CC值基本上都能达到0.75，说明计算值和实测值具有很强的相关性，均方根误差也控制在了较小的范围内.从图6结果来看，K1404站点计算值整体偏大，导致CC值比较低，这可能与观测时候的自然条件变化有关，但是计算结果对涨落潮时刻的模拟符合实际情况.
表2是各站点盐度的验证结果的统计.表中所示位于南支内的03、04、05、08站点WS值接近于0，一方面由于这些站点观测时段本身盐度值接近0且数据量比较少，不便于计算WS值；另一方面由于模型本身在计算该区域时候，由于盐度值比较低，比较难刻画出盐度值的微小变化.剩余的站点WS值基本控制在了0.5以上，通过对这些数据的率定验证，证明了模型对长江口流场、盐度场的模拟具有一定可信性，为预报系统的建立打下了基础.
1.4WRF天气模式简介和设置
WRF（Weather Research and Forecasting Model） 模式是国际上广泛应用的下一代中尺度数值天气预报系统，可用于大气动力学研究和业务化预报需求.其包含两个动力学计算核心、数据同化系统、可扩展的并行软件架构.模型在空间尺度上可以覆盖数十米至数千千米的动力学过程.WRF模式可基于观测和分析数据进行大气动力学模拟，可在不同平台上进行扩展，方便业务化预报系统的使用，提供了先进的物理过程、数学计算和数据同化模块.
为更好地模拟长江口天气风速过程，WRF模型设置为三重嵌套的计算方式，大模型覆盖西北太平洋区域，分辨率81 km；中模型覆盖东中国海及陆架区域，空间分辨率为27 km；小模型聚焦于长江口及临近海域，分辨率较高，为9 km（见图8）.垂向上分为31层，3个区域采用双向数据交换的嵌套方式.
1.5风场模拟验证
采用2006年5月大戢山、滩浒、佘山与嵊山气象站10 m高度风速资料对东海及长江口WRF天气模式进行后报检验，风速后报模拟与观测如图9所示，WRF模式较好地模拟了长江口区域的风速变化过程，可用于长江口的短期（3～7 d）的风速风向预报.
2自动化预报系统介绍及应用
2.1系统架构及流程
系统基于Linux平台进行开发，通过实时下载大通径流数据作为模型河流边界输入参数，并实时下载最新GFS（Global Forecast System）数据作为WRF模式预报的大气驱动.通过基于高速InfiniBand网络支持的MVAPICH2进行并行计算以提高预报时效，输出结果的处理是通过多种可视化软件（NCL，ViSit，Gri）实现的，最后将可视化结果上传至服务器.将上述操作整合到一个shell脚本下，在Linux系统下提交该脚本，便可实现整个流程的自动化.该脚本对整个预报系统作业过程集中控制，包括对预报系统的启动、中断、结束以及执行中间过程的灵活控制，作业序列任何断点和作业子序列的重启，单个作业的重新启动执行，控制文件配置，对预报系统运行情况进行监控和信息获取等功能.
整个系统的运行流程以及各部分的运行时效如图10所示，系统的运行主要包含数据下载处理、数值模式计算、可视化输出及发布四个部分.
数据下载处理主要包括实时径流数据、风场数据、温盐初始场数据的下载处理.
（1） 大通径流数据从长江水文网下载，预报日之前的径流数据使用大通的历史资料，预报日之后的数据用预报前数日的径流数据外插做预测，下载后自动将其转化成适合模型读取的NetCDF格式，每日自动下载的数据自动作为后期预报的历史径流数据.
（2） NCEP GFS风场数据（ftp：//ftpprd.ncep.noaa.gov/pub/data/nccf/com/gfs/prod）每日18时以后更新，因此系统设定每日18时开始自动下载.
（3） 温盐初始场数据是由《渤海、黄海、东海海洋图集》数字化而来，根据起算月份的不同选取温盐资料，插值到计算区域网格上，从而得到温盐初始场输入文件.
数据下载完成后，将NCEP资料作为WRF模式计算的初始条件，计算未来3 d的风场数据，然后转化为适合FVCOM模型读取的NetCDF格式.长江口咸潮入侵FVCOM模型同时调用GOTM湍流模型对整个网格区域进行数值计算，预报天数为3 d.由于盐度的空间分布受径流及大气驱动影响显著，因此只有在基于可靠的径流及风场预测的基础上对盐度的预报才具有实际应用意义.对径流而言，由于采用大通径流，一般大通径流需要5～7 d内影响到长江口，而目前对径流的预测相对没有太大定期预报，因此在采用当前大通径流量的情况下，河口盐度分布的预报时间段约为5～7 d.在本文的预报模型中默认采用3 d 72 h的预报，主要是由于在该时间段内风场的预测精度较高，从而能更好地进行盐度预报.本文采用的WRF模型使用GFS数据进行预测，可以最长选择7.5 d（180 h）的预报.
上述过程通过MVAPICH2在分布式内存架构集群上使用4（节点）×8（CPU）进行并行计算.从天气模式开始运行至咸潮入侵计算完成总共需要耗时3 h.
2.2系统结果的可视化输出
数据可视化部分使用NCL，ViSit，Gri进行自动化作图，输出主要包括三个部分：主要站点的盐度水位数据的时间序列、盐度的平面分布、主要断面的盐度分布，站点和断面位置如图11所示.
2.2.1主要站点的盐度水位时序变化
该系统输出永隆、崇西、南门、堡镇这4个站点以及陈行、东风西沙、青草沙这3个水库取水处的盐度、水位的时序图，受限于篇幅，本文展示4个站点的时序图.如图12所示，表示各站点从2015年1月16日（即预报当日）0时起未来4 d的水位和表底层盐度的时间序列，红线表示表层，黑线表示底层.
2.2.2盐度平面分布
盐度的平面分布主要包括整个长江口的表底层盐度分布以及长江口南支部分的表底层盐度场分布.图13表示预报系统展示的标题时刻长江口盐度的表、底层平面分布情况，系统能够提供今天以及未来3 d各时刻的盐度分布，并整合成gif动态图形式，显示盐度平面分布随时间的动态变化.计算区域网格在南支分辨率比较高，而且上海市三大水库均位于该区域，因此，预报系统细化了长江口南支的盐度分布，如图14所示，整合的动态图能够较好地展示北支倒灌盐水团进入南支后的动态变化.
2.2.3主要断面盐度分布
如图15所示，预报系统提供三条断面的盐度分布：断面1是从南支开始，沿着深槽，一直至北港；断面2也是从南支开始，一直到南槽，一路沿着深槽变化；断面3基本上是沿着深水航道变化的.
预报系统提供这3个断面今天以及未来3 d各时刻的盐度分布图，同时也会展示其动态变化.断面1和断面2的上端均位于南北支分汊口附近， 而断面1的下端处于北港内， 断面2的下端处于南槽中，因此断面1和断面2的盐度分布图能够刻画出随潮汐变化的北支倒灌和外海咸潮入侵的强度变化情况.断面3位于北槽深水航道的位置，能够较好地刻画北槽咸淡水垂向混合的情况.
2.3系统的应用
目前该系统已经在国家海洋环境预报中心安装并处于试运行阶段，系统主要输出长江口区域72 h内咸潮入侵的空间分布和垂向结构，方便业务化预报的相关可视化界面正在制作中.
3实测数据的后报检验
河口海岸学国家重点实验室从2010年开始在堡镇、南门、崇西、永隆等站布放OBS进行长期观测，数据通过GPRS/GSM网络实时回传，之后又增加了青龙港、长兴岛等站点，实时观测这些站点的盐度数据.本文使用2011年1月份的数据进行后报检验，取南北支分汊口附近的崇西站点，以及南北港分汊口附近的堡镇站点做对比.由于预报系统属于后报模拟，没有历史热启动文件，系统自动识别为冷启动，并提前16 d进行计算，采用大通径流历史资料，将WRF模式计算的风场作为风应力驱动，取表层盐度结果和实测值对比，后报检验结果如图16所示.崇西站受北支倒灌的影响比较大，而北支倒灌又与涨落潮有着密切的关系，因此盐度的日变化比较大，该月份实测数据部分时间段有缺失，但是模型能够较好地刻画数据正常的时间段内随着涨落潮交替的盐度变化.堡镇站离南北支分汊口较远，相比崇西站，外海咸潮入侵的影响比较大，整体盐度数值在小潮到大潮的衔接段时间内有较大的增长，这段时间内盐度随着潮汐日变化也比较显著，模型和实测数据均反映了这些规律.从这两个典型站点的后报检验效果来看，该系统能够较好地反映长江口南支以下至口外的盐度变化规律.
4结论
本文详细介绍了一个以长江口FVCOM模型为核心的咸潮入侵模式，在流场、盐度率定验证可靠的基础上，采用经过验证的WRF天气模式进行长江口及临近海域的风场预测，从而能够进行短期（3 d）长江口盐度的三维空间分布以及时间变化预报.该咸潮入侵预报系统集输入文件处理、数值模拟计算、数据可视化一体，输出包括站点的盐度水位时序变化、长江口以及细化南支的表底层盐度分布、断面盐度分布等，显示了长江口盐度的三维空间和时间变化，并且已经在国家海洋环境预报中心试运行.通过2011年1月份后报检验，证明了该系统在受咸潮入侵影响比较大的长江口南支预报效果比较好，具有一定的可信性.综合来讲，该系统为长江口咸潮入侵的预报提供了一个良好的平台.
[参考文献]
[1]顾玉亮， 吴守培， 乐勤.北支盐水入侵对长江口水源地影响研究[J].人民长江，2003，34（4）：1-3.
[2]贺松林， 丁平兴， 孔亚珍.长江口南支河段枯季盐度变异与北支咸水倒灌[J].自然科展，2006，16（5）：585-587
[3]沈焕庭， 茅志昌， 朱建荣.长江河口盐水入侵[M].北京：海洋出版社，2003.
[4]郑晓琴，肖文军，于芸，等.基于径流和潮汐的长江口盐水入侵统计预测研究[J].海洋预报，2014，04：18-23.
[5]吴辉.长江河口盐水入侵研究——北支倒灌、深水航道工程和冬季季风的影响[D].上海：华东师范大学河口海岸科学研究院，2006.
[6]李路.长江河口盐水入侵时空变化特征和机理[D].上海：华东师范大学河口海岸科学研究院，2011.
[7]WU H，ZHU J R，CHOI B H. Links between saltwater intrusion and subtidal circulation in the Changjiang Estuary： A model-guided study [J]. Continental Shelf Research，2010，30：1891-1905.
[8]LI L，ZHU J R，WU H. Impacts of wind stress on saltwater intrusion in the Yangtze Estuary [J].Science China Earth Sciences，2012，55（7）： 1178-1192.
[9]裘诚.长江河口盐水入侵对气候变化和重大工程的响应[D].上海：华东师范大学河口海岸科学研究院，2014.
[10]林唐宇.横沙东滩淡水资源研究及长江盐水入侵预报模型[D].上海：华东师范大学河口海岸科学研究院，2014.
[11]CHEN C，LIU H，BEARDSLEY R C.An unstructured grid，finite-volume，three-dimensional，primitive equations ocean model：Application to coastal ocean and estuaries[J].Journal of Atmospheric and Oceanic Technology，2003，20（1）：159-186.
[12]CHEN C，COWLES G，BEARDSLEY R C.An unstructured grid，finite-volumecoastal ocean model：FVCOM user manual[R].2nd ed. SMAST/UMASSD Technical Report-06-0602，2006：45.
[13]CHEN C，BEARDSLEY R C，COWLES G.An unstructured grid，finite-volumecoastal ocean model （FVCOM） system[J].Oceanography：Special Issue，2006，19（1）：78-89.
[14]UMLAUF L， BURCHARD H， BOLDING K. GOTM. Sourcecode and test case documentation.Version 3.2[R].MeereswissBer， 2005：63.
[15]葛建忠.Multi-Scale FVCOM model system for the East China Sea and Changjiang estuary and its applications[D].上海：华东师范大学河口海岸科学研究院，2010.
[16]海洋图集编委会.渤海黄海东海海洋图集[M].北京：海洋出版社，1993.
[17]WILLMOTT C J. On the validation of models[J].Phys Geogr，1981（2）：184-194.
（责任编辑： 李万会）