Water transit times in the unsaturated zone: Spatio-temporal variation and its application in for the characterization of functional units

Abstract:

Introduction

In catchment hydrology, the pedosphere takes a special role due to filtering, buffering and transformation of water and solutes between the atmosphere, the ground water and the vegetation. However, tracing the fate of rain water over time and space on its way through the soil is challenging. Natural water isotopes - as ideal tracers - have been shown to serve as useful tool for a better process understanding in catchment hydrology in general (Hrachowitz et al., 2011; Kendall and McDonnell, 1998) and in soil hydrology in particular (Eichler, 1966; Garvelmann et al., 2012; Koeniger et al., 2010; Zimmermann et al., 1966). Stable water isotopes in runoff were also successfully used for calibrating and benchmarking catchment models (Birkel et al., 2010; Fenicia et al., 2010; Hartmann et al., 2012; Heidbüchel et al., 2013). In addition, isotope data can help in the parameterization of soil physical models (Stumpp et al., 2012). In general it was shown, that taking isotope data in addition to hydrometric data (e.g. discharge or soil moisture) into account provides additional information that allows estimating not only the response time but also the transit times of the system. The concept of transit times, as an integrative catchment characteristic that describes the storage, flow pathways and source of water (McGuire and McDonnell, 2006) has been increasingly used in the last decade. While the early studies focused on the mean transit time (MTT) and their governing hydrological or landscape characteristics (McDonnell et al.; McGlynn et al., 2003; McGuire et al., 2005; Rodgers et al., 2005; Soulsby et al., 2006), the temporal variability of the transit times are nowadays also taken into account. Hence, the probability of water particles of a precipitation event to pass the hydrological system, as the transit time distribution (TTD), and its variance over time are subject to current research in catchment hydrology (Botter et al., 2011; Heidbüchel et al., 2013; Heidbüchel et al., 2012; Hrachowitz et al., 2013; Rinaldo et al., 2011). Furthermore, various studies showed that the assumption of well mixed conditions, which is usually made in determining MTT via an empirical function, does not hold for most of the hydrological systems (Botter, 2012; Godsey et al., 2010).

However, the role of the unsaturated zone with regard to transport, mixing, and release of water and solutes is not well understood. Especially a better knowledge about the variation of the hydrological processes in the subsurface governed by various characteristics in topography, land use, and pedology across the catchment might help to get insights if landscape entities of similar hydrological functioning exist and if certain areas dominate the transit time distributions at the catchment scale. Therefore, this study aims to work on the following hypothesis: 1.) Catchments consist of elementary functional units (EFU), which can be delineated with regard to the hydrological response in the vadose zone and characterized by transit time distributions on the pedoscale. 2.) The hydrological response in the vadose zone of EFUs is time variant, governed by the catchment state and its forcing (connectivity, storage, precipitation, solar radiation).

To test these hypotheses, we first derived the water and solute transport parameters of a physical-based soil model with the means of site specific soil moisture and pore water isotope data for 46 study sites within the Attert catchment. Afterwards, we applied the model to trace the fate of the water in the unsaturated zone, including percolation through the soil matrix, seepage fluxes at the bottom of the soil profile, and losses due to evaporation and root water uptake. The resulting TTDs for the water leaving the profile via seepage flux or transpiration are then compared with respect to EFU characteristics, type and intensity of driving forces as well as system state/initial conditions.

Study sites

The 46 sites are located in the Attert catchment in Luxembourg. The Attert catchment is characterized by three different geological settings: schist, sandstone, and marls in between the two other areas. The sites were instrumented with sensor clusters during the last 2.5 years within the framework of the CAOS project. At all sites basic meteorological data, soil moisture, soil temperature, matrix potential, sapflow and groundwater table dynamics are monitored. The site locations were chosen to cover a range of possible controls such as the geology, land use, and topography, in order to test whether there are sites of similar hydrological behavior. 23 of the 46 sensor clusters were installed in the schist area, 12 in the sandstone and 11 in the marl areas. The prevailing soils developed on schist are structured loamy Cambisols, but Gleyosols are present in groundwater influenced areas. On marl, clayey Stagnosols have developed and in the sandstone area, sandy entic Podzols dominate. Since the study sites lie within an area of about 30 km2, the climatic controls are assumed to be similar.

