鄭輝， 劉樹(shù)華*， Prabhakar Clement， 劉振鑫， 候旭宏，王姝, 趙靖川， 李源， 繆育聰， 鄭亦佳， 盛黎， 朱琳

1 北京大學(xué)物理學(xué)院大氣與海洋科學(xué)系， 北京　100871 2 Department of Civil Engineering, Auburn University, Auburn, AL 36849 3 中國(guó)科學(xué)院大氣物理研究所， 北京 100029 4 中國(guó)科學(xué)院寒區(qū)旱區(qū)環(huán)境與工程研究所， 蘭州　730000 5 中國(guó)氣象局國(guó)家氣象中心， 北京　100081 6 中國(guó)氣象局衛(wèi)星氣象中心， 北京　100081

?

北京大學(xué)陸面過(guò)程模式PKULM(Peking University Land Model)介紹及檢驗(yàn)

鄭輝1， 劉樹(shù)華1*， Prabhakar Clement2， 劉振鑫3， 候旭宏4，王姝1, 趙靖川1， 李源1， 繆育聰1， 鄭亦佳1， 盛黎5， 朱琳6

1 北京大學(xué)物理學(xué)院大氣與海洋科學(xué)系， 北京100871 2 Department of Civil Engineering, Auburn University, Auburn, AL 36849 3 中國(guó)科學(xué)院大氣物理研究所， 北京 100029 4 中國(guó)科學(xué)院寒區(qū)旱區(qū)環(huán)境與工程研究所， 蘭州730000 5 中國(guó)氣象局國(guó)家氣象中心， 北京100081 6 中國(guó)氣象局衛(wèi)星氣象中心， 北京100081

摘要陸面過(guò)程模式是氣候模式和天氣模式的核心組成部分之一.在土壤—植被—大氣耦合模式(Soil-Plant-Atmosphere Model, SPAM)的基礎(chǔ)上，發(fā)展了新一代北京大學(xué)陸面過(guò)程模式PKULM(Peking University Land Model).本文首先介紹了PKULM的輻射傳輸、湍流輸送、光合作用、土壤水熱輸送等過(guò)程的參數(shù)化方案；采用隱式迭代計(jì)算框架，發(fā)展并應(yīng)用了一個(gè)快速的線性方程組求解算法，提高了模式計(jì)算穩(wěn)定性；提出并使用了二分搜索算法計(jì)算氣孔阻抗，避免了CLM(Community Land Model)等使用的迭代方法在干旱區(qū)不穩(wěn)定的情況，提高了模式的適用性；采用水勢(shì)為基礎(chǔ)的土壤水分?jǐn)U散方程，使模式能夠模擬土壤飽和區(qū)的水分輸送過(guò)程，為進(jìn)一步與水文過(guò)程模式耦合奠定了基礎(chǔ)；還發(fā)展了一個(gè)地表積水與徑流過(guò)程的機(jī)理模型，提高了模式對(duì)地表水分平衡過(guò)程的模擬能力；最后，使用“中國(guó)西北干旱區(qū)陸—?dú)庀嗷プ饔糜^測(cè)試驗(yàn)”平?jīng)稣镜馁Y料對(duì)模式進(jìn)行了檢驗(yàn)并與NOAH(National Center for Environmental Prediction, Oregon State University, Air Force, and Hydrology Lab model)陸面過(guò)程模式的模擬結(jié)果進(jìn)行了比較，結(jié)果表明PKULM能夠較好地模擬西北半干旱區(qū)農(nóng)田下墊面地氣交換過(guò)程.

關(guān)鍵詞陸面過(guò)程模式； 地表能量平衡； 土壤濕度； 氣孔導(dǎo)度

The model included five modules: radiation transfer, turbulence, photosynthesis, soil heat diffusion and soil water transport. PKULM has augmented biophysical and hydrological processes over SPAM and introduced a new numerical solution for coupled nonlinear equation set of land surface processes: (1) Radiation transfer process within canopy is introduced into the model. The two-stream approximation is used to simulate this process. (2) We use the photosynthesis-based Ball-Berry stomatal resistance scheme instead of the Jarvis algorithm for modeling transpiration, which controls canopy carbon and water fluxes consistently and simultaneously. (3) The original algorithm used in Community Land Model (CLM) for solving coupled photosynthesis and stomatal resistance model does not converge in a dry environment. We propose a new and efficient bisection method that can ensure convergence with ambient low environmental vapor pressure. (4) We update the numerical solution of Richards equation from a saturation-based form to a head-based form. It is apt for modeling soil water movement either in the vadose zone or saturated zone. (5) A new mechanism runoff model is adopted to replace the parameterized model. (6) We adapt the traditional Thomas algorithm and use an implicit iterative method for solving the coupled energy balance and water balance equations of soil, canopy and land surface. With the data obtained in NWC-ALIEX (Northwest China Atmosphere-Land Interaction Experiment), we evaluated the performance of PKULM over cropland in Pingliang and compared it with the Noah land surface model.

PKULM outperforms Noah in modeling reflected solar radiation (SR), emitted longwave radiation (LR), sensible heat flux (SH) and latent heat flux (LH) in terms of lower root mean square error (RMSE) and higher correlation coefficients (CC). The RMSEs of PKULM estimations are 4.1, 6.6, 16.3, and 21.8 W·m-2for SR, LR, SH, and LH, respectively. The CCs of PKULM modeled SR, LR, SH, and LH, are 0.999, 0.996, 0.976, and 0.989, respectively. Noah produced top 10 cm soil temperature (ST) with a cold bias of 1.03 K. PKULM performs better than Noah in estimating ST in aspects of reproducing smaller bias (0.289 K), more accurate phase, and amplitude. The RMSE of modeled ST from PKULM in top 5 cm is 0.73 K and CC is 0.978. Both Noah and PKULM reproduced the descending trend of soil moisture in the simulation period, but overestimated the magnitude in the daytime and underestimated it in the nighttime.

