湖泊三维水动力—水温耦合模型及其应用研究

康 玲，靖 争

（华中科技大学水电与数字化工程学院，430074，武汉）

湖泊是重要的城市水体形态，具有调蓄雨水、维持生物多样性、美化环境等功能。水温是影响水生动植物生长、繁衍和迁徙的重要因素之一。然而，随着气候变化和城市化等多因素影响，湖泊的增温极大地改变了水生生物的习性、活动规律和代谢强度，从而影响到水生生物的分布和生长繁殖。因此，对湖泊热学机理的正确认识和水温时空变化过程的准确预测是解决上述问题的关键。

深水型湖库的水温分布是众多学者研究的焦点之一。最新研究成果表明，浅水湖泊也可能出现持续数日甚至更久的热力分层现象。湖泊出现明显热分层会造成一系列生态响应，如湖流分层流动、表底层水质浓度差异、种群结构及富营养化过程发生变化。对湖泊热循环机理及热分层规律的准确认识，有助于更好地理解水体的物理、化学和生物过程，为改善湖泊水环境提供技术支撑。

由于水体与床体的热交换通常远小于水气界面的热交换，在模型中通常被忽略。但是，研究表明某些情况下水体与床底的热交换收支也很重要。Tsay证明湖水层化时，水底热交换收支对垂向剖面的水温分布有显著影响。HydroQual在考虑与沉积物床体的热交换情形下模拟了湖的温度层化现象。虽然少数研究考虑了水体与床底的热交换，但往往只考虑床体对水体的作用，实际上床体和水底的热交换是相互的，床体的温度变化同样受水体影响。再者，水底热交换机理非常复杂，除了不同介质温差引起的热量传递外，还存在太阳辐射对床体的直接加热等热现象，特别是对于城市湖泊这类典型的浅水湖，太阳辐射可能穿透整个水体，加热增益甚至可能深达数十米。因此，亟须提出一个实用且足以反映复杂水底热交换机制的水底热通量计算方法。

传统水温模型通常是一维或垂向二维模型，但是横向或纵向平均处理掩盖了流场、水温分布的空间差异性而造成认识水体状态的歧义，垂向上紊动扩散的强度直接决定了温度和速度的垂向梯度，而不同的温度垂向结构又会抑制或促进紊动的发展，由于紊动本身的三维特性，三维模型最能反映水体在真实空间的状态变化，更适合描述复杂的湍流和对流扩散过程。

在研究湖泊热循环规律的基础上，选择合适的水气热交换通量计算方法，提出一种改进的水底热交换通量计算方法。建立三维水动力—水温耦合模型，采用有限差分法对模型方程进行求解。以典型浅水湖泊东湖为研究对象，根据东湖1978年7月的定点观测资料进行水温模拟，探究东湖水温的日变化、月变化过程。

一、数学模型

1.控制方程

为反映陡峭或水深有较大起伏的地形，本模型在垂向上采用σ坐标系，具有一定计算经济性。模型采用A.F.Blumberg和G.L.Mellor提出的三维不可压缩流控制方程，具体方程如下：

式中，U、V、ω 分别为 x、y、σ 三个方向的速度分量；H 为总水深；η为水位；f=2Ωsinφ为科氏力系数，Ω为地球自转角速度，φ为所处纬度；Fx，Fy为水平扩散系数；ρ为水体密度；ρ0为参考密度；KM为垂向扩散系数。

紊流模型为21/2阶的Mellor-Yamada模型，考虑到水深问题，本文引入了浅水修正形式，其控制方程如下：

式中，q2/2为单位质量流体的紊动能；l为流体的紊动混合长；Fq，Fl为方程中的水平扩散项；k为卡门常数，一般取0.4；GH为理查德森数；z0为粗糙高度；E1、E3为常数，分别取1.8,1.0；cs为水体中的声速；为壁近似函数。

水温模型采用对流扩散方程：

