Modeling sensible flux and latent flux in Heihe and boreal forests based on a 3D ENVI-met model

ENVI-met是由德国美因兹大学地理研究所的MICHAEL BRUSE等开发的三维微气候软件，能通过流体力学和热力学计算模拟小尺度空间内地表-植被-大气之间的相互作用。模拟环境气候参数和热量流动网格可细化到0.50 m，时间分辨率为1~10 s。ENVI-met是通过三维实体建立虚拟现实模型，子模型包括三维主模型、土壤模型及一维边界模型，其中一维边界模型会将三维模型边界按比例扩展到2 500 m的大气边界层来保证模拟的准确性^[15]。在ENVI-met 4中采用强迫边界条件功能可以让用户实时输入背景气象参数。除此之外，植被也细分为树种不同的具有3D结构的模型，植物根部的结构也进行三维模拟，细化下垫面和微环境的能量流动过程^[16]。ENVI-met在通量模拟上的主要原理如下。

下行短波辐射和长波辐射通量通过二流估算模型和一些经验公式结合计算^[14]。在模型中植被辐射吸收模型采用了引入修正系数σ来拟合植被冠层能量传输结构：

以上4个系数依次描述了植被冠层短波辐射直射、短波辐射散射、向下长波辐射和向下长波辐射在冠层不同高度基于叶面积指数（LAI）的修正系数。式（1）~式（4）中：F为消光系数，L及L^*是植被垂直结构中地表（z=0）和冠层顶部（z_p）的叶面积指数。ENVI-met是建筑热环境软件，其引入了天空可视系数σ_svf，由于本研究只考察森林环境，不考虑建筑物的遮挡，所以令σ_svf=1，由此可推导出太阳长短波辐射通量计算公式：

式（5）~式（7）中：$R_{\text{sw, dir}}^{0}$为模型顶部接收的短波直接辐射，$R_{\text{sw, dir}}^{0}$为模型顶部接收的短波反射辐射。${{\varepsilon }_{\text{f}}}$和${{\varepsilon }_{\text{S}}}$为叶片和地表的发射率，$T_{\text{f+}}^{4}$，$T_{\text{f-}}^{4}$和$T_{0}^{4}$分别为叶面上下表面平均温度及地表温度。

地表能量平衡方程是主要的控制方程。模型地表瞬时湍流：

式（8）~式（10）中：k=1时表明与地表相连的第一层，$K_{\text{h}}^{0}$和$K_{\text{q}}^{0}$为热湿传递系数，q₀为土壤第一层的水分含量。R_sw.net及R_lw.net分别为短波及长波净辐射，c_p为空气比热，L为气化潜热。J_h⁰及J_q⁰为地表热流及水分流量，u^*表示正切摩擦速度，${{\theta }^{*}}$为热差，q^*为水分含量，因此${{c}_{\text{p}}}\rho J_{\text{h}}^{0}$可表示为地表显热通量H，$\rho LJ_{\text{q}}^{0}$为地表潜热通量L_E，G为土壤热通量。

式（11）中：E_BR为地表能量闭合率，H_i为显热通量，L_Ei为潜热通量，R_ni为净辐射，G_i为土壤热通量，t为时间。

大野口关滩森林站（以下简称关滩站）位于黑河上游大野口子流域关滩阴坡的森林内（38.53°N，100.46°E），海拔高度为2 835.20 m。该地属于青藏高原垂直结构气候，年平均气温为2.8~7.6 ℃，年降水量为337.1 mm。林内主要是树高15.00~20.00 m的青海云杉Picea crassifolia均匀分布，密度1 000株·hm^-2（图 1A）。土壤以沙土为主，地面覆盖有苔藓，表层土壤含水量达66.24%。植被生长良好。涡动相关系统为Campbell开路式涡动相关通量观测系统附加HMP45C温湿度传感器，架高为20.25 m，采集频率为10 Hz。利用EdiRe软件对原始湍流数据进行质量控制和WPL修正后作为验证数据^[17]。模型模拟场景（图 1B）设置网格尺寸为1 m × 1 m × 1 m，模拟区域场景大小为30 m × 30 m。由于ENVI-met自然地表类型中只有壤土和沙土，不符合现实地表，所以综合考虑选择场景内随机分布V（沙土）：V（壤土）=3：1，更接近真实状态。植被根据样地立木平均形态选择2 m × 2 m网格BS树木模型（树高20.00 m，树冠密集，冠层分明）。

Figure 1.  GUANTAN experiment area image (A) and ENVI-met site scene modeling (B)

