Stochastic parameterization of dust emission and application to convective atmospheric conditions

We develop a parameterization scheme of convective dust emission for regional and global atmospheric models. Convective dust emission occurs in the absence of saltation as large eddies intermittently produce strong shear stresses on the surface and entrain dust particles into the air. This dust emission mechanism has not been included in the traditional dust models. The scheme presented in this study is a new approach which takes account of the stochastic nature of convective dust emission. It consists of the statistical representations of soil particle size, inter-particle cohesion, and instantaneous surface shear stress. A method of determining the probability density function of the latter quantity is proposed. Dust emission is then estimated from the overlap of the probability density functions of the aerodynamic lifting and inter-particle cohesive forces. The new scheme is implemented into the WRF/Chem model and applied to dust modeling in the Taklimakan Desert. A comparison with lidar data shows that the model can reproduce the main features of the dust patterns and their diurnal variations. For the case studied, convective dust emission is typically several μg m −2 s−1 and at times up to 50 μg m −2 s−1.


Introduction
The existing dust emission schemes used in regional and global atmospheric models are mainly concerned with the parameterization of dust emission caused by saltation bombardment (e.g.Marticorena and Bergametti, 1995) and aggregates disintegration (e.g.Shao, 2004).These schemes are fairly effective for the simulation of dust emission during strong dust events (e.g.Shao et al., 2010).However, they are not designed for quantifying dust emission under weak wind con-ditions, which is generated by intermittent turbulence rather than by the mean wind shear and the associated saltation of sand-sized grains.As pointed out in previous studies (Marticorena et al., 1997;Shao et al., 1993), the intensity of dust emission due to direct aerodynamic lifting is much weaker than that due to saltation bombardment and aggregates disintegration, and can therefore be neglected in modeling strong dust events (e.g.dust storms).However, while the intensity of aerodynamic dust emission is weak in general, it may occur frequently, e.g. on daily basis in desert areas.In contrast, strong dust events occur much less frequently, maybe a few times a month during the peak dust season.Therefore, aerodynamic dust emission may constitute a major contributor to the regional and global dust budgets on seasonal, annual or longer time scales.For example, Koch and Renno (2005) reported that convective plumes and vortices contribute to about 35 % of the global budget of mineral dust.Several recent review papers also highlighted the potential importance of micro-scale dust emission (Shao et al., 2011;Knippertz and Todd, 2012).
Aerodynamic dust emission by turbulence is most obvious under convective (atmospheric boundary layer) conditions, as exemplified by dust devils in desert in summer.In this study, we call aerodynamic dust emission due to convective turbulence convective dust emission.A few studies on convective dust emission have been carried out in recent years.Ansmann et al. (2008) investigated the vertical structure of convective dust plumes using lidar observations in Morocco during the Saharan Mineral Dust Experiment (SAMUM).Heintzenberg (2008) and Loosmore and Hunt (2000) carried out wind-tunnel experiments and found that dust emission also occurs in the absence of saltation.Gledzer et al. (2009) considered particles significantly smaller than the thickness Published by Copernicus Publications on behalf of the European Geosciences Union.
of the viscous sublayer and estimated the mass concentration of dispersed fine particles in the viscous thermal boundary layer based on field measurements in the desertified areas near the Caspian Sea.The driving parameters in their approach are friction velocity and temperature drop near the surface.Ito et al. (2010) carried out a large-eddy simulation (LES) to estimate convective dust emission.They considered a particle size range from 1 to 10 µm and computed the dust fluxes according to Loosmore and Hunt (2000).Their simulations show that dust concentration in the mixed layer linearly increases with surface heat flux.
The first aim of this study is to develop a scheme for the parameterization of convective dust emission.An important feature of the scheme is a statistical description of the stochastic variables involved in the process.This is a new approach to dust emission modeling in contrast to the conventional dust emission schemes.The theory for the new scheme is described in Sect. 2. The second aim of the study is to develop the capacity of assessing the contribution of convective dust emission to regional and global dust budgets.To this end, a technique is proposed to implement the new scheme in the framework of the WRF/Chem model.An atmospheric model (e.g.regional) has a typical grid size of a few to a few tens of kilometers.For a given model grid cell, the joint probability density function (joint pdf) of the horizontal and vertical velocity components, i.e. the shear stress, is generated based on the pdfs for the individual components following Manomaiphiboon and Russell (2003).The WRF/Chem model with the new scheme is then implemented for the simulation of a convective dust event in the Taklimakan Desert.The model results are compared with the lidar data collected at Aksu.

