藏東南波堆藏布江流域古鄉冰期冰川重建

覃長雄，許劉兵，謝金明，周尚哲

（1.華南師范大學 地理科學學院，廣東廣州 510631；2.蘭州大學 資源環境學院，甘肅蘭州 730000）

0 引言

冰川是寒冷氣候的產物，對氣溫和/或降水變化十分敏感［1］，古冰川相對現代冰川的物質平衡線高度（equilibrium-line altitude，ELA）、面積、長度、以及體積等特征參數的變化，對了解古氣候和預測未來氣候變化至關重要。其中，現代冰川的各項特征參數可根據實地觀測、測量、遙感影像解譯等方法直接估算得到，對古冰川而言，通常是根據現代冰川與當前氣候環境要素（如氣溫、降水、地形等）的統計關系，再結合冰川地貌和年代學去反演過去不同時期冰川的各項關鍵參數。從目前的研究現狀來看，對古冰川發育環境及特征參數的重建研究主要集中在末次冰期最盛期（Last Glacial Maximum，LGM）以來，主要原因是LGM以來的冰川地貌保存相對完好，地質年代也能準確的測定。

“古鄉冰期”是20世紀80年代由李吉均等［2］提出來的冰期概念，同時提出的還有“白玉冰期”，這兩個冰期名稱在青藏高原第四紀冰川作用眾多地方性冰期名稱中非常具有代表性，數十年來一直被廣泛沿用。其中，古鄉冰期指代的是傳統意義上的“倒數第二次冰期”，已有的宇宙成因核素10Be地表暴露測年結果［3］顯示其發生于海洋氧同位素階段（marine oxygen isotope stage，MIS）6。雖然古鄉冰期概念已被眾人接受，且發生年代也已界定，但古鄉冰期命名地即藏東南波堆藏布江流域的冰川作用規模（如范圍、冰厚及冰量等）和ELA等關鍵參數仍有待進一步考證。本文在多次詳盡野外調查的基礎上，結合高清遙感影像/衛星照片及定量測年結果，對波堆藏布江流域古鄉冰期的冰川作用范圍進行了重建，并進一步嘗試運用冰面縱剖面（glacier surface profile，GSP）模型、積累區面積比率（accumulation area ratio，AAR）法和末端-源頭高度比率（toe-toheadwall altitude ratio，THAR）法等方法對波堆藏布江流域“古鄉冰期”的古冰川特征參數進行重建。

1 古鄉冰期冰川遺跡及其年代

波堆藏布江發源于念青唐古拉山東段南坡（29°89′～30°67′N、95°08′～95°96′E，圖1），流域面積約1.46×104km2，鄰近雅魯藏布江大拐彎，處于南亞夏季風進入高原的水汽通道上。流域內的高山（最高峰則普海拔6 364 m）迎風面降水豐富，使得該區域發育了規模龐大的第四紀和現代冰川，冰川ELA和冰川末端高度基本處于同緯度的最低值［2］。

圖1 波堆藏布江流域位置及流域內現代冰川［4］和古鄉冰期冰磧壟分布圖Fig.1 Map showing location of the Bodui Zangbo River catchment，modern glaciers［4］，and moraines of the Guxiang Glaciation within the catchment

古鄉冰期指的是傳統意義上的“倒數第二次冰期”，由李吉均等［2］于20世紀80年代提出，其命名地位于藏東南波堆藏布江與帕隆藏布江交匯處下游約17 km的古鄉鎮，即古鄉冰期時，流域內支谷冰川匯入主谷（即波堆藏布江谷地），形成復合型山谷冰川，冰舌到達了古鄉鎮附近，遺存有較大規模的終磧-側磧壟。

除古鄉鎮一帶的冰磧壟外，波堆藏布江主谷兩側亦保存有古鄉冰期的沉積物（圖1）。保存在兩側谷坡、形態完好且規模宏大的側磧壟從上游育仁一直延伸至白玉村附近，拔河高度從550 m降至260 m左右［圖2（a）和（b）］；自白玉向下游保存不甚理想，只見部分谷肩上殘存有零星的冰磧物；但在波堆藏布江匯入帕隆藏布江處的卡達橋附近，尚完好地保留著拔河約200 m、長約3.2 km的高大側磧壟，其上分布著大量的花崗巖漂礫［圖2（c）］。周尚哲等［3］對卡達橋附近古鄉冰期冰磧壟進行了10Be暴露測年，測年結果（約130 ka）顯示古鄉冰期發生于MIS 6。除波堆藏布江主谷外，古鄉冰期冰磧壟在支谷亞龍藏布江兩側谷坡上也有保存，其中，亞龍藏布匯入波堆藏布處的側磧壟保存得最為完整［圖2（d）］。亞龍藏布江北岸的側磧壟壟脊非常典型，從谷口往上游延伸約4.5 km，拔河高度由200 m升至320 m，而谷口南岸的側磧壟因受到支溝的侵蝕，壟脊較為平坦。冰磧壟在上游也有斷續保存，直到海拔3 700 m處，此處拔河約500 m。可見，古鄉冰期時，波堆藏布江流域內支谷冰川都匯入了主谷，形成了主干長達100 km左右的復合型山谷冰川，冰川規?？赡転楝F代的數倍［2］。

