周祖昊，劉佳嘉，嚴子奇，王 浩，賈仰文

(1.中國水利水電科學研究院流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038；2.水利部水資源與水生態(tài)工程技術(shù)研究中心,北京 100038)

隨著全球氣候變暖和人類活動影響加劇，流域水資源演變發(fā)生重大變化，這是世界水文水資源研究領(lǐng)域的熱門課題[1-3]。國內(nèi)外學者大多集中于討論氣候和人類活動對實測徑流量演變過程的影響，對天然河川徑流量變化歸因分析相對較少。天然河川徑流量決定了流域可用水資源量以及相關(guān)工程布局，厘清天然河川徑流量衰減原因，有助于為流域水資源高效利用和科學管理提供支撐。

黃河是中國北方重要的淡水資源，河川徑流量約占全國2%，黃河流域承載了全國15%的人口、15%的耕地面積和14%的GDP[4]，是中國水資源短缺非常嚴重的流域，徑流量變化直接影響流域內(nèi)水安全的保障，對流域可持續(xù)發(fā)展有著重要影響。研究表明，20世紀50年代以來，受氣候變化和水利水保工程等人類活動影響，黃河流域?qū)崪y徑流量呈顯著下降趨勢，且下游比上游衰減更顯著，其中人類活動是主要影響因素[5-9]；根據(jù)黃河流域水資源綜合評價相關(guān)成果，黃河流域多年平均地表水資源量也呈下降趨勢[4,10-11]。王國慶等[12]采用模型模擬方法對黃河流域花園口站1986—2018年系列(變化期)實測徑流量相對于1956—1985年系列(基準期)的變化量進行了歸因分析，得出黃河流域人類活動影響對徑流減少的貢獻率為76.2%，氣候變化影響貢獻率為23.8%。Kong等[13]采用雙累積曲線法分析了黃河流域1960s、1970s、1980s、1990s和2001—2012年相對于1951—1959年基準期氣候變化和人類活動對實測徑流變化的貢獻率，得出除1960s外黃河流域90%以上的徑流減少可歸因于人類活動影響。Chang等[14]采用模型模擬分析法對渭河流域各年代相對于基準期(1960s)之間的實測徑流變化進行歸因分析，得出渭河流域1970s、1980s、1990s和2000s人類活動影響貢獻率分別為64%、72%、47%和90%。馮家豪等[15]使用雙累積曲線對1957—2018年系列頭道拐—龍門與龍門—潼關(guān)區(qū)間突變年份(分別為1987年和1999年)前后實測徑流量變化進行歸因分析，得出人類活動對頭道拐—龍門與龍門—潼關(guān)區(qū)間多年平均徑流量變化貢獻率分別為63%和80%。以上研究主要針對實測徑流演變原因開展分析，不能解釋黃河天然徑流量衰減的原因，不利于黃河流域水資源配置與高效利用?？紤]未來氣候變化，黃河流域水資源量可能還會進一步減少[16-17]，在考慮輸沙水量、生態(tài)水量、河口三角洲需水量等多重需求下，未來黃河水資源短缺現(xiàn)象仍較為嚴峻[18]。因此，有必要深入研究氣候變化和人類活動對天然河川徑流量變化的影響，定量計算不同人類活動因素對天然河川徑流量減少的貢獻率，以便提出針對性的應(yīng)對措施，為新時期黃河流域生態(tài)保護與高質(zhì)量發(fā)展提供支持。

本研究采用基于二元水循環(huán)模型的流域天然河川徑流量評價方法[19-20]和多因素歸因分析方法[21-22]，在動態(tài)評價黃河流域天然河川徑流量的基礎(chǔ)上，定量分析氣候變化和人類活動對天然河川徑流量衰減造成的影響。

1 研究區(qū)域概況

黃河是中國第二大河流，流經(jīng)青藏高原、黃土高原、內(nèi)蒙古高原以及華北平原，橫跨中國三大階地，于山東流入渤海，干流全長5 464 km，呈“幾”字形，流域面積79萬km2。流域內(nèi)自西向東氣候分區(qū)為干旱區(qū)、半干旱區(qū)和半濕潤地區(qū)，流域1956—2016年多年平均降水量為452 mm，降水年內(nèi)分配極不均勻，主要集中在6—9月[9]。黃河流域主要水系及二級區(qū)范圍如圖1所示。

