朱子晨胡澤建劉建強張莞君張永強黃炳智

(自然資源部 第一海洋研究所,山東 青島266061)

海灣是海洋資源開發(fā)利用的重要區(qū)域,由于其優(yōu)良的自然環(huán)境和開發(fā)條件,海灣地區(qū)往往是人口聚集、經濟發(fā)展的中心區(qū)域[1]。對多數海灣地區(qū)而言,由于具有較好的波浪掩護條件,潮流常常是泥沙運動的主要動力因素。而海灣內的泥沙運動,影響著地形地貌及灣內港口、能源工程的運營和維護[2-3]。懸浮泥沙能夠吸附重金屬等化學物質,對光具有吸收作用,又對海灣的生態(tài)環(huán)境具有顯著影響[4]。海水中的懸浮泥沙來源于再懸浮、水平輸運等方式,其組成中包含多種粒徑成分。由于不同粒徑的泥沙,其起動、沉降等特征均不相同,各成分的百分含量在一個潮周期內具有較為明顯的變化。王愛軍等[5]2005年對長江口懸沙粒度與濃度關系的研究顯示,特征時段內懸沙中值粒徑與濃度存在指數關系；張一乙等[6]2016年對長江口懸沙組分對再懸浮過程響應的研究顯示,砂質組分是懸沙中對再懸浮過程響應的敏感組分；陳語等[7]2016年對長江口渾濁帶懸沙粒度分布的研究顯示,懸沙粒度的垂向差異在大潮期十分顯著,而小潮期粒徑垂向梯度較低；徐海和李伯根[8]2009年對椒江河口懸沙分選機制研究的結果表明,懸沙粒徑分選主要受物質來源、潮流動力作用下底質再懸浮、絮凝沉降三個因素影響。這些研究表明,我國河口、海岸帶地區(qū)懸沙組成的變化規(guī)律受到廣泛關注。泥沙組分的運動、變化規(guī)律與海洋開發(fā)及環(huán)境保護息息相關,如某些海域開發(fā)后,人工構筑物附近出現的泥化現象等。為探究潮流作用下懸浮泥沙組成成分對動力因素的響應機制,本文采用觀測與模型相結合的方法,選擇潮差較大、波浪掩護條件較好的湄洲灣羅嶼支水道為目標區(qū)域,對半日潮周期的時間尺度內懸沙組成的變化特征進行探討。

1 研究方法

1.1 水動力觀測與泥沙取樣分析

湄洲灣位于福建省莆田、泉州兩市之間,秀嶼站1978-01—1980-12潮汐觀測顯示,該地區(qū)平均潮差5.12 m,最大潮差7.59 m,潮差較大,且海灣狹長,波浪掩護條件較好。本研究選擇在湄洲灣羅嶼東南側支水道設置測站,同步觀測水深(潮汐)、潮流、懸沙,測站處水深9.77 m(以平均海平面計),觀測時段為2007-08-14T19:00—08-15T20:00(大潮期),涵蓋2個完整的半日潮周期,并對底質取樣。研究海域及測站位置見圖1。潮位采用水位計觀測,海流、懸沙采用5層觀測,潮流觀測采用Aanderaa RCM9型海流計,懸沙使用采水器采集后,在實驗室分析,并對漲落潮中間時刻與高低潮時刻4個特征時刻的懸沙粒度組成進行分析。

圖1 研究海域及測站位置Fig.1 Study area and the location of observation site

1.2 一維水深平均懸沙模型

1)控制方程與假設

湄洲灣以M2,S2為主要潮流分量,本次觀測時段內M2分潮流占比59%,S2分潮流占比18%,其他分潮流占比均小于9%,M4分潮流僅占比4%①國家海洋局第一海洋研究所.福建湄州灣電廠二期擴建工程可行性研究工程海域原型觀測專題報告,2007.,因此本文模型水動力因素主要考慮M2,S2分潮,并進行適當簡化。觀測期間,測站處潮流流向幾乎與水道走向平行,往復流態(tài)十分明顯,漲潮流向40°,落潮流向217°,因此,本文使用一維水深平均懸沙控制方程來仿真所觀測到的懸沙序列:

