Influence of soil heterogeneity on evapotranspiration under shallow water table conditions: transient, stochastic simulations

Ensembles of soil column numerical experiments were performed to study the influence of heterogeneity in the saturated hydraulic conductivity on the evapotranspiration or latent heat, LE, under varying water table conditions. In the numerical experiments, a variably saturated groundwater flow model coupled with a land surface model (ParFlow[CLM]) was used. The model was forced at the top with an atmospheric time series from Oklahoma, USA. The heterogeneity was simulated for two different soils (clay and sand) using uncorrelated Gaussian random fields with different variances. Soil heterogeneity has a strong influence on LE during the dry months of the year and negligible influence during months with sufficient moisture availability. The influence is stronger for unstructured soils, such as sand. An increase in shallow water table depth collapses the ensemble onto a single curve, i.e., heterogeneity plays a minor role and LE is limited by the increasing redistribution distance from the water table to the land surface. Accounting for correlations between the hydraulic conductivity and the shape factor in the pressure–saturation function increases the variability in LE. In the case of laterally interconnected soil columns, 3D hydrodynamics homogenize LE fluxes by establishing additional flow paths toward the land surface. Comparison with geometric mean simulations shows good agreement over many time periods for weekly and monthly averaged LE fluxes. Instantaneous daily and hourly differences sometimes exhibit values on the order of 100–101 W m−2 (10−2–10−1 mm d−1) and white noise behavior over extended time periods. In the simulations, perhaps the strongest limitation is the application of a constant atmospheric time series at the top of the simulation domain. In order to assess the impact of this limitation in simulations of subsurface–land surface–atmosphere interactions, it would be necessary to allow for two-way feedbacks with the lower atmosphere. This will require resolving the boundary layer immediately adjacent to the land surface, including the lateral exchange of mass, energy, and momentum.


