气象学报  2013, Vol. 71 Issue (4): 677-691   PDF    
http://dx.doi.org/10.11676/qxxb2013.059
中国气象学会主办。
0

文章信息

王寅钧, 周明煜, 徐祥德, 孙绩华. 2013.
WANG Yinjun, ZHOU Mingyu, XU Xiangde, SUN Jihua. 2013.
MM5和ETA相似理论近地层方案对农田下垫面通量模拟比较研究
Comparative study of the similarity surface layer schemes simulate turbulent flux simulations over cropland between MM5 and ETA
气象学报, 71(4): 677-691
Acta Meteorologica Sinica, 71(4): 677-691.
http://dx.doi.org/10.11676/qxxb2013.059

文章历史

收稿日期:2012-12-06
改回日期:2013-03-18
MM5和ETA相似理论近地层方案对农田下垫面通量模拟比较研究
王寅钧1,2, 周明煜3, 徐祥德2, 孙绩华4    
1. 南京信息工程大学, 南京, 210044;
2. 中国气象科学研究院灾害天气国家重点实验室, 北京, 100081;
3. 国家海洋环境预报中心, 北京, 100081;
4. 云南省气象科学研究所, 昆明, 650034
摘要:近地层湍流通量计算对于中尺度数值模式有重要意义,湍流通量的参数化是当前大气边界层研究的重要课题之一.选择青藏高原东缘大理观象台边界层通量观测系统,离线测试了WRF区域模式中的两种常用的近地层参数化方案(MM5相似理论非迭代方案A和ETA相似理论迭代方案B),并将参数化方案计算结果与边界层铁塔涡动相关法的观测值进行对比分析.在大理观象台观测场不同植被随季节交替的状况下,根据边界层铁塔4层高度风速拟合,发现近地层空气动力学粗糙度随季节变化特征明显.将拟合的空气动力学粗糙度输入模式参数化方案进行通量计算.结果表明:稳定度是影响近地层参数化方案精度的重要因素,在不稳定条件下方案B低估了动量通量,方案A优于方案B,而在稳定条件下方案A低估了动量通量,方案B优于方案A,两种方案总体来看误差不大.对于大理边界层通量观测场农田植被交替的环境条件,不同季节下垫面植被类型的差异,以及植被的稀疏对近地层参数化方案湍流通量计算结果的精度有显着影响.方案B考虑了空气动力学粗糙度z0和热量粗糙度z0h的差异,不稳定条件下感热通量计算结果在裸土或稀少植被条件下明显优于方案A.针对方案B不稳定条件下感热通量计算结果在裸土下垫面仍出现高估的现象,使用了Zeng等1998年提出的用辐射地表温度订正裸土下垫面感热能量方法后,计算结果也有明显改善.
关键词WRF模式     近地层湍流通量     参数化     大理边界层观测资料    
Comparative study of the similarity surface layer schemes simulate turbulent flux simulations over cropland between MM5 and ETA
WANG Yinjun1,2, ZHOU Mingyu3, XU Xiangde2, SUN Jihua4    
1. Nanjing University of Inforrnatiom Science &Technology Nanjing 210044, China;
2. Chinese Academy of Meteorological Sczences (CAMS), Beijing 100081, China;
3. Natiomal Marine Erivironinental Forcastimg Center, Beijing 100081, China;
4. Meteorological Institute of Yunnan Province, Kunrning 650034, China
Abstract:Calculation of surface layer turbulent fluxes is very important for atmosphere numerical models.How to parameterizc the turbulent fluxes is one of the key research questions in current atmosphere boundary layer study.The paper uses two common schemes (MM5 similarity non-iterative Scheme A and ETA similarity iterative Scheme B) in the Weather Research Forecast Model (WRF) to make offline test and intercomparison of the parameterization results with the PBI.eddy-correlation obscrvation.The aerodynamic roughness length (z0)obtained on the Mali boundary layer is determined before calculating the turbulent fluxes.The aerodynamic roughness length (z0)by fitting four different heights wind speed from the Ptil.tower data under neutral condition varies significantly with season due to the obvious changes in underlaying surface conditions during the whole year (horsebcan in winter half year and paddy in summer half year).The results show that the vegetation sparseness degrec has a significant effect on turbulent fluxes calculation in the Mali cropland.In the unstable condition,Schcmc B undcrcstimates the momentum flux, and Scheme A is superior to Scheme B.Conversely, in the stable condition,Scheme A generates lower values with Scheme B is superior to Scheme A.On the whole,both of them produce little error.Scheme B takes account of the difference between aerodynamic roughness length (z0)and heat roughness length (z0).hhe sensible heat flux result of Scheme B is much better than that of Scheme A,especially for the bare soil or sparse vegetation underlaying in the unstable condition.In case of bare soil the calculated sensible heat flux by Scheme B is still larger than the observation.The result is much better, if the (Geng, et al,1998)'s method using the radiometric surface temperature is adopted.
Key words: WRF model     Surface layer turbulent flux     Parameterization scheme     PBL tower observation    
1 引 言

