Projected Wind Impact on Abies balsamea (Balsam fir)-Dominated Stands in New Brunswick (Canada) Based on Remote Sensing and Regional Modelling of Climate and Tree Species Distribution

The paper describes the development of predictive equations of windthrow for five tree species based on remote sensing of wind-affected stands in southwestern New Brunswick (NB). The data characterises forest conditions before, during and after the passing of extratropical cyclone Arthur, July 4–5, 2014. The five-variable logistic function developed for balsam fir (bF) was validated against remote-sensing-acquired windthrow data for bF-stands affected by the Christmas Mountains windthrow event of November 7, 1994. In general, the prediction of windthrow in the area agreed fairly well with the windthrow sites identified by photogrammetry. The occurrence of windthrow in the Christmas Mountains was prominent in areas with shallow soils and prone to localised accelerations in mean and turbulent airflow. The windthrow function for bF was subsequently used to examine the future impact of windthrow under two climate scenarios (RCP’s 4.5 and 8.5) and species response to local changes anticipated with global climate change, particularly with respect to growing degree-days and soil moisture. Under climate change, future windthrow in bF stands (2006–2100) is projected to be modified as the species withdraws from the high-elevation areas and NB as a whole, as the climate progressively warms and precipitation increases, causing the growing environment of bF to deteriorate.


Introduction
Wind plays an important role in defining the functional and structural characteristics of forest ecosystems globally [1][2][3]. Powerful windstorms are particularly notorious for causing widespread ecological and economic losses in wind-impacted landscapes [2]. Wind damage can cause significant challenges for forest managers. These challenges can range from potential fibre shortages in semi-mature forests due to their potential future yield, loss of habitat for some plant and animal species and loss in landscape aesthetics [4]. Salvage operations must contend with increased costs and danger to forest operators [4,5]. An example of such a wind event in New Brunswick (NB), Canada, was the three-hour, Modelling the presence/absence of windthrow for bF and four other tree species was based on windthrow and non-windthrow data from southwestern NB (study area 1, Figure 1) for function development, and northcentral NB (study area 2, Figure 1) for function validation. Windthrow data in study areas 1 and 2 were based on outcomes of extratropical cyclone Arthur, July 4-5, 2014, and the Christmas Mountains windthrow event of November 7, 1994, respectively. Given the predominance of bF in the Christmas Mountains prior to the windthrow event, the function validation was feasible only for bF.

Surface Description and Development
The windthrow functions were constructed from data extracted from: • High-resolution, four-band digital orthophotos taken after the passage of extratropical cyclone Arthur in detection of windthrow-and non-windthrow-affected forests; Modelling the presence/absence of windthrow for bF and four other tree species was based on windthrow and non-windthrow data from southwestern NB (study area 1, Figure 1) for function development, and northcentral NB (study area 2, Figure 1) for function validation. Windthrow data in study areas 1 and 2 were based on outcomes of extratropical cyclone Arthur, July 4-5, 2014, and the Christmas Mountains windthrow event of November 7, 1994, respectively. Given the predominance of bF in the Christmas Mountains prior to the windthrow event, the function validation was feasible only for bF.