Introduction
Soil heterogeneity strongly influences the transport of moisture and solutes in the shallow subsurface (e.g., Harter and Zhang 1999, Maxwell et al 2007, Onsoy et al 2005, Yeh et al 2005. Since the ultimate goal is the prediction of flow and transport processes across various spatial and temporal scales using theoretical models, the influence of heterogeneity must be accounted for in the calculations. This may pose a considerable computational problem especially for large spatial scales, because in order to capture accurately the transport behavior, the heterogeneity must be resolved at all scales, also at high spatial resolution. Additionally, a very limited amount of spatially distributed information may be available on the hydraulic and transport properties. Thus, sophisticated interpolation schemes are required to fill in the information gaps between points of high information content and quantify parameter uncertainty because of original data scarcity (Rubin 2003).
In order to remedy the problem of resolving the heterogeneity across various scales, upscaling of the flow and transport equations and the required parameters has been proposed that focuses on capturing the average behavior of the system under investigation. An excellent review on the history and state-of-the-art of upscaling of soil physical processes was provided by Vereecken et al (2007). Currently, upscaling laws are available only for relatively simple cases including bounded and unbounded domains. A major problem in deriving relationships of upscaling still constitutes highly transient processes, because e.g., constant upscaled parameters may not exist across different timescales and varying flow geometries (Rubin 2003).
An example of a highly transient process is evapotranspiration or latent heat, LE, which is a combination of evaporation from the bare soil and the transpiration of moisture by plants. LE is a key component of the mass and energy balance at the land surface with a distinct daily, monthly and yearly cyclicity overlain by high frequency variability due to e.g. variable cloud cover during the day. It has been suggested that LE strongly depends on the moisture and energy state of the shallow subsurface (e.g., Baldocchi and Xu 2007, Chen et al 2008, Kollet et al 2009, Siqueira et al 2009, which, in turn, is influenced by physical transport processes and thus soil heterogeneity. Since LE is part of the mass and energy balances at the land surface and influences weather generating processes it has been parameterized in land surface models (LSMs) that are applied commonly as lower boundary conditions in atmospheric models (and more recently as upper boundary conditions in subsurface flow and transport models). Again the difficulty arises in representing subgrid scale heterogeneity in LE due to soil heterogeneity, because land surface models are usually applied at very large scales (>10 2 m) from the soil physics perspective. (As a matter of fact, there is no lateral scale explicitly defined in LSMs, because transfer of mass, energy and momentum only occur vertically in laterally isolated, 1D columns. Some lateral scale is only introduced implicitly through the application of the e.g., Monin-Obukhov similarity concept that requires homogeneous land surface conditions at scales >10 2 m (Foken 2006)). In the land surface and atmospheric sciences community, characterization and parameterization of the influence of heterogeneity was mainly limited to lateral variability of soil properties and soil moisture, neglecting vertical heterogeneity of the shallow soil profile (e.g., Intsiful and Kunstmann 2008). It has been shown that this heterogeneity exerts considerable influence on the different components of the land surface mass and energy balances and attempts were made to account for this heterogeneity using stochastic and numerical techniques (Avissar 1991, Entekhabi and Eagleson 1991, Entekhabi and Rodrigueziturbe 1994, Famiglietti and Wood 1995. In these studies, however, vertical heterogeneity of the soil profile has been neglected, which will influence vertical moisture transfer and, thus, moisture availability at the land surface. Additionally, application of a free drainage condition at the bottom of the profile often applied in these models neglects possible boundedness, because of the presence of a shallow water table. Yet simple analytical solutions and stochastic analysis of steady-state evaporation in the presence of lateral and vertical heterogeneity demonstrated combined water table heterogeneity effects on evaporation and infiltration Zhu 2007, Zhu andMohanty 2002). This has also been confirmed by coupled, numerical simulations of a 1440 km 2 watershed using an integrated modeling approach including surface-subsurface hydrodynamics and a free water table (Maxwell et al 2007, Kollet andMaxwell 2008).
This study focuses on a description of the subgrid scale dynamic feedbacks between the subsurface and evapotranspiration at the land surface, which are not accounted for in previous steady-state studies. A Monte Carlo simulation of subsurface heterogeneity is performed using a large number of soil columns for two different soils (clay and sand) in conjunction with an integrated model including a free water table, and forced by a one-year atmospheric time series. In dynamic equilibrium mode (the energy and mass balances are zero over the entire year), instantaneous, hourly mass and energy fluxes for each individual column are calculated. The results from the ensemble of columns is supposed to represent the subgrid scale variability under the given water table and atmospheric forcing conditions and are interrogated in the context of geometric mean simulations, spatiotemporal variability, and degree of heterogeneity.

Methods
In this study, the coupled model ParFlow[CLM] was used to simulated the interactions between land surface processes and variably saturated flow in a heterogeneous subsurface. ParFlow [CLM] consists of the parallel, variably saturated subsurface flow model ParFlow (Jones andWoodward 2001, Kollet andMaxwell 2006) with integrated overland flow and the common land model (CLM; (Dai et al 2001)) that has been implemented modularly into ParFlow indicated by the brackets in ParFlow[CLM] (Kollet and Maxwell 2008, Maxwell and Kollet 2008, Maxwell and Miller 2005. CLM calculates the land surface mass and energy balances where R n is net radiation (W m −2 ); LE is latent heat (W m −2 ); H is sensible heat (W m −2 ); and G is ground heat (W m −2 ). The coupling with ParFlow allows transient, two-way feedbacks with 3D subsurface moisture redistribution processes including heterogeneity in the hydraulic parameters.
Since this study serves only as a starting point to characterize the influence of heterogeneity on the transient mass and energy fluxes at the land surface, some constraining assumptions were made as follows: (1) the subsurface heterogeneity consists of spatially uncorrelated, lognormal Gaussian random fields of saturated hydraulic conductivity K (L/T ) and also simple log-log correlation with the shape factor of the pressure-saturation function; (2) the water table is located close to the land surface (∼10 0 m); (3) subsurface energy transport is simulated by CLM and, thus, includes only the process of conduction (Kollet et al 2009); and (4) the atmospheric variables of the applied time series (temperature, humidity, wind speed, precipitation, barometric pressure, and radiation) are constant and do not change due to transient land surface conditions.
The vegetation type was implemented as crop using the default parameters in CLM. The atmospheric time series used in the simulation stems from a previous study by Kollet and Maxwell (2008) and was constructed from the North American Regional Reanalysis data set based on hourly time steps of the water year 1999. (Note that additional simulations were performed using a mid-latitude atmospheric forcing data set of the year 2002 from the Haarweg Meteorological Station of the Meteorology and Air Quality Section, Wageningen University and Research Center, The Netherlands, (http://www.met.wau. nl/haarwegdata/), which was previously applied by Kollet et al (2009) and is described in detail in Jacobs et al (2008). The results are only discussed briefly here.) The simulations were performed in dynamic equilibrium mode, which means that the atmospheric time series was applied successively over nine years until a dynamic equilibrium in the mass and energy variables was obtained.
The subsurface model was set up using isolated and also laterally interconnected columns with z = 0.025 m. The water table was implemented as a constant head boundary condition at the bottom of the domain with a saturated zone above to allow the root zone to extend below the water table. The shallow water table depth, d, varied between 1 and 2 m. These values are in the critical water table depth range as defined by (Kollet and Maxwell 2008) and reflect water table conditions under which vertical moisture transport is expected to strongly influence the mass and energy fluxes at the land surface. (In case of larger d, the water table is hydraulically disconnected from the land surface whereas for smaller d, moisture can be redistributed readily toward the land surface.) Simulations were performed on an ensemble of 100 soil columns with vertically varying K distributions based on the aforementioned uncorrelated Gaussian random fields. Additional simulations using the geometric mean of the K distribution were carried out for comparison. Two different soils, sand and clay, were approximated using the van Genuchten relationship (Van Genuchten 1980) for the pressure-saturation and relative conductivity functions. The constant van Genuchten parameters are summarized in table 1, implicitly assuming that K , α and n are uncorrelated. This has been supported in the literature, but is the subject of some debate (e.g., Smith and Diekkruger 1996). Table 1 also provides the geometric mean K g of the applied K distributions.
Additional fully 3D simulations were performed for uncorrelated ln K distributions for sand and clay using a water table depth, d = 1 m. The vertical discretization of z =  0.025 m was maintained and the lateral discretization was fixed at x = y = 1 m with nx = ny = 10, which leads to the same number of soil columns compared to 1D ensemble simulations. The rational was to study the influence of 3D variably saturated flow on LE variability influenced by soil heterogeneity.
In order to approximate potential correlation between K and n, a simple log-log regression model without white noise was used in the van Genuchten model Positive correlation was assumed and the regression constants a and b were chosen such that the mean value of n coincided with the value from the uncorrelated heterogeneity simulations. This approach assumes perfect linear correlation between K and n and is of course only a rough approximation of the complexity of the real system. However, it allowed an ad hoc evaluation of correlated K -n parameter distributions on the energy balance at the land surface. The hypothetical regression constants a and b are also listed in table 1. The simulations with the correlated K -n distributions were also performed under dynamic equilibrium and varying water table depth conditions.

Results and discussion
While this section focuses on latent heat (LE, evapotranspiration) as the variable of interest, changes in one variable consequently lead to changes in the other variables of the energy balance following energy conservation (equation (1)) enforced by CLM and the spin-up procedure.

Influence of K-heterogeneity and water table depth
Figures 1 and 2 show the weekly and monthly averages of the LE fluxes for different water table depths and heterogeneity configurations in the case of clay, respectively. The graphs contain the curves from the ensemble simulation, the ensemble mean, and the simulation using K g . For most of the simulated time period, heterogeneity does not influence the average fluxes, because the system is not significantly moisture limited due to adequate precipitation and the shallow water table. For water table depth d = 1.0 m and std ln K = 1.0 the largest ensemble variability occurs in the late summer for August and September. During these months and for this d-std ln K configuration, the K g simulation consistently overestimates the average LE fluxes. Increasing the heterogeneity to std ln K = 2.0 results in a similar overestimation of LE, however the ensemble variability decreases considerably. If the water table depth is increased to d = 2.0 m, heterogeneity in K does not seem to influence the results significantly on a weekly to monthly basis. In these cases it is interesting that the K g simulation coincides very well with the ensemble mean except maybe in May. The increase in water table depth also results in a decrease in LE in the summer months, which is more obvious in the plot of the monthly averages. Figures 3 and 4 show the weekly and monthly average fluxes for the same d-std ln K configurations in the case of sand, respectively. The overall impact of heterogeneity and varying, shallow water table depth on LE is very similar to the simulations of clay. For d = 1 m and std ln K = 1.0 the variability in weekly and monthly averages during the dry months of the year extends over a longer time period from July until October. For d = 1 m the overestimation of the K g simulation is more pronounced, which can be seen more clearly in the monthly averages.
In order to explain better the surprisingly good agreement of the ensemble and K g simulations over long time periods, daily differences between each ensemble mean and the K g simulation for the different d-std ln K configurations were calculated. Figure 5 shows that these differences basically behave like white noise most of the time except in the summer months, which explains the collapse of the longer term averages onto a single curve. Only for d = 1 m and std ln K = 1.0 is there a clear bias during the summer months with maximum values of up to 80 W m −2 (∼2.8 mm d −1 ). In the remaining months, daily differences may slightly exceed 20 W m −2 (∼0.7 mm d −1 ) in the case of sand, which generally produces larger differences for larger d in comparison to the simulation with clay. The behavior is basically the same for instantaneous hourly values (not shown here).
Several of these results are a consequence of the interconnection of soil characteristics, water table depths, and the applied atmospheric forcing. In CLM, LE consists of bare soil evaporation, E s , and transpiration by plants, E tr , which are both a function of soil moisture expressed in more or less complicated parameterizations. For example, E s is parameterized using the thermodynamically based formulation of soil air humidity by Edlefsen and Anderson (1943), which contains an exponential dependence of humidity on the soil matric potential. In the calculation of E tr , assimilation and water vapor transport rates are limited by the stomatal conductance which is a continuous function of the soil moisture content. Thus, soil heterogeneity influencing moisture transport will directly influence LE.
In the months from August until October (and to a limited extent also in May and other shorter time periods), evapotranspiration cannot be maintained by moisture from precipitation and canopy through-fall alone. Moisture must be redistributed from the shallow water table, which is a well-known process and has been discussed previously in great detail (e.g., Butler et al 2007, Chen and Hu 2004, Gardner 1958, Kollet and Maxwell 2008. In the presence of heterogeneity in K , vertical transfer is mainly limited by the layer with the smallest K ; the soil type; and the long term soil moisture profile depending on the location of the water table. This is expressed in the spread in the curves in figures 1-4 for d = 1.0 m and std ln K = 1.0. The large variability in LE stems from the relatively narrow ln K distribution and the shallow water table that limits but does not completely interrupt vertical moisture transfer. An increase in heterogeneity (increasing std ln K ) results in a higher probability that a low-K layer limits vertical transfer, which is reflected in the collapse of the ensemble basically onto a single curve and the strong overestimation by the K g simulation. If the water table depth increases from d = 1 to 2 m, the K g simulation provides a good average behavior of the system for std ln K = 1.0 and 2.0. However, this is not a general rule, because the parameterization in CLM does not lead to a continuous dependence of LE on soil moisture and behaves more like a threshold type formulation. Thus, if the K g simulation cannot maintain a certain soil moisture level with increasing distance from the land surface, instantaneous LE fluxes tend to zero, which may result in the same average fluxes compared to the simulations using a heterogeneous subsurface. The question remains whether this is an accurate representation of the naturally occurring process at the grid scale (>10 2 m). Figures 6 and 7 show the results for the simulations with correlated K -n distributions and 3D hydrodynamics, i.e. lateral moisture redistribution, for clay and sand, respectively. Although the simulations are not comprehensive, they demonstrate that correlation in K -n leads to a slight increase in LE variability and a more continuous distribution of the ensemble during the months of largest variability from July until October. This is reasonable, because the correlation of n with K results in a wider distribution of soil textures, which are now effectively spanning the different soil types from clay to sand and vice versa. However an increase of water table depth to d = 2 m leads again to a collapse of the ensembles onto basically a single curve, which does not need to be shown here because it is almost identical to figures 1 and 3.

Influence of K-n correlation and 3D hydrodynamics
The additional 3D simulations of variably saturated moisture redistribution were performed to juxtapose the results with the aforementioned results from the 1D simulations using isolated columns. Including shallow 3D hydrodynamics of the unsaturated zone apparently homogenizes the LE fluxes at the land surface (figures 6 and 7). This is reflected in a decreasing spread of the ensemble during the dry months of the year. Lateral interconnection of the previously disconnected soil columns at the 10 0 m scale leads to lateral moisture redistribution and the development of additional flow paths from the shallow water table toward the top soil layers. Thus, additional moisture is available for LE to maintain higher values at the land surface, which decreases the ensemble variance and produces a larger ensemble mean value that is quite similar to the K g simulation.

Influence of atmospheric forcing
It is important to keep in mind that the lower atmosphere is fixed in its forcing variables (temperature, humidity, wind speed, precipitation, barometric pressure, and radiation) by the applied atmospheric time series and constitutes the upper boundary condition in the simulations. Thus, there is no twoway feedback between the land surface and the atmospheric conditions in the simulations that may reciprocally influence land surface processes. From a thermodynamical point of view that may lead to inconsistent system responses i.e. physically unrealistic combinations of atmospheric conditions and ambient moisture and energy conditions at the land surface that are purely a function of the lower boundary condition of the soil column and the upper boundary condition represented by the atmospheric time series. In order to study the entire soil-vegetation-atmosphere system including feedbacks with the lower atmosphere, fully coupled simulations over longer time periods would be necessary, which are clearly beyond the scope of this study.
The same simulations were performed using an atmospheric time series from the Haarweg Meteorological Station, Wageningen, The Netherlands, for the year 2002, which is situated in a humid climate. The details of the site and time series are described in more detail in e.g., Jacobs et al (2008) and Kollet et al (2009). The analysis of the results is not shown here, because the influence of heterogeneity on evapotranspiration was negligible. In the year 2002, the soil-vegetation system was never water limited, because of continuous intermittent rainfall also in the summer months. Thus, evapotranspiration never relied on high rates of vertical water transfer from the shallow water table, which may be limited in the case of soil heterogeneity. However, the influence of soil heterogeneity may become more pronounced under conditions with less rainfall as happened for example in the year 2003. Additional simulations are necessary in the future for clarification purposes.

Summary and conclusions
In this study, the influence of subgrid heterogeneity in the shallow hydraulic soil properties on evapotranspiration or latent heat, LE, was studied for two different soil types (sand and clay) and under varying water table conditions. In the case of inadequate moisture availability, soil type and subgrid heterogeneity do exert a major influence on LE.  Ensembles from Monte Carlo simulations showed the largest variability during the dry months of the year, from August until September in the case of clay and from July until October in the case of sand. In the remaining months variability is negligible in terms of weekly and monthly averages, though instantaneous daily and hourly differences can still be on the order of 10 1 -10 2 W m −2 (10 −1 -10 0 mm d −1 ). For shallow water table depths, d, simulations using the geometric mean of the ln K -distribution strongly overestimate the ensemble mean. However, the geometric mean and ensemble basically collapse onto a single curve if d is increased to 2 m for both soils. Additionally, LE values decrease considerably during the dry months reflecting the dependence of LE on d, which has been shown previously. The small influence of heterogeneity on weekly and monthly averages can be explained with the behavior of the instantaneous and daily averaged differences of the LE fluxes that behave very much like white noise and only show a bias during the dry months of the year, when the influence of heterogeneity on LE becomes important.  Correlated distributions of K and the shape factor n in the van Genuchten model result in a more continuous distribution of the ensemble during the dry months of the year, but do not have a pronounced effect on individual LE fluxes and the overall ensemble behavior. On the other hand, 3D hydrodynamics of the shallow unsaturated zone homogenize the ensemble considerably through lateral moisture transport and the development of additional flow paths that are connected with the top soil layers. The latter conclusion is valid only for conditions when topography does not influence moisture transport. Additional hill slope simulations are necessary to account for gravitationally driven 3D shallow moisture transport.
The findings have implications for simulations from the local to the global scale. Under humid climate conditions soil heterogeneity is expected to have an influence on the mass and energy balances only during extended drying periods. However, in (semi-) arid regions characterization and upscaling of subgrid scale soil heterogeneity poses a challenge, because it exerts a large influence on the mass and energy fluxes at the land surface. In the presented simulations, application of the geometric mean of the applied K -distribution resulted in a relatively strong overestimation of LE. This must be kept in mind in the analysis of global results that do not account for subgrid scale soil heterogeneity. Since it will not be computationally feasible in the near future to perform Monte Carlo simulations at the global scale, new generic laws of upscaling must be found that better reflect the mean behavior of the modeled system.
The diagnosed influence of soil heterogeneity on the LE fluxes is also a result of the applied atmospheric time series, which is constant and does not adjust to the spatially changing moisture and energy conditions at the land surface due to soil heterogeneity. That is, the fixed boundary condition determines much of the behavior of the system a priori. In order to account for two-way feedbacks with the atmosphere, detailed simulations of the boundary layer immediately adjacent to the land surface and at the spatial scale of the soil heterogeneity would be required. In the case of simple land surface configurations (e.g., bare soil) large eddy simulations may be an option. However, in the presence of vegetation, novel land surface parameterization schemes need to be developed that relax the assumption of the commonly applied 1D similarity approaches and account for lateral moisture, energy and momentum transport above the land surface.