LETTERS PUBLISHED ONLINE: 5 FEBRUARY 2012 | DOI: 10.1038/NGEO1381

Internal boundary layer model for the evolution of desert dune fields Douglas J. Jerolmack1 *, Ryan C. Ewing2 , Federico Falcini1 , Raleigh L. Martin1 , Claire Masteller1 , Colin Phillips1 , Meredith D. Reitz3 and Ilya Buynevich4 Desert dunes often exhibit remarkable changes in their morphology over short distances. For example, sediment-rich dunes can break up into smaller, isolated features, and then become stabilized by plants, over distances of kilometres1–6 . These pattern transitions often coincide with spatial variations in sediment supply1,3,5 , transport rate6,7 , hydrology8 and vegetation9–11 , but these factors have not been linked mechanistically. Here we hypothesize that the abrupt increase in roughness at the upwind margins of dune fields triggers the development of an internal boundary layer12–18 that thickens downwind and causes a spatial decrease in the surface wind stress. We demonstrate that this mechanism forces a downwind decline in sand flux at White Sands, New Mexico, using a combination of physical theory14–19 , repeated airborne altimetry surveys and field observations. The declining sand flux triggers an abrupt increase in vegetation density, which in turn leads to changes in groundwater depth and salinity—showing that aerodynamics, sediment transport and ecohydrology are tightly interconnected in this landscape. We conclude that, despite the documented complex climatic and geologic history of White Sands20 , internal boundary layer theory explains many of the observed first-order patterns of the dune field. The gypsum dunes of White Sands National Monument were created by Holocene deflation of evaporite deposits from Pleistocene Lake Otero, and align transverse to prevailing winds from the southwest4,20–22 . The dune field abruptly emerges from the upwind deflation plain, the Alkali Flat, to form 5-m-high crescentic dunes, within a few hundred metres (Fig. 1). The large sand ridge at the upwind margin serves as the primary source of sediment to the dune field, and transported sand is trapped within the dunes19,20 , making White Sands essentially a closed sedimentary system. Within the next 8 km downwind, the dune-field pattern evolves from crescentic, to barchan, to vegetated parabolic dunes before terminating abruptly2,6,8,19,20 . Active dunes become more sparse in their areal coverage downwind, indicating a decline in mobile sand availability. Despite arid conditions, the depth to the groundwater table is only h ∼ 1 m below the interdune surface4,8,23 . Recent work using repeated, high-resolution topographic surveys has shown that annually averaged sand flux, hqs i, decreases rapidly in the first 1–2 km downwind of the ridge, and then more slowly for distances up to 10 km (ref. 11). A similar pattern has been observed in other dune fields3,7 . The Alkali-Flat to dune-field transition represents an abrupt change in surface roughness, z0 (smooth → rough). Mesoscale

meteorological studies of boundary layer flow adjustments to such transitions have been undertaken in coastal environments12,13 , but not, to our knowledge, in inland settings such as White Sands. Field studies confirm theoretical predictions15,16,18 and small-scale experiments1,14,17,24 : at the smooth → rough transition there is a sudden increase in the boundary shear stress, τb , exerted by the wind owing to flow convergence and turbulent shear. This is followed by a relaxation of τb toward a new equilibrium value, as the initial momentum perturbation advects downwind and diffuses into the atmosphere above. The rate of thickening of the zone of perturbed flow—the internal boundary layer (IBL)— determines the relaxation rate of τb (Fig. 2). As sand transport is nonlinearly related to τb (see below), IBL development should drive significant spatial variations in sand flux across a dune field. Here we test the idea that growth of an IBL, triggered by dune roughness, is responsible for the observed decline in sand flux downwind of the sediment ridge at White Sands. We consider the idealized case of two-dimensional, steady turbulent flow, a thermally neutral atmospheric boundary layer (ABL) and a large step increase in surface roughness. We verified that these assumptions are consistent with the large-time- and spaceaveraged behaviour of sand transporting winds at White Sands (Supplementary Information). The thickness of the IBL (δi ) is defined by the elevation above the surface at which flow velocity is equal to the free-stream velocity, U . Growth of δi has been found to follow the form16  0.8 x δi = 0.28 (1) z02 z02 where x is the distance downwind of the step and z02 is the downwind roughness (dune) length scale. To derive the downwind profile τb (x), a relation between τb and δi is needed. An empirical expression derived from pipe flow experiments25 is   ν 1/4 τb = 0.0225ρU 2 (2) U δi where for our case ρ = 1.23 kg m−3 and ν = 0.32 × 10−4 m2 s−1 are the density and kinematic viscosity of air, respectively, at 20 ◦ C. This equation is applicable at distances x/z02  1 after the step change, and for length scales larger than that of an individual dune (that is, spatial variation of τb over individual dunes is not considered). From atmospheric soundings, we determined the freestream velocity associated with the most effective wind—the wind whose combination of magnitude and frequency of occurrence