Methods

Data

Pore water stable isotopes information at different depths in the soil were sampled during the instrumentation of the sensor clusters between February 2012 and October 2013 and analyzed for their pore water isotopic composition according to the equilibration method as proposed by Wassenaar et al. (2008). Each isotope profile was determined by taking soil samples in 5 cm intervals from a soil core excavated with a percussion gauge. Profile depths ranged from 120 to 380 cm, depending on depth to water table or depth to bedrock. The soil samples were taken to the laboratory in sealed air tight bags. In addition to the soil samples, standards were prepared, which consisted of oven-dried soil material from the study sites. This material was rewetted to the soil moisture of the sampling moment with three different waters of known isotopic composition. The bags containing the samples or the standards and the sampleswere then filled with dry air and the soil pore water was allowed to equilibrate with the dry atmosphere in the bag for two days under constant temperature (21°C). Then the head space in the bags was directly sampled with a Wavelength-Scanned Cavity Ring Down Spectroscope (Picarro, Santa Clara, USA) for 6 minutes, and the measured concentration of deuterium and 18O during the last 90 seconds was averaged to minimize carryover effects. The isotopic composition of the gas phase was converted to values of the liquid pore water according to the temperature dependent fractionation factor as defined by Majoube (1971). Standards were measured at the beginning, every three hours during, and at the end of the analyses for each profile. The standards were used to account for drift of the laser spectrometer and to calibrate the measurements in order to get values in the δ notation relative to the Vienna Standard Mean Ocean Water (VSMOW).

At the cluster sites, soil moisture was measured continuously in three profiles at 10, 30, and 50 cm soil depth. Time series of daily averages were calculated for each depth and each study site. Precipitation and air temperature data was available since 2005 for the three different geological areas provided by the meteorological stations run by the CRP-GL Lippmann Institute. The isotopic composition of the rainfall is measured at least every 14 days since November 2009. To minimize the influence of the initial conditions of the deuterium concentration in the pore water, the time series of isotope concentration of the precipitation were extended with additional isotope data from GNIP stations for the time period between 2008 and end of 2009.

Modeling

The transient water flow and transport within the unsaturated soil profile was simulated by numerically solving the Richards equation with the finite-element code of Hydrus-1D (Šimůnek et al., 2012). Since solute losses at the upper boundary due to evaporation is not accounted for in the original model code, a modified version of Hydrus introduced by Stumpp et al. (2012) was used which allows the tracer to evaporate. Potential evapotranspiration (PET) was calculated with the Hargreaves Formula, the root water uptake was included according to Feddes et al. (1978) and a snow module based on the degree-day method was applied. The rooting depth and the leaf area index of the vegetation were defined according to field observations, taking the seasonality of the grass and deciduous trees into account.

The required soil physical parameters that describe the water retention and hydraulic conductivity as presented by van Genuchten (1980) and the diffusivity parameter which describes the dispersion in the advection-dispersion model were determined at every study site via inverse modeling. At every site, the depth of the soil profiles for which the simulation were done was set to 200 cm and were discretized into 101 nodes. In order to account for different soil horizons, the profile was divided into two different layers, for each of which six parameters had to be optimized to simulate the water and solute transport of the unsaturated zone: On the one hand the five parameters (θr, θs, α, n, Ks) that describe the water retention and hydraulic conductivity characteristics in accordance to the Mualem-van Genuchten model (MVG) and on the other hand the longitudinal dispersivity (λ), describing the dispersion of a solute were determined.

The global optima of the parameter spaces were determined with the Shuffled-Complex-Evolution algorithm (SCE-UA) (Duan et al., 1992). The modified Kling-Gupta-Efficiency (KGE) (Kling et al., 2012) was applied as the objective function in a multi-objective optimization process: The measured soil moisture time series and isotope profiles were used to simultaneously optimize the parameter for the water and deuterium transport. The KGE for the soil moisture was computed for each soil depth and then averaged (KGE_θ). The two fitting targets, soil moisture and pore water isotope data, were equally balanced, because the KGE was calculated from the average over the efficiencies of the simulated soil moisture series and the isotope profiles (KGE_total = (KGE_θ + KGE_D)/2).

