陳澤峰，王衛(wèi)光，李長(zhǎng)妮，丁一民，傅健宇

(1.河海大學(xué) 水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210098；2.河海大學(xué)水文水資源學(xué)院，南京 210098；3.陜西省河流工程技術(shù)研究中心，西安 710018)

近年來(lái)，伴隨著全球變暖和人類活動(dòng)的影響，降水、陸面蒸散發(fā)、流域下墊面、大氣化學(xué)組分等均發(fā)生顯著變化[1,2]，進(jìn)而影響水文循環(huán)過(guò)程[3]，使得流域徑流過(guò)程發(fā)生改變[4]、極端水文事件頻繁發(fā)生，也加劇了我國(guó)水資源短缺、供需矛盾突出、水環(huán)境惡化等的水資源狀況。因此，精確的評(píng)估陸面可利用水資源量變得至關(guān)重要。

Budyko水熱耦合平衡方程[5,6]由于物理機(jī)制明確、輸入資料少，因此在大時(shí)間尺度上可利用水資源量的評(píng)估和預(yù)測(cè)工作中得到了廣泛的應(yīng)用[7,8]。Budyko理論最初是表示多年平均蒸發(fā)率(實(shí)際蒸散發(fā)AET/降水P)和干旱指數(shù)(潛在蒸散發(fā)PET/降水P)關(guān)系的無(wú)參數(shù)方程。隨后，為了刻畫不同下墊面特征的流域Budyko曲線形狀，含有參數(shù)ω的Budyko方程被逐漸發(fā)展起來(lái)。含參數(shù)的Budyko方程對(duì)Budyko曲線的擬合精度較高，也因此應(yīng)用較為廣泛。雖然，眾多學(xué)者研究過(guò)參數(shù)ω與單個(gè)或多個(gè)流域要素間的關(guān)系[7,9]，但由于參數(shù)ω與觀測(cè)變量間的關(guān)系復(fù)雜，對(duì)于不同氣候狀況、不同流域的ω和各因素間關(guān)系的系統(tǒng)變化仍然沒有一個(gè)確定的結(jié)論。對(duì)此，Greve等[10]用美國(guó)流域的觀測(cè)數(shù)據(jù)，建立了從理論上評(píng)估不同流域特征可用水資源量的概率Budyko方程。隨后，基于上述理論，Wang和Fu[11]結(jié)合流域下墊面特征，改進(jìn)了概率Budyko方程，提出了雙變量概率性Budyko框架，提高了可利用水資源量的準(zhǔn)確性。概率Budyko的發(fā)展，為評(píng)估相同氣候條件下，不同下墊面特征的流域可用水資源及其不確定性提供了新的視角。將基于下墊面因素的參數(shù)ω的概率分布代入Budyko方程，可以定量評(píng)估未來(lái)水資源預(yù)估的不確定性，從而為水資源管理者和決策者提供更精確的依據(jù)。

本文結(jié)合IPCC 第五次報(bào)告 CMIP5 的氣候模式數(shù)據(jù)(CanESM2)，運(yùn)用統(tǒng)計(jì)降尺度方法，預(yù)估我國(guó)不同干濕特征的4個(gè)典型流域(黑河流域、灤河山區(qū)、漢江丹江口以上流域和贛江棟背以上流域)未來(lái)水文氣象變量，分析其變化特征；同時(shí)，基于構(gòu)建的雙變量概率性Budyko方程，預(yù)估典型流域不同排放情景下未來(lái)可用水資源比率的分布，并評(píng)估其不確定性。

1 流域概況及資料

本文選取4個(gè)典型流域(黑河流域、灤河山區(qū)、漢江丹江口以上流域和贛江棟背以上流域)為研究區(qū)域。4個(gè)流域具有不同的干濕特征：黑河流域代表干旱流域、灤河山區(qū)代表半干旱流域、漢江丹江口以上流域代表半濕潤(rùn)流域、贛江棟背以上流域代表濕潤(rùn)流域。研究區(qū)位置及氣象站點(diǎn)分布情況如圖1所示。

圖1 4個(gè)典型流域及其氣象站點(diǎn)的分布Fig.1 Spatial distribution of four typical catchments and locations of meteorological stations

