Abstract

Soil moisture is the essential control of water and energy dynamics at arable sites. Time series of soil moisture reflect the interplay of various processes, each of which influences the overall soil moisture dynamics. In this study we tested an approach to break down observed soil moisture behavior into the respective contributions of individual processes. We applied a principal component analysis to soil moisture time series from a field experiment comprising two crop rotation systems and two different soil tillage practices. We concentrated on 57 soil moisture time series measured over nearly 4 yr at 12 plots and five soil depths, down to 1.5 m. About 77.9% of the variance was reflected by the first component being almost identical to a time series of averaged soil moisture. It described the effect of the meteorological boundary conditions. The second component described the effect of the input signal damping increasing with soil depth and accounted for 7.8% of total variance. The signal transformation over depth proved to be more or less uniform throughout the test site, despite considerable soil heterogeneity. Another 3.6% of the total variance (third component) was unambiguously explained by the different cropping systems. On the contrary, different soil tillage practices had no significant effect. The suggested approach opens up many possibilities to analyze and better understand complex soil system behavior. The data-based approach of time series analysis provides model-independent, quantitative information about the key factors and processes controlling soil water dynamics. Hence, it is especially valuable for model building, calibration, and evaluation.

Soil moisture θ [L3 L−3] is probably the most important measurable status variable to evaluate the water and energy budgets of arable sites (Vereecken et al., 2008). However, soil moisture usually exhibits considerable spatial heterogeneity (Corradini, 2014; Vereecken et al., 2014), making all means of data evaluation more complicated. This heterogeneity can be caused by natural effects, such as the spatial variability of soil texture (e.g., Hohenbrink and Lischeid, 2014; Jawson and Niemann, 2007; Schlüter et al., 2012), as well as by different land use practices, such as soil tillage techniques (e.g., Perfect and Caron, 2002; Schwen et al., 2011). The impacts of land use practices on soil moisture dynamics must be quantified to evaluate the efficacy of specific management decisions. Interactions between different land management effects are usually investigated using simulation models (an overview of which is provided by Boote et al., 2013). However, this always requires a priori knowledge or assumptions about the governing processes. Therefore, model-independent diagnostic tools to disentangle the different sources of soil moisture heterogeneity are urgently needed to analyze existing monitoring data sets. Such tools could even be used to reduce uncertainty in numerical models.

There are several approaches to determine the patterns hidden in moisture data sets (for an overview see Vereecken et al., 2014), each emphasizing different features of soil moisture patterns. Classical geostatistical methods are used to characterize the spatial distribution of soil moisture. The spatial covariance of soil moisture from different locations is usually expressed as a function of distance, illustrated by variograms (e.g., Baroni et al., 2013; Joshi and Mohanty, 2010; Korres et al., 2015). More recently, more sophisticated methods, such as wavelet analysis (e.g., Biswas, 2014; Peng et al., 2013; Rivera et al., 2014), self-organizing maps (e.g., Zou et al., 2012), or fractal analysis (Korres et al., 2015), were used to investigate the spatial distribution of soil moisture. Temporal stability analyses of soil moisture patterns (Vachaud et al., 1985) are designed to identify locations with soil moisture dynamics representative of the whole observation site. This approach is based on the observation that the ranks of soil moisture values from different locations remain almost constant over time (see the overview by Vanderlinden et al., 2012). Standard methods of time series analysis, such as autocorrelation and cross-correlation analyses, have been used to compare soil moisture dynamics from various depths (De Lannoy et al., 2006; Hohenbrink and Lischeid, 2015; Mahmood and Hubbard, 2007) and relate them to other variables such as precipitation and temperature (Mahmood et al., 2012). Low-pass filtering behavior (i.e., changing periodicity of soil moisture oscillations) was investigated for vadose zones (Hohenbrink and Lischeid, 2015; Katul et al., 2007) and entire catchments (Gall et al., 2013). Slopes of the power-spectra were used to characterize and compare the periodicity distributions of time series from different locations. Principal component analysis (PCA), also known as empirical orthogonal functions analysis (see Vereecken et al., 2014), was frequently used to decompose the variance in soil moisture data sets from arable fields into uncorrelated patterns and relate them to specific explanatory variables (e.g., Baroni et al., 2013; Korres et al., 2010; Qiu et al., 2014).

