Estimating spatially distributed soil water content at small watershed scales based on decomposition of temporal anomaly and time stability analysis

Soil water content (SWC) is crucial to rainfallrunoff response at the watershed scale. A model was used to decompose the spatiotemporal SWC into a time-stable pattern (i.e., temporal mean), a space-invariant temporal anomaly, and a space-variant temporal anomaly. The spacevariant temporal anomaly was further decomposed using the empirical orthogonal function (EOF) for estimating spatially distributed SWC. This model was compared to a previous model that decomposes the spatiotemporal SWC into a spatial mean and a spatial anomaly, with the latter being further decomposed using the EOF. These two models are termed the temporal anomaly (TA) model and spatial anomaly (SA) model, respectively. We aimed to test the hypothesis that underlying (i.e., time-invariant) spatial patterns exist in the space-variant temporal anomaly at the small watershed scale, and to examine the advantages of the TA model over the SA model in terms of the estimation of spatially distributed SWC. For this purpose, a data set of near surface (0–0.2 m) and root zone (0–1.0 m) SWC, at a small watershed scale in the Canadian Prairies, was analyzed. Results showed that underlying spatial patterns exist in the space-variant temporal anomaly because of the permanent controls of static factors such as depth to the CaCO3 layer and organic carbon content. Combined with time stability analysis, the TA model improved the estimation of spatially distributed SWC over the SA model, especially for dry conditions. Further application of these two models demonstrated that the TA model outperformed the SA model at a hillslope in the Chinese Loess Plateau, but the performance of these two models in the GENCAI network (∼ 250 km) in Italy was equivalent. The TA model can be used to construct a high-resolution distribution of SWC at small watershed scales from coarseresolution remotely sensed SWC products.


Introduction
Soil water content (SWC) of surface soils exerts a major influence on a series of hydrological processes such as runoff and infiltration (Famiglietti et al., 1998;Vereecken et al., 2007;She et al., 2013a).Soil water content in the root zone is, in many cases, linked to vegetative growth (Wang et al., 2012;Ward et al., 2012;Jia and Shao, 2013).Obtaining accurate information on the spatiotemporal SWC is crucial for improving hydrological prediction and soil water management (Venkatesh et al., 2011;Champagne et al., 2012;She et al., 2013b;Zhao et al., 2010).While remote sensing has advanced SWC measurements of surface soils (<5 cm in depth) at basin (2500-25 000 km 2 ) and continental scales (Robinson et al., 2008), characterization of spatially distributed SWC at small watershed (0.1-80 km 2 ) scales still poses a challenge.A method is needed for estimating spatially distributed SWC in the near surface and root zone at watershed scales.Time stability of SWC, which refers to similar spatial patterns of SWC across different measurement times (Vachaud et al., 1985;Brocca et al., 2009), has been used for estimating spatially distributed SWC (Starr, 2005;Perry and Niemann, 2007;Blöschl et al., 2009).This method is conceptually ap-Published by Copernicus Publications on behalf of the European Geosciences Union.pealing, but assumes completely time-stable spatial patterns of SWC.
The time-stable pattern does not explain all of the spatial variances in SWC, indicating the existence of time-variant components (Starr, 2005).In order to identify underlying patterns of SWC that have time-variant components, the spatiotemporal SWC was decomposed into a spatial mean and a spatial anomaly.The spatial anomaly of the SWC was further decomposed into the sum of the product of time-invariant spatial patterns (EOFs) and temporally varying, but spatially constant coefficients (ECs) using the empirical orthogonal function (EOF) (Fig. 1) (Jawson and Niemann, 2007;Perry andNiemann, 2007, 2008;Joshi and Mohanty, 2010;Korres et al., 2010;Busch et al., 2012).Spatially distributed SWC estimates based on the decomposition of spatial anomaly outperformed those based on time-stable patterns (Perry and Niemann, 2007).
Recently, the spatiotemporal SWC was also decomposed into a temporal mean and a temporal anomaly (Mittelbach and Seneviratne, 2012) (Fig. 1).Previous studies indicated that the contribution of the temporal anomaly to the total spatial variance was notable (Mittelbach and Seneviratne, 2012;Brocca et al., 2014;Rötzer et al., 2015).These studies, however, only focused on surface soils at large scales (>250 km 2 ).Vanderlinden et al. (2012) suggested that the temporal mean may be further decomposed into its spatial mean and residuals, and the temporal anomaly may be further decomposed into space-invariant term (i.e., spatial mean of temporal anomaly) and space-variant term (i.e., spatial residuals of temporal anomaly) (Fig. 1).Note that the spatial variance in the temporal anomaly (Mittelbach and Seneviratne, 2012) equals that of the space-variant term of the temporal anomaly (Vanderlinden et al., 2012).The further decomposition of the temporal anomaly may be physically meaningful, because the space-invariant and space-variant terms in the temporal anomaly may be forced differently.However, the models of Mittelbach and Seneviratne (2012) and Vanderlinden et al. (2012) have not been used for estimating spatially distributed SWC.If the space-variant terms are ignored during the estimation of spatially distributed SWC, their models are equivalent to that based on time-stable patterns.Therefore, estimation of spatially distributed SWC may be improved by incorporating the space-variant term of the temporal anomaly if underlying (i.e., time-invariant) spatial patterns exist in the temporal anomaly.
To our knowledge, the importance of the space-variant term of the temporal anomaly and its physical meaning at small watershed scales is not well-known.Based on previous studies (Perry and Niemann, 2007;Mittelbach and Seneviratne, 2012;Vanderlinden et al., 2012), we assume soil water dynamics at watershed scales can be decomposed into three components (Fig. 1): (1) time-stable pattern (i.e., temporal mean, spatial forcing): the static factors such as soil and topography control the pattern; (2) space-invariant temporal anomaly (temporal forcing): the dynamic factors such as meteorological variables and vegetation change with time, and therefore modify SWC in time, regardless of spatial locations; and (3) space-variant temporal anomaly (interactions between spatial forcing and temporal forcing): this term represents interactions between static and dynamic factors.For example, SWC recharge introduced by a rainfall may be modified by topography through runoff processes; SWC loss triggered by evapotranspiration may be regulated by topography through solar radiation exposure.
The static factors may be persistent in the space-variant temporal anomaly, and their impacts on the space-variant temporal anomaly likely change with time.Thus, we hypothesize that some underlying (i.e., time-invariant) spatial patterns exist in the space-variant temporal anomaly, and their impacts can be modulated by a time coefficient, both of which can be obtained by the EOF method (Fig. 1).If the hypothesis is true, the estimation of spatially distributed SWC utilizing the EOF decomposition may outperform the one suggested by Perry and Niemann (2007).This is because: (1) the spatial anomaly, which was decomposed using the EOF in Perry and Niemann (2007), lumped the time-stable pattern and space-variant temporal anomaly together (Fig. 1); (2) the underlying spatial patterns in the spatial anomaly may not fully capture both time-stable patterns and patterns in the space-variant temporal anomaly due to the possible nonlinear relations between these two terms.
Therefore, the objectives were (1) to test the hypothesis that underlying spatial patterns exist in the space-variant temporal anomaly at small watershed scales and (2) to examine whether the decomposition of the space-variant temporal anomaly using the EOF has any advantages over the decomposition of the spatial anomaly (Perry and Niemann, 2007) for estimating spatially distributed SWC.Two steps were included in the estimation of spatially distributed SWC.First, the spatial mean SWC was upscaled from the SWC measurement at the most time-stable location using time stability analysis.Following this, the spatially distributed SWC was downscaled from the estimated spatial mean SWC.For the purpose of this study, spatiotemporal SWC data sets at depths of near surface (0-0.2 m) and root zone (0-1.0 m) from a Canadian Prairie landscape were used.Spatiotemporal SWC of samples taken 0-0.06 m from a hillslope (100 m) in the Chinese Loess Plateau and 0-0.15 m from the GENCAI network (∼ 250 km 2 ) in Italy were also used to further demonstrate conditions under which the decomposition of the spatial anomaly was beneficial to the estimation of spatially distributed SWC.

