海平面上升对珠江河口东四口门盐水入侵的影响

刘忠辉，宏波

(华南理工大学，广东广州510641)

盐水入侵是由于河口邻近海域的高盐水团随潮汐涨潮流沿着潮汐通道向上游推进，咸淡水混合使河口上游区域水体变咸的现象，会对沿途区域的水环境、农业灌溉和生活供水等产生重要影响。珠江三角洲城市群拥有城镇人口5 500万，已成为全世界最大的城市片区[1]。该区域约占中国国土面积0.5%的土地，但其生产总值却占中国总量的近20%[2]，是中国最重要的经济区域之一。在全球海平面上升的背景下，珠江三角洲地区已成为全球受盐水入侵严重影响的地区之一[3]。从2000年开始，珠江三角洲沿海城市几乎每年都会受到由于盐水入侵导致的淡水供应不足问题[4]。2005年12月，珠海、中山和广州等地均受盐水入侵影响，取水口水闸氯度严重超标，最高达7 500 mg/L[5]。

众多研究表明，海平面上升会强化盐水入侵。Rice等[6]利用HEM-3D模型对美国东海岸切萨皮克湾的盐水入侵长度进行模拟预测，研究发现海平面上升1 m将使得10 ppt等盐线向上游推进10 km，这将改变现有的水质梯度，威胁当地饮用水供应。Hong等[7]利用EFDC模型研究美国东海岸切萨皮克湾的河口盐度与物质输运过程对海平面上升的响应，模拟结果显示，海平面上升将使水体分层加强，同时引起盐度上升，盐水入侵长度增加，并且影响程度与淡水流量有明显关系。孙志林等[8]利用FVCOM模型对钱塘江河口进行模拟，以2012年为基准年，模拟结果表明海平面上升1 m的情况下，大潮期0.45 ppt等盐线向上游推进了6 km，小潮期推进了9.8 km。陈维等[9]预测海平面上升1 m时，1月长江河口的南支受北支倒灌过来的咸水的影响将明显增大，南支水域平均盐度均超过0.45 ppt。

海平面上升对珠江河口盐水入侵的影响也非常显著。何用等[10]利用磨刀门咸潮整体物理模型试验，对海平面上升引起的咸潮运动变化特点进行分析，发现海平面上升使磨刀门河道沿程盐度均增大；以2012年的海平面为基准，每升高10 cm最大咸界的上移距离大致为6.8 km。关帅等[11]利用FVCOM模型预测得海平面平均每上升1 cm，磨刀门水道盐水入侵长度将增加0.85 km。

前人对珠江河口磨刀门盐水入侵的研究报告较多[10]，而对东四口门的研究相对较少，且东四口门的盐水入侵问题将影响到广州、东莞、佛山、中山等重要城市超过3 000万居民的生产生活。鉴于此，本文采用三维水动力模型EFDC(Environmental Fluid Dynamic Code)模拟珠江河口区域的水动力过程，以2007年为基准年，通过设定多种海平面上升幅度的模型计算模式，对珠江河口的盐水入侵情况进行模拟预测，分析珠江河口东四口门(虎门、蕉门、洪奇沥、横门)盐水入侵受海平面上升的响应，为相关部门应对海平面上升背景下该区域淡水资源管理、利用提供有效参考。

1 研究方法

1.1 研究区域

珠江河口是位于中国南海北部的一个大型河口，拥有范围广阔且复杂的河网，共有8个入海口门，分别是：虎门、蕉门、洪奇沥、横门、磨刀门、鸡啼门、虎跳门和崖门，形成“三江汇流，八口入海”的水系特征[11]，而东四口门是指八大口门中向东注入伶仃洋的4个口门[12]。珠江河口位于亚热带季风气候区，夏季高温多雨是洪季，冬季温和少雨是枯季。其承接西江、北江、东江三江的径流量，多年平均总径流量为281.1亿m3[13]。

为了得到真实的水动力过程，本文将模拟区域由珠江河口扩展至其邻近海域，模拟网格见图1。

a)模型网格

b)珠江河口图1 珠江河口及邻近海域的模拟网格

珠江河口的东四口门区域见图2，该区域河网复杂，有多条支流，其中R1—R12是东四口门的12条主要支流，图上还标注了东四口门的4个自来水厂的位置，即新塘水厂、石溪水厂、沙湾水厂和东涌水厂。

图2 东四口门及其主要河网

1.2 研究方法

本文采用的三维水动力模型EFDC是由美国威廉玛丽大学维吉尼亚海洋科学研究所John Hamrick等人开发的三维地表水水质数学模型[14]，可实现河流、湖泊、河口和海洋等水体的水动力和水质模拟，是一个多参数有限差分模型，已经在很多河口海岸地区有过成功应用[15-17]。

