北京大学陆面过程模式PKULM(Peking University Land Model)介绍及检验

引用本文

郑辉, 刘树华, Prabhakar C等. 2016. 北京大学陆面过程模式PKULM(Peking University Land Model)介绍及检验. 地球物理学报, 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). Chinese J. Geophysics, 59(1): 79-92, doi: 10.6038/cjg20160107 复制到剪切板

北京大学陆面过程模式PKULM(Peking University Land Model)介绍及检验

郑辉¹, 刘树华¹, Prabhakar Clement², 刘振鑫³, 候旭宏⁴, 王姝¹, 赵靖川¹, 李源¹, 缪育聪¹, 郑亦佳¹, 盛黎⁵, 朱琳⁶

1. 北京大学物理学院大气与海洋科学系, 北京 100871;
2. Department of Civil Engineering, Auburn University, Auburn, AL 36849;
3. 中国科学院大气物理研究所, 北京 100029;
4. 中国科学院寒区旱区环境与工程研究所, 兰州 730000;
5. 中国气象局国家气象中心, 北京 100081;
6. 中国气象局卫星气象中心, 北京 100081

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

基金项目: 国家科技支撑计划(2012BAH29B03)资助.

作者简介: 郑辉,男,1988年生,博士研究生,主要从事大气边界层物理和区域气候变化研究.E-mail:zhenghui@pku.edu.cn

通讯作者: 刘树华,教授,博士生导师,主要从事大气边界层物理和区域气候变化研究.E-mail:lshuhua@pku.edu.cn

摘要: 陆面过程模式是气候模式和天气模式的核心组成部分之一.在土壤-植被-大气耦合模式(Soil-Plant-Atmosphere Model, SPAM)的基础上,发展了新一代北京大学陆面过程模式PKULM(Peking University Land Model).本文首先介绍了PKULM的辐射传输、湍流输送、光合作用、土壤水热输送等过程的参数化方案;采用隐式迭代计算框架,发展并应用了一个快速的线性方程组求解算法,提高了模式计算稳定性;提出并使用了二分搜索算法计算气孔阻抗,避免了CLM(Community Land Model)等使用的迭代方法在干旱区不稳定的情况,提高了模式的适用性;采用水势为基础的土壤水分扩散方程,使模式能够模拟土壤饱和区的水分输送过程,为进一步与水文过程模式耦合奠定了基础;还发展了一个地表积水与径流过程的机理模型,提高了模式对地表水分平衡过程的模拟能力;最后,使用"中国西北干旱区陆-气相互作用观测试验"平凉站的资料对模式进行了检验并与NOAH(National Center for Environmental Prediction, Oregon State University, Air Force, and Hydrology Lab model)陆面过程模式的模拟结果进行了比较,结果表明PKULM能够较好地模拟西北半干旱区农田下垫面地气交换过程.

关键词: 陆面过程模式地表能量平衡土壤湿度气孔导度

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

ZHENG Hui¹, LIU Shu-Hua¹, Prabhakar Clement², LIU Zhen-Xin³, HOU Xu-Hong⁴, WANG Shu¹, ZHAO Jing-Chuan¹, LI Yuan¹, MIAO Yu-Cong¹, ZHENG Yi-Jia¹, SHENG Li⁵, ZHU Lin⁶

1. Department of Atmospheric and Oceanic Science, School of Physics, Peking University, Beijing 100871, China;
2. Department of Civil Engineering, Auburn University, Auburn, AL 36849, USA;
3. Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China;
4. Cold and Arid Regions Environmental and Engineering Research, Chinese Academy of Sciences, Lanzhou 730000, China;
5. National Meteorological Center, China Meteorological Administration, Beijing 100081, China;
6. National Center for Space Weather, China Meteorological Administration, Beijing 100081, China

Abstract: The 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.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^-2 for 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.

Key words: Land surface model Surface energy balance Soil moisture Stomatal conductance

1 引言

陆地表面约占地球表面积的29.2%，陆面过程是影响地球系统能量水分循环的重要物理过程.一方面，作为大气圈的下垫面，陆地表面控制着太阳净辐射在长波辐射、感热和潜热之间的分配比例，通过湍流输送过程加热大气、驱动大气运动，影响区域和全球气候(周连童和黄荣辉，2008; Yang et al.，2011)；另一方面，作为气候系统中能量和水分的重要储库，土壤温度和湿度与大气环流有长周期的相互作用(Koster et al.，2004)；同时，作为生物活动的主要场所，陆地表面的生物物理、化学特性还对陆气交换过程乃至未来气候变化评估产生着重要的影响(Feddema et al.，2005; Davidson et al.，2006; Jiang et al.，2011).

