Evaporation from undulating soil surfaces under turbulent airflow through numerical and experimental approaches

Evaporation from undulating soil surfaces is rarely studied due to limited modeling theory and inadequate experimental data linking dynamic soil and atmospheric interactions. The goal of this paper is to provide exploratory insights into evaporation behavior from undulating soil surfaces under turbulent conditions through numerical and experimental approaches. A previously developed and verified coupled free flow and porous media flow model was extended by incorporating turbulent airflow through Reynolds‐averaged Navier–Stokes equations. The model explicitly describes the relevant physical processes and the key properties in the free flow, porous media, and at the interface, allowing for the analysis of coupled exchange fluxes. An experiment aiming to simultaneously collect the data of the boundary layer and soil evaporation in the area around the soil–atmosphere interface was conducted using a wind tunnel integrated with a soil tank. The turbulent boundary layer above the undulating soil surface was captured using high‐resolution hot‐wire anemometry, confirming the presence of recirculation zones in the undulating valleys and locally low evaporative flux. Experimental data were used to validate the extended model, and modeling results demonstrate that turbulent airflow enhances evaporation and shortens the duration of Stage I. The surface geometry significantly affects the local evaporative flux by influencing the vapor distribution, concentration gradient, and water availability at the soil surface, especially when recirculation zones form in the valleys. As a joint result of turbulence and surface undulations, the influence of wind speed on both the local and system‐level evaporation rate is restricted.


