湖泊三維水動力—水溫耦合模型及其應(yīng)用研究

康 玲，靖 爭

（華中科技大學(xué)水電與數(shù)字化工程學(xué)院，430074，武漢）

湖泊是重要的城市水體形態(tài)，具有調(diào)蓄雨水、維持生物多樣性、美化環(huán)境等功能。水溫是影響水生動植物生長、繁衍和遷徙的重要因素之一。然而，隨著氣候變化和城市化等多因素影響，湖泊的增溫極大地改變了水生生物的習(xí)性、活動規(guī)律和代謝強度，從而影響到水生生物的分布和生長繁殖。因此，對湖泊熱學(xué)機理的正確認(rèn)識和水溫時空變化過程的準(zhǔn)確預(yù)測是解決上述問題的關(guān)鍵。

深水型湖庫的水溫分布是眾多學(xué)者研究的焦點之一。最新研究成果表明，淺水湖泊也可能出現(xiàn)持續(xù)數(shù)日甚至更久的熱力分層現(xiàn)象。湖泊出現(xiàn)明顯熱分層會造成一系列生態(tài)響應(yīng)，如湖流分層流動、表底層水質(zhì)濃度差異、種群結(jié)構(gòu)及富營養(yǎng)化過程發(fā)生變化。對湖泊熱循環(huán)機理及熱分層規(guī)律的準(zhǔn)確認(rèn)識，有助于更好地理解水體的物理、化學(xué)和生物過程，為改善湖泊水環(huán)境提供技術(shù)支撐。

由于水體與床體的熱交換通常遠小于水氣界面的熱交換，在模型中通常被忽略。但是，研究表明某些情況下水體與床底的熱交換收支也很重要。Tsay證明湖水層化時，水底熱交換收支對垂向剖面的水溫分布有顯著影響。HydroQual在考慮與沉積物床體的熱交換情形下模擬了湖的溫度層化現(xiàn)象。雖然少數(shù)研究考慮了水體與床底的熱交換，但往往只考慮床體對水體的作用，實際上床體和水底的熱交換是相互的，床體的溫度變化同樣受水體影響。再者，水底熱交換機理非常復(fù)雜，除了不同介質(zhì)溫差引起的熱量傳遞外，還存在太陽輻射對床體的直接加熱等熱現(xiàn)象，特別是對于城市湖泊這類典型的淺水湖，太陽輻射可能穿透整個水體，加熱增益甚至可能深達數(shù)十米。因此，亟須提出一個實用且足以反映復(fù)雜水底熱交換機制的水底熱通量計算方法。

傳統(tǒng)水溫模型通常是一維或垂向二維模型，但是橫向或縱向平均處理掩蓋了流場、水溫分布的空間差異性而造成認(rèn)識水體狀態(tài)的歧義，垂向上紊動擴散的強度直接決定了溫度和速度的垂向梯度，而不同的溫度垂向結(jié)構(gòu)又會抑制或促進紊動的發(fā)展，由于紊動本身的三維特性，三維模型最能反映水體在真實空間的狀態(tài)變化，更適合描述復(fù)雜的湍流和對流擴散過程。

在研究湖泊熱循環(huán)規(guī)律的基礎(chǔ)上，選擇合適的水氣熱交換通量計算方法，提出一種改進的水底熱交換通量計算方法。建立三維水動力—水溫耦合模型，采用有限差分法對模型方程進行求解。以典型淺水湖泊東湖為研究對象，根據(jù)東湖1978年7月的定點觀測資料進行水溫模擬，探究東湖水溫的日變化、月變化過程。

一、數(shù)學(xué)模型

1.控制方程

為反映陡峭或水深有較大起伏的地形，本模型在垂向上采用σ坐標(biāo)系，具有一定計算經(jīng)濟性。模型采用A.F.Blumberg和G.L.Mellor提出的三維不可壓縮流控制方程，具體方程如下：

式中，U、V、ω 分別為 x、y、σ 三個方向的速度分量；H 為總水深；η為水位；f=2Ωsinφ為科氏力系數(shù)，Ω為地球自轉(zhuǎn)角速度，φ為所處緯度；Fx，F(xiàn)y為水平擴散系數(shù)；ρ為水體密度；ρ0為參考密度；KM為垂向擴散系數(shù)。

紊流模型為21/2階的Mellor-Yamada模型，考慮到水深問題，本文引入了淺水修正形式，其控制方程如下：

