Toward quantitative prediction of dust storms: an integrated wind erosion modelling system and its applications

Environmental Modelling & Software 16 (2001) 233–249 www.elsevier.com/locate/envsoft

Toward quantitative prediction of dust storms: an integrated wind erosion modelling system and its applications Hua Lu a

a,*

, Yaping Shao

b

CSIRO Land and Water, Canberra Lab, GPO Box 1666, Canberra, ACT 2601, Australia b School of Mathematics, University of New South Wales, Sydney, Australia

Received 7 July 2000; received in revised form 4 September 2000; accepted 25 October 2000

Abstract In this paper, we present an integrated wind-erosion modelling system which couples a physically based wind-erosion scheme, a high-resolution atmospheric model and a dust-transport model with a geographic information database. This system can be used to determine the pattern and intensity of wind erosion, in particular, dust emission from the surface and dust concentration in the atmosphere. The system can also be used for the prediction of individual dust-storm events. We have implemented the system for examining the dust storms over the Australian continent in February 1996. It is shown that over the 1 month period, the total dust emission from the continent was about 1.87 Mt. The dominant dust particles are found in the size range from 0 to 11 µm. Larger particles are only found during dust-storm periods. The major dust emission locations and dust pathways detected by the model are comparable with satellite image and climatology of wind erosion in Australia. A discussion is given on the limitations and uncertainties of the system.  2001 Elsevier Science Ltd. All rights reserved. Keywords: Wind erosion; Dust emission; Dust transport; GIS database; Regional scale

1. Introduction The importance of mineral dust in the atmosphere and the associated impacts have long been recognised (Andreae, 1995; Duce, 1995; Junge, 1979). A prerequisite for studying the impact of aeolian dust on atmospheric radiation and circulation is to determine adequately the spatial distribution of dust concentration and its temporal variations. However, despite numerous studies on the subject (Genthon, 1992; Joussaume, 1990; Tegen and Fung, 1994; Westphal et al., 1988), quantitative estimates of dust-emission rate have not been possible until recently. This is because the processes governing the emission and transport of dust particles are very complex and difficult to model. As shown in Fig. 1, dust emission is determined by an interacting set of physical processes governed by many factors, namely, climate (high wind and low rainfall), soil state (composition, tex-

* Corresponding author. Tel.: +61-2-6246-5923; fax: +61-2-62465965. E-mail address: [email protected] (H. Lu).

ture, aggregation and crusting), and surface roughness (non-erodible elements and vegetation). Since wind erosion is sensitive to a range of environmental factors, erosion events are often spatially variable and temporally highly intermittent. The transport of dust involves particle–turbulence interactions in the atmospheric boundary layer. At large scale, it is often driven by intensive meso-scale to synoptic-scale systems, such as squall lines and cold fronts. This paper constitutes a mathematical and computational modelling effort for the quantitative assessment and prediction of dust-storm events, including the region and the intensity of dust emission and the transport of dust in the atmosphere. For this purpose, we have developed an integrated wind-erosion modelling system (IWEMS), which couples a regional weather-prediction model, a dust-emission model and a dust-transport model. The system links these dynamic models with a geographic information database, which provides the necessary input parameters for the system. We have applied the system to the simulation of the dust-storm events over the Australian continent in the period of February 1996. In this paper, we first describe the system

1364-8152/01/$ - see front matter  2001 Elsevier Science Ltd. All rights reserved. PII: S 1 3 6 4 - 8 1 5 2 ( 0 0 ) 0 0 0 8 3 - 9

234

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

Fig. 1.

An illustration of the physical processes which influence dust emission, transport and deposition.

and then discuss the modelling results in the light of field observations. We shall also discuss the uncertainties of the system.

2. System description Fig. 2 shows the structure and computational procedure of IWEMS which is a modularised three-dimensional system with various options for representing wind erosion and dust transport. The high-resolution weatherprediction model, that drives IWEMS, is a primitive equation model formulated in the s-coordinate system with a Lambert conformal projection. It uses a flexible horizontal resolution with multi-layer structure and can be self-nested. In this study, we have used horizontal resolutions between 5 and 75 km and vertical layers of up to 31 vertical layers. Details of this model can be found in Leslie and Purser (1991) and Shao and Leslie (1997). The atmospheric model provides atmospheric forcing data to both the wind-erosion and the dust-transport models. For each physical time step, relevant atmospheric data (e.g. wind velocity and eddy diffusivity) and land-surface data (e.g. soil moisture) are passed to the wind-erosion scheme through an interface, for the computation of dust-emission rate. This emission rate, together with atmospheric variables, is then passed to the dust-transport model for the calculation of instantaneous

grid-mean dust concentration. This calculation is looped over a number of particle-size classes. The land-surface parameters required for the wind-erosion scheme are prepared by a GIS preprocessor. In the current version of IWEMS, dust is considered as a passive scalar. No chemical reactions are considered in the model. 2.1. The wind-erosion scheme The wind-erosion scheme is an improved version of that developed by Shao et al. (1996) and Shao and Leslie (1997). Its structure is as shown in Fig. 3. It comprises three key parameterisations representing: (1) the erosion threshold friction velocity u∗t; (2) the streamwise saltation flux Q; and (3) the dust-emission rate F(i) for N size classes of dust particles. The treatments of these processes are mainly based on the studies of Raupach (1991) on drag partitioning, Owen (1964) on saltation and Lu and Shao (1999) on dust emission. The effects of surface non-erodible elements are further taken into account by multiplying Q by two factors: the ratio of erodible-to-total surface due to vegetation cover Ev, and the surface fraction of the erodible part of the uncovered soil surface due to the presence of gravels, pebbles and rocks Es. No sediment movement is allowed in the areas occupied by non-erodible elements. The grid values of Ev and Es are estimated from the GIS database. The saltation scheme has been validated against field obser-

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

