Projected Climate-Fire Interactions Drive Forest to Shrubland Transition on an Arizona Sky Island

O’Connor, Christopher D.; Falk, Donald A.; Garfin, Gregg M.

doi:10.3389/fenvs.2020.00137

ORIGINAL RESEARCH article

Front. Environ. Sci., 21 August 2020

Sec. Interdisciplinary Climate Studies

Volume 8 - 2020 | https://doi.org/10.3389/fenvs.2020.00137

This article is part of the Research Topic Climate, Land Use, and Fire: Can Models Inform Management? View all 8 articles

Projected Climate-Fire Interactions Drive Forest to Shrubland Transition on an Arizona Sky Island

$\r\nChristopher D. O&#x;Connor*&#x;$ Christopher D. O’Connor^1*†

Donald A. Falk^1,2

Gregg M. Garfin¹

¹School of Natural Resources and the Environment, University of Arizona, Tucson, AZ, United States
²Laboratory of Tree-Ring Research, University of Arizona, Tucson, AZ, United States

Climate stressors on the forests of the American Southwest are shifting species distributions across spatial scales, lengthening potential fire seasons, and increasing the incidence of drought and insect-related die-off. A legacy of fire exclusion in forests once adapted to frequent surface fires is exacerbating these changes. Reducing stand densities and surface fuel loads has been proposed as a means of moderating fire behavior while reducing competition for water, but it is not established whether thinning treatments and restoration of surface fire regimes will be enough to offset the multiple manifestations of a changing climate. We examined the potential for prescribed fuel treatments and restoration of historical fire frequencies to mitigate the effects of climate on forest species distributions, composition, total biomass, and fire severity. We used an ecosystem process model to simulate the effects of projected climate, fire, and active management interactions along an ecological gradient of shrublands, woodlands, and forests on a mountain range in Arizona in the United States. We used historical climate conditions as a baseline to compare results from projected climate for the period 2005–2055 with and without fire and with no fuel treatments, a single-entry fuel treatment, and a second fuel treatment after 20 years. Simulated desert grassland and shrub communities remained compositionally stable and maintained or expanded their extents while woodland and forest communities lost basal area and total biomass and receded to the coolest and wettest aspects and drainages even without fire. Initial fuel treatments reduced the extent and relative mortality of high-severity patches for the first two decades, and secondary treatments at simulation year 20 extended these effects for the remaining 30 years of simulation. Immediate and future fuel treatments showed potential to mitigate the severity of fire effects under projected conditions and slow the transition from forest to shrubland in some vegetation types, however, a reduction in basal area and spatial extent of some forest species occurred regardless of management actions. Results are being used to inform local land managers and partners of potential landscape changes resulting from climate alone and from climate–fire interactions and to coordinate active management of fuels across ownerships.

Introduction

Projected warming temperatures and increased moisture variability are likely to cause changes to the frequency and severity of disturbances in many forested ecosystems (Bentz et al., 2010; Abatzoglou and Kolden, 2013; Harris et al., 2016; Riley et al., 2019). In semiarid ecosystems, projected changes to vapor pressure deficit and temperature regimes are expected to significantly increase tree mortality, alter forest species distributions, and limit tree size (Allen et al., 2010; Williams et al., 2010, 2013; McDowell et al., 2011, 2016). However, climate-induced changes to fire and insect-outbreak regimes may multiply and accelerate the effects of climate acting alone by causing rapid tree mortality, soil damage, and changes to landscape structure (Dale et al., 2001; Crimmins and Comrie, 2005; Littell et al., 2010; Keane et al., 2015a). Thus, to understand climate-induced vegetation changes on specific landscapes at fine spatial scales under current and future conditions, it is necessary to capture interactions between biophysical landscape conditions that influence the growth of individual species and the disturbance agents that have historically regulated species and assemblage dynamics.

Climate regulates species geographic distributions profoundly across scales of space and time (Turesson, 1925; Pearson et al., 2004; Rehfeldt et al., 2006). Shifts in species ranges are widespread in the paleoecological record, reflecting evolutionary adaptation to changing climate that continues today (Davis and Shaw, 2001; Colwell and Rangel, 2009; Cole et al., 2013). Topoclimatic and edaphic variation across landscapes accounts for a substantial fraction of variation in species distribution at landscape scales (Zimmermann et al., 2009). Consequently, under the influence of changing climate, species ranges are projected to shift across scales from landscapes to entire species ranges (Chen et al., 2011; Notaro et al., 2012). For sessile species, such as plants, these shifts at any scale are the net demographic result of mortality and recruitment failure at the trailing edge of the distribution (extinction debt) and successful recruitment along the leading edge (immigration or colonization credit) (Jackson and Sax, 2010; Evans et al., 2016; Talluto et al., 2017). Species range shifts from climate pressure alone can occur abruptly from transient climate episodes, such as heat waves (Allen et al., 2015; Ruthrof et al., 2018; Law et al., 2019), but broader changes in species distributions are anticipated to occur over multiple decades, even under the accelerated velocity of anthropogenic climate change (Adams et al., 2009; Burrows et al., 2014).

