附錄
本文所建立的廈門灣水動力模型模擬區(qū)域的總海域面積約為1 074 km2。海岸線數據采用美國海洋和大氣管理局(NOAA)的GSHHG數據庫, 數據集的坐標為WGS-84。水深數據來自于美國地球物理中心(NGDC) 2008年發(fā)布的ETOPO1數據庫, 數據集的坐標同樣為WGS-84。為保證模型的精度, 廈門灣近岸區(qū)域和九龍江河口處水深和地形采用小比例尺的海圖數據進行修正, 海圖選取2017年中華人民海事局發(fā)布的圍頭港-廈門港的6517100號海圖, 比例尺大小為1︰75 000。
由于地表水模擬軟件(SMS)中網格編輯器生成的網格和MIKE相比更加規(guī)整, 對于海岸等不規(guī)則邊界的適應性更強, 可以提高模型運算速度, 且其自帶的網格質量驗證功能可以很好地提升網格質量。因此, 將海岸線和水深數據在SMS軟件中預處理, 生成MIKE的網格輸入文件, 調整質量后的網格如附圖1所示。模型水平方向使用漸變的非結構化三角網格, 在近岸和九龍江河口處進行網格加密, 共生成10 644個節(jié)點和19 975個網格, 垂向上使用坐標等比例劃分成5層, 網格的水深值根據插入水深散點的距離進行線性插值, 最后將制作好的地形文件導入MIKE中進行前處理。
附圖1 廈門灣海域計算網格
Supplementary Fig.1 The computational grid of the Xiamen Bay
東部和南部外海開邊界條件設置為潮位強迫條件, 開邊界網格節(jié)點的潮位數據通過13個主要分潮(M2、S2、N2、K2、K1、O1、Q1、P1、MF、MM、M4、MS4、MN4)的調和常數計算得到。而開邊界調和常數的值則取自俄勒岡州立大學(OTIS)潮汐資料中的TPXO7.2數據集。TPXO模型基于二維正壓流體方程, 運用廣義反演法進行實測數據同化, 并采用最小二乘法進行數據擬合。該模型同化數據包括衛(wèi)星測高數據(T/P、Topex Tandem、ERS、GFO)和實測驗潮數據, 模型預測范圍覆蓋全球海平面, 模型的結果已經廣泛應用于我國的沿海地區(qū)(范文藍等, 2018; Wang, 2019)。本文采用TMD (Tide Modal Driver)工具箱計算得到開邊界網格節(jié)點每隔一小時的預報數據, 其他時間的潮位數據通過MIKE模型線性插值確定。同時考慮到九龍江河口徑流會攜帶大量來自上游的微塑料垃圾, 根據駱智斌等(2008)對九龍江河口的水動力潮流潮汐模型研究結果, 在九龍江邊界加入多年平均徑流量371 m3/s作為流量控制邊界, 陸地邊界設置為閉邊界, 即將法向流速設置為0。
風場數據則來自于國家海洋科學數據中心提供的中國臺站觀測數據集, 將廈門測點站每隔1 h的實測數據應用于整個模型區(qū)域, 規(guī)定風摩擦系數為0.001 255。為避免模型產生大的波動, 設置軟啟動間隔為7 200 s, 在此期間內風速從零增加到指定值。在計算中使用可變的時間步長間隔, 使收斂條件判斷數(Courant Friedrich Levy, CFL)在所有計算節(jié)點中小于臨界CFL值。規(guī)定臨界CFL值為0.8, 最小時間步長為0.01 s, 最大時間步長為30 s。
模型通過設置干濕臨界水深和淹沒臨界水深處理移動邊界的問題, 每個單元根據水深值和單元表面是否淹沒分成干燥、部分干燥和濕三種不同的狀態(tài), 其中干燥單元將會在計算中被移除, 部分干燥單元忽略其動量通量, 只計算質量通量。濕單元則計算動量通量和質量通量。本文規(guī)定干水深臨界值為0.005 m, 濕水深臨界值為0.1 m, 淹沒水深臨界值為0.05 m。
在近海潮波模擬中, 測站點多分布在淺水區(qū), 水深一般都不大, 在這種情況下, 底部的摩擦效應較深水而言要重要得多, 本文通過粗糙高度設置底床糙率, 不同位置的粗糙高度取值不同, 范圍為0.01~ 0.05 m。另外, 模型還考慮了渦流黏度和地球自轉產生的科式力的影響。水平渦流黏度通過Smagorinsky公式指定, 規(guī)定Smagorinsky系數為0.28, 渦流黏度的最小值和最大值分別為0.000 018和10 000 000 000 m2/s。垂直渦流黏度選擇-模型計算, 最小值和最大值分別為0.000 018和0.4m2/s。
為了使模型更加平滑地啟動, 本文將水動力模型設置為冷啟動, 將模型初始運行時刻外海開邊界潮位的平均值0.62 m作為初始水位, 初始流速設置為0。
本文將海洋中懸浮塑料顆粒的物理運輸過程簡化為以下四個過程:
(1) 塑料顆粒的大氣沉降: 風會攜帶塑料顆粒進入大氣中, 大氣到海洋表面的塑料顆粒沉降十分常見; (2) 對流擴散: 通過平流和擴散在水柱中運輸塑料顆粒; (3) 沉積: 密度大于海水的塑料顆粒會從水柱中沉積到底床表面; (4) 再懸浮: 當海水流速大于海底沉積塑料顆粒的再懸浮臨界速度時, 塑料顆粒會重新懸浮到水柱中。
將懸浮塑料顆粒(1)、(3)、(4)的三個基本運輸過程在MIKE生物仿真模塊中通過開放式的數學公式定義, 并耦合對流擴散模型則可對懸浮塑料顆粒的運動過程進行較為完整的模擬。
在ECO Lab創(chuàng)建懸浮塑料顆粒遷移模板, 分別定義狀態(tài)變量(state variables, 描述生物個體包括微塑料顆粒對水環(huán)境參數做出的響應的物理量, 如懸浮微塑料的濃度等)、常數(constants, 不隨時間變化的物理量, 如懸浮微塑料的沉降速度等)、過程(processes, 與微塑料遷移運動相關的物理過程, 如前文所述的大氣沉降、對流擴散等)、作用力(forcing, 外力作用, 如水動力模塊輸出的流場、風場信息)、公式表達(expression)、附加輸出(derived outputs)等部分。其中, 塑料顆粒模板定義的狀態(tài)變量為懸浮塑料顆粒的濃度值和沉積塑料顆粒的面密度。記懸浮塑料顆粒的濃度值為MPSS, 沉積物塑料顆粒的面密度為MPSED, 則懸浮塑料顆粒濃度變化可由以下常微分方程描述:
dMPSS/d= prss–sessv+ressv, (B-1)
其中, MPSS為懸浮塑料顆粒的濃度, 單位為mg/L; d為時間步長, 單位為d; prss、sessv、ressv分別為大氣沉降過程塑料顆粒輸入濃度變化率, 懸浮塑料顆粒的沉積率和再懸浮率, 單位均為g/(m3·d)。
式(B-1)等號右邊分別為塑料顆粒的大氣沉降過程、懸浮塑料顆粒的沉積過程和沉積物中塑料顆粒的再懸浮過程的總和, 這三個過程具體的公式表達如下:
prss = papro/d
sessv = vsm×MPSS×86 400 (B-2)
ressv = ressa/d
式中, papro為大氣沉降到海水表面的塑料顆粒沉降速率, 參考Eco Lab關于污染物擴散模型的設置, 假定其主要與空氣中塑料顆粒的濃度有關, 其與風速風向的相關性則需要開展進一步研究; sessv為單位體積沉降速率; d為垂直方向上劃分的層高; vsm為懸浮塑料顆粒的沉降速度, 與塑料顆粒的相對粒徑, 形狀大小(球狀或碎片等)等有關; 本文將沉降過程簡化為當懸浮物濃度達到一定值時, 塑料粒子開始向海底沉降, 其濃度臨界值可以根據經驗或者實測值決定。
ressa為單位面積塑料顆粒懸浮率, 可由下式邏輯表達式求得,
ressa = If bssmx×bresu==1
Then min(resrat, (Max(0, XSED–1))/d×86 400)(B-3)
該邏輯表達式表示了沉積物中塑料顆粒再懸浮的條件既要滿足水柱中懸浮塑料顆粒達到飽和濃度, 同時也要滿足流速大于再懸浮的臨界流速。其中, resrat為最大再懸浮率; bssmx表示水中懸浮是否未飽和, 其中1為飽和, 0為未飽和; bresu表示是否懸浮, 1為懸浮, 0為非懸浮。塑料顆粒的重懸浮率與水流大小及沉積塑料顆粒的沉降速度有關, 相同材質的塑料, 粒徑大的更容易懸浮起來, 因為粒徑大的塑料粒子受到的浮力更大, 更容易在水流的作用下產生擾動。注: 本文將再懸浮過程進行簡化, 只考慮了交換過程中再懸浮量的大小, 而未考慮水流的流向等的影響。
對于沿岸城市向海洋排放塑料垃圾的研究, 本文綜合分析后采用陸源輸入、海濱旅游、排污口排污和船舶垃圾四個參數作為海洋微塑料垃圾的來源, 將廈門灣主要的排污口處、人為活動區(qū)域密集的廈門島西海域和九龍江徑流處設置為微塑料排放的點源。同時在廈門北航道和廈門島近岸分別設置移動點源和面源以考慮濱海旅游和船舶垃圾入海的效應。其中各點源的位置如附圖2所示。
附圖2 廈門灣海洋微塑料垃圾入海點源
Supplementary Fig.2 The sources of microplastic around the Xiamen Bay