基于CLM4.5土壤风蚀模型的河北省土壤释尘量模拟

王小伟 ，王 卫

（1. 河北师范大学资源与环境科学学院，河北省环境演变与生态建设实验室，石家庄 050016；2. 河北女子职业技术学院，石家庄 050051）

0 引 言

随着经济社会不断发展，中国大气污染问题日益严重。研究表明，由地表风蚀而注入到大气中的土壤风蚀尘，是大气污染物的主要组分之一[1-3]。土壤风蚀起尘不但造成大气污染，还破坏了土壤结构，造成土壤物质流失，加剧了地区的荒漠化[4]。悬浮在大气中的土壤风蚀粉尘，不但危害人体健康[5-6]，还改变地面大气辐射平衡和气溶胶的酸碱度[7]，对硫酸盐和硝酸盐的输送和沉降有着重要的影响[8-9]；长距离输送而注入海洋的土壤风蚀粉尘还影响了海洋的生物地球化学过程和生物量，可能引起全球气候系统和生态系统的变化[10-11]。

土壤风蚀粉尘引起一系列环境问题，使很多学者开始关注土壤风蚀的释尘量研究。学术界通过室内风洞试验和野外定点观测，在风沙物理学理论支持下，提出了很多半机理模型，定量模拟土壤的释尘量。这些模型有：Gillette的风蚀粉尘释放模型[12-13]、美国EPA风蚀粉尘释放模式[14-15]、邵亚平的风蚀粉尘释放模型[16-18]、DPM风蚀粉尘释放模型[19-20]、DUSTRAN 模型[21]、Zender的DEAD模型[22]以及全球沙尘循环模型中关于沙尘释放的模拟[23-24]等。相比国外土壤风蚀粉尘模型的研究，国内学者大多是应用上述模型进行释尘量模拟研究[25-32]，也有部分开展自主的方法或模型研究[33-34]。研究结果表明，上述模型都能模拟出土壤风蚀粉尘释放的时空分布特征，但不同模型模拟出的释尘通量存在显著的差异，差异甚至达到几个数量级[35]，主要原因是不同模式在模拟释尘量时对起尘物理过程的简化程度及考虑参数各不相同。尽管上述关于土壤风蚀粉尘的模型模拟研究取得了令人鼓舞的进步[36]，但土壤风蚀粉尘模型仍需更多的区域验证，尽可能多地把地面大气状况、植被覆盖、土壤特性等参数考虑进去，以进一步改善模拟性能。

自全国 74个城市按照新标准开展空气质量监测以来，空气质量相对较差的前 10位城市中，河北省 2013年占7席，2014年占6席，2015年占7席，2016年占4席，2017年1-11月占6席，是中国大气污染最为严重的地区之一。研究表明，土壤风蚀粉尘对大气污染物PM10的贡献率通常在20%～60%之间[28]，是造成大气污染的主要因素之一。在目前观测数据十分有限的情况下，开展河北省土壤风蚀释尘量模拟研究，深入了解土壤风蚀粉尘释放的过程与机理，对防治大气污染具有重要意义。

1 研究区概况与研究方法

1.1 研究区概况

河北省位于 113°27′～119°50′E，36°05′～42°40′N 之间，陆地总面积18.85×106hm2。地貌类型复杂多样，从省域西北向东南依次为坝上高原、燕山和太行山地、河北平原三大地貌单元。属于温带大陆性季风气候，春干、夏雨、秋爽、冬寒；年降水量大体从省域东南向西北由多到少，并在燕山和太行山迎风坡形成降水中心；年大风日数的高值区包括省域西北部的坝上高原、冀西北间山盆地和东部的滨海平原，其他区域属年大风日数的相对低值区。受降水少等地带性因素的影响，省域西北部的坝上高原和冀西北间山盆地的植被覆盖度相对较低；受土壤盐碱化等非地带性因素的影响，省域东部滨海平原特别是沧州滨海平原的植被覆盖度也相对较低；其他区域的植被覆盖度相对较高[37]。

1.2 CLM4.5土壤风蚀模型

CLM4.5土壤风蚀模型是CLM(community land model)的一个子模型，其起尘模式是基于 DEAD（dust entrainment and deposition）模型改进的。沙尘粒子的运动主要受重力和空气动力学拖曳力的影响[38]，当拖曳力大于重力时，沙尘粒子脱离地面上升到某一高度，并随大气运动。垂直起尘通量Fj（kg/(m2·s)）的计算公式[39]如下

