李大鳴, 林 毅, 宋雙霞, 解以揚, 韓素琴, 李培彥

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

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

李大鳴1, 林 毅1, 宋雙霞1, 解以揚2, 韓素琴2, 李培彥2

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

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

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

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

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

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

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

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

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

1.1.1 基本控制方程

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

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

1.1.2 定解條件

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

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

初始條件： 當t=0時,ξ=0,u=v=0

1.2 ADI的有限差分方法

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

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

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

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

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

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

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

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

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

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ū), 其中最大深度可達86 m。渤海的地形等深線圖如圖1所示, 圖中等深線單位為m。

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

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

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

3.1.3 流場模擬與流速驗證

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

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

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

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

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

圖3 邊界點潮位過程Fig. 3 Tide process of border point

圖4 流場圖Fig. 4 Flow field

圖5 測流點流速、流向驗證Fig. 5 Velocity verification of measuring points

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

3.2.1 初始條件

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

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

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

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

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

3.2.3 模型驗證與結果分析

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

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

4 結論

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

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

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

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

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

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

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

[1] 齊雨藻, 黃偉健, 邱璇鴻.大鵬灣夜光藻種群動態(tài)的時間序列模型[J]. 暨南大學學報, 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ù)學模擬[J]. 海洋與湖沼, 1998, 29(2)：163-168.

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

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

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

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

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

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

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

[12] 王壽松, 劉子煌, 夏綜萬, 等. 封閉環(huán)境中赤潮發(fā)生過程的數(shù)學模擬[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-), 男, 河北棗強人, 教授, 博士, 博士生導師, 主要從事水力學及河流動力學、海岸工程, 電話： 022-87401579,E-mail： lidaming@tju.edu.cn

(本文編輯： 劉珊珊)