式中，q2/2為單位質(zhì)量流體的紊動能；l為流體的紊動混合長；Fq，F(xiàn)l為方程中的水平擴散項；k為卡門常數(shù)，一般取0.4；GH為理查德森數(shù)；z0為粗糙高度；E1、E3為常數(shù)，分別取1.8,1.0；cs為水體中的聲速；為壁近似函數(shù)。

水溫模型采用對流擴散方程：

太陽輻射的計算選擇輻射內(nèi)穿透模式，即假設(shè)到達水體表面的太陽輻射部分被表面吸收，剩余部分穿過表層后被中下層水體吸收，采用Beer–Lambert–Bouguer輻射模型：

式中，T為水溫；AHT和AVT為水平溫度擴散系數(shù)和垂向溫度擴散系數(shù)；g為重力加速度；cp為水體比熱容；TC為源匯項；I為水深σ處的太陽輻射；Is為到達水面處（σ=1）的太陽輻射；βf和 βs為快衰減系數(shù)（m-1）和慢衰減系數(shù)（m-1）；r為權(quán)重系數(shù)。

通過聯(lián)立水動力模型方程組式（1）～（10）和水溫模型方程式（11），得到三維水動力—水溫耦合模型控制方程組。

2.水面邊界條件

水面熱交換分為長波輻射（HL）、蒸發(fā)潛熱（HE）和水汽熱傳導(dǎo)顯熱（HC）。 基于 Rosati和Miyakoda（1998年）提出的方法，Hamrick（1992）提出以下模型：

式中，Ts為水面溫度；Ta為大氣溫度；σw是Stefan-Boltzman常數(shù)，為 5.67×10-8；ε 為水面發(fā)射率，取值 0.97；Bc為經(jīng)驗常數(shù)，取值0.8；c為云層覆蓋率，介于0～1之間；ea為蒸汽壓；es為飽和蒸汽壓；ce=ch為湍流交換系數(shù)， 取 1.1×10-3；ρa為大氣密度；L為蒸發(fā)潛熱系數(shù)；w為風(fēng)速；Pa為大氣壓；cpa為空氣比熱容。

3.水底邊界條件

CE-QUAL-W2溫度子模型描述沙質(zhì)床體熱通量的公式為：

式中，HB為水體—床體界面的熱通量；KB為熱交換系數(shù)；T1為最底層水溫；Tb為床體溫度。式（17）的缺點是需要連續(xù)的床體溫度序列資料做輸入，但是床體溫度難以實時監(jiān)測且成本極高，因此需要建立數(shù)學(xué)模型反映床體的熱力學(xué)過程。一個描述床體熱過程的熱平衡方程為：式中，Ib為到達床體的輻射；Hb為床體熱通量厚度；ρb為沙質(zhì)床體密度；cpb為水體—沙質(zhì)床體混合物比熱容；chb為無量綱對流熱傳系數(shù)；和分別為底層水體的水平速度分量。此模型考慮了水體對床體溫度的影響和太陽輻射加熱床體的過程，更適合淺水湖泊。由于水底熱交換的相關(guān)研究較少，式（19）中大量參數(shù)的率定需要考究水底土壤成分和參考詳盡的水情資料，成本過高。本文提出一種改進方法，根據(jù)控制工程線性增益思想，引入兩個變量放大因子MF（Multiplicative Factor）、位移因子 AF（Additive Factor），將式（13）改寫成式（20），通過調(diào)整AF、MF值使模擬值與觀測值吻合。這樣處理將工作重心轉(zhuǎn)移到參數(shù)優(yōu)選上，而不必要在資料收集上花費過多精力，具有實際可行性。

4.定解條件

①水面和水底：

式中,Ts、TB分別為水面處和水底的溫度，HL、HE、HC、HB分別為長波輻射、蒸發(fā)潛熱、水汽熱傳導(dǎo)顯熱、水底—床體熱通量。

②側(cè)邊界采用無滑移邊界條件，出入流邊界根據(jù)實測數(shù)據(jù)指定流量、水位和溫度。

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

④初始條件給定模型起算時水位、流量和水溫的實測值。

二、數(shù)值解法

模型采用C網(wǎng)格布置變量 （圖1），i、j和k分別為x、y和σ方向的網(wǎng)格編號索引；u布置在網(wǎng)格左右界面的中心（i±1/2，j）；v 布置在網(wǎng)格前后界面的中心 （i，j±1/2）；ω 布置在網(wǎng)格前后界面的中心 （i，k±1/2）；T 布置在網(wǎng)格中心（i，j，k）。上標(biāo)n+1和n分別表示下一個時間步和當(dāng)前時間步；Δt為時間步長；Δx、Δy和Δσ分別為x、y和σ方向的網(wǎng)格間距。限于篇幅，有關(guān)水動力控制方程式（1）～（10）的求解，見參考文獻。水溫控制方程（式11）采用有限差分法進行離散，采用物理分步法求解。

