# Water Vapor Transport Across an Arid Sand Surface-Non-Linear Thermal Coupling, Wind-Driven Pore Advection, Subsurface Waves, and Exchange With the Atmospheric Boundary Layer

M. Y. Louge, Alexandre Valance, J. Xu, Ahmed Ould El Moctar, Patrick Chasle

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

M. Y. Louge, Alexandre Valance, J. Xu, Ahmed Ould El Moctar, Patrick Chasle. Water Vapor Transport Across an Arid Sand Surface-Non-Linear Thermal Coupling, Wind-Driven Pore Advection, Subsurface Waves, and Exchange With the Atmospheric Boundary Layer. Journal of Geophysical Research: Earth Surface, 2022, 127 (4), pp.e2021JF006490. 10.1029/2021JF006490. hal-03657541

HAL Id: hal-03657541

https://hal.science/hal-03657541

Submitted on 14 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Water Vapor Transport Across an Arid Sand Surface-Non-Linear Thermal Coupling, Wind-Driven Pore Advection, Subsurface Waves, and Exchange With the Atmospheric Boundary Layer

M. Y. Louge, Alexandre Valance, J. Xu, Ahmed Ould El Moctar, Patrick Chasle

► To cite this version:

M. Y. Louge, Alexandre Valance, J. Xu, Ahmed Ould El Moctar, Patrick Chasle. Water Vapor Transport Across an Arid Sand Surface-Non-Linear Thermal Coupling, Wind-Driven Pore Advection, Subsurface Waves, and Exchange With the Atmospheric Boundary Layer. Journal of Geophysical Research: Earth Surface, 2022, 127 (4), pp.e2021JF006490. 10.1029/2021JF006490. hal-03657541v2

HAL Id: hal-03657541

https://hal.science/hal-03657541v2

Submitted on 24 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Water vapor transport across an arid sand surface -
non-linear thermal coupling, wind-driven pore
advection, subsurface waves, and exchange with the
atmospheric boundary layer

M. Y. Louge<fnref n="1" />, A. Valance<fnref n="2" />, J. Xu<fnref n="1" />, A. Ould el-Moctar<fnref n="3" />, P. Chasle<fnref n="2" />

Key Points:

* We detect tiny spatio-temporal moisture variations below a hyper-arid mobile dune with a new capacitance instrument.

* We predict moisture profiles, internal waves and wind-driven pore advection. Grain evaporation is a kinetic-limited, activated process.

* The moisture surface flux is weaker than expected, and not always proportional to the mass fraction difference between surface and ambient.

<footnote n="1">¹Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY 14853, USA.</footnote>

<footnote n="2">²Institut de Physique de Rennes, Université de Rennes 1, France.</footnote>

<footnote n="3">³Ecole Polytechnique de l'Université de Nantes, France.</footnote>
# JGR Earth Surface

## RESEARCH ARTICLE

10.1029/2021JF006490

**Key Points:**

* We detect tiny spatiotemporal moisture variations below a hyper-arid mobile dune with a new capacitance instrument

* We predict moisture profiles, internal waves and wind-driven pore advection. Grain evaporation is a kinetic-limited, activated process

* The moisture surface flux is weaker than expected, and not always proportional to the mass fraction difference between surface and ambient

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

M. Y. Louge,
MYL3@cornell.edu

**Citation:**

Louge, M. Y., Valance, A., Xu, J., Ould el-Moctar, A., & Chasle, P. (2022). Water vapor transport across an arid sand surface—Non-linear thermal coupling, wind-driven pore advection, subsurface waves, and exchange with the atmospheric boundary layer. *Journal of Geophysical Research: Earth Surface*, 127, e2021JF006490. https://doi.org/10.1029/2021JF006490

Received 19 OCT 2021
Accepted 12 MAR 2022

**Author Contributions:**

**Conceptualization:** M. Y. Louge, A. Valance, A. Ould el-Moctar

**Data curation:** M. Y. Louge

**Formal analysis:** M. Y. Louge, A. Valance, A. Ould el-Moctar

**Funding acquisition:** M. Y. Louge

**Investigation:** M. Y. Louge, A. Valance, J. Xu, A. Ould el-Moctar

**Methodology:** M. Y. Louge, A. Valance, J. Xu, A. Ould el-Moctar, P. Chasle

© 2022. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

## Water Vapor Transport Across an Arid Sand Surface—
Non-Linear Thermal Coupling, Wind-Driven Pore Advection,
Subsurface Waves, and Exchange With the Atmospheric
Boundary Layer

M. Y. Louge¹, A. Valance², J. Xu¹, A. Ould el-Moctar³, and P. Chasle²

¹Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY, USA, ²Institut de Physique de Rennes, Université de Rennes 1, Rennes, France, ³Université de Nantes, CNRS, Laboratoire de Thermique et Energie, UMR 6607, Nantes Cedex, France

**Abstract** Although vapor exchanged across hyper-arid surfaces without free liquid affects the water budget of sand seas, its mechanism is poorly documented for want of accurate instruments with fine spatial resolution. To rectify this, we report bulk density profiles and spatiotemporal variations of vapor mass fraction just below the surface of a mobile dune, acquired with a multi-sensor capacitance probe sensitive to tiny water films adsorbed on sand grains. We also record wind speed and direction, ambient temperature and relative humidity, net radiation flux, and subsurface temperature profiles over 2 days. The data validate a non-linear model of vapor mass fraction. Unlike heat, which conducts through grains, vapor percolates across the interstitial pore space by advection and diffusion. On time scales longer than evaporation, adsorbed films equilibrate with their surroundings and hinder molecular diffusion. Their non-linear coupling with subsurface temperature generates inflections in vapor profiles without counterpart in simpler diffusive systems. Pore advection arises as wind induces subtle pressure variations over the topography. During periods of aeolian transport, flowing sand dehydrates the surface intermittently, triggering evanescent vapor waves of amplitude decaying exponentially downward on a characteristic length implying an adsorption rate governed by a kinetic-limited activated process. Finally, the probe yields diffusive and advective exchanges with the atmospheric boundary layer. During the day, their combined flux is smaller than expected, yet nearly proportional to the difference between vapor mass fraction at the surface and aloft. Under stabler stratification at night, or during aeolian sand transport, this relation no longer holds.

**Plain Language Summary** Deserts inhale and exhale water vapor through their surface. Although this process affects the water balance over vast sand seas, it is poorly understood for want of sensitive instruments. We discover how it operates using a new probe that detects tiny amounts of moisture on sand grains. Our analysis reveals that vapor infiltration is considerably slower in dry sand, and that wind flowing over a dune creates weak internal air currents contributing to the transport of moisture. Their strength depends on dune location, wind speed and direction. When wind is strong enough to let dry sand meander over a dune, the resulting rapid variation in surface moisture sends evanescent waves of humidity downward. An analysis of these waves implies that water evaporation from individual sand grains behaves like a slow chemical reaction. The exchange of moisture with the atmosphere is not always driven by the difference between humidity at the dune surface and in the ambient, as current models assume, and it is weaker than they predict. In future, the new probe can be used as "ground truth" to calibrate satellite observations over deserts, explore extra-terrestrial environments holding scant water, and detect moisture contamination in pharmaceutical products.

## 1. Introduction

Sands of hyper-arid deserts subject to infrequent precipitation experience a diurnal exchange of water vapor between the ambient and a relatively shallow layer just below the surface (Kobayashi et al., 1986). Although free-liquid water is largely absent, this exchange can be intense in fog deserts where morning relative humidity is high (Cereceda et al., 2008; Lancaster et al., 1984). In turn, moisture profiles determine the viability of microorganisms sheltering from high temperatures (Cáceres et al., 2007). If sands imbibe morning dew fast enough (Gambaryan-Roisman, 2014; Johnson & Dettre, 1964), grains can bind with ephemeral capillary bridges forming at their contacts (Mitarai & Nori, 2006), thereby briefly affecting surface cohesion and the rate of
16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

## Abstract

Although vapor exchanged across hyper-arid surfaces without free liquid affects the water budget of sand seas, its mechanism is poorly documented for want of accurate instruments with fine spatial resolution. To rectify this, we report bulk density profiles and spatio-temporal variations of vapor mass fraction just below the surface of a mobile dune, acquired with a multi-sensor capacitance probe sensitive to tiny water films adsorbed on sand grains. We also record wind speed and direction, ambient temperature and relative humidity, net radiation flux, and subsurface temperature profiles over two days. The data validate a non-linear model of vapor mass fraction. Unlike heat, which conducts through grains, vapor percolates across the interstitial pore space by advection and diffusion. On time scales longer than evaporation, adsorbed films equilibrate with their surroundings and hinder molecular diffusion. Their non-linear coupling with subsurface temperature generates inflections in vapor profiles without counterpart in simpler diffusive systems. Pore advection arises as wind induces subtle pressure variations over the topography. During periods of aeolian transport, flowing sand dehydrates the surface intermittently, triggering evanescent vapor waves of amplitude decaying exponentially downward on a characteristic length implying an adsorption rate governed by a kinetic-limited activated process. Finally, the probe yields diffusive and advective exchanges with the atmospheric boundary layer. During the day, their combined flux is smaller than expected, yet nearly proportional to the difference between vapor mass fraction at the surface and aloft. Under stabler stratification at night, or during aeolian sand transport, this relation no longer holds.

38

## Plain Language Summary

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

Deserts inhale and exhale water vapor through their surface. Although this process affects the water balance over vast sand seas, it is poorly understood for want of sensitive instruments. We discover how it operates using a new probe that detects tiny amounts of moisture on sand grains. Our analysis reveals that vapor infiltration is considerably slower in dry sand, and that wind flowing over a dune creates weak internal air currents contributing to the transport of moisture. Their strength depends on dune location, wind speed and direction. When wind is strong enough to let dry sand meander over a dune, the resulting rapid variation in surface moisture sends evanescent waves of humidity downward. An analysis of these waves implies that water evaporation from individual sand grains behaves like a slow chemical reaction. The exchange of moisture with the atmosphere is not always driven by the difference between humidity at the dune surface and in the ambient, as current models assume, and it is weaker than they predict. In future, the new probe can be used as "ground truth" to calibrate satellite observations over deserts, explore extra-terrestrial environments holding scant water, and detect moisture contamination in pharmaceutical products.

54

# 1 Introduction

55

56

57

58

59

60

61

62

63

64

65

Sands of hyper-arid deserts subject to infrequent precipitation experience a diurnal exchange of water vapor between the ambient and a relatively shallow layer just below the surface (Kobayashi et al., 1986). Although free-liquid water is largely absent, this exchange can be intense in fog deserts where morning relative humidity is high (Lancaster et al., 1984; Cereceda et al., 2008). In turn, moisture profiles determine the viability of microorganisms sheltering from high temperatures (Cáceres et al., 2007). If sands imbibe morning dew fast enough (Gambaryan-Roisman, 2014; Johnson & Dettre, 1964), grains can bind with ephemeral capillary bridges forming at their contacts (Mitarai & Nori, 2006), thereby briefly affecting surface cohesion and the rate of aeolian transport (Ravi et al., 2006; Rotnicka, 2013). Understanding how and when hyper-arid sand surfaces exchange water vapor with the atmosphere is therefore crucial to subjects as diverse as zo-
**Resources:** M. Y. Louge, A. Valance, A. Ould el-Moctar, P. Chasle
**Software:** M. Y. Louge, J. Xu
**Supervision:** M. Y. Louge
**Validation:** M. Y. Louge, A. Valance
**Writing – original draft:** M. Y. Louge

aeolian transport (Ravi et al., 2006; Rotnicka, 2013). Understanding how and when hyper-arid sand surfaces exchange water vapor with the atmosphere is therefore crucial to subjects as diverse as zoology (Robinson & Barrows, 2013), microbiology (Heulin et al., 2012; Louge et al., 2013), botany (Wilcox et al., 2004), and desertification management (Fang et al., 2007; T. Wang et al., 2011).

Unfortunately, precise measurements of solid-bound water below the surface of hyper-arid sands are rare. They are usually achieved by weighing and drying samples collected over a relatively large volume (Dincer et al., 1974; Feigenwinter et al., 2020; J. Wang & Mitsuta, 1992; Kobayashi et al., 1991; Ritsema & Dekker, 1994). Using time-domain reflectometry (TDR), Shao et al. (2021) recently reported measurements of water volume fraction in the hyper-arid Taklamakan Desert. They also summarized the approach that is typically adopted to analyze the water retention of arid soils. In general, detailed time-histories of moisture depth profiles are not as well established as the corresponding temperature records (de Félice, 1968), particular near the surface where gradients of water mass fraction are largest.

To address this, Louge, Valance, Babah, et al. (2010) exploited a capacitance technique that Louge et al. (1998) had designed for alpine snowpacks to record solid volume fraction and solid-bound water mass fraction through the first 30 cm of dune sands, where moisture is most likely to be affected by diurnal changes in ambient conditions. However, because those authors did not repeat these measurements throughout the day, it remained unclear how moisture profiles evolved during a complete diurnal period. More significantly, their governing equations for mass transfer ignored water adsorption on the surface of sand grains, therefore overestimating the diffusion rate of moisture into dry sands, and they did not consider the possibility that advection could take place deeper than what they had calculated within sand ripples.

In this context, we designed a new probe to record low water mass fractions versus depth and time on a fine spatial resolution. Our simultaneous development of capacitance instruments for moisture and bulk density measurements in pharmaceutical powders conferred an unprecedented precision to this technique (Louge et al., 2021). The probe is intended for “hygroscopic” conditions, where water either exists as a vapor or is adsorbed on sand grains (Shahraeeni & Or, 2010). We inserted it in a mobile barchan dune of the Qatar desert. Its measurements at 2.7 min intervals revealed hitherto undocumented phenomena that similar instruments lacking fine spatial resolution could not discern.

First, the data set implies the existence of a slow wind-driven “seepage” advection flow in the pore space below the surface that augments water diffusion. Inspired by the theory of P. S. Jackson and Hunt (1975), we suggest that the advection derives from pressure variations on the dune topography. Beside recording stratigraphic layering of the subsurface, the probe also reveals evanescent waves of moisture coinciding with the onset of aeolian transport. The exponential decay of their amplitude with depth implies that evaporation acts as a slow activated reaction. Next, we explain why spatial variations of vapor mass fraction feature multiple moisture profiles inflexions without a counterpart in simpler advection-diffusion systems such as heat transfer through sands (Louge et al., 2013). Finally, the data set provides a unique “ground truth” for the moisture flux through a hyper-arid surface into the atmospheric boundary layer (ABL). As such, it has the potential to complement models of deeper moisture exchange in deserts, such as those of Kamai and Assouline (2018), Assouline and Kamai (2019), Shao et al. (2021) and, when the near-surface holds no liquid water, to revisit the classical formulations of Philip and De Vries (1957), Jury and Letey (1979), Parlange (1980), Brutsaert (1982, 1986), and Cahill and Parlange (1998) for water transport in wetter soils.

We begin with a description of the new probe and its response to the moisture content of sands. We then present the data set from all deployed instruments. We derive governing equations to interpret the resulting observations, and we use this framework to identify the principal mechanisms of hygroscopic subsurface moisture transport and dismiss those that are negligible. In future, the sensitivity of our technique, and the models that the new probe inspired, will make it possible to discern very low levels of moisture and understand its transport in the desert subsurface, or potentially in extra-terrestrial environments (Davidsson & Hosseini, 2021; Honniball et al., 2020).

## 2. Instruments

We designed two bespoke instruments for insertion through sand surfaces: a rectilinear lance with 15 temperature sensors (LM235 National Semiconductor temperature-sensitive Zener diodes), already exploited by Louge et al. (2013), and a capacitance probe with another 15 sensors recording moisture versus depth, shown in Figure 1
ology (Robinson & Barrows, 2013), microbiology (Heulin et al., 2012; Louge et al., 2013), botany (Wilcox et al., 2004), and desertification management (Fang et al., 2007; T. Wang et al., 2011).

Unfortunately, precise measurements of solid-bound water below the surface of hyper-arid sands are rare. They are usually achieved by weighing and drying samples collected over a relatively large volume (Dincer et al., 1974; Kobayashi et al., 1991; J. Wang & Mitsuta, 1992; Ritsema & Dekker, 1994; Feigenwinter et al., 2020). Using time-domain reflectometry (TDR), Shao et al. (2021) recently reported measurements of water volume fraction in the hyper-arid Taklamakan desert. They also summarized the approach that is typically adopted to analyze the water retention of arid soils. In general, detailed time-histories of moisture depth profiles are not as well established as the corresponding temperature records (de Félice, 1968), particular near the surface where gradients of water mass fraction are largest.

To address this, Louge, Valance, Babah, et al. (2010) exploited a capacitance technique that Louge et al. (1998) had designed for alpine snowpacks to record solid volume fraction and solid-bound water mass fraction through the first 30 cm of dune sands, where moisture is most likely to be affected by diurnal changes in ambient conditions. However, because those authors did not repeat these measurements throughout the day, it remained unclear how moisture profiles evolved during a complete diurnal period. More significantly, their governing equations for mass transfer ignored water adsorption on the surface of sand grains, therefore overestimating the diffusion rate of moisture into dry sands, and they did not consider the possibility that advection could take place deeper than what they had calculated within sand ripples.

In this context, we designed a new probe to record low water mass fractions vs depth and time on a fine spatial resolution. Our simultaneous development of capacitance instruments for moisture and bulk density measurements in pharmaceutical powders conferred an unprecedented precision to this technique (Louge et al., 2021). The probe is intended for “hygroscopic” conditions, where water either exists as a vapor or is adsorbed on sand grains (Shahraeeni & Or, 2010). We inserted it in a mobile barchan dune of the Qatar desert. Its measurements at 2.7 min intervals revealed hitherto undocumented phenomena that similar instruments lacking fine spatial resolution could not discern.

First, the data set implies the existence of a slow wind-driven “seepage” advection flow in the pore space below the surface that augments water diffusion. Inspired by the theory of P. S. Jackson and Hunt (1975), we suggest that the advection derives from pressure variations on the dune topography. Beside recording stratigraphic layering of the subsurface, the probe also reveals evanescent waves of moisture coinciding with the onset of aeolian transport. The exponential decay of their amplitude with depth implies that evaporation acts as a slow activated reaction. Next, we explain why spatial variations of vapor mass fraction feature multiple moisture profiles inflexions without a counterpart in simpler advection-diffusion systems such as heat transfer through sands (Louge et al., 2013). Finally, the data set provides a unique “ground truth” for the moisture flux through a hyper-arid surface into the atmospheric boundary layer. As such, it has the potential to complement models of deeper moisture exchange in deserts, such as those of Kamai and Assouline (2018); Assouline and Kamai (2019); Shao et al. (2021) and, when the near-surface holds no liquid water, to revisit the classical formulations of Philip and De Vries (1957); Jury and Letey Jr (1979); Parlange (1980); Brutsaert (1982, 1986); Cahill and Parlange (1998) for water transport in wetter soils.

We begin with a description of the new probe and its response to the moisture content of sands. We then present the data set from all deployed instruments. We derive governing equations to interpret the resulting observations, and we use this framework to identify the principal mechanisms of hygroscopic subsurface moisture transport and dismiss those that are negligible. In future, the sensitivity of our technique, and the mod-
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Fifteen-sensor capacitance probe. (a) Side view of the lance showing the rear guarded connector. (b) Front view with superimposed printed circuit board of 15 independent sensors numbered from the surface downward. Blue rectangles are maintained at sensor voltages. Surrounding red and green surfaces hold guard and ground voltages, respectively. (c) The probe before insertion into sands. In the background, a lance of similar geometry is used to assess sand looseness nearby before insertion. (d) GPS-recorded dune outline and brink (red and blue circles, respectively). The yellow triangle and gray square show insertion spot and weather station, respectively. For scale, the North arrow is 40 m-long.</figcaption></figure>

and described in Appendix A. We inserted these instruments perpendicularly to the sand free surface until their top skimmed it without impeding sand transport, through a relatively soft zone upwind of a barchan sand dune west of Mesaieed, Qatar at 25°00′34.7′′N, 51°20′24.9′′E. We mounted a Kipp and Zonen NR-Lite-2 radiometer nearby, and we buried the electronics, batteries and data acquisition system under a thin reflective emergency blanket approximately 1.5 m away. The system took approximately 20s to scan all temperature, moisture and radiation sensors every 2.7 min for nearly 2 days. In addition, we deployed a weather station above hard ground upstream of the dune at the location shown in Figure 1d to record ambient humidity and temperature at 0.9 and 2 m altitude, as well as wind speed at 0.9 and 2.4 m, and wind direction at 2.4 m, every minute. Michel et al. (2018) later exploited this aeolian data to calibrate optical satellite imagery for dune motion.

Louge et al. (2013) provided additional information on this site, including available long-term rain and wind records, location and speed of barchans in the region, measurements of relative humidity and temperature from a probe buried deeply in our dune over several months, modeling of deep water penetration and thermal energy, near-surface thermal exchanges, particle characterization (size distribution, grain composition), and biology. Arran (2018) also recorded detailed stratigraphy near the dune surface. Meanwhile, Abdul-Majid, Graw, et al. (2016) carried out the first characterization of the microbiota within mobile dunes in this field, including our own, which they dubbed “Nadine,” while Abdul-Majid, Chatziefthimiou, et al. (2016) considered their herpetology.

On long time scales, the capacitance probe is mainly sensitive to the mass fraction Ω of water adsorbed on sand grains, as discussed in Appendix A. When the vapor mass fraction Y (a.k.a. specific humidity) changes rapidly, the probe also discerns its variations versus depth, as Y quickly reaches an equilibrium with the thin solder mask protecting its electrodes. Each of the probe's 15 independent sensors records the real and imaginary parts of the complex effective sand dielectric constant ($K'_e - iK''_e$), where $i^2 = -1$, from which the “loss tangent”
els that the new probe inspired, will make it possible to discern very low levels of moisture and understand its transport in the desert subsurface, or potentially in extra-terrestrial environments (Honniball et al., 2020; Davidsson & Hosseini, 2021).

## 2 Instruments

We designed two bespoke instruments for insertion through sand surfaces: a rec-tilinear lance with 15 temperature sensors (LM235 National Semiconductor temperature-sensitive Zener diodes), already exploited by Louge et al. (2013), and a capacitance probe with another 15 sensors recording moisture vs depth, shown in Fig. 1 and described in Appendix A. We inserted these instruments perpendicularly to the sand free surface until their top skimmed it without impeding sand transport, through a relatively soft zone upwind of a barchan sand dune west of Mesaieed, Qatar at 25°00′34.7″N, 51°20′24.9″E. We mounted a Kipp & Zonen NR-Lite-2 radiometer nearby, and we buried the electronics, batteries and data acquisition system under a thin reflective emergency blanket approximately 1.5 m away. The system took approximately 20 s to scan all temperature, moisture and radiation sensors every 2.7 min for nearly two days. In addition, we deployed a weather station above hard ground upstream of the dune at the location shown in Fig. 1D to record ambient humidity and temperature at 0.9 m and 2 m altitude, as well as wind speed at 0.9 m and 2.4 m, and wind direction at 2.4 m, every minute. Michel et al. (2018) later exploited this aeolian data to calibrate optical satellite imagery for dune motion.

Louve et al. (2013) provided additional information on this site, including available long-term rain and wind records, location and speed of barchans in the region, measurements of relative humidity and temperature from a probe buried deeply in our dune over several months, modeling of deep water penetration and thermal energy, near-surface thermal exchanges, particle characterization (size distribution, grain composition), and biology. Arran (2018) also recorded detailed stratigraphy near the dune surface. Meanwhile, Abdul-Majid, Graw, et al. (2016) carried out the first characterization of the microbiota within mobile dunes in this field, including our own, which they dubbed “Nadine”, while Abdul-Majid, Chatziefthimiou, et al. (2016) considered their herpetology.