发展陆面过程模式，准确模拟陆气间能量物质交换过程，是提高气候系统可预报性的重要途径(Sellers et al.，1997; Pitman，2003).Manabe(1969)发展的“水桶”模型，以地表能量平衡和水分平衡方程作为核心，是建立陆面过程模式的第一次尝试.Deardorff(1978)在他的模式中使用了近地面湍流相似性理论计算陆气间感热和潜热通量，并第一次考虑了植被冠层作用，其工作在陆面过程模式的发展过程中具有划时代的意义.这种考虑了植被冠层作用的陆面过程模式也被称为第二代陆面过程模式；BATS(Dickinson，1984)、SiB(Sellers et al.，1986)、IAP94(Dai and Zeng，1997)是其中的典型代表.20世纪末，考虑植被光合作用的SiB2(Sellers et al.，1996)、CLM(Bonan et al.，2002)、CoLM(Dai et al.，2003)等一大批第三代陆面过程模式相继出现.第三代陆面过程模式以详细考虑光合作用为主要特点，综合考虑了植被冠层内辐射传输，土壤水热输送、植被生化作用、近地面湍流输送等过程，具有较高的模拟精度，是当前陆面过程模式发展的前沿方向之一.

Liu(1996)在Noilhan和Planton(1989)的N89模式的基础上，发展了一个包含植被遮蔽和蒸腾作用的陆面水热交换模式(LSEM，L and Surface Exchange Model).该模式很好地模拟森林和沙漠绿洲下垫面温度变化及能量交换特征.刘和平等(1999)引入Darcy定律描述土壤水分输送过程，并将模式扩展到草原、戈壁及热带森林等多种下垫面，经过实测资料验证，证明了模式能够模拟多种复杂下垫面地气间能量、水分交换过程.Liu等(2004a，b)将陆面过程模式与北京大学大气边界层模式相互耦合，发展了土壤—植被—大气耦合模式(SPAM，Soil-Plant-Atmosphere Model)，准确模拟了沙漠—绿洲系统的“冷岛效应”和“湿岛效应”.Liu等(2006)改进了辐射传输过程，采用了更加适合干旱区半干旱区的通量廓线关系，模拟的地表温度有了大幅度的改进.郑辉和刘树华(2012，2013)和刘树华等(2013)从土壤热传导、湍流输送、植物生理生化过程等多个方面进一步改进模式，经过实测资料验证，取得了较好的模拟结果.

本文在已有工作的基础上，考虑了植被以及地表水分平衡过程，建立了北京大学陆面过程模式PKULM.本文首先介绍PKULM的模式结构、参数化方案以及计算方法，然后利用西北半干旱区农田下垫面观测资料驱动模式，对模式进行了检验，评估模式的模拟能力.

2 模式介绍

如图 1所示，北京大学陆面过程模式包含如下五个相互耦合的过程：(1)辐射传输过程，确定不同波段辐射能量在植被冠层和土壤表面之间的分配，为光合作用、感热输送和潜热输送提供能量来源；(2)湍流输送过程，确定热量、水分和二氧化碳等的湍流输送通量及边界层大气运动的摩擦阻力；(3)土壤热力学过程，确定土壤中水分、热量输送及其与大气相互作用的物理过程；(4)植物生理过程，确定植物蒸腾速率及光合作用强度的生物理化过程；(5)水分平衡过程，确定降水量的分配，估算土壤含水量、径流量、蒸发量、蒸腾量.目前，模式不包含积雪与冻土过程.

图 1 北京大学陆面过程模式示意图 Fig. 1 Diagram of Peking University Land Surface Model

在PKULM中，每一个格点上可以有m种不同的植被冠层覆盖，分别计算每一种冠层所覆盖地表的短波辐射通量、长波辐射通量、感热通量、水汽通量以及二氧化碳通量，然后以相应的植被覆盖率σ_j为权重加总，即可得到陆面与大气之间总的能量、物质交换通量.

2.1 辐射传输过程 2.1.1 太阳辐射

由于冠层的光学特性在不同波段存在显著差异，PKULM将太阳辐射分为0~700 nm和700~2800 nm两个波段Λ分别计算；在每一个波段，还区分对待冠层顶所接收到的直射和散射太阳光.

在植被覆盖的地区，太阳辐射不断在植物叶片间反射、透射和吸收，最终光强趋于稳态分布，该分布可以用二流模型近似为(Dickinson，1984)：