Applying the presented methods to monitoring data sets from various investigation sites revealed that soil moisture patterns can be controlled by topography (e.g., Qiu et al., 2014; Yoo and Kim, 2004), soil textural properties (e.g., Baroni et al., 2013; Jawson and Niemann, 2007; Yoo and Kim, 2004), soil organic carbon content (Korres et al., 2010), vegetation (e.g., Baroni et al., 2013; Korres et al., 2015), land management (e.g., Korres et al., 2010, 2015), and meteorological conditions (e.g., Joshi and Mohanty, 2010; Qiu et al., 2014). However, most researchers only investigated spatiotemporal moisture patterns from a shallow soil depth on a few selected dates. Hupet and Vanclooster (2002) stressed the importance of using measurements from the entire hydrologically active zone to investigate soil moisture spatial variability. However, few studies analyze long-term soil moisture dynamics in vertical profiles deeper than 1 m (De Lannoy et al., 2006; Hu and Si, 2014). Hohenbrink and Lischeid (2015) proposed to apply a PCA on soil moisture time series measured in vertical soil profiles. They investigated the effects of textural heterogeneity on the transformation of hydrological signals (e.g., rainfall, snow melt) propagating through the vadose zone. Furthermore, they stated that the approach has great potential to distinguish independent effects contributing to hydrological behavior observed in soil systems. Their approach was introduced with the example of simulated time series resulting from a numerical experiment. However, the transferability of their results from a “simulated world” to the “real world” still has to be proven based on soil moisture time series measured under field conditions. In fact, that approach has already been successfully applied to groundwater head time series (Böttcher et al., 2014; Lehr et al., 2015; Lischeid et al., 2010; Page et al., 2012).

Our objectives were to identify and describe the specific effects of soil heterogeneity and land management practices on soil moisture dynamics. To achieve this, we decomposed the total variance of measured soil moisture time series by PCA as suggested by Hohenbrink and Lischeid (2015). An adequate monitoring data set of soil moisture time series was required to apply the approach. We used a monitoring data set from a multifactorial field experiment in northeastern Germany as it had the following features: (i) intensive instrumentation, (ii) large number of replicates, (iii) different well-documented soil treatments, (iv) homogeneous precipitation inputs due to the small size of the test site, and (v) large soil heterogeneity spatially uncorrelated with the experimental set-up. In this paper the term hydrological signal designates spatiotemporal changes in the soil moisture which are propagated through the vadose zone (see Hohenbrink and Lischeid, 2015). We use the term functional heterogeneity to express the variability among measured time series (see Basu et al., 2010; Hohenbrink and Lischeid, 2015) bearing information about the hydrological system behavior of the soils investigated.

Methods

Field Experiment

The experimental test site (52°31′01″N, 14° 07′27″E, 62 m a.s.l.) was situated in the lowlands of northeastern Germany close to the city of Müncheberg. It is located in the transition area between a maritime climate and a continental climate. Between 1 Jan. 1981 and 1 Jan. 2012, the annual sums of precipitation and potential evaporation were 529 and 659 mm, and the mean annual temperature was 9.1°C (DWD, 2014). The experimental site is covered by a layer of sediment from the Late Pleistocene, with a flat surface. Sand layers alternate at short distances with glacial till containing clay and marl, resulting in large textural heterogeneity. More information about the depth distribution of soil texture at the experimental site can be found in a data set published by Mirschel et al. (2010). The predominant FAO soil type at the test site is Orthic Luvisol.

Twelve experimental plots with an area of 788 m2 each were arranged on an experimental field (Fig. 1). Six plots were managed using a cultivation system designed to stabilize the organic carbon content of the soil. This featured a 4-yr crop rotation system of winter rye (Secale cereale L.), forage sorghum [Sorghum bicolor (L.) Moench], winter triticale (×Triticosecale Wittm. ex A. Camus [Secale × Triticum]), a mixture of alfalfa, clover and grass (Medicago sativa L., Trifolium pretense L., and Lolium perenne L.), and maize (Zea mays L.) each harvested as whole plants. This system is hereafter referred to as “CropRo4.” The other six plots were cultivated using a 1-yr crop rotation system of maize and winter rye designed to maximize the yield of biomass (hereafter referred to as “CropRo1”). The maize was harvested as a whole plant. The winter rye was only grown for erosion protection in the winter (as a cover crop) and was removed at an early stage of growth. Thus, its soil water consumption was marginal. Half of the plots were managed by plowing and the others were not tilled (direct seeding). Thus, each combination of cultivation system and tillage practices was performed in three replicates. Each experimental plot was equipped with seven FDR soil moisture sensors (ThetaProbe ML2x, Delta-T Devices) at depths of 0.3, 0.6, 0.9, 1.2, 1.5, 2.0, and 3.0 m. All devices and cables were installed at least 30 cm below the ground to ensure that the upper 25 cm of the soil could be plowed. The soil moisture was measured for nearly 4 yr between 1 May 2008 and 23 Apr. 2012 with an hourly time resolution and aggregated to daily data.

