UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Forest recovery from mountain pine beetle attack : synthesis and simulations of stand carbon and water… Meyer, Gesa 2018

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2018_may_meyer_gesa.pdf [ 6.06MB ]
Metadata
JSON: 24-1.0365753.json
JSON-LD: 24-1.0365753-ld.json
RDF/XML (Pretty): 24-1.0365753-rdf.xml
RDF/JSON: 24-1.0365753-rdf.json
Turtle: 24-1.0365753-turtle.txt
N-Triples: 24-1.0365753-rdf-ntriples.txt
Original Record: 24-1.0365753-source.json
Full Text
24-1.0365753-fulltext.txt
Citation
24-1.0365753.ris

Full Text

Forest recovery from mountain pine beetle attack: synthesis and simulations of stand carbon and water balances using a modified version of the 3-PG model by  Gesa Meyer  M.Sc., Rheinische Friedrich-Wilhelms-Universitaet Bonn, Germany, 2010   A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Soil Science)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  April 2018  © Gesa Meyer, 2018 ii  Abstract The most recent mountain pine beetle (MPB) (Dendroctonus ponderosae) outbreak in British Columbia (BC), which began in the late 1990s, killed ~54% of the mature merchantable lodgepole pine volume and was expected to impact gross primary productivity (GPP), ecosystem respiration (R) and thus net ecosystem productivity (NEP), as well as evapotranspiration (E), snow accumulation and melt in infested stands due to tree mortality. To quantify these effects, eddy-covariance (EC) measurements of carbon (C) and water vapour fluxes have been made above two not-salvage-harvested MPB-attacked pine stands, one with little understory (MPB-06) and another with considerable understory for ten and six years, respectively, and for three years in a partial-salvage-harvested stand, complemented with short-term EC measurements in nearby clearcuts. To determine long-term recovery of the C and water balances following attack, I modified the 3-PG (Physiological Principles Predicting Growth) model to simulate the effects of MPB attack on MPB-06. Modifications included a 2-layer canopy with a partly-dying overstory and growing understory, water availability from snowmelt, and a heterotrophic respiration sub-model. Modelled monthly and annual fluxes at MPB-06 agreed well with the respective EC-estimated values during the decade following attack. Modelled annual GPP, R, NEP and E decreased by about 52%, 35%, 126% and 62%, respectively, in the first year following attack compared to pre-attack values in 2005. While modelled GPP and R, as well as EC-estimated GPP, showed a relatively steady increase over the following decade, EC-estimated R changed little in the first eight years after attack and then increased in the last two years. Both modelled and measured NEP increased significantly over the decade with MPB-06 becoming C neutral within three to four years following attack. EC-measured annual E remained remarkably stable iii  for five years after the attack, and then increased in the last five years, whereas the model indicated a relatively steady increase over the decade. Model projections for five climate change scenarios show 2026 average GPP, R, NEP and E being 14%, 1%, 65% and 5%, respectively, lower than pre-attack values. The quick recovery suggests that not-salvage-harvesting can be a beneficial management practice for C sequestration and hydrology.   iv  Lay Summary Across western North America, tree mortality due to insect attacks such as the recent large-scale mountain pine beetle (MPB) outbreak in British Columbia (BC), has impacted forest growth, and hence timber production and carbon (C) sequestration. The goal of this thesis was to synthesize measurements of the C uptake and release, and thus the C balance, of three lodgepole pine-dominated stands, two not-harvested and one partial-harvested (only dead trees were removed), in the BC Interior over the decade following MPB attack. Also, a forest growth model was modified to include the effects of insect attack on stands. Modelled and measured C and water balances of an attacked stand agreed well. Model projections of C and water fluxes for the following decade were made. The remaining trees and the growing understory enabled the not-salvage-harvested stand to recover from insect attack within a decade, which can inform post-attack forest management decisions.  v  Preface I was responsible for the development of the research questions and the modelling aspect of the research project in consultation with my supervisory committee. I performed the data analysis and quality control of the measurements made since 2011 at two forest research sites and for two years at the third site, modified the process-based model 3-PG, performed the model runs, composed the thesis chapters and manuscripts, and revised them for publication in scientific journals. T. Andrew Black, Rachhpal (Paul) S. Jassal and Nicholas C. Coops were the supervisory authors, and throughout the project were involved with the design of the research and the scientific writing. The other co-authors were involved in the data collection, technical design of the field work, collection and analysis of supplementary data such as leaf area index and snowfall measurements, and provided input to the manuscript writing. A version of Chapter 2 has been published in Forest Ecology and Management as: Meyer G., Black T.A., Jassal R.S., Nesic Z., Grant N.J., Spittlehouse D.L., Fredeen A.L., Christen A., Coops N.C., Foord V.N., Bowler R., 2017. Measurements and simulations using the 3-PG model of the water balance and water use efficiency of a lodgepole pine stand following mountain pine beetle attack. Forest Ecology and Management. 393, 89-104. DOI:10.1016/j.foreco.2017.03.019. A version of Chapter 3 has been published in Forest Ecology and Management as: Meyer G., Black T.A., Jassal R.S., Nesic Z., Coops N.C., Christen A., Fredeen A.L., Spittlehouse D.L., Grant N.J., Foord V.N., Bowler R., 2018. Simulation of net ecosystem productivity of a lodgepole pine forest after mountain pine beetle attack using a modified version of 3-PG. Forest Ecology and Management. 412, 41-52. DOI:10.1016/j.foreco.2018.01.034. vi  Table of Contents  Abstract .......................................................................................................................................... ii Lay Summary ............................................................................................................................... iv Preface .............................................................................................................................................v Table of Contents ......................................................................................................................... vi List of Tables ..................................................................................................................................x List of Figures .............................................................................................................................. xii List of Abbreviations ................................................................................................................. xxi List of Symbols ...........................................................................................................................xxv Acknowledgements ................................................................................................................ xxviii Chapter 1: Introduction ................................................................................................................1 1.1 Background and motivation ............................................................................................ 1 1.1.1 MPB attack history in North America ........................................................................ 4 1.1.2 Effects of MPB attack on C and water balances ......................................................... 7 1.1.3 Forest Productivity Models ....................................................................................... 10 1.1.4 Physiological Principles in Predicting Growth Model (3-PG) ................................. 12 1.2 Research objectives ....................................................................................................... 16 1.3 Thesis overview ............................................................................................................ 18 Chapter 2: Measurements and simulations using the 3-PG model of the water balance and water use efficiency of a lodgepole pine stand following mountain pine beetle attack..........19 2.1 Introduction ................................................................................................................... 19 2.2 Methods......................................................................................................................... 23 vii  2.2.1 The model ................................................................................................................. 23 2.2.2 Modifications to the model ....................................................................................... 31 2.2.3 Measurements ........................................................................................................... 34 2.2.4 Uncertainty analysis of EC fluxes............................................................................. 38 2.2.5 Model parameters...................................................................................................... 39 2.2.6 Future projections ..................................................................................................... 39 2.3 Results ........................................................................................................................... 40 2.3.1 Climate ...................................................................................................................... 40 2.3.2 Model results ............................................................................................................. 43 2.3.2.1 Changing stem density and estimated LAI ....................................................... 43 2.3.2.2 Comparison of modelled and measured E ........................................................ 45 2.3.2.3 Monthly and annual GPP and WUE ................................................................. 49 2.3.2.4 Stand water balance .......................................................................................... 52 2.3.2.5 Sensitivity to model parameter values .............................................................. 54 2.3.2.6 Model projections for the next decade .............................................................. 57 2.4 Discussion ..................................................................................................................... 60 2.5 Conclusions ................................................................................................................... 65 Chapter 3: Simulation of net ecosystem productivity of a lodgepole pine forest after mountain pine beetle attack using a modified version of 3-PG ...............................................67 3.1 Introduction ................................................................................................................... 67 3.2 Methods......................................................................................................................... 70 3.2.1 The model ................................................................................................................. 70 3.2.1.1 Model parameter values .................................................................................... 74 viii  3.2.2 Measurements ........................................................................................................... 76 3.2.3 Uncertainty analysis of EC and modelled fluxes ...................................................... 80 3.3 Results and Discussion ................................................................................................. 81 3.3.1 Time series of the climate variables .......................................................................... 81 3.3.2 Comparison of measured and modelled C balance components .............................. 84 3.3.2.1 Treefall and coarse woody debris ..................................................................... 84 3.3.2.2 Heterotrophic respiration .................................................................................. 86 3.3.2.3 Ecosystem respiration ....................................................................................... 87 3.3.2.4 Gross primary productivity ............................................................................... 92 3.3.2.5 Net ecosystem productivity............................................................................... 98 3.4 Conclusions ................................................................................................................. 103 Chapter 4: Synthesis of measurements of the carbon and water balances of lodgepole pine stands in the BC Interior following mountain pine beetle attack ..........................................105 4.1 Introduction ................................................................................................................. 105 4.1.1 MPB attack history in North America .................................................................... 105 4.1.2 Effects of MPB attack on C and water balances ..................................................... 107 4.1.3 Objectives ............................................................................................................... 110 4.2 Methods....................................................................................................................... 110 4.2.1 Site descriptions ...................................................................................................... 110 4.2.2 Measurements ......................................................................................................... 113 4.2.3 Flux quality control and data analysis .................................................................... 116 4.2.4 Uncertainty analysis of EC fluxes........................................................................... 119 4.3 Results and Discussion ............................................................................................... 120 ix  4.3.1 Climate .................................................................................................................... 120 4.3.2 Changing stand structure......................................................................................... 123 4.3.3 Water fluxes following MPB attack ....................................................................... 127 4.3.4 Carbon balances ...................................................................................................... 131 4.3.5 Modelling predictions of R and NEP at MPB-06 with 3-PG .................................. 142 4.4 Conclusions ................................................................................................................. 145 Chapter 5: Conclusions .............................................................................................................147 5.1 Limitations of the research .......................................................................................... 151 5.2 Future outlook ............................................................................................................. 152 Bibliography ...............................................................................................................................155 Appendices ..................................................................................................................................165 Appendix A ............................................................................................................................. 165 Appendix B ............................................................................................................................. 171 Appendix C ............................................................................................................................. 175  x  List of Tables Table 1.1 List of species to which 3-PG has been applied and the respective studies. ................ 14 Table 2.1  List of the 3-PG parameters and their values for lodgepole pine (Pinus contorta var. latifolia) used in this study. .................................................................................................. 24 Table 2.2  Details of climate scenarios for air temperature (Ta), precipitation (P) and vapour pressure deficit (VPD) during 2016 – 2025 used for predicting the effect of MPB attack on future carbon and water fluxes at MPB-06 over the next decade. ....................................... 40 Table 2.3  Measured annual precipitation (P), measured and modelled annual evapotranspiration (E), drainage below the 60-cm depth (Dr) and water use efficiency (WUE) for 2005 - 2015 at MPB-06 and the percentage of half-hourly values of measured E required to be gap-filled per year. ................................................................................................................................ 48 Table 2.4  Annual modelled evapotranspiration (E), gross primary productivity (GPP) and water use efficiency (WUE) at MPB-06 for 2005, 2006, 2015 and 2025. Projections were done for the five climate change scenarios (S1–S5) shown in Table 2.2. The results for S2 represent the predicted trend using the linear regression equation. ..................................................... 59 Table 3.1  List of parameters used in the heterotrophic respiration (Rh) sub-model in the revised 3-PG, and their values for lodgepole pine (Pinus contorta var. latifolia). ........................... 75 Table 3.2 List of value ranges for the model parameters used in the uncertainty calculation of modelled fluxes. ................................................................................................................... 81 Table 3.3 Slopes of the regression, intercepts, coefficients of determination (R2), root mean square errors (RMSE), p-values and sample sizes (n), for EC-estimated monthly and annual values versus modelled monthly and annual values of gross primary productivity (GPP), xi  ecosystem respiration (R) and net ecosystem productivity (NEP) for 2007 - 2016 at MPB-06 shown in Figs. 6, 9 and 12. .............................................................................................. 90 Table 3.4  EC-estimated and modelled annual gross primary productivity (GPP), ecosystem respiration (R) with the two modelled components, autotrophic (Ra) and heterotrophic (Rh) respiration, and net ecosystem productivity (NEP) for 2005 - 2016 at MPB-06. .............. 102 Table 4.1  Stand characteristics at the three main sites MPB-06, MPB-03 and MPB-09. ......... 111 Table 4.2  Annual averages of photosynthetically active radiation (PAR), cumulative precipitation (P), vapour pressure deficit (VPD), air temperature (Ta) and volumetric soil water content (θ) for the 0-10 cm layer for 2007 – 2016 at MPB-06, MPB-03 and MPB-09. Values of VPD were calculated from daytime values. ....................................................... 122 Table 4.3  Measured canopy openness and overstory leaf area index (LAI) at MPB-06 and MPB-03 for 2007-2015 obtained from hemispherical photography. ........................................... 125 Table 4.4  Annual totals of net ecosystem productivity (NEP), gross primary productivity (GPP), ecosystem respiration (R), evapotranspiration (E) and precipitation (P) - E at MPB-06, MPB-03 and MPB-09 for 2007-2016. ................................................................................ 134 Table 4.5  Growing season (GS, defined as 1 May – 30 September) totals of net ecosystem productivity (NEP), gross primary productivity (GPP), ecosystem respiration (R) and evapotranspiration (E) at MPB-06, MPB-03, MPB-09 and MPB-09C for 2007-2016 and average daily NEP and E at CC-97, CC-05 and MPB-06 for 29 June - 23 July 2007 and 24 July - 16 August 2007, respectively. Values shown for MPB-06 for CC-97 and CC-05 are for the same period the clearcut measurements were made. .............................................. 135 Table C.1 Root zone soil water storage (S) at MPB-06, MPB-03 and MPB-09 at the beginning of the growing season (May 1) of the decade following attack. ............................................. 175 xii  List of Figures Figure 1.1  Mountain pine beetle infestation map including observed percentage of pine volume killed in BC by 2015 based on the Provincial Aerial Overview of forest health from 1999 to 2015 and projections using the stochastic model BCMPB version 13 assuming no forest harvesting or beetle suppression activities. Source: https://www.for.gov.bc.ca/ftp/hre/external/!publish/web/bcmpb/year13/BCMPB.v13.2015Kill.pdf. ..................................................................................................................................... 7 Figure 1.2  Locations of the MPB sites in lodgepole pine-dominated stands attacked in 2006 (MPB-06), 2003 (MPB-03) and 2005 and 2006 (MPB-09), respectively, where eddy-covariance flux measurements have been made. MPB-06 and MPB-03 were not harvested, MPB-09 was partial-salvage-harvested in 2009 while MPB-09C, nearby MPB-09, was clearcut-salvage-harvested. CC-05 and CC-97 were clearcuts nearby MPB-06 harvested in 2005 and 1997, respectively. ................................................................................................ 17 Figure 2.1  Basic structure of 3-PG including modifications made (in red) for its application to MPB-attacked lodgepole pine stands. The variables involved in the water balance are shown in blue and the allocation of biomass is shown in green. The subscripts ‘o’ and ‘u’ stand for the overstory and understory, respectively. α is the quantum efficiency of photosynthesis, φ is the physiological modifier applied to the maximum canopy conductance (gcmax) and α, ws is the average tree stem biomass and γn is the tree mortality rate. The visual basic code for the modified version of 3-PG is available on request (gesa.meyer@ubc.ca). .......................................................................................................... 28 xiii  Figure 2.2  Calculated monthly live overstory stem density at MPB-06 using an initial stand density value (900 stems ha-1) and γn from Equation (7) for 2005 – 2015 (curve) (γn0 = 31.6 % month-1 and t = 5.6 months) and observed stand density (circles). ................................ 36 Figure 2.3  Monthly 24-h average photosynthetically active radiation (PAR), cumulative precipitation (P), vapour pressure deficit (VPD) (solid grey line) and specific humidity deficit (D) (dashed black line), air temperature (Ta) (dashed line), soil temperature at the 10-cm depth (Ts) (solid line) and soil water content (θ) for the 0-10-cm layer for 2007 - 2015 at MPB-06. Values of VPD and D were calculated from daytime average values. ................. 42 Figure 2.4  Monthly modelled and measured LAI at MPB-06 for 2005 - 2015. The modelled contributions of the two parts of the stand, the live understory and the live overstory LAI, are shown in green and red, respectively, and the total live LAI is the blue line. The dashed black line is the modelled dead LAI and the dashed light blue line shows the total overstory LAI. The measured total overstory LAI values obtained from hemispherical photography are the black circles. ............................................................................................................. 44 Figure 2.5  Modelled and measured monthly E at MPB-06 for 2005 - 2015. The contributions of the two parts of the stand, the understory and the overstory, are shown in green and red, respectively, and the total is the blue line. The measured values are the black circles. ....... 46 Figure 2.6  Measured monthly values versus modelled monthly values of E at MPB-06 for 2007-2015 including the 1:1 line (solid black line) and regression line (dashed red line). Also shown is the linear regression equation, coefficient of determination (R2), root mean square error (RMSE) and sample size (n). ....................................................................................... 47 xiv  Figure 2.7  Modelled and measured annual E and GPP at MPB-06 for 2005 - 2015. Bars show the 95% confidence interval for annual E and GPP including the uncertainty estimates described in Section 2.2.4. .................................................................................................... 48 Figure 2.8  Modelled and EC-estimated monthly GPP at MPB-06 for 2005 - 2015. The contributions of the two parts of the stand, the understory and the overstory, are shown in green and red, respectively, and the total is the blue line. The EC-estimated values are the black circles. ......................................................................................................................... 51 Figure 2.9  EC-estimated monthly values versus modelled monthly values of GPP at MPB-06 for 2007-2015 including the 1:1 line (solid black line) and regression line (dashed red line). Also shown is the linear regression equation, coefficient of determination (R2), root mean square error (RMSE) and sample size (n). ........................................................................... 52 Figure 2.10  Mean monthly water balance components and annual totals for MPB-06 for 2007-2015: measured precipitation (P), modelled evapotranspiration (E), modelled potential evapotranspiration (Ep), modelled drainage (Dr), modelled change in soil water storage (S), and liquid water (Lw), measured as rain plus modelled snowmelt and equals P on a water year basis. ................................................................................................................... 53 Figure 2.11  Sensitivity of annual E for 2007 to the following model parameters: PAR extinction coefficient (k = 0.6), soil fertility (FR = 0.4), stomatal sensitivity to VPD (sD = 0.055 mbar-1), quantum efficiency (α = 0.047 mol of C (mol)-1 of PAR photons), maximum stand age (ymax = 150 years), boundary layer conductance (ga = 0.2 m s-1), maximum canopy conductance (gcmax = 0.012 m s-1), and maximum interception (MaxIntcptn = 0.15). Absolute values of E are on the left-hand y-axis and relative values of E are on the right-hand y-axis. ........................................................................................................................... 56 xv  Figure 2.12  Annual modelled E (a), GPP (b) and WUE (c) at MPB-06 for 2016-2025 for the five climate change scenarios (S1–S5) shown in Table 2.2. For S2 the black dots are the annual values and the dashed black line is the regression line. Also shown are the linear regression equations for S2. ................................................................................................................... 59 Figure 3.1  Basic structure of 3-PG including modifications made (in red) for its application to MPB-attacked lodgepole pine stands. The subscripts ‘o’ and ‘u’ stand for the overstory and understory, respectively. α is the canopy quantum use efficiency of photosynthesis, φ is the physiological modifier applied to α and the maximum canopy conductance (gcmax), ws is the average tree stem biomass, γn is the tree mortality rate, Soil C and CWD C are the soil and coarse woody debris (CWD) C pools, and Rh soil and Rh cwd are the heterotrophic respiration components originating from the soil and CWD, respectively. In this version of 3-PG, soil C is held constant (see Section 3.2.1.1) with no input from fine root turnover or litter. The visual basic code for the modified version of 3-PG is available on request (gesa.meyer@ubc.ca). .......................................................................................................... 72 Figure 3.2  (a) Calculated monthly treefall rate (black curve) at MPB-06 using a Gaussian fit through the monthly average of the observed annual treefall rate (blue bars) and (b) cumulative monthly modelled coarse woody debris (CWD) C pool content for 2005 - 2016. .............................................................................................................................................. 76 Figure 3.3  Monthly 24-h average photosynthetically active radiation (PAR), cumulative precipitation (P), vapour pressure deficit (VPD), air temperature (Ta) (grey line), soil temperature at the 5-cm depth (Ts) (black line) and soil water content (θ) in the 0-10-cm layer (black line) and the 30-50-cm layer (grey line) for 2007 - 2016 at MPB-06. Values of VPD were daytime average values. ...................................................................................... 83 xvi  Figure 3.4  Modelled monthly Rh at MPB-06 for 2005-2016. The contributions of the two components of Rh, coming from the soil (Rh soil) and coarse woody debris (Rh cwd), are shown in red and green, respectively, and total Rh is the blue line. ................................................. 87 Figure 3.5  Modelled (Rmod) and EC-estimated (REC) monthly R at MPB-06 for 2005 - 2016. The modelled contributions of the two components of R, heterotrophic (Rh) and autotrophic (Ra) respiration, are shown in green and red, respectively, and the total is the blue line. The EC-estimated values are the black circles including the 95% confidence interval represented by the vertical bars. ................................................................................................................... 89 Figure 3.6  EC-estimated monthly (a) and annual (b) values versus modelled monthly (a) and annual (b) values of R at MPB-06 for 2007-2016 with the 1:1 line (solid black line) and regression line (dashed red line). See Table 3.3 for statistics. ............................................. 90 Figure 3.7  Modelled (Rmod) and EC-estimated (REC) annual R at MPB-06 for 2005 - 2016. The red triangles show the modelled values, and the EC-estimated values are the blue filled circles. The red and blue shaded areas show, respectively, the standard deviation for annual modelled R, and annual EC-estimated R with the uncertainty estimates described in Section 3.2.3. ..................................................................................................................................... 91 Figure 3.8  Modelled (GPPmod) and EC-estimated (GPPEC) monthly GPP values at MPB-06 for 2005 - 2016. The blue line shows the modelled values, and the EC-estimated values are the black circles with the 95% confidence interval represented by the vertical bars. ................ 94 Figure 3.9  EC-estimated monthly (a) and annual (b) values versus modelled monthly (a) and annual (b) values of GPP at MPB-06 for 2007-2016 with the 1:1 line (solid black line) and regression line (dashed red line). See Table 3.3 for statistics. ............................................. 95 xvii  Figure 3.10  Modelled (GPPmod) and EC-estimated (GPPEC) annual GPP values at MPB-06 for 2005 - 2016. The red triangles show the modelled values, and the EC-estimated values are the blue filled circles. The red and blue shaded areas show, respectively, the standard deviation for annual modelled and EC-estimated GPP with the uncertainty estimates described in Section 3.2.3. .................................................................................................... 96 Figure 3.11  Modelled (NEPmod) and measured (NEPEC) monthly NEP values at MPB-06 for 2005 - 2016. The blue line shows the modelled values and the measured values are the black circles with the 95% confidence interval represented by the vertical bars. ................ 99 Figure 3.12  Measured monthly (a) and annual (b) values versus modelled monthly (a) and annual (b) values of NEP at MPB-06 for 2007-2016 with the 1:1 line (solid black line) and regression line (dashed red line). See Table 3.3 for statistics. ........................................... 100 Figure 3.13  Modelled (NEPmod) and EC-estimated (NEPEC) annual NEP values at MPB-06 for 2005 - 2016. The red triangles show the modelled values, and the EC-estimated values are the blue filled circles. The red and blue shaded areas show, respectively, the standard deviation for annual modelled and measured NEP with the uncertainty estimates described in Section 3.2.3. .................................................................................................................. 101 Figure 4.1  Monthly 24-h average photosynthetically active radiation (PAR), cumulative precipitation (P), monthly 24-h average vapour pressure deficit (VPD), monthly 24-h average air temperature (Ta) and monthly 24-h average volumetric soil water content (θ) for the 0-10 cm layer for 2007 – 2016 at MPB-06 (black dashed-dotted line), MPB-03 (blue line) and MPB-09 (red dotted line). Values of VPD were daytime average values. .......... 121 Figure 4.2  Mountain pine beetle attack status at MPB-06 for August of 2006 – 2016. Values for 2006-2010 are from Brown et al. (2014) and values for 2010-2016 are from Dale Seip xviii  (personal communication). The bars represent the percentage of healthy (in black), green-attacked (in white) and red-/grey-attacked (in grey) canopy trees. .................................... 124 Figure 4.3  Percentage of monthly incoming photosynthetically active radiation (PAR) reaching below the canopy (3-m height) at MPB-06, MPB-03 and MPB-09 for the decade following attack. ................................................................................................................................. 127 Figure 4.4  Annual evapotranspiration (E) at MPB-06, MPB-03 and MPB-09 for the decade following attack. Bars show the 95% confidence interval for annual E. ........................... 128 Figure 4.5  Cumulative values of daily precipitation (P) minus evapotranspiration (E) and daily change in root zone soil water storage (ΔS) at MPB-06, MPB-03 and MPB-09 for the growing seasons (May 1 – September 30) of the decade following attack. The soil water storage (S) corresponding to ΔS = 0 was 36.34 mm, 7.64 mm and 17.05 mm at MPB-06, MPB-03 and MPB-09, respectively, averaged over the observed years. For values of S for each year see Appendix C. Note that the y-axis ranges for P – E and ΔS are different for each stand. .......................................................................................................................... 130 Figure 4.6  Annual net ecosystem productivity (NEP), gross primary productivity (GPP), and ecosystem respiration (R) at MPB-06, MPB-03 and MPB-09 for the decade following attack. Bars show the 95% confidence interval for annual NEP, GPP and R. The years when MPB attack occurred were 2003 (estimated), 2006 and 2006 (estimated) for MPB-03, MPB-06 and MPB-09, respectively. .................................................................................. 132 Figure 4.7  Monthly, April to October, means of EC-derived gross primary productivity (GPP), evapotranspiration (E) and water use efficiency (WUE) for 2007 - 2016 at MPB-06, for 2007 – 2012 at MPB-03 and for 2010 – 2012 at MPB-09. ................................................ 139 xix  Figure 4.8  Monthly gross primary productivity (GPP) versus evapotranspiration (E) for 2007 - 2016 at MPB-06 (a), for 2007 – 2012 at MPB-03 (b) and for 2010 – 2012 at MPB-09 (c) including the respective regression lines (dashed red line). Also shown are the linear regression equations, coefficients of determination (R2), root mean square errors (RMSE) and sample sizes (n) obtained from EC measurements. ..................................................... 140 Figure 4.9  Annual modelled gross primary productivity (GPP) (a), ecosystem respiration (R) (b), and net ecosystem productivity (NEP) (c) at MPB-06 for 2017-2026 for the five climate change scenarios (S1-S5). For S2 the black circles are the annual values and the dashed black line is the regression line. Also shown are the linear regression equations for S2. .. 145 Figure A.1 Photograph of the eddy-covariance (EC) system at the 26-m height at MPB-06 in August 2007, showing a three-dimensional ultrasonic anemometer (model CSAT3, Campbell Scientific Inc. (CSI), Logan, Utah), an open-path infrared gas analyser (IRGA) (model LI-7500, LI-COR, Inc., Lincoln, Nebraska) and a fine-wire thermocouple. ......... 165 Figure A.2 Photograph of the 32-m tall measurement tower at MPB-09 in July 2011 with the eddy-covariance (EC) and net radiometer booms, three solar panels and the rain gauge (model 52203 R.M. Young Co., Traverse City, MI at MPB-09). ...................................... 166 Figure A.3 Canopy photographs taken from the top of the flux tower at MPB-06 in a northern direction in 2006, 2007, 2010 and 2014. ............................................................................ 167 Figure A.4 Canopy photographs taken from the top of the flux tower at MPB-03 in 2008 and 2014. ................................................................................................................................... 168 Figure A.5 Canopy photograph and view from the ground at MPB-09 in 2009. ....................... 168 Figure A.6 Photographs of the measurement towers at a) CC-97 and b) CC-05 in June and August 2007, respectively. ................................................................................................. 169 xx  Figure A.7 Measurement tower at MPB-09C in July 2010. ....................................................... 170 Figure B.1 Manual portable chamber system with a closed-path infrared gas analyser (LI-800, LI-COR), using an opaque chamber measuring respiration on a log, and a temperature probe. .................................................................................................................................. 171 Figure B.2 Collars on two dead trees on which bole respiration measurements were made. ..... 172 Figure B.3 Collars on two live trees on which bole respiration measurements were made. ...... 172 Figure B.4 Collars on two logs on which bole respiration measurements were made. .............. 173 Figure B.5 Collars in locations where soil respiration measurements were made. .................... 174  xxi  List of Abbreviations  Acronym Units Description 3-PG  Use of Physiological Principles in Predicting Growth aH  Constant in the stem height relationship aS  Constant in the stem mass vs. diam. relationship attackStAge years Stand age at time of beetle attack aV  Constant in the stem volume relationship BC  British Columbia C  Carbon CanCoverMax  Maximum canopy cover CBM-CFS3  Carbon Budget Model of the Canadian Forest Sector CLASS  Canadian Land Surface Scheme CLM  Community Land Model CO2  Carbon dioxide CO2e  Carbon dioxide equivalent CSI  Campbell Scientific Inc. CTEM  Canadian Terrestrial Ecosystem Model CWD  Coarse woody debris DBH cm Diameter at breast height EC  Eddy-covariance fAge  Age modifier FC MPa Field capacity fCalpha700  Assimilation enhancement factor at 700 ppm fCg700  Canopy conductance enhancement factor at 700 ppm fCO2  Atmospheric CO2 modifier FCRN  Fluxnet Canada Research Network fFrost  Frost modifier fN0  Value of 'fNutr' when FR = 0 fNn  Power of (1-FR) in 'fNutr'  fPAR  Fraction of absorbed photosynthetically active radiation FR  Soil fertility xxii  Acronym Units Description fracBB  Branch and bark fraction  fracBB0  Branch and bark fraction at age 0 fracBB1  Branch and bark fraction for mature stands fSW  Soil water modifier fT  Temperature modifier fullCanAge years Age at canopy cover  gDM_mol g mol-1 of DM Molecular weight of dry matter GHG  Greenhouse gas GMAO  Global Modeling and Assimilation Office GPP g m-2 t-1 of Ca Gross primary productivity (sign convention is GPP ≥ 0) GS  Growing season H2O  Water vapour IPCC SRES  Intergovernmental Panel on Climate Change Special report on emissions scenarios IRGA  Infrared gas analyser kF days Days production lost per frost day LAI m2 m-2 Leaf area index LAIgcxu m2 m-2 LAI for maximum conductance for understory LAImaxIntcptn m2 m-2 LAI for maximum rainfall interception m0  Value of 'm' when FR = 0 MAI m3 ha-1 yr-1 Mean annual volume increment MaxAge/ ymax years Maximum stand age used in age modifier MaxCond m s-1 Maximum canopy conductance MaxCondu m s-1 Maximum canopy conductance for understory MaxIntcptn  Maximum interception (proportion of rainfall evaporated from canopy) MERRA  Modern Era Retrospective-Analysis for Research and Applications mF  Fraction mean single-tree foliage biomass lost per dead tree mFu  Fraction mean single-tree foliage biomass lost per dead tree in the understory MinCond m s-1 Minimum canopy conductance MODIS  Moderate Resolution Imaging Spectroradiometer molPAR_MJ mol MJ-1 Conversion of solar radiation to PAR xxiii  Acronym Units Description MPB  Mountain pine beetle (Dendroctonus ponderosae) mR  Fraction mean single-tree root biomass lost per dead tree mRu  Fraction mean single-tree root biomass lost per dead tree in the understory mS  Fraction mean single-tree stem biomass lost per dead tree mSu  Fraction mean single-tree stem biomass lost per dead tree in the understory nAge  Power of relative age in function for fAge NDVI  Normalized Difference Vegetation Index  NE  Northeast NEE g m-2 t-1 of C Net ecosystem exchange of CO2 (positive values indicate C loss to the atmosphere) NEP g m-2 t-1 of C Net ecosystem productivity (positive values indicate C gain by the ecosystem) ngammaN  Shape of mortality response nHB  Power of DBH in the stem height relationship nHN  Power of stocking in the stem height relationship NPP g m-2 t-1 of C Net primary productivity (sign convention is NPP ≥ 0) nS  Power in the stem mass vs. diam. relationship nVB  Power of DBH in the stem volume relationship nVN  Power of stocking in the stem volume relationship NW  Northwest PAR µmol m-2 s-1 Photosynthetically active radiation pFS2  Foliage: stem partitioning ratio @ VPD = 2 cm pFS20  Foliage: stem partitioning ratio @ VPD = 20 cm pRn  Minimum fraction of NPP to roots pRx  Maximum fraction of NPP to roots Qa W m-2 Intercept of net v. solar radiation relationship Qb  Slope of net v. solar radiation relationship rAge  Relative age to give fAge = 0.5 RH % Relative humidity RMSE  Root mean square error S1  First Scenario S2  Second Scenario  xxiv  Acronym Units Description S3  Third Scenario S4  Fourth Scenario  S5  Fifth Scenario  SDM  Synchronous device for measurement SE  Southeast SLA m2 kg-1 Specific leaf area SLA0 m2 kg-1 Specific leaf area at age 0 SLA1 m2 kg-1 Specific leaf area for mature leaves SW  Southwest SWconst  Moisture ratio deficit for fθ = 0.5  SWpower  Power of moisture ratio deficit tBB years Age at which fracBB = (fracBB0+fracBB1)/2 Tg  Teragrams tgammaF months Age at which litterfall rate has median value tgammaN years Age at which mortality rate has median value thinPower  Power in self-thinning rule tRho years Age at which rho = (rhoMin+rhoMax)/2 tSLA years Age at which specific leaf area = (SLA0+SLA1)/2 UNFCC  United Nations Framework Convention on Climate Change VPD kPa Vapour pressure deficit WD mm Water deficit WP MPa Permanent wilting point of the soil WUE g of C (kg of H2O)-1 Water use efficiency a t indicates time (seconds, hours, days or years). xxv  List of Symbols  Symbol Units Description alphau molC molPAR-1 Understory quantum efficiency Amax  g m-2 t-1 of Ca Ecosystem photosynthetic capacity sD / CoeffCond mbar-1 Stomatal sensitivity to vapour pressure deficit Ccwd / Ccwd kg m-2 of C Coarse woody debris carbon storage Cdead avg kg tree-1 of C Average C content per standing dead tree Csoil / Csoil kg m-2 of C Soil (0-60-cm depth) carbon storage D g of H2O (g of air)-1 Specific humidity deficit of canopy air Dr mm t-1 Drainage from the root zone (60-cm depth) S mm t-1 Change in root-zone soil water storage E mm t-1 Evapotranspiration Ep mm t-1 Potential evapotranspiration Fc µmol m-2 s-1 Flux of CO2 fCg  Canopy conductance modifier for atmospheric CO2 concentration G W m-2 Soil heat flux ga / BLcond m s-1 Boundary layer conductance gammaF0 month-1 Litterfall rate at t = 0 gammaFx month-1 Maximum litterfall rate gammaN0 % year-1 Seedling mortality rate (t = 0) gammaNx % year-1 Mortality rate for large t gammaR month-1 Average monthly root turnover rate gc m s-1 Canopy conductance gc limited m s-1 Air temperature-limitation on canopy conductance gcmax m s-1 Maximum canopy conductance gcmin m s-1 Minimum canopy conductance H W m-2 Sensible heat flux k  Photosynthetically active radiation extinction coefficient ku  Extinction coefficient for absorption of PAR by understory Ld m2 m-2 Accumulated dead LAI on dying trees before significant xxvi  Symbol Units Description numbers of these needles start falling Ld t m2 m-2 Current dead leaf area index Lgcx / LAIgcx m2 m-2 Leaf area index for maximum canopy conductance Lw mm Liquid water n  Sample size o  Degree P mm t-1 Precipitation p  Probability value  m3 m-3 Volumetric soil water content Q10 / Q10  Relative increase in base-respiration rate for an increase in temperature of 10 °C a m3 m-3 Available volumetric soil water content a init m3 m-3 Available volumetric soil water content at the beginning of the model run a max m3 m-3 Maximum plant available water a min m3 m-3 Minimum plant available water WP  m3 m-3 Volumetric soil water content at the permanent wilting point of the soil R g m-2 t-1 of C Ecosystem respiration (sign convention is R ≥ 0) r1, r2, r3  Empirical constants from the empirical logistic relationship R10 / R10 g m-2 month-1 of C Base-respiration rate at 10 C R10 cwd / R10cwd g kg-1 month-1 of C Base-respiration rate at 10 C of the coarse woody debris R10 soil / R10soil g kg-1 month-1 of C Base-respiration rate at 10 C of the soil R2  Coefficient of determination Ra g m-2 t-1 of C Autotrophic respiration Rh g m-2 t-1 of C Heterotrophic respiration rhoMax kg m-3 Maximum basic density - for older trees rhoMin kg m-3 Minimum basic density - for young trees Rn W m-2 Net radiation Rs MJ m-2 month-1 Incident solar radiation Rs int MJ m-2 month-1 Intercepted solar radiation rw (t)  Time-varying moving window parameter s  Mixing ratio xxvii  Symbol Units Description S mm t-1 Soil water storage Sc  µmol m-2 s-1 CO2 storage in the air column beneath the eddy-covariance sensors t months Time since attack in months Ta °C Air temperature Tmax °C Maximum temperature for growth Tmin °C Minimum temperature for growth Topt °C Optimum temperature for growth Ts °C Soil temperature u* m s-1  Friction velocity u*th m s-1  Friction velocity threshold v m s-1  Lateral wind speed w m s-1  Vertical wind speed ws kg tree-1 Stem biomass wSx1000 kg tree-1 Max. stem mass per tree @ 1000 trees/hectare wSx1000u kg tree-1 Max. stem mass per tree for understory @ 1000 trees/hectare Y  Ratio NPP/GPP ymax years Maximum stand age α / alpha mol of C (mol of PAR photons)-1  Canopy quantum efficiency γn % month-1 Mortality rate γn0 % month-1 Initial rate of mortality ε  Dimensionless rate of change of saturated vapour pressure with respect to air temperature λ J kg-1 Latent heat of vaporization λE W m-2 Latent heat flux ρ kg m-3 Air density τn months Time constant for the dead needles τt months Time constant for MPB induced tree mortality φ  Physiological modifier a t indicates time (seconds, hours, days or years). xxviii  Acknowledgements I would like to thank my supervisor, Dr. T. Andrew Black, for his commitment and support over the years, as well as my supervisory committee members Dr. Rachhpal (Paul) Jassal and Dr. Nicholas Coops for sharing their knowledge and providing input to my research and the manuscripts. I have learned a lot from you and your input helped improve this thesis.  I am also thankful to the rest of the Biomet group, especially Zoran Nesic, Nicholas Grant, and former member Rick Ketler, who were instrumental in the project and helped me deal with the data analysis. I greatly appreciate the assistance from David Spittlehouse, Vanessa Foord, and Rebecca Bowler at the BC Ministry of Forests, Lands, Natural Resource Operations and Rural Development (BCMFLNRRD), who provided me with their data and made their way out to the sites countless times to collect data for us or to help troubleshoot issues at the sites. Thanks also to Dale Seip, and J.A. Trofymow for providing their data. This research was financially supported by the Pacific Institute for Climate Solutions (PICS), the Terrestrial Research on Ecosystems and World-wide Education and Broadcast (TerreWEB) project of the Natural Sciences and Engineering Research Council of Canada (NSERC) CREATE program, the Canadian Foundation for Climate and Atmospheric Sciences (CFCAS), the BC Forest Science Program, BCMFLNRRD, and a NSERC Strategic Project Grant and Discovery Grant to T.A. Black.  Great thanks to my lab mates Jilmarie Stephens and Hughie Jones, who have provided moral support in tough times, good discussions and fun times together. Thanks for being my friends and convincing me to believe in myself and to not give up! I am also thankful to lots of great and inspiring people, who I have met during the years in Vancouver. To name a few, xxix  thanks to Michael Mitchell, Anna Mittelholz, Anna Grau, Geneviève Savard, Thibaut Astic and Marie-Pier Malouin, who have been great friends, have supported me and taught me lots of things, especially about sports and the outdoors. Thanks for introducing me to several different activities, countless hours in the climbing gym, great dinners and providing a counterpart to work and possibilities to explore parts of BC! Special thanks to my parents, Annette Meyer-Müsch and Bodo Meyer, my sister, Inga, and Raamu! You made this PhD possible with your constant support, morally and financially, convincing me to never give up, as well as introducing me to different countries. Thank you for all the opportunities you have provided me with. Dad, I am sure you are still with me in your thoughts! 1  Chapter 1: Introduction  1.1 Background and motivation Globally, forest ecosystems provide a number of ecosystem services including the regulation of watershed hydrology, air quality, climate and provisioning of timber (Molnar and Kubiszewski, 2012). They also play an important role in carbon (C) sequestration by fixing atmospheric carbon dioxide (CO2) through photosynthesis, thereby mitigating climate change. However, due to natural and anthropogenic disturbances, CO2 loss through ecosystem respiration (R) can exceed CO2 fixation through photosynthesis. Anthropogenic CO2 emissions play an important role in climate change. Largely due to industrial activities including deforestation, global mean annual CO2 concentrations in the atmosphere have increased continuously from 280 ppm before the beginning of the Industrial Revolution (Canadell et al., 2007) to about 403 ppm in 2016. Monthly mean values exceeded 400 ppm for the first time in April 2014 at the Mauna Loa Observatory in Hawaii, and in March 2015 globally (NOAA, 2018). About a third of anthropogenic CO2 emissions are absorbed by the terrestrial biosphere, similar to the amount absorbed by the oceans (Mystakidis et al., 2016; Bellassen and Luyssaert, 2014). The response of vegetation to increased CO2 levels still remains uncertain, raising questions about the future of this C sink with upcoming changes in climate (Mystakidis et al., 2016). Using measurements to constrain water and C fluxes, uncertainties in modelled terrestrial C uptake can be reduced (Mystakidis et al., 2016). An enhancement of terrestrial gross primary productivity (GPP) due to increases in atmospheric CO2 concentration is generally predicted by climate and C cycle models, however the magnitude of this increase ranges from 20 - 60% increases in GPP 2  predicted with a doubling of atmospheric CO2 concentration (Wenzel et al., 2016). Using measurements of atmospheric CO2 at Point Barrow, Alaska, and Cape Kumukahi, Hawaii, Wenzel et al. (2016) estimated a doubling of CO2 to result in GPP being 37 ± 9% higher in high-latitude ecosystems and 32 ± 9% higher in extratropical ecosystems, respectively.  Canada’s managed forests, defined as forests under direct human influence managed for wood production, ecological conservation or protection from natural disturbances, which account for 65% of all forests in Canada with an area of ~226 million ha (Environment Canada, 2017), are generally estimated to be C sinks, e.g., of ~99 Mt of CO2 equivalent (CO2e) in 1990 and of ~72 Mt of CO2e in 2000. Since 2002, the C balance of Canada’s forests has been severely impacted by disturbances with the estimated sink of ~26 Mt of CO2e in 2015 in managed forest areas, that were harvested and regenerated and taking into account the disposal of harvested wood products, being offset by emissions due to natural disturbances of ~247 Mt CO2e, resulting in a net CO2e emission of ~221 Mt CO2e by Canadian managed forests in 2015 (Natural Resources Canada, 2017). Thus, large scale disturbances such as fires, which burned the largest area in the managed forests since 1995 with ~2 million ha burned in 2015, and the ongoing effect of insect outbreaks, significantly impact the ability of Canadian forests to remain a C sink. This makes the prediction of the greenhouse gas (GHG) balance of Canada’s forests under climate change even more challenging, as the increased C source strength due to enhanced forest growth could be offset by increased mortality or fire and insect outbreaks resulting in some regions being C sources (Natural Resources Canada, 2017). GHG balances for unmanaged forests are not reported under the United Nations Framework Convention on Climate Change (UNFCC) adding uncertainty to the GHG balance of Canada’s forests as a whole. 3  The mountain pine beetle (MPB) (Dendroctonus ponderosae) is native to western North America with a range extending from northern Mexico to northwestern BC (Safranyik and Carroll, 2006). The MPB outbreak in BC, which started in the late 1990s, affecting mainly lodgepole pine (Pinus contorta var. latifolia) stands, triggered the loss of an estimated 731 million m3 of commercial timber. Mitigation and restoration measures, ensuring the sustainability of forests, as well as the degree to which forest stands can be harvested in order to produce timber, paper and energy products, are regulated by BC’s policy framework in the Forest Act. Nicholls (2017) describes the determination of the annual allowable cut in the Prince George Timber Supply area in response to insect outbreaks, for example, and the potential changes in forest management to adapt to climate change. There has been an increased interest in renewable sources of energy globally, driven by concerns about energy security and climate change in recent decades. The production of energy from wood, in addition to an increased use of wood as building materials, could help mitigate climate change. Currently, proponents and opponents of increased forest-derived energy production disagree as to whether logging beetle-affected forests will negatively or positively affect the release of GHG emissions and thus BC's ability to meet its GHG emission reduction targets. Turning portions of the affected timber into bioenergy seems a logical solution as the dead trees, which have little remaining value as raw material for lumber and wood pulp products, will otherwise just rot and decay and lose their stored C. Opponents, on the other hand, worry about increased deforestation, declining biological diversity, increased GHG emissions and thus climate change. Management following beetle attack in BC has employed one of the three options: not harvesting, partial-salvage harvesting and clearcut-salvage harvesting (Coates and Hall, 2005). The aim of this study was to collect and evaluate data that can serve as the basis for better-4  informed policy decisions related to forest management by observing C and water balances during the decade following attack and predicting how forest C balances in the beetle-affected area will change over the next decade. It is important for scientists as well as forest managers and policy makers to find out how the affected stands are recovering from the beetle attack and how the recovery times differ depending on the management strategy. To make informed decisions on how to deal with this widespread disturbance, it is important to know the consequences of different management strategies. As Odum (1969) stressed in his classic paper, this requires balancing between young stands or early successional ecosystems, providing high productivity and growth rates, and mature ecosystems, which provide a protective environment and deliver ecosystem services such as wildlife habitat, water purification and nutrient cycling.   1.1.1 MPB attack history in North America The main host species of the MPB in western Canada are lodgepole pine (Pinus contorta Dougl. ex Loud. var. latifolia Engelm.), ponderosa pine (Pinus ponderosa P. Laws. ex C. Laws) and western white pine (P. monticola Dougl. Ex D. Don), although other pine species such as sugar pine (P. lambertiana Dougl.), eastern white pine (P. strobus L.), jack pine (P. banksiana Lamb.), and Scots pine (P. sylvestris L.) can also be infested and killed (Safranyik and Carroll, 2006). Cullingham et al. (2011) have confirmed that jack pine trees have been successfully attacked by the MPB and that beetles were able to complete their larval development in jack pine, which is an important step in possible host-range expansion. Young beetles usually develop by late June to mid-July of the year following infestation, completing one beetle generation per year (Safranyik, 1988). The MPB attack develops in three stages, the first being the green attack, during which adult beetles lay eggs beneath the bark. 5  Larvae, hatching 1-2 weeks later, cut off the tree’s nutrient supply by feeding on the phloem (Taylor et al., 2006). The tree’s water and nutrient transport are also impacted by a beetle-introduced blue-stain fungus killing phloem and xylem cells (Safranyik et al., 1974). A year following infestation, in the red attack stage, the trees have died and their needles turn red. Decaying needles then fall off the trees during the last stage, the grey attack stage (BC Ministry of Forests, Lands and Natural Resource Operations1, 2017). Beetles generally attack trees with a larger diameter at breast height (DBH > 10 cm) and thicker bark, as a minimum bark thickness of ~1.5 mm is needed for gallery construction underneath the bark, where eggs are laid (Safranyik and Carroll, 2006). Tree diameter is also related to phloem thickness with a thicker phloem layer allowing higher brood production with larger beetles (Safranyik and Carroll, 2006). The most recent MPB outbreak in BC reached unprecedented extents. Since 1999, ~731 million m3 of commercial timber in an area of ~18.3 million ha have been killed by the beetle. The current infestation peaked in 2006 and is estimated to have killed ~55% of the mature merchantable lodgepole pine volume in BC by the predicted end of this outbreak in 2020 (BC Ministry of Forests, Lands and Natural Resource Operations, 2016). The extent of the attack in BC as well as the percentages of pine volume killed by 2015 can be seen in Figure 1.1. The rapid expansion of the beetle’s historic range into Northern BC, the Yukon, the Northwest Territories and east of the Rocky Mountains into Alberta during this outbreak (Cooke and Carroll, 2017) has raised concerns about its spread to Eastern Canada through boreal jack pine forests (Nealis and Cooke, 2014). In areas with lower pine volume than in central BC and western Alberta, such as central Alberta and Saskatchewan, climate suitability for MPB range expansion is assumed to be                                                  1 In 2017, the name was changed to BC Ministry of Forests, Lands, Natural Resource Operations & Rural Development. 6  moderate, but farther east, where pine volumes are higher around the Great Lakes, possibilities for the beetle to spread are higher (Safranyik et al., 2010). Thus, when weather conditions allow epidemic populations to develop, the MPB could become persistent outside its historic range (Safranyik et al., 2010). The frequency and extent of insect attacks is expected to increase with rising temperatures due to climate change (Williams et al., 2016).                  Mackenzie Prince George 7   Figure 1.1  Mountain pine beetle infestation map including observed percentage of pine volume killed in BC by 2015 based on the Provincial Aerial Overview of forest health from 1999 to 2015 and projections using the stochastic model BCMPB version 13 assuming no forest harvesting or beetle suppression activities. Source: https://www.for.gov.bc.ca/ftp/hre/external/!publish/web/bcmpb/year13/BCMPB.v13.2015Kill.pdf.  1.1.2 Effects of MPB attack on C and water balances The C balance of forests is impacted by MPB outbreaks, which kill trees, as the photosynthesis rate decreases in attacked stands resulting in lower C uptake (Brown et al., 2010). The net uptake of atmospheric CO2, i.e., the net ecosystem productivity (NEP) of attacked stands also depends on the response of R, which consists of autotrophic respiration (Ra) and heterotrophic respiration (Rh), to the attack. Ra is expected to decrease following beetle attack, but Rh might increase due to decaying foliage, roots, stems, branches and eventually snags. Inclusion of disturbance effects in C budget models is important to avoid overestimation of the forests' ability to offset anthropogenic CO2. The modelling study by Kurz et al. (2008) found that stands that were C sinks turned into C sources immediately after the outbreak. GPP declined more than R, leading to MPB-attacked stands in BC being a C source of 270 teragrams (Tg) between 2000 and 2020 with an additional loss of 50 Tg of C due to subsequent salvage harvesting (Kurz et al., 2008). However, according to Williams et al. (2016), increased atmospheric CO2 concentrations and forest resilience enhanced growth resulting in an increase in 8  C stocks of ~190 Tg year-1 of C in the US and of 10-30 Tg year-1 of C in Canada, despite harvesting, fire, windthrow, insect attacks and drought causing a decline in live biomass of ~200 Tg year-1 of C in the US. Arora et al. (2016) also predicted that the enhancement in C accumulation for 1999-2020 due to changing climate would exceed the beetle-induced loss of 328 Tg of C, leading to a sink of 900-1060 Tg of C depending on the climate change scenario. Using the Community Land Model (CLM) with prognostic C and nitrogen cycling (Version 4), Edburg et al. (2011) found MPB-attacked stands to have reduced stand-level net primary productivity (NPP) for several years following attack. They showed that C fluxes were impacted for decades and C stocks for up to 100 years depending on outbreak severity, delay in snagfall, snagfall rate, and management decisions following disturbance. Mortality rates impacted C fluxes and stocks in the long term, while both C fluxes and stocks were affected in the short term only by outbreak duration (Edburg et al., 2011). On a decadal scale, both GPP and R declined following tree mortality resulting in little change in NEP (Moore et al., 2013). Immediate impacts of MPB attack on C release might have been overestimated by previous studies, as Moore et al. (2013) found the increased amount of litter from needles and roots after the attack caused an increase in R for only 1-2 years, followed by a decline, when C inputs into the soil decreased. Stand structure and management following attack play an important role in determining NEP as indicated by eddy-covariance (EC) measurements conducted in MPB-attacked lodgepole pine-dominated stands in the central interior of BC. Two stands that were not invasively managed and left to naturally regenerate, reached C neutrality 3-5 years after attack due to an increased uptake by the remaining living trees and understory vegetation, whereas a clearcut in the same area remained a growing-season C source 10 years after harvesting (Brown et al., 9  2012). Thus, not harvesting stands with significant secondary structure immediately after attack might prevent them from remaining C sources for as many as 10 years (Brown et al., 2010). EC measurements in a southeast Wyoming beetle-attacked stand showed that C uptake and water use efficiency (WUE) were not impacted when canopy mortality increased from 30 to 78%, suggesting that the remaining trees and sub-canopy vegetation increased their C uptake and water use (Reed et al., 2014). The importance of broadleaf understory and younger coniferous trees in C uptake and evapotranspiration (E) a few years following beetle attack, was also shown using EC measurements (Emmel et al., 2013; Emmel et al., 2014), and component flux measurements using cuvettes and conifer chambers (Bowler et al., 2012), in the BC Interior. Mikkelson et al. (2013), using an integrated hydrologic model, showed that beetle attack also influences stand hydrology, causing reductions in E, less interception of precipitation (P) by the canopy after needle loss, higher snow accumulation, more radiation penetrating through the canopy resulting in earlier and faster snowmelt, and an increase in runoff. Average summertime E over an area of 170,000 km2 of affected forest following MPB attack in BC, determined using Moderate Resolution Spectroradiometer (MODIS) monthly data, decreased by 19% (Maness et al., 2012). This change and an increase in sensible heat flux of 8% may have secondary consequences for climate, such as impacts on atmospheric circulation, cloud cover and P (Maness et al., 2012). Vanderhoof and Williams (2015), also using MODIS data together with MPB outbreak location datasets and Global Modeling and Assimilation Office (GMAO) Modern Era Retrospective-Analysis for Research and Applications (MERRA) data, showed correlations between summertime E and stand density, leaf area index (LAI) and the fraction of absorbed photosynthetically active radiation (PAR). Lodgepole pine stands in the Colorado Rockies following MPB attack showed decreases in E of 6.3% 4-13 years after attack and of 18.7% 14-20 10  years after attack, but for stands 30-40 years after attack, E increased by 21.6% (Vanderhoof and Williams, 2015). Comparing MPB-attacked and mechanically girdled lodgepole pine forests in Colorado, USA, Hubbard et al. (2013) found that transpiration in beetle-attacked trees declined rapidly, while girdled trees showed lower transpiration rates than control trees, but continued to grow during the growing season following treatment, and beetle-attacked trees died due to the additional water supply reduction caused by blue stain fungi. Hubbard et al. (2013) concluded that the stand water balance will likely be influenced by the rapid decline in transpiration following beetle attack and that the surviving trees might benefit from an increase in water availability, which could also result in increased runoff and stream flow.   1.1.3 Forest Productivity Models There are many different models that can be used to simulate ecosystem processes. These models can be divided into six categories according to their applications (Pretzsch, 2009). These categories are (1) short-term, stand scale ecophysiological models, (2) forest growth and yield models, (3) long-term, management-oriented forest-growth and ecosystem dynamics models, (4) structural-functional individual tree models, (5) C budget models, and (6) process-based models with no stand-development component. Models in the first category address ecological processes (e.g., C sequestration, nutrient sustainability) and output biomass development for the stand scale over a few years. In this category, the PnET family of models (Aber and Federer, 1992) and BIOMASS (McMurtrie et al., 1990; McMurtrie and Landsberg, 1992), for example, calculate C and nitrogen balances. The second category includes forest growth models such as EFIMOD (Chertov and Komarov, 1997; Chertov et al., 1999 and 2001) and 3-PG (Landsberg and Waring, 1997), which estimate forest yield in response to environmental variables. They provide output 11  that is useful to forest managers, as they include parameters of interest in forestry such as tree height and stem volume at the stand level, and the models can be run for one generation of trees. Models in the third category are aimed at the long term, and predict forest development over several generations. Models such as FORECAST (Kimmins and Scoullar, 1994), FORCYTE (Kimmins and Scoullar, 1979), EASS-BEPS (Chen et al., 1999), C-CLASS (Arain et al., 2002), and ecosys (Grant et al., 1995) can be used to identify sustainable forest-management strategies. The fourth category includes models such as LIGNUM (Perttunen et al., 1996) and EMILION (Bosc, 2000), which simulate the functional and structural dynamics of trees. They can be applied in the management of individual trees in mixed forests under various site conditions (Pretzsch, 2009). The fifth category consists of C budget models such as CBM-CFS3 (Kurz et al., 2009), which calculate C stocks and stock changes over time. Models in the last category are process-based, short-time-scale models with no stand-development component, such as MAESTRO (Wang and Jarvis, 1990). In forestry, another category of models, empirical models that are based on relationships between inputs and outputs, such as site index curves providing yield estimates of stands (e.g., TIPSY (Table Interpolation Program for Stand Yields) (Mitchell et al., 2000)), are often used, although they are site specific and dependent on the conditions at the time of measurement. Thus, these models cannot be used to simulate the effects of climate change, such as rising atmospheric CO2 concentrations or temperature (Landsberg and Gower, 1997). Process-based models, which are based on the underlying physiological processes determining growth and the physical conditions influencing the trees, are useful in predicting the effects of climate change. To be of practical use, these models have to produce output that is relevant to forest managers and decision- and policy-makers. To achieve this, it is important that model developers work 12  together with the specific target group and keep their interests and capabilities to provide input variables in mind, creating computationally simple, longer time-step models including forest growth, which is relevant to forest managers, and which allow for landscape-level-assessments. Taking into consideration all of the above models, 3-PG was the most appropriate for my research, as it includes stand development, ecophysiological processes, forest yield estimates and is applicable for long-term simulations.   1.1.4 Physiological Principles in Predicting Growth Model (3-PG) 3-PG (use of Physiological Principles in Predicting Growth) is a process-based, stand growth model developed by Landsberg and Waring (1997). The four kinds of inputs required are (i) long-term average monthly or monthly weather data for each year, (ii) initial conditions (e.g., stand age, stocking, available volumetric soil water content (θa)), (iii) site factors (e.g., site latitude, species, soil class) and (iv) species-specific physiological parameter values. The model is run on a monthly time step. Biomass production, allocation of biomass between foliage, roots and stems, stem mortality, the soil water balance and stand characteristics of interest to forest managers are calculated in five sub-models. The output variables include E, NPP, specific leaf area (SLA), canopy LAI, mean annual volume increment (MAI) (since establishment) and mean DBH. The output can be monthly or annual values (Sands, 2004). In 3-PG, GPP is calculated using a light-use efficiency approach as the product of the canopy quantum efficiency (αc) and the intercepted solar radiation (Rs int), which is estimated using Beer's law. This is often referred to as a Monteith approach based on observations that biomass accumulated by a crop has a linear relationship to accumulated Rs (Monteith, 1977). In 3-PG, the environmental conditions reducing αc are taken into account by applying modifiers 13  such as air temperature (Ta), frost occurrence, vapour pressure deficit (VPD), θa, atmospheric CO2 concentration and stand age modifiers with values between 0 (most limiting) and 1 (not limiting). NPP is determined as a constant fraction (47%) of GPP. Growing conditions expressed by θa, VPD, and site fertility influence the allocation of NPP to roots, for example, resulting in a higher allocation of NPP to roots when trees are under stress. Stand age determines biomass allocation to foliage and stems with allocation to foliage decreasing with increasing DBH (Sands, 2004). Both age-dependent and density-dependent tree mortality are included in 3-PG, the latter being described by the -3/2 self-thinning law, which requires an upper limit to the mean single-tree stem mass for the current stem population, resulting in tree mortality when this limit is exceeded (Sands, 2004).  E is calculated using the Penman-Monteith equation (Monteith and Unsworth, 2008) with a fixed value of 0.2 m s-1 for the boundary layer conductance (ga) to avoid the need for wind speed data (Sands and Landsberg, 2002). Factors such as LAI, VPD, θa and the atmospheric CO2 concentration determine the canopy conductance (gc) from the maximum conductance (gcmax) (Landsberg and Sands, 2011). The soil water balance is modelled with the change in root-zone soil water storage (S) calculated as S = P – Dr – E      (1) where P is the precipitation and Dr is the drainage from the root zone. The water holding characteristics of the soil are important for the model as the maximum plant available water in the root zone (θa max) determines Dr, which occurs when θa exceeds θa max (Landsberg and Waring, 1997). The stand characteristics of interest to forest managers such as stem volume, mean stem diameter, basal area and MAI are calculated from the biomass pools and stand density (Sands, 2004).  14  3-PG has been applied to a range of different species which are listed in Table 1.1. It has been shown that 3-PG is successful in predicting the growth of Pinus radiata and Eucalyptus plantations (Landsberg and Waring, 1997). Sands and Landsberg (2002) showed that 3-PG performed well in modelling the growth of even-aged, intensively-managed, fertilized Eucalyptus globulus stands, but overestimated growth when site conditions such as drought were not represented by the long-term average climatic data. By including a temperature-dependent growth modifier reducing GPP at cool sites, an age-dependent SLA and fraction of stem biomass as branch and bark to enable calculation of stem volume or MAI, they improved the original 3-PG model from Landsberg and Waring (1997).   Table 1.1 List of species to which 3-PG has been applied and the respective studies. Species Studies Eucalyptus (Eucalyptus globulus, E. grandis, E. regnans, E. pilularis, E. maculate, E. obliqua) Landsberg and Waring, 1997; Coops et al., 1998; Sands and Landsberg, 2002; Sands, 2004 Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco var. glauca (Beissn.) Franco and Pseudotsuga menziesii (Mirb.) Franco var. menziesii) Coops et al., 2010 Grand fir (Abies grandis) Wei et al., 2014 Birch (Betula platyphylla) Potithep and Yasuoka, 2011 Pinus radiata Landsberg and Waring, 1997; Coops et al., 15  Species Studies 1998 Ponderosa pine (Pinus ponderosa) Coops et al., 2005 Lodgepole pine (Pinus contorta var. latifolia) Coops and Waring, 2011  In a study conducted by Wei et al. (2013), 3-PG was constrained with a new δ13C sub-model, which calculates δ13C in tree rings. They tested the modified 3-PG model using data from two 80-year-old grand fir stands in northern Idaho using as many independent measurements as possible to parameterize 3-PG. The δ13C sub-model was validated using tree-ring δ13C measurements which were then used to parameterize 3-PG. This new sub-model assesses the coupling between the C and water balances by way of GPP and gc. Thus, the δ13C sub-model can be used to make sure that the correct NPP value is passed to the C allocation algorithms and reduce error propagation due to incorrect parameters in the simulation of GPP and transpiration (Wei et al., 2013). Coops et al. (1998) developed a spatial version of 3-PG, called 3-PGS, which uses remotely-sensed data as an input and has been used to study forest productivity across landscapes (Sands and Landsberg, 2002). The fraction of PAR absorbed by forest canopies is estimated from satellite observations, as it is related to the normalized difference vegetation index (NDVI), which is the normalized difference between the reflectances measured in the near-infrared and red wavelengths. Thus, absorbed PAR can be computed from NDVI and Rs measurements (Coops et al., 1998). Estimates of mean annual growth made with 3-PGS over large areas (8 km x 8 km) corresponded well with growth rates derived from the average of forest plot measurements (Coops et al., 1998).   16  1.2 Research objectives The objectives of this research were to synthesize the results of micrometeorological (i.e., EC) flux measurements conducted in MPB-attacked stands in the BC Interior between Prince George and Mackenzie over a decade, and to simulate the C and water balances of an MPB-attacked lodgepole pine stand using 3-PG. The locations of the sites, where EC measurements were made, can be seen in Figure 1.2.        17  Figure 1.2  Locations of the MPB sites in lodgepole pine-dominated stands attacked in 2006 (MPB-06), 2003 (MPB-03) and 2005 and 2006 (MPB-09), respectively, where eddy-covariance flux measurements have been made. MPB-06 and MPB-03 were not harvested, MPB-09 was partial-salvage-harvested in 2009 while MPB-09C, nearby MPB-09, was clearcut-salvage-harvested. CC-05 and CC-97 were clearcuts nearby MPB-06 harvested in 2005 and 1997, respectively.   Specific objectives were to: 1. determine the effects of MPB attack on growing season and annual C and water fluxes in the three main stands, 2. determine the rate of recovery and assess the impacts of different management strategies, specifically not harvesting, partial- and clearcut-salvage harvesting, on the C and water fluxes of lodgepole pine stands recovering from MPB attack, 3. modify 3-PG to simulate the impact of insect attack on C and water balances including water availability from snowmelt, changes in light transmission through the canopy due to canopy mortality following attack as well as treefall, and an Rh sub-model enabling the calculation of NEP, 4. validate the modified model using measured E, NEP, EC-estimated GPP, R and WUE in an 80-year-old lodgepole pine stand (MPB-06) during the first decade following MPB attack, and 5. forecast the effects of beetle attack on E, GPP, R, NEP and WUE for MPB-06 over the next decade.  18  1.3 Thesis overview This dissertation is structured as five chapters, beginning with an introduction including the background and current state of knowledge. This is followed by three main research chapters, which were formulated as the research evolved but are presented as stand-alone documents with their own dedicated introduction, review of literature, methods, results, and discussion sections. Chapter 2 presents measurements and simulations of the water balance and water use efficiency of an MPB-attacked lodgepole pine stand in the BC Interior (MPB-06) using a modified version of 3-PG. The modifications made to 3-PG to enable it to be applied to MPB-attacked stands are explained in this and the following chapter. Chapter 3 focuses on the C balance of an MPB-attacked stand (MPB-06) and presents modelling and measurement results for R, including Rh and Ra, GPP and NEP. Chapter 4 synthesizes the effects of beetle attack on growing season and annual C and water fluxes in the three main stands and nearby clearcuts, and shows the impact of different management strategies on NEP, E, Dr and S. In Chapter 5, the most important findings and conclusions are summarized, limitations of the research are discussed and an outlook on possible future work is presented. The appendices show photographs of the sites and instrumentation.  19  Chapter 2: Measurements and simulations using the 3-PG model of the water balance and water use efficiency of a lodgepole pine stand following mountain pine beetle attack  2.1 Introduction The MPB outbreak, which started in the late 1990s, has killed ~728 million m3 (54%) of merchantable lodgepole pine in BC. The infestation is expected to subside by 2017 (BC Ministry of Forests, Lands and Natural Resource Operations, 2015). Tree mortality due to this outbreak has impacted hydrology through increased snow accumulation due to less interception than in healthy stands, combined with earlier and more rapid snowmelt as well as reduction in late-summer transpiration (Mikkelson et al., 2013). Winkler et al. (2015) found increases in snow accumulation in MPB-attacked forests in BC. Differences in snow water equivalent on April 1 between clearcuts and forests decreased in the red- and especially the grey-attack stages due to canopy loss, which showed that changes in snow accumulation lagged behind beetle attack (Winkler et al., 2015). Thus, clearcut harvesting following attack would impact hydrology more immediately than is observed in not-harvested stands. Biederman et al. (2015) found no changes in stream flow during the decade following tree die-off due to MPB attack in Colorado, and attributed this to increased E by surviving vegetation, and increased snow sublimation and evaporation from the understory. Reed et al. (2014), using EC measurements in a bark-beetle infested lodgepole pine stand in Southeast Wyoming, found no change in CO2 uptake and WUE as tree mortality increased from 30 to 78% over three years indicating increasing CO2 uptake and 20  water use by the remaining trees and sub-canopy vegetation. Emmel et al. (2013), also using EC measurements, showed that summertime E in an MPB-attacked forest in interior BC primarily originated from the understory layer, and that the overstory was not a significant source of water vapour in the late summer. Several studies have examined E in attacked stands using satellite data. In BC, Maness et al. (2012) found that summertime E in affected forests was reduced, on average, by 19% during the first few years following attack. However, this study was based on remote sensing data for a large affected area containing a variety of forest types. Vanderhoof and Williams (2015), however, showed that the effect of MPB attack on summer E of lodgepole pine stands in the Colorado Rockies varied over time. They found that in stands 14-20 years after attack, summer E was reduced by 18.7%, but in stands 30-40 years after attack, E increased by 21.6% relative to non-attacked stands. These changes in E were correlated with changes in stand density, leaf area index (LAI) and the fraction of absorbed PAR (fPAR). Hubbard et al. (2013) compared the effects of MPB attack and mechanical girdling on lodgepole pine forests in Colorado, USA, and found that mortality in beetle-attacked trees was primarily caused by blue stain fungi rather than phloem feeding by beetles. The fungi reduced the water flow from the soil to the canopy and thus played a major role in killing the trees. The rapid decline in transpiration following MPB attack would likely influence the stand water balance. Thus, an increase in available water might benefit the surviving trees, or result in increased runoff to stream flow (Hubbard et al., 2013). Modelling the current MPB attack on lodgepole pine stands in BC using the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3), Kurz et al. (2008) estimated its impact between 2000 and 2020 to be a loss of 270 Tg of C with an additional 50 Tg of C loss due to salvage harvesting beetle-killed and living trees. Williams et al. (2016), reviewing disturbance 21  effects on the C balance of US forests, found that harvesting, fire, windthrow, insect attacks and droughts together decreased live biomass by about 200 Tg year-1 of C across the conterminous US. However, due to forest resilience and enhanced growth, C stocks increased by about 190 Tg year-1 of C in the US at the same time and by about 10-30 Tg year-1 of C in Canada. Insect attacks are expected to become more frequent and extensive due to climate change (i.e., increasing temperatures), which together with management impacts disturbance rates and intensities as well as forest regeneration (Williams et al., 2016). Attacked stands are managed following three different approaches: no harvesting, partial-salvage harvesting and clearcut-salvage harvesting (Coates and Hall, 2005). Edburg et al. (2011), using the CLM with prognostic C and nitrogen cycling (Version 4), studied the impact of several outbreak characteristics, such as severity and duration of insect outbreaks, on C and nitrogen dynamics following different bark beetle outbreak scenarios. They found that while outbreak severity, delay in snagfall, snagfall rate and management decisions all impact C fluxes for decades and C stocks up to 100 years following the attack, outbreak duration only played a major role for C fluxes in the short term. The number of trees killed, on the other hand, impacted C fluxes and stocks in the long term (Edburg et al., 2011). We have been studying the impacts of MPB attack on water and C balances of lodgepole pine forests in northern interior BC, between Prince George and Mackenzie, using EC flux measurements at three tower sites complemented with short-term EC measurements in nearby clearcuts since 2006 (Brown et al., 2010; Mathys et al., 2013; Emmel et al., 2014). In this study, we use direct EC measurements from a flux tower, and modify the ecosystem model 3-PG (use of Physiological Principles in Predicting Growth) (Landsberg and Waring, 1997) to quantify the changes in E and GPP and to study the long-term recovery of water and C balances of a BC 22  interior lodgepole pine stand following MPB attack. 3-PG includes stand development, ecophysiological processes, forest yield estimates, and has been used for long-term simulations in different species such as Eucalyptus (Eucalyptus globulus, E. grandis, E. regnans, E. pilularis, E. maculate, E. obliqua) (Coops et al., 1998; Sands and Landsberg, 2002; Sands, 2004), Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco var. glauca (Beissn.) Franco and Pseudotsuga menziesii (Mirb.) Franco var. menziesii) (Coops et al., 2010), grand fir (Abies grandis) (Wei et al., 2014), birch (Betula platyphylla) (Potithep and Yasuoka, 2011), Pinus radiata (Coops et al., 1998), ponderosa pine (Pinus ponderosa) (Coops et al., 2005) and lodgepole pine (Pinus contorta Dougl.) (Coops and Waring, 2011). Sands and Landsberg (2002) used 3-PG for even-aged, intensively-managed, fertilized stands of Eucalyptus globulus, and showed that it predicted growth and development well at most of the sites. The model overestimated growth in limiting environments such as occurrence of drought or frost damage, which were not accurately represented by the long-term average climate data. However, they did not attempt to parameterize the water balance. Their changes to the original 3-PG model from Landsberg and Waring (1997) included a temperature-dependent growth modifier, which reduces productivity at cool sites relative to warm sites, an age-dependent SLA, and an age-dependent fraction of stem biomass (ws) as branches and bark, which allows the conversion of predicted stem biomass into more practical variables, such as stem volume (excluding bark and branches) or MAI (Sands and Landsberg, 2002). Almeida and Sands (2016) modified the water balance sub-model in 3-PG by dividing the soil profile into two soil layers, separating E into canopy transpiration and soil evaporation, and using time steps shorter than a month (i.e., a daily time step). They found that it markedly improved their estimates of volumetric soil water content (), a and their temporal variation. However, we used the original 3-PG water balance sub-model, as 23  their new sub-model required additional data, e.g., specific root volume, canopy transpiration and evaporation from the soil surface. Also, in many applications, requirement of daily data would be demanding for the stakeholders whereas monthly data are more easily available. For its application to MPB-attacked lodgepole pine forests, it was necessary to modify 3-PG to simulate an MPB-induced overstory mortality and, as a result, increased understory growth due to higher resource availability. E, LAI, GPP and NPP must be calculated for both overstory and understory components. This requires that changes in light transmission through the canopy due to overstory mortality following attack must be taken into account in the model. Specifically, the objectives of this paper were to (i) modify 3-PG to include the effects of insect attack including water availability from altered snow accumulation and snowmelt, (ii) validate the modified model using measured E, EC-estimated GPP and WUE during the first nine years following MPB attack in an 80-year-old lodgepole pine stand, and (iii) forecast the effects of beetle attack on E, GPP and WUE over the next decade. Here we concentrate on the hydrology and the effects of beetle attack on WUE. A subsequent paper will report on the C balance modelling results.  2.2 Methods 2.2.1 The model 3-PG is a process-based stand growth model developed by Landsberg and Waring (1997). The model is run on a monthly time step and consists of five sub-models, which calculate biomass production, allocation of biomass between foliage, roots and stems, stem mortality, soil water balance and stand characteristics that are of interest to forest managers. It requires four 24  kinds of inputs: initial conditions (e.g., stand age, stocking, a in the 0-60-cm layer), site factors (e.g., site latitude, species, soil class), weather data (incident solar radiation (Rs), Ta, VPD, P), and physiological parameter values (see Table 2.1) characterizing the species to be modelled. The output variables include E, NPP, SLA, LAI, MAI, and mean DBH (Sands, 2004). Figure 2.1 shows the basic structure of 3-PG and the causal effects of its variables and processes with the modifications made in this study.   Table 2.1  List of the 3-PG parameters and their values for lodgepole pine (Pinus contorta var. latifolia) used in this study. 3-PG parameter Name Units Lodgepole pine Sourcea 1. Biomass partitioning and turnover       Allometric relationships & partitioning       Foliage: stem partitioning ratio @ VPD = 2 cm pFS2 - 1.53 C & W Foliage: stem partitioning ratio @ VPD = 20 cm pFS20 - 0.55 C & W Constant in the stem mass vs. diam. relationship aS - 0.0073 C & W Power in the stem mass vs. diam. relationship nS - 3.282 C & W Maximum fraction of NPP to roots pRx - 0.8 C & W Minimum fraction of NPP to roots pRn - 0.25 C & W Litterfall & root turnover       Maximum litterfall rate gammaFx month-1 0.015 C & W Litterfall rate at t = 0 gammaF0 month-1 0.01 C & W Age at which litterfall rate has median value tgammaF months 120 C & W Average monthly root turnover rate gammaR month-1 0.015 C & W 2.  NPP & conductance modifiers       Temperature modifier (fT)       Minimum temperature for growth Tmin deg. C -7 C & W Optimum temperature for growth Topt deg. C 15 C & W Maximum temperature for growth Tmax deg. C 30 C & W Frost modifier (fFrost)       Days production lost per frost day kF days 1 Law et al. Soil water modifier (fSW)       Moisture ratio deficit for fθ = 0.5  SWconst - 0.6 Sands 25  3-PG parameter Name Units Lodgepole pine Sourcea Power of moisture ratio deficit SWpower - 2 Sands Atmospheric CO2 modifier (fCO2)       Assimilation enhancement factor at 700 ppm fCalpha700 - 1.4 S & L Canopy conductance enhancement factor at 700 ppm fCg700 - 0.7 S & L Fertility effects       Value of 'm' when FR = 0 m0 - 0 Law et al. Value of 'fNutr' when FR = 0 fN0 - 0.6 R & M Power of (1-FR) in 'fNutr'  fNn - 1 R & M Age modifier (fAge)       Maximum stand age used in age modifier MaxAge/ ymax years 150 C & W Power of relative age in function for fAge nAge - 4 C & W Relative age to give fAge = 0.5 rAge - 0.95 C & W Stem mortality & self-thinning       Mortality rate for large t gammaNx % year-1 0 S & L Seedling mortality rate (t = 0) gammaN0 % year-1 0 S & L Age at which mortality rate has median value tgammaN years 0 S & L Shape of mortality response ngammaN - 1 S & L Stand age at time of beetle attackc attackStAge years 79.3 Brown et al. Initial rate of mortality at time of beetle attackc γn0 % month-1 31.6 This study Time constant for MPB induced tree mortalityc τt months 5.62 This study Time constant for needle fallc τn months 28 This study Max. stem mass per tree @ 1000 trees/hectare wSx1000 kg tree-1 220 C & W Max. stem mass per tree for understory @ 1000 trees/hectare wSx1000u kg tree-1 220 C & Wb Power in self-thinning rule thinPower - 1.5 S & L Fraction mean single-tree foliage biomass lost per dead treec mF - 1 This study Fraction mean single-tree root biomass lost per dead treec mR - 1 This study Fraction mean single-tree stem biomass lost per dead treec mS - 1 This study Fraction mean single-tree foliage biomass lost per dead tree in the understory mFu - 0 S & L Fraction mean single-tree root biomass lost per dead tree in the understory mRu - 0.2 S & L Fraction mean single-tree stem biomass lost per dead tree in the understory mSu - 0.2 S & L 3.  Canopy structure and processes       Specific leaf area       Specific leaf area at age 0 SLA0 m2 kg-1 5.79 Bowler et al. Specific leaf area for mature leaves SLA1 m2 kg-1 3.1 C & W Age at which specific leaf area = (SLA0+SLA1)/2 tSLA years 2 R & M 26  3-PG parameter Name Units Lodgepole pine Sourcea Light interception       Extinction coefficient for absorption of PAR by canopy k - 0.6 This study Extinction coefficient for absorption of PAR by understory ku - 0.6 This study b Age at canopy cover  fullCanAge years 20 This study Maximum proportion of rainfall evaporated from canopy MaxIntcptn - 0.15 S & L LAI for maximum rainfall interception LAImaxIntcptn m2 m-2 3.1 This study Production and respiration       Canopy quantum efficiency alpha molC molPAR-1 0.047 Law et al. Understory quantum efficiency alphau molC molPAR-1 0.03 Law et al. Maximum canopy cover CanCoverMax - 0.71 This study Ratio NPP/GPP Y - 0.47 C & W Heterotrophic respiration       Base-respiration rate at 10 C R10 g m-2 month-1 of C 72.94 Brown et al. 2 Q10 for heterotrophic respirationc Q10   2.8 Brown et al. 2 Conductance       Minimum canopy conductance MinCond m s-1 0 S & L Maximum canopy conductance MaxCond m s-1 0.02 S & L LAI for maximum canopy conductance LAIgcx m2 m-2 3 C & W Maximum canopy conductance for understory MaxCondu m s-1 0.0144 C & W LAI for maximum conductance for understory LAIgcxu m2 m-2 3 C & W Defines stomatal response to VPD CoeffCond/sD mbar-1 0.055 This study Canopy boundary layer conductance BLcond m s-1 0.2 S & L 4.  Wood and stand properties       Branch and bark fraction (fracBB)       Branch and bark fraction at age 0 fracBB0 - 0.5 R & M Branch and bark fraction for mature stands fracBB1 - 0.1 R & M Age at which fracBB = (fracBB0+fracBB1)/2 tBB years 5 R & M Basic Density       Minimum basic density - for young trees rhoMin kg m-3 480 R & M Maximum basic density - for older trees rhoMax kg m-3 480 R & M Age at which rho = (rhoMin+rhoMax)/2 tRho years 4 S & L Stem height       Constant in the stem height relationship aH - 0.6 This study Power of DBH in the stem height relationship nHB - 1.2 This study Power of stocking in the stem height relationship nHN - 0 S & L Stem volume       Constant in the stem volume relationship aV - 0 S & L Power of DBH in the stem volume relationship nVB - 0 S & L 27  3-PG parameter Name Units Lodgepole pine Sourcea Power of stocking in the stem volume relationship nVN - 0 S & L Conversion factors       Intercept of net v. solar radiation relationship Qa W m-2 -60 This study Slope of net v. solar radiation relationship Qb - 0.6 This study Molecular weight of dry matter gDM_mol g mol-1 of DM 24 S & L Conversion of solar radiation to PAR molPAR_MJ mol MJ-1 1.95 This study a C & W stands for Coops and Waring (2011), S & L stands for Sands and Landsberg (2002), Law et al. stands for Law et al. (2000), R & M stands for Raison and Myers (1992), Sands stands for Sands (2010), Bowler et al. stands for Bowler et al. (2012), Brown et al. stands for Brown et al. (2010) and Brown et al. 2 stands for Brown et al. (2012). b This assumes that the same characteristics apply to young trees as to older lodgepole pine trees. c Denotes parameters associated with the modifications of 3-PG for MPB-attacked lodgepole pine stands.  28   Figure 2.1  Basic structure of 3-PG including modifications made (in red) for its application to MPB-attacked lodgepole pine stands. The variables involved in the water balance are shown in blue and the allocation of biomass is shown in green. The subscripts ‘o’ and ‘u’ stand for the overstory and understory, respectively. α is the quantum efficiency of photosynthesis, φ is the physiological modifier applied to the maximum canopy conductance (gcmax) and α, ws is the average tree stem biomass and γn is the tree mortality rate. The visual basic code for the modified version of 3-PG is available on request (gesa.meyer@ubc.ca).  Intercepted Rs (Rs int) is obtained from Rs (MJ m-2 month-1) using Beer’s law:   Rs int = Rs (1 - exp(-k L))     (2)  29  where L (m2 m-2) is the canopy LAI and k is the canopy Rs extinction coefficient. GPP (g m-2 month-1 of C) is determined using a light-use efficiency approach, and is assumed to be proportional to Rs int (MJ m-2 month-1) with the canopy quantum efficiency α (mol of C (mol)-1 of PAR photons) being the proportionality factor, the value of which depends on temperature, VPD, atmospheric CO2 concentration and soil fertility:  GPP = α Rs int  12 g mol-1 of C  1.95 mol of PAR photons MJ-1.  (3)  The last term on the right hand side of Equation (3) converts Rs int into PAR obtained from PAR and Rs measurements at the site in this study. NPP is calculated as a constant fraction (47%) of GPP (Coops and Waring, 2011; Sands, 2004; Campioli et al., 2015).  In the model, a, VPD, and site fertility influence the allocation of NPP to roots. If trees are under stress, e.g., low nutritional status or low a, more NPP is allocated to roots. Biomass allocation to foliage and stems depends on stand age. With larger DBH, allocation to stems increases and allocation to foliage decreases (Sands, 2004). Stem mortality is implemented by using two functions, an age-dependent mortality function and a density-dependent mortality function described by the -3/2 self-thinning law (Drew and Flewelling, 1977). The density-dependent mortality function requires an upper limit to the mean single-tree stem mass for the current stem population as an input. Should this limit be exceeded, trees die and eventually fall down, but remain in the stand, so that the average live tree stem mass does not exceed the maximum mean stem mass. It is assumed that little foliage is lost due to mortality as smaller trees with little foliage usually die first (Sands, 2004). This can, however, be adjusted to the situation, e.g., in the case of an MPB attack, the older and larger trees die first. The soil water balance is modelled using a monthly time step with S (mm) calculated using 30  S = P – Dr – E      (4) where P is in mm, Dr is the drainage from the root zone (mm) and E (mm) is calculated using the Penman-Monteith equation (Landsberg and Gower, 1997; Monteith and Unsworth, 2013):                        ( )    ( 1( ))ccn aag DgR G gEg                                              (5) where λE (W m-2) is the latent heat flux, Rn (W m-2) is the net radiation, λ (J kg-1) is the latent heat of vaporization, ρ (kg m-3) is the air density, ε is the dimensionless rate of change of saturated vapour pressure with respect to Ta (equal to the slope of the saturated vapour pressure curve divided by the psychrometric constant), gc (m s-1) is the canopy conductance, ga (m s-1) is the boundary layer conductance, and D is the specific humidity deficit of canopy air (g of H2O (g of air)-1). D is obtained from the VPD (kPa) by dividing by atmospheric pressure and multiplying by the ratio of the molar mass of water vapour to that of air. ga, which Massman (1999) has shown can be parameterized using LAI and u*, is fixed in 3-PG with a value of 0.2 m s-1 to avoid the need for wind speed data (Sands and Landsberg, 2002; Coops and Waring, 2011). The soil heat flux (G) in Equation (5) is neglected because it is very small relative to Rn at the monthly time scale. The parameter values in Equation (5) used in this study, and their sources, are given in Table 2.1. gc is calculated as gcmax limited by factors including LAI and a as follows:                        max min maxminmin(  ( ) , )   c c c CggccxcLg g g g g fL                   (6) where gcmin (m s-1) is the minimum value of gc, Lgcx (m2 m-2) is the LAI for gcmax (m s-1), fCg is the gc modifier for atmospheric CO2 concentration, and φ is a physiological modifier, which combines the responses of gc to VPD and a, expressed by VPD and a modifiers with values 31  between 0 (most limiting) and 1 (not limiting). fCg equals 1 when the atmospheric CO2 concentration is 350 ppm, and decreases hyperbolically to 0 with increasing CO2 concentration (Landsberg and Sands, 2011). In our case, however, L never exceeded Lgcx. The stand characteristics of interest to forest managers such as stem volume, mean DBH, basal area and MAI are calculated from the biomass pools and stand density (Sands, 2004). a at the beginning of the model run (a init) is determined as  - WP, where WP is  at the permanent wilting point (WP) of the soil. a init is constrained by the minimum (a min) (i.e., zero in our case) and maximum plant available water (a max) in the root zone (0-60-cm layer) with a min ≤ a init ≤ a max. a max corresponds to the difference between  at field capacity (FC) and WP (a max = FC - WP). If a exceeds a max, the excess water is lost as Dr (Landsberg and Waring, 1997).  2.2.2 Modifications to the model In order to apply the model to a mid-latitude MPB-attacked stand, several modifications were required. Since 3-PG has mainly been used in warmer regions, where snowfall does not occur, we modified the water balance sub-model for use in MPB-affected areas by calculating P from measurements of rainfall and snowfall (see Section 2.2.3). The model was driven by the liquid water input, calculated as the sum of rainfall and snow water equivalent in a given month (calculated from snowfall assuming a density of 100 kg m-3). If the monthly mean Ta, rises above 0oC, snowmelt is initiated. The amount of snow melted in each month depends on the monthly mean temperature and is the minimum of the accumulated snow and the potential snowmelt, which is calculated as Ta (C)  melt factor  number of days in the month with Ta > 0C, where the melt factor is 2 mm d-1 C-1 (Buttle, 2009). At our measurement site, all the snow is melted 32  by the end of May because mean monthly Ta for May varies between 8C and 12C. The inclusion of snow in the modified version of 3-PG focused on the liquid water availability due to snowmelt. Other effects of snow cover such as changes in albedo and sublimation of snow occurring during the winter instead of soil evaporation, were not considered in the model and could be included in future versions. The modified version of 3-PG treats the growing understory (lodgepole pine saplings and seedlings) and dying overstory separately, so that stand-level LAI, GPP, NPP, E and Ra are calculated as the sum of the values for the two vegetation layers. The same equations as in the original 3-PG (Equations (2), (3), (5) and (6)) are used for the overstory and understory with their respective parameters and initial values. In addition to the age- and density-dependent mortality functions included in the original model, we defined an insect-induced mortality function with an exponentially declining mortality rate (γn, % month-1):  γn = γn0 exp(-t/ τt) for t ≥ 0    (7)  where t is the time since attack (months), γn0 is the initial rate of mortality, and τt is the time constant for MPB induced tree mortality. The two parameters were obtained from a statistical fit to annual mortality observations (see Section 2.2.3). The decrease in stand density then results in a lower live overstory LAI. The overstory and understory are assumed to have similar root zone depths with the same access to available water. Radiation reaching the understory is calculated using Rs exp(-k L) where k and L are the extinction coefficient (see Table 2.1) and LAI of the overstory, respectively. Total overstory LAI is calculated as the sum of live and dead needles. Arora et al. (2016) used a half-life of 0.35 years (i.e., time constant (τn) = 0.50 years) for the dead needles, which means that the trees have lost close to 90% of their needles within one year after attack. The proportion of dead needles at MPB-06 increased rapidly in the first two years 33  following the beginning of the attack and was calculated from the decrease in live overstory LAI as the stand progressed to the peak of the attack. If the live LAI decreased between the previous and current month, dead LAI was accumulated. Needle fall was assumed to begin two years after the beginning of the beetle attack, as observed in photographs of the stand (no needle fall was measured in this study), and was approximated by an exponential function as follows:   Ld t = Ld exp(-t/τn)      (8)  where Ld t (m2 m-2) is the current dead LAI, Ld (m2 m-2) is the accumulated dead LAI on dying trees before significant numbers of these needles start falling, t is the time since attack and τn is the time constant for needle fall, which, based on observations of the timing of needle fall in this and a nearby stand, was ~2.33 years. We also calculated water deficit (WD) to understand how it varies as the stand recovers from the MPB attack. We estimated potential evapotranspiration (Ep) in the same way as the actual E (Equation (5)) but without soil water limitation. For this purpose, φ in Equation (5) only takes into account the limiting effect of VPD on gc as the soil water modifier was set to 1. WD was obtained as Ep – E. Following Mellander et al. (2004), who showed that stomatal conductance and thus transpiration in a young Scots pine stand was reduced during spring when soil temperature (Ts) was below ~8 C, we introduced a Ta-limitation on gc in the months of March to June as gc limited = min(0.06  Ta, 1)  gc (where Ta is in C).  Further modifications to the model to enable the calculation of NEP, which requires a Rh sub-model, were made, but will be described in detail in a subsequent paper on the effect of beetle attack on the C balance.  34  2.2.3 Measurements To validate the model, we used nine years (2007 to 2015) of EC flux and climate measurements in an 80-year-old (in 2007) MPB-attacked lodgepole pine stand in the BC interior Sub-Boreal Spruce biogeoclimatic zone (Meidinger and Pojar, 1991). The site (MPB-06), described in detail in Brown et al. (2010), is located approximately 35 km southeast of Mackenzie at Kennedy Siding (550642.8 N, 1225028.5 W) where a 32-m-tall scaffold tower (~2.1 m long x ~1.5 m wide) was erected in July 2006. This almost pure 15-m-tall (in 2007) lodgepole pine stand with an understory consisting mainly of pine seedlings and saplings on a gravelly sandy loam soil was attacked in 2006 and not salvage harvested. As shown in Brown et al. (2012), in August 2006 about 50% of the stand was in the green-attack stage with the percentage of trees attacked increasing with time as the beetle continued its attack on this stand so that by August 2012 only 14% of the canopy trees were healthy (D. Seip, BC Ministry of Environment, Prince George, personal communication, 2016). Since then no further infestation has been detected and the number of live trees remained at 14% of the pre-attack value (see Figure 2.2). Flux and climate measurements are described in detail in Brown et al. (2010). Briefly, water vapour, sensible heat and CO2 fluxes were measured using an EC system comprising a three-dimensional ultrasonic anemometer (model CSAT3, Campbell Scientific Inc. (CSI), Logan, Utah) and an open-path infrared gas analyser (IRGA) (model LI-7500, LI-COR, Inc., Lincoln, Nebraska) installed on the tower above the canopy at a height of 26 m. Climate measurements included above-canopy up-welling and down-welling shortwave and longwave radiation (net radiometer, model CNR1, Kipp and Zonen B.V., Delft, The Netherlands) at the 30-m height, and above-canopy up-welling and down-welling (30-m height), and below-canopy 35  down-welling photosynthetic photon flux density (quantum sensor, model LI-190 AS, LI-COR Inc.) at the 3-m height. Rainfall was measured at canopy height with a tipping bucket rain gauge. As the rain gauge was not heated, it did not measure snowfall reliably, so wintertime P was estimated from snowfall data. Ta was measured at the 25-m height. Ts was measured at depths of 5, 10, 20 and 50 cm, and  was measured at the 0-10 cm and 30-50 cm depths. Snow depth was measured acoustically (model SR50A, CSI) in a nearby clearcut and snowfall and its liquid water equivalent was calculated from snow depth and manual measurements of snow water equivalent. Also measured were wind speed, relative humidity (RH) and G. Climate measurements required to run the model from 1 Jan 2005 to 15 July 2006 were obtained from the Mackenzie Airport, supplemented by data from Prince George Airport 130 km south of the site. Comparison of climate data from these airports with measurements at MPB-06 showed good agreement with coefficients of determination, slopes of the regressions and root mean square errors (RMSEs) for Ta of 96%, 0.90 and 2.21 C for July 2006 to February 2007 at Mackenzie Airport, and 99%, 1.02 and 1.08 C for July 2006 to October 2009 at Prince George Airport, respectively. Rs and VPD for 2005-2006 were calculated using linear regression relationships with Ta obtained from measured values from 2007 – 2015, as these variables were not measured at these airports. To test the model sensitivity to the climate data for 2005 – 2006, we also used different climate scenarios as model input for those years when the climate variables were not measured directly at the site. In addition to the procedure of obtaining climate data for 2005 - 2006 described above, we also used the climate scenarios: (1) averages of Rs and VPD measured from August 2006 to December 2015, (2) averages of Rs and VPD measured from January 2007 to December 2015, (3) 2007 and 2008 measurements of Rs and VPD, (4) averages of all climate variables 36  required in 3-PG (Rs, Ta, VPD, P), measured from January 2007 to December 2015, and (5) 2007 and 2008 measurements of all the climate variables.    Figure 2.2  Calculated monthly live overstory stem density at MPB-06 using an initial stand density value (900 stems ha-1) and γn from Equation (7) for 2005 – 2015 (curve) (γn0 = 31.6 % month-1 and t = 5.6 months) and observed stand density (circles).  For the water balance, Dr was obtained using Equation (4) with measured values of P, E and S. Canopy LAI was measured using three optical approaches including hemispherical photography done once a year between 2007 and 2015 (Brown et al., 2010). Hemispherical photographs were taken with a camera with a fisheye converter lens at 1.2-m height in 24 37  locations along 150-m long transects in four directions (NW, NE, SW and SE) from the tower (Foord and Spittlehouse, 2014). Net ecosystem exchange of CO2 (NEE) was calculated half-hourly by adding measured CO2 fluxes to changes in air-column CO2 storage (below the EC sensors). Following Brown et al. (2010) and Morgenstern et al. (2004), the difference between half-hourly averaged CO2 mixing ratios measured at the 26-m height for the following and previous half-hours divided by the difference in time between the measurements and multiplied by the dry air density determined the change in CO2 storage in the whole air column below the EC measurement height (Hollinger et al., 1994). NEP was obtained as –NEE. GPP was then calculated by adding daytime NEP to daytime R where the latter was calculated by extrapolating the nighttime relationship between R (EC-measured NEE) and Ts at the 5-cm depth to daytime (Brown et al., 2010).  Quality control and data analysis details are described in Brown et al. (2010, 2014). Briefly, if more than 30% of values of a trace in a 30-minute period were bad measurements, for example, due to instrument malfunction during rainfall as indicated by an instrument diagnostic warning flag, EC flux data for that half hour were rejected. EC data were also excluded, when water vapour and CO2 mixing ratios measured by the IRGA were beyond bounds set by minimum and maximum values. To avoid times with insufficient turbulent mixing, nighttime EC data when the friction velocity (u*) was less than the u* threshold (u*th) of 0.30 m s-1 (Brown et al., 2010) were excluded. Burba et al. (2008) discovered issues with open-path IRGAs, where wintertime CO2 uptake occurred due to heating of the instrument in cold air conditions. To exclude data with these issues, wintertime (defined as times when Ts at the 5-cm depth < 1 oC and Ta < 5 oC) fluxes with negative NEE were removed. As the occurrence of this instrument malfunction was related to low wind speeds (Burba et al., 2008), we also rejected daytime data 38  when wind speed was < 4 m s-1 during the winter (Brown et al., 2010). In-house software written in MATLAB® (Version 7.5.0, The Math Works Inc., Natick, MA, USA) was used to perform all data calculations and analysis for this paper.  2.2.4 Uncertainty analysis of EC fluxes The uncertainty in annual NEP was obtained as the sum of the systematic error due to the selection of the u*th and the random error with the latter estimated as the square root of the sum of squares of the random error associated with the half-hourly NEP measurements and uncertainties in the gap-filling procedure (Krishnan et al., 2006). The random error of annual NEP was determined by assigning a 20% random error to each half-hourly value of NEP (Morgenstern et al., 2004). Uncertainty in gap-filling was obtained using a Monte Carlo simulation where gaps in measured NEP ranging from one half-hour to 10 days were generated using a uniformly distributed random number generator. Modelled values, using R versus Ts and GPP versus PAR (Michaelis-Menten light response equation) relationships, were used to gap fill following Krishnan et al. (2006). To determine the 95% confidence interval, this was done 500 times. To estimate the systematic error associated with the u*th selection, NEP was recalculated with u*th varied by up to 50% (Morgenstern et al., 2004). The same procedure was used to estimate uncertainties in annual GPP, R and E. The latter was gap-filled using the procedure of Brown et al. (2014). Other uncertainties arising from possible underestimation of vertical wind speed by the CSAT3 sonic anemometer due to transducer shadowing (Horst et al., 2015; Frank et al., 2016) were not included in our calculations.  39  2.2.5 Model parameters For simulations using the modified version of 3-PG, values of the physiological parameters were taken from the literature. Table 2.1 lists the model parameter values, and their sources, used in this study. Newly introduced parameters for the understory trees, and other site-specific parameters, were obtained from the literature, if available, and from our measurements, especially the site-specific parameters such as the SLA, stand age at the time of beetle attack, and the rate of canopy mortality due to beetle attack (e.g., parameters in Equation (7) mentioned above). The proportion of the average live tree biomass (foliage, root and stem) lost per dying tree was higher in our case than the default value used in most studies, because the beetle attacks larger diameter trees instead of suppressed trees, which usually die due to self-thinning (Sands, 2004). Model output sensitivity to input parameter values was studied by varying each parameter value by 25%.  2.2.6 Future projections We used our validated model to make projections up to 2025. For this purpose, we created five different climate scenarios (Table 2.2) appropriate for the Mackenzie region. The first two scenarios (S1 and S2) were, respectively, (i) an average of the climate measured from 2007 - 2015 for each of the subsequent years, and (ii) a reproduction of the monthly values for 2006 - 2015, meaning that the measurements from January 2006 to December 2015 were repeated over the following ten years (2016-2025) assuming that current interannual variability persists for the next decade. Furthermore, to take into account the change in climate that is likely to occur over the next decade in the Mackenzie region, we obtained rates of change in Ta and P 40  for this region from Price et al. (2011), and created two more climate scenarios (S3 and S4) corresponding to the Intergovernmental Panel on Climate Change special report on emissions scenarios (IPCC SRES) scenarios B1 and A2 (Nakicenovic et al., 2000), respectively. The fifth scenario (S5) was similar to the IPCC A2, but also included changes in VPD and Rs as projected by Price et al (2011).  Table 2.2  Details of climate scenarios for air temperature (Ta), precipitation (P) and vapour pressure deficit (VPD) during 2016 – 2025 used for predicting the effect of MPB attack on future carbon and water fluxes at MPB-06 over the next decade. Scenario Ta P VPD Rs S1 Average for 2007 - 2015 S2 Repeat monthly values for 2006 - 2015 S3 (IPCC B1) S1 + 0.02 K per year S1 + 0.1% per year S1 S1 S4 (IPCC A2) S1 + 0.05 K per year ” ” ” S5 (IPCC A2a) ” ” S1 + 0.25% per year Summer: S1 + 0.05% per year, Winter: S1 – 0.03 % per year  2.3 Results 2.3.1 Climate The interannual variability in monthly measured climate variables at MPB-06 for 2007-2015 is shown in Figure 2.3. While 2007 was the wettest year of this period, 2010 and 2012 were the driest. August 2009, January to March 2011, and June to September 2014 had exceptionally 41  low θ. PAR and VPD were low in June and July 2011 as well as in June 2012. Also, July 2014 showed high VPD together with low P and θ. From December 2010 to March 2011, values of θ were very low with values of 0.025 - 0.030 m3 m-3 in the surface 10-cm layer (Figure 2.3) and 0.045 -0.050 m3 m-3 in the 30-50-cm layer (not shown).  This was likely due to soil freezing as indicated by sub-zero Ts down to 20-cm depth and at times down to the 50-cm depth for those 4 months, which was unusual for the years observed. Ts at the 10-cm depth was between -0.45C and -0.70C in December 2010 to March 2011. In other years, monthly Ts only dropped to -0.15C in a few winter months.  42   Figure 2.3  Monthly 24-h average photosynthetically active radiation (PAR), cumulative precipitation (P), vapour pressure deficit (VPD) (solid grey line) and specific humidity deficit (D) (dashed black line), air temperature (Ta) (dashed line), soil temperature at the 10-cm depth 43  (Ts) (solid line) and soil water content (θ) for the 0-10-cm layer for 2007 - 2015 at MPB-06. Values of VPD and D were calculated from daytime average values.   2.3.2 Model results 2.3.2.1 Changing stem density and estimated LAI The modified version of 3-PG was run continuously for 2005 – 2015 with the initial conditions for the pre-attack period in 2005. Tree mortality due to the beetle attack was introduced in May 2006 at a stand age of 79 years. Monthly calculated and annual observed live overstory stem density at MPB-06 for 2005–2015 are shown in Figure 2.2. The beetle-induced mortality function (Equation (7)) included in the modified version of 3-PG resulted in a good fit between calculated and measured live stem density following the attack. There was a sharp decline in stem density in the first two years followed by a very slow decline over the following eight years. The canopy opened up due to overstory mortality as shown by hemispherical photography and the live overstory LAI decreased as shown in Figure 2.4. Averaged over the 24 photograph locations around the tower, canopy openness increased from 30.8% in 2007 to 43.3% in 2015. Measurements below the canopy at the 1-m height showed that the fraction of incoming Rs reaching the understory increased by about 10% between 2007 and 2009 thereby confirming an increase in canopy openness. On the other hand, PAR measurements at the 3-m height at the tower location showed transmitted light to be surprisingly constant during the first nine years after attack with an average of 29% of incoming PAR reaching the understory. In the model, total overstory LAI did not change much in the first two years following attack as the dead needles remained on the trees. When Ld decreased due to dropping of dead needles, the total 44  overstory LAI decreased as the increase in live overstory LAI due to newly forming needles was slower (Figure 2.4). In 2013, modelled total overstory LAI stabilized and then started to increase slightly, in agreement with measured values. The modelled understory LAI increased steadily over the 11 years, only slowing down in the last two years, and the total live LAI decreased following the attack, but recovered to pre-attack levels by 2013.     Figure 2.4  Monthly modelled and measured LAI at MPB-06 for 2005 - 2015. The modelled contributions of the two parts of the stand, the live understory and the live overstory LAI, are 45  shown in green and red, respectively, and the total live LAI is the blue line. The dashed black line is the modelled dead LAI and the dashed light blue line shows the total overstory LAI. The measured total overstory LAI values obtained from hemispherical photography are the black circles.  2.3.2.2 Comparison of modelled and measured E The comparison of the modelled and measured monthly E time series for 2005 – 2015 is shown in Figure 2.5. The model explained 88% of the variance in monthly E with an RMSE of 6.11 mm month-1 and a slope of the regression of 0.86, which was significantly different from 1.00 (p < 0.001) (Figure 2.6). All p-values reported in this paper result from t-tests performed at the 5% significance level. There was also generally reasonable agreement between measured and modelled annual E between 2007 and 2015 with mean values for the nine years of 250 and 245 mm y-1, respectively (R2 = 0.65, p < 0.01) (Figure 2.7 and Table 2.3). The model indicated a significant decrease in E (p < 0.05) of about 62% from 2005 to 2007 following the attack (Table 2.3). Model sensitivity tests using the different climate scenarios for 2005 – 2006 described in Section 2.2.3 resulted in modelled E for 2005 and 2006 being 83 – 95% and 96 – 102%, respectively, of the original modelled value when climate data was obtained from the nearby airport stations. The different climate scenarios resulted in slightly different values for all the years 2005 – 2015. The decrease in E from 2005 to 2007 for the different scenarios ranged from 53 to 60%. Thereafter, modelled annual E increased up to 2013. For 2007 to 2011, monthly and annual measured E, on the other hand, changed little from year to year. Then in 2012, values of measured E began to increase. In the first few years, the understory contribution was small, but it 46  increased more than two-fold over the ten years (since 2005). The model indicated that values of E similar to those in 2005 were reached again in 2015. Even though a large proportion of the main canopy trees died due to the beetle attack, measured E did not decline substantially during the measurement period.    Figure 2.5  Modelled and measured monthly E at MPB-06 for 2005 - 2015. The contributions of the two parts of the stand, the understory and the overstory, are shown in green and red, respectively, and the total is the blue line. The measured values are the black circles.  47   Figure 2.6  Measured monthly values versus modelled monthly values of E at MPB-06 for 2007-2015 including the 1:1 line (solid black line) and regression line (dashed red line). Also shown is the linear regression equation, coefficient of determination (R2), root mean square error (RMSE) and sample size (n).  48   Figure 2.7  Modelled and measured annual E and GPP at MPB-06 for 2005 - 2015. Bars show the 95% confidence interval for annual E and GPP including the uncertainty estimates described in Section 2.2.4.  Table 2.3  Measured annual precipitation (P), measured and modelled annual evapotranspiration (E), drainage below the 60-cm depth (Dr) and water use efficiency (WUE) for 2005 - 2015 at MPB-06 and the percentage of half-hourly values of measured E required to be gap-filled per year.  49  Year P [mm] E [mm] Dr [mm] WUE [g of C (kg of H2O)-1] % gap-filled E Measured Modelled Measured Modelled EC-estimated Modelled 2005 500  385  70  2.10  2006 563  262  223  2.23  2007 854 226 147 546 673 1.94 2.64 50 2008 662 237 173 415 530 2.11 2.54 41 2009 598 233 192 375 412 2.29 2.16 53 2010 507 230 220 307 308 2.48 2.35 43 2011 813 219 236 553 535 2.59 2.42 71 2012 623 265 270 373 318 2.16 2.12 53 2013 668 283 316 373 339 2.31 1.98 63 2014 594 233 288 378 461 2.39 2.00 64 2015 731 328 366 391 343 2.42 1.90 50 2007 – 2015 average 672 250 245 412 435 2.30 2.23 54  2.3.2.3 Monthly and annual GPP and WUE As a large proportion of the main canopy trees died due to the beetle attack in the first two years, modelled annual GPP significantly decreased (p < 0.05) by about 52% from 2005 to 2007, followed by a small increase in 2008, which was very similar to the increase in EC-estimated GPP (Figure 2.7). Model sensitivity tests using the different climate scenarios for 2005–2006 resulted in modelled GPP for 2005 and 2006 being 87–102% and 103–113%, respectively, of the original modelled value when climate data was obtained from the nearby airport stations. Thus, the decrease in GPP from 2005 to 2007 for the different scenarios ranged from 47 to 54%. Similar to modelled E, modelled GPP of the understory increased significantly (p < 0.01), by more than two-fold, in the nine years following attack (Figure 2.5 and Figure 2.8). 50  Monthly EC-estimated and modelled values of GPP, for 2007–2015, are compared in Figure 2.9. The model explained 90% of the variance in GPP with an RMSE of 15.2 g m-2 month-1 of C (slope significantly different from 1 with p < 0.001). Time series of modelled and EC-estimated annual GPP for 2005 to 2015 are shown in Figure 2.7. The model predicted the trend in annual GPP well with the recovery starting in 2008. Monthly WUE, calculated as GPP/E, for the overstory and the understory for the period 2005-2015 (not shown) showed some seasonal and interannual variability, but no clear trend. The linear regression between EC-estimated and modelled monthly values of WUE (not shown) indicated that the model explained 65% of the variance in WUE (highly significant with p < 0.001). Annually, modelled WUE decreased by about 10% over the ten years (Table 2.3). In 2007 and 2008, modelled WUE was higher than in the other years as the attack affected GPP less than E. EC-estimated WUE was lowest in 2007 and showed an increase of about 24% from 2007 to 2015, although not statistically significant (p ≈ 0.12). However, annual modelled WUE showed a significant decrease (p < 0.01) of 28% from 2007 to 2015.   51   Figure 2.8  Modelled and EC-estimated monthly GPP at MPB-06 for 2005 - 2015. The contributions of the two parts of the stand, the understory and the overstory, are shown in green and red, respectively, and the total is the blue line. The EC-estimated values are the black circles.  52   Figure 2.9  EC-estimated monthly values versus modelled monthly values of GPP at MPB-06 for 2007-2015 including the 1:1 line (solid black line) and regression line (dashed red line). Also shown is the linear regression equation, coefficient of determination (R2), root mean square error (RMSE) and sample size (n).  2.3.2.4 Stand water balance The annual course of the monthly averages of the water balance components modelled for 2007 to 2015 is shown in Figure 2.10. Also shown is rain plus snowmelt (i.e., liquid water, Lw), which peaked in April and May, when most of the snowmelt occurred. The seasonal distribution of Dr was bimodal with maxima in spring and fall due to snowmelt and the beginning of the rainy season in October, respectively. On average, Ep exceeded E from March 53  to October indicating that E was limited by θ during those months. The modelled mean growing season (1 June-30 September) E and Dr values for the whole measurement period (2007–2015), 167 mm and 46 mm, respectively, agreed well with the measured values (169 mm E and 49 mm Dr), especially considering that the uncertainties in the three variables (P, E and S) all contribute to the uncertainty in Dr. S was small, but showed that θ was recharged during snowmelt in the spring and at the onset of the rainy season in the fall.   Figure 2.10  Mean monthly water balance components and annual totals for MPB-06 for 2007-2015: measured precipitation (P), modelled evapotranspiration (E), modelled potential 54  evapotranspiration (Ep), modelled drainage (Dr), modelled change in soil water storage (S), and liquid water (Lw), measured as rain plus modelled snowmelt and equals P on a water year basis.  The modelled and measured annual values of the water balance components, E and Dr, for the nine years are compared in Table 2.3. Average annual measured (250 mm) and modelled (245 mm) E agreed well as did measured (412 mm) and modelled (435 mm) Dr. In the first few years, E was underestimated, but it increased along with GPP (Figure 2.7) corresponding to an increase in biomass. As expected, when E was underestimated, Dr was overestimated. The model explained 71% of the variance in Dr (p < 0.05) with the regression equation Dr meas = 0.56 Dr mod + 168.59 and an RMSE of 48 mm year-1. On an annual basis, modelled E and Dr on average accounted for 37% and 65% of annual P, respectively, in close agreement with the respective measured values of 38% and 61%. With the decrease in annual E between 2005 and 2007, modelled Dr increased from 70 mm in 2005 to 673 mm in 2007. The ratio of annual modelled Dr to measured P increased from 14% to 79% between 2005 and 2007, the latter being higher than the measured value of 64% in 2007. In the following years, up to 2015, this ratio decreased for both measured and modelled Dr, and the value of measured Dr in 2015 was 28% lower than in 2007. Measured annual values of S were small (not shown) and modelled values were close to zero (Figure 2.10).   2.3.2.5 Sensitivity to model parameter values The sensitivity of annual E to model parameters for 2007 (as an example) is shown in Figure 2.11. Modelled E was not very sensitive to model parameters, except for stomatal 55  sensitivity to VPD (sD), with E changing by less than 10% for most of the parameters when their values were varied by ±25%. The responses were symmetrical on both sides of the standard parameter value, i.e., linear, except for maximum stand age (ymax) and gcmax. E was negatively correlated with sD and positively correlated with gcmax as expected. Variation in ga and maximum rainfall interception had the least effect on modelled E, varying only by up to 1% and 3% with a ±25% variation in ga and maximum rainfall interception, respectively. GPP was most sensitive to α and sD showing 16 - 20% variation in GPP with a ±25% variation in α (results not shown). GPP was also negatively correlated with sD varying by ∓16% with a ±25% variation in sD. In contrast to E, GPP was negatively correlated with gcmax with a relatively small change of up to ∓2% in GPP resulting from a ±25% change in gcmax. E and GPP both showed a nonlinear sensitivity to ymax, which impacts gc and α through the age modifier.   56   Figure 2.11  Sensitivity of annual E for 2007 to the following model parameters: PAR extinction coefficient (k = 0.6), soil fertility (FR = 0.4), stomatal sensitivity to VPD (sD = 0.055 mbar-1), quantum efficiency (α = 0.047 mol of C (mol)-1 of PAR photons), maximum stand age (ymax = 150 years), boundary layer conductance (ga = 0.2 m s-1), maximum canopy conductance (gcmax = 0.012 m s-1), and maximum interception (MaxIntcptn = 0.15). Absolute values of E are on the left-hand y-axis and relative values of E are on the right-hand y-axis.  57  2.3.2.6 Model projections for the next decade We also made projections up to 2025 using the five different climate scenarios described in Section 2.2.6 assuming no further insect attack after 2015. The projected annual values of E, GPP and WUE for 2016 – 2025 are shown in Figure 2.12. As expected, the second scenario (S2) produced the most variable results because it included the current interannual variability in climate. Over the limited period of a decade, this showed a larger impact than the projected climate change incorporated in the scenarios S3 - S5. The regression lines for S2 show E increasing, GPP staying relatively constant and WUE decreasing from 2016 to 2025, but these trends are not significant. Predictions using the average climate for 2007 - 2015 (S1) showed E slowly increasing up to 2022 and then slightly decreasing, GPP staying relatively constant, and WUE decreasing slightly over the decade. S3 behaved similarly to S1 but resulted in lower E due to the slight increase in T and P in this scenario. With S4, E, GPP and WUE were very similar to S3. The last scenario, S5, which included changes in VPD and Rs as well as changes in T and P, resulted in E being very similar to that in S4. GPP and WUE for S5 were a little lower than for S4. Results show that the effect of the different climate scenarios on E, GPP and WUE would be very small over the next decade.   58   59  Figure 2.12  Annual modelled E (a), GPP (b) and WUE (c) at MPB-06 for 2016-2025 for the five climate change scenarios (S1–S5) shown in Table 2.2. For S2 the black dots are the annual values and the dashed black line is the regression line. Also shown are the linear regression equations for S2.  Table 2.4 shows the annual modelled E, GPP and WUE for 2005, 2006, 2015 and 2025. The estimated 2025 values for the five climate change scenarios except S2 predict reductions of 9-11% in E and of 15-21% in GPP compared to pre-attack values. WUE is expected to be reduced by 7-11% in 2025 compared to 2005. S2 predictions show reductions in E, GPP and WUE of 1%, 16% and 15%, respectively. Based on regression equations shown in Figure 2.12, reductions in E, GPP and WUE in S2 of 8%, 20% and 13% are similar to the reductions for the other scenarios. Modelled annual E and GPP in 2015 were higher than those for 2025, but modelled E compared well to measured E in 2015 (of 328 mm). The projections showed E increasing slightly, but not significantly at the 5% significance level.  Table 2.4  Annual modelled evapotranspiration (E), gross primary productivity (GPP) and water use efficiency (WUE) at MPB-06 for 2005, 2006, 2015 and 2025. Projections were done for the five climate change scenarios (S1–S5) shown in Table 2.2. The results for S2 represent the predicted trend using the linear regression equation. Year E [mm] GPP [g m-2 of C] WUE [g of C (kg of H2O)-1] 2005 385 809 2.10 2006 262 583 2.23 2015 366 696 1.90 60  Year E [mm] GPP [g m-2 of C] WUE [g of C (kg of H2O)-1] 2025 S1 349 685 1.96 2025 S2 354 647 1.83 2025 S3 342 654 1.91 2025 S4 343 655 1.91 2025 S5 344 641 1.87  2.4 Discussion EC measurements at MPB-06 were not made prior to the MPB attack in the summer of 2006 so that the initial reduction in E and GPP due to the attack could not be measured. The model, on the other hand, indicated that the attack and mortality of 70% of the canopy trees in the first year resulted in a reduction of annual E and GPP by 62% and 52%, respectively, in 2007 compared to 2005 (Figure 2.7). Brown et al. (2014) observed little variation in measured E in the first four years following attack and suggested that this was due to an increase in transpiration by the remaining live trees, understory, shrubs, herbs and mosses. The modified 3-PG model, as well as measurements by Emmel et al. (2013) and Bowler et al. (2012), confirmed the importance of the understory in keeping E relatively stable for nine years. Over the eleven years between 2005 and 2015, when mortality occurred in the overstory and more Rs penetrated through the canopy due to needle fall, understory E increased significantly (p < 0.001), more than two-fold, and almost compensated for the decrease in overstory E. The model indicated that the peak monthly E in July 2007 was about 54% lower than in July 2005 and that summertime E decreased by 57% from 2005 to 2007. This is significantly higher than the reduction of 19% in summertime E during the first few years following MPB attack that Maness et al. (2012) found using MODIS data over 170,000 km2 of affected area in BC. This difference is likely due to the 61  variability in tree mortality within the affected areas and the high mortality rate in the first year at MPB-06. Vanderhoof and Williams (2015), studying MODIS-derived E in MPB-affected lodgepole pine dominated stands in the Colorado Rocky Mountains, found that summertime E was only slightly lower (6.3%) for stands 4 - 13 years after attack compared to non-attacked stands, but was reduced by 18.7% for stands at 14 - 20 years after attack with an average mortality rate of 50%, and recovered fully in 20 - 30 years after attack. Our model calculations showed reduction in E and GPP in the first few years following attack, but both measured and modelled values showed the recovery in GPP starting in 2008. In 2012, measured E began to increase, which appeared to be related to the relatively high values of PAR and Ta in 2012 (Figure 2.3) and increased stand biomass, especially in the understory as indicated by its increased GPP (Figure 2.8) resulting in increased transpiration. Furthermore, as overstory mortality continued at a low rate and more Rs penetrated through the overstory due to needle fall, understory growth increased and its contribution to E increased significantly (p < 0.001). In 2014, measured and modelled E decreased when conditions were dry, but the model overestimated E in early summer (Figure 2.5 and Figure 2.7), possibly due, to some extent, to large gaps in our flux measurements in the spring of 2014. Our model calculations and measurements showed that annual E at MPB-06 in 2015, i.e., only 9 years after attack, reached 95% and 85%, respectively, of the estimated pre-attack value in 2005. Our sensitivity analysis showed that modelled and measured values of annual E in 2015 ranged from 100 to 114% and 90 to 103%, respectively, of the modelled value in 2005 obtained for the different 2005-2006 climate scenarios (see Section 2.2.3). Correlation between measured and modelled stand-level monthly E as well as EC-estimated and modelled GPP was high (highly significant with p < 0.001).  62  It was more difficult to model WUE well, since it is a ratio of GPP and E and errors in their estimates propagate, yet the correlation was reasonable and highly significant with p < 0.001. WUE, interestingly, did not show a significant trend in monthly and annual modelled values from 2005 to 2015, but annual modelled WUE decreased significantly (p < 0.01) between 2007 and 2015. EC-estimated WUE, on the other hand, suggested a weak increasing trend during that time (2007-2015). While it is possible that major disturbance can cause change in species composition, EC-estimated WUE not showing a significant trend, likely due to the dominance of lodgepole pine in the understory as well as the canopy, points to the contrary. A change detected in the stand following attack was an increase in canopy openness from 30.8% in 2007 to 43.3% in 2015 (Foord and Spittlehouse, 2014) resulting in an increased fraction of Rs reaching the understory. The degree of mortality and recovery varied spatially within the stand though, as shown by the fact that PAR reaching the 3-m height at the tower location stayed relatively constant between 2007 and 2015 with an average of 29% of incoming PAR reaching the understory, in contrast to the increase in Rs observed at other locations around the tower. In a nearby stand with about 50% of the stand density (trees of height > 10 m) of MPB-06, Emmel et al. (2013, 2014) found that seven years after attack, an average of 60% of the daily total incoming Rs and PAR reached the 1-m height. This agrees with our observed PAR values reaching the understory taking into account the difference in stand density.  We found that including a Ta limitation for gc in the modified version of 3-PG, which reduced E and prevented overestimation of E (and also GPP to some extent) in the spring (April to June), simulated the observed response to changes in Ta very well. In the winter months, on the other hand, E was underestimated by the model. As this occurred mainly in years with low Ta, we speculate that temperature was over-limiting modelled E in winter months. In contrast to 63  E, GPP was negatively correlated with gcmax, which was likely the result of secondary effects of increased E, when gcmax is higher, which leads to a reduced a and possibly water stress. Arora et al. (2016), using the process-based Canadian Terrestrial Ecosystem Model (CTEM, version 1.1) coupled to the Canadian Land Surface Scheme (CLASS, version 3.5), showed that reduced C uptake in BC’s forests of 328 Tg of C for 1999-2020 due to the MPB outbreak was overcompensated by the increase in productivity due to the changing climate and increasing CO2 concentrations, which was estimated to cause a sink of around 900-1060 Tg of C over the same period depending on the climate change scenario. Our model calculations showed that this beetle-attacked stand, which was not harvested, seems to have recovered faster than expected by earlier modelling studies (e.g. Kurz et al., 2008; Edburg et al., 2011), which indicated stand recovery times of 10 – 30 years to reach C neutrality depending on factors such as severity of attack and delay in snagfall with the latter accelerating decomposition. Brown et al. (2012) found that MPB-06 became C neutral within three years following attack, which is consistent with the increase in EC-estimated GPP from 2007 to 2010 (Figure 2.7). The recovery was fast, but not very strong, and GPP plateaued for several years thereafter. Model projections for 2016 - 2025 for all the climate scenarios suggest that GPP and E in this stand have reached their limits and will likely not increase significantly over the next decade, also indicating a low productivity stand. The relatively constant GPP relates to total LAI not changing much between 2020 and 2025, as understory LAI was predicted to peak in 2016 and slightly decrease thereafter with overstory LAI estimated to increase slightly. Our results together with previous studies (e.g. Brown et al., 2010; Edburg et al., 2011) suggest that not clearcut harvesting MPB-attacked stands having sufficient remaining living trees or understory would be an appropriate management practice if the aim is to maintain a positive C 64  balance and minimize the impacts on landscape hydrology. This applies especially to stands on coarse soils with limited production potential and water holding capacity, which offer other ecosystem services, such as winter habitat for Woodland Caribou as in the case of MPB-06 (Seip and Jones, 2009). Our measurements and modelled values indicated a significant decrease (p < 0.05) in E and GPP immediately following the attack, but measured E and EC-estimated GPP reached 85% and 98%, respectively, of their estimated pre-attack values within ten years. This occurred as the canopy opened up due to overstory mortality and on average more light penetrated through the canopy, so that understory E and GPP increased significantly (p < 0.01). As anticipated, Dr increased in the first year following the attack due to the decrease in E and higher P in 2007, and some overestimation of Dr in 2007, but with the recovery of E, Dr decreased again up to 2015 and is projected to remain relatively constant for the following decade. E and GPP in this not-harvested MPB-attacked stand recovered relatively quickly compared to clearcuts in the area, which have been found to remain annual C sources for up to 10 years (Brown et al., 2010). Also, Edburg et al. (2011) in their modelling study found that clearcut-harvested stands remained C sources for a prolonged time, up to 25 years. Mathys et al. (2013), studying a partial-harvested mixed conifer stand, where the dead lodgepole pine trees had been removed following MPB attack, suggest that this stand recovered more slowly than not-harvested stands. This and the fact that E and GPP are not expected to decrease within the next decade, suggest that from a hydrologic and C sequestration perspective the management decision not to salvage harvest immediately following attack can be advantageous depending on the stand characteristics (e.g., proportion of tree species that are not lodgepole pine, stand density, understory density, and soil texture).  65  2.5 Conclusions 1. The modified version of 3-PG performed well in modelling monthly and annual stand-level E during the nine years after MPB attack. Modelled values of annual E decreased by about 62% in the first year after attack compared to pre-attack estimates. Measured annual E remained remarkably stable for the following four years, and then increased somewhat in the last four years.  2. The water balance showed a bimodal Dr distribution during the year with maxima occurring in the spring and fall seasons. Modelled Dr showed good agreement with Dr estimated from measured P, E and S. On average, modelled annual E and Dr accounted for 37% and 65% of P, respectively, in close agreement with measurements.  3. There was good agreement between modelled and EC-estimated monthly and annual GPP. Modelled annual GPP decreased by about 52% in the first year after attack, but then both modelled and EC-estimated GPP showed a slow steady recovery over the nine years. Modelled monthly WUE showed reasonable agreement with EC-estimated monthly WUE, and annual modelled WUE decreased slightly over the ten years. EC-estimated WUE, on the other hand, was lowest in 2007 and increased up to 2015, because EC-estimated GPP increased more than measured E. 4. Parameter sensitivity analysis indicated that modelled E and GPP had low sensitivity to boundary layer conductance and maximum rainfall interception, and high sensitivity to stomatal response to VPD and quantum efficiency of photosynthesis. Similar sensitivities of E and GPP to these parameters showed the overall robustness of the model.  66  5. Over the next ten years, E is expected to increase slightly while GPP is predicted to plateau. WUE is expected to decline slightly.  67  Chapter 3: Simulation of net ecosystem productivity of a lodgepole pine forest after mountain pine beetle attack using a modified version of 3-PG  3.1 Introduction The MPB (Dendroctonus ponderosae) outbreak, which started in the late 1990s, has killed ~731 million m3 of commercial timber over 18.3 million ha, equivalent to 54% of the mature merchantable lodgepole pine (Pinus contorta var. latifolia) volume in BC. The infestation is expected to subside in BC by 2020 with an estimated 55% of the mature merchantable pine volume killed (BC Ministry of Forests, Lands and Natural Resource Operations, 2016). Tree mortality due to the beetle outbreak impacts the C balance of forests (Brown et al., 2010), as C uptake is reduced when the rate of photosynthesis of the attacked stands decreases. Loss of C may increase or decrease depending on how Rh, especially the decomposition of coarse woody debris (CWD) (i.e., decaying roots, stems, branches), and Ra respond to tree mortality. Modelling results suggest that stands that were once C sinks become C sources, i.e., the decline in GPP exceeds the decrease in R. Kurz et al. (2008) estimated the impact of the MPB attack in BC between 2000 and 2020 could result in a loss of 270 Tg of C, equivalent to 36 g m-2 year-1 of C on average over 374,000 km2 of forest, due to a reduction in NPP and an increase in Rh compared to the control scenario without beetle attack, plus an additional 50 Tg of C lost due to salvage harvesting the beetle-killed and living trees. Alternatively, Williams et al. (2016) found that forest resilience, and increased CO2 concentrations in the atmosphere enhanced growth increasing C stocks by about 190 Tg year-1 of C in the US and by about 10-30 Tg year-1 of C in 68  Canada at the same time as disturbances (harvesting, fire, windthrow, insect attacks and drought) caused a decrease in live biomass of about 200 Tg year-1 of C in the US. Using the coupled CLASS-CTEM for 1999-2020, Arora et al. (2016) concluded that the loss of 328 Tg of C due to the MPB attack in BC may be surpassed by the enhanced C accumulation due to changing climate and increasing atmospheric CO2 concentration, resulting in a sink of 900-1060 Tg of C over the 20 years, depending on the future climate change scenario.  Using the CLM with prognostic C and nitrogen cycling (Version 4), Edburg et al. (2011) found that stand-level NPP is reduced and thus the C balance is altered in MPB-attacked stands for several years after the disturbance. The modified CLM, which included snag and dead foliage C pools, showed that the outbreak severity, delay in snagfall, snagfall rate and management decisions all impacted C fluxes for decades, and C stocks up to 100 years after attack. Fluxes of C were impacted by outbreak duration in the short term, but the number of trees killed affected C fluxes and stocks in the long term (Edburg et al., 2011). Moore et al. (2013), using a combination of satellite observations, CO2 observations (EC measurements as well as atmospheric CO2 concentration measurements to obtain the temporal dynamics of R), soil chamber flux and soil organic C content measurements of an uninfested forest and an attacked forested valley about 30 km apart, showed that MPB outbreaks have little effect on NEP on a decadal scale, as both GPP and R decline after the trees die. Previous modelling studies might have overestimated the immediate impact of insect outbreaks on C release, as Moore et al. (2013) found an increase in R for only 1-2 years after the attack due to decomposition of needle and root litter, followed by lower C input into the soil causing a reduction in R. They suggested the need for further studies to improve model estimates. EC measurements in a bark-beetle infested lodgepole pine stand in southeast Wyoming showed that 69  CO2 uptake and WUE did not change with tree mortality increasing from 30% to 78% over three years, suggesting that the remaining trees and sub-canopy vegetation increased their CO2 uptake and water use following attack (Reed et al., 2014). In MPB-attacked pine forests of the BC interior, both component flux (Bowler et al., 2012) and EC measurements (Emmel et al., 2014) showed that the residual vegetation following MPB attack, including broadleaf understory and younger live coniferous trees (secondary structure), were important C sinks in the years following attack, and by contrast, the overstory consisting mainly of dead lodgepole pine trees was a small C source. Clearcut- and partial-salvage harvesting, or no harvesting, have been the most common approach to managing beetle-attacked stands in BC (Coates and Hall, 2005). C and water balances of differently managed MPB-affected lodgepole pine stands in northern interior BC have been studied using EC flux measurements at three tower sites with growing-season EC measurements in nearby clearcuts (Brown et al., 2010; Mathys et al., 2013; Emmel et al., 2013 and 2014). This paper focuses on the long-term effects of MPB attack on GPP, R and NEP of an interior BC lodgepole pine stand and its recovery following attack using flux tower EC data, and the modification of the 3-PG model (Meyer et al., 2017) to include the calculation of Rh for modelling such effects. The original 3-PG (Landsberg and Waring, 1997) has been shown to predict growth and development of even-aged, intensively-managed, fertilized stands of Eucalyptus globulus well (Sands and Landsberg, 2002) and has been applied for different species including Pinus radiata (Coops et al., 1998), ponderosa pine (Pinus ponderosa) (Coops et al., 2005) and lodgepole pine (Pinus contorta Dougl.) (Coops and Waring, 2011). Sands and Landsberg (2002) modified the original 3-PG model and included an age-dependency of the SLA and the fraction of stem biomass (ws) as branches and bark, which enables calculation of 70  variables such as stem volume and MAI in addition to predicted stem biomass. They also adjusted the productivity at cool sites compared to warm sites by including a temperature-dependent growth modifier (Sands and Landsberg, 2002). We have modified 3-PG to represent the effects of beetle attack on the water balance and WUE of a lodgepole pine stand (Meyer et al., 2017). The objectives of this subsequent research are to (i) include a treefall function and Rh sub-model, including the soil and CWD components, in the modified version of 3-PG in order to calculate R and NEP, and (ii) compare modelled values of GPP, R and NEP, calculated using parameter values mostly taken from the literature with some estimated from in-situ measurements, with EC-measured fluxes in an 80-year-old lodgepole pine stand during the first decade following MPB attack.  3.2 Methods 3.2.1 The model Landsberg and Waring (1997) developed the process-based stand growth model, 3-PG, which is run on a monthly time step and determines biomass production, allocation of biomass between foliage, roots and stems, stem mortality, the soil water balance and stand characteristics (such as SLA, LAI, MAI and mean DBH) in five sub-models. Weather data (solar irradiance (Rs), Ta, VPD, P), site factors (e.g., species, site latitude), initial conditions (e.g., a in the 0-60-cm layer, stocking) and physiological parameter values characterizing the species to be modelled (e.g., α, gc, and ga) are required inputs (Sands, 2004). Further details on the original 3-PG model are described in Landsberg and Waring (1997) and Landsberg and Sands (2011). Briefly, a light-use efficiency approach is used to calculate GPP (g m-2 month-1 of C) and NPP is determined as a 71  constant fraction (47%) of GPP (Coops and Waring, 2011; Sands, 2004; Campioli et al., 2015), which results in Ra being 53% of GPP. The allocation of NPP to roots, foliage and stems in 3-PG depends on a, VPD, site fertility and the age of the stand (Sands, 2004). An age-dependent mortality function is implemented in the model and the -3/2 self-thinning law (Drew and Flewelling, 1977) is used to include density-dependent stem mortality. Stem volume, mean DBH and other stand characteristics of interest to forest managers are obtained from stand density and the biomass pools (Sands, 2004). 3-PG includes a simple monthly soil water balance, where P (mm), Dr (mm) and E (mm) determine S (mm). The Penman-Monteith equation (Landsberg and Gower, 1997; Monteith and Unsworth, 2013) is used to calculate actual E with a fixed value for ga of 0.2 m s-1 and gc determined as the maximum gc limited by LAI, a, VPD and atmospheric CO2 concentration (Landsberg and Sands, 2011). We have revised 3-PG for application to MPB-attacked lodgepole pine forests (Meyer et al., 2017). These and our further modifications made in this study to the basic structure of 3-PG are shown in Figure 3.1. Briefly, these included modelling a 2-layer canopy, with the growing understory and partly dying overstory following MPB attack as well as inclusion of snowmelt to account for the total amount of liquid water instead of only rainfall. LAI, GPP, NPP, E and Ra for the understory and overstory are calculated separately using the same set of equations as in the original 3-PG with their respective initial values and parameters. The stand-level values are obtained as the sum of the values for the two canopy components. Overstory death resulting in the reduction of live overstory LAI following beetle attack is expressed by an exponentially declining mortality rate γn in the insect-induced mortality function (as described in Section 2.2.2, Chapter 2), which was added to the age- and density-dependent mortality functions in the original 3-PG. Differences in Rs between the overstory and understory were taken into account 72  and the Rs reaching the understory (Ru) was controlled by the overstory extinction coefficient (k) and the total overstory LAI (L) including live and dead needles, giving Ru = Rs exp(-k L). With the increasing number of dying trees following attack, Ld was accumulated and then reduced beginning two years after the start of the attack due to needle fall, which was described by an exponential function (Meyer et al., 2017).    Figure 3.1  Basic structure of 3-PG including modifications made (in red) for its application to MPB-attacked lodgepole pine stands. The subscripts ‘o’ and ‘u’ stand for the overstory and understory, respectively. α is the canopy quantum use efficiency of photosynthesis, φ is the physiological modifier applied to α and the maximum canopy conductance (gcmax), ws is the 73  average tree stem biomass, γn is the tree mortality rate, Soil C and CWD C are the soil and coarse woody debris (CWD) C pools, and Rh soil and Rh cwd are the heterotrophic respiration components originating from the soil and CWD, respectively. In this version of 3-PG, soil C is held constant (see Section 3.2.1.1) with no input from fine root turnover or litter. The visual basic code for the modified version of 3-PG is available on request (gesa.meyer@ubc.ca).  Further modifications to 3-PG included a Rh sub-model to determine the effect of beetle attack on the C balance, thereby enabling the calculation of NEP, as Ra was already calculated in the original 3-PG. This sub-model has two components of Rh: (1) decomposition of soil organic matter (Rh soil), and (2) decomposition of CWD from the dead trees (Rh cwd), so that Rh = Rh soil + Rh cwd. Similar to the approach used by Arain et al. (2002), the two components of Rh (g m-2 month-1 of C) are calculated as             𝑅h soil =  𝑅10 soil𝐶soil𝑄10𝑇s−1010                   (9)                      𝑅h cwd =  𝑅10 cwd𝐶cwd𝑄10𝑇s−1010        (10) where R10 soil and R10 cwd (g kg-1 month-1 of C) are the base-respiration rates at 10C of the soil and CWD, respectively, Csoil (kg m-2 of C) and Ccwd (kg m-2 of C) are the soil (0-60-cm depth) and CWD C storage, and Q10 is the relative increase in R10 soil and R10 cwd for an increase in temperature of 10C. Temperature of the CWD was approximated by Ts at the 5-cm depth. Shallow Ts was a reasonable approximation of the temperature of CWD lying on the forest floor especially during the winter when Ta was 10-15C less than Ts and CWD temperature beneath the snow. R is obtained by summing Ra and Rh. NEP is then calculated as the difference between GPP and R, which is equivalent to the difference between NPP and Rh. We decided not to 74  include  since in moist-soil ecosystems, Ts and  are generally correlated (Davidson et al. 1998). Further development of the Rh model should include other factors such as the role of enzymes and reaction kinetics (Davidson et al. 2012) but is beyond the scope of this paper.  3.2.1.1 Model parameter values The physiological parameter values used to run the modified version of 3-PG were taken from the literature while site specific parameters such as the stand age, the beetle-caused canopy mortality rate, and the treefall rate were obtained from our measurements (see Meyer et al., 2017 for the complete list of the parameter values). Unlike in the case of self-thinning, where suppressed trees die (Sands, 2004), insect-induced mortality affects larger diameter trees and causes the proportion of the average live-tree biomass lost per killed tree to be higher than the 3-PG default value.  Values of R10 soil and R10 cwd were estimated from in-situ measurements of Rh (see Section 3.2.2). Csoil was considered constant (see Table 3.1). The usual assumption that such systems are in equilibrium with annual C input in root and needle litterfall balanced by decomposition is probably not realistic in heavily disturbed ecosystems, such as MPB-attacked stands. However, based on live needle C measurements made in 2006 at MPB-06 potential litterfall would have amounted to ~0.2 kg m-2 of C (~6% of Csoil) and would likely have caused a negligible increase in Rh in the following years. As standing dead trees were assumed to respire at a negligible rate (Hansen, 2014) and there was a delay between the trees dying and falling, C from the dead biomass was added to the CWD C pool, Ccwd, dependent on the treefall rate starting 1.5 years 75  after the beginning of the attack. Treefall was included in the model using a Gaussian fit of the form:        2 2 exp( ( ) )  exp( ( ) )x b x ey a dc f                    (11) with x being stand age (months), y being treefall rate (stems ha-1 month-1) and coefficients a through f (see Table 3.1) determined using a nonlinear least-squares fit to annual treefall counts (see Section 3.2.2).  This two-term Gaussian function (curve fitting tool in MATLAB®) represented the shape of the treefall observations ( Figure 3.2a) well with a coefficient of determination (R2) of 0.997, a RMSE of 0.422 stems ha-1 month-1 and a sum of squares due to error (SSE) of 0.713 (stems ha-1 month-1)2. The monthly average of the observed annual treefall rate is shown together with the statistical fit in Figure 3.2a. The parameter values used in the Rh sub-model and their sources are listed in Table 3.1.  Table 3.1  List of parameters used in the heterotrophic respiration (Rh) sub-model in the revised 3-PG, and their values for lodgepole pine (Pinus contorta var. latifolia). 3-PG parameter Name Units Lodgepole pine Source Heterotrophic respiration       Base-respiration rate at 10 C of the soil R10soil g kg-1 month-1 of C 7.5 This study Q10 for heterotrophic respiration Q10  - 2 Brown et al. (2012) Soil C content Csoil kg m-2 of C 3.5 Brown et al. (2010) Base-respiration rate at 10 C of the CWD R10cwd g kg-1 month-1 of C 0.6 This study Average C content per standing dead tree Cdead avg kg tree-1 of C 140 NFI1 Gaussian fit coefficients     Treefall function coefficient a  - 12.36  (10.27, 14.44)2 This study Treefall function coefficient b  - 87.08  (86.87, 87.3) This study Treefall function coefficient c  - 0.614  (0.4086, 0.8195) This study 76  3-PG parameter Name Units Lodgepole pine Source Treefall function coefficient d  - 3.958  (2.27, 5.646) This study Treefall function coefficient e  - 86.96  (86.17, 87.75) This study Treefall function coefficient f  - 3.57  (2.069, 5.071) This study 1 National Forest Inventory 2 95% confidence interval in brackets.   Figure 3.2  (a) Calculated monthly treefall rate (black curve) at MPB-06 using a Gaussian fit through the monthly average of the observed annual treefall rate (blue bars) and (b) cumulative monthly modelled coarse woody debris (CWD) C pool content for 2005 - 2016.  3.2.2 Measurements A decade (2007-2016) of EC flux and climate measurements collected at a site (MPB-06), described in detail in Brown et al. (2010), were used to compare with modelled values using the modified version of 3-PG. A 32-m-tall scaffold tower (~2.1 m long x ~1.5 m wide) was 77  erected in July 2006 in this 80-year-old 15-m-tall (in 2007) extensive MPB-attacked almost pure lodgepole pine stand approximately 35 km southeast of Mackenzie at Kennedy Siding (550642.8 N, 1225028.5 W) with a fetch greater than 1 km in all directions from the tower. The stand was not harvested following attack in 2006 and has an understory mainly of pine seedlings and saplings on a gravelly sandy loam soil. About 50% of the canopy trees were green-attacked in August 2006, and by August 2012, about 86% of the canopy trees had been attacked (Brown et al., 2012; D. Seip, BC Ministry of Environment, Prince George, personal communication, 2016) with 14% of the pre-attack canopy trees alive since then, as no further attack has been observed. Together with the observations of stage of attack, numbers of standing dead as well as fallen trees were recorded annually along 30 different 100 m long transects within a 3000-hectare wildlife habitat reserve surrounding the flux tower (D. Seip, BC Ministry of Environment, Prince George, personal communication, 2017). Trees appeared to break off at ground level, when falling, so that the roots did not come out of the ground. The average C content of a standing dead tree (Cdead avg), used to calculate Ccwd in 3-PG as the cumulative of the product of Cdead avg and the treefall rate, and Csoil were obtained from measurements of stand density of dead and alive trees, C contents of standing dead trees, CWD, litter as well as the mineral soil analysis conducted in three ground plots around the MPB-06 tower in September 2006 following Canada’s National Forest Inventory (NFI) ground-plot protocols (NFI, 2008). Hemispherical photography done once a year between 2007 and 2016 was one of three optical approaches used to measure canopy LAI (Brown et al., 2010). A camera with a fisheye converter lens was used to take hemispherical photographs at the 1.2-m height at 24 locations along 150-m long transects in four directions (NW, NE, SW and SE) from the tower (Foord and Spittlehouse, 2014).  78  A detailed description of the flux and climate measurements can be found in Brown et al. (2010). Briefly, an open-path IRGA (model LI-7500, LI-COR, Inc., Lincoln, Nebraska) and a three-dimensional ultrasonic anemometer (model CSAT3, CSI, Logan, Utah) were used to measure CO2, water vapour and sensible heat fluxes above the canopy at a height of 26 m. Half-hourly NEE was calculated as the sum of the measured CO2 fluxes and changes in air-column CO2 storage (below the EC sensors). CO2 storage changes in the whole air column below the EC measurement height (Hollinger et al., 1994) were determined as the product of the dry air density and the difference between half-hourly averaged CO2 mixing ratios measured at the 26-m height for the following and previous half-hours divided by the difference in time between the measurements, following Brown et al. (2010) and Morgenstern et al. (2004). NEP was obtained as –NEE. Daytime R was calculated by extrapolating the nighttime relationship between R (EC-measured NEE) and Ts at the 5-cm depth to daytime and GPP was obtained as the sum of daytime NEP and daytime R (Brown et al., 2010).  Climate measurements included above-canopy down-welling and up-welling shortwave and longwave radiation (net radiometer, model CNR1, Kipp and Zonen B.V., Delft, The Netherlands) at the 30-m height, and above-canopy down-welling and up-welling (30-m height), and below-canopy down-welling photosynthetic photon flux density (quantum sensor, model LI-190 AS, LI-COR Inc.) at the 3-m height. Ta was measured at the 25-m height. Ts was measured at depths of 5, 10, 20 and 50 cm, and  was measured at the 0-10 cm and 30-50 cm depths. A tipping bucket rain gauge (not heated) was used to measure rainfall at canopy height and wintertime P was estimated from snowfall data. In a nearby clearcut (~1.5 km east of the MPB-06 tower), snow depth was measured acoustically (model SR50A, CSI) and snow depth and manual measurements of snow water equivalent were used to calculate snowfall and its liquid 79  water equivalent. Wind speed, RH and G were also measured at the tower site. To run 3-PG for 2005 and 2006, climate measurements from 1 Jan 2005 to 15 July 2006 were obtained from Mackenzie Airport (~28 km northwest of the site), supplemented by data from Prince George Airport (~135 km south of the site). Climate data from these airports agreed well with measurements at MPB-06 with coefficients of determination, slopes of the regressions and RMSEs for Ta of 96%, 0.90 and 2.21 C for July 2006 to February 2007 at Mackenzie Airport, and 99%, 1.02 and 1.08 C for July 2006 to October 2009 at Prince George Airport, respectively. As Rs and VPD were not measured at these airports, they were calculated for 2005-2006 using linear regression relationships with Ta obtained from measured values from 2007-2015.  Brown et al. (2010, 2014) described the quality control and data analysis procedure in detail. Briefly, half-hourly EC flux data were excluded, if more than 30% of a trace in a half-hour were flagged by the quality control or if IRGA measurements were not within bounds set by minimum and maximum values for water vapour and CO2 mixing ratios. A u* threshold (u*th) of 0.30 m s-1 (Brown et al., 2010) was applied to nighttime EC data to ensure sufficient turbulent mixing. To exclude data with wintertime CO2 uptake (Burba et al., 2008), negative NEE fluxes during the wintertime (defined as times when Ts at the 5-cm depth < 1 oC and Ta < 5 oC) as well as winter daytime data when wind speed was < 4 m s-1 were removed (Brown et al., 2010). In August 2014, soil collars were installed in 16 locations in 8 plots, five of which were root-exclusion plots, and soil R measurements were conducted using a portable chamber system equipped with a closed-path IRGA (model LI-800, LI-COR Inc.) and an opaque chamber (chamber volume = 0.0014 m3). Measurements of R were also made on two live and two dead standing trees with DBH ranging from 11 to 15 cm, as well as on two similar sized logs, which would have likely fallen in the previous 5-10 years. Photographs of the chamber system and the 80  measurement locations are shown in Appendix B. The rate of change of CO2 concentration measured in the chamber headspace was used to calculate the flux. The R10 value for soil Rh was estimated from the measurements in the root-exclusion plots and R10 cwd was estimated from measurements on the logs. In-house software written in MATLAB® (Version 7.5.0, The Math Works Inc., Natick, MA, USA) was used to perform all data calculations and analysis for this paper. All p-values were obtained from t-tests performed at the 5% significance level.  3.2.3 Uncertainty analysis of EC and modelled fluxes Uncertainty in EC-measured annual NEP values was obtained as the random error, calculated as the square root of the sum of squares of the random error associated with the half-hourly NEP measurements and uncertainties in the gap-filling procedure, plus the systematic error due to the selection of u*th (Krishnan et al., 2006). A 20% random error was assigned to each half-hourly value of NEP to obtain the random error of annual NEP (Morgenstern et al., 2004). A uniformly distributed random number generator produced gaps in measured NEP ranging from one half-hour to 10 days. Then, a Monte Carlo simulation was used to obtain the uncertainty in gap-filling, which was performed using modelled values from R versus Ts (empirical logistic) and GPP versus PAR (rectangular hyperbola, i.e., the Michaelis-Menten light response equation) relationships (Krishnan et al., 2006). This was repeated 500 times to determine the 95% confidence interval. NEP calculated with a variation in u*th of up to 50% provided an estimate of the systematic error associated with the u*th selection (Morgenstern et al., 2004). Uncertainties in EC-measured annual GPP and R were estimated using the same procedure. In our calculations, uncertainties such as those caused by transducer shadowing in the 81  CSAT3 sonic anemometer possibly leading to underestimation of vertical wind speed (Horst et al., 2015; Frank et al., 2016) were not included (Meyer et al., 2017).  Uncertainties in modelled annual NEP, GPP and R were obtained by defining value ranges for the key model parameters (see Table 3.2), which modelled NEP, GPP and R were the most sensitive to, and then running a Monte Carlo simulation. The model was run 500 times with randomly chosen, by a uniformly distributed random number generator, parameter values within their defined ranges to determine the 95% confidence intervals for the annual fluxes.   Table 3.2 List of value ranges for the model parameters used in the uncertainty calculation of modelled fluxes. 3-PG parameter Name Units Value ranges Canopy quantum efficiency alpha molC molPAR-1 0.045 – 0.055 Understory quantum efficiency alphau molC molPAR-1 0.03 – 0.04 Ratio NPP/GPP Y - 0.30 – 0.72 Extinction coefficient for absorption of PAR by canopy k - 0.5 - 0.6 Extinction coefficient for absorption of PAR by understory ku - 0.5 - 0.6 Defines stomatal response to VPD CoeffCond/sD mbar-1 0.05 – 0.06 Soil fertility FR - 0.3 – 0.5 Base-respiration rate at 10 C R10soil g kg-1 month-1 of C 3 - 15 Q10 for heterotrophic respiration Q10  - 2 - 3 Soil C content Csoil kg m-2 of C 2.2 - 3.5 Base-respiration rate at 10 C of the CWD R10cwd g kg-1 month-1 of C 0.15 – 1.50  3.3 Results and Discussion 3.3.1 Time series of the climate variables The monthly measured climate variables at MPB-06 for the decade 2007-2016 (Figure 3.3) show some interannual variability with 2007 and 2016 being wet years and 2010 and 2012 82  being the driest of this period. θ was generally low in summer months, especially during 2009, 2014, 2015 and 2016, as well as in December 2010 to March 2011. From December 2010 to March 2011, θ was very low in the surface 10-cm layer (Figure 3.3) with values of 0.025 - 0.030 m3 m-3 and 0.045 -0.050 m3 m-3 in the 30-50-cm layer (Figure 3.3). June and July 2011, June 2012, August 2015, and April 2016 to June 2016 had low PAR and VPD, leading to the annual PAR value in 2016 being 92% of the annual average during the decade. VPD was high in July and August 2014, coinciding with low P and θ. Usually, monthly Ts at the 5-cm depth did not drop below -0.4C and at the 10-cm depth not below -0.15C at this site. Unlike the other years, the 2010-2011 winter showed sub-zero Ts down to the 20-cm depth and at times down to the 50-cm depth with Ts at the 5-cm depth falling to between -0.59C and -1.08C, and at the 10-cm depth to between -0.45C and -0.70C from December 2010 to March 2011.  83   Figure 3.3  Monthly 24-h average photosynthetically active radiation (PAR), cumulative precipitation (P), vapour pressure deficit (VPD), air temperature (Ta) (grey line), soil temperature at the 5-cm depth (Ts) (black line) and soil water content (θ) in the 0-10-cm layer (black line) and  84  the 30-50-cm layer (grey line) for 2007 - 2016 at MPB-06. Values of VPD were daytime average values.  3.3.2 Comparison of measured and modelled C balance components Using climate data and initial conditions in 2005 at MPB-06, and independently obtained model parameter values, the modified version of 3-PG was run continuously from 2005 to 2016. The stand was 79 years old, when the beetle attacked in May 2006 and canopy mortality was introduced in the model. Most of the model parameter values used in this study can be found in Meyer et al. (2017) and Table 2.1 with additional parameters required for the simulation of NEP listed in Table 3.1.   3.3.2.1 Treefall and coarse woody debris The monthly average of the annually observed treefall rate at MPB-06 is shown in Figure 3.2a together with the calculated monthly treefall rate in the model using a Gaussian fit to the observations. Treefall at this site started about 1.5 years after the beginning of the main attack in 2006 and both the observed and calculated treefall rates peaked in 2014, 8 years after the attack started, and were markedly lower in the years thereafter. There was no indication in the observed wind data (not shown) that storms caused the higher treefall rates in 2014. The area under the calculated curve is 98% of the observed total treefall. With the trees falling at the site, the CWD C pool content increased up to 2016 reaching a value of 57.5 Mg ha-1 of C (Figure 3.2b), which is more than one and a half times the amount of Csoil, as more C in CWD was added than was 85  respired. Previous studies in BC have shown negligible rates of treefall within the first 3-6 years following mortality due to beetle attack, increasing to 25-50% of the dead trees being down within 9 years and 80-90% of the trees having fallen 15-20 years post-mortality (Lewis and Hrinkevich, 2013; Lewis and Hartley, 2005). Lewis and Hrinkevich (2013) attributed this variability in treefall rates to differences in , wind speed and topography between sites. Mitchell and Preisler (1998) also reported similar treefall rates in Oregon, where they found treefall to begin 3 and 5 years following death in thinned and non-thinned stands, respectively, and 50% of the trees having fallen after 8-9 years. Our treefall rates were slightly higher shortly after attack with 7% and 44% of the dead trees fallen within the first 5 and 8 years after attack, respectively, and 55% of the dead trees down after a decade. Hansen et al. (2015), who simulated trajectories of total, live and dead, C stocks in the central US Rockies, showed a long-term decrease of 1-30% due to MPB outbreaks. In their study, C was redistributed from live to dead pools following outbreak, but aboveground C stocks were resilient, as tree growth recovered and snags decomposed slowly. Hansen (2014) also showed that snags can remain standing for ~10 years or longer and even fallen boles decompose over many decades. Our portable chamber measurements, made under warm and very dry conditions in August 2014, as well as bole respiration measurements in another MPB-attacked stand ~72 km south of MPB-06 made by Emmel et al. (2014) showed that standing dead trees respired very little, only at about a quarter of the rate for live trees, and even logs lying on the ground only respired at about a half of the rate for live trees, agreeing with the slow decomposition of boles found by Hansen (2014).  Our results showing an increase in CWD C pool content (Fig. 2b) as well as fast recovery in GPP (see Section 3.2.4) agrees with the 86  resilience of aboveground C stocks with a redistribution of C from live to dead pools that Hansen et al. (2015) reported.  3.3.2.2 Heterotrophic respiration Modelled monthly Rh at MPB-06 for 2005-2016 is shown in Figure 3.4 together with its two components Rh soil and Rh cwd. The decomposition of soil organic matter was the dominant component of Rh, but Rh cwd increased significantly (p < 0.001), as CWD accumulated from the falling dead trees, making up 11% of Rh in 2016. As both Rh soil and Rh cwd increased significantly (p < 0.01) from 2005 to 2016, total Rh increased significantly by 25%. The increase in Rh soil was due to a small but significant increase in Ts, particularly during the growing season. Edburg et al. (2011) found an increase in Rh for up to 5 years following beetle attack because of the C input from fine and coarse roots and foliage, which decomposed quickly, followed by some decrease in the following five years. However, we did not find the decline in Rh that Edburg et al. (2011) found in the second half of the decade following attack because the CWD pool kept increasing. Kurz et al. (2008) also found an increase in Rh three to seven years after the beginning of beetle attack due to the decomposition of the large addition of biomass to the dead organic matter pools. They estimated that this beetle outbreak resulted in a 6% increase in Rh. This increase is much lower than our estimated increase of 25% from 2005 to 2016, and may be attributable to the high mortality rate at MPB-06. The Kurz et al. (2008) study is based on a large region in BC encompassing stands with a wide range of mortality rates.   87   Figure 3.4  Modelled monthly Rh at MPB-06 for 2005-2016. The contributions of the two components of Rh, coming from the soil (Rh soil) and coarse woody debris (Rh cwd), are shown in red and green, respectively, and total Rh is the blue line.  3.3.2.3 Ecosystem respiration Modelled monthly R values at MPB-06 for 2007-2016 (Figure 3.5) agreed well (R2 of 82%) with monthly EC-estimated R values (Figure 3.6a). The RMSE was 17.21 g m-2 month-1 of C and the slope of the regression was 0.99, which was not significantly different (p = 0.25) from 1.00 (Table 3.3). The two components, Ra and Rh, are also shown in Figure 3.5. On average, 88  during 2007 to 2016, Ra represented 55% of R with the rest being Rh (see Table 3.4). This ratio varied throughout the years with Ra being 48 to 66% of R on an annual basis. Due to a significant decrease (p < 0.05) in annual Ra and little change in Rh during 2005-2007, modelled R decreased by 35% (p < 0.05) (Figure 3.7). In the following decade, both Ra and Rh (p < 0.0005), as well as modelled (p < 0.00001) and EC-estimated R (p < 0.01) showed a significant recovering trend with increases in Ra, Rh, modelled and EC-estimated R of 102%, 26%, 62% and 50%, respectively, by 2016. Annual modelled and EC-estimated R showed good agreement with an R2 of 0.73 (Figure 3.6b and Table 3.3). In the first three years following attack and in 2015, annual modelled R was lower than EC-estimated R, but their uncertainty ranges overlapped. In the other years, modelled R was within the measurement uncertainty range. Looking at the changes from 2005 to 2016, only Rh increased significantly (p < 0.0005) (Figure 3.5). By 2015, both EC-estimated and modelled R exceeded the estimated pre-attack value in 2005 (Figure 3.7) with 2014 being an exception because of decreases in both GPP and R due to drought during the growing season.  89   Figure 3.5  Modelled (Rmod) and EC-estimated (REC) monthly R at MPB-06 for 2005 - 2016. The modelled contributions of the two components of R, heterotrophic (Rh) and autotrophic (Ra) respiration, are shown in green and red, respectively, and the total is the blue line. The EC-estimated values are the black circles including the 95% confidence interval represented by the vertical bars.  90   Figure 3.6  EC-estimated monthly (a) and annual (b) values versus modelled monthly (a) and annual (b) values of R at MPB-06 for 2007-2016 with the 1:1 line (solid black line) and regression line (dashed red line). See Table 3.3 for statistics.   Table 3.3 Slopes of the regression, intercepts, coefficients of determination (R2), root mean square errors (RMSE), p-values and sample sizes (n), for EC-estimated monthly and annual values versus modelled monthly and annual values of gross primary productivity (GPP), ecosystem respiration (R) and net ecosystem productivity (NEP) for 2007 - 2016 at MPB-06 shown in Figs. 6, 9 and 12.  GPP  R  NEP   Monthly Annual Monthly Annual Monthly Annual R2 0.89 0.86 0.82 0.73 0.72 0.61 RMSE (g m-2 of C) 17.49 59.64 17.21 57.94 9.18 40.79 n 120 10 120 10 120 10 Slope of the regression 0.91 0.98 0.99 0.95 0.80 1.01 91   GPP  R  NEP   Monthly Annual Monthly Annual Monthly Annual Intercept 6.59 37.07 2.94 58.14 0.42 0.42 p-value 0.003 0.677 0.251 0.632 0.623 0.977   Figure 3.7  Modelled (Rmod) and EC-estimated (REC) annual R at MPB-06 for 2005 - 2016. The red triangles show the modelled values, and the EC-estimated values are the blue filled circles. The red and blue shaded areas show, respectively, the standard deviation for annual modelled R, and annual EC-estimated R with the uncertainty estimates described in Section 3.2.3.  92  At the current site, there are no EC measurements before the beginning of the MPB attack, but the model indicates a decrease in R of 35% between 2005 and 2007 and recovery within a decade. Moore et al. (2013) found that R in MPB-attacked stands in Colorado did not increase on a decadal scale following tree mortality. First, R declined due to a decrease in new substrate, then it recovered 6-7 years after mortality, when soil organic matter increased due to leaf litter incorporation, followed by a decline in years 8-10 (Moore et al., 2013). This agrees with our modelling results showing little change between the pre-attack value and the value a decade following attack. Our measurements, however, show the increase in R occurring slightly later than observed by Moore et al. (2013), with the biggest increase 9-10 years after attack. Reed et al. (2014), using EC measurements in a beetle-attacked lodgepole pine stand in Southeast Wyoming over three years, also did not find a change in average weekly R with increasing mortality. Borkhuu et al. (2015) found that soil R at the stand level did not change over five years of beetle infestation, suggesting the resilience of the ecosystem as the remaining vegetation, tree growth and microbial respiration compensated for the expected decline in Ra due to tree mortality.   3.3.2.4 Gross primary productivity Monthly modelled GPP values at MPB-06 agreed well with EC-estimated values, generally falling within the uncertainty ranges of the measurements (Figure 3.8 and Figure 3.9a). 89% of the variance in monthly EC-estimated GPP was explained by the model with an RMSE of 17.49 g m-2 month-1 of C and a slope of the regression of 0.91, which was significantly different from 1.00 (p < 0.005) (Table 3.3). Only in the first three years following attack was 93  peak growing season GPP underestimated by the model. However, in 2014, when the summer was dry, the EC-estimated GPP values were lower than the modelled values. Annual modelled and EC-estimated GPP (Figure 3.10) also showed good agreement (R2 = 0.86, slope not significantly different from 1) (Figure 3.9b and Table 3.3) with mean values for the decade of 589 and 617 g m-2 year-1 of C, respectively. Modelled values suggested a significant decrease in GPP of 52% between 2005 and 2007 (p < 0.05). During the following decade, both EC-estimated and modelled GPP increased significantly (p < 0.001), recovering slowly to more than twice as high in 2016 as in 2007 (Figure 3.10). The modelled GPP values were lower than the measurements from 2007 to 2010 when the model seemed to be more sensitive to low θ during late summers, and to the continuing slight canopy mortality.   94   Figure 3.8  Modelled (GPPmod) and EC-estimated (GPPEC) monthly GPP values at MPB-06 for 2005 - 2016. The blue line shows the modelled values, and the EC-estimated values are the black circles with the 95% confidence interval represented by the vertical bars.  95   Figure 3.9  EC-estimated monthly (a) and annual (b) values versus modelled monthly (a) and annual (b) values of GPP at MPB-06 for 2007-2016 with the 1:1 line (solid black line) and regression line (dashed red line). See Table 3.3 for statistics.  96   Figure 3.10  Modelled (GPPmod) and EC-estimated (GPPEC) annual GPP values at MPB-06 for 2005 - 2016. The red triangles show the modelled values, and the EC-estimated values are the blue filled circles. The red and blue shaded areas show, respectively, the standard deviation for annual modelled and EC-estimated GPP with the uncertainty estimates described in Section 3.2.3.  Similar to our findings, Hansen (2014) concluded that recovery of C productivity occurred relatively quickly following MPB outbreaks returning to pre-attack levels within 5-40 years, as the amount of remaining live trees and secondary structure exceeds that left after other 97  disturbances such as wildfires and harvesting. Hansen et al. (2015) simulated GPP to be similar among infested and undisturbed stands within 100 years after MPB outbreak and in some cases even within 60 years. Using Landsat and MODIS satellite data, Bright et al. (2013) found reductions in GPP of 5-26% following tree mortality in northern Colorado, similar to the values of 15-20% and 13-30% reduction in GPP that Coops and Wulder (2010) and Moore et al. (2013) observed in BC and Colorado, respectively. These values are lower than our estimate of 52% reduction in annual GPP due to MPB attack. Edburg et al. (2011) found that GPP was reduced by 10-60% with mortality severities between 25% and 95%. The recovery of GPP following attack can also be observed in an increased mean radial growth of the understory due to increased light levels in central BC (Hawkins et al., 2013). They found that healthy overstory trees can also increase their growth to some extent following outbreak, but younger trees showed more growth response than older trees and it also varied by species. Romme et al. (1986) also found increased growth of surviving trees following MPB attacks in Wyoming and Montana with annual wood production returning to pre-attack levels within 10-15 years. A thinning and fertilization experiment in Oregon showed that increased light penetration caused greater photosynthetic and tree growth efficiency in the remaining trees after MPB attack (Waring and Pitman, 1985). Compensatory responses to beetle outbreaks were also observed in Colorado by Collins et al. (2011), who found that understory trees grew faster and more new seedlings were established following attack. Emmel et al. (2014) showed that the understory and secondary structure contributed significantly to the CO2 uptake during summertime following attack. Simulated aboveground C stocks in MPB-attacked lodgepole pine stands in central Idaho with mortality rates up to 52% decreased by 31-83% depending on number and size of killed and surviving trees, and then continuously increased reaching pre-attack levels within 25 years or 98  less (Pfeifer et al., 2011). In the Southern Rocky Mountains, Caldwell et al. (2013) found that simulated live C pools recovered within 40 years and that lodgepole pine remained the dominant species.  3.3.2.5 Net ecosystem productivity The time series of modelled and measured monthly NEP, which is the difference between GPP and R, is shown in Figure 3.11. Modelled NEP shows reasonable agreement with the measurements, both on a monthly and an annual basis with R2 values of 0.72 and 0.61, respectively (Figure 3.12 and Table 3.3). NEP reached values down to -16 g m-2 month-1 of C during the winter, a phenomenon consistent with continuing small respiratory losses.   99   Figure 3.11  Modelled (NEPmod) and measured (NEPEC) monthly NEP values at MPB-06 for 2005 - 2016. The blue line shows the modelled values and the measured values are the black circles with the 95% confidence interval represented by the vertical bars.  100   Figure 3.12  Measured monthly (a) and annual (b) values versus modelled monthly (a) and annual (b) values of NEP at MPB-06 for 2007-2016 with the 1:1 line (solid black line) and regression line (dashed red line). See Table 3.3 for statistics.  The modelling results suggest a decrease in annual NEP of 126% between 2005 and 2007 (Figure 3.13) with the stand changing from a C sink to a C source. Both modelled and measured annual NEP recovered during the decade following attack showing significant increases (p < 0.01). However, modelled NEP lagged a year behind the EC-measured NEP with the stand becoming a C sink in 2010 instead of 2009. Keeping in mind the relatively high uncertainty in NEP, it is reasonable to say that the stand became a distinct C sink in 2015 as the values in the 95% confidence interval for both measured and modelled NEP are positive. In 2007 and 2008, the modelled NEP suggests that the stand was a weaker C source than the measurements show, but both the model and the measurements show that the stand remained a weak C sink a decade following the attack. Annual EC-estimated and modelled values of GPP, R and NEP for 2005-2016 as well as the respective average values for 2007-2016 at MPB-06 are shown in Table 3.4. 101  The average annual modelled and EC-estimated GPP and R were within 4-5% of each other, and even the modelled and measured NEP agreed well with the average value of 22 g m-2 year-1 of C, respectively.   Figure 3.13  Modelled (NEPmod) and EC-estimated (NEPEC) annual NEP values at MPB-06 for 2005 - 2016. The red triangles show the modelled values, and the EC-estimated values are the blue filled circles. The red and blue shaded areas show, respectively, the standard deviation for annual modelled and measured NEP with the uncertainty estimates described in Section 3.2.3.  102  Table 3.4  EC-estimated and modelled annual gross primary productivity (GPP), ecosystem respiration (R) with the two modelled components, autotrophic (Ra) and heterotrophic (Rh) respiration, and net ecosystem productivity (NEP) for 2005 - 2016 at MPB-06. Year GPP (g m-2 year-1 of C) R (g m-2 year-1 of C) NEP (g m-2 year-1 of C)  EC-estimated Modelled EC-estimated Modelled R (Ra + Rh) Measured Modelled 2005  841  678 (446 + 232)   163 2006  606  567 (321 + 246)  39 2007 440 400 521 443 (212 + 231) -81 -43 2008 500 454 557 475 (240 + 235) -58 -22 2009 535 428 524 469 (227 + 242) 10 -41 2010 571 535 507 521 (283 + 238) 63 13 2011 567 594 559 547 (315 + 232) 8 47 2012 582 593 541 565 (314 + 251) 41 28 2013 654 684 593 633 (363 + 270) 61 52 2014 569 620 565 611 (329 + 282) 3 10 2015 837 778 794 695 (413 + 282) 43 84 2016 915 806 783 718 (427 + 291) 132 88 2007-2016 average 617 589 594 568 (313 + 255) 22 22  Reed et al. (2014) did not find a change in NEP over three years of beetle disturbance, as both GPP and R increased, and concluded that it was due to an increase in diffuse light within the canopy and an increase in the decomposition of soil organic matter, respectively. Hansen (2014), on the other hand, found a decrease in NEP following attack with recovery times of 5-40 years to return to pre-outbreak values, depending on outbreak severity. He concluded that moderately 103  infested stands can remain C sinks despite attack and severely infested stands returned to being C sinks within 5-20 years. This agrees with our observations that MPB-attacked stands became C neutral within 3-5 years after attack due to the remaining live trees and understory (Brown et al., 2012; Emmel et al., 2014; Bowler et al., 2012). As shown in Meyer et al. (2017) and this paper, not salvage harvesting beetle attacked stands with sufficient live trees and understory would be beneficial for C sequestration as well as for stand hydrology, because our observations showed a faster recovery than nearby clearcuts which could remain a C source for as long as a decade (Brown et al., 2010).   3.4 Conclusions 1. The modified version of 3-PG performed well in simulating monthly and annual stand-level R during the decade following MPB attack. On average, Ra represented 55% of R during 2007 to 2016. Ra significantly decreased in the first year following attack compared to pre-attack estimates, but Rh changed little. Thus, R decreased significantly, by about 35% in the first year, but in the following decade Ra, Rh, and hence R increased to estimated pre-attack values. Modelled R showed a slow, continual recovery, while EC-estimated R showed relatively little variation in the first eight years after attack and then increased in the last two years.   2. Treefall following canopy mortality due to MPB attack peaked eight years after the beginning of attack and caused a continual increase in the CWD C pool for the decade after attack. Rh soil was the dominant component of Rh, but Rh cwd increased significantly, making up about 11% of Rh in 2016. 104  3. Modelled and EC-estimated monthly and annual GPP showed good agreement with annual GPP increasing slowly, but significantly, throughout the decade after attack following an estimated initial decrease in modelled GPP of 52%. Its recovery was due to the growth of the remaining trees and of the understory, thought to be due to increased light penetration through the canopy and reduced competition for water and nutrients.  4. Modelled annual NEP decreased by about 126% between 2005 and 2007 with the stand becoming a C source, but both modelled and measured NEP increased significantly in the decade following attack with the model explaining 72% and 61% of the variance in measured monthly and annual NEP, respectively. This is remarkable agreement considering that the calculation of NEP combines the uncertainties of the two large terms, GPP and R. Both the model and measurements showed that the stand remained a weak C sink a decade following attack. 5. Overall, our results indicated that a not-salvage-harvested beetle attacked stand with sufficient live trees and understory showed a marked tendency to become a C sink within 3-4 years after attack thereby making this management practice beneficial for C sequestration.  105  Chapter 4: Synthesis of measurements of the carbon and water balances of lodgepole pine stands in the BC Interior following mountain pine beetle attack  4.1 Introduction 4.1.1 MPB attack history in North America Along with forest fires and harvesting, insect disturbances play an important role in determining the C balance of forest stands. In the BC Interior, the most recent MPB outbreak, which killed an estimated 731 million m3 of commercial timber over 18.3 million ha, equivalent to 54% of the mature merchantable pine volume in BC, raised the question of how to best manage the beetle-affected forests and how these forests would develop after the attack. It is estimated that ~55% of the mature merchantable pine volume will be killed by 2020, when the infestation in BC is expected to have subsided (BC Ministry of Forests, Lands and Natural Resource Operations, 2016). During this beetle outbreak, which started in the late 1990s, the beetle extended its historic range, quickly spreading into Northern BC, the Yukon, the Northwest Territories and east of the Rocky Mountains into Alberta (Cooke and Carroll, 2017). This raised concerns about the beetle spreading east into the boreal jack pine forest, which extends all the way to Eastern Canada (Nealis and Cooke, 2014). Safranyik et al. (2010) expected climate suitability for MPB range expansion to be moderate in central Alberta and Saskatchewan, where the pine volume is lower than in lodgepole pine stands in central BC and western Alberta, but possibilities for continued spread farther east increase again northwest of the Great Lakes with higher pine 106  volumes. Thus, there is a risk that persistent populations will become established outside the beetle’s historic range and spread farther east in the boreal forest under favorable weather conditions, when epidemic infestations can develop (Safranyik et al., 2010).  There are three different stages of attack. In the first stage, the green attack, adult beetles lay eggs beneath the bark which hatch after 1-2 weeks. The tree’s nutrient supply is cut off by the larvae feeding on the phloem (Taylor et al., 2006). A blue-stain fungus, introduced into the tree by the beetle kills living cells in the phloem and xylem and thus impacts the tree’s capacity to transport water as well as nutrients (Safranyik et al., 1974). In the second stage, the red attack, a year after infestation, needles turn red following the death of the trees. The last stage of attack is the grey stage, when the needles are decaying and falling off the trees (BC Ministry of Forests, Lands and Natural Resource Operations, 2017).  Attacked stands are managed following three different approaches, which are no harvesting, partial-salvage harvesting, and clearcut-salvage harvesting (Coates and Hall, 2005). To study the impacts of the beetle on C and water balances of pine forests, EC flux measurements have been made in two MPB-attacked lodgepole pine (Pinus contorta var. latifolia) stands in northern BC, between Prince George and Mackenzie, that were not salvage harvested, and one stand that was partial-salvage-harvested (i.e., only dead lodgepole pine trees were removed), complemented by short-term measurements in nearby clearcuts (Brown et al. 2010; Mathys et al. 2013; Emmel et al. 2013 and 2014). Measurements like these are important to verify model predictions for C losses due to MPB attack.  107  4.1.2 Effects of MPB attack on C and water balances A decrease in the rate of photosynthesis of the attacked stands due to tree mortality following MPB outbreak reduces C uptake thereby impacting the C balance (Brown et al., 2010). How NEP is impacted depends on the responses of Rh and Ra to canopy mortality. Modelling studies suggest that GPP decreases more than R, causing stands to switch from being C sinks to being C sources. Between 2000 and 2020, the MPB attack in BC was estimated to reduce NPP and increase Rh compared to the control scenario without beetle attack, causing a loss of 270 Tg of C over the 374,000 km2 of forest with salvage harvesting the beetle-killed and living trees causing the loss of an additional 50 Tg of C (Kurz et al., 2008). On the other hand, Williams et al. (2016) found that live biomass decreased by ~200 Tg year-1 of C in the US due to disturbances (harvesting, fire, windthrow, insect attacks and drought), while C stocks increased by ~190 Tg year-1 of C in the US and by 10-30 Tg year-1 of C in Canada, as growth was enhanced due to forest resilience and increased CO2 concentrations in the atmosphere. Arora et al. (2016), using the coupled CLASS-CTEM, predicted a sink of 900-1060 Tg of C for 1999-2020, depending on the future climate change scenario, due to changes in climate and increasing atmospheric CO2 concentration causing the enhanced C accumulation to far exceed the loss of 328 Tg of C due to the MPB attack in BC. Edburg et al. (2011), using the modified CLM with prognostic C and nitrogen cycling (Version 4) including snag and dead foliage pools, found that MPB attack reduced stand-level NPP for several years. Fluxes of C were impacted for decades and C stocks up to 100 years after attack, depending on outbreak severity, delay in snagfall, snagfall rate and management. 108  Outbreak duration affected C fluxes only in the short term, while the mortality rate had long-term effects on C fluxes and stocks (Edburg et al., 2011). Comparisons of a non-infested forest stand and an attacked forested valley ~30 km apart, using satellite estimates, CO2 observations (EC and atmospheric CO2 concentration measurements), soil chamber fluxes as well as soil organic C content measurements, confirmed that NEP was affected little by MPB attack on a decadal scale due to a similar decline in GPP and R (Moore et al., 2013). Unlike earlier modelling studies, Moore et al. (2013) found that decomposition of needle and root litter caused only a short-term (1-2 years following attack) increase in R, and that thereafter, R decreased due to reduced soil C inputs, suggesting that model estimates could be improved by further studies. Reed et al. (2014), using EC measurements in a beetle-attacked lodgepole pine stand in southeast Wyoming, found that an increase in canopy mortality from 30 to 78% over three years did not impact CO2 uptake and WUE of the stand, and attributed it to increased CO2 uptake and water use by the remaining trees and sub-canopy vegetation. This was also found in the BC interior, using EC (Emmel et al., 2013 and 2014) and component flux measurements (Bowler et al., 2012), where the broadleaf understory and smaller coniferous trees contributed significantly to the CO2 uptake and E during the summer a few years after the attack, while the canopy was a small source of C and water vapour. Another important concern is the effect of MPB attack on the stand hydrology. Due to tree mortality and a subsequent reduction in transpiration following beetle attack, an increase in runoff would be expected. Using an integrated hydrologic model, Mikkelson et al. (2013) found that MPB outbreaks resulted in decreased E, increased snow accumulation, as needle loss reduced the amount of P being intercepted by the canopy, earlier and faster snowmelt due to higher penetration of radiation through the canopy, and increased runoff. Increases in snow 109  accumulation were also found by Winkler et al. (2015), who showed that the snow water equivalent on April 1 in MPB-attacked forests in BC became more similar to that in clearcuts in the later stages of attack. The lag between MPB attack and needle loss during the red- and especially the grey-attack stage implies that the impact on hydrology in not-harvested beetle-attacked stands lags behind clearcut-harvested stands (Winkler et al., 2015). Maness et al. (2012) used BC Forest Survey annual stand density observations and MODIS monthly data over 170,000 km2 of affected forest in BC to quantify the MPB disturbance impact on the summertime surface energy budget. They found E was reduced, on average, by 19%, the sensible heat flux (H) increased by 8% and the outgoing radiative flux increased by 1% in beetle-attacked stands. Vanderhoof and Williams (2015) used current and historical MPB outbreak location datasets as well as MODIS and GMAO MERRA data to determine summertime E of lodgepole pine stands in the Colorado Rockies up to 60 years following MPB attack. Summertime E was only 6.3% lower 4-13 years after attack, 18.7% lower 14-20 years following attack, but it was 21.6% higher 30-40 years after attack than in control stands, which was correlated with stand density, LAI and the fraction of absorbed PAR (Vanderhoof and Williams, 2015). Hubbard et al. (2013) confirmed that beetle-attacked trees seem to die primarily due to blue stain fungi and decreased water flow from the soil to the canopy, rather than the beetles feeding on the phloem. A comparison of MPB attack and mechanical girdling on lodgepole pine trees in Colorado, USA (Hubbard et al., 2013), showed a rapid decline in transpiration in beetle-attacked trees, but no water supply limitation and continued growth during the growing season following treatment in girdled trees. Surviving trees might benefit from an increase in soil water 110  and nutrient (e.g., nitrate) availability due to the decline in transpiration, or it might lead to an increase in runoff (Hubbard et al., 2013).  4.1.3 Objectives The objective of this study was to synthesize the results of micrometeorological (EC) flux measurements made in MPB-attacked stands in the BC Interior over a decade. Specific objectives were to: (i) determine the effects of MPB attack on growing season and annual NEP, R, GPP, E and WUE in three stands, (ii) determine the rate of recovery, and assess the effects of different management strategies, specifically no harvesting, partial-salvage and clearcut-salvage harvesting, on NEP, E, P - E, as an estimate of runoff, and ΔS, and (iii) forecast the effects of MPB attack on GPP, R and NEP over the next decade in an 80-year-old lodgepole pine stand using a modified version of the 3-PG model.  4.2 Methods 4.2.1 Site descriptions Measurements took place at three locations in the northern interior of BC in the Sub-Boreal Spruce biogeoclimatic zone as described in Brown et al. (2010) and Mathys et al. (2013). Two of the sites, which are located approximately 35 km southeast of Mackenzie (and 120 km north of Prince George) at Kennedy Siding (MPB-06) and 70 km north of Prince George at Crooked River Provincial Park (MPB-03), were not harvested following MPB attack. The third site, which is approximately 46 km north of Prince George near Summit Lake (MPB-09), was 111  partial-harvested. MPB-06, which was attacked by the beetle in 2006 and not harvested, was an almost pure 15-m-tall (in 2007) lodgepole pine stand with an understory consisting mainly of pine seedlings and saplings on a gravelly sandy loam soil. MPB-03, which was attacked in 2003 and also not harvested, consisted of lodgepole pine (92%) and subalpine fir (Abies lasiocarpa) (8%) and had a developed secondary structure consisting of subalpine fir and white spruce (Picea glauca) saplings, pine seedlings and deciduous shrubs. This stand was also on a gravelly sandy loam soil. MPB-09, which was partial-harvested in 2009 following MPB-attack in 2006, originally consisted of 50% lodgepole pine and 50% black spruce (Picea mariana) and hybrid white spruce (Picea engelmannii x glauca) trees with a secondary structure comprising black spruce, hybrid white spruce and subalpine fir on a silty clay loam soil. Table 4.1 lists the stand characteristics of the three sites. As shown in Brown et al. (2012), in 2006 about 50% of the MPB-06 stand was in the green-attack stage and the remaining 50% healthy. As the percentage of trees attacked increased, only 14% of the canopy trees were healthy by August 2011. At MPB-03, the trees were already in the red- and grey-attack stage in 2007 and more than 95% of the canopy pine trees had died. Photographs of the sites, including canopy photographs, and instrumentation are shown in Appendix A.  Table 4.1  Stand characteristics at the three main sites MPB-06, MPB-03 and MPB-09.  MPB-06 MPB-03 MPB-09 Stand age (years) ~ 80 ~ 110 ~ 70 Tower location 55o06’42.8” N 122o50’28.5” W 54o28’24.8” N 122o42’48.4” W 54o13’25.4” N 122o36’53.5” W Elevation (m) 750 710 680 112   MPB-06 MPB-03 MPB-09 Canopy height (m) ~ 15 ~ 17 ~ 16 Stand density (height > 10 m) (stems ha-1) (2007, 2007, 2010) 1275 (204a) 558 (123) 435b Stand basal area (height > 10 m) (m2 ha-1) (2007, 2007, 2010) live: 11.8 – 19.2 dead: 0.2 – 0.9 live: 0.7 – 3.2 dead: 8.1 – 14.7 7.6c Seedling/sapling density  (stems ha-1) (2007 at MPB-06 and MPB-03 and 2010 at MPB-09) Pinus contorta: 7470 Abies lasiocarpa: 100 Picea glauca: 110 Pinus contorta: 2800 Abies lasiocarpa: 2300 Picea glauca: 190 replanted in May 2010: Pinus contorta: 781 Picea engelmannii x glauca: 589  Understory vegetation Alnus tenuifolia, Salix spp., Vaccinium spp. Salix spp., Amelanchier alnifolia, Vaccinium spp., Arctostaphylos uva-ursi Spiraea douglasii, Lonicera involucrata, Salix spp. LAI (overstory) (m2 m-2) 1.24 (in 2007) 0.95 (in 2010) 0.90 (in 2007) 0.60 (in 2010) 1.30 (in 2011) Litter-fibric-humus C content (kg m-2) 1.1 – 3.78 1.88 – 2.81 3.1 – 4.3 Mineral soil C content (0 – 55 cm) (kg m-2) 1.76 – 3.15 1.21 – 2.76 7.6 – 15.9 Fine soil bulk density (kg m-3) 1180 (220) 1160 (323) 1247 - 1495 Soil coarse fragments (% by volume > 2 mm) 34 (11) 70 (7) < 1 Soil texture Gravelly sandy loam Gravelly sandy loam Silty clay loam a Standard deviation in brackets. b Value for the tower footprint area. c Value for the full partial cut block (Nishio, 2010). 113   In addition to these three main sites, short-term measurements were conducted in three nearby clearcuts. CC-97 and CC-05 were clearcut-harvested in 1997 and 2005, respectively, and were located approximately 2 km SW and 1 km E of the MPB-06 flux tower, respectively. Measurements were conducted from 29 June to 23 July 2007 at CC-97 and from 24 July to 16 August 2007 at CC-05. At the time of the measurements, CC-97 was a 10-year-old clearcut, which was left to naturally regenerate. It had a large number of lodgepole pine seedlings (1200 stems ha-1) with a ground cover of lichens and mosses. CC-05 was salvage logged following MPB attack and planted in 2006 with a mixture of lodgepole pine and hybrid white spruce seedlings. The ground cover was similar to that of CC-97, but contained less lichen. Another clearcut stand MPB-09C, located nearby MPB-09, was clearcut-harvested in 2009 with measurements made during the summer of 2010.   4.2.2 Measurements At the three main tower sites MPB-06, MPB-03 and MPB-09 measurements have been made continuously since July 2006 to present, from March 2007 to December 2012 and from October 2009 to December 2012, respectively. At each of the three sites, a 32-m-tall scaffold tower (~2.1 m long x ~1.5 m wide) (see Figure A.2 in Appendix A) was erected on level ground (< 0.5o slope). An EC system consisting of a three-dimensional ultrasonic anemometer (model CSAT3, CSI, Logan, Utah) and an open-path IRGA (model LI-7500, LI-COR, Inc., Lincoln, Nebraska) (see Figure A.1) installed on the tower above the canopy at a height of 26 m measured NEP. High frequency (10 Hz) EC data as well as climate data were measured with two data 114  loggers (model CR1000, CSI) with synchronous-device-for-measurement (SDM) connections and stored on compact flash memory cards, which were exchanged every 4-6 weeks. Calculations of half-hourly covariances and averages performed by the data loggers were received at the UBC Biometeorology and Soil Physics Laboratory daily using a cell-phone connection. The systems were powered by three 100-W solar panels (CTI-130, Carmanah Technologies Corp., Victoria, BC) and an 800-Ah battery unit each as described in Brown et al. (2010).  Climate measurements included above-canopy down-welling and up-welling shortwave and longwave radiation (net radiometer, model CNR1, Kipp and Zonen B.V., Delft, The Netherlands) at the 30-m height and above-canopy down-welling and up-welling, and below-canopy down-welling photosynthetic photon flux density (quantum sensor, model LI-190 AS, LI-COR Inc.) at the 3-m height. At MPB-09, measurements of below-canopy up-welling shortwave radiation were made with a Black & White Pyranometer (model 8-48, Eppley Laboratory Inc., Newport, Rhode Island) as well. P (see List of Symbols) was measured at canopy height at MPB-06 and MPB-03, and at the 5-m height at MPB-09 with tipping bucket rain gauges (model 2501, Sierra Misco, Berkeley, CA at MPB-06, model TE525M, CSI, Logan, Utah at MPB-03 and model 52203 R.M. Young Co., Traverse City, MI at MPB-09). Wind speed (propeller-vane anemometer, model 05103 R.M. Young Inc.), Ta and RH (model HMP45C, Vaisala Oyj, Helsinki, Finland) were measured at the 25-m height. Ts measurements at MPB-06 and MPB-03 (chromel-constantan 30-gauge thermocouple wire, Omega Engineering Stamford, Connecticut) were made at depths of 2 (MPB-03 only), 5, 10, 20 and 50 cm, and at MPB-09, two soil profiles were made for measuring Ts at depths of 3, 10, 20, 50 and 100 cm using soil thermistors (model ST100, Apogee Instruments Inc., Logan, Utah). Also measured were G (three 115  heat-flux plates, model HP01, Hukseflux Thermal Sensors, Delft, The Netherlands at MPB-06 and MPB-03, 4 homemade heat flux plates using Peltier coolers at MPB-09) at a depth of 3 cm corrected for the rate of heat storage change in the 0-3 cm layer, and θ (soil water reflectometer, model CS 616, CSI) at the 0-10 cm and 30-50 cm depths at MPB-06 and (soil moisture probe, model EC-5, Decagon Devices Inc., Pullman, Washington) at the 10, 20 and 50 cm depths at MPB-03, and two soil profiles with depths of 3, 20, 50 and 100 cm at MPB-09 using water content reflectometers (model CS615 and CS616, CSI). At MPB-09, forest floor surface temperature was measured with a downward-facing infrared radiometer (model Apogee SI-111, CSI) and snow depth with an acoustic distance sensor (model SR50A, CSI) mounted at 3.6 m height. Snow depth was also measured at a clearcut weather station located ~ 1 km from MPB-06 (acoustic distance sensor, model SR50A, CSI) and P was calculated from these data and manual measurements of snow liquid water equivalent.  Canopy LAI at the three sites was measured using a LI-COR Plant Canopy Analyser (model LAI-2000, LI-COR Inc.) and hemispherical photography at MPB-06 and MPB-03. At 24 locations along 150-m long transects in four directions from the tower (NW, NE, SW and SE), hemispherical photographs were taken once a year from 2007 to 2015 at MPB-06 and from 2007 to 2013 at MPB-03 using a camera with a fisheye converter lens at the 1.2-m height (Foord and Spittlehouse, 2014). In the summer of 2007, 0.25-m x 3-m tall triangular towers were installed at CC-97 and CC-05 using the same EC instrumentation as on the main towers (see Figure A.6). Climate measurements made in these clearcuts were above-canopy (3-m height) down-welling and up-welling shortwave and longwave radiation (model CNR1, Kipp and Zonen B.V.), Ta and RH (model HMP45C, CSI) at the 3-m height, Ts (chromel-constantan 30 gauge thermocouple wire) 116  at 2, 10, 20 and 50 cm depths, G (3 heat flux plates) at 5-cm depth, and θ (model EC-5, Decagon Devices Inc.) at depths of 5, 15 and 30 cm. Between June 10 and September 22, 2010, EC measurements were made at MPB-09C with the same instruments as above mounted on a tripod at the 2.2-m height (see Figure A.7). A data logger (model CR5000, CSI) recorded EC and climate data at each of these clearcuts and the system was powered by a 100-W solar panel and a 100-Ah battery.   4.2.3 Flux quality control and data analysis As part of the quality control, data were excluded when more than 30% of a half-hour trace showed an instrument diagnostic warning flag. Measured CO2 and water vapour concentrations were limited by minimum and maximum values. As the fetch was greater than 1 km and 600 m, respectively, around the towers at MPB-06 and MPB-03, no data were rejected due to the wind direction. In the rare events, when the wind came through the tower, no effects were discovered. At MPB-09 on the other hand, data were rejected, when the wind came from the NE to SE (45 - 135o) as the footprint included areas outside the partial-harvested stand in those cases (Brown et al., 2010; Mathys et al., 2013). To ensure sufficient turbulent mixing, nighttime EC data were rejected when u* was lower than the threshold u* (u*th) of 0.30 m s-1 for MPB-06 and MPB-03 and of 0.20 m s-1 for MPB-09 (Brown et al., 2010; Mathys et al., 2013). To exclude data when wintertime CO2 uptake occurred due to the open-path IRGA heating up in cold air conditions (Burba et al., 2008), wintertime (defined as times when Ta < 5oC and Ts at the 5-cm depth < 1oC at MPB-06, and Ta < 5oC and Ts at the 5-cm and 3-cm depth < 1.5oC at MPB-03 and MPB-09, respectively) fluxes 117  with negative NEE were removed. As Burba et al. (2008) associated this issue with low wind speeds, we also rejected winter daytime data with wind speeds < 4 m s-1 (Brown et al., 2010). The fluxes of CO2 (Fc) and water vapour (E) were computed as the product of ρ and the covariances of the vertical wind speed (w) and CO2 and water vapour mixing ratios (s), respectively, i.e., CO2, 𝐹c =  𝜌 𝑤′𝑠𝑐′̅̅ ̅̅ ̅̅  and 𝐸 =  𝜌 𝑤′𝑠𝑣′̅̅ ̅̅ ̅̅ ̅, where the subscripts c and v refer to CO2 and water vapour, respectively, the primes are fluctuations from the mean and the overbar indicates half-hourly averaging. The wind data underwent three coordinate rotations, so that the mean vertical (?̅?) and lateral (?̅?) wind vector components became zero (Tanner and Thurtell, 1969). NEE was determined half-hourly as the sum of Fc and the rate of change of CO2 storage in the air column beneath the EC sensors (Sc). Sc was calculated as the product of ρ and the difference between 𝑠?̅? measured at the 26-m height for the following and previous half-hours divided by the difference in time between the measurements (i.e. 60 min) (Morgenstern et al., 2004).  NEP is then defined as –NEE and gaps in NEP of 2 hours or less are filled by linear interpolation. Nighttime R is obtained directly from nighttime NEE. The standard algorithm established by the Fluxnet Canada Research Network (FCRN) of the Canadian C Program (Barr et al., 2004) was used to obtain daytime R using an empirical logistic relationship between nighttime R and shallow Ts (5-cm depth at MPB-06 and MPB-03, and 3-cm depth at MPB-09), applied to daytime. An additional moving window parameter, rw (t), is included in the equation to estimate the seasonal variation of the other parameters giving   12 3ws( ) 1 ex ( )]p[r rr rtRT      (12) 118  where r1, r2 and r3 are empirical constants from the empirical logistic relationship. The slope of the linear regression between R obtained from the annual relationship and measured R within a moving window of 100 data points, moved in an increment of 20 points at a time, determined the time-varying parameter rw (t). The sum of daytime NEP and daytime R provided an estimate of GPP during the day. At nighttime as well as in the winter GPP was set to zero. GPP was modelled using the rectangular hyperbolic Michaelis-Menten relationship between GPP and the down-welling photosynthetic photon flux density fit to the non-zero values of GPP. As in the R calculation, a moving window parameter is included in the equation, thus giving   mamxw ax( )GPPp QAQ At         (13) where α is the quantum use efficiency and Amax is the ecosystem photosynthetic capacity (GPP at light saturation). Modelled GPP from Equation (13) was used to gap-fill GPP and daytime NEP was filled with the difference between GPP and R. As described in Brown et al. (2010), the moving window approach was not applied in the winter because the NEE measurements were sparse following quality control and a 100-point moving window could cover weeks or months instead of a few days.  Similar to the procedure of Amiro et al. (2006), gaps in growing season daytime E and H were filled using orthogonal linear regression equations between E and H and Rn - G obtained in 240 half-hour intervals, moved by 48 half hours at a time. Gaps in E at night were filled with zero. During the winter, E is close to zero and thus does not correlate well with Rn - G. To fill gaps in wintertime E, for each winter the ensemble averaged E was calculated for each half-hour. Remaining gaps in E due to the lack of good measurements for a particular half-hour were filled with linear interpolation. In-house software written in MATLAB® (Version 7.5.0, The Math 119  Works Inc., Natick, MA, USA) was used to perform all data calculations and analysis for this paper. All p-values resulted from t-tests performed at the 5% significance level.  4.2.4 Uncertainty analysis of EC fluxes The sum of the systematic error associated with the u*th selection and the random error, calculated as the square root of the sum of squares of the random error associated with the half-hourly NEP measurements and uncertainties in the gap-filling procedure, determined the uncertainty in annual NEP (Krishnan et al., 2006). Following Morgenstern et al. (2004), each half-hourly NEP value was assigned a 20% random error in order to determine the random error of annual NEP. A Monte Carlo simulation was used to obtain the uncertainty in gap-filling. Using a uniformly distributed random number generator, gaps in measured NEP with sizes of one half-hour up to 10 days were created and then gap-filled using relationships of R versus Ts and GPP versus PAR (Michaelis-Menten light response equation) (Krishnan et al., 2006). The 95% confidence interval was obtained by doing 500 simulations. Estimates of the systematic error due to the selection of u*th were determined by varying u*th by up to 50% and recalculating NEP (Morgenstern et al., 2004). Uncertainties in annual GPP, R and E were estimated the same way. In our calculations, we did not include other uncertainties such as those due to transducer shadowing of the CSAT3 sonic anemometer possibly causing underestimation of the vertical wind speed (Horst et al., 2015; Frank et al., 2016).  120  4.3 Results and Discussion 4.3.1 Climate The measured monthly climate variables at MPB-06, MPB-03 and MPB-09 for the decade 2007-2016 (Figure 4.1) show little variation in PAR, VPD and Ta between the three sites even though MPB-06 was 75 km north of MPB-09. P varied inter-annually, but was generally highest at MPB-03, apart from 2007, when it was higher at MPB-06. Measurements at MPB-03 only began in March of 2007, so the cumulative P does not represent the full year. Accordingly, θ was higher at MPB-03 than at MPB-06, but due to the difference in soil texture, MPB-09 showed values of θ on average 3-4 times as high as at MPB-03. The interannual variability in P showed 2007, 2011 and 2016 being relatively wet years and 2009, 2010 and 2012 being relatively dry (see Table 4.2). At MPB-06, θ was generally low in summer months, especially during 2009, 2010, 2014 and 2016, as well as in December 2010 to March 2011, and from December 2010 to March 2011, it was very low in the surface 10-cm layer (Figure 4.1) with values of 0.025 - 0.030 m3 m-3 and 0.045 -0.050 m3 m-3 in the 30-50-cm layer (not shown). This was likely due to soil freezing, as Ts at the 10-cm depth was between -0.45C and -0.70C from December 2010 to March 2011. Usually, monthly Ts at the 10-cm depth did not drop below -0.15C at this site, but even Ts at the 20-cm depth, and at times down to the 50-cm depth, was below zero during the 2010-2011 winter. At MPB-03 and MPB-09, monthly average Ts did not drop below zero, and only slight reductions in θ were observed during that time. PAR and VPD were low in June-July 2011, June 2012, August 2015 and April-June 2016, leading to annual PAR being below average in 2011 and 2016 (Table 4.2). VPD was high in July-August 2014, coinciding with low P and θ.  121    Figure 4.1  Monthly 24-h average photosynthetically active radiation (PAR), cumulative precipitation (P), monthly 24-h average vapour pressure deficit (VPD), monthly 24-h average air temperature (Ta) and monthly 24-h average volumetric soil water content (θ) for the 0-10 cm 122  layer for 2007 – 2016 at MPB-06 (black dashed-dotted line), MPB-03 (blue line) and MPB-09 (red dotted line). Values of VPD were daytime average values.  Table 4.2  Annual averages of photosynthetically active radiation (PAR), cumulative precipitation (P), vapour pressure deficit (VPD), air temperature (Ta) and volumetric soil water content (θ) for the 0-10 cm layer for 2007 – 2016 at MPB-06, MPB-03 and MPB-09. Values of VPD were calculated from daytime values.  2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 MPB-06           PAR (µmol m-2 s-1) 234.04 240.58 254.24 238.29 224.56 239.06 251.57 251.86 231.87 220.06 P (mm yr-1) 724.47 526.76 464.07 412.92 626.82 454.77 476.16 473.37 515.22 645.42 VPD (kPa) 0.53 0.53 0.59 0.60 0.46 0.60 0.61 0.64 0.58 0.57 Ta (°C) 2.97 2.88 2.43 3.96 3.18 3.19 4.87 4.10 5.09 4.94 Θ (m3 m-3) 0.11 0.10 0.08 0.08 0.08 0.09 0.09 0.07 0.08 0.08 MPB-03           PAR (µmol m-2 s-1) 226.93 253.83 269.94 256.58 238.38 258.25     P (mm yr-1) 442.00 709.43 484.30 472.70 673.60 516.10     VPD (kPa) 0.62 0.54 0.64 0.60 0.49 0.60     Ta (°C) 6.19 3.30 2.95 4.48 3.11 3.82     Θ (m3 m-3) 0.17 0.16 0.16 0.16 0.16 0.16     MPB-09           PAR (µmol m-2 s-1)    244.83 223.59 240.46     P (mm yr-1)    413.00 572.30 501.70     VPD (kPa)    0.65 0.48 0.60     Ta (°C)    4.86 3.06 3.71     Θ (m3 m-3)    0.49 0.56 0.59      123  4.3.2 Changing stand structure The original concept of choosing the three main sites was to have one stand with measurements before MPB attack, one stand following attack, which were not harvested, and one partial-harvested stand following attack. As shown in Figure 4.2, the beetle attacked MPB-06 earlier than expected with about 50% of the canopy being attacked in the year when measurements began. Thus, we do not have pre-attack measurements, but we can compare stands at different stages following attack. At MPB-06, the percentage of canopy trees in the red- or grey-attack stage increased quickly with only 14% of the canopy remaining healthy by August 2011. Since then, no further attack has been detected. At MPB-03, on the other hand, more than 95% of the canopy pine trees were already in the red- and grey-attack stage when measurements started in 2007.  124   Figure 4.2  Mountain pine beetle attack status at MPB-06 for August of 2006 – 2016. Values for 2006-2010 are from Brown et al. (2014) and values for 2010-2016 are from Dale Seip (personal communication). The bars represent the percentage of healthy (in black), green-attacked (in white) and red-/grey-attacked (in grey) canopy trees.  With the increase in mortality at MPB-06, the canopy opened up as can be seen in the canopy openness as well as the overstory LAI values in Table 4.3. The canopy openness increased from 30.8% in 2007 to 43.3% in 2015. Even though there was no new beetle attack since 2007 at MPB-03, its canopy openness increased from 39% in 2007 to 61% in 2013. This was likely due to continued needle fall and then tree fall in the later years.   125  Table 4.3  Measured canopy openness and overstory leaf area index (LAI) at MPB-06 and MPB-03 for 2007-2015 obtained from hemispherical photography.  2007 2008 2009 2010 2011 2013 2015 MPB-06        Canopy openness (%)a 30.8 31.4 33.0 31.3 38.6 39.5 43.3 Overstory LAI (m2 m-2)b 1.24   0.95  0.69 0.72 MPB-03        Canopy openness (%)a 39 47 48  48 61  a Values from Foord and Spittlehouse (2014) and Vanessa N. Foord (personal communication). b Values from Vanessa N. Foord (personal communication).  Rs transmitted through the canopy measured at a height of 1.2-m increased from ~36% and ~48% in 2007 to ~47% and ~60% in 2009 at MPB-06 and MPB-03, respectively (David L. Spittlehouse and Vanessa N. Foord, personal communication). Thus, the opening up of the canopy caused an increase in the fraction of Rs reaching the understory by about 10% between the first and third year following attack. The percentage of incoming PAR reaching below the canopy (3-m height) at the tower location, on the other hand, remained surprisingly constant during the first eight years following attack (Figure 4.3) with averages of 28%, 22% and 60% of incoming PAR reaching the understory at MPB-06, MPB-03 and MPB-09, respectively. It significantly increased in year 9 following attack at MPB-03 and a decade after attack at MPB-06 compared to the first year (p-values < 0.001 and < 0.03, respectively) with averages for the last two years reaching 51% and 33%, respectively. At MPB-09, there was no significant change in PAR below the canopy throughout the three years following harvest. At MPB-03, seven years following attack, the averages of July-August Rs and PAR at the 1-m height at three locations within 50 m from the tower (to capture variability in stand density) were 60% of the incoming Rs and PAR, respectively (Emmel et al. 2013, 2014). The increase in PAR below the canopy at 126  MPB-06 starting in year 9 after attack could be explained by treefall in the stand, which peaked 8 years after the attack (see Section 3.3.2.1). The differences in Rs and PAR reaching the understory at the three sites can be explained by their stand densities with MPB-09 having the lowest stand density. The increase in Rs reaching the 1.2-m height around the tower at MPB-06 in the first few years compared to the relatively constant PAR at the 3-m height on the tower, as well as the differences between the 1-m level PAR observations by Emmel et al. (2014) and our tower measurements at MPB-03, indicate that there was considerable spatial variability in the mortality and recovery rates within the stands. The monthly variation in the percentage of PAR reaching the 3-m height (Figure 4.3) is mainly due to changes in solar elevation.   127  Figure 4.3  Percentage of monthly incoming photosynthetically active radiation (PAR) reaching below the canopy (3-m height) at MPB-06, MPB-03 and MPB-09 for the decade following attack.   4.3.3 Water fluxes following MPB attack Time series of annual E for the decade following attack at MPB-06, MPB-03 and MPB-09 are shown in Figure 4.4. At MPB-06, annual E was conservative during the first five years following attack with up to 5% variation compared to the first year. Then it increased significantly (p < 0.01) over the following five years, being 51% higher a decade following attack than the first year after attack. Annual E at MPB-03 and MPB-09 did not show significant trends (p > 0.67), increasing by 9% and 2%, respectively, between the beginning and end of the respective observation periods with interannual variations up to 9% and 4%, respectively. The relatively constant annual E during the first five years after attack followed by increases in E during the next five years at MPB-06 were likely due to compensating effects of understory transpiration, transpiration by remaining live trees and soil evaporation as the stands recovered (Brown et al., 2014). Measurements by Emmel et al. (2013) and Bowler et al. (2012) confirmed that the understory was a significant contributor to summertime E. Also, our modelling results using 3-PG showed a significant increase in understory E at MPB-06 following attack (see Section 2.3.2.2).   128   Figure 4.4  Annual evapotranspiration (E) at MPB-06, MPB-03 and MPB-09 for the decade following attack. Bars show the 95% confidence interval for annual E.  The growing season cumulative daily P - E and cumulative daily ΔS at the three sites for the decade following attack are shown in Figure 4.5. At MPB-06, growing season P generally exceeded growing season E by 6-74 mm year-1, causing a small surplus of water, although in 2009 and 2014 E exceeded P by 4 mm and 70 mm, respectively. 2011 was a wet and cloudy year resulting in a water surplus of 206 mm at MPB-06. However, as soil water storage (S) decreased at the beginning of the growing season in 2011, and due to drainage below the root zone, this 129  annual water surplus did not cause a positive cumulative ΔS during the growing season. At MPB-03, growing season P was greater than growing season E by 0-23 mm in 2007 and 2010, and by 142-151 mm in 2008 and 2011. In 2009 and 2012, growing season E exceeded growing season P by 33-35 mm. At MPB-09, 2011 was the only year when P exceeded E during the growing season. In 2010 and 2012, growing season E exceeded growing season P by 24-36 mm. Annually, P always exceeded E at all three sites (Table 4.4). As 2007 and 2011 were exceptionally wet years, annual P exceeded annual E by more than 400 mm at MPB-06, so significant runoff would be expected. During the growing season, cumulative ΔS showed a decline in S in most months at all three sites, with recharge usually beginning in August or September. S at MPB-06 showed an especially large decrease in 2014 when E exceeded P for almost the whole growing season. At MPB-06, the decrease in S was generally greater than that at MPB-03. MPB-09 showed an exceptionally large decrease in S during 2010, with a maximum decrease of 354 mm, as P was low and E exceeded P during most of July, August and September. Beginning in September, S recovered again with ΔS at the end of 2010 being close to zero at MPB-09. As can be seen in Figure 4.1, there was a long almost rainless period at all three sites in 2010, resulting in a decrease in S, which was large at MPB-09 and MPB-03 but moderate at MPB-06, with the difference due to P – E during the 2010 growing season being negative for a shorter time at MPB-06 than at MPB-03 and MPB-09.  130   Figure 4.5  Cumulative values of daily precipitation (P) minus evapotranspiration (E) and daily change in root zone soil water storage (ΔS) at MPB-06, MPB-03 and MPB-09 for the growing seasons (May 1 – September 30) of the decade following attack. The soil water storage (S) corresponding to ΔS = 0 was 36.34 mm, 7.64 mm and 17.05 mm at MPB-06, MPB-03 and MPB-09, respectively, averaged over the observed years. For values of S for each year see Appendix C. Note that the y-axis ranges for P – E and ΔS are different for each stand.  131  4.3.4 Carbon balances The annual values of NEP, GPP and R for the three main sites for the decade after attack are shown in Figure 4.6. There was an increase in NEP at the three sites with the not-harvested sites gaining C neutrality within three to five years after attack, and the NEP of the partial-harvested stand, MPB-09, not reaching C neutrality during the three years of measurements. In general, the rate of recovery was slow, as GPP at MPB-06 and MPB-03 plateaued between 4 to 8 years following attack. At MPB-06, annual GPP increased again in the last two years, leading to a significant (p < 0.001) increase in GPP between 2007 and 2016, being twice that one year following the attack. At MPB-03, GPP was 23% higher in 2012 than in 2007. At MPB-09, GPP increased by 14% higher from 2010 to 2012.  132   Figure 4.6  Annual net ecosystem productivity (NEP), gross primary productivity (GPP), and ecosystem respiration (R) at MPB-06, MPB-03 and MPB-09 for the decade following attack. Bars show the 95% confidence interval for annual NEP, GPP and R. The years when MPB attack occurred were 2003 (estimated), 2006 and 2006 (estimated) for MPB-03, MPB-06 and MPB-09, respectively.  Surprisingly, R has remained relatively constant since the attack, being equal in 2007 and 2012 at MPB-03. Only in the last two years has R at MPB-06 increased, resulting in a significant (p < 0.01) increase of 50% by the end of the decade following attack. The number of fallen trees 133  at MPB-06 has continuously increased throughout the decade, but as decomposition of fallen trees occurs over many decades (Hansen, 2014), there has been little evidence of a significant increase in R due to decomposition of these trees. The increase in R nine to ten years following attack at MPB-06 coincides with a significant increase in GPP and thus can be primarily explained by an increase in Ra rather than Rh. While 3-PG modelling results suggest a significant increase in Rh of CWD (Rh cwd), coming from the fallen dead trees, during the decade after attack (see Section 3.3.2.2), it comprises only a small proportion of total R (~4% in 2016) (see Table 3.4). At MPB-09, R was only 6% higher in 2012 than in 2010. Annual 3-year average R at MPB-09 was 81% and 98% higher than at MPB-06 and MPB-03, respectively. Similarly, annual 3-year average GPP was 57% and 73% higher at MPB-09 than at MPB-06 and MPB-03, respectively. The higher GPP and R indicate that MPB-09 is a more productive stand than MPB-06 and MPB-03, which have relatively nutrient-poor, coarse-textured, but well-drained soils. MPB-09 is on a silty clay loam soil, which has higher water holding capacity and thus a higher θ, as well as higher available θ than at the other sites. Reed et al. (2014) also found that average weekly R remained relatively constant over three years when mortality increased in a lodgepole pine stand in Southeast Wyoming. In MPB-attacked stands in Colorado, Moore et al. (2013) found a decrease in R in the first few years, then a recovery 6-7 years after mortality followed by a decline after 8-10 years, resulting in little change in R over a decade. Table 4.4 and Table 4.5 show the annual as well as the growing season totals of NEP, GPP, R and E at all the sites for the decade following attack. During the growing season, the not-harvested stands as well as the partial-harvested stand remained C sinks, despite beetle attack, unlike the clearcuts, which were C sources (Table 4.5). Both annually and for the growing 134  season, MPB-06 showed significant increases in GPP, R, NEP and E from 2007 to 2016 with p-values < 0.001 for GPP, p < 0.01 for NEP and R, and p < 0.01 and p < 0.02 for growing season and annual E, respectively. At MPB-03, only NEP increased significantly from 2007 to 2012 with p < 0.05 and p < 0.04 for the growing season and annual values, respectively. Both annually and for the growing season, none of the fluxes (GPP, R, NEP and E) showed a significant increase at MPB-09.   Table 4.4  Annual totals of net ecosystem productivity (NEP), gross primary productivity (GPP), ecosystem respiration (R), evapotranspiration (E) and precipitation (P) - E at MPB-06, MPB-03 and MPB-09 for 2007-2016.     2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 MPB-06           NEP  (g C m-2 year-1) -81 -58 10 63 8 41 63 3 43 132 GPP (g C m-2 year-1) 440 500 535 571 567 582 656 569 837 915 R (g C m-2 year-1) 521 558 524 507 559 541 593 566 794 783 E (mm year-1) 226 236 233 229 217 254 280 233 346 342 P-E (mm year-1) 498 291 231 183 410 201 196 241 169 303 MPB-03                     NEP  (g C m-2 year-1) -53 5 11 -12 50 48         GPP (g C m-2 year-1) 446 518 514 488 523 548     R (g C m-2 year-1) 499 513 504 499 473 500         E (mm year-1) 284 306 282 287 281 309         P-E (mm year-1) 158 403 202 185 393 207     MPB-09                     NEP  (g C m-2 year-1)       -108 -53 -53         GPP (g C m-2 year-1)    813 962 926     R (g C m-2 year-1)       921 1015 979         E (mm year-1)       316 330 323         P-E (mm year-1)       97 242 178         135   Table 4.5  Growing season (GS, defined as 1 May – 30 September) totals of net ecosystem productivity (NEP), gross primary productivity (GPP), ecosystem respiration (R) and evapotranspiration (E) at MPB-06, MPB-03, MPB-09 and MPB-09C for 2007-2016 and average daily NEP and E at CC-97, CC-05 and MPB-06 for 29 June - 23 July 2007 and 24 July - 16 August 2007, respectively. Values shown for MPB-06 for CC-97 and CC-05 are for the same period the clearcut measurements were made.     2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 MPB-06           NEP  (g C m-2 GS-1) 13 31 89 95 55 60 101 99 84 151 GPP (g C m-2 GS-1) 388 448 477 483 494 510 553 519 693 781 R (g C m-2 GS-1) 375 417 389 389 439 450 452 420 610 631 E (mm GS-1) 170 172 181 171 177 205 222 174 279 274 MPB-03                     NEP  (g C m-2 GS-1) 21 60 50 42 71 86         GPP (g C m-2 GS-1) 386 448 437 389 454 473     R (g C m-2 GS-1) 365 388 387 347 383 387         E (mm GS-1) 230 235 223 203 224 247         MPB-09                     NEP  (g C m-2 GS-1)       12 58 80         GPP (g C m-2 GS-1)    685 855 851     R (g C m-2 GS-1)       673 796 771         E (mm GS-1)       248 270 261         MPB-09Ca           NEP  (g C m-2 GS-1)       -103             GPP  (g C m-2 GS-1)    294       R (g C m-2 GS-1)       397             E (mm GS-1)       197             29 June - 23 July 2007                       CC-97 MPB-06               NEP  (g C m-2 day-1) -0.37 -0.07         136      2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 E (mm day-1) 3.97 3.16                 24 July - 16 August 2007                       CC-05 MPB-06               NEP  (g C m-2 day-1) -0.87 0.40                 E (mm day-1) 4.64 2.87                  aThe growing season values at MPB-09C are for June 10-September 22.   At the two clearcuts nearby MPB-06, CC-97 and CC-05, measurements were made for just over 3 weeks each during the 2007 growing season (Table 4.5). Between 29 June and 23 July 2007, the 10-year-old clearcut had an average daily (24-hour) NEP of -0.37 ± 0.20 g m-2 day-1 of C, while the not-harvested stand MPB-06 was close to C neutral with an NEP of -0.07 ± 0.35 g m-2 day-1 of C. During the period from 24 July to 16 August 2007, average daily NEP at CC-05 and MPB-06 were -0.87 ± 0.13 and 0.40 ± 0.41 g m-2 day-1 of C, respectively (Brown et al., 2010). Thus, both the two- and the ten-year-old clearcut stands were estimated to be growing season as well as annual C sources; and as measurements were made around the time of peak photosynthetic activity, they both were stronger sources than the not-harvested stand during the first year after attack. Average daily E between 29 June and 23 July 2007 at CC-97 and MPB-06 was 3.97 and 3.16 mm day-1, respectively, whereas during 24 July to 16 August 2007, E at CC-05 and MPB-06 was 4.64 and 2.87 mm day-1, respectively. Average daily PAR, Ta and total P at MPB-06 were 460 and 405 µmol m-2 s-1, 16 and 14oC and 35 and 32 mm during the two above noted time periods, respectively. The higher values of E at the clearcuts than at MPB-06 were likely due to a higher soil evaporation rather than transpiration, especially at CC-05, where there was only sparse vegetation consisting of lichen, moss and planted seedlings of lodgepole pine and hybrid white spruce two years following harvest (Brown, 2011). 137  The original hypothesis was that there would initially be a strong negative NEP (i.e., a strong CO2 source) in the not-harvested stands. Instead, our measurements show that the stands were weak C sources following MPB attack, and recovery was faster than expected with the not-harvested stands becoming C neutral only 3-5 years after attack (Figure 4.6). The partial-harvested stand was a C source, but a much weaker source than clearcuts in the area, which remained C sources for up to a decade (Brown et al., 2010). Tree physiological measurements (foliar net assimilation) (Bowler et al. 2012) and within-canopy EC flux measurements at several heights (Emmel et al. 2014) have shown that the secondary structure, particularly seedlings, saplings and immature trees, as well as the understory, consisting mainly of broadleaf shrubs, contributed to the rapid recovery at the not-harvested sites. This is partly due to the response of this vegetation to an increase in within-canopy PAR. Hawkins et al. (2013) observed that increased light levels caused the understory, and, to some extent, the remaining overstory, to increase its mean radial growth following beetle attack in lodgepole pine stands in central BC. Faster understory growth as well as an increased number of natural seedlings established after attack were also found in lodgepole pine stands in Colorado by Collins et al. (2011). Modelling studies such as those by Kurz et al. (2008) and Edburg et al. (2011) had predicted MPB-attacked stands to reach C neutrality after 10-30 years depending on the severity of attack and snagfall timing. Thus, MPB-06 and MPB-03 recovered faster than forecasted by these studies. Hansen (2014) predicted GPP to return to pre-attack levels within 5-40 years, as MPB outbreaks leave more live trees and secondary structure than disturbances like wildfires and harvesting. The monthly means of GPP, E and WUE at the three main sites for April to October during their respective measurement periods are shown in Figure 4.7, showing their interannual 138  variability. In general, GPP and E peaked in July at all three sites with the highest values at MPB-09. The differences between MPB-09 and the two not-harvested sites were larger for GPP than E. In 2010, MPB-03 and MPB-09 were more affected by the drought in July and August than MPB-06. At MPB-03, there were also noticeable effects of drought in 2009. At MPB-06, GPP was the lowest in 2007, especially in June and July, and 2016 had the highest GPP during the first half of the growing season. WUE varied throughout the growing season but was usually lowest at MPB-03 with average values for the April to October period of 2.6, 1.8 and 2.8 g of C (kg of H2O)-1 at MPB-06, MPB-03 and MPB-09, respectively. Annually, the average WUE was 1.7, 1.3 and 2.0 g of C (kg of H2O)-1 at MPB-06, MPB-03 and MPB-09, respectively. At MPB-09, WUE was reduced in June, July and August of 2010 as GPP was affected more by the drought than E.  139   Figure 4.7  Monthly, April to October, means of EC-derived gross primary productivity (GPP), evapotranspiration (E) and water use efficiency (WUE) for 2007 - 2016 at MPB-06, for 2007 – 2012 at MPB-03 and for 2010 – 2012 at MPB-09.  At all three sites, the relationships between monthly GPP and E were strong as can be seen in the linear regression equations in Figure 4.8. E explained 94%, 96% and 96% of the variance in GPP at MPB-06, MPB-03 and MPB-09, respectively, with RMSEs of 12.51, 8.94 and 15.87 g m-2 140  month-1 of C, respectively. With the RMSE lowest at MPB-03, there was the least spread in WUE among different years.   Figure 4.8  Monthly gross primary productivity (GPP) versus evapotranspiration (E) for 2007 - 2016 at MPB-06 (a), for 2007 – 2012 at MPB-03 (b) and for 2010 – 2012 at MPB-09 (c) 141  including the respective regression lines (dashed red line). Also shown are the linear regression equations, coefficients of determination (R2), root mean square errors (RMSE) and sample sizes (n) obtained from EC measurements.  Our study has shown that a faster C recovery and minimal changes in stand hydrology can be achieved, if MPB-attacked stands with sufficient remaining live vegetation are not clearcut-harvested following attack. Depending on the stand density and tree species composition of the remaining live trees, Coates and Hall (2005) found that the expected timber supply shortage following beetle attack can be reduced, if affected stands are not salvage harvested or if non-pine trees remain undamaged after harvest. Focusing on stands with at least 50% of the canopy and sub-canopy basal area being lodgepole pine, they concluded that a mid-term (in 15-50 years) harvest is likely possible in the case of stands with a secondary structure basal area of at least 5-10 m2 ha-1, which was the case in ~20-35% of stands in north-central BC. Coates et al. (2006) determined that in ~60% of the pine-leading stands in the Prince George Forest District, the basal area of the secondary structure is at least equivalent to that of an average 20-year-old spruce plantation. Thus, apart from the possible advantages of increased C sequestration in not-harvested compared to clearcut-harvested stands, there can also be an economic advantage, if stands can be harvested in the mid-term. In addition, not-harvested stands can offer ecosystem services, such as winter habitat for Woodland Caribou at MPB-06 (Seip and Jones, 2009).  142  4.3.5 Modelling predictions of R and NEP at MPB-06 with 3-PG As described in Chapters 2 and 3, I have run a modified version of 3-PG for MPB-06. The model was used to make predictions for the upcoming decade using five different climate scenarios (see Section 2.2.6). These scenarios were: S1, the monthly average of the climate measurements from January 2007 to December 2016 used for 2017-2026; S2, the monthly measurements from January 2007 to December 2016 repeated for 2017-2026; S3, the same as S1 but with an increase in Ta of 0.02 K per year and an increase in P of 0.1% per year; S4, the same as S3 but Ta increased by 0.05 K per year; and S5, the same as S4 but with an increase in VPD of 0.25% per year and an increase of Rs of 0.05% per year during the summer and a decrease of Rs of 0.03% per year during the winter months. Annual modelled GPP, R and NEP at MPB-06 for 2017-2026 for the five climate change scenarios is shown in Figure 4.9. Similar to the projections of E, GPP and WUE shown in Figure 2.12, S2, which accounted for the current interannual climate variability, showed the most variable annual R and NEP. Between 2017 and 2026, S2 showed an increase in GPP of 9%, in R of 10% and in NEP of 6%. The other scenarios, S1 and S3-S5, showed changes in GPP of -1% to +1%, in R of 0% to 1% and in NEP of -14% to +3% between 2017 and 2026. Thus, S2 suggests an increasing trend in GPP and R and little change in NEP, but these trends are not significant. The other scenarios indicate that GPP, R and NEP will remain relatively constant throughout the upcoming decade. Compared with the estimated pre-attack values in 2005, the average values for 2026 of S1 and S3-S5 show a reduction in GPP, R and NEP of 15%, 3% and 65%, respectively. S2 projected values show a reduction in GPP and NEP of 3% and 43%, respectively, but R increases by 7% from 2005 to 2026. Using 2026 values from the regression equations in Figure 4.9 for S2, the projected reductions in GPP and NEP from 143  2005 to 2026 of 12% and 65%, respectively, are similar to the results using the other climate change scenarios. For R, the regression equation still shows a slight increase of 1%. The relatively large reduction in GPP and NEP from 2005 to 2026 could be partly due to a possible overestimation of GPP and NEP just before the attack in 2005. It also indicates, however, that an MPB-attacked stand with only 14% of the canopy remaining healthy following attack might require 30-50 years to grow into a productive stand.  144   145  Figure 4.9  Annual modelled gross primary productivity (GPP) (a), ecosystem respiration (R) (b), and net ecosystem productivity (NEP) (c) at MPB-06 for 2017-2026 for the five climate change scenarios (S1-S5). For S2 the black circles are the annual values and the dashed black line is the regression line. Also shown are the linear regression equations for S2.  4.4 Conclusions This study examined the recovery of MPB attacked stands in the BC interior at different stages after attack and under different management strategies.  1. Not-harvested lodgepole pine stands were relatively resilient following attack, becoming weak annual C sources immediately after the attack but remaining C sinks during the growing season. They recovered to C neutrality within 3-5 years with GPP increasing and R showing little change in the first 5 years after attack, and remained weak annual C sinks thereafter.  2. Salvage harvested stands are less resilient with the partial-harvested stand remaining a C source three years following harvesting. The partial-harvested stand had a lower NEP than the not-harvested stands, but was a weaker C source than the clearcuts, which remained growing season C sources for up to 10 years following harvest. MPB-09 was a more productive stand with higher GPP and R than the not-harvested stands. 3. E was conservative after beetle attack and only showed a significant increase at MPB-06 between year 6 and 10 following attack.  4. The secondary structure and understory in the not-harvested stands significantly contributed to stand recovery by maintaining GPP and E. The understory benefitted from the higher 146  transmission of Rs and PAR due to increasing canopy openness of the stands following attack. Sheltering by the trees might have increased Ta in the understory. 5. Due to interannual climate variability, at all three sites, P exceeded E during some growing seasons and E exceeded P in others, but on an annual basis P was always greater than E, resulting in a generally moderate water surplus. S usually decreased during the growing season, with the decrease generally being greater at MPB-06 than at MPB-03, but recovered beginning in August or September. 6. Projections using the modified version of the 3-PG model for five climate change scenarios for the upcoming decade (2017-2026) suggest no significant trends in GPP, R and NEP at the not-harvested stand MPB-06. The average values in 2026 for the different climate change scenarios show a reduction in GPP, R and NEP of 14%, 1% and 65%, respectively, compared to estimated pre-attack values in 2005, indicating that an MPB-attacked stand with only 14% of the canopy trees remaining alive will need to grow for several decades for stand productivity to return to pre-attack levels.      147  Chapter 5: Conclusions  In this thesis, I examined eddy-covariance (EC)-measured carbon (C) and water fluxes from several mountain pine beetle (MPB)-attacked stands in the BC Interior to determine the rate of recovery of these stands following attack, and the effects of different management strategies on recovery. The studied stands included two not-harvested stands dominated by lodgepole pine at different stages after attack (MPB-06 and MPB-03), a mixed-conifer partial-salvage-harvested stand (MPB-09), and three clearcut-harvested stands (CC-05, CC-97 and MPB-09C) near two of the main stands. My research also involved the use of the process-based 3-PG model, which I modified to include the effects of insect attack on stand growth and ran for an 80-year-old lodgepole pine stand in BC, attacked in 2006 (MPB-06), for 20 years following attack. The modified version of 3-PG, applicable to beetle-attacked lodgepole pine stands, included a two-layer canopy with a partly dying overstory and a growing understory, water availability from snowmelt, changes in solar radiation reaching the understory following attack, and a heterotrophic respiration (Rh) sub-model to estimate soil and coarse woody debris (CWD) Rh, and to obtain stand net ecosystem productivity (NEP). The C and water fluxes in the first decade after the attack were compared with EC measurements at MPB-06. The major conclusions of this thesis are: 1. Not-harvested lodgepole pine stands in the interior of BC with canopy mortality rates above 80%, which became weak annual C sources immediately after the attack, were more resilient than salvage-harvested stands, recovering to C neutrality within 3-5 years and remaining weak C sinks thereafter, showing the benefit of this management practice for C sequestration. 148  The partial-harvested stand, which was monitored for only three years, was more productive with higher gross primary productivity (GPP) and ecosystem respiration (R) than the not-harvested stands, remained an annual C source three years following harvesting with a lower NEP than the not-harvested stands. However, clearcuts in the area remained C sources for up to 10 years after harvest. The non-harvested and partial-harvested stands remained growing season C sinks following MPB attack, while the clearcuts were growing season C sources.  2. In the not-harvested stand, MPB-06, GPP increased slowly, plateauing between years 4 to 8 following attack and then increased significantly (p < 0.001) in the last two years, doubling from one year after attack. R changed little in the first 8 years following attack and then increased in the last two years, resulting in a significant (p < 0.01) increase of ~50% above that observed one year after the attack, coinciding with the increase in GPP.  3. In the three main stands, evapotranspiration (E) was conservative after beetle attack and only showed a significant increase at MPB-06 in the second half of the decade following attack. On an annual basis, precipitation (P) was always greater than E, resulting in a generally moderate water surplus at these sites. During the growing season, P - E, an estimate of runoff, showed interannual variability with P being larger than E in some growing seasons and E being larger than P in others. The decrease in soil water storage (S), generally occurring during the growing season, was greater at MPB-06 than at MPB-03, but S began to recharge in August or September at all three sites. 4. Throughout the decade following attack, treefall at MPB-06 increased, however the decomposition of fallen trees was slow and did not cause a significant increase in R. As the canopy openness of the stands increased following attack, more solar radiation was transmitted through the overstory canopy, resulting in a significant contribution by the 149  secondary structure and understory in the not-harvested stands to stand GPP and E and thus to stand recovery. 5. During the nine years following attack at MPB-06, 3-PG performed well in modelling monthly and annual stand-level E. Measured annual E was remarkably stable in the first five years after attack, and then increased somewhat in the last four years, whereas modelled annual E decreased by ~62% in the first year after attack compared to pre-attack estimates and then recovered slowly. Estimates of root zone drainage (Dr) using measured P, E and S agreed well with modelled Dr, showing a bimodal distribution with maxima in the spring and fall. Similar to the measurements, modelled annual E and Dr accounted for 37% and 65% of annual P on average. 6. There was good agreement between modelled and EC-estimated monthly and annual GPP with annual GPP increasing slowly throughout the decade after attack, due to the growth of the remaining live trees and understory. Modelled GPP showed a decrease of 52% compared to pre-attack estimates in the first year after the attack. There was a slight decrease in annual modelled WUE over the 12 years, especially during the decade from 2007-2016, while EC-estimated WUE increased from its lowest value in 2007 up to its highest value in 2016, as EC-estimated GPP increased more than measured E. 7. The model simulated monthly and annual stand-level R well, suggesting a significant decrease in R of ~35% in the first year after attack, followed by increases in autotrophic respiration (Ra), Rh, and hence R returning to estimated pre-attack values nine years after attack. While EC-estimated R increased significantly in the last two years only, modelled R increased continually over the decade. Modelled Ra, representing 55% of R on average, significantly decreased in the first three years after attack, but Rh changed little. During the 150  decade following attack, the CWD C pool increased due to treefall of the killed canopy trees, which peaked eight years after the attack started at MPB-06. This caused CWD heterotrophic respiration (Rh cwd) to increase significantly, although soil heterotrophic respiration (Rh soil) remained the dominant component being about 89% of Rh in 2016. 8. In the first year after attack, modelled annual NEP for MPB-06 decreased by ~126% compared to the estimated value in the year before the attack (i.e., making it a C source), but increased throughout the following 9 years, making the stand a C sink within 4 years after attack. The model explained 72% and 61% of the variance in measured monthly and annual NEP, respectively, showing remarkable agreement taking into account the uncertainties of the two large terms, GPP and R, which determine NEP. 9. The modified 3-PG parameter sensitivity analysis showed the robustness of the model, because E and GPP showed similar sensitivities to parameters with modelled E and GPP being highly sensitive to quantum use efficiency of photosynthesis (α) and stomatal response to vapour pressure deficit (sD), and not so sensitive to boundary layer conductance (ga) and maximum rainfall interception. 10. The modified 3-PG projections over the following decade (2017-2026) made for five different climate change scenarios showed decreases, averaged for the different climate change scenarios, in GPP, R, NEP, E and WUE of 14%, 1%, 65%, 5% and 10%, respectively, in 2026 compared to pre-attack values in 2005. This suggests that in a stand where 86% of the canopy trees died due to beetle attack, young trees need to grow substantially for stand productivity, and thus GPP and NEP, to return to pre-attack levels.   151  5.1 Limitations of the research The 3-PG model, which is run on a monthly time step and uses monthly data as inputs, has limitations, as it cannot assess shorter time scale hydrological processes, for example. A modified water balance sub-model with a daily time step and a two-layer soil profile, enabling the calculation of canopy transpiration and soil evaporation, was developed by Almeida and Sands (2016), resulting in improvements of modelled volumetric soil water content () and available  (a). As this sub-model required additional data on soil physical properties and daily climate data, which are not generally easily available, the original 3-PG water balance sub-model was used in this research, making the data demand more manageable for forest managers or other users and the model still produced reasonable water balance results. Another limitation of the study is the simplicity of the Rh sub-model, which considers only soil temperature as a controlling factor and not , substrate quality following beetle attack or other factors. To determine the base-respiration rates of the soil and CWD, short-term chamber measurements of Rh were used, which allowed the two components of Rh to be estimated. Longer-term measurements of Rh could improve model estimates and verify Rh from the two substrates, as well as the contribution of Rh to R. In this thesis, only total R, not its components, could be compared to measurements over the decade following attack.  This issue of model verification also applies to the contributions of the overstory and understory, respectively, to ecosystem fluxes. Long-term measurements were only available for the total ecosystem fluxes, but measurements of overstory LAI, compared to the modelled values in Figure 2.4, could provide an estimate of model accuracy for the overstory component. 152  In this study, the modified model has only been applied to a single-species stand and has yet to be tested on mixed species MPB-attacked stands. Some more modifications would be required, if other species were added to the overstory or understory components. Another limitation of this research is that no pre-attack flux measurements were available, as the beetle attacked MPB-06 earlier than expected. Measurements started during the year of attack, but pre-attack values could only be estimated using 3-PG. However, the measurements provide valuable information about the fluxes during the decade following attack.  5.2 Future outlook Given there are still uncertainties regarding the effects of forest disturbances, such as insect outbreaks, on the C and water balances of stands following attack, measurements like those made in this study can contribute to our understanding of forest disturbance effects and the improvement of ecosystem models. MPB attacks are not a local issue, as the beetle’s range extends throughout Western North America, and models such as this modified version of 3-PG could also be applied to stands attacked by other insects causing tree mortality. Under future climate change, the rate and intensity of insect outbreaks is expected to increase due to higher temperatures (Williams et al., 2016) and possibilities exist for higher survival rates of beetles if winter temperatures do not drop low enough to cause the beetles to die. With the likelihood of increasing temperatures, there are also concerns that the MPB might be able to increase its range farther and spread into the boreal forest and eventually into Eastern Canada (Nealis and Cooke, 2014; Safranyik et al., 2010).  153  Further observations or modelling studies, such as this study, informed by EC measurements, can advise scientists as well as forest managers about disturbance effects and provide information about suitable management strategies following attack. It is valuable to obtain measurements in stands with different mortality rates, stand structure and in different climates in order to get an overview of differences in stand responses to MPB attack, and to verify model estimates. Early modelling studies, such as Kurz et al. (2008) and Edburg et al. (2011), predicting stands to return to C neutrality within 10-30 years after attack depending on the severity of attack and timing of treefall, might have overestimated the damaging effect of beetle attack. Williams et al. (2016) and Arora et al. (2016) suggest that increasing atmospheric CO2 concentration together with forest resilience will cause an increase in C stocks and NEP, respectively, between 2000 and 2020 exceeding the reductions due to disturbances in the US and the MPB attack in BC. Since more live trees and secondary structure remain after beetle attack compared to fires and harvesting, Hansen (2014) predicted stand productivity, i.e. GPP, to return to pre-attack levels within 5-40 years. In addition to our study showing the benefit to the C balance of not-harvesting stands immediately following attack instead of clearcut-harvesting, Coates and Hall (2005) suggested that stands, where the basal area of the secondary structure is at least 5-10 m2 ha-1, are likely to provide an opportunity for mid-term (in 15-50 years) harvest, alleviating concerns of expected timber supply shortages following beetle outbreaks. This relates to the importance of the understory and secondary structure in stand recovery shown by Emmel et al. (2014) and Bowler et al (2012) in stands in the BC Interior, as well as by Hawkins et al. (2013) showing increased radial growth of the understory due to higher light levels following attack in central BC. In Colorado, Collins et al. (2011) also found increases in natural seedling establishment and faster understory growth in attacked stands. Moore et al. (2013), similar to our 154  results, found little change in R over the ten years following attack. Thus, no indications of a respiration pulse (i.e., a sudden marked increase in R due to decomposing fallen trees) were discovered in the first decade following beetle attack. Future plans should include running 3-PG for other MPB-attacked stands to determine how well the model performs at those sites compared to MPB-06. Important future work includes the verification of the contributions of Ra and Rh to R. Longer term measurements of the two components of R, as well as more measurements of the decomposition rates of snags, fallen trees and CWD would likely improve modelling results. The Rh sub-model could also be improved by including variables other than temperature, such as  and mortality-induced changes in root exudates, in controlling Rh. It would also be beneficial to continue EC measurements in order to have longer-term observations following attack, instead of only for the first decade after attack. The need to monitor C exchange for at least 2-3 decades was shown by Baldocchi et al. (2018). Continued observations at MPB-09 or other partial-harvested stands would be useful to obtain more information on the effects of this management strategy, as the three years of observations show a recovering trend, but the stand remained a C source. It would be valuable to observe how long it takes for stands like MPB-09 to become C sinks.   155  Bibliography Aber J.D. and Federer C.A., 1992. A generalized, lumped-parameter model of photosynthesis, evapotranspiration and net primary production in temperate and boreal forest ecosystems. Oecologia. 92, 463-474. Almeida A.C. and Sands P.J., 2016. Improving the ability of 3‐PG to model the water balance of forest plantations in contrasting environments. Ecohydrology. 9, 610-630. Amiro B.D., Barr A.G., Black T.A., Iwashita H., Kljun N., McCaughey J.H., Morgenstern K., Murayama S., Nesic Z., Orchansky A.L., Saigusa N., 2006. Carbon, energy and water fluxes at mature and disturbed forest sites, Saskatchewan, Canada. Agric. For. Meteorol. 136, 237-251. Arain M.A., Black T.A., Barr A.G., Jarvis P.G., Massheder J.M., Verseghy D.L., Nesic Z., 2002. Effects of seasonal and interannual climate variability on net ecosystem productivity of boreal deciduous and conifer forests. Canadian Journal of Forest Research. 32, 878-891. Arora V.K., Peng Y., Kurz W.A., Fyfe J.C., Hawkins B., Werner A.T., 2016. Potential near-future carbon uptake overcomes losses from a large insect outbreak in British Columbia, Canada. Geophys. Res. Lett. 43, 2590-2598. Baldocchi D., Chu H., Reichstein M., 2018. Inter-annual variability of net and gross ecosystem carbon fluxes: A review. Agricultural and Forest Meteorology. 249, 520-533. Barr A.G., Black T.A., Hogg E.H., Kljun N., Morgenstern K., Nesic Z., 2004. Inter-annual variability in the leaf area index of a boreal aspen-hazelnut forest in relation to net ecosystem production. Agric. For. Meteorol. 126, 237-255. BC Ministry of Forests, Lands and Natural Resource Operations, 2017. Biology of the Mountain Pine Beetle. 2p. https://www2.gov.bc.ca/assets/gov/farming-natural-resources-and-industry/forestry/forest-health/bark-beetles/biology_of_the_mountain_pine_beetle.pdf. BC Ministry of Forests, Lands and Natural Resource Operations, 2016. Provincial-Level Projection of the Current Mountain Pine Beetle Outbreak (Year 13). https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-health/forest-pests/bark-beetles/mountain-pine-beetle/mpb-projections (Accessed 6 June 2017). BC Ministry of Forests, Lands and Natural Resource Operations, 2015. Provincial-Level Projection of the Current Mountain Pine Beetle Outbreak (Year 12). https://www.for.gov.bc.ca/hre/bcmpb/year12.htm (Accessed 14 February 2016). Bellassen V. and Luyssaert S., 2014. Carbon sequestration: managing forests in uncertain times. Nature. 506, 153-155. Biederman J.A., Somor A.J., Harpold A.A., Gutmann E.D., Breshears D.D., Troch P.A., Gochis D.J., Scott R.L., Meddens A.J.H., Brooks P.D., 2015. Recent tree die-off has little effect on 156  streamflow in contrast to expected increases from historical studies. Water Resources Research. 51, 9775-9789. Borkhuu B., Peckham S.D., Ewers B.E., Norton U., Pendall E., 2015. Does soil respiration decline following bark beetle induced forest mortality? Evidence from a lodgepole pine forest. Agric. For. Meteorol. 214-215, 201-207. Bosc A., 2000. EMILION, a tree functional-structural model: Presentation and first application to the analysis of branch carbon balance. Ann. For. Sci. 57, 555-569. Bowler R., Fredeen A.L., Brown M., Andrew Black T., 2012. Residual vegetation importance to net CO2 uptake in pine-dominated stands following mountain pine beetle attack in British Columbia, Canada. For. Ecol. Manage. 269, 82-91. Brown M., Black T.A., Nesic Z., Foord V.N., Spittlehouse D.L., Fredeen A.L., Grant N.J., Burton P.J., Trofymow J.A., 2010. Impact of mountain pine beetle on the net ecosystem production of lodgepole pine stands in British Columbia. Agricultural and Forest Meteorology. 150, 254-264. Brown M. 2011. The carbon, water and energy balances of two lodgepole pine stands recovering from mountain pine beetle attack in British Columbia. Dissertation. The University of British Columbia, Vancouver, BC, Canada. Brown M.G., Black T.A., Nesic Z., Fredeen A.L., Foord V.N., Spittlehouse D.L., Bowler R., Burton P.J., Trofymow J.A., Grant N.J., Lessard D., 2012. The carbon balance of two lodgepole pine stands recovering from mountain pine beetle attack in British Columbia. Agricultural and Forest Meteorology. 153, 82-93. Brown M.G., Black T.A., Nesic Z., Foord V.N., Spittlehouse D.L., Fredeen A.L., Bowler R., Grant N.J., Burton P.J., Trofymow J.A., Lessard D., Meyer G., 2014. Evapotranspiration and canopy characteristics of two lodgepole pine stands following mountain pine beetle attack. Hydrological Processes. 28, 3326-3340. Burba G.G., McDermitt D.K., Grelle A., Anderson D.J., Xu L., 2008. Addressing the influence of instrument surface heat exchange on the measurements of CO2 flux from open‐path gas analyzers. Global Change Biology. 14, 1854-1876. Buttle J.M., 2009. Using a Temperature-Based Model of Snow Accumulation and Melt to Assess the Long-Term Hydrologic Behaviour of Forested Headwater Basins in South-Central Ontario. In: 66th Eastern Snow Conference, Niagara-on-the-Lake, Ontario, Canada. Caldwell M.K., Hawbaker T.J., Briggs J.S., Cigan P.W., Stitt S., 2013. Simulated impacts of mountain pine beetle and wildfire disturbances on forest vegetation composition and carbon stocks in the Southern Rocky Mountains. Biogeosciences. 10, 8203-8222. Campioli M., Vicca S., Luyssaert S., Bilcke J., Ceschia E., Chapin F.S.I., Ciais P., Fernández-Martínez M., Malhi Y., Obersteiner M., Olefeldt D., Papale D., Piao S.L., Peñuelas J., Sullivan P.F., Wang X., Zenone T., Janssens I.A., 2015. Biomass production efficiency 157  controlled by management in temperate and boreal ecosystems. Nature Geoscience. 8, 843-846. Canadell J.G., Le Quéré C., Raupach M.R., Field C.B., Buitenhuis E.T., Ciais P., Conway T.J., Gillett N.P., Houghton R.A., Marland G., 2007. Contributions to accelerating atmospheric CO₂ growth from economic activity, carbon intensity, and efficiency of natural sinks. Proc. Natl. Acad. Sci. U. S. A. 104, 18866-18870. Chen J.M., Liu J., Cihlar J., Goulden M.L., 1999. Daily canopy photosynthesis model through temporal and spatial scaling for remote sensing applications. Ecol. Model. 124, 99-119. Chertov O.G. and Komarov A.S., 1997. SOMM: A model of soil organic matter dynamics. Ecol. Model. 94, 177-189. Chertov O.G., Komarov A.S., Tsiplianovsky A.M., 1999. A combined simulation model of Scots pine, Norway spruce and Silver birch ecosystems in the European boreal zone. For. Ecol. Manage. 116, 189-206. Chertov O.G., Komarov A.S., Nadporozhskaya M., Bykhovets S.S., Zudin S.L., 2001. ROMUL - a model of forest soil organic matter dynamics as a substantial tool for forest ecosystem modeling. Ecol. Model. 138, 289-308. Coates K.D., DeLong C., Burton P.J., Sachs D.L., 2006. Abundance of Secondary Structure in Lodgepole Pine Stands Affected by the Mountain Pine Beetle. Report for the Chief Forester. Bulkey Valley Centre for Natural Resources Research and Management. Coates K.D. and Hall E., 2005. Implications of Alternate Silvicultural Strategies in Mountain Pine Beetle damaged stands. Technical Report for Forest Science Program Project Y051161. Bulkley Valley Centre for Natural Resources, Research and Management, Smithers, BC. Collins B.J., Rhoades C.C., Hubbard R.M., Battaglia M.A., 2011. Tree regeneration and future stand development after bark beetle infestation and harvesting in Colorado lodgepole pine stands. Forest Ecology and Management. 261, 2168-2175. Cooke B.J. and Carroll A.L., 2017. Predicting the risk of mountain pine beetle spread to eastern pine forests: Considering uncertainty in uncertain times. For. Ecol. Manage. 396, 11-25. Coops N.C., Waring R.H., Landsberg J.J., 1998. Assessing forest productivity in Australia and New Zealand using a physiologically-based model driven with averaged monthly weather data and satellite-derived estimates of canopy photosynthetic capacity. Forest Ecology and Management. 104, 113-127. Coops N.C. and Waring R.H., 2011. A process-based approach to estimate lodgepole pine (Pinus contorta Dougl.) distribution in the Pacific Northwest under climate change. Clim. Change. 105, 313-328. Coops N.C. and Wulder M.A., 2010. Estimating the reduction in gross primary production due to mountain pine beetle infestation using satellite observations. International Journal of Remote Sensing. 31, 2129-2138. 158  Coops N.C., Hember R.A., Waring R.H., 2010. Assessing the impact of current and projected climates on Douglas-Fir productivity in British Columbia, Canada, using a process-based model (3-PG). Canadian Journal of Forest Research. 40, 511-524. Coops N.C., Waring R.H., Law B.E., 2005. Assessing the past and future distribution and productivity of ponderosa pine in the Pacific Northwest using a process model, 3-PG. Ecological Modelling. 183, 107-124. Cullingham C.I., Cooke J.E.K., Dang S., Davis C.S., Cooke B.J., Coltman D.W., 2011. Mountain pine beetle host‐range expansion threatens the boreal forest. Mol. Ecol. 20, 2157-2171. Davidson E.A., Belk E., Boone R.D., 1998. Soil water content and temperature as independent or confounded factors controlling soil respiration in a temperate mixed hardwood forest. Global Change Biol. 4, 217-227. Davidson E.A., Samanta S., Caramori S.S., Savage K., 2012. The Dual Arrhenius and Michaelis–Menten kinetics model for decomposition of soil organic matter at hourly to seasonal time scales. Global Change Biol. 18, 371-384. Drew T.J. and Flewelling J.W., 1977. Some recent Japanese theories of yield-density relationships and their applications to Monterey pine plantations. Forest Science. 23, 517-534. Edburg S.L., Hicke J.A., Lawrence D.M., Thornton P.E., 2011. Simulating coupled carbon and nitrogen dynamics following mountain pine beetle outbreaks in the western United States. Journal of Geophysical Research. 116, 1-15. Emmel C., Paul-Limoges E., Black T., Christen A., 2013. Vertical Distribution of Radiation and Energy Balance Partitioning Within and Above a Lodgepole Pine Stand Recovering from a Recent Insect Attack. Boundary-Layer Meteorol. 149, 133-163. Emmel C., Paul-Limoges E., Bowler R., Black T., Christen A., 2014. Vertical distribution of carbon dioxide sources and sinks in a recovering mountain pine beetle-attacked lodgepole pine stand. Agricultural and Forest Meteorology. 195, 108-122. Environment Canada, 2017. National inventory report 1990–2015: Greenhouse gas sources and sinks in Canada. Foord V.N. and Spittlehouse D.L., 2014. Understory changes and carbon fluxes in MPB-attacked stands in the Central Interior. Northern Silviculture Committee 2014 Winter Workshop, February 18-19, 2014, University of Northern BC, Prince George, BC. Frank J.M., Massman W.J., Ewers B.E., 2016. A Bayesian model to correct underestimated 3-D wind speeds from sonic anemometers increases turbulent components of the surface energy balance. Atmospheric Measurement Techniques. 9, 5933-5953. Grant R.F., Garcia R.L., Pinter P.J., Hunsaker D., Wall G.W., Kimball B.A., LaMorte R.L., 1995. Interaction between atmospheric CO2concentration and water deficit on gas exchange and crop growth: testing of ecosys with data from the Free Air CO2Enrichment (FACE) experiment. Global Change Biol. 1, 443-454. 159  Hansen E.M., Amacher M.C., Van Miegroet H., Long J.N., Ryan M.G., 2015. Carbon dynamics in central US Rockies lodgepole pine type after mountain pine beetle outbreaks. Forest Science. 61, 665-679. Hansen E., 2014. Forest development and carbon dynamics after mountain pine beetle outbreaks. Forest Science. 60, 476-488. Hawkins C.D.B., Dhar A., Balliet N.A., 2013. Radial growth of residual overstory trees and understory saplings after mountain pine beetle attack in central British Columbia. Forest Ecology and Management. 310, 348-356. Hollinger D.Y., Kelliher F.M., Byers J.N., Hunt J.E., McSeveny T.M., Weir P.L., 1994. Carbon Dioxide Exchange between an Undisturbed Old-Growth Temperate Forest and the Atmosphere. Ecology. 75, 134-150. Horst T., Semmer S., Maclean G., 2015. Correction of a Non-orthogonal, Three-Component Sonic Anemometer for Flow Distortion by Transducer Shadowing. Boundary-Layer Meteorol. 155, 371-395. Hubbard R.M., Rhoades C.C., Elder K., Negron J., 2013. Changes in transpiration and foliage growth in lodgepole pine trees following mountain pine beetle attack and mechanical girdling. Forest Ecology and Management. 289, 312-317. Kimmins J.P. and Scoullar K., 1979. FORCYTE: a computer simulation approach to evaluating the effect of whole-tree harvesting on the nutrient budget in Northwest forests. In Forest Fertilization Conference Proceedings. Institute of Forest Resources, Univ. Washington, Seattle. Contribution No. 40. pp. 266-273. Kimmins J.P., Mailly D., Seely B., 1999. Modelling forest ecosystem net primary production: the hybrid simulation approach used in FORECAST. Ecol. Model. 122, 195-224. Krishnan P., Black T.A., Grant N.J., Barr A.G., Hogg E.H., Jassal R.S., Morgenstern K., 2006. Impact of changing soil moisture distribution on net ecosystem productivity of a boreal aspen forest during and following drought. Agricultural and Forest Meteorology. 139, 208-223. Kurz W.A., Dymond C.C., Stinson G., Rampley G.J., Neilson E.T., Carroll A.L., Ebata T., Safranyik L., 2008. Mountain pine beetle and forest carbon feedback to climate change. Nature. 452, 987-990. Kurz W.A., Dymond C.C., White T.M., Stinson G., Shaw C.H., Rampley G.J., Smyth C., Simpson B.N., Neilson E.T., Trofymow J.A., Metsaranta J., Apps M.J., 2009. CBM-CFS3: A model of carbon-dynamics in forestry and land-use change implementing IPCC standards. Ecol. Model. 220, 480-504. Landsberg JJ and Sands PJ. 2011. Physiological ecology of forest production: Principles, processes and models. 1st ed. Amsterdam; Boston: Elsevier/Academic Press. Landsberg JJ and Gower ST. 1997. Applications of physiological ecology to forest management. San Diego: Academic Press. 160  Landsberg J.J. and Waring R.H., 1997. A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning. Forest Ecology and Management. 95, 209-228. Law B.E., Waring R.H., Anthoni P.M., Aber J.D., 2000. Measurements of gross and net ecosystem productivity and water vapour exchange of a Pinus ponderosa ecosystem, and an evaluation of two generalized models. Global Change Biol. 6, 155-168. Lewis K.J. and Hrinkevich K., 2013. Post mortality rate of wood degradation and tree fall in lodgepole pine trees killed by mountain pine beetle in the Foothills and Rocky Mountains regions of Alberta. Final Report for Foothills Research Institute.https://friresearch.ca/resource/post-mortality-rate-wood-degradation-and-tree-fall-lodgepole-pine-trees-killed-mountain(Accessed 11 January 2018). Lewis K.J. and Hartley I., 2005. Rate of deterioration, degrade and fall of trees killed by mountain pine beetle: A synthesis of the literature and experiential knowledge. Mountain Pine Beetle Initiative Working Paper 2005-14, Natural Resources Canada, Canadian Forest Service, Victoria, BC. http://cfs.nrcan.gc.ca/pubwarehouse/pdfs/25483.pdf (Accessed 11 January 2018). Maness H., Kushner P.J., Fung I., 2012. Summertime climate response to mountain pine beetle disturbance in British Columbia. Nature Geoscience. 6, 65-70. Massman W.J., 1999. A model study of kBH−1 for vegetated surfaces using ‘localized near-field’ Lagrangian theory. Journal of Hydrology. 223, 27-43. Mathys A., Black T.A., Nesic Z., Nishio G., Brown M., Spittlehouse D.L., Fredeen A.L., Bowler R., Jassal R.S., Grant N.J., Burton P.J., Trofymow J.A., Meyer G., 2013. Carbon balance of a partially harvested mixed conifer forest following mountain pine beetle attack and its comparison to a clear-cut. Biogeosciences. 10, 5451-5463. McMurtrie R.E. and Landsberg J.J., 1992. Using a simulation model to evaluate the effects of water and nutrients on the growth and carbon partitioning of Pinus radiata. For. Ecol. Manage. 52, 243-260. McMurtrie R.E., Rook D.A., Kelliher F.M., 1990. Modelling the yield of Pinus radiata on a site limited by water and nitrogen. For. Ecol. Manage. 30, 381-413. Meidinger D. and Pojar J., 1991. Ecosystems of British Columbia. Research Branch, Ministry of Forests, 31 Bastion Square, Victoria, BC, 330 pp. Mellander P., Bishop K., Lundmark T., 2004. The influence of soil temperature on transpiration: a plot scale manipulation in a young Scots pine stand. Forest Ecology and Management. 195, 15-28. Meyer G., Black T.A., Jassal R.S., Nesic Z., Grant N.J., Spittlehouse D.L., Fredeen A.L., Christen A., Coops N.C., Foord V.N., Bowler R., 2017. Measurements and simulations using the 3-PG model of the water balance and water use efficiency of a lodgepole pine stand following mountain pine beetle attack. Forest Ecology and Management. 393, 89-104. 161  Mikkelson K.M., Maxwell R.M., Ferguson I., Stednick J.D., McCray J.E., Sharp J.O., 2013. Mountain pine beetle infestation impacts: modeling water and energy budgets at the hill‐slope scale. Ecohydrology. 6, 64-72. Mikkelson K.M., Dickenson E.R.V., Maxwell R.M., McCray J.E., Sharp J.O., 2012. Water-quality impacts from climate-induced forest die-off. Nature Climate Change. 3, 218. Mitchell K.J., Stone M., Grout S.E., Polsson K.R., Di Lucca M., Nigh G.D., Goudie J.W., Stone J.N., Nussbaum A.F., Yanchuk A., Stearns-Smith S., Brockley R., 2000. TIPSY version 3.0. Online. Ministry of Forests, Research Branch, Victoria, BC, Canada.https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-inventory/growth-and-yield-modelling/table-interpolation-program-for-stand-yields-tipsy (Accessed 20 January 2018). Mitchell R.G. and Preisler H.K., 1998. Fall rate of lodgepole pine killed by the mountain pine beetle in Central Oregon. West. J. Appl. For. 13, 23-26. Molnar J.L. and Kubiszewski I., 2012. Managing natural wealth: Research and implementation of ecosystem services in the United States and Canada. Ecosystem Services. 2, 45-55. Monteith J.L. and Moss C.J., 1977. Climate and the efficiency of crop production in Britain. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences. 281, 277-294. Monteith JL and Unsworth MH. 2013. Principles of environmental physics. 4th ed. Amsterdam; Boston: Elsevier/Academic Press. Moore D.J.P., Trahan N.A., Wilkes P., Quaife T., Stephens B.B., Elder K., Desai A.R., Negron J., Monson R.K., Penuelas J., 2013. Persistent reduced ecosystem respiration after insect disturbance in high elevation forests. Ecology Letters. 16, 731-737. Morgenstern K., Andrew Black T., Humphreys E.R., Griffis T.J., Drewitt G.B., Cai T., Nesic Z., Spittlehouse D.L., Livingston N.J., 2004. Sensitivity and uncertainty of the carbon balance of a Pacific Northwest Douglas-fir forest during an El Niño/La Niña cycle. Agricultural and Forest Meteorology. 123, 201-219. Mystakidis S., Davin E.L., Gruber N., Seneviratne S.I., 2016. Constraining future terrestrial carbon cycle projections using observation‐based water and carbon flux estimates. Global Change Biol. 22, 2198-2215. Nakicenovic N., Alcamo J., Davis G., de Vries B., Fenhann J., Gaffin S., Gregory K., Gruebler A., Jung T.Y., Kram T., La Rovere E.L., Michaelis L., Mori S., Morita T., Pepper W., Pitcher H., Price L., Raihi K., Roehrl A., Rogner H., Sankovski A., Schlesinger M., Shukla P., Smith S., Swart R., van Rooijen S., Victor N., Dadi Z., 2000. Special report on emission scenarios. A special report of Working Group III of the Intergovernmental Panel on Climate Change. Cambridge Univ. Press, Cambridge, UK, and New York, NY. 599 pp. Also available from www.grida.no/publications/other/ipcc_sr/?src=/climate/ipcc/emission/index.htm> (Accessed 14 February 2016). 162  Natural Resources Canada, 2017. Indicator: Forest carbon emissions and removals. http://www.nrcan.gc.ca/forests/report/disturbance/16552#new-approach (Accessed 20 January 2018). Nealis V.G. and Cooke B., 2014. Risk assessment of the threat of mountain pine beetle to Canada's boreal and eastern pine forests. Natural Resources Canada, Canadian Council of Forest Ministers, Forest Pest Working Group. NFI, 2008. Canada's National Forest Inventory Ground Sampling Guidelines.https://nfi.nfis.org/resources/groundplot/Gp_guidelines_v5.0.pdf (Accessed 13 September 2017). Nicholls D., 2017. Prince George Timber Supply Area Rationale for Allowable Annual Cut (AAC) Determination. British Columbia Ministry of Forests, Lands, Natural Resource Operations and Rural Development. Nishio G., 2010. Harvesting mountain pine beetle-killed pine while protecting the secondary structure: a comparison of partial harvesting and clearcutting methods. FPInnovations FERIC, Advantage Rep., 12. NOAA, 2018. National Oceanic and Atmospheric Administration, Earth System Research Laboratory, Global Monitoring Division, Trends in Atmospheric Carbon Dioxide. https://www.esrl.noaa.gov/gmd/ccgg/trends/full.html. Odum E.P., 1969. The strategy of ecosystem development. Science. 164, 262-270. Perttunen J., Sievänen R., Nikinmaa E., Salminen H., Saarenmaa H., Väkevä J., 1996. LIGNUM: A Tree Model Based on Simple Structural Units. Annals of Botany. 77, 87-98. Pfeifer E.M., Hicke J.A., Meddens A.J.H., 2011. Observations and modeling of aboveground tree carbon stocks and fluxes following a bark beetle outbreak in the western United States. Global Change Biology. 17, 339-350. Potithep S. and Yasuoka Y., 2011. Application of the 3-PG Model for Gross Primary Productivity Estimation in Deciduous Broadleaf Forests: A Study Area in Japan. Forests. 2, 590-609. Pretzsch H. 2009. Forest growth models. In: Forest dynamics, growth and yield: From measurement to model. Berlin, Heidelberg: Springer Berlin Heidelberg. 423 p. Price D.T., McKenney D.W., Joyce L.A., Siltanen R.M., Papadopol P., Lawrence K.M., 2011. High-resolution interpolation of climate scenarios for Canada derived from general circulation model simulations. Canadian Forest Service Northern Forestry Centre and US Forest Service, Rocky Mountain Research Station. Raison R.J. and Myers B.J., 1992. The Biology of Forest Growth experiment: linking water and nitrogen availability to the growth of Pinus radiata. For. Ecol. Manage. 52, 279-308. Reed D.E., Ewers B.E., Pendall E., 2014. Impact of mountain pine beetle induced mortality on forest carbon and water fluxes. Environmental Research Letters. 9, 1-12. 163  Romme W.H., Knight D.H., Yavitt J.B., 1986. Mountain pine beetle outbreaks in the Rocky Mountains: regulators of primary productivity? The American Naturalist. 127, 484-494. Safranyik L and Carroll AL. 2006. The biology and epidemiology of the mountain pine beetle in lodgepole pine forests. Safranyik L and Wilson B, editors. In: The mountain pine beetle: A synthesis of biology, management, and impacts on lodgepole pine. Victoria, BC, Canada: Natural Resources Canada, Canadian Forest Service, Pacific Forestry Centre. 66 p. Safranyik L., Shrimpton D.M., Whitney H.S., 1974. Management of lodgepole pine to reduce losses from the mountain pine beetle. Forestry Technical Report 1, Natural Resources Canada, Canadian Forest Service, Pacific Forestry Centre, Victoria, BC, Canada. Safranyik L., Carroll A.L., Régnière J., Langor D.W., Riel W.G., Shore T.L., Peter B., Cooke B.J., Nealis V.G., Taylor S.W., 2010. Potential for Range Expansion of Mountain Pine Beetle into the Boreal Forest of North America. The Canadian Entomologist. 142, 415-442. Sands P., 2004. Adaptation of 3-PG to novel species: guidelines for data collection and parameter assignment. Technical Report 141. Cooperative Research Centre for Sustainable Production Forestry. Sands P.J. and Landsberg J.J., 2002. Parameterisation of 3-PG for plantation grown Eucalyptus globulus. Forest Ecology and Management. 163, 273-292. Sands P. 2010. 3PGpjs user manual, pp. 1-27. Seip D. and Jones E., 2009. Response of woodland caribou to partial retention logging of winter ranges attacked by mountain pine beetle. Annual Report. FIA/FSP Project # Y091010. B.C. Ministry of Forests, Prince George, B.C., 28p. Tanner C.B. and Thurtell G.W., 1969. Anemoclinometer Measurements of Reynolds Stress and Heat Transport in the Atmospheric Surface Layer. Research and Development Technical Report ECOM-66-G22F, University of Wisconsin, Madison, Wisconsin. Taylor SW, Carroll AL, Alfaro RI, Safranyik L. 2006. Forest, climate and mountain pine beetle outbreak dynamics in Western Canada. In: The mountain pine beetle: A synthesis of biology, management, and impacts on lodgepole pine. Victoria, BC, Canada: Natural Resources Canada, Canadian Forest Service, Pacific Forestry Centre. 67 p. Vanderhoof M.K. and Williams C.A., 2015. Persistence of MODIS evapotranspiration impacts from mountain pine beetle outbreaks in lodgepole pine forests, south-central Rocky Mountains. Agricultural and Forest Meteorology. 200, 78-91. Wang Y.P. and Jarvis P.G., 1990. Description and validation of an array model - MAESTRO. Agric. For. Meteorol. 51, 257-280. Waring R.H. and Pitman G.B., 1985. Modifying lodgepole pine stands to change susceptibility to mountain pine beetle attack. Ecology. 66, 889-897. 164  Wei L., Marshall J.D., Link T.E., Kavanagh K.L., Du E., Pangle R.E., Gag P.J., Ubierna N., 2014. Constraining 3‐PG with a new δ13C submodel: a test using the δ13C of tree rings. Plant, Cell & Environment. 37, 82-100. Wenzel S., Cox P.M., Eyring V., Friedlingstein P., 2016. Projected land photosynthesis constrained by changes in the seasonal cycle of atmospheric CO2. Nature. 538, 499-501. Williams C.A., Gu H., MacLean R., Masek J.G., Collatz G.J., 2016. Disturbance and the carbon balance of US forests: A quantitative review of impacts from harvests, fires, insects, and droughts. Global Planet. Change. 143, 66-80. Winkler R., Spittlehouse D., Boon S., Zimonick B., 2015. Forest disturbance effects on snow and water yield in interior British Columbia. Hydrology Research. 46, 521. 165  Appendices Appendix A   Photographs of the measurement sites and eddy-covariance instrumentation This appendix shows photographs of the MPB-attacked stands in the BC Interior, where EC-measurements were made. These stands included two not-harvested stands, which were dominated by lodgepole pine and attacked in 2006 (MPB-06) and 2003 (MPB-03), respectively, a mixed-conifer partial-salvage-harvested stand (MPB-09), and three clearcut-harvested stands (CC-05, CC-97 and MPB-09C) near MPB-06 and MPB-09. The photographs were taken by members of the Biometeorology and Soil Physics Group at UBC.    Figure A.1 Photograph of the eddy-covariance (EC) system at the 26-m height at MPB-06 in August 2007, showing a three-dimensional ultrasonic anemometer (model CSAT3, Campbell CSAT3 Thermocouple LI-7500 166  Scientific Inc. (CSI), Logan, Utah), an open-path infrared gas analyser (IRGA) (model LI-7500, LI-COR, Inc., Lincoln, Nebraska) and a fine-wire thermocouple.   Figure A.2 Photograph of the 32-m tall measurement tower at MPB-09 in July 2011 with the eddy-covariance (EC) and net radiometer booms, three solar panels and the rain gauge (model 52203 R.M. Young Co., Traverse City, MI at MPB-09). Net radiometer Rain gauge EC system Solar panel 167      Figure A.3 Canopy photographs taken from the top of the flux tower at MPB-06 in a northern direction in 2006, 2007, 2010 and 2014.  2006 2007 2010 2014 168     Figure A.4 Canopy photographs taken from the top of the flux tower at MPB-03 in 2008 and 2014.    Figure A.5 Canopy photograph and view from the ground at MPB-09 in 2009.  2014 2008 169        Figure A.6 Photographs of the measurement towers at a) CC-97 and b) CC-05 in June and August 2007, respectively. a) b) 170    Figure A.7 Measurement tower at MPB-09C in July 2010.   171  Appendix B   Photographs of respiration instrumentation  This appendix shows photographs of the portable chamber system and collar locations used to make observations of soil, live and dead-tree respiration at MPB-06 on August 26-28, 2014. The photographs were taken by T. Andrew Black.    Figure B.1 Manual portable chamber system with a closed-path infrared gas analyser (LI-800, LI-COR), using an opaque chamber measuring respiration on a log, and a temperature probe.  172    Figure B.2 Collars on two dead trees on which bole respiration measurements were made.    Figure B.3 Collars on two live trees on which bole respiration measurements were made.  173    Figure B.4 Collars on two logs on which bole respiration measurements were made.   174         Figure B.5 Collars in locations where soil respiration measurements were made. 175  Appendix C   Root zone soil water storage at MPB-06, MPB-03 and MPB-09 This appendix shows the root zone soil water storage (S) at the beginning of the growing season (May 1) for the observed years at MPB-06, MPB-03 and MPB-09, respectively. These values correspond to no change in S (ΔS = 0) during the growing season shown in Figure 4.5.  Table C.1 Root zone soil water storage (S) at MPB-06, MPB-03 and MPB-09 at the beginning of the growing season (May 1) of the decade following attack. Year MPB-06 MPB-03 MPB-09 2007 45.22 mm -4.44 mm  2008 44.71 mm 7.65 mm  2009 42.43 mm 15.02 mm  2010 20.70 mm 2.47 mm 46.80 mm 2011 77.15 mm 14.57 mm 15.18 mm 2012 33.08 mm 10.57 mm -10.82 mm 2013 32.97 mm   2014 26.94 mm   2015 21.28 mm   2016 18.95 mm    

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0365753/manifest

Comment

Related Items