曾思棟，夏 軍，杜 鴻，陳向東

（1. 中國科學(xué)院重慶綠色智能技術(shù)研究院，重慶 400714；2. 長江勘測規(guī)劃設(shè)計研究院，湖北 武漢 430010；3. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，湖北 武漢 430072；4. 中南民族大學(xué) 資源與環(huán)境學(xué)院，湖北 武漢 430074；5. 中國水權(quán)交易所，北京 100053）

1 研究背景

水循環(huán)作為地球系統(tǒng)中各過程物質(zhì)交換和能量轉(zhuǎn)換的重要載體和紐帶，與陸地生態(tài)系統(tǒng)過程相互作用并發(fā)生反饋［1］。一方面水循環(huán)通過降水、下滲等過程為生態(tài)系統(tǒng)提供水分，影響植物光合作用碳同化、CO2排放進而影響生態(tài)系統(tǒng)的格局［2］；另一方面陸地生態(tài)系統(tǒng)中植被通過氣孔行為調(diào)控植物的蒸騰，植被動態(tài)變化如葉面積指數(shù)變化影響降水截留和蒸發(fā)、能量分配，對水文過程產(chǎn)生影響［3-4］。研究陸地水循環(huán)與生態(tài)系統(tǒng)的耦合是生態(tài)水文的基本內(nèi)容，準確認識生態(tài)水文的耦合機制是實現(xiàn)流域水資源及生態(tài)管理的重要基礎(chǔ)［5］。目前，生態(tài)水文的耦合研究已經(jīng)成為水文、氣象、生態(tài)等多個學(xué)科共同關(guān)注的課題。

陸地生態(tài)水文研究的主要內(nèi)容體現(xiàn)為植被與水分、能量和物質(zhì)的耦合循環(huán)機理研究［6］，構(gòu)建生態(tài)水文耦合模型是認識這些耦合機理包括水、能量、碳、氮、磷等相互作用的重要基礎(chǔ)。陸地生態(tài)水文耦合從單向耦合向雙向耦合發(fā)展，單向耦合模型主要是從水文模擬角度出發(fā)，考慮植被對水文過程的單向影響，不考慮水文過程對植被動態(tài)的反饋，如DHSVM（Distributed Hydrology-Soil-Vegeta?tion Model）模型［7］等。隨著對生態(tài)水文過程耦合認識的不斷加強，雙向耦合模型的理論框架被提出，如趙風華等［8］根據(jù)生態(tài)水文水、碳循環(huán)要素在植被-土壤-大氣連續(xù)體（SPAC）中的運動過程，將耦合作用分為4 個尺度的關(guān)聯(lián)。在雙向耦合模型的構(gòu)建中，由于傳統(tǒng)水文模型對生態(tài)過程的描述不夠充分，生物地球化學(xué)模型對水循環(huán)過程的描述也欠缺，因此加強分布式水文模型與生物地球化學(xué)模型的耦合可以更好地研究陸地生態(tài)系統(tǒng)與水循環(huán)的相互作用［9］，如Mo等［10］構(gòu)建了基于能量平衡、水循環(huán)和碳氮循環(huán)的生態(tài)水文動力的VIP模型（Vegetation Interface Processes Model）；Park等［11］建立了耦合水文模型與生物地球化學(xué)模型的CHANGE 模型，并應(yīng)用于西伯利亞落葉松林的水、碳通量和能量的模擬；Shen 等［12］將PAWS（Process-based Adaptive Watershed Simulator）模型與CLM 模型（Community Land Model）耦合分析不同土地利用/植被覆蓋以及地下水對蒸散發(fā)及凈初級生產(chǎn)力的影響。