Study area and data collection
This study was mainly conducted in the Canadian Prairie pothole region (hereafter abbreviated as Canadian site) at St. Denis National Wildlife Area (52 • 12 N, 106 • 50 W) with an area of 3.6 km 2 .This area has a humid continental climate (Peel et al., 2007), and had a mean annual air temperature of 1.9 • C and a mean annual precipitation of 402 mm during the study period (Fig. 2).A variety of depressions, knolls, and knobs result in a sequence of undulating slopes (Biswas and Si, 2011).The elevation varies from 554.8 to 557.5 m.The soils are dominated by clay loam textured Mollisols (Soil Survey Staff, 2010) and covered by mixed grass, i.e., smooth brome grass (Bromus inermis) and alfalfa (Medicago sativa L.).The near-surface soil porosity ranges from 38 % (knolls) to 70 % (depressions).Calcium carbonates (CaCO 3 ) derived mostly from fragments of limestone rocks are common in the Canadian Prairies.The CaCO 3 is dissolved by the slightly acidic rainwater moving through the upper horizons and deposited to lower horizons.The heterogeneous amount of infiltrated water resulted in a varying depth of CaCO 3 layer ranging from almost 0 m in the knolls to 2.1 m in the depressions.A 576 m long sampling transect with 128 sampling locations spaced at 4.5 m intervals was established over several rounded knolls and depressions.At each location, a time domain reflectometry probe was used to measure SWC of the near-surface soil (0-0.2 m), and a neutron probe was used to collect SWC measurements at 0.2 m intervals between a depth of 0.2 and 1.0 m.The SWC was measured on a volumetric basis and expressed as a percentage (%) volume of water per unit soil volume.The SWC of the root zone was calculated by averaging the SWC of 0-0.2, 0.2-0.4,0.4-0.6,0.6-0.8, and 0.8-1.0 m.Soil water content was measured on 23 dates from 17 July 2007 to 29 September 2011.The SWC data set was collected in all seasons except winter, and accurately portrays the variations in soil water conditions in the study area.In addition to the SWC data set, the soil, vegetative, and topographical properties were obtained at each sampling location.These properties included soil particle components (clay, silt, and sand contents), bulk density, soil organic carbon (SOC) content for the surface layer, A horizon depth, C horizon depth, depth to the CaCO 3 layer, leaf area index, elevation, cos (aspect), slope, curvature, gradient, upslope length, solar radiation, specific contributing area, convergence index, wetness index, and flow connectivity.Detailed information on the measurements can be found in Biswas et al. (2012).The data sets from the Canadian site were used to demonstrate the following two aspects in detail: (1) different components of spatiotemporal SWC and their contributing factors, and (2) the advantages of the new decomposition method over the method suggested by Perry and Niemann (2007) in terms of the estimation of spatially distributed SWC.
To further test the applicability of the new method, we compared its performance at two other sites, covering both the hillslope and the large watershed scale.Along a hillslope of 100 m in length in the Chinese Loess Plateau, SWC of 0-0.06 m was measured 136 times from 25 June 2007 to 30 August 2008 by a Delta-T Devices Theta probe (ML2x) at 51 locations (Hu et al., 2011).The hillslope was covered by Stipa bungeana Trin.and Medicago sativa L. in sandy loam and silt loam soils.In the GENCAI network (∼ 250 km 2 ) in Italy, SWC of 0-0.15 m was measured by a TDR probe at 46 locations, 34 times from February to December in 2009 (Brocca et al., 2012(Brocca et al., , 2013)).The GENCAI area was dominated by grassland with a flat topography, in silty clay soils.