Surface Description and Development
The windthrow functions were constructed from data extracted from: • Selected forest-state variables from a suite of GIS-thematic variables of stand type, structure, soil depth and texture and other stand indicators (Table 1), and environmental site variables of digital elevation model (DEM)-derived estimates of slope, height above nearest drainage point [17,18] and modelled relative soil water content. Principal component analysis (PCA) was subsequently used to reduce the number of variables from the suite of forest-state and site variables in developing the windthrow functions for bF and the other tree species.
The windthrow functions formulated for bF were later used to investigate the impact that current and future episodes (1950-2005 and 2006-2100) of extreme wind speeds may have on bF-dominated forests in NB. The future wind impact on bF was set in a context of anticipated species shift, as bF and other tree species will respond to changes in local temperature, precipitation, and soil water content over the next nine decades (based on past work by Bourque [22]). Climate-change scenarios for future environmental conditions for NB (i.e., temperature, precipitation, wind speed and direction) were based on 0.44 • × 0.44 • -resolution climate projections from ECCC's most recent RCM (CanRCM4 [23]) developed by the Canadian Centre for Climate Modelling and Analysis (CCCma). CanRCM4 is driven by its parent global model (CanESM2 [24]) for all downscaling applications by employing spectral nudging in CanESM2. In our case, the simulations were based on two representative concentration pathways (RCP's [25]) projected to contribute to a net additional global radiative forcing of 4.5 and 8.5 W m −2 by 2100 (RCP's 4.5 and 8.5, respectively). The following sections describe the analytical procedures and input information (area extractions, for individual forest stands) used in developing the windthrow maps displayed throughout the paper.  (Figure 1). The DEM served as a data layer central to modelling the likelihood of windthrow province-wide under six extreme wind events.

Wind Speed and Direction
Wind speed and direction were modelled spatially according to the full three-dimensional Navier-Stokes equations, incorporating both the effects of atmospheric turbulence and thermal processes [26]. Computational fluid dynamics model [26] calculations of sustained maximum wind speed ( Figure 2) for southwest (study area 1) and northcentral NB (study area 2; Figure 1) were based on a boundary-fitted coordinate system. Initial boundary conditions were specified by the near-surface air temperature and sustained wind speed and direction recorded at the Fredericton Airport and St. Stephen weather stations ( Figure 1) during the peak of extratropical cyclone Arthur (mean air temperature of 14.4 • C and 13.5 • C and wind speed and direction of 23.2 m s −1 and 16.3 m s −1 and 0 • from true north, respectively) and corresponding measurements recorded at the St. Leonard Airport weather station, upwind of the Christmas Mountains, for the Christmas Mountains windthrow event (mean air temperature of 0.5 • C and wind speed and direction of 13.0 m s −1 and 304 • , respectively [16]). In both instances, a wind speed at 500 m AMSL was specified as corresponding to a mean value of 1.3 × near-surface sustained wind speed (after Franklin et al. [27]). Temperature stratification of the atmosphere was set to neutral for both wind events.  [16]). In both instances, a wind speed at 500 m AMSL was specified as corresponding to a mean value of 1.3 × near-surface sustained wind speed (after Franklin et al. [27]). Temperature stratification of the atmosphere was set to neutral for both wind events. Stephen Areas in southwest NB (a; study area 1, Figure 1) and the Christmas Mountains in northcentral NB (b; study area 2, Figure 1). The black solid line and associated curved error bars in the circular wind graphs in subfigures (a) and (b) represent the mean wind direction and ± 1 standard deviation of hourly wind directions during extratropical cyclone Arthur (July 4-5, 2014, subfigure (a)) and the Christmas Mountains windthrow event (November 7, 1994, subfigure (b)). At the peak of cyclone Arthur, winds at the Fredericton and St. Stephen weather stations ( Figure 1) were mostly from true north (i.e., 0 o ). Prevailing winds during the Christmas Mountains windthrow event were mostly from the west-northwest to northwest direction (i.e., ~304 ±7 o from true north). Forests in the Christmas Mountains at the time of windthrow were mostly composed of bF (yellow polygons, middle panel of subfigure (b)); blue polygons indicate nonforested areas (e.g., cutovers, fields, industrial operations, etc.), whereas the red polygons in the lower panel of subfigure (b) indicate areas affected by windthrow, as resolved from photointerpretation.
During a particular windstorm event, windthrow largely results from the destabilising forces generated by large-to-medium size eddies (wind gusts) formed as the prevailing airflow interacts with the underlying irregular terrain and surface cover. However, in this paper, we refer to sustained maximum wind speeds as the main indicator of windthrow, because (i) sustained wind speeds are always weaker then peak wind gusts (PWG's); (ii) they are easier to model spatially with CFD-models and (iii) PWG's, although intermittent, can be related empirically to sustained wind speeds (e.g., [28]). Stephen Areas in southwest NB (a; study area 1, Figure 1) and the Christmas Mountains in northcentral NB (b; study area 2, Figure 1). The black solid line and associated curved error bars in the circular wind graphs in subfigures (a) and (b) represent the mean wind direction and ± 1 standard deviation of hourly wind directions during extratropical cyclone Arthur (July 4-5, 2014, subfigure (a)) and the Christmas Mountains windthrow event (November 7, 1994, subfigure (b)). At the peak of cyclone Arthur, winds at the Fredericton and St. Stephen weather stations ( Figure 1) were mostly from true north (i.e., 0 • ). Prevailing winds during the Christmas Mountains windthrow event were mostly from the west-northwest to northwest direction (i.e.,~304 ±7 • from true north). Forests in the Christmas Mountains at the time of windthrow were mostly composed of bF (yellow polygons, middle panel of subfigure (b)); blue polygons indicate nonforested areas (e.g., cutovers, fields, industrial operations, etc.), whereas the red polygons in the lower panel of subfigure (b) indicate areas affected by windthrow, as resolved from photointerpretation. During a particular windstorm event, windthrow largely results from the destabilising forces generated by large-to-medium size eddies (wind gusts) formed as the prevailing airflow interacts with the underlying irregular terrain and surface cover. However, in this paper, we refer to sustained maximum wind speeds as the main indicator of windthrow, because (i) sustained wind speeds are always weaker then peak wind gusts (PWG's); (ii) they are easier to model spatially with CFD-models and (iii) PWG's, although intermittent, can be related empirically to sustained wind speeds (e.g., [28]). Peak wind gusts during sustained airflow over land were calculated here from sustained maximum wind speeds (η) and a wind-gust function formulated from wind speeds recorded during fifteen major mid-Atlantic hurricanes, including Hurricane Sandy [29], i.e., where η is the sustained maximum wind speed (m s −1 ). The denominator in the quotient in Equation (1), i.e., 0.44704, allows for a unit conversion from m s −1 to miles hour −1 . The bracketed term in Equation (1) is the gust factor [29].

Soil Water Content
The height above nearest drainage point ( Figure 3a) provides a simple static (time-independent) measure of potential drainage, particularly close to drainage channels [17,30], and is described as the vertical separation between raised dry land and a localised estimate of the water table level based on surface water defined by a GIS hydrological layer. For a particular dry cell, the wet cell nearest it is located by means of an iterative search algorithm that minimises the horizontal distance between the subject dry cell and the wet cell downslope, while taking into account DEM-specified flow directions and pathways [17,18]. Here, HNDP=0 m signifies visible surface water, whereas large HNDP's signify low water tables and potentially drier soil conditions. In upland positions of landscapes, HNDP does poorly as an indicator of soil water content, as there is no realistic connection between HNDP and potentially variable soil water content. To avoid the challenge this presents, we introduced a water budget-based calculation of relative soil water content to improve the representation of soil water content in upland positions of landscapes. Because the water budget-based calculation of the variable depends on hydrometeorological inputs and outputs that may conceivably change over time, its calculations are time dependent.
Combining the two variables in a single representation of soil water content is acceptable, as HNDP approximates the influence on soil water content relative to the vertical distance to drainage points, whereas the soil water content module of Landscape Distribution of Soil moisture, Energy, and Temperature (LanDSET) [19,20] does better at representing soil water content relative to flow characteristics of the terrain (i.e., pathways and pooling locations), particularly in upland positions of the landscape. Together, HNDP and LanDSET-based calculations of relative soil water content have the potential to improve description of surface water flow-related processes and features across forested landscapes [18].

7
poorly as an indicator of soil water content, as there is no realistic connection between HNDP and potentially variable soil water content. To avoid the challenge this presents, we introduced a water budget-based calculation of relative soil water content to improve the representation of soil water content in upland positions of landscapes. Because the water budget-based calculation of the variable depends on hydrometeorological inputs and outputs that may conceivably change over time, its calculations are time dependent.  [18,19] and LanDSET-modelled relative soil water content (b) for a small area in the vicinity of the Greater Fredericton Area, NB. Blue to cyan-(a) and black-coloured areas (b) in the maps coincide with wet areas (i.e., exposed surface water and wetlands); white and red polygons in the maps coincide with forest areas unaffected and affected by windthrow, respectively.
Combining the two variables in a single representation of soil water content is acceptable, as HNDP approximates the influence on soil water content relative to the vertical distance to drainage  [18,19] and LanDSET-modelled relative soil water content (b) for a small area in the vicinity of the Greater Fredericton Area, NB. Blue to cyan-(a) and black-coloured areas (b) in the maps coincide with wet areas (i.e., exposed surface water and wetlands); white and red polygons in the maps coincide with forest areas unaffected and affected by windthrow, respectively. Steady-state values of relative soil water content of southwest ( Figure 3b) and northcentral NB were estimated numerically, following methods described in several scientific articles central to the description of LanDSET [19,20,31,32]. LanDSET addresses soil water distribution according to a DEM cell-by-cell assessment of hydrological fluxes, including precipitation, evapotranspiration, percolation and surface runoff. The input to the model includes topography, by way of DEM elevations and cell size specification (20 m), net radiation balance and annual precipitation described spatially by numerical interpolation of annual precipitation recorded at 32 provincial weather stations. Prior to the calculation of relative soil water content, the raw DEM was processed by removing false depressions with the pit-filling algorithm of Planchon and Darboux [33], accessible in SAGA (System for Automated Geoscientific Analyses [34]).
The hydrological inputs to individual DEM gridpoints were precipitation and lateral flow from upslope regions, and outputs were percolation (deep seepage), evapotranspiration and surface runoff flows to regions downslope. Surface water on valley slopes was channelled downslope according to zones of confluence-divergence in the local landscape. The calculation of gridpoint evapotranspiration was based on an application of the Priestley-Taylor equation [35], gridpoint relative soil water content and LanDSET-model evaluation of net all-wave radiation [20]. The values of relative soil water content were in the range of 0.0-1.0, with 0.0 being associated with the drier sites at or below the permanent wilting point, and 1.0 with the wetter sites at or above field capacity. Percolation was based entirely on relative soil water content and occurs whenever soils are above field capacity [20]. Above field capacity, percolation rates were defined as a function of soil property and the saturated conductivity of the soil. The values for the equation parameters have been listed by Bourque et al. [19,20].

Windthrow Detection
The presence/absence of windthrow in study area 1 (Figure 1) was determined by photointerpretation of 30-cm resolution, colour infrared (CIR, surface reflectance-based) digital orthophotos (e.g., Figure 4) taken after the passing of extratropical cyclone Arthur. Individual bands in CIR imagery consist of a matrix of pixels representing colours in the red, green, blue (RGB) or near infrared (NIR) segments of the electromagnetic spectrum.  Figure 1). False colour was used to bring out forest characteristics, e.g., in images (b) through (d), which would have otherwise been difficult to detect.
Living vegetation absorbs blue and red light to produce chlorophyll and drive photosynthesis. Healthy plants possess more chlorophyll and, as a result, reflect more NIR energy than damaged or dying plants. Consequently, evaluating a plant's surface reflectance in terms of both absorbed and reflected visible light (Red Green Blue, RGB) and NIR energy can convey important diagnostic information about the plants' physiological status.
In the false-colour image in Figure 4d (NIR+RG), light pink to greenish colours correspond to damaged and windblown trees and ground debris (e.g., fallen stems, dead detached foliage and branches and exposed forest floor), whereas bright red colours coincide with healthy, undamaged trees and foliage. Colour rendering of images to false colour, as demonstrated in Figures 4c and 4b, helps to locate windthrow-affected areas.

Windthrow Function Development
The windthrow functions were generated with genetic programming (GP)-driven logistic regression. Logistic regression is a statistical procedure widely used in classifying objects as either being present or absent (i.e., binary in nature), given a specific decision threshold value. Regression values below the threshold value signify absence (giving a value of 0.0), and values above the threshold signify presence (giving a value of 1.0). Here, absence refers to forest stands that are observed/predicted to have been left undamaged by wind, and presence refers to stands that are observed/predicted to have been damaged. All wind impact calculations were carried out at the DEM grid cell level and were not carried out for individual stands or soil formations. Information about stands and soil formations was passed on to individual grid cells by rasterising relevant feature classes (shapefiles) at the same resolution as the base DEM (i.e., at 20-m resolution). Living vegetation absorbs blue and red light to produce chlorophyll and drive photosynthesis. Healthy plants possess more chlorophyll and, as a result, reflect more NIR energy than damaged or dying plants. Consequently, evaluating a plant's surface reflectance in terms of both absorbed and reflected visible light (Red Green Blue, RGB) and NIR energy can convey important diagnostic information about the plants' physiological status.
In the false-colour image in Figure 4d (NIR+RG), light pink to greenish colours correspond to damaged and windblown trees and ground debris (e.g., fallen stems, dead detached foliage and branches and exposed forest floor), whereas bright red colours coincide with healthy, undamaged trees and foliage. Colour rendering of images to false colour, as demonstrated in Figure 4b,c, helps to locate windthrow-affected areas.

Windthrow Function Development
The windthrow functions were generated with genetic programming (GP)-driven logistic regression. Logistic regression is a statistical procedure widely used in classifying objects as either being present or absent (i.e., binary in nature), given a specific decision threshold value. Regression values below the threshold value signify absence (giving a value of 0.0), and values above the threshold signify presence (giving a value of 1.0). Here, absence refers to forest stands that are observed/predicted Remote Sens. 2020, 12, 1177 9 of 21 to have been left undamaged by wind, and presence refers to stands that are observed/predicted to have been damaged. All wind impact calculations were carried out at the DEM grid cell level and were not carried out for individual stands or soil formations. Information about stands and soil formations was passed on to individual grid cells by rasterising relevant feature classes (shapefiles) at the same resolution as the base DEM (i.e., at 20-m resolution).
Genetic programming-based regression is used to determine which nonredundant variables are particularly important in explaining variability in a dependent variable from the list of independent variables (Table 1), which, in our case, was the presence/absence of windthrow. Genetic programming-based regression is a procedure founded on evolutionary computation in search of algebraic equations by reducing the difference between target values and values calculated with the equations generated with the procedure [36]. Different from conventional regression, which determines parameter values of known equations, no specific mathematical expression is needed as a starting point with this approach. Rather, primary expressions are formed by randomly combining primitive, base functions of input variables (linear and nonlinear) with algebraic operators. Equations retained by the procedure are those that replicate the target output data better than others, and undesirable solutions are simply discarded. The procedure stops whenever the desired accuracy in data replication or level of machine learning has been reached.
A nontrivial problem in using GP-based methods is the amount of time sometimes required to arrive to an acceptable solution. One method to reduce GP-processing times is to reduce the number of input variables (Table 1) to a few important explanatory variables by means of data mining techniques, e.g., principal component analysis, regression tree analysis, cluster analysis, artificial neural networks and others [37]. Data reduction with any of these methods help to determine which input variables, among those being considered, explain the most variation in a dataset. Involving redundant and nonexplanatory variables in equation development simply helps to slow down solution convergence.
To help reduce problem dimensionality, we employed principal component analysis (PCA) on the suite of forest-state and environmental-site variables addressed in Table 1. Principal component analysis groups redundant variables under a principal component (or primary axis), with each component explaining a portion of the variation in the dataset. Variables grouped together after rotation of the primary axes (based on the VARIMAX-option in SYSTAT ver. 11) are viewed to contain very similar information. As a result, only one of these variables with highest factor loading or easiest to quantify or measure needs to be retained, and all other variables belonging to the same principal component are ignored. For a particular dataset, the number of viable principal components depends on the eigenvalue produced during axis rotation. Principal components, with an eigenvalue > 1.0, collectively explain the majority of the variation in the dataset.

Projected Future Climates
Future climate change scenarios for NB were based on ECCC's RCM projections (i.e., CanRCM4_CanESM2 at 0.44 • × 0.44 • resolution) based on RCP's 4.5 and 8.5 [25]. RCP 4.5 represents a scenario that stabilises a net rise in global radiative forcing to 4.5 W m −2 prior to the end of the 21st century. This scenario assumes that global annual greenhouse gas emissions peak around 2040, and then decline with the implementation of new technology [25], whereas RCP 8.5 represents a significantly more intensive scenario of rising global radiative forcing, reaching 8.5 W m −2 by the end of the 21st century. RCP 8.5 assumes global emissions continue to increase corresponding to the pathway with the greatest greenhouse gas emissions [38]. Figure 5 gives province-wide windspeed frequencies as a function of RCM-modelled historical climatic regime ; inset to Figure 5a) and seasonal wind for two locations in NB (Figure 5b), i.e., one inland location on the southern boundary of study area 2 (gridpoint 61, near the Christmas Mountains) and one coastal location, near Miscou Island, NB (gridpoint 94; see inset to Figure 5a). There was significant deviation (climatology) in overall wind speeds between the two locations. The sustained maximum wind speeds in coastal areas had a greater chance of exceeding 9 m s −1 , attaining wind gusts > 13 m s −1 , particularly during winter, due to stronger atmospheric circulation. The province-wide evaluation of windstorm impact on bF stands for both historical and future climate scenarios ( Figure 6) was based on six one-time events (two historical and four future events) of RCM-generated maximum wind speed over 56 (historical data) to 95 years (future data) at both inland and coastal areas of NB (gridpoints 61 and 94, respectively; Figure 5a). Given the computational grid area and time limitations of CFD models, achieving the level of spatial detail required for windthrow calculations at 20-m resolution for the entire province was enormously challenging. As an effective sidestep to this problem, we instead enhanced the province-wide RCMbased maximum wind speeds using the statistical downscaling of associated corrections (dependent variable) derived by dividing gridpoint CFD-and RCM-based maximum wind speeds occurring within study area 2. The independent variables to the procedure included RCM-based wind speeds and four DEM terrain-based indices, including wind-effect index, terrain ruggedness index, point maximal curvature and relative elevation. The wind-effect index is a nondimensional quantity that defines the areas of the landscape either exposed to or sheltered from high wind speeds by taking into account spatial variation in digital terrain elevations and interpolated RCM wind speeds and directions using SAGA. The statistical relation of windspeed correction is subsequently used to produce province-wide estimates of windspeed correction based on province-wide surfaces of independent variables at 20-m resolution, serving as inputs. Windspeed corrections are then multiplied to the interpolated RCM-generated wind speeds, increasing their overall spatial resolution. Wind directions were predominantly from the northeast and south in proximity to the eastern and southern coast of NB, and west-to-northwest direction in northcentral NB.
The province-wide evaluation of windstorm impact on bF stands for both historical and future climate scenarios ( Figure 6) was based on six one-time events (two historical and four future events) of RCM-generated maximum wind speed over 56 (historical data) to 95 years (future data) at both inland and coastal areas of NB (gridpoints 61 and 94, respectively; Figure 5a). Given the computational grid area and time limitations of CFD models, achieving the level of spatial detail required for windthrow calculations at 20-m resolution for the entire province was enormously challenging. As an effective sidestep to this problem, we instead enhanced the province-wide RCM-based maximum wind speeds using the statistical downscaling of associated corrections (dependent variable) derived by dividing gridpoint CFD-and RCM-based maximum wind speeds occurring within study area 2. The independent variables to the procedure included RCM-based wind speeds and four DEM terrain-based indices, including wind-effect index, terrain ruggedness index, point maximal curvature and relative elevation. The wind-effect index is a nondimensional quantity that defines the areas of the landscape either exposed to or sheltered from high wind speeds by taking into account spatial variation in digital terrain elevations and interpolated RCM wind speeds and directions using SAGA. The statistical relation of windspeed correction is subsequently used to produce province-wide estimates of windspeed correction based on province-wide surfaces of independent variables at 20-m resolution, serving as inputs. Windspeed corrections are then multiplied to the interpolated RCM-generated wind speeds, increasing their overall spatial resolution.
12 Figure 6. RCM-generated wind speed and direction for six extreme wind events corresponding to the historical period  and two future climate scenarios, RCP's 4.5 and 8.5, over the 2006-2100 period. The extreme events were based on searching modelled daily wind data for highest maximum wind speed for gridpoint locations 61 and 94 (Figure 5a). Time of occurrence of projected events is indicated for individual wind distributions.