1 Department

of Earth and Environmental Science, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA, 2 Department of Geological Sciences, University of Alabama, Tuscaloosa, Alabama 35487, USA, 3 Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA, 4 Department of Earth and Environmental Science, Temple University, Philadelphia, Pennsylvania 19122, USA. *e-mail: [email protected]. NATURE GEOSCIENCE | ADVANCE ONLINE PUBLICATION | www.nature.com/naturegeoscience © 2012 Macmillan Publishers Limited. All rights reserved.

1

NATURE GEOSCIENCE DOI: 10.1038/NGEO1381

LETTERS a

35

m 0

0

1

2

3

4 5 6 Distance down dune field (km)

7

8

9

b

Figure 1 | White Sands dune field. a, Lidar-derived topographic map from June 2007. Dune migration is from left to right. Note the abrupt appearance of dunes as a topographic ridge at the upwind (left) margin (x = 0 km). The barchan–parabolic transition begins around x = 7–8 km; the red box shows the area in b. b, Aerial image of the barchan–parabolic transition from Google Earth. Dunes have the following characteristic scales: crest-to-crest spacing, λ ∼ 150 m; height, H ∼ 3 m; migration speed, c ∼ 2 m yr−1 and average migration direction 25◦ N of east.

hqs i = 46x −3/10

(3) 2

δi

h B¬P transition

x

b

3 2 1 0 ¬4

0

4

8

12

x (km)

−1

where x has units of metres and hqs i has units of m yr . Comparing equation (3) with the flux profile estimated from topographic data11 , the agreement is remarkable (Fig. 3). Importantly, the analytical expression was not fit to sand flux data; all aerodynamic parameters were independently estimated (Supplementary Information). Thus, the large-scale spatial pattern of hqs (x)i can be explained completely by the growth of an IBL, owing to an abrupt change in roughness resulting from the dunes. The asymptotic value of hqs (x)i also agrees well with the predicted annual sand flux computed using a meteorologic station just downwind of the dune field (Supplementary Information; Fig. 3), providing independent validation of equation (3). Spatial changes in hqs i drive erosion and deposition through conservation of mass. IBL theory predicts intense scour to occur at the upwind margin of the dune field, and deposition to occur downwind of the sediment ridge (Fig. 2). The resulting transport pattern is like a bulldozer, causing migration of the entire dune field in the dominant wind direction. This pattern is consistent with geologic evidence for dune migration at White Sands4,6,21,22 . 2

a

τ b(x)/τ b∞

moves the most sand—to be U = 20 m s−1 ; from light detection and ranging (lidar) data we computed z02 = 10−2 m (Supplementary Information). Combining equations (1) and (2), we predict that the boundary stress downwind of the ridge at White Sands should decay as τb = 0.41x −1/5 , where stress has units of pascals and distance has units of metres. Applying the Bagnold saltation equation for instantaneous sand flux (qs ) of the form qs ∼ τb3/2 (refs 1,19), and taking into account the intermittency in magnitude and direction of the winds (Supplementary Information), we obtain an analytical expression for the expected decline of annual sand flux with distance downwind of the sediment ridge:

Figure 2 | Definition sketch and boundary stress profile. a, Schematic representation of IBL, vegetation and groundwater depth across White Sands dune field. The barchan–parabolic (B–P) transition is shown, and green symbols represent the relative plant density. All variables are defined in the text. b, Boundary stress profile, normalized by the asymptotic downwind value, across the roughness transition. The solid line is the theoretical profile from equations (1) and (2) for x > 0, valid where x/z02  1; the dashed line is the qualitative trend generalized from numerical and experimental studies14–18 . Red (blue) is the zone of expected erosion (deposition); the colour intensity reflects the magnitude.

