Inference of Soil Hydrologic Parameters from Electronic Soil Moisture Records

Soil moisture is an important control on hydrologic function, as it governs vertical fluxes from and to the atmosphere, groundwater recharge and lateral fluxes through the soil. Historically, the traditional model parameters of saturation, field capacity and permanent wilting point have been determined by laboratory methods. This approach is challenged by issues of scale, boundary conditions and soil disturbance. We develop and compare four methods to determine values of field saturation, field capacity, plant extraction limit and initiation of plant water stress from long term in-situ monitoring records of TDR-measured volumetric water content (. The monitoring sites represent a range of soil textures, soil depths, effective precipitation and plant cover types in a semi-arid climate. The  records exhibit attractors (high frequency values) that correspond to field capacity and the plant extraction limit at both annual and longer time scales, but the field saturation values vary by year depending on seasonal wetness in the semi-arid setting. The analysis for five sites in two watersheds is supported by comparison to values determined by a common pedotransfer function and measured soil characteristic curves. Frozen soil is identified as a complicating factor for the analysis and users are cautioned to filter data by temperature, especially for near surface soils.


INTRODUCTION
Soil moisture is an important control on hydrologic function, as it governs vertical fluxes from and to the atmosphere, groundwater recharge, and lateral fluxes through the soil (Loik et al., 2004;Grayson et al., 2006;Vereecken et al., 2008). Soil moisture is also an excellent metric of hydrologic model performance, as it integrates temporal variation in precipitation and evaporation and is responsive to topography and soil physical properties governing fluxes of water (Wood et al., 1992;Cuenca et al., 1996;Rodriguez-Iturbe et al., 1999Laio et al., 2001;Western et al., 2004;Seneviratne et al., 2010;Bandara et al., 2013;Or et al., 2015;Paniconi and Putti, 2015). As such, the definition and quantification of physical properties that govern the ability of a soil to retain water against gravity drainage, evaporation, and transpiration have long been of interest to agronomists and physical scientists (Briggs and McLane, 1907;Buckingham, 1907). However, measurement of relevant soil physical properties has historically been conducted in laboratories, which may not represent the behavior of water in field conditions. Techniques to measure soil moisture in the field have improved dramatically in recent decades so that near-continuous measurements are now routine. In this paper we demonstrate how multiple years of nearcontinuous measurements of soil water content can improve estimates of hydrologically relevant soil properties.
The physical properties of soils determine several important thresholds in , defined as the volumetric fraction of water in a unit volume of soil (m 3 m −3 ): Soil saturation is conventionally defined as the condition when the entire pore volume is filled with water ( s ), which typically occurs only during subirrigation. The at which gravity drainage effectively ceases has been termed field capacity (FC; Veihmeyer and Hendrickson, 1949), and the at which transpiration ceases has been termed permanent wilting point (PWP; Briggs and Shantz, 1912). The difference in the value of FC and PWP is often used to estimate plant available water (Hansen et al., 1980) and the difference in the value of s and PWP is similarly used to estimate root zone storage (Seneviratne et al., 2010). Unfortunately, most measurements of these important field-scale properties are performed in laboratories on small samples.
The prevailing concepts of FC and PWP are a legacy of efforts to relate soil hydraulic and hydrologic parameters to standardized measurements of soil physical properties. Briggs and McLane (1907), introduced a "moisture equivalent" as a water to soil mass fraction for samples that were air dried, screened through a 2 mm sieve, packed in cylinders, saturated, and spun to achieve a force 1,000 times the force of gravity. These results were later expanded to analysis of PWP of wheat seedlings in a controlled study over a range of soil textures (Briggs and Shantz, 1912). Meinzer (1923) later defined "specific retention" as the volumetric fraction of water retained after long term drainage from saturation, which has been determined by measurement in small cells of disturbed soils (Hazen, 1892), and in situ (Ellis and Lee, 1919). Piper (1933) aptly concluded that experimental scale and depth of the capillary fringe exerted important controls on the results. Subsequently, Richards and Weaver (1944) built on the work of Buckingham (1907) to systematize measurement of matric potential ( )relationships. Using the pressure plate technique, they equated the at −33 kPa to the "moisture equivalent" for 71 near surface samples of irrigated soils (Richards and Weaver, 1944). The use of −10 kPa to determine at field capacity was reserved for coarse and volcanic soils. These advancements led to the development of capillary permeability relations for isotropic porous media (Burdine and Redford, 1952;Brooks and Corey, 1964) and continuous -conductivity models (Mualem, 1976) that remain in use.
The validity of these historic terms as physical constants has been widely criticized (Miller and McMurdie, 1953) resulting in refinements and additional terms for use in fields beyond agronomy. For example, it is understood that saturation may never be achieved in well drained soils under field conditions (Wang et al., 1998). Field saturation, fs , accounts for gas filled soil pore space that may be occluded within the pore matrix at approximately the bubbling pressure (Reynolds et al., 2002) and can occur due to air encapsulation (Fayer and Hillel, 1986), fluctuating groundwater levels (Marinas et al., 2013), and biogenic gases (Morse et al., 2015).
Field capacity has been criticized because the end point for gravity drainage is unclear and possibly method dependent due to differences in boundary conditions and scale (Colman, 1947;Hillel, 1980;Cassel and Nielsen, 1986;Kirkham, 2005). As a result, -values from −33 to −5 kPa have been proposed to define FC (Richards and Weaver, 1944;Salter and Haworth, 1961;Linsley and Franzini, 1972;Romano and Santini, 2002;Kirkham, 2005;Nemes et al., 2011). More recently, Assouline and Or (2014) proposed a method to relate the slope of the soil characteristic curve (n), the air entry value ( ae ) and the residual water of a soil ( r ) (van Genuchten, 1980) to a soil specific FC, based on the balance between capillary and gravitational forces. This dynamical approach seeks to transcend the use of static arbitrary values.
PWP has been criticized because the actual wilting point is species and soil texture dependent and, in the case of many native plants adapted to dry conditions, transpiration can cease with no wilting at in the range of −3 to −5 MPa for both tree and grass species (Scholes and Walker, 1993;Damesin and Rambal, 1995;Sperry et al., 2002;Rambal et al., 2003). To address these issues, Seyfried et al. (2009) used the term "plant extraction limit" (PEL) in place of PWP to describe the below which plants cannot extract water, and evaporation is the primary means of soil water loss. Therefore, at PEL is likely to be equal to or greater than at the limit of capillary continuity Assouline and Or, 2014), since evaporation can continue following the cessation of plant water use.
Furthermore, between FC and PEL, there is range of declining over which plant stress increases and evapotranspiration decreases (Rodriguez-Iturbe et al., 1999). This range of and the related range of is associated with a shift from energy-limited to soil moisture potential-limited evapotranspiration (Budyko, 1956;Koster et al., 2009). The threshold in water content between these states in terms of has been defined as s * Rodriguez-Iturbe et al., 2001), crit (Koster et al., 2004;Seneviratne et al., 2010), or ws (Smith et al., 2011) and conceptualized as an inflection point in evaporation-soil water content relation models. Similarly, Buckingham (1907) proposed that evaporation from soil is initially limited by moisture content, and secondly by vapor diffusion. Various formulations (Idso et al., 1974;Brutsaert and Chen, 1995) and methods have since been used to determine energy limited (Stage I) and supply/transport limited (Stage II) evaporation from soils (Salvucci, 1997). Stage III is vapor diffusion transport (Metzger and Tsotsas, 2005;Lehmann et al., 2008). During the second and third stages, the vaporization plane (e.g., drying front) will move downward from the soil surface as evaporation proceeds (Shokri et al., 2008;Shokri and Or, 2011). Therefore, the annual endpoint of soil drying may range from a value less than FC, to near the hygroscopic point, depending on the characteristics of local climate, soil, and plants .
The current capacity to resolve soil hydrologic properties and at scales sufficient to parameterize and validate distributed hydrologic models remains challenged by the relatively limited data available to define spatial heterogeneity and effective scales of measurement (Njoku et al., 2003;Reichle et al., 2007;Huang et al., 2016). Since transitions in fluxes in agronomic and ecohydrologic models are often indicated by FC, ws , and PEL either indirectly (Liang et al., 1996) or directly (e.g., Wigmosta et al., 1994;Boote et al., 2008;Seneviratne et al., 2010;Paniconi and Putti, 2015) and are bounded by s (or fs ) and either PEL or some lower limit to , we consider these important variables for hydrologic modeling. Recent approaches to determine soil moisture parameters include inverse modeling, remote sensing and other synthetic approaches (Santanello et al., 2007;Montzka et al., 2011;Bandara et al., 2013). However, soil hydrologic properties would ideally be determined by field measurement (Romano and Santini, 2002). Hillel (1980) recommended that properties such as FC must be measured repeatedly in the field, after the entire depth of the soil profile is wetted and that laboratory methods based on equivalent to −1.5 MPa, −33 kPa, or −10 kPa matric potential are arbitrary static measurements that are not representative of systems with external controls on the soil boundary conditions, such as poor drainage or evaporation.
Given the increasing availability of electronic measurements of soil water content and considerable theoretical and conceptual controversy regarding soil water content constants, the primary objective of this paper is to determine if temporal patterns of measured soil water content values are consistent with the concepts of s , FC, and PEL. The proposed methods presented here are intended to determine consistent and physically/ecologically relevant estimates of these soil hydrologic parameters directly from long term field measurements of . We pose the hypothesis: Data attractor values and inflections in time series records represent changes in dQ/dt that are consistent with the concepts of FC, ws , and PEL. This is based on prior observations that FC and PEL exist as quasistable states (Miller and Klute, 1967;Cassel and Nielsen, 1986;Kirkham, 2005) following periods of drainage and evaporation. These states are represented by greater data frequency (data attractors) at intermediate and low values of . That is, following input events, rapidly trends toward FC where it is relatively stable due to capillary forces until decreases by evapotranspiration or increased by input flux greater than the associated unsaturated hydraulic conductivity [K( )]. Then, as soils dry they stabilize near PEL. Thus, FC may be better represented in the absence of strong ET, such as during winter in mid-latitude regions, and PEL is best determined when energy does not limit transpiration. To test the hypothesis, we first compare soil hydrologic property values determined by four methods of analysis to and then assess the accuracy of the estimated values to values determined from a range physical methods at traditional soil water potentials. Finally, we evaluate the applicability of the tested approaches to other environments.

