Soil Surface Sealing Effect on Soil Moisture at a Semiarid Hillslope: Implications for Remote Sensing Estimation

2.1. The Study Site

The LTER (long-term ecological research) Lehavim site (31°20′N, 34°45′E) is located in the Negev Desert in the center of Israel (Figure 1). The climate is characterized as semiarid, with an annual rainfall of 290 mm, usually falling between November and March. Within the LTER, a typical hillslope was chosen (0.075 km²) composed of east- and west-facing slopes. The soils are brown lithosols and arid brown loess, which are prone to surface sealing, with a mean soil texture of 23% sand, 56% silt and 21% clay. The hillslope soil depth ranges from 0 to 55 cm, with a mean value of 22 cm. Rock outcroppings cover 11% of the entire hillslope and are concentrated mainly in the interfluvial area. The vegetation (24% cover) is characterized by scattered dwarf shrubs (the dominant species is Sarcopoterium spinosum) and patches of herbaceous plants [57].

2.2. Methodological Approach

A physically-based modeling approach relying on the numerical solution of the flow equations using Hydrus 1D [58] and measured climatic data was applied to simulate the mean soil moisture dynamics across a hillslope represented by 8240 cells (3 × 3 m). Two cases of bare soil profiles were considered: (1) uniform (unsealed case); (2) topped by a completely formed seal layer (sealed case). An extensive validation of this approach for the upper 15 cm of the soil profile was performed previously [59]. The characterization of the hillslope as an amalgamation of 1D cells was compared with a 2D approach for a subset of the data and was found suitable in reconstructing soil moisture measured in the field (see Appendix A).

To corroborate the model’s ability to reconstruct the soil moisture dynamics at shallower soil depths, the simulated soil moisture was compared with gravimetric soil-moisture data repeatedly collected from the 0–10 cm-deep upper layer at 7 to 13 different locations on the hillslope (Figure 1B), during nine field campaigns spanning the 2010–2011 rainy season (for a total of 82 samples). For each sampling date, each gravimetric sample, in each location, was averaged over three replicates. Details regarding the 2010–2011 rainfall season and sampling dates can be found in Sela et al. [59]. Gravimetric soil moisture data were converted to volumetric soil moisture using soil bulk density calculated following Saxton et al. [60], using soil texture at each sampling site as the input.

To address the objectives of the study, two datasets of soil moisture were created by averaging the simulated data for the hillslope in two different ways: the first dataset represents the soil layer, which corresponds to satellite observations. While in reality, the thickness of this layer is affected by soil moisture, this is neglected at this stage, and the layer thickness is assumed to be 0–2 cm (spatially uniform and constant in time). The second simulated soil moisture dataset represents the soil layer that is usually sampled by soil moisture probes used to validate satellite observations. In this study, the thickness of this layer was assumed to be 0–10 cm and to be spatially uniform. Using a soil moisture dataset created by a hydrological model allows more flexibility in studying near-surface soil moisture dynamics at various temporal scales, from a single rainfall event to the multi-seasonal scale, than if real RS data were used.

The use of model-based soil moisture realizations enables the study of SMTS dynamics at any prescribed spatial scale, both with and without accounting for a seal layer at the surface, thus allowing the quantitative evaluation of its specific effect. In this study, all analyses were applied to a 30 × 30-m spatial scale, corresponding to the scale of synthetic aperture radar (SAR) satellites, such as RadarSat-1 [61]. Consequently, the studied hillslope was mapped to a 30 × 30-m grid, resulting in 69 cells (Figure 1B).

The initial soil moisture conditions at the beginning of each simulation were assumed uniform, both spatially and with depth, and consisted of a soil moisture value of $\theta =0.1\left(\frac{{\text{cm}}^3}{{\text{cm}}^3}\right)$. For each of the simulated seasons, mean soil moisture values for the 0–2 cm (representing the penetration depth of the RS sensor) and 0–10 cm (representing the sampled depth of the ground sensors) layers were computed on an hourly basis in each 3 × 3-m model cell and extracted every 6 h for both the sealed and the unsealed conditions. The known discrepancy between the point scale of the ground validation measurements and the RS sensor footprint scale was not addressed here, and the comparison between the 0–2 cm and 0–10-cm soil depths was conducted directly at the RS sensor footprint scale of 30 × 30 m to emphasize the net effect of the seal layer. Therefore, simulated soil moisture data were upscaled from the 3 × 3-m model cell to the 30 × 30-m RS sensor cells by averaging all respective 3 × 3-m cells contained within each RS sensor cell, for each time step. Finally, to study the factors controlling the bias between mean soil moisture at the 0–2 cm and 0–10 cm depths, the mean difference between the simulated soil moisture at these depths during all seasons was computed for each RS sensor cell and was used as an input for a multi-parameter, stepwise regression (SPSS version 20, IBM, Inc. Armonk, NY, USA).

