Mountain Pine Beetles, Salvage Logging, and Hydrologic Change: Predicting Wet Ground Areas

The mountain pine beetle epidemic in British Columbia has covered 18.1 million hectares of forest land showing the potential for exceptionally large-scale disturbance to influence watershed hydrology. Pine stands killed by the epidemic can experience reduced levels of evapotranspiration and precipitation interception, which can translate into an increase in soil moisture as observed by some forest practitioners during salvage logging in the epicenter of the outbreak. They reported the replacement of summer ground, dry firm soil areas, with winter ground areas identified by having wetter, less firm soils upon which forestry equipment operation is difficult or impossible before winter freeze-up. To decrease the likelihood of soil disturbance from harvesting, a set of hazard indicators was developed to predict wet ground areas in areas heavily infested by the mountain pine beetle. Hazard indicators were based on available GIS data, aerial photographs, and local knowledge. Indicators were selected by an iterative process that began with office-based selection of potential indicators, model development and prediction, field verification, and model refinement to select those indicators that explained most field data variability. Findings indicate that the most effective indicators were lodgepole pine content, understory, drainage density, soil texture, and the topographic index.


Introduction
The mountain pine beetle (MPB) epidemic in British Columbia is the most severe bark beetle infestation recorded in the history of North America.Its origin is the result of climate change, specifically winter minimum temperatures that are too warm to kill beetle larvae in combination with effective forest fire suppression that maintained an abundant mature pine tree population [1,2].The beetle infestation increased dramatically from 2000 to 2006 progressing from an endemic outbreak to an epidemic that has presently affected an area of some 18.1 million hectares.The scale of this infestation is a significant challenge for forestry and watershed management.
The initial response to the infestation was to increase the allowable cut within the most heavily infested areas where the outbreak began to recover the economic value of attacked pine stands and to expedite the regeneration of new forests in those areas.Salvage logging operations reported difficulties due to a loss in summer ground between 2003 and 2005.That is, during those harvest years forest practitioners encountered wet soils that made equipment operation difficult or impossible where they had expected to encounter dry soils capable of supporting heavy equipment.The frequency of these observations over the landscape signaled a possible change in water balance and subsequent ecology of affected areas.
Mountain pine beetles (Dendroctonous pondorosae, Hopkins) can affect forest hydrology by killing pine tree stands.MPB generally burrow into pine trees in the lower bole where they create galleries to mate and lay their eggs.During this process they pass through the phloem and once eggs hatch in the early summer and fall, larvae feed on the phloem, which increases tree stress and lowers transpiration [3].When a tree is attacked by a sufficient number of beetles it dies.Once dead, the tree passes from green to red attack where the pine needles first change from green to red the first year post-infestation, after which they progress to grey attack where trees lose their needles and then fine branches over the next couple of years.During this transition, dead and dying pine tree stands will have lower evapotranspiration rates than stands of live trees [4][5][6].The loss of pine needles and branches decreases precipitation interception, which can vary between 15% and 35% among coniferous forest stands depending upon precipitation amount and form, tree species, and stand characteristics [7][8][9].Transpiration also decreases or ceases, which will also increase soil moisture.The loss of transpiration can be substantial, for example lodegepole pine stand transpiration accounted for 50% and 61% of total evapotranspiration (ET) in pine stands of southeastern Wyoming [10].Although the loss in transpiration may account for lower amounts of water removed from the soil, interception is still considered the important factor accounting for "watering-up", a term used to identify an increase in groundwater table elevation following harvesting [11].As MPB affected trees progress to grey attack they drop needles decreasing the interception of precipitation up to 50% above pre-infestation levels [12].This exceeds the reduction of evapotranspiration between six and 39% observed from sub-boreal watersheds subjected to harvesting alone [13].Increased delivery of precipitation to the ground (i.e., net precipitation) may be stored or moved through the watershed unless remaining live trees or an understory transpire sufficiently as noted by Brown [14].Where increased net precipitation is stored, water table elevation can increase and soils may "water-up" [15].If it moves through the watershed and is exported, the annual water yield will increase [16,17].
The project discussed here developed a watershed level hydrologic hazard assessment procedure to assess the relative hazard of experiencing wet soils during salvage logging due to increased water table elevation and/or delayed surface drainage resulting from watershed characteristics and the MPB infestation.Wet ground during salvage operations increases the likelihood of soil disturbance, which can decrease productivity, increase silviculture costs, and increase surface erosion that can lead to water quality problems for streams receiving run-off.Field observations were used to assess effectiveness of the hazard assessment process as well as to provide insight on the cumulative effect of the MPB infestation and salvage logging on soil hydrology.This project is also distinctive because it examines logging operations under a set of ecological conditions related to massive insect outbreak whereas works most often cited in the literature are the effects of salvage operations following wildfire [18].

