Evaluating multiple historical climate products in ecological models under current and projected temperatures.

Gridded historical climate products (GHCPs) are employed with increasing frequency when modeling ecological phenomena across large scales and predicting ecological responses to projected climate changes. Concurrently, there is an increasing acknowledgement of the need to account for uncertainty when employing climate projections from ensembles of global circulation models (GCMs) and emissions scenarios. Despite the growing usage and documented differences among GHCPs, uncertainty characterization has primarily focused on GCM and emissions scenario choice, while the consequences of using a single GHCP to make predictions over space and time have received less attention. Here we employ average July temperature data from observations and seven GHCPs to model plant canopy cover and tree basal area across central Alaska, USA. We first compare the fit of, and support for, models employing observed temperatures, GHCP temperatures, and GHCP temperatures with an elevation adjustment, finding (1) greater support for, and better fit using, elevation-adjusted vs. raw temperature models and (2) overall similar fits of elevation-adjusted models employing temperatures from observations or GHCPs. Focusing on basal area, we next compare predictions generated by elevation-adjusted models employing GHCP data under current conditions and a warming scenario of current temperatures plus 2°C, finding good agreement among GHCPs though with between-GHCP differences and variation primarily at middle elevations (~1,000 m). These differences were amplified under the warming scenario. Finally, using pooled indices of prediction variation and difference across GHCP models, we identify characteristics of areas most likely to exhibit prediction uncertainty under current and warming conditions. Despite (1) overall good performance of GHCP data relative to observations in models and (2) positive correlation among model predictions, variation in predictions across models, particularly in mid-elevation areas where the position of treeline may be changing, suggests researchers should exercise caution if selecting a single GHCP for use in models. We recommend the use of multiple GHCPs to provide additional uncertainty information beyond standard estimated prediction intervals, particularly when model predictions are employed in conservation planning. The Authors. Ecological Applications published by Wiley Periodicals LLC on behalf of Ecological Society of America.

Introduction

Gridded historical climate products (GHCPs) are composed of one or more continuous surfaces (or fields) of an interpolated meteorological characteristic (e.g., temperature, precipitation) typically averaged over multiple decades and built from weather station data and related inputs. These products are increasingly employed to investigate ecological relationships (Running et al. 1987, Guisan and Thuiller 2005, Abatzoglou 2013) and plan for the ecological consequences of climate change (Carvalho et al. 2011, Littell et al. 2011). These products were developed for a range of applications such as climate trend evaluation (e.g., TopoWx; Oyler et al. 2015), spatial climate analysis (e.g., PRISM; Daly et al. 2008), and to meet ecological research needs (Kittel et al. 2004, Baker et al. 2017, Abatzoglou et al. 2018) such as providing critical input to large‐scale biogeophysical and related process models (Wei et al. 2014, Lawrence et al. 2019). GHCPs continue to improve as observational networks become more representative (Simpson et al. 2005), informed in part by studies in mountainous and remote terrain (Lundquist and Cayan 2007, Dobrowski et al. 2009, Minder et al. 2010, Pike et al. 2013, Sadoti et al. 2018). Recent work capitalizing on this nexus of ecological and climate data, coupled with climate projections from general circulation models (GCMs), has improved confidence in ecological changes anticipated under future greenhouse gas emissions scenarios (Carvalho et al. 2011, Markon et al. 2012, Flato et al. 2013). The literature addressing climatic influences on ecological systems is vast, and our discussion here focuses on new developments driven by the proliferation of relatively high‐resolution GHCPs.

Choices made during data set construction mean that different GHCPs may disagree in their means, trends, and extremes (Hofstra et al. 2010, McAfee 2020). Sources of variation between GHCPs include data inputs (e.g., choice of weather stations or inclusion of reanalysis products; Yin et al. 2015, Behnke et al. 2016), the baseline climatological period (Suggitt et al. 2017, Newman et al. 2019), spatial resolution (Mislan and Wethey 2011, Suggitt et al. 2017), and interpolation algorithm (Hofstra et al. 2008, Newman et al. 2019). Additionally, individual GHCP iterations may include different sets of input data as they are updated, impeding the detection of changes between versions (PRISM Climate Group 2013, Oyler et al. 2015). Although the representation of spatial climate patterns in GHCPs is aided by an increasing density of weather stations at higher elevations and in heterogeneous terrain (McAfee et al. 2019), differences between GHCPs are still likely to be largest where station densities are lower such as high‐elevation areas and in high‐latitude regions (Walton and Hall 2018). Paradoxically, these less accessible areas are those where researchers will have the greatest need for interpolated climate data.

In tandem with a general need for the climate information supplied by GHCPs, ecologists, land managers, and policymakers are increasingly interested in prediction of (1) ecological conditions across large, often unsampled areas under the current (or historical) climate (Elith et al. 2006, Sequeira et al. 2018), and (2) changes in ecological conditions resulting from climate change (Hannah et al. 2002, Araujo et al. 2005, Gogol‐Prokurat 2011). In the pursuit of useable projections of ecological conditions, uncertainty due to variation among GCMs and across emissions scenarios is often assessed (Beaumont et al. 2008, Littell et al. 2011), although additional sources of uncertainty such as modeling methods, the choice of predictors, and ecological data set variation (e.g., density of occurrence records) may be considered (Synes and Osborne 2011, Bagchi et al. 2013). From studies of populations and metapopulations (Ruete et al. 2012, Naujokaitis‐Lewis et al. 2013) to communities (Diniz‐Filho et al. 2009, Cook et al. 2010) and ecosystems (Littell et al. 2011, Sofaer et al. 2017), ecological predictions will vary within ensembles of GCMs and emission scenarios.

