We theoretically investigated the effect of vapor flow on the drying front that develops in soils when water evaporates from the soil surface and on GPR data. The results suggest the integration of the full-wave GPR model with a coupled water, vapor, and heat flow model to accurately estimate the soil hydraulic properties.

## Abstract

We investigated the effects of a drying front that emerges below an evaporating soil surface on the far-field ground-penetrating radar (GPR) data. First, we performed an analysis of the width of the drying front in soils with 12 different textures by using an analytical model. Then, we numerically simulated vertical soil moisture profiles that develop during evaporation for the soil textures. We performed the simulations using a Richards flow model that considers only liquid water flow and a model that considers coupled water, vapor, and heat flows. The GPR signals were then generated from the simulated soil water content profiles taking into account the frequency dependency of apparent electrical conductivity and dielectric permittivity. The analytical approach indicated that the width of the drying front at the end of Stage I of the evaporation was larger in silty soils than in other soil textures and smaller in sandy soils. We also demonstrated that the analytical estimate of the width of the drying front can be considered as a proxy for the impact that a drying front could have on far-field GPR data. The numerical simulations led to the conclusion that vapor transport in soil resulted in S-shaped soil moisture profiles, which clearly influenced the GPR data. As a result, vapor flow needs to be considered when GPR data are interpreted in a coupled inversion approach. Moreover, the impact of vapor flow on the GPR data was larger for silty than for sandy soils. These effects on the GPR data provide promising perspectives regarding the use of radars for evaporation monitoring.

**Soil evaporation is an important land surface process** because it is closely linked to mass and energy exchanges between the land surface and the atmosphere. Soil evaporation rates generally proceed in two different main phases. During the first stage, evaporation is only affected by the amount of energy available to vaporize soil water in the upper layer of the soil given that mass transfer limitations in the air layer above the drying soil surface do not become limiting (Shahraeeni and Or, 2012). During the second phase, transfer of water and vapor in the soil toward the soil surface becomes limiting so that the evaporation rate is determined by the soil hydraulic properties. Dynamic variations of soil moisture in near surface soil layers are key factors to distinguish different stages of evaporation. During Stage I evaporation, the soil surface water content decreases with time and the gradient of soil moisture at the soil surface increases over time, whereas during Stage II, surface soil moisture remains more or less constant but the soil moisture gradient near the surface decreases over time and a drying front develops that penetrates into the soil. In addition, to associate soil surface states to subsurface states and properties, it is important to have knowledge about the dynamic variations of the subsurface soil water content profile.

The one-dimensional Richards equation is often used to simulate vertical soil moisture profiles (Richards, 1931). However, it only considers isothermal liquid water flow, and because of the shape of the hydraulic soil properties, Richards’ equation cannot reproduce S-shaped soil moisture profiles or drying fronts that typically develop during evaporation (Vankeulen and Hillel, 1974). In this respect, Saito et al. (2006) showed that the dynamics of soil water content profiles are highly influenced by vapor transport. As a consequence, it is essential to consider coupled liquid water, water vapor, and heat transport (coupled model) to more accurately estimate soil water content profiles and heat dynamics in field soils. In particular, the coupled model plays an important role in predicting near surface S-shaped soil moisture profiles that develop during the second stage of soil evaporation. In addition, this model also has many agricultural applications (Wuest et al., 1999), engineering applications like landfill management (Khire et al., 1997) and nuclear waste repositories (Spycher et al., 2003), and military applications such as improved detection of buried land mines (Šimùnek et al., 2001; Lopera et al., 2007).