INTRODUCTION
Evaporation is a key component of soil-atmosphere interactions and is influenced by both the soil properties and knowledge of mass and energy exchange between the soil and the atmosphere. One knowledge gap is in our understanding of the influence of the soil surface, which in most cases is not hydrodynamically smooth but rather in the form of undulations at all scales. The undulations may form due to natural processes (e.g., soil erosion, alluvial riverbed) or mechanical manipulation (e.g., agricultural furrow). Unlike a flat surface, topographic variations can induce turbulence adjacent to the soil surface, changing the development of the boundary layer and further affecting the mass and energy exchanges between the soil and the atmosphere (Haghighi & Or, 2015). However, the characteristics and mechanisms of the influence of surface undulations linked with near-surface turbulence on the soil-atmosphere exchange processes during soil evaporation are not clear, being limited by both experimental techniques and appropriate theoretical models.
Particularly, the soil undulations of interest in this study are on the order of several to dozens of centimeters (macroscale surface roughness). There has been some published literature on surface undulations at this scale using laboratory experiments (Athanassiadou & Castro, 2001;Dawkins & Davies, 1981;Dobberschütz, 2014;Poggi, Katul, Albertson, & Ridolfi, 2007;Sugita & Kishii, 2002;Verma & Cermak, 1974a, 1974bZilker & Hanratty, 1979) or numerical simulations (Brutsaert, 1979;Cherukat, Na, Hanratty, & McLaughlin, 1998;Maaß & Schumann, 1996;Ross, 2008;. In terms of the representative experiments, Cermak (1974a, 1974b) used Styrofoam waves placed on a water-filled pan to simulate water evaporation from soil undulations and concluded that the mass transfer rate from an undulating surface was closely linked with the location of the recirculation zones (RZs). Sugita and Kishii (2002) used a lysimeter-wind tunnel system with a low water table to study the effect of macroscale roughness distribution on evaporation. Their work demonstrated that the variation of the roughness intervals could affect the evaporation rate. Most previous experiments focused on only a portion of the overall evaporation process (i.e., atmosphere-dominated Stage I or soil-dominated Stage II) or failed to simultaneously acquire atmospheric information. Other experiments  used wavy impermeable surfaces to explore how the surface geometry affected the atmospheric velocity field while no porous media information was obtained. Alternatively, previous modeling studies of fluid flow over undulating surfaces generally focused on the boundary layer development by different methods such as large eddy simulation (LES) and direct numerical simulation (DNS) (Maaß & Schumann, 1996;Ross, 2008). Some worked on the mathematical expression of the coupling of free flow and porous media flow over a curved interface Dobberschütz, 2014;. Previous modeling studies indicate that boundary layer development is influenced by surface Core Ideas • Turbulent airflow enhances undulating-surface evaporation and shortens Stage I. • Recirculation zones developed in the surface valleys affect the local diffusive flux. • The effect of wind speed is restricted under soil undulations and turbulence. • Hot-wire anemometry captures turbulent boundary layer above undulating soil surfaces.
geometry, yet these studies did not apply their results to evaporation. The boundary layers developed above soil surfaces can be laminar, turbulent, transition, or even more complicated depending on the surface geometry. Low-speed airflow over a surface generally develops a laminar boundary layer. This boundary layer becomes unstable as the downstream distance from the boundary layer origin increases and the boundary layer thickness grows. After a point of instability is reached, the original laminar boundary layer gradually transitions to a fully turbulent boundary layer. This onset of instability and transition is defined by the critical Reynolds number, which accounts for the influence of both the local streamwise position and the free flow speed on boundary layer transition. The addition of roughness to a smooth surface reduces the critical Reynolds number for the development of a turbulent boundary layer (Schlichting & Gersten, 2016;White & Corfield, 2006). This then causes transition to occur at an earlier streamwise position for a constant flow speed, or at a lower flow speed for a constant streamwise position. As a result, turbulent boundary layers commonly exist in nature due to the presence of roughness at a broad range of scales and long boundary layer development lengths in a geophysical environment. If the surface roughness is on a certain scale (e.g., soil undulations), an adverse pressure gradient may occur along the leeward sides of undulating surfaces, reducing the air speed and potentially inducing flow separation (Calhoun & Street, 2001;Finnigan & Belcher, 2004;. Although some studies identified the importance of atmospheric turbulence on subsurface gas transport (Levintal, Dragila, & Weisbrod, 2019;Maier et al., 2012;Poulsen et al., 2018;Pourbakhtiar, Poulsen, Wilkinson, & Bridge, 2017), few studies in hydrology, so far, have considered the interplay of turbulent airflow and undulating surfaces on soil evaporation. Particularly, Haghighi and Or (2015) conducted a study on evaporation from wavy soil surfaces into turbulent airflow. Experimentally, the soil surface temperature was measured by an infrared thermal camera and was used to estimate surface evaporative flux. The difference in the temperature and evaporative flux between the peaks and valleys were identified. However, the corresponding turbulent boundary layer was not measured to connect with the evaporation disparity. Haghighi and Or (2015) used a top-boundary-condition formula to calculate the evaporation rate. The effect of turbulence was incorporated by modifying the aerodynamic resistance term, and an explicit expression was used to describe the mean boundary layer thickness. Although appropriate for system level evaporation estimates (Huntingford, Allen, & Harding, 1995;Swenson & Lawrence, 2014;Villagarcía, Were, Domingo, García, & Alados-Arboledas, 2007), the top-boundary-condition type formulas cannot explicitly quantify the dynamic mass and energy transfer between the soil and the porous media. Hence, the formula proposed in Haghighi and Or (2015) is not necessarily suitable for a fundamental study of soil-atmosphere interactions like the work presented herein.
Given the abovementioned limitations of modeling theory and current experiments, the purpose of this study is to establish a benchmark understanding of the interlinked behavior of the atmosphere and porous media during evaporation from undulating soil surfaces under turbulent conditions. This study builds upon our previous work (Gao, Davarzani, Helmig, & Smits, 2018), in which the mechanisms and characteristics of evaporation from undulating soil surfaces under laminar conditions were investigated. Specifically, there are three main objectives in this paper: (a) testing the application of hot-wire anemometry (HWA) to measure the turbulent atmospheric boundary layer above undulating soil surfaces; (b) extending the coupled free flow and porous media flow model developed in Gao et al. (2018) by including turbulent airflow to simulate the dynamic soil-atmosphere interactions; and (c) analyzing the influence of wind speed, surface geometry, and soil properties on evaporation from undulating surfaces under turbulent conditions.

EXTENSION OF COUPLED FREE FLOW AND POROUS MEDIA FLOW MODEL
In our previous study, Gao et al. (2018) developed a twodomain model at the representative element volume (REV) scale to estimate evaporation by coupled Navier-Stokes free flow and porous media Darcy flow with appropriate coupling conditions at the free flow-porous media interface. This model describes the full nonequilibrium mass, momentum, and heat transfer processes in the whole evaporation system with undulating soil surfaces under laminar conditions and has been validated with corresponding laboratory experimental data (Gao et al., 2018). Compared with the topboundary-condition type formulas, this model avoids many parameterized variables and conditions, thereby excluding much of the uncertainty caused by parameterization. Instead, the physical processes involved in soil evaporation, such as F I G U R E 1 Two-dimensional configuration of subdomains and boundary conditions where C v is vapor concentration, T is temperature, U is mean wind speed, J is flux for T (total mass) and C v , λ is wavelength, and γ is wave amplitude. RANS is Reynolds-averaged Navier-Stokes mass flow and transport, heat transfer, and phase change, are included by the corresponding mass and energy equations. Moreover, the atmospheric airflow states and the development of the boundary layer affected by the surface geometry, the two main focuses in this study, can be explicitly described and obtained. The modeling outputs not only include the evaporation rate, but also the real-time evolution of global and local variables of interest in the system (e.g., moisture, temperature, vapor concentration, mass and thermal fluxes, etc.). Therefore, this model is particularly qualified for the purpose of understanding the soil-atmosphere interactions in this fundamental REV scale study.
In this work, we extend the model developed in Gao et al. (2018) by incorporating turbulent airflow in the free flow through Reynolds-averaged Navier-Stokes (RANS) equations, in which the Reynolds stress term is expressed by a two-equation model for kinetic energy k and energy dissipation ε, respectively (i.e., RANS with k-ε closure model). Among the methods available to simulate turbulent flow (e.g., DNS, LES, and RANS), RANS is widely used in engineering (Allegrini, Dorer, & Carmeliet, 2014;Defraeye, Blocken, & Carmeliet, 2010) because of its high efficiency. Fetzer, Smits, & Helmig (2016) and Mosthaf (2014) used RANS with a zero-equation closure model due to its simplicity, convergence success, and relatively low computational cost. However, zero-equation closure models cannot properly account for the historic effects of turbulence (e.g., diffusion and convection of turbulent energy) (Trautz, 2015). Instead, the RANS with k-ε two-equation closure model is a promising option to simulate turbulent airflow, although extra near-wall flow treatment is necessary. Figure 1 shows the two-dimensional conceptual configuration and boundary conditions for this study. The size of the computation domain and the boundary conditions are decided based on the experimental setup. The whole domain consists of two subdomains: the lower porous media (i.e., soil) and the upper free flow region (i.e., the near-soilsurface atmosphere). The two subdomains are separated by a two-period cosine curve with an aspect ratio (AR) defined by 2γ/λ (γ and λ are the amplitude and wavelength, respectively), representing the undulating soil surface. In the free flow, a single gas phase is assumed, which is composed of two components (dry air and water vapor). The gas flow is assumed incompressible and gravity is neglected. In the porous media, two phases, (gas and liquid water) coexist, and the gas is composed of two components (water vapor and dry air). The dissolution of air in water is neglected, and the two fluids are immiscible. Based on these assumptions, a set of equations describing the fluid flow, vapor transport, and heat transfer is defined separately in the free flow and porous media. This mathematical model was discussed in detail in Gao et al. (2018). Only the extended part (i.e., the description of turbulent airflow and the change of interfacial conditions due to turbulence) is introduced below. The additional equations applied in the porous media, at the interface, and other boundary options are referred to in Gao et al. (2018).
As mentioned above, the turbulent airflow in the free flow region is described by the RANS with k-ε closure model. Besides the mass and momentum balance equations, another two partial differential equations are added. One is for the turbulent kinetic energy, k, and the other is for the energy dissipation, ε.
The mass and momentum balance equations in the free flow are where the superscript "ff" denotes the free flow region and the subscript "g" denotes gas phase. ff g (m s −1 ) and ff g (Pa) are the gas flow velocity and pressure in the free flow region. ρ g (kg m −3 ) and μ g (Pa s −1 ) are the gas density and dynamic viscosity, depending on the temperature and the fraction of water vapor in the gas (their values are updated in real time during calculation). μ g,T is the eddy viscosity in the RANS equations defined by μ g, T = ρ g μ 2 ε (3) where c μ (c μ = 0.09) is a constant. The equations for k and ε are where σ k , σ ε , c ε1 , and c ε2 (σ k = 1.0, σ ε = 1.3, c ε1 = 1.44, c ε2 = 1.92) are model constants and P k is the production term given by Under turbulent flow conditions, the component mass balance equation in the free flow region is modified by including a turbulent diffusion term (Fetzer et al., 2016): where ff v is the mass fraction of water vapor in the gas phase), ff v (m 2 s −1 ) is the diffusion coefficient in the free flow region, and ff v,T (m 2 s −1 ) is the eddy diffusivity caused by turbulence. The eddy diffusivity is usually determined from the turbulent Schmidt number, Sc T , a dimensionless number describing the ratio between the rates of the turbulent transport of momentum and mass: Likewise, the energy balance equation in the free flow region is also revised by including a turbulence-induced term for conductivity (Fetzer et al., 2016): where T ff (K) is the temperature, c p,g [J (kg ⋅ K) −1 ] is the heat capacity of the moist air, and λ g [W (m ⋅ K) −1 ] is the thermal conductivity of gas mixture. λ g,T [W (m ⋅ K) −1 ] is the eddy thermal conductivity formed due to turbulence. Similar to the turbulent mass transfer, the dimensionless Prandtl number, Pr T , defined by the ratio of momentum and heat transfer eddy diffusivity, is used to determine the eddy thermal conductivity. The Sc T and Pr T of gas are about 0.78 and 0.71, Vadose Zone Journal respectively, at atmospheric pressure and temperature of ∼298 K (Mosthaf, 2014): Under turbulent conditions, the near-wall region is different from the mainstream, complicating the calculation of the flow field by turbulent models (Pope, 2001). Although it is possible to modify the k-ε turbulent model so that it can describe the flow in near-wall regions, this is not always desirable because of the accompanying high-resolution requirements. Instead, an analytical expression, known as a wall function, is typically used to describe the near-wall flow (Kuzmin, Mierka, & Turek, 2007;Lacasse, Turgeon, & Pelletier, 2004). One limitation of this approach is that the wall function assumes that the wall is impermeable, whereas porous surfaces allow for mass exchange between the porous media and the free flow through pores at the interface. However, under turbulent flow, a viscous sublayer exists adjacent to the permeable surface where k-ε equations are not applicable. In this case, we assume that the pore-scale mass exchange between the porous media and the free flow at the interface has a negligible influence on the REV-scale evaporation. Hence, the wall function (Equation 11), together with the normal continuity of mass, stresses, and fluxes (Appendix C), is applied at the soil surface. The wall function is described by (Kuzmin et al., 2007): where κ is the von Kármán constant (κ = 0.41), B is a constant set to 0.52 by default, δ + w is the dimensionless thickness of the viscous sublayer, and n is the outward pointing unit vector normal to the boundary.
A summary of the equations to be solved and the primary variables is presented in Table 1.

EXPERIMENTAL SETUP AND RESULTS
To provide a preliminary understanding of the behavior of evaporation from undulating soil surfaces under turbulent airflow conditions, we first conducted an experiment using a wind tunnel and soil tank system (Gao et al., 2018), as shown in Figure 2. Unique in this work, a HWA system was applied to measure the turbulent boundary layer above the undulating soil surfaces. To the authors' knowledge, HWA has not previously been used in such a laboratory experimental system for a soil evaporation study. The main advantage of this method of wind speed measurement compared with sonic anemometers is that the sensor wire is tiny so as to measure the wind speed at one point relatively precisely, and the wire probe can be easily mounted on a high-spatial-resolution traverse. The computer-automated traverse can then accurately locate the sensor to designated positions, allowing for the measurement of wind speed profiles at high spatial resolutions. In this section, the experimental materials, methods, and protocols are discussed.

Development of experimental apparatus
Experiments were conducted using an open-return wind tunnel (Culler & Farnsworth, 2019) in the Experimental Aerodynamics Laboratory at the University of Colorado Boulder (Aerolab, UCB). The wind tunnel has a square test section (76 × 76 cm) and can maintain a steady wind at both low and intermediate speeds. To accommodate the soil tank, the bottom panel of the wind tunnel was redesigned to allow the top of the soil tank to align precisely with the bottom panel. Hence, the floor of the wind tunnel served as the surface of the soil tank.
The soil tank was constructed with acrylic glass with a length of 60 cm, height of 30 cm, and width of 9 cm. Water content was continuously monitored throughout the tank using 22 ECH2O EC-5 moisture sensors (Decagon Devices, accuracy = ±3%). Sensors were inserted into the tank in 10-cm increments horizontally and 5-cm increments vertically ( Figure 2). Data were collected using a series of five-channel continuous data loggers (Decagon Services, ECH2O System). The soil tank was packed with a uniform silica sand, Accusand #50/70 (effective sieve number, particle diameter D 50 = 0.23 mm, Unimin Corporation). This sand has been widely characterized and used in the previous laboratory experiments (Deepagoda, Smits, Ramirez, & Moldrup, 2016;Gao et al., 2018). The porosity of the packed soil tank was ∼0.32.
The speed profiles of the boundary layer were measured by the HWA system. In this paper, "velocity" denotes a vector composed of components in different directions, whereas "speed" denotes the "magnitude" of one velocity component. The HWA system includes a five-channel AN-1003 constant temperature anemometer (CTA) from A.A. Lab Systems, a straight single-sensor wire miniature probe (Model no. 55P11) with the corresponding probe support from Dantec Dynamics, and a software package. The probe has a single 5-μm-diam., 1.25-mm-long plated tungsten wire as a sensor. The probe was mounted on a traverse with the wire perpendicular to the streamwise flow direction, allowing for the time varying measurement of the streamwise speed normal to the sensor. F I G U R E 2 Schematic of experimental setup. HWA is hot-wire anemometry.

Experimental procedures
Soil moisture sensors were tested and calibrated prior to installation into the corresponding locations of the soil tank. The tank was then wet-packed into a two-period cosine surface configuration with #50/70 silica sand in incremental 2-cm layer to achieve a uniform bulk density (Sakaki & Illangasekare, 2007). The packed sand was saturated below the valleys with blue-dyed water to allow for visualization of drying during the experiment. Soil undulations were then built with fully wet sand. The amplitude and length of one-period undulation were 3.6 and 30 cm, respectively. During packing, some free water was left above the valley surface to keep the undulations saturated via capillary rise. Because the characteristic length of the #50/70 sand is ∼12 cm (Lehmann, Assouline, & Or, 2008) based on the sand properties (Gao et al., 2018), the undulating soil surface was assumed to be close to full saturation at the start of the experiment (also confirmed by soil moisture sensors). Prior to the start of the experiment, the soil tank was covered with plastic wrap to prevent evaporation. The side and the bottom walls of the tank did not allow for fluid flow, and the valve at the bottom of the tank was shut off so that no water was supplied at the bottom boundary. After packing, the soil tank was placed on a weight scale (Sartorius Model 11209-95, range = 65 kg, resolution = ±1 g) to continuously monitor the water loss from the system and hence the system-level evaporation rate. The tank and scale setup were precisely fit into the custom opening of the bottom floor of the wind-tunnel test section and adjusted to make sure the tank surface was aligned flush with the floor (Figure 2). The small gap between the opening and the tank was sealed with tape to prevent improper air ventilation. A motorized, computercontrolled traverse was installed from the ceiling of the windtunnel test section. The traverse could move horizontally along the mainstream through the opening at the wind tunnel ceiling. The single-sensor wire probe of the HWA system was F I G U R E 3 Schematic of the coordinate system above the soil surface. The red dashed lines are the test columns then mounted on the traverse, allowing for precise positioning of the sensor to the predefined measurement points. The HWA system was then implemented to measure the mean wind speed field point by point above the undulating soil surfaces. Figure 3 shows the reference coordinate system for the HWA measurement points. The x axis was set at the tank surface. In the horizontal direction (x axis), 12 positions were selected from 0 to 60 cm spaced every 5 cm. In the vertical direction (y axis), the measurement range was from −7 cm (the value of the lowest point of the soil surface on a vertical axis) to 10 cm in the free flow region. Piecewise increments in the vertical direction were chosen in order to guarantee the accuracy of measurement: Δy = 1 mm (h < 1 cm); Δy = 2 mm (1 cm < h < 2 cm); Δy = 5 mm (h > 2 cm), where h is the height above soil surface and Δy is y increment. All the test points were set at the lateral center profile of the tank. For each vertical measurement column (the red dashed lines in Figure 3), the probe sensor was initially placed at the first measurement point, ∼2 mm away from sand surface. The probe was then moved up under the control of the computer and stopped at the next predefined point for 2-3 s to measure the wind speed. In this way, the sensor collected data point by point until reaching the highest test point (y = 10 cm) in each column. Afterwards, the probe shifted to the next column and repeated the same steps until finishing the measurements at all predefined test points. Six hundred and fifty-nine test points in total were collected in this experiment.
To start the experiment, the wind tunnel was set to a high wind speed of 6.0 m s −1 (±0.1 m s −1 ), ensuring the development of a fully turbulent boundary layer over the soil surfaces. The HWA measurements were first taken at the column 108 mm upstream of the tank edge (x = −108 mm in Figure 3) to verify the generated boundary layer. In the vertical direction, the test points were spaced from y = 0 cm to y = 10 cm with the same piecewise increment defined above. The solutions of both laminar and turbulent boundary layer equations were used to fit the data to verify the flow state. The fitting results are presented in Section 3.3. Afterwards, the plastic wrap that covered the top of the tank F I G U R E 4 Measured dimensionless speed profile at 108 mm upstream of the soil tank. Power-law turbulent and Blasius laminar boundary layer solutions are compared by fitting the measured profile. u (m s −1 ) is the measured wind speed at different locations, U (m s −1 ) is the maximum wind speed in the mainstream (i.e., 6 m s −1 in this experiment), and δ* (m) is the displacement thickness was removed, any remaining water located above the valleys was drained by a syringe, and the evaporation experiment commenced. During evaporation, the soil moisture and weight of the soil tank were monitored every minute. The wind speed profiles were measured by HWA point by point. The experiment lasted 6 h due to the limitation of time and laboratory.

Experimental results and model verification
Prior to presenting the boundary layer profiles above the undulating soil surfaces, the profile measured at 108 mm upstream of the tank is given by Figure 4 as compared with different fitting methods. Power-law turbulent boundary layer profiles with different index (β) values and laminar Blasius solutions are compared with the observations, in which β = 4.225 is estimated from the measured data (Appendix D). A detailed mathematical description of the two boundary layers can be found in Appendix D. Figure 4 illustrates that the power law matches with the data much better than the Blasius solution, indicating that a turbulent boundary layer was developed along the wind tunnel floor. According to Cheng (2007), a power-law index (β) varying from 4 to 12 is recommended, and 2 to 6 is suggested in some extreme cases such as macroscale bed roughness. Hence, β = 4.225 estimated from experimental data is qualified for the data fitting, consistent with the theory in Cheng (2007). Figure 5a shows the wind speed profiles above the undulating soil surfaces measured via HWA. The speed decreases with the reduction of the vertical position from 6 m s −1 in the mainstream. In the valley, however, the speed increases again approaching the soil surface, especially near the leeward side of the valley. Since the HWA with a single-sensor wire probe could only measure the average streamwise speed without direction, the changes of the speed with the vertical positions indicate the reverse airflow (i.e., the presence of RZs), thus demonstrating that HWA is able to capture the turbulent boundary layer. The simulated wind speed (Figure 5b) using the numerical model described in Section 2 shows similar trends with the measurements. However, several experimental limitations are responsible for the difference in the profile between the experiment and modeling. First, the manual surface packing results in small surface variations, preventing a perfect cosine shape like that in the model. Second, because of the point-by-point nature of HWA, it took upwards of 6 h to obtain two-dimensional speed profiles in this experiment. Although the wind-tunnel flow field is very precise, small variations in temperature and pressure in the laboratory during this period can affect the wind field development. Finally, the presence of the traverse and sensor itself in the free flow, albeit small, disrupts the wind field development. Despite these limitations, HWA was able to capture the turbulent boundary layer development, demonstrating its promise for such uses. This is especially useful, as many times nonintrusive methods such as particle imaging velocimetry (PIV) are not available or are cost prohibitive. Hot-wire anemometry offers a fairly simple yet precise alternative to gain understanding of the boundary layer. However, to overcome some of these limitations, a nonintrusive method will also be explored for future experiments. Figure 6 shows the observed speed contours obtained by kriging interpolation of the measured speed data. The low-speed area above the leeward surfaces (Surfaces S1 and S3) is caused by the RZs. The presence of the RZs, which will be expanded upon in Section 4.2, results in vapor accumulation at leeward surfaces and thus a locally low evaporation rate. This was also observed visually through the deposition of the blue dye on the soil surface (Shokri & Or, F I G U R E 6 Interpolated speed contours (m s −1 ) using the measured speed profiles in Figure 5a. S1, S2, S3, and S4 represent the four wavy surfaces, respectively. 2011). As seen in Figure A1 in Appendix E, differences in the blue dye deposition are qualitatively very obvious based on the color gradients between the leeward and windward surfaces. The darker blue surfaces correspond to locations of high evaporation rate. Connected with Figure 6, this demonstrates that locally lower evaporation appears at the leeward soil surfaces due to the RZs.
During the 6-h experiment, the weight of the soil tank and soil moisture were monitored. Figure 7 illustrates the experimental and modeled accumulative evaporation and soil moisture 6.5 cm below the top surface of the tank (∼1 cm below the undulating sand surface). The data oscillation at ∼2.5 h is due to an unavoidable stop and restart of the wind tunnel for adjustment purposes. As the soil tank was initially saturated, the accumulative evaporation amount increases almost linearly with time at a high evaporation rate (Stage I evaporation). The modeled soil moisture matches the observations well, further verifying the ability of the model to capture the evaporation behavior under turbulent conditions. This will be further discussed in Section 4.

Vadose Zone Journal
F I G U R E 7 (a) Comparison of observed and modeled evaporation amount; (b) comparison of observed and modeled soil moisture at the first-layer moisture sensor (x = 10 cm, y = −6.5 cm based on the coordinate system showed in Figure 3)  Note: U, maximum wind speed in the mainstream; AR, aspect ratio; γ, amplitude; λ, wavelength. α and n are van Genuchten parameters.

NUMERICAL ANALYSIS OF EVAPORATION FROM UNDULATING SURFACES
This section provides theoretical insights into the characteristics of evaporation from undulating soil surfaces under turbulent conditions. Seven cases, A-G, listed in Table 2 were simulated. Five variables of interest (i.e., flow state, wind speed, AR, size of sand[represented by van Genuchten parameter α and n], and soil permeability) were studied. For each scenario, the simulated porous media was initially saturated. Additionally, the vapor concentration in the free flow and the temperature of the whole system were set equal to the inflow boundary conditions (i.e., 0.006 kg m −3 and 22 • C, respectively), consistent with experimental conditions. The mathematical model was implemented using COMSOL Multiphysics 5.2a (COMSOL, 2016) based on a finite element method. Quadratic Lagrange elements were used for the dependent variables ( pm g , pm w , and ff g ), and linear Lagrange elements were applied for the other variables ( ff g , C v , and T). The whole computational domain was discretized by triangular meshes with 20,506 elements (156,545 degrees of freedom in total). Local mesh refinement was introduced at the inter-face of the two subdomains due to the possible formation of large gradients. Before the calculation, the vapor concentration field and flow field were initialized. A backward Euler time discretization was used to generate an implicit solving system. In every time step, a segregated solver was applied. The flow-concentration-heat variables ( ff g , , v , and T j ; i = w, g; j = ff, pm) and turbulent variables (k and ε) were solved in two subgroups. Each group used a damped version of Newton's method. In each iteration, a linearized version of the nonlinear system was solved using multifrontal massively parallel sparse direct solver (MUMPS) and parallel sparse direct solver (PARDISO) for each subgroup, respectively.

Laminar and turbulent flow conditions
As discussed in Section 1, near-surface turbulence can occur even at a low wind speed due to the presence of roughness at a broad range of scales and long leading edge for boundary layer development. Because an inaccurate assumption of flow state may lead to a great error in the estimation of the evaporation behavior, we first compare the application of laminar (Gao et al., 2018) and turbulent models by Cases A and B for F I G U R E 8 Effect of airflow state, comparison of Cases A and B: modeled evolution of the evaporation rate the same soil at the same low wind speed. This is an important comparison, as many models assume laminar conditions without recognition as to the influence of such an assumption. Figure 8 shows the evolution of the evaporation rate calculated based on the assumption of turbulent and laminar conditions. From a fully saturated condition as presented here, the evaporation process generally progresses in three stages: (1) the atmosphere-dominated Stage I when the evaporation rate is high and constant or slowly decreases; (2) the transition period when the evaporation rate decreases quickly; and (3) the soil-dominated Stage II when the evaporation rate is low and decreases slowly. Figure 8 illustrates that the estimated duration of evaporation stages and evaporation rate by laminar and turbulent models are significantly distinct. Based on the given modeling conditions, the evaporation rate obtained from the turbulent model is obviously higher and stays in Stage I for a shorter time compared with the laminar model. The transition stage occurs at about t = 6 d according to the zoom-in view of the evaporation rate from 0 to 15 d. After a rapid transition, the evaporation then comes to Stage II when the evaporation rate is lower than that under laminar conditions. Hence, the turbulent airflow promotes Stage I evaporation, shortens Stage I duration, and accelerates the transition to Stage II. The comparison illustrates that an inaccurate assumption of airflow states will definitely result in a miscalculation, especially for the fundamental mechanism studies. In the sections below, as an extension of Gao et al. (2018) in which laminar airflow was assumed, the turbulent model is applied to investigate the influence of turbulence on the evaporation from undulating soil surfaces for various atmospheric and soil conditions.

Influence of wind speed on local and system-level evaporation behavior
In this section, the influence of wind speed on the evaporation from undulating surfaces is considered under turbulent conditions by comparing Cases A, C, and D. Figure 9a shows the change of system evaporation rate over time under three different wind speeds from 1 to 6 m/s. Speeds were selected to correspond to low, medium and relatively high speed conditions compared to nature, in which 6 m/s was the speed applied in the above experiment. Consistent with Gao et al. (2018) for laminar airflow, increasing wind speed promotes stage I soil evaporation under turbulent airflow. However, in comparison to laminar conditions, the enhancement with an increase in wind speed is quite small under turbulent conditions, especially for high wind speeds. Interestingly, through laboratory experiments, Davarzani et al. (2014) noticed a threshold wind speed over which an increase in wind speed had a limited effect on the evaporation rate, whereas no explanation was provided. This can be understood by our model in which turbulence is incorporated by further investigation of the diffusive flux. According to Fick's law, diffusive flux is controlled by the thickness of viscous sublayer where vapor diffusion dominates, the vapor concentration difference across the viscous sublayer and the diffusion coefficient related with turbulence under turbulent conditions. Figure 9 (b) shows the distribution of the upward diffusive flux along the undulating soil surface at t = 3 and 30 d. Time t = 3 d corresponds to a time when all three test cases are still in Stage I evaporation. At this time, the peaks and valleys contribute to the evaporation almost equally, as shown by their similar diffusive flux. However, when it comes to the third day after the evaporation starts, the water saturation at the peaks has dropped close to 20%, while the valleys still have a 80% saturation (Figure 9c). This indicates that the peaks have already lost their potential capacity to keep a high diffusive flux due to the water deficiency, and thus the flux at the peaks will decrease quickly (as seen Figure 9b from 3 to 30 d). By contrast, during the fast moisture exhaustion at the peaks, the valleys are still able to continue to provide water to the surface (Figure 9c), and they present a larger diffusive flux than the peaks over time (Figure 9b). It determines that the valleys are the main contributors during most of the evaporation period. Therefore, the influence of wind speed will be limited to the valleys, as the peaks dry quickly.
However, the diffusive flux at the valleys does not increase a lot with the increase in wind speed, as seen in Figure 9b, which can be understood from three aspects. First, RZs are formed in the valleys for the AR set in these three test cases, where the reverse airflow leads to accumulation of water vapor at the leeward surfaces, resulting in locally lower diffusive flux than that at the windward surfaces, especially at the start of RZs (i.e., separation points; Figures 9b and 9d). Thus, the thickness of the viscous sublayer is also related to the RZs, rather than just the wind speed, like in the case of laminar airflow over a flat surface. Second, under turbulent conditions, the diffusion coefficient (Equation 7) consists of two terms, one of which is associated with turbulent  8). In this sense, turbulent variables such as kinetic energy and energy dissipation are highlighted, and the impact of wind speed is somewhat undercut. Third, turbulent airflow promotes the whole-system evaporation and decreases the amount of water reaching the valleys, although slightly (Figure 9c), with the increase in wind speed, resulting in less formation of water vapor (Figure 9d). The vapor concentration difference across the viscous sublayer will decrease to some degree, which is accompanied by the reduction of the viscous sublayer thickness caused by increasing wind speed. Thus, the change in the diffusive flux is unpredictable.
In short, in the case of evaporation from undulating soil surfaces under turbulent conditions, wind speed is not the only key influencing factor for the development of boundary layers, and thus its effect on the local diffusive flux and system-level evaporation rate is somewhat restricted. As for the experiments based on a flat surface in Davarzani et al. (2014), the last two reasons, both caused by turbulence, account for the existence of a wind speed threshold after which the wind speed has reduced impact on system-level evaporation. Further, this also demonstrates that an accurate assumption of airflow states is significant, especially for REV-scale fundamental studies.