研究區(qū)歷史觀測(cè)數(shù)據(jù)包括從國(guó)家氣象局(http:∥cdc.cma.gov.cn/)搜集的4個(gè)流域內(nèi)共21個(gè)氣象站點(diǎn)1961-2000年間逐日降水、最低氣溫、最高氣溫、日照時(shí)數(shù)、平均相對(duì)濕度和平均風(fēng)速資料。潛在蒸散發(fā)通過(guò)具有明確物理機(jī)制的Penman-Monteith方法[12]計(jì)算。各流域上降水、潛在蒸散發(fā)通過(guò)泰森多邊形法[13]計(jì)算流域內(nèi)部氣象站點(diǎn)的面積權(quán)重來(lái)得到。

未來(lái)數(shù)據(jù)本文選擇IPCC 第五次報(bào)告 CMIP5 的CanESM2(the second generation of Canadian Earth System Model)氣候模式網(wǎng)格數(shù)據(jù)，CanESM2氣候模式由加拿大環(huán)境氣候變化部門的氣候模型和分析中心提供，其空間分辨率為2.812 5°×2.812 5°，時(shí)間序列為2006-2100年。CanESM2氣候模式包括3個(gè)排放情景，分別是RCP 26、RCP 45和RCP 85，依次描述的是低排放情景，中等排放情景和高排放情景。

再分析資料(NCEP，National Centers for Environmental prediction)同樣來(lái)自加拿大環(huán)境氣候變化部門的氣候模型和分析中心。該NCEP數(shù)據(jù)集的時(shí)間長(zhǎng)度是1961-2005年，空間尺度與CanESM2氣候模式相同(2.812 5°×2.812 5°)。NCEP數(shù)據(jù)包括26個(gè)大氣因子，主要用來(lái)率定和驗(yàn)證統(tǒng)計(jì)降尺度模型，以預(yù)測(cè)未來(lái)不同排放情景下的氣候變量。

2 研究方法

2.1 概率Budyko方程

1981年，傅抱璞[14]在Budyko假設(shè)的基礎(chǔ)上，根據(jù)流域水文氣候的物理機(jī)制，經(jīng)過(guò)量綱分析和數(shù)學(xué)推導(dǎo)，得到了Budyko假設(shè)的首個(gè)解析表達(dá)式：

(1)

式中：ω為表征流域自然地理特性的無(wú)量綱參數(shù)，且ω≥1。

由于傅抱璞公式的推導(dǎo)過(guò)程為Budyko假設(shè)提供了數(shù)理基礎(chǔ)，具有物理和數(shù)學(xué)意義，所以與其他經(jīng)驗(yàn)的Budyko方程相比，傅抱璞公式的認(rèn)可程度較高并且應(yīng)用較為廣泛。

相比于傳統(tǒng)的Budyko方程給定一個(gè)干旱指數(shù)值只能得到一個(gè)蒸發(fā)率，概率Budyko方程通過(guò)建立參數(shù)ω的概率分布得到的是蒸發(fā)率的一個(gè)分布，繼而可以定量評(píng)估相同氣候條件下，不同下墊面特征流域可利用水資源量的不確定性。

本文以傅抱璞的Budyko水熱耦合平衡方程為依據(jù)，通過(guò)Copula函數(shù)，構(gòu)造Budyko方程參數(shù)ω與流域下墊面因素(植被覆蓋指數(shù)、土壤可用含水量、土壤容重、流域高程、坡度、復(fù)合地形指數(shù))的邊緣分布函數(shù)。經(jīng)過(guò)參數(shù)計(jì)算、擬合優(yōu)度評(píng)價(jià)等選取最優(yōu)的Copula函數(shù)，分別建立ω和各因素的聯(lián)合分布函數(shù)，以植被覆蓋指數(shù)(NDVI)為例，聯(lián)合概率分布見圖2。隨后，計(jì)算ω在流域下墊面條件下的條件分布函數(shù)，反推出在下墊面各因素限制后的參數(shù)ω，進(jìn)而構(gòu)建基于流域下墊面特征的概率Budyko方程。具體方法見Wang和Fu[11]。

圖2 參數(shù)ω與植被覆蓋指數(shù)的聯(lián)合概率分布圖Fig.2 Joint probability distribution of ω and NDVI

2.2 統(tǒng)計(jì)降尺度模型SDSM

