瓯江下游河段污染物质滞留时间数值模拟研究

李谊纯

（河海大学海洋灾害与防护教育部重点实验室，南京 210098；广西科学院，南宁 530000）

瓯江下游河段污染物质滞留时间数值模拟研究

李谊纯

（河海大学海洋灾害与防护教育部重点实验室，南京 210098；广西科学院，南宁 530000）

利用平面二维数学模型对瓯江下游守恒性污染物质滞留时间的空间变化及其对于径流、潮差、初始排放时刻的响应进行了研究。研究将瓯江梅岙至口门河段分为8个子区域分别进行数值试验。数学模型上、下游边界分别采用径流量和潮位控制。上游径流量选取5个代表性流量，下游选取大、中、小潮及混合潮型进行组合研究滞留时间的基本变化规律。研究认为，瓯江下游河段污染物质滞留时间与径流量呈现极好的幂函数关系，与潮差呈现良好的线性函数关系。由于径流量的变幅远大于潮差的变化，综合而言，滞留时间主要是径流量的函数。污染物的排放时刻对其在河口里的平均停留时间影响巨大，低潮位时刻排放，污染物在河口的平均停留时间大为增加。

滞留时间；数值模拟；瓯江

Biography：LI Yi-chun（1977-），male，engineer.

滞留时间是水体内物质更新速率的基本度量，滞留时间可以将河口生态系统研究的许多方面统一起来，因此反映了水动力过程和生态过程间最根本的联系。滞留时间是水体微团或其他要素如盐、污染物等从其进入某一水体至被输运到水体以外滞留在水体中的平均时间，在稳态条件下，可以根据河口内该变量与其在边界处的交换速率的比值来估算。由于空间差异和河流过程的时间依赖性及许多重要物质的非线性行为，使得滞留时间十分复杂［1-2］，所以仅了解滞留时间的平均值是不够的，还应考虑其时空变化等特征。Awaji，Signell等使用水体示踪法对潮汐交换机理进行了有益的探讨［3-5］。漂流物示踪技术［6］的进步为记录河口拉格朗日运动和交换提供了前所未有的机会。随着采样技术的发展，精确的示踪物将为流体交换过程提供更精确的信息。数学模型的优化和更全面的验证使得模型方法成为评价河口净化速率的重要途径。Zimmerman［2］给出的滞留时间的定义为：水体中某一个物质的微团的滞留时间为它到达水体的出口前在水体中的停留时间。可以看出，滞留时间也是时间和空间要素的函数。不同空间位置的物质微团在不同时刻排放将会导致其在水体中具有不同的滞留时间。Takeoka［7］给出了一个滞留时间的定义和计算方法，为通过数学模型手段研究滞留时间问题奠定了基础。随着河口环境问题越来越引起重视，河口物质输运的时间尺度方面的研究在近几年也取得了很大进展［8-10］。本文构建了一个平面二维数学模型并将其应用于瓯江口研究各河段滞留时间及其对于径流、潮汐污染物初始排放时刻等因素的响应。瓯江下游是我国东南沿海重要的经济区，该地区经济高速发展的同时，环境保护问题也日益凸显。因此，为有效控制瓯江口海洋环境污染以及对海洋的其他合法利用，开展瓯河口输运过程及物质交换的时间尺度的研究具有现实和长远的意义。

1 数学模型简介

为能够使建立的数学模型具备良好的河道地形、边界的适应能力，模型方程采用曲线正交坐标系下的控制方程组。控制方程采用结构化的曲线正交网格作为计算网格，变量布置采用交错网格。空间离散采用角输运迎风格式（CTU，Corner-Transport Upwind）并结合TVD限制器（Van leer）进行通量限制，源项采用算子分裂算法处理，时间积分采用可保持TVD性的两步格式计算［11］，方程在时间和空间上精度均基本达到2阶。

瓯江河口数学模型的计算范围上游边界位于瓯江干流圩仁站附近，采用流量控制，外海边界取在飞云江口—南麂—坎门一线，边界给定潮位过程。计算采用曲线正交网格总共包括230×236个网格，网格步长为20～1 000 m。梅岙下游的地形数据为2005年6月的实测数据。验证资料为2005年6～7月的同步测量水文资料。经验证模型能够复演瓯江口水动力过程。

2 瓯江口滞留时间计算

2.1 滞留时间计算方法

