臺風浪風暴潮作用下三維泥沙數值模擬

朱志夏，熊偉

(1.中設設計集團股份有限公司，江蘇 南京 210014; 2.浙江貴仁信息科技股份有限公司，浙江 杭州 310051)

臺風引起動力強大的臺風浪和風暴潮流很容易造成近岸淺灘大量泥沙的起動，在短時間內造成航道或港池的強淤，使堤防與防波堤等港口航道與海岸及近海工程建筑物遭到嚴重損壞，同時，風暴增水也會給沿海地區帶來安全隱患，因此，對工程海域臺風浪、風暴潮以及泥沙輸移、港口航道強淤等的數值模擬和預測提出了更高的要求。為了更加準確地模擬臺風浪風暴潮作用下的泥沙輸移、海床演變、港池航道強淤，有必要進行臺風浪風暴潮作用下三維泥沙數值模擬技術的研究。

有關波浪潮流共同作用下三維泥沙數值模擬，已取得了不少成果。朱志夏[1]基于Navier-Stokes方程和質量傳輸方程，應用多重尺度法、小參數方法及變分原理，導出了波浪作用下的淺水環流方程、不恒定非均勻潮流場中隨機波的折繞射聯合的數值模式及波流聯合作用下的懸沙擴散方程。開發并建立了波流共同作用下二、三維嵌套泥沙數學模型，成功應用于天津新港附近海域疏浚棄土處理的工程。梁丙臣[2]基于國際著名的水動力-生態-懸浮泥沙耦合模式COHERENS-SED和第3代波浪模式SWAN，對COHERENS-SED進行了完善與發展，最終建立了波流聯合作用下的三維水動力-懸浮泥沙耦合數學模型，并應用結構化網格二重嵌套的方法，成功應用于黃河三角洲濱海區的潮流和懸沙輸運規律的研究。王紅[3]采用大、中、小三重結構化網格二、三維嵌套模式，結合第3代波浪模式SWAN和三維水動力-多組分泥沙模型EFDC，將模型應用于濰坊港海域實際工程問題的研究。由于采用結構化網格，工程中關注的防波堤附近水動力的模擬結果有些失真。堵盤軍[4]基于第3代波浪模式SWAN和ECOMSED模式的改進和優化，進行了長江口、杭州灣海域三維懸沙數值模擬的初步研究。研究表明波浪在近岸潮灘區域對泥沙啟動的作用不可忽視。王平[5]根據考慮流場效應的橢圓型緩坡方程，結合水動力泥沙模型FVCOM(an unstructured grid, finite-volume coastal ocean model)，構建了近岸大范圍三維波流耦合模型和波流聯合作用下的近岸物質輸運模型，分別應用于大連灣海域和旅順琥珀灣海域，模擬灣內的潮流和物質輸運。

綜上所述，臺風作用下波浪潮流耦合三維泥沙數值模擬技術仍有待進一步研發，例如高分辨率臺風風場模擬、波-流耦合計算、細顆粒泥沙的沉降與輸移的模擬等等。為此，本文應用中尺度大氣模式WRF、第3代波浪模型MIKE21-SW、二維潮流模型MIKE21-HD和三維水動力泥沙模型FVCOM，針對臺風作用下泥沙輸移、海床演變、港池航道沖淤，通過二次開發，建立了臺風作用下波浪潮流耦合三維泥沙數值模擬系統。著重討論臺風作用下三維泥沙輸移的模擬，主要解決了波流耦合問題和三維模型中增加了波生近岸流、波浪對泥沙與運動的影響、泥沙的絮凝沉降和高含沙量時的制約沉降機制等影響因素，建立了臺風作用下波流耦合作用下三維泥沙數學模型，并將模型應用于連云港港30萬t級淺灘深開挖航道工程的研究。

1 三維泥沙數學模型FVCOM改進

1.1 FVCOM泥沙擴散方程

FVCOM采用三維波流耦合作用下的泥沙模型(ROMS_SED)的泥沙計算模式，將泥沙按粒徑大小分組，考慮黏性沙、非黏性沙，又稱為多組分泥沙模型。其懸沙濃度擴散方程為：

(1)

式中：Ci為第i類泥沙的懸沙濃度；AH為水平渦黏性系數；KH為垂向渦黏性系數；wi為第i類泥沙的沉降速率。

邊界條件：在自由表面處不考慮大氣沉降的顆粒；在海底邊界處，泥沙的沉降和再懸浮過程分別作為水中泥沙的源和匯處理，即：

