山海關(guān)旅游海灘生態(tài)修復(fù)工程水動(dòng)力變化特性

匡翠萍，董博靈，謝華浪，楊燕雄

(1. 同濟(jì)大學(xué) 土木工程學(xué)院，上海 200092；2. 珠江水利科學(xué)研究院，廣東 廣州 510611；3. 河北省地礦局第八地質(zhì)大隊(duì)，河北 秦皇島 066001)

山海關(guān)位于河北省秦皇島市，是我國著名的旅游區(qū).其中，海灘旅游占有重要地位[1-2]，長達(dá)14 km的海岸線是山海關(guān)最重要的旅游資源之一.然而近些年來，灘面下蝕[3]、岸線后退[4]、砂粒粗化[5]等海灘侵蝕問題層出不窮，也因此成為了山海關(guān)海灘旅游發(fā)展的瓶頸.時(shí)間尺度上，海灘侵蝕問題可以分為兩類：一是短暫的、迅速的侵蝕過程，這種侵蝕過程通常是由罕見的強(qiáng)浪和風(fēng)暴潮引起的，相關(guān)研究內(nèi)容包括侵蝕過程[6]、侵蝕量級[7]、累積效應(yīng)[8]等；二是持續(xù)的、緩慢的侵蝕過程，這種侵蝕過程一般是由海域整體的動(dòng)力不平衡引起的，相關(guān)研究內(nèi)容包括侵蝕的地域差異[9]、侵蝕區(qū)域的水動(dòng)力分析[10]、海平面上升引起的侵蝕[11]等.相較于前者，海域動(dòng)力不平衡引發(fā)的海灘侵蝕由于其不可逆性而危害更大.為了改善山海關(guān)海域水動(dòng)力狀況，保護(hù)海灘旅游資源，山海關(guān)旅游海灘修復(fù)工程由此展開，而同時(shí)工程的建設(shè)必然對周邊環(huán)境產(chǎn)生影響[12-14].目前，針對秦皇島市的海灘侵蝕問題，如金夢海灣[15]、北戴河?xùn)|海灘[16]、中海灘和西海灘[17]，采取了拋沙養(yǎng)灘和布置離岸潛堤的措施，使海灘整體趨向穩(wěn)定.Kim等[18]通過野外監(jiān)測證實(shí)了水下人工魚礁對韓國東側(cè)海灘具有保灘促淤的作用.為了探究山海關(guān)旅游海灘生態(tài)修復(fù)工程對海域水動(dòng)力的影響，本研究采用Delft3D數(shù)值模擬軟件建立波流耦合數(shù)學(xué)模型，模擬和預(yù)測工程建設(shè)引起的水動(dòng)力變化特性.

1 數(shù)學(xué)模型

1.1 工程概況

山海關(guān)旅游海灘生態(tài)修復(fù)工程如圖1所示.針對白鷺島至唐子寨岸段海灘，近岸主要采取灘肩補(bǔ)沙，并輔以離岸潛堤阻擋波浪，并在后緣設(shè)計(jì)濱海景觀廊道和人工覆植沙丘等生態(tài)工程.其中，補(bǔ)沙岸線長度為2 km，沙灘寬度在原有基礎(chǔ)上增加30～50 m，沙粒的中值粒徑為0.46～0.69 mm，采用交會(huì)型剖面的形式；兩座離岸潛堤分別距岸240 m和350 m，長300 m，寬8.4 m，堤頂高程為-1.0 m，使用人工魚礁的形式進(jìn)行布置，以盡量減小工程對海域生態(tài)的負(fù)面效應(yīng)；濱海景觀廊道主要包括自行車道、綠化帶和木棧道，形成區(qū)域海岸帶整體建設(shè)的旅游功能與生態(tài)環(huán)境相互平衡.基于Delft3D建立波流耦合數(shù)學(xué)模型，分析研究工程造成的海域水動(dòng)力變化特征，為這類旅游生態(tài)修復(fù)工程的設(shè)計(jì)與建設(shè)提供一定科學(xué)依據(jù).

