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

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%。通過與驗證效果的比較不難看出,這種差別主要由連接斷面處的