北方森林占全球森林总面积的29%，是世界上最大的生物群落^[18-19]。本研究用于检验的数据来源于ORNL DAAC加拿大BOREAS项目南部研究区09号通量塔（以下简称北方森林）（53.99°N，105.12°W）^[20]，海拔为600.80 m，年平均气温为1.0~5.3 ℃。该试验区主要树种为10.00 m高60年生的黑云杉Picea mariana，林木均匀分布（图 2A），冠层高度10.00~11.00 m，密度为500~600株·hm^-2。土壤类型以沙土为主，地表覆盖物为稀疏低矮草本。模型设置模型场景网格和大小同关滩站（图 2B），土壤设置为沙土。ENVI-met中在Biomet界面可以对不同树种的形态进行调整，比如树高、树冠和不同层的叶面积指数等。加拿大地区黑云杉的树高较青海云杉低，所以在Biomet中对T₁树种（树高10.00 m，树冠密集，叶少）进行编辑后采用2 m × 2 m网格放设，50株均匀分布在场景内。通量塔数据经验证符合地表能量守恒^[21]。

Figure 2.  Boreal forest experiment area image (A) and ENVI-met site scene modeling (B)

关滩森林站微气象塔高为24 m，在2，10和24 m高度上分别布设气象传感器，北部6 m处布设土壤温湿传感器^[22]。其中通量塔数据由于仪器运行条件的不足导致部分数据缺失，采用查找表法根据已知的气象数据对通量数据进行插补^[23]。测试点分别选择2010年春季5月1-10日、夏季7月21-30日、秋季9月20-29日共30 d数据进行模拟验证，其中5月4日、7月28日、9月24日有降水。图 3为关滩气象站具体气象数据。

Figure 3.  Air temperature and relative humidity from the 2 m Guantan forest meteorological station

模拟起始时间为当日0：00，共模拟24 h，模拟间隔为30 min，数据输出时间为1 h。模型输入的2 m处气温和相对湿度根据气象站数据设置日变化，风向按当日发生频率最高值。以关滩森林站2010年7月27日为例具体参数输入如表 1。

随机分布、均匀分布和聚集分布为典型的3种森林水平分布形式，为探究森林极端空间分布模式对区域总通量的影响和模型对水平异质性的敏感程度，本研究采用均匀分布和团状分布2种极端的空间分布模式来模拟夏季森林显热通量和潜热通量。场景土壤设置为壤土，植被采用自定义树木模型ZY（树高为10.00 m，冠幅3.00 m，根系深度2.86 m，根系总直径1.00 m），如图 4A所示。场景内设立90株，分布状态分别为均匀分布（图 4B）和团状分布（图 4C）。

Figure 4.  3D single tree view in Albero model with crown and root (A), uniform distributions (B), clumped distribution (C)

能量平衡闭合是地气交换的基本原理。根据能量和物质平衡基本定律，地表可用能量，即净辐射与土壤热通量之差（包括植被冠层热储存量），应等于湍流通量^[24]。但是，由于自然界中能量耗散，且测量手段限制，存在能量不闭合问题。本研究利用三维模型连续模拟2010年7月20日至29日近地面能量闭合情况，探究模型模拟能量闭合能力，设置近地面可利用能量为A，湍流通量为B。

式（12）和式（13）中：R_n为净辐射，G为土壤热通量，H为显热通量，L_E为潜热通量。

为了表明模型对实测值的反映情况，通过有效性验证（validation）对观测数据和模型预测数据进行定量比较^[25]。目前有效性验证方法很多，本研究结合通量数据观测特点，选择拟合度验证法进行模型模拟能力评估：①回归分析：应用线性回归法建立实测值和模拟值的回归方程（假设模拟值和实测值为独立变量），利用相关系数R²判断线性相关程度，并对回归系数进行假设检验判断有无统计学意义。②偏差分析：通常采用均方根误差（R_MSE）作为模拟数据和实测数据偏差的检验标准。均方根误差的计算公式如下：

式（14）中：R_MSE为均方根误差，x_i为实测值，y_i为模拟值，n为样本量，y为模拟平均值。③模型精度：为验证模型模拟精度，采用Nash-Sutcliffe效率系数^[26]（Nash-Sutcliffe coefficient of efficiency, N_S）评判模拟值与观测值之间的方差与观测方差的差异程度，如果观测值与实测值较为接近，且N_S值趋于1，模拟精度越高。计算公式为：