Projected Species Shifts
The discussion of province-wide wind impact on bF stands was done in a context of anticipated species shift, as individual bF trees are expected to respond to changes in local temperature, precipitation and soil water content over the next nine decades (i.e., 2006-2100) for future climate scenarios, RCP's 4.5 and 8.5. Table 2 provides the factor loadings associated with the application of PCA to the full set of variables listed in Table 1. From this analysis, eight principal components (PC's) with eigenvalues > 1.0 were identified. Collectively, the PC's explained 79% of the variance in the dataset. Given group membership, factor loadings ( Table 2) and the ease to which variables can be quantified spatially, mean height above nearest drainage (HNDP_MEAN), relative soil water content (SWC_MEAN), wind speed (WNDSPD_MEAN), density class (DC), tree form (TREEFORM), soil depth (SDEPTH) and mean slope (SLP_MEAN) were judged to be the most important explanatory variables, with HNDP_MEAN explaining the greatest variation (14.2%) in the data among the seven variables, followed by wind speed (13.5%; Table 2). Beyond these seven variables, the remaining variables in Table 1 were discounted from further consideration.

Projected Species Shifts
The discussion of province-wide wind impact on bF stands was done in a context of anticipated species shift, as individual bF trees are expected to respond to changes in local temperature, precipitation and soil water content over the next nine decades (i.e., 2006-2100) for future climate scenarios, RCP's 4.5 and 8.5. Table 2 provides the factor loadings associated with the application of PCA to the full set of variables listed in Table 1. From this analysis, eight principal components (PC's) with eigenvalues > 1.0 were identified. Collectively, the PC's explained 79% of the variance in the dataset. Given group membership, factor loadings ( Table 2) and the ease to which variables can be quantified spatially, mean height above nearest drainage (HNDP_MEAN), relative soil water content (SWC_MEAN), wind speed (WNDSPD_MEAN), density class (DC), tree form (TREEFORM), soil depth (SDEPTH) and mean slope (SLP_MEAN) were judged to be the most important explanatory variables, with HNDP_MEAN explaining the greatest variation (14.2%) in the data among the seven variables, followed by wind speed (13.5%; Table 2). Beyond these seven variables, the remaining variables in Table 1 were discounted from further consideration.