2 方法及數(shù)據(jù)

2.1 基于二元水循環(huán)模型的天然河川徑流量評價方法

天然河川徑流量是衡量流域可用水資源的一個重要指標，指的是流域內(nèi)地表水體中由降水形成的、通過流域出口斷面可逐年更新的全部動態(tài)水量，包括降水直接產(chǎn)流形成的地表水和地下水進入河道的基流部分，并扣除沿途蒸發(fā)、河流滲漏損失的水量，反映了河川徑流的供水能力，也是流域降水及產(chǎn)匯流條件的綜合反映[23-24]。對于無人類活動影響的流域，實測河川徑流量即為天然河川徑流量；對于有人類活動影響的流域，一般需要在實測河川徑流量基礎(chǔ)上還原農(nóng)業(yè)灌溉、工業(yè)和生活用水的耗損量、跨流域引入引出水量、河道分洪水量、水庫蓄變量等；此外，還需要對還原后的天然河川徑流量系列進行一致性修正，反映下墊面條件變化對徑流的影響，得到具有一致性且能反映近期下墊面條件的天然河川徑流量。這就是傳統(tǒng)的“實測+還原+還現(xiàn)”計算方法[23-24]。

傳統(tǒng)天然河川徑流量計算方法主要在實測徑流基礎(chǔ)上還原人類活動直接影響量，沒有考慮人類取用水等活動對地表、地下水循環(huán)過程的間接影響，比如開采地下水會造成地下水位下降、包氣帶加厚、降水入滲補給地下水的路徑延長、土壤水蒸散發(fā)增加，從而導致降水產(chǎn)生的地表、地下產(chǎn)流減少及天然河川徑流量減少。為考慮人類活動對天然河川徑流量的直接和間接影響，以基于二元水循環(huán)模型的流域水資源評價方法[19-20]為基礎(chǔ)，本文通過構(gòu)建具有物理機制的流域分布式二元水循環(huán)模型，對人類活動影響下流域“地表-土壤-地下”水循環(huán)過程進行連續(xù)模擬，根據(jù)天然河川徑流量的概念從模擬得到的各環(huán)節(jié)水循環(huán)變量中分離出因降水而形成的河川徑流量，即“耦合模擬-分離評價”的方法。

WEP-L分布式流域二元水循環(huán)模型(簡稱WEP-L模型)是具有物理機制的流域分布式水循環(huán)模型，可以綜合考慮氣象、下墊面、人類取用水、水利水保工程等因素對水循環(huán)過程的影響，實現(xiàn)“自然-社會”二元水循環(huán)過程耦合模擬和分析，給出水循環(huán)要素時間和空間變化過程以及流域水循環(huán)通量。模型采用子流域套等高帶作為計算單元，同時根據(jù)土地利用不均勻性又進行了細分，各計算單元在垂直方向上分為數(shù)層，用于模擬地表、包氣帶、地下含水層等不同介質(zhì)中的水分運動。自然水循環(huán)過程主要是對各計算單元的降水、入滲、產(chǎn)流、匯流等過程進行模擬；社會水循環(huán)模擬就是對社會水循環(huán)各個子系統(tǒng)的“蓄水-取水-輸水-用水-耗水-排水”6個環(huán)節(jié)進行模擬。社會水循環(huán)同自然水循環(huán)的耦合過程主要是在取水和排水環(huán)節(jié)，取水來源可分為河道取水、水庫取水以及地下水開采3類，取水環(huán)節(jié)主要會導致取水點所在計算單元地表和地下水減少，排水環(huán)節(jié)則會導致排水用戶所在計算單元地表徑流增加。詳細模型原理內(nèi)容參考文獻[19,25-26]。

