張哲茵,衣鵬,c,d,陳鵬
(河海大學,a.水文水資源學院;b.水文水資源與水利工程科學國家重點實驗室;c.長江保護與綠色發(fā)展研究院;d.全球變化與水循環(huán)國際合作聯(lián)合實驗室;e.地球科學與工程學院,南京 210098)
青藏高原作為“世界第三極”,平均海拔超過4 000 m,被稱作亞洲水塔,是亞洲眾多河流的發(fā)源地,為超過14 億人口提供水源[1]。由于青藏高原處于海拔高、溫度低的極端環(huán)境,所以其生態(tài)環(huán)境非常脆弱且自我恢復(fù)能力很差[2]。全球氣候變暖以及愈加頻繁的人類活動對青藏高原本就脆弱的生態(tài)環(huán)境產(chǎn)生了巨大威脅。此外,青藏高原變暖速率是全球平均速率的2倍,加上降水集中在春夏,青藏高原極易發(fā)生土壤侵蝕[1,3-5]。青藏高原東南部是長江、怒江上游源區(qū)所在地[6],相比青藏高原的其他地區(qū),青藏高原東南部水系更多、降水更多、地形復(fù)雜,土壤極易發(fā)生侵蝕。有學者預(yù)測至2050年,青藏高原東南部將是中國潛在土壤侵蝕量最高的地區(qū)之一[7]。
土壤碳循環(huán)是全球碳循環(huán)的重要組成部分,土壤中的碳可能會對全球氣候產(chǎn)生影響[8],有學者認為土壤侵蝕會加速土壤中碳的排放[9]。土壤侵蝕會帶來許多危害,不僅會引起生態(tài)功能的減弱以及土地生產(chǎn)力的消退[10],還會對當?shù)氐纳鐣⒔?jīng)濟、文化方面產(chǎn)生影響[11]。除此之外,還可能產(chǎn)生更嚴重的后果,例如由于土壤侵蝕河流中的含沙量增加,大量泥沙在下游河道淤積,河床抬高,河道的泄洪能力降低,洪澇災(zāi)害發(fā)生的風險增加[12]。由此可知,對土壤侵蝕的研究及對土壤動態(tài)的科學管理極為重要。
目前主要的土壤侵蝕測量方法有傳統(tǒng)測量學方法、地球化學法和GIS 技術(shù)方法三大類。傳統(tǒng)測量法存在著測量周期長、成本高、監(jiān)測范圍小的問題,而GIS 技術(shù)方法也存在著范圍不精確、誤差大的問題[13]。地球化學法即通過同位素來測量土壤的侵蝕量,目前常用的同位素有10Be、137Cs、210Pbex等。與傳統(tǒng)測量法相比,地球化學法測量有著周期短、投入少、不改變原始地貌的優(yōu)點[14,15]。137Cs 產(chǎn)生于20 世紀的核試驗,半衰期為30.08年,對土壤細顆粒有很強的親和力,不易被植被吸收,只隨土壤顆粒運動,因此被作為估算土壤侵蝕量的示蹤劑[16,17]。
目前許多學者對青藏高原開展了土壤侵蝕的研究,各類侵蝕測量方法均被使用,有計算總侵蝕量的,也有單獨計算水力侵蝕、風力侵蝕和凍融侵蝕的[5,18]。但是大部分研究都是針對流域尺度或是一個小范圍的區(qū)域(一個盆地或幾個縣或一段公路途經(jīng)的區(qū)域),很少有學者能針對大范圍的區(qū)域進行侵蝕量的研究[1,19,20]。有部分學者使用GIS 技術(shù)和相關(guān)模型對整個青藏高原的水力侵蝕和風力侵蝕狀況進行了計算[1,21],但模型模擬的結(jié)果誤差較大,不同學者使用同一模型計算同一時間段內(nèi)同一區(qū)域的侵蝕模數(shù)的結(jié)果也存在較大差異[2]。也有許多學者使用同位素示蹤法對青藏高原進行了侵蝕計算,但多位于三江源地區(qū)和青藏高原的中部及北部[5,18,22,23],目前并沒有學者對青藏高原東南部的土壤侵蝕狀況使用同位素手段進行分析。本研究使用137Cs 示蹤法分析青藏高原東南部土壤的侵蝕強度和分布。
研究區(qū)域位于青藏高原東南部(91°27′36″E—103°31′12″E,27°8′24″N—34°6′0″N),涉及西藏的東部、青海省的西南部以及四川省的西北部,見圖1。研究區(qū)域的平均海拔高于4 000 m,多年平均降水量在250~620 mm,多年平均地表溫度在5.5~12.6 ℃。區(qū)域內(nèi)水系發(fā)達,主要河流有雅魯藏布江、瀾滄江、金沙江以及眾多支流河道。區(qū)域內(nèi)的主要土壤類型是高山土和淋溶土,按凍土類型分類則有季節(jié)性凍土和少量的多年凍土。研究區(qū)域內(nèi)植被種類以草甸、灌木、針葉林為主。研究區(qū)域覆蓋了多個生態(tài)保護區(qū)的部分區(qū)域,生態(tài)環(huán)境比較脆弱,尤其是在海拔超過5 000 m 的地區(qū)[6]。
圖1 研究區(qū)域樣點分布
根據(jù)中國的三級流域空間分布(https://www.resdc.cn/data.aspx?DATAID=278),把研究區(qū)域分為5 個區(qū)域進行侵蝕強度和侵蝕分布的分析。從東到西分別是區(qū)域一至區(qū)域五,每個區(qū)域中分別包含1條一級河流,從左到右分別為怒江、扎曲、金沙江、雅礱江、大渡河。區(qū)域一的平均高程為4 643 m,多年平均降水量為400 mm,多年平均地溫為8.9 ℃;區(qū)域二的平均高程為4 475 m,多年平均降水量為344 mm,多年平均地溫為8.2 ℃;區(qū)域三的平均高程為4 088 m,多年平均降水量為398 mm,多年平均地溫為10.2 ℃;區(qū)域四的平均高程為4 147 m,多年平均降水量為448 mm,多年平均地溫為8.5 ℃;區(qū)域五的平均高程為3 841 m,多年平均降水量為486 mm,多年平均地溫為9.7 ℃。
2.1.1 樣點設(shè)置及樣品采集 采樣點多位于地面較平整受人類活動影響較小的區(qū)域,以此來反映土壤在自然狀態(tài)下的侵蝕狀況。此外,所設(shè)置采樣點在研究區(qū)域內(nèi)盡量均勻分布,使其能夠較為準確地反映區(qū)域侵蝕狀況。采樣點設(shè)置情況具體見圖1,共43 個點。2019 年夏季對研究區(qū)域內(nèi)43 個采樣點進行了表層0~5 cm土壤的全樣采集,在每個點采集4~5個土樣進行混合裝樣,土樣的間距在1 m 左右。采樣時去除采樣點表層植物根系后,使用環(huán)刀采樣,裝入塑封袋并帶回實驗室分析。對樣品進行137Cs 的測量(包含圖1 中1~26 號樣點),同時測定有機碳(TOC)(包含圖1 中1~29 號樣點)。此外,對所有點進行土壤粒徑的測量。
2.1.2 測定指標與方法
1)137Cs 比活度的測量。在中國科學院南京土壤研究所進行,采用美國ORTEC 公司的高純鍺低本底伽瑪能譜儀(探頭型號:GWL-120-15;多道型號:DSPEC jr2.0)進行測定。具體步驟為將土樣烘干,去除有機質(zhì)、石塊,研磨過200 目篩子,裝入測試樣品盒,密封10 d,直接放入探測器井,測量時間為80 000 s,137Cs 標準由中國原子能研究院提供,最低檢測限為0.12 Bq/kg。
2)TOC 的測量。在河海大學水文水資源與水利工程國家重點實驗室進行,采用德國耶拿分析儀器股份公司生產(chǎn)的Multi N/C2100 型總有機碳/總氮分析儀測量。具體步驟為將樣品凍干,去除石塊草根,研磨過200 目篩子后,稱取約100 mg 的土樣,加入足量的稀鹽酸,經(jīng)過儀器超過1 000 ℃的燃燒測得樣品中TOC 的含量,每隔15 個樣質(zhì)檢一次,誤差控制在10%以內(nèi)。
3)土壤粒徑的測量。使用LS 13320 型粒徑儀,測量范圍在20 nm 至2 000 μm。將土樣烘干,去除植物根系、石子,經(jīng)簡單研磨后過10 目篩子,加入足量的過氧化氫溶液(30%)進行反應(yīng),靜置10 h 以上,測量前攪拌均勻,加入粒徑儀中測量,測量3 次取平均值。
2.1.3 其他數(shù)據(jù)來源 多年平均降水量、多年平均地表溫度數(shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/data/detail/dataCode/A.0012.0001.html)中研究 區(qū)域附近39 個氣象站2000—2019 年的均值。研究區(qū)域內(nèi)的土地覆蓋數(shù)據(jù)、歸一化植被指數(shù)(NDVI)來自中國科學院資源環(huán)境科學與數(shù)據(jù)中心(https://www.resdc.cn/)2000—2019 年的均值,分辨率均為1 km。數(shù)字高程數(shù)據(jù)(DEM)來自地理空間數(shù)據(jù)云(https://www.gscloud.cn/),分辨率為30 m。
通過與未發(fā)生侵蝕或堆積的點(背景點)的同位素活度比較來計算侵蝕量,若采樣點137Cs 的比活度大于背景值,說明該點發(fā)生堆積;若采樣點137Cs 的比活度小于背景值,說明發(fā)生侵蝕[5]。137Cs 隨大氣沉降后進入土壤,被土壤表層的黏土礦物吸附,隨著土壤細顆粒在重力和雨水淋溶作用向下遷移,最終分布深度一般在10 cm 以內(nèi)[24]。部分學者在青藏高原開展137Cs 的研究,該區(qū)域內(nèi)背景點的137Cs 比活度在1 969~2 376 Bq/m2[18,22,25],與齊永青等[26]得到青藏高原的背景值812~2 571 Bq/m2相吻合。本研究采用的背景值為實際測定值。對背景點處0~10 cm的土樣進行測量,背景點處137Cs 的平均比活度為25.19 Bq/kg,根據(jù)土壤干容重一般在1.1~1.6 g/cm3,折算出背景點處的137Cs 的面積活度為2 771~4 030 Bq/m2,考慮到藏東南降水充沛,137Cs 的沉降量較大,因此該背景點的137Cs 比活度較為合理,可以用于后續(xù)分析。
研究區(qū)域內(nèi)137Cs 的分布很淺,在背景點處的分布深度只有8 cm,0~5 cm 的土樣中的137Cs 占總量的90%以上,故樣點處0~5 cm 土壤中的137Cs 能夠反映土壤的侵蝕程度。參考李元壽等[5]的模型,假設(shè)137Cs在非耕作土中呈指數(shù)分布[27]且侵蝕點的剩余部分137Cs分布與未侵蝕點的對應(yīng)部分的分布相同。
基于上述假設(shè),標準剖面中137Cs 比活度的垂直分布函數(shù)如式(1)所示。
式中,f(z)為z(cm)樣地深度處137Cs 的比活度(Bq/kg);a為地表處(z=0)的137Cs 的比活度(Bq/kg),b為137Cs的深度衰減系數(shù)(cm-1)。
假設(shè)土壤容重沿深度方向不變且各樣點的土壤容重相同,則存在一個面積S1,使得面積為S1、高度為1 cm 的土柱的質(zhì)量為1 kg,則該土柱0~zcm的137Cs比活度的表達式如式(2)所示。
采樣點處面積S1的土柱0~5 cm的137Cs 比活度如式(3)所示。
可得:
式中,C0為采樣點0~5 cm的137Cs 平均比活度(Bq/kg),h為侵蝕點的侵蝕深度(cm),參數(shù)a、b由背景點處137Cs 比活度隨深度的分布確定,背景點處137Cs 隨深度的分布具體為:0~2 cm 的平均比活度為34.2 Bq/kg;2~4 cm 的平均比活度為65.85 Bq/kg;4~6 cm 的平均比活度為17.3 Bq/kg;6~8 cm 的平均比活度為8.6 Bq/kg;8~10 cm 的平均比活度為0。從而得到C(0)=0,C(2)=68.4 Bq,C(4)=200.1 Bq,C(6)=234.7 Bq,C(8)=251.9 Bq,C(9)=251.9 Bq,C(10)=251.9 Bq。使用最小二乘法進行曲線擬合,可得a=96、b=0.37,見圖2。在背景點處,令h=0,可得C0=43.7 Bq/kg。即背景點處0~5 cm 土壤的137Cs 比活度為43.7 Bq/kg。
圖2 C(z)-z 的最小二乘擬合
研究區(qū)域內(nèi)137Cs 的比活度的范圍為0~23.75 Bq/kg,由于低于檢測下限而未測出137Cs 的樣點共6個,分別是樣點4、7、13、18、20、26。測出137Cs 的樣點的137Cs 比活度 范圍為2.6~23.75 Bq/kg,平均為10.48 Bq/kg,小于背景值43.7 Bq/kg,說明研究區(qū)域的土壤整體上呈侵蝕狀態(tài)。
區(qū)域一的137Cs 比活度范圍為0~17.96 Bq/kg,平均為7.31 Bq/kg;區(qū)域二的137Cs 比活度范圍為0~13.67 Bq/kg,平均為8.48 Bq/kg;區(qū)域三的137Cs 比活度范圍為0~5.63 Bq/kg,平均為2.82 Bq/kg;區(qū)域四的137Cs比活度范圍為0~18.94 Bq/kg,平均為8.54 Bq/kg;區(qū)域五的137Cs 比活度的范圍為7.38~23.75 Bq/kg,平均為15.57 Bq/kg。
研究區(qū)域內(nèi),可檢測出137Cs 比活度的樣點處的137Cs 比活度均低于背景點值,表明所有采樣點均存在土壤侵蝕現(xiàn)象。未測出137Cs 的樣點是由于該點處侵蝕強度較大,137Cs 比活度低于監(jiān)測下限,除未能測出137Cs 的樣點,區(qū)域一的土壤侵蝕厚度范圍為2.41~7.63 cm,平均為4.52 cm;區(qū)域二的土壤侵蝕厚度范圍為3.14~4.82 cm,平均為3.90 cm;區(qū)域三的土壤侵蝕厚度為5.54 cm;區(qū)域四的侵蝕厚度范圍為2.26~7.05 cm,平均為4.29 cm;區(qū)域五的侵蝕厚度范圍為1.65~4.91 cm,平均為3.23 cm。區(qū)域一、區(qū)域二、區(qū)域三、區(qū)域四均存在低于137Cs 測定下限的樣點,故實際侵蝕厚度會比計算的結(jié)果大。
對于未測出137Cs 的樣點,由背景點處137Cs 的最大分布深度為8 cm,未測出的樣點由于侵蝕強度過大導(dǎo)致樣點處137Cs 的比活度低于監(jiān)測下限,可得樣點處的侵蝕厚度大于8 cm,這里假設(shè)其侵蝕厚度為8 cm。土壤干容重一般在1.1~1.6 g/cm3,本試驗假設(shè)研究區(qū)域所有采樣點的土壤干容重相同,即都為1.35 g/cm3。通過換算可將137Cs 示蹤法計算的侵蝕厚度轉(zhuǎn)換為1963 年以來的平均侵蝕模數(shù),區(qū)域一的平均侵蝕模數(shù)為13.0 t/(hm2·年),區(qū)域二的平均侵蝕模數(shù)為11.4 t/(hm2·年),區(qū)域三的平均侵蝕模數(shù)為16.3 t/(hm2·年),區(qū)域四的平均侵蝕模數(shù)為12.1 t/(hm2·年),區(qū)域五的平均侵蝕模數(shù)為7.8 t/(hm2·年)。根據(jù)137Cs 示蹤法計算的侵蝕模數(shù),5 個區(qū)域按平均侵蝕模數(shù)從大到小依次為區(qū)域三、區(qū)域一、區(qū)域四、區(qū)域二、區(qū)域五。
研究區(qū)域的年平均溫度大于張建國等[28]提出的西藏地區(qū)凍融侵蝕區(qū)的溫度下界-2.5 ℃,故研究區(qū)域內(nèi)凍融侵蝕力不是主要侵蝕力。青藏高原的風蝕力空間上從東南向西北逐漸增大[29],研究區(qū)域位于青藏高原東南部,風蝕力較弱,加上區(qū)域內(nèi)降水豐富,多年平均降水量為250~620 mm,故本區(qū)域的主要侵蝕類型為水力侵蝕。
水力侵蝕的影響因素有地形因素、氣象因素、植被覆蓋和土壤性質(zhì)[4,30]。將26 個計算侵蝕厚度的點根據(jù)植被類型分為草甸和林地(針葉林、灌木)兩類。去除未測出137Cs 的樣點,余下的樣點進行侵蝕厚度與其影響因子的相關(guān)性分析,對滿足正態(tài)分布的因子使用Pearson 相關(guān)系數(shù),對不滿足正態(tài)分布的因子使用Spearman 相關(guān)系數(shù),曲線擬合計算相關(guān)系數(shù)(R),具體結(jié)果見圖3。
對坡度、高程與土壤侵蝕厚度的關(guān)系進行分析,見圖3a 和圖3b。有研究表明,拉薩河流域的水蝕強度在坡度為15°~25°的區(qū)域最大,其次為25°~35°[31]。在本研究區(qū)域內(nèi)也發(fā)現(xiàn)了類似的規(guī)律,如圖3a,虛線橢圓外的3 個點離點群較遠,去除這3 個點后使用二次多項式進行擬合(R=0.45),在土壤侵蝕厚度最大的點所對應(yīng)的坡度約為18°。也就是說,在研究區(qū)域內(nèi),土壤侵蝕厚度隨著坡度的增加先增大后減小,在坡度約為18°的區(qū)域土壤侵蝕厚度最大。土壤侵蝕厚度與高程之間沒有明顯的相關(guān)性(圖3b)。
對多年平均降水量、地溫、風速與土壤侵蝕厚度進行相關(guān)性分析。多年平均降水量與土壤侵蝕厚度之間的正相關(guān)關(guān)系并不顯著(圖3c),與研究區(qū)域以水力侵蝕為主的情況不符。這是因為并非所有的降水都會產(chǎn)生土壤侵蝕[32],只有產(chǎn)生徑流的降雨才會產(chǎn)生侵蝕[33],所以多年平均降水量與侵蝕量之間沒有顯著的相關(guān)性。在植被類型為草甸的樣點處,土壤侵蝕厚度與多年平均地溫存在顯著的正相關(guān)性(R=0.76,P<0.05);而在植被類型為林地的樣點處,土壤侵蝕厚度與地溫之間的正相關(guān)性并不顯著(圖3d)。這可能與融水有關(guān),也有可能與風化作用相關(guān),但根據(jù)本次研究的數(shù)據(jù)很難判斷。土壤侵蝕厚度與多年平均風速沒有發(fā)現(xiàn)明顯的關(guān)系(圖3e),原因是區(qū)域內(nèi)風力侵蝕作用較弱。
植被覆蓋度的提升能顯著加強土壤的水土保持能力[34],有學者指出在青藏高原的黃河源區(qū),草甸土的土壤侵蝕模數(shù)與植被覆蓋度之間有顯著的負相關(guān)關(guān)系[5]。這里采用2000—2019 年歸一化植被指數(shù)(NDVI)的平均值計算植被覆蓋度(fc),覆蓋度計算見公式(5)。
式中,NDVImin為NDVI的最小值,即純裸土像元的NDVI;NDVImax指NDVI的最大值,即純植被像元的NDVI。
在本研究區(qū)域內(nèi)草甸和林地的侵蝕厚度與植被覆蓋度間的負相關(guān)性都不顯著(圖3f)。其原因可能有:基于分辨率為1 km的NDVI數(shù)據(jù)計算的植被覆蓋度存在誤差;植被分類不夠具體,不同植被的根系分布不同,保持水土的能力也不同。
探討土壤侵蝕厚度與土壤粒徑、有機質(zhì)含量之間的關(guān)系。研究區(qū)域內(nèi),林地處土壤侵蝕厚度與中值粒徑有顯著的負相關(guān)性,而在草甸處,二者的負相關(guān)性不顯著,見圖3g。本研究采用的是0~5 cm 的土壤測量粒徑,這部分土壤會受到侵蝕的影響,侵蝕作用會使得土體顆粒破碎、分散,破壞表層土壤的結(jié)構(gòu)[35]。這部分土不能準確地反映土壤在粒徑方面的差異,因此這里不能判斷土壤粒徑與侵蝕厚度之間的關(guān)系。
土壤中的有機質(zhì)主要存在于黏粒和粉粒中,有機質(zhì)會降低黏粒和粉粒的膠結(jié)能力,從而使土壤團聚體的破壞率增高[4]。因此,土壤有機質(zhì)含量越高,土壤越容易發(fā)生侵蝕。由圖3h 可知,研究區(qū)域內(nèi),在林地,TOC 與侵蝕厚度之間存在顯著的負相關(guān)性(R=-0.62,P<0.05),而在草甸處二者的負相關(guān)性不顯著。在草甸處,TOC 與侵蝕厚度間負相關(guān)性不顯著的原因可能是樣點較少,TOC 的變幅比較小,不容易發(fā)現(xiàn)TOC 與侵蝕厚度間的相關(guān)性。
研究區(qū)域內(nèi)137Cs比活度的范圍為0~23.75 Bq/kg,所有采樣點都發(fā)生土壤侵蝕。研究區(qū)域內(nèi)5 個區(qū)域均發(fā)生輕度侵蝕,按侵蝕模數(shù)從大到小分別為區(qū)域三、區(qū)域一、區(qū)域四、區(qū)域二、區(qū)域五。
研究區(qū)域內(nèi),侵蝕厚度隨坡度的增加先增大后減小,在坡度為18°附近侵蝕厚度最大。在侵蝕厚度與多年平均降水量、風速、植被覆蓋度間沒有觀察到顯著的相關(guān)性。林地土壤的侵蝕厚度與TOC 間存在顯著的負相關(guān)性(P<0.05),草甸土壤處二者的負相關(guān)性不顯著。