圖1 工程總體平面布置圖Fig.1 General engineering layout

1.2 模型方程

Delft3D是由WL Delft Hydraulics(荷蘭Delft水力學(xué)研究所)開發(fā)的一套多功能的數(shù)值模擬軟件，主要應(yīng)用于河流、河口、海岸等自由地表水環(huán)境地區(qū)的水動(dòng)力、水環(huán)境、泥沙運(yùn)動(dòng)等方面的科學(xué)研究.

Delft3D-FLOW潮流模型的基本方程采用淺水假定和Boussinesq近似，由不可壓縮的三維Navier-Stokes方程推導(dǎo)而來.

水深平均連續(xù)方程如下：

(1)

ξ和η為正交曲線坐標(biāo)系的兩個(gè)坐標(biāo)方向，ξ和η方向的動(dòng)量方程為

(2)

(3)

式中：ζ為基準(zhǔn)面以上的水位；t為時(shí)間；Gηη、Gξξ為正交曲線坐標(biāo)系與笛卡爾坐標(biāo)系之間的轉(zhuǎn)換系數(shù)；d為基準(zhǔn)面以下的水深；u、v分別為ξ和η方向的水深平均流速；w為z方向的流速；f為柯氏力系數(shù)；ρ0為水的相對密度；Pξ、Pη分別為ξ和η方向壓力梯度；Fξ、Fη分別為ξ和η方向上紊動(dòng)引起的雷諾應(yīng)力；Mξ、Mη分別為ξ和η方向上泄流或取水、構(gòu)筑物、波浪應(yīng)力等外部作用引起的動(dòng)量變化；vV為垂向渦黏系數(shù)；Q為單位面積的源匯項(xiàng)，包括降水、蒸發(fā)、流入和流出的水流引起的水量變化.

(4)

式中：qin和qout分別為單位體積的源項(xiàng)和匯項(xiàng)；P為降水引起的源項(xiàng)；E表示蒸發(fā)引起的匯項(xiàng).

Delft3D-WAVE波浪模型采用第三代SWAN波浪模型，即采用二維波浪密度譜描述隨機(jī)波浪，其波浪頻譜的變化以譜平衡方程描述，在笛卡爾坐標(biāo)系下，譜平衡方程為

(5)

式中：N為波浪作用密度；σ為相對頻率；θ為波浪傳播方向；cx、cy、cσ和cθ分別是x、y、σ、θ方向上的波浪傳播速度；式中第1項(xiàng)為波浪作用隨時(shí)間的變化，第2、3項(xiàng)為波浪傳播隨空間變化，第4項(xiàng)為地形和水流引起的多普勒頻移效應(yīng)，第5項(xiàng)表示地形和水流引起的波浪折射；S是源項(xiàng)(包括風(fēng)生浪、耗散、非線性波和波浪破碎等過程)，其中，風(fēng)能輸入可表示為

Sin(σ,θ)=A+BE(σ,θ)

(6)

式中：A和B取決于波浪的頻率、角度以及風(fēng)的速度、角度；E(σ,θ)為波浪能量譜密度.

1.3 模型參數(shù)

由于工程區(qū)域較小，為保證足夠的計(jì)算精度，采用渤海大模型、秦皇島中模型和山海關(guān)小模型三重嵌套的方式進(jìn)行計(jì)算，如圖2所示.大模型的網(wǎng)格數(shù)為450×350，網(wǎng)格最小空間步長約為2 km；中模型的網(wǎng)格數(shù)為200×330，網(wǎng)格最小空間步長約為350 m；小模型的網(wǎng)格數(shù)為285×460，網(wǎng)格最小空間步長約為50 m.

圖2 大、中、小模型網(wǎng)格Fig.2 Meshes of big, middle, and small models