模型利用曲线正交网格模拟珠江河口及其邻近海域，垂向上有10个σ层。模型以珠江河口2007年的自然条件为基准，包括地形、风场、径流量等。潮汐谐波数据由美国俄勒冈州立大学潮汐预测软件(OTPS)算得，并用于模型开边界的驱动。开边界的温度和盐度数据取自世界海洋数据库，风力数据来源于香港机场观测站(图1b)采集的2007年风场数据。河流径流量数据采用珠江三条主要支流的上游各站点2007年的径流量观测数据，即石角站(北江)、高要站(西江)和博罗站(东江)。各个淡水源的流量分配按照珠江流域水资源保护局公布的比例进行，见表1。

近几十年来，全球海平面一直呈现上升趋势，并且海平面上升的速度不断增加。根据美国Climate Change Science Program(CCSP) 2009年工作报告指

表1 珠江河口各口门径流量分配比例

出，到2100年海平面可能升高1.0～1.1 m[18]。Intergovernmental Panel on Climate Change(IPCC) 预测到21世纪末海平面将上升0.82 m[19]。从规避风险的角度来说，对未来海平面上升估计的结果应该更偏重于预测上升值的99分位[20]，即区间内至少有99%的概率包含真值。基于典型浓度路径RCP8.5，到2100年全球平均海平面上升的99分位是0.4～1.8 m[21]。故本文以2007年为基准，设定5种模式来研究珠江河口全年对不同海平面上升幅度的响应，模式命名分别为Base、P30、P50、P100、P150，分别代表基准情况，海平面上升30、50、100、150 cm的情况。各模式之间除了海平面高度不同以外，其余设置均相同，各模式从河道上游输入的径流量都按照2007年珠江河口的实际径流量进行。

1.3 模型验证

将珠江河口2007年的模型结果与实际观测数据进行比较，验证模型的可靠性。

1.3.1水位验证

取珠江河口区域的8个站点(三灶、香港、MO1、九洲港、内伶仃、桂山岛、灯笼山、黄埔，图1b)作为水位验证点，所使用的水位验证实测数据均为2007年在各站点的实际观测数据。本文利用Wilmott[22]提出的Skill值来反映模型模拟数据与实测数据之间的匹配程度。Skill按下式计算：

若模型模拟结果与实际观测吻合越好，Skill值越接近于1，反之趋近于0。本文所选的8个验证点的Skill值见表2。结果表明模型的Skill值与1相近，模型模拟结果与实际观测吻合程度较好。

表2 所选验证点的水位验证数据

各站点2007年12月的水位实测值和模拟的水位值随时间的变化过程曲线见图3。各个站点的水位实测值和模拟值都非常接近，模型模拟的水位与实际水位能较好吻合。

a)三灶 b)香港

c)MO1 d)九洲港

e)桂山岛 f)灯笼山

g)黄埔 h)内伶仃

图3 各个站点水位验证曲线

1.3.2潮流验证

以图1b中的MO1站点作为潮流验证点，把验证点处的2007年7月第1周的模拟数据与2007年该站点实测潮流数据进行对比。验证点水深8.3、9.3 m处的潮流矢量变化情况见图4。图中可以看出，8.3、9.3 m处的实测潮流与模拟潮流的流速与流向都吻合良好。

a) 水深为8.3 m

b)水深为9.3 m

图4 MO1站点的潮流验证

1.3.3盐度验证

以图1b中的MO23站点作为潮流验证点，验证该处在珠江典型枯水期模型模拟的盐度变化情况与实际观测结果是否相符。具体验证的时间为2007年第359.5～361.5天，验证点表层、中层和底层盐度变化见图5，验证结果反映模型模拟的盐度在表层、中层和底层都与实际观测的盐度非常吻合，模型模拟效果良好，满足使用要求。

a) 表层图5 MO23站点的盐度验证

b) 中层

c) 底层

续图5MO23站点的盐度验证

2 结果

2.1 盐度变化

珠江河口各个口门的上游就是珠江三角洲城市群，若出现盐水入侵将对周边区域居民的生产生活造成严重影响。模型模拟的结果包含整个模拟区域全年的4D盐度数据。通过对盐度进行筛选，对各口门1年内盐度超过1 ppt的天数进行统计，得到了各个口门在不同海平面上升幅度情况下，盐度超过1 ppt的总天数，见图6。海平面上升对珠江河口各个口门的盐度有明显的影响，海平面上升幅度越大，盐度超过1 ppt的天数越多。海平面上升1 m时，各个口门的高盐天数相比基准情况增多的天数统计见表3。

图6 天数统计

口门增加天数/d虎门63蕉门103洪奇沥60横门39磨刀门110鸡啼门84虎跳门155崖门153

