Estimating aerodynamic resistance of rough surfaces using angular reflectance

Remote Sensing of Environment 114 (2010) 1462–1470

Contents lists available at ScienceDirect

Remote Sensing of Environment j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / r s e

Estimating aerodynamic resistance of rough surfaces using angular reflectance Adrian Chappell a,⁎, Scott Van Pelt b, Ted Zobeck c, Zhibao Dong d a

CSIRO Land and Water, GPO Box 1666, Canberra, ACT 2601, Australia Agricultural Research Service, United States Department of Agriculture, Big Spring, TX, 79720, USA c Agricultural Research Service, United States Department of Agriculture, Lubbock, TX, 79720, USA d Key Laboratory of Desert and Desertification, Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, No. 320 West Donggang Road, Lanzhou 730000, Gansu Province, PR China b

a r t i c l e

i n f o

Article history: Received 23 November 2009 Received in revised form 27 January 2010 Accepted 30 January 2010 Keywords: Dust emission model Wind erosion Sheltering Erodible Flow separation Drag Wake Aerodynamic resistance Aerodynamic roughness length Shadow Illumination Ray-casting Digital elevation model Roughness density Frontal area index Angular reflectance Bi-directional reflectance

a b s t r a c t Current wind erosion and dust emission models neglect the heterogeneous nature of surface roughness and its geometric anisotropic effect on aerodynamic resistance, and over-estimate the erodible area by assuming it is not covered by roughness elements. We address these shortfalls with a new model which estimates aerodynamic roughness length (z0) using angular reflectance of a rough surface. The new model is proportional to the frontal area index, directional, and represents the geometric anisotropy of z0. The model explained most of the variation in two sets of wind tunnel measurements of aerodynamic roughness lengths (z0). Field estimates of z0 for varying wind directions were similar to predictions made by the new model. The model was used to estimate the erodible area exposed to abrasion by saltating particles. Vertically integrated horizontal flux (Fh) was calculated using the area not covered by non-erodible hemispheres; the approach embodied in dust emission models. Under the same model conditions, Fh estimated using the new model was up to 85% smaller than that using the conventional area not covered. These Fh simulations imply that wind erosion and dust emission models without geometric anisotropic sheltering of the surface, may considerably over-estimate Fh and hence the amount of dust emission. The new model provides a straightforward method to estimate aerodynamic resistance with the potential to improve the accuracy of wind erosion and dust emission models, a measure that can be retrieved using bi-directional reflectance models from angular satellite sensors, and an alternative to notoriously unreliable field estimates of z0 and their extrapolations across landform scales. Crown Copyright © 2010 Published by Elsevier Inc. All rights reserved.

1. Introduction Soil-derived mineral dust contributes significantly to the global aerosol load. The direct and indirect climatic effects of dust are potentially large. A prerequisite for estimating the various effects and interactions of dust and climate is the quantification of global atmospheric dust loads (Tegen, 2003). Recent developments in global dust emission models explicitly simulate areas of largely unvegetated dry lake beds as sources of preferential dust emission (Tegen et al., 2002, 2006; Mahowald et al., 2003). In the case of the Earth's largest source of dust (Bodélé Depression; Warren et al., 2007) there are some significant discrepancies between ground measurements of dust emission processes and model assumptions (Chappell et al., 2008). Dust emission is produced by two related processes called

⁎ Corresponding author. E-mail address: [email protected] (A. Chappell).

saltation and sandblasting. Saltation is the net horizontal motion of large particles or aggregates of particles moving in a turbulent nearsurface layer. Sandblasting is the release of dust and larger material caused by saltators as they impact the surface (Alfaro & Gomez, 1995; Shao, 2001). Naturally rough (unvegetated) surfaces usually comprise a heterogeneous mixture (size and spacing) of non-erodible roughness elements that reduce the area of exposed and hence erodible substrate. When such rough surfaces are exposed to the wind, wakes or areas of flow separation (Arya, 1975) are created downwind of all obstacles. These sheltered areas reduce the area of exposed substrate still further and protect some of the roughness elements from the wind (depending on their size and spacing). This heuristic formed the basis for the dimensional analysis of the Raupach (1992) model where dynamic turbulence was replaced by a concept of effective shelter area and was portrayed as a wedge-shaped sheltered area in the lee of the element. The size and shape of the sheltered area is influenced by the wind velocity (speed and direction) and the heterogeneous nature of the surface (Fig. 1). Consequently, the erodible area and the non-

0034-4257/$ – see front matter. Crown Copyright © 2010 Published by Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2010.01.025

A. Chappell et al. / Remote Sensing of Environment 114 (2010) 1462–1470

Fig. 1. Cylinders used to represent non-erodible roughness elements in wind tunnel studies and parameterizations for wind erosion and dust emission models protect a portion of the substrate surface that may include all or part of other roughness elements in a heterogeneous surface (a). A change in wind direction redefines the area of the substrate protected from the wind and may expose previously protected roughness elements (b).

erodible roughness elements that are exposed to, and protected from, drag are an anisotropic function of the heterogeneous surface and wind speed. Central to wind erosion and dust emission models is the turbulent transfer of momentum from the fluid to the bed. The key assumption made by dust emission models (e.g., Marticorena and Bergametti, 1995; p. 16,418) is that the momentum extracted by roughness elements is controlled primarily by their roughness density (λ; Marshall, 1971) and consequently the erodible area is that which is not covered by roughness elements. The λ (also known as lateral cover or the frontal area index) is expressed as λ = nbh / S where n is the number of roughness elements inside an area (or pixel) S and b and h are the breadth and height, respectively of the roughness elements. This assumption forms one of the foundations for the dust production model (Marticorena and Bergametti, 1995) and dust emission scheme (Marticorena et al., 1997) upon which many dust emission models are based (e.g., Tegen et al., 2006). The approach assumes that the roughness elements cover part of the surface, protect it from erosion and that they consume part of the momentum available to initiate and sustain particle motion by the wind. The assumption manifests itself in dust emission models (e.g., Marticorena and Bergametti, 1995; p. 16,422; Eq. 34): Gtot = EC



 ρa 3 2 U ∫ ð1 + RÞ 1−R dSrel Dp dDp g * Dp