Fig. 2.

Fig. 3.

235

Components and computational procedure of IWEMS.

The input and output structure of the wind-erosion scheme.

vations by Shao et al. (1996). The dust-emission scheme has been tested by Lu and Shao (1999), against the windtunnel data of Rice et al. (1996) and the field observations of Gillette (1977).

The wind-erosion scheme receives weather data from the atmospheric prediction model and soil and vegetation parameters from the GIS database. The main outputs of the wind-erosion scheme are the threshold velocity u∗t

236

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

(ms⫺1), horizontal sand flux Q (g m⫺1 s⫺1), and vertical dust flux F (g m⫺2 s⫺1). F then becomes the input to the dust-transport model. 2.1.1. Threshold friction velocity Based on sound physical reasoning, Shao and Lu (2000) found that the threshold friction velocity for uniform particles over bare, dry and loose soil surfaces can be calculated using a simple expression

冪a 冉s gd+rd冊 a2

u∗t0⫽

1

(1)

p

where d is the particle diameter, g is gravity acceleration, r is the air density and sp is the particle-to-air density ratio. The coefficients, a1=0.0123 and a2=3×10−4 kg s⫺2, are determined by fitting the above expression to several wind-tunnel data sets. For a natural soil surface, u∗t not only depends on particle size but also other factors, such as surface roughness elements, soil moisture and soil aggregation. These factors are considered to modify u∗t0 in the following way u∗t(d)⫽u∗t0(d)RHM where R, H, M are the enlargement functions describing the influence of surface roughness elements, soil moisture and soil surface aggregation and crusting, respectively, as proposed by Shao et al. (1996). Raupach (1991) presented an analytical drag partition scheme as R−1(l)⫽

冋

1 (1−srl)(1+mbrl)

册

2.1.2. Horizontal sand flux Owen’s (1964) theoretical equation for transport-limited saltation flux is used. For uniform sand particle soil, the horizontal sediment flux is calculated as ˜⫽ Q

再

(cru3∗/g)[1−(u∗t(d)/u∗)2] u∗ⱖu∗t(d) u∗⬍u∗t(d).

0

where c is Owen’s coefficient which can be calculated as 0.25+0.33wt(d)/u∗ and is of the order of 1. The relative contribution to the total flux of each size range is assumed to be proportional to its weight fraction in the soil particle size distribution. The total horizontal sediment flux is then evaluated as a weighted integral ˜ (d) over each size class defined by the particle size of Q distribution p(d)

冕

˜ (d)p(d)dd Q⫽ Q

(2)

2.1.3. Vertical dust flux For each particle-size class, the vertical dust flux, F, is calculated using the dust-emission model of Lu and Shao (1999) as

冢

冪p

Cagfrb F⫽ 0.24⫹Cbu∗ 2p

rp

冣

Q

(3)

1/2

with br=CR/CS, where CR is the drag coefficient for isolated roughness elements and CS is that for the surface; sr is the basal-to-frontal area ratio and mr is a parameter ⱕ1 accounting for non-uniformity in the surface stress. Raupach et al. (1993) suggested that br⯝90, mr⯝0.5, and s⯝1 are typical values. However, it is difficult to obtain the parameters like br and m for a specific cover condition. Thereafter, those typical values are used for all conditions. The system uses a simple empirical expression derived by Shao et al. (1996) H(w)⫽e22.7w The reason for using this simple empirical equation is due to its robustness and the fact that no additional parameter is required. Lu (2000) has discussed the difference between several available formulae and R and H. A satisfactory expression for M is not available at this stage and hence the values of M are temporarily set to 1 for all soils.

where Q is the horizontal saltation flux, f is the total volumetric fraction of dust in the sediment, rb is the bulk soil density, rp is the particle density, u∗ is the friction velocity, p is the soil plastic pressure exerted by soil on a particle moving through it, Ca and Cb are coefficients of order 1. The effects of surface non-erodible elements are only considered during the calculation of Q. Other effects, such as saltators impacting on non-erodible elements, are not modelled. An integration of Eq. (3) over the sand particle-size range gives the total vertical dust flux from the soil surface for a given u∗. For most natural soils, the plastic pressure, p, is larger than 105 N m⫺2 (Lu and Shao, 1999). This allows a further simplification of Eq. (3) to 0.12Cagfrb F⫽ Q p

(4)

which is used in this study. 2.1.4. Definition of dust The upper limit of dust-particle size, dd, particles with diameters smaller than which can remain in suspension once lifted from the surface, is defined as the solution of

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

wt(d)⫽0.7u∗

(5)

as proposed by Gillette (1974), where wt(dd) is the settling velocity for particles with diameter d. Other definitions of dd can be found in the literature with the ratio of wt(d1)/u∗ varying between 0.2 and 1.25 (Pye, 1987; Shao et al., 1996). Fig. 4 shows the dependence of dd on u∗, with wt(dd) being set to 0.2u∗, 0.7u∗ and 1.25u∗. Particles as large as d=200 µm have been detected in suspension during dust-storm events (Westphal et al., 1987). The choice of wt(dd)=0.7u∗ gives dd=100 µm for u∗=0.8 m s⫺1, which is reasonable for dust-storm conditions (Gillette, 1974). 2.2. Transport model If the concentration of the ith dust-particle group is ci=c(di), then the total dust concentration is given by

冘 N

ctotal=

ci. The transport scheme predicts dust concen-

i⫽1