研究青藏高原近地面层的气象和湍流通量特征一直是青藏高原气象学研究的重大课题之一。第2次青藏高原大气科学试验(TIPEX)于1998年5—8月进行,在改则、当雄、昌都建立了边界层综合观测基地,合作的日本科学家在那曲建立了边界层综合观测基地。在边界层科学试验实施过程中,中国气象科学研究院、中国科学院大气物理研究所、中国国家海洋局环境预报中心、北京大学等10多个研究与业务单位及中国有关省、区气象局投入了各类先进的观测仪器装备(如:风廓线仪、6—7层(20 m)梯度观测塔、鲍恩比系统、超声探测仪(风和温度)、脉动温湿仪、红外辐射温度计、多普勒声雷达、系留气艇、低空探空仪、光学雨量计、各种辐射仪(短波、长波、直接、散射)),成功获取了大量宝贵的高原腹地横贯东西3站(改则、当雄、昌都)大气边界层资料和大量探空、地面与辐射等加密观测资料(徐祥德等,2001)。但是,由于各方面的条件限制,边界层站点主要设置在高原腹地,并且,资料的时间较短。为了弥补边界层观测样本的不足,近年来在中国科学技术部的支持下,实施了中日政府间国际合作“中日气象灾害合作研究中心”项目(简称为JICA,Japan International Cooperation Agency)计划。为增强对“世界屋脊”青藏高原及东部大气三维“立体”气象探测长期综合监测能力,中国气象局联合中国科学院等部门积极推进青藏高原及周边区域新一代气象综合观测系统的建设(Xu et al,2008徐祥德等,2009)。JICA项目在青藏高原东缘大理、林芝、理塘建立了能够进行长期连续观测的边界层探测站,积累了大量青藏高原及其东部周边地区大气边界层的观测资料。

地表与大气的能量、水分交换是影响天气过程和导致气候变化的重要因素,在大气数值模式中用近地层参数化方案来描述这一交换过程。尽管众多学者已经对近地层参数化做了较深入的研究,但是不同学者提出的近地层参数化方法仍有明显的差异,主要表现为以下几点:(1)是否需要循环迭代来计算湍流通量。大致分为迭代方案(Paulson,1970;Businger et al,1971;Dyer,1974;Holtslag et al,1988;Beljaars et al,1991;Lobocki,1993; Högström,1996)和非迭代方案(Louis,1979; Launiainen,1995;Wang et al,2002; Li et al,2010)。迭代方案通过循环方法得到奥布霍夫长度,精度较高,但计算耗时较长;非迭代方案通过近似方法,无需循环即可计算通量,计算效率高,但这要牺牲一定的精度,两者各有优缺点。(2)是否区分空气动力学粗糙度(z0)和热力学粗糙度(z0h),区分z0z0h之后,如何建立两者之间的关系,也是目前近地层湍流通量参数化的热点问题之一。不区分Z0Z0h虽然简化了计算方法,但这样计算感热通量时常出现较大误差,因此,很多学者提出将体积传输方程中计算湍流通量的关键参数Z0Z0h区分开来,这样计算的感热通量更符合实际。(3)不同近地层参数化方案中的动量通量和热量通量的普适函数表达式不同,很多学者提出的普适函数表达式是基于某一地区的观测资料得到的,如经典的Businger等(1971)Dyer(1974)Holtslag等(1988)等,所以,某一种参数化方案是否适用于其他地区需要进一步验证。

由于上述参数化方法的差异,湍流通量计算结果也必然出现明显差异,加之青藏高原东缘的地形复杂性,所以,WRF区域气候模式中常用的近地层参数化方案是否适用于青藏高原以及周边区域是很重要的问题。过去由于缺少青藏高原复杂大地形长期的边界层观测资料,无法验证不同方案的适用性。许多学者利用第2次青藏高原试验的资料已经对青藏高原腹地做了比较细致的研究。刘辉志等(2000)卞林根等(2001)分别分析了青藏高原改则和昌都地区湍流通量特征以及动量输送系数(Cd)和热量输送系数(Ch)在该地区随稳定度ζ=z/L的变化规律。李英等(2009)和杨智等(2010)分析了青藏高原东坡理塘和大理地区近地层湍流特征。近年来,基于对青藏高原东缘的边界层观测数据的分析,对高原东缘的近地层通量时空分布有了更深层次的认识。在对观测认识提高的基础上,如何用于改进数值模式的参数化方案成为问题,为了解决上述难题,首先解读模式的计算方案,采用离线测试的方法,“模拟”不同方案的计算效果,是长期以来比较有效的方法。高志球等(2001)使用了苏北稻田下垫面,将扩展Louis方案计算出的湍流通量与涡动相关法的观测结果进行了比较分析。胡艳冰等(2007)离线测试了6种湍流参数化在不同海面状况下和不同稳定度条件下的适用性。以前学者测试参数化方案时所用数据的时间通常较短,难以了解同一地区由于下垫面植被变化造成的通量计算结果的差异,并且,在青藏高原以及周边地区离线测试当前主流的中尺度WRF模式中的近地层方案的工作也很少。本研究选择了MM5相似理论近地层非迭代方案(以下简称方案A)和ETA相似理论近地层迭代方案(以下简称方案B)这两种WRF模式中常用的近地层参数化方案做了重点研究。通过离线测试的方法,用青藏高原东缘大理边界层近1 a的数据将WRF模式中的近地层参数化方案计算结果与涡动相关法观测结果在不同植被条件和不同稳定度条件下进行对比分析,明确了不同植被下垫面和不同稳定度条件对通量计算的影响,在此基础上对Z0h在感热通量计算中的作用进行了讨论,以评估这两种方案,从而为改进近地层参数化方案寻找可能的途径。

2 观测站、观测仪器、采用的资料和处理方法

大理近地层通量观测系统设置在大理国家气候观象台,海拔1990.5 m,距洱海西岸不到2 km,距观测站东面的苍山山脉约4 km。观测场下垫面较为平坦,观测场内外植被11月至次年4月(干季)种植蚕豆,5—10月(湿季)种植水稻,为大范围均一的下垫面。观测场四周环境基本可代表云贵高原西部山谷盆地地形、地貌和植被特征。