統(tǒng)計(jì)降尺度(Statistical Downscaling Model)的核心是建立大尺度氣候因子與站點(diǎn)觀測(cè)變量之間的關(guān)系。由于運(yùn)行速度快且成本低，統(tǒng)計(jì)降尺度相比動(dòng)力降尺度成為眾多學(xué)者首選的方法[15]，也因此被廣泛應(yīng)用于評(píng)估氣候變化背景下氣象、水文、環(huán)境等的響應(yīng)[16,17]。

統(tǒng)計(jì)降尺度的主要過(guò)程可分為以下4步：預(yù)測(cè)因子的篩選；構(gòu)建模型；模型的率定和驗(yàn)證；未來(lái)氣象變量的預(yù)測(cè)。

3 結(jié)果與分析

3.1 統(tǒng)計(jì)降尺度模型的率定與驗(yàn)證

預(yù)測(cè)因子的篩選是統(tǒng)計(jì)降尺度的核心環(huán)節(jié)，主要依據(jù)觀測(cè)變量和多個(gè)再分析大氣因子間的相關(guān)性和偏相關(guān)性，從而篩選出最相關(guān)的預(yù)測(cè)因子。本文計(jì)算了4個(gè)典型流域日尺度的潛在蒸散發(fā)與降水和26個(gè)氣候因子間的相關(guān)系數(shù)(corr)和偏相關(guān)系數(shù)(p-corr)，并從中選取最相關(guān)的預(yù)報(bào)因子。

因子篩選之后，進(jìn)行模型計(jì)算和率定驗(yàn)證。本文以1961-1990年作為模型率定期，1991-2000年為模型驗(yàn)證期，用相關(guān)系數(shù)R、模擬值標(biāo)準(zhǔn)差與觀測(cè)值標(biāo)準(zhǔn)差的比值RS、確定性系數(shù)Ens和偏差Bias4個(gè)指標(biāo)來(lái)評(píng)價(jià)模型，以灤河山區(qū)的6個(gè)氣象站點(diǎn)為例，評(píng)價(jià)結(jié)果見表1。

可以看出，模型對(duì)潛在蒸散發(fā)的模擬效果較好，R和RS在率定期和驗(yàn)證期均超過(guò)了0.92，多倫站和青龍站的RS甚至等于1；確定性系數(shù)Ens也都大于0.80；偏差Bias大都接近于0。降水的模擬結(jié)果雖然不如潛在蒸散發(fā)好，但R都在0.43以上，RS都大于0.57，也能用來(lái)預(yù)估未來(lái)的氣象變量。

流域尺度上，比較了驗(yàn)證期(1991-2000年)實(shí)測(cè)月潛在蒸散發(fā)和模擬月潛在蒸散發(fā)。黑河流域、灤河山區(qū)、漢江丹江口以上流域和贛江棟背以上流域的R2分別達(dá)到0.993、0.992、0.977和0.975，說(shuō)明整體模擬效果好。同時(shí)比較了驗(yàn)證期(1991-2000年)實(shí)測(cè)月降水和模擬月降水。黑河流域、灤河山區(qū)、漢江丹江口以上流域和贛江棟背以上流域的R2分別達(dá)到0.800、0.912、0.811和0.812，4個(gè)流域的R2都大于0.80，說(shuō)明模擬效果較好。總之，SDSM在構(gòu)建降水和潛在蒸散發(fā)時(shí)具有良好的模擬效果，可用于對(duì)未來(lái)數(shù)據(jù)進(jìn)行降尺度。

3.2 流域未來(lái)水文氣象變量的預(yù)估

根據(jù)建立好的SDSM，輸入不同排放情景的預(yù)報(bào)因子，便可得到未來(lái)不同情景下的日潛在蒸散發(fā)量和日降水量，進(jìn)而累加計(jì)算得到年潛在蒸散發(fā)量和降水量，其結(jié)果如圖3所示。4個(gè)流域的潛在蒸散發(fā)量在未來(lái)RCP 26、RCP 45和RCP 85情景下均呈上升趨勢(shì)，且RCP 85情景下的增加趨勢(shì)更為明顯，尤其在2050年之后。降水的整體增加趨勢(shì)相對(duì)較小，尤其在濕潤(rùn)的贛江棟背以上流域。