Statistical models for decomposing soil water content
Spatiotemporal SWC at small watershed scales was decomposed into three components: time-stable pattern, spaceinvariant temporal anomaly, and space-variant temporal anomaly.This model was compared to the one that decomposed SWC into spatial mean and spatial anomaly (Perry and Niemann, 2007).Both the space-variant temporal anomaly and spatial anomaly were decomposed using the EOF method.The two models are termed the temporal anomaly (TA) model and the spatial anomaly (SA) model.Figure 1 displays the differences between the two models.Each component will be explained in detail later.The explanation of nomenclatures is listed in Table A1.Because we focus on estimating spatial distribution of SWC at any given time, only spatial variances of SWC were taken into account.Therefore, the variance or covariance denotes the quantity in space without specifications.

2.2.1
The SA model Perry and Niemann (2007) expressed SWC at location n and time t (S tn ) as (Fig. 1): where S t n is the spatial mean SWC at time t (temporal forcing) and Z tn is the spatial anomaly of SWC (lumped spatial forcing and interactions).The subscript n (t) indicates a space (time) averaged quantity.
According to Perry and Niemann (2007), S t n can be estimated by remote sensing, water balance models, and in situ soil water measurement at a representative (or time-stable) location.The in situ soil water measurement method was selected because the representative location can be easily de-termined with prior SWC data sets.By measuring SWC only at the most time-stable location (s) and future time t (S ts ), S t n can be estimated using (Grayson and Western, 1998) where the s was identified using the time stability index of mean absolute bias error (Hu et al., 2010(Hu et al., , 2012)).The δ ts is the temporal mean relative difference of SWC at the s, which was calculated with prior measurements.Spatial anomaly (Z tn ) can be reconstructed by the sum of the product of time-invariant spatial structures (EOFs) and temporally varying coefficients (ECs) using the EOF method (Perry and Niemann, 2007;Joshi and Mohanty, 2010;Vanderlinden et al., 2012).The ECs correspond to the eigenvectors of the matrix of spatial covariance of the Z tn , and the EOFs are obtained by projecting the Z tn onto the matrix ECs as EOFs = Z tn ECs.The number of EOF (or EC) series equals the number of sampling dates.Each EOF series corresponds to one value at each location, and each EC series has one value at each measurement time.Each EOF is chosen to be orthogonal to other EOFs, and the lower-order EOFs account for as much variance as possible.The sum of variances of all EOFs equals the sum of variances of Z tn from all measurement times.
Usually, a substantial amount of variance can be explained by a small number of EOFs.Johnson and Wichern (2002) suggested the eigenvalue confidence limits method for selecting the number of EOFs.Once the number of significant EOFs at a confidence level of 95 % is selected, Z tn can be estimated as the sum of the product of significant EOFs and associated ECs as where EOF sig represents the significant EOFs of the Z tn obtained during model development, EC sig is the associated temporally varying coefficient, and the superscript T represents matrix transpose.Following Perry and Niemann ( 2007), the associated significant EC at time t (EC t ), is estimated by the cosine relationship between EC and S t n developed using prior measurements: where a, b, c, and d are the fitted parameters using prior measurements and S t n is estimated from Eq. (2).By using the continuous function, EC t can be estimated at any S t n values, which allows for the estimation of spatially distributed SWC at any soil water conditions.
Hydrol  Mittelbach and Seneviratne (2012) decomposed the S tn into a time-stable pattern (i.e., temporal mean) and a temporal anomaly component (Fig. 1): where M tn is the time-stable pattern (spatial forcing) controlled by static factors such as soil properties and topography; A tn refers to the temporal anomaly (lumped temporal forcing and interactions).The variance of SWC (σ 2 n (S tn )) is the sum of variance of the M tn (σ 2 n (M tn )), variance of the A tn (σ 2 n (A tn )), and two times of covariance between M tn and A tn (2 cov (M tn , A tn )), which can be expressed as: Because the A tn in Mittelbach and Seneviratne ( 2012) is a lumped term, it can be further decomposed into spaceinvariant temporal anomaly (A t n, i.e., temporal forcing) and space-variant temporal anomaly (R tn , i.e., interactions) ( Vanderlinden et al., 2012).At a watershed scale, the A t n is controlled by temporally varying factors such as meteorological variables and vegetation.Positive and negative A t n correspond to relatively wet and dry periods, respectively.The R tn refers to the redistribution of A t n among different locations due to the interactions between spatial forcing and temporal forcing.For example, soil and topography regulate how much rainfall enters soil and how much water runs off or runs on at a location.This, in turn, dictates vegetation growth in a water-limited environment.Therefore, S tn can also be expressed as (Fig. 1) The temporal trends of A t n in Eq. ( 7) and S t n in Eq. ( 1) are the same as both represent temporal forcing.Because the A t n is space-invariant and orthogonal to the M tn and R tn in a space, σ 2 n (S tn ) in Eq. ( 6) can also be written as where cov (M tn , R tn ) is the covariance between the M tn and R tn , and and  σ 2 n (R tn ) out of the σ 2 n (S tn ) are calculated.The cov (M tn , R tn ) can be negative at some conditions, for example, when the depressions correspond to greater M tn and more negative R tn values in the discharge periods.This resulted in percentage of σ 2 n (M tn ) and σ 2 n (R tn ) > 100 % and percentage of 2 cov (M tn , R tn ) < 0 % (Mittelbach and Seneviratne, 2012;Brocca et al., 2014;Rötzer et al., 2015).If R tn is zero at any time or location, there are no interactions between spatial forcing and temporal forcing, σ 2 n (S tn ) and the spatial trends of SWC are consistent over time.Therefore, R tn is directly responsible for temporal change in the spatial variability of SWC.
If some underlying spatial patterns exist in R tn , R tn can be reconstructed by the sum of the product of timeinvariant spatial structures (EOFs) and time-dependent coefficients (ECs) using the EOF method.Note that the number of EOF (or EC) series also equals the number of sampling dates.
For estimation of spatially distributed SWC, R tn is estimated by the same method as Z tn using Eq. ( 3).The M tn is estimated with prior measurements by where m is the number of previous measurement times, and A t n is estimated by: where M t n is the spatial mean of M tn , and S t n is estimated from SWC measurements at the most time-stable location using Eq. ( 2).The Pearson correlation coefficient (R) is used to explore the linear relationships between various spatial components in the two models (i.e., EOF1 of the Z tn in the SA model, M tn , and EOF1 of the R tn in the TA model) and environmental factors (i.e., soil, vegetative, and topographical properties).The multiple stepwise regressions are conducted to determine the percentage of variations in the spatial components which the controlling factors explain.