渤海大模型海域開邊界由大連老虎灘和煙臺(tái)兩個(gè)潮位站所測得的潮位過程所控制，秦皇島中模型的邊界包括西邊界、南邊界和東邊界3條海域開邊界和1條陸地閉邊界，其潮位邊界條件由大模型提供，波浪邊界條件由實(shí)測和模型率定綜合得到.側(cè)固邊界采用不可滑移條件，利用動(dòng)邊界處理干濕交換過程.山海關(guān)小模型邊界的邊界處理與秦皇島中模型相同.大模型和中模型的地形數(shù)據(jù)采用近10年內(nèi)最新的地形數(shù)據(jù)，小模型的地形數(shù)據(jù)采用河北省地礦局第八地質(zhì)大隊(duì)于2015年12月(工程前)和2016年9月(工程后)兩次高密度(測量點(diǎn)間距小于3 m)的水下地形測量數(shù)據(jù).水平紊動(dòng)黏滯系數(shù)取為15 m2·s-1，曼寧系數(shù)取均值為0.015，計(jì)算時(shí)間步長取1.0 min，初始水位為平均海平面.

2 模型驗(yàn)證

通過將當(dāng)前時(shí)刻的Delft3D潮流模型與波浪模型輸出的計(jì)算結(jié)果，分別作為下一時(shí)刻的Delft3D波浪模型與潮流模型的輸入條件，實(shí)現(xiàn)兩個(gè)模型的相互耦合，并利用耦合后的數(shù)學(xué)模型，進(jìn)行模型驗(yàn)證.河北省地礦局第八地質(zhì)大隊(duì)于2015年9月23日、24日在工程海域附近測量得到相關(guān)實(shí)測數(shù)據(jù)，測站分布如圖3所示.根據(jù)2015年9月23日0:00

至9月25日0:00秦皇島和山海關(guān)潮位站的潮位過程資料進(jìn)行潮位驗(yàn)證，驗(yàn)證結(jié)果見圖4.根據(jù)2015年9月23日13:00至9月24日14:00測站BLD01、BLD02的實(shí)測潮流數(shù)據(jù)進(jìn)行潮流驗(yàn)證，驗(yàn)證結(jié)果見圖5.

圖3 測站點(diǎn)分布Fig.3 Positions of measurement points

a 秦皇島

b 山海關(guān)圖4 2015年9月23日0：00至9月25日0：00秦皇島和山海關(guān)潮位站潮位過程驗(yàn)證Fig.4 Tidal level validation at Qinhuangdao and Shanhaiguan

由于研究區(qū)域位于半日潮無潮點(diǎn)和日潮波腹區(qū)，潮汐變化復(fù)雜；實(shí)際測量正值小潮日，潮流流速較小，易受其他環(huán)境因素干擾，因此在個(gè)別時(shí)段，潮流流速的計(jì)算值和實(shí)測值存在一定的誤差.

Wilmott統(tǒng)計(jì)學(xué)方法[19]考慮了實(shí)測值與實(shí)測平均值的偏差、模型計(jì)算值和實(shí)測平均值的偏差這兩者的相關(guān)程度，其計(jì)算方法為

(7)

式中：s>0.65時(shí)結(jié)果為極好；0.5

a BLD01測站流速

b BLD01測站流向

c BLD02測站流速

d BLD02測站流向圖5 2015年9月23日13：00至9月24日14：00測站BLD01和BLD02的垂向平均流速和流向過程驗(yàn)證Fig.5 Depth-averaged velocity and direction validation at BLD01 and BLD02

表1 Wilmott統(tǒng)計(jì)學(xué)方法得到的模型評價(jià)結(jié)果Tab.1 Model evaluation result by descriptive statistic of Wilmott

可以看出，潮位和潮流驗(yàn)證的s值評價(jià)均為極好.由此可知，模型模擬值與實(shí)測值吻合較好.該模型可以用于模擬該區(qū)域?qū)嶋H的水動(dòng)力.

