長江口橫沙淺灘挖入式港池對流場的影響分析Ⅰ：數(shù)值模型和驗(yàn)證

葛建忠， 郭文云， 丁平興

（華東師范大學(xué) 河口海岸學(xué)國家重點(diǎn)實(shí)驗(yàn)室，上海 200062）

0 引 言

長江河口是中等強(qiáng)度的潮汐河口，三級分汊、四口入海；口外為正規(guī)半日潮，口內(nèi)為非正規(guī)半日淺海潮.在徑流和潮流兩股強(qiáng)勁動力作用下，河口段河床沖淤多變，四條入海汊道皆存在攔門沙，均位于口門附近.由于攔門沙的存在，河口航道的自然水深往往比內(nèi)河和海洋較淺，河口攔門沙區(qū)域淺灘灘頂水深常年保持在6.0m左右；這樣的水深再加上灘長、坡緩、變化復(fù)雜等特點(diǎn)，阻礙通海航運(yùn)的發(fā)展，不滿足大型集裝箱船舶的吃水要求［1］.在長江充足的泥沙供給的背景下，長江口及鄰近海域泥沙輸運(yùn)過程顯著，從而形成該區(qū)域水深較淺的整體態(tài)勢；深水岸線不足，難以承接吃水較深的新型集裝箱和石油礦石貨船.為滿足長江流域中長期航運(yùn)發(fā)展需求，提出了橫沙淺灘挖入式港池方案.該方案在橫沙淺灘－5 m等深線內(nèi)建設(shè)挖入式港池，有效避免長江口攔門沙的影響，能夠保證港池內(nèi)具有良好的水深條件以接納新型大噸位船舶［2］.

挖入式港池方案規(guī)劃沿橫沙淺灘－5 m等深線位置建設(shè)圈圍堤，并在南側(cè)與長江口深水航道北導(dǎo)堤連接.此規(guī)劃方案必然對長江口門附近區(qū)域的北港、北槽、南槽的水動力場產(chǎn)生一定的影響.因此，需要對規(guī)劃方案實(shí)施后對周邊流場的可能影響進(jìn)行分析和評估——一般可采用物理模型和數(shù)值模型兩種方法——本文擬從數(shù)值模型角度進(jìn)行分析.

在長江口區(qū)域的有關(guān)研究中，數(shù)值模型的應(yīng)用較為廣泛.顧偉浩等模擬了北槽深水航道對鹽水入侵可能產(chǎn)生的影響［3］.匡翠萍等成功應(yīng)用三維數(shù)值模式模擬了長江口攔門沙地區(qū)水流、鹽度［4］.朱建榮、吳輝等基于改進(jìn)的ECOM模式，進(jìn)行了理想河口鹽水入侵動力過程的模擬分析［5－7］.在此基礎(chǔ)上，吳輝、朱建榮等建立了一個高分辨力的三維水動力、鹽度數(shù)值模式，對長江口水動力模型進(jìn)行了系統(tǒng)性的水動力及鹽度的變化模擬.胡克林等采用Delft3D－Flow也對長江口北槽及鄰近九段沙區(qū)域的水動力和地貌特征進(jìn)行了模擬分析［8］.以上的數(shù)值模型主要采用四邊形曲線網(wǎng)格，而曲線正交網(wǎng)格數(shù)值模型難以良好地?cái)M合圈圍堤、港池內(nèi)碼頭岸線及北槽導(dǎo)堤丁壩等不規(guī)則幾何結(jié)構(gòu)，同時(shí)也難以反映導(dǎo)堤丁壩結(jié)構(gòu)反映對水流的垂向阻隔作用［9］.

橫沙淺灘北接北港航道，南靠北槽深水航道，且與深水航道北導(dǎo)堤結(jié)合，規(guī)劃港池內(nèi)碼頭岸線呈不規(guī)劃走向.而無結(jié)構(gòu)三角形網(wǎng)格可以很好地?cái)M合這些復(fù)雜的平面結(jié)構(gòu)及水下地形變化［10］，保證在所關(guān)注的規(guī)劃工程區(qū)域提供較高的空間分辨率，并能兼顧北槽深水航道導(dǎo)堤丁壩的半淹沒結(jié)構(gòu)［9］.