Convective dust emission
The main mechanisms for dust emission are aerodynamic entrainment, saltation bombardment, and aggregates disintegration (Shao, 2008).In case of strong winds, saltation bombardment and aggregates disintegration are the dominant mechanisms in comparison to aerodynamic entrainment.In the context of dust emission, strong winds refer to the situations when u * > u * t , where u * is friction velocity and u * t is threshold friction velocity for saltation.If wind is too weak to activate saltation (u * < u * t ), then aerodynamic entrainment becomes the prevalent mechanism for dust emission.Convective dust emission is the most important form of aerodynamic entrainment and is thus the focus of the new dust scheme.The applicability of our approach is, however, not limited to convective conditions.It can be easily extended to other turbulent conditions, e.g.shear generated turbulence, by modifying the momentum transport parameterization described in Sect.4.1.
Convective turbulence occurs in regions of strong surface heating, such as desert areas in summer.A super-adiabatic temperature gradient commonly exists in the atmosphere near the surface, which forces the development of thermals and generates strong vertical velocity fluctuations.
Let us consider a unit area covered with dust particles of different sizes.A force τ is exerted by turbulence on the unit area ([τ ] = N m −2 = Pa).If the force is evenly distributed, then the force exerted on a particle with a cross-section of a is The fraction η i of particles of size d i in a given size interval δd i is Here, p (d i ) denotes the particle size distribution function (psd).Therefore, the force exerted on all particles in this size interval is given by The quantity τ is the instantaneous vertical flux of horizontal momentum given by with air density ρ.Note that the current values of u w and v w are used instead of the mean values u w and v w of the Reynolds shear stress.As detailed in Sect.3, τ is parameterized by means of a joint pdf of the horizontal and vertical wind components.Dust emission can be expressed as the number flux n d of dust particles of size d multiplied by the particle mass m p : where N d is the particle number concentration and w p the particle vertical velocity which obeys the equation of particle motion, namely, The first term contains the particle response time T p and the vertical component (w p − w a ) of the particle-to-air relative velocity U r , with w a being the vertical velocity of the air.In general, T p can be expressed as with aerodynamic drag coefficient C d , and U r being the magnitude of U r (Shao, 2008).The particle density ρ p is approximately 2560 kg m −3 .Term I in Eq. ( 6) describes the behavior of a particle in air and is important as soon as the particle is lifted from the surface.Term II reflects the particle acceleration due to gravity.Term III is the most important term concerning the emission process.The force exerted on the particle by wind as described in Eq. ( 3) is f and the cohesive force F i .F i is only active up to a height of the order of particle diameter and is then zero (Fig. 1a).In our scheme, both f and F i are stochastic quantities which obey certain probability distributions.Therefore, dust emission is proportional to the overlap of the two distributions (Fig. 1b).
A dust particle can be considered to be emitted from the surface if it passes through the viscous sublayer adjacent to the surface.Therefore, Eq. ( 6) is integrated over the depth of the laminar layer to determine the vertical velocity of the dust particle motion: with particle terminal velocity w t = gT p .The effective turbulent flux τ (sum of molecular and turbulent fluxes) is approximately constant with height in the surface layer (Stull, 1988).In the viscous sublayer, τ obeys the Newtonian law: with ν being the kinematic viscosity.Suppose U (z = 0) = 0.
Then, an integration of Eq. ( 9) yields Suppose the instantaneous friction velocity is u * = √ τ/ρ and the transition between the laminar and turbulent flows (which defines the top of the viscous sublayer) starts at a friction Reynolds number of u * z/ν = 5 (Schlichting et al., 2003).Then, the thickness δ of the viscous sublayer can be estimated as The first term in Eq. ( 8) represents dust deposition and therefore does not have to be included in the dust emission scheme.The particle vertical velocity w p in Eq. ( 5) can now be substituted by Eq. ( 8).As illustrated in Fig. 2, the particle number concentration N d must be inversely proportional to the depth of the laminar layer, namely, N d = α N /δ with α N in m −2 being the proportionality parameter.It follows that the dust emission flux can be expressed as The parameter α N is an unknown empirical parameter to be determined by comparison of the scheme with observations.Finally, the total convective dust emission for all particles of a given particle size interval is given by The innermost integration accounts for the stochastic behavior of the cohesive force F i (see Sect. 3.1) and the integration over the instantaneous shear stress τ accounts for the stochastic behavior of the lifting force f , which is related to τ by Eq. ( 3).Thus, the two integrals multiplied by the number of particles of size d i describe the amount of dust emitted for this particle size (see also Fig. 1b).Finally, the integration over particle diameter d i of F (d i ) multiplied by the psd yields the total dust emission for all particle sizes.The principal mechanisms of convective dust emission as parameterized by Eq. ( 13) are summarized in Fig. 2.

