珊瑚島礁孤立波爬坡的平面二維數值模擬研究

劉林平，劉維杰, 2，孫志林

(1. 浙江大學 海洋學院，浙江 舟山 316021；2. 南京水利科學研究院 水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210029)

海嘯波通常由海底地震、火山爆發或者山體滑坡等引起，并通過沿岸爬坡過程侵襲陸地，造成極大的破壞，因此研究海嘯波在近岸地區的傳播以及爬坡過程具有重要的科學和工程意義。盡管實際海嘯波通常以瞬態非周期性的形式出現，但孤立波被認為可模擬海嘯首波的一些重要特征，因此學者通常采用孤立波進行海嘯波的物模試驗以及數值模擬研究[1-5]。自印度洋海嘯之后，珊瑚礁在抵御海嘯襲擊中的作用引起了人們的重視，但眾多災后調查發現珊瑚礁對于海嘯災害的防御程度不一[6]。因此，為了更好地了解珊瑚礁在抵御海嘯侵襲過程中的作用，學者們進行了相應的物理模型試驗和數值模擬來探究珊瑚礁不同地形和水動力參數對海嘯災害的影響。

在物理模型試驗方面，Quiroga等[7]通過大比尺水槽試驗研究了礁床粗糙度對孤立波在礁坪上傳播的影響；Yao等[8]在水槽試驗中采用圓柱體陣列近似模擬了珊瑚礁大糙率表面，并提出了預測珊瑚礁海岸孤立波爬高的經驗公式。相比物理模型試驗，數值模擬能更便捷地調整珊瑚礁地形參數(如礁坪寬度、礁前斜坡、礁后斜坡等)和水動力參數(如礁坪水深、入射波高、珊瑚礁糙率等)。過去的幾十年中，Boussinesq方程因其計算精度和效率較高的優勢，在珊瑚礁波浪水動力數值模擬研究中應用最為廣泛[9-14]。對于珊瑚礁孤立波傳播和爬坡，Zhou等[15]利用Boussinesq波浪模型研究了珊瑚礁地形和水動力參數對孤立波在礁緣附近的反射、破碎和透射的影響。Yao等[8]和Ning等[11]利用Boussinesq波浪模型研究了地形和水動力參數對礁后孤立波爬坡高度的影響。Liu等[16]利用Boussinesq波浪模型進一步發現，礁前斜坡的坡度對孤立波爬坡帶內最大動量通量的影響與對爬高高度的影響有所不同。這些研究有助于理解不同因素對緩解珊瑚礁海岸海嘯災害的重要性，也展現了Boussinesq模型模擬珊瑚礁陡變地形上孤立波傳播和爬坡的能力。

然而，目前有關珊瑚礁對海嘯災害影響的研究仍只局限于水平一維(二維地形)的研究，即假設地形在沿岸方向無變化，波浪在沿岸方向無流速，實際更適用于沿大陸分布的岸礁地形上的波浪運動。但對于諸如南海島礁的孤立島嶼，波浪的繞射和折射現象不容忽視，水平一維研究無法充分展示在此情況下珊瑚礁各地形、水動力因素對礁后孤立波最大爬高分布的影響。因此為了彌補現有研究的不足，采用激波捕捉類Boussinesq模型FUNWAVE-TVD[17]，對三維島礁地形上的孤立波傳播和爬坡開展平面二維數值模擬，利用已有物理模型試驗數據驗證了模型模擬孤立波在三維島礁地形上傳播與爬坡的能力，并進行了一系列的數值試驗探究了珊瑚礁對島嶼周圍孤立波最大爬高分布的影響。

1 模型介紹

采用的數值模型是Shi等[17]提出的FUNWAVE-TVD, 控制方程基于Chen等[18]提出的完全非線性Boussinesq方程，并結合Kennedy等[19]提出的參考水深法，其連續性和動量控制方程：

(1)

(2)

(3)

其中，H=h+η表示當地總水深，h為靜水深。

(4)

(5)

(6)

(7)

式中：V1和V2表示Boussinesq方程的色散項，V3表示的是色散項O(μ2)對w×u=wiz×u的貢獻，其中iz表示z方向的單位向量，w表示垂向速度。R表示底部摩擦和網格側向湍流引起的擴散和耗散項，其中底摩擦項為結合曼寧系數的二次摩擦項：

(8)

其中，n表示曼寧系數。

