山海關旅游海灘生態修復工程水動力變化特性

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

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

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

1 數學模型

1.1 工程概況

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

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

1.2 模型方程

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

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

水深平均連續方程如下：

(1)

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

(2)

(3)

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

(4)

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

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

(5)

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

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

(6)

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

1.3 模型參數

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

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

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

2 模型驗證

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

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

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

a 秦皇島

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

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

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

(7)

式中：s>0.65時結果為極好；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統計學方法得到的模型評價結果Tab.1 Model evaluation result by descriptive statistic of Wilmott

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

從動力角度分析，波浪的動力作用與波高的平方成正比.因此，選取以波高平方為依據的波浪統計值作為具有長期代表作用的常浪.由秦皇島波浪測站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°.經過對模型的率定，模擬得到常浪和強浪波況下的波高、周期、波向的特征值均與實測特征值相一致(見表2).

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

3 計算結果

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

3.1 潮流

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

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

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 波浪

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

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

a 常浪

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

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

圖9 波浪監測站點分布Fig.9 Positions of wave observation points

a A區域常浪

b A區域強浪

c B區域常浪

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

4 結論

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

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

(3) 山海關旅游海灘生態修復工程的建設既保護了工程海域水動力環境，又削弱了波浪特別是強浪對海灘的侵蝕，有利于維持當地旅游海灘的生態環境與穩定.