Ecosystem disturbances, such as fire and insect outbreaks, can accelerate changes in species distributions dramatically. Where range shifts driven by climate alone may unfold over years to decades, severe disturbances can trigger rapid and potentially irreversible ecosystem change (O’Connor et al., 2014; Cobb et al., 2017; Stevens-Rumann et al., 2018; Stevens et al., 2019). High-severity wild-land fire commonly creates large areas of overstory tree mortality, extensive soil damage, and vulnerability to extreme hydrologic events (Neary et al., 1999; Yocom-Kent et al., 2015). Tree seedlings of most conifers disperse typically 100–200 m per generation, so recolonization of large contiguous mortality patches can take multiple generations (Haire and McGarigal, 2010). However, even when seeds reach burned areas, soil and climate conditions may prevent successful seedling establishment (Davis et al., 2019). Species adapted to colonizing post-fire environments may become locally abundant and even dominant, leading to abrupt and persistent type conversion (Savage et al., 2013; Barton and Poulos, 2018; Guiterman et al., 2018). As burned areas increase in size and severity, such conversions resulting from wildfire–climate interactions are likely to become more widespread (Picotte et al., 2016; Reilly et al., 2017; Parks et al., 2018; Singleton et al., 2019).

Post-fire states vary widely, and the mechanisms that govern these transitions are incompletely understood. Resilience is an emergent property comprising multiple component processes acting across scales of space, time, and biological organization (Falk, 2017; Falk et al., 2019). When mortality is locally widespread and soils have been severely altered, combinations of climate and disturbance can result in forest-to-shrubland conversion (Tepley et al., 2017, 2018; Serra-Diaz et al., 2018).

Over the next several decades, the southwestern United States is expected to experience a trend of warming annual mean temperatures (Gonzalez et al., 2018), accompanied by decreasing spring season precipitation in the southern part of the region and increasing percentage of heavy precipitation events (Janssen et al., 2014). For the Southwest region, general circulation models (GCMs) project a 4.80°F (2.67°C) increase in mean annual temperature by midcentury (2036–2065) and an 8.65°F (5.00°C) increase by late century (2071–2100) based on assumptions of high rates of greenhouse gas emissions (RCP 8.5) (Vose et al., 2017). Although changes to precipitation patterns are less certain, particularly for the northern hemisphere summer, projected temperature increases are projected to reduce snowfall water equivalent and the number of snow days (Lute et al., 2015), decrease snowpack through sublimation (Bureau of Reclamation, 2016b), and generate a decreasing fraction of snow compared with rain (Klos et al., 2014; Bureau of Reclamation, 2016a), especially in parts of the Southwest region where seasonal temperatures are near freezing. GCM projections suggest that the region along the United States–Mexico border is likely to experience strong temperature increases, including increases in the number of warm days and decreases in the number of days below freezing (Vose et al., 2017). GCMs project the greatest reductions in winter and spring precipitation in the Southwest for the United States–Mexico border region (Easterling et al., 2017) and show cold season aridification of the border region due to decreasing precipitation (Jones and Gutzler, 2016); border region warm season precipitation and aridification projections remain largely within the historic range of variability. The effects of these rapid changes to regional climate on vegetation (Harpold, 2016), water supplies, and forest disturbances, such as wildfire and insect outbreaks, are not well understood, making the information available to landscape managers in the region insufficient for planning decisions or adaptation.

In ecosystems adapted to frequent surface fires, mechanical treatments that modify the abundance, structure, and distribution of surface and canopy fuels have been shown to reduce water stress (McDowell et al., 2007), restore fire resilience (Fulé et al., 2005; Stephens et al., 2009; Kalies and Kent, 2016), and create stable carbon stores (Ager et al., 2010; Hurteau et al., 2014). Simulation modeling of fire–fuel treatment interactions in the Northern Rockies suggest that increased use of fuel-reduction treatments has the potential to maintain the ecological resilience of forested systems even with high levels of fire suppression (Loehman et al., 2018; Keane et al., 2019). However, results from simulations in pine and mixed-conifer forests of the American Southwest under a range of projected climate conditions and thinning intensities suggests that even a high frequency of thinning treatments did little to offset climate-driven ecological reorganizations, and treatments were only partially effective at reducing fire severity (Loehman et al., 2018). Little is known about the efficacy of thinning treatments in more biologically diverse systems in which fire often propagates across bioclimatic zones and climate effects may be more acute on species already residing near the edge of their bioclimatic envelopes.

We parameterized a species-level landscape simulation model to examine the effects of climate–fire interactions on landscape vegetation communities in a Sky Island mountain range in southeastern Arizona. We also evaluated the potential of fuel treatments to mitigate climate–fire interactions and resulting vegetation type change. Simulation goals were twofold: first, to assess the sensitivity of forest ecosystems, fire effects, and carbon stores to 50 years of projected future climate and, second, to test the potential for proposed fuel-reduction treatments to alter the rate and degree of forest changes through modification of competitive interactions and moderation of fire behavior.

Materials and Methods

Study Location

The Huachuca Mountains are a Madrean Sky Island mountain range, approximately eight kilometers (five miles) north of the United States–Mexico border. Vegetation is distributed along gradients of elevation and aspect starting with a mix of Chihuahuan desert scrub and mesquite grasslands near the base elevation of 1199 m (3934 ft). Along the foothills and shoulders of the range, extensive Madrean encinal woodlands along the southern and western slopes are intermixed with pinyon–juniper woodlands on the flats and northern and eastern slopes. This system transitions to Mexican pine woodland with oak understory at mid elevations and then to nearly pure stands of southern Rocky Mountain ponderosa pine and mixed conifer forest types with intermittent pure aspen stands just below the peak elevation of 2885 m (9466 ft). Annual precipitation is 38 cm (15 in) at the base and 51 cm (20 in) at the peak with mean winter temperatures ranging from 2 to 16°C (35°F–60°F) at the base and −7–2°C (20–35°F) at the peak and mean summer temperatures ranging from 18 to 35°C (65–95°F) at the base and 10–24°C (50–75°F) at the peak with ranges representing average low and high temperatures, respectively (Brown and Comrie, 2002; Morehouse et al., 2006). Ownership is split among the USDA Forest Service, United States Army, private lands, National Park Service, and The Nature Conservancy.

