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

Abstract

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.

1. 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, early morning Christmas Mountains blowdown of November 7, 1994, where 3.6 × 10⁶ m³ of timber was lost to windthrow, valued at roughly $100 million in 1994 CAD [5]. On a landscape hierarchy, the degree of windthrow for a specific location arises from localised interactions with prevailing synoptic and small-scale airflow patterns (sustained and turbulent), variation in topography, storm duration and intensity, vegetation cover and on-the-ground management practices [6,7]. Because of the danger windthrow presents, it is important to identify and manage forest areas that are at significant risk of being damaged by strong winds [8].

The goal of management is to reduce risk in a system that is in state of constant flux [9,10]. The risk of windthrow affects not just ecosystem goods, such as maximum sustainable yield, but also ecosystem services and human social, health and economic systems [10]. This risk can be managed using traditional and modern forest management practices. Traditional forest management practices prescribe stand homogeneity [10]. This is done to optimize tree growth, wood production [10] and, thus, economic return. However, this purely economic outlook of forest management does not take into account natural disturbances of forests, nor does it take into account the impact of mechanical extraction of wood on the ecosystem itself [11]. Modern forestry, in contrast, focuses on resilience-based approaches [12]. Resilience-based management of forests is inspired by natural variability, with a particular focus to mimic natural disturbances to ensure ecosystem diversity in maintaining ecological integrity [10].

There is some disagreement over windthrow management solutions. In the boreal forests of northern Canada, the dominant natural disturbance is fire. Therefore, forest operations mimic this disturbance by clearcutting [13]. However, clearcutting is not sustainable in temperate mixed wood forests, like those of southern Quebec and NB, where the dominant natural disturbance is largely driven by defoliation by insects, such as by the spruce budworm (Choristoneura fumiferana Clem.), and windthrow [13]. Windthrow occurrences increase with the proportion of trees harvested, and salvage cutting from insect outbreaks can therefore have a compounding effect [14]. Windthrow and insect outbreaks change species composition and habitat [15] by creating a multicohort forest with an irregular stand structure [13].

This paper describes the intermediary role of high-velocity winds and forest-state variables in producing windthrow in balsam fir (bF; Abies balsamea (L.) Mill.) forests in NB (e.g., growing medium, forest structure and presence of windthrow assessed by remote sensing). The study involves the development of a species-specific windthrow function for bF with genetic programming (GP)-based logistic regression. We subsequently used the species-specific windthrow function to investigate the impact that future climate, based on regional climate model (RCM) simulations, may have on landscape-level windthrow impacts on bF forests. In providing comparisons to the landscape-level behaviour of bF during high-wind events, we generated a secondary set of windthrow functions for four other tree species, including for black spruce (Picea mariana (Mill.) B.S.P.), eastern white cedar (Thuja occidentalis L.) and red (Acer rubrum L) and sugar maple (Acer saccharum Marsh.). The performances of these species at the landscape level were not tested here due to the absence of independent data needed for this activity. However, the equations and their parameterisation provide insight as to the behaviour of these other tree species, in contrast to bF, during episodes of strong wind and species-specific factors that initiate windthrow in these species at the landscape level.

2. Materials and Methods

2.1. Study Areas

New Brunswick falls in the Atlantic Maritime Ecozone of eastern Canada (i.e., 44.59939° to 48.02186°N latitudes and 63.78114° to 69.05178°W longitudes; Figure 1), and is characterised by a fragmented landscape of predominantly forests and agriculture and variable terrain (i.e., elevations range from 0 to 834 m above mean sea level, AMSL). Forests occupy about 85% of the province’s land base, and agriculture occupies about 6%. Provincial climate is largely influenced by the region’s proximity to the Bay of Fundy in the south, the Northumberland Strait in the east and the Bay of Chaleur in the upper northeast corner of the province. On average, winters are milder than most regions of Canada. Springs are usually late and cool, and summers are moderate and breezy. The province experiences, on average, an annual mean air temperature and total precipitation between 2 °C and 7 °C and 1000 mm and 1360 mm, respectively [16].

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.

