Deriving Effective Soil Water Retention Characteristicsfrom Shallow Water Table Fluctuations in Peatlands

The characterization of soil hydraulic properties, including their field variability and scale dependency, is an ongoing research challenge in soil science for decades (Jury et al., 2011). Still, most commonly, soil hydraulic properties are obtained for relatively small soil samples by laboratory measurements imposing steady-state (e.g., hanging water column) (Durner and Lipsius, 2005; Reynolds et al., 2002b) or transient conditions (e.g., evaporation experiments) (Dettmann et al., 2014). However, hydraulic properties obtained in the laboratory may not be representative for field conditions and often disagree with hydraulic properties determined in situ (Basile et al., 2003). One reason for this disagreement is that in the laboratory higher saturation is often achieved compared with field conditions, under which the portion of entrapped air is higher. Another reason for deviation is sample size. For practical reasons and to achieve a sufficient number of replicates, in many studies, soil sample size is chosen too small to sample a representative elementary volume of the soil heterogeneity that is characteristic for the specific soil. Additionally, samples may be disturbed by the sampling procedure, e.g., due to compaction or disturbance of the soil structure. Further, soil samples may behave differently in the laboratory than the soil under field conditions, as is the case for clay and peat that present different shrinkage characteristics once cut from the coherent soil body and root system (Mitchell, 1991; Mitchell and van Genuchten, 1992). In peatlands, shrinkage above and compression below the water table can have considerable influence on how the soil profile stores and releases water when the water table fluctuates (Price and Schlotzhauer, 1999). Further shrinkage and swelling changes the pore size distribution, and soil hydraulic properties are thus temporally variable (Chow et al., 1992; Mapa et al., 1986). Such effects pose limitations on the transfer of laboratory parameters to the field. Therefore, for applications at the field scale, in situ measurements provide important advantages over laboratory measurements because they are more representative (Paquet et al., 1993) and obtain data in the natural environment, including interactions between different soil layers and scales (Wollschläger et al., 2009).

Several methods exist for the hydraulic characterization of a soil under field conditions. The saturated hydraulic conductivity can be measured with bail tests (Hvorslev, 1951) or by infiltration-based methods, e.g., the double-ring infiltrometer (Reynolds et al., 2002a). A tension disk or pressure ring infiltrometer can be used to obtain the unsaturated hydraulic conductivity at specific pressure heads (Ankeny et al., 1991; Basile et al., 2003). This method is applicable for low suctions (approximately more than −25 cm) (Bodhinayake et al., 2004). The direct in situ determination of the soil water retention characteristics (WRC) is limited to simultaneous measurements of water content (θ) and pressure head. Because these are point-like measurements, this approach requires an appropriate number of replicates per horizon. Additionally, measurements of θ often require soil-specific calibrations, especially for soils with high soil organic carbon contents and distinctive shrinkage and swelling characteristics (Nagare et al., 2011; Pepin et al., 1992; Shibchurn et al., 2005). An alternative to direct in situ measurements is the indirect determination of hydraulic parameters with inverse optimization using in situ measured state variables (Jadoon et al., 2012; Wollschläger et al., 2009). However, this approach requires various measurements as input (precipitation, evaporation, surface- and groundwater in- and outputs) besides the observed state variables.

In peatlands, the soil hydraulic properties are crucially influencing water table fluctuations and specific hydrological conditions and therefore physical, chemical, and biological processes (Dimitrov et al., 2010; Holden et al., 2004; Lafleur et al., 2005; McLaughlin and Cohen, 2014; Waddington et al., 2015). Unfortunately, for shallow groundwater systems, knowledge about soil WRC at the field scale is scarce due to the difficulties of common in situ methods. Infiltration methods are problematic due to the influence of the shallow water table that lowers infiltration. Furthermore, accurate in situ measurements of θ are difficult to obtain in wetlands with high soil organic carbon contents (Mortl et al., 2011), which is of concern for both the direct and inverse determination of the soil WRC. On the contrary, water table and standard meteorological data are easy to obtain and widely available for shallow groundwater systems. As mentioned above, water table depth dynamics as a response to boundary fluxes contain information about the WRC of a soil. To our knowledge, this information has not yet been exploited to inversely estimate soil the WRC.

