渤海氮磷營養鹽和葉綠素濃度時空分布數值模擬與富營養化評估*

趙 宇 李克強① 孫 珊 陳 衎 譚光深 張 娟 王修林

(1. 中國海洋大學海洋化學理論與工程技術教育部重點實驗室 中國海洋大學化學化工學院 山東青島 266100; 2. 山東省海洋資源與環境研究院 山東煙臺 264006)

當前, 富營養化是海灣主要的生態環境問題。全球經濟中心大多坐落在環灣地區, 隨著陸源污染物入海排放等人類活動壓力逐漸增大, 海灣富營養化評估和污染防治成為全球陸海環境綜合治理的重點。自1960 年代末特別是1972 年聯合國第一次人類環境會議以來, 美國、日本、歐盟諸國相繼開始近海富營養化綜合治理, 其中, 厘清水體富營養化狀態與變化趨勢是綜合治理的關鍵, 獲取營養鹽和葉綠素等富營養化評估要素長期的演變進程高質量數據是前提和難點。在海洋調查監測數據部分年代缺失條件下,數值模型模擬計算是重構營養鹽和葉綠素等富營養化評估要素數據集的有效手段, 當前逐步發展建立了EFDC、WASP、FVCOM、ECOHAM、ERSEM、COHERENS 等多種數值模型計算方法(任湘湘, 2011;樊嘉蓉, 2022)。按照DPSIR 理論等, 海灣營養鹽、有機質、葉綠素等富營養化評估要素, 主要是環灣地區各種人類經濟社會活動壓力共同作用的結果。這不僅表現在陸源污染物入海排放等改變氮磷營養鹽等化學變量的濃度、結構和組成, 也表現在圍填海、海岸工程等通過改變流速、流向等水動力場物理變量, 間接影響營養鹽、有機質、葉綠素等生地化變量的時空分布的過程。因此, 海灣富營養化評估要素的數值模擬需要物理、化學和生物多種變量的模擬計算, 其中生物生態過程是長時間的演變過程, 需要進行多變量長時序模擬計算。

海灣富營養化評估要素(營養鹽、有機質、葉綠素等)的模擬計算, 相關研究始于20 世紀40 年代(房月英, 2008)。早期研究聚焦富營養鹽、葉綠素等單變量季節或時空分布模擬, 至20 世紀80 年代(Skogenet al, 2000; Molletal, 2003), 海洋三維水動力/生物地球化學模型模擬研究進一步完善和發展, 主要包括海灣富營養化要素控制過程及作用機制研究、模型邏輯架構構建及修正、參數率定、模擬計算等。海灣富營養化要素時空變化是水動力輸運、生物地球化學遷移轉化、浮游植物生長等過程綜合作用的結果。對于海灣水動力輸運過程研究, 通常采用現場觀測結合數值模擬的方法研究流場中潮流和余流變化規律(李高陽, 2015)。對于營養鹽生物地球化學遷移轉化和浮游植物生長過程, 通常采用陸海同步調查、圍隔生態系/培養瓶現場培養實驗和數值模型模擬方法研究(Chenetal, 2022)。對于模型邏輯架構, 海灣水動力模型主要包括ROMS、FVCOM、POM、HAMSOM 等, 目前廣泛采用ROMS 和FVCOM; 生物地球化學和生態動力學模型主要包括WASP、NPZD、ERSEM 等(劉學海, 2009)。對于模型參數率定, 通常采用伴隨矩陣法(呂咸青等, 2002)、模擬退火法(張艷杰等, 2006)、蒙特卡羅法(朱陸陸, 2014)等。對于模擬計算技術, 通常采用粒子追蹤法(翁怡嬋, 2009)、差分法(Kristeket al, 2010)和分步模擬計算方法(齊寧等, 2020)等。