ð1Þ

as the ratio of the erodible area to the total surface area (E) and is set to 1 in the absence of information about non-erodible roughness elements and of vegetation and snow (Tegen et al., 2006). The parameter C is a constant of proportionality (2.61), ρa is the air density, g is a gravitational constant, U3* is the cubic shear stress of the Prandtl–von Karman equation where U* = u(z)(k / ln(z / z0)) and u is the wind speed at a reference height z, k is von Karman's constant (0.4) and z0 is the aerodynamic roughness length. The threshold friction velocity defines R = U*t(Dp, z0, z0s) / U* where the threshold shear stress U*t is a function of particle diameter Dp, z0 and the aerodynamic roughness length of the same surface without obstacles (z0s). The dSrel is a continuous relative distribution of basal surfaces formed by dividing the mass size distribution by the total basal surface and dDp is the particle diameter distribution. This approach includes neither a sheltering effect nor any interaction between the momentum extraction of the roughness elements and the downwind substrate area that they protect (wake). Furthermore, R implicitly assumes homogeneous surface roughness and it does not account for the anisotropy of heterogeneous surface roughness created by changing wind directions i.e., anisotropic z0. Wind erosion and dust emission models should reach a compromise between the realistic representation of the erosion/abrasion processes and the availability of data to parameterize or drive the model (Raupach & Lu, 2004). The requirement here is to reduce the complexity of aerodynamic resistance from an understanding of wake and shelter but capture the essence of the process to make reasonable estimates,

1463

particularly across scales of variation. For example, Shao et al. (1996) provided one of the first physically-based wind erosion models to operate across spatial scales from the field to the continent (Australia). One of the main reasons for its success was its approximation of λ using NDVI (Normalised Difference Vegetation Index) data. To improve this approximation Marticorena et al. (2004) argued that a proportional relationship existed between the protrusion coefficient (PC) derived from a semi-empirical bi-directional reflectance (BRF) model (Roujean et al., 1992) and geometric roughness. Although Roujean et al. (1992) stated the model's limitation for unvegetated situations and Marticorena et al. (2004) recognised this limitation, they developed a relationship between geometric roughness and z0. They retrieved the PC from surface products of the space-borne POLDER (POLarization and Directionality of the Earth's Reflectances) instrument and compared it to geomorphic estimates of z0 (Marticorena et al. 1997; Callot et al. 2000). The authors concluded that z0 could be derived reliably from the PC in arid areas. The main justification for the simplifying assumption of λ in wind erosion and dust emission models appears to be the hypothesis that the configuration and shape of non-erodible (unvegetated) surface roughness elements are unimportant for explaining the drag partition. The concept of drag or shear stress partitioning (Schlichting, 1936) is that the total force on a rough surface Ft can be partitioned into two parts: Fr acting on the non-erodible roughness elements and Fs acting on the intervening substrate surface Ft = Fr + Fs. There is a growing body of evidence that supports this approach. For example, Marshall (1971) studied drag partition experimentally in a wind tunnel and showed no difference between cylinders placed on a regular grid, on a diagonal or at random across the wind tunnel (λ = 0.0002 to 0.2). Raupach et al. (1993) reached a similar conclusion after inspecting Marshall's data and believed that there was only a weak experimental dependence of stress partition on roughness element shape and the arrangement of elements on the surface. Drag balance instrumentation used by Brown et al. (2008) in a wind tunnel, independently and simultaneously measured the drag on arrays of cylinders and the intervening surface, separately. Results were interpreted as confirmation that an increase in surface roughness enhanced the sheltering of the surface, regardless of roughness configuration i.e., irregular arrays of cylinders were analogous to staggered configurations in terms of drag partitioning. The role of flow separation and much-reduced drag in sheltered regions, particularly downwind of roughness elements, is significant for drag partitioning. We posit that the sheltered area is required to account for anisotropic variation in aerodynamic resistance for realistic wind erosion and dust emission models. Furthermore, we posit that current estimates of the erodible area using the area not covered by protruding objects is a poor representation of the erodible substrate exposed to abrasion from mobile material. The aim of the paper is to describe and evaluate the basis for using angular reflectance data to quantify the geometric anisotropy of aerodynamic resistance, account for heterogeneity and estimate the area exposed to abrasion. 2. Estimating aerodynamic roughness length (z0) and erodible area using shadow 2.1. Relationship between reflectance and frontal area index (λ) A new approach is presented here which is based on Chappell and Heritage (2007). The approach is inspired by the dimensional analysis of the Raupach (1992; p. 377–378) model (effective shelter area) and its replacement of dynamic turbulence with the scales controlling an element wake and how the wakes interact (Shao and Yang, 2005) and by the heuristic model of Arya (1975) and hence its similarity with the scheme of Marticorena and Bergametti (1995). In common with Marticorena et al. (2004), we show the relationship between reflectance and aerodynamic resistance estimated by wind tunnel

1464

A. Chappell et al. / Remote Sensing of Environment 114 (2010) 1462–1470