In contrast with the numerous assessments of ecological uncertainty under multiple GCMs and emissions scenarios, it is relatively rare for researchers to assess ecological prediction uncertainty associated with ensembles of GHCPs (Baker et al. 2016). The common assumption, either implicit or explicit, is that ecological responses are less sensitive to variation among historical climate data sources than they are to future climate projections (Baker et al. 2016), despite documented and often substantial differences among GHCPs (Simpson et al. 2005, McAfee et al. 2014, Behnke et al. 2016, Henn et al. 2018). In illuminating the effects of variation among GHCPs, recent studies have found that predictions from ecological models may be sensitive to characteristics of GHCPs such as spatial and temporal scale (Suggitt et al. 2017), input baseline period (Roubicek et al. 2010), and overall approaches used in data creation (Bobrowski and Schickhoff 2017, Lembrechts et al. 2019). Likewise, given the interest in the statistical validity of GHCPs (Johns et al. 2003), the potential consequences of not addressing the within‐product uncertainty of these products when modeling ecological responses have received valuable attention (Fuentes et al. 2006, Stoklosa et al. 2015). Despite this growing awareness, the testing and comparison of ecological predictors across GHCPs remains primarily limited to studies of species distributions at coarse scales (Baker et al. 2016). Quantifying the potential consequences of using individual GHCPs for ecological prediction is critical to their appropriate selection at management‐relevant spatial scales.

To address this deficiency, we selected three large conservation units in a rapidly warming region, interior Alaska, USA, and employed seven GHCPs to investigate patterns of plant canopy cover and tree basal area. Our objectives were to assess (1) overall support for employing GHCPs and the variation in ecological responses to July temperature across GHCPs, (2) differences in predicted ecological responses between current conditions and those under a warming scenario, and (3) uncertainty in ecological prediction across the ensemble of GHCPs employed in modeling. In addition, we assessed support for including elevation adjustment predictors in models to account for temperature error due to the ground resolution of GHCP cells. Our goal was to help inform regional ecologists and land managers in Alaska and beyond of important considerations when selecting one or more GHCPs to utilize in a given study. Our investigation builds on interest in anticipating forest changes (Millar et al. 2007) and expands recent work (Lembrechts et al. 2019) to other high‐latitude areas where climate stations are sparse and gridded climate data are often relied upon (Swanson 2015, Roland et al. 2019, Verbyla and Kurkowski 2019).

Methods

Field measurements

We employed measures of vegetation, temperature, topography, and soils (Table 1) described in Roland et al. (2019) from 83 sampling plots across five transects located within three National Park Service units in interior Alaska, USA (Fig. 1). Transects included between 11 and 30 plots (mean 16.6). Surface (~1.5 m) air temperature values were estimated from locally recorded sensor values using the approach described in Roland et al. (2019). For brevity, we hereafter refer to these estimates as observed temperatures.

Table 1

Temperature and environmental covariates used in models of canopy cover presence or tree basal area in interior Alaska, USA.

Covariate	Description	Mean	SD	Minimum	Maximum
Temperature
Observed temperature	mean estimated July air temperature from sensors (Roland et al. 2019; °C)	12.0	1.4	9.3	14.5
Product temperature	mean July air temperature from climate products (°C)	11.7	1.3	8.0	15.1
Microclimate, topography, and soil
Elevation difference	the difference between elevation of climate source grid cells and plot elevation measured via 5‐m resolution IfSAR data for use as a lapse‐rate adjustment (m)	31	146	−838	605
Slope	measured by clinometer in field (°).	21.7	10.2	0.5	36.0
Equivalent latitude	an index of potential solar radiation (Lee 1962)	58.1	16.5	30.1	86.7
Soil depth	average of deepest four field measurements out of 16 using a metal probe (cm)	67.9	30.0	14.8	129.5
Vegetation
Canopy cover	measured in field on two 16‐m point transects using vertical strata (≥1 m above ground; %)	18.3	30.5	0.0	111.9
Tree basal area	sum of basal area for all live trees (measured at 1.37 m above ground) within 8 m of plot centers (m2/ha)	3.3	6.9	0.0	32.1

Covariates were measured in field plots (n = 83), extracted from high‐resolution elevation data, or extracted from gridded historical climate products (GHCPs) across five transects in three National Park Service units.

Fig. 1

Location of central Alaska, USA study area, three National Park Service units, and five plot sampling transects. Also shown are nearby population centers and major roads. The eastern Denali example area was employed in illustrating gridded historical climate products and model predictions spatially.

Climate products

We employed mean July temperatures from seven GHCPs (Table 2). In addition to having spatial coverage in our Alaska study area, each GHCP had a spatial resolution less than 2 km and provided either (1) a 30‐ or 31‐yr climatology or (2) sufficient data to calculate a 30‐yr normal. GHCPs meeting these criteria were Climatologies at High resolution for the Earth’s Land Surface Areas (CHELSA; Karger et al. 2017), Daymet (Thornton et al. 1997), the Parameter‐elevation Regression on Independent Slopes Model (PRISM; Daly et al. 1997, 2008), Scenarios Network for Alaska and Arctic Planning (SNAP; employing the approach of Gray et al. 2014), and WorldClim (Fick and Hijmans 2017). Two GHCPs were available (employing different baseline periods and sets of stations) for PRISM and WorldClim. We extracted temperature and elevation values from gridded data for each plot for subsequent modeling. Elevation grids were available online for each GHCP with the exception of Daymet, which we received directly (M. Thornton, personal communication). Initial inspection of GHCPs revealed offsets over our study area of one‐half cell (northwest) in WorldClim data, one cell (east) in CHELSA data, and one‐half cell (north) in Daymet data over our study area. These offsets were corrected prior to extracting values.

Table 2

Gridded historical climate products employed in modeling vegetation responses in interior Alaska, USA, National Park units.

