临界理查森数对不稳定状态下边界层模拟的影响

2022-09-29韦俏雅王咏薇孙永

科学技术与工程 2022年22期

韦俏雅，王咏薇*，孙永

(1.南京信息工程大学大气物理学院，南京 210044； 2.南京信息工程大学大气环境中心，南京 210044； 3.南京大学大气科学学院，南京 210044)

大气边界层是受到下垫面影响的大气层。边界层的运动主要为湍流运动。主导着地表与大气之间动量、热量、水汽以及污染物浓度的交换[1-2]。因此，边界层的结构特性和变化规律一直是大气科学界关注的焦点，边界层内通量输送过程的精确模拟对天气气候及污染事件的预测模拟非常重要[3-5]。在数值模拟研究工作中，科学家们发现，不同的边界层参数化方案，对于天气气候以及污染物的预报存在显著的影响[6-7]。当采用同一种边界层参数化方案时，其关键参数(如卡曼常数、临界理查森数等)的设置不同，也会导致湍流运动的模拟存在明显差异，对地气交换产生显著影响，使得边界层结构出现偏差，进而影响到天气、气候、以及污染物扩散的规律发生变化[8-10]。

临界理查森数为基于Richardson[11]和Lyons等[12]的科学实验所获取的经验值，通常取值范围为0.15～1.0。在后期的研究中，研究学者们根据各自实验结果得出了不同值，而数值模式的开发者，在模式中参考不同研究者的研究结果，对理查森数的临界值选取了不同的常数值。例如，Noh等[13]在研究中提出理查森数的临界值设置为0.25；Holtslag等[14]在全球气候模型(GCM)中将理查森数的临界值设置为0.5；Kiehl等[15]和Wu等[16]设定CAM3和BCC_AGCM方案中的理查森数的临界值设置为0.3。

根据之前的研究，临界理查森数是和大气稳定度有关的一个变量。Vogelezang等[17]根据不同条件层结,将临界理查森数的临界值分别进行不同划分,弱中性稳定条件层结在其中的临界值划分为0.23～0.32;强稳定条件层结上的值划分为0.16～0.22。Zhang等[18]发现,临界理查森数在不稳定条件层结中的值是远远大于0.25的。在Esau等[19]的一次大涡模拟中发现,在夜间中性稳定边界层中理查森数的临界值与近中性稳定条件层结中的临界值不同。Richardson等[20]和Basu 等[21]发现，临界理查森数和大气稳定性参数存在线性关系。Zhang等[22]在野外进行观测实验，发现当稳定性的增强时，临界理查森数会减小，在不同条件下得出了不同的临界理查森数。在当前大多数情况下，通常会计算临界理查森数并将其用作每个模式的固定值。

基于以上研究背景，现采用WRF/YSU模式，设定不同临界理查森数，选择2016年8月江苏省一次晴好天气过程作为背景天气，探讨临界理查森数的变化对边界层特征和湍流通量交换的影响。

1 数据及模式设置

1.1 观测数据

选取江苏气象局24个基准气象站点的2 m气温，2 m相对湿度和10 m风速的观测数据对模拟进行了评估。24个站点所在区域下垫面类型分为城市站点(12个)和郊区站点(12个)。并运用南京信息工程大学观测塔的涡度相关系统所观测的感热和潜热通量的观测数据对模拟结果进行了验证。涡度分析相关检测系统主要由两个部分组成，分别为红外超声气旋风速仪(csat3,campbell,美国)和红外气体分析仪(li-7500,li-cor,美国)[23-25]。通过结合测量气体风速、声波中的温度以及测量气体浓度通量来进行计算可以得到气体感热的通量和潜热的通量。仪器架设在南京信息工程大学校园内(118.71°E,32.20°N)的一座铁塔上，安装高度为60 m，仪器采样频率为10 Hz。

1.2 模式及方案

