匡翠萍， 王 潔， 董智超， 劉會欣， 朱 磊

（1.同濟大學 土木工程學院，上海 200092；2.中交天津港灣工程研究院有限公司，天津 300222；3.河北省地質礦產勘查開發(fā)局 第八地質大隊，河北 秦皇島 066001）

潟湖作為陸地水域和海洋水域連接的特殊地貌形態(tài)，受到自然和人類活動的共同影響，擁有豐富的物種多樣性和自然資源，其生態(tài)系統(tǒng)寶貴且敏感［1-2］。隨著沿海經濟的快速發(fā)展，養(yǎng)殖廢水，漁港排污和入湖河流排污使?jié)暫廴締栴}日趨嚴重，藻類等富營養(yǎng)化問題頻發(fā)［3-5］。七里海潟湖是國內僅存的現(xiàn)代潟湖之一，也是昌黎黃金海岸國家級自然保護區(qū)的主要保護對象之一。然而從20世紀70年代開始，由于人類高強度地進行大規(guī)模的圍墾和養(yǎng)殖活動，使得七里海潟湖淤積嚴重，水域面積嚴重萎縮，自然生態(tài)系統(tǒng)遭到嚴重破壞。

探明潟湖水動力和水體交換的變化機制，是保證潟湖生態(tài)修復和生態(tài)保護系統(tǒng)可持續(xù)發(fā)展的基礎［6］。Dias等［7］基于拉格朗日法研究了葡萄牙Riade Aveiro潟湖的水體交換能力，得出潟湖整體的水體交換能力較強，但潟湖內部河道處的水體交換能力較弱。Jeyar等［8］在摩洛哥Nador潟湖的水體交換研究中采用對流擴散模型，結果表明大潮較小潮條件下水體交換的過程更快。Cucco等［9］和Umgiesser等［10］分別針對Venice潟湖和Curonian潟湖采用歐拉法進行水體交換能力的研究，水體交換能力季節(jié)性特征顯著，風和徑流是水體交換能力主要的影響因素。Kuang等［11］基于歐拉法研究了七里海潟湖水體交換能力的季節(jié)性變化特征及其影響因素，認為徑流和潮位與水體交換能力的影響呈明顯的二階指數(shù)函數(shù)分布?？锎淦嫉龋?］探究了清淤疏浚工程對七里海潟湖的影響，表明清淤疏浚工程在提升潟湖納潮量的同時有助于改善潟湖的水體交換能力。本文基于水動力和物質輸運數(shù)學模型，剖析長時間尺度下海平面上升、徑流變化和人類活動作用對七里海潟湖不同時期水體交換能力變化的影響程度，為七里海潟湖濕地生態(tài)保護修復技術體系的建設提供參考依據(jù)，同時也可以為世界范圍內不同區(qū)域的淺水潟湖的水體交換能力研究和生態(tài)修復技術體系的完善提供一定的理論基礎。

1 研究方法

本研究采用歐拉法進行水體沖刷時間的計算。通過分析研究區(qū)域內整體的初始水體質量的變化過程，得到整體的剩余質量的指數(shù)衰減函數(shù)，將整體的初始水體質量衰減至初始質量1/e（約為37%）所需的時間確定為水體沖刷時間。水體沖刷時間越短，表明區(qū)域的水體交換能力越強。故基于變化的地形，綜合考慮氣候變化引起的海平面上升和徑流的變化，同時考慮工程的影響，利用MIKE21軟件建立七里海潟湖二維水動力和物質輸運模型數(shù)學模型［12］，利用有限體積法求解二維不可壓縮Reynolds平均Navier-Stokes方程，考慮了Boussinesq近似、淺水假定和靜水壓力的假定，并耦合保守物質的物質守恒方程建立物質輸運模型，用于評價潟湖的水體交換能力。

數(shù)學模型的具體控制方程如下：

連續(xù)性方程為

式中：t為時間；x和y為笛卡爾坐標；u和v分別為x和y方向垂向平均流速分量；h為總水深，h=d+η，d為靜止水深，η為水位；S為源流量。

物質輸運模型方程為

