2. 广东省广州市气象台,广州 510080
2. Guangzhou Meteorological Observatory of Guangdong Province, Guangzhou 510080
冰雹是广东省较为典型的灾害性强对流天气之一,尤其是大冰雹风暴单体,具有发展迅猛、破坏力大、局地突发性强的特点,常给广东省工农业生产、国防建设及人民生命财产带来严重的危害。广东省的冰雹天气主要发生在春季,3—5月占全年总数的88%;冰雹的形成受地形影响显著,其分布有地域性,广东省以粤北山区和珠江三角洲西北部的河网地带为多发区[1]。2005年5月1日韶关、乐昌、始兴、仁化突遭强冰雹袭击,最大冰雹直径为50~60 mm,受灾人口3万,受损房屋3000多间,直接经济损失近亿元。
天气雷达是冰雹探测、预报和预警的重要工具。与传统雷达相比,多普勒天气雷达的发射功率和接收机灵敏度有了很大提高,因而能探测到传统雷达不能探测到的一些大冰雹特征,如三体散射特征、旁瓣回波特征和有界弱回波区等[2-3]。国内外一些专家学者针对强雹暴的雷达特征进行了研究。Browning等[4]对Fleming雹暴进行了分析,给出了一个超级风暴单体中冰雹增长的3个阶段:弱的上升气流区,小粒子首次上升中增长;小粒子在主上升气流前沿附近的弱上升气流中穿行,增长成为几毫米大小的胚胎,并进入主上升气流核心区;进入主上升气流的雹胚在单次上、下轨迹中增长为冰雹。1987年Zrnic[5]发现了三体散射现象,称其为三体散射特征 (TBSS),1988年Wilson等[6]称其为火焰回波,1998年Lemon[7]称其为三体散射长钉。研究表明:三体散射长度通常小于15 km,其径向速度具有零或朝向雷达的低速度值,谱宽很大,三体散射的出现是大冰雹形成的充分条件。随着我国多普勒天气雷达在业务中逐步投入使用,国内学者也针对大冰雹风暴单体的雷达特征开展了研究工作。朱君鉴等[8]利用CINRAD/SA产品分析了2002年9月27日山东省东阿县的一次大冰雹风暴单体发生、发展过程,观测到了风暴单体低层前端入流缺口、有界弱回波区和中气旋等超级单体特征。廖玉芳等[9]分析了2002年5月14日发生在湖南省常德地区的一次超级风暴单体,在国内首次利用多普勒天气雷达探测到了三体散射特征。朱敏华等[10]分析了强雹暴三体散射的多普勒天气雷达基本反射率因子的形成机制,讨论了影响三体散射观测的主要因素。王令等[11]分析了北京地区32次降雹过程的多普勒天气雷达径向速度特征,认为大风区、中气旋是伴随的常见的速度图像特征。王华等[12]对比分析了2005年北京城区两次强冰雹天气的高低空急流配置、能量和风的垂直切变,以及超级单体特征。廖玉芳等[13]针对6次强对流事件,分析了S波段多普勒天气雷达的旁瓣回波特征,这6次强对流事件均产生了大冰雹。谢健标等[14]分析了广东2005年3月22日一次飑线天气过程中的雹暴雷达特征。王秀明等[15]利用WRF中尺度数值模式模拟了2005年北京地区的一次强冰雹过程,模拟的一个长生命史的强对流雹云的演变及风暴单体结构与北京多普勒天气雷达观测的沿城市中轴线的雹云相似,具有回波墙-有界弱回波区-悬垂回波结构、对峙的倾斜上升和下沉气流等特征。向玉春等[16]应用探空资料和不同时次的降雹点的LAPS分析场作为三维对流云模式的初始场对2008年7月27日、28日湖北西部山区冰雹天气过程进行模拟,结果表明:LAPS输出场具有时空上的优势,能更好地模拟出午后局地降雹。
从2000年开始,广东省逐步建设完成了多普勒天气雷达探测网,对冰雹等中小尺度灾害性天气的发生发展进行了有效监测。本文收集近年来广东省大冰雹个例资料,首先分析了大冰雹风暴单体的垂直积分液态含水量 (VIL)、45 dBZ回波高度和反射率因子垂直递减率等特征;并逐一分析了强雹暴三体散射和旁瓣回波特征,以及0℃和-20℃层高度的回波强度分布情况,对此分析了大冰雹风暴单体与非冰雹风暴单体的反射率因子垂直廓线特征。
1 大冰雹风暴单体个例及资料简介收集广东省近年来有冰雹灾情记录报告和多普勒天气雷达观测资料的大冰雹 (冰雹直径不小于20 mm) 个例,主要包括两个步骤:① 根据冰雹灾情报告,确定近年来广东省大冰雹个例,删去一般冰雹个例 (直径小于20 mm),删去一些冰雹尺寸记录不明确的信息 (如冰雹直径为拇指大小等);② 在上述基础上,利用多普勒天气雷达观测资料进一步分析验证大冰雹可能出现的地点及时间,允许冰雹灾情报告与雷达观测到的风暴单体在时间和空间上存在一定偏差,即在冰雹观测记录地点附近10 km内应存在相对应的风暴单体,且设定1 h的时间窗,删去冰雹灾情报告中时间、地点与雷达观测到的风暴单体出入较大的个例,并将一个风暴单体形成的多个冰雹观测记录进行合并。
经过上述两个步骤处理后,共收集到12个大冰雹风暴单体 (表 1,时间为北京时,下同)。个例较少,主要原因在于广东省冰雹多出现于人烟稀少的粤北山区,一些大冰雹没有被记录到灾情报告中;同时广东省多普勒天气雷达观测网的建设逐渐展开,河源、肇庆、深圳等雷达建设相对较晚,使得一些有记录的大冰雹个例缺少了相应的雷达观测资料。
|
|
表 1 12个大冰雹风暴单体个例 Table 1 12 severe hailstorm cases |
2 大冰雹风暴单体雷达因子统计
依据Fleming超级风暴单体模型和Raymer多单体雹暴模型,发生大冰雹的潜势与风暴单体强度直接相关,也就是在很大程度上取决于上升气流的强度和尺度。当冰雹穿过宽阔的包围在强烈的上升气流核区边缘的中等强度上升气流区时,对形成大冰雹最为有利。因此,从雷达观测角度来讲,大冰雹通常与大片强雷达回波相联系,且与强回波的三维结构、尤其是垂直结构特征密切相关。为此,垂直积分液态含水量、最大反射率因子及其高度、反射率因子核区垂直厚度等因子可以作为判别强雹暴潜势的指标。
一般地,风暴单体内在0℃层高度出现45 dBZ以上强回波时,如果单体内上升气流强盛,强回波将进一步向上发展,达到-20℃层高度,则有利于大冰雹形成。为此,联合最大反射率因子及其高度、回波顶高和45 dBZ回波高度可以较好地分析风暴单体内回波强度及其垂直发展高度,并在一定程度上能反映出风暴单体内垂直上升气流的强度。
计算大冰雹风暴单体在其强盛阶段的最大反射率因子及其高度、回波顶高和45 dBZ回波高度。表 2给出了12个大冰雹风暴单体在其强盛时的上述参数值:风暴单体内最大反射率因子均超过60 dBZ,最大为73 dBZ,最小为61 dBZ (2004年3月30日花都狮岭风暴单体距离雷达站较远,超过150 km),平均为69.5 dBZ,超过65 dBZ的占83.3%。最大反射率因子对应的高度基本上都达到了5 km (个例2除外,为4.6 km),最大高度为9.2 km。10个个例45 dBZ回波高度超过10 km,另外两个个例也非常接近 (为9.7 km),最大值为16.6 km。对于个例11大冰雹风暴单体,最大反射率因子为73 dBZ,对应高度为7.7 km,此时回波顶高为20.2 km,45 dBZ回波高度也达到15.7 km,说明该风暴单体发展极为旺盛。
|
|
表 2 12个大冰雹风暴单体最大反射率因子及其高度、回波顶高、45 dBZ对应高度、VIL和VIL密度 Table 2 The maximum reflectivity and heights, echo tops, 45 dBZ echo heights, VIL and VIL density of 12 severe hailstorms |
另外,垂直积分液态水含量 (VIL) 是在假设所有反射率因子均由液态水滴引起的前提下,对其进行垂直累积,得到在某一确定底面积的垂直柱体内的液态水的总估量,反映了风暴单体的综合强度,对于大冰雹的潜势具有较好的指示作用。Amburn等[17]定义风暴单体VIL与顶高之比为VIL密度,研究表明,如果VIL密度超过4 g·m-3,则风暴单体几乎都会产生大冰雹。为此,表 2还给出了12个大冰雹风暴单体强盛时的VIL及其密度,其中,VIL来自于多普勒天气雷达PUP中的风暴结构 (SS) 产品。结果表明:VIL最大值达到了91 kg·m-2,最小值为32 kg·m-2,平均值为69.3 kg·m-2。仅有1个个例的VIL值小于50 kg·m-2(2004年3月30日,花都狮岭),分析发现该风暴单体距离雷达站约155 km,回波底高为3.2 km,这可能是造成VIL偏小的主要原因。VIL密度普遍为3~5 g·m-3, 平均为4.1 g·m-3,离差系数仅为0.16,VIL密度最小值为2.6 g·m-3,对应上述花都狮岭风暴单体。
对于广东省,一些在海洋上发展的风暴单体移入陆地,虽然其回波也较强,最大反射率因子能够达到50 dBZ以上,但由于其垂直发展不够旺盛,强回波多集中于中低层,尤其是强回波中心垂直向上衰减极为迅速,因此很难形成大冰雹。
针对12个大冰雹风暴单体个例,计算强回波核区垂直伸展厚度和风暴单体内反射率因子垂直递减率。定义50 dBZ以上回波区为风暴单体核心区,强核区垂直伸展厚度是指风暴单体内50 dBZ回波强度所能达到的最高高度与最低高度之间的差值。另外,考虑不同风暴单体之间以及同一风暴单体不同发展阶段之间垂直结构的差异,没有选择某一固定高度层计算风暴单体内反射率因子垂直递减率,而是选择在风暴体内最大反射率因子上方3 km厚度层内计算。表 3给出了12个大冰雹风暴单体强盛时的核心区垂直伸展厚度和核心以上3 km厚度层内的反射率因子垂直递减率。核心区垂直伸展厚度最大值为14.8 km,平均约为10 km,最小值为4.7 km (2004年3月30日花都狮岭)。垂直递减率最大值为5.1 dB·km-1,最小值为0.5 dB·km-1,平均值为2.4 dB·km-1。
|
|
表 3 产生大冰雹的风暴单体的核心区垂直伸展厚度 Table 3 The reflectivity core depths of 12 severe hailstorms |
3 三体散射、旁瓣回波和环境温度层回波特征
对于S波段多普勒天气雷达,出现三体散射特征是风暴单体存在大冰雹的充分条件,但当风暴单体周围有大片雷达回波时,往往无法观测到三体散射现象。另外,按照Fleming超级单体冰雹增长模型,在主上升气流前沿附近的弱上升气流中穿行的小粒子增长为雹胚,并进入主上升气流区,在其上、下轨迹中增长为冰雹。在这个过程中,除上升气流外,合适的环境温度层高度非常重要,一是需要大量过冷水滴的存在,二是大粒子到达穹窿顶部的平衡层时,如果仍处于过冷却水区,则会明显增长。也就是说,风暴单体内强核区向上扩展到0℃等温线以上才能对降雹的潜势有所贡献,当强回波扩展到-20℃高度层之上时,对强降雹的潜势贡献最大。
3.1 三体散射和旁瓣回波特征在12个大冰雹风暴单体中,有6个风暴单体观测到了三体散射特征,有3个风暴单体出现了旁瓣回波现象,其中1个风暴单体既观测到三体散射又观测到旁瓣回波。有4个风暴单体没有观测到三体散射或旁瓣回波特征,主要原因在于这些风暴单体不是孤立的,在其周围有明显的、较大范围的回波 (表 4)。
|
|
表 4 三体散射和旁瓣回波特征 Table 4 Features of three body scatter spike (TBSS) and side-lobe echoes |
进一步分析发现:几乎所有的三体散射特征和旁瓣回波特征均出现在大冰雹降落记录之前,具有一定的预报提前量;三体散射和旁瓣回波特征一般能维持较长一段时间,具备一定的连续性。在12个大冰雹风暴单体中,出现三体散射或旁瓣回波特征至观测到大冰雹这段预报提前时间最大可达40 min,平均为14 min,仅1个风暴单体样本没有预报提前量。三体散射和旁瓣回波特征的持续时间一般为18~36 min,平均为32.2 min。三体散射特征通常先出现在较高的仰角层次上,相应的最大反射率因子均超过了65 dBZ。
3.2 0℃与-20℃层高度回波分布特征一种简单有效的判别大冰雹是否存在的方法是讨论0℃和-20℃层高度回波强度的分布情况。考虑到在获取0℃和-20℃层高度时,气候统计值仅为一种平均态,且广东省仅有4个高空站,每天探测两次,其时空分辨率不理想。为此,利用美国国家环境预报中心 (NCEP) 提供的FNL全球分析资料 (空间分辨率为1°×1°,时间分辨率为6 h,垂直分层28层),针对大冰雹风暴单体实例,计算相近时次0℃和-20℃层的高度,进而利用广东省多普勒天气雷达探测数据,制作0℃和-20℃层高度回波强度分布产品。图 1给出了12个风暴单体在-20℃层高度的回波强度分布。
|
|
| 图 1. 大冰雹风暴单体在-20℃层高度上的回波强度分布 Fig 1. The reflectivity distribution at-20℃ layer of 12 hailstorms | |
表 5给出了12个大冰雹风暴单体对应的0℃和-20℃层高度以及该高度上风暴单体内最大反射率因子。0℃层高度最低为4020 m,最高为4820 m,平均为4502 m;-20℃层高度介于7040~8030 m之间,平均为7682 m。风暴单体在相应等温层高度上,回波强度较强:在0℃层高度上,风暴单体内最大反射率因子为54~67 dBZ,平均值为62 dBZ;在-20℃层高度上,回波强度仍能达到54 dBZ以上,最大为66 dBZ,平均为62.5 dBZ。对于个别风暴单体,尽管其所处的周围环境温度层较高,如个例10(连山) 和个例11(韶关市区),0℃和-20℃层高度分别达到4700 m和8000 m以上,但在该高度风暴单体内仍存在较强的回波,0℃和-20℃高度的最大反射率因子均超过了60 dBZ,因此风暴单体垂直发展非常旺盛,非常有利于大冰雹的生成和增长。
|
|
表 5 大冰雹风暴单体对应的0℃和-20℃层高度及该高度上的最大反射率因子 Table 5 Average heights of 0℃ and-20℃ layers of 12 hailstorms with their maximum echoes |
4 大冰雹风暴单体与非冰雹风暴单体反射率因子垂直廓线特征对比
前面通过分析大冰雹风暴单体的垂直积分液态水含量、45 dBZ回波高度、强核区垂直伸展厚度和垂直递减率等因子,在一定程度上反映了强雹暴的垂直结构信息。为了更清晰地分析风暴单体垂直结构特征,利用广东省多普勒天气雷达基数据进行垂直插值,以获取风暴单体内反射率因子的垂直分布廓线[18]。同时,为进一步将大冰雹风暴单体与普通风暴单体区分开来,还挑选了46个非冰雹风暴单体个例进行对比分析。在选择非冰雹风暴单体时,一方面,为对比分析与强雹暴相类似的环境条件下的普通风暴单体的垂直廓线,选择了冰雹日中的非冰雹风暴单体 (21个个例);另一方面,为分析不利于大冰雹发生的环境条件下的普通风暴单体特征,还随机选择了非冰雹日的风暴单体 (25个个例)。在选择冰雹日的非冰雹风暴单体时,要求风暴单体距离冰雹观测记录位置超过50 km,或时间相差2 h以上。
4.1 风暴单体反射率因子垂直廓线对比图 2、图 3和图 4分别给出了大冰雹风暴单体、冰雹日非冰雹风暴单体和非冰雹日风暴单体的反射率因子垂直分布廓线。
|
|
| 图 2. 大冰雹风暴单体反射率因子垂直分布 Fig 2. The vertical reflectivity distribution of hailstorms | |
|
|
| 图 3. 冰雹日非冰雹风暴单体的反射率因子垂直分布 Fig 3. The vertical reflectivity distribution of non-hail storms in hailstorm days | |
|
|
| 图 4. 非冰雹日风暴单体反射率因子垂直分布 Fig 4. The vertical reflectivity distribution of storms in non-hail days | |
对于大冰雹风暴单体 (图 2),最大反射率因子的极大值超过70 dBZ,极小值也达到60 dBZ以上,且绝大多数的最大反射率因子超过65 dBZ;最大反射率因子对应的垂直高度主要集中在5~10 km。对于所有大冰雹风暴单体,在10 km高度以下,几乎所有观测到的回波强度均超过40 dBZ;一些风暴单体内强度达到30 dBZ的回波可以垂直伸展到15 km以上。表明风暴单体内垂直运动非常强盛,大的水凝物能输送到很高的高度,非常有利于大冰雹的生成和增长。
对于冰雹日中的非冰雹风暴单体 (图 3),最大反射率因子一般为60 dBZ左右,极大值接近70 dBZ,而极小值仅为50 dBZ;最大反射率因子对应的高度相对于大冰雹风暴单体而言较低,多位于2~5 km高度层内。对于这类风暴单体,40 dBZ回波一般很难达到10 km高度以上,主要集中在5~8 km高度范围内,且向上垂直递减率较大,与大冰雹风暴单体垂直发展旺盛有明显差异。即在环境条件与大冰雹风暴单体相似的前提下,由于这类风暴单体的强核区达到的垂直高度不够,且向上垂直伸展有限,因此难以形成大冰雹。
对于非冰雹日风暴单体 (图 4),由于样本选择的随机性,风暴单体内最大反射率因子及其垂直伸展变化范围较大。最大反射率因子为55~68 dBZ,相应的垂直高度主要集中在5 km高度以下。有些风暴单体内40 dBZ的回波强度向上伸展到约15 km高度,有多个高度超过10 km;当然也有半数左右的风暴单体,其体内40 dBZ的回波垂直伸展有限。这类风暴单体没有出现大冰雹,一方面,由于部分风暴单体自身强度及对流发展不够;另一方面,对于那些发展较为强盛的风暴单体,由于所处的周围环境条件不利,也未形成大冰雹。
4.2 不同类型风暴单体的平均垂直廓线图 5给出了大冰雹风暴单体、冰雹日非冰雹风暴单体和非冰雹日风暴单体的平均垂直廓线分布。大冰雹风暴单体的平均廓线与其他风暴单体的平均廓线有明显区别:大冰雹风暴单体的低层回波从55 dBZ迅速增加到65 dBZ,在3~8 km高度层有一个较为深厚的强核心区,在该高度层内,回波垂直递减率很小;在8 km高度以上,回波强度开始迅速减弱,但仍较其他风暴单体偏强。对于冰雹日非冰雹风暴单体和非冰雹日风暴单体,反射率因子平均垂直分布廓线较为接近,低层回波由50 dBZ左右增至最强约60 dBZ,强核心区高度为3~5 km,在6 km高度以上回波强度迅速递减。分析表明,不考虑环境因素,具有较大强度、较高高度和深厚发展的强回波核心区是大冰雹风暴单体与其他两类非冰雹风暴单体的主要差异。
|
|
| 图 5. 3种风暴单体反射率因子平均垂直廓线 Fig 5. Average reflectivity profiles of hailstorms and non-hail storms | |
5 小结
利用广东省多普勒天气雷达资料,对12个大冰雹风暴单体的雷达特征进行分析,得到以下结论:
1) 风暴单体内最大反射率因子均超过60 dBZ,最大为73 dBZ,最小为61 dBZ,平均为69.5 dBZ,超过65 dBZ的占83.3%。最大反射率因子对应的高度基本达到5 km,最大高度为9.2 km。45 dBZ回波高度超过9.7 km。
2) VIL最大值达91 kg·m-2,最小值为32 kg·m-2,平均值为69.3 kg·m-2;VIL密度一般为3~5 g·m-3, 平均为4.1 g·m-3。核心区垂直伸展厚度最大值为14.8 km,平均约为10 km;反射率因子垂直递减率的平均值为2.4 dB·km-1。
3) 在12个大冰雹风暴单体中,有6个观测到了三体散射特征,有3个出现了旁瓣回波现象,其中1个风暴单体既观测到三体散射又观测到旁瓣回波。出现三体散射或旁瓣回波特征的大冰雹预报提前时间最大可达40 min,平均为14 min。
4) 0℃层高度最低为4020 m,最高为4820 m,平均为4502 m;-20℃层平均高度为7682 m。在0℃层高度上,最大反射率因子为54~67 dBZ,平均为62 dBZ;在-20℃层高度上,回波强度仍能达到54 dBZ以上,最大为66 dBZ,平均为62.5 dBZ。
5) 大冰雹风暴单体与非冰雹风暴单体反射率因子垂直廓线存在明显差异:大冰雹风暴单体低层回波从55 dBZ迅速增强到65 dBZ,强核区垂直伸展较厚,位于3~8 km高度层上;非冰雹风暴单体低层回波从50 dBZ, 加强到60 dBZ, 核心区相对较低,为3~5 km, 6 km高度以上雷达回波迅速减弱。
| [1] | 林良勋, 冯业荣, 黄忠, 等. 广东省天气预报技术手册. 北京: 气象出版社, 2006. |
| [2] | 俞小鼎, 姚秀萍, 熊廷南, 等. 多普勒天气雷达原理与业务应用. 北京: 气象出版社, 2006: 3–15. |
| [3] | 刘黎平, 葛润生. 中国气象科学研究院雷达气象研究50年. 应用气象学报, 2006, 17, (6): 682–689. DOI:10.11898/1001-7313.20060606 |
| [4] | Browning K A, Foote G B. Airflow and hail growth in supercell storms and some implications for hail suppression. Quart J R Meteo Soc, 1976, 102: 499–533. DOI:10.1002/(ISSN)1477-870X |
| [5] | Zrnic D S. Three-body scattering produces precipitation signature of special diagnostic value. Radio Sci, 1987, 22: 76–86. DOI:10.1029/RS022i001p00076 |
| [6] | Wilson J W, Reum D. The flare echo:Reflectivity and velocity signature. J Atmos Oceanic Technol, 1988, 5: 197–205. DOI:10.1175/1520-0426(1988)005<0197:TFERAV>2.0.CO;2 |
| [7] | Lemon L R. The radar "three-body scatter spike":An operational large-hail signature. Wea Forecasting, 1998, 13: 327–340. DOI:10.1175/1520-0434(1998)013<0327:TRTBSS>2.0.CO;2 |
| [8] | 朱君鉴, 刁秀广, 黄秀韶. 一次冰雹风暴的CINRAD/SA产品分析. 应用气象学报, 2004, 15, (5): 580–589. |
| [9] | 廖玉芳, 俞小鼎, 郭庆. 一次强对流系列风暴个例的多普勒天气雷达资料分析. 应用气象学报, 2003, 14, (6): 656–662. |
| [10] | 朱敏华, 俞小鼎, 夏峰, 等. 强烈雹暴三体散射的多普勒天气雷达分析. 应用气象学报, 2006, 17, (2): 215–223. DOI:10.11898/1001-7313.20060213 |
| [11] | 王令, 郑国光, 康玉霞, 等. 多普勒天气雷达径向速度图上的雹云特征. 应用气象学报, 2006, 17, (3): 281–287. |
| [12] | 王华, 孙继松, 李津. 2005年北京城区两次强冰雹天气的对比分析. 气象, 2007, 33, (2): 49–56. DOI:10.7519/j.issn.1000-0526.2007.02.008 |
| [13] | 廖玉芳, 俞小鼎, 吴林林, 等. S波段多普勒天气雷达旁瓣回波的特征分析. 热带气象学报, 2008, 24, (2): 183–188. |
| [14] | 谢健标, 林良勋, 颜文胜, 等. 广东2005年"3.22"强飑线天气过程分析. 应用气象学报, 2007, 18, (3): 321–329. |
| [15] | 王秀明, 钟青, 韩慎友. 一次冰雹天气强对流 (雹) 云演变及超级单体结构的个例模拟研究. 高原气象, 2009, 28, (2): 352–365. |
| [16] | 向玉春, 杨军, 唐仁茂, 等. LAPS在冰雹云模拟中的应用. 应用气象学报, 2012, 23, (3): 331–339. |
| [17] | Amburn S A, Wolf P L. VIL density as a hail indicator. Wea Forecasting, 1997, 12: 473–478. DOI:10.1175/1520-0434(1997)012<0473:VDAAHI>2.0.CO;2 |
| [18] | 张家国, 万玉发, 王珏. 风暴生命史雷达特征量反演. 应用气象学报, 2008, 19, (1): 101–105. DOI:10.11898/1001-7313.20080116 |
2015, 26 (1): 57-65