studies of aerodynamic roughness length (z0) and explain the relationship between reflectance and the frontal area index (λ). The λ is the projection of an obstacle's frontal area onto a predefined area or pixel with a flat surface. The projection is defined for a 45° illumination zenith angle i such that tan i is used as a multiplication factor which in this case is restricted to 1 and thus the entire frontal area of the object is projected. If 0°bib45° then tan i reduces the projected frontal area and when 45°bib90° then tan i increases the projected frontal area. If that pixel is viewed at nadir and illuminated for different i, the shadow cast by the object is the same as the projection of the object onto the pixel for the given i. The light reflected from the rough surface is reduced by the proportion of the area that is in shadow and visible (Fig. 2). In the case of many homogeneous objects the shadow area may be reduced if the spacing between the objects is insufficient to allow the shadow cast to reach the underlying surface (Fig. 2). In other words the shadow is projected onto adjacent objects (mutual shadowing). However, if objects with vertical sides (e.g., cylinders) are homogeneous across an area their shadow is truncated because although it is projected on to the adjacent objects it is not visible at nadir. Illumination of natural surfaces demonstrates that a portion of the surface that may otherwise cast shadow may also be under shadow and this effect is dependent on illumination azimuth relative to an arbitrary origin (Fig. 3). Since we believe a priori that the geometry of a rough surface influences aerodynamic resistance and that roughness is also one of the main controls on the proportion of illumination there should be a relationship between the aerodynamic resistance and illumination proportion (viewed at nadir). 2.2. Relationship between shadow and erodible area The convention within current dust emission models is to approximate the erodible area as simply the intervening surface not covered by non-erodible elements. However, saltating soil particles usually strike the soil surface at an elevation angle of approximately 12°–15° (Sorensen, 1985). The point of impact is influenced by the particle jump length, angle of descent and surface roughness (Potter et al., 1990). As the particles bounce, they may jump over obstructions on the surface. If the obstruction is sufficiently tall that the particle cannot jump over it and it is non-erodible, it will shelter a portion of the soil surface from abrasion. Thus, upwind obstructions determine the angle of trajectory shelter angle a particle must achieve in order to strike a given point within the horizontal bounce distance of the saltating particle. The fraction of the soil surface impacted by saltating grains varies with the fraction of the surface sheltered by nonerodible roughness elements. Potter et al. (1990) developed the

cumulative shelter angle distribution and it was included in the Wind Erosion Prediction System to make daily estimates of wind erosion (Hagen, 1990). Here we propose to use the approach by Potter et al. (1990) and developed by Zobeck and Popham (1998) to approximate the erodible area. The curves formed by field measurements by Zobeck and Popham (1998; Figs. 2 and 3) are a function of the surface roughness, the shelter angle is equivalent to the illumination zenith angle described above and the proportion of the surface illuminated and viewed at nadir is equivalent to the surface fraction. Consequently, we propose to approximate the erodible area by predicting the proportion of surface in shadow with an illumination zenith angle of 75° (equivalent to an elevation angle of 15°) and viewed at nadir. 2.3. Evaluation methodology A digital elevation model is used here to reconstruct the surface roughness configuration of previous laboratory wind tunnel and field studies for which shadow/illumination measurements were not available. A ray-casting approach demonstrates the means by which shadow can be estimated remotely. It makes use of a fine resolution of elevation sufficient to discretise the surface obstacles. Such an approach is able to handle heterogeneous bed situations where the object heights and spacings are different across a surface. Variable illumination orientation, to represent wind direction, was accounted for with the same computational procedure and an illumination azimuth angle ϕ relative to a fixed arbitrary origin. A number of models exist to estimate the proportion of reflectance from a rough surface (Cierniewski, 1987; Hapke, 1993; Li & Strahler, 1992; Liang & Townshend, 1996). The most suitable model for describing reflectance of surface roughness across scales is the subject of ongoing research by the authors. An empirical function is used here as an alternative to these models for simplicity and to avoid the modeling becoming a distraction from the retrieval of aerodynamic resistance information. The proportion of illuminated surface, viewed at nadir for given illumination zenith and azimuth angles, was fitted with a Gaussian model with an additional parameter. The function has the desirable quality of resembling a positive exponential or a Gaussian model. Its isotropic form and model parameters are: ( SðiÞ = c exp −

α

i α r

!) ;

ð2Þ

where S is the proportion of illuminated surface for a given illumination zenith angle (i). The function approaches its sill (c) asymptotically and does not have a finite range. Instead, the distance parameter r defines the extent of the model. For practical purposes it can be regarded as having an effective range of approximately (3r)− α,

Fig. 2. Digital elevation models of 38 mm diameter and 24 mm high hemispheroids on a diagonal lattice with a coverage of 5% (a) and 75% (b) using an illumination zenith angle of 75° to cast the shadow (darkest tone) and represent the area sheltered from abrasion (see text for details).

A. Chappell et al. / Remote Sensing of Environment 114 (2010) 1462–1470

1465

Fig. 3. Digital elevation model of a 1 m2 natural soil surface with an overlay of shadow (dark tones) representing illumination azimuth angles of (a) 0°, (b) 45°, (c) 90° and (d) 135° for an illumination zenith angle of 75° representing the abrasion angle (see text for details). The proportion of shadow viewed at nadir covering the surface is 74%, 65%, 63% and 56% for each azimuth direction, respectively.

where it reaches 95% of its sill (Webster & Oliver, 2001). The additional parameter α describes the intensity of variation and the curvature. If α = 1 then the function is a positive exponential and if α = 2 the function is a Gaussian. As i → 90°, c enables S to be greater than 0 to represent the illumination of objects above a rough background or substrate. The proportion of illuminated surface associated with the illumination zenith angle 75° is used to approximate the erodible area (Section 2.2). The goodness of fit was assessed using the RMSE which is here defined as the square root of the mean squared difference divided by the number of degrees of freedom (df). The df is the number of data minus the number of model parameters used. The brightness of the rough surface was assumed to be the same as that of the background (substrate) and in the case of the reconstructed wind tunnel studies to scatter according to the Lambertian distribution. In this case, the Gaussian function was integrated over all illumination zenith (i) and azimuth (ϕ) angles and made relative to the incoming radiation (I) to form a single scattering albedo (SSA): ϕ = 2π i = π = 2

SSA = I =

∫

∫

ϕ=0

i=0

(

α

i c exp − α r

!) :

ð3Þ

