1 引言

沙漠约占陆地表面积的14.2%，是地球上面积最大的陆地系统^[1].沙漠地表反照率大，土壤热容量小，含水量低，是地球系统中重要的感热源，对全球和区域能量平衡及气候变化具有重要的作用^[2].改进沙漠地区陆面参数的获取方法，完善陆面过程模式，提高地气间物质和能量的模拟精度，是提高地球系统模式模拟和预测能力的重要途径^[3-4].

沙漠下垫面地表以及地气能量交换特征与全球其它区域相比有着明显的差异^[5].由于降水量极低，植被稀疏，潜热通常可以忽略不计，感热输送是沙漠地区地气间能量交换的主要形式之一^[6]，也是驱动局地大气环流的主要能量来源^[7].感热与地表反照率、比辐射率、地表粗糙度和土壤热力学性质关系密切.地表反照率是陆面过程的一个重要参数，其值的变化能引起地表能量平衡的显著改变，进而影响地球气候系统.Garratt^[8]利用大气环流模式模拟发现，全球平均地表反照率增加0.13，地面日均降水减少1 mm.在植被稀少的干旱区，地表反照率对降水的影响更为显著^[9]，区域地表反照率的突然增加可能是沙漠形成的主要原因^[10].粗糙度的大小直接决定了地气间能量物质交换的强弱^{[8, 11]}.对全球不同地区地气耦合的研究发现，地表粗糙度的减小可使局地降水量明显减少^[12-14].此外，土壤作为地球系统中重要的热量储库，其热力学性质影响着地表对热量的存储和释放.Richard等^[15]研究发现，土壤热传导系数和土壤热容量是决定土壤热通量和土壤温度变化特征的主要因素.Peters-Lidard等^[16]指出，土壤热性质参数还会通过影响地表能量平衡进一步影响感热通量.

近年来，我国在西北干旱半干旱区进行了若干大型野外观测试验^{[5, 17]}，在获取地表参数、建立陆面过程模式的研究中取得了很多成果.张强等^[18-19]利用敦煌实验观测资料，运用统计分析方法，研究了获取戈壁下垫面陆面过程参数的方法.张宏升等^[20]使用一致的湍流数据采集和处理方法，对比研究了戈壁、草原等多种下垫面的湍流统计特征.Yang等^[21]利用青藏高原西北荒漠区、黑河沙漠及戈壁区域的观测数据，给出了一种计算裸土地表粗糙度的新方法.刘和平等^[22]在简单生物圈模式(Simple Biosphere Model，SiB)的基础上，改进土壤水热计算方法，提出了一个新的干旱区陆面过程模式.在此基础上，刘树华等^[23-25]发展了一套植被-土壤-大气耦合模式，较好地模拟了沙漠以及绿洲地表温度变化及能量交换特征，准确把握和描述了“冷岛效应”和“湿岛效应”.

本文利用中国西北干旱区陆气相互作用观测试验(NWC-ALIEX)巴丹吉林沙漠观测站的数据，研究了沙漠下垫面陆面过程参数化方法，建立了适合沙漠地区的陆面过程模式DLSM.该研究结果对提高沙漠地区地气交换物理过程特征的认识有一定的科学意义，为发展更为全面的陆面过程模式奠定了基础.

2 巴丹吉林沙漠观测资料

巴丹吉林沙漠观测站位于巴丹吉林沙漠边缘，102.37°E，39.47°N.观测站海拔高度1418 m，地面平坦均一，表面为裸露沙土.主要观测仪器为开路涡动协方差(OPEC)快速响应探测系统，观测项目主要有风、温、湿度脉动，四分量辐射通量，多层土壤温度(见表 1)，数据记录时间间隔为1800 s.2009年8月13日至16日，无显著天气系统过站，观测资料完整、连续.本文选用该期间观测数据，作为沙漠地区夏季典型天气的代表.

3 沙漠下垫面陆面过程模式

陆面过程需要准确模拟影响地表能量平衡的诸物理过程，包括辐射传输、湍流输送和土壤热传导过程.在沙漠地区，植被稀疏，可以忽略植物生化过程.此外，由于观测期间无降水，土壤含水量极低，潜热通量远小于净辐射通量，因此水分循环过程也可以不予以考虑.

3.1 辐射传输过程

短波辐射过程是沙漠地区最重要的物理过程之一，地表辐射平衡的改变可能是沙漠区域气候变化的主要原因^[10].由反照率的定义可知，地表向上的短波辐射通量R_s，u (W·m^-2)等于地表短波反照率α与天空向下短波辐射通量R_s，d(W·m^-2)的乘积：

(1)