Windthrow Function Development
Windthrow equation for bF and associated thresholds (δ i , i=1, . . . ,4 and decision thresholds) were derived from GP-based regression (Tables A1 and A2 in Appendix A; these Tables give the windthrow equations and thresholds associated with bF and the other tree species). The windthrow function for bF (Table A1) was able to correctly classify the majority of windthrow-and non-windthrow-affected stands during function development. The area under the receiver operating characteristic (ROC) curve for bF exceeded 94% (Table A2), with a sample size of 518 of windthrow-and non-windthrow-affected stands. An area under the ROC-curve of 100% indicates a perfect match between observed and predicted outcomes. Here, WNDSPD, HNDP, SWC and SLP refer to their respective means ( Table 1). The windthrow-functions for the other tree species (Table A1) performed equally well with determining outcomes of windthrow, despite being based on smaller sample sizes (Table A2).
Relative to the other tree species, bF is particularly vulnerable to blowdown because of the species' shallow rooting habit particularly in shallow, wet soils. This is supported by the relatively low critical wind speed determined for the species during regression, i.e., 5.8 m s −1 compared to 7.0 m s −1 and 9.5 m s −1 for red and sugar maple (Table A2). Estimated PWG's of 8.8 m s −1 , 10.3 m s −1 and 13.5 m s −1 were associated with the calculated critical wind speeds of 5.8 m s −1 , 7.0 m s −1 and 9.5 m s −1 . Red and sugar maple tended to have deeper rooting systems, particularly in well-drained, deep upland soils as a result of being more strongly anchored. These findings are consistent with other researchers' findings based on different methods of evaluation. For example, after a catastrophic wind event in the Charlevoix region of Quebec, Ruel [39] found that 62% of the total windthrow was located in bF stands. For wind speeds between 7.8-10.3 m s −1 , bF stands were more heavily damaged than mixed black spruce stands. In a study of experts' perceptions of windthrow, experts consistently ranked bF, among thirteen northern temperate tree species, including seven hardwood and six softwood species, as especially prone to windthrow [40]. Burns and Honkala [41] also placed bF as not being very resistant to windthrow. Canham et al. [42] have reported that red maple was generally less windfirm than sugar maple. Our results are consistent with that observation (Table A2).