WRF是中尺度预报模式和同化系统[26-27]，自2001年首次公开发布以来，已经成为世界上使用最广泛的天气尺度数值模型[28-29]。YSU边界层参数化方案在考虑大尺度湍流的逆梯度输送同时，并考虑夹卷层与自由大气的相互作用，物理机理较为完善，是目前WRF中最多常用的边界层参数化方案[30-33]。

在数值模式中，理查森数采用两层的位温和风速梯度进行计算，可以被用来判定大气边界层内的湍流状态。但是在本文的研究方案中，近地层方案采用Monin-Obukhov方案[34]，因此，湍流稳定度的判定采用莫宁奥布霍夫长度，没有计算理查森数。模式实际进行计算时,Monin-Obukhov近地层方案中,计算莫宁奥布霍夫长度来直接判定湍流状态稳定度,当湍流为不稳定湍流状态时,YSU方案中所选取临界理查森数模式为0.25。本文主要研究的临界理查森数直接参与边界层高度诊断计算公式的实际计算,对边界层高度诊断有明显影响,并进而影响湍流交换及边界层结构。

1.3 算例及参数设置

图1 d03区域的土地利用类型和站点分布Fig.1 Land-use types for the domain d03 and location of sites

本文研究区域为江苏省，选取2016年8月14—18日晴好天气过程作为研究时段。选择一次夏季的晴好天气过程为背景，在这种天气背景条件下，临界理查森数对大气边界层的影响最为直接。如图1所示，模拟区域采用三重嵌套，水平网格间距分别为27、9、3 km，中心经度和纬度分别为107°E和37°N[35-36]。最内层区域的土地利用类型如图1所示。气象场初始及侧边界条件使用1°×1°的NCEP再分析数据, spin-up时间为72 h，模式在垂直方向分成30层， 2 km以下有15层。参数方法包括YSU边界层方案[37]、Monin-Obukhov 近地层方案[34]、RRTMG长波辐射方案[38]、RRTMG 短波辐射方案[38]、Noah陆面方案[39]、Purdue Lin 微物理方案[40]，Single-layer城市冠层方案[41]。在不稳定层结下，设置了3个算例，分别为CTL算例(模式中，临界理查森数默认值设置为0.25)，High算例(设置为0.5)和Low算例(设置为0.12)。稳定层结条件下，Hong[42]认为，夜间临界里查森数(Ricr)和风速有关，因此设定其为函数形式:

Ricr=0.16(10-7R0)-0.18

(1)

式(1)中：R0为地表罗斯贝数。

(2)

式(2)中：U10为距离地面10 m处的风速；f0为科里奥利参数；z0为地表粗糙长度。

2 模拟评估

2.1 气象要素模拟评估

图2为江苏省城市和郊区站点气温、相对湿度和风速模拟值(CTL)和观测值(Obs)的对比以及统计参数的对比。模拟时间为2016年8月14—18日。

图2 2 m气温、2 m相对湿度和10 m风速模拟和观测值变化对比和参数统计及箱式统计图Fig.2 Diurnal variations of simulated and observed 2 m air temperature, 2 m relative humidity and 10 m wind speed with statistical parametersand box plots

由图2(a)、图2(b)、图2(g)和图2(h)可知，温度和相对湿度的观测值和模拟值的相关系数(R)都大于0.9，城市和郊区站点温度模拟结果类似，较好地体现出了观测值的特征[图2(a)和图2(b)]。平均偏差(MB)均在1.3 ℃左右，均方根误差(RMSE)约为1.5 ℃。城市和郊区站相对湿度的模拟平均值有一定差异，城市站的模拟结果优于郊区站[图2(c)、图2(d)]。城市、郊区站模拟与观测的平均偏差为-2.03%和-4.94%，均方根误差分别为4.11%和6.21%。分析其原因，本次模拟背景天气晴好，风速较小，近地层湿度的变化主要来源于局地蒸发，因此有可能是因为郊区植被的蒸发过程模拟存在一定的偏差。更加准确的下垫面数据，与不同种类植被生长有关的蒸发过程、庄稼地人为浇灌过程以及地表水文过程的准确刻画，对于局地蒸发的准确计算都有所帮助[43]。

