基于河口三维斜压数学模型的珠江河口咸淡水比例研究

黄广灵，邱 静，谭 超，黄本胜

(1.广东省水利水电科学研究院，广州 510635；2.广东省流域水环境治理与水生态修复重点实验室，广州 510635；3.河口水利技术国家地方联合工程实验室，广州 510635)

1 概述

世界很多发达城市均位于河口地区[1]，河口区复杂的水动力交汇过程给河口水资源利用带来与内陆河流不一样的问题，如咸潮对淡水资源的影响等，因此，合理利用及管理河口淡水资源是解决工业及城市居民生活用水的关键问题。珠江河口三角洲网河是粤港澳大湾区重要水源区[2]，然珠江三角洲受潮汐作用影响，下游河口区淡水资源受到近岸高盐陆架水入侵的“污染”[1]，致使湾区局部地区淡水资源的利用受到较大限制，特别是对水质盐度要求的生活、生产取水，每年枯季均会有数天因取水口盐度超标而无法实现取水，严重影响大湾区特别是下游区的取水保障率。因此，从河口区可利用水资源量的角度，需要对珠三角感潮区的咸淡水(咸淡水指河口区含氯度大于250 mg/l的微咸水)进行定量计算及分析。

对于珠江河口咸淡水，以往研究主要关注咸潮的入侵规律[3-6]、影响因素[7]、预报技术[8]以及抑咸措施[9-11]咸潮活动自身的特征，对咸潮活动影响下的河口水资源利用研究目前仍较少，且也仅限于定性分析[12]，缺乏定量的计算。随着最严格水资源管理制度的不断深化，咸淡水作为国家鼓励的非常规水资源，如何计入用水总量控制指标，均直接关系到水资源“三条红线”控制及其考核管理，尤其是取水许可管理、取用水量指标统计、水资源费征收等诸多水资源管理问题。计算咸淡水的比例成为河口区咸淡水管理的核心问题，目前的咸淡水比例计算常规方法是根据水文站长序列盐度监测资料进行统计，但存在受水文站点布设制约，不能全面反映河口区不能水域的咸淡水比例，存在较大的局限性。因此，本文将采用数学模型计算，研究珠江河口年尺度下的咸潮变化情况，以此定量分析河口咸淡水资源量，以为河口区咸淡水资源的管理提供科学依据。

2 三维斜压数学模型

本文采用SELFE[13]模式构建了珠江网河—河口湾—近海陆架整体三维斜压数学模型。

2.1 控制方程

模型基于静压近似及Boussinesq近似求解三维浅水方程和温盐输送方程。 模型求解的主要未知变量包括：自由水面高程、流体速度、温度和盐度矢量。在笛卡尔坐标系统下，模型的主要控制方程如下：

(1)

(2)

(3)

(4)

(5)

(6)

式中：

∇——哈密顿算子；

(x,y)——水平笛卡尔坐标，m；

z——垂向坐标，向上为正，m；

t——时间，s；

η(x,y,t)——自由水面高程，m；

h(x,y)——水深，m；

w——垂向流速，m/s；

f——柯氏力因子，S-1；

g——重力加速度，m/s2；

ψ(φ,λ)——潮汐势，m；

α——有效地球弹性因子；

ρ(x,t)——水的密度，默认参考值ρ0为1 025 kg/m3；

pA(x,y,t)——自由水面大气压强，N/m2；

S,T——水的温度和盐度，psu；

υ——垂向涡动粘性系数，m2/s；

μ——水平涡动粘性系数，m2/s；

κ——温度和盐度的垂向涡动扩散系数，m2/s；

Fs,Fh——输移方程中的水平扩散系数；

Q——太阳辐射吸收率，W/m2；

Cp——水的比热，J/(kg/K)。

模型中所有方程采用半隐格式离散。连续方程及动量方程同时求解，通过底边界条件，将连续方程(1)及动量方程(3)耦合。动量方程(3)中的平流项采用欧拉-拉格朗日(ELM)方法求解，而物质输移方程(4)(5)则可采用ELM或者迎风有限体积法求解。