最近, Gao 等(2020)發展完善了海洋圍隔生態系/實驗瓶現場培養實驗方法, 這不僅揭示了萊州灣浮游植物優勢種與不同形態氮遷移轉化過程之間的關系, 而且從動力學上量化了萊州灣幾種典型浮游植物優勢種演替的營養鹽控制過程和作用機制。同時,Li 等(2017)根據渤海三維水動力/生物地球化學耦合模型, 應用陸海水環境同步監測調查數據, 將渤海溶解無機氮(dissolved inorganic nitrogen, DIN)等濃度時空分布模擬計算結果的準確性平均提高約25%。這樣,圍繞渤海富營養化評估要素的數值模擬計算, 針對渤海三維水動力/生物地球化學耦合模型在邏輯架構、過程耦合、參數地域化等存在的問題, 本文在改進渤海三維水動力/生物地球化學過程耦合模型架構基礎上, 結合陸海同步調查監測數據校正、參數率定、多變量模擬計算等方法, 建立完善了渤海富營養化評估要素長期演變數據重構技術方法, 據此利用修訂的營養狀態質量指數(Nutritional Quality Index,NQI), 科學評估了渤海富營養化狀態和變化趨勢,為進一步提高中國推行的渤海重點海域“深入打好污染防治攻堅戰”行動計劃的實施提供堅實的科技支撐。

1 研究方法

海灣多變量模擬計算準確性不僅取決于模型的邏輯架構和監測數據的質量, 也取決于模型參數的靈敏度。這樣, 海灣富營養化評估要素模擬計算主要包括三維水動力/生物地球化學耦合數值模型構建、參數率定、多變量模擬計算與準確性評價等步驟。

1.1 數值模型構建

王強(2004)和 Dai 等(2015)在 Hamburg Shelf Ocean Model (HAMSOM)(Backhausetal, 1983)基礎上建立了渤海三維水動力/生物地球化學耦合模型,本文在此基礎上, 對相關生物地球化學模塊進行修正。模型中入?？谙嚓P數據參見戴愛泉(2015)研究結果。

根據溶解有機氮(dissolved organic nitrogen, DON)的現場培養實驗結果(戴愛泉, 2015; Lietal, 2017),不同來源的DON 結構組成和降解速率存在顯著差異。本文按照污水處理廠和河流不同來源DON 的降解速率差異性, 將生物地球化學過程模型中DON 狀態變量分為污水處理廠輸入生物易降解DON(BDON1)、河流輸入生物易降解 DON(BDON2)和生物難降解DON(RDON), 方程修正如下:

式中, DON1 為陸源污水處理廠輸入的生物易降解DON 濃度; DON2 為河流輸入的生物易降解DON 濃度; DON3 為各排污口排入的難降解 DON 濃度;DONm為海源DON 濃度;KrDON為難降解DON 組分占比; KB1、KB2、KB3分別為DON1、DON2、DON3的降解速率常數;KtDON為溶解有機氮微生物降解溫度效應系數; temp 為環境溫度。

結合生物地球化學模塊其他狀態變量遷移轉化過程, 建立的模型主要框架及過程如圖1 所示。模型主要包括浮游植物、浮游動物、溶解無機態營養鹽、溶解有機態營養鹽、生物碎屑5 個子模塊。共包含DIN、陸源生物易降解DON (BDON1, BDON2)、海源DON (BDONm)、生物難降解DON(RDON)、無機磷(DIP)、有機磷(DOP)、浮游植物(PHYTO)、浮游動物(ZOO)、碎屑(DPT)10 個狀態變量。物質流統一用氮元素(單位: mmol/m3)表示, 浮游植物生物量可表示為活性顆粒態氮(PPON), chla濃度通過浮游植物葉綠素氮比值關系進行換算(表1)。主要的生物地球化學過程有:浮游植物生長、死亡和營養鹽吸收分泌等過程; 浮游動物生長、死亡和營養鹽釋放等過程; 無機態和有機態營養鹽吸收和轉化過程; 碎屑礦化和沉降過程等。相關生物地球化學過程動力學方程詳見附表Ⅰ。

表1 模型動力學過程所用參數意義及取值Tab.1 Significance and value of parameters used in model dynamics process

圖1 渤海三維水動力/生物地球化學耦合模型架構示意圖Fig.1 Schematic architecture of three-dimensional hydrodynamic/biogeochemical multi-process coupled model for the Bohai Sea