Prior to EuroAmerican settlement in the late 19th century, forests and grasslands of the Huachuca Mountains were shaped by a frequent, typically low-severity fire regime (Danzer, 1998; Barton, 1999; Swetnam et al., 2001). Establishment of a permanent EuroAmerican settlement at Fort Huachuca in 1882 led to the displacement of native tribes and the interruption of several thousand years of naturally occurring and human-augmented wildfires (Danzer et al., 1996). Subsequent expansion of livestock grazing, road building, and resource extraction further reinforced fire exclusion, leading to increasing forest densities and changes to species distributions.

The risk of large, high-severity fires in the Huachuca Mountains has been increasing as human-caused ignitions, fuel loading, changes to forest species and structure, and periods of prolonged drought transition forest ecosystems away from their historical fire-adapted state (Danzer et al., 1996). A series of large and uncharacteristically severe fires burned sections of the range beginning in the early 2000s (the 2002 Ryan fire and 2011 Arlene fire), culminating in the 2011 Monument fire that burned 12,137 ha (29,991 acres) of forest and grassland, of which more than two thirds burned at high severity (75% or more of surface vegetation removed) (MTBS, 2016). Resulting changes to surface cover and soil structure resulted in severe monsoon flooding and secondary damage to homes and other infrastructure following fire containment (Youberg and Pearthree, 2011).

Fort Huachuca has had an active prescribed fire and thinning program in the lower elevation grasslands and woodlands along the eastern front of the range for several decades; however, little has been done to manage fuels in the forested systems that have the greatest potential for damage from wildfires. A 2014 study designed to identify appropriate locations for thinning treatments to mitigate immediate wildfire hazards to listed wildlife species found that targeted fuel-reduction treatments could be used to reduce flame lengths and mitigate other fire behavior characteristics associated with high fire severity and mortality of large trees under assumptions of stable climate (Hollingsworth, 2014). In the analysis detailed in the following sections, we used fuel-treatment locations identified from the 2014 study and thinning prescriptions developed in cooperation with the adjoining National Forest and Fort Huachuca wildlife biologists (Craig Wilcox and Deborah Brewer personal communication).

Simulation Modeling

We used the FireBGCv2 landscape simulation model (Keane et al., 2011) to assess the influence of changing climate on vegetation and fire effects on the Huachuca Mountain landscape. Simulations tracked growth changes in individual trees and shrubs along a 1600 m elevational gradient comprising 10 distinct ecological response units. FireBGCv2 is a tree- to landscape-scale, spatially explicit, ecosystem process model designed for use in montane environments with steep ecological gradients and diverse terrain (Keane et al., 2011). The model tracks the establishment, growth, mortality, and decay of hundreds of thousands of individual trees across a simulated landscape. Disturbance events, such as fire or fuel-management operations, are integrated into the model and influence the growth of trees only on the area of the landscape experiencing the disturbance. On a topographically diverse landscape, such as the Huachuca Mountains, a series of different daily weather streams, modified by elevation, aspect, and topographic index can be applied to adjacent vertically stacked biomes (Figure 1). The model merges vegetation simulation components from FOREST-BGC (Running and Gower, 1991) and BIOME-BGC (Running and Coughlan, 1988; Running and Hunt, 1993; Thornton, 1998), fire initiation and spread outputs from FIRESUM (Keane et al., 1989, 1990), and a series of updated or additional components that simulate weather streams and additional ecosystem processes (Keane et al., 2011).

FIGURE 1

Figure 1. Structure of nested tree, plot, stand, and site layers that make up the FireBGCv2 simulation landscape (Keane et al., 2011). Tree, plot, and stand inputs are developed from field sampling. Site and landscape inputs are developed from the LANDFIRE national vegetation model biophysical setting layer, digital elevation model, and existing vegetation layers (LANDFIRE, 2014). Site-level fire history and fuel data were drawn from a combination of previous data sets (Danzer et al., 1996; Miller et al., 2003) and new field-collected data.

The FIRESUM model in FireBGCv2 uses a simplified, spatially explicit cell percolation algorithm to simulate fire spread, pixel-level fuel parameters to simulate fire intensity, and species-specific physiological traits to determine fire effects on individual trees (Keane et al., 1989). The fire-spread algorithm is simpler and more computationally efficient than that used in FLAMMAP (Finney, 2006), but it still incorporates topographic influences and wind speed and direction to simulate realistic fire progression. Vegetation parameters, fire effects, and allocation of biomass are calculated at an annual time step.

Model Inputs and Species Calibration

Model inputs for species, tree, stand, fuel, and site files were generated from a combination of vegetation and fuel plots, shared databases on southwestern species, and published literature on species-specific ecophysiological parameters and fuel traits. Plot-based data from 156 vegetation and fuel plots (Miller et al., 2003) and 27 supplemental 500 m² forest inventory and age plots were used to validate and adjust biophysical setting maps from LANDFIRE (2014) to develop geo-referenced species and stand databases (Figure 2) and to populate fuel parameters. The network of supplemental plot locations for additional field sampling was targeted for stand types underrepresented by the original vegetation and fuel plot network. The methodology of stand-type determination is detailed in the next section. Supplemental plots were measured over the summer of 2014 and used the same sample protocols as the 2003 vegetation plots to record tree species, diameter, height, canopy base height, and estimated tree ages. Age estimates were based on diameter-age relationships for each species, developed from demographic reconstructions within plots and in the nearby Pinaleño Mountains (O’Connor, 2013).

