CN101865689B - 一种流域咸潮预测方法 - Google Patents

Abstract

本发明涉及海洋灾害监测预警领域，特别是一种流域咸潮预测方法。所述方法的具体步骤为：(11)预测所测量流域的径流，得到径流日平均流量、月平均流量和年平均流量；(12)以步骤(11)所获得的径流日平均流量、月平均流量和年平均流量作为参数代入三维水动力模式，以模拟咸潮动态。通过预先预测流域径流，并以径流流量作为三维水动力模式的输入参数，从而预测出咸潮动态，解决了对咸潮的发生和程度进行预测的技术问题。

Description

一种流域咸潮预测方法

技术领域

本发明涉及海洋灾害监测预警领域，特别是一种流域咸潮预测方法。

背景技术

咸潮是三角洲特别是珠江三角洲最复杂的现象之一，它受到径流、潮汐、地形等诸多因素的驱动、影响和制约，表现出极为复杂的时空变化特征。近十几年来，珠江三角洲冬春季咸潮频繁发生，使当地生产及生活用水受到影响，对国民经济的发展影响极大。广东省位于珠江流域的下游，省内降水的多少直接影响珠江的径流量，而珠江流域径流变化是珠江口咸潮发生的主要原因之一。只有能够精确预测径流量，才能达到精确预测咸潮发生的程度与范围。因此，需要一个能够根据气象预报资料预测珠江径流量的模式，并将预测得到的珠江径流量提供给河网-河口三维一体化数学模式，模拟和预测珠江三角洲冬季咸潮的动态，从而为有关职能部门建立珠江三角洲咸潮预警体系及实施应对措施提供技术支持。

Xue et al曾运用高分辨率无结构网格的有限体积近岸海洋模型FVCOM模拟了长江口的咸潮入侵，并研究了其物理机制(P.Xue，C.Chen，P.Ding，et al，2009：Saltwater intrusion into the Changjiang River：A model-guidedmechanism study.JOURNAL OF GEOPHYSICAL RESEARCH，VOL.114，C02006，doi：10.1029/2008JC004831)。Liu et al曾使用一个三维斜压水动力盐度模型研究了在不同径流条件对台湾Danshuei River estuary余流与盐度入侵的影响(W.-C.Liu，W.-B.Chen，R.T.Cheng，et al，2007：Modeling the influenceof river discharge on salt intrusion and residual circulation in DanshueiRiver estuary，Taiwan.Continental Shelf Research 27，900-921)。这些咸潮研究模型工作原理为提供一定外部驱动条件(外边界、径流、地形等)，根据一系列水动力控制方程模拟出盐水入侵河口的三维动态。陈水森等根据珠江口磨刀门水道长期的含氯度、水位与流量数据推导出磨刀门水道咸潮入侵的经验模型(陈水森，方立刚、李红、李宏丽等，2007：珠江口咸潮入侵分析与经验模型———以磨刀门水道为例。水科学进展18(5)：751-755)，可估算出咸潮入侵的距离。上述模型均以已知径流量为必要条件，可精确描述咸潮入侵的动态或距离，并对其产生机制进行研究，但不能对咸潮的发生和程度进行预测，而降水的多少直接影响了河口咸潮的发生。

发明内容

本发明提供一种流域咸潮预测方法，以解决现有技术不能对咸潮的发生和程度进行预测的技术问题。

本发明采用的技术方案如下：

一种流域咸潮预测方法，所述方法的具体步骤为：

(11)预测所测量流域的径流，得到径流日平均流量、月平均流量和年平均流量；

(12)以步骤(11)所获得的径流日平均流量、月平均流量和年平均流量作为参数代入三维水动力模式，以模拟咸潮动态。

三维水动力模式的径流输入参数并不是已有的数据，而是根据步骤(11)预测得到的径流，因此通过步骤(12)得到的咸潮动态也为经过预测而得到的咸潮动态。