We developed PKULM on the basis of SPAM. PKULM has augmented physical processes and numerical solution that are adapted for arid and semi-arid regions. The modeled energy and water fluxes between land and atmosphere are evaluated and compared to observations. The results show PKULM could reproduce better estimations of SR, LR, SH, and LH than Noah. However, the estimated soil moisture could be improved in the future by adopting more accurate numerical solutions of Richards equation and utilizing more realistic parameterizations of hydrological processes such as surface runoff, evaporation, and drainage.

1引言

陸地表面約占地球表面積的29.2%，陸面過(guò)程是影響地球系統(tǒng)能量水分循環(huán)的重要物理過(guò)程.一方面，作為大氣圈的下墊面，陸地表面控制著太陽(yáng)凈輻射在長(zhǎng)波輻射、感熱和潛熱之間的分配比例，通過(guò)湍流輸送過(guò)程加熱大氣、驅(qū)動(dòng)大氣運(yùn)動(dòng)，影響區(qū)域和全球氣候(周連童和黃榮輝, 2008; Yang et al., 2011)；另一方面，作為氣候系統(tǒng)中能量和水分的重要儲(chǔ)庫(kù)，土壤溫度和濕度與大氣環(huán)流有長(zhǎng)周期的相互作用(Koster et al., 2004)；同時(shí)，作為生物活動(dòng)的主要場(chǎng)所，陸地表面的生物物理、化學(xué)特性還對(duì)陸氣交換過(guò)程乃至未來(lái)氣候變化評(píng)估產(chǎn)生著重要的影響(Feddema et al., 2005; Davidson et al., 2006; Jiang et al., 2011).

發(fā)展陸面過(guò)程模式，準(zhǔn)確模擬陸氣間能量物質(zhì)交換過(guò)程，是提高氣候系統(tǒng)可預(yù)報(bào)性的重要途徑(Sellers et al., 1997; Pitman, 2003).Manabe(1969)發(fā)展的“水桶”模型，以地表能量平衡和水分平衡方程作為核心，是建立陸面過(guò)程模式的第一次嘗試.Deardorff(1978)在他的模式中使用了近地面湍流相似性理論計(jì)算陸氣間感熱和潛熱通量，并第一次考慮了植被冠層作用，其工作在陸面過(guò)程模式的發(fā)展過(guò)程中具有劃時(shí)代的意義.這種考慮了植被冠層作用的陸面過(guò)程模式也被稱(chēng)為第二代陸面過(guò)程模式；BATS(Dickinson, 1984)、SiB(Sellers et al., 1986)、IAP94(Dai and Zeng, 1997)是其中的典型代表.20世紀(jì)末，考慮植被光合作用的SiB2(Sellers et al., 1996)、CLM(Bonan et al., 2002)、CoLM(Dai et al., 2003)等一大批第三代陸面過(guò)程模式相繼出現(xiàn).第三代陸面過(guò)程模式以詳細(xì)考慮光合作用為主要特點(diǎn)，綜合考慮了植被冠層內(nèi)輻射傳輸，土壤水熱輸送、植被生化作用、近地面湍流輸送等過(guò)程，具有較高的模擬精度，是當(dāng)前陸面過(guò)程模式發(fā)展的前沿方向之一.

Liu(1996)在Noilhan和Planton(1989)的N89模式的基礎(chǔ)上，發(fā)展了一個(gè)包含植被遮蔽和蒸騰作用的陸面水熱交換模式(LSEM, Land Surface Exchange Model).該模式很好地模擬森林和沙漠綠洲下墊面溫度變化及能量交換特征.劉和平等(1999)引入Darcy定律描述土壤水分輸送過(guò)程，并將模式擴(kuò)展到草原、戈壁及熱帶森林等多種下墊面，經(jīng)過(guò)實(shí)測(cè)資料驗(yàn)證，證明了模式能夠模擬多種復(fù)雜下墊面地氣間能量、水分交換過(guò)程.Liu等(2004a, b)將陸面過(guò)程模式與北京大學(xué)大氣邊界層模式相互耦合，發(fā)展了土壤—植被—大氣耦合模式(SPAM，Soil-Plant-Atmosphere Model)，準(zhǔn)確模擬了沙漠—綠洲系統(tǒng)的“冷島效應(yīng)”和“濕島效應(yīng)”.Liu等(2006)改進(jìn)了輻射傳輸過(guò)程，采用了更加適合干旱區(qū)半干旱區(qū)的通量廓線關(guān)系，模擬的地表溫度有了大幅度的改進(jìn).鄭輝和劉樹(shù)華(2012, 2013)和劉樹(shù)華等(2013)從土壤熱傳導(dǎo)、湍流輸送、植物生理生化過(guò)程等多個(gè)方面進(jìn)一步改進(jìn)模式，經(jīng)過(guò)實(shí)測(cè)資料驗(yàn)證，取得了較好的模擬結(jié)果.

本文在已有工作的基礎(chǔ)上，考慮了植被以及地表水分平衡過(guò)程，建立了北京大學(xué)陸面過(guò)程模式PKULM.本文首先介紹PKULM的模式結(jié)構(gòu)、參數(shù)化方案以及計(jì)算方法，然后利用西北半干旱區(qū)農(nóng)田下墊面觀測(cè)資料驅(qū)動(dòng)模式，對(duì)模式進(jìn)行了檢驗(yàn)，評(píng)估模式的模擬能力.

圖1　北京大學(xué)陸面過(guò)程模式示意圖

2模式介紹

如圖1所示，北京大學(xué)陸面過(guò)程模式包含如下五個(gè)相互耦合的過(guò)程：(1)輻射傳輸過(guò)程，確定不同波段輻射能量在植被冠層和土壤表面之間的分配，為光合作用、感熱輸送和潛熱輸送提供能量來(lái)源；(2)湍流輸送過(guò)程，確定熱量、水分和二氧化碳等的湍流輸送通量及邊界層大氣運(yùn)動(dòng)的摩擦阻力；(3)土壤熱力學(xué)過(guò)程，確定土壤中水分、熱量輸送及其與大氣相互作用的物理過(guò)程；(4)植物生理過(guò)程，確定植物蒸騰速率及光合作用強(qiáng)度的生物理化過(guò)程；(5)水分平衡過(guò)程，確定降水量的分配，估算土壤含水量、徑流量、蒸發(fā)量、蒸騰量.目前，模式不包含積雪與凍土過(guò)程.