在FUNWAVE-TVD中，上述控制方程被進一步調整為守恒形式且空間上的離散采用有限體積和有限差分的混合數值格式。時間步長采用非線性空間離散化的三階強穩定性(SSP)Runge-Kutta方法，并根據Courant-Friedrichs-Lewy(CFL)標準選擇自適應時間步長。基于混合數值格式，波浪破碎處理采用激波捕捉法，即當波面升高同當地水深的比值達到某一臨界值時，將Boussinesq方程退化為帶有TVD格式的非線性淺水方程進行波浪破碎計算，忽略Boussinesq方程中的高階非線性項和色散項，根據Tonelli等[20]的研究，模型中判斷波浪破碎的臨界值為0.8。干濕界面處通過修正黎曼解的特征速度實現海岸動邊界波浪爬坡的計算。有關模型的具體內容可進一步參考文獻[17]。

2 模型驗證

對于孤立波在二維珊瑚礁陡變地形上傳播和爬坡的模擬，本模型已得到了充分的驗證[5,11]，因此這里僅利用Briggs等[2]的港池試驗數據來驗證模型在模擬三維島嶼地形上孤立波傳播和爬坡的可靠性。基于Briggs等的試驗，數值港池設置的平面和剖面布置如圖1所示。數值港池長為35 m，寬30 m，港池中設置有坡度為1∶4的圓臺來近似實際島嶼，圓臺的中心在x=23 m，y=15 m的位置處，造波帶設置在x=9 m的位置，港池水深為0.32 m，港池兩側設置有海綿層來吸收反射波。物理模型試驗中共設置了27個浪高儀來測量波面高程的變化，選取了其中四個浪高儀的數據用于模型驗證，分別為迎浪面的6號和9號浪高儀，島嶼背浪面的22號浪高儀和島嶼側面的16號浪高儀，其具體位置如圖1所示。同時，試驗也在圓臺周圍布置了爬高計來測量島嶼周圍的孤立波最大爬坡高度。

圖1 基于Briggs試驗的數值模型設置Fig. 1 Numerical setup based on Briggs’ experiments

數值模擬中，網格設置為0.05 m×0.05 m, 曼寧系數設置為0.012以表征島嶼物理模型表面光滑混凝土的糙率，輸出數據的時間步長為0.02 s。圖2展示了入射孤立波波高為0.014 m時四個浪高儀位置處本文模型預測的波面隨時間變化與試驗數據的對比，同時也展示了Lynett等[21]與Fuhrman等[22]所用Boussinesq模型的計算結果。如圖2所示，三種模型計算結果相近且均與試驗數據的吻合程度較好，但在浪高儀9的位置，本文模型對于波峰略有高估，這可能是由于本文模型采用激波捕捉處理波浪破碎造成的。圖3展示了模型預測島嶼周圍的相對最大爬高(最大爬坡高度與入射波高的比值)與試驗數據的對比，同時也展示了Lynett等[21]與Fuhrman等[22]所用Boussinesq模型的計算結果。圖中選取了圓臺一側十個測點，橫坐標表示測點與圓臺圓心的連線與入射波方向的夾角α的度數。如圖3所示，相比Lynett等[21]與Fuhrman等[22]的計算結果，本文模型總體上更為準確地模擬了島嶼周圍的最大爬高分布特征，表明本文模型的混合數值格式在預測波浪爬坡的可靠性和優越性。圖3中島嶼迎浪面和背浪面都出現比較大的爬高，且背浪面出現的爬高比迎浪面的更大，這一現象主要是波浪經過島嶼地形時發生繞射與折射，流向發生改變，波能會在背浪面匯集，使得背浪面爬高急劇增加[3]。

圖2 模型計算的波面高程隨時間變化與試驗結果對比Fig. 2 Comparison of computed and measured time series of surface elevations

圖3 模型預測島嶼周圍相對最大爬高與試驗結果對比Fig. 3 Comparison of computed and measured maximum runup heights around the island

3 數值試驗設置