FVCOM（Finite Volume Coastal and Ocean Model）模型是美國麻州大學(xué)建立的一個三角形網(wǎng)格、有限體積、三維原始方程組的海洋模式，其最大特色是結(jié)合了有限元法易擬合邊界、局部加密與有限差分便于離散計(jì)算海洋原始方程組的優(yōu)點(diǎn)［10－12］.考慮到長江口深水航道工程的實(shí)際情況，導(dǎo)堤、丁壩高于平均海平面約0.2 m，高潮時(shí)被水淹沒，低潮時(shí)露出水面.大部分海洋模型使用緩坡度地形近似處理導(dǎo)堤、丁壩垂向結(jié)構(gòu).葛建忠等通過對FVCOM模式中導(dǎo)堤和丁壩的算法進(jìn)行改進(jìn)，應(yīng)用于航道工程結(jié)構(gòu)引起的水動力及溫鹽變化研究，該模式已成功應(yīng)用于長江口深水航道低鹽水分布和河口環(huán)流模擬［9，13，14］.

基于上述原因，本文擬采用無結(jié)構(gòu)三角形網(wǎng)格海洋數(shù)值模型FVCOM進(jìn)行模擬和分析，建立長江口FVCOM水動力模型，并針對長江口的水文動力環(huán)境進(jìn)行檢驗(yàn)，使之適用于長江口及橫沙淺灘區(qū)域，從而為分析橫沙淺灘挖入式港池方案的可能影響提供合理可靠的模型工具.

1 FVCOM模型的特點(diǎn)與優(yōu)勢

FVCOM是無結(jié)構(gòu)三角形網(wǎng)格架構(gòu)、有限體積、自由表面、三維原始方程海洋數(shù)值模型［10］，其原始方程包含動量方程、質(zhì)量連續(xù)方程以及溫度、鹽度和密度方程，在物理和數(shù)學(xué)上用Mellor－Yamada 2.5階垂向湍流閉合模型及Smagorinsky水平湍流閉合模型對方程組進(jìn)行閉合.在垂向上使用σ坐標(biāo)系或通用垂向s坐標(biāo)系對不規(guī)則底部地形進(jìn)行擬合，在水平上利用無結(jié)構(gòu)三角形網(wǎng)格對水平計(jì)算區(qū)域進(jìn)行空間離散.在數(shù)值計(jì)算上，利用對水平三角形控制體進(jìn)行通量有限體積積分的方式對控制方程進(jìn)行離散求解；該方式結(jié)合了有限元方法的自由幾何擬合特性和有限差分方法的離散結(jié)構(gòu)簡單、計(jì)算高效的特性.利用有限體積積分格式，能更好地保證復(fù)雜幾何結(jié)構(gòu)的河口海灣及海洋計(jì)算中的質(zhì)量、動量、鹽度、溫度及熱量的守恒性.最后，F(xiàn)VCOM模型所采用的無結(jié)構(gòu)三角形網(wǎng)格對復(fù)雜島嶼、岸線及地形具有非常良好的幾何擬合能力.

控制方程組由動量方程、連續(xù)方程、溫度方程、鹽度方程和密度方程組成：

其中，x、y和z分別代表直角坐標(biāo)系中東向、北向和垂向坐標(biāo)，u、v和w 分別是x、y和z三個方向上的流速分量.T是溫度，S是鹽度，ρ是密度，P是壓力，f是科氏力參數(shù)，g是重力加速度，Km是垂向渦動黏性系數(shù)，Kh是熱力垂向渦動摩擦系數(shù)，F(xiàn)u，F(xiàn)v，F(xiàn)θ和Fs代表水平動量、溫度和鹽度擴(kuò)散項(xiàng).總的水深D＝H ＋z，H是平均水深，z是水位變化.

2 區(qū)域模式計(jì)算網(wǎng)格與設(shè)置

長江口模型三角形網(wǎng)格如圖所示，覆蓋長江口、杭州灣、舟山群島、東海內(nèi)陸架及鄰近海域.總?cè)菃卧獢?shù)為94088，節(jié)點(diǎn)數(shù)為49078，該網(wǎng)格在長江口北槽及橫沙淺灘區(qū)域具有較高的空間分辨率，最高達(dá)到200 m左右，網(wǎng)格分辨率最粗位置處于外海開邊界，約14.7 km，能夠較好地刻畫長江口區(qū)域的岸線水深特征，且從局部網(wǎng)格放大圖看出該網(wǎng)格能很好地刻畫北槽深水航道導(dǎo)堤丁壩的平面幾何結(jié)構(gòu)（見圖1）.

本文在長江口內(nèi)及口門海域主要采用1998—2008年期間的實(shí)測地形觀測資料，長江口外水深變化較小，采用海圖數(shù)字化資料.開邊界水位主要采用8個天文分潮M2、S2、K1、O1、N2、K2、P1、Q1合成給出：