The derivation of soil moisture profiles at the field scale from local scale measurements would require interpolation and may overlook variations that are at a smaller scale than the distance between the local soil sensors. Resorting to geophysical methods such as ground-penetrating radar (GPR) that use electromagnetic waves with wavelengths that are smaller than the scales over which soil moisture changes close to the soil surface is therefore promising to resolve drying fronts. Because of wave reflections produced by electromagnetic contrasts, GPR allows for subsurface imaging with a high spatial resolution. Time-lapse GPR has been proven to be a promising method to retrieve soil physical and hydraulic properties in the vadose zone (Hubbard et al., 1997; Hubbard and Rubin, 2000; Lambot et al., 2004a; Moysey, 2010; Moghadas et al., 2010; Andre et al., 2012). Lambot et al. (2004b) proposed an accurate and computationally effective off-ground GPR forward model. The interest in using off-ground GPR stems from its ability to perform measurements in the far-field region of the antenna, resulting in a relatively simple but accurate intrinsic representation of the antenna (Lambot et al., 2004b, 2006a). From a practical point of view, the off-ground antenna minimizes disturbance of the soil surface and its exposure to natural boundary conditions (precipitation, evaporation) (Serbin and Or, 2003, 2004).

In the context of developing integrated hydrogeophysical inverse modeling techniques, Lambot et al. (2006a) integrated the radar model of Lambot et al. (2004b) with the Richards equation, resulting in a coupled electromagnetic–hydrodynamic inverse modeling scheme for the remote estimation of the shallow soil hydraulic properties in unsaturated soils. Jadoon et al. (2008) theoretically investigated the uniqueness and stability of this integrated inversion approach to retrieve soil hydraulic parameters for different soil types. Lambot et al. (2006a, 2009) and Jadoon et al. (2008) showed that enough information may be contained in the off-ground GPR data when it is combined with hydrological models to retrieve time-lapse, near surface soil moisture profiles and hydraulic properties. However, these hydrogeophysical studies were performed mainly for infiltration and drainage conditions but did not consider extended evaporation phases. Moreover, as we already discussed, the dynamics of the soil water content during extended evaporation phases cannot be entirely reconstructed by only considering liquid water flow and neglecting vapor transport. As a result, considering liquid water, water vapor, and heat transport in a coupled electromagnetic–hydrodynamic model is expected to provide additional information about the soil hydraulic properties, particularly regarding the near surface soil evaporation processes.

We evaluated the capability of different soil textures to generate evaporation drying fronts by analytically simulating the soil water content profile at the beginning of the second phase of evaporation. In this respect, the width of the simulated soil moisture profiles (analytical parameter) was considered as a quantitative criterion in our further analysis. First, a sensitivity analysis of the relation between this analytical parameter and soil hydraulic properties was performed. Subsequently, the sensitivity of the GPR data to drying fronts that develop during Stage II of soil evaporation was investigated for different soil textures and correlated with the analytical estimate of the width of the drying front. Therefore, we numerically simulated the vertical soil moisture profiles using HYDRUS-1D (Šimùnek et al., 2008, 2009; Séré et al., 2012; Reading et al., 2012), which numerically solves the Richards equation, and using a modified version of HYDRUS-1D (coupled model) that solves the coupled water, vapor, and heat flow equation (Saito et al., 2006). Then, we modeled the GPR data using a full-waveform GPR forward model (Lambot et al., 2004b, 2006a), taking into account the frequency dependency of dielectric permittivity and apparent electrical conductivity in the GPR operating frequency ranges (800–5200 MHz).

## Materials and Methods

### Richards Equation

*t*is time (d), ψ is the pressure head (cm),

*Z*is the elevation head (cm), and

*K*(ψ) is the soil hydraulic conductivity (cm d

^{−1}) described by van Genuchten (1980) as

^{−1}),

*n*, and

*m =*1 − (1/

*n*) are shape parameters, τ is a tortuosity factor, and

*K*

_{s}(cm d

^{−1}) is the saturated hydraulic conductivity. The volumetric soil water content values are calculated by

_{s}and θ

_{r}are saturated and residual soil water content, respectively, and is the effective saturation.

### Coupled Water, Vapor, and Heat Transport Model

*q*

_{h}is the heat flux density (J cm

^{−2}d

^{−1}), λ is the thermal conductivity (J cm

^{−1}d

^{−1}°C

^{−1}),

*T*is temperature (°C),

*L*

_{0}is the latent heat of vaporization (J g

^{−1}), and

*q*

_{v}is vapor flow (g cm

^{−2}d

^{−1}) described by (Saito et al., 2006)