式中 T是用于平衡模型空间和时间分辨率的调整因子,用来得到一个较合理的气候模拟。S是土壤侵蚀度因子,它与土壤遭受风蚀的敏感性成正比[40]，位于0～1之间。在CLM4.5土壤风蚀模型中，考虑到模型适用于全球范围，为了得到较好的模拟效果，将气候模拟调整因子 T和土壤侵蚀度因子S统一分别赋值为5×10-4和1。而本研究的时间分辨率和空间分辨率均高于CLM4.5模型，为了提高模拟精度，采用Wu等[41]确定的高分辨率的空间调整参数T=1.32×10-2和S=0.02来修正土壤风蚀的释尘量。fm是适宜起尘的裸露土壤占栅格的比例；α是起沙质量功率；Qs是沙尘颗粒的水平跃移通量，kg/(m·s)；Mi,j是第i种尘源模式携带的第j种传输粒径组的质量分数。裸露土壤占栅格的比例fm计算公式如式（2）所示。

式中flake、fwelt、fsno分别代表栅格中湖泊/水库、湿地和积雪覆盖的比例；fv指植被覆盖率，位于0～1之间，具体计算公式见公式（3）。

式中L指叶面积指数，s指茎面积指数，(L+s)t是常数，等于0.3。

公式（2）中的 wliq,1和wice,1分别是表土层中液态水和冰的厚度，m；wliq,1/( wliq,1+wice,1)指地表冰冻比率。土壤中水层厚度等于土层厚度与体积含水量的乘积。冰层厚度等于土层厚度与体积含冰量的乘积。土壤含冰量是基于平衡态热力学理论，按照土壤水势和温度的关系式以及土壤水力学特征关系来确定的[42]。具体见公式（4）和（5）。

式中θmax表示负温状态下的最大液态含水率；θs为土壤孔隙率，Li为熔化潜热，J/kg；Tf为水的冻结温度，K；φs为土壤饱和基质势，Pa；b为Clapp & Hornberger常数，g为重力加速度，m/s2。当土壤冻结时，多余的液态含水率即为土壤含冰率θi。

式中的ρw和ρi分别代表水和冰的密度，kg/m3。

公式（1）中的起沙质量功率 α的计算公式如（6）所示。

式中的Mclay指土壤中黏粒的质量分数。

公式（1）中的水平跃移通量Qs计算公式为

式中跃移常数Cs=2.61，ρatm是大气密度，kg/m3；u*t是摩阻风速临界值，m/s；u*s是基于欧文效应的摩阻风速，m/s。u*s的计算公式为

式中 u*是通常意义上的摩阻风速，m/s；U10是距离地面10 m高的风速，m/s；U10,t是10m高的风速临界值，m/s。u*和U10,t计算公式见公式（9）和公式（10）。

式中K是卡曼常数，等于0.4；Z是高度，m；U是Z高度的平均风速，m/s；Z0是地表粗糙度，m。

公式（7）中的摩阻风速临界值u*t的计算公式为

式中w指土壤的质量湿度，wt是土壤湿度的临界值。

式中mi是第i种尘源模式的质量分数；是第i种尘源模式的质量中位粒径，m；σg,i是第i种尘源模式的几何标准差，具体数值见表1；Dj,min和Dj,max是第j种传输粒径组的最小和最大粒径，m；具体数值见表2。

表1 尘源模式i的质量分数、质量中位粒径和几何标准差Table 1 Mass fraction, mass median diameter, and geometric standard deviation, per dust source mode i

表2 第j种传输粒径组的最小和最大粒径Table 2 Minimum and maximum particle diameters in each dust transport group j

1.3 技术路线

图1 土壤风蚀粉尘模型流程Fig.1 Flow chart of soil wind-blown dust model

2 数据来源及处理

2.1 数据来源

风速和土壤湿度来源于中国气象局陆面数据同化系统（CLDAS），风速是10 m风速；土地覆被分类数据采用2013年的IGBP的全球土地覆被分类方案；植被覆盖度中的叶面积指数采用全球陆表特征参量中 GLASS产品；湖泊、河流、湿地数据来源于河北省土地2010年土地利用图；积雪覆盖比率数据采用 MOD10A2产品；冰冻比率的计算采用中国气象科学共享网的地温数据；土壤数据来源于FAO的世界土壤数据库（HWSD）。