Setting
The analysis exploits the wide range and regular seasonal changes in soil moisture at sites near Boise ID, USA. In semiarid Mediterranean climates, seasonal patterns in soil moisture storage and streamflow follow seasonal changes in precipitation, temperature, and solar radiation (McNamara et al., 2005: Flerchinger et al., 2010Seyfried et al., 2011). In this system, the water year begins with dry soil following the annual summer drought. As autumn rains arrive, the soil wets downward from the soil surface and soil profile water storage increases over a long wet period due to winter rains or snow melt. In spring, warming temperatures drive snowmelt and the annual hydrograph. Ultimately, precipitation decreases and plant water use depletes soil moisture as summer proceeds. In many years, this hydrologic setting imposes different soil surface boundary conditions in each season; intermittent surface wetting of the dry soil in fall, low flux snow melt in winter, high flux snow melt in spring and progressive soil drying over the summer (Figure 1).  Seyfried et al. (2001) and Chauvin et al. (2011). Data for DCEW are available at http://earth.boisestate.edu/drycreek/.

Sites
The Upper Sheep Creek watershed is located within the Reynolds Creek Experimental Watershed, in the Owyhee Mountains of southwestern Idaho. The site ranges in elevation from 1,840 to 2,036 mamsl and is underlain by basalt bedrock. Mean annual precipitation is 426 mm, about 60% of which is snow. The snow is redistributed by wind, resulting in snow accumulations of 3-4 m in drifts and <15 cm in scour zones. Deposition patterns are persistent and result in distinctive snowsoil-vegetation complexes (Flerchinger et al., 1998;Prasad et al., 2001;Chauvin et al., 2011). The low sagebrush complex (Low Sage) is found on shallow (<50 cm deep) soils on mostly westfacing slopes and ridges. Soil textures are loam in the upper 10-20 cm and clay loam and gravel in the subsoil below the argillic horizon. The mountain sagebrush complex (High Sage) is on deeper soils, about 1 m thick, and mostly on north-facing slopes. The aspen complex (Aspen) is on deep (>2 m) soils and are closely associated with snow drifts. Soils in the High Sage and Aspen complexes were formed in aeolian deposits and are highly uniform with depth, mostly categorized as silt loam. Upper Sheep Creek was burned by a prescribed fire in summer 2007. The Big Sage area burned completely, and was replaced with grasses and forbs by 2008. The Aspen was cut down and recovered, in a hydrologic sense, by 2009. The Treeline and Lower Weather sites are located in the Dry Creek Experimental Watershed (DCEW) in the foothills of the Rocky Mountains of southwestern Idaho. The elevation of these sites is 1,620 mamsl at Treeline and 1,150 mamsl at Lower Weather. Both sites are underlain by weathered granite from the Idaho Batholith, although the Lower Weather site is below the shoreline of prehistoric Lake Idaho. Mean annual precipitation is 520 mm at Treeline and 300 mm at Lower Weather. Over the period of measurement, snow accounted for ∼50% of precipitation at Treeline (Kormos et al., 2014) but Lower Weather receives mostly rain. Similar to Upper Sheep Creek, geologic processes including wind scour have formed deeper soils on north facing slopes at Treeline and snow accumulation up to 2 m currently enhances springtime soil water where the sensors for this study are located. The vegetation at the site is typical of a transition between lower elevation grass-lands and higher elevation forests; mountain big sagebrush and ceanothus shrubs, prunus subspecies, forbs, and grasses. The Lower Weather site is not subject to drifting due to the limited snowfall at the lower elevation. Soil texture at DCEW is sandy loam and ranges in depth up to 1.3 m (Gribb et al., 2009;Tesfa et al., 2009).