太阳辐射的计算选择辐射内穿透模式，即假设到达水体表面的太阳辐射部分被表面吸收，剩余部分穿过表层后被中下层水体吸收，采用Beer–Lambert–Bouguer辐射模型：

式中，T为水温；AHT和AVT为水平温度扩散系数和垂向温度扩散系数；g为重力加速度；cp为水体比热容；TC为源汇项；I为水深σ处的太阳辐射；Is为到达水面处（σ=1）的太阳辐射；βf和 βs为快衰减系数（m-1）和慢衰减系数（m-1）；r为权重系数。

通过联立水动力模型方程组式（1）～（10）和水温模型方程式（11），得到三维水动力—水温耦合模型控制方程组。

2.水面边界条件

水面热交换分为长波辐射（HL）、蒸发潜热（HE）和水汽热传导显热（HC）。 基于 Rosati和Miyakoda（1998年）提出的方法，Hamrick（1992）提出以下模型：

式中，Ts为水面温度；Ta为大气温度；σw是Stefan-Boltzman常数，为 5.67×10-8；ε 为水面发射率，取值 0.97；Bc为经验常数，取值0.8；c为云层覆盖率，介于0～1之间；ea为蒸汽压；es为饱和蒸汽压；ce=ch为湍流交换系数， 取 1.1×10-3；ρa为大气密度；L为蒸发潜热系数；w为风速；Pa为大气压；cpa为空气比热容。

3.水底边界条件

CE-QUAL-W2温度子模型描述沙质床体热通量的公式为：

式中，HB为水体—床体界面的热通量；KB为热交换系数；T1为最底层水温；Tb为床体温度。式（17）的缺点是需要连续的床体温度序列资料做输入，但是床体温度难以实时监测且成本极高，因此需要建立数学模型反映床体的热力学过程。一个描述床体热过程的热平衡方程为：式中，Ib为到达床体的辐射；Hb为床体热通量厚度；ρb为沙质床体密度；cpb为水体—沙质床体混合物比热容；chb为无量纲对流热传系数；和分别为底层水体的水平速度分量。此模型考虑了水体对床体温度的影响和太阳辐射加热床体的过程，更适合浅水湖泊。由于水底热交换的相关研究较少，式（19）中大量参数的率定需要考究水底土壤成分和参考详尽的水情资料，成本过高。本文提出一种改进方法，根据控制工程线性增益思想，引入两个变量放大因子MF（Multiplicative Factor）、位移因子 AF（Additive Factor），将式（13）改写成式（20），通过调整AF、MF值使模拟值与观测值吻合。这样处理将工作重心转移到参数优选上，而不必要在资料收集上花费过多精力，具有实际可行性。

4.定解条件

①水面和水底：

式中,Ts、TB分别为水面处和水底的温度，HL、HE、HC、HB分别为长波辐射、蒸发潜热、水汽热传导显热、水底—床体热通量。

②侧边界采用无滑移边界条件，出入流边界根据实测数据指定流量、水位和温度。

③垂向速度 w（0）=w（-1）=0。

④初始条件给定模型起算时水位、流量和水温的实测值。

二、数值解法

模型采用C网格布置变量 （图1），i、j和k分别为x、y和σ方向的网格编号索引；u布置在网格左右界面的中心（i±1/2，j）；v 布置在网格前后界面的中心 （i，j±1/2）；ω 布置在网格前后界面的中心 （i，k±1/2）；T 布置在网格中心（i，j，k）。上标n+1和n分别表示下一个时间步和当前时间步；Δt为时间步长；Δx、Δy和Δσ分别为x、y和σ方向的网格间距。限于篇幅，有关水动力控制方程式（1）～（10）的求解，见参考文献。水温控制方程（式11）采用有限差分法进行离散，采用物理分步法求解。

采用物理分步法对控制方程进行离散，自编程序实现对方程的求解。该方法根据控制方程中物理过程的频率特征进行分裂求解，对流项、水平扩散项等快过程和垂向扩散项等慢过程对于时间步长的要求也不同。因此将温度方程分为两步做时间积分，并定义T1为Tn与Tn+1的中间时间层变量，基本格式框架是时间导数前项差分，空间导数中心差分，垂向扩散项隐式处理。涉及表层和底层的差分参照边界条件。