2.2 数据处理

2.2.1 风 速

研究选取CLDAS的2013年的10 m风速数据，时间分辨率为逐小时，空间分辨率是0.05°，数据为NC格式。在Arcgis10.2的技术支撑下，通过Model builder批量完成数据格式转换、重投影、重采集和研究区裁切等操作，转换成1 km分辨率的河北省风速网格数据。

2.2.2 土壤湿度

研究选取CLDAS土壤湿度数据中的0～5 cm土壤体积湿度数据，在 Arcgis10.2的技术支撑下，通过 Model builder对2013年土壤体积湿度数据进行格式转换、重投影、重采集和研究区裁切等操作，再结合液态水的密度和表层土壤的体积密度，确定河北省1 km分辨率的表层土壤的质量含水率(w)，通过与表层土壤质量含水率的临界值(wt)进行比较，确定公式（11）中的fw。

2.2.3 地表粗糙度

地表粗糙度反映了地表粗糙元对风力的削弱程度，是阻碍土壤风蚀粉尘产生的重要因子。在对 2013年的IGBP的全球土地覆被分类方案进行裁剪的基础上，按照植被/土地覆盖类型与地表粗糙度的关系[43]来确定1 km分辨率的河北省不同土地覆被类型的地表粗糙度网格数据。

2.2.4 植被覆盖度

研究采用我国自主研发的覆盖全球的GLASS LAI数据，空间分辨率为1 km，时间分辨率为8 d，利用MRT

由CLM4.5土壤风蚀模型可知，影响土壤风蚀粉尘释放的因素包括土壤质地、风速、土壤湿度、地表粗糙度、裸露土壤比例等因素。裸露土壤比例由植被覆盖度、地表冰冻比率、积雪覆盖比率、湖和湿地比率决定。土壤质地、地表粗糙度、湖泊、河流和湿地比率属于相对稳定因素。风速、土壤湿度、植被覆盖率、积雪覆盖面积、冰冻比率等数据随时间变化较大。各影响因素对土壤风蚀粉尘的作用见土壤风蚀粉尘模式流程（图1）。进行图形拼接和投影转换，在Arcgis10.2的技术支撑下，将2013年的GLASS LAI数据按照河北区域进行裁剪。茎面积指数的确定按照陆面模型中的动态植被模型（CLM_DGVM）的方法[44]确定，最后根据每8 d的fv数据线性插值成每天的植被覆盖度数据。

2.2.5 湖泊、河流和湿地的比率

湖泊、河流和湿地的比率指栅格中湖泊、河流和湿地面积占栅格总面积的比率。湖泊、河流和湿地的面积随时间变化不大，因此采用2010年MODIS和TM两种遥感数据为主要数据源，以 NDVI作为特征参数，采用CART决策树对影像进行土地覆被分类，确定河北区域内的湖泊、河流和湿地面积，再确定栅格内湖泊、河流和湿地比率，即公式（2）中的flake、fwelt。

2.2.6 积雪覆盖比率

积雪覆盖面积的确定采用MODIS数据MOD10A2产品，空间分辨率为500 m，数据格式为HDF。研究区域在SIN投影方式下每8天有4张图，选取2013年1-3月及10-12月共计96张图，利用MRT软件进行图像拼接与坐标转换，然后利用Arcgis10.2裁剪出河北区域，采用积雪分类编码提取积雪面积，确定栅格的积雪覆盖比率fsno。

2.2.7 地表冰冻比率

选取2013年气象共享网的地温资料，利用克里金插值确定河北省地温网格数据，用每日2点、8点、14点和20点地温监测数据计算地温的日均值，结合依据土壤质地确定的φs，b数据和CLDAS的土壤体积湿度，按照公式（4）和公式（5）计算每日土壤中的最大含水量及含冰量，以确定地表冰冻比率。

2.2.8 土壤质地

土壤质地数据来源于HWSD，数据格式为Grid格式，投影方式为WGS84，空间分辨率为1 km。从世界土壤数据库中裁剪出河北区域，然后在土壤属性表中选取需要的黏粒含量、沙粒含量、有机质比例等因子进行计算，确定区域的起沙质量功率（α）和土壤的质量含水率(w)。

3 结果与分析