More quantitatively, our predicted theoretical deposition rate for the middle of the dune field is within the range of measured accumulation rates for that location determined from optically stimulated luminescence4 (Fig. 3). Although hqs (x)i declines gradually across the dune field, vegetation density changes abruptly at a distance x = 7–8 km downwind of the sediment ridge (Figs 1 and 3). Dunes transition

NATURE GEOSCIENCE | ADVANCE ONLINE PUBLICATION | www.nature.com/naturegeoscience © 2012 Macmillan Publishers Limited. All rights reserved.

NATURE GEOSCIENCE DOI: 10.1038/NGEO1381 a

LETTERS

10

0

4 2

h (m)

2

4

6 x (km)

8

10

8

0.3

6 4 h, GPR h, piezometers h, wells Plants

0.9 1.2

0

2

4

6 x (km)

2

8

10

1.00

0 12

0

1

2 3 E (mm d¬1)

4

5

Figure 4 | Conceptual model of evaporation rate as a function of groundwater depth. E = E0 until Lc , below which E declines rapidly with h (after ref. 28). Inferred bistable states are shown: the blue circle corresponds to the barchan area with no vegetation and high E, whereas the green circle corresponds to the parabolic area with high plant density and low E (Supplementary Information).

Figure 3 | Sand flux, groundwater and vegetation patterns. a, Annual sand-flux profile downwind of the ridge: circles, lidar-derived values; red curve, equation (3) with R2 = 0.69; dashed red line, annual sand flux hqs i = 2.7 m2 yr−1 predicted from wind records. Grey line, predicted deposition rate from mass conservation; triangle (error bars), mean (range) of measured deposition rates4 . b, Groundwater depth and plant density profiles; rapid changes in both at x = 8 km correspond to the barchan–parabolic transition. Depth was estimated using three methods—see text and Supplementary Information. Lines are to guide the eye.

from unvegetated barchans to heavily vegetated parabolics over a distance of only a few wavelengths. Ref. 11 confirmed model predictions9 that this pattern change occurs when dune migration rate drops below a threshold value such that vegetation can colonize. Because plant transpiration will affect the water balance, we anticipate a change in the groundwater depth across the barchan–parabolic transition. The average (equilibrium) depth of the groundwater table, hE , should be determined by the depth at which water losses balance gains26,27 : E(h) + T = P + qw

0.75

12

0

0.6

0.50

Plant density (number per dune)

b

0

Lc

0.25

h (m)

〈

〈qs (m2 yr¬1)

6

0

E0

Deposition rate (mm yr¬1)

8

(4)

where the variables represent annually averaged rates of evaporation, transpiration, precipitation and groundwater efflux, respectively. The only depth-dependent variable in equation (4) is evaporation rate, which should be approximately equal to the surface value (E0 ) for h ≤ Lc , and decline rapidly with h for h > Lc , where Lc is the capillary length determined by the balance between capillary and hydrostatic pressure28 . E0 ≈ 5 mm d−1 at White Sands, which far exceeds inputs (P ≈ 0.5 mm d−1 , qw < 0.1 mm d−1 ); thus, we expect that the groundwater table is drawn down to a depth hE = Lc + ε, where ε is some distance determined by E = f (h) (Fig. 4). We estimate Lc ≈ 0.3 m for White Sands (Supplementary Information). In unvegetated barchan interdunes where T = 0, evaporation rates are expected to be high and thus hE may approach the limit Lc . In densely vegetated parabolic interdunes, transpiration is no longer negligible. The increase in T would draw down the water table until E(h) is small enough that the balance (equation (4)) is again satisfied (Fig. 4).