本研究針對黃河流域梯田、淤地壩、水庫眾多的特點，根據(jù)對水循環(huán)過程調(diào)蓄的機制，對原來的模型進行改進，主要包括：① 將原模型的5類下墊面改進為水域、不透水域、裸地域、林地域、草地域、坡耕地域、灌溉農(nóng)田域、非灌溉農(nóng)田域、梯田域、壩地域共10類進行降水、蒸發(fā)、入滲、產(chǎn)流等水循環(huán)過程模擬。② 在梯田域增加梯田建設(shè)對坡面產(chǎn)匯流過程的模擬。③ 在壩地域增加淤地壩建設(shè)對溝壑產(chǎn)匯流過程的模擬。④ 對原模型水庫調(diào)度模塊進行改進，根據(jù)資料情況采用不同的水庫泄流量計算方法進行模擬，如果有水庫月蓄變量數(shù)據(jù)，則采用月蓄變量數(shù)據(jù)對水庫泄流量進行修正；如果有水庫庫容—下泄量曲線，則根據(jù)水庫庫容計算水庫泄流量；如果上述數(shù)據(jù)資料都缺乏，則采用多年月平均泄流量作為水庫泄流控制。

2.1.2 分離評價

“分離”評價具體步驟如下：

(1) 基于WEP-L模型計算結(jié)果，統(tǒng)計區(qū)域/流域內(nèi)各子流域月平均地表產(chǎn)流量、壤中流產(chǎn)流量、基流產(chǎn)流量、地表用水量(分農(nóng)業(yè)和工業(yè)生活2類)、地下水用水量(分農(nóng)業(yè)和工業(yè)生活2類)、降水量、河道蒸發(fā)量、地表入滲量等。

(2) 計算來自降水的天然地表產(chǎn)流量。首先根據(jù)降水量和農(nóng)業(yè)用水量的比例關(guān)系確定地表產(chǎn)流中降水所占的比重，然后采用地表產(chǎn)流中降水所占的比重和地表產(chǎn)流量確定天然地表產(chǎn)流量。

“教學”主要反映教師的教學水平、教學創(chuàng)新能力，以及教學成果的社會影響。有關(guān)評價要素包括：在各級教學成果評比、教學大賽中獲得的成績，教改課題，教學質(zhì)量評價（大數(shù)據(jù)反映的教師所任教班級學科成績的變化、教師教學情況等）；各級示范課、展示課等。其中，“教改課題”的認可度為中上等（事后訪談得知，有不少調(diào)查對象將該項列入了科研類），其他各評價要素的認可度都為高。

(3) 計算來自降水的天然壤中流產(chǎn)流量。首先根據(jù)降水入滲量、農(nóng)業(yè)灌溉入滲量和用水輸水過程中的滲漏量計算地下產(chǎn)流中降水所占的比重，然后采用地下產(chǎn)流中降水所占的比重和壤中流產(chǎn)流量確定天然壤中流產(chǎn)流量。

(4) 計算來自降水的天然基流產(chǎn)流量。采用地下產(chǎn)流中降水所占的比重和基流產(chǎn)流量確定天然基流產(chǎn)流量。

(5) 計算天然河川徑流量。天然河川徑流量為來自降水的天然地表產(chǎn)流量、天然壤中流產(chǎn)流量與天然基流產(chǎn)流量之和減去河道蒸發(fā)損失量。

2.2 水循環(huán)演變多因素歸因分析方法

本研究采用多因素歸因分析方法[21-22]進行歸因分析，根據(jù)影響水循環(huán)的n個因素分別處于基準期和變化期2種狀態(tài)(基準期和變化期時間序列可重疊)，共設(shè)置2n個模擬情景；使用分布式水循環(huán)模型對每種情景分別進行模擬，求得各情景條件下某個水循環(huán)變量的多年平均值；最后采用以下公式進行貢獻率分解計算：

(1)

(2)

式中：δxi表示因素xi對某個水循環(huán)變量的貢獻量；ei,j是對應(yīng)第i個因素在第j個情景下的權(quán)重系數(shù)(如果該情景下該因素處于變化期則為1，基準期則為-1)；Sj是第j個情景下模擬得到的某個水循環(huán)變量的多年平均值；βi表示第i個因素占總變化的貢獻率，正值表示起到增加作用，負值表示起到減少作用。

