長江口深水航道三維泥沙數(shù)值模型相關(guān)問題探討

顧峰峰，郭 賀

(1. 上海河口海岸科學(xué)研究中心 河口海岸交通行業(yè)重點實驗室，上海 201201; 2. 上海海事大學(xué) 海洋科學(xué)與工程學(xué)院，上海 201306)

長江口作為長江出?？?，受徑潮流相互作用，形成了獨特的水沙鹽輸運和地貌特征：中等潮汐強(qiáng)度作用下的外海驅(qū)動力，年際變化達(dá)6～7倍的上游徑流輸入，攔門沙區(qū)域強(qiáng)勁的鹽水入侵、上溯和鹽淡水交匯，近底層存在形成機(jī)理復(fù)雜的細(xì)顆粒高濃度泥沙帶。長江口歷時13年、投資上百億打造的我國最為繁忙的黃金水道——長江口12.5 m深水航道，這條航道每年花費10～20億的巨資進(jìn)行維護(hù)，疏浚量曾經(jīng)達(dá)到了1億m3/a，目前穩(wěn)定在5×107m3/a左右。在航槽近底處經(jīng)常能觀測到存在約百公斤的高濃度泥沙層，航道回淤的機(jī)理十分復(fù)雜[1-6]。目前面臨著復(fù)雜的水動力—物質(zhì)(水沙)輸運機(jī)理難以用數(shù)學(xué)公式描述的巨大挑戰(zhàn)，如果有一個三維水動力—物質(zhì)輸運耦合計算的數(shù)值模型，能較為全面的揭示其中復(fù)雜的耦合作用機(jī)理，預(yù)測不同工況和假設(shè)條件下的主要水文及航道回淤特征，對于開展長江口航道回淤等基礎(chǔ)理論研究具有重要意義，一個完善和強(qiáng)壯的數(shù)值模型將是研究此類問題的最理想的試驗和驗證載體。

在過去的二三十年里，河口區(qū)域適用的水動力—物質(zhì)輸運耦合數(shù)值模型發(fā)展取得了巨大的成就。首先，研究的幾何空間從平面二維(2D)轉(zhuǎn)向空間三維(3D)，從平面結(jié)構(gòu)網(wǎng)格轉(zhuǎn)向無結(jié)構(gòu)網(wǎng)格；其次，加入了多CPU、GPU計算和機(jī)群計算等不斷發(fā)展的前沿并行計算技術(shù)；另外，在求解方法上從Casulli提出的預(yù)測步求解等方法[7-10]逐步發(fā)展到了可以突破CFL數(shù)的歐拉拉格朗日離散方法等等[11-18]。上述的這些數(shù)值模型的前沿技術(shù)解決或是部分解決了河口區(qū)域數(shù)值模型計算的“跨尺度”問題[16-18]，即從海洋到近岸的不同尺度帶來的計算效率、穩(wěn)定性、模擬精度、貼合復(fù)雜邊界等問題。目前比較常用的三維數(shù)值模型包括結(jié)構(gòu)網(wǎng)格的POM[19]、TRIM[8]、DELFT3D和ROMS[20]，非結(jié)構(gòu)網(wǎng)格的UnTRIM[21]、ELCIRC，F(xiàn)VCOM[22]、MIKE3、SELFE和SCHISM等；同時國內(nèi)也有一些類似的開發(fā)研究工作[12-15]。上述這些模型都得到了很好的應(yīng)用，國內(nèi)幾個主要河口的模型應(yīng)用情況，如表1[23]。

表1 國外模型在國內(nèi)主要河口的應(yīng)用情況表Tab. 1 Application of some models in major estuaries in China