Many kinds of landscapes with shallow groundwater tables are characterized by a distinctive microrelief that has complex effects on hydrological processes (Van der Ploeg et al., 2012). Water table dynamics at high water tables are strongly affected by a partially inundated, uneven soil surface that increases the specific yield of the system and permits the initiation of surface runoff before the whole soil surface is inundated. Also, at lower water tables in the absence of partial inundation, the specific yield is influenced by the microrelief through the non-uniform vertical distribution of the soil volume and the resulting effect on soil water retention (Dettmann and Bechtold, 2016). Dettmann and Bechtold (2016) gave a one-dimensional analytical expression combining the soil WRC and the microrelief effect on the spatially averaged specific yield (S_y). For its application, information about the microrelief as cumulative surface elevation is needed. Besides ground-based surveys, there is an increasing availability of detailed digital elevation models from laser scanning, which can be used to characterize the microrelief.

In this study, we developed an approach to inversely estimate the soil WRC of shallow groundwater systems from frequently available data on water table fluctuations, precipitation (P), and microrelief. Instead of using a continuous soil hydrological model, we focused only on periods of stronger P events. Our approach is built on two main assumptions: (i) for shallow groundwater systems, the soil moisture profile is close to hydrostatic equilibrium before and after rain events; and (ii) during short time periods of high precipitation, the water storage change due to the divergence of lateral fluxes is small compared with the precipitation input. Given these assumptions, water table rises are directly linked to the P amounts, and the height of the water table rise depends on the soil WRC, the initial water table, and the frequency distribution of the microrelief. We applied and evaluated the approach with field data from an ombrotrophic sphagnum bog complex with shallow groundwater tables.

The value of ΔA_soil decreases with shallower water tables. When zl and zu are above ground, i.e., for periods of total inundation, P is equal to Δz = zu − zl, i.e., the height difference of two open water surfaces, further referred as ΔA_surface. Following this, the amount of water received by a system for a depth increment between zl and zu can be separated into ΔA_soil and ΔA_surface, with an abrupt transition from ΔA_soil to ΔA_surface for even surfaces.

The soil moisture profile is vertically extended above the maximum height of the surface elevation by the upper integral bound being infinity, for which the cumulative frequency distribution of the microrelief is 1 and the effective soil moisture profile is 0.

In this study, we did not deal with microrelief issues but with the definition of the soil–vegetation interface at the “point” scale. For many environments, e.g., the peat moss and peat layer in sphagnum bog ecosystems, it is hard to determine a clear position of the soil–vegetation interface due to the continuous vertical transition from soil to vegetation (indicated in Fig. 1). When placing a level staff onto the ground, a position is measured at which the penetration resistance increases to a degree that the level staff is not further penetrating into the soil. In this study, for the modeling, we looked for a position of the soil–vegetation interface that optimally separated the two water storage volumes, ΔA_soil and ΔA_surface. This “optimal” position may not be consistent with the absolute level determined by the survey measurement, as for the water storage modeling the interface is best set where the strongest increase of large easily drainable macropores occurs, which commonly increase from depth to the surface (Moore et al., 2015; Morris et al., 2015). This nearly coincides with the strongest increase in S_y. Because this position is not known beforehand, the position of the interface was a parameter in the inversion. We emphasize that this parameter only adjusts the absolute level of the surface elevations, but a survey measurement is still needed for the frequency distribution of the surface elevations. We refer to the survey-measured mean surface elevation with parameter μ, which is defined to be zero in all model applications. Parameter Δμ (indicated in Fig. 1) is the optimized position of the interface of ΔA_soil and ΔA_surface relative to μ (z = 0), i.e., Δμ is negative when the soil–vegetation interface is lower than the survey elevation. Parameter Δμ can be included in Eq. [2] to [4] by substituting F_s(z) with F_s(z − Δμ). Note that, for clarity, all given and shown water table depths here are relative to μ. Parameter Δμ is given as an additional parameter.

The vG model was chosen because Dettmann et al. (2014) already investigated the sphagnum peat of our study site with laboratory evaporation experiments and inverse optimization. The results showed that the vG model was suitable to describe the peat soil moisture dynamics under laboratory conditions.