tration by solving the dust concentration equation (Eq. (6)) with the boundary conditions given by Eq. (7). In the s-coordinate system, these equations can be written as ∂psci ∂psuci ∂psvci ∂ci ⫹ ⫹ ⫹ (pss˙ ⫹grwti) ∂t ∂x ∂y ∂s

(6)

∂ ∂ ∂ci/r ∂ci/r ∂ci/r g2 ∂ ⫽ps Kphr ⫹ps Kphr ⫹ Kpzr3 ∂x ∂x ∂y ∂y ps ∂s ∂s with boundary conditions g2 ∂ci/r ci(pss˙ ⫹grwti)⫺ Kpzr3 ⫽grFi ps ∂s

at the surface (7)

∂ci/r ⫽0 at the top ∂s where wti and Fi are the particle settling velocity and the surface vertical dust flux for the ith group, respectively; u, v and s˙ are wind velocities, ps is the surface pressure and r is the air density. The horizontal particle diffusivity Kph is assumed to be equal in the x and y directions. The vertical particle diffusivity Kpz is a function of particle size, estimated through a modification of the eddy diffusivity for neutral particles. The relationship between the eddy diffusivities for neutral and heavy particles is complicated in reality, but can be simplified for stationary, homogeneous and isotropic turbulence (Csanady, 1963). While there are more recent models on particle diffusivity, it can be shown that the Csanady (1963) model is sufficiently accurate for our purposes and has therefore been used in this study. Dust particles are removed only by dry deposition at the surface. Although knowing that wet removal is an important process, a large uncertainty exists in available wet removal schemes and the resolution of particle size distribution does not allow us to achieve a proper estimate of the scavenging rate. The value of modelling wet removal is limited and we excluded it from this study. Dry deposition is represented in the model as a boundary condition in the finite difference form of the vertical turbulent transport term. The dry removal flux by gravitational settling and turbulent mixing is modelled as Fd⫽C(r)(wt⫹vd) and is included in Eq. (7). The dry-deposition velocity, vd, is parameterised following Genthon (1992). The settling velocity, wt, can be estimated using

冉

wt(d)⫽

Fig. 4. Upper limit particle diameter dd, under which particles can be suspended in the air, defined as the solutions of wt(d1)/u∗=0.2, 0.7 and 1.25.

237

4rpgd 3rCD(Ret)

冊

1/2

(8)

where Ret=wtd/n is the particle Reynolds number at the 24 settling velocity, and CD(Ret)= (1+0.15Re0.687 ) is the t Ret drag coefficient for a sphere (Durst et al., 1984). Six dust particle size groups are considered in this study. These are, dⱕ2 µm (clay), 2⬍dⱕ11 µm (fine silt), 11⬍dⱕ22 µm (medium silt), 22⬍dⱕ52 µm (large silt), 52⬍dⱕ90 µm (fine sand), and 90⬍dⱕ125 µm (medium sand). The selection of the groups is due to the consideration of observation on atmospheric dust size distribution and the availability of the soil particle size distribution (PSD). It is observed that particles at around 1 and 10 µm are found as dominant classes in the atmosphere (Joussaume, 1990; Westphal et al., 1988). The lower bound of the sand particles is defined by particles of size 20 µm (Gillette, 1977; Gillette and Walker, 1977) and 60 µm particles are used under severe wind-erosion activities (Shao et al., 1996). However, the soil PSDs

238

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

data, which have 38 classes from 2 to 1159 µm do not use 1, 10, 20, and 60 µm as class bounds. Therefore, the nearest values (2, 11, 22, and 52 µm) are used. Two larger upper bounds (90 and 125 µm) are used to account for extreme events. During the simulation, the real upper bound is dynamically determined by the model itself using Eq. (5). The effective settling velocity for each dust section is calculated as

冕

wt(d)f(d)dd

w˜t⫽

冕

p(d)dd

where f(d)=13prr3dN/dd with d being the particle diameter and dN/dd being the number distribution. Following the work of Tegen and Fung (1994) and assuming that dN/dlogd⬀d −3 for dⱕ22 µm, and dN/dd⬀d −3 for 22ⱕ dⱕ125 µm, the effective settling velocity for each class can be calculated. The values are given in Table 1. These values are in good agreement with several measurements (Westphal et al., 1987). Eq. (6) is solved by splitting the advection and diffusion terms. A multidimensional wave-propagation slope-limiter scheme (LeVeque, 1996) is used for horizontal advection. This scheme is second order accurate both in time and space. It eliminates oscillations and undershoots and maintains positivity, which is important for dust concentration. In the vertical direction, Bott’s (1989) advection scheme is used. The scheme is positive definite but not monotonic. It is mass conservative and causes very small numerical diffusion. Second order, area preserving polynomials are used inside the domain. These polynomials were derived assuming variable grid spacing. The order of the polynomials is reduced to one near the domain boundary. The vertical diffusion is solved by using a fully implicit scheme, with the associated algebraic equation system being solved by using the Thomas algorithm. Emission and dry deposition are handled together with the vertical diffusion. 2.3. Data handling 2.3.1. Input parameters: access to the GIS database The land-surface parameters required for the model are: soil type, vegetation type, vegetation height, leafarea index (LAI) and land-use index. These parameters are obtained from the Atlas of Australian Resources,