On long time scales, the capacitance probe is mainly sensitive to the mass fraction Ω of water adsorbed on sand grains, as discussed in Appendix A. When the vapor mass fraction Y (a.k.a. specific humidity) changes rapidly, the probe also discerns its variations vs depth, as Y quickly reaches an equilibrium with the thin solder mask protecting its electrodes. Each of the probe’s 15 independent sensors records the real and imaginary parts of the complex effective sand dielectric constant ($K'_e - iK''_e$), where $i^2 = -1$, from which the “loss tangent” $\tan \varphi \equiv K''_e/K'_e$ can be inferred (J. D. Jackson, 1999; Louge et al., 1997). A sensor of characteristic length $\ell_c$ exposed to air has a capacitance $\epsilon_0\ell_c$, where $\epsilon_0 \approx 8.854 \text{ fF/mm}$ is the permittivity in vacuo. Its introduction in a moist dielectric medium changes its impedance to an equivalent circuit consisting of a resistance $(2\pi f\epsilon_0K''_e\ell_c)^{-1}$ in parallel with a capacitance $\epsilon_0K'_e\ell_c$, where f is frequency of the applied voltage (Louve et al., 2021). Dry sand is a pure dielectric with $K''_e \approx 0$.

As Shahraeeni and Or (2010) showed, van der Waals attraction to the solid and binding liquid forces hold a thin layer of water molecules on the grain surface to a thickness

$$ \ell_w = \frac{\ell_0}{[T^* \ln(1/RH)]^{1/3}}, \qquad (1) $$

with characteristic length

$$ \ell_0 = \left[ \frac{\text{MW}_{\text{H}_2\text{O}}(-H)}{6\pi \rho_w \hat{R} T_{\text{st}}} \right]^{1/3}, \qquad (2) $$

where RH is relative humidity in the interstitial air surrounding a sand grain, MW$_{H_2O}$ is the molar mass of water, $H < 0$ is an effective Hamaker coefficient, $\hat{R} \approx 8.314 \text{ J/mole.K}$ is the fundamental gas constant, $\rho_w \approx 997 \text{ kg/m}^3$ is the density of liquid water, $T_{\text{st}} =$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Loss tangent tan φ ≡ K''/K' versus relative humidity RH_c during consecutive 4 hr plateaus in a moisture chamber calibration at 35.4°C using the capacitance probe of Louge, Valance, Ould el-Moctar, and Dupont (2010) at the base of a 1 cm column of Qatar dune sand. The line is Equation 4 with tan φ₀ ≈ 0.297. The diamond marks the largest stable state observed. Ascending and descending RH_c in the chamber are shown as open and filled circles, respectively, thus revealing negligible hysteresis. The filled triangle is another calibration test at 20.9°C, corrected for temperature with Equation 5, in which the probe of Figure 1 was covered with sand in the lab. This test lasted long enough for both sand and solder mask to reach equilibrium with the ambient RH_c ≈ 0.46.</figcaption></figure>

tan φ ≡ K''/K' can be inferred (J. D. Jackson, 1999; Louge et al., 1997). A sensor of characteristic length ℓ_c exposed to air has a capacitance ε₀ℓ_c, where ε₀ ≈ 8.854 fF/mm is the permittivity in vacuo. Its introduction in a moist dielectric medium changes its impedance to an equivalent circuit consisting of a resistance (2πfε₀K''ℓ_c)⁻¹ in parallel with a capacitance ε₀K'ℓ_c, where f is frequency of the applied voltage (Louge et al., 2021). Dry sand is a pure dielectric with K'' ≈ 0.

As Shahraeeni and Or (2010) showed, van der Waals attraction to the solid and binding liquid forces hold a thin layer of water molecules on the grain surface to a thickness

$$ \ell_w = \frac{\ell_0}{[T^* \ln(1/RH)]^{1/3}}, \quad (1) $$

with characteristic length

$$ \ell_0 = \left[ \frac{\text{MW}_{\text{H}_2\text{O}}(-H)}{6\pi \rho_w \hat{R}_{\text{st}}} \right]^{1/3}, \quad (2) $$

where RH is relative humidity in the interstitial air surrounding a sand grain, MW_H₂O is the molar mass of water, H < 0 is an effective Hamaker coefficient, ˆR ≈ 8.314 J/m K is the fundamental gas constant, ρ_w ≈ 997 kg/m³ is the density of liquid water, T_st = 298.15 K is the absolute temperature T at standard conditions, and T* ≡ T/T_st is its dimensionless counterpart (Iwamatsu & Horii, 1996). This “film adsorption” regime persists until water forms meniscii in small cavities through capillary condensation (Shahraeeni & Or, 2010). As Appendix B shows, the resulting mass fraction of water Ω_v in bulk sand of particle-size distribution f_M at equilibrium with the surrounding RH is given by the “isotherm” relation

$$ \Omega_v = \frac{\Omega_1}{[T^* \ln(1/RH)]^{1/3}}. \quad (3) $$

For sands of the dune under study, Louge et al. (2013) measured Ω₁ ≈ 0.0013 in a moisture chamber.

Because Ω_v, or its instantaneous counterpart Ω, are ratios of masses attached to individual grains, they do not depend on the solid volume fraction ν of the sand bed. Similarly, the loss tangent tan φ is weakly dependent on ν (Louge, Valance, Ould el-Moctar, & Dupont, 2010), but it rises with Ω. To establish their relation, we covered the probe of Louge, Valance, Ould el-Moctar, and Dupont (2010) with a 1 cm-thick layer of Qatar dune sand in a humidity chamber held at T_c ≈ 35.4°C with relative humidity in the range 0.32 < RH_c < 0.9 (Figure 2). We obtained stable signals until RH_c = RH_max ~ 0.821, at which point K'' became too large. We attribute this maximum RH to the onset of capillary condensation. In the chamber, the loss tangent conformed to

$$ \tan \varphi \approx \tan \varphi_0 \times \text{RH}_c / (1 - \text{RH}_c), \quad (4) $$

where RH_c is the relative humidity around grains during calibration.

As sand temperature T varies substantially below the dune surface (Louge et al., 2013), Equation 4 must be extrapolated to values of T other than the calibration temperature T_c. To that end, we use Equation 3 to convert RH_c to RH_v = exp[(T_c/T)ln RH_c]. Substituting Equation 4, we find the equilibrium RH around grains at T,

$$ \text{RH}_v = \exp \left[ - \left( \frac{T_c}{T} \right) \ln \left( 1 + \frac{\tan \varphi_0}{\tan \varphi} \right) \right]. \quad (5) $$

Approximating the saturation pressure of water using “Antoine's law”
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Fifteen-sensor capacitance probe. (A) Side view of the lance showing the rear guarded connector. (B) Front view with superimposed printed circuit board of 15 independent sensors numbered from the surface downward. Blue rectangles are maintained at sensor voltages. Surrounding red and green surfaces hold guard and ground voltages, respectively. (C) The probe before insertion into sands. In the background, a lance of similar geometry is used to assess sand looseness nearby before insertion. (D) GPS-recorded dune outline and brink (red and blue circles, respectively). The yellow triangle and grey square show insertion spot and weather station, respectively. For scale, the North arrow is 40 m-long.</figcaption></figure>
$$p_{\text{sat}} \approx p_{A_0} \exp[-T_{A_0} / (T - T_{A_0})], \qquad (6)$$

with $p_{A_0} \approx 1.24 \ 10^{10}$ Pa, $T_{A_0} \approx 3841.2^{\circ}$K and $T_{0_\alpha} \approx 45.2^{\circ}$K, we convert RH to water vapor mass fraction using

$$Y = \frac{\text{RH}}{M^*(p/p_{\text{sat}}) - \text{RH}(M^* - 1)} \approx \frac{p_{\text{sat}}}{M^*p} \text{RH}, \qquad (7)$$

where $M^* \equiv \text{MW}_{\text{air}}/\text{MW}_{\text{H}_2\text{O}} \approx 1.61$ is the ratio of the molar masses of air and water.

Louge et al. (2013) recorded sand composition for the dune under study using a QEMscan (Gottlieb et al., 2000), reported it in Appendix E of their Supporting Information, and verified that Equation 3 has merit. Other mineral compositions can exhibit qualitatively different isotherms (Shang et al., 1995). However, because the latter have a weaker dependence on temperature than saturation pressure does in Equation 6, our modeling framework has broad applicability.

### 3. Temperature Profiles

In Appendix C, we calculate the derivatives of $\Omega_\alpha$ from the isotherm of Equation 3. This reveals that the amount of water held on sand grains is a strong function of temperature through the exponential in Antoine's law (Equation 6). However, because negligible latent heat is associated with the thin layer in Equation 1, the heat equation is not affected its adsorption or desorption. In other words, Equation 6 “one-way-couples” spatiotemporal variations of $Y$ to $T$, thereby making it crucial to record subsurface temperature, and to evaluate its partial derivatives with time and space accurately. Such evaluation first requires a fit of temperature data that is strictly consistent with the solution of the linear heat equation in time $t$ and depth $x$

$$\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}, \qquad (8)$$

which Louge et al. (2013) satisfactorily integrated subject to the surface boundary condition (BC)

$$-k_s \frac{\partial T}{\partial x}|_{x=0} = \dot{q}''_{\text{rad}} + \dot{q}''_{\text{wind}}, \qquad (9)$$

where $\alpha$ and $k_s$ are, respectively, the effective thermal diffusivity and conductivity of bulk sand. In Equation 9, $\dot{q}''_{\text{rad}}$ is the measured net radiation flux striking the dune and $\dot{q}''_{\text{wind}}$ is the wind-driven convective flux calculated from the Monin-Obukhov model of the coupled momentum and heat ABLs summarized in Appendix D.

Rather than adding Equation 8 to the non-linear set of governing equations to be introduced in Section 6, it is more accurate to access time and space derivatives of $T$ by fitting its record to a sum of four solution harmonics of Equation 8 (Carslaw & Jaeger, 1959),

$$T \approx T_{\infty} + \left[ \sum_{i=1}^{4} T_i \times \mathbb{H}(x, t; i) \right] + T_{1/n} \times \mathbb{H}(x, t; 1/n), \qquad (10)$$

where

$$\mathbb{H}(x, t; i) \equiv \exp\left(-x\sqrt{\frac{i\pi}{J\alpha}}\right) \sin\left[2i\pi\left(\frac{t-t_i}{J}\right) - x\sqrt{\frac{i\pi}{J\alpha}}\right], \qquad (11)$$

and $J = 24$ hr is the diurnal period. The term $\mathbb{H}(x, t; 1/n)$ is a subharmonic capturing the increase in mean temperature on the second day due to slower wind speed. A least squares fit to the data yielded $T_\infty \approx 297.0^{\circ}$K at large depth, amplitudes ($T_{1/n}$, $T_1$, $T_2$, $T_3$, $T_4$) $\approx (19.1, 13.5, -5.5, 0.72, 0.65)^\circ$C, phase lead times ($t_{1/n}$, $t_1$, $t_2$, $t_3$, $t_4$) $\approx (0.789, 0.313, 0.145, 0.441, 0.052) \times J$, $n \approx 822$, and $\alpha \approx 3.53 \ 10^{-7} \text{ m}^2/\text{s}$, which we adopted as the observed sand heat diffusivity. Equations 8–10 suggest that depth and time should be made dimensionless as $x^* \equiv x/\sqrt{\alpha J}$ and $t^* \equiv t/J$, respectively.
298.15 K is the absolute temperature T at standard conditions, and $T^* \equiv T/T_{st}$ is its dimensionless counterpart (Iwamatsu & Horii, 1996). This “film adsorption” regime persists until water forms meniscii in small cavities through capillary condensation (Shahraeeni & Or, 2010). As Appendix B shows, the resulting mass fraction of water $\Omega_e$ in bulk sand of particle-size distribution $f_M$ at equilibrium with the surrounding RH is given by the “isotherm” relation

$$ \Omega_e = \frac{\Omega_1}{[T^* \ln(1/RH)]^{1/3}} \qquad (3) $$

For sands of the dune under study, Louge et al. (2013) measured $\Omega_1 \simeq 0.0013$ in a moisture chamber.

Because $\Omega_e$, or its instantaneous counterpart $\Omega$, are ratios of masses attached to individual grains, they do not depend on the solid volume fraction $\nu$ of the sand bed. Similarly, the loss tangent $\tan \varphi$ is weakly dependent on $\nu$ (Louve, Valance, el Moctar, & Dupont, 2010), but it rises with $\Omega$. To establish their relation, we covered the probe of Louge, Valance, el Moctar, and Dupont (2010) with a 1 cm-thick layer of Qatar dune sand in a humidity chamber held at $T_c \simeq 35.4^\circ C$ with relative humidity in the range $0.32 < RH_c < 0.9$ (Fig. 2). We obtained stable signals until $RH_c = RH_{max} \sim 0.821$, at which point $K''_e$ became too large. We attribute this maximum RH to the onset of capillary condensation. In the chamber, the loss tangent conformed to

$$ \tan \varphi \simeq \tan \varphi_0 \times RH_c / (1 - RH_c), \qquad (4) $$

where $RH_c$ is the relative humidity around grains during calibration.

As sand temperature $T$ varies substantially below the dune surface (Louve et al., 2013), Eq. (4) must be extrapolated to values of $T$ other than the calibration temperature $T_c$. To that end, we use Eq. (3) to convert $RH_c$ to $RH_e = \exp[(T_c/T) \ln RH_c]$. Substituting Eq. (4), we find the equilibrium RH around grains at $T$,

$$ RH_e = \exp \left[ - \left( \frac{T_c}{T} \right) \ln \left( 1 + \frac{\tan \varphi_0}{\tan \varphi} \right) \right]. \qquad (5) $$

Approximating the saturation pressure of water using “Antoine’s law”

$$ p_{\text{sat}} \simeq p_{A_0} \exp[-T_{A_0} / (T - T_{A_0})], \qquad (6) $$

with $p_{A_0} \simeq 1.24 \times 10^{10}$ Pa, $T_{A_0} \simeq 3841.2^\circ K$ and $T_{0a} \simeq 45.2^\circ K$, we convert RH to water vapor mass fraction using

$$ Y = \frac{\text{RH}}{M^*(p/p_{\text{sat}}) - \text{RH}(M^* - 1)} \simeq \frac{p_{\text{sat}}}{M^*p} \text{RH}, \qquad (7) $$

where $M^* \equiv \text{MW}_{\text{air}}/\text{MW}_{\text{H}_2\text{O}} \simeq 1.61$ is the ratio of the molar masses of air and water.

Louve et al. (2013) recorded sand composition for the dune under study using a QEMscan (Gottlieb et al., 2000), reported it in Appendix E of their supporting information, and verified that Eq. (3) has merit. Other mineral compositions can exhibit qualitatively different isotherms (Shang et al., 1995). However, because the latter have a weaker dependence on temperature than saturation pressure does in Eq. (6), our modeling framework has broad applicability.

## 3 Temperature profiles

In Appendix C, we calculate the derivatives of $\Omega_e$ from the isotherm of Eq. (3). This reveals that the amount of water held on sand grains is a strong function of temperature through the exponential in Antoine’s law (Eq. 6). However, because negligible latent heat is associated with the thin layer in Eq. (1), the heat equation is not affected
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Measurements of 19–21 March 2011. Animations, probe depths, and data are provided as Supporting Information. (a) Vapor mass fraction *Y* recorded by capacitance sensors of indices in Figure 1, with matching line colors; black dashed and solid lines are ambient *Y*<sub>b</sub> from the weather station at altitudes *z*<sub>B</sub> = 0.9 m and *z*<sub>T</sub> = 2 m, respectively. Values of *Y* appear as gray symbols in the Movie S1. (b) Subsurface temperatures *T* recorded by the shorter probe of Louge et al. (2013) (increasing indices downward), and appearing as yellow symbols in Movie S4. Black dashed and solid lines are ambient temperatures *T*<sub>B</sub> and *T*<sub>T</sub> at *z*<sub>B</sub> and *z*<sub>T</sub>, respectively. (c) Net thermal fluxes by radiation (red) and wind-driven convection (blue) at the dune surface, positive along depth *x* into the dune. *q*''<sub>wind</sub> is calculated as explained in the Supporting Information of Louge et al. (2013) with relative shear velocity *v*/*U* ≃ 0.036 and roughness *z*<sub>0</sub> ≃ 16 µm recorded at the probe location in January 2017. (d) Wind speeds at the weather station at *z* = 0.9 m (blue) and 2.4 m (red) with 0.3 m/s detection threshold. The common abscissa is time *t*\* dimensionless with *J* = 24 hr, with origin at 00:00 Qatar time (GMT + 3) on 19 March 2011. The vertical solid and dashed orange lines mark sunrise and sunset, respectively. Missing data are due to excessive *K*''\* that induced spurious signals from the capacitance electronics at LF radio frequencies that contaminated both records of *Y* and *T*.</figcaption></figure>

**4. Measurements**

Figure 3 shows the entire data set, which we submit as Supporting Information. It includes vapor mass fractions $Y$ calculated from Equations 5–7, measured subsurface temperatures $T$ used in Equation 10, net radiation flux $q_{\text{rad}}''$ recorded above the probes, the wind-driven heat flux through the surface $q_{\text{wind}}''$ that Louge et al. (2013) modeled, and wind speed $U$ at two altitudes upwind of the dune. We will use this record to validate the model framework in Section 6.

The first unexpected feature of the mass fraction record consisted of multiple inflections versus depth, which have no counterpart in the diffusive heat transfer captured by Equation 10. As Figure 4 and its animation in the Supporting Information illustrate, *Y* exhibited an S-shaped profile with a partial minimum that is co-located with the partial minimum in *T*. This suggests that the peculiar inflections of *Y* were associated with the steep exponential dependence of saturation pressure on *T* modeled by Equation 6 and quantified in Appendix C. The integration of mass transfer partial differential equation (PDE) in Section 7 will confirm this explanation.

A more puzzling observation was sudden jumps in *Y*, detected from $t^* \approx 1.4$ to 1.63 on 20 March 2011. In that interval, *Y* switched from a dry to a humid state in rapid succession. The onset of this apparent signal bistability coincided with an upstream wind speed rising above *U* ~ 7 m/s, and it disappeared as *U* persistently dropped below this value. During this period, wind occasionally uncovered the top of our probes by a few millimeters. When this happened, we replenished the missing sand gently by hand to keep both probes flush with the surface.

Later, no such bistability arose during the corresponding period on 21 March (2.4 $\lesssim$ $t^* \lesssim$ 2.6), when *U* remained below 6 m/s, and $Y_1$ exhibited single values nearly identical to the more humid state of 20 March, implying less efficient surface drying ($Y_1$ > $Y_a$). Acoustic detection of sand impact and anemometry carried out on this dune in January 2017 yielded thresholds for particle motion ~7 m/s. These observations suggest that signal bistability was likely associated with intermittent aeolian sand transport.
its adsorption or desorption. In other words, Eq. (6) “one-way-couples” spatio-temporal variations of Y to T, thereby making it crucial to record subsurface temperature, and to evaluate its partial derivatives with time and space accurately. Such evaluation first requires a fit of temperature data that is strictly consistent with the solution of the linear heat equation in time t and depth x

$$ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}, \qquad (8) $$

which Louge et al. (2013) satisfactorily integrated subject to the surface boundary condition

$$ -k_s \frac{\partial T}{\partial x} |_{x=0} = \dot{q}''_{\text{rad}} + \dot{q}''_{\text{wind}}, \qquad (9) $$

where $\alpha$ and $k_s$ are, respectively, the effective thermal diffusivity and conductivity of bulk sand. In Eq. (9), $\dot{q}''_{\text{rad}}$ is the measured net radiation flux striking the dune and $\dot{q}''_{\text{wind}}$ is the wind-driven convective flux calculated from the Monin-Obukhov model of the coupled momentum and heat atmospheric boundary layers summarized in Appendix D.

Rather than adding Eq. (8) to the non-linear set of governing equations to be introduced in section 6, it is more accurate to access time and space derivatives of T by fitting its record to a sum of four solution harmonics of Eq. (8) (Carslaw & Jaeger, 1959),

$$ T \approx T_{\infty} + \left[ \sum_{i=1}^{4} T_i \times \mathbb{H}(x, t; i) \right] + T_{1/n} \times \mathbb{H}(x, t; 1/n), \qquad (10) $$

where

$$ \mathbb{H}(x, t; i) \equiv \exp \left( -x \sqrt{\frac{i\pi}{J\alpha}} \right) \sin \left[ 2i\pi \left( \frac{t-t_i}{J} \right) - x \sqrt{\frac{i\pi}{J\alpha}} \right], \qquad (11) $$

and $J = 24$ hr is the diurnal period. The term $\mathbb{H}(x, t; 1/n)$ is a subharmonic capturing the increase in mean temperature on the second day due to slower wind speed. A least-squares fit to the data yielded $T_{\infty} \approx 297.0^{\circ}\text{K}$ at large depth, amplitudes $(T_{1/n}, T_1, T_2, T_3, T_4) \approx (19.1, 13.5, -5.5, 0.72, 0.65)^{\circ}\text{C}$, phase lead times $(t_{1/n}, t_1, t_2, t_3, t_4) \approx (0.789, 0.313, 0.145, 0.441, 0.052) \times J$, $n \approx 822$, and $\alpha \approx 3.53 \times 10^{-7} \text{m}^2/\text{s}$, which we adopted as the observed sand heat diffusivity. Equations (8)–(10) suggest that depth and time should be made dimensionless as $x^* \equiv x/\sqrt{\alpha J}$ and $t^* \equiv t/J$, respectively.

# 4 Measurements

Figure 3 shows the entire data set, which we submit as supporting information. It includes vapor mass fractions Y calculated from Eqs. (5)–(7), measured subsurface temperatures T used in Eq. (10), net radiation flux $\dot{q}''_{\text{rad}}$ recorded above the probes, the wind-driven heat flux through the surface $\dot{q}''_{\text{wind}}$ that Louge et al. (2013) modeled, and wind speed U at two altitudes upwind of the dune. We will use this record to validate the model framework in section 6.

The first unexpected feature of the mass fraction record consisted of multiple inflections vs depth, which have no counterpart in the diffusive heat transfer captured by Eq. (10). As Fig. 4 and its animation in the supporting information illustrate, Y exhibited an S-shaped profile with a partial minimum that is co-located with the partial minimum in T. This suggests that the peculiar inflections of Y were associated with the steep exponential dependence of saturation pressure on T modeled by Eq. (6) and quantified in Appendix C. The integration of mass transfer partial differential equations in section 7 will confirm this explanation.

A more puzzling observation was sudden jumps in Y, detected from $t^* \approx 1.4$ to 1.63 on March 20, 2011. In that interval, Y switched from a dry to a humid state in rapid
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Snapshots of Movie S1, showing subsurface depth profiles of temperature in °C (Equation 10, red line) and vapor mass fraction $Y_v$ at equilibrium with the measured Ω (gray symbols) at 16:09 Qatar time ($t^* \approx 1.67$) on 20 March (left) and 09:00 ($t^* \approx 2.38$) on 21 March (right). The blue line is the model implemented in Section 10. The yellow arrow points to co-located partial minima in the Y- and T-profiles. Bars and arrows represent scalars and vectors with lengths proportional to ambient temperature, wind speed, ambient water vapor mass fraction, dimensionless surface fluxes of vapor diffusion $\Psi_s^*/(1-\nu)$ and vapor advection $u_s^* Y_s p^* / T_s^*$ defined in Section 8. On the thin dotted line drawn for $x > 220$ mm, the vapor mass fraction Y is at equilibrium with a constant $\Omega = 1.65 \Omega_1$. The horizontal dotted line marks the sand surface. The horizontal gray and orange timeline indicates progression along the 19–21 March record.</figcaption></figure>

During such transport, mass transfer was enhanced by the presence of saltating particles, which quickly adopted a value of Ω in equilibrium with the ambient $Y_a$. Upon their intermittent deposition, these particles also imposed their lower moisture content on surface grains. Consequently, the drier envelope of water vapor mass fraction $Y_1$ at the sensor closest to the surface (blue line) was nearly identical to its ambient counterpart $Y_a$ recorded independently at the weather station (black line). (The coincidence of $Y_1$ and $Y_a$, which were measured by instruments of widely different principles, also speaks to the accuracy of the capacitance probes.) Without deposition of dry saltating sands, $Y_1$ reverted to the mass fraction $>Y_a$ that grains normally held at the surface. Note that we rarely observed an intermediate mass fraction between high and low states. Instead, the abrupt transition between them likely arose from rapid variations in saltation solid volume fraction, which typically accompany the passage of sharp edges of particles clusters meandering above the dune.

Remarkably, these perturbations in $Y_1$ recorded at the surface propagated quickly and deeply. This suggests that mass transfer was governed by two widely different time scales, namely the diurnal period $J$ that regulates the smooth evolution of $Y$ and $T$, and a much shorter one that allowed the rapid propagation of waves carrying surface disturbances to the subsurface. Section 12 will develop a stability analysis explaining this behavior.

Above the dune, wind speed rose during the day as soon as a net radiation flux stroke the desert floor (Wyngaard, 2010). An unstable ABL then developed, in which ambient temperature grew with elevation, but remained lower than its counterpart on the surface. As the radiation flux turned upward before sunset, so did the gradient of ambient temperature. During the night, surface temperature was lower than the ambient, and both progressively decreased in tandem. Low-level jets intermittently perturbed the nocturnal stable boundary layer (Banta et al., 2007; Klein et al., 2016).

