Obstacles, Interfacial Forms, and Turbulence: A Numerical Analysis of Soil–Water Evaporation Across Different Interfaces

Exchange processes between a turbulent free flow and a porous media flow are sensitive to the flow dynamics in both flow regimes, as well as to the interface that separates them. Resolving these complex exchange processes across irregular interfaces is key in understanding many natural and engineered systems. With soil–water evaporation as the natural application of interest, the coupled behavior and exchange between flow regimes are investigated numerically, considering a turbulent free flow as well as interfacial forms and obstacles. Interfacial forms and obstacles will alter the flow conditions at the interface, creating flow structures that either enhance or reduce exchange rates based on their velocity conditions and their mixing with the main flow. To evaluate how these interfacial forms change the exchange rates, interfacial conditions are isolated and investigated numerically. First, different flow speeds are compared for a flat surface. Second, a porous obstacle of varied height is introduced at the interface, and the effects the flow structures that develop have on the interface are analyzed. The flow parameters of this obstacle are then varied and the interfacial exchange rates investigated. Next, to evaluate the interaction of flow structures between obstacles, a second obstacle is introduced, separated by a varied distance. Finally, the shape of these obstacles is modified to create different wave forms. Each of these interfacial forms and obstacles is shown to create different flow structures adjacent to the surface which alter the mass, momentum, and energy conditions at the interface. These changes will enhance the exchange rate in locations where higher velocity gradients and more mixing with the main flow develop, but will reduce the exchange rate in locations where low velocity gradients and limited mixing with the main flow occur.