给定地表比辐射率σ，地表长波反照率则为1-σ，结合斯蒂芬-玻耳兹曼公式(Stefan-Boltzmann Law)，可以计算地表向上的长波辐射通量R_l，u：

(2)

式中，R_l，d是天空向下的长波辐射通量(W·m^-2)，T_sfc为地表温度(K)，σ₀=5.67×10^-8(W·m^-2· K^-4)为斯蒂芬-玻耳兹曼常数.

3.2 湍流输送过程 3.2.1 动量通量和感热通量

大部分地表吸收的净辐射通量经过湍流输送作用，以感热通量的形式加热大气.地表向上输送至高度为z_atm (m)的动量通量和感热通量分别为^[26]：

(3)

(4)

(5)

式中，τ_x为径向动量通量(kg·m^-1·s^-2)，τ_y为纬向动量通量(kg·m^-1·s^-2)，H_s为地表向上感热通量(W·m^-2)，u_*和T_*分别为摩擦速度(m·s^-1)与温度(K)；ρ_atm为近地面大气密度(kg·m^-3)，C_p为空气定压热容(J·kg^-1·K^-1)；u_atm、v_atm分别为z_atm高度径向与纬向风速(m·s^-1)，T_atm为z_atm高度大气温度(K)；r_am、r_ah分别为大气动力阻抗和热力阻抗(s·m^-1).

3.2.2 Monin-Obukhov相似性理论

根据Monin-Obukhov相似性理论，摩擦速度可以用风速与大气温度的函数表示^[27]：

(6)

(7)

结合式(3)至式(7)，可得到大气阻抗r_am、r_ah的表达式：

(8)

式(8)中，κ=0.4为Karman常数，d为零平面位移(m)，z_0m为地表动力粗糙度(m)，z_0h为地表热力粗糙度(m).L为Monin-Obukhov长度(m)，定义如下：

(9)

g为重力加速度(m·s^-2).ψ_m和ψ_h分别是相似性函数φ_m、φ_h的积分形式，根据Högström^[28]的研究结果，函数形式如下：

(10)

(11)

其中，x=(1-19_zatm/L)^1/4，x₀=(1-19_z0m/L)^1/4，y=(1-11.6_zatm/L)^1/2，y0=(1-11.6_z0h/L)^1/2.

确定摩擦速度u_*是计算大气阻抗的关键.对于给定的一组风速、空气温度等大气强迫，模式首先根据前一时间步所得的摩擦速度，利用Arya^[29]的计算式确定整体理查森数(Bulk Richardson Number)：

(12)

然后利用式(13)，反算出u_*；经过上述两步的若干次循环迭代，最终得到趋于稳定的摩擦速度值.

(13)

3.2.3 地表粗糙度

在贴地层，动量的输送受到粗糙表面附近压力波动的影响，而热量输送则与分子扩散运动有关，因此热力粗糙度与动力粗糙度并不一致.本文采用Zeng ^[30]的热力粗糙度参数化方案：

(14)

式中α=0.13，v=1.5×10^-5 m²·s^-1为空气动力学黏度.

3.3 土壤热传导 3.3.1 热传导方程

土壤是陆面系统中重要的热量储库.白天太阳辐射加热地表，土壤吸收热量；夜晚，地表通过长波辐射向外释放热量.在陆面过程模式中考虑土壤中热量的再分布过程对准确模拟地表温度及地表能量平衡有重要的意义.

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

(15)

(16)

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

3.3.2 土壤格点设置及数值差分格式

模式采用完全离散化的差分方案数值积分热传导方程，土壤从上到下、由密到疏垂直分为若干层，以角标i表示土壤层数.模式中土壤格点设置如图 1所示，虚线表示各层土壤中心，深度为z_i，土壤温度T_s，i和土壤热性质参数C_i、k_i定义在该深度；实线为各层土壤边界，深度z′_i=(z_i+z_i+1)/2，土壤热通量G_i定义在各层界面上，界面处土壤温度以T′_s，i表示.

在土壤层厚度较小的情况下，土壤热通量G_i近似等于从z′_i到z_i的热通量大小，也约等于z_i+1到z′_i间的土壤热通量，即

(17)

消去土壤界面温度T′_i+1，可以得到土壤热通量的差分格式为：

(18)

根据式(17)，采用Crank-Nicolson离散化方法^[31]得到土壤热传导方程的差分形式为：

(19)

积分式(19)，即可模拟土壤温度的分布与变化情况.

3.4 地表能量平衡及模式计算流程

地表向上的净能量通量F_sfc及其对地表温度T_sfc的偏导数分别为：

(20)

(21)

