Gridded Snow Water Equivalent Reconstruction for Utah Using Forest Inventory and Analysis Tree-Ring Data

Snowpack observations in the Intermountain West are sparse and short, making them difficult for use in depicting past variability and extremes. This study presents a reconstruction of April 1 snow water equivalent (SWE) for the period of 1850–1989 using increment cores collected by the U.S. Forest Service, Interior West Forest Inventory and Analysis program (FIA). In the state of Utah, SWE was reconstructed for 38 snow course locations using a combination of standardized tree-ring indices derived from both FIA increment cores and publicly available tree-ring chronologies. These individual reconstructions were then interpolated to a 4-km grid using an objective analysis with elevation correction to create an SWE product. The results showed a significant correlation with observed SWE as well as good correspondence to regional tree-ring-based drought reconstructions. Diagnostic analysis showed statewide coherent climate variability on inter-annual and inter-decadal time-scales, with added geographical details that would not be possible using courser pre-instrumental proxy datasets. This SWE reconstruction provides water resource managers and forecasters with better spatial resolution to examine past variability in snowpack, which will be important as future hydroclimatic variability is amplified by climate change.


Introduction
Water supplies in the U.S. Intermountain West will be affected by a decrease in snowpack compounded by an amplified fluctuation of snowpack variability [1][2][3][4][5]; this tendency may be particularly pronounced in the topographically diverse state of Utah [6]. To further complicate water issues, Utah has a rapidly growing population that is currently straining water resources. Utah's climate can be characterized as warm to hot-humid continental climate type, with hot, dry summers and cold winters [7]. Precipitation amounts are typically highest in the spring and lowest in the summer, but winter storms bring copious amounts of snow to the mountains. This leads to the development of heavy snowpack that is a critical form of water storage for the state. During the winter months (October thru March), nearly 70% of mountain precipitation arrives as snow [8]. This accumulation of snowpack peaks near April 1 and then melts gradually over the summer, and is an important source of runoff. However, the snow water equivalent (SWE) of accumulated snowpack varies widely by location within the state. While 1000 mm is typical, in northern Utah snowpack measurement stations sometimes exceed 1500 mm, while in the southern portion of the state 500 mm is usual. Recent research has suggested that Utah snowpack is declining [6], and understanding the climate controls on past snowpack will help water managers plan for the future. Furthermore, understanding the temporal and spatial variability of SWE prior to instrumentation will improve water resource management and forecasting.
While numerous studies have indicated the possibility for both decreased snowpack and increased variability in snowpack accumulation across the west, subsequent changes in the timing and distribution of peak and annual runoff are also likely [9,10]. However, a relatively short period of instrumental record (starting approximately in 1930) in the Intermountain West limits our understanding of snowpack variability. Fortunately, in the western U.S. tree-ring data have been used to create a broad range of climate reconstructions including precipitation (e.g., [11,12]), temperature (e.g., [13]), and streamflow (e.g., [14][15][16]), as well as snowpack for parts of the Rocky Mountains [17][18][19][20], and could be employed for high resolution reconstructions of SWE.
Tree-ring-based climate reconstructions rely on intensive sampling of trees that have been carefully selected from site(s) and/or species thought to record a climatic variable of interest and can be replicated [21]. Subsequently, statistical reconstruction models rely on the presence and length of tree-ring chronologies across a given region. There have been important efforts to create gridded datasets with broad spatial coverage using tree-ring chronologies, notably for temperature and precipitation [22], as well as the Palmer Drought Severity Index (PDSI) [23]. These gridded products have proven to be useful for the analysis of long-term climate variability at regional to continental scales. These datasets become more limited, however, when it comes to analyzing watershed-scale variability within complex terrain such as that found in the Intermountain West (see Figure 1). The production of multiple tree-ring chronologies that would enable fine-scale and broad-coverage reconstructions is desirable. However, chronologies are extremely time-consuming and cost-prohibitive to develop, thus alternative data are also highly desirable.
Water 2017, 9,403 2 of 13 is usual. Recent research has suggested that Utah snowpack is declining [6], and understanding the climate controls on past snowpack will help water managers plan for the future. Furthermore, understanding the temporal and spatial variability of SWE prior to instrumentation will improve water resource management and forecasting. While numerous studies have indicated the possibility for both decreased snowpack and increased variability in snowpack accumulation across the west, subsequent changes in the timing and distribution of peak and annual runoff are also likely [9,10]. However, a relatively short period of instrumental record (starting approximately in 1930) in the Intermountain West limits our understanding of snowpack variability. Fortunately, in the western U.S. tree-ring data have been used to create a broad range of climate reconstructions including precipitation (e.g., [11,12]), temperature (e.g., [13]), and streamflow (e.g., [14][15][16]), as well as snowpack for parts of the Rocky Mountains [17][18][19][20], and could be employed for high resolution reconstructions of SWE.
Tree-ring-based climate reconstructions rely on intensive sampling of trees that have been carefully selected from site(s) and/or species thought to record a climatic variable of interest and can be replicated [21]. Subsequently, statistical reconstruction models rely on the presence and length of tree-ring chronologies across a given region. There have been important efforts to create gridded datasets with broad spatial coverage using tree-ring chronologies, notably for temperature and precipitation [22], as well as the Palmer Drought Severity Index (PDSI) [23]. These gridded products have proven to be useful for the analysis of long-term climate variability at regional to continental scales. These datasets become more limited, however, when it comes to analyzing watershed-scale variability within complex terrain such as that found in the Intermountain West (see Figure 1). The production of multiple tree-ring chronologies that would enable fine-scale and broad-coverage reconstructions is desirable. However, chronologies are extremely time-consuming and costprohibitive to develop, thus alternative data are also highly desirable.  As has been noted in several studies [16,24,25] there are large gaps in the spatial coverage of publicly available tree-ring chronologies (International Tree-Ring Data Bank, ITRDB). Some of the largest gaps are over central and northern Utah, and southeastern Idaho. Because of the previously mentioned difficulties posed by collecting and building tree-ring chronologies everywhere, DeRose et al. [24] explored the possibility of using increment cores collected by the U.S. Forest Service, Forest Inventory and Analysis (FIA) program as potential proxies for climate reconstructions, and as an alternative means to address the paucity of data in the region. The FIA program monitors and characterizes the nation's forest resources through the use of a comparatively high-density sampling grid (see Section 2.2) and a sampling strategy that resulted in the collection of a large number of increment cores on a geographically unbiased grid (~5 km spacing) throughout the Intermountain West, including Utah.
The feasibility of species-specific subsets of the FIA tree-ring dataset in depicting climate variability has been previously demonstrated in DeRose et al. [24]. In general, increment cores from low elevation species (Psuedotsuga menziesii and Pinus edulis) were significantly correlated with water-year (Aug-Jul) precipitation, and in particular over mountainous areas, and were moderately negatively correlated with temperature. Given that accumulated SWE and water-year precipitation are closely correlated [8], it is reasonable to assume that the FIA tree-ring data will serve as a proxy for SWE as well. Therefore, in this paper we develop a reconstruction of April 1 SWE using the FIA tree-ring data combined with any available tree-ring chronologies in the region. We evaluate the reconstructed grid of SWE based on statistical tests, and by comparison with other gridded products. Climate diagnostics are performed to demonstrate the utility and application of the gridded SWE data to water resources management.