Parameterization
As shown in the previous section, the parameterizations for both the cohesive and lifting forces are necessary, and the parent-soil psd must be specified as an input quantity.As our main concern is convective dust emission under weak mean wind conditions, minimally-disturbed psds are used.These The distributions are calculated using Eq. ( 14).
are approximations to the parent soil psd seen by turbulence, as they are obtained with the minimal mechanical and chemical disturbances to the soil samples.In contrast, fullydisturbed psds are obtained by applying strong mechanical and chemical forces to disaggregate the soil particles, which rarely occurs in real dust emission processes (Shao, 2008).The use of minimally-disturbed psd is thus more appropriate for our scheme.In this study, the minimally-disturbed psds for the parent soils are approximated as the sum of four lognormal distributions as described by Shao (2001).The U.S. Department of Agriculture's (USDA's) soil classification distinguishes 12 soil texture classes based on the percentages of sand, silt, and clay.Due to the lack of psd measurements, these classes are regrouped into four classes in the model, namely, sand, loam, sandy clay loam, and clay.For these soils, the psds are shown in Fig. 3a.

Cohesive force
The particle retarding force includes the cohesive force and the gravity force.The cohesive force is mainly composed of Van der Waals-forces, electrostatic forces, capillary forces, and chemical binding forces (Shao, 2008).For small particles, cohesive force dominates, while for large particles gravity force dominates.Factors such as particle shape, mineral composition, surface roughness, etc., profoundly affect the inter-particle cohesion and consequently, the cohesive force may differ over orders of magnitude for particles in the same size range.In general, the factors contributing to the interparticle cohesion are so various that the cohesive force can be best treated as a stochastic variable with certain probability distributions.Based on the data of Zimon (1982), the retarding force appears to obey a log-normal distribution.In the scheme developed here, the pdf of F i is given by where the mean value Fi and the geometric standard deviation σ F i ( Fi and σ F i in mdyn = 10 −8 N, and d in µm) are The coefficients herein are obtained by fitting the pdf to the data of Zimon (1982).The results for six particle sizes between 1 and 20 µm are shown in Fig. 3b.As seen, Fi increases with particle diameter and the range of F i variation increases with decreasing particle size due to the greater dominance of the stochastic cohesive force and the reduced importance of gravity force.We emphasize, however, the above described parameterization of F i is only provisional.More data is required for improved treatment of F i and tests on the model sensitivity to this treatment is necessary.Nevertheless, it is sufficient to use the data to illustrate our idea of stochastic dust modeling.

Lifting force
To parameterize the shear stress generated by convective turbulence, the joint pdfs of (u , w ) and (v , w ) are required.These are determined by use of the similarity theory.Since the velocity fluctuations u and v behave similarly, they have equal variances which can be combined to σ 2 u h = √ 2σ 2 u .In the remainder of this paper, u denotes the total horizontal wind component.The (u , w ) joint pdf can be constructed on the basis of the pdfs of u and w .The required statistical moments, such as variances and skewnesses can be estimated from the mixed layer similarity theory.The appropriate scaling velocity and length are respectively the convective velocity scale w * and planetary boundary layer (PBL) depth z i .According to the similarity laws (Kaimal and Finnigan, 1994), the variances of u and w are given by while the skewness for w , γ = w 3 /σ 3 w , is defined by The mean value of the momentum transfer, τ R = −ρu w , is also obtained from the similarity theory, i.e.