To study how surface sealing affects the sampling locations that best represent the hillslope mean soil-moisture, a temporal stability analysis [44] was applied using simulated soil moisture data for the 0–2-cm layer, representing the RS sensor output. A long-term dataset of measured climate records was used for this analysis, and both sealed and unsealed conditions were considered. Finally, to study the factors affecting the spatial distribution of temporally stable cells, an analysis of variance (ANOVA) was applied.

2.3. Modeling Approach

For a detailed description of the modeling approach, the reader is referred to Sela et al. [59]. In short, the Hydrus 1D software [58] is used to solve numerically the flow equations in the 8240 independent cells (3 × 3 m) that represent the hillslope study site. MATLAB (Release 2012b, The MathWorks, Inc., Natick, MA, USA) code was used to automate the simulations on an hourly basis during the rainfall season. The 8240 bare soil cells were characterized using an extensive database of topographic and soil-hydraulic parameters collected in the field. Assuming that potential evaporation (PE) is mainly controlled by solar radiation, the model accounts for the spatial variability of PE at the study site using a correction index, E_i [59]. This index adjusts the applied PE values in each cell according to aspect and hillslope location, using the ratio between solar radiation in each model cell and the respective radiation of a flat surface. Using the ratio between the rock and soil covers in each cell, rainfall intensity over all cells was corrected to account for the contribution from local rock outcrops. Boundary conditions consist of an atmospheric condition for the top (open to rainfall and evaporation fluxes) and free drainage for the bottom of the solution domain, reflecting the high level of uncertainty regarding the local bedrock permeability, which was found to be highly disintegrated [62]. In the unsealed case, a uniform soil profile was assumed. When a seal layer was accounted for, a two-layer soil profile was assumed, with the seal having a thickness of 2 cm based on observations of a similar soil type [63,64]. Under field conditions, the thickness of the seal layer is expected to depend on local factors, such as slope inclination and the percentage of fine material in the soil. However, no such relations had been established yet for the study site. Therefore, in this analysis, the seal layer thickness is assumed to be spatially uniform throughout the hillslope.

The hydraulic parameters of the seal layer were directly related to the underlying soil using the Mualem and Assouline [30] approach. The model simulates soil moisture dynamics at the bare soil between vegetation patches and, thus, neglects vegetation effect on infiltration and water redistribution. This approach could be justified at the study site, as the dominant shrubs (Sarcopoterium Spinosum) were found to have mainly vertical roots, extracting water from fissures at the bedrock [65]. Surface water ponding was not allowed in the Hydrus 1D model simulations. Rainfall excess was recorded, but no runoff routing was considered between the cells, based on very limited runoff connectivity, which was measured at the site [66]. Subsurface flow between the cells was not accounted for. Considering all the above-mentioned assumptions, the simulated results were verified extensively in the field. Several sampling campaigns of the 2–15-cm layer for soil gravimetric soil moisture measurements were carried out in the bare soil areas between vegetation patches during the 2010–2011 rainfall season. A comparison between simulated and measured mean soil moisture values has resulted in a satisfactory Nash and Sutcliffe [67] efficiency index of 0.7, indicating that the modeling approach and the assumptions were able to represent soil moisture dynamics in the top soil of the experimental hillslope.

2.4. Climate Data

A 20-y dataset (1993–2013) of daily rainfall measured at the Lahav meteorological station, located 4 km north of the study site, was used for the analysis. All precipitation from a particular rain event was assumed to fall during a 3-h period, the average duration of rainfall events at the study site. Potential evaporation for the years 1993–2009 was measured using a Class A pan at the Lahav station. For the period of 2010–2013, potential evaporation was calculated using the Penman–Monteith equation [68], based on data from the Beer Sheva meteorological station located 10 km south of the study site.

2.5. Time Stability Analysis

