王昆，李楠，宋廣軍，關(guān)春江，胡超魁，杜靜，邵澤偉，吳金浩*

(1.遼寧省海洋水產(chǎn)科學(xué)研究院，遼寧 大連 116023； 2.國家海洋環(huán)境監(jiān)測中心，遼寧 大連 116023)

遼東灣位于渤海的最北部，總水域面積約為10 000 km2，水深較淺，地形多變，潮汐類型復(fù)雜。多年來，人們針對遼東灣海域進(jìn)行了大量的水文動(dòng)力、泥沙輸運(yùn)和水環(huán)境保護(hù)等方面的數(shù)值模擬研究[1-4]，特別是在2014年以來國家開展的渤海綜合治理過程中，相關(guān)部門在遼東灣水環(huán)境保護(hù)研究方面取得了較多的成果。

水體中的浮游植物含量是海洋初級生產(chǎn)力的主要貢獻(xiàn)指標(biāo)，是水體中貝類、牡蠣等生物餌料(單細(xì)胞藻類)的主要貢獻(xiàn)體，在海洋生態(tài)系統(tǒng)和環(huán)境保護(hù)中起著至關(guān)重要的作用。而葉綠素a的濃度往往用來定量表征海水中浮游植物的豐度，在海洋初級生產(chǎn)力的計(jì)算與評價(jià)中，葉綠素a的濃度是主要的決定性因素，故很多研究中常將葉綠素a濃度作為主要水體研究指標(biāo)，并探討其與水動(dòng)力的響應(yīng)關(guān)系[5-7]。相關(guān)研究中，徐文龍等[8]基于光學(xué)觀測數(shù)據(jù)，采用高精度反演和數(shù)理統(tǒng)計(jì)手段對南海東北部海區(qū)夏季葉綠素a濃度的垂向變化特征及其對水動(dòng)力過程的響應(yīng)進(jìn)行了研究；He等[9]基于南海布放生物地球化學(xué)浮標(biāo)數(shù)據(jù)，指出南海中央海盆存在一定深度的葉綠素極大值分布區(qū)。在遼東灣海區(qū)的相關(guān)研究中，Pei等[10]于2013年夏季，在遼東灣進(jìn)行了一項(xiàng)實(shí)地調(diào)查研究，以確定該海灣浮游植物的動(dòng)態(tài)及其初級生產(chǎn)力受環(huán)境制約的程度，結(jié)果顯示，遼東灣的區(qū)域生產(chǎn)受到溫度和光照限制的綜合控制；2016—2017年，田思瑤等[11]在遼東灣中部近岸海域，采用單因子標(biāo)準(zhǔn)指數(shù)法和綜合指數(shù)評價(jià)法對其COD、石油類、葉綠素分布特征及富營養(yǎng)化狀態(tài)進(jìn)行了評價(jià)，結(jié)果顯示，調(diào)查海域內(nèi)COD和石油類指標(biāo)分布呈河口及沿岸濃度較高的特點(diǎn)，存在明顯的季節(jié)變化特征，在監(jiān)測海域中、北部以輕度、中度富營養(yǎng)化狀態(tài)為主，富營養(yǎng)化程度較高的站位為河流入海區(qū)部分站位；2018年，裴紹峰等[12]采用觀測與模擬分析相結(jié)合的手段，對遼東灣海區(qū)的葉綠素a和初級生產(chǎn)力進(jìn)行研究，結(jié)果顯示，葉綠素a低值區(qū)出現(xiàn)在雙臺子河的河口，其原因是河水中過量泥沙懸浮物降低了該海區(qū)水體的透明度，從而導(dǎo)致浮游植物生長受到光的限制。綜上可以看出，前期大部分研究都是基于監(jiān)測數(shù)據(jù)，采用數(shù)理統(tǒng)計(jì)方法來評價(jià)水質(zhì)指標(biāo)在水體中的時(shí)空分布特征，即使考慮多指標(biāo)的水質(zhì)模型，也假定為各指標(biāo)間不發(fā)生任何生化反應(yīng)的保守型問題，但實(shí)際水體存在溫度、鹽度和光照等因素的影響，各個(gè)指標(biāo)間會發(fā)生復(fù)雜的生物化學(xué)反應(yīng)過程，若按單一保守型的變化響應(yīng)來理想化描述，其結(jié)果的呈現(xiàn)往往偏離實(shí)際物理規(guī)律的本質(zhì)。為了解決這一問題，本研究中在遼東灣水文動(dòng)力數(shù)值模型構(gòu)建的基礎(chǔ)上，引入水體中葉綠素a、溶解氧、氨氮、亞硝酸氮、硝酸氮和活性磷酸鹽各指標(biāo)間的生化反應(yīng)過程，建立遼東灣海域以葉綠素a濃度為主要未知量的生態(tài)系統(tǒng)動(dòng)力學(xué)模型。該模型結(jié)合遼東灣海區(qū)的水動(dòng)力條件，基于實(shí)際觀測信息，綜合考慮可能影響到整個(gè)生態(tài)系統(tǒng)的物理、化學(xué)和生物變化過程，模擬預(yù)測多工況下葉綠素a濃度的分布結(jié)果，并用2020年遼東灣的觀測數(shù)據(jù)對該模型進(jìn)行驗(yàn)證，以期為遼東灣海區(qū)的生態(tài)系統(tǒng)動(dòng)力學(xué)過程分析提供合理方案，為水體富營養(yǎng)化和赤潮預(yù)測預(yù)報(bào)研究提供更為科學(xué)有效的技術(shù)手段。