(AAR) (Volumes 1, 3 and 6), except for LAI. Leaf-area indices used for the simulation period are derived from remotely sensed NDVI (Normalised Difference Vegetation Index) data, using empirical relationships (McVicar et al., 1996). The GIS database is of 0.05° spatial resolution. 2.3.1.1. Soil parameters In the Atlas of Australian Resources — Soils, soil types are aggregated into 29 mapping units. Seven of these units are well structured and stabilised soils and three are rocks, peats and saline lakes. These 10 units are non-erodible. The remaining 20 units are reclassified according to the availability of particle-size distributions and chemical similarities. Eight particle-size distributions (PSDs) of typical Australian soils provided by Dr McTainsh are assigned to the 20 erodible mapping units. In these PSDs, the smallest resolvable particle size is 2 µm. No attempt has been made to model the emission of smaller particles. The possible variation of particle sizes with wind velocity has also been neglected. In addition to PSDs, several soil parameters determine the intrinsic ability of a soil to resist wind erosion. These include soil bulk density rb, the plastic pressure p (also known as penetrometer resistance), Ca and Es. These parameters are assigned on the basis of the likelihood of crust formation, the chemical composition, structure and permeability of the soils. The assigned values are as listed in Table 2. 2.3.1.2. Vegetation and roughness Three types of vegetation data are used, namely, LAI, vegetation type (Vt) and vegetation height (Hc). Frontal-area index, which is required for the computation of threshold friction velocity, is estimated from LAI (Shao et al., 1994). The value of Ev is calculated from Ev⫽1⫺exp(⫺LAI/2) assuming a random distribution of foliage above the soil and uniform leaf-angle distribution (Choudhury, 1989). 2.3.2. Preprocessor: GIS converter The horizontal grid spacing used in the atmospheric model is often too coarse to resolve effects of the subgrid scale land-surface heterogeneity on dust emission. The explicit subgrid scale treatment of heterogeneity using the nesting method based on the resolution of the GIS data requires running the wind-erosion model for more than 900×700 nodes over the Australian continent, which

Table 1 Effective settling velocity for each particle-size group Section (µm)

dⱕ2

2⬍dⱕ11

11⬍dⱕ22

22⬍dⱕ52

52⬍dⱕ90

90⬍dⱕ125

wt (m s⫺1)

1.8×10⫺4

2.86×10⫺3

1.876×10⫺2

9.896×10⫺2

0.4259

0.8

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

239

Table 2 Soil type classification and the values of assigned parameters Mapping unit Bb1 Bb2 Bb3 Bb4 Bb5 Bc1 Bc2 Bd1 Bd2 Bd3 Ca1 Cb1 Cb2 Cc1 Cd1 Cd2 Cf1 Cf2 Cf3 Cf4 A1, A2, Ba1, Ba2,

Sample site (soil type) Charleville (massive sesquioxide) Charleville (massive sesquioxide) Charleville (massive sesquioxide) Charleville (massive sesquioxide) Charleville (massive sesquioxide) Buronga (calcareous earth) Buronga (calcareous earth) Betoota (saline soil) Betoota (saline soil) Betoota (saline soil) Diamantina Lakes Dune (sandy) Narrabri (cracking clay) Narrabri (cracking clay) Lark Quarry (loam sandy) Tambo (hard-setting soil) Tambo (hard-setting soil) Diamantina Lakes Dune (sandy) Diamantina Lakes Dune (sandy) Winton (shallow loam soil) Winton (shallow loam soil) Ba3, Ce1, Ce2, Cf5, O1

p (105 N/m2)

rb (kg/m3)

0.5 1 0.5 100 100 25 20 20 300 300 25 200 200 100 25 25 25 25 50 100

is numerically too expensive to implement in the weather-prediction model. A preprocessor is therefore designed to extract land-surface parameters from the GIS database (with 0.05° resolution) into the Lambert conformal projection-based grid of the atmospheric model. The preprocessor minimises the information loss of the original land-surface data and allows computational efficiency of the modelling system. Two important variables display subgrid scale variations, namely, friction velocity u∗, which controls the magnitude of surface wind stress, and threshold friction velocity u∗t, which controls the resistance of the surface to wind erosion. The local values of u∗ and u∗t are probably distributed about their mean values according to certain probability density functions. Thus, for a given grid, u∗ may exceed u∗t locally and temporarily even when u∗¿u∗t. The effect of this phenomenon has been investigated by Westphal et al. (1988). They found that the calculated dust emission can double if the Rayleigh distribution is assumed for u∗. Assuming u∗t to be uniformly distributed between 0.25 and 1.5 m s⫺1, the dust flux can be 10 times larger for the same u∗. Hence, the subgrid scale distribution of u∗t is of central importance to the calculation of wind erosion. Although the enlargement function R is time dependent, it is reasonable to treat u∗td as stationary if a monthly or short time prediction is of interest and the effect of soil moisture on u∗td is considered during an actual model run. For each GIS grid, R and Ev are calculated using the mosaic approach. In this approach, each atmospheric model grid is divided into several fractions according to the soil type. The fractions with the same

1000 1000 800 900 1000 1000 1000 1000 1200 1000 1000 1000 1000 1200 1000 1000 1000 1000 1000 1000

Ca

Es

5 5 2 2 1 5 5 5 1 5 5 5 5 2 5 5 5 5 2 1 Non-wind-erodible

1 1 1 1 0.4 1 1 1 1 1 1 1 1 1 1 1 1 0.8 0.8 0.8

soil index are added together regardless of their location within the grid. Typically, each atmospheric model grid contains 1–7 soil types. For each fraction, the parameters needed for the calculation of R are averaged. For a given atmospheric grid, a set of land-surface parameters is generated as k Is1 Is2 % Isk R1 R2 % Rk Ev1 Ev2 % Evk where Is is soil index and k is the number of fractions. These parameters are input to the wind-erosion model and an array u∗t(i,j) is calculated, with i=1,...,k and j=1,...,np, where np is the number of particle-size groups. The dust flux is calculated for each surface fraction and each particle-size group. The total dust flux for the given atmospheric model grid is given by

冘冘 k

np

F⫽

Fij

i⫽1j⫽1