在图2(e)及图2(f)中能看到，模式能够模拟出较为准确的风速随时间的变化，平均偏差和均方根误差均小于2 m/s。整体来看，无论城市站点还是郊区站点10 m风速的模拟值都偏大，城市的偏大更为明显。这很可能是因为模式采用最内层3 km的模拟分辨率，所选用的单层城市冠层方案无法详细考虑真实城市的高度和密集程度，模式低估了建筑物下垫面的粗糙度，导致模拟出的风速比观测值大。这种效应在城市下垫面更为明显。进一步的研究中，引入更加准确的建筑物资料，对于风速的模拟，可能有较大的帮助。

由箱式图2(g)～图2(j)可以看出，观测的2 m温度和相对湿度与模拟的中值具有相同的变化趋势，模拟结果更好，但最大值和最小值存在一些差异，风速同样呈现出了模拟值大于观测值的状况。模拟温度与观测值的最大误差值出现在中午，最小值出现在凌晨，可以看出模拟值在最高值和最低值处都偏小。对温度低估的原因是单层城市冠层方案低估了城市中实际建筑物高度及密集程度，模拟出了较大的风速[图2(f)、图2(l)]，局地热量被平流项带走，因此2 m气温的模拟相比于观测值要略低。

2.2 感热通量和潜热通量模拟评估

运用南京信息工程大学观测塔的感热和潜热通量的观测值进行模拟评估如图3所示，感热通量的相关系数为0.86，模拟值较大，但与观测值的变化曲线相似，平均偏差为-9.04 W/m2。潜热通量的相关系数为 0.85，模拟值较小。观测塔所处的下垫面分布是非均匀的，包括建筑物，植被和人工湖，然而在模式中，该点下垫面是均一的，忽略了城市水文过程的影响，如绿化灌溉和人为潜热，使得模拟的感热通量呈现增加的趋势[44]。潜热通量的模拟值偏小，平均偏差为38.72 W/m2。绿地灌溉后，随着土壤水分的增加和白天温度的升高，蒸散增强,且观测塔周围的人工湖水份蒸发，这些会使得观测的潜热通量大于模拟的潜热通量[45-46]。但是总体来说，模式较好地模拟出来地表通量的峰值和日变化过程。

图3 南京信息工程大学观测塔的观测和模拟的感热通量和潜热通量Fig.3 Observed and simulated sensible heat flux and latent heat flux of observation tower of Nanjing University of Information Science & Technology

3 结果分析

在WRF/YSU模式中，临界理查森数是重要的参数。用来计算边界层高度。

当边界层高度发生变化后，湍流扩散系数随之产生变化，大气边界层中物质能量的扩散过程也发生变化。

3.1 临界理查森数对边界层高度的影响

YSU方案中，边界层高度诊断公式[47]为

(3)

式(3)中：h为边界层高度；Ricr为临界理查森数。当层结处于不稳定状态时，WRF/YSU方案中临界理查森数默认取值为0.25(CTL算例)；θva为模式最底层虚位温；U(h)为高度h处的风速；θv(h)为高度h处的位温；θs为近地表位温。从式(3)中可以看出，临界理查森数的大小直接影响边界层高度的诊断。

以南京站为例，图4给出了2016年8月14—18日3个算例的边界层高度随时间的变化情况。3个算例模拟的变化趋势一致，白天，地面接收的太阳辐射增强, 导致地面增加温度更快，使地表和大气的温度相差得更大，湍流增强, 促进边界层发展，边界层高度最高值出现在14：00—16：00，临界理查森数对边界层高度的影响主要体现在对边界层高度日变化最大值的影响。白天湍流运动旺盛，模拟出的边界层高度较大。从式(1)中可以看出，白天，边界层高度和临界里查森数成正比，因此大的临界理查森数模拟出来的边界层高度更大。