其中，ξ為潮位，ξ0為余水位，f為節(jié)點(diǎn)因子，Hi為振幅，ωi為角頻率，gi為遲角，Vi＋ui為訂正角.

圖1 長江口FVCOM三角形網(wǎng)格（上）和北槽橫沙淺灘局部網(wǎng)格（下）Fig.1 The unstructured grid for the Changjiang Estuary（upper panel）and its enlarged view around the Hengsha Shoal（lower panel）

由于實(shí)測資料的缺乏，為給出合理的開邊界條件，首先建立了一個東中國海大區(qū)域的模型來計(jì)算各個天文分潮［14］.本文使用的調(diào)和常數(shù)從大區(qū)域在開邊界出的調(diào)和常數(shù)插值得來.上游徑流給大通站歷年實(shí)測流量（長江水文網(wǎng)，http：／／www.cjh.com.cn／）.流速和水位對外力響應(yīng)較快，初始場一般取為零；溫度、鹽度為慢過程，初始場取自《渤海、黃海、東海海洋圖集》［15］數(shù)字化資料.

數(shù)值計(jì)算采用三維模式進(jìn)行，垂向分層總數(shù)為20層，分層方法為隨地形變化的sigma分層.底部糙率在整個計(jì)算域上進(jìn)行分區(qū)表示，從深水向淺水糙率逐漸增大.采用Mellor－Yamada 2.5階湍流閉合模型，其中，水平混合系數(shù)取0.1，垂向混合系數(shù)取10－4.

3 模型的率定與檢驗(yàn)

本文使用的長江口FVCOM模型已經(jīng)過多年實(shí)測數(shù)據(jù)的率定和檢驗(yàn)［9，13，14，16］，主要驗(yàn)證站位如圖2所示.針對本文研究目的，本文采用三次觀測資料進(jìn)行補(bǔ)充驗(yàn)證，分別進(jìn)行長江口外、口內(nèi)主要汊道和北槽區(qū)域的模型驗(yàn)證.

3.1 長江口外水文環(huán)境驗(yàn)證

數(shù)據(jù)資料源自由國家環(huán)?？偩謫拥摹堕L江口碧海行動計(jì)劃》現(xiàn)場調(diào)查，觀測時(shí)間在2005年7月，圖2為此次觀測長江口口門附近地區(qū)測站位置，觀測內(nèi)容為流速、流向和鹽度.本次測量覆蓋面較廣，基本包括了長江口外代表性位置.

圖2 2005年長江河口水文測站位置Fig.2 The observation sites around the Changjiang Estuary in 2005

模式的徑流邊界采用大通站實(shí)測流量.圖3和圖4分別是流速、流向和鹽度驗(yàn)證圖.從圖中可看出，計(jì)算值和實(shí)測值總體吻合良好，除SH2站外各站平均流速偏差小于10%，流向在漲落急時(shí)平均偏差小于10°，滿足“海岸與河口潮流泥沙模擬技術(shù)規(guī)程JTS／T231－2－2010”對潮流模式的誤差要求，同時(shí)各站鹽度的平均誤差在2%～4%，表明水流鹽度模型可較好的應(yīng)用與長江口外海區(qū)流場及鹽度場計(jì)算.

3.2 長江口內(nèi)主要汊道水文環(huán)境驗(yàn)證

選取崇西、南門和堡鎮(zhèn)作為長江口內(nèi)的潮位驗(yàn)證站點(diǎn)，其結(jié)果如圖5所示.水文模擬平均誤差在5.9%～6.1%之間，可見數(shù)值模型對潮位及余水位的變化過程反映良好.

長江河口地區(qū)鹽淡水交匯劇烈，南支、南槽、北槽和北港等主要汊道的流場、鹽度分布規(guī)律復(fù)雜.從2003年起，逐年在自徐六涇以下至長江口近口外海域進(jìn)行洪季、枯季水文、鹽度等綜合測量（見圖6）.

針對2008年洪季模擬上游開邊界條件由大通站實(shí)測的徑流量資料給出，最小流量在月初為25 000 m3／s，16日以后流量增大較快，在23日達(dá)最大42 266 m3／s.

圖7給出靠近攔門沙河段及近口外6個測站流速、鹽度計(jì)算結(jié)果，從中能看出，模式可很好地反應(yīng)河口攔門沙河段流速、鹽度鋒面變化過程：與實(shí)測資料符合良好.