The best performing parameter sets were then used to simulate the fate of water particles in the soil profile. For each day with precipitation, the precipitation water of that day was labeled with a virtual ideal tracer, which allows tracking the water particles in the soil profile. As a result, breakthrough curves (BTC) of water from individual precipitation events leaving the soil profile in 200 cm depth (event transit time distributions (ETTD)) were derived for each rainy day between 2006 and 2011. The median transit time (MTT) of that seepage water was calculated as the median of the BTC for every day with precipitation. Likewise, the MTT of transpired water was derived as the median of the probability of a water particle leaving the soil via root water uptake. The fraction between the precipitation water that leaves the soil via seepage fluxes and the water that is taken up by the roots were also calculated as the maximum of the cumulative distribution function of the ETTD. The ETTDs were multiplied by the event size, superimposed, summed up and divided by the total precipitation to derive site specific master transit time distributions (MTTD) (Heidbüchel et al., 2012). The MTTDs, which represent the average probability of a water particle to pass the -200 cm soil depth plane, were grouped in accordance to the season, to get seasonally specific MTTDs.

Results

The results are shown for three exemplary study sites of the three different geological settings. The above presented methodology will be applied to all 46 sensor cluster sites until the Chapman conference in September.

The results of the optimization process with regard to the performance of the best parameter sets to simulate the soil moisture dynamics and the pore water isotope data are shown in Figure 1 and Figure 2, respectively. The soil moisture dynamics between March 2013 and January 2014 are reproduced well for the Cambisol on schist and the Podzol on sandstone with values of KGE_θof 0.69 and 0.65, respectively. For the Stagnosol on marls, the soil moisture in the upper soil layer, represented by measurements in -10 cm depth, is not simulated in a satisfactory manner. The simulated soil moisture in -30 and -50 cm is in the magnitude of the measurements, but the intense drying during the fall in 2013 is not well represented in the modeling results. The reasons for this mismatch are unclear and will be investigated.

The isotope profiles that were taken at the end of March 2013 are simulated by the soil physical model with some limitations. For example, the seasonal variation, which is preserved in the sampled isotope profile of the Cambisol on schist, is not present in the simulations due to a high dispersivity parameter. A similar effect can be seen for the Stagnosol on marls, where the peak of isotopic enriched precipitation of the last summer is underestimated in the simulation. The best fit is achieved for the Podzol on sandstone, although there are anomalies in -120 to -160 cm soil depth, which might be caused by lateral flow at the hillslope.

The time variant MTTs of the seepage and the transpiration water are presented in Figure 3 for the three exemplary study sites. In addition, the time variant fraction between the amounts of event water that leave the soil profile via seepage and evapotranspiration (ET) is shown in the same figure. The MTTs of seepage water are lowest in the Cambisol ranging between 100 to 500 days. The highest transit times occur in the Stagnosol on marls, where the MTT of seepage water is between 600 and 1150 days. The Podzol on sandstone has medium drainage rates with MTTs of 250 to 750 days. At all the sites, the lowest MTTs occur for water from precipitation events in late fall, because ET is lowest in fall and winter. These circumstances also explain the high fraction of water leaving the soil via seepage and not by ET. High MTTs of seepage are less influenced by seasonality, but more dependent on rainfall dynamics. The root water dynamics are - of course - driven by the physiological processes of the vegetation. As such, the MTT of water to leave the soil via transpiration is highest during the dormant season and approaches sub-weekly values during the summer and beginning of fall.

The MTTDs, which represent the probability of a water particle to pass the -200 cm soil depth plane, are presented in Figure 4. The MTTDs show that the water percolating through the Cambisol on schist has not only the shortest transit times, but also the lowest dispersion. Within the Stagnosol, much more dispersion occurs. Again, the Podzol on sandstone occupies a medium position between the Cambisol and Stagnosol with regard to transit times and dispersion.

Conclusion