2.2 模型范围及网格

建立范围覆盖了珠江三角洲-河口-外海一体的咸潮计算数学模型。模型上边界选取西江高要，北江石角，东江博罗，流溪河及潭江上游，外海下边界最外取100 m等深线处，模型网格如图1所示。模型采用全三角形的非结构性网格，以精确的贴合复杂的河岸边界，模型网格有10万多个节点，17万多个单元，网格单元的尺度跨度为0.01～1.2 km。模型在垂向上采用Sigma坐标，共分20层，均分水深，分层厚度随水深而变化。模型中珠江三角洲网河所用地形采用2005—2010年的大范围组合地形，河口湾及近岸海区地形采用2000—2015年海图数字化成果，外海地区地形采用 ETOPO1全球海洋地形数据。

图1 珠江河口全三维模型网格示意

2.3 模型范围及网格参数设置

模型采用斜压模式，其中外海盐度边界取海水盐度34 psu，上游开边界盐度均设为0。水平网格内部各节点的初始盐度采用以下方法设定：海域节点采用WOA09盐度资料进行插值，网河节点则直接设为0。模型的干湿判定水深选取0.01 m。湍流闭合模型中动量方程和输移方程中的水平及垂向涡扩散系数给常数10-6。在底边界摩擦力的处理上，本文选用二次的阻力公式，网格中各节点的阻力系数分河段设定，并通过模型数次率定结果进行调整。

模型计算的基准面为珠基，计算时间步长为 180 s，网格内各节点初始水位均设为0 m，模型计算50 h后水位稳定，计算约8 d后盐度场达到稳定。模型计算水动力边界条件包括外海及上游，其中外海是潮汐调和常数，调和常数由中国海域潮汐预报系统Chinatide计算得到，共选取了珠江口海区9个主要分潮：Q1，O1，P1，K1，N2，M2，S2，K2、SA。模型上游边界均采用流量边界。

2.4 模型率定和验证

模型主要选取了2001年2月的珠江三角洲河网区同步测量水文资料进行潮位、流速、盐度的率定，率定参数包括底床糙率和扩散系数，糙率的设置对潮位、流速起主要影响作用，而扩散系数则主要影响盐度模拟结果，扩散系数越大则盐水入侵范围越快越大。模型提前于率定数据起算，总计算时长为30 d，选择后8 d(2001年2月7日14:00—2001年2月14日20:00)进行水位、流速及盐度率定。为了使模型更好地模拟珠江三角洲网河的潮汐涨落过程，选择了八大口门的主要测站及网河区的各潮位站进行潮位率定。率定站点共计33个，覆盖整个珠江三角洲。各主要测站的计算潮位与实测潮位的对比如图2所示。率定结果显示，绝大部分测站的潮位误差均在0.1 m以内，模型能够较好的反映真实的原型特征。

图2 潮位率定结果示意(2001年2月)

模型共选择网河区及口门处共17个测站的流速资料进行率定，主要率定结果如图3所示，从流速的率定和验证结果来看，大部分测站的计算流速与实测流速吻合得较好，仅有少数几个站点的流速误差稍大，但总的误差不超过10%，考虑到珠江网河及地形的复杂性，模型的流速计算结果可以较真实地反映水流的实际情况。

图3 流速率定结果示意(2001年2月)

图4为模型率定计算(计算时间为2001年2月)中各口门主要测站的盐度计算值与实测值的对比；从图4中可以看出：盐度的模拟精度相对于潮位及流速有所降低，但总的来说计算盐度的误差在合理范围以内，各站的盐度平均误差均不超过实测范围的20%[14]，能够较好地体现出盐度的大小潮变化特征。考虑到地形、风等因素对盐度的变化有很大影响，而模型地形与率定验证时段地形有所差异，同时也没有把风场加以考虑，故计算结果存在充许范围内的误差符合实际情况。

图4 盐度率定结果示意(2001年2月)

