A fuzzy logic decision support model for climate-driven biomass loss risk in western Oregon and Washington

Dynamic global vegetation model (DGVM) projections are often put forth to aid resource managers in climate change-related decision making. However, interpreting model results and understanding their uncertainty can be difficult. Sources of uncertainty include embedded assumptions about atmospheric CO2 levels, uncertain climate projections driving DGVMs, and DGVM algorithm selection. For western Oregon and Washington, we implemented an Environmental Evaluation Modeling System (EEMS) decision support model using MC2 DGVM results to characterize biomass loss risk. MC2 results were driven by climate projections from 20 General Circulation Models (GCMs) and Earth System Models (ESMs), under Representative Concentration Pathways (RCPs) 4.5 and 8.5, with and without assumed fire suppression, for three different time periods. We produced maps of mean, minimum, and maximum biomass loss risk and uncertainty for each RCP / +/- fire suppression / time period. We characterized the uncertainty due to RCP, fire suppression, and climate projection choice. Finally, we evaluated whether fire or climate maladaptation mortality was the dominant driver of risk for each model run. The risk of biomass loss generally increases in current high biomass areas within the study region through time. The pattern of increased risk is generally south to north and upslope into the Coast and Cascade mountain ranges and along the coast. Uncertainty from climate future choice is greater than that attributable to RCP or +/- fire suppression. Fire dominates as the driving factor for biomass loss risk in more model runs than mortality. This method of interpreting DGVM results and the associated uncertainty provides managers with data in a form directly applicable to their concerns and should prove helpful in adaptive management planning.


EEMS fuzzy logic modeling
The Environmental Evaluation Modeling System (EEMS) [38] is a fuzzy logic [39][40] modeling platform designed to inform answers to management questions. A model is represented by a logic tree, with each node corresponding to a displayable spatial layer or map (e.g. Fig 3). The bottom-most nodes in the tree represent input data layers. Each input layer is first normalized (0 to 1 for this study) to produce a node representing its level of agreement with a user-defined statement. For example, a fuel load metric might be mapped to the statement Simulated Live Biomass is High using user-defined thresholds to characterize High. Normalized values are combined into higher level nodes using fuzzy logic operators that evaluate the relationship between two or more datasets to another statement. For example, data for Simulated Live Biomass is High might be combined with data for Vegetation Stress is High to create a resulting node for Mortality Risk is High. In a complete model, nodes are repeatedly combined to produce a final, top-level node that informs the original management question.
Formally, each node in a fuzzy logic model corresponds to a factual statement, and the values for the node (the normalized values described above) are the values for the statement's fuzzy truth. Fuzzy truths range from 0 for fully false to 1 for fully true. Values between 0.0 and 0.5 are considered partially false, 0.5 is neither true nor false, and values between 0.5 and 1.0 are partially true. Informally, values in the nodes are considered as indices for the attribute associated with the factual statement. For example, a fuzzy value for Vegetation Stress is High might be referred to simply as the level of Vegetation Stress from low (0) to high (1). We use the informal node labels hereafter.
The spatial datasets used in an EEMS model must share the same extent, projection, and reporting units (normally either polygons, or grid cells as in this study). Operations are performed using corresponding reporting units from different data layers (Fig 4A). Reporting units within layers are treated independently of one another and do not influence each other's values.
The EEMS fuzzy logic operators used in this project are And (minimum value of the inputs), Or (maximum value of the inputs), and Union (mean value of the inputs). With the And and Or operators a reporting unit's result value comes from only of the input values (unless multiple input values yield the same minimum value (for And) or maximum value (for Or)). For example, if corresponding cells from nodes A and B have values of 0.3 and 0.5, and these nodes are inputs to the And operator to produce node C, C's value for the corresponding cell would be 0.3 and would come from only the cell in node A. A result of this is that the values in cells of a node produced by And or Or can be attributed to their source node.