FIGURE 2

Figure 2. Distribution of existing ecological response units and sampling plots for model calibration. Each of 10 ERUs received a different daily weather stream based on elevation and solar exposure. Individual tree counts and size distributions from sampled plots were used to develop the simulation landscape.

We developed a database of species parameters for the 16 most common tree, shrub, and grass components in 10 ecological response units (ERUs) representing Chihuahuan desert scrub, semi-desert grassland, Madrean encinal woodland, Madrean pine–oak woodland, pinyon–juniper grassland, ponderosa woodland, mixed-conifer forest, aspen woodland, montane riparian woodland, and non-vegetated/developed (LANDFIRE, 2014). Population-level species parameters (e.g., maximum diameter, maximum age, maximum height) were calculated from field-collected plot measurements and life history descriptions. Physiological parameters (e.g., bark thickness) and physiological tolerances for each species were developed from a series of databases maintained by the United States Forest Service Fire Lab (R.A. Loehman et al. unpublished) and the Ecological Restoration Institute at Northern Arizona University (D. Laughlin unpublished) as well as more general parameters published in Silvics of North America (Burns and Honkala, 1990), and BiomeBGC tables (White et al., 2000; Korol, 2001; Hessl et al., 2004).

Individual species included for modeling were mesquite (Vachellia farnesiana (L.) Wight & Arn.); alligator juniper (Juniperus deppeana Steud.); Mexican pinyon (Pinus cembroides Zucc.); pointleaf manzanita (Arctostaphylos pungens Knuth); a complex of evergreen oaks, including Arizona white (Quercus arizonica Sarg.), silverleaf (Quercus hypoleucoides A. Camus), netleaf (Quercus rugose Née), and associated scrub oak (Quercus turbinella Greene), a complex of broadleaf riparian species, including sycamore (Platanus wrightii S. Watson), walnut (Jugulans major Torr.), bigtooth maple (Acer grandidentatum Nutt.), cottonwood (Populous fremontii S. Watson), and willow (Salix gooddingii C.R. Ball); velvet ash (Fraxinus velutina Torr.); Gambel oak (Quercus gambelii Nutt.); Chihuahua pine (Pinus leiophylla Schiede & Deppe); Apache pine (Pinus engelmannii Carrière); ponderosa/Arizona pine (Pinus ponderosa Engelm. var. arizonica); white fir (Abies concolor (Gord. & Glend.) Lindl. ex Hildebr.), Douglas fir (Pseudotsuga menziesii (Mirb.) Franco var. glauca (Beissn.) Franco), southwestern white pine (Pinus strobiformis Engelm.), aspen (Populus tremuloides Michx.), and mixed grasses. In addition to tree-form vegetation, we developed understory models for shrub-form evergreen sumac (Rhus virens Lindh. ex A. Gray var. choriophylla (Wooton & Standl.) L.D. Benson), scrub oak, manzanita, mesquite, and New Mexico locust (Robinia neomexicana A. Gray).

We populated the simulated Huachuca Mountain landscape with forest, shrubland, and grasslands representative of the 183 sample plots. ERUs were further differentiated into 46 stand types representing differences in height class and aspect. At model initiation, the Huachuca landscape had 3141 unique stands differentiated by ERU, height, and aspect (Supplementary Figure S1).

Model Calibration

The calibration modeling weather stream was drawn from 46 years (1961–2007) of continuous daily weather from Coronado National Monument (Western Regional Climate Center, 2014) (AZ_Coop station ID HRFA3, LAT 31.34550, LON -110.25410, elev. 1604 m), located in the foothills along the southern flank of the Huachuca Mountains. We used the MT-CLIM program (Hungerford et al., 1989; Thornton and Running, 1999) to project the Madrean Enceneal woodlands weather stream onto the nine additional ERUs, accounting for topographic and lapse rate effects. Precipitation for each ERU was calibrated to the 30-year normal at each ERU elevational band (PRISM, 2013).

Initial species calibrations were based solely on vegetation succession dynamics. We simulated 300 years of vegetation growth under 20th century climate without fire to assess simulated species dynamics along gradients of moisture, temperature, and interspecific competition. Species parameters were further adjusted to reflect physiological limits and competitive interactions among species that were observed in sampled plots. Multiple runs of identical initiation conditions yielded a range of results over 300 years of simulation because mature tree seed production and dispersal, seedling survival, and tree mortality are simulated stochastically from an independent probability distribution for each species (Keane et al., 2011). Species parameters were considered stable enough to move to the next calibration phase when 80% or more of modeling runs resulted in species spatial distributions and assemblages representative of late successional development. For example, after 300 years of simulation, stands that were originally ponderosa and Mexican pine dominated with a white fir understory component, matured into nearly pure stands of shade-tolerant white fir with a few old remnant pine; and aspen stands were replaced with a mix of more shade-tolerant Douglas fir, white fir, and southwestern white pine.

Once species parameters were calibrated to the range of moisture and temperature conditions across the landscape, we calibrated fire dynamics based on a 400-year reconstruction of fire history on the modeled landscape (Danzer, 1998). Median fire return intervals and fire sizes were used as initial site file fire parameters. Stand- and site-level fuel depths were generated from plot measurements, and fuel model classifications and initial inputs were drawn from Anderson (1982). Inclusion of fire in the model resulted in a slight reduction in stand biomass and a conversion from dense, shade-tolerant forest types to more open fire-adapted species complexes representative of early 20th century forest conditions (Supplementary Figure S2). Simulated fire behavior and resulting fire effects arise from the fire behavior and linked fire effects modules.