The absence of significant noise in our capacitance record allowed us to differentiate spatiotemporal variations of $Y$ and use Fick's law to calculate the vapor mass flux $\dot{\Psi}_s'' = -\rho(1-\nu)(D/\bar{\omega})(\partial Y/\partial x)_s$ at the sand surface, where $D$ is the diffusion coefficient of water vapor at temperature $T$ and $\bar{\omega}$ is the tortuosity of porous sand (Shen & Chen, 2007). As Section 8 will show, this flux was augmented by a slow “seepage” advection within the pore space driven by wind blowing over topographical variations of the dune. Together, diffusion and advection fluxes provided a BC for the ABL above our instruments.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Loss tangent $\tan \varphi \equiv K_e''/K_e'$ vs relative humidity RH$_c$ during consecutive 4 hr plateaus in a moisture chamber calibration at 35.4 $^{\circ}$C using the capacitance probe of Louge, Valance, el Moctar, and Dupont (2010) at the base of a 1 cm column of Qatar dune sand. The line is Eq. (4) with $\tan \varphi_0 \simeq 0.297$. The diamond marks the largest stable state observed. Ascending and descending RH$_c$ in the chamber are shown as open and filled circles, respectively, thus revealing negligible hysteresis. The filled triangle is another calibration test at 20.9 $^{\circ}$C, corrected for temperature with Eq. (5), in which the probe of Fig. 1 was covered with sand in the lab. This test lasted long enough for both sand and solder mask to reach equilibrium with the ambient RH$_c \simeq 0.46$.</figcaption></figure>

succession. The onset of this apparent signal bistability coincided with an upstream wind speed rising above $U \sim 7$ m/s, and it disappeared as $U$ persistently dropped below this value. During this period, wind occasionally uncovered the top of our probes by a few millimeters. When this happened, we replenished the missing sand gently by hand to keep both probes flush with the surface.

Later, no such bistability arose during the corresponding period on March 21 (2.4 $\lesssim$ $t^* \lesssim 2.6$), when $U$ remained below 6 m/s, and $Y_1$ exhibited single values nearly identical to the more humid state of March 20, implying less efficient surface drying ($Y_1 > Y_a$). Acoustic detection of sand impact and anemometry carried out on this dune in January 2017 yielded thresholds for particle motion $\sim 7$ m/s. These observations suggest that signal bistability was likely associated with intermittent aeolian sand transport.

During such transport, mass transfer was enhanced by the presence of saltating particles, which quickly adopted a value of $\Omega$ in equilibrium with the ambient $Y_a$. Upon their intermittent deposition, these particles also imposed their lower moisture content on surface grains. Consequently, the drier envelope of water vapor mass fraction $Y_1$ at the sensor closest to the surface (blue line) was nearly identical to its ambient counterpart $Y_a$ recorded independently at the weather station (black line). (The coincidence of $Y_1$ and $Y_a$, which were measured by instruments of widely different principles, also speaks to the accuracy of the capacitance probes). Without deposition of dry saltating sands, $Y_1$ reverted to the mass fraction $> Y_a$ that grains normally held at the surface. Note that we rarely observed an intermediate mass fraction between high and low states. Instead, the abrupt transition between them likely arose from rapid variations in saltation solid volume fraction, which typically accompany the passage of sharp edges of particles clusters meandering above the dune.

Remarkably, these perturbations in $Y_1$ recorded at the surface propagated quickly and deeply. This suggests that mass transfer was governed by two widely different time
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Real part of the effective dielectric constant $K'_e$ versus loss tangent $\tan \varphi$ for all valid data points (small circles), shown for sensor 4, 9, and 13. The line is the Böttcher model in Equation 12 obtained with coefficients ($k_0, k_1, k_2, k_3$) $\approx$ (6.7, 20, -7.2, 3.9) from Equation 13 and volume fractions at all 15 sensors that minimize the cost function in Equation 14. Dashed lines above and below the main fit, as well as error bars in inset (b), mark $\pm$ one sample standard deviation in $\nu$ at the corresponding sensor. Inset (a) an excavation conducted in April 2014 reveals a layered stratigraphy after pouring water laced with blue food coloring through the surface. The stratigraphy occasionally appears on the surface after rain partially dries up (inset (c), picture taken on the same dune in January 2014). Inset (b) profile of $\nu$ versus depth $x$ (cm) perpendicular to the dune surface, calculated from Equations 12–14. The dashed line is a visual spline through the data.</figcaption></figure>

## 5. Stratigraphy

By recording both $K'_e$ and $K''_e$, the capacitance probe yielded not only $\Omega$ through $K''_e/K'_e$, but also the time-invariant local “bulk density” $\rho_b \approx \rho_p \nu$ at each sensor, where $\nu$ is the solid volume fraction and $\rho_p$ the material density of sand grains. As Louge, Valance, Ould el-Moctar, and Dupont (2010) showed, $K'_e$ is approximated by the model of Böttcher and Bordewijk (1978),

$$K'_e = \frac{1}{4} \left\{ 2 - K_p + 3\nu(K_p - 1) + \sqrt{8K_p + [2 - K_p + 3\nu(K_p - 1)]^2} \right\}, (12)$$

where $K_p$ is an equivalent material dielectric constant that rises with $\Omega$, and therefore with $\tan \varphi$. Rather than evaluating $K_p$ by staging several volume fractions in the laboratory, we found it more convenient to derive it directly from field data. In this method, we model $K_p$ as

$$K_p \approx k_0 + k_1 \tan \varphi + k_2 \tan^2 \varphi + k_3 \tan^3 \varphi. \quad (13)$$

From Equation 12, $K'_e$ and $K''_e = \tan \varphi \times K'_e$ become, respectively, the functions $f$ and $f'$ of the variables ($\nu_i, \tan \varphi_j; k_0, k_1, k_2, k_3$) recorded at sensor of index $i$ and time $j$. The values of $k_0$ through $k_3$ and $\nu_1$ through $\nu_{15}$ are those yielding the smallest cost function $C$ from the entire data set of $N = 11,557$ valid points acquired with all sensors,

$$C(\nu_1, \nu_2, \dots, \nu_{15}; k_0, k_1, k_2, k_3) = \sum_{j=1}^{N} \sum_{i=1}^{15} [K'_{e}(i, j) - f']^2 + [K''_{e}(i, j) - f'']^2. \quad (14)$$

Adding terms of order higher than $\tan^3 \varphi$ in Equation 13 ($k_4, k_5$, etc.) did not reduce $C$ and thus improve results.

As Figure 5 shows, the first 14 sensors spanned a wide enough range of $\tan \varphi$ to discern variations of $\nu$ with depth. Typically, such variations are associated with layers sorted and deposited by successive avalanches on the downwind face of the mobile dune (Bristow et al., 2007; Fan et al., 2012; Harari, 1996; Kocurek, 1991; Kleinhans, 2004; Qian et al., 2014; Reesink & Bridge, 2009), then sequestered for several months until aeolian transport uncovers them upwind (Louge et al., 2013). On other visits to the same sand dune, the imprint of earlier positions of the lee face was visible for a short time after rain (Figure 5c), suggesting preferential liquid penetration (Al-Shukaili et al., 2019). A complex layered stratigraphy was also evident upon excavation (Figure 5a).

At the specific location of our measurements, Figure 5 reveals relatively light sands with higher bulk density on the surface than below, consistent with the soft region that we selected for easy insertion of the probes. Its inset b also hints at an uneven stratigraphy on a wavelength $\Delta x \sim 3-5$ cm. Because the upwind and downwind faces of the dunes made angles of $\theta_u \approx 7^\circ$ and $\theta_d \approx 31^\circ$ from the horizontal, our capacitance probe inserted approximately normal to the surface along $x$ detected strata separated by a distance $\Delta x \cos(\theta_u + \theta_d) \approx 2-4$ cm. Finally, although the profile in Figure 5b resembled the layering shown in insets a and c, the Böttcher model of Equation 12 likely underestimated $\nu$.

## 6. Modeling Framework

Water evaporation from sands is complicated by capillary forces retaining liquid within its porous structure (Lehmann & Or, 2009; Or et al., 2013; Prat, 2002). Its analysis in the vadose zone therefore requires a model for the hysteretic process of wetting and draining in the unsaturated porous medium (Xu & Louge, 2015), typically obtained by invoking an empirical expression for the water retention curve (Assouline & Kamai, 2019; Kamai & Assouline, 2018; Shao et al., 2021; Vanderborght et al., 2017).
<figure><img src="image_4.png" /></figure>

**Figure 3.** Measurements of March 19-21, 2011. Animations, probe depths and data are pro- vided as supporting information. (A) Vapor mass fraction Y recorded by capacitance sensors of indices in Fig. 1, with matching line colors; black dashed and solid lines are ambient Y_a from the weather station at altitudes z_B = 0.9 m and z_T = 2 m, respectively. Values of Y appear as grey symbols in the movie *MovieDuneSubsurfaceTemperature&VaporMassFraction.mp4*. (B) Subsurface temperatures T recorded by the shorter probe of Louge et al. (2013) (increasing indices downward), and appearing as yellow symbols in *MovieDuneSubsurfaceTemperatur- eRadiation&Wind.mp4*. Black dashed and solid lines are ambient temperatures T_B and T_T at z_B and z_T, respectively. (C) Net thermal fluxes by radiation (red) and wind-driven convection (blue) at the dune surface, positive along depth x into the dune. *q''*<sub>wind</sub> is calculated as explained in the supporting information of Louge et al. (2013) with relative shear velocity *v*<sub>τ</sub>/U ≃ 0.036 and roughness z₀ ≃ 16 µm recorded at the probe location in January 2017. (D) Wind speeds at the weather station at z = 0.9 m (blue) and 2.4 m (red) with 0.3 m/s detection threshold. The common abscissa is time *t*<sup>*</sup> dimensionless with J = 24 hr, with origin at 00:00 Qatar time (GMT+3) on March 19, 2011. The vertical solid and dashed orange lines mark sunrise and sun- set, respectively. Missing data are due to excessive *K''*<sub>e</sub> that induced spurious signals from the capacitance electronics at LF radio frequencies that contaminated both records of Y and T.
The process is simpler near the surface of hyper-arid deserts, where water exists either as a vapor, or as a film directly adsorbed on solid grains (Shahraeeni & Or, 2010). Here, the drying front that Shokri et al. (2008) and Lehmann et al. (2008) documented lies deeper within the dune, and mass transfer near the surface conforms to the last stage of the evaporation process (Shahraeeni et al., 2012).

Such was the case for the data set of Section 4, which suggested a strong spatiotemporal coupling between subsurface variations of vapor mass fraction and temperature, a situation that is difficult to reproduce in laboratory experiments on the drying process (Davarzani et al., 2014; Gao et al., 2018; Z. Li et al., 2020).

In this section, we derive a set of governing equations to interpret these observations in the limit where interstitial water vapor is at equilibrium with water adsorbed on solid grains. We will use their framework to identify the principal mechanisms of moisture transport in the dune subsurface, and dismiss those that are negligible.

In the gas phase, the balance of total mass in an elementary slice of sand at depth x is

$$ (1 - \nu) \frac{\partial \rho}{\partial t} + (1 - \nu) \frac{\partial (\rho u)}{\partial x} = - \rho_p \nu \frac{\partial \Omega}{\partial t}, \quad (15) $$

where the right-hand term represents the rate of evaporation from grains of material density $\rho_p$ on which water was adsorbed, later considered in Section 11 and Appendix G. The variables $\rho$, $u$, $\Omega$, and $\nu$ in this classical mass balance (Bird et al., 2007; Kays & Crawford, 1980) are implicitly coarse-grained on the scale of a few particles, as suggested by T. B. Anderson and Jackson (1967). For clarity of derivations, we ignore relatively small variations of $\nu$ mentioned in Section 5, so $(1 - \nu)$ can be brought out of the spatial derivative. However, prior knowledge of a stratigraphy like that shown in Figure 5 can be substituted in the governing equations and its numerical solutions.

In Equation 15, the density of humid air $\rho = \rho_{H_2O} + \rho_{air}$ is the sum of the partial densities of water vapor $\rho_{H_2O} \equiv \rho Y = MW_{H_2O} \chi_{H_2O} (p/\hat{R}T)$ and dry air $\rho_{air} = MW_{air} \chi_{air} (p/\hat{R}T)$, with complementary mole fractions of water and air $\chi_{H_2O} + \chi_{air} = 1$, so vapor mass fraction is $Y = \chi_{H_2O}MW_{H_2O}/[\chi_{H_2O}MW_{H_2O} + \chi_{air}MW_{air}]$. Meanwhile, dry air density $\rho_{air} = \rho_{st} p^*/T^*$ can be expressed in terms of its counterpart $\rho_{st} \approx 1.185 \text{ kg/m}^3$ under standard conditions at $T_{st} = 298.15^\circ\text{K}$ and $p_{st} = 1.013 \times 10^5 \text{ Pa}$, where $p^* \equiv p/p_{st}$ and $T^* \equiv T/T_{st}$ are dimensionless absolute pressure and temperature, respectively. Using these relations, moist air density is

$$ \rho = \rho_{st} \frac{p^*/T^*}{1 + Y(M^* - 1)}, \quad (16) $$

where $M^* \equiv MW_{air}/MW_{H_2O}$. Because $Y(M^* - 1) \ll 1$, $\rho$ is almost equal to the density of air at the same $p$ and $T$. Therefore, although numerical solutions will adopt Equation 16, we will lighten derivations in this article by simply writing $\rho \approx \rho_{st} p^*/T^*$.

The “seepage” interstitial gas velocity $u$ is given by Darcy's law

$$ u = - \frac{K}{\mu (1 - \nu)} \frac{\partial p}{\partial x} \quad (17) $$

in terms of the local gradient of pressure $p$ and air dynamic viscosity $\mu$. In a bed of equally sided randomly packed spheres at the solid volume fraction $\nu$, the permeability $K$ may be approximated as

$$ \frac{K}{(1 - \nu)} \approx \frac{d^2}{150} \left( \frac{1 - \nu}{\nu} \right)^2, \quad (18) $$

and, because liquid water is not involved, the expression of $K$ need not be corrected for its presence (Brooks & Corey, 1964). In this work, we also ignore variations of $p$ with time and, because we did not record atmospheric pressure during this campaign, we assume $p \approx \rho_{st}$ without significant consequence for our results. Nonetheless, we consider pressure gradients in sand, which drive the seepage velocity $u$ through Equation 17.

Meanwhile, as derived by Bird et al. (2007) and Kays and Crawford, (1980), a balance of water vapor yields

$$ (1 - \nu) \frac{\partial \rho_{H_2O}}{\partial t} + (1 - \nu) \frac{\partial}{\partial x} \left[ \rho_{H_2O} u - \frac{\rho D}{\omega} \left( \frac{\partial Y}{\partial x} + a_T Y \frac{\partial \ln T}{\partial x} \right) \right] = - \rho_p \nu \frac{\partial \Omega}{\partial t}, \quad (19) $$
<figure><img src="image_5.png" /><figcaption>Figure 4. Snapshots of MovieDuneSubsurfaceTemperature&VaporMassFraction.mp4</figcaption></figure>

available as supporting information, showing subsurface depth profiles of temperature in °C (Eq. 10, red line) and vapor mass fraction $Y_e$ at equilibrium with the measured Ω (grey symbols) at 16:09 Qatar time ($t^* \simeq 1.67$) on March 20 (left) and 09:00 ($t^* \simeq 2.38$) on March 21 (right). The blue line is the model implemented in section 10. The yellow arrow points to co-located partial minima in the Y- and T-profiles. Bars and arrows represents scalars and vectors with lengths proportional to ambient temperature, wind speed, ambient water vapor mass fraction, dimensionless surface fluxes of vapor diffusion $\Psi_s^*(1-\nu)$ and vapor advection $u_s^*Y_s p^*/T_s^*$ defined in section 8. On the thin dotted line drawn for $x > 220$ mm, the vapor mass fraction Y is at equilibrium with a constant $\Omega = 1.65 \Omega_1$. The horizontal dotted line marks the sand surface. The horizontal grey and orange timeline indicates progression along the March 19-21 record.

scales, namely the diurnal period *J* that regulates the smooth evolution of *Y* and *T*, and a much shorter one that allowed the rapid propagation of waves carrying surface disturbances to the subsurface. Section 12 will develop a stability analysis explaining this behavior.

Above the dune, wind speed rose during the day as soon as a net radiation flux stroke the desert floor (Wyngaard, 2010). An unstable atmospheric boundary layer (ABL) then developed, in which ambient temperature grew with elevation, but remained lower than its counterpart on the surface. As the radiation flux turned upward before sunset, so did the gradient of ambient temperature. During the night, surface temperature was lower than the ambient, and both progressively decreased in tandem. Low-level jets intermittently perturbed the nocturnal stable boundary layer (Banta et al., 2007; Klein et al., 2016).

The absence of significant noise in our capacitance record allowed us to differentiate spatio-temporal variations of *Y* and use Fick's law to calculate the vapor mass flux $\dot{\Psi}_s'' = -\rho(1-\nu)(D/\varpi)(\partial Y/\partial x)_s$ at the sand surface, where *D* is the diffusion coefficient of water vapor at temperature *T* and $\varpi$ is the tortuosity of porous sand (Shen & Chen, 2007). As section 8 will show, this flux was augmented by a slow “seepage” advection within the pore space driven by wind blowing over topographical variations of the dune. Together, diffusion and advection fluxes provided a boundary condition for the ABL above our instruments.
where

$$ \rho_{H_2O} = \rho Y \qquad (20) $$

is the partial density of water vapor. In the bracket within the spatial derivative, ($\rho_{H_2O}u$) represents vapor advective, and $-(\rho D/\omega)\partial Y/\partial x$ is the Fickian diffusion flux, where $\omega$ is the “tortuosity” of the longer path taken by diffusion around grains (Shen & Chen, 2007). Because desert sands possess strong subsurface temperature gradients, we include the Soret diffusion term $-(\rho D/\omega)a_T Y \partial T/\partial x$ (Chapman & Cowling, 1953; Hudson et al., 2007), and we evaluate the coefficient $a_T$ in Appendix E. We will later see that this term is negligible. However, such diffusion may become important in extra-terrestrial environments.

Another gradient-driven process may be natural convection within the porous medium of permeability $K$ over a typical distance $\Delta x$, characterized by the Rayleigh number $Ra = Kg\Delta x(\Delta T/T)/[\alpha_m(\mu/\rho)]$ (Elder, 1967), where $g$ is gravitational acceleration and $\alpha_m \equiv k_s/(\rho c_p)$ combines the effective thermal conductivity $k_s \approx 0.49 \text{ W/m}\text{°K}$. Taking the highest observed temperature excursion $\Delta T < 20^\circ\text{C}$ over $\Delta x < 20 \text{ cm}$, we find $Ra < 0.006$, which is much too small for natural convection to matter.

Expanding derivatives in Equation 19 and combining with the overall mass conservation (15), we find

$$ \rho(1 - \nu)\frac{\partial Y}{\partial t} + \rho u(1 - \nu)\frac{\partial Y}{\partial x} - (1 - \nu)\frac{\partial}{\partial x}\left[\rho\frac{D}{\omega}\left(\frac{\partial Y}{\partial x} + a_T Y\frac{\partial \ln T}{\partial x}\right)\right] = -\rho_p\nu\frac{\partial \Omega}{\partial t}(1 - Y). \qquad (21) $$

In this equation, the first term on the left hand side captures the rate of change of water vapor mass fraction. The second term models seepage advection at the interstitial velocity $u$. The third term arises from molecular diffusion of water vapor in air with coefficient $D$.

Because temperature varies substantially through depth (Louge et al., 2013), it is important to account for variations of $D$ with $T$, and to keep it within the spatial derivative in Equation 21. Neglecting the dependence of $\rho$ on $Y$ in Equation 16,

$$ \rho D \approx \rho_{st} D_{st} T^{*n_d}, \qquad (22) $$

where $D_{st} \approx 2.576 \text{ 10}^{-5} \text{ m}^2/\text{s}$ and $n_d \approx 0.861$, close to the 1/2 exponent predicted by the kinetic theory of hard spheres (Chapman & Cowling, 1953).

The volumetric rate of water evaporation $\rho_p \nu \partial \Omega / \partial t$ is driven by the difference between the instantaneous vapor mass fraction $Y$ in the interstitial space of density $\rho(1 - \nu)$ surrounding sand grains, and the surface mass fraction $Y_e$ at equilibrium with the actual $\Omega$ that they hold,

$$ \rho_p \nu \frac{\partial \Omega}{\partial t} = \frac{\rho(1 - \nu)}{\tau} (Y - Y_e), \qquad (23) $$

In Section 11, we will derive this expression using the Hertz-Knudsen-Schrage (HKS) theory (Marek & Straub, 2001), and show that the characteristic time $\tau$ is governed by kinetics of an evaporation reaction at the grain scale. On the period $J = 24 \text{ hr}$ of diurnal variations, $\tau$ is $\ll J$, so grains have ample time to reach equilibrium. However, as Section 12 will discuss, rapid variations of $Y$ associated with aeolian transport may be too fast for particles to adjust instantly to a new equilibrium, leading to internal waves.

## 7. Equilibrium Regime

In general, there are two regimes arbitrated by the time scale $\delta t$ of $Y$, $T$, and $\Omega$ variations relative to $\tau$. If $\delta t \gg \tau$, all particles have reached equilibrium irrespective of size, so $\Omega \approx \Omega_e$. In this “equilibrium regime,” $\partial \Omega / \partial t \approx \partial \Omega_e / \partial t$ may be expressed in terms of partial derivatives of $T$, $Y$, and $p$ using functions $f_T$, $f_Y$, and $f_p$ derived in Appendix C,

$$ \frac{\partial \Omega_e}{\partial t} \equiv \Omega_1 \left[ \frac{f_T}{T_{st}} \frac{\partial T}{\partial t} + f_Y \frac{\partial Y}{\partial t} + \frac{f_p}{p_{st}} \frac{\partial p}{\partial t} \right], \qquad (24) $$
## 5 Stratigraphy

By recording both $K'_e$ and $K''_e$, the capacitance probe yielded not only $\Omega$ through $K''_e/K'_e$, but also the time-invariant local “bulk density” $\rho_b \simeq \rho_p\nu$ at each sensor, where $\nu$ is the solid volume fraction and $\rho_p$ the material density of sand grains. As Louge, Valance, el Moctar, and Dupont (2010) showed, $K'_e$ is approximated by the model of Böttcher and Bordewijk (1978),

$$K'_e = \frac{1}{4} \left\{ 2 - K_p + 3\nu(K_p - 1) + \sqrt{8K_p + [2 - K_p + 3\nu(K_p - 1)]^2} \right\}, \quad (12)$$

where $K_p$ is an equivalent material dielectric constant that rises with $\Omega$, and therefore with $\tan \varphi$. Rather than evaluating $K_p$ by staging several volume fractions in the laboratory, we found it more convenient to derive it directly from field data. In this method, we model $K_p$ as

$$K_p \simeq k_0 + k_1 \tan \varphi + k_2 \tan^2 \varphi + k_3 \tan^3 \varphi. \quad (13)$$

From Eq. (12), $K'_e$ and $K''_e = \tan \varphi \times K'_e$ become, respectively, the functions $f'$ and $f''$ of the variables ($\nu_i; \tan \varphi_j; k_0, k_1, k_2, k_3$) recorded at sensor of index $i$ and time $j$. The values of $k_0$ through $k_3$ and $\nu_1$ through $\nu_{15}$ are those yielding the smallest cost function $C$ from the entire data set of $N = 11557$ valid points acquired with all sensors,

$$C(\nu_1, \nu_2, \dots, \nu_{15}; k_0, k_1, k_2, k_3) = \sum_{j=1}^{N} \sum_{i=1}^{15} [K'_{e}(i, j) - f']^2 + [K''_{e}(i, j) - f'']^2. \quad (14)$$

Adding terms of order higher than $\tan^3 \varphi$ in Eq. (13) ($k_4, k_5$, etc) did not reduce $C$ and thus improve results.

As Fig. 5 shows, the first 14 sensors spanned a wide enough range of $\tan \varphi$ to discern variations of $\nu$ with depth. Typically, such variations are associated with layers sorted and deposited by successive avalanches on the downwind face of the mobile dune (Kocurek, 1991; Harari, 1996; Kleinhans, 2004; Bristow et al., 2007; Reesink & Bridge, 2009; Fan et al., 2012; Qian et al., 2014), then sequestered for several months until aeolian transport uncovers them upwind (Louge et al., 2013). On other visits to the same sand dune, the imprint of earlier positions of the lee face was visible for a short time after rain (Fig. 5C), suggesting preferential liquid penetration (Al-Shukaili et al., 2019). A complex layered stratigraphy was also evident upon excavation (Fig. 5A).