Decision support modeling
We created an EEMS decision support model (Fig 3) to evaluate the combined Biomass Loss Risk (formally, Risk of Biomass Loss is High) from fire and climate maladaptation using MC2 ensemble results and aboveground biomass simulated by [34] (hereafter, Hudiburg). For modeled risk to be considered high, the threat to a cell's biomass-either from fire or modeled vegetation type departure from baseline vegetation type-must be high, and the biomass in the cell must also be high. High modeled biomass is insured by input variables in the MC2 Biomass Loss Risk branch of the model. Hudiburg's biomass values are based on observed biomass measurements and their inclusion in the EEMS model serves to adjust Biomass Loss Risk down due to the legacy effects of disturbance and harvest.
Normalization of the datasets to obtain fuzzy values was done by establishing minimum (fully false) and maximum (fully true) thresholds and applying linear interpolation between thresholds, such that where minthresh < = inputval < = maxthresh where inputval > maxthreshwhere fuzzyval is the normalized fuzzy value, inputval is the input (raw) data value, minthresh is the minimum threshold corresponding to the formal node statement, and maxthresh is the maximum threshold.
To normalize MC2 biomass and fire frequency values, we used the distribution of each variable over the study area during the baseline period. The 10 th percentile value for each variable was used as the minimum threshold and the 90 th percentile was used for the maximum threshold (Table 1). Similarly, we normalized Hudiburg's biomass values and used the 10 th and 90 th percentile values from that data set. We calculated MC2 vegetation departure (a shift from the original modeled vegetation type to a new type) by comparing the cell's modal vegetation type for a future period to its vegetation type for the baseline period. A departure value quantifying the level of disparity between past and future vegetation types was obtained from a lookup table based on expert opinion (S1 Table). To normalize the departure values and produce a data layer representing the overall vegetation stress level (MC2 Vegetation Stress) we used departure values of 0 and 3 for minimum and maximum thresholds respectively.

Uncertainty Analysis
We characterized uncertainty by first calculating the variability of Biomass Loss Risk values spatially across each ensemble of results. Results from each ensemble of 20 MC2 runs were combined into 3-dimensional datasets with ensemble members comprising the third dimension ( Fig 4B). Biomass Loss Risk was calculated independently for each ensemble member. An extended version of EEMS was used to produce a data layer for each of the minimum, maximum, and mean fuzzy values (Fig 4B), bracketing the variability. The fuzzy value High Variability (formally, Variability is High) was calculated for each cell in each ensemble by converting standard deviations into fuzzy space using the minimum possible standard deviation (0) as the false threshold and the maximum possible standard deviation (0.5) as the true threshold.
We characterized the non-spatial uncertainty between members of each ensemble using box and whisker plots of their area-weighted means for Biomass Loss Risk. Plots for all scenarios within a time period are displayed together for inter-scenario comparison.
To determine whether climate futures' annual temperature and/or precipitation are tightly coupled with MC2 biomass loss risk, we evaluated the relationships between those 3 variables. First, we compared each ensemble member's area-weighted mean change in temperature from the baseline period against its change in precipitation. Secondly, we compared each ensemble member's contribution to an ensemble's fraction of area matching the maximum MC2 Biomass Loss Risk vs its fraction of area matching the minimum. A visual comparison of an ensemble member's position in the first graph to its position in the second graph illustrates the strength of relationship between these two measures.

Drivers of biomass loss risk
The Or operator in the model node MC2 Biomass Loss Risk (Fig 3) takes the maximum values of the two inputs, one corresponding to the simulated biomass lost by fire from the MC2 model, the other corresponding to the risk of mortality due to vegetation shift (not due to fire) as simulated by MC2. For each ensemble, we took the ensemble mean for each of MC2 Fire Loss Risk, MC2 Mortality Risk, and MC2 Fire Loss Risk minus MC2 Mortality Risk to show which factor most strongly drives MC2 Biomass Loss Risk. Absolute difference values are greatest where one factor produces a high risk and the other produces a low risk. These results reflect the contribution to Biomass Loss Risk from MC2 results without the contribution from Hudiburg Biomass.
We characterized the influence of fire versus that of vegetation shift over the study area. For each ensemble member, we compared the fraction of the area for which mortality due to vegetation shift was the dominant driver of the risk to lose biomass vs the fraction of area where fire was the main driver of risk. Grid cells with a zero risk value were not considered. Level of vegetation departure from historical based on expert opinion (S1 Table).