In addition to the SSA a calculation is made for a specified direction over all i but for only a single ϕ viewed at nadir and made relative to the incoming radiation. This statistic is here defined as the relative directional reflectance (RDR) for use as a measure of the geometric anisotropy of the aerodynamic resistance. Although the ray-casting approach described here is evidently capable of handling the anisotropic nature of heterogeneous surface roughness, it is not intended to provide an operational method for the

retrieval of shadow. Surface illumination/shadow may be readily retrieved using angular sensors on airborne and space-borne platforms. A photometric model can be used to characterise the surface reflectance for given illumination and viewing conditions (cf. Hapke, 1993). The estimation of soil/surface roughness using shadow and geometric models is well established (Cierniewski, 1987) and approximations to radiative transfer theory by Hapke (1993) have provided parameterized models for soil bi-directional reflectance measurements (Pinty et al., 1989; Jacquemoud et al., 1992; Chappell et al., 2006, 2007; Wu et al., 2008). The description above amounts to a hypothesis that the proportion of illumination (viewed at nadir) can approximate the aerodynamic roughness length and the erodible area. That hypothesis is evaluated against existing measurements of aerodynamic resistance. 3. Evaluation data 3.1. Isotropic wind tunnel roughness A wind tunnel study by Marshall (1971) investigated cylinder and hemisphere roughness elements made from solid ‘varnished’ wood. The elements had a uniform height of 2.54 cm with a range of diameters (1.27, 2.54, 5.08, 7.62 and 12.7 cm). The elements were arranged in the working section of a wind tunnel on square and diagonal grids and also using randomly arranged patterns with spacing between elements (2.54–125.73 cm) that produced various coverage (ratio of cylinder base area to the specific cover: approx. 0.01–44.18 %). Each surface configuration was subjected to a single free-stream velocity at 20.3 m s− 1 and the total drag force of the surface roughness and that of the roughness elements separately was measured to deduce the surface drag force in between elements.

1466

A. Chappell et al. / Remote Sensing of Environment 114 (2010) 1462–1470

curve fitting gave a value of R2 b 0.98 at the 5% significance level. Only the data for 20 m s− 1 free-stream velocity was used here to ensure comparison with the results of Marshall (1971). These data are plotted against roughness density in Fig. 4 and the characteristics are summarized in Table 1. 3.2. Geometric anisotropic natural roughness

Fig. 4. The results of wind tunnel measurements by Marshall (1971) and Dong et al., (2002) showing aerodynamic roughness length (z0) standardized by object height (h) against frontal area index (λ; roughness density).

Aerodynamic roughness lengths (z0) of Marshall's surface configurations were derived using the same approach as Raupach et al. (2006; p. 214): " !# z0 δ U = exp k B− δ ; h h u*

ð4Þ

where δ ≈ 3.3+15.0(λϕ)0.43 (cm), k = 0.4, B = 2.5 and Uδ = 20.3 m s− 1 from table 4 of Marshall (1971). These data are plotted against λ in Fig. 4. Their characteristics are summarized in Table 1. An investigation of gravel surface aerodynamic resistance was conducted in a wind tunnel by Dong et al. (2002). Six types of artificial gravel were fabricated in cement to form “parabolic-shaped” elements with circular bases and diameters (19, 29, 38, 47, 57, 65 mm) and heights (12, 19, 24, 31, 37, 43 mm) with a diameter:height of approximately 1.5. The gravels were arranged in the working section of a wind tunnel in a diamond (staggered) pattern so that a range of coverage (1–92%) was provided. When the spacings of the objects with these coverages were evaluated in our model, coverage greater than 85% had a spacing of less than the model resolution (1 mm). Consequently, surface coverages greater than 85% were excluded in all subsequent analyses. Dong et al (2002) used ten free-stream velocities (4, 6, 8, 10, 12, 14, 16, 18, 20, 22 m s− 1) for each type of gravel and coverage. Wind profiles (the distribution of wind speed with height) were measured using an array of 10 Pitot-static probes mounted at ten heights (3, 6, 10, 15, 30, 60, 120, 200, 350 and 500 mm above the surface). Aerodynamic roughness lengths (z0) were derived from measured wind profiles using least squares curve fitting (Dong et al., 2002). Results were eliminated by Dong et al. (2002) if logarithmic

Wind velocity measurements on a circular field (100 m radius; 86,180 m2) at the USDA, Agricultural Research Service experimental farm in Lubbock, Texas were made between March and May 2001 to provide a validation dataset for a larger unpublished study. At the centre of the field a 2 m high mast was erected and equipped with calibrated cup anemometers at heights of 0.5, 1.0, and 2.0 m and with a hot wire anemometer at 0.01 m. Data were logged on a CSI 21X data logger programmed to record the average wind speed at 1 min integrated intervals until the 2 m wind speed exceeded 3.5 m s− 1 for 5 min at which time the integrated sampling interval was reduced to 1 s. Five minute average wind data were available and wind directions were merged into the data files. The parameter values (u* and z0) of the Prandtl equation were optimized against the 5 min average wind velocity using a non-linear weighted least squares routine (PROC NLIN in SAS version 8.2). The characteristics of these data are summarized in Table 2. 4. Digital elevation model reconstructions of wind tunnel and field surface conditions Surface roughness configurations were illuminated using raycasting coded in Matlab for rapid visualization and written to (re) create shapes that could easily be described using geometry. For example, the wind tunnel surface roughness configuration objects constructed by Dong et al. (2002) were recreated by approximating them using oblate spheroids with the equation x2 / a2 + y2 / b2 + z2 / b2 = 1 where x, y and z are the dimensions of a cartesian co-ordinate system and a and b are the semi-major (height) and semi-minor (diameter) axes of an ellipsoid. The shapes were buried in the plane up to their semi-major axis length and their centres were placed with an equidistant spacing that matched the coverages used by Dong et al. (2002). To reconstruct the surface roughness configurations of Marshall's (1971) wind tunnel study we created digital elevation surfaces of portions which were well represented using truncated hemispheroids in the digital reconstruction. 4.1. Isotropic wind tunnel roughness For each simulated surface the number of objects was kept constant and 8 whole shapes were placed on the plane surface. Since spacing, and hence coverage of the objects varied between simulated surfaces, the simulation area also varied (Fig. 2). Each surface was illuminated for a single ϕ and with many i. An example of one half of the configuration surface of hemispheroids for the same i is provided in Fig. 2a and b. Only the downwind half of the DEM surface was used in the calculation of the shadow area to avoid the edge effects of