Data Collection and Conditioning
Soil water content data were collected from representative sites of each of the five soil/vegetation complexes. These sites represent various controls of soil texture, soil depth, vegetation, and annual precipitation (Flerchinger et al., 2010) on (Figure 1). In each instance, a pit was hand-excavated to bedrock or saprolite and two parallel profiles of soil moisture instruments, separated by 2 m, were installed horizontally into the pit face. A summary of sensor depths, soil texture and mean annual precipitation or estimated annual soil water input is provided in Table 1. At Upper Sheep Creek, hourly data were collected with time domain reflectometers (TDR 100, Campbell Sci., Logan UT, USA), with 3-prong, 30 cm long rods. Data are reported on water year basis. At Dry Creek, hourly data were collected with frequency domain reflectometers (CS615, Campbell Sci., Logan UT, USA) calibrated with less frequent samples from co-located TDR waveguides (Chandler et al., 2004). Data were prepared for analysis by removing most electronic noise, out of range data or data from failed sensors (Nayak et al., 2008;Kormos et al., 2014) and are reported by calendar year. Data cleaning did not include removal of records during periods of frozen soil, which were most prominent within 15 cm of the soil surface and affected a limited number of sensors. Frequency distribution of from single sensors was done by constructing histograms from the conditioned data. The frequency analysis of temporally coincident data from paired sensors was done by binning data in two dimensional frequency arrays with increments of 0.01 for 0.0 < < 0.65.
Limited laboratory and field analyses are used to relate the results of the developed methods to traditional capillary pressure-saturation techniques. Measurements from a 66 kPa field tensiometer co-located with the 60 cm TDR sensor at Aspen site represent the range from field-saturated conditions just after snowmelt on May 20, 2004 to the time of removal on August 8, 2004. Laboratory measurements of over the ranges of 1.4-430 kPa (HYPROP, UMS GMBH, Munich Germany), and and (corrected for bulk density) 0.2-115 MPa (dewpoint psychometer, Decagon Devices, Pullman WA, USA) and were made to construct characteristic curves for soil samples collected at depths of 0-3 and 41-43 cm. At Treeline, soil water potential was previously measured in a field experiment to determine FC by surface irrigation from dry initial conditions. The experiment was interrupted by rain during the drainage period and terminated when the tensiometer lost capillary connection to the soil at ∼33 kPa. Soil water potential was monitored by automated tensiometers (Gribb et al., 2009). Values of measured during the natural cycle of wetting and drying provide the basic data used in this analysis. In all cases the soils are well drained, therefore there is no attractor at or near soil saturation. We define attractors as values of that occur at high frequency and develop four analyses that are based on observations made in a climate with distinct annual wet and dry phases. Field Saturation is the maximum observed value of . Since complete soil saturation is infrequent in well drained field soils, the maximum annual value of fs is likely to range between FC and s . Field Capacity is a moisture state that prevails under wet, low flux conditions. In well-drained conditions, saturated soils drain quickly to near FC, then more slowly. Where soil water potential ( ) is seasonally affected by plant water uptake, FC may vary slightly with diel changes in evapotranspiration. Conversely, dry soils wet to a somewhat greater than FC before drainage begins. These two processes enhance the local data frequency near FC and result in a data attractor that we define as the "wet attractor." Figure 2A shows a wet attractor near of 0.20 throughout the winter months when drainage is minimal and transpiration is insignificant. We propose that this attractor is approximately equivalent to FC. Several other descriptive and mechanistic definitions of FC have been developed and are presented by Assouline and Or (2014), who suggest the stability of soil water retention is a function of disruption of liquid phase continuity. Plant Extraction Limit is the at which existing vegetation does not extract soil water by transpiration. Conceptually similar to the PWP, PEL allows for variations among soil-plant combinations and acknowledges that the vegetation doesn't necessarily wilt when transpiration ceases (or becomes very slow). PEL is expected to be an attractor in environments with extended dry periods and perennial vegetation. Under these conditions, declines during the dry The counterclockwise VWC hysteresis trace shows distinct periods of surface wetting and soil profile wetting for 15 and 100 cm sensors and two inflections in the drying limb. crit at 10 cm depth is determined as the initiation of declining dQ/dt. Panels (A,C,D) display 15-min data as semitransparent markers to display data attractors by superposition. "Wet" data attractors are identified by ovals in (C,D).
season until PEL is reached, then remains unchanged until replenished with new rainfall or snowmelt except near the soil surface, which may be subject to direct evaporation and may proceed to values less than the PEL. These processes result in a "dry attractor" in the data. Plant Water Stress Initiation ( ws ) is not an attractor, but is a decrease in plant water use as governed by complex physiological traits that control stomatal aperture (Osakabe et al., 2014), and is an inflection in soil moisture records (Rodriguez-Iturbe et al., 1999;Smith et al., 2011) determined as d 2 Q/dt2 > 0.