圖2 育仁（a）、林瓊（b）、卡達橋（c）和亞龍藏布谷地（d）的古鄉冰期冰磧壟照片Fig.2 Photos of moraines deposited in the Guxiang Glaciation（Moraines in a，b，c，and d are located near Yuren，Linqiong，Kada Bridge，and Yalong Zangbo Valley，respectively）

2 研究方法

2.1 古鄉冰期時冰川范圍的劃定

在多次詳盡野外調查的基礎上，結合已有冰川地貌年代學、高清遙感/衛星影像以及數字高程模型（digital elevation model，DEM），對古鄉冰期冰川作用范圍進行了重建（圖3）。在劃定古鄉冰期冰川范圍時，存在冰磧地貌的地方比較容易確定冰川作用范圍，例如沿著波堆藏布江主谷從則普冰川至古鄉鎮這一段，以及亞龍藏布江與波堆藏布江交匯口附近，可以按照谷地兩側古鄉冰期側磧壟確定古冰川的作用范圍（圖1）。在古鄉冰期冰磧地貌保留不理想甚至完全被侵蝕破壞的地方，如流域內的西南部、西北部以及中部，結合現代冰川的分布范圍來大致劃定古鄉冰期冰川作用范圍。而在流域內的東北部，既無明確的古鄉冰期冰磧地貌，也沒有現代冰川分布，主要依據形態清晰的冰斗和槽谷來劃定古冰川作用范圍。波堆藏布江流域內仍然存在還未定年或無法定年的冰磧壟，要準確地確定古鄉冰期時的冰川作用范圍，仍需進一步實地考察冰川沉積地貌以及進行相關的定年研究。因此，本文采取上述策略來界定古鄉冰期冰川作用范圍，可能與實際情況存在出入，但依據劃定的古冰川范圍估算出的古鄉冰期ELA結果來看，與Zhou等［5］的研究結果基本一致（見3.2節）。

圖3 古鄉冰期時波堆藏布江流域古冰川范圍重建（彩色實線代表35條中流線分布，紅色虛線代表流域邊界）Fig.3 Reconstruction of glacier outline in Guxiang Glaciation in the Bodui Zangbo River catchment（The solid colored lines represent glacier centerlines，and the dotted red line represents the boundary of the Bodui Zangbo River catchment）

2.2 古冰川模擬原理

將根據上文方法得到的古鄉冰期冰川作用范圍導入ArcGIS中，與分辨率為30 m的DEM數據（ALOS World 3D-30 m，https：//www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm）疊加，可得到古冰川的面積。

后運用冰面縱剖面模型模擬古鄉冰期冰川的冰面形態和冰厚等特征參數。冰面縱剖面模型由Schilling等［6］提出，后經Benn等［7］改進，現已被廣泛運用于各地（含青藏高原及周邊山地）古冰川形態的重建［8-12］。

此模型基于的假設為冰川是理想的塑體，其形變由冰的重量及冰表面梯度引起的驅動應力（driving stress，τd）和冰川底部的屈服應力（yield stress，τy）的變化引起。當驅動應力大于屈服應力時，冰川發生運動；當驅動應力小于屈服應力時，就不會發生冰川的運動，但冰川會變厚從而增加驅動應力。為了保持冰川的穩定狀態，即τy=τd，冰川的厚度會不斷進行調整，從而引起驅動應力的變化?；诖说玫搅讼铝斜磉_式。

式中：ρ為冰川密度（900 kg·m-3）；g為重力加速度（9.81 m·s-2）；H為冰川厚度；h為冰表面高程；x為水平坐標（以冰川末端為原點，冰底和冰面高程為縱軸，距冰川末端的長度為橫軸的坐標系）。有時用sinα代替，α為冰川表面坡度［13］。根據以上等式可求出沿著中流線的冰川厚度H，結果為拋物線。

但該式不適用于冰川坡度較大的地區。對于不規則的冰下地形，Benn等［7］指出必須沿著中流線按連續的點（編號1，2，…，i，i+1，i+2）計算，冰的表面坡度表示為

