基于Budyko假設的渭河徑流變化歸因識別

張麗梅,趙廣舉,,*,穆興民,,高 鵬,,孫文義,

1 西北農林科技大學黃土高原土壤侵蝕與旱地農業(yè)國家重點實驗室,楊凌 712100 2 中國科學院水利部水土保持研究所,楊凌 712100

水資源是社會、經濟和生態(tài)環(huán)境可持續(xù)發(fā)展的重要基礎。河川徑流是水資源最重要的存在形式和組成部分之一,是生產生活用水的最主要來源[1]。受到氣候變化和人類活動的共同影響,全球許多河流的河川徑流量發(fā)生了顯著變化,嚴重威脅著區(qū)域水資源狀況[2-3],揭示徑流變化的主要驅動力是全球水資源預測的關鍵,也是維持淡水資源可持續(xù)利用的基礎[4-5]。

渭河流域位于黃河中游,是甘肅和陜西的母親河,也是關中地區(qū)工農業(yè)取水用水的主要來源,是該地區(qū)經濟社會發(fā)展的保障[6]。近年來,渭河徑流量發(fā)生了顯著變化,其變化原因引起了相關管理部門與研究人員的極大關注[7]。郭愛軍等[8]研究發(fā)現(xiàn)人類活動對徑流變化的影響作用處于主導地位,氣候變化的影響呈不明顯的增加趨勢。畢彩霞等[9]采用水量平衡法以渭河流域的華縣站為例研究1958—2011年徑流減少的影響因素,結果表明:降雨變化和人類活動的貢獻率分別為49.0%和51.0%。Guo等[10]研究渭河上游徑流變化的驅動因素,結果表明,人類活動對徑流的貢獻占59%—77%。Huang等[11]分析了渭河不同年代際的徑流變化及其原因,發(fā)現(xiàn)較1960s氣候變化在20世紀70年代和90年代占主要因素,而人類活動在20世紀80年代和2000年之后起主導作用。盡管前人已經開展了許多關于渭河流域徑流變化特征的研究,并識別其變化的主要驅動力,但研究結果仍存在較大的差異。

流域蒸散發(fā)在年尺度上受降水和蒸發(fā)能力兩個因素控制,基于此理論,Bukyko提出了流域水量與能量(太陽輻射產生的能量)耦合平衡方程的構想,即Budyko假設[12]。目前,許多研究人員根據(jù)氣象和徑流資料在全球多個地區(qū)驗證了Budyko假設的成立[13- 15]。研究發(fā)現(xiàn),影響流域水量平衡的因素除了蒸散和降水外,還包括以下兩類因素的影響:一類是下墊面對流域的影響,包括土壤屬性、地形、植被變化等[14-17];另外一類為除蒸散和降水量以外的氣候因子,包括降雨深度、降雨頻次等[18]。因此,近年來許多學者提出了一系列反映下墊面因素的Budyko經驗公式,如Zhou等[19]和楊大文等[20]分別就下墊面參數(shù)進行了推導和驗證,研究表明Budyko解析式方程遵循能量和水分邊界條件?；贐udyko水熱耦合平衡理論的水量平衡法對定量區(qū)分氣候和下墊面變化對流域徑流的影響具有廣闊的應用前景[21-22],國內學者楊大文等[20]和孫福寶等[23]已證實了Budyko水熱耦合平衡方程在黃河流域的適用性。

本研究在深入辨析渭河流域徑流、降水、蒸發(fā)變化的基礎上,采用Budyko水熱耦合平衡方程估算流域年徑流變化的氣候和人類活動的彈性系數(shù),量化降雨、蒸發(fā)以及下墊面變化對徑流變化的影響,識別渭河徑流銳減的原因,以期為渭河流域水資源規(guī)劃以及流域綜合治理提供理論支持。

1 研究區(qū)與數(shù)據(jù)獲取

1.1 研究區(qū)概況

渭河是黃河水量最多、面積最大的一級支流,流域面積13.48萬km2(103°55′—110°20′E,33°40′—37°25′N)(圖1)。渭河發(fā)源于甘肅省渭源縣鳥鼠山,由西向東流經甘肅、寧夏、陜西三省,于陜西省潼關縣附近匯入黃河。渭河干流全長為818 km,支流較多;第一大支流涇河,河長455.1 km,流域面積為4.54萬km2;第二支流為北洛河,河長680 km,流域面積為2.70萬km2。渭河流域多年平均降雨量約為550.8 mm(1958—2015年),在空間分布上呈現(xiàn)南多北少,山區(qū)降雨量大于河谷盆地,且多年平均潛在蒸散量為995.0 mm(1958—2015年)。