Analysis Methods
Four methods of analysis to determine the parameters PEL, FC, and fs are presented. Throughout the remainder of the paper, these analyses will be referred to as time series, frequency, paired, and hysteresis, respectively. Visual inspection of records in time series is a simple approach to identify approximate values for each parameter for individual soil depths. Figure 2A shows this approach for one year of data from a 15 cm sensor at Treeline. At this site the minimal overlap between periods of slow drainage in winter and consistently low in spring and summer. This allows subjective approximation of the wet and dry attractors, and clearly shows the maximum annual , and an inflection in (t) indicating a reduction in plant water use.
Frequency analysis of distributions for single sensors show modal peaks in at dry and wet soil moisture attractors for water year ( Figure 2B) or period of record data. The maximum value of is assumed to be equivalent to fs . Paired sensors at similar depth, are often used to determine an average value of , increase spatial representation of the measurement, or provide redundant data in case of instrument failure. These synchronous data can be plotted as x-y coordinates to determine high frequency attractors and extreme values. Figure 2C shows example data for a water year from paired sensors at 15 cm depth in soil pits separated by 2 m. Ideally, data pairs of from paired sensors at similar depth will fall on a 1:1 line, bounded by the minimum and maximum at the sensor depth. However, differences in the timing of wetting and drying result in divergence and spreading around a 1:1 line. In the example for one water year of data, attractors emerge at dry (≈0.05) and wet (0.19-0.24) attractor values, respectively and the greatest measured values of (0.30, 0.27) for fs on the respective axes of each sensor. Although, similar to the frequency analysis, the paired sensor analysis leverages twice as many sample locations per data point. The use of data pairs decreases the influence of erroneous or extreme values on the frequency of values near the data attractors on the central trend. Additionally, the resulting variance around the trend may provide the user greater guidance on the range of local variability in soil hydrologic properties than would a single sample location.
Hysteresis is often observed in the temporal behavior of over seasonal wetting and drying periods for sensors at two depths in a vertical profile. The lag in both wetting and drying periods between the two depths enhances the relative frequency of wet and dry attractors, near the intersection of the wetting and drying limbs of the hysteresis loop. For example, Figure 2D shows hysteresis between at 15 and 100 cm at Treeline. For the climate of this study, the water year begins near the end of the dry season, with at near minimum values for both the surface and bottom of the soil profile. Accordingly, the trace for a water year begins near the low attractor. As autumn precipitation increases at 15 cm depth, the at 100 cm depth remains nearly constant until the wetting front arrives. For a water year at the example site, the horizontal wetting trace provides strong visual support for the interpretation of low attractor from the small range in the ordinate. Along the right side of Figure 2D, wetting at the deep sensor is shown as the vertical trace up to the wet attractor, found at the intersection of the wetting and drying traces and identified by the black oval. As above, the maximum values on either axis represent fs . Finally, the drying limb terminates at the dry attractor and provides visual support for identification of the dry attractor for depth of the shallow sensor.

RESULTS
We found wet and dry attractors within the expected ranges of FC and PEL by each of the presented methods, equated fs to the maximum value of and interpreted the inflection during the drying phase, d 2 Q/dt2 as ws . Results from the time series, frequency, paired sensors, and hysteresis approaches are in general agreement, although each has a different bias. Below, we first present detailed graphical results from the paired and hysteresis analyses to demonstrate the attractors for the Aspen, Treeline, and Low Sage sites, which span the range of soil texture and climate among the study sites. Then we compare values determined by these methods at three depths for each site with results from the time series and frequency analyses. Example data from soil physics experiments performed in the laboratory and from in situ measurements at Aspen and Treeline sites are then provided to relate the extreme values and data attractors to PEL, FC, and fs . Figure 3 shows two visualizations of paired sensor analyses over the study period for Aspen (90 cm), Low Sage (10 cm), and Treeline (15 cm). Here, and below, we present data that demonstrate analyses over the broadest range in soil texture, soil depth, and vegetation class. For the Aspen the top panel shows semi-transparent data points with makers to represent the data attractors and the bottom panel shows color maps of logtransformed data frequency values. Many of the paired sensors at a depth, such as Aspen, 90 cm show a strong linear correlation between sensors (central trend) with occasional hysteresis traces outside of the central trend. Shallow paired sensors such as at Low Sage and Treeline (Figure 3) show much greater variability in around the central trend. These features result from the temporally local hysteresis between sensors during input events, repeated over the annual period of changes in .