Influence of aspect ratio
This section discusses the influence of the undulating soil surface AR by comparing test Cases A and E. As the AR is defined by the ratio of the depth and length of one soil undulation, a larger AR corresponds to a steeper undulating surface. An AR = 0.27 is selected, as it is close to the experimental design, and AR = 0.15 represents a smoother while natural surface.
The comparison of the evaporation rate for the two ARs ( Figure 10a) shows a rather negligibly faster evaporation rate from the steeper undulating case, but only during Stage I. This can be understood from two aspects. First, large ARs increase the effective surface area in contact with water (Figure 10b) through capillarity and thus enhance evaporation mainly in Stage I. However, due to the fast drying of the peak area under turbulent conditions, the enhancement is not significant. Second, RZs form in the valleys for AR = 0.27, where the reverse airflow causes the local accumulation of water vapor at the leeward surfaces as shown in Figure 11a. By comparison, no RZ is identified for the smoother surface (AR = 0.15), and thus the vapor spreads more evenly along the soil surface as shown in Figure 11b. The vapor distribution in the two cases determines that the surface vapor concentration gradient is Vadose Zone Journal F I G U R E 10 Effect of aspect ratio (AR), comparison of Cases A and E: (a) evaporation rate, (b) water saturation at soil surface at t = 3, 30 d, and (c) upward and (d) horizontal diffusive flux along surface at t = 3, 30 d F I G U R E 11 Vapor concentration (g m −3 ) in the free flow and water saturation in the porous media at t = 3 d lower upward and higher horizontally for the steeper surface than for the smoother surface. In turn, the concentration gradient determines that compared with the smoother surface, the steeper surface has a smaller upward diffusive flux and a larger horizontal flux, as shown by Figures 10c and 10d, respectively. Further, this determines slight difference in the evaporation rate caused by AR, since evaporation rate equals the integration of surface diffusive flux.