观测仪器塔高20 m,分别在铁塔2、4、10、20 m高度处安装风向、风速传感器(METONE 034B)和温、湿度传感器(VAISALA HMP45C-GM),在5.08 m 处安装超声风速温度仪(CSAT3)和CO2/H2O红外气体分析仪(Li-Cor,LI7500)。三维超声风速仪用于测量3个方向脉动风速和声学虚温,采样频率10 Hz,水平、垂直风速分辨率分别为0.001和0.0005 m/s,测量精度分别为±0.004和±0.002 m/s。CO2/H2O红外气体分析仪用于测量大气中CO2和水汽密度,采样频率为10 Hz,分辨率分别为10-3 mg/m3和10-4 g/m3,测量精度分别为±0.01和±0.15 mmol/mol。基于超声风速温度仪和CO2/H2O红外气体分析仪采样频率10 Hz 的资料,采用涡动相关法计算大理观象台近地层湍流通量(感热通量、潜热通量、二氧化碳通量和动量通量)。

湍流资料被处理成30 min平均,首先剔除超合理值、存在明显错误的野点,然后进行线性去趋势和坐标旋转,最后进行平均量、脉动量、方差、协方差、通量等一系列湍流统计运算。对于湍流资料的处理应用了二次旋转修正的方法。

主要使用2008年11月—2009年9月的资料,资料质量总体较好,通过农田下垫面能量平衡闭合分析发现干季(11月—4月)闭合率达到0.83,湿季(5—9月)闭合率达到0.75,这一结果与杨智等(2010)基本一致。由通量印痕分析(Schuepp et al,1990)可知,蚕豆生长后期、裸土和水稻生长后期对通量观测值产生90%贡献的通量贡献区平均长度分别为850、1100和740 m。

3 WRF模式中两种近地层方案的计算方法3.1 MM5相似理论近地层方案

MM5相似理论近地层方案使用Paulson(1970)Dyer等(1970)给出的稳定度函数,Webb(1970)给出的方法计算近地表的动量、热量、水汽输送系数。在该方案中没有考虑Z0h,将稳定度划分为4个范围(Zhang et al,1982)。

该方案通过地表和参考高度(如10 m)的各个变量值来计算通量。需要输入的变量有:参考高度上的风速(U、V)、空气比湿(混合比)、温度和气压,以及地面的温度、气压和空气动力学粗糙度(z0)。此方案中需要对变量摩擦速度(u*)、温度尺度(T*)、稳定度参数(ξ=z/L)赋初值(u*=0.001,T*=0,ζ=0)。该方案可以输出的变量有:热力学粗糙度(z0h)、摩擦速度(u*)、动量通量(τ)、感热通量(H)、稳定度(ξ)等。

计算通量的算法步骤如下:

(1)首先根据Blackadar(1976)方法,计算总体理查森数(Ri)。

(2)由理查森数Ri判断大气层结,然后定出对应层结下的稳定度(ξ)。

(3)计算湍流通量。

3.2 ETA相似理论近地层方案

ETA相似理论近地层方案(Janjic,19962002)基于相似理论(Monin et al,1954),但还包括粘性副层参数化。在水面上,该方案明确地使用Janjic(1994)的方法。在陆面上,粘性副层的影响通过考虑对于温度和湿度的可变粗糙度来解决,这一方法由Zilitinkevich(1995)提出。在近地层不稳定和风速消失情况下,使用Beljaars(1995)订正来避免奇异值的出现。近地层通量通过迭代法计算。

该方案通过地表和参考高度(如10 m)各个变量值来计算通量。需要输入的变量有:参考高度上的风速(U、V)、空气比湿(混合比)、温度和气压,以及地面的温度、气压和空气动力学粗糙度(z0)。此方案需要对变量摩擦速度(u*)、动量通量交换系数(CM)与水汽和热量这两个通量的交换系数(CH,E)赋初值(u*=0.1,CM=Cd×u=0.01,CH,E=Ch×u=0.01)。迭代5次之后的结果,认为其精度满足需求。可以输出的变量有:摩擦速度、动量通量、感热通量、稳定度等。

4 空气动力学粗糙度的确定

用近中性条件下(梯度理查森数-0.03<Ri<0.03)梯度观测资料,采用非线性拟合计算空气动力学粗糙度,结果如图 1所示。将不同下垫面时段的中值作为参数化计算时输入的空气动力学粗糙度。

图 1 通过中性条件下风速廓线数据得到的空气动力学粗糙度(a. 2008年10月20日—2009年3月29日,b. 2009年4月5日—2009年9月13日; 圆圈为每一时段内所有样本按大小排序的中值,上下限为25%和75%处的值) Fig. 1 Estimated aerodynamic roughness lengths from the profile data under neutral conditions(a. from 20 October 2008 to 29 March 2009,and b. from 5 April 2009 to 13 September 2009; the circle represents median value,and the lower and upper limit represents respectively 25% and 75% position value)

空气动力学粗糙度在2008年10月20日—2009年3月28日是增大过程,蚕豆从播种一直到生长成熟,植被逐渐茂密,使得空气动力学粗糙度逐渐增大。10月底蚕豆播种,到11月进入三真叶和分枝阶段,此时植被稀疏,地表近似为裸土,空气动力学粗糙度较小,在0.02 m左右,由于植被生长的高度、密度增加很快,空气动力学粗糙度随时间增大速率较快,至1月以后,蚕豆生长进入中后期,1月中旬开始开花,植被生长的高度、密度增加较慢,空气动力学粗糙度随时间增大速率明显放缓,为0.1 m 左右,2月中旬开始结荚,植被生长的高度、密度基本维持不变,这段时间空气动力学粗糙度基本不变,为0.12 m左右(图 1a)。