2.3 黃河模型數(shù)據(jù)處理

黃河模型輸入數(shù)據(jù)處理主要包括：

(1) 采用30″分辨率的GTOPO30 DEM數(shù)據(jù)作為基礎(chǔ)輸入，并基于此提取模擬河網(wǎng)以及8 485個子流域，并進一步劃分為38 720個等高帶。

(2) 采用407個國家氣象站和1 245個雨量站逐日降水、氣溫、濕度、風速、日照時數(shù)作為基礎(chǔ)輸入，并采用反距離平方法將站點數(shù)據(jù)展布到子流域形心上作為模型驅(qū)動。

(3) 以中科院地理所提供的30 m分辨率的土地利用數(shù)據(jù)為基礎(chǔ)，并重采樣為1 km，共1980年、1990年、2000年、2010年及2015年5期。

(4) 土壤及其特征信息采用全國第二次土壤普查資料。

(5) 主要考慮黃河流域119處10萬畝以上的大型灌區(qū)，數(shù)據(jù)來源黃河水利委員會相關(guān)匯編資料。

(6) 用水數(shù)據(jù)來源于黃河流域水資源公報，并采用向下尺度化將省套三級區(qū)尺度的數(shù)據(jù)展布到等高帶。

3 結(jié)果及討論

3.1 模型率定驗證

黃河模型構(gòu)建過程參見文獻[19]。本次研究將模擬系列延長到2016年，以黃河流域主要干流斷面控制區(qū)間作為參數(shù)分區(qū)，基于各水文站實測月平均流量對模型進行了重新率定和驗證，調(diào)試參數(shù)主要包括洼地儲留深、土壤飽和導水系數(shù)、河床材質(zhì)水力傳導系數(shù)、土壤層厚度、氣孔阻抗等敏感參數(shù)。其中，1956—1980年為率定期，1981—2016年為驗證期。采用2個標準對月徑流量進行評價，模擬月徑流量相對誤差盡可能??；模擬月徑流Nash-Sutcliffe效率系數(shù)盡可能大[27]。計算公式如下所示：

(3)

(4)

各站流量過程模擬結(jié)果和效率系數(shù)如圖2和表1所示。率定期黃河主要干流水文站月流量過程Nash-Sutcliffe效率系數(shù)均在0.7以上，大部分站點相對誤差均在6%以內(nèi)；驗證期Nash-Sutcliffe效率系數(shù)除頭道拐和龍門站外，均在0.6以上，相對誤差除龍門站外，均在3%以內(nèi)。結(jié)果表明，WEP-L模型總體上能夠較好地描述黃河流域的徑流過程。相比驗證期，率定期模型模擬效果稍好，原因在于此階段人類活動較少。

圖2 主要干流站點月徑流過程Fig.2 Process of monthly streamflow of main station in Yellow River

表1 各站月徑流過程效率系數(shù)Table 1 Model efficiency coefficient of monthly streamflow in main station of Yellow River

3.2 不同水平年天然河川徑流量評價

黃河水資源量評價相關(guān)研究成果主要有“八七”分水方案[10]、黃河流域第二次水資源調(diào)查評價[4]、黃河流域第三次水資源調(diào)查評價等。為了保持和已有成果評價口徑一致，本次研究主要設(shè)置3個評價情景，分別為1956—1979年水平年情景、2000年水平年情景、2016年水平年情景，詳細設(shè)置信息如表2所示。

表2 評價情景設(shè)置Table 2 Scene setting of model evaluation