Study Area and Climate Analysis
The study focused on watersheds in the Vanderhoof Forest District of British Columbia (Figure 1).The epidemic began in this forest district and it was particularly affected by the outbreak because forests of the area are more than 80% composed of lodgepole pine (Pinus contorta var.latifolia).To predict the hazard of wet ground within third and fourth order watersheds, each watershed in the district was systematically evaluated using the same indicators to assess how likely it would be to have wet ground relative to other watersheds.Due to the absence of hydrologic information at the scale of our investigation, two hazard assessment approaches were taken.The a priori approach predicted watershed hazard based upon indicators selected from the hydrologic literature as well as professional opinion.The post-hoc approach consisted of an exploratory statistical review of field data to identify indicators most effective at explaining field data variability.In addition to hazard analysis a climate analysis was completed to identify change in precipitation timing and quantity.
To assess climatic trend influence on field observations a review of climate information and an assessment of climatic trends were conducted based on five different Meteorological Services Canada (MSC) weather stations near Vanderhoof.This analysis focuses on data collected between 1980 and 2007 because the weather stations were most similar, being at the same elevation and only 2.2 km apart.Further, this recent period of data provides a relevant climatic context to concerns about the loss of summer ground.Annual and seasonal trend analysis for precipitation was assessed using a t-test of slope.

The a Priori Approach
The a priori approach used two categories of hazard indicators, specifically those that infer the potential for increased net precipitation and those that infer the potential for retention of increased net precipitation in the soil (Table 1).Soil and landform maps of the Vanderhoof Forest District (1:50,000) were coded to identify fine texture soil types (ex.lacustrine) prone to shallow or perched water tables along with organic soils.Fine surface and organic soils will likely have higher ambient soil moisture conditions and will respond more to increased net precipitation than coarser welldrained soil types.Understory Score Multi-storied stands may not see an increase in net precipitation after MPB infestation [14,19].Understory regeneration increases in the sequence of sub-boreal spruce dry cold, dry warm, and moist cool, as well as Engelmann Spruce Subalpine Fir zone to more than 1000 stems/ha [20], which reduces net precipitation.

Drainage Density
Provides an estimate of how efficiently water leaves an area during storms [21] by providing an index of relative distance between where rain falls and flowing channels [22].
Net precipitation indicators targeted forest condition and beetle infestation characteristics such as available rearing habitat (mature pine content), amount of grey attack stage pine trees and the infestation severity data from the 2004 aerial overview survey.Retention indicators focused on watershed characteristics influencing snowmelt and runoff conditions such as understory condition, drainage density, as well as soil texture and moisture.Aspect was excluded because the topography of the area was relatively flat and did not allow adequate watershed differentiation.
The a priori hazard score was calculated as [Equation ( 1)]: Hazard = potential for increased net precipitation × potential for retention of increased net precipitation (1) Scores were normalized to a scale of 100 and hazard categories were assessed as low for the first quartile (0-25), moderate for the second quartile (25-50), and high for the third and fourth quartile (50-100).

