石清, 邢萬秋

(1.中冶華天南京工程技術(shù)有限公司，江蘇 南京 210019； 2.河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098； 3.河海大學(xué) 水文水資源學(xué)院，江蘇 南京 210098)

近幾十年來，眾多研究沿著氣候變化的方向?qū)λ难h(huán)過程和水資源進(jìn)行模擬和預(yù)測(cè)[1-5]。其中廣泛應(yīng)用于氣候變化背景下徑流響應(yīng)預(yù)測(cè)的方法主要有兩類：基于水文模型的方法和基于Budyko假設(shè)的水量平衡方法[6]。在年以上的尺度，基于Budyko假設(shè)的水量平衡方法以水量平衡為基礎(chǔ)，在減少不確定性方面具有優(yōu)勢(shì)，可用于預(yù)測(cè)評(píng)估未來徑流的變化，而目前通過Budyko水熱耦合平衡方程的方法對(duì)區(qū)域水資源進(jìn)行預(yù)測(cè)的研究較少。水文模型在預(yù)測(cè)未來水資源量時(shí)均基于歷史時(shí)期率定和驗(yàn)證后的模型參數(shù)，同樣地，基于Budyko水熱耦合平衡方程在對(duì)未來區(qū)域水資源量進(jìn)行預(yù)測(cè)時(shí)，以往的研究多假設(shè)Budyko方程中的參數(shù)在未來氣候情景下恒定不變。實(shí)際上，隨著未來氣候情景的變化、流域特征的改變和人類活動(dòng)的加劇，區(qū)域水熱耦合平衡狀態(tài)也將隨之發(fā)生改變。在氣候條件和水文情勢(shì)發(fā)生改變的同時(shí)，Budyko方程參數(shù)應(yīng)是變化的。因此，綜合考慮氣候變化和Budyko方程參數(shù)的變化進(jìn)行區(qū)域水資源量預(yù)測(cè)極有必要。

為此，本研究選取了代表我國不同氣候區(qū)的3個(gè)典型流域(松花江流域、黃河流域和珠江流域)，通過研究其水文氣象變量的變化趨勢(shì)，計(jì)算分析Budyko參數(shù)的變化規(guī)律，利用基于Budyko方程的動(dòng)態(tài)水量平衡模型實(shí)現(xiàn)3個(gè)流域未來的水資源量預(yù)測(cè)。

1 研究區(qū)概況

選取分別代表我國不同氣候區(qū)的珠江流域(河源站控制的東江流域)、松花江流域(豐滿水庫站控制的松花江上游流域)和黃河源區(qū)(唐乃亥水文站控制流域)作為典型流域，如圖1所示。

圖1 研究選取的典型流域

河源水文站位于廣東省河源市，控制東江秋香江口以上區(qū)域，流域面積為15 750 km2;東江流域?qū)賮啛釒Ъ撅L(fēng)氣候，四季不甚分明，年平均氣溫為22.1 ℃，多年平均降水量約為2 000 mm。豐滿水庫水文站位于吉林省東部山地的西側(cè)，地處松花江上游，控制面積42 500 km2;松花江流域地處北溫帶季風(fēng)氣候區(qū)，大陸性氣候特點(diǎn)顯著，多年平均氣溫為3～5 ℃，年內(nèi)溫差較大，多年平均降水量約為700 mm。黃河源區(qū)地處青藏高原東部邊緣，以唐乃亥水文站斷面為界，流域面積121 972 km2;屬高原山地氣候，海拔較高，多年平均氣溫僅為1 ℃，多年平均降水量約350 mm。

2 數(shù)據(jù)處理與研究方法

2.1 氣候情景數(shù)據(jù)