参考文献[3]利用敏感度绝对值验证模型敏感性的方法，对ENVI-met模型输入参数进行敏感性分析。气象参数：土壤湿度、土壤温度、2 m处空气相对湿度、风速、风向，由于ENVI-met模型对于下垫面的组成结构需要较为准确的输入，所以对土壤沙壤比也进行敏感性分析，考察场景参数对模型的影响程度。将5个参数选为敏感性参数，变化范围为-50%~50%，对显热通量和潜热通量进行模拟。并以参照日11：00通量模拟结果与改变参数后模拟结果11：00的通量值进行计算，分别得到显热通量的敏感度绝对值（absolute value of sensitive）和潜热通量的敏感度绝对值。

不同生长季通量模拟值与实测值的线性回归分析如图 5所示。实测值来看，显热通量表现出春季最大秋季最小，日平均值分别为103.15和45.84 W·m^-2；潜热通量在夏季最大，秋季最小，日平均值分别为83.10和33.12 W·m^-2；太阳短波辐射夏季最大。从模拟效果看，模型模拟结果与实测值基本一致，显热通量日平均值从高到底分别为春季（109.73 W·m^-2），夏季（78.35 W·m^-2），秋季（70.53 W·m^-2）。对比拟合线与参考线，显热通量春季和夏季拟合方程与1：1参考线较为接近，但高估秋季显热通量。这可能是由于春季近地面可利用能量耗散的主要形式是显热通量。

Figure 5.  Comparison between the observed and stimulated data by ENVI-met at Guantan station in 2010

潜热通量日平均模拟值也表现出夏季最大（95.79 W·m^-2）秋季最小（42.10 W·m^-2）的特点，这与实测情况一致。在夏季由于降水多，气温高，空气湿度增大，水汽蒸散作用强烈，潜热通量高于春季和秋季。太阳短波辐射是地面净辐射的主要来源，模型模拟值和实测值较吻合，其中夏季太阳辐射最大，实测日平均值可达702.00 W·m^-2，模拟值为577.00 W·m^-2。

从模拟有效性上可见（表 2），模型在太阳短波辐射方面模拟效果最佳，决定系数为0.86，R_MSE也最高，为164.81 W·m^-2。偏差是由于模拟时段中包含阴天，且短波太阳辐射平均值较大，R_MSE对其较为敏感造成的。其次是显热通量（R²=0.71），潜热通量模拟效果最差（R²=0.56）。此外，由于季节分析只剔除了降雨日数据，没有剔除多云数据，模型对太阳短波辐射模拟效果具有一定的低估。从N_S值来看，模型对显热通量、潜热通量和太阳短波辐射的模拟精度分别为0.68，0.63和0.80，均为正值，能够解释60%以上的原始数据。由此表明，模型对能量的模拟能力较强，平均精度可达70%。

由图 6可见：ENVI-met能较好地再现具有不同气候条件和土壤类型的地区显热通量日变化特点。关滩森林显热通量日平均实测值（模拟值）为90.18（82.48）W·m^-2（图 6A），日平均潜热通量实测值（模拟值）为68.97（64.54）W·m^-2。关滩站通量实测值在中午（10：00-15：00）并未出现较为明显的峰值，数据波动可能是由于瞬时风或仪器不稳定造成，但模型能够表现出其平均趋势。下午时段（15：00-19：00）模型低估潜热通量44.40%。通过输出模型土壤层数据（未列出），发现土壤含水量在下午时段较低，导致蒸腾作用减少。而关滩森林林下植被丰富，根系密集对土壤保水能力强，持续蒸腾时间长。这说明需细化对地表覆盖物模拟。

Figure 6.  Comparison of sensible heat flux (H) and latent heat flux (L_E) between the observed and stimulated data by ENVI-met at boreal and Guantan station

北方森林湍流通量都为单峰曲线，模拟结果能准确地再现峰值时间（图 6B）。显热通量日平均实测值（模拟值）为154.19（157.75）W·m^-2，模拟日变化结果与实测值较为吻合，除正午12：00，模型对实测值稍有低估。北方森林潜热通量峰值为250.23 W·m^-2，与模拟值相差23.00 W·m^-2。日平均潜热通量实测值和模拟值分别为78.47和80.32 W·m^-2，模拟值在夜间略高于实测值20.00 W·m^-2。模型对北方森林通量估算较为准确。