雖然數(shù)學(xué)模型技術(shù)已經(jīng)取得了巨大的成就，并在近年來開展的長江口水沙鹽、地貌變化特征及其內(nèi)在機(jī)理的深入研究取得了明顯的進(jìn)展[1-2]，但現(xiàn)有的模型技術(shù)在長江口的應(yīng)用時經(jīng)常面臨應(yīng)用困難和鹽度、泥沙輸運模擬精度不高等問題。針對這些問題，很多學(xué)者進(jìn)行了研究，并已取得了不少的研究成果，包括在非結(jié)構(gòu)網(wǎng)格上利用歐拉拉格朗日離散格式、以及并行技術(shù)來克服模型計算效率和復(fù)雜邊界的問題[17-18]。利用有限元方法離散連續(xù)方程，減小網(wǎng)格非正交性計算誤差，以及對物質(zhì)輸運控制方程利用分步歐拉法和動態(tài)時間步長調(diào)整的技術(shù)手段，解決模型計算的“跨尺度”問題和改善物質(zhì)守恒性問題。在長江口區(qū)域的水沙鹽模擬中大量引入鹽度斜壓和紊動制約的作用機(jī)制以提高模擬精度等等[1]。但是，目前這些模型和水沙運動機(jī)理的研究在長江口應(yīng)用中仍然還有不少問題待解決和完善。

1 長江口物質(zhì)輸運模擬計算格式的守恒性和穩(wěn)定性問題

長江口的物質(zhì)輸運包括鹽度、泥沙都是一個長周期的輸運過程，通常至少包括了大中小潮的周期性變化過程(約15天)，因而，對于模型計算物質(zhì)輸運的計算效率及守恒性要求較高，另外攔門沙區(qū)域存在高濃度泥沙層及垂向高度層化特征又要求計算格式兼具極高的穩(wěn)定性，這些要求使得模型的計算步長通常較小，從而使得計算量增大和計算效率降低，通常需要引入并行計算的方法來解決長江口物質(zhì)輸運模擬計算中面臨的守恒性和穩(wěn)定性問題。

另一方面，基于水動力的無條件穩(wěn)定計算格式，研究解決大時間步長條件下的物質(zhì)輸運守恒性和穩(wěn)定性，也是一個重要的研究方向。從目前的技術(shù)條件來看，在無條件穩(wěn)定離散格式下，如對流項采用歐拉拉格朗日離散格式，可以取較大的計算時間步長，使得水動力計算的效率和穩(wěn)定性明顯提升，但較難滿足物質(zhì)輸運的守恒性，通常需要把長時間步長固定或者動態(tài)分解成若干個小時間步長后，利用守恒性控制方程進(jìn)行計算來滿足守恒性要求[17-18]。

根據(jù)上述，長江口物質(zhì)輸運模擬在計算格式還存在一些問題，但總體上來看，并行化計算的模型開發(fā)和無條件穩(wěn)定的守恒性計算格式，也都已經(jīng)取得了不少的研究進(jìn)展，兩者各自或者相結(jié)合都是研究此類問題的模型發(fā)展的重要方向，其研究成果也將大幅提升解決長江口物質(zhì)輸運模擬這類“跨尺度”問題的模擬能力和計算效率。上述問題的示意圖如圖1。

圖1 長江口物質(zhì)輸運模擬計算格式的守恒性和穩(wěn)定性問題Fig. 1 Conservation and stability of the model scheme for the material transport in the Yangtze estuary

2 河口區(qū)域河床切應(yīng)力的準(zhǔn)確計算問題

近底層河床切應(yīng)力是長江口河床泥沙起懸和落淤的關(guān)鍵參數(shù)，其評估方法一般采用對數(shù)流速分布曲線來推算近底河床摩阻和切應(yīng)力(LP方法)[24-29]。然而在河口區(qū)域潮汐作用下，采用LP方法的河床切應(yīng)力評估值的偏差問題十分明顯，誤差甚至可以達(dá)到100%[30-31]，其受非恒定流和垂線層化影響十分明顯。長江口區(qū)域的河床阻力系數(shù)取值通常只有常規(guī)取值的1/4，但實際上它并不是一個常數(shù)。不準(zhǔn)確的河床切應(yīng)力估算往往導(dǎo)致不準(zhǔn)確的河床沖淤變化和近底泥沙通量，同時影響垂線紊流及水動力結(jié)構(gòu)。