Validation and performance parameter
The TA model is more complicated than the SA model.In order to evaluate the two models for parsimony, AICc values are calculated (Burnham and Anderson, 2002) as where k is the number of parameters, n is the sample size, and RSS is the residual sum of squares.Both cross-validation and split sample validation are used to estimate SWC distribution with both models.For the cross-validation, an iterative removal of 1 of the 23 dates is made for model development, and the SWC along the transect corresponding to the removed date is estimated iteratively.For the split sample validation, SWC from 14 dates of the first 2 years (from 17 July 2007 to 27 May 2009) is used for model development, and the SWC distribution of 9 dates in the second 2 years (from 21 July 2009 to 29 September 2011) is estimated.
The Nash-Sutcliffe coefficient of efficiency (NSCE) is used to evaluate the quality of estimation of spatially distributed SWC, which is expressed as where σ 2 measure is the variance of measured SWC, and σ 2 ε is the mean squared estimation error.A larger NSCE value implies a better quality of estimation.A paired samples T test is used to test whether the NSCE values between the TA model and the SA model are statistically significant at P < 0.05.
Many factors may affect the relative performance of spatially distributed SWC estimation between the TA model and the SA model.First, the degree of outperformance of the TA model over the SA model may depend on the amount of R tn variance considered in the TA model.On one hand, the two models are identical if variance of R tn is close to zero or there are negligible interactions between the spatial and temporal components (Fig. 1).On the other hand, if no underlying spatial patterns exist in the R tn or the underlying spatial patterns accounted for little variance of the R tn , the outperformance will also be very limited.Therefore, the greater the variance of R tn considered in the TA model, the more likely the TA model can outperform the SA model.Second, the way of EOF decomposition may also affect the relative performance.In the SA model, EOF decomposition is performed on lumped time-stable patterns (M tn ) and space-variant temporal anomaly (R tn ).In the TA model, however, EOF decom-position is made only on the R tn .In theory, the two models will be identical if the M tn and the first underlying spatial pattern (i.e., EOF1) of the R tn were perfectly correlated.If a nonlinear relationship exists between them, lumping the M tn and R tn together, as in the SA model, would weaken the model performance as compared to the TA model.From this aspect, the greater deviation from a linear relationship between the M tn and EOF1 of the R tn , may lead to a greater outperformance of the TA model over the SA model.Finally, the performances of both models rely on the estimation accuracy of the EC t which depends on both goodness of fit of the cosine function (i.e., Eq. 4) and estimation accuracy of the S t n.Because the same S t n values are used for the two models, the relative performance of the two models is related to the goodness of fit of Eq. (4).

Spatial mean (S t n) and spatial anomaly (Z tn )
The values of spatial mean (S t n) in the SA model varied with the seasons (Fig. 3a).tively great S t n values.In the summer, however, even 1 month after large rainfall events (such as on 19 July 2008 and 21 June 2009), the high evapotranspiration by fast-growing vegetation resulted in small S t n values.The values of S t n also varied between inter-annual meteorological conditions.In 2008, there was less precipitation and higher air temperature than in 2010 (Fig. 2).As a result, S t n was relatively smaller in 2008 than in 2010.
The spatial patterns of spatial anomaly (Z tn ) were similar to those of the original SWC patterns (Fig. 3a).The values of Z tn in wet periods (e.g., 13 May 2011) were much greater than in dry periods (e.g., 23 August 2008) in depressions (e.g., at a distance of 123 and 250 m); at other locations, however, the spatial anomaly was slightly less in wet periods than in dry periods for both soil layers.Moreover, the spatial anomaly in depressions during the wet periods was much greater in the near surface than in the root zone.
When SWCs of all 23 dates were used for model development, only EOF1 was statistically significant (Fig. 4a), which accounted for 84.3 % (0-0.2 m) and 86.5 % (0-1.0 m) of the variances in the Z tn .Correlation analysis indicated that the spatial pattern of EOF1 in the Z tn was identical to the timestable patterns (M tn ) in the TA model (R = 1.0).The controls of EOF1 was therefore the same as those of M tn , and will be discussed later.The relationship between associated EC1 and S t n can be fitted well by the cosine function (R 2 = 0.73 at both the near surface and root zone) (Fig. 4b).For both soil layers, SOC, depth to the CaCO 3 layer, sand content, and wetness index are the dominant factors of M tn ; they together explained 74.5 % (near surface) and 75.6 % (root zone) of the variances in the M tn (Table 1).In addition, the temporal trend of A t n was the same as that of S t n in the SA model (Fig. 3a) as both represent temporal forcing.The R tn varied among landscape positions (Fig. 3b).At a sampling distance of 123 m (in a depression), R tn was negative in dry periods such as 23 August 2008 and positive in wet periods such as 13 May 2011.This was true for all depressions for both the near surface and the root zone.Therefore, topographically lower positions usually corresponded to more positive R tn during the wet periods and more neg- ative R tn during the dry periods.Furthermore, the absolute values of R tn were generally greater in the near surface than the root zone, indicating a greater space-variant temporal anomaly for shallower depths.

