变化环境下万泉河河口区盐分扩散的模拟分析

海来以波，施浩然，刘晓庆，刘 蛟，方明仪，张蓝月

(1.西华大学能源与动力工程学院，四川 成都 610039；2.西南科技大学环境与资源学院，四川 绵阳 621010)

河口区域在内河和外海因素交互影响下形成特殊的水生态环境，是洄游鱼类产卵、育幼和索饵的重要场所[1-2]。随着上游水资源开发带来的内河来水条件变化[3-4]，以及海平面上升等带来的河口潮位波动[5]，导致外海潮波带来的盐分在河口区的扩散发生变化，水体盐度的改变对河口区水生态环境产生了重要的影响。张洪亮等[6]对莱州湾的研究说明盐度升高引起了浮游植物生物量降低、浮游生物的种类减少，同时对海洋生物产卵场造成了不利影响。除了径流和潮位2个主导因素外，河口区盐度尤其是横断面上的分布同样受到河床地形的影响[7]。林一楠等[8]通过对钱塘江河口沿程7个测站的实测数据分析，揭示了盐水主要沿深槽上溯、凹岸盐度及其变幅大于凸岸的客观规律。

河口区的盐分扩散是由内河来流、外海潮位和河床地形三者共同作用下形成的水动力条件决定的。因此，基于水动力学模型的研究成为河口区咸潮入侵防治以及生态环境修复的主要手段。Duong等[9]采用Mike21模型，通过一维和二维的耦合模型模拟了越南湄公河三角洲一条主要支流在未来不同情景下(如海平面上升、上游来流量变化等)的海水入侵变化。Uncles等[10]对强潮作用下英格兰南部Lynher河口径流和海水入侵相互作用的季节变化进行了研究，结果表明海水入侵主要受潮汐和径流的控制。赵高磊等[11]以咸潮上溯和冲淡区为研究对象，分析河口上游修建水库后，在丰水期、平水期、枯水期径流改变的情况下，北门江河口地区咸潮上溯和冲淡区变化规律。Xu等[12]基于Mike21分析了河道排放条件、潮汐范围和海水入侵之间的关系，构建了长江口的咸水入侵函数。陈浩等[13]基于Mike21 FM 建立了研究区域海水入侵模型，对辛集挡潮闸不同下泄流量时原始河道和清淤后河道的海水入侵进行模拟研究，以论证所估算生态基流的合理性。

万泉河为海南省第三大河，近年来由于发展需要，河道水库的修建和运行改变了河口区原有的水文条件[4,14]，同时中国近海区域海平面的上升也影响了万泉河河口潮位变化。径流和潮位使得万泉河河口区盐分扩散的变化，直接影响到了对盐度极为敏感的国家二级保护鱼类花鳗鲡的生境条件[15-16]，同时也对河口区域内的红树林生态林以及农田生态系统造成影响，沿海农田土壤盐渍化问题更加凸显[17- 18]。为此，本文基于Mike21模型，开展径流和潮位变化情景下，不同水文时期河道内盐度分布模拟，分析万泉河特殊的河网结构下，不同季节中径流和潮位变化对海水潮波上行距离、影响范围以及强度等方面的具体作用，阐明不同条件下万泉河河口区的盐分扩散变化规律，为河口区的水环境保护、水生态建设提供依据。

1 研究区域概况

万泉河为海南省第三大河流，发源于琼中黎族苗族自治县五指山，经万宁市至琼海市，在博鳌港汇入南海，干流全长157 km，平均坡降1.12‰，总落差586 m。万泉河流域位于109°37′E～110°38′E和18°46′N～19°31′N，整体形状呈带状，平均宽度15 km，面积3 693 km2。本文对外海潮波上行中盐分在内河河道扩散进行研究，因此根据多年来万泉河河口区潮波影响范围的潮区界确定本次研究的范围为入海口至嘉积镇大坡村，河流长度约为16 km。研究区河道主要流经嘉积镇和博鳌镇，属于国家重点开发区，是海南省生态空间格局中以旅游业为主题的东部核心区，也是亚洲论坛永久性会址所在地(图1)。

图1 研究区域

万泉河河口区的河网结构复杂，尤其是在河流入海口附近的岛屿众多，干流被主要的沙坡岛和东屿岛分为北侧、中侧和南侧三支，与河口联成一体，形成较大水域，而该水域与外海由横亘的海沙与河沙冲积而成的玉带滩隔断，仅在北端有一个小的入海口，为此在万泉河河口形成特殊的葫芦形港湾。

万泉河流域河口的潮汐类型属混合潮，每一太阴日内有时发生2次高潮和2次低潮，有时只出现1次高潮和1次低潮，潮汐日不等现象显著。根据河口处附近港北潮位站的监测数据，按国家85高程基面，河口的多年平均高潮位和低潮位分别为 0.96、0.21 m，潮汐差0.75 m。位于研究区上游约5 km的嘉积水文站为万泉河下游的控制性水文站，监测的多年平均流量为147 m3/s。万泉河河口区上游的天然来流量年内分配极不均匀，汛期6—11月占了全年径流量的80%左右，9—11月的主汛期占60%左右，但是随着上游水资源的开发，经上游水库调节过后，河口区的年平均来流量和年内分布差异均有所减少。