在PKULM中，每一個(gè)格點(diǎn)上可以有m種不同的植被冠層覆蓋，分別計(jì)算每一種冠層所覆蓋地表的短波輻射通量、長(zhǎng)波輻射通量、感熱通量、水汽通量以及二氧化碳通量，然后以相應(yīng)的植被覆蓋率σj為權(quán)重加總，即可得到陸面與大氣之間總的能量、物質(zhì)交換通量.

1.1輻射傳輸過(guò)程

1.2.1太陽(yáng)輻射

由于冠層的光學(xué)特性在不同波段存在顯著差異，PKULM將太陽(yáng)輻射分為0～700 nm和700～2800 nm兩個(gè)波段Λ分別計(jì)算；在每一個(gè)波段，還區(qū)分對(duì)待冠層頂所接收到的直射和散射太陽(yáng)光.

在植被覆蓋的地區(qū)，太陽(yáng)輻射不斷在植物葉片間反射、透射和吸收，最終光強(qiáng)趨于穩(wěn)態(tài)分布，該分布可以用二流模型近似為(Dickinson, 1984)：

(1)

(2)

(3)

IΛ=1-I↑Λ-(1-αg,Λ)I↓Λ，

(4)

冠層和地表吸收的凈輻射通量Sv,net和Sg,net(W·m-2)分別為：

(5)

(6)

2.2長(zhǎng)波輻射

圖2是PKULM長(zhǎng)波輻射傳輸示意圖.Latm↓是大氣向下的長(zhǎng)波輻射強(qiáng)迫(W·m-2).根據(jù)斯特芬-波爾茲曼定律以及比輻射率的定義，植被冠層向下的長(zhǎng)波輻射通量Lv↓(W·m-2)等于Latm↓透過(guò)冠層的部分與冠層熱發(fā)射值之和，即：

(7)

圖2　PKULM長(zhǎng)波輻射傳輸示意圖

其中，εv是冠層的比輻射率，Tv是冠層溫度(K).同理，地表向上和冠層向上的長(zhǎng)波輻射通量Lg↑和L↑(W·m-2)分別為：

(8)

(9)

其中，εg是地表的比輻射率，Tg是地表溫度(K).因此，冠層和地表吸收的凈長(zhǎng)波輻射Lv,net、Lg,net分別為：

(10)

(11)

1.2湍流通量

1.2.3植被冠層感熱、潛熱通量

PKULM的湍流輸送過(guò)程采用阻抗模型計(jì)算.土壤表面、冠層以及大氣之間的感熱、潛熱輸送通道如圖3所示，其湍流阻抗采用Monin-Obukhov相似性理論計(jì)算(Zeng et al., 1998).

假設(shè)冠層間空氣的熱量存儲(chǔ)可以忽略，陸地表面輸送至大氣的感熱通量H(W·m-2)等于冠層輸送至冠層間空氣的感熱通量Hv(W·m-2)與地表輸送至冠層間空氣的感熱通量Hg之和為(W·m-2)：

(12)

(13)

(14)

(15)

圖3　地氣間感熱、潛熱輸送通道及其阻抗示意圖

聯(lián)立式(12)至(15)，可以求出冠層間空氣溫度Ts的表達(dá)式為

(16)

同樣陸地表面輸送至大氣的水汽通量E(kg·m-2·s-1)等于冠層輸送至冠層間空氣的水汽通量Ev(kg·m-2·s-1)與地表輸送至冠層間空氣的水汽通量Eg(kg·m-2·s-1)之和為：

(17)

(18)

(19)

(20)

聯(lián)立式(17)至(20)，可以求出冠層間空氣比濕qs的表達(dá)式為

(21)

1.3光合作用及氣孔導(dǎo)度

由于植物葉片覆蓋了物質(zhì)滲透率極低的角質(zhì)層，氣孔成為了植被與大氣之間進(jìn)行物質(zhì)交換的主要通道，氣孔的開(kāi)閉大小直接影響了植被冠層蒸騰作用的強(qiáng)度.第三代陸面過(guò)程模型均采用Ball-Berry半經(jīng)驗(yàn)?zāi)Ｐ?Collatz et al., 1991)計(jì)算氣孔阻抗rs(s·m-1)：

(22)

其中m、b為經(jīng)驗(yàn)系數(shù)，Patm為大氣壓(Pa)，cs為葉片表面二氧化碳分壓(Pa)，es為葉片表面水汽分壓(Pa)，ei為氣孔內(nèi)水汽分壓(Pa).二氧化碳同化速率A(μmol·m-2·s-1)采用Collatz 等(1991, 1992)方案計(jì)算，它是冠層溫度、土壤水勢(shì)、光合作用有效輻射以及二氧化碳、水氣分壓的函數(shù).

根據(jù)阻抗模型，二氧化碳同化速率A，氣孔內(nèi)二氧化碳和水汽分壓ci、ei，以及大氣中二氧化碳和水汽分壓ca、ea有如下關(guān)系：

(23)

2.4二分查找算法

為了求解氣孔阻抗，令函數(shù)h(ci)=ca-(1.37rb+1.65rs)PatmA.由于氣孔阻抗rs和光合作用強(qiáng)度A都是ci的函數(shù)，因此h(ci)也是ci的函數(shù).利用h(ci)，式(23)等價(jià)于如下的定點(diǎn)問(wèn)題：

(24)

圖4　不同水汽壓下PKULM的氣孔導(dǎo)度定點(diǎn)函數(shù)參數(shù)設(shè)置：ca = 28 Pa, t=25 ℃, rb=1×10-6 s·m2 μmol-1, vcmax = 40 μmol CO2 m-2·s-1.

圖5　計(jì)算氣孔內(nèi)二氧化碳分壓的二分查找法參數(shù)設(shè)置：ca=28 Pa, ea=634 Pa, t=25 ℃, rb=1×10-6 s·m2 μmol-1, vcmax=40 μmol CO2 m-2·s-1.