表3是主要干流斷面和區(qū)間天然河川徑流量評價結(jié)果和黃河水資源調(diào)查評價結(jié)果對比。從表中可以看出，本文評價結(jié)果和黃河水資源評價結(jié)果偏差基本都在10%以內(nèi)，可以認為2種評價結(jié)果總體一致。分析主要干流斷面和區(qū)間2000年水平年和2016年水平年本文評價結(jié)果相對1956—1979年水平年變化量，結(jié)果表明，對黃河流域花園口斷面而言，本次評價的3個情景，天然河川徑流量分別為567.2億m3、486.9億m3和452.6億m3，呈減少趨勢，且2016年水平年情景相對于1956—1979年水平年減少20.2%。對各區(qū)間而言，相比于1956—1979年水平年，2000年水平年和2016年水平年情景除唐乃亥斷面以上區(qū)域有所增加之外，其他區(qū)間天然河川徑流量均呈減少趨勢，其中，龍門—三門峽減少幅度最大，分別減少31.5億m3和37.7億m3；其次是蘭州—頭道拐和頭道拐—龍門區(qū)間；三門峽—花園口區(qū)間減少幅度最小，分別為11.4億m3和15.3億m3。

表3 主要干流斷面和區(qū)間本文評價結(jié)果與調(diào)查評價成果對比 億m3

3.3 天然河川徑流量變化歸因分析

為了支撐黃河流域生態(tài)保護和高質(zhì)量發(fā)展，有必要深入研究氣候變化和人類活動對水資源變化的影響，定量計算不同因素對水資源減少的貢獻率。因此，本研究重點分析相對黃河“八七”分水方案而言，第三次水資源評價情景下黃河流域天然河川徑流量減少的貢獻，即采用多因素歸因分析方法分析氣候、下墊面(包括植被、梯田、壩地)以及社會經(jīng)濟取用水(包括地表、地下取用水)3類因素對2016年水平年相對1956—1979年水平年天然河川徑流量減少的貢獻率。情景設(shè)置如表4所示，由于需要考慮3類因素，因此設(shè)置23共8個情景?；鶞势谂c變化期不同因素的設(shè)置見表5。主要影響因素多年平均值如表6所示。

表4 多因素歸因分析法情景設(shè)置Table 4 Scene setting of multi factor attribution method

表5 基準期與變化期不同因素的設(shè)置

表6 黃河主要干流斷面以上區(qū)域主要水循環(huán)影響因素均值

采用WEP-L模型對各情景分別進行模擬，并采用公式(1)和公式(2)進行貢獻率計算，主要干流斷面2016年水平年情景相比于1956—1979年情景天然河川徑流量歸因分析結(jié)果如表7所示。從表中可以看出，對黃河流域(花園口斷面以上)而言，氣候變化的貢獻率占24.4%，其余均由人類活動影響引起，其中社會經(jīng)濟取用水占50.6%，下墊面變化占25.0%。對各區(qū)間而言，唐乃亥以上區(qū)域降水增加、氣溫升高，水土保持等下墊面變化幅度小、社會經(jīng)濟取用水變化很小，所以氣候變化是導致唐乃亥以上區(qū)域天然河川徑流量增加的主要原因；唐乃亥以下各區(qū)間，天然河川徑流量均在減小，但各種因素的貢獻率各不相同。唐乃亥—蘭州區(qū)間降水增幅小于唐乃亥以上區(qū)域，氣溫增幅大于唐乃亥以上區(qū)域，氣候變化造成天然河川徑流量減少，同時用水量和水土保持面積增量不大，所以氣候變化的貢獻率占比較大，達到60.5%；蘭州—頭道拐和三門峽—花園口區(qū)間用水量增幅較大，社會經(jīng)濟取用水的貢獻率分別達到65.4%和70.5%；頭道拐—龍門和龍門—三門峽區(qū)間水土保持面積增幅為全流域最大，水土保持的貢獻率為全流域最高，達到1/3左右。對比不同分區(qū)下墊面變化貢獻率，可以得出，在影響天然河川徑流量的下墊面各因素中，水土保持(梯田、淤地壩)是貢獻最大的因素，但是涉及到植被、梯田、淤地壩分項因素定量的研究成果，需要在今后的研究中進一步深入論證。

表7 黃河主要干流斷面和區(qū)間天然河川徑流量變化歸因分析Table 7 Attribution result of natural streamflow in the main station and zone of Yellow River