由图 7可见：地面可利用能量模型在上午的上升和下降阶段模拟结果较好，符合日变化趋势，但是模型对能量A产生70.00 W·m^-2左右的高估，且正午模型模拟能量A值较实测值偏高300.00 W·m^-2左右。这是由于夜间地面长波辐射量占主要部分，模型模拟地面向上长波辐射偏高，中午能量以向下短波辐射为主，模型模拟向下短波辐射较实测值偏低200.00 W·m^-2。除此之外，由于关滩站地表除青海云杉以外还覆盖苔藓，土壤热通量所占能量比例低，夏季土壤热通量仅占净辐射的0.3%，然而ENVI-met场景布设并没有考虑地表覆盖物的情况，所以导致太阳辐射出现后土壤热通量模拟值比实测值偏高100.00 W·m^-2左右。

Figure 7.  Simulated and observed surface flux energy balance from Guantan station between 20 July and 29 July in 2010

模型对显潜热通量之和模拟情况较好（图 7B），基本能够模拟出湍流通量日变化特征，模拟相关系数R²为0.90，R_MSE为83.58 W·m^-2。模拟值与实测值基本一致，但是模型在中午高估实测湍流通量值，是由于模型高估了潜热通量导致误差。两者日差值最大在7月21日，这是由降雨造成的。通过计算关滩站夏季实测值地表能量闭合率E_BR为66.4%，模拟值E_BR为92.9%。ENVI-met模型能量闭合情况较好但不完全闭合，可能是模型没有考虑土壤浅层热量储存问题。

图 8可以看出：模型的通量模拟中对气象参数土壤温度、2 m处空气相对湿度、风速、风向在-50%~50%区间变化的敏感度绝对值都低于1，以上4种参数对模型通量模拟影响不大。但是土壤湿度和场景中下垫面沙土壤土的比例对通量模拟的影响较大，特别是土壤湿度在30%时，模型模拟显热通量和潜热通量都对其敏感度最高，分别为2.88和3.88。总体上在土壤沙壤比与土壤湿度变化范围一致的情况下，模型潜热通量的敏感度范围较显热通量大，认为潜热通量较显热通量对土壤湿度和沙壤比更为敏感。

Figure 8.  Variation of sensitivity of ENVI-met parameters

不同林木空间分布格局下近地面太阳净辐射、显热通量、潜热通量模拟结果如图 9所示（7月24日14：00）。辐射异质性是定量描述水平尺度湍流斑块的关键，它直接影响区域能量流动的变化幅度^[9]。太阳净辐射由于树冠的遮挡散射在均匀分布林内地面（红色方框内）形成辐射值不同的光斑，且在东侧形成宽度为2 m的阴影，林内平均太阳净辐射通量为375.14 W·m^-2，辐射透过率为40%。森林团状分布下林木排列均匀无间隙，林内没有太阳辐射直接进入，而是通过树冠折射和散射到达地面，平均太阳净辐射通量为182.19 W·m^-2，辐射通过率为20%。

Figure 9.  ENVI-met energy flux simulation results at 14:00

显热通量方面，林内显热通量水平分布状态与辐射分布有较好的一致性，均匀场景内斑块状分布，团状场景内显热通量较为均匀，且为场景内最低。随着太阳高度角变化，树冠阴影在林缘区域14：00偏东。阴影遮挡部分的显热通量值低于西面地表暴露区域。森林团状分布林缘显热通量较林内高，但明显低于无立木区域。西面显热通量低于东面，且由内而外逐渐增大。可以看出并不是由于太阳辐射直接导致的，原因一方面可能是上午西面为树冠阴影，空气温度低，升温需要一段时间。另一方面根系对土壤水分有固定作用，且根系分布由林分内到外逐渐减少，保水能力也依次减弱，所以地表显热通量的敏感程度有一定差异。

潜热通量方面，森林均匀分布条件下林内水分状态良好，空气湿度大，有太阳辐射进入的像元潜热通量值较高，林缘东面由于上午受到太阳直射蒸腾大，地表湿度低，且下午有阴影遮挡缺少蒸腾动力，在场景中潜热通量值最低。西面受到太阳直射，蒸腾强烈潜热通量最高。团状分布下林内太阳辐射小，近地面潜热通量小，林缘阴影区域潜热通量小于光照区域，西面也有由内而外潜热通量增大趋势，可能是由于植被根系对土壤水分蒸腾的拦截作用使靠近林内的地表水分高于外侧，太阳辐射使水分状态较好的区域蒸腾大。