Observational and Climate Data
SWE observations from a total of 38 snow courses were obtained from the U.S. Department of Agriculture's Natural Resources Conservation Service (NRCS, http://www.wcc.nrcs.usda.gov/snow/, Figure 1). Locations were selected that had a common period of record covering 1930-1989. All courses selected had no more than 10% of the record without observations (i.e., only six or fewer missing values). Many of these courses have longer records, but some FIA collections in Utah were conducted as early as 1989, so a later date for the calibration period could not be used. Details of snow course data and their collection can be found at http://www.wcc.nrcs.usda.gov/factpub/sect_4a.html. The NRCS automated Snow Telemetry (SNOTEL) data were not used here due to the insufficient temporal coverage (no data prior to the early 1980s). We also note that, while SNOTEL observations are considered by many as the "truth" regarding snowpack conditions, these data have drawbacks; known biases include those caused by wind, aspect, and canopy (e.g., [26]). These drawbacks extend to snow course measurement as well; the standardized sampling device used for NRCS snow course measurements has been shown to overestimate SWE by about 10% [27].
As part of the validation of our reconstructions, we used the North American Drought Atlas (NADA) [23] to compare spatial variability between reconstructed SWE and drought. We also used an unsmoothed tripole index for the Inter-decadal Pacific Oscillation (IPO) created from Hadley Center sea ice and sea surface temperature data [28] to test for any relationship between IPO and Utah SWE.

