二維赤潮生態(tài)數(shù)學(xué)模型及其在渤海的應(yīng)用

李大鳴, 林 毅, 宋雙霞, 解以揚(yáng), 韓素琴, 李培彥

(1. 天津大學(xué) 建筑工程學(xué)院暨港口與海洋工程教育部重點(diǎn)實(shí)驗(yàn)室, 天津300072; 2. 天津市氣象科學(xué)研究所, 天津300074)

二維赤潮生態(tài)數(shù)學(xué)模型及其在渤海的應(yīng)用

李大鳴1, 林 毅1, 宋雙霞1, 解以揚(yáng)2, 韓素琴2, 李培彥2

(1. 天津大學(xué) 建筑工程學(xué)院暨港口與海洋工程教育部重點(diǎn)實(shí)驗(yàn)室, 天津300072; 2. 天津市氣象科學(xué)研究所, 天津300074)

采用ADI的有限差分方法對不可壓縮流體二維淺水環(huán)流方程離散和求解, 建立水動力數(shù)學(xué)模型,用迎風(fēng)格式離散赤潮生物動力學(xué)方程, 通過水動力學(xué)和生物動力學(xué)相結(jié)合的方法, 建立了二維赤潮生態(tài)數(shù)學(xué)模型。將所建立的二維赤潮生態(tài)數(shù)學(xué)模型應(yīng)用于渤海, 針對渤海海域2004年6月11～16日發(fā)生的棕囊藻赤潮進(jìn)行了數(shù)值計(jì)算。對EOS/MODIS衛(wèi)星拍攝的2004年6月份的渤海海區(qū)衛(wèi)星遙感圖像進(jìn)行了處理, 提取出海水中的赤潮信息, 并計(jì)算出赤潮面積, 使其與模型計(jì)算出的赤潮面積進(jìn)行對比驗(yàn)證, 結(jié)果基本吻合, 表明該模型能夠較好地模擬赤潮的生消過程, 為渤海地區(qū)的赤潮預(yù)報(bào)提供了科學(xué)依據(jù)。

二維赤潮生態(tài)數(shù)學(xué)模型; ADI法; 迎風(fēng)格式; 棕囊藻; 渤海; 衛(wèi)星遙感

對赤潮生態(tài)的研究迄今為止已有30多年的歷史,從赤潮生物、生理、生化研究到赤潮生物發(fā)生機(jī)制研究以及赤潮預(yù)測預(yù)報(bào)研究, 所涉及的領(lǐng)域越來越深入, 要求也越來越高。越來越多的學(xué)者從赤潮生物的生理、生化特性、毒理學(xué)機(jī)制和發(fā)生機(jī)制、赤潮生物的生態(tài)環(huán)境以及赤潮的預(yù)測預(yù)報(bào)等方面對赤潮生物、赤潮的種類、赤潮發(fā)生與環(huán)境的關(guān)系展開了較為全面的探索。但是, 總的看來, 至今的赤潮研究主要還是集中在赤潮的生理、生態(tài)方面的研究, 因此要達(dá)到預(yù)報(bào)赤潮的發(fā)生還有漫長道路。

生態(tài)數(shù)學(xué)模型是赤潮研究的重要方法之一, 可以通過參數(shù)控制進(jìn)行數(shù)值試驗(yàn), 以便確定哪些因素在赤潮發(fā)生和發(fā)展的過程中起主要作用,既可以確定單個因素的影響, 又能對多因素的效果進(jìn)行研究。齊雨藻、黃偉健等[1]對基礎(chǔ)的統(tǒng)計(jì)學(xué)模型進(jìn)行赤潮研究方面的工作。以微分方程為主要工具的赤潮生態(tài)動力學(xué)模型的研究工作始于Chen[2]、Ditoro[3]開發(fā)的簡單的水質(zhì)動力模型。在我國, 王壽春[4]等建立了這種生態(tài)動力學(xué)模型, 以種群生態(tài)學(xué)和營養(yǎng)動力學(xué)為依據(jù), 提出了赤潮生物動力學(xué)研究對象之間的基本關(guān)系。但這些模型都沒有考慮環(huán)境要素的影響, 特別是水動力環(huán)境的影響。