Data Preparation

The measured time series were initially plotted in a diagram and checked visually for plausibility. They cannot be shown here, due to their large number. Time periods with frozen soil water had to be omitted from the data set, because the FDR probes were calibrated for non-frozen soils only. Thus, an exclusion criterion was defined on the basis of soil temperature time series measured by DWD (2014). Soil moisture values were not considered when the soil temperature at either 0.2 or 0.5 m was smaller than 1°C, yielding four data gaps during the winter. Two additional time periods of 30 d (from 20 May 2009 to 18 June 2009) and 41 d (from 14 July 2011 to 23 Aug. 2011) were omitted due to data gaps in several time series caused by malfunctioning data loggers. The time series measured at 150 cm at Plots 03, 06, and 10 and at 200 cm at Plot 11 had to be omitted totally because those sensors malfunctioned. All measurement gaps in the remaining time series were smaller than 6 d. An autocorrelation analysis of all time series showed that the minimum autocorrelation for time lags of 6 d was 0.80. This means that missing values can be predicted by existing with an accuracy of at least 80%. Thus, all gaps smaller than 6 d were interpolated linearly. The fraction of interpolated data points considered in the analysis was 1.24%.

Each individual moisture time series θ(t) was scaled to zero mean and unit variance (z transformation) 
graphic

Information about the absolute values and amplitudes of θ (t) are lost during z transformation, but it makes the scaled moisture time series θ*(t) comparable and they are weighed equally. In the first step, only data from the upper five depths (z ≤ 150cm) were considered for detailed analyses and discussion, because plant effects were presumably more significant here compared to greater depths. The z-transformed moisture series exhibited annual fluctuations with maxima in winter and minima in summer (Fig. 2). As the soil depth increased, the time series became smoother and lagged in time. However, time lags varied between different time periods. In Fig. 2 time shifts between the shallow and great depths are highlighted for two local minima. Time lags of soil moisture time series were much larger in the autumn of 2008 (marked with an “a”) compared to the summer of 2010 (marked with a “b”), although they were measured at identical positions. All z-transformed time series were organized in an input matrix Θ where each column represented a time series and each row a date.

Principal Component Analysis of Time Series

A PCA was applied to the input matrix of time series. This is an ordination method to analyze the structure of a multivariate data set and to identify common temporal patterns among time series. A PCA finds a new orthonormal basis for the multivariate data space spanned by all input time series. The total variance of the input data is decomposed into independent fractions. This results a set of uncorrelated principal components (PCs). Here, the target is to identify a first component which explains most of the total variance, a second component which explains most of the remaining variance and so on. All the resulting PCs are summarized in a matrix P in the order of the fraction of variance they explain. The transformation of the input matrix Θ into P can be computed by 
graphic

The matrix Λ contains the eigenvectors of the correlation matrix of Θ ordered by their corresponding eigenvalues. The fraction of total variance explained by a PC can be calculated by dividing the corresponding eigenvalue by the sum of all eigenvalues. Often only few PCs are needed to explain much of the total variance. We refer to Jolliffe (2002) for a detailed mathematical description of PCA. One important limitation of PCA is that only linear patterns are detected. Nonlinear patterns in a data set can only be approximated by gradually superposing several linear PCs (Lee and Verleysen, 2007). The input time series should ideally be normally distributed, but this condition is less important when PCA is considered a mainly descriptive technique (Jolliffe, 2002).

The PCs represent time series describing uncorrelated temporal patterns that presumably reflect different effects on the observed dynamics. Pearson correlation coefficients among PCs and the measured time series are known as loadings, L. The fraction of variance forumla of individual time series explained by the m first components can be calculated by the loadings 
graphic