Influence of permeability
Permeability is an intrinsic property of the porous media, which determines the continuous flow ability of liquid and gas through the media. Under laminar conditions, higher permeability results in a longer Stage I evaporation for both undulating and flat surfaces (Gao et al., 2018;Mosthaf, Helmig, & Or, 2014). In comparison, under turbulent conditions, changes in permeability have little to no effect on the Stage I duration (Figure 12a). The reason for the equivalent evaporation rate in Stage I is the adequate water supply at the soil surface for both test cases, which can also be demonstrated by the equal upward diffusive flux at t = 3 d, shown in Figure 12b. Alternatively, as mentioned in Section 4.1, the turbulent airflow shortens the first stage of evaporation. It thus restricts the influence duration of permeability and results in an almost equivalent Stage I period. Therefore, the assumption of laminar conditions may misestimate the influence of soil properties.

Influence of size of sand
Based on the discussion in the above sections, the turbulent condition shortens the first stage of evaporation and accelerates the transition to Stage II. This phenomenon highlights the importance of water capillarity to the soil evaporation, which is closely associated with the size of sand. The drying characteristic of a sand can be described by the van Genutchen parameters (van Genuchten, 1980): α (m −1 ), which is related to the inverse of the air entry suction; and n (dimensionless), which is connected with the pore size distribution. The van Genutchen parameters of the #50/70 sand used in this experiment are α = 5 m −1 and n = 7 (Gao et al., 2018), and for a coarser sand #12/20, they are approximately 8 m −1 and 11 (Wallen, Smits, Sakaki, Howington, & Chamindu Deepagoda, 2016), respectively. The influence of the combined effect of the air suction and the pore size is numerically studied by comparing the evaporation from #50/70 and #12/20 sand. The evaporation rate shown by Figure 13a illustrates that under turbulent conditions, the coarse-sand system evaporates quickly and thus enters transition and Stage II rapidly. Figure 13b also shows that the peak area of a coarse sand loses its contribution early, whereas the valley area plays an important role in the system evaporation due to water availability. Hence, Stage I evaporation is not only shortened by turbulent airflow, but it may even disappear in the case of a coarse sand. An inaccurate assumption of the airflow states may lead to wrong estimation of evaporation stages.