由于赤潮的成因相當(dāng)復(fù)雜, 除了赤潮生物自身的特性外, 還涉及到物理、化學(xué)、水文、氣象等諸多因素。為了對赤潮問題作更深入的研究, 必須借助動力學(xué)方法綜合考慮各種過程以便從整體上探討赤潮發(fā)生、發(fā)展的規(guī)律。如夏綜萬等[5]通過海洋動力學(xué)和赤潮生物動力學(xué)相結(jié)合的方法建立的大鵬灣夜光藻赤潮生態(tài)仿真模型。計(jì)算采取了ADI的方法,模型包括水動力、擴(kuò)散和生物動力學(xué)三部分, 綜合考慮了潮流、營養(yǎng)物質(zhì)等環(huán)境要素的時(shí)空變化對赤潮過程的影響, 并以大鵬灣夜光藻赤潮為例進(jìn)行了數(shù)值模擬。田峰等[6]針對近岸海域赤潮藻類生長及分布的特點(diǎn),將簡化了的赤潮藻類模型與水動力學(xué)中的對流擴(kuò)散方程相耦合, 建立了一個水動力與生態(tài)耦合的赤潮藻類生長的深度模型。

要實(shí)現(xiàn)可靠的赤潮預(yù)報(bào), 就必須將各種動力過程相結(jié)合, 從動態(tài)方面去闡明海洋赤潮高發(fā)區(qū)環(huán)境系統(tǒng)內(nèi)赤潮生物生態(tài)過程、營養(yǎng)物質(zhì)生物化學(xué)過程、赤潮生消的關(guān)鍵物理過程之間的定量關(guān)系等。這些只能通過數(shù)值模擬來實(shí)現(xiàn)。所以本文通過建立赤潮生態(tài)數(shù)學(xué)模型, 將水動力過程與生態(tài)動力過程相結(jié)合, 使較為可靠的赤潮預(yù)報(bào)成為可能。

1 水動力數(shù)學(xué)模型的建立

1.1 水動力數(shù)學(xué)模型的基本方程

1.1.1 基本控制方程

由 N-S方程推導(dǎo)出的不可壓縮流體二維淺水環(huán)流方程為：

式中, 為增水位;Hhξ=+,h為水深;u為x方向的平均流速;v為y方向的平均流速;g為重力加速度;為謝才系數(shù), 其中h是計(jì)算網(wǎng)格圍點(diǎn)的水深平均值,n是粗糙系數(shù), 即糙率, 由試驗(yàn)選取;f為柯氏力系數(shù);HA為水平渦黏系數(shù)。

1.1.2 定解條件

岸邊界條件：vn= 0 (n是邊界法線方向)

ξ=ξ*(ξ*是已知邊界的水位過程)

初始條件： 當(dāng)t=0時(shí),ξ=0,u=v=0

1.2 ADI的有限差分方法

本文對潮流場的數(shù)值模擬采用 ADI的差分方法(Alternating Direction Implicit Method)。它是一種顯隱交替使用的有限差分格式, 其特點(diǎn)是將一個時(shí)間步長分成兩個半步長計(jì)算, 在前半時(shí)間步長內(nèi),聯(lián)立運(yùn)動方程u分量方程和連續(xù)方程, 在x方向采用隱式差分求解 和u,y方向采用顯式差分求解v; 在后半步長內(nèi)聯(lián)立運(yùn)動方程v分量方程和連續(xù)方程,在y方向采用隱式差分求解 和v,x方向采用顯式差分求解u。這樣反復(fù)運(yùn)用顯隱交替的方法運(yùn)算便能算出每個時(shí)間步長上各點(diǎn)的x和y方向的流速和水深。ADI法為一成熟的算法, 關(guān)于應(yīng)用此方法對方程進(jìn)行離散求解的具體過程見文獻(xiàn)[7], 本文不再贅述。