We calculated the arithmetic mean values forumla of the explained variances forumla of all individual time series from each measuring plot. The means of explained variances forumla were used to quantify the prevalence of temporal patterns described by single or several PCs at the individual measuring plots. On the basis of forumla, forumla, L, and additional information (e.g., the measuring depth of the moisture time series), the temporal patterns described by individual components can be related to specific factors causing the observed patterns. We systematically tested whether all PCs were associated with effects of the following factors (significance level: 0.01):

  1. Mean soil moisture behavior: The Pearson correlation was calculated between the daily arithmetic mean values of all input time series and the scores of each PC. A t test was performed to see whether the resulting correlation coefficients r differed from 0.

  2. Soil depth: The Pearson correlation was calculated between the installation depth of single sensors and the loadings of the corresponding time series for each PC. A t test was performed to see whether the resulting correlation coefficients r differed from 0.

  3. Cropping system: A Wilcoxon–Mann–Whitney test was performed for each PC to see whether the loadings from the CropRo1 and CropRo4 plots differ.

  4. Soil tillage: A Wilcoxon–Mann–Whitney test was performed for each PC to see whether the loadings from the tilled and nontilled plots differed.

In cases where the first component describes the mean course of all the considered soil moisture time series and another component represents an effect of soil depth, Hohenbrink and Lischeid (2015) suggested evaluating the input signal transformation behavior of the vadose zone on the basis of both components. A signal damping coefficient D quantifying the extent of smoothing and time lagging in each moisture series can be derived from the ratio of the loadings L1 and L2 as described by Hohenbrink and Lischeid (2015) and Lischeid et al. (2010). The damping coefficient D represents a dimensionless relative measure that can only be interpreted in the context of the data set evaluated by the PCA. Temporal dynamics explained by components of a higher order are neglected by D because it is only based on the first two components. In this study all calculations and statistical analyses were performed with the R software package (R Development Core Team, 2010).

Results

Mean Behavior of Soil Moisture and Effect of Soil Depth

First Principal Component

The first PC explained 77.9% of the data set’s total variance. The mean values of explained variances ranged from forumla = 0.653 at Plot 09 to forumla = 0.841 at Plot 06 (Table 1). All loadings were positive (Fig. 3a). The absolute values of the loadings were between 0.697 and 0.960. Thus, the first PC described a temporal pattern that was apparent in every single time series. At each plot, the largest loadings appeared at medium depths, showing that the pattern described by the first PC was most predominant in the center of the soil profiles. The first PC was nearly identical (r2 > 0.999, p value ≪ 0.01) to a time series representing the daily arithmetic mean values of all input time series. Thus, the first PC described the mean behavior of soil moisture, which was strongly influenced by atmospheric controls.

Second Principal Component

The fraction of total variance explained by the second PC was 7.8%. The loadings of the time series on the second PC varied between −0.446 and 0.439 and were correlated (r2 = 0.884, p value ≪ 0.01) with soil depth. Loadings increased with the depth, from negative values in the topsoil to positive values at greater depths (Fig. 3b). Thus, the second PC reflected deviations from the mean behavior (first component) depending on the soil depth. Combining the first two PCs by linear combination according to the principle of superposition resulted in reconstructed time series. They represented approximations of moisture dynamics that could be expected if soil depth were the only influencing factor. Different desired soil depths were specified using the ratio of combination. Adding the second PC to the first one yielded a time series (blue line in Fig. 4) that was more strongly damped (smoother) and time lagged compared to the mean behavior (black line). Temporal patterns of this kind were observed at large soil depths, where loadings of the second PC were positive. On the contrary, subtracting the second PC from the first one yielded a weaker damped time series (orange line) with temporally preceding dynamics compared to the mean behavior. Such moisture time series were observed at shallow depths and loaded positively on the second PC. The variances of the first and second component were combined in the ratio 70:30 to construct the moisture dynamics shown in Fig. 4. Following Eq. [3] this corresponds to the loadings L1 = forumla = 0.83 and L2 = forumla = ±0.54. Soil depth had a slightly stronger effect in the constructed time series than in the observed ones, since the absolute values of L2 were slightly larger than those of the time series measured at shallow and great depths (Fig. 3).

Time shifts between the moisture dynamics reconstructed for shallow and great soil depths varied between different time periods similarly to the measured moisture series. In Fig. 4, different time lags of local minima in the autumn of 2008 and the summer of 2010 are highlighted with “a” and “b” as in Fig. 2.

Signal Damping