式中：cm為物質濃度；Dcx、Dcy為物質的擴散系數(shù)；cs為源的物質濃度；kp為衰減系數(shù)。

1.1 計算區(qū)域與模型設置

七里海潟湖模型平面坐標系選用北京54高斯克魯格坐標系，采用三重非結構嵌套網格（圖1）。大模型為渤海模型，中模型模擬區(qū)域為秦皇島海域，七里海潟湖模型為加密后的小模型。網格包含6 656個節(jié)點和11 818個單元，空間步長為5～50 m。模型中變化的地形來源于七里海潟湖地貌演變的模擬結果，該模擬基于海平面和徑流數(shù)據(jù)及衛(wèi)星遙感影像，耦合水動力模型，綜合考慮多年來海平面上升過程、徑流變化及人類活動對地貌的影響。潟湖岸邊以岸線即固壁邊界進行考慮，外海開邊界位于新開口潮汐通道口門處，采用秦皇島海域中模型計算得到的Flather邊界進行控制，4條入湖河流流量為實測月均流量。

七里海潟湖由4條河流流入，流量變化采用文獻［13］及全國水文統(tǒng)計年報中搜集的灤河下游灤縣站從1950—2018年河流徑流數(shù)據(jù)河流流量進行換算得到。風應力采用歐洲中期天氣預報中心（ECMWF）給出的海面10 m以上、間隔為6 h的風場資料。底部摩擦力曼寧糙率系數(shù)由水深和底部中值粒徑共同控制，取值范圍為0.011～0.015，時間步長由模型自動調節(jié)為0.001～60 s，CFL（courantfriedrichs-lewy）數(shù)限值為0.8。水平渦黏系數(shù)采用Smagorinsky公式計算。模型采用干濕動邊界處理技術，模型中干點臨界水深取0.005 m，濕點臨界水深取0.05 m。

1.2 數(shù)學模型驗證

模型驗證采用2016年9月23日—2016年10月23日七里海潟湖潮汐通道內的潮位驗證點WL的實測數(shù)據(jù)以及2016年9月23日11：00—2016年9月24日11：00的潮流驗證點C1～C3（測點位置見圖1）的實測數(shù)據(jù)對七里海潟湖水動力模型進行驗證，驗證結果如圖2所示，其中潮流驗證點C1位于潟湖內， C2位于潮汐通道內， C3位于秦皇島海域。

圖1 研究區(qū)域、網格及驗證點位置Fig.1 Study area, computational grid， and positions of validated points

圖2 潮位及潮流驗證Fig.2 Verification of tidal level and tidal current

采用Wilmott［14］提出的效率評價系數(shù)對計算結果進行定量評價，Skill數(shù)的具體計算方程如下：

式中：S為Skill數(shù)；M為模型計算值；D為實測值為實測平均值；N為樣本數(shù)量。當Skill值為1時，代表模型計算值和實測值之間完全一致；Skill值大于0.65時，表示模型計算結果為“極好”；Skill值在0.65至0.50之間時，表示模型計算結果為“非常好”；Skill值在0.50至0.20之間時，表示模型計算結果為“好”；Skill值小于0.20時，表示模型計算結果為“差”。各測站Skill值如表1所示，可見各測站結果均為極好，表示模型能夠較好地模擬實際情況。

表1 模型評價Tab.1 Evaluation of model

1.3 水體交換能力研究

采用歐拉法計算水體沖刷時間進行水體交換能力的研究。具體方法是利用數(shù)學模型中保守物質的輸運模型，求解對流擴散方程，將研究的水體染成單位濃度的示蹤劑，模擬區(qū)域內水體物質濃度（cm）隨時間的變化過程。通過M=Vcm可得到水體質量（M）隨時間的變化過程，其中V為研究區(qū)域水體的體積，剩余初始水體質量隨時間的變化過程可以通過最小二乘法擬合，得到剩余質量的指數(shù)衰減函數(shù)為

式中：Mt和M0分別表示水體在t時刻的質量和初始時刻的質量；e為自然對數(shù)的底；a、b、k分別為通過最小二乘法擬合得到的系數(shù)。

2 1950—2018年七里海潟湖水體交換能力變化特征