式中：Di為懸浮泥沙的沉降通量，Ei為海底泥沙侵蝕通量，Ei的表達式為：

式中:Qi為泥沙侵蝕系數；Pb為海底孔隙率；τb為海底剪切率；τci為泥沙起動的臨界剪切率。Qi依賴于底部泥沙的物理化學特性，可在實驗室中測量得到，根據英國HR實測資料，對于黏性土，Qi=2×10-4～4×10-3kg/m2s，α=1.16。

1.2 泥沙模型的改進

針對連云港港淤泥質海岸開敞海域，模型中需要考慮黏性泥沙的絮凝沉降過程和高含沙量時的制約沉降機制。當泥沙濃度小于制約沉降臨界濃度時，絮凝沉降速度計算公式采用彭潤澤等基于長江口泥沙絮凝實驗的表達式，即：

(2)

當泥沙濃度大于或等于制約沉降臨界濃度時，泥沙沉速計算公式為：

(3)

式中：Cgel為膠凝濃度，根據細顆粒泥沙的相關研究，計算中Cgel取為250 kg/m3。

2 三維泥沙模型的建立與驗證

2.1 模型建立

數值計算采用二維大、中模型和三維小模型三重嵌套的方法，模擬計算連云港港及附近海域“韋帕”臺風作用下波流耦合三維泥沙輸移過程[6]。三維模型計算范圍(D03)北起嵐山頭，南至中山河附近，計算域為北緯34°22.1′ ～35°6.6′, 東經116°0.0′～120°8.8′，模型采用非結構化三角形網格，網格最小尺度為20 m，位于連云港主港區和主航道處，網格最大尺度為2 000 m，位于計算域的深水開邊界處。三維模型的邊界條件由中范圍臺風作用下二維波浪潮流耦合數學模型提供。

2.2 “韋帕”臺風作用下波流耦合三維泥沙模型的驗證

為了準確地模擬連云港港及其附近海域“韋帕”臺風作用下波流耦合的三維泥沙輸移過程，根據天津市氣象科學研究所應用WRF中尺度預報模式和美國環境預報中心NCEP歷史再分析數據，在天河一號大型計算機上模擬計算了整個“韋帕”臺風過程。計算中采用超高分辨率的三重網格嵌套技術，嵌套區域D01、D02、D03的分辨率分別為9、3、1 km，為“韋帕”臺風作用下波浪、潮流、泥沙輸移的模擬計算提供了準確的驅動風場[6]，計算時間為2007年9月17日0時～9月21日20時，共117 h。

1)三維含沙量驗證。

在三維泥沙數值模擬計算中，潮位、潮流、波浪、風場等均由二維嵌套模型提供[6]，三維模型外模計算時間步長取1.0 s，內模計算時間步長選取外模的5倍，垂向平均分為10層。底部摩阻系數范圍取0.002 5～0.003 5。根據“韋帕”臺風期間實測的波浪、水文資料[7]，進行了“韋帕”臺風作用下波流耦合的三維泥沙數學模型的率定與驗證。其中1#點位于-3.0 m等深線，2#位于-5.0 m等深線(如圖1所示)，每個測站分為表下0.5 m、底上0.5 m共3層進行觀測。其中，2#測站的分層含沙量過程的驗證結果如圖2所示，從驗證結果可以看出，計算值與實測值吻合較好。

圖1 水文觀測點位置Fig.1 Hydrological observation point location

圖2 2#觀測點含沙量驗證Fig.2 Verification of sediment concentration at observation point 2#

受“韋帕”臺風影響，2007年9月19日0時“韋帕”臺風遠在浙江東海海面上，距離連云港海域較遠，水體含沙量較小。水體含沙量逐漸增大，到9月20日12點，含沙量增至最大值(如圖2、圖3所示)。然后逐漸減小，主要特征表現為：最大含沙量時刻落后于最大波高時刻。根據三維泥沙數值模擬結果得到，海州灣內南部的高含沙量水體將隨著風暴潮流橫向越過主航道后進入徐圩海域，可成為連云港港主航道泥沙回淤的重要來源之一。因此，在“韋帕”臺風期間，海州灣內南部的高含沙水體對主航道回淤有比較重要的影響。綜上，本文建立的三維泥沙數學模型較好反映了“韋帕”臺風期間連云港海域泥沙分布特征及輸移過程，可為連云港港30萬t級淺灘深開挖航道工程的設計提供技術支撐。

圖3 2007年9月20日12時含沙量分布Fig.3 Sediment concentration distribution at 12 o′clock on September 20, 2007

2)航道回淤驗證。