式中,c表示懸沙濃度；t表示時間；H表示水深；U表示海流流速,羅嶼水道往復流特征顯著,近似假設潮流為一維,以漲潮流方向為正；x為水平距離,正方向與海流方向一致；Qr與Qd則分別代表沉積物再懸浮量和沉降量。〈〉代表沿水深方向平均。由于研究測站處水動力條件較強,本次觀測到的漲急流速達到68 cm/s,落急流速73 cm/s,水平擴散項量級遠小于水平輸沙項,忽略水平擴散項[9-10]。

式(1)中,平流項采用下列分解式表達:

式中,下標“d”代表該量與其垂向平均值的差。式(2)中等號右側第二項代表垂向流速與懸沙濃度剖面特征之間的非線性相互作用,為了對模型進行簡化,將該非線性項忽略[11-12]。最終,平流項可以表示為

式(1)中,再懸浮項由底層海流剪切海床產生,本文采用水深平均模型仿真懸沙變化,因此,為了表達再懸浮項Qr與垂向平均流的關系,采用如下方法計算[13]:

式(4)中,M表示泥沙單位起動量(單位面積、單位時間)；ρ表示海水密度；CD表示拖曳系數；τcr表示臨界起動切應力。假設M,ρ,CD,τcr均為常數,因此影響沉積物再懸浮量的參數B簡化為常數[14]。

懸沙沉降量與沉降速度成正比,并受海水底層懸沙和水深的影響,沉降項采用下式表達[3,14]:

式中,D=Kdws,ws表示泥沙沉降速度；cb表示海水底層懸沙質量濃度。最終,沉降項由垂向平均懸沙質量濃度、水深、沉降速度和沉降參數D表示。根據一維連續(xù)方程:

結合式(6)與式(3),并將式(1)～(6)合并,一維垂向平均懸沙模型控制方程表示為

式(7)表示,懸沙質量濃度的時間序列主要受3個物理過程控制,按等式右側順序依次是平流項、再懸浮項與沉積項。上述3項可使用簡諧波表示[13]。式(7)中,垂向平均流速、總水深(水深與潮位之和)均可利用觀測結果,通過調和分析表達,并可將式(7)以三角傅里葉變換的形式表示為

式中,Ai、ωi與ψi為各個傅里葉分量的振幅、角速度與初相位。將式(8)對時間積分后,懸浮泥沙濃度表示為

式中,u0表示余流；u1、ω1、φ1分別表示 M2分潮流的振幅、角速度與初相位；u2,ω2,φ2分別表示S2分潮流的振幅、角速度與初相位。相比于u0和u2,u1在量級上具有絕對優(yōu)勢,因此做如下變換:

相近地,水深可采用下式表示:

式中,h0為平均水深(以平均海平面計),h1,φ1分別為M2分潮的振幅和初相位；h2,φ2分別為S2分潮的振幅和初相位。由于測站平均水深是潮位振幅的4倍,因此做如下變換:

則式(7)再懸浮項中水深的倒數可表示為

式(7)中平流項受懸沙濃度水平梯度控制,該量十分復雜,且難以觀測,將該梯度進行簡化,設水平濃度梯度為常數[3,9,15],懸沙質量濃度水平梯度表示為

結合式(7)、式(9)、式(11)和工(14),并將振幅中包含u0/u1,u2/u1,h1/h0與h2/h0的高階項(階數＞1)忽略[14],最終將懸沙分解為12個主要的傅里葉分量,其時間序列可表示為

表1 三角傅里葉分量參數Table 1 Parameters of Fourier components

各傅里葉分量中(表1),E1代表在整個潮周期內懸沙的平均水平,不隨時間變化,該分量由M2分潮流與余流共同作用下的水平輸沙、M2分潮流剪切海底引起的再懸浮兩部分組成,是這2個物理過程的疊加；E2也是2個物理過程的疊加,包括M2分潮流作用下的水平輸沙、M2分潮流與余流共同作用剪切海底引起的再懸?。籈3,E4和E8三個分量同源,是M2分潮流剪切海底與M2分潮引起的水深變化共同作用產生的再懸浮項,通過三角變換分解為3個分量；E5是S2分潮流作用下的水平輸沙項；E6,E11和E12三個分量同源,是M2分潮流剪切海底與S2分潮引起的水深變化共同作用產生的再懸浮項,通過三角變換分解為3個分量；E7分量是由M2分潮流剪切海底產生的再懸浮項,且不受其他動力因素控制；E9,E10兩項同源,是M2分潮流與S2分潮流共同作用剪切海底產生的再懸浮項,通過三角變換分解為2個分量,且E10分量的角速度為M2與S2分潮角速度之差,具有半月潮周期,說明在M2與S2分潮的共同作用下,懸沙序列能夠體現出半月潮周期變化。