In this study, SMTS was analyzed using the mean relative difference (MRD) method [44]. The mean relative difference (δ̄₁) $\left(\frac{{\text{cm}}^3}{{\text{cm}}^3}\right)$ for each 30 × 30 m RS cell (i) was calculated using Equation (1):

The combination of these two indices leads to the root mean square error (RMSE) of the relative difference for each cell (Equation (3)) [43].

Lower RMSE values indicate locations with a higher time stability of soil moisture. These locations are considered to represent the best the mean soil moisture of the hillslope [43].

2.6. Statistical Analysis

To study the factors that affect the bias in validation, two statistical methods were applied: linear multiple parameter step-wise regression and analysis of variance (one way ANOVA). Five topographic and soil hydraulic variables (all of them serve also as model inputs) were chosen for the analysis: (1) soil depth; (2) rock cover; (3) the correction index, E_i for PE [59]; (4) the saturated water content value (θ_s) of the soil layer, which is equivalent to soil porosity; and (5) slope inclination. A moderate correlation (R = 0.53) was found between the relations of the soil depth-rock cover and the soil depth-slope inclination. A weak correlation was found for all other combinations of the regression parameters. For the regression analysis, we have used dimensionless standardized regression coefficients, calculated by scaling the sample mean value to zero and the sample variance to one [69]. Using standardized regression coefficients allows us to compare, with necessary caution, the relative importance of each explanatory parameter. One-way ANOVA can be used to test for significant differences between the means of multiple populations [70]. Here, it is applied to study if significant differences are found between populations grouped by topography or soil hydraulic parameters [43], highlighting variables that could be good predictors of SMTS. For this analysis, each of the five variables mentioned above was discretized into four groups, based on an equal interval classification.

3.1. Corroboration of the Model with Field Data

The comparison between simulated and measured mean soil moisture data for the 0–10 cm upper layer is depicted in Figure 2. Accounting for soil surface sealing improved model predictions for the lower soil moisture levels and increased the Nash and Sutcliffe [67] efficiency index from 0.59 for the unsealed case to 0.66 for the sealed case. This latter value is considered a satisfactory agreement according to the guidelines presented in Moriasi et al. [71]. Based on these results, we conclude that the modeling approach accurately simulates soil moisture dynamics in the upper soil layer and that the sealed case is better at representing soil moisture in the field than the unsealed case.

3.2. Differences between Simulated Soil Moisture for the RS Sensor Penetration Depth and Ground Validation Depth

To investigate the bias induced by using soil moisture measured at the 0–10-cm layer to validate soil moisture sensed at the 0–2-cm layer, the model was forced with the climatic data from three representative seasons: 2004–2005 (407 mm rainfall), 2006–2007 (275 mm rainfall) and 2007–2008 (203 mm rainfall), representing above average, average and less than average rainfall seasons, respectively. For each model cell, soil moisture data from all three seasons were combined to a continuous dataset and the difference between the mean soil moisture in the two layers was compared for both the sealed and unsealed systems (Equations (4) and (5), respectively):

3.3. Temporal Variation in the Error between the 0–2 and 0–10 Layers

The computed Diff values for Equations (4) and (5) for the 2006–2007 rainfall season are depicted in Figure 4. This season was chosen, because the seasonal cumulative rainfall depth during that year is similar to the long-term mean value of the rainfall data for the study site. A cyclic pattern for the difference between the mean soil moisture within the two layers was observed under the sealed and the unsealed conditions. Negative values were found to characterize the rainfall events (e.g., higher soil moisture at the 0–2-cm compared to the 0–10-cm layer) and positive values to characterize drying periods (e.g., higher soil moisture at the 0–10-cm compared to the 0–2-cm layer). The bias associated with the validation of remotely sensed soil moisture during rainfall events due to the different RS penetration and ground sampling depths in unsealed systems is a known issue that was addressed previously [21]. The seal layer was found to intensify this bias, especially during low-intensity rainfall events where, in the absence of runoff initiation, the percolation of water through the seal layer is impeded by the seal layer’s lower hydraulic conductivity. Interestingly, the difference in soil moisture values between the 0–2 and 0–10-cm layers during the drying periods was substantially smaller in the sealed system. This is attributed to the seal layer’s suppression of evaporation during the drying periods [33,59], which reduces the difference between the two layers over time. Conversely, in the unsealed system, the progression of a drying front through the soil profile led to an increase in the differences in soil moisture between the 0–2 and 0–10-cm soil layers over time. Figure 4C summarizes the net contribution of the seal layer to the difference in soil moisture between the 0–2 and 0–10-cm soil layers. The seal was found to add up to 10%–15% of the difference between the 0–2 and the 0–10-cm layers during rainfall events. However, immediately following a rainfall event, the seal layer starts to compensate for this increased difference, and by the time the drying period ends, the difference in moisture between the two layers can be reduced by 10%–15%.