其中，I↑和I↓分别是冠层中向上和向下的散射光(入射辐射归一化值)，μ为太阳天顶角的余弦，μ±s是单位叶面积散射光学厚度的倒数，K是单位叶面积直射光的光学厚度；β、β₀分别是散射光和直射光的向上散射系数，ω为叶片散射系数，L和S分别为叶面积指数和茎面积指数(m²·m^-2).对于不同的植被类型，光学性质参数β、β₀以及ω有不同的取值.

根据二流模型，对于Λ波段的太阳辐射，冠层对直射光和散射光的吸收率I_Λ^μ和I_Λ分别为：

其中，α_g，Λ^μ和α_g，Λ分别是在Λ波段地表对直射光和散射光的反射率.

冠层和地表吸收的净辐射通量S_v，net和S_g，net(W·m^-2)分别为：

其中，S_atm↓_Λ^μ为Λ波段入射的太阳直射辐射(W·m^-2)，S_atm↓_Λ为Λ波段入射散射辐射(W·m^-2). 2.1.2 长波辐射

图 2是PKULM长波辐射传输示意图.L_atm↓是大气向下的长波辐射强迫(W·m^-2).根据斯特芬-波尔兹曼定律以及比辐射率的定义，植被冠层向下的长波辐射通量L_v↓(W·m^-2)等于L_atm↓透过冠层的部分与冠层热发射值之和，即：

其中，ε_v是冠层的比辐射率，T_v是冠层温度(K).同理，地表向上和冠层向上的长波辐射通量L_g↑和L↑(W·m^-2)分别为：

其中，ε_g是地表的比辐射率，T_g是地表温度(K).因此，冠层和地表吸收的净长波辐射L_v，net、L_g，net分别为：

图 2 PKULM长波辐射传输示意图 Fig. 2 Diagram of PKULM longwave radiation transfer scheme

2.2 湍流通量 2.2.1 植被冠层感热、潜热通量

PKULM的湍流输送过程采用阻抗模型计算.土壤表面、冠层以及大气之间的感热、潜热输送通道如图 3所示，其湍流阻抗采用Monin-Obukhov相似性理论计算(Zeng et al.，1998).

图 3 地气间感热、潜热输送通道及其阻抗示意图 Fig. 3 Schematic diagram of transfer pathway for sensible heat and moisture/latent heat between land surface and reference height

假设冠层间空气的热量存储可以忽略，陆地表面输送至大气的感热通量H(W·m^-2)等于冠层输送至冠层间空气的感热通量H_v(W·m^-2)与地表输送至冠层间空气的感热通量H_g之和为(W·m^-2)：

其中，ρ_atm为空气密度(kg·m^-3)，C_p为空气的定压热容(J·kg^-1·K^-1)，θ_atm为大气位温(K)，T_s为冠层间空气温度(K)，r_ah、r′_ah分别为大气到冠层间空气、冠层间空气到地表的热力阻抗(s·m^-1)，r_b为叶片边界层阻抗(s·m^-1).

联立式(12)至(15)，可以求出冠层间空气温度T_s的表达式为

其中，系数分别为c_ah=1/r_ah，c_gh=1/r′_ah，c_vh=(L+S)/r_b.

同样陆地表面输送至大气的水汽通量E(kg·m^-2·s^-1)等于冠层输送至冠层间空气的水汽通量E_v(kg·m^-2·s^-1)与地表输送至冠层间空气的水汽通量E_g(kg·m^-2·s^-1)之和为：

其中，q_atm、q_s、q_g分别为大气，冠层间空气以及表层土壤空隙间空气的比湿(kg·kg^-1)，q_sat(T_v)表示温度T_v下的饱和比湿(kg·kg^-1)，β_soil是表层土壤含水量相关的经验函数(Sakaguchi and Zeng，2009)，r′_aw、r_aw分别为地表到冠层间空气以及冠层间空气到大气间的水汽阻抗(s·m^-1)，r_total表示冠层到冠层间空气的总阻抗(s·m^-1)，是气孔阻抗r_s以及叶片边界层阻抗r_b的函数.

联立式(17)至(20)，可以求出冠层间空气比湿q_s的表达式为

其中，系数分别为c_aw=1/r_ah，c_gw=β_soil/r′_ah，c_vw=1/r_total. 2.2.2 光合作用及气孔导度

由于植物叶片覆盖了物质渗透率极低的角质层，气孔成为了植被与大气之间进行物质交换的主要通道，气孔的开闭大小直接影响了植被冠层蒸腾作用的强度.第三代陆面过程模型均采用Ball-Berry半经验模型(Collatz et al.，1991)计算气孔阻抗r_s(s·m^-1)：