Windthrow Function Validation
The windthrow function for bF was subsequently validated against windthrow data acquired from the Christmas Mountains blowdown event (Figure 7). Overall, the prediction of windthrow-affected areas in the Christmas Mountains (study area 2; Figure 1) agreed fairly well with the windthrow sites identified from photointerpretation (Figure 7b), as 93.3% of the impact zone established through photointerpretation corresponded with predictions of the occurrence of windthrow. Unmatched areas outside the impact zone indicated the possibility of windthrow having actually occurred based on the prediction of occurrence, although they were not flagged as containing windthrow. On closer inspection of mostly cloud-free infrared and optical-based Landsat images taken before and after windthrow (i.e., September 21, 1994 andJune 4, 1995) and higher-resolution, coloured aerial film photographs taken after windthrow (August 08, 1995), all 25 additional locations that showed a high likelihood of windthrow possessed some level of windthrow but at varying degrees of destruction. In general, the occurrence of windthrow in the Christmas Mountains was more pronounced in areas (i) with shallow soils (SDEPTH class 2, soils 30-65-cm deep; Figure 7c), even in landscape locations of moderate sustained wind speeds (~6 m s −1 ), and (ii) exposed to localised, terrain-induced accelerations in mean airflow and moderate-to-strong wind gusts (> 9 m s −1 , based on Equation (1)).  Figure 8 gives the projected windthrow areas and associated maximum wind speeds for NB under current and future extreme wind conditions associated with two 56-year-and four 95-yearclassed extreme wind events for historical and future climate scenarios, assuming no net change in tree species distribution from current conditions. Windthrow under these projected wind events are expected to occur most strongly in the northcentral highlands and along the eastern coast of the province (red-coloured areas), where wind speeds are projected to be consistently stronger. Because high winds are common in mountainous terrain and along coastal areas due to atmospheric compression in elevated terrain and strong thermal gradients between land and water, many of the bF in these areas could conceivably become wind-firmed over their lifetime and, as a result, develop resistance to overturning in moderate-to-strong winds. Of course, the level of detail this presents is  distribution from current conditions. Windthrow under these projected wind events are expected to occur most strongly in the northcentral highlands and along the eastern coast of the province (red-coloured areas), where wind speeds are projected to be consistently stronger. Because high winds are common in mountainous terrain and along coastal areas due to atmospheric compression in elevated terrain and strong thermal gradients between land and water, many of the bF in these areas could conceivably become wind-firmed over their lifetime and, as a result, develop resistance to overturning in moderate-to-strong winds. Of course, the level of detail this presents is generally outside our capabilities to model with full confidence. The tendency for elevated windthrow in the NB highlands and regions closer to the coast did not change significantly with prevailing wind direction during the extreme wind events considered in this study. Maximum wind speeds in inland locations of NB are projected to exhibit minor changes over the next nine decades, with an average increase of about 0.037% per year compared to what is being predicted for coastal areas of the province, i.e., +0.066% per year. Although a temporary minor decline in wind speeds is projected for inland locations during the first 30 years under the RCP 8.5 climate scenario, small increases are projected for the remainder of the time.