图4 南京站2016年8月14—18日不同算例的边界层高度日变化曲线Fig.4 Diurnal variation curves of boundary layer height of different cases in Nanjing station from August 14, 2016 to August 18, 2016

在白天不稳定的条件下，大多数时刻边界层高度随临界理查森数增大而增大，有些时刻出现了，临界理查森数大，而边界层高度小的情况，如8月15日10：00；这跟湍流运动有关。对流条件下，湍流大涡的运动会导致边界层内的标量及矢量在时间和空间上出现非周期性和不连续性的，导致边界层高度在空间上分布是不均匀的。因此会出现有些时刻边界层高度随临界理查森数增大而减小的情况。

而在夜间，由于YSU方案中，临界理查森数采用风速的函数进行计算(式1)，本文的算例中对于临界理查森数的设定并无作用，因此在夜间，临界理查森数的变化对边界层高度并无明显的影响。

由图4可知，8月15日的16：00High算例边界层高度出现最大值，约为1 700 m，CTL算例模拟的边界层高度约为1 343 m，Low算例模拟边界层高度最低，约为1 298 m，与Low算例的边界层高度差值最大，达到402 m。18日的16：00High算例边界层高度与Low算例的边界层高度差值也较大，为369 m。本文在后面部分的分析中，以8月15日为例，研究临界理查森数对边界层的影响。

图5给出了2016年8月15日12：00和16：00各算例模式d03区域[图1(a)]边界层高差空间分布图。图5(a)和图5(b)分别显示的是12：00和16：00High算例与CTL算例的边界层高度差值，可以看出，当理查森数的临界值增加时，边界层高度也会增加。图5(c)和图5(d)分别显示的是12：00和16：00Low算例与CTL算例的边界层高度差值，临界理查森数的减少时，边界层高度也减小。这种影响在下午最为明显，边界层高度相差最多达到500 m；而在12：00左右各算例边界层高度差异较小，大部分地区边界层高度相差小于100 m。在大多数地区,在理查森数的临界值持续增加时,边界层会受到很大影响,高度增加,但也有一些地区是随之减少的。这跟一个关于湍流运动的一个空间非均匀性的问题有关,会直接影响导致边界层的高度分布不是均匀的,出现有些情况，在湍流运动区域边界层高度减小,可以说是随临界理查森数不断持续增大而逐渐略微减小的一种异常情况。

3.2 临界理查森数对地气交换的影响

YSU方案的理查森数的临界值主首先影响了边界层高度，当边界层高度发生改变时，边界层内湍流扩散系数也产生了改变，使得边界层内热量、水汽等物质的垂直输送发生了变化，导致地气交换通量的分布存在显著的变化。式(4)给出了YSU中湍流扩散系数的计算方案。从式(4)中可以看出，当边界层高度增大以后，湍流扩散系数有一定程度的增大[31]。

(4)

式(4)中：p为剖面形状指数，取值为2；k为卡曼常数，取0.4；ws为混合层速度尺度；z为距离地面高度；h为边界层高度。

图5 High算例、Low算例与CTL算例边界层高度诊断差值空间分布图Fig.5 Spatial distribution of the diagnostic difference of boundary layer height between High, Low and CTL cases

图6 High算例和CTL算例的水汽混合比(色阶)差值和湍流交换系数(等值线)差值的垂直剖面图Fig.6 Profile of differences of water vapor (shaded) and differences of turbulent exchange coefficient (contour) for the High and CTL cases

图6给出了2016年8月15日16：00High算例与CTL算例水汽和湍流交换系数差值剖面图(位置如图5(b)中AB剖面)。可以看出High算例和CTL算例的湍流交换系数差值的最大值为80 m2/s，水汽混合比差值最大为2.2 g/kg。当临界理查森数增加时，边界层的高度随之增加，湍流交换系数增大，湍流的交换增强，水汽混合比增加。垂直运动增强，水汽的垂直混合增加。