Paired Sensors
Comparison of plots of raw data in the upper panel of Figure 3 with the log transformed frequency colormaps in the lower panel demonstrates the greater visual support of the transformed data for identification of attractor values. Whereas, the point values can be selected directly from the frequency value matrix, as was done to identify attractors for the upper panel, the colormaps simplify the interpretation of the frequency value matrix. In Figure 3 (lower panel) dry attractors for the period of record are clearly shown as dark red points in the colormaps, often near the minimum , and are generally surrounded by a set of high frequency values which may include two or three modes. Similarly, the period of record wet attractor value is often diffuse, with weak frequency modes for the Aspen (silt loam) and Low Sage (clay loam) sites, but with a more distinct attractor for the Treeline (sandy) site. Variability between sensors extends to fs , which varies more by depth and year than the other paired analyses because it is an extreme value. The minimum and maximum of annual estimates (Table 2) bound the period of record estimates, and provide an approximation of the range in variability of the determined values. For the Treeline site, the range in annual estimates for individual sensors within the paired sensor configuration was zero to 0.02 for both the attractors and fs , likely due to the high sand content at this site. At Low Sage and Aspen, the range increased progressively to 0.01-0.05 and 0-0.11 , respectively, primarily due to annual variability in at depth. Spatial differences in (t) resulted in variance from unity slope (1.00) for many sites and depths, but were similar within soil texture class: For the Aspen and High Sage sites the slope of the relation between paired sensors increased from a slope of 1.00-1.03 near the soil surface to as much as 1.33 at depth. At Low sage the slope of the relationship ranged from 1.10 to 1.19 from 10 to 40 cm depth, and 1.36 at 50 cm, in the C horizon, which is very stony. The Treeline and Lower Weather sites have low bulk density surface soil, with frequent macropores the slope of the trend at 5 cm depth (1.17-1.29) is much greater than for the deeper, stony soils (1.12).
Hysteresis plots of for period of study data at Aspen, Low Sage, and Treeline show more complex wet and dry data attractors than are shown for an annual cycle in Figure 2D. Figure 4 presents visualizations of hysteresis  between shallow and mid-profile depths, mid-profile and deep sensors, and log transformed colormaps of selected plots. Despite differences in soil depth, texture, annual water input from rain and snow, and vegetation cover types, the hysteresis plots for most sites are similar in form, with a few important differences. All hysteresis plots with a near surface measurement (e.g., Figures 2D, 4A-C) show a marked inflection in the drying limb, whereas hysteresis plots for lower depths (e.g. Figures 4D,E) are more triangular, with the constant slope of the drying limb extending to near the dry attractor, indicating more uniform drying across depth. For some sites (e.g., Figure 4F) the hysteresis plot can vary annually between these forms, depending on the timing of precipitation. As a result, the range of the wet attractor over the depth of the soil profile is from 0.13 to 0.08 at Aspen and less (0.02-0.07 ) at drier sites such as Low Sage and Treeline (Table 3). Similarly, the distribution of frequency values within the hysteresis loop is quite different among Aspen, Low Sage and Treeline (Figures 4G-I). The differences in the distribution frequency is directly related to differences in depth and timing of precipitation and snowmelt and soil drainage rates for the different soil textures. Figure 5 compares the range of estimated values for fs (Figure 5A), and the wet ( Figure 5B) and dry ( Figure 5C) attractors for selected soil depths at the study sites. Data are presented as maxima and minima of all estimates for time series, frequency, paired, and hysteresis analyses. Representative estimated values for each parameter were determined by comparison of the average of annual estimated values to period of record estimated values. fs values are equated to the maximum found by long term (decadal) frequency analysis ( Figure 5A) and exceed average values from the paired and hysteresis analyses by 0.02-0.16 . For the wet and dry attractors, value estimates from the paired and hysteresis analyses are often similar, despite differences in range, both of which are narrower than the value ranges for time series and frequency analyses (Figures 5B,C). The representative estimated value was determined as the most common value among period of record estimate, the average value of annual estimates of paired analyses, and the average value of six period of record hysteresis estimates. This approach generally resulted in a set of values which differed by <0.02 . Context for the results in Figure 5 is provided by estimates at 0 kPa, 33 kPa, and 1.5 MPa for the various soil textures by FIGURE 4 | Example plots of data from sensors at different depths in one soil pit each at Aspen and Low Sage and for two pits at Treeline, using a range of sensor depth combinations. Data are shown with the deeper sensor on the y-axis to maintain a counterclockwise hysteresis for all plots. Data markers are semi-transparent to highlight high frequency areas by superposition. Circular data markers are placed at points of maximum frequency corresponding to PEL and FC and maximum values at fs . As in Figure 3, the highest frequency data values are displayed as circular markers for the dry and wet attractors and fs . Frames (B,D) are shown as colormaps in frames (G,H). Plot (I) at Treeline pit 4 is displayed for the 30 cm depth analysis due to a failed sensor in pit 3 to compare with pit 3 data in (C,F). Additional data overlays of field tensiometer data (A,C) support later discussion. a commonly used pedotransfer function (Saxton and Rawls, 2006).
The inherent uncertainty in the determining representative attractor values is complicated by physical heterogeneity in soil properties and the somewhat subjective assignment of a single value from within a data attractor that may be weak, have multiple maxima (Figure 4G) or a wide value range ( Figure 4H). Figure 6 provides a summary of data attractor results from Figure 5 with the addition of estimated point values of ws , determined as single values from time series analysis. The ws values ranged from 0.07 to 0.37 among sites, with the values from Treeline and Lower Weather ranging from 0.06 to 0.10, slightly less than or equal to the values determined by Smith et al. (2011) for sites in Dry Creek watershed. The ws values ranged from a seven to 38% of the difference between PEL and FC. Table 4 compares estimated values of wet and dry attractors to data from soil water retention curves, field tensiometry and an intensive soil irrigation, and drainage experiment. We found that the example dry and wet attractors for Aspen and Treeline sites are equivalent to FC and PEL, respectively. Figure 7 shows the soil water retention curves for Aspen for samples from 1-3 cm (blue) and 41-43 cm (brown) and field tensiometer data at 60 cm (black). The wet attractor values for both depths at Aspen (0.41, 0.41) nearly matched the value (0.42, 0.42) at 10 kPa, commonly used as an approximation of FC. Similarly, the dry attractor values (0.13, 0.13) for both depths at Aspen nearly matched the 1.5 MPa values (0.12, 0.12). For Treeline at 60 cm the hysteresis method wet attractor (0.21) compares well with the measured (0.19) corresponding to a tensiometer value of 10 kPa. Although, not exhaustive, this evidence strongly supports the hypothesis that data attractors represent FC and PEL, and that other important soil hydrologic parameters can be determined from records for sites with different soils, wetness and plant cover.

