珠江三角洲的耦合斜压模型

珠江三角洲的耦合斜压模型

1几种类型的耦合计算方法的比较近年来，珠江三角洲的水环境和咸潮入侵日益突出。考虑到河流网络和河流入口区域的原因和影响，我们无法全面分析它们的原因和影响，也无法寻求基本的管理对策。因此，有必要发展珠江三角洲的全球水环境数学模型。针对此问题,国内不少学者将河网区、河口区以及近岸海域结合起来视作一个整体,开展了一维模型与多维模型耦合计算的研究,并已取得一定的进展,但其主要成果集中在一、二维模型方面,一维模型与三维模型连接计算的研究则较为少见。其中,彭静等、徐峰俊等的研究基本实现了珠江口一、二维模型的耦合计算,然而,珠江口由盐度梯度所产生的密度斜压效应明显,尤其是洪季时盐度层化现象异常显著,仅采用二维模型来描述其水动力特征远远不够,无法体现河口中盐度分层结构、密度环流及锋面结构等对物质输送的影响;包芸等的工作虽然实现了一维与三维斜压模型的耦合计算,但实际上并未涉及盐度的连接计算,无法动态反映盐水在河网区与河口区之间的回荡并沿河道上溯的过程;逄勇等开展了一、三维水质模型的耦合计算研究,但其物质输送连接条件简单,存在当一维向三维流动时一维模型无法满足三维模型分层边界条件的问题,这实质上是由于三维模型考虑分层结构而一维模型不考虑分层结构,两个模型信息传递不对称所产生的问题。从研究枯季咸潮入侵的角度来看,需将一维盐度模型与三维斜压模型相耦合,方能模拟盐水的入侵过程;从水环境管理的角度来看,构建整体斜压模型是建立珠江三角洲整体水环境数学模型的前提与基础。因此,本文在前人的基础上将一维盐度模型与三维斜压模型进行耦合计算,在盐度连接过程中,尝试通过延长三维模型上边界河段的方法来获得连接断面处盐度的三维分布信息,以此补充三维模型的分层边界条件,更为合理地构建珠江三角洲一、三维耦合斜压模型,从而为咸潮入侵、污染物通量及其环境效应等问题的研究奠定基础。由于珠江口夏季期间受强烈的太阳辐射影响,温度层化作用并不明显,盐度梯度比温度梯度对密度变化的作用更大,因此本文只考虑盐度梯度对斜压效应的贡献。2维耦合模型模拟范围涵盖整个珠江三角洲感潮河网区以及珠江口海域,如图1所示。河网区采用一维模型进行模拟,共概化河道299条,汊点189个,划分河道断面1726个,上游设置5个流量控制边界,分别为西江的高要、北江的石角、东江的博罗、流溪河的老鸦岗和潭江的石咀;八大入海口门虎门、蕉门、洪奇沥、横门、磨刀门、鸡啼门、虎跳门、崖门作为一维与三维模型计算的连接断面,即为耦合模型的内断面;河口区采用三维模型进行模拟,下边界延伸至外海,其中西边界取至下川岛西部的镇海湾,东边界取至红海湾,南边界取至70m等深线。三维模型采用能与物理边界吻合良好且分辨率较高的正交曲线网格系统,水平方向划分网格为183×186,垂直方向采用sigma坐标,均匀分为6层,网格系统如图1所示。三维模型外海开边界采用由美国NASA喷气推进实验室提供的TPMGDRB数据集高度计观测资料,以M2、S2、K1、O1这4个主要天文分潮进行驱动,以解决外海边界条件因外海潮位资料的缺乏而难以确定的困难。3数值模拟方法3.11河网水动力与盐酸盐的三维模型3.1.1x+zt,t明渠中的非恒定水流运动用一维圣维南方程组来描述,基本方程形式如下。连续性方程1B∂Q∂x+∂Ζ∂t=qL(1)1B∂Q∂x+∂Z∂t=qL(1)动力学方程∂u∂t+u∂u∂x+g∂Ζ∂x+gu|u|C2sR=0(2)式中:Z为断面水位;B为水面宽度;Q为流量;qL为旁侧入流流量;x为沿程距离;t为时间;u为平均流速,u=Q/A;A为过水面积;g为重力加速度;R为水力半径;Cs为谢才系数。对控制方程组采用四点偏心Preissmann隐式差分格式进行离散,采用三级联合解法求解。3.1.2河道断面个数河道方程∂(AS)∂t+∂(QS)∂x-∂∂x(AEx∂S∂x)-Se=0(3)河道汊点方程ΝL∑l=1(QS)l,j=(SΩ)j(dΖdt)j(4)式中:Ex为纵向分散系数;S为盐度;Se为外部源汇项;Ω为汊点水面面积;j为汊点编号;l为与汊点j相联接的河道编号;NL为汊点连接断面个数。一维河网盐度模型采用隐式差分迎风格式离散,方程的求解方法见文献。3.2静力弹性方程连续方程∂u∂x+∂v∂y+∂w∂z=0(5)x方向的运动方程∂u∂t+u∂u∂x+v∂u∂y+w∂u∂z-fv=-1ρ∂p∂x+AΜΔ2u+∂∂z(ΚΜ∂u∂z)+Fx(6)y方向的运动方程∂v∂t+u∂v∂x+v∂v∂y+w∂v∂z+fu=-1ρ∂p∂y+AΜΔ2v+∂∂z(ΚΜ∂v∂z)+Fy(7)垂向静力平衡方程ρg=-∂Ρ∂z(8)盐度守恒方程∂S∂t+u∂S∂x+v∂S∂y+w∂S∂z=AΗΔ2S+∂∂z(ΚΗ∂S∂z)(9)式中:u、v、w分别为x、y、z方向的速度分量;AM、KM分别为水平和垂向涡动黏性系数;f为科氏力系数;P为压力;ρ为海水密度;g为重力加速度;Fx、Fy分别为x、y方向上的外力;S为盐度;AH、KH分别为水平和垂向涡动扩散系数;Δ2为水平拉普拉斯算子。控制方程的具体差分形式及求解步骤参见文献。41.3三维结合斜率模型的构建4.1基于连接条件和耦合计算方法的连接条件4.1.1模型求解与显式耦合一维与三维模型连接断面为八大入海口门。设一维模型在第i个口门的连接断面上的水位、流量、盐度分别为Z(1)i、Q(1)i、S(1)i,三维模型在第i个口门的连接断面上平均水位、流量、平均盐度分别为Z(3)i、Q(3)i、S(3)i,(i=1,2,3,4,5,6,7,8),该断面位于第j列、第k行、第l层的网格点的盐度为S(3)ijkl。考虑联解时,一维和三维的连接断面被视为内断面,因此每个连接断面需补充水位连续、流量连续、盐度输送连续三个关系式,联解的思路如图2所示,采用显式耦合,三维模型将水位传递给一维模型,一维模型将流量传递给三维模型,盐度传递的物质通量包括由移流作用及扩散作用所产生的净通量,其中移流通量传递取决于连接断面水流方向,涨潮时由三维模型向一维模型传递,落潮时由一维模型向三维模型传递,扩散通量传递取决于连接断面处的浓度梯度。由于三维模型在连接断面处含有多个计算水点,因而需获悉连接断面处盐度沿宽度和水深方向的变化,方能补充其分层边界条件。连接的基本条件如下。水位连续Ζ(1)i=Ζ(3)i(10)流量连续Q(1)i+Q(3)i=0(11)盐度输送连续:(1)当水流方向为从三维到一维时S(1)i=(Q(3)i⋅S(3)i+FDΙFi)/Q(1)i(12)(2)当水流方向为从一维到三维时S(3)ijkl=(Q(1)i⋅S(1)i+FDΙFi)/Q(3)i⋅aijkl(13)式中:Q(1)i、Q(3)i以流向连接断面为负;FDIFi为由扩散作用而产生的通量;aijkl为第i个口门的连接断面上盐度沿宽度和水深方向的分布系数。当一维模型断面浓度高于三维模型断面浓度时,扩散通量由一维模型向三维模型传递,此时FDΙFi=Ai⋅Ex⋅S(1)i-S(3)iΔn,反之,FDΙFi=Ai⋅AΗ⋅S(3)i-S(1)iΔn,其中Ai为第i个口门的过水面积,n为法向。4.1.2计算流水管理体制具体耦合计算过程如下:在第n个时层首先进行三维模型斜压计算,将三维模型计算得到的水位赋给一维模型,补充一维模型下边界条件;然后进行一维水动力与盐度模拟计算,将一维模型计算得到的流量赋给三维模型,作为三维模型的上边界条件;最后根据水流方向判断盐度传递方向,补充一维或者三维模型的盐度边界条件,进入第n+1个时层的计算;如此循环重复,直至计算结束。4.1.3计算网格点分布系数由于三维模型描述的是物质在空间三维中的运动和分层现象,在模型计算中要求提供边界上物质浓度的三维分布信息,而一维模型只能获得断面平均浓度,无法满足三维模型的需求,因此如何给出三维模型边界上物质的三维分布信息是一、三维连接模型的难题。本文就这一问题进行尝试性研究,具体分为3个步骤:(1)耦合率定后的一维模型与三维模型,进行诊断模式计算,即不考虑盐度的变化,获得各连接断面的流量过程;(2)将原三维模型计算网格的上边界向上游适当延长(延长后的水平方向网格变为183×208,见图3,垂直方向分层保持不变),使连接断面变为新网格的内部断面,以第一步中诊断模式计算得到的连接断面处流量作为延长网格后的三维模型的上边界条件,单独运行三维模型斜压模式,获取该流量过程控制下连接断面上各网格点的盐度,再分别除以连接断面的平均盐度,即可得到计算时段内各网格点的比例系数aijkl,作为分层的依据;(3)返回原一、三维耦合模型,进行斜压计算,盐度根据关系式(12)、式(13)进行传递与分配,即当水流方向为从一维到三维时,连接断面上各网格点按相应时刻下的分布系数aijkl进行分配。这一方法的依据是设想连接断面处的流量过程提前发生于上游河流,借此获知在该流量过程控制与外海潮汐作用双重影响下连接断面处盐度的分层情况。4.2u3000关于连接计算的验证主要从以下两个角度来考察一、三维耦合斜压模型连接计算的效果,一是比较模拟结果与实测值的拟合程度,即验证效果;二是分析连接计算是否显著改变模型的计算流场与盐度分布。从直观的物理过程考虑,无论是将两个模型单独计算还是连接计算,模型计算的范围和外边界条件都始终一致,假如连接计算保证了连接断面水力要素连续并趋于实际的变化过程,那么,从理论上讲,整个计算域内任意一点连接计算前后的计算结果都应该趋于一致。所以,只要满足了上述两个方面的验证,即可以认为连接计算是成功的。模型采用2001枯季的典型水文条件进行率定,采用1999年洪季的典型水文条件进行验证,具体时间为1999年7月16日～24日(共205h)、2001年2月7日～16日(共213h)。洪季时外海边界的盐度设为表层33.5PSU、底层34.5PSU,枯季时均为34PSU。洪季、枯季以各自预热运行60d后的结果作为正式计算的初始条件,时间步长为40s。糙率是影响一维水动力模型计算精度的主要参数,河网区不同流域河段的糙率分别为:北江流域三水以上取值0.026～0.045,三水以下糙率为0.015～0.030;西江流域糙率为0.020～0.035;口门段糙率为0.010～0.020。一维盐度模型的纵向分散系数取值5.0m2/s。三维模型底部糙率系数Z0取值0.02cm。在连接断面纵剖面上盐度分布系数aijkl的计算中,洪季时上游径流量大,口门附近基本上被淡水所控制,模型计算得到的盐度分布系数aijkl为0,而枯季时上游流量微弱,潮汐作用相对显著,因而盐水向河口上游入侵距离较长,盐度分布系数aijkl有一定变化,具体由模型计算得到。4.2.1洪季枯季盐度计算误差变化水动力计算结果验证采用河网区内的50多个水文站、八大口门水文站、珠江口海区数个潮位站的同步实测资料,图5、图6为虎门、磨刀门的验证结果。分析结果显示,洪季、枯季的水位、流量的相位偏差均在1h以内,洪季时80%河网区站点的水位振幅偏差、流量相对误差分别在20cm、20%以内,枯季时所有站点的水位振幅偏差均在20cm以内;在连接断面处,无论是洪季还是枯季,八大口门的水位计算误差均在12cm以内,总体平均误差为5cm;洪季时八大口门的流量相对误差均在12%以内;表1列出了1999年洪季八大口门分流比计算与实测的对比结果,可以看出,计算结果与实测的基本一致,绝对误差均在2%以内,东四口门与西四口门分流比的绝对误差更是只有0.07%。另外,受资料所限,洪季的盐度计算结果只采用位于内伶仃岛附近的C1站(位置见图4)实测资料进行验证。由图7可以看出,模型模拟的盐度变化趋势与实际过程较为吻合,表层、中层、底层计算偏差依次为0.2、1.7、2.6;底层盐度计算结果较实测的偏低,可能主要是由于模型所采用的垂向分辨率较低,垂向过于均匀化所致。枯季时潮汐作用明显,表2列出八大口门盐度计算与实测的对比结果,计算结果与实测的基本一致,平均偏差为1.3。图8反映了蕉门、麻涌站受潮汐影响的盐度动态变化过程,模型所刻画的变化过程与实测基本一致。综合上述验证结果,表明本文所建立的耦合斜压模型能再现水流、盐度的实际变化过程,模拟的水位、流量、盐度与实测值吻合较好,偏差均控制在合理范围以内。4.2.2连接计算与单独计算的比较为进一步考察连接计算的效果,以1999年洪季实测流量过程作为八大口门的边界条件单独运行三维斜压模型,将其计算结果与耦合斜压模型计算结果进行比较,分析连接计算对三维模型计算流场、盐度分布的影响程度。由于单独计算时八大口门的流量采用实测值,而连接计算时八大口门的流量是通过计算得到,与实测结果不可避免会存在偏差,因此连接计算的结果与单独计算的结果必定存在一定的偏差。但从理论上讲,连接计算与单独计算的偏差越小,说明连接计算越理想。从理论上讲,越靠近连接断面的区域,受连接计算影响越显著,所以将重点讨论连接断面前后一段范围内的计算结果。一维部分选择8个比较断面O1～O8,均与各连接断面相隔一个断面,比较其水位、流量;三维部分选择8个比较点T1～T8,位置分布如图4所示,比较其水位、流速、盐度。受篇幅所限,仅列出部分一维、三维的比较结果,见图9～图11。从水动力比较结果来看,连接计算与单独计算的水位与流量的相位差别均很小,两者一维部分的水位差为5.5cm,流量相对偏差为5.5%,三维部分的水位差为4.0cm,流速相对偏差为3.0%。通过与验证效果的比较不难看出,这种差别主要由连接断面处的