1950—2018 年間，由于海平面、徑流以及人類活動的影響，七里海潟湖地貌發(fā)生變化，對水體交換的影響顯著。20世紀50年代至60年代，七里海潟湖區(qū)域主要屬自然演變。20世紀70年代，七里海潟湖潮汐通道維持自然形態(tài)，對入湖河流進行了綜合治理，增強了排洪和泄洪能力，修建了25 km的圍堤，此時潟湖水域面積為8 km2。同期在潟湖潮汐通道內修建擋潮閘，提高七里海調蓄能力，閘底板高程-1.5 m，閘頂高程3.2 m［15］。由于擋潮閘的修建，減弱了潟湖的水動力及水體交換能力，泥沙在潮汐通道內發(fā)生大量淤積，致使潮汐通道向南遷移。1979年，擋潮閘被停止使用，閘門不再閉合，始終保持過水的狀態(tài)［16］。從1984年的衛(wèi)星圖片（圖3b）可以明顯看出，潟湖大部分岸線修建了圍堤，水域面積縮小至7 km2，擋潮閘的存在使潮汐通道口門寬度縮窄至85 m，潮汐通道的形態(tài)發(fā)生改變［17］。潟湖排水能力的降低和沿岸流的作用，使潮汐通道向南側移動。

1984—1991 年（圖3c），七里海潟湖圍墾區(qū)域主要用于水產養(yǎng)殖，水域面積在1990年縮窄至5.2 km2。1991年，疏浚工程使潮汐通道由彎曲變?yōu)轫樦?，潮汐通道外修建有防波堤?991—2000年（圖3d），由于養(yǎng)殖圍墾的進一步加劇，潟湖水域面積縮小至2.6 km2，2001—2008年（圖3e），圍墾工程顯著減少，水域面積僅減少了0.3 km2，縮小至2.3 km2。2008—2018年（圖3f），擋潮閘上半部分拆除，保留了水下結構，并在原始的位置建橋連接潮汐通道兩岸，橋梁采用墩式結構，橋面高4.7 m。同期，潮汐通道東部靠近秦皇島海域部分拓寬，平均寬度由110 m增加至150 m，外海防波堤區(qū)域建立突堤防止海灘侵蝕。

圖3 七里海潟湖衛(wèi)星圖片F(xiàn)ig.3 Satellite images of Qilihai Lagoon

根據(jù)衛(wèi)星圖片及文獻數(shù)據(jù)，為研究長時間尺度下七里海潟湖水體交換能力變化特征，將該時期按照工況建設情況分為5個階段進行模擬，分析不同工況條件下潟湖水體沖刷時間的變化。不同階段的工況條件及模擬時間見表2。

表2 地貌演變模擬時期、參考衛(wèi)星圖片及主要工況條件Tab.2 Different periods of morphological simulations,reference satellite images, and construction conditions

2.1 第一階段水體交換能力變化特征

圖4為1950—2018年七里海潟湖水體沖刷時間及所對應的海平面、徑流和主要人類活動。第一階段（1950—1969年）七里海潟湖主要為自然演變的過程，潮汐通道入湖口門寬度約300 m，對潮流的阻礙作用較弱，平均水體沖刷時間為79 h。據(jù)圖4顯示，該階段水體交換能力最強的時間為1960年，對應水體沖刷時間為46 h，1962年水體交換能力最弱，對應水體沖刷時間為122 h。但1960年的年均徑流最大，為13.75 m3·s-1，1962年的年均徑流最小，為3.39m3·s-1，故徑流量的增大會提高潟湖的水體交換能力。同時，在比較1964年和1965年的水體沖刷能力時發(fā)現(xiàn)，二者年均徑流均為5.77 m3·s-1，但1964年的海平面-0.12 m，較1965年的-0.20 m更高，但水體沖刷時間為82 h較1965年的86 h更小，說明1964年的水體交換能力更強。因此，在徑流條件相同時，海平面升高會促進潟湖水體交換。該階段平均海平面最高的時間為1962年，但其水體交換能力最弱，故可得出徑流對水體交換能力影響較海平面的升高更為顯著。