*D*

_{v}is the vapor diffusivity in soil (cm

^{2}d

^{−1}), is the saturation vapor concentration (g cm

^{−3}),

*s*is the slope of the saturation vapor concentration versus temperature function,

*z*is the spatial coordinate positive upward (cm),

*h*= exp(

*M*

_{w}

*g*ψ/

*RT*) is the fractional relative humidity in which

*M*

_{w}is the molecular weight of water (= 18.015 g mol

^{−1}),

*g*is the gravitational acceleration (= 980 cm s

^{−2}), and

*R*is the universal gas constant (= 8.31 J mol

^{−1}°C

^{−1}).

Surface precipitation, irrigation, and heat fluxes, including the latent heat flux, are used as boundary conditions for the coupled model. The heat fluxes at the soil surface are obtained from meteorological data (radiation, air temperature, air humidity, and wind speed) by iteratively solving the soil surface energy balance in combination with the coupled water, vapor, and heat transport model. Including vapor transport (Eq. [5]) in the water balance equation (Eq. [1]) and solving it together with the heat flow equation (Eq. [4]) and the surface energy balance provides a very flexible way of linking meteorological information with the soil energy and water balance. The reader is referred to Saito et al. (2006) for additional details on the liquid water, water vapor, and heat transport model.

### Hydrological Simulations

Water content profiles that develop under laboratory conditions in 30-cm thick soil layers were simulated using HYDRUS-1D and the coupled flow model for 12 different soil textures. In the coupled model, the energy balance at the soil surface was solved using a defined incoming longwave radiation, a wind speed of 350 km d^{−1}, air humidity of 50%, and air temperature of 25°C at 2 m above the soil surface to define the upper boundary conditions of water, vapor, and heat fluxes in the soil (Saito et al., 2006). Since we considered an indoor laboratory experiment without shortwave radiation, we calculated the incoming longwave radiation as the radiation emitted by a black body at the temperature of the ambient air multiplied by the short wave emissivity of the soil, which was assumed to be 0.9. The emitted shortwave radiation was calculated from the soil surface temperature, which was due to evaporative cooling lower than the ambient air temperature, and by using the shortwave emissivity. As initial conditions, a hydrostatic equilibrium with ψ equal to 0 at *z* = −30 cm and a uniform soil temperature profile were used. The coupled model simulates the evaporation rate at the soil surface for these conditions. To define the bottom boundary, we assumed a no flow boundary condition for the water, a vapor flow, and a constant temperature that was equal to the ambient temperature in the laboratory. These bottom boundary conditions are representative for a laboratory experiment. For field conditions, the shortwave radiation and the heat and water exchanges with the subsoil also need to be considered.

After a short initialization period, during which the soil surface temperature and the evaporation rate decreased until an equilibrium was reached, the simulated evaporation rate reached a constant value (Phase I evaporation rate) until the soil surface dried out and this evaporation rate could not be sustained by water flow in the soil and the evaporation regime switched to Phase II of evaporation. The simulated evaporation during Phase I by the coupled model was used as flux boundary condition for Richards’ equation until a critical pressure head at the soil surface was reached (−10^{5} cm), after which the pressure head at the soil surface was kept constant and equal to −10^{5} cm.

### Analytical Approach

The magnitude or extent of the drying fronts that develops during the second stage of evaporation depends on soil texture. To assess the extent of a drying front that could develop in a soil, the common approach is to perform simulations with hydrological models. However, hydrological simulations are computationally complex and do not provide a closed-form or direct relation between soil hydraulic parameters and a process property of interest, that is, in our case the width of the drying front. Therefore, we introduced a simple approximate analytical approach to derive the magnitude of the drying front that develops at the beginning of the second stage of the evaporation. We define the magnitude of the drying front as a difference between the soil moisture at a certain depth in the soil, where the evaporation process does not influence the soil water content substantially, and the soil moisture at the soil surface. This method makes it possible to evaluate the capability of different soil textures to develop drying fronts without implementing complex and time consuming numerical simulations. In addition, the effect of different hydraulic parameters on the drying fronts can be simply investigated by this approach.