For global optimization, we used the Markov chain Monte Carlo algorithm DifferRential Evolution Adaptive Metropolis (DREAM) (Vrugt et al., 2008, 2009a, 2009b). The algorithm evolves a posterior density function (pdf) of individual parameters that are treated as probabilistic variables considering the observed data set. Starting with an initial population within the feasible parameter space (prior distribution), the pdfs are evolved in multiple individual Markov chains combining the prior distribution and the data likelihood. Further information about the algorithm can be found in various publications (Vrugt et al., 2008, 2009a, 2009b; Vrugt and Ter Braak, 2011) and is not repeated here.

For all model applications, the prior distributions of the vG parameters were set as a uniform distribution limited by lower and upper parameter boundaries that cover most soil types including peat (Table 1). The upper boundary of θ_s − θ_r was constrained a little more to a value of 0.9, which is lower than reported values for sphagnum peat (0.97; Paavilainen and Päivänen, 1995). The value of 0.9 relies on θ_s values determined by Dettmann et al. (2014) near the central monitoring well and allowed a consistent comparison of the field-derived with the laboratory-derived retention curve from this study, i.e., both retention curves start at a value of 0.9 and could not fall below θ < 0. Due to the spatial heterogeneity around a monitoring well and within the whole study site, θ_s could differ from the determined value. However, we expected differences in θ_s to be small among the three monitoring wells because, in general, the observed differences in peat decomposition were rather small. Furthermore, laboratory-determined values of θ_s are commonly higher than θ_s values that are reached under field conditions due to wetting by capillary rise compared with wetting by precipitation, which commonly results in a higher fraction of entrapped air. Thus, we think that 0.9 is justifiable as the upper boundary of θ_s − θ_r for the field conditions of our site.

The standard deviation (σ) of the normal distribution F_s(z) of the surface elevations was set a priori individually for the three monitoring wells according to the values listed in Table 2. These values relied on the survey data.

DREAM was run with the standard configuration, evolving four parallel chains. The uniform prior distribution was sampled using Latin hypercube sampling. As convergence criteria, the Gelman and Rubin (1992)  convergence diagnostic was set to < 1.1. After the target distribution achieved convergence, DREAM was run for 50,000 additional function evaluations to generate the pdf.

We applied our approach to continuous water table time series from three monitoring wells (south, central, and north) located in a near-natural ombrotrophic bog field site (Schechenfilz, 47°48′ N, 11°19′ E, Germany). All monitoring wells were located at the plateau of the raised bog. The bog has a 5- to 6-m-thick peat layer mainly under permanently water-saturated conditions. The peat substrate and current vegetation is dominated by sphagnum mosses, with spatially variable contributions of sedges, heather meadows, and bog pines (Hommeltenberg et al., 2014). Vegetation around the south monitoring well is only sphagnum moss. The central monitoring well is surrounded by a combination of sphagnum mosses, sedges, heather meadows, and bog pines. The vegetation composition around the north monitoring well is similar to that of the central monitoring well, except the bog pines are missing.

The northern part of the bog was partly affected by peat cutting until the 1950s (Hommeltenberg et al., 2014). Because the peat cutting occurred only in parts, the peat layer is spatially variable, showing different stages of decomposition. However, in most parts the bog is still rather pristine. The peat near the central monitoring well was classified with H2 (0–0.02 m), H1 (0.02–0.12 m), and H7 (0.12–1.1 m) on the von Post scale (Drösler et al., 2015), which classifies the degree of peat humification based on the proportion of visible plant remains and the soil water color (von Post and Granlund, 1926). At the north and south monitoring wells, the soil profiles have not been classified with the von Post scale. Around the north monitoring well, relict and refilled ditches can be observed on satellite pictures. Slightly higher degrees of decomposition cannot be excluded around this monitoring well. The south of the bog complex was not affected by ditch drainage, and the sphagnum peat around the south monitoring well is pristine.

Water tables of all three monitoring wells are close to the surface most of the year, with a mean water table for the available time series (15 Sept. 2010–11 Nov. 2014) of −0.15 m (south monitoring well), −0.18 m (central monitoring well), and −0.04 m (north monitoring well). The minimum water tables are −0.44 m (south and central monitoring wells) and −0.32 m (north monitoring well) (see also Table 2). In our analysis, we used water table rises from zl to zu as a response to P in the following referred as events, with Δz >0.02 m, zu <0 m, mean slope of the rise >0.004 m h⁻¹, and a total sum of P >2 mm.