Analysis Area as a Subset of Total Simulation Area

Dynamics of vegetation and fire were simulated over the entire Huachuca Mountain landscape to allow fire spread and species emigration across the whole of the elevation gradient and among ERUs. To assess the effects of climate, fire, and fuel treatments on species dynamics and fire effects, we limited the results analysis area to a subset of the landscape incorporating the 10 ERUs in the immediate vicinity of fuel treatments (Figure 2).

Selection and Processing of Climate Projections

Modeling and climate model assessment for this research were conducted prior to the release of the fourth National Climate Assessment (NCA4) (Gonzalez et al., 2018; Reidmiller et al., 2018), so all comparisons are made to the third National Climate assessment (NCA3) (Garfin et al., 2014). To simulate changing climate in a region where precipitation patterns are dominated by the North American Monsoon (NAM), we used the subset of three CMIP5 GCMs that had the lowest error rates for NAM prediction from 1975 to 2005 (Sheffield et al., 2013). The second-generation Canadian Earth Systems Model (CanESM2); the Hadley Centre Global Environment Model, version 2–Carbon Cycle (HadGEM2-CC); and the Hadley Centre Global Environment Model, version 2–Earth System (HadGEM2-ES) are run at coarse spatial resolution (on the order of 1–2° latitude and longitude) and required downscaling for application to the study landscape. The high density of local weather stations available to develop GCM transfer functions for downscaling led us to select the Multivariate Adaptive Constructed Analogs (MACA) statistically downscaled product (Abatzoglou, 2013) to project GCMs onto the Huachuca Mountain landscape at a spatial resolution of 4 km.

We used a 50-year time horizon for model projections for several reasons. The planning and funding horizons of federal partners managing the landscape are typically 10 years or less, so their greatest interest was in near-term climate risk. Additionally, the transfer functions used for statistical downscaling that help to constrain spurious results assume constant circulation relationships that are less likely to hold over simulations exceeding 50 years (Pielke and Wilby, 2012).

In consultation with land managers, we selected the IPCC Representative Concentration Pathway (RCP) with a radiative forcing of 8.5 W/m² in the year 2100 for our modeling runs. This pathway represents the radiative forcing effect of no proactive reduction of global greenhouse gas (GHG) emissions (i.e., global emissions policy characterized as “ongoing high rates of GHG emissions”), which may be a conservative estimate of actual greenhouse gas concentrations later in the century if little action is taken globally to reduce GHG emissions. This scenario also was designed to identify climate vulnerabilities and allow for assessment of viable mitigating actions.

We assessed calibration between projected and historical weather streams through a comparison of the temporal overlap between data sets. Downscaled weather streams projected onto the Madrean encinal ERU from 2005 to 2007 were not statistically different from seasonal mean, maximum, and minimum temperatures and average seasonal precipitation measured at Coronado National Monument weather station (Figure 3).

FIGURE 3

Figure 3. Comparison of historical weather station and calibrated downscaled GCM outputs. Boxplot comparisons of (A) daily maximum temperature, (B) minimum temperature, (C) average temperature, and (D) precipitation. Historical data are from Coronado National Monument, Arizona. GCMs are CMIP5 models HadGEM2-ES, HadGEM2-CC, and CanESM2.

Projected climate conditions from the MACA model ensemble from 2005 to 2055 were consistent with projections from the NCA3 for the Southwest United States (Garfin et al., 2014). Modeled average daily winter temperatures increased, approximately 2.5°C (4.5°F) in the first half of the 21st century with a similar trend to that reported in NCA3 in daily minimum and maximum temperatures (Figure 4A). Modeled daily mean summer temperatures exhibited a slightly lower rate of increase by midcentury than NCA3 estimates, rising 1.5°C (2.7°F) on average with greater variability in daily high temperatures and lower variability in daily low temperatures (Figure 4B).

FIGURE 4

Figure 4. Projected temperature and precipitation for the Huachuca Mountains, AZ, from 2005 to 2055 used for landscape model simulations. Winter temperature (A) is the daily average, and precipitation (C) is the daily total for December, January, and February. Summer temperature (B) is the daily average, and precipitation (D) is the daily total for July, August, and September. Projections are generated from the Multivariate Adaptive Climate Analogues (MACA) statistical regional downscaling of an ensemble of three CMIP5 GCMs using the RCP 8.5 scenario (Abatzoglou, 2013). The GCM subset includes the Hadley Centre HadGEM2-ES and HadGEM2-CC, and Canadian CanESM2 models.

Projections of seasonal precipitation produced by the MACA ensemble suggest relatively little change in the total amount of winter precipitation for 2005–2055 and continued inter-annual variability with more potential for winter drought and flooding years by midcentury (Figure 4C). Monsoon precipitation exhibits a slight increase in total volume and variability of precipitation although significant divergence from historical monsoon volume occurred in only one (CanESM2) of the three GCM projections (Figure 4D).

Fuel-Treatment Details

Fuel-treatment scenarios assume thinning of 500 ha per year starting at year 1 and continuing for 10 years for a single-entry thinning and at years 1 and 20 for a double-entry thinning. Treatment locations were identified in Hollingsworth (2014) and reflect stands with road access and pre-thinning basal area of 4–32 m²/ha. Treatments removed 40% of the basal area, preferentially targeting small-diameter stems up to a maximum diameter of 25 cm DBH with an assumption of 20% of slash left on-site (Figure 2).

Model Simulation Scenarios