實測資料表明，長江口北槽深水航道內(nèi)強(qiáng)淤積段近底層的局部落急流速甚至可達(dá)2 m/s以上，且落潮明顯優(yōu)勢，但每年這些區(qū)域的航道淤積強(qiáng)度最大可達(dá)約10 m/a。當(dāng)采用常規(guī)LP計算方法計算切應(yīng)力的數(shù)學(xué)模型模擬計算該區(qū)域的河床變化時，往往表現(xiàn)出強(qiáng)烈的沖刷，這充分說明不考慮河口區(qū)域河床切應(yīng)力估算問題的數(shù)值模型，其計算結(jié)果往往會導(dǎo)致巨大的誤差。

圖2為常規(guī)對數(shù)流速分布計算的河床切應(yīng)力和坐底三腳架基于ADV觀測資料推算的河床切應(yīng)力(TKE法)的比較[1, 24]，兩者的平均比值可達(dá)4倍多，這一結(jié)果說明了兩者差異較大。在實際的模型率定和驗證時，長江口的三維阻力系數(shù)一般取值約為0.000 8～0.001 5，二維計算的糙率取值約為0.011，都明顯小于常規(guī)取值。圖2中TKE法的應(yīng)力計算公式[24]：

圖2 長江口實測河床切應(yīng)力及常規(guī)對數(shù)流速分布估算的結(jié)果比較(LP和TKE法估算的切應(yīng)力比較)Fig. 2 Comparison of the bed shear stress calculated by TKE method and by LP method

(1)

分析形成上述差異的原因，主要是由于長江口區(qū)域的垂線流速分布在非恒定流、近底密度層化影響及垂線制紊作用等綜合作用下偏移了正常對數(shù)分布，這一結(jié)論符合大多數(shù)學(xué)者的研究成果[30]。長江口實測流速曲線和基于TKE法的切應(yīng)力計算得到的常規(guī)對數(shù)分布曲線[1]的對比如圖3，其結(jié)果也驗證了基于實測紊動切應(yīng)力推算得到流速曲線相對于常規(guī)對數(shù)分布曲線已經(jīng)發(fā)生了明顯的偏移。

圖3 基于ADV實測切應(yīng)力的長江口實測流速曲線與對數(shù)分布曲線對比Fig. 3 Deviation of measured velocity curve and calculated logarithmic distribution curve based on the shear stress measured by ADV in Yangtze estuary

在潮汐非恒定流及紊動制約作用下，針對河床切應(yīng)力計算方法進(jìn)行修正，可以修正近底泥沙通量的計算結(jié)果，同時也會影響垂向水動力結(jié)構(gòu)和紊動強(qiáng)度的計算，進(jìn)一步再影響河床切應(yīng)力和垂線流速偏移，因而這里存在一個循環(huán)作用的耦合作用過程。

根據(jù)上述可知：水動力—物質(zhì)輸運耦合模型在長江口的應(yīng)用，需要修正長江口區(qū)域的河床切應(yīng)力估算方法，以解決局部區(qū)域動力、近底泥沙通量計算不合理以及多個影響因子之間耦合、循環(huán)作用的問題。

一般來說非恒定流條件下的河床切應(yīng)力與非恒定流的流速過程密切相關(guān)，因而非恒定流的切應(yīng)力修正可以從流速變化過程來入手。為了說明這個問題，根據(jù)包含近底層切應(yīng)力τb的一維非恒定流動量方程，忽略了對流項：

(2)

式中：U為垂向平均流速，η為水面，H為水深，g為重力加速度，t為時間?；谑?2)可得如下式：

(3)

從式(3)可知非恒定流條件下的河床切應(yīng)力除了和水面坡降有關(guān)，還與流速隨時間的變化梯度有關(guān)，因而河口區(qū)域河床切應(yīng)力計算修正需要考慮徑潮流作用下的加減速流的影響，這和很多學(xué)者已有的研究結(jié)論是一致的。

3 攔門沙區(qū)域近底高濃度泥沙形成機(jī)理研究及數(shù)值模擬問題