文中滞留时间的计算采用Takeoka方法［7］。假设在t=0时刻，水体里的物质的总量为R0，在t=τ时刻仍留在水体里的物质的量为R（τ）。即R（τ）为滞留时间大于τ的物质的量。定义如下函数

那么，平均滞留时间τr可以定义为

分部积分并考虑到lim R（τ）=0，则有

2.2 计算方案

瓯江下游地形复杂，自楠溪江以下被七都岛分为南、北两汊，龙湾下游又被灵昆岛分为南口和北口。滞留时间的研究范围上游至梅岙，下游至口门。为了计算各河段的平均滞留时间，将梅岙以下至口门河段分为8个部分分别予以考虑（图1）。由于南口封堵造陆，南口河段的滞留时间计算未包含在内，但积分区域包括南口河段。分别选取了Q80（85 m3/s）、Q5（1 200 m3/s）代表枯水和丰水流量，多年Qm（470 m3/s）代表平水流量。为能尽量准确拟合流量与滞留时间的关系，增加Q50（235 m3/s）和Q20（560 m3/s）2个流量的计算。各流量均根据上游瓯江圩仁站流量资料得出，下游边界取中潮潮位过程。在关于潮差影响的计算中，上游选取多年平均流量Qm，下游分别组合2005年6月～2005年7月实测大、中、小潮资料，平均潮差分别为4.86 m、3.56 m和2.46 m。污染物质分别在高潮位和低潮位时排放以研究初始排放时刻对于滞留时间的影响。

图2以S2河段为例给出了在平水流量组合中潮情况下的滞留时间及积分区域内的所滞留的物质总量的演变过程。可以看出，随着计算过程的发展，河口内滞留的（污染）物质的总量逐渐减少，相应的，滞留时间曲线的坡度愈趋平缓。至积分临界条件满足，则认为河口内滞留的物质绝大部分已冲刷至口外，由式（3）计算的物质停留在河口内的平均时间即为该情况下的滞留时间。

图1 滞留时间计算河段划分Fig.1 Segmentation of Oujiang estuary in numerical simulation

图2 平均停留时间及水体内滞留的物质总量的演变Fig.2 Cumulative average residence time and residual passive dissolved conservative matter

2.3 径流量的影响

图3为沿程滞留时间，可以看出各方案的滞留时间的沿程变化情况。在多年平均流量下梅岙段（S1）的滞留时间约为 244 h，向下游依次为 215 h（S2）、170 h（S3），乌牛段（S4）为 108 h，状元段（S5）为 162 h，再向下游的龙湾段（S6）为 132 h，七里段（S7）为 34.8 h，黄华段（S8）为 15.2 h。总体而言，滞留时间从上游河段向下游河段逐渐减小。其他流量条件下滞留时间的沿程变化总体趋势与多年平均流量情况下类似。值得指出的是，尽管S6河段位于S4河段的下游，但其滞留时间较S4段滞留时间长。为此在上游边界取Qm，下游组合中潮的情况下将S6段拆分成的南、北2个区域，2个区域的滞留时间分别约为174 h和76 h。可以认为导致该现象的原因可能是S6段的南岸部分受南口动力的影响较大，受同样影响的还有S5段，所以尽管S4和S5距离河口口门距离相近，但S5河段的滞留时间较S4河段大。可以认为S6河段的北部区域受灵昆北口及七都北汊的影响较大，而南部区域则受灵昆南口及七都南汊的影响较大。

随着流量的增加，各河段的滞留时间迅速减小，越向上游，减小幅度越大，梅岙段在Q80时，其滞留时间可达907 h，而在Q5时仅为119 h，二者相差近800 h。而最靠近口门的黄华段，相应的滞留时间分别为52.3 h和6.5 h，二者相差约45 h，由图3可知，越向下游，滞留时间减少越小，受径流的影响也越小；越向上游，减少越大，受径流的影响也越大。

定义参数

式中：Rt5%、Rtm分别为径流量取Q5和Qm的滞留时间。图4为大潮和中潮情况下各河段的R1值。可以看出，一方面不论下游边界采用大潮还是中潮，这一比值均表现为上游大，而下游河段较小，说明径流的影响在上游河段大于下游河段；另一方面，随着潮差的增加，比值相应减小，说明在不同的潮汐动力下，滞留时间对于径流量变化的响应的敏感程度是不同的，大潮情况滞留时间随径流量的变化小，而在中潮情况滞留时间随径流量的变化较大，河流动力和潮流动力对河口的冲刷能力不是简单的线性叠加关系。图5给出了梅岙至口门各河段的滞留时间随流量变化的关系。由其拟合关系式和R2值，可以看出，瓯江口各河段的滞留时间相对于径流量呈现幂函数关系，且拟合极好，R2均大于0.99。