To test these ideas, we conducted a survey of groundwater depth along the downwind profile at White Sands using groundpenetrating radar (GPR) and observation wells, supplemented with previously installed piezometers (Supplementary Information). Data show that h is approximately equal to 0.5 m across the barchan section of the dune field, consistent with our expectation that h ≈ Lc in this region. Groundwater depth then drops across the transition, to an average depth h ≈ 1 m in the parabolics (Fig. 3). We infer that vegetation draws down the water table by about half a metre, within the range of values found in the literature27,29 . Results are in agreement with a less detailed survey conducted previously8 , which also showed that the salinity of the groundwater decreased across the transition. A recent theory27 predicted that two stable states can exist in shallow, saline groundwater systems as a result of feedbacks between salt accumulation and vegetation dynamics: (1) a non-vegetated state with a shallow and salty groundwater table and (2) a vegetated state with a deeper and fresher groundwater table. It seems that both of these states are present at White Sands, and that a threshold in landscape stability forces a shift from (1) to (2). Because the vegetation transition is caused by declining sand flux, the developing IBL facilitates this transition. As none of these theoretical arguments are unique to White Sands, it is likely that this interplay among IBL, sediment transport and (eco)hydrologic dynamics also occurs in other dune fields. Our findings help to illuminate the evolution of the enigmatic White Sands dune field. There is strong evidence that past and ongoing deflation of Lake Otero sediments—which produced the modern dune field—is the result of falling groundwater levels over the past ∼7 kyr due to increased aridity4,20 . However, there is also substantial sediment accumulation within the dune field itself, which has been attributed to a rising groundwater table4,22 . Our work provides an explanation for both of these observations. We propose that climate-driven deflation created the dune field, whereas aerodynamics maintains a zone of deposition within the dunes and erosion at the upwind margin. This transient accumulation zone migrates with velocity of order 100 m yr−1 east-northeast, that is, the group velocity of the dune field. Consider the characteristic timescales of the system: dune migration (wavelength/migration rate), 102 yr; vegetation growth, 101 yr (ref. 11); groundwater response, 10−3 yr (Supplementary Information); aerodynamic response, 10−7 yr (seconds). The rate of landscape change is slow enough that it may be considered quasisteady, in which case the plant, vegetation and air-flow patterns are in equilibrium with the (changing) configuration of the dune field. One interesting implication is that the barchan–parabolic transition—and the

NATURE GEOSCIENCE | ADVANCE ONLINE PUBLICATION | www.nature.com/naturegeoscience © 2012 Macmillan Publishers Limited. All rights reserved.

3

NATURE GEOSCIENCE DOI: 10.1038/NGEO1381

LETTERS accompanying change in groundwater depth and salinity—should occur at a fixed distance from, and should migrate with, the upwind margin of the dune field. Another is that—except for very steep slopes close to palaeo-lake shorelines—the groundwater table depth is slaved to the land surface position, occurring at a fixed depth below the surface that is determined by evaporation rates and vegetation. It is not known, however, how stable this equilibrium is to external perturbations. Proposed groundwater withdrawals in the area, or anthropogenic climate change30 , could drive a rapid change in groundwater depth that would desiccate the surface and facilitate erosion of the currently stable interdunes26 . The IBL model presented here ignores many important details of airflow and sediment transport, for the sake of analytical simplicity and physical insight. That results are robust to such simplifications suggests that this compact model should have broad applicability. Received 16 September 2011; accepted 3 January 2012; published online 5 February 2012

References 1. Bagnold, R. A. The Physics of Blown Sand and Desert Dunes (Methuen, 1941). 2. McKee, E. D. Structures of dunes at White Sands National Monument, New Mexico (and a comparison with structures of dunes from other selected areas). Sedimentology 7, 3–69 (1966). 3. Lancaster, N. Winds and sand movements in the Namib Sand Sea. Earth Surf. Process. Landforms 10, 607–619 (1985). 4. Kocurek, G. et al. White Sands Dune Field, New Mexico: Age, dune dynamics and recent accumulations. Sedim. Geol. 197, 313–331 (2007). 5. Pye, K. & Tsoar, H. Aeolian Sand and Sand Dunes (Springer, 2009). 6. Ewing, R. C. & Kocurek, G. A. Aeolian dune interactions and dune-field pattern formation: White Sands Dune Field, New Mexico. Sedimentology 57, 1199–1219 (2010). 7. Illenberger, W. K. & Ruse, I. C. A sand budget for the Alexandria coastal dunefield, South Africa. Sedimentology 35, 513–521 (1988). 8. Langford, R., Rose, J. & White, D. Groundwater salinity as a control on development of eolian landscape: An example from the White Sands of New Mexico. Geomorphology 105, 39–49 (2009). 9. Duran, O. & Herrmann, H. Vegetation against dune mobility. Phys. Rev. Lett. 97, 188001 (2006). 10. Baas, A. C. W. & Nield, J. M. Modelling vegetated dune landscapes. Geophys. Res. Lett. 34, L06405 (2007). 11. Reitz, M. D., Jerolmack, D. J., Ewing, R. C. & Martin, R. L. Barchan–parabolic dune pattern transition from vegetation stability threshold. Geophys. Res. Lett. 37, L19402 (2010). 12. Hsu, S. A. Boundary-layer meteorological research in the coastal zone. Geosci. Man XVIII, 99–111 (1977). 13. Rasmussen, K. R. Some aspects of flow over coastal dunes. Proc. R. Soc. Edinb. B 96, 129–147 (1989). 14. Bradley, E. F. A micrometeorological study of velocity profiles and surface drag in the region modified by a change in surface roughness. Q. J. R. Meteorol. Soc. 94, 361–379 (1968). 15. Rao, K. S., Wyngaard, J. C. & Coté, O. R. The structure of the two-dimensional internal boundary layer over a sudden change of surface roughness. J. Atmos. Sci. 31, 738–746 (1974).