攔門沙區(qū)域的泥沙問題是長江口三維水沙鹽數(shù)值模型技術(shù)研究中關(guān)注的重點問題，長江口近底高濃度泥沙層的存在是攔門沙區(qū)域的重要特征(圖4)[2]，是導(dǎo)致航道局部超強(qiáng)淤積的重要原因之一，而近底高濃度泥沙層的形成與該區(qū)域的鹽度斜壓力和紊動強(qiáng)度變化密切相關(guān)。其主要作用機(jī)理目前在理論上可簡單描述為：紊動制約導(dǎo)致攔門沙區(qū)域的水、沙、鹽分層加劇，泥沙、鹽度垂線上易于匯聚于近底層，并加強(qiáng)了鹽度斜壓力作用下水沙上溯能力，即在泥沙垂向匯聚近底層和縱向上溯能力增強(qiáng)兩者共同作用下，使得近底層水沙沿航道的縱向向下的輸運能力降低，從而縱向上易于匯聚于滯留區(qū)域。在文獻(xiàn)[1]中，給出了明確的分析結(jié)論：攔門沙區(qū)域的紊動制約與鹽度斜壓力共同作用是縱向上形成近底高濃度泥沙的主要原因。

圖4 長江口實測的攔門沙段近底高濃度泥沙集中分布形態(tài)Fig. 4 The concentrated distribution pattern of high concentration sediment measured in the Yangtze estuary

根據(jù)上述，垂向紊動擴(kuò)散系數(shù)是一個描述近底高濃度泥沙形成機(jī)理的重要指標(biāo)參數(shù)，其計算模型一般采用通用的計算模型，例如FVCOM和SCHISM等中采用了垂向一維GOTM紊流模型計算包[32]，其包括了常用的雙方程模型等。目前不同的紊流模型由于基于不同的假設(shè)，得到的計算結(jié)果有一定的差異，在長江口的適用性需要詳細(xì)比較和論證，也需要進(jìn)行參數(shù)的經(jīng)驗率定。不同紊流模型的數(shù)值試驗的計算結(jié)果的差異性，可以參見如圖5[33]。

圖5 不同紊流模型計算的結(jié)果比較Fig. 5 Comparison of calculation results of different turbulence models

為了進(jìn)一步說明垂向紊動系數(shù)的影響，這里選取零方程模型進(jìn)行實例計算說明，該模型形式簡單，在模擬河口物質(zhì)輸運方面有較好的應(yīng)用；該方程模型中垂向紊動擴(kuò)散系數(shù)和物質(zhì)輸運方程垂向擴(kuò)散系數(shù)僅為Richardson數(shù)的函數(shù)，其中Richardson數(shù)(Ri)定義如下：

Ri=N2/M2

(4)

其中，

(5)

(6)

式中：S和C分別為鹽度和含沙濃度；u，v分別為x和y方在垂線位置z處的流速；ρ為渾水密度；ρ0為清水密度。

Richardson數(shù)計算式中的M2表示水平流速垂向梯度的作用，N2的物理含義為密度分層的影響，即Richardson數(shù)反映了這兩種物理作用的相對強(qiáng)弱；一般認(rèn)為Ri>0.25，即會受到明顯的制紊層化作用影響。

動量方程中的垂向紊動擴(kuò)散系數(shù)(Kmv)的計算式為：

(7)

物質(zhì)輸運方程中的垂向紊動擴(kuò)散系數(shù)(Khv)為：

(8)

式中：v0、vb、Kb均為常數(shù)，根據(jù)Pacanowski和Philander的建議，它們的取值分別為v0=5×10-3、vb=10-4、Kb=10-5m2/s。

根據(jù)上述及式(4)～(8)，紊動擴(kuò)散系數(shù)的變化幅度在紊動制約作用下可以減小幾百倍。由長江口實測固定點的垂向水文資料推算的紊動擴(kuò)散系數(shù)的計算分布參如圖6[1]。根據(jù)圖6可知，長江口深水航道縱向上在中下段紊動擴(kuò)散系數(shù)出現(xiàn)了明顯的大幅度減小。