2 赤潮生態(tài)數(shù)學(xué)模型的建立

赤潮生物動力學(xué)模型的基本方程在二維擴(kuò)散方程基礎(chǔ)上, 參考夏綜萬等[5]的赤潮生物動力學(xué)方程,考慮了速度修正的水平渦動系數(shù)xλ,yλ。

E為海水中營養(yǎng)物質(zhì)的質(zhì)量濃度;N為赤潮藻類的生物量密度;xλ為x方向渦動擴(kuò)散系數(shù);yλ為y方向渦動擴(kuò)散系數(shù); 其他各參數(shù)意義及其取值見表1。

表1 其他各參數(shù)意義及其取值Tab. 1 The significance and value of other parameters

其中,mE為營養(yǎng)物質(zhì)的半飽和常數(shù),是指可作為藻類細(xì)胞能正常生長所需維持水中有效形式營養(yǎng)鹽的臨界濃度, 也可用于比較不同浮游植物吸收營養(yǎng)鹽能力的大小。在光強(qiáng)、水溫及其他條件適宜而營養(yǎng)鹽含量較低時(shí),mE值越小的浮游植物越容易發(fā)展成為優(yōu)勢種;mE值越大的浮游植物會因缺乏營養(yǎng)鹽使生長受到限制。以上的參數(shù)取值, 主要是參考了夏綜萬[5]、王洪禮[8]的取值, 并通過反復(fù)調(diào)試得出的結(jié)果。

迎風(fēng)格式是一種最為常見的差分格式, 采用由下游向上游差分的方法代替微分, 離散方向和赤潮生態(tài)數(shù)學(xué)模型的擴(kuò)散方向一致, 因此本文采用此方法, 對赤潮生物動力學(xué)方程進(jìn)行離散, 具體方法見文獻(xiàn)[7]。

3 二維赤潮生態(tài)數(shù)學(xué)模型在渤海的應(yīng)用

3.1 二維平面水動力數(shù)學(xué)模型在渤海的應(yīng)用

3.1.1 渤海概況[9]

渤海位于 37°07′～41°00′N、117°35′～121°10′E之間。它的東面有渤海海峽與黃海相通, 其余三面均為大陸所圍。在我國行政區(qū)劃上, 它的北界屬于遼寧省, 西界屬河北省和天津市, 南界為山東省, 東以遼東半島的老鐵山西南角與山東半島的蓬萊間的連線——渤海海峽為界。是由遼東灣、渤海灣、萊州灣、中央盆地和渤海海峽五部分組成的。流入渤海的主要河流有遼河、海河、灤河、黃河等。渤海南北長約300 n mile, 東西寬約160 n mile, 沿岸所圍的形狀好似一個葫蘆。渤海海域平均水深約18 m, 大部分海域水深較淺, 深度小于30 m的海域占總面積的95%左右, 海底地形平坦開闊, 5 m與25 m兩級等深線最長, 平均坡度僅0°00′28″。水深大于30 m的海域大部集中在渤海海峽地區(qū), 其中最大深度可達(dá)86 m。渤海的地形等深線圖如圖1所示, 圖中等深線單位為m。

圖1 渤海地形等深線圖Fig. 1 Bathymetric chart of the Bohai terrain

3.1.2 模型網(wǎng)格的剖分

該模型仍采用ADI有限差分方法計(jì)算赤潮發(fā)生時(shí)各點(diǎn)水位及流速, 用正方形網(wǎng)格對計(jì)算區(qū)域進(jìn)行剖分。考慮到模型范圍僅覆蓋了渤海, 因此計(jì)算網(wǎng)格相應(yīng)加密, 采用的時(shí)間步長為60 s,空間步長為3 km, 剖分后為127列×142行。模型網(wǎng)格剖分圖如圖2所示。

3.1.3 流場模擬與流速驗(yàn)證

圖2 渤海網(wǎng)格剖分示意圖Fig. 2 Grid sketch of the Bohai Sea

