Influence of climate variability on water partitioning and effective energy and mass transfer in a semi-arid critical zone

The critical zone (CZ) is the heterogeneous, nearsurface layer of the planet that regulates life-sustaining resources. Previous research has demonstrated that a quantification of the influxes of effective energy and mass transfer (EEMT) to the CZ can predict its structure and function. In this study, we quantify how climate variability in the last 3 decades (1984–2012) has affected water availability and the temporal trends in EEMT. This study takes place in the 1200 km upper Jemez River basin in northern New Mexico. The analysis of climate, water availability, and EEMT was based on records from two highelevation SNOTEL stations, PRISM data, catchment-scale discharge, and satellite-derived net primary productivity (MODIS). Results from this study indicated a decreasing trend in water availability, a reduction in forest productivity (4 g C m per 10 mm of reduction in precipitation), and decreasing EEMT (1.2–1.3 MJ m decade). Although we do not know the timescales of CZ change, these results suggest an upward migration of CZ/ecosystem structure on the order of 100 m decade, and that decadal-scale differences in EEMT are similar to the differences between convergent/hydrologically subsidized and planar/divergent landscapes, which have been shown to be very different in vegetation and CZ structure.


Introduction
The critical zone (CZ) is the surficial layer of the planet that extends from the top of the vegetation canopy to the base of aquifers (Chorover et al., 2011;Brandley et al., 2007).Within its boundaries, complex interactions between air, water, biota, organic matter, soils, and rocks take place that are critical for sustaining life on Earth (Brandley et al., 2007).The CZ has been conceptualized and studied as a weathering engine or reactor where interacting chemical, physical, and biological processes drive weathering reactions (Anderson et al., 2007;Chorover et al., 2011).Over long timescales, the CZ has evolved in response to climatic and tectonic forces and has been recently influenced by human activities (Steffen et al., 2007).Understanding how climate and land use changes affect CZ structure and related processes has become a priority for the scientific community due to the implications it may have on the functioning of life-supporting resources.It has been hypothesized by the researchers from the Jemez River basin (JRB) -Santa Catalina Mountains (SCM) critical zone observatory (CZO) (http://criticalzone. org/catalina-jemez/) that a quantification of the inputs of the effective energy and mass transfer (EEMT) to the CZ can provide insight about its structure and function (Chorover et al., 2011).CZ areas that receive greater EEMT influxes have been shown to have greater structural organization as well as more dissipative product removal (Rasmussen et al., 2011;Zapata-Rios et al., 2015a).The opposite has been observed in regions with less EEMT.

X. Zapata-Rios: Influence of climate variability on water partitioning and EEMT
EEMT is a variable that quantifies energy and mass transfer to the CZ (Rasmussen et al., 2011).EEMT integrates within a single variable the energy and mass associated with water that percolates into the CZ (E ppt ), and reduced carbon compounds resulting from primary production (E bio ) (Rasmussen et al., 2011).It has been demonstrated that other possible energy fluxes to the CZ, such as potential energy from transport of sediments, geochemical potential of chemical weathering, external inputs of dust, heat exchange between soil and atmosphere, and other sources of energy coming from anthropogenic sources, are orders of magnitude smaller (Phillips, 2009;Smil, 1991;Rasmussen et al., 2011;Rasmussen, 2012).Therefore the two dominant terms embodied in EEMT are E ppt and E bio , and only the energy associated with water and carbon is considered in the EEMT quantification.Energy from both water and net primary productivity is essential for CZ processes altering soil genesis, mineral dissolution, solute chemistry, and weathering rates (among others) (Birkeland, 1974;Neilson, 2003).
Previous research has shown that EEMT can become a tool to predict regolith depth, rate of soil production, and soil properties (Rasmussen et al., 2005(Rasmussen et al., , 2011;;Pelletier and Rasmussen, 2009a, b;Rasmussen and Tabor, 2007).For instance, strong correlations were found between EEMT, soil carbon, and clay content in soils on igneous parent materials from California and Oregon (Rasmussen et al., 2005).Furthermore, transfer functions were successfully determined between EEMT and pedogenic indices, including pedon depth, clay content, and chemical indices of soil alteration along an environmental gradient on residual igneous parent material (Rasmussen and Tabor, 2007).EEMT has also been incorporated in geomorphic and pedogenic models on granitic rocks to describe landscape attributes and regolith thickness (Pelletier and Rasmussen, 2009a, b).Rasmussen and Tabor (2007) demonstrated that regolith depth on stable low-gradient slopes increased exponentially with increasing EEMT.Similarly, Pelletier et al. (2013) found that high EEMT values are associated with large above-ground biomass, deeper soils, and longer distance to the valley bottoms across hillslopes in the Santa Catalina Mountains in southern Arizona.More recently, EEMT estimations have been strongly correlated with water transit times, water solute concentrations, and dissolution of silicates on a rhyolitic terrain in northern New Mexico (Zapata-Rios et al., 2015a).In these studies, the main constituents of EEMT (E ppt and E bio ) were quantified as an average value, based on climate records from long-term regional databases as these variables exert first-order controls on photosynthesis and effective precipitation (Rasmussen et al., 2011;Chorover et al., 2011).
It is still uncertain how climate variability influences CZ structure, function, and the timescales of these changes (Chorover et al., 2011;Brooks et al., 2015).Climate variability might directly influence changes in the transfer of mass and energy to the CZ as climate has a direct control on both E ppt and E bio .In the mountains of the southwestern United States, a large percentage of annual precipitation falls as snow, which is stored during the winter and released as snowmelt during the spring (Clow, 2010).The water from the winter snowpack constitutes the main source of regional water supplies and the largest component of runoff (Bales et al., 2006;Nayak et al., 2010).The regional snowpack has been documented to be declining in the southwestern US (Mote et al., 2005;Clow, 2010) and alterations to the snowpack are likely to produce changes in vegetation, impact water availability (Bales et al., 2006;Harpold et al., 2012;Trujillo et al., 2012), and influence inputs of EEMT.For instance, significant increasing trends in air temperature and decreasing trends in winter precipitation in the last decades have been documented in the upper Rio Grande region in northern New Mexico (Harpold et al., 2012).
The objective of this study was to evaluate climate variability and its influence on the temporal changes of water partitioning and EEMT at the catchment scale in a semi-arid CZ over the last few decades.This investigation took place in the upper part of the JRB in northern New Mexico, a basin dominated by a forest cover and limited human infrastructure.Micro-climate variability was studied based on daily records from two SNOTEL stations using records from 1984 through 2012.Water availability and EEMT were estimated during the same time period, based on precipitation and temperature from the Parameter-elevation Regressions on Independent Slopes Model (PRISM), empirical daily observations of catchment-scale discharge, and satellite-derived net primary productivity (MODIS).

Study site
The Jemez River is a tributary of the upper reach of the Rio Grande and is located between Jemez and the Sierra Nacimiento mountains in northern New Mexico (Fig. 1a).Its headwaters originate within the 360 km 2 Valles Caldera National Preserve which contains 30 % of the total basin surface (Fig. 1b).The upper JRB is located at the southern margin of the Rocky Mountain ecoregion between latitudes 35.6 and 36.1 • N and longitudes −106.3 and −106.9W. The basin is characterized by a mean elevation of 2591 m and a gradient in elevation ranging from 1712 to 3435 m.Based on a 10 m digital elevation model, the catchment drains 1218 km 2 above the US Geological Survey (USGS) gauge "Jemez River near Jemez" (35.66 • N and 106.74 • W; USGS 08324000), located at an elevation of 1712 m.The basin has a predominant south aspect and a mean catchment slope of 13.7 • .The geology consists of rocks of volcanic origin with predominant andesitic and rhyolitic compositions that overlie tertiary to Paleozoic sediments along the western margin of the Rio Grande rift (Shevenell et al., 1987)   Inceptisols (Allen et al., 1991(Allen et al., , 2002)).Precipitation has a bimodal pattern with 50 % of annual precipitation occurring during the winter months (primarily as snow) from October to April and originates from westerly frontal systems.The remaining 50 % of precipitation falls as convectional rainfall during the monsoon season between July and September (Sheppard, 2002).According to the National Land Cover Database (NLCD), the basin is a forested catchment with 79 % under evergreen, deciduous, and mixed forest cover and only 0.5 % of area covered by development and agriculture (http://www.mrlc.gov/nlcd06_leg.php)(Table 1).

Climatological stations
There

Climate variability
Climate variability was studied based on 13 variables from the two SNOTEL stations, derived from daily air temperature, precipitation, and maximum SWE, following a similar methodology and data processing procedure as in Harpold et al. (2012).The variables analyzed were winter, summer, and annual air temperature ( (water year day); last day snow cover (water year day); length of snow on the ground (number of days); and SM50, which is the day of the year in which half of the snowpack melts (number of days).Climate records for data analysis were aggregated by water year (from 1 October to 30 September).
Winter season was considered to be between October and April and summer season between May and September.The analysis of climate was conducted from 1984 as a starting year to avoid the anomalous wet years recorded at the beginning of 1980s that were caused by the Pacific Decadal Oscillation (PDO) and El Niño-Southern Oscillation (ENSO) (Harpold et al., 2012;and references therein).The presence of a monotonic increasing or decreasing trend in the 13 climate variables recorded at the two individual stations was evaluated from 1984 through 2012 by applying the nonparametric Mann-Kendall test (MKT) with a α = 0.10 level of significance and the nonparametric Sen's slope estimator of a linear trend (Yue et al., 2012;Sen, 1968).

EEMT estimation
Energy from both water and net primary productivity are essential on CZ processes altering soil genesis, mineral dissolution, solute chemistry, and weathering rates (among others) (Birkeland, 1974;Neilson, 2003;Anderson et al., 2007).In this investigation EEMT was calculated as the sum of E ppt and E bio (Eq.1).We applied two different methods to estimate E ppt and E bio .Following a similar methodology described in Rasmussen and Gallo (2013), the term EEMT emp was empirically estimated at the catchment scale based on baseflow estimations and average basin-scale net primary productivity (NPP) derived from MODIS satellite data.In comparison, EEMT model was estimated at the catchment scale based on long-term climate records from PRISM, developed by the climate group at Oregon State University (http://www.wcc.nrcs.usda.gov/ftpref/support/climate/prism/) and described in Rasmussen et al. (2005Rasmussen et al. ( , 2011)).PRISM is a weighted regression technique that accounts for physiographic factors affecting climate variables, and it has been extensively used in the US (Daly et al., 1994(Daly et al., , 2002)).The assumption of this study is that the 800 m PRISM data will provide a reasonable spatial estimation of basin-scale precipitation.

EEMT emp
Upper JRB precipitation and air temperature from 1984 through 2012 was obtained using PRISM data at an 800 m spatial resolution (Daly et al., 1994(Daly et al., , 2002)).Daily discharge data were available from 1984 through 2012 from the USGS Jemez River near Jemez gauge station (http://waterdata.usgs.gov/nwis).The upper Jemez River has not been subjected to flow regulation, and almost 60 % of the annual discharge occurs during the snowmelt period between March and May.
Daily discharge records were normalized by catchment area, and mean daily discharge was aggregated into water years.Precipitation (P ) on the land surface was partitioned between quickflow (S) and catchment wetting (W ).S represents water that directly contributes to streamflow discharge as a response to precipitation events; thus this amount of water is not transferred to the CZ.W is the total amount of water that infiltrates the soil, of which a portion is available for vaporization (V ) including vegetation uptake.The remaining portion of W flows though the CZ and contributes to baseflow (U ).V was estimated at the annual scale as the difference between P and discharge (Q).Q was separated between S and U using a one-parameter low-pass filter (Lyne and Hollick, 1979;Arnold and Allen, 1999;Eckhardt, 2005;Troch et al., 2009) where a is a filter parameter set to 0.925.This filter was passed twice, backward and forward in time, to improve the partitioning of U and S at the beginning of the time series.After this, daily values of Q, U , and S were integrated to annual timescales.Alterations in snowmelt timing were evaluated with Q 50 , which indicates the day of the water year when 50 % of the total annual discharge is recorded at the catchment outlet (Clow, 2010;Stewart et al., 2004).
The term E ppt_emp (energy input through precipitation) was calculated as stated in Eq. (3) based on estimations of U and mean PRISM-derived air temperature at the catchment scale (Rasmussen et al., 2011;Rasmussen and Gallo, 2013).
In Eq. ( 3), C w is the specific heat of water (4187 J kg −1 K −1 ) and T is the difference in temperature between ambient temperature and 0 • C, calculated as

Net primary productivity
Mean annual NPP at the catchment scale was estimated at a 1 km spatial resolution for the years 2000-2012 using data MOD17A3 from MODIS (Zhao and Running, 2010) (http: //modis-land.gsfc.nasa.gov/npp.html).E bio was calculated as indicated in Eq. ( 4) and presented in Rasmussen et al. (2011) and Rasmussen and Gallo (2013).

EEMT model
E ppt_model was calculated based on estimations of effective precipitation (P eff ) which is defined as the amount of water that enters the CZ in excess of evapotranspiration and is available to flow through the CZ (Rasmussen et al., 2005;Eq. 5).
where P eff(i) is the monthly effective precipitation calculated as the difference between monthly PRISM precipitation and monthly potential evapotranspiration, calculated using the Thornthwaite equation (Rasmussen et al., 2005;Thornthwaite, 1948).P eff , calculated as the difference between monthly precipitation and potential evapotranspiration, has been traditionally used in soil water balances (Arkley, 1963).
C w and T are the same parameters described in Eq. ( 3).E ppt(i) model was calculated on a monthly basis only for the months when precipitation is larger than evapotranspiration (P eff(i) >0), and these values were integrated in water years.E bio_model was estimated as indicated in Eq. ( 4), and NPP was calculated following an empirical relationship based on air temperature (Eq.6; Lieth, 1975).
NPP (i) is the monthly NPP in g m −2 yr −1 , and T a is the monthly air temperature.Days (i) over the number of days in a year is an NPP time correction.Similar to Eq. ( 5), E bio_model was calculated only for the months where P eff(i) >0.For a detailed description of EEMT, see Rasmussen et al. (2005Rasmussen et al. ( , 2011Rasmussen et al. ( , 2015)), Rasmussen and Tabor (2007), and Rasmussen and Gallo (2013).
The EEMT model quantification presented in Chorover et al. ( 2011) has a relative mean prediction error of ∼ 25 % relative to the predicted value.However, we are using mean trends in EEMT at this catchment-scale study, so we believe that even though the EEMT calculation may have errors, the mean trends presented in this investigation are close to the true values.

Water availability, water partitioning, and climate controls on water availability
A trend analysis was conducted using data from 1984 through 2012 on each component of the water partitioning analysis (P , Q, U , S, V , W , Q 50 ), and EEMT using the nonparametric MKT and the Sen's slope estimator of a linear trend with a α = 0.10 level of significance (Yue et al., 2012;Sen, 1968).Relationships between climate, hydrological variables, and EEMT were examined by simple and multiple linear regression analysis with parameters fit through a least-square iterative process (Haan, 1997).

Climate variability
Records from the Quemazon SNOTEL station from 1984 to 2012 indicated a mean annual precipitation of 701 mm, of which 50 % fell during the winter months with a mean maximum SWE of 242.5 mm.The mean annual and winter temperatures at this site were 3.98 and −0.87  2).
During the 3 decades of analysis, 7 out of the 13 climate variables in both stations showed a statistically significant trend (Table 3).Mean winter, summer, and annual air temperatures at the Quemazon station increased significantly by 1.3 • C (p < 0.001), 1.0 • C (p < 0.01) and 1.4 • C decade −1 (p < 0.001), respectively.Similarly, the same variables at the Señorita Divide no. 2 station increased 1.0 • C (p < 0.05), 1.0 • C (p < 0.01), and 1.2 • C (p < 0.001) decade −1 , respec- Mean annual precipitation at the catchment scale correlated significantly with the mean annual precipitation recorded at the Quemazon (R 2 = 0.45; p < 0.0001) and Señorita Divide no. 2 stations (R 2 = 0.73; p < 0.0001).In this same timeframe average, minimum and maximum basin-scale temperatures were 6.1, −1.5, and 13.6 • C, respectively.In general, January was the coldest and July was the warmest month.Basin-scale mean annual and winter temperatures indicated a statistically significant increasing trend of 0.5 and 0.4 • C decade −1 (not shown).Mean annual temperature in the JRB significantly correlated with the mean annual temperature recorded at the Quemazon (R 2 = 0.29; p < 0.006) and Señorita Divide no. 2 stations (R 2 = 0.67; p < 0.0001) (not shown).
Mean river basin discharge during the study period was 0.15 mm day −1 , and the maximum and minimum historical streamflow discharges were 2.97 and 0.008 mm day −1 , respectively.In the 29 years of daily discharge records, 90 % of the time discharge surpassed 0.03 mm day −1 and 10 % of the time exceeded 0.38 mm day −1 .Peak discharge occurred between March and May and 58 % of the annual discharge flowed between these months.From 1984 to 2012, 3 % of annual precipitation became quickflow and contributed directly to the streamflow discharge (3 % P ; standard deviation SD = 1.2 % P ).As a result, 97 % of the annual precipitation (SD = 1.2 % P ) infiltrated and was available for vegetation uptake.This 97 % of annual precipitation is further partitioned between vaporization and baseflow.The amount of water vaporized into the atmosphere represented 91 % of the annual precipitation (SD = 3.4 % P ).Baseflow corresponded to 6.1 % of the annual precipitation (SD = 2.2 % P ) and represented the largest component of discharge (73.2 % Q; SD = 5.4 % Q).Quickflow represented the remaining 26.8 % of annual discharge (SD = 5.4 % Q).

EEMT emp
Using the available 2000-2012 remote-sensing data, mean MODIS NPP was found to be 450 g C m −2 (SD = 57.1 g C m −2 ).Using these 13 years of data, no trend in the mean annual NPP for the upper JRB was found.However, mean annual NPP was positively correlated with basin-scale precipitation (R 2 = 0.56; p = 0.003) and baseflow (R 2 = 0.41; p = 0.02) (Fig. 3).These results indicated that forest productivity in the upper JRB is primarily limited by water availability, since other climate variables recorded at the two SNOTEL stations were not good predictors of NPP.As with any spatial and temporal regression between climate and MODIS data, there are potential errors associated with forest disturbance, interannual lag effects, and interseason variability of water availability and other factors.We also note that the significant relationship, albeit with variability and error, likely captures these effects on this timescale of the study when no large-scale disturbance occurred.
From 1984 through 2012, mean E ppt_emp was 1.03 MJ m 2 yr −1 (SD = 0.49 MJ m 2 yr −1 ) and mean E bio_emp was 9.89 MJ m 2 yr −1 (SD = 1.26 MJ m 2 yr −1 ).Multivariate regression analysis indicated that precipitation at the Quemazon station and the upper JRB were the best predictors of E bio_emp (R 2 = 0.66; p = 0.06).Using this multivariate linear regression model, E bio_emp data were extrapolated for the years 1984-1999.Using the combined data set from extrapolated and measured E bio_emp the mean annual E bio_emp was 10.8 MJ m 2 yr −1 (SD = 1.37 MJ m 2 yr −1 ) for the period from 1984 to 2012.Mean EEMT emp was 11.83 MJ m 2 yr −1 (SD = 1.74 MJ m 2 yr −1 ) and E bio_emp represented 92 % (SD = 0.03 %) of the total EEMT emp during the study period.
EEMT emp was, on average, 1.7 times larger than EEMT model .Both EEMT emp and EEMT model showed a significant linear correlation (R 2 = 0.42; p = 0.0002) and a similar decreasing trend of 1.2 MJ m 2 decade −1 (p ≤ 0.01) and 1.3 MJ m 2 decade −1 (p ≤ 0.05), respectively (Fig. 4).Detailed estimations of EEMT emp and EEMT model and their components can be found in Table S1 (in the Supplement).Figure 5 highlights changes of EEMT in the upper JRB in relation to water availability from 1984 to 2012.EEMT was positively correlated to annual baseflow, increasing during wet years and decreasing during dry years.

Climate variability
Global climate is changing and the instrumental records in the southwestern US for the last 3 decades indicate a decline in precipitation and increasing air temperatures in the region (Hughes and Diaz, 2008;Folland et al., 2001).Global climate models further predict drier conditions and a more arid climate for the 21st century in this region (Seager et al., 2007).For instance, according to low-and high-emission scenarios, global climate models indicate a substantial increase in air temperature between 0.6 and 2.2 • C and 1.3 and 5.0 • C for the period 2021-2050 and by end of the 21st century, respectively (Barnett et al., 2004;Cayan et al., 2013).An increase in winter temperature of about 0.6 • C decade −1 was reported from 1984 to 2012 at a regional level in the upper Rio Grande basin (Harpold et al., 2012).In line with these other studies, we found that mean annual and winter air temperature in the upper JRB have increased 0.5 and 0.4 • C decade −1 , respectively.
Changes in climate have been found to be a predominant influence in snowpack decline as opposed to changes in land use, forest canopy, or other factors (Hamlet et al., 2005;Boisvenue and Running, 2006).There are high-confidence predictions that snowpacks will continue to decline in northern New Mexico through the year 2100, and projections of snowpack accumulation for mid-century (2041-2070) show a marked reduction for SWE of about 40 % (Cayan et al., 2013).Harpold et al. (2012) found a decrease in annual precipitation and maximum SWE for the upper Rio Grande basin of −33 and −40 mm decade −1 , respectively.In this study, a clear decreasing trend in annual, winter precipitation, and maximum SWE was observed in records from 1984 to 2012 in the two high-elevation SNOTEL stations.Records in this study showed approximately twice the rate of decrease in annual precipitation and a smaller decrease in maximum SWE of about 7 mm decade −1 , compared to the regional results from Harpold et al. (2012).Harpold et al. (2012) report that SM50 (−2 days decade −1 ), snow cover length (−4.2 days decade −1 ), day of maximum SWE (−3.31 days decade −1 ), and last day of snow cover (−3.45 days decade −1 ) for the Rio Grande basin showed statistically significant trends.However, based on our analysis from the individual SNOTEL stations, these variables did not show any statistically significant trends.

Changes in discharge and evapotranspiration
Decreasing trends in discharge ranging from 10 to 30 % are expected during the 21st century for the western US (Milly et al., 2005), and maximum peak streamflow is expected to happen 1 month earlier by 2050 (Barnett et al., 2005).Furthermore, it has been reported that streamflow in snowmeltdominated river basins is more sensitive to wintertime increases in temperature (Barnett et al., 2005).In this study, we have found that 50.5 % of annual streamflow occurred between (April) and beginning of the summer (June).This result is congruent with other studies in snowmelt-dominated systems in the region (Clow, 2010).Previous research in the southwest has found that the timing of snowmelt is shifting to early times ranging from a few days to weeks (Stewart et al., 2004;Mote et al., 2005;McCabe and Clark, 2005).For instance, Clow (2010) reports that in southern Colorado rivers, there is a trend toward earlier snowmelt that varied from 4.0 to 5.9 days decade −1 and 1 April SWE decreased between 51 and 95 mm per decade.In this study, it was found that snowmelt timing in the upper JRB occurred 4.3 days decade −1 earlier and 1 April SWE decreased between 54 and 60 mm decade −1 .
The spatial and temporal variability in total evapotranspiration may exhibit significant variability (Tague and Peng, 2013), and contrasting evapotranspiration trends' directions have been reported in different studies around the world (Barnett et al., 2005).In the JRB, a snow-dominated system, the decrease in vaporization (45 mm decade −1 ) is likely a result of the mismatch of the timing of energy and water fluxes.While plant water demand remains relatively low, earlier snowmelt may reduce evapotranspiration by reducing plant/atmospherically available water later during the growing season when demand is higher (Barnett et al., 2005).The decrease in vegetation biomass related to water availability, indicated from the MODIS data at this basin, can also significantly contribute to alter transpiration water losses.An increase in forest water-use efficiency (ratio of water loss to carbon gain) with increasing concentrations of carbon dioxide can also contribute as another cause to the decrease of evapotranspiration fluxes (Keenan et al., 2013).Modeling studies conducted over 100 years support our finding that evapotranspiration has been decreasing in the west arid area of the US (Liu et al., 2013).However, evapotranspiration may increase with temperature in some snow-dominated systems if stored soil or groundwater remains available to plants, either locally or at downslope locations (Goulden et al., 2012;Brooks et al., 2015).

Forest productivity
Reduced carbon compounds resulting from primary production are a fundamental energy component of EEMT (Rasmussen et al., 2011).Modeling and empirical studies indicate that mountain forest productivity in the southwest is sensitive to water and energy limitations (Christensen et al., 2008;Tague et al., 2009;Anderson-Teixeira et al., 2011;Zapata-Rios et al., 2015b;Zapata-Rios, 2015).Trujillo et al. (2012) found that the satellite derived Normalized Difference Vegetation Index (NDVI) greening increased and decreased proportionally to the changes in snowpack accumulation along a gradient in elevation in the Sierra Nevada, while Zapata-Rios et al. (2015b) and Zapata-Rios (2015) found similar results across a gradient of energy created by aspect differences at higher elevations in the Jemez Mountains.Furthermore, energy limitations to productivity have been observed in colder sites at high elevations (Trujillo et al., 2012;Anderson-Teixeira et al., 2011;Zapata-Rios et al., 2015b;Zapata-Rios, 2015).Since the mid-1980s increases have been documented in wildfires and tree mortality rates in highelevation forests due to an increase in spring and summer temperatures and decrease in water availability (Westerling et al., 2006;Van Mantgem et al., 2009).Results from this study indicated that in the upper JRB, forest productivity was primarily responding to water availability (Fig. 3).

EEMT variability
All of the above results indicate that the JRB is highly susceptible to changes in climate that can affect water availability and ecosystem productivity which impacts EEMT.Rasmussen et al. (2005) estimated low rates of EEMT model < 15 MJ m −2 yr −1 for the majority of the continental US and demonstrated that E bio was the dominant component of EEMT, with contributions above 50 % of total EEMT in soil orders associated with arid and semi-arid regions.Regions dominated by E bio corresponded to regions facing water limitation and where E bio accounted for up to 93 % of the total energy and carbon flux to the CZ (Rasmussen et al., 2011;Rasmussen and Gallo, 2013).In semiarid regions vaporization represents over 90 % loss of annual precipitation (Newman et al., 2006), while groundwater recharge accounts for less than 10 % of annual precipitation (Scanlon et al., 2006).Under these conditions, little water remains for CZ processes in semi-arid regions.Other studies have found that the contributions of E bio can be 3-7 orders of magnitude larger than other sources of energy influxes to the CZ (Phillips, 2009;Amundson et al., 2007).In this study, we confirmed that for the upper JRB, E bio was the dominant term from the total EEMT, and E ppt contributions were small.
A comparison of EEMT model and EEMT emp in 86 catchments across the US, characterized by having minimum snow influence, indicated that model and empirical values X. Zapata-Rios: Influence of climate variability on water partitioning and EEMT were strongly linearly correlated (R 2 = 0.75; p < 0.0001) and EEMT model values were larger than EEMT emp (Rasmussen and Gallo, 2013).One limitation of the EEMT model method is that it calculates energy during the months when air temperature is above 0 only and assumes no energy associated with precipitation falling as snow.In a snowmeltdominated system like the upper JRB, where snowmelt is the main source of water availability to ecosystems (Bales et al., 2006), EEMT estimations based only on climate data will likely underestimate the energy transfer to the CZ.Therefore, using EEMT emp methodology may be more suitable for snowmelt-dominated systems.In this study, we found the expected linear correlation between EEMT model and EEMT emp (R 2 = 0.42; p < 0.001); however, EEMT model values were smaller than EEMT emp values.Although the two methods used in this study to calculate EEMT indicated different absolute values of EEMT, the rates of decrease of EEMT per decade are congruent with each other (EEMT emp = 1.2 MJ m 2 decade −1 ; EEMT model = 1.3 MJ m 2 decade −1 ) (Fig. 5).
While the correlation between EEMT and CZ landscape structure does not necessitate causation, previous work has shown that these correlations are widespread and strong, and thus, EEMT has significant predictive ability (Pelletier and Rasmussen, 2009a, b;Rasmussen and Tabor, 2007;Rasmussen et al., 2005Rasmussen et al., , 2011;;Pelletier et al., 2013;Zapata-Rios et al., 2015a).Although we do not know the exact timescale of CZ change (Brooks et al., 2015), we believe the rates of EEMT change found in the upper JRB between 1.2 and 1.3 MJ m 2 decade −1 can be significant for CZ processes.These rates of EEMT change could represent an upward movement of more arid, lower EEMT systems to higher elevations.For instance, in a study conducted in a similar semiarid region in the Santa Catalina Mountains located in southern Arizona, Rasmussen et al. (2015) estimated differences in EEMT of about 25 MJ m 2 yr −1 between the upper elevation (2800 m) covered by mixed conifer forest and low elevation (800 m) covered by a dry semi-arid desert scrub ecosystem.These changes in EEMT along the 2000 m elevation gradient in the SCM are equivalent to a difference of 1.25 MJ m 2 yr −1 per 100 m in elevation change.The rates of EEMT change, every 100 m along the SCM elevation gradient, are similar to observed rates of EEMT change per decade for the entire JRB.Along this elevation gradient contrasting vegetation, soil characteristics, regolith development, chemical depletion, and mineral transformation have been observed between lower and high elevations on similar granitic parent material (Whittaker et al., 1968;Lybrand et al., 2011;Lybrand and Rasmussen, 2014;Holleran et al., 2015).Mollisols and carbon-rich soils have been characterized in convergent areas of higher EEMT versus weakly developed Entisols in lower EEMT landscape positions (Lybrand et al., 2011;Holleran et al., 2015).Furthermore, Rasmussen et al. (2015) determined differences of 3.9 MJ m 2 yr −1 between contrasting north-and south-facing slopes at a similar elevation.In areas with similar EEMT, north-facing slopes have soils characterized by greater clay and carbon accumulation (Holleran et al., 2015).According to topographic wetness, differences of 0.9 MJ m 2 yr −1 were determined between water-gaining and water-losing portions of the landscape (Rasmussen et al., 2015).
It is still uncertain how the CZ evolves over time and how climate, lithology, and biota influence the function of the CZ (Chorover et al., 2011).We postulated that a measure of the energy inputs into the CZ drive CZ evolution, and their quantification can be related to functions and processes within the CZ.The energy inputs and mass transfer have been integrated in a single and transferable metric (EEMT), quantified as water and carbon fluxes that can be easily transferred and quantified in different ecosystems and regions around the world (Rasmussen and Tabor, 2007;Rasmussen et al., 2011).This allows the comparison of energy inputs to the CZ in a broad range of sites, climates, and ecosystems.EEMT can be used as a tool to provide an initial identification of landscape locations subjected to higher energy influx (as a result of water and reduced carbon throughputs) or locations where EEMT is changing over time, as it has been indicated in the present study.Consistent changes in EEMT can be indicators of alteration in the function of the CZ, such as weathering process, hydrochemical, and hydrologic response (among others).In regions where temperature, precipitation, water availability, and vegetation are changing, a quantification of EEMT can provide an initial assessment and a metric to evaluate changes in the CZ.The EEMT model has a limitation in that it does not provide information on how energy is distributed within the CZ and does not provide mechanistic insight into CZ processes.However, it can be used to identify research sites for further instrumentation and measuring of CZ processes.Although the quantification of EEMT using the methodologies applied in this study is suitable for large spatial scales, it is limited in that it does not take into account small-scale variabilities induced by topography in solar energy, effective precipitation, NPP, and redistribution of water by differences in microtopography.Therefore, EEMT estimations at small scales (pedon to hillslopes) need to follow a different approach, as indicated in Rasmussen et al. (2015).

Summary
We investigated how changes in climate in the southwest affect the trends in water availability, vegetation productivity, and the annual influxes of EEMT to the CZ.This investigation took place in the 1200 km 2 upper JRB, a semi-arid basin in northern New Mexico, using records from 1984 to 2012.Results at the two SNOTEL stations indicated clear increasing trends in temperature and decreasing trends in precipitation and maximum SWE.Temperature changes include warmer winters (+1.0-1.3 • C decade −1 ), and generally warmer year-round temperatures (+1.2-1.4 • C decade −1 ).
Precipitation changes include a decreasing trend in precipitation during the winter (−41.6-51.4mm decade −1 ), during the year (−69.8-73.2mm decade −1 ), and maximum SWE (−33.1-34.7 mm decade −1 ).At the upper JRB, all the water-partitioning components showed statistical significant decreasing trends including precipitation (−61.7 mm decade −1 ), discharge (−17.6 mm decade −1 ), and vaporization (−45.7 mm decade −1 ).Similarly, Q 50 , an indicator of snowmelt timing, is occurring −4.3 days decade −1 earlier.Basin-scale precipitation (R 2 = 0.56; p = 0.003) and baseflow (R 2 = 0.41; p = 0.02) were the strongest controls on NPP variability, indicating that forest productivity in the upper JRB is water limited.This study showed a positive correlation between water availability and EEMT.For every 10 mm of change in baseflow, EEMT varies proportionally in 0.6-0.7 MJ m −2 yr −1 .From 1984 to 2012, changes in climate, water availability, and NPP have influenced EEMT in the upper JRB.A decreasing trend in EEMT of 1.2-1.3MJ m −2 decade −1 was calculated in this same time frame.Although we cannot determine the timescales of change, these results suggest an upward migration of CZ/ecosystem structure on the order of 100 m decade −1 , and that decadal-scale differences in EEMT are similar to the differences between convergent/hydrologically subsidized and planar/divergent landscapes, which have been shown to be very different in vegetation and CZ structure.As the landscape moves towards a drier and hotter climate, changes in EEMT of this magnitude are likely to influence CZ processes.
Author contributions.All authors contributed extensively to this research.All authors discussed the methodology, results, and commented on the manuscript at all stages.Xavier Zapata-Rios analyzed data and prepared the manuscript with contributions from all co-authors.

Figure 1 .
Figure 1.(a) Relative location of study area within the northwestern state of New Mexico, (b) upper JRB, ∼ 1200 km 2 , delimited above the USGS gauge station "Jemez River near Jemez" (USGS 08324000) based on a 10 m digital elevation model (DEM).

Figure 2 .
Figure2.Precipitation and water partitioning at the upper Jemez River catchment scale.There was a significant decreasing trend quantified by the Mann-Kendall test (MKT) in the JRB precipitation and all the components of the water partitioning.For instance, precipitation at the catchment scale decreased during the last 3 decades at a rate of 6.17 mm per year and discharge at 1.76 mm yr −1 .Q 50 indicated that discharge is occurring 4.3 days earlier per decade.

Figure 3 .
Figure 3. (a) Positive linear correlation between precipitation in the upper JRB and annual NPP in the upper JRB-derived from MODIS; (b) linear correlation between baseflow and annual NPP in the upper JRB.Forest productivity is water limited in the upper JRB.Other variables such as annual, winter, and summer air temperature did not correlate with NPP.

Figure 4 .
Figure 4. (a) EEMT emp and EEMT model showed similar significant decreasing trends from 1984 to 2012 of 1.2 and 1.3 MJ m −2 yr −1 ; (b) EEMT emp and EEMT model showed a significant linear correlation.

Figure 5 .
Figure 5. Relationship between water availability and EEMT.Baseflow and EEMT showed a positive linear correlation.As water availability in the JRB decreases, indicated by baseflow, EEMT also decreases.

Table 1 .
Land use classification of the JRB area.79.7 % of the total basin is covered by forest, according to the National Land Cover Database (NLCD; http://www.mrlc.gov/nlcd06_leg.php).

Table 2 .
Site and meteorological information for the SNOTEL Quemazon and Señorita Divide no. 2 stations located at high elevations in the upper part of the JRB.
• C, respectively.During the same time period, Señorita Divide no. 2 station had a mean annual precipitation of 686 mm, of which 61 % fell during the winter, with a mean maximum SWE recorded of 239.2 mm.The mean annual and winter temperatures at the Señorita Divide no. 2 site were 4.23 and −0.90 • C, respectively (Table

Table 3 .
Climatic time series trends for the Quemazon and Señorita Divide no. 2 SNOTEL stations from 1984-2012.A trend in the precipitation time series was evaluated with the MKT and Sen's slope estimator.Trends were considered statistically significant at p ≤ 0.1.The results showed an increasing trend in winter, summer, and annual temperature in the two stations.Annual and winter precipitation, maximum SWE, and 1 April SWE decreased in both stations during the 29 years analyzed.The last day of snow cover decreases significantly only at the Quemazon station.No significant trend was observed for the SWE : winter P ratio, duration of snowmelt SM50, and length of snow on the ground.There was no significant trend in the ratio between SWE to winter precipitation at either station.Observed 1 April SWE also decreased −60.5 mm decade −1 (p ≤ 0.05) and −54.4 mm decade −1 (p ≤ 0.1) at the Quemazon and Señorita Divide no. 2 stations, respectively.The day of occurrence of maximum SWE recorded at the Quemazon station showed a significant trend, indicating that maximum SWE is occurring 5.7 days earlier every decade (p ≤ 0.05).However, this same trend was not observed at the Señorita Divide no. 2 station.Variables such as SM50, initiation of snow cover, and snow cover duration did not indicate any trend of change in either station at the 90 % confidence level.
a Statistical significance; * P < 0.1; * * P < 0.05; * * * P < 0.01; * * * * P < 0.001.tively.The rates of increase in winter and annual air temperatures were larger in Quemazon, which was the higher elevation station.Annual precipitation decreased in both stations at similar rates per decade.Quemazon station annual precipitation decreased 69.8 mm decade −1 (p ≤ 0.01) and Señorita Divide no. 2 decreased 73.2 mm decade −1 (p ≤ 0.05).Winter precipitation decreased faster at the Señorita Divide no. 2, the lower elevation station (59.4 mm decade −1 ; p ≤ 0.05) than at the Quemazon station (41.6 mm decade −1 ; p ≤ 0.1).Maximum SWE decreased in both stations at similar rates, −34.7 mm decade −1 at Señorita Divide no. 2 and −33.1 mm decade −1 at the Quemazon station (p ≤ 0.1).In contrast, there is a decreasing trend in the last day of snow cover, which is happening about 6 days sooner per decade in the Quemazon station (p < 0.05).Last day of snow cover at the Señorita Divide no. 2 station did not show a significant trend (Table3).