where Fij is the dust flux for the ith surface fraction and the jth particle-size group. No subgrid distribution is considered for the friction velocity u∗. 3. Modelling results IWEMS has been applied to the simulation of the February 1996 dust-storm events over the Australian conti-

240

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

nent (Fig. 5). A 50 km grid spacing is used (for the atmospheric model) and the output is written in 6-hourly intervals. The integration time step is set to be 120 s. Throughout the simulation, the Australian Eastern Standard Time (EST=UTC+11 h) is used. The simulation began at 23:00 31 January 1996 and ended at 23:00 28 February 1996. The initial dust concentration was set to zero. The initial soil moisture was obtained by running the atmospheric model coupled with the land-surface scheme (without IWEMS) for the entire month of January 1996, starting from a typical summer-time soilmoisture pattern. 3.1. Prediction of dust emission Fig. 6 shows the simulated patterns of daily averaged dust emission. The days free of dust emission were excluded in the graph. Wind erosion occurred at various locations over the Australian continent for the remaining 21 days, including three major dust-storm events. The areas of intensive wind erosion were found mainly in the southern part of the Simpson Desert to the north of Lake Eyre and areas of medium dust uplifting were identified in the Great Victoria Desert, Gibson Desert and the west coast of Western Australia. Relatively weak, intermittent wind erosion activities were predicted in areas around Broken Hill and in Western Australia. To facilitate descriptions, we denote the 6-hourly

Fig. 5.

averages of near surface dust concentration, streamwise ¯ and F¯. sediment flux and dust emission rate as C¯, Q According to the model simulation, the first dust storm occurred on 1–2 February. An extensive area of the Australian continent experienced dust uplifting during this period, with regions in Western Australia being worst affected. However, dust emission for the particular event was not very high and the wind erosion was mainly intermittent. The dust emission in the Simpson Desert was weak with the daily averaged emission rate well below 10 µg m⫺2 s⫺1. During the time, the highest ¯ and F¯ were estimated to be around 7255 µg m⫺3, C¯, Q 200 mg m⫺1s⫺1 and 0.4 mg m⫺2s⫺1, respectively, at (125.9E°, 25.6S°) over the time period between 5:00 and 11:00 on 2 February. The most severe dust storm for the simulation period occurred between 8 and 10 February in Central Australia. The dust storm was associated with the strong winds generated by a deep low-pressure system which was moving from the Southern Ocean toward ¯ and F¯ were estimated to the northeast. The highest C¯, Q be around 9528 µg m⫺3, 1890 mg m⫺1s⫺1 and 0.7 mg m⫺2s⫺1, respectively, at (136.8E°, 27.9S°) during the time period from 5:00 on 9 February to 11:00 on 10 February. The instantaneous near surface dust concentration was as high as 27,060 µg m⫺3. The third, somewhat weaker, dust storm occurred between 14 and 16 ¯ and F¯ February, in the same region. The highest C¯, Q were estimated to be around 3236 µg m⫺3, 573 mg

Map of Australia showing state names, locations of major deserts, lakes and towns referred to in this paper.

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

Fig. 6.

Predicted daily averaged dust emission patterns. F is in µg m⫺2 s⫺1. The head of “010296” reads as 1 February 1996 and so on.

241

242

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

particle-size groups. This suggests that there is a strong correlation between the particle-size distribution of the emitted dust and that of the parent soil. In Fig. 7, no such similarity can be found. The smoother curves for d⬍11 µm suggest that the fine particles can remain in suspension longer than larger particles. Larger particles are found mostly during the wind erosion events. 3.3. Time evolution of variables at given grid points

Fig. 7. Time evolution of the domain averages of total suspended dust for each particle-size group.

m⫺1s⫺1 and 0.2 mg m⫺2s⫺1, respectively, at (136.8E, 27.9S°) at 11:00 on 14 February. The highest instantaneous near surface dust concentration was over 5000 µg m⫺3.

More detailed information can be gained by analyzing various model variables for a single grid point. Fig. 9 shows the evolution of dust concentration C (in log10 µg m⫺3), Q and simulated weather parameters for the grid point located at (126.8E°, 27.9S°). The diurnal variations of wind speed U, temperature T and u∗ can be clearly seen. During the simulation period, two severe dust-storm events occurred at this location, associated with the strong near surface wind.

3.2. Budget of dust concentration The evolution of the 6-hourly total amount of (a) suspended dust for each particle-size group and (b) the emitted dust from each group are shown in Figs. 7 and 8. A comparison of Figs. 7 and 8 reveals that although a large amount of particles with d⬎52 µm were emitted, only a small amount of those particles were suspended in the air due to their large settling velocities. Fig. 8 shows the similar shapes of dust-emission rates for all

Fig. 8. Time evolution of the domain averages of total dust emission for each particle-size group.

Fig. 9. Predicted time evolution of weather parameters and the smallest threshold friction velocity u∗t, streamwise sediment flux Q, and instantaneous near surface layer dust concentration for the grid point (136.8E°, 27.9S°).

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

The dust emission during the period between 9 and 10 February is produced by a south to southwesterly cold air stream driven by a frontal system. Temperature was increasing from 5 February and peaked on 7 and 8 February. The cold front generated a near surface wind as large as U=10 m s⫺1 corresponding to a u∗ of 0.45 m s⫺1 at 5:00 on 9 February. Air temperature during 5–8 February was very high and there was no rainfall, leading to low values of soil moisture. The smallest predicted threshold friction velocity is u∗t=0.33 m s⫺1. It was found that wind erosion started in the early morning of 9 February, peaked at 11:00 of the same day and eased at 17:00 of the next day. For the second dust-storm event, which occurred during the time period between 14 and 16 February, U reached 9.2 m s⫺1 and u∗ reached