1.2 模型參數率定

1.2.1 參數率定 參數設置在主要參考文獻值基礎上, 采用渤海萊州灣海域現場培養實驗結果(李克強等, 2007; Lietal, 2009, 2017), 利用DIN、DIP、DON、chla濃度季節變化平面分布模擬結果與實測結果對比進行校正, 利用年際變化模擬結果與實測結果對比進行驗證, 模型參數結果如表1 所示。

1.2.2 數據來源

模型校驗數據主要取自中國近岸海域環境質量公報(報告書)、中國海洋環境質量(狀況)公報、北海區海洋環境公報、中國海洋生態環境狀況公報和文獻結果, 詳見表2。其中, 2019 年氮磷營養鹽和chla濃度來源于國家自然科學基金委渤黃海共享航次(航次編號: NORC2018-01; NORC2019-01)調查結果。

表2 模型模擬計算結果驗證數據來源Tab.2 Data source for validation on the modelling results

1.3 檢驗評價

本文選用余弦相似性系數(similarity index, SI;方程 4)和相對標準偏差(relative standard deviation,RSD; 方程5)分別評價模擬值和實測值平面分布趨勢相似性以及模擬值偏離實測值的程度(戴愛泉,2015), 從而綜合定量評價模型模擬結果的準確性:

其中,Oi和Pi分別表示第i個水質點上的實測值和模擬值;n表示監測站位數量。

針對渤海DIN、DIP、DON 與chla濃度分布場,同時使用Kappa 系數評價模型模擬計算準確性。采用Cohen 提出的Kappa 系數評價分類標準(Cicchettietal,1990; Feinsteinetal, 1990)(表3)。

表3 Kappa 系數評價分類標準Tab.3 Classification criteria for Kappa coefficient evaluation

1.4 富營養化評價方法

渤海海水中溶解有機氮(DON)的濃度近年來持續升高, 某些海域已成為主要的氮形態(李志林等,2015)。在富營養化評價方法中, 用DON 代替化學需氧量(COD)能更好地評估富營養化狀態。本文應用營養狀態質量指數(NQI)法(陳于望, 1987), 采用DON代替COD 予以修訂, 計算公式如下:

式中,C為下標各要素DON、DIN、DIP 與chla的濃度值, 上標*表示各要素DON、DIN、DIP 與chla的營養狀態評價標準。其中, chla的營養狀態評價標準采用赤潮發生葉綠素濃度閾值, 其他評價要素參照海水水質標準(GB 3097-1997), DON 采用DIN 標準值,各要素營養狀態評價標準分別取值為這樣, 通過各要素濃度與標準的比值可消掉量綱。NQI＞3 為富營養狀態, 3≥NQI≥2 為中等營養狀態, NQI＜2 為貧營養狀態, NQI 指數越高富營養化程度則越高。

2 結果與討論

2.1 渤海氮、磷和chl a 濃度時空變化模擬計算結果

2.1.1 DIN、DIP 和DON 平面分布 應用構建完善的三維水動力/生物地球化學耦合數值模型, 模擬計算了2019 年DIN、DIP、DON 春夏秋冬四個季節的濃度平面分布(圖2)。經模型參數校正, 再與調查站位上的實測結果對比, DIN、DIP、DON 濃度模擬結果在數值大小及變化趨勢上均具有較高的一致性,三者SI 平均值為0.77±0.08, RSD 平均值為26%±5%,Kappa 系數平均值為0.57±0.08。

圖2 2019 年全年渤海DIN、DIP 與DON 濃度平面分布模擬計算結果與陸海同步監測調查結果比較圖Fig.2 Comparison between the simulated and observed (in land sea synchronous monitoring survey) results of DIN, DIP, and DON in concentration plane distribution in the Bohai Sea in 2019