其中m、b为经验系数，P_atm为大气压(Pa)，c_s为叶片表面二氧化碳分压(Pa)，e_s为叶片表面水汽分压(Pa)，e_i为气孔内水汽分压(Pa).二氧化碳同化速率A(μmol·m^-2 ·s^-1)采用Collatz等(1991，1992)方案计算，它是冠层温度、土壤水势、光合作用有效辐射以及二氧化碳、水气分压的函数.

根据阻抗模型，二氧化碳同化速率A，气孔内二氧化碳和水汽分压c_i、e_i，以及大气中二氧化碳和水汽分压c_a、e_a有如下关系：

2.2.3 二分查找算法

为了求解气孔阻抗，令函数h(c_i)=c_a-(1.37r_b+1.65r_s)P_atmA.由于气孔阻抗r_s和光合作用强度A都是c_i的函数，因此h(c_i)也是c_i的函数.利用h(c_i)，式(23)等价于如下的定点问题：

求出使上述等式成立的二氧化碳分压c_i^*，将其带入式(22)，即可计算气孔阻抗r_s.

CLM等陆面过程模式采用迭代算法计算c_i^*，即首先给c_i一个初始值，然后带入函数h(c_i)，将得到的结果更新c_i，如此循环直到收敛.但Sun等(2012)指出该方法在大气水汽压较低的情况下并不收敛，导致全球模式中初级生产力计算结果出现大面积异常.如图 4所示，当水汽压e_a=317Pa(相对湿度10%)时，若c_i初始值取为28 Pa，则下一循环中c_i=h(c_i)为0 Pa，将c_i=0 Pa带入函数h(c_i)，得到c_i=28 Pa，因此c_i的值将在0 Pa与28 Pa之间震荡而不能收敛.

图 4 不同水汽压下PKULM的气孔导度定点函数
参数设置：c_a = 28 Pa, t=25 ℃, r_b=1×10^-6 s·m² μmol^-1, v_cmax = 40 μmol CO₂ m^-2·s^-1.Fig. 4 Fixed-point function for photosynthesis-stomatal conductance model used in PKULM at multiple vapor pressure
Model parameters used to generate these curves are: c_a=28 Pa, t=25 ℃, r_b=1×10^-6 s·m² μmol^-1, v_cmax=40 μmol CO₂ m^-2·s^-1.

PKULM采用二分查找算法，计算定点问题式(24)的解.从图 4可以看到，当c_i=0时，A=0，h(c_i)=c_a>c_i=0；而当c_i→+∞时，A→0，h(c_i)→0<c_i→+∞；所以定点问题式(24)一定有解.并且，在c_i从0增加至正无穷的过程中，函数h(c_i)的导数值一直小于1，因此c_i-h(c_i)单调变化，式(24)的解为唯一解.二分查找算法的具体求解过程如图 5所示：(1)首先给出叶片内二氧化碳分压c_i的初始查找区间[c_i^l，c_i_h]；(2)计算区间中点c_i^m=(c_i^l+c_i^h)/2处定点函数的值h(c_i^m)，若c_i^m<h(c_i^m)，则继续查找左半区间[c_i^l，c_i^m]，否则继续查找右半区间[c_i^m，c_i^h]；(3)重复步骤(2)，直到区间范围c_i^h-c_i^l小于0.01 Pa或者c_i^mi=h(c_i^m)，区间中点c_i^m即可作为方程式(24)的解.

图 5 计算气孔内二氧化碳分压的二分查找法
参数设置：c_a=28 Pa, e_a=634 Pa, t=25 ℃, r_b=1×10^-6 s·m² μmol^-1, v_cmax=40 μmol CO₂ m^-2·s^-1. Fig. 5 Schematic description of the bisection method for calculating CO₂ partial pressure inside leaf
Model parameters used to generate this curve are: c_a=28 Pa, e_a=634 Pa, t=25 ℃, r_b=1×10^-6 s·m² μmol^-1, v_cmax=40 μmol CO₂ m^-2·s^-1.

2.3 冠层及地表能量平衡

PKULM每一个格点上可以有m种不同的植被冠层覆盖，对于第j种植被冠层，根据能量平衡，冠层吸收的太阳净辐射和净长波辐射全部通过湍流输送过程加热大气和地表，即：

对于地表，其吸收的太阳净辐射和净长波辐射等于地表向上的感热、潜热通量与向下的土壤热通量之和，即：

可以将上述两式写为

2.4 土壤热量输送

土壤中热量的存储和释放过程可以用热力学第一定律和梯度热传导方程描述为：

其中，t为时间(s)，z为土壤深度(m)，T为土壤温度(K)，F为土壤热通量(W·m^-2，向上为正)，C为土壤体积热容(J·m^-3·K^-1)，k_T为土壤热传导系数(W·m^-1·K^-1).