作为一种优选方案，所述步骤(11)所述的具体步骤如下：

(21)通过如下公式计算径流日流量D：

D＝r_min×A+r×A×coe，其中r_min为基本流系数，取流域枯水期平均流量除以其集水面积的结果，A为单元面积，coe为水流至注入点过程中的损耗系数；

(22)汇总一个月的径流日流量得到该月的径流月流量，取其平均值得到径流月平均流量；

(23)汇总一年的径流日流量得到该年度径流年流量，取其平均值得到径流年平均流量。

作为进一步的优选方案，所述基本流系数为0.00043m/d，所述损耗系数满足以下不等式：0.45≤损耗系数≤0.8。

作为更进一步的优选方案，所述步骤(12)采用的三维水动力模式为环境流体动力学模式EFDC。

作为更进一步的优选方案，所述环境流体动力学模式的水平网格采用曲线正交网格坐标，垂向采用sigma网格坐标。

作为再进一步的优选方案，所述sigma网络坐标σ定义为其中η是海表面抬升，D＝h+η。h是相对于大地水准面的瞬时水深，z是在笛卡儿坐标系下的垂向坐标，当σ＝0时表示自由表面，而σ＝-1时表示海底。

作为进一步的优选方案，对于环境流体动力学模式EFDC，通过结合有限体积和有限差分法求解运动方程，变量的交错网格参考C网格或者MAC网格方法，模式计算方案是在交错的网格内采用内-外模分离来求解水平动量方程和连续方程。

上述C网格方法采用Arakawa A，V R Lamb，1977：Computational designof the basic dynamical process of the UCLA general circulation model.Methods in Computational Physics 17 Academic Press，173-265.

上述MAC网格方法采用Peyret，R.and T.D.Taylor：ComputationalMethods for Fluid Flow.Springer-Verlag，Berlin，1983.

通过预先预测流域径流，并以径流流量作为三维水动力模式的参数，从而预测出咸潮动态，解决了对咸潮的发生和程度进行预测的技术问题。

附图说明

图1为珠江流域地形高程(单位：m)，最长的一条红色曲线即为珠江主干流；

图2为西江高要站实测流量与水文数字高程模型H-DEM模拟结果的比较(m3/s)，实线为模拟结果，虚线为实测资料。(a)日平均流量(b)月平均流量(c)年平均流量；

图3为咸潮模式模拟范围及模式网格；

图4为2007年12月虎门在大、小潮期间盐度的入侵(上图为大潮期间，下图为小潮期间)；

图5为2007年12月磨刀门在大、小潮期间盐度的入侵(上图为大潮期间，下图为小潮期间)；

图6为预测的2009/2010冬珠三角咸潮入侵情况。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细的说明。

本实施例以广东省珠江口入海口领域作为测量流域。

第一步，预报珠江口各入海口的径流。国际上基于地形表面的流是沿着最速下降方向的，提出了数字地形高度模式(digital elevati on model，DEM)，该模式可以得到较真实的流域径流，从而使得径流可以以较真实的、物理上可信的方式来驱动海洋环流。本项目使用的水文数字高程模型(H-DEM)，区域选定如图1所示为95.00°E-120.05°E，20.00°N-28.04°N，纬向分辨率为0.05°，经向分辨率为0.03°。区域覆盖了整个珠江流域。以水流离开单元的最大坡度方向为水流方向，采用D8单流向方法予以确定。大气强迫主要有降水和表面气温。单元内的日径流量计算如下

D(x，y，t)＝r_min×A+r×A×coe

(5.1.1)

其中r_min为基本流系数，取西江高要站枯水期平均流量除以其集水面积的结果，为0.00043m/d。A为单元面积。coe为水流至注入点过程中的损耗系数，根据实测结果的校正，coe的上下界分别取0.8和0.45。r为单元内的日水流，为与陆面高程有关的量。