1 水動(dòng)力控制方程

利用二維水動(dòng)力模型對遼東灣海域的水動(dòng)力場進(jìn)行模擬，以期數(shù)值再現(xiàn)和預(yù)測該海域水體運(yùn)動(dòng)過程，為浮游植物在水體中的分布模擬提供背景場[13]?；诓豢蓧嚎s流體Reynolds平均Navier-Stokes方程，將其水平動(dòng)量方程和連續(xù)性方程在總水深[0,h]范圍內(nèi)進(jìn)行積分后，可得如下的二維深度平均淺水方程組。

連續(xù)性方程：

(1)

動(dòng)量方程：

(2)

(3)

2 邊界條件

給定兩種邊界條件，即閉邊界和開邊界條件。

1)閉邊界條件。閉邊界即水陸交界邊界，一般由海岸線及海島確定，在閉邊界處法向流速為0。

2)開邊界條件。本模型中開邊界處采用如下的分潮調(diào)和形式計(jì)算：

(4)

其中：m為考慮的分潮個(gè)數(shù)；ωi為第i個(gè)分潮的角速度；fi、ui分別為第i個(gè)分潮的交點(diǎn)因子和遲角修正；Hi和gi為調(diào)和常數(shù)，即分別為各分潮的振幅和遲角；V0i為第i分潮的天文初相位。

3)干濕動(dòng)邊界。對于潮灘，水陸交界的位置隨著潮位漲落而變化，本模型中考慮了動(dòng)邊界內(nèi)網(wǎng)格節(jié)點(diǎn)的干濕變化[14]。

4)初始條件。計(jì)算開始時(shí)采用“冷態(tài)”啟動(dòng)，即各物理量初值均為0。

5)離散方法。模型對計(jì)算區(qū)域的空間離散采用有限體積法，對不同的計(jì)算區(qū)域采用非結(jié)構(gòu)三角形網(wǎng)格剖分形式，大大增強(qiáng)了系統(tǒng)對岸線和地形變化的適應(yīng)性，提高了計(jì)算精度。

3 生態(tài)系統(tǒng)動(dòng)力學(xué)模型

1)描述葉綠素a濃度的生態(tài)系統(tǒng)動(dòng)力學(xué)模型?？刂品匠虨?/p>

(5)

其中：u、v為水動(dòng)力方程中的垂向平均流速；Ex、Ey分別為x、y軸兩方向的紊動(dòng)擴(kuò)散系數(shù)；S為源匯項(xiàng)；F為生化反應(yīng)項(xiàng)。

(6)

其中： Chlacs為浮游植物生長過程對葉綠素a的產(chǎn)生過程項(xiàng)；Chlahx為呼吸過程對葉綠素a的改變過程項(xiàng)；Chlasw為死亡過程對葉綠素a的改變過程項(xiàng)；Chlacj為沉降過程對葉綠素a的改變過程項(xiàng)。

2)影響葉綠素a分布的4個(gè)主要生態(tài)過程。

