Slope and aspect controls on soil climate: Field documentation and implications for large‐scale simulation of critical zone processes

Soil climate, as quantified by soil temperature (TS) and water content (θ), exerts important controls on critical zone processes. It may be sensitive to variations in local slope and aspect (SA), but this attribute remains poorly quantified at the local scale and unresolved in large‐scale models. Estimation of SA effects on soil climate across multiple scales may facilitated using topographically modified, incoming clear‐sky solar radiation (SR,CS,T). We established six paired automated soil climate monitoring stations on opposing north‐facing (NF) and south‐facing (SF) slopes (4 yr) and collected spatial TS and θ data within the hectare surrounding four stations (2 yr) to measure SA effects on soil climate. Results were compared with physically based simulations and evaluated in the context of SR,CS,T. Spatial θ data were more variable than Ts, and both were consistent with values from continuous monitoring stations. On average, the SF TS was much greater (4.7 °C) and the annual summer drought longer (36 d) than on the adjacent NF aspect. Seasonal variations of TS and θ were different from each other and also different from SR,CS,T. Local conditions, including snow cover, precipitation patterns, and soil properties, largely controlled seasonal variations of TS and θ, which cannot be predicted from SR,CS,T. This indicates that realistic simulation of many critical zone processes requires high‐resolution inputs. Simulations captured first‐order SA effects and could be useful for estimating SA effects in lieu of field monitoring.