本研究采用CMIP5公開的4種氣候模式(包括BCCCSM1-1M、GFDL-ESM2M、HadGEM2-ES以及MIROC-ESMCHEM，發(fā)布機(jī)構(gòu)和模式的空間分辨率見表1)，氣候情景采用3種不同的輻射強(qiáng)迫和大氣排放標(biāo)準(zhǔn)，即低排放情景(RCP2.6)、中間穩(wěn)定排放情景(RCP4.5)和高排放情景(RCP8.5)。未來情景數(shù)據(jù)選用降水量、最高氣溫和最低氣溫、潛在蒸發(fā)量4個(gè)變量。

表1 選用的氣候模式發(fā)布機(jī)構(gòu)及其分辨率

2.2 研究方法

2.2.1 基于Budyko方程的動(dòng)態(tài)水量平衡模型

Budyko方程描述了實(shí)際蒸發(fā)依賴于潛在蒸發(fā)量(熱能來源)和可利用水量(水量來源)的程度，利用Budyko框架可以有效評(píng)價(jià)氣候、植被和水文循環(huán)之間的相互影響。其中，經(jīng)量綱分析和數(shù)學(xué)推導(dǎo)并被廣泛應(yīng)用的傅抱璞公式[7]可以較好地表達(dá)流域內(nèi)的水熱耦合狀態(tài)。YANG H等[8]以傅抱璞公式為基礎(chǔ)，將Choudhury公式變換成另一種表達(dá)Budyko曲線的方程，如下：

(1)

式中：E為實(shí)際蒸發(fā)量，mm；E0為潛在蒸發(fā)量，mm；P為降水量，mm；n為參數(shù)，需要率定。

在計(jì)算潛在蒸發(fā)量E0時(shí)，由于具有物理機(jī)理的Penman-Monteith方法對(duì)資料要求較高，容易造成該方法與其它簡(jiǎn)化經(jīng)驗(yàn)公式(如Hargreaves方法)同樣甚至更差的水文循環(huán)模擬效果[9]。為此，本文選擇改進(jìn)的Hargreaves公式，僅需降水量、最高氣溫、最低氣溫等數(shù)據(jù)輸入，具體公式如下：

E0=0.001 3S0(Tmean+17)(TD-0.012 3P)0.76。

(2)

式中：Tmean為某天中最高與最低氣溫均值；TD為最高與最低氣溫差值；S0為太陽能輻射的水當(dāng)量，與測(cè)站所在的地理緯度有關(guān)[10]。

依據(jù)流域水量平衡原理，可建立動(dòng)態(tài)的水量平衡模型：

(3)

式中：Sc(t)為流域總蓄水量；P(t)、E(Sc,t)和R(Sc,t)分別為流域的降水量、蒸發(fā)量和徑流量。多年平均穩(wěn)定狀態(tài)下，Sc(t)=0[11]。

2.2.2 氣候模式數(shù)據(jù)處理方法

1)偏差修正法。在用氣候模式數(shù)據(jù)對(duì)區(qū)域氣候水文預(yù)測(cè)時(shí)，由于氣候模式數(shù)據(jù)分辨率較低以及區(qū)域的不匹配性，首先要對(duì)氣候模式數(shù)據(jù)進(jìn)行降尺度和偏差修正。將氣候模式較粗的網(wǎng)格數(shù)據(jù)降尺度到0.5°×0.5°網(wǎng)格上，然后再求流域平均值。偏差修正法公式如下：

(4)

2)氣候模式數(shù)據(jù)有效性評(píng)價(jià)法。為了較系統(tǒng)地評(píng)價(jià)氣候模式數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)的模擬效果，采用確定性系數(shù)R、效率系數(shù)Ens來分別衡量模型對(duì)量及變化趨勢(shì)的捕捉能力，利用平均變差BIAS來評(píng)價(jià)模型的系統(tǒng)偏差，用標(biāo)準(zhǔn)偏差的比值RS來檢驗(yàn)?zāi)Ｐ蛯?duì)指標(biāo)變率的捕捉能力。各計(jì)算公式如下：

(5)

(6)

(7)

(8)

