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

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

1幾種類(lèi)型的耦合計(jì)算方法的比較近年來(lái)，珠江三角洲的水環(huán)境和咸潮入侵日益突出?？紤]到河流網(wǎng)絡(luò)和河流入口區(qū)域的原因和影響，我們無(wú)法全面分析它們的原因和影響，也無(wú)法尋求基本的管理對(duì)策。因此，有必要發(fā)展珠江三角洲的全球水環(huán)境數(shù)學(xué)模型。針對(duì)此問(wèn)題,國(guó)內(nèi)不少學(xué)者將河網(wǎng)區(qū)、河口區(qū)以及近岸海域結(jié)合起來(lái)視作一個(gè)整體,開(kāi)展了一維模型與多維模型耦合計(jì)算的研究,并已取得一定的進(jìn)展,但其主要成果集中在一、二維模型方面,一維模型與三維模型連接計(jì)算的研究則較為少見(jiàn)。其中,彭靜等、徐峰俊等的研究基本實(shí)現(xiàn)了珠江口一、二維模型的耦合計(jì)算,然而,珠江口由鹽度梯度所產(chǎn)生的密度斜壓效應(yīng)明顯,尤其是洪季時(shí)鹽度層化現(xiàn)象異常顯著,僅采用二維模型來(lái)描述其水動(dòng)力特征遠(yuǎn)遠(yuǎn)不夠,無(wú)法體現(xiàn)河口中鹽度分層結(jié)構(gòu)、密度環(huán)流及鋒面結(jié)構(gòu)等對(duì)物質(zhì)輸送的影響;包蕓等的工作雖然實(shí)現(xiàn)了一維與三維斜壓模型的耦合計(jì)算,但實(shí)際上并未涉及鹽度的連接計(jì)算,無(wú)法動(dòng)態(tài)反映鹽水在河網(wǎng)區(qū)與河口區(qū)之間的回蕩并沿河道上溯的過(guò)程;逄勇等開(kāi)展了一、三維水質(zhì)模型的耦合計(jì)算研究,但其物質(zhì)輸送連接條件簡(jiǎn)單,存在當(dāng)一維向三維流動(dòng)時(shí)一維模型無(wú)法滿足三維模型分層邊界條件的問(wèn)題,這實(shí)質(zhì)上是由于三維模型考慮分層結(jié)構(gòu)而一維模型不考慮分層結(jié)構(gòu),兩個(gè)模型信息傳遞不對(duì)稱(chēng)所產(chǎn)生的問(wèn)題。從研究枯季咸潮入侵的角度來(lái)看,需將一維鹽度模型與三維斜壓模型相耦合,方能模擬鹽水的入侵過(guò)程;從水環(huán)境管理的角度來(lái)看,構(gòu)建整體斜壓模型是建立珠江三角洲整體水環(huán)境數(shù)學(xué)模型的前提與基礎(chǔ)。因此,本文在前人的基礎(chǔ)上將一維鹽度模型與三維斜壓模型進(jìn)行耦合計(jì)算,在鹽度連接過(guò)程中,嘗試通過(guò)延長(zhǎng)三維模型上邊界河段的方法來(lái)獲得連接斷面處鹽度的三維分布信息,以此補(bǔ)充三維模型的分層邊界條件,更為合理地構(gòu)建珠江三角洲一、三維耦合斜壓模型,從而為咸潮入侵、污染物通量及其環(huán)境效應(yīng)等問(wèn)題的研究奠定基礎(chǔ)。由于珠江口夏季期間受強(qiáng)烈的太陽(yáng)輻射影響,溫度層化作用并不明顯,鹽度梯度比溫度梯度對(duì)密度變化的作用更大,因此本文只考慮鹽度梯度對(duì)斜壓效應(yīng)的貢獻(xiàn)。2維耦合模型模擬范圍涵蓋整個(gè)珠江三角洲感潮河網(wǎng)區(qū)以及珠江口海域,如圖1所示。河網(wǎng)區(qū)采用一維模型進(jìn)行模擬,共概化河道299條,汊點(diǎn)189個(gè),劃分河道斷面1726個(gè),上游設(shè)置5個(gè)流量控制邊界,分別為西江的高要、北江的石角、東江的博羅、流溪河的老鴉崗和潭江的石咀;八大入?？陂T(mén)虎門(mén)、蕉門(mén)、洪奇瀝、橫門(mén)、磨刀門(mén)、雞啼門(mén)、虎跳門(mén)、崖門(mén)作為一維與三維模型計(jì)算的連接斷面,即為耦合模型的內(nèi)斷面;河口區(qū)采用三維模型進(jìn)行模擬,下邊界延伸至外海,其中西邊界取至下川島西部的鎮(zhèn)海灣,東邊界取至紅海灣,南邊界取至70m等深線。三維模型采用能與物理邊界吻合良好且分辨率較高的正交曲線網(wǎng)格系統(tǒng),水平方向劃分網(wǎng)格為183×186,垂直方向采用sigma坐標(biāo),均勻分為6層,網(wǎng)格系統(tǒng)如圖1所示。三維模型外海開(kāi)邊界采用由美國(guó)NASA噴氣推進(jìn)實(shí)驗(yàn)室提供的TPMGDRB數(shù)據(jù)集高度計(jì)觀測(cè)資料,以M2、S2、K1、O1這4個(gè)主要天文分潮進(jìn)行驅(qū)動(dòng),以解決外海邊界條件因外海潮位資料的缺乏而難以確定的困難。3數(shù)值模擬方法3.11河網(wǎng)水動(dòng)力與鹽酸鹽的三維模型3.1.1x+zt,t明渠中的非恒定水流運(yùn)動(dòng)用一維圣維南方程組來(lái)描述,基本方程形式如下。連續(xù)性方程1B?Q?x+?Ζ?t=qL(1)1B?Q?x+?Z?t=qL(1)動(dòng)力學(xué)方程?u?t+u?u?x+g?Ζ?x+gu|u|C2sR=0(2)式中:Z為斷面水位;B為水面寬度;Q為流量;qL為旁側(cè)入流流量;x為沿程距離;t為時(shí)間;u為平均流速,u=Q/A;A為過(guò)水面積;g為重力加速度;R為水力半徑;Cs為謝才系數(shù)。對(duì)控制方程組采用四點(diǎn)偏心Preissmann隱式差分格式進(jìn)行離散,采用三級(jí)聯(lián)合解法求解。3.1.2河道斷面?zhèn)€數(shù)河道方程?(AS)?t+?(QS)?x-??x(AEx?S?x)-Se=0(3)河道汊點(diǎn)方程ΝL∑l=1(QS)l,j=(SΩ)j(dΖdt)j(4)式中:Ex為縱向分散系數(shù);S為鹽度;Se為外部源匯項(xiàng);Ω為汊點(diǎn)水面面積;j為汊點(diǎn)編號(hào);l為與汊點(diǎn)j相聯(lián)接的河道編號(hào);NL為汊點(diǎn)連接斷面?zhèn)€數(shù)。一維河網(wǎng)鹽度模型采用隱式差分迎風(fēng)格式離散,方程的求解方法見(jiàn)文獻(xiàn)。3.2靜力彈性方程連續(xù)方程?u?x+?v?y+?w?z=0(5)x方向的運(yùn)動(dòng)方程?u?t+u?u?x+v?u?y+w?u?z-fv=-1ρ?p?x+AΜΔ2u+??z(ΚΜ?u?z)+Fx(6)y方向的運(yùn)動(dòng)方程?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分別為水平和垂向渦動(dòng)黏性系數(shù);f為科氏力系數(shù);P為壓力;ρ為海水密度;g為重力加速度;Fx、Fy分別為x、y方向上的外力;S為鹽度;AH、KH分別為水平和垂向渦動(dòng)擴(kuò)散系數(shù);Δ2為水平拉普拉斯算子?？刂品匠痰木唧w差分形式及求解步驟參見(jiàn)文獻(xiàn)。41.3三維結(jié)合斜率模型的構(gòu)建4.1基于連接條件和耦合計(jì)算方法的連接條件4.1.1模型求解與顯式耦合一維與三維模型連接斷面為八大入海口門(mén)。設(shè)一維模型在第i個(gè)口門(mén)的連接斷面上的水位、流量、鹽度分別為Z(1)i、Q(1)i、S(1)i,三維模型在第i個(gè)口門(mén)的連接斷面上平均水位、流量、平均鹽度分別為Z(3)i、Q(3)i、S(3)i,(i=1,2,3,4,5,6,7,8),該斷面位于第j列、第k行、第l層的網(wǎng)格點(diǎn)的鹽度為S(3)ijkl?？紤]聯(lián)解時(shí),一維和三維的連接斷面被視為內(nèi)斷面,因此每個(gè)連接斷面需補(bǔ)充水位連續(xù)、流量連續(xù)、鹽度輸送連續(xù)三個(gè)關(guān)系式,聯(lián)解的思路如圖2所示,采用顯式耦合,三維模型將水位傳遞給一維模型,一維模型將流量傳遞給三維模型,鹽度傳遞的物質(zhì)通量包括由移流作用及擴(kuò)散作用所產(chǎn)生的凈通量,其中移流通量傳遞取決于連接斷面水流方向,漲潮時(shí)由三維模型向一維模型傳遞,落潮時(shí)由一維模型向三維模型傳遞,擴(kuò)散通量傳遞取決于連接斷面處的濃度梯度。由于三維模型在連接斷面處含有多個(gè)計(jì)算水點(diǎn),因而需獲悉連接斷面處鹽度沿寬度和水深方向的變化,方能補(bǔ)充其分層邊界條件。連接的基本條件如下。水位連續(xù)Ζ(1)i=Ζ(3)i(10)流量連續(xù)Q(1)i+Q(3)i=0(11)鹽度輸送連續(xù):(1)當(dāng)水流方向?yàn)閺娜S到一維時(shí)S(1)i=(Q(3)i?S(3)i+FDΙFi)/Q(1)i(12)(2)當(dāng)水流方向?yàn)閺囊痪S到三維時(shí)S(3)ijkl=(Q(1)i?S(1)i+FDΙFi)/Q(3)i?aijkl(13)式中:Q(1)i、Q(3)i以流向連接斷面為負(fù);FDIFi為由擴(kuò)散作用而產(chǎn)生的通量;aijkl為第i個(gè)口門(mén)的連接斷面上鹽度沿寬度和水深方向的分布系數(shù)。當(dāng)一維模型斷面濃度高于三維模型斷面濃度時(shí),擴(kuò)散通量由一維模型向三維模型傳遞,此時(shí)FDΙFi=Ai?Ex?S(1)i-S(3)iΔn,反之,FDΙFi=Ai?AΗ?S(3)i-S(1)iΔn,其中Ai為第i個(gè)口門(mén)的過(guò)水面積,n為法向。4.1.2計(jì)算流水管理體制具體耦合計(jì)算過(guò)程如下:在第n個(gè)時(shí)層首先進(jìn)行三維模型斜壓計(jì)算,將三維模型計(jì)算得到的水位賦給一維模型,補(bǔ)充一維模型下邊界條件;然后進(jìn)行一維水動(dòng)力與鹽度模擬計(jì)算,將一維模型計(jì)算得到的流量賦給三維模型,作為三維模型的上邊界條件;最后根據(jù)水流方向判斷鹽度傳遞方向,補(bǔ)充一維或者三維模型的鹽度邊界條件,進(jìn)入第n+1個(gè)時(shí)層的計(jì)算;如此循環(huán)重復(fù),直至計(jì)算結(jié)束。4.1.3計(jì)算網(wǎng)格點(diǎn)分布系數(shù)由于三維模型描述的是物質(zhì)在空間三維中的運(yùn)動(dòng)和分層現(xiàn)象,在模型計(jì)算中要求提供邊界上物質(zhì)濃度的三維分布信息,而一維模型只能獲得斷面平均濃度,無(wú)法滿足三維模型的需求,因此如何給出三維模型邊界上物質(zhì)的三維分布信息是一、三維連接模型的難題。本文就這一問(wèn)題進(jìn)行嘗試性研究,具體分為3個(gè)步驟:(1)耦合率定后的一維模型與三維模型,進(jìn)行診斷模式計(jì)算,即不考慮鹽度的變化,獲得各連接斷面的流量過(guò)程;(2)將原三維模型計(jì)算網(wǎng)格的上邊界向上游適當(dāng)延長(zhǎng)(延長(zhǎng)后的水平方向網(wǎng)格變?yōu)?83×208,見(jiàn)圖3,垂直方向分層保持不變),使連接斷面變?yōu)樾戮W(wǎng)格的內(nèi)部斷面,以第一步中診斷模式計(jì)算得到的連接斷面處流量作為延長(zhǎng)網(wǎng)格后的三維模型的上邊界條件,單獨(dú)運(yùn)行三維模型斜壓模式,獲取該流量過(guò)程控制下連接斷面上各網(wǎng)格點(diǎn)的鹽度,再分別除以連接斷面的平均鹽度,即可得到計(jì)算時(shí)段內(nèi)各網(wǎng)格點(diǎn)的比例系數(shù)aijkl,作為分層的依據(jù);(3)返回原一、三維耦合模型,進(jìn)行斜壓計(jì)算,鹽度根據(jù)關(guān)系式(12)、式(13)進(jìn)行傳遞與分配,即當(dāng)水流方向?yàn)閺囊痪S到三維時(shí),連接斷面上各網(wǎng)格點(diǎn)按相應(yīng)時(shí)刻下的分布系數(shù)aijkl進(jìn)行分配。這一方法的依據(jù)是設(shè)想連接斷面處的流量過(guò)程提前發(fā)生于上游河流,借此獲知在該流量過(guò)程控制與外海潮汐作用雙重影響下連接斷面處鹽度的分層情況。4.2u3000關(guān)于連接計(jì)算的驗(yàn)證主要從以下兩個(gè)角度來(lái)考察一、三維耦合斜壓模型連接計(jì)算的效果,一是比較模擬結(jié)果與實(shí)測(cè)值的擬合程度,即驗(yàn)證效果;二是分析連接計(jì)算是否顯著改變模型的計(jì)算流場(chǎng)與鹽度分布。從直觀的物理過(guò)程考慮,無(wú)論是將兩個(gè)模型單獨(dú)計(jì)算還是連接計(jì)算,模型計(jì)算的范圍和外邊界條件都始終一致,假如連接計(jì)算保證了連接斷面水力要素連續(xù)并趨于實(shí)際的變化過(guò)程,那么,從理論上講,整個(gè)計(jì)算域內(nèi)任意一點(diǎn)連接計(jì)算前后的計(jì)算結(jié)果都應(yīng)該趨于一致。所以,只要滿足了上述兩個(gè)方面的驗(yàn)證,即可以認(rèn)為連接計(jì)算是成功的。模型采用2001枯季的典型水文條件進(jìn)行率定,采用1999年洪季的典型水文條件進(jìn)行驗(yàn)證,具體時(shí)間為1999年7月16日～24日(共205h)、2001年2月7日～16日(共213h)。洪季時(shí)外海邊界的鹽度設(shè)為表層33.5PSU、底層34.5PSU,枯季時(shí)均為34PSU。洪季、枯季以各自預(yù)熱運(yùn)行60d后的結(jié)果作為正式計(jì)算的初始條件,時(shí)間步長(zhǎng)為40s。糙率是影響一維水動(dòng)力模型計(jì)算精度的主要參數(shù),河網(wǎng)區(qū)不同流域河段的糙率分別為:北江流域三水以上取值0.026～0.045,三水以下糙率為0.015～0.030;西江流域糙率為0.020～0.035;口門(mén)段糙率為0.010～0.020。一維鹽度模型的縱向分散系數(shù)取值5.0m2/s。三維模型底部糙率系數(shù)Z0取值0.02cm。在連接斷面縱剖面上鹽度分布系數(shù)aijkl的計(jì)算中,洪季時(shí)上游徑流量大,口門(mén)附近基本上被淡水所控制,模型計(jì)算得到的鹽度分布系數(shù)aijkl為0,而枯季時(shí)上游流量微弱,潮汐作用相對(duì)顯著,因而鹽水向河口上游入侵距離較長(zhǎng),鹽度分布系數(shù)aijkl有一定變化,具體由模型計(jì)算得到。4.2.1洪季枯季鹽度計(jì)算誤差變化水動(dòng)力計(jì)算結(jié)果驗(yàn)證采用河網(wǎng)區(qū)內(nèi)的50多個(gè)水文站、八大口門(mén)水文站、珠江口海區(qū)數(shù)個(gè)潮位站的同步實(shí)測(cè)資料,圖5、圖6為虎門(mén)、磨刀門(mén)的驗(yàn)證結(jié)果。分析結(jié)果顯示,洪季、枯季的水位、流量的相位偏差均在1h以?xún)?nèi),洪季時(shí)80%河網(wǎng)區(qū)站點(diǎn)的水位振幅偏差、流量相對(duì)誤差分別在20cm、20%以?xún)?nèi),枯季時(shí)所有站點(diǎn)的水位振幅偏差均在20cm以?xún)?nèi);在連接斷面處,無(wú)論是洪季還是枯季,八大口門(mén)的水位計(jì)算誤差均在12cm以?xún)?nèi),總體平均誤差為5cm;洪季時(shí)八大口門(mén)的流量相對(duì)誤差均在12%以?xún)?nèi);表1列出了1999年洪季八大口門(mén)分流比計(jì)算與實(shí)測(cè)的對(duì)比結(jié)果,可以看出,計(jì)算結(jié)果與實(shí)測(cè)的基本一致,絕對(duì)誤差均在2%以?xún)?nèi),東四口門(mén)與西四口門(mén)分流比的絕對(duì)誤差更是只有0.07%。另外,受資料所限,洪季的鹽度計(jì)算結(jié)果只采用位于內(nèi)伶仃島附近的C1站(位置見(jiàn)圖4)實(shí)測(cè)資料進(jìn)行驗(yàn)證。由圖7可以看出,模型模擬的鹽度變化趨勢(shì)與實(shí)際過(guò)程較為吻合,表層、中層、底層計(jì)算偏差依次為0.2、1.7、2.6;底層鹽度計(jì)算結(jié)果較實(shí)測(cè)的偏低,可能主要是由于模型所采用的垂向分辨率較低,垂向過(guò)于均勻化所致。枯季時(shí)潮汐作用明顯,表2列出八大口門(mén)鹽度計(jì)算與實(shí)測(cè)的對(duì)比結(jié)果,計(jì)算結(jié)果與實(shí)測(cè)的基本一致,平均偏差為1.3。圖8反映了蕉門(mén)、麻涌站受潮汐影響的鹽度動(dòng)態(tài)變化過(guò)程,模型所刻畫(huà)的變化過(guò)程與實(shí)測(cè)基本一致。綜合上述驗(yàn)證結(jié)果,表明本文所建立的耦合斜壓模型能再現(xiàn)水流、鹽度的實(shí)際變化過(guò)程,模擬的水位、流量、鹽度與實(shí)測(cè)值吻合較好,偏差均控制在合理范圍以?xún)?nèi)。4.2.2連接計(jì)算與單獨(dú)計(jì)算的比較為進(jìn)一步考察連接計(jì)算的效果,以1999年洪季實(shí)測(cè)流量過(guò)程作為八大口門(mén)的邊界條件單獨(dú)運(yùn)行三維斜壓模型,將其計(jì)算結(jié)果與耦合斜壓模型計(jì)算結(jié)果進(jìn)行比較,分析連接計(jì)算對(duì)三維模型計(jì)算流場(chǎng)、鹽度分布的影響程度。由于單獨(dú)計(jì)算時(shí)八大口門(mén)的流量采用實(shí)測(cè)值,而連接計(jì)算時(shí)八大口門(mén)的流量是通過(guò)計(jì)算得到,與實(shí)測(cè)結(jié)果不可避免會(huì)存在偏差,因此連接計(jì)算的結(jié)果與單獨(dú)計(jì)算的結(jié)果必定存在一定的偏差。但從理論上講,連接計(jì)算與單獨(dú)計(jì)算的偏差越小,說(shuō)明連接計(jì)算越理想。從理論上講,越靠近連接斷面的區(qū)域,受連接計(jì)算影響越顯著,所以將重點(diǎn)討論連接斷面前后一段范圍內(nèi)的計(jì)算結(jié)果。一維部分選擇8個(gè)比較斷面O1～O8,均與各連接斷面相隔一個(gè)斷面,比較其水位、流量;三維部分選擇8個(gè)比較點(diǎn)T1～T8,位置分布如圖4所示,比較其水位、流速、鹽度。受篇幅所限,僅列出部分一維、三維的比較結(jié)果,見(jiàn)圖9～圖11。從水動(dòng)力比較結(jié)果來(lái)看,連接計(jì)算與單獨(dú)計(jì)算的水位與流量的相位差別均很小,兩者一維部分的水位差為5.5cm,流量相對(duì)偏差為5.5%,三維部分的水位差為4.0cm,流速相對(duì)偏差為3.0%。通過(guò)與驗(yàn)證效果的比較不難看出,這種差別主要由連接斷面處的