具有物理機制的陸面水文模型和生物地球化學(xué)模型分別是描述陸面水文過程及生物地球化學(xué)過程的有效工具，將二者進行完全的緊密耦合對于認識陸面生態(tài)水文相互作用及其反饋機制有著重要的基礎(chǔ)意義。目前仍然需要進一步加強基于物理機制的生態(tài)水文過程“雙向”緊密耦合模型的研究，本研究將分布式時變增益水文模型（Distributed Time Variant Gain Model，DTVGM）與生物地球化學(xué)循環(huán)模型（Carnegie-Ames-Stanford Approach Carbon Nitrogen Phosphorous，CASACNP）進行雙向緊密耦合構(gòu)建DTVGM-CASACNP 模型，提出了耦合模型的架構(gòu)、模擬方法及主要原理和方法，為陸面水文-生物地球化學(xué)過程的耦合模擬提供基礎(chǔ)支撐。

2 生態(tài)水文雙向耦合模型架構(gòu)

生態(tài)水文雙向耦合模型DTVGM-CASACNP 由流域水循環(huán)模塊、能量平衡模塊、光合作用模塊、碳氮磷生物地球化學(xué)循環(huán)模塊等4 大部分組成。其中流域水循環(huán)模塊主要包括垂直方向大氣水-植物水-地表水-土壤水-地下水等的轉(zhuǎn)化以及水平方向的徑流傳輸模擬；能量平衡主要包括冠層輻射傳輸、冠層及地表能量平衡以及土壤熱傳導(dǎo)中能量轉(zhuǎn)化過程的模擬；光合作用模塊主要包括光合有效輻射吸收、光合作用及氣孔導(dǎo)度等碳同化相關(guān)過程的模擬，碳、氮、磷生物地球化學(xué)循環(huán)模塊主要包括植物生長，碳、氮、磷元素在不同生物地球化學(xué)庫包括植物庫、凋落物庫以及土壤庫中的相互轉(zhuǎn)化過程的模擬。

DTVGM-CASACNP 生態(tài)水文雙向耦合模型基于DTVGM 和CASACNP 模型構(gòu)建。DTVGM 模型是夏軍提出并發(fā)展的通過結(jié)合空間數(shù)字化信息將集總水文模型拓展應(yīng)用到流域時空變化模擬的分布式水文模型［13］，該模型將水文非線性系統(tǒng)理論與分布式模型框架進行有機結(jié)合，對于不同流域、水文資料信息不完整或不確定性干擾下的水文模擬具有較強的適應(yīng)能力，已經(jīng)在多個流域得到較好的驗證［14-15］。為了描述生態(tài)水文耦合的相互作用過程，本研究對DTVGM 進行改進，將DTVGM 中基于水量平衡的兩層土壤水模塊改為基于Richards方程的土壤水數(shù)值模擬模型，將基于Bagrov的蒸散發(fā)模塊改為土壤-植被雙源蒸散發(fā)模型，增加了能量平衡模擬以及土壤熱運動模塊。生物地球化學(xué)模型CA?SACNP是由Wang等［16］和Houlton等［17］基于CASA（Carnegie-Ames-Stanford Approach）模型發(fā)展而來的，可以模擬植被動態(tài)過程以及碳、氮、磷的循環(huán)，該模型已經(jīng)用于全球或區(qū)域碳、氮、磷循環(huán)的研究，并得到了較好的驗證［18］。本耦合模型中植物的光合作用采用基于瞬時尺度的Farquhar［19］光合作用模型，氣孔導(dǎo)度模型采用具有一定物理機制的Leuning［20］冠層導(dǎo)度模型。

生態(tài)水文雙向耦合模型各模塊主要運行機制如下：（1）流域水循環(huán)模塊，主要計算各單元內(nèi)的冠層截留、蒸散發(fā)、產(chǎn)流、土壤水運動、坡面匯流以及河網(wǎng)匯流等自然水文過程；（2）能量平衡模塊，主要計算各單元冠層及土壤各項能量收支項包括短波輻射、長波輻射，潛熱通量、感熱通量、土壤熱通量等能量通量以及冠層和地表、土壤溫度變化等；（3）光合作用模塊，主要計算冠層吸收的光合有效輻射、水汽及CO2氣孔導(dǎo)度、光合作用凈光合速率、初級生產(chǎn)力（GPP）等碳同化相關(guān)過程；（4）碳、氮、磷生物地球化學(xué)循環(huán)模塊，主要計算植物自養(yǎng)呼吸、凈初級生產(chǎn)力（NPP）及其分配、葉面積指數(shù)，碳、氮、磷在植物庫、調(diào)落物庫以及土壤庫中的相互轉(zhuǎn)化通量以及異養(yǎng)呼吸等生物地球化學(xué)過程。