把土壤分为n层，采用隐式差分方案，将上述两式离散化，则第i层土壤温度的预报方程为：

其中，z_h，i表示第i层土壤与第i+1层土壤界面的深度(m)，F_s，i是该界面处的土壤热通量(W·m^-2)，k_T[z_h，i]该界面处土壤热传导系数(W·m^-1·K^-1)，Δt为时间步长(s)，Δz_i为第i层土壤的厚度(m)，pT_i表示前一时间步第i层土壤的温度(K，记号p表示状态量在上一时间步的取值，下同).

可以将式(31)写成三对角矩阵形式，以便进一步数值求解，公式为

2.5 冠层降水截留

大气降水q_rain(kg·m^-2·s^-1)在经过植被冠层时，有一部分直接降落至地面q_thru(kg·m^-2·s^-1)，而另一部分被截留q_intr(kg·m^-2·s^-1)，然后部分滴落至地表q_drip(kg·m^-2·s^-1)，没有滴落的部分依附在叶茎表面，成为冠层蓄水w_can(kg·m^-2).冠层截流降水量q_intr、直接落至地表的降水量q_thru以及冠层滴落水量q_drip分别为(Lawrence et al.，2007)：

其中pw_can是前一时间步上冠层蓄水量(kg·m^-2)，w_can，max=0.1(L+S)是最大冠层蓄水量(kg·m^-2)，L和S分别为叶面积指数和茎面积指数(m²·m^-2).

最终落在地表的总水量q_grnd(kg·m^-2·s^-1)为

2.6 地表径流及渗流

PKULM发展了一个地表径流及渗流的机理模型.在降水过程中，地表从水分不饱和状态到产生径流将经历如下三个阶段：

(1)土壤非饱和阶段.该阶段，表层土壤毛细水势小于其饱和水势ψ₁<ψ_sat，1，根据地表水分平衡，降落在地表的总水量q_grnd等于地表蒸发量E_g与渗流量q_infl(kg·m^-2)之和，公式为

(2)土壤饱和、地表出现积水阶段.该阶段，土壤毛细水势大于其饱和水势，且小于地表最大积水深度w_g，max=0.2×10^-3m(Sellers et al.，1996)；降落在地表的总水量与地表蒸发量、地表渗流量之差等于积水的增加量，公式为

(3)径流产生阶段.该阶段，土壤毛细水势等于地表最大积水深度，土壤饱和，积水深度保持w_g，max不变，降落在地表的总水量与地表蒸发量、地表渗流量之差等于径流，公式为

上述三式中，(E_g)_j是第j种植被冠层所覆盖地表的蒸发量(kg·m^-2·s^-1)，ρ_liq是水的密度(kg·m^-3)，q_infl是土壤表面渗流(kg·m^-2·s^-1)，其大小与第一层土壤导水率k_h，1(kg·m^-2·s^-1)，土壤毛细水势ψ₁(m)以及深度z₁(m)有关，公式为

值得注意的是，我们采用统一的变量w_g表示表层土壤毛细水势以及地表积水深度：当w_g≤0时，w_g表示表层土壤毛细水势(m)，而当w_g>0时，w_g表示地表积水深度(m).使用这种表示方法，再利用适当的辅助函数h_g，我们可以将不同阶段的三个地表水分平衡方程式(39)、式(40)和式(41)写成统一的差分形式为：

其中pw_g表示前一时间步上w_g的值.

上式可进一步写为

其中：

2.7 土壤湿度

土壤水分运动可以采用Richards方程描述为(Zeng and Decker，2009)：

其中，w为土壤体积含水量(m·m^-3)，q为土壤水分通量(kg·m^-2·s^-1)，k_h为土壤导水率(kg·m^-2·s^-1)，ψ为土壤毛细水势(m)，q_root为植物根系吸水量(kg·m^-1·s^-1).

为了数值求解上述两个方程，需要利用土壤含水量与土壤毛细水势的函数关系，将方程转化成只含毛细水势或只含土壤含水量的偏微分方程(Ross，2003)：

其中，w_sat是土壤饱和含水量(m·m^-3)，ψ_sat为土壤饱和水势(m)，B为Clapp-Hornberger参数.

CLM、NOAH等陆面过程模式一般将土壤含水量作为预报变量，即模拟土壤含水量变化，然后利用式(51)计算土壤毛细水势.但Ross(2003)指出，由于土壤含水量与土壤毛细水势的一一映射关系只在非饱和时成立，因此这种以土壤含水量为基础的水分运动模式仅能模拟土壤非饱和区(地下水位以上)的水分运动.如图 6所示，在土壤饱和区，土壤含水量维持在饱和含水量不变，根据式(53)，土壤毛细水势保持不变；但实际中，毛细水势随深度增加而增加.