①第一步

仅考虑非线性项、对流项、水平扩散项、太阳辐射项，所有这些项均作显式处理求得中间变量T1。

图1 变量布置

式中：

②第二步

仅考虑垂向扩散项，作隐式处理：

继续整理，得到：

式中，M、S、P、R是已知变量的函数，分别在各个网格写出式（28），得到一个关于三对角、非对称矩阵方程组，采用TMDA（追赶法）求得 Tn+1。 最后将求得的 un+1、vn+1、ωn+1、Tn+1等变量代入新的时间层，循环计算，逐一求得每个时间层的结果。

三、模型验证

东湖位于湖北省武汉市武昌城区东部，现为我国水域面积最广阔的城中湖之一。东湖总面积为33.9 km2，多年平均水深2.61 m（图2）,属于典型浅水湖泊。采用100 m×100 m的规则网格离散东湖，被划分成3 008个网格。#1、#2、#3是东湖水温测量点或观测站。

1.模拟结果

以典型气温年1978年为背景，输入7月的气象资料，时间步长为5 s。图3为计算结果。图4为选取1978年7月14日一天的计算结果。

2.计算结果与误差分析

图2 东湖概况

图3 7—8月水温模型结果（#1观测点）：观测值 (三角形)；模拟值 (实线)

图4 7月14日水温模型结果（#1观测点）：观测值 (圆圈)；模拟值 (实线)

对计算水温与实测水温之间的吻合程度进行量化比较，相对误差（RE）为 1.2%，相对均方根误差（RRE）为16.1%，证明了模拟值与观测值吻合良好。除7月25日7天左右的实测值普遍大于模拟值，其余时间段基本重现了#1站点处实际水温的长期变化，误差原因可能是计算中使用的辐射数据只记录了陆面日太阳直射辐射和散射辐射，反射辐射及其他辐射项目并未观测，加上观测数据的不连续性及插值误差，导致与实际的太阳辐射略有偏差；再者温度不同于水位等参数，昼夜差异明显，做日均后一定程度上消除了白天和黑夜的温差变化。水温与气象数据分别由不同观测站观测，而水温的观测时间及间隔受实际情况约束，与气象数据在采样时间尺度上极不匹配，水温观测数据的日均化抑制了湖体的实际温差，因此造成局部失真的结果。

由图4可看出，从上午8时起，随着太阳辐射的增强，东湖水温迅速提高，至16—17时达到最大；此后逐渐降低，天黑后迅速下降，至次日6时达到最低值。图4的模拟结果与金伯欣在1978年7月对东湖展开的热力学调查成果相符，证明了本文提出的模型的有效性。

四、总 结

本文建立了基于CE-QUAL-W2底部热交换模型方程和一种改进的床体热平衡方程，求得水体—床体热交换通量研究下的三维水动力—水温耦合模型，并根据东湖1978年7月的定点观测资料进行了水温模拟，其预测结果与实测资料中水温变化趋势具有良好的一致性，能有效反映东湖水温的日变化、月变化过程，证明了本文提出的模型及其数值解法的可靠性。

[1]Condie,S.A.,I.T.Webster.Estimating stratification in shallow water bodies from mean meteorological conditions[J].Journal of Hydraulic Engineering,2001,127(4).

[2]Carey,Cayelan C.,et al.Eco-physiological adaptations that favour freshwater cyanobacteria in a changing climate[J].Water research,2012,46(5).

[3]Tsay.Thermal stratification modeling of lakes with sediment heat flux [J]. Journal Hydraulic Engineerring,1992,118(3).

[4]HydroQual.A Primer for SEDZL-3D[M].HydroQual,Inc.,1995.

[5]任华堂，夏建新，陶亚.阿海水库水温数值预测研究[J].长江流域资源与环境,2010,19(7).