Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Towards an understanding of coupled carbon, water and nitrogen dynamics at sand, landscape and regional… Chen, Baozhang 2012

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

Item Metadata

Download

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

Full Text

 TOWARDS AN UNDERSTANDING OF COUPLED CARBON, WATER AND NITROGEN DYNAMICS AT STAND, LANDSCAPE AND REGIONAL SCALES  by Baozhang Chen A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Forestry) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2012    © Baozhang Chen, 2012  ii Abstract One of the critical issues in the prognoses of future climate change is a comprehensive understanding of the global carbon budget. Progress in C balance studies has been achieved either at stand or at continental scales. However, the coupled terrestrial carbon, nitrogen and hydrological dynamics are yet far from well understood and methods to estimate the land- atmosphere carbon fluxes at the landscape and regional scales are notably lacking. The major findings of this dissertation research are as follows: First, this dissertation improves our understanding of the terrestrial C processes, for example, at Douglas-fir stand in the Pacific Northwest: (i) Although the majority of carbon sequestration occurred during March through June, May through August was responsible for about 80% of the inter-annual variability of net ecosystem productivity (NEP). The major drivers of inter-annual variability of annual carbon fluxes were annual and spring mean temperatures (Ta) and water deficiency during late summer to autumn; (ii) Monthly GPP was linearly correlated with photosynthetic active radiation (Q) (r 2  = 0.85) and monthly Re was exponentially correlated with Ta (r 2  = 0.94); (iii) The responses of NEP to changes of Ta and Q were positive during the first and last four months of the year but were negative during the middle four months of the year. (iv) N fertilization increased annual NEP by ~83%, in the first year, resulted from increases in annual GPP by ~8% and from decreases in annual Re by ~5.8%. Secondly, this dissertation develops a pragmatic algorithm with synergy of footprint climatology and geospatial analyses for assessing the spatial representativeness of eddy- covariance flux tower measurements. This algorithm was then applied to the Canadian Carbon Program network. Thirdly, this dissertation develops an innovative up-scaling strategy by integrating ecosystem modeling, footprint climatology modeling, remote sensing, and data-model fusion for the scaling of C fluxes at stand, landscape and regional scales.  iii And fourth, this dissertation develops an analytical scalar concentration footprint model to assess the influences of land surface heterogeneity on tower CO2 concentration measurements.  Summarily, this dissertation research provides a sound basis for shaping future climate change adaptation policy related to carbon management.     iv Preface In this thesis, Chapters 2 through 8 are derived from multi-authored published or in-pressed or submitted peer-reviewed journal papers. The author of this dissertation, Baozhang Chen, is the lead author on all the seven scientific publications. He is the primary person responsible for providing the hypotheses, designing research framework, developing models and methodologies, coordinating data providers, processing data, analyzing results, drawing conclusions and writing and revising all aforementioned works. The author was also involved in data collection for several of the projects. Co-authors, as noted in the footnotes to the relevant Chapters and in the references, were involved with collecting data at many of the field sites, producing ancillary data, and developing the observation frameworks. Colleagues and anonymous reviewers also provided numerous comments during revision of these manuscripts, as noted in the acknowledgements. Copyright permission for reproducing published works has been granted and mentioned in footnotes to the relevant sections with permission letters filed with the thesis office of University of British Columbia. A version of Chapter 2 has been published. [B. Chen] and N.C. Coops (2009), Understanding of coupled terrestrial carbon, nitrogen and water dynamics — an overview. Sensors 9(11), 8624-8657. Impact Factor: 1.771. I reviewed and summarized published research results and conducted all research pertaining to this publication and wrote the manuscript, while my supervisor N.C. Coops helped with some constructive scientific comments and editorial polishing. A version of Chapter 3 has been published. [B. Chen], T.A. Black, N.C. Coops, P. Krishna, R. Jassal, C. Brȕmmer and Z. Nesic (2009), Seasonal controls on interannual variability in carbon dioxide exchange of a Pacific Northwest Douglas-fir forest, 1997 – 2006. Global Change Biology, 15(8), 1962-1981. Impact Factor: 6.346. I performed all the research, data analyses, and interpretation of the results, and wrote the manuscript and finished the final revision. My co-supervisor, T.A. Black leads the data collections and data QC/QA as the principal investigator of this eddy-covariance tower and helped with the direction and interpretation of results; my supervisor, N.C. Coops, provided advices on data analyses  v methodology; all other co-authors helped the long-term (9-years) eddy covariance data collections and data QC/QA. A version of Chapter 4 has been published. [B. Chen], N.C. Coops, T.A. Black, R.S. Jassal, J.M. Chen and M. Johnson (2011), Modeling to discern nitrogen fertilization impacts on carbon sequestration in a Pacific Northwest Douglas-fir forest in the first-post fertilization year. Global Change Biology, 17 (3), 1442–1460. Impact Factor: 6.346. I designed the modeling framework, developed data-model fusion algorithm, conducted model simulations, interpreted and analyzed the results, and wrote the manuscript and finished the final revision. My co-supervisor, T.A. Black helped with the explanation of modeling results; my supervisor, N.C. Coops, provided advice on “Discussion” section; R.S. Jassal provided leaf nitrogen measurements and soil efflux data; all other co-authors made editorial comments. A version of Chapter 5 has been published. [B. Chen], T.A. Black, N.C. Coops, T. Hilker, T. Trofymow, Z. Nesic and K. Morgenstern (2009), Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements. Boundary-Layer Meteorology, 130(2), 137-167. Impact Factor: 1.879. Based on K-theory and an analytical footprint solution of the two dimensional Eulerian advection–diffusion equation by van Ulden (1978) and Horst and Weil (1994), following parameterization methods in Kormann and Meixner (2001), I developed the flux footprint climatology model (SAFE-F: Simple Analytical Footprint model based on Eulerian coordinates for flux). I estimated footprint climatology, interpreted and analyzed the results, and wrote the manuscript and finished the final revision. My co-supervisor, T.A. Black leads the data collections and data QC/QA as the principal investigator of these eddy-covariance towers. My supervisors (T.A. Black and N.C. Coops) helped with the explanation of results and made editorial comments and other co-authors provided ancillary data. A version of Chapter 6 has been published. [B. Chen], N.C. Coops, D. Fu, H.A. Margolis, B.D. Amiro, A.G. Barr, T. A. Black, M.A. Arain, C. P.-A. Bourque, L.B. Flanagan, P.M. Lafleur, J. H. McCaughey and S.C. Wofsy (2011), Assessing eddy-covariance flux tower location bias across the Fluxnet-Canada Research Network based on remote sensing and footprint modeling. Agriculture and Forest Meteorology, 150(1), 87-100. Impact Factor:  vi 3.228. I designed the modeling framework, ran the footprint model, conducted data analyses, interpreted the results, and wrote the manuscript and made the final revisions. My supervisor, N.C. Coops helped with the direction.  Dr. Hank A. Margolis designed the Fluxnet-Canada Research Network as the network leader and other co-authors provided eddy-covariance flux data, metrological data and ancillary data and some of them made editorial comments. A version of Chapter 7 is in press. [B. Chen], N.C. Coops, D. Fu, H.A. Margolis, B.D. Amiro, T.A. Black, M.A. Arain, A.G. Barr, C.P.-A. Bourque, L.B. Flanagan, P. M. Lafleur, J.H. McCaughey, S.C. Wofsy (2011), Characterizing spatial representativeness of flux tower eddy-covariance measurements across the Canadian Carbon Program network using remote sensing and footprint analysis. Remote Sensing of Environment, 156 (in press). Impact Factor: 3.951. I designed the modeling framework, I performed main parts of  the research, data analyses, interpreted the results, and wrote the manuscript and finished the final revision. My supervisor, N.C. Coops provided advices on geospatial analyses methods and other co-authors provided ancillary data and some of them made editorial comments. A version of Chapter 8 has been submitted to a peer-reviewed journal.  [B. Chen], N.C. Coops, D.E.J. Worthy, T.A. Black (submitted), Assessing scalar concentration footprint climatologies and land surface impacts on tower CO2 measurements in the boreal forest of central Saskatchewan, Canada. I conducted all research pertaining to this publication, while other co-authors provided ancillary data and made editorial comments.  vii Table of Contents Abstract .................................................................................................................................... ii Preface ..................................................................................................................................... iv Table of Contents .................................................................................................................. vii List of Tables ........................................................................................................................ xiv List of Figures ....................................................................................................................... xvi List of Abbreviations, Acronyms, and Symbols* ............................................................ xxix Acknowledgements .......................................................................................................... xxxiv Dedication ......................................................................................................................... xxxvi Chapter  1: Introduction ........................................................................................................ 1 1.1 Scientific background ............................................................................................... 1 1.1.1 Terrestrial CO2 cycle............................................................................................. 1 1.1.2 Terrestrial water cycle ........................................................................................... 5 1.1.3 Terrestrial N cycle................................................................................................. 8 1.1.4 Coupling of the C, N and water cycles ................................................................. 8 1.2 State of approaches for studying coupled terrestrial C and water dynamics ............ 9 1.2.1 State of field-scale monitoring and modeling approaches .................................... 9 1.2.2 Arrays of airborne and satellite sensors developed to monitor coupled C and water cycles ................................................................................................................... 10 1.3 Hypothesis and objectives....................................................................................... 12 1.3.1 Knowledge and scale gaps .................................................................................. 12 1.3.2 Hypotheses .......................................................................................................... 13 1.3.3 Objectives ........................................................................................................... 13 1.4 Study areas and data descriptions ........................................................................... 14 1.4.1 Stand-level research ............................................................................................ 14 1.4.2 Landscape-level research .................................................................................... 17 1.4.3 Regional C flux research ..................................................................................... 18 1.5 Thesis structure ....................................................................................................... 21 Chapter  2: Understanding of coupled terrestrial carbon, nitrogen and water dynamics—an overview ........................................................................................................ 23 2.1 Synopsis .................................................................................................................. 23 2.2 Introduction ............................................................................................................. 23 2.3 Scientific background and state of knowledge ....................................................... 25 2.3.1 Overview of terrestrial ecosystems and climate ................................................. 25 2.3.2 Terrestrial C cycling ........................................................................................... 26 2.3.3 Coupling of the C and water cycles .................................................................... 31  viii 2.3.4 Coupling of the C and N cycles .......................................................................... 31 2.4 Monitoring of C and water cycling in terrestrial ecosystems ................................. 33 2.4.1 Global flux tower network (FLUXNET) ............................................................ 33 2.4.2 CO2 concentration measurements, data assimilation and CarbonTracker .......... 34 2.4.3 Stable C isotope measurements .......................................................................... 35 2.4.4 Satellite monitoring ............................................................................................. 36 2.4.5 Other airborne measurements ............................................................................. 37 2.5 Modeling of C and water dynamics in terrestrial ecosystems based on remote sensing ……………………………………………………………………………………..38 2.5.1 Land surface and ecosystem modeling ............................................................... 39 2.5.2 Spatially-distributed hydrological processes modeling ...................................... 41 2.5.3 Modeling dynamics of stable C isotopic exchange between ecosystem and the atmosphere .................................................................................................................... 41 2.5.4 Modeling coupled C and water dynamics --- an eco-hydrological approach ..... 42 2.5.5 Applications of remotely-sensed data in ecohydrological modeling .................. 44 2.6 Research gaps in C and water flux estimates and scaling approaches .................... 46 2.7 Future research directions ....................................................................................... 49 2.7.1 Development of a spatially explicit ecohydrological modeling framework ....... 49 2.7.1.1 Reviewing of the existing EASS model (Chen et al., 2007a) ..................... 50 2.7.1.2 Reviewing of the TerrainLab model (Chen et al., 2005) ............................ 51 2.7.1.3 Development of an ecohydrological model by coupling of EASS with TerrainLab................................................................................................................... 52 2.7.1.4 Model calibration and validation for EASS-TerrainLab ............................ 53 2.7.1.5 Model sensitivity analysis and runs under different scenarios ................... 53 2.7.2 Landscape and regional C and water fluxes estimation: an up-scaling framework  ………………………………………………………………………………….53 2.8 Summary ................................................................................................................. 54 Chapter  3: Seasonal controls on interannual variability in C dioxide exchange of a near end-of rotation Douglas-fir stand in the Pacific Northwest, 1997–2006 .......................... 56 3.1 Synopsis .................................................................................................................. 56 3.2 Introduction ............................................................................................................. 57 3.3 Methods................................................................................................................... 59 3.3.1 Site description.................................................................................................... 59 3.3.2 Instrumentation ................................................................................................... 60 3.3.3 Data processing and analysis .............................................................................. 63 3.3.3.1 Data processing and quality controls .......................................................... 63  ix 3.3.3.2 Gap-filling, time integration and uncertainties ........................................... 64 3.3.3.3 Data analysis methods................................................................................. 67 3.4 Results and discussion ............................................................................................ 67 3.4.1 General weather conditions................................................................................. 67 3.4.2 Diurnal and seasonal CO2 exchange ................................................................... 72 3.4.3 Climatic controls on seasonal CO2 exchange ..................................................... 78 3.4.3.1 Response of CO2 exchange to Q and Ta...................................................... 79 3.4.3.2 Response of CO2 exchange to water stress ................................................. 84 3.4.4 Dependence of the annual NEP on GPP and Re .................................................. 86 3.4.5 Climatic controls on annual and inter-annual CO2 exchange ............................. 89 3.5 Summary and conclusions ...................................................................................... 91 Chapter  4: Modeling to discern nitrogen fertilization impacts on carbon sequestration in a Pacific Northwest Douglas-fir forest in the first-post fertilization year ................... 93 4.1 Synopsis .................................................................................................................. 93 4.2 Introduction ............................................................................................................. 94 4.3 Material and methods .............................................................................................. 96 4.3.1 Site description.................................................................................................... 96 4.3.2 Stand fertilization ................................................................................................ 97 4.3.3 Needle mass and N concentration and growth increment measurements ........... 98 4.3.4 Climate and EC measurements ........................................................................... 98 4.3.5 Soil CO2 efflux data used in this study ............................................................... 99 4.3.6 Modeling approach ............................................................................................. 99 4.3.6.1 Model description ....................................................................................... 99 4.3.6.2 Model experimental design ....................................................................... 100 4.3.6.3 Model parameter optimization algorithm ................................................. 101 4.3.6.4 Parameter selection and ensemble generation .......................................... 101 4.3.7 Data analysis methods....................................................................................... 102 4.4 Results ................................................................................................................... 103 4.4.1 Comparison of measured environmental variables between pre- and post- fertilization year .......................................................................................................... 103 4.4.2 Comparison of measured C component fluxes between pre- and post- fertilization year .......................................................................................................... 106 4.4.3 Response of needle mass, N concentration and growth increment to fertilization  …………………………………………………………………………………107  x 4.4.4 Discerning N fertilization effects on C component fluxes using the optimized BEPS model ................................................................................................................ 107 4.4.4.1 Calibration and verification of the optimized BEPS model...................... 108 4.4.4.2 Discerning N fertilization impacts ............................................................ 114 4.5 Discussion ............................................................................................................. 117 4.5.1 Forest fertilization effects on photosynthesis ................................................... 117 4.5.2 Forest fertilization effects on respiration .......................................................... 119 4.5.3 N-use efficiency for C sequestration ................................................................. 122 4.5.4 Uncertainty and limitation of this modeling approach ..................................... 123 4.6 Conclusions ........................................................................................................... 124 Chapter  5: Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements ...................................................................... 126 5.1 Synopsis ................................................................................................................ 126 5.2 Introduction ........................................................................................................... 127 5.3 Materials and methods .......................................................................................... 131 5.3.1 Research area and site descriptions................................................................... 131 5.3.2 EC flux and meteorological measurements ...................................................... 133 5.3.3 Land surface remote sensing measurements ..................................................... 135 5.3.4 Footprint modelling .......................................................................................... 136 5.3.4.1 Half-hourly footprint modelling ............................................................... 136 5.3.4.2 Footprint climatology estimates ................................................................ 139 5.3.5 Modelling spatially-distributed GPP using land surface remote sensing data .. 141 5.3.6 Scaling of remotely sensed GPP to EC-flux derived GPP ................................ 142 5.3.6.1 Overlaying pure footprint map on remotely-sensed GPP field ................. 142 5.3.6.2 Sensor location bias .................................................................................. 142 5.4 Results ................................................................................................................... 143 5.4.1 Wind roses ........................................................................................................ 143 5.4.2 Seasonal variations in C fluxes ......................................................................... 145 5.4.3 Half-hourly footprints ....................................................................................... 146 5.4.4 Footprint climatology........................................................................................ 150 5.4.4.1 Variations in seasonal footprint climatology ............................................ 150 5.4.4.2 Annual footprint climatology and its interannual variations .................... 151 5.4.5 Overlaying footprint climatology on the remotely sensed GPP field ............... 158  xi 5.4.6 Comparing estimated values of remotely sensed and EC derived GPP ............ 159 5.5 Discussion ............................................................................................................. 160 5.5.1 Advantages and disadvantages using an analytical footprint model for assessing footprint climatology ................................................................................................... 160 5.5.2 Uncertainty in long-term EC flux caused by footprint climatology variations 163 5.5.3 Spatial scaling using footprint climatology ...................................................... 164 5.6 Conclusions and implications ............................................................................... 166 Chapter  6: Assessing eddy-covariance flux tower location bias across the Fluxnet- Canada Research Network based on remote sensing and footprint modeling.............. 167 6.1 Synopsis ................................................................................................................ 167 6.2 Introduction ........................................................................................................... 168 6.3 Site characteristics, tower measured data and Landsat images ............................. 170 6.4 Theory and methodology ...................................................................................... 175 6.4.1 Point-to-area representativeness ....................................................................... 175 6.4.2 Footprint and footprint climatology modelling................................................. 176 6.4.2.1 Footprint PDF modelling .......................................................................... 176 6.4.2.2 Footprint climatology PDF modelling ...................................................... 177 6.4.3 Analysis of Landsat images .............................................................................. 178 6.4.3.1 Normalized Difference Vegetation Index (NDVI) ................................... 178 6.4.3.2 Enhanced Vegetation Index (EVI) ............................................................ 179 6.4.3.3 Sensor location bias (SLB) ....................................................................... 179 6.5 Results ................................................................................................................... 180 6.5.1 Footprint climatology........................................................................................ 180 6.5.2 Sensor location bias .......................................................................................... 185 6.6 Discussion ............................................................................................................. 192 6.7 Summary and conclusions .................................................................................... 195 Chapter  7: Characterizing spatial representativeness of flux tower eddy-covariance measurements across the Canadian Carbon Program network using remote sensing and footprint analysis ................................................................................................................. 197 7.1 Synopsis ................................................................................................................ 197 7.2 Introduction ........................................................................................................... 198 7.3 Methods and data .................................................................................................. 199 7.3.1 Site characteristics ............................................................................................ 199 7.3.2 Footprint and footprint climatology modelling................................................. 200 7.3.3 Digital elevation model (DEM) data................................................................. 201 7.3.4 Land cover classification and analysis .............................................................. 201  xii 7.3.5 Analysis of Landsat images .............................................................................. 203 7.3.6 Semivariogram .................................................................................................. 203 7.3.7 Window size analysis ........................................................................................ 205 7.4 Results ................................................................................................................... 205 7.4.1 Footprint climatology........................................................................................ 205 7.4.2 Topography ....................................................................................................... 208 7.4.3 Land cover ........................................................................................................ 210 7.4.4 Semivariograms ................................................................................................ 215 7.4.5 Window size analysis ........................................................................................ 219 7.5 Discussion ............................................................................................................. 220 7.5.1 Limitations and uncertainties ............................................................................ 220 7.5.2 Implications....................................................................................................... 223 7.5.3 Applications and future directions .................................................................... 224 7.6 Conclusions ........................................................................................................... 225 Chapter  8: Assessing scalar concentration footprint climatologies and land surface impacts on tower CO2 concentration measurements in the boreal forest of central Saskatchewan, Canada ....................................................................................................... 227 8.1 Synopsis ................................................................................................................ 227 8.2 Introduction ........................................................................................................... 228 8.3 Materials ............................................................................................................... 230 8.3.1 Study site measurements ................................................................................... 230 8.3.2 Study site descriptions ...................................................................................... 230 8.3.3 Land cover information over the tower concentration footprint area ............... 231 8.4 Methodology ......................................................................................................... 234 8.4.1 Daily scalar concentration footprint modeling ................................................. 234 8.4.2 Footprint climatology determination ................................................................ 238 8.4.3 Land cover classification and analysis .............................................................. 239 8.5 Results ................................................................................................................... 241 8.5.1 Daily concentration footprints .......................................................................... 241 8.5.2 Variations in seasonal concentration footprint climatologies ........................... 242 8.5.3 Annual concentration footprint climatologies .................................................. 246 8.5.4 Land surface impacts on tower CO2 measurements ......................................... 248 8.6 Discussion ............................................................................................................. 253 8.7 Conclusions and implications ............................................................................... 255 Chapter  9: Conclusions of thesis ...................................................................................... 257 9.1 Summary of findings............................................................................................. 257 9.2 Future work ........................................................................................................... 260 9.2.1 Development of a spatially explicit ecohydrological modeling framework ..... 260  xiii 9.2.2 Application of the algorithm for assessing the spatial representativeness of flux tower measurements in the global network of EC towers (FLUXNET) ..................... 261 9.2.3 Reducing uncertainties in regional C flux estimates by mutual constraint of “top- down” and “bottom-up” approaches based on tall-tower CO2 concentration measurements and remote sensing .............................................................................. 261 References ............................................................................................................................ 263 Appendix A: Photosynthesis and respiration calculations in the BEPS –EASS model 304 Appendix B: Sensitivity of stomatal conductance to soil water variability ................... 307   xiv List of Tables Table 1.1 Major stocks of water on the Earth .......................................................................... 6 Table 1.2 Major components of annual terrestrial hydrological fluxes ................................... 7 Table 1.3 Characteristics and tower measurement information for the CCP sites a)  .............. 20 Table 3.1 Climate statistics for the research site .................................................................... 69 Table 3.2 Climate statistics for a long-term climate station near the research site *  ............... 72 Table 3.3 Seasonal and interannual variability of C fluxes and associated climate variables expressed as the coefficient of variation (CV), 1998 - 2006* ................................................ 78 Table 3.4 Characteristics of linear regression analysis [y = ax + b] of monthly mean ecosystem respiration (Re, in μmol m -2  s -1) and gross primary productivity (GPP, in μmol m-2 s -1 ) vs. monthly mean air temperature (Ta, in ºC), and daytime mean photosynthetically active radiation (Q, in μmol m-2 s-1) for individual month and annual clusters. Data are from Nov 1997 to Dec 2006* .................................................................................................................. 80 Table 3.5 Characteristics of linear regression analysis [y = ax + b] of monthly mean net ecosystem productivity (NEP, in μmol m-2 s-1) vs. monthly mean air temperature (Ta, in ºC) and daytime mean photosynthetically active radiation (Q, in μmol m-2 s-1) for each 4-month period. Data are from Nov 1997 to Dec 2006* ....................................................................... 83 Table 3.6 Annual sums of net ecosystem productivity (NEP), gross primary productivity (GPP), and ecosystem respiration (Re) in g C m -2  and number of C sink days (CSD) during the year and the dry-season months (Apr - Sep) when the daily GPP to Re ratio was greater than 1 ....................................................................................................................................... 87 Table 4.1 Model parameters optimized for a 58-year-old Douglas-fir stand ....................... 102 Table 4.2 Comparison of measured annual C component fluxes between non-fertilization year (1998-2006 and the fertilization year (2007)  *  .............................................................. 106  xv Table 4.3 Comparisons of measured and modeled annual gross primary productivity (GPP), net ecosystem productivity (NEP), ecosystem respiration (Re) and  soil respiration (Rs) * . 112 Table 5.1 Selected stand, topography, and EC system characteristics for the three sites .... 133 Table 5.2 Monthly median half-hourly meteorological values (July, 2004) during daytime and nighttime for the three sites, selected to illustrate typical half-hourly flux footprints under unstable and stable conditions* ............................................................................................ 146 Table 5.3 Comparison of selected mean annual cumulative footprints weighted by NEP, Re and GPP and pure (unweighted) footprint with their corresponding areas (and radii) for the three sites. ............................................................................................................................. 153 Table 6.1 List of abbreviations, key attributes of vegetation at FCRN sites*...................... 172 Table 6.2 List of abbreviations, site characteristics and tower measurement information for FCRN sites a)  ......................................................................................................................... 173 Table 6.3 Acquired Landsat imagery information for FCRN sites used in this study ......... 175 Table 6.4 Comparison of annual cumulative footprints with their corresponding areas and SLB (δ, in %) for all the 12 sites........................................................................................... 183 Table 7.1 Land cover classes of EOSD and aggregation schemes....................................... 202 Table 7.2 Percentages of observed dominant land cover type in a reference footprint area: 50%, 80%, 90% and 99% annual footprint climatology areas and the circular areas centred at each tower location with radius of 1, 2, and 3 km ................................................................ 214 Table 8.1 Land cover type classification of the GLC2000 products and two regrouping (8 classes and 4 classes, respectively) schemes ........................................................................ 233 Table 8.2 Parameters for characterizing atmospheric stability for two arbitrarily selected days (a sunny and cloudy days in summer) and for describing daily footprint features for measurements at four heights as shown in Figure 8.1 * ....................................................... 243  xvi List of Figures Figure 1.1 Study areas.  A, Stand scale: (a) pictures for the DF49 EC flux tower site of British Columbia flux station, with Douglas-fir dominated forest. In 2010 the stand ages adjacent to the tower was approximately 62 years old.  DF49 was harvested in December 2010; B, Landscape scale: (b) a land cover map of a 6 × 6 km area centered at the DF40 EC tower at a 30-m resolution, the white line is the 90% annual footprint climatology contour at the DF49 tower in 2004; (c) British Columbia flux station, which includes three EC towers (DF49, HDF88, and HDF00) with different ages of Douglas-fir dominated forest and in 2010 the stand ages adjacent to each tower were approximately 62, 22, and 11 years old, respectively; (d) the eco-region map of Canada’s landmass and with the locations of high- precision CO2 concentration towers in the long term observations network for atmospheric measurements of CO2; C, Regional scale: (e) a 1000 × 1000 km region surrounding two high-precision CO2 concentration towers: East Trout Lake (ETL) 106-m tall tower, Saskatchewan, located at (54°21'N, 104°59'W) and  Candle Lake 28-m high tower (CL), located at (53°59'N, 105°07'W); while triangles are the cluster of EC towers within the tower concentration footprint area. Land cover data at 1 km resolution are from the Global Land Cover Map 2000 (GLC2000); (f) picture of the ETL 106-m tall tower with CO2 concentration measurements at four levels. ................................................................................................... 16 Figure 1.2 Locations of flux towers in the Fluxnet Canada Research Network (FCRN) / Canadian Carbon Program (CCP) (2002-2011) and eco-region distributions of Canada’s landmass. For description of flux sites, see Table 1. 3. .......................................................... 19 Figure 2.1 Schematic view of the components of the climate system, their processes and interactions (Treut, 2007). ....................................................................................................... 26 Figure 2.2 Global CO2 budget from 1959 to 2006. Upper panel (A): CO2 emissions to the atmosphere (sources) as the sum of fossil fuel combustion, land-use change, and other emissions. Lower panel (B): The fate of the emitted CO2, including the increase in atmospheric CO2 plus the sinks of CO2 on land and in the ocean (Canadell et al., 2007). .... 27 Figure 2.3 Temporal and spatial scales of different approaches. ........................................... 47  xvii Figure 2.4 Structure of the EASS model. Three components (soil, vegetation and the atmosphere) are considered in EASS, which are integrated with two interfaces. The right panel illustrated energy fluxes between these three components. LE, H, Rs, Rl and G are the latent heat flux, sensible heat flux, shortwave radiation, longwave radiation, and soil conductive heat flux, respectively; the subscripts g and c present the energy fluxes at soil- canopy and canopy-atmosphere interfaces, respectively. The left panel describes soil water fluxes. The symbol F represents conductive water flux between soil layers, and F0 represents the incoming water flux from the surface to the top soil layer (i.e. the actual infiltration rate ), and Fb is the water exchange (drainage or capillary rise) between the bottom soil layer and the underground water (Chen et al., 2007a). ........................................................................... 50 Figure 2.5 An up-scaling framework synthetically integrating eco-hydrological and footprint modeling, remote sensing and land surface measurements. ................................................... 55 Figure 3.1 Seasonal and interannual variability of climatic and environmental variables measured at the research site for 1997-2006: (a) Monthly totals of precipitation (P), (b) air temperature (Ta) above the canopy, (c) soil temperature (Ts) at 5-cm depth, (d) soil water content of 0-60 cm soil layer (θ was weighted by the distribution of roots for the top 0-60 cm), (e) vapor pressure deficit (D) above the canopy, and (f) photosynthetically active radiation (Q) above the canopy.  The dotted and dashed lines in panel (e) are the estimated soil water content at field capacity (soil water matric potential = -0.01 MPa) and permanent wilting point (soil water matric potential = -1.5 MPa), respectively (Coops et al. 2007b). ... 70 Figure 3.2 Climatic features for the period 1965-2004, observed at Campbell River airport climate station (Latitude: 49º 57.000’N, Longitude: 125º16.200’W, Elevation: 105.5 m), which is 10 km NE of the research site (DF49), showing interannual and seasonal trends and anomalies. (a) Air temperatures were calculated on the basis of calendar year: the dry season extends from April - September while the remainder of the year is the wet season. (b) The winter, spring, summer and autumn seasons were classified following Mote (2003) as December-January-February, March-April-May, June-July-August and September-October- November, respectively. Hydrological-year precipitation was calculated from October of the preceding year to September of the nominal year. (c) The hydrological year was broken I  xviii down into the wet and the dry seasons. The dotted lines show the 40-year mean values for each season and the thick lines show measured air temperatures above the canopy and precipitation at DF49. ............................................................................................................. 71 Figure 3.3 Nine-year averaged monthly composite diurnal cycles of net ecosystem productivity (NEP, thin solid line) and gross primary productivity (GPP, thick solid line). Vertical bars show ±1 standard deviation of annual variation. ............................................... 73 Figure 3.4 Monthly totals of net ecosystem productivity (NEP), gross primary productivity (GPP) and ecosystem respiration (Re) from November 1997 to the end of 2006. .................. 74 Figure 3.5 Seasonality of C fluxes and the corresponding interannual variability for 1998- 2006. (a) Time series of the 9-year-average value of the ratio of the 4-month (current month and the 3 following months) total C flux (NEP, GPP and Re) to the annual total C flux for each month of the year. (b) Time series of the 9-year-average value of the ratio of the 4- month (current month and the 3 following months) average value of σ for monthly C fluxes to the annual value of sigma for each month of the year. (c) Standard deviation (σ) of monthly C fluxes. (d) Standard deviation of monthly averages of environmental variables. NEP is net ecosystem productivity, GPP is gross primary productivity, Re is ecosystem respiration, Ta is the air temperature above the canopy, θ is volumetric soil water content of 0-60 cm soil layer weighted by the distribution of roots and D is vapor pressure deficit above the canopy. The gray area (May - September) in panels (b) and (d) shows the period with high σ values of C fluxes. ....................................................................................................... 75 Figure 3.6 Monthly gross primary productivity to ecosystem respiration (GPP /Re) ratios for 1998 - 2006 and the 9-year mean. ........................................................................................... 76 Figure 3.7 Monthly departures of gross primary productivity (GPP), ecosystem respiration (Re), and net ecosystem productivity (NEP) fluxes from the 1998 - 2006 means. Linear fits are shown as dashed lines for each panel. ............................................................................... 77 Figure 3.8 Photosynthetically active radiation (Q) and air temperature dependences of gross primary productivity (GPP) and ecosystem respiration (Re) rates over annual cycles, at DF49, Nov. 1997 - Dec. 2006. (a) Relationship between monthly means of gross primary  xix productivity (GPP) and daytime (Q > 5 μmol m-2 s-1) Q, and the solid and dashed lines are the linear regression fits of EC-derived GPP (sum of gap-filled daytime net ecosystem productivity and calculated daytime Re and modeled GPP, respectively; (b) Relationship between residuals of monthly GPP and monthly mean soil water content (θ, at the 0-60-cm depth), and the solid and dashed lines are the linear regression trends for dry-soil months and wet-soil months, respectively; (c) Responses of monthly mean Re to Ta , and the solid and dashed lines are the exponential regression fits of gap-filled measured Re and modeled Re using, respectively; (d) Relationships of residuals of monthly Re with monthly mean θ, and the solid and dashed lines are the linear regression trends for dry-soil months and wet-soil months, respectively;  (e) Responses of monthly mean soil respiration (Rs) to soil temperature (Ts) at the 5-cm depth. The residuals of GPP or Re are the differences between the observed values and the corresponding values that are predicted by the fitted models (GPP - Q model or Re-Ta relationships, respectively) and thus they represent the variance that is not explained by the relationships (see text). The dry-soil months are July-October and the remainder of months is the wet-soil months................................................................................................. 81 Figure 3.9 Net ecosystem productivity (NEP) vs. air temperature (Ta) above the canopy (a) and daytime photosynthetically active radiation (Q) above the canopy (b) on monthly time scales. The data shown are for Nov 1997 - Dec 2006. The legend is the same as in Figure 3.8. The results of linear regression analyses for the same data are shown in Table 3.5. ............. 83 Figure 3.10 Relationships of monthly mean daytime net ecosystem productivity (NEPd) vs. monthly mean soil water content (θ30-50cm) at the 30-50-cm depth (a) and vapor pressure deficit (D) above the canopy (b). Linear regression fits applied to the dry season months (April - September). The data shown are for DF49, Nov. 1997 - Dec. 2006. ........................ 84 Figure 3.11 Relationships of the monthly mean gross primary productivity (GPP), daytime ecosystem respiration (Red), and daytime net ecosystem productivity (NEPd) with soil water content (θ30-50cm) at the 30-50-cm depth under severely dry conditions (θ30-50cm < 0.1 m 3  m -3 ). Solid, long-dashed, and short-dashed lines are linear fits for GPP, Red, and NEPd vs. θ30-50cm, respectively. ............................................................................................................................ 86  xx Figure 3.12 Departures in annual gross primary productivity (GPP) and ecosystem respiration (Re) fluxes from their respective 1998 - 2006 means. The linear fit is shown as a dashed line while the 1:1 line is shown as a dotted line. ........................................................ 88 Figure 3.13 Responses of annual C fluxes (GPP: gross primary productivity, Re: ecosystem respiration, and NEP: net ecosystem productivity) to mean spring temperature (a, 1998 - 2006) and annual mean diffuse photosynthetically active radiation (Qdif) above the canopy (b, 2000 - 2006). ........................................................................................................................... 90 Figure 4.1 Comparison of monthly mean values of climatic and environmental variables measured at the 58-year-old West Coast Douglas fir stand (DF49) for the fertilized year 2007 and for the 1998-2006 means. (a) Precipitation (P), (b) air temperature above the canopy (Ta), (c) 0-60-cm soil water content (θ), (d) water vapor saturation deficit above the canopy (D) and (e) daytime total downward photosynthetically active radiation (Q). The error bars are ± 1 standard deviation for all the 9-year data (1998-2006). ............................................ 104 Figure 4.2 Comparison of monthly mean values of component C fluxes measured at the 58- year-old West Coast Douglas-fir stand (DF49) for the fertilized year 2007, for the previous year 2006 and for the 1998-2006 means. (a) Gross primary productivity (GPP), (b) ecosystem respiration (Re), (c) belowground ecosystem (soil) respiration (Rs), and (d) net ecosystem productivity (NEP).  The error bars are ± 1 standard deviation for all the 9-year data (1998-2006). .................................................................................................................. 105 Figure 4.3 BEPS predicted and EC tower measured daily gross primary productivity (GPP), ecosystem respiration (Re) and net ecosystem productivity (NEP) at DF49 for1998 to 2007, (a) for the EnKF applied years1998-2005, (b) for the model validation year 2006 and (c) for the first-post-fertilization year 2007. .................................................................................... 108 Figure 4.4 Comparison of BEPS predicted and EC tower measured daily gross primary productivity (GPP), ecosystem respiration (Re) and net ecosystem productivity (NEP) at DF49 for non-fertilized years, (a) - (c) for the EnKF applied years 1998-2005 and (d) - (f) for the validation year 2006. ....................................................................................................... 109  xxi Figure 4.5 BEPS predicted and measured daily soil respiration (Rs) using an auto-chamber at DF49, (a) for the ENKF applied years 2003-2005, (b) for the model validation year 2006 and (c) for the first-post-fertilization year 2007. ......................................................................... 110 Figure 4.6 Comparisons of measured and modeled annual gross primary productivity (GPP), ecosystem respiration (Re), soil respiration (Rs), and net ecosystem production (NEP). The error bars show the standard errors (s.e.). The s.e. values for the EC measured annual C fluxes were estimated using the method of Schwalm et al. (2007), while the s.e. values for modeled annual values were estimated using the ±1 standard deviation of optimized model parameters. (a) for the EnKF applied years, (b) for the model validation year 2006 and (c) for the first-post-fertilization year 2007. .................................................................................... 111 Figure 4.7 BEPS predicted and EC tower measured net ecosystem productivity (NEP) in the early (a), middle (b) and late (c) stages of the growing season in 2006 at DF49. ................ 113 Figure 4.8 Comparison of BEPS predicted (no fertilization) and EC tower measured daily gross primary productivity (GPP), ecosystem respiration (Re) and net ecosystem productivity (NEP) at DF49 for the fertilized year 2007. ......................................................................... 115 Figure 4.9 Seven-day moving averages of BEPS predicted and measured C component fluxes for the fertilization year 2007 at DF49. The C component fluxes are (a) ecosystem respiration (Re), (b) soil respiration (Rs), (c) gross primary productivity (GPP), and (d) net ecosystem respiration (NEP). The N fertilization effects on the C component fluxes were estimated as the differences between the measured C fluxes and their corresponding modeled values. ................................................................................................................................... 116 Figure 4.10 Cumulative N fertilization effects on C component fluxes with ±1 standard errors in 2007 at DF49, which were estimated as the differences between the measured C fluxes and the corresponding modeled values. Standard errors were estimated using the ±1 standard deviation of optimized model parameters. The C component fluxes are gross primary productivity (GPP), ecosystem respiration (Re), soil respiration (Rs) and net ecosystem respiration (NEP)................................................................................................. 117  xxii Figure 5.1 The location of the three flux tower sites of the coastal BC Flux Station on Vancouver Island, part of the Canadian Carbon Program. ................................................... 132 Figure 5.2 Land surface heterogeneites of vegetation (tree canopy density maps or orthophoto) and topography for the three study tower sites.  Continuous canopy tree densities for 25  25 m pixels were determined from individual canopy tree crowns identified on 0.5- m-resolution orthophotos (as of 1999) for the DF49 and HDF88 sites, while for HDF00 only the orthophoto is shown.  The thin isolines show elevation contours while the thick isolines show 80% (outer) and 50% (inner) contours of the multi-year cumulative NEP-weighted footprint-climatology contours. Star symbols show the tower locations. And the locations of permanent sample plots (numbered square) are also shown. ................................................ 134 Figure 5.3 Annual mean daytime and nighttime wind roses (with 10º intervals) for the three sites; for 1998-2006, 2001-2006 and 2002-2006 at DF49, HDF00 and HDF88, respectively. The dashed circles with numbers (5, 10, 15, …) represent the percentages of probability of wind direction for each 10º interval. Left hand panels (A-1, B-1 and C-1) and right hand panels (A-2, B-2 and C-2) show data without and with applying fraction threshold ( ) to nighttime, respectively. The values of for DF49, HDF88 and HDF00 are 0.3, 0.16 and 0.08 m s -1 , respectively. 71%, 70% and 66% data were removed by excluding the data when < for DF49, HDF88 and HDF00, respectively. .......................................................... 144 Figure 5.4 Five-day ensemble means of CO2 component fluxes for the three sites for 1998 – 2006, 2002 – 2006, and 2001 – 2006 for DF49, HDF88, and HDF00, respectively. Measured net ecosystem productivity (NEP) was partitioned into gross primary productivity (GPP) and ecosystem respiration (Re) using procedures described in Morgenstern et al. (2004). ......... 145 Figure 5.5 Crosswind-integrated footprints for a typical day (unstable conditions: panel A) and night (stable conditions: panel B). Model inputs are monthly median half-hourly meteorological values (July, 2004) as listed in Table 5.2. .................................................... 147 Figure 5.6 Cumulative footprint contours for typical half-hourly unstable (daytime, solid line) and stable (nighttime, dotted line) conditions for the three sites, showing the footprint areas for given percentages (20%, 40%, 60%, 80%, 90%, 95%, and 99%) of the total flux  thu* thu* *u thu*  xxiii measured by the EC flux tower. The triangle and square are the location of footprint peak during daytime and nighttime, respectively. ......................................................................... 148 Figure 5.7 Nine-year (1998-2006) mean monthly footprints ( ) weighted by NEP (left panel), GPP (middle panel) and Re (right panel) for January, April, July and October for DF49. The corresponding cumulative footprint contours overlie on their corresponding footprint maps. ...................................................................................................................... 149 Figure 5.8 Relationship between nine-year (1998-2006) mean monthly cumulative NEP weighted and pure (unweighted) footprints and their corresponding areas for DF49. ......... 151 Figure 5.9 Multi-year mean of annual C-fluxes-weighted and pure (unweighted) footprint climatologies for the three sites: DF49, 1998-2006 (left hand panel); HDF88, 2002-2006 (middle panel); and HDF00, 2001-2006 (right hand panel). The corresponding cumulative footprint contours overlie on their corresponding footprint maps. ....................................... 154 Figure 5.10 Relationship between the multiple-year mean cumulative C-fluxes-weighted and pure (unweighted) footprints and their corresponding areas for the three sites. ................... 156 Figure 5.11 Multiple-year standard variation ( ) and coefficient of variation ( ) of the NEP-weighted footprint climatology for the three sites, showing the overall patterns and interannual variations in the spatial structure of the weighted annual footprint climatology during the measurement periods. The measurement periods for DF49, HDF88 and HDF00 are 1998-2006, 2002-2006, and 2001-2006, respectively. was calculated as the ratio of the standard deviation ( ) to the multiple-year mean ( ) (see text for detail). ............................................................................................... 157 Figure 5.12 Remotely sensed monthly mean GPP maps at 10 m resolution surrounding the DF49 tower for January 2007 (A), April 2006 (B), July 2006 (C) and October 2006 (D), overlaid by corresponding cumulative monthly daytime pure (unweighted) footprint climatology contours. ............................................................................................................ 159  ),( yx ),( yxCV ),( yxCV ),( yx ),( yx  xxiv Figure 5.13 Comparisons in ecosystem-scale GPP estimates between integrated from remotely sensed GPP maps and directly derived from EC measurements.  Both footprint integration and “Equal” integration were calculated over the same area of ΩΠ. The source area ΩΠ is the area contains all the grids (cells) whose footprint function values are larger than the values at cutoff point, where a given cumulative fraction Π of 99% is achieved (see text for detail). (A) 1:1 plot of “equally”-integrated versus EC-derived half-hourly GPP; (B) 1:1 plot of footprint-integrated versus EC-derived half-hourly GPP; (C) temporal variations in estimates of monthly sums of GPP. The inset shows the cumulative curves of the three estimations of GPP from April 1 st  2006 to March 31 st  2007. The monthly mean root sensor location bias δ (circle symbols with thick solid line) is also shown in (C)........................... 162 Figure 6.1 False color composite images (bands 2, 3 and 4) of Landsat at 30-m resolution for each area (6 km  6 km) centred at individual towers. The scene path and row and the acquired dates are shown in Table 6.3. The 90% contour of cumulative annual footprint climatology and a circle with 1-km radius centred at the tower overlie their corresponding Landsat imagery maps. For description of abbreviation of vegetation fuctional type and site ID, see Table 6.1. .................................................................................................................. 174 Figure 6.2 Annual footprint climatology PDF maps in 2006 for the 12 CCP sites. The cumulative footprint contours overlie their corresponding footprint PDF maps. For description of abbreviation of vegetation fuctional type and site ID, see Table 6.1............. 182 Figure 6.3 Exponential relationship between the annual cumulative footprint climatology percentages and their corresponding areas for the 12 sites. For description of abbreviation of vegetation fuctional type and site ID, see Table 6.1. ............................................................ 184 Figure 6.4 Seasonal variations in the size of 90% cumulative footprint climatology for the 12 sites, presented by the maximum and minimum of monthly footprint climatology sizes. For description of abbreviation of vegetation fuctional type and site ID, see Table 6.1............. 185 Figure 6.5 Variations of annual sensor location biases using NDVI (a, c) and EVI (b, d) as flux surrogates with corresponding varying cumulative footprint climatology percentages (a, b) and footprint climatology areas (c, d) for the 12 sites. The reference area of interest in Eq.   xxv (6.10) is assumed to be 90% of annual footprint climatology (A90) for individual sites, where A90 is the footprint area for 90% cumulative footprint percentage. For description of abbreviation of vegetation fuctional type and site ID, see Table 6.1. ................................... 187 Figure 6.6 Variations of annual sensor location biases using NDVI (a, c) and EVI (b, d) as flux surrogates with corresponding varying cumulative footprint climatology percentages (a, b) and footprint climatology areas (c, d) for the 12 sites. The reference area of interest in Eq. (6.10) is assumed to be a circular area with 1-km radius centred at the individual towers. For description of abbreviation of vegetation fuctional type and site ID, see Table 6.1............. 189 Figure 6.7 Variations of annual sensor location biases using NDVI (a) and EVI (b) as flux surrogates for the 12 FCRN sites with corresponding varying reference circular areas of interest in Eq. (6.10) centred at the individual towers from radii of 0.1 km to 3 km, while the cumulative footprint percentage in Eq. (9) is fixed to be 90%. For description of abbreviation of vegetation fuctional type and site ID, see Table 6.1. ........................................................ 190 Figure 6.8 Seasonal variations in sensor location biases (SLB) using NDVI and EVI as flux surrogates with 90% of annual footprint climatology for the 12 sites, presented by the maximum and minimum of monthly SLB. The reference areas of interest in Eq. (10) are 90% of annual footprint climatology area and a circular area with 1-km radius centred at the individual towers, respectively. For description of abbreviation of vegetation fuctional type and site ID, see Table 6.1. ..................................................................................................... 191 Figure 7.1 Coefficients of a typical semivariogram model. The x axis shows the scale of variation, while the y axis shows the amount of variation. ................................................... 204 Figure 7.2 NDVI (Normalized Difference Vegetation Index) maps at a 30-m resolution for each area (6  6 km) centred at individual CCP towers. The scene path and row and the acquisition dates are shown in Table 6.3. The contours from inner to outer are the 50%, 80%, 90% and 99% annual cumulative footprint climatology, respectively. For description of abbreviation of vegetation type and site name, see Table 6.1. ............................................. 207   xxvi Figure 7.3 Footprint areas of 50%, 80%, 90% and 99% annual footprint climatology for the twelve CCP sites in 2006. For description of abbreviation of vegetation type and site name, see Table 6.1. ........................................................................................................................ 208 Figure 7.4  GDEM maps at a 30-m resolution for each area (6  6 km) centred at individual CCP sites. The contours from inner to outer are the 50%, 80%, 90% and 99% annual cumulative footprint climatology, respectively. For description of abbreviation of vegetation type and site name, see Table 6.1. ........................................................................................ 209 Figure 7.5 Frequency distribution of DEM for the twelve CCP sites within the areas of 50%, 80%, 90% and 99% annual footprint climatologies and the circular areas centred at individual towers with radius of 1, 2, and 3 km. For description of abbreviation of vegetation type and site name, see Table 6.1. ....................................................................................................... 210 Figure 7.6 Land cover maps of the ten data-available CCP sites at a 30-m resolution for each area (6  6 km) centred at individual CCP sites. The contours from inner to outer are the 50%, 80%, 90% and 99% annual cumulative footprint climatology, respectively. For description of abbreviation of vegetation type and site name, see Table 6.1. ...................... 212 Figure 7.7 Frequency distribution of land cover types for the ten data-available CCP sites in the areas of 50%, 80%, 90% and 99% annual footprint climatologies. For description of abbreviation of vegetation type and site name, see Table 6.1. ............................................. 213 Figure 7.8 Directional semivariograms of NDVI over the 90% annual footprint climatology area. For description of abbreviation of vegetation type and site name, see Table 6.1. ....... 217 Figure 7.9 Omnidirectional semivariogram parameters of NDVI over the areas of 50%, 80%, 90% and 99% annual footprint climatologies and the circular areas centred at individual towers with radius of 1, 2, and 3 km. For description of abbreviation of vegetation type and site name, see Table 6.1. ....................................................................................................... 218 Figure 7.10 The window-averaged NDVI of varying size with the tower location at the center of the window. Vertical lines show the window width (d) matching the 90% annual    xxvii footprint climatology area (A90), where      . For description of abbreviation of vegetation type and site name, see Table 6. 1. ....................................................................................... 220 Figure 8.1 Daily footprints for two arbitrarily selected days (a sunny and cloudy days in summer) at four heights (95, 55, 33 and 22 m) on the ETL tower. Simulated concentration footprints (a, b) along the centerline of mean flow (upwind of the tower) and (c, d) across the wind direction from the the centreline of the mean flow. The parameters for characterising the two day’s atmospheric stability are listed in Table 8.2. .................................................. 242 Figure 8.2 Monthly footprint areas for 90% cumulative footprint climatologies of the four measurement heights (95, 55, 33 and 22 m) on the ETL tower and for the CL tower at 28-m height for the period 2006-2008. Asterisks, circles and triangles represent the years: 2006, 2007 and 2008, respectively. The lines connect the points for 2007. ................................... 244 Figure 8.3 Monthly footprint climatologies ( ) for January, April, July and October for measurements at the 95-m height (left panels) and at 22-m height on the ETL tower (middle panels) and for measurements at 28-m height on the CL tower (right panels). The corresponding cumulative footprint contours (from inner to outer: 50%, 65%, 80% and 95%) are also shown. ...................................................................................................................... 245 Figure 8.4 Annual footprint climatologies ( ) for the year 2008 for measurements on the ETL tower ((a) 95 m, (b) 56 m, (c) 33 m and (d) 22 m) and (e) for measurements at 28-m height on the CL tower. The corresponding cumulative footprint contours (from inner to outer: 50%, 65%, 80% and 95%) are also shown. ................................................................ 247 Figure 8.5 Relationships between annual footprint percentage and area for measurements at four eights (95, 56, 33 and 22 m) on the ETL tower for the years 2006-2008. .................... 248 Figure 8.6 LC map for a 1000 × 1000 km region surrounding the ETL 106-m tall tower, in Saskatchewan (54°21'N, 104°59'W) at 1 km resolution. The circles are the ETL and CL high- precision CO2 concentration towers, while the triangles are the cluster of EC towers. LC data are from the GLC2000, which are 1-km resolution products based on the Food and Agriculture Organization Land Cover Classification System (Bartholome  ́and Belward,  xxviii 2005). The LC classification schemes are shown in Table 8.1.  The original 26 LC type classification (panels (c) and (d)) are regrouped into 8 types (panel (b)) and further into 4 types (panel (a)). 95% cumulative footprint climatology contours for the years 2006-2008 for the ETC tower at heights of 22 m (panel (a)), 33 m (panel (b), 56 m (panel (c)) and 95 m (panel (d)). The contours for the CL tower at the 28-m height are shown in panel (a). ....... 250 Figure 8.7 The percentages of land cover types calculated with footprint weighting factors (using Eqs. (8.7-8.8))  over 95% footprint climatology areas for measurements at four levels on the ETL tower (95, 56, 33 and 22 m above the ground) and at 28-m height on the CL tower, (a) eight types LC classification and (b) four types LC classification. Each error bar represents the standard deviation of calculated percentages during the years of 2006-2008 for the ETL tower and the years of 2003-2009 for the CL tower. .............................................. 251 Figure 8.8 The differences in percentages of land cover types over 95% footprint climatology areas between footprint weighting (using Eqs. (8.7-8.8)) and no weighting (using Eqs. (8.7-8.8) but setting ,C i to 1) for measurements at four heights on the ETL tower (95, 56, 33 and 22 m) and at the 28-m height on the CL tower, (a) eight types LC classification and (b) four types LC classification. Each error bar represents the standard deviation of calculated percentages during the years of 2006-2008 for the ETL tower and the years of 2003-2009 for the CL tower. ........................................................................................................................ 252   xxix List of Abbreviations, Acronyms, and Symbols* Symbol Units Description  AB n/a Alberta Ac µmol C m -2  s -1  Rubiso-limited photosynthesis rate Aj µmol C m -2  s -1  light-limited gross photosynthesis rate Amax µmol C m -2  s -1  light-saturated asymptote for canopy photosynthetic capacity AMMA n/a African Monsoon Multidisciplinary Analysis AMSPEC n/a automated multi-angular spectral radiometer platform Anet µmol C m -2  s -1  net photosynthetic rate ASTER n/a Advanced Spaceborne Thermal Emission and Reflection Radiometer AVHRR n/a Advanced Very High Resolution Radiometer BATS n/a Biosphere-atmosphere Transfer  Scheme (model name) BC n/a British Columbia BOREAS n/a Boreal Ecosystems-Atmosphere Study bsc dimensionless ratio of the molecular diffusivity of water to that of CO2 (bsc =1.6) c μmol CO2 mol -1  CO2 molar mixing ratio C n/a carbon C/N  dimensionless ratio of carbon to nitrogen Ca mol CO2 mol -1  CO2 mole fraction of ambient air in the canopy CBL n/a convective boundary layer Cc mol CO2 mol −1  intercellular CO2  mole fraction CC n/a crown closure CCP n/a Canadian Carbon Program Ci mol CO2 mol -1  CO2 mole fraction of ambient air in the intercellular spaces CL n/a Candle Lake  (CO2 concentration tower’s name) CLASS n/a Canadian LAnd Surface Scheme (model name) CO2 n/a carbon dioxide COST n/a cosine approximation model CR n/a Campbell River Cs µmol CO2 mol -1  CO2 mole fraction at the leaf surface Csoil, n/a soil C pool CV dimensionless coefficient of variation CWH n/a coastal western hemlock D kPa vapor pressure deficit DAS n/a data-acquisition system DEM n/a digital elevation model DF49 n/a Douglas fir (established in 1949) (flux tower’s name) Do kPa an empirical parameter determining the sensitivity of stomatal conductance to water vapor saturation deficit, DOC mol C dissolved organic carbon Ds kPa water vapor saturation deficit at the leaf surface, ds/dt mm s -1  change of water storage EASS n/a ecosystem atmosphere simulation scheme (model name)  xxx EC n/a eddy-covariance EOBS n/a Eastern Old Black Spruce (flux tower’s name) EOS n/a Earth Observing System EOSD n/a Earth Observation for Sustainable Development of forest EPL n/a Eastern peat land (flux tower’s name) ESRL n/a the Earth System Research Laboratory ET mm evapotranspiration ETL n/a East Trout Lake (CO2 concentration tower’s name) ETM+ n/a Enhanced TM Plus FCRN n/a Fluxnet Canada Research Network F MA -1 t -1  vertical flux cf m -2  concentration footprint function Ff  m -2  flux footprint function y Ff m -1  crosswind-integrated footprint rsGPPF ,  MA -1 t -1  Remotely-sensed spatially-distributed gross primary production fPAR dimensionless fraction of photosynthetically active radiation fw dimensionless bulk water availability factor G Wm -2  soil conductive heat flux gc mol m -2  s -1  conductance to intercellular spaces for CO2 diffusion from the ambient air GCM n/a General Circulation Model GDEM n/a Global DEM go mol m -2  s -1  stomatal conductance GOES n/a Geostationary Operational Environmental Satellites GOSAT n/a Greenhouse Gases Observing Satellite GPP MA -1 t -1  Gross primary production GRS n/a Grassland (flux tower’s name) gs mol m -2  s -1  bulk stomatal conductance H Wm -2  sensible heat flux hc m canopy height He m elevation above sea level hm m EC sensor height IBIS n/a Integrated BIosphere Simulator (model name) J  µmol C m -2  s -1  electron transport rate 25 maxcJ  µmol C m -2  s -1  maximum electron transport rate at 25 °C 0k  dimensionless effective eddy diffusivity coefficient Kc Ms -1  the eddy diffusivity Kc mol mol −1  Michaelis-Menten constants for CO2 Ko mol mol −1  Michaelis-Menten constants for O2 L m Obukhov length LAI m 2 m -2  leaf area index LC n/a land cover LE Wm -2  the latent heat flux LiDAR n/a Light Detection and Ranging LSM n/a land surface model LUE µmol CO2 µmol -1  quanta light use efficiency  xxxi m dimensionless exponent of the wind velocity power law in the SAFE model (Eq. 5.7) m dimensionless a dimensionless slope of the Ball-Berry model for predicting stomatal conductance, which is related to the intercellular CO2 mole fraction by Ci/Cs = 1 -1/m at maximal stomatal opening (when both Ds and go are zero and fw =1) in the EASS model (Eq. A3) MB n/a Manitoba METI n/a Ministry of Economy, Trade, and Industry MODIS n/a Moderate Resolution Imaging Spectroradiometer MSS n/a Multispectral Scanner System N n/a Nitrogen n dimensionless exponent of the eddy diffusivity power law in Eq. 5.7 NACP n/a North American Carbon Program NASA n/a United States National Aeronautics and Space Administration NB n/a New Brunswick NCAR n/a National Center for Atmospheric Research NCEP n/a National Centers for Environmental Prediction NDVI n/a normalized difference vegetation index (vegetation index), NEE MA -1 t -1  net ecosystem exchange of CO2 NEP MA -1 t -1  net ecosystem production (GPP – Re) NOBS n/a Northern Old Black Spruce (flux tower’s name) NPP MA -2 t -1  Net primary production (GPP – Ra) Nr, M all nitrogen species, such as NH3 and NOx Nr0 g N m −2  leaf area the leaf Rubisco-N  in the top canopy OA n/a Old Aspen (flux tower’s name) OBF n/a Old Balsam Fir (flux tower’s name) Oc  mol mol -1  intercellular O2 mole fractions OJP n/a Old Jack Pine (flux tower’s name) OMW n/a Old mixed wood (flux tower’s name) ON n/a Ontario, P mm precipitation p n/a statistically significant level PAR µmol quanta m -2  s -1  photosynthetically active radiation PBL n/a planetary boundary layer PILPS n/a Inter-comparison of Land-surface Parameterization Schemes PRI  dimensionless photochemical reflectance index PVM n/a potential vegetation model (model name) pw dimensionless a time-varying parameter that quantifiers the fractional departure of GPP from the mean  GPP vs. Q relationship Q µmol quanta m -2  s -1  Photosynthetically active radiation Q10 dimensionless factor by which Re increases for a 10  o C increase in temperature QC n/a Quebec R mm s -1  percolation and surface runoff RADAR n/a RAdio Detection And Ranging Rb dimensionless the bulk Richardson number Rd µmol m -2  s -1  daytime leaf dark respiration Re MA -2 t -1  ecosystem respiration  xxxii Red MA -2 t -1  daytime Re Rg µmol C m -2  s -1  growth respiration Rh µmol C m -2  s -1  heterotrophic respiration Rl Wm -2  longwave radiation Rm µmol C m -2  s -1  maintenance respiration RMSE n/a root mean square error Rn Wm -2  net radiation Rref µmol C m -2  s -1  respiration rate at a reference temperature refT (here 10 o C was used) Rref µmol C m -2  s -1  respiration rate at a reference temperature (here 10 o C was used) Rs Wm -2  shortwave radiation rw dimensionless time-varying parameter that quantifiers the fractional departure of Re from the mean ‘wet’ Re vs. Ts relationship SAFE-C n/a Simple Analytical Footprint model based on Eulerian coordinates for the scalar concentration (model name) SAFE-F n/a Simple Analytical Footprint model based on Eulerian coordinates for the scalar flux (model name) SAR n/a Synthetic-aperture RADAR SAT n/a sonic anemometer-thermometer Sc µmol C mol -1  CO2 storage SiB n/a Simple Biosphere model (model name) SiB2 n/a Simple Biosphere model version 2.0 (model name) SK n/a Saskatchewan, SLB n/a sensor location bias SOBS n/a Southern Old Black Spruce (flux tower’s name) SOC n/a soil organic carbon SSA n/a Southern Study Area SVATS n/a soil-vegetation-atmosphere transfer schemes Ta °C land surface air temperature Ta °C air temperature TBM n/a terrestrial biogeochemistry model TDL n/a tunable diode laser Te Year A.C. stand established time TIR n/a remotely-sensed thermal-infrared TIROS n/a Television Infrared Observation Satellite TM n/a the Thematic Mapper Tref °C a reference temperature (here 10 o C was used) for soil respiration calculation Ts, °C soil temperature u m s -1  wind speed 0u  dimensionless effective speed of plume advection u  m s -1  friction velocity thu*  m s -1  friction wind speed (u*) threshold for EC flux calculation )(xup  m s -1  effective velocity of the plume. Vcmax µmol m -2  s -1  maximum carboxylation rate 25 maxcV  µmol m -2  s -1  maximum carboxylation rate at 25 °C VFT  n/a vegetation function type  xxxiii WD ° wind direction WMO n/a World Meteorological Observation WPL n/a Western peat land (flux tower’s name) WPP39 n/a White Pine Plantation (established 1939) (flux tower’s name) fxmax  m the distance of the peak value of the footprint function from the tower z0 m surface roughness length Zh m the maximum height of the PBL zm  m height of the sensor θ m3 m-3 soil water content Ψ MPa soil water matric potential *  mol CO2 mol −1  CO2 compensation point without dark respiration _ a  μmol CO2 mol -1  mean molar density of dry air,  ____ ''cw  m s-1mol mol-1 covariance of vertical velocity ( w , in m s-1) and CO2 molar mixing ratio ( c , in μmol mol-1) φmax  % max is the peak value of the footprint function w  m s-1 vertical velocity   n/a quantum yield v  m s -1  lateral wind speed y  m s -1  standard deviation of the plume in the y dimension ),( zxc  μmol CO2 mol -1  crosswind integrated concentration,   n/a Gamma function * n/a represents not applicable; MA -2 t -1 : M - mass unit, A - area unit, and t - time unit   xxxiv Acknowledgements I am gratefully indebted to my mentors, advisors and friends, Dr. Nicholas C. Coops and Dr. T. Andy Black, who have guided me over the years at the University of British Columbia during which I was pursuing my Ph.D. and before as a research associate. I also extend special appreciation to the other two members of my supervisory committee, Dr. Rob Guy and Dr. Mark Johnson for sharing their wisdom and enthusiasm. Equally important has been the patience, support and understanding of my family and especially my wife, confidant and best friend, Jie Teng, who provided support and encouraged me to never give up or lose focus. This dissertation is dedicated to my two lovely sons, Derek Yuteng Chen and Alex Jia Teng Chen, who have to leave their friends in Toronto and move to Vancouver accompanying me and faced the difficulty of starting a new school and daycare and had to make new friends in Vancouver. They have grown up during my Ph.D research. This research was supported by Alexander Graham Bell Canada Graduate Scholarship (CGS, $100,500 for three years) by the Natural Sciences and Engineering Research Council of Canada (NSERC), and Four-Year Fellowships (FYF) for PhD Students, Asa Johal Graduate Fellowship, George S. Allen Memorial, and Donald S. Mcphee Fellowship by UBC. I watched Nicholas C. Coops’ Lab and T. Andy Black’s Lab, evolve over time and I would like to recognize and thank all past and current members for their support and camaraderie. I am also appreciative of the entire Department of Forest Resources Management and the Faculty of Forestry, UBC, for providing an outstanding academic environment, and the administrative staff, for supporting efforts to make this the best graduate program in the country. Numerous other colleagues and co-authors were instrumental in my research and I am thankful to everyone for discussion, manuscript comments, data, support and diversions. Those that have provided contributions to me for research related to this dissertation include:  xxxv Hank A. Margolis (Faculté de Foresterie et de Géomatique, Université Laval, Québec); Jing M. Chen (Department of Geography and Program in Planning, University of Toronto, Ontario); Alan G. Barr (Climate Research Branch, Meteorological Service of Canada, Saskatoon, Saskatchewan);  Brian D. Amiro (Department of Soil Science, University of Manitoba, Winnipeg, Manitoba); M. Altaf Arain (School of Geography and Earth Sciences and McMaster Center For Climate Change, McMaster University, Hamilton, Ontario); Charles P.-A. Bourque (Faculty of Forestry and Environmental Management, University of New Brunswick, Fredericton, New Brunswick); Lawrence B. Flanagan (Department of Biological Sciences, University of Lethbridge, Lethbridge, Alberta); Peter M. Lafleur (Department of Geography, Trent University, Peterborough, Ontario); J. Harry McCaughey (Department of Geography, Queen’s University, Kingston, Ontario); Steven C. Wofsy (Department of Earth and Planetary Science, Harvard University, Cambridge); J. A. (Tony) Trofymow (Canadian Forest Service (Pacific Forestry Centre), Natural Resources Canada, Victoria); Douglas E.J. Worthy (Air Quality Research Branch, Meteorological Service of Canada). I also wish to acknowledge the help of numerous field crew, technicians, engineers and students involved in installation, maintenance, trouble shooting and data collection at all the study sites of the Canadian Carbon Program (CCP) network (previous Fluxnet Canada Research Network) and land owners for allowing access to field locations.  xxxvi Dedication  I dedicate this thesis to my family for their endless love and support and to the great teachers and professors I had.  1  Chapter  1: Introduction 1.1 Scientific background The Earth climate is a complex, interactive system, determined by a number of connected physical, chemical and biological processes occurring in the atmosphere, land and ocean. The terrestrial biosphere plays many pivotal roles in the coupled Earth system providing both positive and negative feedbacks to climate change (Treut et al., 2007). Terrestrial vegetation converts carbon dioxide (CO2) into organic carbon (C) by absorbing solar radiation through photosynthesis that would otherwise reside in the atmosphere as CO2, thereby moderating climate. Vegetation also transfers water between belowground reservoirs and the atmosphere to maintain precipitation and surface water flows through transpiration. The terrestrial C cycle is also closely linked to hydrological and nutrient controls (e.g. the most limiting nutrient: nitrogen (N)) on vegetation (Betts et al., 2000; Cox et al., 2000). 1.1.1 Terrestrial CO2 cycle The exchange of C between the biosphere, hydrosphere, geosphere and atmosphere is recognized as the global C cycle, which can be viewed as a series of reservoirs of C in the Earth system. Atmospheric CO2 has increased by a factor 1.4 since the beginning of the Industrial Revolution around 1750 and its growth rate is increasing rapidly (Conway and Tans, 2011; Honisch et al., 2009; Luthi et al., 2008; Petit et al., 1999). Two processes contributing to this rapid increase are fossil fuel emission and cement manufacturing (Andres et al., 2011) and land use change (deforestation) (Friedlingstein and Prentice, 2010;  Friedlingstein  et al., 2010; LeQuere et al., 2009). Average fossil fuel and cement manufacturing emissions were 7.7 ± 0.5 PgC yr –1  in the decade 2000–2009 with an average growth rate of 2.9% yr –1 (Andres et al., 2011; Friedlingstein  et al., 2010). This rate of increase of fossil fuel emissions is higher than during the 1990’s (1.0% yr–1). Emissions from land use change over the same decade are dominated by tropical deforestation, and are estimated at 0.9 ± 0.5 PgC yr –1 (Houghton, 2010). The human- 2  caused release of CO2 to the atmosphere is absorbed partly by the ocean, and partly sequestered by the terrestrial ecosystems, where it is stored in the vegetation and soils. Atmospheric CO2 concentrations grew by 4.0 ± 0.2 PgC yr –1  in the 2000s. The estimated mean ocean and land CO2 sinks are 2.3 ± 0.5 PgC yr –1  and 2.3 ± 0.9 PgC yr –1 , respectively for 2000–2009. Terrestrial ecosystems mediate a large part of CO2 exchange between the atmosphere and the land, by uptake due to photosynthesis and by release back to the atmosphere due to respiration (Prentice et al., 2001). Imbalances between gross primary productivity (GPP) and ecosystem respiration (Re) lead to positive or negative net ecosystem exchange of CO2 (NEE), i.e.                                                                                                                               (   ) with land surfaces being either CO2 sinks (NEE < 0) or sources (NEE > 0). Terrestrial ecosystems have accumulated 124 ± 59 Pg C since 1750 (Beer et al., 2010; Sarmiento et al., 2010, 2011). The increasing gain of C by the land is estimated to take place mainly due to enhanced photosynthesis at higher CO2 levels and N deposition, longer growing seasons in high latitudes, and the expansion and recovery of forests from past land use (Kleinen et al., 2010; Lemmen, 2009; Olofsson and Hickler, 2008; Jung et al., 2011). However, significant spatial and temporal variability in terrestrial CO2 fluxes limits our ability to reasonably diagnose and predict C cycling in the Earth system. The magnitudes of sinks and sources of terrestrial ecosystems fluctuate on diurnal, seasonal and longer time scales due to variable climate, land use change, disturbance, and changes in the age distribution and species composition of ecosystems (Battle et al., 2000; Arain et al., 2002; Law et al., 2002; Morgenstern et al., 2004; Humphreys et al., 2005, 2006; Urbanski et al., 2007). Terrestrial ecosystems modify the atmospheric C balance through many mechanisms and a detailed understanding of the interactive relationships involved in atmosphere–biosphere exchange is highly relevant to stand- 3  scale analysis and is needed to improve our knowledge of the global C cycle (Falk et al., 2008). In recent decades, scientists have become increasingly aware that terrestrial ecosystems’ vegetation, soil (Melillo et al., 1989; Knapp et al., 1993) and animals (Naeem et al., 1995; Hattenschwiler and Bretscher, 2001) play key roles in mediating the terrestrial C cycle. Plants being the primary producers provide mass and energy which gets transformed to other living organisms (Engel and Odum, 1999) within an ecosystem. The process of photosynthesis fixes atmospheric CO2 into the terrestrial ecosystem. Atmospheric CO2 enters the plant through stomatal openings that are controlled by a variety of environmental factors (Jarvis, 1976; Griffis et al., 2003). These factors include ambient temperature, atmospheric CO2 concentration, nutrient availability, soil water availability and forest age (Schimel, 1995; Prentice et al., 2001). Changes in atmospheric CO2 concentration and the corresponding changes in the climate have altered the magnitudes of terrestrial C balance component. For example, climate change induced increase in vegetation growth due to earlier spring and lengthened growing season was detected from the phase shift in the seasonal atmospheric CO2 cycle by Keeling et al. (1996) and satellite-based vegetation index analysis by Myneni et al. (1997). Studies indicate that an increase in atmospheric CO2 enhances photosynthesis (e.g. Woodward and Friend, 1988) and hence increases assimilation of atmospheric CO2 by the terrestrial vegetation. N availability to plants is another factor that can affect photosynthesis. This is because N is a limiting nutrient for plant growth in most ecosystems. In recent years, variations in plant N availability, occurring mainly due to natural and anthropogenic N deposition, have altered the trends in the terrestrial C cycle. Based on modeling studies, researchers have demonstrated that N deposition is responsible for about 0.1 - 2.3 Pg C yr -1  fixed by terrestrial vegetation (Townsend et al., 1996; Asner et al., 1997; Holland et al., 1997), which is almost half of the magnitude of C flux due to fossil fuel emission. Another factor that determines the nature of terrestrial C balance of an ecosystem is the age of the vegetation. Schimel et al. (1995) have demonstrated that forest re-growth can account for part of terrestrial C uptake by 0.5 ± 0.5 Pg C yr -1 , especially in northern mid and high 4  latitudes. This is because younger vegetation actively grows and hence sequesters more atmospheric CO2 than mature forest stands. There are many other processes that directly and indirectly affect photosynthesis and thus, the C cycle. They include land use and land cover change (Caspersen et al., 2000; Houghton and Hackler, 2006; Easter et al., 2007), reforestation (House et al., 2002; Paul et al., 2002), agricultural and grazing activities (Cerri et al., 2005), insect attack (Chapman et al., 2003; Throop et al., 2004) and invasive species (Szlavecz et al., 2006). Respiration is a process by which C is added to the atmosphere from the biosphere. There are studies that indicate that total ecosystem respiration is a major determinant of terrestrial C balances (Valentini et al., 2000). Total ecosystem respiration includes respiration by aboveground plant parts (boles, branches, twigs, and leaves) and soil respiration, which is the sum of the heterotrophic respiration, and root respiration including respiration of symbiotic microorganisms. The temporal variability of respiratory metabolism is influenced mostly by temperature and humidity conditions (Davidson and Janssens, 2006). Although ecosystem respiration has received considerable attention in recent decades, much less is known about the relative contributions of its sub- components (Jassal et al., 2007), and our understanding of how they will respond to global climate change is poor. Soil respiration (root + heterotrophic respiration) is a dominant component of C exchange in terrestrial ecosystems which accounts for more than half of the total ecosystem respiration (Black et al. 2005). This is because soils of terrestrial ecosystems contain more C than the atmosphere and live biomass together (Eswaran et al., 1993). Components of respiration can have different responses to temperature and soil water content (Boone et al., 1998; Lavigne et al., 2004), thus the effects of these environmental controls need to be understood in order to fully comprehend the soil C cycling mechanism. There are many other mechanisms that can release terrestrial C to the atmosphere. This includes both natural and anthropogenic reasons. Emission of large amounts of C to the atmosphere from vegetation can occur during forest fires (Amiro et al., 2002; Soja et al., 2004; Amiro et al., 2004) or biomass burning (Fernandez et al., 1999; Tanaka et al., 2001). These C emissions are about 1–3% of the annual C uptake by forests although their duration is very short. Forest fires and 5  biomass burning also affect the nutrient status of the soil which could have positive effects on the succeeding vegetation (Prietofernandez et al., 1993; Deluca and Sala, 2006). Another form of C flux in almost all terrestrial ecosystems is the import and export of dissolved organic C (DOC) (Neff and Asner 2001; Hornberger et al. 1994). DOC fluxes include C in the form of simple amino acids to large molecules that are transported through water flows. Fluxes of DOC into the ocean via runoff from terrestrial ecosystems are estimated to be 0.2 (Harrison et al., 2005) to 0.4 PgC yr -1  (IPCC, 2001). The strength of the terrestrial C sink was estimated to be 0.5-2.0 PgC yr -1  (Schimel et al., 1995). By sequestering atmospheric C, terrestrial ecosystems help decrease the rate of accumulation of anthropogenic CO2 in the atmosphere and its associated radiative forcing which causes anthropogenic climate change (Cihlar, 2007). Terrestrial C sinks may be responsible for taking up about one-third of all the CO2 that is released by anthropogenic activities into the atmosphere (Canadell et al., 2007). The terrestrial C sink, inferred on the basis of our current understanding, may not be permanent (Luo et al., 2003; Cox et al., 2000; Friedlingstein et al., 2003). Over the last few years there have been several studies suggesting that the size of this terrestrial C sink is vulnerable to rising near- surface temperatures (Martin et al., 1998; Nemani et al., 2002; Canadell et al., 2007). The metabolism of terrestrial ecosystems is complex and highly dynamic because ecosystems consist of coupled, non-linear processes that possess many positive and negative feedbacks (Levin, 2002; Ma et al. 2007). How the C budget of major ecosystems will respond to changes in climate is not quantitatively well understood (Baldocchi & Meyers 1998, Goulden et al., 1998; Black et al., 2000; Baldocchi et al., 2001a; Baldocchi & Wilson, 2001; Law et al., 2002; Barr et al., 2004, 2007). 1.1.2 Terrestrial water cycle There are more than one billion km 3  of water on Earth (Chahine, 1992; Schmitt, 1995; Schlesinger, 1997), covering most of the surface area (Shiklomanov, 1993).The vast majority of that water, however, is in forms unavailable to land-based or freshwater ecosystems. Globally, freshwater lakes and rivers hold 100,000 km 3 , less than one ten- 6  thousandth of all water on Earth (Jackson et al, 2001). Major stocks of water on the Earth are summarized in Table 1.1. Water molecules are transferred continually between oceans and the land surface back into the atmosphere by evaporation, falling to the land as precipitation, and transferred back to the sea by rivers and groundwater. This continual circulation is known as the water cycle. The terrestrial water cycle includes precipitation, runoff, streamflow, soil water storage and evapotranspiration (ET), the major components of which are summarized in Tables 1.1 & 1.2. Table 1.1 Major stocks of water on the Earth  Distribution Area (10 3  km 2 ) Volume (10 3  km 3 ) Percent of Total Water (%) Percent of Total Fresh Water (%) Total water 510,000 1,386,000 100 N/A World oceans 361,300 1,340,000 96.5 N/A Saline groundwater N/A 13,000 1 N/A Saline lakes 822 85.4 .006 N/A Total freshwater 149,000 35,000 2.53 100 Arctic islands 226 84 .006 .24 Antarctic glaciers 13,980 21,600 1.56 61.7 Greenland glaciers 1,800 2,340 .17 6.7 Mountain glaciers 224 40.6 .003 .12 Ground ice/permafrost 21,000 300 .022 .86 Fresh groundwater N/A 10,500 .76  30 Freshwater lakes 1,240 91 .007 .26 Wetlands 2,680 11.5 .0008 .03 In biological matter N/A 1.12 .0001 .0003 In the atmosphere  N/A 12.9 .0001 .04 Rivers (as flows on average) N/A 2.12 .0002 .006 Reservoirs* N/A 6.6 .0005 0.019 7  Source: Shiklomanov (1993) except for * from Postel (1996).  The terrestrial water balance is a key component in the global water cycle. The mass conservation equation for the land surface water balance can be written as:                                                                                                                          (   ) where ds/dt is the change of water storage with respect to time, P is precipitation, ET is evapotranspiration, and R is percolation and surface runoff. Most of the variables in the terrestrial water cycle are now observable at varying spatial and temporal resolutions and accuracy via remote sensing and in situ measurements. However attempts to close the land surface water balance equation using remotely sensed data alone have generally been unsuccessful (Tang et al., 2009). As a result, remote sensing observations need to be merged with process-based land surface models using data assimilation techniques to reasonably close the surface water balance (Eq.1.2) within an acceptable accuracy (Reichle, 2008). Table 1.2 Major components of annual terrestrial hydrological fluxes  Volume(10 3  km 3  yr -1 ) References Precipitation to the land 110 Schwarz et al. 1990 Evaporation from the oceans 425 Jackson et al, 2001 Evapotranspiration from the land 65±3 Jung et al., 2010 Evapotranspiration from the land 65.5 Oki and Kanae, 2006 Evapotranspiration from the land* 58-85 Dirmeyer et al., 2006 Annual rivers (as flows) 2.12 Shiklomanov, 1993 Annual rivers (as flows) 15 Jackson et al, 2001 Annual accessible runoff 12.5 Jackson et al, 2001 Annual total runoff 38.75 Jackson et al, 2001 8  *Multi-model comparison range 1.1.3 Terrestrial N cycle In most terrestrial ecosystems, reactive nitrogen (Nr, defined as all nitrogen species (such as NH3 and NOx) except molecular atmospheric N2) constitutes a limiting element for growth, hence the cycles of Nr and C are closely coupled (Goldstein et al., 2003). In the pre-industrial world, Nr was made available to living organisms dominantly by biological N fixation and the total amount of naturally cycled Nr between the atmosphere and the biosphere was quite small (Ayres et al., 1994). Thus, most life forms experienced intensive competition under N-limited conditions to form the natural biodiversity. By about the mid-1970s, human systems became more important than natural systems in creating Nr (Galloway et al., 1995, 2004, 2008). Currently food production accounts for ~75% of the Nr created by humans, with fossil fuel combustion and industrial uses accounting for ~13% each. Global inputs to the terrestrial N cycle have doubled in the past century due to anthropogenic activities, particularly fertilizer use and fossil fuel burning (Vitousek et al., 1997; Galloway and Cowling, 2002). 1.1.4 Coupling of the C, N and water cycles The cycling of C and N is strongly coupled to water flux through the patterns of plant growth and microbial decomposition, and this coupling creates additional feedbacks between vegetation and climate. Nr addition to terrestrial ecosystems has shown variable effects, including increases, decreases, or unchanged rates in C sequestration (Salonius and Mahendrappa, 1975; Brumme and Beese, 1992; Castro et al., 1994; Haynes and Gower, 1995; Vose et al., 1995; Bowden et al., 2000, 2004; Burton et al., 2004; Micks et al., 2004; Cleveland and Townsend, 2006; Phillips and Fahey 2007; Mo JM et al., 2006, 2008). With increasing rates of anthropogenic N deposition (Vitousek et al., 1997; Galloway and Cowling, 2002), there is a strong need to understand links between ecosystem N inputs and C sequestration. Understanding the coupled terrestrial C, N and water cycle is required to build models that project future climates and investigate the 9  role that terrestrial ecosystems play in global climate change. The importance of coupled C and water dynamics for the climate system has been increasingly recognized (Cox et al., 2000; Pielke Sr, 2001; Friedlingstein et al., 2003; Seneviratne et al., 2006; Betts et al, 2000, 2007a,b), and much progress has also been made in gaining insight into the coupling between C, water and N cycles across a range of time and spatial scales (Pielke Sr, 2001; Friedlingstein et al., 2003; Seneviratne et al., 2006; Betts et al, 2007a,b; Baldocchi, 2008), however, the mechanisms behind these coupled cycles are not well understood, especially at the landscape scale. 1.2 State of approaches for studying coupled terrestrial C and water dynamics 1.2.1 State of field-scale monitoring and modeling approaches A variety of approaches are being used in C, N and water cycles studies, including direct field-scale observations (e.g. biometric measurement, soil chamber efflux measurements), eddy-covariance (EC) flux tower and high-precision CO2 concentration tower measurements, targeted field campaigns (e.g. using aircraft or tethered balloon measurement platforms), and satellite remote sensing. These datasets are then coupled using a wide variety of modeling approaches including: (i) bottom-up models (i.e. process-based ecosystem models, hybrid models that adopt a combination of empirical and process-based algorithms); (ii) atmospheric inversion modeling (Gurney et al., 2002); and (ii) data assimilation systems which combine diverse data and models into a unified description of a physical/biogeochemical system consistent with observations (Levy et al., 1999). An example of a recent data assimilation system is CarbonTracker, which uses a state-of- the-art atmospheric transport model coupled to an Ensemble Kalman filter (Peters, et al., 2009). The major monitoring and modeling methods are evaluated and compared in Chapter 2; while the arrays of airborne and satellite sensors developed for studying the coupled terrestrial C and water cycles are reviewed and evaluated below. 10  1.2.2 Arrays of airborne and satellite sensors developed to monitor coupled C and water cycles Remote sensing is the observation of a phenomenon from a distance, using devices that detect electromagnetic radiation (Sass and Greed, 2011). Satellite-borne remote sensing offers unique opportunities to parameterize land surface characteristics over large spatial extents at variable spatial and temporal resolutions. Since 1960s, there has been a substantial increase in the number of satellite and airborne sensors for Earth observations that cover a large range of the electromagnetic radiation spectrum. Scientists began applying these observations to vegetation studies from the 1970s (Rouse et al., 1974; Tucker et al., 1979). For example, Tucker et al. (1986) exploited the properties of chlorophyll pigments to absorb wavelengths in the red spectral region and structural properties of leaves to reflect near-infrared spectra using imagery acquired by the Advanced Very High Resolution Radiometer (AVHRR) sensor onboard TIROS (Television Infrared Observation Satellite). This pioneering study for the first time allowed a synoptic view (1 – 5 km) of the coupled atmosphere-biosphere system (Tuker et al., 1986) and opened possibilities for global perspectives in ecology. The first Landsat satellite, launched in 1972, carried the Multispectral Scanner System (MSS) sensors which were specifically designed to map land resources with finer spatial resolution (68 m × 82 m) than the AVHRR. The program was the first civil, non-meteorological satellite program and provided observations for any place on Earth once every 18 days, offering a wide range of studies on terrestrial vegetation and C and water cycles. Landsat Thematic Mapper sensors carried onboard the later series of Landsat satellites, acquire images at a 30-m spatial resolution with a 16-day interval. The acquired data have been the backbone for land-cover, vegetation and C cycle studies. The National Aeronautics and Space Administration’s (NASA’s) Earth Observing System (EOS) launched in 1999 (Tilford S. 1984), brought new capabilities for monitoring vegetation productivity and other properties with near-daily and global coverage. Multispectral sensors---Moderate Resolution Imaging Spectroradiometer (MODIS), onboard the Terra and Aqua platforms, have produced valuable global observation datasets for C and water cycles research since 11  the early 2000s. MODIS provides a global coverage every 1-2 days with 36 bands. The spatial resolution of MODIS (pixel size at nadir) is 250 m for channels 1 and 2 (0.6µm - 0.9µm), 500 m for channels 3 to 7 (0.4µm - 2.1µm) and 1000 m for channels 8 to 36 (0.4µm - 14.4µm), respectively. Data from the satellite-borne MODIS are currently used in the calculation of global weekly GPP and ET at 1-km spatial resolution (Running et al., 2004). Sensors that have potential applications in C and hydrology studies fall into two groups: optical and microwave. Optical sensors cannot penetrate vegetation or clouds.  In contrast, microwave sensors are able to penetrate vegetation and can collect data independently of cloud cover and solar illumination. There are two types of microwave sensors: active sensors and passive sensors. The former send and receive their own energy; while the latter detect the microwaves emitted by the Earth’s surface. The microwave bands, useful for vegetation and C and water cycles, are K, X, C, and L, ranked in increasing wavelengths. K- and X-bands are useful for detecting surface temperature, snow density, and rainfall rates, whereas C- and L-bands are sensitive to soil moisture (Sass and Greed, 2011). Microwave wavelengths penetrate greater depths into plant canopies than optical sensors (Kasischke et al. 1997). The potential for using RADAR (RAdio Detection And Ranging) to study terrestrial C and water cycles, particularly for assessing standing woody biomass is promising. The sensitivity of RADAR to vegetation biomass strongly depends on wavelength. Single-band RADAR is able to detect aboveground biomass up to approximately 100 Mg per hectare (Dobson et al. 1992, Luckman et al. 1998). In addition, multiband RADAR enables researchers to separate biomass into component fractions (e.g. stem and canopy) (Saatchi and Moghaddam 2000). Synthetic aperture RADAR (SAR) is also sensitive to vegetation structure and to the amount of biomass, including both photosynthetic (green) and nonphotosynthetic vegetation components (Turner et al., 2004). 12  LiDAR (Light Detection and Ranging) is a remote sensing technology that determines distances to an object or surface using laser pulses, which is a relatively new technology compared to optical sensors, and has the added capability of characterizing the distribution of foliage with height in the canopy (Lefsky et al. 2002, Treuhaft et al. 2002, 2004; Turner et al., 2004). LiDAR data have proved to be highly effective for determining three dimensional forest attributes including leaf area index (LAI) and the probability of canopy gaps within different layers of canopy (Coops et al., 2004, 2007c). Interpreted LiDAR data have been further used for landscape C modeling and scaling studies (Hilker et al. 2008b; Chen et al., 2009c). 1.3 Hypothesis and objectives 1.3.1 Knowledge and scale gaps Progress in C balance studies has been achieved at the extreme ends of the spatial scale spectrum, either at local scales (less than 1-3 km 2 , e.g. EC-measurements) or at large continental scales (larger than 10 6  km 2 , e.g. global inverse modeling). However, methods to estimate and assess the land-atmosphere C fluxes at the intermediate scales (i.e. landscape to regional scales, finer than those used in global inversions and larger than local EC flux measurements and roughly defined as the range between 10 and 10 6  km 2 ) are notably lacking and there are significant uncertainties in C balance estimates at these scales. It is not appropriate extremely unreliable to upscale stand-level fluxes (i.e. EC measurements) to a region by simple spatial extrapolation and interpolation because of the heterogeneity of the land surface and the nonlinearity inherent in ecophysiological processes (Levy et al., 1999). It is also challenging to apply atmospheric inversion techniques to regional scales for quantifying annual C budgets because at such intermediate scales the atmosphere is often poorly constrained (Gloor et al., 1999; Matross et al., 2006). Moreover, aggregation errors and errors in atmospheric transport, both within the planetary boundary layer (PBL) and between the PBL and free troposphere, can also be obstacles to obtaining quantitative estimates of regional C fluxes (Lin et al., 2006). Hence, there is a strong motivation to develop methods to quantify and 13  validate estimates of the C balance at these intermediate scales (Lin et al., 2006; Bakwin et al., 2004; Matross et al., 2006, Chen et al., 2007c). Reliable estimates of terrestrial C sources and sinks at landscape to regional spatial scales are fundamental to improving our understanding of the C cycle (Crevoisier et al., 2006). As discussed above, to date, the two major weakness areas in the study of coupled terrestrial C, N and water cycles are: (i) in nature, hydrological and biogeochemical processes are closely linked; however, the coupled processes involved in terrestrial C, N and hydrological cycles are far from well understood; and (ii) there are significant uncertainties in estimates of the components of the C balance, especially at landscape and regional scales. 1.3.2 Hypotheses The two fundamental hypotheses of this work are: (i) a synthesis aggregation methodology --- integrating remotely sensed data, local-scale measurements, N fertilization data, process-based ecosystem modeling and data-model assimilation, is a pragmatic approach towards a better understanding of the coupled C, N and water dynamics at landscape/watershed scales; and (ii) landscape- and regional-scale C fluxes can be reasonably estimated using a multi-tiered approach involving direct land surface and remote-sensing measurements, high-precision CO2 concentration observations and ecosystem-, footprint- and inversion- modeling. 1.3.3 Objectives The overarching goal of this thesis is to test the above-mentioned hypotheses, that is, to provide an improved quantitative understanding of coupled C, N and water cycles in temperate and boreal forests at intermediate scales (local to regional). The specific objectives in this research are fourfold: 1) to identify key processes and controls on within-site and cross-site variability of C, water and energy fluxes in space and time (from instantaneous to inter- 14  annual); 2) to develop state-of-the-art inverse analyses and data-model fusion techniques to advance parameter optimization methods in a common ecosystem modeling framework and to study the effects of N addition on C sequestration; 3) to assess and characterize spatial representativeness of flux network measurements based on  footprint and geospatial analyses and to upscale local EC C fluxes to the landscape level; and 4) to assess land surface impacts on tower CO2 concentration measurements and to derive regional C fluxes using an integrated approach including concentration footprint modeling and PBL budgeting. The first two objectives were performed in a coastal temporal forest landscape encompassing the British Columbia (BC) Flux Station (three flux towers); the third objective was achieved using data obtained from Fluxnet Canada Research Network (FCRN) and the Canadian Carbon Program (CCP) (2002-2011); while the last objective was achieved using data obtained in a boreal region centred at Canada’s only tall tower, which is located at East Trout Lake, Saskatchewan. 1.4 Study areas and data descriptions This thesis research was conducted at three scales associated with the different objectives, which are local (stand), landscape, and regional scales (Figure 1.2). 1.4.1 Stand-level research The stand-level research associated with the first two objectives, took place at the BC flux station area (49º30 Ń - 49º55 Ń, 124º50 Ẃ - 125º30 Ẃ), located on the east coast of Vancouver Island, BC between Denman Island and Campbell River. This area is in the very dry maritime subzone (CWHxm) of the coastal western hemlock (CWH) biogeoclimatic zone (Meidinger and Pojar 1991). The BC Flux Station contains three EC towers (DF49, HDF88, and HDF00) with different ages of Douglas-fir (Pseudotsuga menzesii) dominated forest with similar soil and weather, which are expected to represent 15  the typical ages of forest found in this region after harvesting, site preparation, and planting. In 2010, stand age adjacent to each tower was approximately 62 (DF49), 22 (HDF88), and 11 (HDF00) years old, respectively. The three sites, DF49, HDF88 and HDF00 are on a 5-10º north-east slope, a 2-5º south-east slope, and on a level gravelly fluvial terrace, respectively, which results in spatial variation in soil moisture and plant associations within each site. The EC tower and chamber flux measurements and a suite of climate measurements made at these sites available from 1997 were used in this thesis study (Morgenstern et al. 2004; Humphreys et al. 2003, 2006, Chen et al. 2009b, 2010). Objectives 1 & 2 were mostly implemented at the DF49 flux tower site located 10 km southwest of Campbell River and 9 km west of Georgia Strait. The site is on glaciated terrain with a 5-10º northeast-facing slope.  LiDAR data were acquired for a 12 km2 area around the DF49 tower covering EC footprint area on June 8th 2004, by Terra Remote Sensing (Sidney, British Columbia, Canada) with a spacing density of 0.7 hits per m2 and a footprint (spot size) of 0.19 m. Separation of vegetation and terrain was carried out using Terrascan v.4.006 (Terrasolid, Helsinki, Finland) to iteratively classify LiDAR data into either ground or non-ground returns (Hilker et al., 2008a). The elevation of the DF49 tower area covered by LiDAR varies from about 470 m in the southwest to 111 m in the northeast, with the flux tower at 330 m. 16   Figure 1.1 Study areas.  A, Stand scale: (a) pictures for the DF49 EC flux tower site of British Columbia flux station, with Douglas-fir dominated forest. In 2010 the stand ages adjacent to the tower was approximately 62 years old.  DF49 was harvested in December 2010; B, Landscape scale: (b) a land cover map of a 6 × 6 km area centered at the DF40 EC tower at a 30-m resolution, the white line is the 90% annual footprint climatology 17  contour at the DF49 tower in 2004; (c) British Columbia flux station, which includes three EC towers (DF49, HDF88, and HDF00) with different ages of Douglas-fir dominated forest and in 2010 the stand ages adjacent to each tower were approximately 62, 22, and 11 years old, respectively; (d) the eco-region map of Canada’s landmass and with the locations of high-precision CO2 concentration towers in the long term observations network for atmospheric measurements of CO2; C, Regional scale: (e) a 1000 × 1000 km region surrounding two high-precision CO2 concentration towers: East Trout Lake (ETL) 106-m tall tower, Saskatchewan, located at (54°21'N, 104°59'W) and Candle Lake 28-m high tower (CL), located at (53°59'N, 105°07'W); while triangles are the cluster of EC towers within the tower concentration footprint area. Land cover data at 1 km resolution are from the Global Land Cover Map 2000 (GLC2000); (f) picture of the ETL 106-m tall tower with CO2 concentration measurements at four levels. 1.4.2 Landscape-level research To upscale local EC C fluxes to the landscape level, this study developed an approach for evaluating the representativeness of EC tower flux measurements and assessing sensor location bias based on footprint modeling and remote sensing (objective 3). This approach was applied to the BC Flux Station (three sites: DF49, HDF88 and HDF00, Figure 1.1b & 1.1c) and then to the 12 main sites of the FCRN/CCP located along an east-west continental-scale transect, covering grassland, forest, and wetland biomes (Figure 1.3).  The site characteristics and tower measurement information for the12 CCP main flux towers are shown in Table 1.3. Year-round half-hourly flux and meteorological measurements made in 2006, which was generally considered a normal weather year over most of the sites, were used in this study. Missing flux and meteorological data were filled using the gap-filling method of Chen et al. (2009b). Those gap-filled datasets were used to drive the footprint model for assessing the monthly and annual footprint climatology for all the sites.  Landsat imagery acquired in 2006, or as close as possible, for each of the 12 sites within a 6   6 km area centred at each tower location (downloaded from the U.S. Geological Survey at http://glovis.usgs.gov/).   18  1.4.3 Regional C flux research The regional scale research (objective 4) was achieved in a 1000 × 1000 km region surrounding the East Trout Lake (ETL) 106-m-tall high-precision CO2 concentration tower, Saskatchewan (54°21'N, 104°59'W) (Figure 1.1). This study area is centred at the ETL tower (Figure 1.1). The ETL monitoring site is located near the mouths of the rivers draining to the East Trout and Nipekamew Lakes. The terrain in the surrounding area is slightly rolling with vegetation primarily consisting of black spruce and a small amount of deciduous larch and trembling aspen (Figure 1.1). Meteorology (air temperature, humidity, wind speed and direction) and CO2 concentration have been measured at 4 levels (95, 55, 33 and 22 m) on a 106-m high SaskTel communication tower.  Another high-precision CO2 concentration tower --- the Candle Lake  (CL, 53°59'N, 105°07'W) monitoring site, was located 40 km southwest of the ETC tower, near the southern edge of the boreal forest approximately 100 km NE of Prince Albert, Saskatchewan.  The CL CO2 concentration sensors are mounted 28 m above the ground on the same EC flux tower, which is referred to as Southern Old Black Spruce (SOBS) and was established in 1994 as part of the Boreal Ecosystems-Atmosphere Study (BOREAS) (Sellers et al., 1997). The greenhouse gas measurement program began at Candle Lake in June 2002. The landscape in the region is predominantly flat, with slight topographical undulations. This site has an elevated water table and is generally wet. Ground level vegetation consists mostly of moss and Labrador tea. Further site details are found in Jarvis et al. (1997), Griffis et al. (2003), and Kljun et al. (2006). 19   Figure 1.2  Locations of flux towers in the Fluxnet Canada Research Network (FCRN) / Canadian Carbon Program (CCP) (2002-2011) and eco-region distributions of Canada’s landmass. For description of flux sites, see Table 1. 3.  Table 1.3 Characteristics and tower measurement information for the CCP sites a) a)  Site names: First two letters indicate the province (MB – Manitoba, QC – Quebec, SK – Saskatchewan, NB – New Brunswick, BC – British Columbia, ON – Ontario, AB – Alberta); NOBS – Northern Old Black Spruce, EOBS – Eastern Old Black Spruce, OJP – Old Jack Pine, SOBS – Southern Old Black Spruce, OBF – Old Balsam Fir, DF49 – Douglas fir (established 1949), WPP39 – White Pine Plantation (established 1939), OA – Old Aspen, GRS – Grassland, OMW – Old mixed wood, WPL – Western peat land, EPL – Eastern peat land).  Site a)  Vegetation type Eco-region Latitude ( o N  ) Longitude ( o W) References MB- NOBS Coniferous boreal forest Boreal shield 55.87956 98.48084 Dunn et al. (2007) and Turner et al .(2003) QC-EOBS Coniferous boreal forest Boreal shield  49.69247  74.34204 Richardson et al. 2006 SK-OJP Coniferous boreal forest Boreal plains 53.91634 104.69203 Barr et al. (2006) SK-SOBS Coniferous boreal forest Boreal plains 53.98717  105.11779 Barr et al. (2006) NB-OBF Coniferous maritime forest Atlantic maritime 46.47388 67.09937 Xing et al. (2007) BC-DF49 Coniferous temperate rainforest Pacific maritime 49.86883  125.33508 Chen et al (2009b) ON- WPP39 Coniferous temperate forest Mixedwood plains 42.71222  80.3572 Arain and Restrepo-Coupe (2005) and Peichl et al. (2010) SK-OA Deciduous boreal forest Boreal plains 53.62889  106.19779 Barr et al. (2006) AB-GRS Short/mixed grass prairie (C3); grassland Prairies 49.70919 112.94025 Flanagan and Johnson (2005), Flanagan et al. (2002, 2009) and  Richardson et al. (2006) ON-OMW Mixedwood temperate forest Boreal shield 48.21738 82.15553 McCaughey et al. (2006) AB-WPL Treed fen; wetland Boreal plains 54.95384 112.46698 Flanagan (2009) and Syed et al. (2006) ON-EPL Ombrotrophic bog; wetland Mixed wood plains 45.40940  75.51870 Admiral et al. (2007) and Peichl et al. (2010) 21  1.5 Thesis structure The body of the thesis consists of 7 Chapters as outlined below, covering three phases of research, which are (i) improvement of our mechanistic understanding of the coupled C, N and water at the stand level based on decadal EC data and the data-model fusion approach (Chapters 3 & 4); (ii) landscape up-scaling strategy for C and water fluxes based on EC tower measurements and remote sensing (Chapters 5, 6 & 7); and (iii) retrieving regional C flux from tall-tower CO2 concentration measurements and remote sensing-based footprint integration (Chapter 8). Chapter 2 distils, synthesises and extracts information from the rapidly growing literature on C, N and water cycles across local to global spatial scales and over a range of time scales. Based on perceived current research gaps in this field, this Chapter highlights recent advances, discusses current research trends and near-future directions and proposes an up- scaling framework for landscape and regional C and water fluxes estimates. The whole thesis research framework and content are designed on the basis of the results presented in Chapter 2. Chapter 3 presents research on the dynamics of net C dioxide exchange between the biosphere and atmosphere (triggered by such critical features as switches, pulses, lags, and acclimation) based on the EC data carried out in a Pacific Northwest Douglas-fir forest (60- year old in 2009) on the east coast of Vancouver Island, Canada. This Chapter advances our mechanistic understanding of the biophysical controls of the ecosystem C budget and on how trends and inter-annual variations in climate affect C and water exchange between terrestrial ecosystems and the atmosphere on a decadal time scale. Chapter 4 presents an approach using a data-model fusion technique which encompasses both model parameter optimization and data assimilation, minimizes the effects of inter- annual climatic perturbations, and focuses on the biotic and abiotic factors controlling seasonal C fluxes using a pre-fertilization 9-year-long time series of EC data (1998-2006) to discern N fertilization impacts on C sequestration. 22  Chapter 5 presents a pragmatic and reliable method to examine the influence of patch-scale heterogeneities on the uncertainty in long-term EC flux data and to upscale the EC C fluxes to the landscape scale. A footprint climatology model was developed and was applied to three Douglas-fir stands (BC Flux Station). Chapters 6 & 7 extend the methodology from Chapter 5 to the Fluxnet-Canada Research Network for assessing the location bias and characterizing spatial representativeness of the EC flux tower measurements. Chapter 8 presents a case study on concentration footprint climatology analyses.  A Simple Analytical Footprint model based on Eulerian coordinates for the scalar concentration (SAFE-C) was applied to the two high-precision CO2 mixing ratio towers in Saskatchewan, Canada.The concentration footprint climatologies for these two towers were analyzed and the impacts of the complex heterogeneous distribution of land cover types over regions surrounding the towers on the CO2 mixing ratio measurements were assessed.  23  Chapter  2:  Understanding of coupled terrestrial carbon, nitrogen and water dynamics—an overview2 2.1 Synopsis Coupled terrestrial carbon, nitrogen and hydrological processes play a crucial role in the climate system, providing both positive and negative feedbacks to climate change. Chapter 2 summarizes published research results to gain an increased understanding of vegetation - atmosphere dynamics. A variety of methods, including monitoring (e.g. EC flux tower, remote sensing, etc.) and modeling (i.e. ecosystem, hydrology and atmospheric inversion modeling) of terrestrial carbon (C) and water fluxes are evaluated and compared. Chapter 2 highlights two major research areas where additional research could be focused: (i) conceptually, the hydrological and biogeochemical processes are closely linked; however, the coupled terrestrial C, nitrogen (N) and hydrological processes are not well understood; and (ii) there are significant uncertainties in estimates of the components of the C balance, especially at landscape and regional scales. To address these two issues, a synthetic research framework is needed which includes both bottom-up and top-down approaches integrating scalable (footprint and ecosystem) models and a spatially-nested hierarchy of observations including multispectral remote sensing, inventories, existing regional clusters of EC flux towers and CO2 mixing ratio towers and chambers. 2.2 Introduction The terrestrial biosphere plays a crucial role in the climate system providing both positive and negative feedbacks to climate change (Treut et al., 2007). The terrestrial C cycle is closely linked to hydrological and nutrient controls on vegetation (Betts et al., 2000; Cox et  2  Reprinted from Sensors, Vol. 9(11), 8624-8657, [B. Chen] and N.C. Coops (2009), Understanding of coupled terrestrial carbon, nitrogen and water dynamics — an overview. Sensors, 9(11), 8624-8657, (doi:10.3390/s91108624) Copyright (2009), with permission from MDPI Publishing, Basel, Switzerland. 24  al., 2000). Understanding the coupled terrestrial C, N and water cycles is required to gain a comprehensive understanding of the role that terrestrial ecosystems play in the global climate change. Much progress has been made in gaining insight of the coupled processes between C, N and water cycles across a range of time and spatial scales (Pielke et al., 2001; Friedlingstein et al., 2003; Seneviratne et al., 2006; Betts, 2007a; Baldocchi, 2008). Since the early 1990s, there has been an increased interest in monitoring of the CO2, water vapor and energy exchanges between the atmosphere and terrestrial ecosystems by a variety of methods, such as the eddy-covariance techniques, satellite and airborne remote sensing, CO2 concentration and isotope measurements. Meanwhile, there are various kinds of models have been developed to better understanding of these processes and for large-scale C and water budgeting. The large number of papers published since the 1980s on the terrestrial C and water cycles have resulted in the publication of several major reviews from different perspectives. For example, Running et al. (2004) described a blueprint for more comprehensive coordination of the various flux measurement and modeling activities into a global terrestrial monitoring network by reviewing the literature published before the middle 1990s. Baldocchi (2008) recently provided a comprehensive review of research results associated with a global network of C flux measurements systems. The topics discussed by this review include history of the network, errors and issues related with the EC method, and a synopsis of how these data are being used by ecosystem and climate modellers and the remote-sensing community (Baldocchi, 2008). Kalma et al. (2008) reviewed satellite-based algorithms for estimating evepotranspiration (ET) and land surface temperatures at local, regional and continental scales, with particular emphasis on studies published since the early 1990s; while Verstraeten et al. (2008) provided a comprehensive review of remote-sensing methods for assessing ET and soil moisture content across different scales based on the literature published after 1990s. Marquis and Tans (2008) reviewed satellite-based instruments on CO2 concentration measurements. 25  In this Chapter, we distil and synthesize the rapidly growing literature on C and water cycles across local to global spatial scales and over a range of time scales. To give the reader a perspective of the growth of this literature, a search of Web of Science produced over 2000 papers with the key words ‘ecosystem carbon, water and N cycles’ published since 1990 which is indicative of the large amount of research recently being undertaken on these topics. In order to filter through this large body of literature, we concentrate on papers discussing on the coupling processes between C, water and N cycles and we extract information from a database of published results that we have collated during the past decade (available on request). In terms of content, the Chapter covers the state of knowledge, monitoring and modeling of the coupled terrestrial C and water cycles. Our aim is to highlight the recent advances in this field, and propose areas of future research based on perceived current gaps in the literature. The review is divided into several inter-connected sections. First, we review the scientific background of the linkage between terrestrial ecosystems and climate, and revise the state of knowledge on terrestrial C cycling and coupling of C and N cycles. Second, we discuss the ground-based and satellite-based monitoring methods and observation networks associated with measuring C and water fluxes, CO2 concentration and C isotopes. Third, we report on the recent advances in modeling approaches associated with the terrestrial biochemical and hydrological studies. Fourth, we discuss research gaps in C sinks/sources estimates and finally, we discuss current research trends and near-future directions in this field and propose an up-scaling framework for landscape and regional C and water fluxes estimates. 2.3 Scientific background and state of knowledge 2.3.1 Overview of terrestrial ecosystems and climate The climate system is controlled by a number of complex coupled physical, chemical and biological processes (Figure 2.1). The terrestrial biosphere plays a crucial role in the climate system, providing both positive and negative feedbacks to climate change through biogeophysical and biogeochemical processes (Treut, 2007). Couplings between the climate 26  system and biogeochemistry are mainly through tightly linked dynamics of C and water cycles. The importance of coupled C and water dynamics for the climate system has been increasingly recognized (Cox et al., 2000; Pielke Sr, 2001; Friedlingstein et al., 2003; Seneviratne et al., 2006; Betts et al, 2000, 2007a, b); however the mechanisms behind these coupled cycles are still far from well understood.   Figure 2.1 Schematic view of the components of the climate system, their processes and interactions (Treut, 2007). 2.3.2 Terrestrial C cycling One of the crucial issues in the prognosis of future climate change is the global budget of atmospheric CO2. The growth rate of atmospheric CO2 is increasing rapidly. Three processes contribute to this rapid increase: fossil fuel emission, land use change (deforestation), and ocean and terrestrial uptake. As shown in Figure 2.2, terrestrial C budgets have large uncertainties and inter-annual variability.  27   Figure 2.2 Global CO2 budget from 1959 to 2006. Upper panel (A): CO2 emissions to the atmosphere (sources) as the sum of fossil fuel combustion, land-use change, and other emissions. Lower panel (B): The fate of the emitted CO2, including the increase in atmospheric CO2 plus the sinks of CO2 on land and in the ocean (Canadell et al., 2007).  Terrestrial ecosystems mediate a large part of CO2 flux between Earth’s surface and the atmosphere, with ~120 Pg C yr -1  taken up by photosynthesis and roughly the same amount released back to the atmosphere by respiration annually (Treut et al., 2007; Prentice et al., 2001). Imbalances between gross ecosystem photosynthesis or gross primary productivity (GPP) and ecosystem respiration (Re) lead to land surfaces being either CO2 sinks or sources. The magnitudes of sinks and sources have fluctuated on annual and longer time scales due to variable climate, land use change, disturbance, and changes in the age distribution and species composition of ecosystems (Battle et al., 2000; Arain et al., 2002; Law et al., 2002; Morgenstern et al., 2004; Humphreys et al., 2005, 2006; Urbanski et al., 2007). Terrestrial ecosystems modify atmospheric C balance through many mechanisms. A detailed understanding of the interactive relationships in atmosphere–biosphere exchange is relevant 28  to ecosystem-scale analysis and is needed to improve our knowledge of the global C cycle (Falk et al., 2008). In recent years, scientists have learnt that vegetation, soil (Melillo et al., 1989; Knapp et al., 1993) and animals (Naeem et al., 1995; Hattenschwiler and Bretscher, 2001) play key roles in mediating the terrestrial C cycle. Plants being the primary producers, it is from them that mass and energy gets transformed to other living organisms (Engel and Odum, 1999) within an ecosystem. The process of photosynthesis fixes atmospheric C into the biosphere. Atmospheric C enters the plant through stomatal openings that are controlled by a variety of environmental factors (Jarvis, 1976; Griffis et al., 2003). These factors include ambient temperature, atmospheric CO2 concentration, nutrient availability, soil water availability and forest age (Schimel, 1995; Prentice et al., 2001). Changes in the atmospheric CO2 concentration and the corresponding changes in the climate have altered the magnitudes of terrestrial C cycling. For example, climate change induced increases in vegetation growth due to earlier springs and lengthened growing seasons were detected by the phase shift of seasonal atmospheric CO2 cycle by Keeling et al. (1996) and satellite-based vegetation index analysis by Myneni et al. (1997). Studies indicate that an increase in atmospheric CO2 enhances photosynthesis (e.g. Woodward and Friend, 1988) and hence increases assimilation of atmospheric CO2 by the terrestrial vegetation. N availability to plants is another factor that can affect photosynthesis. This is because N is a primary nutrient for plant growth. In the recent years, variations in plant N availability have also altered the trends in the terrestrial C cycles. Variations in plant N availability occur mainly due to natural and anthropogenic N- deposition. Based on modeling studies, Townsend et al. (1996), Asner et al. (1997) and Holland et al. (1997) have demonstrated that N deposition is responsible for about 0.1-2.3 Pg C yr -1  fixed by terrestrial vegetation. Another factor that determines the nature of terrestrial C balance of an ecosystem is the age of the vegetation. Schimel et al. (1995) have demonstrated that forest re-growth can account for part of terrestrial C uptake as much as 0.5 ± 0.5 Pg C yr - 1 , especially in northern mid and high latitudes. This is because younger vegetation actively grows and hence sequesters more atmospheric CO2 as opposed to mature forest stands. There are many other processes that directly and indirectly affect photosynthesis and thus, the C 29  cycle. They include land use and land cover change (Caspersen et al., 2000; Houghton and Hackler, 2006; Easter et al., 2007), reforestation (House et al., 2002; Paul et al., 2002), agricultural and grazing activities (Cerri et al., 2005), insect attack (Chapman et al., 2003; Throop et al., 2004) and invasive species (Szlavecz et al., 2006). Respiration is a process by which C is added to the atmosphere from the biosphere. There are studies that indicate that total ecosystem respiration is a major determinant of terrestrial C balances (Valentini et al., 2000). Total ecosystem respiration includes respiration by aboveground plant parts (boles, branches, twigs, and leaves) and soil respiration, which is the sum of the heterotrophic respiration, and root respiration including respiration of symbiotic microorganisms. The temporal variability of respiratory metabolism is influenced mostly by temperature and humidity conditions (Davidson and Janssens, 2006). Although ecosystem respiration has received considerable attention in recent decades, much less is known about the relative contributions of its sub- components (Jassal et al., 2007), and our understanding of how they will respond to global warming is poor. Soil respiration (root + heterotrophic respiration) is a dominant component of C exchange in terrestrial ecosystems which accounts for more than half of the total ecosystem respiration (Black et al. 2005). This is because soils of terrestrial ecosystems contain more C than the atmosphere and live biomass together (Eswaran et al., 1993). Components of respiration can have different responses to temperature and soil water content (Boone et al., 1998; Lavigne et al., 2004), thus the effects of these environmental controls need to be understood in order to fully comprehend the soil C cycling mechanism. There are many other mechanisms that can release terrestrial C to atmosphere. This includes both natural and anthropogenic reasons. Emission of large amounts of C to the atmosphere from vegetation can occur during forest fires (Amiro et al., 2002; Soja et al., 2004; Amiro et al., 2004) or biomass burning (Fernandez et al., 1999; Tanaka et al., 2001). These C emissions are of very high magnitudes although their duration is very short. Forest fires and biomass burning also affect the nutrient status of the soil which could have positive effects on the succeeding vegetation (Prietofernandez et al., 1993; Deluca and Sala, 2006). Another form of C flux in almost all terrestrial ecosystems is the import and export of dissolved organic C (DOC) (Neff and Asner 2001; Hornberger et al. 1994). DOC fluxes include C in 30  the form of simple amino acids to large molecules that are transported through water flows. Fluxes of DOC into the ocean via runoff from terrestrial ecosystems are estimated to be 0.2 (Harrison et al., 2005) to 0.4 Pg C yr -1  (IPCC, 2001,Thompson et al., 1996; Chen et al., 2003; Canadell et al., 2007; Schulze, 2006). The strength of the terrestrial C sink was estimated to be 0.5-2.0 Pg C yr -1  (Schimel et al., 1995). By sequestering atmospheric CO2, the terrestrial ecosystems help decrease the rate of accumulation of anthropogenic CO2 in the atmosphere, and its associated climate change (Cihlar, 2007). Terrestrial C sinks may be responsible for taking up about one-third of all the CO2 that is released into the atmosphere (Canadell et al., 2007). The terrestrial C sink, inferred based on our current understanding, may not be permanent (Luo et al., 2003; Cox et al., 2000; Friedlingstein et al., 2003). Over the last few years there have been several studies suggesting that the size of this terrestrial C sink is vulnerable to global warming (Martin et al., 1998; Nemani et al., 2002; Canadell et al., 2007). The metabolism of terrestrial ecosystems is complex and highly dynamic because ecosystems consist of coupled, non-linear processes that possess many positive and negative feedbacks (Levin, 2002; Ma et al. 2007). How the C budget of major ecosystems will respond to changes in climate is not quantitatively well understood (Baldocchi & Meyers 1998, Goulden et al., 1998; Black et al., 2000; Baldocchi et al., 2001a; Baldocchi & Wilson, 2001; Law et al., 2002; Barr et al., 2004, 2007). A detailed understanding of the interactive relationships in atmosphere-biosphere exchange is relevant to ecosystem-scale analysis and is needed to improve our knowledge of the global C cycle (Falk,M et al., 2008). The metabolism of terrestrial ecosystems is complex and highly dynamic because ecosystems consist of coupled, non-linear processes that possess many positive and negative feedbacks (Levin et al., 2002; Ma et al., 2007). Complex features of ecosystem metabolism are relatively unknown and how the C budget of major ecosystems will respond to changes in climate is not quantitatively well understood (Black et al., 2000; Baldocchi et al., 2001a,b; Barr et al., 2004; Law et al., 2002).   31  2.3.3 Coupling of the C and water cycles Thermodynamically, a terrestrial ecosystem is an open system. Therefore, hydrological and C cycles are closely coupled at various temporal and spatial scales (Betts et al., 2007a; Ball et al., 1987; Levis et al., 1999; Rodriguez-Iturbe and Ecohydrology, 2000; Joos et al., 2001; Arain et al., 2006; Balanken et al., 2004; Snyder et al., 2004). C uptake for example, is closely coupled to water loss by ecosystems mainly through leaf stomatal pathway governed principally through leaf conductance (Javis, 1976; Harris et al., 2004; Rodriguez-Iturbe et al., 2001). Soil organic C decomposition is very sensitive to soil moisture content via microbial activity and other processes (Betts, 2007a; Levis et al., 1999; Snyder et al., 2004; Parton et al., 1993a, b; D'odorico et al., 2004). The flux of terrestrial organic C by river runoff to the ocean and wetland discharge is an important component of the global organic C cycle (Hedges, 1992; Wang et al., 2004). It is estimated that 0.25 Pg C DOC is discharged to the ocean by the rivers each year (Meybeck, 1982). The land surface hydrological processes (in particular the terrestrial river systems) play an important role in transport of dissolved and particulate organic C from terrestrial to marine ecosystems (Wang et al., 2004). However, the interactions between C and water cycles and the mechanisms of how these interactions will shape future climatic and biosphere conditions are far from well understood. 2.3.4 Coupling of the C and N cycles N is an essential element that can limit the growth of living organisms across a wide range of ecosystems (Vitousek and Howarth 1991; Allison et al., 2008). However, global inputs to the terrestrial N cycle have doubled in the past century due to anthropogenic activities, particularly fertilizer use and fossil fuel burning (Vitousek et al., 1997; Galloway and Cowling, 2002). N depositions have changed the dynamics of N and C in many ecosystems of developed countries (e.g. in Europe and North America) (Brumme and Khanna, 2008), with the impact in other Asian developing countries projected to increase over the next few decades (Galloway and Cowling, 2002). N addition to forest and other ecosystems will affect the health and vitality of ecosystems and the terrestrial C cycle (Aber et al., 1989, 1998). There is a continuing discussion on whether N deposition has positive (e.g. increases in forest 32  growth and C sequestration) or negative (e.g. increases in nitrate leaching, reduction in forest growth and C sequestration) effects on an ecosystem which depends on the N status of the system and the rate and duration of N deposition (Aber et al., 1989, 1998). Bauer et al. (1996) summarized that (1) if the remaining available N does not exceed the capacity for N uptake by vegetation net primary productivity (NPP), C sequestration may be enhanced; (2) if deposition rates exceed the capacity for N uptake, nutrient imbalances can lead to forest decline due to N saturation (Agren and Bosatta, 1988; Nihlgard (1985) and Aber  et al. (1989,1998) hypothesized that forest ecosystems positively respond to N addition in the short-term, while negatively respond to N addition over the long-term. N fertilization in northern temperate zones has been estimated to enhance C storage by 0.3–0.5 Pg C year-1 (Townsend, 1996). Pregitzer et al. (2005) reported that simulated chronic N deposition has increased C storage in northern temperate forests. Olsson et al. (2005) found that fertilization of a boreal Norway spruce stand led to a three-fold increase in aboveground productivity, possibly due to decreased C allocation to roots in response to higher nutrient availability. Leggett and Kelting (2006) found that N fertilization of Loblolly pine plantations not only increased aboveground and belowground biomass but also increased soil C pools. However, other estimates suggest that N loading on ecosystems only makes a minor contribution to C sequestration (Nadelhoffer, 1999) or does not likely account for significant C storage (Korner, 2000; De Viries et al., 2006) or in some cases may actually reduce ecosystem productivity and C storage (Aber, 1989; Schulze, 1989). N fertilisation may affect the pool of soil organic C (SOC) through the changes of litterfall. Both negative (Gundersen et al., 1998) and positive (Majdi and Anderson, 2005) effects of N fertilisation on root biomass have been observed. A detailed review by Nadelhoffer (2000) found it is likely that N deposition will decrease fine-root biomass but stimulate fine-root turnover and production. N addition or fertilization to soils has also shown variable effects on SOC storage by enhanced, reduced or unchanged heterotrophic and/or autotrophic respiration as a response to litter C/N ratio (Bowden et al., 2000, 2004; Burton, et  al., 2004; Cleveland and Townsend, 2006; Phillips et al., 2007; Mo et al., 2006, 2008; Hyvönen et al., 2008). With increasing rates of anthropogenic N deposition (Vitousek et al., 1997; Galloway and Cowling, 33  2002), there is a strong need to understand links between ecosystem N inputs and C sequestration. 2.4 Monitoring of C and water cycling in terrestrial ecosystems 2.4.1  Global flux tower network (FLUXNET) The EC technique is commonly used to directly measure the CO2, water vapor and energy exchange between the atmosphere and terrestrial ecosystems (Baldocchi, 2008). A global network of collaborating regional networks, FLUXNET, began in 1997 (Baldocchi, 2008) and today, there exist more than 400 EC-flux towers across the globe. EC measurements are a rich source of information on temporal variability and environmental controls of CO2 exchange between the atmosphere and terrestrial ecosystems (Law et al., 2000). These global EC datasets allow us to (1) explore emergent-scale properties by quantifying how the metabolism of complex ecosystems respond to perturbations in climate variables on diurnal, seasonal, inter-annual and decadal time scales and elucidate physical and biological controlling factors (Baldocchi, 2008; Law et al., 2000); (2) examine carry-over effects that may be introduced by either favorable or deleterious conditions during antecedent years (Barford et al., 2001); (3) observe a disturbance and the recovery from it or to span a natural sequence of ecological development coupled with fluctuations in climate (Amiro et al., 2006; Stoy et al., 2006); and (4) test and validate ecosystem process models (Chen et al., 2007a; Urbanski et al., 2007), since most of these models span timescales from hours to decades. Although the available EC data have been rapidly accumulating, there are some issues associated with its use due to difficulties/uncertainties in (i) assessing/interpreting the associated measuring biases of EC data, and (ii) upscaling of the EC fluxes at the ecosystem (typically less than 1-3 km 2  for each site) to larger scales, e.g. landscape and regional scales.   34  2.4.2  CO2 concentration measurements, data assimilation and CarbonTracker Observations of CO2 over the continent within the atmospheric boundary layer reflect exchange processes occurring at the surface at a regional scale (10 2  – 105 km2). The flux information contained in CO2 concentration data represents footprints of up to 10 5  km 2  (Gloor et al., 2001; Lin et al., 2004; Chen et al., 2008), which is several orders of magnitude larger than the direct EC-flux footprint. These measurements form a record of integrated net CO2 exchange from multiple processes, geographic areas, and times (Peters, 2007). CO2 concentration measurements are made from tall towers (Bakwin et al., 1995, 1998, 2004) or balloons that reach the top of the planetary boundary layer (PBL) of the Earth. These measurements therefore provide significant information to upscale from site to region. Moreover, the number of CO2 mixing ratio measurements above the land surface, made by either tower or aircraft, is steadily increasing. Data are collected by numerous agencies around the world, for instance, the National Oceanic and Atmospheric Administration’s (NOAA’s) Earth System Research Laboratory (ESRL) monitors CO2 in the atmosphere as a contribution to the North American Carbon Program (NACP) (Wofsy and Harriss, 2002). In addition, Peters et al. (2007) suggest that direct satellite observations of CO2 are available already for the upper troposphere, whereas near-surface CO2 from space will become available within several years to augment the current efforts. Previous efforts to interpret the signal of regional CO2 exchange making use of tower concentration data have focused on simple one-dimensional PBL budgets that rely on gradients in CO2 concentrations between the PBL and the free troposphere (Bakwin et al., 2004; Chen et al., 2008; Betts et al., 2004; Helliker et al., 2004). These methods are limited to monthly resolution because of the need to smooth and average over several synoptic events (Matross et al., 2006). The atmosphere integrates surface fluxes over many temporal and spatial scales and links scalar sources and sinks with concentrations and fluxes. This principle has been successfully used to develop inverse models to estimate annual C budgets (Tans et al., 1990; Gurney et al., 35  2002; Stephens et al., 2007). However, due to model limitations and paucity of continental CO2 observations these studies have yielded C fluxes only at coarse resolution, over large spatial regions (i.e. at continental scale, Rodenbeck  et al., 2003). Recently, Stephens (2007) and Yang (2007a) showed that a large set of atmospheric inverse model results were inconsistent with free troposphere CO2 concentration. Transport biases are one of the largest unknown sources of error in flux inversions (Gurney et al., 2002). A powerful way to use all these CO2 data is in a data assimilation system, which combines diverse data and models into a unified description of a physical/biogeochemical system consistent with observations (Marquis and Tans, 2008). Such a new data assimilation system called CarbonTracker, has been built in the NOAA’s ESRL, using a state-of-the-art atmospheric transport model coupled to an ensemble Kalman filter (Peters et al., 2007). The first release of CarbonTracker marks a significant step in our ability to monitor month-by- month surface sources and sinks of CO2 (Peters et al., 2007), with all results, data and code of CarbonTracker freely available from NOAA ERSL’s CarbonTracker web site, http://carbontracker.noaa.gov. 2.4.3 Stable C isotope measurements The information on the biological and physical processes that exchange CO2 between terrestrial ecosystems and the atmosphere is recorded by the signals of 13 C / 12 C ratio in the atmosphere CO2. The stable isotope ratio of CO2 ( 13 C) in the atmosphere contains unique information to study the overall balance of surface CO2 fluxes (Tans. 1980; Tans et al., 1993). As reviewed by Suits et al. (2005), C isotopes can be helpful in investigations of the following four aspects at ecosystem and local scales: (i) plant water-use efficiency and the response of plants to changes in precipitation and relative humidity (Faquhar and Richards, 1984; Farquhar and Hubick, 1998; Condon et al., 1993; Hall et al., 1993; Ekblad and Ho g̈berg, 2001; Bowling et al., 2002), (ii) variation in light distribution and stand structure (Berry et al., 1997; Buchmann et al., 1997 a,b; Le Roux et al., 2001), (iii) recycling of respired CO2 (Keeling, 1961; Schleser and Jayasekera, 1985; Lloyd et al., 1996; Sternberg et al., 1989), and (iv) determining the relative contributions of photosynthesis and respiration to 36  the total net ecosystem exchange (Yakir  et al., 1996; Bowling et al., 2001; Oge é et al., 2003; Lai et al., 2003; Knohl and Buchmann, 2005). Multiple efforts to measure stable C isotopes at both flux towers and flask stations around the world have been achieved. Several new techniques (e.g. automated measurement systems, tunable diode laser (TDL) spectrometer and pulsed quantum cascade laser spectrometer) have also been applied to this investigation (Bowling et al., 2003; Lai et al., 2004; Griffis et al., 2005). Isotope measurements, however, are still lacking considering the land surface diversity/heterogeneity. This shortage of long-term measurements and of sampling frequency still limits isotopic studies and applications to various spatial/temporal scales. Efforts to expand CarbonTracker to assimilate observations of C isotopes ( 13 CO2, 14 CO2) and other observation types (eddy-flux measurements, satellite radiances) are underway. Such observations could facilitate attribution of C fluxes to specific processes such as fossil fuel burning, biomass burning, or agricultural food and biofuel production (Peters et al., 2007). 2.4.4 Satellite monitoring Satellite-borne remote sensing offers unique opportunities to parameterize land surface characteristics over large spatial extents at variable spatial and temporal resolutions. For example, the Moderate Resolution Imaging Spectroradiometer (MODIS) provides a global dataset every 1-2 days with 36 bands. The spatial resolution of MODIS (pixel size at nadir) is 250 m for channels 1 and 2 (0.6µm - 0.9µm), 500 m for channels 3 to 7 (0.4µm - 2.1µm) and 1000 m for channels 8 to 36 (0.4µm - 14.4µm), respectively. Data from the satellite-borne MODIS are currently used in the calculation of global weekly primary productivity (GPP) at 1-km spatial resolution (Running et al., 2004). Other sensors, such as the Landsat Thematic Mapper sensors carried onboard the Landsat series of satellites, acquire images at a 30-m spatial resolution with a 16-day interval. In general, water evapotranspired from ecosystems into the atmosphere will reduce the land surface air temperature (Ta). Reduction in soil moisture will decrease plant transpiration and evaporation from soil and plant surfaces. Reduction in ET will increase Ta. Ta can be derived 37  from remotely-sensed thermal-infrared (TIR) band (8-14 microns) from various operational satellites. Based on the relationship between Ta and ET, remotely sensed Ta has been used to estimate regional ET (Gillies et al., 1997; Kite et al., 2000; Su, 2000). The existing thermal imaging sensors provide adequate coverage of thermal dynamics that are useful for operational monitoring applications of ET. For example, thermal images at 15 minute intervals and at a spatial resolution of 5 kilometers can be obtained from the NOAA Geostationary Operational Environmental Satellites (GOES), and TIR data at a fine spatial resolution (60 m or 120 m) with a much longer time interval (16 days) have been provided by the Thematic Mapper (TM) and the Enhanced TM Plus (ETM+) instruments on Landsat 5 and Landsat 7. ET, the largest component of water loss from ecosystems, plays an important role in affecting soil moisture, vegetation productivity, C cycle, and water budgets in terrestrial ecosystems (Dirmeyer, 1994; Pielke et al., 1997; Betts and Ball, 1997). Verstraeten et al. (2008) provided a comprehensive review of remote sensing methods for assessing ET and soil moisture content across different scales and Kalma et al. (2008) reviewed satellite-based algorithms for estimating ET and land surface temperatures at local, regional and continental scales, with particular emphasis on studies published since the early 1990s. In addition, as Marquis and Tans  (2008) reviewed, satellite-based instruments can also provide information about CO2 concentration in the atmosphere, but no current satellite- borne instrument comes close to providing the accuracy, precision, and continuity required to determine regional CO2 concentrations and local fluxes. Future satellites, including the Greenhouse Gases Observing Satellite (GOSAT) (Basu et al., 2011), are expected to provide more accurate CO2 measurements than do today’s satellites (Marquis and Tans, 2008; Crisp et al., 2004). 2.4.5 Other airborne measurements Besides satellite monitoring, other airborne observation techniques (e.g. aircraft, airplane and land surface remote sensing) have been developed rapidly since the latest decade. For 38  instance, a new approach is LiDAR (Light Detection and Ranging), which is a remote sensing technology that determines distances to an object or surface using laser pulses. LiDAR data have proved to be highly effective for the determination of three dimensional forest attributes. The suitability of airborne LiDAR for the determination of forest stand attributes including leaf area index (LAI) and the probability of canopy gaps within different layers of canopy has been widely acknowledged by various studies (Coops et al., 2004, 2007). The interpreted LiDAR data have been further used for landscape C modeling and scaling (Hilker et al., 2008; Chen et al., 2009c). 2.5 Modeling of C and water dynamics in terrestrial ecosystems based on remote sensing The land surface of the Earth represents significant sources, sinks, and reservoirs of C, heat and moisture to the atmosphere. C and energy fluxes and water cycles at soil-atmosphere and plant-atmosphere interfaces are therefore important land surface processes. Due to the complexity and non-linearity of C, N and water dynamics in terrestrial ecosystems, various modeling tools are needed for better understanding of these biogeochemical and hydrological processes and their feedback mechanisms with the land surface climate system (Rannik et al., 2006). The rapidly proliferating volume of spatial data generated by remote sensing has created a significant challenge in terms of designing model algorithms. A spatially distributed process-based model uses spatial data for computing eco-hydrological and biophysical processes. The model algorithms represent hypotheses that can be assessed and potentially revised after confrontation with remote sensing and land surface-based observations. It is well known that realistic simulations of C and water dynamics in terrestrial ecosystems are of critical importance, not only for the surface microclimate, but also for the large-scale physics of the atmosphere (Cox et al., 1999; Gedney et al., 2006; Dickinson et al., 2002). Depending on the scientific objectives or applications, C and water cycle models have been designed with varying degrees of aggregation with respect to ecosystem processes, components, and remote sensing data as model inputs. Such models can be flagged by land surface, ecosystem and hydrological models based on their objectives and emphases. The 39  former focus on ecosystem processes and the interactions between ecosystems and the atmosphere; while the latter place emphasis on the land surface hydrology processes, including lateral flow resulting from catchment topography. 2.5.1 Land surface and ecosystem modeling Global climate and the global C cycle are controlled by exchanges of water, C, and energy between the terrestrial biosphere and atmosphere. Thus land surface models (LSMs) are essential for the purpose of developing predictive capability for the Earth's climate on all time scales (Sellers et al., 1997). Most current LSMs can be associated with three broad types (Seth et al., 1994): soil-vegetation-atmosphere transfer schemes (SVATS), potential vegetation models (PVMs), and terrestrial biogeochemistry models (TBMs). The first generation of SVATS evolved from simple bucket schemes focusing on soil water availability (Manabe et al., 1969), through the schemes of Deardorff et al. (1978). Marked improvements of the second generation from the first generation are the separation of vegetation from soil and the inclusion of multiple soil layers for dynamic heat and moisture- flow simulations (Chen et al., 2007a). The following models, for example, Biosphere- atmosphere Transfer Scheme (BATS) (Seth et al., 1994), Simple Biosphere model (SiB) (Sellers et al., 1986, 1997) and Canadian LAnd Surface Scheme (CLASS) (Verseghy et al., 1993, 1999), are typical second-generation SVATS models. The second generation SVATS firstly modeled plant physiology in an explicit manner in GCMs (General Circulation Model or Global Climate Model) (Henderson et al., 1993). For most second-generation SVATS, land cover was fixed, with seasonally-varying prescriptions of parameters such as reflectance, LAI or rooting depth (Wang et al., 2002a,b; Kickert et al., 1999; Kley et al., 1999; Schwalm et al., 2001). Some SVATS incorporated satellite data to characterize more realistically the seasonal dynamics in vegetation function (Kickert et al., 1999; Bonan et al., 1994). The latest (third generation)  SVATS used more recent theories relating photosynthesis and plant water  relations to provide a consistent description of energy exchange,  ET, and C exchange by plants (Chen et al., 2007a; Sellers et al., 1996). In our effort in understanding the impact of climate change on terrestrial ecosystems, energy, water, and C cycles need to 40  be modeled simultaneously (Sellers et al., 1996; Williams et al., 2001). Recently, most of SVATS have thus been enhanced to include the CO2 flux between the land surface and the atmosphere, such as Simple Biosphere model version 2.0 (SiB2) (Sellers et al., 1996), Integrated BIosphere Simulator (IBIS) (Foley et al., 1996), The National Center for Atmospheric Research LSM (NCAR-LSM) (Bonan et al., 1995), BATS (Dickinson et al., 2002), CLASS-C (Wang et al., 2002a,b) and Ecosystem-Atmosphere Simulation Scheme (EASS)  (Chen et al., 2007a). The earlier generation of PVMs comprised a suite of schemes that focus on modeling distributions of vegetation as a function of climate (Holdridge et al., 1947; Prentice et al., 1990) without influences of anthropogenic or natural disturbance. The second generation of PVMs included more sophisticated modules to account for factors controlling vegetation distributions, such as competition, varying combinations of plant functional types, and physiological and ecological constraints (Prentice et al., 1992). TBMs developed from scaling up local ecological models, are process-based models that simulate dynamics of energy, water, and C and N exchange among biospheric pools and the atmosphere (Seth et al., 1994). Few of the existing TBMs incorporate PVMs. These models are not applicable to transient climate change experiments without coupling with PVMs. In recent decades, the interactions among soil, vegetation and climate have been studied intensively and modeled successfully on the basis of water and energy transfer in the soil- vegetation-atmosphere system (Seth et al., 1994; Sellers et al., 1986; Verseghy et al., 1999; Verseghy et al., 1993; Zhang et al., 2003). Also the construction and refinement of LSMs have received increasing attention (Sellers et al., 1996; Viterbo et al., 1995; Christopher et al., 2004). Combination of these three different LSMs and utilization of remotely sensed land surface parameters are critical in the future LSM development, because of (1) the tight coupling of exchanges of water, energy and C between the land surface and the atmosphere; (2) the sophisticated impact/feedback mechanisms between climate change and terrestrial ecosystems; and (3) increasingly strong anthropogenic alterations to land cover. On-line coupling of a LSM with a GCM is needed for studying interannual to multi-decadal climate variations. 41  Several model inter-comparisons have focused on evaluating SVATS and TBMs with particular objectives. For instance, the Project for Inter-comparison of Land-surface Parameterization Schemes (PILPS) was initiated to evaluate an array of LSMs existing in GCMs (General Circulation Model or Global Climate Model) (Henderson et al., 1993); while the AMMA (African Monsoon Multidisciplinary Analysis) Land Surface Model Inter- comparison Project (ALMIP) was conducted to get a better understanding of the role of soil moisture in land surface processes in West Africa (de Rosnay et al., 2009). Coordinated land surface modeling activities have improved our understanding of land surface processes (de Rosnay et al., 2009). 2.5.2 Spatially-distributed hydrological processes modeling Hydrology and ecosystem have, for the most part, been studied independently. Most LSMs and ecosystem models make an assumption of “flat Earth” with the absence of lateral redistribution of soil moisture. On the other hand, hydrological models have mostly been concerned with runoff production. Spatially-distributed models are needed, especially for hydrological simulation objective, because of heterogeneity of land surface and non-linearity of hydrological processes. Spatially-distributed hydrological models are not only able to account for spatial variability of hydrological processes, but enable computation of internal fluxes and state variables. Such kinds of models are increasingly applied to simulate spatial variability of forcing variables (e.g. precipitation), physiographic characteristics, detailed processes and internal fluxes within a catchment (Liang et al., 1994; Liang et al., 2004; Beldring et al., 2003; Brath et al., 2004; Christensen et al., 2007; Reed et al., 2004). 2.5.3 Modeling dynamics of stable C isotopic exchange between ecosystem and the atmosphere It is recognized that the atmospheric measurements are still too sparse, relative to its spatial variability, to be used for inferring the surface flux at high spatial resolution (Ciais et al., 1995). The use of the isotope ratio as an additional constraint to identify various C sources and sinks can contribute to a significant reduction in the uncertainty. Though available 42  isotopic datasets are being accumulated quickly (Griffis et al., 2005; Ponton et al., 2006; Lai et al., 2006; Lai et al., 2005) isotope measurements are still lacking considering land surface diversity and heterogeneity. This shortage of long-term measurements and of sampling frequency still limits C isotopic studies. Mechanistic ecosystem models that couple micrometeorological and eco-physiological theories have the potential to shed light on how to extend efforts and applications of stable isotopes of CO2 to global C budgeting, because biophysical models have the capacities of simulating isotope discrimination in response to environmental perturbations and can produce information on its diurnal, seasonal and inter-annual dynamics. Few biophysical models, however, have been developed to assess stable C discrimination between a plant canopy and the atmosphere (Suits et al., 2005; Oge’e et al., 2003; Baldocchi et al., 2003). Most existing biophysical models are based on individual leaf level discrimination equations given by Farquhar et al. (Farquhar et al., 1989; Farquhar et al., 1982) and only focus on the land surface layer (ignoring vertical and horizontal advection effects beyond 50~100 m above the ground (Baldocchi et al., 2003b). However, in nature, the convective boundary layer (CBL) integrates the effects of photosynthesis, respiration, and turbulent transport of CO2 over the landscape (Lloyd et al., 1996; Pataki et al., 2003). The influence of the CBL cannot be ignored when using isotope composition of CO2 to investigate biological processes (Bowling et al., 1999), because the effect of atmospheric stability on turbulent mixing/diffusion has an important impact on scalar fluxes and concentration fields within and above canopies (Baldocchi et al., 1995; Leuning et al., 2000). Few such models considering the CBL effects on isotope fractionation have been developed to date (Lloyd et al., 1996, 2001; Chen et al., 2006a, b; Chen et al., 2007a). 2.5.4 Modeling coupled C and water dynamics --- an eco-hydrological approach C and N dynamics and hydrological processes are closely linked. The stomatal conductance (gs) is the key linkage between C assimilation (photosynthesis) and transpiration. An empirical equation is used in the second-generation LSMs to calculate  gs, which is 43  hypothesized to be controlled by the environmental conditions  (Jarvis et al., 1976). While field and laboratory studies have documented that leaf photosynthesis also affects gs. Therefore, Ball et al. (1987) proposed a semi-empirical stomatal conductance formulation (Ball-Woodrow-Berry model), in which gs is controlled by both photosynthesis and the environmental conditions. Most of third-generation LSMs (Ecological models, e.g. SiB2 (Sellers et al., 1996, 1997); CN-CLASS (Arain et al., 2006); Ecosys (Grant et al., 1999, 2007), and EASS (Chen et al., 2007a) fully couple photosynthesis and transpiration processes by employing the Ball-Woodrow-Berry stamatal conductance formulation. In addition to the coupling of hydrological condition and C assimilation through the linkage of gs, C assimilation is also coupling with N dynamics through another biochemical parameter, 25maxcV --- maximum carboxylation rate at 25 °C. In the photosynthesis model proposed by Farquhar et al. (1980), the net photosynthetic rate Anet at leaf level is a function of two tightly-correlated parameters 25maxcV and 25 maxcJ (the maximum electron transport rate at 25 °C), and is calculated as, djcnet RAAA  ),min(                                                                                                 (2.1) where Ac and Aj are Rubiso-limited and light-limited gross photosynthesis rates, respectively, and Rd is the daytime leaf dark respiration and computed as Rd = 0.015 Vc max. Ac and Aj are expressed as, )/1( * max occc c cc KOKC C VA     (2.2a) and, )2(4 * * max    c c j C C JA  (2.2b) where Cc and Oc are the intercellular CO2 and O2 mole fractions (mol mol −1 ), respectively; * is the CO2 compensation point without dark respiration (mol mol −1 ); Kc and Ko are 44  Michaelis-Menten constants for CO2 and O2 (mol mol −1 ), respectively. In the nutrient-limited stands, Anet is generally limited by Ac, while Ac is dominantly controlled by a parameter maxcV (see Eq. 2.2a). Many research results showed 25maxcV is very sensitive to leaf N status (more specifically leaf Rubisco-N) (Dickinson et al., 2002; Wilson et al., 2000a,b; Wilson et al., 2001; Warren et al., 2001). As a result in some ecosystem models (i.e. C&N-CLASS (Arain et al., 2006)), 25maxcV  is calculated as a nonlinear function of Rubisco-N following observations made by Warren and Adams (Dickinson et al., 2002):  0 25 max 8.1exp(1)( rc NNV    (2.3) where α is the maximum value of 25maxcV  and Nr0 is the leaf Rubisco-N (g N m −2  leaf area) in the top canopy. The coupled C, N and water processes have been carefully considered in most of the third- generation LSMs (e.g. SiB2: Sellers et al., 1986, 1997, 1997; CN-CLASS: Arain et al., 2006; and Ecosys: Grant et al., 1999, 2007), the models’ grids, however, are isolated from their neighboring grids mainly due to the availability of input data. Vertical soil hydrological processes are hard to be realistically simulated if the lateral flows are ignored by assuming that the Earth is “flat”. However, simulations of the topographically-driven lateral water flows are important components in most of spatially-distributed models, while the detailed eco-physiological processes are weakly represented (Govind et al., 2009). Much effort to bridge these two different models has been increasingly made (Rodriguez et al., 2001; D’odorico et al., 2004; Govind et al., 2009; Creed et al., 1998; Band et al., 2001; Porporato et al., 2002; Porporato et al., 2003; Daly et al., 2004; Chen et al., 2005). However, a model coupling approach --a full combination of ecosystem model and hydrological model, i.e. eco- hydrological modeling, is still lacking. 2.5.5 Applications of remotely-sensed data in ecohydrological modeling Remote sensing techniques, which inherently have the ability to provide spatially comprehensive and temporally repeatable information of the land surface, may be the only 45  feasible way to obtaining data needed for land surface and ecological modeling (Sellers et al., 1986; Gurney et al., 2003; Kite et al., 1996; Engman et al., 1996; Melesse et al., 2008). The most common rationale for interfacing remote sensing and land surface-ecosystem models is using remotely sensed data as model inputs (Plummer et al., 2000). These input data, corre- sponding to forcing functions or state variables in ecological modeling, include land cover (LC), LAI, normalized difference vegetation index (NDVI), and the fraction of photosynthetically active radiation (fPAR) (Sellers et al., 1986; Running et al., 1998; Chiesi et al., 2002; Loiselle et al., 2001). Another effort is the direct estimation of GPP and NPP (Goetz et al., 1999; Seaquist et al., 2003) of biomass (Seaquist et al., 2003; Bergen et al., 1999) and of plant growth (Maas et al., 1988; Kurth et al., 1994), by making use of fPAR and NDVI. It has been shown that the direct estimation has lower accuracy than the integration of remotely sensed data with process-based models (Goetz et al., 1999). Remotely sensed data have also been used to parameterize hydrological models (Chen et al., 2005; Kite et al., 1996; Boegh et al., 2004). For instance, a hydrological model (TerrainLab) was further developed using remotely sensed data as inputs (Chen et al., 2005). TerrainLab is a spatially distributed, process-oriented hydrological model using the explicit routing scheme of Wigmosta et al. (2004). This model has been applied to flat areas such as boreal and wet land region (Chen et al., 2005; Govind et al., 2009), but it has not yet been applied to mountainous areas. Different from traditional hydrological models, which have coarse spatial resolutions, the grid-based-distributed eco-hydrological models have a high demand for spatial data (Kite et al., 1996; Montzka et al., 2008). Some researchers highlight that the main obstacles in current distributed eco-hydrological modeling is the lack of sufficient spatially distributed data for input and model validation (Stisen et al., 2008). Remote sensing can potentially fill in some of the gaps in data availability and produce means of spatial calibration and validation of distributed hydrological models. As a result the application of remote-sensing techniques in hydrological studies and water resources management has progressed in the past decades (see review by Kite et al., 1996). 46  In general, the applications of remotely sensed data in eco-hydrological modeling can be in the two ways (Kite et al., 1996; Chen et al., 2005;  Boegh et al., 2004; Montzka et al., 2008; Stisen et al., 2008; Ritchie et al., 1996; Schultz et al., 1996; Melesse et al., 2007; Schmugge et al., 2002; Jain et al., 2004; Pietroniro et al., 2005; French et al., 2006): (i) multispectral remote sensing data are used to quantify surface parameters, such as vegetation types and density. Although the usefulness of remotely sensed data is widely recognized, there remain few cases where remotely sensed data have been actually used in eco-hydrological simulations. Difficulties still exist in choosing the most suitable spectral data for studying hydrological processes as well as in interpreting such data to extract useful in formation (Chen et al., 2005; Kite et al., 1996; Engman et al., 1996); and (ii) processed remotely sensed data are used to provide fields of hydrological parameters for calibration and validation of eco-hydrological models, such as precipitation (Kite et al., 1996; Wang et al., 2001), and soil moisture (Jackson et al., 1993; Hollenbeck et al., 1996., 1996; Kim et al., 2002; Koster et al., 2006). Koster and others (2006) pointed out that remotely sensed data take the form of emitted and reflected radiances and thus are not the type of data traditionally used to run and calibrate models. Hence, it is important to understand and develop relationships between the electromagnetic signals and hydrological parameters of interest (Chen et al., 2005). Kite and Pietroniro (1996) stated that the use of remote sensing in hydrological modeling was limited. Even though a number of new sensors have been launched since then and research has documented that remote sensing data have promising perspectives, operational uses of satellite data in hydrological modeling still appear to be in its infancy (Stisen et al., 2008). 2.6 Research gaps in C and water flux estimates and scaling approaches A variety of methods are being used in the C and water cycles studies. As shown in Figure 2.3, different approaches have different temporal and spatial scales. The most direct measurements of the terrestrial C flux are made either at the plot scale (10 −2 -10 1  m 2 ), e.g. using biometric methods and various forms of chamber, or at the ecosystem (patch) scale (10 4 - 10 6  m 2 ), using the EC technique. Eco-hydrological and ecosystem modeling and 47  remote sensing estimations are generally available across variable spatiotemporal scales. These estimates are normally available within a nested framework  that permits a progressive comparison of measurements made by  surface instrumentation (scale: 1 to 10 m), surface flux equipment  (10 m to 1 km), airborne remote sensing equipment (100 m to several km), satellite remote sensing (30 m to global scale) and EC tower (1-3  km).  Figure 2.3 Temporal and spatial scales of different approaches. The atmosphere integrates surface fluxes over many temporal and spatial scales and links scalar sources and sinks with concentrations and fluxes. This principle has been successfully used to develop inverse models to estimate annual C budgets (Tans et al., 1990; Enting et al., 1995; Fan et al., 1998; Bousquet et al., 1999; Gurney et al., 2002, 2003). However, due to model limitations and paucity of continental CO2 observations these studies have yielded C fluxes only at coarse resolution, over large spatial regions (Gurney et al., 2004, 2005, 2008). Progress in C balance studies has been achieved at both ends of the spatial scale spectrum, either large continents (larger than 10 6  km 2 , e.g. global inverse modeling) or small vegetation stands (less than 1-3 km 2 , e.g. EC-measurements). Methods to estimate CO2 sources and sinks at the intermediate scale (i.e. landscape to regional scales) between continental and 48  local scales are less well advanced. Moreover, the C cycle in different regions can vary markedly in response to changing climate (Friedlingstein et al., 2003). Reliable estimates of terrestrial C sources and sinks at landscape to regional spatial scales (finer than those used in global inversions and larger than local EC flux measurements and roughly defined as the range between 10 2  and 10 6  km 2 ) are required to quantitatively account for the large spatial variability in sources and sinks in the near-field of a measurement location (Gerbig et al., 2003), as well as fundamental to improving our understanding of C cycle (Crevoisier et al., 2006). It is generally considered unreliable to upscale stand-level fluxes (i.e. EC measurements) to a region by simple spatial extrapolation and interpolation because of the heterogeneity of the land surface and the nonlinearity inherent in ecophysiological processes (Levy et al., 1999). It is also challenging to apply atmospheric inversion technique to regional scales for quantifying annual C budgets because at such intermediate scales the atmosphere is often poorly constrained (Matross et al., 2006; Gloor et al., 1999). Moreover, aggregation errors and errors in atmospheric transport, both within the PBL and between the PBL and free troposphere, can also be obstacles to using these approaches to obtain quantitative estimates of regional C fluxes (Lin et al., 2004). Hence, there is a strong motivation to develop methods to quantify and validate estimates of the C balance at these intermediate scales (Lin et al., 2004; Chen et al., 2008; Bakwin et al., 2004; Matross et al., 2006). Observations of CO2 over the continent within the PBL reflect exchange processes occurring at the surface at a regional scale (10 2  – 105 km2). The flux information contained in CO2 concentration data represents footprints of up to 10 5  km 2  (Gloor et al., 2001; Lin et al., 2004), which are several orders of magnitude larger than the direct EC-flux footprint. This information is therefore needed in our effort to upscale from site to region. Moreover, the number of CO2 mixing ratio measurements above the land surface, made by either tower or aircraft, is steadily increasing. Previous efforts to interpret the signal of regional CO2 exchange making use of tower concentration data have focused on simple one-dimensional PBL budgets that rely on gradients in CO2 concentrations between the PBL and the free troposphere (Bakwin et al., 49  2004; Helliker et al., 2004). These methods are limited to monthly resolution because of the need to smooth and average over several synoptic events (Matross et al., 2006). 2.7 Future research directions A synthetic research framework is needed to strengthen the less well researched areas as reviewed in Section 2.6: bottom-up and top-down approaches integrating scalable (footprint and ecosystem) models and a spatially nested hierarchy of observations which include multispectral remote sensing, inventories, existing regional clusters of EC flux towers and CO2 mixing ratio towers and chambers. The current research trends and the future directions in this field include: (i) A synthesis aggregation method --- integrating eco-hydrological and isotopic models, remote sensing and component flux data, is becoming a pragmatic approach towards a better understanding of the coupled C, N and water dynamics at landscape/watershed scales; and (ii) The landscape and regional scale C fluxes are being estimated using an integrated approach involving direct land surface measurements, remote sensing measurements, and ecosystem-, footprint- and inversion- modeling. 2.7.1 Development of a spatially explicit ecohydrological modeling framework Coupled modeling will help refine the experimental and instrumental design and generate cross-disciplinary hypotheses that can be tested in the experiment. Eco-hydrological models are powerful tools for quantitative and predictive understanding the coupled C, N and water mechanisms. Spatially-explicit eco-hydrological modeling can be used to infer aspects of the land surface system that are difficult to measure by mass and energy balance, and will be critical to improving the accuracy of forecasts of landscape change and C dynamics in the real world. While developing a process-based, coupled-system model is a significant task, the coupling to existing models may provide a relatively easy way forward. For example, a spatially explicit, process-based eco-hydrological model, EASS-TerrainLab, is developed to 50  improve the representation of the coupled C, N and hydrological processes by integrating of two existing models (an ecosystem LSM model---EASS and a distributed hydrological model--TerrainLab). 2.7.1.1 Reviewing of the existing EASS model (Chen et al., 2007a)  Figure 2.4 Structure of the EASS model. Three components (soil, vegetation and the atmosphere) are considered in EASS, which are integrated with two interfaces. The right panel illustrated energy fluxes between these three components. LE, H, Rs, Rl and G are the latent heat flux, sensible heat flux, shortwave radiation, longwave radiation, and soil conductive heat flux, respectively; the subscripts g and c present the energy fluxes at soil- canopy and canopy-atmosphere interfaces, respectively. The left panel describes soil water fluxes. The symbol F represents conductive water flux between soil layers, and F0 represents the incoming water flux from the surface to the top soil layer (i.e. the actual infiltration rate ), and Fb is the water exchange (drainage or capillary rise) between the bottom soil layer and the underground water (Chen et al., 2007a). I 51   EASS is based on a single layer vegetation canopy overlying a seven-layer soil, and includes physically-based treatment of energy and moisture fluxes from the vegetation canopy and through it. It also incorporates explicit thermal separation of the vegetation from the underlying ground (Dickinson et al., 2002). Moreover, EASS includes a scheme with stratification of sunlit and shaded leaves to avoid shortcomings of the “big leaf” assumption (Melesse et al., 2008; De Pury et al., 1997). It has been referred as a “two-leaf” canopy model (Norman, 1980; Chen et al., 2007a). The structure of EASS is shown in Figure 2.4. With spatially explicit input data on vegetation, meteorology and soil, EASS can be run pixel by pixel over a defined domain, such as Canada’s landmass, or any of its parts, or the globe. EASS has flexible spatial and temporal resolutions, as long as the input data of each pixel are defined. In short, EASS has the following characteristics (Chen et al., 2007a): (i) satellite data are used to describe the spatial and temporal information on vegetation; (ii) energy and water exchanges and C assimilation in the soil-vegetation-atmosphere system are fully coupled and are simulated simultaneously. The Ball-Woodrow-Berry stamatal conductance formulation (Ball and Berry, 1987) was employed (see Section 2.5.4); (iii) the energy and C assimilation fluxes are calculated with stratification of sunlit and shaded leaves to avoid shortcomings of the “big-leaf” assumption; and (iv) the lateral movement of water is ignored. 2.7.1.2 Reviewing of the TerrainLab model (Chen et al., 2005) TerrainLab (Chen et al., 2005) was adopted concepts from Wigmosta et al. (1994) and it has been tested in the Southern Study Area (SSA) of the BOREAS region to study the spatio- temporal variation of ET and soil moisture (Chen et al., 2007a). The horizontal boundary of the simulated area is a watershed delineated from a digital elevation model (DEM) where divides between neighboring watersheds are identified. Vertically, the simulation extends from the saturated zone in the soil to the top of vegetation canopy. Within a watershed, the forest ecosystem is divided into basic spatial units, or pixels in remote sensing. Each pixel is 52  treated as a unique vegetation-soil system, except for the ground and runoff water exchanges. Basic model simulations of the physical and biological processes are made at the pixel scale. According to the need of simulating hydrological processes, a pixel is vertically divided into five strata, i.e. overstory, understory, litter or moss layer, soil unsaturated zone, and soil saturated zone. Precipitation, solar radiation, topographic parameters, LC, LAI, and soil properties are the major inputs to the model. The soil water balance component of TerrainLab allows for the spatially explicit simulation of topographically driven lateral subsurface flow and its influence on water table depth and soil water content dynamics based on a raster grid DEM (Chen et al., 2007a). TerrainLab uses the explicit routing scheme of Wigmosta et al. (1994). In this approach the soil profile is subdivided into an unsaturated and a saturated zone. Saturated hydraulic conductivity is assumed to be depth-dependent, while other hydraulic parameters are considered to be vertically homogenous. However, the assumption of vertical homogeneity of permanent wilting point, field capacity, and porosity represents a major simplification in peat lands. The lateral movement of water between each raster grid cell (pixel) of the modeling domain and its maximum eight neighboring pixels based on a multi-predilection-flow algorithm occurs in the saturated zone, i.e. as groundwater flow. Groundwater follows the local hydraulic gradient (3 x 3 pixels) that is assumed to be approximated by local ground surface slopes (Wigmosta et al., 1994; Chen et al., 2005) 2.7.1.3 Development of an ecohydrological model by coupling of EASS with TerrainLab The coupled model (EASS-TerrainLab) explicitly describes ecohydrological processes (coupled C, N and water dynamics). Each process will be explicitly calculated separately for overstory and understory with sunlit and shaded leaf groups. Vertical hydrological processing will be adequately described using detailed water balance equations, the topographically driven lateral water flows will be considered for both saturated and unsaturated soil layers, and their influence on soil water content and water table depth are calculated based on a  53  raster grid DEM (Chen et al., 2005) using the explicit routing scheme of Wigmosta et al. (1994). 2.7.1.4 Model calibration and validation for EASS-TerrainLab Model calibration and validation will conducted in three ways: (i) Modeled surface runoff will be compared with the watershed discharge measurements; (ii) Simulated soil water storage and its changes will be verified with several soil water measurement profiles; and (iii) Simulated ET will be validated with EC-flux tower measured ET based on footprint averaging and the remotely sensed spatially-distributed ET. 2.7.1.5 Model sensitivity analysis and runs under different scenarios A set of sensitivity analysis will conducted to test the coupled capacity. The model will be run under several different scenarios to explore: (i) the effects of surface and sub-surface base flow on GPP and ET; (ii) the influence of mesoscale topography on hydrological processes (runoff and soil storage) and C exchange (water use efficiency and light use efficiency); (iii) the influence of photosynthesis on ET by comparing simulation results using different stomatal conductance schemes; and (iv) the effects of representations of spatial variability on simulations of photosynthesis and ET. 2.7.2 Landscape and regional C and water fluxes estimation: an up- scaling framework Figure 2.5 schematically shows a synthetic up-scaling framework for estimating landscape and regional C and water budgets with a reasonable accuracy. Using a data-model fusion in conjunction of footprint analysis to optimize the eco-hydrological model’s parameters is strongly suggested because the land surface is normally heterogeneous. This is a key component in the up-scaling framework. Model parameters are optimized by minimizing the ‘cost’ function (Tarantola, 2005): 54          b1bTb1oT xxPxxY(x)OCY(x)O 2 1 J(x)    (2.4) where x is the vector of unknown parameters and xb is a vector of a priori values of x; O is the vector of observations and Y is the nonlinear model results, Co is the covariance matrix of observations and Pb is the covariance matrix of a priori parameters. Different from general data-fusion method, Y is weighting averaged by footprint fj for each pixel (j) within the footprint domain as Y = ∑Yjfj. We also note that the dynamics of the PBL is tightly coupled with the land surface processes (e.g. energy / water / C fluxes). But for simplification, we proposed off-line modeling instead of the on-line coupling with regional climate model or the GCMs. The meteorological driving inputs for the eco-hydrological model, therefore, are required for model runs. These inputs can either acquire from spatially extrapolation from climate station (including EC flux towers) measurements, or from the re-analysis data (e.g. the data provided by the National Centers for Environmental Prediction (NCEP)), or from the regional climate modeling outputs. 2.8 Summary After comprehensive reviewing of a variety of approaches being used in research on the C/water cycles, the concluding remarks are as following: Research gaps in this field are (i) The coupled terrestrial C, N and hydrological dynamics are far from well understood, especially at landscape (watershed) and regional scales; (2) Much progress has been achieved at the extreme ends of the spatial-scale spectrum, either large regions/continents or small vegetation stands. Because of the heterogeneity of the land surface and the nonlinearity inherent in eco-physiological and eco-hydrological processes in response to their driving forces, it is difficult to upscale stand level results to regions and the globe by extrapolation. Budgets of C and water at landscape intermediate regional scales (10 2–105 km2) have large uncertainties. 55  A coupled spatially-explicit eco-hydrological model is a powerful tool for quantitative and predictive understanding of the coupled C, N and water mechanism. This modeling framework can be used to infer aspects of the land surface system that are difficult to measure, and will be critical to improving the accuracy of forecasts of landscape change and C dynamics in the real world.  Figure 2.5 An up-scaling framework synthetically integrating eco-hydrological and footprint modeling, remote sensing and land surface measurements. Combining and mutually constraining the bottom-up and top-down methods to reduce their uncertainties using data assimilation techniques is a practical and effective means to derive regional C and water fluxes with reasonably high accuracy. In the proposed up-scaling framework, spatially nested hierarchy of observations, including multispectral remote sensing, inventories, existing regional clusters of EC flux towers and CO2 concentration towers and chambers, are able to be integrated using scalable (footprint and ecosystem and eco-hydrological) models and data-model fusion techniques. 56  Chapter  3: Seasonal controls on interannual variability in C dioxide exchange of a near end-of rotation Douglas-fir stand in the Pacific Northwest, 1997–20063 3.1 Synopsis This study analysed 9 years of EC data obtained in a Pacific Northwest Douglas-fir (Pseudotsuga menzesii) forest (58-year old in 2007) on the east coast of Vancouver Island, Canada and characterizes the seasonal and interannual variability in net ecosystem productivity (NEP), gross primary productivity (GPP) and ecosystem respiration (Re) and primary climatic controls on these fluxes. The annual values (mean ± standard deviation) of NEP, GPP and Re were 357 ± 51, 2,124 ± 125, and 1,767 ± 146 g C m -2  yr -1 , with ranges of 267 - 410, 1,592 - 2,338, and 1,642 - 2,071g C m -2  yr -1 , respectively. Spring to early summer (March - June) accounted for more than 80% of annual NEP while late spring to early autumn (May - August) was mainly responsible for its inter-annual variability (~80%). The major drivers of inter-annual variability in annual C fluxes were annual and spring mean air temperatures and water deficiency during late summer and autumn (July - October) when growth in this Douglas-fir forest was often water-limited. Photosynthetically active radiation (Q) and the combination of Q and soil water content (θ) explained 85% and 91% of the variance of monthly GPP, respectively; and 91% and 96% of the variance of monthly Re was explained by air temperature (Ta) and the combination of Ta and θ, respectively. Annual net C sequestration was high during optimally warm and normal precipitation years, but low in unusually warm or severely dry years.  Excluding results for1998 and 1999, the two years strongly affected by an El Niño/La Niña cycle, annual NEP significantly decreased with increasing annual mean air temperature. Annual NEP will likely decrease whereas both  3  Reprinted from Global Change Biology, Vol. 15(8), 8624-8657, [B. Chen], T.A. Black, N.C. Coops, P. Krishna, R. Jassal, C. Brȕmmer, Z. Nesic (2009), Seasonal controls on interannual variability in carbon dioxide exchange of a Pacific Northwest Douglas-fir forest, 1997 – 2006. Global Change Biology, 15(8), 1962-1981, (DOI: 10.1111/j.1365-2486.2008.01832.x) Copyright (2009), with 57  annual GPP and Re will likely increase if the future climate at the site follows a trend similar to that of the past 40 years. 3.2 Introduction One of the crucial issues in the prognosis of future climate change is the global budget of atmospheric CO2. Terrestrial ecosystems mediate a large part of CO2 flux between the Earth’s surface and the atmosphere, with approximately 120 Pg C as CO2 taken up by photosynthesis and roughly the same amount released back to the atmosphere by respiration annually (Prentice et al., 2001). A detailed understanding of the interactive relationships in atmosphere–biosphere exchange is relevant to ecosystem-scale analysis and is needed to improve our knowledge of the global C cycle (Falk et al., 2008). The metabolism of terrestrial ecosystems is complex and highly dynamic because ecosystems consist of coupled, non-linear processes that possess many positive and negative feedbacks (Levin, 2002; Ma et al. 2007). Complex features of ecosystem metabolism at longer, multi-year time scales are relatively unknown. Imbalances between gross ecosystem photosynthesis or gross primary productivity (GPP) and ecosystem respiration (Re) lead to land surfaces being either CO2 sinks or sources. The magnitudes of sinks and sources have fluctuated on annual and longer time scales due to variable climate, land use change, disturbance, and changes in the age distribution and species composition of ecosystems (Law et al., 2002; Morgenstern et al., 2004; Urbanski et al., 2007). How C budget of major ecosystems will respond to changes in climate is not quantitatively well understood (Black et al., 2000; Baldocchi et al., 2001a; Baldocchi and Wilson, 2001; Law et al., 2002; Barr et al., 2004, 2007). Continuous and long-term measurements of net ecosystem exchange of CO2 (NEE) between ecosystem and the atmosphere provide investigators opportunities and information to (1) explore emergent-scale properties by quantifying how the metabolism of complex ecosystems respond to perturbations in climate variables on diurnal, seasonal, inter-annual and decadal time scales and elucidate physical and biological controlling factors (Black et al., 58  2000; Law et al., 2002; Barr et al., 2004); (2) examine carry-over effects that may be introduced by either favorable or deleterious conditions during antecedent years (Barford et al., 2001); (3) observe a disturbance and the recovery from it or to span a natural sequence of ecological development coupled with fluctuations in climate (Amiro et al. 2006; Stoy et al., 2006); and (4) test and validate ecosystem process models (Grant et al.,1999; Amthor et al., 2001; Law et al., 2000; Chen et al. 2007a,b; Chen & Chen 2007; Urbanski et al., 2007), since most of these models span timescales from hours to decades. To date, only a few decadal-scale CO2 flux studies have been published (e.g. a boreal deciduous aspen site: Barr et al. 2007, a boreal conifer forest site: Dunn et al., 2007; a temperate mixed deciduous forest site at Harvard forest: Urbanski et al., 2007). In the Pacific Northwest, Falk et al. (2005, 2008) reported studies on the temporal dynamics of soil respiration and NEE for an old-growth Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco)– western hemlock (Tsuga heterophylla (Raf.) Sarg.) forest based on 6-years of eddy- covariance (EC) measurements, but no decadal-long CO2 flux studies have been reported. Pacific Northwest coastal forests of United States and Canada cover approximately 10 5  km 2 between Oregon and Alaska and likely play a significant role in the global C cycle (Paw U et al., 2004; Falk et al., 2008). The climate of the Pacific Northwest is characterized by relatively mild winters with high precipitation and moderately warm, dry summers. The low summer rainfall distinguishes this coastal climate from most other temperate regions in Europe, eastern Asia, and the east coast of the United States, where precipitation is relatively constant throughout the year (Waring and Franklin, 1979).  To date, a number of studies have provided estimates of net ecosystem production (NEP, i.e. NEE) and the components of the C balance (i.e. GPP and Re) across different stand ages in this region (Price & Black, 1990, 1991; Drewitt et al., 2002; Paw U et al., 2004; Chen et al., 2004; Harmon et al., 2004; Unsworth et al., 2004; Morgenstern et al., 2004; Humphreys et al., 2006; Falk et al. 2005, 2008). During the 20th century, globally averaged surface air temperature increased by 0.6 ºC (Folland et al., 2001). The increase in temperature (0.91 ºC) over the last 100 years in the Pacific Northwest exceeded that in global temperature with the largest warming rates 59  occurring in winter and spring (Mote, 2003). Little is known about the interaction between climate change and ecosystems in the Pacific Northwest region. There is a need to investigate what are the climate change effects on and the feedbacks from the coastal forests of the Pacific Northwest. In this Chapter, we analyze 9 years of EC data (1997-2006) obtained in a Douglas-fir forest established in 1949 on the east coast of Vancouver Island, Canada. We address three major questions: (i) What are the patterns and variability of NEP over both seasonal and multiyear time scales, and how do climatic factors control the variability in NEP, Re and GPP? (ii) How does low summer rainfall affect NEP, Re and GPP?  (iii) What is the potential C sequestration response of this Douglas-fir forest to future climate change? 3.3 Methods 3.3.1 Site description The study site (DF49) is located at 4952 7́.8"N, 12520 6́.3"W (flux tower location) on the east of Vancouver Island, British Columbia, Canada. The site is 10 km south west of Campbell River (CR) and 9 km west of Georgia Strait. The site is on a glaciated terrain with a 5-10º northeast-facing slope. The elevation of a rectangular area (2.5 × 2.5 km, centred at the flux tower), which covers most of the flux footprint area (Chen et al. 2009c), varies from about 500 m in the southwest to 100 m in the northeast, with the flux tower at 330 m. The site was harvested and slash-burned in 1943, followed by planting of Douglas-fir seedlings in 1949 resulting in a relatively homogeneous stand. The stand covers an area of 130 ha with slightly younger trees extending for at least a km in all directions. It consisted of 80% Douglas-fir (Pseudotsuga menziesii var menziesii (Mirb.) Franco), 17% western redcedar (Thuja plicata Donn ex D. Don) and 3% western hemlock (Tsuga heterophylla (Raf.) Sarg.). In 2007, the height of dominant trees was about 33.5 m; the mean stand density and tree diameter at breast height for live trees >1.3 m tall were about 1100 stems ha -1  and 40 cm, respectively. The overstorey leaf area index (LAI) was estimated to be 7.3 m 2 m −2  (Chen J 60  et al. 2006). The understory is sparse (LAI was ~1 m 2  m 2 ), mainly consisting of salal (Gaultheria shallon Pursh.), Oregon grape (Berberis nervosa Pursh), vanilla-leaf deer foot (Achlys triphylla (Smith) DC), and various ferns and mosses. The soil is a humo-ferric podzol (Quimper sandy loam) underlain by a compacted till at a depth of 1 m (Jungen, 1986). The thickness of the surface organic layer ranges from 1 to 10 cm with a mean of 3 cm.  Below the organic layer, the soil texture gradually changes from a gravelly sandy loam in the upper 40 cm to gravelly sand with increasing depth. The total soil C content to the 1-m depth is 11.5 kg C m −2 , of which 2.5 kg C m −2 is located in the upper 10 cm of soil. More detailed descriptions of the site characteristics can be found in Morgenstern et al. (2004) and Humphreys et al. (2006). 3.3.2 Instrumentation EC instrumentation and measurements at this site were described in detail in Jassal et al. (2007), Humphreys et al. (2006) and Morgenstern et al. (2004). Briefly, EC instruments were mounted on a 45-m-tall, 51 cm triangular open-lattice type tower at a height of 43 m. The EC instruments included a three-dimensional sonic anemometer-thermometer (SAT) (model 1012R2A up to March 2001, and then a model 1012R2 till May 2005, and  thereafter a model R3-50, Gill Instruments, Lymington, UK) and a closed-path infrared gas (CO2/H2O) analyzer (IRGA) (model LI-6262, LI-COR Inc., Lincoln, NE, USA). Different SAT models were overlapped for a certain periods (~ a month) for calibration and comparison when one model was changed to another. Calculated values of EC-measured CO2 flux (see below) using these three different SAT models were identical. The IRGA was maintained in a temperature- controlled box mounted 2 m lower than the SAT. It was operated in absolute mode with N gas scrubbed to remove water vapour and CO2 flowing through the reference cell at 60-80 cm 3  min −1 . Air was drawn from 30 cm below the center of the SAT array into the IRGA sample cell through a 4.2 m-long heated sampling tube (Dekoron type 1300, 4.0mm inner diameter (i.d.), Cable USA, Cape Coral, FL, USA) and a filter located in the temperature controlled box (model ACRO 50, Gelman Sciences, Ann Arbor, MI, USA). A flow rate of 9.5 l min −1  was maintained using a linear AC pump (model SPP-15EBS-101 Gast 61  Manufacturing Inc., Benton Harbor, MI, USA) located at the base of the tower. Delay times for the CO2 and water vapour mixing ratios calculated using the SAT temperature were about 0.8 and 1 s, respectively. A solenoid-valve system was used to automatically zero and calibrate the gas-analysis system once a day at midnight. This was done by releasing, first, N and then CO2 in dry air calibration gas through a tee at the entrance of the sampling tube, causing almost no pressure change in the sample cell (Chen et al., 1999). The calibrations proved to be stable under field conditions (typically less than 0.5% variation in the span and less than 1 μmol μmol−1 variation in the zero-offset measured in the field over the course of a year). The calibration gas cylinders were calibrated in the laboratory using standards provided by the Canadian Greenhouse Gases Measurement Laboratory, Meteorological Service of Canada, Downsview, ON, Canada. The analog IRGA signals required for EC calculations were filtered using a low-pass R-C analog filter (cut-off frequency 50 Hz) to minimize aliasing and then sampled by a data- acquisition system (DAS) (model DaqBook/200, IOtech Inc., Cleveland, OH, USA), equipped with an analog-to-digital converter (16 channels, 16 bit), at a frequency of 125 Hz. To match the non-adjustable sampling speed of the digital SAT output, the analog signals were subsequently down-sampled to 20.83 Hz after applying a digital low-pass filter (cut-off frequency 10 Hz) to the signals. The digital output from the DAS and the SAT were transferred to a PC via the parallel and serial ports, respectively, and written into output files every 30 min. On-line flux calculations were performed by the PC at the site using software written in MATLAB (The MathWorks Inc., Natick, MA, USA). Half-hourly flux of CO2, water vapor and sensible heat fluxes at this height have been continuously measured since September 1997. Change in the storage of CO2 in the air column below the EC measurement level was measured using a four-level profile sampling system (z = 2, 12, 27, 43 m) equipped with variable lengths of 4 mm i.d. Air was drawn sequentially down sample tubes (Dekoron type 1300, 20 mm i.d.) at a flow rate of 22 l min −1 for 7.5 min, providing a single measurement for 62  each sampling height every half-hour. Dekoron tube gives the same pressure drop for all sampling tubes. A sample was diverted into the LI-800 at a flow rate of 0.8 l min −1 . After flushing the IRGA for 60 s an average CO2 mole fraction was calculated from measurements in the following 390 s. The IRGA was calibrated daily at midnight by flushing N and calibration gas through the 2 m height sampling tube. A suite of climate measurements were also made at the site (Morgenstern et al. 2004; Humphreys et al. 2003, 2006). Initially, net radiation (Rn) was measured by a model S-1 net radiometer (Swissteco Instruments, Oberriet, Switzerland) and, since August 1999, by a four- way radiation sensor (model CNR 1, Kipp and Zonen B.V., Delft, The Netherlands) consist- ing of two pyranometers (model CM3) and two pyrgeometers (model CG3) measuring upwelling and downwelling solar and longwave radiation, respectively,  at the 41-m height. Vertical profile of air temperature (Ta) was measured at 8 levels (z = 0.2, 2, 9, 15, 21, 27, 33, 44 m) with 75 μm chromel-constantan thermocouples. Relative humidity was measured at 3 levels (z = 4, 27, 40 m) using shielded humidity probes (model HMP45CF, Vaisala Oyj, Helsinki, Finland). Wind speed and direction were measured using the SAT. Precipitation (P) was measured using two tipping-bucket rain gauges (model 2501, Sierra Misco, Berkeley, CA, USA) mounted on the flux tower at the 25-m height and a precipitation gauge (model I-200B, Geonor A/S Oslo, Norway) (for determining water equivalent of snowfall in winter) installed in a young plantation about 3 km from the site. The vertical profile of soil volumetric water content (depth: 1-2, 10-12, 35-48, 70-100 cm) was measured at two locations using water content reflectometers (model CS 615, Campbell Scientific Inc. (CSI), Logan, UT, USA). Because the root density in Douglas-fir stands tends to decrease exponentially with depth, thus extracting water non-uniformly with depth (Nnyamah et al., 1977), we followed the method developed by Baldocchi et al. (2004) to compute a representative metric of soil water content by weighting soil water content by an exponential distribution of root density with depth. Measurements of the vertical profile of soil temperature (Ts, depth: 2, 5, 10, 20, 50 cm) were made with copper-constantan thermo- couples, which were located about 5 m from one of the CS-615 profiles. Soil respiration 63  measurements were made using automated closed-type CO2 efflux chambers started in June, 2003 (Drewitt et al., 2002; Jassal et al. 2007). 3.3.3 Data processing and analysis 3.3.3.1 Data processing and quality controls The half-hourly NEE was computed as the sum of EC-measured CO2 flux (Fc, positive upward) and the rate of change in CO2 storage (Sc) in the air column from the ground to the EC measurement height, i.e. cc SFNEE  . Fc was calculated as ______ ''cwF ac  , where _ a  is the mean molar density of dry air, ____ ''cw  is the covariance of vertical velocity ( w , in m s-1) and CO2 molar mixing ratio ( c , in μmol mol -1 ), where the over bar and prime indicate mean and fluctuation, respectively. A coordinate rotation was applied to bring the mean vertical and transverse wind speed components to zero for each half hour. Time lags, due to sample travel time and adsorption in the sample line, were determined by maximizing the covariance between w and c . Sc was computed as dzzc t a    _____ )( based on the profile measurement of CO2 molar mixing ratio ( c  at z = 2, 12, 27, 43 m) using the c  measured in the half hours before and after the current half hour. Net ecosystem productivity (NEP) was calculated as NEP = -NEE. Measurements made during rain or snow fall periods were not explicitly excluded in the processing. However, half-hourly data were filtered to exclude high rates of errors indicated by the sonic and IRGA error flags typically attributable to heavy rainfall. After instrument calibration, quality control procedures were implemented. Fluxes with statistics that did not fall within reasonable limits and/or occurred during instrumental function were removed. We only accepted the half-hourly NEE measurements under well-mixed conditions indicated by the friction velocity ( u ) being larger than 0.3 m s -1  during nighttime (see Morgenstern et al. 2004). The overall half-hourly NEE coverage ratios after applying data quality criteria and removing the data when u < 0.3 m s -1  during night time for the measurement period (15 64  October 1997 to 31 December 2006) were 97.3%, 29.3% and 54.4% for daytime, night time and 24-hours, respectively. The average coverage ratio of our NEE measurements is slightly higher than typical values (~50%) for continuous EC-measurement towers (Falge et al., 2001). To evaluate long-term (e.g. seasonal, annual and multi-year) EC CO2 flux data and relate fluxes to flux chamber measurements, an analytical footprint model based on Eulerian advection-diffusion equations was used for multi-year flux footprint analysis (Chen et al., 2009b). The results showed that the cumulative annual footprint area accounting for 80% of the EC flux was a relatively small area around the tower (< 1 km 2 ), within which the canopy was relatively homogeneous. 3.3.3.2 Gap-filling, time integration and uncertainties The missing half-hourly meteorological and Ts and soil water content (θ) data were filled using (1) appropriately adjacent measurements from adjust levels in the measurement profile (e.g. air temperature profile), and (2) the mean diurnal method with a 15-day window. A 15- day window was moved in 5-day increments. When necessary, the 15-day window was expanded to ensure a minimum of 3 valid observations for each half-hour of the day for the averaging period. The Fluxnet-Canada Research Network procedure for NEE gap-filling and partitioning of NEE into Re and GPP (Barr et al., 2004, 2007) was followed. Briefly, short gaps (< 2h) were filled using the linear interpolation method. NEP was calculated as –NEE. Re was estimated as NEE during periods when GPP was known to be zero, i.e. at night and during the cold season (i.e. Ta < -1ºC). Missed Re at night and daytime Re (Red) were estimated using an empirical ( )( se TfR  ) model, fit to the measured data. An exponential function (i.e. a Q10 model) was used in this study (Black et al. 1996; Goulden et al., 1996; Morgenstern et al., 2004):   10/)(10),( refs TTrefwse QRrtTfR   ,                                                                                (3.1) 65  where Ts is the soil temperature at 5-cm depth; refR is the respiration rate at a reference temperature refT (here 10 o C was used) and 10Q is the factor by which Re increases for a 10  o C increase in temperature (Lloyd and Taylor, 1994); and rw is a time-varying parameter that quantifiers the fractional departure of Re from the mean ‘wet’ Re vs. Ts relationship (the term inside the square bracket in Eq. (3.1), fit to high soil water data only). The value for rw in Eq. (3.1) is independent of Ts at 5-cm depth and thus quantifies the effect of other time-varying factors on Re, including (potentially) Ta, Ts at other depths, soil moisture, LAI, root and shoot growth, soil microbial dynamics and the size of the labile C pool. We fit the parameters in two steps. First, the constants refR and 10Q  were fit to all measured (night-time) Re and Ts data from periods with high soil water (soil water content of the 0-60 cm soil layer, θ0-60cm > 0.2 m 3  m -3 ), with rw fixed at 1.0. Second, the temporal variation in parameter rw was fit using a flexible, moving window of 100 measured (not-missing) data points, corresponding to a period of 3–14 days. The window was moved in increments of 20 points. For each window, rw was evaluated by linear regression, forced through the origin, of measured Re vs. the term in the square bracket of Eq. (3.1). The use of Eq. (3.1) has two advantages over other empirical Re models, such as annual exponential or logistic relationships: it is responsive to short-term and seasonal variations in Re that are not driven by temperature and it gives conservative estimates of Re at high temperatures that may lie outside the range over which the model is fit. GPP then was estimated as NEP + Red (warm season periods) or zero (cold season periods). Gaps in GPP were filled using an empirical light response curve, fit to the measured data. The rectangular hyperbolic model (Hollinger et al., 1999; Lee et al., 1999; Griffis et al., 2003) was used, max max ( , ) w A Q GPP f Q t p Q A        ,                                                                                        (3.2) where Q is photosynthetically active radiation, is the quantum yield, maxA  is the light- saturated asymptote for canopy photosynthetic capacity, and wp  is a time-varying parameter that quantifiers the fractional departure of GPP from the mean  GPP vs. Q relationship. We 66  evaluated the seasonal variations in wp , and maxA  by fitting Eq. (3.2) to measured, daytime NEP and Q data using a flexible, moving window of 240 measured (not-missing) data points, moved in increments of 48 points. The width of the moving window was typically 7-14 days but increased during periods with significant data gaps. Gaps in NEP during daytime were filled as the difference between modeled Re using Eq. (3.1) and modeled GPP using Eq. (3.2). Although the methodology used to fill gaps in NEP and partitioning NEP into Re and GPP introduced a small degree of uncertainty in the absolute magnitudes of Re and GPP, sensitivity testing showed that the effects on annual NEP, Re and GPP were relatively small. We also compared our gap-filled dataset with one produced using the neural network method of Papale and Valentini (2003) and found that the two datasets agreed well. To assign statistical uncertainties to the annual NEP, Re and GPP, we first estimate an error distribution for each half-hour, consisting of (i) instrumentation error calculated via daily difference method of Hollinger and Richardson (2005) and Richardson et al. (2006), (ii) error associated with the model fit for Re and GPP, and (iii) error associated with filling Ts and Q during gaps.  These error distributions were then used to generate different 500 simulations of the NEP time series, from which confidence intervals were bootstrapped using the Monte Carlo method (Efron and Tibshirani, 1993). The range of annual NEP estimation was determined using one-sampling percentile method and at 95% confidence intervals (C.I.). Random sampling of NEP error populations was used to determine the uncertainties in the annual sums of NEP. Each year was divided into 24 periods (12 seasons, day and night). The uncertainties for gap- filled half-hourly NEP were estimated by randomly drawing from the error population defined by half-hours with valid NEP observations (NEPobserved  - NEPestimated). The mean difference between observations and estimates was 0 and such random error in the estimates of annual NEP was found to be within 30 g C m-2 (see also Morgenstern et al., 2004; Schwalm et al., 2007). The biases for valid half-hourly NEP observations were also estimated by randomly drawing errors from a double exponential 67  probability distribution with shape factors for daytime and nighttime measurements from Richardson et al. (2006). We found no relative biases in annual sums for each year. 3.3.3.3 Data analysis methods Regression analysis between C fluxes and environmental variables at a monthly time scale was chosen to detect ecosystem responses to changes in driving variables because of the pronounced spectral gap on the time scale of a month (Baldocchi, et al. 2001b). We used least squares, single factor linear regression and residual analysis to evaluate the dependence of NEP, Re and GPP on several climatic and biophysical variables. For each dependent variable, a residual analysis with other independent factors was followed by one-single factor regression analysis. The multiple regressions were run for all possible combinations of the independent variables, and the best-fit model was selected as the one that produced the lowest value (the probability of a type-one error). For the regression analysis, we determined the probability of significance, the p-value, using the F-test or t-test, at the significance levels of 0.05 or 0.01. We report the p-value and the coefficient of determination (r 2 ). If the p-value is less than 0.001, we show it as p < 0.001. 3.4 Results and discussion 3.4.1 General weather conditions The general annual weather pattern on the west coast of Canada includes a warm dry season (April – September) and a mild wet season (October – March) (Figure 3.1 and Table 3.1). Monthly mean Ta varied between 1 and 9 ºC in the wet season and 10 and 18 ºC in the dry season with relatively few occurrences of freezing or extremely high temperatures. The annual sums of hydrological-year P, which was calculated from October 1 of the preceding year to September 30 of the nominal year, ranged from 935 mm (2001) to 1777 mm (1999) with most of it (~77%) occurring in the wet season (ranging from 59% (2005) to 87% (2006)) with < 60 mm occurring during June through September. Low P and high Ta in the dry wseason resulted in high D (Figure 3.1e) and low soil water content of the 0-60 cm soil layer (θ0-60cm < 0.2 m 3  m -3 , Figure 3.1d) during July to October (this period is hereafter 68  referred as the dry-soil months and the remainder of the year is referred as the wet-soil months). Generally most of the available water in the 0-60 cm soil layer was removed during July-October (Figure 3.1d). The beginning of the study period coincided with an El Niño/La Niña cycle during the northern hemispheric winters of 1997/1998 and 1998/1999. As a result, the warmest and coldest spring and summer temperatures occurred in 1998 and in 1999, respectively. The 2002-2003 El Niño event did not significantly affect temperatures and precipitation at the study site. As shown in Figure 3.2 and Table 3.2, Ta in the region has gradually warmed over the past 40 years (1965-2004), with a slightly higher warming rate in the dry season than in the wet season (0.040 vs. 0.037 ºC yr -1 ). During the study period 1999-2006 (excluding the 1998 El Niño year), spring and the dry season continually warmed following the past 40-year trend, while winter and the wet season showed no trend.   69  Table 3.1 Climate statistics for the research site  1998 1999 2000 2001 2002 2003 2004 2005 2006 Avg. Mean temperature (ºC)   Annual 9.07 7.62 8.20 8.06 8.45 8.44 8.75 8.28 8.36 8.36   Dry 14.08 11.73 12.62 12.16 12.81 13.05 13.72 12.72 13.39 12.92 Daily minimum temperature (ºC)   Annual 6.37 5.10 5.71 5.53 5.68 5.77 6.10 5.72 5.57 5.73   Dry 10.35 8.43 9.26 8.93 9.12 9.52 10.12 9.42 9.71 9.43 Hydrological-year precipitation (mm)   Annual 1432 1777 1145 935 1249 1277 1349 1352 1699 1357   Dry 216 290 284 278 217 313 375 554 225 305 Soil water content (0-60 cm, m 3  m -3 )   Annual 0.182 0.220 0.208 0.209 0.198 0.205 0.218 0.223 0.206 0.21   Dry 0.166 0.203 0.187 0.185 0.172 0.173 0.188 0.202 0.173 0.18 Dry season length (days)   Annual 104 87 100 117 158 118 76 65 124 105 Dry (dry season) extends from April to September. The hydrological-year is measured from October of preceding year to September of the nominal year. Dry season length is the number of days when daily maximum soil water content of 0-60 cm soil layer (θ0-60cm) was consistently below 0.20 m3 m-3. Avg.—is the 9-year average value  The 40-year normal hydrological-year precipitation was 1418 mm with a standard deviation of 239 mm. No significant trend in hydrological-year precipitation was measured over the past 40 years and during the study period. However, autumn and the dry season showed a significantly increasing trend during the late 1970s to 1998 (Figure 3.2 and Table 3.2).  70   Figure 3.1 Seasonal and interannual variability of climatic and environmental variables measured at the research site for 1997-2006: (a) Monthly totals of precipitation (P), (b) air temperature (Ta) above the canopy, (c) soil temperature (Ts) at 5-cm depth, (d) soil water content of 0-60 cm soil layer (θ was weighted by the distribution of roots for the top 0-60 cm), (e) vapor pressure deficit (D) above the canopy, and (f) photosynthetically active radiation (Q) above the canopy.  The dotted and dashed lines in panel (e) are the estimated soil water content at field capacity (soil water matric potential = -0.01 MPa) and permanent wilting point (soil water matric potential = -1.5 MPa), respectively (Coops et al. 2007b). 71   Figure 3.2 Climatic features for the period 1965-2004, observed at Campbell River airport climate station (Latitude: 49º 57.000’N, Longitude: 125º16.200’W, Elevation: 105.5 m), which is 10 km NE of the research site (DF49), showing interannual and seasonal trends and anomalies. (a) Air temperatures were calculated on the basis of calendar year: the dry season extends from April - September while the remainder of the year is the wet season. (b) The winter, spring, summer and autumn seasons were classified following Mote (2003) as December-January-February, March-April-May, June-July-August and September-October- November, respectively. Hydrological-year precipitation was calculated from October of the preceding year to September of the nominal year. (c) The hydrological year was broken down into the wet and the dry seasons. The dotted lines show the 40-year mean values for each season and the thick lines show measured air temperatures above the canopy and precipitation at DF49. 72   Table 3.2 Climate statistics for a long-term climate station near the research site * * The station is at the Campbell River airport located at (49º 57.000’N, 125º16.200’W), 10 km northeast of the research site. Air temperatures were calculated on the basis of calendar year: Dry season extends from April – September while the remainder of the year is the wet season. The winter, spring, summer and autumn seasons were classified following Mote (2003) as December-January-February, March-April-May, June-July-August and September- October-November, respectively. Precipitation was calculated on the basis of “hydrological- year” from October 1 of preceding year to September 30 of the nominal year. σ is the ±1 standard deviation calculated using yearly values. Linear correlation coefficient (r) and statistically significant level (p) values are for the temperature or precipitation trends with time at yearly time steps. The bold numbers indicate those are statistically significant (p < 0.05) 3.4.2 Diurnal and seasonal CO2 exchange The interannual variability in the 9-year ensemble monthly composite diurnal cycles of NEP and GPP for alternative 6 months is shown in Figure 3.3. Daytime CO2 uptake occurred throughout the year. High interannual variability in GPP and daytime NEP occurred during  Annual Dry Wet Winter Spring Summer Autumn Air temperatures, 1965-2004 Mean (ºC) 8.674 13.462 3.886 2.757 9.883 15.904 6.129 Stdv (ºC) 0.691 0.766 0.817 1.658 2.111 0.831 1.842 Slope(ºC yr -1 ) 0.039 0.040 0.037 0.083 0.151 0.021 -0.098 R 0.656 0.615 0.534 0.585 0.837 0.297 -0.619 P <0.001  <0.001 <0.001 <0.001 <0.001 0.062 <0.001 Precipitation,  hydrological-year 1965-2004 Mean (m) 1.418 0.381 1.038 0.509 0.235 0.137 0.539 σ (m) 0.239 0.127 0.211 0.172 0.071 0.054 0.147 Slope (m yr -1 ) 0.004 0.005 -0.002 -0.002 -0.001 0.001 0.007 R 0.170 0.471 -0.091 -0.142 -0.214 0.112 0.535 P 0.294 0.002 0.577 0.383 0.185 0.493 <0.001 73  the mid-day hours for all months. As expected, interannual variability in nighttime NEP (i.e. Re) in the dry-season months was much higher than that in the wet-season months.  Figure 3.3 Nine-year averaged monthly composite diurnal cycles of net ecosystem productivity (NEP, thin solid line) and gross primary productivity (GPP, thick solid line). Vertical bars show ±1 standard deviation of annual variation. The general seasonal pattern of monthly NEP for the study period can be summarized as follows (Figure 3.4a): small C losses (5-10 g C m -2  month -1 ) occurred during November to January; C uptake started in February (20-30 g C m -2  month -1 ), peaked either in April or May (~100 gC m -2  month -1 ) and gradually decreased to August; small C sinks or sources occurred in September and October. Monthly GPP quickly increased beyond the annual monthly average of GPP (~180 g C m -2  month -1 , Figure 3.4b) in spring in response to high Q and sufficient soil water supply (March to June, Figures 3.1f & d) while monthly Re in spring was less than its annual monthly average of Re (~150 gC m -2  month -1 , Figure 3.3c) as a result of low Ta (Figure 3.1b). Both 74  monthly GPP and Re reached their maxima in July with values of 350 and 315 gC m -2  month - 1 , respectively.  Figure 3.4 Monthly totals of net ecosystem productivity (NEP), gross primary productivity (GPP) and ecosystem respiration (Re) from November 1997 to the end of 2006. Figure 3.5a shows that spring to early summer (March-June) accounted for more than 80% of the annual NEP while late spring to early autumn accounted for more than 60% of annual GPP and Re (May-August and June-September for GPP and Re, respectively). However, the highest interannual variability in NEP for a 4-month period was observed during late spring to summer (the standard deviation σ for the period of May-August was ~80% of σ for the full year, Figure 3.5b) while high interannual variability in GPP and Re for a 4 or 3 month period was observed during spring to early summer (the σ value for the period of March-July was ~90% of the σ value for the full-year course, Figure 5b). High inter-annual variability of 75  monthly NEP, GPP and Re over the 9 years (high monthly values of σ, which was calculated from the monthly flux values) was observed during May-September with Re showing higher inter-annual variability than GPP (Figure 3.5c).  Figure 3.5 Seasonality of C fluxes and the corresponding interannual variability for 1998- 2006. (a) Time series of the 9-year-average value of the ratio of the 4-month (current month and the 3 following months) total C flux (NEP, GPP and Re) to the annual total C flux for each month of the year. (b) Time series of the 9-year-average value of the ratio of the 4- month (current month and the 3 following months) average value of σ for monthly C fluxes to the annual value of sigma for each month of the year. (c) Standard deviation (σ) of monthly C fluxes. (d) Standard deviation of monthly averages of environmental variables. NEP is net ecosystem productivity, GPP is gross primary productivity, Re is ecosystem respiration, Ta is the air temperature above the canopy, θ is volumetric soil water content of 0-60 cm soil layer weighted by the distribution of roots and D is vapor pressure deficit above the canopy. The gray area (May - September) in panels (b) and (d) shows the period with high σ values of C fluxes. The GPP/Re ratio is a useful diagnostic in studies of seasonal and inter-annual variability (Valentini et al., 2000; Law et al., 2002). Figure 3.6 shows the seasonal variations in GPP/Re. An unusually low March ratio was observed in 2002, whereas higher than normal ratios were observed in the spring and early summer of 1999. To project the inter-annual variability, the 76  departures in monthly and annual C fluxes were calculated by removing the 9-year mean. Large positive anomalies in monthly GPP (ΔGPP) and Re (ΔRe) were observed in the springs and summers of 2004 and 2005 and the autumn of 1999, whereas low departures were observed in the springs of 1999 and 2002 (Figure 3.7). Departures in monthly NEP (ΔNEP) from the 9-year mean negatively were significantly negatively correlated with ΔRe and ΔGPP (Figures 3.7b, c). This was similar to observations made in an old-growth forest by Falk et al. (2008).  Figure 3.6 Monthly gross primary productivity to ecosystem respiration (GPP /Re) ratios for 1998 - 2006 and the 9-year mean. 77   Figure 3.7 Monthly departures of gross primary productivity (GPP), ecosystem respiration (Re), and net ecosystem productivity (NEP) fluxes from the 1998 - 2006 means. Linear fits are shown as dashed lines for each panel. 78  The coefficient of variation (CV) was calculated from the monthly mean values to assess seasonal variability of C fluxes and several associated climate variables (Table 3.3). The seasonal variability of NEP (9-year average value of CV = 1.42) was significantly greater than that of GPP (9-year average value of CV = 0.67) and Re (9-year average of CV = 0.73). The interannual variability of NEP, GPP, Re and all associated climate variables was remarkably less than their seasonal variability (Table 3.3). Table 3.3 Seasonal and interannual variability of C fluxes and associated climate variables expressed as the coefficient of variation (CV), 1998 - 2006*  1998 1999 2000 2001 2002 2003 2004 2005 2006 Avg IAV NEP 1.18 1.42 1.31 1.36 1.70 1.15 1.76 1.47 1.39 1.42 0.21 GPP 0.69 0.67 0.66 0.64 0.67 0.66 0.70 0.67 0.67 0.67 0.06 Re 0.74 0.82 0.68 0.74 0.66 0.73 0.74 0.71 0.74 0.73 0.09 Ta 1.04 0.82 0.59 0.81 0.91 0.93 0.54 0.65 0.93 0.80 0.19 Ts 1.15 0.82 0.71 0.57 0.84 0.96 0.68 0.55 0.88 0.80 0.18 P 0.66 0.70 0.67 0.64 0.67 0.70 0.68 0.63 0.74 0.68 0.05 Phydro 0.82 0.77 0.74 0.62 0.81 0.85 0.92 0.70 0.75 0.78 0.15 D 0.62 0.74 0.65 0.60 0.62 0.57 0.56 0.54 0.67 0.62 0.08 Ψ30-50cm 0.94 1.09 1.00 0.93 0.88 1.01 1.08 1.16 1.02 1.01 0.24 θ0-60cm 0.29 0.19 0.21 0.21 0.27 0.25 0.19 0.16 0.27 0.23 0.07 *IAV was calculated from annual means, in dictating interannual variability. Seasonal CV was calculated using monthly average values and the 9-year average values (Avg) of seasonal CV are also shown. Variables include net ecosystem productivity (NEP), gross primary productivity (GPP), ecosystem respiration (Re), calendar year precipitation (P), hydrological-year precipitation (Phydro), air temperatures (Ta) above canopy, air vapor pressure deficit (D) above the canopy, soil temperatures (Ts) at 5-cm depth, soil water matric potential of 30-50-cm layer (Ψ30-50cm) and soil water content of  0-60-cm layer (θ0-60cm). Phydro was calculated from October of preceding year to September of the nominal year 3.4.3 Climatic controls on seasonal CO2 exchange To understand the mechanisms behind the relationship between ecosystem C sequestration and climatic factors, we examined how GPP, Re and NEP were related to changes in Q, Ta, water stress and other factors over the course of the year.  79  3.4.3.1 Response of CO2 exchange to Q and Ta As shown in Figure 3.8a, monthly mean GPP was linearly related to Q even though the light response on the half-hourly time scale was well described by a non-linear Michaelis-Menten model (Humphreys et al., 2006). Although monthly mean GPP was significantly correlated with Q (r 2 = 0.85, p < 0.001), the residuals of the GPP-Q response (GPP, residual), which are the differences between the observed GPP and the corresponding values predicted by the linearly-fitted GPP-Q model, were still considerable (see the scatter around the regression line in Figure 3.8a). Figure 3.8b shows that the correlation of GPP, residual to θ0-60cm was statistically significant; it was positive during the dry-soil months (July - October) but it was negative during the wet-soil months. GPP, residual was not correlated with either Ta or vapor pressure deficit (D). This suggests that during dry-soil months soil drought stress altered the response of photosynthetic assimilation to Q. This also implies that forest growth was water- limited during the dry-soil months (see Section 3.4.3.2). The combination of Q and θ explained 91% of the variance of monthly GPP over the 9-year period. As shown in Figure 3.8c, the temperature response of Re was exponential even when data were integrated at a monthly scale (see also Morgenstern et al., 2004). Does the exponential relationship between the monthly mean of Re and Ta (Figure 3.8c) really reflect their inherent functional relationship? To answer this question, we examined the relationship between chamber-measured soil respiration (Rs) and Ts since it is a large component of Re (Humphreys et al., 2006). This relationship was also exponential even after averaging the data to the monthly time scale (Figure 3.8e). This supports the exponential relationship between Re and Ta at monthly time scales. The variance shown by residuals of the Re-Ta response (Re, residual) was greater during the dry- season months than during the wet-season months (Figure 8c). Most of positive values of Re, residual occurred when θ0-60cm was 0.12  between 0.21 m 3  m -3 . The relationship of Re,residual to monthly mean θ0-60cm during the dry-soil months (July - October) was positive whereas during the wet-soil months it was negative, and both were statistically significant (Figure 3.8d). This suggests the optimal soil water content for Re was  between 0.12 m 3  m -3 and 0.21 80  m 3  m -3 , and Re, residual significantly decreased with increasing θ0-60cm when θ0-60cm > 0.21 m 3  m -3  or with decreasing θ0-60cm when θ0-60cm < 0.12 m 3  m -3 . The combination of Ta and θ explained 96% of variance of monthly Re. Table 3.4 Characteristics of linear regression analysis [y = ax + b] of monthly mean ecosystem respiration (Re, in μmol m -2  s -1 ) and gross primary productivity (GPP, in μmol m-2 s -1 ) vs. monthly mean air temperature (Ta, in ºC), and daytime mean photosynthetically active radiation (Q, in μmol m-2 s-1) for individual month and annual clusters. Data are from Nov 1997 to Dec 2006*   Jan Feb Mar Apr May Jun July Aug Sep Oct Nov Dec Annual Re vs. Ta r 0.19 0.60 0.89 0.82 0.82 0.88 0.75 0.09 -0.43 0.56 0.59 0.62 0.95 a 0.019 0.131 0.309 0.581 1.064 0.854 0.665 0.350 -0.491 0.403 0.103 0.215 0.569 b 1.31 1.16 0.87 -0.89 -5.04 -3.01 -0.59 3.73 13.64 0.67 1.90 1.08 -0.14 GPP vs. Ta r 0.73 0.70 0.93 0.76 0.81 0.64 0.70 0.05 -0.64 0.36 0.24 0.56 0.92 a 0.165 0.177 0.525 0.617 0.606 0.494 0.454 -0.111 -0.531 0.260 0.037 0.084 0.611 b 0.86 1.79 1.82 1.93 2.40 3.26 3.61 11.44 14.00 1.75 1.80 0.97 0.44 GPP vs. Q r 0.41 0.12 0.51 0.21 0.18 0.48 0.23 -0.63 -0.11 0.89 0.58 0.24 0.92 a 0.005 0.001 0.012 0.002 0.005 0.005 0.003 -0.009 -0.001 0.006 0.004 0.001 0.0145 b 0.42 2.06 -1.17 4.96 5.77 6.69 8.96 15.70 7.64 1.94 1.24 0.97 -1.10 *The bold numbers indicate those are statistically significant (p < 0.05) and r is the correlation coefficient  None of the relationships between the C fluxes and climate driving variables derived from the individual month clusters (Table 3.4) followed the annual relationships (Figure 3.8) probably because of the different variance and covariance at different time scales in the records of C fluxes (Baldocchi, et al. 2001b). The monthly clusters showed very different sensitivities of Re and GPP to Ta and Q (Table 3.4): (1) Both Re and GPP were significantly correlated with Ta during the spring (starting from winter for GPP) to middle summer months. (2) The sensitivities (regression slopes) of both Re and GPP to Ta increased steadily from winter to late spring and decreased in early summer, with maxima of 1.604 μmol m-2 s-1 ºC -1  in May and 0.617 μmol m-2 s-1 ºC-1 in April, respectively; (3) the response of GPP to Q was not statistically significant for almost all months in contrast to the annual relationship. 81   Figure 3.8 Photosynthetically active radiation (Q) and air temperature dependences of gross primary productivity (GPP) and ecosystem respiration (Re) rates over annual cycles, at DF49, Nov. 1997 - Dec. 2006. (a) Relationship between monthly means of gross primary productivity (GPP) and daytime (Q > 5 μmol m-2 s-1) Q, and the solid and dashed lines are the linear regression fits of EC-derived GPP (sum of gap-filled daytime net ecosystem productivity and calculated daytime Re and modeled GPP, respectively; (b) Relationship 82  between residuals of monthly GPP and monthly mean soil water content (θ, at the 0-60-cm depth), and the solid and dashed lines are the linear regression trends for dry-soil months and wet-soil months, respectively; (c) Responses of monthly mean Re to Ta , and the solid and dashed lines are the exponential regression fits of gap-filled measured Re and modeled Re using, respectively; (d) Relationships of residuals of monthly Re with monthly mean θ, and the solid and dashed lines are the linear regression trends for dry-soil months and wet-soil months, respectively;  (e) Responses of monthly mean soil respiration (Rs) to soil temperature (Ts) at the 5-cm depth. The residuals of GPP or Re are the differences between the observed values and the corresponding values that are predicted by the fitted models (GPP - Q model or Re-Ta relationships, respectively) and thus they represent the variance that is not explained by the relationships (see text). The dry-soil months are July-October and the remainder of months is the wet-soil months.  The difference in the relationship between GPP and Q and between Re and Ta, along with the seasonal variation patterns of Q and Ta (Figures. 3.1f, b) led to a delicate balance in their influences on NEP. The responses of NEP to changes of Q and Ta showed a seasonal cycle with three distinctive stages: January to April, May to August, and the remainder of the year (Figure 3.9, Table 3.5). The responses of NEP to changes of Ta and Q were positive during the first and last four months of year while they were negative during the middle four months of year, and NEP was more sensitive to changes in Ta and Q during January to April than that during September to December. These seasonal differences in the relationship between NEP and Ta and Q possibly resulted from: (1) the differences in the tree phenology of photosynthesis and autotrophic respiration in the first and the second halves of the year; (2) the difference in seasonal distributions of Ta and Q: Ta was generally higher during January - April than during September - December (Figure 3.1b), whereas Q was generally lower during January - April than during September - December (Figure 3.1f); (3) the difference in the sensitivities of GPP and Re to changes in Ta (linear and exponential, respectively); and (4) the different decreases of GPP and Re resulting from water deficiency in summer and autumn (see Section 3.4.3.2). 83   Figure 3.9 Net ecosystem productivity (NEP) vs. air temperature (Ta) above the canopy (a) and daytime photosynthetically active radiation (Q) above the canopy (b) on monthly time scales. The data shown are for Nov 1997 - Dec 2006. The legend is the same as in Figure 3.8. The results of linear regression analyses for the same data are shown in Table 3.5. Table 3.5 Characteristics of linear regression analysis [y = ax + b] of monthly mean net ecosystem productivity (NEP, in μmol m-2 s-1) vs. monthly mean air temperature (Ta, in ºC) and daytime mean photosynthetically active radiation (Q, in μmol m-2 s-1) for each 4-month period. Data are from Nov 1997 to Dec 2006* Seasonal periods NEP vs. Ta NEP vs. Q a b r p a b r p Jan - Apr 0.49 -0.5 0.88 <0.001 0.007 -1.25 0.94 <0.001 May - Aug -0.351 6.52 -0.85 <0.001 -0.004 4.57 -0.34 0.04 Sep - Dec 0.035 -0.38 0.43 <0.005 0.001 -0.43 0.52 <0.001 *r is the linear correlation coefficient and the p-value using F-test indicates the probability of significance. The bold numbers indicate the values of slopes are negative. The observation numbers n = 37 or 39 for the sampling group of Sep - Dec or for the groups of Jan - Apr and May - Aug, respectively)  84   Figure 3.10 Relationships of monthly mean daytime net ecosystem productivity (NEPd) vs. monthly mean soil water content (θ30-50cm) at the 30-50-cm depth (a) and vapor pressure deficit (D) above the canopy (b). Linear regression fits applied to the dry season months (April - September). The data shown are for DF49, Nov. 1997 - Dec. 2006. 3.4.3.2 Response of CO2 exchange to water stress A decrease in midday latent heat flux was observed when D was high and the soil was dry during the dry season (not shown). This decrease coincided with low values of canopy conductance calculated by inverting the Penman-Monteith equation (Monteith and Unsworth, 1990), indicative of a decrease in midday stomatal conductance in the dry season. Photosynthesis rate was also depressed during midday hours likely as a result of stomatal closure and/or a directed effect of the drop in the leaf water potential due to soil water stress in the dry season. This was also found in an old-growth Douglas-fir stand in south western Washington: photosynthesis rates were attenuated by high D in the dry season (Paw U et al., 2004; Chen et al., 2004). Coops et al. (2007a) reported the limitations on growth of coastal Douglas-fir in Vancouver Island mainly caused by suboptimal temperature and high D based on a physiology-based forest-growth model analysis and remotely-sensed data. On the other hand, dry-soil conditions can also result in a decrease of Re possibly by affecting enzyme activity, diffusion of solutions and O2, growth of root tissue, and microbial populations 85  (Keith et al. 1997; Davidson et al., 1998; Lavigne et al. 2004; Hutyra et al., 2007). Falk et al. (2005, 2008) reported that summer soil-water stress attenuated annual Re by about 13% in an old-growth Douglas-fir stand in southwestern Washington. Similar attenuation of Re was also observed in this Douglas-fir stand (Jassal et al., 2007). In theory and many process-based models (e.g. Leuning et al., 1998; Grant et al., 1999; Chen et al. 2007a,b), water deficiency (drought) would simultaneously provide a strong constraint on both GPP and Re. Declines in GPP and Re were found to occur simultaneously with the heat wave and drought of 2003 in Europe (Reichstein et al., 2006). But how does drought influence NEP? Is GPP or Re more greatly depressed by moisture limitation in the dry season? As shown in Figure 3.10b, in this study, daytime NEP significantly decreased with increasing D in the dry season (p < 0.001), suggesting GPP was constrained by a decline of stomatal conductance when the ambient air in the canopy was dry (high D). The increase of NEP with increasing D in the wet season (Figure 3.10b) likely reflects the increase of NEP is due to increasing Ta as D simultaneously increased with Ta in the wet season.  Data also showed that daytime NEP significantly decreased with decreasing θ30-50cm in the dry season (Figure 3.10a, p < 0.001). This suggests that moisture limitation in the dry season depressed GPP greater than Re. Figure 3.11 shows that both GPP and Re decreased with increasing soil water stress (i.e. decreasing θ) under drought conditions (monthly mean θ30-50cm < 0.1 m 3  m -3 ), however the former had a steeper slope (a t-test indicated that the regression slope of GPP with θ30-50cm significantly differed from that of Re, p < 0.001). Daytime NEP significantly decreased with soil water stress (p < 0.025).  Strong moisture limitation depressed or even disrupted the relationships (between GPP and Ta, between GPP and Q, and between Re and Ta: the coefficients of determination declined considerably under severe drought conditions, not shown). 86   Figure 3.11 Relationships of the monthly mean gross primary productivity (GPP), daytime ecosystem respiration (Red), and daytime net ecosystem productivity (NEPd) with soil water content (θ30-50cm) at the 30-50-cm depth under severely dry conditions (θ30-50cm < 0.1 m 3  m -3 ). Solid, long-dashed, and short-dashed lines are linear fits for GPP, Red, and NEPd vs. θ30-50cm, respectively. 3.4.4 Dependence of the annual NEP on GPP and Re Average annual values (± stand deviations) of NEP, GPP and Re for the 9 years were 357 ± 51, 2,124 ± 125, and 1,767 ± 146 g C m -2  yr -1 , respectively (Table 3.6). The departures of annual NEP (ΔNEP) from the 9-year mean ranged from -90 to 53 g C m-2 yr-1 (Table 3.6). The two years that can be considered to be anomalies are 2002 and 2004 (|ΔNEP| > 80 g C m -2  yr -1 ). Departures of GPP (ΔGPP) and Re (ΔRe) from the 9-year respective means are shown in Figure 3.12. Suppressed photosynthesis in 2002 (ΔGPP = -173 g C m-2 yr-1, ΔRe = - 92 g C m -2  yr -1) and enhanced respiration in 2004 (ΔRe = 304 g C m -2  yr -1, ΔGPP = 215 g C m -2  yr -1 ) caused these two years to have the lowest NEP of the 9 years (Table 3.6, Figure 3.12). 2002 was a very dry year with the longest dry-season length (158 days), while 2004 had sufficient precipitation but was a very warm year. Both years had the fewest days with a positive NEP both annually and during the dry season (Table 3.6). 2002 had the warmest wet season of the 9 years (Figure 3.1b), resulting in high wet season Re (Figure 3.4c) while GPP was still close to the 9-year average (Figure 3.4b) due to the moderately low Q (Figure 3.1f) 87  which limited photosynthesis even though Ta was higher than the 9-year average. This resulted in the lowest wet season NEP (-12 g C m -2 ), the only wet season losing C of the 9 years. Table 3.6 Annual sums of net ecosystem productivity (NEP), gross primary productivity (GPP), and ecosystem respiration (Re) in g C m -2  and number of C sink days (CSD) during the year and the dry-season months (Apr - Sep) when the daily GPP to Re ratio was greater than 1   1998 1999 2000 2001 2002 2003 2004 2005 2006 Avg NEP Annual 379 382 400 410 277 353 267 355 386 357 Dry 319 303 358 326 289 278 206 291 298 296 GPP Annual 2131 2024 2091 2077 1952 2078 2338 2310 2112 2124 Dry 1694 1575 1643 1592 1541 1634 1873 1832 1659 1671 Re Annual 1752 1642 1693 1667 1676 1725 2071 1955 1727 1768 Dry 1375 1272 1286 1265 1253 1356 1667 1541 1361 1375 CSD Annual 250 248 245 260 226 260 226 230 253 244 Dry 150 141 149 145 137 143 127 138 136 141  What are the sources of the variation in the annual C budget ---as reflected in GPP and Re? There is debate about whether photosynthesis or ecosystem respiration contributes more to inter-annual variability of NEP (Ma et al., 2007). Valentini et al. (2000) found that Re was the main determinant of variation of NEE in European forests. Similarly, Hutyra et al (2007) found that climate anomalies exerted a strong influence on NEE in Amazonian rain forest principally through effects on Re according to 5 years (2002 - 2006) of EC CO2 flux data. On the other hand, 13 years (1992 - 2004) of CO2 flux measurements in a temperate mixed deciduous forest site at Harvard forest showed that NEE exhibited significant inter-annual variability due to variations in GPP (Urbanski et al., 2007). Reported anomalies in Re and GPP in the literature are generally positively correlated for the majority of ecosystems (e.g. Ciais et al. 2005; Falk et al. 2008). ΔGPP and ΔRe were highly correlated on both annual (Figure 3.12) and monthly (Figure 3.7a) time scales, indicating that Re were coupled to GPP, i.e. they both increased or decreased depending on climate perturbations. However, the high correlation between Re and GPP was not true for all years 88  (see the outliers of 2002 and 2004 in Figures 3.7a, 3.12). The difference between increasing or decreasing rates of Re and GPP due to climate perturbations resulted in anomalies in annual C budget, e.g. 2002 and 2004.  Figure 3.12 Departures in annual gross primary productivity (GPP) and ecosystem respiration (Re) fluxes from their respective 1998 - 2006 means. The linear fit is shown as a dashed line while the 1:1 line is shown as a dotted line. In the very dry year of 2002 (Figures 3.1d, e & Table 3.1), both GPP and Re were depressed to the lowest values over the 9 years of 1,952 and 1,676 gC m -2 (Table 3.6),  respectively. The greater depression of GPP by drought than Re in this year led to the second lowest C sequestration year in the study period. In the very warm year of 2004, both GPP and Re were enhanced by high temperature to the largest values over the 9 years of 2,338 and 2,071 gC m -2  (Table 3.6), respectively. The increase of Re due to high Ta was greater than that of GPP in this year because Re was more sensitive to Ta than GPP, especially at high Ta (see Figure 3.8, GPP linearly correlated with 89  Ta vs. Re exponentially correlated with Ta).  As a result, 2004 was the lowest NEP year in the study period (Table 3.6). These 9-year observations confirm that both respiratory and photosynthetic fluxes were enhanced by increasing air temperature in years with sufficient water supply (normal precipitation years) and Re was enhanced more than GPP (e.g. 2004); whereas both GPP and Re were suppressed by water stress in drought years and GPP was suppressed more than Re (e.g. 2002). Annual net C sequestration was maximized in optimally warm and normal precipitation years (e.g. 2001 and 2000), but minimized in unusually warm (e.g. 2004) or severely dry (e.g. 2002) years. 3.4.5 Climatic controls on annual and inter-annual CO2 exchange Although responses of plant and ecosystem physiological processes to weather are expected, our data showed poor correlations of C fluxes with single climatic factors (i.e. Ta, Q, P, θ, , D, etc) on yearly time scales. Our results, however, show that the changes in C fluxes depended more on seasonal climatic variability (especially in the dry season). NEP was negatively correlated with Ta for annual (r = -0.60) and all seasons (|r| < 0.50) but not significantly. Excluding the two years strongly affected by an El Niño/La Niña cycle (1998-1999), annual NEP was significantly correlated with mean annual Ta (r = -0.96, p <= 0.001). Although autumn NEP was significantly correlated with autumn Ta (r = -0.77, p = 0.04), annual NEP was not correlated with it. This finding is inconsistent with that of Piao et al. (2008) that autumn warming led to net C losses in northern ecosystems. Spring Ta explained 89% and 80% of the variance in annual GPP and Re, respectively, but annual NEP was not correlated with spring Ta (Figure 3.13a). Annual GPP and Re were not correlated with annual P but were correlated with hydrological- year P (r = 0.58, p = 0.08 and r = 0.63, p = 0.07, respectively), and with dry season P (r = 0.71, p < 0.05, and r = 0.50, p = 0.14, respectively). Autumn NEP was significantly correlated with summer P in (r  = 0.80, p < 0.01). This finding of a lag of C fluxes with 90  respect to P is similar to that of Falk et al. (2008) that NEP was significantly correlated with precipitation that fell 100-105 days earlier in an old-growth Douglas-fir stand.  Figure 3.13 Responses of annual C fluxes (GPP: gross primary productivity, Re: ecosystem respiration, and NEP: net ecosystem productivity) to mean spring temperature (a, 1998 - 2006) and annual mean diffuse photosynthetically active radiation (Qdif) above the canopy (b, 2000 - 2006). 91  Multiple linear regression analysis revealed that annual NEP was significantly (r 2 = 0.82, p < 0.001) determined by the combination of  mean March Ta, May-June Ta, and mean soil water content of 0-60 cm soil layer for August-October. Canopy light use efficiency (LUE, which is the slope of the NEP vs. Q response curve) has been reported to be sensitive to whether the sky is clear or cloudy (Price and Black 1990; Hollinger et al. 1994; Baldocchi 1997a, 2008; Schwalm et al., 2007). Using the diffuse Q measurements made at the site since the beginning of 2000 we found that the correlations of annual GPP and NEP to annual mean diffuse Q (Qdif) were statistically significant (Figure 3.13b). This is similar to the result reported by Gu et al. (2002, 2003) and Niyogi et al (2004) that LUE nearly doubles when incident sunlight is diffuse, compared with when it is mostly direct. The general consensus (Baldocchi 1997a, 2008) is that under cloudy skies, incoming solar radiation is more isotropic, so that it penetrates deeper into the canopy. This process allows more light to reach and illuminate shaded leaves. On sunny days, the photosynthetic rate of upper sunlit leaves tend to be light-saturated. They also experience a higher heat load, which enhances respiration, and thus lowers their photosynthesis rates. Annual and spring temperatures, and hydrological-year and the dry-season precipitation have been steadily increasing over the past half century in the study area (see Figure 3.2 and Table 3.2) and in the Pacific Northwest region (Mote, 2003). Assuming these trends persist in the near future, dry-season NEP will likely increase, annual NEP will likely decrease, whereas both annual GPP and Re will likely increase. 3.5 Summary and conclusions (1) Daytime C uptake occurred throughout the year. The observed overall seasonal patterns of monthly total NEP for the study years were: (i) small C losses occurring during November to January; (ii) large C gains starting in February and continuing through August with a maximum value of ~100 gC m -2  month -1 occurring in either April or May; and (iii) small C gains or losses occurring in September and October. The annual average sums (± stand 92  deviations) of NEP, GPP and Re were 357 ± 51, 2,124 ± 125, and 1,767 ± 146 gC m -2  yr -1 , with ranges of 267 - 410, 1,592 – 2,338, and 1,642 – 2,071gC m-2 yr-1, respectively. (2) Although the majority of C sequestration (more than 80% of the annual total) occurred during March through June, May through August was responsible for about 80% of the inter- annual variability of NEP. Both Re and GPP were enhanced in warm years with sufficient water supply with Re being more enhanced than GPP; whereas both GPP and Re were suppressed by water stress in dry years with GPP being more suppressed. The major drivers of inter-annual variability of annual C fluxes were annual and spring mean temperatures, and water deficiency during late summer to autumn when productivity of this Douglas-fir stand was water limited. (3) Monthly GPP was linearly correlated with Q (r 2  = 0.85) and monthly Re was exponentially correlated with Ta (r 2  = 0.94). The combination of Q and θ explained 91% of variance of monthly GPP and the combination of Ta and θ explained 96% of variance of monthly Re. The responses of NEP to changes of Ta and Q were positive during the first and last four months of the year but were negative during the middle four months of the year. NEP was more sensitive to changes of Ta and Q during January to April than that during September to December. (4) Assuming future climate at the site follows a trend to that of the past 40 years, the dry season NEP will likely increase, annual NEP will likely decrease, whereas both annual GPP and Re will likely increase.     93  Chapter  4: Modeling to discern nitrogen fertilization impacts on carbon sequestration in a Pacific Northwest Douglas-fir forest in the first-post fertilization year 4  4.1 Synopsis This study investigated how nitrogen (N) fertilization with 200 kg nitrogen ha -1  of urea affected ecosystem carbon (C) sequestration in the first-post fertilization year in a Pacific Northwest Douglas-fir (Pseudotsuga menziesii) stand on the basis of multi-year EC and soil- chamber measurements before and after fertilization in combination with ecosystem modeling.  The approach uses a data-model fusion technique which encompasses both model parameter optimization and data assimilation, and minimizes the effects of interannual climatic perturbations and focuses on the biotic and abiotic factors controlling seasonal C fluxes using a pre-fertilization 9-year-long time series of EC data (1998-2006). A process- based ecosystem model was optimized using the half-hourly data measured during 1998- 2005, and the optimized model was validated using measurements made in 2006 and further applied to predict C fluxes for 2007 assuming the stand was not fertilized. The N fertilization effects on C sequestration were then obtained as differences between modeled (unfertilized stand) and EC or soil-chamber measured (fertilized stand) C component fluxes. Results indicate that annual net ecosystem productivity (NEP) in the first-post N fertilization year increased by ~83%, from 302 ±19 g m -2  yr -1  to 552 ± 36 g m -2  yr -1 , which resulted primarily from an increase in annual gross primary productivity (GPP) of ~8%, from 1938 ± 22 g m -2  yr -1  to 2095 ± 29 g m -2  yr -1  concurrent with a decrease in annual ecosystem respiration (Re) of ~5.7%, from 1636 ± 17  g m -2  yr -1  to 1543 ± 31 g m -2  yr -1 . Moreover, with respect to  4  Reprinted from Global Change Biology, Vol. 17(3), 1442–1460, [B. Chen], N.C. Coops, T.A. Black, R.S. Jassal, J.M. Chen, M. Johnson (2011), Modeling to discern nitrogen fertilization impacts on carbon sequestration in a Pacific Northwest Douglas-fir forest in the first-post fertilization year. Global Change Biology, 17(3), pages 1442–1460, (DOI: 10.1111/j.1365-2486.2010.02298.x) Copyright (2011), with permission from John Wiley & Sons.  Minor updates, rewording and corrections have been applied. 94  respiration, model results showed that the fertilizer-induced reduction in Re (~93 g m -2  yr -1 ) principally resulted from the decrease in soil respiration Rs (~62 g m -2  yr -1 ). 4.2 Introduction Global inputs to the terrestrial nitrogen (N) cycle have increased two- to five-fold (Vitousek et al., 1997; Galloway et al., 2002; Janssens et al., 2010) in the past century mainly due to human activities, particularly fertilizer use and fossil fuel burning. N addition to forest and other ecosystems has been hypothesized to have a positive or negative effect on the health and vitality of ecosystems and the terrestrial C cycle depending on the degree of ecosystem N saturation (Aber et al., 1989, 1993, 1998). Temperate and boreal forest ecosystems are generally N-limited and the N addition has been hypothesized to result in increasing carbon (C) sequestration (Aber et al., 1998). A range of  studies has shown that the effects of N addition are positive on forest growth and also mostly positive on soil C sequestration in the Northern Hemisphere (e.g. Vitousek and Howarth, 1991; Aber et al., 1995; Bergh et al., 1999; Franklin et al., 2003; Raey et al., 2008; De Vries et al., 2009; Liu and Greaver, 2009). Earlier modeling results suggested that N addition could account for an increased C sequestration of 0.44 - 0.74 Pg yr -1  by simply assuming that most (80%) of the deposited  N would be stored in plant tissues (Townsend et al., 1996; Holland et al., 1997) and these estimates  have been questioned to be overestimated (Nadelhoffer et al., 1999). Olsson et al. (2005) found that fertilization of a boreal Norway spruce stand led to a three-fold increase in aboveground net primary productivity (NPP), possibly due to decreased C allocation to roots in response to higher nutrient availability. Leggett and Kelting (2006) found that N fertilization of Loblolly pine plantations not only increased living biomass but also increased the size of soil C pools. The contribution of an average additional N deposition on C sequestration in European forests and soils in the period 1960–2000 was estimated to be 0.0118 G ton yr-1 by de Vries et al (2006), being equal to ~10% of their net C sequestration in that period (0.117 Pg C yr -1 ). Pregitzer et al. (2008) also reported that chronic N deposition increased C storage in northern temperate forests.  Magnani et al. (2007) even reported N deposition to be the dominant driver of C sequestration in forest ecosystems, which generated an intense debate about the 95  magnitude and sustainability of the N-induced C sink and its underlying mechanisms (de Vries et al., 2008; Sutton et al., 2008; Janssens and Luyssaert, 2009; Jamssens et al, 2010). However, other research results have shown that N addition to ecosystems had only a slightly positive effect on C sequestration (Nadelhoffer et al., 1999) or had almost no contribution to C storage (Korner, 2000), or in some cases reduced NPP and C storage at very high N inputs (hypothesized to be due to N saturation, e.g. Aber et al., 1989) or by other pollutants (e.g. acidification, Schulze, 1989). Bauer et al. (2004) reported that NPP and C sequestration may be enhanced when available C does not exceed the vegetation capacity for N uptake; otherwise, nutrient imbalances can lead to a decrease in C sequestration due to N saturation (Agren and Bosatta, 1998; Nihlgard, 1985).  Aber et al. (1989, 1998) hypothesized that forest ecosystems respond positively to N addition in the short-term and negatively in the long- term. Variable effects of N addition/fertilization to soils on soil organic C storage have been also observed as a result of enhanced, reduced or unchanged soil respiration (Rs) as a response to litter C/N ratio and/or a reduction in belowground allocation (Bowden et al., 2000, 2004; Burton et al., 2004; Cleveland and Townsend, 2006; Phillips and Fahey, 2007; Mo JM et al., 2006, 2008; Hyvönen et al., 2008). Pacific Northwest coastal forests of the United States and Canada cover approximately 10 5  km 2 between Oregon and Alaska and play a significant role in the global C cycle (Paw U et al., 2004; Falk et al., 2008; Chen et al., 2009b). In this region, there is very little N deposition from air pollution sources owing to their remote location, and the soils are generally considered deficient in N (Hanley et al., 1996). As a result, N fertilization is a common management practice in British Columbia (Brix, 1981; Fisher and Binkley, 2000; Chapin et al., 2002; Brooks et al., 2009) with a standard forest fertilization application rate of 200 kg N ha -1  from prilled urea at mid-rotation (Hanley et al., 1996). N fertilizer-induced decreases in Rs would increase soil C storage (Johnson and Curtis, 2001; Phillips and Fahey, 2007). N fertilization has been reported to be the only forest management activity having positive effects on the soil C pools, presumably as a result of a reduction in Rs (Johnson and Curtis, 2001). Since additional merchantable timber volumes can result from fertilization of stands of mid-rotation trees (i.e. of commercial thinning size, 20-40 year-old), N fertilization just 96  prior to harvest of near-end-of-rotation (50- to 60-year-old) Douglas-fir stands provides an attractive financial return with an average increase in bole volume growth of 20% (Hanley et al., 1996). This effect must result from imbalance between N-fertilization induced changes in the two ecosystem processes, gross primary productivity (GPP) and ecosystem respiration (Re). The objective of this study is to investigate how N fertilization affects ecosystem C component fluxes at daily, monthly and annual time scales in the first year following N fertilization, in a coastal Douglas-fir forest in British Columbia, Canada. While EC allows estimation of net C sequestration, it is different from traditional experiments directed to fertilization effects on stand growth. How to distinguish fertilization and climatic effects on C fluxes on the basis of EC measurements is a challenge. It is hard to find an area with similar land surface properties to set up a control EC tower because of high spatial heterogeneity of land surface. With absence of a control EC tower, a modeling approach to predict C fluxes for the fertilized year 2007 assuming the stand was not fertilized  is necessary for discerning the impact of fertilization on C sequestration. In this study, we apply a data- model fusion technique to optimize model parameters. This type of approach has recently been shown to be a powerful tool for minimizing the uncertainties in the land surface C flux estimates due to model parameter biases and measurement errors (Raupach et al., 2005; Sacks et al., 2006, Mo X et al., 2008) and allows us to integrate EC measurements, available soil-chamber measurements and a previously published process-based ecosystem model (BEPS: Boreal Ecosystem Productivity Simulator, see Liu et al., 2002, Ju et al., 2006, Chen et al., 2007a). 4.3 Material and methods 4.3.1 Site description This investigation occurred at a near-end-of-rotation costal Douglas-fir (Pseudotsuga menziesii) stand on the eastern side of Vancouver Island, British Columbia, Canada.  This area is in the very dry maritime subzone (CWHxm) of the coastal western hemlock (CWH) 97  biogeoclimatic zone (Meidinger and Pojar 1991). The CWH covers three million hectares in coastal and interior British Columbia, as well as parts of Alaska, Oregon, and Washington, USA. The EC tower (namely DF49, 4952 7́.8"N, 12520 6́.3"W) was established in 1997. The previous Douglas-fir stand at this site was harvested and slash-burned in 1943, followed by planting of Douglas-fir seedlings in 1949 resulting in a relatively homogeneous stand. The current second-growth stand surrounding the tower covers an area of 130 ha ranging from 300 – 400 m above sea level, with 80% Douglas-fir, 17% western red cedar (Thuja plicata), and 3% western hemlock (Tsuga heterophylla). Further details on soil and vegetation characteristics can be found in Morgenstern et al. (2004), Humphreys et al. (2006) and Chen et al. (2009b, c). 4.3.2 Stand fertilization About 110 ha surrounding the DF49 tower (including the EC flux footprint area of about 70 ha, which contributed >80% of the cumulative fluxes observed at the tower; see Chen et al., 2009b) was aerially fertilized with urea at 200 kg N ha -1  on Jan 13, 2007 using a Eurocopter SA315B helicopter (Western Aerial Applications Ltd., Chilliwack, British Columbia, Canada) with an in-house engineered hydraulic-driven spreader bucket and a GPS-assisted guidance system (Jassal et al., 2008). A non-fertilized area of about 17 ha (200 m × 850 m) on the southeast side of the fertilized area (500 m from the flux tower) served as a control for comparing differences in tree growth, C stocks and C fluxes. The location of the control area was chosen to be outside of the 95% cumulated tower footprint area (Chen et al., 2009b,c). Some of urea was retained in the snow-laden foliage on the day of fertilizer application, which was washed down to the ground surface with the melting of intercepted snow in the following days.  98  4.3.3 Needle mass and N concentration and growth increment measurements To assess the response to fertilization, the increase in size of the individual needles and the needle N content were measured and a preliminary tree ring analysis was made.  In the end of 2007, the foliage samples and bole cores at DBH (diameter at breast (1.3 m) height) from five representative healthy dominant and co-dominant trees in both the fertilized and the control areas were collected. These trees had a mean DBH of 39-cm. Samples were collected from the upper (5/6th), middle (3/6th) and lower (1/6th) sections of the live crown of each tree.  Current-year needles were carefully removed. The mass of 100 needles was determined for each of the trees. Dried needle samples were ground in an electric coffee grinder, digested in concentrated sulphuric acid and analyzed for total N content using an autoanalyzer (Autoanalyzer II, Pulse Instrumentation Ltd., SK, Canada). An image analysis system (model WinDENDRO, Regent Instruments Canada Inc., Nepean ON) was used to measure the annual ring widths to the nearest 0.01 mm. Allometric equations for coastal BC Douglas-fir were used to calculate mean annual (bole) increment (MAI) (Feller, 1992; Marshall and Turnblom, 2005). 4.3.4 Climate and EC measurements Details on the  climate and EC measurements and analysis for the site can be found in Morgenstern et al. (2004), Humphreys et al. (2006) and Chen et al. (2009a). Half-hourly net ecosystem exchange (NEE) was computed as the sum of EC-measured CO2 flux (Fc, positive upward) and the rate of change in CO2 storage (Sc) in the air column from the ground to the EC measurement height, i.e. cc SFNEE  . Net ecosystem productivity (NEP) was calculated as NEP = -NEE. The Fluxnet-Canada Research Network procedure for NEE gap-filling and partitioning of NEE into Re and GPP (Barr et al., 2004; Chen et al., 2009a) was followed except for using a logarithmic transformation  of an exponential (Q10) ReTs model rather than using a logistic model (see Chen et al., 2009a for additional details and rationale). The nighttime relationship 99  between log Re vs. Ts obtained using half hours with friction velocity greater than 0.35 m s -1  (i.e. the threshold value) was used in gap filling and estimating daytime Re. Random sampling of NEE error populations was used to determine the uncertainties in the annual sums of NEE. Each year was divided into 24 periods (12 months, day and night). Uncertainties in gap-filled half-hours were estimated by randomly drawing from the error population defined by half-hours with valid NEE observations. The mean difference between observations and estimates (NEEobserved  - NEEestimated) was less than 1 g C m 2  month 1  and the random error in the estimates of annual NEP was found to be within 30 g C m-2 (see also Morgenstern et al., 2004; Schwalm et al., 2007, Chen et al. 2009b). Two systematic biases of the analytical procedure due to variation of the threshold of friction wind speed and correction for the lack of energy balance closure were investigated (Morgenstern et al., 2004), the results showed no relative biases in annual sums of C fluxes for a single method between years. Standard errors (s.e.) of measured annual C component fluxes were calculated to be less than 50 g C m -2  yr -1  for GPP and 30 g C m -2  yr -1  for NEP and Re following the method of Schwalm et al. (2007). 4.3.5 Soil CO2 efflux data used in this study As part of an larger, ongoing soil respiration experiment (Jassal et al, 2010b), consecutive series of half-hourly soil respiration (Rs) data (2003-2007) measured using an automated closed-type CO2 efflux chamber were available and used for model calibration and validation in this study. 4.3.6 Modeling approach 4.3.6.1 Model description The process-based ecosystem model BEPS used in this study is an updated version by coupling with a land surface scheme (EASS: Ecosystem Atmosphere Simulation Scheme, see Chen et al., 2007a). It includes modules for photosynthesis, autotrophic respiration and live biomass allocation and soil C pools dynamics, soil biogeochemical and hydrological processes modules, and a scheme for the computation of energy balance, sensible and latent 100  heat fluxes, soil water and soil temperature status (Ju et al., 2006; Chen et al., 2007a, b). In the model framework, the canopy is stratified into overstory and understory layers, each of which is separated into sunlit and shaded leaf groups. A foliage clumping index, in addition to leaf area index (LAI), is used to characterize the effects of 3-dimensional canopy structure on radiation, water, and C fluxes.  In the soil hydrological processes module, snow packing and melting, rainfall infiltration and runoff, and soil vertical percolation are modelled. To estimate the vertical distribution of soil moisture and temperature, the soil profile is divided into seven layers and the thickness of the layers increases exponentially from the top layer to the sixth layer (0.05, 0.1, 0.2, 0.4, 0.8, 1.6 m, respectively). The first 6 soil layers with a total depth of 3.15 m are designed to ensure the complete simulation of energy dissipation in the soil column. The depth of the bottom soil layer is adjusted according to water table depth. Canopy fluxes of C, water, and energy are computed through stratification of the sunlit and shaded leaf components. A brief description of the processes directly related to photosynthesis and respiration is given in Appendix A. The time step of the model simulations is 30 min and model input data include atmospheric variables (temperature, relative humidity, wind speed, precipitation, and solar irradiance), vegetation type and stand age, canopy clumping index, soil texture and physical properties, and initial size of C pools. 4.3.6.2 Model experimental design As discussed, our approach utilises a model-data synthesis, which encompasses both model parameter optimization and data assimilation. To undertake the simulations and validate model predictions a data-model fusion technique was combined with an 8-years of EC- measured fluxes acquired at the DF49 site. We hypothesized that the underlying processes and multi-year seasonal patterns of C fluxes can be retrieved by assimilating a long time series of EC data into a process-based ecosystem model. The BEPS model was parameterized using the EC data obtained during 1998-2005. We assume that the optimized BEPS model captures the multi-year seasonal variations of C fluxes which we then tested by predicting C fluxes in 2006 and comparing them to actual observed data. Once we had confidence in the 101  BEPS model parameterization we simulated 2007 fluxes with differences between the observations and simulations expected to be principally the result by the N fertilization. 4.3.6.3 Model parameter optimization algorithm The Ensemble Kalman filter (EnKF) data-model assimilation technique, known as a stochastic-dynamic system, was used in this study. The basis of the EnKF technique is that a previous measurement can provide information about the state at the current time (Mo X et al., 2008). The optimum values of the model parameters, therefore, are assumed to correspond to the minima of the cost function J(x) (Tarantola, 1987),         b1bTb1oT xxPxxY(x)OCY(x)O 2 1 J(x)    ,                                            (4.1) where x is the vector of unknown parameters and xb is the a priori values of x, O is the vector of observations, Y is the nonlinear model (BEPS), Co is the covariance matrix of observations and Pb is the covariance matrix of a priori parameters. Following Diego et al. (2007), we adopted a gradient-based algorithm which converges more rapidly than standard Monte-Carlo methods to solve Eq. (4.1) for optimizing model parameters, typically converging to a minimum of J(x) within 100 iterations. 4.3.6.4 Parameter selection and ensemble generation The ensemble size (the number of parameters and the size of the moving window) is an important parameter in EnKF, which represents the number of model states predicted and analyzed concurrently. The size should be large enough to ensure the correct estimate of the error variance in the predicted model state (Williams et al., 2005). However, the very large ensemble size may be a heavy computation burden. Careful selection of which parameters are to be inversely optimized is required.  To identify which parameters are most sensitive to photosynthesis and respiration, we analyzed the responses of parameters to predicted C component fluxes by random sampling of parameters within their possible ranges. Seven parameters (Table 4.1), which are significantly sensitive to photosynthesis and respiration, were selected to be optimized. 102   Table 4.1 Model parameters optimized for a 58-year-old Douglas-fir stand Parameter symbols Parameter Description Units Prior Value Optimized values Range Mean s.d.  *  25 maxcV Maximum carboxylation rate at 25 °C (μmol m-2 s-1) 35 15-100 56.8 6 25 maxJ Maximum electron transport rate at 25 °C (μmol m-2 s-1) 70 30-180 89.6 8 m A dimensionless slope of the Ball-Berry model for predicting stomatal conductance - 6 3-14 7.3 0.5 D0 Sensitivity of stomatal conductance to water vapor saturation deficit kPa 1.6 1.1-2.5 1.9 0.16 fstress Soil water stress dependency of the sunlit/shaded leaf stomatal conductance - 0.5 0.05-1 0.44 0.27 Q10,s Temperature dependency of microbial respiration (K -1 ) 2 1-6 3.6 0.5 Q10,m Temperature dependency of maintenance respiration (K -1 ) 2 1-6 3.2 0.7 *  Standard deviation Based on a nine-year dataset analysis, Chen et al. (2009b) found that the seasonal variability of C component fluxes was significantly greater than their interannual variability. In addition, most ecological model parameters have been found to vary seasonally (Mo X et al. 2008). We therefore binned the 8-year dataset (1998-2005) into each calendar month with the ensemble size being the number of 8 days in a month   48 half-hourly data points per day. The selected parameters were optimized by minimizing the difference between observations and predictions, considering model and data uncertainties, and prior information on parameters.The model was then continuously run from the beginning of 1998 through to the end of 2007 at half-hourly time steps. 4.3.7 Data analysis methods Data analysis was conducted using Matlab (the Mathworks Company). Linear regression analysis between the measured C component fluxes and the corresponding model outputs across half-hourly to monthly time scales was chosen to test model behavior. Root mean square error (RMSE) was used to estimate the model errors. The significance levels of 103  differences between the two data series were detected using t-test. We determined the probability of significance, the p-value, using the F-test or t-test, at the significance levels of 0.05 or 0.01. We report the p-value and the coefficient of determination (r 2 ). If the p-value was less than 0.001, we only show it as p < 0.001. 4.4 Results 4.4.1 Comparison of measured environmental variables between pre- and post- fertilization year The climate variables (P, Ta, θ, D, and Q) in 2007 followed seasonal patterns similar to those observed during the previous non-fertilized 9-years (Figure 4.1). Compared to the 9-year means, annual total P was ~100 mm higher; annual Q was about ~7% lower, and annual average Ta was ~0.5 C lower in 2007, which was mainly due to a higher than normal occurrence of cloudy weather in the second half of the year (Figures 4.1a, e). As a result,  in the 0-60-cm soil layer in 2007 was much higher than the previous 9-year means (Figure 4.1c). According to Morgenstern et al. (2004) and  Chen et al. (2009b), the weather conditions in 2007 would suggest a possible tree growth limitation due to low Q, and enhancement or reduction of Re and Rs due to high θ or lower Ta, respectively. 104   Figure 4.1 Comparison of monthly mean values of climatic and environmental variables measured at the 58-year-old West Coast Douglas fir stand (DF49) for the fertilized year 2007 and for the 1998-2006 means. (a) Precipitation (P), (b) air temperature above the canopy (Ta), (c) 0-60-cm soil water content (θ), (d) water vapor saturation deficit above the canopy (D) and (e) daytime total downward photosynthetically active radiation (Q). The error bars are ± 1 standard deviation for all the 9-year data (1998-2006). 105   Figure 4.2 Comparison of monthly mean values of component C fluxes measured at the 58- year-old West Coast Douglas-fir stand (DF49) for the fertilized year 2007, for the previous year 2006 and for the 1998-2006 means. (a) Gross primary productivity (GPP), (b) ecosystem respiration (Re), (c) belowground ecosystem (soil) respiration (Rs), and (d) net ecosystem productivity (NEP).  The error bars are ± 1 standard deviation for all the 9-year data (1998-2006). 106  4.4.2 Comparison of measured C component fluxes between pre- and post- fertilization year Figure 4.2 compares the variations in the monthly values of C component fluxes in 2006 and 2007 with the mean values for 1998-2006. Overall, the monthly values of GPP in 2007 were close to the previous 9-year average values (Figure 4.2a). However, the monthly Re values in 2007 were similar to the previous 9-year averaged values before May but lower for May and the later months (Figure 4.2b). The values of below-ground ecosystem respiration (i.e. soil respiration, Rs) values in 2007 were higher for January through May but lower for June, September and October compared with the2003-2006 mean values (Figure 4.2c). As a result, the monthly values of NEP in 2007 were higher than the 1998-2006 mean values during April to December (Figure 4.2d). However, it is critical to recognize that besides N addition, environmental factors play an important role in C cycling processes. Climate perturbations may result in anomalies in C component fluxes. Table 4.2 Comparison of measured annual C component fluxes between non-fertilization year (1998-2006 and the fertilization year (2007)  *  NEP GPP Re Re/GPP 1998-2006 2007 1998-2006 2007 1998-2006 2007 1998-2006 2007 356  ± 51 552 ± 19 2124 ± 125 2095 ± 22 1768 ± 146 1543 ± 17 0.8 ±  0.03 0.74 The ± 1 standard deviation (s.d.) of annual values over 1998-2006 and the standard errors (s.e.) for the year 2007 estimated using the method of Schwalm et al. (2007) are also shown. * Urea at 200 kg N ha -1  applied on 13 January 2007  As shown in Table 4.2, despite 2007 being wetter and cooler than previous years, annual NEP in 2007 was much higher than the average for the previous 9 years, while the Re/GPP ratio for 2007 was lower than the 1998-2006 mean value because there was a much larger decrease in Re than GPP.  107  4.4.3 Response of needle mass, N concentration and growth increment to fertilization The averaged N contents in current-year needles (averages of lower, middle and upper sections of the live crown of each tree crown) in the unfertilized and fertilized trees were respectively 1.15 and 1.61% (dry needle basis). Foliar N analysis showed that the unfertilized trees were severely deficient in N (Ballard and Carter, 1996). As expected improved nutrition resulted in that foliar biomass increased from 445 to 559 mg per100 dry needles on average. Tree ring analysis showed that fertilization resulted in higher annual growth increment. The annual bole diameter increments at the 1.3-m height during the first year following fertilization was 1.15 mm for the fertilized trees, which were significantly higher (p < 0.05) than 0.78 mm for the trees in the control area, which is consistent with the modeled increases in GPP at this site. The fertilized-induced additional MAI was calculated to be 5.4 and 4.4 m 3  ha -1 , respectively, following the method of Feller (1992) and using a growth-and-yield model of Marshall and Turnblom (2005). By assuming that Douglas-fir dry wood density is 450 kg m -3  (Gartner et al., 2002), the increase in stem wood biomass due to fertilization was further approximated to be 1980 to 2430 kg ha -1 . We routinely measured LAI every growing season at 8 plots at the research site. The measured LAI was about 7.3 (± 15%) over the past 10 years (Chen et al., 2006) and no significant difference was found between before and after fertilization likely owing to the stage of near-end-of-rotation. We also found no significant change in fraction of absorbed photosynthetically active radiation (fAPAR) after fertilization. 4.4.4 Discerning N fertilization effects on C component fluxes using the optimized BEPS model The optimized BEPS model with inversed parameters was first calibrated using the EnKF applied 8-year data and further verified using the reserved non-fertilized data in 2006. Finally, the validated model was used to discern N fertilization impacts on C sequestration in 2007.  108  4.4.4.1 Calibration and verification of the optimized BEPS model To evaluate the accuracy of the optimized BEPS model in the prediction of the C component fluxes, the modeled results were calibrated and verified against the measurements made during 1998-2005 and during 2006, respectively.  We also found the optimized parameters using 8-year-long dataset (1998-2005) and using 9-year-long dataset (1998-2006) are identical through a model experiment, implying that the 8-year-long dataset is long enough to reveal the inherent seasonal relationship between the modeled C fluxes and the driving variables.  Figure 4.3 BEPS predicted and EC tower measured daily gross primary productivity (GPP), ecosystem respiration (Re) and net ecosystem productivity (NEP) at DF49 for1998 to 2007, (a) for the EnKF applied years1998-2005, (b) for the model validation year 2006 and (c) for the first-post-fertilization year 2007. 109   Figure 4.4 Comparison of BEPS predicted and EC tower measured daily gross primary productivity (GPP), ecosystem respiration (Re) and net ecosystem productivity (NEP) at DF49 for non-fertilized years, (a) - (c) for the EnKF applied years 1998-2005 and (d) - (f) for the validation year 2006. 110  Comparisons of measured and BEPS-simulated daily GPP, Re and NEP for 1998-2005 are shown in Figure 4.3a and Figure 4.4a -c.  BEPS simulations explained about 89% of the variance of daily GPP, 90% of the variance of daily Re, and 71% of the variance of daily NEP (Figure 4.4a-c). Linear regression analysis of simulated vs. measured daily values (Figure 4a-c) indicated that the  p value was less than  0.001 for all C component fluxes; the slopes were 1.01, 1.02 and 1.03 for GPP, Re and NEP, respectively, and RMSE values were 1.12, 1.05 and 1.26 g C m -2  day -1  for GPP, Re and NEP, respectively. Figure 4.5a compares the modeled and chamber-measured values of Rs. Regression analysis for the available measured period (2003-2005) showed that p value < 0.001; slope = 0.98, r 2  = 0.91 and RMSE = 0.46 g C m -2  day -1 . These model results suggested that the BEPS model reasonably simulated both photosynthesis and respiration processes at a daily time step using the calibrated parameters.  Figure 4.5 BEPS predicted and measured daily soil respiration (Rs) using an auto-chamber at DF49, (a) for the ENKF applied years 2003-2005, (b) for the model validation year 2006 and (c) for the first-post-fertilization year 2007. Figure 4.6 and Table 4.3 compares measured and modeled annual sums of the C component fluxes. Linear regression results showed p <0.001; r 2  = 0.9, 0.98 and 0.95 for GPP, Re and 111  NEP, respectively. Model biases (i.e. measured – modeled values) are significantly smaller (p <0.01) than the departures in C fluxes from their respective 1998-2006 means, suggesting BEPS reasonably captured interannual variations of observed C component fluxes.  Figure 4.6 Comparisons of measured and modeled annual gross primary productivity (GPP), ecosystem respiration (Re), soil respiration (Rs), and net ecosystem production (NEP). The error bars show the standard errors (s.e.). The s.e. values for the EC measured annual C fluxes were estimated using the method of Schwalm et al. (2007), while the s.e. values for modeled annual values were estimated using the ±1 standard deviation of optimized model parameters. (a) for the EnKF applied years, (b) for the model validation year 2006 and (c) for the first-post-fertilization year 2007. 112  Table 4.3 Comparisons of measured and modeled annual gross primary productivity (GPP), net ecosystem productivity (NEP), ecosystem respiration (Re) and  soil respiration (Rs) * * Urea at 200 kg N ha -1  applied on 13 January 2007. The unit for GPP, NEP, Re and Rs is g m -2  yr -1 . dep = departures in C fluxes from their respective1998-2006 mean; bias = measured – modeled; percentage of bias (bias%) = (measured - modeled)/model*100%; Σ = the average of the absolute values of the deviations of data points from their mean for 1998-2005; σ = standard deviation of annual values over 1998-2005. The bias and bias% for 2007 are approximated as absolute and relative N fertilization effects on respective C fluxes, respectively, and positive values indicate positive effects, visa versa  Figure 4.7 compares modeled and measured half-hourly NEP for three separate weeks in 2006, which were selected to show the model’s performance in different phases of the growing season (early, middle and late) and under different weather conditions (clear, cloudy and rainy). Modeled half-hourly NEP followed the observations reasonably except the spikes in observations during some nights. Linear regression analysis for half-hourly NEP showed that p < 0.001; r 2  = 0.76 and RMSE = 4.56 μmol m-2 s-1. Data-model plots of daily C fluxes and linear regression results are shown in Figure 4.4 d-f, which are similar to data-model comparisons for the model calibration period 1998-2005 (Figure 4.4). Relative annual model biases for 2006 were -0.9%, -1.7% and 2.9% for GPP, Re and NEP are similar to the respective values for 1998-2005 (Table 4.3). After verification against measurements at half-  GPP Re Rs NEP  dep bias bias% dep bias bias% bias bias% dep bias bias% 1998 7 29 1.4 -16 11 0.6   23 18 5.0 1999 -100 -11 -0.5 -126 -36 -2.1   26 25 7.0 2000 -33 -31 -1.5 -75 -17 -1.0   42 -14 -3.4 2001 -47 26 1.3 -101 16 1.0   54 10 2.5 2002 -172 14 0.7 -92 31 1.9   -80 -17 -5.8 2003 -46 -16 -0.8 -43 -38 -2.2   -3 22 6.6 2004 214 22 0.9 303 36 1.8 16 1.7 -89 -14 -5.0 2005 186 -17 -0.7 187 -31 -1.6 -8 -0.8 -1 14 4.1 Σ 101 21 1.0 120 27 1.5 12 1.2 41 15 4.6 σ 134 23 1.1 155 31 1.7 17 1.8 54 18 5.3 2006 -12 -19 -0.9 -41 -30 -1.7 13 1.4 29 11 2.9 2007 -29 157 8.1 -225 -93 -5.7 -62 -6.3 196 250 82.8 113  hourly to annual time scales, we have confidence that the BEPS model with the optimized parameters performed well in predicting C component fluxes.  Figure 4.7 BEPS predicted and EC tower measured net ecosystem productivity (NEP) in the early (a), middle (b) and late (c) stages of the growing season in 2006 at DF49. 114  4.4.4.2 Discerning N fertilization impacts After optimization, BEPS was applied to predict C fluxes for 2007 assuming the stand was not fertilized.  The N fertilization impacts in 2007 were then obtained as differences in C component fluxes between measured (fertilized) and the predicted values (unfertilized). As shown in Figures. 4.3-4.6, 4.8, the differences between modeled and measured C component fluxes at different time scales in 2007 were significantly larger than that in other unfertilized years (student t-test showed the regression slopes for 2007 were significantly different from all other years with p <0.01. BEPS modeled daily values were lower than measurements for GPP (Figure 4.8a, slope of predicted vs. observed = 0.88) and higher for Re (Figure 4.8b, slope of predicted vs. observed = 1.06), and as a result, it simulated daily NEP values were lower than measured values (Figure 4.8c, slope of predicted vs. observed = 0.69). BEPS predicted daily Rs values were also higher than measured values (Figure 4.5c, slope of predicted vs. observed = 1.07, not shown). Figure 4.9 shows 7-day running averages of measured and modeled C fluxes, and N fertilization effects on C component fluxes, which were quantified by the differences between measured and modeled values. For most days in 2007, N fertilization had negative effects on Re and Rs, while having positive effects on GPP and NEP. Large N effects were found in summer and fall. Figure 4.10 shows the cumulative plot for each of the C component fluxes for 2007 with error bars, which were calculated using ±1standard deviation of the optimized parameters.  The overall N effects on C fluxes in the first-post fertilization year are shown in Table 4.3 and Figure 4.6 and summarized as follows: (i) annual GPP increased by ~8.1% from 1938 ± 22 g m -2  yr -1  to 2095 ± 29 g m -2  yr - 1 , (ii) annual Re decreased by 5.7% from 1636 ± 17  g m -2  yr -1  to 1543 ± 31 g m -2  yr -1 , and (iii) as a result, NEP increased by 82.8% from 302 ±19 g m -2  yr -1  to 552 ± 36 g m -2  yr -1 . Moreover, in terms of respiration, the model results showed that fertilizer-induced reduction in total Re (~93 g m -2  yr -1 ) mostly resulting from the decrease in belowground Re (i.e. Rs, ~62 g m -2  yr -1 ).   115   Figure 4.8 Comparison of BEPS predicted (no fertilization) and EC tower measured daily gross primary productivity (GPP), ecosystem respiration (Re) and net ecosystem productivity (NEP) at DF49 for the fertilized year 2007. 116   Figure 4.9 Seven-day moving averages of BEPS predicted and measured C component fluxes for the fertilization year 2007 at DF49. The C component fluxes are (a) ecosystem respiration (Re), (b) soil respiration (Rs), (c) gross primary productivity (GPP), and (d) net ecosystem respiration (NEP). The N fertilization effects on the C component fluxes were estimated as the differences between the measured C fluxes and their corresponding modeled values.  117   Figure 4.10 Cumulative N fertilization effects on C component fluxes with ±1 standard errors in 2007 at DF49, which were estimated as the differences between the measured C fluxes and the corresponding modeled values. Standard errors were estimated using the ±1 standard deviation of optimized model parameters. The C component fluxes are gross primary productivity (GPP), ecosystem respiration (Re), soil respiration (Rs) and net ecosystem respiration (NEP). 4.5 Discussion 4.5.1  Forest fertilization effects on photosynthesis Fertilization application in such N-limited stands likely stimulates above-ground NPP (Fisher and Binkley, 2000; Chapin et al., 2002). After fertilization, needle mass at the end of the first year increased by ~26%. Our model results also showed fertilization increased GPP by ~8.1% (= 157 g m -2  yr -1 ) in this stand. These results are consistent with, but lower than, the results for this same site obtained by Jassal et al. (2010a) that N-fertilization increased GPP by 11% (=203 g m -2  yr -1 ). However, in their study, NEE was partitioned by calculating daytime respiration using the relationship between daytime NEP and photosynthetically active radiation and the monthly values of unfertilized GPP for 2007 were calculated using a simple linear regression analysis based on monthly values in the non-fertilized years 1998- 118  2006. These findings are also consistent with other studies: Brix (1991) reported an increase in net photosynthesis rate in response to improved nutrition in Douglas-fir; Canary et al. (2000) report that N fertilization of Douglas-fir plantations in western Washington State added an average of 26.7 Mg ha -1  to the live tree component over a 16-year-period comparing with adjacent unfertilized control sites;  Adams et al. (2005) reported a similar result for the same forest as Canary et al. (2000) reported that the N-fertilized sites (161 Mg C ha -1 ) had an average of 20% more C in the tree biomass compared to unfertilized sites (135 Mg C ha -1 ) over an 8-year-period, suggesting enhancement of  tree growth due to N fertilization; Footen et al. (2009) found that N fertilization increased site productivity of young Douglas-fir stands on low quality sites in the Pacific Northwest 15-22 years after application by a carryover effect; Hyvönen et al. (2007) found the positive effects of fertiliser N on C stocks in trees (stems, stumps, branches, needles, and coarse roots) and these effects were quantified by analysing data from 15 long-term (14-30 years) experiments in Picea abies and Pinus sylvestris stands in Sweden and Finland; Olsson et al. (2005) reported a three-fold increase in above-ground productivity in response to fertilization of a boreal Norway spruce stand; Magill et al. (2004) found that long-term fertilization at Harvard Forest resulted in an increase in aboveground NPP of the high N hardwood stand relative to the control plot, and Leggett and Kelting (2006) found N fertilization increased both aboveground and belowground biomass of Loblolly pine plantations. Brooks et al. (2009) found that latewood Δ13C sharply decreased by 1.4‰ after fertilization and was significantly lower than controls for four years, but no differences existed between fertilization levels, and the effect disappeared after four years in an 85-year-old Douglas-fir stand based on tree ring and isotope analyses. They further concluded that these findings indicate intrinsic water use efficiency increased in response to fertilization (Brooks et al., 2009). Jassal et al., (2010a), on contrast, showed that there was no first-year response of evapotranspiration and water use efficiency to N fertilization in a 58-year-old Douglas fir stand. These findings may imply the N-fertilization had little first-year response on leaf stomatal conductance. The fertilizer-induced increases in photosynthesis may be related to the effect of leaf N on maxcV . The net photosynthetic rate Anet at the leaf level is a function of 119  two tightly-correlated parameters maxcV and maxcJ (see Eqs. (A1) - (A2)).  In nutrient-limited forests, such as West Coast Douglas-fir stands, Anet is generally limited by Ac, while Ac is dominantly controlled by maxcV  (e.g. Arain et al., 2006; see also Eq. A1a).  As many research results (e.g. Arain et al., 2006) have shown maxcV is a nonlinear function of leaf N status,  rc NNV 8.1exp(1)(max  ,                                                                                           (4.2) where α is the maximum value of maxcV  and Nr is leaf Rubisco-N content (in g N m -2  leaf area) of the canopy. Substituting the above values of foliar N of unfertilized trees (1.15%) and fertilized trees (1.61%) into Eq. (4.2) indicates that the fertilizer-induced increase in maxcV is 8%, which is similar to modeled increase in GPP (7.9%) due to N fertilization (Table 4.3 & Figure 4.10). 4.5.2 Forest fertilization effects on respiration The results of our model simulation indicate that Re decreased in the first year after fertilization by ~93 g m -2  yr -1  (~ 5.7%), in which, the aboveground Re declined by ~31 g m -2  yr -1  (~4.4%) and belowground Re (i.e. Rs) declined ~62 g m -2  yr -1 (~6.3%, Figures 4.5c, 4. 6 & Table 4.3). This contrasts to an increase of 35 g m -2  yr -1  calculated by Jassal et al. (2010b) following a different type of NEE partitioning and empirical modelling as explained above, though, however, the disagreement can partly be due to a reported uncertainty of ~30 g m -2  yr -1  in the EC-measured Re values (e.g. Schwalm et al., 2007). Given a significant increase in GPP (by ~157 g m -2  yr -1 ) and an increase in aboveground NPP (by ~126 g m -2  yr -1 ), it is an unexpectedly weak response of above-ground Re to fertilization (reduction of only ~31 g m -2  yr -1 ). Fertilizer-induced increases in aboveground productivity (Olsson et al., 2005), increases in both aboveground and belowground biomass (Leggett and Kelting, 2006), and increases in aboveground biomass with decreases in belowground biomass (Teskey et al., 1995) have also been reported for other forests.  Reich et al. (1998) reported that mass-based leaf dark respiration rate was positively related to leaf N content for 120  69 species from four functional groups (forbs, broad-leafed trees and shrubs, and needle- leafed conifers) in six biomes traversing the Americas (alpine tundra/subalpine forest, Colorado; cold temperate forest/grassland, Wisconsin; cool temperate forest, North Carolina; desert/shrubland, New Mexico; subtropical forest, South Carolina; and tropical rain forest, Amazonas, Venezuela). Based on their findings, we may deduce that a reduction in above ground Re found in this study might result mainly from stem respiration. Comparatively weak effects of fertilization on aboveground Re with considerable increases in aboveground biomass may be related to the following mechanisms: N fertilization may (i) depress the sensitivity of aboveground Re (i.e. both growth and maintenance respiration) to temperature (i.e. a decline in Q10), (ii) alter the relationship between GPP and growth respiration, and (iii) change patterns of photosynthate allocation which is related to Re, i.e. as summarized by Janssens and Luyssaert (2009) that C allocation tends to shift from fine roots and mycorrhizal symbionts, with a relatively low C:N ratio, to woody biomass with a high ratio. For instance, Koehler et al. (2009) found a shift in C partitioning from below- to aboveground in the N- addition plots in a tropical montane forest in which stem diameter growth was promoted. Our findings that N fertilization decreased annual Rs are consistent with those studies for temperate forests (e.g. Haynes & Gower, 1995; Fahey et al., 1998; Maier & Kress, 2000; Bowden et al., 2000, 2004; Butnor et al., 2003; Burton et al., 2004; Phillips and Fahey, 2007); however, the magnitude of the impact of N fertilization on Rs in this study (a 6.3% decrease in Rs) is at the low end of the range (5 - 40%) of the previous findings. The mechanisms behind N-fertilizer-induced Rs reduction have not been clearly identified because it is usually difficult to separate the contributions of root respiration (Rr) and microbial decomposition to Rs (Hanson et al., 2000; Baggs, 2006; Phillips and Fahey, 2007). We simply assumed that a reduction in Rs attributed to both below ground Ra and Rh. N- fertilizer-induced decreases in Ra are most likely associated with reduction in fine-root biomass after fertilization as fine-root biomass is significantly correlated with Rs (e.g. Davidson et al., 2004; Mo JM et al., 2008). Decreases in fine-root biomass in fertilized forest soils have been reported in numerous studies (Aber and Melillo, 1991; Haynes and Gower, 1995; Boxman et al., 1998; Bowden et al., 2004), and are consistent with plant C allocation 121  theory which states that trees decrease C allocation to roots when nutrient availability is high (Bloom et al., 1985). Janssens et al. (2010) revealed that the average response of Rh to N addition is much more pronounced than that of leaf-litter decomposition alone on the basis a statistical meta-analysis. Averaged over 36 N-manipulation studies in forest ecosystems, N addition decreased Rh by 15%, with responses ranging from a reduction of 57% to a stimulation of 63% (Janssens et al., 2010). N-fertilizer-induced Rh reduction is thought to be related to the soil C:N ratio, it has been found to be associated with a decrease in soil microbial biomass (Arnebrandt et al., 1990; Wallenstein, 2003; Bowden et al., 2004; Compton et al., 2004; Frey et al., 2004) or with an increase in microbial C use efficiency if an increased proportion of C is assimilated into new biomass (Phillips and Fahey, 2007). The findings in this study which show that N fertilization decreased Rs, however, differ from the field observations at Harvard Forest reported by Micks et al. (2004), and manual and automated chamber measurements made at this site by Jassal et al. (2010b). Using manual soil CO2 efflux measurements, made at 2-4-weekly intervals in 16 plots (eight root-exclusion plots (Rh) and eight adjacent control plots (Rs)), Jassal et al. (2010b) showed that fertilization resulted in a significant (20%) increase in Rs during the first 3-4 months due mainly to an increase in autotrophic (or rhizospheric) soil respiration (Rs). In the following 3 months, they found little effect on Rs but a small (~6%) decrease in Rh. They attributed the increase in Ra to fertilization causing an increased production of fast-turnover fine roots (Pan et al., 2009, Cleveland and Townsend 2006), which would stimulate higher C allocation belowground to balance high N concentration in these tissues (Field et al., 1992; Chapin et al., 1990). Thus, a definitive answer on the short-term responses of Rs and its components to N addition is not readily apparent, however regardless of sign, changes in Rs induced by N were small and within uncertainties or biases in both modeling and/or measurement approaches. Future research is therefore required to help resolve the question of N-fertilization effects on Rs and its components.  122  4.5.3 N-use efficiency for C sequestration As Janssens and Luyssaert (2009) summarized, three main mechanisms for a N–C link in forests are:  (i) accelerated photosynthesis, (ii) a change in tree C allocation from roots and root symbionts to woody biomass, and (iii) slowing of the microbial degradation of soil C (these were discussed at a recent workshop ‘The Impact of N on the Forest C Cycle’ held in Stockholm, Sweden, in February, 2009). The response of C sequestration to N addition, in term of N-use efficiency (kg C (sequestrated) kg -1 (N added)) relates to the spatial variations in each of the three mechanisms, in other words, it depends on geographical situation, soil status, tree species, stand age, fertilizer composition and dose (Hyvönen et al., 2008). Magnani et al. (2007) reported a strong positive relation between mean lifetime C sequestration (in terms of net ecosystem production, NEPav) and N deposition for mid- latitude forests based on an analysis of 20 European forest stands using an Arrhenius function. Such a finding differs markedly from other estimates (De Vries et al., 2008; Sutton et al., 2008), and induces a major debate on the relationship between atmospheric N deposition and forest C sequestration (de Vries et al., 2008; De Schrijver et al., 2008; Sutton et al., 2008). Sutton et al. (2008) found that the response of net C sequestration to total N was about 50 - 75:1, by re-analysis of 22 European forest stands and accounting for the effects of inter-site climatological differences. A most common range of C sequestration per kg N addition in above-ground biomass and in soil organic matter for forests was estimated to be about 20-40 kg C/kg N using multiple approaches (de Vries et al., 2009).  By analysing data from 15 long-term (14-30 years) experiments in Picea abies and Pinus sylvestris stands in Sweden and Finland, Hyvönen et al (2008) found low application rates (30–50 kg N ha-1 year -1 ) were always more efficient per unit of N than high application rates (50–200 kg N ha-1 year -1 ), suggesting the size of dose will affect the N-use efficiency. N fertilization with the amount of 200 kg N ha -1  in one single dose is different from N deposition with characteristics of long-term and low-dose N addition.  N-use efficiency in the near-end-of-rotation Douglas fir stand (58-year old in 2007) in this study was 12.5 kg C (sequestrated) kg -1 (N added) in the first year following fertilization, which is expected to be lower than the common range of N- use efficiency (20-40 kg C/kg N) for forests due to atmospheric N deposition (de Vries et al., 123  2009). N fertilization with the amount of 200 kg N ha -1  in one single dose is a common forest management activity in N-limited forests, e.g. Pacific Northwest coastal forests of the United States and Canada covering approximately 10 5  km 2 between Oregon and Alaska. Therefore, to assess the effects of this common N fertilization activity on net C sequestration is also important.  N-use efficiency in the 58-year old Douglas fir stand is lower than that in the stand at middle rotation stage (19-year old) with value of 18 kg C (kg N)  -1 (Jassal et al., 2010a), implying the stand age effect on N-use efficiency. 4.5.4 Uncertainty and limitation of this modeling approach The 9-year EC data (1998-2006) at DF49 before fertilization shows the apparent interannual variability of NEP, GPP and Re were (± s.d.) 357 ± 51, 2,124 ± 125, and 1,767 ± 146 g C m -2  yr -1 , with ranges of 267 - 410, 1,592 – 2,338, and 1,642 – 2,071g C m-2 yr-1, respectively (see Chen et al., 2009a). The major drivers of interannual variability in annual C fluxes were interannual climatic perturbations (especially annual and spring mean air temperatures and water deficiency during late summer and autumn; Chen et al., 2009a). How to distinguish fertilization and climatic effects on C fluxes is a challenge. One may either design measurement or modeling experiments. Given the high spatial heterogeneity of land surface in nature, it is hard to find an EC flux area with similar land surface properties ( e.g. tree age, density, leaf area index and soil properties etc.) to set up a control EC tower as in this case. The EC sensor location bias (a measure of uncertainties of EC tower data owing to the land surface heterogeneity) for DF49 based on footprint climatology and remote sensing analyses was about 15% for GPP estimation (Chen et al., 2009b), which is larger than the magnitude of N fertilization effects on GPP and Re but lower than of the effects on NEP. The measuring experiment with control tower, therefore, is not a pragmatic and reliable approach. For example, Hollinger et al. (2004) had shown the considerable difference in C fluxes between the paired measurements from two flux towers (e.g. Richardson et al., 2006). Soil chamber measuring experiments face a similar spatial representativeness challenge. The responses of forest ecosystems to N addition are likely different among years after fertilization (Aber et al., 1989, 1998). The modeling approach is necessary to explore the N fertilization effects on 124  C sequestration in first-post year in this research. A desirable process-based ecosystem model with high accuracy (model uncertainty must be less than the N fertilization effect on C sequestration, i.e. 5-10%) is required in a modeling approach. This type of model, one with high accuracy, may not exist. In this study, we designed a modeling experiment by taking advantage of the data-model fusion technique to eliminate the effects of interannual variations on C fluxes based on an 9-year-long series of EC data (1998-2006) before fertilization (see Section 4.2.4). A set of key model parameters was optimized  based on the binned month dataset with the ensemble size of ~8  30   48 in order to account for the effects of interannual climatic perturbations but reveal seasonal controlling factors of C fluxes including both biotic and abiotic factors. The objective of model runs is to distinguish the N fertilization effects on C exchanges rather than to obtain good agreement between model outputs and measurements for each modeling year. One may argue about the uncertainty of the N fertilization effects on net and component C fluxes derived from the data-model approach as it is not easy to assess. Alternative methods are worthwhile to explore for future studies for assessing the effects of fertilization, e.g. using a process-based model with a full description of N dynamics. 4.6 Conclusions Fertilization of a 58-year-old Pacific Northwest Douglas-fir forest with 200 kg N ha -1  was found to increase annual NEP by ~83%, from 302 ±19 g m -2  yr -1  to 552 ± 36 g m -2  yr -1  in the first year. On the basis of a data-model fusion approach, this N fertilizer-induced increased in NEP was found to have resulted primarily from increases in annual GPP by ~8%, from 1938 ± 22 g m -2  yr -1  to 2095 ± 29 g m -2  yr -1  and secondly from decreases in annual Re by ~5.7%, from 1636 ± 17  g m -2  yr -1  to 1543 ± 31  g m -2  yr -1 . Modeling results indicated that fertilizer-induced decreases in belowground Re (i.e. Rs) mainly contributed to the reduction in total Re while aboveground Re had less response to N fertilization. The fertilizer-induced enhancement of photosynthesis (~157 g m -2  yr -1 ) was more pronounced than the depression of soil respiration (~62 g m -2  yr -1 ). The former is likely related to the functional relationship of the maximum carboxylation rate to leaf N content. 125  The uncertainties or biases in estimated N effects, stemming from errors of both measurements and modeling results are still considerable but less than 10% of their quantities.  These results suggests that N fertilization as a forest management activity in such N-limited forests and soils, may significantly increase C sequestration and may have potential consequences for feedbacks to global change.   126  Chapter  5: Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements 5  5.1 Synopsis In this Chapter, we describe pragmatic and reliable methods to examine the influence of patch-scale heterogeneities on the uncertainty in long-term EC C flux data and to scale between the C-flux estimates derived from land surface optical remote sensing and directly derived from EC flux measurements on the basis of assessment of footprint climatology. Three different aged Douglas-fir stands with EC flux towers located on Vancouver Island and part of the Fluxnet Canada Research Network were selected. Monthly, annual and inter- annual footprint climatology, unweighted or weighted by C fluxes, were produced by a simple model based on an analytical solution of the Eulerian advection-diffusion equation. The dimensions and orientation of the flux footprint depended on the height of the measurement, surface roughness length, wind speed and direction, and atmospheric stability. The weighted footprint climatology varied with the different C flux components and was asymmetrically distributed around the tower and its size and spatial structure significantly varied monthly, seasonally and inter-annually. The gross primary productivity (GPP) maps at 10 m resolution were produced using a tower mounted multi-angular spectroradiometer, combined with the canopy structural information derived from airborne laser scanning (Light Detection and Ranging --LiDAR data).  The horizontal arrays of footprint climatology were superimposed with the 10-m-resolution GPP maps. Monthly and annual uncertainties in EC flux caused by variations in footprint climatology of a 59-year-old Douglas-fir stand were  5  Reprinted from Boundary-Layer Meteorology, Vol. 130(2), 137–167, [B. Chen], T.A. Black, N.C. Coops, T. Hilker, T. Trofymow, Z. Nesic, K. Morgenstern (2009), Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements. Boundary- Layer Meteorology, 130(2), pages 137-167, (DOI: 10.1007/s10546-008-9339-1) Copyright (2009), with permission from SpringerLink.  Minor updates, rewording and corrections have been applied. 127  estimated to be ~15 - 20% based on comparison of GPP estimates derived from EC and remote sensing measurements and on sensor location bias analysis. The footprint-variation induced uncertainty in long-term EC flux measurements was mainly dependent on the site spatial heterogeneity. The bias in C-flux estimates from spatially-explicit ecological models or tower-based remote sensing at finer scales was estimated by comparing the footprint weighted and EC derived flux estimates. This bias is useful for model parameter optimization. Optimization of parameters in remote-sensing algorithms or ecosystem models using satellite data, will in turn, increase the accuracy in the up-scaled regional C-flux estimation. 5.2 Introduction Growing interest in climate change has stimulated recent research that aims to quantify components of the terrestrial C cycle. Many of these studies have focused on quantifying land-surface-atmosphere fluxes of the primary naturally occurring radiatively active trace gases (Fox et al., 2008), i.e. CO2 and methane (CH4). Although a variety of methods are used in such studies, the most direct measurements of the terrestrial C flux are made either at the plot scale (10 -1 – 102 m2), e.g. using biometric methods and various forms of chamber, or at the ecosystem (patch) scale (10 4 - 10 6  m 2 ), using the eddy covariance (EC) technique. The number of EC flux tower measurements has been rapidly increasing over the world (Black et al., 1996; Baldocchi et al., 2001a). There are over 400 stations globally using the EC method to measure the exchange of CO2, water vapour and energy between the land and the atmosphere. Although the available EC data have been rapidly accumulating, much of this information is of limited use, unless methods of processing and interpretation of such data are available (Sogachev et al., 2004). The EC method is based on measurements of turbulent fluctuations of the vertical velocity and the concentration of a passive tracer. Knowledge of the area’s soil and vegetation that impacts the EC flux is clearly important both in planning the site tower location and in the interpretation of measured fluxes (Finnigan, 2004). The adoption of the EC technique to estimate surface exchange is based on the assumption that certain meteorological conditions (e.g. horizontal homogeneity, steady-state, and non- 128  advection) are satisfied (Göckede et al., 2004). Since such conditions are often violated in complex terrain, e.g. at flux monitoring sites in forests, correct interpretation of EC-data is still a matter of some difficulty (Sogachev et al., 2004). In particular, the spatial variability of vegetation density influences the lower atmospheric circulation and surface exchange of energy, water and C over a wide range of scales (e.g. Shen and Leclerc, 1995; Buermann et al., 2001; Cosh and Brutsaert, 2003). As a result, evaluation of the spatial representativeness of long-term accumulated EC-flux measurements is still challenging. Footprint analysis is a recognized part of the establishment and siting of flux towers and the analysis of their output (Finnigan, 2004). The footprint concept in its modern form is attributed to van Ulden (1978) and Pasquill and Smith (1983), which has been widely used. Horst and Weil (1992, 1994) refined the simple analytical model of Schuepp et al. (1990) by making use of realistic velocity profiles and stability dependence. Schmid (1994, 1997) extended the model to two-dimensions (i.e. FSAM) which has been widely adopted, however only applied during period of neutral and moderately unstable atmospheric stability and with a limited range of crosswind turbulence intensity. Kormann and Meixner (2001) developed a simple analytical model to describe the crosswind integrated and crosswind-distributed footprint under all conditions of atmospheric stability. The Lagrangian stochastic approach was first applied to footprint estimation by Leclerc and Thurtell (1990). Baldocchi (1997b) studied the qualitative effects of canopy turbulence on the footprint using the Lagrangian simulation method with Kljun et al. (2002, 2004) developing a footprint model for a wide range of atmosphere boundary layer conditions. The interpretation of EC flux measurements over a heterogeneous surface depends largely on the footprint over which the fluxes are sampled. The dimensions and orientation of the flux footprint depend on the height of the sensor (zm), surface roughness length (z0), wind speed (u) and direction, and atmospheric stability (e.g. Leclerc and Thurtell 1990; Schmid 2002). In recent years, a number of flux footprint models have been applied to quantitatively estimate the upwind distribution of source-weighting factors for a given short-time (e.g. half hour) (e.g. Schuepp et al., 1990; Leclerc and Thurtell, 1990; Horst and Weil, 1992,1994; Schmid 129  1997; Baldocchi, 1997b; Schmid and Lloyd, 1999; Kormann and Meixner, 2001; Schmid 2002; Leclerc et al., 2003a, b; Soegaard et al., 2003). Despite many current studies on detailed footprint modeling and experimental validation (e.g. Leclerc and Thurtell, 1990; Wilson and Swaters, 1991; Horst and Weil, 1994; Finn et al., 1996; Leclerc et al., 1997, 2003a, b; Cooper et al., 2003), the temporal and spatial variability of footprints and the associated influence of varying site heterogeneities on tower flux measurements has yet to be fully been investigated. One of the practical problems in using a footprint model as an operational tool is that the source contribution in the area of a prospective measurement site is not known a priori (Kim et al., 2006). Recently, it has been demonstrated that long-term (e.g. seasonal, annual and multi-years)  patterns of source contributions (i.e. ‘footprint climatology’) provide essential information about the vegetation sampled when measuring long-term fluxes especially over heterogeneous landscapes (e.g. Amiro, 1998; Schmid and Lloyd, 1999; Stoughton et al., 2000; Kim et al., 2006). However, footprint climatology has not been applied in a practical way to quantitatively illustrate the upwind distribution of source-weighting factors over long time periods. Currently researchers are assessing the agreement between micrometeorological (usually using EC) measurements and biometric estimates of ecosystem-atmosphere C exchange (e.g. Barford et al., 2001; Curtis et al., 2002). The latter has finer spatial scales (10 1  versus 10 6  m 2 ) and has much coarser temporal scales than the former (yearly or multi-years versus half- hourly or hourly).  To validly compare them, correct weighting is critical; therefore the footprint climatology is needed. In order to use EC flux measurements to provide estimates of components of the natural C cycle at landscape (10 1 - 10 2  km 2 ), regional (10 3  - 10 6  km 2 ), or hemispheric to global land surface (10 7  - 10 8 km 2 ) scales, they must be reasonably ‘‘upscaled’’ using either models and/or remote sensing measurements (e.g. Earth observation (EO) data). Upscaling models usually use some form of vegetation or land cover map, often derived from EO data, to assess the areal extent of each surface type. The overall fluxes at larger spatial scales are then estimated by calculating the area-weighted sums of fluxes measured for the various land surface types (e.g. Heikkinen et al., 2004; Soegaard et al., 2000). If such models do not explicitly take into acc