−u w w
where u w represents the covariance σ 2 uw between u and w .L is the Obukhov length.
Following Manomaiphiboon and Russell (2003), the pdf of τ , and thus that of f , is obtained in three steps.First, the pdfs of u and w , p u and p w , and the corresponding cumulative distribution functions (cdfs), P u and P w , are computed.Second, the (u , w ) joint pdf, p j , is computed using the above pdfs and cdfs.Last, the pdf of τ , p τ , is obtained from p j with τ = −ρu w .Note that since different combinations of u and w can yield the same τ , p τ is the sum of p j for the same u w product as follows To determine p u , P u , etc., the method proposed by Manomaiphiboon and Russell (2003) is used.While p u is Gaussian, p w is approximated with a Bi-Gaussian pdf as convective turbulence exhibits a negative skewness.Manomaiphiboon and Russell (2003) parameterized the joint pdf of u and w according to the technique of Koehler and Symanowski (1995).The essence of this technique is to derive the joint pdf from the predefined marginal distributions.The shapes of the marginal distributions are conserved during the transformation to a joint pdf.
Figure 4 shows the process of determining p τ step by step.The mean value of τ , τ R , as well as the variances of u and w , estimated by means of the mixed layer similarity theory, are herein used to match the parameterized pdf of τ to the environmental conditions given by the model simulations.The above described dust emission scheme and the parameterizations have been integrated in the WRF model (Advanced Research version) (Wang et al., 2009) with chemistry (WRF/Chem, Grell et al., 2005).The GOCART (Georgia Tech/Goddard Global Ozone Chemistry Aerosol Radiation Transport) aerosol scheme (Chin et al., 2000) is used as the basis for the implementation of our scheme.The convective dust emission module, together with the other dust emission modules, can be chosen in combination with the GOCART simple chemistry option (Kang et al., 2011).We refer to the above model as WRF/Chem Dust.

Offline tests
To examine the performance of the new dust scheme, offline tests are first carried out.In these tests, α N is set to 1.The scheme sensitivity to w * and soil attributes is first examined.For this purpose, p τ is computed for various w * values in the range between 0.5 and 4 m s −1 .The PBL depth z i and the Obukhov length L are set to z i = 1000 m and L = −10 m.The results for p |τ | are shown in Fig. 5a.
Analysis shows that τ R is almost constant and always positive, indicating that the net momentum flux is directed downwards.In contrast to the small variation of τ R , the standard deviation of τ increases with w * corresponding to the rising intensity of turbulence, i.e. the stronger variations in u and w also result in stronger variations in τ .Due to the greater variance of τ , the overlap of the pdfs of τ and F i also increases (Fig. 1b).As a result, dust emission increases with w * .The small variation of τ R is understandable, because the shear stress due to the mean wind is not included.Clearly, the inclusion of u * would lead to increased variations in τ R .The characteristics of p τ for different wind and stability conditions as well as the implications to dust emission will be described in a future study.
Tests are performed to investigate the dependency of dust emission F on w * for different soil types (Fig. 5b).As α N is not yet known, F /α N is shown (a preliminary estimate of α N is given in Sect.4.2.1).For constant w * , F is weak for soils with large particles (e.g.sand), and clearly increases for soils rich in small particles (e.g.clay).Figure 5b also shows that F substantially increases with w * .

