干旱與濕潤區流域時變水熱耦合參數的歸因分析

趙香桂,黃生志,*,趙 靜,韓知明,魏曉婷,黃 強,鄧銘江, 2

1 西安理工大學 西北旱區生態水利國家重點實驗室,西安 710048 2 新疆寒旱區水資源與生態水利工程研究中心,烏魯木齊 830000

隨著全球氣候變暖,黃河和長江流域地區氣候發生了明顯變化；同時,兩區域人類活動(如城市化、工業化、退耕還林還草和農田灌溉等)也在不斷加劇,使得流域下墊面條件發生改變,進而影響降水在蒸發、下滲和徑流之間的分配,從而影響水循環過程[1- 3]。無定河流域是黃河中游的典型一級支流,自1999年以來實施的水土保持和退耕還林還草等生態修復工程實施以來。流域植被覆蓋條件明顯改善,但也造成徑流量急劇減少,從而使得生態環境和經濟發展之間供需水矛盾愈發突出[4-6]。同樣地,地處長江流域的漢江上游地區經濟的快速發展及引水工程的逐年增多,也加劇了徑流的衰減[3],水矛盾問題日益顯著。因此,在復雜的變化環境下,研究不同氣候區流域水文過程的變化并探究其影響因素,有助于深刻認識不同氣候區水文過程的變化特征。

氣候和下墊面條件的變化是影響流域水文循環的主要因素,定量分析兩者對徑流和流域水熱耦合動態變化的影響是當前水文研究的熱點和難點。目前較為常用的研究方法主要是水文模型模擬和Budyko水熱耦合理論[3,7]。水文模型雖然能有效模擬流域水文過程,但對輸入數據要求高,模型的不確定性較大,故模擬結果存在較大誤差[8- 9]。而Budyko水熱耦合理論因包含表征流域下墊面特征的參數,使模型具有一定物理意義,且計算過程相對簡單,被廣泛應用于水文領域的研究中[10-11]。Patterson等[12]應用Budyko理論評估了氣候和直接的人類活動對美國南大西洋徑流量的影響；Xu等[13]基于Budyko理論對海河流域徑流的減小進行了歸因分析；Huang等[14]使用Budyko假設和SVM模型量化了氣候和人類活動對徑流變化的貢獻率；楊大文等[15]基于Budyko理論,以黃河流域38個典型子流域為研究對象對徑流變化進行了歸因分析；張麗梅等[16]基于Budyko框架,估算了渭河徑流變化對各驅動因素的彈性系數,進而定量分解了氣候變化和人類活動對徑流變化的貢獻率；上述基于Budyko框架的研究主要聚焦于徑流的定量歸因分析,側重于水熱平衡理論的應用,而關于Budyko框架本身的研究甚少,故本研究擬開展對Budyko框架中水熱耦合參數n的定量歸因研究。此外,以往關于Budyko框架理論與應用的研究中,表征下墊面特征的水熱耦合參數均被視作恒定值,即認為其在多年尺度上保持不變。但流域下墊面條件往往存在年際或年內尺度上的變化,故基于恒定參數所模擬的流域水熱狀況不能精準的反應流域的水文演變規律。因此,本文探究了變化環境下流域水熱耦合參數n的動態變化規律,并定量分析了氣候變化和下墊面條件的改變對其動態變化的驅動機制。此外,干旱區和濕潤區不僅氣候條件不同,植被覆蓋度和農田灌溉情況也大不相同,以往開展的相應研究也未在不同氣候區進行歸因對比分析。

綜上所述,水熱耦合狀況對氣候變化和下墊面條件改變的響應是一個較復雜的過程。鑒于此,本研究選取氣候類型差異顯著且下墊面條件變化較大的兩個流域——無定河和漢江上游為研究對象,估算干旱與濕潤區流域的時變水熱耦合參數,并通過逐步多元回歸模型、敏感性和貢獻率分析,定量揭示水熱耦合控制參數演變的歸因,并將干旱與濕潤區流域作對比分析,揭示不同氣候區水文循環過程的聯系與區別,研究結果可為深入探究流域水文過程的演變、指導流域水資源合理開發利用和生態環境的保護提供科學依據。

