珠江河口水质模拟的二维水流盐度场分析

珠江河口水质模拟的二维水流盐度场分析

近年来，随着河口地区经济的快速发展，河口淡水资源开发利用的不断发展，世界各国的一体化、海拔高度、降水时间和空间变化、流域建设等自然因素的综合影响，河口淡水资源的入侵呈现出普遍加剧的趋势。由于盐水入侵对社会生产和生活带来的影响和经济损失愈显严重,其已经成为河口进一步发展的主要障碍。盐水入侵还对河口地区泥沙絮凝沉降和最大浑浊带具有极其重要的影响,影响航道的正常营运和维护。因此,面对这种愈演愈烈的盐水入侵现象,全面深入了解其变化规律以探讨影响其变化的主要动力因素具有重要的理论和实际意义。盐水入侵是潮汐河口特有的自然现象。国外早在20世纪50年代就开始了关于河口盐水入侵问题的研究,相关研究在早期多应用河口环流和盐淡水混合机制分析,随后在此基础上逐步发展采用数值模拟方法和其他较先进的方法。国内对于盐水入侵的研究始于20世纪80年代初对南水北调工程的影响论证,以及上海市饮用水和工业用水的需要,因此相关研究主要集中在长江口。事实上,在珠江三角洲地区,由于近年来频繁的河道挖沙等人类活动的影响,盐水入侵已日趋严重,尤其是近10年来,珠江河口咸潮更显其高频度、长时间和远距离上溯的特点,对珠江流域地区的社会、经济和环境的发展产生了极大影响。另外,由于珠江河网纵横交错,水源地密布,咸潮侵袭显得尤为复杂。盐通量是构成地球系统中物质循环的重要环节,也是陆海交互作用过程的重要影响因子。对于“通量”的概念,沈焕庭将其定义为:在一定时间内通过某一面积的某种物质的质量或者体积,其中物质可以是水、泥沙或化学元素等。随着水库拦沙的作用日益显著,许多学者对珠江的水沙通量进行了不少重要研究,而对于珠江口盐通量的研究并不多见。珠江三角洲水系错综复杂,动力因子众多,为全面掌握盐水入侵的具体过程及盐通量的大小,需借助于数学模型,其不但可弥补统计上的不足,而且能够提高通量估算的精度。近年来,珠江三角洲盐水入侵的机理分析和规律研究工作已经展开了不少,其多数侧重于动力条件与咸潮上溯特征的分析,而对盐水入侵的数值模拟计算极其少见;再者,研究多偏重于珠江口的部分河口区(磨刀门水道),而面向整个珠江河口区,尤其是将珠江河网区与河口区视作整体的盐水入侵模拟研究则极少。实际上,珠江河口河道密布,八大出海口门的盐水运动规律各不相同。河网区、河口区与近岸海域之间水流相互贯通,多种动力交互作用剧烈,难以将其截然分开;此外,河网区的物质输送形式非常独特,其对物质的滞留、积蓄作用对于通量的确定以及对河口海岸环境的影响至关重要。因此,准确描述与预测盐分在河口中的输移及其通量过程是全面认识河口动力学过程、污染物迁移转化以及河口生态环境问题的前提与基础。龚政指出,对于不需要深入了解流场和盐度场垂向结构情况下的数值模拟,可采用平面二维模型。本文旨在模拟计算珠江口的盐通量变化特征,不需要深入研究其垂向变化,因此采用珠江整体平面二维模型对珠江口盐水入侵的盐通量进行计算与分析,并展示咸潮在不同时刻对珠江流域的影响程度。具体地,利用二维数学模型模拟珠江三角洲河口及网河区的枯季咸潮上溯过程,主要计算珠江八大口门盐通量的大小及分配情况,并以狮子洋为例初步探讨盐通量的沿程变化,分析大、小潮期间不同咸潮对3种不同类型用水在珠江三角洲可用范围的影响。希望该成果能为珠江枯季河口区的咸潮防治工作、珠江河口区的水资源规划与利用提供参考。1数值模型1.1平均水流速度、平均流速笛卡尔坐标系下的平面二维水流运动方程包括连续方程和动量方程。水流连续方程:水流运动方程:盐度扩散方程:(1)―(4)式中,t为时间;x、y为右手Cartesian坐标系;η为水面相对于计算基面的水位;h为计算平面下的总水深;u、v分别是x、y方向的垂线平均流速;Pa为当地大气压;ρ为水密度,ρ0为水密度;f=2Ωsinφ为Coriolis参量(其中Ω=0.729×10-4rad/s,为地球自转角速率,φ为地理纬度);为地球自转引起的加速度;S为源汇项,(us,vs)源汇项水流流速,为平均水深盐度;Fs水平扩散项。1.2模型的计算范围与网格基于目前所收集的珠江三角洲网河区及口门区水文基础资料,将模型上游控制边界分别取在石咀(潭江)、高要(西江)、石角(北江)、老鸭岗(广州上游)和博罗(东江),同时考虑到洪枯季节外海潮汐作用的影响,模型外海开边界设定在珠江口外南海-30m等深线附近,将珠江八大口门、整个伶仃洋及黄茅海包括在内,八大口门分别为虎门、蕉门、洪奇沥、横门、磨刀门、鸡啼门、虎跳门和崖门,前4个口门统称为东四口门,后4个口门统称为西四口门。模型的计算范围与网格见图1。网格采用三角形与四边形相结合的混合形式,网格步长为50～2000m不等,网格数共计56911个,节点数共计55454个。1.2.1模型的计算方法模型驱动力主要包括上游流量、外海潮波和风应力,还应综合考虑径流、潮流以及季风等因素对盐水入侵过程的作用。其中径、潮流是最主要的控制因子,而季风虽然对盐通量的计算有一定影响,但由于缺少风速、风向资料,而且本文的计算时段并非在夏季,出现对咸潮上溯影响较大的台风暴的可能性很小。因此,本文的计算过程未考虑风的影响,只考虑上游径流和外海潮波的驱动作用。该模型计算时间步长取10s,初始条件水位取为0;盐度初始条件的给定对盐度模型的计算结果影响很大,一般情况下盐度扩散在3个月以后基本达到稳定,因此本模型中的盐度初始场采用根据循环计算3个月所得的盐度结果,以保证计算结果具有一个较为可靠的依据。图1-a显示了虎门处某点在计算过程中最后1169h的盐度变化情况,可见经过3个月的模拟计算,盐度扩散已经趋于稳定。图1-b为计算所得盐度初始场。模型上边界条件为给定实测水位,外海开边界潮位由南中国海潮波模型提供。模型上游边界盐度取为0;根据外海观测的统计资料,模型在西南开边界从陆地到外海的盐度值为20‰～33‰,东北开边界为25‰～33‰,东南开边界为30‰～33‰。由于模型较大,需考虑柯氏力影响,为了便于计算,模型中纬度统一取为23°N。曼宁系数是主要体现河床阻力的参数,由于网河区与口门处地形相差较大,所以曼宁系数也需分区给定。本模型中曼宁系数取值范围是0.015～0.030,从上游到口门逐渐增大,口外最大;伶仃洋内平均为0.016,网河区上游马口段为0.028,北江上游三水段为0.030,网河区中段平均为0.025。等盐度扩散均采用水流涡黏方程乘以尺度参数的形式,比例参数取为1.1。经调试,涡黏系数取为0.28。1.2.21年枯季补充验证考虑到盐水入侵主要发生在枯季,因此模型主要对枯季情况进行率定和验证,外海流速流向采用2007年中水进行补充验证。(1)2001年枯季:2001-02-07T14:00―02-16T10:00;(2)2007年中水:2007-08-16T10:00―08-17T16:00。验证站点包括网河区及口门区共34个站点(图2),部分验证结果见图3。该数值模式验证良好,能较好地反应盐水入侵的时空分布规律。2计算2.1小潮期间的盐通量在该数值模式验证良好的基础上,可计算得到指定断面连续潮周期内的盐通量。其断面盐通量情况能较好地反映珠江口的盐量变化情况,能为以后在珠江流域拟实施的“挡咸蓄淡”工程或运用水库调度“以淡压咸”等方案提供重要的科学依据。本文模拟2001-02-07―16河网与河口区的水流盐度输移过程,统计分析了大潮(02-08)与小潮(02-14)期间珠江各口门的盐通量。盐通量的计算公式为:式中,Fs为盐通量,Q为断面流量(m3/s),S为断面平均盐度,t1为涨(落)潮流起始时刻,t2为涨(落)潮流终止时刻,ρ为水体体积与重量的换算系数(kg/m3)。本文以某断面涨潮起止时刻来计算,因此,下文所述盐通量是指从某断面某次涨潮起始时刻到该次涨潮终止时刻通过该断面的盐分总量。虽然珠江各口门涨潮起止时刻不一致,但为了更直白地展示咸潮对整个珠江口的影响,本文将八大口门的盐通量相加近似得到珠江口的总盐通量。各口门的盐通量计算结果见表1,从表中可以看出,虎门和崖门为最主要的盐分输入通道。大潮间,珠江三角洲的八大口门总盐通量为1954×104t,其中大虎和官冲的盐通量最大,分别为811×104t和577×104t,各占到进入珠江口总盐通量的41.53%和29.51%。而在小潮期间珠江口总盐通量为512×104t,仅为大潮期间的26.2%,充分说明盐水入侵主要发生在大潮期间。小潮期间各口门的盐通量相对于大潮期间都显著减小,各口门的盐通量所占比例也发生了变化,其中虎门与崖门的盐通量之和占同时段总盐通量的94.97%,而磨刀门的盐通量几乎为零。可见,在小潮期间,虎门和崖门输入盐分的主导性仍然非常明显。这是因为盐水入侵主要是由于外海海水盐度大,内河径流盐度小,当潮波向里推进时,盐分随着水流向内陆扩散。由于不同口门水道的径流量与涨潮流量的相对大小不同,即山潮比不同,导致盐水入侵程度不一。山潮比越小,表明径流作用相对越小,潮流量相对作用则越大,盐通量也就越大;反之,山潮比越大,表明径流的影响相对越大,潮波作用被径流所压制,盐通量自然就越小。根据赵焕庭对珠江八大口门的特征分析可知,虎门和崖门皆以潮流动力为主,是八大口门中山潮比最小的两个口门,全年各月平均山潮比均小于1,其盐通量相对最大。而磨刀门作为西江最主要的泄洪输沙出口,山潮比最大,各月均大于1,多年平均为5.5,其中8月高达21.33,说明磨刀门以径流动力为主导,潮流动力很弱,故其对应的盐通量所占比例极小。2.2虎的盐通量随咸潮增长情况为了了解盐分从外海进入口门后在网河沿程的变化情况,本文对比分析了狮子洋沿程的7个断面(图4)的盐通量情况。从各断面盐通量变化数据(表2)可知,盐水从虎门进入狮子洋后,在一个连续潮周期内通过断面的盐通量沿程持续降低:一号断面取在和大虎平行的位置,其盐通量等于大虎的盐通量;二号断面在大虎上游5.89km处,其盐通量为大虎的94.78%;而当上溯至离大虎上游29.75km处的七号断面时,其盐通量减少至大虎的7.88%。可见盐通量沿程的减小趋势非常迅速而明显。盐通量随咸潮上溯距离的变化情况见图5-a。盐通量的沿程减少部分是由于珠江网河的分岔,使得部分盐分随水体进入支流,从而减少了狮子洋主要通道断面的盐通量(图5-a)。在涨潮期间,从大虎往上游方向,支流的分水作用使得水通量出现沿程减少的趋势,七号断面的水通量仅为大虎断面的16.26%。除此之外,盐通量沿程急剧减少也是断面盐度沿程降低所致。潮流在从河口上溯过程中,由于河床阻力的作用,其涨潮量会越来越小,混合结果使得盐度越来越低。图5-b是大潮和小潮期间各断面中心点盐度值的沿程变化,可知在任意断面大潮盐度都要高于小潮盐度,而且不论大潮还是小潮,盐度都随着离口门距离的增加而明显减少。2.3生活饮用水盐度限值及盐度分布根据世界卫生组织(WHD)标准和国家生活饮用水标准,结合国内水资源和用水状况,居民饮用水的氯化物含量不得高于250mg/L。盐度(S)与含氯度(cl)有如下关系:S‰=1.80655cl‰,换算可知生活饮用水的盐度限值为0.45‰。同理,由于灌溉用水和生态用水的含氯度分别不得高于1100mg/L和3150mg/L,因此,这两种类型的用水盐度限值分别为1.99‰和5.69‰。图6是大、小潮期间0.45‰、1.99‰、5.69‰3条盐度等值线的分布情况,说明盐水入侵已经严重影响到了澳门、珠海、中山、广州、香港、深圳和东莞等地的正常供水。对比图6-a、b可以发现大潮期间盐水入侵明显比小潮期间严重,3条盐度等值线都较小潮期间有明显上移。大潮期间,相应盐度限制等值线上溯至番禺、广州境内,而在小潮期间,盐度等值线则明显后退。3盐水入侵主要发生在枯季、潮涌、盐液、用水、生态用水等区域的通道综上所述,在枯季大潮期间,由于珠江入海径流量较小,使得盐水入侵较小潮期间明显严重。珠江河口枯季大潮时的总盐通量为1954×104t,而小潮时的总盐通量为512×104t,仅为大潮时的26.2%。无论是在大潮还是小潮期间,大虎和官冲的盐通量所占比例最大,说明盐水入侵主要是通过虎门与崖门这两个通道进入的。这是由于在珠江八大口门中,虎门和崖门的山潮比最小,潮波影响相对最大,