4個(gè)流域未來(lái)時(shí)段2021-2050年3個(gè)排放情景的多年平均降水量(P)和潛在蒸散發(fā)量(PET)，繼而計(jì)算得到相應(yīng)的干旱指數(shù)(AI)。結(jié)果表明，在RCP26情景下，黑河流域、灤河山區(qū)、漢江丹江口以上流域以及贛江棟背以上流域的AI分別為8.944、1.651、1.144和0.884；在RCP45情景下，黑河流域、灤河山區(qū)、漢江丹江口以上流域以及贛江棟背以上流域的AI分別為8.797、1.727、1.145和0.867；在RCP85情景下，黑河流域、灤河山區(qū)、漢江丹江口以上流域以及贛江棟背以上流域的AI分別為8.745、1.710、1.131和0.879。

表1 統(tǒng)計(jì)降尺度模型在構(gòu)建灤河山區(qū)各氣象站點(diǎn)日潛在蒸散發(fā)和日降水時(shí)的模擬結(jié)果評(píng)價(jià)Tab.1 Evaluation of daily potential evapotranspiration and precipitation at each at each meteorological station of Luanhe Mountain that simulated by statistical downscaling model

圖3 流域歷史實(shí)測(cè)潛在蒸散發(fā)、降水和CanESM2模式不同氣候情景下的潛在蒸散發(fā)、降水Fig.3 Measured potential evapotranspiration and the potential evapotranspiration under the different emission scenarios based on CanESM2 for four catchments

3.3 流域2021-2050年可用水資源比率的預(yù)估

在構(gòu)建概率Budyko方程的過(guò)程中，通過(guò)對(duì)比下墊面各因素限制后的參數(shù)ω′和ω的概率密度函數(shù)，得到了各下墊面因素對(duì)參數(shù)ω的影響程度。結(jié)果表明，黑河流域和贛江棟背以上流域的最大影響因素為植被覆蓋指數(shù)，漢江丹江口以上流域的最大影響因素為土壤可用含水量，而灤河山區(qū)的最大影響因素為高程。因此，根據(jù)建立的概率Budyko方程，假定流域2021-2050年參數(shù)ω分布與歷史時(shí)期(1982-2008年)ω的分布保持一致，即黑河流域和贛江棟背以上流域的參數(shù)是植被覆蓋指數(shù)限定后的ω′，灤河山區(qū)和漢江丹江口以上流域的參數(shù)分別是由高程和土壤可用含水量限定后的ω′。將4個(gè)典型流域?qū)?yīng)的ω′及得到的未來(lái)不同排放情景下干旱指數(shù)(AI)分別代入Budyko方程，便可得到未來(lái)可用水資源比率的分布，從而評(píng)估其不確定性。

圖4展示了4個(gè)流域不同排放情景下2021-2050年多年平均可用水資源比率的分布。可發(fā)現(xiàn)，黑河流域可用水資源比率值分布范圍較窄，多在0～0.06之間，不同排放情景的可用水資源比率值較為接近[圖4(a)]。灤河山區(qū)可用水資源比率范圍相對(duì)較寬，多集中在0～0.5之間，且有雙峰現(xiàn)象，說(shuō)明有兩段值出現(xiàn)的頻率較大，0.1～0.25、0.25～0.4為兩個(gè)雙峰的大概位置；RCP 26 情景與另外兩個(gè)排放情景下的可用水資源比率略有差異[圖4(b)]。漢江丹江口以上流域可用水資源比率大都集中在0.2～0.6間，雙峰現(xiàn)象略不明顯，不同排放情景的可用水資源比率值較為接近[圖4(c)]。贛江棟背以上流域可用水資源比率多分布在0.3～0.8的范圍內(nèi)，雙峰現(xiàn)象明顯，且在峰值以前的曲線斜率相對(duì)峰值以后的曲線斜率較小，不同排放情景的可用水資源比率值差異較小[圖4(d)]。

圖4 流域2021-2050年不同排放情景下可用水資源比率的分布Fig.4 Water availability ratio under the different emission scenarios during 2021-2050 for four catchments