运用CLM4.5土壤风蚀模型，对2013年河北省的土壤风蚀释尘量进行模拟估算：2013年河北省土壤风蚀粉尘 PM10的释尘通量为1.02 t/hm2，全年共向大气中释放PM101 922.87×104t 。

3.1 释尘量时空分布特征

3.1.1 释尘强度分级特征

2013年河北省土壤风蚀释尘强度存在明显差异，按照增强时空分布差异区分度和各级别面积比例适度的原则，将河北省土壤风蚀的释尘强度划分为 0～0.20、0.20～0.50、0.50～1.00、1.00～5、5～10 及＞10 t/(hm2·a) 6个等级（表 3），其中Ⅵ级为释尘强度最高的等级。河北省释尘强度最高的Ⅵ级区域总释尘量为571.21× 104t，占全省总释尘量的29.71%，面积占全省面积的1.42%；释尘强度Ⅴ级区域总释尘量为622.70×104t，占总释尘量的32.28%，面积仅占4.82%；释尘强度第3高的Ⅳ级区域总释尘量为610.24×104t，占释尘总量的31.74%，面积占比为12.76%；释尘强度Ⅲ级区域总释尘量为61.45×104t，占释尘总量的3.20%，面积占全省面积的4.52%；释尘强度Ⅱ级区域总释尘量为 34.33×104t，占全省总释尘量的1.79%；释尘强度最低的Ⅰ级区域总释尘量仅为22.94×104t，占释尘总量的1.19%，面积却占全省总面积的71.10%。其中Ⅴ、Ⅵ级区域总释尘量为1 193.91×104t，占到总释尘量的62%，而面积占比仅为6.24%，由此可见，河北省土壤风蚀粉尘释放属于小区域、高强度的排放模式。

表3 河北省土壤风蚀粉尘PM10释尘强度分级Table 3 Emissions intensity grades of wind-blown dust PM10 in Hebei Province

3.1.2 时间分布特征

1）月际分布特征

河北省2013年1-12月土壤风蚀粉尘PM10的释尘通量 按 时 间 顺 序 分 别 为 ： 28.31×10-3、 130.98×10-3、281.14×10-3、 238.50×10-3、 26.20×10-3、 10.60×10-3、2.29×10-3、2.38×10-3、6.59×10-3、14.69×10-3、157.17×10-3、121.66×10-3t/hm2（图2）。全年12个月中，3月份释尘量最高，占全年的27.6%。植被覆盖度低、冻土和积雪融化导致的裸露土壤比例增大，再加上土壤湿度低、风力大是3月份土壤风蚀粉尘释尘量高的原因。7月份的土壤风蚀粉尘最低，仅占全年释尘量的0.2%，土壤湿度大、植被覆盖度高是释尘量低的主要原因。

图2 2013年河北省各月和各季PM10的释尘通量Fig.2 Monthly and seasonal average emission flux of PM10 in Hebei province at 2013

2）季节分布特征

2013年春季河北省大部分地区的土壤风蚀PM10的释尘通量低于1 t/hm2，坝上高原和沧州滨海平原释尘通量较高，局部地区高于10 t/hm2，区域平均释尘通量为0.55 t/hm2，占到全年的53.6%；夏季全省土壤风蚀PM10的释尘通量较低，只有沧州滨海平原和冀东滨海平原地区释尘通量高于 0.50 t/hm2，区域平均释尘通量为 0.015 t/hm2，仅占全年的1.5%；秋季全省释尘通量逐渐回升，但大部分地区的土壤风蚀 PM10的释尘通量低于 0.50 t/hm2，沧州滨海平原和冀东滨海平原地区释尘通量相对较高，但低于5 t/hm2，区域平均释尘通量为0.18 t/hm2，占全年的17.5%，远低于春季；冬季释尘通量在秋季基础上进一步升高，区域平均释尘通量为0.28 t/hm2，占全年的27.4%，是释尘通量次高的季节（图2，图3）。河北省土壤风蚀粉尘PM10的释尘量呈现明显的春高夏低季节变化。春季由于植被稀少，气温回升快，裸露地表冻土和积雪融化，再加上风力较大，使土壤风蚀释尘量最高；夏季雨水增多，植被旺盛，土壤湿度大，不利于起尘；秋季植被减少，土壤湿度下降，使秋季的释尘量回升；冬季由于土壤湿度进一步降低、植被覆盖度下降，使地面起尘能力上升，释尘量仅次于春季。