Case study
The WRF/Chem Dust model is implemented to the Taklimakan Desert to the simulation of a weakly convective dust event.The numerical results are then compared with the measurements of a ground-based lidar at the Aksu Water Balance Experimental Station, Xinjiang Institute of Ecology and Geography of the Chinese Academy of Sciences (Jin et al., 2010).The station is located in the northern part of the Taklimakan Desert (40.62 • N, 80.83 • E, 1028 m above sea level).The measurements were taken with a Mie-scattering polarization lidar, which continuously determines the vertical distribution of aerosols from the PBL through the troposphere up to the stratosphere (Kai et al., 2008).The three-day period, 23-25 March 2009, is chosen for comparison with the simulations.During these days, the noon time (14 LST) surface heat flux fell between 100 and 250 W m −2 in much of the simulation domain (not shown).At the fringes of the Taklimakan desert, it exceeded 300 W m −2 on occasions.The ratio of z i /L was mostly less than −45.Thus, the case studied is convective, and satisfies the requirement for testing the scheme.Although convective turbulence is more prevalent in summer, we have at this stage no other suitable data for model verification.
The model run is set up with a horizontal resolution of 25 km.The domain extends over 1500 km × 750 km, which corresponds to 60 grid points in x-direction and 30 in ydirection.The domain as well as the topographical height can be seen in Fig. 7.In the vertical direction, 28 model levels are used up to a pressure level of 50 hPa.The Yonsei University (YSU) PBL scheme (Hong et al., 2006) was applied to estimate w * and z i .The source areas for dust emission calculation have been defined as suggested by Shao and Dong (2006).The authors calculated the dust concentration on the basis of synoptic visibility reports using an empirical relationship.A location is classified as a potential dust source area if the average dust concentration exceeds a threshold value and additional criteria regarding erodibility and vegetation cover are satisfied.The potential dust source area is shown in Fig. 7 (background,dotted).
The geographical data are interpolated from terrestrial data based on the default 24-category land use classification and 16-category soil classification in WRF with a 10 m resolution.The vegetation cover data used in this study is combined from vegetation type data of the State Key Laboratory of Resources and Environment Information System (LREIS) of the Institute of Geography, Chinese Academy of Sciences, and NDVI (Normalized Differenced Vegetation Index) data derived from NOAA/NASA (National Oceanic and Atmospheric Administration/National Aeronautics and Space Administration) Pathfinder AVHRR (Advanced Very High Resolution Radiometer) land dataset (Shao and Dong, 2006).The meteorological initial and lateral boundary conditions are specified by the 6-hourly Final Analysis data (FNL) of the NCEP Global Forecast System (GFS) with 1 • resolution.Four sequential model runs have been made to enable a full update of the meteorological conditions every 24 h.Thus, the simulation covers 22 March 2009 00:00 UTC to 26 March 2009 00:00 UTC, which corresponds to 22 March 05:30 LST to 26 March 05:30 LST, including one day spin up time before the period of comparison.The resulting dust concentration of each one-day simulation is passed to the consecutive simulation as initial condition.Four particle size bins are currently used in the model: d ≤ 2.5, 2.5 < d ≤ 5, 5 < d ≤ 10, and 10 < d ≤ 20 µm respectively for bins 1-4.

Comparison to lidar measurements
A lidar measures the backscattering of aerosols, water droplets, and other scattering objects in the atmospheric column through which the lidar beam passes.The backscattering coefficient of aerosols, β a , can be used to calculate quantities such as backscattering ratio R, which is the ratio of total (molecular plus aerosol) to molecular backscattering, with β m being the backscattering coefficient of molecules (Kai et al., 2008).R is measured in 60 m height intervals and the lidar overlap effect is empirically corrected before the lidar data are used for the analysis.Dust concentration, c, can be approximated based on the backscattering ratio, R, according to where a = 0.04 mg m −3 is a coefficient determined by fitting R to c derived from the near-surface dust concentration (observed using high volume sampler) and a prespecified dust concentration profile.
The column dust load in turn can be derived by integrating the dust concentration in the vertical direction.We assume that all aerosols below the cloud base are dust particles and the column dust load is determined by emission, advection, and deposition.In reality, not all aerosols are dust particles and a background aerosol concentration is present in the atmosphere, which has to be removed from the measurements for comparison with the model simulation.Since no initial conditions of aerosol concentration for the dust-free situations are available, a mean profile is calculated and removed from the concentration data.
We first compared the model-simulated and lidar-observed dust load below the lowest cloud base (about 4300 m for the study period).This comparison turned out to be less meaningful, because for reasons yet to be clarified, the highest dust concentration occurred in heights above the boundary layer.The dust there was unlikely to be related to local convective dust emission.For this reason, we used the model-simulated boundary layer height (pblh) as the reference level and calculated the PBL dust load by integrating the dust concentration up to pblh.To find the most appropriate comparison, we tested several options by comparing the lidar data with the model data (I) from the Aksu grid cell, (II) averaged over 9 grid cells surrounding Aksu, and (III) averaged over 25 grid cells surrounding Aksu.The α N parameter is calculated for each of these options by fitting the model data to the lidar data for the 12-h period between 12:00 and 24:00 LST, 24 March 2009, by using a MATLAB ® robust curve fitting technique with a Cauchy weighting function.Figure 6a, c, and e show the scatter plots of the model versus the lidar PBL dust load, together with the linear regressions.Colors indicate the time of the data points.From the slope of the straight lines, we found α N = 1912.9,785.2, and 685.1 m −2 for comparison options I, II and III, respectively.Figures 6b, d, and e show the PBL dust load for the days of 23-25 March 2009 estimated from the lidar data and the model simulations using the calibrated α N values.Note that the lidar PBL dust load also varies slightly as the pblhs estimated for options I, II, and III are somewhat different.
Figure 6a and b present the results for option I.The diurnal cycle of the dust loading is reproduced, but some problems exist.The model overestimated the PBL dust load for the mornings of 24 and 25 March and the evening of 25 March.Figure 6c and d present the results for option II and III.The degree of agreement between the lidar and model data is similar for both options, except for 23 March.Again, the model over predicted the PBL dust load for the evening of 25 March.Only in option III (Fig. 6f), the model shows a decreasing PBL dust load at this time.No substantial differences between options II and III can be seen for other times.The calibrated value of α N decreased from 1912.9 to 685.1 m −2 from option I to III.By comparing the model results for the Aksu grid cell and the adjacent grid cells, it is found that the model-simulated dust emission for the Aksu grid cell is smaller.To reduce the model uncertainties, it seems reasonable to accept the α N values obtained from option II and III.The difference between the α N values is relatively small and therefore, α N = 785.2m −2 is used for subsequent model runs.The coefficient of determination r 2 (see Fig. 6) underlines this choice as option II also has the highest value of r 2 = 0.90.