Tree-Ring Data
At each FIA sampling plot, a suite of measurements are made periodically in order to characterize forest attributes and change over time [29,30]. As part of the Interior West FIA program, increment cores were collected from a subset of trees on each plot to establish plot ages and periodic growth rates. These trees represent the dominant species on the plot in terms of size and forest type. Because increment cores were sampled to represent dominant forest types, many tree species are represented in the FIA tree-ring data set. For this study, we used cores collected from P. menziesii and P. edulis, all of which grow at relatively low elevations and exhibit ring-width variability known to reasonably reflect moisture availability. These species are potentially susceptible to drought and have been shown to be appropriate for use as climate proxies (e.g., [11,31]). The FIA increment cores were subsequently measured and cross-dated, and each time-series was then detrended. In an attempt to remove the effects of stand dynamics and biological growth patterns, we used a cubic smoothing spline with a frequency response of 0.5 at a wavelength of 2/3 times the tree-ring series length [32]. In our region of interest (36.5-43.5 N, 114.5-107.5 W), there was a total of 967 FIA samples. Some of these samples did not provide sufficient temporal coverage (i.e., young trees) for our intended reconstruction period of 1850-1990, and thus were excluded; after our selection process, which included temporal coverage and species, the final number of FIA plots included in our study as potential predictors was 110.
The FIA tree-ring ring-width indices were supplemented by (a) 31 P. menziesii and P. edulis chronologies available from the ITRDB and (b) 31 additional tree-ring chronologies developed by the Wasatch Dendroclimatology Research group, WADR [33]. These chronologies were developed for the species P. menziesii, P. edulis, P. ponderosa, Juniperus scopulorum, and J. osteosperma, and have been used in reconstructions of streamflow for the Logan River [25] and Weber River [16], and of Great Salt Lake (GSL) water level [34]. The WADR chronologies in particular provided added sample depth for northern Utah, where both FIA samples and ITRDB samples were sparse ( Figure 1).

Reconstruction Method
The FIA tree-ring data were evaluated for temporal consistency of their relationship to 1 April SWE using correlations. Low correlations indicated possible instability, and high correlations were necessary to be included in the models. For each snow course, the tree-ring data with the highest correlations (both concurrent with observations and lagged one year) were then used in a multiple linear regression model, of the form: built in a stepwise fashion by adding predictors one at a time, ranging from one predictor to a maximum of 12, objectively determined to maintain model parsimony. We set no spatial control on potential model predictors such that any snow course site could be modeled by any combination of tree-ring data. The predictive potential of each model was tested using the predictive error sum of squares (PRESS; [35]), a leave-one-out validation method taking the form: where n = number of observations, Y i is the i th observation, andŶ i,−i is the predicted value of Y i when Y i was excluded from the regression. For each snow course, the regression model with the smallest PRESS and thus the highest predictive skill was selected as the final model. As metrics for model performance, we report variance explained (R 2 ), as well as adjusted R 2 , calculated as: where n = sample size and k = number of predictors. This test penalized the addition of predictors and thus reduced the inflation of R 2 due simply to additional regression parameters, and thereby provided a stricter test [36]. We also used a cross-validated R 2 based on the PRESS statistic: In the case of poor predictive models, R 2 cv can take negative values if PRESS is greater than total sum of squares (TSS), that is if model predictions are worse than using the observational data.