At the specific location of our measurements, Fig. 5 reveals relatively light sands with higher bulk density on the surface than below, consistent with the soft region that we selected for easy insertion of the probes. Its inset B also hints at an uneven stratigraphy on a wavelength $\Delta x \sim 3$ to $5$ cm. Because the upwind and downwind faces of the dunes made angles of $\theta_u \simeq 7^\circ$ and $\theta_d \simeq 31^\circ$ from the horizontal, our capacitance probe inserted approximately normal to the surface along $x$ detected strata separated by a distance $\Delta x \cos(\theta_u + \theta_d) \simeq 2$ to $4$ cm. Finally, although the profile in Fig. 5B resembled the layering shown in insets A and C, the Böttcher model of Eq. (12) likely underestimated $\nu$.

## 6 Modeling framework

Water evaporation from sands is complicated by capillary forces retaining liquid within its porous structure (Prat, 2002; Lehmann & Or, 2009; Or et al., 2013). Its analysis in the vadose zone therefore requires a model for the hysteretic process of wetting and draining in the unsaturated porous medium (Xu & Louge, 2015), typically obtained by invoking an empirical expression for the water retention curve (Vanderborght et al., 2017; Kamai & Assouline, 2018; Assouline & Kamai, 2019; Shao et al., 2021).

The process is simpler near the surface of hyper-arid deserts, where water exists either as a vapor, or as a film directly adsorbed on solid grains (Shahraeeni & Or, 2010).
where Equation 3 is the sand isotherm averaged over the particle size distribution. Because these functions all depend strongly on temperature, either explicitly or through $p_{\text{sat}}$, the term $\partial\Omega_e/\partial t$ couples spatial and temporal variations in $T$ and $Y$, ultimately leading to multiple inflexions in the $Y$-profile (Figure 4). With $\partial\Omega_e/\partial t \approx \partial\Omega_e/\partial t$, the right-hand source term in Equation 21 augments its unsteady counterpart to the left,

$$ \rho(1-\nu)\frac{\partial Y}{\partial t} + \rho_p\nu\frac{\partial\Omega}{\partial t}(1-Y) \\ \approx \rho(1-\nu)\left[1 + \frac{\rho_p\nu\Omega_1}{\rho(1-\nu)}f_Y(1-Y)\right]\frac{\partial Y}{\partial t} + \rho_p\nu\Omega_1(1-Y)\frac{f_T}{T_{\text{st}}}\frac{\partial T}{\partial t}. \quad (25) $$

Defining $\mathbb{R} \equiv \rho_p\nu/[\rho_{\text{st}}(1-\nu)]$ and $u^* \equiv u\sqrt{J/\alpha}$, recalling Equation 22 and $\rho_{\text{st}}/\rho \approx (T^*/p^*)$, dividing by $(1-\nu)\rho_{\text{st}}$ and reintroducing $t^* = t/J$, and $x^* = x/(\alpha J)^{1/2}$, Equation 21 has the dimensionless form

$$ \frac{p^*}{T^*} \left\{ \left[ 1 + \mathbb{R}\Omega_1 \frac{T^*}{p^*} f_Y(1-Y) \right] \frac{\partial Y}{\partial t^*} + \mathbb{R}\Omega_1 \frac{T^*}{p^*} f_T(1-Y) \frac{\partial T^*}{\partial t^*} \right\} + \frac{p^*}{T^*} u^* \frac{\partial Y}{\partial x^*} - \mathbb{L} \frac{\partial}{\partial x^*} \left\{ T^{*n_d} \left( \frac{\partial Y}{\partial x^*} + a_T Y \frac{\partial \ln T^*}{\partial x^*} \right) \right\} = 0, \quad (26) $$

where $\mathbb{L} \equiv D_{\text{st}}/(\varpi\alpha)$ is a hybrid Lewis number comparing interstitial vapor diffusivity with thermal diffusivity through bulk sand. With a tortuosity estimate for spheres $\varpi \sim \pi/2$ (Shen & Chen, 2007), $\mathbb{L} \approx 46.5$. With $\rho_p \approx 2,630$ kg/m$^3$ and $\nu \approx 0.545$, $\mathbb{R} \approx 2660$. For the entire data set, $-51 < f_T < -1.2$ and $19 < f_Y < 210$. Therefore, in the unsteady term in straight brackets of Equation 26, the expression

$$ \mathbb{R}\Omega_1 \frac{T^*}{p^*} f_Y(1-Y) \approx \mathbb{R}\Omega_1 \frac{T^*}{p^*} f_Y \gg 1 \quad (27) $$

is very large. To illustrate the significance of this inequality, consider an isobaric and isothermal system without seepage at $Y \ll 1$. In this hypothetical situation, the dimensional version of Equation 26 could be approximated as

$$ \frac{\partial Y}{\partial t} \approx \left[ \frac{D/\varpi}{\mathbb{R}\Omega_1 (T^*/p^*) f_Y} \right] \frac{\partial^2 Y}{\partial x^2}. \quad (28) $$

If such system, unlike sand, also exhibited a linear isotherm $\Omega = \Omega_{p_0} + \Omega_{p_1} \times \text{RH}$, as is the case in pharmaceutical powders (Klinzing & Zavaliangos, 2016; Louge et al., 2021), which absorb, rather than adsorb water, differentiating the isotherm (Appendix C) would yield an approximately invariant $f_Y = \partial\text{RH}/\partial Y \approx pM^*/p_{\text{sat}} \gg 1$, and Equation 28 would govern an unsteady diffusion with apparent coefficient $D_{\text{eff}} \approx (D/\varpi)/(\mathbb{R}\Omega_{p_1}(p_{\text{st}}/p_{\text{sat}})T^*M^*)$, much smaller than its counterpart without water adsorption on grains (or absorption within them). In that case, the mass fraction in a powder initially at $Y = Y_i$ subject to $Y = Y_s$ at the surface of a half-space would evolve as $Y = Y_s - (Y_s - Y_i)\text{erf}(x/2\sqrt{D_{\text{eff}}t})$. In short, grain adsorption (or absorption) would result in much slower diffusive behavior.

For sands, the isotherm in Equation 3 is more complicated, with $f_Y \propto 1/Y$, thereby precluding a simple closed-form solution. However, on long time scales $\delta_i \gg \tau$, the source term $\partial\Omega_e/\partial t$ also acts as an inhibitor of diffusion, as successive layers in an initially dehydrated sand must adsorb water before letting the vapor diffuse deeper. If we ignored the subtlety of $f_Y \propto 1/Y$ for our own data set, the range of $f_Y$ would roughly imply that the apparent $D_{\text{eff}}$ is 70–730 times weaker than $D$. Similar observations are found in moisture-absorbing pharmaceutical powders (Louge et al., 2021). Louge, Valance, Babah, et al. (2010) did not consider such hindered diffusion. However, because they derived the steady distribution of moisture through sand ripples, their model needs no repair, unless sand is exposed to rapid variations in surface moisture.

While temperature variations described in Section 3 conformed to the linear PDE (8), the corresponding model of vapor mass fraction in Equation 26 is highly non-linear, as it strongly couples $Y$ to subsurface variations in $T$. Another complication is a variable advection flow within the dune subsurface, which we consider next.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Real part of the effective dielectric constant $K'_e$ vs loss tangent $tan\phi$ for all valid data points (small circles), shown for sensor 4, 9 and 13. The line is the Böttcher model in Eq. (12) obtained with coefficients ($k_0, k_1, k_2, k_3$) $\simeq (6.7, 20, -7.2, 3.9)$ from Eq. (13) and volume fractions at all 15 sensors that minimize the cost function in Eq. (14). Dashed lines above and below the main fit, as well as error bars in inset B, mark $\pm$ one sample standard deviation in $\nu$ at the corresponding sensor. Inset A: an excavation conducted in April 2014 reveals a layered stratigraphy after pouring water laced with blue food coloring through the surface. The stratigraphy occasionally appears on the surface after rain partially dries up (inset C, picture taken on the same dune in January 2014). Inset B: profile of $\nu$ vs depth $x$ (cm) perpendicular to the dune surface, calculated from Eqs. (12)–(14). The dashed line is a visual spline through the data.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Snapshots of Movie S2. Left: best fits (lines) of recorded available data (symbols) of Y versus x* to the seepage velocity u* using ordinary differential equation (ODE) (30) at 08:03 on 20 March. Dotted lines mark Y ± 7% Y, from which we estimate uncertainties in u*. The integration of ODE (30) without seepage, shown as a dashed line, fails to capture data. Blue, yellow, and red vectors are proportional to dimensionless fluxes of, respectively, advection u_s*Y_sP*/T_s*, diffusion Ψ_s*/(1 - v) from Equation 32, and their sum at the dune surface. Later in the movie (not shown here), the signal splits at 13:25 on 20 March in two bistable envelopes analyzed in Section 12. During the subsequent period of high winds, dryer conditions prevail with aeolian transport (pink) than without (blue). Right: measured unsteady terms U_M (blue line), U_T (red line), and U_M + U_T (black dashed line) defined in Equation 29.</figcaption></figure>

8. Seepage

It would not be straightforward to measure the small interstitial velocity *u* in the field. However, the recorded mass fraction *Y* (Figure 3) exhibits so little noise that its temporal derivative can be reliably evaluated using a smoothing spline. Then, it is possible to fit the value of *u*\* at each recorded time step by transforming the governing PDE (26) into an ordinary differential equation (ODE) where the unsteady term *U*<sub>M</sub> + *U*<sub>T</sub> is calculated from data,

$$
\begin{align}
U_M &\equiv \left[ \frac{p^*}{T^*} + \mathbb{R}\Omega_1 f_Y(1-Y) \right] \frac{\partial Y}{\partial t^*} \approx \mathbb{R}\Omega_1 f_Y(1-Y) \frac{\partial Y}{\partial t^*}, \tag{29} \\
U_T &\equiv \mathbb{R}\Omega_1 f_T(1-Y) \frac{\partial T^*}{\partial t^*}. \nonumber
\end{align}
$$

Rearranging Equation 26 and replacing $\partial^2 T^*/\partial x^{*2}$ by $\partial T^*/\partial t^*$ using the heat equation, we find the ODE

$$
\begin{equation}
\begin{split}
& \frac{d^2 Y}{dx^{*2}} + \frac{dY}{dx^*} \left[ n_d \frac{\partial \ln T^*}{\partial x^*} - \frac{p^* u^*}{\mathbb{L} T^{*n_d} + 1} \right] \\
& \qquad + a_T Y \left[ \frac{d \ln Y}{dx^*} \times \frac{\partial \ln T^*}{\partial x^*} + \frac{\partial \ln T^*}{\partial t^*} + (n_d - 1) \left( \frac{\partial \ln T^*}{\partial x^*} \right)^2 \right] = \frac{U_M + U_T}{\mathbb{L} T^{*n_d}},
\end{split}
\tag{30}
\end{equation}
$$

which we integrate subject to the measured *Y* at depths *x*<sub>1</sub><sup>*</sup> and *x*<sub>10</sub><sup>*</sup> of sensors 1 and 10 using the two-point boundary value toolbox bvp4c of MATLAB. For consistency, we limit the integration to the interval [*x*<sub>1</sub><sup>*</sup>, *x*<sub>10</sub><sup>*</sup>], since moisture on sensors 11 to 15 is sometimes too large for the probe to return a stable measurement of *Y* (Figure 3). Using the simplex search method of Lagarias et al. (1998), which is implemented as fminsearch in MATLAB, we repeat the integration of ODE (30) for several values of *u*\* until the root-mean-square (RMS) difference between signal and prediction is minimized.

Figure 6 illustrates this procedure. It is best viewed in the animation Movie S2. In this figure, the left graph shows the profile of least RMS error, as well as fits of *u*\* to ±7% excursions in *Y*, which we use as uncertainty estimates in *u*\*. These excursions enclose spatial oscillations of relatively small amplitude in *Y*, which may be associated with the mild stratigraphy shown in Figure 5. The corresponding uncertainties in *u*\* (dashed lines in Figure 7) are
Here, the drying front that Shokri et al. (2008) and Lehmann et al. (2008) documented lies deeper within the dune, and mass transfer near the surface conforms to the last stage of the evaporation process (Shahraeeni et al., 2012).

Such was the case for the data set of section 4, which suggested a strong spatio-temporal coupling between subsurface variations of vapor mass fraction and temperature, a situation that is difficult to reproduce in laboratory experiments on the drying process (Davarzani et al., 2014; Gao et al., 2018; Z. Li et al., 2020).

In this section, we derive a set of governing equations to interpret these observations in the limit where interstitial water vapor is at equilibrium with water adsorbed on solid grains. We will use their framework to identify the principal mechanisms of moisture transport in the dune subsurface, and dismiss those that are negligible.

In the gas phase, the balance of total mass in an elementary slice of sand at depth $x$ is

$$ (1 - \nu) \frac{\partial \rho}{\partial t} + (1 - \nu) \frac{\partial (\rho u)}{\partial x} = - \rho_p \nu \frac{\partial \Omega}{\partial t}, \qquad (15) $$

where the right-hand term represents the rate of evaporation from grains of material density $\rho_p$ on which water was adsorbed, later considered in section 11 and Appendix G. The variables $\rho$, $u$, $\Omega$ and $\nu$ in this classical mass balance (Kays & Crawford, 1980; Bird et al., 2007) are implicitly coarse-grained on the scale of a few particles, as suggested by T. B. Anderson and Jackson (1967). For clarity of derivations, we ignore relatively small variations of $\nu$ mentioned in section 5, so $(1 - \nu)$ can be brought out of the spatial derivative. However, prior knowledge of a stratigraphy like that shown in Fig. 5 can be substituted in the governing equations and its numerical solutions.

In Eq. (15), the density of humid air $\rho = \rho_{H_2O} + \rho_{air}$ is the sum of the partial densities of water vapor $\rho_{H_2O} \equiv \rho Y = MW_{H_2O}\chi_{H_2O}(p/\hat{RT})$ and dry air $\rho_{air} = MW_{air}\chi_{air}(p/\hat{RT})$, with complementary mole fractions of water and air $\chi_{H_2O} + \chi_{air} = 1$, so vapor mass fraction is $Y = \chi_{H_2O}MW_{H_2O}/[\chi_{H_2O}MW_{H_2O} + \chi_{air}MW_{air}]$. Meanwhile, dry air density $\rho_{air} = \rho_{st}p^*/T^*$ can be expressed in terms of its counterpart $\rho_{st} \approx 1.185 \text{ kg/m}^3$ under standard conditions at $T_{st} = 298.15 \text{ K}$ and $p_{st} = 1.013 \times 10^5 \text{ Pa}$, where $p^* \equiv p/p_{st}$ and $T^* \equiv T/T_{st}$ are dimensionless absolute pressure and temperature, respectively. Using these relations, moist air density is

$$ \rho = \rho_{st} \frac{p^*/T^*}{1 + Y(M^* - 1)}, \qquad (16) $$

where $M^* \equiv MW_{air}/MW_{H_2O}$. Because $Y(M^* - 1) \ll 1$, $\rho$ is almost equal to the density of air at the same $p$ and $T$. Therefore, although numerical solutions will adopt Eq. (16), we will lighten derivations in this article by simply writing $\rho \approx \rho_{st} p^*/T^*$.

The “seepage” interstitial gas velocity $u$ is given by Darcy’s law

$$ u = - \frac{K}{\mu(1 - \nu)} \frac{\partial p}{\partial x} \qquad (17) $$

in terms of the local gradient of pressure $p$ and air dynamic viscosity $\mu$. In a bed of equally-sided randomly-packed spheres at the solid volume fraction $\nu$, the permeability $K$ may be approximated as

$$ \frac{K}{(1 - \nu)} \approx \frac{d^2}{150} \left( \frac{1 - \nu}{\nu} \right)^2, \qquad (18) $$

and, because liquid water is not involved, the expression of $K$ need not be corrected for its presence (Brooks & Corey, 1964). In this work, we also ignore variations of $p$ with time and, because we did not record atmospheric pressure during this campaign, we assume $p \approx p_{st}$ without significant consequence for our results. Nonetheless, we consider pressure gradients in sand, which drive the seepage velocity $u$ through Eq. (17).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Time-history of dimensionless seepage velocity $u^*$ at the dune surface versus $t^*$ (solid line) with surrounding uncertainties (dashed lines) associated with the measured $Y \pm 7\% Y$ shown in Figure 6. The gray regions enclose values of $u^*$ calculated in Appendix F for depths $x^* < 0.87$ above sensor 10. Splitting of the $u^*$-record in the period $1.4 < t^* < 1.63$ arises as signal envelopes separate during aeolian transport. Solid lines in insets (a–c) are typical profiles of $u^* - u_s^*$ calculated through these depths at times $t^* = 0.98, 1.54$ and $2.47$; dashed lines illustrate the role of vapor release on $u^*$ by artificially making $f_Y = f_Y = 0$ in these calculations. Insets and gray regions reveal that $u^*$ varies little with depth in this case.</figcaption></figure>

small. The integration of ODE (30) fails with $u^* = 0$ (dashed lines), thereby confirming that seepage advection is essential.

Shao et al. (2021) agreed with these observations. They suggested that a small advection on the order of 100 µm/s is necessary to capture data, and that thermal processes are crucial. However, by ignoring the soil's isotherm, Shao et al. (2021) had to adopt heuristic closures for the local evaporation rate, and they could not relate the apparent diffusion coefficient $D_{eff}$ to its known counterpart $D$ for water vapor in air. In addition, the spatial resolution and accuracy of their TDR was not sufficient to discern subtle variations in equilibrium vapor mass fraction that are needed to evaluate seepage time-history. They reported $5.4 \times 10^{-3} < \theta_w < 8.2 \times 10^{-3}$, where $\theta_w \approx (\rho_p \nu / \rho_w) \Omega$ is water volume fraction. Combining Equations 3 and 7, and using properties of our sands as a rough estimate, their smallest detected $\theta_{w_{min}} \approx 5.4 \times 10^{-3}$ translates to a vapor mass fraction $Y_{min} \approx 0.02$, which is well over an order of magnitude larger than what our own capacitance probe discerned ($Y_{min} \approx 0.0012$).

The right graph of Figure 6 also shows relative magnitudes of the two contributions $\cup_M$ and $\cup_T$ to the total unsteady term. At long diurnal time scales for which grains achieve equilibrium with the surrounding $Y$, these terms are dominated by the rate $\partial\Omega/\partial t^*$ of the moisture exchange between solid and gas, ultimately captured by the functions $f_T$ and $f_Y$ derived in Appendix C. Because $\cup_M$ and $\cup_T$ often have opposite signs but similar magnitudes, temporal variations of mass fraction and temperature are both crucial, and they compete to govern vapor transfer through the dune subsurface.

A subtlety is that the seepage velocity $u^*$ varies along $x^*$ in response to density changes arising from temperature, pressure or the release of water vapor. Appendix F provides the corresponding ODEs (F1)–(F2), which we solve using MATLAB's ODE45 4-th order Runge-Kutta algorithm, subject to dimensionless seepage velocity $u^* = u_s^*$ and pressure $p^* = p_s^*$ at the surface (subscript $s$), once $Y$ and its derivatives have been obtained from the solution of Equation 30.

Insets in Figure 7 show the resulting profiles as solid lines, while dashed lines indicate calculations with $f_Y = f_T = 0$. Contrasting predictions with and without these functions underscores how the rate of water vapor release $\partial\Omega/\partial t^*$ affects $u^*$. Under planetary conditions where ambient pressure is small or gradients are large, ODEs (30) and (F1)–(F2) may be coupled strongly enough that the evolution of $Y$, $\rho$, and $u^*$ with $x^*$ would require their simultaneous solution. However, for this dune, depth variations in $u^*$ (solid lines) are small, and so it is legitimate to assume $u^* \approx u_s^*$. A consequence is that seepage is not directly affected by internal gradients of temperature or mass fraction. Instead, as the next section articulates, it is likely due to wind-driven variations of surface pressure.

In short, integration of ODE (30) yields our best estimate of the dimensionless instantaneous velocity

$$v_s^* = +(1 - \nu)u_s^* \quad (31)$$

counted >0 into sand from the ABL. By fitting the recorded depth profile of $Y$ with a smoothing spline and differentiating the result at the surface, we also evaluate the diffusive vapor flux $\Psi_s'' = -\rho(1 - \nu)(D/\sigma)(\partial Y / \partial x)_s$ across the surface, counted >0 into sand. Made dimensionless with $\rho_{st} \sqrt{\alpha/J}$, this flux

$$\Psi_s^* = -(1 - \nu)\mathbb{L}T_s^{*nd}(\partial Y / \partial x)_s, \quad (32)$$

augmented by the advection $u_s^* Y_s p^* / T_s^*$, contributes to the surface BC for the ABL, which we will examine in Section 13.

## 9. Origin of Seepage

Louge, Valance, Babah, et al. (2010) recognized that wind-driven seepage advects moisture and dust through sand ripples by letting the Bernoulli effect induce a higher static pressure at their troughs than at their crests. Gao et al. (2018) staged a wavy surface to investigate its role in drying. By deploying a porous plastic ripple with
Meanwhile, as derived by (Kays & Crawford, 1980; Bird et al., 2007), a balance of water vapor yields

$$ (1 - \nu) \frac{\partial \rho_{\text{H}_2\text{O}}}{\partial t} + (1 - \nu) \frac{\partial}{\partial x} \left[ \rho_{\text{H}_2\text{O}} u - \frac{\rho D}{\varpi} \left( \frac{\partial Y}{\partial x} + a_T Y \frac{\partial \ln T}{\partial x} \right) \right] = - \rho_p \nu \frac{\partial \Omega}{\partial t}, \quad (19) $$

where

$$ \rho_{\text{H}_2\text{O}} = \rho Y \quad (20) $$

is the partial density of water vapor. In the bracket within the spatial derivative, ($\rho_{H_2O}u$) represents vapor advection, and ($\rho D/\varpi \partial Y/\partial x$) is the Fickian diffusion flux, where $\varpi$ is the “tortuosity” of the longer path taken by diffusion around grains (Shen & Chen, 2007). Because desert sands possess strong subsurface temperature gradients, we include the Soret diffusion term ($\rho D/\varpi a_T Y \partial \ln T/\partial x$) (Chapman & Cowling, 1953; Hudson et al., 2007), and we evaluate the coefficient $a_T$ in Appendix E. We will later see that this term is negligible. However, such diffusion may become important in extra-terrestrial environments.

Another gradient-driven process may be natural convection within the porous medium of permeability $K$ over a typical distance $\Delta x$, characterized by the Rayleigh number $Ra = Kg\Delta x(\Delta T/T)/[\alpha_m(\mu/\rho)]$ (Elder, 1967), where $g$ is gravitational acceleration and $\alpha_m \equiv k_s/(\rho c_p)$ combines the effective thermal conductivity $k_s \approx 0.49 \text{ W/m}\cdot{}^{\circ}\text{K}$ of the sand bed with the density of air and its specific heat $c_p \approx 1000 \text{ J/kg}\cdot{}^{\circ}\text{K}$. Taking the highest observed temperature excursion $\Delta T < 20 \text{ }^{\circ}\text{C}$ over $\Delta x < 20 \text{ cm}$, we find $Ra < 0.006$, which is much too small for natural convection to matter.

Expanding derivatives in Eq. (19) and combining with the overall mass conservation (15), we find

$$ \rho(1-\nu)\frac{\partial Y}{\partial t} + \rho u(1-\nu)\frac{\partial Y}{\partial x} - (1-\nu)\frac{\partial}{\partial x}\left[\rho\frac{D}{\varpi}\left(\frac{\partial Y}{\partial x} + a_T Y \frac{\partial \ln T}{\partial x}\right)\right] = -\rho_p \nu \frac{\partial \Omega}{\partial t}(1-Y). \quad (21) $$

In this equation, the first term on the left hand side captures the rate of change of water vapor mass fraction. The second term models seepage advection at the interstitial velocity $u$. The third term arises from molecular diffusion of water vapor in air with coefficient $D$.

Because temperature varies substantially through depth (Louge et al., 2013), it is important to account for variations of $D$ with $T$, and to keep it within the spatial derivative in Eq. (21). Neglecting the dependence of $\rho$ on $Y$ in Eq. (16),

$$ \rho D \approx \rho_{\text{st}} D_{\text{st}} T^{*n_d}, \quad (22) $$