*L*(cm) from which water is uniformly extracted by evaporation. The water flux

*J*

_{w}at a certain depth

*z*> −

*L*(

*z*= 0 at the soil surface) is hence

*E*is the evaporation rate. Neglecting gravity, the water flux at

*z*is

^{5}cm at the soil surface, that is, at

*z*= 0, and integrate left- and right-hand sides, we obtain

*z*. Using the van Genuchten (1980) relation for θ(ψ), the soil water content profile can be derived. This profile is an approximation of the profile that is simulated by solving Richards’ equation at the end of Phase I. The proposed analytical method provides a simple and rapid approach to estimate the soil water content profile at the end of Stage I of evaporation using evaporation rate, soil layer thickness, and van Genuchten–Mualem parameters as input. The difference between the soil water content at the surface and at the bottom of the profile or the width of the soil water content profile, Δθ

^{anl}(see Fig. 1), is considered as a criterion to assess the extent of a drying front that could develop in a soil. We expected that the analytical parameter Δθ

^{anl}can be considered as a proxy for the impact that a drying front could have on the radar data.

### GPR Forward Modeling

*S*

_{11}is the quantity measured by the VNA,

*b*(ω) and

*a*(ω) are, respectively, the backscattered and incident waves at the VNA reference calibration plane, ω is the angular frequency,

*H*

_{i}(ω) is the antenna return loss,

*H*(ω) is the transmitting–receiving transfer function accounting for antenna gain and phase delay,

*H*

_{f}(ω) is the feedback loss accounting for the interactions occurring between the antenna and the ground, and is the spatial domain Green’s function of the air–subsurface system modeled as a planar three-dimensional multilayered medium formulated as (Lambot et al., 2004b, 2006a, 2007)

In this expression, the subscripts denote layer indexes, *R*^{TM} and *R*^{TE} are, respectively, the transverse magnetic (TM) and transverse electric (TE) global reflection coefficients accounting for all reflections and multiples from inferior interfaces, *k*ρ is the spectral domain counterpart of ρ and ρ(=0) is the distance between the source and the receiver, Γ is the vertical wavenumber defined as , while with the magnetic permeability μ, dielectric permittivity ɛ, and electrical conductivity σ, *j* is equal to , , and . For the free-space layer 1, we have with *c* being the free-space wave velocity. We simulated GPR data from estimated soil water content profiles for the 800 to 5200 MHz frequency range (6 MHz frequency step) with 27-cm antenna height above the soil surface (far-field condition). The antenna transfer functions used in the radar model inherently account for the frequency-dependent phase center position of the antenna (Jadoon et al., 2011).

### Frequency Dependency

_{e}) is calculated as

_{w}, ɛ

_{a}, and ɛ

_{s}are the permittivities of the water, air, and soil particles, respectively, and γ is the empirical coefficient representing for the soil structure and its orientation with respect to the electromagnetic field. We considered γ = 0.5 for this equation, which is widely known as Complex Refractive Index Model (Shutko and Reutov, 1982). The permittivity of water ɛ

_{w}is a function of the frequency and expressed by Debye’s equation (Debye, 1929):

*ɛ*

_{∞}= 4.9 is the dielectric constant at infinite frequency,

*ɛ*

_{0}= 8.854 × 10

^{−12}is the permittivity of the free space,

*ɛ*

_{sp}= 88.045 − 0.4147

*T*+ 6.295 × 10

^{−4}

*T*

^{2}+ 1.075 × 10

^{−5}

*T*

^{3}is the static dielectric constant (Klein and Swift, 1977), τ

_{0}= (1/2π) (1.1109 × 10

^{−10}− 3.824 × 10

^{−12}

*T*+ 6.938 × 10

^{−14}

*T*

^{2}− 5.096 × 10

^{−16}

*T*

^{8}) (Stogryn, 1971) is the relaxation time, and σ

_{DC}(S m

^{−1}) is the ionic conductivity of water.