Introduction
Throughout our natural and engineered surroundings, fluids flow in a multitude of different fashions, transporting mass, momentum, and energy to form the dynamic environment we know. In specific flow domains, mathematical descriptions of these flows have been developed. For free-flow domains (e.g., atmosphere, surface water, or mechanically convected flows), the Navier-Stokes equations were developed, describing the momentum transport in viscous fluids. For flows within a porous media flow domain (e.g., groundwater flow, filtered flow), Darcy's law was developed and adapted to describe many complex cases from nature and industry. Unfortunately, many realistic applications, spanning spatial scales and fields of interest, cannot be described as one isolated domain, as flow dynamics in adjacent domains often exhibit a coupling to the dynamics of their neighbors. For example, within proton-exchange membrane (PEM) fuel cells, supply channels and porous diffusers should be designed such that reactive gases efficiently supply the porous catalyst layers with reactant, maintain a certain level of humidity, and remove any excess of liquid water produced (Siegel 2008). In urban environments, impermeable surfaces reduce evaporation rates, where heat and air pollution can stagnate due to disturbed flow paths. This causes the urban heat island (UHI) effect, where cities see larger increases in temperature than rural areas, reducing livability and adding economic strain (Allegrini et al. 2015). Further applications, such as the salinization of agricultural land (Jambhekar et al. 2015), flow through oil filters (Iliev and Laptev 2004), rocket cooling (Dahmen et al. 2014), and nuclear waste storage (Masson et al. 2016), all require an understanding of the coupled behaviors between porous media flow and free-flow domains. In a majority of these applications, the domain interface will not be flat and may include flow obstacles or interfacial forms. These interfacial shapes will cast different flow structures along the interface, modifying the flow conditions and the exchange processes across the surface. In this work, we perform a numerical investigation of the non-flat interface between a porous media flow and a free-flow domain, focusing specifically on the exchange processes between a partially saturated bare soil and an atmospheric turbulent flow.
Mass, momentum, and energy transport at the soil-atmosphere interface depend on the dynamics of flow both above, below, and at the surface interface. In the free-flow domain, the development of boundary layers depends on the momentum, energy, and compositional content of the atmosphere, as well as the shape of the domain interface. Viscous boundary layers, developing at the interface, define the diffusive distance across which water vapor diffuses to the ambient atmospheric conditions (Shahraeeni et al. 2012). In the subsurface, the supply of liquid water held by capillary forces will also limit the location of the vaporization front (Shahraeeni and Or 2010). Considering soil-water evaporation from a flat surface, the development of boundary layers is fairly well understood considering both laminar and turbulent free-flow regimes (Haghighi et al. 2013;Davarzani et al. 2014;Trautz et al. 2015), heterogeneous soils (Fetzer et al. 2017b;Vanderborght et al. 2017), and soil grain roughness (Fetzer et al. 2016). If a surface exhibits curvature and undulations or is covered in porous or impermeable obstacles, the formation of sublayers at the surface will change significantly as flow structures will develop, complicating energy and mass transport. Evaporation from surface undulations and around flow obstacles has been investigated experimentally in Verma and Cermak (1974), Haghighi and Or (2015a), Haghighi and Or (2015b), Gao et al. (2018), Trautz et al. (2018), and Gao et al. (2020). These experimental works have further compared their results to numerical models based on surface renewal theory (Haghighi and Or 2015b), or by evaluating either a single flow domain (Trautz et al. 2019), or both (Gao et al. 2020). Each of these works has drawn conclusions about evaporation from a non-flat interface, but has investigated specific experimental geometries. With the results and model concepts presented in these works in mind, individual geometrical conditions are isolated in this work and their effects on exchange rates analyzed numerically.
Numerical investigations of flow systems can be performed using many approaches. For applications with coupled flow regimes, various methods exist, in particular, onedomain methods (Brinkman 1949), or mixed-theoretic methods (Varsakelis and Papalexandris 2011). Further, Yang et al. (2018Yang et al. ( , 2019 each investigate flow conditions in turbulent channel flows adjacent to porous structures. These works resolve the full geometry of all structures, simulating turbulent heat and momentum transport within both domains using Reynolds stress models (RSM) as outlined in Speziale et al. (1991). Due to the geometric and computational limitations of these methods, two-domain models are used in this work. Both flow domains are modeled using individual model concepts, and these domains are coupled using coupling conditions (Ochoa-Tapia and Whitaker 1995;Layton et al. 2002;Baber et al. 2016;Weishaupt et al. 2019). In this work, a numerical investigation is performed on the representative elementary volume (REV) scale, using the two-domain concepts introduced in Mosthaf et al. (2011Mosthaf et al. ( , 2014 and expanded to turbulent conditions in Fetzer et al. (2016).
This compositional, non-isothermal, and two-domain approach is further developed using turbulence models and higher-order discretization methods in the free-flow domain. In addition, the interface is expanded to include shapes and obstacles. These obstacles and shapes will cast different flow structures depending on their size, spatial parameters, proximity to other obstacles, and shape. The flow structures that develop around these obstacles and shapes will alter the mass, momentum, and energy conditions along the interface, which will modify the exchange processes across the interface. In areas where the flow structures provide high velocity gradients at the interface and additional mixing with the main flow, higher evaporation rates will occur. In areas where the flow structures collect evaporated water vapor, due to low velocity gradients and limited mixing with the main flow, evaporation rates will decrease. The aim of this study is to analyze how these strongly coupled processes function under different isolated topographical conditions.
Using the model concepts and numerical methods outlined in Sects. 2 and 3, evaporation is first compared for different flow conditions along a flat surface to outline how nearinterface free-flow conditions alter exchange rates. Then, an analysis of interfacial obstacles and forms is performed. The height of a single rectangular partially saturated porous obstacle is first investigated and changes to the evaporation rate, up and downstream of the obstacle, are evaluated. Following this, a square obstacle containing no liquid water is introduced with varied flow properties. With different flow parameters in the porous obstacle, the free-flow flow structures will change, and the changes in evaporation rate upstream and downstream are evaluated. Next, a second obstacle is introduced. Due to their proximity, the flow structures cast by the upstream obstacle will merge with those of the downstream obstacle. The flow structures between the obstacles will then change their shape according to the spacing between the obstacles, and the evaporation rate across the surface between the obstacles is evaluated. Lastly, the shapes of these obstacles are modified from rectangular forms to triangular, sinusoidal, and sawtooth forms. Each of these forms casts a different set of flow structures up and downstream of their scope. These flow structures are compared, and their effect on evaporation rates discussed.

Model Concept
Using a fully implicit monolithic scheme, continuum-scale models describing turbulent free flow are coupled to REV Darcy scale models simulating multiphase flow and transport in a porous medium using interface coupling conditions. This two-domain concept is based on the work developed in Mosthaf et al. (2011) andFetzer et al. (2016) and is coupled using the sharp interface techniques described in Fetzer et al. (2017a). In this section and the following, the model concepts used in each domain, their coupling, and their numerical implementation are described.

Free-Flow Domain
To simulate the mean effects of turbulence in the turbulent free-flow domain ( ff ), the nonisothermal compositional Navier-Stokes equations are modified using Reynold's decomposition and time averaging to develop the RANS equations. The RANS equations are able to simulate the mean effects of turbulence in a free flow; each transported quantity is split into its fluctuating value ( ′ g ) and its time-averaged mean value ( g ). Fluctuating terms alone average to zero over time ( p � g = � g = 0 ), where products of these fluctuating terms do not ( � g � g ≠ 0) (Versteeg and Malalasekra 2009). This decomposition and time averaging create an additional term, the Reynolds stress tensor, t , and replace each transported value with its time-averaged mean value ( , p g , T, w g ). For simplicity, each time-averaged primary variable will be written plainly throughout this work. The Reynolds stress tensor, g , arising from the convective term, contains fluctuation terms that cannot be resolved on this scale of analysis. Using the theory posed by Boussinesq (1877) and discussed in Schmitt (2007), this term is assumed to act in strict analogy to viscous stresses, adding a so-called turbulent eddy viscosity ( g,t ) to the viscosity in the diffusive term. For non-isothermal and compositional models, the same addition to the diffusive term needs to be made for the mass fraction and the temperature, added in the form of the turbulent eddy diffusivity ( D g,t ) and the turbulent eddy thermal conductivity ( g,t ). Although the turbulent dynamics in the mean flow are important in this work, the viscous near-wall flow conditions must be accurately resolved in order to evaluate the interfacial exchange processes. The free-flow model used must be able to resolve both of these flow conditions.
The RANS equations used in this work are developed according to the following assumptions: (i) single-phase gaseous flow (g), consisting of components ( ), water (w), and air (a), (ii) the simulated gas phase can be described as a Newtonian fluid, neglecting dilation, (iii) local thermodynamic equilibrium, specifically Fickian diffusion, (iv) turbulent stresses act isotropically, (v) turbulent pressure fluctuation does not significantly modify gas density, and (vi) turbulent stresses can be treated similarly to viscous stresses via the Boussinesq assumption (Schmitt 2007).
Here, the time-averaged velocity vector and pressure, g and p g , are the primary variables and the turbulent shear stress, g,ff,t , is defined, neglecting dilatation, as: The variable g,t , or the turbulent eddy viscosity, is solved using turbulence closure models. is the identity matrix, and g and g are the density and the viscosity in the gas phase, respectively, both of which depend on the compositional make up of the gas and the temperature.
To evaluate the transport of water (w) in the gas phase, the following component mass balance equation (4) is used: Here, the mass fraction of water in the gas, X w g , is the primary variable, and the turbulent diffusive flux, w Diff,g,ff,t , is defined here as: The term w g is the mole fraction of component ( w ) in the gas phase ( g ), or the quotient of the mass fraction and the molar mass ( X w g ∕M w ), and g,mol is the molar density of the gas phase. D w g , is the binary diffusion coefficient, dependent on the temperature, and D w g,t is the turbulent eddy diffusivity, determined as a function of the turbulent eddy viscosity: Values for the turbulent Schmidt number, Sc t , are found in Wilcox (2006). Finally, to model energy transport, the following energy balance equation (7) is used: Here, the temperature T is evaluated as the primary variable. The variables u e g , h g , and g are the specific internal energy, the enthalpy per component ( ) in the gas phase, and the thermal conductivity of the gas phase, respectively. The diffusive term, h g Diff,g,ff,t , represents the enthalpy transport due to component diffusion. The turbulent conductive flux, Cond,g,ff,t , is defined here as: g,t is the turbulent eddy conductivity, determined as a function of the turbulent eddy viscosity: Values for the turbulent Prandtl number, Pr t , are found in Flesch et al. (2002).
In order to calculate the turbulent eddy viscosity ( t ) , turbulence models of varying complexity have been developed, typically classified by the number of additional partial differential equations they require. Within the coupled model environment, zero-equation models, such as those proposed by Baldwin and Lomax (1978) and based on mixing length theory, have been used in the past with adjustments for interface soil grain roughness (Fetzer et al. 2016). For this investigation, these models and the one-equation models proposed by Spalart and Allmaras (1994) are insufficient as they rely on geometrical and mean flow assumptions (Schlichting and Gersten 2006). Two-equation models are more complete in that they incorporate a spatial and temporal scale, and do not rely on geometrical assumptions made to define the mixing length (Pope 2006). These models aim to simulate turbulence on these scales with two of three turbulent variables: k, the turbulent kinetic energy, , the dissipation, or , the turbulence frequency. There are an abundance of variations of these models, but within the coupled environment used in this work, the k-and k-models have been implemented. Kmodels, such as those developed in Launder and Sharma (1974), are popular and have been used in the analysis of coupled systems (Gao et al. 2020), but alone are suitable only for higher Reynolds numbers ( t >> g ). To counter this, wall functions, describing flow dynamics in a viscous sublayer, are introduced, but these wall functions are often not valid in complex geometries. The k-model, originally proposed in 1942 by Kolomogorov (Wilcox 2006), and expanded to the format used in this work in Wilcox (2008), does not use wall functions or other damping functions, but uses grid specific boundary conditions to simulate near-wall effects. This model was also used in Yang et al. (2019) when comparing coupled two-domain model results to a fully resolved coupled channel and porous media flow simulation using a Reynolds stress model.
In this work, the k-model will be used for all comparisons. Although all turbulence models have trouble predicting the precise dynamics of separated flows in arbitrary geometries, this model has been tested, evaluated, and adapted against turbulence validation cases, such as the backward facing step (Driver and Seegmiller 1985), which analyze separated flows similar to those seen in this work (Vescovini 2019). Due to the k-model's flexibility in comparison with the other models, and its tested resolution of near-wall viscous flow, it will be used for the remaining tests.
The additional balance equations making up the k-model are as follows: . These two variables then resolve the turbulent eddy viscosity, t , using the following relations and coefficients: These models incorporate the following wall specific boundary conditions, k wall = 0 and wall,P = 6 g y P , where y P is a discretization specific distance to wall length. For more information about this turbulence model, please refer to Pope (2006) and Wilcox (2006).

Porous Media Flow Domain
The porous media flow domain ( pm ) is modeled using a non-isothermal compositional multiphase Darcy's law and is based on the following assumptions (Helmig 1997): (i) gas (air) and liquid (water) as mobile phases (g, w), (ii) Newtonian fluids, (iii) creeping flow, (iv) a non-deformable rigid solid matrix, (v) unique capillary pressure saturation relationships (no hysteresis), and (vi) local thermodynamic equilibrium.
To balance mass in a porous medium, the following general mass balance is written as the following: where Adv, pm, , and Diff, pm, are the advective and diffusive fluxes of each component in each phase : Here, is the porosity, S the saturation, X , the mass fraction, the phase density, D pm, the effective binary diffusion coefficient for a porous medium, and the fluid phase velocity as determined by Eqs. 16 and 17. To evaluate the momentum, the multiphase Darcy's law is used, written here for isotropic porous media considering gravitational effects: where K is the intrinsic permeability, k r, the relative permeability, the phase viscosity, P the pressure, and the gravity vector. For more on the development of the mass and momentum balances, please refer to the model concept developed in Fetzer et al. (2016). In addition, for flow through porous media with higher Reynolds numbers, an additional law can be introduced to describe inertial stresses, the Darcy-Forchheimer law (Nield and Bejan 2006): where c F is the Forchheimer coefficient, set here as 0.55.
Evaporation is an energy intensive process, meaning that heat transfer within the porous media will affect the final exchange rates. To evaluate heat transfer within the porous media, we assume a local thermodynamic equilibrium, such that the temperature across phases is the same within one control volume. The energy balance used in this work is as follows: Here, the subscript (s) describes the solid phase, or matrix material. The subscript (pm) on the other hand describes the full partially saturated porous medium, as conducted fluxes are not separated into phase specific quantities under a local thermodynamic equilibrium.

Coupling Conditions
The conditions used to couple the free flow and porous media flow in this work are outlined in Mosthaf et al. (2011) and are based on the assumption of local thermodynamic equilibrium from Moeckel (1975), Layton et al. (2002).
In the free-flow domain, a tangential momentum condition is specified, representing a nonzero slip condition at the interface. This is based on the work of Beavers and Joseph (1967) and Saffman (1971) and was originally defined for the interface between a laminar flow and single phase flow in a porous medium. Its validity has been investigated experimentally in Terzis et al. (2019) and mathematically in Eggenweiler and Rybak (2020), where opportunities for improvement are outlined. In a numerical investigation regarding turbulent flows, converging and parallel, Yang et al. (2019) also investigate the validity of this condition. As converging flows in this work are, at most, limited to low Reynolds numbers, and as alternative conditions for this interfacial condition have not been fully developed, this condition is used in this work. Further, as each model concept defines momentum differently, there exists a potential pressure jump between the two domains, , which can lead to minor differences in the mole fraction, X w , and the temperature, T , at the interface. Instead, the normal stresses are held to be continuous across the interface.
The development of these conditions is outlined in Mosthaf et al. (2011) and expanded to turbulent conditions in Fetzer et al. (2016). The options for implementation over a sharp interface with no further degrees of freedom are investigated in Fetzer et al. (2017a). Each condition is shown in Table 1.

Free-Flow Domain
The free-flow domain is discretized using finite volumes aligned via the marker and cell (MAC) scheme (Harlow and Welch 1965), here referred to as a staggered grid. As depicted in Fig. 1, scalar variables are associated with cell centers. Components of the vector variable velocity are located in centers of cell edges, such that they point along the normals to the edges. For the finite volume scheme, scalar equations, such as the continuity equation (4), are integrated over the grid cells. The control volumes for the components of the vectorial momentum equation (2) are located around the respective vector variable component, as shown by the rectangles in Fig. 1.
Central difference approximations are applied to the diffusive terms in the momentum equations (2). Transporting velocities in convective terms are arithmetically averaged when integrals are built along lines in the centers of which there is no velocity degree of freedom.
The transported velocities in convective terms are treated by total variation diminishing (TVD) methods, which allow second-order accuracy while preserving monotonicity and avoiding oscillation in the solutions, see, e.g., Harten (1983) and Sweby (1984). Specifically, the formulation of Hou et al. (2012) is chosen for the transported velocities on nonuniform grids,  For r < 0 , i.e., at a local maximum or minimum, the flux limiter becomes zero and, hence, first-order upwinding is applied. A full investigation of this discretization scheme using various flux limiters is performed in Vescovini (2019). For all non-porous wall boundaries, a no-flow no-slip condition is established. Symmetry boundaries, used throughout this work as an upper boundary condition, establish a zero-gradient condition for normal fluxes, while maintaining the tangential velocity. The inflow conditions use a Dirichlet condition for a uniform velocity field, and the outflow conditions set a Dirichlet condition for the total mass balance.

Porous Media Flow Domain
The porous media flow domain is discretized using cell-centered finite volumes. The simulations are performed on a rectangular structured non-uniform grid, using a two-point flux approximation. The multiphase Darcy's law is solved using the pressure saturation formulation, with the primary variables pressure, p g , liquid phase saturation, S w , and temperature, T, which are all located at the finite volume cell center. After a minimal saturation S w,r has been reached, the cell has numerically dried out, and a primary variable switch is performed, replacing S w with X w g (Helmig 1997;Class et al. 2002). All non-interface boundaries in this work specify no-flow boundary conditions.

Interface Conditions
Within the scope of this work, the interface has been modeled using a sharp interface approach; at the interface there are no extra degrees of freedom, meaning no mass, momentum, or energy can be stored between the domains. Further, the grids must be conforming. For a further discussion of the implementation of these coupling conditions, please refer to Fetzer et al. (2017a) and Grüninger et al. (2017).
In order to provide better initial conditions for the coupled simulation runs, simulations of each domain are run with coarse initial conditions until realistic equilibrium conditions are reached. The results from these separated domains are then used as initial conditions within the coupled model context.

Implementation
The concepts presented above are implemented using DuMu x Koch et al. 2020), an open-source simulation environment based on the open-source numerical toolbox Dune (Bastian et al. 2008a, b). In addition to the standard Dune modules, Dunesubgrid (Gräser and Sander 2009) is used to split the full discretization into two conforming domains coupled across an irregularly shaped interface. The grid used is limited to rectangular elements, but allows for conforming non-uniformity. Using zonal grading, the grid is refined near the interface and in areas of interest to fully resolve the flow structures, while areas further away were left coarser. For simulations with simpler geometries, around 6000 grid elements were used, whereas for larger obstacles around 13,000 elements were used. For specifics on the grids used, please refer to the below-referenced publication module.
In this work, the models defining each domain are coupled using a monolithic approach, meaning that all balance equations are assembled into one system of equations and solved at once. Newton's method is used to linearize the nonlinear system of equations, and the system is solved using the direct linear solver UMFPACK (Davis 2004). Convergence is established, stopping the iterative Newton method, if the current iteration step n (≥ 1) holds that: with N being the number of degrees of freedom, x i , and the tolerance, r tol , ranging from 10 −5 to 10 −6 . Other publications using this two-domain approach include ) where a Navier-Stokes free-flow model is coupled to a pore-network model, and Schneider et al. (2020) where a Navier-Stokes free-flow model is coupled to a single-phase porous media flow model discretized with a multi-point flux approximation for unstructured cell-centered finite volumes.
All of the source code used to perform the simulations referenced below can be found in a DuMu x publication module (Coltman 2020) together with installation instructions. Upon further questions regarding the source code, please contact the authors.

Numerical Analysis
Using the model concept described in Sect. 2 and the discretization described in Sect. 3, numerical experiments are developed to investigate soil-water evaporation from a surface with different obstacles and interfacial forms. An illustration of these experiments is shown in Fig. 2. To begin, the free-flow velocity is varied above a flat surface and the effects on evaporation are outlined. Next, an obstacle, hydraulically connected to the porous media is introduced. The height of the obstacle is varied, and the flow field and evaporation rates upstream and downstream of the obstacle are evaluated. Following this, an dry porous obstacle is introduced. The permeability of this porous obstacle is varied, and the flow field and evaporation rates upstream and downstream of the obstacle are again evaluated. As a further expansion, a second obstacle is introduced. The spacing between these obstacles is varied, and both the flow field and the evaporation rates between the obstacles are evaluated. Finally, the effect the shape of the interfacial obstacle, formed to five different shapes, has on the surrounding flow field and the evaporation rate is investigated. Each of these obstacles casts different flow structures in the free flow, altering the conditions near the surface. These changes are investigated, and their effect on the exchange rates across the interface is explained.
In each of the following simulations, the default properties, parameters, and conditions listed in Table 2 are used unless otherwise specified.

Flat Surfaces: Varied Free-Flow Velocities
First, we investigate evaporation from a flat surface of 1m length (Fig. 3a). Similar configurations have been performed experimentally in Davarzani et al. (2014) and evaluated numerically in Fetzer et al. (2016). Evaporation rates from a porous medium are defined in this work as the mass exchange rate across the domain interface. The rates shown here represent a unit evaporation, or the height of water (mm) removed over time (day), averaged across a length of interface ( 1 m 2 ). This evaporation is split into two phases, phase I and phase II. Phase I evaporation can be identified as the initial higher evaporation rate, governed by atmospheric demand and characterized by a continuous liquid water supply at the surface . As the rate of phase I evaporation is dependent on atmospheric demand, the distance across which water vapor diffuses within the free flow, will contribute to determining the evaporation rate. This distance is directly correlated to the thickness of the mass, momentum, and energy boundary layers (Haghighi et al. 2013). As soon as the total mass loss has advanced such that air replaces the first layer of pores, phase II begins, characterized by vapor diffusion rates within the porous media, warmer surface temperatures, and a falling evaporation rate (Lehmann et al. 2008).
As shown in Fig. 3b, the evaporation rate during phase I increases and the duration of this phase decreases with increasing Reynolds number, defined here using a characteristic length of 1m, or the free-flow height given a mirrored vertical domain across the symmetry boundary. This increase in evaporation rate is due to the reduction in evaporative diffusive distance, in this case, the thickness of the viscous sublayer of the boundary layer, as well as turbulence-enhanced mixing in the free-flow. With faster mass removal, the duration of this phase is reduced.
For non-flat interfaces, the diffusive distances along the interface will vary depending on the flow structures that develop adjacent to the interface. Additionally, the available water saturation at the surface will vary according to the height of any elevated interfacial forms. With these two factors in mind, evaporation rates across non-flat interfaces may Finally, the shape of these interfacial forms is changed. Each additional interface condition will produce different flow structures adjacent to the interface. These flow structures will alter the mass, momentum, and energy sublayer thicknesses and mix differently with the main flow, altering the exchange rates across the interface include a mixture of varied phase I evaporation rates, as well as both phase I and II evaporation rates in different locations at one time.

Obstacle Height: One Partially Saturated Obstacle of Varied Height
To begin, a rectangular obstacle is introduced at the interface. This obstacle is fully hydraulically connected to the underlying porous medium and begins with a initial water saturation, S w , of 0.6. Upwind and downwind of the obstacle, there is a flat coupled interface of 1m length, and the length of the obstacle is held at 10 cm. The height of the obstacle, H, is then varied: H = 2 cm, 4 cm, 6 cm, 8 cm, and 10 cm, as shown in Fig. 4.
Introducing an obstacle to the interface will disrupt the boundary-layer conditions at the interface on both sides of the obstacle. Upwind of the obstacle, the flow will separate and converge upward to flow around the obstacle, and a small corner eddy at the upwind base of the obstacle will develop. Downstream of this obstacle, a long recirculation zone develops, along with a small corner eddy at the base of the obstacle. With increasing obstacle height, the size of these flow structures grows, casting different flow conditions along the interface. An illustration of these flow structures is shown in Fig. 5, where the upstream detached length D, the reattachment length R, and the surface length of the eddy L are shown. Two examples of these regions are displayed in Fig. 6 where the free-flow direction is drawn with arrows and the flow magnitude plotted in the background.
The coefficient of skin friction, C f , a metric used to qualitatively describe near-wall flow conditions used in turbulence validation tests cases such as Driver and Seegmiller (1985), is defined as the follows:  Downstream of the obstacle, in Fig. 7b, there is a short initial length of positive recirculation with a low velocity gradient, the downstream corner eddy. This is then followed by a length of negative recirculation with a higher velocity gradient, the downstream recirculation zone. The length from the obstacle to the point where the skin friction becomes positive after recirculation zone is called the reattachment length, R. The size of the corner eddy and the full reattachment length both increase with obstacle height, but both the corner eddy and the recirculation zone do not see significant changes in velocity gradient with increased obstacle height.
These flow structures upwind and downwind of the obstacle will alter atmospheric demand at the interface. In areas where no separation occurs, but the velocity gradient is reduced, the evaporative diffusive distance increases, reducing the evaporative demand. Within the small corner eddies with low velocity gradients, cool water vapor is circulated Depending on their location, these eddies will add a surface area for diffusive mixing with the mean flow, enhancing the demand, but also involve lower flow velocity gradients, reducing the demand. Downstream of the obstacle, the larger recirculation zones will both return air of main flow content to the surface, due to the additional mixing area across the length L, and maintain a higher velocity gradient at the surface. Both of these factors will enhance the evaporation rate.    Fig. 8, the evaporation rates across the surface length upstream of the obstacle (Fig. 8a) and downstream of the obstacle (Fig. 8b) are shown. In addition, the average evaporation rate across the surface at 4 days, or ∫ E(x)dx , is calculated where is the length of the coupled interface (1m) and E(x) is the location specific evaporation rate. These values are shown in comparison with the obstacle height. Upstream of the obstacle, the velocity gradient decreases over an increasing portion of the interface with increasing obstacle height. This corresponds to a decrease in evaporation rate with increasing obstacle height. Closer to the obstacle, within the detached distance, D, the corner eddy will cast a decreased evaporation rate due to the low velocity gradient, but as obstacle height increases, this reduction in evaporation rate will wane, due to the additional mixing area provided by the surface of the eddy. This can be seen in the local maximum near the obstacle; as the corner eddy turns clockwise, a peak is located at the right side of this eddy, where the circulated air meets the surface.
Downstream of the obstacle, a corner eddy, roughly proportional to the obstacle height, develops, where low velocity gradients and limited mixing with the main flow reduce the evaporation rate. The additional diffusive mixing area, provided by the surface of the eddy mixing with the surface of the recirculation zone, slightly increases the evaporation rate within this area, as seen by the local maximum adjacent to the obstacle. In comparison with the upstream corner eddy, evaporation rates are much lower, mainly as mixing with the main flow is limited downstream of the obstacle by the recirculation zone. Further downstream, the larger recirculation zone develops, characterized by higher near-wall velocity gradients and a large additional diffusive mixing area, L. At the reattachment length R, the highest evaporation rate is reached, where the atmospheric demand is enhanced by the recirculation of dry mean flow air to the interface.
Across the 1-m interface downstream of the obstacle, the evaporation rate is both enhanced in the recirculation zone and reduced beneath the corner eddy. For smaller obstacles (2 cm, 4 cm, and 6 cm), the area of reduction is outweighed by the enhanced effects in the recirculation, and the net evaporation rate across the interface increases with the obstacle height. After a height of around 6 cm, the area of reduction continues to increase, and Fig. 8 (a) Evaporation rates upwind of an obstacle of height H, plotted across the upstream interface length at 4 days. (b) Evaporation rates after an obstacle of height H, plotted across the downstream interface length at 4 days. Also shown in both figures is the average upstream/downstream evaporation rate at 4 days plotted against obstacle height the peak associated with the recirculation zone shifts further downstream. For this coupled area of analysis, the trend reverses; for higher obstacles, the net evaporation rate will still be enhanced, but will decrease with increased obstacle height.

Obstacle Permeability: One Dry Obstacle of Varied Permeability
Similar to the tests performed in Sect. 4.2, one obstacle is placed at the soil-atmosphere interface. This obstacle is identical in geometry to that of Fig. 4, but the saturation and flow resistance parameters within the porous obstacle are changed. This obstacle is not hydraulically connected to the underlying porous medium, separated here by eliminating any capillary action ( vg = 0.1, n vg = 100 ) and contains no liquid water ( S w = 0 ). The porosity is held at a large 0.75, and the permeability, K, of the obstacle is varied: K = 10 −7 , 10 −8 , 10 −9 , 10 −10 (m 2 ) . These obstacles do not share the flow parameters of common natural soils, but could describe other natural obstacles located at the soil-atmosphere interface such as shrubs or engineered boundaries. Momentum is evaluated within this obstacle using Forchheimer's law, as described in Sect. 2.2,Eq. (17).
In comparison with the previous section, the obstacle is much less resistant to flow; no liquid water blocks pathways within the obstacle, and the permeability is in some cases reduced. As a result, more momentum from the free-flow domain will enter the obstacle, rather than flowing around it. To visualize how this momentum transfer from the free flow to the porous obstacle affects the flow structures surrounding the obstacle, the C f [Eq. (22)], is plotted upstream and downstream of the obstacle in Figure 9. Upstream of the obstacle, increasing permeability will have only a minor effect on the flow structures. This corner eddy will have slightly higher velocity gradients and the detached distance D will reduce slightly with increasing permeability. As more air can flow through the porous obstacle, the shape of the downstream flow structures will change more noticeably. The corner eddy will increase in length with increasing obstacle permeability, therefore extending the distance to reattachment.
Upstream of the obstacle, these minor changes to the flow structures will create slight changes in the evaporation rate. As shown in Fig. 10a, the slight reductions in detached correspond to slight increases in the evaporation rate. Additionally, beneath the corner eddy, the increasing velocity gradients with increasing permeability will enhance evaporation rates. Downstream of the permeable obstacle, the changes to the flow conditions create more considerable changes in the evaporation rate. As shown in Fig. 10b, extensions to the length of the corner eddy reduce evaporation rates.
When these results are compared to the evaporation rates from the partially saturated obstacle of the same size from Sect. 4.2 (H = 10 cm), the evaporation rates within both corner eddies are higher for the dry obstacle case ( s w = 0 ). Along the surface of the dry obstacle, no evaporation will occur, maintaining a temperature and water vapor concentration closer to the main stream. In comparison with the partially saturated obstacle, the recirculated air will be warmer and dryer as it returns to the interface, increasing the evaporation rate. If a saturated obstacle were to dry out over time, the air recirculated to the surface would increase in evaporation potential, increasing the evaporation rate at the base of the obstacle (Fig. 11).

Obstacle Spacing: Two Partially Saturated Obstacles of Varied Spacing
To further analyze the effects of interfacial forms, a series of two obstacles is introduced. The two obstacles, hydraulically connected to the subsurface, of height and width H = 10 cm, and identical in flow parameters and saturation to those in Sect. 4.2, are added to the interface. The spacing, S, between these two obstacles is then varied: S = 25 cm, 50 cm, 75 cm, 100 cm, 125 cm, and 150 cm, as shown in Fig. 12. The free-flow domain is coupled with the porous media flow domain for 50 cm upstream of the first obstacle and for 50 cm downstream of the second obstacle for each test.
As given in Sects. 4.2 and 4.3, flow structures develop upstream and downstream of flow obstacles and their scope depends on the size and flow properties of the obstacle. In this example, two obstacles are placed near to each other, such that downstream and upstream flow structures within the spacing will overlap. This overlap of structure lengths will modify the form of the structures that develop between them, and the Fig. 9 The coefficient of friction C f (22) along the domain interface upstream of (a), and downstream of (b) a dry obstacle of permeability K Fig. 10 (a) Evaporation rates upstream of an obstacle of permeability K, plotted across the interface length at 4 days. (b) Evaporation rates downstream of an obstacle of permeability K, plotted across the interface length at 4 days. Also shown in both figures is the average upstream/downstream evaporation rate at 4 days against the obstacle permeability near-surface flow characteristics. This is shown in Fig. 13a, where the flow direction and magnitude are shown for two of the test cases. Additionally, Fig. 13b shows the total vertical exchange in and out of the spacing at the height of the obstacle. The following analysis of flow properties will be evaluated across the unit distance between the obstacles, or the distance, x, from the first obstacle divided by the spacing S.
To analyze the near-interface flow conditions, the C f [Eq. (22)], is plotted along the interface between the obstacles in Fig. 14a. The flow along the interface can be classified here into three sections from left to right: (1) a positive, low velocity gradient corner eddy, (2) a negative recirculation zone with a velocity gradient that decreases with increased spacing, and (3) a positive corner eddy with a velocity gradient that decreases with increased spacing. The proportional size of the second negative recirculation zone increases with increasing spacing, and the velocity gradients within this zone decrease with increasing spacing. The third, positive recirculation zone decreases in proportional length with increasing spacing, and for two spacing cases, 125 cm and 150 cm, does not occur. The velocity gradient in this recirculation zone also decreases with obstacle spacing. For the smallest spacing case, 25 cm, this last recirculation zone splits into a further negative flow zone directly adjacent to the second obstacle.
These flow conditions affect the evaporation rates across the coupled interface, as shown in Fig. 14b. As discussed in the previous sections, corner eddies, characterized by low velocity gradients and proximity to obstacles, are characteristic of lower evaporation rates, with local maximums dependent on their rotation and mixing with the mainstream. This can be seen in the first positive recirculation zone. Within the second zone, characterized  by higher velocity gradients and more mixing with the mainstream, the evaporation rate increases and reaches a maximum at the point of flow direction change. For the two larger spacing cases where the flow direction does not change, this local maximum occurs at the local minimum in C f , where the two negative recirculation structures merge. In the last positive corner eddy, the evaporation rate decreases toward the obstacle. Mixing with the mainstream can further be evaluated using the vertical velocities shown in Fig. 13b. Higher velocities in and out of the spacing between obstacles will correspond to more mixing to the mainstream, while larger obstacle spacings will correspond to larger mixing surfaces (Fig. 15).

Obstacle Shape: Two Partially Saturated Obstacles of Varied Shape
As the shape of the interface has a considerable effect on the resulting flow structures and transport, the interface's form can also alter the atmospheric demand at the soil-atmosphere Evaporation rates between two obstacles separated by a spacing S, plotted across the interface length at 4 days. Also shown is the average evaporation rate at 4 days plotted against obstacle spacing interface. To investigate this, five different obstacle shapes, or interfacial forms, are introduced: sinusoidal, rectangular, triangular, sawtooth, and reverse sawtooth (shown in Fig. 16a). To first evaluate the flow structures that develop from these forms, we consider these objects as isolated non-porous flow obstacles along the base of a free-flow domain.
In Fig. 15, the C f (22) is plotted along the base wall upstream and downstream of the obstacle. Upstream of the obstacle, as shown in Fig. 15a, the rectangular obstacle disrupts the flow the most, casting the largest flow structure upstream, followed by the reverse sawtooth shape where a similar upstream shape is cast. The triangular and sinusoidal obstacles, both similar in shape and profile, have similarly shaped upstream flow structures, where the smoother sinusoidal obstacle casts a slightly smaller structure in comparison with the sharper triangular form. The sawtooth shape, with its most gradual flow interruption, casts the smallest upstream flow structure. Downstream of the obstacle, the sinusoidal and the triangular forms both cast similarly shaped flow structures with the reattachment length for the sinusoidal form again slightly shorter due to its rounded shape. The rectangular shape casts the largest downstream flow structure. Downstream, the sawtooth form and the reverse sawtooth form switch roles, with the reverse sawtooth form casting the smallest flow structure and the sawtooth form a larger structure.
To test these shaped obstacles in a more complex coupled environment, we introduce a coupled interface of 1 m length. Upon this interface, there are two shaped obstacles, each 10 cm in height and length, spaced 20 cm away from each other. There are additionally 30cm of coupled interface upstream of the first obstacle and downstream of the second. A conceptual diagram of this test is shown in Fig. 16. Here, the shapes of the flow structures surrounding each obstacle will be affected by the shape of the obstacle, as well as the proximity to the other obstacle.
The flow structures cast upstream, downstream, and between these shaped obstacles all alter the atmospheric demand at the interface. To visualize the atmospheric demand and the impact on the evaporation rate, we look at Fig. 17b, where two interfacial forms are shown after 2 days of evaporation. In the free-flow domain, the water vapor content is shown, overlain with flow streamlines. In the porous media flow domain, the temperature is shown, where colder areas correspond to areas of higher evaporation and warmer areas to lower evaporation. For the rectangular obstacles, we can see that the low velocity gradient corner eddies in the free-flow domain directly adjacent to the obstacles will collect higher water vapor contents, reducing the atmospheric demand in this area. The upstream corner eddies will both mix more with mean flow quality air, creating higher evaporation rates in comparison with the downstream eddies. The effect this has on the evaporation rate can be seen in the porous medium with temperature. Areas in the porous medium that are adjacent to collection zones are warmer, corresponding to reduced evaporation. Areas adjacent to the recirculation zones where higher velocity gradients and more mixing with the main flow correspond to colder porous media flows.
For the reverse sawtooth obstacles, only two corner eddies develop, upstream of each obstacle, where the eddy upstream of the second obstacle is reduced in scope due to its distance from the first obstacle. Downstream of the second obstacle, a recirculation zone does not develop, and no corner eddy develops. The two corner eddies that do develop upstream of the obstacles collect cooler water vapor, reducing the atmospheric demand, but still mix with main stream quality air, raising the evaporation rate. Looking at the porous media, only the area adjacent to the first obstacle shows a higher temperature, meaning evaporation is largely enhanced across this surface.
Dry out can also be observed within the obstacles heights where higher temperatures develop. When the saturation in these peaks is completely reduced, air with main flow properties invades the pore space, and as no evaporation will occur, leading to higher temperatures. In these locations, phase II evaporation, where these higher temperatures are standard (Haghighi and Or 2015a), has already begun.
To compare the exchange rates across these formed surfaces, the average evaporation rate across the full coupled surface is calculated per time step and plotted over time, as shown in Fig. 17a. As first mentioned in Sect. 4.1, the evaporation rates seen here are a mixture of various phase I evaporation rates, depending on the flow conditions at the surface, as well as some phase II evaporation rates in locations where dry out has begun at the peaks of the obstacles.
The flow structures developing about the rectangular obstacles reduce the atmospheric demand at the surface the most in comparison with the other interfaces. Following this is the sawtooth form, that casts downstream flow conditions similar to those of the rectangular domain, but develops smaller upstream flow structures. The reverse sawtooth forms maintain the highest atmospheric demand along the interface, as the only low velocity flow structures they cast are the smaller corner eddies upstream of the obstacles. The sinusoidal

Conclusions
A review: The flow dynamics in a free-flow domain play a large role in the exchange processes across the soil-atmosphere interface. Interfaces that include non-flat geometries, such as obstacles or interfacial forms, create flow structures adjacent to the domain interface that alter the transport of mass and energy in the free flow. These flow structures modify the shape of the mass, momentum, and energy sublayers at the interface, and therefore the exchange rates across the interface.
After investigating the various evaporation rates that can occur across a flat surface, the model concept and implementation outlined in Sects. 2 and 3 are used to investigate the coupled system under different geometries. First, single unsaturated obstacles of varied height are investigated. Dependent on their height, these obstacles cast different flow structures upstream and downstream. In locations where the velocity gradient decreases, such as corner eddies, the evaporation rates will decrease. If these corner eddies provide a surface where mixing with the main flow can occur, evaporation rates will remain higher, whereas eddies that only mix with other recirculation zones will decrease further due to the collection of water vapor within their scope. Within the recirculation zones with higher velocity gradients, evaporation rates will remain high, and the turbulence enhanced mixing with the main flow, along the eddy's surface, will further enhance the evaporation rate.
The flow structures that develop around a dry highly permeable obstacle will also change their form given different obstacle flow parameters. Higher obstacle permeabilities will thin the flow structures upstream of the obstacle and extrude those downstream, enhancing evaporation upstream and reducing it downstream. Further, when comparing these dry obstacles with partially saturated obstacles, the evaporation rates within the Fig. 17 (a) The average evaporation rate plotted over time due with respect to changes in obstacle shape. (b) Atmospheric conditions around two different interfacial form configurations. In the free-flow domain, the water vapor mass fraction is shown overlain with flow streamlines. In the porous media flow domain, the temperature (K) throughout the domain is shown corner eddies adjacent to the obstacles are higher, as temperatures will remain higher, and no evaporated water vapor from the surface of the obstacle will collect in the eddies.
In a series of two partially saturated obstacles, the spacing between the obstacles will change the evaporation rates across the interface that separates them. The flow structures that develop between these obstacles will merge and change the flow dynamics depending on their spacing. Depending on the velocity gradients along the surface separating the obstacles and the overall mixing with the flow above this spacing, evaporation rates will be either reduced or enhanced.
The shape of these obstacles will also interrupt the flow in different manners. Blunt sharp obstacles will reduce the evaporation rate adjacent to them the most, where more gradually sloping obstacles will produce smaller flow structures, maintaining a higher evaporation rate. As the flow structures that develop downstream of the obstacles are larger, obstacles with more gradual slopes on their downstream side will reduce the evaporation rate the least.
These findings are consistent throughout each investigation; flow structures that create higher velocity gradients and greater mixing with the main flow provide more atmospheric demand at the interface, thus enhancing the evaporation rate. Flow structures with lower velocity gradients that circulate and collect evaporated water vapor at the surface without mixing with the main flow will reduce the atmospheric demand and the exchange rates across the interface.
An outlook: In the future, an expansion of the current model to include solar radiation effects to further develop our knowledge of larger scale energy and water balances would be of interest. For flow around obstacles and interfacial forms, the effects of shading and direct radiation would affect the evaporation rate around these forms, as well as contribute to the flow dynamics via thermal buoyancy and altered flow parameters. In addition to this, enhancements to the porous media model would also allow for the simulation of related environmental problems such as subsurface gas migration or evaporation driven salt precipitation, as an expansion to the work described in Shokri-Kuehni et al. (2020). Modifications to this model including multi-component gas transport or mineralization could be performed to allow for this. Further, although this model builds on model concepts developed on past experimental work, a further detailed experimental evaluation done on the same scale as this work, evaluating the formation of flow structures and local evaporation rates as presented in this article, would be beneficial. In addition, experimental work on a larger scales could help transfer the information learned in this work to reduced climatic heat and water balance models.