Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Comparison of carbon and energy balances of a Douglas-fir forest from pre- to post-harvest Paul-Limoges, Eugenie 2013

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

Item Metadata

Download

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

Full Text

COMPARISON OF CARBON AND ENERGY BALANCES OF A DOUGLAS-FIR FOREST FROM PRE- TO POST-HARVEST   by Eugenie Paul-Limoges  B.Sc., The University of British Columbia, 2011  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Geography)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  September 2013  ? Eugenie Paul-Limoges, 2013  ii  Abstract Stand-replacing disturbances, such as harvesting, have a major impact on the exchanges of carbon (C) and energy between forested land and the atmosphere. The former forest CO2 sinks become net CO2 sources due to the continued respiratory losses and to the significantly reduced photosynthetic uptake. Chronosequence studies, where current different-aged stands are used to reconstruct the development of an older stand, have been widely used to quantify the influence of harvesting on C and energy exchanges of forested stands. Almost no replicated measurements have been made within the Fluxnet community for same-age stands within an ecozone. Chronosequence studies assume that all sites differ only in age, and have had the same history in their abiotic and biotic components; this main assumption has been shown to be invalid in several ecological studies using chronosequences and replications are needed to explain these differences.  This study used data from the well-studied Fluxnet-Canada Douglas-fir chronosequence on Vancouver Island, where the most mature site recently reached harvesting age and was commercially harvested. C and energy balances were measured using the eddy-covariance technique and other micrometeorological instruments at the recently harvested site (HDF11) for two years following the harvest. These measurements were then compared to pre-harvest measurements at the same site (DF49) and to post-harvest measurements from another previously harvested stand (HDF00) 3 km away in the chronosequence.  The results from this study showed that the net radiation decreased from pre- to post-harvest due to the increase in albedo and surface temperature. The average annual Bowen ratio increased slightly due to the reduction in evapotranspiration following the harvest. From pre- to post-harvest, the site transitioned from being a moderate sink to being a strong source of CO2. In iii  comparison, the previously harvested stand (HDF00) was a weaker source of CO2 due to lower respiratory rates and faster vegetation recovery. The results show the importance of replicated measurements to characterize the C and energy exchanges for an ecosystem-specific stand age following a stand-replacing disturbance.   iv  Preface  The footprint analysis in this thesis is part of a different project I published during my MSc (see Appendix B). I did the LiDAR and micrometeorological analysis and wrote most of the manuscript. Dr. Andreas Christen did the footprint analysis. Paul-Limoges, E., Christen, A., Coops, N.C., Black, T.A. and Trofymow, J.A., 2013. Estimation of aerodynamic roughness of a harvested Douglas-fir forest using airborne LiDAR. Remote Sensing of Environment, 136: 225-233. v  Table of Contents Abstract .......................................................................................................................................... ii?Preface ........................................................................................................................................... iv?Table of Contents ...........................................................................................................................v?List of Tables .............................................................................................................................. viii?List of Figures ............................................................................................................................... ix?List of Symbols and Acronyms ................................................................................................. xiv?Acknowledgements .................................................................................................................... xvi?Dedication .................................................................................................................................. xvii?1? Introduction ............................................................................................................................1?1.1? Previous studies on carbon and energy exchanges following harvesting ....................... 2?1.1.1? Carbon exchange ......................................................................................................... 2?1.1.2? Energy exchange ......................................................................................................... 6?1.2? Issues with current studies .............................................................................................. 8?1.2.1? Main chronosequence assumption .............................................................................. 8?1.2.2? Lack of replication ...................................................................................................... 9?1.3? Vancouver Island harvest chronosequence ................................................................... 11?1.4? Research objectives ....................................................................................................... 13?2? Methods .................................................................................................................................14?2.1? Site descriptions ............................................................................................................ 14?2.2? EC measurements.......................................................................................................... 16?2.3? Energy balance measurements ...................................................................................... 19?2.4? Chamber measurements ................................................................................................ 20?vi  2.5? Weather measurements ................................................................................................. 21?2.6? Flux analysis and gap filling ......................................................................................... 22?3? Results and Discussion .........................................................................................................26?3.1? Forest stand characteristics and regrowth at HDF11 .................................................... 26?3.2? Weather at HDF11 ........................................................................................................ 26?3.3? Footprints for HDF11 ................................................................................................... 28?3.4? Energy exchange at HDF11 .......................................................................................... 31?3.4.1? Annual radiation and energy fluxes .......................................................................... 31?3.4.2? Diurnal and seasonal radiation and energy fluxes .................................................... 36?3.4.3? Energy balance closure ............................................................................................. 38?3.5? Comparison of radiation and energy fluxes in the chronosequence ............................. 42?3.5.1? Pre- to post-harvest monthly temperature, albedo and Bowen ratio ......................... 42?3.5.2? Annual radiation and energy fluxes .......................................................................... 45?3.6? CO2 exchange at HDF11 ............................................................................................... 50?3.6.1? Annual NEP, GEP and Re from EC measurements .................................................. 50?3.6.2? Gross ecosystem photosynthesis ............................................................................... 51?3.6.3? Ecosystem respiration ............................................................................................... 53?3.6.4? Respiration from chambers and comparison with EC measurements ...................... 56?3.7? Comparison of CO2 fluxes in the chronosequence ....................................................... 60?3.7.1? Annual NEP, GEP and Re for DF49, HDF00 and HDF11 ........................................ 60?3.7.2? Comparison of DF49 and HDF11 ............................................................................. 62?3.7.3? Comparison of HDF00 and HDF11 .......................................................................... 64?3.7.4? Implications for chronosequence studies .................................................................. 65?vii  4? Conclusions ...........................................................................................................................69?References .....................................................................................................................................71?Appendices ....................................................................................................................................79?Appendix A : Site pictures ........................................................................................................ 79?A.1? Pre- to post-harvest surface cover ............................................................................. 79?A.2? Instrumentation at HDF11 ........................................................................................ 80?A.3? Year 2 in pictures ...................................................................................................... 81?Appendix B : Paper on estimation of aerodynamic roughness ................................................. 82? viii  List of Tables Table 1: CO2 flux results for chronosequences in the first 10 years following harvesting ............ 5?Table 2: Site characteristics for DF49/HDF11 and HDF00 (Humphreys et al., 2006; Krishnan et al., 2009). ...................................................................................................................................... 15?Table 3: Mean air temperature and total precipitation for Campbell River airport climate station (49? 57?N, 125?16?W, 109 m.a.s.l., station ID 1021261). ............................................................. 28?Table 4: Annual radiation and energy fluxes for the first two years following harvesting at HDF11........................................................................................................................................... 33?Table 5: Annual radiation and energy flux densities for the first two years following harvesting at HDF00........................................................................................................................................... 48?Table 6: Annual radiation and energy flux densities for the last four years before harvesting at DF49. ............................................................................................................................................ 49?Table 7: Comparison of annual radiation and energy flux densities for DF49, HDF11 and HDF00. The values shown are the mean with the minimum and maximum in parentheses. ....... 50?Table 8: Annual values of GEP, Re, NEP, E and P for the first two years after harvesting at HDF11........................................................................................................................................... 51?Table 9: Parameters for the rectangular hyperbolic relationship between GEP and Q?. .............. 53?Table 10: Parameters for the logistic equation fit between Re and Ts. .......................................... 54?Table 11: Comparison of annual GEP, Re, NEP, E and P for DF49, HDF11 and HDF00. The values for DF49 are the average of the last four years before harvesting, whereas the values for HDF11 and HDF00 are the averages of the first two years following harvesting. ....................... 62? ix  List of Figures Figure 1: GEP, Re and NEP exchanges as a function of stand age. Source: Humphreys (2004) ... 4?Figure 2: NEP as a function of the age of stand for the Douglas-fir harvest chronosequence on Vancouver Island. Source: Black et al. (2008) ............................................................................. 12?Figure 3: Location of the three Douglas-fir stands on the east coast of Vancouver Island. Source: Chen et al. (2009) .......................................................................................................................... 16?Figure 4: Climate variables at HDF11 for the first two years following harvesting. The solid line represents the first year from 5 May 2011 to 4 May 2012, and the dashed line represents the second year from 5 May 2012 to 4 May 2013. Both lines are 5-day averages. ............................ 27?Figure 5: Daytime cumulative flux footprint contours in 2011 for HDF11 overlaying the map of log10(z0, in m) estimated from LiDAR analysis. The red cross represents the micrometeorological tower location. The contour lines are in % of the cumulative probability assuming a unit source strength. ......................................................................................................................................... 29?Figure 6: Nighttime cumulative flux footprint contours in 2011 for HDF11 overlaying the map of log10(z0, in m) estimated from LiDAR analysis. The red cross represents the micrometeorological tower location. The contour lines are in % of the cumulative probability assuming a unit source strength. ......................................................................................................................................... 30?Figure 7: Monthly-averaged radiation components at HDF11 for the first two years following harvesting (May 2011 to May 2013). ........................................................................................... 31?Figure 8: Monthly albedo (calculated as the ratio of monthly total S? / S?) at HDF11 for the first two years following harvesting (May 2011 to May 2013). .......................................................... 32?Figure 9: Monthly-averaged energy flux densities terms for the first two years following harvesting at HDF11 (May 2011 to May 2013). ........................................................................... 34?x  Figure 10: Monthly-averaged Bowen ratio for the two years following harvesting at HDF11 (May 2011 to May 2013). ............................................................................................................. 35?Figure 11: Storage components for the two years following harvesting at HDF11. The solid line represents the first year of measurements and the dashed line represents the second year of measurements. ............................................................................................................................... 35?Figure 12: Average diurnal radiation components for the first two years following harvesting at HDF11 for a) summer, b) fall, c) winter and d) spring. Seasons refer to the normal definition based on equinoxes and solstices. ................................................................................................. 36?Figure 13: Average diurnal energy balance components for the first two years following harvesting at HDF11 for a) summer, b) fall, c) winter, and d) spring. Seasons refer to the normal definition based on equinoxes and solstices. ................................................................................ 37?Figure 14: Diurnal energy storage components for the first two years following harvesting at HDF11 for a) summer, b) fall, c) winter and d) spring. The solid line represents the storage in soil while the three other overlapping lines (dots, dashes and dot-dashes) represent the storage in biomass, latent heat and sensible heat, respectively. Seasons refer to the normal definition based on equinoxes and solstices. ........................................................................................................... 37?Figure 15: Diurnal energy balance closure on a 30 min resolution for the first two years following harvesting at HDF11 for a) summer, b) fall, c) winter, and d) spring. ......................... 40?Figure 16: EBC with ? ? filtering by season for the first two years following harvesting at HDF11.The dots represent half-hourly measurements and the percents shown in each panel are the slopes of the regression lines. The line is the 1: 1 line. ........................................................... 41?xi  Figure 17: EBC ? ? filtering for all points for both years. The dots represent half-hourly measurements and the percent shown is the slope of the regression line. The line is the 1: 1 line........................................................................................................................................................ 41?Figure 18: Average monthly Ta (heavy) and Ts (light) for HDF11 (solid lines) and DF49 (dashed lines). The data for HDF11 represent the average of the first two years of measurements following harvesting, and the data for DF49 represent the average of the last four years of measurements before the harvest. ................................................................................................. 42?Figure 19: Average monthly albedo for HDF11 (solid line) and DF49 (dashed line). The data for HDF11 are averages for the first two years of measurements following harvesting, and the data for DF49 are averages for the last four years of measurements before the harvest. ..................... 44?Figure 20: Average monthly Bowen ratios for HDF11 (solid line) and DF49 (dashed line). The data for HDF11 are averages of the first two years of measurements following harvesting, and the data for DF49 are averages of the last four years of measurement before harvesting. ........... 45?Figure 21: NEP, GEP and Re at HDF11. The solid line represents the first year from 5 May 2011 to 4 May 2012, and the dashed line represents the second year from 5 May 2012 to 4 May 2013. Both lines are 5-day averages. ...................................................................................................... 51?Figure 22: GEP as a function of downwelling PAR (Q?). The solid line represents year 1 and the dashed line represents year 2 rectangular hyperbolic fits to the binned values by 50 ?mol m-2 s-1 and the parameters are in  Table 9. ............................................................................................... 52?Figure 23: Ecosystem respiration versus soil temperature at the 5-cm depth. The solid line (circles) represents year 1 and the dashed line (triangles) represents year 2. The soil temperature values represent binned values by 1? C. The lines are logistic equation fits to the binned values and the parameters are shown in Table 10. ................................................................................... 54?xii  Figure 24: Respiration residuals, after normalization using soil temperature, as a function of soil water content. ................................................................................................................................ 55?Figure 25: Daytime respiration derived from the monthly intercept of the rectangular-hyperbolic relationship between NEP and PAR. ............................................................................................ 56?Figure 26: Ecosystem (Re) and soil (Rs) respiration during the two years after harvesting at HDF11. The solid line represents the respiration from the EC system and the dashed line represents the average soil respiration measured by the three automated chambers. Both are 5-day averages. Also shown are average manual soil respiration measurements (empty circles) and average manual respiration measurements for logs and stumps (filled circles). .......................... 59?Figure 27: Orthogonal linear regression between respiration derived from the average of the three automated chambers and the EC system (5-day averages). The equation of the regression line is Rs = 0.88Re-0.83 ?mol m-2 s-1. ........................................................................................... 60?Figure 28: Comparison of CO2 fluxes at DF49, HDF11 and HDF00. DF49 represents the average of the last 4 years. HDF11 and HDF00 represent the first year following harvesting. ................ 61?Figure 29: Relationship of NEP to stand age for all sites in the Douglas-fir chronosequence. .... 66?Figure 30: Photographs showing the forest cover at DF49 and the surface cover immediately following harvesting, which became HDF11. The photographs were taken from the top of the DF49 tower at 40-m height and the view is towards east-southeast. ............................................ 79?Figure 31: Post-harvest surface residue and micrometeorological tower in June 2011 looking downslope towards the east-northeast. ......................................................................................... 80?Figure 32: Eddy-covariance sensors at 4.5- m height at HDF11. Photograph shows the Gill Instruments R3 sonic anemometer, the fine wire thermocouple (right) and the air intake to the LI-COR Inc. LI-7000 infrared gas analyzer (bottom). ....................................................................... 80?xiii  Figure 33: Year 2 in pictures. These 12 photographs represent one picture per month over the second year of measurements (in order from May 2012 to April 2013). The photographs are taken from the top of the DF49 tower at 40-m height and the view is towards the northwest. .... 81? xiv  List of Symbols and Acronyms Symbol/Acronym  Units Definition BC  British Columbia C    carbon CO2    carbon dioxide CSI  Campbell Scientific Inc. E mm evapotranspiration EBC  energy balance closure EC     eddy-covariance Fc ?mol m-2 s-1 CO2 flux density G W m-2 or MJ m-2 day-1 or GJ m-2 yr-1 soil heat flux density GEP   g C m-2 time -1or  ?mol m-2 s-1 gross ecosystem photosynthesis H W m-2  or MJ m-2 day-1 or GJ m-2 yr-1 sensible heat flux density IRGA    infrared gas analyzer LAI   m2 m-2 leaf area index  L? W m-2 or MJ m-2 day-1 or GJ m-2 yr-1 downwelling longwave  L? W m-2 or MJ m-2 day-1 or GJ m-2 yr-1 upwelling longwave  NEE   ?mol m-2 s-1 net ecosystem exchange NEP   g C m-2 time -1 or  ?mol  m-2 s-1  net ecosystem productivity  P mm precipitation PAR ?mol m-2 s-1 photosynthetically active radiation Px ?mol m-2 s-1 photosynthetic capacity Q? ?mol m-2 s-1 downwelling photosynthetically active radiation Ra    g C m-2 time -1 or   ?mol m-2 s-1 autotrophic respiration xv  Symbol/Acronym  Units Definition Rd g C m-2 time -1 or   ?mol m-2 s-1 daytime respiration Re   g C m-2 time -1 or   ?mol m-2 s-1 ecosystem respiration Rh   g C m-2 time -1 or  ?mol m-2 s-1 heterotrophic respiration Rn W m-2 or MJ m-2 day-1 or  or GJ m-2 yr-1 net radiation Rs   g C m-2 time -1 or  ?mol m-2 s-1 soil respiration r1,2,3 ?mol m-2 s-1 or ?C empirical parameters sc ?mol CO2 mol-1 dry air CO2 molar mixing ratio ?S   W m-2 or MJ m-2 day-1 or  or GJ m-2 yr-1 energy storage in the air column and biomass below the EC sensors SD    standard deviation S?  W m-2 or MJ m-2 day-1  or GJ m-2 yr-1 downwelling shortwave  S? W m-2 or MJ m-2 day-1  or GJ m-2 yr-1 upwelling shortwave  Ta   ?C air temperature  Ts   ?C soil temperature  t hours time of day (PST) u m s-1 horizontal wind speed u*   m s-1 friction velocity ??,?? m s-1 threshold friction velocity v m s-1 lateral wind speed w m s-1 vertical wind speed zm m height of EC system ? mol C mol-1 photons quantum yield ?E W m-2 or MJ m-2 day-1  or GJ m-2 yr-1 latent heat flux density ?a mol dry air m?3 dry air density ? m3 m?3 volumetric soil water content xvi  Acknowledgements I am incredibly grateful for having been co-supervised by two extremely enthusiastic and knowledgeable supervisors, Dr. Andy Black and Dr. Andreas Christen. Their passion and excitement for their science have been contagious and provided me with continuous motivation to get through the challenges of my project and life during my MSc.  Enormous thanks to Zoran Nesic, Rick Ketler and Nicholas Grant for their curiosity, patience, excitement and care relating to the technical aspects of my project. I am forever thankful for everything you taught me; my project would not have been the same success without your help. I also had the pleasure to work with the best research assistant, Trevor Baker. I very much appreciated your autonomy, problem-solving skills and hard-work in the field and laboratory. Thanks to Carmen Emmel for being a supportive friend and colleague. Sincere thanks to the Natural Sciences and Engineering Research Council of Canada (NSERC) and to the Canadian Meteorological and Oceanographic Society (CMOS) for having believed in me and supported me financially during my MSc. Also, thanks to the Department of Geography for having taken me as a candidate and funded part of my MSc. I also appreciated the financial support provided by NSERC discovery grants to Dr. Andy Black and Dr. Andreas Christen, which supported my research project. Finally, special thanks to my family and friends for encouraging me through the different steps of my life, I would not be where I am today without you all. Guillaume, Jacques, Madeleine and Sebastian, your unwavering support made this thesis possible; thank you for always being there for me.   xvii  Dedication      To Guillaume, Jacques, Madeleine and Sebastian for their continuous support.       To Andy, Andreas and Zoran for being inspiring mentors.    1  1 Introduction Disturbance and succession in terrestrial ecosystems have been of interest to ecologists for more than a century (Clements, 1916). After a disturbance, the functional characteristics of ecosystems including the carbon (C) and energy exchanges vary with the structural characteristics (Odum, 1969; Sprugel, 1985). The net carbon dioxide (CO2) exchange of an ecosystem depends on the balance between the amount fixed by photosynthesis and the amount released through respiration. A stand-replacing disturbance, such as harvesting, results in the forest stand becoming a net source of CO2 due to the continued respiratory loss of CO2 and to the significant reduction of photosynthetic uptake of CO2. On a global scale, deforestation is responsible for approximately 20% of the rise in atmospheric CO2 (Denman et al., 2007; Schimel et al., 2001). An accurate quantification of C fluxes following disturbance is therefore important in order to model regional C budgets and improve management strategies toward greater C sequestration. Chronosequence studies, where current different-aged stands are used to reconstruct the development of an older stand, have been used to characterize the C and energy exchanges at different stand ages following disturbances. However, the issues with the main chronosequence assumption and the lack of replication make the use of these numbers uncertain in regional C models and in management strategies.  One of these chronosequence studies is located on Vancouver Island, BC, Canada. The mature site of this chronosequence recently reached harvesting age and was commercially logged, thereby providing a unique and unprecedented opportunity to study the influence of a disturbance at a specific site. This study tested the main assumption of the Vancouver Island chronosequence by quantifying the effect of harvesting on the C and energy balances at a 2  specific site and determining whether the results are consistent with those of the chronosequence study.   1.1 Previous studies on carbon and energy exchanges following harvesting 1.1.1 Carbon exchange Net CO2 exchange between terrestrial ecosystems and the atmosphere is measured using the eddy-covariance (EC) technique (Baldocchi, 2003) at more than 500 sites worldwide through the FLUXNET project, which started in 1998 (Baldocchi et al., 2001; Baldocchi et al., 1996). This global network of long-term tower sites has allowed for an improved understanding of C fluxes in different biomes and climates (e.g., Baldocchi, 2008; Law et al., 2002). At first, most multi-year studies focused on the influence of climate variables on the net CO2 exchange of mature forest ecosystems (e.g., Barford et al., 2001; Carrara et al., 2003; Morgenstern et al., 2004; Pilegaard et al., 2001). Efforts were then made to quantify C dynamics at different stages of succession using chronosequence studies (Amiro, 2001; Anthoni et al., 2002; Chen et al., 2002; Clark et al., 2004; Goulden et al., 2011; Goulden et al., 2006; Humphreys et al., 2006; Kolari et al., 2004; Law et al., 2001; Litvak et al., 2003; Mkhabela et al., 2009; Rannik et al., 2002; Schulze et al., 1999), where current different-aged stands are used to reconstruct the development of an older stand, thereby bringing considerable information on net CO2 exchange trajectories following disturbances. The magnitude and direction of the net CO2 exchange of an ecosystem is highly influenced by the time since the last disturbance (Sprugel, 1985). Figure 1 illustrates C exchange trajectories following a stand-replacing disturbance, as hypothesized by Humphreys (2004) for the growth of a Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirbel) Franco) stand. Gross 3  ecosystem photosynthesis (GEP), also known as gross primary productivity (GPP), represents the uptake of CO2 through photosynthesis. Stand-replacing disturbances such as harvesting result in the removal of the photosynthesizing vegetation; the low leaf area index (LAI) following a disturbance limits C sequestration (Humphreys et al., 2005), resulting in a significant reduction in CO2 being fixed by the stand. As the vegetation recovers, so does GEP. Ecosystem respiration (Re) includes the respiration from the living biomass (autotrophic) and from the decomposition of dead organic matter by microbes (heterotrophic). Autotrophic respiration (Ra) often decreases with the loss of the canopy due to the loss of respiring roots, shoots, boles and leaves. However, an exception has been found in trembling aspen (Populus tremuloides) stands that reproduce from basal shoots, and thus the roots keep living and metabolising after harvest (Amiro, 2001). Heterotrophic respiration (Rh) has usually been found to be invariant with stand age (Amiro et al., 2010; Law et al., 2003; Luyssaert et al., 2008) or greater in younger stands due in part to their warmer soil temperatures (Noormets et al., 2008), to the increased root decay, and to the decomposition of the large detrital mass after a disturbance (Harmon et al., 1990). Reduced decomposition rates due to drier surface conditions following harvesting have also been found in a Douglas-fir chronosequence (Addison et al., 2003; Trofymow, 1998).  Net ecosystem productivity (NEP) represents the balance between GEP and Re, i.e.    NEP = GEP - Re. After a disturbance, NEP is usually negative, meaning that the stand is a C source due to the heterotrophic respiration and lack of photosynthesis. Stand-replacing disturbances therefore have a major impact on net CO2 exchange, changing CO2 sinks into large sources. Many studies have now quantified C loss from a stand after a disturbance; the magnitude of this source varies based on the biome and climate (Table 1).  4   Figure 1: GEP, Re and NEP exchanges as a function of stand age. Source: Humphreys (2004)            5  Table 1: CO2 flux results for chronosequences in the first 10 years following harvesting Authors and year Forest types Type of study Harvest site results Amiro, 2001 Canadian boreal forest: jack pine forest in the Northwest territories, aspen forest in Alberta, and mixed forest in Saskatchewan -EC with paired towers for the different forest types -No stand age replication Source of 1.6 g C m-2 day-1  Rannik et al., 2002 Scots pine forest in southern Finland -EC at a 5yr-old site compared to a 38yr-old -No stand age replication Source of 2.5 to 4.0 g C m-2 day-1 Kowalski et al., 2003 Maritime pine stands in Les Landes, Southwestern France -EC measurement at a harvest site compared to a mature site -No stand age replication Source of 200 to 340 g C m-2 yr-1 Clark et al., 2004 Slash pine plantations in North Florida -EC and chambers along a chronosequence of 3 sites -No stand age replication Source of 889 to 1269 g C m-2 yr-1 Kolari et al., 2004 Scots pine forest in southern Finland -EC along a chronosequence of 4 sites -No stand age replication Source of 386 g C m-2 yr-1 Humphreys et al., 2005 Humphreys et al., 2006 Krishnan et al., 2009  Douglas-fir stands on Vancouver Island -EC and chambers along a chronosequence of 3 sites -No stand age replication  Source of 564 to 623 g C m-2 yr-1 Mkhabela et al., 2009  Canadian boreal forest: jack pine forest -EC along a chronosequence of 3 sites -No stand age replication Source of 139 g C m-2 yr-1 Zha et al., 2009 Boreal jack pine forest in Saskatchewan -EC along a chronosequence of 4 sites -No stand age replication Source of 126 to 148 g C m-2 yr-1 Goulden et al., 2011 Canadian boreal forest at various stages of succession: Jack pine forest and black spruce forest. -EC with biometry and biomass harvests along a chronosequence of 7 sites  -No stand age replication Source of 50 to 150 g C m-2  yr-1  The stand remains a CO2 source for several years after harvest, and the source strength usually attenuates over time as plants and seedlings grow and photosynthesize. Many studies have quantified that the stand is a source for 10 to 20 years after a stand-replacing disturbance, depending on the biome and climate (Amiro et al., 2010). The stand shifts from being a source of C to a sink when GEP starts to exceed respiration (Rannik et al., 2002).  6  1.1.2 Energy exchange Studies have found a 10 to 22% reduction in daily net radiation at harvested sites when compared to mature sites (Amiro, 2001; Brown, 1972; Chen et al., 2004; McCaughey, 1981; McCaughey, 1985). The higher net radiation of mature coniferous stands is usually explained by their lower albedo (Baldocchi et al., 2000; Brown, 1972; Clark et al., 2004; Gholz and Clark, 2002; Hornbeck, 1970; McCaughey, 1981; McCaughey, 1985) in part due to the effect of multiple reflections and radiation trapping within the 3-dimensional forest canopy compared to a flat harvested surface. Coniferous forests are also optically darker (Betts and Ball, 1997; Jarvis et al., 1976; Sellers et al., 1995) and the conical structure of conifer trees allows them to shed snow and intercept shortwave radiation efficiently (Stenberg et al., 1995).These attributes allow them to absorb more solar radiation and gives them a greater potential to evaporate water and heat the air and soil.  When surface temperature and vapour pressure do not limit the fluxes, mature coniferous forests have a greater ability to transfer latent and sensible heat to the atmosphere due to their greater aerodynamic roughness, resulting in greater turbulence and consequently low aerodynamic resistance (Baldocchi et al., 2000; Jarvis et al., 1976; Jarvis and Stewart, 1979; Jarvis and McNaughton, 1986; Lee et al., 2011; McCaughey, 1981). In contrast, harvested areas have an aerodynamic resistance about an order of magnitude larger (Jarvis and Stewart, 1979; Oke, 1987). The differences in aerodynamic resistances of mature forests and harvested areas result in warmer daytime surface temperatures in harvested areas, and consequently, in greater emitted longwave radiation (Lee et al., 2011; McCaughey, 1981; McCaughey, 1985). The situation is opposite at night, with cooler surface temperatures in harvested areas than in forested areas. Lee et al. (2011) observed that harvesting resulted in a local cooling effect due to the 7  higher surface albedo during daytime and cooler surface air temperatures at night. They hypothesized that surface air temperatures in forests are warmer at night under stable and stratified conditions due to the presence of trees causing turbulent mixing, thus bringing heat from aloft to the surface. During daytime, the smaller surface roughness in harvested areas causes air temperatures to rise faster than in forests. However, they observed that this roughness effect for harvested areas north of the 45? N was offset by cooling associated with albedo and Bowen ratio changes, resulting in almost identical daily maximum temperatures for forests and harvested areas. As a result of the lower daily net radiation in harvested areas, the sensible heat flux is often reduced (Amiro, 2001; Gholz and Clark, 2002). The latent heat flux in harvested areas has also been found to be lower due to the reduced leaf area resulting in decreased evapotranspiration  (Amiro, 2001; Gholz and Clark, 2002). However, the effect of harvesting on the latent heat flux has been found to be lower than expected in a slash pine (Pinus elliottii) chronosequence study due to greater evaporation from the soil (Gholz and Clark, 2002). Jassal et al. (2009) reported that after clearcut harvesting of Douglas-fir on Vancouver Island, evapotranspiration dropped to about 70% of that of a 60-year-old stand while ecosystem evapotranspiration was fully recovered when stand age was 12 years. Soil heat fluxes have been found to be greater in harvested stands due to the lack of canopy and consequent warmer soil surface temperatures (Amiro, 2001; Kowalski et al., 2003). In general, the soil heat flux is the smallest of the energy flux terms, but it becomes important where the leaf area index has been significantly reduced such as after harvesting (Gholz and Clark, 2002). The greater soil heat flux following stand-replacing disturbances often affects the 8  C budgets, as below-surface soil temperatures influence the rate of soil respiration (Gholz and Clark, 2002).  1.2 Issues with current studies 1.2.1 Main chronosequence assumption To date, all studies investigating the influence of a stand-replacing disturbance on the C and energy exchanges of a stand have used chronosequences (Amiro, 2001; Anthoni et al., 2002; Chen et al., 2002; Clark et al., 2004; Goulden et al., 2011; Goulden et al., 2006; Humphreys et al., 2006; Kolari et al., 2004; Law et al., 2001; Litvak et al., 2003; Mkhabela et al., 2009; Rannik et al., 2002; Schulze et al., 1999). Chronosequence studies are a useful tool to obtain long-term age-related data in a short period of time. The option of quantifying the effect of a stand-replacing disturbance for a specific stand is often unavailable and would require a much longer period to study ecosystem development, which is often not available with recent EC measurements (Amiro, 2001). However, at some locations, EC measurements have now been made for more than a decade (Barr et al., 2007; Dunn et al., 2007; Krishnan et al., 2009; Urbanski et al., 2007) . The issue with chronosequence studies is that they rely on the assumption that differences in C fluxes measured at sites with similar characteristics can be linked to age or disturbances, despite other differences including site location, soil texture, soil moisture dynamics and management history. It assumes that all sites differ only in stand age, and have had the same history in their abiotic and biotic components; this main assumption has been shown to be invalid for many chronosequence studies (Johnson and Miyanishi, 2008). A study by Goulden et al. (2011) underlined the large source of uncertainty in chronosequence studies due to landscape 9  heterogeneity. The study was based on a fire chronosequence, and one of the stands showed anomalously low productivity. Later remote sensing analyses showed lower enhanced vegetation indices for this stand when compared to other similarly aged stands (Goulden et al., 2006). Although in this case it was apparent that one of the sites in the sequence was different in more than just age, this type of uncertainty is present in most chronosequence studies.    However, the full dismissal of chronosequence studies would impede the improvement in our understanding of ecological processes occurring over the long term (Walker et al., 2010).  Inferences from chronosequences must therefore be validated using a different method of study; unfortunately, this is rarely done and chronosequences are assumed to hold true without any validation (Johnson and Miyanishi, 2008). Many studies simply justify the use of a chronosequence based on the justification of a similar soil, topography or climate (e.g., Humphreys et al., 2006; Zha et al., 2009) and, although these conditions are necessary, they may not be sufficient (Johnson and Miyanishi, 2008).  1.2.2 Lack of replication In addition, due to monetary and technical limitations, current studies have usually used only one site to characterise a species-specific stand age after a certain disturbance, without any replication (Amiro, 2001; Anthoni et al., 2002; Chen et al., 2002; Clark et al., 2004; Goulden et al., 2011; Goulden et al., 2006; Humphreys et al., 2006; Kolari et al., 2004; Law et al., 2001; Litvak et al., 2003; Mkhabela et al., 2009; Rannik et al., 2002; Schulze et al., 1999). A particular site might differ from other same age stands within a certain forest type and replications are needed to account for that. Site-to-site variations can be expected due to the influence and interaction of many factors including climate and site microclimate; soil carbon, nutrients and 10  water dynamics; historical and current management practices; belowground and aboveground respiration processes; and type and rate of the re-colonizing vegetation (Amiro et al., 2006; Humphreys et al., 2005).  Multiple trajectories are therefore possible following stand-replacing disturbances, and the addition of same age sites within a certain forest type would help account for the possible spatial and temporal variations (Amiro et al., 2010). Significant spatial variability in ecosystem C fluxes has been found for sites of less than 20 years of age in the same geographical area in Wisconsin (Amiro et al., 2010), where intersite variability has been found to be of the same magnitude as interannual variability (Noormets et al., 2009). For stands with rapid recovery, interannual and spatial variability between similar sites has been found to be lower in the early than in the later years after a disturbance, but this depends on the length of the disturbance recovery (Amiro et al., 2010). At sites where disturbance recovery is slower, such as the Vancouver Island harvest sites (Humphreys et al., 2006), the first 20 years are important and need to be accurately accounted for to get the net rotational C exchange (Amiro et al., 2010). Krishnan et al. (2009) also found that intersite variability for the Vancouver Island harvest sites was greater than the interannual variability, thus showing the importance of stand age in determining the C exchange.  Evaluating variations in C and energy exchanges over the landscape is therefore difficult and many processes, such as heterotrophic respiration, are likely to differ among the sites because of the different variables influencing those processes (Amiro et al., 2006). Heterotrophic respiration is the main component of the C balance after a stand-replacing disturbance and its response to climate variations has been found to be greater in the first 20 years following a disturbance (Noormets et al., 2007; Noormets et al., 2009). Heterotrophic respiration is highly 11  influenced by the nature of the site and disturbance, and thus replicated post-disturbance chronosequences are required to accurately quantify the pattern of heterotrophic respiration and the associated NEP over time (Harmon et al., 2011).  1.3 Vancouver Island harvest chronosequence Canada?s Pacific coastal temperate rainforest is a very productive ecosystem, holding the highest amount of biomass per unit area in North America (Turner et al., 1995). Smithwick et al. (2002) showed that these forests could store more C than they currently do, if management practices were improved. Continuous measurements of CO2, water vapour and energy exchanges of this ecosystem have been made for more than a decade through a harvest chronosequence of three Douglas-fir stands: DF49 planted in 1949, HDF88 planted in 1988, and HDF00 planted in 2000 (Humphreys et al., 2006; Krishnan et al., 2009). Previous studies have found the near-end-of-rotation stand, DF49, to be a moderate C sink and the harvested stand, HDF00, to be a large C source (Figure 2) (Black et al., 2008; Humphreys, 2004; Humphreys et al., 2006; Krishnan et al., 2009; Morgenstern et al., 2004).  12   Figure 2: NEP as a function of the age of stand for the Douglas-fir harvest chronosequence on Vancouver Island. Source: Black et al. (2008)  The three sites in the chronosequence were located within 50 km of each other to minimize the confounding effect of weather (Humphreys, 2004). However, other site characteristics including ?elevation, soil nutrients, soil texture and available water content, and harvesting and management practices varied slightly between stands? (Humphreys, 2004). These characteristics have been shown to have an important influence on C fluxes following disturbances (e.g., Thornton et al., 2002). Humphreys (2004) also mentioned that the C fluxes were measured in unreplicated aged-stands and, consequently, these results were not meant to definitely state the effects of harvesting on the C dynamics of Douglas-fir stands and should be used carefully in regional C budgets. She suggested replicated experiments in recently harvested sites to fully understand the C dynamics (Humphreys, 2004). In January 2011, DF49 reached harvesting age and a large section (77 ha) of it was commercially clearcut harvested. The clearcut was replanted in April 2011 and the site was renamed from DF49 to HDF11. The harvesting of this well-studied stand provides a unique and 13  unprecedented opportunity to quantify the effect of harvesting at a specific site, and determine whether the results are consistent with those of the chronosequence study.   1.4 Research objectives The overarching goal of this thesis is to test the underlying assumption of the Vancouver Island chronosequence study and evaluate its validity when studying C and energy exchanges following harvesting. In order to fulfill this goal, this thesis has three specific objectives: 1) To quantify the C and energy balances for two years following harvesting at a recently harvested site. 2) To compare the C and energy balances from pre- to post-harvest at a specific site.  3) To compare the C and energy balances in two clearcut harvested stands 3 km apart, and assess whether they are similar.  14  2 Methods 2.1 Site descriptions The Douglas-fir harvest chronosequence is located near Campbell River on the east coast of Vancouver Island, BC, Canada (Figure 3). This chronosequence is located in the dry maritime Coastal Western Hemlock biogeoclimatic subzone, with an average annual precipitation of 1,500 mm and mean annual temperature of 9.1?C (Pojar et al., 1991). This subzone experiences a maritime climate with generally cool summers and warm winters (Pojar et al., 1991).  Prior to harvesting in 2011, DF49 (49?52?N, 125?20?W, 300 m.a.s.l.; Table 2) was a dense (1100 stems ha-1) coniferous stand of 62 years-of-age with tree heights varying between 30 and 35 m (Hilker et al., 2010). The leaf area index was 7.3 m2 m-2 (Chen et al., 2006), and the mean tree diameter at breast height was 29 cm (Morgenstern et al., 2004). This 130-ha area was planted in 1949 after the eastern half of the original old-growth stand was logged and slashed-burned in 1937, and the remainder logged and slashed-burned in 1943 (Humphreys et al., 2006). DF49 comprised 80% Douglas-fir, 17% western redcedar (Thuja plicata Donn), and 3% western hemlock (Tsuga heterophylla (Raf.) Sarg.) with only a sparse understory (Humphreys et al., 2006; Morgenstern et al., 2004). Continuous measurements of CO2, water vapour and energy exchanges at this site started in 1998, as part of a project funded by the federal-provincial Forest Renewal and Development Agreement. It became part of the Fluxnet-Canada Research Network in 2002. Detailed descriptions of EC and climate measurements at DF49 can be found in Humphreys et al. (2003) and Morgenstern et al. (2004). DF49 was fertilized with 200 kg N (urea) ha-1 in January 2007. A large section (77 ha) of DF49 was harvested from January to March 2011, and replanted with 1-year-old seedlings in April 2011 with 97% Douglas-fir and 15  3% Sitka spruce (Picea sitchensis). Continuous measurements at HDF11 started in April 2011 as part of this thesis. This thesis also makes use of previous measurements made at HDF00 (49?52?N, 125?17?W, 175 m.a.s.l.; Table 2), located 3 km ESE of DF49. HDF00 was clearcut harvested in the winter of 1999/2000 with roadside debris piles burnt, and replanted in 2000 with 93% Douglas-fir and 7% western redcedar (Humphreys et al., 2006). Prior to harvesting, this 32 ha area was a second growth Douglas-fir stand established in 1940 (Humphreys et al., 2006). Continuous measurements at HDF00 started in September 2000. Detailed descriptions of EC and climate measurements can be found in Humphreys et al. (2005) and Humphreys et al. (2006).  Table 2: Site characteristics for DF49/HDF11 and HDF00 (Humphreys et al., 2006; Krishnan et al., 2009).  DF49/HDF11 HDF00 Location 49?52?N, 125?20?W 49?52?N, 125?17?W Elevation (m.a.s.l.) 300 175 Slope 5-10? facing ENE 0-2?  Soil type Humoferric podzol Humoferric podzol Soil mineral fraction (0 to 1 m depth) Texture FC?  and WP?  (m3 m-3) Available water content (mm)  Gravelly loamy sand 0.21, 0.06 150  Gravelly loamy sand to sand0.19, 0.05 140  Soil mineral fraction (0 to 15 cm) C (mg g-1 dry soil) N (mg g-1 dry soil)   18-74 0.4-4.4   32-105 1.0-5.1  Surface organic horizons C (mg g-1 dry soil) N (mg g-1 dry soil) Average thickness (cm)   244-257 6.8-16.5 0.9-19.3   401-496 4.0-16.3 0.6-17.2  16   Figure 3: Location of the three Douglas-fir stands on the east coast of Vancouver Island. Source: Chen et al. (2009)  2.2 EC measurements Continuous turbulent fluxes of CO2, water vapour and sensible heat were measured using the EC technique (Baldocchi, 2003) at a height of 4.5 m above the soil surface. The lower height of the measurements following harvesting was to ensure that the vertical flux from the site was being measured and that the footprint did not include any substantial part of the forested area (Rannik et al., 2012). The EC instrumentation consisted of a three-dimensional sonic anemometer-thermometer (R3, Gill Instruments, Lymington, UK) and a closed-path, temperature-controlled infrared gas analyzer (IRGA) (model LI-7000, LI-COR Inc., Lincoln, NE, USA). These two instruments were chosen because they perform better under wet conditions than the CSAT3 three-dimensional sonic anemometer-thermometer (model CSAT3, Campbell Scientific Inc. (CSI), Logan, UT, USA) and open-path IRGA (model LI-7500, LI-COR Inc., 17  Lincoln, NE, USA) (Humphreys et al., 2005). The two chosen instruments were also used before harvesting (Humphreys et al., 2006) and, therefore, allowed for an accurate comparison. The LI-7000 was housed within an insulated box and air was sampled through a 3-m long insulated tube and filter at an unregulated flow of 6-9 L min-1, providing turbulent flow in the tube. The LI-7000 was run in absolute mode by flowing scrubbed nitrogen gas through the reference cell at 90 cm3 min-1. The sampling tube inlet for the closed-path IRGA was located 20 cm beneath the centre of the sonic transducer array in order to minimize sensor separation errors (Kristensen et al., 1997). The sampling tube and filter were replaced frequently to ensure a fast time response. An external stainless steel filter was added to the sample tube for one month over the 2011 winter to see if it helped improve energy balance closure; however, the stainless steel filter increased the response time for H2O fluxes and decreased the IRGA sample-cell pressure, so it was removed after the one-month period. Delays were calculated from the covariance maximization at high flux and applied to the vertical velocity time series.  EC measurements were made at a frequency of 20 Hz; using such a high frequency for measurements was particularly important on a short tower above a clearcut because near the ground there are more high frequency fluctuations in vertical wind velocity and scalar concentrations than further away from the ground (Foken et al., 2012). The sonic anemometer-thermometer and IRGA were connected to a computer using an RS232 serial connection, and half-hourly means and covariances were transmitted daily by cell phone to UBC for quality checks. High-frequency data were retrieved from the site every 3 to 6 weeks. The EC system was calibrated once a day using nitrogen gas free of CO2 and water vapour, and a span gas with a known CO2-in-air concentration. Half-hourly fluxes of CO2 (Fc, ?mol m-2 s-1) were calculated from the covariance of the CO2 mole mixing ratio and vertical velocity (Webb et al., 1980) after 18  it had undergone a three-axis coordinate rotation (Tanner and Thurtell, 1969); this three-axis coordinate rotation was particularly important at HDF11 because the site is located on a slight slope. Similarly, half-hourly fluxes of latent heat and sensible heat were calculated from the covariance of the vertical velocity with the water vapour mole mixing ratio (Webb et al., 1980) and with the air temperature derived from the sonic temperature (Kaimal and Finnigan, 1994), respectively. In addition, a fine-wire (75 ?m), home-made type E (constantan-chromel) thermocouple measured temperature within the sonic array at a frequency of 20 Hz to compare with the sonic temperature. A half-hourly averaging period was chosen because it was long enough to capture most of the large eddies (i.e., low frequency) involved in turbulent transport (Finnigan et al., 2003; Kidston et al., 2010), and stationarity was assumed to hold true for this length of time (Aubinet et al., 2000; Kidston et al., 2010). Linear detrending was not used as it likely results in some loss of the low frequency contribution to the flux.  Net ecosystem exchange, NEE (?mol m-2 s-1), was calculated as the sum of Fc and the rate of change of CO2 mixing ratio (storage) in the air column below the height of the eddy-flux measurements ( mz , m) (Humphreys et al., 2005; Morgenstern et al., 2004),  NEE /c m a cF z ds dt?? ? ,                    (1) where cs  is the half-hourly average CO2 mole mixing ratio (?mol mol dry air-1) at EC sensor height, t is the time (s), and a?  is the average molar density of dry air (mol dry air m-3). /cds dtwas calculated for each half-hour using the difference between cs for the following and previous half-hours. Positive NEE values indicate a release of CO2 to the atmosphere while negative values indicate CO2 uptake by the ecosystem. Assuming that advection, fluxes of methane and 19  other volatile organic compounds, and drainage of dissolved organic and inorganic C are negligible, NEP (?mol m-2 s-1) is equal to -NEE.  2.3 Energy balance measurements Energy balance closure (EBC) was used to assess the accuracy of EC measurements. Based on the first law of thermodynamics, the energy entering the stand volume (ground surface to canopy top) must equal the energy stored in the volume plus the energy leaving the volume. Net radiation (Rn) must therefore equal the sum of the sensible heat flux (H), latent heat flux (?E), soil-heat flux (G) and rate of energy storage change in the air column and above-ground biomass (?S), i.e., Rn = H + ?E + G + ?S. Sensible and latent heat fluxes were measured using the EC technique. Net radiation was measured using a non-ventilated CNR 1 four-way net radiometer (Kipp & Zonen B.V., Delft, The Netherlands), for which the four individual components had been calibrated against Environment Canada standards. The CNR 1 was at a height of 10 m above the harvested surface pointing south of the tower.  Soil-heat flux density was measured using eight soil-heat flux plates (four CN3, Middleton Solar, Victoria, Australia and four waterproofed Peltier coolers HP-127-1.0-1.3-71P, TE Technology Inc., Traverse City, MI, USA) to account for spatial variability in the soil following harvesting. The soil-heat flux plates were calibrated in the laboratory following the procedure outlined in Emmel et al. (2013). Soil-heat flux measurements were made at a depth of 3 cm below the soil surface and corrected for heat storage change above the measurement depth using soil temperature measurements at the 1.5- cm depth, soil moisture measurements at the 2-cm depth, and the volumetric heat capacity of organic matter, minerals and water in the soil following Blanken et al. (1997). The soil-heat flux plates were in the mineral layer while the 20  layer above comprised mainly organic matter. Starting in April 2012, temperatures of logs, stumps and debris were measured using 30-gauge, home-made insulated type E thermocouples inserted to various depths through north-south transects (north aspect below the bark, centre, south aspect below the bark) for one 1943 stump, one 2011 stump, one log and some branch debris. For each thermocouple a small (3.2-mm) hole was drilled to the required depth and sealed with silicone caulk after the thermocouple had been inserted. Surface temperatures of the soil, the 2011 stump and the log were measured using infrared thermometers (IRTS-PC, Apogee Instruments, Logan, UT, USA). The net radiometer, soil-heat flux plates, thermocouples and infrared thermometers were all connected to a CR3000; measurements were made every 3 s and average values output every 30 min. The rate of change of energy stored in debris piles, stumps and logs was then calculated (in the same way as for air column CO2 storage change) in addition to the rates of energy storage change in sensible and latent heat in the air column following Blanken et al. (1997).   2.4 Chamber measurements Three non-steady-state automated chambers, each covering an area of 0.216 m2, were used to continuously measure soil CO2 efflux (i.e., soil respiration) and assess its temporal variation. The automated chambers were fabricated in the laboratory and have been described in previous studies (Jassal et al., 2005; Jassal et al., 2007; Jassal et al., 2012). Chamber headspace concentration was measured at 1 s intervals using an IRGA (LI-840, LI-COR Inc.) for a minute and a half after lid closure once every half hour. CO2 fluxes were calculated using the average rate of change over 60 s (Jassal et al., 2012). Any photosynthesizing vegetation was physically removed from the chambers to ensure that only soil respiration was measured. The chambers 21  were located within the tower flux footprint (source area) so that the soil CO2 effluxes could be compared with the CO2 fluxes measured using the EC system. Soil respiration measurements from 30 soil collars distributed systematically over the source area footprint were also made periodically with a portable chamber system, each chamber covering 0.0079 m2 (see Jassal et al., 2005; Jassal et al., 2007), to quantify the spatial variability in soil respiration and scale the measurements from the automated chambers accordingly (Irvine et al., 2008; Vickers et al., 2012). The sampling design extended 70 m upslope and 40 m downslope from the tower with transects at every 20-m increment from the tower. In addition, respiration from 10 logs and stumps was periodically measured and the fluxes weighted by the fractional area covered by logs and stumps at the site, i.e. 0.20 m2 m-2.  2.5 Weather measurements Air temperature and relative humidity were measured at a height of 5 m using a platinum resistance temperature detector and a Vaisala HUMICAP capacitive polymer chip (HMP45C, CSI). Rainfall was measured using two tipping-bucket rain gauges (TR525M, Texas Electronics Inc., Dallas, Texas, USA); one of the rain gauges had a snow adaptor (CS705, CSI) for snowfall measurements enabling it to measure precipitation year round and a sonic ranging  sensor (SR50A, CSI) measured the depth of the snow. The incoming and outgoing photosynthetically active radiation (PAR) were measured using upward- and downward-facing quantum sensors (LI-190, LI-COR Inc., Lincoln, NE, USA), respectively. The PAR sensors were mounted 10 m above the harvested surface on a boom pointing south of the tower. Soil temperature and volumetric water content were measured at depths of 5, 10, 20, 50 and 80 cm at two locations (about 15 m apart) using home-made type T (copper-constantan) thermocouples and water 22  content reflectometers (CS616, CSI). All weather sensors were connected to a data logger (CR3000, CSI); measurements were made every 3 s and output every 30 min. For an accurate comparison of the climate during the different site years, data from the nearby Campbell River airport weather station (station ID 1021261; 49? 57?N, 125?16?W, 109 m.a.s.l.) was used to compare between the different site-years used in this thesis, as the forested and harvested land covers can influence the measured weather variables.  2.6 Flux analysis and gap filling  Quality control procedures were applied to the data through a series of stages. Data quality checks were performed daily to ensure that any malfunctioning instruments were repaired or replaced quickly, thereby ensuring almost continuous data records (largest gap in EC data over the two years was 4 days). In the first stage, EC measurements went through a quality control procedure that involved despiking and removing portions of the high frequency time series associated with calibrations before computing EC statistics (Aubinet et al., 2000; Humphreys et al., 2005). Fluxes with statistics that did not fall within reasonable limits and/or occurred during instrument malfunctions were removed.  In the second stage, the best traces were selected when there were multiple instruments measuring the same variable, and derived variables such as energy storage in soil, biomass and air were calculated. Nighttime measurements made under calm conditions were removed from all CO2 analyses (Baldocchi, 2003; Falge et al., 2001). These periods were selected based on the half-hours when the measured friction velocity (??) was less than a threshold value (??,??) of 0.15 m s-1. After sorting nighttime NEE by ??, this threshold value was determined as the minimum ?? above which NEE was no longer correlated with ???(Barr et al., 2006; Goulden et 23  al., 1997; Morgenstern et al., 2004). The relationship between ??and NEE was tested for different months and seasons, and was found to be consistent over the two years. In the third stage, gap filling was performed for any gaps in the measurement time series based on the standard procedures for the Fluxnet-Canada Research Network (Barr et al., 2004), as continuous traces over the entire year were required to output annual NEP. The gap-filling procedure made use of two simple annual empirical relationships determined from measured data. One was between Re and soil temperature at 2-cm depth and the other was between GEP and downwelling PAR. NEP was assumed to equal GEP ? Re. Measured Re was estimated as ?NEP when GEP was known to be zero; that is, at night and during the cold season (periods when both air temperature and soil temperature at the 2-cm depth were below 0?C). Measured Re was then filled based on the annual logistic relationship between Re and soil temperature using the following equation,  ? ?1e 2 31 exp srR r r T? ? ?? ?? ? ,                 (2) where Ts is the soil temperature at the 2-cm depth and 1r , 2r and 3r are model fitted empirical constants. GEP was estimated as NEP + Re during daytime or zero during nighttime and periods when both the air and soil temperatures at the 2-cm depth were below 0?C. A rectangular hyperbolic relationship between GEP and downwelling PAR was then fitted to the non-zero GEP data using the equation GEP xxQ PQ P???? ? ?  ,               (3) where ?  is the quantum yield, Q? is the downwelling PAR and xP  is the photosynthetic capacity. For each relationship, parameters were first obtained for the annual relationships using 24  Equations (2) or (3). One additional parameter per relationship was then allowed to vary over time in order to account for changes in other environmental variables or phenological stage over a short period of time. The time-varying parameters were determined using a moving window of 100 consecutive half-hour measurements moving through the year in increments of 20 half hours (for more details on this procedure see Barr et al. (2004)). Finally, gaps in NEP were modeled using the difference between the modeled GEP and Re.  The uncertainty associated with annual estimates of NEP was estimated as follows. An annual error estimate associated with a 20% random error on each half-hour value of NEP was computed using a propagation of errors procedure following Morgenstern et al. (2004). The half-hourly fluxes within the ? 20% variation were re-sampled using a bootstrap Monte Carlo method and annual sums were calculated; this procedure was repeated 500 times and the 95% confidence intervals were calculated. The uncertainty associated with the derivation of the empirical relationships used for gap-filling was assessed using a resampling technique with replacement described by Humphreys et al. (2005).   In addition to estimating Re from the relationship between nighttime measurements and soil temperature at 2-cm depth, Re was estimated following two other independent methods. Previous studies have found that daytime leaf respiration is reduced compared to nighttime values (Brooks and Farquhar, 1985; Jassal et al., 2011; Villar et al., 1995) and consequently, when using nighttime estimates, daytime Re can be overestimated (Griffis et al., 2003; Janssens et al., 2001; Reichstein et al., 2005; Suyker and Verma, 2001; Wohlfahrt et al., 2005).  Janssens et al. (2001) suggested that respiration during daytime can be overestimated by as much as 15% if photoinhibition is neglected. Stoy et al. (2006) found that using methods with daytime EC data gave more accurate gap-filled Re. Lee et al. (1999) found that using daytime Re rather than 25  nighttime Re refined their correlation, as daytime measurements were less prone to measurements issues which occurred at night under stable conditions. For this reason, the intercept of the monthly light response curves between daytime NEP and PAR was used to estimate daytime respiration (Rd) i.e., dNEP xxQ P RQ P???? ?? ? .              (4) Soil respiration (Rs) measured with the automated chambers was also compared to Re obtained from the EC measurements.  LAI was measured at the site in August 2013. Eight quadrats of 0.5 m2 were destructively sampled and the leaf area in each quadrat was then measured using an LAI-3100C Area Meter (LI-COR Inc., Lincoln, NE, USA). LAI was calculated as the total area found by the area metre for one quadrat in m2 divided by 0.5 m2 and the average LAI from the eight quadrats was taken to yield one value for the site. Footprint analysis based on Kormann and Meixner (2001) was performed to assess for the probable source area contribution of the fluxes measured by the EC system, and characterize the source area footprint (Rannik et al., 2012) and its diurnal and seasonal variations (Chen et al., 2009). More details on the footprint analysis can be found in Paul-Limoges et al. (2013).    26  3 Results and Discussion 3.1 Forest stand characteristics and regrowth at HDF11 During the first year after harvesting, the only vegetation consisted of small 1-year-old Douglas-fir seedlings and sparse understory vegetation, mostly left from the previous 61-year-old stand, resulting in an LAI of about 0.2 m2 m-2. The Douglas-fir seedlings had on average a height of 40 cm and a basal diameter of 1.02 cm. The LAI increased to a maximum of 1.78 m2 m-2 in the second year. By the end of the second year, the average height of the seedlings was 43.7 cm and their basal diameter was on average 1.16 cm. The establishment and growth of pioneer species increased over the two years. With the use of transects, the more common species were found to be purple-leaved willowherb (Epilobium ciliatum), vanilla leaf (Achlys triphylla), dagger-leaved rush (Juncus ensifolius), Oregon grape (Mahonia nervosa), foamflower (Tiarella trifoliata) and various grasses. Using the logarithmic wind profile equation and independent LiDAR- based measurements, the aerodynamic roughness of the site was found to be 0.13 m in late 2011 (Paul-Limoges et al., 2013).  3.2 Weather at HDF11 The climate at the study site is typically cool and wet during winter and warm and dry during summer. For the first two years following harvesting, HDF11 experienced similar weather (Figure 4) and typical climate (Table 3). Table 3 shows the mean air temperature and total precipitation at the Campbell River airport climate station from 1999 to 2012 and Environment Canada climate normals for 1971-2000. The mean air temperature at the study site was 7.9 ?C for the first year and 8.5 ?C for the second year, both years being within one SD of the 30-year mean for the area, while the mean soil temperature (5-cm depth) was 8.1 ?C and 8.9 ?C for the first and 27  second years, respectively. The annual precipitation was 1266 mm for the first year and 1136 mm for the second year. The annual precipitation measured at the site was lower than that measured at Campbell River airport for the 2011 and 2012 calendar years (1378 and 1644 mm for 2011 and 2012, respectively). These values were within 1 SD of the 30-year mean of 1452 mm. The mean wind speed at the site was 1.55 m s-1 in the first year and 1.72 m s-1 in the second year. The dominant wind directions at the site are downslope towards the east-northeast during the night and upslope towards the west-southwest during the day. Figure 4: Climate variables at HDF11 for the first two years following harvesting. The solid line represents the first year from 5 May 2011 to 4 May 2012, and the dashed line represents the second year from 5 May 2012 to 4 May 2013. Both lines are 5-day averages.      28  Table 3: Mean air temperature and total precipitation for Campbell River airport climate station (49? 57?N, 125?16?W, 109 m.a.s.l., station ID 1021261). Year Mean air temperature (?C) Total precipitation (mm)2012 9.2 1644 2011 8.5 1378 2010 9.6 1904 2009 n/a n/a 2008 n/a n/a 2007 n/a n/a 2006 9.4 1823 2005 9.6 1422 2004 10.2 1440 2003 9.6 1497 2002 9.3 1249 2001 8.7 1295 2000 8.7 1121 1999 8.8 1986 Mean 1971-2000 8.6 1452 SD 1971-2000 0.7 267  3.3 Footprints for HDF11 Figures 5 and 6 show the results from the turbulent source areas (footprints) of the EC system overlaid on the roughness length distribution of the site (for more details on this analysis see Paul-Limoges et al. (2013)) for daytime and nighttime, respectively. Unstable conditions during daytime allowed for a more constrained footprint surrounding the tower in a circular fashion (Figure 5). The 85% contour line (enclosing 85% of the cumulative probability for a unit source) was entirely on the clearcut harvested surface, and covered a similar area upslope and downslope from the tower. However, stable conditions at night resulted in a much larger footprint, where a significant proportion (approximately 25%) of the fluxes came from the forest upslope due to katabatic flows at night (Figure 6). Figure 6 demonstrates the importance of 29  filtering based on the ??,?? for low turbulence conditions at night, as the footprint measured under nighttime stable conditions was not as representative of the harvested surface.    Figure 5: Daytime cumulative flux footprint contours in 2011 for HDF11 overlaying the map of log10(z0, in m) estimated from LiDAR analysis. The red cross represents the micrometeorological tower location. The contour lines are in % of the cumulative probability assuming a unit source strength.   30   Figure 6: Nighttime cumulative flux footprint contours in 2011 for HDF11 overlaying the map of log10(z0, in m) estimated from LiDAR analysis. The red cross represents the micrometeorological tower location. The contour lines are in % of the cumulative probability assuming a unit source strength.   31  3.4 Energy exchange at HDF11 3.4.1 Annual radiation and energy fluxes  The net radiation (Rn) was greater in the second year (1.83 GJ m-2 yr-1) after harvesting than in the first year (1.36 GJ m-2 yr-1) due to the greater amount of downwelling shortwave (S?) and longwave (L?) radiation in year 2, despite the greater amount of upwelling shortwave (S?) and longwave (L?) radiation (Table 4 and Figure 7). The annual albedo increased from year 1 to year 2, due to the greater snow cover resulting in a maximum albedo of almost 0.70 for most of January (Figure 8). The albedo of the site was consistently at about 0.13 for both growing seasons. The albedo increased over both winters due to snowfall, but this change was greater and more sustained in year 2.  Figure 7: Monthly-averaged radiation components at HDF11 for the first two years following harvesting (May 2011 to May 2013).  32   Figure 8: Monthly albedo (calculated as the ratio of monthly total S? / S?) at HDF11 for the first two years following harvesting (May 2011 to May 2013).              33  Table 4: Annual radiation and energy fluxes for the first two years following harvesting at HDF11. Energy terms First year1 Second year2 Average ? SD Annual sums (GJ m-2 yr-1)    S?  3.70 4.08 3.88 ? 0.28   S? 0.60 0.70 0.65 ? 0.07   L? 9.56 9.93 9.74 ? 0.26   L? 11.29 11.48 11.38 ? 0.14   Rn 1.36 1.83 1.59 ? 0.33   H  0.56 0.61 0.59 ? 0.04   ?E 0.56 0.75 0.65 ? 0.13   G -0.01 0.03 0.01 ? 0.03   ?S  -0.01 0.03 0.01 ? 0.03 Ratios    Albedo 0.16 0.17 0.17 ? 0.01   Bowen ratio 1.00 0.82 0.91 ? 0.13 Other climate variables    P (mm) 1266 1136 1201? 92   E (mm) 228 305 267 ? 55   Ta (?C)3 7.9 8.5 8.2 ? 0.4   Ts (?C)4 8.1 8.9 8.5 ? 0.61 5 May 2011 to 4 May 2012 2 5 May 2012 to 4 May 2013 3 Ta at a height of 5 m 4 Ts at a depth of 5 cm    The sensible heat flux density (H) and latent heat flux density (?E) increased from the first to the second year from 0.56 to 0.61 GJ m-2 yr-1 and from 0.56 to 0.75 GJ m-2 yr-1, respectively (Table 4). During the first growing season, H strongly dominated ?E due to a drier period with less precipitation, and thus lower volumetric water content from June to September over the first year (Figure 9). However, during the second growing season, H and ?E followed a similar pattern except in June 2012, when abundant precipitation led to ?E being greater than H (Figure 9). During the first year, the Bowen ratio was positive and approximately 2 during the summer, and then became negative (approximately -3) during November, December and January 34  (Figure 10). It was then fairly constant at about 1 over the second growing season until snowfall in November, which resulted in first a Bowen ratio of -4 and then a very high Bowen ratio of 70 in January (out of scale in Figure 10). Overall, the first year had an average Bowen ratio of 1.00 whereas the second year had an average of 0.82 (Table 4). The soil heat flux density (G) (Figure 9) and the rate of energy storage change (?S) (Figure 11) were positive in the summer and negative in the winter, but the trends remained relatively constant over the two years. The energy storage in soil was found to be by far the largest of the energy storage terms, and the energy storage in sensible and latent heat in the air column between ground and measurement level and energy storage in biomass was found to be negligible (Figure 11).  Figure 9: Monthly-averaged energy flux densities terms for the first two years following harvesting at HDF11 (May 2011 to May 2013).    35   Figure 10: Monthly-averaged Bowen ratio for the two years following harvesting at HDF11 (May 2011 to May 2013).  Figure 11: Storage components for the two years following harvesting at HDF11. The solid line represents the first year of measurements and the dashed line represents the second year of measurements.   36  3.4.2 Diurnal and seasonal radiation and energy fluxes Summer had the highest Rn due to the greater amount of S? from the higher solar altitude, longer daylight and clearer skies, and due to the amount of greater L? from warmer atmospheric temperatures (Figure 12). Fall and winter at HDF11 were characterized by continuous and thick cloud cover, resulting in lower Rn for these two seasons (Figure 12). On a diurnal basis, H strongly dominated ?E during daytime during the summer and spring months, whereas at night ?E slightly exceeded zero while H was negative (Figure 13). H and ?E were both very small during fall and winter (Figure 13). G was greatest in the summer due to greater soil warming (Figure 13). Energy storage in soil was important during daytime during the summer, but this energy was then lost at night (Figure 14). All other energy storage terms were found to be negligible (Figure 14).  Figure 12: Average diurnal radiation components for the first two years following harvesting at HDF11 for a) summer, b) fall, c) winter and d) spring. Seasons refer to the normal definition based on equinoxes and solstices.   37   Figure 13: Average diurnal energy balance components for the first two years following harvesting at HDF11 for a) summer, b) fall, c) winter, and d) spring. Seasons refer to the normal definition based on equinoxes and solstices.   Figure 14: Diurnal energy storage components for the first two years following harvesting at HDF11 for a) summer, b) fall, c) winter and d) spring. The solid line represents the storage in soil while the three other overlapping lines (dots, dashes and dot-dashes) represent the storage in biomass, latent heat and sensible heat, respectively. Seasons refer to the normal definition based on equinoxes and solstices.  38  3.4.3 Energy balance closure Energy balance closure (EBC) is commonly used to assess the quality of the CO2 fluxes. The failure of EC CO2 flux measurements at low ?? at night is well-documented (Black et al., 1996; Goulden et al., 1996; Wofsy et al., 1993) and, for that reason, it is common practice to filter nighttime CO2 fluxes based on a ??,??, as the CO2 flux should not be controlled by turbulence (Barford et al., 2001; Black et al., 1996; Goulden et al., 1996). Barr et al. (2006) found that EC measurements of H and ?E were equally affected by the low turbulence conditions at night, and thus showed the necessity to extend the low ?? exclusion to H and ?E. Figure 15 shows that the diurnal EBC dropped during nighttime in the summer, thus suggesting that turbulent energy fluxes were affected by the same lack of turbulence as the CO2 fluxes. In addition, the nighttime source area footprint (Figure 6) showed that the turbulent source areas were too large during nighttime. Since the purpose of the EBC in this study is to assess the quality of the CO2 flux measurements, H and ?E were also filtered using the same threshold ?? as the CO2 fluxes. After filtering using the ??,??,?EBC was 83% in the summer (81% no filtering), 61% in the fall (56% no filtering), 49% in the winter (41% no filtering) and 78% in the spring (77% no filtering) (Figure 16). EBC was found to be best in the summer due to greater turbulent mixing allowing for a better measurement of turbulent fluxes. High and constant precipitation during the fall and winter contributed to the lower EBC. The lowest EBC was linked to snow cover and it increased drastically as soon as snow melted. Issues related to EBC when there is snow could possibly be due to the energy used in snowmelt; however, the heat flux required to melt snow has been estimated to only about 2.2% (Barr et al., 2012). The presence of a snow cover can also possibly have resulted in a smoother surface and consequently, in a larger footprint. A larger footprint 39  would include a disproportional fraction from outside the clearcut and invalidate EBC, as radiative and turbulent signals would originate from different systems. Overall, the EBC for all seasons over the two years was 81% (78% no filtering) (Figure 17). As this lack of closure can be attributed to errors in the measurement of the available energy flux, as much as to errors in the turbulent energy flux densities, no correction for EBC was applied in this study. The diurnal pattern in EBC found during summer at HDF11 (Figure 15) was similar to that found by Kidston et al. (2010) for a harvested jack pine site during June. They observed lower EBC at night, large fluctuations around sunrise and sunset, and a constant EBC during daytime. They found that at sunrise and sunset, the available energy flux was positive while the turbulent energy fluxes were negative or close to zero. Rn was negative with a larger magnitude than G and the sun was above the horizon. The low solar altitude during those times resulted in the amount of energy penetrating the above-ground biomass and residue being decoupled from G, as a disproportionately greater heat flux density went into the above-ground biomass than the soil surface.   Kidston et al. (2010) also pointed out that the fact that the storage terms were approximately 90? out of phase with the other energy fluxes, i.e., the storage fluxes were positive in the morning and negative in the afternoon, suggested that the diurnal variation in EBC can be used as a qualitative assessment of the accuracy of the storage fluxes. Assuming that an underestimation of the turbulent energy fluxes was not dependent on time of day, EBC would systematically decrease throughout the day if storage terms were underestimated or systematically increase if overestimated. The fact that EBC remained constant during daytime in the summer at HDF11 is therefore an indication that the storage terms were well measured. 40   The lack of energy balance closure is a common and well-documented issue among EC studies (Barr et al., 2006; Foken, 2008; Kidston et al., 2010; Leuning et al., 2012; Twine et al., 2000; Wilson et al., 2002). Wilson et al. (2002) found a mean EBC of 80% for 50 years of data at 22 Fluxnet sites and, therefore, the values for EBC in this study were consistent with previous studies.  Figure 15: Diurnal energy balance closure on a 30 min resolution for the first two years following harvesting at HDF11 for a) summer, b) fall, c) winter, and d) spring.  41   Figure 16: EBC with ?? filtering by season for the first two years following harvesting at HDF11.The dots represent half-hourly measurements and the percents shown in each panel are the slopes of the regression lines. The lines are the 1: 1 lines.   Figure 17: EBC with ?? filtering for all points for both years. The dots represent half-hourly measurements and the percent shown is the slope of the regression line. The line is the 1: 1 line.  42  3.5 Comparison of radiation and energy fluxes in the chronosequence 3.5.1 Pre- to post-harvest monthly temperature, albedo and Bowen ratio   The average monthly air temperature (Ta) was very similar from pre- to post-harvest (Figure 18). The average monthly Ta warmed up faster in the spring at HDF11, but then peaked later in August as opposed to in July at DF49. The removal of the canopy resulted in consistently warmer average monthly soil temperatures (Ts) at 5-cm depth at HDF11, with the largest differences being seen in the winter (Figure 18).   Figure 18: Average monthly Ta (heavy) and Ts (light) for HDF11 (solid lines) and DF49 (dashed lines). The data for HDF11 represent the average of the first two years of measurements following harvesting, and the data for DF49 represent the average of the last four years of measurements before the harvest.     43  The albedo was consistently greater at HDF11 compared to DF49 (Figure 19). Highest albedo at HDF11 was found in the winter due to the greater snow cover on the open harvested surface, compared to the forest at DF49. The albedo at HDF11 was lowest in April, May and June due to saturated and dark surface organic layer leading to more S? being absorbed, and it increased with plant growth and drier site conditions over the growing season. In contrast, the albedo at DF49 stayed relatively constant throughout the year thereby demonstrating the lack of snow cover on the trees and the smaller changes in surface cover over the year. The albedo values measured at DF49 in the last four years (about 0.08) were similar to those usually reported for coniferous forests, i.e., about 0.083 (Baldocchi et al., 2000; Jarvis et al., 1976). Jassal et al. (2009) reported a dry foliage albedo of 0.085 before the fertilization, although this value would possibly have been lower if the wet foliage values would have been included. The decrease in albedo following fertilization might possibly indicate the influence of fertilization on this very productive stand, resulting in more reflections within the very tall canopy.   44   Figure 19: Average monthly albedo for HDF11 (solid line) and DF49 (dashed line). The data for HDF11 are averages for the first two years of measurements following harvesting, and the data for DF49 are averages for the last four years of measurements before the harvest.    The Bowen ratio for HDF11 and DF49 was similar for the period from February to July (Figure 20). The Bowen ratio was lower at HDF11 during January and December due to the presence of a thicker and longer lasting snow cover resulting in lower ?E from sublimation at HDF11 than from the sublimation, evaporation and transpiration from the forest canopy at DF49. The Bowen ratio at HDF11 increased in August and September, as dry surface conditions resulted in less moisture available for evapotranspiration at the site and early senescence of the broadleaf herbs and shrubs. In comparison, the forest was able to maintain a more humid environment within the forest canopy even during the drier summer period, perhaps due to roots reaching soil moisture at deeper depths.   45   Figure 20: Average monthly Bowen ratios for HDF11 (solid line) and DF49 (dashed line). The data for HDF11 are averages of the first two years of measurements following harvesting, and the data for DF49 are averages of the last four years of measurement before harvesting.   3.5.2 Annual radiation and energy fluxes From pre- to post-harvest, annual Rn decreased by 35% at HDF11 and by 25% at HDF00. The average Rn was 2.48 GJ m-2 yr-1 at DF49 (Table 6), 1.59 GJ m-2 yr-1 at HDF11 (Table 4) and 1.82 GJ m-2 yr-1 at HDF00 (Table 5). From pre- to post-harvest, the albedo increased in both clearcuts and therefore, so did S?. The average pre-harvest albedo was 0.08 at DF49 (Table 6), whereas the average post-harvest albedo was 0.17 at HDF11 (Table 4) and 0.18 at HDF00 (Table 5). Other studies have also found a reduction in Rn at harvested sites when compared to mature sites (Amiro, 2001; Brown, 1972; Chen et al., 2004; McCaughey, 1981; McCaughey, 1985), and this has often been explained by the increase in albedo following harvesting (Baldocchi et al., 2000; Brown, 1972; Clark et al., 2004; Gholz and Clark, 2002; Hornbeck, 1970; McCaughey, 46  1981; McCaughey, 1985). The removal of the DF49 canopy resulted in less incident shortwave radiation getting trapped and absorbed by the land-surface, and thus, more reflected. In addition, the smoother surface at HDF11 resulted  in a deeper and longer lasting snow cover from pre- to post-harvest at HDF11, thereby resulting in a greater shortwave reflectivity.   Annual H and ?E decreased from the forest at DF49 (0.67 GJ m-2 yr-1 and 1.03 GJ m-2 yr-1, respectively) to the harvested stand at HDF11 (0.59 GJ m-2 yr-1 and 0.65 GJ m-2 yr-1, respectively) (Table 7). Annual H and ?E were greater at HDF00 (0.82 GJ m-2 yr-1 and 0.70 GJ m-2 yr-1, respectively) than at HDF11 (Table 7). DF49 therefore had a greater evapotranspiration (E) (420 mm for 1198 mm of precipitation) compared to HDF11 (267 mm for 1201 mm of precipitation). This is consistent with previous results for HDF00 from Jassal et al. (2009) that reported that after clearcut harvesting of Douglas-fir on Vancouver Island, E dropped to about 70% of that of the 60-year-old stand. As indicated earlier, when surface temperature and vapour pressure do not limit the fluxes, mature coniferous forests have a greater ability to transfer latent and sensible heat to the atmosphere due to their greater aerodynamic roughness, resulting in greater turbulence and consequently low aerodynamic resistance (Baldocchi et al., 2000; Jarvis et al., 1976; Jarvis and Stewart, 1979; Jarvis and McNaughton, 1986; Lee et al., 2011; McCaughey, 1981). In contrast, harvested areas have an aerodynamic resistance about an order of magnitude larger (Jarvis and Stewart, 1979; Oke, 1987). Other factors influencing are the lower daily Rn resulting in less energy available for H and ?E, and the reduced leaf area resulting in decreased E (Amiro, 2001; Gholz and Clark, 2002). The partitioning of the available energy flux as H and ?E changed from pre- to post- harvest. The annual Bowen ratio was on average greater at HDF11 (0.91) and HDF00 (1.17) compared to DF49 (0.65), thereby indicating a greater proportion of the turbulent energy fluxes leaving as sensible heat following the harvest.  47  Average annual air temperatures were very similar for the three sites (DF49: 7.9?C, HDF11: 8.2?C, and HDF00: 8.3?C). Annual soil temperatures at the 5-cm depth were on average higher at HDF11 (8.5?C) and at HDF00 (9.7 ?C) compared to DF49 (7.3?C). L? was on average slightly greater at HDF11 compared to DF49 (11.38 GJ m-2 yr-1 vs. 10.83 GJ m-2 yr-1, respectively).  L? was influenced by higher soil surface temperatures at HDF11, whereas at DF49 it was to a large extent controlled by leaf temperatures, which are moderated by the cooling effect of transpiration. The greater aerodynamic roughness at DF49 resulted in the leaf temperature being coupled with the air temperature. The differences in aerodynamic resistance of mature forests and harvested stands have been shown to result in warmer daytime surface temperatures in harvested areas, and consequently, in greater L? (Lee et al., 2011; McCaughey, 1981; McCaughey, 1985). The opposite has been found at night, with cooler surface temperatures in harvested areas than in forested areas. As mentioned earlier, Lee et al. (2011) observed that harvesting resulted in a local cooling effect due to the greater surface albedo during daytime and cooler surface air temperatures at night. They hypothesized that surface air temperatures in forests are warmer at night under stable and stratified conditions due to the presence of trees causing turbulent mixing, thus bringing heat from aloft to the surface. During daytime, the smaller surface roughness in harvested areas causes air temperatures to rise faster than in forests, but this roughness effect for harvested sites north of the 45? N was found to be offset by cooling associated with albedo and Bowen ratio changes, resulting in almost identical daily maximum temperatures for forests and harvested areas. However, on an annual basis in this study, no significant difference was found in average air temperatures.   48  Table 5: Annual radiation and energy flux densities for the first two years following harvesting at HDF00. Energy terms  First year1 Second year2 Average ? SD Annual Sums (GJ m-2 yr-1)    S?  3.83 4.17 4.00 ? 0.23   S? 0.67 0.79 0.73 ? 0.08   L? n/a n/a n/a   L? n/a n/a n/a   Rn  1.78 1.86 1.82 ? 0.06   H  0.82 0.82 0.82 ? 0.00   ?E  0.72 0.68 0.70 ? 0.03   G 0.03 0.02 0.02 ? 0.01   ?S 0.03 0.02 0.03 ? 0.01 Ratios    Albedo 0.18 0.19 0.18 ? 0.01   Bowen ratio 1.14 1.20 1.17 ? 0.05 Other climate variables    P (mm) 973 1224 1099 ? 178   E (mm) 295 278 286 ? 12   Ta (?C)3 8.2 8.4 8.3 ? 0.1   Ts (?C)4 9.7 9.6 9.65 ? 0.11 1 September 2000 to 31 August 2001 2 1 September 2001 to 31 August 2002 3 Ta at a height of 3 m 4 Ts at a depth of 5 cm           49  Table 6: Annual radiation and energy flux densities for the last four years before harvesting at DF49. Energy terms  2007 2008 2009 2010 Average ? SD   Annual Sums (GJ m-2 yr-1)     S?  4.28 3.70 3.98 4.44 4.10 ? 0.33    S? 0.34 0.28 0.30 0.34 0.31 ? 0.03    L? 9.99 9.71 9.15 9.26 9.53 ? 0.39    L? 11.35 10.87 10.39 10.69 10.83 ? 0.40    Rn 2.57 2.26 2.44 2.67 2.48 ? 0.18    H 0.72 0.45 0.69 0.83 0.67 ? 0.16    ?E 1.00 1.07 1.00 1.06 1.03 ? 0.04    G -0.03 -0.03 -0.02 -0.02 -0.02 ? 0.00    ?S  -0.03 -0.03 -0.02 -0.02 -0.03 ? 0.00   Ratios     Albedo 0.08 0.08 0.07 0.08 0.08 ? 0.00    Bowen ratio 0.72 0.42 0.69 0.78 0.65 ? 0.16   Other climate variables     P (mm) 1711 1106 884 1091 1198 ? 357    E (mm) 407 437 407 433 421 ? 16    Ta (?C)1 8.4 7.7 7.3 8.1 7.9 ?0.5    Ts (?C)2 7.6 7.3 7.1 7.2 7.3 ? 0.2 1 Ta at a height of 41 m 2 Ts at a depth of 5 cm             50  Table 7: Comparison of annual radiation and energy flux densities for DF49, HDF11 and HDF00. The values shown are the mean with the minimum and maximum in parentheses.  Energy terms  DF49 HDF11 HDF00     Annual Sums (GJ m-2 yr-1)       S?  4.10 (3.70-4.44) 3.88 (3.69-4.08) 4.00 (3.83-4.17)    S? 0.31 (0.28-0.34) 0.65 (0.60-0.70) 0.73 (0.67-0.79)    L? 9.53 (9.15-9.99) 9.74 (9.56-9.93) n/a    L? 10.83 (10.39-11.35) 11.38 (11.29-11.48) n/a    Rn 2.48 (2.26-2.67) 1.60 (1.37-1.84) 1.82 (1.78-1.86)    H 0.67 (0.45-0.83) 0.59 (0.56-0.61) 0.82 (0.82-0.82)    ?E 1.03 (1.00-1.07) 0.65 (0.56-0.75) 0.70 (0.68-0.72)    G -0.02 (-0.03- -0.02) 0.01 (-0.01-0.03) 0.02 (0.02-0.03)    ?S  -0.02 (-0.03- -0.02) 0.01 (-0.01-0.03) 0.03 (0.02-0.03)     Ratios       Albedo 0.08 (0.07-0.08) 0.17 (0.16-0.17) 0.18 (0.18-0.19)    Bowen ratio 0.65 (0.42-0.78) 0.91 (0.82-1.00) 1.17 (1.14-1.20)     Other climate variables       P (mm) 1198 (884-1711) 1201 (1136-1266) 1099 (973-1224)   E (mm) 420 (407-437) 267 (228-305) 286 (276-295)    Ta (?C) 7.9 (7.3-8.4) 8.2 (7.9-8.5)  8.3 (8.2-8.4)    Ts (?C) 7.3 (7.1-7.6) 8.5 (8.1-8.9)  9.65 (9.6-9.7)  3.6 CO2 exchange at HDF11 3.6.1 Annual NEP, GEP and Re from EC measurements Figure 21 and Table 8 show the annual courses and values of net ecosystem productivity (NEP), gross ecosystem photosynthesis (GEP) and ecosystem respiration (Re) for the first two years following harvesting at HDF11. HDF11 was a strong source of carbon for both of these years with an NEP of -1000 g C m-2 yr-1 and -700 g C m-2 yr-1 for the first and second years, respectively. Over the first year, HDF11 was a much stronger source in summer than in winter, whereas it was more constant source during the second year. During the first year, GEP was almost nonexistent due to the slow recovery of vegetation (130 g C m-2 yr-1), while it increased 51  considerably (almost tripled) in the second year with vegetation recovery during the growing season (385 g C m-2 yr-1). Re was much greater than GEP during the entire year, thus accounting for most of the measured NEP. Re was slightly greater in the first year (1130 g C m-2 yr-1) than in the second year (1085 g C m-2 yr-1).  Figure 21: NEP, GEP and Re at HDF11. The solid line represents the first year from 5 May 2011 to 4 May 2012, and the dashed line represents the second year from 5 May 2012 to 4 May 2013. Both lines are 5-day averages.    Table 8: Annual values of GEP, Re, NEP, E and P for the first two years after harvesting at HDF11.   Year 1 Year 2 AverageGEP (g C m-2 yr-1) 130 385 258Re (g C m-2 yr-1) 1130 1085 1108NEP (g C m-2 yr-1) -1000 -700 -850P (mm) 1266 1136 1201E (mm) 228 305 267 3.6.2 Gross ecosystem photosynthesis Photosynthesis was almost nonexistent during the first year (Figure 22) due to the small amount of living vegetation at the site, with an LAI of approximately 0.2 m2 m-2. The few plants 52  growing at the site following harvesting were mostly pre-harvest understory species in the forest and showed clear signs of water stress over the first year during the peak of the growing season, as understory plant species are not used to warmer temperatures and higher light levels. Unlike the relationship between GEP and downwelling PAR (Q?) observed in the second year, there was virtually no response to Q? in the first year. Previous studies have found that the net photosynthetic rate of most understory species at DF49 was about 3.0 to 3.5 ?mol m-2 s-1 (Leitch, 2010). The light response curve for year 2 (Figure 22) is consistent with pre-harvest understory measurements, as GEP leveled off at about 3 ?mol m-2 s-1.   Figure 22: GEP as a function of downwelling PAR (Q?). The solid line represents year 1 and the dashed line represents year 2 rectangular hyperbolic fits to the binned values by 50 ?mol m-2 s-1 and the parameters are in  Table 9.    53  Table 9: Parameters for the rectangular hyperbolic relationship between GEP and Q?.  ?  (mol C mol-1 photons) Px (?mol m-2 s-1) R2Year 1 0.01 0.72 0.00Year 2 0.03 3.57 0.18 3.6.3 Ecosystem respiration  For the same soil temperature at the 5-cm depth, Re was slightly greater in the first year than in the second year (Figure 23). The decomposition of labile materials and fine roots over the first growing season likely contributed to the steeper slope of Re vs. soil temperature (Ts) in the first year of measurements. The logistic relationship between Re and soil temperature at HDF11 contrasts with that found at HDF00, where an exponential relationship was found to be a better fit (Humphreys et al., 2006). This suggests that other factors were limiting Re at HDF11 at higher Ts. To investigate whether soil water content was a limiting factor, residuals were obtained by subtracting the best-fit values obtained using Eq. (2) from the original values. Figure 24 shows the relationship between the residuals and volumetric soil water content at the 10-cm depth. However, as Figure 24 demonstrates, no relationship was found. This is consistent with previous results from Humphreys et al. (2006) and Drewitt et al. (2002) indicating that, despite the variations in soil water content during the growing season, no clear relationship between Re and soil water content could be established for all three stands on Vancouver Island. Drewitt et al. (2002) suggested that the poor relationships between soil water content and Rs for the Vancouver Island sites may have been due to the highly correlated nature of soil water content and temperature at these sites. Other studies have also found that the influence of soil water content on Re was hard to quantify due to the confounding effect of a negative correlation 54  between soil water content and temperature (Davidson et al., 1998). Other possible factors influencing Re relate to the interactions between the climate variables, the physical and chemical characteristics of the substrate, and the decomposer organisms, which are all influenced by the nature of the disturbance (Harmon et al., 2011).  Figure 23: Ecosystem respiration versus soil temperature at the 5-cm depth. The solid line (circles) represents year 1 and the dashed line (triangles) represents year 2. The soil temperature values represent binned values by 1? C. The lines are logistic equation fits to the binned values and the parameters are shown in Table 10.  Table 10: Parameters for the logistic equation fit between Re and Ts.  r1 (?mol m-2 s-1) r2 (?C-1) r3 (?C) R2Year 1 5.18 0.34 6.54 0.42Year 2 5.59 0.24 8.61 0.28 55   Figure 24: Respiration residuals, after normalization using soil temperature, as a function of soil water content.   Figure 25 shows the intercepts of the monthly rectangular hyperbolic relationships (see Eq. 4) between NEP and PAR. Although these points should be considered as only estimates of daytime respiration, these results are consistent with the results found from the nighttime relationship with soil temperature (Figure 21). As mentioned previously, some studies have found that daytime leaf respiration is reduced compared to nighttime values (Brooks and Farquhar, 1985; Janssens et al., 2001; Villar et al., 1995) and consequently, when using nighttime estimates, daytime Re can be overestimated. However, in our case, there is only a minimal amount of biomass photosynthesizing at the site, so photoinhibition had only a small effect. The role of photoinhibition is likely to be greater where photosynthesis is more important, such as for a forest canopy. Janssens et al. (2001) found that photoinhibition reached a reduction of up to 15% in European forests, while Suyker and Verma (2001) found a 22% reduction in a 56  tallgrass prairie. Cai (2006) found annual Re of DF49 obtained using the light response relationship was approximately 75% of that obtained using the nighttime NEE vs. Ts relationship. However, Griffis et al. (2004) found that the reduction in boreal mature aspen, black spruce and jack pine was about 10% and was within their uncertainty in the nighttime Re estimates.  Figure 25: Daytime respiration derived from the monthly intercept of the rectangular-hyperbolic relationship between NEP and PAR.  3.6.4 Respiration from chambers and comparison with EC measurements The chambers measured soil respiration (Rs) only, whereas the EC system measured ecosystem respiration (Re). Rs measurements derived from the average of the three automated chambers were well correlated but significantly less than Re especially in winter (Figure 26). Rs and Re followed similar trends annually, with much greater respiratory rates during the summer than winter.  57  Figure 26 and Figure 27 show the constant offset (about 1 ?mol m-2 s-1) between Rs and Re, representing the respiration from the decomposing stumps, logs and other debris (not included in the automated chambers), as well as autotrophic respiration. As there are possibly measurement errors in both the automated chamber and EC measurements, an orthogonal linear regression was used to assess the relationship between the two. On an annual basis, Rs was 0.88 of Re (Figure 27). Similarly, the average Rs was found to be about 0.88 of the average Re at HDF00 (Humphreys et al., 2006).  The post-harvest ratio of Rs to Re was greater than the pre-harvest ratio at DF49. Jassal et al. (2007) found that, on an annual basis, Rs accounted for 0.62 of Re at DF49. In that study, Rs :  Re was a minimum in the spring with 0.52 and gradually increased until a maximum of 0.86 in winter. The post-harvest ratio at HDF11 (0.88) was similar to the winter ratio at DF49 (0.86), thereby showing the reduced influence of the autotrophic component during the winter season. In the winter, low air temperatures results in low aboveground respiration, even in forests, but some heterotrophic respiration still occurs from the warmer and deeper soil layers. With the increase in temperature, both Rs and Re increase, but the increase in Re is greater as autotrophic respiration increase with photosynthesis. Point measurements using the manual respiration chamber are also shown in Figure 26. These measurements represent different temporal and spatial scales than those of the automated chambers; the manual chamber measurements represented one point measurement on a day (i.e., at one specific time in the daytime), whereas the automated chambers and EC system measured fluxes continuously which were then averaged. Nevertheless, these point measurements help to quantify for the spatial variability in Rs, as they cover a larger area than the automated chambers, and the similarity in the measurements show that the three automated chambers well represented 58  the spatial variability of the site. Respiratory fluxes from logs and stumps were also measured. The fluxes measured from the logs and stumps were multiplied by their estimated fractional area at the site, which was 0.20 m2 m-2, whereas soil was assumed to cover the entire site. For specific measurement periods, the ratio of respiration from logs and stumps to that of the soil was on average 7%. The average respiratory flux from the logs and stumps was 0.15 ?mol m-2 s-1 and the standard deviation was 0.6 ?mol m-2 s-1. Diurnal patterns of respiration from the chambers and the EC system were also examined following Vickers et al. (2012) in order to investigate whether advection accounted for CO2 loss from the site. They used the nocturnal ratio of NEE from EC to chamber-measured ecosystem respiration to identify the occurrence of advection and found advection occurred in stable conditions when atmospheric mixing was suppressed. In most EC studies where there is adequate fetch, horizontal and vertical flux divergences are assumed to be negligible and advection terms are often omitted in the NEE equation. Advection can lead to the horizontal transport of C and energy below the EC sensors, and thus can result in systematic measurement errors (Federer, 1970). HDF11 is located on a slight slope. Comparing fluxes measured by the chambers and the EC system permits the investigation of the effects of advection; if, for example, advective fluxes were present during nighttime, flux densities measured with the chambers would be greater (or lower) than those measured with the EC system. However, preliminary results suggested that the difference between the measurements from the automated-chamber and EC systems was not dependent on time of day, and thus suggested no significant advection. Leitch (2010) found non-negligible horizontal and vertical advection of C at DF49, and the advection correction was similar to the low friction velocity????) filtering correction, thus suggesting that the ??? criterion can successfully filter for cases when advection is relevant. 59   Figure 26: Ecosystem (Re) and soil (Rs) respiration during the two years after harvesting at HDF11. The solid line represents the respiration from the EC system and the dashed line represents the average soil respiration measured by the three automated chambers. Both are 5-day averages. Also shown are average manual soil respiration measurements (empty circles) and average manual respiration measurements for logs and stumps (filled circles). 60   Figure 27: Orthogonal linear regression between respiration derived from the average of the three automated chambers and the EC system (5-day averages). The equation of the regression line is Rs = 0.88Re-0.83 ?mol m-2 s-1.  3.7 Comparison of CO2 fluxes in the chronosequence 3.7.1 Annual NEP, GEP and Re for DF49, HDF00 and HDF11 The CO2 fluxes in Figure 28 represent the first year of measurements following harvesting at HDF11 and HDF00, and the average of the last four years of measurements at DF49. From pre- to post-harvest, the stand transitioned from being a carbon sink of 560 g C m-2 yr-1 to being a strong carbon source of 1000 g C m-2 yr-1. With the removal of the canopy, GEP significantly declined from 2020 g C m-2 yr-1 to 130 g C m-2 yr-1. Re only slightly declined from 1460 g C m-2 yr-1 to 1130 g C m-2 yr-1. In comparison to HDF11, the previously harvested stand, HDF00, was a weaker source in the first year following harvesting of 620 g C m-2 yr-1 due to a lower Re (840 g C m-2 yr-1) and a greater GEP (220 g C m-2 yr-1).  61   Figure 28: Comparison of CO2 fluxes at DF49, HDF11 and HDF00. DF49 represents the average of the last 4 years. HDF11 and HDF00 represent the first year following harvesting.  The difference in CO2 fluxes between HDF11 and HDF00 remained similar in the second year. In the second year, HDF11 was a source of 700 g C m-2 yr-1 while HDF00 was a source of 520 g C m-2 yr-1. The average NEP for the first two years following harvesting at HDF11 was      -850 g C m-2 yr-1, whereas it was -570 g C m-2 yr-1 for HDF00 (Table 11). In the second year, GEP at HDF00 increased to 530 g C m-2 yr-1, while it stayed low at HDF11 with 385 g C m-2    yr-1. The average GEP for the two years was 258 and 375 g C m-2 yr-1  for HDF11 and HDF00, respectively (Table 11). Over the second year, Re increased to 1050 g C m-2 yr-1 at HDF00, while it decreased to 1085 g C m-2 yr-1 at HDF11. The average Re for the two years following harvesting was 1108 g C m-2 yr-1 at HDF11 and 945 g C m-2 yr-1 for HDF00 (Table 11).   62  Table 11: Comparison of annual GEP, Re, NEP, E and P for DF49, HDF11 and HDF00. The values for DF49 are the averages of the last four years before harvesting, whereas the values for HDF11 and HDF00 are the averages of the first two years following harvesting. Annual total DF49 HDF11 HDF00 GEP (g C m-2 yr-1) 2020 258 375 Re (g C m-2 yr-1) 1460 1108 945 NEP (g C m-2 yr-1) 560 -850 -570 E (mm) 420 267 286 P (mm) 1198 1201 1099  The uncertainty in annual estimates of NEP results from the random and systematic errors in the measurements and their analysis. The bootstrap Monte Carlo procedure indicated that a random error of 20% for each of the half-hour measurements of NEP resulted in errors of only 2 g C m-2 yr-1 for HDF00 and HDF11, and 5 g C m-2 yr-1 for DF49. The uncertainty associated with the derivation of the empirical relationships used for gap-filling was found to be less than 25 g C m-2 yr-1 for NEP for all three sites (Humphreys et al., 2006). Previous studies have shown that given a consistent treatment of the data, differences among forests (Humphreys et al., 2006) and among years (Humphreys et al., 2005; Morgenstern et al., 2004) were similar regardless of the gap-filling strategy.  3.7.2 Comparison of DF49 and HDF11 The removal of the canopy at DF49 and subsequent low LAI resulted in a large change in GEP, whereas there was only a small change in Re. Similar results have been found by Noormets et al. (2012). However, the partitionning of Re into its components changed significantly. Re decreased from pre- to post-harvest due to the significant reduction in autotrophic respiration (Ra), as the removal of the canopy resulted in the loss of respiring roots, shoots, boles and leaves. Since Re values stayed high from pre- to post-harvest, this suggests that heterotrophic respiration 63  (Rh) stayed either constant or increased. Previous studies have shown that Rh can be invariant with age (Amiro et al., 2010; Law et al., 2003; Luyssaert et al., 2008) or greater in younger stands due in part to their warmer soil temperatures (Noormets et al., 2008), to the increased root decay, and to the decomposition of the large amount of residue after a disturbance (Harmon et al., 1990).   Jassal et al. (2007) found that pre-harvest at DF49, 0.54 of GEP was respired back to the atmosphere as Ra. This Ra fraction is similar to the 0.53 in the empirical equation by Waring and Running (1998), where Rh = Re - 0.53 GPP. Before fertilization at DF49, this equation gave somewhat reasonable results for Rh in a study by Grant et al. (2010) compared to the measurement results for Rh from Jassal et al. (2007) (758 ? 127 g C m-2 yr-1 vs. 585 ? 32 g C m-2 yr-1, respectively). However, the fertilization at DF49 increased GEP while Re remained relatively constant (Jassal et al., 2010a). Using the average values for the last four years at DF49 (Table 11), the equation by Waring and Running (1998) gives an Rh value of 389 g C m-2 yr-1 and suggests that Rh decreased with fertilization as was found in the study by Grant et al. (2010). Although this Rh value seems very low for DF49 compared to the 6% decrease found by Jassal et al. (2010b) and possibly indicates that the equation does not apply as well to fertilized stands, it strongly suggests that Rh did not decrease from pre- to post-harvest like for other Vancouver Island sites, where drier surface conditions following harvesting had been associated with reduced decomposition rates (Addison et al., 2003; Trofymow, 1998). Using the Waring and Running?s equation, Rh would be 1061 at HDF11 g C m-2 yr-1 and 723 g C m-2 yr-1 at HDF00.  Noormets et al. (2012) clearly showed that the Re partitioning in the first years immediately after a disturbance is different than in the later years, as a large fraction of Re comes from the decomposition of coarse woody debris and Ra is a small fraction. Special care should be 64  taken when using models in the first few years following a disturbance as the global ratios, Rh : Rs (e.g., Bond-Lamberty et al., 2004) or Ra : GEP (e.g., Waring and Running, 1998) might not hold or might vary between stand-specific disturbances. Noormets et al. (2012) found that Ra : GEP increased significantly with age, while Rh : Rs was significantly higher immediately after the disturbance. Harmon et al. (2011) found that variations in Rh : Re depended largely on the forest type and climate, although after a disturbance Rh could account for 100% of Re (Wang et al., 2002).   3.7.3 Comparison of HDF00 and HDF11  There can be a few reasons for the differences in C fluxes following harvesting at HDF11 and HDF00. Vegetation recovery was slower at HDF11 than at HDF00. Humphreys et al. (2005) reported considerable growth of the vegetation and seedlings over the three years of measurements at HDF00, whereas in the case of HDF11 the trees grew very slowly. Both sites started with an LAI of 0.2 m2 m-2 from small 1-year-old Douglas-fir seedlings (average height of 40 cm and basal diameter of 1.02 cm) and sparse understory vegetation left from the previous Douglas-fir forest. By the end of the second year, the seedlings at HDF11 had an average height of 43.7 cm and basal diameter of 1.16 cm, whereas the seedlings at HDF00 had an average height of 62 cm and basal diameter of 1.53 cm (Humphreys et al., 2005). The establishment and growth of pioneer species increased over the two years for both sites, but more at HDF00 where an LAI of 2.17 m2 m-2 was reached over the second summer whereas a maximum LAI of 1.78 m2 m-2 was reached at HDF11. HDF11 is located on a slope and drainage at the site resulted in a patchy surface of wet and dry areas following harvesting, which seems to have significantly affected vegetation recovery and tree growth. Other possible factors influencing the regrowth are 65  the higher elevation of 300 m.a.s.l. at HDF11 compared to 175 m.a.s.l. at HDF00 and the slope facing ENE at HDF11, which could have resulted in less solar irradiance at the surface.  Lower respiratory fluxes at HDF00 can possibly have resulted from the coarser soil texture with a greater mineral fraction and less organic matter compared to HDF11 (Table 2). In addition, the forest at DF49 was fertilized in 2007. However, Jassal et al. (2010a) showed that fertilization had very little effect on Re, so this is not likely to explain the difference in respiratory rates between HDF00 and HDF11. It is important to note that the first year of measurements at HDF00 started in September rather than in May due to postponed instrument installation at the site (Humphreys et al., 2005). The first growing season was the biggest source of C at HDF11, and this growing season was missed at HDF00. However, even during the second year, the C source stayed consistently greater at HDF11, so this cannot entirely explain the difference. Taking the same time period from 1 September to 31 August as at HDF00, HDF11 would have been a source of 780 g C m-2  yr-1 over the first year and thus, still a much larger source than HDF00.  3.7.4 Implications for chronosequence studies The solid line in Figure 29 represents the fit for the Douglas-fir chronosequence (Humphreys et al., 2006; Jassal et al., 2010a) before measurements at HDF11. As in most Fluxnet studies, this chronosequence study was based on only one replicate for a particular stand age. Chronosequence studies such as this one and many others (Amiro, 2001; Anthoni et al., 2002; Chen et al., 2002; Clark et al., 2004; Goulden et al., 2011; Goulden et al., 2006; Humphreys et al., 2006; Kolari et al., 2004; Law et al., 2001; Litvak et al., 2003; Mkhabela et al., 2009; Rannik et al., 2002; Schulze et al., 1999) have been used to investigate the influence of a 66  stand-replacing disturbance on the C balance of a stand. Chronosequence studies are a useful tool to obtain long-term age-related data in a short period of time; however, the fundamental assumption of chronosequence studies has been shown to be invalid for some of those studies (Howard et al., 2004; Johnson and Miyanishi, 2008; Yanai et al., 2003) while in most other cases, it has not been tested.   Figure 29: Relationship of NEP to stand age for all sites in the Douglas-fir chronosequence.  With only one replicate for a specific stand-age, it is dangerous to generalize on the forest C balance following a stand-replacing disturbance, as there can be considerable spatial and temporal variability in vegetation recovery and respiratory fluxes among same-age sites, as the results from this study have indicated. Similarly, significant spatial variability in C fluxes has been found for sites of less than 20 years of age in the same geographical area in Wisconsin (Amiro et al., 2010). Spatial variability in the early years following a disturbance has been shown 0 10 20 30 40 50 60 70-1000-800-600-400-2000200400600800NEP (g C m-2  yr-1 )Age of stand (yr)HDF00HDF88DF49SinkSourceHDF1167  to be particularly important for sites with slower recovery (Amiro et al., 2010), such as the Vancouver Island sites.  Differences between sites are related to the complex interactions between site microclimate, soil characteristics, nutritional status, management history, belowground and aboveground respiration processes, and the nature and speed with which the vegetation returns (Harmon et al., 2011; Thornton et al., 2002). These factors influence the forest growth and carbon cycling and, despite the fact that some of these differences are often considered negligible in chronosequence studies, this study showed that they can significantly affect the C balance of same age sites within an ecozone. Multiple trajectories are possible following a stand-replacing disturbance. At the seedling age, competition and/or facilitation by the understory and pioneer species can either decrease or increase the establishment and regrowth of the seedlings and this has an important effect on the long-term C sequestration of the stand. HDF11 will likely never be as productive as HDF00 as the seedlings have had a high mortality in the first years and regrowth has been minimal. Rh is the main component of the C balance after a stand-replacing disturbance and its response to climate variations has been found to be greater in the first 20 years following a disturbance (Noormets et al., 2007; Noormets et al., 2009). Rh is highly influenced by the nature of the site and disturbance (Harmon et al., 2011). Immediately after a disturbance, changes in Rh depend on the type and severity of the disturbance, the degree to which the disturbance removes the carbon legacy, and changes the relative size of the carbon pools contributing to Rh (Harmon et al., 2011). In this study, different respiratory fluxes were found between HDF00 and HDF11 despite the similarity of the sites and disturbances.  68  In order to account for and predict the consequences of stand-replacing disturbances such as harvesting on regional and global C budgets, special care should be taken to not trust blindly data from a single site and generalize it to all other same-age sites within that ecozone. Replicated C flux measurements are required to validate the measurements and ensure accurate C budgets, as errors at the measurement level then propagate to remote sensing indices and process- and inventory-based models. 69  4 Conclusions This study examined the use of a chronosequence on Vancouver Island to study C and energy balances from pre- to post-harvest. C and energy balances were measured using the EC technique for two years following harvesting at a recently harvested site. The C and energy balances measured were then compared to the pre-harvest measurements at the same site, and to post-harvest measurements at another previously harvested site in the chronosequence.  Net radiation decreased from pre- to post-harvest due to an increase in albedo and emitted longwave radiation. The albedo increased because of the deeper and longer lasting snow cover, whereas the longwave emitted increased due to higher surface temperatures. The annual Bowen ratio increased following harvesting due to reduced evapotranspiration following the removal of the forest canopy. Air temperatures were found to be similar between mature and harvested sites, whereas soil temperatures were found to be consistently greater at the two harvested sites. The recently harvested stand, HDF11, was a strong source of C during the first two years following harvesting due to high respiratory rates and slow recovery of vegetation, although this source of C was weaker during the second year due to the regrowth of plant cover. In the first year, HDF11 was a strong source during the summer and a weak source during the winter, as NEP was driven mainly by Re. HDF11 was a more constant C source during the second year due to GEP balancing Re during the summer. No relationship between GEP and PAR could be found in the first year of measurements as photosynthesis was almost nonexistent, whereas the recovery was observed in the second year with a well-defined hyperbolic relationship. Re followed a logistic relationship with soil temperature in both years and no relationship to soil water content was found.   70  From pre- to post-harvest, the stand transitioned from being a moderate sink of C to being a strong source of C. GEP markedly declined due to the removal of the canopy, while Re decreased only slightly. The two post-harvest sites were very different C sources. The previously clearcut harvested site (HDF00) was a weaker source of C due to lower respiration rates and faster vegetation recovery following the disturbance.  This study showed significant differences between two clearcut harvested sites and showed the importance of replicating measurements in same-age sites within an ecozone. However, this study only had two replicates for a same-age site and thus, more same-age sites would be required to properly quantify the variability which can occur within the ecosystem. Furthermore, despite the fact that the climate in the years compared in this study was similar, ideally C and energy exchanges should be measured simultaneously in different aged-stands to ensure that climate differences are not influencing the results.  Considering the global impact of disturbances on forest ecosystems, this work contributes to a better understanding of the impacts of harvesting on C and energy balances. In order to account for and predict the consequences of stand-replacing disturbances such as harvesting on regional and global C and energy budgets, replicated measurements are crucial to ensure accurate budgets, as errors at the measurement level then propagate to remote sensing indices, and process- and inventory-based models.  An accurate quantification of C and energy balances following disturbances is important in order to improve forest management strategies towards increased C sequestration and climate change mitigation.    71  References Addison, J.A., Trofymow, J.A. and Marshall, V.G., 2003. Functional role of Collembola in successional coastal temperate forests on Vancouver Island, Canada. Applied Soil Ecology, 24(3): 247-261. Amiro, B.D., 2001. Paired-tower measurements of carbon and energy fluxes following disturbance in the boreal forest. Global Change Biology, 7(3): 253-268. Amiro, B.D. et al., 2010. Ecosystem carbon dioxide fluxes after disturbance in forests of North America. Journal of Geophysical Research-Biogeosciences, 115. Amiro, B.D. et al., 2006. Carbon, energy and water fluxes at mature and disturbed forest sites, Saskatchewan, Canada. Agricultural and Forest Meteorology, 136(3-4): 237-251. Anthoni, P.M. et al., 2002. Seasonal differences in carbon and water vapor exchange in young and old-growth ponderosa pine ecosystems. Agricultural and Forest Meteorology, 111(3): 203-222. Aubinet, M. et al., 2000. Estimates of the annual net carbon and water exchange of forests: the EUROFLUX methodology. Advances in Ecological Research, 30: 113-175. Baldocchi, D., 2008. Breathing of the terrestrial biosphere: lessons learned from a global network of carbon dioxide flux measurement systems. Australian Journal of Botany, 56(1): 1-26. Baldocchi, D. et al., 2001. FLUXNET: A new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities. Bulletin of the American Meteorological Society, 82(11): 2415-2434. Baldocchi, D., Kelliher, F.M., Black, T.A. and Jarvis, P., 2000. Climate and vegetation controls on boreal zone energy exchange. Global Change Biology, 6: 69-83. Baldocchi, D., Valentini, R., Running, S., Oechel, W. and Dahlman, R., 1996. Strategies for measuring and modelling carbon dioxide and water vapour fluxes over terrestrial ecosystems. Global Change Biology, 2(3): 159-168. Baldocchi, D.D., 2003. Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: past, present and future. Global Change Biology, 9(4): 479-492. Barford, C.C. et al., 2001. Factors controlling long- and short-term sequestration of atmospheric CO2 in a mid-latitude forest. Science, 294(5547): 1688-1691. Barr, A.G. et al., 2007. Climatic controls on the carbon and water balances of a boreal aspen forest, 1994-2003. Global Change Biology, 13(3): 561-576. Barr, A.G. et al., 2004. Inter-annual variability in the leaf area index of a boreal aspen-hazelnut forest in relation to net ecosystem production. Agricultural and Forest Meteorology, 126(3-4): 237-255. Barr, A.G., Morgenstern, K., Black, T.A., McCaughey, J.H. and Nesic, Z., 2006. Surface energy balance closure by the eddy-covariance method above three boreal forest stands and implications for the measurement of the CO2 flux. Agricultural and Forest Meteorology, 140(1-4): 322-337. Barr, A.G., van der Kamp, G., Black, T.A., McCaughey, J.H. and Nesic, Z., 2012. Energy balance closure at the BERMS flux towers in relation to the water balance of the White Gull Creek watershed 1999?2009. Agricultural and Forest Meteorology, 153(0): 3-13. 72  Betts, A.K. and Ball, J.H., 1997. Albedo over the boreal forest. Journal of Geophysical Research-Atmospheres, 102(D24): 28901-28909. Black, T.A. et al., 1996. Annual cycles of water vapour and carbon dioxide fluxes in and above a boreal aspen forest. Global Change Biology, 2(3): 219-229. Black, T.A., Jassal, R.S. and Fredeen, A.L., 2008. Carbon sequestration in British Columbia?s forests and management options: a white paper on forestry, Pacific Institute for Climate Solutions, Victoria, BC. Blanken, P.D. et al., 1997. Energy balance and canopy conductance of a boreal aspen forest: Partitioning overstory and understory components. Journal of Geophysical Research: Atmospheres, 102(D24): 28915-28927. Bond-Lamberty, B., Wang, C. and Gower, S.T., 2004. A global relationship between the heterotrophic and autotrophic components of soil respiration? Global Change Biology, 10(10): 1756-1766. Brooks, A. and Farquhar, G.D., 1985. Effect of temperature on the CO2/O2 specificity of ribulose-1,5-bisphosphate carboxylase/oxygenase and the rate of respiration in the light. Planta, 165(3): 397-406. Brown, J.M., 1972. Effect of clearcutting a black spruce bog on net radiation. Forest Science, 18(4): 273-277. Cai, T.B., 2006. Analysis of the net ecosystem exchange of CO2 in a 56-year-old coastal Douglas-fir stand: its relation to temperature, soil moisture and photosynthetically active radiation. Unpub. PhD Thesis, University of British Columbia. Carrara, A. et al., 2003. Net ecosystem CO2 exchange of mixed forest in Belgium over 5 years. Agricultural and Forest Meteorology, 119(3-4): 209-227. Chen, B. et al., 2009. Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements. Boundary-Layer Meteorology, 130(2): 137-167. Chen, J.M. et al., 2006. Leaf area index measurements at Fluxnet-Canada forest sites. Agricultural and Forest Meteorology, 140(1?4): 257-268. Chen, J.Q. et al., 2002. Biophysical controls of carbon flows in three successional Douglas-fir stands based on eddy-covariance measurements. Tree Physiology, 22(2-3): 169-177. Chen, J.Q. et al., 2004. Net ecosystem exchanges of carbon, water, and energy in young and old-growth Douglas-fir forests. Ecosystems, 7(5): 534-544. Clark, K.L., Gholz, H.L. and Castro, M.S., 2004. Carbon dynamics along a chronosequence of slash pine plantations in north Florida. Ecological Applications, 14(4): 1154-1171. Clements, F.E., 1916. Plant succession: an analysis of the development of vegetation. Carnegie Institution of Washington, Washington, 512 pp. Davidson, E.A., Belk, E. and 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 Biology, 4(2): 217-227. Denman, K.L. et al., 2007. Couplings between changes in the climate system and biogeochemistry. In: S. Solomon et al. (Editors), Climate change 2007: the physical science basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK, pp. 500-584. 73  Drewitt, G.B. et al., 2002. Measuring forest floor CO2 fluxes in a Douglas-fir forest. Agricultural and Forest Meteorology, 110(4): 299-317. Dunn, A.L., Barford, C.C., Wofsy, S.C., Goulden, M.L. and Daube, B.C., 2007. A long-term record of carbon exchange in a boreal black spruce forest: means, responses to interannual variability, and decadal trends. Global Change Biology, 13(3): 577-590. Emmel, C., Paul-Limoges, E., Black, T. and 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 Meteorology: 1-31. Falge, E. et al., 2001. Gap filling strategies for defensible annual sums of net ecosystem exchange. Agricultural and Forest Meteorology, 107(1): 43-69. Federer, C.A., 1970. Measuring forest evapotranspiration. Theory and problems. (Ne-165): 25. Finnigan, J.J., Clement, R., Malhi, Y., Leuning, R. and Cleugh, H.A., 2003. A re-evaluation of long-term flux measurement techniques - Part I: Averaging and coordinate rotation. Boundary-Layer Meteorology, 107(1): 1-48. Foken, T., Aubinet, M. and Leuning, R., 2012. The eddy covariance method. In: M. Aubinet, T. Vesala and D. Papale (Editors), Eddy Covariance. Springer Atmospheric Sciences. Springer Netherlands, pp. 1-19. Gholz, H.L. and Clark, K.L., 2002. Energy exchange across a chronosequence of slash pine forests in Florida. Agricultural and Forest Meteorology, 112(2): 87-102. Goulden, M.L. et al., 1997. Physiological responses of a black spruce forest to weather. Journal of Geophysical Research-Atmospheres, 102(D24): 28987-28996. Goulden, M.L. et al., 2011. Patterns of NPP, GPP, respiration, and NEP during boreal forest succession. Global Change Biology, 17(2): 855-871. Goulden, M.L., Munger, J.W., Fan, S.-M., Daube, B.C. and Wofsy, S.C., 1996. Measurements of carbon sequestration by long-term eddy covariance: methods and a critical evaluation of accuracy. Global Change Biology, 2(3): 169-182. Goulden, M.L. et al., 2006. An eddy covariance mesonet to measure the effect of forest age on land-atmosphere exchange. Global Change Biology, 12(11): 2146-2162. Grant, R.F. et al., 2010. Net ecosystem productivity of temperate and boreal forests after clearcutting?a Fluxnet-Canada measurement and modelling synthesis. Tellus B, 62(5): 475-496. Griffis, T.J. et al., 2004. Seasonal variation and partitioning of ecosystem respiration in a southern boreal aspen forest. Agricultural and Forest Meteorology, 125(3?4): 207-223. Griffis, T.J. et al., 2003. Ecophysiological controls on the carbon balances of three southern boreal forests. Agricultural and Forest Meteorology, 117(1?2): 53-71. Harmon, M.E., Bond-Lamberty, B., Tang, J.W. and Vargas, R., 2011. Heterotrophic respiration in disturbed forests: A review with examples from North America. Journal of Geophysical Research-Biogeosciences, 116. Harmon, M.E., Ferrell, W.K. and Franklin, J.F., 1990. Effects on carbon storage of conversion of old-growth forests to young forests. Science, 247(4943): 699-702. Hilker, T. et al., 2010. Remote sensing of photosynthetic light-use efficiency across two forested biomes: Spatial scaling. Remote Sensing of Environment, 114(12): 2863-2874. Hornbeck, J.W., 1970. Radiant energy budget of clearcut and forested sites in West-Virginia. Forest Science, 16(2): 139-145. 74  Howard, E.A., Gower, S.T., Foley, J.A. and Kucharik, C.J., 2004. Effects of logging on carbon dynamics of a jack pine forest in Saskatchewan, Canada. Global Change Biology, 10(8): 1267-1284. Humphreys, E.R., 2004. Net ecosystem production of three coastal Douglas-fir stands at different stages of development after harvesting. Unpub. PhD Thesis, University of British Columbia, 122 pp. Humphreys, E.R. et al., 2003. Annual and seasonal variability of sensible and latent heat fluxes above a coastal Douglas-fir forest, British Columbia, Canada. Agricultural and Forest Meteorology, 115(1-2): 109-125. Humphreys, E.R. et al., 2006. Carbon dioxide fluxes in coastal Douglas-fir stands at different stages of development after clearcut harvesting. Agricultural and Forest Meteorology, 140(1?4): 6-22. Humphreys, E.R., Black, T.A., Morgenstern, K., Li, Z. and Nesic, Z., 2005. Net ecosystem production of a Douglas-fir stand for 3 years following clearcut harvesting. Global Change Biology, 11(3): 450-464. Irvine, J., Law, B.E., Martin, J.G. and Vickers, D., 2008. Interannual variation in soil CO2 efflux and the response of root respiration to climate and canopy gas exchange in mature ponderosa pine. Global Change Biology, 14(12): 2848-2859. Janssens, I.A. et al., 2001. Productivity overshadows temperature in determining soil and ecosystem respiration across European forests. Global Change Biology, 7(3): 269-278. Jarvis, P., James, G. and Landsberg, J., 1976. Coniferous forest. In: J.L. Monteith (Editor), Vegetation and the Atmosphere. Academic Press, pp. 171-240. Jarvis, P. and Stewart, J., 1979. Evaporation of water from plantation forest. In: E. Ford, D. Malcolm and J. Atterson (Editors), The Ecology of Even-Aged Forest Plantations. Institute of Terrestrial Ecology, Cambridge UK, pp. 327-349. Jarvis, P.G. and McNaughton, K.G., 1986. Stomatal control of transpiration - scaling up from leaf to region. Advances in Ecological Research, 15: 1-49. Jassal, R. et al., 2005. Relationship between soil CO2 concentrations and forest-floor CO2 effluxes. Agricultural and Forest Meteorology, 130(3?4): 176-192. Jassal, R.S. et al., 2010a. Impact of nitrogen fertilization on carbon and water balances in a chronosequence of three Douglas-fir stands in the Pacific Northwest. Agricultural and Forest Meteorology, 150(2): 208-218. Jassal, R.S. et al., 2007. Components of ecosystem respiration and an estimate of net primary productivity of an intermediate-aged Douglas-fir stand. Agricultural and Forest Meteorology, 144(1?2): 44-57. Jassal, R.S., Black, T.A., Nesic, Z. and Gaumont-Guay, D., 2012. Using automated non-steady-state chamber systems for making continuous long-term measurements of soil CO2 efflux in forest ecosystems. Agricultural and Forest Meteorology, 161(0): 57-65. Jassal, R.S., Black, T.A., Roy, R. and Ethier, G., 2011. Effect of nitrogen fertilization on soil CH4 and N2O fluxes, and soil and bole respiration. Geoderma, 162(1?2): 182-186. Jassal, R.S., Black, T.A., Spittlehouse, D.L., Brummer, C. and Nesic, Z., 2009. Evapotranspiration and water use efficiency in different-aged Pacific Northwest Douglas-fir stands. Agricultural and Forest Meteorology, 149(6-7): 1168-1178. 75  Jassal, R.S., Black, T.A., Trofymow, J.A., Roy, R. and Nesic, Z., 2010b. Soil CO2 and N2O flux dynamics in a nitrogen-fertilized Pacific Northwest Douglas-fir stand. Geoderma, 157(3?4): 118-125. Johnson, E.A. and Miyanishi, K., 2008. Testing the assumptions of chronosequences in succession. Ecology Letters, 11(5): 419-431. Kaimal, J.C. and Finnigan, J.J., 1994. Atmospheric boundary layer flows: their structure and measurement. Oxford University Press, 289 pp. Kidston, J. et al., 2010. Energy balance closure using eddy covariance above two different land surfaces and implications for CO2 flux measurements. Boundary-Layer Meteorology, 136(2): 193-218. Kolari, P. et al., 2004. Carbon balance of different aged Scots pine forests in Southern Finland. Global Change Biology, 10(7): 1106-1119. Kormann, R. and Meixner, F., 2001. An analytical footprint model for non-neutral stratification. Boundary-Layer Meteorology, 99(2): 207-224. Kowalski, S., Sartore, M., Burlett, R., Berbigier, P. and Loustau, D., 2003. The annual carbon budget of a French pine forest (Pinus pinaster) following harvest. Global Change Biology, 9(7): 1051-1065. Krishnan, P., Black, T.A., Jassal, R.S., Chen, B.Z. and Nesic, Z., 2009. Interannual variability of the carbon balance of three different-aged Douglas-fir stands in the Pacific Northwest. Journal of Geophysical Research-Biogeosciences, 114. Kristensen, L., Mann, J., Oncley, S.P. and Wyngaard, J.C., 1997. How close is close enough when measuring scalar fluxes with displaced sensors? Journal of Atmospheric and Oceanic Technology, 14(4): 814-821. Law, B.E. et al., 2002. Environmental controls over carbon dioxide and water vapor exchange of terrestrial vegetation. Agricultural and Forest Meteorology, 113(1-4): 97-120. Law, B.E., Sun, O.J., Campbell, J., Van Tuyl, S. and Thornton, P.E., 2003. Changes in carbon storage and fluxes in a chronosequence of ponderosa pine. Global Change Biology, 9(4): 510-524. Law, B.E., Thornton, P.E., Irvine, J., Anthoni, P.M. and Van Tuyl, S., 2001. Carbon storage and fluxes in ponderosa pine forests at different developmental stages. Global Change Biology, 7(7): 755-777. Lee, X., Fuentes, J.D., Staebler, R.M. and Neumann, H.H., 1999. Long-term observation of the atmospheric exchange of CO2 with a temperate deciduous forest in southern Ontario, Canada. Journal of Geophysical Research: Atmospheres, 104(D13): 15975-15984. Lee, X. et al., 2011. Observed increase in local cooling effect of deforestation at higher latitudes. Nature, 479(7373): 384-387. Leitch, A., 2010. Summertime horizontal and vertical advective carbon dioxide fluxes measured in a closed-canopy Douglas-fir forest on a slope. Unpub. MSc Thesis, University of British Columbia, 107 pp. Litvak, M., Miller, S., Wofsy, S.C. and Goulden, M., 2003. Effect of stand age on whole ecosystem CO2 exchange in the Canadian boreal forest. Journal of Geophysical Research-Atmospheres, 108(D3). Luyssaert, S. et al., 2008. Old-growth forests as global carbon sinks. Nature, 455(7210): 213-215. 76  McCaughey, J.H., 1981. Impact of clearcutting of coniferous forest on the surface radiation balance. Journal of Applied Ecology, 18(3): 815-826. McCaughey, J.H., 1985. A radiation and energy-balance study of mature forest and clear-cut sites. Boundary-Layer Meteorology, 32(1): 1-24. Mkhabela, M.S. et al., 2009. Comparison of carbon dynamics and water use efficiency following fire and harvesting in Canadian boreal forests. Agricultural and Forest Meteorology, 149(5): 783-794. Morgenstern, K. et al., 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(3?4): 201-219. Noormets, A., Chen, J. and Crow, T.R., 2007. Age-dependent changes in ecosystem carbon fluxes in managed forests in northern wisconsin, USA. Ecosystems, 10(2): 187-203. Noormets, A., Chen, J., Gu, L. and Desai, A., 2009. The phenology of gross ecosystem productivity and ecosystem respiration in temperate hardwood and conifer chronosequences. In: A. Noormets (Editor), Phenology of Ecosystem Processes. Springer, New York, pp. 59-85. Noormets, A. et al., 2008. Drought during canopy development has lasting effect on annual carbon balance in a deciduous temperate forest. New Phytologist, 179(3): 818-828. Noormets, A. et al., 2012. The role of harvest residue in rotation cycle carbon balance in loblolly pine plantations. Respiration partitioning approach. Global Change Biology, 18(10): 3186-3201. Odum, E.P., 1969. Strategy of ecosystem development. Science, 164(3877): 262-270. Oke, T.R., 1987. Boundary layer climates. Routledge, New York, 435 pp. Paul-Limoges, E., Christen, A., Coops, N.C., Black, T.A. and Trofymow, J.A., 2013. Estimation of aerodynamic roughness of a harvested Douglas-fir forest using airborne LiDAR. Remote Sensing of Environment, 136: 225-233. Pilegaard, K., Hummelshoj, P., Jensen, N.O. and Chen, Z., 2001. Two years of continuous CO2 eddy-flux measurements over a Danish beech forest. Agricultural and Forest Meteorology, 107(1): 29-41. Pojar, J., Klinka, K. and Demarchi, D.A., 1991. Coastal western hemlock zone. In: D. Meidinger, Pojar, J., (Editor), Ecosystems of British Columbia. BC Special Report Series No 6 Victoria: BC Ministry of Forests, pp. 95-111. Rannik, ?. et al., 2002. Fluxes of carbon dioxide and water vapour over Scots pine forest and clearing. Agricultural and Forest Meteorology, 111(3): 187-202. Rannik, ?. et al., 2012. Footprint Analysis, Eddy Covariance. Springer, pp. 211-261. Reichstein, M. et al., 2005. On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm. Global Change Biology, 11(9): 1424-1439. Schimel, D.S. et al., 2001. Recent patterns and mechanisms of carbon exchange by terrestrial ecosystems. Nature, 414(6860): 169-172. Schulze, E.D. et al., 1999. Productivity of forests in the Eurosiberian boreal region and their potential to act as a carbon sink - a synthesis. Global Change Biology, 5(6): 703-722. Sellers, P. et al., 1995. The Boreal Ecosystem-Atmosphere Study (BOREAS)-an overview and early results from the 1994 field year. Bulletin of the American Meteorological Society, 76(9): 1549-1577. 77  Smithwick, E.A.H., Harmon, M.E., Remillard, S.M., Acker, S.A. and Franklin, J.F., 2002. Potential upper bounds of carbon stores in forests of the Pacific Northwest. Ecological Applications, 12(5): 1303-1317. Sprugel, D., 1985. Natural disturbance and ecosystem energetics. In: S. Pickett and P. White (Editors), The Ecology of Natural Disturbance and Patch Dynamics. Academic Press, New York, pp. 335-352. Stenberg, P., DeLucia, E.H., Schoettle, A.W. and Smolander, H., 1995. Photosynthetic light capture and processing from the cell to canopy. In: W. Smith and T. Hinckley (Editors), Resource Physiology of Conifers. Academic Press, New York, pp. 3-38. Stoy, P.C. et al., 2006. An evaluation of models for partitioning eddy covariance-measured net ecosystem exchange into photosynthesis and respiration. Agricultural and Forest Meteorology, 141(1): 2-18. Suyker, A.E. and Verma, S.B., 2001. Year-round observations of the net ecosystem exchange of carbon dioxide in a native tallgrass prairie. Global Change Biology, 7(3): 279-289. Tanner, C.B. and Thurtell, G.W., 1969. Anemoclinometer measurements of reynolds stress and heat transport in the atmospheric surface layer: final report. University of Wisconsin, Madison, WI. Thornton, P.E. et al., 2002. Modeling and measuring the effects of disturbance history and climate on carbon and water budgets in evergreen needleleaf forests. Agricultural and Forest Meteorology, 113(1-4): 185-222. Trofymow, J.A., 1998. Detrital carbon fluxes and microbial activity in successional Douglas-fir forests. Northwest Science, 72: 51-53. Turner, D.P., Koerper, G.J., Harmon, M.E. and Lee, J.J., 1995. A carbon budget for forests of the conterminous United States. Ecological Applications, 5(2): 421-436. Urbanski, S. et al., 2007. Factors controlling CO2 exchange on timescales from hourly to decadal at Harvard Forest. Journal of Geophysical Research-Biogeosciences, 112(G2). Vickers, D., Irvine, J., Martin, J.G. and Law, B.E., 2012. Nocturnal subcanopy flow regimes and missing carbon dioxide. Agricultural and Forest Meteorology, 152: 101-108. Villar, R., Held, A.A. and Merino, J., 1995. Dark leaf respiration in light and darkness of an evergreen and a deciduous plant species. Plant Physiology, 107(2): 421-427. Walker, L.R., Wardle, D.A., Bardgett, R.D. and Clarkson, B.D., 2010. The use of chronosequences in studies of ecological succession and soil development. Journal of Ecology, 98(4): 725-736. Wang, C., Bond-Lamberty, B. and Gower, S., 2002. Environmental controls on carbon dioxide flux from black spruce coarse woody debris. Oecologia, 132(3): 374-381. Waring, R.H. and Running, S., 1998. Forest ecosystems: analysis at mutiple scales. Academic Press, London UK. Webb, E.K., Pearman, G.I. and Leuning, R., 1980. Correction of flux measurements for density effects due to heat and water-vapor transfer. Quarterly Journal of the Royal Meteorological Society, 106(447): 85-100. Wilson, K. et al., 2002. Energy balance closure at FLUXNET sites. Agricultural and Forest Meteorology, 113(1-4): 223-243. Wofsy, S.C. et al., 1993. Net exchange of CO2 in a mid-latitude forest. Science, 260(5112): 1314-1317. 78  Wohlfahrt, G., Bahn, M., Haslwanter, A., Newesely, C. and Cernusca, A., 2005. Estimation of daytime ecosystem respiration to determine gross primary production of a mountain meadow. Agricultural and Forest Meteorology, 130(1?2): 13-25. Yanai, R.D., Currie, W.S. and Goodale, C.L., 2003. Soil carbon dynamics after forest harvest: an ecosystem paradigm reconsidered. Ecosystems, 6(3): 197-212. Zha, T. et al., 2009. Carbon sequestration in boreal jack pine stands following harvesting. Global Change Biology, 15(6): 1475-1487.   79  Appendices Appendix A  : Site pictures A.1 Pre- to post-harvest surface cover  Figure 30: Photographs showing the forest cover at DF49 and the surface cover immediately following harvesting, which became HDF11. The photographs were taken from the top of the DF49 tower at 40-m height and the view is towards east-southeast.            80  A.2 Instrumentation at HDF11  Figure 31: Post-harvest surface residue and micrometeorological tower in June 2011 looking downslope towards the east-northeast.    Figure 32: Eddy-covariance sensors at 4.5- m height at HDF11. Photograph shows the Gill Instruments R3 sonic anemometer, the fine wire thermocouple (right) and the air intake to the LI-COR Inc. LI-7000 infrared gas analyzer (bottom). EC system Net radiometer PAR sensors 81  A.3 Year 2 in pictures     Figure 33: Year 2 in pictures. These 12 photographs represent one picture per month over the second year of measurements (in order from May 2012 to April 2013). The photographs are taken from the top of the DF49 tower at 40-m height and the view is towards the northwest.     May 2012 June 2012 July 2012August 2012 September 2012 October 2012 November  2012 December 2012 January  2013 February  2013 March  2013 April  201382  As Figure 33 shows, plants started to reestablish at HDF11 in May and grew until July. Drier surface conditions due to reduced precipitation in July and August led to plant senescence. Heavy and constant precipitation during the fall led to very wet and saturated soil conditions, resulting in a darker soil surface. HDF11 was then mostly or partly covered in snow in January and February, and then started drying out in March for the beginning of the new growing season in April.   Appendix B  : Paper on estimation of aerodynamic roughness As mentioned in the preface, the following paper was published during my MSc:  Paul-Limoges, E., Christen, A., Coops, N.C., Black, T.A. and Trofymow, J.A., 2013. Estimation of aerodynamic roughness of a harvested Douglas-fir forest using airborne LiDAR. Remote Sensing of Environment, 136: 225-233.  83   Estimation of aerodynamic roughness of a harvested Douglas-fir forestusing airborne LiDARE. Paul-Limoges a,b,?, A. Christen a, N.C. Coops c, T.A. Black b, J.A. Trofymow da Department of Geography, University of British Columbia, 1984 West Mall, Vancouver, BC V6T 1Z2, Canadab Faculty of Land and Food Systems, University of British Columbia, 2424 Main Mall, Vancouver, BC V6T 1Z4, Canadac Faculty of Forestry, University of British Columbia, 2357 Main Mall, Vancouver, BC V6T 1Z4, Canadad Canadian Forest Service, Natural Resources, Pacific Forestry Centre, 506 West Burnside Road, Victoria, BC V8Z 1M5, Canadaa b s t r a c ta r t i c l e i n f oArticle history:Received 18 September 2012Received in revised form 1 May 2013Accepted 5 May 2013Available online 2 June 2013Keywords:LiDARMicrometeorologyDigital terrain modelRoughness lengthThe aerodynamic roughness length (z0) is a key variable for the parameterization of momentum, mass andheat exchanges between land surfaces and the atmosphere. Its estimation however is complicated due tothe large number of input variables such as height and arrangement of roughness elements on the surface,and its measurement relies on complex micrometeorological instrumentation that is typically unavailable.One remote sensing technology well suited to measuring the height of objects is light detection and ranging(LiDAR). This study demonstrates the use of pre- and post-harvest LiDAR data to quantify the aerodynamicroughness of a post-harvest forest surface. LiDAR data were acquired before and after clearcut harvestingof a 77-ha Douglas-fir dominated site on Vancouver Island, for which a micrometeorological tower provideddirect year-long measurements of shear or Reynolds stress (i.e., momentum flux) and wind speed, thus per-mitting the independent assessment of z0 using the logarithmic wind profile equation.The LiDAR data were used to estimate z0 based on the standard deviation of roughness element heights withinthe source areas of themicrometeorological tower. Estimated z0 from the LiDAR analysis comparedwell to z0 cal-culated using the micrometeorological measurements. The standard deviation of roughness element height es-timated from the LiDAR analysis resulted in z0 = 0.13 ? 0.01 m (mean ? SD) for neutral atmosphericstability conditions and z0 = 0.13 ? 0.01 m for all stability conditions. The value of z0 calculated using the log-arithmicwind profile equationwas 0.13 ? 0.13 m for neutral conditions and 0.12 ? 0.30 m for all stability con-ditions after applying diabatic profile corrections. The results from this study demonstrate the potential of usingLiDAR data to estimate z0 across large areas and in complex situations where direct measurements of z0 areimpossible.? 2013 Elsevier Inc. All rights reserved.1. IntroductionThe aerodynamic roughness length (z0) is a key variable for theparameterization of momentum, mass and heat exchange betweenland surfaces and the atmosphere, and is defined as the height atwhich the vertical profile of mean horizontal wind speed extrapolatesto zero (Kaimal & Finnigan, 1994). Many variables influence z0 includ-ing the height, geometry, density and arrangement of surface rough-ness elements (Garratt, 1992; Lettau, 1969; Raupach, 1992, 1994;Shaw & Pereira, 1982) and, as a result, it cannot be quantified easily.Micrometeorological wind profile measurements can be used todirectly determine z0 (Stull, 1988). However, suchmeasurements are ex-pensive, site-specific and limited to the source area of themeasurements,thereby not allowing for the determination of z0 at a greater spatial scale.Wind tunnel and computational fluid dynamics experiments are used tounderstand the influence of variables such as spacing, density, heightand layout on z0 (e.g. Crago et al., 2012; Xie et al., 2008) andmorphomet-ricmethods have been used to relatemeasured dimensions of roughnesselements to z0 using algorithms (e.g. Grimmond & Oke, 1999). Such ex-periments are fundamental to the understanding of the influence of dif-ferent variables on z0, but they rely on laboratory simulations usingsimplified geometries, and thus are limitedwhen describing actual, com-plex surfaces, especially on larger spatial scales, and require validation inthe field. In recent years, an increasing number of models have been de-veloped to estimate the exchange between the surface and the atmo-sphere using passive remote sensing (e.g. Bastiaanssen et al., 1998;Colin et al., 2006a; Roerink et al., 2000; Su, 2002); many of these modelsrequire an accurate estimate of z0 that is often not available, in whichcase it is calculated as a simple ratio of spectral channels such as theNor-malizedDifference Vegetation Index (NDVI). Different formulations usedto parameterize z0 frompassive remote sensing are often applied beyondtheir intended range and on heterogeneous land surfaces, often resultingin significant errors in turbulent flux estimates (Colin et al., 2006b).Remote Sensing of Environment 136 (2013) 225?233? Corresponding author at: Department of Geography, University of British Columbia,249-1984 West Mall, Vancouver, BC V6T 1Z2, Canada. Tel.: +1 604 822 5654.E-mail address: eugeniepl@hotmail.com (E. Paul-Limoges).0034-4257/$ ? see front matter ? 2013 Elsevier Inc. All rights reserved.http://dx.doi.org/10.1016/j.rse.2013.05.007Contents lists available at SciVerse ScienceDirectRemote Sensing of Environmentj ourna l homepage: www.e lsev ie r .com/ locate / rse84   Advances in remote sensing over the past 20 years have providednew tools to directly measure the three dimensional structure of landsurfaces. Light detection and ranging (LiDAR) is capable of characteriz-ing the size and spatial distribution of surface roughness elements withan accuracy in the decimeter range. Despite this advantage, only a fewstudies have used LiDAR to characterize z0. Menenti and Ritchie(1994) were the first to use LiDAR to estimate z0 based on the geomet-rical regularity of vegetation canopies by multiplying the ratio of thestandard deviation of vegetation height to vegetation height in seg-ments of a transect by the average height of the vegetation along thetransect. They found reasonable agreement between their estimates ofz0 derived from LiDAR measurements (corrected for instrument noise)and calculated using the Monin?Obukhov similarity theory applied tomeasurements of horizontal wind speed profiles (Kustas et al., 1994).Brown and Hugenholtz (2011) adapted Menenti and Ritchie's methodfor use in a mixed grassland prairie, and found that the variation inroughness heights explained 76% of the variation in z0 derived from awind speed profile with five measurement heights. De Vries et al.(2003) extended the use of LiDAR data to estimate z0 for complex ter-rain consisting of coppice dunes covered with honey mesquite withbare interdunal areas. They also found good agreementwhen they com-pared their LiDAR-derived estimates of z0 to those derived from a verti-cal profile of horizontal wind speed measured at six heights. Colin andFaivre (2010) estimated z0 for a heterogeneous landscape using a com-bination of LiDAR and computational fluid dynamics models. Theydiscussed the need for footprint definitions and ground measurementsto validate their results. Tian et al. (2011) used a combination of LiDARand SPOT-5 spectral data with micrometeorological measurements totest fourmodels to parameterize z0. They showed that themodel gener-ated maps of z0 from LiDAR are superior to those from satellite opticalremote sensing data and suggested the use of high density LiDAR com-bined with source areas for eddy covariance (EC) towers to improvetheir validation.This paper reports the evaluation of the potential of using pre-harvest and post-harvest LiDAR data to estimate values of z0 for aharvested Douglas-fir forest in the Pacific Northwest. The evaluationwas achieved by comparing values of LiDAR-derived z0 with values cal-culated using micrometeorological measurements of shear or Reynoldsstress and wind speed.2. Background on LiDAR accuracy for terrain mappingLiDAR is widely used for mapping terrain elevations due to its highprecision and accuracy (e.g. Hodgson et al., 2005). The last returnpulse is often considered the best estimate of the ground surface; how-ever, this is not always the casewhen other obstacles such as vegetationcanopy, understory or shrub cover totally obscure the ground (Hodgsonet al., 2005). LiDAR systems can achieve an accuracy of 15 cmover openflat surfaces (Gomes Pereira & Janssen, 1999). However, for roughersurfaces such as dense canopies and understories, the accuracy is usual-ly lower due to the uncertainties in the reflections coming from the sur-face or vegetation. In dense canopies, there are also fewer pulsesreaching the ground, making it more difficult to extract a surface fromthese points. Greater data processing must be performed in such casesto classify the points as ground or non-ground.Reutebuch et al. (2003) tested the accuracy of digital terrainmodels (DTMs) generated in heavily forested 70-year-old Douglas-firforest by producing a high-resolution DTM from high density bareground LiDAR returns. In that study, they found that the DTM errorwas 0.16 ? 0.23 m (mean ? SD) for clearcuts, 0.18 ? 0.14 m forheavily thinned, 0.18 ? 0.18 m for lightly thinned, and 0.31 ?0.29 m for uncut stands. Although they found the DTM error to in-crease with stand density, the differences were small and the accuracyunder a dense forest canopy remained very high. They found that themean and standard deviation of the differences between the DTMand 347 ground points were 0.22 m and 0.24 m, respectively (rootmean square error (RMSE) = 0.32 m). These values were similar to,but less than, those of Kraus and Pfeifer (1998) who compared theirDTM results for a wooded area in Austria to 466 ground points(RMSE = 0.57 m). Kraus and Pfeifer (1998) also found a systematicoverestimation of elevation of 0.20 m for the DTM.3. Methods3.1. Study siteThe study site (49?52?9.54? N, 125?20?9.00? W, 300 m.a.s.l.) is locat-ed near Campbell River on the east coast of Vancouver Island, Canadaand it has a downward slope of about 5? towards the east-northeast.The site is located in the dry maritime Coastal Western Hemlockbiogeoclimatic subzone, with an average annual precipitation of1500 mm and a mean annual temperature of 9.1 ?C (Pojar et al.,1991). These climate conditions contribute significantly to making theforests in this subzone the most productive in North America (Turneret al., 1995).Prior to harvesting in 2011, the standwas a dense (1100 stems ha?1)coniferous stand of 62 years-of-agewith tree heights varying between 30and 35 m (Hilker et al., 2010). The pre-harvest stand was composed of80% Douglas-fir (Pseudotsuga menziesii), 17% western redcedar (Thujaplicata Donn), and 3% western hemlock (Tsuga heterophylla (Raf.) Sarg.)with only a sparse understory (Ferster et al., 2011; Humphreys et al.,2006; Morgenstern et al., 2004). The leaf area index was 7.3 m2 m?2(Chen et al., 2006). This second growth stand established followingclearcut harvesting (1937, 1938 and 1943) of the original old-growthstand (Ferster et al., 2011). A large section of the site (77 ha) washarvested from January to March 2011, and replanted in April 2011with 97% Douglas-fir and 3% Sitka spruce (Picea sitchensis). The post-harvest surface roughness is now composed of old stumps from the1943 clearcut, recent stumps from the 2011 clearcut, harvest residuepiles, logs and holes from the harvesting machinery (Fig. 1).3.2. Light detection and ranging (LiDAR)3.2.1. LiDAR data acquisitionPre-harvest LiDAR data were acquired on 14 August 2008, using aLeica ALS50-II able to record up to 4 returns per laser pulse. The flyingaltitude was approximately 900 m and the sensor pulse rate was110 kHz. The estimated global positioning system (GPS) accuracy ofthe sensor was 0.02, 0.03 and 0.05 m in x, y and z, respectively. TheFig. 1. Post-harvest surface residue and micrometeorological tower in June 2011 lookingdownslope towards east-northeast. The instrument on the boom to the left is the ultra-sonic anemometer?thermometer used to measure horizontal wind speed and Reynoldsstress.226 E. Paul-Limoges et al. / Remote Sensing of Environment 136 (2013) 225?23385   raw data point cloud had an average density of 3.74 points m?2.Post-harvest LiDAR data were acquired on 18 August 2011. The flyingaltitude was approximately 500 m. The static GPS accuracy was lessthan 2 cm and the kinematic GPS accuracy was less than 5 cm. Thepixel resolution was 15 cm and the raw data point cloud had an aver-age point density of 15 points m?2.3.2.2. Derivation of DTMs from LiDAR dataThe software program FUSION v 2.90 from theUnited States Depart-ment of Agriculture Forest Service was used to produce the DTMs(McGaughey, 2010). DTMs were produced by filtering the LiDAR datapoints to eliminate any non-ground points, including vegetation andbuildings, and by interpolating the bald forest floor surface. The filteringalgorithmwas adapted from Kraus and Pfeifer (1998) and is based on alinear prediction (Kraus & Mikhail, 1972) with an individual accuracyfor each LiDAR point. A surface is first computed with equal weightsfor all points; there is a greater probability for ground points to have anegative vertical residual and for vegetation points to have a small neg-ative or positive vertical residual. Each vertical residual (vi) is then usedtoweight each point based on their distance and direction to the surfaceusing the following function,pi ?1 vi ? g11? a vi?g? ?b  g b vi ? g ? w;0 g ? w b vi8><>:?1?where a and b determine the steepness of the weighting function, g is abelow ground shift value (in m), and w is the above ground offset pa-rameter (in m). In general, a value of 1 for a, and 4 for b has beenshown to produce accurate results if the LiDAR dataset is of high enoughdensity. The value of g is computed in three different ways based on ahistogram of the residuals (Pfeifer et al., 1998) and plausibility rulesare applied to choose the best value for g. Points that are below the sur-face bymore than the shift value g get a value of pi = 1,whilew sets anupper limit for points to have an influence on the surface. Large positiveresiduals (g + w) are set to pi = 0. As a result, points with large nega-tive residuals influence the computed surface themost, points withme-dium residuals influence the surface less, and points with large positiveresiduals are eliminated.3.2.3. Derivation of roughness element height from LiDARStudies have found that the accuracy of DTMs is impacted by nearground features due to the small differences in height with respect tothe true ground surface, making it difficult to separate these nearground features from the ground (Estornell et al., 2011; Streutker &Glenn, 2006). Meng et al. (2010) reported that low vegetation is oftenignored by ground filters and included in the ground surface. For thatreason, DTMs were generated from the pre- and post-harvest LiDARdatasets. The roughness of the surface for the post-harvest DTM waspreserved by not undertaking any filtering, whereas the pre-harvestDTM was filtered to represent only the smooth ground surface withthe influence of plants and other near-surface, non-ground obstaclesremoved.The LiDAR-derived DTMs were exported to ArcGIS v 10 (ESRI Inc.,Redlands, CA, USA) for further analysis. To quantify the post-harvestroughness of the stand, the pre-harvest DTM was subtracted from thepost-harvest DTM. The surface following harvesting is covered withwoody residues and was also changed by the harvesting machinery(i.e. compaction and holes); subtracting a filtered post-harvest DTMwould have resulted in the loss of those post-harvest roughness fea-tures. The pre-harvest ground surface was covered with little plantcover and a more accurate surface could be interpolated from thesemeasurements due to the greater height difference between the vegeta-tion canopy and ground surface.3.2.4. Verification of LiDAR roughness element heightField measurement of roughness element heights were taken at thesite by recording the coordinate locations with a hand-held GPS andmeasuring the vertical height of roughness elements above or belowthe ground surface with a tape measure to minimize error due tosmall errors in elevations measured with the GPS. The 24 randomlydistributed check points included different roughness elements includ-ing soil surfaces, 1943 stumps, 2011 stumps, residue piles and holes.These ground points were overlaid on the LIDAR-derived roughness el-ement map to assess the vertical accuracy of the residual approach.3.2.5. Estimations of z0 from LiDAREstimations of z0 from the LiDAR residual were performed for eachgrid cell of the entire clearcut using an equation proposed by Menentiand Ritchie (1994). In this equation, z0 is estimated from the variabilityin roughness heights asz0;l ?1N? ?hi;j =hi;jh ih; ?2?where z0,l is the LiDAR estimate of z0;N is the number of 1 m ? 1 mgridcells within the source area; ?hi;j is the standard deviation of theLiDAR-derived roughness element height in grid cell i, j; hi,j is theLiDAR-derived roughness element height in grid cell i, j; andh is the av-erage LiDAR-derived roughness element height for the entire harvestedarea. The pixel size used in this analysis must be smaller than the char-acteristic width of the roughness elements in order to resolve themproperly. In order to calculate the standard deviation for a cell i, j, thecell was compared to the neighboring cells within a circle of 1 m (F1),2 m (F2), 3 m (F3) and 4 m (F4) radii to determine the influence of ra-dius size. As the post-harvest surface is characterized by negative andpositive roughness element heights, the absolute value of themost neg-ative roughness element height was added to all cells in order to makethem all positive. This approach mimics the influence of differentroughness element heights in a forest canopy, for example, where thesoil and the different roughness elements influence the effective z0; al-though the spaces between the different tree canopies can be seen asholes, they are still positive roughness elements.3.3. Direct micrometeorological measurements3.3.1. Aerodynamic measurementsMean horizontal wind speed u (in m s?1) and Reynolds stress ? (inN m?2) were measured at the site using an ultrasonic anemometer?thermometer (R3, Gill Instruments, Lymington, UK) at a height of4.5 m above the soil surface. Longitudinal (u), lateral (v) and vertical(w) components of the wind velocity vector were measured by the ul-trasonic anemometer?thermometer at 20 Hz for the period of 1 May2011 to 31 December 2011. Half-hour values of Reynolds stress (i.e.,momentum flux) ? were calculated from the covariance of the u andw components, after coordinate rotation to make mean v and w equalto zero, using ? ? ?u?w? , where ? is the mean air density in kg m?3,the overbar indicates the mean value and the primes indicate fluctua-tions from the mean (Stull, 1988). The wind stream line was approxi-mately parallel to the surface.3.3.2. Estimation of z0 from micrometeorological measurementsThe effective aerodynamic roughnesswas calculatedusing half-hourlymeasurements of mean wind speed and Reynolds stress by rearrangingthe logarithmic wind profile equation under near neutral conditions (sta-bility conditions are defined later in this section) to obtain an expressionfor z0 (Stull, 1988):z0;a ? z?d? ? exp?kuu ?3?227E. Paul-Limoges et al. / Remote Sensing of Environment 136 (2013) 225?23386   where z0,a is the aerodynamic estimate of z0, u is the mean vector windspeed, u* is the friction velocity, k is the von Karman constant (0.41), z isthe measurement height, and d is the zero-plane displacement height.Near neutral conditions are used because empirical diabatic corrections(Garratt, 1992) are not required in this case. Friction velocitieswere cal-culated from the Reynolds stress measurements using u?ffiffiffiffiffiffiffiffi?=?p. Cal-culations performed on half-hourly data with wind speeds less than1.5 m s?1 were removed as they did not provide sufficient mixing,and the low Reynolds stress caused larger errors in estimated z0 values.The zero-plane displacement height d was estimated to be 2/3 of themedian roughness element height (Monteith & Unsworth, 2008). Inthis study where d ? z, a small error in d would not have a big effecton the calculated z0. Sensitivity analysis indicated that a 1 m changein d resulted in a 0.02 m change in z0; an error as great as 0.5 m in dis very unlikely. Large harvest residue piles were burnt at the site 19 Oc-tober 2011, so the calculations were performed on the data collectedafter the burning of the piles, as they were also removed in the LiDARdata analysis.In addition, the effective aerodynamic roughness was also calcu-lated for stable and unstable conditions using the diabatic profileequation (Campbell & Norman, 1998; Yasuda, 1988),z0;a ? z?d? ? exp?kuu??M ; ?4?where?M is the profile diabatic correction factor. The diabatic correc-tion factor is zero for neutral conditions, but negative for unstable andpositive for stable conditions. It can be expressed in terms of the at-mospheric stability parameter (? = (z ? d)/L) where L, the Obukhovlength, is given byL ???u3kgw????5?in which g is the gravitational acceleration (m s?2), ? is themean abso-lute air temperature (K), and w??? the kinematic sensible heat flux(K m s?1) measured by the sonic anemometer?thermometer. The ex-pression for ?M in unstable flow (? b 0), i.e., enhanced turbulence, is?M ??1:2 ln1? 1?16?? ?1=22" #?6?while in stable flow (? > 0), i.e., suppressed turbulence, it is?M ? 6 ln 1? ?? ?: ?7?Near neutral conditions are defined as?0.005 b ? b 0.005 for calcu-lations in our analysis.3.3.3. Turbulent source areasIf the flow is in equilibrium with the surface momentum exchange,the roughness length upwind of the tower (the ?source area? of theultrasonic anemometer?thermometer) will affect the ratio of u to ?.This source area changes constantly with changing wind direction andturbulent state of the atmosphere. To estimate the instantaneous areathat influences the sampled ratio of wind to Reynolds stress, a2-dimensional gradient diffusion and crosswind dispersion model(Kormann & Meixner, 2001) was run for all half-hour periods for the2011 measurements made after harvesting at a 1-m grid resolutionover a domain of 1000 m by 1000 m. The inputs to the model were thedirectly measured wind direction at a height of 4.5 m, the standarddeviation of the lateral wind speed (?v) at 4.5 m and L. An a priori rea-sonable homogeneous value of z0 of 0.12 m (i.e., one tenth of the heightof the main roughness elements; Monteith & Unsworth, 2008) was as-sumed. Ensemble turbulent source areas were calculated for differentstability conditions and wind directions. The integrated source areaswere calculated as the average of all individual half-hour source areasfor the characteristics of interest (i.e., stability, wind direction, etc.) sim-ilar to Chen et al. (2009).3.4. Comparing roughness lengths derived from LiDAR andmicrometeorological measurementsWhile the LiDAR analysis provides an estimate of z0 for each grid cell,the micrometeorological method provides a single estimate that is aweighted spatial average over the entire instantaneous source area ofthe ultrasonic anemometer?thermometer. To allow a direct comparisonof the two datasets, the estimated z0,l for all grid cells was weighted bymultiplying them by the gridded footprint probability, or vertical fluxper unit point source ? (m?2) of the integrated source area of interest(Kormann & Meixner, 2001). A small fraction of the source area is pre-dicted to be outside the 1000 m by 1000 m domain for which thesource area model was run (typically b15%, except under stable condi-tions). This fraction was set to the average z0,l within the entire studyarea following Christen et al. (2011). This approach assumes that theharvested area continues in a similar pattern outside the modeled area.4. Results4.1. Quality of the LiDAR-derived DTMs and their subtractionsFig. 2 demonstrates that the pre- and post-harvest LiDAR-derivedDTMs of the ground surfaces were in good agreement; this can beseen by the similarity in the 10-m elevation contour lines as well as inthe larger relief features. The pre-harvest DTM characterized thesmooth ground surface to remove any plants, shrubs and other nearsurface obstacles, which would otherwise have been incorrectlyaccounted for in the DTM, thereby resulting in a higher ground surface.As discussed, in order to maintain the roughness elements of the site,the post-harvest DTM was not filtered. These results show the high ac-curacy of LiDAR-derived DTMs, even under dense Douglas-fir canopies.Fig. 3 shows the post-harvest residual roughness element height fol-lowing the subtraction of the pre-harvest DTM from the post-harvestDTM. This residual roughness element height is composed of the differ-ent roughness element heights for 1 ? 1 m cells with the influence oftopography removed due to the DTM subtraction. The highest rough-ness element areas occur near the roads and correspond to large pilesof logs that remained after harvesting. The lowest roughness elementareas are also located near the roads and often correspond to drainagechannels that were created during road construction. Except for nearthe roads, roughness element heights varied mostly between plus andminus 1 m. The machinery used to harvest the trees highly disturbedthe soil, resulting in large holes (about 1 m3) being created at the sitefrom erosion and drainage following harvesting. Above-ground rough-ness elements are composed of the 1943 stumps, 2011 stumps, residuepiles and logs. Fig. 3 illustrates that the DTM subtraction was successfulin removing the topographical features, as the non-existent 10-mcontours (i.e., there is only the 0-m contour separating the above- andbelow-ground features) show that there was no residual slope left.Heights measured at 24 locations at the site were used to validatethe roughness element heights obtained from the LiDAR data. These lo-cations included undisturbed soil surfaces, holes, 1943 stumps, 2011stumps and residue piles. The average absolute error for all locationheights was 20 cm and the median absolute error was 8 cm (Table 1).The largest errors were found for the old 1943 stumps that were alsopresent in the pre-harvest LiDAR dataset, and consequently, were likelyincluded in the pre-harvest DTM and lost after filtering and subtraction.The smallest errors were found for the soil that was undisturbed frompre- to post-harvest. Heightmeasurements need to bemade atmore lo-cations in a later study in order to better quantify the accuracy of theresults.228 E. Paul-Limoges et al. / Remote Sensing of Environment 136 (2013) 225?23387   4.2. Comparing roughness lengths derived from LiDAR andmicrometeorological measurementsFig. 4 shows the distribution of z0,a calculated from the micromete-orological measurements for neutral (?0.005 b ? b 0.005), stable(0.005 b ? b 1) and unstable (?1 b ? b ?0.005) conditions. Therewas very good agreement between the values ofmean z0,a for the differ-ent stability conditionswith 0.127 m for neutral conditions, 0.122 m forstable conditions, and 0.119 m for unstable conditions. The analysis fo-cused primarily on neutral and unstable conditions due to the smallerdomain covered by the source area compared to stable conditions. Alarge percentage of the source area was outside the domain under sta-ble conditions (on average about 38%) and, as a result, it did not permitan accurate z0,a analysis. Fig. 5 illustrates the cumulative source areas forneutral (?0.005 b ? b 0.005) and unstable (?1 b ? b ?0.005) condi-tions overlaying the logarithm of the estimated z0,l, i.e., from theLiDAR analysis. Only 10.7% and 12.2% of the turbulent source areaswere outside the domain under neutral and unstable conditions,respectively.The influence of wind direction under unstable and neutral condi-tions was also explored to assess the homogeneity of the site (Table 2).Although the average values of z0,a and z0,l for all wind direction sectorswere similar, there were significant differences among the differentwind direction sectors. As Table 2 indicates, northerly winds tended tohave a larger z0,a, which can be explained by the close proximity, about200 m away, to the non-harvested forest (30 m tall) in that direction(see Fig. 5b), so the flow has not yet readjusted to the smoother clearcutsurface. The influence of unstable and neutral stability conditions on thecalculated z0,a is shown in Table 3. Larger values of z0,a were found underunstable conditions due to the greater number of occurrences of unstableconditions when the wind was from the north.Spatially averaged values of z0,l obtained using a filter radius of2-m (F2) matched z0,a best. The filter's radius needs to capture thecharacteristic width of the roughness elements and their influenceon their surroundings. The larger the filter's radius, the higher theoverall variability in the DTM that influences ?h and the higher thecalculated z0,l. From our comparison we conclude that a filter radiusof about twice the size of the individual roughness elements gave400 mN0240260280300320340360240260280300320340360a)b)Fig. 2. Digital terrain models derived from LiDAR data for a) pre-harvest and b) post-harvest. The gray shading represents a hillshade showing the roughness of the DTMs and theyellow lines represent the 10-m contour lines.229E. Paul-Limoges et al. / Remote Sensing of Environment 136 (2013) 225?23388   the best agreement between z0,l and z0,a. Likely, the 2-m radius waslarge enough to include the morphometric dimension of the individ-ual roughness elements (i.e., the element and their influence on thesurrounding), as the dominant roughness elements were ~1 m inwidth. The 1-m filter (F1) was inadequate in accounting for the influ-ence of the roughness elements as in some cases the filter radius wasentirely covered by a roughness element. The 3-m and 4-m filters (F3and F4) included clusters of several roughness elements, so the arealweight given to the tallest elements in the cluster probably led to theoverestimation of z0,l. The values of the mean and standard deviationof z0,l for neutral and all stability conditions obtained using the F2 fil-ter were identical: 0.13 ? 0.01 m, while the respective values of z0,awere 0.13 ? 0.13 m and 0.12 ? 0.30 m (Table 4).5. Discussion5.1. Feasibility of generating DTMs under dense forest canopiesThe results from this study showed that high accuracy could beachieved when using LiDAR data to generate DTMs, even under denseDouglas-fir canopies. The subtraction of a pre-harvest DTM from apost-harvest DTM gave highly accurate results, with no topographicalfeatures remaining. Differences could have been expected between thetwo different LiDAR datasets, as many studies have found under-predictions (e.g. Adams & Chandler, 2002; Hodgson et al., 2005) orover-predictions (e.g. Kraus & Pfeifer, 1998) of elevations when usingLiDAR to generate DTMs. The subtraction of the two DTMs in this caseaccounted for this under- or over-prediction as there was no resultingslope. The accuracy of LiDAR datasets and their derived products is af-fected by many factors including sensor, aircraft platform, navigation,LiDAR point processing and the characteristics of the surface (Hodgsonet al., 2005). All products derived from the DTM such as slope areinfluenced by the same factor affecting elevation inaccuracies(Hodgson et al., 2005). However, it is possible that while the elevationdata are less accurate, the surface represented by the LiDAR observationscan be accurate if the over- or under-estimations are internally consis-tent (Hodgson et al., 2005). In this study, the surface was highlyRoughness (m)High: 3.95Low: ?2.16400 mN0Fig. 3. Roughness element heights (in m) following harvesting. The black lines separating the greener features (b0 m) and yellower features (>0 m) represent the 0-m contour line.Table 1Mean, median, maximum and minimum differences (in m) between LiDAR-derivedroughness element heights and measurements using a GPS and a measuring tape.# of points Mean Median Maximum MinimumSoil 8 0.026 0.010 0.086 0.000New stumps 5 0.187 0.115 0.385 0.000Old stumps 4 0.697 0.589 0.830 0.000Holes 4 0.085 0.080 0.180 0.000Residue 3 0.085 0.080 0.086 0.000Total 24 0.203 0.078 0.830 0.00005101510002000 1200FrequencyNeutralStableUnstablelog(z0 in m)0.122Fig. 4. Frequency distribution of log10(z0,a in m) for neutral, stable and unstable condi-tions. The black arrow indicates the mean value.230 E. Paul-Limoges et al. / Remote Sensing of Environment 136 (2013) 225?23389   consistent between the two LiDAR-derivedDTMs, despite the lower pointdensity for thepre-harvest DTM, and the under- or over-estimations of el-evations had to be similar to give such accurate results. This shows thatcomparisons can be made between different LiDAR datasets for thesame land area obtained at different timeswith sufficiently high accuracy.5.2. Feasibility of quantifying roughness element heights from differencesbetween two DTMsRoughness element heights derived from LiDAR were very accuratewhen compared to GPS and ground measurements. As low-height,near-ground features were studied, larger errors could have beenexpected due to the time response of LiDAR sensors. Previous studieshave found that the largest height errors for LiDAR were found whencharacterizing near-ground features due to the small height differencebetween the ground surface and roughness elements (e.g. Estornell etal., 2011; Meng et al., 2010; Streutker & Glenn, 2006). In this case,subtracting the pre-harvest DTM from the post-harvest DTM helpedin detecting more accurate roughness heights, as most roughness ele-ments were nonexistent before harvest.5.3. Feasibility of estimating aerodynamic roughness length from LiDARThis study showed also that standard deviation of roughness elementheightswas a good estimate of z0 of the clearcut surface. This is importantbecause it confirms that any method capable of representing the spatialroughness element height distribution could possibly be used to estimatez0. Other studies have also found good agreement between the variabilityin roughness element heights and z0 (Brown & Hugenholtz, 2011;Menenti & Ritchie, 1994). LiDAR can therefore be used to improve the ac-curacy of passive remote sensing formulations requiring z0 estimates inorder to correctly quantify surface-atmosphere heat and mass exchange.In addition, LiDAR-derived variations in roughness element heights canbe used for areas with no micrometeorological measurements or wherethe logarithmic wind profile does not apply due to limitations in fetchor complex topography. In this study, we were able to compare theLiDAR-derived roughness length to an independent estimate based on552620055264005526600552680055270005527200North (m)East (m)a) Neutral331800332000332200332400332600East (m)b) Unstable33180033200033220033240033260085%25%85%75%60%80%50%70%85%85%80%70%75%60%50%0 1log(z0)Fig. 5. The cumulative turbulent source area for a) neutral and b) unstable atmospheric stability conditions overlaying the map of log10(z0,l in m) estimated from the LiDAR analysis.Weighted z0,l is 0.13 m for both a) and b). The fraction of the turbulent source area outside the domain was 10.7% for a) and 12.2% for b). The yellow cross represents the micro-meteorological tower. The contour lines show the flux footprint areas corresponding to cumulative percentages of the probability of ?, i.e., the proportion of the flux measured bythe EC flux tower.Table 2Comparison of z0,a and z0,l for different filter sizes (F1 to F4) separated by wind direction. Selection criteria: wind speed > 1.5 m s?1, fraction outside the domain b15%, stability ?10 b ? b 0.01 (unstable and neutral) and successful determination of z0,a in the range 0.01 m b z0,a b 1 m. The sample count (n) represents the number of half-hour values used inthe aerodynamic calculations.Wind direction z0,a (m) z0,l [F1] (m) z0,l [F2] (m) z0,l [F3] (m) z0,l [F4] (m) nSector Median (25P?75P) Median (25P?75P) Median (25P?75P) Median (25P?75P) Median (25P?75P)NE (0??90?) 0.176 (0.07?0.39) 0.103 (0.10?0.10) 0.138 (0.13?0.14) 0.165 (0.16?0.17) 0.181 (0.18?0.18) 909SE (90??180?) 0.091 (0.05?0.17) 0.103 (0.10?0.10) 0.138 (0.14?0.14) 0.164 (0.16?0.17) 0.179 (0.18?0.18) 866SW (180??270?) 0.091 (0.04?0.18) 0.092 (0.09?0.09) 0.119 (0.12?0.12) 0.138 (0.13?0.14) 0.149 (0.14?0.15) 64NW (270??360?) 0.241 (0.12?0.39) 0.093 (0.09?0.09) 0.121 (0.12?0.12) 0.142 (0.14?0.14) 0.154 (0.15?0.16) 263All sectors 0.133 (0.06?0.27) 0.102 (0.10?0.10) 0.137 (0.13?0.14) 0.164 (0.16?0.17) 0.179 (0.17?0.18) 2102231E. Paul-Limoges et al. / Remote Sensing of Environment 136 (2013) 225?23390   the logarithmic wind profile equation, but this equation only applies tohomogeneous surfaces. Studies have found that only a small fraction ofwind profiles are suitable for the determination of z0 (e.g. Hatfield,1989; Kustas et al., 1989; Matthias et al., 1990). LiDAR is currently theonly technology able to estimate z0 over large, heterogeneous surfaces.5.4. Limitations and potential for future researchAlthough the current study was successful in estimating z0 using theequation adapted from Menenti and Ritchie (1994), this equation reliesonly on the variability in roughness elements and thus does not takeinto consideration other variables influencing z0. As shown in other stud-ies (e.g. Cheng & Castro, 2002; Jiang et al., 2008), the importance ofheight variability is likely to change depending on the spatial densityand distribution of the roughness elements and this is not reflected inthe formula applied in this paper. Furthermore, our study showed a con-siderable influence of the radius of the filter on the resulting estimate ofz0. A filter with a 2-m radius, which corresponded to about twice the di-mension of individual elements, was in best agreement with the micro-meteorological determination of z0. Yet, there is no objective set ofcriteria for choosing a particular filter radius. More studies need to bedone to examine the limitations of this equation and the influence of fil-ter size. In addition, for cases where most cell heights are similar to themean height of the cells, it is possible that the normalization, wherethe standard deviation is multiplied by the ratio of mean height to cellheight, is not needed as using the standard deviation alone would givesimilar results. In our case a subsample of 75,000 points suggested thatthe normalization had a relatively small effect (b5%).In addition, the source areamodel used an apriori estimate of z0 basedon initial calculations of z0 from measurements of u* and u. Another pos-sible approach would be to use an iterative process where z0 is adjustedbased on the LiDAR results for each half-hour calculation. However, thechange in the calculated z0 from slightly different source areas shouldbe very small when doing such iterations and the processing timewould be significantly increased. Also the source area model relies on ahomogeneous z0 upwind of themeasurement site, so the method is any-way not applicable to situations with substantial spatial changes in z0.Nevertheless, when validating LiDAR-derived z0 estimates, it is importantto compare the same area as for themicrometeorological measurements.Some previous studies have used micrometeorological measurementswithout appropriate source area calculations. The clearcut in our studywas surrounded by forest edges and is a good example of how surround-ing surfaces can influence the micrometeorological measurements. Bycomparing results under differentwind direction and stability conditions,we were able to discard measurements where the surrounding rough-ness features influenced our measurements. Without doing this, therewould have been an overestimation in the calculated z0.6. ConclusionsThis study showed that a high accuracy could be achieved whenusing LiDAR to generate DTMs, even under dense Douglas-fir canopies.It also demonstrated that subtractions can be performed using differentLiDAR datasets without significant residual errors. The standard devia-tion of roughness element height was found to give a good estimate ofz0 when compared to the results from the logarithmic wind profile cal-culation. These results suggest that LiDAR has the potential to be used tocharacterize z0 on larger spatial scales to improve surface flux estimatesmade with passive remote sensing technologies.AcknowledgmentsThis project was funded by the Natural Sciences and Engineering Re-search Council of Canada (NSERC) Discovery grants # 342029-07 (Chris-ten) and # 6123-07 (Black). Many thanks to our collaborators atTimberWest Forest Corp. for their assistance and use of their land. Weare grateful for assistance in the field and in the lab from Z. Nesic, N.J.Grant, T.D. Baker, and R. Ketler. E. Paul-Limoges was supported byNSERC through an Alexander Graham Bell Canada Postgraduate Scholar-ship and through the Canadian Meteorological and Oceanographic Soci-ety (CMOS) ? Weather Research House NSERC Scholarship supplement.ReferencesAdams, J., & Chandler, J. (2002). Evaluation of lidar and medium scale photogrammetryfor detecting soft-cliff coastal change. The Photogrammetric Record, 17, 405?418.Bastiaanssen, W. G. M., Menenti, M., Feddes, R. A., & Holtslag, A. A. M. (1998). A remotesensing surface energy balance algorithm for land (SEBAL). 1. Formulation. Journalof Hydrology, 212?213, 198?212.Brown, O. W., & Hugenholtz, C. H. (2011). Estimating aerodynamic roughness (zo) in mixedgrassland prairie with airborne LiDAR. Canadian Journal of Remote Sensing, 37, 422?428.Campbell, G. S., & Norman, J. M. (1998). An introduction to environmental biophysics(2nd ed.): Springer-Verlag New York, IncChen, B., Black, T. A., Coops, N. C., Hilker, T., Trofymow, J. A., & Morgenstern, K. (2009).Assessing tower flux footprint climatology and scaling between remotely sensedand eddy covariance measurements. Boundary-Layer Meteorology, 130, 137?167.Chen, J. M., Govind, A., Sonnentag, O., Zhang, Y., Barr, A., & Amiro, B. (2006). Leaf area indexmeasurements at Fluxnet-Canada forest sites. Agricultural and Forest Meteorology, 140,257?268.Cheng, H., & Castro, I. P. (2002). Near wall flow over urban-like roughness.Boundary-Layer Meteorology, 104, 229?259.Christen, A., Coops, N. C., Crawford, B. R., Kellett, R., Liss, K. N., Olchovski, I., et al. (2011).Validation of modeled carbon-dioxide emissions from an urban neighborhood withdirect eddy-covariance measurements. Atmospheric Environment, 45, 6057?6069.Table 3Comparison of z0,a and z0,l for different filter sizes (F1 to F4) separated by atmospheric stability. Selection criteria: wind speed > 1.5 m s?1, fraction outside the domain b30%, suc-cessful determination of z0,a in the range 0.01 m b z0,a b 1 m, and restricted to the wind directions between 70? and 280?. The sample count (n) represents the number of half-hourvalues used in the aerodynamic calculations.Stability range z0,a (m) z0,l [F1] (m) z0,l [F2] (m) z0,l [F3] (m) z0,l [F4] (m) nMedian (25P?75P) Median (25P?75P) Median (25P?75P) Median (25P?75P) Median (25P?75P)Extremely unstable(?10 b z/L b ?1)0.115 (0.05?0.25) 0.100 (0.09?0.10) 0.135 (0.12?0.14) 0.162 (0.15?0.17) 0.179 (0.17?0.18) 646Unstable(?1 b z/L b ?0.5)0.174 (0.08?0.38) 0.104 (0.10?0.10) 0.140 (0.14?0.14) 0.165 (0.16?0.17) 0.180 (0.18?0.18) 387Weakly unstable(?0.5 b z/L b ?0.2)0.119 (0.05?0.28) 0.104 (0.10?0.10) 0.138 (0.14?0.14) 0.164 (0.16?0.17) 0.179 (0.18?0.18) 356Near-neutral(?0.2 b z/L b ?0.1)0.083 (0.05?0.15) 0.104 (0.10?0.10) 0.139 (0.14?0.14) 0.165 (0.16?0.17) 0.181 (0.18?0.18) 202Neutral(?0.1 b z/L b 0)0.107 (0.06?0.19) 0.101 (0.10?0.10) 0.136 (0.13?0.14) 0.163 (0.15?0.17) 0.179 (0.17?0.18) 419Table 4Mean and standard deviation for estimated z0.l using the F2 filter and for calculated z0,aunder neutral and all stability conditions.Mean z0,l ? SD (m) Mean z0,a ? SD (m)Neutral stability conditions 0.13 ? 0.01 0.13 ? 0.13All stability conditions 0.13 ? 0.01 0.12 ? 0.30232 E. Paul-Limoges et al. / Remote Sensing of Environment 136 (2013) 225?23391   Colin, J., & Faivre, R. (2010). Aerodynamic roughness length estimation from veryhigh-resolution imaging LIDAR observations over the Heihe basin in China. Hydrol-ogy and Earth System Sciences, 14, 2661?2669.Colin, J., Menenti, A., Rubio, E., & Jochum, A. (2006a). A multi-scales surface energy bal-ance system for operational actual surface evapotranspiration monitoring. In G.Durso, M. A. O. Jochum, & J. Moreno (Eds.), Earth observation for vegetation monitor-ing and water management (pp. 178?184).Colin, J., Menenti, M., Rubio, E., & Jochum, A. (2006b). Accuracy vs. operability: A casestudy over barrax in the context of the DEMETER project. In G. Durso, M. A. O.Jochum, & J. Moreno (Eds.), Earth observation for vegetation monitoring and watermanagement (pp. 75?83).Crago, R., Okello, W., & Jasinski, M. (2012). Equations for the drag force and aerodynamicroughness length of urban areas with randombuilding heights. Boundary-Layer Mete-orology, 145, 423?437.DeVries, A. C., Kustas,W. P., Ritchie, J. C., Klaassen,W.,Menenti,M., Rango, A., et al. (2003).Effective aerodynamic roughness estimated from airborne laser altimeter measure-ments of surface features. International Journal of Remote Sensing, 24, 1545?1558.Estornell, J., Ruiz, L. A., & Velazquez-Marti, B. (2011). Study of shrub cover and heightusing LIDAR data in a Mediterranean area. Forest Science, 57, 171?179.Ferster, C. J., Trofymow, J. A., Coops, N. C., Chen, B. Z., Black, T. A., & Gougeon, F. A. (2011).Determination of ecosystem carbon-stock distributions in the flux footprint of aneddy-covariance tower in a coastal forest in British Columbia. Canadian Journal of For-est Research-Revue Canadienne De Recherche Forestiere, 41, 1380?1393.Garratt, J. R. (1992). The atmospheric boundary layer. Cambridge: Cambridge UniversityPress (316 pp).Gomes Pereira, L. M., & Janssen, L. L. F. (1999). Suitability of laser data for DTM gener-ation: A case study in the context of road planning and design. ISPRS Journal of Pho-togrammetry and Remote Sensing, 54, 244?253.Grimmond, C. S. B., & Oke, T. R. (1999). Aerodynamic properties of urban areas derivedfrom analysis of surface form. Journal of Applied Meteorology, 38, 1262?1292.Hatfield, J. L. (1989). Aerodynamic properties of partial canopies. Agricultural and ForestMeteorology, 46, 15?22.Hilker, T., Hall, F. G., Coops, N. C., Lyapustin, A., Wang, Y., Nesic, Z., et al. (2010). Remotesensing of photosynthetic light-use efficiency across two forested biomes: Spatialscaling. Remote Sensing of Environment, 114, 2863?2874.Hodgson, M. E., Jensen, J., Raber, G., Tullis, J., Davis, B. A., Thompson, G., et al. (2005). Anevaluation of lidar-derived elevation and terrain slope in leaf-off conditions. Photo-grammetric Engineering and Remote Sensing, 71, 817?823.Humphreys, E. R., Black, T. A., Morgenstern, K., Cai, T., Drewitt, G. B., Nesic, Z., et al. (2006).Carbon dioxide fluxes in coastal Douglas-fir stands at different stages of developmentafter clearcut harvesting. Agricultural and Forest Meteorology, 140, 6?22.Jiang, D. H., Jiang, W. M., Liu, H. N., & Sun, J. N. (2008). Systematic influence of differentbuilding spacing, height and layout on mean wind and turbulent characteristicswithin and over urban building arrays. Wind and Structures, 11, 275?289.Kaimal, J. C., & Finnigan, J. J. (1994). Atmospheric boundary layer flows: Their structureand measurements. New York: Oxford University Press.Kormann, R., & Meixner, F. (2001). An analytical footprint model for non-neutral strat-ification. Boundary-Layer Meteorology, 99, 207?224.Kraus, K., & Mikhail, E. M. (1972). Linear least-squares interpolation. PhotogrammetricEngineering, 38, 1016?1029.Kraus, K., & Pfeifer, N. (1998). Determination of terrain models in wooded areas withairborne laser scanner data. ISPRS Journal of Photogrammetry and Remote Sensing,53, 193?203.Kustas, W. P., Blanford, J. H., Stannard, D. I., Daughtry, C. S. T., Nichols, W. D., &Weltz, M.A. (1994). Local energy flux estimates for unstable conditions using variance datain semiarid rangelands. Water Resources Research, 30, 1351?1361.Kustas, W. P., Choudhury, B. J., Kunkel, K. E., & Gay, L. W. (1989). Estimate of the aero-dynamic roughness parameters over an incomplete canopy cover of cotton. Agri-cultural and Forest Meteorology, 46, 91?105.Lettau, H. (1969). Note on aerodynamic roughness-parameter estimation on the basisof roughness-element description. Journal of Applied Meteorology, 8, 828?832.Matthias, A. D., Kustas, W. P., Gay, L. W., Cooper, D. I., Alves, L. M., & Pinter, P. J., Jr.(1990). Aerodynamic parameters for a sparsely roughened surface composed ofsmall cotton plants and ridged soil. Remote Sensing of Environment, 32, 143?153.McGaughey, R. J. (2010). FUSION/LDV: Software for LIDAR data analysis and visualization.United States Department of Agriculture, Forest Service.Menenti, M., & Ritchie, J. C. (1994). Estimation of effective aerodynamic roughness ofwalnut gulch watershed with laser altimeter measurements. Water Resources Re-search, 30, 1329?1337.Meng, X., Currit, N., & Zhao, K. (2010). Ground filtering algorithms for airborne LiDARdata: A review of critical issues. Remote Sensing, 2, 833?860.Monteith, J. L., & Unsworth, M. H. (2008). Principles of environmental physics. Oxford:Academic Press.Morgenstern, K., Black, T. A., Humphreys, E. R., Griffis, T. J., Drewitt, G. B., Cai, T., et al.(2004). Sensitivity and uncertainty of the carbon balance of a Pacific NorthwestDouglas-fir forest during an El Ni?o/La Ni?a cycle. Agricultural and Forest Meteorol-ogy, 123, 201?219.Pfeifer, N., K?stli, A., & Kraus, K. (1998). Interpolation and filtering of laser scanner data ?Implementation and first results. International archives of photogrammetry and remotesensing. Columbus, OH, USA.Pojar, J., Klinka, K., & Demarchi, D. A. (1991). Coastal western hemlock zone. In D.Meidinger, & J. Pojar (Eds.), Ecosystems of British Columbia. BC special report series,No 6. (pp. 95?111). Victoria: BC Ministry of Forests.Raupach, M. R. (1992). Drag and drag partition on rough surfaces. Boundary-Layer Me-teorology, 60, 375?395.Raupach, M. R. (1994). Simplified expressions for vegetation roughness length andzero-plane displacement as functions of canopy height and area index. Boundary-LayerMeteorology, 71, 211?216.Reutebuch, S. E., McGaughey, R. J., Andersen, H. E., & Carson, W. W. (2003). Accuracy ofa high-resolution lidar terrain model under a conifer forest canopy. Canadian Jour-nal of Remote Sensing, 29, 527?535.Roerink, G. J., Su, Z., & Menenti, M. (2000). S-SEBI: A simple remote sensing algorithmto estimate the surface energy balance. Physics and Chemistry of the Earth, Part B:Hydrology, Oceans and Atmosphere, 25, 147?157.Shaw, R. H., & Pereira, A. R. (1982). Aerodynamic roughness of a plant canopy ? A nu-merical experiment. Agricultural Meteorology, 26, 51?65.Streutker, D. R., & Glenn, N. F. (2006). LiDAR measurement of sagebrush steppe vege-tation heights. Remote Sensing of Environment, 102, 135?145.Stull, R. (1988). An introduction to boundary-layer meteorology. The Netherlands: KluwerAcademic Publishers.Su, Z. (2002). The Surface Energy Balance System (SEBS) for estimation of turbulentheat fluxes. Hydrology and Earth System Sciences, 6, 85?99.Tian, X., Li, Z. Y., van der Tol, C., Su, Z., Li, X., He, Q. S., et al. (2011). Estimating zero-plane dis-placement height and aerodynamic roughness length using synthesis of LiDAR andSPOT-5 data. Remote Sensing of Environment, 115, 2330?2341.Turner, D. P., Koerper, G. J., Harmon, M. E., & Lee, J. J. (1995). A carbon budget for forestsof the conterminous United States. Ecological Applications, 5, 421?436.Xie, Z. -T., Coceal, O., & Castro, I. (2008). Large-eddy simulation of flows over randomurban-like obstacles. Boundary-Layer Meteorology, 129, 1?23.Yasuda, N. (1988). Turbulent diffusivity and diurnal variations in the atmosphericboundary layer. Boundary-Layer Meteorology, 43, 209?221.233E. Paul-Limoges et al. / Remote Sensing of Environment 136 (2013) 225?233

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-0074280/manifest

Comment

Related Items