2.2. 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;
• Spatially explicit reconstructions of wind fields during the peak of cyclone Arthur obtained with a computational fluid dynamics model (CFD) and input wind data (speed and direction) from two Environment and Climate Change Canada (ECCC) weather stations in southwestern NB (i.e., Fredericton Airport and St. Stephen weather stations; Figure 1);
• 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. Suite of 23-predictor variables prior dimensionality reduction with principal component analysis (PCA). WNDSPD (3), SLP (3), SWC (3) and HNDP (3) took into account the minimum, mean and maximum values of the variable for a given stand and were treated as independent (predictor) variables.

  Grouping	Variable	Significance	Data Source
  Wind-related variable	WNDSPD (3)	Wind speed (in m s−1); meteorological agent needed for windthrow	CFD-calculations based on input data relevant to the wind event
  Forest-state and tree-related variables	TREEFORM	Relates to the sail area of the crown; grouped according to dominant stand type, i.e., hardwood or softwood (1 or 2)	GIS-thematic forest data, based on dominant first species code (L1S1 1)
  	DEVS	Dominant tree layer development stage class (1–6)	GIS-thematic forest data (L1DS 1)
  	FRAC	Dominant tree layer first species % ratio (2–10)	GIS-thematic forest data (L1PR1 1)
  	CC	Dominant tree layer crown closure code (1–5)	GIS-thematic forest data (L1CC 1)
  	DC	Dominant tree layer density class (1–5)	GIS-thematic forest data (L1DC 1)
  	SC	Dominant tree layer size class (1–3)	GIS-thematic forest data (L1SC 1)
  Terrain- and soil moisture-related variables	SLP (3)	Stand slope (in degrees); windthrow on slopes varies depending on the prevailing airflow and secondary air circulation	DEM-based calculation, based on finite differencing
  	SWC (3)	Stand relative soil water content (unitless); high soil water content tends to constrain root development and, as a result, soils of high SWC can encourage windthrow	DEM-based water-budget calculation with LanDSET 2
  	HNDP (3)	Terrain height above nearest drainage point in metres; functions similar to SWC in defining windthrow especially for wet areas in the landscape	DEM-based height difference with minimising-distance function to locate nearest drainage point
  Soil-related variables	SFERT	Soil fertility class (1–4); potentially defines growing potential of the soil and species wind firmness	GIS-thematic soil data, derived from soil parent material nutrient content and degree of weatherability
  	SDEPTH	Depth to contrasting layer (1–4); shallow soils function to restrict vertical root growth and reduces the anchoring potential of the soil	GIS-thematic soil data 3
  	CONTLYR	Contrasting layer description (1–7)	GIS-thematic soil data 3
  	STEXT	Soil texture class (1–3); addresses mechanical and water-drainage properties of the soil, i.e., anchoring potential	GIS-thematic soil data 3
  	SDRAIN	Soil drainage class (1–7); a field indicator of soil drainage taking into account soil texture and topographic position	GIS-thematic soil data 3

  ¹ NB DNR Data Dictionary, 2004; ² Landscape Distribution of Soil Moisture, Energy and Temperature [19,20]; ³ after Fahmy et al. [21].

  ), 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⁻² 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.

2.3.1. Digital Elevation Model

The DEM of NB was based on the 30-m resolution Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model ver. 2.0 (GDEM; http://asterweb.jpl.nasa.gov/gdem.asp, accessed on March, 2015) resampled at a 20-m resolution with bicubic interpolation (Figure 1). The DEM served as a data layer central to modelling the likelihood of windthrow province-wide under six extreme wind events.

2.3.2. 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⁻¹ and 16.3 m s⁻¹ 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⁻¹ 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.

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.,

$$\mathrm{PWG}’\mathrm{s}=\left[1.9506-0.173\cdot \ln \left(\frac{\eta }{0.44704}\right)\right]\cdot \eta ,$$

(1)

where η is the sustained maximum wind speed (m s⁻¹). The denominator in the quotient in Equation (1), i.e., 0.44704, allows for a unit conversion from m s⁻¹ to miles hour⁻¹. The bracketed term in Equation (1) is the gust factor [29].

2.3.3. 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].

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].

2.3.4. 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.

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.

2.3.5. 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).

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.