3.3 北槽水文環(huán)境驗(yàn)證

長江口深水航道工程建設(shè)對北槽流場分布影響顯著.利用2003—2007年洪枯季期間，北槽W3站點(diǎn)流速、鹽度資料，對建立的模式進(jìn)行檢驗(yàn).圖8給出北槽W3站點(diǎn)位置.

圖3 2005年7月表層（上）、中層（中）、下層（下）流速、流向計(jì)算值（實(shí)線）和實(shí)測值（點(diǎn)）比Fig.3 Model－data comparison of velocity and direction at site JS1，JS2，SH1，SH2，SH4 and SH5

圖4 2005年7月各測站表、0.2、0.4、0.6、0.8和底層鹽度計(jì)算值（實(shí)線）和實(shí)測值（點(diǎn)）比較Fig.4 Model－data comparison of six－layer salinity process at sites JS1，SH2，SH2，SH3，SH4 and SH5

2003—2005年北槽二期工程實(shí)施過程中，北槽水流鹽度處于變化過程，同時(shí)長江河口地形變化劇烈.模型水深采用2003—2007年工程階段長江河口實(shí)測地形，導(dǎo)堤和丁壩結(jié)構(gòu)采用FVCOM模型中改進(jìn)算法模塊［17］.

根據(jù)FVCOM模擬結(jié)果，利用模型流速數(shù)據(jù)計(jì)算北槽下斷面大潮期間落潮分流比［18］，與實(shí)測分流比對比.如表1所示，模型計(jì)算的分流比與實(shí)測值基本一致，說明該模型能較好底把握北槽及鄰近區(qū)域主要汊道的水量分配.

圖5 長江口內(nèi)堡鎮(zhèn)、南門及崇西站2010年4月潮位及余水位驗(yàn)證Fig.5 Model－data comparison of hourly water levels at Baozhan，Nanmen，Chongxi stations

從圖9可看出，W3測站落潮流速大于漲潮流速，表層最大流速達(dá)2.5 m／s，底層最大流速可達(dá)1.5 m／s左右.漲潮時(shí)，表層流向略有分散；落潮時(shí)流向比較規(guī)律，約在120°左右.計(jì)算結(jié)果較好的模擬了歷年流速、流向的一個潮周期過程變化特征.

圖10給出了歷年洪枯季，大小潮期間W3測站鹽度過程變化，計(jì)算結(jié)果較好地模擬了鹽度變化過程：2003—2007年四次模擬觀測對比的流速平均誤差分別是6.4%，8.2%，8.9%和7.7%，流向在主流時(shí)段誤差基本在10°以內(nèi)，滿足“海岸與河口潮流泥沙模擬技術(shù)規(guī)程JTS／T231－2－2010”對潮流模式的誤差要求；鹽度模擬誤差為5.7%，也具有良好的精度.

圖7－1 0804、0805、0806和0807測站流速、流向和鹽度計(jì)算值（實(shí)線）和實(shí)測值（點(diǎn)）比較Fig.7－1 Model－data comparison of velocities，directions and salinities at sites 0804，0805，0806 and 0807

圖7－2 0812和0813測站流速、流向和鹽度計(jì)算值（實(shí)線）和實(shí)測值（點(diǎn)）比較Fig.7－2 Model－data comparison of velocities，directions and salinities at sites 0812 and 0813

圖8 北槽W3測站位置Fig.8 Location of W3 site in the North Passage

表1 下斷面落潮分流比驗(yàn)證結(jié)果Tab.1 Comparisons of ebb flow diversion ratio between modeled results and observed data

圖9 W3站點(diǎn)表層（上）、中層（中）、底層（下）流速、流向計(jì)算值（實(shí)線）和實(shí)測值（點(diǎn)）比較Fig.9 Model－data comparisons of velocity and direction of surface，middle and bottom layers at W3 point

圖10 W3站點(diǎn)表層（上）、中層、底層（下）鹽度計(jì)算值（實(shí)線）和實(shí)測值（點(diǎn)）比較Fig.10 Model－data comparison of three－layer salinity at W3 point

4 結(jié) 語