空气动力学粗糙度在2009年4月15日—2009年9月3日是先减小后增大过程。蚕豆在4月底成熟,5月1日收割,收割前的半个月由于叶片枯萎,空气动力学粗糙度逐渐减小,收割后下垫面为裸土,空气动力学粗糙度为0.05 m左右,比图 2a中2008年11月时0.03 m要大一些,这可能与这一时段翻地有关,5月下旬水稻移栽后,空气动力学粗糙度

总体呈缓慢增大的趋势,从0.10 m左右增大到0.20 m左右(图 1b)。

图 1中7—8月水稻成熟期的空气动力学粗糙度比3—4月蚕豆成熟期大很多,两者的植被高度基本相同,都为1 m左右,产生的原因可能是: 6月以后,大理进入雨季,近地层风速减小,相对较弱的风使水稻弯伏的波动减小,覆盖层变得更为粗糙,也会使粗糙度增大。

图 1得出不同下垫面时段粗糙度的均值如表 1所示。

表 1 估算出的各下垫面的空气动力学粗糙度 Table 1 Estimated aerodynamic roughness lengths for the different underlying surfaces
下垫面类型空气动力学粗糙度(m)
蚕豆的生长初期(10—11月)0.01—0.03
蚕豆的生长后期(1—4月)0.08—0.12
裸土(4月下旬—5月上旬)0.03—0.07
水稻移栽后的生长初期(5—6月)0.08—0.12
水稻移栽后的生长后期(7—8月)0.12—0.20
5 不同下垫面通量测试结果分析5.1 生长后期蚕豆下垫面的湍流通量

在不稳定条件下,蚕豆下垫面动量通量观测值与模式计算值大小的分布总体相差不大,观测值和模式计算值的中值和上下限为25%和75%处的值都很接近,方案A高估10%左右,方案B低估10%左右。但在稳定条件下,两方案都有明显的不足,并且出现明显的差异,方案A计算值偏小,方案B计算值偏大,且层结越稳定方案B误差越明显。方案A偏小的原因可能是稳定层结模式计算值出现较多u*=0.1 m/s的点(此方案认为陆面u*最小为0.1 m/s),比实际观测值明显偏小(图 2)。

图 2 蚕豆下垫面不同稳定度条件下动量通量

(a.方案A,b.方案B; 圆圈或者叉点为中值,上下限为25%和75%处的值,黑色为观测值,灰色为WRF模式计算值)
Fig. 2 Momentum flux for the horsebean underlying surface under the different stability conditions

(a. Scheme A,and b. Scheme B; circle or cross represent median value,and the lower and upper limit represent respectively 25% and 75% position value; blue represent observed value,and red represent WRF model-calculated value)

方案B和方案A观测值与模式计算拟合的斜率分别为0.88和0.62,方案B更优,但方案B点的离散似乎更大一些(图 3)。

图 3 蚕豆下垫面动量通量观测值和WRF计算值的散点图

(a. 方案A,b. 方案B)
Fig. 3 Measured and simulated variations of the momentum flux for the horsebean underlying surface

(a. Scheme A,b. Scheme B)

蚕豆下垫面在不稳定条件下,两种方案计算的感热通量与观测值结果相符的较好,稳定条件下,考虑到蚕豆下垫面感热通量的观测值存在-30—20 W/m2的系统偏差,故方案A的高估比图中所示要小,方案B出现一定程度的低估,总体两种方案效果接近(图 4)。

图 4 蚕豆下垫面不同稳定度条件下感热通量

(a. 方案A,b. 方案B; 圆圈或者叉点为中值,上下限为25%和75%处的值,黑色为观测值,灰色为WRF模式计算值)
Fig. 4 Heat flux for the horsebean underlying surface under the different stability conditions

(a. Scheme A,b. Scheme B; circle or cross represent median value,and the lower and upper limit represent respectively 25% and 75% position value,blue represent observed value,and red represent WRF model-calculated value)

拟合的斜率方案B为0.75,而方案A为0.58,方案B更接近1,更接近观测值,效果更好。相关系数方案B为0.82,而方案A为0.88,方案A略好一些(图 5)。

图 5 蚕豆下垫面感热通量观测值和WRF计算值的散点图

(a. 方案A,b. 方案B)
Fig. 5 Measured and simulated variations of the heat flux for the horsebean underlying surface

(a. Scheme A,b. Scheme B)
5.2 生长后期水稻下垫面的湍流通量

图 6可见,与生长后期的蚕豆下垫面类似,在不稳定条件下,方案A高估,方案B低估,稳定条件下,方案A误差较小,方案B高估。从图 7散点图中可以看到,方案A在稳定条件(u*=0.1 m/s为最小值时)的个别时段出现模式计算偏小的情况,不稳定条件下计算值偏大,方案B不稳定条件下误差较小,稳定条件下高估,从拟合的斜率来看,方案A分别为1.15,方案B为0.81,相关系数方案A为0.80,方案B为0.75,总体来看,两种方案效果比较接近。

图 6 水稻下垫面不同稳定度条件下动量通量

(a. 方案A,b. 方案B; 圆圈或者叉点为中值,上下限为25%和75%处的值,黑色为观测值,灰色为WRF模式计算值)
Fig. 6 Momentum flux for the paddy underlying surface under the different stability conditions

