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

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

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

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

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

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

1 模型介紹

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

(1)

(2)

(3)

其中，H=h+η表示當(dāng)?shù)乜偹睿琱為靜水深。

(4)

(5)

(6)

(7)

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

(8)

其中，n表示曼寧系數(shù)。

在FUNWAVE-TVD中，上述控制方程被進(jìn)一步調(diào)整為守恒形式且空間上的離散采用有限體積和有限差分的混合數(shù)值格式。時(shí)間步長采用非線性空間離散化的三階強(qiáng)穩(wěn)定性(SSP)Runge-Kutta方法，并根據(jù)Courant-Friedrichs-Lewy(CFL)標(biāo)準(zhǔn)選擇自適應(yīng)時(shí)間步長?；诨旌蠑?shù)值格式，波浪破碎處理采用激波捕捉法，即當(dāng)波面升高同當(dāng)?shù)厮畹谋戎颠_(dá)到某一臨界值時(shí)，將Boussinesq方程退化為帶有TVD格式的非線性淺水方程進(jìn)行波浪破碎計(jì)算，忽略Boussinesq方程中的高階非線性項(xiàng)和色散項(xiàng)，根據(jù)Tonelli等[20]的研究，模型中判斷波浪破碎的臨界值為0.8。干濕界面處通過修正黎曼解的特征速度實(shí)現(xiàn)海岸動(dòng)邊界波浪爬坡的計(jì)算。有關(guān)模型的具體內(nèi)容可進(jìn)一步參考文獻(xiàn)[17]。

2 模型驗(yàn)證

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

圖1 基于Briggs試驗(yàn)的數(shù)值模型設(shè)置Fig. 1 Numerical setup based on Briggs’ experiments

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

圖2 模型計(jì)算的波面高程隨時(shí)間變化與試驗(yàn)結(jié)果對比Fig. 2 Comparison of computed and measured time series of surface elevations

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

3 數(shù)值試驗(yàn)設(shè)置

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

圖4 數(shù)值試驗(yàn)設(shè)置Fig. 4 Numerical experiment setup

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

表1 數(shù)值試驗(yàn)工況

4 結(jié)果與討論

圖6展示了組1珊瑚島礁四周孤立波最大爬高隨入射波高的變化。如圖6所示，珊瑚島礁四周的最大爬高隨著入射波高H的增大而不斷增大，入射波高為6 m的迎浪面(α=0°)最大爬高相比入射波高為1 m的計(jì)算結(jié)果增大約250%，背浪面(α=180°)最大爬高增大約125%。圖中結(jié)果也顯示入射波高越大爬高增大的幅度越小，這主要是因?yàn)椴ǜ咴酱蟛ɡ说钠扑槌潭仍綇?qiáng)烈，能量損耗更大，一定程度上減弱了爬坡高度的增長，該結(jié)果也與前人水平一維研究[3,11]的成果類似。

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 語

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