從動(dòng)力角度分析，波浪的動(dòng)力作用與波高的平方成正比.因此，選取以波高平方為依據(jù)的波浪統(tǒng)計(jì)值作為具有長期代表作用的常浪.由秦皇島波浪測站2011年3月至9月的資料，可得常浪條件為波高H1/3=0.62 m，周期T=4.68 s，波浪方向?yàn)?52.3°；另外，取2011年的最大有效波高作為強(qiáng)浪，故強(qiáng)浪條件為H1/3=1.50 m，T=5.50 s，波浪方向?yàn)?68.5°.經(jīng)過對模型的率定，模擬得到常浪和強(qiáng)浪波況下的波高、周期、波向的特征值均與實(shí)測特征值相一致(見表2).

表2 常浪和強(qiáng)浪波況下的實(shí)測值與計(jì)算值對比Tab.2 Comparison of measured value and model result of normal and strong wave

3 計(jì)算結(jié)果

利用經(jīng)過驗(yàn)證的波流耦合模型，得到了山海關(guān)旅游海灘修復(fù)工程建設(shè)前后的潮流場和波浪場，對工程前后的水動(dòng)力進(jìn)行分析和比較，得到旅游海灘生態(tài)修復(fù)工程手段的成效.

3.1 潮流

山海關(guān)海域位于M2和S2分潮的無潮點(diǎn)附近[20]，潮差較小，水動(dòng)力較弱，潮汐為正規(guī)全日潮，潮流為正規(guī)半日潮.漲潮流基本為WSW向，落潮流基本為ENE向，潮流為典型的往復(fù)流.以2015年6月夏季典型潮型為例，工程前山海關(guān)海域大小潮漲、落急時(shí)刻潮流場如圖6所示(圖中的坐標(biāo)系統(tǒng)為北京54坐標(biāo)系統(tǒng)).總體上，山海關(guān)海域大小潮漲、落急流速大小相差不大，漲潮流向大致為SW向，落潮流向大致為NE向，為周期性往復(fù)流.其中，漲潮潮流歷時(shí)約為6.1 h，落潮潮流歷時(shí)約為6.4 h，落潮潮流歷時(shí)大于漲潮潮流.流速平面分布上，流速由外海向岸邊遞減，外海流速約為0.3 m·s-1，近岸流速不足0.05 m·s-1.由于受西側(cè)樂島、東側(cè)白鷺島形成的雙岬角掩護(hù)作用，潮流繞過岬角進(jìn)入工程海域之后流速明顯減小，距離岸邊約600 m范圍內(nèi)漲、落急流速均在0.01～0.15 m·s-1范圍內(nèi).

以夏季典型大潮為例，潛堤建造前后工程海域的流速差值場如圖7所示.由圖可知，離岸潛堤附近的水體運(yùn)動(dòng)方向與離岸潛堤的布置方向平行，潛堤建設(shè)后潮流場僅有微小的改變.因此，潛堤不會(huì)減小掩護(hù)區(qū)海域內(nèi)的水體交換能力，對海域環(huán)境沒有負(fù)面影響，保證了旅游海灘海域水質(zhì).

a 大潮漲急

b 大潮落急

c 小潮漲急

d 小潮落急圖6 工程前大小潮漲、落急時(shí)刻流速場Fig.6 Depth-averaged velocity distribution at the moment of maximum flood and ebb of spring and neap tide before the project

a 漲急

b 落急圖7 工程前后的漲、落急時(shí)刻流速差值場(工程后減工程前)Fig.7 Current velocity changes induced by offshore submerged breakwaters

3.2 波浪

利用經(jīng)過驗(yàn)證的波流耦合模型，分別計(jì)算了平均海平面下常浪(波高0.62 m)和強(qiáng)浪(波高1.50 m)下的波浪場，在此基礎(chǔ)上研究工程對波浪的影響.