2)參數設置

測站處底質粒徑分布見圖2a,低潮、漲急、高潮、落急四個特征時刻懸沙的粒徑組成見圖2b。測站處底質包含淤泥、粉砂與細砂各個成分,最大粒徑2 mm,但懸沙成分僅包含淤泥、粉砂兩種成分,最大粒徑0.063 mm。觀測結果說明細砂并未起動,因此模型將懸沙分為淤泥、粉砂兩種組分進行模擬,以粒徑0.016 mm代表粉砂,粒徑0.002 mm代表淤泥組分,懸沙濃度為兩者之和。

模型自由參數B,D,k與泥沙特性及海水特性等有關,如泥沙臨界起動切應力、沉降速度等,由于其難以測定,在同類的模型研究中通常采用觀測結果進行率定[9,16-17]。在率定前,部分參數可根據現有研究成果限定合理的取值范圍。

圖2 底質與懸沙粒徑組成Fig.2 Grain diameter compositions of sediment and suspension

k取值范圍在-2.0×10-5～2.0×10-5kg/m4[14],本研究取k=-1.0×10-5kg/m4。單位起動量M取值范圍在2.0×10-5～3.5×10-3kg·m-2·s-1[18],粉砂組分取3.5×10-3kg·m-2·s-1,淤泥組分取2.48×10-3kg·m-2·s-1。粉砂臨界起動切應力τcr在0.05～0.13 N/m2[19],本文研究取0.05 N/m2；淤泥臨界起動切應力τcr在0.05～0.2 N/m2[19-21],本研究取τcr=0.14 N/m2。CD取值5.0×10-3。海水密度ρ取1 025.0 kg/m3。粒徑為0.01 mm的泥沙沉速大于0.07 mm/s[22],絮凝體的沉降速度小于3 mm/s[23]。對于各組分沉降速度ws,粉砂組分取1.06 mm/s[19]；淤泥組分沉降速度在0.003～0.2 mm/s[24],本文取0.2 mm/s。沉降參數Kd則較為復雜,影響因素包括垂向平均懸沙與近底層懸沙之比、水深等,經過率定,參數Kd粉砂組分取值1 800,淤泥組分取值4 680。由式(5)計算得到,粉砂組分D=0.195,淤泥組分D=0.096。

2 結果與驗證

本期觀測在大潮期內進行,觀測到的最大潮差6.32 m,漲潮過程最大流速0.68 m/s,落潮過程最大流速0.73 m/s,懸沙質量濃度最高達到94.5 mg/L,最低質量濃度39.3 mg/L。對應時段內,經調和分析計算到的最大潮差6.01 m,漲潮過程最大流速0.73 m/s,落潮過程最大流速0.60 m/s,模擬得到的最高懸沙質量濃度81.9 mg/L,最低質量濃度42.0 mg/L。

將計算結果與觀測資料比對驗證,水深、潮流及懸沙質量濃度驗證結果見圖3a～3c,粉砂、淤泥組分占比驗證見圖3d,驗證結果顯示,計算結果與觀測資料基本吻合,能夠反映懸沙質量濃度隨時間的變化趨勢,粉砂、淤泥組分隨時間變化的計算結果也與觀測結果吻合。本文由于對模型進行簡化,潮流僅考慮了M2,S2兩個分量及余流,因此,潮汐、潮流計算結果與觀測結果存在誤差,且水動力計算誤差可傳遞至懸沙計算,最終引發(fā)懸沙計算誤差。除上述原因外,懸沙計算誤差的產生原因還包括:1)雖然湄洲灣內波浪掩護條件較好,但依然受到波浪影響,因此懸沙的部分高頻率的變化特征沒有被仿真出來；2)測站處底質與懸沙的組成成分復雜,含有粉砂、淤泥的各個粒徑成分,本文采用2種粒徑的泥沙來代替成分復雜的混合物,因此存在計算誤差。經準調和分析后,垂向平均流速及水深近似以式(10)及式(12)表示為

圖3 計算結果驗證Fig.3 Validadtion of model simulations