Following model calibration, we set up a series of climate change risk scenarios for the 50-year period from 2005 to 2055 to assess potential effects of changing climate conditions, fire, and thinning treatments on forest conditions. We assessed the spatial distributions of dominant forest species, total basal area, and total ecosystem carbon under scenarios of no fire, with fire, no thinning, single-entry thinning, and double-entry thinning. The experimental design included 12 runs of each climate scenario–fire factor–treatment combination, resulting in 288 total landscape simulations. Historical climate results were summarized from 12 replicates of each scenario, and projected climate ensemble results were summarized from 36 replicates (12 from each of three GCMs).

For change analysis, landscape simulation outputs were compared to a baseline case of 50 years of landscape simulation with no fire exclusion under historical climate conditions (1960–2010). We used annual total ecosystem carbon to track gross carbon dynamics as a general summary of climate, fire, and thinning effects on woody biomass over the whole of the simulation area and then used decadal summaries constrained to the analysis area to track fine-scale interactions between climate, vegetation, fire, and fuel treatments.

For each pixel of the analysis area, categorical variables, such as dominant species by biomass, were summarized from the mode of model replicates, whereas continuous variables, such as stem basal area and fire-caused mortality, were calculated from the median value of model replicates. To assess changes in median basal area along the ecological gradient of vegetation types, we grouped the 10 initial site simulations into four general vegetation zones based on the majority of vegetative forms at the start of simulations. We report changes in basal area at decadal time steps in shrubland, woodland, forest, and riparian zones.

Fire-caused mortality, a proxy for burn severity, is calculated at an annual time step in Fire BGCv2. We summarize these results across modeling runs by standardizing to a decadal rate of mortality (median percent mortality) over the 50-year simulation period. Results were further summarized by vegetation zone to assess the interactions between fire, climate, vegetation type, and fuel treatments along the ecological gradient. We present median and maximum fire-caused mortality rates within vegetation zones.

Results

At the landscape scale, total ecosystem carbon (TEC), a proxy for biomass, was influenced strongly by climate and fire interactions (Figure 5). In simulations under historical climate and without fire, TEC gradually increased over the first four decades before plateauing near 7 kg/m². Inclusion of fire in simulations of historical climate resulted in a dynamic equilibrium between fire-related losses and new growth at ∼3.5 kg/m². In climate ensemble projections without fire, TEC increased slowly for the first decade and then leveled off with some minor variability at around 4 kg/m² for the remainder of the simulation. Inclusion of fire resulted a similar trend of dynamic equilibrium in TEC observed under historical climate for the first 20 years of simulations, followed by a slow, steady decline in TEC, reaching what may have been a new equilibrium point for the last two decades of simulation near 2.8 kg/m². The effect of thinning treatments on TEC under projected climate and without fire was negligible with no discernable differences between simulations with or without fuel treatments. Thinning treatment effects on TEC when fire was included with projected climate were somewhat counterintuitive. Under a single fuel treatment, TEC values were similar to those of historical and projected climate values for the first 20 years before declining at a faster rate than in the ensemble TEC scenarios without thinning. The reduction in TEC slowed near the end of the simulation period at a value of 2.2 kg/m². This effect was further amplified with the inclusion of a second thinning treatment, for which TEC values diverged below those of historical and projected climate scenarios within the first 15 years but reached a relative equilibrium at ∼2.1 kg/m² starting around simulation year 35 and continuing to year 50.

FIGURE 5

Figure 5. Change in total ecosystem carbon over 50 years of simulation under historical (Hist) and projected model ensemble (Ens) climate with and without fire (Fire, No Fire) and with no fuel treatments, a single treatment (T1x), and double-entry fuel treatment (T2x).

Species Response to Climate, Fire, and Fuel Treatments

Under historical climate conditions and without thinning treatments, fire had little effect on the composition of tree-dominated communities. However, there were shifts in shrubland communities away from evergreen oak and toward manzanita species as well as a general reduction in the distribution of pinyon–juniper communities (Figure 6A). In climate ensemble projections, climate alone was a potent driver of vegetation changes with fire promoting species diversity and thinning treatments prolonging the presence of a suite of fire-adapted tree species.

FIGURE 6

Figure 6. Species dynamics (relative abundance) over a 50-year time series with varying climate, disturbance, and management actions. Historical climate with fire (A) is used as a baseline for comparison to projected climate ensemble without fire and (B) no fuel treatments, (C) single-entry thinning, and (D) double-entry thinning as well as projected climate with fire and (E) no fuel treatments, (F) single-entry thinning, and (G) double-entry thinning.

Changes to Dominant Species Without Fire

Under projected climate alone, significant changes to forest species composition began to occur at simulation year 20 when evergreen oak and manzanita communities began to encroach on formerly tree-dominated forests. By simulation year 30, Mexican pine, ponderosa pine, pinyon, juniper, and white fir forests were functionally extirpated; southwestern white pine retained a small population; and the distribution of Douglas fir and riparian species increased nominally for a decade before declining as well (Figure 6B). A single-entry thinning treatment in climate-only projections prolonged the retention of Douglas fir, white fir, and ponderosa pine for an additional decade, but other species trends were unchanged (Figure 6C). A second-entry thinning treatment prolonged retention of white fir but reduced the abundance of Douglas fir and ponderosa pine (Figure 6D).

Changes to Dominant Species With Fire