where $D_{\text{st}} \approx 2.576 \times 10^{-5} \text{ m}^2/\text{s}$ and $n_d \approx 0.861$, close to the 1/2 exponent predicted by the kinetic theory of hard spheres (Chapman & Cowling, 1953).

The volumetric rate of water evaporation $\rho_p \nu \partial \Omega / \partial t$ is driven by the difference between the instantaneous vapor mass fraction $Y$ in the interstitial space of density $\rho(1-\nu)$ surrounding sand grains, and the surface mass fraction $Y_e$ at equilibrium with the actual $\Omega$ that they hold,

$$ \rho_p \nu \frac{\partial \Omega}{\partial t} = \frac{\rho(1-\nu)}{\tau} (Y - Y_e), \quad (23) $$

In section 11, we will derive this expression using the Hertz-Knudsen-Schrage (HKS) theory (Marek & Straub, 2001), and show that the characteristic time $\tau$ is governed by kinetics of an evaporation reaction at the grain scale. On the period $J = 24 \text{ hr of diurnal variations}$, $\tau$ is $\ll J$, so grains have ample time to reach equilibrium. However, as section 12 will discuss, rapid variations of $Y$ associated with aeolian transport may be too fast for particles to adjust instantly to a new equilibrium, leading to internal waves.
harmonic amplitude $h_0$ and wavelength $\lambda_0$ in a wind tunnel, Musa et al. (2014) confirmed that the pore pressure gradient driving this advection scales as $\rho U^2 h_0 / \lambda_0^2$, where U is the bulk wind speed. With origin above the trough at mid-depth between trough and crest, they showed that a ripple with ondulations (counted positive downward like x)

$$h = h_0 \cos(\omega_{x_0} \mathcal{X}) \qquad (33)$$

along the wind direction of unit vector $\hat{\mathbf{X}}$ produces a pore pressure field in the subsurface half-space that is a solution of the Laplace equation governing quasi-steady Darcy flow,

$$p = p_a + \rho v_r^2 \times \delta p^* \exp(-\omega_{x_0}x), \qquad (34)$$

where $\omega_{x_0} = 2\pi/\lambda_0$ is the single-mode wavenumber of the harmonic ripple along $\mathcal{X}$, $v_r$ is the turbulent shear velocity, and

$$\delta p^* = 2\pi (h_0/\lambda_0) \sqrt{C^2 + D^2} \cos[(\omega_{x_0} \mathcal{X}) - \phi_0] \qquad (35)$$

is the dimensionless excursion from the ambient pressure $p_a$ at the surface, with phase lag $\phi_0 = \arctan(D/C) > 0$, Fourrière (2009) and Fourrière et al. (2010) calculated the coefficients C and D in terms of the parameter $\ln(2\pi z_0 / \lambda_0)$, where $z_0$ is the hydrodynamic roughness (Claudin et al., 2016). Through Darcy's Equation 17, the gradient of the pressure field in Equation 34 then induces the seepage velocity at $x = 0$,

$$u_s = U^2 \left[ \frac{K \rho (v_r / U)^2}{\mu (1 - v)} \right] \omega_{x_0} \delta p^* \equiv U^2 \tau_K \omega_{x_0} \delta p^*, \qquad (36)$$

which, like surface pressure, features a peak lagging behind the harmonic ripple profile by a distance $\lambda_0 \phi_0 / (2\pi)$ along the wind. In Equation 36, $\tau_K \equiv K \rho (v_r / U)^2 / [\mu(1 - v)]$ has units of time and incorporates the ratio $v_r / U$, which typically varies on the sand surface.

For our dune, we estimate $\lambda_0 \approx 220$ m as approximately four times its toe to crest distance, and $h_0 \approx 4.5$ m as its crest height (Figure 8). Anemometry carried out near the capacitance probe location in January 2017 also found $z_0 \approx 16$ µm. With $\ln(2\pi z_0 / \lambda_0) \approx -14.6$, Fourrière (2009) calculated $C \approx -1$, 150 and $D \approx 49$. Because C and D has a weak logarithmic dependence on $z_0$ and modal wavelengths, we ignore variations of these quantities on the dune and, for simplicity, treat C and D as constants.

Because the single-mode wavenumber $\omega_{x_0} = 2\pi/\lambda_0$ is very small, Equation 36 would predict a much smaller seepage than a ripple of comparable aspect ratio, if the dune possessed a harmonic profile. However, the dune shape is not harmonic, but instead features gentle undulations in its surface elevation $z(\mathcal{X}, \mathcal{Y})$ along the wind unit direction $\hat{\mathbf{X}}$ and its cross-wind counterpart $\hat{\mathbf{Y}}$. The small ratio of dune height on wavelength $h_0 / \lambda_0 \sim 0.02$ further suggests that $z(\mathcal{X}, \mathcal{Y})$ may be expanded in Fourier modes that contribute linearly to the overall pressure excursion (Musa et al., 2014). Then, as the following analysis shows, some of these modes may possess large enough wave-numbers to induce significant seepage.

Inspired by the theory of P. S. Jackson and Hunt (1975) and Hunt et al. (1988), Kroy et al. (2002) derived the transfer function relating the two-dimensional Fourier transforms of dune elevation $\tilde{z}$ and perturbation in turbulent shear stress $\delta\tilde{\tau}$ in terms of two parameters A and B. The dimensionless pressure excursion $\delta p^*$ is subject to a similar transfer function, in which A and B are respectively replaced by C and D,

$$\frac{\delta \tilde{p}^*}{\tilde{z}} = \frac{\omega_X (C \omega_X + iD |\omega_X|)}{\sqrt{\omega_X^2 + \omega_Y^2}}, \qquad (37)$$

where $\omega_X \equiv 2\pi/\lambda_X$ and $\omega_Y \equiv 2\pi/\lambda_Y$ are wavenumbers along $\hat{\mathbf{X}}$ and $\hat{\mathbf{Y}}$, respectively. With Equation 37, the two-dimensional Fourier transform of $z$ obtained from a dune survey therefore yields $\delta \tilde{p}^*$, which is then inverted to recover the pressure excursion everywhere on the dune surface. To evaluate $\tilde{z}$ at each time step, we tilt the dune slightly so it rises above an artificially flattened desert floor, rotate it around the origin of $(\mathcal{X}, \mathcal{Y})$ at the barycenter of its outline to align $\hat{\mathbf{X}}$ with the instantaneous wind direction, interpolate survey points using the C1-continuous scattered Interpolant of MATLAB, subdivide the square central domain circumscribing the rotated
7 Equilibrium regime

In general, there are two regimes arbitrated by the time scale $\delta t$ of Y, T and $\Omega$ variations relative to $\tau$. If $\delta t \gg \tau$, all particles have reached equilibrium irrespective of size, so $\Omega \simeq \Omega_e$. In this “equilibrium regime”, $\partial\Omega/\partial t \simeq \partial\Omega_e/\partial t$ may be expressed in terms of partial derivatives of T, Y and p using functions $f_T$, $f_Y$ and $f_p$ derived in Appendix C,

$$ \frac{\partial \Omega_e}{\partial t} = \Omega_1 \left[ \frac{f_T}{T_{\text{st}}} \frac{\partial T}{\partial t} + f_Y \frac{\partial Y}{\partial t} + \frac{f_p}{p_{\text{st}}} \frac{\partial p}{\partial t} \right], \quad (24) $$

where Eq. (3) is the sand isotherm averaged over the particle size distribution. Because these functions all depend strongly on temperature, either explicitly or through $p_{\text{sat}}$, the term $\partial\Omega_e/\partial t$ couples spatial and temporal variations in T and Y, ultimately leading to multiple inflexions in the Y-profile (Fig. 4). With $\partial\Omega/\partial t \simeq \partial\Omega_e/\partial t$, the right-hand source term in Eq. (21) augments its unsteady counterpart to the left,

$$ \begin{aligned} \rho(1-\nu)\frac{\partial Y}{\partial t} &+ \rho_p\nu\frac{\partial \Omega}{\partial t}(1-Y) & (25) \\ &\approx \rho(1-\nu)\left[1+\frac{\rho_p\nu\Omega_1}{\rho(1-\nu)}f_Y(1-Y)\right]\frac{\partial Y}{\partial t} + \rho_p\nu\Omega_1(1-Y)\frac{f_T}{T_{\text{st}}}\frac{\partial T}{\partial t}. \end{aligned} $$

Defining $\mathbb{R} \equiv \rho_p\nu/[\rho_{\text{st}}(1-\nu)]$ and $u^* \equiv u\sqrt{J/\alpha}$, recalling Eq. (22) and $\rho_{\text{st}}/\rho \simeq (T^*/p^*)$, dividing by $(1-\nu)\rho_{\text{st}}$ and reintroducing $t^* = t/J$, and $x^* = x/(\alpha J)^{1/2}$, Eq. (21) has the dimensionless form

$$ \frac{p^*}{T^*} \left\{ \left[ 1 + \mathbb{R}\Omega_1 \frac{T^*}{p^*} f_Y(1-Y) \right] \frac{\partial Y}{\partial t^*} + \mathbb{R}\Omega_1 \frac{T^*}{p^*} f_T(1-Y) \frac{\partial T^*}{\partial t^*} \right\} + \frac{p^*}{T^*} u^* \frac{\partial Y}{\partial x^*} - \mathbb{L} \frac{\partial}{\partial x^*} \left\{ T^{*n_d} \left( \frac{\partial Y}{\partial x^*} + a_T Y \frac{\partial \ln T^*}{\partial x^*} \right) \right\} = 0, \quad (26) $$

where $\mathbb{L} \equiv D_{\text{st}}/(\varpi\alpha)$ is a hybrid Lewis number comparing interstitial vapor diffusivity with thermal diffusivity through bulk sand. With a tortuosity estimate for spheres $\varpi \sim \pi/2$ (Shen & Chen, 2007), $\mathbb{L} \simeq 46.5$. With $\rho_p \simeq 2630 \text{ kg/m}^3$ and $\nu \simeq 0.545$, $\mathbb{R} \simeq 2660$. For the entire data set, $-51 < f_T < -1.2$ and $19 < f_Y < 210$. Therefore, in the unsteady term in straight brackets of Eq. (26), the expression

$$ \mathbb{R}\Omega_1 \frac{T^*}{p^*} f_Y(1-Y) \simeq \mathbb{R}\Omega_1 \frac{T^*}{p^*} f_Y \gg 1 \quad (27) $$

is very large. To illustrate the significance of this inequality, consider an isobaric and isothermal system without seepage at $Y \ll 1$. In this hypothetical situation, the dimensional version of Eq. (26) could be approximated as

$$ \frac{\partial Y}{\partial t} \approx \left[ \frac{D/\varpi}{\mathbb{R}\Omega_1(T^*/p^*)f_Y} \right] \frac{\partial^2 Y}{\partial x^2}. \quad (28) $$

If such system, unlike sand, also exhibited a linear isotherm $\Omega = \Omega_{p0}+\Omega_{p1}\times\text{RH}$, as is the case in pharmaceutical powders (Louge et al., 2021), which absorb, rather than adsorb water, differentiating the isotherm (Appendix C) would yield an approximately invariant $f_Y = \partial\text{RH}/\partial Y \simeq pM^*/p_{\text{sat}} \gg 1$, and Eq. (28) would govern an unsteady diffusion with apparent coefficient $D_{\text{eff}} \simeq (D/\varpi)/[\mathbb{R}\Omega_{p1}(p_{\text{st}}/p_{\text{sat}})T^*M^*]$, much smaller than its counterpart without water adsorption on grains (or absorption within them). In that case, the mass fraction in a powder initially at $Y = Y_i$ subject to $Y = Y_s$ at the surface of a half-space would evolve as $Y = Y_s-(Y_s-Y_i)\text{erf}(x/2\sqrt{D_{\text{eff}}t})$. In short, grain adsorption (or absorption) would result in much slower diffusive behavior.

For sands, the isotherm in Eq. (3) is more complicated, with $f_Y \propto 1/Y$, thereby precluding a simple closed-form solution. However, on long time scales $\delta_t \gg \tau$, the source term $\partial\Omega_e/\partial t$ also acts as an inhibitor of diffusion, as successive layers in an initially dehydrated sand must adsorb water before letting the vapor diffuse deeper. If we ignored
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Seepage predictions of Equation 39 calculated for the dune shape in January 2017, the wind record *U* in Figure 3, and the instantaneous wind direction $\hat{\mathbf{X}}$. (a) Relative surface seepage velocity $u_s/\tau_K$ (m/s² > 0 into sand) at the position shown by the gray dot in inset (e) versus diurnal time *t*/*J*. (b) Instantaneous wind bearing $\theta$ (°) versus *t*/*J*, and mean historical bearing from Louge et al. (2013) (horizontal line). In insets (a and b), the line is a visual fit, and vertical solid and dashed lines mark sunrise and sunset, respectively. (c) Solid line: predicted relative pressure excursion $\delta p^*$ versus distance $\mathcal{X}$ (m) along the transect marked by the dashed line in inset (e); dashed line: dune elevation *z* (m) along this transect. (d) Dimensionless spectral power $a_i^2 + b_i^2$ of pressure excursions in Equation 38 versus wavelength $\lambda_i = \Delta\mathcal{X}/i$ with $i = 1, \dots, 2^{N+1}$ for the wind bearings shown. (e) Snapshot of $u_s/\tau_K$ on the dune surface at 11:14 on 21 March, marked by gray dots in insets (a and b). For the whole sequence, see Movie S3. Black isolines are placed every 0.5 m. For best contrast, colors are obtained by converting $[1 + \tanh(u_s/(\tau_K\gamma_0))]/2$ to the jet color map of MATLAB with $\gamma_0 = 10^3$ m/s².</figcaption></figure>

dune outline in $2^N$ cells with $N = 9$, and buttress this domain with eight identical squares holding zero elevation to avoid aliasing. The resulting physical domain of width $\Delta\mathcal{X}$ therefore holds $2^{N+2}$ cells along both directions, leading to $2^{N+2}$ discrete wavenumbers $\omega_X = (2\pi/\Delta\mathcal{X}) \times (0, \dots, 2^{N+1}, -2^{N+1} + 1, \dots, -1)$, organized in the peculiar order with which MATLAB optimizes its discrete Fourier transforms (DFT). Comparing predictions of Equation 37 from the two-dimensional fft2 DFT and its one-dimensional counterpart fft with $\omega_Y = 0$, we find that transverse wavenumbers $\omega_Y$ play a minor role. Therefore, for simplicity, we reconstruct the pressure excursion $\delta p^*$ on any wind-directed transect using the first $1 + 2^{N+1}$ coefficients of the one-dimensional DFT $\delta\tilde{p}^*$ in Equation 37 with $\omega_Y = 0$,

$$ \delta p^* = a_0 + \sum_{i=1}^{2^{N+1}} \delta p_i^* = a_0 + \sum_{i=1}^{2^{N+1}} a_i \cos(2\pi\mathcal{X}/\lambda_i) + b_i \sin(2\pi\mathcal{X}/\lambda_i), \quad (38) $$

where $a_0 = \delta\tilde{p}_0/2$ corresponds to $\omega_X = 0$, while $a_i = \mathfrak{R}(\delta\tilde{p}_i^*)/2^{N+1}$ and $b_i = -\mathfrak{F}(\delta\tilde{p}_i^*)/2^{N+1}$ are for $\omega_X = \omega_{X_i} = (2\pi/\lambda_i) = (2\pi i/\Delta\mathcal{X})$ with $i = 1, \dots, 2^{N+1}$, and $\mathfrak{R}$ and $\mathfrak{F}$ denote real and imaginary parts. To evaluate the surface seepage velocity, we then superimpose modal contributions as
the subtlety of $f_Y \propto 1/Y$ for our own data set, the range of $f_Y$ would roughly imply that the apparent $D_{eff}$ is 70 to 730 times weaker than $D$. Similar observations are found in moisture-absorbing pharmaceutical powders (Louge et al., 2021). Louge, Valance, Babah, et al. (2010) did not consider such hindered diffusion. However, because they derived the steady distribution of moisture through sand ripples, their model needs no repair, unless sand is exposed to rapid variations in surface moisture.

While temperature variations described in section 3 conformed to the linear partial differential equation (PDE) (8), the corresponding model of vapor mass fraction in Eq. (26) is highly non-linear, as it strongly couples $Y$ to subsurface variations in $T$. Another complication is a variable advection flow within the dune subsurface, which we consider next.

## 8 Seepage

It would not be straightforward to measure the small interstitial velocity $u$ in the field. However, the recorded mass fraction $Y$ (Fig. 3) exhibits so little noise that its temporal derivative can be reliably evaluated using a smoothing spline. Then, it is possible to fit the value of $u^*$ at each recorded time step by transforming the governing PDE (26) into an ordinary differential equation (ODE) where the unsteady term $U_M+U_T$ is calculated from data,

$$ U_M \equiv \left[ \frac{p^*}{T^*} + \mathbb{R}\Omega_1 f_Y(1-Y) \right] \frac{\partial Y}{\partial t^*} \approx \mathbb{R}\Omega_1 f_Y(1-Y) \frac{\partial Y}{\partial t^*}, \quad (29) \\ U_T \equiv \mathbb{R}\Omega_1 f_T(1-Y) \frac{\partial T^*}{\partial t^*}. $$

Rearranging Eq. (26) and replacing $\partial^2 T^*/\partial x^{*2}$ by $\partial T^*/\partial t^*$ using the heat equation, we find the ODE

$$ \frac{d^2 Y}{dx^{*2}} + \frac{dY}{dx^*} \left[ n_d \frac{\partial \ln T^*}{\partial x^*} - \frac{p^* u^*}{\mathbb{L} T^{*n_d+1}} \right] + a_T Y \left[ \frac{d \ln Y}{dx^*} \times \frac{\partial \ln T^*}{\partial x^*} + \frac{\partial \ln T^*}{\partial t^*} + (n_d - 1) \left( \frac{\partial \ln T^*}{\partial x^*} \right)^2 \right] = \frac{U_M + U_T}{\mathbb{L} T^{*n_d}}, \quad (30) $$

which we integrate subject to the measured $Y$ at depths $x_1^*$ and $x_{10}^*$ of sensors 1 and 10 using the two-point boundary value toolbox bvp4c of MATLAB. For consistency, we limit the integration to the interval $[x_1^*, x_{10}^*]$, since moisture on sensors 11 to 15 is sometimes too large for the probe to return a stable measurement of $Y$ (Fig. 3). Using the simplex search method of Lagarias et al. (1998), which is implemented as fminsearch in MATLAB, we repeat the integration of ODE (30) for several values of $u^*$ until the root-mean-square (rms) difference between signal and prediction is minimized.

Figure 6 illustrates this procedure. It is best viewed in the animation MovieCalculationDuneSeepageVelocity.mp4 available as supporting information. In this Fig., the left graph shows the profile of least rms error, as well as fits of $u^*$ to $\pm 7\%$ excursions in $Y$, which we use as uncertainty estimates in $u^*$. These excursions enclose spatial oscillations of relatively small amplitude in $Y$, which may be associated with the mild stratigraphy shown in Fig. 5. The corresponding uncertainties in $u^*$ (dashed lines in Fig. 7) are small. The integration of ODE (30) fails with $u^* = 0$ (dashed lines), thereby confirming that seepage advection is essential.

Shao et al. (2021) agreed with these observations. They suggested that a small advection on the order of 100 µm/s is necessary to capture data, and that thermal processes are crucial. However, by ignoring the soil's isotherm, Shao et al. (2021) had to adopt heuristic closures for the local evaporation rate, and they could not relate the apparent diffusion coefficient $D_{eff}$ to its known counterpart $D$ for water vapor in air. In addition, the
$$u_s = U^2 \tau_K \sum_{i=1}^{n} \omega_{x_i} \delta p_i^*. \quad (39)$$

To avoid spurious noise associated with survey errors at small wavelengths, we truncate the series at $n < 2^{N+1}$, such that remaining modes account for 97% of the total power $\sum a_i^2 + b_i^2$ of the DFT.

Because we did not conduct a detailed survey of the dune in March 2011, we derive $z$ from theodolite measurements carried out on the same dune in January 2017, and thus our calculations are only suggestive. However, they explain the seepage variations shown in Figure 7. As inset (d) of Figure 8 illustrates, variations of dune elevation along the wind direction, once converted to excursions in surface pressure, have their greatest Fourier power content at wavelengths on the order of the dune size ~100 m. Unlike sand ripples (Louge, Valance, Babah, et al., 2010), such wavelengths are too large to contribute appreciably to seepage. However, much power remains at smaller wavelengths ~10–30 m, particularly when wind flows toward the mean 159° historical bearing of this mobile dune. Because the gradient in Darcy's Equation 17 effectively amplifies surface pressure variations, these modes produce a seepage advection with significant longitudinal oscillations. As colored zones in Figure 8e suggest, parts of the dune inhale (red, $u_s > 0$), while others simultaneously exhale (blue, $u_s < 0$). These patterns change as wind changes direction, and as aeolian transport alters the dune's surface morphology.

The resulting predictions of $u_s/\tau_K$ in Figure 8 are consistent with magnitudes of $u_s$ calculated in Section 8. At the capacitance probe location (gray circle), inset a has a minimum of $u_s/\tau_K \approx -3,300 \text{ m/s}^2$ on the first day, while Figure 7 has $u^* \approx u_s^* \approx -100$. This corresponds to $\tau_K \approx 6 \ 10^{-8} \text{ s}$, which agrees well with $\tau_K \approx 5 \ 10^{-8} \text{ s}$ derived from parameters adopted in this article: $d \approx 351 \ \mu\text{m}$, $\nu \approx 0.545$, $K \approx 2.6 \ 10^{-10} \text{ m}^2$ (Equation 18), $\rho \approx 1.19 \text{ kg/m}^3$, $\mu \approx 1.8 \ 10^{-5} \text{ kg/ms}$, and $v/U \approx 0.036$.

However, a surprising observation from Figure 7 is the larger seepage advection velocity magnitude on 21 March ($t^* \sim 2.5$) than 20 March ($t^* \sim 1.5$) despite a smaller wind speed on the second day ($U \sim 5 \text{ m/s}$) than the first ($U \sim 7.5 \text{ m/s}$). A possible reason for this paradox may be a change in surface morphology from one day to the next, which affects $\tilde{z}$ and therefore shifts seepage patterns. Another could be associated with sand loading during aeolian transport on 21 March. Here, the energy required to lift sand grains on the first day may have lowered static pressure on the surface, thereby reducing the gradient between ambient and porous sand beneath it.

Finally, this calculation revealed that seepage advection varies greatly with location. Therefore, it is difficult to predict it for an arbitrary point using Equation 39, unless topography is known precisely for the entire dune at every instant, along with wind strength and direction. Fortunately, we showed in Section 8 that seepage can be inferred locally from measured spatiotemporal variations of vapor mass fraction.

## 10. Diurnal Predictions

In this section, we use MATLAB's pdepe to integrate PDE (26) with the seepage velocity inferred in Section 8 and predict the observed spatiotemporal evolution of vapor mass fraction $Y$ on relatively long diurnal time scales. Deferring a discussion of rapid variations until Section 12, we view the record of $Y$ during the onset of aeolian transport as two coinciding envelopes of a bistable signal. Integration is achieved in two overlapping periods. The first one begins at 19:12 Qatar time ($t^* \approx 0.80$) on 19 March and ends at 15:02 ($t^* \approx 1.63$) on 20 March when aeolian transport stops. Here, the PDE is subject to a surface BC $Y_1$ equal to the low continuous signal envelope at sensor 1 closest to the surface (Figure 3). The second period starts at the onset of aeolian transport on 20 March (12:29 or $t^* \approx 1.52$) and it persists until data acquisition is over on 21 March at 15:02 ($t^* \approx 2.63$). Its surface BC is the high envelope at sensor 1.

Unlike the procedure outlined in Section 8 to determine seepage velocity (Figure 6), we integrate the PDE over the entire probe depth, even if the probe does not return data at the deepest sensors. When it does, we note that the resulting values of $Y$ align on the thin dotted line marking the adsorbed mass fraction $\Omega \approx 1.65 \Omega_1$. This suggests that a suitable bottom BC is an invariant $\Omega$ at depth, $\partial\Omega/\partial x^*|_{\text{depth}} = 0$. In this case, using equations derived in Appendix C, such BC may be written
spatial resolution and accuracy of their time-domain reflectometry (TDR) was not sufficient to discern subtle variations in equilibrium vapor mass fraction that are needed to evaluate seepage time-history. They reported $5.4 \cdot 10^{-3} < \theta_w < 8.2 \cdot 10^{-3}$, where $\theta_w \simeq (\rho_p\nu/\rho_w)\Omega$ is water volume fraction. Combining Eqs. (3) and (7), and using properties of our sands as a rough estimate, their smallest detected $\theta_{w_{min}} \simeq 5.4 \cdot 10^{-3}$ translates to a vapor mass fraction $Y_{min} \simeq 0.02$, which is well over an order of magnitude larger than what our own capacitance probe discerned ($Y_{min} \simeq 0.0012$).