Table 1 Characteristics of the surfaces used to test the shadow model.

Marshall (1971) Cylinders Dong et al. (2002) Hemispheroids Van Pelt unpublished Natural surface

Layout

Lateral cover (λ)

Object height (mm)

Object diameter (mm)

Average height (mm)

Street

0.0002–0.22

25.4

12.7–127

0.002–11

Diagonal

0.01–0.61

19–65

12–43

0.07–19

N/A

N/A

N/A

N/A

0.0008

N/A = not applicable because the surface is natural.

Wind velocity (m s− 1) At 0.6 m 20.3 At 0.6 m 4–22 At 2 m 0.7–16

Wind direction (°)

Fixed Fixed 1–358

A. Chappell et al. / Remote Sensing of Environment 114 (2010) 1462–1470

1467

Table 2 Aerodynamic roughness lengths (z0) measured in a field for which 1 m2 was used to predict geometric anisotropic z0 using the new model. Wind direction (azimuth angle)

Number of data

z0 median

z0 half interquartile range

Predicted z0

z0 RMSE

mm

mm

mm

mm

16 129 152 182 285 183 46 28 4786

0.058 3.641 0.529 0.033 0.053 0.254 0.043 0.012 0.079

0.500 1.981 0.566 0.032 0.311 0.938 0.038 0.004 0.459

0.540 0.516 0.522 0.477 0.530 0.529 0.542 0.534 0.530

0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.002

± 5° 0 45 90 135 180 225 270 315 All data

absent upstream objects that would contribute shadow to the surface (cf., Chappell & Heritage, 2007). 4.2. Geometric anisotropic natural roughness There is a dearth of available data that combines digital elevation models (DEMs) of natural surfaces with wind speed data and aerodynamic roughness length (z0) observations. In their absence we used a DEM taken from a circular field from which wind velocity measurements were available at its centre (Section 3.2). The field was prepared consistently with machinery to produce a minimal roughness (approx. 1 cm) to satisfy a working hypothesis of that study: there was no discernible preferential orientation of surface roughness to influence aerodynamic resistance. Prior to making wind velocity measurements at the site, the DEM was produced using a laser on a horizontal frame that measured the elevation relative to an arbitrary datum over an area of 1 m2 with a horizontal resolution of 5 mm. The natural surface elevation data were filtered to remove a few spurious measurements. The DEM of the surface is shown in Fig. 3 using examples overlayed of the proportion of illuminated surface (viewed at nadir) that represents four azimuth directions (0°, 45°, 90° and 135°). The shadow overlay is described below (Section 8). Note that an illumination zenith angle of 75° was used and therefore Fig. 3 also represents the variability of erodible area with azimuth angle.

Fig. 5. The sill variance (c) divided by the exponent (α) parameter values of the Gaussian model fitted to the illuminated surface roughness of wind tunnel studies by Dong et al. (2002) and Marshall (1971) plotted against the frontal area index (λ; roughness density).

reconstructed wind tunnel studies of Dong et al. (2002) and Marshall (1971) due to the nature of the objects (discussed below). This relationship is sufficient to enable predictions of z0 / h for other surface roughness configurations. However, the relationship is restricted to this particular model and its parameter values. To investigate a more general approach involving other appropriate models (kernel) a relationship is sought initially between the single scattering albedo (SSA) and z0 / h (Fig. 7). There are strong separate relationships between the SSA and z0 / h for the reconstructed surface roughness of the two wind tunnel studies. The separation between the studies is similar to that evident in the model parameter values (Fig. 6) and is due to the bluff wall nature of the cylinders used in Marshall's study as opposed to the hemispheroids used by Dong et al. (2002). The former does not enable mutual shadowing whilst the latter does. Consequently, when the spacing between objects is smaller than the projected area of the objects (cast shadow) the mutual shadowing of the hemispheroids reduces the SSA but the vertical

5. Evaluation of the relationship between shadow and aerodynamic roughness length (z0) Perhaps the most fundamental assumption made in the current understanding of the relationship between roughness configuration and aerodynamic resistance is that frontal area index (roughness density; λ) adequately captures the characteristics. This is evident in Fig. 4 which shows the relationship between λ and aerodynamic roughness length (z0) standardized by object height (h) using the data from Marshall (1971) and Dong et al. (2002). Using the same surface configurations for both studies, a range of illumination conditions was simulated using ray-casting and these data were fitted with Eq. (2) and then integrated over all illumination zenith angles using Eq. (3). All model fits achieved an RMSE value of less than 0.0045 of shadow proportion and consequently the model parameters were accepted. For brevity, the results of the model fitting and the integrated values are not shown. The values of the variance model parameter (c) were divided by the values of the exponent (α) and are plotted against λ to demonstrate the relationship between λ and model parameters representing illumination/shadow (Fig. 5). Since there is a linear relationship between λ and the illumination function parameters, the relationship between the same parameter values and the aerodynamic roughness length (z0) standardized by object height (h) (Fig. 6) is very similar to the relationship between λ and z0 / h (Fig. 4). The results show a small separation between the

Fig. 6. The sill variance (c) divided by the exponent (α) parameter values of the Gaussian model fitted to the illuminated surface roughness of wind tunnel studies by Dong et al. (2002) and Marshall (1971) plotted against the aerodynamic roughness length (z0) standardized by object height (h).

1468

A. Chappell et al. / Remote Sensing of Environment 114 (2010) 1462–1470

6. Geometric anisotropic aerodynamic resistance (z0) of a heterogeneous surface