The Schechenfilz bog complex has a considerable microrelief consisting of hummocks and hollows (Nungesser, 2003). Because water table changes are affected by the microrelief, the surface elevation was measured around the three investigated monitoring wells. Along six 8-m transects, each beginning at the monitoring well, elevations were measured every 0.25 m in relation to the elevation of the monitoring well. However, as detailed above, sphagnum bogs are characterized by a continuous transition from peat soil to vegetation (Fig. 1), and it is hardly possible to define and measure a soil–vegetation interface. Thus, to obtain objective and reproducible elevation data, we placed a level staff (locating surface: 17 cm²; weight: 1.8 kg) onto the ground without manually adding any extra pressure.

Figure 2 shows the surface elevation distributions as an example around the south monitoring well and the fitted normal frequency distribution and F_s used in the inversion. Usually water table depths are determined with respect to the surface height at the location of the monitoring well. However, this makes water table depths dependent on the specific microtopographic position of the monitoring well, and the comparison of water table depths with other monitoring wells in the study site is often more influenced by this arbitrary position and less by real differences in hydrology. Therefore, in this study, all water tables were related to the mean surface elevation (μ) of zero for the area around the monitoring well.

Frequency distributions of surface elevations were similar and “normal-like” around all three monitoring wells. Table 2 gives the standard deviations (σ) of the fitted normal distributions for the south, central, and north monitoring wells.

The predefined criteria were fulfilled for 145 events for the south monitoring well, 113 events for the central well, and 114 events for the north well. Differences in the number of detected events among the monitoring wells were caused by differences in the values of Δz, zu, and the mean slope of the water table rise.

Figure 3a illustrates, as in Fig. 2 also for the south monitoring well, the water table rises from zl to zu vs. S_y, calculated from the detected events (S_y = Δz/P; Logsdon et al., 2010). Most of the detected events were close to the surface. As seen in Fig. 3a, S_y values increased for near-surface water table rises. This clearly indicates the contribution of the surface storage for shallow water tables and demonstrates the necessity of accounting for the microrelief when predicting shallow water table changes. For lower water tables, the influence of the WRC of the soil increased.

Figure 3a further indicates that a macroporous layer (with high S_y values) below the measured mean surface elevation contributed to the surface storage (see also above). This can be seen by the S_y values for the events with water tables between approximately −0.2 and −0.1 m. Applying the normally distributed surface elevations of the microrelief (σ = 0.071 m, μ = 0) at the south monitoring well to a water table change from −0.2 to −0.1 m results into a S_y,surface value of 0.02504 (after division by Δz). Thus, the S_y values shown in Fig. 3a have an offset to the expected behavior of about 0.10 m. This offset can be explained by a topsoil layer that acts like a porous system that releases water in the range of the occurring matrix potential fluctuations. How deep the top surface soil acts like surface storage is shown by the optimization of Δμ as detailed above. The water table rises from zl to zu vs. S_y at the central monitoring well (not shown) showed a similar behavior as at the south monitoring well. At the north monitoring well (not shown), the water table rises vs. S_y indicated a lower contribution of the top surface soil to the surface storage.

For all three investigated monitoring wells, the 95% confidence intervals (CI₉₅) of the simulated water table changes, obtained from the pdfs, matched the observed ones fairly accurately for both deep and shallow water tables (Fig. 3b; again only shown as an example for the south monitoring well).

The pdfs of the estimated vG parameters and Δμ are depicted in Fig. 4 for the south and central monitoring wells (horizontal axes for parameters α, n, and Δμ are narrowed from the uniform prior). The corresponding CI₉₅ quantiles of the optimized parameters n, α, and Δμ are given in Table 3 for all three monitoring wells.