2.4 潮差的影响

图6为上游选取多年平均流量，下游组合大、中、小潮的情况下，沿程的各河段的滞留时间的变化。可以看出，随着潮差的增大，各河段的滞留时间明显减小。最上游的梅岙段在外海边界为大潮情况下，其滞留时间约为208 h，而小潮情况下约为272 h。龙湾段大潮时为105 h，小潮时约为160 h。灵昆以下的七里段和黄华段变化的绝对值较小，约为2～13 h。自龙湾段以上除乌牛段外，滞留时间均超过120 h。一般来讲，不论大、中、小潮，任一潮型均不会恒定持续很久时间，为了更好的研究实际情况下的各区段的滞留时间状况，选取外海边界为大、中、小潮组合做为外海的计算边界条件，上游取多年平均流量进行数值试验。可以看出各区段的滞留时间与外海边界取为中潮情况下非常接近。

图3 不同流量下的滞留时间Fig.3 Longitudinal variance of residence time for different river discharge

图4 大、中潮情况下各河段的R1值Fig.4 R1of different segments for spring and middle tide

图5 滞留时间与流量关系Fig.5 Relationship between residence time and river discharge of different segments

图6 不同潮差下滞留时间的沿程变化Fig.6 Longitudinal variance of residence time for different tidal range

图7 滞留时间与潮差的关系（Qm）Fig.7 Relationship between residence time and tidal range of difference segments for mean discharge

图8 各河段滞留时间与Ur/Ut关系Fig.8 Relationship between residence time and Ur/Utof different segments

图7为径流取Qm时各区段的滞留时间与潮差的关系，从各河段进行的数值试验的计算结果来看，滞留时间相对于潮差基本上呈线性关系，潮差越大，滞留时间越小。由于每个河段仅为3个数据点（大、中、小3个潮型），滞留时间与潮差是否存在更复杂的关系，有待进一步研究。但仅就现有的数值模拟结果来看，滞留时间相对于潮差二者呈线性关系是可接受的。图8为各河段的滞留时间随Ur/Ut的变化。其中Ur、Ut分别为径流速度和潮流最大速度，均取各河段中间位置的平均值。可以看出滞留时间与Ur/Ut仍呈现良好的幂函数关系，由于径流量的变幅远大于潮差的变化，瓯江口径流洪枯比可达2 000，而潮差变化一般不会超过3倍，综合而言，滞留时间主要是径流量的函数。

2.5 初始排放时刻的影响

前述滞留时间计算中，污染物质的排放均选在高潮位时刻，即河口内水体体积最大的时刻。瓯江河口内潮波以驻波为主，兼有前进波的性质，呈现一种混合型式。潮流转流时刻均在高、低平潮时刻附近。若是在高潮位时排放，所排放物质将很快随落潮流向口外输运。相反，如果在低潮位时排放，则污染物将会先随涨潮流上溯，其在河口中停留的时间必长于高潮位时刻排放的情况。为了研究瓯江口各河段污染物排放时刻对其在河口内停留的平均时间的影响，对比了水文控制条件采用多年平均流量组合中潮的情况下污染物在高潮位及低潮位排放后在河口内的平均停留时间，由图9可以看出，低潮位时排放将导致污染物质在河口里的平均停留时间大为增加。由梅岙段至最下游的黄华段，增加的时间从35 h逐渐增加至120 h。定义参数

式中：Rth、Rt1分别为污染物质在高、低潮位时排放情况下物质在河口内的平均停留时间。图10为各河段的R2值。可以看出，排放时刻对于下游区段的影响要大于上游区段。影响最大的黄华段增加约为9倍，而最上游的梅岙段仅增加了15%。

图9 不同排放时刻平均停留时间Fig.9 Relationship between average residence time and the release time

图10 不同排放时刻的平均停留时间比Fig.10 Sensitivity of average residence time to the release time

3 结语