圖6 長江口航道沿程垂線紊動系數(shù)(基于實測資料和經(jīng)驗公式推算)和鹽度場分布(注：圖中實線為鹽度等值線，2012年8月17日05:00)Fig. 6 The vertical turbulence coefficient(based on measured data and empirical formula)and salinity along the Yangtze estuary channel (the solid line is salinity isoline, 05:00 17/8/2012)

通常來說，紊動制約作用將大大降低攔門沙區(qū)域垂線上水、沙、鹽的擴(kuò)散，這將明顯增加垂向上泥沙、鹽度的層化現(xiàn)象，并加劇沿航道縱向的鹽度斜壓上溯作用，兩者都將導(dǎo)致航道近底層縱向向外海的輸運能力降低，從而促進(jìn)近底泥沙和鹽度更易于匯聚于底部，這又進(jìn)一步導(dǎo)致垂線密度分層加劇和紊動制約加強(qiáng)，從而形成了一個循環(huán)促進(jìn)的過程和耦合作用的機(jī)制，這是形成目前攔門沙河段深水航道近底高濃度泥沙層存在的主要原因之一[1-2, 7]。闡述這一問題的示意如圖7所示。

圖7 近底高濃度泥沙形成機(jī)理示意Fig. 7 Turbulence restriction and formation mechanism of high concentration sediment near the bottom

實際上長江口長期受到風(fēng)浪的影響(長江口固定水文站——牛皮礁站的年平均有效波高約0.7 m)，風(fēng)浪條件下的垂向紊動擴(kuò)散系數(shù)將發(fā)生明顯變化，尤其是航道邊灘水深較小使得泥沙垂向分布受到的紊動系數(shù)變化的影響加劇，進(jìn)而影響邊灘泥沙輸運，主要表現(xiàn)為邊灘泥沙紊動加強(qiáng)，泥沙活動性增強(qiáng)。

根據(jù)上述，在長江口進(jìn)行水沙鹽的模擬，采用和研發(fā)一個能夠反映上述近底高濃度泥沙形成機(jī)制的模擬技術(shù)是十分必要的。

4 長江口近底高濃度泥沙橫向入槽模式問題

在長江口人工開挖的航槽內(nèi)，經(jīng)常能在橫斷面上觀測到明顯的高濃度泥沙大量堆積于主槽的現(xiàn)象(圖8)。圖8中的近底高濃度泥沙可以在航道內(nèi)維持相當(dāng)長的時間，但航道邊灘上幾乎觀測不到。一般認(rèn)為密度流作用下泥沙歸槽輸運的影響是重要的成因之一，這一點推論目前尚需要通過更多的數(shù)模試驗和進(jìn)一步的現(xiàn)場觀測來驗證。這里需要補(bǔ)充指出值得關(guān)注的有兩點：

1) 垂線密度分層導(dǎo)致的紊動制約作用對泥沙橫向輸運的影響

在中等強(qiáng)度潮汐作用下長江口北槽航道內(nèi)的垂線紊動強(qiáng)度理論上是較強(qiáng)的。如果沒有前述分析得到的結(jié)論，即近底層存在很強(qiáng)的制紊作用(紊動系數(shù)減小1～2個量級)導(dǎo)致劇烈的水動力及物質(zhì)密度垂線分層，這樣的泥沙匯聚也會瞬時被破壞，無法持續(xù)。

2) 近底層高泥沙濃度的制約沉降作用對泥沙橫向輸運的影響

當(dāng)近底層泥沙以正常沉速和概率落淤時，理論上分析可知近底層高濃度泥沙形成時間也是較為短暫的，大部分泥沙將在邊灘上落淤；然而當(dāng)匯聚到一定程度時，在形成的高濃度泥沙層內(nèi)由于制約沉降作用，會使得沉降速度迅速降低(根據(jù)室內(nèi)試驗成果可知減小為原來的約1/6[1])，使得更多的近底高濃度泥沙參與到上述由密度差異導(dǎo)致的橫向“匯聚”中來，加劇入槽的泥沙量。這一沉降制約的影響同樣適用于縱向近底泥沙匯聚的機(jī)理分析。