Product	Climatological baseline(s)	Spatial resolution	Data sources(s) and methods	References
CHELSA	1979–2013†	30‐arc second‡	Elevation‐based interpolation of ERA‐Interim (reanalysis) temperatures were based on linear free‐air lapse rates estimated by the product.	Karger et al. (2017)
Daymet§	1981–2010†	1 km	Local lapse rates are calculated from available Global Historical Climatology Network—Daily stations and applied to a high‐resolution DEM.	Thornton et al. (1997)
PRISM	1971–2000, 1981–2010	30‐arc second‡	Large numbers of stations were used to define local regressions relating temperature to elevation, slope, aspect, and their distance from the coast.	Daly et al. (1997, 2008)
SNAP¶	1961–1990†	1 km	Interpolated temperature anomalies from the CRU TS4.0 data set (Harris and Jones 2017) were added to the 1961–1990 PRISM and then bilinearly interpolated to 1‐km2 resolution.	Gray et al. (2014)
WorldClim	1960–1990 (v1.4), 1970–2000 (v2.0)	30‐arc second‡	Meteorological station data from multiple networks was interpolated with thin‐plate spines accounting for latitude, longitude, and elevation.	Hijmans et al. (2005), Fick and Hijmans (2017)

We calculated the 1981–2010 30‐yr climatology from the available data.

30‐arc second resolution is ~900 m (north‐south) by 400 m (east‐west) at the study area center.

Daymet calculates daily temperatures directly for 1980 to near present, without the use of a climatology, so we calculated a 1981–2010 climatology from the available data.