Quantitative analysis
The above described case is rerun with α N set to 785.2 m −2 .The model simulation shows that dust emission in the Taklimakan Desert is primarily limited to the desert fringes.
As example, Fig. 7 shows the predicted dust emission for 14:00 LST on 23, 24 and 25 March 2009, when convective turbulence is expected to be the strongest.The pattern of dust emission clearly shows its dependency on soil type.In the interior of the Taklimakan Desert, where sand is the dominant soil type, there is little dust emission.For areas where dust emission occurred, the total (all particle size) dust emission fell between 1 and 30 µg m −2 s −1 at noon time, reaching occasionally a maximum of 50 µg m −2 s −1 .The dust concentration in the lowest model layer is up to 150 µg m −3 in the areas of dust emission, reaching on occasions a maximum of 300 µg m −3 .
To study the influence of soil composition on dust emission in greater detail, four locations are selected from the domain representing the four soil groups used for the simulation as follows: The time series of dust emission, column dust load, and w * are analyzed for the four locations.Figure 8 shows the time series of dust emission and dust load for all particle sizes.Clay consists a high proportion of small particles and reveals a clear connection between dust emission and convection, which is most prevalent during noon (Fig. 8a).Sandy clay loam also contains a high fraction of dust and an obvious relation between dust emission and w * (for this location, w * is relatively small, not shown).In contrast, loam contains only a small fraction of dust (Fig. 3a) and hence gives less dust emission.The relatively large dust load (Fig. 8b) found at the location is due to advection rather than local emission, as revealed by the fact that the high dust loads occurred during night.Sand is predominantly composed of large particles and hence produces little convective dust emission.

Dust budget
A dust budget is calculated to estimate the total convective dust emission over the study domain.The dust budget equation, integrated over the 3-D study domain, can be written as where δD is the change with time of the dust load over the domain, AF the total dust emission, AF D the total dust deposition and AF A u and AF A v the transport through the lateral boundaries.Figure 9 shows their time series.Up to about 1500 kg s −1 of dust are emitted at noon time each day by convective turbulence.This makes a total emission of 86.5 kt for the three-day period in the study domain.Deposition reaches about 500 kg s −1 during the time of maximum dust concentration and the total amount of deposited dust for the study period is 70.3 kt.Advection is relatively small during the study period due to the weak winds.δD reaches up to 1000 kg s −1 at noon time.The accumulated dust load over the study period is 7.3 kt.