The Post-Hoc Approach
Hazard indicators were identified in the post-hoc approach using a coarse and fine filtering approach [23].The coarse filtering approach consisted of a principle components analysis (PCA) of field data from seven of the 17 sites to identify groups of correlated hazard indicators that explained a high proportion of data variability [24].Hazard indicators include those identified for the a priori approach as well as the topographic index [25,26] and relief ratio [24].Once groups of indicators were identified, they were fine filtered using a stepwise general linear model, which identified those indicators that had the highest predictive power for identifying field verified "wet sites".Pos-hoc model verification was completed by comparing predicted and observed conditions at the ten sites not used for model-building

Site Classification and Study Design
A total of 17 watersheds were chosen for study in 2005-2007 using the first hazard assessment results.These include six low hazard, four moderate hazard, and seven high hazard watersheds as identified by the a priori method (Table 2).
The prominence of high and low hazard watersheds was intentional to maximize the amount of data gathered to identify field differences between these classifications.Study sites were in pine dominated stands with dry to average soil moisture conditions and they were located along hill slopes in lower watershed reaches to ensure similar sampling environments between watersheds.Within these 17 watersheds, seven detailed assessment study areas were established in 2005 while the remaining ten qualitative assessment areas were chosen in 2005 and 2006 (Table 2).Field information from 2006 and 2007 was used for the a priori and post-hoc assessment approach because those sampling seasons provided the most consistent datasets across all the detailed and qualitative assessment sites.Volumetric soil moisture content (θ) was measured in each watershed at a grey attack stand as well as a bordering harvested area along the same slope at the toe, mid-slope and summit positions at 10, 20, 40, and 60 cm depth.Average soil moisture from the four depths was measured in the late summer and fall, when summer ground issues should be observed, using a Thetaprobe connected to a hand-held reader (Delta-T Devices Ltd., Cambridge, UK).Soil-specific calibration of the Thetaprobe followed the general procedure for calibrating capacitance sensors and provided accuracy of ±1%-3% for our soils [27,28].Volumetric soil moisture content above 30% was considered to be at levels where harvesting equipment may cause soil disturbance and concurrently experience operational difficulties.

Detailed Assessment Approach
Detailed assessment sites have nine wells located along MPB killed forested pine-leading stands and harvested slopes (Figure 2).Specifically, there are three wells installed along a transect at the level summit, middle slope and toe slope of the harvested and forested sites to study the range of variability in soil hydrologic properties [9,29].Soil structural-textural conditions were confirmed in proximal pedons to ensure within-site characteristics are as uniform as possible.
Field measurements were gathered at two to three week intervals in spring and summer as well as early fall.These periods were chosen to examine seasonal water table fluctuations.This sampling frequency also allowed for observation of surface ponding, soil saturation, and surface flow due to precipitation as measured by a proximal rain gauge.Shallow wells (<1 m) were excavated by auger and lined with a 4 cm interior diameter PVC pipe [30].Water table depths were measured at each site using a dipper or electrical buzzer probe [31].Water table depths less than 60 cm below surface were considered "shallow" and to be an operational concern because the capillary fringe may be as high as 30 cm above the water table indicating a reduced amount of soil available for storage during rainfall events.The amplitude of the capillary fringe was determined by steel rods driven vertically in mineral soil [32] weekly over four growing seasons under similar soil conditions in another as yet unreported study.
To examine potential changes as to how water moves under saturated conditions, field saturated hydraulic conductivity (Kfs) was measured in five watersheds using a simplified falling-head technique [33].A value of α*-parameter = 12 m −1 was used to calculate Kfs [34].The five watersheds chosen cover the range of hydrological risks i.e. low, moderate, and high hazard systems as determined by the a priori approach.Kfs data was measured in both forested areas and harvested areas at the summit position to assess harvesting effect (i.e., soil disturbance) on soil drainage.Six randomly chosen 24 m 2 grid sampling areas were located within each of the forested and harvested areas.Within each grid, samples were drawn from 12 points.
As soil texture influences hydraulic conductivities, surface soil particle size of the fine soil fraction (<2 mm) was measured at each Kfs sampling point (Table 3) using the hydrometer method [35].Soils were predominantly coarse and relatively uniform in texture within the top layer (0-10 cm depth) across the sites except at Belisle Creek, which had high clay content (Table 3).Within the sites, forested and clear-cut areas had similar particle size distribution.
Bulk density (Db) information was also gathered to provide an index of soil compaction [36].Db was calculated on wet volume basis and determined by the core technique [37]), two cores (4 cm long × 5 cm diameter) collected at the mid-point of the 10 and 20 cm depths.