wall of the cylinders does not allow it and SSA reduces at a slow rate. When the surface roughness configuration is widely spaced the values of SSA for both studies are very similar and their different geometry makes little difference. Natural particles, aggregates of particles and non-erodible roughness elements (e.g., reg deflation surface) tend towards hemispheroids because of the inherent strength of the shape and because of the weathering/abrasion processes. The hemispheroids represent a more realistic behaviour in mutual shadowing than the cylinders. Hemispheroids have also been found by other workers to adequately represent the light scattering behaviour of a range of unvegetated (Cierniewski, 1987) and vegetated surfaces (Li & Strahler, 1992). This provides a compelling basis to adopt Dong et al. (2002) results using hemispheroids to establish a general relationship between SSA and z0 / h (Fig. 8). A positive variant of the Gaussian model (Eq. (2)) is shown in Fig. 8 (c = 0.2116, α = 10.9940, r = 0.8892) and has a RMSE= 0.0026 for z0 / h. It is this model which is used in subsequent predictions described below.

The final stage of testing the shadow model is to make predictions of z0 over the heterogeneous surface with changing wind directions. Recall that a digital elevation model (DEM) was measured over a 1 m2 area of a 100 m radius circular field prepared with minimal roughness. The DEM has already been presented (Fig. 3) but the shadow overlay has not been described. The pattern of shadow is created by an illumination azimuth angle from an arbitrary origin (0°) that represents the wind direction and an illumination zenith angle of 75° to represent the abrasion angle. The most noticeable aspect of the shadow overlays (Fig. 3) is the large proportion of the DEM that they cover. More than 70% of the surface is sheltered with a wind direction at 0°. This finding has important implications for the estimation of erodible area in dust emission models (see Section 7). Secondly, differences amongst the shadow overlays include the changing patterns of shadow and changing coverage of shadow with different azimuth angle. These changes are responses to the variability in the surface elevation and the systematic decrease in elevation towards the origin of the DEM surface. The patterns evident in Fig. 3 can be generalized as a shadow function using illumination zenith angles and illumination azimuth angles (Fig. 9). In this case, the shadow function is the proportion of shadow visible at nadir; when the shadow function is close to 1 most of the surface is in shadow and when close to 0 most of the surface is directly illuminated. Each curve shows an increase in the illumination function with increasing illumination zenith angle. In general, the curves show only small differences in their shape as a consequence of the prepared surface roughness. Nevertheless, the angular anisotropy of the roughness remains evident in some curves. For a given azimuth angle (e.g., 0°) and its opposing angle (180°) the rate of change is different and is readily evident as a contrast in the magnitude of the shadow at the largest illumination angles. These patterns may reverse as is the case with the remaining azimuth angles. Such complicated angular anisotropic responses are characterised readily by shadow with changing illumination zenith angle offering the potential to model the geometric anisotropic aerodynamic resistance using sensors which measure angular reflectance. Using the previously established isotropic relationship between SSA and z0 (Fig. 8) we predicted the z0 of the heterogeneous surface for a range of wind directions (Table 2). Since the prediction requires a value for object height (h) this was estimated when the average

Fig. 8. Gaussian model fitted to the single scattering albedo (SSA) of Dong et al. (2002) reconstructed wind tunnel surface roughness and its relationship with the measured aerodynamic roughness length (z0) standardized by object height (h).

Fig. 9. Shadow function for several illumination zenith (i) and azimuth (ϕ) angles representing the wind directions over the natural surface digital elevation model shown in Fig. 3.

Fig. 7. Single scattering albedo (SSA) of the reconstructed wind tunnel studies hemispheroid surface roughness and its relationship with the measured aerodynamic roughness length (z0) standardized by object height (h).

A. Chappell et al. / Remote Sensing of Environment 114 (2010) 1462–1470

predicted z0 was the same as the median measured z0 (0.056 mm) for all observations and resulted in a value for h = 2.59 mm. This h value is consistent with the small surface roughness of the prepared surface. Predicted z0 changed only slightly with wind direction. This pattern is also consistent with the prepared nature of the soil surface. However, it is clear that some anisotropy remains since the differences in predicted z0 exceed the uncertainty of the model estimate. Only those measured z0 with approximately the same wind direction (±5°) were used to make predictions. The observed z0 values were highly skewed (not shown) and consequently parametric statistics were avoided. Instead, the median z0 and half the z0 interquartile range are presented here to provide uncertainty about the estimate (Table 2). The measured z0 values are highly variable and some measurement ranges do not contain the values of the predicted z0. The interpretation of these data is difficult because field measurements of z0 are notoriously unreliable (Lancaster et al., 1991; Bauer et al., 1992). Nevertheless, there is good general agreement between the measurements and predictions. 7. Erodible area used in dust emission models The proportion of a surface sheltered from abrasion by saltating material is not included in regional and global dust emission estimates despite for example, a measure (the cumulative shelter angle distribution) being included in the Wind Erosion Prediction System. Dust emission models assume that the relative contribution to the total flux of each size range (e.g., Marticorena and Bergametti, 1995; p. 16,422) is “… proportional to the relative [area] it occupies on the total surface.” There are several weaknesses with this assumption: 1. non-erodible roughness elements protect a. the underlying substrate from the wind and abrasion/sandblasting, b. other non-erodible roughness elements from the wind and hence reduce the momentum extracted by the roughness elements (i.e., z0). 2. the area protected will be a geometric anisotropic function of wind velocity (speed and direction) and configuration of the roughness elements or surface roughness heterogeneity. 3. the spatial organisation of the surface will play a significant role in the sheltering of substrate and non-erodible roughness elements (e.g., ripples). The implications for dust emission models of the geometric anisotropic sheltering of non-erodible surface roughness are illustrated using the model of Marticorena and Bergametti (1995). The vertically integrated horizontal flux (Fh) was calculated to consider the way in which the exposed erodible surface contributes dust emission. Two configurations of the same hemispheroid (38 mm high and 24 mm diameter) used in Dong's et al. (2002) wind tunnel study were included with different spacings (Table 3). The obstacle-free surface was estimated following the approach of most dust emission models and the alternative erodible area was calculated using the zenith abrasion angle (75°) to represent exposure to sandblasting. The Fh used a single lognormal mode at 100 μm and σ=2, measured z0 and a smooth roughness z0s ≈0.026 mm, following the wind tunnel measurements of Dong et al. (2002).