选择了2005年1月的实测潮位及流速对模型进行验证，另选择了2005年1月、2008年1月、2009年12月和2011年1月共4组盐度实测资料对模型的盐度计算进行验证。在4组验证计算中，验证结果均显示模型参数选择是合理的，各计算要素的误差均在允许范围之内(见图5)。综上可知，利用本次建立的模型进行长时间序列咸潮入侵计算合理可行。

图5 盐度验证结果示意(2005年1月)

2.5 计算边界选择

珠江河口咸淡水比例变化受上游径流影响较大，受枯季流量的影响尤为显著，因此，上游流量过程采用1960—2009年每年10月—次年3月的平均流量进行排频，计算枯季平均流量保证率，本文珠江河口咸淡水比例计算的主要典型年选取西江下游控制站高要站+北江下游控制站石角站50%、75%、90%枯季来水保证率、东江下游控制站博罗50%、75%、90%枯季来水保证率。

采用西江高要水文站、北江石角水文站及东江博罗水文站1960—2009年共50 a的流量资料，分别将西江下游控制站高要站+北江下游控制站石角站历年枯季平均流量(以下简称“高+石”)、博罗站历年枯季平均流量以P-Ⅲ适线法进行频率计算，得到高+石及博罗站不同枯水保证率对应流量及典型年，在此基础上选择年内洪枯季流量差异最大的典型年作为计算边界，最终选取的典型年见表1所示，其余上边界流溪河老鸦岗及潭江石咀流量典型年的选取与高+石一致。

表1 不同枯水保证率对应流量及典型年 m3/s

3 结果分析

珠江河口水资源管理对象主要为各取水户，其取水规模相对较稳定，因此，其取水中咸淡水的占比可以转化成河道咸淡水时间占比计算。珠江河口河道内咸淡水时间比例(简称咸淡水比例)数学模拟计算采用方法如下：本次计算的上边界水文条件共采用3组，分别是上游枯季来水保证率为50%、75%和90%的年份的全年逐日实测流量。计算时长为1 a，从典型年的前1 a的10月到该典型年份的9月，并提前30 d起算，以保证盐度场达到稳定状态。计算结果的处理：计算结果可得到每一个点的全年逐时含氯度过程，从而可统计每个点全年咸淡水的比例，按全年含氯度大于250 mg/L的小时数除以全年总小时数得到结果百分率统计，如式(7)所示。各点咸淡水比例经空间线性插值后在绘出其咸水比例的等值线。

(7)

式中:

Sper——计算点的咸淡水比例；

Ts——该点计算年中含氯度大于250 mg/L的小时数；

Ty——计算年的总小时数。

3.1 计算咸淡水比例与实测数据的对比

3.1.1磨刀门水道

表2为枯季平均流量保证率为50%、75%和90%的水文条件下磨刀门水道挂定角、大涌口、广昌泵站、灯笼山站及平岗泵站的年咸淡水比例计算结果。与根据实测资料的统计结果对比可知，在同样的枯季平均流量保证率下通过模型计算得到的咸淡水比例更高。如根据马口水文站枯水平均流量枯水保证率的统计结果，2005年高+石的流量保证率为74%，2005年的挂定角、大涌口水闸和平岗泵站的实测年咸淡水比例分别为51%、41%和13%，而高+石枯季平均流量保证率为75%水文条件模型计算得到挂定角、大涌口水闸和平岗泵站3站咸淡水比例分别为62.0%、48.0%和9.0%。

表2 磨刀门水道监测站点咸淡水比例计算结果与实测结果 %

模型计算结果与实测值相比存在一定的误差，主要有下几点原因：① 河口咸淡水混合受诸多因素影响，而模型模拟过程未能考虑所有影响因子，只考虑了其中最重要的，故计算结果存在一定误差是合理的；② 与数据的处理方法有关，本文中所出现的咸淡水比例结果代表的是站点所在的断面水体的垂向平均盐度超标小时数除以1 a小时数(365 d×24 h)的算术结果，实测值则是某个取水口或者水闸所在位置的某一水深处的检测值的统计结果，而在磨刀门水道的同一个断面的不同位置(如主槽和边滩、支汊)，或者在同一个位置的不同水深处，其盐度值都有着较大的变化，因此，模型计算结果和实测统计值之间会存在一定范围内的误差；③ 上游流量边界不完全一致。