2.2 第二階段水體交換能力變化特征

第二階段（1970—1979年），在地貌演變的模擬中考慮了潮汐通道口門處擋潮閘的建設（圖3），且該階段擋潮閘處于關閉階段，圖5為典型年份七里海潟湖漲急時刻的流場圖，可以看出在1973年潟湖區(qū)域的流速較第一階段的自然狀態(tài)下顯著降低，說明擋潮閘的建設對潟湖水體流入和流出有明顯的阻礙作用，潟湖水動力大幅減弱，水體交換能力明顯減弱。第二階段平均徑流量為6.45 m3·s-1，與第一階段6.99 m3·s-1的徑流量近似，但平均水體沖刷時間為259 h，是第一階段平均水體沖刷時間的3.28倍，說明擋潮閘的建設顯著降低了水體交換能力。同時，根據(jù)圖4可以看出，1970年徑流為6.65 m3·s-1，較1969年徑流（3.68 m3·s-1）有大幅增加，但1970年的水體沖刷時間（247 h）較1969年（104 h）卻增加了138%，水體交換能力總體呈現(xiàn)減弱趨勢，說明擋潮閘的建設對水體交換能力的影響較徑流更為明顯。1973年的水體沖刷時間最大，為436 h，但年徑流最小，為2.75 m3·s-11；1979年的水體沖刷時間最小，為189 h，但年徑流最大，為10.17 m3·s-11，說明徑流增加會增強潟湖的水體交換能力，此結論與第一階段一致。該階段的平均海平面變化范圍為-0.18～-0.08 m，但是海平面的變化與水體沖刷時間的變化趨勢相關性較差。因此，該階段由于擋潮閘的建設對水體流入與流出的阻礙作用較大，對水體交換能力的影響遠大于徑流和海平面變化的影響。

圖5 1960、1973以及2010年七里海潟湖流速分布Fig.5 Distribution of current velocity in Qilihai Lagoon in 1960,1973, and 2010

2.3 第三階段水體交換能力變化特征

第三階段（1980—1989年），由于環(huán)境及泥沙淤積問題，擋潮閘停止關閉，始終保持開啟狀態(tài)，且潟湖周圍修建圍堤使得水域面積從8.0 km2縮小至5.2 km2，平均水體沖刷時間由第二階段的259 h降低為175 h，水體交換能力增強32%。根據(jù)圖4，1980年徑流和1979年徑流接近，但水體沖刷時間從189 h降低至93 h，水體交換能力增加超過50%，說明擋潮閘的開啟能夠顯著改善潟湖的水體交換能力。與自然狀態(tài)（第一階段）進行對比時發(fā)現(xiàn)，1952年和1988年徑流均為4.58 m3·s-1，但水體沖刷時間由1952年的88 h增加至141 h，說明圍墾工程使?jié)暫乃w交換能力較自然條件下有明顯的減弱。為了研究圍墾面積與水體交換能力的關系，參考陳紅霞等［18］給出的納潮量計算公式（公式（5）），進行了該階段單位面積納潮量的計算。

式中：W為納潮量；S1、S2分別為平均高、低潮潮位的水域面積；h1、h2分別為S1、S2所對應的潮高。

計算得到該階段潟湖單位面積納潮量由0.526 m3·m-2降低為0.519 m3·m-2，降低幅度僅為1%。納潮量與圍墾工程所造成的水域面積減小及海平面變化有關。據(jù)圖4可知，該階段的平均海平面變化微弱，納潮量的改變主要來源于圍墾面積的變化，也就是該階段潟湖圍墾對水體交換能力的影響有限，而擋潮閘的開啟是水體交換能力改善的關鍵。該階段的平均水體沖刷時間的變化范圍為93～228 h，其中1980年水體交換能力最強，對應的年徑流最大，為9.99 m3·s-1；1983年水體交換能力最弱，對應的年徑流最小，為1.33 m3·s-1，故在相同的工況條件下，徑流越大，潟湖水體交換能力越強，此結論與第一階段一致。

2.4 第四階段水體交換能力變化特征