243

0.4 m s⫺1. The predicted air temperature was below 25°C. The light rainfall associated with the cold front in the previous days did not increase soil moisture substantially. The threshold friction velocity u∗t was predicted to be around 0.35 m s⫺1. The maximum streamwise sediment flux is about one-third of that in the previous duststorm event. The wind erosion lasted for 2 days, starting in the early morning, peaking at noon time, then weakening at midnight. Although U at this grid was not the largest (the largest U in the computational domain was over 20 m s⫺1), severe wind erosion occurred at this grid point because of its low vegetation cover and low soil moisture, which provide ideal surface conditions for wind erosion. Again, this shows that wind erosion is controlled by a combination of surface properties and wind stress.

Fig. 10. Predicted monthly mean dust concentration in the first model layer for the six different size sections.

244

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

Fig. 11. Predicted vertical profile of dust concentration for six continuous days at 11:00 at 136.8E°.

3.4. Long-range transport Fig. 10 shows the modelled near surface (s=0.999) monthly mean concentration of the individual particlesize groups. Particles with d⬍11 µm were found to travel far beyond the surrounding oceans, while particles with d⬎52 µm were found to be deposited in areas surrounding the source region. The relative amount of smaller particles to larger particles increases with increasing distance from the sources, as the larger particles are removed more quickly due to their larger settling velocities. The horizontal distribution pattern of the small particles is more diffusive than of the larger particles. In reality, part of the fine particles will be washed out by rain and the actual travelling distance of those

fine particles might be smaller than simulated, as wet deposition is not included in the model. Fig. 11 shows the modelled vertical distribution of the instantaneous dust concentration at 11:00 for the 6 days from 9 to 14 February. For times when dust storms occurred, dust plumes reached higher levels on hot days (9th) than on cool windy days (14th). This confirms that convection plays an important role in the transport of particles into higher levels of the troposphere. For the same location and same altitude, dust concentration may change dramatically with time. For times with dust emission, dust concentration was higher than 1000 µg m⫺3. For times without dust emission, just 1 day after the storm, dust concentration was smaller than 1 µg m⫺3, because a large proportion of the large dust particles

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

lifted during the storm deposited back to the surface. They were removed within several hours from the atmosphere due to their large settling velocities. This effect can be seen more clearly in Fig. 12 which shows the instantaneous profile of dust concentration for the individual particle-size groups at grid point (136.8E°, 27.9S°) for two different times. The larger particles can only be found for the wind-erosion day and only reach the height with s=0.95. The change in the magnitude of dust concentration for particles with d⬍11 µm is in the order of 1–2, while for particles with d⬎11 µm, the change can be as large as several orders of magnitude. Therefore, particles suspended in the atmosphere for a long time are mainly those with d⬍11 µm.

245

If the residence time for each particle-size group is defined as tr⫽te⫺td where te is the time of particle emission and td is the time when all particles from this group are deposited, then tr is estimated to be 6–10 h for particles with 90⬍dⱕ125 µm, 12–18 h for 52⬍dⱕ90 µm, 1–2 days for 22⬍dⱕ52 µm, 10 days for 11⬍dⱕ22 µm and more than 1 month for dⱕ11 µm, if no wet removal is considered. These findings agree well with the various measurements mentioned by Westphal et al. (1987) and the model results of Tegen and Fung (1994). As shown previously, the size distribution of airborne

Fig. 12. Predicted vertical profile of concentration of six particle-size groups for a wind-erosion day (11:00 9 February 1996) and non-erosion day (11:00 21 February 1996) at the grid point (136.8E°, 27.9S°).

246

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

dust particles can change dramatically both in time and space. Near the surface, the dominant mode is 52⬍dⱕ 90 µm during wind-erosion and 2⬍dⱕ11 µm during periods free of wind-erosion. The dominant mode shifts gradually to the particle-size group of 2⬍dⱕ11 µm with increasing height during wind-erosion episodes, which is also the overall long-term dominant mode as can be seen in Fig. 7. These rapid and dramatic changes of size distribution of suspended particles in space and time imply that the assumption of static particle-size distribution in the atmosphere may lead to poor assessment of the effect of mineral aerosols on atmospheric radiation balance.

4. Comparison with observations Very few measurements are available for a quantitative verification of the model. Nevertheless, the predicted dust concentration can be compared with measurements from two sites. Fig. 13 shows that the model has produced reasonable estimates of dust concentration for Birdsville for 9 and 16 February, although the predicted values were smaller than the measured ones by a factor of two. However, the comparison for the Mildura site is rather unsatisfactory. The predicted dust concentration was 500 times smaller than the measured values though the model has predicted correctly the timing of high dust concentrations, except for 1 February. Several factors may have contributed to such a quantitative disagreement. The first one is the coarse resolution of the model. Although the subgrid variation of u∗t arising from the

heterogeneity of soil type and vegetation cover has been modelled, that arising from subgrid variations of friction velocity u∗, soil moisture and other quantities has not been taken into consideration. At the Mildura site, human-induced erosion, which happens at a much finer scale, is the dominant cause compared with a desert dust source area, such as Birdsville. The land surface conditions at Birdsville are more uniform compared with the Mildura site. In Australia, wind-erosion prone agricultural areas are relative small. They spread in the south part of the country with paddock sizes smaller than 1 km2. The 5 km GIS database is not capable of detecting human-induced wind erosion activities even though the subgrid variation of other variables are modelled. The second factor is that the representation of surface crusts, roughness elements and soil moisture may still be too crude. The third factor is the parameters, especially soil parameters, such as horizontal plastic pressure p, particle size distribution, and bulk density are insufficiently accurate for the entire continent. The predictions can also be compared with satellite images. Fig. 14 shows the comparison between the simulated dust emission rate and the visible light picture from the Geostationary Meteorological Satellite (GMS) for 00Z 9 February (11:00 EST). From the satellite image, dust clouds are visible in the region around 137°E and 28°S, which are identified by comparing visible light and near infrared satellite images. Fig. 15 shows the simulated results of several variables for the same time at continent scale. The simulated dust clouds, in general, show good qualitative agreement when compared with the visible features from the satellite picture and meteorological records during the episode. The shape, location and extent of the dust cloud coincide well with the wind-erosion areas predicted by the model. The most intensive wind-erosion areas are well identified.