It is interesting to note that the application and implementation of the power law model (Eq. [13]) is subject to a priori knowledge on soil porosity and permittivity of the mineral matrix (*ɛ*_{s}). However, no direct measurement technique has been introduced for *ɛ*_{s}, and the estimated value of 5 is commonly used for soil particle permittivity (Robinson and Friedman, 2003). In this paper, we assumed that the porosity is equal to the saturated soil water content and that *ɛ*_{s} is equal to 5 for the different soil textures. The uncertainties in soil matrix permittivity are expected to introduce minor errors that are not considered in our numerical simulations.

## Numerical Simulations and Results

### Analytical Method

*L*= 30 cm vertical thickness that is exposed to an evaporation rate of

*E*= 0.2 cm d

^{−1}. Table 1 shows the van Genuchten–Mualem parameters and the corresponding analytical results for different soil textural groups. As can be clearly seen, the sandy soil provides the lowest values of Δθ

^{anl}(= 0.09), while silty soil and clay loam soil exhibit a higher potential to generate wide drying fronts (Δθ

^{anl}= 0.3). For all soil textural groups, water content profiles were simulated in a 30-cm-thick soil layer using HYDRUS-1D and the coupled flow models. The evaporation periods of sand, loamy sand, and sandy loam were 34, 34, and 47 d, respectively (the simulations were stopped as the simulated water contents were close to the residual water contents), and for the other soil textures, in which the Phase I evaporation period lasted longer, evaporation periods of 80 d were considered. The GPR signals were simulated from estimated soil water content profiles, and subsequently, the RMS errors between simulated GPR data for each soil texture were calculated using the following equation:

*N*is the total number of Green’s function elements. Table 1 shows the RMSE values for all soil textures. A plot of RMSE values versus Δθ

^{anl}(Fig. 2) demonstrates that the deviation between GPR signals derived from moisture profiles simulated with HYDRUS-1D and with the coupled model increases with increasing Δθ

^{anl}or increasing extent of the drying front. As a consequence, the analytical parameter Δθ

^{anl}can be considered as a proxy for the impact that vapor flow could have on soil moisture profiles in drying fronts and consequently on radar data. For Δθ

^{anl}larger than 0.15, huge RMSE are observed in Fig. 2. Consequently, the analytical parameter larger than 0.15 is considered as a quantitative criterion for a significant impact of vapor flow and soil moisture profiles in drying fronts and on GPR data.

### Sensitivity Analysis

In addition to the relation between soil texture and Δθ^{anl}, we investigated the relation between Δθ^{anl} and the van Genuchten–Mualem parameters. In general, θ_{r} and θ_{s} do not vary over large ranges between different soil textures. The parameter θ_{r} is regarded as an empirical parameter and can be fixed either to a value which yields the best fit to the experimental water retention data or to zero. Furthermore, θ_{s} and θ_{r} can be directly derived from GPR measurements, respectively, under saturated and dry conditions at the soil surface (Lambot et al., 2009; Jadoon et al., 2012). The parameter τ is generally considered to have a small effect on the hydrodynamic events and can be either fixed to the average value of 0.5 (Mualem, 1976) or, preferably, estimated from bulk density and hydraulic conductivity with which it is highly correlated (Vereecken, 1995). Given these considerations, we assumed θ_{r}, θ_{s}, and τ to be known in the sensitivity analysis and the response function Δθ^{anl} was calculated considering three parameter pairs α–*n*, α–*K*_{s}, and *n*–*K*_{s}. For all calculations, the parameter ranges were relatively large so that they represent the prevailing soil hydraulic properties of most soils, namely, 0.001 < α < 0.2 cm^{−1}, 1.02 < *n* < 3, and 5 < *K*_{s} < 1000 cm d^{−1}. The range of each parameter was divided into 50 discrete values, resulting in 2500 calculations for each contour plot. The sensitivity analysis was performed for a soil column with 30-cm thickness and evaporation rate of 0.2 cm d^{−1}. The response surfaces of the analytical parameter are presented in Fig. 3.