Gridded Interpolation
After regression modeling, the reconstructed SWE for all snow courses were interpolated to a regular rectangular grid. Such a grid extended the usefulness of data by enabling the integration of SWE estimates over space, and resulted in estimates of total snowpack water content over a region. Previous, gridded SWE datasets have been created with ground observations [37] and satellite measurements [38], but are less common than other gridded climate variables such as temperature and precipitation and have limited temporal coverage (post-1978).
In this study, the combined SWE reconstructions were interpolated to a regular grid using Cressman objective analysis [39,40] to 4-km resolution, similar to the grid spacing used by the PRISM dataset [41]. Cressman analysis used an iterative process of considering station data within successively smaller radii of influence to determine the interpolated value. We used four radii of influence: 3 • , 2 • , 1 • , and 1 2 • for our grids. A data mask was applied to any grid points with an elevation below 1800 m (~6000 ft.), where there is generally not snowpack accumulation. Because Cressman analysis did not account for orography in the gridding process, we adjusted the grid point values using: where E g = elevation of grid point and E re f = average elevation of snow courses used in the reconstruction. This took into account the relationship between SWE observations and elevation in the data set ( Figure 2). Cressman analysis also assumed that the data being interpolated have spatial and temporal stationarity; however, this assumption might not hold true with our data. To address the potential for spatial non-stationarity, we added dummy observations in our region of interest filled with zeroes, essentially zero-padding the data. This was accomplished by adding 250 locations, which were assigned using a random number generator, provided that the locations were at least 75 km from any snow course location, and yielded a density of dummy locations that was similar to the density and distribution of the snow courses stations. The zero-padding approach was intended only to minimize extrapolation problems at the edges of the domain, and was not necessarily an assumption that the point was actually snow-free ( Figure 3).

Gridded Interpolation
After regression modeling, the reconstructed SWE for all snow courses were interpolated to a regular rectangular grid. Such a grid extended the usefulness of data by enabling the integration of SWE estimates over space, and resulted in estimates of total snowpack water content over a region. Previous, gridded SWE datasets have been created with ground observations [37] and satellite measurements [38], but are less common than other gridded climate variables such as temperature and precipitation and have limited temporal coverage (post-1978).
In this study, the combined SWE reconstructions were interpolated to a regular grid using Cressman objective analysis [39,40] to 4-km resolution, similar to the grid spacing used by the PRISM dataset [41]. Cressman analysis used an iterative process of considering station data within successively smaller radii of influence to determine the interpolated value. We used four radii of influence: 3°, 2°, 1°, and ½° for our grids. A data mask was applied to any grid points with an elevation below 1800 m (~6000 ft.), where there is generally not snowpack accumulation. Because Cressman analysis did not account for orography in the gridding process, we adjusted the grid point values using: where = elevation of grid point and = average elevation of snow courses used in the reconstruction. This took into account the relationship between SWE observations and elevation in the data set ( Figure 2). Cressman analysis also assumed that the data being interpolated have spatial and temporal stationarity; however, this assumption might not hold true with our data. To address the potential for spatial non-stationarity, we added dummy observations in our region of interest filled with zeroes, essentially zero-padding the data. This was accomplished by adding 250 locations, which were assigned using a random number generator, provided that the locations were at least 75 km from any snow course location, and yielded a density of dummy locations that was similar to the density and distribution of the snow courses stations. The zero-padding approach was intended only to minimize extrapolation problems at the edges of the domain, and was not necessarily an assumption that the point was actually snow-free ( Figure 3).

Climate Diagnostics
To test the added value of including the FIA data in these reconstructions, we applied our methodology to two sets of data: (1) all tree-ring data, which included the FIA plots, and (2) by including only the ITRDB and WADR tree-ring chronologies. We tested for model improvement in the all tree-ring data model by conducting a modified F-test to compare models with a potentially different number of parameters. We further investigated the spatial and temporal dynamics of the gridded SWE reconstruction by the application of an empirical orthogonal function (EOF) analysis, which indicated the major modes of variability in space and time in a particular dataset. Specifically, we correlated the first and second principal components (PC) derived from the SWE reconstruction and the NADA, and we performed a spectral decompostion of both PCs. Finally, we used the IPO index, low-pass filtered at four years to eliminate high-frequency noise, and performed a lagged correlation analysis on the first and second PCs. We chose IPO because it characterizes pan-Pacific tropical sea surface temperature anomalies [42] that represent both inter-annual events (i.e., El Niño/Southern Oscillation, ENSO) as well as the decadal to multi-decadal variability inherent in such teleconnections as the Pacific Decadal Oscillation or the Quasi-Decadal Oscillation.