圖8為常浪和強(qiáng)浪下潛堤建成前后波高差值分布.由圖可知，常浪下，潛堤對其背后區(qū)域具有一定的掩護(hù)作用，波浪通過潛堤后發(fā)生波能損耗使波高略微減小，在離岸潛堤后方100 m范圍內(nèi)的波高削減程度約為0.05 m；強(qiáng)浪下，潛堤對波浪的削減作用較強(qiáng)，波浪通過潛堤后波高明顯減小，但同時(shí)由于波高增大，波浪的繞射現(xiàn)象也變得更為明顯，使得離岸潛堤對較遠(yuǎn)的近岸區(qū)域的波浪削減作用較小，在離岸潛堤后方約100 m范圍內(nèi)對強(qiáng)浪的波高削減度達(dá)0.2～0.4 m，往近岸方向潛堤對波高的削減作用越來越小.

a 常浪

b 強(qiáng)浪圖8 常浪和強(qiáng)浪下潛堤建成前后波高差值分布(工程后減工程前)Fig.8 Normal and strong wave height changes induced by offshore submerged breakwaters

為了統(tǒng)計(jì)分析潛堤掩護(hù)區(qū)不同區(qū)域的波高變化情況，根據(jù)波流耦合模型，在潛堤掩護(hù)區(qū)內(nèi)布置了30個(gè)波浪監(jiān)測站點(diǎn)，分為A和B兩個(gè)區(qū)域，分別位于兩個(gè)潛堤的波影區(qū)，見圖9.圖10為工程前后波浪監(jiān)測站點(diǎn)處的波高直方圖.由圖10可知，常浪下，A區(qū)東側(cè)區(qū)域幾乎不受影響，A區(qū)西側(cè)岸邊區(qū)域受養(yǎng)灘工程影響波高減小了0.3 m，在潛堤掩護(hù)區(qū)波高的減小均約為0.05 m；B區(qū)岸邊區(qū)域波高減小了0.3～0.4 m，在潛堤掩護(hù)區(qū)的減小幅度大于A區(qū)，為0.03～0.08 m.而在強(qiáng)浪下，在A區(qū)域A15測站距離潛堤最近，波高減幅達(dá)0.38 m；在B區(qū)域B11測站位于潛堤掩護(hù)區(qū)，波高減幅達(dá)0.22 m.因此，潛堤掩護(hù)區(qū)在常浪下總體的減小幅度約為10%，在強(qiáng)浪下總體的減小幅度約為20%.

圖9 波浪監(jiān)測站點(diǎn)分布Fig.9 Positions of wave observation points

a A區(qū)域常浪

b A區(qū)域強(qiáng)浪

c B區(qū)域常浪

d B區(qū)域強(qiáng)浪圖10 工程前后波浪監(jiān)測站點(diǎn)處的波高對比Fig.10 Comparison of wave height on observation points before and after the project

4 結(jié)論

(1) 生態(tài)修復(fù)工程實(shí)施后潮流場僅有微小的改變，離岸潛堤幾乎沒有減小掩護(hù)區(qū)海域內(nèi)的水體交換能力，保持了海域內(nèi)水動(dòng)力環(huán)境.

(2) 對于常浪(入射波高0.62 m)，離岸潛堤的存在使其后方100 m內(nèi)的波高減小約0.05 m，而對于強(qiáng)浪作用(入射波高1.50 m)，堤后100 m內(nèi)波高減小達(dá)0.2～0.4 m，離岸潛堤對強(qiáng)浪的防護(hù)作用更大，可大大減少風(fēng)暴天氣下的海灘侵蝕.

(3) 山海關(guān)旅游海灘生態(tài)修復(fù)工程的建設(shè)既保護(hù)了工程海域水動(dòng)力環(huán)境，又削弱了波浪特別是強(qiáng)浪對海灘的侵蝕，有利于維持當(dāng)?shù)芈糜魏┑纳鷳B(tài)環(huán)境與穩(wěn)定.