It was possible to calculate signal damping coefficients D quantifying the damping status of each measured time series, since the loadings on the second PC were correlated with the soil depth. The calculation of D was based on the loadings of the first two components only. Thus, the determination of D made use of 85.7% of the data set’s total variance. The mean value of variances explained by the first two PCs used to calculate D varied between forumla = 0.743 at Plot 09 and forumla = 0.910 at Plot 02 (Table 1). At each plot D increased along with the soil depth; in other words, time series measured at a greater soil depth were more strongly damped compared to those from shallow depths (Fig. 5). The relation between D and soil depth was almost linear. A regression line (intersect: −84.5 cm; slope: −129.4 cm per unit of D; r2 = 0.88) described the field-averaged signal damping status of the measured time series. Its slope indicated the averaged signal transformation intensity and thus characterized the signal transformation behavior of the vadose zone. Moisture dynamics showing the averaged damping state can be reconstructed for any soil depth by combining the first two components. That way, the measured soil moisture time series can be functionally averaged. This means averaging the signal transformation behavior of the vadose zone at the experimental site (see Hohenbrink and Lischeid, 2015). The damping coefficients of individual time series differed only slightly from the field average, showing that hydrological signals were similarly transformed at all plots. Mean absolute errors (MAEs) and correlation coefficients (R2) between individual and averaged damping profiles were used to evaluate deviations from the averaged behavior. Plot 12 (MAE = 0.031, R2 = 0.989) and Plot 03 (MAE = 0.040, R2 = 0.994) behaved most similarly to the average while the largest deviations from the mean behavior occurred at Plot 10 (MAE = 0.147, R2 = 0.990) and Plot 05 (MAE = 0.132, R2 = 0.815). All time series from Plot 10 were more strongly damped than the field average. Hence, the damping profile ran parallel to the averaged one. Note that the mean values of explained variances forumla were correlated neither with MAE (r2 < 0.005) nor with R2 (r2 < 0.04). This shows that moisture dynamics deviating from the field-averaged damping status can be approximated by the first two components to the same extent as those representing the average behavior.

Effect of Cropping System

The third PC covered 3.6% of the data set’s total variance. In contrast to the first two components, the mean value of explained variances (Table 1) was largest at Plot 09 (forumla = 0.094) and smallest at Plot 02 (forumla = 0.009). The loadings were clustered in two groups (Wilcoxon–Mann–Whitney Test, p value ≪ 0.01). The time series measured at the CropRo1 plots were positively correlated with the third PC (Fig. 6), and those measured at the CropRo4 plots were negatively correlated. Consequently, the third PC discriminated between the cropping systems with respect to differences in the soil moisture dynamics between the CropRo4 and the CropRo1 plots (Fig. 7). Positive scores of the third PC indicate that the z-transformed soil moisture values were greater at the CropRo1 plots than at the CropRo4 plots. The opposite was true during time periods with negative scores. Color bars at the top and bottom of Fig. 7 indicate the plants grown at the CropRo1 (top) and CropRo4 (bottom) plots during different time periods. In May 2008 the rye plants at the CropRo4 plots consumed more water than the very young maize plants at the CropRo1 plots. Hence, soil moisture at the CropRo4 plots decreased compared to the CropRo1 plots, as indicated by increasing scores of this component. This development promptly reversed after the rye harvest (marked with an “a” in Fig. 7). The scores decreased until the maize was harvested at the CropRo1 plots (marked with a “b”). During this time period, the maize plants built up much of their biomass, consuming more water than the millet plants at the CropRo4 plots. The same interplay of plant growth and water consumption was observed in the following year, 2009 (turning points marked “c” and “d”). The course of the third PC also showed a response to the alfalfa–clover–grass mixture being mowed (marked “e” and “f”). The scores increased before the mowing dates, showing that water consumption was higher at the CropRo4 plots. After mowing, the water consumption of the alfalfa–clover–grass mixture was reduced, leading to decreasing scores. This effect was not clearly visible when the mowing date was shortly before or after maize harvesting at the CropRo1 plots (marked “d” and “g”). In the last observed season in 2011 and 2012, all plots were uniformly cultivated with maize followed by rye in the winter. During this time span, the scores of the third component were close to zero. Thus, there were no major differences in soil moisture dynamics between the CropRo4 and CropRo1 plots. It is worth mentioning that the third component described relative soil moisture dynamics instead of providing information about absolute soil moisture values, since each input time series was z transformed (Eq. [1]).