(a. Scheme A,b. Scheme B; circle or cross represent median value; and the lower and upper limit represent respectively 25% and 75% position value,blue represent observed value,and red represent WRF model-calculated value)
图 7 水稻下垫面动量通量观测值和WRF计算值的散点图

(a. 方案A,b. 方案B)
Fig. 7 Measured and simulated variations of the momentum flux for the paddy underlying surface

(a. Scheme A,b. Scheme B)

在稳定条件下,考虑到水稻下垫面感热通量的观测值存在-20—10 W/m2的系统偏差,方案A和方案B感热通量计算值都与观测值非常接近,在不稳定条件下,方案A感热通量计算值比观测值有一定程度的偏大,方案B误差较小,效果较好(图 8)。

图 8 水稻下垫面不同稳定度条件下感热通量

(a. 方案A,b. 方案B; 其他同图 4)
Fig. 8 The heat flux for the paddy underlying surface under the different stability conditions

(a. Scheme A,b. Scheme B; the others as in Fig. 4)

方案A的拟合斜率和截距分别为1.17和19.27,拟合斜率大于1的主要原因是在感热通量为高值时(50—250 W/m2),模式计算值偏大比较明显,方案B的拟合斜率和截距分别为0.77和18.48,拟合斜率小于1的主要原因是在感热通量为低值时(-100—0 W/m2),模式计算值偏大。相关系数方案A为0.86,方案B为0.79,前者更优。两种方案计算效果总体上看比较接近(图 9)。

图 9 水稻下垫面感热通量观测值和WRF计算值的散点图

(a. 方案A,b. 方案B)
Fig. 9 Measured and simulated variations of the heat flux for the paddy underlying surface

(a. Scheme A,b. Scheme B)
5.3 裸土下垫面的湍流通量

观测值在不稳定条件下动量通量在0.08 N/m2左右,稳定条件下为0.03 N/m2左右。在不稳定条件下,方案A高估,方案B低估,稳定条件下,方案A低估,方案B高估,两种方案的误差随稳定度变化呈现相反分布的情况。两种方案模拟的这一误差规律与蚕豆下垫面有一定的相似(图 10)。

图 10 裸土下垫面不同稳定度条件下动量通量

(a. 方案A,b. 方案B; 其他说明同图 6)
Fig. 10 Momentum flux for the bare soil underlying surface under the different stability conditions

(a. Scheme A,b. Scheme B; the others as in Fig. 6)

方案A相关系数为0.71,拟合斜率为0.73,方案B相关系数为0.62,拟合斜率为0.68,方案A优于方案B(图 11)。

图 11 裸土下垫面动量通量观测值和WRF计算值的散点图

(a. 方案A,b. 方案B)
Fig. 11 Measured and simulated variations of the momentum flux for the bare soil underlying surface

(a. Scheme A,b. Scheme B)

裸土下垫面时,在稳定条件下,方案A计算的感热通量值与观测值非常接近,方案B计算的感热通量值比观测值略微偏小,但在不稳定条件下两种方案均出现模式计算值比观测值明显偏大的现象(图 12)。

图 12 裸土下垫面不同稳定度条件下感热通量

(a. 方案A,b. 方案B; 其他同图 4)
Fig. 12 Heat flux for the bare soil underlying surface under the different stability conditions

(a. Scheme A,b. Scheme B; the others as in Fig. 4)

方案A和方案B拟合斜率分别达到1.77和1.33,远大于1,方案B计算结果更接近观测。两种方案模拟结果与观测值相关系数都较高,方案A为0.94,方案B为0.89(图 13)。

图 13 裸土下垫面感热通量观测值和WRF计算值的散点图

(a.方案A,b.方案B)
Fig. 13 Measured and simulated variations of the heat flux for the bare soil underlying surface

(a.Scheme A,b.Scheme B)

由于裸土下垫面时间较短(4月下旬—5月上旬),不足一个月,为了证明方案A误差较大这一结果并非偶然,使用2008年裸土下垫面的资料,发现方案A拟合斜率为1.54,高估仍然明显(图略)。

5.4 不同植被下垫面和不同稳定度条件下误差分析

图 14中横纵坐标为标准化后的标准差f。方位角的大小代表了相关系数R。绿色虚线表示了各红色圆点到参考点的距离,为中心形式的均方根误差。对于有植被的蚕豆和水稻下垫面,方案A和B的R相差不大,但方案B的f比方案A更接近于1,计算结果与观测数据更接近,效果更好一些。对于裸土下垫面,情况有所不同,方案A计算的动量通量f更接近1,效果更好,但感热通量方案A的计算误差过大,方案B明显优于方案A。总体来看,方案B效果更好一些。

图 14 方案A和B对于不同下垫面湍流通量的泰勒统计图 Fig. 14 Talor diagram of the simulation results with A and B schemes for the different underlaying surfaces

表 2可见,方案A和方案B计算动量通量在稳定和不稳定条件下均方根误差总体而言差异不大,都在0.05 N/m2左右;方案B裸土下垫面计算感热通量在不稳定条件下均方根误差明显小于方案A,其他条件下,两者差异较小。不稳定条件下,方案A计算的动量通量效果优于方案B,方案B的结果明显低估,稳定条件下,方案B计算的动量通量效果较好,误差在10%左右,方案A结果低估较多,误差较大。