DTVGM-CASACNP 以流域水循環(huán)為紐帶，耦合了伴隨的能量通量以及碳、氮、磷的生物地球化學(xué)過程，通過本架構(gòu)實現(xiàn)水、能量通量、碳、氮、磷過程的緊密耦合以及反映其相互作用及反饋的實時動態(tài)關(guān)系，可以實現(xiàn)陸面水文-生物地球化學(xué)過程的耦合模擬，為陸地水、碳、氮、磷循環(huán)研究及水資源及生態(tài)管理提供基礎(chǔ)支撐。圖1為耦合模型的結(jié)構(gòu)示意圖。

圖1 生態(tài)水文雙向耦合模型DTVGM-CASACNP結(jié)構(gòu)示意圖

3 雙向耦合模擬方法

陸地生態(tài)水文的雙向緊密耦合模擬體現(xiàn)在水、能量通量、碳、氮、磷通量在每個計算時段內(nèi)的相互作用及實時反饋，這些相互作用及實時反饋的描述通過關(guān)鍵的耦合變量實時傳遞實現(xiàn)。水、能量通量及碳、氮、磷變量實時傳遞主要包括：（1）能量過程為水文過程提供冠層和地表輻射傳輸、土壤熱通量，直接影響水文過程蒸散發(fā)、土壤水運動等；水文過程為能量過程提供蒸散發(fā)、土壤含水量，直接影響能量傳輸過程的潛熱通量、土壤熱通量等；（2）能量過程為碳、氮、磷過程提供冠層和地表輻射，冠層、地表及土壤溫度，直接影響植物光合作用、呼吸作用等；碳、氮、磷過程為能量過程提供植物葉面積指數(shù)等，直接影響能量的分配；（3）水文過程為碳、氮、磷過程提供土壤含水量、葉面濕潤比例等，直接影響植物光合作用、植物干物質(zhì)分配等；碳、氮、磷過程為水文過程提供葉面積指數(shù)、氣孔導(dǎo)度等，直接影響蒸散發(fā)、冠層截留等。水、能量通量及碳、氮、磷通量主要變量實時傳遞關(guān)系與模擬流程，如圖2、圖3所示。

圖2 耦合模型主要變量傳遞示意圖

圖3 雙向耦合模擬基本方法流程圖

實時動態(tài)反饋的雙向耦合模擬具體算法如下：（1）將模擬時段劃分為M個時段，對每個時段分別進行水、能量通量、碳、氮、磷通量的耦合模擬；（2）設(shè)定每個時段初始冠層溫度、地表溫度，進行冠層和地表輻射傳輸計算；（3）進行植物光合作用-氣孔導(dǎo)度-蒸騰作用的耦合模擬計算；（4）計算地表產(chǎn)流，土壤水、熱運動；（5）計算植物、凋落物和土壤中碳、氮、磷循環(huán)過程；（6）根據(jù)能量平衡迭代求冠層、地表溫度，滿足收斂條件結(jié)束本時段計算，不滿足收斂條件回到步驟（2）直到能量平衡收斂為止；（7）進行下時段耦合模擬計算，直到所有計算時段完成。

4 模型主要原理

4.1 地表產(chǎn)流過程地表水產(chǎn)流根據(jù)時變增益非線性產(chǎn)流概念計算，其理論基礎(chǔ)是降雨-徑流系統(tǒng)的非線性關(guān)系，表征為與前期影響雨量指數(shù)（API）有關(guān)的時變增益因子［13］。計算公式如下：

式中：Rs為地表產(chǎn)流，mm；Ke為反映土壤平均含水量變化的參數(shù)；Pe為有效降雨，mm；g1、g2為時變增益因子。

