中国气象学会主办。
文章信息
- 全利红, 袁春红, 王丙兰. 2016.
- QUAN Lihong, YUAN Chunhong, WANG Binglan. 2016.
- 应用印痕函数研究测风塔资料代表性
- Application of footprint function in characterizing the representation of wind tower data
- 气象学报, 74(3): 442-449
- Acta Meteorologica Sinica, 74(3): 442-449.
- http://dx.doi.org/10.11676/qxxb2016.030
-
文章历史
- 2015-10-22 收稿
- 2016-03-16 改回
风况观测数据是桥梁、建筑等工程抗风设计、风能资源开发、电力、铁路调度运行等必不可少的基础依据。换言之,风况实地观测数据的效果影响着许多工程建设设计和运行维护的全过程。因此,提高测风数据的准确性和可靠性对工程建设和运行具有非常重要的意义。为了得到准确的观测数据,需要对风况观测进行科学、合理的设计,以及对观测结果的代表性进行分析。
测风塔是风况观测的重要手段。尽管目前已有相关规范对测风塔的布设和安装做出了规定和要求,但是,针对复杂地形的相关规定过于笼统。在复杂地形上,气流随地形的变化而变化,具有很强的局地性,单个测风塔的观测结果往往不能满足要求。因此,在设计测风方案时,测风塔的布设需要结合地形和气流运动特征来确定,要使得测风塔的观测结果既能够反映整个关心区域的特征,又不能过多布设测风塔,增加建设成本。解决这个问题,首先需要了解测风塔观测结果的代表性,即单个测风塔的观测结果能够反映多大范围的气流特征。
目前,复杂地形上测风塔位置的确定主要还是依靠经验分析,根据地形地貌、下垫面粗糙度变化、障碍物的分布等大致估计。包小庆等(2008)通过对不同地形气流特征的分析,提出了风电场测风塔选址原则。在复杂地形上确定测风塔位置还可以通过数值模式计算发电量来反推(杨冬雷,2011)。然而,实际上测风塔观测结果综合反映了多种因素的特征,包括下垫面粗糙度、大气稳定度、风向风速、天气系统等。为了更合理地确定测风塔观测结果能够代表的区域范围,需要通过适当的方法综合考虑这些因素。
印痕函数方法(或源权重方法)是一个综合考虑风速、大气稳定度、下垫面粗糙度等因素对观测量影响的方法,可以用来分析单点观测结果的代表性(蔡旭辉,2008)。印痕函数f(r,r′)是一个将单点观测结果η(r)与上风向的“源”项或强迫项Qη(r+r′)联系起来的传递函数η(r)=∫RQη(r+r′)f(r,r′)dr′,其中,r是观测位置,r′是观测点与“源”区的距离,它反映了地-气相互作用对观测结果的影响(Schmid,2002)。Schmid(2002)曾给出一个印痕函数的概念模型。在这个模型中,距观测点很近的范围源权重比较小,随着距观测点距离的增大,源权重增大,在某一点达到最大值,然后随着距观测点距离的增大而减小。早期在研究湍流通量低估问题时,为了研究观测资料的代表性引出了印痕函数。目前,印痕函数方法已用于湍流通量和标量浓度的代表性研究中(Vesala et al,2008)。由于印痕函数方法在一定程度上反映了观测信息与地面状况的关系,在对湍流大气中进行的各类观测的分析和解释都是十分有用的(蔡旭辉,2008)。
目前几乎所有的印痕函数模型(或公式)都是基于湍流脉动量建立的(Leclerc et al,1997;Hsieh et al,2000,2009;Schmid,2002; Lee,2003;Kljun et al,2004;Vesala et al,2008)。然而,测风塔上主要是常规梯度观测资料,湍流观测资料很少。根据Horst(1999),梯度观测资料与涡相关法得到的印痕函数的上风向距离相近。因此,本研究将尝试结合常规梯度观测资料和印痕函数分析测风塔观测资料的代表性。
2 测风塔资料
所用测风塔资料来自中国风能资源观测网。中国风能资源观测网共有400多座测风塔,自2007年开始建设。此资料主要用于第4次全国风能资源详查。在《全国风能资源详查和评价报告》(2014)中对各测风塔的布设情况进行了详细介绍。在此选取其中大风较多的44座测风塔2010年的资料进行分析。这44座测风塔分布在新疆、内蒙古、西藏、广东、湖南等省区,下垫面类型包含了平坦地形以及山脉丘陵等复杂地形(图 1)。其中,1—6号塔在海拔较高的草原上,17—23号塔在沿海地区或海岛上,39、40号塔在天山地区三十里风区附近、41号塔在天山地区百里风区附近,9号塔周围有农田和城镇,15号塔在城镇附近,7、30号塔在山坡上,24、43号塔在湖边,27、28号塔在西南山区较高的山脊上(海拔高度高于3000 m)。
44座测风塔有70和100 m两种高度。70 m高测风塔在10、30、50、70 m各层装有测风仪器;100 m高测风塔在10、30、50、70、100 m各层装有测风仪器;各塔10和70 m高度层有温、湿度传感器,8.5 m高度层装有气压传感器。部分测风塔上还装有超声风速仪。其中,测风塔上风速传感器型号为EL-15-1A,风向传感器型号为EL-15-2D,温、湿度传感器型号为HMP45D,强风传感器型号为XFY3型,仪器性能指标见表 1。观测资料每10 min保存一次,包括各层平均风速、风向、风速方差、最大风速、温度、气压等。
风速 | 0—60 m/s(强风仪:0—90 m/s) | 0.1 m/s | ±(0.5+0.03V) m/s | 3 s1 min | 1 s-1 |
风向 | 0°—360° | 3° | ±5° | 2 和10 min | |
温度 | -50—+50℃ | 0.1℃ | ±0.2℃ | 1 min | 6 min-1 |
湿度 | 0—100% | 1% | ±4%(≤80%),±8%(>80%) | 1 min | 6 min-1 |
印痕函数方法中一个非常重要的变量是空气动力学粗糙度。中性层结条件下,粗糙度可由下式得到
式中,u是z高度的风速,u*是摩擦速度,k是冯卡门常数,z0是粗糙度。通过拟合中性条件下风廓线得到z0。拟合风廓线与观测风廓线的相关系数达到99.9%(图 2)以上的数据才用来进行粗糙度计算,取其平均值作为输入印痕函数中的粗糙度值。根据Smeets等(2008),层结稳定性由10和70 m两层的风速u和温度T确定
其中,,g是重力加速度, T是两层的平均温度。这里取|Lp|<0.2为中性层结。按此方法挑选出来的资料大部分风速较大,70 m高度平均风速大于6 m/s的测风塔占90%以上。
3.2 印痕函数印痕函数的模型(或公式)很多,有欧拉型、拉格朗日型以及基于大涡模式的模型(Schmid,2002;Vesala et al,2008)。相关的参数也可能存在差异,除了风速、稳定度和下垫面粗糙度参数外,很多印痕函数模型中还包含了其他参数,如垂直速度方差、摩擦速度等(Leclerc et al,1997;Hsieh et al,2000,2009; Lee,2003;Kljun et al,2004)。此外,不同的印痕函数模型使用条件也不同。大涡模式和一些拉格朗日模型需要较多的计算时间,在此不予考虑。由于测风塔梯度观测资料中没有包含湍流量和垂直速度等的观测信息,不能直接得到速度方差和摩擦速度等参数值,只能利用梯度观测结果计算通量,因此,包含垂直速度或方差的模型也不能使用。最后选定了两种简单易实现的印痕函数模型,分别为Schuepp等(1990)和Hsieh等(2000)的印痕函数模型。
Schuepp印痕函数模型是一个欧拉型的解析模型,首先由Schuepp等(1990)通过近似求解平流-扩散方程得到,用来分析中性层结条件下标量通量的源权重。中性条件下印痕函数微分模型可以表示为
式中,f(x)是印痕函数,d是零平面位移,x是观测点距上风向的距离,U是地表至观测高度z的平均风速,。通过对风廓线的拟合,得到z0。尽管Schuepp模型假设下垫面均匀,但是,根据Schmid(2002),Schuepp模型与Leclerc等(1990)的考虑了垂直不均匀性的拉格朗日模型计算结果非常一致。Hsieh印痕函数模型是Hsieh等(2000)利用实测资料对一个拉格朗日随机扩散模型结果进行回归分析得到的,可用于各种层结。函数的微分形式为
式中,L为奥布霍夫长度,zu=z(ln(z/z0)-1+z0/z),D和P是两个随稳定度变化的参数,中性条件下,D=0.97,P=1;不稳定条件下,D=0.28,P=0.59;稳定条件下,D=2.44,P=1.33。 4 结果分析由于大部分测风塔都安装在复杂地形上,各个方向下垫面的状况不同,测风塔的代表性范围也不同,需要分不同来流方向来分析其代表范围。为了确保中性层结的样本数,选取70 m高度主导风方向和风向频率较高的其他方向的资料进行计算分析。
4.1 粗糙度地形坡度和风速对粗糙度影响很大。粗糙度最大值出现在7号塔西-西北西方向。7号塔位于山坡上,与山顶海拔高度差约400 m,水平距离约23 km。40号塔也在山坡上,与山顶海拔高度差约800 m,水平距离约38 km,粗糙度比7号塔西-西北西方向的粗糙度小了一个量级。30号塔同样也位于山坡上,与山顶海拔高度差约800 m,但是,水平距离约150 km,因此,粗糙度比40号塔西北-北西北方向粗糙度又小了一个量级。 由此可见,对于山坡上的塔和沿山坡的风向,坡度越大,则粗糙度越大。虽然有不少测风塔在山地丘陵地区,但是,由于山坡较缓或者在山顶,粗糙度并不大。例如27和28号测风塔,虽然处在西南山区,但是,均在山脊上,海拔高度(分别高于3200和3600 m)较附近地区高,山坡非常陡,因此,与山脊垂直的方向粗糙度值较小。39和40号塔平均风速都很大,分别为15.2和17.1 m/s,但是,39号塔位于狭谷中,沿狭谷方向地形比较平坦,而且,由于风速很大,沙石也很容易被吹走,所以下垫面很光滑,粗糙度很小。同样位于狭谷地区的43号塔,平均风速只有6.9 m/s,所以粗糙度比39号塔大得多(0.0067 m)。9号测风塔西南-西和东南方向主要是农田和村庄,但是,西南-西方向平均风速为5.4 m/s,而东南方向平均风速为13.7 m/s,相应地,西南-西方向的粗糙度(0.5142 m)比东南方向的粗糙度(0.0014 m)大得多。
4.2 印痕函数在研究通量或标量浓度的源权重时,通常取印痕函数积分为90%的范围作为代表性范围,根据印痕函数的概念模型(Schmid,2002),认为对单点观测的通量或标量浓度有贡献的源90%都来自于这个范围。对平均风速而言,由于湍流扩散的影响,如果也选取印痕函数积分90%的范围,则意味着在这个距离的风速到达观测点时90%的风速特征已发生改变。对于较大的风速,如系统性的大风、台风等,即使风速改变了90%可能对观测点的影响还不能忽略时,需要考虑印痕函数积分90%以上的距离作为测风塔的代表性范围。对于较小的风速,90%的改变可能对观测点的影响已经很小,可以忽略,这种情况需要考虑印痕函数积分值较小时对应的距离作为测风塔的代表性范围。在此选择两种阈值进行分析,其中一种是对测风塔观测结果影响最大的位置,即印痕函数最大值对应的距离,作为其代表性范围。这个距离通常在印痕函数积分为30%的范围内,有些甚至只是印痕函数积分的10%左右。图 3是印痕函数分布的一个例子(第16号测风塔),其中,印痕函数最大值对应积分值为13.5%。另一种是印痕函数积分90%的范围。
表 2给出了两种方法计算的印痕函数最大值对应距离及印痕函数积分90%对应的距离。从表中可见,两种方法得到的距离比较接近。印痕函数最大值对应距离大部分在200—3000 m;印痕函数积分90%对应的距离大部分在4—60 km。
测风塔编号 | 风向 | z0(m) | Schuepp方法(m) | Hsieh方法(m) | 测风塔编号 | 风向 | z0(m) | Schuepp方法(m) | Hsieh方法(m) | ||||
xmax | x90% | xmax | x90% | xmax | x90% | xmax | x90% | ||||||
1 | NW-N | 0.0037 | 1936 | 36742 | 1880 | 35640 | 23 | NE-E | 0.00092 | 2240 | 42520 | 2170 | 41240 |
2 | WNW-NW | 0.0016 | 2119 | 40222 | 2060 | 39010 | 24 | N-NNE | 0.0047 | 1883 | 35749 | 1830 | 34670 |
3 | NW-N | 0.0056 | 1845 | 35022 | 1790 | 33970 | 25 | N-NE | 0.0038 | 1930 | 36632 | 1870 | 35530 |
4 | SW-W | 0.0011 | 2201 | 41778 | 2130 | 40520 | 26 | S-SSW | 0.002 | 2070 | 39296 | 2010 | 38120 |
4 | SE | 0.00018 | 2597 | 49294 | 2520 | 47810 | 27 | SW-WSW | 0.000065 | 2820 | 53523 | 2740 | 51920 |
5 | SW-W | 0.0011 | 2201 | 41778 | 2130 | 40520 | 28 | SSW-WSW | 0.000043 | 2910 | 55239 | 2820 | 53580 |
5 | WNW-NNW | 0.0023 | 2040 | 38716 | 1980 | 37550 | 29 | WSW-WNW | 0.0032 | 1967 | 37345 | 1910 | 36220 |
6 | WNW-NNW | 0.00089 | 2247 | 42657 | 2180 | 41380 | 30 | NNE-ENE | 1.2395 | 680 | 12899 | 650 | 12290 |
7 | W-WNW | 118.37 | 40 | 670 | 31 | E-ESE | 0.0127 | 1666 | 31626 | 1620 | 30670 | ||
8 | WSW-WNW | 0.8056 | 769 | 14603 | 740 | 14000 | 31 | W-WNW | 0.0127 | 1666 | 31626 | 1620 | 30670 |
9 | SW-W | 0.5142 | 864 | 16402 | 830 | 15790 | 32 | N-NNE | 0.2236 | 1042 | 19785 | 1010 | 19130 |
9 | SE | 0.0014 | 2148 | 40777 | 2080 | 39550 | 33 | NW-N | 0.001 | 2222 | 42174 | 2160 | 40910 |
10 | S-SSW | 0.00012 | 2685 | 50977 | 2600 | 49450 | 33 | WSW-W | 0.0113 | 1692 | 32110 | 1640 | 31140 |
11 | W-NW | 0.0073 | 1787 | 33922 | 1730 | 32900 | 34 | NW-NNW | 0.5769 | 840 | 15939 | 810 | 15330 |
12 | WSW-WNW | 0.0053 | 1857 | 35251 | 1800 | 34190 | 34 | SSE-S | 0.0017 | 2106 | 39971 | 2040 | 38770 |
13 | WSW-WNW | 0.0026 | 2013 | 38207 | 1950 | 37060 | 35 | WNW-NW | 0.00038 | 2433 | 46191 | 2360 | 44810 |
14 | WNW-NNW | 0.0013 | 2164 | 41084 | 2100 | 39850 | 36 | E-ESE | 0.0107 | 1703 | 32336 | 1650 | 31360 |
14 | SSE | 0.0115 | 1688 | 32037 | 1640 | 31070 | 36 | W-WNW | 0.005 | 1870 | 35493 | 1810 | 34430 |
15 | WNW-NNW | 1.4108 | 653 | 12393 | 620 | 11780 | 37 | WSW-W | 0.0022 | 2049 | 38900 | 1990 | 37730 |
16 | S-SW | 0.00018 | 2597 | 49294 | 2520 | 47810 | 38 | WSW | 0.0208 | 1558 | 29581 | 1510 | 28680 |
16 | NW-NNW | 0.00029 | 2492 | 47313 | 2420 | 45890 | 38 | WNW-NW | 0.00093 | 2238 | 42475 | 2170 | 41200 |
17 | NNE-ENE | 0.002 | 2070 | 39296 | 2010 | 38120 | 39 | WNW-NW | 1.22E-05 | 3186 | 60470 | 3090 | 58660 |
18 | NNE-NE | 0.000059 | 2841 | 53925 | 2760 | 52310 | 40 | NW-NNW | 12.3787 | 242 | 4587 | 190 | 3660 |
19 | NNE-NE | 0.000011 | 3208 | 60900 | 3110 | 59070 | 41 | NNW-N | 3.44E-05 | 2959 | 56165 | 2870 | 54480 |
19 | SW | 0.0227 | 1539 | 29218 | 1490 | 28330 | 41 | WSW-W | 0.0038 | 1930 | 36632 | 1870 | 35530 |
20 | N-NE | 0.002 | 2070 | 39296 | 2010 | 38120 | 42 | E-ESE | 0.2887 | 987 | 18742 | 950 | 18110 |
20 | SSW | 0.0011 | 2201 | 41778 | 2130 | 40520 | 42 | NW-NNW | 0.000098 | 2730 | 51818 | 2650 | 50260 |
21 | N-NNE | 0.0707 | 1292 | 24519 | 1250 | 23760 | 43 | SE-SSE | 0.0067 | 1806 | 34278 | 1750 | 33250 |
22 | NNE-ENE | 0.0054 | 1853 | 35173 | 1800 | 34120 | 44 | WNW-NW | 0.0028 | 1997 | 37899 | 1940 | 36760 |
印痕函数不同取值对应的距离与下垫面的变化有关。例如,16号测风塔在西北方向约2500 m是农田和村庄的分界线,两种方法计算得到的印痕函数最大值对应的距离分别为2492和2420 m,与测风塔距下垫面分界线距离接近;20号测风塔北偏东方向约2000 m是海陆分界线,两种方法计算得到的印痕函数最大值对应的距离分别为2070和2010 m,与测风塔距下垫面分界线距离接近。又如10号测风塔,在南方约50 km是戈壁和农田的分界线,两种方法计算得到的印痕函数积分为90%时对应的距离分别是50.977和49.45 km,与测风塔距分界线的距离接近;43号测风塔东南偏南方向约33 km是湖陆交接处,两种方法计算得到的印痕函数积分为90%时对应的距离分别是34.278和33.254 km,与测风塔距交接处的距离接近。
根据Wiernga(1993),均匀下垫面粗糙度与风浪区长度可以表示为
从Hsieh方法计算的印痕函数不同积分值对应距离与粗糙度的关系(图 4)可见,随着粗糙度的增大,印痕函数不同积分值对应的上风向距离减小,粗糙度小于10 m时,粗糙度与印痕函数对应的上风向距离近似为对数-线性关系,这种变化趋势与式(5)一致,但是,测风塔得到的印痕函数积分值对应
的上风向距离随粗糙度的变化速率与式(5)不一致,而且粗糙度大于10 m时,印痕函数对应的上风向距离随粗糙度的变化开始变缓,而式(5)印痕函数对应的上风 向距离随粗糙度的变化在粗糙度大于100 m时才变缓。因此,测风塔的印痕函数积分值对应的距离与粗糙度的关系不能用式(5)来描述。图 5给出了粗糙度小于10 m时与印痕函数不同积分值对应的距离的拟合结果。粗糙度与印痕函数最大值对应距离的拟合结果分别为:xmax=-217.0lnz0+723.8(Scheupp方法),xmax=-211.6lnz0+693.8(Hsieh方法)。粗糙度与印痕函数积分90%对应的距离的拟合结果分别为:x90%=(-4.1lnz0+13.7)×1000(Scheupp方法),x90%=(-4.0lnz0+13.2)×1000(Hsieh方法)。
比较两种印痕函数计算方法得到的结果(表 1和图 5),随着粗糙度的增大,两者的差异减小。Scheupp方法是基于均匀下垫面假设得到的,在复杂山地则不适用,在粗糙度大于100 m时Scheupp方法不能算出xmax和x90%。
5 结 论本研究分析了44个测风塔2010年的梯度观测资料,根据10和70 m的风速和温度确定中性层结,选取中性层结条件下满足对数风廓线分布的数据,引入印痕函数方法对主导风方向和风向频率较大方向的测风塔资料的代表性范围进行了分析。
根据各种印痕函数方法的应用条件和测风塔的资料,选择两种方法:Scheupp方法和Hsieh方法,进行印痕函数的分析。粗糙度小于10 m范围,两种方法得到的代表性范围距离比较接近。但是,Scheupp方法适用于比较均匀的下垫面,在复杂山地则建议采用Hsieh方法。同时,Hsieh方法还可用于不同层结。
粗糙度小于10 m时,印痕函数不同积分值对应的上风向距离随粗糙度增大成对数-线性关系减小,与Wiernga(1993)的结果相比,测风塔得到的印痕函数积分值对应的上风向距离随粗糙度增大而减小得快。粗糙度大于10 m时印痕函数积分值对应的上风向距离随粗糙度增大减幅变小。Wieginga公式中系数2是由摩擦速度u*与垂直速度方差σw的比值u*/kσw估算得到的,在复杂地形,u*/kσw的值可能不为2,内边界层厚度可能也不是10z,因此,在复杂地形测风塔代表性范围不能简单由Wiernga(1993)公式(式5)计算得到。
借鉴通量和浓度研究的结果,测风塔的代表性范围也可取为Hsieh印痕函数积分值的90%对应的距离,但是,如果在这个范围内下垫面有显著变化,则测风塔的代表性范围应为测风塔所在下垫面变化之前的范围。由于印痕函数方法很多,这里只是选取两种较简单的方法进行了比较分析,进一步的研究需要借助于示踪实验、实验室实验或者数值模拟来检验印痕函数和积分阈值确定的合理性(Leclerc et al,1997;Foken et al,2004;Göckede et al,2005;Masseroni et al,2014)。
包小庆, 张国栋. 2008. 风电场测风塔选址方法. 资源节约与环保, 24(6): 55-56. Bao X Q, Zhang G D. 2008. The method to determine wind tower sites for wind plant. Resour Econom Environ Prot, 24(6): 55-56(in Chinese) |
蔡旭辉. 2008. 湍流微气象观测的印痕分析方法及其应用拓展. 大气科学, 32(1): 123-132. Cai X H. 2008. Footprint analysis in micrometeorology and its extended applications. Chinese J Atmos Sci, 32(1): 123-132(in Chinese) |
杨冬雷. 2011. 复杂地形风场测风塔位置的选择. 风能, (1): 69-71. Yang D L. 2011. Selecting the location of wind tower for wind field on complex terrain. Wind Energy, (1): 69-71(in Chinese) |
Foken T, Leclerc M Y. 2004. Methods and limitations in validation of footprint models. Agri Forest Meteor, 127(3-4): 223-234 |
Göckede M, Markkanen T, Mauder M, et al. 2005. Validation of footprint models using natural tracer measurements from a field experiment. Agri Forest Meteor, 135(1-4): 314-325 |
Horst T W. 1999. The footprint for estimation of atmosphere-surface exchange fluxes by profile techniques. Bound Layer Meteor, 90(2): 171-188 |
Hsieh C I, Katul G, Chi T. 2000. An approximate analytical model for footprint estimation of scalar fluxes in thermally stratified atmospheric flows. Adv Water Resour, 23(7): 765-772 |
Hsieh C I, Katul G. 2009. The Lagrangian stochastic model for estimating footprint and water vapor fluxes over inhomogeneous surfaces. Int J Biometeor, 53(1): 87-100 |
Kljun N, Calanca P, Rotach M W, et al. 2004. A simple parameterisation for flux footprint predictions. Bound Layer Meteor, 112(3): 503-523 |
Leclerc M Y, Thurtell G W. 1990. Footprint prediction of scalar fluxes using a Markovian analysis. Bound Layer Meteor, 52(3): 247-258 |
Leclerc M Y, Shen S H, Lamb B. 1997. Observations and large-eddy simulation modeling of footprints in the lower convective boundary layer. J Geophys Res, 102(D8): 9323-9334 |
Lee X. 2003. Fetch and footprint of turbulent fluxes over vegetative stands with elevated sources. Boundary Layer Meteor, 107(3): 561-579 |
Masseroni D, Corbari C, Mancini M. 2014. Validation of theoretical footprint models using experimental measurements of turbulent fluxes over maize fields in Po Valley. Environ Earth Sci, 72(4): 1213-1225 |
Schmid H P. 2002. Footprint modeling for vegetation atmosphere exchange studies: A review and perspective. Agri Forest Meteor, 113(1-4): 159-183 |
Schuepp P H, Leclerc M Y, MacPherson J I, et al. 1990. Footprint prediction of scalar fluxes from analytical solutions of the diffusion equation. Bound Layer Meteor, 50(1-4): 355-373 |
Smeets C J P P, van den Broeke M R. 2008. Temporal and spatial variations of the aerodynamic roughness length in the ablation zone of the Greenland ice sheet. Bound Layer Meteor, 128(3): 315-338 |
Vesala T, Kljun N, Rannik Ü, et al. 2008. Flux and concentration footprint modelling: State of the art. Environ Pollute, 152(3): 653-666 |
Wiernga J. 1993. Representative roughness parameters for homogeneous terrain. Bound-Layer Meteor, 63(4): 323-363 |