The right graph of Fig. 6 also shows relative magnitudes of the two contributions $U_M$ and $U_T$ to the total unsteady term. At long diurnal time scales for which grains achieve equilibrium with the surrounding $Y$, these terms are dominated by the rate $\partial\Omega/\partial t^*$ of the moisture exchange between solid and gas, ultimately captured by the functions $f_T$ and $f_Y$ derived in Appendix C. Because $U_M$ and $U_T$ often have opposite signs but similar magnitudes, temporal variations of mass fraction and temperature are both crucial, and they compete to govern vapor transfer through the dune subsurface.

A subtlety is that the seepage velocity $u^*$ varies along $x^*$ in response to density changes arising from temperature, pressure or the release of water vapor. Appendix F provides the corresponding ODEs (F1)–(F2), which we solve using MATLAB’s ODE45 4-th order Runge-Kutta algorithm, subject to dimensionless seepage velocity $u^* = u_s^*$ and pressure $p^* = p_s^*$ at the surface (subscript $s$), once $Y$ and its derivatives have been obtained from the solution of Eq. (30).

Inserts in Fig. 7 show the resulting profiles as solid lines, while dashed lines indicate calculations with $f_Y = f_T = 0$. Contrasting predictions with and without these functions underscores how the rate of water vapor release $\partial\Omega/\partial t^*$ affects $u^*$. Under planetary conditions where ambient pressure is small or gradients are large, ODEs (30) and (F1)–(F2) may be coupled strongly enough that the evolution of $Y$, $p$ and $u^*$ with $x^*$ would require their simultaneous solution. However, for this dune, depth variations in $u^*$ (solid lines) are small, and so it is legitimate to assume $u^* \simeq u_s^*$. A consequence is that seepage is not directly affected by internal gradients of temperature or mass fraction. Instead, as the next section articulates, it is likely due to wind-driven variations of surface pressure.

In short, integration of ODE (30) yields our best estimate of the dimensionless instantaneous velocity

$$v_s^* = +(1 - \nu)u_s^* \qquad (31)$$

counted > 0 into sand from the atmospheric boundary layer (ABL). By fitting the recorded depth profile of $Y$ with a smoothing spline and differentiating the result at the surface, we also evaluate the diffusive vapor flux $\dot{\Psi}_s'' = -\rho(1 - \nu)(D/\varpi)(\partial Y/\partial x)_s$ across the surface, counted > 0 into sand. Made dimensionless with $\rho_{st}\sqrt{\alpha/J}$, this flux

$$\Psi_s^* = -(1 - \nu)\mathbb{L}T_s^{*n_d}(\partial Y/\partial x^*)_s, \qquad (32)$$

augmented by the advection $u_s^*Y_s p^*/T_s^*$, contributes to the surface boundary condition for the ABL, which we will examine in section 13.

## 9 Origin of seepage

Louge, Valance, Babah, et al. (2010) recognized that wind-driven seepage advects moisture and dust through sand ripples by letting the Bernoulli effect induce a higher static pressure at their troughs than at their crests. Gao et al. (2018) staged a wavy surface to investigate its role in drying. By deploying a porous plastic ripple with harmonic amplitude $h_0$ and wavelength $\lambda_0$ in a wind tunnel, Musa et al. (2014) confirmed that the pore pressure gradient driving this advection scales as $\rho U^2 h_0 / \lambda_0^2$, where $U$ is the bulk wind speed. With origin above the trough at mid-depth between trough and crest, they
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Detail of Figure 3 showing two envelopes of the vapor mass fraction *Y* during a period of aeolian transport. Circles are data of sensors shown versus dimensionless time *t*\*. Thin lines indicate when the signal jumps from one envelope to the other. For each sensor, the lower *Y* corresponds to aeolian sand drying the surface. The superimposed lines are smoothing spline fits used to evaluate *∂Y/∂t*\* in the analysis. The solid and dashed lines marked "ambient" are *Y*<sub>a</sub> measured at the top and bottom hygrometers of the weather station, respectively. Inset: the resulting mean difference

$$ 
\Delta Y = (t_1^* - t_0^*)^{-1} \int_{t_0^*}^{t_1^*} (Y_u - Y_i) dt^* \text{ between the upper } Y_u \text{ and lower } Y_i 
$$

envelopes of *Y* versus dimensionless depth *x*\* (circles) with *t*<sub>0</sub><sup>*</sup> ≈ 1.52 and *t*<sub>1</sub><sup>*</sup> ≈ 1.63. The line is the model of Equation 51 with *ℓ*<sup>*</sup> = 0.28 ± 0.05.</figcaption></figure>

$$
\left. \frac{\partial Y}{\partial x^*} \right|_{\text{depth}} = Y \left[ 1 + Y (M^* - 1) \right] \left. \frac{\partial T^*}{\partial x^*} \right|_{\text{depth}} \left[ \frac{1}{T^*} \ln \left( \frac{1}{\text{RH}} \right) + \frac{T_{A_a} T_{\text{st}}}{(T - T_{A_0})^2} \right] . (40)
$$

We impose it as the bottom of our integration domain, set for simplicity at the deepest sensor. To show predicted variations with depth and time, and to compare measurements and predictions, we produced a movie available in the Supporting Information, with typical frames in Figure 4. This analysis also underscores the coupling of spatiotemporal variations of *T* and *Y*. Making *f*<sub>*T*</sub> ≡ 0 to eliminate it from PDE (26) would suppress the multiple observed inflections in the depth profile of *Y*, and altogether fail to capture data.

Finally, deeply buried water from antecedent rains, reported in Table 1 of Louge et al. (2013), matters to spatiotemporal variations of adsorbed moisture on sand grains nearer the surface. Its effects are captured by BC (40) at the bottom of our near-surface domain without liquid water. This BC can also be used to match our model of the near-subsurface to solutions for the deeper, moister domain below, for example, those of Kamai and Assouline (2018), Assouline and Kamai (2019), and Shao et al. (2021), or our own integration of Richards' equation (Louge et al., 2013; Richards, 1931).

However, because unsaturated porous media are governed by a history-dependent and hysteretic relation between capillary pressure and liquid volume fraction (Xu & Louge, 2015), definitive models of deeper liquid-infused regions are more challenging than our own. Conversely, their measurement techniques of water content are well established (Kidron & Kronenfeld, 2020; Kizito et al., 2008).

In short, until now, spatiotemporal variations of adsorbed water had not been measured in the liquid-free near-surface region of hyper-arid dunes with enough accuracy or spatial resolution. Doing so allowed us to show that this region could be analyzed without arbitrary closures. The principal remaining challenge is to model the deeper unsaturated region holding free liquid, where capillary effects continue to defy rigorous predictions.

11. Evaporation at the Grain Scale

The capacitance probe detects phenomena that other instruments with thicker measurement volume, such as microlysimeters (Boast & Robertson, 1982) or frequency domain reflectometry sensors (Kizito et al., 2008), cannot discern. In particular, during periods of aeolian sand transport, the capacitance probe returned values of *Y* belonging to two distinct envelopes switching rapidly from one to the other.

Figure 9 shows that the difference ΔY between these envelopes decreases exponentially with depth, that the lowest envelope coincides with the ambient vapor mass fraction Y_a, and that the probe hardly ever discerned any value of Y in between the envelopes. These observations suggest that strands of saltating particles meandering over the dune dried the surface intermittently, thereby sending rapid waves of Y downward on time scales too small for water vapor to come to equilibrium with the thin film of water adsorbed on sand grains.

To support this explanation, we briefly summarize the underlying physics of evaporation at the grain scale, and we show that the characteristic time τ for the film's return to an equilibrium is governed by slow first-order reaction kinetics, rather than by diffusion to the surface of individual grains. Then, we conduct a linear stability analysis of Equation 21 in the limit of short time scales. The calculations are consistent with the observed exponential decrease of the wave amplitude with depth. The corresponding length scale then lets us estimate τ, thus providing insight into the drying process for these sands.

During our two consecutive days of field data collection, the top 10 sensors of our capacitance probe detected no free liquid water within the first 15 cm below the surface. In this case, the exchange of moisture conformed to the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Snapshots of MovieCalculationDuneSeepageVelocity.mp4 available as supporting information. Left: best fits (lines) of recorded available data (symbols) of Y vs $x^*$ to the seepage velocity $u^*$ using ODE (30) at 08:03 on March 20. Dotted lines mark $Y \pm 7\% Y$, from which we estimate uncertainties in $u^*$. The integration of ODE (30) without seepage, shown as a dashed line, fails to capture data. Blue, yellow and red vectors are proportional to dimensionless fluxes of, respectively, advection $u_s^*Y_s p^*/T_s^*$, diffusion $\Psi_s^*/(1 - \nu)$ from Eq. (32), and their sum at the dune surface. Later in the movie (not shown here), the signal splits at 13:25 on March 20 in two bistable envelopes analyzed in section 12. During the subsequent period of high winds, dryer conditions prevail with aeolian transport (pink) than without (blue). Right: measured unsteady terms $U_M$ (blue line), $U_T$ (red line) and $U_M + U_T$ (black dashed line) defined in Eq. (29).</figcaption></figure>

showed that a ripple with ondulations (counted positive downward like $x$)

$$h = h_0 \cos(\omega_{X_0} \chi) \qquad (33)$$

along the wind direction of unit vector $\hat{X}$ produces a pore pressure field in the subsurface half-space that is a solution of the Laplace equation governing quasi-steady Darcy flow,

$$p = p_a + \rho v_{\tau}^2 \times \delta p^* \exp(-\omega_{X_0} x), \qquad (34)$$

where $\omega_{X_0} = 2\pi/\lambda_0$ is the single-mode wavenumber of the harmonic ripple along $\chi$, $v_{\tau}$ is the turbulent shear velocity, and

$$\delta p^* = 2\pi(h_0/\lambda_0)\sqrt{C^2 + D^2} \cos[(\omega_{X_0} \chi) - \phi_0] \qquad (35)$$

is the dimensionless excursion from the ambient pressure $p_a$ at the surface, with phase lag $\phi_0 = \arctan(D/C) > 0$. Fourrière (2009) and Fourrière et al. (2010) calculated the coefficients $C$ and $D$ in terms of the parameter $\ln(2\pi z_0/\lambda_0)$, where $z_0$ is the hydrodynamic roughness (Claudin et al., 2016). Through Darcy's Eq. (17), the gradient of the pressure field in Eq. (34) then induces the seepage velocity at $x = 0$,

$$u_s = U^2 \left[ \frac{K \rho (v_{\tau}/U)^2}{\mu(1-\nu)} \right] \omega_{X_0} \delta p^* \equiv U^2 \tau_K \omega_{X_0} \delta p^*, \qquad (36)$$

which, like surface pressure, features a peak lagging behind the harmonic ripple profile by a distance $\lambda_0\phi_0/(2\pi)$ along the wind. In Eq. (36), $\tau_K \equiv K\rho(v_{\tau}/U)^2/[\mu(1-\nu)]$ has units of time and incorporates the ratio $v_{\tau}/U$, which typically varies on the sand surface.
“hygroscopic” regime, in which capillary forces are insignificant, and water molecules either form a thin layer on the surface of grains, or exist as a vapor in air.

The amount of moisture bound to grains in the size class of diameter *d* is measured by the mass fraction *Ω*<sub>d</sub> of water that the class holds. In the vapor phase within the interstitial space between packed grains, the mass fraction of water is *Y*. Given enough time, *Ω*<sub>d</sub> and *Y* reach an equilibrium, denoted by the subscript *e*, in which the rate of evaporation from grains is balanced by condensation at a given temperature. Appendix B provides the form of the corresponding “isotherm” of size class *d*. Once integrated over the particle-size distribution, Equation B1 yields the isotherm for the average mass fraction *Ω*<sub>e</sub> that Louge et al. (2013) recorded in terms of the relative humidity RH<sub>e</sub> in a humidity chamber.

In general, Ω is a macroscopic measure of the thin coating of water arranged as molecular layers on grain surfaces. As Equations 1 and 2 indicate (Shahraeeni & Or, 2010), his coating has mean thickness

$$ \ell_w \approx \left(\frac{d_{2,3}}{6}\right) \left(\frac{\rho_p}{\rho_w}\right) \frac{\Omega_1}{T^{*1/3} \ln^{1/3}[p_{\text{sat}}/(pM^*Y_e)]} \quad (41) $$

where ρ<sub>p</sub> ≈ 2,630 and ρ<sub>w</sub> ≈ 997 kg/m<sup>3</sup> are, respectively, sand and water material densities, d<sub>2,3</sub> ≈ 300 µm is the Sauter mean diameter that arises from the integration over the particle size distribution, and Ω<sub>1</sub> ≈ 0.0013 (Louge et al., 2013). In this field campaign, the probe recorded 0.0011 < Y<sub>e</sub> < 0.017, or 120 nm < ℓ<sub>w</sub> < 330 nm, corresponding to 400 to 1,100 layers of water molecules of individual 0.3 nm thickness (Yeşilbaş & Boily, 2016).

With so many layers, it is reasonable to neglect interactions with the solid as water molecules escape from grains or adsorb on them. In this case, grain evaporation can be regarded as a non-equilibrium process releasing water molecules at the interface between interstitial air and the thin water coating, mitigated by a reverse reaction of vapor condensation (Koffman et al., 1984). Even under controlled laboratory conditions, the coupling of vapor and liquid phases, or possibly local variations of temperature near the gas-liquid interface (Jafari et al., 2018), make it challenging to model this process (Persad & Ward, 2016). Because we must also contend with grains of varying size, shape, surface geometry, impurities and composition, immersed in a random interstitial space, we resort to a rudimentary model that is inspired from the Hertz-Knudsen-Schrage (HKS) theory (Marek & Straub, 2001). In this view, the net vapor mass flux at the saturated gas-liquid interface,

$$ \dot{W}_{\text{HKS}}'' = \frac{2\kappa_e}{2 - \kappa_c} \sqrt{\frac{\text{MW}_{\text{H}_2\text{O}}}{2\pi\hat{R}}} \left[ \left( \frac{\kappa_c}{\kappa_e} \right) \frac{p_{\text{sat}}}{T_{\ell}^{1/2}} - \frac{p_{\text{H}_2\text{O}}}{T_v^{1/2}} \right] \quad (42) $$

is built upon the prediction of the kinetic theory of gases for the water mass flux p<sub>H₂O</sub>√MW<sub>H₂O</sub>/(2πRT) crossing a control surface under a Maxwellian velocity distribution function (Vincenti & Kruger, 1965). In this expression, p<sub>H₂O</sub> ≈ pYM* is the partial pressure of water vapor in the interstitial space, R ≈ 8.314 J/m K is the fundamental gas constant, and the “accommodation coefficients” of evaporation κ<sub>e</sub> and condensation κ<sub>c</sub> represent the probability that a water molecule colliding with the interface undergoes its respective first-order phase transition (Y. Q. Li et al., 2001).