Model uncertainties
In conventional dust emission schemes, the threshold friction velocity, u * t , is a key parameter.An ideal u * t is defined, which depends only on particle diameter, and it is then corrected to account for the influences of the environmental factors, such as surface roughness, soil moisture, salt crust, etc., by multiplying the ideal u * t with correction functions (Shao, 2008).In our scheme, the concept of threshold friction velocity is not used, and the influences of the environmental factors are reflected in the probability distributions of F i and |τ |.As pointed out in Sect.3.1, our estimate of p (F i ) is provisional.The parameters used to estimate p (F i ) are likely to have considerable uncertainties and vary from case to case depending on the prevailing environmental conditions.To investigate the uncertainties arising from the p (F i ), we conducted a set of sensitivity experiments in which we (a) increased and (b) reduced the mean value Fi by 20 %, and (c) increased and (d) reduced the standard deviation σ F i by 20 %.The perturbed p (F i ) functions for the sensitivity experiments are shown in Fig. 10a.We repeated the model simulations for 24 March with the new p (F i ) functions and compared the herewith obtained dust emissions with respect to the reference run. Figure 10b shows the domain integrated emission for 24 March 2009 for the sensitivity experiments together with the reference dust emission ±10 %.As can be seen, a 20 % increase/decrease in Fi leads to about 5 % decrease/increase in dust emission.In comparison, the changes in σ F i only affect the predicted dust emission to a small degree.This is because the changes in p (F i ) are rather small.Spatially, the changes in Fi resulted in about 5 % changes in dust emission in most regions where dust emission occurred, whereas the changes in σ F i produced significant changes only in the regions of strong emission (up to 5 %, not shown).

Concluding remarks
In this paper, a parameterization scheme for convective dust emission is presented.We have pointed out that convective dust emission may play an important role in the global dust budget and the mechanisms for convective dust emission differs profoundly from that for dust emission generated by the saltation of sand-sized grains driven by the mean wind.The construction of the new scheme is based on two important observations: (1) convective eddies generate intermittently large shear stresses on fractions of the aeolian surface; and (2) due to the stochastic nature of inter-particle cohesion, there always exists a fraction of free dust that can be entrained into the air by turbulence or weak winds without saltation.The fundamental difference between our scheme and the conventional dust emission schemes is that in our scheme the stochastic nature of dust emission has been taken into consideration, in terms of the probability distributions of the inter-particle cohesive forces and the turbulent shear stress.
We have developed the WRF/Chem Dust model by integrating the new dust emission scheme, together with several other conventional schemes (Kang et al., 2011), into the WRF/Chem model, which is then applied to the simulation of a convective dust event in the Taklimakan Desert.The model results are compared with the lidar data obtained at Aksu.The model is found to be able to reproduce the basic spatial and temporal features of dust patterns in the atmospheric boundary layer.We have used the lidar data for a 12-h period to calibrate the model parameter α N and found α N to be around 785.2 m −2 .Based on this choice of α N , convective dust emission is found to be of the order of magnitude 1 to 10 µg m −2 s −1 up to a maximum of 50 µg m −2 s −1 at noon time.During the three-day study period, a total of 86.5 kt dust are emitted from, and 70.3 kt are deposition to, the study area (1500 × 750 km 2 ), resulting in a net dust emission of 7.3 kt.
The model-lidar comparison must be viewed in light of several limitations.The study area and the selected dates   (23-25 March 2009) for the model validation are by no means optimal.In the Taklimakan Desert, convective turbulence during this time of the year is not as pronounced as in summer.Also, the presence of clouds during the study period reduced the incoming radiation and prevented surface heating.The topography around Aksu is rather complicated, as it is located close to the Tianshan Mountains.The flow there can be strongly influenced by the mountain-valley winds, especially in the morning and evening.The flow in the Taklimakan Desert, surrounded by the Tibetan Plateau, the Pamir Plateau and the Tianshan Mountains, is also very complex.Kim et al. (2009) studied the dust layer height at Aksu in April 2002 and found its diurnal variation is strongly influenced by the local circulation which can lead to increased dust load in the morning and night.This phenomenon is also embedded in the lidar measurements used in this study.Aksu is therefore not an ideal site for the calibration of our model.The Sahara would be a more suitable reference for testing our model, but we have so far no available data for the area.
A more reliable calibration of the α N parameter is also desirable.The lidar data used in this study have considerable uncertainties in determining dust load in the atmospheric boundary layer due to the overlapping effect.Further, as indicated in Sect.3.1, the parameterization of the cohesive force is only provisional.
Given the above restrictions, we consider the model and the lidar data to be in reasonable agreement.With the calibrated α N value, typical convective dust emission is found to be about 3 µg m −2 s −1 for clay and 1 µg m −2 s −1 for sandy clay loam for w * = 1 m s −1 .The typical magnitude of convective dust emission is consistent with that of the dust emission observed under weak wind conditions (u * ≤u * t ), as reported in the literature (Fig. 7.2 of Shao, 2008).Nickling and Gillies (1993) measured the near-surface vertical dust flux for non-dust storm periods in Mali and reported values around 1 µg m −2 s −1 for u * <0.2 m s −1 .Dust flux measurements made under weak wind conditions for other regions show the same order of magnitude (Nickling et al., 1999).
Several improvements to the proposed scheme are planned.In this study, we concentrated on dust emission by convective turbulence and assumed turbulence is buoyancy driven.We have therefore used w * as the scaling velocity for the variances of turbulent velocity according to the mixedlayer similarity theory.In fact, turbulence can be both sheardriven and buoyancy-driven.Moeng and Sullivan (1994) introduced a more general velocity scale w m for turbulence, which combines buoyancy and shear production of turbulence, For the generalization of the proposed scheme, such a combination of u * and w * would provide a more adequate parameterization for the pdf of the instantaneous shear stress.
In conventional dust emission schemes, threshold friction velocity u * t is used, which is first calculated for ideal conditions (dry and bare soil) and then corrected to account for the effects of surface roughness, soil moisture, salt content, crust, etc. (Shao, 2008).In this study, u * t is not used, but the impact of the above-mentioned environmental factors must be reflected in the probability distributions of the cohesive forces.The stochastic nature of the cohesive force and its statistical quantification will be the most formidable, yet unavoidable, problem for the application of the new scheme.