In proportion to the prior distribution of n (1.01–30), α (0.001–0.5 cm⁻¹), and Δμ (−20 to 20 cm), the pdfs of the south and central monitoring wells are well constrained. In contrast, the results for the north monitoring well demonstrate the applicability limits of our approach. For the north monitoring well, the inversion failed to constrain parameter n at the lower boundary of the prior distribution, and the pdf of parameter α is wide compared with the south and central monitoring wells. This leads to unrealistic soil WRC within the CI₉₅. We interpret this as a failure of our approach for the north monitoring well. Water tables at the north monitoring well are approximately 0.1 m shallower than at the other monitoring wells. For the south and central monitoring wells, S_y values increased for events with zl < −0.3 m, reaching S_y values of approximately 0.2 (indicated in Fig. 3a for the south monitoring well). The increasing S_y values at comparatively low water tables indicate that the soil substantially releases water at equivalent pressure heads in the upper soil. For constraining the vG parameters α and n, this information is crucial because these parameters define at which pressure heads the soil dewaters substantially and how abruptly this occurs. The lowest water table in the time series for the north monitoring well was −0.32 m. Below −0.3 m, three events were detected with S_y values between 0.07 and 0.08. However, these S_y values are not increased compared with events at shallower water tables. Following this, at the north monitoring well, the deepest water tables were not sufficient, i.e., the extent of the unsaturated zone was too limited, so that a substantial dewatering of the topsoil did not take place. Therefore, the vG parameters n and α could not be constrained.

For all monitoring wells, the inversion failed to constrain θ_s − θ_r. The information for the lowest water tables was not sufficient to constrain θ_s − θ_r.

The uncertainty bounds (CI₉₅) of the optimized soil WRC are visualized in Fig. 5. Also depicted is one of the soil WRC derived by Dettmann et al. (2014) near the central monitoring well. The figure shows a steeper decline of the soil WRC for decreasing pressure heads (increasing suctions) for the south monitoring well than for the central monitoring well. This indicates that the peat at the south monitoring well has a narrower pore size distribution (higher vG parameter n) and thus dewaters more abruptly. This may be an effect of different peat degradation and plant composition of the peat substrates, with pure sphagnum at the south monitoring well and the combination of sphagnum, sedges, heather meadows, and bog pines at the central monitoring well. According to Bartels and Kuntze (1968), pristine sphagnum peat is characterized by a more abrupt dewatering characteristic than less pristine sphagnum peat. In former centuries, the central monitoring well was probably partly affected by nearby ditch drainage. Hence the peat around the south monitoring well was expected to be more pristine, although this was not determined directly at the site. Following this, the steeper soil WRC at the south monitoring well (shown in Fig. 5) may be related to the more pristine peat than at the central monitoring well. The pressure heads at which the soils start to dewater substantially are in the same range for both monitoring wells, which is indicated by only slight differences in the characteristic pore size (1/α).

There are also methodological errors that may have contributed to the differences in parameters between the monitoring wells. They are discussed below.

Dettmann et al. (2014) studied the WRC of the sphagnum peat at the study site near the central monitoring well with evaporation experiments and inverse optimization. Laboratory values from their study (α = 0.0456 cm⁻¹, n = 1.72) indicate an earlier (at lower suctions) dewatering compared with the WRC determined in situ (Fig. 5). This difference between the WRC indicates that WRC determined in the laboratory may not be transferrable to the field. To investigate the transferability, we predicted water table changes with Eq. [7] by applying the laboratory WRC and the surface storage parameters, i.e., σ from the survey data and Δμ of the optimization (median of the pdf). Figure 6 shows the observed water table changes for the central monitoring well plotted vs. the predicted values. The figure indicates a bias of the prediction toward underestimated water table changes.

Besides methodological errors (discussed below) that may influence the parameters derived in situ, it can be speculated about different hydraulic behavior of the peat under laboratory and field conditions. As mentioned above, shrinkage and swelling of peat might be more distinctive under laboratory conditions. This affects the θ–h relation derived in the laboratory. Note that we defined the reference volume to calculate the volumetric moisture, θ, to be constant (cracks still belong to the total peat volume). It is known that shrinkage, as a result of increasing soil water suction, partly releases water stresses (Baumgartl and Köck, 2004), i.e., in the laboratory more water is effectively released at low suctions. In other words, with the laboratory WRC, more water can be stored in the unsaturated peat profile when it rains, i.e., the predicted water table changes underestimate the observed changes.

A second reason for the underestimation of the water table changes might be a scale effect. In the laboratory, the soil samples had a height of 0.2 m and width of 0.3 m. They were taken from the upper part of the peat profile, just below the vegetation layer. In contrast, the in situ WRC represents a spatial integration around the monitoring well while assuming a uniform soil profile with a vertical extent of about 0.65 m, which is the distance between the highest surface elevation and the lowest water table of the time series. As for rather pristine peatlands, where porosity and pore sizes commonly decrease with depth, the observed water table changes are higher than the ones predicted from laboratory parameters that were only based on samples of the upper 0.2 m of the peat profile.