以往針對黃河流域徑流變化歸因分析的研究中，均為對實測徑流變化的歸因分析，且基準期和變化期不重疊?？紤]到本研究的變化期(1956—2016年)與基準期(1956—1979年)存在包含關(guān)系，因此，如果對比1980—2016年相對于1956—1979年人類活動對天然徑流量變化的貢獻，應(yīng)比本文研究得到的貢獻率75%還要大。另一方面，考慮到人類活動對實測徑流既有直接影響，又有間接影響，而對天然河川徑流只有間接影響，人類活動對實測徑流影響的貢獻應(yīng)比天然河川徑流還要大一些。因此，按照本文分析結(jié)果，對于花園口斷面，人類活動對實測徑流變化的影響應(yīng)大于75%。參照前人對黃河流域?qū)崪y徑流變化影響的研究結(jié)果[12-15,28]，人類活動的影響均占主導作用，大多數(shù)研究得到貢獻率在70%以上，有的甚至在90%以上，這與本文結(jié)果基本一致。

4 結(jié) 論

為進一步區(qū)分氣候變化和不同人類活動對天然河川徑流量變化的貢獻，采用基于二元水循環(huán)模型的流域天然河川徑流量評價方法和水循環(huán)演變多因素歸因分析方法進行分析。主要結(jié)論如下：

(1) 采用1956—2016年長系列資料，構(gòu)建黃河流域分布式二元水循環(huán)模型，并對模型進行率定驗證。結(jié)果表明，對黃河主要干流水文站月流量過程而言，率定期Nash-Sutcliffe效率系數(shù)均在0.7以上，大部分站點相對誤差均在6%以內(nèi)；驗證期Nash-Sutcliffe效率系數(shù)除頭道拐和龍門站外，均在0.6以上，相對誤差除龍門站外，均在3%以內(nèi)。說明模型能夠較好地描述黃河流域?qū)崪y徑流過程。

(2) 采用基于分布式二元水循環(huán)模型的“耦合模擬-分離評價”的天然河川徑流量計算方法，評價分析黃河流域主要干流斷面天然河川徑流量。結(jié)果表明，在1956—1979年水平年、2000年水平年、2016年水平年情景下，花園口斷面天然河川徑流量分別為567.2億m3、486.9億m3和452.6億m3，呈減少趨勢，2016水平年情景相對于1956—1979年水平年減少20.2%。2000年水平年和2016年水平年情景除唐乃亥斷面以上區(qū)域天然河川徑流量有所增加之外，其他區(qū)間均呈減少趨勢，其中龍門—三門峽區(qū)間減少幅度最大，分別減少31.5億m3和37.7億m3；其次是蘭州—頭道拐和頭道拐—龍門區(qū)間；三門峽—花園口區(qū)間減少幅度最小，分別為11.4億m3和15.3億m3。

(3) 采用水循環(huán)演變多因素歸因分析方法分析氣候、下墊面和社會經(jīng)濟取用水3個因素對天然河川徑流量減少的貢獻率。結(jié)果表明，花園口斷面2016年水平年情景(現(xiàn)狀年情景)相比于1956—1979年水平年情景(“八七”分水情景)，天然河川徑流量減少114.6億m3，其中氣候變化、下墊面變化和社會經(jīng)濟取用水的影響貢獻率分別為24.4%、25.0%和50.6%。從分區(qū)來看，蘭州以上區(qū)域氣候變化是主導因素，蘭州以下各區(qū)間人類活動是主導因素，受水土保持和社會經(jīng)濟取用水活動程度不同，各種因素的貢獻率各不相同。說明人類活動對天然河川徑流量的間接影響(因用水新增區(qū)域蒸散發(fā)和滯留包氣帶)要大于氣候變化的直接影響(降水量減少和因氣溫升高增加蒸散發(fā))。

(4) 從研究結(jié)果來看，黃河流域天然河川徑流量衰減是大趨勢，人類取用水和水土保持等活動是主要因素，應(yīng)加強深度節(jié)水、剛性控水、適度增水、強化管水和立法護水，以遏制天然河川徑流量衰減的趨勢，促進流域生態(tài)保護與高質(zhì)量發(fā)展。

致謝:感謝黃河水利委員會水文局的張學成正高級工程師和李東正高級工程師在研究過程中給予的幫助。