图6为磨刀门水道咸淡水比例实测值及计算值的沿程变化示意。由磨刀门水道咸淡水比例沿程变化可知：计算结果和实测资料统计结果均表现为从下游向上游咸淡水比例逐渐降低；而从不同水文条件下的咸淡水比例变化可知：计算结果和实测资料统计结果均表明咸淡水比例变化与上游来流量有一定的负相关关系，枯水平均流量保证率越高则相应的年咸淡水比例也越大。

图6 磨刀门水道咸淡水比例实测值、计算值沿程变化示意

3.1.2东四口门

表3为东四口门黄埔站、中大站、三沙口站、南沙站和冯马庙站的年淡水比例计算结果，图7为珠江狮子洋水道咸水比例实测值与计算值的沿程变化。根据与实测资料的统计结果对比，三沙口站、南沙站和冯马庙站的咸水比例计算与实测值之间的差异较小，如三沙口站的咸水比例模型计算结果为26%～37%，实测统计结果为33%～38%；南沙站的模型计算结果为33%～45%，实测统计结果为35%～40%；冯马庙站淡水比例计算结果为12%～18%，实测统计结果为25%。黄埔站和中大站的咸水比例模型计算结果要小于统计结果，而且误差较大，与实测资料的局限和统计方法有关。该站实测资料的咸淡水比例计算方法是采用每日高低潮特征时刻盐度值计算，1 a中1 d高低潮2次盐度监测资料，超标1次为0.5 d、超标2次为 1 d，咸水比例=含氯度大于250 mg/L /365 d。由此可知：三沙口站、南沙站和冯马庙站的咸水比例计算结果与实测资料统计结果相符，而黄埔站和中大站的咸水比例模型计算结果与实测资料统计值存在较大误差是合理的。因为三沙口站、南沙站和冯马庙站位于河口咸水区的中段位置，在枯水期的大部分时间内都不在咸水界移动范围内，因此，按超标1次为0.5 d、超标2次为1 d是比较符合实际的；而黄埔站及中大站位于河口咸水区的上段位置，枯水期的大部分时间内都在咸水界移动范围内，假如其1 d盐度超标1次，极有可能该天咸水界刚好在此位置，而导致实际其也只超标1～2 h，却按0.5 d(12 h)计，最终咸水比例的统计结果将会比实际大很多。由此可知：模型计算得到的年咸水比例比较符合实际。珠江咸水比例从沿程变化可知：计算结果和实测资料统计结果均表现为从下游向上游咸水比例逐渐降低。

表3 珠江监测站点咸淡水比例计算结果与实测结果 %

图7 珠江咸淡水比例实测值、计算值沿程变化示意

3.2 典型年水文条件下珠江河口咸淡水比例的空间分布特征

图8为计算得到的3个典型水文年条件下的珠江河口区全年咸淡水比例等值线分布示意。从图8中可知：不同保证率典型年水文条件下的咸淡水比例等值线的分布规律基本一致。咸水比例0%等值线系计算典型年中咸水上溯的最远距离，其位置主要受年最枯流量的影响，3组典型年水文条件下，咸水比例0%等值线的位置变化不大。从咸水比例等值线的疏密可知：0%～10%等值线较稀疏，且咸淡水比例越小各等值线之间的距离就越大，主要是咸水最大上溯距离在枯水期对流量的变化最为敏感(在较枯流量范围内，流量的较小变化能引起咸水界的较大距离的变化)；同时也因最枯流量等级流量出现的频率很小。伶仃洋的咸水比例等值线分布与其盐度场的分布一样具有“东高西低”的特征，咸水比例100%的等值线在淇澳岛与内伶仃岛之间。磨刀门河口咸水比例100%的等值线位于交杯滩外2～3 km处，形状由河口内向外突出。黄茅海的咸水比例100%等值线向西南上游凸起，咸水比例呈“西高东低”之势，主要是由于径流经东部深槽下泄，减小了东部的咸水比例。