Projected Windthrow under Current and Future Climatic Conditions
It is very unlikely that the projections illustrated in Figure 8 are reliable, especially when viewed at the individual stand level. Statistical downscaling of RCM-generated wind speeds may capture some macroscale features of the landscape (e.g., mountainous and near-coastal features) as shown in Figure 8, but the smaller-scaled features within these areas are generally too small to distinguish, especially considering that input wind speeds used in statistical downscaling ( Figure 6) are from computations performed at coarse spatial resolutions, i.e., ~50-km resolution, as compared to the 20m target resolution. CFD-based wind models can capture some of this spatial variability as demonstrated in Figure 7a, but its grid-area limitations, due to its exorbitant need of computer random access memory and processing time, render these types of models unsuitable for large DEMs at subhectare resolutions. As a compromise, we can view the probability of windthrow as a function of normalised ranking of statistically downscaled wind speeds, similar to the logistic density function addressed in Figure 7. Stands with a value of 1.0, indicating highest possible sustained maximum The tendency for elevated windthrow in the NB highlands and regions closer to the coast did not change significantly with prevailing wind direction during the extreme wind events considered in this study. Maximum wind speeds in inland locations of NB are projected to exhibit minor changes over the next nine decades, with an average increase of about 0.037% per year compared to what is being predicted for coastal areas of the province, i.e., +0.066% per year. Although a temporary minor decline in wind speeds is projected for inland locations during the first 30 years under the RCP 8.5 climate scenario, small increases are projected for the remainder of the time.
It is very unlikely that the projections illustrated in Figure 8 are reliable, especially when viewed at the individual stand level. Statistical downscaling of RCM-generated wind speeds may capture some macroscale features of the landscape (e.g., mountainous and near-coastal features) as shown in Figure 8, but the smaller-scaled features within these areas are generally too small to distinguish, especially considering that input wind speeds used in statistical downscaling ( Figure 6) are from computations performed at coarse spatial resolutions, i.e.,~50-km resolution, as compared to the 20-m target resolution. CFD-based wind models can capture some of this spatial variability as demonstrated in Figure 7a, but its grid-area limitations, due to its exorbitant need of computer random access memory and processing time, render these types of models unsuitable for large DEMs at subhectare resolutions. As a compromise, we can view the probability of windthrow as a function of normalised ranking of statistically downscaled wind speeds, similar to the logistic density function addressed in Figure 7. Stands with a value of 1.0, indicating highest possible sustained maximum wind speeds > 20 m s −1 with PWG's > 26 m s −1 , are expected to produce windthrow compared to stands with a lower probability of windthrow under lower wind speeds (captured in blue).
As bF habitat and bF retreats northward because of net increases in growing-season degree-days and increases in soil water content, the presence of bF and, therefore, windthrow in bF is anticipated to decline over the next nine decades (Figure 9). The rate of bF decline province-wide is expected to vary subject to differences in topography (low (warmer) vs. high-elevation areas (cooler)), growing degree-day sums and soil moisture. This decline is particularly pronounced for the RCP 8.5 climate scenario (associated with a net increase of 5.4 • C in summer over the 2081-2100 period, with respect to the 1986-2005 reference period), whereby 95.6% of viable bF habitat is expected to be lost by 2100 under RCP 8.5 ( Figure 9) compared to 50.9% under RCP 4.5. Increases in soil water content in already shallow soils, particularly in saturated-prone areas of the landscape (e.g., channel depressions), can further increase the risk of windthrow by restricting tree root growth (and anchoring potential) even more. The distribution of red and sugar maple habitat is expected to expand with a province-wide amelioration of growing conditions for the two species ( Figure 9). The replacement of bF with windfirm hardwood species, like the maples, has the potential to elicit substantial shifts in forest composition and wind firmness in future NB forests.
to decline over the next nine decades (Figure 9). The rate of bF decline province-wide is expected to vary subject to differences in topography (low (warmer) vs. high-elevation areas (cooler)), growing degree-day sums and soil moisture. This decline is particularly pronounced for the RCP 8.5 climate scenario (associated with a net increase of 5.4 °C in summer over the 2081-2100 period, with respect to the 1986-2005 reference period), whereby 95.6% of viable bF habitat is expected to be lost by 2100 under RCP 8.5 ( Figure 9) compared to 50.9% under RCP 4.5. Increases in soil water content in already shallow soils, particularly in saturated-prone areas of the landscape (e.g., channel depressions), can further increase the risk of windthrow by restricting tree root growth (and anchoring potential) even more. The distribution of red and sugar maple habitat is expected to expand with a province-wide amelioration of growing conditions for the two species ( Figure 9). The replacement of bF with windfirm hardwood species, like the maples, has the potential to elicit substantial shifts in forest composition and wind firmness in future NB forests. In situations where frozen soils provide additional anchoring for many of the trees, an incremental shift toward higher winter temperatures (net increase of 6.4°C in winter over the 2081-2100 period; Figure 10) may result in greater windthrow in winter, especially in areas of the province that are particularly susceptible to strong westerly and northwesterly winds common to the interior highlands of the province and northeastly winds for remote landward locations close to the eastern coast of NB ( Figure 6). With an abundance of snow on the ground in Atlantic Canada, soils in highly- In situations where frozen soils provide additional anchoring for many of the trees, an incremental shift toward higher winter temperatures (net increase of 6.4 • C in winter over the 2081-2100 period; Figure 10) may result in greater windthrow in winter, especially in areas of the province that are particularly susceptible to strong westerly and northwesterly winds common to the interior highlands of the province and northeastly winds for remote landward locations close to the eastern coast of NB ( Figure 6). With an abundance of snow on the ground in Atlantic Canada, soils in highly-to-moderately dense coniferous forests are seldom frozen to any great depths during the snow-accumulation period of the year (Figure 10) because of the insulating properties of both the overstorey vegetation and snow layer and the upward transfer of heat in seasonally cooling and snow-covered soils (Figure 10b). Consequently, coniferous forests, including bF forests in NB, rarely benefit from added anchoring in winter. Future lowering of the snowpack as winter air temperatures shift upward (Figure 10c,d) may cause midwinter thaw-refreeze events and cavitation in stemwood-porous species to occur, particularly in the maples and birches, which may predispose these species to dieback [43] and, thus, to stem-breakage and windthrow. This effect, however, will gradually dissipate in the future (e.g., 2070-2100) as air temperatures continue to climb and freezing of the soil complex is potentially confined to a very shallow depth of the forest floor ( Figure 10d). In their study, Saad et al. [44] found an overall increased risk of windthrow in bF forests in eastern Canada for the future (by 3% to 30%) mostly as a result of an increased duration of unfrozen soil conditions by the end of the 21st century under an RCP 8.5 climate scenario. Their windthrow projections for Atlantic Canada, however, do not account for the ongoing and projected erosion of bF growing conditions (habitat) in the region with continued climatic warming.
vegetation and snow layer and the upward transfer of heat in seasonally cooling and snow-covered soils (Figure 10b). Consequently, coniferous forests, including bF forests in NB, rarely benefit from added anchoring in winter. Future lowering of the snowpack as winter air temperatures shift upward (Figures 10c and 10d) may cause midwinter thaw-refreeze events and cavitation in stemwood-porous species to occur, particularly in the maples and birches, which may predispose these species to dieback [43] and, thus, to stem-breakage and windthrow. This effect, however, will gradually dissipate in the future (e.g., 2070-2100) as air temperatures continue to climb and freezing of the soil complex is potentially confined to a very shallow depth of the forest floor ( Figure 10d). In their study, Saad et al. [44] found an overall increased risk of windthrow in bF forests in eastern Canada for the future (by 3% to 30%) mostly as a result of an increased duration of unfrozen soil conditions by the end of the 21 st century under an RCP 8.5 climate scenario. Their windthrow projections for Atlantic Canada, however, do not account for the ongoing and projected erosion of bF growing conditions (habitat) in the region with continued climatic warming. Simulations with the ForHyM-model [45] were conducted with a soil complex involving a forest floor (7-cm thick), A-(11-cm), B-(50-cm), C1-(50-cm), C2-horizon (50-cm) and 11 subsoil layers, each 100cm thick. Soil temperature at each depth was calculated at the centre of each layer.