采用物理分步法對控制方程進行離散，自編程序?qū)崿F(xiàn)對方程的求解。該方法根據(jù)控制方程中物理過程的頻率特征進行分裂求解，對流項、水平擴散項等快過程和垂向擴散項等慢過程對于時間步長的要求也不同。因此將溫度方程分為兩步做時間積分，并定義T1為Tn與Tn+1的中間時間層變量，基本格式框架是時間導(dǎo)數(shù)前項差分，空間導(dǎo)數(shù)中心差分，垂向擴散項隱式處理。涉及表層和底層的差分參照邊界條件。

①第一步

僅考慮非線性項、對流項、水平擴散項、太陽輻射項，所有這些項均作顯式處理求得中間變量T1。

圖1 變量布置

式中：

②第二步

僅考慮垂向擴散項，作隱式處理：

繼續(xù)整理，得到：

式中，M、S、P、R是已知變量的函數(shù)，分別在各個網(wǎng)格寫出式（28），得到一個關(guān)于三對角、非對稱矩陣方程組，采用TMDA（追趕法）求得 Tn+1。 最后將求得的 un+1、vn+1、ωn+1、Tn+1等變量代入新的時間層，循環(huán)計算，逐一求得每個時間層的結(jié)果。

三、模型驗證

東湖位于湖北省武漢市武昌城區(qū)東部，現(xiàn)為我國水域面積最廣闊的城中湖之一。東湖總面積為33.9 km2，多年平均水深2.61 m（圖2）,屬于典型淺水湖泊。采用100 m×100 m的規(guī)則網(wǎng)格離散東湖，被劃分成3 008個網(wǎng)格。#1、#2、#3是東湖水溫測量點或觀測站。

1.模擬結(jié)果

以典型氣溫年1978年為背景，輸入7月的氣象資料，時間步長為5 s。圖3為計算結(jié)果。圖4為選取1978年7月14日一天的計算結(jié)果。

2.計算結(jié)果與誤差分析

圖2 東湖概況

圖3 7—8月水溫模型結(jié)果（#1觀測點）：觀測值 (三角形)；模擬值 (實線)

圖4 7月14日水溫模型結(jié)果（#1觀測點）：觀測值 (圓圈)；模擬值 (實線)

對計算水溫與實測水溫之間的吻合程度進行量化比較，相對誤差（RE）為 1.2%，相對均方根誤差（RRE）為16.1%，證明了模擬值與觀測值吻合良好。除7月25日7天左右的實測值普遍大于模擬值，其余時間段基本重現(xiàn)了#1站點處實際水溫的長期變化，誤差原因可能是計算中使用的輻射數(shù)據(jù)只記錄了陸面日太陽直射輻射和散射輻射，反射輻射及其他輻射項目并未觀測，加上觀測數(shù)據(jù)的不連續(xù)性及插值誤差，導(dǎo)致與實際的太陽輻射略有偏差；再者溫度不同于水位等參數(shù)，晝夜差異明顯，做日均后一定程度上消除了白天和黑夜的溫差變化。水溫與氣象數(shù)據(jù)分別由不同觀測站觀測，而水溫的觀測時間及間隔受實際情況約束，與氣象數(shù)據(jù)在采樣時間尺度上極不匹配，水溫觀測數(shù)據(jù)的日均化抑制了湖體的實際溫差，因此造成局部失真的結(jié)果。

由圖4可看出，從上午8時起，隨著太陽輻射的增強，東湖水溫迅速提高，至16—17時達到最大；此后逐漸降低，天黑后迅速下降，至次日6時達到最低值。圖4的模擬結(jié)果與金伯欣在1978年7月對東湖展開的熱力學(xué)調(diào)查成果相符，證明了本文提出的模型的有效性。

四、總 結(jié)

本文建立了基于CE-QUAL-W2底部熱交換模型方程和一種改進的床體熱平衡方程，求得水體—床體熱交換通量研究下的三維水動力—水溫耦合模型，并根據(jù)東湖1978年7月的定點觀測資料進行了水溫模擬，其預(yù)測結(jié)果與實測資料中水溫變化趨勢具有良好的一致性，能有效反映東湖水溫的日變化、月變化過程，證明了本文提出的模型及其數(shù)值解法的可靠性。

[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]任華堂，夏建新，陶亞.阿海水庫水溫數(shù)值預(yù)測研究[J].長江流域資源與環(huán)境,2010,19(7).