图 6 土壤毛细水势和土壤体积含水量廓线
参数设置：饱和含水量为0.3731 m³·m^-3，饱和土壤毛细水势为-0.0473 m，Clapp-Hornbeger参数为 3.387，地下水位为 -1.737 m. Fig. 6 Schematic diagram of soil capillary pressure head and volumetric water content profile
Model parameters used to generate these curves: Soil volumetric water content at saturation is 0.3731 m³·m^-3, saturated soil capillary pressure head is -0.0473 m, Clapp-Hornberger parameter is 3.387, water table depth is -1.737 m.

PKULM采用土壤水势为基础的Richards方程.使用隐式差分方案，将式(49)和式(50)离散化，得到第i层土壤毛细水势的预报方程为

其中，pψ_i为前一时间步第i层土壤的毛细水势(m)；k_h[z_h，i]为第i层与第i+1层土壤界面处的导水率，公式为

式(52)和式(53)可以进一步写成三对角矩阵形式为

其中：

3 控制方程及隐式Picard迭代计算算法

PKULM的核心是能量与水分平衡控制方程组.冠层能量平衡方程式(27)、地表能量平衡方程式(28)以及土壤热扩散方程式(34)构成了PKULM的能量平衡控制方程组.方程组可以写成如下的矩阵形式为

地表水分平衡方程式(53)以及Richards方程式(54)构成了PKULM的水分平衡控制方程组.方程组可以写成如下的矩阵形式为

求解上述两个非线性矩阵方程，即可获得PKULM全部的状态参量：土壤温度T_i和水势ψ_i、地表温度T_g和积水深度w′_g以及冠层温度T_v，j，在此基础上，可以计算地气间全部的能量与物质交换通量.

为了方便讨论，将式(59)和式(60)记作Aφ=B，其中φ是全部状态参量所构成的向量，A是全部系数构成的矩阵，B是全部常数项构成的矩阵.由于矩阵A、B都是关于φ的非线性函数，我们采用Picard迭代法(Hipsey et al.，2004)将方程线性化并求解，具体做法为：(1)使用上一时刻状态量的值pφ作为φ的初值；(2)将上一迭代步状态量的值φ^l，带入系数矩阵，求解线性方程组A(φ^l)φ^l+1=B(φ^l)，以方程的解φ^l+1作为当前迭代步状态量φ的值；(3)经过若干步迭代，当φ^l+1≈φ^l时，可认为φ^l+1就是非线性方程组A(φ)φ=B(φ)的解.

对于形如式(60)的线性方程组，可以直接采用Thomas算法快速求解.而对于形如式(59)的线性方程组，需要对Thomas算法进行修改，具体算法如下：(1)依次用矩阵第j(1≤j≤m)行消去第m+1行的各项系数h_v，j(1≤j≤m)；(2)采用标准Thomas算法求解方程组第m+1行至第m+n+1行；(3)最后，用矩阵第m+1行依次消去第m+1列各项系数v_v，j(1≤j≤m).

4 观测资料介绍

本文使用“我国西北典型干旱半干旱区能量和水分循环观测试验与分析”研究项目(黄荣辉等，2013)平凉观测站2010年加强观测期资料作为模式输入，进行模拟试验.平凉站位于甘肃省平凉市(35°35′N，106°42′E)，海拔为1417 m.地表平坦均一，土壤为黄土，下垫面作物为玉米.平凉站的主要观测仪器包括：四分量辐射计(CNR1，Kipp & Zonen)、三维超声风速仪(CSAT3，Campell)、温度和相对湿度传感器(HMP45C，Vaisala)、大气压传感器(CS105，Campbell)、CO₂/H₂O开路气体分析仪(LI-7500，LI-COR)、4层(5，10，20，40 cm)土壤温度测量系统(CS616-L，Campbell)和4层(5 cm，10，20，40 cm)土壤湿度测量系统(109-L，Campbell).2010年8月14日19时至2010年8月16日09时，无天气系统过境，观测记录完整，且根据稳态测试(SST，Steady State Test)和整体湍流特征检验(ITC，Integral Turbulence Characteristics)对湍流资料进行了质量控制(Foken，2008).

5 模拟结果