為了準(zhǔn)確模擬 2004年 6月渤海所發(fā)生的赤潮,作者運(yùn)用該模型對渤海2004年5月30日7時(shí)到6月30日23時(shí)的水流潮汐進(jìn)行了模擬計(jì)算。水邊界條件采用邊界點(diǎn) P(97,60)、Q(84,25)(如圖 2所示)在2004年5月30日7時(shí)到6月30日23時(shí)的潮位過程,其中括號中的數(shù)值表示網(wǎng)格位置坐標(biāo), 潮位過程如圖 3所示。并在渤海取一點(diǎn)進(jìn)行流速的大小及方向的驗(yàn)證。2004年5月30日7時(shí)開始后的第127 小時(shí)和第603 小時(shí)的流場示意圖如圖4所示。渤海測流點(diǎn)的流速和流向驗(yàn)證曲線如圖5所示。

從以上驗(yàn)證曲線可以看出, 模型計(jì)算結(jié)果與實(shí)際基本吻合。盡管用于模型計(jì)算的實(shí)測資料有限, 以及計(jì)算中參數(shù)的概化與近似帶來不可避免的誤差,使得計(jì)算結(jié)果與實(shí)際有所偏差, 但模型計(jì)算所得的潮流場及潮位過程與實(shí)際情況的符合程度都比較令人滿意, 可以作為進(jìn)一步計(jì)算赤潮模型的基礎(chǔ)。

3.2 二維赤潮生態(tài)數(shù)學(xué)模型在渤海的應(yīng)用

2004年 6月 11日, 渤海海域出現(xiàn)棕囊藻赤潮,16日赤潮消失。棕囊藻隸屬定鞭藻綱, 為廣溫廣鹽性藻類, 對環(huán)境有出色的適應(yīng)能力。在世界各地的海洋中, 棕囊藻赤潮一般于春季爆發(fā)。從專業(yè)角度講,1L海水中的棕囊藻密度要達(dá)到 1×107個, 才能稱為赤潮。從2004年5月30日開始, 我們對模型進(jìn)行了1個月的計(jì)算, 得到各時(shí)間段的藻類密度分布圖, 圖6為6月11日、16日的藻類密度分布圖, 模擬了這 次渤海赤潮的生消過程。

圖3 邊界點(diǎn)潮位過程Fig. 3 Tide process of border point

圖4 流場圖Fig. 4 Flow field

圖5 測流點(diǎn)流速、流向驗(yàn)證Fig. 5 Velocity verification of measuring points

圖6 藻類密度分布圖Fig. 6 Distribution map of algae density

3.2.1 初始條件

初始時(shí)刻, 渤海海域內(nèi), 每一點(diǎn)的營養(yǎng)物質(zhì)濃度都取為4× 1 0μg/L, 藻類密度都取為 0.1個/L。即：當(dāng)t=0時(shí), 取Ei,j= 4 × 1 0μg/L ,Ni,j=0.1 個/L。

3.2.2 從衛(wèi)星遙感圖像上提取赤潮信息

赤潮水體與非赤潮水體光譜相比具有明顯的吸收峰和反射峰, 這些特征在衛(wèi)星遙感探測上則會表現(xiàn)出在相應(yīng)波段的離水輻射率(透射入水的輻射光經(jīng)水分子、浮游生物、懸浮物等散射, 其中一部分離開水面反射出來, 這部分稱為離水輻射率。)發(fā)生不同程度的增強(qiáng)或減弱, 而非赤潮水體在相應(yīng)波段的變化較小。因此, 可以利用赤潮水體的這一特征作為赤潮是否發(fā)生的一個條件。

MODIS第三通道(459～479 nm)遙感信息中含有赤潮水體的吸收信息, 第四通道(545～565nm)遙感信息中含有赤潮水體的反射信息, 因此, 利用上述兩通道反射率的比值建立如下算法用于提取海水中赤潮水體的信息。