Calibration and Model Validation
Reconstructed SWE model skill ( 2 ) varied across the 38 snow coreses used in the state from as low as 0.25 to as high as 0.71, but averaged 0.49 when all tree-ring data were used. Model verification also varied across the state, 2 ranged from 0.23 to 0.67, and averaged 0.46, and 2 ranged from 0.17 to 0.59, and averaged 0.36 (Table S1). Limiting the tree-ring data to ITRDB and WADR data reduced the range of skill to 0.18-0.62, and averaged 0.43, when compared to the full tree-ring data set. Model verification on the reduced tree-ring data set resulted in 2 that ranged from 0.15 to 0.57, and averaged 0.4, and 2 ranged from 0.09 to 0.49, and averaged 0.34 (Table S2). Fully 60% of the snow courses had significantly higher skill when FIA data were included, than with tree-ring chronologies alone (Table S1). Therefore, overall the full tree-ring data set produced the gridded reconstruction with the highest skill, where all but five of the snow courses yielded an R 2 of at least 0.4, and over half exceeded 0.5. Depicted spatially, the gridded results for the individual

Climate Diagnostics
To test the added value of including the FIA data in these reconstructions, we applied our methodology to two sets of data: (1) all tree-ring data, which included the FIA plots, and (2) by including only the ITRDB and WADR tree-ring chronologies. We tested for model improvement in the all tree-ring data model by conducting a modified F-test to compare models with a potentially different number of parameters. We further investigated the spatial and temporal dynamics of the gridded SWE reconstruction by the application of an empirical orthogonal function (EOF) analysis, which indicated the major modes of variability in space and time in a particular dataset. Specifically, we correlated the first and second principal components (PC) derived from the SWE reconstruction and the NADA, and we performed a spectral decompostion of both PCs. Finally, we used the IPO index, low-pass filtered at four years to eliminate high-frequency noise, and performed a lagged correlation analysis on the first and second PCs. We chose IPO because it characterizes pan-Pacific tropical sea surface temperature anomalies [42] that represent both inter-annual events (i.e., El Niño/Southern Oscillation, ENSO) as well as the decadal to multi-decadal variability inherent in such teleconnections as the Pacific Decadal Oscillation or the Quasi-Decadal Oscillation.

Calibration and Model Validation
Reconstructed SWE model skill (R 2 ) varied across the 38 snow coreses used in the state from as low as 0.25 to as high as 0.71, but averaged 0.49 when all tree-ring data were used. Model verification also varied across the state, R 2 adj ranged from 0.23 to 0.67, and averaged 0.46, and R 2 cv ranged from 0.17 to 0.59, and averaged 0.36 (Table S1). Limiting the tree-ring data to ITRDB and WADR data reduced the range of skill to 0.18-0.62, and averaged 0.43, when compared to the full tree-ring data set. Model verification on the reduced tree-ring data set resulted in R 2 adj that ranged from 0.15 to 0.57, and averaged 0.4, and R 2 cv ranged from 0.09 to 0.49, and averaged 0.34 (Table S2). Fully 60% of the snow courses had significantly higher skill when FIA data were included, than with tree-ring chronologies alone (Table S1). Therefore, overall the full tree-ring data set produced the gridded reconstruction with the highest skill, where all but five of the snow courses yielded an R 2 of at least 0.4, and over Water 2017, 9, 403 7 of 13 half exceeded 0.5. Depicted spatially, the gridded results for the individual reconstructions indicated the relationships (R 2 ) between the reconstructed and observed grids. Almost the entire grid scored above 0.4 despite the predictor data set (Figure 4). However, the use of the FIA tree-ring data enhanced the spatial representation of model skill when compared to the reduced data set (Figure 4a versus Figure 4b), where the model skill for almost the entire state exceeded 0.5.
Water 2017, 9, 403 7 of 13 reconstructions indicated the relationships (R 2 ) between the reconstructed and observed grids. Almost the entire grid scored above 0.4 despite the predictor data set ( Figure 4). However, the use of the FIA tree-ring data enhanced the spatial representation of model skill when compared to the reduced data set (Figure 4a versus Figure 4b), where the model skill for almost the entire state exceeded 0.5.

