To project how big sagebrush populations in British Columbia may migrate to cope with climate warming, ArcGIS’s Habitat Suitability Model tool was utilized in conjunction with climate model data. Data on plant occurrences were collected from GBIF (GBIF, 2021). Elevation data used in the analysis were collected from the Government of British Columbia (specifically, GeoBC), as was data used for map aesthetics. Current climate data and future climate projection data were collected from AdaptWest. Five variables that are assumed to relate to habitat suitability were examined. These were elevation, slope, aspect, temperature, and precipitation. Using these variables, a suitability model was developed from the current distribution of big sagebrush in the province. The model was then altered by inputting future temperature and precipitation projections to determine how suitability will change throughout the province by the year 2080 under two climate scenarios: low-warming conditions and high-warming conditions.

First, the data listed in Table 1 were imported and cleaned within ArcGIS (version 2.8.3). GBIF occurrences outside of B.C. were deleted. This is a justifiable step for this analysis, as the occurrences within B.C. represent the plant’s most northern populations. It is anticipated by some researchers that populations at the northern range edge of a species’ distribution are most evolutionarily prepared for migrating northward under climate warming (Miller et al., 2020). Artemisia tridentata occurrences across the contiguous United States experience very different climatic conditions, and are not representative of the conditions that the northern edge populations are locally adapted to. Additionally, focusing on a single province for this analysis allows for the use of data at a finer scale (such as high-resolution elevation data) that could otherwise not be utilized at a continental-scale due to processing limitations. In total, this deletion left 705 GBIF occurrences. This also removed some occurrences that are likely erroneous, such as lone observations in Alaska, Ontario, and New York city. 

Digital Elevation Models (DEMs) covering the entirety of British Columbia were imported and combined using the ‘Merge Rasters’ tool. The ‘Slope’ and ‘Aspect’ tools were used to calculate slope and aspect from the DEM. Rasters for the present climate and two future climate projections were imported. Specifically, data were extracted for two important climate variables: MAT (mean annual temperature) and MAP (mean annual precipitation). Present-day MAT and MAP data are representative of the average values across the 1991-2020 period. Future climate projection data utilized the CNRM-ESM2-1 Atmosphere-Ocean General Circulation Model (AOGCM). The climate projection data was extracted for the year 2080 under two climate warming scenarios: SSP2-4.5 (from here on, ‘low warming’) and SSP5-8.5 (from here on, ‘high warming’). SSPs, or ‘Shared Socioeconomic Pathways’, represent different hypothetical future scenarios that exhibit varying amounts of climate warming depending on how much society adjusts its greenhouse gas emissions. All six climate rasters (MAT and MAP rasters for the present climate and both future climate scenarios) were clipped to the extent of the British Columbia shapefile. 

The distributions of big sagebrush across these five environmental variables (elevation, slope, aspect, MAT, and MAP) were determined using the ‘Extract Values to Points’ tool. The distribution of occurrences across each of these variables were then plotted. Values of each variable’s mean, maximum, and minimum occurrence value were used to inform the Habitat Suitability Model. For example, the mean MAT value across all big sagebrush occurrences on GBIF was found to be 7.8˚C. The highest value (the maximum) was 10.6˚C, while the lowest (the minimum) was 0.3˚C (Figure 1d). It was also found that there were several outlier occurrences in the elevation data (Figure 1a). Upon examination of these data points, it was determined that most were from a small number of outlying high-elevation populations. It was decided that the mean of the elevation data used for the Suitability Modeller should exclude these outliers so that fuzzy Gaussian distribution used in the modeller better represents the data. To do this, all elevation values above 1400 metres were excluded, changing the mean from 693 metres to 600 metres. However, the maximum threshold used for the Suitability Modeller remained at 2,216 metres, as this value represents the highest-observed occurrence within the dataset. Another alteration to the acquired values was for the aspect data. The mean aspect value was found to be 177˚, but this was changed to 180˚ for simplicity of the calculations. It is likely that the plant does prefer ‘south-facing’ slopes and that the 3˚ difference between direct south (180˚) and the observed mean value (177˚) is due to natural variability in the dataset.

These plots in Figure 1 were also used to estimate the relative ‘importance’ of each environmental variable. While this method is only an estimation, it was deemed that variables with ‘tighter’ correlations that have less variation/noise are more important for dictating the distribution of big sagebrush. By examining the plots in Figure 1, it was decided that slope and MAT had the tightest and least noisy data distribution. These variables are followed by elevation and MAP, which are followed by the least important variable, aspect. Using the AHP (analytical hierarchy process) calculator from the AHP Online System website (AHP Online System, 2021), the hierarchy of these variables was used to calculate their importance (as a percent-weighting) for the Habitat Suitability Modeller. Slope and MAT were defined as ‘moderately more important’ than elevation and MAP, which were defined as ‘moderately more important’ than aspect. This resulted in priority rankings shown in Table 2.

The five variables were used in the Suitability Modeller to determine the spatial habitat suitability of big sagebrush under present climate conditions, future low-warming conditions, and future high-warming conditions. As mentioned previously, each variable’s mean, maximum, and minimum occurrence value were extracted from the data. These three values for each variable were utilized in the Habitat Suitability Model to represent the midpoint and the maximum/minimum thresholds for the fuzzy Gaussian curve in the transformation pane. Table 3 shows all of the values utilized in the fuzzy Gaussian transformations. For all variables except aspect, the values above and below the upper and lower thresholds (respectively) were set to 0. For aspect, the value below the lower threshold (0) was set to 50, as aspect values of -1 represent flat surfaces that are neither north- or south-facing. The weightings from the AHP model (Table 2) rounded to the nearest whole number were used as the weightings of each variable in the Suitability Model. Once all fuzzy Gaussian transformations were performed and the weightings were assigned, the Habitat Suitability Model was run. It should be noted that this was initially attempted with the original high-resolution data, but it was found that the process was too intensive for ArcGIS Pro to perform in a reasonable amount of time. For this reason, each of the five environmental rasters were aggregated using the ‘Aggregate’ tool. The MAT and MAP rasters were aggregated by a cell factor of 4, while the elevation, slope, and aspect rasters were aggregated by a cell factor of 25. The original high-resolution data was still used to acquire the mean, maximum, and minimum values for each variable in order to prevent a loss of accuracy in the Suitability Modeller. The Suitability Model was run three times, once for each climate scenario (present, low-warming, and high-warming). The weightings of each variable remained the same between each model, as did the elevation, slope, and raster values (as these will be mostly unchanged 60 years in the future). MAT and MAP rasters were changed between each of the model runs to reflect the climate scenario being modelled.

After the three models were run, the final Habitat Suitability Model rasters were cleaned up. The rasters were clipped to the exact British Columbia shapefile using the ‘Extract by Mask’ tool. Additionally, the change in habitat suitability was calculated for both low climate warming and high climate warming conditions using the ‘Raster Calculator’ tool. To find the change in suitability, the present-day suitability raster was subtracted from each future climate scenario’s suitability. This generated a raster that showed values that could range from -100 to 100. The values represent how much the suitability of a pixel changes under climate warming, with positive numbers indicating an increase in the habitat suitability and negative numbers indicating a decrease in habitat suitability.