图3 河北省土壤风蚀粉尘PM10释尘通量季节变化Fig.3 Emission flux of wind-blown dust PM10 in different seasons in Hebei province

3.1.3 空间分布特征

2013年河北省土壤风蚀释尘量呈现明显的地区差异（图 4），全省高等级释尘强度（Ⅴ、Ⅵ级）区域主要位于沧州滨海平原东北部、坝上高原西北部和冀东滨海平原南部地区，中等释尘强度（Ⅳ级）区域主要位于坝上高原、冀西北间山盆地、沧州滨海平原、中部平原的东北部地区、太行山山地南部与平原交接地带；冀北山地、燕山山前平原、太行山山地、太行山山前平原、中部平原的南部地区的释尘强度都较低。按照全省自然区域统计，释尘通量最高的为沧州滨海平原，年释尘通量为5.365 t/hm2；其次为坝上高原，年释尘通量为 3.542 t/hm2；冀东滨海平原、冀西北间山盆地、中部平原、太行山山地、燕山山前平原、冀北山地的年释尘通量依次为：3.382、2.369、0.550、0.423、0.321、0.240 t/hm2。释尘量最低的为太行山山前平原，年释尘通量为0.198 t/hm2（图5）。

图4 2013年河北省土壤风蚀粉尘PM10的释尘通量Fig.4 Emission flux of PM10 from wind-blown dust in Hebei province at 2013

图5 河北省土壤风蚀粉尘PM10释尘通量区域分布Fig.5 Emission flux of wind-blown dust PM10 in different regions

3.2 释尘量成因特征

2013年河北省的年均风速、土壤湿度和裸露土壤比例如图 6所示，年均风速的高值区域分布在滨海平原和坝上高原两地，太行山山地和冀北山地的风速较小。土壤湿度的高值区域主要是沿海区域，低值区域主要是坝上高原、冀北山地西北部-冀西北间山盆地-太行山山地一线。裸露土壤比例的高值区域主要分布在冀西北间山盆地、燕山山前平原和冀东滨海平原、山前平原和太行山山地南部与平原交接地带。对比土壤风蚀释尘通量（图4）与风速、土壤湿度和裸露的土壤比例的分布规律可以发现，风速高值区域、土壤湿度低值区域和裸露土壤比例的高值区域是河北省土壤风蚀释尘强度的高等级地区，与河北省土壤风蚀释尘通量空间分布具有一致性。

图6 2013年河北省风速、土壤湿度、裸露土壤比例年均值Fig.6 Average wind speed, soil moisture and proportion of bare soil of Hebei province at 2013

为了定量分析风速、土壤湿度、裸露土壤的比例对土壤风蚀粉尘释放的影响，本研究对风速、土壤湿度、裸露土壤比例的年均值及风蚀粉尘的年释尘通量进行分区抽样，研究区共取样400个点，去除异常值后剩余371个点。利用SPSS19.0建立多元回归模型，模型的t检验和F检验均达到0.01的显著性，调整后的R2为0.230。回归模型如下：

式中X1为风速；X2为土壤湿度；X3为裸露土壤比例。由模型结果可知，土壤风蚀释尘量与风速、裸露土壤比例成正比，风速和裸露土壤比例越大，释尘量越高；与土壤湿度成反比，土壤湿度越大，释尘量越低。风速、土壤湿度和裸露土壤比例这 3个影响因素中，土壤风蚀释尘量对裸露土壤比例变化最为敏感，对风速的变化最不敏感。因此，应从减少裸露土壤比例方面入手，通过植树造林、退耕还林、实施农田保护性耕作等方式降低土壤风蚀的释尘量。

4 讨 论

随着观测手段的不断进步，沙尘模式中输入参数的精度逐渐提高，对于沙尘排放的测算和计算也更为准确。但目前针对沙尘模式的模拟验证只局限在小区域内，在较大范围内仍无法进行。在本研究区这种自然起尘与人为起尘并存的区域如何进行土壤风蚀粉尘模拟结果的定量验证，有待于进一步探索。

通过与相同区域不同模式模拟的释尘通量进行对比，本研究年均土壤风蚀释尘通量为 3.23 μg/(m2·s)，与郭俊等[31]采用GOCART模型模拟的2000～2007年的年均释尘通量 0～4.6 μg/(m2·s)之间的模拟结果相近，能够较好的模拟河北区域土壤风蚀PM10的释尘量。但由于本研究和对比研究的研究区范围和研究基期存在差异，风速、土壤湿度等影响土壤风蚀粉尘产生的因子不同，会导致模拟结果对比的随机性。因此，需要进一步深化不同模拟成果对比的规范化。