For the tight porous medium of the dune with grains of relatively high heat capacity, it is reasonable to assume that the temperature T<sub>ℓ</sub> at the gas-liquid interface and its interstitial vapor counterpart T<sub>v</sub> are equal to the local value *T* of the “heat bath” that our instrument measured. Then, because the vapor is saturated at the interface, RH ≡ 1, an expression similar to Equation 7 sets the mass fraction there, *Y*<sub>sat</sub> ≈ *p*<sub>sat</sub>/*(pM*)*.

Meanwhile, at equilibrium, evaporation and condensation balance (Ẇ''<sub>HKS</sub> = 0), and molecular diffusion makes the interstitial vapor fraction *Y* uniform and equal to the value *Y*<sub>e</sub> recorded by the isotherm of Equation B1 for each size class. Therefore, the expression (κ<sub>c</sub>/κ<sub>e</sub>)Y<sub>sat</sub> is identified with *Y*<sub>e</sub>, and Equation 42 becomes

$$ \dot{W}_{\text{HKS}}'' \approx \frac{2\kappa_e}{2 - \kappa_c} \sqrt{\frac{\text{MW}_{\text{H}_2\text{O}}}{2\pi\hat{R}T}} pM^* (Y_e - Y) \equiv k_m (Y_e - Y) \quad (43) $$

which serves as definition of the rate constant *k*<sub>m</sub> of a balanced first-order reaction of evaporation returning the system to equilibrium.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time-history of dimensionless seepage velocity $u^*$ at the dune surface vs $t^*$ (solid line) with surrounding uncertainties (dashed lines) associated with the measured $Y \pm 7\% Y$ shown in Fig. 6. The grey regions enclose values of $u^*$ calculated in Appendix F for depths $x^* < 0.87$ above sensor 10. Splitting of the $u^*$-record in the period $1.4 < t^* < 1.63$ arises as signal envelopes separate during aeolian transport. Solid lines in insets A, B and C are typical profiles of $u^* - u_s^*$ calculated through these depths at times $t^* = 0.98, 1.54$ and $2.47$; dashed lines illustrate the role of vapor release on $u^*$ by artificially making $f_T = f_Y = 0$ in these calculations. Insets and grey regions reveal that $u^*$ varies little with depth in this case.</figcaption></figure>
Because the interstitial space may also possess local gradients of Y, vapor diffusion on the pore scale could also play a role in setting the flux at the film's interface. For example, Philip (1964) and Shahraeeni and Or (2010) assumed that the transport of vapor to grain asperities is dominated by diffusion. (Section 12 will show that such assumption is not consistent with the depth scale of subsurface waves that we observed).

In general, as in the combustion of a solid particle (Kanury, 1975), there are two limiting regimes. In the “kinetic limit,” $k_m$ is small relative to $\rho D/d$, so diffusion of coefficient D erases gradients of Y in a pore with interstitial air of local density $\rho$. Conversely, in the “diffusion limit,” the evaporation rate of Equation 43 is much more intense than the local diffusive flux.

As Appendix G shows, the dominant process setting the net vapor flux $\dot{W}_s''$ at the grain surface is the slowest of the two. Once it is identified, one can infer the rate of change of the mass fraction $\Omega_d$ for each size class. Assuming for simplicity that grains of material density $\rho_p$ are spherical,

$$ \rho_p \nu \frac{\partial \Omega_d}{\partial t} = - \frac{6\nu \dot{W}_s''}{d} \equiv \frac{\rho(1-\nu)}{\tau_d} (Y - Y_e), \quad (44) $$

which defines the characteristic reaction time $\tau_d$ for the size class of diameter d. Upon integration over the particle size distribution in Appendix B, Equation 44 yields the ODE that governs the rate of change of the overall moisture mass fraction $\Omega$ adsorbed on bulk sand of solid volume fraction $\nu$,

$$ \frac{\partial \Omega}{\partial t} = \frac{\rho(1-\nu)}{\rho_p \nu} \left( \frac{Y - Y_e}{\tau} \right). \quad (45) $$

In the next section, we conduct a stability analysis to determine the origin of the subsurface waves of Y revealed by our instrument, and to calculate the characteristic depth of their decay. In turn, this depth sets bounds for possible values of $\tau$ and, as Appendix G shows, conclusively places evaporation at the grain scale in the kinetic limit.

## 12. Evanescent Waves

In Section 7, we considered the limit when variables change slowly enough that the water film of mass fraction $\Omega$ remains in equilibrium with the interstitial Y. When instead Y varies on a time scale shorter than $\tau$, $\Omega$ does not change appreciably, so the right side of Equation 21 becomes small. This situation arises during aeolian transport (Dupont et al., 2013), when clustered strands of dry particles meander over the dune, alternately forcing the surface Y to adopt the ambient value $Y_a$, or letting it quickly readjust to its higher original level without them. As Figure 9 shows, recorded mass fractions then switch between these two extreme states, and the resulting difference $\Delta Y$ decreases exponentially with depth. Toward explaining this, we write Equation 21 in a dimensionless form denoted by asterisks,

$$ \frac{p^*}{T^*} \frac{\partial Y}{\partial t^*} + \frac{p^*}{T^*} u^* \frac{\partial Y}{\partial x^*} - \mathbb{L} \frac{\partial}{\partial x^*} \left[ T^{*nd} \left( \frac{\partial Y}{\partial x^*} + a_T Y \frac{\partial \ln T^*}{\partial x^*} \right) \right] = - \frac{p^*}{T^*} \frac{(Y - Y_e)}{\tau^*} (1 - Y), \quad (46) $$

where $\tau^* \equiv \tau/J$. Here, the diffusion term involves a Soret effect with coefficient $a_T$ coupling variations of T and Y, which Appendix E dismisses as unimportant.

In this equation, $Y_e$ is known in terms of the instantaneous $\Omega$ held by grains from Equations 3 and 7. On short time scales, it is therefore invariant as equilibrium cannot be restored fast enough. In this analysis, we consider waves of Y traveling about its mean value $\tilde{Y}$ at the dimensionless frequency $f^* \equiv fJ$,

$$ Y = \tilde{Y} + \delta Y \exp \left[ t \left( 2\pi f^* t^* - k^* x^* \right) \right], \quad (47) $$

where $t^2 = -1$. Substituting this expression in Equation 46, expanding to leading order in $\delta Y/\tilde{Y} \ll 1$, simplifying $(1 - Y) \simeq 1$, subtracting the time-averaged equation, applying derivatives of Equation 47, and dividing by the exponential term and $\delta Y$, Equation 46 yields the characteristic equation

$$ A_c k^{*2} - \iota B_c k^* + C_c = 0, \quad (48) $$
For our dune, we estimate $λ_0 \simeq 220$ m as approximately four times its toe to crest distance, and $h_0 \simeq 4.5$ m as its crest height (Fig. 8). Anemometry carried out near the capacitance probe location in January 2017 also found $z_0 \simeq 16$ µm. With $\ln(2\pi z_0/λ_0) \simeq -14.6$, Fourrière (2009) calculated $C \simeq -1150$ and $D \simeq 49$. Because $C$ and $D$ has a weak logarithmic dependence on $z_0$ and modal wavelengths, we ignore variations of these quantities on the dune and, for simplicity, treat $C$ and $D$ as constants.

Because the single-mode wavenumber $ω_{X_0} = 2π/λ_0$ is very small, Eq. (36) would predict a much smaller seepage than a ripple of comparable aspect ratio, if the dune possessed a harmonic profile. However, the dune shape is not harmonic, but instead features gentle undulations in its surface elevation $z(\mathcal{X}, \mathcal{Y})$ along the wind unit direction $\hat{\mathbf{X}}$ and its cross-wind counterpart $\hat{\mathbf{Y}}$. The small ratio of dune height on wavelength $h_0/λ_0 \sim 0.02$ further suggests that $z(\mathcal{X}, \mathcal{Y})$ may be expanded in Fourier modes that contribute linearly to the overall pressure excursion (Musa et al., 2014). Then, as the following analysis shows, some of these modes may possess large enough wavenumbers to induce significant seepage.

Inspired by the theory of P. S. Jackson and Hunt (1975) and Hunt et al. (1988), Kroy et al. (2002) derived the transfer function relating the two-dimensional Fourier transforms of dune elevation $\tilde{z}$ and perturbation in turbulent shear stress $\delta\tau$ in terms of two parameters $A$ and $B$. The dimensionless pressure excursion $\delta p^*$ is subject to a similar transfer function, in which $A$ and $B$ are respectively replaced by $C$ and $D$,

$$ \frac{\delta \tilde{p}^*}{\tilde{z}} = \frac{\omega_X (C\omega_X + \iota D|\omega_X|)}{\sqrt{\omega_X^2 + \omega_Y^2}}, \qquad (37) $$

where $ω_X \equiv 2π/λ_X$ and $ω_Y \equiv 2π/λ_Y$ are wavenumbers along $\hat{\mathbf{X}}$ and $\hat{\mathbf{Y}}$, respectively. With Eq. (37), the two-dimensional Fourier transform of $z$ obtained from a dune survey therefore yields $\delta p^*$, which is then inverted to recover the pressure excursion everywhere on the dune surface. To evaluate $\tilde{z}$ at each time step, we tilt the dune slightly so it rises above an artificially flattened desert floor, rotate it around the origin of $(\mathcal{X}, \mathcal{Y})$ at the barycenter of its outline to align $\hat{\mathbf{X}}$ with the instantaneous wind direction, interpolate survey points using the C1-continuous scattered Interpolant of MATLAB, subdivide the square central domain circumscribing the rotated dune outline in $2^N$ cells with $N = 9$, and buttress this domain with eight identical squares holding zero elevation to avoid aliasing. The resulting physical domain of width $\Delta\mathcal{X}$ therefore holds $2^{N+2}$ cells along both directions, leading to $2^{N+2}$ discrete wavenumbers $ω_X = (2π/\Delta\mathcal{X}) \times (0, \dots, 2^{N+1}, -2^{N+1} + 1, \dots, -1)$, organized in the peculiar order with which MATLAB optimizes its discrete Fourier transforms (DFT). Comparing predictions of Eq. (37) from the two-dimensional fft2 DFT and its one-dimensional counterpart fft with $ω_Y = 0$, we find that transverse wavenumbers $ω_Y$ play a minor role. Therefore, for simplicity, we reconstruct the pressure excursion $\delta p^*$ on any wind-directed transect using the first $1+2^{N+1}$ coefficients of the one-dimensional DFT $\delta \tilde{p}^*$ in Eq. (37) with $ω_Y = 0$,

$$ \delta p^* = a_0 + \sum_{i=1}^{2^{N+1}} \delta p_i^* = a_0 + \sum_{i=1}^{2^{N+1}} a_i \cos(2\pi\mathcal{X}/\lambda_i) + b_i \sin(2\pi\mathcal{X}/\lambda_i), \qquad (38) $$

where $a_0 = \delta p_0/2$ corresponds to $ω_X = 0$, while $a_i = \Re(\delta \tilde{p}_i^*)/2^{N+1}$ and $b_i = -\Im(\delta \tilde{p}_i^*)/2^{N+1}$ are for $ω_X = ω_{X_i} = (2π/λ_i) = (2πi/\Delta\mathcal{X})$ with $i = 1, \dots, 2^{N+1}$, and $\Re$ and $\Im$ denote real and imaginary parts. To evaluate the surface seepage velocity, we then superimpose modal contributions as

$$ u_s = U^2 \tau_K \sum_{i=1}^{n} \omega_{X_i} \delta p_i^*. \qquad (39) $$

To avoid spurious noise associated with survey errors at small wavelengths, we truncate the series at $n < 2^{N+1}$, such that remaining modes account for 97% of the total power $\sum a_i^2 + b_i^2$ of the DFT.
where

$$
\begin{align}
A_c &\equiv \mathbb{L} T^{*n_d} \nonumber \\
B_c &\equiv \frac{p^*}{T^*} u^* - \mathbb{L} (n_d + a_T) T^{*n_d} \frac{\partial \ln T^*}{\partial x^*} \tag{49} \\
C_c &\equiv -a_T \mathbb{L} T^{*n_d} \left[ (n_d - 1) \left( \frac{\partial \ln T^*}{\partial x^*} \right)^2 + \frac{\partial \ln T^*}{\partial t^*} \right] + \frac{p^*}{T^*} \left[ 12\pi f^* + 1/\tau^* \right], \nonumber
\end{align}
$$

and $T^*$ and $u^*$ are sampled in the midst of the aeolian transport period at $t^* \simeq 1.57$. There are two solutions $k^{*\pm} = k_r^{*\pm} - t k_l^{*\pm}$ to quadratic Equation 48, with real and imaginary parts $k_r^{*\pm}$ and $k_l^{*\pm}$,

$$k^{*\pm} = t \frac{B_c}{2A_c} \pm \frac{1}{2A_c} \sqrt{-B_c^2 - 4A_c C_c} \quad (50)$$

The two resulting traveling waves in $Y$ have the form $\delta Y \exp(k_l^{*\pm}x^*) \times \exp[t(2\pi f^*t^* - k_r^{*\pm}x^*)]$. Of these, only the $k^{*\pm}$ wave has $k_l^{*\pm} < 0$, and therefore vanishes as $x^* \to \infty$. In the regime of rapid oscillations and slow return to equilibrium (Figure 9), we match its decreasing amplitude along $x^*$

$$\Delta Y = \delta Y \exp(-|k_l^{*\pm}|x^*) \quad (51)$$

to the exponential decay with characteristic length $\ell^* \equiv -1/k_l^{*\pm} = 0.28 \pm 0.05$, or $\ell = 4.9 \pm 0.9$ cm. This relatively short $\ell$ justifies why it is crucial to measure water mass fraction with a probe of narrow vertical spatial resolution.

Meanwhile, Equations 49 and 50 yield $k_l^{*\pm}$ in terms of the known parameters $\mathbb{L}$, $n_d$, and $a_T$, the measured variables $T^*$, $u^*$, and $p^*$, as well as $\tau^*$ and $f^*$. Therefore, they can be recast as the relation between $\tau^*$ and $f^*$ that yields the measured $\ell^*$. For our conditions, this relation is conveniently fitted as $\tau^* \simeq \tau_l^* + \delta\tau^* / [(1/f^* - (1/f_l^*)^*)]$, where $(1/f_l^*) \simeq 0.0084 (12\text{min})$, $\delta\tau^* \simeq 4.2 \times 10^{-5}$ and $\tau_l^* \simeq 0.01 (14\text{min})$. Although we did not sample $Y$ fast enough to establish its high frequency content, an inspection of oscillations in the bistable $Y$-record (thin lines in Figure 9) suggests that $0.013 \lesssim 1/f^* \lesssim 0.033 (19-48\text{min})$, thus yielding $\tau^* = 0.015 \pm 0.004 (22 \pm 5\text{min})$.

Such time is much larger than the diffusion time $\tau_{\text{DIFF}} \sim 240 \ \mu\text{s}$ estimated in Appendix G. Therefore, if the diffusion limit prevailed, sand grains would have ample time to return everywhere to an equilibrium with the surrounding vapor mass fraction. In that case, as Section 7 showed, the right side of Equation 46 would become $-\mathbb{R}\Omega_1 f_Y \partial Y / \partial t^*$, where temperature changes negligibly, $f_Y \propto 1/Y$, and $\mathbb{R} \equiv \rho_p \nu / [\rho_{\text{st}}(1-\nu)]$. Then, the procedure leading to characteristic Equation 48 would leave $A_c$ and $B_c$ unchanged, while replacing $C_c$ in Equation 49 with

$$C_c' = \frac{p^*}{T^*} t 12\pi f^* + \mathbb{R}\Omega_1 \bar{f}_Y \left( t 12\pi f^* - \frac{\partial \ln \bar{Y}}{\partial t^*} \right), \quad (52)$$

an expression that no longer involves $\tau$. Using recorded $T^*$, $u^*$, $\bar{f}_Y$, and $\partial \ln \bar{Y} / \partial t^*$ at $t^* \simeq 1.6$, we find that $\ell^*$ would never exceed 0.045, which is much smaller than the observed $\ell^*$. Therefore, we conclude that evaporation is a kinetic-limited process in our case.

A common alternative approach to the kinetics of evaporation assumes that water in the adsorbed film is at thermodynamic equilibrium with the surrounding RH, while time-history is reintroduced through mass conservation and diffusion (Equation 21). Such equilibrium approach to the evaporation rate, which is suitable for the drying of porous media such as wood or biomass, see for example, Nasrallah and Perré (1988), Zhang and Datta (2004), and Borujerdi et al. (2019), would fail for these hyper-arid sands.

Lastly, our determination of $\tau^* \simeq \tau_{\text{HKS}}^* = 0.015 \pm 0.004$ yields the accommodation coefficient of evaporation. Substituting the mean Sauter diameter $d_{2,3}$ for $d$ in Equation G1, and taking $\kappa_c$ to be on the same order, we find $\kappa_c = (2.3 \pm 0.5) \times 10^{-10}$. Such value is very small by the standards of pure, smooth surfaces. However, microscopic impurities, which are certain to contaminate desert sands, and interface coverage (Rubel & Gentry, 1984) are known to thwart interfacial mass transfer and reduce $\kappa_c$ by several orders of magnitude (Marek & Straub, 2001).

Coincidentally, this low probability is also consistent with an activation energy $\hat{E} = -\hat{R}T \ln \kappa_c = 54.6 \pm 0.6 \text{ kJ/mole}$, another common interpretation of the evaporation process (Cappa et al., 2007; Rubel & Gentry, 1984; Smith et al., 2006; Tsukahara et al., 2009; Y. Q. Li et al., 2001). This value is typical of the energy needed to release
Because we did not conduct a detailed survey of the dune in March 2011, we derive $z$ from theodolite measurements carried out on the same dune in January 2017, and thus our calculations are only suggestive. However, they explain the seepage variations shown in Fig. 7. As inset (D) of Fig. 8 illustrates, variations of dune elevation along the wind direction, once converted to excursions in surface pressure, have their greatest Fourier power content at wavelengths on the order of the dune size ~100 m. Unlike sand ripples (Louge, Valance, Babah, et al., 2010), such wavelengths are too large to contribute appreciably to seepage. However, much power remains at smaller wavelengths ~10 to 30 m, particularly when wind flows toward the mean 159° historical bearing of this mobile dune. Because the gradient in Darcy's Eq. (17) effectively amplifies surface pressure variations, these modes produce a seepage advection with significant longitudinal oscillations. As colored zones in Fig. 8E suggest, parts of the dune inhale (red, $u_s > 0$), while others simultaneously exhale (blue, $u_s < 0$). These patterns change as wind changes direction, and as aeolian transport alters the dune's surface morphology.

The resulting predictions of $u_s/\tau_K$ in Fig. 8 are consistent with magnitudes of $u_s$ calculated in section 8. At the capacitance probe location (grey circle), inset A has a minimum of $u_s/\tau_K \simeq -3300 \text{ m/s}^2$ on the first day, while Fig. 7 has $u^* \simeq u_s^* \simeq -100$. This corresponds to $\tau_K \simeq 6 \times 10^{-8} \text{ s}$, which agrees well with $\tau_K \simeq 5 \times 10^{-8} \text{ s}$ derived from parameters adopted in this article: $d \simeq 351 \text{ µm}$, $\nu \simeq 0.545$, $K \simeq 2.6 \times 10^{-10} \text{ m}^2$ (Eq. 18), $\rho \simeq 1.19 \text{ kg/m}^3$, $\mu \simeq 1.8 \times 10^{-5} \text{ kg/m.s}$ and $v_\tau/U \simeq 0.036$.

However, a surprising observation from Fig. 7 is the larger seepage advection velocity magnitude on March 21 ($t^* \sim 2.5$) than March 20 ($t^* \sim 1.5$) despite a smaller wind speed on the second day ($U \sim 5 \text{ m/s}$) than the first ($U \sim 7.5 \text{ m/s}$). A possible reason for this paradox may be a change in surface morphology from one day to the next, which affects $\tilde{z}$ and therefore shifts seepage patterns. Another could be associated with sand loading during aeolian transport on March 21. Here, the energy required to lift sand grains on the first day may have lowered static pressure on the surface, thereby reducing the gradient between ambient and porous sand beneath it.

Finally, this calculation revealed that seepage advection varies greatly with location. Therefore, it is difficult to predict it for an arbitrary point using Eq. (39), unless topography is known precisely for the entire dune at every instant, along with wind strength and direction. Fortunately, we showed in section 8 that seepage can be inferred locally from measured spatio-temporal variations of vapor mass fraction.

## 10 Diurnal predictions

In this section, we use MATLAB's `pdepe` to integrate PDE (26) with the seepage velocity inferred in section 8 and predict the observed spatio-temporal evolution of vapor mass fraction $Y$ on relatively long diurnal time scales. Deferring a discussion of rapid variations until section 12, we view the record of $Y$ during the onset of aeolian transport as two coinciding envelopes of a bistable signal. Integration is achieved in two overlapping periods. The first one begins at 19:12 Qatar time ($t^* \simeq 0.80$) on March 19 and ends at 15:02 ($t^* \simeq 1.63$) on March 20 when aeolian transport stops. Here, the PDE is subject to a surface boundary condition $Y_1$ equal to the low continuous signal envelope at sensor 1 closest to the surface (Fig. 3). The second period starts at the onset of aeolian transport on March 20 (12:29 or $t^* \simeq 1.52$) and it persists until data acquisition is over on March 21 at 15:02 ($t^* \simeq 2.63$). Its surface boundary condition is the high envelope at sensor 1.

Unlike the procedure outlined in section 8 to determine seepage velocity (Fig. 6), we integrate the PDE over the entire probe depth, even if the probe does not return data at the deepest sensors. When it does, we note that the resulting values of $Y$ align on the thin dotted line marking the adsorbed mass fraction $\Omega \simeq 1.65 \Omega_1$. This suggests that
solid-bound water in solid fuels (Borujerdi et al., 2019) or powders (Prado & Vyazovkin, 2011) and, as expected, it is above, —but on the order of—, the latent heat ≈43.9 kJ/m to vaporize unbound water.

In short, the characteristic decay depth of the evanescent waves suggests that evaporation from these sand grains is a kinetic-limited process. Although the latter may be interpreted as having a very low accommodation coefficient, it is more intuitive to regard it as an activated reaction.

## 13. Atmospheric Boundary Layer

The instantaneous evaporation flux exchanged with the ABL from the vadose zone is a crucial BC for understanding the hydrology (Assouline & Kamai, 2019; Shao et al., 2021) and microbiology (Kidron & Starinsky, 2019) of hyper-arid regions. Because it is important to distinguish evaporation from transpiration when vegetation is involved (R. G. Anderson et al., 2017), flux techniques deployed in the ABL, such as eddy-covariance, are supplemented by carbon dioxide measurements (Scanlon & Kustas, 2010), sometimes involving stable isotopes (Griffis, 2013).

Measurements of the evaporation flux from within the sand surface are equally challenging. For example, in a fog desert, Feigenwinter et al. (2020) used an automated microlysimeter (ML) that recorded the mass of non-rainfall water (NRW) deposited in a sampling cup mounted flush with the surface. However, Kidron and Kronenfeld (2020) noted how the ML technique can overestimate the NRW amount, especially at night.

Meanwhile, transport in the ABL faces instabilities arising from natural convection during the day (Wyngaard, 2010), and it is complicated at night by stratified, weak turbulence (Mahrt, 2014; Optis et al., 2014) and by low-level jets intermittently perturbing the stable nocturnal boundary layer (Banta et al., 2007; Klein et al., 2016). For flat expanses of land, the resulting transports of momentum and heat are typically modeled using the Monin-Obukhov similarity (Monin, 1970). However, in the case of isolated sand dunes surrounded by a hard, more impermeable desert floor, it is unclear how such similarity can be extended to mass transfer of water vapor, or whether a boundary layer of vapor mass fraction developing from the dune leading edge may dominate the exchange. Therefore, it is opportune to validate mass transfer BCs at the base of the ABL at a realistic field scale with direct flux measurements from within the dune.

In this context, the capacitance probe described in Section 2 provides an alternative to traditional evaporation instruments by capturing spatiotemporal variations of the water vapor mass fraction Y just below the sand surface. As Section 8 showed, an integration of the governing equations then returns a best estimate of the dimensionless instantaneous wind-driven seepage velocity $u^*$. At the surface, denoted by the subscript s, this yields the superficial gas velocity

$$v_s^* = +(1 - v)u_s^* \qquad (53)$$

counted >0 into sand from the ABL.

In addition, by virtue of their accuracy and low noise, the recorded depth profiles of Y can be smoothed with a spline and differentiated to evaluate the vapor flux $\Psi_s'' = -\rho(1 - v)(D/\omega)(\partial Y/\partial x)_s$ across the surface, again counted >0 into sand. Made dimensionless with $\rho_{st}\sqrt{\alpha/J}$, the diffusive flux in Equation 32 augmented by its advection counterpart $u_s^*Y_s p^*/T_s^*$, constitutes a unique measurement of the mass transfer BC at the base of ABL. As Appendix H shows from the integral equations of the boundary layer, the proper way is to count it >0 into the ABL aloft as

$$\Sigma^* = -\Psi_s^* - \frac{p^*}{T_s^*}v_s^*(Y_s - Y_a). \qquad (54)$$

In this determination, the principal uncertainty is associated with the tortuosity $\omega \sim \pi/2$ (Shen & Chen, 2007), which we estimated for the nearly spherical and monodisperse grains of this dune.

As the panels of Figure 10 show, the diffusion flux into the ABL and its sign are, as expected, nearly proportional to the “driving potential” ($Y_s - Y_a$) in the second day, that is, to the difference between the surface mass fraction $Y_s$ and its ambient counterpart $Y_a$. However, two surprising observations stand out. First, the driving potential nearly vanishes during aeolian transport in the first day (Figure 10a), while the flux maintains roughly the
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Seepage predictions of Eq. (39) calculated for the dune shape in January 2017, the wind record *U* in Fig. 3, and the instantaneous wind direction $\hat{X}$. (A) Relative surface seepage velocity $u_s/\tau_K$ (m/s$^2$ > 0 *into* sand) at the position shown by the grey dot in inset (E) vs diurnal time *t*/*J*. (B) Instantaneous wind bearing $\theta$ (deg) vs *t*/*J*, and mean historical bearing from Louge et al. (2013) (horizontal line). In insets (A) and (B), the line is a visual fit, and vertical solid and dashed lines mark sunrise and sunset, respectively. (C) Solid line: predicted relative pressure excursion $\delta p^*$ vs distance *X* (m) along the transect marked by the dashed line in inset (E); dashed line: dune elevation *z* (m) along this transect. (D) Dimensionless spectral power $a_i^2 + b_i^2$ of pressure excursions in Eq. (38) vs wavelength $\lambda_i = \Delta X/i$ with $i = 1, \dots, 2^{N+1}$ for the wind bearings shown. (E) Snapshot of $u_s/\tau_K$ on the dune surface at 11:14 on March 21, marked by grey dots in A and B. The whole movie MovieDuneSurfaceSeepagePrediction.mp4 is available as supporting information. Black isolines are placed every 0.5 m. For best contrast, colors are obtained by converting $[1 + \tanh(u_s/(\tau_K\gamma_0))]/2$ to the jet color map of MATLAB with $\gamma_0 = 10^3$ m/s$^2$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Time-histories in the atmospheric boundary layer (ABL). (a) Driving potential ($Y_s - Y_a$), nearly vanishing during aeolian transport from $t^* \simeq 1.4$ to 1.6, thus imposing a Dirichlet boundary condition on the ABL. (b) Thick line: measured dimensionless net flux of Equation 54. Thin line: turbulent boundary layer (TBL) model of Equation 57 with $a_m \simeq 0.0010$ and $a_m \simeq 0.0015$ under stable ($L > 0$) and unstable ($L < 0$) conditions, respectively, $\Sigma^*_{\infty} \simeq 0.43$ in the first night, $\Sigma^*_{\infty} \simeq 0.22$ in the second night, and $\Sigma^*_{\infty} \simeq 0.075$ in the second day. Dotted line: TBL model of the second day extrapolated to the period of stability after sunrise. (c) Inverse Monin-Obukhov length, dimensionless with $\sqrt{\alpha J}$; dotted lines at $\pm$ the inverse dimensionless altitude $z_T$ of the top RH sensor. (d) Re$_{\chi}$ and Gr$_L^{1/2}$. The horizontal dashed line marks Re$_{x,crit}$ (Kays & Crawford, 1980).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Small difference $\Delta Y_a = Y_a(z_T) - Y_a(z_B)$ in the ambient water vapor mass fraction time-history between the top and bottom sensor elevations at $z_T \simeq 2$ m and $z_B \simeq 0.9$ m measured by the weather station above hard ground with much higher aerodynamic roughness $z_0 \simeq 0.011$ m (Louge et al., 2013) than found over the sand surface. Vertical lines mark sunrise and sunset. The horizontal dashed line is the mean $\Delta Y_a$. (b) Corresponding change in ambient temperature $\Delta T = T(z_T) - T(z_B) \equiv T_T - T_B$ (°K). (c) Ratio $U(z_T)/U(z_B)$. The horizontal dashed and dotted lines show $\ln(z_T/z_0)/\ln(z_B/z_0)$ and $z_T/z_B$, respectively consistent with a turbulent log-law and a viscous sublayer.</figcaption></figure>

same magnitude (Figure 10b). As discussed in Section 12, we attribute this phenomenon to enhanced mass transfer by saltating particles, which effectively produce a Dirichlet BC imposing $Y_s \simeq Y_a$ at the surface.

The second paradox is more puzzling. When the ABL is stable at night, the dry surface keeps losing moisture to the wetter ambient, in an apparent inversion of the proportionality between flux and driving potential. In addition, the net surface flux changes sign in the last quarter of the night, while the driving potential does not.

Profile measurements at our short weather station are not sufficient to resolve this paradox. However, data support the following comments. As Figure 3c showed, the ABL is stable at night, when a colder dune at $\sim 290^\circ$K gives up heat at long-wavelength infrared $\sim 10$ µm, and unstable during the day, when warm air rises from sand stricken by a solar irradiation that exceeds wind-driven losses. Therefore, turbulent momentum in the ABL is coupled with natural convection during the day (Wyngaard, 2010). To capture this effect, Monin and Obukhov proposed a similarity between the transports of momentum, heat and mass in the turbulent ABL (Monin, 1970), which we summarize in Appendix D. The similarity is captured by corrections to the Prandtl-von Kármàn log-law of the wall using heuristic functions $\phi_m(z/L)$, $\phi_h(z/L)$, and $\phi_y(z/L)$ for momentum, heat and moisture. In these functions, the altitude $z$ is measured against the Monin-Obukhov length

$$L = \frac{(\rho c_p T) v_{\tau}^3}{\kappa g \dot{q}_{\text{wind}}''}, \qquad (55)$$

which is positive when the cold dune is warmed by a stable ABL at night ($\dot{q}_{\text{wind}}'' > 0$), and negative when an unstable ABL cools off sand during the day. In Equation 55, $g$ is the gravitational acceleration, $\kappa \simeq 0.41$ is von Kármán's constant, $T$ is absolute temperature, and $\dot{q}_{\text{wind}}''$ is the convective thermal flux at the surface. (Note that, because we count $\dot{q}_{\text{wind}}'' > 0$ into sand, the expression of $L$ in Equation 55 has opposite sign than most conventions for the ABL.) For example, field tests in Kansas in 1968 conformed to

$$\phi_m = \begin{cases} 1 + a_{mN} z/L & , L > 0 \\ 1/(1 - a_{mD} z/L)^{1/4} & , L < 0 \end{cases} \qquad (56)$$

with coefficients $a_{mN} \simeq 4.8$ and $a_{mD} \simeq 19.3$ (Wyngaard, 2010). At low altitudes such that $z \ll L$, all three functions tend to 1, and the log-law remains unchanged. Such is the case for our relatively short weather station (Figure 10c).

The Monin-Obukhov analogy was meant for a fully developed ABL established over vast, uniform, flat expanses, like the horizontal desert floor surrounding our sparsely populated barchan dunes. However, it is questionable whether its mass transfer predictions apply to the dune itself during the day, or whether it is altogether valid under stratified, weak turbulence at night (Mahrt, 2014; Optis et al., 2014).

As Figure 11a shows, ambient moisture changes negligibly with altitude at the weather station, $\Delta Y_a \simeq -0.01 Y_a$, in contrast with the noticeable difference in ambient temperature at the same location (Figure 11b). Therefore, the moisture flux through the hard, relatively more impermeable desert floor should be very small. However, our measurements reveal that it is significant
a suitable bottom boundary condition is an invariant $\Omega$ at depth, $\partial\Omega/\partial x^*|_{\text{depth}} = 0$. In this case, using equations derived in Appendix C, such boundary condition (BC) may be written

$$ \left. \frac{\partial Y}{\partial x^*} \right|_{\text{depth}} = Y[1 + Y(M^* - 1)] \left. \frac{\partial T^*}{\partial x^*} \right|_{\text{depth}} \left[ \frac{1}{T^*} \ln \left( \frac{1}{\text{RH}} \right) + \frac{T_{A_a} T_{\text{st}}}{(T - T_{A_0})^2} \right]. \quad (40) $$

We impose it as the bottom of our integration domain, set for simplicity at the deepest sensor. To show predicted variations with depth and time, and to compare measurements and predictions, we produced a movie available in the supporting information, with typical frames in Fig. 4. This analysis also underscores the coupling of spatio-temporal variations of $T$ and $Y$. Making $f_T \equiv 0$ to eliminate it from PDE (26) would suppress the multiple observed inflections in the depth profile of $Y$, and altogether fail to capture data.

Finally, deeply buried water from antecedent rains, reported in Table 1 of Louge et al. (2013), matters to spatio-temporal variations of adsorbed moisture on sand grains nearer the surface. Its effects are captured by BC (40) at the bottom of our near-surface domain without liquid water. This BC can also be used to match our model of the near-subsurface to solutions for the deeper, moister domain below, for example those of Kamai and Assouline (2018); Assouline and Kamai (2019); Shao et al. (2021), or our own integration of Richards' equation (Richards, 1931; Louge et al., 2013).