在驗證了模型預測島嶼周圍孤立波傳播和爬坡的能力后，設置了一系列現場尺度的數值港池試驗用于探究珊瑚礁地形和水動力因素對島礁周圍孤立波最大爬高的影響。數值港池試驗設置如圖4所示，港池長17 km，寬16 km，網格大小為10 m×10 m。港池中設置有概化的圓形珊瑚島礁，島礁原型為印度尼西亞(南緯8°25′32″，東經122°30′30″)的Babi島(如圖5所示)，該島四周發育有成熟珊瑚礁。利用Google earth得到Babi島的海岸線周長約為6 400 m，珊瑚礁礁緣的周長約為9 000 m，將其概化為圓形得到海岸線直徑D約為2 044 m，礁坪寬約為400 m。島礁的中心位于x=8 km，y=8 km的位置，內部造波帶位于x=3 km。在上述島礁參數的基礎上，共設置了4組數值試驗，分別研究水動力(入射波高H、礁坪水深hr、珊瑚礁糙率n2)及珊瑚礁地形(礁坪寬度w、礁前斜坡坡度Cotθ、礁后斜坡坡度Cotβ)參數對島礁周圍孤立波爬坡的影響。各組數值試驗工況如表1所示，每組工況參數的取值變化范圍均根據現有文獻報道設定，其中礁坪寬度的變化范圍(50～1 000 m)、礁坪水深的取值(3 m)和變化范圍(1～5 m)、礁前(后)斜坡坡度取值(4)和變化范圍(1～20)根據Quataert等[23]總結的珊瑚礁參數設定，珊瑚礁糙率的取值(0.09)和變化范圍(0.02～0.09)根據Gelfenbaum等[6]的研究確定，入射波高取值(2 m)根據Ford等[24]的文獻報道確定，入射波高變化范圍(1～6 m)以及礁前水深hd的取值(60.5 m)根據Briggs等[2]物理模型試驗中采用的入射孤立波的非線性系數范圍(入射波高與水深的比值)而確定。此外，港池底床的曼寧系數設置為n1=0.02(沙質底床)，礁后斜坡的曼寧系數設置為n3=0.04(礫石)，數值模擬中數據輸出的時間步長為1 s。

圖4 數值試驗設置Fig. 4 Numerical experiment setup

圖5 巴比島 (Google Earth)Fig. 5 Babi island (Google Earth)

表1 數值試驗工況

4 結果與討論

圖6展示了組1珊瑚島礁四周孤立波最大爬高隨入射波高的變化。如圖6所示，珊瑚島礁四周的最大爬高隨著入射波高H的增大而不斷增大，入射波高為6 m的迎浪面(α=0°)最大爬高相比入射波高為1 m的計算結果增大約250%，背浪面(α=180°)最大爬高增大約125%。圖中結果也顯示入射波高越大爬高增大的幅度越小，這主要是因為波高越大波浪的破碎程度越強烈，能量損耗更大，一定程度上減弱了爬坡高度的增長，該結果也與前人水平一維研究[3,11]的成果類似。

圖6 島礁四周孤立波最大爬高隨入射波高的變化Fig. 6 Variation of maximum runup height around the reef-fringed island with different incident wave heights

圖7展示了組2珊瑚島礁四周孤立波相對最大爬高隨礁坪水深的變化，同時也展示了無珊瑚礁(僅圓臺島嶼)時，島嶼四周的相對最大爬高分布情況。如圖7所示，珊瑚島礁周圍的相對最大爬坡高度均隨礁坪深水的增大而不斷增大，當礁坪水深為1 m時，相比于無珊瑚礁島嶼，迎浪面(α=0°)最大爬高可降低70%，背浪面(α=180°)最大爬高可降低56%，而礁坪水深增大至5 m時，相比于無珊瑚礁島嶼，迎浪面(α=0°)最大爬高仍可降低約36%，但背浪面波能匯集區域的最大爬高已經與無珊瑚礁時相差無幾，表明礁坪水深增大到一定程度時，島礁背浪面的孤立波爬坡高度將不受島嶼周圍珊瑚礁的影響。

圖7 島礁四周孤立波相對最大爬高隨礁坪水深的變化Fig. 7 Variation of relative maximum runup height around the reef-fringed island with different reef flat water depths

圖8展示了組3中珊瑚島礁四周孤立波相對最大爬高隨礁坪寬度的變化。如圖8所示，珊瑚島礁周圍的相對最大爬坡高度均隨礁坪寬度的增大而不斷減小。當礁坪寬度為800 m時，相比于無珊瑚礁地形，迎浪面(α=0°)最大爬高降低約71%，背浪面 (α=180°) 最大爬高降低約56%。當礁坪寬度為100 m，整個珊瑚島礁四周的相對最大爬高已與無珊瑚礁時十分接近，甚至當礁坪寬度只有50 m時，珊瑚島礁的相對最大爬高會高于無珊瑚礁時的爬高，該結果也與部分學者[6]的觀點一致，珊瑚礁的存在未必總是能降低波浪災害，其陡變地形引起的波浪強非線性作用有時可能會加劇波浪災害，尤其是珊瑚礁本身的能量損耗作用較弱的時候。