表 2 不稳定条件和稳定条件下不同下垫面的误差 Table 2 The error under unstable and stable conditions for the different underlying surfaces
下垫面类型方案动量通量均方根误差/百分比偏差感热通量均方根误差/百分比偏差
不稳定条件z/L<0(N/m2)稳定条件z/L>0(N/m2)不稳定条件z/L<0(W/m2)稳定条件z/L>0(W/m2)
蚕豆A0.071/ -7% 0.077/-48% 36.45/-21%38.41/-63%
B0.070/-12% 0.076/4%45.73/-27%36.39/-23%
裸土A0.060/-10% 0.068/-59%139.56/101% 24.23/-57%
B0.075/-43%0.059/-15% 85.20/33%36.22/17%
水稻A0.050/28%0.031/-29%51.36/35%25.76/-60%
B0.034/-21%0.035/-16% 39.25/-5%27.03/-87%
5.5 空气动力学粗糙度对湍流通量计算结果的影响

由于两种参数化方案计算湍流通量前需要输入z0,本研究采用梯度观测资料按对数律拟合的方法得到z0,但是计算z0有多种方法,不同的计算方法结果有一定的差异,所以有必要做z0对湍流通量计算结果的敏感性分析。图 15中,随着z0的增大,两种方案计算出的湍流通量数值也增大,方案A计算的动量通量增大的速度更为明显,两种方案计算出的感热通量增大的速度基本一致。根据前面z0的计算结果,裸土下垫面的z0≈0.05±0.02 m,生长后期的水稻下垫面的z0≈0.16±0.04 m。从图 15a和b都可以看出由于z0计算结果的不确定性引起了湍流通量结果也存在一定的不确定性,由此带来的湍流通量波动范围大约在20%,如果取z0等于估算范围的中间值,则由z0引起的相对误差在±10%以内。方案A计算的动量通量相对误差略小一些,但是方案B计算的感热通量远优于方案A。

图 15 空气动力学粗糙度对湍流通量计算结果的影响

(a. 裸土下垫面,b. 生长后期的水稻下垫面; 虚线为方案A,实线为方案B)
Fig. 15 Effect of aerodynamic roughness length z0 on the turbulent flux calculation

(a. bare soil underlying surface,b. paddy underlying surface; dash line represents Scheme A,solid line represents Scheme B)
6 讨论与结论

关于农田下垫面近地层湍流通量的参数化研究,许多学者已经做了大量研究,并得出很多有意义的结论。高志球等(2001)使用了苏北稻田下垫面2001年6月10日—7月20日的数据,将扩展Louis方案计算出的湍流通量与涡动相关法的观测结果进行了比较分析,发现总体来看扩展Louis方案模拟动量通量较为准确,低估感热通量10%左右,低估潜热通量11%左右。郭建侠(2006)采用Businger(1971)Dyer(1974)Louis(1982)Beljaars(1991)Uno(1995)5种方案计算华北玉米地下垫面的近地层湍流通量,与涡动相关法观测资料结果进行对比分析。结果表明: Businger(1971)方案估计的动量通量与实测值的平均偏差最小,为0.002 N/m2,而后依次为Louis(1982)方案、 Uno(1995)方案、 Beljaar(1991)方案、 Dyer(1974)方案;各参数化方案对感热通量的估计偏差在玉米生长前期和后期存在明显的不同,各方案对生长前期一致地高估感热通量,但在生长后期模拟的感热通量普遍较好,误差减小,高估和低估均可出现。本研究采用WRF模式目前使用的两种方案,与以前学者的工作略有不同,从模拟结果来看,动量通量在高值时比较离散,稳定条件下偏差较大是目前这两种方案的明显不足之处,此外,感热通量在植被稀疏或无植被时普遍高估也是各方案需要改进的。

裸土下垫面两种方案计算感热通量结果都偏大较多,并简要分析了造成这一现象的可能原因。由于无植被叶片反射太阳辐射以及植被的蒸腾作用,地面白天增温加快。另外,对于裸土下垫面,人为翻耕的土壤,使得土壤孔隙度增大,在土壤湿度没有明显变化的情况下,土壤热容量降低,这也使得地面在白天吸收太阳辐射后增温加快。地温增温快,而气温响应偏慢,近地层出现强不稳定层结,稳定度ξ=z/L负值较大,以热力湍流为主。参数化方案中的感热通量与地-气位温之差(θs-θa)成正比,裸土状态下,晴天中午θs-θa比蚕豆下垫面大2—3倍(蚕豆为5—10 K,裸土为10—30 K),而从实际涡动观测结果来看,感热通量在裸土下垫面与蚕豆下垫面中午峰值无明显差异,都在100 W/m2左右,由此可见,在风速变化不大的情况下,裸土下垫面时热量输送系数(Ch)是蚕豆下垫面的,所以计算感热通量的关键问题是得到较为准确的Ch

图 13中方案A的误差表现得更为明显,方案A不区分动力学粗糙度z0和热力学粗糙度z0h,虽然简化了计算,但这样计算的感热通量出现较大误差,显然不适合应用于裸土下垫面,方案B将z0z0h区分开来,这样计算的感热通量更符合实际。做离线测试前,z0可以通过中性条件下多个高度风速进行对数拟合得到,而z0h一般来源于参数热传输附加阻尼kB-1=ln(z0/z0h)确定,这里需要进行参数化。不同学者提出的计算kB-1的方法有所差异,但大体的思想都是建立粗糙度雷诺数Re*kB-1的关系,只是具体表达式不同,如Brutsaert(1982)Zilitinkevich(1995)等。Chen等(1997)在测试方案B性能时,发现使用Zilitinkevich订正确定z0h比假设z0hz0之间为固定比值在计算感热通量时有明显的改进。目前WRF中方案B通过假设关键参数Zilitinkevich系数(Czil=0.1)来确定z0hkB-1的大小,从