式中：Δx為每個點的間隔距離，本文取Δx=200 m。

因此，由式（1）和式（3）得

式中：H=h-B，B為冰川底部高程，此式可以解出冰面高程h。

為了解決冰舌處冰厚為0的問題，van der Veen［14］進一步推導得到

但該模型是以冰蓋出發進行推導的，冰蓋不受山谷兩側阻力的約束，其流動需要的驅動應力等于冰川底部剪切應力。而山谷冰川則需要考慮山谷兩側阻力的影響，因此引入形狀因子f。

式中：τB為底部剪切應力。

在冰川中流線處，τB=τd，f為

式中：A為冰川橫截面面積；H為冰厚度；p為橫截面的“冰川周長”［15］；f的取值范圍一般為0.4～1［8］。此時的驅動應力大于同等條件下冰蓋流動所需要的驅動應力，在計算時，需要把式（4）中的τy用τy/f替換。

2.3 古冰川模擬過程

使用ArcGIS軟件并結合DEM數據以及研究區的谷歌地球（Google Earth）影像圖，在古冰川范圍上提取出中流線，為了使中流線上插值點分布范圍更全面，減小插值誤差，人工選取了23條分布在面積較大的現代冰川上的中流線和12條分布在主谷位置的中流線（合計35條冰川中流線，圖3）。在Arc-GIS中，按等距離（Δx=200 m）對中流線進行分割，把分割得到的線段轉化為點，最后提取每個點的高程數據，古冰川模擬流程如圖4。

圖4 古冰川模擬流程圖Fig.4 Flow chart of paleo-glacier simulation

在此ExcelTM表格程序中進行古冰川表面高程模擬時，步長（Δx=200 m）、冰的密度（ρ=900 kg·m-3）、重力加速度（g=9.81 m·s-2）是已知的，而B為冰底地形高程，存在現代冰川的地方需要減掉其厚度。輸入的剪切應力τy一般取值100 kPa，根據谷歌地球影像圖，可求得終磧-側磧壟的高程數據作為“目標高程”，將其與中流線上相應點的高程對應起來。下面將以中流線E4（圖3，古鄉至關星冰川北部頂端）為例子，模擬結果見圖5。

模擬過程中，從下游（古鄉）至上游（關星）沿著中流線作出幾道橫剖面，計算每一段的形狀因子［圖5（d）］并插值，之后通過不斷調整中流線上各點的剪切力，使冰面高程與“目標高程”（終磧-側磧壟高度）接近一致［圖5（c）］。

圖5 根據冰川中流線E4模擬的冰川表面高程和其他參數Fig.5 Characteristic parameters of the glacier reconstructed from the glacier centerline E4（Figure a，b，c，and d represents basal shear stress，ice thickness，surface elevation，and shape factor，respectively）

使用上述模擬過程，對其他34條中流線進行模擬，將模擬結果分別用反距離加權（inverse distance weight，IDW）、泛克里金（universal Kriging，UK）和普通克里金（ordinary Kriging，OK）插值法進行插值，最終得到古冰面高程，而古冰川厚度可由模擬古冰面高程減去古冰川冰底高程，結果見圖6。

2.4 冰川ELA估算方法

由Meier［17］提出的AAR法被廣泛用于計算古冰川和現代冰川的ELA，這一方法的關鍵是選擇合適的AAR值。研究表明，分布于不同緯度的冰川，有不同的AAR值。中高緯度冰川的AAR值介于0.5～0.8［18］或0.55～0.65之間［19］，而低緯地區則可能高達0.8［20］。AAR法假定冰川物質平衡隨高度的變化趨近于一條直線，并且考慮了冰川的面積-高程信息。Kern等［21］認為冰川面積S與AAR值之間正相關，且呈對數變化關系：AAR=0.0648×lnS+0.483。對于一些數據缺乏的地區，當S為0.1～1 km2時，AAR=0.44±0.07；當S為1～4 km2時，AAR=0.54±0.07；當S大于4 km2時，AAR=0.64±0.04。Kaser等［22］認為熱帶地區無積雪覆蓋的冰川最佳AAR范圍在0.65～0.70之間，且這個范圍對于熱帶有積雪覆蓋的冰川同樣適用，地貌證據也表明了這個取值范圍同樣適用于古冰川。