Relation of Inferred Values to Physical Measurements
Additional inspection of Figure 7 and results from the frequency analysis of are presented here to clarify the interpretation of s , and fs . For the Aspen the laboratory saturated (0 kPa) ranges from s of 0.64 at 1-3 cm to 0.50 at 41-43 cm depth. The frequency method estimated value of fs at 10 cm depth (0.55) falls between these values and is substantially less than s for the soil surface measurement. This difference is likely due to either a decrease in soil porosity between 2 and 10 cm depth or a difference between at s and fs . At greater depth the difference between the frequency method value for fs at 60 cm (0.51) and the value at for 0 kPa at 41-43 cm (0.50) is negligible, indicating fs is equivalent to s at these depths. For sandy soils at Treeline, constant irrigation achieved s at 60 cm

DISCUSSION
The objectives of this paper are to evaluate four approaches to determine if temporal patterns of measured soil water values are consistent with the concepts of s , FC, PEL, and ws . The method exploits the increasing availability of electronic soil water  Saxton and Rawls (2006).  FIGURE 7 | Soil moisture characteristic curve measurements for Aspen site (60 cm), samples from similar silt loam sites in RCEW, and sandy loam at Treeline site (60 cm). Soil water tension values are shown as pF, (log hPa) and identified on the curves at 10, 33, and 1,500 kPa as vertical gray lines in both panels. Measurements for the silt loam include Aspen site 60 cm field tensiometer (black) and laboratory measurement by psychrometer (violet) and HYPROP at 0-2 cm (brown) and 41-43 cm (blue) for silt loam samples from a similar site nearby Aspen. Measurements for Treeline site sandy loam are shown for a prior experiment to determine FC that was interrupted by a period of rain. Laboratory and model estimates for r are soil samples from 24 to 25 cm depth (Gribb et al., 2009). Values determined by frequency analyses are shown as black dashed lines in both panels. A range of relatively high frequency data greater than FC are shown to indicate the range of estimates in this value by various methods.
data to address the considerable controversy, both theoretical and conceptual, regarding soil water content constants. We found good agreement between soil hydrologic parameter values determined from attractors found in in situ sensor records and from conventional laboratory techniques, for a limited number of samples representing sandy loam and silt loam soils (Figure 7). These results support the hypothesis that data attractors represent FC and PEL, and that other important soil hydrologic parameters can be determined from records for sites with different soils, wetness, and plant cover. We demonstrate the results of the analyses by overlaying the representative parameter values on the time series data from which the values were derived (Figure 8) to support discussion of the determination of values for various soil hydrologic properties, applicability of this approach to other sites and caveats for use where frozen soils occur. In the following sections we comment on the value of the approaches for estimating each property ( s , FC, PEL, and ws ), followed by a discussion of uncertainty and errors.

Saturation ( s )
Saturation was uncommon in our data. All of the study sites are well-drained and did not show evidence of a water table. Therefore, it is not surprising that we found no attractor around saturation. However, near surface soil pipes and macropores may intermittently fill and appear as data outliers. This may explain why the maximum measured values of can be quite different from the pedotransfer function estimates of saturation, especially near the soil surface. Results from our analyses for Aspen and Treeline sites show maximum ( fs ) values of 0.55, 0.51, and 0.34 (Table 4) from controlled experiments and values of 0.64, 0.50, and 0.39 from field measurements. The difference in peak saturation for shallow soil at Aspen may arise from differences in porosity between the surface 3 cm and at 10 cm depth, or the greater capacity for occluded soil gas below the very dense mat of fine roots within above the mineral soil. We note that complete saturation of soils requires extensive irrigation, either in the field (Ellis and Lee, 1919;Hillel et al., 1972) or laboratory (Hazen, 1892;Watson, 1966) to ensure near complete removal of gas from the soil pores. Analysis of both laboratory and field measurements should be repeated for multiple depths in the soil profile.
At most sites, field saturation, fs , is a transient state achieved during infrequent periods of high flux and can range from FC to s in annual records (Figure 8). Assuming a consistent pore size distribution over time, it is logical to select the greatest recorded in situ value as fs either as a tail of the frequency analysis over the period of record. This selection requires judgment to determine if the extreme value is a valid measurement point or an artifact of measurement. This is likely an appropriate proxy value ( Figure 5A) for most modeling purposes as long as the site remains freely draining. fs decreases with depth for all sites except Low Sage (which had a strong textural contrast), presumably in response to decreased porosity, decreased organic matter increased occluded volume by gravel or attenuation of flux by storage. These effects are lumped by the approach presented here. Due to the dependence of on flux for well-drained sites, determination of fs is likely improved with longer data series and shorter time step data, especially for semi-arid or drier sites.

Field Capacity (FC)
The most consistent representation of FC often corresponded with the intersection of the wetting and drying limbs of the hysteresis loops ( Figure 2D) and near the low edge of the wet data attractor in the paired analysis, as shown for Aspen and Treeline (Figure 3). These constraints on selection of FC recognize that this attractor requires that soil wetting is followed by a period of several days with no evaporation. Therefore, FC is most evident after cool season precipitation events, which are dominant in the study environment. During snowmelt, radiation and temperature patterns drive diel fluctuations of melt water flux to the soil that can to lead to positive bias in estimation of FC.
We found very similar values for FC at all depths within Aspen, Treeline, and Lower Weather but variability in FC with depth at High Sage and Low Sage (Figure 6). Although, texture is a first-order control on FC, the values from our analyses were often quite different from the pedotransfer function predictions, and even the same texture: FC at High Sage ranges from 74 to 82% of FC at Aspen for the surface and deep sensors, respectively, which is remarkable, given the close proximity and nearly identical soils at these sites. The difference in attractor values between Aspen and High Sage is due to the different soil water input environments. Aspen site is located beneath a seasonal snow drift that is typically 3 m deep. This results in extended periods during which daily water inputs from snowmelt maintain soil water contents at values greater than FC. Winter snow cover at High Sage site is typically 50 cm. It therefore melts much earlier and is subject to periodic cool season rainfall events that allow for drainage to FC, thus shifting the attractor downward to a value more consistent with the concept of FC. The pedotransfer function predictions appear reasonable (with hindsight) for sites other than Aspen (Figure 6), but appear to correspond better to fs values presented for Low Sage and Treeline. These differences in representation of the mobile water fraction in soil are likely attributable to soil plant interactions such as litter cover and secondary soil structure and are important for careful parameterization for ecohydrological models.