对比不同典型水文年下的咸淡水比例情况，可知不同典型年水文条件下相同的咸淡水比例等值线的位置变化不会太大，且一般情况下均是枯季平均流量保证率越高，各咸淡水比例等值线越往上游。不同频率水文年下的咸淡水比例等值线的位置变化不足以说明珠江网河区的咸淡水分布具有相对稳定性，在非极端水文条件下，网河区各处咸淡水比例均在一定的幅度内变化。从图8中亦可知：网河区偶有局部地区在枯季平均流量保证率高时候咸淡水比例等值线却更往下游，主要是因为典型年的年内流量分布并不一致，同时也从另一方面反映了珠江网河区的径潮动力的复杂性。图9为不同频率水文年高要+石角水文站年内日均流量分布，从图9中知：枯季平均流量保证率90%的水文年的年内日均流量并不都小于枯季平均流量保证率75%的水文年的年内日均流量，如从0到100 d，保证率90%的日均流量便大于保证率为75%的日均流量，也从而导致网河区咸淡水比例等值线位置分布与对应枯季平均流量保证率不完全一致。

图8 枯季平均流量不同保证率水文条件下的咸淡水比例示意

图9 不同频率水文年高+石日均流量年内分布示意

3.3 河口区主要取水口咸淡水比例

统计珠江河口区23家主要取水户的不同保证率下的咸水比例可知，基本上所有取水户的咸水比例均是在90%保证率的水文条件下最大，咸水比例受流量大小影响比较明显，越是枯水年，各取水口的咸水比例越大。从河道沿程的取水口的年咸水比例变化来看，越往下游其咸水比例越大，而且越往下游，咸水比例增大的速度越快。其中，枯水期平均流量为50%、75%和90%保证率的水文条件下，取水口年咸水比例超过30%的取水户分别为8户、10户和11户。

将2019年珠江河口区主要电厂实际取水量及咸淡水比例折算应计入用水总量控制部分水量，珠江河口主要电厂2019年实际取水量为319 807万m3，按50%、75%和90%来水保证率计，取用咸淡水量分别为76 829万m3、97 920万m3和109 816万m3。

4 结语

1) 建立基于非结构网格的珠江河口三维盐水数学模型，模型主要选取了2001年2月的珠江三角洲河网区同步测量水文资料进行潮位、流速、盐度的率定。选择2005年1月、2008年1月、2009年12月和2011年1月共4组盐度实测资料对模型的盐度计算进行验证。水动力及盐度验证结果良好，能够较好地反映原体实际情况。

2) 利用珠江河口三维盐水数学模型，选取上游50%、75%、90%枯季来水保证率典型年，下游选取河口典型潮位过程，计算河口区不同水域(取水口)全年盐度超标小时数并统计作图，实现不同水文条件下珠江河口各水域盐度值及不同水域盐度超标时间的快速简易查询。提取珠江河口主要取水口咸淡水比例的计算结果与实测数据进行对比，发现计算结果与实测值基本相符。

3) 根据珠江河口主要取水户咸淡水比例结果，在较典型水文条件下，取水口年咸水比例超过30%的取水户超过8户；主要电厂实际取水量中，咸淡水比例约为24.0%～34.3%。

4) 河口咸淡水区水资源利用具有独特的区域特点，因此，建议利用本文提出的咸淡水比例计算方法，对河口大型取水户的取水中的淡水和咸淡水进行量化，并将其量化结果应用于河口区用水总量控制、水资源管理日常统计、水资源公报统计、水资源征收等方面，从而推进珠江河口区水资源精细化管理，为最严格水资源管理制度提供科学支撑。