The topography of the response function provides valuable insights into the sensitivity of the Δθ^{anl} to the different hydraulic parameters and parameter correlations. Regarding the α–*n* plane (Fig. 3a), for 1.04 < *n* < 1.7 and α > 0.015 cm^{−1}, the response function presents values higher than our threshold (= 0.15). In these ranges of parameters, α is positively correlated with Δθ^{anl}. Concerning the α–*K*_{s} pair (Fig. 3b), *K*_{s} is negatively correlated with Δθ^{anl}, and for *K*_{s} < 200 cm d^{−1} and α > 0.001 cm^{−1}, the response function has values higher than our threshold. For the *n*–*K*_{s} pair (Fig. 3c), the response function presents values higher than 0.15 for *K*_{s} < 200 cm d^{−1} and 1.04 < *n* < 3. This analysis demonstrates that, theoretically, soils with *n* < 1.7, α > 0.015 cm^{−1}, and *K*_{s} < 200 cm d^{−1} present higher capability than other soil textures to develop significantly large drying fronts (Δθ^{anl} > 0.15) during the second phase of evaporation and consequently present more eminent impacts of vapor flow on GPR data. These ranges of parameters are mainly attributed to medium textured soils that are characterized by a wide pore size distribution (Leij et al., 1996). Given the van Genuchten–Mualem hydraulic parameters presented in Table 1, this analysis shows that soil textures such as loam, silt, silt loam, sandy clay loam, clay loam, and sandy clay present more considerable influences of vapor flow on GPR data than other soils, which is consistent with the results presented in Fig. 2.

### Water Content Profiles

In the following, we discuss the results of numerical simulations for sandy and silty textures in detail as these two soil textures provide the lowest and the highest values of Δθ^{anl}, respectively. The simulated evaporation rates during Phase I amounted to 0.22 cm d^{−1} for the sandy soil and 0.25 cm d^{−1} for the silty soil (the higher evaporation rate in the silty soil is a result of the higher heat conductivity of the silty soil due to its higher soil water content and the constant temperature at the bottom of the profile). Figure 4 shows the soil water content profiles for the whole simulated evaporation period. The soil water content profiles derived using the analytical approach are also shown by dashed lines. We calculated the width of the soil water content profiles for the transition time for HYDRUS-1D (Δθ^{HD}) and also for coupled model (Δθ^{CM}) data to compare with the analytical results. The soil moisture profiles simulated by the coupled model are shown in Fig. 4c and 4d. We replotted these two subplots for smaller depths (Fig. 4e and 4f) to better visualize the drying fronts during the second phase of evaporation.

For the silty soil, the transition between Phase I and Phase II occurs at Day 17, and the analytical result is almost in agreement with the simulations (Δθ^{HD} ≈ Δθ^{CM} ≈ Δθ^{anl} ≈ 0.3); that is, the analytical model estimates the width of soil water content profile at the end of Phase I of the evaporation, close to those simulated by the hydrological models. For the sandy soil, the Δθ^{HD} and Δθ^{CM} at Day 7 (transition time) show almost the same values (≈0.3), which is completely different from the analytical results (Δθ^{anl} = 0.09). This difference between numerically simulated and analytically derived width of the soil moisture profile in the sandy soil may be attributed to the hydrostatic initial conditions, which still have an important impact on the moisture distribution deeper in the soil profile. However, when looking at the soil moisture profiles closer to the soil surface, which is of relevance for the development of drying fronts and also is the region of interest for GPR, a better agreement between the numerical and analytical profiles can be observed. This indicates that the analytical method may still be useful for providing a proxy of the drying front width.

### GPR Forward Modeling Results