結果表明, 模擬計算值再現了渤海的DIN、DIP與DON 濃度平面分布春夏秋冬的季節變化特征。在春、夏兩季, DIN 濃度因浮游植物生長被消耗利用, 渤海灣與遼東灣南部濃度都比較低, 萊州灣西部受黃河徑流輸入的影響, 成為夏季渤海DIN 濃度最高的海區(圖2b, 2d); 秋季渤海DIN 濃度高值區不斷擴大,渤海中部濃度相對較低, 但也高于夏季水平; 冬季渤海DIN 濃度達到最高峰, 濃度高值區較冬季繼續擴大, 渤海中部平均濃度達到秋季的兩倍水平, 僅有渤海海峽東部濃度較低(圖2h)。由于浮游植物在春、夏季生長的消耗, 海水中的DIP 濃度較低(圖2j, 2l); 秋季渤海灣南岸、萊州灣以及遼東灣北岸DIP 濃度較高,渤海中部濃度較低, 夏秋季是我國北方降雨期, 三大灣沿岸DIP 濃度上升與陸地徑流輸入有關, 從而形成了由沿岸高值區向渤海中部逐漸遞減的趨勢(圖2l,2n)(王燕等, 2021); 冬季渤海中DIP 含量上升, 主要因為冬季時淺海風力增大, 導致躍層失去了穩定性,從而促使海水底層顆粒的有機磷分解所產生的DIP隨著垂直混合被帶到了表層(圖2p)(黃爽, 2012)。在春、秋和冬季, DON 高值區主要集中于萊州灣西部、渤海灣、秦皇島近岸、遼東灣北部及金州灣沿岸地區(圖2r, 2v, 2x); 而夏季DON 受高值區出現在渤海中部(圖2t), 模型沒有很好的模擬再現這一特征, 這與DON 等有機質受浮游植物和微生物調控有關, 由于浮游植物夏季較高的生物量, DON 的源匯過程較為復雜, 模型模擬很難實現精準刻畫, 但在調查站位上與調查監測結果比較, 31%的RSD 值還是符合海上調查監測平均誤差范圍(約30%)的。

在季節變化上, 2019 年渤海夏季DIN、DIP 濃度顯著低于其他季節, 而秋、冬季DON 濃度高于其他季節(圖2), 這與王婷(2009)在2002～2004 年的研究結果在變化規律上基本一致, 構建的模型很好的體現了黃河調水調沙引起的營養鹽輸入通量變化對DIN 和DIP 濃度平面分布產生的影響(王婷, 2009; Dingetal, 2020)。

2.1.2 chla濃度平面分布 應用模型模擬計算了2019 年chla春夏秋冬四個季節的濃度平面分布(圖3)。在調查站位上與調查監測結果對比, chla濃度模擬結果在數值大小及變化趨勢上均具有較高的一致性, SI平均值為0.77±0.04, RSD 平均值為23%±2%, Kappa系數平均值為0.63±0.13。

圖3 2019 年渤海chl a 濃度平面分布模擬計算結果與陸海同步監測調查結果比較圖Fig.3 Comparison between the simulated and observed (in land sea synchronous monitoring survey) results of chl a concentration plane distribution in the Bohai Sea in 2019

渤海的chla濃度平面分布模擬計算結果體現了浮游植物春夏秋冬的季節變化特征。春季chla濃度的高值區出現在黃河口附近海域、萊州灣南部、渤海灣與遼東灣北部、秦皇島外海與金州灣, 渤海中央以及遼東灣南部海區為低值區(圖3b); 夏季chla濃度高值區主要集中于黃河與遼河入?？诟浇Ｓ? 中央海區濃度較春季增加一倍(圖3d); 秋季由于光照、溫度等因素影響, chla濃度急劇下降, 僅在遼東灣沿岸出現高值區(圖3f); 冬季整個渤海的表層chla濃度繼續下降, 只在渤海灣沿岸與秦皇島沿岸地區濃度較高(圖3h)。

在季節變化上, 2019 年渤海夏季chla濃度顯著高于其他季節, chla與DIN 和DIP 濃度的季節性分布趨勢呈負相關, 與DON 濃度季節性分布趨勢呈正相關(圖3)。主要的濃度高值區分布在黃河口臨近海域,構建的模型很好的體現了黃河調水調沙引起的營養鹽輸入通量變化對浮游植物生物量平面分布產生的影響(Dingetal, 2020)。