除上述方法之外，本研究還采用了時(shí)間序列趨勢(shì)檢驗(yàn)法對(duì)降水、徑流、氣溫等水文氣象要素以及Budyko參數(shù)的變化規(guī)律加以分析；采用了滑動(dòng)平均法結(jié)合最小二乘法推求Budyko方程參數(shù)。

3 結(jié)果分析

3.1 水文氣象變量變化趨勢(shì)

對(duì)歷史時(shí)期1956—2000 年間的降水量、潛在蒸發(fā)量、最高和最低氣溫4個(gè)水文氣象變量的趨勢(shì)變化加以分析，結(jié)果如圖2所示。圖2中，P為降水量，E0為潛在蒸發(fā)量，Tmax為最高氣溫，Tmin為最低氣溫。1956—2000年水文氣象變量變化及其顯著性結(jié)果見表2。表2中坡度和p分別為線性回歸的坡度值和顯著性；Z為Mann-Kendall顯著性檢測(cè)值；“*”、“**”和“***”分別代表顯著性檢測(cè)通過90%、95%和99%置信水平。

由圖2可見：黃河流域的降水量、最高氣溫、最低氣溫和潛在蒸發(fā)量均呈增加或升高趨勢(shì)，其中，最高氣溫和最低氣溫上升趨勢(shì)顯著(見表2)；松花江流域上的降水量呈現(xiàn)非顯著性的降低趨勢(shì)，氣溫顯著增加，潛在蒸發(fā)量非顯著上升；珠江流域上的降水量呈增加趨勢(shì)，最低氣溫呈顯著上升趨勢(shì)，而最高氣溫和潛在蒸發(fā)量卻呈下降趨勢(shì)，其中潛在蒸發(fā)量下降趨勢(shì)顯著(通過T參數(shù)檢驗(yàn)90%置信區(qū)間)。

圖2 1956—2000年水文氣象變量變化過程

表2 1956—2000年水文氣象變量變化及其顯著性

3.2 Budyko方程參數(shù)的變化

對(duì)實(shí)測(cè)1956—2000 年降水量、潛在蒸發(fā)量和徑流深用11 a窗口移動(dòng)對(duì)其進(jìn)行平滑和濾波處理，再通過最小二乘法求得每個(gè)滑動(dòng)窗口的Budyko方程參數(shù)，Budyko方程參數(shù)變化趨勢(shì)及其顯著性結(jié)果見表3，參數(shù)結(jié)果如圖3所示。

表3 Budyko方程參數(shù)變化趨勢(shì)及其顯著性

圖3 11 a滑動(dòng)平均后Budyko方程參數(shù)變化

由圖3和表3可看出：除了黃河源參數(shù)滑動(dòng)平均值序列呈顯著增加趨勢(shì)(通過0.01 顯著性檢驗(yàn))外，其他兩個(gè)流域(松花江流域和珠江流域)上的Budyko 方程參數(shù)均呈非顯著性的減小趨勢(shì)；珠江流域上參數(shù)的滑動(dòng)平均值序列在歷史時(shí)期變化幅度較大，而黃河流域和松花江流域上的參數(shù)值序列變化較為平穩(wěn)。

3.3 區(qū)域水資源量預(yù)測(cè)

3.3.1 氣候模式數(shù)據(jù)處理結(jié)果

根據(jù)歷史實(shí)測(cè)降水量、最高氣溫、最低氣溫計(jì)算潛在蒸發(fā)量E0，作為歷史實(shí)測(cè)E0的標(biāo)準(zhǔn)值。氣候模式數(shù)據(jù)經(jīng)降尺度和偏差修正后，歷史時(shí)期的氣候模式數(shù)據(jù)和歷史實(shí)測(cè)月值序列相比較，結(jié)果見表4(以黃河流域?yàn)槔?。表4中：4個(gè)變量中最高氣溫的模擬效果最好，最低氣溫的次之，降水量的模擬效果最差；4個(gè)氣候模式中，以GFDL模式對(duì)黃河流域上的水文氣候變量模擬精度最高，以HadGEM2模式對(duì)黃河流域上的水文氣候變量的模擬精度最低；4種氣候模式均在一定程度上較好地模擬了區(qū)域水文氣候變量，本文采納4種氣候模式數(shù)據(jù)進(jìn)行研究。