5. Comparison with other model results

Fig. 13. Measured and modelled dust concentrations for Birdsville and Mildura.

Shao and Leslie (1997) have modelled the dust storms during 8–11 February. In their study, u∗t was calculated for each atmospheric model grid using the dominant soil type and the average value of LAI with no subgrid treatment. Their prediction of wind erosion was done off-line using a predicted wind field at half hourly intervals. While there is a general agreement between Shao and Leslie (1997) and this work, the former predicted a stronger wind-erosion event both in affected areas and dust-emission rates. The largest vertical dust flux F predicted by Shao and Leslie (1997) is over 3 mg m⫺2 s⫺1, about five times that of the present work. Overall, the total dust emission calculated by Shao and Leslie was about 6 Mt for four intensive dust-storm days. For the same period, we found the total dust emission was about 1.113 Mt, again about five times smaller than their pre-

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

Fig. 14.

247

Predicted dust emission locations at the Central Australia and visible satellite image at the time 00Z 9 February 1996.

diction. There are two possible reasons for the differences between the two studies. The first one is that the dust-emission schemes used in the two studies differ. The dust-emission scheme used by Shao and Leslie (1997) is based on the study of Shao et al. (1993), while the one used here is based on the study of Lu and Shao (1999). Both models have quantitative uncertainties as discussed in Lu and Shao (1999). The scheme of Shao et al. (1993) was derived from wind tunnel experiments using a loosely packing dust bed as the saltating target. Compared with the field data of Gillette (1977), the scheme of Shao et al. (1993) appears to overestimate dust emission for large values of u∗ (Lu, 2000). The second possible reason lies in the fact that subgrid variations of u∗t were not considered by Shao and Leslie (1997). The values for LAI and vegetation cover used in their study may have resulted in a smaller u∗t and

hence larger streamwise sediment flux and dust-emission rate. The total amount of dust emitted from the continent is estimated to be 1.87 Mt in the simulated month. Comparing with the estimated global dust uplifting rates (ranging from 200 to 5000 Mt yr⫺1) (Pye, 1987; Tegen and Fung, 1994), it suggests that Australia is not a strong source of aeolian dust though it occupies over 10% of global desert area. This is consistent with other observations. Maximum dust storm frequencies in Australia are 15 day yr⫺1, compared with 30–60 days in Asia (Pye, 1987). Satellite data shows that the optical thickness over the oceans surrounding Australia is often less than 0.1 compared with 0.6–2 over the oceans near the Sahara. One of the causes of these differences in emission rate is that the level of human disturbance in the Saharan and Asian dust-source regions is much higher compared

248

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

Fig. 15. Predicted wind erosion intensity, near surface dust concentration, vertically integrated dust concentration, surface wind, friction velocity, and threshold friction velocity for the time 00Z 9 February 1996 dust-storm events.

with that in Australia. It was found that about 60–90% contemporary atmospheric dust is generated from disturbed regions (Tegen and Fung, 1995). The predicted wind-erosion areas coincide with the areas of largest dust-storm frequencies derived by McTainsh and Pitblado (1987) based upon meteorological records at 149 stations. This study also gives a good indication of the dust paths and deposition regions. Fig. 15 shows that dust particles emitted from Central Australia were transported toward the Indian Ocean in the west, passing offshore to the south-east. The dust plume was confined to a narrow region of the frontal area. For

the whole simulation month, it is found that the largest fraction of dust particles was transported toward the northwest. All these features are consistent with the climatology of wind erosion in Australia (Sprigg, 1982; McTainsh and Leys, 1993).

6. Discussion and conclusion An integrated wind-erosion prediction system has been briefly described with emphasis on the physically based wind-erosion scheme and its coupling with a GIS

H. Lu, Y. Shao / Environmental Modelling & Software 16 (2001) 233–249