模式中地面感热通量和潜热通量的计算在陆面过程模块进行。图7为2016年8月15日16：00High算例与CTL算例地面热通量差值(感热和潜热)空间分布情况。临界理查森数增加，地表感热通量和潜热通量在空间上总体增加。地表潜热通量与感热通量的变化趋势相同，但是潜热通量变化幅度比地面感热通量变化幅度更大。临界理查森数会对水汽和地表热通量产生影响。临界理查森数增大，湍流扩散系数也增大，水汽混合比增加，并且地面热通量总体增大，又会进一步使边界层高度增大。

图7 8月15日16：00High算例和CTL算例的地面热通量差值空间分布图Fig.7 Spatial distribution of surface heat flux difference of High and CTL cases at 16:00 on August 15

3.3 临界理查森数对边界层通量影响分析

目前，不同尺度的数值天气预报模型和空气质量预报模型无法辨别小于模型网格尺度的边界层湍涡，通常将湍流通量参数化，可用于计算大气边界层中的热量、水汽、动量和其他物质的输送。YSU边界层参数化方案使用K理论来进行湍流参数化，并假设湍流输送方式相同。湍流扩散项的方程为

(5)

(6)

夹卷通量计算时，YSU方案并不直接求解湍流通量项，而是直接使用了参数化方案，设定夹卷通量是地面热通量的0.15倍。夹卷过程控制了边界层内物质和自由大气的交换，当地表热量越强，夹卷通量越大，自由大气内更多的大气被卷入大气边界层中，边界层高度又进一步增长。

从前面的分析可以看出，16：00边界层高度差异最大，此时边界层高度受到的影响最大。为了讨论临界理查森数对边界层湍流的影响，本节将以16：00为例，图8给出了2016年8月15日16：00边界层内各算例水汽通量、热通量和动量通量的垂直分布情况。正值代表向上的通量，负值代表向下的通量。High算例和CTL算例边界层下部向上水汽通量和边界层上部向下水汽通量的值相当，而Low算例模拟出的向下水汽通量大于向上的水汽通量。随着临界理查森数的增加，边界层上部向下的水汽通量逐渐减小，边界层下部向上的水汽通量变化不明显。对于热通量，High算例和Low算例的热通量变化相差不大，CTL算例向上和向下的热通量是最小的。CTL算例边界层下部向上的动量通量和边界层上部向下的动量通量的值相当，High算例和Low算例的动量通量主要是在边界层下部向上传输。这种变化可由局地通量项，非局地通量项和夹卷通量项的变化共同解释。

3.3.1 临界理查森数对局地通量项的影响

由式(5)可见，临界理查森数对局地通量项的影响，主要表现在湍流扩散系数与局地梯度的变化。当临界理查森数变化，边界层高度发生变化，湍流扩散系数必然随之变化，进而影响局地通量项的计算。边界层高度越高，与之对应的边界层内的湍流扩散系数就越大。边界层内最大湍涡所影响的高度即为边界层高度，边界层高度增加意味着边界层内最大湍涡的增加，湍流输送能力增强，能将物质输送到更高的高度。

图8 8月15日16：00不同算例边界层水汽通量、热通量和动量通量廓线Fig.8 Profiles of water vapor and heat flux and momentum flux in the boundary layer with different cases at 16:00 on August 15

图9 8月15日16：00不同算例边界层湍流扩散系数、局地热通量项、局地水汽通量项和局地动量通量项的垂直廓线Fig.9 Vertical distribution of turbulent diffusion coefficient, local heat flux term, local water vapor flux term and local momentum flux term for the boundary layer in different cases at 16:00 on August 15