zone such as stream flow generation, plant productivity, C cycling, and mineral weathering (Riveros-Iregui & McGlynn, 2009;Seyfried et al., 2009;Stielstra et al., 2015). Improved quantification of C, water, and nutrient fluxes across the landscape is increasingly required to better understand and manage resources under changing climatic conditions and increasing resource demands. Earth system models currently simulate these processes at continental to global scales (Fan et al., 2019), implying effective quantification of the spatial and temporal distribution of θ and T S . These large-area simulations are built on relatively large model grid cells of about 20-200 km (Fan et al., 2019).
Largely due to limited availability of soil data, T S and θ have historically been estimated from more easily obtained measurements of precipitation and air temperature (T a ). These are obviously important determinants of soil climate and have been shown to be useful for regional-and continentalscale predictions of soil climate (Jenny, 1980;Minasny et al., 2013;Rasmussen, 2006). In mountainous regions, climatic variations are often linked to elevation, enabling estimation of high-elevation conditions from low-elevation data. Other factors, such as vegetative cover and topography, may exert important controls on soil climate at smaller spatial scales than can be resolved with 20-to-200-km grid cells (Bond-Lamberty et al., 2005;Ebel, 2012;Flerchinger et al., 2010;Ivanov et al., 2008). The cumulative magnitude of these smallscale effects may negatively affect the accuracy of larger-scale modeling that generally does not resolve processes at finer spatial scales (Freund & Kirchner, 2017). Where these smallscale effects are deterministic in the sense of being directly related to known, measurable parameters, it may be possible to reasonably represent them in earth system models (Fan et al., 2019).
In this paper, we focus on topography, in particular upland processes associated with the orientation of the land surface relative to the sun, which is locally controlled by land surface slope and aspect (SA). Although direct measurement of SA effects on soil climate are limited, there is considerable indirect evidence that SA controls on soil climate are both substantial and extensive. Differences in plant species, density, and vegetative productivity on slopes with contrasting aspects have been widely reported (Birkeland et al., 2003;Burnett et al., 2008;Cerda, 1998;Desta et al., 2004;Florinsky et al., 1994;Jenny, 1980;Swetnam et al., 2017;Whittaker & Niering, 1965). Similarly, differences in soil development and C content have been documented on contrasting slopes (Anderson et al., 2014;Beaudette & O'Geen, 2009;Birkeland et al., 2003;Kunkel et al., 2011;McNamara et al., 2018;Patton et al., 2019;Rech et al., 2001). These observations, which imply differences in soil climate, are often described in somewhat qualitative terms, such as mesic or xeric, which are difficult to quantify in terms of soil climate.
The fundamental driver of SA effects on soil climate is the differential incoming shortwave solar radiation (S R ) incident to slopes with contrasting aspects, which is affected by season and latitude as well as SA. The widespread availability of topographic information and simple calculation of clearsky solar radiation (S R,CS ) suggests the potential for nearuniversal estimation of SA effects on soil climate (Chorover et al., 2011). Although this approach cannot yield a local energy balance or calculate soil climate directly, it may prove

Core Ideas
• Spatial variability of soil water content was much greater than that of soil temperature. • Soil climate was very different on slopes with contrasting aspects. • Local conditions control effect of solar radiation on soil climate in complex terrain. • Slope and aspect effects can be well represented by simulation and/or point-scale monitoring.
useful if the measured SA effects on soil climate are consistent with S R, CS corrected for topography (S R,CS,T ). Within a given (temperate) latitude, equator-facing slopes receive more S R on an annual basis than pole-facing slopes. For example, the annual S R,CS,T incident on a 20˚equatorfacing slope at a latitude of 40˚is approximately twice that on an adjacent, pole-facing slope, as calculated by the methods of Tian et al. (2001). This difference between slopes with contrasting aspects varies seasonally, with the maximum difference around the winter solstice and minimum difference around the summer solstice (Tian et al., 2001;Zou et al., 2007). Based strictly on S R,CS,T , we expect that 1. On an annual basis, T S is greater and θ is lower on equatorfacing than on pole-facing slopes because the greater equator-facing S R, S,T differentially warms soils and generates greater evaporative demand than on pole-facing slopes. 2. On a seasonal basis, T S and θ differences between contrasting equator-facing and pole-facing slopes are greatest in winter, when the effects of lower solar angles result in the greatest difference in S R,CS,T , and least in summer, when S R, CS, T is nearly equal on both slopes.
The seasonality, or timing, of SA effects is important because both T S and θ control most of the relevant processes (Pelletier et al., 2018). For example, soil respiration may be limited by either very dry or cool conditions (Moyano et al., 2013;Wood et al., 2013). For S R,CS,T to be an effective basis for estimating SA effects on soil climate, the expected trends described above should be consistent with field observations. Otherwise, other factors, probably driven by specific local conditions, dominate SA effects and more site-specific simulation is required.
The limited field data available support the first expectation with regard to annual T S , with annual measured T S on equator-facing slopes being substantially greater (∼5˚C) than on pole-facing slopes (Burnett et al., 2008;Ebel, 2012;Gutiérrez-Jurado et al., 2013;Radcliffe & Lefever, 1981; Vadose Zone Journal Seyfried et al., 2016). However, the reported seasonality of T S does not always support the second expectation. In some cases, T S differences between contrasting aspects were consistent with S R,CS,T , with a winter maximum and summer minimum (Burnett et al., 2008;Gutiérrez-Jurado et al., 2013). In contrast, Radcliffe and Lefever (1981) found that T S differences between contrasting slopes were virtually constant throughout the year, and Seyfried et al. (2016) found a reversal in seasonality, with a late summer maximum and winter minimum.
Generalizations regarding θ on contrasting aspects are not clear. In some cases, it appears that spring transpiration is greater on equator-facing slopes (Ebel, 2012;Ebel, et al., 2012;Langston et al., 2015;Reid, 1973) but with little difference in measured θ over time. In contradiction with expectations, Radcliffe and Lefever (1981) found that θ on a pole-facing slope was greater than the equator-facing slope throughout the growing season but similar during the winter. Where snow is significant, θ values on opposing slopes may be determined by the timing of snowmelt (Hinkley et al., 2012;Reid, 1973) and/or spatial distribution of snow (Flerchinger et al., 2010;Kormos et al., 2015;Seyfried et al., 2009). Differences in θ between slopes may also be confounded by soil differences (Geroy et al., 2011).
One explanation for the lack of consistency among reported data is that the spatial and temporal variability of T S and θ is so great that deterministic patterns of differences between aspects are obscured by random variability. In general, the data described above were collected with little or no spatial replication for a limited time duration. The underlying assumptions are that the measurement locations are representative of the slope they occupy and that interannual variations in T S and θ are relatively unimportant. For T S , the only relevant data we are aware of does support the assumption that point data may effectively represent a slope, but no interannual data were provided . Regarding θ, there is a wealth of information on spatial variability (Grant et al., 2004;Seyfried, 1998;Western et al., 2002;Wilson et al., 2004), but this has not been evaluated in the context of SA. It is, however, clear that θ generally exhibits a high degree of spatial variability (Kutilek & Nielsen, 1994;Wilding & Drees, 1983), so that distinctions between slopes may be problematic. An alternative explanation is that local conditions may variably affect the exchange of heat and/or water between the atmosphere and soil in ways that alter or contradict expectations based on S R,CS,T . At present, there is little empirical basis for accepting either explanation.
The general objective of this research is to advance the understanding of how soil climate varies in complex terrain using a combination of replicated continuous monitoring stations, spatially extensive periodic sampling, and physically based modeling to determine the magnitude and seasonality of differences that result from contrasting slope and aspect. The investigation is focused on the following four questions: 1. Are SA effects on T S and θ subsumed by background random spatial variability, or can point data effectively represent hillslopes? 2. Is the magnitude of SA effects on T S and θ sufficient to warrant specific consideration in larger scale ecohydrologic and biogeochemical models, or can those effects be ignored as insignificant? 3. Can observed SA effects on soil climate, as quantified by seasonal T S and θ differences between contrasting slopes, be directly attributed to S R,CS,T , or is additional, local, sitespecific information required? 4. Can current, process-based modeling adequately simulate SA impacts on soil climate, or is further model development required?
We then consider the implications that the answers to those questions have for modeling important critical zone fluxes related to T S and θ.

Site description
The study site, Johnston Draw, is in the Reynolds Creek Experimental Watershed (RCEW) and Critical Zone Observatory. The 1.8-km 2 subbasin ranges from 1,490 to 1,850 m in elevation, which spans the current rain-snow transition zone during most winter storms. The climate is semiarid, with a mean annual precipitation of 550 mm and mean annual temperature of 7.4˚C at Station 125 ( Figure 1). The site experiences a pronounced dry period in July and August with little to no precipitation. Streamflow is intermittent, ceasing during the annual summer dry period and resuming during the winter. The lower (eastern) half of Johnston Draw flows nearly due east. Our study focuses on this lower part of the watershed, where the polar, or north-facing (NF) slope, has an average gradient of 21˚and aspect of 16˚, and the equatorial or south-facing (SF) slope has an average gradient of 17˚and aspect of 146˚. These steep slopes and nearly north-south facing aspects were chosen to represent a strong expression of SA effects on soil climate. Soils of the study area are underlain by granite and are generally classified as Haploxerolls on the NF slope and Xeropsamments and Xerorthents on the SF slope. Typical soil profile characteristics for NF and SF soils are reported in Patton et al. (2019). The lower soil boundary separating relatively hard granite and "loose" soil (or mobile regolith) is typically more abrupt on the SF slope. Soils on both slopes are sandy, F I G U R E 1 Johnston Draw (JD) is located within the heavily instrumented Reynolds Creek Experimental Watershed and Critical Zone Observatory. Stations jdt2, jdt3, jdt4, jdt2b, jdt3b, and jdt4b (referred to a 2, 3, 4, 2b, 3b, and 4b, respectively) were instrumented to measure depth profiles of soil temperature (T S ) and water content (θ). Spatial data were collected in 100-m-long transects centered on Stations 2, 4, 2b, and 4b. Weather data were measured at Site 125, and streamflow at 125b usually loamy sands or sandy loams, with a large component of small coarse fragments.
Mountain big sagebrush [Artemisia tridentata Nutt. ssp. vaseyana (Rydb.) Beetle] is a large component of the vegetative cover on both slopes. On the SF slope, bitterbrush [Purshia tridentata (Pursh) DC.] with mixed western Juniper (Juniperus occidentalis Hook.) are also common. A substantial part of the ground cover is composed of annual grasses, forbs, and bare ground. On the NF slope, snowberry (Symphoricarpos oreophilus A. Gray) is also common, and there is a dense ground cover of perennial grasses and forbs. Alder trees (Acer glabrum Torr.) are found in depressions on the NF slope.
Peak vegetative cover on the NF slope is denser and more vigorous than on the SF slope and varies strongly by season. The 16-yr average normalized difference vegetation index (NDVI), as measured with moderate resolution imaging spectroradiometer (MODIS) determined from a representative, 250-m pixel from each slope (Figure 2), illustrates this seasonality (note that winter data are excluded due to complications with snow cover). Both slopes follow the annual cycle of increasing NDVI during the spring and decreasing NDVI during summer. In early spring, NDVI on both slopes is similar but increases more rapidly on the NF slope with significantly greater NDVI values throughout late spring and into fall. By F I G U R E 2 Average normalized difference vegetation index (NDVI) on the north-facing (NF) and south-facing (SF) slopes calculated from moderate resolution imaging spectroradiometer (MODIS) data collected over a 16-yr period. Difference (Diff) is calculated as the 16-yr average NF NDVI -SF NDVI, with significant (α = .10) differences (Sig) indicated by an asterisk. The 90% confidence interval around both curves is shaded. Periods of consistent snow accumulation are omitted late fall, the vegetation on both slopes is either dead or dormant and NDVI returns to values similar to early spring.

Continuous station measurements
The soil climate data network at Johnston Draw was built upon a previously established set of meteorological stations designed to monitor snow accumulation and melt across the rain-snow transition zone and established at 50-m elevation intervals from near the bottom to the top of the watershed (Marks et al., 2013). The soil climate data network consists of six stations (Figure 1). Stations 2, 3, and 4 are on the NF slope at elevations of 1,600, 1,650 and 1,700 m asl. We added three stations, 2b, 3b, and 4b, at the same elevations on the opposing SF slope (i.e., 2b corresponds to 2, etc.). Soil depth profiles of soil water and temperature sensors were installed at each station. On the NF slope (Stations 2, 3, and 4), the depth to relatively hard bedrock was about 100 cm and the sensors were installed at depths of 5, 20, 50, 75, and 100 cm. On the SF slope, relatively hard bedrock was encountered at shallower depths, and instruments were installed at depths of 5, 20, 35, and 50 cm at all three (2b, 3b, and 4b) stations with an additional sensor at 2b at 75 cm. Soil temperature and θ, were measured hourly using Hydra Probe II sensors (Stevens, 2007). Soil water content was calculated from the real dielectric permittivity measured at 50 MHz using the "general" calibration equation recommended by Seyfried et al. (2005). Soil temperature was measured using thermistors within each Hydra Probe.

Spatial soil data collection
Due to the expense of instrumenting numerous sites for continuous monitoring, we made periodic, spatially extensive measurements that provided temporal "snapshots" of T S and θ variability. We measured T S at a depth of 30 cm (T S,30 ) with a thermistor (Omega Engineering HH806AU, THSS-18G-RSC-12) and θ from 0-30 cm (θ 30 ) with a Soil Moisture Mini Trase time-domain reflectometer (Jones et al., 2002) using 30-cm waveguides inserted vertically. Field measurements were made surrounding Stations 2, 4, 2b, and 4b in a 10-m grid of 22 measurements centered on the instrument stations and oriented to cardinal directions. The result was field measurements spanning 50 m in each cardinal direction from the corresponding station ( Figure 1). The 30-cm depth was used because T S at that depth exhibits very little diurnal fluctuation allowing for comparisons throughout the day. Similarly with θ, barring a major rainfall event, the depth integrated 0-to-30-cm water content changes only slightly in the course of several hours. Field measurements were taken approximately monthly from June to August 2011, and from April to September 2012, for a total of eight samplings.

Hydrometeorological data
Meteorological and snow depth data from Johnston Draw during this time are described by and available through Godsey et al. (2018). We used hourly snow depth measured with sonic depth sensors (Judd Communications) at all six stations, and T a measured at Stations 2, 4, 2b, and 4b to analyze SA effects on those variables. For simulation, T a wind speed and relative humidity from Stations 4 and 4b were used with precipitation and incoming solar radiation from Station 125 ( Figure 1). Stream discharge was measured at a weir installed in 2003. Johnston Draw typically starts flowing in early winter and ceases around mid-July during the very dry summer. Stage height was converted to stream discharge using a rating curve and frequent field measurements to ensure high-quality flow records (Pierson et al., 2001). Average stream discharge over the period of record is approximately 0.007 m 3 s −1 (0.33 mm d −1 ) with the largest discharge of 1.63 m 3 s −1 (78.2 mm d −1 ) on 14 Feb. 2014 during a rain-on-snow and frozen soil event.

Model
We used the Simultaneous Heat and Water (SHAW) model (Flerchinger, 1991;Flerchinger & Saxton, 1989;Flerchinger et al., 1996) to simulate T S and θ. The SHAW model is a onedimensional model that includes a vegetation canopy, snow cover (if present), plant residue, and the soil profile. The surface energy balance is determined from weather inputs and surface conditions and then the interrelated heat, liquid water, and vapor fluxes are calculated within and from the soil to the atmosphere and lower model boundary. These processes are fully coupled in the model, which also calculates soil freezing processes and snow accumulation and melt. The model has been tested and applied extensively over a range of vegetation types in semiarid environments including different parts of the RCEW (Chauvin et al., 2011;Flerchinger et al., 2016). In SHAW, incoming solar radiation incident to a sloping surface (S R,T ) is computed from the measured total incoming solar radiation (S t ) on a horizontal plane, which consists of direct (or beam, S b ), and diffuse (S d ) components. Total incoming solar radiation is separated into the two components by the following equation developed by Bristow et al. (1985): where τ d is the atmospheric diffuse transmission coefficient where β is the angle which the sun's rays make with the sloping surface and ϕ s is computed based on the latitude of the site, the time of year, and the hour of the day. Total radiation incident on the sloping surface is the sum of S s and S d , assuming that there is negligible adjacent terrain shading. Model simulations for the NF and SF slopes were initialized with measured soil temperature and moisture profiles on 1 Oct. 2010 and run through 31 Dec. 2014. Soil texture and moisture release curves were parameterized based on field measurements. Saturated hydraulic conductivity was estimated from pedotransfer functions (Saxton & Rawls, 2006).

Spatial variability
The standard deviation of T S,30 around each station was similar and ranged from 1 to 2˚C (Figure 3). There is a positive correlation between the standard deviation and average T S,30 for each site that explains much of this range. The average T S,30 values at the two SF stations were within about 1˚C of each other on every measurement date, and differences were not significant (α = .05) on any of the measurement dates.
On the NF slope, T S,30 at Station JD4 was always cooler than at Station JD2, about 1.5˚C on average, and significantly (α = .05) cooler on four of the eight measurement dates. Comparing across opposing slopes, Station JD2b was significantly warmer than Station JD2 on seven of the eight measurement dates, and Station JD4b was significantly warmer than Station JD4 on all measurement dates. The average of the two NF stations was lower than the average of SF stations on every measurement date with an overall average difference of 4.1˚C. These data also demonstrate that the continuous station data are representative of the larger adjacent hillslope. For each spatial measurement date, the continuous data are within one standard deviation of the spatial data. In fact, values are nearly identical, with a high correlation coefficient (.98), a slope near one (0.98), and a small offset (0.80˚C) when the two data types are regressed. Though no spatial data were collected over winter, the range of measured T S,30 range of more than 15˚C encompasses much of the total annual T S,30 range, and it is reasonable to expect the spatial data to follow the continuous data for the entire year.

Continuous measurements
We illustrate the continuous T S data from each of the six stations in Figure 4 using measured values at a depth of 50 cm (T S,50 ). This depth was chosen because it is common to all stations and does not exhibit the pronounced diurnal fluctuations observed at more shallow depths that obscure seasonal trends, making visual comparison difficult. Comparisons among stations at this depth represent other measurement depths because seasonal patterns, with a slight damping and time lag, are similar at all depths and because long-term statistics, such as mean annual soil temperature, are virtually identical at all depths within a given profile. Clear T S,50 differences between the two aspects are obvious in Figure 4, with the three SF stations warmer than the three NF stations throughout the year. Consistent with the spatial data, T S,50 for the three SF stations was nearly identical. Annual maxima were about 24˚C in mid-August and minima were about 2˚C in January. Soils on the SF slope did not freeze at 50 cm during the study, although there was some limited freezing closer to the soil surface. On the NF slope, T S,50 followed the same seasonal pattern but was always cooler, with maximum values between 15 and 20˚C and minima between 1 and −5˚C. Unlike the SF stations, there were consistent differences among the three NF stations. During the summer months, T S,50 values were such that JD2 > JD3 > JD4. During autumn, those differences abated and there was a slight tendency for a reversal of order. At the onset of snowmelt in spring, T S,50 was nearly equal for the three stations at about 0˚C.

Vadose Zone Journal
F I G U R E 4 Continuous station T S,50 (soil temperature at a depth of 50 cm) data over the 4-yr study period. South-facing (SF) data are represented with warmer colors, north-facing (NF) stations with cooler colors. Snow depths are the average measured at three sites on the NF slope (dark gray) and of the three sites on the SF slope (light gray) F I G U R E 5 Time series illustrating the seasonality of slope and aspect (SA) effects on soil temperature at 50 cm (T S, 50 ) and topographically modified, incoming clear-sky solar radiation (S R,CS,T ). ΔT S,50 is the difference between the south-facing (SF) T S,50 and the north-facing (NF) T S,50 (SF T S,50 − NF T S,50 ) averaged for each day of the year over the 4-yr study period, ΔT S,sim is the corresponding simulated value, and ΔS R,CS,T is SF S R,CS,T -NF S R,CS,T , as calculated from Tian et al. (2001), for each day of the year. Values of ΔT S,50 and ΔT S,sim were "smoothed" with a 5-d moving average Snow depth differed among slopes (Figure 4). The NF slope developed a seasonal snowpack each year, whereas the snow cover at the SF stations was ephemeral as snow that fell on the SF slope quickly melted. On average, the duration of snow cover was 119 d on the NF slope compared with 37 d on the SF slope. A rapid increase in T S,50 followed the snow ablation date each year on the NF slope. Those dates ranged from 13 Mar. 2014 to 12 Apr. 2011.
To quantify SA effects on T S , we calculated the difference between the average of the three SF stations, SF T S,50 , and that of the three NF stations, NF T S,50 , such that ΔT S,50 = SF T S,50 -NF T S,50 . The mean annual ΔT S,50 was almost 5˚C (4.7˚C), with a mean annual SF T S,50 of 11.9˚C and NF T S,50 of 7.2˚C. This result is consistent with initial expectations based on S R,CS,T and similar to the results of other studies referred to previously. These annual differ-ences were measured while the difference in mean annual T a between JD4 and JD4b was negligible (0.3˚C), supporting the importance of S R as a driver of T S differences across slopes.
Unlike the average annual ΔT S , seasonal values contradict some reported results and are quite different from the patterns expected from S R,CS,T . The four year average daily minimum ΔT S,50 was in winter at just over 2˚C, when the analogous difference in S R,CS,T , ΔS R,CS,T (SF S R,CS,T -NF S R,CS,T ) is greatest ( Figure 5). The measured maximum ΔT S,50 of almost 8˚C occurred in fall, somewhat prior to the maximum ΔS R,CS,T.
Peak spring values were about 2˚C less than the corresponding fall values, whereas ΔS R,CS,T is the same at those times. Consistent with ΔS R,CS,T , there was a local minimum at the summer solstice, but it was about 1.5˚C greater than the winter minimum. F I G U R E 6 Slope-average measured soil temperature at a depth of 50 cm (T S,50 , three stations on each slope) compared with the simulated value (T S,sim )

Simulation
To evaluate the simulation of T S , we compare the aspectaverage (three station) 50-cm measured soil temperature, NF T S,50 , and SF T S,50 , with the corresponding simulated values of NF T S,sim and SF T S,sim ( Figure 6). This is reasonable because the model parameterization of each individual site is nearly identical and thus the simulation effectively weights each site on a given slope equally. Visually, the annual cycle of NF T S,sim is in close agreement with NF T S,50 . Major freezing events and snowmelt timing, as indicated by the date of rapid T S,50 increase in spring, were also in close agreement with measured values. This is supported by the regression of NF T S,50 and NF T S,sim , which had a high r 2 (.98), slope near one (1.10), and an intercept near zero (0.63˚C). Differences are roughly symmetric around the measured values such that NF T S,sim tends to be greater in summer and lower in winter than NF T S,50 . The SF T S,50 was almost always slightly warmer (about 1˚C) than SF T S,sim , but the overall agreement was also excellent, with a high r 2 (.99), slope near one (1.01), and a modest intercept of 1.19˚C, which was the primary difference between SF T S,50 and SF T S,sim . Annual average values of T S,sim were similar those measured, with SF T S,sim = 10.6˚C and NF T S,sim = 7.4˚C. Consistent with these overall observations, ΔT S,sim (SF T S,sim -NF T S,sim ) captures the basic pattern of minimum values in winter, a local maximum in early spring, a local minimum around the summer solstice, and a fall peak ( Figure 5). However, the model did not capture the full extent of this latter trend. This can also be seen in Figure 6, where, each year, the largest discrepancy between simulated and measured values occurred in the late summer and early fall, with a slight (≈1˚C) underestimation of SF T S,50 and similar slight overestimation (≈3˚C) of NF T S,50 , resulting in a reduced ΔT S,sim relative to ΔT S,50 .

Spatial variability
In general, the spatial variability of θ 30 around the two SF stations, as indicated by the standard deviation, was about twice that around the NF stations (0.29 m 3 m −3 for SF stations versus 0.15 m 3 m −3 for NF stations). The standard deviation of θ 30 increased with the average θ 30 at all stations, roughly doubling over the range of measured values. Unlike T S , there is no clear θ difference between the slopes (Figure 7). On some dates, the two SF stations (2b and 4b) have both the highest and lowest average θ 30 among the four stations. Rather, the primary difference among stations is that 2b always had a higher average θ 30 than the other three stations, which were similar. Statistically, θ 30 at station 2b was significantly (α = .05) greater than the other three stations on seven of the eight measurement dates that all four were collected, whereas the three other stations were not significantly different from each other. On average, spatial variability around a given station, as indicated by the coefficient of variation, was about three times greater for θ 30 (27%) than for T S,30 (9%) (note that we use coefficient of variation here because it is the usual metric in soil variability studies and because the spatial T S,30 values were always greater than 0˚C). As with T S , the continuous station θ data were highly correlated with the spatial data, although correlation was weaker (r 2 = .74) and the spatial data were somewhat lower than the continuous data, with a regression slope of 0.76. Similarly, the number of dates the continuous data fell within one standard deviation of the average was much lower than for T S . The greater variability of θ relative to T S is evident comparing Figures 7 and 3. Although these data also support the use of station data to represent a much larger slope area, they do not indicate a difference between slopes.

Vadose Zone Journal
F I G U R E 7 Spatial average and standard deviation of spatial water content (θ) measured at 0-to-30-cm depth and corresponding continuous station data calculated as the weighted average of θ measured at 5-, 20-, and 35-cm depths, represented as solid lines. Spatial data agree closely with the continuous and do not show distinct cross-slope differences

Continuous measurements
Comparisons of soil water across sites is less straightforward than for T S for two reasons. The first is soil freezing. Electronic sensors, like the Hydra Probe or TDR, measure, to a close approximation, the liquid water content of the soil when ice is present (Seyfried & Murdock, 1996), so that direct comparison among frozen and unfrozen soils is not possible. Soil freezing also has the effect of dramatically altering the soil water potential, creating gradients that may result in redistribution (generally upward) of soil water within the soil profile (Flerchinger, 1991;Flerchinger & Saxton, 1989;Miller, 1980). Second, for many ecohydrological and biogeochemical processes, the measured θ at a given depth is not necessarily the value of interest and may, in fact, be misleading because differences in soil texture may obscure important soil water potential differences. In terms of soil water measurements, parameters of field capacity (θ FC ) and the plant extraction limit, θ PEL , are useful for hydrologic purposes (Finzel et al., 2016;Seyfried & Wilcox, 2006) and can be identified from field data (Chandler et al., 2017). Similarly, levels of microbial activity are also better described by θ values normalized to the same parameters of θ PEL , θ FC in addition to the saturated value of θ SAT (Moyano, 2013;Paul et al., 2003), and not the actual θ value. In addition, plants respond to all water within the rooting zone, as opposed to water from a particular depth in the profile.
For these reasons we evaluate soil water relations in terms of the effective soil water storage, S e , expressed in centimeters of water stored and calculated as where θ PAW is the plant available water content, θ PAW = θ − θ PEL , for a given depth increment i, multiplied by the thickness of that depth increment (THK i ). In this case, soil layers are defined by the vertical space between soil water sensors in the profile. Note that, when all soil layers have θ of θ PEL , S e = 0, or S e , PEL , and that S e when all layers are equal to θ FC , S e is equal to the effective soil storage capacity, S e,FC . Soil water trends, expressed in terms of effective storage to a depth of 55 cm (S e,55 ) at all six stations (Figure 8), generally follow temporal trends similar to those described by Grant et al. (2004) and Rawls et al. (1973) at Reynolds Creek and more recently noted by others in a wet winter, dry summer environment (Bales et al., 2011;McNamara et al., 2005;Oroza et al., 2018;Seyfried et al., 2011). During the warm dry summer, S e,55 values are low, less than 2 cm of water at all sites, as evaporative demand greatly exceeds precipitation. In winter, S e , 55 remains near S e , FC , as precipitation exceeds evapotranspiration (ET) and excess soil water drains to greater depths. These two stable periods are linked by often abrupt transitions we refer to as wet-up in autumn transitioning from summer to winter conditions, and dry-down in spring when S e,55 returned to values near 0 cm. Consistent with the spatial data, the three NF stations track each other closely throughout the period while one SF station, 2b, clearly retains more water than the other SF stations.
There are three key differences between the NF and SF slopes which we illustrate in terms of slope-average S e,55 for the NF (NF S e,55 ) and SF (SF S e,55 ) slopes and slope-average NF S e,105 in Figure 9. The first is that soil freezing had a greater impact on NF S e,55 than on SF S e,55 . Periods when T S measured at one or more depths on the slope were less than 0˚C are highlighted as shaded regions in Figure 9. This difference was especially clear during the winters of 2011-2012 and 2013-2014. In late 2011, soils on both slopes froze, but the SF slope thawed in January while the NF soil remained F I G U R E 8 Continuous station effective storage to a depth of 55 cm (S e,55 ) over the 4-yr study period. South-facing (SF) data are represented with warmer colors, and the north-facing (NF) stations are represented with cooler colors F I G U R E 9 Slope-average effective storage to a depth of 55 cm (S e,55 ) for the north-facing (NF) and south-facing (SF) slopes and slope-average effective storage to a depth of 105 cm (S e,105 ) for the NF slope. No corresponding S e,105 is calculated for the SF slope due to the relatively shallow soils on that facet. Periods when soil freezing was detected at any of the sensors (all depths) are shaded. The horizontal line at 3.5 cm demarks the condition of "dry" soil frozen and appeared to be dry until mid-March of 2012 ( Figure 4) (note that values of S e less than 0 are possible when the soil is frozen because soil freezing can result in water potentials lower than those corresponding to θ PEL ). At that time, NF S e , 55 quickly rose to the storage capacity of about 8 cm and was similar to SF S e,55 . Also, from early December 2013 to early March 2014, frozen soil at one or more of the NF station caused a drop in NF S e , 55 . During the same period, none of the SF stations froze and SF S e , 55 rose in response to winter rain. When the NF soil thawed in early March, NF S e , 55 values quickly rose to about 8 cm as in other years.
The second key difference is that the timing of spring drydown and autumn wet-up is different on the NF and SF slopes. The dry-down at the SF stations preceded that at the NF stations each year. Using 3.5 cm of S e , 55 as a somewhat arbitrary metric for "dry" soil, the SF stations, on average, became dry 7 d earlier than the NF stations in 2011, 28 d earlier in 2012, 26 d earlier in 2013, and 12 d earlier in 2014. This implies that, given equal soil depths, vegetation on the SF slope will experience about 2 wk longer summer drought than on the NF slope. The wet-up is much less regular due to a combination of frost effects in November and variable precipitation in October. However, in October of 2013, SF S e , 55 was clearly less than NF S e , 55 (Figure 9). The third key difference is that soils on the NF slope are deeper, by about 50 cm, than those on the SF slope, with a correspondingly greater storage capacity. Note that NF S e , 55 and NF S e,105 are similar and near 0 cm each summer as the vegetation effectively extracts soil water throughout the profile to near θ PEL . Values diverge each autumn when the storage capacity of about 8 cm was exceeded so that NF S e , 55 drained to greater depths and NF S e,105 continued to increase to about 14 cm. As a result, NF S e,105 greatly exceeded NF S e , 55 on both slopes in spring, and the NF dry-down date was delayed resulting in slope-average dry-down dates 24, 40, 43, and 28 d after the SF S e , 55 for 2011, 2012, 2013, and 2014, F I G U R E 1 0 Slope-average effective storage to a depth of 55 cm (S e,55 ), three stations on each slope, compared with the simulated value (S e,sim ) for both slopes. The south-facing comparison is in warm colors, and the north-facing comparison is in cooler colors respectively, for an average of 36 d. Thus, accounting for greater soil depth, the summer dry period is about a month longer on the SF slope.

Simulation
The simulated SF S e , 55 (SF S e,sim ) was overall in close agreement with the measured values both in terms of the timing of the annual wet-up and dry-down and the winter-spring field capacity-like values around 8 cm ( Figure 10). The most obvious differences were 12 Jan. 2011 and 4 Feb. 2013, when SF S e,sim decreased sharply while measured values remained relatively constant. In both cases, the model simulated freezing to a depth of 20 cm when, in fact, little or no freezing was observed. T S differences of less than 1˚C were responsible, demonstrating that very small simulation errors near 0˚C can have a major impact on simulated liquid water contents. Regression statistics of r 2 = .89, slope = 1.06, and offset = 0.12 cm support the general appearance of excellent agreement between measured and simulated values on the SF slope. Comparison of measured and simulated soil water storage on the NF slope is more problematic (Figure 10). As with the SF results, the timing of the dry-down is very closely simulated, but the NF S e,sim tended to be 1-2 cm greater in spring and more variable over short time intervals. In addition, NF S e,sim was much lower than the measured value during the wetup in 2013 and 2014. The main differences are again associated with soil freezing. The model tended to underestimate the liquid water content of frozen soil compared with the measured values. In the winters of 2011, 2013, and 2014, both NF S e,sim and NF S e,55 decreased dramatically, but in all cases NF S e,sim dropped to substantially lower values than NF S e,55 . This led to an exaggerated response "spike" immediately after thaw as relatively large quantities of stored (frozen) soil water were rapidly released. This again illustrates the high sensitivity of θ (and therefore S e ) simulations when T S is near 0˚C. These problems are reflected in the relatively low regression statistics comparing the measured and simulated S e,55 on the NF slope with r 2 = .58, slope less than 1 (0.832), and intercept of 0.4 cm.

Spatial variability
The within-slope spatial variability data lend insight into the natural variability of T S relative to θ due to stochastic or other relatively subtle factors not easily considered. The coefficient of variation of θ was generally between 0.15 and 0.4, which is larger than that reported by Western and Grayson (2000), but not inconsistent with typical values (>0.05) suggested by Wilding (1985) and generally in the range reported by Seyfried (1998) and Williams et al. (2009) of between 0.05 and 0.35. The corresponding values for T S were always lower than those for θ, by a factor of about three, and are also evident comparing Figures 3 and 7. The T S variability reported here is similar to that observed by Seyfried et al. (2016) in a nearby subwatershed in the RCEW. This observation of greater θ than T S variability can probably be generalized. The soil hydraulic parameters (e.g., hydraulic conductivity and water holding capacity) that control θ are more sensitive to variations in soil properties than the soil thermal parameters that control T S (e.g., soil heat capacity and thermal conductivity). This was exemplified in our spatial dataset by one point near Site 4b that happened to intersect an intermittent clay layer. As expected, θ at that point was always greater than the surrounding, coarse-textured points while T S at that point was not distinguishable from neighboring points. In addition, thermal fluxes tend to be more isotropic than hydraulic fluxes which have a strong vertical component owing to gravity.
The greater variability of θ indicates that it is generally more difficult to detect SA effects on θ than on T S . This partially explains the lack of difference observed between slopes in this study and probably contributes to the mixed results reported in other studies. We discuss other critical factors in the section below. In either case (T S or θ), data from the continuous stations fit well within the spatial data representing the overall slope, supporting the use of point station data as representative of a given slope-aspect combination.

Soil temperature
The mean annual ΔT S is in the direction predicted from S R,CS,T (greater on the SF slope), substantial (about 5˚C), and generally consistent with other reported research (Burnett et al., 2008;Ebel, 2012;Gutiérrez-Jurado et al., 2013;Radcliffe & Lefever, 1981;Seyfried et al., 2016). Within the context of the larger landscape, the 4.7˚C cross-slope difference is about the same as the 900 m elevation difference within the RCEW (see Seyfried et al., 2016). If we use T a as a surrogate for T S , this is roughly equivalent to a 5˚latitude shift in the central United States or to the 6˚latitude shift estimated by Radcliffe and Lefever (1981) in New Zealand. Both are much greater than the typical cell size for continental-scale models Fan et al. (2019). The fact that these differences were observed while T a differences across slopes were negligible is notable because models often assume that T S varies directly with T a , which is clearly not the case in Johnston Draw. This result is probably generalizable to areas where T a is measured above the plant canopy, thus allowing for free circulation of air across slopes . The seasonality of ΔT S was very different from ΔS R,CS,T . The primary factor responsible for this is the differential snow cover on the opposing slopes, which affected the magnitude and interannual variability of ΔT S,50 . Snow cover affects T S because it is an excellent insulator, reducing loss of heat from the soil and thus effectively "warming" the soil during winter. This reduced soil cooling under snow cover has been widely reported (Hardy et al., 2001), but not in the context of SA.
During periods of significant snow cover, near surface diurnal T S fluctuations are strongly damped and T S at all depths trended slowly towards 0˚C. In most winters, a snowpack developed on the NF slope prior to soil freezing at 50 cm and T S,50 gradually decreased towards 0˚C during winter while SF T S,50 and hence ΔT S,50 also decreased. In the winter of 2011-2012, a snow pack failed to develop until relatively late in the season and T S,50 below −5˚C were recorded, resulting in much greater ΔT S,50 values for that year. After the development of a snow pack T S,50 warmed toward 0˚C until the snow ablated.
A second impact of snow cover is that it provides a 0˚C boundary condition at the soil-snow interface during warming conditions in spring, effectively maintaining T S near 0˚C until snow ablation. The effect is that, while ΔS R,CS,T decreased during March, ΔT S,50 increased because the snowfree SF soils warmed while NF T S,50 remained near 0˚C. Once the snow completely ablated, NF T S,50 also increased, ΔT S,50 stabilized and began to follow the trends of ΔS R,CS,T in that it declined until the summer solstice and increased through the fall equinox.
The observed ΔT S,50 in late summer of nearly 8˚C may have been enhanced by the denser plant cover on the NF slope during summer ( Figure 2). As the annual grass or forb cover on the SF slope senesced, bare soil was exposed to direct solar radiation while the NF slope was consistently covered by a perennial sod. Although the effect of vegetative cover in general on T S is well documented (Bond-Lamberty et al., 2005;Flerchinger & Pierson, 1991), little has been reported in the context of SA. Research in central New Mexico (USA) has shown that vegetative cover may have a large impact on annual T S dynamics (Gutiérrez-Jurado et al., 2013;Zou et al., 2007) and serve to accentuate SA effects.
The close agreement between T S,50 and T S,sim in winter and spring indicates that the impacts of S R,CS,T , snow cover and soil freezing are well represented in the SHAW model and supports the explanation of ΔT S dynamics presented above. After 1 July, the simulations continue to accurately portray the dynamics of T S,50 with modest differences in maximum temperatures. We note that these differences are consistent with differing vegetation effects described above that may not have been accounted for in the model. These relatively subtle and dynamic effects of vegetative cover may warrant further research.
In general, we expect that the effects of SA on T S will follow S R,CS,T most closely in environments with sparse vegetative cover and little or no snow cover. Where there is significant snow cover, those effects will likely be altered and even reversed. Vegetative cover, especially where it changes during the year, may be an important factor controlling ΔT S in many environments.

Soil water
The expectation of drier soils on the SF slope was not supported by the spatial θ data. This is partly due to the relatively high spatial variability of θ but more fundamentally due to the interplay between soil properties, climate, and vegetation.
The primary mechanism creating drier soils on SF slopes is greater ET from those slopes due to greater S R . For this differential ET to be evident from measured S e , the following three conditions must be met: 1. S e < S e,FC on the SF slope (otherwise, drainage dynamics dominate). 2. S e > S e,PEL on both slopes (otherwise, ET has little effect). 3. ET from SF slope > ET from the NF slope. These conditions are met at different times depending on the local climate, soil properties, and vegetation in addition to ΔS R,T .
At Johnston Draw all three conditions consistently applied only during the spring dry-down. In winter, the first condition was not met despite maximum ΔS R,CS,T . In summer, the second condition was not met. The autumn wet-up was complicated by precipitation patterns, soil freezing, and vegetation phenology. In spring, ET from the SF slope exceeded that on the NF slope as expected from ΔS R,CS,T but did not result in drier soil until spring rain was insufficient to maintain S e near S e,FC . That occurred on different dates each year depending on the amount of spring precipitation. In the relatively wet spring of 2011, when the April to June total precipitation was 217 mm, high S e was maintained on both slopes and the drydown date on the SF slope was 29 June, when ΔS R,CS,T was minimal, followed by the NF slope only 7 d later. In contrast, the following 3 yr were much drier (and closer to average), with precipitation inputs of 123, 110, and 91 mm for the same period. Dry-down dates for the SF slope were about six weeks earlier, in mid-May, about 22 d before from the NF slope drydown date.
The soil water simulations generally captured these seasonal trends on both slopes including the dry-down window ( Figure 10). An exception is the apparent model overestimation of ET on the NF slope in the fall of 2013 and 2014, which is probably due to differential vegetation phenology. That is, the vegetation on the NF slope may not revive quickly following the summer drought. These issues require further investigation. It is clear, however, that the simulations advance our understanding of soil water dynamics on opposing slopes well beyond simple conceptual drier vs. wetter or mesic versus xeric, or greater effects in winter than summer and that highresolution simulation is good tool for estimating SA effects.

Soil depth
Soil depth, or more precisely, the soil profile water-holding capacity, may be an important factor determining the magnitude of SA effects on soil water. Given two opposing slopes with equivalent soil hydraulic properties, relatively shallow, low-S e,FC soils will be more difficult to distinguish because there will be less time that S e is between S e,PEL and S e,FC . It has often been observed that soils on NF slopes tend to be deeper and store more soil water than on SF slopes (Jenny, 1980;Geroy et al., 2011;Kutiel & Lavee, 1999), as we found at Johnston Draw. In this case, greater storage enhances the effects of SA, resulting in more transpiration and shorter summer droughts on NF slopes. On the other hand, Rech et al. (2001) found that soils on SF slopes were deeper and more highly weathered than those on NF slopes, which would tend to reduce SA effects. We should also note that soil water storage does not necessarily define the total amount of water available to the native vegetation. Recent research confirms Arkley's (1981) finding that that vegetation commonly exploits water stored well below the traditional soil boundary in the subjacent saprolite and/or weathered bedrock that is difficult to measure (Bales et al., 2011;Klos et al., 2018). This certainly is a factor at Johnston Draw, but we suspect that it is small relative to the examples from California cited above. First, based on our excavations, the rock at Johnston Draw appears to be of sufficient integrity to absorb little water and prevent root penetration except in fractures, so that the water holding capacity below the soil boundary is small. Second, use of deep water implies considerable excess of water entering the soil surface, which is clearly the case in Californian forests, but much less the case in the 500 mm precipitation regime of Johnston Draw. The cessation of streamflow each August supports this. In addition, recent ET data collected using eddy covariance has shown that transpiration rates are very low in late summer at other locations in the RCEW (Flerchinger et al., 2020).

CONCLUSIONS
Soil climate controls a number of Critical Zone processes. Although there is considerable evidence of slope-aspect impacts on soil climate, there has been relatively little quantification of the two primary parameters, T S and θ, in this context. The effects of slope-aspect on soil-related processes are expressed at spatial scales that are not resolved by large-scale, regional, or continental models but may be estimated from knowledge of incoming, clear sky solar radiation corrected for topography, S R,CS,T , without detailed, explicit simulation. We measured slope-aspect effects on soil climate and evaluated them in the context of S R,CS,T . Using 22 spatially distributed measurements in the hectare surrounding four soil climate monitoring stations on opposing slopes north-facing, or pole-facing and south-facing, or equator-facing slopes, T S and θ were measured on eight dates. We found that continuous station data were highly correlated with, and of similar magnitude to, the spatially distributed data for both T S and θ, indicating that station data are useful for characterizing larger landscape conditions. These data also showed that 1. θ was more variable that T S , with a coefficient of variation three times greater than T S , probably due to greater sensitivity of θ to other soil properties such as texture. 2. There were clear T S differences between aspects, with the SF slope always warmer than the NF slope (average difference of 4.1˚C), while no trends with aspect were apparent from the spatial θ data. 3. T S differences within sites on the same aspect were always smaller than those between aspects, which was not true for θ.
Using continuous station data from six paired stations, we found that, on an annual basis, slope-aspect effects were consistent with the general expectation that pole-facing slopes be wetter and cooler than equator-facing slopes and that the magnitude of those differences probably effects the accuracy of large-scale simulations. Regarding T S , the 4.7˚C difference between two slopes separated by 300-500 m was roughly equivalent to a north-south latitudinal shift of about 500 km (multiple simulation cells), or a 900-m elevation difference. We analyzed soil water data in terms of the effective water storage, S e , as a better descriptor of soil climate than θ at a given depth. The annual summer drought was about 36 d longer, on average, on the SF than the NF slope.
We also used the station data to analyze the seasonality of slope-aspect effects on T S and S e with the expectation that they would coincide with S R,CS,T . Seasonality is critical because most processes of interest are simultaneously controlled by both T S and S e . We found that the seasonal variation of T S was different from S e and that both were quite different from S R,CS,T . For T S , differences between opposing slopes were smallest in winter, when the S R,CS,T effect is maximal, and greatest in early fall, when it is intermediate between a summer minimum and winter maximum. Soil water differences were consistently apparent only during relatively short time periods after the spring equinox.
Accurate simulation of soil climate and associated processes in complex terrain therefore requires high-resolution, local inputs because neither T S nor θ varies directly with S R,CS,T . Both are controlled by local conditions. For T S , snow cover has an overriding effect, causing a homogenization of values across slopes in winter. For S e , factors that control the local water balance, such as soil water holding capacity, timing and amount of soil water input, and vegetation phenology, largely determine slope-aspect effects.
These findings highlight the importance of smaller scale, process-based simulations, for assessing the impacts of slopeaspect. Using the SHAW model, we found that model simulations generally agreed with continuously measured data on both aspects. For T S , close agreement visually is supported by regression statistics for both slopes and the simulated seasonal pattern of slope differences. Similarly, S e,55 was generally well simulated with the complication of soil freezing, which resulted in extreme sensitivity of liquid water contents near the freezing point and the model-simulated liquid water contents lower than we estimate from the measured values. For nonfreezing periods, the model was in close agreement with measured values and the seasonal patterns were preserved. Importantly, the spring dry-down was reasonably simulated on both slopes. We suspect that some of the differences between simulated and measured T S and θ were related to the parametrization of vegetation, but this requires further investigation.
Our results at Johnston Draw demonstrate the dramatic effects of slope-aspect on soil climate. Although differential incoming solar radiation drives these effects, local conditions in the atmosphere, land cover, and soil modify those effects so that the seasonality of slope-aspect effects is very different from that predicted from S R,CS,T . We expect that slopeaspect effects at other locations will be similarly affected by local conditions, and we have shown that field monitoring and application of existing models can provide valuable information for determining the impacts of slope-aspect at those locations.

A C K N O W L E D G M E N T S
The research described in this paper was conducted in the Reynolds Creek Experimental Watershed, which is managed by the Northwest Watershed Research Center (NWRC) and part of the USDA-ARS. The NWRC supplied logistical support. Much of the field work and half of the instrumentation was supplied via grant, NSF CBET-0854553 and NSF EPS-0814387. In addition, we would like to thank Mark Murdock, Pat Kormos, and Ricki Loughridge for installation assistance.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.

D A T A AVA I L A B I L I T Y S T A T E M E N T
We are currently preparing data used in this study to be available via our website, www.ars.usda.gov/pacific-west-area/ boise.