根據2007年9月13日(“韋帕”臺風前)和2007年9月21日(“韋帕”臺風后)連云港港主航道斷面水深測量[8]，2次實測的航道水深沿程分布得到(如圖4所示)，風前航道水深沿程分布在0～1.5 km和1.5～3.0 km呈階梯狀變化，3.0 km以外水深逐漸增大；由于臺風期間航道產生淤積，風后航道水深沿程分布在0～3.0 km內基本一致，向外水深逐漸加大。

圖4 “韋帕”臺風前后主航道沿程水深Fig.4 Water depth along main channel before and after typhoon Wipha

對比2次水深測量結果可知，在航道里程0～1.3 km內，航道內的淤積強度基本一致，在0.3～0.5 m；航道里程1.3～3.0 km的淤積強度較大，約為前者的2倍。模擬計算得到的“韋帕”臺風作用下主航道回淤與實測較為一致，模擬計算的最大實質回淤強度為0.40 m。

3 “韋帕”臺風浪作用下三維泥沙數學模型應用

3.1 航道工程方案

連云港港30萬t級航道工程由連云港區航道、連云港區外航道、徐圩進港航道組成，其中，連云港區外航道外段寬350 m，底標高-23.0 m，外航道內段寬310 m，底標高-22.5 m；徐圩進港航道寬370 m，底標高-22.0 m，徐圩港區內航道寬210 m，底標高-13.3 m；徐圩二港池航道寬170 m，底標高-11.0 m。徐圩港區防波堤口門分為“平口”、“八字口”，結合不同港內圍墾工程方案布置，分為4個方案，本文以方案2、方案4為例進行分析，具體布置如圖5所示。

圖5 方案2、方案4布置圖Fig.5 Layout of plan 2 and 4

3.2 含沙量計算分析

圖6～圖11給出了工程方案2、方案4連云港海域在“韋帕”臺風期間代表時刻的表層和底層含沙量分布情況。2007年9月19日0時，“韋帕”臺風遠在浙江東海海面上，距離連云港海域甚遠，除灌河口含沙量較大外，連云港海域內的含沙量均不大(如圖6、圖9所示)。隨著臺風中心逐漸往北移動，臺風中心距離連云港越近，臺風浪也隨之加大，導致近岸水體含沙量逐漸增加，波浪掀沙、波流共同輸沙作用明顯，9月20日0時表層、底層含沙量均有明顯增大(如圖7、圖10所示)。至9月20日12時，連云港海域含沙量達到最大(如圖8、圖11所示)，滯后最大臺風浪出現時刻約5 h[6]。值得注意的是，由于該時刻的風暴潮流為由北向南的沿岸流，海州灣內的高含沙水體會越過連云港區主航道進入徐圩海域，與實測資料較為相符，即-5.0 m測點的含沙量從最高濃度逐漸減小后，仍然有大幅度增大的趨勢。

圖11 2007年9月20日12時方案4含沙量分布Fig.11 Plan 4 sediment concentration distribution at 12 o′clock on September 20, 2007

圖10 2007年9月20日0時方案4含沙量分布Fig.10 Plan 4 sediment concentration distribution at 0 o′clock on September 20, 2007

圖8 2007年9月20日12時方案2含沙量分布Fig.8 Plan 2 sediment concentration distribution at 12 o′clock on September 20, 2007

圖7 2007年9月20日0時方案2含沙量分布Fig.7 Plan 2 sediment concentration distribution at 0 o′clock on September 20, 2007

由于工程方案2、方案4港區圍填面積的不同，港池內部的水體含沙量也有差異，但與口門外相比，含沙量都較小。一方面，由于各方案中防波堤均已建成，防波堤的防浪效果較好，港池內的波高均不大，口門處波高最大值約為2.2 m，波浪自口門處逐漸遞減，港池內最大時刻波高介于1～2 m，和港池外相比小得多，難以引起港池內部泥沙大面積啟動和底部高濃度含沙層的產生。另一方面，圍填區域面積的不同，造成港池內的環流強度均有不同，其特征表現為港池圍填面積越小，環流強度越大。如工程方案2環流強度較大的時刻，港池內的環流區含沙量濃度稍高(圖6～8)；工程方案4港區圍填全部形成，環流強度最弱，港內含沙量較小(圖9～11)。因此，徐圩大環抱防波堤建成后，港池內的水體泥沙主要來源于港外輸入和港池內大強度環流造成的淺灘沖刷。

圖9 2007年9月19日0時方案4含沙量分布Fig.9 Plan 4 sediment concentration distribution at 0 o′clock on September 19, 2007