模型对不同森林空间分布格局显热通量与潜热通量的模拟结果如图 10所示。由于太阳辐射是引起能量变化的主要因素，所以只取日照时段分析。可以看出：均匀分布和团状分布在能量峰值方面造成很大差异，团状分布下显热通量峰值高于均匀分布155.52 W·m^-2，日变化幅度分别为254.33和100.19 W·m^-2。潜热通量团状分布峰值与均匀分布相差110.58 W·m^-2，日变化幅度分别为388.03和276.56 W·m^-2。团状分布下显热通量和潜热通量峰值都高于均匀分布，日变化幅度也大于均匀分布。从日变化趋势来看，团状分布对太阳辐射引起的通量日变化反映更剧烈，随时间变化，能量上升和下降过程迅速。均匀分布下能量日变化幅度较为平缓。

Figure 10.  Simulated result of sensible heat flux (H) and latent heat flux (L_E) in different distribution pattern

陆-气能量交换是研究区域生态环境乃至气候变化的基础。陆面过程模式为在地表类型一致情况下模拟通量交换的首选，但是由于植被形态各异，空间分布不均等使陆面过程模式的参数化过程困难^[27]。本研究利用小气候三维模型ENVI-met对森林地表能量流动进行模拟，探究三维模型对森林空间结构通量的模拟能力。利用黑河关滩森林以及北方森林地区森林实测通量数据进行模拟验证，并对黑河地表通量日变化和季节变化特征进行分析。此外，还探究了空间异质对模型的影响对地表能量平衡和模型敏感性进行了分析。主要结论如下：①ENVI-met模型在春季对显热通量模拟最佳，在夏季模拟潜热通量的效果较好，表明模型在不同季节对占主导地位的能量解释能力较强，这一点与蔡福等^[27]的研究结果一致。在模拟短波辐射方面，模拟值对实测值解释能力为80%。由模型在不同季节模拟结果可以看出，模型在模拟显热通量与太阳辐射的能力较强，且可以很好地反映季节变化规律。②模型可以较好地模拟出北方森林和关滩森林显热通量和潜热通量的日变化趋势。但是在关滩森林下午时段潜热通量模拟值偏低，有可能是模拟潜热通量峰值提前。有研究表明，在净辐射一定情况下，潜热通量与地表气温、地表湿度、相对湿度显著相关，峰值出现的时间也与深层水含量有极大关系^[28-29]。所以这种结果可能是因为，一方面，模型中添加土壤类型的模块尚未开发完整，导致非单一土壤类型区域的通量模拟有较大误差。另一方面，模拟关滩林下植被较为简单，与土壤之间的水分流动不够细化造成土壤持水力弱，上午蒸腾强烈导致下午土壤相对缺水，蒸腾弱。由此可见，三维模型虽然对通量具有一定的模拟能力，但是其三维模型库需要较为多种的，且具有精确参数植被。③在湍流能量（H+L_E）的模拟方面模拟值接近于实测值，但是在可利用能量方面模型造成低估，模拟地表能量闭合率E_BR为92.9%。若不考虑原数据实测值闭合率较低的情况，模型可以较好地估算地表可用能量（R_n-G）。LIU等^[30]指出显热通量与潜热通量分配比例受地表植被特点的影响较大。综上所述，模型在显热通量和潜热通量方面造成的误差是由于模型在模拟复杂自然场景，特别是复杂地表的能力有限，造成净辐射分配不合理造成的。④对ENVI-met模型通量模拟敏感度分析，发现模型对地表类型和土壤湿度较为敏感，而对风速、风向、空气湿度、土壤温度不敏感。这与LARSEN等^[31]的研究结果一致，地表土壤类型决定了地下可利用水量，而这是能量传输的关键。可以看出在通量模拟方面，准确的地表参数对模型参数输入极为重要。⑤森林分布方式不同对通量模拟有较大影响，实质上是由于不同的空间分布方式造成植被在地面投影面积不同，导致地表利用太阳辐射能量的方式不同^{[7, 32]}。团状聚集分布由于在地面造成阴影面积比均匀分布情况小，地表升温速度相对快，显热通量对空间异质较为敏感。而水平空间异质性对潜热通量在土壤水分条件较好的森林影响不大^[7]。由此可见，三维模型可以更为准确地模拟包含空间信息的地表通量。

ENVI-met模型不仅能够定性的描述森林分布状态对地表能量交换方式的影响，也能够通过像元值输出定量地表现森林空间异质性在森林通量方面的差异。引入具有较强空间三维模拟能力的小气候三维模型可以改善空间异质性对实验结果的误差，对具有不同结构的空间都能进行较为合理的参数化。将来还需进一步加强对其他不同森林类型和站点的模拟分析，定量评价水平异质性带来的误差。