綜上, DIN、DIP、DON 與chla濃度模擬值與實測值的相對標準偏差、相似性系數和Kappa 系數平均值分別為24%、0.77 和0.60, 說明通過該模型能夠較真實的反映渤海海域DIN、DIP、DON 與chla濃度分布情況, 可以用于渤海DIN、DIP、DON 與chla濃度長期變化模擬計算和數據重構。

2.1.3 年際變化 以1980 年至2020 年部分年代調查監測的DIN、DIP、DON 和chla濃度年際變化對模型進行驗證。圖4 中模擬值為模型對全渤海逐月模擬結果的平均值, 監測值為對應某一月份的所有監測值的平均值。模擬結果與調查監測結果比較表明,模擬結果在數值大小與變化趨勢上均能夠很好的再現渤海近40 年營養鹽和chla的年際變化, DIN、DIP、DON 和chla年際變化模擬結果與調查結果的RSD平均值小于20%, SI 平均值高于0.80(圖4)。

模擬結果表明, 20 世紀80 年代至90 年代中期DIN 濃度顯著增加, 1995 年冬季DIN 濃度達到最大值(24.4 μmol/L), 在1996 年至2004 年間, DIN 濃度逐漸下降, 在2005 年至2013 年期間, DIN 濃度再次呈現上升趨勢, 近幾年DIN 濃度又呈下降趨勢(圖4a), 構建的模型能很好地再現這兩次DIN 濃度峰值變化(王修林等, 2006; Wangetal, 2019)。分析表明, 農業化肥使用、污水排放等人類活動導致總氮河流輸入和大氣沉降通量增加, 是造成渤海DIN 濃度升高的主要原因, 自2002 年起水利部黃河委員會于夏季實施的調水調沙工程, 是導致渤海夏秋季營養鹽輸入和濃度增加的原因(Dingetal, 2020)。DIP 濃度在20 世紀80年代至90 年代初逐漸下降, 在1992 年夏季達到最小值(0.1 μmol/L), 然后增加至90 年代后期, 在1998 年冬季DIP 濃度達到最大值(1.1 μmol/L), 之后呈下降趨勢(圖4b), 模擬結果在年際變化和季節變化上均與文獻結果基本一致(Ningetal, 2010; Xinetal, 2019)。20 世紀90 年代后期DIP 濃度的下降, 與含磷洗滌劑的禁用有關。在20 世紀80 年代到90 年代初DON 濃度總體處于緩慢上升趨勢, 直至1994 年濃度有所下降;21 世紀10 年代末有所上升, 在2014 年又降低, 之后呈上升趨勢(圖4c), 這與陸源長期輸入變化與海源浮游植物、微生物長期演變過程有關(唐永等, 2017)。

Chla濃度自20 世紀80 年代至21 世紀初呈波動緩慢升高趨勢, 在年際變化上與營養鹽濃度, 特別是與DIN 濃度變化呈現很好的相關性, 這與孟慶輝等(2022)的研究結果一致。Chla濃度季節變化主要是夏季峰(圖4d), 與DIN 濃度呈反向變化關系, 這與夏季藻華對營養物質的消耗和冬季浮游植物死亡營養鹽恢復的生態學規律一致; 與DON 濃度呈復雜源匯變化關系, 說明DON 濃度主要受陸源和海源的影響,也是支持浮游植物生長的重要物質(唐永等, 2017)。

2.2 模擬結果準確性評價

為驗證模型改進后對模擬計算結果準確性的提高程度, 將改進前后的模擬計算量化驗證結果進行比較, 結果如圖5 所示。氮、磷營養鹽與chla濃度的相似性系數平均值分別提高了0.07 和0.13, Kappa系數分別提高了0.09 和0.10, 說明改進后的模型對氮、磷營養鹽與chla濃度的變化趨勢的模擬要顯著比模型改進前與實測結果更加吻合。二者的相對標準偏差也均有明顯降低, 氮、磷營養鹽與chla濃度分別降低了約11%和6%, 說明在數值大小上, 模型改進后的結果與實測結果也更為接近。綜上所述, 通過參數優化和優質數據選取對模型進行改進, 提高了模型模擬計算的準確性。