浮游植物生長過程引起的葉綠素a濃度變化的表達(dá)式：

Chlacs=Gf×a1×a2×GN。

(7)

其中：Gf為光合作用過程函數(shù)；a1、a2為浮游植物生長速率[1](表1)；GN為營養(yǎng)限制因子[15]，且

GN=Amax×e(-1.6dz/θgs)×Φd/dz。

(8)

式中：Amax代表海區(qū)的正午最大產(chǎn)氧量[16]；dz為水層厚度；θgs為光衰減系數(shù)[17]；Φd為太陽輻射強(qiáng)度的日變化函數(shù)[18]。

呼吸過程引起的葉綠素a濃度變化的表達(dá)式：

Chlahx=Of×a1×a2×GN。

(9)

其中，Of為浮游植物呼吸過程的需氧量函數(shù)，且

(10)

式中：CDO為浮游植物呼吸過程中所需的溶解氧半飽和濃度；sfz為浮游植物呼吸速率[19]；a5為溫度依賴系數(shù)；T為水溫(℃)。

死亡過程引起的葉綠素a濃度變化的表達(dá)式：

在相同作業(yè)溫度下，更高的催化劑活性能夠提升苯乙烯單體的產(chǎn)量。高選擇性則能生成更高比例的苯乙烯單體，降低副產(chǎn)物的生成，如苯和甲苯等。當(dāng)前可用的苯乙烯單體催化劑只能在高活性和高選擇性之間二選一。這是因?yàn)樵诔蚐HR作業(yè)條件下，催化劑活性的上升通常會導(dǎo)致選擇性下降。

Chlasw=a3×Chla。

(11)

其中，a3為浮游植物的死亡率[1]。

沉降過程引起的葉綠素a濃度變化的表達(dá)式：

Chlacj=a4/dz×Chla。

(12)

其中，a4為沉降速率[20]。

3)主要生態(tài)過程求解需要的6個(gè)子系統(tǒng)[21]。

生化需氧量子系統(tǒng)1：

(13)

其中：JBOD為BOD的降解函數(shù)，與BOD的降解速率、溫度依賴系數(shù)和溶解氧半飽和濃度有關(guān)。

氨氮子系統(tǒng)2：

(14)

其中：SN為氨氮釋放BOD的衰減函數(shù)，與BOD衰減釋放的氨氮比例與降解有關(guān)；FfN、FxN分別為浮游植物和細(xì)菌對氨氮的吸收函數(shù)，與水深、吸收的氨氮量及溶解氧半飽和濃度有關(guān)；Sx為硝化過程中氨氮轉(zhuǎn)化為亞硝酸鹽的硝化速率函數(shù)，同樣與相應(yīng)過程的衰減速率、溫度依賴系數(shù)和溶解氧半飽和濃度有關(guān)。

溶解氧子系統(tǒng)3：

(15)

其中：Fat為大氣復(fù)氧函數(shù)，與大氣復(fù)氧速率和溶解氧濃度有關(guān)；Fx為需氧量的硝化速率函數(shù)，與單位質(zhì)量的氨氮到亞硝酸鹽、單位質(zhì)量的亞硝酸鹽到硝酸鹽硝化作用的需氧量，硝化作用中20 ℃時(shí)的一級衰減速率，以及亞硝酸鹽到硝酸鹽衰減過程的溫度依賴系數(shù)等有關(guān)；Fsed為沉積物的需氧量函數(shù)，由底層碎屑的半飽和濃度、單位面積沉積物的需氧量及沉積物衰變的溫度依賴系數(shù)控制。

亞硝酸鹽子系統(tǒng)4：

(16)

硝酸鹽子系統(tǒng)5：

(17)

其中：K2、K3分別為溫度20 ℃時(shí)水體的一級硝化和反硝化速率；Cxo為硝化作用中所需的溶解氧半飽和濃度；θx、θfx分別為硝化和反硝化的溫度依賴系數(shù)。

活性磷酸鹽子系統(tǒng)6：

(18)

其中：PrB為每克溶解態(tài)BOD中的磷含量；Pfz為磷系統(tǒng)中每克浮游植物磷的吸收量。