1 研究區域與數據

1.1 研究區概況

無定河流域是黃河中游河口至龍門區間最大的支流,位于陜西省北部,發源于定邊縣白于山,流經定邊、靖邊、米脂、綏德等縣,于清澗縣河口村匯入黃河。無定河全長491 km,流域面積3.03萬km2,流域地理位置處于37°02′—39°00′N,107°47′—110°34′E,流域出口控制站為白家川水文站,集水面積為2.97萬km2,河道平均比降1.8‰。流域氣候屬于溫帶大陸性干旱半干旱季風氣候類型,多年平均氣溫為8.9℃,多年平均降水量為369.9 mm,且降水量空間分布為自東南向西北遞減,年內分布極其不均。

漢江是長江最大的支流,發源于秦嶺南麓,流經陜西和湖北,于武漢市漢口龍王廟匯入長江。漢江流域面積15.9萬km2,干流全長1577 km,通常分上、中、下游3段。丹江口以上為上游,河長925 km；丹江口至鐘祥為中游,河長270 km；鐘祥以下為下游,河長約382 km。該流域位于北緯31°41′—34°11′,東經106°05′—109°22′之間,本文選取安康站以上的漢江上游為研究區域,該區域全長426 km,集水面積3.86萬km2。流域氣候屬于亞熱帶濕潤性季風氣候類型,四季分明,雨量充沛。流域多年平均氣溫為13.4℃,多年平均降水量912.4 mm,70%左右的降水集中在5—9月。

圖1 研究區地理位置Fig.1 Location of the study area

1.2 數據來源