database. Comparing with some existing dust-transport models, this system approach has the ability to detect the dust sources and to calculate dust-emission rate using both meteorological and surface information. The simulations of the February 1996 dust-storm events over the Australian continent have been carried out to demonstrate the suitability of the approach. The simulated evolution of dust storms has been found to be in qualitative agreement with the observations. Quantitative agreement has also been partially achieved, although still less than satisfactory. Some subjectivity has been inevitably involved in the derivation of soil parameters. Future work will be directed at improving the quality of the GIS database, enhancing the reliability of the model, studying the impact of dust concentration on the atmospheric radiation budget, and eventually creating a satisfactory modelling system for quantitative dust-storm prediction. Acknowledgements This work is supported by Australia Research Council. The authors are grateful to M.R. Raupach and P.A. Coppin at CSIRO Land and Water for insightful comments which greatly improved the original manuscript. References Andreae, M.O., 1995. Climate effects of changing atmospheric aerosol levels. In: Henderson-Sellers, A. (Ed.), World Survey of Climatology, Future Climates of the World, vol. 16. New York pp. 341–392. Bott, A., 1989. A positive definite advection scheme obtained by nonlinear renormalization of the advective fluxes. Mon. Wea. Rev. 117, 1006–1015. Choudhury, B.J., 1989. Estimating evaporation and carbon assimilation using infrared temperature data: vistas in modeling. In: Asrar, G. (Ed.), Theory and Applications of Political Remote Sensing. New York, pp. 628–690. Csanady, G.T., 1963. Turbulent diffusion of heavy particles in the atmosphere. J. Geophys. Sci. 20, 201–208. Duce, R.A., 1995. Sources, distributions and fluxes of mineral aerosols and their relationship to climate. In: Charlson, R.J., Heintzenberg, J., (Eds.), Aerosol Forcing of Climate. New York, pp. 43-72. Durst, F., Milojevic, D., Scho¨nung, B., 1984. Eulerian and Lagrangian predictions of particulate two-phase flows: a numerical study. Appl. Math. Modelling 8, 101–115. Genthon, C., 1992. Simulations of desert dust and sea salt aerosols in Antarctica with a general circulation model of the atmosphere. Tellus 44, 371–389. Gillette, D.A., 1974. On the production of soil wind erosion aerosols having the potential for long range transport. J. Rech. Atmos. 8, 735–744. Gillette, D.A., 1977. Fine particulate emission due to wind erosion. Trans. of the ASAE 20 (5), 890–897. Gillette, D.A., Walker, T.R., 1977. Characteristics of airborne particles produced by wind erosion on sandy soil, high plains of West Texas. Soil Sci. 123, 97–110. Joussaume, S., 1990. Three-dimensional simulation of the atmospheric cycle of desert dust particles using a general circulation model. J. Geophys. Res. 95, 1909–1941.

249

Junge, C., 1979. The importance of mineral dust as an atmospheric constituent. In: Morales, C. (Ed.), Saharan Dust. New York, pp. 49-60. Leslie, L.M., Purser, R.J., 1991. High-order numerics in a three-dimensional time-split semi-Lagrangian forecast model. Mon. Wea. Rev. 119, 1612–1632. LeVeque, R.J., 1996. High-resolution conservative algorithms for advection in incompressible flow. SIAM J. Numer. Anal. 33, 627–665. Lu, H., 2000. An integrated wind erosion modelling system with emphasis on dust emission and transport. Ph.D. thesis, The University of New South Wales, Australia. Lu, H., Shao, Y.P., 1999. A new model for dust emission by saltation bombardment. J. Geophys. Res. 104 (D14), 16827–16841. McTainsh, G.H., Leys, J.F., 1993. Soil erosion by wind. In: McTainsh, G.H., Broughton, W. (Eds.), Land Degradation in Australia. Longman Cheshire, Sydney, pp. 188–230. McTainsh, G.H., Pitblado, J.R., 1987. Dust storms and related phenomena measured from meteorological records in Australia. Earth Surface Processes and Landforms 12, 415–424. McVicar, T.R., Walker, J., Jupp, D.L.B., Pierce, L.L., Byrne, G.T., Dallwitz, R., 1996. Relating AVHRR vegetation indices to in situ measurements of leaf area index. Technical Report 96.5. CSIRO, Division of Water Resources. Owen, R.P., 1964. Saltation of uniform grains in air. J. Fluid Mech. 20, 225–242. Pye, K., 1987. Aeolian Dust and Dust Deposits. Academic Press, San Diego, CA. Raupach, M.R., 1991. Saltation layer, vegetation canopies and roughness lengths. Acta Mech. Suppl. 1, 135–144. Raupach, M.R., Gillette, D.A., Leys, J.F., 1993. The effect of roughness elements on wind erosion threshold. J. Geophys. Res. 98, 3023–3029. Rice, M.A., Willetts, B.B., McEwan, I.K., 1996. Wind erosion of crusted soil sediments. Earth Surface Processes and Landforms 21, 279–293. Shao, Y.P., Leslie, L.M., 1997. Wind erosion prediction over the Australian continent. J. Geophys. Res. 102, 30091–30105. Shao, Y.P., Lu, H., 2000. A simple expression for wind erosion threshold friction velocity. J. Geophys. Res. 105 (D17), 22437–22443. Shao, Y.P., Raupach, M.R., Findlater, P.A., 1993. The effect of saltation bombardment on the entrainment of dust by wind. J. Geophys. Res. 98, 12719–12726. Shao, Y.P., Raupach, M.R., Leys, J.F., 1996. A model for predicting aeolian sand drift and dust entrainment on scales from paddock to region. Aust. J. Soil Res. 34, 309–342. Shao, Y.P., Raupach, M.R., Short, D., 1994. Preliminary assessment of wind erosion patterns in the Murray–Darling basin. Aust. J. Soil Water Conserv. 7, 46–51. Sprigg, R.C., 1982. Alternating wind cycles of the Quaternary era, and their influence on aeolian sedimentation in and around the dune deserts of south eastern Australia. In: Wasson, R.J. (Ed.), Quaternary Dust Mantles of China, New Zealand and Australia. Proceedings INQUA Loess Commission Workshop, Canberra. Tegen, I., Fung, I., 1994. Modelling of mineral dust in the atmosphere: sources, transport, and optical thickness. J. Geophys. Res. 99, 22897–22914. Tegen, I., Fung, I., 1995. Contribution to the atmospheric mineral aerosol load from land surface modification. J. Geophys. Res. 99, 22897–22914. Westphal, D.L., Toon, O.B., Carson, T.N., 1987. A two-dimensional numerical investigation of the dynamics and microphysics of Saharan dust storms. J. Geophys. Res. 92 (D3), 3027–3049. Westphal, D.L., Toon, O.B., Carson, T.N., 1988. A case study of mobilisation and transport of Saharan dust. J. Atmos. Sci. 45, 2145–2175.