Fig. 1 .
Fig. 1.(a) Illustration of the lifting force affecting a particle in the viscous sublayer; (b) illustration of the probabilistic distributions used for the description of the cohesive and lifting forces.

Fig. 2 . 19 Fig. 2 .
Fig. 2. Illustration of a large eddy above the viscous sublayer (flow indicated by arrows), which exerts a on the particles, reduces the depth of the viscous sublayer δ (dotted line) and increases the particle nu concentration.

Fig. 3 .Fig. 3 .
Fig. 3. (a) Parameterizations of particle size distribution composed from four log-normal distributions to fit experimental data.The results for the four soil types used in this study are shown; (b) probabilistic distributionof the cohesive force for particle diameters of 2, 3, 5, 10, and 20 µm plotted in F i • p(F i ).The distributions are calculated using Eq.(14).

Fig. 4 .Fig. 4 .
Fig. 4. Procedure of determining the pdf of τ , p τ , by using the pdfs of u and w , p u and p w (a), their cdfs P u and P w (b), as well as their joint pdf (c).Shown in (d) is the pdf of |τ |, p |τ | .

Fig. 5 .
Fig. 5. (a) Sensitivity of p |τ | to w * ; (b) Sensitivity of dust emission flux F to w * for different soil types.The insert shows the F and w * relationship for w * ≈ 1 m s −1 .

Fig. 6 .
Fig. 6.(a) Scatter plot of model versus lidar PBL dust load for option I (one grid cell) for the 12-h period (12:00 to 24:00 LST, 24 March 2009) for determination of α N as the slope of the linear relationship, with r 2 denoting the coefficient of determination.(b) PBL dust load derived from lidar data (black) and simulations (red) for 23, 24, and 25 March 2009.(c) as (a) and (d) as (b), but for option II (9 grid-cell average); (e) as (a) and (f) as (b) but for option III (25 grid-cell average).Lidar data from K. Kai and Y. Jin, with acknowledgment.

Fig. 8 .Fig. 8 .
Fig. 8. Time series of (a) dust emission and (b) dust load of all particle sizes for the days 23, 24, an 2009 for the four locations representing clay, sandy clay loam, clay, and sand.

Fig. 10
Fig. 10.a) Distribution of F i × p(F i ) for Fi ± 20 % and σ F i ± 20 % with respect to the reference values.Domain integrated dust emission [kg s −1 ] for the sensitivity experiments shown in a). 27

Fig. 10 .
Fig. 10.(a) Distribution of F i × p (F i ) for Fi ± 20 % and σ F i ± 20 % with respect to the reference values.(b) Domain integrated dust emission [kg s −1 ] for the sensitivity experiments shown in (a).