本文以朱西冰川、白玉冰川、則普冰川、日母冰川和關星冰川為例，在已有冰川編目數據的基礎上，計算這5條冰川AAR的平均值，并將得到的AAR平均值用來估算古鄉冰期時冰川的ELA。具體操作如下：①根據第一次冰川編目［23］這5條冰川的積累區面積和總面積計算AAR值，取此AAR的平均值來估算古鄉冰期時的冰川ELA；②根據得到的AAR值和第二次編目的冰川范圍來計算各冰川的ELA（表1），其平均值用于下文中選取THAR值?？紤]到這些冰川上部寬、下部窄，粒雪盆（積累區）與冰舌（消融區）都發育正常［24-25］，因此用這5條冰川AAR的平均值（0.73）來恢復古鄉冰期時的冰川ELA。

THAR法也常用于估算古冰川和現代冰川的ELA，THAR值為ELA與冰川最低點高程之差和冰川最高點與最低點高程之差的比值。但該方法并未涉及到冰川的物質平衡或者冰川的面積-高程信息，THAR值的選取有一定的主觀性［26-27］。Charlesworth［28］、Manley［29］曾使用冰川的中值高度來代替ELA。THAR值小于0.5時的ELA小于中值高度［30-33］。Meierding［32］在研究科羅拉多冰川時，認為在0.35～0.40之間的THAR值比較合適。使用該方法估算古冰川ELA時，古冰川的末端高程主要根據冰川終磧的高度計算，即采用沿終磧頂端的最低高度作為冰舌的最低高度［24］；但最高高程不容易確定。為了選取合適的THAR值估算現代冰川和古冰川ELA，本文以第二次冰川編目［4］時觀測的5條現代冰川（分別為朱西冰川、白玉冰川、則普冰川、日母冰川和關星冰川）為例（圖1），將每條冰川上緣處的最高點作為最大高程，使用不同THAR值估算其ELA，估算結果見表2。結果顯示，相同的THAR值下這5條冰川的ELA有著類似的變化趨勢，即隨著冰川海拔的增高其ELA增大，本文取這5條冰川ELA的平均值來代表整個流域的ELA。繼而與第二次冰川編目下這5條冰川ELA（表1）的平均值4 451 m比較，發現THAR=0.35時二者比較接近（表2），因此，在重建古鄉冰期ELA時，THAR值取0.35。

表1 第一次冰川編目［23］下波堆藏布江流域內5條現代冰川的AAR值及估算的第二次冰川編目［4］下各ELA值Table 1 The five modern glaciers AARvalues in the Bodui Zangbo River catchment in glacier inventory of China［23］and their estimated ELA values in second glacier inventory of China［4］

表2 波堆藏布江流域不同THAR值下的5條現代冰川ELA估算結果Table 2 ELA estimates of five modern glaciers under different THAR values in the Bodui Zangbo River catchment

3 結果與討論

3.1 冰川規模變化

從中流線E4上的剪切力分布［圖5（a）］發現，下游（距冰川末端的長度為0～45 km）剪切力變化較小，中游（距離冰川末端的長度為45～75 km）剪切力變化較大，上游（距離冰川末端的長度為75～120 km）剪切力較小，總體分布在7～120 kPa之間。經過計算，剪切力τy平均值為61.3 kPa，平均形狀因子為0.75，等效平均剪切力τy/f為81.7 kPa，在剪切力合理取值范圍內（50～150 kPa）［13］。圖5（b）是冰面高程與“目標高程”相吻合時的冰川厚度分布，由圖5（c）中冰面高程與冰底地形高程相減所得?？梢詮膱D5（b）和5（c）看出，從下游到上游，冰川厚度總體上呈增大的趨勢，其最大值分布在上游，但臨近冰川陡峭的后壁時厚度迅速減小。另外，下游地形較為平緩而中上游地形起伏較大。一般來說，地形平緩的下游更容易積累冰川，但在下游主谷內冰川厚度較小，從白玉村往上游［圖5（c），距離末端50 km以上］冰川厚度逐漸增大，則普村之后冰川厚度基本維持著很大的值，這種趨勢在圖6（d）中也很好地體現了出來。從圖6（a）～（c）中可以得知當時的ELA大致分布在則普村（圖1則普冰川末端）附近，因此冰厚的這種分布是合理的。

圖6 古冰面高程（a，b，c）和古冰川厚度（d）模擬結果以及不同古冰面高程下的ELA重建結果（a，b，c中的黑線）（a，b，c分別為根據IDW、UK和OK插值法得到的古冰面高程）Fig.6 The modeled glacier surfaceelevation（a，b，c）and glacier thickness（d）and reconstructed ELAs（black linesin a，b，c）based on different glacier surface elevation（Figures a，b，c show the glacier surface elevation using IDW，UK，and OK interpolation methods，respectively）