從可用水資源比率的箱線圖(圖5)可看出，各流域可用水資源比率在不同排放情景下的分布極為相似，表明不同的排放情景對(duì)可用水資源比率的影響不大；與常理一致，可用水資源比率值隨著流域濕潤(rùn)程度的增加而增加，黑河流域、灤河山區(qū)、漢江丹江口以上流域和贛江棟背以上流域可用水資源比率的中位數(shù)依次為0.10、0.20、0.38和0.57；箱子的高度(即四分位距IQR)在黑河流域最矮(0.022)，漢江丹江口以上流域次之(0.144)，接著是贛江棟背以上流域(0.151)，灤河山區(qū)最高(0.194)，說(shuō)明可用水資源比率的不確定性在灤河山區(qū)最大，黑河流域最小，其余兩個(gè)流域介于之間。此外，黑河流域和贛江棟背以上流域的下墊面最大影響因素同為植被覆蓋指數(shù)，兩個(gè)流域的Budyko參數(shù)均為由植被覆蓋指數(shù)限定之后的參數(shù)ω，但其可用水資源比率的大小及不確定性卻相差甚大，說(shuō)明流域氣候特征對(duì)流域水資源狀況影響強(qiáng)烈，且加入植被覆蓋指數(shù)后可用水資源比率的不確定性在干旱區(qū)域(黑河流域)比濕潤(rùn)區(qū)域(贛江棟背以上流域)小。

箱子最底端和最頂端的短線分別表示可用水資源比率的最小值和最大值，箱子從下往上的三條線依次表示25%、50%、75%處的值，實(shí)心圓代表可用水資源比率的平均值圖5 流域2021-2050年不同排放情景下可用水資源比率的箱線圖Fig. 5 Boxplot of water availability ratio under the different emission scenarios during 2021-2050 for four catchments

總之，不同排放情景對(duì)流域可用水資源比率的影響較小，流域可用水資源比率值隨著流域濕潤(rùn)程度的增加而增加，可用水資源比率的不確定性從高到低依次是灤河山區(qū)、贛江棟背以上流域、漢江丹江口以上流域和黑河流域。此外，濕潤(rùn)區(qū)贛江棟背以上流域的可用水資源比率的不確定性略高于半濕潤(rùn)區(qū)漢江丹江口以上流域的這一結(jié)果也反映出研究不同下墊面因素對(duì)參數(shù)ω影響程度的必要性。

4 結(jié) 語(yǔ)

本文運(yùn)用統(tǒng)計(jì)降尺度模型(SDSM)，結(jié)合CanESM2氣候模式、實(shí)測(cè)降水、由Penman-Monteith公式計(jì)算的潛在蒸散發(fā)，預(yù)估未來(lái)3種不同輻射強(qiáng)迫和大氣排放標(biāo)準(zhǔn)(RCP 26、RCP 45、RCP 85)下的降水量和潛在蒸散發(fā)量；基于構(gòu)建的雙變量概率Budyko方程，對(duì)我國(guó)不同干旱特征的典型流域(黑河流域、灤河山區(qū)、漢江丹江口以上流域和贛江棟背以上流域)未來(lái)可用水資源比率進(jìn)行預(yù)估，并分析其不確定性。主要結(jié)果如下：

(1)4個(gè)典型流域的潛在蒸散發(fā)量在未來(lái)RCP 26、RCP 45和RCP 85情景下均呈上升趨勢(shì)，且在2050年之后，RCP 85情景下的潛在蒸散發(fā)量最大，RCP 45情景下的潛在蒸散發(fā)量次之，RCP 26情景下的潛在蒸散發(fā)量最??；降水的整體增加趨勢(shì)則相對(duì)較小，尤其在濕潤(rùn)的贛江棟背以上流域。

(2)黑河流域、灤河山區(qū)、漢江丹江口以上流域和贛江棟背以上流域可用水資源比率的中位數(shù)依次為0.10、0.20、0.38和0.57，即4個(gè)流域2021-2050年多年平均可用水資源比率值隨著濕潤(rùn)程度的增加而增加；不同排放情景下的可用水資源比率分布的差異較小；灤河山區(qū)和贛江棟背以上流域的可用水資源比率分布有雙峰現(xiàn)象。

(3)不同排放情景對(duì)流域可用水資源比率的不確定性影響較?。?個(gè)流域可用水資源比率的不確定性從高到低依次是灤河山區(qū)(IQR為0.194)、贛江棟背以上流域(IQR為0.151)、漢江丹江口以上流域(IQR為0.144)和黑河流域(IQR為0.022)。