2.4冠層及地表能量平衡

PKULM每一個(gè)格點(diǎn)上可以有m種不同的植被冠層覆蓋，對(duì)于第j種植被冠層，根據(jù)能量平衡，冠層吸收的太陽(yáng)凈輻射和凈長(zhǎng)波輻射全部通過(guò)湍流輸送過(guò)程加熱大氣和地表，即：

(25)

對(duì)于地表，其吸收的太陽(yáng)凈輻射和凈長(zhǎng)波輻射等于地表向上的感熱、潛熱通量與向下的土壤熱通量之和，即：

(26)

可以將上述兩式寫(xiě)為

(27)

(28)

2.5土壤熱量輸送

土壤中熱量的存儲(chǔ)和釋放過(guò)程可以用熱力學(xué)第一定律和梯度熱傳導(dǎo)方程描述為：

(29)

(30)

其中，t為時(shí)間(s)，z為土壤深度(m)，T為土壤溫度(K)，F(xiàn)為土壤熱通量(W·m-2，向上為正)，C為土壤體積熱容(J·m-3·K-1)，kT為土壤熱傳導(dǎo)系數(shù)(W·m-1·K-1).

把土壤分為n層，采用隱式差分方案，將上述兩式離散化，則第i層土壤溫度的預(yù)報(bào)方程為：

(31)

(32)

(33)

可以將式(31)寫(xiě)成三對(duì)角矩陣形式，以便進(jìn)一步數(shù)值求解，公式為

(34)

2.6冠層降水截留

大氣降水qrain(kg·m-2·s-1)在經(jīng)過(guò)植被冠層時(shí)，有一部分直接降落至地面qthru(kg·m-2·s-1)，而另一部分被截留qintr(kg·m-2·s-1)，然后部分滴落至地表qdrip(kg·m-2·s-1)，沒(méi)有滴落的部分依附在葉莖表面，成為冠層蓄水wcan(kg·m-2).冠層截流降水量qintr、直接落至地表的降水量qthru以及冠層滴落水量qdrip分別為(Lawrence et al., 2007)：

(35)

(36)

(37)

其中pwcan是前一時(shí)間步上冠層蓄水量(kg·m-2)，wcan,max=0.1(L+S)是最大冠層蓄水量(kg·m-2)，L和S分別為葉面積指數(shù)和莖面積指數(shù)(m2·m-2).

最終落在地表的總水量qgrnd(kg·m-2·s-1)為

(38)

2.7地表徑流及滲流

PKULM發(fā)展了一個(gè)地表徑流及滲流的機(jī)理模型.在降水過(guò)程中，地表從水分不飽和狀態(tài)到產(chǎn)生徑流將經(jīng)歷如下三個(gè)階段：

(1)土壤非飽和階段.該階段，表層土壤毛細(xì)水勢(shì)小于其飽和水勢(shì)ψ1<ψsat,1，根據(jù)地表水分平衡，降落在地表的總水量qgrnd等于地表蒸發(fā)量Eg與滲流量qinfl(kg·m-2)之和，公式為

(39)

(2)土壤飽和、地表出現(xiàn)積水階段.該階段，土壤毛細(xì)水勢(shì)大于其飽和水勢(shì)，且小于地表最大積水深度wg,max=0.2×10-3m(Sellers et al., 1996)；降落在地表的總水量與地表蒸發(fā)量、地表滲流量之差等于積水的增加量，公式為

ψ1≥ψsat,1,wg

(40)

(3)徑流產(chǎn)生階段.該階段，土壤毛細(xì)水勢(shì)等于地表最大積水深度，土壤飽和，積水深度保持wg,max不變，降落在地表的總水量與地表蒸發(fā)量、地表滲流量之差等于徑流，公式為

(41)

上述三式中，(Eg)j是第j種植被冠層所覆蓋地表的蒸發(fā)量(kg·m-2·s-1)，ρliq是水的密度(kg·m-3)，qinfl是土壤表面滲流(kg·m-2·s-1)，其大小與第一層土壤導(dǎo)水率kh,1(kg·m-2·s-1)，土壤毛細(xì)水勢(shì)ψ1(m)以及深度z1(m)有關(guān)，公式為

(42)

值得注意的是，我們采用統(tǒng)一的變量wg表示表層土壤毛細(xì)水勢(shì)以及地表積水深度：當(dāng)wg≤0時(shí)，wg表示表層土壤毛細(xì)水勢(shì)(m)，而當(dāng)wg>0時(shí)，wg表示地表積水深度(m).使用這種表示方法，再利用適當(dāng)?shù)妮o助函數(shù)hg，我們可以將不同階段的三個(gè)地表水分平衡方程式(39)、式(40)和式(41)寫(xiě)成統(tǒng)一的差分形式為：

(43)

(44)

其中pwg表示前一時(shí)間步上wg的值.

上式可進(jìn)一步寫(xiě)為

(45)

其中：

(46)

(47)

(48)

2.8土壤濕度

土壤水分運(yùn)動(dòng)可以采用Richards方程描述為(Zeng and Decker, 2009)：

(49)

(50)

其中，w為土壤體積含水量(m·m-3)，q為土壤水分通量(kg·m-2·s-1)，kh為土壤導(dǎo)水率(kg·m-2·s-1)，ψ為土壤毛細(xì)水勢(shì)(m)，qroot為植物根系吸水量(kg·m-1·s-1).

為了數(shù)值求解上述兩個(gè)方程，需要利用土壤含水量與土壤毛細(xì)水勢(shì)的函數(shù)關(guān)系，將方程轉(zhuǎn)化成只含毛細(xì)水勢(shì)或只含土壤含水量的偏微分方程(Ross, 2003)：

(51)

其中，wsat是土壤飽和含水量(m·m-3)，ψsat為土壤飽和水勢(shì)(m)，B為Clapp-Hornberger參數(shù).