上述泥沙在槽內(nèi)匯聚的多種影響因子示意圖見圖8。

圖8 長江口深水航道近底高濃度泥沙斷面分布實測及橫向入槽模式示意Fig. 8 Section distribution and transverse channel model of the highly concentrated sediment near the bottom in the Yangtze estuary deep water channel

長江口北槽航道內(nèi)的高濃度泥沙在洪季時的實測沿程縱向分布參見圖9。從圖8和圖9及資料分析來看這種槽內(nèi)高濃度泥沙的存在有長期性和普遍性，這種近底泥沙分布特征是該區(qū)域航道局部回淤強(qiáng)度(最大可超過10 m/a)明顯大于非航槽區(qū)域的主要原因之一。因而圖8所示的多種影響因子的作用機(jī)制，也是長江口水沙鹽三維數(shù)值模型及水沙運動機(jī)理研究中需要重點考慮和深化研究的重要問題之一。

圖9 長江口深水航槽近底高濃度泥沙沿程分布(高低頻水深測量)Fig. 9 The high concentration sediment distribution near the bottom along the Yangtze estuary deep water channel measured by high and low frequency sounder

5 結(jié) 語

基于近年來長江口三維水沙鹽數(shù)值模型及水沙運動機(jī)理研究的發(fā)展，系統(tǒng)的梳理了在長江口深水航道數(shù)學(xué)模型計算中實際所面臨的若干個重要問題，主要包括以下四個方面：

1) 長江口三維水沙鹽模型離散求解面臨“跨尺度”問題，即從海洋到近岸的不同尺度帶來的計算效率、穩(wěn)定性、模擬精度、貼合復(fù)雜邊界等問題，針對這一問題，需要在現(xiàn)有研究成果的基礎(chǔ)上，進(jìn)一步提升長江口物質(zhì)輸運模擬計算格式的守恒性和穩(wěn)定性。

2) 近底層河床切應(yīng)力是長江口河床泥沙起懸和落淤的關(guān)鍵參數(shù)，目前河床切應(yīng)力估算的不準(zhǔn)確往往導(dǎo)致河床沖淤變化和近底泥沙通量的不準(zhǔn)確，同時影響垂線紊流及水動力結(jié)構(gòu)模擬，因而在長江口應(yīng)用時，需要修正河床切應(yīng)力的估算方法，以解決局部區(qū)域動力、近底泥沙通量計算不合理以及多個影響因子之間耦合、循環(huán)作用的問題。

3) 長江口近底高濃度泥沙層的存在是攔門沙區(qū)域重要特征，是導(dǎo)致深水航道局部強(qiáng)淤積的主要原因之一，因而長江口的三維水沙鹽數(shù)值模型中必須引入近底層高濃度泥沙形成機(jī)制，闡明這一機(jī)制需要選擇和率定一個長江口適用的紊流模型，并需要完善反映紊動制約、鹽度斜壓力等耦合作用機(jī)理和模擬技術(shù)。

4) 對長江口深水航道內(nèi)的高濃度泥沙來說，一般認(rèn)為密度流作用下泥沙歸槽輸運的影響是重要的成因之一，這一點推論需要通過理論分析、模擬試驗和進(jìn)一步的現(xiàn)場觀測來驗證。這里的泥沙歸槽運動需要考慮垂向密度分層導(dǎo)致的紊動制約作用和近底層高泥沙濃度的制約沉降作用等耦合作用的影響。

當(dāng)然也不僅限于本文所列的問題，在長江口還存在其他一些問題也是值得和亟需進(jìn)一步研究，如風(fēng)浪對于水沙運動、地形沖淤問題、航道內(nèi)的浮泥和疏浚影響下航道淤積等問題，希望這些問題能隨著研究的深入而逐步得到解決，從而使得數(shù)值模型得到更多和更廣泛的應(yīng)用。