瓯江口各河段滞留时间随径流量的增加而减小。愈向上游河段，滞留时间受径流影响越大；越向下游河段，滞留时间受径流影响越小。各河段的滞留时间与流量呈幂函数的关系。瓯江口各河段滞留时间随潮差的增大而减小，滞留时间随潮差呈线性变化。不同的潮汐动力下，滞留时间对于径流量变化的响应的敏感程度不同；不同径流量情况下，潮差对于瓯江口滞留时间的影响亦不同。流量较小时相对大流量情况下滞留时间随径流量的变化更为显著，河流动力和潮流动力对河口的冲刷能力并非简单的线性叠加关系。潮差对下游河段滞留时间的影响大于上游河段。瓯江口滞留时间主要是径流的函数，潮汐动力作用次之。瓯江口属强潮河口，污染物的排放时刻对其在河口里的平均停留时间影响巨大。低潮位时刻排放，污染物在河口的平均停留时间大为增加。

［1］Oliveria A，Baptista A M.Diagnostic Modeling of Residence Time in Estuaries［J］.Water Resources Research，1997，33：1 935-1 946.

［2］Zimmerman J T F.Estuarine Residence Times［C］//Kjerfve B.Hydrodynamics of Estuaries，Volume I：Estuarine Physics.Florida：CRC Press，1988.

［3］Awaji T，Imasato N，Kunishi H.Tidal Exchange through a Strait：A Numerical Experiment Using a Simple Model Basin［J］.Journal of Physical Oceanography，1980，10（10）：1 499-1 508.

［4］Signell R P，Geyer W R.Numerical Simulation of Tidal Dispersion around a Coastal Headland［C］//Cheng R.Proceedings of the International Conference on the Physics of Shallow Estuaries and Bays.New York：Springer-verlag，1988.

［5］Signell R P，Butman B.Modeling tidal exchange and dispersion in Boston harbor［J］.Journal of Geophysical Research，1992，97：15 591-15 606.

［6］ Hitchcock G L，Olson D B，Cavendish S L.A GPS-tracked Surface Drifter with Cellular Telemetry Capabilities［J］.Marine Technology Society Journal，1996，30：40-49.

［7］Takeoka H.Fundamental Concepts of Exchange and Transport Time Scales in a Coastal Sea［J］.Continental Shelf Research，1984，3（3）：322-326.

［8］Shen J，Lin J.Modeling Study of the Influences of Tide and Stratification on Age of Water in Tidal James River［J］.Estuarine，Coastal and Shelf Science，2006，68：101-112.

［9］Delhez E J M，Heemink A W.Residence Time in a Semi-enclosed Domain from the Solution of an Adjoint Problem［J］.Estuarine，Coastal and Shelf Science，2004，61：691-702.

［10］Sheldon J E，Alber M.A Comparison of Residence Time Calculations Using Simple Compartment Models of the Altamaha River Estuary，Georgia［J］.Estuaries，2002，25（6）：1 304-1 317.

［11］Leveque R J.双曲问题用的有限元法［M］.北京：世界图书出版社，2003.

Numerical simulation of residence time in Oujiang estuary

LI Yi-chun
（Key Laboratory of Coastal Disaster and Defence，Ministry of Education，Hohai University，Nanjing210098，China；Guangxi Academy of Sciences，Nanning530000，China）

A numerical modeling study to predict the residence time of a conservative tracer in Oujiang estuary was given in this paper.An integrated hydrodynamic-dispersion model was used to predict the average residence time in the estuary for various tidal range and river discharge conditions.The estuary from Meiao to the river mouth was divided into eight segments.Five representative river discharges and three tidal ranges were selected to give the computational conditions in the numerical simulations.Respectively，the residence time of each segment under different conditions was calculated to investigate the influences of drainage，tidal range and initial condition.The numerical results show that in Oujiang estuary the general regression relationships between the residence time and river discharge are power function，and linear function between the residence time and tidal range.Because of much greater variance of river discharge than tidal range，synthetically，the residence time is mainly affected by drainage.The initial condition affects the average residence time significantly.The average residence time for low tidal level is much greater than that for high tidal level.

residence time；numerical simulation；Oujiang estuary

X 522；O 242.1

A

1005-8443（2011）06-0434-06

2010-11-03；

2011-01-14

广西科学院基本科研业务费资助项目（10YJ25HY03）

李谊纯（1977-），男，河北省人，工程师，博士研究生，主要从事河口动力学研究。