Conclusions
In this paper, we present a set of windthrow equations used to determine the presence/absence of windthrow in NB forests assessed from remote sensing following episodes of hurricane-strength  [45] were conducted with a soil complex involving a forest floor (7-cm thick), A-(11-cm), B-(50-cm), C1-(50-cm), C2-horizon (50-cm) and 11 subsoil layers, each 100-cm thick. Soil temperature at each depth was calculated at the centre of each layer.

Conclusions
In this paper, we present a set of windthrow equations used to determine the presence/absence of windthrow in NB forests assessed from remote sensing following episodes of hurricane-strength winds. The system integrates the binary calculation of windthrow to modelled biophysical surfaces of wind speed, forest state (density and size classes) and site soil conditions (soil depth, soil water content, slope) through genetic-programming-driven logistic regression. Results from a validation of the windthrow function for bF showed reasonable agreement between modelled and photointerpreted areas of windthrow in the Christmas Mountains following a major windstorm. From a wind firmness point of view, bF is the least resistant to windthrow when confronted with moderate-to-strong winds. This is illustrated by the relatively low critical wind speeds needed to initiate windthrow in bF stands, compared to other tree species, such as the maples.
Maximum wind speeds are expected to increase throughout NB over the next nine decades, but at nominal rates, i.e., 0.066 vs. 0.037% per year under climate-scenario RCP 8.5 for coastal and inland locations. Because coastal winds are generally stronger most of the time, due to strong thermal gradients between land and water, the net increase in maximum wind speeds by the end of the 21st-century is expected to be greatest for coastal areas of the province. Inland areas of NB will see much lower increases overall. The predominance of windthrow in currently existing bF stands in NB is expected to be moderated in the future, as the species responds to an overall decline in viable growing conditions with increased temperatures and soil moisture, and a gradual physiological and structural weakening of residual bF. Windfirm hardwood species are likely to replace bF as it progressively withdraws from NB, potentially causing future, transitioning forests in NB to be more resistant to windthrow.
Given the range of natural climate variability and uncertainties regarding future greenhouse-gas emission pathways and climate response, changes projected by one climate model, or one emission scenario, should not be used in isolation (e.g., [46]). A range of projections from multiple climate models (ensembles) and emission scenarios is therefore a way to consider (a good practice to identify the main sources of uncertainties), especially when regional or local climate variables as maximum wind speed and wind gusts are concerned. There is low confidence in the projection of regional storm track changes, and the impact of such changes on regional surface climate from global climate models (e.g., [46]). There is a clear need to assess more regional and local scale windstorm changes using high-resolution information (i.e., finer than 50 km) from RCM's. This is true for both extratropical and tropical cyclones that affect the maritime areas over the major part of the year. This will help confirm the outcomes of our study and provide improved risk assessments of windthrow in all tree species.

