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

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

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

1 數(shù)學模型

1.1 工程概況

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

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

1.2 模型方程

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

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

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

(1)

ξ和η為正交曲線坐標系的兩個坐標方向，ξ和η方向的動量方程為

(2)

(3)

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

(4)

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

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

(5)

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

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

(6)

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

1.3 模型參數(shù)

由于工程區(qū)域較小，為保證足夠的計算精度，采用渤海大模型、秦皇島中模型和山海關(guān)小模型三重嵌套的方式進行計算，如圖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

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

2 模型驗證

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

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

圖3 測站點分布Fig.3 Positions of measurement points

a 秦皇島

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

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

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

(7)

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

a BLD01測站流速

b BLD01測站流向

c BLD02測站流速

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

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

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

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

表2 常浪和強浪波況下的實測值與計算值對比Tab.2 Comparison of measured value and model result of normal and strong wave

3 計算結(jié)果

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

3.1 潮流

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

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

a 大潮漲急

b 大潮落急

c 小潮漲急

d 小潮落急圖6 工程前大小潮漲、落急時刻流速場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 工程前后的漲、落急時刻流速差值場(工程后減工程前)Fig.7 Current velocity changes induced by offshore submerged breakwaters

3.2 波浪

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

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

a 常浪

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

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

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

a A區(qū)域常浪

b A區(qū)域強浪

c B區(qū)域常浪

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

4 結(jié)論

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

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

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