CONCLUSIONS
With the goal of understanding the influence of atmospheric turbulent airflow on the evaporation from undulating soil surfaces, we initiated an explorative study using both modeling and experimental approaches. The coupled free flow and porous media flow model that was developed and testified in Gao et al. (2018) was modified by incorporating turbulent airflow in the free flow through RANS equations, in which the Reynolds stress term is expressed by a two-equation model for kinetic energy k and energy dissipation ε, respectively (i.e., RANS with k-ε closure model). The prominent feature of this coupled model compared with the commonly used parameterized top-boundary-condition-type formulas is that it uses the most basic REV-scale equations to explicitly describe the physical processes during soil evaporation (e.g., mass and heat transfer, phase change, etc.), as well as the atmospheric airflow states and their development as influenced by the surface geometry. In that sense, the uncertainty caused by simplified parameterization can be mostly excluded, and it makes the coupled model particularly appropriate for fundamental studies. In particular, the evaporation rate, as well as other global and local information in the evaporation system (e.g., moisture, temperature, vapor concentration, mass and thermal fluxes, etc.), can be obtained in real time by this method.
A preliminary laboratory experiment was conducted by integrating a soil tank within a wind tunnel to initiate this study. Hot-wire anemometry was for the first time applied in the soil tank-wind tunnel system to measure the wind field for the laboratory study of soil evaporation. The development of the boundary layer above the undulating soil surfaces was successfully measured at high spatial resolutions and RZs were clearly captured in the valley area. Recirculation zones caused locally low-level evaporation at the leeward soil surfaces, which was observed by less dye deposition at leeward surfaces compared with the windward surfaces. Besides, the collected soil moisture and accumulative evaporation amount in the experiment, as well as the HWA-measured wind speed profiles, matched well with the simulation results from the extended coupled model, allowing us to have an initial understanding of the influencing factors using numerical modeling.
From the perspective of modeling, the influence of wind speed, AR, permeability, and size of sand on evaporation from undulating soil surfaces under turbulent conditions was investigated. Compared with laminar conditions, atmospheric turbulent airflow enhances the system-level evaporation rate and shortens the duration of the Stage I due to fast loss of water, especially at the undulation peaks. The undulation valleys play a major role during most of the evaporation period due to their vicinity to the receding water table. First, increasing wind speed has negligible influence on the local evaporation at the peaks while slightly affecting the valleys. The limited influencing area leads to limited contributions of increases in wind speed to the evaporation rate of the whole system. In the meantime, RZs are prone to develop in the undulating surface valleys with certain ARs, which can affect the development of the viscous sublayer and further influence the local diffusive flux. Especially under turbulent conditions, the vapor diffusion within the viscous sublayer is also associated with turbulent features such as kinetic energy and energy dissipation. Under the intertwined circumstances including both soil undulations and atmospheric turbulent airflow, wind speed is no longer particularly important, as is commonly assumed for laminar airflow over a flat surface. Therefore, the influence of wind speed on undulating-surface evaporation under turbulent conditions will be restricted. Second, the surface AR directly determines the presence and extent of the RZs, which play a significant role in the vapor distribution and the diffusive flux at the soil surfaces. Thus, the AR affects evaporation through influencing both the vertical and horizontal diffusive flux out of the soil surface, as well as the surface water availability. Third, permeability shows little impact on the duration of Stage I under turbulent conditions, and further, Stage I may even vanish in the case of a coarse sand.
Additionally, in terms of the experiment, although it is possible to capture the turbulent boundary layer using HWA, the single-sensor wire probe used in this study could only measure the mean streamwise wind speed. Therefore, further follow-on study is warranted, using particle image velocimetry to capture a full instantaneous boundary layer profile with flow directions. Also, the soil tank used in this experiment was originally designed for a single-period wavy surface presented in Gao et al. (2018). Hence, the locations of the soil moisture sensors were not ideal (i.e., close to the surface of peaks and valleys). A new soil tank will be designed based on the locations of interest to measure moisture and temperature. The soil surface undulations will be built at the tank surface rather than inside the tank to avoid the influence of the side walls of the tank on the airflow development in the valleys. A series of experiments will be carried out considering various variables (e.g., wind speed, undulation space) to obtain more data to further improve the numerical model and understanding of the characteristics of evaporation from undulating surfaces under turbulent airflow.
We expect the planned follow-on work to be useful for improving the models for evaporation prediction on the engineering scale and similar research problems involving interactions between the free flow and porous media flow.