Statistical Analysis
Water table data collected from well sites were log-normally distributed and were transformed prior to Analysis of Variance (ANOVA) and used to identify differences in water table depth across slope location, hazard class, treatment (MPB and cutblock), season, and year using SYSTAT 11 ® software.Soil moisture measurements were subjected to similar analyses but did not require transformation.ANOVA for soil moisture data was completed using average soil moisture collected from the 10, 20, 40, and 60 cm depths.Significance for all tests was determined at a level of 0.05.
The Kfs values were also log-normally distributed and were transformed accordingly prior to a group t test (PROC TTEST) to compare transformed Kfs and normally distributed soil bulk density data (Db) within sites (MPB vs. clear-cut).All analyses were based on a significance level of p = 0.05.

Climate
Analysis identified that annual precipitation has increased over historical levels but not significantly (Table 4).The summers of 2005 and 2007 were the wettest on record.The ratio of summer to winter precipitation has increased significantly since 1997 (Table 4).Between 2001 and 2003, summer precipitation totals (June-September) were within 4% of the 1971-2000 normal of 191 mm.Summer precipitation levels for 2004 and 2005, and 2007 were considerably greater, ranging from 250 to 329 mm or 30%-75% higher than the 30-year normal.Given that summer months are generally wetter now than they were in the earlier period of the climatic record, there may be an effect on trafficability during harvesting activities in areas prone to poor drainage.

