王仁政, 單正垛, 孟思雨, 宮響
1. 青島科技大學數(shù)理學院, 山東 青島 266061;
2. 中國海洋大學環(huán)境科學與工程學院, 山東 青島 266100
次表層葉綠素最大值(subsurface chlorophyll maxima, SCMs)是全球海洋生物剖面的顯著共性之一,是層化水體中葉綠素濃度最大值發(fā)生在次表層的一種現(xiàn)象。在熱帶、亞熱帶海域SCMs 所貢獻的初級生產力約占水體總初級生產力的30%~70%(Takahashi et al, 1984), 并且對海洋新生產力也有很大貢獻。Hodges 等(2004)研究發(fā)現(xiàn)在SCMs 發(fā)生條件下新生產力的f比(新生產力在總初級生產力中的比例)增大了一倍, Hanson 等(2007)在澳大利亞附近寡營養(yǎng)海域發(fā)現(xiàn)這一最大值層中f比約為7%~32%。而海洋新生產力反映了海洋真光層對大氣中CO2的凈吸收能力。因此研究SCMs 現(xiàn)象, 特別是其年際變化對深入理解海洋生態(tài)系統(tǒng)在氣候調節(jié)中的作用具有重要意義。
中國南海北部海域擁有十分復雜的環(huán)流體系、生態(tài)系統(tǒng)和水文特征(樂鳳鳳 等, 2006), 導致南海北部海域次表層葉綠素時空分布特征的多樣性。眾多學者開展了南海北部葉綠素垂直分布的研究, 如,張釩 等(2000)基于觀測數(shù)據(jù)分析了1997 年8 月臺灣海峽南部SCMs 特征因子的周日變化; 柯志新等(2013)研究了2008 年8、9 月份南海北部SCMs 的空間分布特征及其影響因素; Lu 等(2010)模擬研究了南海北部夏季SCMs 深度、強度與上升流、生物平衡的關系; Gong 等(2014)基于數(shù)值模擬研究了南海北部海盆區(qū)SCMs 的季節(jié)變化。
盡管對南海北部海域SCMs 現(xiàn)象已有一定認識,但由于現(xiàn)場觀測數(shù)據(jù)不足且SCMs 發(fā)生在遙感可探測深度以外, 目前對SCMs 年際變化的研究較少。數(shù)值模擬是研究SCMs 的重要工具。本文采用一維物理—生態(tài)耦合模型模擬南海北部海盆區(qū)1994—2019 年海溫、鹽度、葉綠素等要素, 并基于模型結果, 采用3 種統(tǒng)計方法(創(chuàng)新趨勢分析、經驗模態(tài)分解和Mann-Kendall 突變檢驗)分別給出總體趨勢、不同時間尺度變化過程及突變節(jié)點, 探討了海水溫度與SCMs 特征因子的關系及影響機制, 最后利用時間序列模型對2020 年3—11 月SCMs 特征因子的變化趨勢進行預測。
本文采用的數(shù)值模型為一維物理-生態(tài)耦合模型(Modular Ecosystem Model-1D, MEM-1D)(Ricci et al, 2000)。生態(tài)模式采用的是歐洲區(qū)域性海洋生態(tài)系統(tǒng)模型(European Regional Sea Ecosystem Model,ERSEM), 詳細過程見圖1; 該模型已被成功應用于我國黃海與南海北部海洋生態(tài)系統(tǒng)的研究中(張沖等, 2011; Gong et al, 2014; Butensch?n et al, 2016;Pan et al, 2017)。
圖1 European Regional Sea Ecosystem Model 食物網(wǎng)結構Fig. 1 Schematic of the trophic structure of European Regional Sea Ecosystem Model
物理模式采用的是一維普林斯頓海洋模型(Princeton Ocean Model, POM), 其中湍混合模型采用傳統(tǒng)的 2 階 Mellor 和 Yamada 湍流閉合模型(Mellor et al, 1982), 并采用夏潔等(2006)研究中關于POM 模式的垂直渦動動量Kh和熱混合系數(shù)Km的改進方案, 見公式1、2。
式中:W為風速;k為卡曼常數(shù), 取0.4;δ為波陡, 取0.1;β為波齡, 取1.0; 重力加速度g取9.8m·s-2;P為與Richardson 數(shù)有關的無量綱系數(shù), 取0.1;Z為水深;Kh′、Km′為考慮海浪作用后的垂直混合系數(shù)。
模擬區(qū)域為南海北部海盆區(qū)(南海18°N 以北且水深大于1000m 的區(qū)域)。模型水深取南海的平均深度1200m, 垂直方向上劃分為上密下疏的40 層。時間步長為216s。模型模擬時間為1990—2019 年(其中1990—1993 年為模型spin-up 時間)。模型采用南海北部SEATS(South East Asian Time-series Study)站(18°N, 116°E)的溫度、鹽度、營養(yǎng)鹽濃度和葉綠素濃度多年平均觀測數(shù)據(jù)作為初始條件, 外強迫主要考慮風與太陽輻射, 均采用 SEATS 站附近1990—2019 年數(shù)據(jù), 其中海表風速、凈長波輻射通量與凈短波輻射通量來自NCEP 再分析數(shù)據(jù)集, 潛熱通量與顯熱通量通過COARE 算法(Fairall et al,1996)計算得到。本文模型參數(shù)主要來自ERSEM 中的常用參數(shù), 部分選自南海及鄰近海域的研究結果,具體參數(shù)見表1—4。
表1 ERSEM 模型基本參數(shù)的取值Tab. 1 Values of basic parameters of ERSEM model
表2 浮游植物參數(shù)的取值Tab. 2 Values of phytoplankton parameters
表3 浮游動物參數(shù)的取值Tab. 3 Values of zooplankton parameters
表4 細菌和碎屑參數(shù)的取值Tab. 4 Values of bacteria and debris parameters
本文采用多種統(tǒng)計方法, 從不同角度分析SCMs 特征因子的年際變化特征。
1) 創(chuàng)新趨勢分析(innovative trend analysis, ITA)是?en(2012)提出的一種新型趨勢檢驗方法。該方法的步驟為: 將一個時間序列等分為第一、第二兩個子序列(若序列長度為奇數(shù), 等分時不考慮第一項),將各子序列按照升序排列后將其分別放置于笛卡爾坐標系x,y軸, 通過觀察數(shù)據(jù)點在1∶1 線附近的分布來判斷數(shù)據(jù)趨勢。趨勢斜率nx可用于判斷變化趨勢, ITA 指數(shù)D(Wu et al, 2017)可用于判斷趨勢變化幅度(計算如式3)。
2) 經驗模態(tài)分解(empirical mode decomposition,EMD)由Huang 等(1998)提出, 被廣泛應用于趨勢項提取、海洋數(shù)據(jù)分析(姜浩 等, 2018)、降雨量時間序列分析(陳國鼎 等, 2018)等問題。該方法的實質是不斷地從數(shù)據(jù)中分解出高頻分量, 并把剩下的低頻分量作為新數(shù)據(jù)繼續(xù)分解的過程。被分解出的高頻分量被稱作本征模態(tài)函數(shù)(intrinsic mode function,IMF), 其包含了原數(shù)據(jù)在不同時間尺度上的局部特征, 而最后無法繼續(xù)分解的數(shù)據(jù)被稱為殘差(residual, RES), 通常情況下代表原始時間序列數(shù)據(jù)的趨勢項。原始時間序列數(shù)據(jù)x(n) 可以表示為:
式中ci(n)為數(shù)據(jù)在第i時間尺度上的本征模態(tài)函數(shù)(IMF),rt(n) 為數(shù)據(jù)分解后的殘差(RES)。
對于各不同時間尺度的本征模態(tài)函數(shù), 其變化周期T可通過計算 IMF 平均頻率ω(式 5)求得(Flandrin et al, 2004)。
式中:N為IMF 分量與x軸交點個數(shù),L為第一交點與最后交點之間的距離。T為ω的倒數(shù)。
3) Mann-Kendall(M-K)突變檢驗是世界氣象組織推薦并已被廣泛使用的時間序列分析方法, 可以用于判斷氣候序列中是否存在氣候突變, 如果存在,可確定出突變發(fā)生的時間。Mann-Kendall 檢驗法也經常用于氣候變化影響下的降水、干旱頻次趨勢檢測。該方法最初由H. B. Mann 和M. G. Kendall 提出(Mann, 1945; Kendall, 1975), 其原理為: 對于時間序列數(shù)據(jù)xn, 構造一個秩序列ri(i=2,3,…,n-1), 定義ri為: 針對xi所有滿足j≤i且xj<xi的數(shù)據(jù)xj的個數(shù), 由此定義統(tǒng)計量Sk:
在滿足時間序列獨立且隨機的前提下定義統(tǒng)計量(式7), 再將原序列數(shù)據(jù)逆序重新排列, 重復上述步驟構建新統(tǒng)計量。當時間序列為正序時, 統(tǒng)計量用UFk表示; 逆序時用UBk表示。根據(jù)UFk和UBk繪圖曲線來判斷數(shù)據(jù)趨勢變化與突變點
差分整合移動平均自回歸(auto regressive integrated moving average, ARIMA)模型是由Box 等(1970)提出的一種隨機時間序列預測分析方法。該方法利用外推機制描述時間序列的變化, 是最小方差意義下一種精度較高的時序短期預測方法。為了處理具有周期性或季節(jié)性時間序列數(shù)據(jù), 在ARIMA 模型的基礎上拓展出了季節(jié)性差分自回歸滑動平均(seasonal autoregressive integrated moving average, SARIMA)模型(Box 等, 1970)。SARIMA 模型是在假定季節(jié)相關和普通相關的交互作用下建立的乘法模型, 可以表示為SARIMA(p,d,q)×(P,D,Q)s。 模型一般形式如下式:
式中yt為時間序列,μt為隨機項;φp(B)、ΦP(BS)分別為非季節(jié)與季節(jié)的自回歸部分; (1-B)d、(1-Bs)D分別為非季節(jié)與季節(jié)的差分部分;θq(B)、ΘQ(BS)分別為非季節(jié)與季節(jié)的移動平均部分。
Lewis 等(1983)最早提出了通過高斯函數(shù)擬合將海洋中葉綠素濃度垂直分布數(shù)據(jù)參數(shù)化的方法,這為SCMs 的解析計算奠定了理論基礎。由于上混合層內葉綠素垂直分布較為均勻, 而SCMs 多發(fā)生在上混合層以深, 本文采用分段函數(shù)擬合真光層中葉綠素濃度ρ的垂直分布(Gong et al, 2014):式中h是葉綠素濃度在區(qū)間[zs,ze]的垂直積分值,z代表水深(非負值),zs與0ρ均為常數(shù), 分別為葉綠素濃度均勻分布處的最大水深和葉綠素濃度,σ為高斯分布的標準差,ze為真光層深度。SCMs 深度定義為葉綠素濃度最大值所處水深(zm)。SCMs 強度為次表層葉綠素濃度最大值(h/σ)。SCMs 厚度為SCMs 上下邊界(滿足 ?2ρ(z)/?z2= 0條件處)之間的水深(一般情況下為2σ, 當上混合層深度低于上邊界時, 計算時需減去兩者之間的差值)。
為了驗證MEM-1D 模型, 將模擬結果與南海SEATS 站觀測數(shù)據(jù)進行比較(如圖2—4)?;谀虾1辈亢S蚨維CMs 現(xiàn)象不明顯, 本文不分析南海北部12 月、1 月和2 月的SCMs 情況。
由圖2、3 可知, 溫度和鹽度的模擬數(shù)據(jù)整體與觀測數(shù)據(jù)較為吻合, 偏差主要集中在表層鹽度, 這可能與所用模式沒有考慮海表輸入有關。MEM-1D模式較好地再現(xiàn)了南海北部海盆區(qū)的SCMs 現(xiàn)象(圖4), 反映了SCMs 的季節(jié)變化和不同年份間差異, 特別是SCMs 深度與實際基本一致, 主要位于60~80m范圍內。但模式對SCMs 強度的模擬存在一定偏差,如2000 年7 月和10 月SCMs 強度較實際偏小約0.15mg·m-3, 而 2004 年春季則略偏大, 這可能與SCMs 在不同站位的差異有關。本文所采用的分段函數(shù)能有效擬合上混合層以深的葉綠素垂直分布(圖4), 得到SCMs 特征值。
圖2 海溫模擬觀測數(shù)據(jù)對比Fig. 2 Comparison of simulated and observed ocean temperature
圖3 鹽度模擬觀測數(shù)據(jù)對比Fig. 3 Comparison of simulated and observed ocean salinity
圖4 葉綠素濃度模擬觀測數(shù)據(jù)對比Fig. 4 Comparison of simulated and observed chlorophyll concentrations
2.2.1 創(chuàng)新趨勢分析
對1994—2019 年SCMs 特征因子數(shù)據(jù)進行創(chuàng)新趨勢分析(ITA), 結果如圖5。
由圖5 可知, 1994—2019 年期間SCMs 的深度和厚度均呈上升趨勢, 而強度的散點在y=x直線上下分布為6∶7, 呈較弱的減小趨勢; 值得注意的是圖中厚度的低值區(qū)表現(xiàn)出十分明顯的上升趨勢(高低區(qū)域以橫坐標平均線為界, 以垂直虛線表示)。特征因子的總趨勢與趨勢斜率(強度斜率為-9.5×10-5、深度斜率為3.5×10-3、厚度斜率為1.8×10-3)相吻合。SCMs 強度、深度和厚度的ITA 指數(shù)分別為0.03、0.05 和0.07, 故在1994—2019 年SCMs 厚度變化幅度最大, 深度次之, 強度最小。
圖5 SCMs 特征因子趨勢的創(chuàng)新趨勢分析Fig. 5 Trend of SCM characteristics by innovative trend analysis
2.2.2 經驗模態(tài)分解分析
據(jù)國家旅游部門統(tǒng)計,2016年我國出境游人數(shù)達到1.22億人次,旅游消費1098億美元,出境游客人次連續(xù)13年以兩位數(shù)的速度增長,這兩項指標連續(xù)四年位居全球第一,對全球旅游收入的貢獻年均超過13%。與之相伴,中國留學大軍也遍布五大洲,2015年總人數(shù)首次突破50萬,2016達到54.45萬人。中國也因此摘得世界第一大留學生輸出國的桂冠。
對1994—2019 年SCMs 特征因子數(shù)據(jù)進行經驗模態(tài)分解分析, 結果如圖6—8。
圖6 SCMs 強度變化的經驗模態(tài)分解Fig. 6 Variation of SCM intensity by empirical mode decomposition
由式5 計算可得, 強度IMF 分量的時間尺度分別為2、5、10、18、44 個月。圖6 中IMF3 近似體現(xiàn)了強度年際變化過程: 在1994—2006 年趨勢變化較為平緩, 2006—2011 年趨勢變化較大,2011—2018 年趨勢變化再次減緩。IMF5 近似體現(xiàn)了強度以9 年為周期的年代際變化過程。由RES部分可知, SCMs 強度在1994—2004 年有一定幅度的下降, 在 2004—2012 年趨勢逐漸上升, 在2012—2019 年趨勢再次下降。
由式5 計算可得, SCMs 深度的4 個IMF 分量時間尺度分別為4、11、16、48 個月。圖7 中IMF2近似體現(xiàn)了深度年際變化過程: 在1994—2004 和2009—2015 年呈波動變化, 在2004—2009 和2015—2019 年趨勢較為平緩。IMF4 近似體現(xiàn)了深度年代際變化過程: 在1994—2005 年不斷變深, 在2005—2019 年呈現(xiàn)周期波動。由RES 部分可知: SCMs深度先變深(1994—2011 年), 后變淺(2011—2019年)。
圖7 SCMs 深度變化的經驗模態(tài)分解Fig. 7 Variation of SCM depth by empirical mode decomposition
由式5 計算可得, 厚度IMF 分量的時間尺度分別約為2、4、10、25、43、115 個月。圖8 中IMF3 近似表現(xiàn)了厚度的年際變化過程: 1994—2009、2011—2019 年呈現(xiàn)較穩(wěn)定的周期波動, 但在2009—2011 年維持在較低水平。IMF5 近似表現(xiàn)了厚度年代際變化過程: 1994—2002 年不斷變大, 2002—2008 年逐漸變小, 2008—2019 年維持平穩(wěn)。由RES 部分可知, 1994—2019 年SCMs 厚度呈增大趨勢。
圖8 SCMs 厚度變化的經驗模態(tài)分解Fig. 8 Variation of SCM thickness by empirical mode decomposition
2.2.3 M-K 突變檢驗分析
對 1994—2019 年 SCMs 特征因子數(shù)據(jù)進行Mann-Kendall 突變檢驗, 結果如圖9。
圖9 SCMs 特征因子變化趨勢的M-K 法分析Fig. 9 Trend of SCM characteristics by M-K mutation test
由圖9 可知, 厚度在1996 年產生突變(突變點位于95%置信區(qū)間內), 整體呈上升趨勢(UF>0), 且于1997 年起上升趨勢顯著(通過顯著水平α=0.05 的檢驗); 強度于1994 年起點處存在交點, 經改變時間尺度(1994— 2003 年)后做M-K 突變檢驗發(fā)現(xiàn), 該交點不成為突變點, 1994—2019 年SCMs 強度呈下降趨勢(UF<0), 且1999—2004 年明顯減小; SCMs 深度變化不存在突變點, 整體呈上升趨勢(UF>0),1994—1997 年上升趨勢不明顯, 之后呈顯著上升趨勢。
綜合上述結果可以看出, 在 1994—2019 年SCMs 深度、厚度整體均為上升趨勢, 而強度下降趨勢不明顯; SCMs 強度在2006—2011 年有相對較大的減小趨勢, 深度在1994—2004 和2009—2015 年變化相對較大, 厚度在2009—2011 年維持在較低水平; 厚度在1996 年產生突變并于1997 年變?yōu)轱@著上升趨勢。
為理解氣候變化相關因素對SCMs的影響, 本文討論海水溫度(表面與深層海溫)與SCMs 特征因子的關系。年際變化上, 相關分析結果顯示, 海表面溫度SST 與 SCMs 三個特征因子均不存在相關關系(P>0.05)。由SST 與SCMs 三個特征因子標準化后的年際變化序列(圖10)可見, 在厄爾尼諾影響下, 1998年和2016 年南海北部海盆區(qū)SST 顯著升高, 但SCMs三個特征因子變化較小。這表明氣候變化導致的海表面升溫可能對次表層浮游植物的生長影響不大。
圖10 SCMs 特征因子與海表溫度標準化值的年際變化Fig. 10 Interannual variation of standardized value of SCM characteristics and SST
月份變化上, SCMs 深度、強度均與海表面溫度SST 呈顯著的正相關關系(r分別為 0.79、0.49,P<0.001, 圖11), 而SCMs 厚度與SST 不存在相關關系(P>0.05)。特別是春季和秋季, 隨 SST 增加,SCMs 深度顯著變深(r分別為0.86、0.81,P<0.001,表1), SCMs 強度也顯著變大(r分別為0.51, 0.69,P<0.001, 表5)。這可能是由于海表溫度上升, 次表層海水密度變小(Midhun shah et al, 2020), 更多表層浮游植物沉降到次表層(Fennel et al, 2003), 導致SCMs 強度變大, 深度變深; 同時由于海表升溫使得密度分層得到加強, 上混合層深度變淺, 營養(yǎng)鹽垂直向上混合層的輸送減少, 表層浮游植物生長受到限制(Lewandowska et al, 2014; Meng et al, 2021),這使得次表層水體中光合有效輻射增強, 因此SCMs 深度和強度隨之增加。夏季SCMs 強度和深度與SST 不存在顯著相關性或相關性較小, 這可能是由于夏季南海北部海盆區(qū)溫度較高, 當溫度超過一定范圍時, 對浮游植物生長的影響較小, 同時由于夏季光輻射較強, 從而導致光抑制作用。而SCMs厚度主要受浮游植物種群結構、生物特性及次表層水體湍混合強度等因素影響(Beckmann et al, 2007;Gong et al, 2015), 與SST 不存在相關性(Midhun shah et al, 2020)。
表5 SST 與SCMs 特征因子的相關性Tab. 5 Correlations between SST and SCM characteristics
圖11 海表面溫度與SCMs 強度、深度的線性相關關系Fig. 11 Linear correlations of SST with the intensity and depth of SCMs
為了進一步探討海水溫度與SCMs 特征因子的相關關系, 選取深層海溫(SCMs 上邊界處、深度處和下邊界處海水溫度, 各處水深取歷年平均值為59、75和91m)與特征因子進行相關性分析, 結果見表6。
表6 深層海溫與SCMs 特征因子的相關性Tab. 6 Correlations between deep-sea temperature and SCM characteristics
由表6 可知, SCMs 特征因子僅與SCMs 上邊界處海水溫度存在顯著相關性(P<0.05), 但相關系數(shù)較小, 強度與其呈負相關, 深度和厚度與其呈正相關, 這表明混合層以深的水溫對SCMs 影響不大。南海北部海盆區(qū)混合層深度除冬季外均較淺, 約20~40m(Tseng et al, 2005), 而葉綠素峰值多發(fā)生60~80m, SCMs 厚度約30m, 其上邊界多位于混合層以深15m 處, 上邊界處海溫增加, 意味著上混合層對營養(yǎng)鹽的向上卷挾較多, 次表層浮游植物可在更大范圍內快速生長, SCMs 厚度變大, 峰值減小(Beckman et al, 2007)。
本文首先利用1994-2018 年SCMs 特征因子的月平均值, 根據(jù)貝葉斯信息準則(Bayesian information criterion, BIC)得到SCMs 特征因子的最優(yōu)SARIMA模型, 結果顯示: SCMs 強度 SARIMA 模型為
(0,0,1)×(0,1,1)9、深度SARIMA 模型為(0,0,1)×(0,1,1)9、厚度為SARIMA(0,0,0)×(0,1,1)9。
利用所得的SARIMA 模型預測2019 年SCMs特征因子變化趨勢, 并采用平均絕對百分比誤差(mean absolute percentage error, MAPE)這一指標來評估模型結果(式 10), 可得 SCMs 特征因子的MAPE 分別為5.33%(強度)、0.62%(深度)、2.49%(厚度), 表明SARIMA 模型可用于SCMs 特征因子的預測分析。
式中:T為預測時間,為t時刻的預測值,ty為t時刻的實際值。
為預測2020 年3—11 月SCMs 特征因子的變化情況, 本文利用1994—2019 年SCMs 特征因子的月平均值, 再次根據(jù)貝葉斯信息準則確定三個特征因子的最優(yōu)SARIMA 模型(與1994—2018 年模型相同)。預測結果如圖12。如圖12 所示, 2020 年3—11月SCMs 強度、深度和厚度的預測值與歷史值相差不大, 其極大值與極小值發(fā)生的時間基本與歷史同期, 表明SARIMA 模型可用于對SCMs 特征因子的預測。
圖12 1994-2020 年SCMs 強度、深度與厚度的回歸擬合及預測Fig. 12 Regression fitting and prediction of the intensity, depth and thickness of SCMs from 1994 to 2020
本文利用MEM-1D 模型模擬了1994—2019 年南海北部海盆區(qū)SCMs 現(xiàn)象, 并對模擬結果進行了統(tǒng)計分析, 得到了SCMs 特征因子的年際變化特征,并探討了SCMs 特征因子變化與海水溫度的關系。
在1994—2019 年, 深度、厚度整體均為上升趨勢而強度為下降趨勢, 厚度在1996 年開始突變并于1997 年變?yōu)轱@著上升趨勢。具體表現(xiàn)為, SCMs 強度在1994—2004 年逐漸變小, 2004—2012 年不斷增加,而在2012—2019 年強度再次減小; 深度在1994—2011 年不斷變深, 在2011—2019 年略微變淺; SCMs厚度在1994—2019 年呈持續(xù)上升趨勢。年際變化上,海表面溫度與SCMs 特征因子間不存在相關性, 表明在南海北部海盆區(qū), 氣候變化導致的海表面升溫對SCM 層的浮游植物生長影響不大。但在月份變化上, SCMs 強度和深度與海表溫度之間均為顯著正相關關系, 特別在春季和秋季, 當海表溫度上升時SCMs 強度增加、深度變深。深層海溫與SCMs 特征因子基本不具有相關性。最后采用SARIMA 模型對2020 年3—11 月的SCMs 特征因子變化做出了預測, 模型評估誤差小于10%且預測值與歷史同期相近, 表明SARIMA 模型可用于預測SCMs 特征因子變化趨勢。