The CI₉₅ of the derived pdfs of Δμ (Fig. 4) are between −8.9 and −10.4 cm for the south monitoring well, between −10.4 and −12.7 cm for the central monitoring well, and between −2.2 and −0.1 cm for the north monitoring well (also listed in Table 3). Regarding the south monitoring well, the pdf of Δμ confirms the visual assumption from Fig. 3a that approximately the top 0.1 m of the soil contains a high amount of macropores, i.e., presents high S_y values and thus is better described by ΔA_surface. The same effect was observed for the central monitoring well.

Figure 3a shows that there is considerable scatter in the depth dependency of S_y obtained from the analyzed rain events and related water table rises. This scatter can be attributed to several potential error sources, which can be separated into precipitation- and soil-related errors and can be both statistical and systematic.

About 35% of the P events were measured with a gauge within the bog, 30 to 325 m away from the three monitoring wells. Given the proximity, the error due to spatial variability of P can be considered low for these data. However, data gaps were filled with data from the German weather forecast station Eberfing, which is 10 km away. The spatial variability of P can be high for such distances. Further, the distance leads to a temporal offset of the rainfall. For the gap-filled data, thus, we expect higher errors for the total P sums of the events. Comparing the P sums of the two rain gauges for the overlap period did not show any systematic difference. The P sums were either higher or lower for the same event. This can be an effect of different wind directions, wind speeds, or other climate conditions. Given the short overlap period, a correction of the error by the distant rain gauge was not possible.

Measurements with P gauges are adversely affected by undercatch bias, induced by wind turbulence, evaporative losses, wetting losses, or mechanical plugging (Groisman and Legates, 1994). The P undercatch is highly variable depending on the wind characteristic (Neff, 1977). On average, measured P values are lower than the actual P. In contrast, measured P values are higher than the amount of water that actually infiltrated due to interception in the vegetation layer. Regarding Δz, several model assumptions for the soil processes are not entirely fulfilled in reality. Macropore flow, lateral flow, shrinkage and swelling, hydraulic gradients within the soil profile before and after P, and water repellency effects are variable in their presence and intensity among the different events and thus lead to statistical and systematic errors.

The pore size distribution of the peat soil matrix and surface elevations change due to shrinkage and swelling. Therefore, hydraulic properties can be different for different P events. This can affect the corresponding Δz and contributes to the scatter in Fig. 3.

During macropore flow, water is bypassing the soil matrix (Dekker and Ritsema, 1996). At our site, such flow might be channeled through root channels and along coarse peat-forming plant material. As result, water tables rise abruptly, with a successive decline after the maximum when water redistributes into the soil matrix. We observed immediate and steep water table declines directly after the water table maximum for a variety of events for both deep and shallow water table rises but not for all. When water tables showed such a decline, the length and steepness of the decline was very variable. This shows that macropore flow is variable and the intensity may depend on several conditions, such as, e.g., the P intensity or water repellency condition of the sphagnum moss and peat before P. Water repellency is a time-dependent physical property of the soil (Dekker and Ritsema, 1996). It depends on the initial moisture at an event and can induce macropore flow. In dependence on the degree of decomposition and plant origin, sphagnum mosses and peats have been reported to be either water repellent (Michel et al., 2001) or hydrophilic (Kettridge et al., 2014). For the investigated study site, no water repellency measurements have been done. It is possible that water repellency conditions occur in the top centimeters for deep water tables during dry periods. Because the decline behavior after water table rises was very variable and for longer time periods and lateral water fluxes may become important, we decided not to correct Δz for the decline. Thus, if macropore flow and water repellency are relevant for our site, the measured Δz values are systematically higher than the prediction of a model that ignores these processes.