1469

Table 3 shows that there are large differences in z0 and U* for each roughness configuration caused by the spacing of the hemispheroids. In the case of widely-spaced hemispheroids the area not covered was 95% of the surface area but using the abrasion angle only approximately 79% was erodible. Consequently, estimates of Fh for the area not covered were 17% larger than the erodible area. In the case of closely-spaced hemispheroids the area not covered was 25% but the erodible area was approximately 4% of the surface area. Using the same z0 and U* values the Fh for the erodible area was 85% smaller than that of the area not covered (Table 3). These results show that substantial over-estimates of Fh and hence dust emission can arise from the use of area not covered to estimate the erodible area. The zenith abrasion angle and sheltered area offers an alternative rapid and consistent estimate of the erodible area for use in dust emission models.

8. Conclusions Predictions of aerodynamic resistance were made using the relationship between single scattering albedo (SSA) and aerodynamic roughness length (z0) provided by wind velocity profile measurements in a wind tunnel study. This new model provided the ability to account for the anisotropy of aerodynamic resistance caused by changing wind direction, which is particularly important over heterogeneous surfaces. These achievements were possible because the model used illumination and shadow to represent the geometric projection of frontal area index in any orientation. The estimation of shadow cast at a 75° illumination zenith angle was equivalent to the sheltering from abrasion by saltating particles. The large area sheltered from abrasion was in contrast to the area not covered used in current dust emission models. Consequently, the findings above suggested strongly that wind erosion and dust emission models that do not account for geometric anisotropic sheltering of the surface may considerably over-estimate the horizontal flux and hence dust emission. The estimation of z0 using measured wind velocity profiles is notoriously unreliable particularly over large areas or multiple dunes. That unreliability would perhaps be more evident if more studies included the uncertainty associated with their estimates. Notwithstanding that omission, much of the uncertainty in z0 field estimates arise from the difficulty in accounting for the inevitable geometric anisotropy in z0 as a consequence of varying wind directions over heterogeneous surfaces. Consequently, field z0 is likely to be neither precise nor accurate. The results of this study suggest that the new model provides an alternative measure that is suitable to investigate geometric anisotropic aerodynamic resistance across scales of variation (e.g., grain, ripple, dune etc.).

Acknowledgements The first author is grateful to M. Ekström and N. Chappell for their enduring support and critical observations. He is grateful to P. Hairsine for his incisive comments on an early version of this manuscript. The authors would like to thank B. Marticorena and G. Bergametti for making the dust emission code available and the two anonymous referees who provided constructive comments on the manuscript.

Table 3 Vertically integrated horizontal flux (Fh) for two configurations of hemispheres (38 mm height and 24 mm diameter) exposed to simulated sandblasting using the proportion of the wind tunnel not covered and the erodible proportion using the zenith abrasion angle (see text for details). Not covered

Erodible proportion

Lateral cover

Measured z0

U*

Area

Fh

Area

Fh

Difference Fh

λ

mm

mm s−1

%

g mm−1 s−1

%

g mm−1 s− 1

%

0.04 0.60

0.65 0.21

884.0 707.3

95.0 25.0

4.0 0.9

78.8 3.7

3.4 0.1

17 85

1470

A. Chappell et al. / Remote Sensing of Environment 114 (2010) 1462–1470

References Alfaro, S. C., & Gomez, L. (1995). Improving the large-scale modelling of the saltation flux of soil particles in the presence of non-erodible elements. Journal of Geophysical Research, 100, 16357−16366. Arya, S. P. S. (1975). A drag partition theory for determining the large-scale roughness parameter and wind stress on the Arctic Pack Ice. Journal of Geophysical Research, 80, 3447−3454. Bauer, B. O., Sherman, D. J., & Wolcott, J. F. (1992). Sources of uncertainty in shear stress and roughness length estimates derived from velocity profiles. Professional Geographer, 44 (4), 453−464. Brown, S., Nickling, W. G., & Gillies, J. A. (2008). A wind tunnel examination of shear stress partitioning for an assortment of surface roughness distributions. Journal of Geophysical Research. doi:10.1029/2007JF000790. Callot, Y., Marticorena, B., & Bergametti, G. (2000). Geomorphologic approach for modeling the surface features over arid environments in a model of dust emissions: Application to the Sahara desert. Geodinamica Acta, 13, 245−270. Chappell, A., & Heritage, G. L. (2007). Using illumination and shadow to model aerodynamic resistance and flow separation: An isotropic study. Atmospheric Environment, 41(28), 5817−5830. Chappell, A., Strong, C., McTainsh, G., & Leys, J. (2007). Detecting induced in situ erodibility of a dust-producing playa in Australia using a bi-directional soil spectral reflectance model. Remote Sensing of Environment, 106, 508−524. Chappell, A., Warren, A., O'Donoghue, A., Robinson, A., Thomas, A., & Bristow, C. (2008). The implications for dust emission modeling of spatial and vertical variations in horizontal dust flux and particle size in the Bodélé Depression, Northern Chad. Journal of Geophysical Research, 113, D04214. doi:10.1029/2007JD009032. Chappell, A., Zobeck, T. M., & Brunner, G. (2006). Using bi-directional soil spectral reflectance to model soil surface changes induced by rainfall and wind-tunnel abrasion. Remote Sensing of Environment, 102(3–4), 328−343. Cierniewski, J. (1987). A model for soil surface roughness influence on the spectral response of bare soils in the visible and near-infrared range. Remote Sensing of Environment, 23, 97−115. Dong, Z., Liu, X., & Wang, X. (2002). Aerodynamic roughness of gravel surfaces. Geomorphology, 43(1–2), 17−31. Hagen, L. J. (1990). A wind erosion prediction system to meet user needs. Journal of Soil and Water Conservation, 46, 106−111. Hapke, B. (1993). Theory of reflectance and emittance spectroscopy. Cambridge, UK: Cambridge Univ. Press. Jacquemoud, S., Bater, F., & Hanocq, J. F. (1992). Modeling spectral and bidirectional soil reflectance. Remote Sensing of Environment, 41, 123−132. Lancaster, N., Greeley, R., & Rasmussen, K. R. (1991). Interaction between unvegetated desert surfaces and atmospheric boundary layer: A preliminary assessment. Acta Mechanica. Supplement, 2, 89−102. Li, X., & Strahler, A. H. (1992). Geometric-optical bidirectional reflectance modeling of the discrete crown vegetation canopy: Effect of crown shape and mutual shadowing. IEEE Transactions on Geoscience and Remote Sensing, 30(2), 276−291. Liang, S., & Townshend, J. R. G. (1996). A modified Hapke model for soil bidirectional reflectance. Remote Sensing of Environment, 55, 1−10. Mahowald, N., Bryant, R. G., del Corral, J., & Steinberger, L. (2003). Ephemeral lakes and desert dust sources. Geophysical Research Letters, 30(2), 1074. doi:10.1029/2002GL016041. Marshall, J. K. (1971). Drag measurements in roughness arrays of varying density and distribution. Agricultural Meteorology, 8, 269−292. Marticorena, B., & Bergametti, G. (1995). Modeling the atmospheric dust cycle: 1. Design of a soil-derived dust emission scheme. Journal of Geophysical Research, 100 (D8), 16415−16430.