表4 歷史的氣候模式數(shù)據(jù)和歷史實(shí)測(cè)月值數(shù)據(jù)的比較

圖4給出了黃河流域上歷史實(shí)測(cè)數(shù)據(jù)和氣候模式RCP2.6情景下數(shù)據(jù)的11 a滑動(dòng)平均值長(zhǎng)期序列。

圖4 氣候模式數(shù)據(jù)和實(shí)測(cè)歷史數(shù)據(jù)長(zhǎng)期序列11 a滑動(dòng)平均值

由圖4可見，歷史時(shí)期4種氣候模式對(duì)潛在蒸發(fā)量的模擬結(jié)果較為一致，同時(shí)未來情景下潛在蒸發(fā)量呈現(xiàn)增加的趨勢(shì)。

3.3.2 2021—2050年多年平均徑流預(yù)測(cè)結(jié)果

采用4種氣候模式、3種氣候情景，基于黃河流域、松花江流域和珠江流域重新構(gòu)建的區(qū)域Budyko水熱耦合平衡方程，計(jì)算未來時(shí)段(2021—2050年)區(qū)域水資源量。3個(gè)流域的降水量、潛在蒸發(fā)量、最高氣溫和最低氣溫的變化如圖5所示。圖5中，在3個(gè)典型流域，最高氣溫和最低氣溫均隨著輻射強(qiáng)迫和大氣排放標(biāo)準(zhǔn)的升高而變幅增大。

圖5 未來時(shí)段(2021—2050年)降水量、潛在蒸發(fā)量、最高氣溫、最低氣溫相對(duì)基準(zhǔn)期(1971—2000年)的變化率或變化

從圖5中可以看出：相較于各氣候模式的歷史時(shí)期均值，黃河流域在未來時(shí)段(2021—2050年)內(nèi)4個(gè)水文氣象變量均增加，同時(shí)MIROC模式的降水量、最高氣溫和最低氣溫在4個(gè)模式中增幅最大，而潛在蒸發(fā)量在HadGEM2模式下的變化率最大；潛在蒸發(fā)量、最高氣溫和最低氣溫隨著輻射強(qiáng)迫和大氣排放標(biāo)準(zhǔn)的升高而變化率或變幅增大，而降水量則相反；松花江流域2021—2050年內(nèi)降水量在GFDL模式下相較于歷史時(shí)期將減少，潛在蒸發(fā)量在GFDL模式下隨著排放情景的升高而變化率變小，在BCCCSM1模式下隨著排放情景的升高而變化率增大；珠江流域在MIROC模式下的降水量隨著排放情景的升高而減小幅度變大，而在HadGEM2模式下剛好相反；相較于黃河流域和松花江流域，珠江流域除MIROC模式外，潛在蒸發(fā)量在未來時(shí)段(2021—2050年)的變化率較小。

圖6顯示出未來時(shí)段(2021—2050年)3個(gè)典型流域上Budyko水熱耦合平衡方程參數(shù)與歷史1971—2000年同時(shí)期參數(shù)的對(duì)比。

由圖6可以看出，黃河流域的水熱耦合平衡方程參數(shù)自1971—2000年間的0.609增加到2021—2050年間的0.771，松花江流域和珠江流域的參數(shù)分別從1.032、1.560減小到0.921、1.494。說明未來時(shí)段(2021—2050年)相較于歷史時(shí)期，黃河流域下墊面改變將變大，流域特征對(duì)流域水循環(huán)過程的影響將加劇；而松花江流域和珠江流域的下墊面改變將減弱，對(duì)流域水循環(huán)過程的影響將減弱。