Under projected climate with fire, the composition of forest, woodland, and shrub communities changed significantly, but conversion from forest to shrubland was dependent on management actions. Without thinning treatments, Mexican pine and pinyon–juniper species were lost; however, ponderosa pine was retained into the fifth decade of simulations, and the proportion of the landscape dominated by Douglas fir and white fir increased over the simulation period. Shrub and evergreen oak communities expanded modestly from 60 to 80% of landscape area for the first 30 years before reaching an equilibrium with tree-dominated communities (Figure 6E). Single-entry thinning treatments further reduced shrub and evergreen oak encroachment over the simulation period, resulting in a balance of approximately 70% shrub- and evergreen oak–dominated communities and 30% tree-dominated communities (Figure 6F). Forest composition changes with a second-entry fuel treatment mirrored those of no treatment for the first four decades of simulation, followed by a precipitous decline in tree-dominated area and increase in both evergreen oak and manzanita shrublands (Figure 6G).

Climate, Fire, and Treatment Effects on Basal Area Zones

Under historical climate with fire and without thinning treatments, the basal area of upper-elevation forests increased by 10% while vegetation in riparian, shrubland, and woodland zones decreased by 50, 63, and 32%, respectively (Figure 7A and Supplementary Table S1).

FIGURE 7

Figure 7. Effects of climate, fire, and fuel treatments on median basal area within vegetation zones. Stacked bar charts summarize changes to shrubland, woodland, forest, and riparian vegetation zones at decadal time steps. Historical climate with fire (A) is used as a baseline for comparison to projected climate ensemble without fire and (B) no fuel treatments, (C) single-entry thinning, and (D) double-entry thinning as well as projected climate with fire and (E) no fuel treatments, (F) single-entry thinning, and (G) double-entry thinning. Basal area values within each vegetation zone are presented in Supplementary Table S1.

Fuel-management treatments were not effective for retaining basal area under projected climate conditions in any of the four vegetation zones. By simulation year 30, stem basal area within the study area was reduced by more than half in simulations without fire and by more than two thirds in simulations with fire. Climate effects were the primary driver of reductions in basal area in riparian, woodland, and shrubland vegetation zones. Dramatic reductions in riparian (89–91%), woodland (88–100%), and shrubland (73–82%) basal area were consistent with and without fire and across thinning treatments (Figures 7B–G). In the forested zone, fire had the greatest effect on stem basal area followed by marginal effects of thinning treatments. The reduction in forest basal area over the simulation period with fire (40% with no treatments, 32% with single entry, 43% with second entry) was more than double that of simulations without fire (13% with no treatments, 16% with single or double entry) (Figures 7B–G and Supplementary Table S1).

Climate, Vegetation, and Fuel-Treatment Effects on Fire Severity

Projected climate effects on fire-caused mortality (relative fire severity) varied by vegetation zone and fuel-treatment frequency. Although median decadal mortality rates of 5 to 8% were consistent across shrubland, woodland, forest, and riparian systems regardless of treatment, maximum mortality rates within vegetation zones were more sensitive to fuel treatments. Fire-caused mortality declined by 6–7% of the no-treatment rate with a single-entry thinning treatment and by 16–29% of the no-treatment rate with a second-entry thinning treatment over 50 years of simulation (Figure 8).

FIGURE 8

Figure 8. Effects of thinning treatments on fire-caused mortality under projected climate ensemble conditions. Maps depict median decadal mortality rate within the study area under scenarios of no thinning (A), single-entry thinning (B), and double-entry thinning (C). Histograms display median and maximum mortality rates within shrubland, woodland, forest, and riparian vegetation zones.

Fire-caused mortality under historical climate conditions was most strongly associated with vegetation zones in which median decadal mortality rates in woodland and riparian zones were insensitive to treatments, and shrubland and forest rates ranged from 6 to 10% depending on treatments (Supplementary Figure S3). Changes to maximum severity rates within vegetation zones showed no trend.

Discussion

The combined effects of climate, disturbance, and land management are reshaping many ecosystems in North America. Each of these three primary factors influences forest composition, structure, and distribution individually, and their interacting effects constitute a powerful influence on current and future forests (Cobb et al., 2017; Schoennagel et al., 2017).

On the simulated landscape, forests of the Huachuca Mountains are projected to undergo significant shifts in biomass, species distributions, and patterns of fire over 50 years of projected future climate. Total landscape biomass, summarized as ecosystem carbon, is expected to remain relatively stable over the next decade. However, by 2030, carbon is less likely to be recovered following future fires even with thinning treatments, suggesting that the resilience of current vegetation communities will be moderated by the interaction of fire and climate.

By midcentury, the expansive mid-elevation forests, historically dominated by a multilayered canopy of pine and understory oak, are projected to convert to shrublands even in the absence of fire, and upper-elevation pine and mixed conifer forests are expected to lose more than a third of their basal area and species diversity. These rapid and lasting changes to basal area and ecosystem carbon suggest that forests of the Huachuca Mountains may transition from a carbon-neutral system to a significant carbon source over the coming decades.

Fuel treatments demonstrate the potential to limit the extent and severity of fire-induced mortality in a range of vegetation types. Simulation of a secondary thinning treatment at year 20 further reduced the level of fire-induced mortality. The combination of thinning treatments and fire may also have reduced competition among trees, allowing larger, older trees to persist on the landscape longer than in forest without fire. However, changes to fire dynamics, either through fire exclusion or fuel reduction treatments, did not slow the rate of landscape-scale biomass loss or changes to species distributions, which were still driven inexorably by climate.

Considerations for Interpreting Ecosystem Responses to Projected Climate and Disturbance