圖5 模型改進前后渤海氮、磷營養鹽與chl a 濃度模擬計算評價結果的Taylor 圖Fig.5 The Taylor diagram of the simulation results of N-P nutrients and chl a concentrations in Bohai Sea before and after model improvement

2.3 渤海富營養化評估

2.3.1 渤海富營養化狀態年際變化 渤海富營養化狀態指數(NQI)年際變化計算結果表明, 自1980 年至2020 年, NQI 整體呈先上升后下降的倒“U”型變化特征(圖6)。其中, 1980 年至1995 年, NQI 上升達到極大值2.27, 1995 年至2013 年NQI 緩慢下降, 2013 年至2020 年迅速下降; 從數值上看, 1993 年至2017 年, NQI基本都大于2, 富營養狀態處于中等水平, 當前處于貧富營養狀態。利用渤海復合富營養化指數(CEI)(Linetal,2020)年際變化計算結果進行后驗, 結果表明, NQI 與CEI 在變化趨勢上相當一致(SI=0.92) (圖6)。

圖6 1980 至2020 年渤海復合富營養化指數(CEI)與富營養化狀態(NQI)年際變化圖Fig.6 Interannual changes of Bohai composite eutrophication index (CEI) and eutrophication status (NQI) from 1980 to 2020

2.3.2 渤海富營養化狀態季節變化 在渤海全海域富營養化狀態季節變化規律上, NQI 呈倒“U”型變化特征, 夏季富營養化程度最高; NQI 在5～10 月份超過2.0, 處于中等富營養狀態。利用渤海復合富營養化指數(CEI)(Linetal, 2020)季節變化結果對比后驗,NQI 與CEI 在變化規律上也較為一致(SI=0.74)(圖7)。

圖7 渤海復合富營養化指數(CEI)與富營養化狀態(NQI)季節變化圖Fig.7 Seasonal change in composite eutrophication index (CEI)and eutrophication status (NQI) in Bohai Sea

2.3.3 渤海富營養化狀態平面分布 渤海富營養化狀態指數(NQI)平面分布模擬計算結果如圖8 所示。渤海富營養化狀態平面分布總體呈現近岸較高、中部較低的趨勢, 這與Lin 等(2020)的研究結果一致。近岸較高的富營養化狀態主要與人類活動有關, 且海灣近岸地區水體交換能力差, 二者綜合導致近岸海域富營養化程度顯著高于中部海域。在季節變化上,富營養化狀態與chla濃度模擬結果(圖3)變化規律基本一致, 即隨chla濃度春夏季升高, 富營養化水平升高, 隨chla濃度秋冬季下降, 富營養化水平下降,同時還受海洋中的環境因子(光照、溫度、營養鹽、混合層深度等)季節變化調控。

圖8 2019 年渤海富營養化狀態(NQI)平面分布圖Fig.8 Horizontal distribution of eutrophication status (NQI) in 2019 in the Bohai Sea

3 結論

本文通過改進生物地球化學模型架構, 建立完善了用于重構富營養化評價要素長期演變數據的渤海三維水動力/生物地球化學耦合數值模型。利用構建的模型成功模擬了渤海近四十年來DIN、DIP、DON與chla濃度時空分布[(相對標準偏差平均值)為24%,(余弦相似性系數平均值)為0.77,(Kappa系數平均值)為0.60], 說明模型相關動力學方程與參數能夠反映渤海海域的地域化特征。應用模型模擬結果計算了營養狀態質量指數(NQI), 與復合富營養化指數(CEI)具有較好的一致性(=0.83), 說明該模型能夠反映渤海富營養化狀況。

本文構建的三維水動力/生物地球化學耦合模型,對富營養化狀態的模擬計算結果能夠較真實的反映渤海海域富營養化狀態的時空變化情況, 建立的富營養化評估方法具有較高可靠性, 可以用于渤海富營養化評估, 并支撐渤海富營養化綜合治理。