EOF Analysis and Climate Diagnostics
EOF 1 indicated statewide coherent variability in the FIA tree-ring-based SWE reconstruction, and was by far the largest mode of variability (~87%). EOF 2 showed a northwest-southeast dipole, indicated by a much weaker mode (~9%). In comparison, the first leading mode for both the FIAbased reconstruction and the NADA explained 87.2% and 43.7% of the variability, respectively, and were depicted similarly over space (Figure 5a,b). Likewise, the second leading mode for FIA and the NADA explained an additional 8.7% and 21.4% of the variability, respectively, and had reasonably similar spatial representation (Figure 5e,f). While the first two modes of variability in the FIA-based reconstruction explained over 95% of the total variance, the first two modes of NADA explained ~67%. Furthermore, both principle component (1 & 2) time-series for the SWE reconstruction compared to NADA were highly significantly correlated (p < 0.01, Figure 5c,g). Significant periodicity at 2-4 years and ~16-40 years was apparent in PC1 of the SWE reconstruction (Figure 5d). For PC2 there was significant periodicity at 4, 12 and ~50 year frequencies (Figure 5h). Finally, lagged correlation between PC1 the low pass-filtered IPO was significant (p < 0.01) at a lag of 6 years, whereas the relationship for PC2 was significant at a lag of 1 year (Figure 6).

EOF Analysis and Climate Diagnostics
EOF 1 indicated statewide coherent variability in the FIA tree-ring-based SWE reconstruction, and was by far the largest mode of variability (~87%). EOF 2 showed a northwest-southeast dipole, indicated by a much weaker mode (~9%). In comparison, the first leading mode for both the FIA-based reconstruction and the NADA explained 87.2% and 43.7% of the variability, respectively, and were depicted similarly over space (Figure 5a,b). Likewise, the second leading mode for FIA and the NADA explained an additional 8.7% and 21.4% of the variability, respectively, and had reasonably similar spatial representation (Figure 5e,f). While the first two modes of variability in the FIA-based reconstruction explained over 95% of the total variance, the first two modes of NADA explained~67%. Furthermore, both principle component (1 & 2) time-series for the SWE reconstruction compared to NADA were highly significantly correlated (p < 0.01, Figure 5c,g). Significant periodicity at 2-4 years and~16-40 years was apparent in PC1 of the SWE reconstruction (Figure 5d). For PC2 there was significant periodicity at 4, 12 and~50 year frequencies (Figure 5h). Finally, lagged correlation between PC1 the low pass-filtered IPO was significant (p < 0.01) at a lag of 6 years, whereas the relationship for PC2 was significant at a lag of 1 year (Figure 6).