4

16. Wood, D. H. Internal boundary layer growth following a step change in surface roughness. Bound.-Layer Meteorol. 22, 241–244 (1982). 17. Stevenson, D. C. & Lee, W. L. Proc. 9th Australasian Fluid Mech. Conf. 379–382 (University of Auckland, 1986). 18. Garratt, J. R. The internal boundary layer: A review. Bound.-Layer Meteorol. 50, 171–203 (1990). 19. Jerolmack, D. J., Reitz, M. D. & Martin, R. L. Sorting out abrasion in a gypsum dune field. J. Geophys. Res. 116, F02003 (2011). 20. Langford, R. The Holocene history of the White Sands dune field and influences on eolian deflation and playa lakes. Quat. Int. 104, 31–39 (2003). 21. Fryberger, S. Geological overview of White Sands National Monument Preprint at http://www.nps.gov/archive/whsa/Geology%20of%20White% 20Sands/GeoHome.html (2009). 22. Szynkiewicz, A. et al. Origin of terrestrial gypsum dunes—implications for Martian gypsum-rich dunes of Olympia Undae. Geomorphology 121, 69–83 (2010). 23. Newton, T. 2010 Annual Report to the National Park Service (National Park Service, 2010). 24. Antonia, R. A. & Luxton, R. E. The response of a turbulent boundary layer to a step change in surface roughness. Part 1. Smooth to rough. J. Fluid Mech. 48, 721–761 (1971). 25. Glauert, M. B. The wall jet. J. Fluid Mech. 1, 625–643 (1956). 26. Chen, J. S. et al. Water resources: Groundwater maintains dune landscape. Nature 432, 459–460 (2004). 27. Runyan, C. & D’Odorico, P. Ecohydrological feedbacks between salt accumulation and vegetation dynamics: Role of vegetation–groundwater interactions. Wat. Resour. Res. 46, W11561 (2010). 28. Lehmann, P., Assouline, S. & Or, D. Characteristic lengths affecting evaporative drying of porous media. Phys. Rev. E 77, 056309 (2008). 29. Gribovszki, Z., Szilágyi, J. & Kalicz, P. Diurnal fluctuations in shallow groundwater levels and streamflow rates and their interpretation—a review. J. Hydrol. 385, 371–383 (2010). 30. Maxwell, R. M. & Kollet, S. J. Interdependence of groundwater dynamics and land-energy feedbacks under climate change. Nature Geosci. 1, 665–669 (2008).

Acknowledgements Research partially supported by the National Science Foundation through EAR-PF-0846233 to R.C.E., and EAR-0810038 to D.J.J. We thank K. Litwin, A. Bhattachan, A. Boles and B. McNutt for their assistance in collecting these data. We are especially grateful to D. Bustos of the National Park Service (NPS) for facilitating this study. Lidar data supplied by G. Kocurek, funded by a grant from NPS as part of the Chihuahuan Desert Network Inventory and Monitoring Program. Ideas and suggestions by G. Kocurek and D. Mohrig improved the paper.

Author contributions All authors contributed equally to data collection and analysis. D.J.J. designed the study with input from R.C.E., F.F., R.L.M., C.P. and M.D.R. The paper was written by D.J.J. and edited by all authors. I.B. and C.M. led the groundwater component.

Additional information The authors declare no competing financial interests. Supplementary information accompanies this paper on www.nature.com/naturegeoscience. Reprints and permissions information is available online at http://www.nature.com/reprints. Correspondence and requests for materials should be addressed to D.J.J.

NATURE GEOSCIENCE | ADVANCE ONLINE PUBLICATION | www.nature.com/naturegeoscience © 2012 Macmillan Publishers Limited. All rights reserved.