2.3.6. 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⁻² 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⁻² 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 (1950–2005; 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⁻¹, attaining wind gusts > 13 m s⁻¹, particularly during winter, due to stronger atmospheric circulation.

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.

2.3.7. 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.

3. Results and Discussion

3.1. Principal Component Analysis

Table 2. Factor loadings. Values in bold are the factor loadings (large positive or negative loadings) used in grouping variables to individual principle components (PCs). Percentages given are the proportion of total variance explained by individual PCs. Variables in bold are used as independent variables in the development of the windthrow function for bF and the other tree species.

Variable 1	PC 1 (14.2%)	PC 2 (13.5%)	PC 3 (8.8%)	PC 4 (11.7%)	PC 5 (7.7%)	PC 6 (8.4%)	PC 7 (8.0%)	PC 8 (6.7%)
WNDSPD_MIN	0.026	0.954	−0.039	0.068	−0.039	0.021	−0.112	0.102
WNDSPD_MAX	0.019	0.940	−0.022	0.111	−0.051	−0.059	0.174	−0.144
WNDSPD_MEAN	0.024	0.981	−0.038	0.097	−0.043	−0.033	0.052	−0.037
TREEFORM	−0.119	−0.061	0.078	0.073	0.901	0.006	−0.047	0.072
DEVS	−0.025	0.023	0.283	0.031	−0.843	0.137	0.024	−0.022
FRAC	0.114	−0.009	0.010	−0.054	0.430	0.194	0.088	−0.077
CC	−0.143	0.096	−0.805	0.004	−0.040	−0.132	0.039	0.125
DC	−0.054	0.023	−0.931	−0.063	0.028	0.002	−0.015	0.005
SC	−0.208	0.034	0.644	0.021	−0.179	−0.135	0.056	0.166
SLP_MIN	0.209	−0.040	0.002	−0.011	0.005	0.049	0.106	0.912
SLP_MAX	0.142	0.033	0.015	0.017	0.017	−0.172	0.937	−0.047
SLP_MEANb	0.285	0.026	0.014	0.003	0.024	−0.099	0.7662	0.5252
SWC_MIN	−0.190	−0.111	0.008	0.014	0.016	0.798	−0.340	0.216
SWC_MAX	−0.520	0.026	−0.026	0.151	0.050	0.541	0.158	−0.451
SWC_MEAN	−0.413	−0.060	0.007	0.097	0.058	0.836	−0.116	−0.122
HNDP_MIN	0.880	0.071	−0.025	−0.082	0.031	−0.085	−0.004	0.193
HNDP_MAX	0.853	−0.083	0.049	−0.011	0.018	−0.219	0.290	−0.003
HNDP_MEAN	0.929	−0.010	0.009	−0.048	0.022	−0.178	0.168	0.102
SFERT	0.334	0.508	−0.003	−0.438	−0.075	0.111	0.095	−0.074
SDEPTH	−0.075	0.000	−0.013	−0.867	0.089	−0.079	−0.008	0.088
CONTLYR	0.078	−0.185	0.026	−0.9223	0.039	−0.022	0.040	0.042
STEXT	0.238	−0.064	−0.066	−0.815	−0.033	0.070	0.034	−0.081
SDRAIN	0.104	0.116	0.049	0.396	0.083	0.249	0.174	−0.000

¹ For variable definition, consult Table 1 in the main body of the manuscript; ² SLP_MEAN had a strong association with both PC’s 7 and 8 (columns 8 and 9) and was easiest to estimate spatially. As a result, we used SLP_MEAN as the only descriptor of slope (i.e., SLP_MIN and SLP_MAX were not considered further); ³ The description of the soil contrasting layer was not very meaningful in comprehending the role of soils on windthrow. As a result, we used SDEPTH instead.

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.

3.2. Windthrow Function Development

Windthrow equation for bF and associated thresholds (δ_i, i=1,…,4 and decision thresholds) were derived from GP-based regression (

Table A1. Windthrow functions, p(ϕ_k), for five native tree species of southwest New Brunswick based; ${\mathsf{\delta }}_1=\{\begin{array}{l}1,\hspace{1em}\mathrm{WS}>{\mathsf{\eta }}_{\mathrm{crit}}\\ 0,\hspace{1em}\mathrm{otherwise}\end{array},$ where η_crit is the critical sustained wind speed (m s⁻¹), 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.

Species	p(ϕk)1,2	δ23	δ3	δ4
Balsam Fir	$\frac{{\delta }_1}{1.0+0.002\exp \left\{9.2\cdot \left[\left(\frac{\sin \left(69.0\cdot \mathrm{HNDP}\right)-9.2\cdot {\mathrm{SWC}}^2}{\mathrm{SLP}-23.4}\right)+{\delta }_2\cdot {\delta }_3-{\delta }_4\right]\right\}}$	$\mathrm{SWC}<0.7$	$\begin{array}{l}\left[\begin{array}{l}\mathrm{HNDP}+9.2\cdot {\mathrm{SWC}}^2\\ +\cos \left(69.0\cdot \mathrm{HNDP}\right)\end{array}\right]\\ <\mathrm{SLP}\end{array}$	$\begin{array}{l}4.3-\max \left(\begin{array}{l}\mathrm{SWC},(\mathrm{SLP}\\ -23.4)\end{array}\right)\\ <\mathrm{SD}\end{array}$
Black Spruce	$\frac{{\delta }_1}{1.0+\exp \left\{-829-1.5\times {10}^3{\delta }_2-\min \left[\frac{\mathrm{DC}}{\min \left(0.3\cdot \mathrm{SLP},\mathrm{HNDP}\right)},\cos \left[\cos \left(0.3\cdot \mathrm{SLP}\right)\right]\right]\right\}},$	$\begin{array}{l}\min \left(\mathrm{SLP}-7.6,8\right)\\ <\mathrm{SWC}\cdot \left(\mathrm{SLP}-7.6\right)\end{array}$
Eastern White Cedar	$\frac{{\delta }_1}{1.0+\exp \left\{385-395\cdot \max \left[\sin \left(\mathrm{HNDP}\left(-47.8-2.1\cdot \mathrm{HNDP}\right)\right),\cos \left(-0.3\cdot \mathrm{HNDP}\cdot \left(\mathrm{SLP}-4.3\cdot \mathrm{HNDP}\right)\right)\right]\right\}}$
Red Maple	$\frac{{\delta }_1}{1.0+\exp \left\{81.0-120.3\cdot \sin \left[\mathrm{HNDP}\cdot \min \left(14.9,\max \left(8.2,\mathrm{SLP}\right)\right)-{\delta }_2+\mathrm{HNDP}+\sin \left(1.5\cdot \mathrm{SLP}\right)\cdot {\delta }_3\right]\right\}}$	$\mathrm{SWC}<0.5$	$\mathrm{HNDP}<(\mathrm{DC}-1.1)$
Sugar Maple	$\frac{{\delta }_1}{1.0+\exp \left\{\frac{430\cdot \left(\mathrm{SWC}+\mathrm{HNDP}\right)}{10.4+\sin \left(5.9+\mathrm{SWC}+\mathrm{HNDP}\right)-\sin \left(\mathrm{SD}\right)\cdot {\delta }_2}-751\right\}}$	$\begin{array}{l}{\left(0.3+\mathrm{HNDP}\right)}^{\sin \left(\mathrm{HNDP}\right)}\\ <\mathrm{DC}\end{array}$

¹ p(ϕ_k)-values < threshold (Table A2) denote stands predicted to be unaffected by windthrow, whereas p(ϕ_k)-values > threshold, stands predicted to have sustained some level of windthrow; ² 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 5000 (unmerchantable), 2 = number of stems up to 601–1200 (merchantable) or up to 5001–10,000 (unmerchantable), 3 = number of stems >1200 (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 Table 1 and Table 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 (°); 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⁻¹; ³ δ_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.

and

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.

Tree Species	Accuracy 1 (%)	ηcrit (m s−1)	Peak Wind Gust Associated with ηcrit 2 (m s−1)	Threshold (unitless)	Number of Observations
Balsam Fir	94.6	5.8	8.8	0.03	518
Black Spruce	~100	6.4	9.6	0.49	129
Eastern White Cedar	93.4	5.8	8.8	0.02	332
Red Maple	96.3	7.0	10.3	0.02	243
Sugar Maple	~100	9.5	13.5	0.50	112

¹ % accuracy is defined as the area under the receiver operating characteristic curve; ² calculations of peak wind gusts are based on eqn. (1) in the main text.

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⁻¹ compared to 7.0 m s⁻¹ and 9.5 m s⁻¹ for red and sugar maple (Table A2). Estimated PWG’s of 8.8 m s⁻¹, 10.3 m s⁻¹ and 13.5 m s⁻¹ were associated with the calculated critical wind speeds of 5.8 m s⁻¹, 7.0 m s⁻¹ and 9.5 m s⁻¹. 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⁻¹, 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).

3.3. 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 and June 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⁻¹), and (ii) exposed to localised, terrain-induced accelerations in mean airflow and moderate-to-strong wind gusts (> 9 m s⁻¹, based on Equation (1)).

3.4. Projected Windthrow under Current and Future Climatic Conditions

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-year-classed 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 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.

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⁻¹ with PWG’s > 26 m s⁻¹, 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.

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.

4. 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.