Our study demonstrates the potential of site specific transit time distributions for comparison between the responses of different soils developed under different conditions. Therefore, this approach will be expanded to the other cluster sites with different slope positions, land-use and soil properties to further explore how these sites differ in relation to their hydrological functioning. We have shown that the inclusion of pore water isotope data in the calibration of physical soil hydrological models allows simulating the transition of the water through the soil. The first results for the three different geologies support the existing conceptual models and field measurements (van den Bos et al., 2006): High percolation rates within the Cambisols on schist due to its structured characteristics with high rock contents and preferential flow paths; low water flows through the clayey Stagnosols on marl, where dispersion plays a major role; relatively high percolation rates in the homogenous sandy Podzols on the sandstone. The consideration of different catchment states point out the temporal variability of the hydrological functioning. However, it was shown that it is important to distinguish between seasonal variation (governed by radiation forcing, evaporation and dormant season of the vegetation) and non-seasonal variation due to forcing of the system via precipitation.

The above presented methodology and its application for the exemplary study sites show that transit time distributions provide promising opportunities to compare hydrological processes across space. An extension of the study to all 46 sites within the CAOS framework will show general differences and their governing site specific characteristics. Similarity measures such as the Earth Movers Distance to compare probability density functions with each other and the visualization via non-metric multi dimensional scaling will contribute to the search and possible delineation for EFUs. Acknowledgements

First author was funded by the DFG Research Group: From Catchments as Organised Systems to Models based on Functional Units (FOR 1598). The isotope data in the precipitation was provided by FNR/CORE/SOWAT project of the Département 'EnVironnement et Agro-biotechnologies' CRP-GL Lippmann Institute. Sampling of the isotope profiles was made possible by the support of the CAOS-Team and Begona Lorente Sistiaga, Benjamin Gralher, Andre Böker, Marvin Reich, Andrea Popp. Special thanks to Britta Kattenstroth Jean Francois Iffly for their technical support in the field and Barbara Herbstritt for her support in the laboratory.

References

Birkel, C., Dunn, S.M., Tetzlaff, D., Soulsby, C., 2010. Assessing the value of high-resolution isotope tracer data in the stepwise development of a lumped conceptual rainfall-runoff model. Hydrol Process 24, pp. 2335-2348.

Botter, G., 2012. Catchment mixing processes and travel time distributions. Water Resour. Res. 48.

Botter, G., Bertuzzo, E., Rinaldo, A., 2011. Catchment residence and travel time distributions: The master equation. Geophys. Res. Lett. 38, pp. L11403.

Eichler, R., 1966. Deuterium-Isotopengeochemie des Grund- und Oberflächenwassers. Geol Rundsch 55, pp. 144-159.

Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of field water use and crop yield. Centre for Agricultural Publishing and Documentation, Wageningen. Simulation monographs.

Fenicia, F., Wrede, S., Kavetski, D., Pfister, L., Hoffmann, L., Savenije, H.H.G. et al, 2010. Assessing the impact of mixing assumptions on the estimation of streamwater mean residence time. Hydrol Process 24, pp. 1730-1741.

Garvelmann, J., Külls, C., Weiler, M., 2012. A porewater-based stable isotope approach for the investigation of subsurface hydrological processes. Hydrol. Earth Syst. Sci. 16, pp. 631-640.

Godsey, S.E., Aas, W., Clair, T.A., Wit, H.A. de, Fernandez, I.J., Kahl, J.S. et al, 2010. Generality of fractal 1/f scaling in catchment tracer time series, and its implications for catchment travel time distributions. Hydrol Process 24, pp. 1660-1671.

Hartmann, A., Kralik, M., Humer, F., Lange, J., Weiler, M., 2012. Identification of a karst system’s intrinsic hydrodynamic parameters: upscaling from single springs to the whole aquifer. Environ Earth Sci 65, pp. 2377-2389.

Heidbüchel, I., Troch, P.A., Lyon, S.W., 2013. Separating physical and meteorological controls of variable transit times in zero-order catchments. Water Resour. Res. 49, pp. 7644-7657.