如图2所示，模拟结果得到珠江口入海径流日平均流量，月平均流量从日平均流量计算而来，年平均流量从月平均流量计算而来。

第二步，结合珠江口河网水动力过程，驱动珠江河口三维水动力模式模拟咸潮动态。本技术运用的The Environmental Fluid Dynamics Code(EFDC)模式，是由美国国家环保署资助开发，由美国弗吉尼亚州海洋研究所Hamrick教授等人于1992年开发并适时进行发展和维护，用于模拟湖泊、水库、海湾、湿地和河口等地表水的三维数值模型。从1997年开始被美国国家环保局指定为WASP(水质分析与模拟程序)水质模型的耦合水动力模型部分。模式水平网格采用曲线正交网格坐标，垂向采用sigma网格坐标。动量和连续方程以及盐度和温度的输运方程可以表达为以下的形式：

$h_1{h}_2\frac{\∂\mathrm{\η}}{\∂t}+\frac{\∂}{\∂\mathrm{\ξ}}\left(h_2\mathrm{uD}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(h_1\mathrm{vD}\right)+h_1{h}_2\frac{\∂\mathrm{\ω}}{\∂\mathrm{\σ}}=0---\left(\mathrm{5.2.1}\right)$

$h_1{h}_2\frac{\∂\mathrm{\η}}{\∂t}+\frac{\∂}{\∂\mathrm{\ξ}}\left(h_{2}D{\∫}_{-1}^0\mathrm{ud\σ}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(h_{1}D{\∫}_{-1}^0\mathrm{vd\σ}\right)=0---(5.2.2)$

$\frac{\∂\left(h_1{h}_2\mathrm{uD}\right)}{\∂t}+\frac{\∂}{\∂\mathrm{\ξ}}\left(h_2{u}^{2}D\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(h_1\mathrm{uvD}\right)+h_1{h}_2\frac{\∂\left(\mathrm{\ωu}\right)}{\∂\mathrm{\σ}}+\mathrm{Dv}(-v\frac{{\∂h}_2}{\∂\mathrm{\ξ}}+u\frac{{\∂h}_1}{\∂\mathrm{\ζ}}-h_1{h}_{2}f)$

$=-h_{2}D\frac{\∂}{\∂\mathrm{\ξ}}(\mathrm{g\η}+p)-h_2(\frac{\∂h}{\∂\mathrm{\ξ}}-\mathrm{\σ}\frac{\∂D}{\∂\mathrm{\ξ}})\frac{\∂p}{\∂\mathrm{\σ}}+\frac{\∂}{\∂\mathrm{\σ}}\left(\frac{h_1{h}_2{K}_m}{D}\frac{\∂u}{\∂\mathrm{\σ}}\right)+Q_u---\left(\mathrm{5.2.3}\right)$

$\frac{\∂\left(h_1{h}_2\mathrm{vD}\right)}{\∂t}+\frac{\∂}{\∂\mathrm{\ξ}}\left(h_2\mathrm{uvD}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(h_1{v}^{2}D\right)+h_1{h}_2\frac{\∂\left(\mathrm{\ωv}\right)}{\∂\mathrm{\σ}}+\mathrm{Du}(v\frac{{\∂h}_2}{\∂\mathrm{\ξ}}-u\frac{{\∂h}_1}{\∂\mathrm{\ζ}}+h_1{h}_{2}f)$

$=-h_{1}D\frac{\∂}{\∂\mathrm{\ζ}}(\mathrm{g\η}+p)-h_1(\frac{\∂h}{\∂\mathrm{\ζ}}-\mathrm{\σ}\frac{\∂D}{\∂\mathrm{\ζ}})\frac{\∂p}{\∂\mathrm{\σ}}+\frac{\∂}{\∂\mathrm{\σ}}\left(\frac{h_1{h}_2{K}_m}{D}\frac{\∂v}{\∂\mathrm{\σ}}\right)+Q_v---\left(\mathrm{5.2.4}\right)$

$\frac{\∂p}{\∂\mathrm{\σ}}=-\mathrm{gD}(\mathrm{\ρ}-{\mathrm{\ρ}}_0){{\mathrm{\ρ}}_0}^{-1}---\left(\mathrm{5.2.5}\right)$

$Q_u=\frac{\∂}{\∂\mathrm{\ξ}}\left(2A_m\frac{h_2}{h_1}D\frac{\∂u}{\∂\mathrm{\ξ}}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(A_m\frac{h_1}{h_2}D\frac{\∂u}{\∂\mathrm{\ζ}}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(A_{m}D\frac{\∂v}{\∂\mathrm{\ξ}}\right)---\left(\mathrm{5.2.6}\right)$

$Q_v=\frac{\∂}{\∂\mathrm{\ζ}}\left(2A_m\frac{h_1}{h_2}D\frac{\∂v}{\∂\mathrm{\ζ}}\right)+\frac{\∂}{\∂\mathrm{\ξ}}\left(A_m\frac{h_2}{h_1}D\frac{\∂v}{\∂\mathrm{\ξ}}\right)+\frac{\∂}{\∂\mathrm{\ξ}}\left(A_{m}D\frac{\∂u}{\∂\mathrm{\ζ}}\right)---\left(\mathrm{5.2.7}\right)$

ρ＝ρ(p，S，T)

(5.2.8)

$\frac{\∂}{\∂t}\left(h_1{h}_2\mathrm{DS}\right)+\frac{\∂}{\∂\mathrm{\ξ}}\left(h_2\mathrm{DuS}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(h_1\mathrm{DvS}\right)+\frac{\∂}{\∂\mathrm{\σ}}\left(h_1{h}_2\mathrm{\ωS}\right)=\frac{\∂}{\∂\mathrm{\σ}}\left(\frac{h_1{h}_2{K}_h}{D}\frac{\∂S}{\∂\mathrm{\σ}}\right)+Q_S---\left(\mathrm{5.2.9}\right)$

$\frac{\∂}{\∂t}\left(h_1{h}_2\mathrm{DT}\right)+\frac{\∂}{\∂\mathrm{\ξ}}\left(h_2\mathrm{DuT}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(h_1\mathrm{DvT}\right)+\frac{\∂}{\∂\mathrm{\σ}}\left(h_1{h}_2\mathrm{\ωT}\right)=\frac{\∂}{\∂\mathrm{\σ}}\left(\frac{h_1{h}_2{K}_h}{D}\frac{\∂T}{\∂\mathrm{\σ}}\right)+Q_T---\left(\mathrm{5.2.10}\right)$

$Q_S=\frac{\∂}{\∂\mathrm{\ξ}}\left(A_h\frac{h_2}{h_1}D\frac{\∂S}{\∂\mathrm{\ξ}}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(A_h\frac{h_1}{h_2}D\frac{\∂S}{\∂\mathrm{\ζ}}\right)---\left(\mathrm{5.2.11}\right)$

$Q_T=\frac{\∂}{\∂\mathrm{\ξ}}\left(A_h\frac{h_2}{h_1}D\frac{\∂T}{\∂\mathrm{\ξ}}\right)+\frac{\∂}{\∂\mathrm{\ζ}}\left(A_h\frac{h_1}{h_2}D\frac{\∂T}{\∂\mathrm{\ζ}}\right)+h_1{h}_2\frac{\∂R}{\∂\mathrm{\σ}}---\left(\mathrm{5.2.12}\right)$

在这些方程中，h₁和h₂是方程从笛卡儿坐标转化为曲线正交坐标使时的矩阵转化系数。垂向sigma坐标σ定义为其中η是海表面抬升，D＝h+η。h是相对于大地水准面的瞬时水深。z是在笛卡儿坐标系下的垂向坐标。因此，当σ＝0时表示自由表面，而σ＝-1时表示海底。u和v表示在曲线正交坐标系ξ和ζ下的两个水平速度分量。ω表示垂直于σ面的速度，它可以用笛卡尔坐标系下的垂向速度w来表示：

$\mathrm{\ω}=w-\frac{1}{h_1{h}_2}[h_{2}u(\mathrm{\σ}\frac{\∂D}{\∂\mathrm{\ξ}}+\frac{\∂\mathrm{\η}}{\∂\mathrm{\ξ}})+h_{1}v(\mathrm{\σ}\frac{\∂D}{\∂\mathrm{\ζ}}+\frac{\∂\mathrm{\η}}{\∂\mathrm{\ζ}})]-(\mathrm{\σ}\frac{\∂D}{\∂t}+\frac{\∂\mathrm{\η}}{\∂t})---\left(\mathrm{5.2.13}\right)$

表示科氏力参数。Ω＝7.29×10^-5s^-1表示地球自转的角速度。其中

表示纬度，g是重力加速度，ρ₀是平均密度，压力p是超出参考密度的静力压力ρ₀gD(1+σ)部分与参考密度的比值。Q_S和Q_T分别表示温度和盐度的源汇项。R表示太阳辐射进入水体的大小。K_m和K_h分别表示动量、温度和盐度混合方程中的垂向湍流粘性项和扩散项。它们是通过由Mellor & Yamada(1982)发展的一个2.5层湍流闭合模型来计算的。这个方程被Galperin et al.(1988)所简化。在这个闭合模型中，湍流强度和长度尺度运用分析学上的决定稳态的方法和输运方程的解来计算。水平湍流粘性系数A_m和扩散系数A_h通过Smagorinsky(1963)方程来计算：

其中A_C是经验系数。

上层风应力和底层摩擦应力运用下面的公式来计算：

τ_w＝C_Dρ_sU|U|

(5.2.15)

τ_b＝C_fρ_bU_b|U_b|

(5.2.16)

其中τ_w和τ_b表示风应力和底层摩擦应力。C_D和C_f是表层和底层的摩擦系数。U和U_b分别表示风速和底层流速。

通过结合有限体积和有限差分法可以求解运动方程(5.1.1-5.1.5)，变量的交错网格位置通常参考C网格或者MAC网格方法，模式计算方案是在交错的网格内采用内-外模分离来求解水平动量方程和连续方程。外模态与正压的长波运动有关，它是通过半隐式三次时间标准方案来计算，并伴随着两次时间标准的订正。内模态与水平速度分量的垂向剪切有关，它的垂向剪切项是通过运用分步算法结合垂向剪切项的隐式步骤来求解，其他部分是通过一个显式步骤来计算。

模型网格区域包括珠江所有河网，东江，西江和北江河网区域，以及广海湾(Guanghai Bay)，黄茅海，大亚湾和广东沿岸到60-100米水深区域。经纬度为东经112.6度到115.5度，北纬21.1到23.1度。西江和北江八个口门上边界分别为：虎门最上界到广州，磨刀门最上界到天河水位站(包括磨刀门，鸡啼门和虎跳门)，横门最上界达到小揽，焦门和洪奇门最上界到达沙湾，崖门最上界到达谭江。东江河网最上界达到博罗水位站。整个网格区域水平方向为263×246个网格，垂向分9层。如图3所示，在河网区域网格水平分辨率可达到100米，而在广东近岸海域分辨率为3-5公里。模型的时间步长设为30秒。模型开边界分为东、南和西边界三部分，模型河流上边界在进行咸潮机制研究时由北江石角、西江高要两站日均流量根据多年统计平均的分流比来给出，在进行咸潮发生预测时由第一步数字高程模式预测的径流量作为输入条件。

模拟结果输出网格区域内的三维流场和温盐场，根据咸潮预警盐度0.25psu判断咸潮在珠江河口与河网上溯的程度和距离。

如图4-5所示，在咸潮机制研究方面，采用2007.12-2008.1水文站日径流量数据为输入条件，模拟了该段时间内珠三角的咸潮。模拟结果与同期实测结果相比较吻合较好，水位验证相关系数都在0.93以上，潮位与盐度、流速验证结果如下表(表1，表2)，证实模式可用。通过设计的数值实验，对影响咸潮上溯距离的径流、潮汐、局地风和海平面上升等因素进行了机制性的研究。模式结果和观测结果都表明大潮期间，潮汐作用使得水体垂向混合良好；而在小潮期间随着潮汐作用的减弱，水体表现出较明显的层化；此外虎门水道大潮期间咸潮入侵较强，而磨刀门水道小潮期间咸潮入侵较强。在探讨风对珠江口盐度的影响的敏感性数值实验中，我们发现对虎门和磨刀门而言，在风速大小相同的条件下，西北风能增加盐水线上溯的距离，即盐度入侵的强度达到最大，东北风和西南风对盐度入侵的影响与无风时相比影响很小，而东南风能使盐水线往下推移，也就是东南风能减小咸潮入侵的强度。此外数值结果表明海平面上升会加剧咸潮的危害，比如增加水道的平均盐度以及盐度入侵的距离；冬季大潮时，当海平面上升0.8m时，咸端可能到达广州附近。

表1为珠江口三个潮位站(九洲港、内伶仃、桂山岛)观测(O)与模拟(M)的四个分潮(M2、S2、N2、K1)的比较：

表2为大、小潮时四个实时监测站观测与模拟结果的比较

当模式使用于预测时，基于气象部门2009年华南区域冬春气候趋势会作出的2009/2010年冬春连旱的天气预测，运用本技术对2009年12月至2010年1月咸潮发生的可能状况进行应急预估，如图6所示。

Claims (6)

1.一种流域咸潮预测方法，其特征在于，所述方法的具体步骤为：

(11)预测所测量流域的径流，得到径流日平均流量、月平均流量和年平均流量；

(12)以步骤(11)所获得的径流日平均流量、月平均流量和年平均流量作为参数代入三维水动力模式。

所述步骤(11)所述的具体步骤如下：

(21)通过如下公式计算径流日流量D：

D＝r_min×A+r×A×coe，其中r_min为基本流系数，取流域枯水期平均流量除以其集水面积的结果，A为单元面积，coe为水流至注入点过程中的损耗系数，r为单元内的日水流；

(22)汇总一个月的径流日流量得到该月的径流总量，取其平均值得到径流月平均流量；

(23)汇总一年的径流日流量得到该年度径流总量，取其平均值得到径流年平均流量。

2.根据权利要求1所述的咸潮预测方法，其特征在于，所述基本流系数为0.00041-0.00045m/d，所述损耗系数满足以下不等式：0.45≤损耗系数≤0.8。

3.根据权利要求1所述的咸潮预测方法，其特征在于，所述步骤(12)采用的三维水动力模式为环境流体动力学模式EFDC。

4.根据权利要求3所述的咸潮预测方法，其特征在于，所述环境流体动力学模式的水平网格采用曲线正交网格坐标，垂向采用sigma网格坐标。

5.根据权利要求4所述的咸潮预测方法，其特征在于，所述sigma网络坐标σ定义为

其中η是海表面抬升，D＝h+η。h是相对于大地水准面的瞬时水深，z是在笛卡儿坐标系下的垂向坐标，当σ＝0时表示自由表面，而σ＝-1时表示海底。

6.根据权利要求3所述的咸潮预测方法，其特征在于，对于环境流体动力学模式EFDC，通过结合有限体积和有限差分法求解运动方程，变量的交错网格参考C网格或者MAC网格方法，模式计算方案是在交错的网格内采用内-外模分离来求解水平动量方程和连续方程。

Images