CLM、NOAH等陸面過(guò)程模式一般將土壤含水量作為預(yù)報(bào)變量，即模擬土壤含水量變化，然后利用式(51)計(jì)算土壤毛細(xì)水勢(shì).但Ross(2003)指出，由于土壤含水量與土壤毛細(xì)水勢(shì)的一一映射關(guān)系只在非飽和時(shí)成立，因此這種以土壤含水量為基礎(chǔ)的水分運(yùn)動(dòng)模式僅能模擬土壤非飽和區(qū)(地下水位以上)的水分運(yùn)動(dòng).如圖6所示，在土壤飽和區(qū)，土壤含水量維持在飽和含水量不變，根據(jù)式(53)，土壤毛細(xì)水勢(shì)保持不變；但實(shí)際中，毛細(xì)水勢(shì)隨深度增加而增加.

圖6　土壤毛細(xì)水勢(shì)和土壤體積含水量廓線參數(shù)設(shè)置：飽和含水量為0.3731 m3·m-3，飽和土壤毛細(xì)水勢(shì)為-0.0473 m，Clapp-Hornbeger參數(shù)為 3.387，地下水位為 -1.737 m.

PKULM采用土壤水勢(shì)為基礎(chǔ)的Richards方程.使用隱式差分方案，將式(49)和式(50)離散化，得到第i層土壤毛細(xì)水勢(shì)的預(yù)報(bào)方程為

(52)

其中，pψi為前一時(shí)間步第i層土壤的毛細(xì)水勢(shì)(m)；kh[zh,i]為第i層與第i+1層土壤界面處的導(dǎo)水率，公式為

(53)

式(52)和式(53)可以進(jìn)一步寫(xiě)成三對(duì)角矩陣形式為

(54)

其中：

(55)

(56)

(57)

(58)

3控制方程及隱式Picard迭代計(jì)算算法

PKULM的核心是能量與水分平衡控制方程組.冠層能量平衡方程式(27)、地表能量平衡方程式(28)以及土壤熱擴(kuò)散方程式(34)構(gòu)成了PKULM的能量平衡控制方程組.方程組可以寫(xiě)成如下的矩陣形式為

(59)

地表水分平衡方程式(53)以及Richards方程式(54)構(gòu)成了PKULM的水分平衡控制方程組.方程組可以寫(xiě)成如下的矩陣形式為

(60)

為了方便討論，將式(59)和式(60)記作Aφ=B，其中φ是全部狀態(tài)參量所構(gòu)成的向量，A是全部系數(shù)構(gòu)成的矩陣，B是全部常數(shù)項(xiàng)構(gòu)成的矩陣.由于矩陣A、B都是關(guān)于φ的非線性函數(shù)，我們采用Picard迭代法(Hipsey et al., 2004)將方程線性化并求解，具體做法為：(1)使用上一時(shí)刻狀態(tài)量的值pφ作為φ的初值；(2)將上一迭代步狀態(tài)量的值φl(shuí)，帶入系數(shù)矩陣，求解線性方程組A(φl(shuí))φl(shuí)+1=B(φl(shuí))，以方程的解φl(shuí)+1作為當(dāng)前迭代步狀態(tài)量φ的值；(3)經(jīng)過(guò)若干步迭代，當(dāng)φl(shuí)+1≈φl(shuí)時(shí)，可認(rèn)為φl(shuí)+1就是非線性方程組A(φ)φ=B(φ)的解.

對(duì)于形如式(60)的線性方程組，可以直接采用Thomas算法快速求解.而對(duì)于形如式(59)的線性方程組，需要對(duì)Thomas算法進(jìn)行修改，具體算法如下：(1)依次用矩陣第j(1≤j≤m)行消去第m+1行的各項(xiàng)系數(shù)hv,j(1≤j≤m)；(2)采用標(biāo)準(zhǔn)Thomas算法求解方程組第m+1行至第m+n+1行；(3)最后，用矩陣第m+1行依次消去第m+1列各項(xiàng)系數(shù)vv,j(1≤j≤m).

4觀測(cè)資料介紹

本文使用“我國(guó)西北典型干旱半干旱區(qū)能量和水分循環(huán)觀測(cè)試驗(yàn)與分析”研究項(xiàng)目(黃榮輝等, 2013)平?jīng)鲇^測(cè)站2010年加強(qiáng)觀測(cè)期資料作為模式輸入，進(jìn)行模擬試驗(yàn).平?jīng)稣疚挥诟拭C省平?jīng)鍪?35°35′N(xiāo)，106°42′E)，海拔為1417m.地表平坦均一，土壤為黃土，下墊面作物為玉米.平?jīng)稣镜闹饕^測(cè)儀器包括：四分量輻射計(jì)(CNR1，Kipp&Zonen)、三維超聲風(fēng)速儀(CSAT3，Campell)、溫度和相對(duì)濕度傳感器(HMP45C，Vaisala)、大氣壓傳感器(CS105，Campbell)、CO2/H2O開(kāi)路氣體分析儀(LI-7500，LI-COR)、4層(5，10，20，40cm)土壤溫度測(cè)量系統(tǒng)(CS616-L，Campbell)和4層(5cm，10，20，40cm)土壤濕度測(cè)量系統(tǒng)(109-L，Campbell).2010年8月14日19時(shí)至2010年8月16日09時(shí)，無(wú)天氣系統(tǒng)過(guò)境，觀測(cè)記錄完整，且根據(jù)穩(wěn)態(tài)測(cè)試(SST,SteadyStateTest)和整體湍流特征檢驗(yàn)(ITC,IntegralTurbulenceCharacteristics)對(duì)湍流資料進(jìn)行了質(zhì)量控制(Foken, 2008).

5模擬結(jié)果

模式以2010年8月14日19時(shí)觀測(cè)得到的各層土壤溫度和濕度資料作為初始狀態(tài)，以平均風(fēng)速、平均氣溫、相對(duì)濕度、大氣壓、二氧化碳分壓以及向下的長(zhǎng)短波輻射作為輸入量，以半小時(shí)為步長(zhǎng)向前積分，從2010年8月15日06時(shí)起輸出模擬結(jié)果.同時(shí)使用相同的大氣強(qiáng)迫資料和初始土壤溫度、濕度場(chǎng)驅(qū)動(dòng)NOAH陸面過(guò)程模式，將模擬結(jié)果與PKULM進(jìn)行比較.