Discussion
The SWE reconstruction indicated particularly good skill in the southern and central Wasatch and central Uinta Mountains, which are crucial areas for snowpack water storage as they feed into the largest and most important tributaries in the region, the Great Basin and the Colorado Basin. The reconstruction skill is consistent with, but shows improvement over, the results in DeRose et al. [24] where water-year precipitation was correlated with FIA tree-ring data and showed the strongest relationships in mountainous regions. In comparison, the gridded reconstruction created without the FIA tree-ring data showed reduced variance explained throughout the entire state. One important similarity between both reconstruction grids is the failure to adequately explain a few snow courses in the extreme northern portion of the state. Increasing predictive skill in this region should be the subject of future inquiry. It should be noted that while the variance explained by these regression models is lower than that reported in water-year or streamflow reconstructions (which can produce

Discussion
The SWE reconstruction indicated particularly good skill in the southern and central Wasatch and central Uinta Mountains, which are crucial areas for snowpack water storage as they feed into the largest and most important tributaries in the region, the Great Basin and the Colorado Basin. The reconstruction skill is consistent with, but shows improvement over, the results in DeRose et al. [24] where water-year precipitation was correlated with FIA tree-ring data and showed the strongest relationships in mountainous regions. In comparison, the gridded reconstruction created without the FIA tree-ring data showed reduced variance explained throughout the entire state. One important similarity between both reconstruction grids is the failure to adequately explain a few snow courses in the extreme northern portion of the state. Increasing predictive skill in this region should be the subject of future inquiry. It should be noted that while the variance explained by these regression models is lower than that reported in water-year or streamflow reconstructions (which can produce

Discussion
The SWE reconstruction indicated particularly good skill in the southern and central Wasatch and central Uinta Mountains, which are crucial areas for snowpack water storage as they feed into the largest and most important tributaries in the region, the Great Basin and the Colorado Basin. The reconstruction skill is consistent with, but shows improvement over, the results in DeRose et al. [24] where water-year precipitation was correlated with FIA tree-ring data and showed the strongest relationships in mountainous regions. In comparison, the gridded reconstruction created without the FIA tree-ring data showed reduced variance explained throughout the entire state. One important similarity between both reconstruction grids is the failure to adequately explain a few snow courses in the extreme northern portion of the state. Increasing predictive skill in this region should be the subject of future inquiry. It should be noted that while the variance explained by these regression models is lower than that reported in water-year or streamflow reconstructions (which can produce R 2 values as high as~0.7), April 1 is an arbitrarily selected time to measure SWE that is nonetheless accepted as a standard for water resource management because in many years this date provides a reasonable estimate of peak snowpack. Ideally, peak SWE is a better estimate of water availability, and therefore potential tree growth, but such data are not available prior to the SNOTEL era (early 1980s).
As in any regression modeling context, the amount of variance explained decreases from R 2 to R 2 adj due to the differing number of parameters. Similarly, our chosen metric of validation, R 2 cv , will be necessarily lower than statistics that describe model fit. The relationship among fit statistics resulted in some interesting patterns that bear further scrutiny. In a number of models, the addition of FIA predictors resulted in better R 2 (e.g., King's Cabin-Upper, Silver Lake-Brighton, Panguitch Lake), but the improvement was not statistically significant. In contrast, there were a number of sites where the addition of FIA data did not noticeably change the model fit (e.g., Mosby Mountain, Garden City Summit, Mammoth-Cottonwood). In most of these situations the stations were located at very high elevations (greater than approximately 2900 m), suggesting that our modeling approach, which was composed of almost entirely low-elevation P. menziesii and P. edulis as predictors in both datasets, might be enhanced by incorporating other less-conventional species [43]. Although these species have been shown to be moisture-sensitive in this region, and therefore some tree growth would be associated with moisture that arrived other than from snowpack, there was~30-50% variance left unexplained within our SWE reconstructions. At least part of this may be due to the indirect linkage between tree growth and previous winter snowpack. Regardless, given the close relationship between accumulated snowpack and water-year precipitation in the region, we are confident in the validity of our reconstruction. Finally, in some situations the no-FIA model actually produced higher skill (e.g., Fish Lake), which was a result of applying our stepwise linear models independently. Although, it is not unreasonable to expect tree-ring chronologies, in the absence of more local FIA data, to produce high skill models.
While our analyses suggested that the gridded SWE data in Utah were broadly consistent with drought variability at the regional scale, it should be noted that some of the predictors (i.e., ITRDB chronologies) were likely used in the creation of both our reconstruction and the NADA, and so some correspondence should be expected. It is also worth noting that the NADA was calibrated for the June, July, August season, while April 1 SWE is the result of winter precipitation, so we suggest that the similarity of the PCs between reconstructions is more than simply a statistical artifact, but rather reflects the effect of broad-scale climatic drivers on regional-scale hydroclimatology (i.e., accumulated snowpack is the primary source of moisture for the region, and therefore will have a strong influence on drought indices).

Comparison of Spatial Variability
The periodicity seen in PC1 is consistent with the prominent quasi-decadal variability recorded by the PDSI in Utah [44] and Great Salt Lake levels [45]. These studies proposed that this unique climate variability is caused by a type of teleconnection excited during the transitional phases of the Pacific quasi-decadal oscillation (QDO), rather than the extreme warm/cold phases that appeared in the central tropical Pacific in previous research on Great Salt Lake levels [46,47]. This connection was also found with the IPO, which occurs on a~30 year frequency; this finding, along with the PDSI analysis of Herweijer et al. [48], lends support to the multidecadal variability found in PC1 of reconstructed SWE. Additionally, although PC1 had a much stronger effect (i.e., explained more variance) than PC2, when both occur simultaneously such as in 1934 (Figure 5c-g), the reduction in SWE can be severe.
The year 1934 was one of the strongest regional drought years in Utah, often used as the drought year-of-record in water management, which coincided with the Dust Bowl era. This result indicates that ENSO forcing can be modulated by the larger-scale and stronger effects linked with PC1 or the Pacific QDO's transition phases, as was the case proposed for the Dust Bowl [42], and record floods in northern Utah [49] and the Missouri River Basin [50]. Indeed, the significant correlation of PC1 to IP at a lag of six years is consistent with the transition-phase teleconnection presented by Wang et al. [47].
Variability in the second mode was consistent with ENSO-driven variability in terms of spatial patterns (i.e., the north/south dipole), the relatively weak loading, and the lower percentage of variance explained. This second mode corroborates previous observations that Utah lies along the fulcrum of the regional-scale, north-south ENSO dipole [44,[51][52][53]. The correlation between PC2 and the IPO at a lag of one year was consistent with ENSO forcing. These modes of variability were strikingly consistent with the first two modes found in the EOF analysis of NADA [23], in which EOF1 exhibited domain-wide coherent behavior with a regional center over Utah. By comparison, EOF2 of the continental-scale PDSI depicted the ENSO dipole with the wet/dry division cutting across Utah. However, the NADA provides only a coarse representation of the location of the dipole fulcrum; EOF2 of gridded SWE depicted much more clearly the spatial representation of this fulcrum, which bisected the northern Wasatch Mountains from the central Utah plateau region, and the Uinta Mountains from the Uinta Basin. The dipole indicated that in some years the response to ENSO within Utah can be opposed in extreme northern or southern domains, and that little correspondence to ENSO can be expected for central portions of Utah.

Conclusions
The gridded, high resolution (4 km), SWE reconstruction represents a useful proxy of past mountain snowpack at the watershed level. Observational data in the Intermountain West are sparse and have limited temporal coverage at primarily low elevations (except the more recent SNOTEL), and therefore are difficult for use in depicting past climate extremes. Our SWE reconstruction for Utah more than doubled the temporal coverage available from instrumental records (i.e., snow courses). Because the understanding of historical controls on snowpack variability provides insight into snowpack conditions, any extended record of SWE like the one presented here can help improve our understanding of hydroclimatic variability, and specifically on water availability and climate change impacts in the region. The results of the analyses here help indicate which areas of the state are likely to respond to relative changes in climatic forcing (e.g., ENSO versus IPO-driven) in the future, which can influence water availability forecasting. In an era of decreasing snowpack, understanding the watershed-level detail of past snowpack will be important for water managers and water management. Furthermore, understanding the potential repercussions of the simultaneous occurrence of the primary and secondary modes of variability in the region will help water managers anticipate and plan for extreme droughts (e.g., 1934) or extreme pluvials (e.g., early 1980s), in particular, given that the primary modes of SWE variability statewide is predominantly driven by decadal to multi-decadal climate forcing. The FIA tree-ring data and the reconstructed SWE grid data are provided to the public at http://cliserv.jql.usu.edu/FIAdata/.

Supplementary Materials:
The following are available online at www.mdpi.com/2073-4441/9/6/403/s1, Table S1: Summary statistics for snow course reconstructions using all tree-ring data. Table S2: Summary statistics for snow course reconstructions using only ITRDB and WADR tree-ring data.