CONFLICT OF INTEREST
The authors declare no conflict of interest.

A. Equations for porous media flow
The mass balance equations of these two phases are (Bear, 2013) where the superscript "pm" denotes the porous media. S i , ρ i , and u i (i = w, g) are the saturation, density, and velocity of the i phase ("w" for water and "g" for gas), respectively. ϕ is the porosity of the porous media, p c (Pa) is the capillary pressure between the water and gas, and f vw (kg m −3 s −1 ) is the watergas phase change rate during evaporation, defined by Bixler (1985) as where b is a fitting parameter dependent on equilibrium time of phase change (Gao et al., 2018); S rw is the residual water saturation; R (J mol −1 K −1 ) is the universal gas constant; T (K) is the temperature; M w (kg mol −1 ) is molecular mass; ρ v,eq (kg m −3 ) is the density of vapor at equilibrium; and ρ v (kg m −3 ) is the vapor density at any condition. The vapor density at equilibrium is defined by according to Kelvin's equation, where ρ v,s (kg m −3 ) is the saturated vapor density. The saturated vapor density is estimated empirically (Campbell, 1985): Detailed explanation of the phase change rate term can be found in Bixler (1985), Gao et al. (2018_ Halder, Dhall, and Datta (2011), and Smits, Cihan, Sakaki, and Illangasekare (2011.
The velocity in Equation A1 is determined by Darcy's law: where K int is the intrinsic permeability of the soil; k ri (i = w, g) is the relative permeability of the i phase. The van Genuchten model is used for the water retention curve and relative permeability (van Genuchten, 1980).
where S ew is the normalized water saturation defined by S ew = (S w − S rw )/(1 − S rw − S rg ); S rw and S rg are residual saturation of water and gas, respectively; and α and n are van Genuchten fitting parameters (m = 1 − 1/n).
where l is a van Genuchten fitting parameter. The mass balance equation for the two components in the gas is described by Bear (2013) as where pm v and pm d are the effective diffusion and dispersion coefficient in the porous media, respectively. The former is determined by where D v is the gas phase diffusion coefficient (Campbell, 1985) and τ is the tortuosity of the porous media, which is estimated according to Millington and Quirk model (Millington & Quirk, 1961). The latter is considered as a function of product of dispersity and velocity. The energy balance equation is expressed by ( ρ p ) ef f = ϕ w ρ w p,w + ϕ g ρ g p,g + (1 − ϕ) ρ s p,s (A11) where c p,i (J kg −1 K −1 ) is the heat capacity of the i phase (i = w, g, s); λ i (J kg −1 K −1 ) is the thermal conductivity of the i phase (i = w, g, s); and λ eff (J kg −1 K −1 ) is the effective thermal conductivity of all the phases including water, gas, and solid grain, which is defined by (Nield & Bejan, 2006): λ ef f = ϕ w λ w + ϕ g λ g + (1 − ϕ) λ s (A12) The subscript "s" in Equations A11 and A12 denotes the solid phase. Q s and L in Equation A10 describe the energy losses from the soil tank and latent heat during water vaporization, respectively (Moradi, Smits, Massey, Cihan, & McCartney, 2015).