而得到相应的Ch。方案B区分z0hz0,计算的热量输送系数Ch比方案A更接近实际情况,有效解决了裸土下垫面不稳定条件下方案A热量输送系数Ch明显偏大的不足,感热通量的模拟效果有明显改善。

但最新的研究表明,系数Czil对于不同植被下垫面可能不为常数。Rosero等(2010)通过敏感性试验发现方案B中系数Czil对下垫面为裸土或稀疏植被时的感热通量计算起到关键性作用,但对茂密植被下垫面的影响相对较小。

Chen等(2009)使用美国中西部12个不同下垫面边界层站数据拟合得到了Czil与植被高度(h)的函数关系近似为Czil=10(-0.4h)。改进裸土下垫面感热模拟的一条可能的思路是参数Czil不设为默认值0.1,增大Czil,从而减小z0h,使得地表热量输送强度减弱,湍流减弱,Ch减小,最终达到改进感热通量计算的目的。

针对方案B在裸土下垫面仍然会出现感热通量计算偏大较多的现象,本研究引进Zeng等(1998)有关使用辐射地表温度时裸土下垫面订正公式

用上式订正过后,方案B计算感热值有较大改进,不稳定条件下计算值与观测值更为接近,拟合斜率下降到1.13,散点图中偏大的点数量明显减少(图 16)。

图 16 使用Zeng等(1998)订正方法后的方案B效果

(a. 裸土下垫面不同稳定度条件下感热通量,b. 裸土下垫面感热通量观测值和WRF计算值散点图)
Fig. 16 Revised Scheme B which adopts the correction method in Zeng et al.(1998)

(a. the heat flux for the bare soil underlying surface under the different stability conditions,b. measured and simulated variations of the heat flux for the bare soil underlying surface)

本研究认为在数值模式在线运行时,对于植被稀疏或裸土下垫面,计算感热通量是可采用订正后的方案B,这样可以有效减小误差。

将WRF模式中常用的两种近地层参数化方案(MM5相似理论非迭代方案A和ETA 相似理论迭代方案B)使用大理边界层铁塔资料进行了离线测试,得到以下主要结论:

(1)两种方案在不同稳定度条件下计算动量通量的效果呈现出不同程度的偏差。在不稳定条件下方案B低估动量通量,方案A优于方案B,而在稳定条件下方案A低估动量通量,方案B优于方案A,两种方案的误差随稳定度变化呈现相反分布的情况,两种方案的稳定度函数的形式都较适用于大理农田下垫面。

(2)下垫面的植被类型以及植被的稀疏对于感热通量的计算精度有明显影响。由泰勒统计图和误差分析表可知,对于有植被的蚕豆和水稻下垫面,方案B的标准化后的标准差f比方案A更接近于1,计算效果更好,相关系数R和中心形式的均方根误差方案A和方案B效果差异不大;在无植被或植被较为稀疏情况下,区分空气动力学粗糙度z0和热量粗糙度z0h的方案B在不稳定条件下效果明显优于不区分的方案A,并且植被越稀少,差异越明显。

(3)参数z0对湍流通量的计算有重要的影响。随着z0的增大,两种方案计算出的湍流通量数值也增大,方案A计算的动量通量增大的速度更为明显,两种方案计算出的感热通量增大的速度基本一致。由于估算出的z0存在一定的不确定性,由此会引起湍流通量计算结果大约±10%的误差。

(4)对于裸土下垫面,方案B虽优于方案A,但仍然出现较大误差,这与观测资料使用辐射地表温度可能有一定关系,使用Zeng等(1998)订正公式可以有效减小误差,提高计算精度。

(5)上述离线测试的研究表明,数值模式在线运行近地层参数化方案B计算感热通量时,对于植被稀疏或裸土下垫面,可采用订正后的方案B,这样可以有效减小误差。