4.2 蒸散發(fā)過程蒸散發(fā)是陸地水汽返回大氣的主要形式，直接影響地表能量的收支平衡及反饋大氣的潛熱通量。DTVGM 原版模型中將土壤和冠層作為一個整體并通過潛在蒸散發(fā)與土壤水分脅迫來計算實際蒸散發(fā)，而陸地生態(tài)水文耦合的關(guān)鍵是基于植被氣孔行為的蒸騰作用與光合作用的耦合。因此本研究中將蒸散發(fā)進行分層細化，分為土壤蒸發(fā)、冠層截留蒸發(fā)和植被蒸騰，并引入莫興國等［21］對雙源蒸散發(fā)模型的改進方法，計算公式如下：

4.3 土壤水運動過程土壤水運動不僅是水循環(huán)的重要環(huán)節(jié)，對于植被的蒸騰、土壤蒸發(fā)、植被生長以及碳、氮、磷的生物地球化學(xué)循環(huán)也起著關(guān)鍵性的作用。為細致反映土壤水的運動規(guī)律，本研究土壤水運動基于Richards方程［22］，該方程微分形式如下：

式中：θ為土壤含水率，m3/m3；ψ為土壤基質(zhì)勢，m；K為土壤水力傳導(dǎo)率，mm/d；R為根系吸水項，mm；t為時間，s；z為土壤深度，m。采用van Genuchten 公式［23］描述具有高度非線性的水分特征曲線θ(ψ) 和水力傳導(dǎo)函數(shù)K(ψ) 如下：

式中：θr、θs分別為土壤殘余含水率和土壤飽和含水率，m3/m3；α、n為水分特征曲線參數(shù)。

Richards 方程的隱式差分離散格式見方程（9），模型求解采用變時間步長法，即當?shù)螖?shù)大于設(shè)定次數(shù)時，將模型的時間步長變?yōu)樯弦淮蔚?/2，如此反復(fù)直到前后兩次土壤負壓滿足設(shè)定精度為止。本次研究采用Celia［24］方法以減少計算過程中的水量平衡誤差，用追趕法求解該方程。

式中：j為時間；m為迭代次數(shù)；z為單元中心節(jié)點的縱坐標。

本次研究上、下邊界條件分別為：

4.4 植物根系吸水植被主要通過根系從土壤中吸收水分，滿足其植被的生長以及代謝等生理活動和葉片蒸騰。本研究假設(shè)根系吸水R(z,t) 是蒸騰Ec(t)以及有效根分布Le(z,t)的函數(shù)，在此基礎(chǔ)上對根系吸水加入土壤水分約束［25］。

4.5 冠層和地表輻射傳輸?shù)竭_地表的輻射包括太陽總輻射（短波輻射直射和散射）以及天空的長波輻射，同時地表又會反射太陽短波輻射，也會發(fā)出長波輻射。對于冠層凈輻射Rnc（單位為W/m2，以下同）平衡包括太陽短波輻射Rsc(↓,表示輻射傳輸方向向下,以下同)、冠層反射的短波輻射Rscr(↑,表示輻射傳輸方向向上,以下同)、天空長波輻射Rlc(↓)，冠層向天空Rlca(↑)及向地表Rlcs(↓)兩個方向發(fā)射的長波輻射以及地表長波輻射Rlsc(↑)。對于地表凈輻射Rns平衡包括冠層截留剩下的太陽短波輻射Rss(↓)、天空長波輻射Rls(↓)、冠層長波輻射Rlcs(↓)，地表反射的短波輻射Rssr(↑)以及地表發(fā)射的長波輻射Rlsac(↑)。

4.6 冠層和地表能量平衡冠層和地表能量平衡方程可表示為如下式：

式中：λEc、λEs分別為冠層蒸潛熱通量及土壤蒸發(fā)潛熱通量，W/m2；Hc、Hs分別為冠層感熱通量以及土壤表面感熱通量，W/m2；G為土壤熱通量，W/m2；Cc和Cgs分別是冠層及地表的容積熱容量，J/（m2·℃）；Tc和Tgs分別是冠層和地表溫度，℃；

其中，冠層和土壤表面的感熱通量［26］采用下式計算：

土壤熱通量通過對兩層土壤的熱量估計來計算：