We simulated GPR signals using soil water content profiles presented in Fig. 4. The GPR data are represented by the Green’s function in both frequency and time domains. Figure 5 shows the phase of Green’s functions in the frequency domain. For better illustration and comparison, we plotted the phase data in frequency ranges between 4000 and 5000 MHz for five different evaporation times including the transition days. The drying fronts that developed during the second phase of evaporation give rise to small phase shifts in GPR data, while significant changes can be observed in the amplitude of Green’s functions in the frequency domain (Fig. 6). Phase changes are small as the layer thicknesses are small compared to the wavelengths. For the coupled model, the simulated Green’s function amplitudes at the higher frequencies do not change monotonically with time whereas for Richards’ equation, they decrease monotonically over time. The Richards equation leads to moisture (dielectric permittivity) gradients at the soil surface that decrease monotonically with time and therefore leads to a decrease of the reflection coefficient with time. The soil water content gradients at the soil surface produced by the Richards equation are larger than those simulated by the coupled model. This explains the larger Green’s function amplitudes corresponding to larger reflection coefficients for the Richards equation. The development of drying fronts first leads to a decrease of the soil water content gradients at the soil surface, and concurrently a decrease in GPR signal amplitude. During the establishment of the drying front, the soil moisture gradient at the drying front below the soil surface increases which explains the increase of the GPR signal amplitude over time again.

Figure 7 shows the Green’s function in time domain for five different evaporation times including the transition days. In time domain, the amplitude of the Green’s function decreases monotonically during the evaporation period for the soil moisture profiles simulated by the Richards equation. For the coupled model simulations, the amplitude of the Green’s function initially decreases, and after the establishment of the drying front, it increases again. When there is a constant soil water content profile (no layering), the amplitude of the reflection depends directly on the electromagnetic contrast between the air and the soil surface. In the presence of layering, the reflections at the different interfaces may interfere in a constructive or destructive way, thereby causing a larger or smaller reflection amplitude (Kowalsky et al., 2001; Lunt et al., 2005; Lambot et al., 2006b). During the evaporation process, the permittivity at the soil surface continuously decreases and therefore the surface reflection continuously decreases as well. However, the amplitude of the Green’s function increases again for the coupled model simulations that originate from a second reflector. In that case, either the decrease comes from a transition in the profile that becomes less sharp in a first stage and then sharper in a second stage (thereby producing a larger reflection), or the transition in the profile is quite sharp in both cases, leading to significant reflection; however, depending on the depth of the transition, it leads to constructive or destructive interferences with respect to the surface reflection (depending on the wavelength and transition depth) and therefore to a larger or smaller amplitude of the Green’s function, respectively. This effect is more prominent for the silty soil (Fig. 7b and 7d) than for the sandy soil (Fig. 7a and 7c).

## Conclusions

We theoretically investigated the effect of vapor flow on the drying front that develops in soils when water evaporates from the soil surface and on GPR data. Interpretation of the GPR data with a model that does not include vapor flow or with a simpler model that does not resolve soil moisture profile and provides an average of depth profile soil water content in a sensitive soil layer hampers the accurate retrieval of the soil hydrogeophysical properties.

The results suggest that the description of the liquid water, water vapor, and heat transport is important for GPR modeling, since the GPR method is sensitive to the soil moisture profile close to the soil surface. This sensitivity not only leads to a change in signal strength but also to a fundamentally different temporal evolution of the GPR data. The numerical results presented here suggest the integration of the full-wave GPR model with a coupled water, vapor, and heat flow model to accurately estimate the soil hydraulic properties. Moreover, the analytical approach introduced in this paper permits a simple investigation of the capacity of different soil types to develop a drying front and provides a proxy for the impact that an S-shaped soil moisture profile could have on GPR data. Our results indicate that the largest effects can be expected in medium textured soils that are characterized by a wide pore size distribution. We believe that the proposed integrated approach constitutes a valuable alternative to the commonly used methods and represents the basis of a new integrated tool with great promise for in situ characterizations of the shallow soil hydraulic properties with a high spatial resolution. Moreover, since the GPR method allows for electrical conductivity measurement, the proposed approach is promising for solute transport modeling and solute accumulation at evaporating surfaces. In addition, this study can be further extended to cases where remote sensing data with active radars have to be interpreted in an integrated hydrogeophysical modeling scheme. Our intention is to focus next on the validation of the proposed approach under laboratory and field conditions.

D. Moghadas is funded by the research unit Multi-scale Interfaces in Unsaturated Soil (MUSIS) of the DFG FOR 1083 (Deutsche Forschungsgemeinschaft).