將以上微分方程及相關(guān)數(shù)學(xué)表達(dá)式(6)～(18)聯(lián)立起來組成一個(gè)七元一階的微分方程組，利用Euler法或四階Runge-Kutta法[22-25]求解。然后，將生化反應(yīng)項(xiàng)代入方程(5)，以方程(1)～(3)求得的水動(dòng)力場為背景場，便可得到葉綠素a濃度的時(shí)空分布。

4 遼東灣海域浮游植物分布特征模擬

浮游植物的豐度分布用葉綠素a的濃度來表征，采用上述生態(tài)系統(tǒng)動(dòng)力學(xué)模型對遼東灣海域的水質(zhì)和葉綠素a濃度的分布狀況進(jìn)行模擬分析。

4.1 計(jì)算域及參數(shù)設(shè)置

計(jì)算域范圍為整個(gè)渤海范圍，整個(gè)模擬區(qū)域內(nèi)由25 888個(gè)節(jié)點(diǎn)和47 197個(gè)三角單元組成，最小空間步長約為80 m，時(shí)間步長為30 min。地形取自自然資源部網(wǎng)站，插值得到研究海域的水深，岸界取自海域水深地形的平面圖。計(jì)算中所需常量的取值見表1。

表1 計(jì)算中所需的常量取值Tab.1 Constant value required in calculation

4.2 模型驗(yàn)證與模擬結(jié)果

為了獲取準(zhǔn)確可靠的數(shù)值模型，需要利用觀測值對模型進(jìn)行驗(yàn)證。分別于2020年5月和8月，在遼東灣海域布設(shè)4個(gè)水文和14個(gè)生態(tài)監(jiān)測站位，對遼東灣海域的水動(dòng)力和生態(tài)狀況進(jìn)行了監(jiān)測分析。水文觀測站位置坐標(biāo)見表2，站位位置見圖1。

表2 潮位和海流觀測站位位置Tab.2 Location of tide level and current observation stations

圖2給出了監(jiān)測站位H1、H2實(shí)測潮位值和計(jì)算值的比較，可以看出，監(jiān)測站位處潮汐屬于不規(guī)則半日潮，即一天24 h內(nèi)出現(xiàn)兩次漲潮和兩次落潮過程，漲落潮歷時(shí)大致相同，計(jì)算和實(shí)測潮位過程的高、低潮位及過程線均吻合良好。說明本研究中的水動(dòng)力模型模擬的海域潮波運(yùn)動(dòng)與天然潮波運(yùn)動(dòng)基本相似，采用的邊界條件合理，能夠真實(shí)地重演和反映海域內(nèi)潮波的傳遞和變形。監(jiān)測站位v1、v2的流速和流向模擬值與實(shí)測值對比見圖3、圖4，兩個(gè)監(jiān)測點(diǎn)的流速和流向模擬過程線與實(shí)測值吻合良好，呈現(xiàn)出海區(qū)典型的NE-SW往復(fù)流分布特征，v1站位略顯旋轉(zhuǎn)流特征。綜上可見，本研究中水動(dòng)力模型的模擬結(jié)果具有較高的精確性，能夠滿足后期生態(tài)系統(tǒng)指標(biāo)濃度值模擬計(jì)算的需要。

本圖基于自然資源部標(biāo)準(zhǔn)地圖服務(wù)網(wǎng)站GS(2019)3266號標(biāo)準(zhǔn)地圖為底圖制作，底圖邊界無修改。The figure is based on the standard map GS(2019)3266 in the standard map service website of Ministry of Natural Resources of the People’s Republic of China, with no modifications of the boundaries in the standard map.圖1 監(jiān)測站位位置圖Fig.1 Location map of surveyed station

圖2 H1和H2站位的潮位驗(yàn)證Fig.2 Tidal level verification of H1 and H2 stations

為進(jìn)一步增強(qiáng)模型的可靠性，本研究中對生態(tài)系統(tǒng)動(dòng)力學(xué)模型的生態(tài)指標(biāo)濃度模擬結(jié)果進(jìn)行了驗(yàn)證，并基于2020年5月和8月的生態(tài)指標(biāo)監(jiān)測值來進(jìn)行分析。表3給出了5月和8月2個(gè)航次14個(gè)站位葉綠素a濃度監(jiān)測值和模擬值的對比，可以看出，模型的模擬結(jié)果與實(shí)測值吻合良好。同時(shí)，為了完整表達(dá)本研究中生態(tài)系統(tǒng)動(dòng)力學(xué)模型的模擬結(jié)果，以及與保守水質(zhì)模型[25]相比，體現(xiàn)該模型的精確度與優(yōu)勢度，圖5、圖6分別給出了2020年5月3個(gè)生態(tài)指標(biāo)的濃度變化過程及一個(gè)潮期內(nèi)4個(gè)典型潮時(shí)的葉綠素a濃度場的平面分布狀況。