Plant Extraction Limit (PEL)
PEL, like ws is an ecohydrologic parameter that represents an interaction among plant phenology, the surface energy balance, and soil water potential. For this study, the land surface cover is dominated by perennial vegetation with varying extent of litter and bare soil across sites, soil water input is primarily derived from spring snowmelt resulting in an extended dry growing season, and peak PET is in July . These conditions, in conjunction with the low mean annual precipitation outside of snow drifts all facilitate very dry soil states required for development of an attractor at PEL. The dry attractor values in this study were generally slightly less than those predicted by the pedotransfer function, and converse from fs , tended to increase with depth. Our initial premise was that PEL at any depth is represented by a constant value of for an extended period during the growing season, as shown for Treeline and Low Sage (Figure 8).
For other sites, PEL may not be consistently represented across years for either the shallow or deep sensors for two reasons: First, Aspen is often energy limited and retains water at depth in most years, only achieving PEL in 2003, 2007, and 2012 (Figure 8).
Similarly, other sites with insufficient energy to evaporate the mean annual precipitation are unlikely to reach PEL at depth. Second, less than PEL commonly occurs in near surface soils. Records from Treeline and Lower Weather show as low as 0.01 for sensors at 5 cm depth, but a PEL of 0.06 to 0.08 at depth (Figure 8). Evaporative drying near the soil surface can decrease below PEL. In this case, when the evaporative potential at the leaf surface is balanced or exceeded by soil water tension, plant water use is negligible and vapor diffusion governs evaporation (Buckingham, 1907;Salvucci, 1997). This transition from viscous capillary (S1) to vapor diffusion (S2) control on evaporation (Brutsaert and Chen, 1995;Lehmann et al., 2008) is apparent in data, and complicates the interpretation of the dry attractor and PEL.
Thus, interpretation of the range of the dry attractor requires consideration, as with the wet attractor. Lower Weather receives the greatest energy and least precipitation among the study sites, and has coarse textured soils with a nearly bare surface, making it a good case study for dry soils. Initial estimates of the dry attractor range from 0.04 at 15 cm to 0.08 at 100 cm (Figure 8). The much lower values of the dry attractor for 5 cm sensors at 0.00 < < 0.02 (Figure 8) likely represents the limit of evaporation via vapor diffusion near the surface, for the available seasonal energy. We found a gradient 0.00 < < 0.08 develops from 5 to 100 cm depth in response to the vapor diffusion gradient above PEL at 100 cm depth. These observations lead to the conclusion that near the soil surface, the shift from evapotranspiration to evaporation below PEL is difficult to distinguish using a data attractor.

Plant Water Stress ( ws )
Important features of soil drying below FC include the initiation of plant water stress, ws , the endpoint of plant water use, PEL, and the transition to evaporative drying by vapor diffusion. We found that the semi-arid climate of the study sites provided ideal conditions for soil drying through spring and summer, as shown by a nearly constant drying trend over approximately the middle half of the range in for all sites (except Low Sage). Yet, a notable difference is evident between drying patterns for shallow and deep sensor locations: Soils at depths of 15 cm or less dried first and fastest at near steady state from below FC to near PEL (Figure 8). For soils below 15 cm the drying function (t) was quasi-sinusoidal between FC and PEL, with a shift to faster drying when the shallow sensor reached ws . We attribute this initial increase in the rate of deep soil drying to consistent (energylimited) plant water use from a smaller storage volume following surface drying below ws . For all sites, ws occurred annually for the shallow sensor depths, but was less consistent for the deep sensors. For example, records from the deep sensors Aspen and High Sage regularly show linear declines until approximately Oct 1, when the water year ends with the onset of freezing nighttime air temperatures and abrupt termination of plant water use at these high elevation sites, e.g., WY 2008-2010 (Figure 8). Thus, as s may not be consistently achieved due to surface infiltration rates that seldom exceed soil profile drainage flux, ws may not be consistently achieved due to insufficient seasonal energy. Yet, the occasional dry years (e.g., WY 2006) at Aspen and High Sage are clearly water limited and are sufficient to estimate ws for the high elevation sites. Low Sage and the lower, elevation Treeline and Lower Weather sites are annually water limited and show consistent inflections indicating ws during July and August with values ranging from 0.20 for the deep clay soils at Low Sage to 0.07 for deep soils at Treeline (Figures 6, 8).