由能量守恒定律有：

(22)

自此，式(1)至式(22)构成完整的闭合方程组.对于每一时刻的大气强迫，利用能量平衡式(22)和式(21)采用牛顿迭代方法计算出地表温度T_sfc，进而得到该时刻地表向上的辐射、感热通量和地表热通量.模式整体的计算流程如图 2所示.

4 陆面过程参数化

陆面过程参数反映了不同下垫面的物理属性，是准确模拟陆面物理过程的先决条件.利用现有的陆面过程观测资料，可以获取地表反照率、长波发射率、地表粗糙度、土壤热容量及土壤热扩散系数，作为陆面过程模式的输入参量.

4.1 地表反照率

地表反照率一般由表层土壤颜色、湿度决定^[18].由于沙漠地区极端干燥且观测时段无降水过程，可以不考虑土壤湿度对地表反照率的影响.图 3为2009年8月13日至16日地表反照率的日变化情况.采用加权平均方法^[18]获取地表反照率：

(23)

式中R_s，d，i为向下短波辐射的观测时间序列，α_i为地表反照率的观测时间序列.计算得到巴丹吉林沙漠加权平均地表反照率为0.273，高于黑河实验在沙漠中确定的地表反照率0.246^[32]，接近于敦煌戈壁地区的0.255^[18].

4.2 地表比辐射率

干旱地区地表比辐射率并没有一个公认大小，不同作者或模式取值均有不同.Stewart等^[33]假设土壤比辐射率σ=1，然而Verhoef等^[34]认为σ的取值范围是0.91~0.94，在NOAH陆面过程模式中，沙漠地区该值设为0.88^[35].

利用式(2)建立向上长波辐射通量与地表比辐射率的线性关系如下：

(24)

本文使用5 cm、10 cm两层土壤温度T_d5、T_d10计算地表温度T_sfc.结合式(15)和式(16)，有土壤温度扩散方程及其差分形式：

(25)

(26)

式中变量上标表示时间步，Δt为时间步长(s)，Δz=5 cm为土壤层厚度，经过等式变换，得到地表温度T_sfc的计算式如下：

(27)

在干旱半干旱区，向下长波辐射的测量误差较大^[36]，且CNR1型辐射计的测量值受到云量的影响^[37].为了保证计算准确，本文舍弃白天多云时的资料，也不考虑可能出现边界层云的夜间.

综合式(24)至式(27)，做出净长波辐射通量R_l，u -R_l，d与σ₀T_sfc⁴ -R_l，d的散点关系图 4.使用最小二乘法拟合，可以得到二者比值即地表比辐射率为0.950.

4.3 地表粗糙度

由式(6)可以得到动力粗糙度的计算式^[21]：

(28)

利用快速响应涡动相关系统观测资料，计算z_0m，并获得其概率密度分布函数(PDF，Possibility Distribution Function)见图 5.结果符合正态分布，使用最小二乘法拟合，得到ln (z_0m)的最佳估计值为-6.47，即z_0m=1.55×10^-3 m.结果接近敦煌戈壁地区的1.9×10^-3 m^[18]，低于黑河实验沙漠地区测得值2.67×10^-3 m^[38].

4.4 土壤热容量

土壤体积热容量与热传导率由土壤成分、土壤密度、土壤孔隙度和土壤含水量等性质决定，在沙漠地区，可忽略土壤水分的影响.根据热力学第一定律，有如下恒等关系：

(29)

式中，Δz为土壤层厚度(m)，ΔT_s为土壤层温度变化(K)，ΔG为上下两界面土壤热通量之差(W·m^-2)，Δt为时间间隔(s).利用5 cm与20 cm土壤热通量观测资料，做出方程(25)两边的散点关系图 6及其最小二乘拟合线，可得到巴丹吉林沙漠土壤体积热容量为1.08×106 J·m^-3·K^-1，该结果接近敦煌实验^[18]获得的戈壁地区土壤体积热容量值1.12× 106 J·m^-3·K^-1.

4.5 土壤热扩散系数

土壤热扩散系数λ (m²·s)定义为^[39]：

(30)

如图 7所示，采用正弦曲线拟合5 cm和10 cm土壤温度观测资料，得到两层土壤温度变化振幅分别为6.18 ℃和3.60 ℃.由振幅法^[40]：

(31)

计算出巴丹吉林沙漠土壤热扩散系数为3.342× 10^-7 m²·s.其中，P为基波周期(s)，z₁、z₂和A₁、A₂分别为两层土壤深度(m)和温度振幅(K).该值高于敦煌地区^[18]的2.52×10^-7 m²·s.

5 模拟结果与分析 5.1 模拟方案