Appendix A (Follows below, in Landscape Format)
where ηcrit is the critical sustained wind speed (m s -1 ), below which windthrow is very unlikely to occur. The critical wind speed is determined numerically from windthrow and CFD-modelled wind data associated with extratropical cyclone Arthur, July 4-5, 2014.  (Table A2) denote stands predicted to be unaffected by windthrow, whereas p(ϕk)-values > threshold, stands predicted to have sustained some level of windthrow; 2 independent variables in windthrow-function definition: DC is the dominant-layer density class from 1-5, where 1=number of stems up to 600 (merchantable) or up to 5,000 (unmerchantable), 2=number of stems up to 601-1,200 (merchantable) or up to 5,001-10,000 (unmerchantable), 3=number of stems >1,200 (merchantable) or 10,001-20,000 (unmerchantable), 4=number of stems 20,001-30,000 (unmerchantable) and 5=number of stems > 30,000 (unmerchantable; NB DNR Data Dictionary, 2004); HNDP is the height above nearest drainage point (m); SD is the soil-depth class (same as SDEPTH in Tables 1 and 2 in the main body of the manuscript) from 1-4, where 1= < 30 cm, 2=30-65 cm, 3=65-100 cm and 4= >100 cm; SLP is the mean slope (SLP_MEAN) of the local terrain ( o ); SWC is the mean relative soil water content (SWC_MEAN) that ranges from 0-1, where 0 represents the permanent wilting point and 1, at field capacity; and WS is the mean wind speed (WNDSPD_MEAN) in m s -1 ; 3 δi's (i=2,…,4) are several variable functions, evaluated at exactly 1.0, when conditions in columns 3-5 are satisfied and 0.0, otherwise.

Species
where η crit is the critical sustained wind speed (m s −1 ), below which windthrow is very unlikely to occur. The critical wind speed is determined numerically from windthrow and CFD-modelled wind data associated with extratropical cyclone Arthur, July 4-5, 2014.  Table A2. Input records and results from windthrow-function generation with nonlinear, logistic regression. Results include function accuracy, critical sustained wind speeds, decision thresholds and number of observations used in function formulation (Table A1); η crit (species' critical windspeed requirement for windthrow) and decision thresholds are set during regression.