圖6 2021—2050年與歷史1971—2000年的Budyko方程參數(shù)變化

圖7給出了未來時(shí)段(2021—2050年)內(nèi)3個(gè)典型流域的模擬徑流深相較于歷史實(shí)測(cè)徑流深的變化。

圖7 氣候模式下2021—2050年相較于1971—2000年實(shí)測(cè)平均徑流深的變化率

由圖7可以看出：黃河流域未來時(shí)段(2021—2050年)的徑流均值相較于基準(zhǔn)期的均將減少，3種氣候排放情景下徑流的變化率差別不大，在HadGEM2模式下徑流減少的百分率最大，在RCP 8.5情景下達(dá)到40%；松花江流域在2021—2050年的徑流深多年平均值相較于基準(zhǔn)期基本均將增加，僅RCP2.6情景下的GFDL模式、RCP8.5情景下的MIROC模式預(yù)測(cè)的徑流深減少，減少率均在10%以內(nèi)；珠江流域在2021—2050年預(yù)測(cè)徑流深將減少，僅RCP2.6情景下的BCCCSM1模式、RCP4.5和RCP8.5情景下的GFDL模式的徑流深將增加。

綜合圖5—7得到以下結(jié)論：預(yù)測(cè)的黃河流域未來時(shí)段(2021—2050年)降水量增加而徑流深卻將減少，可能是由于Budyko方程參數(shù)變大導(dǎo)致；松花江流域未來時(shí)段(2021—2050年)的降水量增加、潛在蒸發(fā)量增加、Budyko方程參數(shù)減小，未來徑流深也將增加；珠江流域未來時(shí)段(2021—2050年)的降水量減小，未來徑流深也將減少。

3.3.3 與恒定Budyko方程參數(shù)預(yù)測(cè)結(jié)果的比較

為評(píng)價(jià)Budyko方程參數(shù)對(duì)區(qū)域未來水資源量的影響，假設(shè)未來氣候情景下Budyko方程參數(shù)維持在基準(zhǔn)期1971—2000年的恒定不變，以此參數(shù)預(yù)測(cè)2021—2050年的多年平均徑流深，用變化的Budyko方程參數(shù)預(yù)測(cè)未來徑流深的變化率，結(jié)果見表5。

表5 基于基準(zhǔn)期Budyko方程參數(shù)預(yù)估2021—2050年的徑流深變化率 %

從表5中可以看出:僅氣候發(fā)生改變而流域下墊面不變的情況下，4種氣候模式預(yù)測(cè)的3個(gè)典型流域上的徑流深在未來30 a的均值大多將減少，和基準(zhǔn)期實(shí)測(cè)徑流深相比，BCCCSM1氣候模式下以珠江流域上RCP4.5情景下減少量最大，為-18.6%，該期間降水量的變化率為-9.9%;GFDL模式下預(yù)測(cè)的2021—2050多年平均徑流深相較于其他氣候模式的變化率最小，在松花江流域和珠江流域的中高排放情景下(RCP4.5和RCP8.5)的徑流深微弱增加；HadGEM2模式下預(yù)測(cè)的3個(gè)流域上的徑流深均將減少，僅在松花江流域RCP4.5情景下的徑流深增加10.9%(預(yù)測(cè)的降水量將增加19.9%)；MIROC模式下預(yù)測(cè)的徑流深也都將減少，僅松花江流域上RCP4.5情景下的徑流深將增加2.5%。

結(jié)合圖7和表5可知:盡管未來氣候情景下預(yù)測(cè)的黃河流域徑流深在未來期間(2021—2050年)都減少，但相比基于變化的Budyko方程參數(shù)預(yù)測(cè)的徑流深減少得更多。