图9给出了2016年8月15日16：00湍流扩散系数、局地热通量项、局地水汽通量项和局地动量通量项垂直变化情况。局地通量项为湍流扩散系数与局地梯度的乘积，局地通量的变化受二者共同影响。由图9(a)可以看出，临界理查森数越大，湍流扩散系数越大。但是局地热量通量[图9(b)]及水汽通量[图9(c)]随高度向上的变化趋势受到局地梯度的影响，垂直分布与湍流扩散系数的分布不一致。随着临界理查森数的增加，局地通量项的量值呈现增加的趋势。热的局地通量差异主要体现在边界层中上部向下的通量，这主要是由于热量以及水汽的局地梯度在这个高度上比较显著。High算例模拟出的向下热通量最大，CTL算例次之，Low算例模拟出的向下热通量最小。临界理查森数增大，边界层内模拟出的热通量方向向下的区域增大。水汽的局地通量在整个边界层内都为正，代表局地项向上输送水汽，且临界理查森数越大，向上输送的水汽通量可以达到更高的高度。与局地热通量相比，水汽的局地通量向上输送的高度略低，说明局地通量对热量的输送作用要比对水汽的输送作用强。High算例和Low算例的局地动量通量变化趋势相近，在边界层下部向上输送的通量比上部输送的通量值更大，CTL算例中边界层下部向上的动量通量和边界层上部向下的动量通量的值相差不大。

3.3.2 临界理查森数对非局地通量项的影响

由于边界层总体梯度的存在，边界层内最大尺度的湍涡(即湍涡影响范围为整个边界层)会将地面动量、热量、水汽及其他物质向上输送，即使在某些高度局地梯度为正。例如，热力在对流边界层的中上部占主导地位，热量可以从冷区转移到热区，违反温度梯度的限制。或者仅有很小的温度梯度，仍然有显著的热流存在，这两种现象为逆梯度输送现象。

YSU边界层参数化方案考虑了逆梯度输送的作用，在局地混合所描述的梯度扩散项中加入了考虑非局地扩散的逆梯度修正项。图10给出了热量非局地项、水汽非局地项和动量非局地项通量廓线图。从图10中可以看出，热量、水汽和动量的非局地通量，在整个边界层内非局地通量都始终为正值。与CTL算例相比，High算例的非局地通量变化更为明显。较低的临界理查森数对非局地通量的输送影响不大。以热量通量为例，临界理查森数增加，边界层高度增高，地表热通量减小，根据式(6)可得出，逆梯度项γc减小，则热量的非局地通量项增大，非局地热通量向上输送的高度增高。非局地通量项对热量的输送作用与局地项相当，而非局地通量对水汽和动量的输送作用却远远小于局地通量项。

3.3.3 临界理查森数对边界层夹卷过程中的影响

夹卷过程不仅直接影响大气边界层的发展，同时也会对大气中的其他物理过程产生一定的影响。图11给出了热量、水汽和动量夹卷通量项在边界层内分布情况，可知夹卷通量项对热量、水汽和动量在边界层内垂直分布作用相同，在整个边界层内夹卷通量项方向均向下，夹卷通量项主要影响体现在边界层中上部区域，对边界层下部影响较小。由4.2节分析可知，临界理查森数增加，湍流扩散增强，水汽混合增加，夹卷通量项呈现出减小的趋势。对于热量通量而言，夹卷通量项对总热通量的贡献与局地项、非局地项相当。水汽和动量的夹卷通量与热量夹卷通量的变化趋势相同，也是随着临界理查森数增加而减小。与热量的夹卷通量项不同，局地水汽和动量通量大小与夹卷通量大小相差不大，边界层中水汽和动量的局地输送和夹卷作用同样重要。高临界理查森数和较低的夹卷通量对应，使边界层和自由大气之间的交换减弱，更多的水汽保留在边界层上部，导致边界层上部湿度更大。临界理查森数越低，反而使夹卷通量越高。

图10 8月15日16：00不同算例的非局地热通量项、非局地水汽通量项和非局地动量通量项垂直廓线Fig.10 Vertical profiles of non-local heat flux term，non-local water vapor flux term and non-local momentum flux term of different cases at 16:00 on August 15

图11 8月15日16：00不同算例边界层内夹卷热通量、夹卷水汽通量和夹卷动量通量垂直分布Fig.11 Vertical profiles of entrainment heat flux,entrainment water vapor flux and entrainment momentum flux for the boundary layer of different cases at 16:00 on August 15