On the observed hysteresis in field-scale soil moisture variability and its physical controls

The spatiotemporal variability of soil moisture (θ) has rarely been studied at the field scale across different seasons and sites. Here, we utilized 9 months of θ data in two semiarid ecosystems of North America to investigate the key relationship between the spatial mean (〈θ〉) and standard deviation (σθ) at the field-scale (∼100 m). Analyses revealed a strong seasonal control on the σθ versus 〈θ〉 relation and the existence of hysteretic cycles where wetting and dry-down phases have notably different behavior. Empirical orthogonal functions (EOFs) showed that θ variability depends on two dominant spatial patterns, with time-stable and seasonally varying contributions in time, respectively. Correlations between EOFs and land surface properties also indicated that θ patterns are linked to vegetation (terrain and soil) factors at the site with higher (lower) vegetation cover. These physical controls explained the observed hysteresis cycles, thus confirming interpretations from previous modeling studies for the first time.


Introduction
Soil moisture (θ) is a critical state variable controlling the feedback between the water, energy and carbon cycles (Entekhabi 1995). As a result, quantifying its spatiotemporal variability and identifying its associated physical controls at different scales is crucial for diagnostic and prognostic studies in the areas of hydrology, climate, atmospheric science, ecology, and agronomy, among others (Tao et al 2003, Entekhabi et al 2010, Seneviratne et al 2010, Mascaro and Vivoni 2012, Dillon et al 2016. A simple way to characterize the spatiotemporal variability of θ is through the relation between its spatial mean (〈θ〉) and variability (measured by the standard deviation, σ θ , or the coefficient of variation, CV=σ θ /〈θ〉). This relation has been investigated using observations and models over different spatial extents, ranging from field (∼100 m) to regional (>100 km) scales, and diverse settings (e.g., Famiglietti et al 2008, Brocca et al 2010, Li and Rodell 2013. Previous studies have shown that the σ θ versus 〈θ〉 relation has an upward convex shape independently of the spatial extent, while the range of 〈θ〉 depends on local climate. Moreover, Rosenbaum et al (2012) recently showed the presence of hysteretic cycles in the relation by analyzing one year of θ data collected in a humid basin of 0.27 km 2 in Germany. Prior to this study, the occurrence of hysteretic patterns was only identified and discussed through modeling experiments (Ivanov et al 2010, Vivoni et al 2010, Fatichi et al 2015, Ji et al 2015. While the general shape of the σ θ versus 〈θ〉 relation is explained by the bounded nature of θ, the physical factors (climate, vegetation, soil and topography) controlling the range of σ θ values observed for the same 〈θ〉, the possible existence of hysteresis, and the shape of the relation within different 〈θ〉 intervals are not fully understood (Vanderlinden et al 2012). A number of hypotheses have been formulated to explain the controls on these patterns based on synthetic numerical experiments, which highlighted the greater importance of vegetation (soil) properties in temperate and semiarid (humid) climates, where 〈θ〉 Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. tends to be low (high) (e.g., Lawrence andHornberger 2007, Fatichi et al 2015). However, the validation of these hypotheses has been prevented by the lack of high-resolution soil moisture observations, such as those provided by dense automated networks. This is especially true in semiarid areas, where distributed measurements would be required to capture the high spatial heterogeneity of soil and vegetation properties (Ochsner et al 2013).
We address this gap by conducting a comparative study in two semiarid ecosystems that are representative of the Sonoran and Chihuahuan Deserts (∼622 000 km 2 in North America). Both sites have transitioned from grasslands to shrublands as a consequence of woody plant encroachment (Van Auken 2000). Regional climate is characterized by a marked seasonality due to North American monsoon occurring from July to September, which accounts for 40%-70% of the annual precipitation (Adams and Comrie 1997) and leads to the greening of droughtdeciduous plants. Numerical studies have shown that the seasonal variation of climatic forcings and vegetation lead to complex soil moisture dynamics, with possible hysteretic patterns in the σ θ versus 〈θ〉 relation and whose physical controls change in time (e.g., Vivoni et al 2010, Yetemen et al 2015). Here, we use direct observations over a 9 month period from dense sensor networks installed at the two sites to: (i) investigate the field-scale spatiotemporal variability of θ, its seasonality in the σ θ versus 〈θ〉 space, and the occurrence of hysteresis; and (ii) identify the associated physical controls in the two ecosystems by means of an empirical orthogonal function (EOF) analysis.

Study area and datasets
The two study sites are located in the Santa Rita (SRER) and Jornada (JER) experimental ranges in southern Arizona and New Mexico, respectively ( figure 1(a)). Climate at SRER is characterized by hotter conditions than JER (mean annual temperature of 22°C and 18°C, respectively), mostly because of differences in elevation (1170 and 1470 m). SRER also receives a larger mean annual precipitation than JER (370 and 280 mm, respectively). As a consequence of these differences, vegetation is characterized as a woody savanna at SRER and a mixed shrubland at JER. Recently, an eddy covariance (EC) tower was installed at each site to measure hydrometeorological variables and surface turbulent fluxes, along with a distributed sampling network to capture the spatial variability of θ around the tower. The networks have a similar design consisting of twenty locations distributed over a fiveby-four grid (∼150 m by 120 m) with a regular spacing of ∼30 m and sensors at two depths (5 and 15 cm). Details on the instrument networks are provided by Anderson and Vivoni (2016).
Here, we use observations of θ, as well as precipitation (P) and evapotranspiration (ET), from 1 July 2013 to 31 March 2014 at 30 min resolution. Soil moisture data were averaged over the two sensor depths, thus representing the first 20 cm of soil. Missing θ data were reconstructed through temporal interpolation for gaps smaller than 12 h or by averaging concurrent data of the five neighboring locations with high correlations (always >0.84). To characterize local terrain, soil and vegetation, we computed the following factors in a 5 m by 5 m area centered on each sensor location: mean slope and curvature, derived from aerial imagery at 1 m resolution; percentage of gravel, sand, silt and clay obtained from soil samples collected at comparable depths; and percentage of the vegetation classes through the analysis of orthoimages at 50 cm resolution (grass, prickly pear, mesquite and bare soil at SRER; bare soil, grass, mesquite, creosote and other shrubs at JER) (Anderson and Vivoni 2016). As shown in figures 1(b) and (c), SRER is mostly characterized by sandy soil and grass vegetation, while JER is mainly bare soils with sand and gravel textures. More complete information on the land surface factors is reported in figures S1 and S2. Finally, 16-day composites of the normalized difference vegetation index (NDVI) from the moderate resolution imaging spectroradiometer (MODIS) were acquired in the two 250 m by 250 m pixels containing each network to quantify vegetation greening.

Methods
We investigated the relation between σ θ and 〈θ〉 by (i) dividing the study period into summer (July-September, 2013), fall (October-December, 2013), and winter (January-March, 2014) periods, with the aim of exploring differences due to the seasonality; and (ii) tracking the time evolution of the relation to identify the occurrence of hysteretic cycles. To identify the physical controls on the spatiotemporal variability, we conducted an EOF analysis (Hannachi et al 2007). This technique decomposes a space-time dataset into time-invariant orthogonal spatial patterns (the EOFs) and corresponding space-invariant time series of coefficients (the principal components, PCs) that control the importance of each EOF in time. The EOF analysis has been previously applied at catchment and regional scales using θ fields from instrument networks collected on a limited number of dates (up to 15 d) (e.g., Korres et al 2010, Busch et al 2012). Here, we focus on a smaller spatial extent (field scale) that has received much less attention (Vereecken et al 2014) and, for the first time, we use seasonally varying θ data covering a wider range of wetness conditions and climate forcing.
The EOF technique was applied as follows: (i) for each time step, we computed the spatial anomalies by subtracting the spatial mean from θ values observed at each location; (ii) we calculated eigenvalues and eigenvectors (the PCs) of the covariance matrix of the spatial anomalies; (iii) we projected the original anomaly dataset into the new space defined by the eigenvectors, obtaining the EOFs; and (iv) we identified the EOFs that explain most of the spatial variability and verified their statistically significance through the Kaiser's rule (Wilks 2006). Next, we provided a physical interpretation of the dominant EOFs by computing the correlation coefficients (CCs) between each EOF and the land surface factors previously defined. Finally, to explore how the importance of the EOFs (and the associated physical controls) varies in time, we calculated the percentage of the spatial variance explained by each EOF at each time step. For the kth EOF and jth time step, this is computed as: where s q j , 2 is the spatial variance of the soil moisture network at time j;s f k , 2 is the spatial variance of the kth EOF; and e jk is the PC of the kth EOF at time j. The derivation of equation (1) is described in text S1 in the supporting information. Additional mathematical details on the EOF analysis applied to θ fields are available in Jawson and Niemann (2007).

Results and discussion
4.1. Seasonality and Hysteresis in the relations between σ θ and 〈θ〉 Figures 2(a) and (d) present the seasonal relations between σ θ and 〈θ〉 at the two sites. To reduce noise, θ data at 30 min resolution have been aggregated at the daily scale. We first analyzed the overall behavior by fitting all points to the analytical equation of Famiglietti et al (2008). The general shape of the relations ('Overall Fit') is upward convex at SRER and asymptotically increasing at JER. Nevertheless, a closer inspection reveals the importance of seasonality in establishing the general shape at each site. In summer, a large number of storm events (N) and high values of ET lead to a large range of 〈θ〉, and a greater scatter in σ θ , as compared to fall and winter seasons (detailed values of N, P and ET are reported in table S1 of the supporting information). During the study period, the convex upward relation at SRER is caused primarily by a subset of fall and winter days with relatively high 〈θ〉 and low σ θ , whereas JER exhibits an asymptotic increase due to large number of summer days with high 〈θ〉 and high σ θ . These differences can be attributed to the high P at SRER (JER) during the fall (summer) season, exceeding the 75th (90th) percentile of the long-term records (table S1). This analysis reveals that using overall (single) relations obscures the superposition of complex spatial patterns that vary substantially across different seasons.
We further explored the θ spatiotemporal variability by tracking the time-evolution of the σ θ versus 〈θ〉 relations. In both networks and for each season, we identified the occurrence of consecutive hysteretic cycles. Specifically, we found six (seven) cycles at SRER (JER), with three (four) in summer, one (one) in a transition between summer and fall, and two (two) in the fall and winter (lumped as fall/winter). Figure 2 presents examples of the observed hysteresis in summer (b), (e) and fall/winter (c), (f) seasons for each site, with all hysteretic cycles reported in figures S3 and S4. Following Ji et al (2015), we identified three phases in each hysteretic cycle: (i) a wetting phase with sharp increases in 〈θ〉 and σ θ due to precipitation; (ii) a redistribution phase prompted by a precipitation hiatus of 2-5 d followed by a short sequence of storm/interstorm periods; and (iii) a dry-down phase with gradual decreases in 〈θ〉 and σ θ due to ET. For all hysteretic cycles, we computed statistics of their duration (D), average P and ET, and the change in 〈θ〉 (Δ〈θ〉) and σ θ (Δσ θ , table S2). While a large diversity of hysteretic cycles is observed, shorter durations occur in summer (∼15-17 d) as compared to fall/winter (∼65-68 d) at both sites. This is due to a large (small) number of storms and high (low) ET rates in summer (fall/winter) and their respective quick (slow) wetting and drydown phases. Wetting phases usually induce positive Δ〈θ〉 and Δσ θ , while dry-down phases have the opposite effect, thus closing the hysteretic cycles. Both clockwise and counterclockwise rotations are possible, depending on the spatiotemporal dynamics during redistribution phases. The redistribution (wetting) phase induces the largest changes in spatial variability (Δσ θ ) at SRER (JER), related to the sequence of storm/interstorm periods. These different pathways indicate that spatial variations during wetting and drydown are affected in varying ways by land surface factors, as explored next.

EOF analyses of θ and links to land surface factors
EOF analyses were applied on daily θ observations at 14 locations at each site with greater than 50% of nonmissing data. The first two EOFs (EOF1 and EOF2, figure S5) were found to be statistically significant (text S2) and together explain more than 80% of the variability of both datasets (table 1). As expected, EOF1 matches the spatial pattern of the time-averaged θ (figure S6), while EOF2 resembles the pattern of the time standard deviation of θ (the CCs are reported in table S3). Similar EOFs were also obtained with 30 min data, indicating that diurnal variations do not affect the dominant patterns. At SRER, EOF1 (76% of the explained variance) is primarily correlated with the  percentage of mesquite cover, while EOF2 (7%) does not have a statistically significant correlation with any factors. At JER, EOF1 (53%) is controlled by slope and sand content, while EOF2 (29%) is related to clay content. Table S3 details  Associated with the EOFs are PCs (e j1 and e j2 ) and their percentage of spatial variance explained (PV 1,j and PV 2,j ). Figure 3 relates these to time series of climatic forcings (P and ET), vegetation greenness (NDVI), 〈θ〉 and σ θ . As a result of the different shape of the σ θ versus 〈θ〉 relation, the strength of the correlation between the time series of e j1 and e j2 and those of 〈θ〉 and σ θ changes with the season (table S4). As noted earlier, the summer season is characterized by sharp changes in 〈θ〉 and σ θ due to storm events and high ET rates, while more gradual variations occur in fall and winter. At both sites, the PCs nicely show that θ variability results from the superposition of EOF1 with a stable contribution during the year (i.e., e j1 is positive with limited variability) and EOF2 whose influence depends strongly on seasonality (e.g., e j2 is positive in wet conditions and negative for dry periods). Sign changes in e j2 also occur after each storm event, suggesting an important role of EOF2 in the response to precipitation forcing. This is further illustrated through the relative importance of each EOF (PV 1,j and PV 2,j ), which can be linked to the land surface factors exhibiting high correlations to each EOF (table 1). For example, EOF1 is almost always dominant at both sites (i.e., PV 1,j larger than PV 2,j ), implying that spatial patterns of mesquite cover (SRER) and slope and sand content (JER) dictate soil moisture variability. Time periods when EOF2 increases in importance can be linked to interesting dynamics. For instance, storm events in all seasons lead to increases in PV 2,j followed by recessions during subsequent hiatuses, such that the effects of a mixed combination of land surface factors (SRER) and clay content (JER) on soil moisture variability operate during these periods (table 1). In addition, green-up periods with high NDVI in the summer and prolonged dry periods in the fall and winter seasons also induce a larger role for EOF2, indicating that soil moisture is impacted by vegetation seasonality.

Physical controls explaining the hysteresis between σ θ and 〈θ〉
The influence of land surface factors on the σ θ versus 〈θ〉 relations is analyzed in figure 4 using the proportion of the spatial variance explained by EOF2 relative to the total variance of the first two EOFs (PV 2,j /(PV 1,j +PV 2,j )). In both semiarid ecosystems, EOF1 generally dominates (blue colors) in the summer season for drier states (〈θ〉 less than 0.08 m 3 m −3 ), while EOF2 increases in importance (green to red colors) for wetter states (〈θ〉 greater than 0.08 m 3 m −3 ). Inspection of hysteretic cycles (figure 4 insets and mean values of PV 2,j /(PV 1,j +PV 2,j ) in table S3) for the summer season reveals that wetting (dry-down) phases are associated with higher (lower) contributions from EOF2, while the redistribution phase exhibits the maximum contribution of EOF2. Furthermore, under these summer conditions, the various physical controls on EOF1 and EOF2 act together to substantially increase σ θ , since both e j1 and e j2 are positive, as is clearly seen at JER during storm/interstorm sequences (i.e., high σ θ and 〈θ〉, dark red). In fall and winter, a remarkably different behavior is observed, EOF1 is generally more important at intermediate values of 〈θ〉 and the contributions of EOF2 increase towards the driest and wettest states, in particular at SRER. Based on the hysteretic cycles, the wetting and dry-down phases exhibit the influence of physical factors associated with EOF2 that act to substantially decrease σ θ , since e j2 is negative in these cases. At intermediate values of 〈θ〉, linked to redistribution (SRER) or a transition between wetting and drydown phases (JER), the dominance of EOF1 indicates that vegetation (soil and terrain) factors at SRER (JER) lead to high σ θ . These novel insights from the soil moisture observations are consistent with the numerical experiments of Fatichi et al (2015) indicating semiarid vegetation has a major control on σ θ , but that its sparseness and seasonality may lead to conditions where abiotic factors also become important.

Conclusions
Field-scale spatiotemporal variability of θ from dense observation networks in two semiarid ecosystems revealed that the overall shape of the σ θ versus 〈θ〉 relation results from the superposition of seasonal patterns and the relation evolves in time according to consecutive hysteretic cycles with shorter (longer) duration in summer (fall/winter). The observed spatiotemporal variability is largely explained by two dominant EOF patterns representing time-stable and seasonally varying contributions which act together to increase (reduce) the field-scale spatial variability in summer (fall/winter). Correlations between EOFs and land surface properties measured at a commensurate highresolution also indicated that θ patterns and hystereses are mainly controlled by vegetation (terrain and soil) factors at the site with higher (lower) vegetation cover. These findings (i) confirm mechanistic interpretations only previously available from modeling studies (e.g., Ivanov et al 2010, Fatichi et al 2015, and (ii) significantly advance the understanding of θ variability and its underlying physical controls at the field scale, which can aid soil moisture retrievals from multiple techniques. . Daily relations of σ θ and 〈θ〉 with the proportion of spatial variance (PV 1,j +PV 2,j ) explained by EOF2 (PV 2,j ) as colored symbols at (a) SRER and (b) JER for summer (left) and fall/winter (right). Insets show relations between EOFs and the hysteretic cycles of figure (D-D is dry-down, W is wetting, and R is redistribution).