08美國(guó)數(shù)學(xué)建模比賽A題.doc

Team # 3694 Page 30 of 30不可忽視的事實(shí)：一個(gè)海平面上升預(yù)測(cè)模型杜克大學(xué)：Jason Chen，Joonhahn Cho，Brian Choi目 錄目錄.1問(wèn)題介紹.2II模型建立.4 海平面上升模型.4 溫度數(shù)據(jù)分析.5 冰原模型.5 物質(zhì)平衡-積累模型.6物質(zhì)平衡-消融模型.7物質(zhì)平衡與海平面上升模型.9熱膨脹模型.9局域化說(shuō)明.9III. 結(jié)果.11海平面上升模型輸出數(shù)據(jù).11淹沒(méi)模型仿真結(jié)果.12IV. 討論與結(jié)論.17V建議.18參考文獻(xiàn).20附錄A 海平面上升仿真代碼.21附錄B 地質(zhì)圖點(diǎn)矩陣腳本.25附錄C 淹沒(méi)模型仿真代碼.27附錄D 佛羅里達(dá)城市數(shù)據(jù)初值.29I問(wèn)題介紹據(jù)有力證據(jù)表明全球變暖趨勢(shì)已經(jīng)存在，目前已經(jīng)建立了很多有效模型來(lái)預(yù)測(cè)將來(lái)的氣候變化。在近15年內(nèi)，全球氣溫上漲了大約0.5，全球溫度達(dá)到了過(guò)去一千年以來(lái)的最高水平。雖然變暖的趨勢(shì)相當(dāng)明顯，但人們對(duì)這種大規(guī)模的氣候變化仍然是知之甚少。全球變暖帶來(lái)的最令人擔(dān)心的后果之一便是海平面的上升。TOPEX/Poseidon衛(wèi)星測(cè)高儀觀測(cè)出，在1993-1998年期間，海平面每年上漲3.2 0.2毫米。泰特斯等人估計(jì)，僅在美國(guó)，海平面上升1米，就可能會(huì)導(dǎo)致270-475億美元的賠償金。若干復(fù)雜的因素促成了海平面的上升。人們一直認(rèn)為，因溫度變化而導(dǎo)致的水的熱膨脹是海平面上升的重要組成部分；然而，近期研究表明，單是熱膨脹并不能解釋已觀察到的大多數(shù)海平面的上升。大塊冰原的物質(zhì)平衡，特別是格陵蘭冰原的物質(zhì)平衡，現(xiàn)在被認(rèn)為是海平面上升的主要因素。物質(zhì)平衡是由2個(gè)主要過(guò)程控制的，積累（冰塊匯集到冰原）和消融（冰塊脫離冰原）。積累主要是由于降雪，消融主要是由于水分的升華和冰體的融化。與大家普遍的認(rèn)識(shí)恰恰相反的是，浮冰對(duì)海平面的上升并沒(méi)有發(fā)揮重大的作用。根據(jù)阿基米德原理，重量為W（假設(shè)淡水的密度為water）的浮冰融化成密度為ocean 的水體，其體積的增量為V，即：（1） 海水的密度大約是1024.8kg/m3 ；北極海冰的重量大約是2x103 kg。因此，如果北極所有的海冰都融化，體積變化將由下式給出： （2）近360噸水可以導(dǎo)致海平面上升1毫米。（3）海平面如此微小的變化對(duì)我們的模型來(lái)說(shuō)是無(wú)關(guān)緊要的，因?yàn)榫纫呀?jīng)低于千分之一毫米。我們也忽略了南極冰原的影響，因?yàn)槠淇傮w影響海平面上升的量是極少的，難以量化。在1978到1987年之間，衛(wèi)星微波輻射的數(shù)據(jù)表明，北極的冰體減少了3.5%，然而同時(shí)，南極的冰體卻沒(méi)有可統(tǒng)計(jì)的顯著變化?？ㄍ吡欣热艘呀?jīng)預(yù)測(cè)出未來(lái)50年內(nèi)，南極冰體融化的最低限度。出于這樣的考慮，我們?cè)谀Ｐ椭袃H考慮格陵蘭冰原的影響。關(guān)于物質(zhì)平衡和熱膨脹理論，已經(jīng)出現(xiàn)了一些模型。但是這些模型考慮了很多方面的變量，非常復(fù)雜，并且經(jīng)常彼此不統(tǒng)一。我們希望提出一個(gè)基于簡(jiǎn)單的物理過(guò)程的模型，僅僅體現(xiàn)溫度和時(shí)間的影響。利用這種方法，分析影響氣候變暖的過(guò)程將被簡(jiǎn)化，并且使海平面上升對(duì)溫度的依賴體現(xiàn)得更為明顯。此外，我們建立了一個(gè)可以擴(kuò)展到在若干不同溫度影響下的計(jì)算的模型，允許我們直接比較碳的排放量對(duì)海平面上升的影響。模型概述深入了解冰原融化能夠給海平面上升提供寶貴的見(jiàn)解。通過(guò)構(gòu)造一個(gè)可以整合冰原融化和熱膨脹作用的框架，我們可以估計(jì)50年期限的全球平均海平面。這個(gè)模型實(shí)現(xiàn)了幾個(gè)重要的目標(biāo)：1. 正確的吻合了過(guò)去海平面上升的數(shù)據(jù)2. 為預(yù)測(cè)50年期限的海平面提供了充分的一般性3. 作為一個(gè)獨(dú)立的全球溫度和時(shí)間函數(shù)計(jì)算佛羅里達(dá)州海平面的上升最后，這一模型從邏輯上預(yù)測(cè)了人類的人口。我們尤其分析了佛羅里達(dá)州海平面上升的影響，很多人認(rèn)為由于它的海拔低而且臨近大西洋而使它極易受到攻擊。從分析中我們?cè)u(píng)估了將全球變暖導(dǎo)致海平面上升帶來(lái)的損失降低到最小的可能策略。假定為精簡(jiǎn)模型進(jìn)行一些主要的假設(shè)：1. 海平面上升主要由于兩個(gè)因素：格陵蘭冰原的積累/消融平衡和大洋的熱膨脹。忽略了例如冰裂和直接的人類干預(yù)等作用的影響，而這些影響很難正確的模擬并且對(duì)海平面上升影響很小。2. 空氣是融化冰原的唯一熱源。格陵蘭大陸是永久凍結(jié)帶，而且由于它表面覆蓋的大量冰原，所以假設(shè)它的溫度相對(duì)穩(wěn)定。這方便我們使用對(duì)流作為熱傳遞的模型。3. 冰原內(nèi)部穩(wěn)態(tài)的溫度線性變化。這一假設(shè)可以使我們用紐曼條件解決熱量方程式。通過(guò)減少熱量方程式的穩(wěn)態(tài)條件，我們可以解決齊次邊界條件。4. 升華和融化過(guò)程互不干涉。這一假設(shè)徹底簡(jiǎn)化模型中分別考慮升華和融化帶來(lái)的計(jì)算。假設(shè)又是合理的，升華發(fā)生在零度以下，而融化在零度以下根本不會(huì)發(fā)生。因此在我們的模型中兩個(gè)過(guò)程是時(shí)間上獨(dú)立的。5. 冰原的表面在溫度、壓力和化學(xué)組成上是均勻的，因?yàn)樵谖覀兊目蚣苤懈窳晏m大陸高分辨率空間溫度數(shù)據(jù)是不能得到的，所以這一假設(shè)是必要的。另外我們?nèi)鄙儆?jì)算的資源和時(shí)間來(lái)仿真如此一個(gè)要求用有限元法和網(wǎng)格劃分的復(fù)雜拓?fù)鋵W(xué)的變量。確定問(wèn)題M表示格陵蘭冰原的物質(zhì)平衡，我們?cè)跀?shù)量上估計(jì)海平面上升結(jié)果出現(xiàn)的W衛(wèi)星激光測(cè)距，給出一個(gè)以溫度為變量的函數(shù)，這些上升是M和熱膨脹TE影響的總和根據(jù)當(dāng)?shù)刳厔?shì)做的校正。進(jìn)一步，我們必須定量和定性地分析全球變暖對(duì)佛羅里達(dá)州主要城市和區(qū)域長(zhǎng)期（50年）的影響，結(jié)果是高SLR。這一分析可以用來(lái)為如何最好的應(yīng)對(duì)和減少SLP效應(yīng)提出建議。II. 模型建立建立海平面上升模型海平面的上升大多是因?yàn)楦窳晏m冰原的物質(zhì)平衡及由于氣候變暖帶來(lái)的熱膨脹。為了表示海平面的上升以及其他后處理的影響，我們用了一個(gè)物質(zhì)平衡模型和熱膨脹模型。邏輯模擬過(guò)程詳細(xì)見(jiàn)圖1 。模擬開(kāi)始IPCC監(jiān)測(cè)大氣中的碳含量模型熱膨脹格陵蘭冰原物質(zhì)平衡升華積累融化區(qū)域增長(zhǎng)趨勢(shì)海平面上升淹沒(méi)模擬提出建議EdGCM臨時(shí)程序調(diào)用圖1：模擬流程圖溫度數(shù)據(jù)分析溫度數(shù)據(jù)在我們的模型中是唯一使用的數(shù)據(jù)，所以必須仔細(xì)分析考慮。因?yàn)樾枰⒉煌榫跋碌哪Ｐ停覀兊臏囟葦?shù)據(jù)就必須包含這些受控制的且只有一個(gè)不同變量的情景。另外，溫度數(shù)據(jù)必須是高質(zhì)量的，并且能保證為我們的模擬提供正確的臨時(shí)分析。為了到達(dá)這個(gè)目的，我們決定使用一個(gè)全球氣候模型（GCM），通過(guò)輸入便于我們控制的變量來(lái)處理我們自己的溫度數(shù)據(jù)。由于計(jì)算能力和時(shí)間的限制，我們選擇EdGCM模型。EdGCM模型是一個(gè)為教學(xué)設(shè)計(jì)的快速模型。該模型是基于美國(guó)航天局GISS氣候變化模型設(shè)計(jì)出來(lái)的。此程序可以滿足我們所有的需要，特別是在模擬的速度方面（大約10小時(shí)模擬一個(gè)50年的氣候階段），使我們可以分析若干個(gè)不同的溫度情景。溫度情景分析中納入了IPCC第三次評(píng)估報(bào)告得到的關(guān)于碳排放量的三個(gè)因素在IS92系列中的結(jié)果，低、高和中等。IS92e(高)，IS92a(中)以及IS92c(低)情景都是在EdGCM中非常近似的結(jié)果。這些近似的碳作用在下圖2中表示出來(lái)。其他所有的參量根據(jù)美國(guó)航天局GISS模型，被默認(rèn)保存。用這種方式，在三個(gè)時(shí)間序列中，獲得了全球地表空氣溫度。圖2：二氧化碳在EdGCM模型中的影響對(duì)于EdGCM模型，一個(gè)不利的影響是它只能輸出全球氣溫變化。雖然能計(jì)算區(qū)域氣溫變化，但是卻很難到達(dá)并且所得結(jié)果空間精度較低。然而，根據(jù)Chylek等人的研究，格陵蘭島的氣溫和全球氣溫之間的關(guān)系可以近似的表示為：（4）Chylek等人用未受NAO影響的區(qū)域證明，并用氣候模型的輸出預(yù)測(cè)了這一結(jié)果。冰原模型將冰原簡(jiǎn)化假設(shè)為一個(gè)矩形，且冰原上層的每一處溫度假設(shè)為恒定的Ta。這是因?yàn)槲覀兊臍夂蚰Ｐ蜎](méi)有達(dá)到格陵蘭島的空間精度，所以微小的溫度變化被忽略，不考慮。冰原下層的凍土層有恒定的溫度Tl。冰原模型的描述如圖3所示。圖3：冰原模型的側(cè)面圖為了計(jì)算冰原融化和升華的熱通量，我們假設(shè)其為無(wú)限多的微分單元，如下圖4。圖4：冰原微分單元最初，高度h的計(jì)算是利用Williams等人提供的數(shù)據(jù)。 我們的模型中主要的計(jì)算方法是依據(jù)物質(zhì)平衡，物質(zhì)平衡即計(jì)算減去消融量后所積累的量。積累，是指冰原上增加冰塊，其主要形式是降雪。冰體的消融主要是兩個(gè)過(guò)程的結(jié)果，即升華和融化。物質(zhì)平衡積累首先，我們建立積累的模型。Huybrechts等人證明格陵蘭島的溫度還沒(méi)有高達(dá)可以融化大面積的冰雪。此外，以Knight的經(jīng)驗(yàn)推斷，冰原的積累率和時(shí)間存在非常近似的線性關(guān)系，覆蓋格陵蘭大陸的冰雪積累量是0.3m/year。因此，積累率是0.025m/month。根據(jù)物質(zhì)平衡得出： （5）LD的乘積指冰原的表面積。物質(zhì)平衡消融其次，我們對(duì)消融的升華和融化這兩個(gè)過(guò)程進(jìn)行建模。升華率（質(zhì)量變化量）由下式給出： （6）Mw 為水的分子質(zhì)量。這個(gè)表達(dá)式可以從理想氣體法則和麥克斯韋玻爾茲曼分布導(dǎo)出。帶入esat 的巴克表達(dá)式，我們可以得出： （7）巴克方程適用于范圍較大的溫度和壓力的變化，適合格陵蘭島的環(huán)境。雖然近似地舍去了極端溫度和壓力，但是這樣一來(lái)計(jì)算相對(duì)簡(jiǎn)單。為了將質(zhì)量變化量轉(zhuǎn)換成冰塊的厚度變化量，我們根據(jù)冰塊的密度來(lái)劃分質(zhì)量變化量。由此我們可以表述冰塊厚度的變化，如下所示： （8）d表示焊敷系數(shù)(克/小時(shí)安)，即d=(1-deposition rate)=0.01。此條件是必要的，因?yàn)樯A和沉積是保持平衡的。利用此升華率表達(dá)式，現(xiàn)在可以找出冰原每一時(shí)刻在計(jì)算模型中的微小的厚度變化。由消融中的升華過(guò)程引起的冰原厚度的變化可以由下式給出： （9）公式中的h是冰原的厚度，而t是一個(gè)時(shí)間段后所經(jīng)歷的時(shí)間。我們用已經(jīng)導(dǎo)出的表達(dá)式代替并且替代了關(guān)于水分子重力場(chǎng)的已知數(shù)值 （10）這個(gè)方程規(guī)定了冰的升華。對(duì)于消融的第二個(gè)組成部分，融化模型，我們應(yīng)用熱力方程。熱力方程規(guī)定了如下關(guān)系： （11）其中k=0.0104是冰的熱擴(kuò)散率。為了解決在諾伊曼條件下的熱力方程，我們以與相同的邊界條件并且獨(dú)立于時(shí)間的條件假設(shè)一個(gè)常量。剩余溫度V有齊次邊界條件以及在-中發(fā)現(xiàn)的起始提條件。因此我們可以重新寫(xiě)出如下熱力方程： （12）熱力方程的恒定量解決被如下給出: (13)其中滿足0xS(t)與0tk，在t不變時(shí)U(x,t)0，因此當(dāng)kx Sx_val = S;endif x_val Sx_val = S;endif x_val Sx_val = S;endif x_val 0error(X is less than 0);endtimecourseL = timecourseL x_val;h3 = x_val + accu_rate;if mcount = 13ann_tempchg = ann_tempchg + Ta; %temp change recordingmcount = 0;SLR3 = SLR3 + (6.89*ann_tempchg); %thermal expansion effect in centimetersycount = ycount + 1;SLR3 = SLR3 + (inih-h3)*700/h; %melting effect in centimetersSLRtot3(ycount,1) = SLR3 - (unifrnd(flomin,flomax)*(ycount) +(floavg*ycount); %record the rise(local) ; %record the rise ann_tempchg = 0;endendfunction out = U(Tl, Ta, t, S, x)%heat equationn = 1:100;Us = Tl+1e-2*(Ta-Tl)*x/S;V = (Tl+Ta)/2+sum(2.*(Ta-Tl).*(-1).n-1).*exp(-n.2.*pi2.*t./S.2).*cos(n.*pi.*x./S)/(n.*pi).2);out = Us+V;附錄B 建立地質(zhì)圖點(diǎn)矩陣腳本% Create elevation data from GTOPO30 datadatagrid refvec = gtopo30(W100N40, 2, 24 31, -90 -78); %gather topological data for 1 resolutiondatagrid(isnan(datagrid) = -100;datagridt = flipud(datagrid);image(datagridt);% % Create an indexing matrix to quickly determine latitude and longitude for% % lat = 31:(-1/60):24;% % lon = -90:(1/60):-78;% % latrow = 0;% % for latind = 1:length(lat)-1% % for lonind = 1:length(lon)-1% % latlonmat(latind, lonind).coord = lat(latind) lon(lonind)