分別使用平均偏差(MeanBiasError,MBE)、均方根誤差(RootMeanSquareError,RMSE)和相關(guān)系數(shù)(CorrelationCoefficient,CC)分析模擬誤差，定義如下：

(61)

(62)

(63)

表1　PKULM和NOAH模擬能量通量的平均誤差(MBE)、標(biāo)準(zhǔn)差(RMSE)和相關(guān)系數(shù)(CC)

圖7是PKULM和NOAH模擬的反射短波輻射及其與觀測(cè)值的比較.中午12時(shí)，平?jīng)稣居^測(cè)得到的反射短波輻射達(dá)到最大值，約為140 W·m-2，此時(shí)PKULM模擬的結(jié)果略高于觀測(cè)值0.39 W·m-2，而NOAH模擬結(jié)果則低于觀測(cè)值約4.0 W·m-2.總體而言，兩個(gè)模式均能較好地模擬短波輻射傳輸過(guò)程，模擬的反射短波輻射均方根誤差分別為4.061 W·m-2和11.34 W·m-2.

圖8是PKULM和NOAH模擬的地表發(fā)射長(zhǎng)波輻射及其與觀測(cè)值的比較.可以看到，PKULM的模擬結(jié)果與觀測(cè)值接近，模擬均方根誤差為6.608 W·m-2.NOAH的模擬結(jié)果偏低，這可能與其感熱通量模擬偏高(圖9)，導(dǎo)致地表溫度模擬結(jié)果偏低(圖11)造成的.

圖9和圖10分別給出了感熱通量和潛熱通量的模擬值與觀測(cè)值.從圖中可以看到，由于植被蒸騰作用，平?jīng)鲛r(nóng)田下墊面上潛熱通量占主導(dǎo)地位，潛熱通量大小約為感熱通量的兩倍.兩個(gè)模式均能比較好地模擬潛熱通量大小及日變化過(guò)程，PKULM和NOAH模擬值的均方根誤差分別為21.81 W·m-2和24.05 W·m-2.而對(duì)感熱通量而言，PKULM的模擬結(jié)果更接近觀測(cè)值，均方根誤差為16.29 W·m-2，而NOAH模擬結(jié)果的均方根誤差為29.95 W·m-2.

圖11和圖12分別是5 cm深度土壤溫度、濕度的模擬值及其與觀測(cè)值的比較.PKULM的模擬值均比NOAH的模擬值更接近實(shí)測(cè)值.對(duì)于土壤溫度，PKULM模擬的溫度日變化振幅與相位均與實(shí)測(cè)值相符，平均偏差為0.289 K，均方根誤差為0.732 K，相關(guān)系數(shù)為0.978，而NOAH模擬的土壤溫度偏低，平均偏差為-1.03 K，均方根誤差為1.79 K，相關(guān)系數(shù)為0.903.對(duì)于土壤濕度，PKULM模擬結(jié)果的平均偏差為-7.4×10-3m3·m-3，均方根誤差為7.5×10-3m3·m-3，相關(guān)系數(shù)為0.906；NOAH模擬結(jié)果的平均偏差為-2.7×10-2m3·m-3，均方根誤差為3.0×10-2m3·m-3，相關(guān)系數(shù)為0.792.兩個(gè)模式模擬的土壤濕度變化趨勢(shì)與實(shí)際一致，但大小與實(shí)測(cè)有一定偏差.兩個(gè)模式模擬的白天土壤濕度降低速度高于實(shí)測(cè)，而夜晚低于實(shí)測(cè)，模式對(duì)土壤水分運(yùn)動(dòng)的參數(shù)化方案仍然有待改進(jìn).

圖7　向上短波輻射的觀測(cè)值與PKULM和NOAH模擬值的比較圖

圖8　向上長(zhǎng)波輻射的觀測(cè)值與PKULM和NOAH模擬值的比較圖

圖9　感熱通量的觀測(cè)值與PKULM和NOAH模擬值的比較圖

圖10　潛熱通量的觀測(cè)值與PKULM和NOAH模擬值的比較圖

圖11　5 cm深度土壤溫度觀測(cè)值與PKULM和NOAH模擬值的比較圖

圖12　5 cm深度土壤濕度觀測(cè)值與PKULM和NOAH模擬值比較圖

6結(jié)論

本文介紹了北京大學(xué)陸面過(guò)程模式PKULM的參數(shù)化方案和計(jì)算算法，并利用平?jīng)鲛r(nóng)田下墊面觀測(cè)資料對(duì)模式進(jìn)行了檢驗(yàn).主要結(jié)論如下：

(1) 通過(guò)使用二分查找算法計(jì)算氣孔導(dǎo)度，避免了CLM等模式使用的迭代算法在低水汽壓環(huán)境下不穩(wěn)定的問(wèn)題，提高了模式在干旱和半干旱地區(qū)的適用性.

(2) 本文發(fā)展了一套求解陸面過(guò)程非線性方程組的Picard隱式迭代算法，該算法計(jì)算穩(wěn)定，且有效地將各物理過(guò)程參數(shù)化方案與模式總體框架隔離開(kāi)來(lái)，有利于模式的進(jìn)一步發(fā)展.

(3) PKULM能夠模擬半干旱區(qū)農(nóng)田下墊面陸氣交換過(guò)程.模式模擬的短波輻射、長(zhǎng)波輻射、感熱和潛熱通量與實(shí)測(cè)值吻合，模擬誤差略低于NOAH陸面過(guò)程模式.

(4) 模式對(duì)土壤濕度的模擬誤差較大，與土壤濕度密切相關(guān)的根系過(guò)程、土壤水分輸送過(guò)程仍然需要進(jìn)一步改進(jìn).

References

Bonan G B, Oleson K W, Vertenstein M, et al. 2002. The land surface climatology of the community land model coupled to the NCAR community climate model.J.Climate, 15(22): 3123-3149, doi: 10.1175/1520-0442(2002)015〈3123:TLSCOT〉2.0.CO;2.