模式以2010年8月14日19时观测得到的各层土壤温度和湿度资料作为初始状态，以平均风速、平均气温、相对湿度、大气压、二氧化碳分压以及向下的长短波辐射作为输入量，以半小时为步长向前积分，从2010年8月15日06时起输出模拟结果.同时使用相同的大气强迫资料和初始土壤温度、湿度场驱动NOAH陆面过程模式，将模拟结果与PKULM进行比较.

分别使用平均偏差(Mean Bias Error，MBE)、均方根误差(Root Mean Square Error，RMSE)和相关系数(Correlation Coefficient，CC)分析模拟误差，定义如下：

其中，X_obs表示观测值，X_sim表示模拟值，E(X)=[∑_i=1^NX(i)]/N为序列X的期望.表 1列出了PKULM和NOAH模拟的向上短波辐射、向上长波辐射、感热通量和潜热通量的误差.

表 1 PKULM和NOAH模拟能量通量的平均误差(MBE)、标准差(RMSE)和相关系数(CC) Table 1 Mean bias error (MBE), root mean square error (RMSE) and correlation coefficient (CC) of energy flux of PKULM and NOAH

图 7是PKULM和NOAH模拟的反射短波辐射及其与观测值的比较.中午12时，平凉站观测得到的反射短波辐射达到最大值，约为140 W·m^-2，此时PKULM模拟的结果略高于观测值0.39 W·m^-2，而NOAH模拟结果则低于观测值约4.0 W·m^-2.总体而言，两个模式均能较好地模拟短波辐射传输过程，模拟的反射短波辐射均方根误差分别为4.061 W·m^-2和11.34 W·m^-2.

图 7 向上短波辐射的观测值与PKULM和NOAH模拟值的比较图 Fig. 7 Comparison of observed, PKULM modeled and NOAH modeled upward solar radiation

图 8是PKULM和NOAH模拟的地表发射长波辐射及其与观测值的比较.可以看到，PKULM的模拟结果与观测值接近，模拟均方根误差为6.608 W·m^-2.NOAH的模拟结果偏低，这可能与其感热通量模拟偏高(图 9)，导致地表温度模拟结果偏低(图 11)造成的.

图 8 向上长波辐射的观测值与PKULM和NOAH模拟值的比较图 Fig. 8 Comparison of observed, PKULM modeled and NOAH modeled upward longwave radiation

图 9 感热通量的观测值与PKULM和NOAH模拟值的比较图 Fig. 9 Comparison of observed, PKULM modeled and NOAH modeled sensible heat flux

图 10 潜热通量的观测值与PKULM和NOAH模拟值的比较图 Fig. 10 Comparison of observed, PKULM modeled, and NOAH modeled latent heat flux

图 9和图 10分别给出了感热通量和潜热通量的模拟值与观测值.从图中可以看到，由于植被蒸腾作用，平凉农田下垫面上潜热通量占主导地位，潜热通量大小约为感热通量的两倍.两个模式均能比较好地模拟潜热通量大小及日变化过程，PKULM和NOAH模拟值的均方根误差分别为21.81 W·m^-2和24.05 W·m^-2.而对感热通量而言，PKULM的模拟结果更接近观测值，均方根误差为16.29 W·m^-2，而NOAH模拟结果的均方根误差为29.95 W·m^-2.

图 11和图 12分别是5 cm深度土壤温度、湿度的模拟值及其与观测值的比较.PKULM的模拟值均比NOAH的模拟值更接近实测值.对于土壤温度，PKULM模拟的温度日变化振幅与相位均与实测值相符，平均偏差为0.289 K，均方根误差为0.732 K，相关系数为0.978，而NOAH模拟的土壤温度偏低，平均偏差为-1.03 K，均方根误差为1.79 K，相关系数为0.903.对于土壤湿度，PKULM模拟结果的平均偏差为-7.4×10^-3 m³·m^-3，均方根误差为7.5×10^-3 m³·m^-3，相关系数为0.906；NOAH模拟结果的平均偏差为-2.7×10^-2 m³·m^-3，均方根误差为3.0×10^-2 m³·m^-3，相关系数为0.792.两个模式模拟的土壤湿度变化趋势与实际一致，但大小与实测有一定偏差.两个模式模拟的白天土壤湿度降低速度高于实测，而夜晚低于实测，模式对土壤水分运动的参数化方案仍然有待改进.

图 11 5 cm深度土壤温度观测值与PKULM和NOAH模拟值的比较图 Fig. 11 Comparison of observed, PKULM modeled, and NOAH modeled soil temperature at 5 cm depth

图 12 5 cm深度土壤湿度观测值与PKULM和NOAH模拟值比较图 Fig. 12 Comparison of observed, PKULM modeled, and NOAH modeled soil moisture at 5 cm depth

6 结论