参考文献
卞林根.2001.青藏高原东南部昌都地区近地层湍流输送的观测研究.应用气象学报,12(1):1-13
郭建侠.2006.华北玉米下垫面湍流输送特征及参数化方案比较研 究[D].南京:南京信息工程大学,93-98
胡艳冰.2007.不同海况下六种湍流通量参数化方案的对比分析. 海洋预报,24(2):9-16
李英.2009.青藏高原东坡理塘地区近地层湍流特征研究.高原气 象,28(4):745-753
刘辉志.2000.青藏高原改则地区近地层湍流特征.大气科学,24 (3):289-300
徐祥德,周明煜,陈家宜等.2001.青藏高原地-气过程动力热力结 构、综合物理图象.中国科学(D辑),31(5):428-440
徐祥德.2009.青藏高原“敏感区”对我国灾害天气气候的影响及其 监测.中国工程科学,11(10):96-107
杨智,刘劲松,朱以维等.2010.云贵高原西部大理地区近地层湍流 特征分析.大气科学学报,33(1):117-124
张人禾,徐祥德.2012.青藏高原及东缘新一代大气综合探测系统 应用平台_中日合作JICA项目.中国工程科学,14(9):101-112
周明煜,徐祥德,卞林根等.2000.青藏高原大气边界层观测分析与动力学研究.北京:气象出版社,17-25
Arya S P. 1988. Introduction to Micrometeorology. New York: Academic Press,307pp
Beljaars ACM,Holtslag A A M. 1991. Flux parameterization over land surfaces for atmospheric models. J Appl Meteor, 30 (3):327-341
Beljaars ACM. 1995. The parameterization of surface fluxes in large-scale models under free convection.Quart J Roy Meteor Soc,121 : 255-270
Blackadar A K. 1976. Modeling the nocturnal boundary layer//Proc. Third Symposium on Atmospheric Turbulence Diffnsion and Air Quality. Raleigh,Amer Meteor Soc,46-49
Brutsaert W A. 1982. Evaporation into the Atmosphere. Reidel Publishing,299pp
Businger J A,Wyngaard J C,Izumi Y,et al. 1971. Flux-profile relationships in the atmospheric surface layer. J Atmos Sci,28(2) 181-189
Chen F, Mitchell K,Schaake J. 1996. Modeling of land-surface evaporation by Lour schemes and comparison with FIFE observations J Ueophys Res,101,7251-7268
Chen F, Janjic Z I,Mitchell K. 1997. Impact of atmospheric surface-layer parameterizations in the new land-surface scheme of the NCEP mesoscale Eta Model. Bound Layer Meteor,85;391-421
Chen F,Zhang Y. 2009. On the coupling strength between the land surface and the atmosphere: From viewpoint of surface exchange coellicients. Ueophys Res Lett, 36, L10404, doi: 10. 1029/2009gl037980
Dyer AJ,Hicks B B. 1970. Flux-gradient relationships in the constant flux layer. Quart J Roy Meteor Soc,96:715-721
Dyer A J. 1974. A review of flux-profile relationships. Bound-Layer Meteor,7(3):363-372
Uao Z Q, Bian L U, Zhou X J. 2003. Measurements of turbulent transfer in the near-surface layer over a rice paddy in China. J Ueophys Res,108(D13):4387
Uarratt J R, Francey R J. 1978. Bulk characteristics of heat transfer in the unstable, baroclinic atmospheric boundary layer. Bound Layer Meteor,15:399-421
Högström U LF. 1996. Review of some basic characteristics of the atmospheric surlace layer. Bound Layer Meteor,78(3-4):215-246
Holtslag A A M,de Bruin H A R. 1988. Applied modeling of the nighttime surlace energy balance over land. J Appl Meteor, 27 (6);689-704
Janjic Z I. 1994. The step-mountain Eta coordinate model further developments of the convection viscous sublayer and turbulence closure schemes. Mon Wea Rev,122:927-945
Janjic Z I. 1996. The surlace layer in the NC}EP Eta Model//Eleventh Conlerence on Numerical Weather Prediction. Norlolk, VA,19-23 August;Amer Meteor Soc,Boston,MA,354-355
Janjic Z I. 2002. Nonsingular Implementation of the Mellor-Yamada Level 2. 5 Scheme in the NCEP Meso model. NCEP Office Note. No. 437.6100
Jiménez P A,Dudhia J,et al. 2012. A revised scheme for the WRF surlace layer formulation. Mon Wea Rev,140(3):898-918
Monin A S, A M Obukhov. 1954. Basic laws of turbulent mixing in the surlace layer of the atmosphere. Contrib Ueophys Inst Acad Sci, USSR, (151):163-187
Launiainen J. 1995. Derivation of the relationship between the Monin-Obukhov stability parameter and the bulk Richardson number for flux-prolile studies. Bound Layer Meteor,76(1-2);165-179
Lobocki L. 1993. A procedure for the derivation of surface-layer bulk relationships from simplified second-order closure models.J Appl Meteor,32;126-138
Louis J F. 1979. A parametric model of vertical eddy Iluxes in the atmosphere. Bound Layer Meteor,17(2):187-202
Louis J F,TiedtleM,Ueleyn J F. 1982. A short history of the PBL Parameterization at ECMWF. Proceedings of the 1981 ECMWF Workshop on Planetary Boundary Layer Parameterization, 59-80
Paulson C A. 1970. The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. J Appl Meteor,9:857-861
Rosero E, Yang Z L, Wagener T,et al. 2010. Quantifying parameter sensitivity, interaction, and transferability in hydrologically enhanced versions of the Noah land surface model over transition zones during the warm season. J Ueophys Res, 115(D03106):1-21
Schuepp P H, Leclerc M Y, MacPherson J I, et al. 1990. Footprint prediction of scalar Iluxes Irom analytical solution of the dillution equation. Bound Layer Meteor,50:355-373
Uno I, Cai X M, Steyn D U, et al. 1995. A simple extension of the Louis method for rough surlace layer modeling. Bound Layer Meteor,76;395-409
Xu X D,Zhang R H,Koike T,et al. 2008. A new integrated obser-vational system over the Tibetan Plateau. Amer Meteor Soc, 1492-1496
Wang S P, Wang Q, Doyle J. 2002. Some improvements of Louis surlace Ilux parameterization刀Preprints, 15th Symposium on Boundary Layer and Turbulence. Amer Meteor Soc, Wageningen,Netherlands,547-550
Webb E K. 1970. Prolile relationships: The log-linear range, and extension to strong stability. Quart J Roy Meteor Soc,96:67-90
Zeng X,R E Dickinson. 1998. EIIect of surlace sublayer on surlace skin temperature and Iluxes. J Climate,11;537-550
Zhang D L,Anthes R A. 1982. A high-resolution model of the planetary boundary layer-Eensitivity tests and comparisons with SESAME-79 data. J Appl Meteor,21:1594-1609
Zilitinkevich S S. 1995. Non-local turbulent transport:pollution dispersion aspects of coherent structure of convective fiiow//Power H,Moussiopoulos N, Brebbia C A. Air Pollution 111-Volume I. Air Pollution Theory and Simulation. Southampton Boston Computational Mechanics Publ ications, 53-60