Collatz G J, Ball J T, Grivet C, et al. 1991. Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer.Agr.ForestMeteor., 54(2-4): 107-136, doi: 10.1016/0168-1923(91)90002-8.Collatz G J, Ribas-Carbo M, Berry J A. 1992. Coupled photosynthesis-stomatal conductance model for leaves of C4 plants.Aust.J.PlantPhysiol., 19(5): 519-538.

Dai Y J, Zeng Q C. 1997. A land surface model (IAP94) for climate studies part I: Formulation and validation in off-line experiments.Adv.Atmos.Sci., 14(4): 433-460, doi: 10. 1007/s00376-997-0063-4.

Dai Y J, Zeng X B, Dickinson R E, et al. 2003. The Common Land Model.Bull.Amer.Meteor.Soc., 84(8): 1013-1023, doi: 10.1175/BAMS-84-8-1013.

Davidson E A, Janssens I A. 2006. Temperature sensitivity of soil carbon decomposition and feedbacks to climate change.Nature, 440(7081): 165-173, doi: 10.1038/nature04514.

Deardorff J W. 1978. Efficient prediction of ground surface temperature and moisture, with inclusion of a layer of vegetation.J.Geophys.Res.-Atmos., 83(C4): 1889-1903, doi: 10.1029/JC083iC04p01889.

Dickinson R E. 1984. Modeling evapotranspiration for three-dimensional global climate models. ∥ Hansen J E, Takahashi T, eds. Climate Processes and Climate Sensitivity. Washington D. C.: American Geophysical Union, 58-72.

Feddema J J, Oleson K W, Bonan G B, et al. 2005. The importance of land-cover change in simulating future climates.Science, 310(5754): 1674-1678.

Foken T. 2008. Micrometeorology. Berlin: Springer, 1-308.

Hipsey M R, Sivapalan M, Clement T P. 2004. A numerical and field investigation of surface heat fluxes from small wind-sheltered waterbodies in semi-arid Western Australia.Environ.FluidMech., 4(1): 79-106, doi: 10.1023/A:1025547707198.

Huang R H, Zhou D G, Chen W, et al. 2013. Recent progress in studies of air-land interaction over the arid area of Northwest China and its impact on climate.ChineseJournalofAtmosphericSciences(in Chinese), 37(2): 189-210, doi: 10.3878/j.issn.1006-9895.2012.12303.

Jiang D B, Zhang Y, Lang X M. 2011. Vegetation feedback under future global warming.TheoreticalandAppliedClimatology, 106(1-2): 211-227.

Koster R D, Dirmeyer P A, Guo Z C, et al. 2004. Regions of strong coupling between soil moisture and precipitation.Science, 305(5687): 1138-1140, doi: 10.1126/science.1100217.

Lawrence D M, Thornton P E, Oleson K W, et al. 2007. The partitioning of evapotranspiration into transpiration, soil evaporation, and canopy evaporation in a GCM: Impacts on land-atmosphere interaction.J.Hydrometeor., 8(4): 862-880.

Liu H P, Liu S H, Sang J G. 1999. A modified SiB to simulate momentum, heat and water transfer over various underlying surfaces.ChineseJournalofAtmosphericSciences(in Chinese), 23(4): 449-460. Liu S H. 1996. A parameterized model for the influence of vegetation on the moisture-heat exchange at near-ground layer.ActaSci.Nat.Univ.Pek., 32(1): 69-77.

Liu S H, Liu Z X, Zheng H, et al. 2013. Multi-scale atmospheric boundary-layer and land surface physics process models.Sci.Sin.Phys.Mech.Astron. (in Chinese), 43(10): 1332-1355.

Liu S H, Pan Y, Deng Y, et al. 2006. Numerical simulation experiment of land surface physical processes and local climate effect in forest underlying surface.Acta.Meteor.Sin., 20(1): 72-85.

Liu S H, Yue X, Hu F, et al. 2004a. Using a Modified Soil-Plant-Atmosphere Scheme (MSPAS) to simulate the interaction between land surface processes and atmospheric boundary layer in semi-arid regions.Adv.Atmos.Sci., 21(2): 245-259, doi: 10.1007/BF02915711.

Liu S H, Yue X, Liu H Z, et al. 2004b. Using a Modified Soil-Plant-Atmosphere Scheme (MSPAS) to study the sensitivity of land surface and boundary layer processes to soil and vegetation conditions.Adv.Atmos.Sci., 21(5): 717-729, doi: 10.1007/BF02916369.Manabe S. 1969. Climate and the ocean circulation 1. The atmospheric circulation and the hydrology of the earth′s surface.Mon.Wea.Rev., 97(1): 739-774.

Noilhan J, Planton S. 1989. A simple parameterization of land surface processes for meteorological models.Mon.Wea.Rev., 117(3): 536-549, doi: 10.1175/1520-0493(1989)117〈0536:ASPOLS〉2.0.CO;2.

Pitman A J. 2003. The evolution of, and revolution in, land surface schemes designed for climate models.Int.J.Climatol., 23(5): 479-510, doi: 10.1002/joc.893.

Ross P J. 2003. Modeling soil water and solute transport-fast, simplified numerical solutions.Agron.J., 95(6): 1352-1361.

Sakaguchi K, Zeng X B. 2009. Effects of soil wetness, plant litter, and under-canopy atmospheric stability on ground evaporation in the Community Land Model (CLM3. 5).J.Geophys.Res., 114(D1): D01107, doi: 10.1029/2008JD010834.

Sellers P J, Dickinson R E, Randall D A, et al. 1997. Modeling the exchanges of energy, water, and carbon between continents and the atmosphere.Science, 275(5299): 502-509, doi: 10.1126/science.275.5299.502.

Sellers P J, Mintz Y, Sud Y C, et al. 1986. A Simple Biosphere Model (SiB) for use within general circulation models.J.Atmos.Sci., 43(6): 505-531, doi: 10.1175/1520-0469(1986)043〈0505:ASBMFU〉2.0.CO;2.