本文介绍了北京大学陆面过程模式PKULM的参数化方案和计算算法，并利用平凉农田下垫面观测资料对模式进行了检验.主要结论如下：

(1)通过使用二分查找算法计算气孔导度，避免了CLM等模式使用的迭代算法在低水汽压环境下不稳定的问题，提高了模式在干旱和半干旱地区的适用性.

(2)本文发展了一套求解陆面过程非线性方程组的Picard隐式迭代算法，该算法计算稳定，且有效地将各物理过程参数化方案与模式总体框架隔离开来，有利于模式的进一步发展.

(3)PKULM能够模拟半干旱区农田下垫面陆气交换过程.模式模拟的短波辐射、长波辐射、感热和潜热通量与实测值吻合，模拟误差略低于NOAH陆面过程模式.

(4)模式对土壤湿度的模拟误差较大，与土壤湿度密切相关的根系过程、土壤水分输送过程仍然需要进一步改进.

参考文献

[1]	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.
[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. Forest Meteor., 54(2-4):107-136, doi:10.1016/0168-1923(91)90002-8.
[3]	Collatz G J, Ribas-Carbo M, Berry J A. 1992. Coupled photosynthesis-stomatal conductance model for leaves of C4 plants. Aust. J. Plant Physiol., 19(5):519-538.
[4]	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.
[5]	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.
[6]	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.
[7]	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.
[8]	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.
[9]	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.
[10]	Foken T. 2008. Micrometeorology. Berlin:Springer, 1-308.
[11]	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. Fluid Mech., 4(1):79-106, doi:10.1023/A:1025547707198.
[12]	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. Chinese Journal of Atmospheric Sciences(in Chinese), 37(2):189-210, doi:10.3878/j.issn.1006-9895.2012.12303.
[13]	Jiang D B, Zhang Y, Lang X M. 2011. Vegetation feedback under future global warming. Theoretical and Applied Climatology, 106(1-2):211-227.
[14]	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.
[15]	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.
[16]	Liu H P, Liu S H, Sang J G. 1999. A modified SiB to simulate momentum, heat and water transfer over various underlying surfaces. Chinese Journal of Atmospheric Sciences(in Chinese), 23(4):449-460.
[17]	Liu S H. 1996. A parameterized model for the influence of vegetation on the moisture-heat exchange at near-ground layer. Acta Sci. Nat. Univ. Pek., 32(1):69-77.
[18]	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.
[19]	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.
[20]	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.
[21]	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.
[22]	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.
[23]	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.
[24]	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.
[25]	Ross P J. 2003. Modeling soil water and solute transport-fast, simplified numerical solutions. Agron. J., 95(6):1352-1361.
[26]	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.
[27]	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.
[28]	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.
[29]	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.
[30]	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.
[31]	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.
[32]	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.
[33]	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.
[34]	Zheng H, Liu S H. 2012. Effects of difference and grid scheme in soil temperature simulation. Chinese J. Geophys.(in Chinese), 55(8):2514-2522, doi:10.6038/j.issn.0001-5733.2012.08.004.
[35]	Zheng H, Liu S H. 2013. Land surface parameterization and modeling over desert. Chinese J. Geophys.(in Chinese), 56(7):2207-2217, doi:10.6038/cjg20130708.
[36]	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. Chinese Journal of Atmospheric Sciences(in Chinese), 32(6):1276-1288.
[37]	黄荣辉, 周德刚, 陈文等. 2013. 关于中国西北干旱区陆-气相互作用及其对气候影响研究的最近进展. 大气科学, 37(2):189-210, doi:10.3878/j.issn.1006-9895.2012.12303.
[38]	刘和平, 刘树华, 桑建国. 1999. 不同下垫面水分与能量传输模式. 大气科学, 23(4):449-460.
[39]	刘树华, 刘振鑫, 郑辉等. 2013. 多尺度大气边界层与陆面物理过程模式的研究进展. 中国科学:物理学力学天文学, 43(10):1332-1355.
[40]	郑辉, 刘树华. 2012. 数值差分格式及格点设置对土壤温度模拟结果的影响. 地球物理学报, 55(8):2514-2522, doi:10.6038/j.issn.0001-5733.2012.08.004.
[41]	郑辉, 刘树华. 2013. 沙漠陆面过程参数化与模拟. 地球物理学报, 56(7):2207-2217, doi:10.6038/cjg20130708.
[42]	周连童, 黄荣辉. 2008. 中国西北干旱、半干旱区感热的年代际变化特征及其与中国夏季降水的关系. 大气科学, 32(6):1276-1288.


地球物理学报 2016, Vol. 59 Issue (1): 79-92	PDF