Field Findings
There was a significant slope location and seasonal effect (F 2,189 = 55.9, p < 0.001 and F 2,189 = 5.9, p = 0.003 respectively) on depth to water table across all sites.Toe slope locations had shallower water table levels than the other slope positions (Figure 3a) and were most often above the 60 cm threshold used to identify wet locations in 2006 and 2007 combined.Summer rainfalls were higher in 2007 (309 mm) than in 2006 (154 mm), however water table trends were similar between years.As expected, spring months had shallower water table levels than the summer months (Figure 3b).
The variability of Kfs values was high within and across sites, with coefficients of variation (CV's) ranging between 62% and 161% (Table 5).Kfs typically exhibits large spatial variability to which texture (e.g., highest Kfs in Belisle Creek-silty clay loam) and structure of soils is most directly related [38,39].Average Kfs results were particularly high overall compared to K values published elsewhere [33,40] near the low end of saturated hydraulic range for coarse-textured soils [39].This can be attributed to an overestimation of the α*-parameter [41].Minimum K results appear more like actual values (Table 5).The abundant presence of silt in the sand matrix may disperse and clog up the conductive pores upon wetting.There was no clear salvage logging effect on drainage patterns.Saturated soil infiltration did not show consistently lower values in cutblock areas (Table 5), possibly due to careful logging practices and large sampling error.Where compaction was evident such as at Targe Creek (i.e., lower Kfs and higher Db between 0 and 5 cm depth in the clear-cut) soil disturbance from harvesting led to poor drainage.Interestingly, the finer-texture soils at Belisle Creek had the highest rates of Kfs, which cannot be explained by texture alone since lower rates would have been expected (Table 3).At this site, the surface horizon had a loose crumb structure, which provided a very porous medium.Some differences in volumetric soil moisture were observed at the qualitative assessment sites across years, treatments, slope locations, and seasons (Figure 4).Not all differences were statistically significant.Although there appears to be differences between treatments, it is not statistically significant (F 1,186 = 3.86, p = 0.05).There were differences in seasons and locations (F 2,186 = 11.14, p < 0.001 , F 1,186 = 952.05,p < 0.001, respectively) with both the spring season and toe location having higher soil moisture than other seasons or slope positions.There was also a significant difference between years, with 2007 having generally higher soil moisture than 2006.This was expected because 2007 received more rainfall than 2006.

Post-Hoc Model Evaluation
The principal components analysis of the water table and average volumetric soil moisture content identified two groups of indicators that explained 90% and 80% of field data variability, respectively (Table 6).The general linear model (GLM) for water table and soil moisture data respectively refined this list of indicators, identifying lodgepole pine content, understory, drainage density, sensitive soils, and the topographic index as the most significant indicators.Although each GLM analysis provided an equation to predict specific values for water table or soil moisture, these formulae are not presented here because water table elevations and soil moisture cannot be predicted at the watershed scale.Instead, the equations were used to develop a new hazard prediction formula based upon the coefficient's scale and sign (i.e., positively or negatively correlated to depth to water table or soil moisture) for each indicator.Hazard rankings were considered correct when high hazard sites were wet in both the forest and cutblock locations, moderate sites could be wet in the cutblock due to the loss in transpiration, and low sites were dry in both locations.Table 6.Hazard indicators that were most effective at explaining data variability as identified by the principal components analysis.Note the same indicators were identified for both measurements.
This formula is more hydrologically relevant than that presented for the a priori approach because it emphasizes watershed characteristics that have direct influence on net precipitation and its retention such as understory, soil type, and the relative slope of the watershed.For example, understory can lower the increase in net precipitation and also transpire [14,19], areas with less sensitive soils may have better drainage than those with sensitive soils, and the area based slope of the watershed indicates retention time of water on the soil surface [22].
Scores generated by the post-hoc formula were then ranked from one to 100 with ties receiving the same rank (i.e., 50, 51, 51, 52 were ranked 50, 51.5, 51.5, 53).High hazard sites were those with the upper 25th percentile of ranked scores (i.e., , moderate hazard watersheds were the middle 50% , and the low hazard watersheds were the lower 25th percentile of scores (75-100).In contrast to the a priori approach, the post-hoc assessment correctly identified all sites (Table 7).
High hazard watersheds had significantly shallower depth to water table at the summit across years (Figure 5, F 2,57 = 5.61, p = 0.006).Harvesting effects on water table depth were not detectable as dead and dying pine stands had low transpiration due to dead pine trees as well as increased water delivery to soil more comparable to cutblock areas than to non-infested stands at toe and summit locations.High hazard sites had shallow water tables that were on average 25 cm closer to the soil surface than moderate and low hazard sites (Figure 5).Mid-slope water table was not affected by risk, season or treatment because midslope drainage is mostly controlled by gravity.All toe slope sites were wet regardless of site condition and hazard.Toe slopes had shallower water table levels in spring compared to summer whereas water table levels were the same across seasons for mid-slope and summit positions.The deepest water table was recorded in the fall suggesting effects of spring runoff on toe-slope receiving areas diminished during the summer until fall precipitation replenishes the soil moisture (Figure 6).The analysis of ln-transformed Kfs from all sites indicated a statistically significant effect on watershed hazard (F 3,256 = 4.10, p = 0.007) and site condition (F 5,254 = 3.71, p = 0.003).Highest measured hydraulic conductivities were found in the high hazard forested sites.In contrast, moderate hazard clear-cut area had the lowest Kfs but soil disturbance was not observed.Great variability in Kfs in part due to very high spatial variability in soil properties may explain this result.Based on our sampling, harvesting did not lead to a significant reduction in Kfs across watershed hazard.There was faster surface drainage in the high hazard MPB areas than in the low hazard MPB areas.This indicates that differences in Kfs may not be explained by high water table levels, which are less likely to occur where surface drainage is fast.Hard almost cemented layers less than 60 cm deep were observed during soil pit excavation at some high and moderate hazard sites (e.g., Targe Creek, watersheds 10557 and 10426) that may impede drainage similar to that observed in Ortstein layers [42].Although not impervious to water, the naturally compacted layer has a slower percolation rate that may be inadequate to drain large quantities of water reaching the soil in stands with a dead pine overstory or large salvage harvested areas.Under these conditions, soil saturation persists longer after spring runoff and a large summer storm can quickly fill up available storage in the soil profile raising the water table quickly, which may impede forest management activities.
The influence of pre-existing conditions in the soil profile such as a moist and soft layer lying over a dry or hard subsurface layer can be exacerbated by compaction [43] and may result in a higher hazard for salvage-logged areas.For example, there was a statistically significant relationship between site condition and Kfs (p < 0.01) at Targe Creek (Table 5).Compaction was evident in the clear-cut sampling areas and sample sites showed a significant increase in bulk density in the top 5 cm of soil following skidding (Table 5).These compacted areas were characterized by a platy structure and loss of original structure [44].The reduction in large pore space in the clear-cut, which is responsible for most of the saturated flow, produced an average Kfs rate of 115 mm h −1 ), which represented a reduction of 57% in Kfs from the MPB forest.Harvesting operations and subsequent soil compaction can decrease field saturated hydraulic conductivity [45].

Conclusions
This project developed a field-verified watershed level hydrologic hazard assessment procedure to assess the relative hazard of experiencing wet soils during salvage logging of MPB-affected forests in the Vanderhoof Forest District.Watering-up occurred in response to precipitation in watersheds with characteristics that increase net precipitation and retention of that precipitation such as dominant forest stand cover of dead pine trees, lack of understory, low watershed slope, and fine textured soils.Findings indicate that salvage logging did not influence water table elevations or soil moisture when compared to standing beetle-killed stands but it did increase soil compaction, which can alter drainage pattern and efficiency.Although the model presented here was developed for the Vanderhoof Forest District, some model components may be transferrable to other areas experiencing a forest pest outbreak or other watershed-level disturbance.Similarly, the model-development process presented here is suitable for transfer to other areas requiring a watershed-level hazard analysis for increase in soil moisture.

Figure 1 .
Figure 1.Location of the Vanderhoof Forest District in British Columbia Canada.

Figure 2 .
Figure 2. Orthophoto image of the Belisle Creek detailed assessment site showing site design (Inset shows Thetaprobes at 10, 20, 40, and 60 cm depth).

Figure 3 .
Figure 3. (a) Average water table depth at each slope location during 2006 and 2007 (least squares mean and 95% CI, n = 69) and (b) average seasonal water table depths at all locations for 2006 and 2007 (least square means and 95% CI, spring and summer n = 84, fall n = 39).

Figure 5 .
Figure 5. Average depth to water table at the summit slope location for each hazard class (error bars represent 95% CI, high n = 24, moderate n = 20, low n = 25).

Figure 6 .
Figure 6.Seasonal average depth to water table at toe slope locations (error bars represent 95% CI, n = 27 spring and summer n = 13 for fall).

Table 1 .
Hazard indicators used to assess likelihood of wet ground areas in the Vanderhoof Forest District through the a priori approach.

Table 2 .
The a priori predicted hydrologic hazard for studied watersheds, their harvest level and the assessment approach used in 2006-2007.

Table 3 .
Particle-size distribution in top 10 cm soil in the MPB and clearcut for the five selected sites (n = 4).

Table 4 .
Precipitation trends between 1980 and 2007 at the Vanderhoof climate station.

Table 5 .
Comparison of saturated hydraulic conductivity (Ksat) between forested and clear-cut areas for the five selected sites.Kmean is geometric mean Kfs value, Kmax is maximum Kfs value, Kmin is minimum Kfs value, CV is coefficient of variation, and Db is soil bulk density.

-5 cm) Db (5-10 cm) (mm h −1 ) (%) (Kg/m −3 ) (Kg/m −3 )
Number of measurements; ‡ Different letters following geometric mean Kfs indicate significant differences between Forested and clear-cut within the same site at P < 0.01; § Different letters following soil bulk density indicate significant differences between Forested and clear-cut within the same site at P < 0.01.

Table 7 .
Watershed hazard prediction for a priori and post-hoc assessment, field data verification summary for volumetric soil moisture and water table elevation along with comments on whether the prediction was correct, over-or underestimated.Table values identify whether the average condition was wet or dry during the summer months of 2007, the wettest summer during the sample period.Post-hoc hazard prediction and hazard scores for 2007 are not included for those watersheds used to generate the post-hoc model.