Marticorena, B., Bergametti, G., Aumont, B., Callot, Y., N'doume´, C., & Legrand, M. (1997). Modeling the atmospheric dust cycle: 2 — Simulation of Saharan sources. Journal of Geophysical Research, 102, 4387−4404. Marticorena, B., Chazette, P., Bergametti, G., Dulac, F., & Legrand, M. (2004). Mapping the aerodynamic roughness length of desert surfaces from the POLDER/ADEOS bidirectional reflectance product. International Journal of Remote Sensing, 25, 603−626. Pinty, B., Verstraete, M. M., & Dickinson, R. E. (1989). A physical model for predicting bidirectional reflectances over bare soil. Remote Sensing of Environment, 27, 273−288. Potter, K. N., Zobeck, T. M., & Hagen, L. J. (1990). A microrelief index to estimate soil erodibility by wind. Trans. ASAE, 33, 151−155. Raupach, M. R. (1992). Drag and drag partition on rough surfaces. Boundary-Layer Meteorology, 60, 374−396. Raupach, M. R., & Lu, H. (2004). Representation of land-surface processes in aeolian transport models. Environmental Modelling & Software, 19, 93−112. Raupach, M. R., Gillette, D. A., & Leys, J. F. (1993). The effect of roughness elements on wind erosion thresholds. Journal of Geophysical Research, 98, 3023−3029. Raupach, M. R., Hughes, D. E., & Cleugh, H. A. (2006). Momentum absorption in roughwall boundary layers with sparse roughness elements in random and clustered distributions. Boundary-Layer Meteorology, 120, 201−218. Roujean, J. L., Leroy, M., & Deschamps, P. Y. (1992). A bidirectional reflectance model of the Earth's surface for the correction of remote sensing data. Journal of Geophysical Research, 97, 20455−20468. Schlichting, H. (1936). Experimental investigation of the problem of surface roughness. Ingenieur-Archiv, 7, 1−34 680. Shao, Y. (2001). A model for mineral dust emission. Journal of Geophysical Research, 106 (D17), 20239−20254. Shao, Y., & Yang, Y. (2005). A scheme for drag partition over rough surfaces. Atmospheric Environment, 39, 7351−7361. Shao, Y., Raupach, M. R., & Leys, J. F. (1996). A model for predicting aeolian sand drift and dust entrainment on scales from paddock to region. Australian Journal of Soil Research, 34, 309−342. Sorensen, M. (1985). Estimation of some aeolian saltation transport parameters from transport rate profiles. Proceeding of International Workshop on the Physics of Blown Sand. Aarhus, May 28–31 (pp. 141−190). Denmark: University of Aarhus. Tegen, I. (2003). Modeling the mineral dust aerosol cycle in the climate system. Quaternary Science Reviews, 22, 1821−1834. Tegen, I., Heinold, B., Todd, M. C., Helmert, J., Washington, R., & Dubovik, O. (2006). Modelling soil dust aerosol in the Bodélé Depression during the BoDEx campaign. Atmospheric Chemistry and Physics, 6, 4345−4359. Tegen, I., Harrison, S. P., Kohfeld, K. E., Prentice, I. C., Coe, M. C., & Heimann, M. (2002). The impact of vegetation and preferential source areas on global dust aerosol: Results from a model study. Journal of Geophysical Research, 107. doi:10.1029/2001JD000963. Warren, A., Chappell, A., Todd, M. C., Bristow, C., Drake, N., Engelstaedter, S., et al. (2007). Dust-raising in the dustiest place on earth. Geomorphology, 92(1–2), 25−37. Webster, R., & Oliver, M. A. (2001). Geostatistics for environmental scientists.Statistics in Practice Series: Wiley 271 pp. Wu, Y., Gong, P., Liu, Q., & Chappell, A. (2008). Retrieving photometric properties of desert surfaces in China using the Hapke model and MISR data. Remote Sensing of Environment. doi:10.1016/j.rse.2008.09.006. Zobeck, & Popham (1998). Wind erosion roughness index response to observation spacing and measurement distance. Soil & Tillage Research, 45, 311−324.