實(shí)際上在預(yù)測(cè)區(qū)域未來(2021—2050年)水資源量時(shí)，僅考慮氣候情景的變化是不夠的，本文基于前面的研究對(duì)代表流域下墊面的Budyko方程參數(shù)進(jìn)行了預(yù)測(cè)，基于氣候變化和流域下墊面改變預(yù)測(cè)的未來(2021—2050年)徑流的變化將更準(zhǔn)確;以RCP4.5 情景為例，圖8給出了Budyko 方程參數(shù)變化前、后預(yù)測(cè)的徑流變化率的對(duì)比。

圖8 RCP4.5 情景下Budyko方程參數(shù)變化前后預(yù)估徑流變化率的對(duì)比

由圖8結(jié)合圖5可以看出：雖然黃河流域2021—2050年期間降水預(yù)測(cè)量將增加，徑流深卻將減少(在HadGEM2模式下的減少率最高)，而未來(2021—2050年)流域下墊面的改變更將加劇該區(qū)水資源的緊張(圖8(a))。由圖8(b)可以看出，松花江流域上隨著未來(2021—2050年)降水量的增加，徑流深也將增加(在HadGEM2模式下的增加率最高)，同時(shí)，下墊面的改變更將造成該區(qū)未來水資源量的增加，在對(duì)該區(qū)進(jìn)行未來水資源規(guī)劃時(shí)，應(yīng)防止該區(qū)發(fā)生洪澇災(zāi)害。圖8(c)顯示出，在珠江流域上除了GFDL 模式下預(yù)測(cè)的徑流在2021—2050年期間呈增加趨勢(shì)外，其他3 個(gè)氣候模式預(yù)測(cè)的未來徑流深均將減少(在BCCCSM1模式下的減少率最大)。同時(shí)，從圖8(c)也可看出，流域下墊面對(duì)該區(qū)水資源的影響不大。

4 結(jié)語

本研究構(gòu)建了基于Budyko方程的動(dòng)態(tài)水量平衡模型，選用CMIP5中的4種氣候模式(BCCCSM1-1M、GFDL-ESM2M、HadGEM2-ES 和MIROC-ESMCHEM)、3種不同的輻射強(qiáng)迫和大氣排放標(biāo)準(zhǔn)(RCP2.6、RCP4.5和RCP8.5)，對(duì)我國3個(gè)代表不同氣候區(qū)的典型流域(黃河流域、松花江流域和珠江流域)的未來徑流進(jìn)行預(yù)測(cè)，得到結(jié)論如下：

1)3個(gè)流域的最低氣溫均呈顯著上升趨勢(shì)，黃河流域與松花江流域最高氣溫呈顯著上升趨勢(shì)，珠江流域的潛在蒸發(fā)量下降趨勢(shì)顯著(通過T參數(shù)檢驗(yàn)90%置信區(qū)間)；除了黃河流域的滑動(dòng)平均值序列呈顯著增加趨勢(shì)外，松花江流域和珠江流域的區(qū)域Budyko方程參數(shù)呈現(xiàn)非顯著性的減小趨勢(shì)。

2)預(yù)測(cè)黃河流域2021—2050年間降水量、潛在蒸發(fā)量增加，Budyko方程參數(shù)增大，徑流深減?。凰苫ń饔蛭磥斫邓?、潛在蒸發(fā)量增加，Budyko方程參數(shù)減小，徑流深也增加；珠江流域未來在RCP4.5和RCP8.5情景下僅在GFDL模式下的徑流深增加，其余3個(gè)氣候模式預(yù)測(cè)的徑流深減少。

3)與基于不變的Budyko模型參數(shù)預(yù)測(cè)的徑流深相比，黃河流域由于未來下墊面改變會(huì)加劇該區(qū)水資源量減少，松花江流域下墊面的改變更將造成該區(qū)未來水資源量的增加，珠江流域上下墊面改變對(duì)Budyko模型預(yù)測(cè)該區(qū)2021—2050年間的徑流深改變不明顯。