圖8 島礁四周孤立波最大爬高隨礁坪寬度的變化Fig. 8 Variation of relative maximum runup height around the reef-fringed island with different reef flat widths

圖9展示了組4珊瑚島礁四周孤立波相對最大爬高隨礁前斜坡坡度的變化。如圖9所示，礁前斜坡坡度對島礁周圍的爬高分布影響幾乎沒有影響，該結果也與前人水平一維研究[5,11]的成果類似。盡管礁前斜坡變緩會降低波浪反射，但礁前斜坡變緩也會使波浪破碎點向外海方向移動，使得波浪破碎提早發生，能量損耗更多。

圖9 島礁四周孤立波最大爬高隨礁前斜坡坡度的變化Fig. 9 Variation of relative maximum runup height around the reef-fringed island with different fore-reef slopes

圖10展示了組5珊瑚島礁四周孤立波相對最大爬高隨珊瑚礁糙率的變化。如圖10所示，珊瑚島礁周圍的相對最大爬坡高度均隨珊瑚礁糙率的減小而不斷增大，整個珊瑚島礁周圍，相比于無珊瑚礁地形，珊瑚礁糙率為n2=0.09時迎浪面(α=0°)最大爬高降低約36%，背浪面(α=180°)最大爬高降低約56%。當珊瑚礁糙率減小至n2=0.03時，除迎浪面部分區域(α<30°)外，島礁周圍的相對最大爬高已經與無珊瑚礁時相差無幾，這與礁坪寬度的影響相似。而當珊瑚礁糙率減小至n2=0.02時，島礁周圍的相對最大爬高甚至會高于無珊瑚礁時的爬高，該現象的原因應與圖8中礁坪寬度為50 m的情況類似，當珊瑚礁本身的能量損耗作用較弱的時候，其陡變地形引起的波浪強非線性作用反而會加劇波浪災害。值得一提的是，n2=0.02通常用于近似表征珊瑚礁因白化而徹底死亡以后的表面糙率[6]，由此可見，氣候變化引起的珊瑚礁白化不僅會降低珊瑚礁對海嘯災害的防御作用，甚至會加重致災程度。

圖10 島礁四周孤立波最大爬高隨珊瑚礁糙率的變化Fig. 10 Variation of relative maximum runup height around the reef-fringed island with different reef roughnesses

圖11展示了組6珊瑚島礁四周孤立波相對最大爬高隨礁后斜坡坡度的變化。如圖11所示，珊瑚島礁周圍的相對最大爬坡高度隨礁后斜坡坡度的變緩而有所減小，礁后斜坡坡度為20的迎浪面(α=0°)最大爬高相比礁后斜坡坡度為2的計算結果降幅約38%，背浪面(α=180°)最大爬高降幅約33%，這主要是由于隨著礁后斜坡的變緩，波浪在水平方向上傳播的距離變大，波浪在整個爬坡過程中會有更多的摩擦損耗。

圖11 島礁四周孤立波最大爬高隨礁后斜坡坡度的變化Fig. 11 Variation of relative maximum runup height around the reef-fringed island with different back-reef slopes

5 結 語

基于驗證后的激波捕捉類Boussinesq模型FUNWAVE-TVD對孤立波在三維珊瑚礁島嶼附近的傳播和爬坡進行了現場尺度的水平二維數值模擬，分別研究了入射波高、礁坪水深、礁前斜坡、礁坪寬度、珊瑚礁糙率、礁后斜坡對島礁周圍孤立波最大爬坡高度的影響。結果表明：總體上珊瑚礁的存在可有效降低島嶼四周孤立波的最大爬坡高度，入射波高、礁坪水深、礁坪寬度、珊瑚礁糙率是影響珊瑚島礁四周孤立波爬坡分布的主要因素，島礁四周最大爬坡高度會隨入射波高和礁坪水深的增大、礁坪寬度和珊瑚礁糙率的減小而不斷增大。當礁坪水深增大到一定程度時，島礁背浪面的孤立波爬坡高度將不受島嶼周圍珊瑚礁的影響，而當礁坪寬度和珊瑚礁糙率減小至一定程度時，會出現島礁四周最大爬高高于無珊瑚礁時爬高的現象。礁后斜坡的變緩會使島礁周圍的最大爬高有所減小，而礁前斜坡坡度對珊瑚島礁周圍的最大爬高幾乎沒有影響。