The first three PCs can be combined to reconstruct time series approximating the soil moisture dynamics that could be expected if only the soil depth and cropping system affected soil moisture. Such reconstructed time series aggregate the most pronounced effects of soil depth and the cropping system. They can be used both for (i) detailed analyses of the specific effects and (ii) upscaling purposes. The time series shown in Fig. 8 were composed by superposing the first three PCs. The third PC was considered both positively (CropRo1) and negatively (CropRo4). The ratios for combination were determined by their fraction of explained total variance, where the sum of explained total variance of 89.2% was set to 100%. Thus, the third PC accounted for 4.03%. The remaining 95.97% were then assigned to the first (95.80%) and second PC (0.17%) in a ratio determined by the field-averaged value of D for a depth of 90 cm. Note that the second PC only accounted for 0.17% because the pattern represented by this PC was almost nonexistent in the middle of the soil profile (Fig. 3b). The time series show reconstructed soil moisture dynamics that can be expected for both cropping systems at a soil depth of 90 cm. The temporal patterns emerging from mean behavior, soil depth, and cropping system were functionally averaged for the whole experimental field. The moisture dynamics are controlled by precipitation and plant transpiration. The most significant differences between the cropping systems appeared in the summer months when transpiration rates were highest and soil moisture deficits were not rapidly compensated by precipitation.

Effect of Soil Tillage

The design of the field experiment included two soil tillage options (Fig. 1). The loadings of the time series from tilled and nontilled plots differed significantly for the 27th component (Wilcoxon–Mann–Whitney Test, p value ≪ 0.01), indicating that it provides information about soil tillage effects. However, it only accounted for 0.046% of the total variance. This fraction is too small for reasonable process interpretations. Nevertheless, it provides an estimation of the maximal fraction of total soil moisture variance below the plowing zone that could be caused by soil tillage effects.

Considering Data from Different Soil Depths

The previous sections showed the results from analyses of the upper five measuring depths between 30 and 150 cm. The described effects of mean behavior, signal damping, and plant water consumption were also detected using PCA when different numbers of measuring depths were considered. However, their relevance differed. The fraction of variance explained by the mean behavior decreased from 84.9% (upper three depths) to 72.0% (all seven depths) when time series from more soil depths were considered (Table 2). Conversely, the effect of signal damping was strongest (explained total variance: 9.0%) when all measuring depths were considered. The effect of water consumption of plants slightly became weaker with increasing number of considered measuring depths. The total amount of variance that could clearly be assigned to specific influencing factors decreased from 93.9% (upper three depths) to 84.0% (all seven depths).

Discussion

Decomposing the Soil Moisture Variance

Applying PCA to the z-transformed input data set yielded three meaningful components accounting for 89.3% of the total variance in the upper five soil depths. A large share of the observed temporal dynamics recurs in each time series, showing that soil moisture dynamics exhibited only a small degree of functional heterogeneity at the experimental site. This becomes particularly obvious when looking at the first component. More than three-quarters of the total variance could already be aggregated by one single temporal pattern that was almost identical to a time series of mean values from all measuring locations. These findings are in good agreement with the results of Korres et al. (2010), who analyzed spatiotemporal soil moisture patterns at shallow soil depths using empirical orthogonal functions. They also found a PC reflecting the mean soil moisture dynamics at a pasture site and an arable field. Similarly to our results, this main pattern explained more than 70% of soil moisture variance from the arable site, although Korres et al. (2010) considered only 10 measuring dates.

The explanatory power of a PCA to identify specific effects strongly depends on the information content of the analyzed data set. Thus, the experimental design already determines the applicability of a PCA. Our results have shown that the fraction of variance explained by specific effects was determined by the number of soil depths considered. However, the nature of the resulting PCs and, thus, the specific effects identified by the PCA were similar for all the numbers of soil depths considered. This underlines how robustly PCA identifies first-order controls.

A PCA is often used as a dimensionality reduction tool. In that case it is necessary to test how many components contain significant information and can thus be considered for interpretation (Peres-Neto et al., 2005). Our purpose, however, was to identify and extract the effects of specific factors known beforehand. Thus, we systematically tested whether each PC contained information about the factor of interest, notwithstanding the fraction of variance it explained.

Signal Transformation Behavior

The second component described the part of the deviation from the mean behavior (first PC) that can be explained by effects of soil depth. The fraction of total variance explained by the first two components (85.7%) was very similar to the results of Hohenbrink and Lischeid (2015). In their study both components explained 88.7% of soil moisture dynamics simulated in a heterogeneous soil profile. It is noteworthy that the measured dynamics are no more complex than the simulated ones, although the latter were generated for strongly idealized and simplified model cases. However, it is worth mentioning that Hohenbrink and Lischeid (2015) considered model cases covering nearly the whole range of possible soil textures, while the soils at the monitoring site were dominated by sandy textures.