由于實際懸沙在海水中的運動是三維的,而本文采用一維水深平均模型進行簡化,因此模型本身存在誤差。觀測期間,漲落潮流主流向相差177°,在平面上,本文模型簡化為一維,這是造成模型從動力方面產生誤差的一個方面。另一方面,本文選擇垂向平均懸沙濃度作為研究指標,因此忽略了懸沙的垂向結構,沒有能夠反映海流作用下懸沙的垂向輸運和垂向擴散兩個物理過程,這會造成模型計算結果與實際觀測結果之間存在濃度和相位兩方面的誤差。從模型驗證結果來看,簡化模型基本反映了實際情況,說明簡化后的模型抓住了主要的物理過程。

圖4 潮流矢量圖Fig.4 Tidal current vectors

圖5 懸沙組分質量分數對潮流的響應Fig.5 Responses of grain diameter composition to current

3 討 論

3.1 各組分起動、沉降能力比較

圖5所示為模型計算得到的粉砂、淤泥組分在懸沙中所占質量分數隨潮流變化的過程,由于模型沒有考慮波浪作用,剪切作用強弱僅與潮流流速有關,因此對流速取絕對值進行對照分析。在所模擬的2個半日潮周期內,粉砂所占質量分數在62.0%～69.4%,淤泥占比在30.6%～38.0%；在觀測到的4個特征時刻,粉砂占比在62%～71%,淤泥占比在29%～38%。模擬與觀測結果均顯示,粉砂組分占優(yōu)勢。而底質分析結果顯示,細砂占比21.6%,粉砂占比40.6%,淤泥占比37.8%,底質中粉砂與淤泥組分較為接近。從沉降角度考慮,粉砂組分由于粒徑較粗,沉降速度高于淤泥。本文模型中沉降參數D經率定后,粉砂組分D=0.195,淤泥組分D=0.096,粉砂沉降能力明顯高于淤泥。粉砂組分在沉降能力明顯高于淤泥組分、且在底質中占比與淤泥組分接近的條件下,在懸沙中的占比明顯高于淤泥,這說明,粉砂組分的起動能力明顯高于淤泥組分,淤泥組分起動條件較高與黏性泥沙黏結力等特性有關[19]。由于粉砂組分的起動和沉降能力均高于淤泥組分,因此,懸沙中粉砂百分含量的極大值出現在動力條件較強時刻附近,約在漲落急時刻1.5 h后,淤泥質量分數的極大值則出現在動力較弱的漲落憩時刻1.5 h后,這說明,漲落急或漲落憩等潮流特征時刻所觀測到的懸沙組成,并不能代表懸沙組成的極值,但采集到的懸沙樣本基本包含了本地懸浮體中各粒徑成分。

3.2 懸沙組分傅里葉分量

本文模型分解得到的12個傅里葉分量(表1)中,E1分量代表懸沙序列的平均水平,不隨時間變化；E10分量的角速度為M2分潮與S2分潮角速度之差,具有半月潮周期,說明在M2與S2分潮共同作用下,懸沙能夠表現出大小潮變化趨勢,這是因為M2與S2分潮的角速度接近卻不重合,但由于本文研究沒有觀測到完整的半月潮周期,因此主要討論E2～E9,E11,E12分量。

圖6所示為M2,S2分潮流及余流共同作用下,對粉砂、淤泥兩種組分進行傅里葉分析得到的主要諧波分量,結果顯示,2種懸沙組分均以E7為主要傅里葉分量,E7具有M2分潮2倍的角速度(表1),呈1/4日潮周期,是在M2分潮流的剪切作用下產生,揭示了半日潮海區(qū)在潮流作用下,懸沙呈1/4日潮周期的原因。其次,主要傅里葉分量還包括E2,E9,E3分量。E2由水平輸沙、余流、M2分潮流共同作用產生,具有與M2分潮相同的角速度及半日潮周期；E9由M2,S2分潮共同作用產生,角速度為M2,S2分潮角速度之和,具有1/4日潮的周期；E3主要由M2分潮剪切海床產生,具有M2分潮的角速度及半日潮周期?？倯腋∧嗌秤筛鞣至刊B加而成,而上述4個主要傅里葉分量分別具有1/4日潮周期與半日潮周期,因而懸沙濃度在一個潮周期內呈現出不對稱特征,即一個半日潮周期內兩次峰值具有明顯差異(圖3c),這與觀測結果吻合。