Uncertainty and Errors
One approach to the uncertainty in evaluating contemporaneous soil moisture states is to visualize (z) at a daily to weekly time step. Figure 9 shows (z) for wetting and drying periods over one water year at Aspen and Treeline. Although, this approach is somewhat difficult to implement over multiple years, it provides insight on the gradation of soil hydrologic properties with depth and the different dynamics of wetting and drying across sites. This type of diagram is common in texts, but typically focuses attention on the dynamics of wetting and drying, rather than the predominant soil moisture states. For both sites, a vertical profile at FC emerges after the wetting front arrives at the base of the soil. Yet (z) may not reach fs , as determined from the long term record for any depth in a given year. During soil drying, similar increasing gradients in (z) emerge slightly below FC at both sites in response to plant water use, and indicate a rapid decline in drying associated with ws , near PEL. The much sharper inflection in (z) near the soil surface at Treeline than Aspen reflects the difference in surface evaporation between sites. Finally, the much broader band of overlap among (z) profiles at Treeline reflects the longer period and greater depth of plant water stress at that site.
Each approach to determining soil hydrologic properties has a different bias, but all are in general agreement. The applicability of measured data for this model parameterization depends on the nature of soil water dynamics under consideration. Selection of specific values for the wet and dry attractors by time series analysis is more subjective than by frequency analysis, as it depends on selection of a single value from a range of near the extreme values. The two depth hysteresis method provides estimates that closely match the at standard soil water potential values commonly used to define FC and PEL. The improvement over the other methods for determining these wet and dry attractors is related to the reduced overlap between valid data, and the invalid data and data errors which can skew the frequency distribution of the attractor. In particular, the juxtaposition of for a sensor depth (a) with stable to a sensor at a different depth (b) with either increasing or decreasing develops a linear trace toward a point of inflection where the rate of change in sensor a increases and sensor b decreases. The relatively slow rate of change near the inflection enhances the frequency of (a,b) data at the attractor. Error values tend to fall at the margins of the hysteresis loop which reduces error distortion of the main attractors (e.g., Figure 4). Unlike the hysteresis analysis, the paired sensor approach can also be biased by any temporal lag in wetting front depth between the sensors. Soil depth and cover effects are clearly represented in the presented analyses. We found the sites with very well drained soils (Treeline, Lower Weather) showed less uncertainty in fs and FC near the bottom of the soil profile (Figure 5), likely as a result of faster transit times and less retention. Similarly, the control of surface cover on evaporation is clearly represented in general across sites. PEL is not strictly a soil property, but also depends on vegetation where roots are active. The soil at High Sage has the same texture as Aspen, but the surface cover is dense shrub canopy 0.5-1.5 m over with 2-5 cm of litter and bare interspace and associated decreasing values of PEL at 60 cm (0.10) and 10 cm (0.08). This cover appears to effectively limit surface evaporation. Low Sage has a soil textural gradient from silt loam to clay loam and a bare soil surface and results in the strongest gradient in PEL (0.07-0.14) over the least depth of soil, 50 cm. Treeline and Lower Weather are both sandy sites with shrub cover and low bulk density surface, which supports extensive drying. Difference in PEL following vegetation removal by prescribed fire in 2007 are clear at the Aspen and High Sage sites.
There are two primary effects. First, electronic soil water sensors are based on relating measured permittivity to liquid water content, which is generally a robust relationship because the permittivity of liquid water (80) is so much greater than soil (5) or air (1). However, the permittivity of ice (3.1) is very similar to that of air, so frozen soil appears to have less water content and effectively approximates the liquid water content of the soil (Seyfried and Murdock, 1996). This effect is apparent when dramatic "drying" events during the winter appear as noise and stray from the FC attractor. Freezing effects are apparent in all four panels in Figure 2 and can be clearly seen by comparison of the soil moisture trends at 15 cm and 30 cm in Figure 9. Second, freezing effectively reduces the soil water potential. When the near-surface soil freezes, the normal downward hydraulic gradient is reversed and water moves upward toward the freezing front. This can cause either an increase or decrease in the measured soil water content, depending on the proximity of the sensor to the freezing front (Hillel, 1980). These are most pronounced near the soil surface, which commonly freezes in this environment. In this study the freezing front rarely exceeded 30 cm. There is evidence of occasional upward gradient effects apparent deeper in the soil profile. These artifacts can introduce erroneous data into the analysis. This is more problematic for the frequency and paired analyses, in which the surface sensors are likely to be similarly affected, but less important for the hysteresis analysis. We suggest that soil temperature records, which are commonly collected along with TDR or other soil moisture measurements, could be used to censor data during periods of soil frost prior to conducting the presented analyses.

CONCLUSIONS
We demonstrate a novel approach to extract important soil hydrologic parameters directly from field data. The approach takes advantage of the increasing availability of continuously monitored . Because it is based on in-situ data, results are directly related to the local climate, soils and vegetation and does not rely on assumed pedotransfer functions or proxy measurements such as . The approach builds primarily on concepts of s , FC, and PEL inherited from irrigated agriculture, but regards these quantities in terms of the amount of water retained in the soil, rather than in terms of soil water potentials. This approach is consistent with mass-based modeling schemes and provides detail for estimation of other additional transitional and extreme states fs , ws . However, the applicability of measured data for this model parameterization depends on the nature of soil water dynamics under consideration. The implicit assumptions are that drainage from s or fs to FC is relatively fast in the absence of soil water input, and below that threshold 1D unsaturated flow responds primarily to plant water uptake, which decreases progressively below ws and ceases at PEL. The endpoint of drying in any annual cycle may range from a value greater than PEL to r , depending on plant type and soil texture and depth. These assumptions are supported by observations and data from the sites in this study, which indicate a downward wetting front progression for every major rainfall or snow melt and the regular and extensive periods of soil wetting and drying. Application of the developed approach in more humid climates, or in the presence of a persistent phreatic surface is expected to improve estimation of s , but may provide little guidance on PEL or ws . Preferential flow is not specifically addressed, but the non-unity slopes in the paired analysis indicate that variable rates of wetting front advance can complicate interpretation of this approach. Similarly, vertical bypass flow in individual sensor profiles may complicate interpretation of FC and ws , but not PEL.
We found that the frequency of measured values tended to be greatest around observable "attractors" that correspond to the specific hydrologic parameters of interest. Of the four approaches analyzed, we found the hysteresis analysis approach is the most robust predictor of data attractors and frequency analysis is the simplest approach to determine the extreme values. Nevertheless, we consider the construction of a two dimensional frequency matrix as the most effective and practical single approach to parameter value identification. This approach enables analysis from a single profile of sensors and provides the greatest visual support to estimate parameter values. The data and analyses presented here clearly demonstrate how near surface often differs from the deeper soil at many sites, due to differences in macroporosity, soil freezing, hysteresis in wetting. These controls on soil profile storage may confound attempts to measure soil profile storage by remote sensing methods, and opportunities to overcome some of these obstacles will be addressed in a successive paper.

AUTHOR CONTRIBUTIONS
DC: Conceived research concept and analyses, conducted fieldwork, was primary author. MS: Contributed to research concept, conducted fieldwork and laboratory experiments, was second author. JM: Contributed to research concept, conducted fieldwork, contributing author. KH: Conceived research concept and analyses.