SNAP used the 1961–1990 PRISM climatology as a base for interpolation, but they serve monthly data for 1901–2015. The SNAP data were constructed by adding interpolated temperature anomalies from the CRU TS4.0 data set (Harris and Jones 2017) to the 1961–1990 PRISM climatology (http://ckan.snap.uaf.edu/dataset/historical‐monthly‐and‐derived‐temperature‐products‐downscaled‐from‐cru‐ts‐data‐via‐the‐delta‐m).

Statistical analysis

Using models identified in Roland et al. (2019) as a baseline, we modeled the relationships of canopy cover (percent cover of vegetation ≥1 m tall) and tree basal area (m2/ha of all live stems ≥1.37 m tall) with observed and GHCP‐derived temperatures using generalized linear mixed models in the glmmTMB package (Brooks et al. 2017) in R (R Development Core Team 2019). Each model included random effects (intercepts) of transect nested within park to account for potential non‐independence of error. We controlled for slope angle and soil depth in canopy cover models and an index of incident solar radiation, slope angle, and soil depth in basal area models (Roland et al. 2019). To address model convergence issues when using GHCP‐derived temperature covariates, we (1) modified the tree basal area model used in Roland et al. (2019) by removing soil temperature as a predictor and (2) modified canopy cover models by modeling the presence vs. absence of canopy cover (i.e., 0% vs. >0%) rather than the continuous range. Additionally, we found more realistic predictions of basal area were achieved by including logit‐transformed scaled (mean = 0, SD = 1) rather than raw (scaled) temperature values in models.

We built two models using temperature data from each GHCP using the set of control variables and either (1) temperature alone or (2) temperature plus a covariate representing an elevation adjustment. The elevation adjustment covariate was the difference between the plot elevation and the GHCP grid cell elevation. A hypothetical linear relationship between an ecological response (y) in location i and temperature values from a GHCP cell (Tgrid) would bewhere β0 is the estimated intercept, β1 is the estimated effect of GHCP‐derived temperature, and ε is random error. In mountainous areas where sub‐cell temperature varies with elevation, it is common to adjust for lapse rate. For example, the approach of Roberts et al. (2019) employed a local lapse rate L (°C/km) and the difference between Eobs and Egrid:

where L is calculated from data in a local network of weather stations, Eobs is a local, higher‐accuracy measure of elevation (e.g., 10‐m digital elevation model), and Egrid is the elevation associated with the GHCP (e.g., average elevation across a 1‐km cell).

Our approach does not assume the availability of local lapse rate information and instead allows the estimation of an indirect effect of lapse rate on y via the addition of an estimated parameter (β2)

In our approach, the elevation adjustment covariate (Eobs – Egrid) is specific to each GHCP, and Eobs is the plot elevation extracted from 5‐m IfSAR data (data available online).2 An illustration of this approach is shown in Fig. 2. If Eobs – Egrid is positive, the actual elevation is higher than the surrounding grid cell elevation. As both canopy cover and basal area tend to decrease with increasing elevation in our study area we would expect β2 to be negative.

Fig. 2

Illustration of how gridded climate data were extracted and processed for inclusion in models predicting canopy cover and basal area in central Alaska, USA. Area shown is the Nabesna transect (Wrangell‐St. Elias National Park and Preserve) and circles indicate the location of sampling plots. Example elevation (elev.) data (A and B), average July temperature (temp., T avg; °C) and an elevation adjustment covariate (D) are shown (diff., difference). High‐resolution (res.) elevation data (A) is 5 m IfSAR data. Thirty arc second elevation (B) and temperature (C) data are from PRISM (1981–2010). The elevation adjustment covariate (D) was calculated from the difference between IfSAR (A) and PRISM (B) elevation.

We identified the relative support among models for both canopy cover presence and basal area using the Bayesian information criterion (BIC; Schwarz 1978) and calculated both (1) marginal (population‐level, no random effects) and (2) k‐fold cross‐validated measures of model fit or accuracy for each model. We calculated pseudo‐R 2 values and root mean‐square error (RMSE) as measures of fit among basal area models. Pseudo‐R 2 was calculated from the ordinary least squares relationship between log‐transformed (+1) observed basal area and model‐predicted values. We calculated the area under the receiver operating curve (AUC, via the AUC package; Ballings and Van den Poel 2013) and the Brier score (mean‐square error) in assessing canopy cover models. K‐fold cross‐validation entailed iteratively withholding each transect of sampling plots, refitting models, and predicting canopy cover presence or basal area values at sampling plots on the withheld transect. K‐fold validation on spatially aggregated samples can also be considered an assessment of model transferability (Wenger and Olden 2012). We assessed variance inflation factors (VIF) of fitted models to ensure multicollinearity was not an issue, finding all models had VIF < 4. In addition, we assessed the independence of model residuals via Moran’s I correlograms (Legendre and Legendre 1998) with bin increments of 200 m, finding no significant correlation (α = 0.05; progressive Bonferroni adjustment) among 1,000 Monte Carlo permutations per bin.

Inter‐product comparisons and variation in model predictions

To assess how models employing different GHCPs may vary in their predictions under current and future climates, we considered one of our ecological responses (tree basal area) and predictions from both current (or near contemporary) July temperatures and a scenario in which July temperatures increased by 2°C above current conditions. This 2°C increase is on the conservative end of 2070–2099 interior Alaska climate projections using RCP 4.5 (Markon et al. 2012). By employing a fixed increase in temperature, we focus the analysis on the consequences of GHCP choice, limiting the uncertainty associated with projections (see Beaumont et al. 2008, Littell et al. 2011, Sofaer et al. 2017). Using models with the elevation adjustment covariate, we assessed inter‐pair similarities across predictions from the seven models employing GHCP data. We considered predictions from the 83 sampling plots and, to illustrate patterns spatially, an example area of 45 × 45 km (2,025 km2) encompassing eastern Denali National Park and Preserve that contains the Mt. Healy and Igloo Creek transects (Fig. 1). We generated predictions for this area at a 100‐m resolution (n = 202,500 cells) standardized across all GHCPs. We calculated three statistics based on predictions of basal area from pairs of models employing different GHCPs (n = 21 pairs): (1) Spearman correlation coefficients (r s), (2) mean absolute differences (MAD), and (3) percentage of plots or cells for which 95% prediction intervals did not overlap.

In addition to comparisons of predicted basal area between paired models, we calculated two aggregate uncertainty statistics from predictions at sampling plots or cells under current and warming‐scenario July temperatures: (1) standard deviation of predictions (log‐transformed values + 1) and (2) the number of inter‐pair model comparisons in which 95% prediction intervals did not overlap. Standard deviation or coefficient of variation are commonly used measures of uncertainty within GCM and emission scenario ensembles (Fuentes et al. 2006), while prediction intervals provided a more robust measure of difference accounting for sometimes large prediction errors. We did not weight models when calculating either statistic (e.g., using BIC model weight) but this may be warranted in other studies.

We assessed spatial patterns of these aggregate uncertainty statistics qualitatively among cells in the eastern Denali example area and, for sampling plots, employed these statistics in generalized linear mixed models to help illuminate their relationships with biophysical characteristics previously considered (elevation, solar radiation, and slope) while accounting for correlated error (e.g., within transects). All models employed random effects as in canopy cover and basal area models and Poisson error. We drew inference from the most supported (lowest BIC) model for each response.

Results

Observed and GHCP temperature

Mean July temperatures at sampling plots differed between GHCPs by an average 0.44°C (range = 0.04–1.00; Fig. 3). Observed temperatures were on average 0.28°C higher than GHCP‐derived temperatures (range = −0.21–0.79). July temperatures from GHCPs were moderately to highly correlated (r s range = 0.41–0.93) and all exhibited strong agreement with observed temperature (r s ≥ 0.69), with the exception of SNAP (r s = 0.31).

Fig. 3

Mean July temperatures (°C; T avg) from observations and gridded climate data products in interior Alaska, USA, National Park Service units. The upper panel shows temperatures from gridded historical climate products (GHCPs) in an example area of eastern Denali National Park and Preserve while the above‐diagonal area of the lower panel shows differences between pairs of GHCPs. The diagonal and below‐diagonal areas of the lower panel show the distribution, correlation (r s), and differences (Δ°C) between values measured at sampling plot locations in which vegetation and other characteristics were measured. Below‐diagonal axes are in °C. Elevation‐adjusted temperature (top panel) is presented to illustrate the corrective effects of employing information from climate products and a 5‐m resolution elevation data set (Fig. 2), though this differs from the implementation in models. Cells are aggregated to a resolution of 500 m for improved clarity.

Vegetation models

Standardized effects of temperature varied across models employing observed and GHCP‐derived temperature data with parameter estimates in the range of 1.0–3.0 in canopy cover models, 0.5–1.2 in conditional basal area model components, and −6.4 to −1.1 in zero‐inflation basal area model components (Fig. 4). Parameter 95% confidence intervals did not bound zero in canopy cover models but did bound zero in basal area models employing 1960–1990 WorldClim. Despite these differences, most 95% confidence intervals associated with gridded and observed temperatures overlapped in both canopy cover and basal area models. The only exception was in the conditional component of basal area models; in this model set, the effect of 1981–2010 PRISM temperature was greater than that of temperature from several other GHCPs.

Fig. 4

Standardized effects (means and 95% confidence intervals) of average July temperature in generalized linear mixed models of canopy cover presence and tree basal area in interior Alaska, USA, National Park Service units. We compared July 2016 temperature (observed; estimated employing temperature sensors) to seven gridded historical climate products in models with and without an elevation adjustment (EA). Each covariate was scaled (mean = 0, SD = 1). Basal area models (zero‐inflated negative binomial error) included both conditional and zero‐inflated temperature effects.

The classification accuracy of canopy cover models was good (AUC ≥ 0.8) to excellent (AUC ≥ 0.9) and pseudo‐R 2 of basal area models was generally over 0.5 (Table 3). K‐fold cross‐validated AUC values of canopy cover models were moderately lower than internal AUC values (mean change = −0.04, range = −0.02 to −0.08), indicating adequate transferability between transects. While AUC and pseudo‐R 2 are not directly comparable statistics, differences between internal and cross‐validated pseudo‐R 2 of basal area models were more substantial (mean change = −0.27, range = −0.10 to −0.40), indicating poorer transferability between transects.

Table 3

Comparison of models employing mean July temperature estimated from sensors (observed) or from gridded historical climate products (GHCPs) in predicting the presence of canopy and tree basal area on 83 plots across five transects in three National Park Service units in interior Alaska, USA.

Product and model	Canopy cover						Tree basal area
	k	ΔBIC	AUC		Brier score		k	ΔBIC	Pseudo‐R 2		RMSE a
			Int.	CV	Int.	CV			Int.	CV	Int.	CV
Null (no temperature) a	5	40.1	0.77	0.72	0.19	0.22	7	66.3	0.29	0.19	6.00	6.31
Observed (2016)	6	13.3	0.90	0.85	0.11	0.17	10	12.6	0.69	0.50	4.70	5.64
CHELSA (1981–2010)
Raw a	6	12.1	0.91	0.89	0.12	0.13	8	30.2	0.60	0.45	4.95	10.01
EA	7	3.3	0.94	0.92	0.10	0.11	11	16.3	0.67	0.44	4.68	9.92
Daymet (1981–2010)
Raw	6	24.0	0.87	0.85	0.15	0.16	10	42.5	0.53	0.20	4.76	8.78
EA	7	3.6	0.94	0.92	0.10	0.12	11	4.7	0.74	0.56	3.72	6.42
PRISM (1971–2000)
Raw	6	15.2	0.89	0.86	0.13	0.15	10	15.6	0.74	0.51	4.71	7.63
EA	7	9.6	0.93	0.90	0.11	0.13	11	0.0	0.77	0.56	4.10	6.99
PRISM (1981–2010)
Raw a	6	9.1	0.92	0.90	0.11	0.13	8	20.2	0.64	0.49	6.94	16.20
EA	7	0.0	0.95	0.92	0.09	0.11	11	21.1	0.68	0.46	6.26	17.36
SNAP (1981–2010)
Raw	6	37.0	0.81	0.79	0.18	0.19	10	61.6	0.42	0.30	5.81	12.75
EA	7	10.4	0.93	0.90	0.11	0.13	11	23.4	0.64	0.49	4.69	11.84
WorldClim (1960–1990)
Raw	6	27.2	0.86	0.80	0.15	0.18	10	43.3	0.56	0.25	5.08	7.09
EA	7	17.3	0.90	0.84	0.11	0.17	11	22.3	0.60	0.37	5.03	8.25
WorldClim (1970–2000)
Raw	6	29.6	0.85	0.79	0.15	0.19	10	42.5	0.55	0.15	5.16	9.37
EA	7	20.7	0.89	0.81	0.12	0.19	11	20.5	0.64	0.35	4.37	7.25

Two models (raw and elevation‐adjusted [EA] temperature) were employed for each GHCP. EA models included an elevation difference covariate (Table 1, Fig. 2). The number of estimated parameters (k) and relative Bayesian information criteria (ΔBIC) is indicated for each model. Area under the operating receiver curve (AUC; a measure of classification accuracy) and Brier score (a measure of predicted probability error) are indicated for canopy cover presence models. Pseudo‐R 2 (calculated from ordinary least squares) and root mean square error (RMSE) are indicated for basal area models. Both internal (int.) and k‐fold cross‐validated (CV; where each transect was a fold) measures are provided. The most supported model in each set is shown in boldface type.

To ensure convergence, these models did not include zero‐inflation parameters.

The inclusion of an elevation adjustment covariate reduced BIC in all canopy cover models (mean = 12.9, range = 6.3–26.8; Table 3) and basal area models (mean = 21.5, range = 13.8–39.4). Classification accuracy (AUC) improved in canopy cover models (mean increase = 0.05, range = 0.03–0.12) with the inclusion of this covariate, while fit (pseudo‐R 2) improved in basal area models (mean increase = 0.11, range = 0.02–0.23). The largest model improvements (reduced BIC, increased AUC or pseudo‐R 2) were observed in models employing SNAP data. A pattern of larger temperature effects in canopy cover models was observed when adjusting for elevation (mean increase = 0.66, range = 0.41–1.01), but a pattern was less evident among basal area models. Model improvements can be attributed primarily to an effective increase in the spatial resolution of predictions due to elevation adjustment covariates (Fig. 3).

With the exception of models employing CHELSA or 1981–2010 PRISM data, models of canopy cover employing raw (no elevation adjustment) GHCP data were less supported than the model employing observed temperature. All basal area models employing raw GHCP data were less supported than the model employing observed temperature. After inclusion of elevation adjustment covariates, five of seven (71%) canopy cover models employing GHCP data and two of seven (29%) basal area models were more supported than models employing observed temperature. The most supported model for both vegetation responses employed a PRISM GHCP.

Inter‐product comparisons of basal area predictions

Under current (1960–2010) climate conditions, predictions of tree basal area at sampling plots from models employing GHCPs were strongly correlated with observed temperatures (mean r s = 0.87, range = 0.84–0.90; Table 4) and with predictions based on other GHCPs (mean r s = 0.87, range = 0.80–0.95; Table 4). Mean absolute differences (MAD) between model predictions were overall small in these comparisons (mean = 1.7 m2/ha, range = 0.9–2.8; Table 4). The mean percentage of plots with non‐overlapping 95% prediction intervals was 3.2% (range = 0.0–8.4%; Table 4, Fig. 5). Within the eastern Denali example area, correlation among model predictions was high (mean r s = 0.97, range = 0.92–0.97; Table 4), MAD between model predictions was higher than in plots (mean = 7.7 m2/ha, range = 3.7–17.9; Table 4), and the percentage of area with non‐overlapping 95% prediction intervals was 1.9% (range = 0–18.3%; Table 4, Fig. 6).

Table 4

Comparisons of predicted tree basal area (m2/ha) from pairs of models employing mean July temperatures from observations (2016) and seven different gridded historical climate products (see Table 2) in interior Alaska, USA National Park Service units.

Product	Sampling plots (n = 83)								Eastern Denali example area
	Obs.	CH	DM	PR1	PR2	SN	WC1	WC2	CH	DM	PR1	PR2	SN	WC1	WC2
Spearman rank correlation (r s)
Observed (2016)		0.93	0.93	0.94	0.77	0.91	0.90	0.90
CHELSA (1981–2010; CH)	0.87		0.95	0.92	0.85	0.95	0.91	0.91		0.94	0.97	0.97	0.84	0.97	0.96
Daymet (1981–2010; DM)	0.89	0.91		0.93	0.71	0.93	0.96	0.95	0.99		0.92	0.91	0.76	0.93	0.93
PRISM (1971–2000; PR1)	0.84	0.85	0.85		0.80	0.90	0.92	0.90	0.99	0.99		0.99	0.79	0.98	0.97
PRISM (1981–2010; PR2)	0.84	0.95	0.89	0.86		0.81	0.66	0.66	0.98	0.97	0.98		0.78	0.97	0.97
SNAP (1981–2010; SN)	0.87	0.86	0.88	0.85	0.87		0.89	0.88	0.94	0.94	0.94	0.92		0.80	0.75
WorldClim (1960–1990; WC1)	0.90	0.89	0.89	0.89	0.82	0.84		0.99	0.99	0.99	0.99	0.97	0.93		0.99
WorldClim (1970–2000; WC2)	0.88	0.87	0.83	0.83	0.82	0.80	0.95		0.97	0.98	0.98	0.96	0.92	0.98
Mean absolute difference (m2/ha)
Observed (2016)		5.9	4.8	5.6	5.8	6.4	9.4	6.7
CHELSA (1981–2010)	2.1		5.5	3.7	4.9	4.2	7.3	5.1		8.3	8.0	9.2	17.2	7.6	29.3
Daymet (1981–2010)	2.0	1.2		5.5	4.9	5.3	9.4	6.6	2.7		13.7	8.1	16.2	12.4	33.2
PRISM (1971–2000)	1.5	1.7	1.7		4.7	4.3	7.0	4.5	3.6	5.6		11.9	20.8	7.9	27.8
PRISM (1981–2010)	2.8	1.2	1.6	2.0		4.2	10.4	7.6	5.6	4.5	7.3		13.2	14.4	37.2
SNAP (1981–2010)	2.3	1.2	1.5	1.5	1.3		8.1	6.0	6.8	6.1	8.5	4.2		23.0	43.5
WorldClim (1960–1990)	1.5	1.6	1.4	1.1	2.2	1.9		3.8	2.8	3.1	4.5	4.0	6.0		23.2
WorldClim (1970–2000)	1.1	1.6	1.5	1.2	2.4	2.0	0.9		12.3	13.4	11.2	17.2	17.9	14.4
Significant difference (%)
Observed (2016)		1.2	0.0	1.2	4.8	4.8	6.0	1.2
CHELSA (1981–2010)	1.2		1.2	1.2	1.2	1.2	9.6	4.8		0.0	5.2	0.0	2.2	0.0	1.1
Daymet (1981–2010)	6.0	1.2		0.0	10.8	6.0	2.4	2.4	0.0		11.2	0.1	8.7	0.1	1.2
PRISM (1971–2000)	4.8	1.2	6.0		6.0	4.8	3.6	4.8	0.0	0.5		1.9	3.1	17.8	26.8
PRISM (1981–2010)	4.8	0.0	2.4	2.4		2.4	20.5	15.7	0.0	0.0	1.7		5.3	0.4	9.0
SNAP (1981–2010)	2.4	1.2	3.6	2.4	0.0		25.3	14.5	0.1	0.5	1.0	1.6		12.0	23.6
WorldClim (1960–1990)	4.8	0.0	4.8	1.2	6.0	2.4		0.0	0.0	0.0	0.0	0.1	0.4		0.1
WorldClim (1970–2000)	8.4	0.0	6.0	8.4	4.8	2.4	0.0		2.0	0.2	18.3	1.5	9.2	3.3

Shown are statistics for sampling plots (n = 83) and a 45 × 45 km area in eastern Denali National Park and Preserve (see Fig. 1). Matrices of 8 × 8 or 7 × 7 correspond to panels in Figs. 5 and 6, respectively. Values below the diagonal in each matrix are for predictions using current mean July temperatures, while those above the diagonal are for a warming scenario of current temperatures plus 2°C. The percent significant statistic is the percentage of plots or area in which 95% prediction intervals did not overlap.

Fig. 5

Comparisons of predicted tree basal area (m2/ha) in plots sampled in interior Alaska, USA, National Park Service units. Predictions are from models built using observed average July temperatures (x‐axes) and temperatures extracted from seven gridded historical climate products. Predictions are presented both for models employing current temperatures (below diagonal) and a warming scenario of current temperatures plus 2°C (above diagonal). Symbols indicated park units. Black outlined symbols indicate 95% confidence intervals of observed and gridded product temperature model predictions did not overlap. The difference in predictions is indicated by color. Descriptive statistics associated with these comparisons are presented in Table 4.

Fig. 6

Predictions and prediction differences of basal area (m2/ha) using seven gridded climate data products across a 45 × 45 km area in eastern Denali National Park and Preserve. Predictions are from models employing July average air temperature, an elevation adjustment covariate, an index of solar radiation, slope, and soil depth (values were held at 67.9 cm). Panels below the diagonal illustrate differences in predictions between models employing different climate products (row minus column). Panels above the diagonal illustrate differences in predictions between models employing different climate products under a warming scenario (product plus 2°C; column minus row). Areas in which 95% prediction intervals did not overlap between models are outlined in black. Statistics associated with these between‐product comparisons are presented in Table 4. Black areas are masked as they are situated below the range of elevation values employed in models. Cells are aggregated to a resolution of 500 m for improved clarity.

Predictions at sampling plots were similarly correlated among those generated from models using observed temperature and GHCP data under the warming scenario (mean r s = 0.88, range = 0.66–0.99; Table 4). MAD between model predictions at sampling plots was larger (mean = 6.0 m2/ha, range = 3.7–10.4; Table 4), as was the percentage of plots with non‐overlapping 95% prediction intervals (mean = 5.6%, range = 0.0–25.3%; Table 4, Fig. 5). Within the eastern Denali example area, correlation among model predictions remained overall high (mean r s = 0.91, range = 0.75–0.99; Table 4), MAD between model predictions more than doubled over current conditions (mean = 18.4 m2/ha, range = 7.6–43.5; Table 4), and the average percent area with non‐overlapping 95% prediction intervals increased threefold to 6.2% (range = 0.0–26.8%; Table 4, Fig. 6).

Basal area prediction uncertainty

Quantifying the uncertainty among predictions of basal area from models employing different GHCPs in eastern Denali (Fig. 7) revealed patterns of lowest variation and fewer inter‐pair model differences at higher elevations, moderate variation and few inter‐pair model differences at low elevations, and greatest variation and the most inter‐pair model differences at middle elevations. Variation and the number of inter‐pair model differences between predictions were overall greater under the warming scenario than current conditions, with largest positive differences in middle‐elevation areas.

Fig. 7

Variation and inter‐pair differences in predicted tree basal area (m2/ha) in the eastern Denali National Park and Preserve example area among models employing mean July temperatures from seven gridded historical climate products. Predictions were made using current climate conditions and a warming scenario of current temperatures plus 2°C. Variation between predictions was calculated as 10 times the log standard deviation (plus 1) across seven models while the measure of pairwise differences was the total number of climate product model combinations (out of 21 possible) in which 95% prediction intervals did not overlap. Each panel is an area of 45 × 45 km with cells aggregated to a resolution or 500 m for improved clarity. Black areas are masked as they are situated below the range of elevation values employed in models. See Fig. 1 for location of the example area.

The assessment of basal area prediction uncertainty at sampling plots under current climate conditions revealed similar associations with biophysical characteristics for variation and inter‐pair differences (Table 5). Prediction variation was lower and inter‐pair differences were fewer in sampling plots on steeper slopes and with higher equivalent latitudes (i.e., cooler, more north‐facing slopes). Variation and the probability of at least one inter‐pair difference exhibited a quadratic relationship with elevation, with both peaking at ~1,000 m. This elevation is approximately the level of treeline in our study area (Roland et al. 2013). Variation and the number of inter‐pair differences in predicted basal area declined with equivalent latitude under the warming scenario but only the probability of at least one inter‐pair difference increased with elevation. Model fits (pseudo‐R 2 and RMSE) varied across models and most exhibited poor transferability (lower k‐fold cross‐validated pseudo‐R 2, higher RMSE; Table 5).

Table 5

The most supported (lowest Bayesian information criterion [BIC]) models of variation and inter‐pair differences in predictions of tree basal area at sampling plots in interior Alaska, USA National Park Service units.

Response and climate conditions	Conditional	Zero‐inflation	k	ΔBIC	Pseudo‐R 2		RMSE
					Internal	CV	Internal	CV
Variation
Current	elevation + elevation2 + slope + EQ	–	7	53.1	0.63	0.54	1.71	1.84
Warming scenario	EQ	–	4	23.2	0.14	0.03	2.53	2.86
Pairwise differences
Current	slope + EQ	elevation + elevation2	8	6.7	0.20	0.06	1.05	1.83
Warming scenario	EQ	elevation	6	24.9	0.31	<0.01	2.09	3.65

Predictions are from models employing seven different climate products (see Table 2) and two climate conditions (current vs. warming [current plus 2°C]). Variation between predictions was calculated as ten times the log standard deviation (plus 1) across seven models while the measure of inter‐pair differences was the total number of climate product model combinations (out of 21 possible) in which 95% prediction intervals did not overlap. Information is presented for each model including covariates employed in the conditional (Poisson) and zero‐inflation model components, the number of estimated parameters (k), the difference between the BIC of a model employing only random effects and the given model, and the pseudo‐R 2 and RMSE of each model both internally and using a k‐fold cross validation (CV) in which each transect was withheld. Elevation, slope, and equivalent latitude (EQ; a measure of incident solar radiation) were extracted or calculated from 5‐m IfSAR elevation data. The most supported model for each response was determined from all combinations of the three covariates.

Discussion

With few exceptions, models of canopy cover and basal area confirm that similar temperature influences are detected using observations or GHCPs. In models employing GHCP data, those with an elevation adjustment received more support (lower BIC) and had improved model fit and accuracy (higher AUC, lower Brier score, higher pseudo‐R 2, and lower RMSE). This suggests elevation adjustment via our approach or others may be applied effectively across a broad range of analyses in mountainous terrain employing gridded climate data when both the elevation data used in GHCP development and higher‐resolution elevation data are available. Despite the overall similarity of the models employing different GHCPs, there were variations in model estimates and fits that we cannot fully explain but that are likely driven in part by the sets of stations and interpolation algorithms employed in developing these climate products.

While temperatures differed between GHCPs, there were overall strong correlations between all pairs. Given these correlations, it is not surprising that canopy cover model accuracy and basal area model fit were similar across GHCPs. Among models without an elevation adjustment, the canopy cover models employing CHELSA and 1981–2010 PRISM temperature had more support than the model employing 2016 observed temperature. While both models had excellent accuracy (AUC ≥ 0.9), this suggests that data from GHCPs may be more predictive of some ecological responses than higher‐precision short‐term data even though observations effectively capture the spatial variation in local microclimates (Roland et al. 2019). At present, monthly PRISM time series data are not available for Alaska to allow a direct comparison for July 2016. Because some ecological responses, such as tree basal area, are not reflective of short‐term weather conditions, but integrate response to climate conditions over longer periods, it is reasonable that gridded data representing a 30‐yr climatology may be more predictive of ecological responses than observations from a single growing season, which may not be representative of average conditions (Lembrechts et al. 2019).

Despite the growing availability of gridded climate data, ecologists typically employ a single GHCP when modeling an ecological response. Multiple GHCPs, in concert with elevation adjustment and other controlling covariates such as soil depth, provided overall similar predictions of our vegetation responses of interest. However, our study illustrates how the choice to employ a single GHCP may result in biased or lower‐precision effects or predictions relative to other available GHCPs. Further, the unknown causes of differences in model fit and precision between GHCPs lend further support to employing multiple GHCPs to better capture the potential range and uncertainty in model predictions (Lembrechts et al. 2019). This is particularly critical in mountainous terrain due to errors in lapse rate assumptions, heterogeneity in radiative effects, and other factors (Strachan and Daly 2017).

We expected the addition of elevation adjustment covariates to improve model fit or accuracy over unadjusted models as they described elevation variation at a scale smaller than the climate cell (~1 km). Models employing SNAP data had greatest improvement with elevation adjustment. While SNAP (based on CRU TS 4.0) is available at 1‐km2 resolution downscaled from 4‐km2, it exhibits overall low variation across cells within 4‐km2 blocks (data available online).3 However, even in GHCPs such as Daymet with clear temperature variation among neighboring ~1‐km2 cells, there were notable increases in accuracy of canopy cover predictions and fit of basal area predictions following elevation adjustment. As models based on other GHCPs (e.g., PRISM, WorldClim) showed more modest accuracy gains with elevation adjustment, specific reasons for large improvements in the model based on Daymet are not immediately clear, but may be due to the density of input data, interpolation methods, and decisions about which physiographic controls on temperature to account for during GHCP creation (Daly et al. 2008).

In comparing predictions among basal area models employing different GHCPs, several compelling patterns emerged in both our sampling plots and our example area in eastern Denali National Park and Preserve. Both (1) the significant differences in predictions (non‐overlapping 95% prediction intervals) and (2) between‐model variation in average predictions made across multiple GHCPs under the range of current climate (1960–2010) were most prominent in sampling plots or grid cells with intermediate elevations (~1,000 m). In some instances, these differences in predicted basal area exceeded 50 m2/ha or were associated with radiative effects. The greater variation among model predictions at the approximate elevation where much research is focused, treeline, underscores the need for caution when interpreting model results employing one GHCP. The magnitude of differences, the overall area of significant differences, and between‐model variation in predictions tended to be larger under the warming scenario. These differences in predicted basal area between models employing different GHCPs under current conditions vs. the warming scenario highlight the complex behavior of multivariate ecological models and resulting variation in the estimated effects of temperature and other biophysical characteristics (e.g., solar radiation).

Most comparisons of potential ecological change have focused on the uncertainties generated when employing multiple general circulation models and emissions scenarios (Baker et al. 2016, 2017). Likewise, the development of ecological projections tends to focus on understanding how the characteristics of simulated climate interact with ecosystem model sensitivities (Sofaer et al. 2017). While less commonly considered, the variation induced by employing multiple GHCPs during modeling may have conservation consequences when a model response such as species occurrence is sensitive to baseline period or other characteristics of a historical climate data set (Roubicek et al. 2010, Lembrechts et al. 2019). In our study, the choice of GHCP and resulting variation in predicted canopy cover or basal area may be important in addressing the conservation and management of tundra ecosystems either (1) over space under current temperatures or (2) projected for a modest warming scenario (Myers‐Smith et al. 2011). For example, with increasing temperatures, understanding the variation in predicted basal area increases and expansion of forests into tundra may assist in management planning for forest and shrub–tundra species assemblages (Mizel et al. 2016) or tundra obligates such as the white‐tailed ptarmigan (Lagopus leucura; Jackson et al. 2015).

By examining temperature‐vegetation relationships using a suite of gridded historical climate products and accounting for local elevation effects on temperature, we bolstered confidence in predicting vegetation attributes across a large, sparsely instrumented area of interior Alaska under both current climatic conditions and a warming scenario. However, despite correlation among GHCPs and overall agreement among model predictions, the disagreement at the elevation of regional treeline ecotones, and to a greater degree in a scenario of modest warming, highlights the potential pitfalls of employing any single GHCP in modeling. As interest in understanding ecological responses to climate will continue to dictate a need for GHCPs, more holistic approaches in future investigations will better inform best practices. Regardless of application and location, given the increasing availability of (1) multiple GHCPs and (2) elevation data at GHCP‐associated and higher resolutions, we recommend the more widespread adoption of approaches such as ours and, when appropriate, the incorporation of uncertainty resulting both from the creation of GHCPs and from the use of GCM ensembles and multiple emissions scenarios. These practices will improve inference from, and confidence in, ecological models employing historical climate data, and thus better inform management and policy‐making.