2 数值模拟

2.1 Mike21模型

由于本次研究河口区域水深较浅，加上纵向水流对盐度分布的扰动，盐度分层不明显，因此本次研究根据收集的数据，采用丹麦DHI公司开发的Mike21模型对万泉河河口区在2017年3个典型水文月中日尺度上的盐度分布开展模拟。Mike21主要包含二维水动力模块、波浪模块、简单对流扩散模块、水生态模块以及泥沙运移模块等，在河口、海湾和海岸地区的水动力及其相关现象的平面二维仿真模拟得到广泛使用[11,13,19-20]。Mike21水动力模型遵循纳维耶-斯托克斯(Navier-Stokes)方程，并服从Boussinesq假设，即流体低速流动中，密度的变化不考虑由压强变化所引起的，仅考虑温度对密度的影响，同时服从静水压力。Mike21的控制方程组如下：

(1)

(2)

Mike21模型采用有限体积法，对研究区进行三角网格化离散处理。本次模型划分为7 690个网格。网格最小角度29°，最大角度59.9°，共计4 388个计算节点，计算的最小时间步长取1.5 s，最大时间步长取200 s，CFL(Courant-Friedrichs-Lewy)数为0.8。

2.2 模型数据

按流量大小，对万泉河河口区2017年5、2、4月3个不同水文月的盐度分布开展模拟，模型上边界条件的驱动数据为模拟期内嘉积水文站的逐日流量数据；下边界条件的驱动数据为模拟期内港北潮位站的逐日潮位数据。降水、蒸发、风速风向等气象数据来自于国家气象科学数据中心(http://data.cma.cn/)的琼海气象站实测数据。

模拟河道的河床组成主要为中粗砂，床面不平整，采砂形成的局部深坑散乱分布，滩地较宽阔生长有杂草和低矮灌木，两侧边岸一般为茂密灌木和椰子、槟榔等高大乔木林带，河道两岸为较为光滑的堤防工程，根据河道现状，结合河床地形分布，研究区域的糙率取为0.025～0.030。

研究河段内无实测的水文数据，且考虑不同径流和潮位条件下的盐度分布，而盐度分布是盐分在水动力学作用下的结果，因此直接根据汀州水质监测站的监测断面在各月份的实测盐度数据开展模型验证。

2.3 模拟情景

首先根据2017年的实测数据，对代表性水文月的盐度模拟结果进行率定和验证，并将其作为参照模型。在此基础上，考虑上游来流变化和河口潮位变化2种情景，通过改变参照组的输入条件，设置对比模型。具体情景设置如下。

a)参照组模型。万泉河河口区多年平均流量为147 m3/s，该流量对河道的塑造以及稳定生态环境的形成具有重要的意义。同时根据水文的季节性变化，结合上游的水库调节作用，选择了2017年5、2、4月，开展模拟与验证。3个月的平均流量分别为178、148、106 m3/s，平均高潮位分别为0.96、0.93、0.96 m，平均低潮位分别为0.15、0.12、0.18 m。

b)流量对比组模型。考虑研究区域上游在水资源开发中的取水增加后，在参照组的基础上，分别假设河口区的来流减少10%和20%，而潮位不变的2种情景，对河口区盐分扩散进行模拟，分析内河来流对万泉河河口区盐分扩散的影响。

c)潮位对比组模型。基于中国近海近50年海平面变化的周期性、趋势性等规律，预测研究河口区的海平面可能具有较大幅度的上升[21]。因此，在参照组的基础上，假设潮位增加0.1、0.3 m，上游来流量不变的2种情景，对河口区盐分扩散进行模拟，量化分析外海潮位对万泉河河口区盐分扩散的影响。

3 模型验证与结果分析

3.1 模型验证

根据3个月份模拟的平均盐度与实测数据的对比(表1)看出，3个月的盐度误差不大，最大的5月为-13.9%，4月的误差仅为3.3%。该结果说明构建Mike21模型能够较好地反映出万泉河河口区在不同径流和潮位作用下的盐分扩散情况。

表1 汀州断面盐度的模拟与实测结果对比

3.2 盐分沿河扩散距离变化

本文参照相关研究以河道内盐度浓度为0.5‰位置为外海盐分上溯扩散顶点对上溯距离进行统计[11]，得到了不同情景下万泉河河口区咸潮沿河上溯的影响范围，结果见表2。总体上来看，上游来流量的减少以及外海潮位的上升均会对盐分沿河扩散距离产生正向影响，且来流变化的影响表现更加明显。在不同月份，当流量发生同等幅度的减小时，流量偏丰月份的咸潮上溯距离变化更大。当潮位增加0.1 m时，就会对来流量为106 m3/s的4月份盐分上溯距离产生影响，而在来流量为148 m3/s的2月时，潮位需要增加0.3 m才会对盐分上溯距离产生比较明显的影响；当来流量增大至178 m3/s的5月时，假设情景下潮位的增加对盐分上溯距离基本没有影响。通过不同情景下，盐分上溯距离的变化可以看出，当上游来流量偏枯的时候，外海潮汐变化的作用更加明显；而当上游来流量偏丰的时候，内河来流量变化的作用更加明显。