1.2 數(shù)據(jù)來源與處理

本文選取渭河流域咸陽、張家山、狀頭3個水文站點(表1),年徑流量數(shù)據(jù)來源于中國水文水資源科學數(shù)據(jù)共享網和《黃河泥沙公報》等。本文選擇了渭河流域22個氣象站(圖1),數(shù)據(jù)采用中國氣象數(shù)據(jù)共享網(http://data.cma.cn)1958—2015年逐月氣象資料,包括降水量、平均氣溫、最高和最低氣溫、日照時數(shù)、平均風速、相對濕度等。本文采用FAO修正的Penman-Monteith公式[24]計算潛在蒸散發(fā),計算公式如下:

(1)

式中,ET0為潛在蒸散量(mm/d);為飽和水汽壓曲線斜率(kPa/℃);Rn為太陽凈輻射(MJ m-2d-1);G為土壤熱通量(MJ m-2d-1);γ為干濕常數(shù)(kPa/℃);U2為2 m高處的風(m/s);T為平均氣溫(℃);es為平均飽和水汽壓(kPa),ea為實際水汽壓(kPa),(es-ea)為飽和水汽壓差(kPa)。根據(jù)月氣象資料計算得到潛在蒸發(fā)量,并統(tǒng)計不同季節(jié)和年潛在蒸發(fā)量。季節(jié)的劃分采用氣象季節(jié)即3—5月為春季,6—8月為夏季,9—11月為秋季,12—次年2月為冬季。

表1 水文站基本情況

圖1 渭河流域水文站點和氣象站點分布圖Fig.1 Locations of hydro-climatic stations in the Wei River basin

2 方法

2.1 Mann-Kendall趨勢檢驗

Mann-Kendall趨勢檢驗法,在氣候和水文序列變化趨勢分析中應用廣泛,可以確定氣候及水文時間序列變化趨勢的顯著性[25]。

2.2 雙累積曲線

利用徑流量與降雨量的雙累積曲線,研究直線斜率的變化過程,直線斜率發(fā)生明顯偏離的點對應徑流量開始發(fā)生顯著變化的年份[26]。

2.3 徑流變化歸因識別

基于Budyko水熱耦合平衡理論的水量平衡法是區(qū)分氣候和下墊面變化對徑流變化貢獻的有效方法。該方法與傳統(tǒng)的數(shù)理統(tǒng)計方法相比,其物理意義更加明顯,所用參數(shù)相對易獲取,計算方法較簡單,是研究河川徑流變化年序列尺度上一種有效且理想的分析方法,在研究氣候和下墊面變化對徑流影響程度時被廣泛應用[27]。

2.3.1 水量平衡方程

流域水量平衡方程為:

R=P-ET-ΔS

(2)

式中,R為徑流深(mm);P為降水量(mm);ET為實際蒸散量(mm);ΔS為儲水量變化(mm)。

在流域尺度上,徑流深和降雨量可通過實際觀測獲得,實際蒸散量可采用Budyko假設計算獲取。Choudhury[13]和Yang[28]等基于Budyko假設,推導出流域水熱耦合平衡方程,表達式如下:

(3)

式中,ET0為年平均潛在蒸散量(mm);ω為下墊面特征參數(shù)。

下墊面特征參數(shù)(ω)是控制Budyko曲線形狀的參數(shù),其主要反映下墊面特征對流域水平衡的綜合影響[29-30],在本流域中主要表現(xiàn)為植被變化對流域蒸散發(fā)的影響。分析長時間水文序列,儲水量的變化(ΔS)一般假定為零值[27],即在假定儲水量不變的情況下,對多年平均尺度及年際尺度的彈性系數(shù)進行計算。結合(2)、(3)式,水量平衡方程可以表示為如下公式:

(4)

式中,已知R,P和ET0可求算ω值。

2.3.2 敏感性分析

由于水量平衡方程(4)式可表示為R=f(P,ET0,ω),徑流對特定獨立變量x的彈性系數(shù)可采用下式表示:

(5)

式中,εxi是徑流對特定獨立變量xi的彈性系數(shù)(Elasticity coefficent),xi表示P,ET0或ω。假設:

(6)

各變量的彈性系數(shù)計算如下[27]:

(7)

(8)

(9)

某個變量彈性系數(shù)的正值表示徑流深(R)隨著該變量的增加而增加,負值則表示徑流深(R)隨著該變量的增加而減小,使用(7)—(9)式分別求得,可以獲得徑流對降水的彈性系數(shù)(εP)、徑流對潛在蒸散發(fā)的彈性系數(shù)(εET0)和徑流對下墊面的彈性系數(shù)(εω)。

2.3.3 下墊面變化對徑流變化的貢獻

根據(jù)突變點將研究時期劃分為不同時段,時段1的多年平均徑流深為R1,時段2的多年平均徑流深為R2,從時段1到時段2的年徑流的變化可以用徑流前后兩時段的多年平均徑流深之差(dR)表示,即:

dR=R2-R1

(10)

同理,降水量(dP)、潛在蒸散發(fā)(dET0)和下墊面(dω)的變化表示為:

dP=P2-P1

(11)

dET0=ET02-ET01

(12)

dω=ω2-ω1

(13)

由一定因子引起的徑流變化可以通過因子變化及其偏導數(shù)的乘積來估計。因此,每個因子對徑流變化的貢獻可以使用以下微分方程來計算:

(14)

式中,dR′為計算求得的徑流深變化。上式可化簡為:

dR′=dRP+dRET0+dRω

(15)

dRP、dRET0和dRω分別是氣候變化(P和ET0)和下墊面(ω)的變化引起的徑流變化。

將(5)式代入,可表示為:

(16)

每個因素對徑流變化的相對貢獻可以計算如下:

(17)

式中,xi表示P,ET0或ω,Cxi表示各因子對徑流變化的貢獻率(Contribution rate)。

3 結果與分析

3.1 渭河徑流、降雨及潛在蒸發(fā)變化趨勢分析

表2為渭河流域咸陽站、張家山站、狀頭站1958—2015年徑流深的統(tǒng)計特征。由表可知,咸陽站多年平均徑流深較張家山和狀頭站高,且極值差異顯著,各站年徑流深的極大值均出現(xiàn)在1964年。采用非參數(shù)Mann-Kendall趨勢檢驗對各站徑流深進行趨勢分析得出:渭河流域各站徑流深均呈顯著減少趨勢,且均達到0.01的顯著性水平。圖2為渭河流域各水文站1958—2015年徑流深變化。通過線性趨勢分析發(fā)現(xiàn),咸陽站、張家山站、狀頭站的徑流深減少速率分別為:-1.520、-0.501、-0.322 mm/a。

表2 各水文站1958—2015年徑流深變化特征

顯著性水平為0.05時,統(tǒng)計檢驗臨界值為±1.96;顯著性水平為0.01時,統(tǒng)計檢驗臨界值為±2.58;顯著性水平為0.1時,統(tǒng)計檢驗臨界值為±1.28

圖2 咸陽站、張家山站、狀頭站1958—2015年徑流深變化趨勢Fig.2 Variation of runoff depth at three stations during 1958—2015

圖3為渭河各流域1958—2015年降雨量和潛在蒸散量的年際變化過程線,并采用MK檢驗分析渭河干流區(qū)域、涇河流域、北洛河流域的降雨量和潛在蒸散量1958—2015年變化趨勢,結果表明(表3):各流域降雨量的檢驗統(tǒng)計量為負,但顯著性水平均未達到0.1,減少趨勢不顯著;潛在蒸發(fā)的檢驗統(tǒng)計量呈正值,但顯著性水平均未達到0.1,增加趨勢不顯著。

3.2 渭河徑流量變化的歸因分析

3.2.1 渭河干支流徑流量突變特征

渭河流域降雨-徑流雙累積曲線如圖4所示。從圖中可以發(fā)現(xiàn),渭河干流降水徑流累積曲線于1971年以及1993年發(fā)生明顯偏移,涇河和北洛河流域分別在1997、1993年發(fā)生偏移。降水-徑流的雙累積曲線在發(fā)生偏移前的階段一般視作河川徑流量未受人類活動干擾或干擾較小的時期,即基準期[31]。因此,可以確定渭河干流、涇河和北洛河流域的基準期分別為1958—1970年、1958—1996年和1958—1992年,圖中基準期后的累積曲線的斜率發(fā)生了一定程度的改變,曲線斜率變小,說明徑流相對于降雨有了減少趨勢,并以此將渭河流域徑流序列劃分為不同階段,具體各個站點的劃分時段見表4。

表3 渭河流域降雨量和潛在蒸發(fā)量MK檢驗分析結果

圖4 渭河流域徑流-降雨雙累積曲線Fig.4 Double mass curve analysis between cumulative runoff and precipitation in the Wei River basin

流域 Watershed時期 PeriodET0/mmR/mmP/mmωR/PET0/P彈性系數(shù)Elasticity coefficientεPεET0εω渭河干流1958—1970891.9 129.7 674.6 2.08 0.19 1.32 2.51-1.51-1.321971—1992852.9 84.5 637.7 2.67 0.13 1.34 3.06-2.06-1.531993—2015912.7 50.0 603.7 2.97 0.08 1.51 3.51-2.51-1.99涇河1958—1996964.8 42.2 508.1 2.33 0.08 1.90 3.02-2.02-2.251997—20151039.2 22.1 495.0 2.72 0.04 2.10 3.50-2.50-2.83狀頭1958—19921085.7 36.6 488.0 2.14 0.08 2.23 2.89-1.89-2.471993—20151147.7 24.9 467.2 2.26 0.05 2.46 3.06-2.06-2.83

ET0:年平均潛在蒸散量,Annual average potential evapotranspiration;R:徑流深,Runoff depth;P:降水量,Precipitation;ω:下墊面特征參數(shù),Underlying surface feature parameters;R/P:徑流系數(shù),Runoff coefficient;ET0/P,干旱指數(shù),Drought index;εP:降雨量的彈性系數(shù),Elasticity coefficent of precipitation;εET0:潛在蒸散發(fā)的彈性系數(shù),Elasticity coefficents of potential evapotranspiration;εω:下墊面特征參數(shù)的彈性系數(shù),Elasticity coefficents of underlying surface feature parameters

3.2.2 徑流對氣候要素和下墊面變化的敏感性分析

表4列出了各流域各時段水文氣候特征值和徑流對潛在蒸散發(fā)、降水以及地表參數(shù)的彈性系數(shù)。各流域變化期的年平均降水、年均徑流深均較基準期減少;除渭河干流區(qū)1971—1992年外,其余時段各流域年平均潛在蒸散量較基準期增加;干旱指數(shù)(ET0/P)較基準期增大,徑流系數(shù)(R/P)較基準期減小。

總體而言,徑流與ET0、ω呈負相關,但與P呈正相關。3個彈性系數(shù)的絕對值最大的為P,中間值為ω,最小的為ET0。ET0的彈性系數(shù)范圍為-2.51—-1.51,P為2.51—3.51,ω為-2.83—-1.32。這些范圍表明,ET0,P或ω增加1%將導致1.51%—2.51%的下降,2.51%—3.51%的增加或1.32%—2.83%的徑流減少。渭河干流區(qū)域的干旱系數(shù)均較小,表明渭河干流氣候較北洛河和涇河相對濕潤。北洛河流域干旱指數(shù)大于2.20,徑流系數(shù)較小,年徑流深約32.0 mm,屬溫帶大陸性季風氣候,日照充足,且雨量偏少。由表4知,彈性系數(shù)的絕對值在時間上均有增大趨勢。

圖5為1958—2015年各流域徑流對氣候和地表變化的年彈性系數(shù),其反映了氣候與地表參數(shù)對徑流影響不同時期的影響。由圖可知,渭河干流和涇河流域εp和εET0的絕對值顯著增加(P<0.01),北洛河流域的增加趨勢不顯著,這表明渭河干流和涇河流域徑流深對氣候變化的敏感性較高,但北洛河流域的敏感性較差。3個流域的εω的絕對值均顯著增加(P<0.01),這表明徑流深對下墊面條件的變化更為敏感。

圖5 潛在蒸散量、下墊面參數(shù)、降雨量的彈性系數(shù)年際變化Fig.5 Relationship between runoff depth and PE, land surface parameters and precipitation

3.2.3 渭河徑流變化的歸因識別

氣候變化(降雨和潛在蒸散發(fā))和下墊面變化(參數(shù)ω)對徑流變化的影響程度如表5所示。由表知,計算求得的徑流深變化(dR′)與實際徑流深的變化(dR)相差很小,這表明本文在評估相關環(huán)境因素對徑流變化的貢獻時所用的方法有效。降水量減少(或增加)和潛在蒸散發(fā)的增加(或減小)以及下墊面特征參數(shù)ω值增大對徑流的減少有正(負)貢獻。

由表5看出,渭河流域在不同時期、不同流域降雨、潛在蒸散量的變化和人類活動對徑流變化影響程度不同。人類活動期與基準期相比,年降水量和徑流量均顯著減少,下墊面參數(shù)增加,除渭河干流1971—1992年均潛在蒸發(fā)量減少外,其他流域的年均潛在蒸發(fā)量人類活動期較基準期均增加;不同流域的人類活動期下墊面的變化均為徑流變化的主導因素;渭河干流、涇河流域下墊面的變化對徑流的影響均大于60%,即渭河干流和涇河的年徑流深對下墊面的變化更為敏感;涇河、北洛河流域對潛在蒸發(fā)因子變化的敏感度較高,潛在蒸散量的變化對徑流變化也起到了重要作用。

綜上,下墊面變化是渭河流域徑流減少的主要影響因素,降雨量變化次之,潛在蒸散發(fā)的影響較小,下墊面變化對徑流量減少的影響程度在渭河干流和涇河較劇烈。

表5 渭河徑流變化的歸因識別

dRP:降雨量P引起的徑流變化,Runoff changes caused by precipitation;dRET0:潛在蒸散量ET0引起的徑流變化,Runoff changes caused by potential evapotranspiration;dRω:下墊面ω引起的徑流變化,Runoff changes caused by underlying surface;dR:徑流深之差,Runoff depth difference;dR′:計算求得的徑流深變化,Runoff depth change calculated;δ:dR′與dR的差值,Difference between dR′ and dR;CP:降雨對徑流變化的貢獻率,The contribution of precipitation to runoff changes;CET0:潛在蒸散量對徑流變化的貢獻率,The contribution of potential evapotranspiration to changes in runoff;Cω:下墊面對徑流變化的貢獻率,The contribution of underlying surface to changes in runoff

4 討論

氣候變化與人類活動是影響河川徑流量變化的兩個最重要的因素。分析渭河流域22個氣象站點氣象要素,發(fā)現(xiàn)降雨量在1958—2015年期間呈非顯著性減少,而潛在蒸散發(fā)顯著增加,這表明流域近60年來呈現(xiàn)暖干化趨勢。降雨、潛在蒸發(fā)變化并不顯著,二者對渭河徑流量的減少的貢獻占比約為18.28%—62.12%,且在時間和空間上存在顯著差異。

近幾十年,黃土高原經過大規(guī)模的水土流失治理,土壤侵蝕量已顯著降低,大部分河流的徑流量和輸沙量顯著減少[32]。渭河流域除受水土保持措施的影響,土地利用、水庫、工農業(yè)取水等對徑流的減少起到很大的作用。表6為渭河流域水土保持措施的統(tǒng)計信息[7]。由表可知,截至2006年,流域累積修建梯田206萬hm2,淤地壩壩控面積1.41萬hm2。大規(guī)模的水土保持措施大大改變了地表產匯流特征與路徑。1999年,國家實施“退耕還林草”,渭河流域林草地面積大幅度上升,林草地的持水性能和截留量不斷上升,導致徑流量不斷減少[33]。此外,流域內共有水庫302座,總庫容達到27.3億m3,蓄、引、提工程有效灌溉面積121萬hm2[9]。水利工程建設和工農業(yè)耗水量的不斷上升導致人類活動對渭河徑流的影響越來越劇烈,致使徑流量不斷減少。壩庫工程建成初期,其蓄水對河川徑流的年內分配影響較大,然而在多年平均狀況下,導致整個流域蒸發(fā)量增加?；贐udyko水熱耦合平衡理論的基本假設是流域在某一段時間內儲水變化量ΔS為零,但該計算方法忽略了流域水利工程尤其是水庫和淤地壩對徑流的攔蓄。盡管相關研究將人類活動的影響作為一個整體考慮[20,27],但部分研究表明下墊面參數(shù)(ω)的物理意義和數(shù)學表達尚不準確,未來研究仍需剔除壩庫工程攔蓄的影響。

表6 渭河流域水土保持措施信息表

5 結論

本文選取渭河流域為研究區(qū),采用Mann-Kendall檢驗、雙累積曲線等分析渭河流域1958—2015年的氣候、水文要素變化,應用基于Budyko水熱耦合平衡方程量化氣候變化和人類活動對徑流變化的貢獻,主要結論如下:

(1)Mann-Kendall趨勢檢驗法對渭河流域咸陽、張家山、狀頭3個控制水文站的徑流和氣候資料趨勢分析得出,徑流深呈顯著減少趨勢,減少速率分別為:-1.520、-0.501、-0.322 mm/a。降雨量和潛在蒸散量分別呈減少和微弱增加趨勢,但均不顯著。

(2)采用徑流—降雨雙累積曲線識別渭河各流域徑流階段性變化特征:渭河干流區(qū)域、涇河流域和北洛河流域徑流變化的基準期分別為1958—1970年、1958—1996年和1958—1992年。依據(jù)Budyko水熱平衡方程計算徑流對各因子彈性系數(shù),結果顯示徑流深對下墊面條件變化更為敏感,人類活動期渭河干流和涇河流域下墊面條件對徑流變化的貢獻率均在60%以上,北洛河流域在1993—2015年下墊面條件對徑流變化的貢獻率為37.88%。總體上,下墊面變化是渭河流域徑流銳減的主要因素,降雨因素次之,潛在蒸散發(fā)的影響較小。