式中：Tgs、Ts1、Ts2分別為在地表、土壤深度d1、d2處的土壤溫度，℃；k1、k2分別是土壤深度d1、d2處的土壤的熱傳導(dǎo)系數(shù)，W/（m·℃）。

4.7 土壤熱運動土壤溫度直接影響到能量平衡的各分量，土柱中熱量傳輸?shù)姆匠蹋?7］如下：

式中：C為土壤體積熱容量，J/（m3·℃）；T為土壤溫度，℃。初始條件可以是來自站點的實測數(shù)據(jù)，或者通過氣溫進行估計。上邊界條件為地表溫度，地表溫度由地表能量平衡計算中得出。下邊界條件根據(jù)實際情況可以選擇固定溫度Tb或者零通量邊界。

4.8 光合作用與氣孔導(dǎo)度本研究采用基于Furquhar 和Collatz［28］的光合作用模型計算植被碳同化，該模型具有一定的物理機制，得到了廣泛的應(yīng)用。該模型中葉子的凈光合速率A（umol C/（m2·s））采用RuBP 再生限制的凈光合作用率Aj、羧化速率限制的凈光合作用率Ac以及最大凈同化速率As的最小平滑值進行計算，即將凈光合速率看作是式（18）二次多項式的最小的根：

式中，Rubisco活性限制的光合速率為：

RuBP（光限制率）再生限制的光合速率為：

最大凈同化速率為：

式中：Vc,max為最大羧化作用率，umol C/（m2·s）；Kc和Ko分別為CO2和O2的Michaelis-Menten 系數(shù)，umol/mol；Γ*為CO2補償點，umol/mol；ci和oi分別為葉肉細胞間隙的CO2和O2濃度，umol/mol；J為電子傳輸速率，umol electrons/（m2· s）；β1、β2為 經(jīng)驗參數(shù)，C3 和C4 分別 指C3 植 物 和C4 植物。

蒸騰和光合作用均受到植物氣孔的開關(guān)控制，本次采用Leuning 在Ball-Berry 模型［20］基礎(chǔ)上改進的冠層導(dǎo)度模型，該模型基于CO2氣孔導(dǎo)度和水汽的氣孔導(dǎo)度成線性關(guān)系如下：

式中：gs為水汽的冠層導(dǎo)度，m/s；a1為經(jīng)驗系數(shù)；cs為葉片表面的CO2濃度，umol/mol；D為葉片表面的飽和水汽壓差，kPa；Dx為表征氣孔對敏感性的經(jīng)驗系數(shù)；g0為最小的冠層導(dǎo)度，m/s。

4.9 植物生長與物候模擬植物的自養(yǎng)呼吸包括維持呼吸Rm和生長呼吸Rg。維持呼吸［29］計算如下：

式中：cni為植物各庫的碳氮比；Ci為植物各庫的含碳量，umol C/m2；ri為根莖葉呼吸速率。生長呼吸通常與扣除維持呼吸的初級生產(chǎn)力（GPP，umol C/（m2·s））成正比，則凈初級生產(chǎn)力（NPP，umol C/（m2·s）及其分配如下：

式中：ac,i為NPP分配到根、莖、葉的分配系數(shù)；τP,i為植物各庫到凋落物庫的轉(zhuǎn)化速率，s-1。耦合模型中植被的物候分為4個階段，第1階段從發(fā)芽到穩(wěn)定生長，第2階段從穩(wěn)定生長到葉子開始凋落，第3階段從葉子凋落開始到凋落結(jié)束，第4階段從葉子凋落結(jié)束到發(fā)芽。

4.10 碳的生物地球化學(xué)過程植物凋落后進入凋落物庫，凋落物一部分被吸收到土壤有機碳庫中，一部分轉(zhuǎn)化為CO2。凋落物碳庫的控制方程如下：

式中：bj,i為碳從植物庫i到凋落物庫j的分配比例；mn為土壤無機氮限制消減率；τL,j為沒有被土壤無機氮限制的凋落物分解速率，s-1。