表2 万泉河河口区不同情景下盐分沿河扩散距离统计

3.3 河道内盐度分布变化

2017年4月万泉河河口区以及4种变化情景下河道内的盐度分布见图2，可以看出不同情境下，河道的盐度分布变化规律基本一致，在3条支流上都表现为盐度在0.8 km范围内快速地从4‰降低到0.5‰以下。这主要是受4月份上游来流量较少的影响，内河径流对入海口附近盐分扩散的影响不大。因此，在径流量减小的2种情景下，1‰以上的盐度分布区域基本没有变化，只是上溯距离分别增加了0.26、0.28 km。而潮波通过入海口进入内河后，首先较大的水域起到了缓冲作用，降低了潮汐作用；同时万泉河河口区在入海口附近河床平缓，而3条支流的河床地形往上游抬升明显，这进一步降低了潮汐对盐分的上行扩散作用，因此在潮位升高的2种情景下，1‰以上的盐度分布区域也基本没有变化，只是上溯距离分别增加了0.21、0.38 km。

图2 2017年4月万泉河河口区不同情景盐度分布

相对于4月来水量，5月万泉河的流量偏丰，其不同情景下河道内的盐度分布见图3，可以看出不同来流减少以及潮位升高后盐分上溯顶点(0.5‰)的上移变化，原有流量178 m3/s减少10%和20%后，外海盐分的影响范围分别增加了0.25、0.52 km；而随着潮位增加0.1、0.3 m，咸水上溯距离无明显变化。从盐度的空间分布变化上来看，总体上不同情景的分布规律较为类似，表现为在狭小的入海口附近区域盐度最高，均高于4‰，在继续沿河扩散中随着内河水域的扩大，盐度比较均匀地往上游扩散后逐步减小，然而继续向上游扩散进入分叉河道后，受弯道水动力学塑造的河床地形影响，在三支分叉河道中表现出不同的规律。首先在北侧支流河岸凸点下游，河道的流量大，流速大，该区域的水深较其他地区浅，盐度扩散快，因此盐度较低，为0.1‰～0.5‰，而在上游，受弯曲河道顶托作用，沿河流速减少，盐分在此处累积，呈现出比下游更高的盐度分布。同样在南侧支流，东屿岛与玉带滩对中间河道的束窄，盐分在潮汐作用下通过束窄区进入上游后，受到增强的内河径流影响，盐分在此处累积，盐度达到4‰左右。在经过分叉河道，进入汀州断面上游的主河道之后，河道内的盐度扩散又呈现出比较均匀的扩散，直至上溯顶点。

图3 2017年5月万泉河河口区不同情景盐度分布

4 讨论

通过对5种情景下万泉河河口盐度分布的对比模拟分析，可以看出，外海盐分在河口区的扩散影响范围受到了潮汐增加和来流减小的正向影响，这与目前基本的研究成果是一致的[10-11,22]。但是在万泉河口，不同程度的流量减小都会使得盐分扩散范围增加，但是潮位升高的影响还受到了内河径流的制约。当万泉河上游来流为106 m3/s，潮位增加0.1 m就能使得盐分扩散距离上移，而来流量增大到148 m3/s，潮位需要增加0.3 m才能使得盐分扩散距离上移。这说明径流和潮汐对盐分扩散范围的影响是相互影响、相互制约的。

研究发现，万泉河河口区特殊的水网形状以及河床地形等对径流和潮汐相互作用下形成的水动力条件具有决定性的作用，这在很大程度上决定了河口区的盐度分布。这与Parsa等[23]得出的河口的盐度入侵在很大程度上取决于河口的形态、地形、水动力和潮汐特性的结论相一致。万泉河在出海口形成的水域对潮汐起到了缓冲作用，使得大部分盐分聚集在此区域，盐度分布迅速降低；加之岛屿的阻隔与河床抬升的影响，这都削弱了潮汐对盐分的扩散影响作用。这就可能造成了模拟结果中潮位变化情景下的盐度分布改变弱于径流变化。同时由于河道的弯曲和束窄变化，在局部形成特殊的水流特征，出现盐分在局部累积的特殊分布情况。

本文从对河道环境等具有重要意义的流量出发，考虑一定的水文变化，在不同的来流和潮位变化情景下，对万泉河河口在2017年3个典型水文月份上的盐度分布情况进行了模拟分析，阐明了在特殊的水系结构及河道地形条件下，径流和潮位变化对万泉河河口盐度分布的影响作用，可为河口区的水生态环境建设提供科学依据。