長江口深水航道三維泥沙數值模型相關問題探討

顧峰峰，郭 賀

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

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

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

表1 國外模型在國內主要河口的應用情況表Tab. 1 Application of some models in major estuaries in China

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

1 長江口物質輸運模擬計算格式的守恒性和穩定性問題

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

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

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

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

2 河口區域河床切應力的準確計算問題

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

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

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

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

(1)

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

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

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

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

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

(2)

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

(3)

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

3 攔門沙區域近底高濃度泥沙形成機理研究及數值模擬問題

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

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

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

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

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

Ri=N2/M2

(4)

其中，

(5)

(6)

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

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

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

(7)

物質輸運方程中的垂向紊動擴散系數(Khv)為：

(8)

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

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

圖6 長江口航道沿程垂線紊動系數(基于實測資料和經驗公式推算)和鹽度場分布(注：圖中實線為鹽度等值線，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)

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

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

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

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

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

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

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

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

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

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

上述泥沙在槽內匯聚的多種影響因子示意圖見圖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

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

圖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 結 語

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

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

2) 近底層河床切應力是長江口河床泥沙起懸和落淤的關鍵參數，目前河床切應力估算的不準確往往導致河床沖淤變化和近底泥沙通量的不準確，同時影響垂線紊流及水動力結構模擬，因而在長江口應用時，需要修正河床切應力的估算方法，以解決局部區域動力、近底泥沙通量計算不合理以及多個影響因子之間耦合、循環作用的問題。

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

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

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