Decision support modeling
We used normalized biomass values from Hudiburg ( The area weighted mean of Biomass Loss Risk increases with time, is lower for RCP 4.5 than for RCP 8.5, and is slightly higher for NFS scenarios than for FS (Fig 7, Table 2). The range of values increases for all scenarios through time (Fig 7, Table 2).

Uncertainty
Uncertainty is 0 where Hudiburg Biomass is 0, and is generally lower or higher corresponding to lower and higher values for Biomass Loss Risk (Fig 6 and S1-S3 Figs D, H, L). In the Olympic Peninsula, uncertainty is generally higher overall except near the end of the century for the RCP 8.5 scenarios. In the southeastern portion of the study area, uncertainty is low across all scenarios. Area weighted mean uncertainty is similar overall and increases with time ( Table 2).
Between RCP 4.5 and RCP 8.5, uncertainty ranges from 0.01 to 0.09, increasing through time (Table 3). Between FS and NFS, uncertainty ranges from 0.00 to 0.03, with the lowest values for the early 21 st c. (Table 3).

Drivers of Results
MC2 Fire Loss Risk (Fig 8 and S4 The general spatial separation of high MC2 Fire Loss Risk from high MC2 Mortality Risk (Fig 8 and S4-S6 Figs A, B, D, E, G, H) is reflected in the driver difference maps (Fig 8 and S4-S6 Figs C, F, I). The southern portion and northeastern corner of the study area are driven by fire whereas the coast and Olympic Peninsula are driven by mortality. Under RCP 4.5, mortality drives MC2 Biomass Loss Risk more strongly in the Cascades (S4-S5 Figs C, F, I). However, under RCP 8.5, the stronger driver shifts from mortality to fire in the Cascades in the late 21 st c (Fig 8 and S6 Fig. C, F, I).
The fraction of area with 0 MC2 Biomass Loss Risk is somewhat greater for RCP 4.5 scenarios than for RCP 8.5, and declines through time, with a minimum of 11% for any ensemble member (Table 4). Fire (MC2 Fire Loss Risk) contributes more to the risk of losing biomass (Biomass Loss Risk) than does vegetation shift (Mortality Risk) for all ensembles with the exception of RCP 8.5, FS, 2071-2099 (Table 5, Fig 9). The difference is generally smaller for FS scenarios than for NFS scenarios.
Ensemble members with the greatest (or least) change in annual temperature generally do not correspond to climate futures responsible for the largest area of maximum (minimum) risk of biomass loss (Fig 10 and S7 Fig.). One exception is under HadGEM2-ES365 (model 9) which drives the greatest change in temperature for 2071-2099 under both RCP 4.5 and RCP 8.5 and causes the greatest simulated area at risk in the node Maximum Biomass Loss Risk. Over time, the ratio of the number of climate models causing the largest vs the smallest areas at risk of losing biomass (Biomass Loss Risk) increases substantially (2-5/x models for 2011-2030 to 17-19/x for 2071-2099).

Context
We compared our method of presenting uncertainty with those from seven studies covering our study area using inputs from climate models. Climate-based uncertainty was handled in a variety of ways. [15] used the average results from two GCMs, predicating their results on the correctness of those GCMs and the CO 2 projections driving them. [16] used an ensemble mean of results from 17 GCMs as a means of defining a consensus future climate before using that climate in their model. Many studies presented graphs and/or tables of region-summarized values of selected drivers and results [14, 19-20, 22, 25, 33]. Several presented sets of maps allowing visual comparisons of result variation [14, 19-20, 25, 33], but only two presented spatial uncertainty. [19] classified risk to Douglas fir in terms of the percentage of models agreeing or disagreeing on its occurrence, and [20] mapped the number of models agreeing on the direction of change in ecosystem carbon, burn area, and vegetation type. To our knowledge, ours is the first study in this region to provide a detailed, quantified, spatial measure of climate-based uncertainty for modeled future vegetation.

Limitations
An ensemble mean provides a single measure of risk for the ensemble, however climate models driving the results may not be completely independent of one another [41]. Weighting results based on the similarities of the underlying climate models could adjust for this but understanding the provenance of many climate models can be onerous. We did not account for the uncertainty resulting from assumed ignitions or the built-in CO 2 fertilization effect in MC2. The consequences of these assumptions on MC2 fire and carbon dynamics results were found to be substantial [25], but including those uncertainties was beyond the goal of this study. Likewise, we did not incorporate the uncertainty in Hudiburg's [34] data due to the study's limited scope.

General implications
Our results show the risk of biomass loss generally increasing with time in current high biomass areas within the study region. The pattern of increased fire-driven risk through time is generally south to north and upslope as fires become more frequent due to increasing temperatures. Mortality-driven risk increases along the coast where vegetation becomes maladapted to warming and where coastal climate influences reduce fire risk.
Changes in biomass are directly related to ecosystem services such as timber production, carbon sequestration, wildlife habitat provision, recreational opportunities, and fresh water  quality [42]. Thinning, prescribed fire, and suppression can mitigate fire risk, however each of these actions has associated economic and other costs [43][44]. Thinning may increase forest resistance and resilience to drought, however, it may make forests less resistant and resilient as forests age [45]. When physiological processes cannot be buffered against environmental variability, maladaptation leads to mortality [18]. Maintaining biomass in forested areas under climate change-induced maladaptation may depend on management strategies such as sourcing seeds and species from better climatically suited sources (i.e. assisted migration) [46][47]. Tools to help mangers implement such strategies have been developed (e.g. Seedlot Selection Tool, https://seedlotselectiontool.org/sst/).

Management and planning
Our results provide managers with spatial datasets representing three aspects of the Biomass Loss Risk metric. The mean value of biomass lost across all climate futures provides an overall idea of potential magnitude of loss while minimum and maximum values bracket the range, suggesting limits for management alternatives. Uncertainty quantifies the variability of the model results. Land management takes place at multiple scales, with planning and assessment at the national or regional level and implementation at more local levels [48]. Using terminology appropriate for managers such as risk and uncertainty, as well as using a spatial scale appropriate for local information makes our results appropriate for local planning. Environmental models are often found to be insufficiently accurate to use as forecasts [31]. However, they provide insights and scenarios useful for scenario planning [31] and are useful for decreasing uncertainty rather than making predictions [28]. Our work is intended to be viewed within this context, providing one set of results with as much clarity as possible regarding uncertainty, sources of uncertainty, and drivers of risk.
Process-based models, such as MC2, are considered more limited than empirical models in quantifying uncertainty [28], thus limiting their usefulness in management planning. Our method alleviates this limitation, making it easier for managers to use process-based models in their decision making.
It has been suggested that when a range of future possibilities is needed for planning, selecting the most extreme climate projections (e.g. warmest, coolest, wettest, driest) as inputs to ecological models provides brackets for the needed answers [28]. The lack of correspondence we found between the most extreme climate futures and their influence on minimum or maximum risk indicate that simple metrics for climate extremes are not sufficient for bracketing our model results. In a process-based model such as MC2, seasonal patterns and extreme events that are not reflected in annual values or averages over multi-year time periods have the potential to strongly affect fire and vegetation trajectories. Finding climate metrics that predict the most extreme results would be challenging, if not impossible, due to complex interactions within the model. While culling input datasets from GCMs and ESMs that perform poorly in the region may be required to reduce uncertainty [49], culling less extreme climate futures may inadvertently reduce the desired range of results.
Our model may prove useful to managers by itself, but it has the potential to provide greater utility when combined with other metrics reflecting landscape condition, status, and value. The modular nature of the EEMS framework would allow our model to be easily combined into new models. Ignition probabilities, fire spread probabilities, and fire refugia data [56][57][58] could be added to provide greater detail for fire risk. Submodels of climate refugia related to microclimate and enduring landscape features [59][60] could provide more realism for mortality-based risk. Combining our risk model with submodels for current habitat quality (e.g. [61]), connectivity corridors [62], and species presence or absence could help guide management conservation decisions. Similarly, incorporating risk with submodels for economic, social, and cultural values could help managers with biocultural approaches to conservation [63]. Stakeholder input and expert opinion can be used to parameterize these models so that they precisely reflect management concerns. It is our hope that this model and this methodology can contribute to sound decision making for a wide variety of purposes in our study region and beyond.
Supporting information S1 Table. Lookup table for