However, because unsaturated porous media are governed by a history-dependent and hysteretic relation between capillary pressure and liquid volume fraction (Xu & Louge, 2015), definitive models of deeper liquid-infused regions are more challenging than our own. Conversely, their measurement techniques of water content are well established (Kizito et al., 2008; Kidron & Kronenfeld, 2020).

In short, until now, spatio-temporal variations of adsorbed water had not been measured in the liquid-free near-surface region of hyper-arid dunes with enough accuracy or spatial resolution. Doing so allowed us to show that this region could be analyzed without arbitrary closures. The principal remaining challenge is to model the deeper unsaturated region holding free liquid, where capillary effects continue to defy rigorous predictions.

## 11 Evaporation at the grain scale

The capacitance probe detects phenomena that other instruments with thicker measurement volume, such as microlysimeters (Boast & Robertson, 1982) or frequency-domain reflectometry sensors (Kizito et al., 2008), cannot discern. In particular, during periods of aeolian sand transport, the capacitance probe returned values of $Y$ belonging to two distinct envelopes switching rapidly from one to the other.

Figure 9 shows that the difference $\Delta Y$ between these envelopes decreases exponentially with depth, that the lowest envelope coincides with the ambient vapor mass fraction $Y_a$, and that the probe hardly ever discerned any value of $Y$ in between the envelopes. These observations suggest that strands of saltating particles meandering over the dune dried the surface intermittently, thereby sending rapid waves of $Y$ downward on time scales too small for water vapor to come to equilibrium with the thin film of water adsorbed on sand grains.

To support this explanation, we briefly summarize the underlying physics of evaporation at the grain scale, and we show that the characteristic time $\tau$ for the film's return to an equilibrium is governed by slow first-order reaction kinetics, rather than by diffusion to the surface of individual grains. Then, we conduct a linear stability analysis of Eq. (21) in the limit of short time scales. The calculations are consistent with the observed exponential decrease of the wave amplitude with depth. The corresponding length scale then lets us estimate $\tau$, thus providing insight into the drying process for these sands.
above the porous sand surface (Figure 10b), which typically traps more moisture during long periods without rain (Louge et al., 2013). Consequently, the windward toe of the dune likely constitutes the leading edge of a turbulent mass transfer boundary layer (TBL) driven by forced horizontal convection, developing below an ABL subject to vertical natural convection when unstable, and stratified otherwise.

We consider such boundary layer in Appendix H. To establish its Reynolds number $Re_X \equiv \rho U \chi / \mu$ at the probe, where $\mu \approx 1.8 \times 10^{-5}$ kg/ms is the air dynamic viscosity, we calculate the distance $\chi$ from the toe to the buried probe along the instantaneous wind direction. As the dashed line in Figure 10d shows, $Re_X$ almost always exceeds the critical value $Re_{X,crit} \approx 3 \times 10^5$ for transition to a TBL, except from sunset to midnight, when the lower ABL intermittently behaves as a viscous sublayer or as log-law turbulence (Figure 11c).

To make the TBL analysis tractable, we invoke the integral equations of the boundary layer (Kays & Crawford, 1980), and adopt the analogy between turbulent heat and mass transfer to estimate a Stanton number capturing wind-driven mass transfer. We account for exchanges at the upper boundary of the TBL by including a vapor flux $\dot{\Sigma}_{\infty}''$ counted >0 into the ABL aloft. Equation H6 then predicts the net dimensionless flux, which we equate to what we measured in Equation 54,

$$ \Sigma^* = -\Psi_s^* - \frac{p^*}{T_s^*} v_s^* (Y_s - Y_a) = a_m \frac{\rho^* U^* (Y_s - Y_a)}{\mathrm{Re}_X^{1/5} \mathrm{Sc}^{0.4}} + \Sigma_{\infty}^*, \quad (57) $$

where $\rho^* \equiv \rho / \rho_{st} \approx p^* / T_T^*$ is the relative air density in the TBL, and $\mathrm{Sc} \equiv (\mu / \rho) / D \approx 0.60$ is the Schmidt number. As Figure 10 shows, this prediction requires a value of $a_m$ about an order of magnitude smaller than a typical laboratory correlation $a_m \approx 0.0307$ for the TBL over a flat plate (Kays & Crawford, 1980), suggesting a much thicker TBL than expected. However, this disagreement may not be surprising, since Reynolds et al. (1958) fixed $a_m$ for heat transfer (not mass transfer) at Reynolds numbers within $9.2 \times 10^4 < Re_X < 1.27 \times 10^6$, which are below values in Figure 10d, except during quiet periods at night. Another possible reason for the discrepancy may be a possible gradual evolution of $Y_s$ along the dune, which we did not measure. Lastly, to gauge whether the apparent TBL thickening may be associated with the presence of strong natural convection during the day, Figure 10d also compares $Re_X$ with the square root of a Grashof number $Gr_L^{1/2} \equiv \rho [gL(T_s - T_T)/T_T]^{1/2} L / \mu$, which amounts to a Reynolds number based on $L$ and speed driven by buoyant air rising from the hot surface at $T_s$ and the colder measured ambient temperature $T_T$. However, because $Gr_L^{1/2} < Re_X$ during much of the day, natural convection is unlikely to dominate mass transfer in the TBL.

At night, fitting our data to Equation 57 also suggests a strong upward vapor flux $\Sigma_{\infty}^*$ at the top edge of the TBL, perhaps due to a moisture inversion above the dune, whereby dry air aloft lurks above a more humid layer closer to the surface (Mahrt, 2014). However, because the turbulent behavior of the stable night-time ABL is complex (Mortarini et al., 2016; Optis et al., 2014), we do not know whether this interpretation has merit. In future, detailed profiles of $Y_a$, $T$, and $U$ versus elevation should be simultaneously established above a capacitance instrument recording surface moisture flux to inform the detailed behavior of the mass transfer ABL, especially at night.

## 14. Conclusions

In this work, we deployed a unique capacitance instrument of fine vertical resolution to measure spatiotemporal variations of water mass fraction adsorbed on sand grains in the hygroscopic regime without free liquid, as well as stratigraphy of the solid volume fraction, within the first 30 cm from the surface of a hyper-arid dune. We simultaneously recorded temperature profiles there, as well as ambient conditions, net solar radiation, and wind speed and direction. We also characterized relevant sand properties, including the van der Waals isotherm that relates water mass fraction in the bulk solid to relative humidity around grains at equilibrium.

Unlike heat, which propagates mainly through the sand contact network by conduction, water vapor diffuses and advects through the tortuous interstitial pore space between grains. As a result, its governing equations are more complicated and lead to richer physics. Accordingly, we derived a PDE describing the unsteady advection-diffusion of water vapor in that space without free liquid, coupled with an ODE capturing evaporation kinetics of water adsorbed on solid grains.

We used this framework to interpret data. First, we considered long diurnal time scales, during which grains have ample opportunity to achieve thermodynamic equilibrium with the surrounding relative humidity. Here, the
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Detail of Fig. 3 showing two envelopes of the vapor mass fraction Y during a period of aeolian transport. Circles are data of sensors shown vs dimensionless time $t^*$. Thin lines indicate when the signal jumps from one envelope to the other. For each sensor, the lower Y corresponds to aeolian sand drying the surface. The superimposed lines are smoothing spline fits used to evaluate $\partial Y/\partial t^*$ in the analysis. The solid and dashed lines marked “ambient” are $Y_a$ measured at the top and bottom hygrometers of the weather station, respectively. Inset: the resulting mean difference $\Delta Y = (t_1^* - t_0^*)^{-1} \int_{t^*=t_0^*}^{t_1^*} (Y_u - Y_l) dt^*$ between the upper $Y_u$ and lower $Y_l$ envelopes of Y vs dimensionless depth $x^*$ (circles) with $t_0^* \approx 1.52$ and $t_1^* \approx 1.63$. The line is the model of Eq. (51) with $\ell^* = 0.28 \pm 0.05$.</figcaption></figure>
source term in the solid phase ODE effectively augments the unsteady term of the water vapor PDE. The result is an apparent vapor diffusion that is nearly three orders of magnitude slower than without any solid. In simpler terms, each dry sand layer hinders vapor diffusion by adsorbing moisture before letting it proceed to the next layer.

Another attribute of the long-term equilibrium behavior is a strong coupling between depth profiles of water vapor mass fraction and temperature. As the latter evolves on a diurnal basis, the non-linear coupling produces multiple inflections in moisture mass fraction that have no counterpart in linear diffusive systems. It originates in the dependence of the amount of adsorbed water on interstitial relative humidity, which varies exponentially with temperature through the saturation pressure. We considered two other transport processes driven by strong temperature gradients, namely the Soret effect and natural convection. However in our case, both proved negligible.

The porous matrix created by well-rounded packed sand also allows a slow “seepage” advection to compete with diffusion, or to augment it. To evaluate its velocity at each measurement time step, we transformed the coupled governing equations into a single ODE along depth, which we solved as a two-point boundary value problem subject to the measured vapor mass fraction just below the surface and at the deepest probe sensor that consistently yielded a stable signal. We then selected the seepage velocity that presented the smallest least squares error between ODE prediction and data.

Separately, we showed that seepage is caused by gentle variations in dune topography that induce streamwise surface static pressure oscillations decaying exponentially with depth, thereby creating Darcy pressure gradients that effectively amplify small topographical deformations. As such, seepage changes from vapor inhalation to exhalation along the surface, and it depends on wind direction.

Then, we applied the governing equations to the other limit where water mass fraction varies on short time scales. This let us elucidate a peculiar moisture wave phenomenon, whereby disturbances in vapor mass fraction on the dune surface propagate rapidly downward into its sands. A linear stability analysis of the governing equations identified the origin of these evanescent waves. Whereas subsurface moisture transport on the decimetric scale operated in an effective diffusion limit, we showed that the microscopic size of the interstitial space between grains was too small for the evaporation process itself to be diffusion-limited. Accordingly, we derived an estimate for the characteristic time of first-order drying kinetics from the exponential decay of the moisture wave amplitude with depth. Here, data suggested a drying process at the grain scale that is limited by the kinetics of a reaction with activation energy larger than, but on the order of, the latent heat of vaporization.

Finally, the depth profiles of water mass fraction *Y* allowed us to determine the time-history of the moisture flux into the ABL. To our knowledge, this constitutes the first direct measurement of this quantity on a hyper-arid surface. We registered three unexpected observations. First, when the flux was, as anticipated, nearly proportional to the difference (*Y*<sub>s</sub> − *Y*<sub>a</sub>) between vapor mass fraction at the surface and in the ambient, it was about an order of magnitude smaller than predicted for a turbulent mass transfer boundary layer developing from the dune toe, or by a Monin-Obukhov correlation of the ABL. Second, during aeolian sand transport, the flux was effectively decoupled from a vanishing (*Y*<sub>s</sub> − *Y*<sub>a</sub>). Third and most intriguing, although the flux slowly decreased in unisson with (*Y*<sub>s</sub> − *Y*<sub>a</sub>) during the stable ABL at night, it was augmented by a nearly invariant upward moisture flux aloft. To elucidate these observations in future, detailed profiles of vapor mass fraction *Y*, wind speed *U* and temperature *T* in the ABL should be recorded just above probes measuring *Y* and *T* in the subsurface.

The capacitance instrument and modeling framework that we presented in this article have wider relevance than hyper-arid dunes. In industry, they can be deployed, for example, to predict vapor contamination through pharmaceutical powders (Louge et al., 2021). Because the probe can reliably detect very low moisture, it could also be used as “ground-truth” for satellite remote sensing over sand seas (Bürgi & Lohman, 2021; Myeni et al., 2019; Zribi et al., 2014), or to look for scant water in future space applications (Davidsson & Hosseini, 2021; Honniball et al., 2020).

Studies of the Earth life-sustaining Critical Zone have underscored the importance of water exchange at interfaces among the zone's compartments (Sprenger et al., 2019). In temperate regions, the top soil often behaves as an unsaturated porous medium, where liquid water content is reliably measured (Kidron & Kronenfeld, 2020; Kizito et al., 2008), albeit on a relatively coarse spatial resolution. There, the principal challenge is to specify the
During our two consecutive days of field data collection, the top ten sensors of our capacitance probe detected no free liquid water within the first 15 cm below the surface. In this case, the exchange of moisture conformed to the “hygroscopic” regime, in which capillary forces are insignificant, and water molecules either form a thin layer on the surface of grains, or exist as a vapor in air.

The amount of moisture bound to grains in the size class of diameter $d$ is measured by the mass fraction $\Omega_d$ of water that the class holds. In the vapor phase within the interstitial space between packed grains, the mass fraction of water is $Y$. Given enough time, $\Omega_d$ and $Y$ reach an equilibrium, denoted by the subscript $e$, in which the rate of evaporation from grains is balanced by condensation at a given temperature. Appendix B provides the form of the corresponding “isotherm” of size class $d$. Once integrated over the particle-size distribution, Eq. B1 yields the isotherm for the average mass fraction $\Omega_e$ that Louge et al. (2013) recorded in terms of the relative humidity $RH_e$ in a humidity chamber.

In general, $\Omega$ is a macroscopic measure of the thin coating of water arranged as molecular layers on grain surfaces. As Eqs. (1)-(2) indicate (Shahraeeni & Or, 2010), his coating has mean thickness

$$ \ell_w \approx \left(\frac{d_{2,3}}{6}\right) \left(\frac{\rho_p}{\rho_w}\right) \frac{\Omega_1}{T^{*1/3} \ln^{1/3}[p_{\text{sat}}/(pM^*Y_e)]}, \quad (41) $$

where $\rho_p \approx 2630 \text{ kg/m}^3$ and $\rho_w \approx 997 \text{ kg/m}^3$ are, respectively, sand and water material densities, $d_{2,3} \approx 300 \text{ µm}$ is the Sauter mean diameter that arises from the integration over the particle size distribution, and $\Omega_1 \approx 0.0013$ (Louge et al., 2013). In this field campaign, the probe recorded $0.0011 < Y_e < 0.017$, or $120 \text{ nm} < \ell_w < 330 \text{ nm}$, corresponding to 400 to 1100 layers of water molecules of individual 0.3 nm thickness (Yeşilbaş & Boily, 2016).

With so many layers, it is reasonable to neglect interactions with the solid as water molecules escape from grains or adsorb on them. In this case, grain evaporation can be regarded as a non-equilibrium process releasing water molecules at the interface between interstitial air and the thin water coating, mitigated by a reverse reaction of vapor condensation (Koffman et al., 1984). Even under controlled laboratory conditions, the coupling of vapor and liquid phases, or possibly local variations of temperature near the gas-liquid interface (Jafari et al., 2018), make it challenging to model this process (Persad & Ward, 2016). Because we must also contend with grains of varying size, shape, surface geometry, impurities and composition, immersed in a random interstitial space, we resort to a rudimentary model that is inspired from the Hertz-Knudsen-Schrage (HKS) theory (Marek & Straub, 2001). In this view, the net vapor mass flux at the saturated gas-liquid interface,

$$ \dot{W}_{\text{HKS}}'' = \frac{2\kappa_e}{2 - \kappa_c} \sqrt{\frac{\text{MW}_{\text{H}_2\text{O}}}{2\pi\hat{R}}} \left[ \left(\frac{\kappa_c}{\kappa_e}\right) \frac{p_{\text{sat}}}{T_{\ell}^{1/2}} - \frac{p_{\text{H}_2\text{O}}}{T_v^{1/2}} \right], \quad (42) $$

is built upon the prediction of the kinetic theory of gases for the water mass flux $p_{\text{H}_2\text{O}} \sqrt{\text{MW}_{\text{H}_2\text{O}}/(2\pi\hat{R}T)}$ crossing a control surface under a Maxwellian velocity distribution function (Vincenti & Kruger, 1965). In this expression, $p_{\text{H}_2\text{O}} \approx pYM^*$ is the partial pressure of water vapor in the interstitial space, $\hat{R} \approx 8.314 \text{ J/mole.K}$ is the fundamental gas constant, and the “accommodation coefficients” of evaporation $\kappa_e$ and condensation $\kappa_c$ represent the probability that a water molecule colliding with the interface undergoes its respective first-order phase transition (Y. Q. Li et al., 2001).

For the tight porous medium of the dune with grains of relatively high heat capacity, it is reasonable to assume that the temperature $T_\ell$ at the gas-liquid interface and its interstitial vapor counterpart $T_v$ are equal to the local value $T$ of the “heat bath” that
water retention curve, a history-dependent relation between water volume fraction and capillary pressure subject to “return-point” memory hysteresis, which is predictable for known pore networks (Xu & Louge, 2015), but is typically modeled instead with the Van Genuchten (1980) empirical correlation, see for example, Assouline and Kamai (2019), or with local closures involving meniscus curvature (Gray & Miller, 2005) that elicit controversy (Baveye, 2013).

For hyper-arid soils without free liquid, measurements of small amounts of water adsorbed on solids lacked until now the spatial resolution and accuracy that could reveal the essential physics of water transport. With such instrument at hand, we have shown that a model free of arbitrary closures can capture observations and record the exchange of moisture at the soil-atmosphere interface. Recent improvements in the stability of the instrument's processing electronics suggest that its range can be extended to greater water content (Louge et al., 2021), thereby enabling measurements into the more humid vadose zone with a single instrument.

## Appendix A: Multiple Capacitance Probe

Louge et al. (1997) developed capacitance instruments to measure density and dielectric signature of snow. Recently, Louge et al. (2021) improved their data reduction to reach unprecedented accuracy in pharmaceutical applications. Our design is based on a similar principle (Figure 1). The probe records the impedance of the medium between any one of 15 sensors and a “target” held at a constant reference voltage. To achieve this, conductors carrying a sensor voltage are surrounded by a “guard” driven at precisely the same potential by an independent circuit. By guiding electric field lines that penetrate sand at the probe face, the guard focuses the extent of the probe's measurement volume originating from the sensor, thereby avoiding external interference known as “stray” capacitance. Because the guard is also connected to the outer conductor of a high-quality coaxial cable connecting probe and processing electronics, it also shields the wire and all electrical components carrying the sensor voltage. Therefore, because there can be no charge accumulation in the cable, the latter's capacitance does not perturb the measurement, unlike conventional bridge circuits. Consequently, this technique allows precise detection of extremely small capacitances.

For the invasive design in Figure 1, symmetry produces circular electric field lines that are shed from sensor to target in planes perpendicular to the long axis of the probe. Their center lies on the target conductor at a distance $x_c$ from the centerline and with radius $R$ satisfying $x_c^2 - R^2 = a^2$, where $a \approx 12.6$ mm is the half guard width (Louge et al., 1996). These field lines delimit a measurement volume of $W \approx 5$ mm vertical spatial resolution set by the sensor height $W$, small enough to resolve steep gradients of $\Omega$ perpendicular to the free surface. In the probe of Figure 1, each sensor has respective outer and inner distances $b \approx 11.8$ mm and $c \approx 4.8$ mm from the lance's centerline. It is separated from the grounded target by a thin guarded strip that absorbs the singular voltage jump from AC guard to DC target. As Louge et al. (1996) calculated, this makes the measurement volume penetrate sand no farther than $(a^2 - c^2)/2c \approx 14$ mm from the probe. It also confers each sensor a capacitance $C_0 = \epsilon_0 \ell_c = \epsilon_0 (W/\pi) \ln\{(b+a)(c-a)/[(b-a)(c+a)]\} \approx 38 \ 10^{-15}$ F when the probe is exposed to air or, equivalently, $\ell_c \approx 4.3$ mm.

Such tight dimensions are achieved by deploying all conductive surfaces in a thin printed-circuit board (PCB) consisting of four layers. The top one is where sensor, guard and target are exposed to sand. To limit wear, it is covered with a thin TA1YO PSR-4000BN solder mask. The next layer is a copper film held at guard voltage, designed to shield the back of all sensors. The third layer contains wires of resistance $\lesssim 1 \Omega$ bringing sensor voltages to a location where each can be individually soldered to one of 15 wires assembled into a guarded coaxial cable ~1.5 m linking probe and processing electronics. Finally, the bottom PCB layer is another guarded copper film shielding these sensor wires. Our PCB was printed by ADVANCED CIRCUITS.

Coincidentally, the thin solder mask, which equilibrates quickly to the local interstitial vapor mass fraction $Y$, also lets the probe discern rapid changes in that quantity, even if sand lags in achieving an adsorbed mass fraction $\Omega_c$ at equilibrium with $Y$. Although $\Omega$ is the primary variable determining the recorded complex effective dielectric constant, adsorbed moisture on the film, which is closest to sensor and ground, can affect the electric field, thereby acting as an impedance in series with the sand's, regardless whether $\Omega$ has reached $\Omega_c$. We noticed this behavior when the phase lead of the bare probe in moist air disappeared after it was briefly dried in an oven. While it is unclear whether the probe detects rapid variations in $Y$ quantitatively, we exploited such capability to reveal the signal bistability shown in Figure 3, which we analyzed in Section 12.
our instrument measured. Then, because the vapor is saturated at the interface, RH ≡ 1, an expression similar to Eq. (7) sets the mass fraction there, $Y_{sat} \approx p_{sat}/(pM^*)$.

Meanwhile, at equilibrium, evaporation and condensation balance ($\dot{W}_{HKS}'' = 0$), and molecular diffusion makes the interstitial vapor fraction Y uniform and equal to the value $Y_e$ recorded by the isotherm of Eq. (B1) for each size class. Therefore, the expression ($\kappa_c/\kappa_e$) $Y_{sat}$ is identified with $Y_e$, and Eq. (42) becomes

$$ \dot{W}_{HKS}'' \approx \frac{2\kappa_e}{2 - \kappa_c} \sqrt{\frac{\mathrm{MW}_{H_2O}}{2\pi \bar{R}T}} pM^* (Y_e - Y) \equiv k_m (Y_e - Y), \quad (43) $$

which serves as definition of the rate constant $k_m$ of a balanced first-order reaction of evaporation returning the system to equilibrium.

Because the interstitial space may also possess local gradients of Y, vapor diffusion on the pore scale could also play a role in setting the flux at the film's interface. For example, Philip (1964) and Shahraeeni and Or (2010) assumed that the transport of vapor to grain asperities is dominated by diffusion. (Section 12 will show that such assumption is not consistent with the depth scale of subsurface waves that we observed).

In general, as in the combustion of a solid particle (Kanury, 1975), there are two limiting regimes. In the “kinetic limit”, $k_m$ is small relative to $\rho D/d$, so diffusion of co-efficient D erases gradients of Y in a pore with interstitial air of local density $\rho$. Conversely, in the “diffusion limit”, the evaporation rate of Eq. (43) is much more intense than the local diffusive flux.

As Appendix G shows, the dominant process setting the net vapor flux $\dot{W}_s''$ at the grain surface is the slowest of the two. Once it is identified, one can infer the rate of change of the mass fraction $\Omega_d$ for each size class. Assuming for simplicity that grains of material density $\rho_p$ are spherical,

$$ \rho_p \nu \frac{\partial \Omega_d}{\partial t} = - \frac{6\nu \dot{W}_s''}{d} \equiv \frac{\rho(1-\nu)}{\tau_d} (Y - Y_e), \quad (44) $$

which defines the characteristic reaction time $\tau_d$ for the size class of diameter d. Upon integration over the particle size distribution in Appendix B, Eq. (44) yields the ODE that governs the rate of change of the overall moisture mass fraction $\Omega$ adsorbed on bulk sand of solid volume fraction $\nu$,

$$ \frac{\partial \Omega}{\partial t} = \frac{\rho(1-\nu)}{\rho_p \nu} \left( \frac{Y - Y_e}{\tau} \right). \quad (45) $$

In the next section, we conduct a stability analysis to determine the origin of the subsurface waves of Y revealed by our instrument, and to calculate the characteristic depth of their decay. In turn, this depth sets bounds for possible values of $\tau$ and, as Appendix G shows, conclusively places evaporation at the grain scale in the kinetic limit.

## 12 Evanescent waves

In section 7, we considered the limit when variables change slowly enough that the water film of mass fraction $\Omega$ remains in equilibrium with the interstitial Y. When instead Y varies on a time scale shorter than $\tau$, $\Omega$ does not change appreciably, so the right side of Eq. (21) becomes small. This situation arises during aeolian transport (Dupont et al., 2013), when clustered strands of dry particles meander over the dune, alternately forcing the surface Y to adopt the ambient value $Y_a$, or letting it quickly readjust to its higher original level without them. As Fig. 9 shows, recorded mass fractions then switch between these two extreme states, and the resulting difference $\Delta Y$ decreases exponentially