Time series with different degrees of signal transformation were reconstructed by combining the scores of the first two components (Fig. 4) as proposed by Hohenbrink and Lischeid (2015). This was a prerequisite to derive D from the loadings of the first two components. The depth profiles of D (Fig. 5) provided evidence that hydrological input signals were continuously transformed into a more strongly damped status while propagating through the soil profiles. This might be a trivial result known before (e.g., Mahmood et al., 2012; Pan et al., 2011), but the novelty of this approach is that D represents a robust measure to quantify the signal transformation behavior of soil systems. The fact that D mostly increased linearly with the soil depth easily enables the functional averaging of soil moisture dynamics, that is, averaging the signal transformation behavior irrespective of small scale heterogeneities.

The damping coefficients of single time series only deviated slightly from the field-averaged damping profile. Furthermore, the slopes of the damping profiles from all measuring plots were very similar. Thus it can be concluded that the general signal transformation behavior of the vadose zone proceeds relatively homogeneously at the field scale. Similarly to temporal stability analyses (Vachaud et al., 1985; Vanderlinden et al., 2012), the concept of signal transformation analysis can be used to identify measuring locations that reflect the mean behavior of the whole investigation site.

At some plots, D increased linearly with the soil depth, indicating that water flow was relatively uniform. At other plots, D was scattered around the mean damping profile, as indicated by larger values of the mean absolute errors MAE and the correlation coefficients R2 between the individual and averaged damping profiles (Fig. 5). This shows that heterogeneous water flow fields occurred at plot scale; that is, hydrological signals were preferentially propagated along distinct water pathways. This effect can be caused by small-scale textural heterogeneity. In heterogeneous flow fields, the signal transformation along preferred water pathways is less intensive than in surrounding regions with a reduced water flow (Hohenbrink and Lischeid, 2015). Thus, time series measured at a certain depth can be less strongly damped than others from shallower depths. However, the prevalence of both temporal patterns captured by D was not reduced on the occurrence of heterogeneous flow fields, since forumla was neither correlated with MAE nor with the correlation coefficients R2. Soil moisture dynamics emerging in heterogeneous flow fields were still composed of the first two components to the same extent as under uniform flow conditions. Similarly to Hohenbrink and Lischeid (2015), this shows that textural heterogeneity does not necessarily cause functional heterogeneity, since no additional, more complex temporal patterns were generated by heterogeneous water flow. Soil heterogeneity merely caused different signal transformation intensities per depth unit instead of inducing completely different kinds of signal transformation. Thus, the effect of textural heterogeneity is much more systematic than, for example, random noise added to the data (Hohenbrink and Lischeid, 2015).

Effects of local heterogeneities, that is, different intensities of signal transformation, can balance each other out. This can be visualized using the example of layered soil profiles. Input signals entering any soil layer correspond to the output signal of the layer above. An input signal entering the soil surface is propagated in this way through several layers with different thicknesses, textures, and small-scale heterogeneities. Given that (i) the general kind of signal transformation does not differ between different layers and that (ii) signals can only be transformed to a more strongly damped status, it is possible to find a corresponding homogeneous soil with the same average signal transformation behavior as the layered soil. This is a basic condition for applying the concept of temporal averaging for different purposes, such as for upscaling soil moisture dynamics by averaging the signal transformation behavior of the vadose zone.

The calculation of D makes use of the first two PCs, representing linear approximations of the two most prevalent uncorrelated temporal patterns, that is, the mean soil moisture behavior and the effect of soil depth. If all depth effects occurring cannot fully be described by the second principal component, there is a risk of neglecting some parts of relevant variance that are also induced by soil depth. This typically happens if signal transformation induces nonlinear structures that can only be approximated by stepwise superposition of several linear components (Hohenbrink and Lischeid, 2015). There is further need for research into the nonlinear behavior of signal transformation. We could not find Pearson correlations between the soil depth and the loadings on higher-order PCs. However, loadings of the fourth component (data not shown) seemed to have a curve-shaped dependency on depth. Thus, more sophisticated measures for nonlinear dependence, such as mutual information—as defined by Shannon (1948), for its estimation see Fraser and Swinney (1986) or Kraskov et al. (2004)—would be needed to detect further nonlinear relations between loadings and soil depth. However, this is beyond the scope of our study.

Effects of Cropping System and Soil Tillage