Sellers P J, Randall D A, Collatz G J, et al. 1996. A revised land surface parameterization (SiB2) for atmospheric GCMs. Part I: Model formulation.J.Climate, 9(4): 676-705, doi: 10.1175/1520-0442(1996)009〈0676:ARLSPF〉2.0.CO;2.

Sun Y, Gu L H, Dickinson R E. 2012. A numerical issue in calculating the coupled carbon and water fluxes in a climate model.J.Geophys.Res.Atmos., 117(D22): D22103, doi: 10.1029/2012JD018059.

Yang T, Wang X Y, Zhao C Y. 2011. Changes of climate extremes in a typical arid zone: Observations and multimodel ensemble projections.J.Geophys.Res.Atmos., 116(D19): D19106, doi: 10.1029/2010JD015192.

Zeng X B, Decker M. 2009. Improving the numerical solution of soil moisture-based Richards equation for land models with a deep or shallow water table.J.Hydrometeorol., 10(1): 308-319,doi: 10.1175/2008JHM1011.1.

Zeng X B, Zhao M, Dickinson R E. 1998. Intercomparison of bulk aerodynamic algorithms for the computation of sea surface fluxes using TOGA COARE and TAO data.J.Climate, 11(10): 2628-2644.

Zheng H, Liu S H. 2012. Effects of difference and grid scheme in soil temperature simulation.ChineseJ.Geophys. (in Chinese), 55(8): 2514-2522, doi: 10.6038/j.issn.0001-5733.2012.08.004.Zheng H, Liu S H. 2013. Land surface parameterization and modeling over desert.ChineseJ.Geophys. (in Chinese), 56(7): 2207-2217, doi: 10.6038/cjg20130708.

Zhou L T, Huang R H. 2008. Interdecadal variability of sensible heat in arid and semi-arid regions of Northwest China and its relation to summer precipitation in China.ChineseJournalofAtmosphericSciences(in Chinese), 32(6): 1276-1288.

附中文參考文獻(xiàn)

黃榮輝, 周德剛, 陳文等. 2013. 關(guān)于中國(guó)西北干旱區(qū)陸—?dú)庀嗷プ饔眉捌鋵?duì)氣候影響研究的最近進(jìn)展. 大氣科學(xué), 37(2): 189-210, doi: 10.3878/j.issn.1006-9895.2012.12303.

劉和平, 劉樹(shù)華, 桑建國(guó). 1999. 不同下墊面水分與能量傳輸模式. 大氣科學(xué), 23(4): 449-460.

劉樹(shù)華, 劉振鑫, 鄭輝等. 2013. 多尺度大氣邊界層與陸面物理過(guò)程模式的研究進(jìn)展. 中國(guó)科學(xué): 物理學(xué) 力學(xué) 天文學(xué), 43(10): 1332-1355.

鄭輝, 劉樹(shù)華. 2012. 數(shù)值差分格式及格點(diǎn)設(shè)置對(duì)土壤溫度模擬結(jié)果的影響. 地球物理學(xué)報(bào), 55(8): 2514-2522, doi: 10.6038/j.issn.0001-5733.2012.08.004.

鄭輝, 劉樹(shù)華. 2013. 沙漠陸面過(guò)程參數(shù)化與模擬. 地球物理學(xué)報(bào), 56(7): 2207-2217, doi: 10.6038/cjg20130708.

周連童, 黃榮輝. 2008. 中國(guó)西北干旱、半干旱區(qū)感熱的年代際變化特征及其與中國(guó)夏季降水的關(guān)系. 大氣科學(xué), 32(6): 1276-1288.

(本文編輯張正峰)

基金項(xiàng)目國(guó)家科技支撐計(jì)劃(2012BAH29B03)資助.

作者簡(jiǎn)介鄭輝，男，1988年生，博士研究生，主要從事大氣邊界層物理和區(qū)域氣候變化研究. E-mail: zhenghui@pku.edu.cn *通訊作者劉樹(shù)華，教授，博士生導(dǎo)師，主要從事大氣邊界層物理和區(qū)域氣候變化研究. E-mail:lshuhua@pku.edu.cn

doi:10.6038/cjg20160107 中圖分類(lèi)號(hào)P404

收稿日期2014-02-25，2015-05-16收修定稿

Description and evaluation of the Peking University Land Model (PKULM)

ZHENG Hui1，LIU Shu-Hua1*， Prabhakar Clement2， LIU Zhen-Xin3， HOU Xu-Hong4，WANG Shu1, ZHAO Jing-Chuan1，LI Yuan1，MIAO Yu-Cong1，ZHENG Yi-Jia1， SHENG Li5， ZHU Lin6

1DepartmentofAtmosphericandOceanicScience,SchoolofPhysics,PekingUniversity,Beijing100871,China2DepartmentofCivilEngineering,AuburnUniversity,Auburn,AL36849,USA3InstituteofAtmosphericPhysics,ChineseAcademyofSciences,Beijing100029,China4ColdandAridRegionsEnvironmentalandEngineeringResearch,ChineseAcademyofSciences,Lanzhou730000,China5NationalMeteorologicalCenter,ChinaMeteorologicalAdministration,Beijing100081,China6NationalCenterforSpaceWeather,ChinaMeteorologicalAdministration,Beijing100081,China

AbstractThe land surface model is a key component of climate and weather models. It is intractable to accurate modeling land-atmosphere interaction in arid and semi-arid regions. Evolved from the Soil-Plant-Atmosphere Model (SPAM), a new Peking University Land Model (PKULM) has been developed with several augments in physical processes and numerical solution for application in arid and semi-arid environments.

KeywordsLand surface model； Surface energy balance； Soil moisture； Stomatal conductance

鄭輝， 劉樹(shù)華， Prabhakar C等. 2016. 北京大學(xué)陸面過(guò)程模式PKULM(Peking University Land Model)介紹及檢驗(yàn).地球物理學(xué)報(bào)，59(1):79-92,doi:10.6038/cjg20160107.

Zheng H, Liu S H, Prabhakar C, et al. 2016. Description and evaluation of the Peking University Land Model (PKULM).ChineseJ.Geophys. (in Chinese),59(1):79-92,doi:10.6038/cjg20160107.