通過對其余34條中流線的模擬，可以重建古鄉冰期時古冰川的表面高程和冰川厚度（圖6），大致估算出古鄉冰期時波堆藏布江谷流域內冰川的形態參數，即流域內平均冰川厚度約360 m，最大厚度約893 m，面積約2 648 km2，體積約953 km3。古鄉冰期時，流域內約63%的區域被冰川所覆蓋，而根據第二次冰川編目統計［4］的數據，現代冰川的覆蓋率約為22.4%。

3.2 古冰川ELA

冰面縱剖面模型只能模擬中流線上的冰面高程和冰川厚度，因此重建古冰川整個冰面高程時需要進行插值。本文使用了IDW、UK和OK等三種常用空間插值方法，得到了三個不同的古冰面高程分布（圖6），結合AAR工具［12］和THAR法計算了相應的冰川ELA，結果顯示：當采用AAR法時，古鄉冰期時冰川的ELA分別為4 021 m、3 867 m、3 847 m。采用THAR法時，估算出古鄉冰期時冰川ELA分別為3 720 m、3 649 m、3 649 m（表3）。從估算的結果來看，AAR法與THAR法得到的古冰川ELA相差較大，可能有以下幾個方面的原因：①使用AAR法計算ELA時需要考慮冰川形狀、類型和規模等特征［26］，而THAR法沒有考慮冰川形狀或不同谷地的面積-高程變化［34］；②采用AAR法時，可能會因為重建的古冰川表面面積不準確而產生誤差；③使用THAR法時，若雪崩對冰川的積累影響較大，則會難以界定冰川的最高高程。因此，本文取兩種估算方法得到的古冰川ELA平均值作為最終結果，這樣能有效避免使用單一計算方法所帶來的誤差。

表3 波堆藏布江流域冰川ELATable 3 ELAs in the Bodui Zangbo River catchment

不同的空間插值方法得到的古冰面高程對于重建古冰川冰量、ELA等重要參數影響較大。為選擇合適空間插值方法，本文在中流線上的數千個插值點中隨機抽取10個點，對比其不同插值結果與模擬值（表4中的實測高程）之間的誤差，結果見表4。從誤差指標可以看出，IDW插值的古冰面高程值的絕對誤差小于5 m，與實測值之間的相對誤差在0.01%～0.1%之間。相對地，UK和OK插值的絕對誤差分別可達347 m和277 m，相對誤差分別在0.002%～12%和0.002%～10%之間。由此可見，IDW插值法的古冰面高程插值精度明顯好于UK和OK插值法。主要原因是IDW插值法對已知高程點設置權重值，并搜索插值點與已知點之間的最短距離來進行空間插值，而UK和OK插值法要求已知高程點的密度分布較大，對于密度較小的區域，空間插值精度較低。本文選取中流線上每個點之間的間隔為200 m，密度較小，已知點較多，因此UK和OK插值法的精度要低于IDW插值法。

表4 不同空間插值方法得到的古冰面高程精度對比Table 4 Precision comparison of paleo-glacier surface elevation estimated from different spatial interpolation methods

最終選擇IDW插值得到的古冰面高程，結合AAR法和THAR法計算古鄉冰期時波堆藏布江流域內冰川ELA，結果取平均值3 871 m（圖6和表3），相較于現代冰川ELA（4 455 m）下降了約584 m（表3）。Zhou等［5］曾對古鄉冰期時波堆藏布江流域內冰川ELA進行估算，結果顯示，古鄉冰期時冰川ELA約為3 800 m，相較于現代冰川ELA（4 600 m）下降約800 m，本文古鄉冰期時冰川ELA的平均值與此前研究之間僅相差約71 m，下降值相差約216 m。

4 結論

本文在冰川地貌外野實地考察和已有年代學為基礎上，輔以高清遙感影像/衛星照片，對波堆藏布江流域古鄉冰期時的冰川作用范圍進行了重建；后進一步運用冰面縱剖面模型對波堆藏布江流域古鄉冰期時古冰川進行了模擬；此外運用AAR和THAR方法對古鄉冰期和現代冰川ELA進行了估算。古鄉冰期時，波堆藏布流域冰川的平均厚度約為360 m，冰川總面積約為2 648 km2，體積約為953 km3，流域內的冰川覆蓋率由古鄉冰期的63%縮減至現在的22.4%；現代冰川ELA約4 455 m，古鄉冰期時ELA約為3 871 m，較現代冰川下降了約584 m。

“古鄉冰期”是李吉均先生主導提出的第四紀冰川重要冰期概念，謹以此文，紀念李吉均院士！

致謝：感謝華東師范大學地理科學學院許可芃在數據處理中提供了幫助。