無定河流域1970—2013年白家川水文站的徑流觀測數據來源于黃河流域水文年鑒,氣象數據來源于中國氣象科學數據共享服務平臺(http://data.cma.cn/)提供的4個氣象站(榆林、橫山、綏德、靖邊)的地面氣候數據集,主要包括日降雨量(P)、日平均氣溫(T)、日最高氣溫(Tmax)、日最低氣溫(Tmin)、相對濕度(RH)、平均風速(u2)以及日照時數(S)資料。用ArcGIS 10.2軟件中泰森多邊形工具將點數據轉換為面數據。

漢江上游徑流控制站點安康水文站,觀測數據來源于長江流域水文年鑒；選取流域內9個氣象站點(略陽、留壩、太白、漢中、佛坪、石泉、鎮巴、安康、鎮坪),并于中國氣象科學數據共享服務平臺獲得P、T、Tmax、Tmin、RH、S及u2資料,同樣采用泰森多邊形將其處理為面數據。

本研究使用的植被覆蓋指數(NDVI)是由美國國家海洋大氣局發布的Global Inventory Modeling and Mapping Studies Normalized Difference Vegetation Index 3rd generation(GIMMS NDVI3g)數據集,數據覆蓋時段為1982—2015年,時間分辨率和空間分辨率分別為15 d和0.083°×0.083°。有效灌溉面積(EIA,1000 hm2)來源于農業部統計局官網。首先將下載的各省份數據除以該省行政區面積得到灌溉密度,再乘以研究流域在各省的分布面積,最后將各省的計算結果求和,即得到研究流域的有效灌溉面積。

2 研究方法

2.1 潛在蒸散發的計算

目前用于估算潛在蒸散發(ET0)的方法主要有溫度法、輻射法、質量傳輸法以及綜合法[17]。其中,Penman-Monteith(PM)法因具有明確物理機制而被聯合國糧農組織(FAO)推薦使用[18]。本文采用考慮多種氣象要素,并經Shuttleworth于1993年修正的PM公式計算ET0[19],具體公式如下：

(1)

式中:λ為潛熱(MJ/kg)；Δ為飽和水汽壓與溫度關系曲線斜率(kPa/℃)；γ為干濕常數(kPa/℃)；Rn為凈輻射(MJ m-2d-1)；G為土壤熱通量(MJ m-2d-1)；u2為2m高處的風速(m/s)；es為空氣飽和水汽壓(kPa)；RH為相對濕度(%)。

2.2 基于Budyko框架計算參數n

Budyko框架是指實際蒸散發(E)受水分供應條件(降水,P)和能量供給條件(潛在蒸散發,ET0)的共同限制,流域P和ET0之間存在耦合平衡關系[20]。其中,經量綱分析和數學推導的傅抱璞公式可較好地反映流域內的水熱耦合狀態,故被廣泛應用[21],其公式為：

(2)

式中：n為流域水熱耦合控制參數,受下墊面變化和氣候變化的影響,決定著流域的水熱分配。

公式(2)中E可通過水量平衡公式計算,其表達式為：

E=P-R-ΔS

(3)

式中：E為多年平均實際蒸散發量(mm)；R為多年平均徑流量(mm)；ΔS為流域蓄水變化量(在多年時間尺度上,閉合流域近似為0[12,20],故可選取5 a作為滑動窗口使其為0來計算E)。

選取5 a為滑動窗口,將滑動窗口內的P、E、ET0值作為輸入數據,即可根據最小二乘法擬合出一個n值來代表滑動窗口中間年份的水熱耦合控制參數值[3,20]。以此類推,將所有窗口的n值全部擬合,即可得到流域時變水熱耦合控制參數序列。

2.3 Mann-Kendall趨勢檢驗法

Mann-Kendall趨勢檢驗簡稱MK檢驗,是一種非參數檢驗方法。由于此方法適用性廣、不受異常值干擾,因而在水文、氣象以及農業等領域被廣泛應用。本研究采用MK法對水熱耦合控制參數n、氣候要素和下墊面因子進行趨勢分析。MK法通過統計量Z來確定趨勢,當顯著性水平為0.05時,若|Z|>1.96,被觀測序列具有顯著性變化趨勢,Z值的正負表示觀測序列呈上升或下降趨勢,具體計算公式見參考文獻[22-23]。

2.4 逐步多元線性回歸模型(SMLR)

為確定氣候和下墊面因子對時變流域水熱耦合控制參數的影響,本研究采用逐步多元回歸模型(Stepwise multiple linear regression model, SMLR)建立時變參數n與氣候和下墊面因子的響應關系。首先建立一個包含所有可能的因變量的模型,然后將其逐步剔除,在保持參數顯著性的同時保持決定系數最高的模型,即最后建立“最優”方程的回歸分析[24]。該回歸方法有效避免了自變量的多重共線性,使建立的模型更加可靠。

選取氣候要素(P、T、Tmax、Tmin、RH、S、u2)和下墊面因子(NDVI和EIA)作為自變量,時變參數n作為因變量,進行逐步回歸擬合方程：

n=α0+(α1xc1+…+αkxck)+(b1xh1+…+bkxhk)

(4)

式中：α0為逐步回歸方程的截距；α1,…,αk及b1,…,bk為自變量的模型系數；xc1,…,xck代表氣候因子,xh1,…,xhk代表下墊面因子。最后將公式(4)計算得到的n值作為模擬值,將最小二乘法確定的n值作為實際值,并基于模擬值和實際值的確定性系數(R2)來確定SMLR模型的擬合效果。

2.5 敏感性分析法

不同氣候區,時變水熱耦合控制參數n對氣候因子和下墊面變化的響應不同。本文選用敏感性分析法來估算干旱與濕潤區流域氣候和下墊面變化對時變參數n的影響。假設時變參數n對某影響因子xi的敏感性系數為Sxi,則Sxi為n對該影響因子的偏導除以n與該影響因子的比值,其計算公式為：

(5)

2.6 貢獻率計算

為定量評估各因子對參數n變化的影響,可將時變參數n對影響因子的敏感性系數與研究時段內該因子的相對變化率相乘,即可得到該影響因子對n變化的貢獻率[25- 26]。具體計算公式如下所示：

Cxi=Sxi×Rxi

(6)

(7)

所有影響因子對時變參數n的貢獻率之和則為氣候因子和下墊面因子對參數n變化的總貢獻率。

3 結果分析

3.1 影響因子的演變規律分析

繪制研究區氣候和下墊面因子在1970—2013年間的年際變化,如圖2所示。對于干旱區的無定河流域,徑流R呈下降趨勢,氣候因子中P、E0、RH、S與u2在1970—2013年間均呈現下降趨勢,下墊面因子NDVI和EIA均呈上升趨勢,從圖2中可以明顯看出,無定河降水呈下降趨勢,而近50年來氣溫呈上升趨勢,這使得流域暖干化現象愈發加劇,而無定河流域既是氣候變化的敏感區,又是環境脆弱區,暖干化這一現象對該地區的水資源可持續利用與工農業發展造成巨大威脅[27]。對于濕潤的漢江上游,徑流R呈下降趨勢,氣候因子中P、E0、S與u2在1970—2013年間均呈現下降趨勢,下墊面因子NDVI和EIA均呈上升趨勢,EIA從1975年開始增長逐漸變緩。

圖2 研究區氣候和下墊面因子變化趨勢示意圖Fig.2 Changes in climatic and underlying surface variables in the study area

采用MK法對無定河流域和漢江上游流域代表氣候和下墊面因子的時間序列進行趨勢檢驗,各因子的趨勢檢驗結果如表1所示。

由表1結果可知,對于干旱區的無定河流域,R、E0、u2通過了顯著性水平為0.05的趨勢檢驗,呈現顯著下降趨勢,P、RH、S呈不顯著的下降趨勢,T、Tmax、Tmin呈現顯著的上升趨勢,這與周園園等[28- 29]所得出的結論一致；下墊面因子NDVI和EIA在1970—2013年呈明顯的上升趨勢,這與Hao等[30]得出的結果一致。對于濕潤區的漢江上游,R和u2呈現顯著下降趨勢,P、E0和S呈不顯著的下降趨勢,T、Tmax、Tmin均呈明顯的上升趨勢,這與李紫妍[32]所得出的結論一致。下墊面因子中NDVI和EIA均呈顯著的上升趨勢。EIA從1975年開始增長逐漸變緩。總體而言,兩流域的氣溫均呈顯著上升趨勢,降水呈下降趨勢,說明所研究的干與濕流域都存在暖干化現象。

表1 MK統計量Z值結果表

3.2 時變水熱耦合參數的趨勢分析

圖3 水熱平衡參數的動態變化Fig.3 The dynamic change of hydrothermal equilibrium parameter

基于水量平衡方程和傅抱璞公式,以11a為滑動窗口來估算時變水熱耦合參數n,估算結果如圖3所示。由圖3可知,干旱區無定河流域與濕潤區漢江上游參數n均呈現上升趨勢。采用MK法對時變參數n進行趨勢檢驗,無定河和漢江上游的Z統計量分別為4.98和4.24,表明兩流域水熱耦合參數n的上升趨勢均顯著。

干旱區無定河流域動態參數n的變化范圍為2.45—2.82其值符合孫福寶等[31]所給的參數范圍,且從1996年開始n值持續上升,最大值出現在2008年,表明在干燥指數(ET0/P)保持不變時,無定河流域自1996年以來蒸發率持續增加,即在同一降水條件下,隨著無定河流域水熱耦合狀況不斷變化,流域的蒸發量持續增大,而下墊面因子中NDVI在1982—1999年呈波動的緩慢增加趨勢,在1999年以后NDVI呈顯著的增加趨勢,說明1999年前該流域生態工程建設相對緩慢,1999年開始施行大規模退耕還林還草和水土保持等生態修復工程,同時,隨著灌溉面積的增加,共同驅動著徑流量顯著減小,流域干旱程度持續加劇。

濕潤區漢江上游動態參數n的變化范圍為1.28—1.65,這與Li等[3]的結果一致,其中1983—1998年參數n持續增大,自1999年以后,參數值明顯減小。表明在相同降水條件下,漢江上游區域1983—1998年間蒸發量增大,流域干旱程度增強。自1999年以后,流域蒸散發有減小的趨勢,流域干旱程度有所緩解,主要是由于1999年以后漢江上游地區EIA呈顯著下降趨勢,徑流有增加的趨勢。

3.3 水熱耦合參數變化的歸因分析

3.3.1時變參數n對影響因子的響應關系

為進一步研究參數n對各影響因子的響應,采用SMLR模型來擬合參數n和氣候因子(P、T、Tmax、Tmin、S、u2、RH)及下墊面因子(NDVI、EIA)的關系式。在逐步回歸的過程中,已剔除與參數n變化相關性弱的影響因子,保留了與時變參數密切相關的因子。模擬結果如下：

無定河：

n1=2.238+0.0476P+0.1127T+0.0577u2+1.8212NDVI

(8)

漢江上游：

n2=2.266-0.0403P+0.0474T-1.218NDVI+0.0351EIA

(9)

由式(8)可知,對干旱區無定河流域來說,P、T、u2和NDVI對水熱耦合參數n的變化具有顯著意義；由式(9)可知,對濕潤區漢江上游來說,P、T、NDVI和EIA與水熱耦合參數n的變化最為相關。表明氣候因素和下墊面因素共同影響著無定河流域和漢江上游水熱耦合參數。此外,將干旱區和濕潤區影響因子的實測值代入公式(8)和(9),計算參數n的模擬值,并將其和參數n的實測值進行對比,結果如圖3所示。SMLR模型對干旱區和濕潤區流域參數n擬合的確定性系數R2分別為0.85和0.98,說明SMLR模型的模擬效果較好。

圖4 參數n實際值與模擬值的比較Fig.4 The comparison between the actual value of parameter n and the simulated value of parameter n0

3.3.2參數n對影響因子的敏感性分析

采用敏感性系數法進一步分析參數n對影響因子的敏感程度,結果如表2所示。

表2 驅動因子對參數n變化的敏感性系數

對干旱區無定河流域,參數n對P的敏感系數最大,為0.485,對T和NDVI的敏感系數次之,對u2的敏感系數最小。表明在各影響因子的變化量相同時,無定河流域水熱耦合參數n對降水變化的響應最為強烈,對氣溫和NDVI次之,對風速的響應最弱,且風速對參數n的變化有負向影響。因為作為水量限制型區域(E0/P>1),無定河降水的變化會直接影響該區域的水熱狀況,而溫度的變化會直接影響潛在蒸散發,進而間接影響流域的水熱狀況。

對濕潤區漢江上游,參數n對NDVI和T的敏感系數最大,對EIA和P的敏感系數次之。表明漢江上游水熱耦合參數對植被覆蓋度和日平均氣溫變化的響應較為強烈,對有效灌溉面積和降水變化的響應次之。其中,植被覆蓋度的改變對參數n的變化有負向影響。總體來說,干旱區流域與濕潤區流域的時變水熱耦合參數對氣候和下墊面變化的敏感程度不同。

3.3.3影響因子對參數n變化的貢獻率

為了進一步定量評估氣候變化和下墊面條件的改變對參數n變化的影響,基于研究區域SMLR模型的結果,計算氣候因子和下墊面因子對水熱耦合參數變化的貢獻率,結果如表3和圖4所示。

表3 驅動因子對參數n變化的貢獻率

圖5 氣候和下墊面因子對參數n變化的貢獻率 Fig.5 Contributions of climatic and underlying surface characteristic factors to parameter n changes

在干旱區無定河流域,P對水熱耦合參數n的貢獻率為-54.7%,T、u2和NDVI對參數n的貢獻率分別為29.5%、35.7%和89.5%,表明P對水熱耦合參數的增加具有負向驅動作用,而T、u2和NDVI的改變對參數的增加具有正向驅動作用。由于降水對參數增加的負向驅動作用小于氣溫、風速和植被覆蓋度的正向驅動作用,故參數呈現增加趨勢。其中,NDVI對參數增加的貢獻率最大,表明植被覆蓋是引起干旱區無定河流域水熱耦合參數增加的主導因子。總體來說,氣候因素對干旱區無定河流域參數n變化的貢獻率為10.5%,而下墊面條件的貢獻率為89.5%,故下墊面條件的改變主導著干旱區無定河流域水熱耦合參數的變化。自1999年黃土高原施行退耕還林(草)生態修復措施以來,黃河流域逐年變綠[32]。無定河位于黃土高原西北部地區,同樣存在變綠的趨勢。無定河作為水量限制型區域(E0/P>1),大規模的植樹造林會通過蒸騰消耗大量的水資源,會對流域水文循環過程產生強烈干擾作用,進一步改變流域水熱狀況,導致水熱耦合參數n的變化,從而使流域徑流減少。因此,植被覆蓋度的增加是無定河流域水熱狀況變化的主導因素。

在濕潤區漢江上游,P、T、NDVI和EIA對參數n變化的貢獻率分別為31.3%、14.4%、-28.8%和83.1%,其中,P、T和EIA的變化對參數n的增加具有促進作用,而NDVI的變化對參數n的增加具有抑制作用。由于降水、氣溫和有效灌溉面積對參數增加的正向驅動作用大于植被覆蓋度的負向驅動作用,故參數在1970—2013年間呈現增加趨勢。從貢獻率的數值大小來看,有效灌溉面積是濕潤區流域參數n變化的主導因子。總體來看,氣候因素對濕潤區水熱耦合參數n變化的貢獻率為45.7%,而下墊面因子的貢獻率為54.3%。因為漢江上游是能量限制型區域(E0/P<1),該區域的植被變綠并未消耗大量的水資源,故植被覆蓋的增加并非導致漢江上游流域水熱耦合參數變化的主要因素。此外,漢江上游流域為陜南的“糧倉”,是重要的農業基地,有效灌溉面積較大[3]。經貢獻率計算結果可知,有效灌溉面積主導著該地區水熱耦合參數的變化。

綜上所述,不同氣候帶的異同主要體現在干燥度指數(E0/P),干旱區無定河流域和濕潤區漢江上游水熱狀況不同,故主導其水熱耦合參數變化的因子也不同。

4 結論

本文選取干旱區無定河流域與濕潤區漢江上游為研究區域,基于Budyko框架,采用MK趨勢檢驗、SMLR模型、敏感性系數法等,對不同水熱條件下,水熱耦合參數演變規律及影響因素進行了分析,并定量評估了干旱區流域與濕潤區流域氣候與下墊面條件的改變對水熱耦合參數變化的貢獻程度,主要結論如下：

(1)在1970—2013年間,兩流域水熱耦合參數整體呈顯著上升趨勢,蒸發量增加而徑流量下降,流域干旱程度總體呈加劇狀態。而對漢江上游而言,自1999年以后,有效灌溉面積顯著減少,蒸發量減少而徑流量增加,其水熱耦合參數出現下降趨勢。

(2)影響干旱流域與濕潤流域水熱平衡參數n的敏感性因子不同。無定河參數n對降水變化具有較高的敏感性,對風速的敏感性最弱；漢江上游參數n對NDVI和氣溫變化均較為敏感。

(3)對于水量限制型流域(無定河),植樹造林通過蒸騰消耗大量的水資源,導致降水在蒸發、下滲及徑流之間的分配發生改變。故植被覆蓋度的增加主導著無定河流域水熱耦合狀況的變化,其貢獻率為89.5%；對于能量限制型流域(漢江上游),有效灌溉面積的增加主導著該流域水熱耦合狀況的變化,其貢獻率為83.1%。

本研究對比分析了不同氣候區氣候和下墊面條件變化對流域水熱耦合狀況的影響,研究結果有助于揭示

變化環境下流域水文循環的演變規律。此外,本研究采用的流域水熱耦合平衡動態變化的驅動力分析框架可以推廣至其他流域。由于數據的限制,本研究在人類活動因子的選取上還不全面,由于下墊面變化對水熱平衡參數變化的影響較為復雜,接下來的研究中會進一步考慮土壤濕度、基流以及遙相關因子(太陽黑子、厄爾尼諾南方濤動等)對流域水熱耦合狀況的影響。