5 结 论

1）2013年河北省土壤风蚀粉尘PM10的年释尘通量值为1.02 t/hm2，全年共向大气中释放PM101 922.87×104t。春季释尘通量为0.55 t/hm2，占全年的53.6%，是一年中土壤风蚀粉尘排放量最大的季节，夏季释尘通量最低，为0.015 t/hm2。CLM4.5土壤风蚀模型可以很好的模拟出河北省土壤风蚀粉尘发生的时间变化特征。

2）2013年河北省土壤风蚀释尘量呈现明显的地区差异，释尘量最高的为沧州滨海平原，其次为坝上高原，冀东沿海平原、冀西北间山盆地、中部平原、太行山山地、燕山山前平原、冀北山地和太行山山前平原的释尘通量依次降低，其中，排名前四的区域是河北省土壤风蚀释尘较严重的区域。

3）CLM4.5土壤风蚀模型特定的参数体系空间分辨率较低，普遍为0.05°，对于小区域的模拟研究并不适用。因此，本研究充分利用卫星遥感和地理大数据的优势，采用了空间分辨率为1 km的高分辨率的植被覆盖度、风速、土壤湿度、土地覆被分类、湖泊、河流和湿地数据等本地化参数，提高了模型输入参数的精度，提高了模拟结果的准确性，在小区域土壤风蚀释尘量模拟方面进行了新的尝试。

[参 考 文 献]

[1] 陈璐，陈莉，石梦双. 天津市郊区土壤细粒径空间分布及变异性分析[J]. 环境科学与技术，2013，36(7)：29－34.Chen Lu, Chen Li, Shi Mengshuang. Spatial distribution and variability of soil fine diameter in suburbs of Tianjin [J].Environmental Science & Technology, 2013, 36(7): 29－34.(in Chinese with English abstract)

[2] Waldman J M, Lioy P J, Zelenka M, et al. Winter time measurements of aerosol acidity and trace elements in Wuhan,acity incentral China[J].Atmospheric Environment, 1991,25(1): 113－120.

[3] 陈杰，赵素平，殷代英，等. 沙尘天气过程对中国北方城市空气质量的影响[J]. 中国沙漠，2015，35(2)：423－430.Chen Jie, Zhao Suping, Yin Daiying, et al. Effect of dust process on air quality in cities of Northern China [J]. Journal of Desert Research, 2015, 35(2): 423－430. (in Chinese with English abstract)

[4] Tegen I, Fung I. Modeling of mineral dust in the atmosphere:Sources, transport, and optical thickness [J]. Journal of Geophysical Research, 1994, 99(D11): 22897－22914.

[5] Goudie A S. Desert dust and human health disorders [J].Environment International, 2014, 63: 101－113.

[6] Crooks J L, Cascio W E, Percy M S, et al. The association between dust storms and daily non-accidental mortality in the United States, 1993－2005[J]. Environmental Health Perspectives, 2016, 124(11): 1735－1743.

[7] Losno R, Bergametti G, Carlier P, et al. Majorions in marine rainwater with attention to sources of alkaline and acidic species[J]. Atmospheric Environment, 1991, 25(3－4): 763－770.

[8] Maki T, Puspitasari F, Hara K, et al. Variations in the structure of airborne bacterial communities in a downwind area during an Asian dust (Kosa) event [J]. Science of the Total Environment, 2014, 488: 75－84.

[9] Li H, Qian X, Hu W, et al. Chemical speciation and human health risk of trace metals in urban street dusts from ametropolitan city, Nanjing, SE China [J] Science of the Total Environment, 2013, 456: 212－221.

[10] Ren Y, Ichinose T, He M, et al. Aggravation of ovalbumin-induced murine asthma by co-exposure to desert-dust and organic chemicals: An animal model study[J].Environment Health, 2014, 13: 83.

[11] Choobari O A, Zawar-Reza P, Sturman A. The global distribution of mineral dust and its impacts on the climate system: A review [J]. Atmosphere Research, 2014, 138: 152－165.

[12] Gillette D A, Passi R. Modeling dust emission caused by wind erosion [J]. Journal of Geophysical Research, 1988,93(D11): 14233－14242.