圖3 v1站位的流速驗(yàn)證Fig.3 Tidal current verification of v1 station

圖4 v2站位的流速驗(yàn)證Fig.4 Tidal current verification of v2 station

表3 葉綠素a濃度實(shí)測值與模擬值對比Tab.3 Comparison of chlorophyll a concentration results between measured and simulated values

圖5 遼東灣灣頂海區(qū)生態(tài)指標(biāo)濃度模擬值隨時(shí)間的變化Fig.5 Change in simulated values of ecological index concentrations with time at top sea area in Liaodong Bay

圖6 遼東灣4個(gè)典型潮時(shí)的葉綠素a濃度模擬場分布Fig.6 Distribution of simulated field of chlorophyll a concentration at four typical tides in Liaodong Bay

模擬結(jié)果表明：葉綠素a的分布具有明顯的時(shí)空分布特征，遼東灣靠近遼河入?？凇⒑J島和長興島附近海區(qū)的葉綠素a濃度明顯高于其他海區(qū)；8月的平均濃度明顯高于5月，二者的比值在1.6左右(表3)。葉綠素a濃度的模擬值在5月中旬達(dá)峰值，為13.4 μg/L(圖5)，之后由于藻類呼吸和沉降過程，又引起濃度的顯著下降；而在8月，隨著氣溫和溶解氧的升高，藻類的量會快速累加，葉綠素a的濃度則達(dá)到一年內(nèi)的另一個(gè)峰值，結(jié)果與歷年觀測的浮游植物豐度呈現(xiàn)雙峰結(jié)構(gòu)的分布趨勢一致，與實(shí)際情況基本吻合。

5 討論

5.1 生態(tài)指標(biāo)濃度的時(shí)間分布特征

近年來，以葉綠素a為主要表征指標(biāo)的水體中浮游植物分布密度的研究廣泛開展。本研究在借鑒前人基于監(jiān)測數(shù)據(jù)數(shù)理統(tǒng)計(jì)分析的基礎(chǔ)上，綜合考慮實(shí)際水體的流動(dòng)性特征，真實(shí)展現(xiàn)可能影響到整個(gè)生態(tài)系統(tǒng)的物理、化學(xué)和生物變化過程，在水動(dòng)力數(shù)值模型基礎(chǔ)上引入各生態(tài)指標(biāo)之間的耦合作用，建立了遼東灣海域的生態(tài)系統(tǒng)動(dòng)力學(xué)模型。本模型的數(shù)值模擬結(jié)果表明，其指標(biāo)濃度分布具有明顯的時(shí)間分布特征。提取遼東灣灣頂海區(qū)葉綠素a濃度的時(shí)間序列曲線圖發(fā)現(xiàn)，在5月中旬濃度達(dá)到峰值13.4 μg/L，到5月15日基本是平穩(wěn)期，而后由于藻類呼吸和沉降過程，5月16日之后濃度出現(xiàn)了顯著的下降趨勢，5月26日到6月初濃度趨于穩(wěn)定。8月葉綠素a又出現(xiàn)另外一個(gè)峰值，平均濃度明顯高于5月，二者的比值在1.6左右，表明葉綠素a濃度即浮游植物分布存在明顯的季節(jié)性分布特征。而無機(jī)氮濃度呈現(xiàn)出5月初到6月初持續(xù)遞減的趨勢；活性磷酸鹽濃度在5月6日達(dá)到最高值，然后隨時(shí)間變化至5月18日(12 d左右)時(shí)濃度大幅度下降，到5月23日(5 d左右)，活性磷酸鹽的濃度開始緩慢降低，最后基本以極低的濃度值穩(wěn)定在6月初。與2013年建立的水質(zhì)模型[1，25]模擬結(jié)果相比較，本模型在整個(gè)模擬過程中，生態(tài)指標(biāo)濃度量值變化明顯，規(guī)律趨勢分辨率更高，模型計(jì)算結(jié)果在詳細(xì)反映實(shí)際的內(nèi)部生態(tài)機(jī)理變化特征方面更有優(yōu)勢，準(zhǔn)確性和可靠性更高。