Further, deviations from the hydrostatic equilibrium assumption within the soil profile before and after P and lateral fluxes during the event may generate further uncertainty. If drying-induced deviations from hydrostatic equilibrium within the soil profile occurred due to an evapotranspiration period before the event, parts of P would wet up the soil instead of causing Δz. Deviations from hydrostatic equilibrium have been observed for peatlands, e.g., by Lukenbach et al. (2016), during extended periods without rainfall and for water levels lower than −0.6 m in the Utikuma Lake Research Study Area of the Boreal Plain ecozone. However, at our site we expected the soil moisture profile to be close to hydrostatic equilibrium before and after P events, as pressure observations from unpublished tensiometer data at −10-, −20-, and −30-cm depths showed for this peat soil. For the top vegetation–soil layer, however, deviations from hydrostatic equilibrium may occur due to evapotranspiration. Furthermore, it should be noted that in peatlands with less humid climates and deeper water tables, deviations from hydrostatic equilibrium may be higher for the whole profile.

A water storage change due to the local difference in lateral in- and outflow (i.e., local flow divergence) represents another error source. In particular, for shallower water tables, overland flow directed out of the bog system may occur during the time span of the event and generate significant local differences of lateral in- and outflows within the peatland. We emphasize that we do not refer here to the magnitude of the lateral fluxes, which can be very high during precipitation events, but to the local differences in them that generate water storage changes. These lateral flux differences are neglected in our approach for the time period of an event. Both vertical hydraulic gradients during evapotranspiration periods and lateral fluxes lead to measured Δz values that are systematically lower than the prediction of a model that neglects these influences. It can be noted that for both input variables, P and Δz, there are statistical errors and opposing systematic errors. Because sufficient information about the error sources is not available for a reliable quantification, our strategy in this study was to not attempt to correct for any bias and to apply the data as they are.

Above, we discussed the differences between the obtained soil WRC at the south and central monitoring wells and how they could be explained by field variability. From the uncertainty discussion here, it becomes obvious that also methodological errors that deviate among the monitoring wells may have contributed to these differences. The observed water table rises at the south and central monitoring wells showed different Δz for the same P events. For most events, Δz was higher at the south monitoring well than at the central monitoring well. This led to different n and Δμ values at the south monitoring well.

Most of the error sources above may have led to higher Δz at the south monitoring well compared with the central monitoring well. The most evident influence may have had the different compositions of the vegetation at the two monitoring wells. At the central monitoring well, a combination of sphagnum mosses, sedges, heather meadows, and bog pines is present, while the south monitoring well is surrounded by sphagnum mosses only. We expect substantially higher interception for the higher vegetation at the central monitoring well, which would directly decrease Δz because less water would have actually infiltrated.

Further, one could argue that water repellency and macropore flow are more relevant at the south compared with the central monitoring well. The lower vegetation at the south monitoring well may have led to an enhanced drying of the uppermost soil–vegetation layer, which may have provoked flow channeling due to water repellency and thus macropore flow. In contrast, higher vegetation at the central monitoring well additionally decreased the rainfall intensity due to the retardation effect of the canopy, which further increased the potential macropore flow at the south compared with the central monitoring wells.

Summarizing, both statistical and systematic uncertainties influence the shape and scatter of the observed depth distribution of S_y (indicated in Fig. 3a). The obtained pdfs are less constrained and potentially biased due to these uncertainties. Systematic uncertainties can have a different extent at different monitoring wells and therefore may result in apparently different soil WRC within a study site. The magnitude of these uncertainties relative to the actual field variability was beyond the scope of this study.

We have provided a new approach for the determination of the effective soil WRC, obtained in situ in the natural environment and as a spatial average around a monitoring well. The proposed approach requires only water table and P measurements, and for field sites with uneven surfaces, additional information about the surface elevation distribution of the microrelief around the monitoring well. Therefore, labor-intensive determination in the laboratory on small soil samples and the typical uncertainties of such methods (e.g., compaction or disturbance of the soil structure) can be avoided. Furthermore, the soil WRC of our approach are obtained for more representative soil volumes compared with laboratory soil samples. A disadvantage of our approach is that time series of a few years are required to obtain sufficient events that can be analyzed.