[13] Gillette D A, Admas J, Muhs D, et al. Threshold friction velocitie sand rupture module for crusted desert soils for the input of soil particles into the air [J]. Journal of Geophysical Research, 1982, 87: 9003－9015.

[14] Foroutan H, Young J, Napelenok S, et al. Development and evaluation of a physics-based windblown dust emission scheme implemented in the CMAQ modeling system[J].Journal of Advances in Modeling Earth Systems, 2017, 9(1):585－608.

[15] Jeon W, Choi Y, Percell P, et al. Computationally efficient air quality forecasting tool: Implementation of STOPS v1.5 model into CMAQ v5.0.2 for a prediction of Asian dust[J].Geoscientific Model Development,2016, 9 (10): 3671－3684.

[16] Shao Y P, Raupach M R, Leys J F. A model for predicting Aeolian sand drift and dust entrainment on scales from paddock to region [J]. Soil Research, 1996, 34(3): 309－342.

[17] Lu H, Shao Y P. A new model for dust emission by saltation bombardment [J]. Journal of Geophysical Research, 1999,104(D14): 16827－16842.

[18] Shao Y P. Physics and Modeling of Wind Erosion[M].London: Kluwer Academic Publishers, 2000: 195－197.

[19] Alfaro S C, Gomes L. Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distribution in source areas [J]. Journal of Geophysical Research, 2001, 106(D16): 18075－18084.

[20] Alfaro S C, Rajot J L, Lafon S, et al. Does soil composition have an influence on the sandblasting process [R]. Paris: The 2nd Workshop on Mineral Dust, 2003: 10－12.

[21] Shaw W J, Jerry A K, Fritz B G, et al. An evaluation of the wind erosion module in DUSTRAN [J]. Atmospheric Environment, 2008, 42(8): 1907－1921.

[22] Zender C S, Bian H S, Newman D. Mineral dust entrainment and deposition (DEAD) model: Description and 1990s dust Climatology [J]. Journal of Geophysical Research, 2003,108(D14): 4416.

[23] Joussaume S. Three-dimensional simulations of the atmospheric cycle of desert dust particles using a general circulation model [J]. Journal of Geophysical Research, 1990,95(D2): 1909－1941.

[24] Ginoux P, Chin M, Tegen I, et al. Sources and distribution of dust aerosols simulated with GOCART model[J]. Journal of Geophysical Research, 2001,106(D17): 20225－20273.

[25] 张小曳. 亚洲粉尘的源区分布、释放、输送、沉降与黄土堆积[J]. 第四纪研究，2001，21(1)：29－40.Zhang Xiaoye. Source distributions, emission, transport,deposition of Asian dust and loess accumulation[J].Quaternary Sciences, 2001, 21(1): 29－40. (in Chinese with English abstract)

[26] 宿兴涛，王宏，许丽人，等. 沙尘气溶胶直接气候效应对东亚冬季风影响的模拟研究[J]. 大气科学，2016，40 (3)：551－562.Su Xingtao, Wang Hong, Xu Liren, et al. A numerical simulation study on the impacts of dust aerosol direct climatic effects on the East Asian winter monsoon[J].Chinese Journal of Atmospheric Sciences, 2016, 40 (3): 551－562. (in Chinese with English abstract)

[27] 车慧正，张小曳，李杨，等. DPM模型计算中国北方沙漠地区粉尘释放通量[J]. 干旱区资源与环境，2005，19(5)：49－55.Che Huizheng, Zhang Xiaoye, Li Yang, et al. Calculations of dust emission flux in northern deserts of China by DPM model[J]. Journal of Arid Land Resources and Environment,2005, 19(5): 49－55. (in Chinese with English abstract)

[28] 陈莉，李涛，韩婷婷，等. WEPS模型下天津郊区风蚀尘对城区空气质量的影响[J]. 中国环境科学，2012，32(8)：1353－1360.Chen Li, Li Tao, Hang Tingting, et al. Effect of wind-blown dust derived from Tianjin suburban on air quality of central district based on WEPS model[J]. China Environmental science, 2012, 32(8): 1353－1360. (in Chinese with English abstract)

[29] 元天刚，陈思宇，康丽泰，等. 1961～2010年中国北方沙尘源区沙尘强度时空分布特征及变化趋势[J]. 干旱气象，2016，34 (6)：927－935.Yuan Tiangang, Chen Siyu, Kang Litai, et al. Temporal and spatial distribution characteristics and change trends of dust intensity in dust source regions of northern China during 1961-2010[J]. Journal of Arid Meteorology, 2016, 34 (6):927－935. (in Chinese with English abstract)