本文建立了無結(jié)構(gòu)三角形網(wǎng)格高分辨率長江口及鄰近區(qū)域的FVCOM模型，綜合考慮多種物理過程，包括徑流、天文潮汐、潮流、溫度、鹽度，并利用長江口主要汊道、口門及口外的水文觀測資料對數(shù)值模型進(jìn)行水文環(huán)境檢驗(yàn).其無規(guī)則三角形網(wǎng)格很好地?cái)M合了長江口門區(qū)域的島嶼岸線及復(fù)雜地形；采用三維模式，模型能合理地模擬從口內(nèi)到口外的水動力變化情況，并且在橫沙淺灘區(qū)域及鄰近的北槽和北港水道具有良好的精度；其主要站位的流速和流向模擬滿足“海岸與河口潮流泥沙模擬技術(shù)規(guī)程（JTS／T231－2－2010）”對潮流模式的誤差要求；可以作為下一步橫沙淺灘挖入式港池規(guī)劃方案對周邊流場影響分析的數(shù)值模型工具，以計(jì)算不同的工程方案布置對北港，特別是北槽深水航道的直接和間接影響，以及工程后港池內(nèi)及外航道的水動力變化情況.

［1］ 陳吉余，沈煥庭，惲才興.長江河口動力過程和地貌演變［M］.上海：上?？茖W(xué)技術(shù)出版社，1988.

［2］ 交通部第三航務(wù)工程勘察設(shè)計(jì)研究院.上海國際航運(yùn)中心橫沙淺灘挖入式港池規(guī)劃方案研究報(bào)告［R］.上海：交通部第三航務(wù)工程勘察設(shè)計(jì)研究院，2012.

［3］ 顧偉浩.長江口南北槽咸水入侵—兼談開挖北槽為深水航槽［J］.水運(yùn)工程，1985（2）：1－3.

［4］ 匡翠萍.長江口鹽水入侵三維數(shù)值模擬［J］.河海大學(xué)學(xué)報(bào)1997，25（54）：54－60.

［5］ 朱建榮，胡松.河口形狀對河口環(huán)流和鹽水入侵的影響［J］.華東師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2003（2）：68－73.

［6］ 朱建榮，胡松，傅德健，等.河口環(huán)流和鹽水入侵I：模式及控制數(shù)值試驗(yàn)［J］.青島海洋大學(xué)學(xué)報(bào)，2003，33（2）：180－184.

［7］ 吳輝.長江河口鹽水入侵研究—北支倒灌、深水航道工程和冬季季風(fēng)的影響［D］.上海：華東師范大學(xué)，2006.

［8］ HU K，DING P.The Effect of Deep Waterway Construction on Hydrodynamics and Salinities in Yangtze Estuary，China［J］.Journal of Coastal Research，2009：961－965.

［9］ GE J，CHEN C，QI J，et al.A dike－groyne algorithm in a terrain－following coordinate ocean model（FVCOM）：Development，validation and application［J］.Ocean Modelling，2012，47（C）：26－40.doi：10.1016／j.ocemod.2012.01.006

［10］ CHEN C，LIU H，BEARDSLEY R C.An unstructured grid，finite－volume，three－dimensional，primitive equations ocean model：application to coastal ocean and estuaries［J］.Journal of atmospheric and oceanic technology，2003，20（1）：159－186.

［11］ CHEN C，COWLES G，BEARDSLEY R C.An unstructured grid，finite－volume coastal ocean model：FVCOM user manual［R］.2nd ed.SMAST／UMASSD Technical Report－06－0602，2006：45.

［12］ CHEN C，BEARDSLEY R C，COWLES G.An unstructured grid，finite－volume coastal ocean model（FVCOM）system［J］.Oceanography：Special Issue，2006，19（1）：78－89.

［13］ GE J，DING P，CHEN C，et al.Low－salinity plume in the Changjiang and adjacent coastal regions：a model－data comparison［C］.Proceedings of the 31st International Conference of Coastal Engineering Hamburg，Germany，2008：4471－4481.

［14］ GE J，DING P，CHEN C.Imapcts of deep waterway project on local circulations and salinity in the Changjiang estuary，China［C］.Proceedings of the 32st International Conference of Coastal Engineering Shanghai，China，2010.

［15］ 海洋圖集編委會.渤海、黃海、東海海洋圖集［M］.北京：海洋出版社，1993.

［16］ 范中亞，葛建忠，丁平興，等.長江口深水航道工程對北槽鹽度分布的影響［J］.華東師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2012（4）：181－189.

［17］ 葛建忠.Hulti－scale FVCOM model system for the East China Sea and Changjiang Estuary and its applications［D］.上海：華東師范大學(xué)，2010.

［18］ 高敏，范期錦，談澤煒，等.對長江口北槽分流比的分析研究［J］.水運(yùn)工程，2009（5）：82－86.