模式使用平均风速、平均气温、大气压以及向下的长短波辐射作为输入量，以半小时为时间步长向前积分.模拟期间，地表参数保持不变.各层土壤及其界面深度如表 2所示.使用观测资料循环反复驱动模式50次，以最后模式中各变量的值作为初始状态.

同时使用相同的大气强迫资料及陆面过程参数驱动NOAH陆面过程模式^{[34, 38]}，并采用相同的初始化方法.输出结果与DLSM模拟结果及观测数据进行比对.

分别使用标准差(RMSE，Root Mean Square Error)和判定系数(R²，Coefficient of Determination)分析模拟误差、判定模拟质量，定义如下^[39]：

(32)

(33)

式中，X_obs表示观测值，X_sim表示模拟值，Var为方差.

5.2 模拟结果分析

图 8比较了巴丹吉林沙漠地表能量通量的观测值与两个模式的模拟值，表 3为模拟值与观测值的回归分析结果.可以看到，DLSM和NOAH模拟结果接近，均能较好地模拟出沙漠地区地表能量通量的日变化特征，其中DLSM略好于NOAH.由于两个模式使用相同的地表反照率，因此短波辐射的模拟结果完全一致，模拟值与观测值之间的判定系数为0.99，标准差为7.98 W·m^-2.长波辐射模拟结果接近，两个模式结果与观测值之间的判定系数分别为0.99和0.98，标准差分别为6.14，7.72 W·m^-2.感热通量的模拟结果略差，其判定系数分别为0.92和0.89，标准差分别为33.9，46.6 W·m^-2，模拟精度还有待进一步改进.

图 9为5 cm和10 cm两层土壤温度的测量与模拟值时间序列.两模式的模拟结果接近，随着土壤深度的增加，均模拟出了温度振幅衰减、相位滞后的特征，表明模式对土壤温度的模拟是合理的；但在5 cm深度，两个模式均低估了最低温度，有待进一步改进.通过对5 cm深度的土壤温度的回归分析，两个模式的模拟结果与观测值之间的判定系数分别为0.838和0.881，相应的标准误差分别为1.93，1.60 K；对10 cm深度土壤温度，两模式的判定系数分别为0.943和0.821，标准误差分别为0.48，0.79 K.DLSM的模拟精度略高于NOAH.

5.3 敏感性试验

定义归一化标准差NSEE (Normalized Standard Error of Estimate)为^[43]：

(34)

衡量从X_original调整至X_adjusted，时间序列X的相对变化大小.

为了确定地表参数变化对能量通量模拟结果的影响程度，将参数分别调增5%和调减5%输入DLSM进行模拟，将结果与原始结果进行比较(见表 4).我们定义参数调整后，NSEE大于2.5%为影响显著，大于1%认为是有明显影响.

可以看到，地表辐射特征参数，包括地表反照率、地表比辐射率，是影响沙漠地区陆面过程的主要参数.增加地表反照率，地面反射的太阳辐射显著增加，净辐射通量随之减小，而感热通量减小明显.增大地表比辐射率，向上长波辐射明显增加，导致地表温度降低，地气温差减小，感热通量显著减小.另外，由于贴地层动量与感热通过不同的机制输送，所以地表动力学粗糙度对地表能量平衡没有明显影响.在所有敏感性试验中，地表反照率变化所引起的能量平衡的改变最为显著，这是影响沙漠地区地表能量平衡最重要的参数.

6 结论

本文根据沙漠地区气候环境特点，分析和确定了地表反照率等关键参数的参数化方案；基于地表能量平衡方程，建立陆面过程模式；最后，利用巴丹吉林沙漠观测资料，模拟了该地区地表能量平衡的日变化过程.主要结论如下：

(1) 辐射传输过程、感热输送过程和土壤热传导过程是影响沙漠地区地表能量平衡的主要物理过程，通过对这三种过程的合理参数化与准确模拟，可以得到与国际知名的陆面过程模式相当的模拟结果.

(2) 地表反照率是沙漠地区最重要的陆面过程参数，其变化显著改变了地表能量平衡.地表比辐射率对地表能量平衡也有明显的影响.

(3) 巴丹吉林沙漠地表反照率为0.273，比辐射率为0.950，地表粗糙度为1.55×10^-3 m，土壤热容量和热扩散系数分别为1.08×106 J·m^-3·K^-1和3.34×10^-7 m²·s，地表参数接近于敦煌等干旱地区地表参数值.

致谢

巴丹吉林沙漠观测资料为中国科学院寒区旱区环境与工程研究所大气研究室候旭宏高级工程师提供，深表感谢.