B. Coupling conditions at the interface
Since different model concepts are used in the free flow and porous media, extra interfacial conditions are necessary to couple the two subdomains (Mosthaf et al., 2011) based on the assumption of local mechanical, chemical, and thermal equilibrium.
1. Continuity of total mass flux Considering there is only a gas phase that exchanges between the free flow region and the porous media, the continuity of total mass flux in the normal direction of the interface should be ( ρ g ff g where n denotes the normal vector of the interface.

Continuity of normal stress
The normal stress at the interface is continuous:

C. Other boundary conditions
As shown in Figure 1, the left, right, and bottom boundaries of the porous media are closed where Neumann conditions are applied − pm pm = 0 ( = w, g) We assume that the center line of the free flow region has the largest wind speed, which is considered as the top boundary, and thus the following no viscous stress and normal flow conditions should be satisfied: The left boundary of the free flow region is set as the inlet, where a constant average velocity, constant vapor concentration, and temperature are specifically defined. The right boundary of the free flow region is the outlet, and the manometer pressure is zero.

D. Turbulent and laminar boundary layer solution 1. Power-law turbulent boundary layer
where u (m s −1 ) is the streamwise time-mean flow speed, and U (m s −1 ) is the maximum flow speed in the mainstream. y is the surface-normal distance. The power-law index β is determined by β = 2 − 1 ; = δ * δ * * where H is shape factor; δ * and δ ** are is displacement thickness and momentum thickness, respectively, calculated from experimental data in this paper. 2. Blasius laminar boundary layer (Sanders, 2014) Blasius transformed the second-order partial differential equations governing the growth of the laminar boundary layer on a flat plate under two-dimensional steady and incompressible assumption to a nonlinear, third-order ordinary differential equation by introducing a similarity variable η and dimensionless stream function ψ: where y is the surface-normal distance, x is the streamwise position, υ is the kinetic viscosity, and U is the average freestream speed. This nonlinear, third-order ordinary differential equation can be solved precisely by numerical methods, resulting in a series of η, f(η), f ′ (η) = u/U, and f ′′ (η). The calculation shows that η ≈ 4.9 for u/U = 0.99 when y = δ, where δ is the overall boundary layer thickness. The displacement thickness δ* is related with the overall thickness δ and for the flat surface is δ* = 0.35δ. Therefore, the surface-normal distance y normalized by displacement thickness δ* can be expressed as based on η: δ * =