[30] 陈思宇，黄建平，李景鑫，等. 塔克拉玛干沙漠和戈壁沙尘起沙、传输和沉降的对比研究[J]. 中国科学：地球科学，2017，47(08)：939－957.Chen Siyu, Huang Jianping, Li Jingxin, et al. Comparison of dust emissions, transport, and deposition between the Taklimakan Desert and Gobi Desert from 2007 to 2011[J].Science China Earth Sciences, 2017, 47(08): 939－957. (in Chinese with English abstract)

[31] 郭俊，银燕，王咏薇，等. 东亚沙尘分布、源汇及输送特征的模拟研究[J]. 中国环境学，2017，37(3)：801－812.Guo Jun, Yin Yan, Wang Yongwei, et al.Numerical study of the dust distribution, source and sink, and transport features over East Asia [J]. China Environmental Science, 2017, 37(3):801－812. (in Chinese with English abstract)

[32] 陈广善，刘晓东，陈葆德. 亚洲大陆2000～2002年春季大气沙尘时空特征的数值模拟[J]. 环境科学，2006，27(1)：1－8.Chen Guangshan, Liu Xiaodong, Chen Baode. Numerical simulations of spatial and temporal characteristics of airborne dust over Asia during springs of 2000 to 2002[J].Environmental Science, 2006, 27(1): 1－8. (in Chinese with English abstract)

[33] 王仁德，常春平，彭帅，等. 基于粒度对比法的坝上农田风蚀与粉尘释放量估算[J]. 农业工程学报，2013，29(21)：108－114.Wang Rende, Chang Chunping, Peng Shuai, et al. Estimation on farmland wind-erosion and dust emission amount in Bashang of Hebei province by grain composition contrast[J].Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(21): 108－114. (in Chinese with English abstract)

[34] 康永德，夏开伟，杨兴华，等. 近地层沙尘水平通量与集沙效率野外观测分析[J]. 农业工程学报，2017，33(13)：168－174.Kang Yongde, Xia Kaiwei, Yang Xinghua, et al. Analysis of horizontal flux and sand collecting efficiency of sand dust near ground surface based on field observations [J].Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(13): 168－174. (in Chinese with English abstract)

[35] Shao Y P, Dong C H. A review on east Asian dust storm climate, modelling and monitoring [J]. Global and Planetary Change, 2006, 52(1－4): 1－22.

[36] Harrison S P, Kohfeld K E, Roelandt C, et al. The role of dust in climate changes today, the last glacial maximum and in the future [J]. Earth-Science Reviews, 2001, 54(1-3): 43－80.

[37] 王卫，冯忠江，陈辉. 河北地理[M]. 北京：北京师范大学出版社，2012：3－9.

[38] Clark I, Assamoi P, Bertrand J, et al. Characterization of potential zones of dust generation at eleven stations in the southern Sahara[J]. Theoretical and Applied Climatology,2004, 77(3): 173－184.

[39] Oleson K W, Lawrence D M. Technical description of version 4.5 of the Community Land Model[R]. NCAR Earth System Laboratory Climate and Global Dynamics Division,USA, 2013.

[40] Webb N P, Strong C L. Soil erodibility dynamics and its representation for wind erosion and dust emission models [J].Aeolian Research, 2011, 3(2): 65－179.

[41] Wu C L, Lin Z H, He J X, et al. Aprocess-oriented evaluation of dust emission parameterizations in CESM: Simulation of a typical severe dust storm in East Asia[J]. Journal of Advances in Modeling Earth Systems, 2016, 8(3): 1432－1452.

[42] Niu G Y, Yang Z L. Effects of frozen soil on snowmelt runoff and soil water storage at a continental scale[J]. Journal of Hydrometeorology, 2006, 7(5): 937－952.

[43] Zeng X B, Shaikh M, Dai Y J, et al. Coupling of the Common Land Model to the NCAR Community Climate Model [J]. Journal of Climate, 2002, 15(14): 1832－1854.

[44] Levis S, Bonan G B, Vertenstein M, et al. The community land model's dynamic global vegetation model(CLM-DGVM): Technical description and user's guide[R].National Center for Atmospheric Research, Colorado, 2004.