圖6 傅里葉分量的時間序列Fig.6 Time series of Fourier components

對比圖6a與圖5b,粉砂組分各分量的振幅均顯著高于淤泥組分,這是由于淤泥組分起動和沉降能力均弱于粉砂組分造成的。觀測懸沙總質量濃度波動范圍在39～95 mg/L,特征時刻觀測到的粉砂組分波動范圍在29～41 mg/L,淤泥組分波動范圍在17～20 mg/L,觀測結果也顯示,粉砂組分的振幅明顯高于淤泥組分。本次模型計算結果顯示,總懸沙波動范圍在42～82 mg/L,粉砂組分波動范圍在26～57 mg/L,淤泥組分波動范圍在16～25 mg/L,與觀測結果基本吻合,趨勢一致。作為主要傅里葉分量,粉砂組分的E7分量振幅為6.1 mg/L,而淤泥組分的E7分量振幅為1.5 mg/L；此外,粉砂組分的E2,E9,E3分量振幅分別為4.9,3.7和2.7 mg/L,淤泥組分的E2,E9,E3分量振幅分別為1.2,0.9和0.7 mg/L,均說明粉砂組分的波動幅度高于淤泥組分。

3.3 余流對懸沙組分的影響

為討論余流對海區(qū)懸沙組分變化特征的影響,假定其他因素固定,對比真實余流(u0=0.06 m/s,最適率定值)、余流消失(u0=0)、余流增強一倍(u0=0.12 m/s)與余流反向(u0=-0.06 m/s)四種情況下懸沙組分百分含量的變化特征,如圖7所示。對比結果顯示,由于余流方向與漲潮流方向相同,當余流消失時,漲潮過程粉砂組分所占百分比將下降,落潮過程粉砂組分將升高,而淤泥組分則相反,且由于余流對主要半日潮傅里葉分量(E2)的影響,如果余流消失,半日潮周期內懸沙相鄰的兩個波峰之間的不對稱性(漲潮過程波峰高于落潮過程波峰)減弱；當余流增強一倍時,半日潮周期內懸沙相鄰波峰之間的不對稱性將增強,漲潮過程粉砂組分占比將上升,淤泥組分將下降；當余流反向時,半日潮周期內懸沙相鄰波峰之間的不對稱性發(fā)生反轉,變?yōu)闈q潮過程波峰低于落潮過程波峰,漲潮過程粉砂組分將下降,淤泥組分將提高。

圖7 余流對懸沙組分的影響Fig.7 Influence of residual current on suspended sediment components

4 結 論

以泥沙的起動、沉降和水平輸運為主要物理過程,根據所研究海域動力環(huán)境,以M2,S2為主要潮汐、潮流分量,建立了一維水深平均懸沙模型,采用傅里葉分析方法,將懸沙時間序列分解為12個主要傅里葉分量,模擬了湄洲灣海域2個完整半日潮周期內懸沙中粉砂、淤泥兩種組分的運動過程。模型能夠與水動力、懸沙同步觀測吻合,并能仿真出各粒級組分對潮流的響應規(guī)律。

模型分解得到的主要傅里葉分量具有M2分潮兩倍的角速度與1/4日潮周期,該分量粉砂組分振幅6.1 mg/L,淤泥組分振幅1.5 mg/L；次主要傅里葉分量具有M2分潮的角速度,振幅受水平輸沙、余流、M2分潮流共同影響,粉砂組分振幅4.9 mg/L,淤泥組分振幅1.2 mg/L；振幅第三的分量由M2,S2分潮共同作用產生,角速度為M2,S2分潮角速度之和,具有1/4日潮的周期,粉砂組分振幅3.7 mg/L,淤泥組分振幅0.9 mg/L。觀測與模型主要分量均顯示,粉砂組分振幅高于淤泥組分,這是由于粉砂組分單位起動能力強、沉降速度高造成的,且淤泥組分由于黏結力等因素,起動條件較高。

本文測站位置余流與漲潮流方向一致,余流的存在致使?jié)q潮過程中粉砂組分含量所占百分比上升,而落潮過程下降,淤泥組分則相反。此外,傅里葉分析說明,具有半日潮周期的傅里葉分量中振幅最大者,其振幅受余流影響,因此余流影響了半日潮周期內懸沙波峰之間的不對稱性,由于余流方向同漲潮流,漲潮時段懸沙波峰高于落潮時段內的波峰。