As an application example, we want to highlight the relevance of our study to improve the methodology of an evapotranspiration estimation approach. Evapotranspiration is frequently estimated by water table depth fluctuations following the method of White (1932) (Fahle and Dietrich, 2014; Loheide et al., 2005; McLaughlin and Cohen, 2014; Mould et al., 2010; Wang and Pozdniakov, 2014). The method of White (1932) uses diurnal patterns consisting of a water table decline in the daytime and a recovery phase overnight. Following the theory of White (1932) and above, evapotranspiration can be calculated directly with Eq. [4] for shallow water table changes if knowledge about the soil WRC is available. Using the soil WRC obtained by the approach of this study improves the estimation substantially. First, the soil WRC are directly estimated at the monitoring well, and the vG parameters are available as pdfs. Thus, uncertainties can be propagated to the evapotranspiration estimates. In our opinion, this is essential because evaporation estimates usually are highly uncertain. Second, surface storage plays a major role in the estimation at shallow groundwater systems. However, if the top vegetation at the soil surface contributes to the surface storage, as is shown by the optimization of Δμ, surface storage is substantially underestimated. The optimization of Δμ gives the possibility for a better estimation of the surface storage. This leads to more realistic evapotranspiration estimates.

Further, the optimization of Δμ can substantially contribute to research on the microrelief of shallow groundwater systems. For these systems, a separation into vegetation and soil layer is often difficult. The particular height of this interface cannot be defined objectively, regardless whether data are derived by surveying or by laser scanning. The optimization of Δμ gives an objective estimation about the depth of the soil–vegetation interface from a hydrological perspective.

The in situ determination of soil hydraulic properties by inverse modeling of field observations is difficult for peatlands when trying to use the full time series because there are large uncertainties in lateral flow and evapotranspiration. In this study, we thus focused only on the periods of water table rises caused by P. Our study indicated that these observations can contain sufficient information to constrain the wet range of the soil WRC around monitoring wells in peatlands. For two of the three investigated monitoring wells, the obtained soil WRC are in a plausible range and comparable to soil WRC determined in the laboratory for the site and to those of previous studies for similar sites (Bartels and Kuntze, 1968; Dettmann et al., 2014; Letts et al., 2000; Price et al., 2008). At one monitoring well, water tables were not low enough and the upper soil did not dewater enough to derive information about the soil WRC. In principle, our approach should also be applicable to other shallow groundwater systems, which however needs further field evaluation.

The major advantages of our approach over laboratory methods is the determination of soil WRC at more representative scales and under typical dynamic field conditions. Other site-specific factors, as for example shrinkage and swelling, which is typical for peatlands and leads to a significant fluctuating soil surface, are lumped together in the description of the effective soil WRC obtained from our approach.

For low pressure heads (high suctions), the indicated uncertainty was high. Because most of the unsaturated zone of peat soils remains very wet throughout the year, the effective soil moisture variations occur in the wet range of the soil water retention curve. Thus, water table depth observations do not provide the necessary information to constrain the dry range, i.e., θ_s− θ_r could not be constrained.

Besides this inevitable limitation in the dry range, uncertainty can be reduced with high-quality and continuous P data observed directly at the field site. Further, there are several potential systematic errors due to simplified model assumptions. Given the lack of knowledge about detailed peat properties, their relevance can hardly be evaluated. Non-equilibrium phenomena such as hydrophobicity and macropore flow may be investigated by the water table declines after the maximum of the events in combination with soil moisture probes. More knowledge about such processes could help to improve the method developed in this study. A quantitative consideration was beyond the scope of this study and should be a topic of future research.

Finally, we emphasize that our approach did not account for vertical soil heterogeneity but has been used to estimate effective parameters assuming a uniform soil profile. In particular for peat soils, vertical heterogeneity often occurs and can be high (Morris et al., 2015), and functions describing typical vertical trends could be used to advance our approach.

We acknowledge the financial support by the joint research project Organic Soils funded by the Thünen Institute. We appreciate support from Bärbel Tiemeyer (Thünen Institute, Braunschweig), Enrico Frahm (PTB, Braunschweig), Eike Maurer, Annette Freibauer (LFL, Freising), Benedikt Scharnagl (UFZ, Halle), Marika Bernrieder, and Janina Hommeltenberg (KIT, Garmisch-Partenkirchen). We further thank the District Government of Upper Bavaria for making access to the site possible. Special thanks to Jan Vanderborght as associate editor for his excellent comments and suggestions to improve the manuscript. We are grateful to Martine van der Ploeg for her review that included several helpful comments and suggestions, to a second anonymous reviewer, and to Insa Neuweiler as co-editor.