5.2 生態(tài)指標(biāo)濃度的空間分布特征

諸多的研究成果中，采用斷面觀測和大面觀測的方式來反映水體中指標(biāo)濃度的變化狀況，其全面性尚待提高。數(shù)值模型的典型特征就是在重演歷史時(shí)期有限數(shù)據(jù)的基礎(chǔ)上，通過模型驗(yàn)證和調(diào)參過程，進(jìn)而預(yù)測未來時(shí)間段整個(gè)空間平面或立體的數(shù)據(jù)信息，使結(jié)果分析更加具有說服力。為了較全面地展現(xiàn)生態(tài)系統(tǒng)動(dòng)力學(xué)模型的研究結(jié)果，本研究中，在給出各生態(tài)指標(biāo)時(shí)間序列曲線結(jié)果的同時(shí)，也給出各生態(tài)指標(biāo)模擬結(jié)果在遼東灣海域整個(gè)空間的分布狀況。葉綠素a濃度的模擬平面圖表明，遼東灣灣頂海區(qū)的葉綠素a濃度明顯偏高，而且經(jīng)統(tǒng)計(jì)分析發(fā)現(xiàn)，高、低潮時(shí)的濃度高值區(qū)范圍均大于漲急、落急時(shí)刻的濃度高值區(qū)，遼東灣遼河口入海區(qū)、葫蘆島和長興島附近海區(qū)的葉綠素a濃度明顯高于其他海區(qū)，這也與歷史上遼東灣夏季赤潮易發(fā)區(qū)分布位置相一致。分析原因，可能是由于夏季遼河口處于豐水期，排入河口入海區(qū)的營養(yǎng)鹽比較豐富，使灣頂海區(qū)的藻類大量繁殖，而且葉綠素a濃度高值區(qū)的位置分布也與遼東灣灣頂海區(qū)及葫蘆島一帶海區(qū)為產(chǎn)卵場與索餌場的漁場特征相吻合[26-27]，體現(xiàn)出遼東灣葉綠素a濃度明顯的空間分布特征。

6 結(jié)論

1)遼東灣海區(qū)葉綠素a的濃度分布具有明顯的時(shí)空分布與季節(jié)性分布特征，本研究中建立的浮游植物分布生態(tài)系統(tǒng)動(dòng)力學(xué)模型，對生態(tài)指標(biāo)濃度的變化規(guī)律性展現(xiàn)出較以往保守型模型具有更高的分辨率，模擬結(jié)果與歷年觀測的浮游植物豐度分布趨勢的雙峰結(jié)構(gòu)特征相吻合，說明模型具有較高的精確性。

2)本模型模擬的生態(tài)指標(biāo)濃度分布規(guī)律與遼東灣歷年的赤潮易發(fā)區(qū)和漁場分布特征一致，說明模型中考慮的整個(gè)生態(tài)系統(tǒng)的物理、化學(xué)和生物變化過程，能夠反映出模擬時(shí)間段海水中的客觀自然變化規(guī)律，證明文中建立的模型在反映內(nèi)部生態(tài)指標(biāo)間相互影響機(jī)理方面表達(dá)更準(zhǔn)確，具有靈敏性，說明模型具有較高的可靠性。

3)對遼東灣海域浮游植物分布模型的模擬結(jié)果，可為根據(jù)葉綠素a高值區(qū)的時(shí)空分布位置科學(xué)推測和跟蹤漁業(yè)資源密集區(qū)，以及提高捕撈效率或選擇增殖放流海區(qū)提供重要參考；葉綠素a模擬結(jié)果對與浮游植物密度相關(guān)性極強(qiáng)的浮游動(dòng)物，如毛蝦集群效應(yīng)等研究，也具有預(yù)測和預(yù)警作用，進(jìn)而也可為規(guī)避電廠取水口堵塞風(fēng)險(xiǎn)、保障產(chǎn)業(yè)的安全穩(wěn)定運(yùn)行提供技術(shù)支撐，說明該模型具有較好的實(shí)用性。