The cyclic pattern in the differences in soil moisture between the 0–2 and 0–10-cm soil layers, driven by the rainfall events and drying periods (Figure 4), suggests that there might be a characteristic duration until the difference in soil moisture between these layers vanishes. To test this hypothesis, this duration was computed following each of the rainfall events of the three representative seasons (a total of 40 rainfall events). The results are summarized in Figure 5, which shows a histogram of the computed durations. In most of the cases, the difference in soil moisture between the layers vanished 10 h for the unsealed system and 30 h for the sealed system after rainfall, representing an ideal duration for sampling. In reality, validation ground sampling schemes is designed to coincide with the timing of satellite overpass and not according to the local climatic conditions. However, since the bias between the layers was found to follow a cyclic pattern, using these calculated durations can assist validation schemes by indicating whether a positive or negative bias in soil moisture validation could be expected during the actual satellite overpass.

3.4. Effects of a Seal Layer on SMTS and the Determination of Optimal Hillslope-Scale Sampling Locations

To study the effect a seal layer has on SMTS patterns at the hillslope scale, the model was forced with 20 years of measured climate data (1993–2013), and MRD and RMSE values (Equations (1)–(3)) were calculated for each 30-m cell based on the entire period. Figure 6 presents the MRD and RMSE values of the sealed and unsealed systems, ranked from the smallest to the largest, according to the MRD values. Sixteen percent and 14% of the cells in the sealed and unsealed systems, respectively, estimated the hillslope mean soil moisture within ±2%, thus exhibiting a relatively high SMTS. Fifty five percent and 39% of the cells captured up to ±5% of the mean hillslope soil moisture, while 87% and 78% of the cells captured up to ±10% of the mean hillslope soil moisture (for the sealed and unsealed systems, respectively). Overall, soil surface sealing improves SMTS, reduces the RMSE by 26%, the mean absolute MRD by 20% and its variance by 46%. In both the sealed and unsealed systems, the variance of the MRD values was greater in cells with soil moisture that is consistently higher than the mean value (positive values of the MRD) than in cells that were consistently drier then the hillslope mean. This result agrees with those of several previous studies [1,2,43,52,72]. The MRD values calculated in this study were lower than the values reported in other studies conducted on similar semiarid hillslopes or small catchments [47,72,73]. This could be related to the fact that we considered soil moisture data averaged over the 30-m footprint, while those studies used data originating from a network of point measurements.

Following the methodology presented in Jacobs et al. [43] and Joshi et al. [1], to further study whether topographic and soil hydraulic variables could be good predictors of the MRD values in a quantitative manner, the results of the ANOVA analysis are presented in Table 2. The results indicate that for both the sealed and unsealed systems, the difference in MRD between the groups is most significant for the saturated water content value and the slope inclination. Figure 7A depicts the mean MRD in each of the four slope inclination groups, together with the spatial distribution of each zone along the hillslope (Figure 7B). For the sealed system, both the higher and lower backslope (10–15 and 15–20 degrees) emerge as optimal locations for sampling, estimating the hillslope mean soil moisture within ±2%. For the unsealed system, only sampling from the lower backslope (15–20°) allowed the estimation of the hillslope mean soil moisture within ±2%. However, for both the sealed and unsealed systems, all four slope-inclination zones had highly variable MRD levels (e.g., relatively large standard deviations), indicating that slope inclination alone may not be a good predictor of optimal locations for sampling soil moisture. Figure 8 depicts the relation between MRD and the slope inclination and the saturated water content value. While the relation between the MRD and the slope inclination is scattered across the entire range of MRD values, a strong linear relation was found between the saturated water content and the MRD (R² values of 0.97 and 0.93 for the sealed and unsealed systems, respectively). This indicates that temporally stable hillslope locations (e.g., MRD of ±2%) could only be found within a narrow range of saturated water content (or porosity) values (θ_s = 0.482 − 0.486), closely representing the hillslope mean value (θ_s = 0.484).