土壤有機碳分解時，一部分被其他有機碳庫吸收，一部分通過異氧呼吸轉(zhuǎn)化為CO2進入大氣，控制方程如下：

式中：ck,j為碳從凋落物碳庫到土壤碳庫的分配比例；dk,kk為碳從土壤碳庫kk到土壤碳庫k的比例；τs,k是各碳庫的分解速率。

4.11 氮的生物地球化學(xué)過程植物有機氮庫分為根、莖、葉三個庫，植物將從根系吸收的氮按比例分配到各庫。植物衰老過程中，一部分氮被重新吸收到活的組織中，剩下的轉(zhuǎn)化到凋落物庫中［30］。

式中：NP,i為植物各庫有機氮含量，umol N/m2；Fn,up為根系氮吸收率，umol N/（m2·s）；rn,i是重新吸收系數(shù)。

凋落物氮庫的控制方程如下。

式中：nL,str為STR庫的氮碳比；τL,i為各庫的分解速率，s-1；mn是氮限制因子。

土壤有機氮庫的控制方程如下：

式中：ck,j為氮從凋落物庫j到土壤有機氮庫k的比例；dk,kk為氮從土壤有機氮庫k到土壤有機氮庫kk的比例；τs,k為土壤有機氮各庫的分解速率，s-1。

無機氮庫被認為是植物根系吸氮的主要來源，其控制方程如下：

式中：Ns,min為土壤無機氮庫的量，umol N/m2；Fn,dep、Fn,fert、Fn,fix、Fn,net、Fn,up、Fn,loss分別為大氣氮沉降率、施肥添加率、固氮率、凈礦化率、氮吸收率、氮損失率，umol N/（m2·s）。當根系吸氮量小于氮的需求量時，NPP將會減少，根系吸氮計算如下：

式中：nPmax,i、nPmin,i分別為各植物庫的最大和最小氮碳比；KN,up為經(jīng)驗系數(shù)。

4.12 磷的生物地球化學(xué)過程磷通過大氣沉降、巖石風化進入生態(tài)系統(tǒng)，通過淋失和侵蝕輸出系統(tǒng)。有機磷植物各庫的控制方程如下［16］：

式中：PP,i為植物各庫中有機磷的量，umol P/m2；Fp,up為植被根系的磷吸收率，umol P/（m2· s）；rp,i為凋落前磷的被植物組織重新吸收系數(shù)。

凋落物氮庫的控制方程如下：

式中：pL,str為STR庫的磷碳比；τL,i為各庫的分解速率，s-1。

土壤有機磷庫的控制方程如下：

式中：ck,j為磷從凋落物庫j到土壤有機磷庫k的比例；dk,kk為磷從土壤有機磷庫k到土壤有機磷庫kk的比例；τs,k為土壤有機磷各庫的分解速率，s-1；Fp,tase為磷的生物化學(xué)礦化率。

土壤無機磷庫按照植被和細菌對無機磷吸收的有效程度分為3個部分，分別為不穩(wěn)定磷、吸附磷和強吸附磷。以不穩(wěn)定磷Ps,lab為例，其控制方程如下：

式中：Fp,dep、Fp,wea、Fp,fert、Fp,net、Fp,tase、Fp,up、Fp,loss分別代表磷沉降率、固定率、施肥添加率、凈礦化率、生物化學(xué)礦化率、吸收率、損失率，umol P/（m2· s）；Spmax為吸附磷的最大量；Kplab為吸附系數(shù)。

5 結(jié)論

（1）本研究提出了生態(tài)水文雙向緊密耦合的模擬方法，通過耦合的關(guān)鍵變量實時傳遞實現(xiàn)水、能量通量、碳、氮、磷通量在每個計算時段內(nèi)的相互作用與實時反饋。

（2）本研究提出的耦合模型能夠較好地刻畫水、能量及碳循環(huán)等陸面水文物理、生物地球化學(xué)過程及其相互作用，為陸面水文-生物地球化學(xué)過程的耦合模擬及應(yīng)用提供基礎(chǔ)支撐。

該生態(tài)水文雙向耦合模型DTVGM-CASACNP的應(yīng)用案例在下篇中給出。