圖6 2007年9月19日0時方案2含沙量分布Fig.6 Plan 2 sediment concentration distribution at 0 o′clock on September 19, 2007

3.3 回淤計算分析

連云港主航道內外段的沿程回淤強度如圖12所示。由圖可知，航道最大回淤發生在航道里程3～5 km，相當于旗臺北防波堤外1～3 km內(旗臺北防波堤距離主航道拐點約2.2 km)，回淤強度約為0.40 m。總體來說，主航道外段回淤小于內段，航道里程5～40 km測點隨著距離口門越遠，回淤強度迅速減小，然后趨于平緩。由于受防波堤的阻擋作用，旗臺港區內部回淤較小，如航道里程0～3 km隨著離口門越遠回淤越小。

圖12 主航道內外段回淤強度Fig.12 Siltation intensity of the main channel

徐圩進港航道方案2、方案4沿程回淤強度如圖13所示。可見，徐圩進港航道回淤最大范圍出現在航道里程3～7 km，該段航道2個方案回淤強度差別不大；其中，最大回淤位于距離徐圩大環抱口門約5 km處，其風后最大回淤強度約為0.30 m。位于口門內的0、1 km處回淤比口門小，平均回淤約0.15 m。總體來說，徐圩進港航道最大回淤強度比主航道略小，但沿程回淤強度變化較為平緩。

圖13 方案2、4徐圩進港航道回淤強度Fig.13 Plan 2、4 siltation intensity of Xuwei entry channel

徐圩一港池航道方案2、方案4沿程回淤強度如圖14所示。其中，方案4與方案2相比，港區陸域圍填全部完成，港池內水域面積最小，一港池航道回淤強度明顯減小，并隨著與口門之間距離的增大而減小，尤其是一港池內航道里程0～4 km回淤強度明顯小于方案2。

圖14 方案2、4徐圩一港池航道回淤強度Fig.14 Plan 2、4 siltation intensity of Xuwei port #1

總之，與徐圩進港航道回淤相比，徐圩港區港池內回淤強度較小，主要是因為在臺風期間，防波堤掩護作用較好，港內水體含沙量較小；另外，在臺風風速較大期間，港池口門處的潮流流向幾乎平行于港池口門，使得直接進入港池內的泥沙較少。雖然徐圩航道靠近南部的灌河口海域，由于“韋帕”臺風期間，潮流的主流是由西北向東南沿岸流動，所以，灌河口海域的高含沙水體不易將泥沙輸移至徐圩航道工程海域，進一步說明波浪掀沙、波流共同輸沙的重要性。

4 結論

1)基于三重網格嵌套的超高分辨率(分辨率分別為9、3、1 km)“韋帕”臺風風場，應用第3代波浪模型MIKE21-SW、二維潮流模型MIKE21-HD和三維水動力泥沙模型FVCOM，以二維、三維大、中、小模型三重嵌套的方式，主要解決了波流耦合問題和三維模型中增加了波生近岸流、波浪對泥沙運動、泥沙的絮凝沉降和高含沙量時的制約沉降機制等影響因素，建立了臺風作用下波浪潮流耦合三維泥沙數學模型。

2)應用臺風作用下波浪潮流耦合三維泥沙數學模型，模擬了2007年“韋帕”臺風作用下連云港海域三維泥沙輸移過程。2#測站的分層含沙量過程的計算值與實測值吻合較好。“韋帕”臺風作用下連云港區主航道的回淤強度與實測較為一致，最大實質回淤強度為0.40 m。較好地反映在“韋帕”臺風作用下連云港海域的泥沙運動特征。

3)計算分析了方案2、方案4連云港港海域“韋帕”臺風作用下波流耦合三維泥沙運動特征，得到：連云港主航道航道最大回淤發生在航道里程3～5 km，相當于旗臺北防波堤外1 km～3km內回淤強度約為0.40 m。受防波堤的阻擋作用，旗臺港區內部回淤較小。徐圩進港航道最大回淤值比主航道略小，最大回淤強度約為0.30 m，但其沿程回淤強度變化較為平緩。徐圩港區港池內回淤強度較小，主要是因為臺風期間，防波堤掩護作用較好，港內水體含沙量較小。另外，由于“韋帕”臺風期間，風暴潮流的主流是由西北向東南沿岸流動，所以，灌河口海域的高含沙水體不易將泥沙輸運至徐圩航道工程海域，進一步說明波浪掀沙、波浪潮流共同輸沙的重要性。