图6和表3中的结果表明，在相同海平面上升幅度的影响下，不同口门的响应情况是不同的，其中虎跳门和崖门对海平面上升最敏感，海平面上升1 m时高盐天数相比基准情况分别增加155、153 d。通过研究各口门盐度超过1ppt情况的出现时间，发现在海平面上升的情况下，各口门高盐现象仍是多发生于春、冬两季，也是珠江的枯水期，这与基准情况相似，但高盐现象的持续时间会随海平面上升幅度的增加而延长。

当自来水厂取水口处水体的盐度超过自来水取用标准(0.5 ppt)时，自来水厂将可能采取停止自来水供应的措施。珠江上游区域4个自来水厂，水厂位置见图2，全年最高盐度和超标天数的统计见表4。将海平面上升1 m的模拟结果与基准情况进行对比可以看出，各个水厂盐度受海平面上升影响非常明显，盐度超标现象开始出现，尤其是东涌水厂的最高盐度将增加7.84 ppt，同时盐度超标天数增加了114 d。

表4 天数统计

2.2 盐水入侵长度变化

珠江河口各口门上游区域河流分叉较多，单以某一条支流的盐水入侵长度来代表一个口门的盐水入侵并不充分，故将各个口门的主要支流(图2)都作为研究对象进行分析。模拟的珠江河口在不同海平面上升幅度情况下，2007年1月各主要支流平均的0.5 ppt等盐线位置见图7。海平面上升幅度越大，0.5 ppt等盐线向上游推进越多，盐水会入侵到越上游区域，加剧盐水入侵。

图7 平均0.5 ppt 等盐线所在位置

海平面上升对珠江河口的盐水入侵有明显的强化作用。基准情况下，洪季各口门都不出现盐水入侵，但当海平面上升150 cm时，即使在洪季，各口门也会出现不同程度的盐水入侵问题。把图2中标定的各口门位置作为计算盐水入侵长度的起算点，以各支流0.5 ppt等盐线沿河道至对应口门的距离作为盐水入侵长度。海平面上升100、150 cm时，东四口门各支流盐水入侵长度相比基准情况向上游的最大推进距离见表5。由于蕉门(R7、R8)、洪奇沥(R9、R10)和横门(R11、R12)各处盐水都未入侵到各自的河流分叉点(图2中A、B、C为分叉点)，故这3个口门各自的支流的入侵情况相同。表 5的结果可以看出海平面上升1 m时，各支流的0.5 ppt等盐线相较于基准情况都向上游推进了较大距离，其中R2的推进距离最大，达21.371 km；蕉门最大推进9.644 km；洪奇沥最大推进9.753 km，横门最大推进4.819 km。

表5 0.5 ppt等盐线相比基准情况最大向上游推进的距离 km

3 讨论

受Yang等[23]和Chen等[24]将海平面上升幅度与盐度变化通过数据回归的方式得到拟合度较好的直线或曲线的启发，并考虑到本文模型设定的几个计算模式之间，仅有海平面上升幅度一个变量，便将各条支流盐水入侵长度进行月平均，与海平面上升幅度进行线性回归，得到海平面上升幅度与各支流盐水入侵长度的拟合直线。选择入侵情况最严重的1月份作为线性回归分析对象，得到的回归直线(图8)可以用于预测未来海平面上升某一特定幅度时，珠江河口各口门上游区域盐水入侵的距离。每条回归直线的误差范数见表6。误差范数都很小，则回归直线可靠，说明珠江河口各支流由于海平面上升幅度变化引起的盐水入侵长度变化与海平面上升幅度之间存在着较好的线性关系。得到的线性回归关系可用于预测特定海平面上升幅度下某支流的盐水入侵长度，指导工程实践，如自来水厂迁址、选址等。

表6 R1-R12线性回归直线误差范数

a)虎门R1 b)虎门R2 c)虎门R3

d)虎门R4 e)虎门R5 f)虎门R6

g)蕉门R7-R8 h)洪奇沥R9-R10 i)横门R11-R12

图8 海平面上升幅度与盐水入侵长度的回归直线

4 结论

海平面上升会加剧珠江河口东四口门的盐水入侵问题。在海平面上升1 m的情况下，东四口门高盐天数(盐度≥1 ppt)都有明显增加，其中蕉门增加最多，达3个多月。而建在东四口门上游的自来水厂也会出现由于海平面上升引起的取水口处盐度超标的问题，其中新塘水厂和东涌水厂受到的影响最大，海平面上升1 m情况下，全年盐度超标天数将增加3个多月，所以在未来各个水厂必须要考虑迁址。

模型模拟结果还表明海平面上升会使东四口门的盐水入侵长度增加，其中虎门各支流盐水入侵长度增加最多，海平面上升150 cm情况下，0.5 ppt等盐线最多将向上游推进27 km。通过对珠江河口东四口门主要支流在一月份的平均盐水入侵长度与海平面上升幅度进行线性回归，得到了可靠的回归直线。借助回归直线，可以预测未来海平面上升特定幅度时，珠江河口各支流的盐水入侵长度，可以为未来水厂迁址提供依据。