Modeling complex landscape interactions of climate and disturbance remains a challenging frontier in ecological research (Keane and Finney, 2003; Keane et al., 2015b). The model ensemble in this study was calibrated to decade-scale trends over the three decades prior to the period of model projections under the specified GCMs and regional downscaling method. Although general trends are similar among GCMs, the degree of uncertainty in the projections increases greatly at annual or shorter time steps (Hawkins and Sutton, 2009). Model agreement is highest for decade-scale trends in seasonal temperature. Climate scientists are less confident in trends in seasonal precipitation due to the wide array of estimates among multiple GCM projections, especially in the region of the North American Monsoon (Easterling et al., 2017). Results from this modeling simulation, although relying on the best available suite of GCMs for the Southwest monsoon region, should be interpreted with caution. The simulation is not a forecast, but rather a projection based on particular assumptions about future global greenhouse gas emissions and is constrained by the limited statistical robustness associated with using a small suite of GCMs to make projections. Nevertheless, the changes to vegetation and fire effects simulated here may be useful for understanding trends in landscape change.

Assumptions about temperature and precipitation inherent in the simulation weather stream underlie many of the physiological stressors and fuel-curing conditions that initiated shifts in forest species and vegetation structure over the simulation period. Warmer winter temperatures projected in the GCM ensemble would reduce the number of days with snowpack and shorten the critical snow melt season. Although snowpack is not simulated directly in the FireBGCv2 model, temperature-driven changes to soil moisture and species-specific drought stress thresholds within each of the 10 topographically distributed biomes capture some proportion of the negative effects of drought stress on high-elevation conifer species, which are dependent upon snow melt for spring bud break, cambial division, and wood formation and were responsible for a significant amount of the basal area loss from the modeled system.

Under the specified climate, fire, and treatment scenarios, the effects of changing climate overwhelmed any benefits of fuel-reduction treatments for reducing water stress although these did reduce the potential for high-severity fire (Hurteau, 2017; Loehman et al., 2018). Although the species physiological parameters developed for this landscape performed well in the calibration weather stream, actual species responses to future climate conditions are inherently uncertain because limited information is available regarding field- or laboratory-quantified drought or heat thresholds of Madrean forest species.

The FireBGCv2 model was designed for use as a research tool and for guiding landscape-management strategies but not individual management decisions, such as treatment locations or specific interventions. The identification of novel changes to landscapes, such as the conversion from conifer to evergreen oak dominance of the mid-elevation forests approximately 20–30 years into a future climate scenario, is clearly a model result that warrants further study.

The statistical methods used for regional downscaling of the GCM model ensemble are considered most appropriate for short-term climate projections. The series of cascading changes to forest species and basal area in this series of simulations occur near the maximum threshold of realistic use of statistical downscaling. Additional modeling runs using a suite of different GCMs and different emissions scenarios as well as different modeling methods, such as dynamic downscaling (Chang et al., 2015; Shamir et al., 2019) would be useful for comparing trends in projected climate effects on species distributions and fire under a range of treatment options.

Using simulation modeling, Keane et al. (2018) used comparisons between historical and future variability to assess ecosystem resilience to climate and different levels of fire suppression on landscapes of the northern Rocky Mountains. Application of this type of modeling for assessing management options under a changing climate is increasing as managers recognize that the historical range of variability may no longer apply to current or future conditions (Crausbay et al., 2017; Falk, 2017; Loehman et al., 2018). Additional work of this kind will help identify the environmental threshold responses of individual species in this landscape and the relative probability of conditions likely to surpass these thresholds. Such inquiries will be valuable to both scientists and managers to reduce the number of ecological surprises associated with future climate and disturbance.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

All authors conceived the project and contributed to the article and approved the submitted version. CO’C performed the modeling, analysis, and manuscript writing. DF and GG provided specific sections of the text and helped with edits, structure, and interpretation of results.

Funding

This research was funded by the Strategic Environmental Research and Development Program (SERDP), an environmental science and technology program supported by the United States Department of Defense and co-administered by the Department of Energy and the Environmental Protection Agency, SERDP Project RC-2232.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank the Natural Resources staff at Fort Huachuca, Sheridan Stone, Debbie Brewer, and Shane Hall, for their expertise in Madrean forest ecology and allowing access to the base; the director of the Nature Conservancy Ramsey Canyon Reserve, Brooke Gebow, for providing accommodations and logistical support; LaWen Hollingsworth for sharing fire modeling results as well as Jeremy Weiss for assistance with climate model file manipulation; Elizabeth Bercel for field assistance; and Rachel Loehman, Bob Keane, and Lisa Holsinger for support with parameterizing, running and summarizing results from FireBGCv2.

Dedication

This article is dedicated to the memory of Harry Sheridan Stone, a passionate ecologist, conservationist and colleague who spent decades in service to the ecosystems, stewards, and visitors to the Huachuca Mountains.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2020.00137/full#supplementary-material

FIGURE S1 | Initial stand locations for model initiation under all scenarios. Color differences represent variability in vegetation type, height, and aspect.

FIGURE S2 | Distribution of dominant vegetation (by biomass) after 50 years of landscape simulation under historical climate with fire.

FIGURE S3 | Effects of thinning treatments on fire-caused mortality under historical climate. Maps depict median decadal mortality rate within the study area under scenarios of no thinning (A), single-entry thinning (B), and double-entry thinning (C). Histograms display median and maximum mortality rates within shrubland, woodland, forest, and riparian vegetation zones.

TABLE S1 | Effects of climate, fire, and thinning treatments on median basal area (m²/ha) within vegetation zones. Scenarios include historical climate with fire and without thinning (Hist No Thin) for reference against climate ensemble projections (Ens) with and without fire and without thinning, single-entry thinning (Thin 1x), and double-entry thinning (Thin 2x).

References

Abatzoglou, J. T. (2013). Development of gridded surface meteorological data for ecological applications and modelling. Int. J. Climatol. 33, 121–131. doi: 10.1002/joc.3413