第四階段（1990—2007年）七里海潟湖內潮汐通道由彎曲變?yōu)轫樦保瑵暫M一步圍墾，水域面積減小了2.9 km2，該階段平均徑流為2.55 m3·s-1，平均水體沖刷時間由第三階段的175 h降低至152 h，水體交換能力提升13%。據(jù)圖4可以看出，雖然1990年的徑流為2.28 m3·s-1， 略小于1989年的2.86 m3·s-1，但水體沖刷時間從179 h降低至137 h，水體交換能力增加23%，說明潮汐通道的變化顯著增強潟湖水體交換能力。同時，由于該階段繼續(xù)進行潟湖圍墾，海平面微弱上升，根據(jù)公式（5）計算該階段單位面積納潮量由0.519 m3·m-2降低至0.474 m3·m-2，降低幅度為9%，但水體交換能力并未隨之減弱，說明潮汐通道的順直較圍墾工程來說對水體交換能力的改善更為顯著。

綜合第三和第四兩個階段，潟湖大規(guī)模的圍墾均使得潟湖單位面積納潮量減小，減小幅度分別為1%和9%，對水體交換能力的影響不及潮汐通道工況變化的影響，因此潮汐通道工況的變化是該階段影響水體交換能力的關鍵因素。

2.5 第五階段水體交換能力變化特征

第五階段（2008—2018年）七里海潟湖潮汐通道拓寬，平均寬度由110 m增加至150 m，并在入湖口門處建橋，其中4個直徑為2 m的圓形橋墩與水體接觸，由于橋墩直徑較潟湖口門寬度很小，故建橋對動力的影響較弱。根據(jù)圖5c可知，該階段水體的流入與流出有明顯改善，水體交換能力增強。此階段的徑流條件趨于穩(wěn)定，但海平面有了明顯的上升，增加約0.15 m（圖4）。由于潮汐通道的拓寬以及海平面的變化，潟湖的單位面積納潮量為0.517 m3·m-2，較第四階段的0.474 m3·m-2提升了9%，平均水體沖刷時間由第四階段的152 h減小至126 h，水體交換能力提升17%。根據(jù)第一階段得出的結論，海平面的變化對潟湖水體交換能力的影響較弱，故該階段沖刷時間的降低主要來自于納潮量的變化，即潮汐通道的拓寬提升了潟湖的納潮量，有利于潟湖發(fā)生水體交換。對比2001年和2013年的情況，徑流條件相同，均為0.82 m3·s-1，但2013年的水體沖刷時間較2001年更小，水體交換能力更強，也就是潮汐通道拓寬增強了潟湖內的水體交換。

3 結論

本研究基于1950—2018年月徑流量、海平面變化與七里海潟湖地貌演變的模擬結果，利用MIKE21建立七里海潟湖水動力模型，根據(jù)水體沖刷時間指標分析長時間尺度下潟湖的水體交換能力。研究結果表明：

在自然條件下，徑流增加與海平面升高均可增強潟湖的水體交換能力，但徑流對水體交換能力的影響更為顯著；在相同工況的條件下，徑流越大，水體交換能力越強。

受擋潮閘建設與潟湖圍墾工程的影響，七里海潟湖的水體交換能力發(fā)生變化。擋潮閘建設阻礙了潟湖的流入與流出，降低了潟湖水體交換能力，但當其開啟時水體交換能力得到大幅提升，因而擋潮閘的啟閉成為影響潟湖水體交換能力的關鍵因素；同時，潟湖圍墾降低了潟湖的納潮量，但由于納潮量變化幅度較小，減少水體交換能力有限。

潮汐通道的順直和拓寬使得潟湖水體交換能力進一步增強。其中，第四階段由于圍墾工程所帶來的納潮量變化并未對潟湖水體交換能力造成顯著影響，因而潮汐通道由彎曲變?yōu)轫樦焙屯貙捠窃鰪娝w交換能力的主要因素。

作者貢獻聲明：

匡翠萍：方法提出，論文撰寫與修改。

王 潔：數(shù)值模擬，論文撰寫。

董智超：數(shù)值模擬，圖像繪制。

劉會欣：現(xiàn)場工作。

朱 磊：現(xiàn)場工作。