Time-stable pattern (M
The SWC variances and associated components (Eq.8) also varied with time (Fig. 5).Often, wetter conditions corresponded to greater σ 2 n (S tn ), as further indicated by moderate correlation between σ 2 n (S tn ) and S t n (R 2 of 0.51 and 0.38 for the near surface and the root zone, respectively).This was in agreement with others (Gómez-Plaza et al., 2001;Martínez-Fernández and Ceballos, 2003;Hu et al., 2011).Furthermore, there were greater σ 2 n (S tn ) values at the near surface than in the root zone, indicating greater variability of SWC in the near surface.
The time-invariant σ 2 n (M tn ) accounted for the σ 2 n (S tn ) with percentages ranging from 25 to 795 % for the near surface and from 40 to 174 % for the root zone (Fig. 5).The σ 2 n (M tn ) exceeded the σ 2 n (S tn ) mainly under dry conditions, such as July-October in 2008 and 2009.This excess was offset by the σ 2 n (R tn ) and 2 cov (M tn , R tn ), with the latter accounting for the σ 2 n (S tn ) negatively with mean absolute percentages of 210 % for the near surface and 17 % for the root zone.In the dry period, the absolute percentage of 2 cov (M tn , R tn ) was up to 1327 % for the near surface and 122 % for the root zone.These values are comparable to those in Mittelbach and Seneviratne (2012) and Brocca et al. (2014).
The σ 2 n (R tn ) accounted for less percentage of the σ 2 n (S tn ) than other components did (Fig. 5).The percentages of σ 2 n (R tn ) ranged from 11 to 632 % (arithmetic average of 118 %) for the near surface and from 6 to 48 % (arithmetic average of 19 %) for the root zone; the percentage of σ 2 n (R tn ) tended to be greater in drier periods.This indicates that the space-variant temporal anomaly cannot be ignored, particularly in dry conditions.Furthermore, the percentage of σ 2 n (R tn ) was greater in the near surface than in the root zone, confirming stronger temporal dynamics of soil water at the near surface.Compared with larger-scale studies (Mittelbach and Seneviratne, 2012;Brocca et al., 2014), the percentage of σ 2 n (R tn ) out of the σ 2 n (S tn ) at the near surface was greater, with a mean percentage of 118 %, versus 9-68 % in the other, larger-scale studies.This indicates that interactions between spatial and temporal forcing were stronger, resulting in relatively more intensive temporal dynamics of soil water in our study area than at larger scales.
Three significant EOFs of R tn for both soil layers were identified when SWC of all 23 dates were used for model development.The first three EOFs explained 61.1, 13.4, and 8.1 %, respectively, of the total R tn variance for the near surface, and 44.3, 20.2, and 12.4 %, respectively, of the total R tn variance in the root zone.Therefore, our hypothesis that underlying spatial patterns exist in the R tn was supported.Due to the negligible contribution of EOF2 and EOF3 to the estimation of spatially distributed SWC, only EOF1 is shown in Fig. 6a.The associated EC1 changed with soil water conditions (S t n) (Fig. 6b).When SWC was close to average levels, the EC1 was close to 0, resulting in negligible R tn .This was in accordance with Mittelbach and Seneviratne (2012) and Brocca et al. (2014), who showed that the spatial variance of the temporal anomaly was the smallest when water contents were close to average levels.The cosine function (Eq.4) explained a large amount of the variances in EC1 for both soil layers (R 2 = 0.76 at the near surface and 0.88 in the root zone).
The contribution of EOF1 to the space-variant temporal anomaly can be examined through the product of the EOF1 and the associated EC1.The EC1 values tended to be positive during wet periods and negative during dry periods (Fig. 6b); more positive EOF1 values were usually observed at locations with greater M tn values (Figs.3b and 6a).Therefore, the product of EOF1 and EC1 led to greater temporal SWC dynamics at wetter locations of both layers in both the wet and dry periods.
Depth to the CaCO 3 layer and SOC had significant, positive correlations with EOF1 for both soil layers (R ranging from 0.76 to 0.88; Table 1).They jointly accounted for 81.6 % (near surface) and 81.0 % (root zone) of the variances in EOF1.This implies that locations with a greater depth to the CaCO 3 layer and SOC, which correspond to wetter locations such as depressions, usually have greater temporal SWC dynamics during both wet and dry periods.

Estimation of spatially distributed SWC
When all 23 data sets were used and only EOF1 was considered, the TA model had an AICc value of 4093 for the near surface and 562 for the root zone, while the corresponding values for the SA model were 6370 and 3460.This indicated that even when penalty for complexity was given, the TA model was better than the SA model.The two models in terms of spatially distributed SWC estimation are compared below.

The TA model
The R tn terms and associated EOFs differed slightly with each validation.The number of significant EOFs varied between one (accounting for 60 % of the total cases) and three for both soil layers.A paired samples T test indicated that more EOFs did not result in a significant increase of NSCE in the estimation of spatially distributed SWC for both validation methods.This is also supported by the increasing AICc values with the increasing number of parameters resulting from more EOFs (data not shown).This indicates that higherorder EOFs, even if they are statistically significant, are negligible for SWC prediction.Therefore, SWC distribution was estimated with EOF1 only.Estimated SWCs generally approximated those measured at different soil water conditions during the cross-validation (Fig. 7).However, on 27 October 2009, there were unsatisfactory overestimates at the 100-140 and 220-225 m locations near the surface (Fig. 7a).Unsatisfactory NSCE values  The poor performance obtained with the TA model on those dates (Fig. 8a) was a result of overestimation in depressions, which is shown for example on 27 October 2009 (Fig. 7a).These dates also corresponded to a high percentage of σ 2 n (R tn ) to the σ 2 n (S tn ) (203-439 %).For 23 August and 17 September in 2008, which were in dry periods, the percentage of σ 2 n (R tn ) at the near surface was also high (580 and 630 %).Because a fair amount of σ 2 n (R tn ) was accounted for with the TA model, the TA model performed satisfactorily (NSCE of 0.43 and 0.60).For the remaining 20 dates, the resulting NSCE value ranged from 0.38 to 0.90 in the near surface and from 0.65 to 0.96 in the root zone (Fig. 8).This suggests that the TA model was generally satisfactory, with better performance in the root zone than in the near surface.
During the split sample validation, the TA model resulted in SWC estimations with NSCE values ranging from 0.61 to 0.85 near the surface and from 0.32 to 0.92 in the root zone, with exception of 2 days (27 August 2009 and 27 October 2009 with NSCE values of −2.63 and −5.12, respectively) at 0-0.2 m (Fig. 8).This suggested that the TA model performed well in estimating spatially distributed SWC patterns except on 27 August 2009 and 27 October 2009 at 0-0.2 m.The estimation in the root zone was also generally better than in the near surface.

Comparison with the SA model
One significant EOF of Z tn was identified for both soil layers, irrespective of the validation method.The SA model with only EOF1 produced reasonable SWC estimations for both validations in all dates in the root zone and in every date except five dates (23 August 2008, 17 September 2008, 22 October 2008, 27 August 2009, and 27 October 2009) in the near surface (Fig. 8).Similarly, when more EOFs were included, NSCE values did not increase significantly (data not shown) and consequently, estimation of spatially distributed SWC was not improved.This was because EOF2 and EOF3 together explained a very limited (< 10 %) amount of variability of Z tn and thus had low predictive power in terms of variance.
The difference in NSCE values between the TA and SA models for both validations are presented in Fig. 9. Generally, the difference decreased as A t n increased, and then slightly increased with a further increase in A t n.A paired samples T test indicated that the NSCE values of the TA model were significantly (P < 0.05) greater than those of the SA model for both soil layers, irrespective of validation methods.This indicates that the TA model outperformed the SA model, particularly in dry conditions.This was because when the soil was dry, there was a high percentage of σ 2 n (R tn ), and thus strong variability in the space-variant temporal anomaly.

A hillslope in the Chinese Loess Plateau
On average, the σ 2 n (M tn ), σ 2 n (R tn ), and 2 cov (M tn , R tn ) accounted for 53, 74, and −27 % out of the σ 2 n (S tn ), indicating that both time-stable pattern and temporal anomalies were the main contributors to the σ 2 n (S tn ).The EOF analysis showed that only the EOF1 was statistically significant for both the R tn and Z tn , and the EOF1 explained 23 and 47 % of the total variances of R tn and Z tn , respectively.This illustrated that underlying spatial patterns exist in the R tn on the hillslope.Cross-validation was used to estimate the spatially distributed SWC along the hillslope.The results showed that the NSCE varied from −4.25 to 0.83 (TA model) and from −4.30 to 0.81 (SA model), with a mean value of 0.25 and 0.19, respectively (Fig. 10a).A paired samples T test showed that the NSCE values for the TA model were significantly (P < 0.05) greater than those for the SA model, indicating that the TA model outperformed the SA model.As Fig. 10a shows, the outperformance was greater when SWC deviated from intermediate conditions, especially for dry conditions, which was similar to the Canadian site.

The GENCAI network in Italy
The σ 2 n (M tn ), σ 2 n (R tn ), and 2 cov (M tn , R tn ) accounted for 38, 68, and −7 % out of the σ 2 n (S tn ) (Brocca et al., 2014), indicating the dominant role of temporal anomalies in SWC variability.The first three EOFs of the R tn explained 19, 16, and 8 % of the total σ 2 n (R tn ), and no EOFs were statistically significant, indicating that no underlying spatial patterns exist in the R tn .The EOF1 of the Z tn was significant and accounted for 37 % of the variances in the Z tn .Although the EOF1 of the R tn was not significant, it was considered in the TA model for estimating spatially distributed SWC.The cross-validation indicates that the NSCE varied from −0.79 to 0.50 (TA model) and from −0.87 to 0.56 (SA model), with mean values of 0.09 and 0.08, respectively (Fig. 10b).The SWC estimation based on these two models was not satisfactory except for a few days.As Fig. 10b shows, the differences in NSCE values between the two models were scattered around 0. A paired samples T test showed that the NSCE values between the TA model and the SA model were not significant (P < 0.05), indicating no differences in estimating spatially distributed SWC between these two models.The R tn played an important role in the temporal change in spatial patterns of the SWC.The underlying spatial patterns and physical meaning in the R tn were examined in our study for the first time.Although three significant EOFs of the R tn existed in some cases, only EOF1 rather than higherorder EOFs of the R tn should be considered for the spatially distributed SWC estimation.Among many factors influencing the EOF1 of the R tn , depth to the CaCO 3 layer followed by the SOC, were the most important factors.Depressions have deeper CaCO 3 layers than knolls, and the shallow CaCO 3 layer on knolls limited water infiltration during rainfall or snowmelt, resulting in less water recharge on knolls than in depressions.The depth to CaCO 3 layer and SOC were negatively correlated with elevation (R = −0.54,P < 0.01).Therefore, the influence of depth to CaCO 3 layer and SOC partially reflected the role of topography in driving snowmelt runoff along slopes in the spring, which contributes to increasing water recharge in depressions.As already demonstrated, topographically lower positions corresponded to more negative R tn during the dry periods.This implies that depressions lost more water during discharge.This is because depressions usually corresponded to vegetation with a larger leaf area index, which would result in higher evapotranspiration and more water loss during discharge periods.
As Table 1 shows, both the depth to the CaCO 3 layer and SOC controlled the M tn .This was because deeper CaCO 3 layers and higher SOC were observed in depressions where soils were usually wetter in most of the year because of the snowmelt runoff in the spring and rainfall runoff in the summer and autumn (van der Kamp et al., 2003).Therefore, the roles of soil and topography were two-fold: On one hand, they were highly correlated with the time-stable patterns and thus the time stability of SWC (Gómez-Plaza et al., 2000;Mohanty and Skaggs, 2001;Grant et al., 2004); on the other hand, soil and topography, interplaying with temporal forcing, triggered local-specific soil water change and destroyed time stability of SWC.Their roles in protecting time stability persisted, but their roles in destroying time stability varied with time.Greater σ 2 n (R tn ) implies greater contribution of these factors in soil water dynamics, resulting in less time stability of SWC.

Model performance for spatially distributed SWC estimation
The outperformance of the TA model for estimating spatial SWC at the Canadian site and Chinese site can be partly explained by the high percentages (average of 19-118%) of the σ 2 n (R tn ) out of the total variance.When SWC is close to average levels, R tn is also close to zero, resulting in negligible percentage of σ 2 n (R tn ).In this case, the soil water patterns are stable in time, the SA model performs well, and there will be little difference between these two models.As is well known, the spatial patterns in soil water content are inherently time unstable.For example, when evapotranspiration becomes the dominant process at the small watershed scale, more water will be lost in depressions due to the denser veg-etation than on knolls (Millar, 1971;Biswas et al., 2012), effectively diminishing the spatial patterns and increasing temporal instability.In this case, the σ 2 n (R tn ) accounts for more percentage of the total variance (e.g., high up to 632 %) and the TA model may outperform the SA model.This explained why the outperformance of the TA model was more obvious in the dry conditions.For the GENCAI network in Italy, although the σ 2 n (R tn ) accounted for 68 % of the total variance, the performance of the TA model was identical to the SA model.This was because there were no underlying spatial patterns in the R tn .Similarly, because the first underlying spatial pattern (i.e., EOF1) explained greater percentages of the σ 2 n (R tn ) at the Canadian site (44-61 %) than the Chinese site (23 %), the outperformance of the TA model over the SA model was more obvious at the former site (Figs. 9 and 10a).Therefore, the TA model is advantageous only if the percentage of σ 2 n (R tn ) out of the total variance is substantial and underlying spatial patterns exist in the R tn .
The existence of underlying spatial patterns in the R tn is related to the controlling factors, which may be scale specific.At small scales, static factors such as the depth to the CaCO 3 layer and SOC at the Canadian site may affect not only the time-stable patterns but also the R tn .The persistent influence of static factors on the R tn resulted in significant underlying spatial patterns in the R tn .Thus, the TA model outperformed the SA model at the small scales.At large scales such as the basin scale or greater, time-stable patterns may be controlled by, in addition to soil and topography (Mittelbach and Seneviratne, 2012), the climate gradient (Sherratt and Wheater, 1984); at those scales, R tn is more likely to be controlled by the meteorological anomaly (i.e., spatially random variation) (Walsh and Mostek, 1980), and the effects of soil and topography may be reduced.Consequently, spatial patterns in the R tn may be weakened and the TA model may have no advantages over the SA model such as for the Italian site.
The M tn and the underlying spatial patterns (EOF1) in the R tn were controlled by the same spatial forcing (e.g., depth to CaCO 3 layer and SOC) at the Canadian site (Table 1), and they were correlated with an R 2 of 0.83 for the near surface and 0.42 for the root zone.Although the relationships between M tn and R tn were strong, they were not strictly linear, suggesting that M tn and R tn were affected differently by these factors.Therefore, the nonlinear relationship between M tn and R tn partially contributed to the outperformance of the TA model over the SA model.
The relationship between the S t n and EC1 was better fitted by the cosine function in the TA model than the SA model (Figs.4b and 6b), with R 2 of 0.76 versus 0.73 in the near surface and 0.88 versus 0.73 in the root zone.The reduced scatter in the S t n and EC1 relationship for the TA model may also partly explain the outperformance of the TA model over the SA model.
Therefore, the outperformance of the TA model over the SA model depends on counterbalance among the variance of R tn explained in the TA model, the linear correlation between the M tn and EOF1 of the R tn , and the goodness of fit for the S t n and EC1 relationship.For example, the variance of EOF1 in the R tn for the near surface (i.e., 264 % 2 ) was much greater than that for the root zone (i.e., 43 % 2 ).However, M tn and underlying spatial patterns (EOF1) in the R tn in the root zone deviated more from a linear relationship, and the reduced scatter in the S t n and EC1 relationship in the TA model was more obviously in the root zone than in the near surface.As a result, the outperformance of the TA model was comparable between the near surface and root zone at the Canadian site (Fig. 9).
In the real world, the relations between the M tn and underlying spatial patterns in the R tn may rarely be perfectly linear.Therefore, when underlying spatial patterns exist in the R tn and the R tn has substantial variances, the TA model is preferable to the SA model for the estimation of spatially distributed SWC.On the other hand, when underlying spatial patterns do not exist in the R tn or the R tn has negligible variances, the SA model may be selected although these two models yield the same quality of SWC estimation.This is because the TA model needs one more spatial parameter (i.e., M tn ) than the SA model.
Previous studies on SWC decomposition mainly focus on near-surface layers (Jawson and Niemann, 2007;Perry andNiemann, 2007, 2008;Joshi and Mohanty, 2010;Korres et al., 2010;Busch et al., 2012).This study decomposed spatiotemporal SWC using the TA model for both the near surface and the root zone.The results showed that the estimation of spatially distributed SWC at small watershed scales was improved by the TA method that considers the R tn .The σ 2 n (M tn ) was greater than the σ 2 n (R tn ) (Fig. 5), indicating that time stability was more important than time instability for SWC estimation.For the three dates in the fall (i.e., 22 October 2008, 27 August 2009, and 27 October 2009), strong evapotranspiration and deep drainage in depressions resulted in a much lower SWC at the near surface than in the spring.This resulted in reduced time stability of SWC patterns and poor performance of both models in terms of SWC evaluation (Fig. 8a).Because of the stronger time stability of SWC in deeper soil layers (Biswas and Si, 2011), SWC evaluation was more accurate for soil layers extending from the surface to greater depth.This is particularly important because SWC data for deeper soil layers in a watershed is more difficult to collect than that of surface soil.

Conclusions
The TA model was used to decompose spatiotemporal SWC into time-stable patterns M tn , space-invariant temporal anomaly A t n, and space-variant temporal anomaly R tn .This study indicated that underlying spatial patterns may exist in the R tn at small scales (e.g., small watersheds and hillslope) but may not exist at large scales such as the GENCAI network (∼ 250 km 2 ) in Italy.This was because the R tn at small scales was driven by static factors such as depth to the CaCO 3 layer and SOC at the Canadian site, while the R tn at large scales may be dominated by dynamic factors such as meteorological anomaly.Compared to the SA model, estimation of spatially distributed SWC was improved with the TA model at small watershed scales.This was because the TA model considered a fair amount of spatial variance in the R tn , which was ignored in the SA model.Furthermore, the improved performance was observed mainly when there was less or more soil water than the average level, especially in drier conditions due to the high σ 2 n (R tn ) value.
This study showed that outperformance of the TA model over the SA model is possible when σ 2 n (R tn ) accounts for substantial variance of SWC, and significant spatial patterns (or EOFs) exist in the R tn .Further application of the TA model for the estimation of spatially distributed SWC at different scales and hydrological backgrounds is recommended.If the TA model parameters (i.e., M tn , EOF1 of the R tn , and relationship between EC and S t n) are obtained from historical in situ SWC data sets, a detailed spatially distributed SWC of near-surface soil at watershed scales can be constructed from remotely sensed SWC.Note that both models rely on in situ SWC measurements for model parameters.Therefore, future research should be conducted to estimate spatially distributed SWC in un-gauged watersheds based on the estimation of the model parameters using pedotransfer functions.The codes for decomposing SWC with the SA and TA models and related EOF analysis were written in Matlab and are freely available from the authors upon request.

Figure 1 .
Figure 1.Decomposition of spatiotemporal soil water content (SWC) in different models.

Figure 2 .
Figure 2. Daily mean air temperature and precipitation during the study period.

Figure 3 .
Figure 3. Components of soil water content in (a) the SA model (spatial mean soil water content S t n and spatial anomaly Z tn ) and in (b) the TA model (time-stable pattern M tn , space-invariant temporal anomaly A t n, and space-variant temporal anomaly R tn ) for 0-0.2 and 0-1.0 m.Also shown is the elevation.
Figure3bdisplays the three components in the TA model.The first component M tn fluctuated along the transect, with high values in depressions and low values on knolls; the M tn also had greater spatial variability in the near surface (variance = 36.7 % 2 ) than in the root zone (variance = 19.5 % 2 ).For both soil layers, SOC, depth to the CaCO 3 layer, sand content, and wetness index are the dominant factors of M tn ; they together explained 74.5 % (near surface) and 75.6 % (root zone) of the variances in the M tn (Table1).In addition, the temporal trend of A t n was the same as that of S t n in the SA model (Fig.3a) as both represent temporal forcing.The R tn varied among landscape positions (Fig.3b).At a sampling distance of 123 m (in a depression), R tn was negative in dry periods such as 23 August 2008 and positive in wet periods such as 13 May 2011.This was true for all depressions for both the near surface and the root zone.Therefore, topographically lower positions usually corresponded to more positive R tn during the wet periods and more neg-

Figure 4 .
Figure 4. (a) The EOF1 of the spatial anomaly Z tn and (b) relationships of associated EC1 versus spatial mean soil water content Z tn fitted by the cosine function (Eq.4).

Figure 5 .
Figure 5. Spatial variances of different components in Eq. (8) expressed in % 2 (upper panel) and as percentage (lower panel) for (a) 0-0.2 and (b) 0-1.0 m.Spatial mean soil water content S t n on each measurement day is also shown.

Figure 6 .
Figure 6.(a) The EOF1 of the space-variant temporal anomaly R tn and (b) relationships of associated EC1 versus spatial mean soil water content S t n fitted by the cosine function (Eq.4).

Figure 7 .
Figure 7.Estimated soil water content (SWC) versus measured SWC for three dates at different soil water conditions (23 August 2008, 27 October 2009, and 13 May 2011 are associated with relatively dry, medium, and wet days, respectively) using the TA model for (a) 0-0.2 and (b) 0-1.0 m.

Figure 8 .
Figure 8.The Nash-Sutcliffe coefficient of efficiency (NSCE) of soil water content estimation using the TA and SA models for (a) 0-0.2 and (b) 0-1.0 m for both cross-validation (CV) and split sample validation (SV).At 0-0.2 m, three dates (22 October 2008, 27 August 2009, and 27 October 2009) as indicated by green lines present negative NSCE values (−4.05, −1.83, and −3.81, respectively, for the CV on the three dates; −2.63 and −5.12, respectively, for the SV on the latter two dates).Spatial mean soil water content S t n on each measurement day is also shown.

Figure 9 .
Figure 9. Nash-Sutcliffe coefficient of efficiency (NSCE) difference between the TA and SA models in terms of soil water content estimation using both cross-validation (CV) and split sample validation (SV) as a function of space-invariant temporal anomaly A t n for (a) 0-0.2 and (b) 0-1.0 m.

Figure 10 .
Figure 10.Nash-Sutcliffe coefficient of efficiency (NSCE) difference between the TA and SA models in terms of soil water content estimation using cross-validation as a function of space-invariant temporal anomaly A t n for (a) 0-0.06 m of the Chinese Loess Plateau hillslope and (b) 0-0.15 m of the GENCAI network in Italy.The NSCE values for both models are also shown.

Hu and B. Si: Estimating spatially distributed soil water content at small watershed scales 577Table 1 .
Pearson correlation coefficients between time-stable pattern M tn , EOF1 of space-variant temporal anomaly R tn and various properties.
1 Percent of variance explained by the controlling factors obtained by the multiple stepwise regressions.* Significant at P < 0.05; * * significant at P < 0.01.