其中,R3,R4分別是MODIS第三通道和第四通道的反射率,Cr是常數(shù), 其值的大小與發(fā)生赤潮的海區(qū)和赤潮藻種有關(guān)。上式的比值在一定程度上表明海水表層浮游藻類的聚集程度, 也是水體中葉綠素a濃度的一種反映, 它隨著水體藻類細(xì)胞密度的增加而增大。因此, 可以據(jù)此提取海水中的赤潮信息。

3.2.3 模型驗(yàn)證與結(jié)果分析

對EOS/MODIS衛(wèi)星拍攝的2004年6月份的渤海海區(qū)衛(wèi)星遙感圖像(如圖 7)進(jìn)行了處理, 提取出海水中的赤潮信息, 并計(jì)算出赤潮面積, 使其與模型計(jì)算出的赤潮面積進(jìn)行對比(如圖8所示), 通過衛(wèi)星遙感照片算出的赤潮面積大約是2 862 km2, 數(shù)值模擬算出的赤潮面積為2 637 km2, 面積大小基本相符,結(jié)果比較令人滿意。

選擇渤海內(nèi)幾個典型格點(diǎn)A(15, 44), B(46, 12),C(68, 63), 它們分別位于渤海灣、萊州灣及中央盆地(圖2)。將這幾點(diǎn)處的營養(yǎng)物質(zhì)濃度E和藻類密度N的時(shí)間序列繪于圖 9, 可以看出, 赤潮的發(fā)生與海洋中的營養(yǎng)物質(zhì)濃度有著密切地關(guān)系, 即： 營養(yǎng)物質(zhì)濃度大量增加后導(dǎo)致藻類急劇增殖, 而藻類的繁殖吸收了大量的營養(yǎng)物質(zhì)又致使?fàn)I養(yǎng)濃度轉(zhuǎn)為下降,這和前人研究的結(jié)果一致[10～12]。從這幾點(diǎn)的情況來講, 處于渤海灣、萊州灣的這些地方, 似乎比渤海中開闊區(qū)域的部分, 如 C點(diǎn)更容易聚集營養(yǎng)物質(zhì)和浮游植物, 因而更容易形成赤潮。這也說明了在這些地方潮流量級小, 擴(kuò)散能力差, 營養(yǎng)物質(zhì)容易聚集, 一旦產(chǎn)生了赤潮, 也不容易分散。這種趨勢和這一帶經(jīng)常容易產(chǎn)生赤潮的情況相符[10]。

4 結(jié)論

圖7 2004年6月11日MODIS渤海赤潮遙感圖Fig. 7 Remote sensing map of the red tide in the Bohai Sea on June 11, 2004

圖8 赤潮擴(kuò)散范圍對比圖Fig. 8 Comparison of the proliferation of red tide

圖9 各個典型格點(diǎn)E和N的時(shí)空變化序列Fig. 9 Sequential changes of E and N in time and space of points

本文綜合考慮了潮、流、營養(yǎng)物質(zhì)對赤潮過程的影響, 建立了一個水動力學(xué)與生物動力學(xué)相結(jié)合的二維生態(tài)數(shù)學(xué)模型, 并以渤海棕囊藻赤潮為例進(jìn)行了數(shù)值模擬, 對赤潮發(fā)生的數(shù)值模擬問題進(jìn)行了初步嘗試。

1)以不可壓縮流體二維淺水環(huán)流方程為基本控制方程, 采用有限差分 ADI方法對方程進(jìn)行離散,建立了渤海水流潮汐模型, 計(jì)算潮流場, 驗(yàn)證了測流點(diǎn)流速、流向, 取得了較好的模擬效果。

2)采用有限差分思想, 應(yīng)用迎風(fēng)格式對赤潮生物動力學(xué)方程進(jìn)行離散, 建立二維赤潮生態(tài)數(shù)值模型并把該模型應(yīng)用于渤海海區(qū)。

3)用該赤潮生態(tài)數(shù)值模型模擬驗(yàn)證了 2004年 6月11～16日發(fā)生在渤海的棕囊藻赤潮, 并與渤海海區(qū)赤潮衛(wèi)星遙感圖對比, 模擬結(jié)果與實(shí)測資料較吻合。