Heidbüchel, I., Troch, P.A., Lyon, S.W., Weiler, M., 2012. The master transit time distribution of variable flow systems. Water Resour. Res. 48.

Hrachowitz, M., Bohte, R., Mul, M.L., Bogaard, T.A., Savenije, H.H.G., Uhlenbrook, S., 2011. On the value of combined event runoff and tracer analysis to improve understanding of catchment functioning in a data-scarce semi-arid area. Hydrol Earth Syst Sc 15, pp. 2007-2024.

Hrachowitz, M., Savenije, H., Bogaard, T.A., Tetzlaff, D., Soulsby, C., 2013. What can flux tracking teach us about water age distribution patterns and their temporal dynamics? Hydrol. Earth Syst. Sci. 17, pp. 533-564.

Kendall, C., McDonnell, J.J., 1998. Isotope tracers in catchment hydrology. Elsevier, Amsterdam.

Koeniger, P., Leibundgut, C., Link, T., Marshall, J.D., 2010. Stable isotopes applied as water tracers in column and field studies. Stable Isotopes in Biogeosciences (III) 41, pp. 31-40.

Majoube, M., 1971. Fractionnement en oxygene-18 et en deuterium entre l’eau et sa vapeur. J. Chim. phys 68, pp. 1423-1436.

McDonnell, J., Rowe, L., Stewart, M.K.1999. A combined tracer-hydrometric approach to assess the effect of catchment scale on water flow path, source and age, in: Leibundgut, C., McDonnell, J., Schultz, G. (Eds.). IAHS, General Mean Residence Time, pp. 265-273.

McGlynn, B., McDonnell, J., Stewart, M., Seibert, J., 2003. On the relationships between catchment scale and streamwater mean residence time. Hydrol Process 17, pp. 175-181.

McGuire, K.J., McDonnell, J.J., 2006. A review and evaluation of catchment transit time modeling. J Hydrol 330, pp. 543-563.

McGuire, K.J., McDonnell, J.J., Weiler, M., Kendall, C., McGlynn, B.L., Welker, J.M. et al, 2005. The role of topography on catchment-scale water residence time. Water Resour Res 41.

Rinaldo, A., Beven, K.J., Bertuzzo, E., Nicotina, L., Davies, J., Fiori, A. et al, 2011. Catchment travel time distributions and water flow in soils. Water Resour Res 47.

Rodgers, P., Soulsby, C., Waldron, S., 2005. Stable isotope tracers as diagnostic tools in upscaling flow path understanding and residence time estimates in a mountainous mesoscale catchment. Hydrol Process 19, pp. 2291-2307.

Šimůnek, J., M. Sejna, H. Saito, M. Sakai, M. Th. van Genuchten, 2012. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media, Version 4.15, Riverside, California.

Soulsby, C., Tetzlaff, D., Rodgers, P., Dunn, S., Waldron, S., 2006. Runoff processes, stream water residence times and controlling landscape characteristics in a mesoscale catchment: An initial evaluation. J Hydrol 325, pp. 197-221.

Stumpp, C., Stichler, W., Kandolf, M., Šimůnek, J., 2012. Effects of Land Cover and Fertilization Method on Water Flow and Solute Transport in Five Lysimeters: A Long-Term Study Using Stable Water Isotopes. Vadose Zone Journal 11

van den Bos, R., Hoffmann, L., Juilleret, J., Matgen. P., Pfister, L.2006. Conceptual modelling of individuals as a trade-off between bottom-up and top-down modelling, a case study, iEMSs Third Biennial Meeting. "Summit on Environmental Modelling and Software".

van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44, pp. 892-898.

Wassenaar, L., Hendry, M., Chostner, V., Lis, G., 2008. High Resolution Pore Water δ 2 H and δ 18 O Measurements by H 2 O (liquid) −H 2 O (vapor) Equilibration Laser Spectroscopy. Environ. Sci. Technol. 42, pp. 9262-9267.

Zimmermann, U., Munnich, K.O., Roether, W., Kreutz, W., Schubach, K., Siegel, O., 1966. Tracers Determine Movement of Soil Moisture and Evapotranspiration. Science 152, pp. 346-347.