The effect of the cropping systems represented by the third PC accounted for 3.6% of the total variance. This fraction was smaller than we would have expected beforehand, since both cropping systems comprised plants with different water demands, varying harvest times, and nonparallel growing seasons. The third component nevertheless described a clear and unambiguously interpretable temporal pattern. This holds true for both (i) the clustering of loadings in the two groups of CropRo1 and CropRo4 and (ii) the interpretation of the temporal dynamics described by the scores of the third component. The clarity of these results highlights the great potential of PCA in detecting minimal but significant patterns hidden in data sets of time series (Thomas et al., 2012). Soil moisture temporal patterns accounting for only small fractions of variance can still be interesting and can even have important impacts on dominant soil water processes. This can be illustrated by the example of threshold-controlled processes such as macropore flow occurring under particular moisture conditions. It can be the main reason for solute leaching (reviewed by Jarvis, 2007), although it often accounts for only small fractions of the total soil moisture variance.

Korres et al. (2015) used various analyzing tools to reveal spatiotemporal soil moisture patterns emerging under different land uses. With regard to arable sites, they found large spatial soil moisture variability among neighboring fields. They explained this effect by strongly varying evaporation rates due to shifted periods of maximum water uptake by different crops and different agricultural management (e.g., planting dates, harvesting dates, field sizes). It is worth noting that we identified very similar effects of crops on soil moisture patterns by time series analysis, such as Korres et al. (2015) found by spatial analysis. Baroni et al. (2013) performed a PCA on a data set containing soil moisture measured at a shallow depth in an agricultural field. Among other PCs they identified a component describing a plant effect. It was spatially correlated with the crop factors leaf area index and plant height. This vegetation-controlled pattern prevailed when the soil became dryer. Our third component also revealed that soil moisture differences between the cropping systems were the greatest in dry conditions, because the scores showed the highest absolute values during the dry summer months (Fig. 7).

There is agreement in the literature that in comparison to plowing, nontillage causes (i) greater bulk density, (ii) more stable soil aggregates, (iii) increased soil biota abundance generating macropores, and (iv) higher soil carbon contents (e.g., Arai et al., 2014; Holland, 2004; Soane et al., 2012). Tillage-induced alterations of these substantial soil characteristics can have various effects on soil hydraulic properties (e.g., Schwen et al., 2011; Strudley et al., 2008). Thus, it might be expected that soil moisture dynamics would also differ between different tillage practices. However, we could not detect significant differences between tilled and nontilled plots. Some authors reported tillage effects on spatial soil moisture patterns at shallow soil depths (e.g., Korres et al., 2010; Perfect and Caron, 2002). Similar effects could possibly also be found at very shallow depths at our investigation site. However, we only analyzed time series from greater soil depths, since it is not possible to measure long-term data with sensors placed at fixed positions when the surrounding soil is regularly tilled.

Conclusions

By applying a principal component analysis to measured soil moisture time series, we achieved our objectives of identifying, describing, and evaluating the specific effects of soil heterogeneity and land management on soil moisture dynamics. The method turned out to be powerful as long as relative temporal dynamics are of interest rather than absolute values. Based on the results that three meaningful components accounted for 89.2% of the total soil moisture variance, we can draw the following conclusions.

First, contrary to common assumption, the interactions of infiltration, soil heterogeneity, and different land management practices do not necessarily induce complex soil moisture dynamics in deeper soil layers. Nearly 78% of the observed soil moisture variance was identical at all measuring locations. Thus, functional heterogeneity, that is, variability among all soil moisture time series, only accounted for the remaining 22% of variance. About 35% of that was unambiguously attributed to the deterministic effect of input signal transformation with increasing soil depth. The large textural heterogeneities at the test site had no effect on the general nature of signal transformation, but did affect its intensity, which varied at different sites. Land management only slightly affected soil moisture dynamics, since the different cropping systems induced 16.3% of functional heterogeneity. Soil tillage was not found to have any significant effect.

Second, the suggested approach opens up new possibilities to analyze and better understand complex soil system behavior. Functional averaging, that is, averaging the signal transformation behavior of the vadose zone, provides time series representing the most relevant characteristics of the system behavior. The approach does not require a priori assumptions about the nature of physical processes, since it is solely based on information provided by the data. Thus, it provides model-independent information on how individual effects contribute to the observed dynamics, making it especially valuable for model building, calibration, and evaluation.

This study was funded by the German Research Foundation DFG (project Li 802/3-1). We would like to thank Wilfried Hierold for providing soil description data, Ralf Dannowski for preparing information about climatic conditions, and Dietmar Barkusky for providing land management data. We are also grateful to Anne-Kathrin Schneider, Christian Lehr, Steven Böttcher, Florian Reverey, and Björn Thomas for proofreading and fruitful discussions, as well as two anonymous reviewers for their valuable comments and suggestions.

This is an open access article distributed under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).