[1] 齊雨藻, 黃偉健, 邱璇鴻.大鵬灣夜光藻種群動態(tài)的時(shí)間序列模型[J]. 暨南大學(xué)學(xué)報(bào), 1991, 12(3)：96-103.

[2] Chen C W. Concepts and utilities of ecologic model[J].Sanitary Engineering Division, 1970, 96： 1085-1097.[3] Di Toro D M, Oconnor D J, Thomann R V. A dynamic model of phytoplankton population in the Sacramento san Joaquin Delta[J]. Advances in Chemistry, 1971,106： 131-180.

[4] 王壽松, 劉子煌, 夏綜萬, 等. 封閉環(huán)境中赤潮發(fā)生過程的數(shù)學(xué)模擬[J]. 海洋與湖沼, 1998, 29(2)：163-168.

[5] 夏綜萬. 大鵬灣的赤潮生態(tài)仿真模型[J]. 海洋與湖沼, 1997, 28(5)： 468-474.

[6] 田峰, 葛根, 楊晨. 水動力與生態(tài)耦合的赤潮藻類生長模型研究[J]. 海洋技術(shù), 2007, 26(2)： 34-37.

[7] 李萬平. 計(jì)算流體力學(xué)[M]. 武漢： 華中科技大學(xué)出版社, 2004. 160-182.

[8] 王洪禮. 渤海赤潮藻類生態(tài)動力學(xué)模型的非線性動力學(xué)研究[J]. 海洋技術(shù), 2002, 21(3)： 8-12.

[9] 中國科學(xué)院海洋研究所海洋地質(zhì)研究室. 渤海地質(zhì)[M]. 北京： 科學(xué)出版社, 1985. 1-44.

[10] 齊雨藻. 中國沿海赤潮[M]. 北京： 北京出版社, 2003.1-302.

[11] 錢宏林. 南海北部沿海夜光藻赤潮的生態(tài)模式研究[J]. 生態(tài)科學(xué), 1994, 1： 39-46.

[12] 王壽松, 劉子煌, 夏綜萬, 等. 封閉環(huán)境中赤潮發(fā)生過程的數(shù)學(xué)模擬[J]. 海洋與湖沼, 1998, 29(2)： 163-168.

Two-dimensional ecological mathematical model of red tide and its application in the Bohai Sea

LI Da-ming1, LIN Yi1, SONG Shuang-xia1, XIE Yi-yang2, HAN Su-qin2, LI Pei-yan2
(1. Tianjin University School of Civil Engineering & Key Laboratory of Harbor & Ocean Engineering Ministry of Education, Tianjin 300072, China; 2. Institute of Tianjin Meteorological Bureau, Tianjin, 300074, China)

Mar., 9, 2009

two-dimensional ecological mathematical model of red tide; ADI method; upwind scheme; Phaeocystis; the Bohai Sea; remote sensing

The two-dimensional ecological mathematical model of red tide was established by combining water dynamics with biological dynamics. The water dynamics mathematical model was established by dispersing and solving the incompressible two-dimensional circulation equations in shallow water by the ADI (Alternating Direction Implicit) Method of the finite difference method, and the red tide biodynamic equations was dispersed by upwind scheme. The established 2D ecological mathematical model of red tide was applied to calculate the red tide of Phaeocystis, which occurred in the Baohai Sea during the period of June 11～16th, 2004. Satellite remote sensing data of the red tide area and calculation results based on the model were in good agreement. It shows that the model can simulate the process of growth and death of red tide well and provide scientific basis for red tide forecast in the Bohai Sea.

X55

A

1000-3096(2010)09-0087-07

2009-03-09;

2010-07-05

李大鳴(1957-), 男, 河北棗強(qiáng)人, 教授, 博士, 博士生導(dǎo)師, 主要從事水力學(xué)及河流動力學(xué)、海岸工程, 電話： 022-87401579,E-mail： lidaming@tju.edu.cn

(本文編輯： 劉珊珊)